考虑火电机组储热改造的电力系统低碳经济调度

这几年做电力系统优化调度方向的课题，有一个绕不开的痛点：新能源占比越来越高，但火电依然承担着调峰和供热的主力角色。尤其是北方地区，热电联产机组占比极大，"以热定电"的刚性约束让机组在冬季深调能力非常有限，弃风弃光问题在供暖期集中爆发。储热改造就是针对这一矛盾的核心手段之一——在火电机组侧加装储热罐，把热负荷和电出力解耦，给机组腾出更多调节空间。而"低碳经济调度"则是在传统经济调度框架之上，进一步叠加碳配额、碳交易成本、碳流追踪等因素，把碳排放变成调度决策中的可量化成本，从而引导系统走向更低碳的运行方式。

这篇文章总结的是我用Matlab实现"考虑火电机组储热改造的电力系统低碳经济调度"的全过程，从模型拆解、约束构建、代码实现到调试排坑。面向的读者是正在做电力系统优化调度、储能/储热配置、低碳经济调度相关毕业设计或课题研究的学生和工程师。如果你手上正缺一套能跑通的思路框架，或者对Yalmip建模、Cplex求解MILP问题还不太熟，这篇文章应该能帮你省不少时间。

1. 项目核心思路拆解

1.1 储热改造到底在改什么

很多人第一次接触"储热改造"这个概念时，容易把它和电网侧的储能电站混淆。这里要先把物理对象搞清楚。

火电机组储热改造，改造对象是热电厂。传统热电联产机组（CHP）在采暖期必须保证供热出力，而供热出力又和电出力强耦合——典型的抽凝式机组，电出力范围被热出力限制在一个狭窄区间内，这就是所谓"以热定电"。你想多发电，热也得跟着多产；你想少发电，热又供不上。风电大发需要火电压出力时，热电机组往往压不下去，只能眼睁睁弃风。

加装储热罐之后，热负荷不再是"必须由锅炉产热即时满足"，而是可以拆成两部分：一部分直接由机组供热，另一部分由储热罐放热供给。反过来，机组在需要的时候也可以多产热，把热量存进罐里，等热负荷高峰再放出来。相当于给热力系统加了一个缓冲池，热和电的强耦合关系被打破了。

从数学上讲，改造前机组电出力P和热出力H之间存在一个刚性可行域，改造后可行域在储热罐SOC的调节下发生了变化。调度模型里需要考虑的不再是单个时段的电热耦合，而是一个跨时段的热量转移约束。这对模型的动态建模能力提出了要求，也是整个项目建模的核心难点之一。

1.2 低碳经济调度解决了什么问题

传统经济调度的目标很单纯：让系统总运行成本最低，主要就是煤耗成本。但"双碳"目标下，只算煤耗账是不够的。碳排放要进入成本函数，才会影响调度决策。

低碳经济调度做的其实就是两件事的折中：保证供电和供热的经济性，同时尽可能降低系统碳排放。常见的做法有几种：

• 引入碳交易机制，给每个机组分配碳排放配额，实际排放超过配额的部分需要购买配额，低于配额的部分可以出售获利；
• 将碳捕集设备的能耗和收益纳入调度模型；
• 在目标函数中加入碳税或碳价的显式项。

这个项目采用的是第一种思路，即基于碳配额与碳交易成本的建模方式。每个火电机组根据出力水平计算实际碳排放量，再与政府分配的配额对比，差值乘以碳价计入目标函数。配额分配方式常用基准线法，也就是基于机组类型和发电量给出一个排放基准。

这样的建模方式有一个明显优势：目标函数仍然是线性或者分段线性的，配合混合整数线性规划（MILP）求解器可以直接求解，不需要引入复杂的非线性求解器，计算效率和收敛性都有保障。

1.3 为什么用Matlab来做

我见过不少人在Python和Matlab之间反复纠结，客观讲两者都能做，但Matlab在这个项目中确实有不可替代的便利性。

首先是Yalmip工具箱。它提供了一个极其简洁的建模层，让你几乎不需要关心求解器底层的接口细节。写约束就是一个个直观的公式，比如Constraints = [Constraints, Pgen >= Pmin];这种写法，跟数学表达式几乎没有区别。对于做电力系统调度研究的同学来说，Yalmip几乎是最短路径。

其次是Cplex求解器的Matlab接口相当成熟，MILP求解能力非常稳定。电力系统调度问题本质上是带0-1变量的混合整数线性规划，机组启停、储热罐充放状态是最典型的0-1决策变量。Cplex处理这类问题的速度值得信赖，我用它跑过24时段、几十台机组的算例，完全不担心性能。

另外Matlab的可视化能力和矩阵运算能力也省了很多事。调度结果要画机组出力曲线、风电出力曲线、SOC曲线、碳排放量曲线，Matlab几行代码就搞定，数据分析和结果诊断的效率非常高。

2. 低碳经济调度模型的关键建模细节

2.1 目标函数的结构设计

目标函数是调度的"指挥棒"，写得好不好决定模型行为是否符合预期。这个项目的目标函数有四个关键部分。

第一项是煤耗成本。机组煤耗量和出力基本呈二次函数关系，但Cplex只能处理线性和二次规划，MILP要求目标函数是线性的，所以我们对二次煤耗曲线做分段线性化近似。实际操作中我常用F = a*P^2 + b*P + c的形式，将P的区间切分成若干段，每段用线性函数代替，分段越多精度越高，但变量也越多。10段以内精度已经完全够用。

第二项是碳交易成本。COST_carbon = P_carbon * (E_actual - E_allowance)，其中E_actual是机组实际碳排放量，E_allowance是分配的配额，P_carbon是碳价。如果实际排放低于配额，这项成本为负，相当于系统通过"卖碳"获得收益。碳价的取值对调度结果影响很大，我在实验中常用的基准范围是30~80元/吨，灵敏度分析时在这个区间内扫值，观察机组出力和风电消纳的变化趋势。

第三项是弃风惩罚成本。风电的边际成本几乎为零，调度中原则上应全额消纳，但实际受限于系统调峰能力和线路传输能力，必须弃掉一部分。给弃风设置惩罚成本，本质上是给"少弃风"一个优先级，让模型在多个可行方案中优先选择风电消纳更多的方案。惩罚系数的设置需要远高于煤耗和碳成本单价，否则模型会把弃风当作理性的"最优解"。

第四项是储热系统的运行成本，这部分相对较小，主要是储放热过程中的损耗、维护成本等。如果暂不考虑这些，也可以设为0，让模型自然地为了降低煤耗和碳排而使用储热。

2.2 储热系统建模

储热罐的建模是整个模型的核心，也是这个项目和普通经济调度模型最大的区别所在。

储热罐本质上是能量在时间维度上的移动装置。它不产生能量，只扮演"时间搬运工"的角色。模型里我用一个状态变量SOC(t)表示t时段结束时储热罐的储热量（单位MWh），它的动态方程非常简单：

SOC(t) = SOC(t-1) + eta_charge * Q_ch(t) - Q_dis(t) / eta_discharge - Q_loss(t)

其中Q_ch(t)是t时段的充热功率，即机组向储热罐注入的热功率；Q_dis(t)是放热功率，即储热罐向外输出供热的热功率；eta_charge和eta_discharge分别是充放热效率，一般取0.9左右。注意放热等式里除法方向不能搞反，因为放热时储热罐输出给用户的是有效热量，罐内实际消耗的热量会大于输出热量。

约束方面，储热罐有三个关键边界条件：储热量上下限SOC_min <= SOC(t) <= SOC_max，充放热功率上下限0 <= Q_ch(t) <= Q_ch_max、0 <= Q_dis(t) <= Q_dis_max，以及调度周期始末储热量相等的要求，保证一个调度周期内储能状态守恒。

这里有个非常容易踩坑的点：充热和放热不能同时进行。物理上当然可以同时充放，但无意义的对充对放会浪费能量并且让结果变得不经济。要建模这个约束，需要引入0-1变量b，写成：

0 <= Q_ch(t) <= b(t) * Q_ch_max 0 <= Q_dis(t) <= (1-b(t)) * Q_dis_max

这才是典型的MILP建模手法。如果没有这对互补约束，模型可能给出同时充放热的"伪最优解"。

2.3 火电机组运行约束

火电部分的约束分两类，一类是传统机组约束，一类是热电联产机组约束。

传统凝汽式机组的约束比较基础：出力上下限约束P_min <= P(t) <= P_max、爬坡约束-DR <= P(t) - P(t-1) <= UR、最小启停时间约束。其中爬坡约束在含风电的调度模型中尤其重要，风电出力的波动性直接考验机组的爬坡能力。

热电联产机组的约束需要格外注意。对抽凝式机组，电出力P和热出力H的可行域是一个多边形的凸包，并不是简单的矩形。典型的可行域约束如下：

P_min + 0.5 * H(t) <= P(t) <= P_max - 0.1 * H(t) H_min <= H(t) <= H_max

第一行是关键——随着热出力增加，电出力可行区间会整体收窄。这就是"以热定电"数学化的核心表达。加装储热罐之后，机组的实际热出力 H(t) 会被拆分为两部分：直供热出力 H_direct(t) 和储热充热 H_storage(t)，即H(t) = H_direct(t) + H_storage(t)。而系统满足的热负荷则是H_direct(t) + H_discharge(t) = H_load(t)。

这个拆分让热负荷不再直接锁定机组的电出力：当风电大发需要机组压电出力时，可以让储热罐承担更多放热任务，机组热出力降低，电出力可行域顺势扩大，风电消纳空间自然打开。

2.4 碳配额与碳交易建模

碳交易建模要注意两个细节，一个是配额的计算方式，一个是碳成本如何进入目标函数。

配额的分配方式我采用的是行业基准线法，即对每台机组设定一个单位电量的碳排放基准值，配额等于基准值乘以发电量：

E_allowance = E_base * P(t)

实际碳排放量的计算基本是线性的：煤耗量乘以排放系数，再考虑机组运行效率。由于煤耗量本身就是出力的函数，实际排放量也可以表达为出力的线性函数。

当实际排放量大于配额时，系统需要在碳市场购买差额；小于配额时，多余的配额可以出售。不管哪种情况，碳成本都是线性的，放进目标函数没有任何问题。我个人在实验中更关注的是碳价变化对调度结构的敏感性——碳价提高时，系统倾向于降低高碳机组出力、增加低碳机组出力和风电消纳、加大储热罐利用以优化机组运行点。这套逻辑在灵敏度分析中可以作为检验模型正确性的一个重要依据。

3. Matlab实现过程与代码架构解析

3.1 程序整体框架

拿到项目第一步不是急着敲代码，而是先把程序架构想清楚。市面上的电力调度Matlab代码五花八门，但好代码的共性是一致的：数据结构清晰、模型与求解分离、结果输出规范。

我采用的框架分为五个模块：

1. 数据输入模块：机组参数、负荷曲线、风电出力曲线、热负荷曲线、碳价、配额基准、储热罐参数，全部集中在数据文件里；
2. 参数预处理模块：把原始数据整理成Yalmip建模需要的格式，比如时序数据、机组索引、分段线性化的断点坐标；
3. 模型构建模块：定义决策变量、目标函数、约束条件；
4. 求解与结果提取模块：调用Cplex求解，检查求解状态，提取各变量结果；
5. 结果分析与可视化模块：绘制出力曲线、SOC曲线、碳排放柱状图等。

模块之间尽量解耦，数据参数在预处理时统一校验，避免后面模型构建时因为维度不匹配来回返工。

3.2 Yalmip变量定义与约束构建

定义变量是整个建模环节最核心的一步。Yalmip里需要区分连续变量和二进制变量。

连续变量：机组电出力、热出力、储热罐充放热功率、SOC、风电实际出力（或者弃风量）、碳排放量等。

二进制变量：机组启停状态、储热罐充放热状态。

举个例子，定义24时段、N台机组的连续变量矩阵：

Pgen = sdpvar(N, 24, 'full'); % 电出力 Hgen = sdpvar(N, 24, 'full'); % 热出力 Qch = sdpvar(N, 24, 'full'); % 充热功率 Qdis = sdpvar(N, 24, 'full'); % 放热功率 SOC = sdpvar(N, 24, 'full'); % 储热罐热储状态 Pw = sdpvar(1, 24, 'full'); % 风电实际出力 delta = binvar(N, 24, 'full'); % 储热充放互斥状态

注意sdpvar第三个参数写成'full'，明确告诉Yalmip这是普通稠密变量矩阵。如果不写，Yalmip默认建立稀疏结构，后续索引时报错会让人抓狂。

约束的构建遵循"先等式后不等式、先单机后系统"的顺序，方便排错。系统功率平衡约束一定不能遗漏，它连接了电源和负荷：

Constraints = [Constraints, sum(Pgen,1) + Pw == P_load(1:24)];

这里的sum(Pgen,1)是在时间维度上对所有机组出力求和，得到每个时段的总电出力。热功率平衡类似：

Constraints = [Constraints, sum(Hgen,1) - sum(Qch,1) + sum(Qdis,1) == H_load(1:24)];

时序索引要小心，P_load如果是从Excel读进来的行向量，维度是1×24，和sum(Pgen,1)维度一致，但如果你读进来的是24×1的列向量，两个向量维度就错开了，Yalmip会直接报错或者隐式广播，最后结果完全不对。

3.3 核心约束的Yalmip写法

储热罐SOC动态方程在Yalmip里写起来非常直观，循环上每个时段：

for t = 2:24 Constraints = [Constraints, SOC(:,t) == SOC(:,t-1) + eta_ch * Qch(:,t) - Qdis(:,t)/eta_dis]; end Constraints = [Constraints, SOC(:,1) == SOC_init + eta_ch * Qch(:,1) - Qdis(:,1)/eta_dis];

注意SOC初值SOC_init的设定很关键。如果把初值设得太低，系统第一时段会约束过紧；设得太高，末时段可能被迫多放热。更合理的做法是设置周期始末SOC相等约束：

Constraints = [Constraints, SOC(:,24) == SOC(:,1)];

同时为了避免第一时段就不可行，可以适当放宽SOC初值范围，允许它在一定区间内自由选择。

充放热互斥约束用二进制变量实现：

Constraints = [Constraints, Qch <= delta * Qch_max]; Constraints = [Constraints, Qdis <= (1-delta) * Qdis_max];

注意这里delta是N×24的矩阵，每台机组每个时段都有一个独立的0-1变量。如果你想让整个系统共享一个储热罐（比如只对一台参考机组做储热改造），那么delta就是1×24的向量，索引方式对应调整。

火电机组的热电耦合可行域约束也是一组不等式：

Constraints = [Constraints, Pgen >= P_min + 0.5 * Hgen]; Constraints = [Constraints, Pgen <= P_max - 0.1 * Hgen];

如果机组类型不一致，这些系数需要分别配置。我在数据输入时用一个2×n的矩阵存每台机组的热电比上下限，而不是写死在公式里，方便后续换算例改参数。

3.4 求解器配置与结果提取

Yalmip调用Cplex的核心代码非常简洁：

ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'showprogress', 1); ops.cplex.mip.tolerances.mipgap = 0.0001; ops.cplex.mip.tolerances.integrality = 1e-06; optimize(Constraints, Objective, ops);

几个参数值得展开说。

mipgap是MILP的停机判据，默认0.0001已经够用，追求快速得到工程可行解可以放宽到0.001，求解时间会明显下降。integrality是0-1变量整数容差，默认1e-06即可，调大了可能导致结果中出现明显的"伪0-1"数值。

求解完成之后，不要急着画图，先检查optimize的返回值和Yalmip的求解状态：

if optimize(Constraints, Objective, ops) == 0 disp('求解成功'); Pgen_value = value(Pgen); SOC_value = value(SOC); else disp('求解失败'); disp(optimize(Constraints, Objective, ops)); end

初学者最容易犯的错就是不管求解状态直接画图。如果问题是不可行的，value(Pgen)返回的是NaN，图倒是画出来了，全图都是洞，根本没法用。而且这种情况下最怕的是你以为是模型正确，但实际上约束已经被Yalmip悄悄修改或者松弛掉了。

4. 实际调试过程中的常见问题与排查心得

4.1 求解器报"Infeasible problem"怎么办

这是做MILP调度模型最常见的噩梦，几乎每个入门的同学都会碰到。维度、约束、参数有一处不对，模型就可能直接变成不可行。

我自己的排查顺序是从松到紧、从内到外。

第一步，先做"无约束测试"。把目标函数设成常数，去掉储热罐SOC的动态约束，只保留所有变量的上下限约束。如果这样还不可行，一定是变量的边界条件冲突。最常见的原因是SOC初值设置和容量边界有冲突，或者某台机组的出力范围下限高于系统负荷，任何一个时段功率平衡都无法满足。

第二步，加上功率平衡，逐步恢复约束。每次只加一组约束，求解一次，直到问题变成不可行，就在最近一次加入的约束里找冲突。这个方法虽然繁琐，但定位问题非常精准。不要一次把所有约束都灌进去，然后把报错信息发给别人看，那排查起来太痛苦了。

第三步，检查储热约束的时序耦合。储热SOC方程天然带有从初始时段到末尾的累计效应，任何一端的能量不平衡都会在末端集中爆发。比如初始SOC设了80%，24小时内储热罐用掉的热量超过总容量，最后时段SOC会突破下限。这种情况下要检查初值选择，或者确认始末SOC约束是否合理。

第四步，检查配额的归属和符号。碳配额如果设置不当，会导致目标函数为负但约束不可行。虽然不常见，但我在测试不同碳价时确实遇到过因为截图问题导致约束范围被排挤的情况。

4.2 求解时间太长怎么优化

MILP求解时间跟二进制变量的数量以及约束的稠密度高度正相关。24时段、每时段多台机组、每台机组又有启停和储热状态两个0-1变量，这种配置的模型规模非常可观。

我实测的经验是，先跑一个小规模算例验证模型逻辑正确，再从6节点系统逐步扩大到几十台机组。不要一上来就跑全规模算例，等模型逻辑完全通了再放开规模。

如果大规模算例确实求解太慢，有几个实用技巧：

• 问题本质上是凸MILP，可以通过设置MIP gap阈值提前终止，比如设置mipgap=0.005，得到的次优解和最优解的偏差在0.5%以内，对调度研究完全够用；
• 把一些布尔状态转化为连续区间约束，比如储热充放互斥，如果不影响方案质量，可以尝试去掉0-1变量，只靠目标函数的经济性引导；
• 热负荷曲线和电负荷曲线的相关性如果比较强，可以适当合并时段，从24时段合并为12时段进行预测试，灵敏度分析阶段再细分。

4.3 常见问题速查表

4.4 结果合理性校验的三个技巧

模型能跑通只是第一步，结果合不合理才是真正见功底的地方。分享三个我在实践中反复使用的校验技巧。

技巧一：观察储热罐SOC曲线的形态。正常的SOC曲线应该是平滑变化的，类似于电池的充放电曲线，不会出现锯齿状剧烈波动。如果SOC曲线在相邻时段剧烈波动，说明充放热功率上下限设置不当或者互斥约束未生效。这个检查能在30秒内发现大量代码Bug。

技巧二：对比改造前后的风电消纳量。同一个风电出力场景下，有储热改造的系统应该比无储热的系统消纳更多风电。如果两个方案的弃风量完全一样，大概率是储热约束没生效，或者热负荷太低导致储热罐没有用武之地。我常用的检查方法是分别跑"有储热"和"无储热"两版模型，对比弃风电量和碳排总量。

技巧三：检查碳排放配额和实际排放的相对关系。合理情况下，系统会倾向于让部分机组处于配额基准附近、部分机组超排购买配额、部分机组减排出售配额。如果所有机组都远低于配额，说明碳价太低，减排没有激励；如果都远高于配额，说明碳价设置过高导致系统过度反应。

5. 参数调试与场景扩展的做法

5.1 碳价的灵敏度测试

碳价是整个低碳调度模型中最值得做的灵敏度参数。我通常的做法是设定一组梯度，比如0、30、50、80、120元/吨，分别求解并记录风电消纳率、系统总煤耗量、碳排放总量、储热罐利用率这几个关键指标。

从实测趋势来看，碳价从0升到30元/吨时，系统碳排放总量下降最明显，因为相对低成本的减排手段（储热利用、负荷转移、机组间出力重新分配）在这个区间内被激活。碳价继续升高到80元/吨以上后，减排效果增速放缓，说明剩余减排空间需要靠更大规模的改造（比如机组灵活性改造、新增储能）才能实现。这种"S型曲线"是符合预期的，如果测试结果没有这个趋势，就要回头检查模型是否真的把碳成本传达到了调度决策。

5.2 两个可以扩展的方向

如果你的课题想做深一点，以下几个方向在现有框架上扩展非常顺滑。

一是多场景鲁棒性分析。当前模型用的是确定性风电出力和负荷曲线。实际中风电出力预测必然有误差，可以引入场景生成和缩减技术，把多个风光出力场景加权后放入目标函数，或者采用鲁棒优化的思路，让调度方案在"最坏风电场景"下依然可行。Yalmip对这类扩展有很好的支持，只需要把变量维度从单场景扩展到多场景再求和平均即可。

二是系统级储热容量规划。当前模型中储热罐参数是给定的，能不能反过来？把储热罐容量也作为决策变量，目标函数里加上储热改造的投资成本，模型就变成了"储热容量优化配置+调度运行"两层决策问题。虽然问题规模变大，但本质上还是一个MILP，Cplex照样能解。这是很多期刊论文惯用的扩展思路，适合在课题中后期往这个方向深入。

6. 一些补充经验和最终建议

最后聊几个我在这类项目上积累的实操细节，都是踩过坑之后总结出来的。

第一个是关于热负荷曲线数据的处理。很多公开数据集里热负荷的数据精度不高，甚至有缺失值。做调度之前一定要先把热负荷曲线可视化一遍，看看是否有明显的不合理跳变。如果热负荷出现尖峰但储热罐放热上限不够，模型的调度结果无论怎么调都会在对应时段出现功率不平衡，最后表现为不可行。处理办法是适当平滑数据或者调大储热罐功率上限，而不是强行更改约束。

第二个是关于Matlab内存和计算环境的建议。大规模MILP求解时，Cplex占用内存非常大。我遇到过求解过程中Matlab直接卡死的情况，后来排查发现是循环里不停地往约束集里添加约束，导致Yalmip内部数据结构反复重建。解决办法是在建模前预先分配Constraints变量，或者每加完一段约束就调用optimize验证，而不是把所有约束堆到最后一次性求解。

第三个是关于结果可视化的细节。对应的调度结果除了画常规的电源出力堆叠图之外，强烈建议再画一张"热力系统平衡图"，把每个时段的直供热量、储热充热量、储热放热量、热负荷量四条曲线画在一起，一眼就能看到储热罐在哪些时段充热、哪些时段放热、什么时候在"搬移热量"。这张图对于论文写作和答辩展示的说服力非常强，效果远超一张密密麻麻的出力表格。

如果你正在做这个方向，我的建议是不要急着追求大算例、复杂模型，先把一个简单的场景跑通、把结果校验清楚，理解清楚储热改造到底给调度模型带来了什么样的自由度变化，再逐步往里加细节。跟储热相关的课题，不管最后是写毕业论文还是发期刊，核心卖点都在于"热-电-碳"三者之间的耦合关系。把这个关系表达清楚，模型的框架就立住了，剩下的都是工程实现层面的问题。