做综合能源系统优化这几年我接手的项目里出现频率最高的关键词基本就是综合能源系统、电转气、碳捕集系统、热电联产这四个词的任意组合。原因不复杂大家都在找一条既能消化富余可再生能源、又能压低系统碳排放的可行技术路径而把这几个词串起来的思路恰好是——用P2G电转气把多余的电变成天然气用碳捕集系统把热电联产机组排放的CO₂捕回来再作为P2G甲烷化的原料形成电-热-气-碳的耦合闭环。这篇内容不打算复现某篇具体论文而是根据这套技术路线从机理建模到Matlab优化求解完整梳理一遍我实际做过的建模框架、约束构建方式和踩坑经验适合正在做综合能源系统调度的研究生、工程师以及想入门低碳优化方向的朋友参考。1. P2G与CCS耦合到底在解决什么现实问题先说清楚这套系统要解决的真实痛点不然很容易陷入为了建模而建模的状态。综合能源系统的核心枢纽通常是一台热电联产机组CHP它同时供应电和热能量利用率看起来很高但实际运行中有两个长期被人诟病的毛病一是发电和供热强耦合二是碳排放量大。P2G和CCS就是在这个基础上被引入的。1.1 CHP的以热定电困境为什么必须给它开一扇新门传统供热期里CHP的运行方式是以热定电热负荷定了电出力就被热电比锁定在很窄的区间内。比如一台抽凝式机组当采暖抽汽量固定后低压缸发电量基本没多少调节空间。结果就是在夜晚风电大发、热负荷又高的时候CHP必须维持较高电出力风电只能弃掉反过来白天热负荷降低CHP电出力又跟着往下走可能不够满足电负荷。这背后的本质是CHP的电-热可行域是一个刚性比较强的区域而不是一条可自由伸缩的曲线。新能源渗透率越高“以热定电”和“可再生能源消纳”之间的矛盾就越尖锐。光靠电储能去消纳风电成本高而且容量往往不够光靠提高CHP调峰能力又受制于汽轮机本身的机械和热力特性。这时接入P2G相当于给系统开了一个“电→气”的泄压阀电多了不弃掉而是转化为天然气存起来需要的时候再去发电或供气。1.2 P2G与CCS的技术路线氢、碳、气如何闭环P2G过程通常分两步第一步电解水制氢把电能转化为氢气的化学能第二步氢气与二氧化碳在甲烷化反应器中合成甲烷也就是人工合成天然气。这一步很关键因为纯氢的储存和运输成本太高而甲烷可以直接注入天然气网或者用常规模拟的储气罐存储。CCS的引入则补齐了甲烷化的原料短板。CHP燃烧天然气或煤时烟气和尾气中含有大量CO₂通过胺法吸收、再生等工艺把CO₂从烟气中分离出来一部分封存或利用另一部分直接送到甲烷化反应器与电解水产生的氢气反应生成甲烷。这样一来链条就变成一个闭环CHP发电产生CO₂ → CCS分离CO₂ → CO₂与P2G的氢气甲烷化 → 甲烷回到气网/储气罐 → 再次供给CHP发电。这个闭环最有意思的地方在于CCS本身也是一个耗电设备它在白天多耗电、在夜间给风电让出更多消纳空间到了低谷时段P2G再大功率投入。两个原本看似独立的新技术因为“碳-氢-气”在分子层面的转化被绑定在一起也让整个系统的运行自由度变大了。1.3 耦合建模的整体思路从孤岛设备到统一能量流做优化模型时最容易犯的错是把每个设备单独建模最后凑一堆互不相干的方程。正确的做法是先画出统一能量流框架电母线连接CHP电出力、风电、P2G耗电、CCS耗电、电负荷、电储能。热母线连接CHP热出力、燃气锅炉补热、热负荷、储热。气母线连接P2G产气、天然气网购气、储气、CHP耗气、气负荷。碳流线连接CHP排碳、CCS捕碳、P2G耗碳、碳封存/碳利用。有了这个框架每个设备的输入输出只是在这几条母线之间“搬运能量和物质”。Matlab建模时也按这个思路去定义变量和平衡约束代码结构会清晰很多。后面所有数学方程都围绕这套统一框架展开。2. 从物理过程到数学方程四大核心设备的建模要点这一部分是整个优化模型的地基。我在实际操作中体会是设备模型不要一味追求精度够用就行——优化调度关心的是“边界范围内的最优工况”而不是微观反应动力学。因此下面每个设备我都采用工程上常用、又能方便写入优化模型的数学表达。2.1 抽凝式CHP用可行域顶点描述热电之间的柔性区间抽凝式CHP是最常见的一类热电联产机组它的电出力P_chp和热出力H_chp不是简单线性关系而是一个凸多边形可行域。工程上最通用的做法是把可行域用几个顶点来表示然后通过一组线性不等式约束描述P_chp_min ≤ P_chp(t) ≤ P_chp_max H_chp(t) ≤ H_max_at_P_chp(P_chp(t)) P_chp(t) c_h · H_chp(t) ≤ c_some实际上我会直接用顶点枚举把可行域写成矩阵形式A_chp * [P_chp(t); H_chp(t)] ≤ b_chp比如某台机组可行域可以近似为四个顶点(A点)最小电出力最大热出力、(B点)最大电出力最大热出力、(C点)最大电出力零热出力、(D点)最小电出力零热出力。每个时段机组运行状态必须落在这个四边形内部。顶点约束的好处是不管机组实际特性多复杂只要能抽出若干关键运行工况点就可以线性化表示后面交MILP求解器处理几乎无压力。注意如果机组是背压式热出力与电出力是固定比例那就更简单了直接一个线性比例约束搞定。但抽凝式的通用性更强我建议模型直接按抽凝式写一个函数背压作为抽凝的特例传入参数即可。2.2 P2G电解槽变工况效率与甲烷化的化学计量约束P2G的耗电功率P_p2g和产气量之间有一个效率曲线。电解槽在额定负荷附近效率最高低负荷和过载效率都会明显下降。这个特性如果简单用固定效率表示会严重高估低谷时段的产气量。我通常的做法是把P2G整体效率曲线作分段线性化G_h2(t) η_elec_seg * P_p2g(t) 按分段区间取对应效率然后甲烷化环节再乘以一个反应效率η_meth最终产甲烷功率G_ch4(t) η_meth * G_h2(t)化学计量关系也要进模型。简单写就是产生1MW当量甲烷大约需要消耗一定质量的CO₂。从反应方程 CH₄ 2O₂ ⇌ CO₂ 2H₂O 走反向甲烷化来看每公斤甲烷需要约2.75公斤二氧化碳。优化模型里我一般把这个系数换算成 t/MW当量得到一个常数系数k_co2M_co2_meth(t) k_co2 * G_ch4(t)这个约束把P2G和CCS焊在了一起——没有足够的CO₂P2G产气量就会被卡住这也是很多人在代码里漏掉的关键耦合约束。2.3 CCS捕集率、再生能耗与0-1运行状态建模CCS系统在模型中主要体现为两个量捕集到的CO₂量M_capture(t)以及为此付出的运行电耗P_ccs(t)。通常用一个线性关系近似能耗中包含固定能耗和可变能耗两部分。固定能耗来自溶剂循环、泵等设备只要系统开启就在消耗可变能耗则随捕集量上升而上升。为了刻画这一点我引入0-1变量u_ccs(t)0 ≤ M_capture(t) ≤ M_capture_max * u_ccs(t) P_ccs(t) P_ccs_fixed * u_ccs(t) p_rate * M_capture(t)u_ccs(t)同时限制了最大捕集率。捕集率不能超过进入CCS的烟气中CO₂总量于是还要写一条来源上限约束M_capture(t) ≤ co2_content * F_gas(t) * capture_ratio_max。这里的co2_content与CHP耗气量对应F_gas(t)是天然气进气速率。很多人会忽略CCS的另一个角色它是电力负荷。在优化模型里P_ccs(t)出现在电平衡方程右侧一方面增加系统电负荷拉低整体电平衡另一方面也给风电低谷提供了消纳出口。这一点后面在仿真结果里表现得很明显。2.4 储氢/储气/储热动态SOC方程和充放约束储能设备让优化模型从“静态平衡”升级为“跨时段协调”。我至少要建三类储能储气罐存储P2G产出的甲烷、储热罐存储CHP多余热量和蓄电/储氢设备。它们的统一形式是SOC(t1) SOC(t) η_charge * P_charge(t) - P_dis(t)/η_discharge - loss * SOC(t)对储热和储气罐充放功率上限、SOC上下限约束必须写全对储气罐如果气网可以接收P2G产气也可以把它算作系统外出口。储氢罐由于氢可以直接用于后续甲烷化也可以单独建模但为了模型简洁我更倾向于将氢视为中间变量不单独设储能产出的氢气直接进入甲烷化反应器除非研究目标明确涉及氢能市场再单独展开。储热罐对CHP运行有直接影响白天热负荷大、电负荷也不小CHP全开可能电多热少夜间热负荷大、电负荷小CHP如果按热平衡跑电又会过剩。有了储热罐就可以在白天让CHP多产热存储夜间放热从而减少白天购电或夜间弃风。3. 优化调度模型的搭建目标函数、约束顺序和线性化技巧设备方程都有了接下来就是把这些方程组织成一个可求解的优化问题。我采用的建模框架是混合整数线性规划MILP。这一节直接给出一个典型的24小时调度模型的完整数学结构以及我在整理约束时的一套顺序法则。3.1 24h调度的目标运行成本碳排放的量化权衡目标函数我习惯写成四部分之和min ∑_{t1}^{24} [ C_buy(t) * P_buy(t) C_gas(t) * F_buy(t) C_om * (P_chp(t)P_p2g(t)P_ccs(t)) C_carbon * E_net(t) ]C_buy(t)时段电价峰谷电价下是动态的。C_gas(t)网购天然气价格。C_om运维成本系数按出力线性折算。C_carbon碳排放成本系数可以理解为碳价或碳税E_net(t)是系统的净碳排放量。净碳排放量的计算是这条模型的灵魂。它不能简单等于CHP总排放而是要扣掉捕集量以及甲烷化消耗的碳E_net(t) E_chp_raw(t) - M_capture(t) - M_co2_meth(t)注意M_co2_meth取决于P2G产气量而M_capture由CCS决定。这两笔“扣减”分别代表“捕集后封存/利用”和“碳循环再利用”两种减碳路径。目标函数里加上C_carbon * E_net(t)C_carbon设得越高模型就越倾向于开启CCS、增加P2G产气用高成本换低碳排放。这就是经济性和环保性的显式权衡。3.2 约束的组装顺序与线性化处理要点约束多了以后代码能不能一次跑通很大程度取决于你按什么顺序写。我总结为以下顺序第一优先级母线平衡约束。电、热、气三条母线的守恒方程先写它们决定了系统能量是否闭合。第二优先级设备出力上下限与可行域约束。每台设备的P_min/P_max、爬坡、最小启停时间等这些约束独立于其他设备单独一块写。第三优先级设备间耦合约束。比如P2G产气量与CHP耗气的关联、CCS捕碳量与P2G耗碳量的关联、CCS电耗对电母线的影响。这部分是模型最容易出错的地方建议逐条校验量纲和方向。第四优先级储能动态约束。包含SOC更新、储充放上下限、初末SOC相等约束。初末SOC相等是为了保证调度结果可以循环执行否则模型会像“吃老本”一样把初始储能全部放空。线性化方面唯一真正棘手的是P2G效率曲线。解决方案有两种一是用SOS2约束做分段线性化MILP求解器原生支持二是预先计算各离散区间的线性系数用凸组合方式写。MatlabYALMIP里我通常用replacer或pw类工具或者干脆自己手写SOS2辅助变量这两种在工程上都稳定。3.3 求解范式为什么选择了MILP而不是启发式或纯非线性规划你可能要问为什么非要用MILP直接调fmincon不行吗我实测对比下来原因有三点。第一MILP可以保证全局最优性。多时段调度涉及储能跨时段耦合目标函数和约束中存在大量非线性关联用启发式算法容易陷入局部最优而且每次跑的结果还不一样论文里没法交代。第二设备状态的自然表达是整数变量。CCS起停、CHP启停、分段线性化选择段位这些天然是0-1变量。fmincon不能直接处理离散变量只能靠罚函数或分支割稳定性很差。第三现代MILP求解器的性能远远超出普通认知。Gurobi和Cplex对于这种规模规模通常在几千个变量、几千条约束以内的24小时调度模型通常几秒到几十秒内就能收敛到0.1%以内的最优gap完全满足研究和实时调度的需求。如果你的模型里确实存在不可线性化的非凸部分——比如P2G动态学中某些热力学细节——我会建议把这一部分拆出主优化循环作为一个“查表/内环评估器”处理主问题仍然保持MILP结构。4. Matlab代码实现从参数表到求解结果的可复现链路理论模型再漂亮最后也得落成能跑的Matlab代码。下面给出一个我在多个项目中反复使用的代码骨架所有代码都基于YALMIPGurobi的组合稍加修改也可以换成Cplex、intlinprog或开源求解器。4.1 基础数据定义把24个时段的负荷切面组织成参数结构体写优化代码之前我习惯先把所有外部数据装进一个结构体避免后面散落一堆全局变量。一个典型的参数定义如下%% 初始化参数结构体 dt 1; % 每个时段1小时 T 24; % 调度周期24小时 params.T T; params.dt dt; % 负荷曲线示例实际数据从Excel或预测模型读入 params.P_load [20 19 18 18 17 18 22 26 30 32 34 33 ... 31 30 32 35 37 38 36 33 31 28 25 22]; % 电负荷MW params.H_load [45 44 43 42 41 40 39 42 45 48 50 49 ... 48 47 46 45 44 45 46 47 48 47 46 45]; % 热负荷MW params.P_wind [10 12 14 16 18 20 18 15 12 10 8 7 ... 6 5 6 8 10 12 11 9 8 7 6 5]; % 风电预测出力MW params.price_elec [0.4 0.4 0.4 0.4 0.35 0.35 0.5 0.8 ... 0.9 1.0 1.0 0.9 0.8 0.8 0.7 0.6 ... 0.7 0.9 1.0 1.0 0.9 0.7 0.5 0.4]; % 分时电价元/kWh params.price_gas 0.35; % 天然气价格元/kWh当量这里的核心思路是所有随时间变化的外部输入全部按“1×T向量”存放时间索引统一到第1~24行后面写约束时直接params.P_load(t)引用出错率最低。4.2 YALMIP变量定义与约束写入的实操片段%% 定义决策变量 % CHP P_chp sdpvar(1, T); % 电出力MW H_chp sdpvar(1, T); % 热出力MW u_chp binvar(1, T); % 启停状态 % P2G与CCS P_p2g sdpvar(1, T); % P2G耗电MW G_ch4 sdpvar(1, T); % 甲烷产量MW当量 u_p2g binvar(1, T); % P2G运行状态 M_capture sdpvar(1, T); % 捕集CO2量t/h P_ccs sdpvar(1, T); % CCS耗电MW u_ccs binvar(1, T); % CCS运行状态 % 储能 SOC_soc sdpvar(1, T1); % 储气SOC P_cha_gas sdpvar(1, T); % 储气充入功率 P_dis_gas sdpvar(1, T); % 储气放出功率 % 同理还有储热SOC和蓄电SOC此处省略以保持篇幅 % 网购天然气、购电 F_buy sdpvar(1, T); % 网购天然气功率当量MW P_buy sdpvar(1, T); % 购电功率MW %% 约束集合 Constraints [];注意我把设备启停状态也建模为0-1变量。这样CCS的固定电耗就可以写成0-1状态×固定功率物理意义明确。代码写到这里变量的数量级和含义都应该心中有数。接下来是典型的约束写入方式我逐条列出%% 电功率平衡约束 Constraints [Constraints, ... P_chp P_wind P_dis_elec P_buy ... P_load P_p2g P_ccs P_cha_elec]; %% 热功率平衡约束 Constraints [Constraints, ... H_chp H_boiler P_dis_heat ... H_load P_cha_heat];这里有一个容易忽略的细节CCS耗电要放在平衡式的负荷侧。我刚接触时把它漏了导致每次结果里CCS都是满负荷运行因为模型认为它免费耗电直到一个月后才检查出来。4.3 求解器调用与结果解包optimize之后的一堆坑%% 目标函数 Cost_buy_elec params.price_elec * P_buy; % 购电成本 Cost_buy_gas params.price_gas * F_buy; % 购气成本 Cost_om 0.01 * (P_chp P_p2g P_ccs) * ones(T,1); % 运维费用 Cost_carbon params.carbon_price * sum(E_net); % 碳排放成本 Objective Cost_buy_elec Cost_buy_gas Cost_om Cost_carbon; %% 求解 ops sdpsettings(solver, gurobi, verbose, 1, ... gurobi.MIPGap, 1e-4); optimize(Constraints, Objective, ops); %% 解包结果 P_chp_opt value(P_chp); H_chp_opt value(H_chp); P_p2g_opt value(P_p2g);跑完之后第一个要检查的不是目标值是多少而是optimize返回的problem标志if optimize(Constraints, Objective, ops).problem ~ 0 error(求解失败请检查模型); endYALMIP的problem0表示正常1表示不可行2表示无界×表示数值问题。我强烈建议在跑大规模算例之前先把这个检查写在脚本开头否则后面出了莫名其妙的负数结果排查成本高得离谱。4.4 可视化电热平衡、碳排放与P2G运行曲线的画法拿到结果后我通常第一件事是画出三条曲线电平衡堆叠图、热平衡图、P2G与CCS运行状态图。堆叠图画起来很简单figure; bar(1:T, [P_chp_opt, P_wind, P_buy_opt], stacked); hold on; plot(1:T, P_load, r-o, LineWidth, 2); % 同时画一条P2GCCS耗电的下限看电平衡是否真的闭合但我更推荐画“平衡偏差检查图”把电平衡左右两侧相减理论上应全为0。如果有一个时段不为0说明约束漏了一条或变量序号错了。这种检查图在开发阶段比任何调试器都管用balance_check P_chp_opt P_wind P_dis_elec_opt P_buy_opt ... - P_load - P_p2g_opt - P_ccs_opt - P_cha_elec_opt; figure; plot(1:T, balance_check, k-); ylim([-1e-6 1e-6]); title(电平衡残差);锅炉、CHP热出力与热负荷之间也应画同样的残差图。做完这一步模型是否“闭合”就不言自明了。5. 典型日仿真结果三种方案对比揭示的耦合价值模型和代码都有了接下来用一个典型日算例说明结果。我强调一下以下数值来自我测试用的一组公开类型参数不是某篇论文的精确复现但趋势和结论具有普遍性可以作为你跑自己数据时的参考对照。5.1 方案设置与参数基线我设计了三个方案做横向对比方案A纯CHP无P2G、无CCS作为基线。方案BCHPP2G但CCS只捕集不参与甲烷化捕集到的CO₂直接封存外部碳汇。方案CCHPP2GCCS耦合捕集的CO₂优先供给P2G甲烷化剩余封存。CHP容量设为100MW电额定最大热出力60MWP2G额定功率25MWCCS最大捕集量20t/h碳价从0到200元/t以20元/t为步长变化。所有方案都允许购电购气储能设备保持一致。5.2 能流走向拆解多余的绿电最终去了哪里以方案C在中等碳价、夜间风电大发时段的能流为例此时电负荷约20MW热负荷约35MWCHP电出力被压低到45MW风电出力18MW电侧总供给63MW电负荷20MWP2G耗电15MWCCS耗电6MW储电充电2MW合计43MW——多余20MW体现为购电减少或CHP进一步压出力。P2G的15MW电投入后按全链效率55%折算产甲烷约8.3MW当量这里需要约0.75t/h的CO₂而CCS在此时段捕集量约1.2t/h完全够用剩余0.45t/h进入封存/利用。这条链路最直观的价值是原本被迫弃掉的18MW风电在方案C中超过一半被转化成了可储存的甲烷而这部分甲烷在第二天晚高峰又可以被CHP用来发电供热实现了“风光储”在时间维度上的平移。5.3 碳价敏感性成本与碳排放之间的权衡曲线碳价的变化会引导系统做出非常清晰的行为切换。我实测的趋势如下表所示碳价区间系统行为表现净碳排放运行成本0~20元/tCCS基本不开启P2G仅在谷电时段运行最高接近CHP原始排放的85%以上最低40~80元/tCCS开启但与P2G耦合度低捕集量逐步上升中等比基线下降30%~40%成本小幅上升100元/t以上CCS满负荷P2G全时段尽量高产气储气罐积极套利比基线下降50%以上成本明显上升但上升斜率放缓这个趋势说明一个问题CCS与P2G的耦合不是“装上就一定省钱”而是要取决于碳价信号是否足够强。碳价低时捕集的CO₂没有经济价值CCS就是个纯耗电设备碳价高到一定程度捕集的CO₂变成P2G的“免费原料”耦合系统才有动力满负荷运转。6. 建模和编程中容易踩的坑代码跑不通时先查这四类问题这节单独拎出来写因为我在这个方向上的时间至少有一半浪费在调试上。很多问题在理论推导时完全看不出来但一到代码层面就原形毕露。6.1 单位不一致是最隐蔽的无声错误这是所有错误里最可恨的一种因为Matlab不会报错但结果完全失真。典型场景CHP发电量用MW天然气购气量用m³/h碳排放用t/h目标函数的电价用元/kWh——四个不同的量纲放在同一个约束里直接加减得到的解根本不是任何物理情形下的最优解。我的做法是全模型统一采用能量当量所有电、热、气功率统一为MW所有物质的量统一换算为t/h或按化学反应系数折算为MW当量。引入的每个系数比如P2G的CO₂消耗系数都单独注释它的量纲来源。代码开头写一段自检assert(abs(k_co2 - 2.75/1.36) 1e-6, ... CO2消耗系数可能错误);量纲统一后再进模型很多“看起来合理但实际荒谬”的结果会自动消失。6.2 非线性项的线性化别让MILP悄悄变成MIQP我在模型里写过这样一个经验性约束储能SOC变化率与当前SOC成正比于是SOC²出现在约束里。YALMIP不会拒绝它会默默把问题类型改成MIQP混合整数二次规划求解时间从5秒变成50秒甚至更久gap还不一定收敛。排查方法很简单求解前打印YALMIP识别的问题类型确认它是MILP而不是MIQP。如果哪里突然变成二次了去搜一遍代码里的乘法——最常见的是两个sdpvar相乘比如sdpvar(1)*sdpvar(1)。这种隐性类型变化是性能杀手也是我在代码评审时第一眼看的地方。6.3 不可行解与无界解的排查从松弛变量到约束注释法一旦返回problem1不可行不要一上来就改模型。我有一套固定流程第一步用YALMIP的feasrelax去掉一部分约束看是哪些约束在冲突[primalfeas,feasflag,recoverdata,diagnostic] ... feasible(Constraints, ops);第二步采用“注释法”定位把约束分成几大组比如平衡组、设备组、储能组、耦合组先只保留平衡设备组求解能行再加储能组再能行加耦合组。第一次不可行的那一组就是问题所在。第三步检查物理量的方向。最经典的错误是变量符号定义相反导致平衡方程永远不可能成立比如把储能充电功率定义在负荷侧却用了流入方向。这套流程虽然“土”但比任何高级调试器都有效。我用它解决了不下二十个不可行问题。6.4 求解器选型Gurobi、Cplex、intlinprog与HiGHS的取舍最后说说求解器。学术环境最常见的是Gurobi和Cplex两家MILP性能都很强对几千变量规模的模型差异不大选哪个主要看是否有授权。没有商业授权时我的建议是小规模模型变量2000约束3000Matlab自带的intlinprog就够用完全免费。中大规模模型可以考虑开源求解器HiGHSYALMIP直接支持solver,highsMILP性能在开源里属于第一梯队。需要绝对稳定和良好WL松弛、指示变量、回调支持时再上Gurobi。另外无论用哪个求解器都要设置好MIPGap和TimeLimit。研究场景我通常设MIPGap1e-4如果只是做敏感性分析MIPGap1e-2可以显著提速误差成本完全可以接受。这套模型的后续扩展空间也很大。我在最新一版代码里已经加进了风速/负荷的不确定性场景把确定性MILP改成了两阶段鲁棒优化内层用强对偶把min-max问题转化为单层MILP求解。P2G和CCS的耦合正好为这类随机调度提供了天然的灵活性——低谷时段可以多储气高峰时段再多放气抵抗新能源波动全靠这个“气氢碳”缓冲池。如果这篇内容对你有用后面我可以把随机化的那部分约束也整理出来分享给大家。