简介：面向能源系统优化方向的研究人员与具备 MATLAB 基础的开发者，这份 PDF 文档围绕阶梯型碳交易与电制氢的综合能源系统优化调度展开，借助 YALMIP 调用 CPLEX 求解混合整数线性规划模型，处理低碳约束下机组出力与氢能多向转换的协同问题。内容覆盖传统机组、风电、电解槽、储氢罐与燃料电池五个单元的建模，阶梯碳价的分段线性化（big-M 与辅助变量写法）、储氢罐时间耦合状态方程，以及兼顾发电成本、弃风惩罚与碳成本的目标函数设计，并附结果可视化图示便于对照复现。压缩包内为 1 个 PDF 文件，约 129KB，篇幅紧凑。已有 163 人学习。读者可获得可直接迁移的建模框架、分段计价与时间索引对齐等易错点的处理思路，以及风电大发时段优先制氢、碳排放接近阈值时火电出力出现拐点等现象的分析视角。

1. 24小时调度里最贵的不是燃料：阶梯碳价和电制氢把模型推成了MILP

凌晨两点风电顶到装机上限，弃风惩罚和电解槽启停成本同时压在目标函数上，此时决定机组出力顺序的往往不是燃料价格，而是碳排放总量落在第几档。用MATLAB配YALMIP搭一套含阶梯碳交易与电制氢的综合能源系统优化调度，难点从来不是把约束一行行敲出来，而是让这套模型被CPLEX识别成一个结构良好、松弛够紧的混合整数线性规划——同样一份物理模型，分段写法差一行，求解时间就可能从十几秒变成半小时。这份代码包把工程拆成参数初始化、模型构建、结果可视化三块，覆盖传统机组、风电、电解槽、储氢罐、燃料电池五个单元，碳价按三段超额累进，氢侧带储罐动态和日启停次数限制。适合已经会写sdpvar和optimize、但一碰到分段函数和双向效率折算就卡住的人。

2. 系统骨架与参数标定：五类单元、24个时段和那组最容易拍脑袋的系数

2.1 决策变量怎么摆：矩阵还是向量，先定下来再动手

综合能源系统调度的变量天然是二维的——设备 × 时段。新手最常见的翻车点是把机组出力和电解槽功率混在一个矩阵里，后面写约束时又得反复切片，代码可读性直接塌掉。我一般按「同类设备打包成矩阵、异构设备各自成向量」来摆：火电有三台机组，就开一个 24×3 的矩阵；风电、电解槽、燃料电池、储氢量各开一个 24×1 的列向量。这样写循环约束时下标语义统一，调试时value(P_gt(:,2))一眼就能看出是二号机组全天出力。

%% 决策变量定义（T = 24，dt = 1 h） T = 24; Ngt = 3; P_gt = sdpvar(T, Ngt, 'full'); % 火电机组出力 MW P_wind = sdpvar(T, 1, 'full'); % 风电实际消纳 MW P_curt = sdpvar(T, 1, 'full'); % 弃风功率 MW P_ele = sdpvar(T, 1, 'full'); % 电解槽输入电功率 MW P_fc = sdpvar(T, 1, 'full'); % 燃料电池发电功率 MW H_tank = sdpvar(T, 1, 'full'); % 时段末储氢量 kg u_ele = binvar(T, 1, 'full'); % 电解槽运行状态 0/1

逻辑说明：'full'这个参数容易被忽略。YALMIP 默认会对变量做稀疏处理，如果模型里出现大量零系数约束，稀疏结构反而会让 CPLEX 的预求解多做一轮无用功；显式声明 full 可以避免后面P_gt(t,:)这类整行切片出错。

参数说明：P_curt单独设变量而不是用wind_avail - P_wind隐式表达，是为了让目标函数里的弃风惩罚项线性可加，也方便出图时直接对比可用风功率和实际消纳。u_ele必须是binvar，写成intvar再配0<=u<=1虽然数学上等价，但 CPLEX 会走整数分支而不是二进制分支，速度差一截。

2.2 参数标定与量纲对齐：吨、千克、兆瓦时别混着用

碳交易段、储氢罐、燃料成本三套单位体系撞在一起，是这类模型最隐蔽的错误来源。碳排放按吨算，储氢按千克算，功率按兆瓦算，时段长度 1 小时折算成兆瓦时。任何一处没对齐，模型照样能解出来，只是结果离谱到没法解释。我习惯在sys_init.m里把所有换算系数显式命名，禁止在约束里写裸数字。

提示：碳排放因子的单位是「吨每兆瓦时」，而碳价常被写成「元每吨」。如果你手上的碳价数据是「元每千克」，别忘了乘 1000，否则碳成本项会小三个数量级，最后看到火电全天满发还以为是求解器的问题。

2.3 功率平衡与时间耦合：索引对齐是第一个必过的坎

24 时段模型的索引约定必须全局统一：t代表第 t 个时段「结束时刻」的状态，功率变量代表第 t 个时段「内」的平均功率。这两者混用，储氢罐方程就会在首尾两个时段错位，表现为最后时段库存对不上初始值。

%% 功率平衡（逐时段） cons = []; P_load = load_profile; % 24×1，MW for t = 1:T cons = [cons, ... sum(P_gt(t,:)) + P_wind(t) + P_fc(t) == ... P_load(t) + P_ele(t)]; cons = [cons, P_wind(t) + P_curt(t) == wind_avail(t)]; cons = [cons, 0 <= P_curt(t) <= wind_avail(t)]; end %% 机组出力上下限 for t = 1:T cons = [cons, 0 <= P_gt(t,:) <= P_gt_max]; end %% 储氢罐状态方程：H(t) 为第 t 时段末库存 for t = 2:T cons = [cons, ... H_tank(t) == H_tank(t-1) ... + k_ele * P_ele(t) * 1.0 ... % 电解制氢 kg - P_fc(t) / (eta_fc * 33.3) * 1000]; % 耗氢 kg，33.3 kWh/kg 为低位热值 end cons = [cons, H_tank(1) == H0 + k_ele*P_ele(1) - P_fc(1)/(eta_fc*33.3)*1000]; cons = [cons, H_tank(T) == H0]; % 日周期闭合 cons = [cons, 0.1*H_cap <= H_tank <= 0.9*H_cap];

逻辑说明：储氢方程最容易写错的是首时段。如果循环从t = 2开始，H_tank(1)就没有任何方程约束，YALMIP 会把它当成自由变量，求解器随便给个值，最后H_tank(T) == H0这一条又强行拉回初始值，整个氢侧的时序就失真了。补一条首时段方程，或者干脆把H_tank(0)作为参数写进t=1的等式，两种写法都比留空强。

参数说明：燃料电池耗氢项除以效率而不是乘，因为效率是「氢→电」的转换比，同样的发电量，效率越低耗的氢越多。这一点写反了，模型会解出「越高效的燃料电池越费氢」这种自相矛盾的结果。分母里的 33.3 是氢低位热值，单位 kWh/kg，乘 1000 把兆瓦时换算成千瓦时，保证结果落在千克量纲上。如果你的氢数据来自高压气瓶按体积计量，还要再叠一层理想气体状态方程的换算，那就别在约束里硬写，单独封一个函数。

3. 阶梯碳交易的分段线性化：别在目标函数里写 if-else

3.1 为什么 if-else 一定会翻车

把手写的if emission <= 1000 ... elseif直接塞进目标函数，是这类模型里出现频率最高的错误。YALMIP 在做符号解析时，遇到以决策变量为条件的if分支，会直接抛出「非凸」或者「无法确定分支」的警告；更麻烦的是有些版本不报错，而是默默把你的表达式当成常量处理，于是碳成本恒等于某个固定值，模型解出来的调度方案看上去「合理」，但碳价机制形同虚设。这个坑我在早期项目里踩过一次，直到对账时发现碳成本占比恒定在 3% 才察觉。

正确的姿势是把阶梯结构拆成「分段 + 权重」。这里要先区分两种阶梯计费语义，它们对应的写法完全不同：

• 超额累进：每段只对超出上一档阈值的部分按本档单价计价，和居民阶梯电价一致。此时分段变量d_k表示落在第 k 段的排放量，总排放等于各段之和。
• 全额跳档：总排放一旦超过某阈值，全部排放按该档单价重新计价。这种语义才是真正需要 0-1 变量的。

绝大多数能源系统文献用的是第一种。而第一种在「成本最小化」的目标下，连续松弛就能得到正确解，不需要任何二进制变量。

3.2 辅助变量构造三段超额累进碳价

%% 阶梯碳交易：三段超额累进 segW = [200, 300, 500]; % 各段宽度 tCO2 segP = [ 80, 120, 160]; % 各段边际碳价 元/tCO2 d = sdpvar(3, T, 'full'); % d(k,t)：第 t 时段落在第 k 段的排放量 E_t = sdpvar(T, 1, 'full'); % 第 t 时段总排放 tCO2 for t = 1:T cons = [cons, E_t(t) == ef_gt * sum(P_gt(t,:)) * dt]; cons = [cons, sum(d(:,t)) == E_t(t)]; for k = 1:3 cons = [cons, 0 <= d(k,t) <= segW(k)]; end end carbon_cost = sum(sum(segP' .* d)); % 3×24 与 3×24 逐元素乘再求和 total_emission = sum(E_t);

逻辑说明：这段代码把每个时段的排放拆成三个非负变量，各自带一个上界（段宽）。因为segP严格递增，而目标函数是最小化碳成本，求解器为了让总成本更小，会本能地把d(1,t)填满再填d(2,t)，最后才动d(3,t)。这相当于自动实现了「先按低价段计费」的语义，不需要任何 0-1 变量。

参数说明：segW三段相加等于 1000 吨，意味着超过 1000 吨的部分在模型里没有变量承接，sum(d(:,t)) == E_t(t)会直接导致不可行。所以要么把最后一段上界设成一个足够大的数（但不能是 1e6，那样 LP 松弛会松得离谱），要么按最大可能排放量算一个紧上界：segW(3) = ef_gt * sum(P_gt_max) * dt * T。我一般取历史最大排放的 1.2 倍，既有余量又不至于毁掉松弛质量。

注意：如果换了求解器，或者把目标函数改成了最大化（比如做碳排放权交易的收益模型），上面「自动填满低段」的性质会失效，必须显式加二进制变量约束分段顺序。判断依据很简单：看你的目标函数和碳价方向是否一致。

3.3 SOS2 与二进制两种写法的取舍

除了辅助变量法，YALMIP 还支持sos2特殊有序集，适合处理连续分段线性函数（比如后面要讲的电解槽效率曲线）。两种方式的差别在实际跑的时候很明显：

工程上我基本只用前两种，二进制只留给真正需要离散决策的地方（设备启停）。把碳价这种纯经济分段也用 0-1 表达，等于白白把 MILP 的难度往上抬一档。

3.4 碳成本回代校验：算完必须手算一遍

模型能解出来不等于碳成本算对了。我习惯在求解后单独把value(d)拿出来，按段汇总，再和value(carbon_cost)对一遍，两分钟的事，能挡掉后面九成的诡异结论。

%% 碳成本回代校验 d_val = value(d); % 3×24 E_val = value(E_t); cost_chk = sum(sum(segP' .* d_val)); fprintf('模型目标中的碳成本: %.2f 元\n', value(carbon_cost)); fprintf('手工回代碳成本: %.2f 元\n', cost_chk); fprintf('总排放: %.1f tCO2，其中各段占比 %.1f%% / %.1f%% / %.1f%%\n', ... sum(E_val), 100*sum(d_val(1,:))/sum(E_val), ... 100*sum(d_val(2,:))/sum(E_val), 100*sum(d_val(3,:))/sum(E_val));

如果两个数对不上，八成是把d的维度搞反了，或者segP'的转置用错了位置——3×1转成1×3再和3×24逐元素乘，MATLAB 会做隐式扩展，但如果写成了segP .* d，维度不匹配会直接报错，反而是好事。

4. 电制氢与储氢罐：效率折算方向、索引偏移和启停次数

4.1 电解槽效率不是常数，分段线性化比二次拟合更稳

电解槽在低负载率区间的效率下降很明显，额定点附近反而平缓。有人喜欢用二次函数拟合，写起来短，但目标函数一旦带上二次项，整个问题就从 MILP 变成 MIQP，CPLEX 的求解时间可能翻好几倍。我一般用 3 到 4 个断点的分段线性，误差控制在 2% 以内，换来的是纯线性结构。

%% 电解槽效率分段线性化（SOS2） bp = [0, 0.3, 0.6, 1.0]; % 归一化负载率断点 eta = [0.60, 0.72, 0.71, 0.66]; % 对应电解效率 lam = sdpvar(4, T, 'full'); % 各断点的权重 for t = 1:T cons = [cons, sos2(lam(:,t))]; cons = [cons, sum(lam(:,t)) == 1]; cons = [cons, P_ele(t) == P_ele_max * sum(bp' .* lam(:,t))]; cons = [cons, H_prod(t) == k_ele * eta_curve_val(t) * P_ele(t)]; end

逻辑说明：sos2保证权重向量里最多两个相邻元素非零，这样P_ele必然落在两个断点之间，效率由这两个断点的eta线性插值得到。YALMIP 会把sos2转成求解器原生的 SOS 约束，CPLEX 原生支持，比手写 big-M 干净得多。

参数说明：断点必须包含 0 和 1，否则低负载段会外推。eta的取值取决于你的电解槽类型，碱性电解槽和质子交换膜电解槽的效率曲线形状不同，前者在低负载时衰减更明显。如果手头没有实测数据，按额定点效率上下浮动 10% 构造一条平滑曲线，比直接用常数效率更接近真实调度行为。

如果求解器不支持 SOS2，退而求其次用二进制选段，但断点数要压到 3 个以内，否则 0-1 变量数量乘以 24 个时段，规模会很可观。

4.2 储氢罐的容量约束与「负储能」现象

储氢罐的约束表面上看只有一条状态方程加一条上下限，但实际踩的坑集中在两处：一是容量单位，二是充放不同时。

%% 储氢罐容量与充放互斥 for t = 1:T cons = [cons, 0.1*H_cap <= H_tank(t) <= 0.9*H_cap]; end % 同一时段不允许既制氢又发电（可选，视设备配置） for t = 1:T cons = [cons, P_ele(t) <= P_ele_max * u_ele(t)]; cons = [cons, P_fc(t) <= P_fc_max * (1 - u_ele(t))]; end

逻辑说明：上下限取 10% 到 90% 而不是 0 到 100%，是因为真实储氢罐需要留缓冲压力，同时这也是个建模技巧——把边界收紧能显著改善 LP 松弛。互斥约束是否需要，要看你的系统里电解槽和燃料电池是不是共用一套管路，如果是分立的设备，这条约束可以去掉，模型会宽松很多。

参数说明：H_cap的单位必须和状态方程里的千克一致。见过有人把储氢罐按标方（Nm³）标定，直接填进以千克为单位的方程，结果容量约束形同虚设，储氢量能冲到物理上不可能的数值。换算关系大致是 1 kg 氢气约 11.1 Nm³，标定前先统一。

4.3 启停次数约束：不加这条，调度方案在工程上没法执行

纯经济优化会给出一个「电解槽全天在 0 和满功率之间反复横跳」的方案，成本确实低，但设备寿命撑不住。加启停约束的写法有两种，一种限制总启停次数，一种加最小运行/停机时间。前者简单，后者更贴近实际，我一般先用前者把模型跑通。

%% 启停次数约束 startup = sdpvar(T, 1, 'full'); shutdown = sdpvar(T, 1, 'full'); for t = 2:T cons = [cons, startup(t) >= u_ele(t) - u_ele(t-1)]; cons = [cons, shutdown(t) >= u_ele(t-1) - u_ele(t)]; cons = [cons, startup(t) >= 0, shutdown(t) >= 0]; end cons = [cons, sum(startup) <= 3]; % 24 小时内最多启停 3 次 cons = [cons, P_ele >= r_ele_min * P_ele_max * u_ele]; % 最小出力耦合

逻辑说明：startup(t) >= u(t) - u(t-1)这种写法是「大等于」而不是「等于」，配合目标函数里给startup一个正的小系数（比如 1 元/次），求解器会自动把它压到最紧，不会虚高。如果目标函数里完全不给startup系数，也没有sum(startup)上界，这些变量就是自由的，约束等于没加。

参数说明：r_ele_min取 0.3 是工程经验值，低于这个负载率电解效率下降快，产氢纯度也受影响。这个系数直接决定了P_ele的可行域，如果设成 0.5，模型的可行解空间会明显收窄，某些极端风光场景下可能直接不可行——这时候要么放宽系数，要么允许弃风。

5. 避坑与排查：五个让模型跑不出结果或结果离谱的现场

5.1 求解器与数值层面的坑

坑一：约束全对，但 CPLEX 报 infeasible。现象是optimize返回problem = 1，diagnostics显示不可行，逐条检查却看不出问题。原因通常出在储氢罐的日闭合约束上——H_tank(T) == H0要求一天的净制氢量等于零，而风光出力曲线和负荷曲线本身不自洽，物理上就做不到。解决方法是把日闭合改成软约束，加一个松弛变量进目标函数并给高惩罚，或者把H0放开成可优化变量，只在年终做一次总量平衡。我一般选后者，因为日调节本来就不该强求每日归零。

坑二：big-M 拍得太大，求解时间爆炸。现象是模型规模看着不大，CPLEX 却跑了二十分钟还没收敛，日志里分支定界节点数疯涨。原因是P_ele <= 1e6 * u_ele这种写法让 LP 松弛几乎无效，求解器只能靠分支硬搜。解决方法是把 M 换成物理紧上界：P_ele_max，如果还有时段耦合，就再和该时段的负荷上界取小。M 每缩小一个数量级，节点数往往能掉一大截。

坑三：目标函数里出现min/max函数，模型被判定为非线性。现象是 YALMIP 提示切换到了非线性求解器，或者直接警告Non-convex。原因是min(total_emission, 200)里的比较对象是决策变量，符号引擎无法把它展成线性形式。解决方法是统一改写成辅助变量的形式，用上下界和等式约束替代，让表达式在符号层面保持仿射。

坑四：mipgap 设成 0.01，两次运行结果不一样。现象是同样的模型和数据，重新跑一遍总成本差了几百块。原因很简单，gap 是相对最优界定义的，成本在十万量级时 1% 就是一千块，求解器在达到这个精度后就停了。解决方法是出报告时把实际 gap 打印出来，涉及方案对比的场景压到 0.001，同时在sdpsettings里固定随机种子相关的分支策略参数。

5.2 建模与结果校核层面的坑

坑五：储氢量出现负值。现象是value(H_tank)里有负数，但约束明明白白写了H_tank >= 0.1*H_cap。原因有两个，一是效率折算方向反了——燃料电池那侧本该除以效率，写成了乘，导致耗氢被低估，库存被高估后又被别的约束拉扯；二是上下限约束写在循环外面，只约束了最后一个时段。解决方法是在求解后立刻对H_tank做一次min()检查，同时在建模阶段就用向量化写法一次性约束所有时段，别在循环里漏。

坑六：碳成本占比恒定不变。现象是不同风光场景下，碳成本在总成本里的比例几乎一模一样，看起来「很稳定」。原因大概率是分段变量被写成了常量，或者segP' .* d里的转置写错导致广播成了标量乘。解决方法就是前面那段回代校验，两行fprintf就能暴露。

坑七：电解槽启停次数远小于约束上限，但调度方案依然抖。现象是加了sum(startup) <= 3之后，P_ele曲线还是在 30% 和 100% 之间来回跳。原因是u_ele是连续松弛出来的——如果你不小心用了intvar而不是binvar，并且没配0<=u<=1，求解器可能给出 0.5 这样的值，P_ele >= 0.3*P_max*u就变成了一个空约束。解决方法很直接，检查value(u_ele)的取值范围，出现非整数就是这个问题。

6. 结果怎么读与怎么验：三条曲线、一次扫描

6.1 三条必须亲手画一遍的曲线

模型跑完之后，出图不是走流程，而是校核。我固定画三张：第一张是各单元出力堆叠图，横轴 24 小时，重点看风电大发时段电解槽有没有跟上，以及火电曲线在碳价跳档时刻有没有出现拐点；第二张是储氢量曲线，看它是不是在容量区间内平滑波动，如果出现「贴着上下限走」的锯齿，说明容量约束太紧或者互斥约束没生效；第三张是分段碳排放的堆叠柱状图，直接看每一档用了多少吨，第三段一旦被触发，就说明碳价机制真的开始起作用了。

%% 碳价敏感性扫描：找阶梯被触发的临界点 mult_list = 0.5:0.25:2.5; res = zeros(numel(mult_list), 3); for i = 1:numel(mult_list) segP_i = segP * mult_list(i); carbon_cost_i = sum(sum(segP_i' .* d)); obj_i = sum(sum(gen_cost .* P_gt)) + sum(pen_curt .* P_curt) + carbon_cost_i; sol = optimize(cons, obj_i, ops); if sol.problem == 0 res(i,:) = [value(obj_i), value(total_emission), value(carbon_cost_i)]; end end

这段扫描把碳价整体放大 0.5 到 2.5 倍，每档记录总成本、总排放和碳成本。mult_list的步长取 0.25 是我试出来的经验值，再细就只是在重复同样的解，再粗会跨过拐点。跑完把res画成双纵轴图，通常能看到一条明显的分界：碳价倍数低于某个值时，排放基本不变，因为减排成本高于碳成本；跨过这个点之后，排放陡降，火电出力被压到最低，这就是碳价开始真正改变调度行为的临界区。工程上拿到这个临界点，比拿到一个具体的最优解有用得多。

ops里的mipgap在这一步建议压到 0.002，因为扫描的目的是比较不同碳价下的相对变化，不是求单个场景的精确最优，但如果 gap 太大，相邻两个倍数之间的差异会被求解误差淹没，曲线上看不到台阶。扫描耗时大概是单次求解的 9 倍，40 MW 规模的系统单次 30 秒左右，整体五分钟能跑完，值得等。

6.2 从那以后的习惯

早期做这类模型，我最容易犯的错是「解出来了就当对了」，直到有一次拿结果去对账，发现碳成本比手算少了两个数量级，回头查才发现是分段变量的维度广播出了问题。从那以后我每次跑完阶梯碳交易模型，都强制走一遍「回代校验 — 储氢量最小值检查 — 分段占比打印」这三步，缺一步不出图。这套流程听起来笨，但它挡掉的返工比任何调试技巧都多。希望帮到你。

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