企业官网建设流程全解析

简介：面向能源系统优化研究者及具备MATLAB基础的技术人员，这份PDF聚焦阶梯型碳交易与电制氢综合能源系统的优化调度，帮助读者理解传统机组、风电、电解槽、储氢罐和燃料电池五单元协同下的经济性与减排平衡。内容以YALMIP调用Cplex求解为主线，覆盖阶梯碳成本分段线性化与big-M处理、储氢罐动态状态方程、电-氢-电双向效率折算、目标函数设计以及结果可视化，并结合案例分析风电大发时段优先制氢、碳排放接近阈值时火电出力调整等现象。整份文档围绕模型搭建、参数初始化与结果解读组织，关键约束和易错细节有较充分说明，便于对照复现。压缩包共1个PDF文件，大小129KB；已有162人学习，适合参考约束写法、求解设置与数值稳定性处理细节，也可用于能源优化课题复现、课程设计及论文实验参考。

1. 一套会算碳账的调度器：MATLAB 里跑通阶梯碳交易与电制氢

如果只按固定碳价算成本，综合能源系统的调度结果往往很温和，火电多出一点力，交点碳费就过去了。可一旦把碳价改成阶梯制，再塞进一整套电制氢链路，模型的行为就会突然变得很不线性：某个时段的碳排放只要跨过阈值，后面的边际成本立刻跳档，而电解槽、储氢罐、燃料电池又要求你在时间轴上前后对齐。这份 MATLAB + YALMIP 的代码包干的就是这件事：传统机组、风电、电解槽、储氢罐、燃料电池五个单元，24 点逐时调度，目标是最小化燃料成本、弃风惩罚、运维与阶梯碳成本之和，约束交给 CPLEX 求解。它适合已经能看懂优化模型、想把碳机制和氢储能真正落进代码里的开发者，也适合要改造现有调度程序、把碳成本从常数项变成分段函数的人。

2. 系统架构拆解：五单元拓扑与三个 .m 文件的分工

2.1 能量流拓扑：母线、弃风变量与氢的往返路径

整个系统挂在一条电母线上，供电侧有三台可调火电、一片风电场和一个燃料电池，用电侧是固定负荷加一台电解槽。风电按预测曲线出，但不一定全被消纳，所以模型里专门留了一个弃风变量；电解槽是可控负荷，耗电产氢，氢进储氢罐，储氢罐再给燃料电池供氢，燃料电池又回到母线上发电。这条电极到氢、氢再到电的往返路径，就是整个模型区别于普通机组组合问题的地方。功率平衡约束写出来是：

% 母线功率平衡：发电侧 = 用电侧（每个时段一条约束） for t = 1:T cons = [cons, ... sum(P_gen(t,:)) + (P_wind(t) - P_curt(t)) + P_fc(t) ... == P_load(t) + P_ele(t)]; end

这里的逻辑是「火电总和 + 实际消纳风电 + 燃料电池出力 = 负荷 + 电解槽耗电」。参数说明：P_wind(t)是 t 时段风电可用出力预测值，属于常量；P_curt(t)是弃风功率，是决策变量，上行被P_wind(t)卡住，不能凭空多弃；P_fc(t)是燃料电池发电功率。把弃风写成变量而不是简单丢弃，是为了在目标函数里给它一个惩罚价，模型才知道该尽量消纳。设备参数统一放在初始化结构体里：

需要强调单位：储氢罐里存的不是 kg，而是「等效电量 MWh」，也就是罐里那部分氢如果拿去发电最多能出多少电。这样一来电解槽乘效率、燃料电池除效率就能直接加减，不用来回换算热值。建模时所有量统一到 MW 和 MWh 上，是后面不出现诡异负储氢的前提。

2.2 三文件分工：sys_init、build_model、plot_results 各管一段

代码拆成三个文件，边界划得很清楚。sys_init.m只做一件事：把上面那张表拼成一个units结构体，再生成 24 点负荷曲线和风电预测曲线，顺手算出几个派生量，比如总装机、M 值的上界、储氢容量对应的最大连续制氢时长。build_model.m负责全部建模和求解，输入units，输出解结构体sol和一个诊断结构体diag（求解状态、耗时、各分项成本）。plot_results.m只读不写，拿sol画三张图：出力堆叠图、储氢量与碳排分档曲线、成本构成饼图。调用骨架就这么几行：

clear; clc; yalmip('clear'); % 清空 YALMIP 内部缓存，避免上次的模型残留 units = sys_init(); % 生成参数结构体 + 24 点预测数据 [sol, diag]= build_model(units); % 建模、求解、成本拆分 plot_results(sol, units); % 出力曲线 / 氢罐曲线 / 碳档位 if diag.exitflag ~= 1 warning('求解未正常收敛，code = %d，先查第 5 章的排查清单', diag.exitflag); end

这样拆的好处是二次开发时改动面可控：换数据只动sys_init.m，加约束只动build_model.m，出图样式完全不碰。diag里我习惯把碳成本、燃料成本、弃风成本、启停成本四项单独存一份，因为调试阶段最常见的错误就是总成本算错了却不知道错在哪一项，有了分解就能一眼看出是不是碳成本被放大了 1000 倍这种低级事故。

2.3 决策变量清单：维度错一次，后面全乱

变量定义写在build_model.m开头，维度必须一次定对。火电用T×N_gen的矩阵，其余时序量用T×1列向量，启停状态用binvar。YALMIP 里最容易翻车的地方是sum的默认行为：对矩阵做一次sum(P_gen)得到的是1×3行向量而不是标量，直接扔进目标函数会报 nonscalar 错误。

T = units.T; % 24 N_gen = numel(units.gen_cost); % 3 台火电 P_gen = sdpvar(T, N_gen, 'full'); % 火电出力 MW，列对应机组 P_ele = sdpvar(T, 1); % 电解槽耗电功率 MW P_fc = sdpvar(T, 1); % 燃料电池发电功率 MW H2 = sdpvar(T, 1); % 储氢量（等效电量）MWh P_curt = sdpvar(T, 1); % 弃风功率 MW on_ele = binvar(T, 1); % 电解槽开停机状态 0/1 u_ele = binvar(T, 1); % 启动动作 0/1 E_seg = sdpvar(1, 3); % 阶梯碳排放的三段量，见第 3 章

配套的基础约束也一并给出，机组上下限、弃风上限、氢罐容量这三组是每个时段都要满足的：

cons = []; cons = [cons, repmat(units.gen_pmin, T, 1) <= P_gen <= repmat(units.gen_pmax, T, 1)]; cons = [cons, 0 <= P_curt <= P_wind]; % 弃风不能超过可用风资源 cons = [cons, units.h2_min <= H2 <= units.h2_tank_max]; % 留 20% 安全余量

repmat那一行是把1×3的上下限扩成T×3，YALMIP 支持这种逐元素双边不等式，写起来比循环干净。氢罐下限设成容量的 20% 而不是 0，是为了避免模型为了省钱把罐子抽干、第二天开局无氢可用，这个细节在第 4 章讲终端约束时还会再提一次。变量都定成连续量还是整数，直接决定求解时间：只有启停和分段选择用二进制，功率和储量一律连续，24 点的规模下 CPLEX 通常几秒内出解。

3. 阶梯碳交易的建模：把 if-else 变成 MILP 的分段线性化

3.1 为什么不能直接写 if-else 和 max

很多人第一版会写成if emission <= 200, cost = 80*emission; ...这种形式，跑起来发现要么报错，要么结果完全不对。原因是emission是sdpvar，它的值在建模阶段根本不存在，MATLAB 的if拿到一个未知对象只能按对象真假走固定分支，等于把你的分段逻辑整个跳过了。另一种写法是用min、max拼成本函数，比如80*min(E,200) + 120*max(0,E-200)，读起来漂亮，但max是分段线性的非光滑算子，YALMIP 会把它解释成非线性算子，CPLEX 的 MILP 分支根本不接受，求解器直接告诉你模型类型不匹配。

正确做法是在建模阶段就把分段结构展开成「若干非负变量 + 线性等式 + 上下界」，让整个模型保持纯 MILP。这也是为什么这份代码里碳排放量不是一个标量变量，而是一个长度为 3 的分段向量。碳排总量本身是线性的，由各机组出力和排放强度相乘求和得到，加一个固定排放基数（比如供热分摊部分）便于对齐实际园区的口径：

% 碳排总量：火电出力 × 排放强度，逐机组逐时段求和 total_emission = sum(P_gen * units.co2_int', 'all') + units.base_emission;

注意'all'这个参数，不加的话sum会对矩阵按列求和返回行向量，total_emission就不是标量，后面所有跟它相关的等式约束都会维度不匹配。这是我在别人代码里见过最频繁的低级错误之一，报错信息还特别不友好，只说约束维度不一致，不说到底哪一行。养成习惯：凡是求和得到标量的地方都显式写'all'，或者套两层sum(sum(...))。

3.2 三段填充法：递增碳价下的最小写法

代码包默认用一套三段递增碳价：前 200 吨每吨 80 元，200 到 500 吨每吨 120 元，超过 500 吨每吨 160 元。实现上引入三个非负变量代表落在各档里的排放量，加起来等于总排放，各档再卡上自己的容量上限，成本就是三个量分别乘单价再求和：

% 阶梯碳交易：三档递增价格的分段线性化 co2_seg = [200, 500]; % 档位分界点，单位 tCO2 co2_price = [80, 120, 160]; % 各档单价，元/tCO2 E_seg = sdpvar(1, 3); % E_seg(1) 落第一档，以此类推 cons = [cons, sum(E_seg) == total_emission]; % 三段之和 = 总排放 cons = [cons, 0 <= E_seg(1) <= co2_seg(1)]; % 第一档容量 200 t cons = [cons, 0 <= E_seg(2) <= co2_seg(2) - co2_seg(1)]; % 第二档容量 300 t cons = [cons, E_seg(3) >= 0]; % 第三档无上限 carbon_cost = co2_price(1)*E_seg(1) + co2_price(2)*E_seg(2) + co2_price(3)*E_seg(3);

这段代码里没有二进制变量，靠的是两个性质共同成立：单价严格递增，且目标函数是最小化。求解器为了让成本最小，会本能地先把便宜的第一档填满，再填第二档，绝不会出现「第一档没满却先用第三档」的浪费解，所以线性松弛本身就等于整数最优，不需要额外加 0-1 变量。一旦政策设计成第二档比第一档便宜（某些减免区间），或者第三档有封顶返还，这个性质立刻失效，模型会钻空子把排放全塞进便宜档，这时候就必须上二进制变量强制顺序填充。判断标准很朴素：把各档单价排成一列，看是否单调递增，是就用填充法，不是就加整数变量。

3.3 SOS2、二进制与动态 big-M：三种写法的取舍

除了填充法，工程上还有两种常见写法。一种是特殊有序集 SOS2，用一组权重 λ 表示当前解落在哪两个相邻断点之间，CPLEX 原生支持这种结构，分支效率比二进制高，缺点是 YALMIP 里要手工把 λ 和断点绑定，写错索引就退化成普通连续变量，约束形同虚设。另一种是最直白的二进制选择，给每一档配一个 0-1 变量表示「是否进入该档」，再用 big-M 把前后档位串起来：

% 二进制顺序填充：保证低档填满才允许用高档 z = binvar(3, 1); cons = [cons, sum(z) == 1]; cons = [cons, E_seg(1) >= co2_seg(1) * z(2)]; % 用第二档则第一档必满 cons = [cons, E_seg(2) >= (co2_seg(2)-co2_seg(1)) * z(3)]; % 用第三档则第二档必满 cons = [cons, E_seg(1) <= co2_seg(1)]; cons = [cons, E_seg(2) <= (co2_seg(2)-co2_seg(1))];

这里z(2)=1表示已经进入第二档，于是强制第一档的下界等于它的容量。真正的坑在 big-M 的取值：如果模型里还有别的逻辑用 M 表示「足够大的数」，很多人顺手写1e6，CPLEX 的对偶边界立刻变得极端病态，分支定界在节点上反复震荡，跑十几分钟给不出可行解。我的做法是让 M 由上界推出来，而不是拍脑袋：

% 动态 M：一天里火电可能产生的最大碳排，作为任何分段逻辑的紧上界 M_big = T * sum(units.gen_pmax .* units.co2_int) + units.base_emission;

这个值通常只有几千吨量级，跟真实解的差距在两三倍以内，数值条件好得多。三种写法的对比如下，选哪个主要看碳价是否单调、模型规模有多大：

4. 电制氢链路：电解槽效率、储氢罐状态方程与启停约束

4.1 电解槽的功率-产氢关系：二次拟合好看，进 MILP 就翻车

电解槽的实际效率不是常数，低负载时效率明显下滑，中高负载区间反而最稳。文献里常见的做法是拟合一条二次曲线，比如产氢量关于功率写成0.45*P - 2e-4*P^2这种形式。这条曲线本身没错，问题在于它带着平方项，一旦写进目标函数或约束，整个模型就从 MILP 变成 MIQCP，CPLEX 要动用二阶锥或者空间分支定界，24 点三机组的小模型还能撑住，扩到一周或者加进更多设备，求解时间就不是翻倍而是翻十倍。

代码包的处理方式是离散化：取负载率的三个断点 0、0.5、1.0，对应效率 0.60、0.72、0.68，用一组权重 λ 和一组段选择变量 z 把「当前工作在哪个区间」表达清楚。这种写法本质上是把曲线掰成两段折线，误差控制在 2 个百分点以内，但模型保持纯线性：

% 电解槽效率分段线性化：断点 (0, 0.5, 1.0)，对应效率 (0.60, 0.72, 0.68) bp = [0, 0.5, 1.0]; eta = [0.60, 0.72, 0.68]; lam = sdpvar(T, 3); % 每个时段一组断点权重 z = binvar(T, 2); % 两个线性段的选择变量 for t = 1:T cons = [cons, sum(lam(t,:)) == 1, sum(z(t,:)) == 1]; cons = [cons, lam(t,1) <= z(t,1)]; cons = [cons, lam(t,2) <= z(t,1) + z(t,2)]; % 相邻断点才允许非零 cons = [cons, lam(t,3) <= z(t,2)]; cons = [cons, P_ele(t) == units.ele_pmax * (bp * lam(t,:)')]; cons = [cons, H2_prod(t) == units.ele_pmax * ((bp .* eta) * lam(t,:)')]; end

逻辑说明：lam表示当前工作点落在两个相邻断点之间的权重组合，z负责把非相邻的权重掐成 0，这就保证了「同一时刻只在一段折线上」。P_ele(t)由断点值加权得到，H2_prod(t)用断点上的产氢量加权，两者用同一组 λ，所以功率和产氢严格一一对应。参数上ele_pmax是 40 MW，意味着满负载时段产氢约 28.8 MWh 等效电量；中途效率峰值在 0.72，模型会自然倾向于让电解槽跑在半载以上，这也解释了为什么结果图里电解槽很少出现在 20% 负载以下的区间。这套写法比二次曲线的精度差一点，但换来了可预测的求解时间，工业调度里我觉得这笔交易划算。

4.2 储氢罐状态方程：索引偏移是负储氢的元凶

储氢罐是整条链路里唯一带时间耦合的环节，也是 bug 高发区。它的物理含义很直白：这一时刻罐里的氢 = 上一时刻的余量 + 本时段电解槽制氢 − 本时段燃料电池耗氢。用等效电量表达之后，耗氢那项要除以发电效率，因为燃料电池发 1 MWh 电需要消耗 1/0.58 ≈ 1.72 MWh 的氢。这个「乘还是除」的方向如果搞反，模型会认为燃料电池越发电罐里氢越多，跑出来的结果看着收敛，实际物理上完全不成立。

% 储氢罐动态：初值 + 制氢（乘效率） - 耗氢（除效率） cons = [cons, H2(1) == units.h2_init ... + units.ele_eff * P_ele(1) - P_fc(1) / units.fuel_cell_eff]; for t = 2:T cons = [cons, H2(t) == H2(t-1) ... + units.ele_eff * P_ele(t) - P_fc(t) / units.fuel_cell_eff]; end % 日终闭合：末时段储量不低于日初，防止吃空库存 cons = [cons, H2(T) >= units.h2_init];

这里我把每个时段的功率都按「本时段」写，而不是原文里那种「t 时刻储量由 t-1 时刻功率决定」的滞后写法。两种离散方式都能用，但滞后写法有个致命副作用：最后一个时段的功率不会影响任何储量，模型会把最便宜的时段全堆在第 24 点，图上看起来末时段制氢功率爆表、罐子却没变化。统一成本时段的写法，再补一条日终闭合约束，问题就消失了。h2_init一般取容量的 30% 到 40%，再配合第 2 章那条 20% 下限，罐子全天都在 20% 到 100% 之间游走，不会出现负储量。

4.3 启停次数与最小出力：0.3 这个系数从哪来

电解槽是贵设备，不能让它每小时在 0 和满功率之间来回横跳。模型里用两个二进制变量表达这件事：on_ele表示当前时段是否开机，u_ele表示当前时段是否发生了「从停到开」的动作。最小出力约束靠on_ele撑起来，启动次数约束靠u_ele累加：

% 电解槽开停逻辑与最小出力 cons = [cons, P_ele <= units.ele_pmax * on_ele]; % 停机时功率必须为 0 cons = [cons, P_ele >= units.ele_pmin * on_ele]; % 开机则不低于 12 MW cons = [cons, u_ele(2:T) >= on_ele(2:T) - on_ele(1:T-1)]; % 状态跃迁即一次启动 cons = [cons, u_ele(1) >= on_ele(1) - units.on_init]; % 首时段也要算启动 cons = [cons, sum(u_ele) <= units.max_start]; % 全天不超过 3 次

ele_pmin取 40 MW 的 0.3 倍，也就是 12 MW，这个系数不是拍脑袋来的：低于 30% 负载时电解效率会明显掉到 0.5 以下，单位氢的耗电量飙升，与其这样不如让设备停着。max_start默认 3，是设备厂家给的日均启停次数建议值。u_ele(1)那行很容易漏，漏了之后模型会把「凌晨先开一次机」当成免费的，启动成本少算一次。启动成本在目标函数里按次计价，通常设成 200 到 500 元一次，具体看设备容量，数值不需要很准，它的作用是把频繁启停的收益压下去。

5. 求解配置与踩坑排查：从 MIP gap 到负储氢的五条血泪记录

5.1 目标函数组装与 CPLEX 选项：分项成本必须能拆开看

目标函数一共四项：燃料成本、运维成本、弃风惩罚、阶梯碳成本，再加上电解槽启动成本这一小块。写的时候逐项算完再相加，不要图省事写成一个巨大的表达式，否则调参时根本看不出是哪一项在主导结果。

% 分项成本：每一项单独成变量，方便事后拆解 fuel_cost = sum(sum(P_gen .* repmat(units.gen_cost, T, 1))); om_cost = units.om_coef * sum(P_gen, 'all') + units.h2_om * sum(P_ele); curt_cost = units.curtail_penalty * sum(P_curt); start_cost = units.start_cost * sum(u_ele); total_cost = fuel_cost + om_cost + curt_cost + start_cost + carbon_cost; % 求解器选项：MILP 三条常用设置 ops = sdpsettings('solver', 'cplex', 'verbose', 1); ops.cplex.timelimit = 600; ops.cplex.mip.tolerances.mipgap = 0.01; ops.cplex.mip.tolerances.absmipgap = 50; sol = optimize(cons, total_cost, ops); diag.exitflag = sol.problem; diag.fuel = value(fuel_cost); diag.carbon = value(carbon_cost);

参数说明：mipgap = 0.01表示允许 1% 的相对最优间隙，调度问题里这个精度完全够用，实际总成本几百万元的时候误差几千元，低于风电预测误差带来的不确定性；absmipgap = 50是绝对间隙，作用是当目标函数量级很小的时候不让求解器为了 1% 死磕；timelimit必须设，没有时间上限的 MILP 在极端数据下可能跑一整夜。verbose = 1会打印分支节点数、当前最好解和下界，调试的时候看这两个值的差距缩得快不快，比盯着进度条有用得多。

5.2 五条踩坑记录：现象、原因、解决

现象一：结果里储氢量出现负值。原因几乎总是效率折算方向写反，或者状态方程的时序错位导致末时段没闭合。解决方式是先检查H2(t)是不是被h2_min这类下界约束真正包住了，再核对制氢乘效率、耗氢除效率的方向；如果模型仍然给出负值，把H2的下界临时收紧到 0 跑一遍，对比两次结果就能定位是约束漏写还是求解精度问题。

现象二：模型显示 infeasible，但怎么看参数都合理。最常见的原因是储氢容量和制氢功率不匹配：40 MW 电解槽满发一小时产 28.8 MWh 等效氢，如果罐子只给了 30 MWh 容量，早上四个小时连发就把罐子撑爆，模型只能报不可行。解决办法是把h2_tank_max和ele_pmax联动设计，经验上罐容至少能装下电解槽满发 6 到 8 小时的产量，或者干脆给弃风留出更大空间。

现象三：求解时间从几秒变成十几分钟。十有八九是某处 big-M 取了1e6这类夸张值，或者碳排总量用了一个很宽的人工上界。把 M 换成由设备容量推出来的紧上界之后，同一个模型经常在一分钟内收敛。这条血泪经验我踩过不止一次，后来养成习惯：模型里每出现一个 M，都要在注释里写清它是怎么算出来的。

现象四：碳成本高得离谱，比燃料成本还大一个数量级。单位混用。有的文献碳价按元/kg 写，有的按元/t 写，两者差 1000 倍；还有些参数表里排放强度给的是 kgCO2/MWh 而代码按 t 处理。排查方法很简单，把value(total_emission)打出来，看量级是不是在几百到几千吨，如果只有零点几，说明单位错了一千倍。

现象五：电解槽全天频繁启停，启动次数约束好像没生效。检查u_ele的定义，很多人只写了u_ele(t) >= on_ele(t) - on_ele(t-1)，却忘了首时段，也忘了把u_ele放进目标函数。如果启动成本系数是 0，模型会认为启停免费，约束虽然限制了次数上限但 3 次之内随便跳，图上就显得很乱。补上首时段约束、把启动成本调到几百元每次，曲线立刻变干净。

6. 二次开发：从单日 24 点扩到多日滚动，我一般这么改

单日模型最大的局限是储氢罐必须在当天闭合，跨天的能量搬移被硬生生切断。要扩成一周或者一个月的滚动调度，改动集中在外层：把build_model.m包进一个按天推进的循环，每天调用一次，把前一天的末态作为当天的初值传进去，同时把日终约束从「不低于日初」改成「不低于某个保底值」，只在最后一天强制闭合。

% 多日滚动：末态承接，只有最后一天强制闭合 h2_state = units.h2_init; for d = 1:N_day units.h2_init = h2_state; units.is_last = (d == N_day); % 最后一天才要求罐子回到初始水平 [sol, diag] = build_model(units); h2_state = value(sol.H2(end)); rec.cost(d) = diag.total; end

改成滚动之后，验证方法也要跟着变。我一般做三件事：一是拿固定碳价版本跑同一组数据做基准，看阶梯碳价下的总成本是不是落在 8% 到 15% 的下降区间，超出这个范围多半是参数写错；二是检查能量守恒残差，把每个时段的功率平衡左右两边相减取绝对值，最大值应该小于 1e-6，否则就是有约束没写全；三是把碳档位曲线画出来，看火电出力在跨档时刻有没有出现预期的拐点，没有拐点说明碳价传导到调度决策的链路断了。

再往深一层，可以加碳捕集设备把排放压回第一档，或者把风电预测换成多场景鲁棒版本。鲁棒版不用重写模型，只要把P_wind从常量改成场景集合，再把功率平衡约束对每个场景各写一份，用binvar选最坏场景即可，代价是二进制变量翻倍，求解时间也会跟着上去。改造氢能部分时，我每次都会先把储氢状态方程单独拿出来，用固定功率输入跑一遍开环仿真，确认罐子曲线单调合理了再放进优化模型，不然一旦出现负储量，你根本分不清是物理模型错了还是求解器给了个奇怪的解。这套排查顺序希望帮到你。

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