简介面向综合能源系统优化调度方向研究者与技术开发者的MATLAB实战文档聚焦阶梯型碳交易机制与电制氢模块的联合建模求解。内容围绕传统机组、风电、电解槽、储氢罐与燃料电池五个单元展开讲解阶梯碳成本的分段线性化处理、氢能系统状态方程的时间耦合写法以及兼顾经济性与碳排放的目标函数构建并配有案例结果与可视化输出读者可据此理解风电大发时段优先制氢、碳排放逼近阈值时火电出力出现拐点等现象的成因。压缩包内仅含1个PDF文件共129KB篇幅精炼适合打印或按章节精读。目前已有163人学习可供具备MATLAB与YALMIP基础、希望把碳交易与氢能环节接入调度模型的读者参考其中对big-M取值与索引对齐等实现细节的讨论也有助于排查模型不收敛或储能异常等常见问题。1. 阶梯碳价与电制氢为什么分段计价会把调度问题从 LP 逼成 MILP把 MATLAB 和 YALMIP 摆到综合能源系统优化这道题前面多数人的第一反应是把火电、风电、电解槽的出力变量一列目标函数一写optimize一跑就完事。真跑起来才会发现一旦把阶梯碳交易和电制氢塞进去模型的性质就变了原本一个线性规划能秒解的调度问题会因为碳价的分段跳变和电解槽的启停逻辑变成必须靠分支定界才收敛的混合整数线性规划。更反直觉的是阶梯碳价并不只是给模型加了一个环保指标它改变的是每个时段的边际成本符号于是火电出力曲线会在阈值附近出现拐点而氢能设备在特定参数下反而变成了省钱的那一个。搞懂这条从线性到分段线性的分界线在哪比会调几个求解器参数有用得多。2. 阶梯碳交易为什么 if-else 写不出可解的分段成本2.1 阶梯碳价的数学结构本质是一个凸分段线性函数先把这个函数写清楚。设某时段总排放量为 $E$碳价分档阈值为 $T_1 T_2 \dots$对应单价 $p_1 p_2 \dots$那么碳成本是$$ C(E) p_1\min(E,T_1) p_2\min\big(\max(E-T_1,0),,T_2-T_1\big) \dots $$这是个分段线性的凸函数。凸这个字是全局的关键——它意味着边际成本随排放量单调不减多排一吨永远比少排一吨贵。很多人写代码时第一反应是照抄一个if-elseif-else函数像这样% 反例直接写成分段判断看起来直观实际上不可用 function cost carbon_cost(emission) if emission 200 cost emission * 80; elseif emission 600 cost 200*80 (emission-200)*120; else cost 200*80 400*120 (emission-600)*160; end end问题在于emission在优化模型里是决策变量的表达式各机组出力乘以排放因子累加而来不是数值。YALMIP 遇到if分支时无法判断走哪条路只能退化成非线性约束甚至直接报错说表达式不支持。这类写法在调试阶段可能还能跑通一旦把 24 个时段串起来就彻底罢工。2.2 辅助变量拆段把跳变摊平成一组线性约束正确姿势是引入一组非负辅助变量让它们分别表示落在第 k 档里的那部分排放量。总量是各分量之和成本是各分量乘以对应单价再求和阶梯的上限用不等式约束表达。% 阶梯碳交易的分段线性建模 tier_upper [200, 600, 1200]; % 前三档的累计上限吨最后一档开口 price [80, 120, 160, 220]; % 各档单价元/吨 width diff([0, tier_upper]); % 每档的可用宽度[200 400 600] nTier numel(price); delta sdpvar(nTier, 24); % 每个时段、每一档上的排放分量 cons []; for t 1:24 % 分解恒等式各档分量之和必须等于该时段真实排放量 cons [cons, sum(delta(:,t)) E_total(t)]; cons [cons, delta(:,t) 0]; % 前 nTier-1 档有宽度上限最后一档不设限吸收全部超出部分 for k 1:nTier-1 cons [cons, delta(k,t) width(k)]; end end % 碳成本各档分量 × 对应单价跨时段求和 carbon_cost sum(price * delta, all);这段代码里没有任何if全部是线性等式与不等式。E_total(t)是由机组出力拼出来的线性表达式delta是求解器要决定的变量二者通过分解恒等式绑定。参数含义上tier_upper是累计阈值而不是每档宽度所以要先做一次diff最后一档故意不设上限否则当排放超预期时模型会直接不可行这是新手最容易踩的一个坑。2.3 凸性为什么能省掉二进制变量以及什么时候省不掉上面的模型能成立靠的是一个不太显眼但非常重要的事实当价格递增时线性规划松弛会自动给出贪心解——先填满便宜的档再填贵的。因为目标是最小化成本求解器没有任何动机把排放量塞进高价档而让低价档空着。也就是说不需要 0-1 变量来标记当前落在哪一档模型自动就落对了位置。反过来如果哪天政策设计成用得越多单价越低价格递减这个性质立刻失效。LP 松弛会把全部排放量塞进最便宜的最高档得到的成本远低于真实值求出来的调度方案是假的。这时候必须上 SOS2 或者二进制变量。SOS2Special Ordered Set of type 2的写法是把分段函数表示成断点与权重的凸组合% SOS2 写法用断点 权重表示分段线性函数不依赖凸性 E_bp [0, 200, 600, 1200, 3000]; % 排放量断点 C_bp [0, 16000, 64000, 160000, 556000]; % 对应的累计碳成本 lam sdpvar(numel(E_bp), 1); % 每个断点的权重 cons [cons, sum(lam) 1, lam 0, sos2(lam)]; % 权重和为 1且最多两个相邻非零 cons [cons, E_t E_bp * lam, carbon_t C_bp * lam];SOS2 的约束语义是权重向量中至多两个分量非零且必须相邻它强制解落在某一段的端点上天然绕开凸性要求。代价是分支定界树的规模变大——CPLEX 对 SOS 约束做的是特殊分支节点数通常比纯 LP 放松多。所以取舍很清楚价格递增就用辅助变量法纯 LP 松弛最快价格非递增或者要做任意形状的分段曲线才动 SOS2。我一般会先按辅助变量写跑出来发现成本偏低得离谱再回头检查是不是价格序列本身不满足递增。3. 电制氢链路乘效率还是除效率决定了储氢量会不会变成负数3.1 能量流的方向性电到氢是乘氢到电是除电制氢这条链路有三个设备、两个方向。电解槽吃电产氢燃料电池吃氢发电储氢罐夹在中间做时间上的搬运。效率的定义都是有效输出 / 输入所以两个方向的写法天然不对称电解槽输入电功率 $P_{ele}$输出氢能量 $ \eta_{ele} \cdot P_{ele}$是乘法燃料电池输出电功率 $P_{fc}$消耗氢能量 $ P_{fc} / \eta_{fc}$是除法。写成代码就是下面这样注意符号和时间戳的对齐% 储氢罐状态递推注意索引从 t2 开始功率从 t-1 取值 for t 2:24 cons [cons, H2(t) H2(t-1) ... eta_ele * P_ele(t-1) * dt ... % 电解槽功率 × 效率产氢 - P_fc(t-1) / eta_fc * dt]; % 燃料电池功率 ÷ 效率耗氢 end % 日周期闭合首末状态相等否则模型会免费消耗或囤积氢气 cons [cons, H2(1) H2_init, H2(24) H2_init]; cons [cons, H2_min H2 H2_max];dt是时段长度小时P_ele单位是 kW乘完得到 kWh如果储氢量以 kg 计量还要再除以氢的低位热值常见取值约 33.3 kWh/kg做单位统一。这一步在论文里经常被一笔带过实际写代码时漏掉就会导致量纲对不上——储氢罐容量设成 500结果算出来是几万模型看着可行调度曲线却完全没有物理意义。3.2 索引偏移那个让储氢量在最后一小时崩掉的细节H2(t)依赖H2(t-1)和t-1时刻的功率这个错位是刻意的。第 t 小时结束时的储氢量等于第 t 小时开始时的储量加上第 t 小时内电解槽产出的氢减去第 t 小时内燃料电池消耗的氢。如果写成P_ele(t)和H2(t-1)配对等价于假设这一小时的产出在小时开始前就已经到账时间上穿越了。后果不是报错而是末端状态对不上由于循环从t2走到t24实际只约束了 23 个时段的状态转移第 24 小时的功率完全没有进入任何状态方程模型可以放心地在最后一小时满功率制氢或者满功率发电而储量不变。表现出来就是末时段储氢曲线突然跳变、或者日周期闭合约束莫名其妙不可行。稳妥做法是在写完后加一条自检把所有时段的H2解出来手动按物理公式递推一遍看能不能对上对不上就一定是索引或者符号写错了。3.3 启停次数、最小负载率与 0.3 这个系数电解槽不是想开就开、想调多低就多低。两个工程约束必须进模型最小负载率和启停次数上限。u_on binvar(24, 1); % 电解槽运行状态 start binvar(24, 1); % 启动动作标记 for t 1:24 cons [cons, P_ele(t) 0.3 * P_ele_max * u_on(t)]; % 最小负载率 30% cons [cons, P_ele(t) 1.0 * P_ele_max * u_on(t)]; % 最大功率 end % 启动动作只在状态由 0 变 1 的时刻为 1其余时刻为 0 cons [cons, start(2:end) u_on(2:end) - u_on(1:end-1)]; cons [cons, start 0, sum(start) 3]; % 24 小时内最多启停 3 次0.3 这个下限系数的来源是电解槽的极化曲线负载率过低时单位产氢能耗急剧上升效率会掉到一个不划算的区间同时低负载运行对电极寿命也不友好。把它写成P_ele 0.3 * P_max * u_on而不是P_ele 0.3 * P_max前者允许停机后者强制设备 24 小时都开着是完全不同的两个模型。启停次数的约束有个隐蔽写法问题。很多资料里直接写sum(abs(diff(u_on))) 3abs和diff作用在决策变量上会让 YALMIP 生成非线性表达式MILP 结构当场被破坏求解器要么拒绝要么慢到不可接受。正确做法就是上面那个start辅助变量约束它不小于状态变化量再由目标函数或者上界把它压住。这条属于典型的看起来能写、跑起来翻车。4. 目标函数与求解器MILP 结构怎么决定收敛速度4.1 四项成本的量纲统一与权重设计目标函数长这样% 四项成本燃料、运维、弃风惩罚、阶梯碳成本 cost_fuel sum(C_fuel .* P_gen, all); % 元机组燃料成本 cost_om sum(C_om .* P_ele C_om_fc .* P_fc); % 元氢能设备运维 cost_curt sum(C_curt .* P_curtail); % 元弃风惩罚 cost_carbon sum(price * delta, all); % 元前面算好的阶梯碳成本 total_cost cost_fuel cost_om cost_curt cost_carbon;四个量都是元这是能直接相加的前提。弃风惩罚系数C_curt的取值其实是个隐含的政策变量设得太小模型宁愿弃风也不肯多开电解槽设得够大制氢消纳弃风才有经济性。我一般会把它设成接近上网电价的量级否则结果里凌晨风电高峰优先制氢这个现象根本出不来。另外要注意碳成本是整个调度周期的全局量必须放在时段循环外面累加写进循环里会导致同一笔排放被重复计费。4.2 CPLEX 参数mipgap 到底在放松什么ops sdpsettings(solver,cplex,verbose,2); ops.cplex.timelimit 3600; % 墙钟时间上限秒 ops.cplex.mip.tolerances.mipgap 0.01; % 相对 gap 1% ops.cplex.mip.tolerances.absmipgap 1e2; % 绝对 gap防止相对 gap 失效 ops.cplex.threads 1; % 单线程保证结果可复现mipgap的语义是(UB - LB) / |UB|即当前最好整数解与线性松弛下界之间的相对差距。设成 0.01 意味着允许 1% 的最优性损失——如果总成本是 10 万元那就是 1000 元的差距工程调度完全能接受。但要留个心眼当目标函数值本身接近 0 时比如做了归一化处理相对 gap 会变得极其苛刻求解器可能跑到时间上限都停不下来这时候必须补一个absmipgap。threads 1这条不是性能优化是可复现性保障。CPLEX 的多线程分支定界在节点选择上有非确定性同样的模型跑两遍可能给出不同的次优解论文里画出来的曲线每次都不一样会被审稿人问死。牺牲一点速度换确定性在做对比实验时是值得的。4.3 big-M 的数值敏感性为什么 1e6 是个坏习惯线性化要么为零、要么不低于某个值这类逻辑标准手法是 big-M。问题出在 M 的取值上。教科书写法喜欢用一个大数比如 1e6% 反例M 远大于物理上界 cons [cons, P_ele(t) 1e6 * u_on(t)]; % 正解M 取该设备的容量上限这是最紧的合法上界 M_ele P_ele_max; cons [cons, P_ele(t) M_ele * u_on(t)];M 过大有三个后果线性松弛几乎不起约束作用下界极松分支定界树规模指数级膨胀系数矩阵条件数变差CPLEX 报数值警告甚至给出错误结论变量量级跨度太大会让单纯形法的容差判断失准。而 M 取物理容量上限是最紧的选择——它恰好覆盖所有可行解又不会多放一丝冗余空间。原文提到用动态计算 M 值才稳定收敛本质就是同一个道理每个约束的 M 应该由它涉及的那几个变量自身的物理上界推出来而不是一个全局常数。比如功率平衡约束里的 M 用总装机储氢约束里的 M 用罐容各算各的。5. 排查这几类翻车几乎每次建模都会撞上5.1 储氢量出现负值或者日末状态回不到初值现象求解成功但value(H2)里有负数或者H2(24)和H2(1)差了一大截。原因通常是两个一是状态递推里的时间索引没有和功率对齐导致最后一小时的能量流没被记账二是只写了H2 0却忘了首末状态相等模型就有了跨日免费借氢的口子。解决办法是把H2(1) H2(24)显式写进约束并在求解后手动按物理公式向前递推一遍做交叉验证。5.2 模型报不可行但明明存在显然的可行解现象CPLEX 返回 infeasible把碳价档位调宽一点又能跑了。原因往往是分段上限写死导致的最后一档也设了宽度上限一旦排放量超出所有档的总宽度分解恒等式无解。解决是把最后一档的宽度约束去掉让它成为开口的溢出档同时给个足够大的物理上界比如所有机组满发时的最大排放量。5.3 电解槽在 0 和满功率之间反复跳变现象调度曲线上电解槽功率呈现方波一小时开一小时停。原因是模型里只有最小负载率、没有启停成本或启停次数约束而电价或碳价在相邻时段有小幅波动求解器就利用跳变去贴每一个局部最优点。解决是加start变量并把启停次数限制在合理范围或者在目标函数里加一项启动惩罚让频繁切换付出代价。5.4 碳成本算出来比手工核算低一截现象模型收敛目标函数里的碳成本看着挺合理但拿解出来的排放量代回分段函数一算对不上。原因多半是价格序列本身不满足递增LP 松弛钻了空子把排放塞进了低价档。解决是先确认diff(price) 0不满足就换 SOS2 或者引入 0-1 变量标记所在档位。5.5 同样的代码跑两遍结果不一样现象参数没动两次运行的调度曲线和总成本都有差异。原因是求解器多线程导致的非确定性或者 YALMIP 内部变量命名顺序变化影响了分支顺序。解决是把cplex.threads设为 1并在每次运行前yalmip(clear)保证模型构建的起点干净。做论文对比实验时这一条几乎是硬性要求。6. 设计取舍氢能省下的那笔钱到底是谁付的跑出一版结果看到引入电制氢后总成本下降 8%、最高碳价档从未被触发很容易得出一个乐观结论氢能投资能靠碳交易机制赚回来。这个结论对参数极度敏感换一组数就可能完全反转。真正决定结论的是下面几个量参数基准取值变化后结论走向最高档碳价220 元/吨降到 150氢能替代收益基本消失火电照常出力电解槽投资成本未计入日运行成本按容量摊销计入8% 的节省直接被摊薄吃光弃风惩罚系数接近上网电价降为 0制氢失去消纳弃风的附加价值电解效率模型恒定 0.72分段0.6/0.72/0.68拐点位置移动满功率时段缩短表格里第一行和第二行是最容易被忽略的。很多模型只算运行成本把电解槽、储氢罐的资本开支当成沉没成本排除在外这样当然能算出氢能省钱——省的是碳成本付钱的是没被计进来的设备折旧。公平的比较应该把日运行成本里加一项容量摊销再去看总账。第三行则揭示了另一个机制制氢的经济性有相当一部分来自它替系统消化了弃风如果弃风惩罚系数设得很低模型根本不在乎弃多少风制氢的这块隐性收益就归零了。还有一个更微妙的取舍在时间尺度上。储氢罐的容量约束天然是跨时段的它把 24 小时耦合成一个整体问题求解规模比纯发电调度大一圈而如果为了求解速度做时段解耦比如滚动优化、每小时独立求解储氢的状态就断了氢能的核心价值——把凌晨的风电搬到晚高峰用——也就没了。我一般会在做敏感性分析时先跑一遍全日耦合的版本作为基准再用滚动窗口跑一遍对比两者的成本差这个差值本身就是时间耦合值多少钱的直接度量。从那以后我每次搭这类含阶梯碳价和储氢的模型都会强制走一遍固定流程先检查价格序列是否单调递增再用物理公式手工递推一遍储能状态最后把 big-M 全部换成物理容量上限。这三步看着笨但能挡掉绝大部分求解器给出的假收敛。希望帮到你。本文还有配套的精品资源点击获取