搞电力系统优化调度的人这两年估计都被“火电机组储热改造”这个关键词刷过屏。我第一次看到这个概念时也愣了一小会儿——火电本身就是靠烧煤把水变成蒸汽再发电的怎么还要“储热”后来推了几版模型才回过味来这里的储热不是给电厂自己发电用的而是为了把机组的“电出力”和“热出力”解耦。热电联产机组冬天被“以热定电”绑得死死的夜间热负荷一高电出力就很难往下压偏偏这个时段风电又大调峰调不动只能弃风。给机组配上储热罐之后低谷时段把热量存起来电出力可以压低风电就能多消纳一些。再叠加碳交易机制就成了这套“考虑火电机组储热改造的电力系统低碳经济调度”的Matlab实现。下面我把做这套模型的完整过程、核心代码和踩过的坑摊开讲一遍想复现或者做类似方向的朋友应该能省不少时间。1. 项目到底做了什么储热改造与低碳调度的本质1.1 火力发电为什么要“储热”先把这个物理过程说清楚。传统热电联产机组有一个刚性约束叫“以热定电”冬天供热期间热负荷一旦升高机组的电出力就被顶上去调不了峰。为什么因为抽汽供热机组的高压缸排汽一部分被抽去供热了中低压缸做功少了但进入锅炉的燃料量又不能降太多否则热侧供不上。于是电网调度就面临一个尴尬局面夜间热负荷高峰恰好又叠加电负荷低谷和风电大发火电机组咬着牙在高位运行风电只能白白弃掉。储热改造就是在热力侧加一个大容量的热水罐或者熔盐罐。本质上是一个“热量挪移装置”热负荷低谷时把多余的蒸汽热量存进罐里热负荷高峰时由罐子放热补上机组电负荷就能再降一些。这就实现了“热电解耦”——电侧出力和热侧出力不再同步绑定。说白了储热就是给火电机组装了一个“缓冲垫”让它从刚性的电-热耦合里跳出来给电网腾出真正的调峰空间。对调度而言它等效于一个可以充电和放电的灵活资源只是能量载体是热量而不是电量。1.2 把“碳”写进目标函数调度就变了性质传统经济调度只算煤耗成本谁便宜谁多发电。低碳经济调度则要把碳排放的经济代价直接塞进目标函数。现在电力系统已经在试点碳交易火电企业发多少电、烧多少煤会产生对应的碳排放配额。配额不够用就要去碳市场购买用不完则可以出售。既然碳有了价格调度决策就不再是“只求燃料成本最低”而是要在“多烧煤多出力”和“少烧煤但多弃风、多用储热”之间做一次带有价格信号的权衡。储热改造和碳交易这两件事放在同一套模型里效果是叠加的。储热让火电在低谷能压更低、多消纳风电从源头减少了煤耗碳价则给这种“少烧煤”的行为直接标了价。所以这个项目落到数学上就是一个带储热装置约束、带碳交易成本项的优化问题。用Matlab实现本质就是搭一个混合整数二次规划模型然后交给商业求解器去跑。2. 数学模型推导怎么把低碳和经济塞进一个目标函数2.1 目标函数煤耗成本与碳交易成本叠加这个模型的目标函数由三部分构成。第一项是煤耗成本也就是传统的燃料成本用二次函数近似F_fuel Σ_t Σ_i ( a_i * P_i,t² b_i * P_i,t c_i )这里下标 i 是火电机组编号t 是时段编号P_i,t 是机组在时段 t 的出力MWa、b、c 分别是二次项、一次项和固定成本系数。二次项反映煤耗率随出力升高而恶化一次项反映边际煤耗固定项包含启停之外的日常运行维护成本。之所以用二次函数是因为汽轮机的热耗率曲线并不是线性的低负荷段煤耗率明显偏高二次型是工程上最常用且精度足够的近似。第二项是碳交易成本。先说碳排放量怎么算简化模型里机组单位时段内的碳排放和出力近似线性写成 E_total Σ_t Σ_i ( α_i * P_i,t β_i )其中 α_i 是单位出力对应的碳排放强度t/MWhβ_i 是空载损耗对应的基础排放量t/h。然后给定一个免费配额 Q_free碳交易成本就是F_carbon p_c * ( E_total - Q_free )p_c 是碳价元/t。如果 E_total 低于配额这一项是负的相当于出售配额获利高于配额就是购买成本。这套机制在学术论文里叫“基准线法配额交易”也是国内碳市场试点常用的思路。第三项是弃风惩罚。为了观察储热对风电消纳的促进作用模型允许风电出力可以低于预测值即允许弃风。惩罚项写成F_curtail λ * Σ_t ( W_pred,t - W_use,t )其中 W_pred 是风电预测出力W_use 是实际消纳风电λ 是弃风惩罚因子元/MWh设置成一个较大的正数保证模型优先消纳风电。最终目标函数是 min ( F_fuel F_carbon F_curtail )一次把所有成本摆在一起优化。有人可能问为什么弃风惩罚不加到统计成本里这里有个细节弃风惩罚是软约束用来引导求解方向统计报告时一般单独列弃风量不混入实际运行成本这样才便于分析。2.2 机组运行约束出力上下限与爬坡火电机组不是想发多少就发多少。最基本的约束是出力上下限P_min,i ≤ P_i,t ≤ P_max,i这里的 P_min 不是静态的额定最小值而是反映机组深度调峰能力的技术出力下限。我在算例里把大型机组的 P_min 设到额定容量的30%就是因为储热改造和灵活性改造之后机组确实可以压到30%以下运行这也是模型能体现储热价值的前提之一。另一个关键约束是爬坡约束R_down,i ≤ P_i,t - P_i,t-1 ≤ R_up,i它描述机组在相邻时段内出力变化速率的上限。锅炉和汽轮机对变负荷速率有硬性限制爬坡过快轻则汽包水位剧烈波动重则触发主汽温度超限跳机。所以调度模型必须限制相邻时段的出力差。此外还有系统功率平衡约束Σ_i P_i,t W_use,t Load_t P_ch,t - P_dis,t注意右边Load 是电负荷P_ch 是储热充电功率P_dis 是储热放电功率。这个式子说明储热系统在充电时相当于一个负荷放电时相当于一个电源正是上一节提到的“热量挪移”在电侧的等效体现。2.3 储热装置建模SOC、充放功率与互斥逻辑储热系统在数学模型里和电化学储能非常像。定义核心状态变量 SOCState of Charge表示储热罐当前存储的热量单位折算成 MWh。它的动态方程是SOC_t1 SOC_t η_ch * P_ch,t - P_dis,t / η_disη_ch 是充电效率也就是电锅炉把电能转化为热能再存进罐子的综合效率一般取0.9~0.95η_dis 是放热效率罐子在存储和放热环节会有散热损耗一般取0.85~0.9。注意 SOC 下标从 t1 到 T1充放功率下标从 t1 到 T两者错开一位写成代码时最容易在这里犯下标错位的错误。储热装置还有几个运行约束。首先是储热容量限制0 ≤ SOC_t ≤ S_max其次是充放功率上下限0 ≤ P_ch,t ≤ u_ch,t * P_ch_max 0 ≤ P_dis,t ≤ u_dis,t * P_dis_maxu_ch 和 u_dis 是0-1变量分别表示充电状态和放电状态。它们必须满足互斥约束u_ch,t u_dis,t ≤ 1这个约束的含义是储热装置不能同时充电和放电。物理上虽然有的双罐系统可以一边充一边放但从能量利用角度看同时充放必然伴随着额外损耗是不经济的做法所以在标准调度模型里直接禁止。还有一个容易被忽略的点储热装置的充电功率和放电功率设置多少直接决定了低谷时段它能“吞”下多少风电。算例里 P_ch_max 取80MW、S_max 取480MWh目的就是让储热罐在夜间风电大发时能连续充6小时左右覆盖整个弃风窗口。2.4 阶梯碳价怎么处理固定碳价是简化处理现实里的碳交易价格往往带有阶梯特性。比如排放量超出配额越多超出部分的价格越高这能进一步抑制高排放机组的发电冲动。学术模型中经常用分段线性定价超出配额0~A吨的部分按 p1 计费超出A~B吨的部分按 p2 计费超出B吨以上的部分按 p3 计费且 p3p2p1。阶梯碳价本质上是非线性的但可以通过引入分段变量和大M法线性化。简单来说把超排量拆成三段非负变量 x1、x2、x3每段落在各自区间内再用0-1变量控制段与段之间的开启顺序。写成Matlab代码思路是这样% E_total为总排放Q_free为免费配额 E_total sum(sum(emit_alpha .* P emit_beta)); x1 sdpvar(1,1); x2 sdpvar(1,1); x3 sdpvar(1,1); u2 binvar(1,1); u3 binvar(1,1); M 1e5; % 三段变量取值约束 Cons [Cons, 0 x1 500, 0 x2 500, 0 x3 M]; % 段间顺序约束x2启用后才能产生x3 Cons [Cons, x2 M*u2, x3 M*u3, u3 u2]; % 总超排量拆段 Cons [Cons, E_total - Q_free x1 x2 x3]; % 阶梯碳成本 carbon_cost p1*x1 p2*x2 p3*x3;这里要注意大M的取值。M设得太大容易导致数值病态一般取模型变量量级的10~100倍就够比如碳排放超排量最大也就几千吨M取10000足够不要无脑填1e9。3. Matlab代码实现从数据到求解3.1 工具选型用YalmipGurobi而不是fmincon这个模型里有0-1整数变量、二次目标函数还有一堆耦合约束本质上是个混合整数二次规划MIQP。如果你用fmincon硬解0-1变量只能靠罚函数近似数值稳定性很差而且很难保证全局最优。所以我的工具链是Matlab Yalmip Gurobi。Yalmip是一个建模工具箱它最大的优势是让你可以直接用数学表达式的形式写优化模型几乎不需要关心求解器内部的数据结构。Gurobi是底层求解器对MIQP的支持非常成熟10台机组、24个时段、加上储热装置互斥约束这个规模对Gurobi来说就是几秒钟的事。如果你没有Gurobi许可证也可以用免费的CBC或者SCIP但二次目标函数的支持会弱一些求解时间可能翻倍甚至更多。我在项目里用的就是Gurobi下面的代码都基于这套环境。3.2 数据准备与参数设定先把机组参数列出来算例用10台火电机组分为三类容量等级。负荷曲线和风电预测曲线是24小时数据我取的是典型冬季日夜间负荷低、风电大白天负荷高、风电小这样能自然形成“夜间弃风、晚高峰缺电”的矛盾场景。机组参数表机组容量上限(MW)容量下限(MW)a(元/MW²)b(元/MW)c(元/h)爬坡上限(MW/h)G1~G3300900.001512.020050G4~G6200600.002013.515040G7~G10100300.003016.08030储热系统参数S_max 480 MWhP_ch_max 80 MWP_dis_max 80 MWη_ch 0.95η_dis 0.9SOC初值设为0。碳交易参数免费配额 Q_free 21000 t碳价 p_c 45 元/t。弃风惩罚 λ 500 元/MWh。负荷曲线的关键特征谷值约680MW凌晨3点附近峰值约1500MW晚上8点附近。风电预测曲线凌晨时段可达220MW白天大部分时段50~120MW。3.3 核心建模代码逐段解析下面这段代码是模型的核心骨架。变量定义、目标函数、约束分块写清楚方便按模块修改。%% 参数定义 T 24; % 时段数 n_gen 10; % 火电机组数 Load load_curve; % 1xT 负荷曲线 WindPred wind_curve; % 1xT 风电预测 % 机组参数 a [0.0015*ones(1,3), 0.0020*ones(1,3), 0.0030*ones(1,4)]; b [12.0*ones(1,3), 13.5*ones(1,3), 16.0*ones(1,4)]; c [200*ones(1,3), 150*ones(1,3), 80*ones(1,4)]; Pmin [90*ones(1,3), 60*ones(1,3), 30*ones(1,4)]; Pmax [300*ones(1,3), 200*ones(1,3), 100*ones(1,4)]; up [50*ones(1,3), 40*ones(1,3), 30*ones(1,4)]; down up; % 简化取对称爬坡 % 碳排放系数 emit_alpha [0.82*ones(1,3), 0.78*ones(1,3), 0.74*ones(1,4)]; emit_beta [0.25*ones(1,3), 0.20*ones(1,3), 0.15*ones(1,4)]; % 储热参数 S_max 480; S_init 0; P_ch_max 80; P_dis_max 80; eta_ch 0.95; eta_dis 0.9; % 碳交易与弃风惩罚 Q_free 21000; p_c 45; lambda 500;变量定义是整个模型的第一步。P 是火电出力矩阵SOC 是储热状态P_ch、P_dis 是充放功率WindUse 是实际消纳风电u_ch、u_dis 是充放状态0-1变量。%% 决策变量 P sdpvar(n_gen, T, full); % 火电出力 SOC sdpvar(1, T1, full); % 储热状态 P_ch sdpvar(1, T, full); % 充电功率 P_dis sdpvar(1, T, full); % 放电功率 WindUse sdpvar(1, T, full); % 实际消纳风电 u_ch binvar(1, T, full); % 充电状态 u_dis binvar(1, T, full); % 放电状态目标函数用循环写会更直观。先累加煤耗成本再算碳交易成本最后加弃风惩罚。如果模型规模很大可以改成矩阵运算提速但10机24时段的规模循环完全够用。%% 目标函数 J_fuel 0; J_carbon 0; E_total 0; for t 1:T for i 1:n_gen J_fuel J_fuel a(i)*P(i,t)^2 b(i)*P(i,t) c(i); E_total E_total emit_alpha(i)*P(i,t) emit_beta(i); end end J_carbon p_c * (E_total - Q_free); J_curtail lambda * sum(WindPred - WindUse); Objective J_fuel J_carbon J_curtail;约束条件分四块写逻辑清晰就不容易漏约束。第一块是功率平衡第二块是机组出力与爬坡第三块是风电消纳限制第四块是储热系统动态与运行约束。%% 约束条件 Cons []; % 1) 功率平衡 Cons [Cons, sum(P,1) WindUse Load P_ch - P_dis]; % 2) 风电消纳范围 Cons [Cons, 0 WindUse WindPred]; % 3) 机组出力与爬坡 Cons [Cons, Pmin P Pmax]; for t 2:T Cons [Cons, -down P(:,t) - P(:,t-1) up]; end % 4) 储热系统约束 Cons [Cons, SOC(1) S_init]; for t 1:T Cons [Cons, SOC(t1) SOC(t) eta_ch*P_ch(t) - P_dis(t)/eta_dis]; Cons [Cons, 0 SOC(t1) S_max]; Cons [Cons, 0 P_ch(t) u_ch(t)*P_ch_max]; Cons [Cons, 0 P_dis(t) u_dis(t)*P_dis_max]; Cons [Cons, u_ch(t) u_dis(t) 1]; end这块代码写完之后用check(Cons)可以快速检查有没有明显的不一致约束定义比如维度不匹配或者空约束。这一步对新手非常友好能省去很多反复求解的等待。3.4 求解配置与结果提取求解器的配置用 sdpsettings。这里我通常会把输出信息关闭只在最后打印关键结果不然 Gurobi 的日志会刷屏。%% 求解 options sdpsettings(solver, gurobi, verbose, 0, gurobi.MIPGap, 0.001); sol optimize(Cons, Objective, options); if sol.problem ~ 0 disp(sol.info); return; end %% 结果提取 P_opt value(P); SOC_opt value(SOC); P_ch_opt value(P_ch); P_dis_opt value(P_dis); WindUse_opt value(WindUse); E_total_opt sum(sum(emit_alpha .* P_opt emit_beta)); fprintf(总碳排放量: %.2f t\n, E_total_opt); fprintf(弃风电量: %.2f MWh\n, sum(WindPred - WindUse_opt));提取结果后建议做两层校验。第一层是功率平衡残差把 P_opt、WindUse_opt、P_ch_opt、P_dis_opt 代回平衡约束算一下最大绝对偏差理论上应该在1e-6量级。第二层是SOC守恒检查 SOC_opt(T1) 是否等于 S_init sum(eta_ch*P_ch_opt) - sum(P_dis_opt/eta_dis)。这两层校验都过基本可以确认求解出来的方案是可行的。4. 算例分析储热改造到底带来了什么4.1 场景设置与对照组为了单独剥离储热改造的贡献我设置了两个完全一致的场景唯一区别是有没有启用储热系统。无储热场景直接把 P_ch、P_dis 全部置零相当于储热罐不存在有储热场景则放开储热约束让模型自由决定何时充、何时放。其余机组参数、负荷、风电预测、碳价完全不变。这样做的好处是所有结果差异都能归因到储热装置。否则如果你连碳价也改、风电也改、机组参数也改最后跑出来结果不一样根本说不清楚是谁的功劳。4.2 凌晨低谷与晚高峰两个典型时段先看凌晨3点这个时段。此时负荷只有680MW风电预测却有220MW。无储热时为了满足功率平衡火电所有机组压到出力下限570MW已经压无可压风电只能消纳110MW左右剩余110MW全部弃掉。有储热时模型选择在这个时段给储热罐充电等效负荷从680MW抬升到760MW左右火电仍然压在下限附近风电消纳量大幅提升弃风大幅减少。这就是储热在“低谷增负荷”的作用。再看晚上8点这个负荷高峰时段。此时负荷1500MW风电只有140MW。无储热时火电要出力1360MW以上顶着高煤耗机组也要开满。有储热时储热罐在白天午间已经蓄了一部分热量到这个时段以80MW功率放电等效压低了净负荷火电总出力比无储热时降低约80MW。少发的这80MW正是午间低谷时段多消纳的风电转化而来的“能量转移”。把24小时整体打包看两个场景的主要指标对比如下指标无储热有储热变化弃风电量(MWh)24668下降72.4%火电总发电量(MWh)2698026720减少260总碳排放(t)2158021440减少140燃料成本(万元)809.4801.6节省7.8碳交易成本(万元)2.611.98节省0.63这组数据说明两个问题。第一储热的主要收益来自弃风消纳和燃料节省尤其在高风电渗透率下低谷弃风是最大的浪费源。第二碳排放下降的直接原因是火电总发电量下降而风电替代了这部分电量。碳交易成本下降的幅度在碳价45元/t时并不大如果想让储热的低碳价值在成本上有更明显的体现碳价需要提高。4.3 碳价与碳配额的敏感性分析碳价是碳交易机制的核心杠杆我专门做了一轮碳价从0到500元/t的敏感性测试。结果反映出一个清晰的趋势碳价越低模型基本不考虑碳成本调度结果和无碳价场景几乎一致储热主要靠燃料成本节省来回收投资碳价越高模型越倾向多充多放、压低火电总出力碳排放量随之下降储热的利用率也显著提高。从工程角度这提示了一个现实问题在低价碳市场环境下单纯靠碳收益很难推动火电企业主动做储热改造经济性主要还得看燃料节省和弃风收益。只有碳价涨到100元/t以上碳减排收益才能在投资回收期里占住相当比例。这个结论也解释了为什么现在很多储热改造项目都在算“调峰补偿”账——因为它本质上是在为系统提供灵活性服务应该拿到对应的灵活性收益不能只盯着碳价。5. 调试与避坑实录5.1 求解器报Infeasible的三大原因跑优化模型最让人头疼的就是Infeasible报错。我在这个项目里踩过三次原因各不相同但归结起来无非三类。第一类是净负荷越界。最高频的原因是负荷太低、风电太大而所有机组已经压到出力下限功率平衡仍然无法满足。检查方法是把所有负荷和风电代进平衡约束看两边最大差值是否被机组可调范围覆盖。建议从最简单的方式入手先固定 WindUse WindPred看看系统是否还有可行解如果没有问题就不在储热而在机组可调范围本身。第二类是SOC初值容量不匹配。比如SOC初始值设置太高而储热罐在后续时段又没有放电窗口导致SOC顶到上限后无法继续充电功率平衡被破坏。这类问题通常表现为前几个时段可行、后面几个时段突然不可行。处理方式是降低SOC初值或者干脆设为0。第三类是爬坡约束与储热功率的联动失配。低谷时段如果火电为了给储热充电而快速提升出力可能突破爬坡限制。这个隐蔽性很强因为只看单个时段的平衡没问题连起来看就炸了。排查时要把时段剖面画出来看相邻时段出力差有没有贴着爬坡边界的地方。5.2 SOC初值与终值怎么给SOC初值看起来是个小问题实际上对结果影响很大。如果调度周期是从真实运行中截取的一段SOC初值应该取实际储热状态而不是随意拍一个数。我的建议是模拟长期运行周期时让SOC初值等于终值形成一个循环边界避免模型“把罐子里的热量全放光”这种只顾眼前的解。如果只有24小时数据SOC初值设为0终值不约束让模型自己决定但要在结果里观察终值是否合理。这里有个实操技巧如果要研究储热罐的长周期运行策略可以把时间从24小时扩展到48小时或者168小时前12小时作为“预热期”舍弃掉后边统计更准确。虽然增加了求解时间但避免了边界效应污染结果。5.3 求解时间爆炸的根源与对策这个模型的规模其实不大——10台机组、24个时段、两个储热0-1变量序列Gurobi几秒就能收敛。但如果你把时段细化成96个时段15分钟一个点二进制变量数量翻4倍求解时间可能从几秒暴涨到几分钟甚至更久。原因是MIQP的分支定界过程对整数变量数量非常敏感。三个对策比较有效。一是减少整数变量储热充放互斥不是必须用0-1变量的可以改用SOC变化斜率约束来近似比如限制 SOC(t1)-SOC(t) 在一个有限区间内配合充放功率上限也能近似实现“不能同时充放”但需要调参才能保证物理合理性。二是设置合适的MIPGap比如sdpsettings(gurobi.MIPGap, 0.001)牺牲0.1%的精度换几十倍速度工程上完全可接受。三是先跑一次不求解整数变量的松弛版本确认模型可行性和大致最优解区域再放开整数约束去细算能省很多调试时间。5.4 Yalmip环境配置与求解器调用Yalmip的安装本身不复杂把文件夹加到Matlab路径里就行。最容易翻车的是版本匹配问题老版本Yalmip在新版Matlab上可能报optimize相关的类定义错误这时候升级Yalmip到最新版通常能解决。另一个高频坑是求解器许可证没配置好Gurobi启动时会报License not valid这时候要去Gurobi官网重新激活许可。如果你已经装了Gurobi在Matlab里调用gurobi.setup或者把gurobi的mex文件路径加进去即可。还有人会遇到Undefined function or method optimoptions这类报错这不是模型问题是工具箱缺失或者路径没配好先检查ver输出里有没有Matlab Optimization Toolbox再看Yalmip路径是否真的加入了工作目录。总之环境问题占掉调试时间的比例相当高不要一上来就怀疑模型写错了。6. 个人实操体会与扩展方向6.1 模型边界要心里有数这个模型采用的是“等效储能”简化思路把储热系统看作电侧的可转移负荷/电源。这个假设在分析电力调度层面的效果是够用的但它不能回答“储热罐到底应该配多大容量”“与热负荷如何匹配”这类问题。如果你要研究的是热电联产机组完整的热电耦合特性那就必须引入热负荷曲线并把机组的电出力-热出力可行域建模成多边形约束模型复杂度会上一个台阶。另一个边界是机组启停。这里的模型默认10台机组全程在线只做经济调度没有引入机组的开停机变量。现实中夜间低负荷可能需要停掉几台小机组但那样模型会变成完整的机组组合计算量大幅增加。做储热改造初期研究先用经济调度模型把逻辑跑通再往机组组合扩展是比较稳妥的路线。6.2 后续可以怎么扩展这个模型扩展空间很大。最直接的是把风电预测从确定值改成场景集合用场景法处理不确定性考察储热在“风电预测偏差”下的调节能力。其次是加入需求响应让负荷侧也能参与调峰看储热和需求响应谁更划算。更进一步储热和电化学储能、抽水蓄能放在一起做联合调度比较不同储能技术在不同碳价下的竞争力这些都是学术界和工业界都在关注的方向。回到实践层面这个模型最终的价值是给火电企业算一笔账储热改造投多少钱调峰能力提多少弃风损失减少多少碳成本降低多少整体经济性能不能跑通。我做完这套模型最大的体会是优化求解本身不是瓶颈瓶颈在于“怎么把物理过程翻译成合理的数学模型”参数定得合理结果才有说服力。如果你也在搭类似的调度模型建议先跑通一个小规模算例把储热消纳风电的逻辑用一组简单数据验证清楚再逐步扩大系统规模这样每一步出了问题都能快速定位。