考虑火电机组储热改造的电力系统低碳经济调度这几年做电力系统优化调度方向的课题有一个绕不开的痛点新能源占比越来越高但火电依然承担着调峰和供热的主力角色。尤其是北方地区热电联产机组占比极大以热定电的刚性约束让机组在冬季深调能力非常有限弃风弃光问题在供暖期集中爆发。储热改造就是针对这一矛盾的核心手段之一——在火电机组侧加装储热罐把热负荷和电出力解耦给机组腾出更多调节空间。而低碳经济调度则是在传统经济调度框架之上进一步叠加碳配额、碳交易成本、碳流追踪等因素把碳排放变成调度决策中的可量化成本从而引导系统走向更低碳的运行方式。这篇文章总结的是我用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代码五花八门但好代码的共性是一致的数据结构清晰、模型与求解分离、结果输出规范。我采用的框架分为五个模块数据输入模块机组参数、负荷曲线、风电出力曲线、热负荷曲线、碳价、配额基准、储热罐参数全部集中在数据文件里参数预处理模块把原始数据整理成Yalmip建模需要的格式比如时序数据、机组索引、分段线性化的断点坐标模型构建模块定义决策变量、目标函数、约束条件求解与结果提取模块调用Cplex求解检查求解状态提取各变量结果结果分析与可视化模块绘制出力曲线、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阈值提前终止比如设置mipgap0.005得到的次优解和最优解的偏差在0.5%以内对调度研究完全够用把一些布尔状态转化为连续区间约束比如储热充放互斥如果不影响方案质量可以尝试去掉0-1变量只靠目标函数的经济性引导热负荷曲线和电负荷曲线的相关性如果比较强可以适当合并时段从24时段合并为12时段进行预测试灵敏度分析阶段再细分。4.3 常见问题速查表现象可能原因排查方法求解状态InfeasibleSOC初值和末值冲突去掉始末SOC约束测试模型是否可解求解状态Infeasible功率平衡约束和负荷曲线不匹配检查P_load维度、单位是否转化正确结果中充放热同时非零未添加0-1互斥约束引入binvar并添加互补约束风电出力曲线几乎为0弃风惩罚系数设置过低增大弃风惩罚系数到煤耗成本的10倍以上储热罐SOC长期保持上限或下限储热罐容量过小或热负荷缺乏波动性调整储热罐容量检查热负荷曲线碳交易成本对结果无影响碳价设置过低或配额过宽松提高碳价做灵敏度分析结果严重不连续机组启停0-1变量数量过少增加时段数或给启停变量添加过渡约束求解时间过长0-1变量过多、约束矩阵过密设置mipgap、削减时段数、变量降维4.4 结果合理性校验的三个技巧模型能跑通只是第一步结果合不合理才是真正见功底的地方。分享三个我在实践中反复使用的校验技巧。技巧一观察储热罐SOC曲线的形态。正常的SOC曲线应该是平滑变化的类似于电池的充放电曲线不会出现锯齿状剧烈波动。如果SOC曲线在相邻时段剧烈波动说明充放热功率上下限设置不当或者互斥约束未生效。这个检查能在30秒内发现大量代码Bug。技巧二对比改造前后的风电消纳量。同一个风电出力场景下有储热改造的系统应该比无储热的系统消纳更多风电。如果两个方案的弃风量完全一样大概率是储热约束没生效或者热负荷太低导致储热罐没有用武之地。我常用的检查方法是分别跑有储热和无储热两版模型对比弃风电量和碳排总量。技巧三检查碳排放配额和实际排放的相对关系。合理情况下系统会倾向于让部分机组处于配额基准附近、部分机组超排购买配额、部分机组减排出售配额。如果所有机组都远低于配额说明碳价太低减排没有激励如果都远高于配额说明碳价设置过高导致系统过度反应。5. 参数调试与场景扩展的做法5.1 碳价的灵敏度测试碳价是整个低碳调度模型中最值得做的灵敏度参数。我通常的做法是设定一组梯度比如0、30、50、80、120元/吨分别求解并记录风电消纳率、系统总煤耗量、碳排放总量、储热罐利用率这几个关键指标。从实测趋势来看碳价从0升到30元/吨时系统碳排放总量下降最明显因为相对低成本的减排手段储热利用、负荷转移、机组间出力重新分配在这个区间内被激活。碳价继续升高到80元/吨以上后减排效果增速放缓说明剩余减排空间需要靠更大规模的改造比如机组灵活性改造、新增储能才能实现。这种S型曲线是符合预期的如果测试结果没有这个趋势就要回头检查模型是否真的把碳成本传达到了调度决策。5.2 两个可以扩展的方向如果你的课题想做深一点以下几个方向在现有框架上扩展非常顺滑。一是多场景鲁棒性分析。当前模型用的是确定性风电出力和负荷曲线。实际中风电出力预测必然有误差可以引入场景生成和缩减技术把多个风光出力场景加权后放入目标函数或者采用鲁棒优化的思路让调度方案在最坏风电场景下依然可行。Yalmip对这类扩展有很好的支持只需要把变量维度从单场景扩展到多场景再求和平均即可。二是系统级储热容量规划。当前模型中储热罐参数是给定的能不能反过来把储热罐容量也作为决策变量目标函数里加上储热改造的投资成本模型就变成了储热容量优化配置调度运行两层决策问题。虽然问题规模变大但本质上还是一个MILPCplex照样能解。这是很多期刊论文惯用的扩展思路适合在课题中后期往这个方向深入。6. 一些补充经验和最终建议最后聊几个我在这类项目上积累的实操细节都是踩过坑之后总结出来的。第一个是关于热负荷曲线数据的处理。很多公开数据集里热负荷的数据精度不高甚至有缺失值。做调度之前一定要先把热负荷曲线可视化一遍看看是否有明显的不合理跳变。如果热负荷出现尖峰但储热罐放热上限不够模型的调度结果无论怎么调都会在对应时段出现功率不平衡最后表现为不可行。处理办法是适当平滑数据或者调大储热罐功率上限而不是强行更改约束。第二个是关于Matlab内存和计算环境的建议。大规模MILP求解时Cplex占用内存非常大。我遇到过求解过程中Matlab直接卡死的情况后来排查发现是循环里不停地往约束集里添加约束导致Yalmip内部数据结构反复重建。解决办法是在建模前预先分配Constraints变量或者每加完一段约束就调用optimize验证而不是把所有约束堆到最后一次性求解。第三个是关于结果可视化的细节。对应的调度结果除了画常规的电源出力堆叠图之外强烈建议再画一张热力系统平衡图把每个时段的直供热量、储热充热量、储热放热量、热负荷量四条曲线画在一起一眼就能看到储热罐在哪些时段充热、哪些时段放热、什么时候在搬移热量。这张图对于论文写作和答辩展示的说服力非常强效果远超一张密密麻麻的出力表格。如果你正在做这个方向我的建议是不要急着追求大算例、复杂模型先把一个简单的场景跑通、把结果校验清楚理解清楚储热改造到底给调度模型带来了什么样的自由度变化再逐步往里加细节。跟储热相关的课题不管最后是写毕业论文还是发期刊核心卖点都在于热-电-碳三者之间的耦合关系。把这个关系表达清楚模型的框架就立住了剩下的都是工程实现层面的问题。