去年冬天我在做一个区域电网日前经济调度项目最头疼的不是电负荷而是供暖期火电机组“以热定电”导致的风电弃电。白天勉强能平衡一到夜间风电大发、热负荷也跟着踩上高峰机组为了保住供热必须顶着电出力下限运行风电场几乎每隔两天就被调度压一次出力。当时和我合作的工程师提了一个思路在热网侧加蓄热罐让机组的热出力不必被热负荷完全绑死。我把这个想法写进了低碳经济调度模型用Matlab跑通了一整套代码也踩了不少坑。这篇文章就把这套“火电机组储热改造低碳经济调度”的建模思路、Matlab代码框架和调试经验完整整理出来给正在做灵活性改造、新能源消纳和低碳调度方向研究或工程项目的朋友一个可复现的参考。1. “以热定电”的困局为什么储热改造会成为调度的破局点1.1 供热期机组调峰受限的根源抽汽供热机组在中压缸某个级后抽取一部分蒸汽去加热热网水这部分抽汽不再继续参与发电膨胀做功所以机组的电出力范围和供热抽汽量是强耦合的。供热抽汽量越大能够发电的蒸汽通流比例越低机组电出力可调范围也越窄。更麻烦的是为了保证供热参数机组往往存在一个由热负荷决定的电出力下限这个下限在严寒期会被抬得很高。一句话概括电负荷低谷和热负荷高峰同时出现在夜间而机组为了“保热”不得不维持一个较高的电出力电网调度只能让火电硬给风电“让路”。这就是“以热定电”限制调峰能力的直接体现。生活一点理解相当于你家暖气片必须一直烧着但夜里没人开灯发电机却不能因为你不用电就停下来——它在夜里不但要发电还要按最小的那部分发电量去“保底”。从电网整体看这个保底出力挤占了新能源的上网空间。1.2 热网侧蓄热罐把热负荷从低谷搬到高峰储热改造的方式有很多比如机组旁布置大型热水蓄热罐、电极锅炉配合蓄热、固体蓄热等。我这次用的建模方式是“热网侧集中式蓄热罐”也就是在热力站和机组之间加一个带斜温层的大型热水罐利用供回水温差储存热量。夜间低谷时段蓄热罐从热网侧抽取热量存储起来相当于机组面对的热负荷暂时变小白天电负荷上升后蓄热罐再放热把这部分热量补回热网。这个过程本质上是把“热负荷的时序曲线”做了一次平移——将低谷时段的热需求搬到高峰时段。机组不再需要严格跟随热负荷实时出力而是可以跟随“净热负荷”调度。由于机组热出力下限和电出力下限绑在一起热负荷曲线被削平之后电出力下限也会跟着降低风电在低谷时段的消纳空间就出来了。这就是所谓“热电解耦”的通俗解释用储热这条柔性路径在电和热之间加一个缓冲打破“热电必须同时平衡”的强约束。1.3 从物理改造到调度模型到底要新增什么变量物理上加了蓄热罐调度模型里就要多出三类东西储热罐的充热功率和放热功率变量储热罐的储热状态变量代表当前罐里存了多少热量充、放热工况对应的0-1标志变量用来保证同一时段不能边充边放。同时题目里“低碳”二字意味着碳排放相关变量也要进模型比如实际碳排放量、配额量、碳交易量。原本的“纯电调度模型”升级成“电-热-碳联合调度模型”机组出力和热出力同时被决策优化结果要在一个目标函数里平衡燃料成本、碳排放成本和弃电风险。理解了新增变量的意义后面代码才不会写乱。2. 储热改造机组的数学刻画可行域、储热罐与碳成本建模2.1 热电联产机组的电热可行域抽汽式热电机组在电-热二维平面上存在一个凸多边形可行域工程上常用顶点法描述。比如一台300MW抽汽机组电出力范围大约在150300MW供热出力范围在0250MWth但这两个约束不是独立的顶点表给出的是一组线性不等式供热出力越大电出力下限和上限都会相应移动机组存在最大进汽量约束对应电出力上限随热出力增加而下降最小凝汽流量约束对应电出力下限随热出力增加而上移。实际建模时把顶点坐标转换成行列式或线性不等式即可。本文算例里我用的是简化参数两条机组曲线如下表所示方向性比较清晰机组类型电出力上限 MW电出力下限 MW热出力上限 MWth排放强度 kg/MWh电排放强度 kg/GJ热G1300MW抽汽机组300120250820110G2200MW抽汽机组20080180820110G3600MW纯凝机组60024008600这里要先说清楚我采用的是热网侧集中式储热机组本身的电热可行域不变变化的是机组需要满足的净热负荷。两种建模思路区别很大——机组侧储热会直接改变机组可行域形状而热网侧储热只改变热平衡方程右侧的负荷曲线。工程中热网侧蓄热罐实施更灵活调度模型改起来也更可控所以我优先推荐这种方式。2.2 集中式储热系统的动态约束方程蓄热罐的数学模型核心是能量递推关系。假设调度步长 dt1h第 t 个时段初的储热状态为 S_t则S_{t1} S_t η_ch · C_ch,t - C_dis,t / η_dis其中 C_ch,t 是充热功率C_dis,t 是放热功率η_ch 和 η_dis 分别是充放热效率本文取 0.95。这个递推方程说明储热罐是有“记忆”的——当前状态影响后续所有时段的可用热量不能当成简单的静态设备。其余必要约束包括储热容量上下限0 ≤ S_t ≤ C_max比如本文算例取 C_max1000 MWh充放功率上下限0 ≤ C_ch,t ≤ C_ch_max·u_ch,t0 ≤ C_dis,t ≤ C_dis_max·u_dis,t同时性约束u_ch,t u_dis,t ≤ 1防止边充边放这种无意义的能量空转周期约束S_{T1}S_0保证一个调度周期结束时储热状态回到初始值。热平衡方程也要改。原来的形式是“机组热出力总和热负荷”加入储热后变成Σ_j H_j,t C_dis,t - C_ch,t H_load,t也就是说机组总供热量可以大于或小于当前热负荷差额由蓄热罐充放热来补。注意符号方向放热是给系统加热量充热是从系统取热量这一步一旦写错整个结果都是反的。2.3 碳排放成本模型与目标函数里的“低碳”分量碳排放成本最常见的写法是“配额-实排”差额交易模式。系统每个调度时段获得一定免费碳排放配额 E_alloc同时按实际发电量和供热量计算实际碳排放 E_actualE_actual Σ_i (e_P,i · P_i,t e_H,i · H_i,t)碳交易成本为C_carbon,t γ · (E_actual - E_alloc)γ 是碳价单位为元/吨。配额可以按机组容量分配也可以按历史排放水平设定。研究性仿真里常将配额设为 0相当于对全部碳排放征收碳成本工程测算则按既有配额规则取值。我这里做敏感性分析时会把 γ 从 0 扫到 400 元/吨观察它对弃风率和碳排放量的影响。完整的目标函数写成min Σ_t [ Σ_i F_i(P_i,t, H_i,t) Σ_i SU_i,t γ·(E_actual - E_alloc) λ·(P_w_fore,t - P_w,t) ]其中 F_i 是煤耗成本SU_i,t 是启动成本最后一项是弃风惩罚成本λ 取 300 元/MWh代表风电被压出力时给予一个机会成本惩罚。这样“低碳”就不仅仅是口号而是通过碳价进入目标函数与经济性在同一框架内做权衡。2.4 模型汇总约束的类型与含义把全套模型列清楚方便对照代码逐条检查约束类型表达式含义电功率平衡ΣP_i,t P_w,t P_load,t发电量等于电负荷热功率平衡ΣH_j,t C_dis,t - C_ch,t H_load,t机组供热和储热共同满足热负荷机组出力约束P_i,min ≤ P_i,t ≤ P_i,max电出力上下限热出力约束H_j,minleq H_j,t ≤ H_j,max热出力上下限电热耦合约束顶点线性不等式描述抽汽供热可行域爬坡约束-R_down ≤ P_i,t - P_i,t-1 ≤ R_up机组爬坡速率限制风电出力约束0 ≤ P_w,t ≤ P_w_fore,t风电预测值作为上限旋转备用约束ΣP_i,max·u_on,i,t ≥ P_load,t R_t系统留有备用容量储热递推约束S_{t1}S_tη_ch·C_ch,t-C_dis,t/η_dis储热状态时序耦合储热容量约束0 ≤ S_t ≤ C_max储热罐容量限制充放同时性约束u_ch,tu_dis,t≤1防止边充边放碳交易约束C_carbon,tγ(E_actual-E_alloc)碳成本入目标函数这是一个典型的混合整数线性规划MILP问题。连续变量是各机组的电出力、热出力、储热功率、储热状态、风电出力和碳排放量0-1变量是机组启停和储热充放标志。规模不大用常规商用求解器可以秒解。3. Matlab实现路线Yalmip建模与求解器选择3.1 为什么我优先推荐Yalmip而不是手写矩阵用MATLAB做这类调度优化两条路摆在面前一是用优化工具箱的linprog/intlinprog把模型展开成标准矩阵形式二是用Yalmip建模工具箱直接以变量和约束语句描述数学模型。我强烈推荐后者原因是电热耦合、跨时段递推这类约束手写矩阵既容易算错索引又极难排查。Yalmip的写法接近数学表达式变量声明、约束拼接、目标函数书写都很自然调试时还能查看diagnostics信息。更重要的是Yalmip底层可以无缝切换CPLEX、Gurobi、intlinprog等求解器前期开发用一个免费求解器就能跑通正式算例再切换到商用求解器。3.2 代码框架与关键片段整体代码目录分四块参数输入区、变量定义区、约束生成区、求解与结果导出区。核心片段如下我已做过简化便于复现%% 参数区 T 24; % 调度时段数h nG 3; dt 1; % 步长1h % 机组参数Pmax, Pmin, Hmax, Hmin, ramp, 排放强度... G struct(Pmax,[300;200;600], Pmin,[120;80;240], ... Hmax,[250;180;0], Hmin,[0;0;0], ... c1,[230;230;210], c2,[0.15;0.15;0.18], ... eP,[820;820;860], eH,[110;110;0]); % 负荷与预测 P_load load_data(electric_load.xlsx); % 1xT H_load load_data(heat_load.xlsx); % 1xT P_w_fore load_data(wind_forecast.xlsx); % 1xT % 储热参数 HST.Cmax 1000; % MWh HST.cmax 200; % MWth HST.eta_ch 0.95; HST.eta_dis 0.95; % 碳价与配额 gamma 200; % 元/t E_alloc zeros(1,T);% 配额设为0研究性假设 lambda_w 300; % 弃风惩罚元/MWh %% 变量区 P sdpvar(nG, T, full); % 电出力 H sdpvar(nG, T, full); % 热出力 u_on binvar(nG, T, full); % 启停标志 P_w sdpvar(1, T, full); % 风电并网功率 C_ch sdpvar(1, T, full); % 充热功率 C_dis sdpvar(1, T, full); % 放热功率 u_ch binvar(1, T, full); % 充热标志 u_dis binvar(1, T, full); % 放热标志 S sdpvar(1, T1, full); % 储热状态含初末值 %% 约束区 Cons []; % 电功率平衡 Cons [Cons, sum(P,1) P_w P_load]; % 热功率平衡 Cons [Cons, sum(H,1) C_dis - C_ch H_load]; % 机组电出力上下限与热出力上下限 for i 1:nG Cons [Cons, G.Pmin(i).*u_on(i,:) P(i,:) G.Pmax(i).*u_on(i,:)]; Cons [Cons, G.Hmin(i).*u_on(i,:) H(i,:) G.Hmax(i).*u_on(i,:)]; end % 储热递推方程 Cons [Cons, S(1) 500]; % 初始储热态 for t 1:T Cons [Cons, S(t1) S(t) HST.eta_ch*C_ch(t) - C_dis(t)/HST.eta_dis]; end Cons [Cons, S(T1) S(1)]; % 周期约束 Cons [Cons, 0 S HST.Cmax]; % 容量约束 Cons [Cons, 0 C_ch HST.cmax*u_ch]; Cons [Cons, 0 C_dis HST.cmax*u_dis]; Cons [Cons, u_ch u_dis 1]; % 同时性约束 % 风电出力 Cons [Cons, 0 P_w P_w_fore]; %% 目标函数 % 简化线性煤耗成本正式模型中把二次项做分段线性化 fuel_cost sum(sum( G.c1 .* P G.c2 .* H )); carbon_cost gamma * sum(sum( G.eP.*P G.eH.*H ) - E_alloc); wind_penalty lambda_w * sum(P_w_fore - P_w); Objective fuel_cost carbon_cost wind_penalty; %% 求解 ops sdpsettings(solver,gurobi,verbose,2,gurobi.MIPGap,0.001); sol optimize(Cons, Objective, ops);因为电热耦合可行域在简化参数里暂时没有展开实际项目里你需要在“机组电出力上下限”和“热出力上下限”之外加入顶点法生成的线性不等式。我当时处理是用一个小函数把机组顶点表转换成行向量并批量添加到约束集里逻辑清晰且不易出错。3.3 求解器配置CPLEX、Gurobi与intlinprog的取舍如果你有Gurobi或CPLEX许可证直接在sdpsettings里指定solver即可。几个参数值得注意verbose 建议设成 2能看到求解过程定位卡在哪一拍MIPGap 设 0.1% 就够工程上没必要追求 0 对偶间隙TimeLimit 设为 600 秒防止复杂算例卡死。没有商业求解器时Yalmip会自动调用MATLAB自带的intlinprog。模型不大时完全能跑但求解速度比Gurobi慢不少。我实测过同样一个24时段的算例intlinprog要十几秒Gurobi基本零秒解决。所以正式研究还是想办法搞到Gurobi许可证学生和科研机构一般都能申请。3.4 数据组织和参数化的习惯一个容易忽略的点代码千万别写成一坨。我是把所有机组参数、储热参数、负荷曲线都封装成结构体或表格再把主程序写成函数function res lowcarbon_dispatch(gens, hst, load_data, gamma)这样的好处是参数扫描极其方便——比如要做碳价敏感性一个for循环调用这个函数就行。调度周期 T 和步长 dt 也建议做参数从典型日扩展到一周甚至全年时不需要改任何方程逻辑。4. 算例验证与结果解读改造前后调度结果的对比4.1 算例系统与对照组设计我用前面那张表里的3台机组搭了一个简化单节点系统电负荷峰值1200MW热负荷峰值约400MWth风电装机500MW且夜间大风时段预测出力较高。设了三组场景S0不改造机组原生“以热定电”约束生效S1加装储热罐1000MWh容量200MWth充放功率S2加装储热罐同时碳价取200元/吨。仿真结果汇总如下场景弃风率 %总运行成本 万元系统碳排放 t夜间最低电出力 MWS016.8182.45420780S13.9171.65015640S22.8176.24728632先看 S0 和 S1 的对比储热罐带来的变化立竿见影弃风率从16.8%降到3.9%总成本下降约10万元碳排放减少约8%。这说明储热改造在系统层面带来的收益是双重的——既省钱又减碳。4.2 低谷时段的风电消纳变化我把夜间几个典型时段的出力数据拉出来看机理非常清楚。改造前晚上22:00到次日4:00机组为了维持热负荷必须保持较高电出力风电即使预测出力很高也只能上网一小部分剩下的全部被弃。改造后储热罐在2:00-5:00之间持续充热相当于把热负荷从晚上搬到了白天机组电出力下限整体降了约140MW这部分空间完全让给了风电。一个更直观的说法储热罐在低谷时段“吃掉”了火电的多余热量让火电在电侧能够更大幅度地压出力风电才真正有机会多发。如果没有这个缓冲单纯靠优化算法“调整权重”火电的物理约束摆在那里再怎么调也挤不出多少空间。4.3 碳价变化对调度决策的影响碳价是“低碳经济调度”里最容易做敏感性分析的参数。我分别设0、100、200、400元/吨扫描结果见表碳价 元/t弃风率 %总运行成本 万元碳排放 t储热罐平均充热量 MWh03.9171.650156101003.4173.948806502002.8176.247286724002.1183.54460720碳价从0涨到400弃风率从3.9%降到2.1%碳排放减少约11%。原因是碳价越高系统就越倾向于用风电替代火电而储热罐的存在让这种替代有地方“使力”——低谷时段机组压得下出力风电才能多发。注意没有储热罐的S0场景在碳价上升时碳排放下降非常有限因为火电可行域卡住了这就是“热电解耦能力”决定低碳转型空间的核心逻辑。5. 代码调试中的坑与应对不可行解、非线性项与求解效率5.1 不可行解排查顺序这类模型最常遇到的就是find 后报“infeasible”尤其是第一次加储热约束时。我的排查顺序是固定的先把储热相关约束全部注释掉跑基准模型确认原机组模型本身没有缺陷加入储热递推和容量约束单独检查周期约束看S(1)初值和S(T1)S(1)是否冲突检查热平衡方程符号最容易错的是把C_ch和C_dis的正负号搞反导致蓄热罐一边充一边放才算出平衡用Yalmip的诊断输出区分是primal infeasible还是dual infeasible前者说明约束自相矛盾后者往往是目标函数存在无界项。还有一个常见问题0-1变量和连续变量相乘造成的非凸约束。Yalmip里如果直接写P.*u_on会生成双线性项这是无法用MILP求解器处理的。正确写法是像前面代码那样用不等号分离变量或者引入Big-M辅助变量。这一步是我刚开始频繁踩坑的地方。5.2 非线性项的线性化处理煤耗成本在实际中往往是二次函数储热调度模型即便只保留最关键的乘积项也需要做线性化。我通常的做法是将 P_i,t 的可行范围分割成35段每段用线性函数近似引入分段标志的0-1变量用Yalmip的implies或Big-M约束实现碳价如果采用阶梯碳价超过配额部分分段涨价同样用分段线性化。Big-M的值要谨慎取值。M太大会导致数值病态太小又会错误地削减可行域。通用经验是取该变量物理上限的1.1倍比如电出力上限300MWM取330即可。还有一点目标函数里的弃风惩罚项 λ·(P_w_fore - P_w) 本身就是线性的不要画蛇添足加max或绝对值。风电预测值作为上限约束后压缩风电只会因为经济性和碳排放考虑被自动优化不需要额外引入惩罚逻辑。5.3 求解效率与长时间尺度扩展24小时或168小时的MILP对Gurobi来说非常轻松但如果你要把模型扩展到全年8760小时直接一次性求解会非常慢且内存占用夸张。我的经验是分周滚动求解以一周为窗口前一周的最后储热状态作为下一周的初始值并保留跨周周期约束的松紧程度。这样做虽然会损失一点全局最优性但工程上是完全可接受的。另外约束生成时尽量用矩阵运算批量添加而不是逐时段写for循环。Yalmip虽然会自动做转换但你自己组织的约束集如果本身是稀疏拼接转换效率会高很多。以 T8760 为例储热递推约束如果用循环一列一列加光建模时间就能等上几分钟矩阵化后秒级完成。5.4 让我反复踩坑的细节调试过程中有几个细节说出来供你避雷功率和能量单位别混。步长1h时MW和MWh数值相等但换了15min步长就完全不同储热递推方程里不乘dt结果必然不是你想的那样初始储热状态不要拍脑袋设。我用S(1)500的时候前几个时段充放都在围绕这个初值做调整初值不合理会让储热设备“被迫”先调整状态导致前几个小时的结果明显偏离最优风电预测序列不要全填同一个数。我测试时图省事填了常数结果优化结果完美得不像话浪费了半天检查才发现是数据太理想碳配额设成0是一种研究假设真实系统里配额和基准线直接决定碳交易成本的大小做分析时一定要注明否则结果很容易被误读。如果你正准备复现这套模型我的建议是先拿一个最简算例跑通流程确认储热罐的充放热功率和状态量落在合理区间再扩展完整案例。调度层面的结论也别直接外推到改造项目的投资决策——储热罐单位造价、设备寿命、改造停机时间都没有进入这个模型真要评估项目收益还得叠加投资回收期的测算。对我来说这套代码最有价值的地方是把“热-电-碳”三个时间尺度不同的物理量放进了同一个优化框架里看着储热罐自动选择凌晨充热、白天放热会比读一百遍教科书都直观。