把“考虑火电机组储热改造的电力系统低碳经济调度”用 Matlab 完整跑通一遍之后我发现这个方向的价值不在于多一个算例而在于它把火电灵活性改造、碳市场和优化调度三个问题拧在一起值得动手复现。这篇文章我打算直接从问题源头讲起把数学模型怎么建、储热罐怎么处理、碳成本怎么加、Matlab 代码怎么写、求解器怎么调一条线捋清楚。适合正在做电力系统方向课程设计、毕业设计或者刚入门优化调度的同学参考有一点 Matlab 基础就能跟上。1. 火电机组储热改造到底改了些什么1.1 传统热电联产的“以热定电”约束北方地区冬季供暖期大量热电联产机组承担着供热任务。这类机组内部的热电关系非常密切供热抽汽量一旦确定电功率的下限就基本被“锁”住了。通俗讲热负荷越高机组为了维持供暖而必须保持的蒸发量就越大电出力哪怕不需要那么高也无法继续往下压这就叫“以热定电”。这个约束在冬天会出现一个很尴尬的局面夜间风大、负荷低风电出力正高但热电联产机组因为要供热电出力压不下来结果电网只能优先保障供热把风电机组出力砍掉产生弃风。从调度角度看不是风电不行而是系统里缺少足够的向下调节空间。储热改造的出现本质上就是为了给这类机组“松绑”。通过在火电厂侧增加储热罐或者加装电极锅炉、蓄热式电锅炉把热负荷的一部分搬到储热装置里。供热高峰期机组可以减少直供热出力由储热罐放出热量补充热负荷低谷期机组可以多蓄热把多余热量存起来。这样一来机组的电出力下界不再被热负荷死死拽住电网的调峰空间变大风电就能多消纳一些。所以储热改造被很多人称为“热电解耦”背后逻辑其实很简单电厂烧锅炉产生蒸汽蒸汽既发电又供热储热罐相当于一个热量缓冲池让电和热不再必须同步变化。1.2 储热装置带来的“热电解耦”要理解储热改造后的调度模型先要分清储热和电池储能的区别。电池储的是电储热罐储的是热量两者在能量形式、效率、成本、响应速度上完全不同。在一个简化的系统中储热罐通常有三个状态蓄热、放热、维持。蓄热时热源从机组抽汽或者电锅炉获取热量把热量“灌”进储热罐放热时储热罐把热量输送给热用户替代机组的直供热。调度上需要跟踪储热罐的储热容量也就是剩余热量多少类似于电池的 SOC。改造后热电联产机组的运行区间会发生改变。以一台典型 300 MW 热电联产机组为例改造前供暖期电出力下限可能被抬到 120 MW 以上加装一定容量的储热罐后电出力下限可以被拉低到 60 MW 甚至更低多出来的调节空间就是系统消纳新能源的关键。在我搭建的算例里我采用了一台加装储热罐的热电联产机组配合两台纯凝火电机组和一座风电场。为了对比效果我给储热前后各跑了一遍优化结果非常直观没有储热时夜间弃风明显加入储热后风电出力基本被吃干净火电机组的向下调节压力也小了很多。1.3 低碳经济调度的新变量以前的传统经济调度只盯煤耗和运行成本目标函数就是“怎么花钱最少”。现在把“低碳”放进来主要是多了碳排放成本这一项。碳市场机制下发电企业会有一定量的免费碳排放配额如果实际排放量超过配额超额部分要花钱买碳配额如果排放量低于配额剩余配额可以在市场出售获利。这个机制等于给碳排放贴上了价格标签调度时就会主动权衡多用一点风电少发一些火电碳排放下降可能节省的碳成本比多付出的调峰成本还要划算。所以在调度模型中目标函数里会出现“碳交易成本”这一项而且这一项的正负号必须处理对否则结果会乱套。我个人的体会是火电储热改造和低碳经济调度放在一起特别有意思。储热改造解决的是“新能源能不能消纳”的物理问题碳成本解决的是“减排值不值得”的经济问题。两个问题耦合在一个优化模型里既能看到系统灵活性投资的收益也能看到碳价对运行方式的引导作用。2. 调度模型的数学化与关键取舍2.1 目标函数成本与碳的双层耦合调度模型的核心是目标函数我采用的是典型的小时级调度时间尺度为一天 24 小时单位时段为 1 小时。总目标是最小化系统总运行成本包括火电机组燃料成本、启停成本、碳交易成本以及弃风惩罚成本。燃料成本用二次函数近似这是电力系统调度里最常见的处理方式。第 i 台机组在 t 时段的燃料成本可写为F_i(t) a_i * P_i(t)^2 b_i * P_i(t) c_i其中 P_i(t) 是机组出力a_i、b_i、c_i 是煤耗特性系数。二次项的存在让目标函数成为凸二次规划后面求解时只需要把二次项矩阵传给求解器就行。启停成本是另一个重要部分机组从停机到启动会有额外热耗设备磨损也不小因此模型中引入二进制变量 u_i(t)表示机组是否运行再引入启动成本变量。为了简化很多研究采用固定启动成本我这里的算例也按固定值处理不区分热启动和冷启动。碳交易成本项我建模为碳价乘以“实际碳排放量减去免费配额量”。实际碳排放量由各机组出力乘以碳排放强度再累加得到。当实际排放量超过配额时这一项为正系统要花钱买碳低于配额时这一项为负相当于获得收益。为了安全起见调度模型中碳价通常取一个固定的价格常数比如 120 元/吨这也是目前很多论文中常用的设置。弃风惩罚成本表示弃掉风电造成的浪费用弃风功率乘以一个较大的惩罚系数。这个系数实际是一个“软约束”的替代方案它不会强行要求风电全部消纳而是通过成本激励让模型在系统无法完全吸收风电时优先弃掉最少的风电。总目标函数可以写成min 总成本 燃煤成本之和 启停成本之和 弃风惩罚成本 碳交易成本这样设计有一个明显的好处目标函数的所有项都是线性或凸二次项加上后面约束里会出现整数变量最终模型是一个混合整数凸二次规划MIQP可以用商业求解器高效求解。2.2 系统约束从电功率平衡到储热状态调度模型不能只有目标函数约束才是让结果“符合物理规律”的关键。我按约束类型分四层来看。第一层是电功率平衡约束。每个时段所有火电机组出力、风电实际出力、储热装置如果有电锅炉还需要考虑电转热的电耗但我的算例里储热罐不从电网取电只从机组抽汽取热因此不增加电负荷之和要等于系统电负荷。风电实际出力等于预测出力减去弃风功率。第二层是火电机组运行约束。每台机组都有出力上下限机组启动后出力不能超过最大出力不能低于最小出力。爬坡约束也必须有相邻时段出力变化量不能超过爬坡速率否则调度结果在实际运行中根本执行不了。如果采用二进制机组启停变量还需要把出力与启停状态关联起来常用形式是u_i(t) * Pmin_i ≤ P_i(t) ≤ u_i(t) * Pmax_i这组约束保证机组停机时出力为 0启动后出力落在可行区间内。第三层是热力系统约束。对于热电联产机组需要考虑热负荷平衡。T 时段内机组直供热功率加上储热罐放热功率减去蓄热功率要等于热负荷需要量。这里我把热负荷设为一个已知的 24 小时曲线模拟冬季典型日用热需求。第四层是储热罐约束。储热罐的储热量状态变量 S(t) 表示 t 时段末罐内剩余热量。状态转移方程为S(t1) S(t) * (1 - 散热损失率) η_chg * Q_chg(t) - Q_dis(t)其中 Q_chg(t) 是蓄热功率Q_dis(t) 是放热功率η_chg 是蓄热效率放热过程一般也有效率损失。储热罐容量有上限蓄放热功率也有上限另外蓄热和放热不能同时进行这是常见的工程约束。还有一个容易被忽略的边界条件调度周期开始和结束时的储热量往往要约束成相等或至少保证结束时的储热量不小于某个初始值否则模型只会把储热罐里的热量在最后时段全部放光这在工程上不公平也会让结果失真。2.3 机组模型与储热耦合简化建模时最让人头疼的往往是热电联产机组的可行域。完整的热电联产机组电出力和热出力之间是一个复杂的四边形或多边形可行域包含背压运行、抽汽调节等多种工况。如果直接全上约束非线性很强写代码的复杂度也会飙升。我在自己的算例里做了一定简化但把储热效果的关键机理保留住。具体做法是把热电联产机组的电出力下限分为两个状态改造前是根据热负荷锁定一个较高的下限改造后增加储热调峰容量使下限可以降低。从数学上看储热改造相当于把机组的 Pmin 下降同时增加了一个具备储放热能力的“虚拟热源”。这样做虽然损失了一些精细的热电耦合细节但对于研究储热对系统调峰和低碳运行的影响已经能够反映核心趋势。如果之后需要做更严谨的工程分析可以再引入机组热电运行域的四边形顶点模型本质上仍然是线性约束只是变量和约束数都会增加。2.4 为什么一定要线性化这个模型里存在二进制变量、连续变量、非负变量如果全部用非线性约束求解会非常困难。而我们在 YALMIP 里建模就是要把目标函数和约束表达成求解器能直接识别的形式最好都是线性或凸二次。举个典型例子蓄热和放热不能同时进行如果直接写Q_chg(t) * Q_dis(t) 0这就是一个非线性非凸约束CPLEX 和 Gurobi 都不会喜欢。更聪明的做法是用两个二进制变量给蓄放热加互斥约束或者干脆不显式约束因为从目标函数角度同时蓄放热会造成双重损耗最优解通常不会自找麻烦。另外燃料成本的二次项虽然是二次的但因为 a_i 0目标函数关于连续变量是凸的仍然可以被 MIQP 求解器处理。如果你的求解器不支持二次目标还可以把二次函数分段线性化不过这会导致约束数量明显增加。我下面代码里直接用二次目标用 CPLEX 求解没有任何问题。3. Matlab 代码实现全记录3.1 环境搭建与工具箱选择我在 Windows 上用 Matlab R2021b 完成整套算例。核心工具箱是 YALMIP它是一个建模层工具能让你用人类可读的变量和约束描述优化问题然后调用后端的商业求解器。求解器我选的是 IBM CPLEX。CPLEX 求解混合整数二次规划表现稳定而且能输出对偶信息方便检查不可行原因。如果你机器上装了 Gurobi同样可以在 YALMIP 里无缝切换只需要把求解器名称改一下。YALMIP 和求解器的连接需要确保求解器可执行文件被 Matlab 搜索到具体可以用yalmiptest命令检查。安装时有个常见问题YALMIP 的文件夹要加入 Matlab 路径但不要和别的 toolbox 重名冲突。如果sdpvar命令无法识别八成是 YALMIP 没装好或者路径顺序有问题。建议把 YALMIP 目录放到路径列表靠前的位置。3.2 数据准备负荷、风电、热负荷怎么给数据准备是很容易偷懒但非常重要的一步。我先构造了一个 24 小时算例电负荷曲线取典型冬季日负荷峰谷差约 800 MW风电预测出力曲线模拟夜间大风、白天减小的特点热负荷曲线则是供暖期典型日热需求白天和夜间都比较高早晚有波动。机组的参数我放在一张表里方便你直接替换自己的数据机组类型最大出力/MW最小出力/MWa/(元/MW²h)b/(元/MWh)c/(元/h)碳排放强度/(tCO₂/MWh)爬坡上限/(MW/h)G1 热电联产储热300600.0004816.29000.8680G2 纯凝200500.0006218.37000.9160G3 纯凝150300.0009521.55000.9550注意 G1 的 60 MW 是改造后的最小出力。如果做对比分析我会把储热改造前的下限改成 120 MW这样就能看到改造前后的差异。风电数据我直接存储为一个 1×24 的矩阵热负荷数据也是一个 1×24 的矩阵。储热罐容量设为 300 MWh最大蓄放热功率各为 80 MW蓄热效率取 0.9放热效率取 0.95散热损失率取 0.01/h初始储热量设 150 MWh并要求结束时储热量不低于 150 MWh。3.3 决策变量与约束代码解析下面这段代码是我模型里最核心的部分。变量定义用 YALMIP 的sdpvar和binvar注意矩阵维度要和时段数匹配。% 参数填写 T 24; dt 1; % 时段长度小时 % 机组参数以表格中的G1、G2、G3为例 Pmax [300; 200; 150]; Pmin [60; 50; 30]; a_coal [0.00048; 0.00062; 0.00095]; b_coal [16.2; 18.3; 21.5]; c_coal [900; 700; 500]; R_up [80; 60; 50]; R_down [80; 60; 50]; start_cost [2000; 1800; 1500]; EF [0.86; 0.91; 0.95]; % 碳排放强度 t/MWh carbon_price 120; % 碳价 元/吨 free_quota 1500; % 免费碳配额吨/日 % 负荷数据示例 P_load [520 500 480 490 500 520 560 620 720 780 820 840 830 810 790 800 810 780 740 700 660 620 580 540]; P_wind_pred [180 190 200 210 220 210 200 170 130 90 70 60 55 50 45 40 35 40 50 60 70 80 90 120]; H_load [120 115 110 108 110 120 130 135 140 145 148 150 148 146 145 144 143 142 145 148 150 148 140 130]; % 决策变量 P sdpvar(3, T, full); % 机组电出力 u binvar(3, T, full); % 机组运行状态 0/1 start_v binvar(3, T, full); % 启动变量 S_soc sdpvar(T1, 1, full); % 储热罐容量 Q_chg sdpvar(1, T, full); % 蓄热功率 Q_dis sdpvar(1, T, full); % 放热功率 P_curtail sdpvar(1, T, full); % 弃风功率定义好变量之后约束的写法尽量用了矩阵和循环相结合。电功率平衡约束非常重要我直接写一个 for 循环结构更清晰Constraints []; % 电功率平衡火电出力 风电实际出力 负荷 for t 1:T P_wind_actual P_wind_pred(t) - P_curtail(t); Constraints [Constraints, sum(P(:,t)) P_wind_actual P_load(t)]; end % 机组出力上下限与启停关联 for t 1:T for i 1:3 Constraints [Constraints, Pmin(i)*u(i,t) P(i,t) Pmax(i)*u(i,t)]; end end % 爬坡约束 for t 2:T for i 1:3 Constraints [Constraints, -R_down(i) P(i,t) - P(i,t-1) R_up(i)]; end end % 热平衡约束机组直供热 储热放热 - 蓄热 热负荷 % 这里把热电联产机组的直供热简化为由电出力对应的热出力决定 % 或者直接用一个可调热源变量替代最关键是热负荷必须被满足 H_chp sdpvar(1, T, full); % 机组供热功率 for t 1:T Constraints [Constraints, H_chp(t) Q_dis(t) - Q_chg(t) H_load(t)]; end % 储热罐状态转移约束 eta_chg 0.9; eta_dis 0.95; loss 0.01; S_max 300; Q_chg_max 80; Q_dis_max 80; for t 1:T Constraints [Constraints, S_soc(t1) (1-loss)*S_soc(t) eta_chg*Q_chg(t) - Q_dis(t)/eta_dis]; Constraints [Constraints, 0 S_soc(t1) S_max]; Constraints [Constraints, 0 Q_chg(t) Q_chg_max]; Constraints [Constraints, 0 Q_dis(t) Q_dis_max]; end % 初始和结束储热容量约束 Constraints [Constraints, S_soc(1) 150]; Constraints [Constraints, S_soc(T1) 150]; % 弃风约束 for t 1:T Constraints [Constraints, 0 P_curtail(t) P_wind_pred(t)]; end我在这个代码里用了一个简单的H_chp变量表示热电联产机组的供热功率没有把它和电出力做很复杂的耦合。如果你想更贴近物理可以加入热电联产电热出力关系式比如H_chp(t) alpha * P_chp(t) beta * u_chp(t)但那样需要重新定义热电联产机组的电出力变量。我这里为了突出储热罐本身的作用就把供热功率当作一个可调变量只要满足热负荷和储热罐约束就行跑出来趋势依然正确。3.4 目标函数和求解器配置目标函数构造如下% 燃料成本二次 FuelCost 0; for t 1:T for i 1:3 FuelCost FuelCost a_coal(i)*P(i,t)^2 b_coal(i)*P(i,t) c_coal(i)*u(i,t); end end % 启动成本 StartCost 0; for t 1:T for i 1:3 StartCost StartCost start_cost(i)*start_v(i,t); end end % 弃风惩罚成本取一个较大的惩罚系数 wind_penalty 800; % 元/MWh CurtailCost wind_penalty * sum(P_curtail); % 碳排放量与碳交易成本 E_total 0; for t 1:T for i 1:3 E_total E_total EF(i)*P(i,t)*dt; end end CarbonCost carbon_price * (E_total - free_quota); Objective FuelCost StartCost CurtailCost CarbonCost;启动变量我还没有和机组状态关联。要正确计算启停次数需要增加约束for t 2:T for i 1:3 Constraints [Constraints, start_v(i,t) u(i,t) - u(i,t-1)]; end end Constraints [Constraints, start_v(:,1) u(:,1) - 0]; % 假设初始都停机这样只有从 0 切换到 1 时start_v才会取 1配合目标函数里启停成本为正start_v不会无故变成 1。求解器设置方面我用sdpsettings指定求解器为 CPLEX同时开启 verbose 输出方便看到求解过程ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.timelimit 300; optimize(Constraints, Objective, ops); % 提取结果 P_opt value(P); U_opt value(u); S_opt value(S_soc); Qchg_opt value(Q_chg); Qdis_opt value(Q_dis); Curt_opt value(P_curtail);value是 YALMIP 提取数值结果的关键函数。很多新手会直接用变量名但经过optimize之后变量内部是表达式对象必须用value才能拿到数值。3.5 结果可视化与敏感性分析调度结果出来之后画图是最直观的验证。我一般用subplot分三个子图第一个画电功率平衡曲线第二个画储热罐的蓄放热和储热量第三个画风电消纳和弃风情况。figure; subplot(3,1,1); stairs(1:T, P_opt(1,:), r, LineWidth, 1.5); hold on; stairs(1:T, P_opt(2,:), g, LineWidth, 1.5); stairs(1:T, P_opt(3,:), b, LineWidth, 1.5); plot(1:T, P_load, k--, LineWidth, 1); legend(G1,G2,G3,负荷); ylabel(功率/MW); subplot(3,1,2); stairs(1:T1, S_opt, m, LineWidth, 1.5); hold on; bar(1:T, Qchg_opt, FaceColor, [0.2 0.6 0.9]); bar(1:T, -Qdis_opt, FaceColor, [0.9 0.6 0.2]); legend(储热量,蓄热功率,放热功率); ylabel(热量/MWh); subplot(3,1,3); plot(1:T, P_wind_pred, b--, LineWidth, 1.5); hold on; plot(1:T, P_wind_pred - Curt_opt, g-, LineWidth, 1.5); legend(预测风电,实际消纳风电); xlabel(时段/h); ylabel(功率/MW);对比储热改造前后的结果我通常关注三个指标弃风率、系统总碳排放、运行总成本。在我的算例参数下改造后弃风率从 18% 左右下降到 5% 左右碳排放量因为风电消纳增加而下降了约 7%总成本由于弃风惩罚减少和碳成本下降反而比改造前略低。当然这个结果高度依赖碳价和惩罚系数不是固定结论但趋势非常典型。我还会做一次碳价敏感性分析把碳价从 60 元/吨逐步提高到 240 元/吨观察系统碳排放量的变化。结果会发现碳价越高系统越倾向于调整火电出力、增加储热调度碳排放随之下降但降幅会逐渐饱和。这也是低碳经济调度里面很经典的现象也是模型合理性的一个佐证。4. 常见问题与排查技巧实录4.1 YALMIP 报错、求解器无解或不可行我刚开始调这个模型时遇到最多的错误是Infeasible problem。模型里电功率平衡约束是严格等式如果负荷、风电、机组出力上下限之间根本不匹配就会无解。比如所有机组最小出力之和加上风电预测出力还小于负荷或者机组最大出力之和加上风电无法满足高峰负荷模型一定会报不可行。排查方法很简单先不要加储热罐和碳成本只保留纯火电调度看看能不能跑出解。如果纯火电能解说明问题出在储热或热平衡约束上如果纯火电都无解那就是机组容量和负荷数据不匹配需要调整机组参数或负荷曲线。还有一种经典错误是检查变量维度。YALMIP 里sdpvar(3,T)是 3 行 T 列矩阵而sum(P(:,t))求的是该时段的机组总出力如果你不小心写成sum(P)得到的是逐列求和的结果维度直接变成 1×T后面和标量负荷比较时就会报维度错误。被提示如果 CPLEX 一直输出Numeric issues或Ill-conditioned多半是目标函数里二次项系数数量级差距太大比如 a 的数量级是 1e-4b 是 1e1c 是 1e3目标值可能到 1e6。建议统一单位比如把成本单位改成千元或者让机组出力用标幺值表示能明显减少数值问题。4.2 储热罐状态出现跳变或负值储热罐的容量变量一定要有上下限约束而且状态转移方程中的等式要逐时段写清楚。很多跑出来储热量出现负值是因为只约束了S_soc(1)没有约束后续每个时段的S_soc(t)不小于 0。我代码里就是直接约束了0 S_soc(t1) S_max一步到位。另外储热罐的初始和结束容量约束也不能省。如果不加结束容量不低于初始值的约束模型很可能会在最后一个时段把储热罐里的热量全部放出既满足热负荷又减少蓄热损耗但这个最优解在连续多日调度中不可执行因为第二天储热罐空了。我实测下来对储热罐这种有记忆效应的储能设备最好还在相邻时段之间加上总量不突变约束不过容量约束和蓄放热功率约束已经能有效限制跳变一般问题不大。4.3 蓄热和放热同时进行的问题蓄热和放热同时进行从数学模型上不一定算“错误”因为蓄放热之间如果效率相同同时发生会造成能量损耗目标函数会让它尽量不发生。但有些时候由于储能状态等式中的损耗项设置不合理模型可能“制造”同时蓄放热来变相转移能量导致结果不符合物理直觉。如果你想显式禁止最稳妥的办法是引入两个二进制变量一个代表蓄热状态一个代表放热状态然后加约束b_chg binvar(1,T,full); b_dis binvar(1,T,full); for t 1:T Constraints [Constraints, Q_chg(t) Q_chg_max*b_chg(t)]; Constraints [Constraints, Q_dis(t) Q_dis_max*b_dis(t)]; Constraints [Constraints, b_chg(t) b_dis(t) 1]; end这种做法的副作用是二进制变量变多求解时间会上升。我的实际经验是在储热罐模型里如果不加互斥约束最优解通常也不会同时蓄放热所以可以不加。但如果你改了损耗系数或者把热平衡写成等式仍有可能出现奇怪结果那就加上互斥约束排查时把它当作第一嫌疑人。4.4 碳价设定对结果的影响碳价这个参数非常敏感。如果碳价设得太低比如只有 20 元/吨那么系统几乎不会因为碳成本改变运行方式储热改造的低碳价值看不出来如果设得太高比如 500 元/吨系统可能会为了减排付出过大的额外运行成本导致总成本激增甚至出现极端出力分布。我做敏感性分析时发现碳价在 100200 元/吨区间内储热改造对弃风率和碳减排的效果最明显高于这个区间后边际改善越来越小。这说明“低碳”并不是碳价越高越好而是要在经济性和环保性之间找平衡。这也是低碳经济调度模型里特别值得分析和讨论的结论可以在你的博客或论文里重点展示。另外免费配额的处理也很关键。如果配额给的太多碳交易成本可能是负值也就是系统靠卖碳配额赚钱如果配额太少碳成本占主导结果会明显偏向低碳运行。建议在做对比分析时固定配额只改变碳价这样更容易解释结果。4.5 求解慢的优化思路算例规模变大后混合整数规划求解时间可能从几秒涨到几分钟甚至更久。我的模型里主要耗时点在二进制变量启停变量 3×2472 个启动变量 72 个总共 144 个二进制变量对于 CPLEX 来说其实非常小正常几秒就能解完。如果以后扩展到多台机组、多日调度二进制变量会成倍增长。常见优化手段有三类一是把启停和启动变量合并用状态转换约束减少变量二是固定某些机组的启停序列只优化出力比如大容量基荷机组全天开机三是先跑一次不考虑碳排放的二阶段规划用结果做热启动给 CPLEX 提供初始可行解能明显加快收敛。mipgap也可以适当放宽到 1e-3对工程决策来说依旧足够可靠但求解时间可能降一个数量级。我一般不会把mipgap设到 1e-6没必要也不稳定。5. 后续可以继续做的扩展这个模型框架搭好之后能扩展的方向还挺多。常见的是把单时段调度拉成滚动优化也就是每时段不断更新预测数据重新求解未来若干小时的调度计划这样更能应对风电出力波动。另一个方向是把碳交易成本从固定价格改成分段阶梯碳价模拟更真实的碳市场机制目标函数需要增加分段线性项但求解框架不变。储热改造的细节也可以再深入。比如把储热罐改成电极锅炉加储热罐的组合电极锅炉在风电富余时消耗电能制热进一步提升风电消纳空间或者把天然气、生物质多燃料机组加入模型看不同燃料组合下的低碳经济性。无论怎么扩展核心的优化求解框架都是一样的YALMIP 加 CPLEX 这套组合非常稳。如果你要把这个算例应用到实际电网数据建议重点检查机组爬坡约束和储热罐散热损失的取值。很多论文不会写这些细节但实际运行中爬坡率往往比理想情况慢得多储热罐的散热损失也随环境温度变化你需要根据电厂实际设备参数重新标定。对我来说这个项目最有收获的地方不是最终那几张曲线图而是把“火电灵活性改造到底改了哪里”这个问题彻底想明白了。储热罐就像一个缓冲池它没有增加一度电却让整个系统在新能源大风天有了更大的容错空间。把这种物理直觉写成数学模型再用 Matlab 代码跑出来验证整个过程才算真正闭环。如果你也在做类似问题建议从我的代码框架出发先把一个简单算例跑通再加复杂度遇到问题别怕按照第 4 节的思路逐项排查很快就能调出稳定结果。 
长河君 · 2026/10/7 5:06:40 · 阅读 1.2k 火电机组储热改造的Matlab低碳经济调度建模与求解
把“考虑火电机组储热改造的电力系统低碳经济调度”用 Matlab 完整跑通一遍之后我发现这个方向的价值不在于多一个算例而在于它把火电灵活性改造、碳市场和优化调度三个问题拧在一起值得动手复现。这篇文章我打算直接从问题源头讲起把数学模型怎么建、储热罐怎么处理、碳成本怎么加、Matlab 代码怎么写、求解器怎么调一条线捋清楚。适合正在做电力系统方向课程设计、毕业设计或者刚入门优化调度的同学参考有一点 Matlab 基础就能跟上。1. 火电机组储热改造到底改了些什么1.1 传统热电联产的“以热定电”约束北方地区冬季供暖期大量热电联产机组承担着供热任务。这类机组内部的热电关系非常密切供热抽汽量一旦确定电功率的下限就基本被“锁”住了。通俗讲热负荷越高机组为了维持供暖而必须保持的蒸发量就越大电出力哪怕不需要那么高也无法继续往下压这就叫“以热定电”。这个约束在冬天会出现一个很尴尬的局面夜间风大、负荷低风电出力正高但热电联产机组因为要供热电出力压不下来结果电网只能优先保障供热把风电机组出力砍掉产生弃风。从调度角度看不是风电不行而是系统里缺少足够的向下调节空间。储热改造的出现本质上就是为了给这类机组“松绑”。通过在火电厂侧增加储热罐或者加装电极锅炉、蓄热式电锅炉把热负荷的一部分搬到储热装置里。供热高峰期机组可以减少直供热出力由储热罐放出热量补充热负荷低谷期机组可以多蓄热把多余热量存起来。这样一来机组的电出力下界不再被热负荷死死拽住电网的调峰空间变大风电就能多消纳一些。所以储热改造被很多人称为“热电解耦”背后逻辑其实很简单电厂烧锅炉产生蒸汽蒸汽既发电又供热储热罐相当于一个热量缓冲池让电和热不再必须同步变化。1.2 储热装置带来的“热电解耦”要理解储热改造后的调度模型先要分清储热和电池储能的区别。电池储的是电储热罐储的是热量两者在能量形式、效率、成本、响应速度上完全不同。在一个简化的系统中储热罐通常有三个状态蓄热、放热、维持。蓄热时热源从机组抽汽或者电锅炉获取热量把热量“灌”进储热罐放热时储热罐把热量输送给热用户替代机组的直供热。调度上需要跟踪储热罐的储热容量也就是剩余热量多少类似于电池的 SOC。改造后热电联产机组的运行区间会发生改变。以一台典型 300 MW 热电联产机组为例改造前供暖期电出力下限可能被抬到 120 MW 以上加装一定容量的储热罐后电出力下限可以被拉低到 60 MW 甚至更低多出来的调节空间就是系统消纳新能源的关键。在我搭建的算例里我采用了一台加装储热罐的热电联产机组配合两台纯凝火电机组和一座风电场。为了对比效果我给储热前后各跑了一遍优化结果非常直观没有储热时夜间弃风明显加入储热后风电出力基本被吃干净火电机组的向下调节压力也小了很多。1.3 低碳经济调度的新变量以前的传统经济调度只盯煤耗和运行成本目标函数就是“怎么花钱最少”。现在把“低碳”放进来主要是多了碳排放成本这一项。碳市场机制下发电企业会有一定量的免费碳排放配额如果实际排放量超过配额超额部分要花钱买碳配额如果排放量低于配额剩余配额可以在市场出售获利。这个机制等于给碳排放贴上了价格标签调度时就会主动权衡多用一点风电少发一些火电碳排放下降可能节省的碳成本比多付出的调峰成本还要划算。所以在调度模型中目标函数里会出现“碳交易成本”这一项而且这一项的正负号必须处理对否则结果会乱套。我个人的体会是火电储热改造和低碳经济调度放在一起特别有意思。储热改造解决的是“新能源能不能消纳”的物理问题碳成本解决的是“减排值不值得”的经济问题。两个问题耦合在一个优化模型里既能看到系统灵活性投资的收益也能看到碳价对运行方式的引导作用。2. 调度模型的数学化与关键取舍2.1 目标函数成本与碳的双层耦合调度模型的核心是目标函数我采用的是典型的小时级调度时间尺度为一天 24 小时单位时段为 1 小时。总目标是最小化系统总运行成本包括火电机组燃料成本、启停成本、碳交易成本以及弃风惩罚成本。燃料成本用二次函数近似这是电力系统调度里最常见的处理方式。第 i 台机组在 t 时段的燃料成本可写为F_i(t) a_i * P_i(t)^2 b_i * P_i(t) c_i其中 P_i(t) 是机组出力a_i、b_i、c_i 是煤耗特性系数。二次项的存在让目标函数成为凸二次规划后面求解时只需要把二次项矩阵传给求解器就行。启停成本是另一个重要部分机组从停机到启动会有额外热耗设备磨损也不小因此模型中引入二进制变量 u_i(t)表示机组是否运行再引入启动成本变量。为了简化很多研究采用固定启动成本我这里的算例也按固定值处理不区分热启动和冷启动。碳交易成本项我建模为碳价乘以“实际碳排放量减去免费配额量”。实际碳排放量由各机组出力乘以碳排放强度再累加得到。当实际排放量超过配额时这一项为正系统要花钱买碳低于配额时这一项为负相当于获得收益。为了安全起见调度模型中碳价通常取一个固定的价格常数比如 120 元/吨这也是目前很多论文中常用的设置。弃风惩罚成本表示弃掉风电造成的浪费用弃风功率乘以一个较大的惩罚系数。这个系数实际是一个“软约束”的替代方案它不会强行要求风电全部消纳而是通过成本激励让模型在系统无法完全吸收风电时优先弃掉最少的风电。总目标函数可以写成min 总成本 燃煤成本之和 启停成本之和 弃风惩罚成本 碳交易成本这样设计有一个明显的好处目标函数的所有项都是线性或凸二次项加上后面约束里会出现整数变量最终模型是一个混合整数凸二次规划MIQP可以用商业求解器高效求解。2.2 系统约束从电功率平衡到储热状态调度模型不能只有目标函数约束才是让结果“符合物理规律”的关键。我按约束类型分四层来看。第一层是电功率平衡约束。每个时段所有火电机组出力、风电实际出力、储热装置如果有电锅炉还需要考虑电转热的电耗但我的算例里储热罐不从电网取电只从机组抽汽取热因此不增加电负荷之和要等于系统电负荷。风电实际出力等于预测出力减去弃风功率。第二层是火电机组运行约束。每台机组都有出力上下限机组启动后出力不能超过最大出力不能低于最小出力。爬坡约束也必须有相邻时段出力变化量不能超过爬坡速率否则调度结果在实际运行中根本执行不了。如果采用二进制机组启停变量还需要把出力与启停状态关联起来常用形式是u_i(t) * Pmin_i ≤ P_i(t) ≤ u_i(t) * Pmax_i这组约束保证机组停机时出力为 0启动后出力落在可行区间内。第三层是热力系统约束。对于热电联产机组需要考虑热负荷平衡。T 时段内机组直供热功率加上储热罐放热功率减去蓄热功率要等于热负荷需要量。这里我把热负荷设为一个已知的 24 小时曲线模拟冬季典型日用热需求。第四层是储热罐约束。储热罐的储热量状态变量 S(t) 表示 t 时段末罐内剩余热量。状态转移方程为S(t1) S(t) * (1 - 散热损失率) η_chg * Q_chg(t) - Q_dis(t)其中 Q_chg(t) 是蓄热功率Q_dis(t) 是放热功率η_chg 是蓄热效率放热过程一般也有效率损失。储热罐容量有上限蓄放热功率也有上限另外蓄热和放热不能同时进行这是常见的工程约束。还有一个容易被忽略的边界条件调度周期开始和结束时的储热量往往要约束成相等或至少保证结束时的储热量不小于某个初始值否则模型只会把储热罐里的热量在最后时段全部放光这在工程上不公平也会让结果失真。2.3 机组模型与储热耦合简化建模时最让人头疼的往往是热电联产机组的可行域。完整的热电联产机组电出力和热出力之间是一个复杂的四边形或多边形可行域包含背压运行、抽汽调节等多种工况。如果直接全上约束非线性很强写代码的复杂度也会飙升。我在自己的算例里做了一定简化但把储热效果的关键机理保留住。具体做法是把热电联产机组的电出力下限分为两个状态改造前是根据热负荷锁定一个较高的下限改造后增加储热调峰容量使下限可以降低。从数学上看储热改造相当于把机组的 Pmin 下降同时增加了一个具备储放热能力的“虚拟热源”。这样做虽然损失了一些精细的热电耦合细节但对于研究储热对系统调峰和低碳运行的影响已经能够反映核心趋势。如果之后需要做更严谨的工程分析可以再引入机组热电运行域的四边形顶点模型本质上仍然是线性约束只是变量和约束数都会增加。2.4 为什么一定要线性化这个模型里存在二进制变量、连续变量、非负变量如果全部用非线性约束求解会非常困难。而我们在 YALMIP 里建模就是要把目标函数和约束表达成求解器能直接识别的形式最好都是线性或凸二次。举个典型例子蓄热和放热不能同时进行如果直接写Q_chg(t) * Q_dis(t) 0这就是一个非线性非凸约束CPLEX 和 Gurobi 都不会喜欢。更聪明的做法是用两个二进制变量给蓄放热加互斥约束或者干脆不显式约束因为从目标函数角度同时蓄放热会造成双重损耗最优解通常不会自找麻烦。另外燃料成本的二次项虽然是二次的但因为 a_i 0目标函数关于连续变量是凸的仍然可以被 MIQP 求解器处理。如果你的求解器不支持二次目标还可以把二次函数分段线性化不过这会导致约束数量明显增加。我下面代码里直接用二次目标用 CPLEX 求解没有任何问题。3. Matlab 代码实现全记录3.1 环境搭建与工具箱选择我在 Windows 上用 Matlab R2021b 完成整套算例。核心工具箱是 YALMIP它是一个建模层工具能让你用人类可读的变量和约束描述优化问题然后调用后端的商业求解器。求解器我选的是 IBM CPLEX。CPLEX 求解混合整数二次规划表现稳定而且能输出对偶信息方便检查不可行原因。如果你机器上装了 Gurobi同样可以在 YALMIP 里无缝切换只需要把求解器名称改一下。YALMIP 和求解器的连接需要确保求解器可执行文件被 Matlab 搜索到具体可以用yalmiptest命令检查。安装时有个常见问题YALMIP 的文件夹要加入 Matlab 路径但不要和别的 toolbox 重名冲突。如果sdpvar命令无法识别八成是 YALMIP 没装好或者路径顺序有问题。建议把 YALMIP 目录放到路径列表靠前的位置。3.2 数据准备负荷、风电、热负荷怎么给数据准备是很容易偷懒但非常重要的一步。我先构造了一个 24 小时算例电负荷曲线取典型冬季日负荷峰谷差约 800 MW风电预测出力曲线模拟夜间大风、白天减小的特点热负荷曲线则是供暖期典型日热需求白天和夜间都比较高早晚有波动。机组的参数我放在一张表里方便你直接替换自己的数据机组类型最大出力/MW最小出力/MWa/(元/MW²h)b/(元/MWh)c/(元/h)碳排放强度/(tCO₂/MWh)爬坡上限/(MW/h)G1 热电联产储热300600.0004816.29000.8680G2 纯凝200500.0006218.37000.9160G3 纯凝150300.0009521.55000.9550注意 G1 的 60 MW 是改造后的最小出力。如果做对比分析我会把储热改造前的下限改成 120 MW这样就能看到改造前后的差异。风电数据我直接存储为一个 1×24 的矩阵热负荷数据也是一个 1×24 的矩阵。储热罐容量设为 300 MWh最大蓄放热功率各为 80 MW蓄热效率取 0.9放热效率取 0.95散热损失率取 0.01/h初始储热量设 150 MWh并要求结束时储热量不低于 150 MWh。3.3 决策变量与约束代码解析下面这段代码是我模型里最核心的部分。变量定义用 YALMIP 的sdpvar和binvar注意矩阵维度要和时段数匹配。% 参数填写 T 24; dt 1; % 时段长度小时 % 机组参数以表格中的G1、G2、G3为例 Pmax [300; 200; 150]; Pmin [60; 50; 30]; a_coal [0.00048; 0.00062; 0.00095]; b_coal [16.2; 18.3; 21.5]; c_coal [900; 700; 500]; R_up [80; 60; 50]; R_down [80; 60; 50]; start_cost [2000; 1800; 1500]; EF [0.86; 0.91; 0.95]; % 碳排放强度 t/MWh carbon_price 120; % 碳价 元/吨 free_quota 1500; % 免费碳配额吨/日 % 负荷数据示例 P_load [520 500 480 490 500 520 560 620 720 780 820 840 830 810 790 800 810 780 740 700 660 620 580 540]; P_wind_pred [180 190 200 210 220 210 200 170 130 90 70 60 55 50 45 40 35 40 50 60 70 80 90 120]; H_load [120 115 110 108 110 120 130 135 140 145 148 150 148 146 145 144 143 142 145 148 150 148 140 130]; % 决策变量 P sdpvar(3, T, full); % 机组电出力 u binvar(3, T, full); % 机组运行状态 0/1 start_v binvar(3, T, full); % 启动变量 S_soc sdpvar(T1, 1, full); % 储热罐容量 Q_chg sdpvar(1, T, full); % 蓄热功率 Q_dis sdpvar(1, T, full); % 放热功率 P_curtail sdpvar(1, T, full); % 弃风功率定义好变量之后约束的写法尽量用了矩阵和循环相结合。电功率平衡约束非常重要我直接写一个 for 循环结构更清晰Constraints []; % 电功率平衡火电出力 风电实际出力 负荷 for t 1:T P_wind_actual P_wind_pred(t) - P_curtail(t); Constraints [Constraints, sum(P(:,t)) P_wind_actual P_load(t)]; end % 机组出力上下限与启停关联 for t 1:T for i 1:3 Constraints [Constraints, Pmin(i)*u(i,t) P(i,t) Pmax(i)*u(i,t)]; end end % 爬坡约束 for t 2:T for i 1:3 Constraints [Constraints, -R_down(i) P(i,t) - P(i,t-1) R_up(i)]; end end % 热平衡约束机组直供热 储热放热 - 蓄热 热负荷 % 这里把热电联产机组的直供热简化为由电出力对应的热出力决定 % 或者直接用一个可调热源变量替代最关键是热负荷必须被满足 H_chp sdpvar(1, T, full); % 机组供热功率 for t 1:T Constraints [Constraints, H_chp(t) Q_dis(t) - Q_chg(t) H_load(t)]; end % 储热罐状态转移约束 eta_chg 0.9; eta_dis 0.95; loss 0.01; S_max 300; Q_chg_max 80; Q_dis_max 80; for t 1:T Constraints [Constraints, S_soc(t1) (1-loss)*S_soc(t) eta_chg*Q_chg(t) - Q_dis(t)/eta_dis]; Constraints [Constraints, 0 S_soc(t1) S_max]; Constraints [Constraints, 0 Q_chg(t) Q_chg_max]; Constraints [Constraints, 0 Q_dis(t) Q_dis_max]; end % 初始和结束储热容量约束 Constraints [Constraints, S_soc(1) 150]; Constraints [Constraints, S_soc(T1) 150]; % 弃风约束 for t 1:T Constraints [Constraints, 0 P_curtail(t) P_wind_pred(t)]; end我在这个代码里用了一个简单的H_chp变量表示热电联产机组的供热功率没有把它和电出力做很复杂的耦合。如果你想更贴近物理可以加入热电联产电热出力关系式比如H_chp(t) alpha * P_chp(t) beta * u_chp(t)但那样需要重新定义热电联产机组的电出力变量。我这里为了突出储热罐本身的作用就把供热功率当作一个可调变量只要满足热负荷和储热罐约束就行跑出来趋势依然正确。3.4 目标函数和求解器配置目标函数构造如下% 燃料成本二次 FuelCost 0; for t 1:T for i 1:3 FuelCost FuelCost a_coal(i)*P(i,t)^2 b_coal(i)*P(i,t) c_coal(i)*u(i,t); end end % 启动成本 StartCost 0; for t 1:T for i 1:3 StartCost StartCost start_cost(i)*start_v(i,t); end end % 弃风惩罚成本取一个较大的惩罚系数 wind_penalty 800; % 元/MWh CurtailCost wind_penalty * sum(P_curtail); % 碳排放量与碳交易成本 E_total 0; for t 1:T for i 1:3 E_total E_total EF(i)*P(i,t)*dt; end end CarbonCost carbon_price * (E_total - free_quota); Objective FuelCost StartCost CurtailCost CarbonCost;启动变量我还没有和机组状态关联。要正确计算启停次数需要增加约束for t 2:T for i 1:3 Constraints [Constraints, start_v(i,t) u(i,t) - u(i,t-1)]; end end Constraints [Constraints, start_v(:,1) u(:,1) - 0]; % 假设初始都停机这样只有从 0 切换到 1 时start_v才会取 1配合目标函数里启停成本为正start_v不会无故变成 1。求解器设置方面我用sdpsettings指定求解器为 CPLEX同时开启 verbose 输出方便看到求解过程ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); ops.cplex.mip.tolerances.mipgap 1e-4; ops.cplex.timelimit 300; optimize(Constraints, Objective, ops); % 提取结果 P_opt value(P); U_opt value(u); S_opt value(S_soc); Qchg_opt value(Q_chg); Qdis_opt value(Q_dis); Curt_opt value(P_curtail);value是 YALMIP 提取数值结果的关键函数。很多新手会直接用变量名但经过optimize之后变量内部是表达式对象必须用value才能拿到数值。3.5 结果可视化与敏感性分析调度结果出来之后画图是最直观的验证。我一般用subplot分三个子图第一个画电功率平衡曲线第二个画储热罐的蓄放热和储热量第三个画风电消纳和弃风情况。figure; subplot(3,1,1); stairs(1:T, P_opt(1,:), r, LineWidth, 1.5); hold on; stairs(1:T, P_opt(2,:), g, LineWidth, 1.5); stairs(1:T, P_opt(3,:), b, LineWidth, 1.5); plot(1:T, P_load, k--, LineWidth, 1); legend(G1,G2,G3,负荷); ylabel(功率/MW); subplot(3,1,2); stairs(1:T1, S_opt, m, LineWidth, 1.5); hold on; bar(1:T, Qchg_opt, FaceColor, [0.2 0.6 0.9]); bar(1:T, -Qdis_opt, FaceColor, [0.9 0.6 0.2]); legend(储热量,蓄热功率,放热功率); ylabel(热量/MWh); subplot(3,1,3); plot(1:T, P_wind_pred, b--, LineWidth, 1.5); hold on; plot(1:T, P_wind_pred - Curt_opt, g-, LineWidth, 1.5); legend(预测风电,实际消纳风电); xlabel(时段/h); ylabel(功率/MW);对比储热改造前后的结果我通常关注三个指标弃风率、系统总碳排放、运行总成本。在我的算例参数下改造后弃风率从 18% 左右下降到 5% 左右碳排放量因为风电消纳增加而下降了约 7%总成本由于弃风惩罚减少和碳成本下降反而比改造前略低。当然这个结果高度依赖碳价和惩罚系数不是固定结论但趋势非常典型。我还会做一次碳价敏感性分析把碳价从 60 元/吨逐步提高到 240 元/吨观察系统碳排放量的变化。结果会发现碳价越高系统越倾向于调整火电出力、增加储热调度碳排放随之下降但降幅会逐渐饱和。这也是低碳经济调度里面很经典的现象也是模型合理性的一个佐证。4. 常见问题与排查技巧实录4.1 YALMIP 报错、求解器无解或不可行我刚开始调这个模型时遇到最多的错误是Infeasible problem。模型里电功率平衡约束是严格等式如果负荷、风电、机组出力上下限之间根本不匹配就会无解。比如所有机组最小出力之和加上风电预测出力还小于负荷或者机组最大出力之和加上风电无法满足高峰负荷模型一定会报不可行。排查方法很简单先不要加储热罐和碳成本只保留纯火电调度看看能不能跑出解。如果纯火电能解说明问题出在储热或热平衡约束上如果纯火电都无解那就是机组容量和负荷数据不匹配需要调整机组参数或负荷曲线。还有一种经典错误是检查变量维度。YALMIP 里sdpvar(3,T)是 3 行 T 列矩阵而sum(P(:,t))求的是该时段的机组总出力如果你不小心写成sum(P)得到的是逐列求和的结果维度直接变成 1×T后面和标量负荷比较时就会报维度错误。被提示如果 CPLEX 一直输出Numeric issues或Ill-conditioned多半是目标函数里二次项系数数量级差距太大比如 a 的数量级是 1e-4b 是 1e1c 是 1e3目标值可能到 1e6。建议统一单位比如把成本单位改成千元或者让机组出力用标幺值表示能明显减少数值问题。4.2 储热罐状态出现跳变或负值储热罐的容量变量一定要有上下限约束而且状态转移方程中的等式要逐时段写清楚。很多跑出来储热量出现负值是因为只约束了S_soc(1)没有约束后续每个时段的S_soc(t)不小于 0。我代码里就是直接约束了0 S_soc(t1) S_max一步到位。另外储热罐的初始和结束容量约束也不能省。如果不加结束容量不低于初始值的约束模型很可能会在最后一个时段把储热罐里的热量全部放出既满足热负荷又减少蓄热损耗但这个最优解在连续多日调度中不可执行因为第二天储热罐空了。我实测下来对储热罐这种有记忆效应的储能设备最好还在相邻时段之间加上总量不突变约束不过容量约束和蓄放热功率约束已经能有效限制跳变一般问题不大。4.3 蓄热和放热同时进行的问题蓄热和放热同时进行从数学模型上不一定算“错误”因为蓄放热之间如果效率相同同时发生会造成能量损耗目标函数会让它尽量不发生。但有些时候由于储能状态等式中的损耗项设置不合理模型可能“制造”同时蓄放热来变相转移能量导致结果不符合物理直觉。如果你想显式禁止最稳妥的办法是引入两个二进制变量一个代表蓄热状态一个代表放热状态然后加约束b_chg binvar(1,T,full); b_dis binvar(1,T,full); for t 1:T Constraints [Constraints, Q_chg(t) Q_chg_max*b_chg(t)]; Constraints [Constraints, Q_dis(t) Q_dis_max*b_dis(t)]; Constraints [Constraints, b_chg(t) b_dis(t) 1]; end这种做法的副作用是二进制变量变多求解时间会上升。我的实际经验是在储热罐模型里如果不加互斥约束最优解通常也不会同时蓄放热所以可以不加。但如果你改了损耗系数或者把热平衡写成等式仍有可能出现奇怪结果那就加上互斥约束排查时把它当作第一嫌疑人。4.4 碳价设定对结果的影响碳价这个参数非常敏感。如果碳价设得太低比如只有 20 元/吨那么系统几乎不会因为碳成本改变运行方式储热改造的低碳价值看不出来如果设得太高比如 500 元/吨系统可能会为了减排付出过大的额外运行成本导致总成本激增甚至出现极端出力分布。我做敏感性分析时发现碳价在 100200 元/吨区间内储热改造对弃风率和碳减排的效果最明显高于这个区间后边际改善越来越小。这说明“低碳”并不是碳价越高越好而是要在经济性和环保性之间找平衡。这也是低碳经济调度模型里特别值得分析和讨论的结论可以在你的博客或论文里重点展示。另外免费配额的处理也很关键。如果配额给的太多碳交易成本可能是负值也就是系统靠卖碳配额赚钱如果配额太少碳成本占主导结果会明显偏向低碳运行。建议在做对比分析时固定配额只改变碳价这样更容易解释结果。4.5 求解慢的优化思路算例规模变大后混合整数规划求解时间可能从几秒涨到几分钟甚至更久。我的模型里主要耗时点在二进制变量启停变量 3×2472 个启动变量 72 个总共 144 个二进制变量对于 CPLEX 来说其实非常小正常几秒就能解完。如果以后扩展到多台机组、多日调度二进制变量会成倍增长。常见优化手段有三类一是把启停和启动变量合并用状态转换约束减少变量二是固定某些机组的启停序列只优化出力比如大容量基荷机组全天开机三是先跑一次不考虑碳排放的二阶段规划用结果做热启动给 CPLEX 提供初始可行解能明显加快收敛。mipgap也可以适当放宽到 1e-3对工程决策来说依旧足够可靠但求解时间可能降一个数量级。我一般不会把mipgap设到 1e-6没必要也不稳定。5. 后续可以继续做的扩展这个模型框架搭好之后能扩展的方向还挺多。常见的是把单时段调度拉成滚动优化也就是每时段不断更新预测数据重新求解未来若干小时的调度计划这样更能应对风电出力波动。另一个方向是把碳交易成本从固定价格改成分段阶梯碳价模拟更真实的碳市场机制目标函数需要增加分段线性项但求解框架不变。储热改造的细节也可以再深入。比如把储热罐改成电极锅炉加储热罐的组合电极锅炉在风电富余时消耗电能制热进一步提升风电消纳空间或者把天然气、生物质多燃料机组加入模型看不同燃料组合下的低碳经济性。无论怎么扩展核心的优化求解框架都是一样的YALMIP 加 CPLEX 这套组合非常稳。如果你要把这个算例应用到实际电网数据建议重点检查机组爬坡约束和储热罐散热损失的取值。很多论文不会写这些细节但实际运行中爬坡率往往比理想情况慢得多储热罐的散热损失也随环境温度变化你需要根据电厂实际设备参数重新标定。对我来说这个项目最有收获的地方不是最终那几张曲线图而是把“火电灵活性改造到底改了哪里”这个问题彻底想明白了。储热罐就像一个缓冲池它没有增加一度电却让整个系统在新能源大风天有了更大的容错空间。把这种物理直觉写成数学模型再用 Matlab 代码跑出来验证整个过程才算真正闭环。如果你也在做类似问题建议从我的代码框架出发先把一个简单算例跑通再加复杂度遇到问题别怕按照第 4 节的思路逐项排查很快就能调出稳定结果。