★ FEATURED ARTICLE 简介面向具备一定MATLAB基础、从事综合能源系统优化的研究人员与工程技术人员这份PDF围绕阶梯型碳交易与电制氢的联合调度建模展开重点解决碳成本分段跃迁、电—氢—电多向转换及24小时时序耦合下的优化求解问题。文档以传统机组、风电、电解槽、储氢罐和燃料电池五个单元构成系统架构拆解阶梯碳交易成本计算、储氢罐状态方程、弃风惩罚与目标函数设计并借助YALMIP调用CPLEX完成混合整数线性规划求解。压缩包仅含1个PDF文件大小约129KB篇幅紧凑但覆盖模型构建、代码实现与结果可视化。已有162人学习适合想理解分段线性化、big-M处理、时间索引对齐和求解器数值稳定性等细节的读者通过算例还能观察风电大发时段优先制氢、碳排放接近阈值时火电出力拐点等现象便于将方法迁移到自身调度模型。1. 阶梯碳价与电制氢为什么一个成本拐点能把线性规划逼成混合整数规划把阶梯碳交易塞进综合能源系统调度模型直觉会说不过是给成本函数加了个拐点。真跑起来才发现这个拐点会把原本用单纯形法几十毫秒解完的线性规划变成必须靠分支定界才能收敛的混合整数规划再叠上电制氢这条电—氢—电的双向链路MATLAB 里一行写得下的效率折算会在第 24 个时段把储氢量算成负数。用 YALMIP 建模、CPLEX 求解这套组合做过 24 时段调度的开发者多半都卡在同样的三个地方阶梯段怎么线性化、时间索引怎么对齐、big-M 取多大。这三个问题的答案直接决定模型是半分钟出结果还是跑一小时 gap 还挂在 15%。2. 阶梯碳交易为什么一个 if-else 会把整个模型毁掉2.1 阶梯计价的经济学直觉本质是一个凸分段线性函数固定碳税和阶梯碳价的区别不在“单价高低”而在“边际成本是不是常数”。固定碳税下每多排一吨的代价永远是同一个数调度模型只会做一个权衡清洁能源的边际成本低于碳税就用清洁能源。阶梯碳价则让边际成本随排放量跃迁前 200 吨可能 80 元/吨跨过阈值瞬间变成 120 元/吨再往上 160 元/吨。这意味着系统在排放接近阈值时会主动改变机组组合去“躲”那一跳火电出力曲线因此出现拐点而不是平滑下压。从建模角度看这个成本函数写出来是C(E) 80*min(E,200) 120*max(0,E-200) ...斜率单调不减属于凸函数。凸性非常关键在最小化问题里凸性等价于“优化器会自动从低价档往高价档填”这个性质后面会直接换掉一堆二进制变量。反过来如果阶梯设计成第二档比第一档便宜补贴退坡或者负碳价场景凸性就没了模型必须老老实实引入整数变量求解规模立刻上一个台阶。还有一个容易忽略的点total_emission不是常数它是决策变量P_gt的线性组合。所以碳成本不能像后处理那样用一个函数算完再塞进目标它必须在建模阶段就展开成决策变量的表达式否则 YALMIP 拿到的只是一个没有导数的黑箱。2.2 辅助变量写法用上界约束替代分段判断最常见的实现思路是给每一档碳价配一个独立的排放量变量各档之间靠容量上界互相约束% 碳排放分段变量E(1) 落在第一档E(2) 第二档E(3) 第三档 E sdpvar(1,3); cap [200, 200, 400]; % 各档可容纳的排放增量吨 price [ 80, 120, 160]; % 各档单价元/吨 % total_emission 是火电出力的线性组合必须先算出来 total_emission sum(e_factor * P_gt) * dt; cons [cons, sum(E) total_emission]; % 分段量之和等于总排放 cons [cons, 0 E(1) cap(1)]; cons [cons, 0 E(2) cap(2)]; cons [cons, 0 E(3) cap(3)]; carbon_cost price * E; % 直接线性加权这段代码的关键不在那三行上界而在“没有写任何东西”的地方——没有二进制变量、没有 big-M、没有 SOS2。原因是price严格递增而目标是最小化优化器在把sum(E)顶到total_emission的同时会本能地把E(1)填满再去动E(2)。填不满低档就跳到高档等于自己给自己加钱最优解不会这么干。这就是凸性带来的红利分段逻辑被目标函数的方向“自动”实现了。参数上有一处极易写错cap存的是每档的增量容量不是累计阈值。如果按cap [200, 400, 600]来写等于第二档能装 400 吨总容量被放大到 1200 吨模型会算出一个明显偏低的碳成本而且不会报错——这是最难查的一类错误。2.3 效率曲线用 SOS2连续分段函数的另一种打开方式阶梯碳价的分段点落在“成本”上电解槽的分段点落在“效率”上。后者是连续折线不是跳变处理手法不同。YALMIP 里比较干净的写法是 SOS2特殊有序集让一组权重变量最多只有两个相邻分量非零% 电解槽分段效率负载率-效率折线 bp [0, 0.5, 1.0]; % 负载率断点归一化 eta [0.60, 0.72, 0.68]; % 对应效率实测值 lam sdpvar(3,1); % 凸组合权重 P_emax 30; % 电解槽额定功率MW cons [cons, lam 0, sum(lam) 1]; cons [cons, P_elec P_emax * (bp * lam)]; % 功率落在断点凸包内 cons [cons, eta_eff eta * lam]; % 效率同步插值 cons [cons, sos2(lam)]; % 至多两个相邻权重非零 H2_prod eta_eff * P_elec * dt; % 时段产氢量kg 或 Nm3sos2是 CPLEX 原生支持的约束类型求解器在分支时只需要沿折线相邻段二选一分支树比“3 个二进制变量 每个配一组 big-M”小一个量级。代价是它只能表达连续折线表达不了 2.2 节那种真正的跳变——阶梯碳价如果非要上 SOS2得先把它改写成“斜率累加”的形式反而绕远。这里还有一个物理细节值得强调eta在负载率 0.5 处最高、1.0 处反而下降说明电解槽并不是越满负荷越划算。这条曲线一旦进入模型调度结果会出现一个反直觉现象——风电大发时段电解槽不一定顶到额定功率而是停在 0.5 到 0.8 负载率之间。如果偷懒用常数效率这个拐点直接消失模型给出的“最优”调度在现场根本复现不了。2.4 三种分段写法的取舍写法是否要求凸性新增变量量级求解表现典型场景辅助变量 上界必须凸价格递增每段 1 个连续变量几乎无额外开销阶梯碳价、阶梯电价SOS2 凸组合不要求但只能连续每断点 1 个连续变量分支浅收敛快效率曲线、成本曲线拟合二进制 big-M无要求每段 1 个二进制 1 个连续分支树最大对 M 敏感非凸阶梯、启停逻辑选型顺序很清楚先判断函数是不是凸的凸就用辅助变量不是凸但连续用 SOS2既不凸又是跳变才动二进制。很多模型跑不动不是 CPLEX 不行是一开始就把一个凸问题写成了整数问题。3. 电制氢电—氢—电的双向折算为什么会算出负储能3.1 效率的方向性该乘的时候乘该除的时候必须除氢能系统最容易被写错的一行是效率到底乘在哪一侧。电解槽把电能转成氢效率打折打在“产物”上所以产氢量等于效率乘功率燃料电池把氢转回电能效率打折同样打在“产物”电能上于是反推氢气消耗量时必须除以效率而不是再乘一次。物理上这很直白——70% 的效率意味着要烧掉 1/0.7 份氢才能发出 1 份电——但写成代码时两个方向都用乘号是高频翻车点而且模型不会报错只会给出一个“凭空多出氢气”的可行解。% 储氢罐动态注意两个方向效率的处理方式不同 for t 2:T cons [cons, H(t) H(t-1) ... eta_e * P_elec(t-1) * dt ... % 电解产氢乘效率 - P_fc(t-1) / eta_f * dt]; % 燃料电池耗氢除效率 endeta_e取 0.72、eta_f取 0.58 是常见的工程区间两个数字一乘只剩 0.42也就是说“电—氢—电”走一圈要损失近六成能量。这个损耗率直接决定了氢能在系统里的定位它不是一个高效的储能装置而是一个跨时段、跨品类的调节手段用来把凌晨的低价风电搬到晚高峰同时顺带躲开阶梯碳价的第二档。理解这一点才不会对“为什么模型只在特定时段开电解槽”感到困惑。3.2 索引偏移为什么最后一个时段总是状态不闭合for t 2:T这条循环只产生了T-1个方程而H有T个分量。数学上这意味着H(1)和H(T)都处在自由状态模型有两处“漏洞”开头可以凭空造出一罐氢结尾可以把氢直接丢掉。两种行为都能降低总成本所以最优解一定会往这个方向走。% 日循环闭合防止跨日“造氢”或“丢氢” cons [cons, H(1) H_init]; % 初始储量 cons [cons, H(T) H_init]; % 日末回到初值 cons [cons, 0.2*H_max H 0.9*H_max]; % 安全区间不是 0 到满罐这条约束看起来只是一行等式但它决定了结果的可用性。漏写时的典型症状是储氢曲线末端一路下探到 0总成本比补齐后低 5% 到 8%模型“显得”特别经济实际调度时第二天开局就没氢可用。另一个细节是储氢上下界不能取 0 和满罐罐压过低或过高都会触发保护取 0.2 到 0.9 的区间既符合运行规程也顺手避开了边界上的数值不稳定。3.3 电解槽的物理约束下限、启停与 big-M 绑定电解槽不是想开多大就开多大。低负载率下电解效率断崖式下跌所以工程上有一条“要么停机要么不低于额定 30%”的硬约束。这条约束本身是非凸的——P 0或者P 0.3*P_max中间的区间不允许——只能靠二进制变量表达u binvar(T,1); % 电解槽启停状态 start binvar(T,1); % 启动动作标记 P_emax 30; for t 1:T cons [cons, P_elec(t) P_emax * u(t)]; % 上限绑定 cons [cons, P_elec(t) 0.3 * P_emax * u(t)]; % 运行下限 end % 启动动作仅在 u 由 0 变 1 时为 1 cons [cons, start(2:T) u(2:T) - u(1:T-1)]; cons [cons, sum(start) 3]; % 24 小时内最多启停 3 次这里的P_emax同时扮演了两个角色功率上限以及 big-M 的取值。这是有意为之。big-M 最理想的取值就是“物理上不可能超过的那个界”取大了会让线性松弛过松、分支定界剪不动枝取小了会直接切掉可行域。用设备铭牌值当作 M等于让模型自己给出最紧的松弛。启停次数约束的写法也值得看一眼start(t) u(t) - u(t-1)是必要非充分条件优化器在最小化启停成本时会把start压到最小于是不等式自动收紧成等式不需要额外的逻辑约束。这是一个很典型的 MILP 建模技巧——用目标函数的压力替代显式的逻辑等式。3.4 储氢罐的容量换算压力不是质量别直接当变量用储氢罐的物理上限是压力而调度模型需要的是质量或标准状态下体积两者之间隔着一个状态方程。常见做法是用理想气体方程把压力上限折算成等效储氢量上限再对温度项做线性化处理% 压力上限折算为等效储氢质量上限等温假设下的简化 % m_max P_max * V * M / (R * T_amb) H_max P_max * V_tank * M_h2 / (R_gas * T_amb); % 低压段对压缩因子做一阶修正避免大范围泰勒展开误差 H_eff H_max * (1 - beta * (P_max - P_op) / P_max);这一步能把求解时间压下来三到四成因为原本需要迭代校核的非线性关系被换成了一个常数上限。代价是精度在压力偏离工作点较远时等效质量会有几个百分点的误差。如果项目对储氢量的绝对值敏感比如要做容量配置这个简化就得慎用宁可保留分段线性也不要一口气线性化。4. 目标函数与求解器MILP 的规模是怎么被撑大的4.1 四块成本的量纲对齐目标函数看着简单实际最容易出问题的是量纲。燃料成本按元/MWh 乘功率乘时长弃风惩罚按元/MWh 乘弃风电量碳成本按元/吨乘吨数启停成本按元/次乘次数——四块必须全部归到“元”否则优化器会优先牺牲量纲最小的那一块结果看起来“成本很低”实际上是把排放或弃风甩出去了。% 目标函数四块成本量纲统一为元 fuel_cost sum(a_gen .* P_gt b_gen .* u_gt) * dt; % 燃料 空载 curtail_pen c_curt * sum(P_curt) * dt; % 弃风惩罚 carbon_cost price * E; % 阶梯碳成本 start_cost c_start * sum(start); % 电解槽启停 total_cost fuel_cost curtail_pen carbon_cost start_cost;有一个顺序问题必须注意carbon_cost依赖的E是由total_emission推出来的而total_emission又是全天P_gt的求和所以碳成本的计算必须放在时段循环外面。放进循环里逐时段算等于把阶梯阈值按小时切碎了200 吨的分档变成每小时 8.3 吨阶梯机制彻底失效算出来的曲线会平滑得像固定碳税。c_curt这个弃风惩罚系数是另一个隐形杠杆。它设得太低模型宁愿弃风也不肯启电解槽设得太高又会逼着电解槽在低效区间硬撑。工程上一般把它设在“风电边际价值 少量惩罚”的量级让弃风和制氢在成本上真正可比而不是一边倒。4.2 求解器参数1% 的 gap 为什么是工程上的合理选择CPLEX 解 MILP 的本质是分支定界不断用线性松弛给出下界用可行整数解给出上界两者夹逼到足够近就停。mipgap控制的就是这个“足够近”的门槛。ops sdpsettings(solver,cplex,verbose,2); ops.cplex.timelimit 3600; % 硬性时限防止卡死 ops.cplex.mip.tolerances.mipgap 0.01; % 上下界相对差 1% ops.cplex.mip.strategy.search 1; % 传统分支定界规模小时更稳1% 听起来很粗但要看基数。如果总成本在百万元量级1% 就是一万块而调度指令本身的分辨率是 MW 级、时段是小时级再往下抠小数点后三位对现场没有意义。反倒是不设timelimit更危险——碰到退化严重或对称性强的算例分支定界可能跑几个小时还在原地打转而你要的只是一条能用的出力曲线。verbose设成 2 是为了看收敛过程。真正有用的信息不是最终的 gap 数值而是“下界提升的速度”如果前 100 秒下界就爬到接近最优说明松弛够紧如果跑了十分钟下界几乎没动问题多半出在 big-M 太松或者约束之间线性相关改参数没用得回去改模型。4.3 结果自检三条必须跑的校验优化器返回Problem solved不等于结果可信。求解完之后至少要跑三条反算把“数学上可行”翻译成“物理上说得通”。% 校验一功率平衡残差数值误差应在 1e-5 以内 res value(P_gt P_wind_used P_fc - P_load - P_elec); assert(max(abs(res)) 1e-5, 功率平衡不闭合); % 校验二碳排放反算确认与分摊变量一致 E_check value(sum(e_factor * P_gt) * dt); assert(abs(E_check - sum(value(E))) 1e-4, 碳分摊量与实际排放不符); % 校验三储氢日循环闭合 assert(abs(value(H(1)) - value(H(T))) 1e-4, 储氢未回到初值);第三条尤其重要。前面提过漏写日循环约束会让成本虚低但还有一种情况约束写了、求解也收敛了可因为H_init设得比H的下界还低模型在第一个时段就“违规”造氢。这种解不会触发 infeasible只会让校验二、校验三同时飘红所以三条断言最好一起跑。4.4 两个现象背后的机制跑通之后有两个结果值得单独解释因为它们恰好验证了前面所有建模选择的正确性。第一个是凌晨风电高峰时段电解槽优先启动而不是弃风。这不是因为制氢“更环保”而是成本账算下来的结果弃风要付c_curt的惩罚制氢要付启停成本和运维成本但制氢产出的氢可以在晚高峰替代火电、顺带少排碳。当弃风惩罚高于制氢的边际成本时电解槽自然被顶起来。第二个是碳排放接近阈值时火电出力出现拐点。这背后是边际碳价从 80 跳到 120 的那一跳——在阈值左侧多排一吨的代价是 80跨过去之后新增排放的代价变成 120等效于火电的边际成本突然上涨了 40 元/吨乘以排放因子。系统在阈值附近“刹车”正是阶梯机制在数学上的必然表现而不是数值噪声。5. 常见问题排查数值、索引与约束松紧的五次翻车5.1 求解器报 infeasible但约束逐条看都满足现象CPLEX 直接返回 infeasible把碳交易那几行注释掉就立刻可行说明问题出在分段部分。原因有两个高发点。一是cap总和小于实际排放分段变量装不下总量sum(E) total_emission自然无解二是cap被误写成累计阈值[200, 400, 600]第二档容量算成了 400虽然看起来更宽松但和后续的单价加权一起把成本结构扭曲模型在别的约束上被迫违约。解决求解前先打印sum(cap)和“火电全开时的最大排放”比一比再用 YALMIP 的诊断输出定位是哪一条约束被松弛掉。养成习惯把每个分段变量的上界单独assert一次比事后调试省时间。5.2 big-M 取大了跑不动取小了直接切掉可行域现象同一个模型M 取 1e4 时三十秒收敛取 1e6 时跑一小时 gap 还挂在 15%。原因big-M 的本质是线性松弛的“松紧度”。M 越大松弛后的可行域越宽分支定界拿到的下界越松剪枝越难同时系数矩阵的条件数变差CPLEX 在数值上还会额外做缩放进一步拖慢。解决不要拍一个大数用设备物理上限当 M。电解槽用额定功率储氢罐用容量上限火电用最大出力。这些值天然是“不可能超过”的界既不会切掉可行解又能给求解器最紧的松弛。% 动态 big-M全部来自设备铭牌不出现魔数 M_elec P_emax; % 电解槽功率上限 M_h2 H_max; % 储氢容量上限 M_gt P_gt_max; % 机组最大出力5.3 结果里储氢量跑成负数或者超过罐容现象value(H)里有小于 0 的分量或者某几个时段顶到了远大于H_max的值。原因效率方向写反是最常见的。燃料电池那一项写成- eta_f * P_fc等于把耗氢量缩小了将近一半罐子里的氢怎么用都用不完模型就会放心地超额放氢最后在闭合约束上崩掉。另一个原因是漏了H 0的下界只写了上界。解决先把效率量纲单独拿出来核对一遍物理单位再显式补上下界。调试时可以临时把eta_f设成 1.0如果负储能消失基本可以确认是效率方向的锅。5.4 碳成本算出负数或者高档变量先被填满现象carbon_cost是负值或者结果里E(3) 0但E(1)没满。原因一是price不单调比如第二档手滑写了 60优化器发现“多排更便宜”就会把排放往第二档猛塞二是给E加了不该有的等式约束强行指定了某档的取值破坏了自动分层的性质。解决建模前assert(all(diff(price) 0))把单调性写成硬检查检查E上是否只挂了上界和总量等式没有多余约束。5.5 参数微调一点点出力曲线完全变样现象碳价从 80 改到 82机组组合和电解槽启停序列整个换了一套看起来像模型不稳定。原因这不是 bug是阶梯阈值附近的退化解现象。当最优解恰好落在分档边界上存在多个成本几乎相同但结构完全不同的解MILP 求解器返回哪一个取决于分支顺序。再加上机组之间的对称性同容量机组可互换解的不唯一性会被放大。解决如果只是要一条可执行曲线接受它现实调度本来就有裕度如果要对比不同碳价场景加一个极小的二次正则项比如1e-4 * sum(P_gt.^2)把解“拉平”代价是模型从 MILP 变成 MIQP求解时间上升一个档次折中做法是先固定机组组合再对出力做二次优化。6. 设计取舍把政策曲线当成模型的一部分而不是一个可调参数阶梯碳交易的段数是个容易被低估的决策。三段够用五段更贴近真实政策但每多一段就多一组连续变量和一条上界约束凸性还在时开销很小一旦某一段设计成非凸就得配二进制分支树立刻膨胀。我一般的判断标准是先用三段跑通、验证机制正确再逐步加密每加一段看一次mipgap收敛曲线如果段数翻倍带来的是求解时间翻三倍而成本只降 0.3%就该停手。阶梯段数新增连续变量是否需二进制典型求解时间量级适用阶段3 段3否凸秒级机制验证、方案比选5 段5否凸十秒级政策敏感性分析非凸分段每段 1 连续 1 二进制是分钟到小时特殊政策设计验证阶梯机制有没有真的生效有个很朴素的办法把碳价参数扫一遍比如 60 到 200步长 10画总成本随碳价变化的曲线。如果阶梯确实起作用曲线会在阈值附近出现斜率突变如果画出来是一条直线说明分段逻辑被写成了固定碳税多半是cap写成了累计值或者碳成本被放进了时段循环。效率曲线也是同样的取舍逻辑。三点折线能抓住“低负载效率低、中负载最高、满负载略降”这个形状五点更准但边际收益递减用二次函数拟合看着优雅可一旦进了 MILP 就得再做一次分段线性化绕一圈回到原点。我倾向于三点折线加 SOS2精度够用求解器友好。从那以后我每次建这类模型都强制先做一件事把成本函数单独画出来确认它是凸的、确认分段点是增量不是累计、确认最后一个时段的储氢能闭合三条都过了再交给 CPLEX。这三步花不到十分钟但能省掉后面几个小时的“结果看着对、其实全是漏洞”。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
