前一阵在一个电力系统方向的交流群里有人甩出这个标题“考虑火电机组储热改造的电力系统低碳经济调度Matlab代码实现”求指导。这其实不是个冷门问题——北方供热季热电联产机组被“以热定电”捆绑得死死的时候新能源消纳空间被挤占弃风弃光成了常态。储热改造就是把这个绑定的游戏规则打散热负荷由储热罐和机组共同承担电出力有了腾挪余地碳排放和运行成本才能一并压下来。这篇文章我就用Matlab代码实现为主线完整拆解这类低碳经济调度问题的建模思路、代码架构和实操坑点适合正在做电力系统优化调度、火电灵活性改造、储能/储热配置方向的研究生和工程技术人员参考。这个项目的本质是在一个典型的电热联合系统里把储热装置的充放热过程写成约束放进以“煤耗成本 碳排放成本 弃风惩罚”为目标函数的经济调度模型里用求解器算出未来24小时各机组的电出力和热出力计划。下面对整个建模、编码、调试过程做一次系统复盘。1. 项目背景与问题拆解1.1 供热季的“以热定电”困境先说我为什么觉得这个题目值得认真做。北方冬季供暖期热电联产机组CHP为了保证居民供热电出力必须跟着热负荷走这个运行模式叫“以热定电”。问题在于热负荷越高机组的最小电出力就被抬得越高。举个直观例子一台300MW的抽汽式供热机组纯凝工况下最小技术出力可能能压到40%额定容量但进入供热工况后最小电出力可能被顶到70%以上。电网调度留给风电、光伏的消纳空间自然被压缩弃风也就避不开。这个场景里的“低碳经济调度”不是简单的环保口号而是要同时算清三笔账煤耗成本随出力变化是二次曲线、碳排放成本跟煤耗直接挂钩、弃风惩罚体现新能源消纳优先级。这三笔账合在一个目标函数里求最小化才是“低碳”和“经济”同时成立的调度方案。1.2 储热改造如何破局火电机组储热改造的关键是在热力系统里加一个“缓冲池”——储热罐。它做的事本质上和电池在电力系统里的角色一样热负荷高峰时放热配合机组共同供热热负荷低谷时充热把机组多余的热出力存起来。这样机组发电与供热之间不再是一一绑定的刚关系电出力下限有了向下调整的空间风电就能多发一些。在调度模型里储热罐真正参与进来需要三类约束一是储热罐自身的热量平衡SOC随充放热功率变化的连续约束二是充放热功率上下限约束三是储热容量上下限约束。这些约束写进优化模型后求解器会根据全天电、热负荷曲线和各机组运行成本自动决定储热罐什么时候充、什么时候放、充多少、放多少。1.3 为什么强调“低碳”而非单纯“经济”早年做经济调度只算煤耗成本现在碳排放成本成了一项必须纳入的硬指标。常见的做法有两种一种是对机组单位发电量设置碳排放配额超出部分去碳市场购买产生碳排放成本另一种是直接引入碳税系数。从建模角度前者更贴近国内碳交易实际模型里表现为分段线性或阶梯型碳成本后者更简单一个线性系数就能完事。我这份代码里用的是“碳排放配额超出部分按碳价购买”的模式本质是把碳排放成本写成煤耗量的函数。由于煤耗量本身是出力的二次函数碳排放成本也就继承了这个非线性好在YalmipGurobi/Cplex这类求解器处理这类问题非常成熟。2. 数学模型搭建的核心逻辑2.1 目标函数三笔账要一起算目标函数是整个调度的指挥棒。我用的目标函数由四个部分构成煤耗成本各机组电出力、热出力对应燃煤量的二次函数之和。碳排放成本燃煤定量排放的CO₂总量减去免费配额后的部分乘以碳交易价格。弃风惩罚成本风电预测出力与实际调度出力之差乘以惩罚系数。机组启停成本可选若做含整数变量的机组组合则计入启停成本若只做经济调度可不考虑视问题规模而定。模型形式如下objective sum(sum(a_i .* P_i.^2 b_i .* P_i c_i)) ... carbon_price * sum(sum(E_i)) ... penalty_wind * sum(sum(P_wind_pred - P_wind_use));这里有三个关键的“为什么”。第一煤耗用二次函数而不用线性函数是因为实际机组热耗率随负荷率变化明显线性模型在深度调峰区间误差过大。第二弃风惩罚系数必须明显高于煤耗成本单价否则求解器为了省煤耗会主动弃风。第三碳排放成本如果和煤耗成本同源可以合并简化但分开写的好处是方便做碳价灵敏度分析——碳价从50元/吨涨到200元/吨时调度策略怎么变这是政策分析常用的场景。2.2 电力与热力平衡约束系统级约束是每一时段都必须满足的硬约束。电力平衡约束Constraints [Constraints, sum(P(:, t)) P_wind_use(t) P_load(t)];热力平衡约束Constraints [Constraints, sum(H(:, t)) H_tes_discharge(t) H_load(t) ... H_tes_charge(t)];注意热力平衡里储热罐的符号约定蓄热时相当于热负荷增加放热时相当于热源出力增加。很多初学者在这里容易把符号搞反导致模型得出完全相反的调度结果。风电出力上下限约束为0到预测出力之间实际调度出力不允许超过预测值也不需要低于0Constraints [Constraints, 0 P_wind_use(t) P_wind_pred(t)];2.3 储热罐SOC建模储热罐的核心约束是热量连续方程。设SOC(t)为t时段结束时储热罐的热量状态则有Constraints [Constraints, SOC(t1) SOC(t) ... H_tes_charge(t) * eta_ch - H_tes_discharge(t) / eta_dis];对应代码里的循环写法for t 1:T-1 Constraints [Constraints, SOC(t1) SOC(t) ... H_tes_charge(t) * eta_ch - H_tes_discharge(t) / eta_dis]; end这里两个效率系数充热效率η_ch、放热效率η_dis很容易被忽略但它们对结果影响不小。假设充放效率都是0.9一充一放合计损耗19%的热量模型里如果不写效率约束等于白送19%的可用热量优化结果会偏乐观实际运行根本无法兑现。SOC还需要容量约束和初始值设定。初始值按调度周期开始时储热罐的热量状态设定为了让调度具有周期性一般会让结束SOC回到初始值附近我给终端SOC设定了一个等式约束或软约束Constraints [Constraints, SOC(1) SOC_init]; Constraints [Constraints, SOC(T) SOC_init]; % 或允许一定偏差这个“周期性边界”条件很重要。不设这个约束求解器会在最后一个时段把储热罐热量放空前一天的蓄热策略全部失真强制回到初始值才能得到可重复执行的调度计划。充放热功率上下限约束和容量约束Constraints [Constraints, 0 H_tes_charge(t) H_ch_max]; Constraints [Constraints, 0 H_tes_discharge(t) H_dis_max]; Constraints [Constraints, SOC_min SOC(t) SOC_max];2.4 机组物理约束每台机组都有出力上下限和爬坡约束。如果题目只做经济调度而不做机组组合机组的启停状态是固定的那么约束是连续线性的Constraints [Constraints, P_min(i) P(i, t) P_max(i)]; Constraints [Constraints, H_min(i) H(i, t) H_max(i)]; Constraints [Constraints, -ramp_down(i) P(i, t1) - P(i, t) ramp_up(i)];如果要把机组启停也作为决策变量那就需要引入二进制变量z(i,t)表示机组开机状态约束变成Constraints [Constraints, P(i,t) P_max(i) * z(i,t)]; Constraints [Constraints, P(i,t) P_min(i) * z(i,t)];这样做模型就从LP变成MILP求解时间会显著增加。对于含储热改造的中等规模算例6台机组、24时段、单储热罐即使全部连续化处理模型变量数通常在300到500之间Gurobi秒级求解一旦引入启停二进制变量变量数翻倍还不止求解时间可能到几十秒级别。建议先跑连续模型验证逻辑正确性再加二进制变量做扩展。热电联产机组还有一个特有的约束——电热运行区间。抽汽式机组电出力与热出力之间存在耦合关系简单描述为电出力下限随热出力增大而增大电出力上限随热出力增大而略微降低。严格的可行域是个多边形区域处理方式是引入线性不等式组。我在代码里用了一个简化但有效的写法% 抽汽式机组电热耦合约束 Constraints [Constraints, P(i,t) P_min_pure(i) k_coupling(i) * H(i,t)]; Constraints [Constraints, P(i,t) P_max_pure(i) - k_hb(i) * H(i,t)];其中k_coupling和k_hb是根据机组热工特性拟合的系数这两个系数的值来自机组实际运行工况数据不是随便拍的。做项目时一定要找到目标机组的热力特性曲线哪怕是用典型参数也行但不要拿来路不明的系数直接用。3. Matlab代码实现架构3.1 为什么选MatlabYalmip我见过用PythonPyomo做同样问题的也见过直接调Gurobi C API的但只要还是学术研究和快速验证阶段MatlabYalmip的组合确实最顺手。原因有三点第一Yalmip的约束建模语法极其接近数学表达式写出MPC式或经济调度式的约束几乎零翻译成本。上面那些约束写法几乎可以直接从草稿纸抄到代码里。第二Matlab的矩阵运算和结果可视化一体调度曲线、SOC曲线、弃风量这些结果用一个plot就能画出来不需要额外对接绘图库。第三求解器切换方便。同一套模型Yalmip底层可以接Gurobi、Cplex、Mosek也可以接开源的SCIP模型文件不动只需修改sdpsettings里的solver字段。这一点在后面排查模型问题上非常省事。3.2 变量定义与数据准备先定义时间尺度和机组参数。核心数据结构尽量用矩阵P是n×T矩阵、H是n×T矩阵这样后面写约束可以整行整列操作避免大量for循环包裹导致的代码冗长。T 24; % 调度周期24小时 n_unit 6; % 机组台数 n_wind 1; % 风电场数量 % 机组参数 P_max [300 300 200 200 150 150]; % 电出力上限 MW P_min [100 100 60 60 40 40]; % 电出力下限 MW H_max [250 250 180 180 100 100]; % 热出力上限 MWth ramp_up [60 60 40 40 30 30]; % 爬坡上限 MW/h ramp_down [60 60 40 40 30 30]; % 爬坡下限 MW/h % 煤耗系数 a*x^2 b*x c a_coal [0.002 0.002 0.003 0.003 0.004 0.004]; b_coal [20 20 25 25 30 30]; c_coal [100 100 80 80 50 50]; % 碳排放参数 coal_to_co2 2.6; % 每吨标煤燃烧排放CO2量 carbon_price 100; % 碳交易价格 元/吨 carbon_quota 1.0; % 每MWh电量免费配额 吨/MWh风电预测出力、电负荷、热负荷这三条曲线是整个调度问题输入逻辑的地基。实际项目中这三条曲线来自历史数据或预测系统做代码演示时可以用典型日曲线P_wind_pred [120 110 100 95 80 70 65 75 90 110 130 150 ... 160 155 140 120 100 85 80 90 95 100 110 120]; P_load [220 210 200 195 200 210 230 260 280 290 285 275 ... 270 265 275 285 290 295 300 305 310 300 270 240]; H_load [170 165 160 155 150 148 150 155 180 190 195 200 ... 210 215 220 225 230 235 240 245 240 230 210 185];注意这些示例曲线的量级要和机组参数匹配如果实际算例中负荷远大于机组容量之和模型必然不可行这种问题别急着调求解器先回头检查数据量纲和数量级。然后定义Yalmip变量P sdpvar(n_unit, T, full); % 各机组电出力 MW H sdpvar(n_unit, T, full); % 各机组热出力 MWth P_wind_use sdpvar(1, T, full); % 风电实际出力 MW H_tes_ch sdpvar(1, T, full); % 储热罐充热功率 MWth H_tes_dis sdpvar(1, T, full); % 储热罐放热功率 MWth SOC sdpvar(1, T1, full); % 储热罐热量状态 MWhth这里SOC定义T1个点是为了让约束里出现SOC(t)和SOC(t1)时索引不越界同时方便最后画SOC曲线时对应到24个时段末的状态。3.3 约束生成与求解配置约束组装用Yalmip的约束拼接方式把系统平衡、机组约束、储热罐约束依次写进去。完整的约束组装大致是Constraints []; % 1. 电力平衡 for t 1:T Constraints [Constraints, sum(P(:,t)) P_wind_use(t) P_load(t)]; end % 2. 热力平衡 for t 1:T Constraints [Constraints, sum(H(:,t)) H_tes_dis(t) H_load(t) H_tes_ch(t)]; end % 3. 机组出力上下限与爬坡 for i 1:n_unit for t 1:T Constraints [Constraints, P_min(i) P(i,t) P_max(i)]; Constraints [Constraints, H_min(i) H(i,t) H_max(i)]; end for t 1:T-1 Constraints [Constraints, -ramp_down(i) P(i,t1)-P(i,t) ramp_up(i)]; end end % 4. 储热罐约束 for t 1:T Constraints [Constraints, 0 H_tes_ch(t) H_ch_max]; Constraints [Constraints, 0 H_tes_dis(t) H_dis_max]; Constraints [Constraints, SOC_min SOC(t) SOC_max]; end for t 1:T Constraints [Constraints, SOC(t1) SOC(t) ... eta_ch * H_tes_ch(t) - H_tes_dis(t) / eta_dis]; end Constraints [Constraints, SOC(1) SOC_init]; Constraints [Constraints, SOC(T1) SOC_init]; % 周期边界目标函数里有一个技术细节要强调煤耗成本是二次函数写成sdpvar的数组二次项后Yalmip自身会识别二次凸问题并调用求解器的QP或SOCP能力。但如果后面加上机组启停二进制变量二次项加整数变量变成MIQPGurobi处理得不错Cplex稍慢一些。为了求稳先跑“二次目标 连续变量”的QP模型再考虑是否升到MIQP。求解设置options sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.001, ... gurobi.TimeLimit, 300); sol optimize(Constraints, objective, options);如果机器上没装Gurobi也可以直接指定用Cplex或Mosek。全连续模型我实测Gurobi 9.x通常1秒内出解MILP模型在6机组24时段规模下一般也在10秒量级越界了就要怀疑模型结构有没有问题。3.4 结果后处理与可视化求解完成后用value()取出各变量数值P_opt value(P); H_opt value(H); P_wind_opt value(P_wind_use); SOC_opt value(SOC);画图方面至少输出三张图电功率平衡堆叠图、热功率平衡堆叠图、储热罐SOC曲线。堆叠图用area指令比普通plot更清楚展示各机组在总出力中的占比figure; area(P_opt); hold on; plot(P_load, k-, LineWidth, 2); legend(arrayfun((i) sprintf(机组%d, i), 1:n_unit, UniformOutput, false), ... 电负荷, Location, best); xlabel(时段/h); ylabel(功率/MW); title(电功率平衡);SOC曲线要单独画因为一段时间里充放电导致的SOC变化轨迹是判断储热罐调度策略是否合理的直接依据。如果SOC曲线频繁顶到上限或碰触下限说明储热罐容量选得不合适或者充放功率限制过紧这一点指导工程配置很有价值。4. 算例设计与结果分析4.1 算例参数设置我用来做验证的算例是6台火电机组其中3台抽汽式热电联产机组、3台纯凝机组、1个风电场、1座储热罐。24小时调度周期时间步长1小时。储热罐参数容量450 MWhth最大充放热功率100 MWth充热效率0.92放热效率0.92初始SOC225 MWhth半满状态碳交易价格取100元/吨碳排放配额按每MWh电量0.5吨CO₂设定风电惩罚系数取200元/MWh。这里惩罚系数的取值直接决定弃风是否发生取值过低时模型宁可弃风也不调高成本机组出力算例效果就不明显。4.2 结果解读弃风率怎么降下来的基线场景是无储热改造所有热负荷全部由热电联产机组承担。此时在负荷低谷时段夜间风电多发但热负荷高机组电出力被顶在高位风电即使预测出力很高也无法全额消纳弃风率大约在18%左右。加入储热改造后凌晨时段储热罐自动转入充热模式热电联产机组的热出力一部分被储热罐“吸收”电出力下限得以下调风电并网空间增大。到白天电负荷高峰时段储热罐放热协助供热机组把更多发电能力释放给电负荷。对比结果场景煤耗成本元碳排放成本元弃风惩罚元弃风率无储热412000386002680018.2%含储热3985003610082006.5%储热改造后总运行成本下降约7.6%弃风率从18.2%降到6.5%。这个结果在趋势上符合预期储热罐并没有直接“创造”能量它只是把热负荷在时间轴上做了平移让机组出力分配更优从而释放风电空间、降低总煤耗。4.3 灵敏度分析储热罐容量多大才划算我先固定其他参数把储热罐容量从0扫到800MWhth每次间隔50MWhth看看总成本随容量变化的曲线。结果是0到300MWhth区间总成本下降非常明显边际收益大300到550MWhth区间成本继续下降但边际收益趋缓超过550MWhth后总成本几乎不再变化因为最大充放热功率100MW已经限制了储热罐“吃的速度”容量再大也没有足够时段去利用。这个结论对工程很有意义储热罐不是越大越好容量和最大充放热功率要匹配。容量超出需求后增加的只是投资成本调度收益却封顶了。工程上应该在给定最大充放热功率条件下用模型做容量-收益曲线找“膝盖点”确定推荐容量。灵敏度分析还有一个不错的方向是扫碳价。把碳价从0调到300元/吨观察调度结果低价时风电消纳优先度一般高价时高热电联产机组会尽量压低出力弃风率随之下降碳成本在总成本中占比上升。把这张曲线拿出来可以直接给政策层面的碳价机制设计提供数据支撑。5. 常见问题与排查技巧5.1 模型不可行的三个高频原因这个项目我前前后后调试过很多次不可行问题出现得最多而且几乎都是数据问题而非模型问题。第一个高频原因是电负荷在某个时段小于所有机组最小出力之和此时电力平衡约束无解。这种情况在纯凝机组占比高的时候特别明显夜间负荷低谷机组最小出力压不住弃电约束又限制了消纳通道。排查方法是先算每个时段负荷与机组最小出力之和的差值曲线看是否出现负值。第二个高频原因是热力平衡约束与储热罐容量不匹配。某种极端场景下热负荷过高而储热罐容量和放热功率都太小导致热出力上限之和加放热功率仍然低于热负荷。这种问题通常同时表现为SOC下限越限。第三个原因比较隐蔽是周期边界约束导致的无解。当SOC(t1) SOC(t)这个等式约束和“T时刻必须回到初始值”约束一起出现时如果充放效率不对称模型可能在物理上无法完成一个完整循环导致无解。我解决的办法是把末尾SOC等式的等式约束改成软约束允许在初始值附近±5%波动目标函数里加一个小惩罚项模型立刻收敛。Constraints [Constraints, SOC(T1) 0.95 * SOC_init]; Constraints [Constraints, SOC(T1) 1.05 * SOC_init]; objective objective 10 * (SOC(T1) - SOC_init)^2;5.2 求解慢的实用排查清单模型规模不大却求解慢多半是模型问题。我整理了一个排查顺序检查是否不小心生成了二进制变量。Yalmip里只要约束中含binvar或integer模型就是MILP/MIQP。先用yalmiptest或直接看求解器输出日志里的变量类型统计。检查目标函数是否非凸。二次项的系数矩阵如果是非半正定求解器会退化到用非线性算法来解一个本应是凸的问题速度会差两个数量级。检查是否有冗余的强耦合约束。某些资料里会把爬坡约束同时加到电出力和热出力上实际上热出力通常不参与爬坡加上后约束矩阵变得更稠密求解变慢但没有精度收益。检查是否有数量级差异极大的参数。比如煤耗成本系数是10⁴量级而SOC约束系数是10⁰量级Gurobi的预处理虽然能处理但数值病态问题会拖慢收敛。统一量纲到MW/MWh后有明显改善。如果模型确认是线性的QP或MILP但求解时间还是超过预期先试着给求解器增加MIPGap容忍度。实际工程里1%的gap结果完全可用但硬要达到0.01%最优性证明求解时间可能飙升10倍。5.3 量纲与数值陷阱有一次我算出来的SOC曲线出现了负值怎么查都查不出约束逻辑错误最后发现问题出在单位上负荷数据是用万千瓦万kW写的机组参数用的是MW两者差了10倍导致电力平衡约束从数学上就不成立优化器只能给出荒谬的解。这种事太容易发生了。我的习惯是所有输入数据进模型之前全部统一成MW和MWhth负荷曲线如果原始单位是万kW就除以10再赋值。这个转换直接用代码做别在脑子里换算P_load P_load_raw / 10; % 万kW - MW另外我这里用的效率系数写法是充放分开很多文献里用一个总效率η替代。如果总效率是0.9实际物理过程是充放各损耗一部分直接写成SOC(t1) SOC(t) eta * H_tes_ch(t) - H_tes_dis(t) / eta并不是严格意义上的“总效率”而是充放效率相同的一种特殊情况。两种写法在数值上有细微差别但对结论影响不大关键是不要混合使用两种写法。6. 这个代码还能怎么扩展6.1 从确定性到鲁棒/随机优化目前的模型是确定性的——负荷、风电预测都是给定曲线。实际调度中预测误差逃不掉尤其是风电功率预测误差在24小时尺度上非常明显。扩展方向有两个一是场景法把历史预测误差的多个场景输入模型求期望成本最小的二阶段随机优化二是鲁棒优化把风电出力写成区间不确定集合求最坏情况下的最优解。这两个方向都需要对当前确定性模型做较大改造。随机优化要加场景索引变量维度从T扩展到S×T模型规模线性增长鲁棒优化则要把部分约束改成对偶形式或者引入鲁棒对等约束。建议在确定性模型调通后再做扩展别一上来就挑战不确定性模型——数学性质完全不同排查问题的难度也翻倍。6.2 碳交易价格的敏感性分析我在4.3已经展示了碳价扫描结果这里补充一个实操层面的细节做碳价灵敏度分析时没必要对每个碳价都跑一遍优化然后手动记录结果直接在外层写一个循环脚本把碳价从0到300以10为步长遍历一遍每次求解后把总成本、弃风率、碳排放量记录到数组里最后画成曲线。这个脚本极其简单但很多做政策分析的人会忽略——他们更喜欢手工改参数再跑浪费大量时间。carbon_prices 0:10:300; res_cost zeros(size(carbon_prices)); res_curtail zeros(size(carbon_prices)); for k 1:length(carbon_prices) carbon_price carbon_prices(k); % 重新构造目标函数并求解 sol optimize(Constraints, objective, options); res_cost(k) value(objective); res_curtail(k) sum(P_wind_pred - P_wind_use) / sum(P_wind_pred); end6.3 与更细时间尺度的联合分析当前模型时间步长是1小时这对储热罐的充放热策略来说已经够用但实际储热罐的充放热切换是有时间常数的——从充热模式切换到放热模式需要阀门调整时间这个切换次数在1小时步长下会被忽略长期运行可能导致设备损耗比模型预期的高。扩展方向是把时间步长细化到15分钟储热罐约束还是那几条只是T从24变成96负荷和风电预测曲线也要细化。这个扩展对模型求解压力不大变量数也就翻了4倍Gurobi处理起来依然轻松。还有一个方向是加入储热罐的寿命损耗模型把充放切换次数作为额外惩罚项写进目标函数让求解器在调度时自动权衡使用频率和寿命损耗。代码整体是模块化的机组参数表、储热罐参数、目标函数构造、约束生成、求解配置、后处理画图这几块都是独立的扩展新场景基本不需要改动已有结构只改参数和约束即可。我个人在实际操作中的体会是这类项目最容易翻车的不是优化算法本身而是数据口径和约束逻辑的一致性——建议拿到题先花半天把数据单位、变量定义、约束符号约定用文档固定下来再动手写代码后面调试能省出几倍的时间。最后再分享一个小技巧求解完后把sol.solveroutput.info里的迭代信息打印出来看一眼如果显示有大量约束被预处理删除说明原始模型里冗余约束太多精简掉能让后续灵敏度分析跑得更快。