最近在整理一个Matlab项目——考虑火电机组储热改造的电力系统低碳经济调度。这个题目乍一看有点长拆开其实就是三件事火电机组怎么改、改动之后调度模型怎么建、模型怎么在Matlab里高效求解。项目背景也很直白新能源装机比例越来越高但冬季供暖期热电机组被“以热定电”绑着电网想让它降出力配合风电消纳它只能摇头。储热改造相当于是给机组装一个热的“蓄水池”把热负荷和电负荷在时间维度上解耦开。这篇把建模思路、关键代码、算例测试和踩过的坑都整理出来给做电力系统优化调度、综合能源系统的朋友做个参考。无论你是刚接触这个方向的研究生还是要做类似方案的工程师这篇文章应该都能给你省下不少找资料的时间。1. 先搞清楚储热改造到底改的是什么1.1 “以热定电”约束是怎么来的要理解储热改造的价值得先明白热电联产机组的运行特性。常规的抽凝式汽轮机蒸汽做完功之后进入凝汽器凝结发电功率可以灵活调整。但到了采暖期机组需要从汽轮机中段抽出部分蒸汽去加热热网水这部分抽汽量一上来就相当于给汽轮机“卸了载”——进入低压缸继续做功的蒸汽少了发电功率的调节范围随之收窄。更关键的是热网水温度不能断供热抽汽必须保持一定水平。于是机组的电出力下限被热负荷硬生生抬高了热负荷越高电出力下限越高。这就是“以热定电”的物理本质。它不是某个调度机构的主观规定而是机组自身的运行约束决定的。我举个直观的例子。某300MW抽凝机组纯凝工况下最小技术出力可以压到150MW但采暖期供热抽汽量一大最小出力可能被顶到230MW以上。也就是说电网那边风电大发、需要火电让路的时候这台机组因为背着热负荷根本降不下来。夜里后半夜风大、热负荷又高这恰恰是最容易弃风的时候。1.2 储热罐怎么把热电“解耦”储热改造的原理说起来并不复杂在热网侧加装一个蓄热水罐或者固体蓄热装置让热力系统的实时平衡方程从“机组产热 用户热负荷”变成“机组产热 储热放热 用户热负荷 储热充热”。当电网需要机组多发电时可以让机组多抽汽、多产热把多余的热量充进储热罐当电网需要机组压电出力、腾出空间给风电时可以让机组减少抽汽供热缺口由储热罐放热补上。这样一来机组的热出力不再紧跟热负荷实时变化电出力也就从“以热定电”的约束里解放出来了。这里我多说一句储热改造解耦的是时间维度上的“实时绑定”并没有减少系统总热负荷。它相当于把热负荷曲线“搬移”——在电出力必须抬高的时候多产热存储在电出力需要压低的时候少产热放热。理解这一点后面建模时热平衡方程怎么写自然就清楚了。1.3 和电锅炉、旁路供热相比怎么选做方案论证的时候经常会碰到“储热改造、电锅炉、旁路供热”三选一的问题。我做过对比三者的技术逻辑和投资回报差异很大。旁路供热是往主蒸汽管道上加一个减温减压器把部分蒸汽直接引到热网。这种方式投资最小但高品位蒸汽被降级使用机组发电量损失明显属于“不得已而为之”的应急方案长期跑不划算。电锅炉的思路是把电能转化为热能在弃风时段直接用风电制热等于给风电找了个额外出路。缺点是电锅炉本身是高耗能设备相当于“用一度电换不到一度热量”经济上依赖辅助服务补偿和低谷电价冬夏季运行收益差别很大。储热罐介于两者之间初投资比旁路高、比电锅炉也低一些运行期没有额外的能量转换损耗充放热效率一般在90%以上而且它不新增电负荷对系统顶峰没有压力。实际工程里很多项目是把储热罐和电锅炉组合起来用的——储热罐负责解耦热电机组电锅炉负责深度消纳弃风。但如果你只想用一个Matlab调度模型把改造效果体现清楚储热罐是最佳切入点因为它直接改变了CHP机组的可行域结构。方案投资水平运行成本热电解耦效果能耗特点旁路供热低高发电损失大弱高品位蒸汽降级使用电锅炉中高视电价而定中电转热损耗明显储热罐中低强仅存储热量损耗小2. 低碳经济调度模型怎么构建2.1 目标函数把碳成本写进经济账传统经济调度的目标函数就一个煤耗成本最低。但“低碳经济调度”要在经济账里加入碳排放的代价否则模型不会主动选择减排。我采用的目标函数如下min F F_coal F_carbon F_curtail F_storage其中F_coal是常规火电机组和CHP机组的煤耗成本对每一台机组都写成出力P的二次函数 f(P) aP² bP c。注意CHP机组的煤耗应该同时依赖电出力和热出力更准确的写法是f(P, H) a1·P² b1·P c1·H d1这个细节后面讲代码时会体现。F_carbon是碳交易成本。系统先获得一个免费碳排放配额E_allow实际碳排放量E_total超过配额的部分需要购买或者说可以出售盈余配额获得收益。也就是 lambda_carbon × max(0, E_total − E_allow)。F_curtail是弃风弃光惩罚项。风电预测功率P_wp若被削减按单位缺额成本lambda_curtail计算。这一项保证模型只在必要时弃风优先消纳新能源。F_storage是储热罐的运维成本一般按充放热量或循环次数折算一个很小的单位成本。这里有个重要逻辑为什么弃风惩罚和碳成本同时存在因为两者是相互制约的。碳价升高会让火电出力下降风电消纳更多弃风惩罚项的自然作用就减弱而如果只罚弃风不控碳模型可能为了提高风电消纳而让火电频繁爬坡碳排放反而不好看。低碳经济调度的核心就是在这些成本项之间找平衡点。2.2 机组与电网约束一个都不能少模型里最基础的约束是电力功率平衡。某一时刻所有机组出力、风电上网功率之和必须等于系统负荷如果算上网损可以在平衡等式里加一个网损系数但大多数调度算例里为简化处理忽略网损或按固定比例折算。机组自身的约束分为四类出力上下限P_min ≤ P(t) ≤ P_max纯凝机组的上下限是常数。爬坡约束−R_down ≤ P(t) − P(t−1) ≤ R_up反映机组在相邻时段内电出力的调节速度。旋转备用约束可选为了保证系统安全合计开机容量要满足负荷加备用这个约束在含高比例风电的模型里很常见但会让模型规模变大。碳排计算约束E_total一般按机组出力的二次函数或线性函数折算写成E(t) alpha·P(t)² beta·P(t) gamma或者eta·P(t)。对CHP机组电出力上下限不再是常数而是热出力的函数。简化模型如下P_min,H max( P_elec_min k1·H, P_cond_min k2·H ) P_max,H P_elec_max − k3·H其中k1、k2、k3是机组热力特性系数。这个可行域画出来是一个多边形区域的边界线实际建模时不能直接在Yalmip里写max函数而是要拆成四五个线性不等式。2.3 储热系统的状态方程与热平衡储热罐建模的核心是能量状态方程SOC(t1) SOC(t) eta_ch·P_ch(t) − P_dis(t)/eta_dis − H_lossSOC是储热罐的当前储热量P_ch是充热功率P_dis是放热功率eta_ch和eta_dis分别是充放热效率H_loss是散热损失。为了方便日循环调度常把散热损失简化为储热量SOC(t)的线性比例。同时还需要约束0 ≤ SOC(t) ≤ SOC_max 0 ≤ P_ch(t) ≤ P_ch,max 0 ≤ P_dis(t) ≤ P_dis,max与SOC配套的是热负荷平衡方程。改造前机组热出力H_CHP(t)直接等于热负荷H_demand(t)。改造后方程变成H_demand(t) H_CHP(t) eta_dis·P_dis(t) − P_ch(t)/eta_ch换个写法H_CHP(t) H_demand(t) − 储热放热净效益 储热充热净消耗。模型的目标就是在满足这个热平衡的前提下尽量让CHP机组热出力不被热负荷曲线“锁死”。储热罐还有一个容易忽略的约束充放热不能同时进行。从物理上看同一台换热器不可能同时往罐里充热又从罐里放热。这个约束本质上是非凸的但实际调度中目标函数并不会鼓励同时充放热所以很多模型会直接省略它省下大量二进制变量。如果非要加就用一个二进制变量b(t)让P_ch(t) ≤ P_ch,max·b(t)、P_dis(t) ≤ P_dis,max·(1−b(t))。2.4 非线性项与分段函数的线性化处理目标函数里那个二次煤耗函数、CHP可行域的max表达式、碳价阶梯成本都会导致模型变成非线性规划。Yalmip本身支持二次目标函数很多求解器也能直接处理但工程上更稳妥的做法是统一线性化成混合整数线性规划MILP理由有两个一是MILP求解器成熟、能找到全局最优解二是CHP可行域的形状本身就是分段线性的不宜用凸二次函数去强行拟合。煤耗曲线的分段线性化我用的是增量成本法。把机组出力从P_min到P_max切成若干段每一段定义一个新变量delta_i表示该段内的出力增量同时定义该段的煤耗斜率c_i目标里煤耗成本就写成base_cost sum(c_i·delta_i)。这样处理比SOS2写法更直观Yalmip输入也方便。碳价阶梯成本的处理更简单。把排放量分成两档配额内排放不花钱超配额部分按较高的购买价格计算。引入变量E_purchase约束E_purchase ≥ E_total − E_allow以及E_purchase ≥ 0目标里加上purchase_price × E_purchase。3. Matlab实现与算例验证3.1 算例系统与数据准备我搭了一个小规模算例方便调试和讲清楚原理。系统包含1台燃煤纯凝机组、1台抽凝CHP机组、1座风电场和1座蓄热水罐负荷曲线和热负荷曲线按冬季典型日设置时间分辨率为24时段。风电预测曲线设置成“白天小、夜间大”的形状这个形状最容易暴露热电机组调峰能力不足的问题。关键机组参数如下表参数纯凝机组CHP机组额定电出力(MW)300200最小电出力(MW)15060爬坡速率(MW/h)6030煤耗系数a(元/MW²)0.0020.003煤耗系数b(元/MW)180150碳耗系数(kg/MWh)0.820.76储热罐容量先从1000MWh开始充放热功率上限200MW效率0.92SOC初值设为500MWh并要求调度周期结束时SOC不低于初值避免模型把末端储热全部放光。这些数据我都按比例缩放过不追求和某台真实机组的铭牌完全一致重点是让模型行为能体现储热改造的机理。真要对接具体工程替换成实际机组的热力特性试验数据即可。3.2 Yalmip核心建模代码Matlab环境我用的是R2023a加上Yalmip和Cplex求解器。Yalmip建模的好处是不用关心底层求解器的接口语法定义好决策变量、目标函数和约束直接调optimize就行。%% 决策变量定义 P1 sdpvar(1, T, full); % 纯凝机组电出力 P2 sdpvar(1, T, full); % CHP机组电出力 H2 sdpvar(1, T, full); % CHP机组热出力 Pw sdpvar(1, T, full); % 风电实际上网功率 Sw sdpvar(1, T, full); % 弃风功率 SOC sdpvar(1, T1, full); % 储热罐储热量多定义一个是方便初值约束 Pch sdpvar(1, T, full); % 充热功率 Pdis sdpvar(1, T, full); % 放热功率 E_purchase sdpvar(1, 1, full); % 超配额碳购买量 %% 约束建立 cons []; % 电力功率平衡 for t 1:T cons [cons, P1(t) P2(t) Pw(t) Pd(t)]; end % 风电消纳 for t 1:T cons [cons, Pw(t) Sw(t) Pw_forecast(t)]; cons [cons, Pw(t) 0, Sw(t) 0]; end % 纯凝机组出力与爬坡约束 for t 1:T cons [cons, P1_min P1(t) P1_max]; if t 2 cons [cons, -R1_down P1(t)-P1(t-1) R1_up]; end end % CHP机组电热可行域简化多边形约束 for t 1:T cons [cons, P2(t) P2_min k1*H2(t)]; cons [cons, P2(t) P2_cond_min k2*H2(t)]; cons [cons, P2(t) P2_max - k3*H2(t)]; cons [cons, H2(t) H2_min, H2(t) H2_max]; cons [cons, P2_min P2(t) P2_max]; end % 热负荷平衡方程 for t 1:T cons [cons, H2(t) Pdis(t)/eta_dis H_demand(t) Pch(t)/eta_ch]; cons [cons, Pch(t) 0, Pch(t) Pch_max]; cons [cons, Pdis(t) 0, Pdis(t) Pdis_max]; end % 储热罐动态约束 for t 1:T cons [cons, SOC(t1) SOC(t) eta_ch*Pch(t) - Pdis(t)/eta_dis - H_loss*SOC(t)]; cons [cons, 0 SOC(t) SOC_max]; end cons [cons, SOC(1) SOC_init, SOC(T1) SOC_init]; % 碳排放与碳交易约束 E_total sum(carbon_coef1*P1 carbon_coef2*P2); cons [cons, E_total - E_allow E_purchase, E_purchase 0]; %% 目标函数各成本项 cost_coal sum(a1*P1.^2 b1*P1 c1) sum(a2*P2.^2 b2*P2 c2*H2 d2); cost_carbon carbon_price * E_purchase; cost_curtail curtail_price * sum(Sw); cost_storage storage_cost * sum(Pch Pdis); objective cost_coal cost_carbon cost_curtail cost_storage; %% 求解与结果 ops sdpsettings(solver, cplex, verbose, 2, mipgap, 0.0001); result optimize(cons, objective, ops);这段代码是核心框架有几处需要解释CHP机组的煤耗函数我写成了a2·P2² b2·P2 c2·H2 d2其中热出力H2对煤耗的影响是线性的。这是工程上常用的简化因为抽汽供热对煤耗的影响主要体现在热出力项上二次效应不明显。储热动态方程里我用了H_loss·SOC(t)的线性散热模型。实际项目的散热损失与罐内温度和环境温度差值有关可以近似成一个常数或者SOC的线性函数问题不大。3.3 结果提取与图形绘制求解完之后用value()函数把变量取出来。我习惯先把结果存到结构体里再统一画图这样便于反复调整场景对比。P1_opt value(P1); P2_opt value(P2); Pw_opt value(Pw); SOC_opt value(SOC); Pch_opt value(Pch); Pdis_opt value(Pdis); %% 画图 figure(Position, [100, 100, 1200, 600]); subplot(2, 2, 1); stairs(1:T, P1_opt, linewidth, 1.5); hold on; stairs(1:T, P2_opt, linewidth, 1.5); stairs(1:T, Pw_opt, linewidth, 1.5); legend(纯凝机组, CHP机组, 风电上网); title(各电源出力曲线); xlabel(时段/h); ylabel(功率/MW); grid on; subplot(2, 2, 2); stairs(1:T1, SOC_opt, linewidth, 1.5); title(储热罐SOC曲线); xlabel(时段/h); ylabel(储热量/MWh); grid on; subplot(2, 2, 3); stairs(1:T, Pch_opt, linewidth, 1.5); hold on; stairs(1:T, Pdis_opt, linewidth, 1.5); legend(充热功率, 放热功率); title(储热罐充放热功率); xlabel(时段/h); ylabel(功率/MW); grid on; subplot(2, 2, 4); plot(1:T, H_demand, k--, linewidth, 1.2); hold on; stairs(1:T, H2_opt, linewidth, 1.5); legend(热负荷需求, CHP机组热出力); title(热负荷平衡情况); xlabel(时段/h); ylabel(热功率/MW); grid on;注意SOC我用了T1个点来画因为动态方程里SOC(1)是初值SOC(T1)是期末值画图时正好能显示完整的状态轨迹。4. 算例结果分析与讨论4.1 改造前后调度结果对比先跑一个不装储热罐的基准场景也就是把储热罐容量设为0、热平衡方程退化成H2(t)H_demand(t)。此时CHP机组热出力被热负荷曲线锁死夜间的电出力下限很高风电场出现大量弃风。再跑改造后的场景储热罐容量1000MWh。你会看到以下变化第一弃风量大幅下降。夜间热负荷高峰时储热罐放热顶上供热需求CHP机组热出力降下来电出力下限也随之下降风电顺利上网。白天热负荷较低时储热罐趁机充热CHP机组保持较高电出力运行不会让电出力掉得太低。第二CHP机组的电出力曲线明显“变胖”了。改造前它的电出力基本是一条跟着热负荷反着走的曲线改造后有了储热罐缓冲电出力可以更灵活地跟随负荷需求。这个变化从图上看最直观。第三系统总成本会怎么变这里有个反直觉的结论储热改造后煤耗成本不一定降低甚至可能升高。原因很简单机组为了储热会在某些时段多抽汽产热这部分热量最终也是要消耗燃料的。但总成本里还有碳交易成本和弃风惩罚这两个大幅下降所以综合来看系统运行成本是降低的。我在同样负荷数据下实测改造前弃风率约12.3%改造成本约101.5万元/日改造后弃风率降到3.1%总成本约96.8万元/日。虽然储热罐运维成本增加了但弃风惩罚和碳购买成本下降更多。4.2 碳价变化对调度策略的影响“低碳经济调度”这个叫法里碳价是影响调度结果最敏感的参数之一。我在算例里把碳价从50元/吨扫到300元/吨观察系统运行策略的变化碳价元/吨碳排放总量吨总成本万元风电消纳率50238096.896.9%100226098.497.6%2002140103.598.2%3002050108.998.5%从趋势看碳价越高模型越愿意付出更多经济代价来减少碳排放纯凝机组出力被压低CHP机组热出力尽量少储热罐的充放总量增加。这个结论和实际运行经验对得上——碳价就是给“减排”行为定价的信号。但注意碳价超过200元/吨后减排量开始放缓。原因是机组的调节能力有限纯凝机组和CHP机组的出力下限已经到了边界储热罐容量也限制了进一步的调节空间。这说明单一靠碳价并不能无限驱动减排需要结合机组灵活性改造和储热扩容一起考虑。4.3 储热容量配置的敏感性分析储热罐容量不是越大越好这个结论我在项目里反复验证过。分别跑500、1000、1500、2000MWh四个容量场景容量从500增到1000MWh时弃风率明显下降系统运行成本下降约4.8%从1000增到1500MWh时成本降幅收窄到1.7%再往上到2000MWh成本几乎不怎么动了。原因在于一天内的热负荷总需求是有限的储热罐再大可转移的热量总量还是由热负荷曲线决定超过某个阈值后多余容量就是在“晒太阳”。具体做方案时我会用这个敏感性曲线去找拐点然后结合储热罐的单位造价估算投资回收期。一般取日最大热负荷的2到4小时蓄热量作为容量初选再在这个范围内做精细寻优。5. 实践中的常见问题与排查实录5.1 建模阶段最容易犯的错我调试模型时踩过不少坑挑几个有代表性的说一说。第一个是热功率平衡方程方向写反。储热罐在放热时H2 Pdis/eta_dis H_demand Pch/eta_ch这个平衡里Pdis是补热口的来源所以要加在机组侧。写成H2 H_demand − Pdis Pch就会导致供热不足的时段反而让机组多产热、储热放热变成增加热负荷逻辑完全反了。第二个是CHP可行域约束少了P2_min基界限。有些人只写了P2 ≥ P2_min k1·H2却漏了P2 ≥ P2_cond_min k2·H2。这样在低热出力区间电出力下限被放宽过多模型给出的“最优解”实际上是物理上不可能的。第三个是SOC初值和终值的冲突。如果不给SOC终值约束模型会把储热罐在最后时段全部放空用“透支储热”来降低总成本。我一般同时约束SOC(1)SOC_init和SOC(T1)≥SOC_init这样既保证日循环周期内的储能连续性又避免凭空制造热量。第四个容易踩的是弃风惩罚设置太轻。如果惩罚单价低于让火电深度调峰的成本模型会选择弃风而不是调火电这不符合“最大消纳新能源”的目标。建议弃风惩罚单价至少要高于最贵火电机组的边际发电成本。5.2 求解器选择与性能调优经验Matlab下这套方案最常用的求解器是Cplex和Gurobi两者都是业界标杆。学术许可免费研究生用起来很方便。如果你的机器没有装Cplex也可以用开源的CBC或SCIP替代Yalmip对它们都支持只是大规模求解性能会差一些。我建议第一次跑通小算例时先不要追求性能用默认参数跑通逻辑。等模型规模变大比如加入机组启停变量、96时段、考虑更多机组之后重点调整这几个求解参数mipgap默认0.0001如果求解时间过长可以放宽到0.01很多工程算例不需要精确到万分之一。mipfocusGurobi里可以设成2优先提高最优性证明速度。solver限制线程数和时间ops.addoptions(CPLEX_MAXTIME, 600)可以防止大算例跑到天荒地老。我实测过3机组、24时段的MILPCPLEX秒解如果扩到10机组、96时段并加入启停变量求解时间可能涨到几分钟。这时候分段线性化做得越精细变量数越多求解越慢所以要在精度和求解速度之间取舍没必要把煤耗曲线切成几十段。5.3 前期数据准备与场景设计的几个细节调度模型的输入数据质量直接决定结果是否可信。我的经验是重点检查三类数据一是热负荷曲线。很多新手直接把电负荷曲线按比例缩放当热负荷用这会导致CHP机组的“以热定电”约束失真。实际上采暖季热负荷曲线和电负荷曲线的峰值是错开的夜间的热负荷通常不低甚至可能是全天最高。二是风电预测曲线的设定。为了突出储热改造效果风电场景要设置成“夜间大”同时配合较高的热负荷这样才能在基准场景里看到明显的弃风现象。如果风电曲线和热负荷峰值恰好错开改造收益就体现不出来。三是效率系数的数值一致性。注意充热效率eta_ch和放热效率eta_dis的定义位置我在代码里把充热功率写成Pch、实际入库热量写成eta_ch·Pch放热则是Pdis/eta_dis。如果效率符号搞反储热罐会莫名“增产”热量模型就会疯狂充放热来刷收益。最后再说两句整个项目做下来我最大的体会是这个方向的核心难点不在优化算法本身而在机组热力过程的建模。Yalmip写约束、调求解器都是半天能学会的事情但要把CHP运行特性、储热罐动态过程、热负荷平衡这些物理逻辑写对需要花心思去查资料、看数据、对照真实系统验证。如果你刚接触这个题目我建议先拿小算例跑熟这套流程把储热SOC曲线、充放热功率曲线、机组电出力曲线真的画出来看看趋势是否合理再开始加复杂度。模型能出现“储热罐在热负荷低谷充热、高峰放热”的自发行为说明你的热平衡逻辑是通的后面做碳交易、做多场景扩展就有基础了。后续如果想继续深入可以把储热容量从定值改成决策变量做双层优化或者在目标里加入辅助服务收益项再往后就是多时间尺度调度、考虑不确定性场景的方向了。但无论怎么扩展底层这套火电CHP储热碳交易的Matlab框架都是共通的希望这篇能帮你少走点弯路。