做综合能源系统调度的时候很多人会把碳交易和需求响应当成两个孤立的功能模块碳交易嘛就是在目标函数里加个碳价乘排放量需求响应嘛就是让负荷曲线削个峰填个谷。但实际上这两个机制一旦同时进入优化模型它们会通过设备出力、购电策略、储能充放电互相影响甚至可能改变整个系统的运行方式。这篇文章就来拆解我是怎么在Matlab里把碳交易和需求响应同时放进综合能源系统优化调度模型并且真正跑出可解释的结果。这篇内容适合几类人正在做综合能源、微电网、园区多能互补相关毕设或课题的同学需要复现碳交易需求响应优化调度代码的工程师以及想了解从数学建模到Matlab求解全流程的入门者。我会把碳配额的核算方式、阶梯碳价的线性化处理、价格型与激励型需求响应的建模思路以及YalmipCplex/Gurobi这套求解组合的关键代码都过一遍。文末还会聊聊我实际调试中踩过的一些坑——这些东西光看论文是学不到的。1. 为什么要把碳交易和需求响应放进同一个优化模型1.1 传统经济调度只算用能成本不算碳排放成本过去做综合能源系统调度目标函数通常只有三类成本购能成本买电、买气、设备运行维护成本以及可能的弃风弃光惩罚。这个框架本身没问题但在双碳背景下它有一个明显的盲区——碳排放变成了一个外部性系统不会主动去约束它。举个例子如果天然气价格很低、而电网购电价格相对较高传统调度模型会让燃气轮机CHP满发因为燃气发电的单位成本更低。但如果把碳交易成本算进来天然气燃烧的直接排放和电网购电的间接排放都要承担碳成本这时候CHP的低价优势可能就被碳价抵消了。更复杂的是碳配额的计算方式会影响这个判断——如果系统获得了免费配额只要实际排放低于配额甚至还能通过出售富余配额赚钱。这就不是简单加个碳价的问题了。1.2 碳交易与需求响应如何互相影响调度结果需求响应进入模型后系统端的负荷曲线不再是固定的刚性需求而是可调节的。这里有个容易忽略的联动关系价格型DR通过分时电价引导用户把峰时负荷挪到谷时系统侧的购电曲线随之变化。谷时负荷增加意味着原本在谷时被压低的CHP出力可以适当上调或者储能可以在谷时充入更多电能。碳交易机制会改变设备的边际成本排序。碳价高的时候高排放设备燃气锅炉、没装CCUS的CHP的出力会被压降低排放设备电锅炉、热泵的出力会上升。如果DR把负荷转移到了这些低排放设备的可用时段碳排放和运行成本就能同时降低。储能的作用也会被重新定义。传统调度里储能主要是套利——低充高放引入碳交易后储能的充放电策略还会影响购电量和自发电量的配比进而影响间接排放和碳配额盈亏。所以两件事必须放到同一个优化模型里联立求解分开做等价于先定负荷曲线再定设备出力而实际上负荷和设备出力是同时被决策出来的。2. 碳交易机制建模配额核算、差额结算与阶梯碳价线性化2.1 免费配额怎么算基准线法下的排放差额碳交易建模的第一步是确定配额。国内目前电力行业用基准线法比较多综合能源系统里常见做法是对供电量和供热量分别给定单位配额系数E_quota δ_e × ΣP_e,t δ_h × ΣH_t其中δ_e是单位电量的免费配额单位tCO2/MWhδ_h是单位热量的免费配额单位tCO2/GJ。这个系数的取值通常参考行业基准比如供电配额0.45~0.7 tCO2/MWh供热配额0.2~0.3 tCO2/GJ。如果是做算例设计也可以自己设定但要注意合理性——配额太紧会导致碳成本占比过高太松则碳交易形同虚设。实际排放量这边主要包含两个来源一个是天然气消耗的直接排放一个是外购电力的间接排放E_emission Σ F_gas,t × EF_gas Σ P_buy,t × EF_grid天然气排放因子EF_gas按低位热值折算大约是0.2 kgCO2/kWh实际项目里用0.185~0.21都正常电网排放因子EF_grid国内不同区域差异很大从0.3到0.7 kgCO2/kWh都有算例里取0.4~0.6比较常见。然后定义排放差额ΔE E_emission - E_quotaΔE0说明配额不足需要在碳市场购买ΔE0说明有富余配额可以出售。要注意符号方向——很多初次写代码的人在这里搞反最后目标函数里碳成本变成了负数越大越好整个调度逻辑就崩了。2.2 阶梯碳价的线性化0-1变量与大M约束实际碳市场往往采用阶梯价格目的是让排放量越高的主体边际惩罚越大。常见设定是若 ΔE ≤ 0则碳收益为 p_c × ΔE此时为负成本也就是收益其中p_c是基准碳价若 0 ΔE ≤ λ1碳成本为 p_c × ΔE若 λ1 ΔE ≤ λ2碳成本为 p_c × λ1 1.1p_c × (ΔE - λ1)若 ΔE λ2则超出部分按1.2p_c计即 p_c × λ1 1.1p_c × (λ2 - λ1) 1.2p_c × (ΔE - λ2)。这个函数是分段线性的不能在Matlab里直接用if-else写进目标函数——因为优化求解器需要的是显式的数学表达式而不是程序逻辑。分段线性函数的标准做法是引入0-1变量把定义域切分成几个区间再用大M约束把每个区间的成本和排放差额对应起来。具体实现时可以定义三个连续非负变量E1、E2、E3分别表示落在三个梯度区间内的排放差额部分以及对应的0-1变量z1、z2、z3E1代表第一档区间 [0, λ1] 内的排放量约束 E1 ≤ λ1 × z1E2代表第二档区间 (λ1, λ2] 内的排放量约束 λ1 × z2 ≤ E2 ≤ λ2 × z2E3代表第三档区间 (λ2, ∞) 内的排放量约束 E3 ≥ 0且用到时 z31。同时令 ΔE_pos E1 E2 E3且 z1 z2 z3 1只有当ΔE_pos0时才需要这些变量ΔE≤0的部分另算收益。目标函数里碳成本部分就可以写成C_carbon p_c × E1 1.1 × p_c × E2 1.2 × p_c × E3 - p_c × E_sell其中E_sell是富余配额出售量约束 E_sell ≥ -ΔEΔE为负时成立。这套东西看起来繁琐但它是把不可导的分段函数变成线性约束MILP可解形式的关键。我在Yalmip里通常用binvar定义z变量再用大M系数M取远大于物理排放量上限的数比如1e4来写互斥条件。2.3 配额制度下的一个特殊边界负差额是否允许出售有些算例会假设富余配额可以自由出售有些则不允许设一个最低持有量。这个设定会显著影响系统行为。如果允许出售系统可能因为电锅炉供热替代燃气锅炉而获得额外的碳收益如果不允许出售排放低于配额的部分就只是不亏系统的减排动力会弱一些。我的建议是算例里把两种情况都跑一遍对比结果这样论文或报告里能多一个有意义的敏感性分析。实现方式也很简单出售变量E_sell的增加一个上限约束比如不超过0就是在代码里一行的事。3. 需求响应建模价格型弹性矩阵与激励型可中断负荷3.1 价格型DR用弹性系数矩阵描述负荷转移价格型需求响应的经典模型是弹性矩阵。它的思路是用户对电价的反应可以用自弹性和交叉弹性来描述。自弹性是当前时段价格变化对当前时段负荷的影响通常为负即电价越高负荷越低交叉弹性是其他时段价格变化对本时段负荷的影响通常为正即其他时段涨价会把负荷挤到本时段。用数学式表达L_dr,t L_base,t Σ_t ε(t,t) × (π_t - π_0,t) / π_0,t × L_base,t其中π_t是优化后或者分时电价政策的价格π_0,t是基准价格。这个公式里负荷转移量和价差、基准负荷、弹性系数三者直接挂钩非常直观。不过在实际算例中弹性矩阵法有一个麻烦价格变量π本身往往也是调度模型的决策变量尤其在引入需求响应后有些模型把电价设计成实时电价这时候公式里出现π_t × L_base,t这种双变量乘积会让问题变成非线性。我的处理办法是避开来解决——把需求响应当作已知的分时电价方案比如峰谷平时段各设几档固定价格弹性矩阵计算出来的负荷变化率直接在优化前算好再代入调度模型。这样既保留了DR对负荷曲线的重塑效果又不引入非线性。如果一定要做电价内生化那就要用二元乘积线性化的手段复杂度会上升一个档次。3.2 激励型DR可中断负荷的0-1变量处理激励型DR相对更好建模。它的核心逻辑是系统在高峰时段可以呼叫用户中断一部分负荷作为交换用户获得补偿像签订一个可中断负荷合同一样。数学上每个时段定义可中断负荷量ΔL_DR,t配有0-1状态变量u_DR,t0 ≤ ΔL_DR,t ≤ ΔL_DR,max × u_DR,t∑_t u_DR,t ≤ N_max全天最多中断N次∑ ΔL_DR,t ≤ E_DR,max全天累计中断电量上限目标函数里加上补偿成本C_DR Σ c_DR × ΔL_DR,t关键是要把中断电量从该时段的电负荷里扣掉——它等于系统少供的那部分电。这个负向负荷会直接影响功率平衡约束P_supply,t L_base,t - ΔL_DR,t P_storage_net ...。我在初版代码里就吃过一次亏DR变量定义好了目标函数也加了成本但忘了更新功率平衡方程结果DR形同虚设。3.3 两类DR联合使用的协调约束价格型DR的变化量受弹性矩阵约束是连续但总量守恒的——用户转移负荷总量不减少只是挪到别的时段。激励型DR则是真减少中断的负荷就是少用了。当两者同时存在时要注意价格型DR把峰荷转移去谷时会抬高谷时负荷激励型DR在峰时中断部分负荷会进一步压低谷时需要充入的功率。如果储能策略跟不上可能出现谷时负荷依然很高、储能没空间充电、光伏被弃的尴尬结果。所以调度模型里这两类DR的变量都要参与功率平衡不能先算好价格型再叠加激励型否则会重复计算负荷削减量。4. Matlab建模与求解Yalmip与求解器选型4.1 为什么不用Matlab内置优化工具箱很多人一上来用fmincon或者linprog做调度优化遇到0-1变量就只能用intlinprog。小规模算例比如单设备、24时段intlinprog勉强能跑但综合能源系统一旦加上CHP热电耦合、储能SOC时序约束、碳交易阶梯区间变量和DR中断变量整个模型的决策变量数量和约束矩阵规模会迅速膨胀intlinprog的求解效率会变得很难看而且数值稳定性一般。我更推荐Yalmip Cplex/Gurobi这套组合。Yalmip是一个Matlab下的建模层它把优化问题用人类友好的方式写出来sdpvar、binvar、constraints、objective底层调用商业求解器。Gurobi在MILP上几乎是当前最快的学术license申请也很方便。Cplex如果拿不到新license老版本配合旧版Matlab也能用。我自己现在主力是Gurobi备胎是Cplex。4.2 核心代码骨架从变量定义到optimize调用下面给一个能跑通基本框架的MatlabYalmip片段覆盖了我们前面讨论的关键要素%% 基础数据 T 24; dt 1; % 时段长度1h load_base [...]; % 基础电负荷1x24 heat_load [...]; % 热负荷1x24 price_buy [...]; % 分时购电价 price_DR [...]; % 价格型DR套用的电价 %% 决策变量 P_chp sdpvar(1, T); H_chp sdpvar(1, T); P_gb sdpvar(1, T); % 燃气锅炉 P_eb sdpvar(1, T); % 电锅炉 P_pv sdpvar(1, T); % 光伏实际出力预测值 P_wt sdpvar(1, T); P_buy sdpvar(1, T); P_sell sdpvar(1, T); % 向电网售电 P_es_c sdpvar(1, T); % 电储能充电 P_es_d sdpvar(1, T); % 电储能放电 SOC sdpvar(1, T1); % 电储能SOC u_es binvar(1, T); % 储能充放互斥 u_chp binvar(1, T); % CHP启停 u_dr binvar(1, T); % 激励型DR中断状态 L_dr_inc sdpvar(1, T); % 激励型DR中断量 L_dr_price sdpvar(1, T); % 价格型DR变化量可为负表示负荷转移到该时段然后是约束的骨架Constraints []; % 电功率平衡 Constraints [Constraints, P_buy - P_sell P_pv P_wt P_chp ... - P_eb - P_es_c P_es_d load_base L_dr_price - L_dr_inc]; % 热功率平衡 Constraints [Constraints, H_chp P_gb P_eb heat_load]; % CHP热电耦合简化定热电比 k_chp 0.9; Constraints [Constraints, H_chp k_chp * P_chp]; Constraints [Constraints, P_chp 0.3 * P_chp_max * u_chp]; Constraints [Constraints, P_chp P_chp_max * u_chp];目标函数就比较直白了C_fuel_chp sum(P_chp / eta_chp * price_gas); % CHP耗气成本 C_fuel_gb sum(P_gb / eta_gb * price_gas); C_grid sum(price_buy .* P_buy) - sum(price_sell .* P_sell); C_om sum(om_chp * P_chp) sum(om_gb * P_gb) ...; C_carbon ...; % 按2.2节分段线性 C_dr sum(c_dr * L_dr_inc) ...; % 如果有DR补偿 Objective C_fuel_chp C_fuel_gb C_grid C_om C_carbon C_dr; ops sdpsettings(solver, gurobi, verbose, 2); Diagnostics optimize(Constraints, Objective, ops); if Diagnostics.problem ~ 0 disp(求解失败); end需要提醒的是L_dr_price是通过弹性矩阵在优化前计算好的固定值还是决策变量要分清楚。如果作为决策变量就必须额外加上弹性约束矩阵如果提前算好直接当作已知量代入平衡方程即可。两种都有人用但混用会出错。4.3 目标函数里非线性项的线性化处理这个模型里最容易出现非线性项的地方有三处一是CHP热电耦合如果是运行域多边形可能引入出力点是否在多边形内部这类约束处理方式是顶点组合法或线性不等式组近似不建议用二次约束。二是1.1节提到的碳交易阶梯价格用分段线性化。三是储能充放电的功率损耗建模。如果写成P_es_c × η_es这样的形式还好一旦要表达充电时的损耗和放电时的损耗不对称就必须用两个0-1变量互斥再加辅助变量避免出现P_es_c × P_es_d这种产品项。凡是出现两个sdpvar相乘的地方都要警惕。Yalmip会尝试自动处理部分情况但代价是模型变成非线性半定规划SDPGurobi就不认了。5. 算例设计与结果解读怎么让数字说话5.1 测试系统构成与参数设置我不建议一上来就搭一个巨型系统先把一个中型园区模型跑通更有价值。我常用的算例配置如下设备容量/参数说明CHP机组2 MW电出力热电比0.9电效率0.35主要电源燃气锅炉4 MW热效率0.9调峰热源电锅炉1 MW热效率0.95低碳热源光伏风电3 MW 2 MW可再生出力曲线给出电储能1 MW / 2 MWh效率0.95SOC范围0.1~0.9热储能2 MWh热力时移电价采用峰平谷三段比如峰时1.2元/kWh、平时0.8元/kWh、谷时0.4元/kWh天然气价格取2.8元/m³并折算成单位热值价格。碳价基准设为50元/tCO2排放因子按0.6 kgCO2/kWh电网和0.2 kgCO2/kWh天然气来设。这些参数的选取没有放诸四海皆准的标准关键是保持内部一致性。我见过不少论文把天然气价和碳价取到某个比例之后结果里燃气设备从头到尾都不开机这明显是参数失配。5.2 三个场景怎么设置才有对比价值我会固定所有设备参数和负荷曲线只改动三个变量场景A不引入碳交易、不考虑DR即传统经济调度场景B引入碳交易机制但没有DR场景C碳交易和两类DR同时启用。这个三场景对比的好处是能拆分两个机制的贡献。B和A对比看碳交易的减排效果和成本影响C和B对比看DR在碳约束下的进一步优化能力。如果想把需求响应的效果单独拆出来还可以加一个场景D只有DR没有碳交易。四个场景一起跑能展现的结论链条就非常完整了。5.3 结果里最值得关注的三个指标跑完优化后我习惯先看三样东西而不是直接贴出一堆曲线第一是总成本和成本构成变化。碳交易引入后总成本不一定上升——如果系统通过电锅炉替代燃气锅炉获得了碳收益反而可能下降。DR引入后购电成本会降但补偿成本会增加要看净效果。第二是碳排放总量的变化。这是碳交易建模的验收指标如果考虑碳交易后排放不降反升那多半是配额系数给得太松或者排放因子设置有问题。第三是负荷曲线形态。把场景A和场景C的L_baseL_dr_price-L_dr_inc画出来对比峰谷差缩小比例。这个指标能直观反映DR有没有起作用也是报告里最有说服力的图之一。我还会把每个时段的CHP出力、储能SOC、购电量画在一张图上检查是否存在不合理的跳变——比如SOC突然从0.9掉到0.1再冲回0.9那大概率是储能约束写错了。6. 实际调试中的坑与排查思路6.1 求解器报infeasible时的排查顺序Yalmip/Gurobi返回不可行时第一反应不要去看数学公式先做三件事一是检查功率平衡约束的方向对不对。等号约束最容易因为符号习惯出错尤其是把售电、储能放电、DR削减这些带方向性的量写反一个符号模型立刻不可行。二是检查储能SOC约束的初始条件。如果SOC(1)没赋初值或者SOC(T)被限定为0.5而T时刻前后出现了奇怪的强制关系也很容易不可行。我习惯先松开末端SOC约束跑一遍确认模型是真不可解还是约束太紧。三是检查0-1变量的互斥约束。比如储能充放互斥如果两个0-1变量之和≤1约束写成了≥1系统要求每个时段必须充或放在某个零边界时段就会冲突。另一个定位技巧是逐组注释约束先把全部约束注释掉只跑目标函数然后逐步放开约束组看哪一组放开之后才变得不可行问题就出在哪里。这个方法笨但极有效。6.2 阶梯碳价0-1变量导致的非线性陷阱碳交易阶梯成本如果实现不当很容易引入sdpvar的乘积项。举个例子有人在Yalmip里直接写C_carbon 0.5 * p_c * lambda1 * z1 1.1 * p_c * (delta_E - lambda1) * z2;这里z2是binvardelta_E是sdpvar两者相乘就是一个双线性项。Gurobi不能直接处理Yalmip要么把它转为非凸二次规划要么报错。解决方案就是我2.2节说的把每个梯度区间的排放量拆成独立的非负连续变量E1、E2、E3让0-1变量只和这些切片变量的上限绑定而不是直接参与乘法。此外大M量的选择也容易出问题。M取太小会导致可行域被错误截断M取太大会引起数值振荡。我通常按该变量的物理上限×10来定E的物理上限就是最大排放量算一下单位数量级再给M不要随手写个1e6。6.3 DR参数不合理导致的伪优化结果需求响应参数如果设得太激进会出现一种看起来漂亮但不真实的调度结果系统通过DR把大量负荷挪到谷时然后谷时电价和碳排放同时很低设备全开总成本大幅下降峰谷差几乎抹平。现实中这不可能因为用户的负荷转移意愿有限。我的经验是价格型DR的弹性系数绝对值不要超过0.5激励型DR的中断容量不要超过峰值负荷的15%中断次数限制在2~3次以内。另外DR削减的负荷总量要监控削减时段负荷下降转移时段负荷上升但如果转移后的负荷超过了设备总出力上限又会出现不可行。这时候要回头检查弹性矩阵配置是否合理而不是盲目加设备容量。最后分享一个我自己的调试习惯每跑完一组场景我会把SOC曲线、DR削减量曲线和碳成本曲线同时画出来看。这三个量的行为能反映模型里90%的逻辑错误。调度优化这东西模型能跑通只是万里长征第一步结果能不能解释得通、物理上是否自洽才是最花时间的地方。希望这篇分享能帮你少走一段弯路。