MATLAB/YALMIP综合能源调度：阶梯碳交易与电制氢建模避坑

简介：面向能源系统优化方向的研究人员与具备MATLAB基础的工程技术人员，这份PDF文档围绕阶梯型碳交易机制与电制氢模块的综合能源系统优化调度展开，重点解决碳成本跃迁与电—氢—电多向转换下的经济环保协同调度问题。包内仅含1个PDF文件，约129KB，篇幅紧凑，便于直接查阅模型思路、约束构建与算例结果。文档结合YALMIP建模与CPLEX求解，拆解传统机组、风电、电解槽、储氢罐和燃料电池五单元架构，讲清阶梯碳交易分段线性化、储氢罐状态方程、目标函数设计及可视化输出，并给出完整代码片段与关键注释。内容还提示big-M取值、时间索引对齐等数值稳定性与建模细节，适合用于论文复现、课程设计或项目方案比选。已有163人学习，可作为理解综合能源系统低碳调度建模的参考材料。

1. 从一条 Nonlinear operator 报错说起：阶梯碳价这样写必翻车

用 MATLAB 加 YALMIP 搭综合能源系统优化调度，模型里只要同时放进阶梯碳交易和电制氢，头一天大概率会撞上一条不太起眼的报错：Nonlinear operator ... not supported，或者更隐蔽的Warning: Solver not applicable (cplex)。模型建得挺顺，变量、约束一条没漏，optimize一跑，CPLEX 直接退场。阶梯碳交易的分段计价、电解槽效率曲线、燃料电池的双向折算，这三处都能把本该是 MILP 的问题悄悄变成非线性，等你发现时已经调了两小时参数。这篇写给正在跑 24 时段调度、模型能建起来却解不出结果的人。

2. 阶梯碳交易报错：min/max 一上手，MILP 就变成了非线性

2.1 报错现场：min/max 一出现，CPLEX 就退场

最典型的一行代码长这样，本意是“200 吨以内按 80 元、超出部分按 120 元”，写法上很像普通 MATLAB 函数，语义也对，但在 YALMIP 世界里它已经不是线性表达式了：

% 反例：看起来人畜无害，实际把模型拖成了非线性 E_total = sum(emission_coef .* P_gen) * dt; % 总排放量，sdpvar 线性表达式 carbon_cost = 80*min(E_total, 200) + 120*max(0, E_total - 200); cons = [cons, P_gen >= 0, P_gen <= P_max]; sol = optimize(cons, fuel_cost + carbon_cost, sdpsettings('solver','cplex'));

跑下去大概率看到Warning: Solver not applicable (cplex)，或者更直白的Nonlinear operator (max) not supported。原因不复杂：min和max在 YALMIP 里对sdpvar属于非光滑算子，它不会替你展开成分段线性，只会把这个表达式标记成“非线性”。CPLEX 是 MILP/MIQP 求解器，面对非线性项只有两条路，要么换成 IPOPT 之类的非线性求解器（这类问题它是凸的吗？不一定），要么求解器匹配失败直接退出。所以我一般看到这个 warning，第一反应不是去换求解器，而是回去看目标函数里有没有偷懒写的min、max、abs。

顺带提一句，if-else版本的阶梯函数更隐蔽：

% 反例二：if-else 只对数值生效，对 sdpvar 完全无效 function cost = carbon_cost(emission) if emission <= 200 cost = 80 * emission; else cost = 80*200 + 120*(emission-200); end end

sdpvar传进去，emission <= 200返回的是一个约束对象而不是true/false，MATLAB 的判断分支会走进一个你完全没预期的方向，成本项可能直接变成常数 0 或者报维度错误。这类 bug 最坑的地方在于它不报错，只是让碳成本悄悄消失，最后结果“看起来很合理”，成本还特别低。

2.2 分段变量法：把阶梯写成一组线性约束

正确做法是把“排放量落在哪一档”这件事交给辅助变量和 0-1 变量表达，保持整个模型是 MILP。三档写法如下，档位数和阈值都可以改：

% 阶梯碳交易：三档分段，用非负辅助变量拆分排放量 tier = [200, 500, 1e4]; % 各档区间上限（吨），最后一档给足够大的上界 price = [80, 120, 160]; % 各档碳价（元/吨），必须递增 nTier = numel(price); E_seg = sdpvar(nTier, 1); % 各档对应的排放量 z = binvar(nTier, 1); % 档位开关，1 表示该档已被占用 cons = [cons, E_seg >= 0]; cons = [cons, sum(E_seg) == E_total]; % 分段之和等于总排放 cap = [tier(1), tier(2)-tier(1), tier(3)-tier(2)]; % 每档的容量 for k = 1:nTier cons = [cons, E_seg(k) <= cap(k) * z(k)]; % 大 M 用档位容量本身 end % 档位顺序：低档没占满，高档不许开 for k = 2:nTier cons = [cons, z(k) <= z(k-1)]; end carbon_cost = price * E_seg; % 线性加权，CPLEX 认得

逻辑上一共三件事。第一，E_seg把总排放量切成互不重叠的三段，sum(E_seg) == E_total保证不丢量也不多算。第二，E_seg(k) <= cap(k)*z(k)里的 M 直接用该档容量，而不是随手写 1e6，这一点决定了后面数值稳不稳。第三，z(k) <= z(k-1)强制阶梯按顺序爬升。

关于第三条要多说一句：在价格递增、目标是极小化的前提下，求解器本来就会优先把排放量塞进便宜档，这条顺序约束看起来多余。但实际跑下来，缺少它时解的退化很常见——某些时段z的组合不唯一，E_seg在两个等价解之间跳，导致你对比两天结果时发现“同样的输入成本差了 3 块钱”，排查半天以为是随机性。加上顺序约束，解唯一，曲线才稳定。这也是我踩过一次之后固定加进去的习惯。

参数上还有一点：tier最后一档给1e4这种“足够大但不夸张”的值，不要写inf，YALMIP 对inf系数的处理会让 CPLEX 直接报数值错误；也不要写1e8，量级差太远会让矩阵条件数变得很难看。

2.3 阈值和单价怎么定：单位统一比数值大小重要

比写对分段更容易翻车的是单位。我见过同一份代码里，碳价一处写80（元/吨），另一处写1.2（元/kg），两处都“看起来合理”，合起来算出的碳成本差了三个数量级，模型当然还是能解，只是碳交易完全不起作用——因为碳成本相对燃料成本可以忽略，优化器压根不理会它。

最省事的验证办法是打印量纲：disp(value(E_total))，如果 24 小时总排放跑出来是 0.05 这种数，说明某处少乘了功率量级；如果跑出来是 5e5，多半是重复乘了时段数。碳成本占燃料成本的比例落在 5%~25% 之间通常比较正常，低于 1% 基本可以怀疑单位没对齐。

2.4 用固定碳价做一次对照验证

阶梯模型写完后，别急着看漂亮的调度曲线，先跑一次固定碳价对照：

% 对照组：把阶梯拍平成单一碳价，验证模型主体没错 carbon_cost_flat = 100 * E_total; sol_flat = optimize(cons, fuel_cost + carbon_cost_flat, ops); % 实验组：阶梯碳价 sol_tier = optimize(cons, fuel_cost + carbon_cost, ops); fprintf('固定碳价总成本: %.2f\n', value(fuel_cost + carbon_cost_flat)); fprintf('阶梯碳价总成本: %.2f\n', value(fuel_cost + carbon_cost));

正常情况下阶梯碳价的总成本不会低于固定碳价（同等排放下阶梯只会更贵或持平），但机组出力会往低碳方向偏，排放量应该更低。如果跑出来阶梯版的成本反而更低、排放还更高，那基本可以确定分段变量写反了或者price没递增。这个对照跑一次只要几十秒，比后面盯着诡异曲线猜半天划算得多。

3. 电制氢模型翻车：储氢量解出负值、效率折算方向反了

3.1 负储氢量从哪来：效率折算方向与索引错位

原文那行状态方程我照着敲过一遍：

% 有问题的写法：t=1 被跳过，t=24 状态不闭合 for t = 2:24 cons = [cons, H2_storage(t) == H2_storage(t-1) + ... units.ele_eff*P_electro(t-1) - P_fuelcell(t-1)/units.fuel_cell_eff]; end

两个毛病叠在一起。一是t=1的储氢量完全没被约束，求解器可以随便给它一个负值，只要t=2之后能圆回来；二是右边用的是t-1时刻的功率，意味着第 24 个时段的电解槽和燃料电池功率压根没进状态方程，最后一个时段的状态是“悬空”的，H2_storage(24)只受前面影响，末端想放空就放空。跑出来的曲线在最后几个点突然砸到负值，看起来像“负储能”，实际是约束没闭合。

效率折算方向是另一个高频错误。电解槽是电转氢，产氢量 = 效率 × 用电量；燃料电池是氢转电，耗氢量 = 发电量 / 效率。一个乘一个除，写反了不会报错，只会让整个氢能环节的能量凭空多出来或者凭空少掉，最直观的表现就是“加了氢能之后总成本降了 20%”，看着像好事，实际是能量不守恒。

% 正确写法：统一用 t 时刻功率，首末两端都闭合 T = 24; dt = 1; H2 = sdpvar(T,1); P_el = sdpvar(T,1); % 电解槽用电功率 MW P_fc = sdpvar(T,1); % 燃料电池发电功率 MW eta_el = 0.72; % 电解槽：电 -> 氢 eta_fc = 0.58; % 燃料电池：氢 -> 电 H2_init = 200; % 初始储氢量 kg H2_min = 20; H2_max = 900; cons = [cons, H2(1) == H2_init + eta_el*P_el(1)*dt - P_fc(1)*dt/eta_fc]; for t = 2:T cons = [cons, H2(t) == H2(t-1) + eta_el*P_el(t)*dt - P_fc(t)*dt/eta_fc]; end cons = [cons, H2 >= H2_min, H2 <= H2_max]; % 末端闭合：调度周期结束时储氢量回到初值附近 cons = [cons, H2(T) >= 0.95*H2_init, H2(T) <= 1.05*H2_init];

几点说明。eta_el*P_el(t)*dt得到的才是这一小时实际产氢量，漏掉dt在dt=1时看不出问题，换成 15 分钟粒度立刻错四倍。H2_min不要设成 0，储氢罐留一点安全余量，也顺便避免求解器把罐子榨干到恰好 0 这种物理上不可行的边界解。末端闭合约束是最容易被忽略的一条，不加的话优化器一定会在第 24 小时把氢气全放掉换成电，因为“反正后面没有时段了”，这在数学上最优，在工程上纯属耍赖。

3.2 启停约束与 0.3 倍运行下限

电解槽不能想开就开、想停就停，低负载运行时效率掉得厉害，所以通常给一个最小出力比例。原文里那个 0.3 的系数就是干这个的，但要配合 0-1 变量才能表达“要么不运行，要么至少 30% 额定功率”：

u_el = binvar(T,1); % 电解槽运行状态 P_el_max = 800; % 额定功率 MW P_el_min = 0.3 * P_el_max; % 最小技术出力 for t = 1:T cons = [cons, P_el(t) >= P_el_min * u_el(t)]; cons = [cons, P_el(t) <= P_el_max * u_el(t)]; end

这两条不等式组合起来就是一个大 M 形式的开关：u_el(t)=0时功率被压到 0，u_el(t)=1时功率被夹在[240, 800]之间。中间那段 0 到 240 的区间被主动排除掉了，这正是“最小技术出力”在 MILP 里的标准表达。

启停次数限制要引入上升沿变量，别直接对u_el求和：

v_el = binvar(T,1); % 启动动作标识 for t = 2:T cons = [cons, v_el(t) >= u_el(t) - u_el(t-1)]; cons = [cons, v_el(t) <= u_el(t)]; cons = [cons, v_el(t) <= 1 - u_el(t-1)]; end cons = [cons, sum(v_el) <= 3]; % 24 小时内最多启动 3 次

v_el(t) >= u_el(t) - u_el(t-1)只能保证“启动时v_el不小于 1”，但对“不启动时v_el能不能等于 1”没约束，目标函数里如果没有和v_el相关的成本项，求解器会随便给值，导致sum(v_el) <= 3这条约束要么过紧要么白给。补上后两条上界约束，v_el才真正等于上升沿。如果模型里本来就加了启动成本系数，前一条就够用，但把三条都写上更保险，代价只是多两行约束。

3.3 电解槽效率非线性：分段线性化别用双线性

电解槽效率随负载率变化，通常是先升后降的一条曲线，满载和低载都不是最优点。直接把eta(P)*P塞进模型就是非线性，只能用分段线性逼近。常见做法是用断点加权重：

% 电解槽效率曲线分段线性化：4 个断点切成 3 段 load_bp = [0.30, 0.50, 0.80, 1.00]; % 负载率断点（标幺） eta_bp = [0.62, 0.72, 0.70, 0.64]; % 对应效率，按实测曲线填 K = numel(load_bp) - 1; lam = sdpvar(K+1, T); % 断点权重 zseg = binvar(K, T); % 段选择 for t = 1:T cons = [cons, lam(:,t) >= 0, sum(lam(:,t)) == 1]; cons = [cons, sum(zseg(:,t)) == 1]; cons = [cons, lam(1,t) <= zseg(1,t)]; for k = 2:K cons = [cons, lam(k,t) <= zseg(k-1,t) + zseg(k,t)]; end cons = [cons, lam(K+1,t) <= zseg(K,t)]; % 功率与产氢量都用断点的凸组合表达 cons = [cons, P_el(t) == P_el_max * (load_bp * lam(:,t))]; cons = [cons, H2_prod(t) == P_el_max * ((load_bp .* eta_bp) * lam(:,t))]; end

关键在最后两行：P_el和H2_prod都是断点值的凸组合，两者共用一个lam，所以效率关系被自动保持在这条折线上，而且全程线性。zseg保证同一时刻只有相邻两个断点权重非零，否则lam可以在三个断点上同时取正，相当于把效率曲线“拉平”，得到一条并不存在的效率。这个约束一漏，模型照样能解，只是效率永远取到折线凸包上界，产氢量被系统性高估。

有个细节要避开：不能写成P_el(t) == P_el_max * (load_bp*lam(:,t)) * u_el(t)。u_el(t)是 0-1 变量，乘上去就是双线性项，CPLEX 立刻就懵了。实际上不需要乘——前面P_el(t) <= P_el_max*u_el(t)已经保证u_el=0时功率为 0，lam自然也只能落在零点附近。我第一次写的时候顺手乘了一下，换来的是一条Nonlinear operator报错和半小时的自我怀疑。

4. 求解器罢工：CPLEX 报 Infeasible 与 big-M 数值敏感

4.1 Infeasible 的八成原因在 big-M

模型能建、能提交，跑出来CPLEX Error 1217: No solution exists，或者 YALMIP 层面给你一句Infeasible problem。这种情况我一般先怀疑三件事：约束互相矛盾、储氢末端闭合给死了、big-M 把可行域切掉了。

第三点最隐蔽。big-M 取太小，物理上可行的运行方式被数学上排除；取太大，系数矩阵量级跨度拉开，CPLEX 的分支定界在数值容差附近反复横跳，最后报“不可行”或者给出一个明显违反约束的解。

% 动态计算 M：按每个约束里变量的真实量程取，不用拍脑袋的 1e6 M_el = P_el_max; % 电解槽功率上界 M_h2 = H2_max - H2_min; % 储氢量可变化范围 M_curt = wind_cap * ones(T,1); % 弃风量按时段风电出力上限取 % 弃风量约束：弃风 <= 该时段可用风电 for t = 1:T cons = [cons, P_curt(t) <= P_wind_avail(t)]; cons = [cons, P_curt(t) <= M_curt(t) * z_curt(t)]; % 只有允许弃风时才放开 end

原则就一条：M 取“这个约束里变量真正能取到的最大差值”，不要全局统一一个 1e6。上面M_curt按每个时段的风电可用出力单独给，比全局 150 这种常数更紧，求解器搜起来轻松很多。判断标准也简单：如果你的 M 和变量典型值的量级差超过 1e4 倍，基本可以认为 M 取松了。

4.2 用调试模式和松弛变量定位不可行约束

硬猜没意义，让求解器自己说。

% 第一步：拿到求解器的原始返回，而不是只有一句 Infeasible ops = sdpsettings('solver','cplex','verbose',2,'debug',1); sol = optimize(cons, total_cost, ops); if sol.problem ~= 0 disp(sol.info); % CPLEX 原始输出 yalmiperror(sol.problem) % 转成人话 end

sol.problem的取值要记几个：0成功，1不可行，2无界，3超时但有可行解，4数值问题。看到1就往下走松弛法，看到3说明模型没写错，只是规模或 M 的问题。

% 第二步：给关键约束加非负松弛，看哪条被撑破 s_bal = sdpvar(T,1); % 功率平衡松弛 s_h2 = sdpvar(T,1); % 储氢上下限松弛 cons_soft = [cons, s_bal >= 0, s_h2 >= 0]; cons_soft = [cons_soft, P_gt + P_wind - P_curt == P_load + P_el - P_fc + s_bal]; cons_soft = [cons_soft, H2 + s_h2 >= H2_min, H2 - s_h2 <= H2_max]; sol_soft = optimize(cons_soft, total_cost + 1e4*sum(s_bal + s_h2), ops); idx_bal = find(value(s_bal) > 1e-6); % 这些时段功率平衡不上 idx_h2 = find(value(s_h2) > 1e-6); % 这些时段储氢越界

惩罚系数1e4的取法有讲究：要比正常成本项大一个量级以上，才能保证求解器优先消除松弛；但也不能到 1e9，否则又是数值问题。跑出来idx_bal非空，说明那段时间的功率根本配不平，回头查负荷数据或风电出力上限；idx_h2非空，说明储氢容量或初始值给得不合理，比如初值 200 却把H2_max设成 150。

4.3 求解参数：mipgap、timelimit 和分支策略

MILP 规模稍大一点，默认参数往往跑不完，或者跑完了但 gap 很大。这几个参数是我固定会调的：

ops = sdpsettings('solver','cplex','verbose',2); ops.cplex.timelimit = 600; % 单次求解上限，秒 ops.cplex.mip.tolerances.mipgap = 0.01; % 相对 gap 1% ops.cplex.mip.tolerances.absmipgap = 1e-4; % 绝对 gap ops.cplex.mip.strategy.bbinterval = 1; % 分支定界更积极 ops.cplex.mip.limits.nodes = 50000; % 节点上限，防内存爆

mipgap设 0.01 对调度类问题完全够用，1% 的成本误差落在工程可接受范围里，而且能把求解时间压下来一大截。timelimit一定要设，不然遇到写错的模型能把机器挂一整晚——这种血泪经验一次就够。设了之后要习惯性检查返回码：

if sol.problem == 3 warning('达到时间上限，当前解 gap = %.4f', sol.info); end

超时返回的解不是不能用，但要明确知道自己拿到的是次优解，别拿着它去论证“阶梯碳价能省 8.6%”这种结论。真要下结论，把mipgap收紧到 0.001 再跑一遍确认趋势没变。

5. 常见问题排查：五个把我卡了半天的报错现场

5.1 现象：储氢罐曲线末端突然掉到负值

现象是最后两三个时段储氢量直接穿到 0 以下，前面的曲线都正常，唯独收尾崩了。原因是状态方程从t=2开始写，第 24 时段的电解槽、燃料电池功率没有进方程，末端状态完全自由。解决办法是把首端H2(1)显式写出来，末端加闭合约束H2(T)落在初值 ±5% 区间内。加完再跑，末端那条断崖就消失了。

5.2 现象：报 Nonlinear operator，但目标函数里看不到 min/max

现象是明明只写了线性成本项，仍然报非线性算子错误。原因是某个自定义函数里藏了abs()或者if-else，比如启停成本写成abs(u(t)-u(t-1))*cost。abs对sdpvar同样是非线性。解决办法是把abs拆成两个方向的辅助变量，或者直接用 5.1 那套上升沿、下降沿变量v_up、v_dn分别计价，全程保持线性。

5.3 现象：碳成本算出来是个常数，跟排放量完全无关

现象是改排放系数、改机组出力，碳成本一分钱不变。原因是阶梯函数里的判断条件对sdpvar失效，走进了某一个固定分支，或者分段变量E_seg的求和约束写反了方向。解决办法是先disp(value(E_total))确认排放量确实在变，再检查sum(E_seg) == E_total这条等式有没有被后面的cons覆盖掉。我遇到过一次是把cons重新赋值而不是拼接，前面的约束全丢了。

5.4 现象：启停次数约束形同虚设，曲线疯狂开关

现象是电解槽在相邻时段反复启停十几次，sum(v_el) <= 3完全没起作用。原因是v_el只写了下界v_el(t) >= u_el(t) - u_el(t-1)，没有上界约束，求解器可以随意把它压到 0 来满足总次数限制。解决办法是补上v_el(t) <= u_el(t)和v_el(t) <= 1 - u_el(t-1)两条上界，v_el才真正等于上升沿。

5.5 现象：同样输入跑两次，结果差一点点

现象是数据没改，重跑一次成本差几块钱，曲线在几个时段有细微抖动。原因是模型存在多解（比如并列最优的机组组合），加上mipgap放到 0.01，求解器在不同分支上停下来的点不一样。解决办法是给阶梯碳价的档位加顺序约束消除退化、把mipgap收到 0.001、必要时给目标函数加一个极小的二次正则项1e-4*sum(P_gt.^2)来挑选唯一解。抖动本身不影响工程判断，但会让你在对比两版模型时误以为改出了效果。

6. 快速定位：拿到一条报错后我的三步走顺序

报错信息通常只有一行，能提供的线索有限，所以顺序比技巧重要。我现在固定按这三步走，基本十分钟内能判断问题出在建模层还是求解层。

第一步看sol.problem返回码，别急着看sol.info那一大段英文。1不可行、2无界、3超时、4数值问题，这四个码对应四条完全不同的排查路径。判断成不可行就直接上松弛变量，判断成超时就去查规模和大 M，别在错误的路径上浪费时间。

第二步看量纲。把value()打出来的关键变量逐个过一遍：总排放量、总成本、储氢量峰值、电解槽利用小时数。这四个数只要有一个明显不合常理（比如 24 小时总排放 0.05 吨，或者储氢罐峰值超过容量），问题一定在建模层，跟求解器没关系。这一步能拦掉我遇到的大半问题。

第三步才是看解的形状。把各单元出力画出来，重点看三个位置：风电大发时段的弃风量和电解槽功率是否对得上、碳排放在阈值附近的时段火电出力有没有拐点、储氢曲线的首末两端是否闭合。形状对了，再回头核对目标函数的各项占比，碳成本占比落在 5%~25% 之间我才会相信这个结果能拿去汇报。

从那以后我每次搭这类模型，都会在optimize之前强制走一遍“打印量纲 + 固定碳价对照 + 松弛变量备用”这三件事，宁可多花五分钟，也不想再经历一次跑到凌晨三点发现碳成本根本没进目标函数。希望帮到你。

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