MATLAB+YALMIP阶梯碳交易与电制氢调度建模

简介：这份PDF面向能源系统优化研究者及具备MATLAB/YALMIP基础的技术人员，聚焦含阶梯型碳交易与电制氢模块的综合能源系统优化调度问题。内容围绕传统机组、风电、电解槽、储氢罐和燃料电池五类单元展开，讲解阶梯碳成本分段线性化与big-M处理、氢能系统状态方程、目标函数设计及案例结果可视化，并给出完整代码实现思路与文献参考。压缩包仅含1个PDF文件，大小约129KB，便于集中阅读和存档。已有163人学习，适合用来理解电—氢—电多向转换、弃风制氢与碳排阈值附近火电调整等建模细节，也可为搭建MILP调度模型、处理时间耦合约束和调试数值稳定性提供参考。

1. 凌晨两点的那条曲线，让我确认了模型没写错

有人发来一张 24 小时出力图问我，为什么火电曲线在某个时段突然往下压了一截，是不是约束写漏了。我看了十几分钟，回他：没写漏，那是阶梯碳价在起作用。碳排放一旦跨过分档边界，边际成本从 80 元一吨跳到 120 元，优化器自己就会把火电往下压，这个拐点恰恰说明模型是活的。

这份 MATLAB 程序干的就是这一件事：用 YALMIP 搭模型、CPLEX 求解，把传统机组、风电、电解槽、储氢罐、燃料电池五类单元塞进同一个 24 小时调度框架，目标函数里既有燃料成本和弃风惩罚，也有按阶梯计价的碳成本。它适合已经会写基本机组组合、想往综合能源系统方向迈一步的人，也适合手上有个调度项目临时要加碳约束、正被分段线性化卡住的开发者。阶梯碳交易和电制氢这两块，一个卡在"非线性"，一个卡在"时间耦合"，恰好是这类模型最容易翻车的两个位置。

2. 模型骨架与文件分工：五类设备、四个文件怎么拼成一次求解

拿到一份 MATLAB 优化代码，我第一件事不是看目标函数，而是看文件怎么切。切得乱的项目，后面加设备、加约束基本只能推倒重来。这份代码的切法比较规矩，四个文件各管一件事，改了参数不用碰模型，改了模型不用碰画图。

这种切法不新鲜，但真正做到"参数和模型隔离"的项目不多。我见过太多脚本把所有东西写在一个 600 行的.m文件里，clear all一跑到底，客户想换个碳价得在一堆约束里找常量。

2.1 sys_init.m：参数集中营

sys_init.m返回一个结构体，把机组成本、装机容量、效率、碳价档位全部收进去。写法上我习惯用结构体而不是散落的工作区变量，原因是build_model.m只依赖这一个入参，函数签名清清楚楚。

function units = sys_init() % 系统参数初始化：所有可调参数集中在这里，换场景只改这个文件 units.T = 24; % 调度时段数，1 小时一步 units.gen_cost = [28 31 25]; % 三台火电机组燃料成本系数，元/MWh units.gen_pmax = [200 150 120]; % 机组出力上限，MW units.gen_pmin = [60 50 40]; % 机组出力下限，MW units.gen_emiss = [0.95 1.05 1.20]; % 单位出力碳排放强度，tCO2/MWh units.wind_cap = 150; % 风电装机容量，MW units.wind_curve = [ ... ]; % 24 点风电可用出力，MW units.load_curve = [ ... ]; % 24 点电负荷，MW units.ele_eff = 0.72; % 电解槽额定效率 units.ele_pmax = 60; % 电解槽额定功率，MW units.ele_pmin_ratio = 0.30; % 最低负载率，低于此值效率骤降 units.fc_eff = 0.58; % 燃料电池发电效率 units.fc_pmax = 30; % 燃料电池功率上限，MW units.h2_cap = 3000; % 储氢罐容量上限，kg units.h2_init = 600; % 初始储氢量，kg units.h2_min_ratio = 0.15; % 储氢罐最低安全储量比例 units.carbon_tier = [2000 4000]; % 碳排放分档边界，tCO2 units.carbon_price = [80 120 160]; % 各档边际碳价，元/tCO2 units.curtail_pen = 200; % 弃风惩罚系数，元/MWh end

参数说明上有个容易忽略的细节：ele_pmin_ratio单独抽出来做了变量。很多代码直接把这 0.3 硬编码在约束表达式中，等到客户说"这台设备最低负载率是 0.25"的时候，得全局搜索 0.3，一搜几十处，全是血泪经验。凡是带物理含义的系数，一律进结构体。

carbon_tier和carbon_price一起定义了阶梯结构：0 到 2000 吨按 80 元，2000 到 4000 吨按 120 元，超过 4000 吨部分按 160 元。注意这里存的是边际价而不是累进总价，两者的差别在下一章会展开，写混了碳成本能差出一倍。

2.2 build_model.m 的三段式结构

模型构建文件按"变量、约束、目标"三段排，中间用注释分隔，这是 YALMIP 脚本最经得起改的结构。三段之间不该有交叉依赖，否则加一条约束就得回头看变量定义在哪。

function [ops, cons, obj, vars] = build_model(units) % 变量定义：维度统一按 (T, N) 排，避免后面到处转置 T = units.T; Ngen = numel(units.gen_cost); P_gen = sdpvar(T, Ngen, 'full'); % 火电机组出力，MW u_gen = binvar(T, Ngen, 'full'); % 机组启停状态，1 表示出力大于 0 P_wind = sdpvar(T, 1); % 风电实际消纳，MW P_curt = sdpvar(T, 1); % 弃风量，MW P_elec = sdpvar(T, 1); % 电解槽用电功率，MW u_elec = binvar(T, 1); % 电解槽启停状态 H2_prod = sdpvar(T, 1); % 制氢量，kg H2_use = sdpvar(T, 1); % 燃料电池耗氢量，kg H2_soc = sdpvar(T, 1); % 储氢罐存量，kg P_fc = sdpvar(T, 1); % 燃料电池发电功率，MW E_seg = sdpvar(1, 3); % 三档碳排放量分解，tCO2 lam = sdpvar(T, 3); % 电解槽效率曲线的凸组合权重 vars = struct('P_gen', P_gen, 'u_gen', u_gen, 'P_wind', P_wind, ... 'P_curt', P_curt, 'P_elec', P_elec, 'H2_soc', H2_soc, ... 'P_fc', P_fc, 'E_seg', E_seg);

维度约定值得单独说一句：P_gen用(T, Ngen)而不是(Ngen, T)。原因很简单，24 小时是行，机组是列，画图的时候plot(value(P_gen))直接出多条时间曲线，sum(P_gen,1)是全时段各机组发电量，sum(P_gen,2)是各时段总出力。反过来排的话，代码里到处是'转置，多写几次就会有一处漏掉，然后报维度不匹配——这种错不致命，但一晚上能耗进去。

变量定义完成之后是约束段，按设备分块：

cons = []; % ---- 火电机组：出力上下限 + 启停逻辑 ---- for i = 1:Ngen cons = [cons, units.gen_pmin(i)*u_gen(:,i) <= P_gen(:,i) <= units.gen_pmax(i)*u_gen(:,i)]; end cons = [cons, P_gen(1,:) == 0.5*(units.gen_pmin + units.gen_pmax)]; % 首时段给定初值，防退化 % ---- 风电消纳与弃风 ---- cons = [cons, P_wind + P_curt == units.wind_curve]; cons = [cons, P_wind >= 0, P_curt >= 0]; % ---- 功率平衡：源侧等于荷侧加电解槽 ---- cons = [cons, sum(P_gen, 2) + P_wind + P_fc == units.load_curve + P_elec];

功率平衡这一行是整个模型的骨架。sum(P_gen,2)是各时段的火电总出力，加上风电消纳、燃料电池发电，等于负荷加上电解槽吃掉的电。这里把电解槽当作负荷挂在平衡式右边，是电制氢模型最自然的写法——电解槽对电网而言就是个可控负荷。

u_gen这个二进制变量容易被省略。只在 pmin 和 pmax 之间做约束、不要启停变量，模型也能跑通，解出来也是可行的，但机组会在 pmin 上挂着运行一整晚。实际调度里机组停机是常态，加上u_gen之后模型才能表达"这台机组今天不发电"。代价是多了 24×3 个二进制变量，MILP 求解时间会涨，这个取舍在后面讲求解器配置时还会遇到。

2.3 目标函数与结果自检

目标函数分四块拼，写法上尽量拆成独立的中间变量，方便出问题的时候逐项打印。

% 燃料成本：gen_cost 是行向量，需要转置成 (1,Ngen) 才能和 (T,Ngen) 对齐 fuel_cost = sum(sum(P_gen .* units.gen_cost')); % 弃风惩罚 curtail_cost = units.curtail_pen * sum(P_curt); % 氢能运维：制氢按 kg 计、燃料电池发电按 MWh 计 h2_om = 12 * sum(H2_prod) + 8 * sum(P_fc); % 阶梯碳成本（具体构造见第 3 章） carbon_cost = units.carbon_price * E_seg'; obj = fuel_cost + curtail_cost + h2_om + carbon_cost; ops = sdpsettings('solver', 'cplex', 'verbose', 2);

P_gen .* units.gen_cost'这一行的转置是必须的。units.gen_cost是1×3行向量，P_gen是24×3，直接点乘会触发 MATLAB 的隐式扩展——从 R2016b 起这个行为是合法的，结果也正确，但很多人的环境还是老版本，跑出来直接报维度不匹配。我一般会把参数在sys_init.m里就存成列向量，从源头掐掉这个问题。

plot_results.m里除了常规的堆叠面积图，我建议加一段自检代码。优化模型最怕的不是报错，是安安静静跑出一个看起来合理的错误结果。

% 结果自检：分段碳排放量之和必须等于实际排放量 E_real = sum(sum(value(P_gen) .* units.gen_emiss')); E_seg_v = sum(value(E_seg)); assert(abs(E_real - E_seg_v) < 1e-3, ... '碳排放分段量与实际排放不一致：%.2f vs %.2f', E_seg_v, E_real); % 功率平衡残差检查 resid = sum(value(P_gen),2) + value(P_wind) + value(P_fc) ... - units.load_curve - value(P_elec); assert(max(abs(resid)) < 1e-3, '功率平衡存在残差，检查约束是否被覆盖');

这两段断言平时不显眼，出问题的时候能省掉一整个下午。尤其是第一条，如果碳成本构造写错（比如只算了第一档），模型照样收敛、图照样好看，只有自检能把它拦住。

提示：assert放在run_main.m里，而不是build_model.m里。建模函数不该依赖求解结果，职责要分清。

3. 阶梯碳交易建模：把 if-else 改成分段线性约束

阶梯碳价是这套模型里最容易被做错的部分。做错的方式主要有三种：拿if-else直接对着sdpvar写、用min/max套在决策变量上、把边际价当累进价用。这三种写法里前两种会让模型跑不动或者跑出线性化错误，第三种模型能跑但结果是错的，更危险。

3.1 三段边际计价与辅助变量法

先把计价规则说清楚。阶梯电价那种算法是"分段累进"：前 200 度按 0.5 元，超过部分按 0.8 元，这是对总电量分段。碳交易这里用的是同一套逻辑，对总排放量分段，每段的单价是边际价。

用辅助变量把总排放量拆成三段：

% 总排放量由火电出力决定，先算出来 E_total = sum(sum(P_gen .* units.gen_emiss')); % 标量 sdpvar % 三档分解，E_seg = [E1 E2 E3] E1 = E_seg(1); E2 = E_seg(2); E3 = E_seg(3); t1 = units.carbon_tier(1); % 2000 t2 = units.carbon_tier(2); % 4000 p = units.carbon_price; % [80 120 160] cons = [cons, E1 + E2 + E3 == E_total]; % 分解守恒，这一条最关键 cons = [cons, 0 <= E1 <= t1]; % 第一档最多 2000 吨 cons = [cons, 0 <= E2 <= t2 - t1]; % 第二档最多 2000 吨 cons = [cons, E3 >= 0]; % 第三档无上限 carbon_cost = p(1)*E1 + p(2)*E2 + p(3)*E3;

四行约束加一行成本表达式，整个阶梯逻辑就表达完了。这里没有一处if，也没有min/max，进模型的是标准 LP 结构。

参数含义再明确一遍：t1、t2是排放量的档位边界，单位 tCO2；p是各档的边际碳价，单位元/tCO2。carbon_cost算出来的是总碳成本，单位元。三档累加就等于"前 2000 吨 × 80 + 中间 2000 吨 × 120 + 超出部分 × 160"，与分段函数完全等价。

很多人第一次看这段会觉得少了点什么——为什么不用限制E2必须填满才能填E3？答案是目标函数是最小化方向。单价递增的情况下，优化器会本能地把便宜的档位先填满，绝不可能出现E1 = 100、E3 = 500这种"跳过便宜档直接买贵的"的解，因为那样成本更高。这一点在下一节展开。

3.2 为什么价格递增时不需要二进制变量

分段线性化有个经典权衡：用辅助变量法省二进制、求解快，但只在"目标函数方向与分段单调性一致"时才正确；用 big-M 加二进制变量，逻辑严格但求解慢。这份代码选的是前者。

判断标准很干脆：如果各档边际价单调递增，且目标是最小化，就可以省掉二进制变量。原因在于，假设某个可行解没有填满低价档就动用了高价档，那么把高价档的量挪一部分回低价档，总排放量不变（分解守恒约束仍然满足），成本一定下降。所以任何非"低价优先"的解都是劣解，最优解自动满足顺序填充。

这个性质一旦被打破，辅助变量法就失效了。两种常见情况：

一是碳价倒挂，比如第二档因政策优惠反而比第一档便宜。这时候优化器会故意把排放量塞进第二档，成本算出来偏低，但和现实不符。

二是计价规则不是边际累进，而是"整段跳价"——比如超过 4000 吨之后，全部排放量都按 160 元重算。这种规则是非凸的，必须上二进制变量：

% 整段跳价模型：用二进制变量 z(k) 表示落在第 k 档 z = binvar(3, 1); cons = [cons, sum(z) == 1]; cons = [cons, E_total <= t1*z(1) + t2*z(2) + 1e5*z(3)]; cons = [cons, E_total >= t1*z(2) + t2*z(3)]; % 档位下界 carbon_cost = E_total * (p(1)*z(1) + p(2)*z(2) + p(3)*z(3));

第三种容易踩的坑是拿min/max写成本：

% 错误写法：对 sdpvar 用 min/max carbon_cost = 80*min(E_total, 2000) + 120*max(0, E_total - 2000);

这段代码在 YALMIP 里看起来能跑，实际上会触发非线性算子处理，求解器要么直接报不支持，要么被转成一堆大 M 约束，最后慢得莫名其妙。min/max是给数值计算用的，不是给决策变量用的。

3.3 档位参数怎么改，改完看什么

阶梯参数经常要跟着政策口径调，改动位置在sys_init.m里，改完有三处要核对。

档数改动要同步改三处：sys_init.m里的两个数组、build_model.m里E_seg的维度和上下限约束、plot_results.m里成本分解的标签。我一般会写个小循环根据carbon_tier自动生成长度，而不是手写E1、E2、E3三个变量名，这样从三档变四档只改一行数组。

% 通用写法：档数由参数决定，加档不用改模型 K = numel(units.carbon_price); E_seg = sdpvar(1, K); edges = [0, units.carbon_tier, inf]; % 长度 K+1 for k = 1:K cons = [cons, 0 <= E_seg(k) <= edges(k+1) - edges(k)]; end cons = [cons, sum(E_seg) == E_total]; carbon_cost = units.carbon_price * E_seg';

改完参数之后别急着看成本降了多少，先看拐点位置。如果碳价提高之后火电曲线完全没变化，八成是carbon_price的单位和gen_emiss对不上——排放强度按 tCO2/MWh 存，碳价按元/t 存，乘出来是元/MWh，和燃料成本的量级要能对得上。我见过一次把碳价写成元/kg 的，成本里碳的部分直接贡献了九成，火电全时段压到最低，图看起来"效果很好"，实际上是量纲错了。

4. 电制氢链路：电解槽、储氢罐、燃料电池的时间耦合

电制氢这部分比碳交易更考验建模细心程度。碳交易只是单个标量的分段，氢能是一条带时间耦合的链：电进电解槽、产出氢、存进罐子、罐子放氢给燃料电池、燃料电池再发电。任何一个环节的索引或效率折算写反，结果都会以"储能曲线诡异"的形式暴露出来。

4.1 储氢状态方程与索引对齐

状态方程本身不复杂，写成一行就行：

% 储氢罐动态：本期存量 = 上期存量 + 本期制氢 - 本期耗氢 for t = 2:T cons = [cons, H2_soc(t) == H2_soc(t-1) + H2_prod(t) - H2_use(t)]; end

注意这里H2_prod(t)和H2_soc(t)的时段下标是一致的：第 t 时段制出来的氢，计入第 t 时段结束时的存量。写成H2_prod(t-1)就变成"上一时段的氢存到这一时段"，模型也能解出来，但整条曲线会整体后移一小时，储氢罐的响应看起来总是慢一拍。

首末时段的处理是另一个高频翻车点：

% 第 1 时段：接入初始存量 cons = [cons, H2_soc(1) == units.h2_init + H2_prod(1) - H2_use(1)]; % 第 24 时段：日末存量不低于初始值，防止模型"把罐子掏空换成本" cons = [cons, H2_soc(T) >= units.h2_init]; % 全时段容量与安全下限 cons = [cons, units.h2_min_ratio*units.h2_cap <= H2_soc <= units.h2_cap];

末端的H2_soc(T) >= h2_init这条约束看着不起眼，删掉之后总成本能降一大截，而且解完全可行——模型会在最后几个时段把罐子里的氢全拿去发电，等于白赚了一批存量。这种解在数学上没错，在物理上也没错，但在调度业务上是不接受的：今天把库存清空，明天开机就没氢用了。多日滚动调度的时候，这条约束必须保留，否则逐日优化会累积出明显的偏差。

如果要做多日连续调度，末端处理要换成软约束或者跨日传递：

% 多日调度：把最后一天的末端存量作为下一天的初始值传出去 h2_end = value(H2_soc(T)); save('h2_state.mat', 'h2_end'); % 下一天调用时：units.h2_init = h2_end;

4.2 电解效率分段线性化：别让 eta 和 P 相乘

电解槽的效率不是常数，负载率低的时候效率掉得厉害。常见做法是给一条三点曲线，负载率 0.3 对应效率 0.60，0.6 对应 0.72，满负载对应 0.66。制氢量等于效率 × 功率 × 时间，这里效率是决策变量的函数，直接乘就变成双线性，MILP 求解器处理不了。

第一个坑就出在这：很多人先线性化效率，再把它和功率相乘。

% 错误做法：eta 和 P_e 都依赖 lam，相乘后是双线性 P_e = (bp .* lam)' * units.ele_pmax; eta = (eta_bp .* lam)'; % eta 也是 lam 的线性组合 H2 = eta * P_e * dt / 1000; % 两个 lam 相乘，非线性

正确做法是在 (功率, 产氢量) 平面上取断点，让两者都成为同一组权重的凸组合，乘积关系被吸收进断点坐标里：

% 断点：负载率、对应效率 bp = [0.30 0.60 1.00]; eta = [0.60 0.72 0.66]; % 先在 (P, H2) 平面上算出三个断点坐标 P_bp = bp * units.ele_pmax; % MW H_bp = eta .* P_bp * 1 / 1000; % t/h，1 小时一步 % 凸组合权重，SOS2 保证只有相邻两个权重非零 lam = sdpvar(T, 3); for t = 1:T cons = [cons, sum(lam(t,:)) == 1, lam(t,:) >= 0]; cons = [cons, sos2(lam(t,:))]; cons = [cons, P_elec(t) == P_bp * lam(t,:)']; cons = [cons, H2_prod(t) == H_bp * lam(t,:)']; end

sos2是特殊有序集约束，语义是"至多两个相邻元素非零"，CPLEX 原生支持。它把凸组合限制在相邻断点之间，插值结果就落在那条折线段上，不会跑到三角形内部去。如果不加sos2，模型会在三个断点之间做凸包，效率曲线被"拉直"，制氢量被高估。

H_bp的单位要盯紧。功率是 MW，时间是 1 小时，效率无量纲，乘出来是 MWh，除以 1000 换成吨。储氢罐容量按 kg 存的话，这里就别除 1000，保持量纲一致比什么都重要。我在一个项目里因为储氢量一半用 kg 一半用 t，导致容量约束形同虚设，罐子能存设计容量的 1000 倍，求解器毫无怨言。

4.3 启停约束与双向效率折算

电解槽频繁启停会伤设备，所以加一条日启停次数约束：

% 24 小时内启停不超过 3 次 cons = [cons, sum(abs(diff(u_elec))) <= 6]; % 一次启停涉及两次状态跳变

diff和abs直接作用在决策变量上同样会引入非线性。稳妥写法是引入辅助变量：

d = binvar(T-1, 1); % 状态变化指示 for t = 1:T-1 cons = [cons, d(t) >= u_elec(t+1) - u_elec(t)]; cons = [cons, d(t) >= u_elec(t) - u_elec(t+1)]; end cons = [cons, sum(d) <= 6]; % 总计最多 6 次跳变，约合 3 次启停

功率上下限和状态变量的绑定也要写全：

% 运行下限按额定功率的比例给，低于 0.3 倍额定效率掉太快 cons = [cons, units.ele_pmin_ratio*units.ele_pmax*u_elec <= P_elec]; cons = [cons, P_elec <= units.ele_pmax*u_elec];

燃料电池侧是反过来的：给定发电功率，耗氢量要除以效率，而不是乘。

% 燃料电池：发 P_fc 的电，需要消耗 P_fc/eta 的氢能 cons = [cons, H2_use == P_fc ./ units.fc_eff * 1 / 1000]; % t/h cons = [cons, 0 <= P_fc <= units.fc_pmax];

这一乘一除是整个氢能链路最容易写反的地方。电解槽耗电产氢，效率乘在功率上；燃料电池耗氢发电，效率除在功率上。写反之后模型仍然可解，但会出现"充进去 1 度电、放出来 1.7 度电"的永动机式结果，储氢曲线一路向上不回头。我第一版就栽在这，盯着图看了半小时才反应过来效率位置错了。

注意：如果氢能还要外送或者供给氢负荷，H2_use那一项就要拆成"发电用氢"和"外送氢"两部分，状态方程里减去的是两者之和。预留这个拆分接口，后面接氢负荷不用改结构。

5. 常见问题排查：五条从报错日志里抄下来的记录

模型能跑通和模型能用之间隔着一堆细节。下面这五条是我在这套结构里反复遇到的，每条都按现象、原因、解决排一遍，遇到同类问题可以直接对号。

5.1 求解层面的三类异常

第一条：求解器返回 infeasible，但约束看起来都对。现象是 CPLEX 报Infeasible，diagnostics一看，冲突集中在功率平衡那一行。原因通常是风电可用出力曲线和负荷曲线的量级不匹配，比如风电存的是标幺值（0～1），负荷存的是 MW（几百），平衡式永远成立不了。解决办法是先做一次量纲体检：把wind_curve、load_curve、各设备pmax打印出来看数量级，再逐条删约束做二分定位。YALMIP 的optimize带'debug'选项能导出冲突集，比重头读代码快得多。

第二条：解出来了，但 gap 一直是 20% 以上，跑一小时也不收敛。现象是终端刷屏，verbose=2下能看到 incumbent 和 best bound 差距巨大。原因多半是 big-M 取得太大，LP 松弛极弱。手写M = 1e6是最常见的做法，也是最要命的做法。正确姿势是让 M 等于约束里那个变量的物理上界：启停约束P_elec <= P_max*u_elec里的 M 就是P_max本身，写成 1e6 等于告诉求解器"这个变量可以到一百万"，分支定界的下界自然松得离谱。把 M 逐个替换成物理上界之后，同一个模型求解时间从几十分钟掉到几十秒。

第三条：换了一组碳价参数，模型突然报数值错误或者返回 NaN。现象是optimize返回Numerical problems，或者解出来的成本是 NaN。原因是分段变量的上下界和档位参数之间出现了矛盾，比如把carbon_tier改成递减数组[4000 2000]，第二档的t2 - t1变成负数，约束直接不可行。解决办法是在sys_init.m末尾加一行校验：

assert(all(diff([0 units.carbon_tier]) > 0), '碳档位边界必须严格递增'); assert(all(diff(units.carbon_price) >= 0), '碳价需单调不减，否则分段线性化失效');

这两行断言花不了五秒钟，能拦住一整类问题。

5.2 建模逻辑层面的两个隐形错误

第四条：储氢量出现负值，或者曲线一路单调向上。现象是H2_soc里有负数，或者整条曲线从早到晚只涨不跌。负值一般是末端约束或容量下界漏写，导致模型可以在某些时段"借"氢。单调向上则基本可以确定是效率折算写反了，燃料电池那侧用了乘法而不是除法，等于凭空造能量。解决办法是先检查H2_use == P_fc/eta/1000这一行的除号，再检查H2_soc是否同时有上下界约束。加一条cons = [cons, H2_soc >= 0]能让问题立刻暴露成 infeasible 而不是悄悄跑出错误结果。

第五条：碳成本算出来了，但和手工核算对不上。现象是成本分解饼图里碳的部分和按排放量手算的差了一截。原因有三种：分段变量之和不等于总排放（漏了守恒约束）、碳价用了累进总价而不是边际价、排放强度矩阵转置错了。解决办法就是第 2 章那段自检代码，把E_real和sum(E_seg)对比，差值超过 1e-3 就断言失败。另外一个笨办法但很有效：把最优解代回原始分段函数，用if-else手算一遍碳成本，和模型输出的value(carbon_cost)逐位对比。前一种快，后一种能查出逻辑错。我一般两个都做。

提示：plot_results.m里把碳成本的三个分档量分别画成柱状堆叠图，一旦某一档长期为 0 或者长期顶格，都能一眼看出来是参数问题还是模型问题。

6. 二次开发：从 24 小时单日模型改成能交付的调度脚本

单日模型跑通只是起点。真要做成能交给别人用的东西，第一件事是把单日循环改成多日滚动。外层加一个日期循环，每一天读当天的风光负荷预测，把前一天的末端储氢量作为当天的初始值传进去，跑完再把结果写回。这套结构改起来不复杂，但有个坑必须提前处理：yalmip('clear')要放在循环内部，不放在外面。放到外面的话，第二天的模型会带着第一天的变量一起求解，规模翻几十倍，而且约束互相污染。我吃过这个亏，第二天直接跑到内存溢出。

第二个改动点是求解器配置。单日模型 gap 设 1% 够用，多日滚动里如果每天都等 1%，跑七天就是一小时。常见的处理方式是首日设严一点（0.5%），后续日期放宽到 2%～3%，同时给每天设一个时间上限，超时就用当前可行解。工程上没人关心第 4 天的最优解比次优解便宜 3 块钱。

第三个改动点是目标函数的层次。现在的写法是加权求和，碳成本、弃风惩罚、氢能运维各给一个系数。真要交付的时候，客户往往问"如果我想把碳排放再压 10%，成本要多花多少"。这时候加权求和就不够用了，得改成 ε-约束：固定一个排放上限，最小化成本，把上限从紧到松扫一遍，画出一条权衡曲线。代码上只是加一条cons = [cons, E_total <= eps_k]，外层套一个eps_k的循环。这个改动很小，但能回答的问题完全不一样，交付时说服力强很多。

最后一个建议关于数据接口。现在曲线是硬编码在sys_init.m里的定长数组，换数据要手改。我一般会把它改成从csv读，加一层校验：长度是否等于T、是否有负值、风电可用出力是否超过装机。这些检查写在sys_init.m里，比在模型里出问题之后再回头查要省事得多。

从那以后我每次改完参数，都会先跑一遍自检断言，再看图和成本分解，最后才把结果发出去。这套流程走下来多花五分钟，但能把"图看着对、数是错的"这种最难查的问题挡在门外。希望帮到你。

本文还有配套的精品资源，点击获取