企业官网建设流程全解析

简介：面向能源系统优化研究人员与具备MATLAB基础的技术人员，这份PDF聚焦阶梯型碳交易与电制氢综合能源系统优化调度，解决碳价分段跃升、电-氢-电多向转换下机组出力安排、成本控制与碳减排问题。内容涵盖传统机组、风电、电解槽、储氢罐、燃料电池五单元建模，阶梯碳交易分段线性化与big-M处理、储氢动态方程、弃风惩罚及目标函数设计，并用24小时算例展示风电大发时优先制氢、碳排放接近阈值时火电调整等现象；同时提供CPLEX+YALMIP求解设置、结果可视化、关键注释、参考文献及big-M数值调试思路。包内共1个PDF，约129KB，便于集中阅读与对照建模。已有163人学习，适合希望理解综合能源系统混合整数线性规划建模细节、复现阶梯碳价与氢能耦合调度案例的读者参考。

1. 阶梯碳价加电制氢：为什么调度表会在凌晨两点翻脸

凌晨两点风电突然抬高，调度表却把电解槽压到最低，碳成本还跳到了第二档——这是我在调 MATLAB+YALMIP 综合能源系统模型时最常遇到的“假合理”结果。不是求解器算错了，而是阶梯碳交易的分段变量或电制氢的时间索引有一处没对齐，模型依然能给出光滑曲线，但成本已经悄悄偏了。这份代码包围绕传统机组、风电、电解槽、储氢罐、燃料电池五个单元，把阶梯碳交易和电制氢塞进 24 小时混合整数线性规划里，交付的是sys_init.m、build_model.m、plot_results.m三个模块。适合已经写过 YALMIP 小模型、准备把碳成本与氢储能耦合进调度的人；新手可以照着参数表改数，熟手重点看分段线性化和索引闭合那几段。

2. 从设备参数到决策变量：24小时调度模型的骨架怎么搭

2.1 五个单元的变量与功率平衡

综合能源系统调度最怕一上来就写目标函数，变量没定义清楚，后面碳排和储氢全对不上。我一般先把五个单元的功率、储氢状态、启停状态和碳排分段变量一次性摆出来，再写功率平衡。传统机组按多台机组处理，风电拆成“可用出力—实际消纳—弃风”三部分，电解槽和燃料电池分别作为电负荷和电源进入平衡式。这样写的好处是弃风惩罚有明确载体，后续阶梯碳交易也能直接挂到机组出力上。

% sys_init.m 参数入口（示例值，按本地数据替换） T = 24; % 调度时段数 N_gen = 2; % 传统机组台数 P_load = [120 115 110 108 112 125 150 170 180 185 190 188 ... 186 182 178 180 190 205 215 210 195 175 150 130]; % MW P_wind_avail = [95 110 120 125 118 100 85 70 60 55 50 48 ... 52 58 65 72 80 88 92 90 85 80 75 70]; % MW P_gt_max = [180 150]; % 机组出力上限 MW P_gt_min = [60 50]; % 机组出力下限 MW em_factor = [0.85 0.92]; % 机组排放因子 t/MWh

代码块里P_load和P_wind_avail是典型日曲线，单位统一为 MW；em_factor决定后续碳排计算，不要和燃料成本系数混用。P_gt_min不能设成 0，否则机组会频繁在极小出力下“蹭”平衡，碳排曲线看起来平滑但机组效率失真。实际项目里我习惯把负荷和风电放在独立 CSV，sys_init.m只负责读表和做单位检查。

% build_model.m 变量定义 P_gt = sdpvar(T, N_gen, 'full'); % 传统机组出力 P_wind_use = sdpvar(T, 1); % 风电实际消纳 P_curt = sdpvar(T, 1); % 弃风功率 P_ele = sdpvar(T, 1); % 电解槽耗电 P_fc = sdpvar(T, 1); % 燃料电池发电 H = sdpvar(T, 1); % 储氢罐等效氢量 u_ele = binvar(T, 1); % 电解槽运行状态 u_fc = binvar(T, 1); % 燃料电池运行状态 v_ele = binvar(T, 1); % 电解槽启动动作

变量定义阶段要守住两个边界：P_curt必须非负，H必须给上下限，否则求解器会在可行域边缘给出负储能或负弃风。u_ele、u_fc是设备启停状态，v_ele只表示启动动作，后面用v_ele(t) >= u_ele(t) - u_ele(t-1)关联。很多翻车案例不是约束写错，而是变量没设'full'或维度少一列，YALMIP 自动扩展后索引悄悄错位。

功率平衡式要同时覆盖电平衡和风电场内部分解。电解槽耗电放在负荷侧，燃料电池发电放在电源侧，储氢罐不直接参与电平衡，只通过电解槽和燃料电池间接耦合。

cons = []; for t = 1:T cons = [cons, sum(P_gt(t,:)) + P_wind_use(t) + P_fc(t) == ... P_load(t) + P_ele(t)]; cons = [cons, P_wind_use(t) + P_curt(t) == P_wind_avail(t)]; cons = [cons, P_curt(t) >= 0, P_wind_use(t) >= 0]; end

第一行是电功率平衡，sum(P_gt(t,:))表示两台机组同时刻出力之和；第二行把风电可用出力拆成消纳和弃风，保证风电不会凭空多出来。P_curt只设非负还不够，最好再给P_curt(t) <= P_wind_avail(t)，否则极端参数下会出现“弃风大于可用风电”的数学解。这个约束在正常数据里几乎不触发，但一旦风电预测文件有负值或单位写成 kW，它就会暴露问题。

2.2 目标函数的五项成本与线性化边界

目标函数不要只写燃料成本，阶梯碳交易和电制氢进来后，成本构成会改变调度顺序。我通常拆成五项：燃料成本、机组运维、弃风惩罚、碳成本、启停成本。燃料成本用线性系数，机组出力乘成本系数；运维成本按各设备出力线性累加；弃风惩罚要足够大，让模型优先消纳风电，但也不能大到让碳成本失去意义。碳成本单独由分段变量生成，启停成本只挂在电解槽启动动作上。

% 成本系数（示例值） c_fuel = [28 31]; % 元/MWh c_om_gt = 12; % 机组运维 元/MWh c_om_ele = 8; % 电解槽运维 元/MWh c_om_fc = 10; % 燃料电池运维 元/MWh c_curt = 120; % 弃风惩罚 元/MWh c_start_ele = 300; % 电解槽单次启动成本 元/次 fuel_cost = sum(sum(P_gt .* c_fuel)); om_cost = c_om_gt*sum(P_gt(:)) + c_om_ele*sum(P_ele) + c_om_fc*sum(P_fc); curt_cost = c_curt*sum(P_curt); start_cost = c_start_ele*sum(v_ele);

c_fuel按机组列向量乘到P_gt上，YALMIP 会自动按列广播；sum(P_gt(:))把所有时段所有机组拉平，避免维度不一致。弃风惩罚c_curt我一般设在燃料成本的 3 到 5 倍之间，太小模型会宁愿弃风也不让电解槽多耗电，太大又会让碳成本分段变量被压制。启停成本不要设得比单次启动实际损耗低太多，否则电解槽会在凌晨反复启停，图形上看着“灵活”，设备寿命账算不回来。

碳成本后面单独用分段变量拼，不要在这里直接写if-else。目标函数最终形式是五项相加，其中碳成本是一个标量，启停成本也是标量，求和时注意不要重复计入储氢罐的运维，储氢罐通常只算在电解槽和燃料电池的运维系数里。

2.3 约束的索引与初始末端条件

储氢罐动态是这套模型最容易翻车的地方。电解槽在t时段耗电产氢，燃料电池在t时段耗氢发电，储氢状态从t-1过渡到t。如果循环写成for t = 1:T然后直接引用H(t-1)，MATLAB 在t=1会报错；如果写成for t = 2:T但功率用了P_ele(t)，又会出现最后一个时段状态不闭合。我的习惯是固定H(1)为初始储氢，循环从t=2开始，方程右侧统一用t-1时段的功率。

H_init = 80; % 初始储氢量 kg H_min = 20; H_max = 300; % 储氢上下限 kg k_ele = 21.6; % 电解槽产氢系数 kg/MWh k_fc = 0.0193; % 燃料电池耗氢系数 MWh/kg cons = [cons, H(1) == H_init]; for t = 2:T cons = [cons, H(t) == H(t-1) + k_ele*P_ele(t-1) - P_fc(t-1)/k_fc]; cons = [cons, H_min <= H(t) <= H_max]; end cons = [cons, H(T) == H_init]; % 日末回到初始，避免跨日套利

k_ele的量纲是 kg/MWh，k_fc的量纲是 MWh/kg，所以燃料电池项是P_fc/k_fc，不是乘效率。这里用等效氢量把 LHV 和效率折进去，代码可读性好，但参数表必须写清折算依据，否则换一组氢负荷数据就会对不上。H(T) == H_init是日调度常用闭合条件，防止模型把储氢罐放空来降低当日成本。如果做连续多日滚动，末端条件要改成软约束或终端价值函数，不要硬套日末相等。

机组出力上下限和爬坡约束也要一起给。P_gt_min和P_gt_max按机组列向量约束，爬坡限制用相邻时段差值。没有爬坡约束时，阶梯碳价会让机组在阈值附近瞬间切换出力，曲线出现不现实的垂直跳变。

ramp_up = [60 50]; ramp_down = [60 50]; for t = 1:T cons = [cons, P_gt_min <= P_gt(t,:) <= P_gt_max]; end for t = 2:T cons = [cons, P_gt(t,:) - P_gt(t-1,:) <= ramp_up]; cons = [cons, P_gt(t-1,:) - P_gt(t,:) <= ramp_down]; end

ramp_up和ramp_down单位是 MW/h，时段步长 1 小时时可直接用差值。若数据是 15 分钟粒度，要乘 0.25 或把系数改成 MW/15min。这个细节不写进参数表，后面换数据的人一定会踩坑。

3. 阶梯碳交易的分段线性化：big-M、SOS2与碳成本跃迁

3.1 先判断凸性：递增阶梯不需要big-M

阶梯碳交易最容易被写复杂。很多人一看到“阶梯”就上 big-M 和二进制变量，结果模型规模翻倍，CPLEX 数值稳定性下降。先判断阶梯结构：如果碳价随排放量递增，且每一段只对超出部分计价，这是凸分段线性成本，直接用非负辅助变量表示各段排放量即可，不需要二进制。优化目标是最小化成本，价格低的段会被自然填满，价格高的段只在总排放超过阈值时启用。

% 碳交易参数：分段阈值与单价（示例） E_cap = [200; 300; 1e4]; % 各段容量 t，第三段给足够大上界 p_carbon = [80; 120; 160]; % 元/t d = sdpvar(3,1); % 各段碳排放量 E_total = sum(sum(P_gt .* em_factor)); % 调度周期总排放 t cons = [cons, d >= 0]; cons = [cons, d(1) <= E_cap(1)]; cons = [cons, d(2) <= E_cap(2)]; cons = [cons, d(3) <= E_cap(3)]; cons = [cons, sum(d) == E_total]; carbon_cost = p_carbon' * d;

E_cap不是阈值本身，而是每段可容纳的排放量。第一段 200 t，第二段 300 t，第三段给 10000 t 上界，防止高排放场景不可行。p_carbon递增，所以d(1)会先被填满，再填d(2)，最后才用d(3)，不需要额外顺序约束。E_total是所有时段所有机组排放之和，是标量；如果按小时分别计阶梯，d要变成3×T，成本再按小时求和，那是另一套机制，不要混。

这里的关键不是代码长短，而是选型理由：凸分段线性化保持模型为 MILP，甚至可以是 LP（如果没有启停二进制），求解速度快，碳价跃迁通过分段变量自然体现。big-M适合处理逻辑条件，不适合处理这种“价格递增、按段计价”的成本函数。把凸分段硬写成 big-M，等于给求解器增加不必要的分支，血泪经验是模型能跑但 gap 降不下来。

3.2 非凸阶梯与SOS2的适用边界

如果碳交易规则不是“超出部分按高价”，而是“一旦超过阈值，全部排放按高价重算”，成本函数出现跳变，凸性被破坏，这时才需要二进制或 SOS2。比如排放 201 t 时，前 200 t 也按第二档价格计费，总成本不是连续分段函数。YALMIP 里可以用 SOS2 构造分段线性函数，也可以用二进制变量选择工作区间。代码包的主模型采用凸分段，非凸版本只保留接口，避免默认模型过重。

% 非凸阶梯示意：用二进制选择区间，仅作边界说明 z = binvar(3,1); M = 1e4; cons = [cons, sum(z) == 1]; cons = [cons, E_total <= 200 + M*(1-z(1))]; cons = [cons, E_total >= 200*z(2), E_total <= 500 + M*(1-z(2))]; cons = [cons, E_total >= 500*z(3)]; carbon_cost = 80*200*z(1) + 120*E_total*z(2) + (120*500 + 160*(E_total-500))*z(3);

这段只说明非凸阶梯需要区间选择，M不能随手写 1e6。M应取E_total的物理上界，比如最大机组出力乘排放因子乘时段数。M过大是 CPLEX 数值敏感的主要来源之一，后面避坑章会展开。实际项目里，如果政策规则是“超额累进”，优先用 3.1 的凸分段；如果规则是“整段跳价”，再考虑 SOS2 或二进制，并且先做小规模测试。

3.3 碳排计算与单位换算

碳排计算本身不复杂，难的是单位。机组出力是 MW，时段是小时，排放因子是 t/MWh，乘出来才是 t。如果负荷数据是 kW，排放因子是 kg/MWh，最后会差 1000 倍，碳成本看起来像 0，阶梯永远不触发。我在sys_init.m里加了一个单位检查：负荷峰值、风电峰值、机组容量三者量级要一致，碳排总量要落在碳段阈值附近，否则先查单位。

碳成本算完后，我习惯用value(E_total)和value(d)复算一遍：value(d(1))+value(d(2))+value(d(3))必须等于value(E_total)，p_carbon'*value(d)必须等于目标函数里的碳成本项。这两个等式不成立，说明分段变量或目标函数拼错了。图形上碳成本曲线在阈值处有拐点，但拐点不等于跳变；如果看到碳成本垂直跳，先检查是不是把非凸规则误写进了凸分段。

4. 电制氢与储氢罐：效率折算、时间索引和启停约束

4.1 电解槽与燃料电池的效率折算

电制氢模块最容易写成“效率乘功率”就完事，实际单位换算比代码长。电解槽耗电产氢，产氢速率等于电功率乘电解效率除以氢热值；燃料电池耗氢发电，耗氢速率等于电功率除以燃料电池效率再除以氢热值。代码里把热值和效率折成k_ele、k_fc两个系数，储氢方程里一个乘、一个除，方向不能反。k_ele用 21.6 kg/MWh 时，1 MWh 电大约产 21.6 kg 氢；k_fc用 0.0193 MWh/kg 时，1 kg 氢大约发 0.0193 MWh 电，效率损失体现在系数里。

% 电解槽与燃料电池约束 P_ele_max = 80; P_ele_min_ratio = 0.3; P_fc_max = 50; for t = 1:T cons = [cons, 0 <= P_ele(t) <= P_ele_max*u_ele(t)]; cons = [cons, P_ele(t) >= P_ele_min_ratio*P_ele_max*u_ele(t)]; cons = [cons, 0 <= P_fc(t) <= P_fc_max*u_fc(t)]; end

P_ele_min_ratio = 0.3表示电解槽运行下限为额定功率的 30%，低于这个值效率下降明显，常见做法是用实验曲线或厂家数据替换常数。u_ele为 1 时电解槽才允许出力，为 0 时功率强制为 0。燃料电池没有设最小负荷，因为它的调节范围通常更宽；如果实际设备有最小出力，同样加P_fc(t) >= P_fc_min*u_fc(t)。这里不要用P_ele(t) <= P_ele_max代替二进制约束，否则启停次数统计会失真。

4.2 储氢罐动态与末端闭合

储氢罐是时间耦合的核心。H(t)等于上一时段储氢量加电解产氢、减燃料电池耗氢，再减氢负荷。氢负荷可以设为固定工业用氢或加氢站需求，没有氢负荷时置 0。代码里H(1)固定为初始值，循环从t=2开始，右侧用t-1的功率，保证每个时段的状态只由上一时段决定。末端H(T) == H_init是日调度闭合条件，防止模型把储氢放空。

H_demand = zeros(T,1); % 氢负荷 kg/h，示例置零 for t = 2:T cons = [cons, H(t) == H(t-1) + k_ele*P_ele(t-1) ... - P_fc(t-1)/k_fc - H_demand(t-1)]; cons = [cons, H_min <= H(t) <= H_max]; end

如果直接把H(t)写成H(t) = H(t-1) + k_ele*P_ele(t) - P_fc(t)/k_fc，功率和状态同一时段，会导致“先产氢再储氢再发电”的代数环，求解器可能给出同时充放的非物理解。H_demand的单位是 kg/h，和k_ele*P_ele一致；如果氢负荷是 Nm³/h，要先换算成 kg/h。末端闭合不是绝对必要，但在单日调度里能避免“今天放空、明天补回”的跨日套利，做滚动调度时改成终端软约束更合理。

4.3 启停次数与最小负荷

电解槽启停次数约束是电制氢模型里最容易被忽略的寿命成本。没有启停约束时，模型会在电价或碳价拐点处频繁开关电解槽，调度曲线看起来“灵活”，实际设备受不了。我一般用二进制u_ele表示运行状态，v_ele表示启动动作，约束v_ele(t) >= u_ele(t) - u_ele(t-1)，再限制 24 小时内启动次数不超过 3 次。启动成本挂在v_ele上，目标函数自然会权衡。

for t = 2:T cons = [cons, v_ele(t) >= u_ele(t) - u_ele(t-1)]; cons = [cons, v_ele(t) >= 0]; end cons = [cons, v_ele(1) >= u_ele(1)]; cons = [cons, sum(v_ele) <= 3]; cons = [cons, start_cost == c_start_ele*sum(v_ele)];

v_ele(t) >= u_ele(t) - u_ele(t-1)只保证启动时v_ele为 1，不保证停机时v_ele为 0，但目标函数有启动成本，优化会尽量让v_ele取小，所以不需要额外等式。sum(v_ele) <= 3是 24 小时启动次数上限，按设备手册或运维经验改。v_ele(1)单独处理初始状态，如果初始时刻电解槽已运行，v_ele(1)应设为 0，否则会多算一次启动。这个小细节在日滚动调度里经常导致首时段启动成本虚高。

电解槽和燃料电池同时运行在物理上并非绝对禁止，但效率损失会让它在经济上不划算。如果结果里出现两者同时满出力，先查氢负荷和末端储氢约束，再考虑加u_ele(t) + u_fc(t) <= 1的互斥约束。互斥约束会增强整数性，可能降低求解速度，所以不要默认加，先看结果是否合理。

5. 求解器避坑排查：CPLEX数值敏感、索引错位与负储能

5.1 先看可行性：功率平衡与储氢闭合

现象：模型返回Infeasible，但去掉碳成本后又能求解。
原因：碳段容量上界设得太小，总排放超过E_cap总和，或者H(T) == H_init与氢负荷、设备容量冲突。
解决：先跑无碳价基准，确认功率平衡和储氢闭合可行，再把E_cap(3)放到物理上界，末端储氢改成软约束试算。我习惯在build_model.m里加一个checks结构体，求解前打印总负荷、总风电、最大排放和储氢容量，数字对不上就不进求解器。

现象：求解成功，但H出现负值或最后时段不闭合。
原因：循环索引从t=1开始引用H(t-1)，或者功率用t而状态用t-1，导致最后一个时段功率没进状态方程。
解决：固定H(1)=H_init，循环t=2:T，右侧统一t-1；求解后检查value(H)最小值和value(H(T))-H_init。这两个数不满足就回头查索引，不要先调求解器参数。

5.2 再看数值：big-M与MIP gap

现象：CPLEX 长时间不收敛，日志里出现数值困难或分支定界树爆炸。
原因：非凸阶梯用了过大的M，或者凸分段误加了二进制变量；也有可能是机组爬坡约束和启停约束同时过紧。
解决：凸分段去掉二进制，非凸阶梯的M用物理上界动态计算，比如M = max(P_gt_max)*max(em_factor)*T；先设mipgap=0.01跑通，再逐步收紧到 0.005。不要一上来追求 1e-4，工程调度 1% 的 gap 通常够用。

现象：碳成本为 0 或阶梯完全不触发。
原因：排放因子单位写成 kg/MWh，而碳段阈值按 t 设，量级差 1000 倍；或者E_total用了单时段排放而非全周期总排放。
解决：打印value(E_total)和E_cap对比，量级不对先查单位；全周期阶梯用sum汇总，不要放在时段循环里逐时判断。

5.3 最后看设备逻辑：启停与同时运行

现象：电解槽 24 小时内启停十几次，启动成本却不高。
原因：只约束了功率上下限，没有约束v_ele和启动次数；或者启动成本系数远低于实际损耗。
解决：加v_ele(t) >= u_ele(t)-u_ele(t-1)和sum(v_ele)<=3，启动成本按设备手册或运维统计设，不要用 0 或 1 元这种象征值。

现象：电解槽和燃料电池同一时段同时满出力，储氢量却平稳。
原因：效率系数方向写反，或者k_ele、k_fc单位不一致，导致“耗氢发电”被算成“产氢”。
解决：用手算校验：1 MWh 电产多少 kg 氢，1 kg 氢发多少 MWh 电，代进方程看储氢变化方向；必要时加u_ele(t)+u_fc(t)<=1做互斥测试，确认不是效率参数问题后再决定是否保留。

6. 结果验证与参数扫掠：用成本构成图判断模型是否可信

6.1 三个复算指标

求解器返回exitflag=1不代表结果可信。我每次先复算三个指标：功率平衡最大误差、储氢闭合误差、碳成本复算误差。功率平衡误差用value(sum(P_gt,2)+P_wind_use+P_fc-P_load-P_ele)的最大绝对值；储氢闭合误差看value(H(T))-H_init；碳成本复算用p_carbon'*value(d)和目标函数里的碳成本项对比。三个误差都在 1e-6 量级，才继续看曲线。

% 结果复算（放在 optimize 之后） P_balance_err = max(abs(value(sum(P_gt,2) + P_wind_use + P_fc - P_load - P_ele))); H_close_err = abs(value(H(T)) - H_init); carbon_recalc = p_carbon' * value(d); fprintf('功率平衡最大误差: %.2e MW\n', P_balance_err); fprintf('储氢闭合误差: %.2e kg\n', H_close_err); fprintf('碳成本复算: %.2f 元\n', carbon_recalc);

P_balance_err大于 1e-4 时，先查单位，再查是否漏了某台机组；H_close_err大于 1e-4 时，查末端约束和索引；carbon_recalc与目标函数碳成本项不一致时，查d的维度和E_total的汇总方式。这三个检查不通过，后面的成本构成图没有意义。

6.2 参数扫掠表

单点结果容易让人误判。我一般做三组扫掠：碳价整体上浮 20%、电解效率从 0.65 到 0.78、弃风惩罚从 80 到 160。每组只改一个参数，记录总成本、碳成本占比、弃风率和电解槽启动次数。下面是一组示例结果，用来判断趋势，不是标准答案。

如果碳价上浮后碳成本占比不升反降，先查碳段变量是否被目标函数正确引用；如果电解效率提高但总成本没降，查k_ele是否同步更新；如果弃风惩罚加大后启动次数暴涨，说明启停约束或启动成本太松。参数扫掠的目的不是找“最优参数”，而是看模型响应是否符合物理直觉。

6.3 一个具体技巧：先跑无碳价基准

我现在每次拿到新的综合能源系统数据，先注释掉碳成本项，跑一个无碳价基准，记录机组出力、弃风、储氢曲线；再把碳价设为常数跑第二遍；最后才切到阶梯碳交易。这样一旦阶梯版本出问题，能快速判断是碳分段写错，还是电制氢模块本身不可行。这个习惯来自一次翻车：阶梯碳价和储氢末端约束同时改，结果模型不可行，查了半天才发现是E_cap(3)给得太小。从那以后我每次改模型都强制走一遍“无碳价—固定碳价—阶梯碳价”的三步对照，先确认每个模块的边际影响，再调参数。希望帮到你。

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