做电力系统优化调度的人估计对“以热定电”这四个字都不陌生。北方冬季热电联产机组一开发电出力被供热需求绑得死死的电负荷低的时候也没法降太多风电只能在夹缝里求生存。“考虑火电机组储热改造的电力系统低碳经济调度”这个题目说白了就是给火电机组加一个储热罐把热和电解耦然后在满足负荷的前提下让系统的碳排放和总运行成本一起降下来。我最近把这个课题从建模到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后面适合百变量量级的小模型练手。各方案对比如下求解器许可证速度表现适用场景Gurobi学术免费/商业付费快千变量秒级论文级算例、工业规模CPLEX学术免费快与Gurobi接近企业项目、教学研究CBC开源免费中等几百变量尚可入门验证、个人学习intlinprogMatlab内置偏慢小规模跑通流程安装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); % 储热罐存储量T1对应初值 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的维度设成T1而不是T是为了把初始储热量S(0)放进变量里循环约束可以直接写S(1) S(T1)代码逻辑更干净。约束逐条添加时我用一个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(t1) (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(T1)让储热罐回到初值。这一段看起来简单但写的时候一定要逐时段检查我曾经把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); endmipgap设0.001还是0.01要根据算例规模权衡。小规模算例Gurobi能直接给到0.001甚至更低但机组和时段多了以后可以适当放宽到0.01节约的求解时间远大于误差带来的影响。timelimit是保险丝防止算例组合爆炸后无限跑下去。结果导出和数据分析我习惯再用一段后处理脚本从Yalmip的变量对象中抽取value(P)转成数值矩阵存到results结构体里方便后续画图和对比。抽取时一定要确认求解状态是稳定可行解否则value返回的是空值或NaN画出来的图基本没法看。4. 算例验证储热改造前后的低碳经济对比4.1 典型日数据与参数设置为了验证模型和代码逻辑我搭了一套小型算例2台抽凝式热电机组额定容量均为300MW电出力下限为其额定值的40%1台纯凝火电机组容量同为300MW1座风电场装机200MW。热负荷曲线用北方冬季典型日的形状大致是早晚高、午间和夜间维持中等水平峰值约180MW。储热罐容量设300MWh最大充放热功率60MW充放热效率均取0.9散热损失率0.001/h。碳价要是设得太低模型里低碳只是个摆设设太高又会压倒煤耗成本出现不切实际的“疯狂压煤电”调度结果。我调了几轮80元/吨左右在目前的算例里平衡得最好。弃风惩罚单价取600元/MWh这个值比煤电边际成本高但又不至于高到让储热罐为了消纳风电做过度投资式的充放。4.2 改造前后调度结果对比我把模型设成两个对照组方案A不接入储热罐直接去掉储热相关变量和约束方案B完整包含储热模型其余参数保持一致。运行完成后我把弃风率、总煤耗、碳排放指标摆在同一个表格里对比指标方案A无储热方案B有储热弃风率15.6%3.8%总煤耗成本万元86.482.1碳排放吨39203650储热罐日均循环次数-1.2最直观的变化是弃风率从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文档上花的任何功夫都值。如果你也正在复现类似调度模型建议从这个小系统起步把储热罐的昼夜循环曲线调顺了再往里面加风电、加碳价、加更多机组每一步都有参照就不会乱。