考虑储热改造的热电联产电力系统低碳经济调度Matlab实现

做电力系统优化调度的人，估计对“以热定电”这四个字都不陌生。北方冬季热电联产机组一开，发电出力被供热需求绑得死死的，电负荷低的时候也没法降太多，风电只能在夹缝里求生存。“考虑火电机组储热改造的电力系统低碳经济调度”这个题目，说白了就是给火电机组加一个储热罐，把热和电解耦，然后在满足负荷的前提下让系统的碳排放和总运行成本一起降下来。我最近把这个课题从建模到Matlab代码实现完整跑了一遍，这篇文章把过程拆开讲讲，从模型设计到代码框架再到算例验证都覆盖到，适合正在做调度优化、或者想了解储热改造实际价值的同学参考。

1. 项目在解决什么问题：从“以热定电”到热电解耦

1.1 “以热定电”是怎么拖累风电消纳的

传统热电联产机组在冬季供暖期，运行方式基本是“以热定电”：热负荷有多高，机组电出力就被约束在一个相应范围内。具体来说，抽凝式机组的电出力下限会随着供热抽汽量的增加而抬高，背压式机组干脆是“发多少电就供多少热”，几乎没有调节弹性。这套逻辑在以前电网结构简单、负荷波动不大时没什么问题，但风电大规模接入之后就麻烦了。

风电的特点大家都知道，夜间出力往往比白天高，而夜间恰恰是电负荷低谷、热负荷却很高的时段。此时热电机组为了满足夜间供热需求，必须维持较高的最小电出力，比如一台300MW的机组可能最低只能压到180MW。如果当晚的电网负荷只有400MW，风电预测出力有150MW，热电机组再占走300MW以上的空间，风电就被挤出去了，只能弃掉。弃风不仅浪费清洁能源，还抬高了系统整体的碳排放水平，因为被挤掉的风电得靠煤电补上。

储热改造的核心价值就在这里：在热电机组和热网之间嵌入一个储热罐，相当于给“供热”加了缓冲池。热负荷高的时候可以由储热罐放热补充，机组不必为供热强撑高电出力；热负荷低或电负荷高的时候，机组多发电并给储热罐充热。这样一来，“热”和“电”不再是硬绑定的关系，机组发电出力有了更大的调节区间，风电消纳空间自然就出来了。

1.2 低碳经济调度到底是一个什么样的优化问题

低碳经济调度本质上是一个多目标优化问题，但工程实现上通常把它揉成一个单目标函数来处理。低碳和经济并不总是同方向的：压低煤电出力可以降碳排放，但可能引起备用不足或需要频繁启停机组，经济上未必划算；反过来，让机组在高效区满发虽然煤耗低，碳排放总量可能更高。所以目标函数里需要既有煤耗成本、又有碳排放成本，还要有弃风惩罚成本，通过权重系数（比如碳价、惩罚价格）把它们统一到“总成本”这一个量纲里。

约束条件方面，要比常规经济调度多出不少东西：热功率平衡约束、热电耦合约束、储热罐容量和充放热状态约束，以及储热罐的“初始储热量 = 终态储热量”循环约束。决策变量也多了储热罐的充热功率、放热功率和相应的0-1状态变量，问题从线性规划变成了混合整数规划。这也直接决定了Matlab代码实现的复杂度。

2. 模型构建：储热装置物理建模与目标函数设计

2.1 储热罐和热电机组的物理约束怎么写

储热罐的模型看起来很复杂，其实能量守恒一句话就能说清楚：相邻两个时段之间，罐内储存的热量变化量，等于流入热量减去流出热量，再减去散热损失。写成数学表达式就是：

S(t) = S(t-1) + η_ch * Q_ch(t) - Q_dis(t) / η_dis - μ * S(t-1)

其中S是储热罐当前储存的热量，Q_ch是充热功率，Q_dis是放热功率，η_ch和η_dis分别是充放热效率，μ是散热损失率。充热和放热不能同时进行，这是储热罐物理特性决定的，后续要靠两个0-1变量来约束。

热电机组部分，完整的抽凝式机组可行域是一个四边形围成的区域，手工建模比较繁琐。实际课题中更常用的是简化模型：认为电出力上下限都和供热出力线性相关，即

P_i,min + α_i * H_i(t) <= P_i(t) <= P_i,max - β_i * H_i(t)

其中H_i是热电机组供热量，α_i表示供热对电出力下限的提升系数，β_i表示供热对电出力上限的压低系数。现实中α_i通常为正，意味着供热量越大，机组为保供热的最低电出力越高，这正是“以热定电”的数学来源。储热改造后，供热平衡变成：

Σ H_i(t) + Q_dis(t) - Q_ch(t) = H_load(t)

也就是说，热负荷可以由热电机组直接供热，也可以由储热罐放热来满足，机组的一部分供热义务被“转移”到了储热罐身上。

2.2 目标函数：煤耗、碳排放、弃风惩罚怎么加权

调度模型的目标函数我一般写成四项加和：

min Σ_t [ Σ_i (a_i * P_i(t)^2 + b_i * P_i(t) + c_i) + λ_carbon * Σ_i ε_i * P_i(t) + λ_curtail * P_curtail(t) + λ_storage * (Q_ch(t) + Q_dis(t)) ]

第一项是煤耗成本，二次项描述机组煤耗随出力增加而上升的非线性关系；第二项是碳排放成本，用机组碳排放强度ε_i乘以出力再乘碳价λ_carbon得到；第三项是弃风惩罚，P_curtail是弃风功率，λ_curtail设得够高才能体现“优先消纳风电”的导向；第四项是储热罐的充放热运行维护成本，系数通常很小，主要是为了防止模型在无关紧要的地方随意动作。

这里有一个容易踩的坑：煤耗二次项如果直接保留，目标函数就是二次约束二次规划，求解器处理起来慢。我在代码里常用的处理是把二次成本做分段线性化，比如把0到额定出力的区间切成4段，每段用线性函数逼近。这样问题退化成混合整数线性规划，Gurobi和CPLEX求解速度会快很多，而且调度结果差别很小。分段线性化的实现方法是引入分段插值变量和对应的0-1区间指示变量，Yalmip里有内置的pwl相关接口，但直接手动写线性不等式组更可控。

2.3 约束条件分层：从功率平衡到充放热逻辑

约束我习惯分成四层来组织，逐层加到模型中，能有效避免漏写或者重复冲突。

第一层是系统平衡约束。电功率平衡要求任意时段的机组总出力加上风电实际出力等于电负荷，值得提醒的是风电实际出力等于预测出力减去弃风功率。热功率平衡在上面已经写了，它是储热改造模型区别传统调度的最关键一项。

第二层是机组运行约束。包括火电机组出力上下限、爬坡约束和启停逻辑。爬坡约束要区分升爬坡和降爬坡速率，对于热电机组还要同时考虑供热变化带来的约束，如果做了简化的热电耦合线性模型，爬坡约束可以直接对电出力和热出力分别施加。

第三层是储热罐运行约束，这层最容易出错。首先容量约束要写成0到S_max的上下界；充放热功率边界要乘上对应的0-1状态变量，保证不充热时该项强制为0；同时u_ch(t) + u_dis(t) <= 1防止同时充放。其次循环约束S(0) == S(T)很关键，表示储热罐在一个调度周期结束时回到初始储能状态，否则求解器可能会把罐里热量“用光”再“偷来用”，结果虽然可行但没法实际运行。

第四层是风电场出力约束，实际出力不得大于预测出力，弃风功率为0到预测值之间的连续变量。

3. Matlab实现全流程：从变量定义到求解部署

3.1 工具箱选型与求解器配置

我用的环境是Matlab R2023a搭配Yalmip工具箱加Gurobi求解器。Yalmip的魅力在于建模语法几乎和数学公式一一对应，变量定义、约束拼接、目标函数赋值都直观得很，免去手工写大规模矩阵的折磨。Gurobi面向学术有免费license，求解混合整数线性规划的效率和稳定性在同类工具里数一数二。

没有Gurobi的时候也能跑，Matlab自带的intlinprog可以直接解MILP，代价是求解速度慢、对大算例容易卡住。CBC这类开源求解器也可以接在Yalmip后面，适合百变量量级的小模型练手。各方案对比如下：

安装Yalmip之后要跑一次yalmiptest，确认诊断信息里Gurobi或CPLEX旁边显示“found”。我遇到过太多“明明装了求解器但Yalmip提示找不到”的情况，多半是环境变量没配对或求解器版本太旧不兼容。别嫌这一步麻烦，它能为后面省大量调试时间。

3.2 决策变量与约束的Yalmip写法

假设系统里有2台热电机组、1台纯凝火电机组和1座风电场，调度周期T取24小时。先定义基础变量：

P = sdpvar(n_gen, T, 'full'); % 机组出力 H = sdpvar(n_chp, T, 'full'); % 热电机组供热出力 Q_ch = sdpvar(1, T, 'full'); % 储热罐充热功率 Q_dis = sdpvar(1, T, 'full'); % 储热罐放热功率 S = sdpvar(1, T + 1, 'full'); % 储热罐存储量，T+1对应初值 P_wind = sdpvar(1, T, 'full'); % 风电实际出力 P_curtail = sdpvar(1, T, 'full'); % 弃风功率 U = binvar(n_gen, T, 'full'); % 机组启停状态 u_ch = binvar(1, T, 'full'); % 充热状态 u_dis = binvar(1, T, 'full'); % 放热状态

这里有个细节：S的维度设成T+1而不是T，是为了把初始储热量S(0)放进变量里，循环约束可以直接写S(1) == S(T+1)，代码逻辑更干净。约束逐条添加时，我用一个Constraints = []数组不断往尾部拼接。比如功率平衡和热电耦合：

Constraints = [Constraints, ... sum(P, 1) + P_wind == P_load]; Constraints = [Constraints, ... sum(H, 1) + Q_dis - Q_ch == H_load]; Constraints = [Constraints, ... P(i, :) >= P_min(i) + alpha(i) * H(i, :)]; Constraints = [Constraints, ... P(i, :) <= P_max(i) - beta(i) * H(i, :)];

注意热电机组的供热变量H只在热电机组下标范围内定义，如果直接用矩阵运算需要小心维度匹配。爬坡约束按时间维度展开：

Constraints = [Constraints, ... P(:, 2:end) - P(:, 1:end-1) <= ramp_up(:, ones(1, T-1))]; Constraints = [Constraints, ... P(:, 1:end-1) - P(:, 2:end) <= ramp_down(:, ones(1, T-1))];

储热罐那组约束是整个模型最“有存在感”的部分：

Constraints = [Constraints, ... S(t+1) == (1 - mu) * S(t) + eta_ch * Q_ch(t) - Q_dis(t) / eta_dis]; Constraints = [Constraints, 0 <= S <= S_max]; Constraints = [Constraints, 0 <= Q_ch <= Q_ch_max * u_ch(t)]; Constraints = [Constraints, 0 <= Q_dis <= Q_dis_max * u_dis(t)]; Constraints = [Constraints, u_ch(t) + u_dis(t) <= 1];

最后加上S(1) == S(T+1)让储热罐回到初值。这一段看起来简单，但写的时候一定要逐时段检查，我曾经把sdpvar维度设错，导致约束矩阵尺寸对不上，报错信息看得人一头雾水。后来学乖了，每个变量定义后先size()检查一遍再往下写。

3.3 求解器参数整定与结果导出

调用求解器的过程用optimize函数：

ops = sdpsettings('solver', 'gurobi', ... 'verbose', 2, ... 'mipgap', 0.001, ... 'timelimit', 1800); diagnostics = optimize(Constraints, Objective, ops); if diagnostics.problem ~= 0 disp(diagnostics.info); end

mipgap设0.001还是0.01要根据算例规模权衡。小规模算例Gurobi能直接给到0.001甚至更低，但机组和时段多了以后可以适当放宽到0.01，节约的求解时间远大于误差带来的影响。timelimit是保险丝，防止算例组合爆炸后无限跑下去。

结果导出和数据分析我习惯再用一段后处理脚本：从Yalmip的变量对象中抽取value(P)，转成数值矩阵存到results结构体里，方便后续画图和对比。抽取时一定要确认求解状态是稳定可行解，否则value返回的是空值或NaN，画出来的图基本没法看。

4. 算例验证：储热改造前后的低碳经济对比

4.1 典型日数据与参数设置

为了验证模型和代码逻辑，我搭了一套小型算例：2台抽凝式热电机组，额定容量均为300MW，电出力下限为其额定值的40%；1台纯凝火电机组，容量同为300MW；1座风电场装机200MW。热负荷曲线用北方冬季典型日的形状，大致是早晚高、午间和夜间维持中等水平，峰值约180MW。储热罐容量设300MWh，最大充放热功率60MW，充放热效率均取0.9，散热损失率0.001/h。

碳价要是设得太低，模型里低碳只是个摆设；设太高又会压倒煤耗成本，出现不切实际的“疯狂压煤电”调度结果。我调了几轮，80元/吨左右在目前的算例里平衡得最好。弃风惩罚单价取600元/MWh，这个值比煤电边际成本高但又不至于高到让储热罐为了消纳风电做过度投资式的充放。

4.2 改造前后调度结果对比

我把模型设成两个对照组，方案A不接入储热罐（直接去掉储热相关变量和约束），方案B完整包含储热模型，其余参数保持一致。运行完成后，我把弃风率、总煤耗、碳排放指标摆在同一个表格里对比：

最直观的变化是弃风率从15.6%降到3.8%。原因也好解释：方案A里夜间热负荷高，热电机组被迫抬高电出力，风电场基本处于被压制状态；方案B里夜间由储热罐放热承担一部分供热，热电机组电出力可以压低，风电在夜间可以全额上网。

总煤耗成本下降大概5万元，一方面因为风电替代了部分煤电出力，另一方面储热让机组可以更稳定地运行在高效区间，没有出现频繁深调导致的高煤耗状态。碳排放下降约7%，趋势和煤耗一致，说明低碳和经济在这个场景下并不矛盾。

4.3 储热罐的工作逻辑与经济性解读

画储热罐容量曲线时能清晰看到它的昼夜循环规律：白天电负荷高、风电一般，机组发电充足，多余热量存入罐中，储热量持续上升；傍晚到夜间热负荷上升，同时风功率也大，热电机组压低电出力，储热罐放热满足热网需求，储热量持续下降。第二天早上回到初始值附近，完成一次完整循环。

这个规律和直觉完全一致，但真正跑出结果前我差点以为储热罐会“不分青红皂白地整天充放”。后来检查目标函数才发现，如果没有给充放热次数加任何惩罚，储热罐确实会出现频繁小幅动作的现象，所以我在目标函数里加的小系数运维成本是有实际意义的，它给储热罐的动作频率上了“紧箍咒”，让调度策略更贴近工程实际。

5. 实操中的坑与排查经验

5.1 Yalmip和求解器的衔接问题

Toolbox装好了、路径也添加了，但一跑optimize就提示找不到求解器，这个问题我遇到不下三次。排查顺序是这样的：先输入yalmiptest，看输出表格里Gurobi对应行是不是“found”；如果不是，到Yalmip官网重新下载最新版，覆盖旧文件；再检查系统环境变量里有没有包含Gurobi的bin目录。还有一种情况是Matlab路径里同时存在多个版本的Yalmip文件夹，新旧版本互相干扰，这种必须彻底删除一个版本才能解决。

5.2 模型不可行时的排查顺序

diagnostics.problem == 1代表模型不可行，新手往往一头雾水，其实有固定的排查路径。第一优先检查电功率平衡是否漏项，比如风电实际出力、弃风功率和负荷项是否同时存在；第二检查储热罐的循环约束和容量约束是否相容，如果S_max设置得太小，储热罐就帮不上热负荷的忙；第三检查热电耦合系数是否合理，α_i太大会导致热电机组电出力下限过高，与风电消纳目标直接冲突。

我有一个非常实用的排查技巧：要么把储热罐容量S_max放大到极大，要么把充放热功率上限设成无穷大，如果模型从不可行变成可行，问题基本锁定在这几个约束上。然后逐步缩回参数，找到出问题的最小设置点。

5.3 数值缩放和效率优化

这个问题在Matlab里比在GAMS里更突出。煤耗成本动辄几十万到上百万元，而废弃惩罚和碳价可能只有几百元，多个数量级混在一个目标函数里会让求解器的数值算法吃尽苦头。我的做法是统一量纲：功率单位用GW，成本单位用万元或百万元。比如500元/MWh改写为0.05万元/MWh乘以1000MW得到100万元，量级就不会太离谱。做完单位统一后Gurobi的求解时间和稳定性都会有肉眼可见的改善。

5.4 结果不合常理时先查这四个环节

调度结果中出现“储热罐白天疯狂放热、夜间疯狂充热”这类反直觉现象，或者风电明明可以消纳却人为弃风，优先从四个环节一一排查：弃风惩罚系数是否设得太低、碳价是否设得过高或过低、储热罐初末储能是否闭合、热负荷平衡约束的变量符号是否写反。符号方向写反是Yalmip代码里最隐蔽的错误之一，我把Q_ch和Q_dis在热平衡式里的位置写反过一次，结果储热罐变成了“供热用户”，出力的物理意义全是反的。

个人经验是，拿到一套新数据先跑一个只有24时段、1台热电机组加1个储热罐的最小系统，把逻辑理顺后再扩展到完整算例。这个习惯帮我省下的调试时间，比我在Yalmip文档上花的任何功夫都值。如果你也正在复现类似调度模型，建议从这个小系统起步，把储热罐的昼夜循环曲线调顺了，再往里面加风电、加碳价、加更多机组，每一步都有参照就不会乱。