尧图网页设计 
做综合能源系统优化调度的朋友这两年应该没少在论文里碰到“信息间隙决策理论IGDT”和“阶梯碳交易”这两个词。尤其是做双碳方向、园区多能互补、以及电热气氢耦合系统研究的基本绕不开这套组合。这个题目“基于信息间隙决策理论的多能系统-阶梯碳交易优化调度Python代码实现”本质上就是在解决一个问题风光出力不确定、碳价政策又在变这时候你怎么让多能系统既省钱又稳得住还能算得出来。这篇文章就把这套东西从头到尾拆开讲清楚。我会把IGDT的核心思想、多能系统的建模套路、阶梯碳交易的数学处理以及Python代码怎么一步步实现整个优化调度模型全部梳理出来。不管你是刚入门的硕士生还是已经在做园区综合能源项目的工程师只要按着这篇文章的思路走一遍都能自己把代码跑起来并且理解每个环节为什么要这么写。1. 整体思路拆解为什么是IGDT加阶梯碳交易1.1 多能系统优化调度的本质是什么先聊清楚多能系统到底在算什么账。所谓多能系统简单说就是把电、气、热、冷这几类能源放在一个平台里统一调度常见的设备包括燃气轮机、燃气锅炉、余热回收、电锅炉、储能、光伏、风电有时候还带P2G电转气和氢能。系统内部通过能源转换设备把一种能源变成另一种比如天然气进去出来的是电和热光伏发出来的电既可以供给电负荷也可以驱动电锅炉制热。这种系统的优化调度核心目标就是在满足各类负荷需求的前提下决定每个时段每台设备的出力使得总成本最低。这里面有几个关键约束能量平衡约束电、热、气各自守恒、设备出力上下限约束、爬坡约束、储能充放能约束以及碳排放相关约束。写成一个标准的数学优化问题就是一个带连续变量和0-1变量的混合整数线性规划MILP问题。但问题是现实场景里光伏和风电的出力不是精确知道的值负荷预测也是存在误差的。传统做法是假设这些不确定量服从某种概率分布然后用随机优化或者场景法处理但概率分布本身往往很难准确获得。这里IGDT的优势就体现出来了。1.2 信息间隙决策理论到底解决什么问题信息间隙决策理论英文叫Info-Gap Decision Theory第一次看到这个名字会觉得挺唬人其实思想非常朴素。它假设你手里的预测数据比如明天光伏出力预测是100MW和真实值之间可能存在误差但这个误差有多大你不知道也不需要你去猜一个具体的概率分布。你只需要定义一个不确定半径α表示预测值在多大范围内浮动比如α0.1就表示真实值落在预测值上下10%的区间内。IGDT一般分两种决策模式。一种是风险规避型就是我做调度决策时要保证即使不确定量往最坏的方向变化系统的运行成本也不会超过某个我可以接受的阈值在这个前提下系统能承受多大的不确定性也就是α能取到多大越大说明系统抗风险能力越强。另一种是风险追求型也叫机会型是反过来利用不确定性期望在实际值往好的方向偏差时能获得额外收益此时关注的是最小可能成本。把IGDT用到多能系统调度里步骤也很清晰。先求一个确定性模型的最优解得到一个基准成本C0。然后引入决策者的风险偏好系数β表示你愿意在基准成本的基础上多付出β比例的费用来换取抵御不确定性的能力。接下来构建一个鲁棒优化模型在不确定参数的区间范围内找到使得系统成本恰好不超过(1β)C0的最大不确定性半径α。这个α值越大说明系统越“抗造”。通过对比不同β下的α就能得到一组决策者可以根据自身风险态度来选择的帕累托前沿。1.3 碳交易为什么要做成阶梯式至于碳排放这块传统的碳交易模型大多采用固定碳价或者单一碳价每排放一吨碳就付一笔固定费用。这种做法的缺陷在于它对碳排放量大的主体约束力不够强排得多和排得少在边际成本上没有区别。阶梯碳交易机制就聪明在这里它像阶梯电价一样把碳排放配额按区间划分超出基准配额的排放量越多对应的碳价越高碳排放的边际成本是递增的。具体建模时先用基准线法确定系统的免费碳排放配额这个配额通常跟机组的发电量或者装机容量挂钩。然后计算系统实际碳排放量包括外购电力的等效碳排放和燃气设备燃烧产生的直接碳排放。实际排放量减去免费配额差额如果为正就要去碳市场买配额按照阶梯价格付费如果为负说明有多余的配额可以卖出获利。阶梯碳交易函数是一个分段线性函数在优化模型里需要引入0-1变量来分段处理这会让模型变成一个混合整数规划。但这部分在Python里用Gurobi或者Cplex求解并不困难关键是分段点的设置和松弛处理后面代码部分我会详细讲。1.4 三块拼起来为什么能自洽把这三件事拼在一起的逻辑其实是闭环的。多能系统因为风光不确定导致运行成本波动IGDT给出了一个量化这种风险的工具碳排放约束通过阶梯碳交易嵌入到成本函数中让碳排成本成为调度决策的一部分两者合在一起目标函数就成了“燃料成本购电成本设备运维成本碳交易成本”然后在这个目标下寻找既经济又稳健的运行方案。从学术角度看这种组合能解决单一方法解决不了的问题只有IGDT没有碳交易模型不够完整不符合双碳政策导向只有碳交易没有IGDT无法处理风光不确定性两者结合恰好覆盖了多能系统优化调度领域最核心的两个痛点。这也是为什么近两年这个方向的文章特别多你如果做研究用这个思路很容易出成果而且代码跑出来的结果会非常有说服力——对比确定性模型和IGDT模型的结果画几张图就能直观展示鲁棒性和经济性之间的权衡。2. 数学建模全流程从设备模型到目标函数2.1 典型多能系统物理架构与设备建模我假设系统中的设备包括燃气轮机CHP、燃气锅炉GB、电锅炉EB、光伏PV、风电WT、蓄电池BESS、蓄热罐TESS以及外网购电。这样一个架构覆盖了大部分园区级综合能源系统的典型配置代码写出来也方便你后续往里面加电转气、氢燃料电池、冷负荷等模块。设备建模时需要注意几个关键量的单位统一这是个特别容易踩坑的地方。功率类变量统一用MW时间尺度选1小时能量类变量用MWh成本类金额用元碳排量用吨。只要单位统一了后面写约束才不会出数值问题。燃气轮机的建模核心是热电联产关系。抽凝式CHP的可运行域是一个四边形区域电出力和热出力之间存在耦合关系通常用几个不等式来刻画# CHP可运行区域约束示例电量P_chp热量H_chp # 区域顶点单位MW由机组参数决定此处为四个约束限制 m.addConstr(P_chp 0) m.addConstr(H_chp 0) m.addConstr(P_chp - 0.6 * H_chp 20) m.addConstr(P_chp 0.4 * H_chp 120)如果你用的是背压式CHP那热电比是一个固定值模型更简单热出力就等于电出力乘上热电比。燃气锅炉就是典型的线性模型天然气输入热量转换成热出力效率一般取0.85到0.9。电锅炉也类似效率取0.95左右。储电和储热模型则采用荷电状态SOC递推公式SOC[t1] SOC[t] P_charge[t] * eta_charge - P_discharge[t] / eta_discharge这里有几个防止模型吹牛的约束必须写不然求出来的解根本没有实际意义。储能设备同一时段不能同时充放电需要用0-1变量限制初始和末尾SOC要相等这是调度周期闭环的要求尤其是做日前调度时一天结束和一天开始的储能量应当一致否则第二天没法继续。2.2 目标函数燃料购电运维碳交易四本账目标函数的设计是整个模型的灵魂。我推荐把总成本写成四项之和第一项是燃料成本主要包括天然气的消耗成本。燃气轮机的气耗量可以近似表示为其出力的线性函数# 燃料成本单位天然气价格gas_price元/MWh # CHP和燃气锅炉消耗的天然气功率MWh/h fuel_cost gas_price * sum(GT_gas[t] GB_gas[t] for t in range(T))第二项是购电成本从上级电网购电的价格按照分时电价处理峰谷平时段价格不同electricity_cost sum(price_buy[t] * P_buy[t] for t in range(T))第三项是设备运维成本一般取各设备出力的比例系数光伏和风电的运维成本很低CHP和锅炉稍高储能取充放电功率的成本。这一项虽然单项不大但加上去会让模型更真实。第四项就是阶梯碳交易成本这部分是重头戏下面单独拆解。目标函数最终写成total_cost fuel_cost electricity_cost maintenance_cost carbon_cost2.3 碳交易成本的计算逻辑与数学表达式碳排放量的计算分两部分。燃气轮机和燃气锅炉燃烧天然气产生直接碳排放排放因子取单位热值天然气的CO2排放系数大约是0.2吨CO2/MWh不同文献取值略有差异你按自己的参考标准就行。外购电力虽然不在本地直接排放但在全生命周期口径下电网端的火电占比决定了它的等效碳排放业内一般取0.58到0.997吨CO2/MWh看你所在区域的电网排放因子。免费碳排放配额的确定方式比较多常见做法是按照系统的总输出能量乘以一个配额系数。我代码里用的是按CHP电出力和热出力、锅炉热出力分别乘以各自的配额系数再加总。这样做的合理性在于多能系统里不同设备的能源利用效率不同给效率高的CHP更高的配额系数相当于政策上鼓励热电联产。计算完实际碳排放量E_actual和免费配额E_quota之后需要购买的碳排放额度E_purchase等于E_actual减去E_quota的差额。阶梯碳交易的过程就是针对这个差额分段计价。我设置三个阶梯第一阶梯是0到E1碳价为c1第二阶梯是E1到E2碳价为c2第三阶梯是超过E2的部分碳价为c3。c2比c1高c3比c2高体现碳排放边际成本递增的政策意图。在Python代码中实现这个分段计价我推荐用Gurobi自带的分段线性约束addGenConstrPWL这样最稳妥不需要手动引入0-1变量代码简洁且不容易出错。如果你用的是其他求解器那就老老实实引入0-1变量做混合整数线性化逻辑也不复杂就是把每个阶梯段的选择用二进制变量表示并加相应的约束限制该段的取值范围。2.4 IGDT鲁棒模型的构建与转换过程IGDT模型的构建需要先做一次确定性求解。以光伏预测值PV_forecast、风电预测值WT_forecast和负荷预测值Load_forecast为输入求出基准成本C0。然后定义不确定参数的实际值落在预测值的(1-α)到(1α)区间内PV_actual[t] (1 - alpha * delta_PV[t]) * PV_forecast[t] # 风险规避下取最坏情况 WT_actual[t] (1 - alpha * delta_WT[t]) * WT_forecast[t] # 同上 Load_actual[t] (1 alpha * delta_Load[t]) * Load_forecast[t] # 负荷往大偏风险规避型IGDT的逻辑是在最坏的不确定性情景下优化调度方案使得系统成本不超过(1β)C0并且在此约束下最大化α。注意这里最坏情况的方向是有讲究的。光伏和风电往小偏是坏方向因为发电量少了系统得多买电或者多烧气负荷往大偏是坏方向因为要多供应能量。所以鲁棒模型中不确定参数取极值时光伏取(1-α)倍风电取(1-α)倍负荷取(1α)倍。这样IGDT模型就变成了一个两层问题。外层最大化α内层在给定α下求最优调度。实际代码处理时可以采用迭代法从α0开始每次增加一个步长比如0.02在固定α下求解一次带鲁棒约束的优化模型看最小成本是否突破了(1β)C0。如果没突破就继续增大α如果突破了说明系统能承受的最大不确定性半径就是上一个α值。这个迭代法虽然慢一点但非常直观而且每步求解的都是一个MILPGurobi处理起来没有压力。固定α下单次求解的模型就是把确定性模型中的光伏、风电、负荷替换成受α影响的值然后加一条总成本不超过(1β)C0的约束即可。目标函数此时可以设成就地最小化成本也可以直接设为一个常数反正只是判断可行性。我习惯用最小化成本作为目标这样还能顺便拿到该α下的调度方案。3. Python代码实现从环境配置到跑通全流程3.1 编程环境与依赖库安装代码用的是Python 3.8以上版本求解器是Gurobi。Gurobi的学术许可申请很快用校园网邮箱注册就能拿到。如果你拿不到Gurobi也可以把模型改写成用Cbc或者HiGHS求解但MILP求解性能会差一些尤其是当系统规模扩大后差距会非常明显。需要安装的库只有这么几个pip install gurobipy pandas numpy matplotlibpandas用来读Excel数据文件numpy处理数值计算matplotlib画结果图。如果电脑里还没装Python直接装Anaconda发行版最省事自带大部分科学计算库Gurobi装完后把安装目录下的gurobipy路径加进去或者直接用pip安装gurobipy就行现在Gurobi 10以上的版本pip安装已经非常成熟了。3.2 数据准备与参数配置数据的组织方式是整个代码能不能跑通的关键。我习惯把所有系统参数集中放在一个配置模块里这样改起来方便不用在模型代码里到处翻找某个常数的定义。核心参数包括负荷数据、光伏预测数据、风电预测数据、分时电价这些我放在三个DataFrame里另外设备参数单独定义成字典。分时电价这块尤其要注意不同地区的峰谷时段划分不同而且冬夏可能不一样。代码里用一个列表price_buy表示24个时段的购电价格峰段时间设为1.2元/kWh平时段0.7元/kWh谷时段0.3元/kWh1元/kWh等于1000元/MWh换算关系别搞混了。IGDT的参数配置就是两个核心值决策者风险偏好系数beta和不确定性步长alpha_step。beta一般取0.05到0.3之间分别代表愿意多花5%到30%的成本来换取抗不确定性能力。alpha_step取0.02这样算出来的α精度到0.02对大多数分析已经够用了你要是想画更细致的曲线可以取0.01但求解时间会翻倍。这里我给一个示例的参数字典你运行代码前需要把这个模块替换成你自己系统的数据# 系统基本参数单位统一见注释 T 24 # 调度周期24小时 params { GT_max: 100, # 燃气轮机最大电出力 MW GT_min: 20, # 燃气轮机最小电出力 MW GT_eta_e: 0.35, # 燃气轮机发电效率 GT_eta_h: 0.45, # 燃气轮机余热回收效率 GB_max: 80, # 燃气锅炉最大热出力 MW GB_eta: 0.88, # 燃气锅炉效率 EB_max: 50, # 电锅炉最大热出力 MW EB_eta: 0.95, # 电锅炉效率 BESS_cap: 80, # 蓄电池容量 MWh BESS_pmax: 20, # 蓄电池最大充放电功率 MW BESS_eta: 0.92, # 蓄电池充放电效率 TESS_cap: 120, # 蓄热罐容量 MWh TESS_pmax: 30, # 蓄热罐最大充放热功率 MW TESS_eta: 0.90, # 蓄热罐充放热效率 gas_price: 2.8, # 天然气价格 元/m3换算后约等于元/MWh CO2_factor_gas: 0.2, # 天然气燃烧CO2排放系数 吨/MWh CO2_factor_grid: 0.58, # 电网购电等效CO2排放系数 吨/MWh carbon_base_price: 60, # 基础碳价 元/吨 carbon_step_price1: 60, # 第一阶梯碳价 元/吨 carbon_step_price2: 80, # 第二阶梯碳价 元/吨 carbon_step_price3: 100, # 第三阶梯碳价 元/吨 carbon_interval: 50, # 每个阶梯的碳排放区间宽度 吨 }3.3 确定性模型的核心代码实现确定性模型是整个代码框架的地基IGDT模型只是在这个基础上做循环修改。模型主体用gurobipy来写先创建模型对象然后创建变量再写约束和目标函数。变量的定义方式我给出关键部分import gurobipy as gp from gurobipy import GRB m gp.Model(MultiEnergySystem) # 连续变量 P_gt m.addVars(T, lb0, ub100, nameP_gt) # CHP电出力 H_gt m.addVars(T, lb0, ub150, nameH_gt) # CHP热出力 H_gb m.addVars(T, lb0, ub80, nameH_gb) # 燃气锅炉热出力 H_eb m.addVars(T, lb0, ub50, nameH_eb) # 电锅炉热出力 P_buy m.addVars(T, lb0, ub150, nameP_buy) # 购电功率 P_bess_c m.addVars(T, lb0, ub20, nameP_bess_c) # 蓄电充电 P_bess_d m.addVars(T, lb0, ub20, nameP_bess_d) # 蓄电放电 H_tess_c m.addVars(T, lb0, ub30, nameH_tess_c) # 蓄热充热 H_tess_d m.addVars(T, lb0, ub30, nameH_tess_d) # 蓄热放热 SOC_bess m.addVars(T1, lb0, ub80, nameSOC_bess) # 蓄电池SOC SOC_tess m.addVars(T1, lb0, ub120, nameSOC_tess) # 蓄热罐SOC # 二进制变量 u_gt m.addVars(T, vtypeGRB.BINARY, nameu_gt) # CHP启停 u_gb m.addVars(T, vtypeGRB.BINARY, nameu_gb) # 锅炉启停 u_bess_c m.addVars(T, vtypeGRB.BINARY, nameu_bess_c) # 蓄电充电状态 u_bess_d m.addVars(T, vtypeGRB.BINARY, nameu_bess_d) # 蓄电放电状态 u_tess_c m.addVars(T, vtypeGRB.BINARY, nameu_tess_c) u_tess_d m.addVars(T, vtypeGRB.BINARY, nameu_tess_d)电功率平衡约束是系统所有能量流的核心约束它表达的是电源出力等于电负荷加电锅炉消耗加储能净充电for t in range(T): m.addConstr( P_gt[t] P_wind_actual[t] P_pv_actual[t] P_buy[t] P_bess_d[t] P_load_actual[t] H_eb[t] / EB_eta P_bess_c[t], namefpower_balance_{t} )热功率平衡约束是热源出力等于热负荷加储热净充热for t in range(T): m.addConstr( H_gt[t] H_gb[t] H_eb[t] H_tess_d[t] H_load_actual[t] H_tess_c[t], namefheat_balance_{t} )这里有个细节要注意电锅炉消耗的是电功率产出的却是热功率要不是在电平衡约束里除以EB_eta要不是在热平衡约束里乘EB_eta方向千万别搞反了。我习惯在电平衡里除以效率这样H_eb变量的物理意义是产出的热功率和热平衡约束的写法保持一致。很多初学者第一次跑模型结果报不可行排了半天查不出问题结果就是电锅炉效率这里除法和乘法用反了。CHP的燃气消耗量和电出力、热出力之间的关系我直接用经验公式。对于抽凝式机组在运行区域内气耗量近似等于电出力除以发电效率减去一个修正项但更稳妥的做法是把CHP的燃料消耗表示成电出力和热出力的线性函数。你可以根据具体机组的性能曲线拟合我这里用简单模型for t in range(T): m.addConstr(GT_gas[t] P_gt[t] / GT_eta_e 0.15 * H_gt[t], namefGT_fuel_{t})燃气锅炉的气耗量就简单了热出力除以锅炉效率for t in range(T): m.addConstr(GB_gas[t] H_gb[t] / GB_eta, namefGB_fuel_{t})储能约束除了SOC递推还要限制同一时段不能同时充放。这个用两个二进制变量来做for t in range(T): m.addConstr(P_bess_c[t] BESS_pmax * u_bess_c[t], namefbess_c_ub_{t}) m.addConstr(P_bess_d[t] BESS_pmax * u_bess_d[t], namefbess_d_ub_{t}) m.addConstr(u_bess_c[t] u_bess_d[t] 1, namefbess_state_{t}) # 初始和末尾SOC一致 m.addConstr(SOC_bess[0] SOC_bess[T], namebess_cycle)3.4 阶梯碳交易的代码实现碳交易这部分是模型的亮点也是最容易出Bug的地方。我分三步写第一步计算系统的总碳排放量。燃气设备的碳排放量等于气耗量乘以排放因子购电的碳排放量等于购电量乘以电网等效排放因子total_carbon sum( (GT_gas[t] GB_gas[t]) * CO2_factor_gas P_buy[t] * CO2_factor_grid for t in range(T) )第二步计算免费配额。配额系数按设备区分CHP由于热电联产效率高配额系数给得也相对宽松carbon_quota sum( P_gt[t] * quota_gt H_gt[t] * quota_gt_h H_gb[t] * quota_gb for t in range(T) )第三步计算需要购买的配额量并分段计算阶梯碳价。在使用Gurobi的情况下我推荐用addGenConstrPWL但那个写法是给表达式用的不能直接用于自定义变量。所以更通用的做法是创建一个新的连续变量carbon_buy表示需要购买的碳配额然后把total_carbon减去carbon_quota赋值给它再用分段函数计算成本carbon_buy m.addVar(lb0, ubGRB.INFINITY, namecarbon_buy) m.addConstr(carbon_buy total_carbon - carbon_quota, namecarbon_balance) # 阶梯碳交易分段线性函数 carbon_cost_var m.addVar(lb0, ubGRB.INFINITY, namecarbon_cost) points [0, 50, 100, 200] # 碳排放购买量分段点 values [0, 50*60, 50*6050*80, 50*6050*80100*100] # 各分段点的累计成本 m.addGenConstrPWL(carbon_buy, carbon_cost_var, points, values, namecarbon_pwl)这段代码的含义是碳购买量在0到50吨之间时按60元/吨计价50到100吨之间这50吨按80元/吨再往上按100元/吨。如果你想把碳交易成本直接加进目标函数用carbon_cost_var变量就行。把carbon_cost_var加进目标函数的时候要和燃料成本、购电成本、运维成本放一起Gurobi会自动处理这个分段线性函数的线性化。这个addGenConstrPWL是Gurobi非常强大的功能内部会做MIP线性化不需要你自己去写一堆0-1变量非常推荐。3.5 IGDT鲁棒模型的迭代求解框架IGDT的核心迭代逻辑我单独写成函数。先跑一次确定性模型拿到基准成本C0。然后写一个固定beta下的最大α搜索循环def solve_igdt_robust(beta, alpha_max0.5, alpha_step0.02): C0, _ solve_deterministic() threshold_cost (1 beta) * C0 alpha_feasible 0.0 alpha alpha_step while alpha alpha_max: status, cost solve_robust_at_alpha(alpha, threshold_cost) if status optimal: alpha_feasible alpha alpha alpha_step else: break return alpha_feasiblesolve_robust_at_alpha这个函数做的事情就是把确定性模型里光伏、风电、负荷的实际值用受α影响的最坏情景值替换再加一条总成本不超过threshold_cost的约束然后重新求解。调用Gurobi求解时如果模型可行且最优说明在当前α下系统还能承受如果返回Infeasible说明α太大系统扛不住了。注意这里有个小技巧虽然名义上叫“最大化α”但实际迭代法里是在逐步试探可行的α上界。如果你嫌迭代速度慢可以用二分法搜索先设alpha_low等于0alpha_high等于1每次取中点求解根据可行性收缩区间。理论上二分法效率更高但迭代法画出来的那个阶梯状曲线更直观。我两个都试过实际工程里参数没那么多迭代法多跑几次也就几十次求解Gurobi单次求解也就几秒钟完全能接受。3.6 结果可视化与对比分析模型跑完之后光看数值没感觉画图才是最能说明问题的。我一般画四张图。第一张是多能系统的电功率平衡图用堆叠面积图展示每个时段光伏、风电、CHP、购电、储能放电分别承担了多少电负荷和电锅炉消耗能非常直观地看出不确定情景下系统的电源结构变化。第二张是热功率平衡图展示燃气轮机余热、燃气锅炉、电锅炉、蓄热罐在热负荷侧的分配。第三张是碳交易成本的阶梯变化图横轴是碳购买量纵轴是碳交易费用把阶梯碳价的台阶画出来然后在图上标出确定性方案和IGDT方案对应的点一眼就能看出IGDT方案因为不确定性裕度导致的碳排放增加和费用增加。第四张是风险评估曲线横轴是风险偏好系数β纵轴是系统能承受的最大不确定性半径α。这条曲线是整个模型的核心产出相当于一台机器的“抗风险能力说明书”。当β从0.05增加到0.3时α会随之增大但增速会逐渐放缓说明边际收益递减。决策者看这张图做决策清楚得很。画图代码用matplotlib即可注意中文字体问题在Windows上要设置SimHei在Linux服务器上要设置中文字体路径否则标签全部乱码import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei] plt.rcParams[axes.unicode_minus] False4. 实操中的常见问题与排查技巧4.1 模型求解速度慢、一直转圈的排查方案多能系统优化调度模型规模并不大24个时段几十个变量Gurobi求解一般都在几秒内。如果你发现求解时间异常长多半是模型本身出现了一些“病态”结构。最常见的是储能约束里0-1变量和连续变量之间的Big-M取值太大导致求解器在分支定界时搜索空间巨大。解决办法是把Big-M紧到和设备最大出力一致不要随便给个一万。比如蓄电池最大充放电功率是20MW那约束里的M就写20写1000就是给自己挖坑。还有一种情况是阶梯碳交易分段数太多或者分段点设置不合理。我曾经试过把碳交易分成五个阶梯每个阶梯宽度只有10吨结果模型求解时间从两秒涨到四十秒。后来发现阶梯分得太细MIP的凸包松弛效果变差分支数量暴增。工程上三个阶梯完全够用甚至两个阶梯都能也能表达核心思想。另外Gurobi的MIPGap参数也是提速利器。学术论文追求精确解默认的MIPGap是1e-4但工程场景下设置成0.01完全够用就是把最优性 gap 放宽到1%求解速度能提升好几倍。在代码里加一行m.setParam(MIPGap, 0.01)4.2 IGDT模型求解出来的方案过于保守怎么办IGDT模型因为取的是最坏情景调度结果必然比确定性模型更保守表现出来就是CHP和锅炉出力增大、购电量增加、储能留更多裕量总成本上升。有些刚接触IGDT的同学会觉得方案太不经济了怀疑是不是代码写错了。首先检查不确定参数的方向对不对。风险规避模型下光伏和风电取(1-α)倍负荷取(1α)倍这三者都是往让系统更紧张的方向偏离。如果你弄反了把光伏取了(1α)倍那你求出来的所谓“最坏情景”其实是最好情景α值会虚高模型完全失真。其次β的选择要结合实际需求。β设得越大模型预留的抗风险裕度就越大但方案也就越保守。如果做的是日前调度光伏预测的准确度相对高β取0.05到0.1就够如果做的是实时调度且天气变化剧烈β取0.2甚至0.3也合理。还有一个常见的认知误区鲁棒模型的α并不是越大越好。α大说明系统能扛的极端情况多但代价是正常情况下的运行成本更高。IGDT的输出不是让你无脑选α最大的方案而是给你一条α随β变化的曲线让决策者根据风险偏好自己拍板。这点在论文里务必写清楚否则评审会质疑你的方案缺乏经济性。4.3 碳排放数据对不上免费配额和实际排放算错这个坑我印象很深。早期版本我把燃气轮机的碳排放直接按热值效率折算结果实际排放量算出来高得离谱后来发现是把天然气热值单位搞混了。天然气的热值大约是35.6 MJ/m3折合成能量单位就是9.89 kWh/m3约等于0.00989 MWh/m3。如果你直接用立方米数和排放因子相乘必须先统一单位。再有就是免费配额的计算口径。不同文献给出的配额系数差异很大有的按设备额定容量给有的按实际发电量给。我推荐按实际发电量给因为这样碳排放强度低的CHP在实际运行时能获得更多配额相当于变相奖励。如果你的模型结果是碳交易成本为负数即卖配额获利也不是不可能但你要检查一下是不是配额系数给得太宽松导致模型疯狂发电去薅碳配额的羊毛这种方案虽然成本低但失去了物理合理性和政策约束的本意。4.4 不同beta下结果对比结果不单调需要检查哪些地方理论上β越大可承受的α应当越大也就是曲线应当是单调递增的。如果你画出来的曲线出现回退也就是β从0.1增加到0.15α反而变小了那基本可以断定是求解精度问题。我遇到过的一次情况是固定α下求解时MIPGap设置太大导致每次求解出来的成本是近似最优判断是否超过threshold_cost时出现了误判。解决方案是IGDT内层求解时把MIPGap收紧到1e-3外层迭代再适当放宽。另外在某些α值下模型刚好处于可行与不可行的边界附近数值计算误差会引起抖动这种情况可以通过在判断约束里加一个很小的松弛来处理比如成本允许超过threshold_cost的0.1%以内这也符合工程上对决策偏差的容忍度。4.5 模拟结果画图时中文字体乱码和坐标轴不美观matplotlib中文字体乱码是老生常谈的问题我建议在项目开始时就创建好一个style配置文件统一设置字体、字号、线宽、颜色否则后面论文配图会改到崩溃。坐标轴的标签要写清楚单位功率单位写成MW时间单位写成h成本单位写成元碳排量单位写成吨这些细节直接决定你的图能不能被评审一眼看懂。还有一个小技巧画堆叠面积图时要把储能放电和充电的净功率分开处理否则会出现某一时段储能的净功率导致面积图出现负值面积图就不会正确堆叠。我在代码里是把P_bess_d和P_bess_c分别放在图的两侧充电作为负荷的一部分向下画放电作为电源的一部分向上画这样视觉效果最清晰。5. 模型扩展方向与工程落地建议5.1 从园区系统扩展到区域多能系统我上面这套模型处理一个园区的电热耦合是够用的。但如果你要扩展到一个区域加入天然气网、氢能、碳捕集等环节模型规模会呈指数级增长。这时有几个扩展方向一是加入P2G设备让富余风电转成天然气储存起来进一步消纳可再生能源但这个设备模型会引入电-气双向耦合约束复杂度提升不少二是加入碳捕集与封存装置让CHP的碳排放量变成可控变量这时候碳交易模型的逻辑会从“被动买配额”变成“主动选择捕集还是买配额”优化空间更大。这两个方向都是目前的热门研究方向代码框架不用改加设备模型和对应约束即可。5.2 从IGDT扩展到其他不确定性处理方法的对比IGDT的好处是不需要概率分布但它的缺陷也在这里——它等概率看待区间内所有情景不考虑真实概率密度。如果系统所在地的天气数据积累充分你能准确知道光伏出力预测误差的概率分布那随机优化方法的成本会更低。我在实际项目中做过的做法是用IGDT得到一套保守的基准方案再用蒙特卡洛模拟做后验评估发现IGDT方案在90%以上的场景下成本都控制在预算内而传统确定性方案只有不到50%的场景满足要求。这种“不确定性方法加蒙特卡洛验证”的组合方案既具备学术上的说服力又在工程上有实践价值。5.3 代码工程化过程中的几个建议最后补充几个代码工程化的建议。第一所有的参数尽量外部化不要写死在模型代码里用Excel或者CSV文件维护这样换一个园区案例只需要改数据文件不用动代码逻辑。第二模型求解结果一定要做可行性校验尤其是能量平衡约束脚本里加一段断言检查for t in range(T): lhs P_gt[t].X P_wind[t].X P_pv[t].X P_buy[t].X P_bess_d[t].X rhs P_load[t].X H_eb[t].X / EB_eta P_bess_c[t].X assert abs(lhs - rhs) 1e-4, ft{t} 电功率不平衡偏差{abs(lhs-rhs):.5f}MW第三多跑几组参数做敏感性分析把β、碳价、天然气价格都扫一遍画出敏感性热力图这个在写报告时非常有用能直接告诉决策者哪个因素的变动对成本影响最大。我在实际项目中还碰到过一个有意思的情况某个园区用户看了分析报告后问既然IGDT能保证最坏情况下的成本不超标能不能反过来用让运行人员在负荷和光伏比较稳定的时段t切换回确定性调度模式只在预测误差大的时段启用鲁棒方案。这个思路其实非常合理相当于做了一个“自适应的不确定性裕度”把IGDT从整体调度策略粒度细化到时段粒度。代码上实现也不难就是把原来的单一α改成按时段变量α_t再对不同时段的α分别设约束。这个扩展我给好几个朋友都推荐过效果很好发小论文或者做毕业设计都是加分项。做优化调度的项目最怕的就是模型建得花团锦簇代码一跑就报Infeasible。遇到这种情况千万别慌先检查单位再查功率平衡约束然后逐条注释约束看哪条出了问题。多能系统模型因为设备多、耦合多排查过程确实繁琐一些但一旦跑通后面扩展其他设备、换其他场景都非常方便。这套代码我前前后后迭代了四个版本从最初只有电热两个环节到后来加了碳交易、IGDT、P2G和碳捕集核心框架一直没变过。你照着这篇文章把基础版搭起来后续往哪个方向加设备、加约束、加政策机制都不会觉得无从下手。 基于信息间隙决策理论的多能系统阶梯碳交易优化调度Python实现
做综合能源系统优化调度的朋友这两年应该没少在论文里碰到“信息间隙决策理论IGDT”和“阶梯碳交易”这两个词。尤其是做双碳方向、园区多能互补、以及电热气氢耦合系统研究的基本绕不开这套组合。这个题目“基于信息间隙决策理论的多能系统-阶梯碳交易优化调度Python代码实现”本质上就是在解决一个问题风光出力不确定、碳价政策又在变这时候你怎么让多能系统既省钱又稳得住还能算得出来。这篇文章就把这套东西从头到尾拆开讲清楚。我会把IGDT的核心思想、多能系统的建模套路、阶梯碳交易的数学处理以及Python代码怎么一步步实现整个优化调度模型全部梳理出来。不管你是刚入门的硕士生还是已经在做园区综合能源项目的工程师只要按着这篇文章的思路走一遍都能自己把代码跑起来并且理解每个环节为什么要这么写。1. 整体思路拆解为什么是IGDT加阶梯碳交易1.1 多能系统优化调度的本质是什么先聊清楚多能系统到底在算什么账。所谓多能系统简单说就是把电、气、热、冷这几类能源放在一个平台里统一调度常见的设备包括燃气轮机、燃气锅炉、余热回收、电锅炉、储能、光伏、风电有时候还带P2G电转气和氢能。系统内部通过能源转换设备把一种能源变成另一种比如天然气进去出来的是电和热光伏发出来的电既可以供给电负荷也可以驱动电锅炉制热。这种系统的优化调度核心目标就是在满足各类负荷需求的前提下决定每个时段每台设备的出力使得总成本最低。这里面有几个关键约束能量平衡约束电、热、气各自守恒、设备出力上下限约束、爬坡约束、储能充放能约束以及碳排放相关约束。写成一个标准的数学优化问题就是一个带连续变量和0-1变量的混合整数线性规划MILP问题。但问题是现实场景里光伏和风电的出力不是精确知道的值负荷预测也是存在误差的。传统做法是假设这些不确定量服从某种概率分布然后用随机优化或者场景法处理但概率分布本身往往很难准确获得。这里IGDT的优势就体现出来了。1.2 信息间隙决策理论到底解决什么问题信息间隙决策理论英文叫Info-Gap Decision Theory第一次看到这个名字会觉得挺唬人其实思想非常朴素。它假设你手里的预测数据比如明天光伏出力预测是100MW和真实值之间可能存在误差但这个误差有多大你不知道也不需要你去猜一个具体的概率分布。你只需要定义一个不确定半径α表示预测值在多大范围内浮动比如α0.1就表示真实值落在预测值上下10%的区间内。IGDT一般分两种决策模式。一种是风险规避型就是我做调度决策时要保证即使不确定量往最坏的方向变化系统的运行成本也不会超过某个我可以接受的阈值在这个前提下系统能承受多大的不确定性也就是α能取到多大越大说明系统抗风险能力越强。另一种是风险追求型也叫机会型是反过来利用不确定性期望在实际值往好的方向偏差时能获得额外收益此时关注的是最小可能成本。把IGDT用到多能系统调度里步骤也很清晰。先求一个确定性模型的最优解得到一个基准成本C0。然后引入决策者的风险偏好系数β表示你愿意在基准成本的基础上多付出β比例的费用来换取抵御不确定性的能力。接下来构建一个鲁棒优化模型在不确定参数的区间范围内找到使得系统成本恰好不超过(1β)C0的最大不确定性半径α。这个α值越大说明系统越“抗造”。通过对比不同β下的α就能得到一组决策者可以根据自身风险态度来选择的帕累托前沿。1.3 碳交易为什么要做成阶梯式至于碳排放这块传统的碳交易模型大多采用固定碳价或者单一碳价每排放一吨碳就付一笔固定费用。这种做法的缺陷在于它对碳排放量大的主体约束力不够强排得多和排得少在边际成本上没有区别。阶梯碳交易机制就聪明在这里它像阶梯电价一样把碳排放配额按区间划分超出基准配额的排放量越多对应的碳价越高碳排放的边际成本是递增的。具体建模时先用基准线法确定系统的免费碳排放配额这个配额通常跟机组的发电量或者装机容量挂钩。然后计算系统实际碳排放量包括外购电力的等效碳排放和燃气设备燃烧产生的直接碳排放。实际排放量减去免费配额差额如果为正就要去碳市场买配额按照阶梯价格付费如果为负说明有多余的配额可以卖出获利。阶梯碳交易函数是一个分段线性函数在优化模型里需要引入0-1变量来分段处理这会让模型变成一个混合整数规划。但这部分在Python里用Gurobi或者Cplex求解并不困难关键是分段点的设置和松弛处理后面代码部分我会详细讲。1.4 三块拼起来为什么能自洽把这三件事拼在一起的逻辑其实是闭环的。多能系统因为风光不确定导致运行成本波动IGDT给出了一个量化这种风险的工具碳排放约束通过阶梯碳交易嵌入到成本函数中让碳排成本成为调度决策的一部分两者合在一起目标函数就成了“燃料成本购电成本设备运维成本碳交易成本”然后在这个目标下寻找既经济又稳健的运行方案。从学术角度看这种组合能解决单一方法解决不了的问题只有IGDT没有碳交易模型不够完整不符合双碳政策导向只有碳交易没有IGDT无法处理风光不确定性两者结合恰好覆盖了多能系统优化调度领域最核心的两个痛点。这也是为什么近两年这个方向的文章特别多你如果做研究用这个思路很容易出成果而且代码跑出来的结果会非常有说服力——对比确定性模型和IGDT模型的结果画几张图就能直观展示鲁棒性和经济性之间的权衡。2. 数学建模全流程从设备模型到目标函数2.1 典型多能系统物理架构与设备建模我假设系统中的设备包括燃气轮机CHP、燃气锅炉GB、电锅炉EB、光伏PV、风电WT、蓄电池BESS、蓄热罐TESS以及外网购电。这样一个架构覆盖了大部分园区级综合能源系统的典型配置代码写出来也方便你后续往里面加电转气、氢燃料电池、冷负荷等模块。设备建模时需要注意几个关键量的单位统一这是个特别容易踩坑的地方。功率类变量统一用MW时间尺度选1小时能量类变量用MWh成本类金额用元碳排量用吨。只要单位统一了后面写约束才不会出数值问题。燃气轮机的建模核心是热电联产关系。抽凝式CHP的可运行域是一个四边形区域电出力和热出力之间存在耦合关系通常用几个不等式来刻画# CHP可运行区域约束示例电量P_chp热量H_chp # 区域顶点单位MW由机组参数决定此处为四个约束限制 m.addConstr(P_chp 0) m.addConstr(H_chp 0) m.addConstr(P_chp - 0.6 * H_chp 20) m.addConstr(P_chp 0.4 * H_chp 120)如果你用的是背压式CHP那热电比是一个固定值模型更简单热出力就等于电出力乘上热电比。燃气锅炉就是典型的线性模型天然气输入热量转换成热出力效率一般取0.85到0.9。电锅炉也类似效率取0.95左右。储电和储热模型则采用荷电状态SOC递推公式SOC[t1] SOC[t] P_charge[t] * eta_charge - P_discharge[t] / eta_discharge这里有几个防止模型吹牛的约束必须写不然求出来的解根本没有实际意义。储能设备同一时段不能同时充放电需要用0-1变量限制初始和末尾SOC要相等这是调度周期闭环的要求尤其是做日前调度时一天结束和一天开始的储能量应当一致否则第二天没法继续。2.2 目标函数燃料购电运维碳交易四本账目标函数的设计是整个模型的灵魂。我推荐把总成本写成四项之和第一项是燃料成本主要包括天然气的消耗成本。燃气轮机的气耗量可以近似表示为其出力的线性函数# 燃料成本单位天然气价格gas_price元/MWh # CHP和燃气锅炉消耗的天然气功率MWh/h fuel_cost gas_price * sum(GT_gas[t] GB_gas[t] for t in range(T))第二项是购电成本从上级电网购电的价格按照分时电价处理峰谷平时段价格不同electricity_cost sum(price_buy[t] * P_buy[t] for t in range(T))第三项是设备运维成本一般取各设备出力的比例系数光伏和风电的运维成本很低CHP和锅炉稍高储能取充放电功率的成本。这一项虽然单项不大但加上去会让模型更真实。第四项就是阶梯碳交易成本这部分是重头戏下面单独拆解。目标函数最终写成total_cost fuel_cost electricity_cost maintenance_cost carbon_cost2.3 碳交易成本的计算逻辑与数学表达式碳排放量的计算分两部分。燃气轮机和燃气锅炉燃烧天然气产生直接碳排放排放因子取单位热值天然气的CO2排放系数大约是0.2吨CO2/MWh不同文献取值略有差异你按自己的参考标准就行。外购电力虽然不在本地直接排放但在全生命周期口径下电网端的火电占比决定了它的等效碳排放业内一般取0.58到0.997吨CO2/MWh看你所在区域的电网排放因子。免费碳排放配额的确定方式比较多常见做法是按照系统的总输出能量乘以一个配额系数。我代码里用的是按CHP电出力和热出力、锅炉热出力分别乘以各自的配额系数再加总。这样做的合理性在于多能系统里不同设备的能源利用效率不同给效率高的CHP更高的配额系数相当于政策上鼓励热电联产。计算完实际碳排放量E_actual和免费配额E_quota之后需要购买的碳排放额度E_purchase等于E_actual减去E_quota的差额。阶梯碳交易的过程就是针对这个差额分段计价。我设置三个阶梯第一阶梯是0到E1碳价为c1第二阶梯是E1到E2碳价为c2第三阶梯是超过E2的部分碳价为c3。c2比c1高c3比c2高体现碳排放边际成本递增的政策意图。在Python代码中实现这个分段计价我推荐用Gurobi自带的分段线性约束addGenConstrPWL这样最稳妥不需要手动引入0-1变量代码简洁且不容易出错。如果你用的是其他求解器那就老老实实引入0-1变量做混合整数线性化逻辑也不复杂就是把每个阶梯段的选择用二进制变量表示并加相应的约束限制该段的取值范围。2.4 IGDT鲁棒模型的构建与转换过程IGDT模型的构建需要先做一次确定性求解。以光伏预测值PV_forecast、风电预测值WT_forecast和负荷预测值Load_forecast为输入求出基准成本C0。然后定义不确定参数的实际值落在预测值的(1-α)到(1α)区间内PV_actual[t] (1 - alpha * delta_PV[t]) * PV_forecast[t] # 风险规避下取最坏情况 WT_actual[t] (1 - alpha * delta_WT[t]) * WT_forecast[t] # 同上 Load_actual[t] (1 alpha * delta_Load[t]) * Load_forecast[t] # 负荷往大偏风险规避型IGDT的逻辑是在最坏的不确定性情景下优化调度方案使得系统成本不超过(1β)C0并且在此约束下最大化α。注意这里最坏情况的方向是有讲究的。光伏和风电往小偏是坏方向因为发电量少了系统得多买电或者多烧气负荷往大偏是坏方向因为要多供应能量。所以鲁棒模型中不确定参数取极值时光伏取(1-α)倍风电取(1-α)倍负荷取(1α)倍。这样IGDT模型就变成了一个两层问题。外层最大化α内层在给定α下求最优调度。实际代码处理时可以采用迭代法从α0开始每次增加一个步长比如0.02在固定α下求解一次带鲁棒约束的优化模型看最小成本是否突破了(1β)C0。如果没突破就继续增大α如果突破了说明系统能承受的最大不确定性半径就是上一个α值。这个迭代法虽然慢一点但非常直观而且每步求解的都是一个MILPGurobi处理起来没有压力。固定α下单次求解的模型就是把确定性模型中的光伏、风电、负荷替换成受α影响的值然后加一条总成本不超过(1β)C0的约束即可。目标函数此时可以设成就地最小化成本也可以直接设为一个常数反正只是判断可行性。我习惯用最小化成本作为目标这样还能顺便拿到该α下的调度方案。3. Python代码实现从环境配置到跑通全流程3.1 编程环境与依赖库安装代码用的是Python 3.8以上版本求解器是Gurobi。Gurobi的学术许可申请很快用校园网邮箱注册就能拿到。如果你拿不到Gurobi也可以把模型改写成用Cbc或者HiGHS求解但MILP求解性能会差一些尤其是当系统规模扩大后差距会非常明显。需要安装的库只有这么几个pip install gurobipy pandas numpy matplotlibpandas用来读Excel数据文件numpy处理数值计算matplotlib画结果图。如果电脑里还没装Python直接装Anaconda发行版最省事自带大部分科学计算库Gurobi装完后把安装目录下的gurobipy路径加进去或者直接用pip安装gurobipy就行现在Gurobi 10以上的版本pip安装已经非常成熟了。3.2 数据准备与参数配置数据的组织方式是整个代码能不能跑通的关键。我习惯把所有系统参数集中放在一个配置模块里这样改起来方便不用在模型代码里到处翻找某个常数的定义。核心参数包括负荷数据、光伏预测数据、风电预测数据、分时电价这些我放在三个DataFrame里另外设备参数单独定义成字典。分时电价这块尤其要注意不同地区的峰谷时段划分不同而且冬夏可能不一样。代码里用一个列表price_buy表示24个时段的购电价格峰段时间设为1.2元/kWh平时段0.7元/kWh谷时段0.3元/kWh1元/kWh等于1000元/MWh换算关系别搞混了。IGDT的参数配置就是两个核心值决策者风险偏好系数beta和不确定性步长alpha_step。beta一般取0.05到0.3之间分别代表愿意多花5%到30%的成本来换取抗不确定性能力。alpha_step取0.02这样算出来的α精度到0.02对大多数分析已经够用了你要是想画更细致的曲线可以取0.01但求解时间会翻倍。这里我给一个示例的参数字典你运行代码前需要把这个模块替换成你自己系统的数据# 系统基本参数单位统一见注释 T 24 # 调度周期24小时 params { GT_max: 100, # 燃气轮机最大电出力 MW GT_min: 20, # 燃气轮机最小电出力 MW GT_eta_e: 0.35, # 燃气轮机发电效率 GT_eta_h: 0.45, # 燃气轮机余热回收效率 GB_max: 80, # 燃气锅炉最大热出力 MW GB_eta: 0.88, # 燃气锅炉效率 EB_max: 50, # 电锅炉最大热出力 MW EB_eta: 0.95, # 电锅炉效率 BESS_cap: 80, # 蓄电池容量 MWh BESS_pmax: 20, # 蓄电池最大充放电功率 MW BESS_eta: 0.92, # 蓄电池充放电效率 TESS_cap: 120, # 蓄热罐容量 MWh TESS_pmax: 30, # 蓄热罐最大充放热功率 MW TESS_eta: 0.90, # 蓄热罐充放热效率 gas_price: 2.8, # 天然气价格 元/m3换算后约等于元/MWh CO2_factor_gas: 0.2, # 天然气燃烧CO2排放系数 吨/MWh CO2_factor_grid: 0.58, # 电网购电等效CO2排放系数 吨/MWh carbon_base_price: 60, # 基础碳价 元/吨 carbon_step_price1: 60, # 第一阶梯碳价 元/吨 carbon_step_price2: 80, # 第二阶梯碳价 元/吨 carbon_step_price3: 100, # 第三阶梯碳价 元/吨 carbon_interval: 50, # 每个阶梯的碳排放区间宽度 吨 }3.3 确定性模型的核心代码实现确定性模型是整个代码框架的地基IGDT模型只是在这个基础上做循环修改。模型主体用gurobipy来写先创建模型对象然后创建变量再写约束和目标函数。变量的定义方式我给出关键部分import gurobipy as gp from gurobipy import GRB m gp.Model(MultiEnergySystem) # 连续变量 P_gt m.addVars(T, lb0, ub100, nameP_gt) # CHP电出力 H_gt m.addVars(T, lb0, ub150, nameH_gt) # CHP热出力 H_gb m.addVars(T, lb0, ub80, nameH_gb) # 燃气锅炉热出力 H_eb m.addVars(T, lb0, ub50, nameH_eb) # 电锅炉热出力 P_buy m.addVars(T, lb0, ub150, nameP_buy) # 购电功率 P_bess_c m.addVars(T, lb0, ub20, nameP_bess_c) # 蓄电充电 P_bess_d m.addVars(T, lb0, ub20, nameP_bess_d) # 蓄电放电 H_tess_c m.addVars(T, lb0, ub30, nameH_tess_c) # 蓄热充热 H_tess_d m.addVars(T, lb0, ub30, nameH_tess_d) # 蓄热放热 SOC_bess m.addVars(T1, lb0, ub80, nameSOC_bess) # 蓄电池SOC SOC_tess m.addVars(T1, lb0, ub120, nameSOC_tess) # 蓄热罐SOC # 二进制变量 u_gt m.addVars(T, vtypeGRB.BINARY, nameu_gt) # CHP启停 u_gb m.addVars(T, vtypeGRB.BINARY, nameu_gb) # 锅炉启停 u_bess_c m.addVars(T, vtypeGRB.BINARY, nameu_bess_c) # 蓄电充电状态 u_bess_d m.addVars(T, vtypeGRB.BINARY, nameu_bess_d) # 蓄电放电状态 u_tess_c m.addVars(T, vtypeGRB.BINARY, nameu_tess_c) u_tess_d m.addVars(T, vtypeGRB.BINARY, nameu_tess_d)电功率平衡约束是系统所有能量流的核心约束它表达的是电源出力等于电负荷加电锅炉消耗加储能净充电for t in range(T): m.addConstr( P_gt[t] P_wind_actual[t] P_pv_actual[t] P_buy[t] P_bess_d[t] P_load_actual[t] H_eb[t] / EB_eta P_bess_c[t], namefpower_balance_{t} )热功率平衡约束是热源出力等于热负荷加储热净充热for t in range(T): m.addConstr( H_gt[t] H_gb[t] H_eb[t] H_tess_d[t] H_load_actual[t] H_tess_c[t], namefheat_balance_{t} )这里有个细节要注意电锅炉消耗的是电功率产出的却是热功率要不是在电平衡约束里除以EB_eta要不是在热平衡约束里乘EB_eta方向千万别搞反了。我习惯在电平衡里除以效率这样H_eb变量的物理意义是产出的热功率和热平衡约束的写法保持一致。很多初学者第一次跑模型结果报不可行排了半天查不出问题结果就是电锅炉效率这里除法和乘法用反了。CHP的燃气消耗量和电出力、热出力之间的关系我直接用经验公式。对于抽凝式机组在运行区域内气耗量近似等于电出力除以发电效率减去一个修正项但更稳妥的做法是把CHP的燃料消耗表示成电出力和热出力的线性函数。你可以根据具体机组的性能曲线拟合我这里用简单模型for t in range(T): m.addConstr(GT_gas[t] P_gt[t] / GT_eta_e 0.15 * H_gt[t], namefGT_fuel_{t})燃气锅炉的气耗量就简单了热出力除以锅炉效率for t in range(T): m.addConstr(GB_gas[t] H_gb[t] / GB_eta, namefGB_fuel_{t})储能约束除了SOC递推还要限制同一时段不能同时充放。这个用两个二进制变量来做for t in range(T): m.addConstr(P_bess_c[t] BESS_pmax * u_bess_c[t], namefbess_c_ub_{t}) m.addConstr(P_bess_d[t] BESS_pmax * u_bess_d[t], namefbess_d_ub_{t}) m.addConstr(u_bess_c[t] u_bess_d[t] 1, namefbess_state_{t}) # 初始和末尾SOC一致 m.addConstr(SOC_bess[0] SOC_bess[T], namebess_cycle)3.4 阶梯碳交易的代码实现碳交易这部分是模型的亮点也是最容易出Bug的地方。我分三步写第一步计算系统的总碳排放量。燃气设备的碳排放量等于气耗量乘以排放因子购电的碳排放量等于购电量乘以电网等效排放因子total_carbon sum( (GT_gas[t] GB_gas[t]) * CO2_factor_gas P_buy[t] * CO2_factor_grid for t in range(T) )第二步计算免费配额。配额系数按设备区分CHP由于热电联产效率高配额系数给得也相对宽松carbon_quota sum( P_gt[t] * quota_gt H_gt[t] * quota_gt_h H_gb[t] * quota_gb for t in range(T) )第三步计算需要购买的配额量并分段计算阶梯碳价。在使用Gurobi的情况下我推荐用addGenConstrPWL但那个写法是给表达式用的不能直接用于自定义变量。所以更通用的做法是创建一个新的连续变量carbon_buy表示需要购买的碳配额然后把total_carbon减去carbon_quota赋值给它再用分段函数计算成本carbon_buy m.addVar(lb0, ubGRB.INFINITY, namecarbon_buy) m.addConstr(carbon_buy total_carbon - carbon_quota, namecarbon_balance) # 阶梯碳交易分段线性函数 carbon_cost_var m.addVar(lb0, ubGRB.INFINITY, namecarbon_cost) points [0, 50, 100, 200] # 碳排放购买量分段点 values [0, 50*60, 50*6050*80, 50*6050*80100*100] # 各分段点的累计成本 m.addGenConstrPWL(carbon_buy, carbon_cost_var, points, values, namecarbon_pwl)这段代码的含义是碳购买量在0到50吨之间时按60元/吨计价50到100吨之间这50吨按80元/吨再往上按100元/吨。如果你想把碳交易成本直接加进目标函数用carbon_cost_var变量就行。把carbon_cost_var加进目标函数的时候要和燃料成本、购电成本、运维成本放一起Gurobi会自动处理这个分段线性函数的线性化。这个addGenConstrPWL是Gurobi非常强大的功能内部会做MIP线性化不需要你自己去写一堆0-1变量非常推荐。3.5 IGDT鲁棒模型的迭代求解框架IGDT的核心迭代逻辑我单独写成函数。先跑一次确定性模型拿到基准成本C0。然后写一个固定beta下的最大α搜索循环def solve_igdt_robust(beta, alpha_max0.5, alpha_step0.02): C0, _ solve_deterministic() threshold_cost (1 beta) * C0 alpha_feasible 0.0 alpha alpha_step while alpha alpha_max: status, cost solve_robust_at_alpha(alpha, threshold_cost) if status optimal: alpha_feasible alpha alpha alpha_step else: break return alpha_feasiblesolve_robust_at_alpha这个函数做的事情就是把确定性模型里光伏、风电、负荷的实际值用受α影响的最坏情景值替换再加一条总成本不超过threshold_cost的约束然后重新求解。调用Gurobi求解时如果模型可行且最优说明在当前α下系统还能承受如果返回Infeasible说明α太大系统扛不住了。注意这里有个小技巧虽然名义上叫“最大化α”但实际迭代法里是在逐步试探可行的α上界。如果你嫌迭代速度慢可以用二分法搜索先设alpha_low等于0alpha_high等于1每次取中点求解根据可行性收缩区间。理论上二分法效率更高但迭代法画出来的那个阶梯状曲线更直观。我两个都试过实际工程里参数没那么多迭代法多跑几次也就几十次求解Gurobi单次求解也就几秒钟完全能接受。3.6 结果可视化与对比分析模型跑完之后光看数值没感觉画图才是最能说明问题的。我一般画四张图。第一张是多能系统的电功率平衡图用堆叠面积图展示每个时段光伏、风电、CHP、购电、储能放电分别承担了多少电负荷和电锅炉消耗能非常直观地看出不确定情景下系统的电源结构变化。第二张是热功率平衡图展示燃气轮机余热、燃气锅炉、电锅炉、蓄热罐在热负荷侧的分配。第三张是碳交易成本的阶梯变化图横轴是碳购买量纵轴是碳交易费用把阶梯碳价的台阶画出来然后在图上标出确定性方案和IGDT方案对应的点一眼就能看出IGDT方案因为不确定性裕度导致的碳排放增加和费用增加。第四张是风险评估曲线横轴是风险偏好系数β纵轴是系统能承受的最大不确定性半径α。这条曲线是整个模型的核心产出相当于一台机器的“抗风险能力说明书”。当β从0.05增加到0.3时α会随之增大但增速会逐渐放缓说明边际收益递减。决策者看这张图做决策清楚得很。画图代码用matplotlib即可注意中文字体问题在Windows上要设置SimHei在Linux服务器上要设置中文字体路径否则标签全部乱码import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei] plt.rcParams[axes.unicode_minus] False4. 实操中的常见问题与排查技巧4.1 模型求解速度慢、一直转圈的排查方案多能系统优化调度模型规模并不大24个时段几十个变量Gurobi求解一般都在几秒内。如果你发现求解时间异常长多半是模型本身出现了一些“病态”结构。最常见的是储能约束里0-1变量和连续变量之间的Big-M取值太大导致求解器在分支定界时搜索空间巨大。解决办法是把Big-M紧到和设备最大出力一致不要随便给个一万。比如蓄电池最大充放电功率是20MW那约束里的M就写20写1000就是给自己挖坑。还有一种情况是阶梯碳交易分段数太多或者分段点设置不合理。我曾经试过把碳交易分成五个阶梯每个阶梯宽度只有10吨结果模型求解时间从两秒涨到四十秒。后来发现阶梯分得太细MIP的凸包松弛效果变差分支数量暴增。工程上三个阶梯完全够用甚至两个阶梯都能也能表达核心思想。另外Gurobi的MIPGap参数也是提速利器。学术论文追求精确解默认的MIPGap是1e-4但工程场景下设置成0.01完全够用就是把最优性 gap 放宽到1%求解速度能提升好几倍。在代码里加一行m.setParam(MIPGap, 0.01)4.2 IGDT模型求解出来的方案过于保守怎么办IGDT模型因为取的是最坏情景调度结果必然比确定性模型更保守表现出来就是CHP和锅炉出力增大、购电量增加、储能留更多裕量总成本上升。有些刚接触IGDT的同学会觉得方案太不经济了怀疑是不是代码写错了。首先检查不确定参数的方向对不对。风险规避模型下光伏和风电取(1-α)倍负荷取(1α)倍这三者都是往让系统更紧张的方向偏离。如果你弄反了把光伏取了(1α)倍那你求出来的所谓“最坏情景”其实是最好情景α值会虚高模型完全失真。其次β的选择要结合实际需求。β设得越大模型预留的抗风险裕度就越大但方案也就越保守。如果做的是日前调度光伏预测的准确度相对高β取0.05到0.1就够如果做的是实时调度且天气变化剧烈β取0.2甚至0.3也合理。还有一个常见的认知误区鲁棒模型的α并不是越大越好。α大说明系统能扛的极端情况多但代价是正常情况下的运行成本更高。IGDT的输出不是让你无脑选α最大的方案而是给你一条α随β变化的曲线让决策者根据风险偏好自己拍板。这点在论文里务必写清楚否则评审会质疑你的方案缺乏经济性。4.3 碳排放数据对不上免费配额和实际排放算错这个坑我印象很深。早期版本我把燃气轮机的碳排放直接按热值效率折算结果实际排放量算出来高得离谱后来发现是把天然气热值单位搞混了。天然气的热值大约是35.6 MJ/m3折合成能量单位就是9.89 kWh/m3约等于0.00989 MWh/m3。如果你直接用立方米数和排放因子相乘必须先统一单位。再有就是免费配额的计算口径。不同文献给出的配额系数差异很大有的按设备额定容量给有的按实际发电量给。我推荐按实际发电量给因为这样碳排放强度低的CHP在实际运行时能获得更多配额相当于变相奖励。如果你的模型结果是碳交易成本为负数即卖配额获利也不是不可能但你要检查一下是不是配额系数给得太宽松导致模型疯狂发电去薅碳配额的羊毛这种方案虽然成本低但失去了物理合理性和政策约束的本意。4.4 不同beta下结果对比结果不单调需要检查哪些地方理论上β越大可承受的α应当越大也就是曲线应当是单调递增的。如果你画出来的曲线出现回退也就是β从0.1增加到0.15α反而变小了那基本可以断定是求解精度问题。我遇到过的一次情况是固定α下求解时MIPGap设置太大导致每次求解出来的成本是近似最优判断是否超过threshold_cost时出现了误判。解决方案是IGDT内层求解时把MIPGap收紧到1e-3外层迭代再适当放宽。另外在某些α值下模型刚好处于可行与不可行的边界附近数值计算误差会引起抖动这种情况可以通过在判断约束里加一个很小的松弛来处理比如成本允许超过threshold_cost的0.1%以内这也符合工程上对决策偏差的容忍度。4.5 模拟结果画图时中文字体乱码和坐标轴不美观matplotlib中文字体乱码是老生常谈的问题我建议在项目开始时就创建好一个style配置文件统一设置字体、字号、线宽、颜色否则后面论文配图会改到崩溃。坐标轴的标签要写清楚单位功率单位写成MW时间单位写成h成本单位写成元碳排量单位写成吨这些细节直接决定你的图能不能被评审一眼看懂。还有一个小技巧画堆叠面积图时要把储能放电和充电的净功率分开处理否则会出现某一时段储能的净功率导致面积图出现负值面积图就不会正确堆叠。我在代码里是把P_bess_d和P_bess_c分别放在图的两侧充电作为负荷的一部分向下画放电作为电源的一部分向上画这样视觉效果最清晰。5. 模型扩展方向与工程落地建议5.1 从园区系统扩展到区域多能系统我上面这套模型处理一个园区的电热耦合是够用的。但如果你要扩展到一个区域加入天然气网、氢能、碳捕集等环节模型规模会呈指数级增长。这时有几个扩展方向一是加入P2G设备让富余风电转成天然气储存起来进一步消纳可再生能源但这个设备模型会引入电-气双向耦合约束复杂度提升不少二是加入碳捕集与封存装置让CHP的碳排放量变成可控变量这时候碳交易模型的逻辑会从“被动买配额”变成“主动选择捕集还是买配额”优化空间更大。这两个方向都是目前的热门研究方向代码框架不用改加设备模型和对应约束即可。5.2 从IGDT扩展到其他不确定性处理方法的对比IGDT的好处是不需要概率分布但它的缺陷也在这里——它等概率看待区间内所有情景不考虑真实概率密度。如果系统所在地的天气数据积累充分你能准确知道光伏出力预测误差的概率分布那随机优化方法的成本会更低。我在实际项目中做过的做法是用IGDT得到一套保守的基准方案再用蒙特卡洛模拟做后验评估发现IGDT方案在90%以上的场景下成本都控制在预算内而传统确定性方案只有不到50%的场景满足要求。这种“不确定性方法加蒙特卡洛验证”的组合方案既具备学术上的说服力又在工程上有实践价值。5.3 代码工程化过程中的几个建议最后补充几个代码工程化的建议。第一所有的参数尽量外部化不要写死在模型代码里用Excel或者CSV文件维护这样换一个园区案例只需要改数据文件不用动代码逻辑。第二模型求解结果一定要做可行性校验尤其是能量平衡约束脚本里加一段断言检查for t in range(T): lhs P_gt[t].X P_wind[t].X P_pv[t].X P_buy[t].X P_bess_d[t].X rhs P_load[t].X H_eb[t].X / EB_eta P_bess_c[t].X assert abs(lhs - rhs) 1e-4, ft{t} 电功率不平衡偏差{abs(lhs-rhs):.5f}MW第三多跑几组参数做敏感性分析把β、碳价、天然气价格都扫一遍画出敏感性热力图这个在写报告时非常有用能直接告诉决策者哪个因素的变动对成本影响最大。我在实际项目中还碰到过一个有意思的情况某个园区用户看了分析报告后问既然IGDT能保证最坏情况下的成本不超标能不能反过来用让运行人员在负荷和光伏比较稳定的时段t切换回确定性调度模式只在预测误差大的时段启用鲁棒方案。这个思路其实非常合理相当于做了一个“自适应的不确定性裕度”把IGDT从整体调度策略粒度细化到时段粒度。代码上实现也不难就是把原来的单一α改成按时段变量α_t再对不同时段的α分别设约束。这个扩展我给好几个朋友都推荐过效果很好发小论文或者做毕业设计都是加分项。做优化调度的项目最怕的就是模型建得花团锦簇代码一跑就报Infeasible。遇到这种情况千万别慌先检查单位再查功率平衡约束然后逐条注释约束看哪条出了问题。多能系统模型因为设备多、耦合多排查过程确实繁琐一些但一旦跑通后面扩展其他设备、换其他场景都非常方便。这套代码我前前后后迭代了四个版本从最初只有电热两个环节到后来加了碳交易、IGDT、P2G和碳捕集核心框架一直没变过。你照着这篇文章把基础版搭起来后续往哪个方向加设备、加约束、加政策机制都不会觉得无从下手。