简介面向能源系统优化研究者及具备MATLAB基础的技术人员这份PDF聚焦阶梯型碳交易与电制氢综合能源系统的优化调度帮助读者理解传统机组、风电、电解槽、储氢罐和燃料电池五单元协同下的经济性与减排平衡。内容以YALMIP调用Cplex求解为主线覆盖阶梯碳成本分段线性化与big-M处理、储氢罐动态状态方程、电-氢-电双向效率折算、目标函数设计以及结果可视化并结合案例分析风电大发时段优先制氢、碳排放接近阈值时火电出力调整等现象。整份文档围绕模型搭建、参数初始化与结果解读组织关键约束和易错细节有较充分说明便于对照复现。压缩包共1个PDF文件大小129KB已有162人学习适合参考约束写法、求解设置与数值稳定性处理细节也可用于能源优化课题复现、课程设计及论文实验参考。1. 一套会算碳账的调度器MATLAB 里跑通阶梯碳交易与电制氢如果只按固定碳价算成本综合能源系统的调度结果往往很温和火电多出一点力交点碳费就过去了。可一旦把碳价改成阶梯制再塞进一整套电制氢链路模型的行为就会突然变得很不线性某个时段的碳排放只要跨过阈值后面的边际成本立刻跳档而电解槽、储氢罐、燃料电池又要求你在时间轴上前后对齐。这份 MATLAB YALMIP 的代码包干的就是这件事传统机组、风电、电解槽、储氢罐、燃料电池五个单元24 点逐时调度目标是最小化燃料成本、弃风惩罚、运维与阶梯碳成本之和约束交给 CPLEX 求解。它适合已经能看懂优化模型、想把碳机制和氢储能真正落进代码里的开发者也适合要改造现有调度程序、把碳成本从常数项变成分段函数的人。2. 系统架构拆解五单元拓扑与三个 .m 文件的分工2.1 能量流拓扑母线、弃风变量与氢的往返路径整个系统挂在一条电母线上供电侧有三台可调火电、一片风电场和一个燃料电池用电侧是固定负荷加一台电解槽。风电按预测曲线出但不一定全被消纳所以模型里专门留了一个弃风变量电解槽是可控负荷耗电产氢氢进储氢罐储氢罐再给燃料电池供氢燃料电池又回到母线上发电。这条电极到氢、氢再到电的往返路径就是整个模型区别于普通机组组合问题的地方。功率平衡约束写出来是% 母线功率平衡发电侧 用电侧每个时段一条约束 for t 1:T cons [cons, ... sum(P_gen(t,:)) (P_wind(t) - P_curt(t)) P_fc(t) ... P_load(t) P_ele(t)]; end这里的逻辑是「火电总和 实际消纳风电 燃料电池出力 负荷 电解槽耗电」。参数说明P_wind(t)是 t 时段风电可用出力预测值属于常量P_curt(t)是弃风功率是决策变量上行被P_wind(t)卡住不能凭空多弃P_fc(t)是燃料电池发电功率。把弃风写成变量而不是简单丢弃是为了在目标函数里给它一个惩罚价模型才知道该尽量消纳。设备参数统一放在初始化结构体里参数字段含义参考取值单位gen_cost三台火电线性成本系数[28 31 25]元/MWhgen_pmax/gen_pmin单机出力上下限[80 60 50]/[10 8 5]MWco2_int三台机组碳排放强度[0.85 0.92 0.78]t/MWhwind_cap风电装机容量150MWele_eff电解槽额定效率0.72—fuel_cell_eff燃料电池发电效率0.58—ele_pmax/ele_pmin电解槽功率上下限40/12MWh2_tank_max储氢容量折算为等效电量300MWhcurtail_penalty弃风惩罚单价200元/MWh需要强调单位储氢罐里存的不是 kg而是「等效电量 MWh」也就是罐里那部分氢如果拿去发电最多能出多少电。这样一来电解槽乘效率、燃料电池除效率就能直接加减不用来回换算热值。建模时所有量统一到 MW 和 MWh 上是后面不出现诡异负储氢的前提。2.2 三文件分工sys_init、build_model、plot_results 各管一段代码拆成三个文件边界划得很清楚。sys_init.m只做一件事把上面那张表拼成一个units结构体再生成 24 点负荷曲线和风电预测曲线顺手算出几个派生量比如总装机、M 值的上界、储氢容量对应的最大连续制氢时长。build_model.m负责全部建模和求解输入units输出解结构体sol和一个诊断结构体diag求解状态、耗时、各分项成本。plot_results.m只读不写拿sol画三张图出力堆叠图、储氢量与碳排分档曲线、成本构成饼图。调用骨架就这么几行clear; clc; yalmip(clear); % 清空 YALMIP 内部缓存避免上次的模型残留 units sys_init(); % 生成参数结构体 24 点预测数据 [sol, diag] build_model(units); % 建模、求解、成本拆分 plot_results(sol, units); % 出力曲线 / 氢罐曲线 / 碳档位 if diag.exitflag ~ 1 warning(求解未正常收敛code %d先查第 5 章的排查清单, diag.exitflag); end这样拆的好处是二次开发时改动面可控换数据只动sys_init.m加约束只动build_model.m出图样式完全不碰。diag里我习惯把碳成本、燃料成本、弃风成本、启停成本四项单独存一份因为调试阶段最常见的错误就是总成本算错了却不知道错在哪一项有了分解就能一眼看出是不是碳成本被放大了 1000 倍这种低级事故。2.3 决策变量清单维度错一次后面全乱变量定义写在build_model.m开头维度必须一次定对。火电用T×N_gen的矩阵其余时序量用T×1列向量启停状态用binvar。YALMIP 里最容易翻车的地方是sum的默认行为对矩阵做一次sum(P_gen)得到的是1×3行向量而不是标量直接扔进目标函数会报 nonscalar 错误。T units.T; % 24 N_gen numel(units.gen_cost); % 3 台火电 P_gen sdpvar(T, N_gen, full); % 火电出力 MW列对应机组 P_ele sdpvar(T, 1); % 电解槽耗电功率 MW P_fc sdpvar(T, 1); % 燃料电池发电功率 MW H2 sdpvar(T, 1); % 储氢量等效电量MWh P_curt sdpvar(T, 1); % 弃风功率 MW on_ele binvar(T, 1); % 电解槽开停机状态 0/1 u_ele binvar(T, 1); % 启动动作 0/1 E_seg sdpvar(1, 3); % 阶梯碳排放的三段量见第 3 章配套的基础约束也一并给出机组上下限、弃风上限、氢罐容量这三组是每个时段都要满足的cons []; cons [cons, repmat(units.gen_pmin, T, 1) P_gen repmat(units.gen_pmax, T, 1)]; cons [cons, 0 P_curt P_wind]; % 弃风不能超过可用风资源 cons [cons, units.h2_min H2 units.h2_tank_max]; % 留 20% 安全余量repmat那一行是把1×3的上下限扩成T×3YALMIP 支持这种逐元素双边不等式写起来比循环干净。氢罐下限设成容量的 20% 而不是 0是为了避免模型为了省钱把罐子抽干、第二天开局无氢可用这个细节在第 4 章讲终端约束时还会再提一次。变量都定成连续量还是整数直接决定求解时间只有启停和分段选择用二进制功率和储量一律连续24 点的规模下 CPLEX 通常几秒内出解。3. 阶梯碳交易的建模把 if-else 变成 MILP 的分段线性化3.1 为什么不能直接写 if-else 和 max很多人第一版会写成if emission 200, cost 80*emission; ...这种形式跑起来发现要么报错要么结果完全不对。原因是emission是sdpvar它的值在建模阶段根本不存在MATLAB 的if拿到一个未知对象只能按对象真假走固定分支等于把你的分段逻辑整个跳过了。另一种写法是用min、max拼成本函数比如80*min(E,200) 120*max(0,E-200)读起来漂亮但max是分段线性的非光滑算子YALMIP 会把它解释成非线性算子CPLEX 的 MILP 分支根本不接受求解器直接告诉你模型类型不匹配。正确做法是在建模阶段就把分段结构展开成「若干非负变量 线性等式 上下界」让整个模型保持纯 MILP。这也是为什么这份代码里碳排放量不是一个标量变量而是一个长度为 3 的分段向量。碳排总量本身是线性的由各机组出力和排放强度相乘求和得到加一个固定排放基数比如供热分摊部分便于对齐实际园区的口径% 碳排总量火电出力 × 排放强度逐机组逐时段求和 total_emission sum(P_gen * units.co2_int, all) units.base_emission;注意all这个参数不加的话sum会对矩阵按列求和返回行向量total_emission就不是标量后面所有跟它相关的等式约束都会维度不匹配。这是我在别人代码里见过最频繁的低级错误之一报错信息还特别不友好只说约束维度不一致不说到底哪一行。养成习惯凡是求和得到标量的地方都显式写all或者套两层sum(sum(...))。3.2 三段填充法递增碳价下的最小写法代码包默认用一套三段递增碳价前 200 吨每吨 80 元200 到 500 吨每吨 120 元超过 500 吨每吨 160 元。实现上引入三个非负变量代表落在各档里的排放量加起来等于总排放各档再卡上自己的容量上限成本就是三个量分别乘单价再求和% 阶梯碳交易三档递增价格的分段线性化 co2_seg [200, 500]; % 档位分界点单位 tCO2 co2_price [80, 120, 160]; % 各档单价元/tCO2 E_seg sdpvar(1, 3); % E_seg(1) 落第一档以此类推 cons [cons, sum(E_seg) total_emission]; % 三段之和 总排放 cons [cons, 0 E_seg(1) co2_seg(1)]; % 第一档容量 200 t cons [cons, 0 E_seg(2) co2_seg(2) - co2_seg(1)]; % 第二档容量 300 t cons [cons, E_seg(3) 0]; % 第三档无上限 carbon_cost co2_price(1)*E_seg(1) co2_price(2)*E_seg(2) co2_price(3)*E_seg(3);这段代码里没有二进制变量靠的是两个性质共同成立单价严格递增且目标函数是最小化。求解器为了让成本最小会本能地先把便宜的第一档填满再填第二档绝不会出现「第一档没满却先用第三档」的浪费解所以线性松弛本身就等于整数最优不需要额外加 0-1 变量。一旦政策设计成第二档比第一档便宜某些减免区间或者第三档有封顶返还这个性质立刻失效模型会钻空子把排放全塞进便宜档这时候就必须上二进制变量强制顺序填充。判断标准很朴素把各档单价排成一列看是否单调递增是就用填充法不是就加整数变量。3.3 SOS2、二进制与动态 big-M三种写法的取舍除了填充法工程上还有两种常见写法。一种是特殊有序集 SOS2用一组权重 λ 表示当前解落在哪两个相邻断点之间CPLEX 原生支持这种结构分支效率比二进制高缺点是 YALMIP 里要手工把 λ 和断点绑定写错索引就退化成普通连续变量约束形同虚设。另一种是最直白的二进制选择给每一档配一个 0-1 变量表示「是否进入该档」再用 big-M 把前后档位串起来% 二进制顺序填充保证低档填满才允许用高档 z binvar(3, 1); cons [cons, sum(z) 1]; cons [cons, E_seg(1) co2_seg(1) * z(2)]; % 用第二档则第一档必满 cons [cons, E_seg(2) (co2_seg(2)-co2_seg(1)) * z(3)]; % 用第三档则第二档必满 cons [cons, E_seg(1) co2_seg(1)]; cons [cons, E_seg(2) (co2_seg(2)-co2_seg(1))];这里z(2)1表示已经进入第二档于是强制第一档的下界等于它的容量。真正的坑在 big-M 的取值如果模型里还有别的逻辑用 M 表示「足够大的数」很多人顺手写1e6CPLEX 的对偶边界立刻变得极端病态分支定界在节点上反复震荡跑十几分钟给不出可行解。我的做法是让 M 由上界推出来而不是拍脑袋% 动态 M一天里火电可能产生的最大碳排作为任何分段逻辑的紧上界 M_big T * sum(units.gen_pmax .* units.co2_int) units.base_emission;这个值通常只有几千吨量级跟真实解的差距在两三倍以内数值条件好得多。三种写法的对比如下选哪个主要看碳价是否单调、模型规模有多大写法需要整数变量求解速度适用场景三段填充否最快单价严格递增绝大多数政策场景SOS2 分段求解器内部特殊结构较快断点多于 3 段、需要精确拟合曲线二进制 big-M是每档一个较慢单价非单调、存在减免或返还区间4. 电制氢链路电解槽效率、储氢罐状态方程与启停约束4.1 电解槽的功率-产氢关系二次拟合好看进 MILP 就翻车电解槽的实际效率不是常数低负载时效率明显下滑中高负载区间反而最稳。文献里常见的做法是拟合一条二次曲线比如产氢量关于功率写成0.45*P - 2e-4*P^2这种形式。这条曲线本身没错问题在于它带着平方项一旦写进目标函数或约束整个模型就从 MILP 变成 MIQCPCPLEX 要动用二阶锥或者空间分支定界24 点三机组的小模型还能撑住扩到一周或者加进更多设备求解时间就不是翻倍而是翻十倍。代码包的处理方式是离散化取负载率的三个断点 0、0.5、1.0对应效率 0.60、0.72、0.68用一组权重 λ 和一组段选择变量 z 把「当前工作在哪个区间」表达清楚。这种写法本质上是把曲线掰成两段折线误差控制在 2 个百分点以内但模型保持纯线性% 电解槽效率分段线性化断点 (0, 0.5, 1.0)对应效率 (0.60, 0.72, 0.68) bp [0, 0.5, 1.0]; eta [0.60, 0.72, 0.68]; lam sdpvar(T, 3); % 每个时段一组断点权重 z binvar(T, 2); % 两个线性段的选择变量 for t 1:T cons [cons, sum(lam(t,:)) 1, sum(z(t,:)) 1]; cons [cons, lam(t,1) z(t,1)]; cons [cons, lam(t,2) z(t,1) z(t,2)]; % 相邻断点才允许非零 cons [cons, lam(t,3) z(t,2)]; cons [cons, P_ele(t) units.ele_pmax * (bp * lam(t,:))]; cons [cons, H2_prod(t) units.ele_pmax * ((bp .* eta) * lam(t,:))]; end逻辑说明lam表示当前工作点落在两个相邻断点之间的权重组合z负责把非相邻的权重掐成 0这就保证了「同一时刻只在一段折线上」。P_ele(t)由断点值加权得到H2_prod(t)用断点上的产氢量加权两者用同一组 λ所以功率和产氢严格一一对应。参数上ele_pmax是 40 MW意味着满负载时段产氢约 28.8 MWh 等效电量中途效率峰值在 0.72模型会自然倾向于让电解槽跑在半载以上这也解释了为什么结果图里电解槽很少出现在 20% 负载以下的区间。这套写法比二次曲线的精度差一点但换来了可预测的求解时间工业调度里我觉得这笔交易划算。4.2 储氢罐状态方程索引偏移是负储氢的元凶储氢罐是整条链路里唯一带时间耦合的环节也是 bug 高发区。它的物理含义很直白这一时刻罐里的氢 上一时刻的余量 本时段电解槽制氢 − 本时段燃料电池耗氢。用等效电量表达之后耗氢那项要除以发电效率因为燃料电池发 1 MWh 电需要消耗 1/0.58 ≈ 1.72 MWh 的氢。这个「乘还是除」的方向如果搞反模型会认为燃料电池越发电罐里氢越多跑出来的结果看着收敛实际物理上完全不成立。% 储氢罐动态初值 制氢乘效率 - 耗氢除效率 cons [cons, H2(1) units.h2_init ... units.ele_eff * P_ele(1) - P_fc(1) / units.fuel_cell_eff]; for t 2:T cons [cons, H2(t) H2(t-1) ... units.ele_eff * P_ele(t) - P_fc(t) / units.fuel_cell_eff]; end % 日终闭合末时段储量不低于日初防止吃空库存 cons [cons, H2(T) units.h2_init];这里我把每个时段的功率都按「本时段」写而不是原文里那种「t 时刻储量由 t-1 时刻功率决定」的滞后写法。两种离散方式都能用但滞后写法有个致命副作用最后一个时段的功率不会影响任何储量模型会把最便宜的时段全堆在第 24 点图上看起来末时段制氢功率爆表、罐子却没变化。统一成本时段的写法再补一条日终闭合约束问题就消失了。h2_init一般取容量的 30% 到 40%再配合第 2 章那条 20% 下限罐子全天都在 20% 到 100% 之间游走不会出现负储量。4.3 启停次数与最小出力0.3 这个系数从哪来电解槽是贵设备不能让它每小时在 0 和满功率之间来回横跳。模型里用两个二进制变量表达这件事on_ele表示当前时段是否开机u_ele表示当前时段是否发生了「从停到开」的动作。最小出力约束靠on_ele撑起来启动次数约束靠u_ele累加% 电解槽开停逻辑与最小出力 cons [cons, P_ele units.ele_pmax * on_ele]; % 停机时功率必须为 0 cons [cons, P_ele units.ele_pmin * on_ele]; % 开机则不低于 12 MW cons [cons, u_ele(2:T) on_ele(2:T) - on_ele(1:T-1)]; % 状态跃迁即一次启动 cons [cons, u_ele(1) on_ele(1) - units.on_init]; % 首时段也要算启动 cons [cons, sum(u_ele) units.max_start]; % 全天不超过 3 次ele_pmin取 40 MW 的 0.3 倍也就是 12 MW这个系数不是拍脑袋来的低于 30% 负载时电解效率会明显掉到 0.5 以下单位氢的耗电量飙升与其这样不如让设备停着。max_start默认 3是设备厂家给的日均启停次数建议值。u_ele(1)那行很容易漏漏了之后模型会把「凌晨先开一次机」当成免费的启动成本少算一次。启动成本在目标函数里按次计价通常设成 200 到 500 元一次具体看设备容量数值不需要很准它的作用是把频繁启停的收益压下去。5. 求解配置与踩坑排查从 MIP gap 到负储氢的五条血泪记录5.1 目标函数组装与 CPLEX 选项分项成本必须能拆开看目标函数一共四项燃料成本、运维成本、弃风惩罚、阶梯碳成本再加上电解槽启动成本这一小块。写的时候逐项算完再相加不要图省事写成一个巨大的表达式否则调参时根本看不出是哪一项在主导结果。% 分项成本每一项单独成变量方便事后拆解 fuel_cost sum(sum(P_gen .* repmat(units.gen_cost, T, 1))); om_cost units.om_coef * sum(P_gen, all) units.h2_om * sum(P_ele); curt_cost units.curtail_penalty * sum(P_curt); start_cost units.start_cost * sum(u_ele); total_cost fuel_cost om_cost curt_cost start_cost carbon_cost; % 求解器选项MILP 三条常用设置 ops sdpsettings(solver, cplex, verbose, 1); ops.cplex.timelimit 600; ops.cplex.mip.tolerances.mipgap 0.01; ops.cplex.mip.tolerances.absmipgap 50; sol optimize(cons, total_cost, ops); diag.exitflag sol.problem; diag.fuel value(fuel_cost); diag.carbon value(carbon_cost);参数说明mipgap 0.01表示允许 1% 的相对最优间隙调度问题里这个精度完全够用实际总成本几百万元的时候误差几千元低于风电预测误差带来的不确定性absmipgap 50是绝对间隙作用是当目标函数量级很小的时候不让求解器为了 1% 死磕timelimit必须设没有时间上限的 MILP 在极端数据下可能跑一整夜。verbose 1会打印分支节点数、当前最好解和下界调试的时候看这两个值的差距缩得快不快比盯着进度条有用得多。5.2 五条踩坑记录现象、原因、解决现象一结果里储氢量出现负值。原因几乎总是效率折算方向写反或者状态方程的时序错位导致末时段没闭合。解决方式是先检查H2(t)是不是被h2_min这类下界约束真正包住了再核对制氢乘效率、耗氢除效率的方向如果模型仍然给出负值把H2的下界临时收紧到 0 跑一遍对比两次结果就能定位是约束漏写还是求解精度问题。现象二模型显示 infeasible但怎么看参数都合理。最常见的原因是储氢容量和制氢功率不匹配40 MW 电解槽满发一小时产 28.8 MWh 等效氢如果罐子只给了 30 MWh 容量早上四个小时连发就把罐子撑爆模型只能报不可行。解决办法是把h2_tank_max和ele_pmax联动设计经验上罐容至少能装下电解槽满发 6 到 8 小时的产量或者干脆给弃风留出更大空间。现象三求解时间从几秒变成十几分钟。十有八九是某处 big-M 取了1e6这类夸张值或者碳排总量用了一个很宽的人工上界。把 M 换成由设备容量推出来的紧上界之后同一个模型经常在一分钟内收敛。这条血泪经验我踩过不止一次后来养成习惯模型里每出现一个 M都要在注释里写清它是怎么算出来的。现象四碳成本高得离谱比燃料成本还大一个数量级。单位混用。有的文献碳价按元/kg 写有的按元/t 写两者差 1000 倍还有些参数表里排放强度给的是 kgCO2/MWh 而代码按 t 处理。排查方法很简单把value(total_emission)打出来看量级是不是在几百到几千吨如果只有零点几说明单位错了一千倍。现象五电解槽全天频繁启停启动次数约束好像没生效。检查u_ele的定义很多人只写了u_ele(t) on_ele(t) - on_ele(t-1)却忘了首时段也忘了把u_ele放进目标函数。如果启动成本系数是 0模型会认为启停免费约束虽然限制了次数上限但 3 次之内随便跳图上就显得很乱。补上首时段约束、把启动成本调到几百元每次曲线立刻变干净。6. 二次开发从单日 24 点扩到多日滚动我一般这么改单日模型最大的局限是储氢罐必须在当天闭合跨天的能量搬移被硬生生切断。要扩成一周或者一个月的滚动调度改动集中在外层把build_model.m包进一个按天推进的循环每天调用一次把前一天的末态作为当天的初值传进去同时把日终约束从「不低于日初」改成「不低于某个保底值」只在最后一天强制闭合。% 多日滚动末态承接只有最后一天强制闭合 h2_state units.h2_init; for d 1:N_day units.h2_init h2_state; units.is_last (d N_day); % 最后一天才要求罐子回到初始水平 [sol, diag] build_model(units); h2_state value(sol.H2(end)); rec.cost(d) diag.total; end改成滚动之后验证方法也要跟着变。我一般做三件事一是拿固定碳价版本跑同一组数据做基准看阶梯碳价下的总成本是不是落在 8% 到 15% 的下降区间超出这个范围多半是参数写错二是检查能量守恒残差把每个时段的功率平衡左右两边相减取绝对值最大值应该小于 1e-6否则就是有约束没写全三是把碳档位曲线画出来看火电出力在跨档时刻有没有出现预期的拐点没有拐点说明碳价传导到调度决策的链路断了。再往深一层可以加碳捕集设备把排放压回第一档或者把风电预测换成多场景鲁棒版本。鲁棒版不用重写模型只要把P_wind从常量改成场景集合再把功率平衡约束对每个场景各写一份用binvar选最坏场景即可代价是二进制变量翻倍求解时间也会跟着上去。改造氢能部分时我每次都会先把储氢状态方程单独拿出来用固定功率输入跑一遍开环仿真确认罐子曲线单调合理了再放进优化模型不然一旦出现负储量你根本分不清是物理模型错了还是求解器给了个奇怪的解。这套排查顺序希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/28 14:50:25 · 阅读 1.2k MATLAB+YALMIP实现阶梯碳交易与电制氢调度
简介面向能源系统优化研究者及具备MATLAB基础的技术人员这份PDF聚焦阶梯型碳交易与电制氢综合能源系统的优化调度帮助读者理解传统机组、风电、电解槽、储氢罐和燃料电池五单元协同下的经济性与减排平衡。内容以YALMIP调用Cplex求解为主线覆盖阶梯碳成本分段线性化与big-M处理、储氢罐动态状态方程、电-氢-电双向效率折算、目标函数设计以及结果可视化并结合案例分析风电大发时段优先制氢、碳排放接近阈值时火电出力调整等现象。整份文档围绕模型搭建、参数初始化与结果解读组织关键约束和易错细节有较充分说明便于对照复现。压缩包共1个PDF文件大小129KB已有162人学习适合参考约束写法、求解设置与数值稳定性处理细节也可用于能源优化课题复现、课程设计及论文实验参考。1. 一套会算碳账的调度器MATLAB 里跑通阶梯碳交易与电制氢如果只按固定碳价算成本综合能源系统的调度结果往往很温和火电多出一点力交点碳费就过去了。可一旦把碳价改成阶梯制再塞进一整套电制氢链路模型的行为就会突然变得很不线性某个时段的碳排放只要跨过阈值后面的边际成本立刻跳档而电解槽、储氢罐、燃料电池又要求你在时间轴上前后对齐。这份 MATLAB YALMIP 的代码包干的就是这件事传统机组、风电、电解槽、储氢罐、燃料电池五个单元24 点逐时调度目标是最小化燃料成本、弃风惩罚、运维与阶梯碳成本之和约束交给 CPLEX 求解。它适合已经能看懂优化模型、想把碳机制和氢储能真正落进代码里的开发者也适合要改造现有调度程序、把碳成本从常数项变成分段函数的人。2. 系统架构拆解五单元拓扑与三个 .m 文件的分工2.1 能量流拓扑母线、弃风变量与氢的往返路径整个系统挂在一条电母线上供电侧有三台可调火电、一片风电场和一个燃料电池用电侧是固定负荷加一台电解槽。风电按预测曲线出但不一定全被消纳所以模型里专门留了一个弃风变量电解槽是可控负荷耗电产氢氢进储氢罐储氢罐再给燃料电池供氢燃料电池又回到母线上发电。这条电极到氢、氢再到电的往返路径就是整个模型区别于普通机组组合问题的地方。功率平衡约束写出来是% 母线功率平衡发电侧 用电侧每个时段一条约束 for t 1:T cons [cons, ... sum(P_gen(t,:)) (P_wind(t) - P_curt(t)) P_fc(t) ... P_load(t) P_ele(t)]; end这里的逻辑是「火电总和 实际消纳风电 燃料电池出力 负荷 电解槽耗电」。参数说明P_wind(t)是 t 时段风电可用出力预测值属于常量P_curt(t)是弃风功率是决策变量上行被P_wind(t)卡住不能凭空多弃P_fc(t)是燃料电池发电功率。把弃风写成变量而不是简单丢弃是为了在目标函数里给它一个惩罚价模型才知道该尽量消纳。设备参数统一放在初始化结构体里参数字段含义参考取值单位gen_cost三台火电线性成本系数[28 31 25]元/MWhgen_pmax/gen_pmin单机出力上下限[80 60 50]/[10 8 5]MWco2_int三台机组碳排放强度[0.85 0.92 0.78]t/MWhwind_cap风电装机容量150MWele_eff电解槽额定效率0.72—fuel_cell_eff燃料电池发电效率0.58—ele_pmax/ele_pmin电解槽功率上下限40/12MWh2_tank_max储氢容量折算为等效电量300MWhcurtail_penalty弃风惩罚单价200元/MWh需要强调单位储氢罐里存的不是 kg而是「等效电量 MWh」也就是罐里那部分氢如果拿去发电最多能出多少电。这样一来电解槽乘效率、燃料电池除效率就能直接加减不用来回换算热值。建模时所有量统一到 MW 和 MWh 上是后面不出现诡异负储氢的前提。2.2 三文件分工sys_init、build_model、plot_results 各管一段代码拆成三个文件边界划得很清楚。sys_init.m只做一件事把上面那张表拼成一个units结构体再生成 24 点负荷曲线和风电预测曲线顺手算出几个派生量比如总装机、M 值的上界、储氢容量对应的最大连续制氢时长。build_model.m负责全部建模和求解输入units输出解结构体sol和一个诊断结构体diag求解状态、耗时、各分项成本。plot_results.m只读不写拿sol画三张图出力堆叠图、储氢量与碳排分档曲线、成本构成饼图。调用骨架就这么几行clear; clc; yalmip(clear); % 清空 YALMIP 内部缓存避免上次的模型残留 units sys_init(); % 生成参数结构体 24 点预测数据 [sol, diag] build_model(units); % 建模、求解、成本拆分 plot_results(sol, units); % 出力曲线 / 氢罐曲线 / 碳档位 if diag.exitflag ~ 1 warning(求解未正常收敛code %d先查第 5 章的排查清单, diag.exitflag); end这样拆的好处是二次开发时改动面可控换数据只动sys_init.m加约束只动build_model.m出图样式完全不碰。diag里我习惯把碳成本、燃料成本、弃风成本、启停成本四项单独存一份因为调试阶段最常见的错误就是总成本算错了却不知道错在哪一项有了分解就能一眼看出是不是碳成本被放大了 1000 倍这种低级事故。2.3 决策变量清单维度错一次后面全乱变量定义写在build_model.m开头维度必须一次定对。火电用T×N_gen的矩阵其余时序量用T×1列向量启停状态用binvar。YALMIP 里最容易翻车的地方是sum的默认行为对矩阵做一次sum(P_gen)得到的是1×3行向量而不是标量直接扔进目标函数会报 nonscalar 错误。T units.T; % 24 N_gen numel(units.gen_cost); % 3 台火电 P_gen sdpvar(T, N_gen, full); % 火电出力 MW列对应机组 P_ele sdpvar(T, 1); % 电解槽耗电功率 MW P_fc sdpvar(T, 1); % 燃料电池发电功率 MW H2 sdpvar(T, 1); % 储氢量等效电量MWh P_curt sdpvar(T, 1); % 弃风功率 MW on_ele binvar(T, 1); % 电解槽开停机状态 0/1 u_ele binvar(T, 1); % 启动动作 0/1 E_seg sdpvar(1, 3); % 阶梯碳排放的三段量见第 3 章配套的基础约束也一并给出机组上下限、弃风上限、氢罐容量这三组是每个时段都要满足的cons []; cons [cons, repmat(units.gen_pmin, T, 1) P_gen repmat(units.gen_pmax, T, 1)]; cons [cons, 0 P_curt P_wind]; % 弃风不能超过可用风资源 cons [cons, units.h2_min H2 units.h2_tank_max]; % 留 20% 安全余量repmat那一行是把1×3的上下限扩成T×3YALMIP 支持这种逐元素双边不等式写起来比循环干净。氢罐下限设成容量的 20% 而不是 0是为了避免模型为了省钱把罐子抽干、第二天开局无氢可用这个细节在第 4 章讲终端约束时还会再提一次。变量都定成连续量还是整数直接决定求解时间只有启停和分段选择用二进制功率和储量一律连续24 点的规模下 CPLEX 通常几秒内出解。3. 阶梯碳交易的建模把 if-else 变成 MILP 的分段线性化3.1 为什么不能直接写 if-else 和 max很多人第一版会写成if emission 200, cost 80*emission; ...这种形式跑起来发现要么报错要么结果完全不对。原因是emission是sdpvar它的值在建模阶段根本不存在MATLAB 的if拿到一个未知对象只能按对象真假走固定分支等于把你的分段逻辑整个跳过了。另一种写法是用min、max拼成本函数比如80*min(E,200) 120*max(0,E-200)读起来漂亮但max是分段线性的非光滑算子YALMIP 会把它解释成非线性算子CPLEX 的 MILP 分支根本不接受求解器直接告诉你模型类型不匹配。正确做法是在建模阶段就把分段结构展开成「若干非负变量 线性等式 上下界」让整个模型保持纯 MILP。这也是为什么这份代码里碳排放量不是一个标量变量而是一个长度为 3 的分段向量。碳排总量本身是线性的由各机组出力和排放强度相乘求和得到加一个固定排放基数比如供热分摊部分便于对齐实际园区的口径% 碳排总量火电出力 × 排放强度逐机组逐时段求和 total_emission sum(P_gen * units.co2_int, all) units.base_emission;注意all这个参数不加的话sum会对矩阵按列求和返回行向量total_emission就不是标量后面所有跟它相关的等式约束都会维度不匹配。这是我在别人代码里见过最频繁的低级错误之一报错信息还特别不友好只说约束维度不一致不说到底哪一行。养成习惯凡是求和得到标量的地方都显式写all或者套两层sum(sum(...))。3.2 三段填充法递增碳价下的最小写法代码包默认用一套三段递增碳价前 200 吨每吨 80 元200 到 500 吨每吨 120 元超过 500 吨每吨 160 元。实现上引入三个非负变量代表落在各档里的排放量加起来等于总排放各档再卡上自己的容量上限成本就是三个量分别乘单价再求和% 阶梯碳交易三档递增价格的分段线性化 co2_seg [200, 500]; % 档位分界点单位 tCO2 co2_price [80, 120, 160]; % 各档单价元/tCO2 E_seg sdpvar(1, 3); % E_seg(1) 落第一档以此类推 cons [cons, sum(E_seg) total_emission]; % 三段之和 总排放 cons [cons, 0 E_seg(1) co2_seg(1)]; % 第一档容量 200 t cons [cons, 0 E_seg(2) co2_seg(2) - co2_seg(1)]; % 第二档容量 300 t cons [cons, E_seg(3) 0]; % 第三档无上限 carbon_cost co2_price(1)*E_seg(1) co2_price(2)*E_seg(2) co2_price(3)*E_seg(3);这段代码里没有二进制变量靠的是两个性质共同成立单价严格递增且目标函数是最小化。求解器为了让成本最小会本能地先把便宜的第一档填满再填第二档绝不会出现「第一档没满却先用第三档」的浪费解所以线性松弛本身就等于整数最优不需要额外加 0-1 变量。一旦政策设计成第二档比第一档便宜某些减免区间或者第三档有封顶返还这个性质立刻失效模型会钻空子把排放全塞进便宜档这时候就必须上二进制变量强制顺序填充。判断标准很朴素把各档单价排成一列看是否单调递增是就用填充法不是就加整数变量。3.3 SOS2、二进制与动态 big-M三种写法的取舍除了填充法工程上还有两种常见写法。一种是特殊有序集 SOS2用一组权重 λ 表示当前解落在哪两个相邻断点之间CPLEX 原生支持这种结构分支效率比二进制高缺点是 YALMIP 里要手工把 λ 和断点绑定写错索引就退化成普通连续变量约束形同虚设。另一种是最直白的二进制选择给每一档配一个 0-1 变量表示「是否进入该档」再用 big-M 把前后档位串起来% 二进制顺序填充保证低档填满才允许用高档 z binvar(3, 1); cons [cons, sum(z) 1]; cons [cons, E_seg(1) co2_seg(1) * z(2)]; % 用第二档则第一档必满 cons [cons, E_seg(2) (co2_seg(2)-co2_seg(1)) * z(3)]; % 用第三档则第二档必满 cons [cons, E_seg(1) co2_seg(1)]; cons [cons, E_seg(2) (co2_seg(2)-co2_seg(1))];这里z(2)1表示已经进入第二档于是强制第一档的下界等于它的容量。真正的坑在 big-M 的取值如果模型里还有别的逻辑用 M 表示「足够大的数」很多人顺手写1e6CPLEX 的对偶边界立刻变得极端病态分支定界在节点上反复震荡跑十几分钟给不出可行解。我的做法是让 M 由上界推出来而不是拍脑袋% 动态 M一天里火电可能产生的最大碳排作为任何分段逻辑的紧上界 M_big T * sum(units.gen_pmax .* units.co2_int) units.base_emission;这个值通常只有几千吨量级跟真实解的差距在两三倍以内数值条件好得多。三种写法的对比如下选哪个主要看碳价是否单调、模型规模有多大写法需要整数变量求解速度适用场景三段填充否最快单价严格递增绝大多数政策场景SOS2 分段求解器内部特殊结构较快断点多于 3 段、需要精确拟合曲线二进制 big-M是每档一个较慢单价非单调、存在减免或返还区间4. 电制氢链路电解槽效率、储氢罐状态方程与启停约束4.1 电解槽的功率-产氢关系二次拟合好看进 MILP 就翻车电解槽的实际效率不是常数低负载时效率明显下滑中高负载区间反而最稳。文献里常见的做法是拟合一条二次曲线比如产氢量关于功率写成0.45*P - 2e-4*P^2这种形式。这条曲线本身没错问题在于它带着平方项一旦写进目标函数或约束整个模型就从 MILP 变成 MIQCPCPLEX 要动用二阶锥或者空间分支定界24 点三机组的小模型还能撑住扩到一周或者加进更多设备求解时间就不是翻倍而是翻十倍。代码包的处理方式是离散化取负载率的三个断点 0、0.5、1.0对应效率 0.60、0.72、0.68用一组权重 λ 和一组段选择变量 z 把「当前工作在哪个区间」表达清楚。这种写法本质上是把曲线掰成两段折线误差控制在 2 个百分点以内但模型保持纯线性% 电解槽效率分段线性化断点 (0, 0.5, 1.0)对应效率 (0.60, 0.72, 0.68) bp [0, 0.5, 1.0]; eta [0.60, 0.72, 0.68]; lam sdpvar(T, 3); % 每个时段一组断点权重 z binvar(T, 2); % 两个线性段的选择变量 for t 1:T cons [cons, sum(lam(t,:)) 1, sum(z(t,:)) 1]; cons [cons, lam(t,1) z(t,1)]; cons [cons, lam(t,2) z(t,1) z(t,2)]; % 相邻断点才允许非零 cons [cons, lam(t,3) z(t,2)]; cons [cons, P_ele(t) units.ele_pmax * (bp * lam(t,:))]; cons [cons, H2_prod(t) units.ele_pmax * ((bp .* eta) * lam(t,:))]; end逻辑说明lam表示当前工作点落在两个相邻断点之间的权重组合z负责把非相邻的权重掐成 0这就保证了「同一时刻只在一段折线上」。P_ele(t)由断点值加权得到H2_prod(t)用断点上的产氢量加权两者用同一组 λ所以功率和产氢严格一一对应。参数上ele_pmax是 40 MW意味着满负载时段产氢约 28.8 MWh 等效电量中途效率峰值在 0.72模型会自然倾向于让电解槽跑在半载以上这也解释了为什么结果图里电解槽很少出现在 20% 负载以下的区间。这套写法比二次曲线的精度差一点但换来了可预测的求解时间工业调度里我觉得这笔交易划算。4.2 储氢罐状态方程索引偏移是负储氢的元凶储氢罐是整条链路里唯一带时间耦合的环节也是 bug 高发区。它的物理含义很直白这一时刻罐里的氢 上一时刻的余量 本时段电解槽制氢 − 本时段燃料电池耗氢。用等效电量表达之后耗氢那项要除以发电效率因为燃料电池发 1 MWh 电需要消耗 1/0.58 ≈ 1.72 MWh 的氢。这个「乘还是除」的方向如果搞反模型会认为燃料电池越发电罐里氢越多跑出来的结果看着收敛实际物理上完全不成立。% 储氢罐动态初值 制氢乘效率 - 耗氢除效率 cons [cons, H2(1) units.h2_init ... units.ele_eff * P_ele(1) - P_fc(1) / units.fuel_cell_eff]; for t 2:T cons [cons, H2(t) H2(t-1) ... units.ele_eff * P_ele(t) - P_fc(t) / units.fuel_cell_eff]; end % 日终闭合末时段储量不低于日初防止吃空库存 cons [cons, H2(T) units.h2_init];这里我把每个时段的功率都按「本时段」写而不是原文里那种「t 时刻储量由 t-1 时刻功率决定」的滞后写法。两种离散方式都能用但滞后写法有个致命副作用最后一个时段的功率不会影响任何储量模型会把最便宜的时段全堆在第 24 点图上看起来末时段制氢功率爆表、罐子却没变化。统一成本时段的写法再补一条日终闭合约束问题就消失了。h2_init一般取容量的 30% 到 40%再配合第 2 章那条 20% 下限罐子全天都在 20% 到 100% 之间游走不会出现负储量。4.3 启停次数与最小出力0.3 这个系数从哪来电解槽是贵设备不能让它每小时在 0 和满功率之间来回横跳。模型里用两个二进制变量表达这件事on_ele表示当前时段是否开机u_ele表示当前时段是否发生了「从停到开」的动作。最小出力约束靠on_ele撑起来启动次数约束靠u_ele累加% 电解槽开停逻辑与最小出力 cons [cons, P_ele units.ele_pmax * on_ele]; % 停机时功率必须为 0 cons [cons, P_ele units.ele_pmin * on_ele]; % 开机则不低于 12 MW cons [cons, u_ele(2:T) on_ele(2:T) - on_ele(1:T-1)]; % 状态跃迁即一次启动 cons [cons, u_ele(1) on_ele(1) - units.on_init]; % 首时段也要算启动 cons [cons, sum(u_ele) units.max_start]; % 全天不超过 3 次ele_pmin取 40 MW 的 0.3 倍也就是 12 MW这个系数不是拍脑袋来的低于 30% 负载时电解效率会明显掉到 0.5 以下单位氢的耗电量飙升与其这样不如让设备停着。max_start默认 3是设备厂家给的日均启停次数建议值。u_ele(1)那行很容易漏漏了之后模型会把「凌晨先开一次机」当成免费的启动成本少算一次。启动成本在目标函数里按次计价通常设成 200 到 500 元一次具体看设备容量数值不需要很准它的作用是把频繁启停的收益压下去。5. 求解配置与踩坑排查从 MIP gap 到负储氢的五条血泪记录5.1 目标函数组装与 CPLEX 选项分项成本必须能拆开看目标函数一共四项燃料成本、运维成本、弃风惩罚、阶梯碳成本再加上电解槽启动成本这一小块。写的时候逐项算完再相加不要图省事写成一个巨大的表达式否则调参时根本看不出是哪一项在主导结果。% 分项成本每一项单独成变量方便事后拆解 fuel_cost sum(sum(P_gen .* repmat(units.gen_cost, T, 1))); om_cost units.om_coef * sum(P_gen, all) units.h2_om * sum(P_ele); curt_cost units.curtail_penalty * sum(P_curt); start_cost units.start_cost * sum(u_ele); total_cost fuel_cost om_cost curt_cost start_cost carbon_cost; % 求解器选项MILP 三条常用设置 ops sdpsettings(solver, cplex, verbose, 1); ops.cplex.timelimit 600; ops.cplex.mip.tolerances.mipgap 0.01; ops.cplex.mip.tolerances.absmipgap 50; sol optimize(cons, total_cost, ops); diag.exitflag sol.problem; diag.fuel value(fuel_cost); diag.carbon value(carbon_cost);参数说明mipgap 0.01表示允许 1% 的相对最优间隙调度问题里这个精度完全够用实际总成本几百万元的时候误差几千元低于风电预测误差带来的不确定性absmipgap 50是绝对间隙作用是当目标函数量级很小的时候不让求解器为了 1% 死磕timelimit必须设没有时间上限的 MILP 在极端数据下可能跑一整夜。verbose 1会打印分支节点数、当前最好解和下界调试的时候看这两个值的差距缩得快不快比盯着进度条有用得多。5.2 五条踩坑记录现象、原因、解决现象一结果里储氢量出现负值。原因几乎总是效率折算方向写反或者状态方程的时序错位导致末时段没闭合。解决方式是先检查H2(t)是不是被h2_min这类下界约束真正包住了再核对制氢乘效率、耗氢除效率的方向如果模型仍然给出负值把H2的下界临时收紧到 0 跑一遍对比两次结果就能定位是约束漏写还是求解精度问题。现象二模型显示 infeasible但怎么看参数都合理。最常见的原因是储氢容量和制氢功率不匹配40 MW 电解槽满发一小时产 28.8 MWh 等效氢如果罐子只给了 30 MWh 容量早上四个小时连发就把罐子撑爆模型只能报不可行。解决办法是把h2_tank_max和ele_pmax联动设计经验上罐容至少能装下电解槽满发 6 到 8 小时的产量或者干脆给弃风留出更大空间。现象三求解时间从几秒变成十几分钟。十有八九是某处 big-M 取了1e6这类夸张值或者碳排总量用了一个很宽的人工上界。把 M 换成由设备容量推出来的紧上界之后同一个模型经常在一分钟内收敛。这条血泪经验我踩过不止一次后来养成习惯模型里每出现一个 M都要在注释里写清它是怎么算出来的。现象四碳成本高得离谱比燃料成本还大一个数量级。单位混用。有的文献碳价按元/kg 写有的按元/t 写两者差 1000 倍还有些参数表里排放强度给的是 kgCO2/MWh 而代码按 t 处理。排查方法很简单把value(total_emission)打出来看量级是不是在几百到几千吨如果只有零点几说明单位错了一千倍。现象五电解槽全天频繁启停启动次数约束好像没生效。检查u_ele的定义很多人只写了u_ele(t) on_ele(t) - on_ele(t-1)却忘了首时段也忘了把u_ele放进目标函数。如果启动成本系数是 0模型会认为启停免费约束虽然限制了次数上限但 3 次之内随便跳图上就显得很乱。补上首时段约束、把启动成本调到几百元每次曲线立刻变干净。6. 二次开发从单日 24 点扩到多日滚动我一般这么改单日模型最大的局限是储氢罐必须在当天闭合跨天的能量搬移被硬生生切断。要扩成一周或者一个月的滚动调度改动集中在外层把build_model.m包进一个按天推进的循环每天调用一次把前一天的末态作为当天的初值传进去同时把日终约束从「不低于日初」改成「不低于某个保底值」只在最后一天强制闭合。% 多日滚动末态承接只有最后一天强制闭合 h2_state units.h2_init; for d 1:N_day units.h2_init h2_state; units.is_last (d N_day); % 最后一天才要求罐子回到初始水平 [sol, diag] build_model(units); h2_state value(sol.H2(end)); rec.cost(d) diag.total; end改成滚动之后验证方法也要跟着变。我一般做三件事一是拿固定碳价版本跑同一组数据做基准看阶梯碳价下的总成本是不是落在 8% 到 15% 的下降区间超出这个范围多半是参数写错二是检查能量守恒残差把每个时段的功率平衡左右两边相减取绝对值最大值应该小于 1e-6否则就是有约束没写全三是把碳档位曲线画出来看火电出力在跨档时刻有没有出现预期的拐点没有拐点说明碳价传导到调度决策的链路断了。再往深一层可以加碳捕集设备把排放压回第一档或者把风电预测换成多场景鲁棒版本。鲁棒版不用重写模型只要把P_wind从常量改成场景集合再把功率平衡约束对每个场景各写一份用binvar选最坏场景即可代价是二进制变量翻倍求解时间也会跟着上去。改造氢能部分时我每次都会先把储氢状态方程单独拿出来用固定功率输入跑一遍开环仿真确认罐子曲线单调合理了再放进优化模型不然一旦出现负储量你根本分不清是物理模型错了还是求解器给了个奇怪的解。这套排查顺序希望帮到你。本文还有配套的精品资源点击获取