做综合能源系统优化调度绝大多数人绕不开MATLAB。这段时间我刚把一个基于绿证交易与综合需求响应的综合能源系统优化调度模型从建模到求解完整跑通趁着热乎把思路、公式、代码结构和调试时踩过的坑都梳理一遍。这个方向现在很热不管你是做园区级电热气联供、微电网调度还是研究生开题做IES优化这套模型都能直接用代码层面也给出了完整的YALMIP实现思路照着改数据就能出结果。这篇博文重点解决三个问题绿证交易这个政策工具怎么量化进目标函数和约束里综合需求响应电、热、气负荷的灵活性怎么建模不失真以及整个优化模型用MATLAB落地的代码架构和求解经验。我尽量把每一步为什么这么做讲明白不只是给一段能跑的代码。1. 建模之前先想清楚两件事1.1 绿证交易到底在模型里扮演什么角色绿证的全称是绿色电力证书1张绿证对应1000千瓦时可再生能源上网电量是证明可再生能源发电量的唯一凭证。国内实行可再生能源电力配额制之后承担配额义务的主体售电公司、电力用户、自备电厂等每年必须持有一定比例的非水可再生能源绿证没达标的部分就得去绿证市场买富余的可以卖。放到综合能源系统里园区或者微网运营商就是一个配额的义务主体。它自己装了风机、光伏这部分绿证计入自产发电量不够配额要求时就需要额外购买绿证如果自产绿证超过配额还能出售赚一笔。这个逻辑落到优化调度模型里就形成了两个层面的影响目标函数层面绿证购买支出是成本项绿证出售收入是收益项直接影响调度方案的经济性。约束层面配额要求是一个硬性不等式约束它把新能源出力、外购电量、绿证交易量三者耦合在一起。很多初学的人容易忽略的一点是绿证配额制和碳排放权交易不是一回事。碳交易约束的是二氧化碳排放总量绿证约束的是可再生能源电力消纳比例。两者可以同时存在但这一篇我们聚焦绿证。实际构建模型时绿色证书的购买量通常与外购电量正相关因为外购火电占比越高你需要补的绿证就越多。1.2 综合需求响应解决什么问题传统需求响应主要关注电力负荷用户在电价高的时候削峰填谷赚补偿。但综合能源系统的优势是多能互补需求响应不只可以做电的还能做热的、气的。这就是“综合需求响应”Integrated Demand ResponseIDR的核心思路。举一个具体的场景某园区有电锅炉和燃气锅炉电价尖峰时段如果能让电负荷减少、同时把一部分供暖需求转移到燃气锅炉侧那电网侧压力下降系统运行成本也下降。类似的可削减负荷、可转移负荷、可替代负荷这三类资源都可以拿出来调度。综合需求响应在模型里本质上是对负荷侧弹性的数学描述它把原来固定的负荷曲线变成了决策变量的一部分给了优化器更大的可行域。IDR并不是“用户随便用”它是有代价的。削减负荷要给用户补偿转移负荷可能存在舒适度成本替代负荷涉及设备转换效率损失。这些代价都需要量化进入目标函数。如果补偿成本设得太低模型会一股脑地把所有可响应负荷全砍掉结果虽然成本好看但实际不可行设得太高则响应量小起不到调节作用。这是一个需要根据实际场景反复标定的参数。1.3 系统的整体架构与数据流我用的系统架构是一个典型的电-热-气耦合综合能源系统包含风力发电、光伏发电、热电联产机组CHP、燃气锅炉、电锅炉、电储能和蓄热罐外部与配电网、天然气网相连。负荷侧包含电负荷、热负荷和气负荷其中部分电、热负荷具备需求响应能力。整个模型的数据流是这样的系统输入光伏和风电的预测出力曲线、各类负荷预测曲线、分时电价、天然气价格、绿证交易价格、配额比例、设备参数以及需求响应补偿价格决策变量是各类设备的出力、储能充放功率、与外网交互功率、绿证购买量、负荷响应量输出是最优运行成本、各设备24小时的出力计划、购电购气计划、绿证交易计划和负荷响应计划。调度周期我设置为一天24小时时间尺度为1小时。这个时间粒度对大多数园区级综合能源系统的分析够用了。如果你需要更精细的控制可以改成15分钟粒度但要注意求解规模会成倍增长MILP的求解时间可能从几十秒涨到几分钟甚至更久。2. 数学模型目标函数与核心约束2.1 目标函数搭建模型的目标函数是系统一天的运行总成本最小化包含外购电力成本、天然气购置成本、绿证交易成本可以是支出也可以是收入、需求响应补偿成本、设备运维成本。公式可以写为min F C_elec C_gas C_gc C_dr C_om各项展开来看外购电成本C_elec sum_t (price_elec(t) * P_buy(t) * dt)其中 P_buy(t) 是t时段从电网购电的功率price_elec(t) 是分时电价。一般不考虑售电如果光伏出力非常大允许倒送电的话可以加一项售电收益但要在约束里限制倒送功率否则模型会为了卖电收益而虚增出力。天然气成本C_gas sum_t (gas_price(t) * V_gas(t) * dt)V_gas(t) 是t时段消耗的天然气体积CHP和燃气锅炉各自消耗的天然气要分开计算再求和。绿证交易成本C_gc gc_price * N_buy - gc_price * N_sellN_buy是购买的绿证数量N_sell是出售的绿证数量。需要注意的是同一时段模型不应该既买又卖绿证这需要额外的约束或者巧妙设置价格参数来避免套利。需求响应补偿成本C_dr sum_t (DR_price_e * E_cut(t) DR_price_h * Q_cut(t) ...)不同类型负荷的响应补偿单价可以相同也可以不同。运维成本C_om sum( k_c * P_chp(t) k_b * P_eb(t) ... )运维成本系数一般按启停或按出力比例给一个很小的常数。2.2 CHP与多能设备的运行约束CHP机组是电热耦合的关键。最简单的模型是定热电比模型即热出力Q_chp ratio * P_chp电出力P_chp在上下限范围内。这种模型优点是线性、简洁适合做MILP缺点是不够精确因为它忽略了CHP实际运行的可行域形状。更精确的做法是用可行域约束P_chp_min ≤ P_chp(t) ≤ P_chp_max Q_chp_min ≤ Q_chp(t) ≤ Q_chp_maxP_chp(t) alpha * Q_chp(t) ≥ M_chp_min这个约束定义了电热出力的组合可行域形式可以根据机组类型调整电锅炉的约束很简单Q_eb(t) eta_eb * P_eb(t)且 0 ≤ P_eb(t) ≤ P_eb_max。电锅炉的作用就是把便宜的电比如夜间谷电转化为热能是电热耦合的又一个支点。燃气锅炉约束0 ≤ Q_gb(t) ≤ Q_gb_max消耗天然气量 V_gb(t) Q_gb(t) / (eta_gb * LHV)。LHV是天然气低位热值一般取9.7 kWh/Nm³也就是每立方米天然气完全燃烧产生的低位热量。电储能模型是经典的SOC荷电状态递推SOC(t) SOC(t-1) (P_ch(t) * eta_ch - P_dis(t) / eta_dis) * dt / E_cap同时满足充放电功率限制0 ≤ P_ch(t) ≤ P_ch_max * u_ch(t) 0 ≤ P_dis(t) ≤ P_dis_max * u_dis(t) u_ch(t) u_dis(t) ≤ 1这里的0-1变量u_ch、u_dis避免储能同时充电和放电。这个约束初看很合理但实际运行中容易导致求解速度变慢因为引入了大量二进制变量。如果系统模型已经很大可以考虑用松弛掉互斥约束的方式先求一个参考解。2.3 综合需求响应的建模方式综合需求响应我采用了分类建模的方式分电负荷和热负荷两类处理每类各考虑可削减和可转移两种响应手段。可削减负荷比较简单每个时段允许削减一定比例不能超过上限0 ≤ E_cut(t) ≤ alpha_cut * E_load(t) 0 ≤ Q_cut(t) ≤ beta_cut * Q_load(t)alpha_cut、beta_cut是削减比例上限一般取10%到20%比较合理。削减的部分乘以补偿单价进入目标函数。可转移负荷比较麻烦。它要求负荷的总量在一段时间内保持不变只是在不同时段之间搬移。比如说某工厂有一条连续生产线每天耗电1000千瓦时它可以灵活安排到凌晨低谷时运行但不能取消。建模时可以设定转移比例gamma让原有时段的负荷减少gamma*原始负荷同时在别的时段增加同样的量E_shift(t) gamma * E_load_shiftable(t) * y(t) sum_t E_shift(t) 0这里的y(t)可以直接作为连续变量处理因为转移总量守恒约束已经保证了不会出现无中生有的情况。热负荷的转移通常和建筑热惯性相关联更精细的建模要用到建筑物的热动态方程但作为规划级调度模型用上述简化手法足够。还有一个容易被忽略的约束是削减量和转移量叠加之后不能让剩余负荷变成负数。即E_load(t) - E_cut(t) E_shift(t) ≥ 0否则优化器会钻空子把响应量设得超过实际负荷结果莫名其妙多拿补偿。2.4 绿证配额与交易约束绿证的配额约束是整个模型里最容易写错的部分。基本原则是系统当年或当日持有的绿证数量不得低于配额要求。由于我们是日调度模型按日折算配额。绿证持有量 自产绿证数量 购置绿证数量 - 出售绿证数量自产绿证数量由风电和光伏的实际出力决定N_self sum_t (P_wind(t) P_pv(t)) * dt / 1000这里dt是小时数功率单位kW转换为绿证时要除以1000得到MWh。配额要求的计算公式为N_quota quota_rate * total_supplytotal_supply是总供电量包括外购电和新能源发电。这就形成了内生耦合外购电越多需要持有的绿证也越多。如果系统内新能源占比不足就得去买绿证补足缺口如果新能源出力很充足超过配额的部分可以卖出。建模时要注意一个陷阱直接用max(0, N_quota - N_self)这种表达是非线性、非凸的CPLEX和Gurobi都不吃这一套。正确的做法是引入两个非负变量N_buy和N_sell加上约束N_self N_buy - N_sell ≥ N_quota因为目标函数里买绿证是成本、卖绿证是收益模型自然会选择在满足该约束的前提下尽量少买多卖。但不排除某种极端价格情况下模型既买又卖进行套利所以需要再加一条逻辑约束。这可以用Big-M方法N_buy ≤ M * z N_sell ≤ M * (1 - z)其中z是一个0-1变量M是一个足够大的常数这个约束确保同一时段只能买或者只能卖。3. MATLAB代码实现从数据到求解3.1 程序结构设计我习惯把程序拆成四个文件参数文件、数据文件、主程序和结果绘图脚本。参数文件放设备参数、价格参数、约束上下限数据文件放负荷曲线、新能源出力曲线、分时电价等时序数据主程序负责建模和求解绘图脚本输出结果。Project/ ├── setup_params.m % 设备参数与价格参数 ├── load_data.m % 读取负荷、新能源预测曲线 ├── main_scheduling.m % 主程序建模、求解、保存结果 ├── plot_results.m % 绘制调度结果图主程序的整体流程是加载参数和数据 - 定义变量 - 写目标函数 - 写约束 - 设置求解器选项 - 求解 - 提取结果 - 保存。用YALMIP的话代码层次会非常清晰。3.2 关键代码片段与实现技巧变量定义部分T是24所有功率类的决策变量都是T维连续变量P_chp sdpvar(1, T, full); % CHP电出力 Q_chp sdpvar(1, T, full); % CHP热出力 P_buy sdpvar(1, T, full); % 购电功率 P_eb sdpvar(1, T, full); % 电锅炉电功率 Q_gb sdpvar(1, T, full); % 燃气锅炉热出力 SOC sdpvar(1, T, full); % 储能荷电状态 E_cut sdpvar(1, T, full); % 电负荷削减量 n_buy sdpvar(1, 1, full); % 绿证购买量 n_sell sdpvar(1, 1, full); % 绿证出售量 z binvar(1, 1, full); % 绿证买卖互斥电功率平衡约束Constraints [Constraints, P_wind(t) P_pv(t) P_chp(t) P_buy(t) P_dis(t) ... P_load(t) - E_cut(t) P_eb(t) P_ch(t)];热功率平衡约束Constraints [Constraints, Q_chp(t) Q_gb(t) Q_eb(t) Q_hdis(t) Q_load(t) - Q_cut(t) Q_hch(t)];CHP约束和电锅炉、储能、需求响应、绿证约束逐个追加。最后设置求解器和求解ops sdpsettings(solver, cplex, verbose, 1, showprogress, 1); result optimize(Constraints, Objective, ops);3.3 求解器选型与参数配置模型规模上来后求解器选型很关键。我分别用CPLEX和Gurobi跑过同一套算例速度差异不大都远快于MATLAB内置的linprog和intlinprog。YALMIP里只要写成约束和目标的集合求解器是透明切换的不需要改模型代码。针对MILP模型几个关键参数值得说明MIP gap相对最优间隙CPLEX里用cplex.mip.tolerances.mipgap设置Gurobi用OutputFlag配合MIPGap。我一般设到0.01%也就是1e-4。设得太小会白白增加求解时间特别是在有上千个整型变量的大规模算例中。时间限制设到600秒。如果600秒内没找到最优解至少能得到一个可行解用来分析趋势足够。数值容差CPLEX的EpOpt默认1e-6一般不用动。但如果你发现结果里功率值出现0.0001这种脏数据可以考虑把容差调大一点到1e-5或1e-4换一个干净的解。对称性MILP求解器非常怕对称性。如果你有多个同类型机组并联可以加简单的排序约束比如P_chp1 P_chp2能显著加速求解。这个技巧很实用但对单机组模型没有帮助。3.4 典型参数设置参考我把一套能跑出合理结果的典型参数列在下面方便你直接作为初始值使用参数数值说明风电装机200 kW预测出力曲线按最大出力的0-90%波动光伏装机150 kW白天呈钟形夜间为0CHP容量300 kW电出力上限热电比取1.3电锅炉容量100 kW效率0.95燃气锅炉容量200 kW效率0.9电储能容量200 kWh最大充放电功率50 kW蓄热罐容量300 kWh最大充放热功率80 kW绿证配额比例15%对应非水可再生能源消纳责任权重绿证价格50元/张每张对应1 MWh可再生能源电量需求响应补偿0.4元/kWh略高于分时电价谷段低于峰段这套参数下模型一天的总运行成本大约在1.2万到1.8万元量级具体取决于外购电价格和新能源出力情况。绿证成本约占总成本的5%到12%需求响应能降低3%到8%的总成本。4. 场景对比与结果分析4.1 实验场景设置为了说清楚绿证交易和需求响应各自的作用我设计了四个对比场景场景编号绿证配额约束综合需求响应说明S1不启用不启用基准场景只做常规经济调度S2启用不启用只看绿证配额的影响S3不启用启用只看需求响应的影响S4启用启用绿证 需求响应完整模型4.2 结果怎么看成本构成与出力曲线完整的模型跑完之后结果文件会包含各时段出力、储能SOC曲线、购电曲线、绿证交易量、各成本项明细。我建议重点看三张图各设备24小时出力堆叠图、电负荷平衡图负荷基线、响应后负荷、各电源供给关系、成本构成饼图。根据我的实际运行数据几个典型结论如下开启绿证配额后S2相比S1的风电、光伏消纳量会提升因为新能源出力不足时需要购买绿证而自产绿证可以抵配额等于变相给了新能源更高的“影子价格”。具体数值上如果配额比例从0调到15%弃风弃光量下降约40%外购电量下降约10%到15%。单独启需求响应的S3效果最直接的是系统购电曲线的峰谷差变小。基线场景的峰谷差大约是180kW加入需求响应后降到140kW左右因为一部分高峰负荷被削减或转移到了凌晨低谷时段。电价峰段购电量明显减少谷段购电量略有增加。S4完整模型在成本上最低相比S1能降低约7%的运行成本。这里有个微妙之处绿证机制和需求响应不是简单的相加关系。绿证提高了新能源的价值需求响应又增加了系统消纳新能源的灵活性两者配合在一起新能源出力占比可以从42%提升到53%绿证购买量也随之减小。绿证价格对调度结果影响很大。我做了一组敏感性分析绿证价格分别取20、50、80元/张时系统策略有明显变化。价格低的时候绿证支出几乎可以忽略系统倾向于多买少生产价格高的时候模型会尽量多出力、多卖绿证甚至出现为了获取更多绿证而在午间光伏大发时反而降低储能充电功率的情况。这个结果初看反直觉但细想是对的储能充电再放电会有损耗如果绿证价格足够高直接卖绿证的收益超过了储能套利的收益。还有一个有意思的现象需求响应的价值在绿证价格越高时越明显。这是因为高绿证价格下削减负荷带来的“机会成本”更大——少用电意味着减少外购电力外购电少了需要配的绿证也少了双重的经济激励会推动模型更积极地去调用负荷侧弹性。5. 调试过程中踩过的坑与排查技巧5.1 求解不可行的排查路径刚写完模型的人最常见的问题是求解器直接报infeasible。我的排查顺序是固定的先查功率平衡约束。很多时候是各设备出力上限加总小于负荷峰值这种不可行是物理意义上的不可行模型再聪明也没用。我的排查手段是把每时段所有电源的出力上限求和与负荷曲线对比看看是否存在负的裕度。再查储能SOC约束。储能初始SOC设得不对也会导致不可行。比如把SOC(1)初始值设为0同时约束SOC(t) ≥ 0.1那前几个时段如果有放电需求就必然不可行。这个问题的排查方法是先把SOC约束松弛掉再跑如果可解说明问题出在SOC相关约束。最后查绿证配额约束。配额比例太高而新能源装机太小系统无论如何都买不起足够的绿证这时候模型也会不可行。这种情况不应该用调大Big-M糊弄而应该思考是不是参数本就不合理比如配额超出了系统物理能力。一个实用的技巧把模型里的整数约束全部注释掉跑一个LP松弛解。LP不可行说明约束本身有冲突LP可行但MILP不可行说明整数约束导致了问题这时候重点检查Big-M约束和逻辑充分条件类约束。5.2 非线性约束的处理细节绿证约束里max(0, quota - self)是最典型的非线性表达。还有几个常见的非线性项也要小心效率与出力成非线性关系如CHP热效率随出力变化曲线这种可以用数据包络法做分段线性化。把出力区间分成5到7段每段用线性系数近似精度已经非常好。天然气流量与压力关系的约束这是典型的非凸约束在日调度模型里建议直接按等压简化处理即认为气网恒压只考虑气量平衡。要精确到管网压力的话需要另建天然气网络潮流模型代码复杂度会翻倍。需求响应成本如果是阶梯价格前50kW削1元/kWh超出部分削2元/kWh需要引入多个0-1变量分段表示。这个分段线性化可以通过YALMIP的implies命令方便地表达但求解速度会受影响。5.3 结果不合理的常见原因求解返回成功但结果一看就离谱这种情况在调试中更常见。我总结出三类高频问题第一类是“储能一直满充满放”。如果储能几乎每个周期都是0到满载来回切换大概率是没设合理的充放电次数限制或者没有维护成本项。储能频繁深度充放会严重缩短寿命模型不知道这一点所以它只盯着套利收益。解决办法是在目标函数里加入储能充放电的损耗惩罚项比如每次充放一个周期加0.02元/kWh的折算损耗费用效果立竿见影。第二类是“需求响应把某个时段的负荷几乎切光”。这说明削减比例上限或者补偿成本有问题。检查E_cut上限约束是否生效如果约束写成了E_cut(t) alpha * P_load(t)而P_load(t)取的是没响应前的基线这个约束就有问题——因为模型可能同时用了可转移负荷导致当前实际负荷已经变小再去削减就会击穿下限。必须在约束右侧使用响应操作发生之前的原始负荷值而不是决策之后的负荷。第三类是“绿证购买量为负”。如果你忘了把n_buy和n_sell定义为非负变量或者用了sdpvar但没加≥0约束模型就会通过“负购买”来白嫖绿证收益。这个问题非常隐蔽因为YALMIP默认是自由变量不会自动加非负约束。排查绿证结果时第一步就是检查这两个变量的符号。5.4 运行时间与内存优化模型规模一大求解时间就开始放飞。我的经验是先跑LP松弛看目标函数值和约束状态再用MILP跑。如果LP松弛解和MILP解的差距很小其实可以直接采用LP结果做分析省下大把时间。整套模型用24时段、单台CHP、单台储能、两三种需求响应大约有300到500个连续变量、几十个0-1变量CPLEX求解时间一般在半分钟到两分钟内完全可以接受。但如果扩展到多台机组、多能源枢纽或者96时段整型变量会急剧增加求解时间可能直接爆炸。两个有效的提速技巧一是给所有二进制变量设置合理的初值。YALMIP里可以用assign(z, 1)这种形式给变量一个初始可行解。CPU时间在复杂算例里能节省30%以上。二是给MILP设置求解策略。Gurobi的MIPFocus设为2CPLEX的MIPEmphasis设为1重点找可行解可以让模型更快返回一个质量较好的解。6. 从模型到论文或工程应用的一些体会这个模型我前前后后改了三版从最初的纯电力调度一步步扩展到电热耦合、加绿证、加需求响应。个人最大的体会是综合能源系统优化调度最花时间的不是求解而是把工程问题翻译成数学问题的过程。具体来说绿证交易机制的政策条文翻译成约束条件时要反复核对交易规则需求响应的负荷侧数据获取比想象中困难得多——你说用户负荷可削减10%这个10%从哪来不是拍脑袋定的要依据历史负荷曲线的波动特性、用户问卷意愿、以及实际控制系统能力三者综合确定。从应用角度看这套模型即时性很强。现在很多园区的实际运营确实需要同时考虑绿证配额和负荷侧调节缺一不可。绿证配额让系统有动力最大化消纳新能源需求响应则给了系统实际消纳新能源的手段。两者联合优化的价值大于单独考虑任何一个机制。做这个模型时还有一个细节值得提醒代码注释和命名一定要规范。优化调度模型的变量动辄几十个Power、P、E、Q这几个符号在不同论文里含义经常打架。我在代码里统一用P代表电功率Q代表热功率所有变量名都带单位后缀比如P_chp_kw这样过两个月再翻代码不至于看不懂。目前这个模型还可以往几个方向扩展加入碳排放配额和碳交易机制、考虑新能源出力和负荷预测的不确定性用场景法或者鲁棒优化的方式、扩展到多个园区之间的协同调度、加装氢能系统做电氢热耦合。每一层扩展都会带来新的建模挑战和更丰富的结论但也意味着复杂度再上一个台阶。建议想上手的朋友先把基础版调试通理解每一个约束在做什么再逐层加复杂度。最后分享一个实用的小建议调试绿证约束时把绿证价格设得极高或者极低看模型的买卖行为是否符合直觉这是验证逻辑正确性最快速的方法。比如价格极高时模型应该尽量多生产绿证来卖价格极低时模型应该直接放弃超配额生产转而买证。如果结果不符合这个逻辑说明约束里一定有哪里写错了。这个思路不只适合绿证模型任何带有交易变量的调度模型都可以用类似的极限参数法来验证。