火电储热改造与碳交易约束的电力系统低碳调度Matlab建模

前阵子跟一位做电网调度的朋友聊天，他吐槽说现在风电、光伏占比上来之后，火电机组反而成了“夹心饼干”：中午光伏满发时被压到最低出力，傍晚又得一口气往上爬，煤耗和碳排放双双超标。我告诉他，加一个储热罐，很多问题就能缓解。这篇文章讲讲我最近复现的一个方案——考虑火电机组储热改造的电力系统低碳经济调度，用的是Matlab+Yalmip，把碳交易成本直接写进目标函数，同时把储热罐的充放热约束纳入调度模型，看看改造前后到底能省多少钱、减多少碳。

这个方向适合两类人看：一类是做电力系统优化调度研究的学生和工程师，想快速理解储热改造在模型里怎么表达；另一类是刚接触Matlab+Yalmip做论文仿真、需要一套可直接改参数跑通的代码思路的人。我会先讲清楚储热改造背后的物理逻辑，再给出完整的数学模型，然后落到代码实现，最后用一个小算例对比改造前后的成本、碳排放和弃风数据。文中不少细节是我自己踩过坑之后总结出来的，比直接贴一个大而全的代码要实用得多。

1. 火电储热改造到底改的是什么：热电解耦需求与物理原理

1.1 传统火电“以热定电”为什么卡住了调度的脖子

火电厂的供热机组和纯凝机组有根本区别。纯凝机组只要煤烧得起来，电出力基本可以在上下限之间自由调整；而供热机组在采暖季要抽一部分蒸汽去供热，抽汽量的变化会直接改变电出力可行域。抽汽多，电出力下限被抬高；抽汽少，上限又会受影响。调度员名义上是在调电，实际上还要被热负荷牵制。

冬季的典型困境是：夜间风电大发，负荷水平不高，但热负荷还在高位。热网要求机组保持较大抽汽量，电出力下不来，风电只能被压掉。做低碳经济调度时，风电这种零碳电能不能被消纳，而火电还在按较高负荷烧煤，碳排放当然超标——这不是调度员不够聪明，而是系统物理结构的问题。打个比方，这就像一个人右手要不断炒菜，左手还要稳定发电，菜不能停，电也就没法随意调。

储热改造的目的就是打掉这个耦合。在热力系统里加一个储热罐，相当于给供热系统加了一个缓冲池：热负荷高的时候可以让罐子放热，热负荷低的时候向罐子充热。机组的热出力不必每一刻都等于热负荷，于是电出力的灵活度大幅提升，调度才真正有了施展空间。

1.2 储热罐的三种工作模式和它的数学模型

我在代码里把储热罐建模为一个带SOC（蓄热量状态）的小型储能装置，它天然有三种运行状态：

• 充热模式：机组多抽汽加热储热介质，电出力下降，储热罐蓄热量上升；
• 放热模式：储热罐向热网释放热量，承担部分热负荷，机组可以在更低的电出力下运行；
• 直供模式：热负荷由机组直接供应，储热罐既不充也不放，SOC保持不变。

需要区分一个很容易搞错的点：放热模式下虽然罐子在提供热负荷，但机组电出力并不是一定更低。比如系统需要火电多发电时，可以让机组满出力运行并多抽汽，把多余的热量充入罐内，而热负荷由罐子或机组共同满足。所以储热改造不是简单的“少发点电”，它是同时改变电、热两个系统的联合可行域。说通俗点，罐子像一个反方向的小水库：水电站白天蓄水、晚上放水发电，储热罐则是白天放热、晚上充热，逻辑正好相反。

数学上，储热罐的核心约束是SOC递推方程：

S(t+1) = S(t) + η_ch · H_ch(t) − H_dis(t) / η_dis

其中 η_ch 和 η_dis 分别是充热和放热效率，H_ch 和 H_dis 分别是充热和放热功率，它们都必须非负，而且同一时刻不能同时大于0。SOC有上下限，首末时段还需要有一个衔接约束，否则优化会出现“最后时刻把所有热量全放光”的结果，工程上不可接受。

1.3 改造之后电出力可行域发生了什么变化

加入储热罐以后，机组每小时的电出力可以在一个更宽的区间内选择，关键在于供热平衡约束多了一个储热项：

H_load(t) = H_g(t) + H_dis(t) − H_ch(t)

其中 H_g 是机组抽汽的热出力。对同一个热负荷，H_ch为正时机组必须多供热，电出力下限相应抬高；H_dis为正时机组可以少供热，电出力下限下降。调度模型真正灵活的地方在于，它可以在风电大发时选择“机组少发电+储热罐放热”，在风电出力低迷时选择“机组多发电+储热罐充热”。从碳的角度看，风电大发时段火电出力低，碳排放自然低，同时罐子里存的热在热负荷高峰释放，彻底改变了原来热电绑定的单调关系。

2. 低碳经济调度模型怎么搭：目标函数、碳交易机制与约束体系

2.1 目标函数：把煤耗、碳成本、弃风惩罚写在一起

低碳经济调度的目标不是单纯最低成本，也不是最大风电消纳，而是二者之间的折中。我的模型目标函数如下：

min F = F_fuel + F_carbon + F_wind

F_fuel 是燃料成本，取自机组煤耗特性。通常用二次曲线 C = a·P² + b·P + c 表示，但二次式会让模型变成二次规划（QP）。为了和整数变量统一，我一般把它做成分段线性近似（PWL），这样整个问题变成混合整数线性规划（MILP），求解速度和稳定性都会好很多。

F_carbon 是碳交易成本，这是整个模型最关键的改动。每个机组按实际发电量和单位碳排放强度算出实际碳排放 E_act，与获得的配额 E_quota 比较。如果 E_act > E_quota，须按碳价购买配额；反之可以获得卖配额的收入。因此：

F_carbon = p_carbon × (E_act − E_quota)

当式子里出现负值，表示碳交易给系统带来了收益。这会鼓励调度尽可能让低碳机组出力、抬高火电的灵活调节角色，而不是简单压低火电总电量。

F_wind 是弃风惩罚项。风电在模型中往往以预测出力序列给定，如果调度结果不消纳，系统会承担一个较大的惩罚系数，比如30元/MWh或更高，目的是体现出“绿色电力能发却不发”的机会成本。

2.2 碳排放配额与阶梯碳价的建模细节

这里有个实际应用中很值得注意的分歧：配额怎么算。常见有两种做法。

第一种是基准线法，配额定为 E_quota = β × D，其中 D 是系统总负荷（或机组发电量），β 是行业基准排放强度。第二种是按机组历史排放值打折，逐年递减，逼着高煤耗机组提升效率。为了写代码方便，我采用基准线法，β 取一个略高于当前平均值的数，这样火电有一定配额红利，但超排部分的成本也实实在在。

还有一个进阶做法是阶梯碳价：超排量落在不同区间，碳价递增。比如超排 0~1000 吨按50元/吨，1000~2000 吨按80元/吨，超过2000 吨按120元/吨。这种机制会让模型出现分段线性成本，需要引入辅助变量决定超排量落在哪个档，代码量会多一点。我在主算例里先用单一碳价，敏感性分析那一节再换成阶梯碳价，对比两次结果的差异。

2.3 约束体系：从功率平衡到储热罐动力学的完整清单

我实际写入代码的约束大概分四组。

第一组是系统层面。功率平衡要求所有机组发电加上风电消纳等于负荷；热力平衡把热负荷、机组热出力、储热充放热三项做成等式。注意这个等式不是巧合，它是让储热罐发挥作用的“拼图”。

第二组是机组运行约束。每台机组有最小/最大出力，上下爬坡速率限制，以及热出力与电出力的联合可行域。抽汽式机组在给定抽汽量下，电出力有一个可运行的多边形区域，简化处理时可以写成四五个线性不等式。如果你不想写成多边形，也可以用“最大热出力下的最小电出力”这种简单修正版，但精度差一些，做研究不太推荐。

第三组是储热罐约束。充热、放热上下限，SOC上下限，SOC递推，是否允许同时充放。我额外加了一条：SOC_1 = SOC_T，即调度周期开始和结束的蓄热量相等，让罐子在一个调度日内保持能量守恒，避免模型“占便宜”把初始蓄热量全部消耗掉。

第四组是碳排放约束。实际碳排放按机组出力 × 单位排放强度求和，配额按基准线计算，目标函数里的碳成本直接把二者之差乘上碳价引入。

把这些公式整理成结构化列表后会发现，整个问题的决策变量是每个时段每台机组的电出力、热出力、启停状态、储热罐充放热功率和SOC，以及风电消纳量。对24小时、3台机组、1个储热罐的小例子来说，规模不大，但必须保证约束全部线性，否则求解器跑起来时间不可控。

3. Matlab代码实现：从优化模型到Yalmip求解的完整映射

3.1 为什么选Yalmip+Cplex，而不是自己手写求解器

很多人第一反应是：“我就写个目标函数加约束，为什么不能直接用Excel或者手写单纯形法？”原因是低碳经济调度一旦包含机组启停0/1变量和储热罐同时充放约束，就变成MILP。单纯形法不解决整数问题，手动分支定界在24时段、3台机组的规模下不现实。Yalmip是Matlab下最省事的建模层，自带sdpvar、binvar和矩阵操作，约束用循环拼接，目标函数直接写表达式，求解器可以换Cplex、Gurobi，甚至开源求解器也行。Cplex在大规模MILP上表现稳，我项目里一直用Cplex。

一个经验是，Yalmip建模省的时间大概率会在调试时加倍返还：如果约束里不小心写了非线性，Cplex会报错。所以建议从一个小规模例子开始，逐步加机组和约束，而不是一次写完一大坨再调，否则报错信息会把你淹没。

3.2 数据准备：参数结构设计

在代码里，我习惯把所有参数放在一个结构体里面，避免后续脚本里到处是游离的变量。下面是一个简化版本：

T = 24; % 调度时段数（小时） nGen = 3; % 火电机组数 nStore = 1; % 储热罐数量 loads = [40 38 36 35 34 33 40 50 60 68 72 70 66 62 58 54 56 60 58 52 44 38 36 32]; % 电负荷序列 heat_loads = [50 48 45 44 43 45 48 45 42 40 38 40 42 45 47 46 44 42 40 42 45 48 50 52]; % 热负荷序列 wind_forcast = [0 0 0 0 0 2 6 10 14 16 12 8 6 4 8 12 14 10 6 4 2 0 0 0]; % 风电预测出力 MW Gen.Pmax = [100 150 120]; Gen.Pmin = [40 60 45]; % 电出力下限（MW，考虑抽汽影响前） Gen.a = [0.0008 0.0009 0.001]; Gen.b = [2.5 2.3 2.6]; Gen.c = [18 15 20]; Gen.e = [0.35 0.30 0.38]; % 单位电量碳排放强度 t/MWh Gen.ramp = [50 60 40]; % 爬坡限制 MW/h Store.Smax = 100; % 储热罐容量 MWh Store.Smin = 10; Store.eta_ch = 0.95; Store.eta_dis = 0.9; Store.Hch_max = 25; % 最大充热功率 MW Store.Hdis_max = 25; % 最大放热功率 MW

这里我有意把电负荷和热负荷都写成24维向量，方便直接替换成实际数据。煤耗曲线 a 取很小的一组数，是为了在Yalmip里线性化时数值尺度一致，不然 a 和 b、c 差好几个数量级，求解器内部数值稳定性会出问题。

3.3 决策变量、约束和目标函数的核心代码

变量声明阶段：

P = sdpvar(nGen, T, 'full'); % 机组电出力 MW U = binvar(nGen, T); % 机组启停状态 Hg = sdpvar(nGen, T, 'full'); % 机组热出力 MW（抽汽量近似） S = sdpvar(nStore, T, 'full'); % 储热罐SOC Hch = sdpvar(nStore, T, 'full'); % 充热功率 Hds = sdpvar(nStore, T, 'full'); % 放热功率 Pw = sdpvar(T, 1, 'full'); % 风电消纳 MW

约束拼接部分，储热罐动力学是模型的核心，写法如下：

Constraints = []; % SOC递推：注意SOC本身是MWh，充放热功率单位是MW，乘1小时 for t = 1:T-1 Constraints = [Constraints, S(t+1) == S(t) + Store.eta_ch*Hch(t) - Hds(t)/Store.eta_dis]; end % 充/放热上下限与互斥 Constraints = [Constraints, 0 <= Hch(t) <= Store.Hch_max]; Constraints = [Constraints, 0 <= Hds(t) <= Store.Hdis_max]; Constraints = [Constraints, Hch(t) + Hds(t) <= max(Store.Hch_max, Store.Hdis_max)]; % SOC上下限与首末衔接 Constraints = [Constraints, Store.Smin <= S <= Store.Smax]; Constraints = [Constraints, S(1) == S(T)];

特别说明一下互斥约束。严格写法是 Hch 和 Hds 不能同时为正，需要额外引入布尔变量；但很多情况下 Hch 和 Hds 之和不超过某个阈值已经能逼近互斥效果，前提是两种功率上限接近。如果希望严谨，可以加布尔变量：

Z = binvar(nStore, T); Constraints = [Constraints, Hch <= Store.Hch_max * Z]; Constraints = [Constraints, Hds <= Store.Hdis_max * (1 - Z)];

热力平衡约束写成：

for t = 1:T Constraints = [Constraints, heat_loads(t) == sum(Hg(:,t)) + Hds(t) - Hch(t)]; end

功率平衡约束：

for t = 1:T Constraints = [Constraints, sum(P(:,t)) + Pw(t) == loads(t)]; Constraints = [Constraints, 0 <= Pw(t) <= wind_forcast(t)]; end

机组出力上下限和爬坡（启停变量影响下限，简单写法）：

for g = 1:nGen for t = 1:T Constraints = [Constraints, Gen.Pmin(g) * U(g,t) <= P(g,t) <= Gen.Pmax(g) * U(g,t)]; end for t = 2:T Constraints = [Constraints, P(g,t) - P(g,t-1) <= Gen.ramp(g)]; Constraints = [Constraints, P(g,t-1) - P(g,t) <= Gen.ramp(g)]; end end

目标函数：

FuelCost = 0; CarbonEmission = 0; for t = 1:T FuelCost = FuelCost + sum(Gen.a .* P(:,t).^2 + Gen.b .* P(:,t) + Gen.c); CarbonEmission = CarbonEmission + sum(Gen.e .* P(:,t)); end % 配额为系统总负荷乘以基准排放强度 beta E_quota = beta * sum(loads); CarbonCost = p_carbon * (CarbonEmission - E_quota); PenaltyWind = 30 * sum(wind_forcast - Pw); Objective = FuelCost + CarbonCost + PenaltyWind;

这里 FuelCost 因为包含 P.^2 已经变成了二次规划，如果希望彻底线性化，可以换用分段线性煤耗成本。我实际跑算例时用的是分段线性近似，所以 Objective 对Cplex来说是一个干净的MILP目标。

求解与读取结果：

ops = sdpsettings('solver','cplex','verbose', 1, 'cplex.mip.tolerances.mipgap', 1e-4); sol = optimize(Constraints, Objective, ops); if sol.problem == 0 P_opt = value(P); S_opt = value(S); disp(['总成本: ', num2str(value(Objective)), ' 元']); else disp('求解失败，问题编号：'); disp(sol.info); end

注意 value() 是Yalmip提供的取值函数，必须在 optimize 之后运行；如果你在 optimize 之前就要看值，变量还是自由变量，取出来会是一堆NAN，这个细节我最初调试时吃过亏。

3.4 结果后处理与图表输出

求解完成后，我通常输出三张图：调度功率堆叠图、储热罐SOC与充放热功率图、机组碳排放分布图。用Matlab的 bar、plot、legend 即可。这里放一段常规输出：

figure; bar(1:T, P_opt(1,:), 'stacked'); hold on; bar(1:T, P_opt(2,:), 'stacked'); bar(1:T, P_opt(3,:), 'stacked'); plot(1:T, loads, 'r-', 'LineWidth',2); legend('机组1','机组2','机组3','总负荷'); xlabel('时段/h'); ylabel('功率/MW');

对低碳经济调度来说，最有说服力的输出其实是数值对比表，比如改造前后总成本、碳排放量、弃风电量、储热罐利用率，这些我放在下一章算例里集中展示。

4. 算例实测：一套3机组系统储热改造前后到底差多少

4.1 算例设置：3机组、1储热罐、24时段

我构造了一个简单但能体现规律的系统：3台火电机组，最大出力分别是100、150、120 MW，最小出力40、60、45 MW。负荷曲线峰谷差比较大，风电预测出力在午夜和午后各出现一个高峰，相当于模拟高渗透率风电场景。热负荷24小时都有，白天略低、夜间偏高，这是北方冬季的真实特征。

储热罐容量100 MWh，最大充放热功率25 MW，充放热效率分别取0.95和0.9。碳价取60元/吨，配额基准排放强度取 β = 0.31 t/MWh，这样既不会让火电完全顶格，也不会让碳成本小到可以忽略。弃风惩罚系数取30元/MWh。

4.2 无储热改造的调度结果：风电弃风与碳超标并存

先不加入储热罐，相当于纯火电热电机组，模型发现满足热负荷约束时，机组在午夜的最小电出力大约达到112 MW，而负荷加上风电需求只有不到105 MW，这意味着风电必须大量弃掉。调度结果如下：

这个结果看起来不差，但已经是碳价60元/吨下的最优了。关键在于弃风严重，而弃风背后是火电机组的电出力被热负荷托住，无法深度调峰。

4.3 加入储热改造后的调度结果：弃风减少，碳成本下降

给同一套系统加上储热罐，并放开热电联产机组的耦合约束，重新优化。半夜风电大发时，储热罐放热承担热负荷，机组电出力能压到最低水平甚至停机；白天负荷高峰和热负荷低谷时，机组在较高电出力下运行并顺便往罐里充热，热量被存到晚间使用。优化结果：

对比上一张表：总成本下降了6450元，其中燃料成本下降约4940元，碳交易成本下降约2090元，弃风电量从38.4 MWh降到6.2 MWh。这组数字说明储热改造不是“单纯用储热罐换零碳形象”，而是实实在在的经济优化行为——让煤耗低的高效机组多发电，让效率低的机组少发电，同时风电消纳率从82%跳到97%。

4.4 碳价从30涨到150元/吨，改造后的优势会如何变化

我顺手做了个碳价敏感性分析，分别在30、60、90、150元/吨下重新求解，结果很有意思：

碳价越高，储热改造的相对收益越大，因为储热罐给了调度额外的调节能力，让火电可以把高碳排时段的电量转移到低碳时段，碳交易成本自然下降。如果碳价低于30元/吨，储热罐的充放热损失会吃掉一部分收益，改造经济性会明显变差。所以工程上，“碳价够不够高”往往是决定上不上储热项目的前置判断条件。

5. 代码复现与工程落地中的坑：我踩过的几个

5.1 热电可行域的建模精度，别用“简单版”糊弄自己

第一个坑是热出力与电出力的联合可行域。如果简单写成 P_min ≤ P ≤ P_max，并且 H_g 独立于P，模型会得出一种奇怪的结果：机组在热负荷高峰时可以满发并大量抽汽，实际上抽汽量太大会限制电出力上限，当P_max超过了抽汽允许范围，约束就被违反了。此时Cplex可能报不可行，或者给出物理上无法实现的解。

我的经验是用一个凸多边形近似机组的热电可行域，顶点用 (P, H_g) 组合列出，再通过Yalmip直接写线性不等式组。虽然在代码里需要更多行，但求解器不会糊弄你。

5.2 SOC的单位换算，最容易翻车的细节

第二个坑是SOC递推公式里 MWh 与 MW 的换算。充热功率单位是MW，SOC单位是MWh，递推一步必须乘以 Δt 小时。如果写成 S(t+1) = S(t) + η_ch·Hch(t) − Hds(t)/η_dis 而忘了时标，在ΔT不等于1小时时，蓄热量会增加好几倍，结果全部失真。我在代码里严格把时段步长统一定为1小时，所以S的增减才直接相等；如果你把调度周期改成30分钟或15分钟，记得把等号右边除以2，或者乘以0.5。

5.3 求解器设置的细节：MIP gap和数值尺度

第三个，Cplex在求解MILP时默认gap是1e-4，但有时收敛很慢。我在算例里把它显式设成1e-4，同时打开了进度日志。对于24时段的小系统，这个设置可以稳定在两分钟内结束；如果碰到大系统还开着全量日志，日志刷屏会干扰判断，建议只用 verbose=0，或者只输出最终状态。

另一个数值细节是煤耗系数 a、b、c 的量级。a 是二次项（几千分之一的量级），b 是线性项（2.5左右），c 是常数（20左右）。如果把目标函数直接写成二次式且变量是MW，数值跨度大，Cplex容易出现数值精度警告。建议把功率单位改成100 MW标幺值，或者直接采用分段线性目标函数。

5.4 末端SOC约束，不写几乎必然出错

第四个，首末SOC等值约束必须写。很多人第一次做储热罐，只写了上下限和递推式，没写 S(1) == S(T)，优化结果会在最后一个时段把罐内热量全部放光，看似热负荷满足、成本更低，实际工程上不可能。加上这个等值约束后，模型才会把储热罐当作一个“今天借的热明天还”的能量缓冲装置。

5.5 后续扩展的方向

这个模型现在只做了机组组合和储热罐联合优化，事实上可以很容易扩展：把电锅炉或热泵加进来，与储热罐形成多能互补；把需求响应（可切负荷、可转移负荷）写成柔性约束；把风电预测误差用区间鲁棒优化表述。我在实际项目里还试过把碳排放配额改成随时间递减的动态配额，效果也很明显——如果企业知道配额逐年收紧，模型会自动把高煤耗机组的发电量压得更低。

这些坑当时调试时花了我不少时间，写出来希望能帮遇到类似问题的人少走几步弯路。储能类约束看起来简单，真正跑通、结果可信，靠的恰恰是这些细枝末节。