1. 项目背景与优化问题定义1.1 为什么把P2G和碳捕集设备放在同一个系统里研究我做综合能源系统优化调度这块也有几年了坦白说早期看到P2G和碳捕集设备的组合模型时第一反应是“这不就是把两个模块拼一起嘛”但真正动手复现之后才发现这两个设备放在一起会产生很强的耦合效应远不是简单叠加。先说P2GPower to Gas电转气。它的原理分两步走第一步是电解水制氢2H2O电解生成2H2和O2第二步是氢气和二氧化碳发生甲烷化反应CO2加4H2生成CH4和2H2O。所以P2G设备本质上是一个“耗电耗碳、产出天然气”的环节。很多人做系统优化时把P2G简化成一个“电转气”的效率模型忽略了对CO2的消耗这在单独研究P2G消纳风电时问题不大但一旦和碳捕集设备联合建模就必须把这个CO2需求明确写出来因为碳捕集捕下来的CO2正好可以作为P2G的原料。碳捕集设备Carbon Capture SystemCCS项目里具体做的是燃烧后捕集则是另一个方向它从CHP机组或锅炉的烟气中把CO2分离出来。主流做法是化学吸收法用单乙醇胺MEA溶液吸收烟气中的CO2再通过加热再生塔把CO2解析出来。这个过程需要消耗大量蒸汽如果蒸汽来自CHP机组的抽汽那就会影响机组的热电出力分配这个耦合关系直接改变了系统的运行可行域。两个设备合在一起之后系统里就形成了一条很有意思的内部循环链路CHP燃烧天然气发电产热产生烟气烟气进碳捕集装置CO2被分离出来一部分卖给外部或者封存一部分送到P2G做甲烷化原料P2G产出的天然气再回到天然气网络里供气负荷使用。用我的话说这就等于在系统内部造了一条“碳的传送带”而整个传送带的驱动力来自风电和电网买来的电。碳排放成本和经济运行成本的目标对冲就藏在这条传送带的每一个环节里。1.2 双目标优化为什么选epsilon约束法而不是加权求和这个项目标题里写得很清楚要同时优化碳排放成本和运维成本。这是一个典型的多目标优化问题两个目标之间存在冲突想多捕碳、多产气就要多耗电多耗热运维成本就上去了想控制运行成本可能就要牺牲一部分碳减排效益。很多初学者拿到双目标问题第一反应是加权求和把两个目标乘上权重加起来变成一个单目标。我早年也这么干过但用在实际复现里容易碰到三个麻烦。第一个是量纲问题碳排放成本是吨乘以碳交易价格运维成本是设备运维加购能费用两个数量级可能相差好几倍权重系数很难解释清楚。第二个是权重本身的主观性拍脑袋定0.5和0.5整个系统的最优运行点就变了别人复现你的结果时没有办法还原你的决策逻辑。第三个也是更本质的问题加权法在Pareto前沿是非凸的情况下会漏掉一部分真实的Pareto最优解画出来的前沿曲线形状是残缺的。epsilon约束法完全是另一种思路。它把其中一个目标作为主目标另一个目标转成不等式约束。拿这个项目来说就是把碳排放成本当成主目标最小化把运维成本写成约束条件要求它小于等于某个给定的ε值。每给一个不同的ε就得到一个不同的Pareto最优解把ε从最小允许值一路扫到最大允许值就能得到一条完整的Pareto前沿。相比加权法epsilon约束法有几个突出的优点。它不需要预先确定权重避免了主观性它每次求解的都是单目标优化问题可以直接用成熟的求解器处理比如Gurobi或者CPLEX它生成的Pareto前沿点分布更均匀可控可以通过设定等间距的ε序列来控制解的密度。这个方法对非凸问题同样有效这也是为什么很多SCI一区文章在综合能源系统多目标调度上都采用它。1.3 复现前需要理清的系统边界和调度尺度在动手写代码之前我强烈建议先把系统边界画清楚。这个项目考虑的是典型日前调度场景时间尺度是24小时步长1小时也就是一共24个调度时段。系统内部的设备包括CHP热电联产机组、燃气锅炉、P2G设备、碳捕集装置、风电机组以及向上级电网和天然气网的购能接口。电负荷和热负荷是已知的预测值风电出力是已知的预测序列这些作为模型的输入参数传入。这里有一个细节值得注意。标题里强调了“热电联供”所以CHP在这个系统里承担的是供电压舱石和供热主力军的双重角色。它的电出力和热出力不是独立的这个耦合特性会在模型里体现为热电联产可行域约束。为了方便建模我在复现时把CHP按抽凝式机组处理也就是说在一定的电出力范围内热出力可以在一定区间内调节而不是固定热电比。这种处理方式能更真实地反映机组运行灵活性也是论文里比较常见的设定。系统与外部网络的交互也需要注意。向上级电网购电是有上下限的不能无限买向天然气网购气也有流量上限。这些边界条件在优化模型里都是硬约束如果漏掉任何一个求解出来的“最优方案”可能在工程上根本不可行。我在最初复现时漏掉了电网交互功率的爬坡限制结果算出来的24小时购电曲线在相邻时段出现了剧烈跳变后来补上约束之后曲线才变得合理。2. 综合能源系统建模与关键设备数学模型2.1 电-热-气三种能量流的总体架构整个系统的能量流动可以用三条主线来梳理。电流方面风电出力、CHP电出力、上级电网购电共同满足电负荷、P2G设备耗电和碳捕集设备耗电热流方面CHP余热和燃气锅炉的出力满足热负荷和碳捕集装置再生塔的耗热需求气流方面天然气网购气和P2G产气共同供应气负荷同时CHP和燃气锅炉消耗天然气。这三条能量流通过CHP、P2G和碳捕集装置实现耦合。用公式化的语言来写系统的电功率平衡约束$$P_{wt}(t) P_{chp}(t) P_{grid}(t) P_{load}(t) P_{p2g}(t) P_{ccs}(t)$$热功率平衡约束$$H_{chp}(t) H_{boiler}(t) H_{load}(t) H_{ccs}(t)$$气功率平衡约束$$G_{buy}(t) G_{p2g}(t) G_{load}(t) F_{chp}(t) F_{boiler}(t)$$三个公式看着简单但每一项背后都有设备模型在支撑。尤其是P2G产气量和碳捕集捕碳量之间不是独立的P2G产气需要消耗CO2而CO2的来源主要就是碳捕集设备这就在气平衡和碳排放平衡之间建立了一条隐含的耦合通道。我复现时的体会是建模阶段最忌讳的就是把每个设备孤立地写成一个效率模型然后简单叠加那样算出来的结果看似合理但实际上忽略了设备间的物质流约束换一组参数可能就崩了。2.2 CHP机组的运行可行域建模CHP在这个系统里是最核心的设备它的建模质量直接决定了整个优化结果的工程合理性。如果简单地写一个电出力上下限和热出力上下限那就等于把热电联产机组当成了两个独立设备完全没有体现“联产”的意义。我采用的是抽凝式CHP的线性可行域模型。这类机组的特点是在一定的电出力范围内可以通过调整凝汽流量和抽汽流量来灵活调节热出力但电出力和热出力之间存在一个耦合关系这个关系通常表示为一个多边形可行域。简化的写法是一组线性不等式约束$$P_{chp}(t) - c_v \cdot H_{chp}(t) \geq P_{chp}^{min}$$$$P_{chp}(t) - c_v \cdot H_{chp}(t) \leq P_{chp}^{max}$$$$0 \leq H_{chp}(t) \leq H_{chp}^{max}$$其中 cv 是电热特性系数。这组约束的物理意义是在纯凝工况下机组的电出力有上下限当从汽轮机抽汽供热时凝汽发电能力下降所以电出力和热出力之间存在一个反向的线性折减关系。虽然我在复现时对原论文的可行域做了一定简化没有考虑背压工况段但对日前调度这个尺度来说精度已经足够。另外CHP机组的爬坡约束也不能忽略$$-R_{chp}^{down} \leq P_{chp}(t) - P_{chp}(t-1) \leq R_{chp}^{up}$$如果不加爬坡约束优化器可能会让机组从一个极端工况瞬间跳到另一个极端工况这在数学上可行但在物理上不可实现。我在最初复现时偷懒跳过爬坡约束结果Pareto前沿上的某些解出现了很不合理的功率跳变后来补上之后情况立刻就好了。2.3 P2G设备的能量转化与CO2消耗关系P2G设备的建模看起来简单就是一个输入电功率、输出天然气功率的效率模型但里面有几个坑值得单独说。第一个坑是效率的层级。P2G的电转气全过程包含电解和甲烷化两个环节每个环节都有各自的效率。在日前优化调度中为了保持模型线性通常把两个环节合并成一个总效率η_p2g输入电能到输出天然气功率的关系写成$$G_{p2g}(t) \eta_{p2g} \cdot P_{p2g}(t)$$这个处理方式是论文里最常见的简化方式它隐含假设了电解槽和甲烷化反应器在中低负载率下效率基本不变。实际操作中这个假设是可接受的因为P2G在这个系统里主要承担的是消纳多余风电和利用捕集CO2的角色不是作为主力供气源。第二个坑是CO2消耗约束。P2G甲烷化反应需要CO2根据化学反应计量关系每生成1单位的CH4按能量计大约需要消耗固定比例的CO2质量。这个比例在模型中通常写成$$M_{co2, p2g}(t) \beta_{p2g} \cdot G_{p2g}(t)$$其中 βp2g 是CO2消耗系数。这笔CO2需求从哪里来一部分来自碳捕集设备不足的部分可能需要外购。如果系统里碳捕集的量不够P2G使用就会产生额外的碳源成本这个关系在做双目标优化的时候会自动体现出系统内部的权衡。我在代码里把CO2平衡写成了单独的约束块并把它与碳排放成本目标函数关联起来这样做的好处是后续如果想研究“碳捕集容量配置”或“P2G规模扩展”这类问题只需要改设备上限参数不需要动主体模型框架。第三个坑是P2G设备的运行上限和爬坡限制。P2G的电解槽虽然调节速度快但也不能在1小时内从0跳到满负荷我在复现时给它加了简化的爬坡约束。另外P2G的启停会影响电解槽寿命有论文做了启停惩罚项但在这个项目里为了保持模型简洁没有加入整数变量来精确建模启停状态而是通过最小出力限制来近似避免频繁启停。如果你的场景对设备寿命损耗很敏感可以考虑加入二进制变量代价是求解时间会明显增加。2.4 碳捕集设备的能耗与捕碳量模型碳捕集设备的建模核心要搞清楚三件事捕了多少碳、消耗了多少热、消耗了多少电。捕碳量取决于进入碳捕集装置的烟气量和捕集率。烟气量跟CHP机组和锅炉的燃料消耗量直接相关燃料烧得越多烟气CO2越多。简化模型里系统产生的CO2总量可以近似写成$$E_{total}(t) \alpha_{gas} \cdot (F_{chp}(t) F_{boiler}(t)) E_{grid}(t)$$其中 αgas 是天然气燃烧的碳排放因子Egrid(t) 是购电对应的间接碳排放。碳捕集装置捕集CO2的量是$$M_{capture}(t) \eta_{ccs} \cdot k_{capture} \cdot E_{total}(t)$$这里 kcapture 是碳捕集系统的CO2捕集率通常取0.85到0.9。被捕获的CO2分两路走一路作为P2G原料一路可以外送或封存。在本文模型里我优先让被捕集的CO2去满足P2G需求富余量才算作系统的碳减排量。碳捕集装置的能耗是模型里最容易被忽略、但影响又很大的部分。溶剂再生塔需要蒸汽这个蒸汽在本文模型中折算为热功率消耗 Hccs(t)而吸收塔的风机、泵等辅助设备需要消耗电功率 Pccs(t)。这两部分能耗都跟处理的CO2量成正比写成线性模型$$H_{ccs}(t) \lambda_h \cdot M_{capture}(t)$$$$P_{ccs}(t) \lambda_e \cdot M_{capture}(t)$$λh 和 λe 分别是单位捕碳量的热耗和电耗。这就是我之前说的耦合来源——碳捕集设备自己会“吃掉”一部分热和电而这部分能量供给又分流了CHP和燃气锅炉的出力对电负荷平衡和热负荷平衡都会产生压力。我复现时曾经想过把捕集设备的能耗设成常数拉倒结果发现捕碳量大时能耗占比相当可观低负荷时又显得浪费。后来坚持用与捕碳量线性相关的模型整个系统的运行结果合理多了。2.5 储能设备和其他系统约束原项目标题里没有提到储能设备但我在复现时发现如果不加任何储能系统的灵活调节能力会明显受限尤其是当夜间风电出力大而电负荷小时多出来的电只能通过P2G消化P2G容量不够就得弃风。如果你参考的论文里带了电储能和热储能那模型会更丰富但要做的工作量也会大不少。考虑到原项目是复现文章基础版配置不含储能我在主体代码里没有加入电储能和热储能而是用P2G作为柔性负荷来提升风电消纳率。在讨论部分如果后面想做扩展把储能加进去是比较自然的方向。其他系统约束还包括电网交互功率上下限、气网购气量上下限、风电出力上限即预测可用功率、P2G和碳捕集设备的容量上限以及各设备出力的非负约束。这些约束不是凑数的每一个都对应着实际工程中的物理边界。写模型的时候我建议把所有变量统一一个备用命名规则比如P_开头的都是电功率H_开头的都是热功率G_开头的都是气流量M_开头的都是CO2质量流这样调试代码的时候能省下大量时间。3. 双目标优化模型与epsilon约束算法实现3.1 目标函数一碳排放成本的定义与计算方式碳排放成本是本项目的核心目标之一。它不是简单地把系统总碳排放量乘以碳价而是要引入碳配额的概念这是很多论文里都在用的机制。系统实际碳排放量超过配额的部分需要花钱购买碳排放权低于配额的部分则可以出售获利。碳配额的设定方式直接影响了碳排放成本的计算结果。常见的做法是按机组出力和发电量的比例分配免费配额CHP机组发电越多分到的配额越多这是对供热为主的机组的一种政策倾斜。我在复现时采用的配额模型是$$E_{quota}(t) \gamma_{elec} \cdot P_{chp}(t) \gamma_{heat} \cdot H_{chp}(t)$$其中 γelec 和 γheat 分别是单位电出力和单位热出力对应的碳配额系数。碳捕集系统捕集的CO2从总排放中扣除所以系统实际需要购买配额的碳排放量为$$E_{net}(t) E_{total}(t) - M_{capture}(t)$$最终碳排放成本目标函数$$f_1 \sum_{t1}^{24} c_{co2} \cdot (E_{net}(t) - E_{quota}(t))$$c_co2 是碳交易价格。这个模型算出来的碳排放成本可以为正也可以为负。如果碳捕集量足够大E_net 小于配额那系统实际上是卖碳配额获利碳排放成本就是负数这是完全合理的代表低碳运行带来了经济收益。这一点在结果讨论时很值得讲因为它直观展示了碳捕集设备的经济价值。3.2 目标函数二系统运行维护成本的定义第二个目标函数是系统运行维护成本在优化模型里它通常包含三块向上级电网购电的费用、向天然气网购气的费用、各设备运行维护费用。购电费用是分时电价结构我采用峰谷平三段式电价这在复现案例里比较常见。购气费用则按照天然气价格乘以购气量计算。设备运维费用采用单位出力的运维成本率乘以出力的形式包括CHP、燃气锅炉、P2G和碳捕集设备各自的运维费用。综合起来$$f_2 \sum_{t1}^{24} [c_{buy}(t) \cdot P_{grid}(t) c_{gas} \cdot G_{buy}(t) \sum_{i \in \Omega} c_{om,i} \cdot P_i(t)]$$其中 Ω 是所有设备的集合。注意我这里用P_i表示设备出力对于产热设备在代码里要换成热量单位否则量纲会出问题。写代码时我建议把目标函数二单独封装成一个函数方便调试和对比不同参数下的成本结构。有一点要提醒碳排放成本里已经包含了碳交易费用所以运维成本里不要再重复加任何碳相关费用否则两个目标的物理含义就说不清了。这也是双目标建模时常见的边界混淆陷阱。3.3 epsilon约束法的原理与完整伪代码epsilon约束法的本质是将多目标问题转化为一系列带约束的单目标问题具体到这个项目就是保留碳排放成本作为主目标把运维成本加上一个上限约束$$\min f_1(x)$$$$\text{s.t.} \quad f_2(x) \leq \varepsilon_i$$$$x \in \Omega_{feasible}$$其中 Ω_feasible 是原问题的所有约束集合。每取一个不同的 εi 值就求解一次上面的单目标问题得到一组最优解和对应的两个目标函数值。把所有的 εi 遍历完就得到了Pareto前沿的离散近似。写代码的时候ε序列怎么取是需要认真考虑的。通常的做法是先单独最小化f1得到最优值f1_min同时记录这个解对应的f2值记为f2_max因为在这个解下f2一般是最大的再单独最小化f2得到f2_min同时记录对应的f1值。然后就可以确定ε的取值范围是[f2_min, f2_max]把这一整个区间均匀分成N段就得到了N1个ε值。我做的算法伪代码如下输入系统参数、负荷数据、风电数据 输出Pareto前沿点集合 Step 1: 求解约束化问题 求解 min f1仅含原始约束得到f1_min和对应的f2_max 求解 min f2仅含原始约束得到f2_min和对应的f1_max Step 2: 确定epsilon序列 设区间步长delta (f2_max - f2_min) / N 生成序列 epsilon_i f2_min i * delta, i 0, 1, ..., N-1 Step 3: for i 0 to N-1: 求解 min f1 约束条件 原始约束 f2 epsilon_i 记录解 x_i计算f1_i和f2_i 将(f1_i, f2_i, x_i)加入Pareto集合 Step 4: 对Pareto集合做非支配排序剔除被支配解 Step 5: 计算折中解采用模糊隶属度法这里有一个细节值得多说一句。第5步非支配排序不是可选项。由于整数变量的存在某些ε值下求解出来的点可能不是严格最优的或者出现两个不同ε对应同一个解的情况导致前沿上出现冗余点。做一遍非支配排序可以把这些冗余点剔除让Pareto前沿更干净。我在第一次写代码时跳过了这一步画出来的前沿曲线有毛刺加上之后平滑了很多。3.4 折中解的选择方法与工程意义Pareto前沿上的每个点都是数学意义上的最优解但工程上最终只能选一个调度方案执行。怎么选这里需要引入一个“折中解”的概念最常用的方法是模糊隶属度函数法。对Pareto前沿上的每一个解k计算它在第m个目标上的满意度$$\mu_{k,m} \frac{f_{m}^{max} - f_{k,m}}{f_{m}^{max} - f_{m}^{min}}$$然后把每个解对所有目标的满意度取平均$$\mu_k \frac{1}{M} \sum_{m1}^{M} \mu_{k,m}$$满意度最大的解就是折中解它的物理意义是两个目标的综合表现最均衡。在本文的双目标问题中折中解对应的调度方案既不会太追求碳减排而让运维成本飙升也不会太抠成本而让碳排放失控。我用这个方法选出的折中解两个目标值大概都处于各自极值区间的中间偏低碳一侧。这也是符合直觉的因为碳交易价格给定的情况下适当的碳捕集投入会带来碳配额的收益相当于给系统增加了一个收入来源所以折中解不会极端偏向单纯的经济最优。4. Matlab代码实现与调试要点4.1 开发环境与工具箱配置整个项目我用的是Matlab R2023a版本建模用YALMIP工具箱求解器用Gurobi 10.0。YALMIP是Matlab环境下非常好用的优化建模语言它的设计思路是让你用接近数学公式的方式写优化模型然后自动调用后端求解器求解。安装配置方面有几点经验分享。第一Gurobi安装后需要获取license学术用户可以直接申请免费license个人用完全够。第二Matlab里要用Gurobi需要运行gurobi_setup脚本把路径加入Matlab环境中然后在YALMIP里用solvesdp或者optimize时指定solvergurobi。第三YALMIP版本不要太老老版本对Gurobi新版本的支持有问题我遇到过接口不兼容导致Gurobi一直报错的情况升级YALMIP后就好了。如果你的电脑没装Gurobi用YALMIP自带的求解器比如sedumi或SDPT3也能算但对大规模MILP问题求解速度会慢很多。这个项目是24时段的混合整数线性规划变量规模不算特别大用sedumi也许能跑但我建议直接用Gurobi效率高一个量级。4.2 决策变量与主程序框架设计我把整个程序分成四个模块数据参数模块、模型定义模块、epsilon求解模块、结果分析模块。模块化的好处是后续调整碳价、负荷曲线或者设备参数时不需要动模型主体代码。数据参数模块用结构体存储所有参数包括系统参数、设备参数、电价气价、负荷曲线、风电预测数据等。我专门用一个init_parameters.m文件来写方便对比不同场景时直接修改参数文件。模型定义模块的核心代码如下% 决策变量定义以24时段为例 P_chp sdpvar(1, 24, full); % CHP电出力kW H_chp sdpvar(1, 24, full); % CHP热出力kW P_boiler sdpvar(1, 24, full); % 燃气锅炉热出力kW P_p2g sdpvar(1, 24, full); % P2G输入电功率kW G_p2g sdpvar(1, 24, full); % P2G输出天然气功率kW P_grid sdpvar(1, 24, full); % 电网购电功率kW G_buy sdpvar(1, 24, full); % 气网购气流量kW M_capture sdpvar(1, 24, full); % 碳捕集量kg P_wt sdpvar(1, 24, full); % 风电出力kW注意这里我用sdpvar而不是binvar意味着这些变量默认是连续变量。CHP的启停状态在这个简化模型里没有显式建模因为阶梯式最小出力可以通过可行域约束来近似。如果你参考的论文里考虑了启停你需要额外引入binvar变量模型的求解复杂度会上升不少。4.3 约束条件的代码实现细节约束条件在YALMIP里就是一堆等式和不等式逐条写上去就行。但有几个细节要注意。第一个是气平衡约束里的单位统一。天然气有体积立方米和能量kWh两种计量方式模型里我用的是能量单位kW气体负荷和P2G产气都用kW表示。如果你从外部数据源拿到的气负荷是立方米需要乘以天然气热值换算1立方米天然气约9.7kWh再传入模型。第二个是电功率平衡约束里碳捕集设备的电耗不是独立变量它是M_capture的函数所以不需要单独定义变量直接在约束里写成lambda_e乘以M_capture就行。同理热平衡约束里碳捕集热耗也是M_capture的函数。第三是P2G的CO2消耗约束这里要连接M_capture和G_p2g。代码里我写的是一个等式约束M_co2_p2g beta_p2g * G_p2g; % CO2消耗量 Constraints [Constraints, M_co2_p2g M_capture]; % 捕集CO2优先给P2G直接写成等式的话系统会强制所有P2G消耗的CO2都来自碳捕集这在建模里有点太紧了万一捕集量不够P2G就不能满发导致弃风加剧。我在复现时给P2G留了一条外购CO2的口子M_co2_buy sdpvar(1, 24, full); Constraints [Constraints, M_co2_p2g M_capture_to_p2g M_co2_buy];然后在目标函数里对外购CO2加一个惩罚成本。这样做的好处是模型更贴合实际工程缺点是目标函数里多了一项碳排放成本的含义会稍微变化。如果只想忠实复现论文不加外购CO2的口子问题也不大因为大多数论文里都默认捕集CO2足够P2G使用。4.4 epsilon循环求解的代码框架epsilon循环是整个程序的核心它把单目标优化问题包在循环里反复求解。我写的代码框架大致是这个样子% 第一步求两个单目标极值 optimize(Constraints, f1, ops); f1_min value(f1); f2_at_f1min value(f2); optimize(Constraints, f2, ops); f2_min value(f2); f1_at_f2min value(f1); % 第二步确定epsilon序列 N 20; epsilon_seq linspace(f2_min, f2_at_f1min, N); % 第三步循环求解 pareto_f1 []; pareto_f2 []; for i 1:N epsilon epsilon_seq(i); Constraint_eps [Constraints, f2 epsilon]; optimize(Constraint_eps, f1, ops); pareto_f1 [pareto_f1, value(f1)]; pareto_f2 [pareto_f2, value(f2)]; end有一个小坑要注意linspace生成的epsilon序列如果N取得很大比如50以上相邻两个epsilon对应的最优解差异可能很小求解时间成倍增加但Pareto前沿的改善很有限。我试过N30和N15前沿曲线差别肉眼几乎看不出所以一般取N15到20就够用了。另外每次求解前不需要清理工作区变量但建议在循环体内把上一次的optimize结果用assign函数重新赋值避免YALMIP里面有些变量的初始值残留影响求解。我在写第一版代码时因为没注意这个问题出现了偶发的不收敛排查了很久才发现是上一次求解的初始点没有清理。5. 算例分析与结果讨论5.1 案例系统参数设置与数据说明复现论文时算例参数设置是最让人头疼的部分。原论文的详细参数往往藏在附录里公众号或者复现包里一般只给一个大概的配置。我的做法是参考原论文里给出的系统结构图和数据表格尽量还原一个典型冬季日的调度场景。我这里提供一组基准测试参数供没有原论文数据的读者参考参数项数值单位CHP容量2000kWCHP电效率0.35-CHP热电比范围0.5~1.2-燃气锅炉容量1500kWP2G容量800kWP2G综合效率0.60-碳捕集率0.88-风电装机1200kW碳交易价格120元/t天然气价格3.2元/m³电网购电上限1500kW碳排放因子气2.16t/MWh电网碳排放因子0.58t/MWh电负荷和热负荷数据我采用的是典型冬季日的24小时曲线峰值电负荷约2800kW峰值热负荷约2200kW风电出力夜间高白天低有明显的反调峰特性。这些数据都放在init_parameters.m里一行一个数组方便替换。5.2 Pareto前沿的形态与关键解读跑完epsilon循环之后把得到的(f1, f2)点画在二维坐标上就是Pareto前沿。在我的基准算例下前沿曲线呈现典型的下降趋势碳排放成本从低到高变化时运维成本从高到低变化两个目标呈现明显的冲突关系。曲线左侧碳排放成本低、运维成本高对应的是碳捕集设备全力运行、P2G尽可能多地消纳风电并制气的状态。这个时候CHP可能会降低出力由电网购电和风电来平衡负荷因为CHP燃烧天然气会产生大量CO2增加碳排放成本。但同时电网购电成本高P2G运行也要耗电所以运维成本上去了。曲线右侧碳排放成本高、运维成本低则是相对保守的运行状态。碳捕集设备可能低负荷运行P2G也不太开系统主要靠CHP和燃气锅炉供能从电网购电尽量少这样设备运维和购能成本最低但碳排放量最大碳排放成本自然就高。折中解落在前沿曲线的中间偏左位置。用模糊隶属度法选出来的折中解碳排放成本比最小值大约高百分之二十左右运维成本比最小值大约高百分之十左右这个比例在不同碳价下会变化。碳价越高折中解越偏向低碳侧这说明碳交易机制对系统低碳转型有明确的引导作用。5.3 不同碳价下的系统运行策略演变我额外做了一组敏感性分析把碳交易价格从60元/吨一直加到200元/吨观察折中解对应的各设备出力变化。碳价低的时候系统倾向多开CHP因为碳排放的惩罚成本低CHP的效率优势能充分体现P2G和碳捕集设备的利用率都不高。碳价升高之后系统开始减少CHP出力增加电网购电同时碳捕集设备和P2G的利用率明显上升。碳价在150元/吨以上的时候P2G基本处于满发状态碳捕集设备也接近上限运行系统几乎把所有能捕的碳都捕下来了。这个结果从工程角度看很有参考价值。它说明碳捕集设备和P2G的投资是否能收回成本很大程度上取决于碳交易价格水平。碳价太低的时候这套设备就是纯成本项碳价够高它们才真正变成盈利单元。在做方案规划时这些设备配置的边界条件就是碳价而不是单纯的技术可行性。5.4 风电消纳与碳捕集的联动效应另一个值得关注的结果是风电消纳率。由于P2G的存在夜间风电大发时P2G会大量消纳风电来制气这直接提升了系统整体的风电消纳率。在折中解下风电消纳率可以达到95%以上接近完全消纳。更有意思的是风电消纳和碳捕集之间存在一个双向促进的联动效应。P2G制气需要CO2而碳捕集设备捕碳越多P2G的原料就越充足。当碳价较高时系统愿意投钱让碳捕集设备满负荷运行捕下来的CO2一路送进P2GP2G制气再替代一部分外部购气形成了“高碳价→多捕碳→多产气→少购气→更低碳排放”的正向循环。这条循环链在结果图上看得很清楚碳价从80元/吨升到160元/吨的过程中购气量下降的幅度非常可观。这个联动效应在单目标优化模型里也能体现但双目标Pareto前沿把它呈现得更直观——前沿上的每一个点实际上是这个循环链在不同“强度”下的运行结果。从低碳端到低成本端对应的是这条循环链从“全力拉满”到“基本关停”的渐变过程。6. 复现踩坑记录与后续扩展方向6.1 建模阶段的高频错误我在复现过程中踩了不少坑有些问题折腾了大半天才找到原因。这里整理一个速查表给后来者省点时间。问题现象根因解决办法求解结果出现负的P2G产气量变量没有加非负约束给所有出力变量加P0约束CHP热出力极端跳变缺少热出力爬坡约束在约束组中补充H_chp的爬坡限制Pareto前沿上出现大量重复点没有做非支配排序输出前增加去冗余步骤目标函数值出现NaN单位不统一导致模型病态统一kW、kg、元三种单位体系epsilon序列超出可行范围f2_max取的是f1最小值解对应的f2而不是f2上界先验证f2_max的可行性再生成序列Gurobi求解失败代码10005YALMIP版本与Gurobi接口不兼容升级YALMIP到最新版第一个问题其实是最常见的新手写约束时容易忘记“出力非负”这种默认常识。YALMIP里如果变量没有加非负约束求解器可能会让变量取负值来降低目标函数结果就是负产气量这在物理上完全说不通。6.2 求解性能优化的实际经验这个项目的MILP模型规模不算大变量数量大约在两三百个量级约束条件几百条Gurobi求解单目标问题通常几秒到几十秒就能收敛。但epsilon循环要跑20次加上两个极值点的求解整体耗时大约几分钟到十几分钟。如果想进一步提高求解效率我有几个实践经验。第一打开Gurobi的MIP gap容差限制设到1%甚至2%求解时间能缩短一半以上对日前调度结果的影响很小。第二给整数变量提供好的初始解可以先用一次不加整数约束的松弛求解结果作为初始点MIP的搜索会快很多。第三如果机器内存不够或者求解时间太长可以适当减少epsilon序列的个数N从20减到12前沿曲线形态变化很小但时间能省将近一半。6.3 从复现到扩展几个值得尝试的方向完成了基础版的复现读者完全可以在这个模型基础上做扩展。我列几个我认为有价值的扩展方向。第一个是加入储能设备。电储能可以平抑风电波动热储能可以解耦CHP电出力和热出力让CHP在碳价高时少发电、多产热配合储能满足热负荷。加入储能后风电消纳率和整体经济性都会有明显改善这是一个值得做的方向。第二个是把确定性优化换成鲁棒优化或随机优化。风电出力和负荷预测都是有误差的把不确定性描述成场景集或模糊集用分布鲁棒优化来处理能让调度结果在真实运行中更可靠。这也是目前综合能源系统研究的前沿热点。第三个是考虑碳捕集设备的动态特性。本文用的是线性静态模型但实际碳捕集设备的溶剂循环、吸收塔负载变化都有动态过程在更细的时间尺度比如15分钟下动态建模会显著影响碳捕集量和能耗的精确度。第四个是增加需求侧响应。电负荷和热负荷不完全是刚性的引入可平移负荷、可削减负荷让负荷侧参与系统调节能在不大幅增加运行成本的前提下进一步降低碳排放。6.4 关于论文复现心态的一点建议最后说点题外话。复现一篇SCI论文的代码很多时候最难的并不是优化算法本身而是弄清楚论文里没有明说的那些前提条件。我复现这个Energy一区文章时光是理解原系统的能量流框架和碳流框架就花了不少时间参数调校又花了不少时间。这个过程确实煎熬但一旦模型跑通Pareto前沿画出来你会有一种“终于拿到钥匙”的感觉。我的建议是不要指望一次性把整个系统全部跑通先从最基本的版本开始比如先用单目标优化只优化运维成本把模型跑通确认结果合理再加入碳排放目标、再把epsilon循环加上。每一步都验证一下输出结果的合理性这样出问题时你能快速定位是哪个环节的问题。我自己吃过这个亏一上来就写双目标加epsilon循环结果怎么调都出不来合理的Pareto前沿后来回退到单目标排查才发现是CHP可行域约束的符号写反了。这个项目后续还可以做很多花样但我觉得最核心的收获不是代码本身而是建模思路当你面对一个有多个设备、多股能量流、多种目标的系统时怎么把它分解成可求解的数学问题怎么判断哪些细节要精细建模、哪些地方可以简化处理这种判断力才是做研究工作真正需要积累的东西。