最近一直在折腾一个综合能源系统的项目核心是把电转气Power-to-Gas, P2G和碳捕集系统Carbon Capture System, CCS硬塞进热电联产Combined Heat and Power, CHP机组的优化模型里然后用Matlab把整套程序撸出来。说实话这个方向不是那种能让你“跑个demo截两张图交差”的简单任务而是要把电力、热力、天然气、碳四套流动在同一个时间断面上做联合优化光是把物理过程翻译成数学语言的功夫就不小。这篇文章就是我这次项目的完整复盘从建模思路、数学推导、代码实现到调试踩坑全都有。如果你是做综合能源系统、低碳调度或者碳捕集与电转气协同优化方向的研究生或者工程师并且打算用Matlab复现一套可运行的仿真模型那这篇内容基本就是你需要的“抄作业”模板。先说清楚这套东西能解决什么问题。传统热电联产机组有一个著名的“以热定电”约束也就是冬季供暖需求高的时候机组必须保证一定的电出力下限这直接挤压了风电、光伏的上网空间弃风和碳排放两个问题一起冒出来。引入电转气之后盈余电力可以电解水制氢进一步甲烷化注入天然气网络引入碳捕集之后燃煤燃气机组排放的二氧化碳有机会被捕集下来而P2G甲烷化刚刚好需要二氧化碳作为原料。这样一来CHP、P2G、CCS三者形成了一个“电-热-气-碳”闭环既降低了机组的电出力下限又让捕集下来的碳有了去处。我这次用Matlab把这三套设备在一个统一的优化框架下建模目标函数是最小化系统运行总成本约束条件包括功率平衡、爬坡、储罐容量、碳捕集运行范围等最终求解结果让我觉得这条路是真实可行的。1. 项目从哪来为什么把P2G、CCS和CHP放在一起建模1.1 综合能源系统里的“硬骨头”综合能源系统Integrated Energy System, IES这几年已经是学术圈和工程圈的高频词了核心思想很简单电、热、气、冷这些能源形式别再各干各的而是通过耦合元件把它们串起来让能量在多能流之间梯级利用、互相补充。理想很丰满但真到了建模优化这一步就发现IES的复杂度是指数级上升的。电网有节点潮流和频率约束热网有传输延迟和温度梯度气网有压力和管存动态碳流还要考虑捕集、运输、利用全链条任何两个网络耦合在一起约束条件就会交叉感染。我在项目里选择以CHP机组作为核心枢纽原因是它在综合能源系统里太典型了。CHP同时输出电和热往上是燃气供应和电网连接往下是热力管网和电力负荷天然就是一个“多点耦合”的设备。但正是这种耦合让调度变得扭捏——冬天热负荷上来了机组电出力不得不顶着下限运行碰上凌晨风电大发的时候系统就只能被迫弃风。这是真实电网里的痛点也是我这篇文章里整套模型的核心动机。1.2 P2G和CCS为什么是一对“天作之合”单独看P2G和CCS其实各自都有一些“尴尬点”。P2G的全链条是“电解水制氢→氢甲烷化”甲烷化需要二氧化碳如果二氧化碳来自外部购买或者空气直捕成本高、碳源不稳定经济性就很差。CCS的尴尬则在于捕集下来的二氧化碳如果只是封存在地下监管成本和技术风险都不小而且纯封存没有任何经济回报只能靠碳价硬撑着。这两者凑到一起就变了CCS捕集的CO2直接就地供给P2G的甲烷化反应器省掉了碳运输环节同时P2G生产的人造天然气还能回注气网替代一部分天然气采购等于把“碳排放”变成了“碳资源”。这个耦合关系不是硬凑的它是化学上的天然匹配。甲烷化反应的反应式是CO2 4H2 → CH4 2H2O二氧化碳和氢气按比例消耗温度压力条件在常规工业范围内是可以实现的。所以从系统角度看P2G和CCS在能量流上还多了一层碳元素流建模的时候必须把每一单位甲烷产气对应的CO2消耗量和H2消耗量都精确算清楚。这也是我在这套模型里最注意的地方一个反应式系数写错整个碳平衡就崩了。1.3 为什么我用Matlab而不是Python这个问题的答案可能要得罪一部分Python党但我还是想说句公道话。单论生态Python的Pyomo、Gurobipy在优化建模上确实很能打但Matlab在“矩阵化约束建模快速原型验证”这个组合上依然有独特优势。Matlab的YALMIP工具箱是学术界做调度的老朋友了它允许你直接用符号变量描述决策变量和约束然后底层调用Cplex、Gurobi、甚至Matlab自带的intlinprog求解这对“快速验证一个优化想法”来说是极其顺手的。而且对于理工科背景的研究生来说Matlab的矩阵操作习惯和论文里的数学表达式天然契合很多公式几乎可以直接翻译成代码。我这次用的具体版本是R2023b不过核心代码并没有依赖任何新特性YALMIP和求解器配置部分从R2019a到R2026b都可以跑通所以版本适配基本不用太操心。下面的分享我一律按“能用就行”的原则来写不故意炫技保证你拿过去改改参数就能复现。2. 核心数学模型电-热-气-碳四流耦合框架怎么搭2.1 CHP机组的可行域与运行约束热电联产机组的建模有几个流派。一种是最精细的非线性模型把汽轮机抽汽工况的热电特性曲线直接写进去精度高但求解困难不适合放在小时级调度模型里。另一种是工程上常用的“线性可行域”模型把机组电出力和热出力的关系描述为一个凸多边形所有运行点都落在多边形范围内。后者虽然精度略低但对于一个以“系统级优化”为目标的项目来说完全够用而且线性约束能让模型保持为混合整数线性规划MILP求解稳定性好得多。以我项目里的抽凝式CHP机组为例电出力P和热出力H之间的关系可以写成一组线性不等式。简单来说电出力有个上下限但上限会随着热出力增加而下降因为抽汽多了做功少了下限也可能和热出力呈耦合关系。把这些边界写成约束就是P_min a * H ≤ P ≤ P_max - b * H其中a、b是根据机组特性曲线拟合出来的系数。同时还要加上爬坡约束也就是相邻两个时段电出力变化量的上限-P_ramp ≤ P(t1) - P(t) ≤ P_ramp如果模型里需要区分开停机状态那就再加一组整数变量u(t)配合开机成本、最小运行时间、最小停机时间约束。我这次并没有加上最小启停时间约束因为算例只需要看连续48小时的调度结果机组启停次数本身不频繁加不加对成本影响不大。2.2 P2G的能耗与物料平衡P2G系统内部其实包含两个核心反应电解水制氢和二氧化碳甲烷化。电解水部分消耗电能产物是氢气和氧气甲烷化部分消耗氢气和二氧化碳产物是甲烷和水。整个P2G系统的输入能量是“电能二氧化碳”输出能量是“天然气甲烷”同时还伴随一部分热损失。建模时需要用如下几个关键参数η_elec电解效率kW·h电对应多少Nm³氢气η_meth甲烷化效率氢气和CO2转化为CH4的转化率r_CO2每生产1单位甲烷所需CO2量化学反应计量比r_H2每生产1单位甲烷所需H2量那么某一时段t内P2G产气量V_P2G(t)与耗电量E_P2G(t)、CO2消耗量C_P2G(t)、H2消耗量的关系是E_P2G(t) V_P2G(t) / η_elec / η_meth C_P2G(t) r_CO2 * V_P2G(t) H2_cons(t) r_H2 * V_P2G(t)这里的H2可能还要考虑中间储氢罐因为电解制氢和甲烷化的运行不一定完全同步。如果不想引入过多的中间变量可以直接让P2G作为一个整体设备电解和甲烷化同步运行。我项目中为了简化就采用了这种同步方案H2只在系统内部传递不对外销售也不外购因此H2平衡约束自动满足不需要额外建模。但如果你后期想研究“氢储能”或“掺氢天然气”方向那就要把H2储罐的动态约束单独拿出来建。2.3 碳捕集系统的能耗与“碳流”出口CCS系统的本质是“把烟气里的CO2分离出来”技术路线有燃烧后捕集、燃烧前捕集、富氧燃烧等建模时最常用的是燃烧后化学吸收法。它的关键参数有两个捕集率η_ccs和单位捕集能耗β。单位捕集能耗通常由两部分组成再生塔需要的热耗通常来自CHP抽汽和吸收塔、压缩机需要的电耗。这就有意思了因为CCS不仅仅消耗电力还消耗热力所以在热电联产系统里CCS相当于一个“电热双耗”负荷这恰恰是CHP的优势所在——热负荷需求增大时机组可以多抽汽供热电出力下限也能被CCS的耗电“吃掉”一箭双雕。捕集下来的CO2去向无外乎三种封存、利用、驱油。在我的模型里CO2的主要去向是供给P2G甲烷化剩余部分按碳交易价格出售或者封存抽象成一个经济项。因此碳流约束是C_capture(t) C_P2G(t) C_sell(t)其中C_capture(t)是CCS捕集总量C_P2G(t)送入P2GC_sell(t)外售或封存。另外CCS装置自身还有一个运行范围约束捕集量不能超过吸收塔的处理上限也不能低于最小稳定负荷。这个约束在极端工况下挺容易踩到后面调试部分我会细说。2.4 目标函数与约束条件清单整套优化模型的目标函数是系统在一个调度周期通常24小时或48小时内的总运行成本最小化主要组成项包括购电成本从主网购入电力的费用P_grid(t) * price_grid(t)燃料成本天然气消耗费用CHP机组耗气量V_gas(t) * price_gas(t)碳交易成本净碳排放量总排放 - 捕集量* 碳价启停成本机组启动和停机产生的额外费用运维成本P2G、CCS设备的单位维护成本和运行量成正比完整的约束条件我用一个清单来整理方便你对号入座电力平衡CHP出力 P2G耗电 CCS耗电 风电光伏出力 主网购电 电负荷热力平衡CHP热出力 燃气锅炉如果有 热负荷 CCS再生热耗天然气平衡CHP耗气 天然气负荷 气网购气 P2G产气CHP运行约束可行域、爬坡、启停逻辑P2G运行约束转化效率、产气上下限、储罐容量CCS运行约束捕集量上下限、能耗线性关系、CO2储存容量碳流平衡捕集量 P2G耗碳 外售碳这里每一项都能在Matlab里翻译成一行或几行约束代码。我建议你别把约束写成一坨而是像下面这样分模块组织后续排查问题的时候会很舒服。3. Matlab实现从零搭建仿真框架的实操记录3.1 代码总体架构与变量定义整套代码我分成了五个模块case_data数据输入、build_model模型构建、solve_model求解调用、post_process结果处理、plot_results画图展示。这种分法的好处是想换一个算例或者修改某个参数时只需要动case_data模型和求解部分都不用大改。变量定义这一块我用的是YALMIP的sdpvar和binvar语法% 时间粒度1小时调度周期T24 T 24; % 连续变量 P_Chp sdpvar(1, T, full); % CHP电出力 H_Chp sdpvar(1, T, full); % CHP热出力 V_P2G sdpvar(1, T, full); % P2G产气量 C_capture sdpvar(1, T, full); % CCS捕集CO2量 P_grid sdpvar(1, T, full); % 主网购电 P_wind sdpvar(1, T, full); % 风电上网出力假设可弃风 % 整数变量 u_Chp binvar(1, T, full); % CHP开停机状态 % 辅助变量可选 C_sell sdpvar(1, T, full); % CO2外售量刚开始接触YALMIP的人容易在维度匹配上翻车。sdpvar(1,T,full)是行向量如果你某条约束不小心写出了列向量Matlab会直接报维度错误或者更隐蔽地做自动广播导致约束偷偷丢掉了。我的经验是所有变量统一用行向量遇到循环赋值时也要保证维度和索引方向一致。3.2 约束条件怎么写成YALMIP能消化的形式YALMIP写约束非常直白就像把数学公式抄下来一样。我举两个最容易出问题的例子。第一个是CHP可行域约束。假设抽凝式机组的可行域由四个不等式围成% CHP可行域P_min a * H P P_max - b * H Constraints [Constraints, P_Chp P_Chp_min a * H_Chp]; Constraints [Constraints, P_Chp P_Chp_max - b * H_Chp]; Constraints [Constraints, H_Chp 0]; Constraints [Constraints, H_Chp H_Chp_max];注意这里的a和b如果是向量比如不同时段有不同的最大值系数那就直接写成向量乘法如果是标量YALMIP会自动广播这个不会出问题。真正容易出问题的是“if-else”逻辑。很多初学者习惯写循环for t 1:T if t 1 Constraints [Constraints, ...]; else Constraints [Constraints, ...]; end end这在YALMIP里其实是合法的因为t是Matlab的数值变量不是sdpvar循环可以展开。但一旦你试图把if作用于sdpvar比如检查“P(t) 0”那就完蛋了这等于问一个符号变量“你现在大于0吗”它没法回答。如果确实需要条件约束请用逻辑变量或者用大M法改写绝对不要尝试“if sdpvar”。第二个例子是P2G和CCS的物料平衡约束。这里我偷懒用了同步模型所以CO2平衡可以在每个时段单独闭合成% CCS捕碳 P2G耗碳 外售碳 Constraints [Constraints, C_capture r_CO2 * V_P2G C_sell];再加上P2G耗电约束% P2G耗电 产气量 / 总效率 E_P2G V_P2G / (eta_elec * eta_meth);如果你后续还要考虑储气罐容量那就要引入一个状态变量S_tank并加上状态转移约束S_tank sdpvar(1, T1, full); Constraints [Constraints, S_tank(1) S_tank_init]; Constraints [Constraints, S_tank(2:T1) S_tank(1:T) V_P2G(1:T) - V_out(1:T)]; Constraints [Constraints, 0 S_tank S_tank_max];这里的时间索引错位是很多朋友容易搞晕的地方。我的习惯是把S_tank定义成T1长度用“t1时段结束后的储气量”作为S_tank(2)这样循环递推写起来不容易串位。3.3 求解器配置与热启动技巧模型搭好之后求解是另一道坎。我用YALMIP的solvesdp或者optimize命令都行关键在sdpsettings里的配置options sdpsettings(solver, gurobi, verbose, 2, debug, 1); options.gurobi.MIPGap 1e-4; % MIP间隙 options.gurobi.TimeLimit 300; % 单次求解限时5分钟 options.gurobi.NumericFocus 2; % 数值稳定性 diagnosis optimize(Constraints, Objective, options);如果没有Gurobi或者Cplex的许可证也可以把求解器改成Matlab自带的intlinprogoptions sdpsettings(solver, intlinprog, verbose, 2);不过intlinprog对MIP问题的处理能力相比商业求解器有明显差距模型规模一大比如连续变量超过2万个整数变量超过500个就很吃力。所以我强烈建议在学生版或者项目预算允许的情况下用Gurobi或Cplex。如果都没有至少装个SCIPMIT的免费求解器也能解不少MIP问题。还有一个特别实用的经验热启动。如果你的模型要做很多组不同参数比如不同碳价、不同P2G容量的重复求解每次从头算很亏信息。YALMIP支持设置初始解assign(P_Chp, x0.P_Chp); % 用上一次的最优解作为初值 assign(H_Chp, x0.H_Chp); options.gurobi.MIPStart 1;这样做的好处是在参数变化不大的情况下求解器能在几秒内收敛到新最优解整体计算时间能下降50%以上。这个套路在很多学术论文里不会写但实际做项目是真的香。3.4 结果可视化不只画曲线还要画“多能流平衡”画图这块常规做法是把电出力、热出力、P2G产气量各画一张图但我觉得更有价值的是画“平衡图”和“碳流图”。平衡图是把每个时段各个能源的供给和需求用正负柱状图展示出来一眼就能看出哪个时段缺电、哪个时段弃风、哪个时段CCS捕碳量最大。碳流图则展示碳捕集、碳利用、碳外售三项的时序曲线。核心绘图代码大概长这样figure; bar(t_axis, P_Chp_opt, g); hold on; bar(t_axis, P_wind_opt, b); bar(t_axis, -E_P2G_opt, r); bar(t_axis, -E_CCS_opt, m); plot(t_axis, Load_elec, k-, LineWidth, 2); legend(CHP电量, 风电, P2G耗电, CCS耗电, 电负荷);用堆叠柱状图去观察系统平衡是最直接的。你会发现在午夜时段风电出力大而电负荷小传统模型里只能弃风但在P2G和CCS的模型里弃风被P2G和CCS消耗掉一部分整个系统的弃风率大幅下降。这个结果比任何成本数字都更能直观说明模型的价值。4. 优化求解与结果分析实战中的关键问题4.1 MILP求解时间上不去的常见瓶颈调度模型的数学本质是一个混合整数线性规划。T24小时、3台CHP机组、1台P2G、1台CCS整数变量的数量大概是3×2472个连续变量几百个这个规模Gurobi几秒钟就能搞定。但一旦把时间粒度从1小时改成15分钟或者把机组数量扩展到10台整数变量立刻涨到几千个问题就变成难啃的骨头了。我这次做的算例是48小时T48CHP机组3台整数变量144个Gurobi默认参数下通常需要30-60秒才能把MIPGap压到1e-4级别。如果时间也按15分钟细化我试过把T变成96求解时间直接飙到20分钟以上。这时候就需要调整策略关闭无关输出verbose降到1将MIPGap放宽到1e-3工程上已经完全够用使用分支优先级对启停变量赋予更高的优先级适当增加预求解强度减少冗余约束这些参数在sdpsettings里都能配。我不建议一开始就追求零MIPGap因为最后0.01%的成本差异在工程上毫无意义但耗时可能是数量级的差距。4.2 典型算例结果解读弃风降了碳排放降了成本未必涨我在一个典型冬季算例热负荷高、风电大里跑了优化模型对比“无P2G/CCS”和“有P2G/CCS”两种配置。结果非常有代表性无P2G/CCS时系统弃风率约12.6%碳排放量为基准值有P2G/CCS时弃风率降到4.1%碳排放量下降约22.8%系统总运行成本上升约3.5%主要来自P2G和CCS的设备运维成本和额外的能耗成本上升这个现象要正确看待。如果碳价够高比如超过某个阈值那么捕集碳带来的碳交易收益能够平衡设备成本总成本甚至可能更低。这就是碳价灵敏度分析的用处找出“碳价多高时P2GCCS划算”这个临界点。我实际算下来在这个算例参数下临界碳价大约是65元/吨左右高于这个值系统配置就有经济动力。这个数字虽然依赖具体参数但量级上可以给同行做一个参考。4.3 灵敏度分析别把参数拍脑袋定跑一圈才知道对工程报告或者学术论文来说单点最优结果的说服力远远不够评审人几乎必然要问“参数变化时结论还立不立得住”。我这次做了两组灵敏度分析第一组是碳价从30元/吨逐步升到120元/吨观察P2G产气量、CCS捕集量、总成本的变化。结果符合物理直觉碳价越高CCS捕集量越大P2G产气量也随之上升因为更多的CO2被循环利用到了甲烷化反应里。这个趋势曲线可以直接画出来作为“碳价驱动碳循环”的有力证据。第二组是P2G转化效率从50%扫到80%观察它对弃风消纳的影响。效率越高相同电量能产出更多的甲烷P2G运行的经济性越强系统最终会更多地把盈余电力转向P2G。这个分析能用来论证“电解水效率提升是P2G经济性的命门”也算是对未来技术路线的一个量化判断。在Matlab里做灵敏度扫描很简单包一层for循环改参数重新求解结果存成矩阵最后统一画图。注意每组参数求解前都做一次热启动整体耗时能控制在可接受范围。5. 常见问题与调试实录5.1 模型无可行解先把约束“松”到极致再逐步加上“Infeasible problem”是我在这类项目里被折磨最多的情况。模型第一次跑出个“无可行解”人很容易慌然后开始瞎改参数其实大部分时候是你的约束之间存在矛盾。我的调试方法是第一步把目标函数改成0只看约束是否有解第二步把储罐容量、CCS捕集量上下限这些“软约束”放宽到离谱的水平第三步保留核心的功率平衡约束去掉所有的耦合项比如把CO2平衡断开一旦找到可行区域再逐步把约束加回来在YALMIP里排查约束很简单当optimize返回infesible时可以用check(Constraints)查看每条约束的残差一般残差为负或者非数的那条就是问题所在。我遇到过最经典的问题是把P2G产气量上限设得太低同时又把CCS捕集量下限设得很高导致送不出去的CO2堆在系统里碳平衡跑不通。5.2 求解时间长三招就能解决大部分痛点求解时间过长是另一个高频问题。除了前面提到的调MIPGap和热启动还有一个非常实用的技巧明确指定整数变量的分支优先级。对于调度模型启停变量是最重要的整数变量可以给它们设置高优先级变量序号% 假设u_Chp是主要的整数变量 options.gurobi.VarBranch 1; % 手动分支规则不过说实话对于YALMIP用户来说这个操作有点绕一般人不太需要。我更推荐的做法是降粒度。如果你要跑30天的长期调度把1小时粒度放宽到2小时甚至4小时模型规模直接缩小一半以上结果趋势基本不变。如果你要研究日内精细调度再用1小时粒度跑一周内典型日两者是互补关系没必要在一个模型里把所有尺度都吃掉。5.3 几个让人哭笑不得的编码细节这节是我觉得最有价值的部分记录几个看着不起眼但真能“爆雷”的细节问题。第一YALMIP变量和非YALMIP变量混用。比如你有一个针对某些时段的常系数数组直接用普通Matlab数组算出来的结果和sdpvar放在同一个表达式里没问题但如果你不小心把某个普通数组与sdpvar做了点乘.*且数组里有NaN位置那约束就悄悄丢掉了。所以我的习惯是涉及sdpvar的表达式所有系数全都先用full()和double()强制转成普通数值避免出现符号变量悄悄变成了ssstate space对象的情况。第二维度陷阱。sdpvar(1,T)和sdpvar(T,1)长得像但完全不同。所有约束都要统一方向我干脆约定所有变量都是行向量所有系数矩阵也都是行向量只使用水平维度。第三Gurobi的许可证问题。在2024年之后Gurobi对学术使用是提供免费license的需要你在官网申请绑mac地址。Matlab的YALMIP调用Gurobi时还要设置path建议把gurobi的mex文件路径手动添加到Matlab搜索路径里。不然你solver明明配置了gurobiYALMIP还是会悄悄回退到默认求解器结果你看了半天结果才发现根本没调用到Gurobi。关于Matlab本身的版本问题我不建议一上来就追最新版稳定够用就好。R2023b和R2024a都可以YALMIP对这两个版本的兼容性都不错。如果你用的是R2026b这种未来版本大概率也没问题但建议先跑一下YALMIP自带的testyalmip脚本确认环境正常。5.4 快速问题速查表问题现象可能原因排查/解决办法返回infeasible功率平衡或碳平衡约束矛盾先移除耦合约束逐步加回用check()查残差求解时间暴涨整数变量太多、MIPGap过严降时间粒度、开放MIPGap、热启动结果出现NaNsdpvar维度不对或系数含Inf检查变量尺寸用size()逐个打印Gurobi不生效路径配置错误或license问题检查sdpsettings(solver,gurobi)看verbose输出碳平衡对不上甲烷化反应系数或捕集效率设置错手算几个时段物料平衡逐项核对YALMIP报“无合适求解器”缺少求解器或选项写错安装Gurobi/SDPT3或改用intlinprog写在最后我的几点实在体会这项目里最让我印象深刻的一条经验是别急着上复杂算法先把“潮流平衡”这种最基本的物理约束校核清楚。很多人在建模时过度关注目标函数和求解器调参忽视了碳元素守恒、能量守恒这些底层逻辑结果模型跑得很“顺利”但结果一看就是瞎编的。我在调试阶段就发现P2G和CCS如果没有严格按化学反应计量比挂钩系统会多出“免费”的甲烷和“凭空消失”的二氧化碳成本低得离谱却完全不符合现实。所以不管你是不是要用Matlab先把平衡关系想清楚建模就成功了一大半。另外一个小建议参数不一定要完全按文献里拿遇到不确定的地方先给定一个“合理值”跑通全流程再做灵敏度分析看结果敏感性。比如P2G的电解效率文献里从60%到80%都有你与其纠结哪一个更准确不如把模型里的参数设成变量跑一组扫描曲线这样审稿人或项目经理要什么口径你都能给出来。如果后续想做更深入的工作可以把单目标优化扩展成多目标优化比如“成本碳排放”双目标也可以把确定性优化升级为分布鲁棒优化或随机优化加入风电出力的不确定性集合。模型框架和Matlab代码的基础我已经搭好扩展不会太痛苦但那是下一个项目的话题了。这次就先分享到这里希望这篇记录能帮你少走几个我走过的弯路。