企业官网建设流程全解析

做了大半年方案对比和代码调试，终于把含P2G-CCS耦合、燃气掺氢和阶梯碳交易的虚拟电厂优化调度模型在Matlab里完整跑通了。这个项目解决的不只是“省成本”，而是把电、气、氢、碳四条能量流捏在一起做协同调度，让虚拟电厂既省燃料钱，又省碳配额钱，还能顺带消纳富余风光。整套代码用YALMIP建模、CPLEX求混合整数规划，跑一个典型日的调度计划只要一两分钟。这篇就把我的建模思路、关键约束、Matlab实现细节和踩过的坑都整理出来，给打算做碳交易、P2G、掺氢方向调度研究的朋友做个参考。

这篇东西适合两类人看：一类是做电力系统优化调度研究的研究生，想快速搭出含碳交易机制的VPP调度模型；另一类是在做综合能源系统可行性分析、需要算经济账的工程师。我尽量把从数学模型到代码实现的每一步都拆开讲，包括阶梯碳价怎么线性化、P2G和CCS耦合时的物料平衡怎么写、掺氢比例对燃气轮机出力范围的约束，以及这些约束在Matlab里怎么表示才不会让求解器报错。

1. 模型设计的底层逻辑：为什么要加碳交易和掺氢

1.1 虚拟电厂调度为什么难做

虚拟电厂（VPP）说到底就是把分散的风机、光伏、燃气轮机、储能、气网负荷聚合起来，作为一个整体参与电网调度。它跟传统火电厂调度最大的区别是：不确定性电源占比高、响应速度快的小型机组多、能量转换环节多。一个典型的VPP可能同时涉及电、气、热、氢四条能量回路，彼此之间还有耦合关系。

比如燃气轮机一边发电一边产热，余热还能补电热负荷；P2G设备用电制氢，氢气既能卖给加氢站，又能掺到天然气管网里给燃气轮机烧。这些耦合关系让调度问题的复杂度从单纯的“机组组合”升级成了“多能互补优化”，决策变量从十几个涨到上百个，约束也从线性等式变成了一堆非线性、含0-1变量的混合整数规划问题。

传统做法是只盯着运行成本做经济调度，把气、氢、碳的成本当成固定参数或者干脆忽略。但在双碳目标逐渐落到具体考核指标的背景下，碳成本变成了一个随排放量变化的分段函数，如果不用阶梯价格引导减排，模型算出来的“最优”方案很可能让他实际碳排放超标，后边碳配额的罚款比省下的燃料费还多。所以我在设计目标函数时，把“运行成本+碳交易成本”当成一个整体来优化，而不是分开两步走。

1.2 阶梯碳交易和传统碳交易的本质区别

传统碳交易机制里，碳价是一个常数，企业拿到免费配额后，每超出配额一吨就按固定价格买配额。这种机制的缺点是对深度减排的激励不够：如果碳价不高，企业宁愿花钱买配额也不肯改造设备，因为减排投入回收期太长。

阶梯碳交易机制相当于把碳价按照“超额排放量”分档：排放量超过配额越多，超出部分适用的碳价越高。举个例子，假设免费配额是100吨，阶梯设置了三个档次：

注意，阶梯碳价是“分档累进”的，不是“超额部分统一按最高档计价”。每一段的钱是分段计算的，加起来才是总碳交易成本。这样做的好处是让调度模型在每一档边界上产生明显的成本跳变，从而逼着优化器去寻求低碳运行方式。

在Matlab的优化模型里，我喜欢用一个叫carbon_tax的分段线性函数来计算碳交易成本，然后用0-1变量来标志当前所处的排放区间。模型求解后你会看到一个典型的“逼到墙角”的效果：只要风电出力够用，优化器宁可让机组低负荷运行也绝不让排放跨入第二档，这就是阶梯价格的杠杆作用。

1.3 P2G-CCS耦合是干什么的

P2G（Power to Gas）是用电能电解水制氢，再把氢气和捕集到的二氧化碳反应生成甲烷，这样电能被转化成了化学能储存在天然气里。CCS（Carbon Capture and Storage）是捕集燃气轮机和其它排放源产生的CO2，防止它直接排入大气。

单独看P2G和CCS各有各的问题：P2G制出来的氢气不好储存，直接掺入天然气管网有比例限制；CCS把CO2捕集后要运输到封存地，又是一笔费用。但把两者耦合起来，思路就打开了：

• P2G制取的氢气一部分直接用于燃气掺氢，一部分用于和CCS捕集的CO2反应生成甲烷；
• 甲烷可以注入天然气管网，替代一部分外购天然气；
• 燃气轮机燃烧甲烷时产生的CO2又被CCS捕集，捕集下来的CO2再次进入P2G的甲烷化单元。

这个循环本质上就是用“电力+捕集的碳”来制造“燃料”，实现了碳的循环利用而非简单封存。在模型里最直接的变化是：燃气轮机的碳排放不再是“排放多少就交多少碳税”，而是扣除CCS捕集量后的净排放才算碳账，这一下就改变了碳交易成本函数的结构。

我在代码里用一个eta_p2g参数表示P2G电解制氢效率，用eta_meth表示甲烷化反应效率，碳捕集量则根据燃气轮机出力乘一个捕集率来算。耦合的物料平衡是整个模型里最容易写错的地方，后面专门讲。

1.4 燃气掺氢为什么值得做

燃气掺氢是把P2G制得的氢气按一定体积比例混入天然气管网，供给燃气轮机燃烧。这么做有三个直观收益：

• 清洁能源占比提升：氢气燃烧产物主要是水，不会直接产生CO2；
• 为富余风光提供消纳路径：当风电光伏出力过剩而电网又限制出力时，可以用多余的电制氢，把弃电转成气存储；
• 降低对外购天然气的依赖：掺了氢气以后，同等的能量需求只需要更少的天然气。

但掺氢不是想掺多少就掺多少。燃气轮机对氢气的适应性和管网安全限制了掺氢比例，工程上一般体积比不超过20%，更高比例会导致燃烧回火、NOx排放升高、管道材料氢脆等问题。在调度模型里，我用一个h2_blend_ratio参数来约束掺氢上限，同时还要注意：掺氢后天然气-氢混合气的热值变了，燃气轮机的出力范围也会跟着变。

之前我处理这个环节时吃过亏，直接用固定热值算燃气轮机的耗气量，后来发现掺氢比例提高后，同样的功率输出对应的混合气流量其实增加了，如果不做热值折算，气网平衡的约束就会出错。正确的做法是：把掺氢比例映射到等效天然气流量，即天然气的等效消耗量=混合气流量×混合气热值÷纯天然气热值。

2. 数学模型搭建：目标函数、约束和碳交易线性化

2.1 目标函数：运行成本加碳交易成本

调度目标函数我写成最小化总成本，包含T个调度时段的四项成本：

• 购电成本和售电收益：从主网购电是正的，向主网售电是负的；
• 燃气轮机的燃料成本：基于天然气的等效耗量来算；
• 运维成本：P2G设备、CCS设备、储能设备的单位运维成本乘运行功率；
• 碳交易成本：由阶梯碳价函数计算。

用公式表达大概是：

[ \min \sum_{t=1}^{T} \left[ c_{buy}^t P_{buy}^t - c_{sell}^t P_{sell}^t + c_{gas} V_{gas}^t + k_{om} P_{om}^t + C_{carbon}^t \right] ]

在Matlab里，我先把每个设备的逐时功率设成sdpvar变量，然后将所有成本项累加成Objective对象，交给YALMIP处理。这种写法的好处是逻辑直观，每一类成本对应一段代码，排查错误时很方便。

这里有一个值得留意的细节：碳交易成本本身是分段线性函数，把它放到目标函数里时，我们不能直接写“if 排放>配额 then 碳价=80 else 碳价=120”这种判断语句，必须把分段函数转换成一组线性约束和0-1变量。下一节专门讲这个线性化处理。

2.2 阶梯碳交易成本的分段线性化

把阶梯碳交易成本写成数学表达式，设总CO2排放量为E，免费碳配额为E0，阶梯碳价为c1、c2、c3，对应的超额区间是[0, Δ1]、[Δ1, Δ2]、[Δ2, ∞)。那么碳交易成本F(E)可以定义为：

[ F(E) = \begin{cases} c_1 (E - E_0), & E_0 \le E \le E_0 + \Delta_1 \ c_1 \Delta_1 + c_2 (E - E_0 - \Delta_1), & E_0 + \Delta_1 \le E \le E_0 + \Delta_2 \ c_1 \Delta_1 + c_2 \Delta_2 + c_3 (E - E_0 - \Delta_1 - \Delta_2), & E \ge E_0 + \Delta_1 + \Delta_2 \end{cases} ]

因为F(E)是凸分段函数，可以直接用YALMIP自带的分段线性化工具yalmip('clear')配合pwa函数，或者手动引入三个0-1变量来建模。我习惯手动建模，因为这样生成的线性约束更精简，求解速度更快。

手动建模的思路是：把E-E0的超额排放量拆解为三段Δ1、Δ2、Δ3，每段的值域受0-1变量控制，保证低档区间先被填满。用一组大M约束就能实现。具体写出来是：

% 排放超额量拆成三段 e1 = sdpvar(1, 1); % 第一档超额量 [0, d1] e2 = sdpvar(1, 1); % 第二档超额量 [0, d2] e3 = sdpvar(1, 1); % 第三档超额量 [0, inf] % 大M常量，取一个明显大于所有可能排放量的数 M = 5000; % 0-1变量表示每档是否被“占用” z1 = binvar(1, 1); z2 = binvar(1, 1); % 第一档：要么整体不超过d1，要么等于d1 Constraints = [Constraints, e1 >= 0, e1 <= d1]; % 第二档：只有第一档满时才启用 Constraints = [Constraints, 0 <= e2 <= d2*z1]; % 第三档：只有第二档满时才启用 Constraints = [Constraints, 0 <= e3 <= M*z2]; % 触发条件 Constraints = [Constraints, e1 >= d1*z2 - M*(1-z2)]; Constraints = [Constraints, e1 >= d1*z1 - M*(1-z1)]; Constraints = [Constraints, e2 >= d2*z2 - M*(1-z2)]; % 总超额量 = e1 + e2 + e3 Constraints = [Constraints, total_carbon_over = e1 + e2 + e3];

这里面最关键的是确定大M的取值。M设得太大，求解器的数值稳定性会变差；设得太小，有可能把可行域错误地剪掉。一般我取“预测最大总排放量”的2~3倍，对典型的24小时调度问题，5000吨左右足够了。

2.3 关键约束一：电功率平衡与储能

电平衡是所有调度模型的地基。约束表达式是：

[ P_{buy}^t + P_{wt}^t + P_{pv}^t + P_{gt}^t + P_{dis}^t = P_{sell}^t + P_{load}^t + P_{p2g}^t + P_{ccs}^t + P_{ch}^t ]

其中P2G和CCS设备在这里都是耗电设备，储能充电是负荷，放电是电源。燃气轮机的出力是决策变量，风电、光伏、负荷在确定性模型中是预测值，在鲁棒模型里可以是区间变量。

储能约束按标准写法加：SOC递推方程、充放电功率上限、SOC上下限、爬坡限制。注意充放电状态不能同时为1，这个需要加互斥约束，否则求解器会为了“既要充电又要放电”而滥用二进制变量。

在Matlab里，储能SOC状态用一个24维的sdpvar连续向量表示，递推关系写成：

for t = 2:T Constraints = [Constraints, SOC(t) == SOC(t-1) - P_dis(t)/eta_dis*dt + P_ch(t)*eta_ch*dt]; end

这里dt是调度时段时长，取1小时。纯电储能和氢储能我都处理成类似结构，只是氢储能还要考虑向燃气轮机供氢的路径。

2.4 关键约束二：气网平衡与氢平衡

气网平衡是建模最容易出错的地方，因为它涉及“等效天然气流量”的概念。VPP用气的来源有三个：外购天然气、P2G甲烷化生成的甲烷、储气罐释放的天然气。用气去向有两个：燃气轮机消耗的天然气能量、外部气负荷。

我把所有用气量都用“等效纯天然气体积”来计量，这样掺氢的混合气流量才能和纯天然气流量做加减运算。具体写法是：

[ V_{buy}^t + V_{meth}^t + V_{tank,out}^t = V_{gt}^t + V_{load}^t ]

燃气轮机的等效天然气耗量由出力、效率、燃料热值共同决定。在Matlab里我用一个向量V_gt表示逐时耗气量，用等式约束关联到燃气轮机出力：

Constraints = [Constraints, V_gt(t) == P_gt(t) / eta_gt(t) / LHV_gas];

掺氢比例约束则写为：氢气体积流量 ≤ 掺氢比例上限 × 混合气总体积流量。

2.5 关键约束三：P2G-CCS耦合的物料平衡

这里我把装置之间的耦合逻辑展开讲。P2G单元用电生产氢气H2_prod，这部分氢气有三个去向：

• 直接并入燃气管网用于掺氢燃烧；
• 输送到加氢站对外销售；
• 进入甲烷化反应器，和CCS捕集的CO2反应生成甲烷。

甲烷化反应的简化计量关系是：每1摩尔CO2消耗4摩尔H2，生成1摩尔CH4和2摩尔H2O。实际计算时我把它折算成质量或体积比例关系：生成1立方米甲烷需要约4立方米氢气和1立方米二氧化碳。这个比例参数co2_to_ch4、h2_to_ch4直接影响模型里气网平衡与碳平衡。

CCS捕集的CO2来源于燃气轮机燃烧排放。燃气轮机的总碳排放量，先乘以捕集率η_ccs得到捕集量，捕集量中一部分送入甲烷化单元参与反应，剩余部分封存。因此碳交易成本的“净排放”等于总排放量减捕集量再减去P2G甲烷化消耗掉的CO2，最后和配额比较算阶梯价。

我当时在这个环节优化了一个做法：把P2G和CCS耦合引起的CO2自循环做成一个“节点”，要求进节点的CO2流量等于出节点的CO2流量，这让模型的可解释性变强了。代码长这样：

% 燃气轮机CO2排放量 CO2_gt = P_gt * emission_factor; % 捕集量 CO2_captured = CO2_gt * eta_ccs; % 捕集量去向：甲烷化 + 封存 Constraints = [Constraints, CO2_captured == CO2_meth + CO2_stored]; % 甲烷化消耗CO2量 = 甲烷化产天然气量 * co2_to_ch4 Constraints = [Constraints, CO2_meth == V_meth * co2_to_ch4]; % 净排放量 Constraints = [Constraints, CO2_net == CO2_gt - CO2_meth - CO2_stored];

2.6 关键约束四：燃气轮机的掺氢运行区间

燃气轮机的运行受最小技术出力限制。掺氢比例变了以后，混合气的热值下降，同样的流量对应的能量输出发生变化。为了简化，我假设燃气轮机按照“等效天然气能量输入”来发电，那么实际混合气的体积流量可以用等效天然气耗量除以（1-掺氢体积比）来推算。

此外，掺氢比例越高，燃气轮机可调范围越窄，我在模型中把出力上限做了一个随掺氢比例下降的修正，即：

[ P_{gt}^{max} \rightarrow P_{gt}^{max} \times (1 - k_{derate} \times \phi_{h2}) ]

这里的k_derate取0.05~0.1之间，表示掺氢比例达到上限时燃气轮机最大出力会降低5%~10%。这个修正不是标准教科书内容，是我在实际调试中发现掺氢比例太高时燃烧稳定性确实会变差才加进去的，虽然不是严谨的热力学模型，但对调度结果更接近工程实际。

3. Matlab代码实现：从建模到求解

3.1 为什么用YALMIP而不是手写线性规划矩阵

如果你只用Matlab自带的linprog或intlinprog，那么写约束矩阵会非常痛苦。100多个决策变量、几十个交叉耦合约束，手写矩阵等于自找麻烦。

YALMIP作为一个建模层，最大的好处就是符号化建模、自动转换成求解器需要的标准形式。你可以直接用sdpvar创建变量，用加法减法写目标函数和约束，然后调用optimize函数指定求解器。兼容性也很关键：同一个模型既能用CPLEX跑MILP，也能无痛切换到Gurobi，方便做求解器对比和解耦验证。

我的求解环境是Matlab R2023a + YALMIP R20230622 + CPLEX 12.10。Gurobi我也试过，速度比CPLEX快大概15%~20%，但CPLEX在Academic环境里更常见，所以我最终用CPLEX出结果。

3.2 数据准备：典型日负荷曲线和分时电价

调度周期取24小时，单位时间间隔1小时。模型输入数据包括：风电出力预测曲线、光伏出力预测曲线、电负荷预测曲线、气负荷预测曲线、分时购售电价、碳配额总量、燃机参数、P2G参数、储能参数。

下表是我用来复现论文结果的一组典型数据（部分），注意这只是示意，具体数值要跟你的研究对象绑定：

把数据写成Matlab的列向量后，模型代码就不需要再引用具体数值了，这样要切换算例，只需要改数据文件的数值即可。

风力光伏的波动性我是直接采用预测值，不考虑预测误差的场景树。这样做的好处是模型只验证调度思路，不掺入随机优化的额外复杂度。如果你要做鲁棒优化，可以把这个基础模型当成确定性子问题去嵌入。

3.3 核心代码框架：YALMIP变量定义与约束组装

代码结构和变量定义直接关系到可读性和调试效率。我习惯把一个大型调度模型的代码拆成4个区块：参数定义、变量定义、约束组装、目标函数与求解。

%% 1. 参数定义 T = 24; dt = 1; % 风电、光伏、负荷数据 P_wt = [18.5; 16.2; 15.1; ...]; P_pv = [0; 0; 0; ...]; P_load = [39.2; 35.8; 33.4; ...]; % 分时电价 c_buy = [320; 320; ...]; c_sell = [260; 260; ...]; % 设备参数 eta_p2g = 0.7; % P2G电解效率 eta_meth = 0.6; % 甲烷化效率 eta_ccs = 0.85; % 碳捕集率 eta_gt = 0.45; % 燃气轮机发电效率 h2_blend_max = 0.20; % 最大掺氢体积比 % 碳交易 E0 = 100; % 免费配额(吨) c1 = 80; c2 = 120; c3 = 200; % 阶梯碳价(元/吨) delta1 = 50; delta2 = 100; % 阶梯分档区间 M = 5000; % 大M常数 %% 2. 变量定义 P_buy = sdpvar(T, 1); % 从主网购电 P_sell = sdpvar(T, 1); % 向主网售电 P_gt = sdpvar(T, 1); % 燃气轮机出力 V_gt = sdpvar(T, 1); % 燃气轮机等效天然气耗量 V_buy = sdpvar(T, 1); % 外购天然气量 V_meth = sdpvar(T, 1); % P2G甲烷化产气量 P_p2g = sdpvar(T, 1); % P2G耗电功率 CO2_gt = sdpvar(T, 1); % 燃气轮机排放CO2量 CO2_captured = sdpvar(T, 1); % 捕集的CO2量 CO2_meth = sdpvar(T, 1); % 用于甲烷化的CO2量 CO2_stored = sdpvar(T, 1); % 封存的CO2量 CO2_net = sdpvar(T, 1); % 净排放量 % 储能变量 SOC_bess = sdpvar(T+1, 1); P_ch = sdpvar(T, 1); P_dis = sdpvar(T, 1); u_ch = binvar(T, 1); u_dis = binvar(T, 1); % 掺氢比例变量 phi_h2 = sdpvar(T, 1); % 逐时掺氢体积比 % 碳交易分段变量 e1 = sdpvar(1, 1); e2 = sdpvar(1, 1); e3 = sdpvar(1, 1); z1 = binvar(1, 1); z2 = binvar(1, 1);

变量定义阶段一定要把边界和单位写清楚。我见过不少人在sdpvar里忘记限定非负，结果求解器把购电量“优化”成负数，这种错误排查起来特别痛苦。所以连续变量定义后立刻加上非负约束，养成习惯。

%% 3. 约束组装 Constraints = []; % 非负约束 Constraints = [Constraints, P_buy >= 0, P_sell >= 0, P_gt >= 0, ...]; % 电平衡 for t = 1:T Constraints = [Constraints, P_buy(t) + P_wt(t) + P_pv(t) + P_gt(t) + P_dis(t) == ... P_sell(t) + P_load(t) + P_p2g(t) + P_ch(t)]; end % P2G耗电与产氢的关系 H2_prod = eta_p2g * P_p2g / LHV_h2; % 甲烷化产气量与耗氢量、耗CO2量的关系 V_meth = eta_meth * H2_to_meth / h2_LHV_meth; Constraints = [Constraints, CO2_meth == V_meth * co2_to_ch4]; Constraints = [Constraints, H2_to_meth == V_meth * h2_to_ch4]; % 气网平衡 Constraints = [Constraints, V_buy + V_meth == V_gt + V_load]; % 燃气轮机CO2排放 Constraints = [Constraints, CO2_gt == P_gt * emission_factor]; % 碳捕集 Constraints = [Constraints, CO2_captured == CO2_gt * eta_ccs]; % 捕集量分配 Constraints = [Constraints, CO2_captured == CO2_meth + CO2_stored]; % 净排放 Constraints = [Constraints, CO2_net == CO2_gt - CO2_stored - CO2_meth]; % 掺氢比例约束 Constraints = [Constraints, phi_h2 <= h2_blend_max]; % 掺氢量与V_gt的关联（按体积比） Constraints = [Constraints, phi_h2 == V_h2 / (V_gt + V_h2)]; % 修正后的燃气轮机出力上限 Constraints = [Constraints, P_gt <= P_gt_max * (1 - 0.1*phi_h2)];

3.4 约束组装这个环节，最容易出Bug

这里我要特别提醒：sdpvar变量的维度、for循环的索引、约束的拼接方式，只要有一个细节不对，YALMIP就会报维度不匹配或者逻辑错误。

常见的报错分三类：

• “Nonlinear constraints are not supported”——这说明你把逻辑运算符“或”“且”直接写进了约束，或者用了if判断；
• “Unbounded objective”——目标函数里漏了某个成本项，导致优化器认为某些变量取无穷大也能让成本最小，检查是不是有变量忘了进目标函数；
• “Infeasible problem”——约束本身冲突，常见于储能互斥约束或碳捕集量与排放量之间的耦合矛盾。

调试时我通常的做法是：先把所有约束注释掉，只保留目标函数和几个关键变量，然后迭代打开约束块。每打开一大块就求解一次，看是哪个环节把可行域剪空的。

3.5 求解设置与结果展示

求解代码很简单，核心是设置求解器参数，特别是MIP gap和求解时限：

%% 4. 目标函数 % 购电成本 objective = sum(c_buy .* P_buy) - sum(c_sell .* P_sell); % 燃料成本 objective = objective + c_gas * sum(V_buy + V_meth); % 运维成本 objective = objective + k_p2g * sum(P_p2g) + k_ccs * sum(CO2_captured) + k_bess * (sum(P_ch) + sum(P_dis)); % 碳交易成本 objective = objective + c1*e1 + c2*e2 + c3*e3; %% 5. 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'cplex.mip.tolerances.mipgap', 1e-4); diagnosis = optimize(Constraints, objective, ops); if diagnosis.problem == 0 disp('求解成功'); else disp(diagnosis.info); end

求解成功之后，我喜欢一次性把结果全部画出来：

figure; subplot(3, 1, 1); bar(1:T, [value(P_wt), value(P_pv), value(P_gt), value(P_buy)], 'stacked'); legend('风电', '光伏', '燃气轮机', '购电'); % 继续画电负荷、P2G耗电等

注意，value()函数把YALMIP变量转换为数值，这是出图前必须做的一步。有些人直接用sdpvar变量画图，结果Matlab报错，就是这个原因。

结果曲线上重点看三点：

• 燃气轮机的出力是否跟着碳价波动走。碳价走高时，燃气轮机大概率会降低出力，把负荷让给风电光伏；
• P2G设备的运行时段是否集中在电价低的时段。电价低意味着制氢成本低，P2G会倾向于半夜或午间充电制氢；
• 掺氢比例是否突破了上限。如果某个时段phi_h2卡在上限，说明该时段氢气供给充足，甲烷化反应器吃掉了大部分CO2、净排放量降下来了。

3.6 代码复现时要注意的细节

• YALMIP里sdpvar(T,1)创建的是列向量，如果你用sdpvar(1,T)，后面所有涉及sum和矩阵乘法的表达式都要对应调整，否则维度不匹配报错够你调一个小时；
• CPLEX求解后要检查diagnosis.problem字段，它等于0才代表找到最优解，等于1时不收敛或不可行，等于2时数据问题。不要只看优化器里打印的“Optimal solution found”就认为成功，因为那是CPLEX的输出，和YALMIP返回的诊断信息要分开看；
• 如果你的模型规模大，可以加上'cplex.mip.strategy.search', 1强制深度优先搜索，对含大量0-1变量的机组组合问题有奇效；
• 我建议开写代码之前，先在纸面上把设备和系统框架图画出来，标明每条能量流的流向和装置之间的耦合关系。不然直接上手写把两个无关变量的约束连在一起，排查逻辑错误时会疯掉。

4. 典型算例分析与对比：不同掺氢比下的调度策略

4.1 算例A：不掺氢，纯天然气 + CCS

我把这个作为基准算例。此时P2G产生的氢气全部卖给加氢站或者直接储存，甲烷化反应没有足够的CO2输入，基本上只有少量CO2封存。结果来看，总碳排放较高，碳交易成本占总成本的比例最大，达到18%左右。这在数学上很好理解，因为燃气轮机烧的是纯天然气，排放CO2全计入碳账，P2G生成的甲烷占比小，自循环碳量少。

4.2 算例B：掺氢比例上限10%

当掺氢比例设成10%以后，P2G的一部分氢气进入燃气管网，燃气轮机燃烧的混合气里氢气的体积比上升，等效天然气量下降。效果是外购天然气量下降约3%，碳排放量下降约4%，碳交易成本下降明显。

但注意：不是所有时段都适合掺氢。在晚间风电出力高峰且电价低廉时，P2G制氢成本低，产出的氢气可以有效掺入，燃气轮机会适当提高出力，因为在碳账上“烧氢”这个过程不产生碳配额成本，比外购电更划算。而在日间光伏出力大、电价高的时段，燃气轮机不太出力，掺氢需求也低，P2G更倾向于停机或只给加氢站供氢。

4.3 算例C：掺氢比例上限20%

掺氢比上限提升到20%之后，模型对P2G产能的利用率明显提高，外购气量和碳交易成本继续下降。但此时我设置的燃气轮机出力修正系数k_derate=0.1开始起作用，在某些高负荷时段，燃气轮机因出力上限下调不得不减少供电，缺口需要靠购电填补。

这就是掺氢高比例下的一个经济学权衡：省下的碳交易成本和燃料成本，与增加的购电成本和由于出力受限导致的备用损失，谁大谁小。我实际跑出来的结果是：在给定的典型日数据下，20%掺氢比的总成本反而比10%略高一点，原因就是高峰期燃气轮机的出力限制带来了额外的购电费用。

这个结果说明，掺氢不是越多越好，要结合燃气轮机运行区间、碳价水平和电价结构综合决定。调度模型的魅力就在于此：它不会给你一个拍脑袋的“掺氢越多越环保”，而是帮你算出每个掺氢比例对应的最优运行方式和经济性边界。

4.4 一个值得深挖的现象：阶梯碳价的“跳档”行为

用阶梯碳价时，模型会出现一种有意思的行为：某些方案的总排放量会恰好卡在配额附近，或者卡在某一档的边缘。这是因为跨过档位的边际碳价会突然上涨，导致优化器宁可增加运行成本也不让排放量越过档位阈值。

我在结果里就观察到，算例B的总排放量是94.6吨，刚好压在配额95吨附近。你再逼近一步就会发现，如果把碳配额从100改成94吨，那么最优解会立刻发生跳跃，燃气轮机减小出力、P2G加足马力多制氢，用成本更高的方法把排放压下车。这种非连续的响应，就是阶梯碳价设计者想要的效果。

5. 常见问题与排查技巧

5.1 求解器报“Infeasible problem”怎么办

先检查储能相关约束，尤其是SOC的初始值和终值。我经常忘记设置SOC终端约束，导致储能系统在调度周期结束时还保留大量电量，这从优化角度看是可用的，但如果你人为要求SOC_T=SOC_0，就可能在风机光伏出力低时出现无解。

再检查互斥约束的写法。我见过有人用u_ch + u_dis <= 1约束，但因为u_ch和u_dis分别在充电、放电约束里被>=强制成1，导致两者永远不可能同时满足。逻辑上应该允许u=0时充放电都为0，用大M约束放开就行。

5.2 P2G-CCS耦合部分总是报非线性错误

把phi_h2 == V_h2 / (V_gt + V_h2)这种分式约束原样写进YALMIP，大概率报非线性错误。正确做法是把分母乘到等式另一边：

[ \phi_{h2} \cdot (V_{gt} + V_{h2}) = V_{h2} ]

这仍然是个双线性项，因为phi_h2 * V_gt是变量乘积。我最终用线性化方案替代：引入辅助变量V_h2_eq = phi_h2 * (V_gt + V_h2)，然后加约束V_h2_eq == V_h2，麻烦的是这个等式本身也是双线性的。

实际上，更工程化的做法是：不把掺氢比例设成逐时自由变量，而设成固定参数。比如直接设phi_h2=0.1或者phi_h2=0.2，这样约束变成一个线性等式，求解器毫无压力。如果你想做掺氢比例优化，可以用0-1变量离散化几个档位：10%、15%、20%，再通过大M约束选出最优点。这是工程上更常规的处理思路，学术上也可以用分段线性化和二进制展开去逼近连续变量，但代价是求解复杂度急剧上升。

5.3 结果里出现“购电和售电同时为正”

如果目标函数里购电价大于售电价，理论上优化器绝不会同时购和售。如果结果出现了，说明某个约束把购电或售电变量固定住了。检查是不是有“最大购电量不超过负荷”之类的硬约束，或者因为电平衡表达式里变量符号写反，导致优化器用购电加售电来“平衡”多余变量。

5.4 数值稳定性问题

大M取值过大，或者P_gt上限设到1000MW而没有实际量纲，就容易导致求解器报“Numerical issues”。解决办法：一是把所有变量量纲统一到MW、MWh、吨，别有的用kW、有的用W；二是给每个变量设定合理的上下界，避免变量取值范围跨好几个数量级；三是如果非要用大M，尽量让M只比对应变量上限大一点，别为了省事给个1e9。

5.5 和文献结果对不上，误差在哪

多数对不上的情况出在碳交易成本的分段计算上。文献里阶梯碳价可能有两种口径：一种是“超额部分全部按最高档价格计算”，另一种是“分档累计累进”。这两者的碳交易成本差很多，我在做结果对比时就吃过亏，一开始按第一种口径算，后来发现总成本比文献结果高出一截，换成第二种才吻合。

另一个常见差异是免费碳配额的分配对象。有的文献把配额分给燃气轮机，有的分给虚拟电厂整体排放，还有的按机组逐台分配。对照公式时，你一定得确认文献里的排放源边界覆盖了哪些设备。风电光伏没排放，但P2G耗电如果是从主网购的，火电厂的间接排放算不算进配额，不同文献口径差异很大。建议单独用表格记录每个文献的配额口径和阶梯价格设计，然后再决定自己的模型采用哪一种。

6. 实际项目中的一点体会

整个项目做下来，我最大的感触是：调度模型难的不在求解算法，而在把“碳-氢-气-电”这几种能量流的物理逻辑翻译成数学约束的过程。每个等式都要问一句：它真的描述了物理世界里的平衡吗？

比如燃气掺氢，工程上和模型里的“掺氢比例上限”不是同一个概念。工程上，20%是安全极限，但在实际运行里你还要考虑管道末端的组分波动、燃烧器的燃料适应性、甚至环境温度对火焰温度的影响。我在模型里用出力修正系数来模拟这种限制，虽然比不上CFD模拟精确，但足以让调度结果比“理想化掺氢无限额”贴近现实。

另外，阶梯碳交易机制对调度的影响远超我最初预期。它不像固定碳价那样只拉高成本，而是改变了最优解的结构。不同档位之间的边际碳价差异，导致燃气轮机会在某一个排放水平附近出现“驻足不前”的调度行为。这种非连续的策略切换，恰恰是双碳政策设计者想要的激励效果。

代码我后来整理成了一个带GUI的版本，输入风电光伏负荷曲线，选择掺氢上限和碳交易参数，点一下就能跑出24小时的机组出力和碳交易成本报表。有这个工具之后，领导要算“掺氢到20%一年能省多少碳配额”这种问题，都是几分钟出结果。后续我还打算把鲁棒优化和分布鲁棒优化加进去，替代现有的确定性模型，应对风电光伏的不确定性。如果你在复现这个模型时遇到具体的报错或对不上的结果，建议先从碳交易分段函数的数学表述查起，那是最隐蔽也最影响结果的地方。