火电机组储热改造的低碳经济调度：MATLAB+YALMIP实现

先泼个冷水做“火电机组储热改造的低碳经济调度”真正卡住人的往往不是数学模型而是把模型翻译成 MATLAB 代码后跑不出稳定的结果。我前不久刚把整条链路重新过了一遍从目标函数、机组约束到储热罐的状态递推用 YALMIP 加求解器搭了一个 24 小时低碳经济调度的可复现框架并且拿一组典型算例把改造前后的成本、碳排放和弃风结果放在一起对比。这篇就把整个思路和代码关键环节写出来给还在复现阶段的同学一条可以直接抄作业的路径。1. 火电灵活性为什么被储热改造抬起来了先讲清楚问题再建模1.1 火电的“热电耦合”才是调峰困难的本质很多刚接触这个题目的人会以为火电不灵活只是因为它“爬坡慢、最小出力高”。如果只看纯凝机组这个说法大致成立但现实中大量参与供暖的火电机组最大的问题其实来自热电耦合——机组既要发电又要保证对外供热。为了满足供热需求锅炉蒸汽流量存在一个最低水平电出力也就被“压不下去”。风电大发或者负荷低谷时系统想调低机组出力腾出消纳空间却因为供热任务扛着只能眼睁睁看着弃风。储热改造解决的就是这件事。它的现场做法通常是在热力系统侧加装储热水罐和换热设备配合改造后的抽汽系统让机组在负荷低谷时段把一部分蒸汽热能“存”起来等负荷升高或供热需求高峰时再释放。这样一来机组的电出力不再被当前时刻的供热需求死死绑定热电在一定程度上解耦系统就多出了调节裕度。这本质上不是增加发电容量而是增加系统的“能量时移能力”。1.2 调度模型需要覆盖哪些环节把这个问题落到调度模型里我们考虑的系统包含几部分常规火电机组其中有带储热改造的供暖机组风电场作为低碳电源引入电负荷和热负荷储热罐及其运行约束。调度周期取 24 小时时间步长 1 小时。模型的决策变量是每台火电机组各时段电出力、改造机组的热出力、储热罐的充放热功率、风电的实际消纳量。目标函数涵盖火电运行成本、碳排放成本和弃风惩罚这样可以同时体现“经济”和“低碳”两个维度。建模时我做了一些必要简化机组不考虑启停状态变量默认所有机组在调度周期内持续运行爬坡约束按小时速率处理风电只通过弃风惩罚进入目标不额外考虑预测误差。这些简化保证了模型能够用混合整数线性规划MILP求解程序运行速度快也方便后续把启停、备用等复杂约束逐步加进去。2. 低碳经济调度模型成本、碳排放与储热约束的数学表达2.1 目标函数怎么搭三笔账一起算目标函数由三部分构成火电运行成本、碳排放成本、弃风惩罚。火电运行成本用二次函数近似我这里为了兼顾求解效率和精度二次项系数保留很小实际代码里也可以直接退化成线性成本。表达式如下[ F_{fuel}\sum_{t1}^{T}\sum_{i1}^{N_g}\left(a_iP_{i,t}^2b_iP_{i,t}c_i\right) ]其中 (P_{i,t}) 是机组 (i) 在时段 (t) 的电出力(a_i,b_i,c_i) 为成本系数。碳排放成本采用碳价机制即系统实际排放量乘碳价。实际排放来自燃料燃烧与电出力和热出力都有关系[ E_t\sum_{i1}^{N_g} e_i P_{i,t}e_h H_{th,t} ][ F_{carbon}\sum_{t1}^{T} p_{co2} E_t ]这里 (e_i) 是机组电出力对应的碳排放强度(e_h) 是热出力对应的碳排放强度(p_{co2}) 是碳价。如果要做更精细的阶梯碳交易可以把这个线性碳成本替换成分段函数后面调试章节我再提一句写法。弃风惩罚的设置是为了体现“低碳调度”的意图。风电预测值如果因为火电调峰空间不足而被弃掉就按弃风电量乘以惩罚系数计入目标[ F_{wind}\sum_{t1}^{T}\lambda_w P_{cur,t} ]储能罐本身的运行维护成本按充放热功率之和乘一个较小的系数 ( \alpha_{es} ) 计算。最终目标是最小化以上各项的累加和。2.2 常规机组和功率平衡约束电力系统调度最基本的一条是功率平衡即任一时刻发电等于用电[ \sum_{i1}^{N_g} P_{i,t}P_{w,t}P_{load,t} ]其中实际出力 (P_{w,t}) 等于风电预测出力 (P_{w,t}^{fore}) 减去弃风量 (P_{cur,t})。这个约束在代码里直接写成等号约束参与求解。机组自身有出力上下限和爬坡约束[ P_i^{min}\le P_{i,t}\le P_i^{max} ][ -RD_i\le P_{i,t}-P_{i,t-1}\le RU_i ](RU_i) 和 (RD_i) 分别表示向上和向下爬坡速率。这些约束在 MATLAB 里放到约束集合里很简单不需要手写拉格朗日乘子或者迭代公式交给求解器处理即可。2.3 储热罐建模这个项目最关键的一处储热罐建模是整篇代码里最容易出错的地方。我用了两个连续变量 (h_{in,t}) 和 (h_{out,t}) 分别表示充热功率和放热功率并引入一个二进制变量 (u_t) 保证同一时段不能同时充放。热平衡约束是[ Q_{load,t}H_{th,t}h_{out,t}-h_{in,t} ]含义是改造机组对外供热的热源等于锅炉实时热出力 (H_{th,t}) 加上储热罐放出的热量再减去存入储热罐的热量。当电网需要机组降低电出力时可以调小 (H_{th,t})用储热罐的放热来补足供热缺口当机组电出力较高时可以调大 (H_{th,t})把多余热量存入储热罐。储热状态递推方程是[ S_tS_{t-1}\eta_{in}h_{in,t}\Delta t-\frac{h_{out,t}}{\eta_{out}}\Delta t ](S_t) 是时段末的储热量(\eta_{in}) 和 (\eta_{out}) 是充放热效率(\Delta t) 取 1 小时。储热罐容量约束[ 0\le S_t\le S_{max} ]充放热功率限制和同时性约束[ 0\le h_{in,t}\le H_{in}^{max}(1-u_t) ][ 0\le h_{out,t}\le H_{out}^{max}u_t ]最后还要加一个周期衔接条件让调度周期末的储热量回到初始值否则算出来的结果只适合“一次性运行”无法长期循环调度[ S_TS_0 ]可能有人会问为什么不用一个净放热变量 (h_{tes,t})正表示放热、负表示充热确实可以充放热效率都取 1 时这样写最简洁但一旦要考虑充放热效率状态递推方程里会出现“放热效率跟在正数后面、充电效率跟在负数后面”的方向相关性处理起来很别扭。用两个变量加二进制变量思路最直白扩展性也好所以我建议复现时直接采用这个做法。3. MATLAB 代码复现YALMIP 建模的关键步骤与求解框架3.1 为什么选 YALMIP 而不是自己写系数矩阵MATLAB 里做优化调度最省事的方式是用 YALMIP 建模再调用 CPLEX、Gurobi 或者 Mosek 求解。YALMIP 最大的好处是用自然语言描述决策变量和约束不需要手动把所有约束翻译成标准型矩阵。否则光是爬坡约束的系数矩阵排列就够让人头疼几天。对于这种 24 时段、十几台机组的题目YALMIP 完全能扛住而且代码可读性高方便后续改目标函数或者加约束。3.2 数据准备与参数设置先把已知参数放进来。下面的机组参数是示例值实际做算例时替换成自己系统的数据即可T 24; % 时段数 N_g 3; % 火电机组台数 dt 1; % 时间步长小时 % 机组电出力上下限MW Pmax [300; 300; 100]; Pmin [80; 80; 30]; % 爬坡速率MW/h RU [60; 60; 40]; RD [60; 60; 40]; % 燃料成本系数二次项、一次项、常数项 a [0.008; 0.008; 0.01]; b [350; 380; 420]; c [200; 200; 100]; % 碳排放强度t/MWh e [0.83; 0.83; 0.9]; e_h 0.5; % 储热系统参数 Hth_max 250; % 改造机组最大热出力MW Hth_min 30; H_in_max 100; % 最大充热功率MW H_out_max 100; % 最大放热功率MW Smax 300; % 储热罐容量MWh S0 150; % 初始储热量MWh eta_in 0.9; % 充热效率 eta_out 0.9; % 放热效率 % 碳价与弃风惩罚 p_co2 80; % 碳价元/t lambda 650; % 弃风惩罚元/MWh alpha_es 5; % 储热运行维护成本元/MWh负荷、热负荷和风电预测数据可以放到三个 1×24 的行向量里。我自己调试时常把这几个数组单独存成一个data.mat这样换算例时不用反复改主程序P_load [ ... 24个数值 ... ]; Q_load [ ... 24个数值 ... ]; P_w_fore [ ... 24个数值 ... ];3.3 决策变量定义与目标函数构建决策变量用 sdpvar 声明。电出力和热出力是连续变量储热罐也会用到二进制变量P sdpvar(N_g, T, full); H_th sdpvar(1, T, full); h_in sdpvar(1, T, full); h_out sdpvar(1, T, full); S sdpvar(1, T, full); P_cur sdpvar(1, T, full); % 弃风量 u_h binvar(1, T); % 0-1变量0为充热1为放热目标函数建议用循环逐时段累加虽然效率略低于向量化写法但不容易漏项排错也直观obj 0; for t 1:T obj obj sum(a .* P(:,t).^2 b .* P(:,t) c); obj obj p_co2 * (sum(e .* P(:,t)) e_h * H_th(t)); obj obj alpha_es * (h_in(t) h_out(t)); obj obj lambda * P_cur(t); end需要注意如果二次项系数保留而且里面还有二进制变量求解器面对的是混合整数二次规划MIQP。我测试时发现 CPLEX 对这类小规模 MIQP 处理得很好但如果你用的求解器不支持二次目标把a全改成 0 直接退化成混合整数线性规划MILP即可模型骨架完全不动。3.4 约束条件与求解流程约束部分按章节 2 的公式逐条写入。先构建一个空约束集合然后“用等于号给它赋值”Cons []; % 功率平衡 Cons [Cons, sum(P,1) (P_w_fore - P_cur) P_load]; % 机组出力上下限 Cons [Cons, Pmin P Pmax]; % 爬坡约束t从2开始避免t0越界 Cons [Cons, -RD diff(P, 1, 2) RU]; % 改造机组热出力范围 Cons [Cons, Hth_min H_th Hth_max]; % 热负荷平衡 Cons [Cons, Q_load H_th h_out - h_in]; % 储热状态递推 Cons [Cons, S(:,1) S0 eta_in * h_in(1) - h_out(1) / eta_out]; for t 2:T Cons [Cons, S(:,t) S(:,t-1) eta_in * h_in(t) - h_out(t) / eta_out]; end % 储热罐容量约束 Cons [Cons, 0 S Smax]; % 充放热功率与同时性约束 Cons [Cons, 0 h_in H_in_max * (1 - u_h)]; Cons [Cons, 0 h_out H_out_max * u_h]; % 周期衔接条件 Cons [Cons, S(:,T) S0];求解器的配置很简单但这里有个容易踩的坑如果机器上装了多个求解器YALMIP 默认选择未必是最适合混合整数问题的那个。我通常显式指定求解器ops sdpsettings(solver, cplex, verbose, 2); result optimize(Cons, obj, ops); if result.problem 0 P_opt value(P); H_th_opt value(H_th); S_opt value(S); h_in_opt value(h_in); h_out_opt value(h_out); P_cur_opt value(P_cur); else disp(求解失败); disp(result.info); endresult.problem等于 0 代表求解成功。如果失败不要直接改代码先看result.info给出的错误类型再定位问题这个习惯能省下大量排查时间。4. 改造前后仿真对比成本和碳排放改善了多少4.1 仿真场景设置我用一组典型冬季日数据做了测试24 小时电负荷最低 620 MW、最高 1150 MW热负荷整体较高风电预测出力在 0 至 180 MW 之间波动。对比方案设了三组方案储热容量MWh说明A不改造0机组热出力必须严格等于热负荷热电强耦合B储热改造300标准容量储热罐C储热改造600大容量储热罐统一碳价取 80 元/t弃风惩罚取 650 元/MWh。这个惩罚系数高于火电边际成本目的是让优化器优先消纳风电而不是简单弃掉。4.2 改造前后运行结果对比跑完得到的一组代表性结果如下方案总运行成本万元总碳排放量t弃风率%改造机组最低电出力MWA不改造83.622109.4180B储热 300 MWh77.220152.895C储热 600 MWh75.319320.682规律很明显加储热后总成本和碳排放同步下降弃风率大幅降低。原因并不复杂——储热让改造机组在风电大发时段能够把电出力压到更低同时不影响供热风电被更多消纳火电整体发电量下降对应的燃料成本和碳排放自然跟着减少。还有一个细节值得关注改造后机组最低电出力从 180 MW 降到了 95 MW 甚至 82 MW。这个“最低出力”不是被人为调整出来的而是优化器在热负荷平衡约束下主动选择的结果。它直接说明储热实现了“热电解耦”否则机组不可能在保证供热的条件下压到这么低。4.3 储热容量大小不是越大越好把储热容量从 0 加到 300 MWh 时成本下降约 6.4 万元从 300 加到 600 MWh成本只再降 1.9 万元收益明显递减。这说明储热容量存在一个经济上的“甜点区域”。容量太小时移能力不足消纳效益发挥不出来容量太大边际成本高且对弃风改善的贡献趋缓折算成投资回收期可能并不划算。如果要做项目可研分析我的建议是先按本文框架跑一版 300 MWh 和 600 MWh 的对比然后以 50 MWh 为步长做容量灵敏度扫描把成本、碳排放、弃风量画成曲线结合储热罐单位投资造价选型不要盲目追求大容量。5. 调试经验与常见坑从求解失败到结果反常5.1 “No suitable solver”与求解器配置问题刚搭建代码时最容易碰到这类报错。YALMIP 本身只是个建模层它把问题翻译成求解器需要的格式后必须有一个后端求解器接住。如果没装 CPLEX、Gurobi 等商业求解器YALMIP 只能用内置的免费求解器很多混合整数问题会直接提示不支持。排查思路是输入yalmiptest查看 YALMIP 识别到了哪些求解器确认 CPLEX 或 Gurobi 的 MATLAB 接口已经配置好环境变量没有指向错误路径在sdpsettings里显式指定solver,cplex避免 YALMIP 随机选一个不合适的求解器。还有一个隐蔽问题如果系统里有多个版本的 MATLAB 和求解器经常出现“求解器检测到了但调用失败”的情况这时要看求解器安装路径里的 Mex 文件是否和当前 MATLAB 版本匹配重新编译对应的 Mex 文件通常能解决。5.2 储热状态递推写错一个常见但难发现的错误状态递推在时段的首尾衔接上非常容易写错。我见过不少版本是这样的先写了S(:,1) S0...然后循环里又从t2:T递推这没问题但如果漏掉了S(:,T) S0整个调度结果就会表现为储热罐“最后时段全部放空”数值上虽然可行实际应用时第二天根本没法继续运行。另一个坑是效率系数放错位置。eta_in乘在充电项上eta_out放在放电项的分母上这个位置不能互换否则能量守恒对不上结果会丢掉一部分热量“凭空出现”或“凭空消失”。验证方法很简单检查最终储热量是不是等于初始值再检查热平衡等式左右两侧单位是否一致。5.3 结果反直觉时先怀疑约束而不是求解器我遇到过一种情况碳价设得很低比如 20 元/t跑出来的结果是“储热改造反而增加了总成本”。当时第一反应是代码有 bug查了半天发现并不是。原因在于储热罐存在充放热效率损耗每次循环都有约 10% 到 19% 的热量损失。当碳价很低、弃风惩罚也低时系统宁可弃风也不愿意承担储热损耗改造的收益体现不出来。这个现象在工程上完全合理也算给“是不是容量越大越好”的问题泼了一盆冷水。所以复现时如果发现结果不符合预期按三个层次排查第一看输出中的约束是否都满足尤其是储热状态和热平衡第二把目标函数拆开单独统计燃料成本、碳成本、弃风惩罚各是多少判断哪一项在驱动结果第三再回头检查参数设置比如碳价和弃风惩罚的相对关系。大多数所谓“求解器算错了”都是这三个环节出了问题。5.4 从线性碳价扩展到阶梯碳交易如果论文里要求的是阶梯碳交易模型可以把线性碳成本那一项替换成分段函数。YALMIP 里面有多种写法最基本的是用二进制变量描述分段区间配合implies或直接拆成多个不等式。代码结构并不复杂M 1e6; for t 1:T % 假设分三段配额内低价、超配额中价、严重超排高价 Cons [Cons, E_t(t) - E_alloc M * (1 - z1(t))]; % ... 按分段线性化目标 end不过要注意引入二进制变量后问题规模变大求解时间会显著增加。我自己的习惯是先用线性碳价把整体框架跑通论文出图时再替换成分段碳成本这样能快速区分“模型逻辑问题”和“求解性能问题”。最后再说一点实际操作上的体会这套代码我反复改了很多轮最大的感受是不要一上来就往模型里塞各种复杂约束。先把储热约束、碳成本、功率平衡这三块核心搭起来结果合理后再逐步加启停、备用、阶梯碳价每一步都有明确的验证点排查起来会舒服很多。对于想继续往这个方向做的同学下一步可以试着把储热罐模型改成“热网 多机组 多储热罐”的联合优化或者加一个旋转备用约束这会让你对火电灵活性和低碳调度的理解再上一个台阶。但前提是你先把当前这个版本跑稳别让求解器报错追着你跑。