简介：面向从事综合能源系统优化、且具备一定MATLAB基础的科研与技术人员，这份PDF围绕阶梯型碳交易与电制氢模块的协同调度展开，重点解决传统机组、风电、电解槽、储氢罐和燃料电池五单元联合运行时的经济性与碳排平衡问题。包内共1个PDF文件，大小约129KB，以图文与代码片段结合的方式呈现系统建模、阶梯碳成本分段处理、氢能状态方程、目标函数设计及案例可视化结果。已有162人学习，适合用来理解风电大发时段优先制氢、碳排放接近阈值时火电出力调整等现象。读者可借此梳理YALMIP/Cplex求解混合优化模型的实现思路，掌握分段线性化、时间耦合约束和结果分析的关键细节，为综合能源系统调度研究与方案复现提供参考。
1. 一个能跑出解、但储氢量是负数的调度模型
上周帮 A 同学看一个 MATLAB + YALMIP 的综合能源系统优化调度模型，加进阶梯碳交易和电制氢模块之后，CPLEX 不再报错，求解状态显示成功，可画出来的储氢曲线在凌晨三点掉到 −80 kg，燃料电池偏偏在风电满发的时段满功率发电。这类“能求解、结果却是错的”故障比直接报错难缠得多，因为它不给你任何提示，线索全埋在约束的索引和分段逻辑里。阶梯碳价随手写成 if-else 会把模型变成非凸，电制氢的双向效率折算写反会让氢气凭空产生，储氢罐递推少算一拍会让末端状态不闭合——这三处基本覆盖了这类模型九成以上的失败现场。下面按报错现象、根因、可以直接照抄的改法依次拆开，MATLAB + YALMIP + CPLEX 这套组合的参数怎么调、结果怎么验，都落到具体代码上。
2. YALMIP 报 Non-convexity detected：阶梯碳价写成 if-else 的三种后果
2.1 报错长什么样，哪几种写法必然触发
先看现场。用分段线性思路写碳成本，很多人第一版是这样：
% 错误写法 A：直接对 sdpvar 用 min / max carbon_cost = 80*min(total_emission,200) + 120*max(0,total_emission-200);求解时 YALMIP 会抛出一串警告，大意是检测到非凸结构，随后要么退化成一个很慢的分支流程，要么干脆把 min 项当成常数处理。根因在于min是凹函数，最小化一个凹函数不满足凸性，CPLEX 的 LP/MILP 引擎吃不下，YALMIP 只能尝试自动引入二进制变量做 big-M 展开，而自动展开的 M 值往往取不到合适的量级，最后报CPLEX Error 5002: Objective is not convex或者长时间卡在根节点。
更隐蔽的是错误写法 B——把碳成本单独写成普通 MATLAB 函数：
% 错误写法 B：用 if-else 写分段，对符号变量完全不成立 function c = calc_carbon_cost(E) if E <= 200 c = 80*E; else c = 80*200 + 120*(E-200); end end carbon_cost = calc_carbon_cost(total_emission); % total_emission 是 sdpvarsdpvar进不了if的逻辑判断：MATLAB 会调用重载的比较运算，返回一个约束对象，if拿到非逻辑值后要么直接报“无法转换为逻辑值”，要么在部分版本里悄悄取当前数值（默认全 0）走分支，于是carbon_cost恒等于 0。这种失败最坑人——模型能解、不发警告，但改tier、改price、把碳价翻三倍，总成本曲线一模一样。判断方法很简单：把碳价参数全置 0 跑一遍，再恢复跑一遍，如果两次总成本差异小于千分之一，碳成本根本没进目标函数。
还有第三种：分段逻辑里用了value()、double()这类取值函数，等于把符号变量在求解前固化成常数，模型里留下一个常数项，同样静默失效。
2.2 增量分段法：价格单调递增时最省事
真正能进 MILP 的写法是引入辅助变量，把“总量落在哪一段”交给约束表达，而不是交给 if 判断。假设碳排总量按吨计，分三段计价，第一段 0～2000 吨单价 100 元/吨，第二段 2000～4000 吨单价 150，第三段 4000～6000 吨单价 200：
% 阶梯碳交易：增量分段线性化 tier = [2000, 4000, 6000]; % 各分段的右端点，单位：吨 price = [100, 150, 200]; % 各分段单价，单位：元/吨，必须单调递增 seg = sdpvar(3,1); % 落在每一段区间内的排放量 E = sum(emission_factor .* P_gen) * dt; % 总排放量，吨 cons = [cons, seg >= 0]; cons = [cons, sum(seg) == E]; % 分段量之和等于总量 cons = [cons, seg(1) <= tier(1)]; % 第一段容量 cons = [cons, seg(2) <= tier(2) - tier(1)]; % 第二段容量 cons = [cons, seg(3) <= tier(3) - tier(2)]; % 第三段容量 carbon_cost = price * seg; % 纯线性表达式，MILP 友好这段代码能成立的关键在于单价递增。目标函数是最小化，如果某个解让seg(2)有值而seg(1)没填满，把它挪回第一段一定更便宜，所以最优解会自动从低价段往高价段依次填满，不需要任何 0-1 变量。这也是为什么增量法比 big-M 法快——它保持了 LP 松弛的紧致性，分支定界不受额外二进制拖累。
代价是它有前提：单价必须严格单调递增。如果政策设计成“用量越大单价打折”，或者分段里出现持平甚至下降，增量法会给出错误的成本，此时必须换 sos2 或二进制。这条限制几乎没有教程会写，但它是实际项目里最容易出错的地方。
参数上有两个坑要盯住。一是量纲：tier用吨、price用元/吨，而很多文献给的是元/kg，直接抄过来会差 1000 倍，模型会把碳成本算成可以忽略的小量，阶梯永远不会被触发。二是E的时间口径：P_gen是 MW，dt是小时，乘出来是 MWh，再乘排放因子（吨/MWh）才是吨。24 个时段各 1 小时，dt=1可以省；换成 15 分钟粒度就必须写上，否则排放量直接差 4 倍。
2.3 价格非单调或需要严格分段：sos2 与二进制的写法
sos2（第二类特殊有序集）的思路是把成本曲线表示成断点的凸组合，YALMIP 里直接有对应函数：
% sos2 写法：先手算各断点的累计碳成本 pts_E = [0, 2000, 4000, 6000]; % 断点排放量，吨 pts_cost = [0, 200000, 500000, 900000]; % 断点累计成本，元 lambda = sdpvar(4,1); cons = [cons, sum(lambda) == 1, lambda >= 0, sos2(lambda)]; cons = [cons, E == pts_E * lambda]; % 用凸组合表示当前排放量 carbon_cost = pts_cost * lambda; % 成本随之插值pts_cost必须手算，这一步很多人算错。0～2000 吨按 100 元/吨是 20 万；2000～4000 吨按 150 是 30 万，累计 50 万；4000～6000 吨按 200 是 40 万，累计 90 万。如果只写各段单价乘长度，第二段就变成 300000 而不是累计值，成本曲线整体错位，碳价灵敏度分析全废。sos2约束保证lambda里最多两个相邻分量非零，从而落在一条线段上，价格单调与否都能正确表达。
二进制写法适合分段很少、又需要额外逻辑（比如“超过阈值才启动碳捕集”）的场景：
% big-M 写法：M 用动态上界，别拍 1e6 M1 = max(1, E_max - tier(1)); % 见第 4 章，E_max 为理论上界 z = binvar(2,1); % z(1)=1 表示已越过第一段 cons = [cons, z(1) <= z(2)]; % 单调性：先越第一段再越第二段 cons = [cons, E - tier(1) <= M1 * z(1)]; % 未越界时 E 不超过 tier(1) cons = [cons, E - tier(1) >= -M1 * (1 - z(1))]; % 越界后 E 可超过 tier(1)这一版只展示第一段跃迁的骨架，第二段同理复制。要提醒的是implies()内部同样展开成 big-M，24 个时段、5 台设备如果有 120 条implies，等价于 120 组 M，数值条件会明显变差。能用增量法就别上二进制。
2.4 求解后必须做的两步校验
模型里怎么算是一回事，算得对不对是另一回事。建议把校验函数和模型彻底分开：模型里只用符号表达式，校验时用纯数值函数复算。
% 纯数值分段函数，只用于事后校验，不进模型 function c = tiered_ref(E, tier, price) edges = [0, tier]; c = 0; for k = 1:numel(price) lo = edges(k); hi = edges(k+1); c = c + price(k) * max(0, min(E, hi) - lo); end end % 求解后对比 E_val = value(E); c_val = value(carbon_cost); c_ref = tiered_ref(E_val, tier, price); fprintf('模型碳成本 %.2f / 复算 %.2f\n', c_val, c_ref); assert(abs(c_val - c_ref) <= 1e-4 * max(1, c_ref), '碳成本与分段函数不一致');第二步是看阶梯触发情况。把每个时段的E打出来，标出它落在第几段，再对照火电出力曲线。正常结果应该出现“碳排接近阈值时火电主动压出力”的拐点；如果 24 个时段全落在第一段，要么碳价给低了，要么量纲错了；如果一上来就顶在第三段，检查tier是不是写成了 kW 口径的数值。
3. CPLEX 返回 Infeasible：储氢递推的索引偏移与效率折算方向
3.1 末时段状态不闭合：递推式里到底丢了哪一段功率
原始版本里储氢状态方程是这样写的：
% 错误写法：功率索引比状态索引慢一拍 for t = 2:24 cons = [cons, H2(t) == H2(t-1) + eta_ele*P_ele(t-1) - P_fc(t-1)/eta_fc]; endH2(t)描述的是 t 时段结束时的储量，它应该由 t 时段自身的充放功率决定，而不是 t−1 时段的。这个偏移带来三个连锁反应。第一，最后一个时段（t=24）的电解槽和燃料电池功率对任何状态变量都不产生影响，模型的“免费时段”出现了：电解槽可以在 24 时段满功率运行却不占用储氢空间，燃料电池也能凭空发电，功率平衡莫名其妙被满足，总成本虚低。第二，第一个时段的功率被算了两次或完全丢失，具体取决于H2(1)是固定值还是变量。第三，储氢曲线末端会出现台阶式跳变，这正是 A 同学看到的“负储量”的来源——末时段被抽干，越过了下限约束，而求解器在容差范围内没报不可行。
正确写法是索引对齐，并显式写出步长：
% 正确写法：状态与功率同索引，dt 为时段步长（小时） H2 = sdpvar(T,1); cons = [cons, H2(1) == H2_init]; % 初值固定，不作为决策变量 for t = 2:T cons = [cons, H2(t) == H2(t-1) ... + eta_ele * P_ele(t) * dt ... % 耗电产氢 - P_fc(t) / eta_fc * dt]; % 耗氢发电 end另一种常见做法是把状态定义成T+1维：H2(1)是初始值，H2(t+1) = H2(t) + ...，闭合约束写成H2(T+1) == H2(1)。两种写法等价，但千万别把两种索引混在一个模型里，那是比偏移更难查的错误——它不会报错，只会让储氢曲线在中间某处对不上。
除了索引，还要加容量上下限和周期闭合：
cons = [cons, H2_min <= H2, H2 <= H2_max]; % 例如 50 ~ 500 kg cons = [cons, H2(T) == H2_init]; % 日周期闭合闭合约束看着可有可无，实际影响很大。不加它，模型会把初始储氢当成免费燃料在当天用光，总成本能虚低 5% 以上；加上之后，储能调度才是一个可循环的方案。这是能源调度模型最经典的“后悔药”约束，漏写一次，后面所有结论都要重跑。
3.2 燃料电池效率写反：氢气凭空多出来的三种症状
电解槽和燃料电池的效率折算方向相反，这一点写错的概率极高。制氢侧：耗电P_ele，产氢按效率打折，m_H2 = eta_ele * P_ele，效率 0.72 意味着产氢能量小于耗电，合理。发电侧：耗氢m_H2_in，发电P_fc = eta_fc * m_H2_in，效率 0.58 意味着耗氢能量大于发电，所以状态方程里必须写成- P_fc / eta_fc，除以效率，而不是乘。
写反的第一种症状是氢量单调上涨：如果写成- eta_fc * P_fc，耗氢被缩小到 58%，模型发现“制氢—储氢—发电”这条链路几乎不亏，于是让电解槽和燃料电池同时满功率运行，储氢量长期顶在上限，火电出力被压得很低。第二种症状是碳成本几乎不触发阶梯，因为火电被氢能替掉了，结论变成“引入电制氢可降碳 30%”，实际上全是效率折算错误导致的假象。第三种是往返效率大于 1：如果两个方向都乘效率，电—氢—电的净效率变成eta_ele / eta_fc ≈ 1.24，模型会靠“充电再放电”白赚钱，这已经违反热力学，但 MIP 求解器不在乎物理，它只在乎约束。
一条硬性自检足够挡住大部分错误：
% 往返效率必须小于 1，否则模型必然套利 assert(eta_ele * eta_fc < 1, '电-氢-电往返效率异常，检查效率折算方向'); fprintf('往返效率 = %.3f\n', eta_ele * eta_fc); % 0.72*0.58 ≈ 0.418求解后再验一遍氢平衡：把每个时段的“储量变化 − 制氢量 + 耗氢量”算成残差，最大值应该在 1e-6 量级。
resid = value(H2(2:end) - H2(1:end-1) ... - eta_ele*P_ele(2:end)*dt + P_fc(2:end)/eta_fc*dt); fprintf('储氢方程最大残差 %.3e\n', max(abs(resid)));把怀疑的约束先写成显式残差变量放进模型里,再在求解后打印value(resid)，比事后猜哪条约束有问题快得多。
3.3 不可行的定位顺序：先分组，再逐时段
CPLEX 报Infeasible时，sol.problem通常返回非 0，sol.info里带一行简短说明。第一件事是把约束分组，逐组松弛，看去掉哪一组后可行：
% 把约束按物理含义分组保存，排查时逐组剥离 groups = {cons_balance, cons_storage, cons_device, cons_carbon, cons_bound}; for k = 1:numel(groups) cons_test = []; for j = 1:numel(groups) if j ~= k cons_test = [cons_test, groups{j}]; % 跳过第 k 组 end end s = optimize(cons_test, [], sdpsettings('solver','cplex','verbose',0)); fprintf('剥离第 %d 组后 problem = %d\n', k, s.problem); end去掉哪组约束后变成可行，冲突就在那组里。实际项目里最常见的三处：储氢上限设得比单时段最大制氢量还小（物理上直接冲突）、某时段负荷大于所有机组出力上限之和、末端闭合约束和初值不匹配。第二种情况在 24 时段模型里很隐蔽，因为负荷峰值往往只出现在一两个时段，把sum(P_max)和max(load)打出来看一眼就能排掉。
第二步是逐时段截断。把约束按t分组，只保留前 6、12、18 个时段分别求解，找到第一个不可行的时段数，问题就锁定在那一个小时附近。这个方法比读求解器日志快，也不依赖任何高级功能。
4. CPLEX 求解超时与 big-M 敏感：把 M 从拍脑袋改成算出来
4.1 M 取值不当的三种典型症状
原始版本里踩过的坑是 M 值反复调不稳。症状分三类。第一类是求解慢：LP 松弛被过大的 M 撑得极松，分支定界下界提升缓慢，24 个时段、几十个二进制变量也能跑半小时以上，日志里 Gap 长时间停在 30% 不动。第二类是数值不稳：系数矩阵里 1e6 和 1 混在一起，CPLEX 在日志里提示数值困难，返回的解违反约束却声称最优，回代之后残差在 1e−3 量级。第三类是偶发不可行：M 给小了，某些物理上可行的解被错误切掉，换一组负荷数据就报 Infeasible，换回来又好。
根因都一样：M 的作用是“大到足以让约束在另一种逻辑状态下失效”，它只需要覆盖该约束在物理上可能超出的最大幅度，而不是覆盖整个变量的全局上界。很多人习惯M = 1e6图省事，等于给求解器递了一个几乎无信息的松弛。
4.2 动态计算 M：用物理上界替代经验值
M 是可以算出来的。以碳排分段为例，系统 24 小时的理论最大排放量由各机组容量和排放因子决定：
% 动态上界：各机组满发时的最大排放量 P_gen_max = [200; 150; 100]; % 机组最大出力，MW ef = [0.85; 0.90; 0.75]; % 排放因子，吨/MWh T = 24; dt = 1; % 时段数与步长 E_max = sum(P_gen_max .* ef) * T * dt; % 理论上界，吨 % 每条跃迁约束只需覆盖“可能超出该断点”的幅度 M1 = max(1, E_max - tier(1)); M2 = max(1, E_max - tier(2));假设算出来E_max = 7000吨，tier(1) = 2000，那M1 = 5000，比 1e6 小两个数量级。如果某条约束只涉及储氢罐，M 就取罐容上限减阈值，量级更小。一个可操作的规则：M 的量级不应超过该约束里最大物理量的两倍，超了就是没算。
碳成本用增量法时其实不需要 M，这也是为什么第 2 章推荐增量法——少一组二进制、少一组 M，求解时间和数值稳定性都会好一截。只有在需要额外逻辑（启停、阈值触发、分段互斥）时才不得不用 big-M。
4.3 求解器参数：该动的和不该动的
参数不是越多越好，下面这几项是实际调过的：
ops = sdpsettings('solver','cplex','verbose',2); ops.cplex.timelimit = 1800; % 求解时限，秒 ops.cplex.mip.tolerances.mipgap = 0.01; % 相对间隙 1%，工程够用 ops.cplex.tolerances.feasibility = 1e-7; % 可行性容差，别放大到 1e-4 ops.cplex.exportmodel = 'model.lp'; % 导出模型，便于离线排查| 参数 | 作用 | 建议取值 | 说明 |
|---|---|---|---|
timelimit | 硬性时限 | 600～3600 | 调度类问题 24 时段通常分钟级，超时说明建模有问题 |
mip.tolerances.mipgap | 相对间隙 | 0.01 | 精确到小数点后三位对调度无意义，1% 足够 |
feasibility | 可行性容差 | 1e-7 | 放到 1e-4 会让“负储量”这类错误藏进容差里 |
mip.strategy.bbinterval | 分支策略 | 1 | 加强分支，变量规模小时收益明显 |
verbose | 日志级别 | 2 | 调试期开 2，正式跑改 0 |
verbose=2的价值在于能看到根节点松弛、每次 incumbent 更新和 Gap 变化。如果根节点 Gap 就超过 50%，问题多半在 big-M 或约束线性化，不在参数；如果 Gap 卡在 5% 十几分钟不动，才考虑调mipgap或给初始解。exportmodel导出的 lp 文件可以离线用求解器的冲突检测功能定位不可行约束，比在 MATLAB 里反复加日志省事。
单位统一是另一件容易被忽略的事。功率统一用 MW、能量用 MWh、碳排用吨、成本用元，整张系数矩阵的量级就压在 1e0～1e3。混用 kW 和 MW 会让矩阵出现 1e3 的量级差，求解器日志里的数值警告基本都来自这里，而不是算法本身。
5. 避坑清单：阶梯碳交易与电制氢建模里最容易踩的五个坑
5.1 建模层的三个坑
坑一：碳成本函数用 if-else 或 value() 实现。现象是模型能解、没有警告，但改碳价参数结果不变，灵敏度分析画出一条水平线。原因是对sdpvar做逻辑判断时，MATLAB 要么报错、要么取默认数值走固定分支，碳成本变成常数项。解决办法是把分段逻辑改写成增量法或 sos2 的约束形式，纯数值的分段函数只保留在求解后的校验环节，两套代码物理隔离。
坑二：储氢递推索引比功率慢一拍。现象是末时段储氢量跳变、出现负储量，或者储氢曲线整体滞后一个小时。原因是H2(t)由P_ele(t-1)决定，末时段的功率不进状态、首时段功率被重复计算或丢失。解决办法是把状态方程统一成H2(t)对应P_ele(t)、P_fc(t)，并显式乘上步长dt；采用T+1维写法时，闭合约束写成H2(T+1) == H2(1)，两种索引不要混用。
坑三：燃料电池效率乘反，或末端储量不闭合。现象是氢气用量长期顶在上限、火电出力异常低、碳成本几乎不触发阶梯，结论变成“电制氢大幅降碳”。原因是耗氢量被写成eta_fc * P_fc，等于把 0.58 的效率当成增益；或者模型把初始储氢当免费燃料用光，成本虚低。解决办法是发电侧统一用P_fc / eta_fc，加eta_ele * eta_fc < 1的断言，并补上日周期闭合约束。
5.2 求解层的两个坑
坑四：弃风与制氢没有硬性耦合，出现“一边弃风一边制氢”的解。现象是最优解里P_curtail > 0的同时电解槽还在满功率吃电，看起来自相矛盾。原因是功率平衡允许多种来源，而弃风惩罚系数小于制氢带来的碳收益，模型自然选择“多发电、多弃、多制氢”这条更划算的路径。解决办法是给弃风惩罚系数设一个下限（大于单位弃风的机会成本），并补一条P_ele(t) <= P_wind_avail(t) - P_curtail(t) + P_gt(t) + P_fc(t) - P_load(t)的耦合约束，把逻辑堵死。
坑五：big-M 拍 1e6，或者每条 implies 各带一个 M。现象是求解时间忽长忽短、日志出现数值警告、换一组数据就偶发不可行。原因是 M 远大于物理可能幅度，LP 松弛失效，同时系数矩阵量级差被放大。解决办法是按第 4 章的方式用机组容量和排放因子算出动态 M，能用增量法就不用二进制，implies的数量控制在个位数。
6. 快速定位：拿到一份跑不通的调度代码，我先看哪三处
第一步永远先看sol.problem和sol.info，别急着读模型。problem非 0 时先分清类型：不可行、无界、数值问题，三者的排查路径完全不同。不可行走第 3 章的分组松弛，无界通常是某个变量缺了上界或者成本系数符号写反，数值问题先查单位量纲和 big-M。花三十秒确认类型，比盲改两小时约束有效。
| 求解返回 | 大概率原因 | 第一个动作 |
|---|---|---|
Infeasible | 储氢上下限冲突、末端闭合与初值矛盾、负荷超总容量 | 约束分组剥离，逐组求解 |
Unbounded | 变量缺上界、成本系数符号错、/eta写成*eta | 检查所有sdpvar是否有界 |
| 数值警告 | 单位混用、big-M 过大 | 统一 MW/MWh/吨，动态算 M |
| 能解但结果怪 | 碳成本被常数化、效率折算反、索引偏移 | 回代残差校验 |
第二步只跑可行性，不带目标函数：optimize(cons, [], ops)。这一步能过，说明约束体系自洽，问题在目标函数或者线性化方式；过不了，就老老实实做分组剥离。很多人跳过这步直接调求解器参数，纯属浪费时间。
第三步是把解代回约束算残差，这一步能抓住所有“能解但错”的问题：
% 关键约束逐条回代，残差超过阈值就报警 tol = 1e-6; % 1) 功率平衡 r_bal = value(P_gt + P_wind_used + P_fc - P_load - P_ele - P_curtail); fprintf('功率平衡最大残差 %.3e\n', max(abs(r_bal))); % 2) 储氢动态 r_h2 = value(H2(2:end) - H2(1:end-1) ... - eta_ele*P_ele(2:end)*dt + P_fc(2:end)/eta_fc*dt); fprintf('储氢方程最大残差 %.3e\n', max(abs(r_h2))); % 3) 阶梯碳成本：与纯数值函数复算对比 E_val = value(E); c_val = value(carbon_cost); c_ref = tiered_ref(E_val, tier, price); fprintf('碳成本 模型 %.2f / 复算 %.2f\n', c_val, c_ref); assert(max(abs(r_bal)) < tol && abs(c_val - c_ref) < 1e-3*max(1,c_ref), ... '结果未通过回代校验，别急着画图');三条残差各管一类故障：功率平衡残差超限说明约束漏写或索引错位；储氢残差超限说明效率折算方向或步长有问题；碳成本对不上说明分段线性化写错或者量纲差 1000 倍。这套校验我一般封装成一个check_model.m，每次改完模型先跑它再画图，能省掉大量“图看起来还行、结论全错”的返工。
还有一个习惯值得养成：改完任何一处参数，先把碳价全置零跑一遍作为基线，再恢复参数跑对照。两次结果的差异如果小于千分之一，说明碳交易模块根本没生效，后面的成本对比、降碳比例都没有意义。这个基线对照花不了两分钟，但挡住过我好几次“参数改了、模型没读到”的乌龙。
从那以后我每次接手这类带分段价格和储能耦合的调度模型，都会强制走一遍“分组可行性 → 残差回代 → 参数基线对照”的流程，哪怕代码是自己上周刚写的。故障从来不是因为模型复杂，而是因为某个索引偏移了 1、某个效率乘除反了，而这些东西不会自己冒出来报警。希望帮到你。
本文还有配套的精品资源，点击获取
