我的实际经验是这类“火电储热改造 低碳经济调度”的题目看起来像一篇典型的电力系统优化论文选题但如果只看摘要和公式很容易陷入“模型写出来但跑不动、跑通了但结果不合理”的尴尬境地。真正动手用Matlab实现时难点往往不在优化理论本身而在目标函数怎么线性化、约束怎么避免冗余、以及求解器怎么在可接受时间内给出稳定解。这篇文章我就以“考虑火电机组储热改造的电力系统低碳经济调度”为题把完整的建模思路、Matlab代码实现架构和我在调试中踩过的坑拆开讲。内容偏工程实操适合正在做电力系统优化调度课题的研究生以及需要做灵活性改造方案预研的规划工程师参考。1. 火电储热改造为什么值得做从“以热定电”到热电解耦1.1 热电联产机组为什么消纳不了风电先讲一个很多人都绕不过去的背景。传统火电里的热电机组CHPCombined Heat and Power在供热季必须同时满足电负荷和热负荷而它的电出力下限并不是一个固定值而是由当前热出力决定的热出力越高最小电出力就越高这就是所谓的“以热定电”约束。画一条简化的CHP电热运行特性曲线来看横轴是热出力纵轴是电出力机组的可行域是一个梯形区域。左上边界是最大电出力右下边界是最小电出力而最小电出力会随热出力的增大而抬升。夜间风电大发时电网希望把火电出力压下去但热电机组因为还在供热电出力压不到纯凝机组的深度调峰水平风电就得被限制。我做过一个简单算例一个区域电网里有2台300MW级热电机组、1台600MW纯凝机组和300MW风电装机供热季夜间热负荷约为500MWth。此时两台热机最小电出力合计超过240MW加上纯凝机组按旋转备用要求至少运行120MW风电上网空间只剩不到90MW而夜间风电可能出力到250MW以上弃风率轻松超过15%。1.2 储热改造的工程本质与经济逻辑储热改造解决的就是上面这个矛盾。在热电机组的供热管网上并联一个储热罐相当于在“热源—热用户”之间加了一个缓冲池。当电网需要压火电电出力时热电机组保持锅炉和汽机供热抽汽不变但把一部分热出力送到储热罐而不是直接送向热网从而降低对外供热出力电出力下限也随之下降。这么说可能抽象打个比方原本热电机组必须“即发即用”地满足热负荷就像你工资一到账就必须花完加了储热罐之后你可以把一部分钱存起来当下个月要交房租时再取出来。储热罐让热电机组从“刚性供热”变成了“柔性供热”。从经济角度看储热改造的收益来自三个方面一是减少弃风带来的发电量补偿二是通过降低煤耗率取得燃料节约三是在碳交易机制下减少碳排放配额支出。以北方某300MW级热电机组的典型改造为例储热罐容量按机组最大热出力放4小时设计总投资大约在3000万至5000万元级别如果供热季每天能增加4小时以上的风电消纳空间静态回收期普遍在4到6年。但要注意储热罐本身也有损耗。充放热过程存在热损失储罐放热温度越高、保温越差损耗越大。建模时通常用一个充热效率和放热效率来表征一般在0.9到0.95之间。这个参数直接决定了储热罐“存一度热、放出能用的热”到底有多少很多人初期建模时忽略它导致调度结果过于乐观。2. 低碳经济调度的数学模型目标函数与约束的逐层拆解2.1 目标函数五类成本怎么算低碳经济调度的“低碳”和“经济”两个目标我建议不要做成双目标加权而是直接把碳排放成本化纳入统一的目标函数。这样既符合碳交易的现实机制求解时也不用处理Pareto前沿简化了模型复杂度。目标函数我通常写成这样min C C_fuel C_carbon C_wind C_storage C_start各部分含义如下C_fuel燃料成本常规做法是用二次函数 aP² bP c 拟合机组煤耗特性再乘以标煤价格。二次函数是经典的工程近似拟合精度在工程允许范围内。C_carbon碳排放成本对纯凝机组和热电机组分别按电出力和热出力折算碳排放量再乘以碳价。碳价可以直接设为一个常数也可以做成阶梯碳价来模拟碳市场不同配额约束下的边际价格变化。C_wind弃风惩罚成本风电预测出力与实际上网电量的差值乘以惩罚系数。这个惩罚系数本质上不是真实成本而是为了让模型优先消纳风电而设置的一个“影子价格”取值一般要高于单位煤电成本否则模型宁可弃风也不压火电。C_storage储热运行维护成本按储热罐充放热量计算的维护费用通常取值很小主要是防止模型毫无意义的频繁充放热。C_start机组启停成本风电波动剧烈时模型可能倾向于频繁启停机组来消纳。加启停成本就是为了避免这种不实际的调度行为。2.2 约束体系从电网侧到储热侧约束是模型能跑出合理结果的关键。我粗略地把约束分为四层每一层都有容易漏掉的细节。第一层电功率平衡每个调度时段所有机组电出力之和加上风电上网功率必须等于系统电负荷。这里注意如果模型里包含储能如电化学储能或抽蓄还要把储能充放电功率也放进平衡方程。本文讨论的储热属于热侧设备不直接出现在电功率平衡方程中这一点很容易搞混。第二层机组运行约束包括火电机组出力上下限、爬坡速率约束、最小启停时间约束。对于热电机组还要加入电热耦合可行域约束。最常见的简化建模方式是P_min_chp(t) P_min_conventional beta * H_out(t) / H_max其中 H_out(t) 是热电机组在 t 时段的对外热出力beta 是电出力下限随热出力变化的斜率。H_out 越小电出力下限就越低这就是储热改造调节能力的体现。第三层储热系统约束储热罐的核心约束是能量递推关系S(t1) S(t) eta_ch * Q_ch(t) - Q_dis(t) / eta_disS(t) 是 t 时段储热罐的储热量Q_ch 和 Q_dis 分别为充热和放热功率eta_ch 和 eta_dis 分别是充、放热效率。还要限制储热罐容量上下限、充放热功率上限以及防止同时充放热。第四层热负荷平衡每个时段热电机组对外热出力加上储热罐放热功率、减去充热功率必须等于热网热负荷H_chp(t) Q_dis(t) - Q_ch(t) H_load(t)这个方程是储热改造能够“解耦热电”的关键储热罐放热时热电机组对外热出力可以减少电出力下限随之降低从而腾出更多空间给风电。2.3 从非线性到线性MILP化处理目标函数里有二次项、约束里有双线性项如果直接用非线性求解器处理求解速度和稳定性都很难保证。工程上最常用的做法是把模型转化为混合整数线性规划MILP问题。燃料成本二次函数的线性化我一般用分段线性逼近。把机组出力范围分成若干区间在每个区间内用一条直线近似原成本曲线。分段数取3到5段即可再增加段数对精度的提升很有限但求解时间会显著增长。这里给出一个我自己调试时的参考经验一台600MW机组燃料成本曲线分4段线性化时相对二次函数模型的误差通常在0.5%以内但求解时间分段的整数变量增加带来的影响不大因为分段线性化用的是连续变量的特殊有序集表示不会引入大量0/1变量。对于热电机组可行域的线性化需要注意“电出力下限随热出力变化”这个关系。热出力本身也是连续性变量所以直接用一个线性表达式约束即可不需要额外分段。3. Matlab实现架构与核心代码逻辑3.1 整体框架数据、建模、求解三阶段Matlab代码我建议按三个文件组织而不是把所有代码堆到一个脚本里。这样后续改参数、换场景、做敏感性分析都方便。第一个文件是参数初始化脚本负责定义所有基础数据机组参数、负荷曲线、风电出力场景、储热罐参数、碳价、成本系数。第二个文件是模型构建脚本用Yalmip定义决策变量、目标函数和约束。第三个文件是求解与后处理脚本调用求解器并输出结果表格和图形。这种分离的好处是做敏感性分析时只需要修改参数文件模型文件完全不动做不同场景对比时直接循环调用模型文件即可。我在实际项目里甚至会把这套逻辑再封装成一个函数输出每个时段的机组出力、储热罐SOC和总成本方便批量跑64个典型场景。3.2 核心代码Yalmip建模的关键片段下面给出模型构建脚本中最核心的Yalmip代码片段变量定义和约束构造的逻辑我都写了注释。%% 决策变量定义 % P: 火电机组电出力维度为 n_gen x T % H_chp: 热电机组对外热出力维度为 n_chp x T % Q_ch, Q_dis: 储热罐充、放热功率维度为 1 x T % S: 储热罐储热量维度为 1 x (T1)首末时段储热量固定 P sdpvar(n_gen, T, full); H_chp sdpvar(n_chp, T, full); Q_ch sdpvar(1, T, full); Q_dis sdpvar(1, T, full); S sdpvar(1, T 1, full); P_wind sdpvar(n_wind_farm, T, full); %% 目标函数燃料成本 碳排放成本 弃风惩罚 储热运维成本 objective 0; for t 1:T for g 1:n_gen objective objective coeff_fuel(g,1) * P(g,t)^2 ... coeff_fuel(g,2) * P(g,t) ... coeff_fuel(g,3); end objective objective carbon_price * sum(emission_intensity .* P(:,t)) ... wind_penalty * sum(P_wind_pred(:,t) - P_wind(:,t)) ... storage_om * (Q_ch(t) Q_dis(t)); end %% 约束条件 Constraints []; % 电功率平衡 for t 1:T Constraints [Constraints, sum(P(:,t)) sum(P_wind(:,t)) P_load(t)]; end % 火电出力上下限和爬坡约束 for t 1:T for g 1:n_gen Constraints [Constraints, P_min(g) P(g,t) P_max(g)]; if t 1 Constraints [Constraints, P(g,t) - P(g,t-1) Ramp_up(g)]; Constraints [Constraints, P(g,t-1) - P(g,t) Ramp_down(g)]; end end end % 风电上网功率约束 for t 1:T Constraints [Constraints, zeros(1, n_wind_farm) P_wind(:,t) P_wind_pred(:,t)]; end % 储热罐SOC递推 Constraints [Constraints, S(1) S_init]; Constraints [Constraints, S(T1) S_init]; for t 1:T Constraints [Constraints, S(t1) S(t) eta_ch * Q_ch(t) - Q_dis(t) / eta_dis]; Constraints [Constraints, 0 Q_ch(t) Q_ch_max]; Constraints [Constraints, 0 Q_dis(t) Q_dis_max]; Constraints [Constraints, S_min S(t1) S_max]; end %% 求解设置 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 0.001); % 如果目标函数包含P的平方项Yalmip会自动将其视为二次规划 % 若想用纯线性MILP则需预先对燃料成本做分段线性化。3.3 求解器选型与Yalmip常用坑Matlab环境下我首选Yalmip做建模层因为它屏蔽了不同求解器API的差异换求解器只需要改一行sdpsettings配置。求解器方面如果学校或单位有Gurobi/Cplex许可证优先选它们。MILP求解速度和稳定性最好。没有商业求解器时可以用Matlab自带的intlinprog。但Yalmip对intlinprog的优化力度不如Gurobi大规模算例可能出现求解极慢的情况。或者用HiGHS这是开源求解器里表现比较好的Yalmip新版已经支持。关于储热罐SOC的初值设置这里有个容易忽略的细节如果固定S(1)和S(T1)都等于某个给定初值模型会强制储热罐在调度周期结束时回到初始状态。这样做的合理性在于调度是滚动进行的不能让储热罐在一个周期内“用完即弃”否则下一个周期的调度结果会失真。另外我调试时经常遇到Yalmip报错“Nonlinear constraints not supported”或者“Currently Yalmip does not support ...”。多数情况下是因为某处变量相乘产生了双线性项。排查方法很简单用check(Constraints)检查每一个约束或者直接用yalmip(clear)清空变量空间再重新声明。4. 算例设计与结果解读改造容量、碳价、弃风率的关键影响4.1 算例基础数据为了便于说明我构造了一个典型的区域电网测试系统。基础参数如下项目数值热电机组2台额定电出力300MW最大热出力400MWth纯凝机组1台额定电出力600MW风电场装机300MW系统峰值负荷1200MW热负荷曲线夜间高800MWth白天低500MWth碳价基准100元/吨储热罐容量基准600MWh充放热功率上限150MWth负荷曲线和风电出力曲线取的是典型的冬季供热日数据夜间负荷低但热负荷高风电出力夜间高、白天低。4.2 三个方案的结果对比我分别跑三个方案完全不装储热、储热容量300MWh、储热容量600MWh。在相同负荷和风电场景下结果如下指标无储热储热300MWh储热600MWh弃风率18.6%6.8%3.1%系统总运行成本万元126.5119.2116.7碳排放总量吨7,8427,1756,952需要说明的是这是算例演示的典型结果具体数值会随系统参数变化。但从趋势上可以清楚看到两件事第一储热容量从0增加到300MWh时弃风率大幅下降成本下降明显。这说明在原本“以热定电”约束很强的系统里储热改造的边际效益很高。第二从300MWh增加到600MWh时弃风率只下降了不到4个百分点成本下降幅度也变小了。这说明储热容量存在边际效益递减的特征。原因在于弃风主要发生在夜间几个时段而储热罐最多只能把夜间的热量搬到白天当罐容量已经能覆盖整个夜间热负荷时再扩容就不能增加更多的风电消纳空间了。4.3 储热容量与碳价的敏感性分析我把碳价从50元/吨逐步提高到200元/吨同时观察模型在固定储热容量下的碳排放量变化碳价元/吨碳排放量吨系统总成本万元508,214112.31007,175116.71506,763119.82006,581122.6碳价提升会让模型更积极压低火电出力、增加风电上网但碳排放下降的幅度会逐渐放缓。这是因为系统里的灵活性资源已经用尽继续提高碳价也无法进一步减排只能体现为成本增加。做实际项目时这个拐点很重要它对应着“碳价高到一定程度后单纯靠调度优化无法实现更多减排必须从电源侧或储能侧做增量投资”的结论。5. 从仿真到现场边界条件、调试心得与工程化建议5.1 模型假设在真实工程中的边界仿真模型跑通了只是第一步真正做工程方案预研时必须清楚模型的假设边界。储热罐的散热损失在模型里通常简化成固定效率但实际中散热速率与罐内温度、环境温度、保温层状态都有关。长期调度模拟时如果按固定效率外推会低估冬季连续运行时的热损失。热电机组的电热耦合特性在不同抽汽工况下差异很大。我在建模时用的是一条线性斜率来描述最小电出力与热出力的关系但实际机组的可行域可能是分段折线而且受主蒸汽压力、凝汽器背压影响。如果要做可行性研究级别的方案建议找热机专业人员要机组实际运行特性曲线把可行域改成分段线性约束代码并不需要大幅改动只是多几组约束行。还有风电出力的不确定性。本文静态优化用的是预测曲线但在实际调度中风电预测误差很大。做工程方案时至少要叠加一个“预测误差场景集合”来做鲁棒优化或随机优化否则只看单场景结果容易过度乐观。5.2 调试与排错经验我在Matlab代码调试中遇到过几个比较典型的问题逐个说下我的处理方式。第一个问题是变量维度不匹配导致的约束错误。比如P_wind_pred是一个行向量而P_wind用sdpvar定义了多行T列矩阵直接做减法时Yalmip不报错但会生成一张维度不对的约束矩阵最终结果完全错误。我的经验是每构造完一组约束就调用size(Constraints(i))检查或者直接对变量用相同维度声明。第二个问题是分段线性化后求解变慢。分段数加得太多、或者每个机组的断点选择不恰当会导致MILP模型规模失控。我试过把煤耗曲线分15段线性化模型变量数翻了一倍求解时间从2分钟涨到30分钟而精度只提升了不到0.1%。最终的折中方案是把每台机组分4段并在波动较大的低出力区段加密断点。第三个问题是储热罐的SOC波动异常。我曾经遇到过一个情况SOC曲线在各个时段剧烈振荡看起来像是求解器在“抖”。排查后发现是充热效率和放热效率没有区分或者说我用了同一个变量同时做充放热模型钻了“同时充放热”的空子把热损失当成了收益来源。解决办法是加一个防止同时充放热的约束Q_ch(t) M * u(t)、Q_dis(t) M * (1-u(t))其中u(t)是0/1变量。如果不想引入整数变量也可以直接限制Q_ch(t) * Q_dis(t) 0但这个约束是非线性的并不推荐。第四个问题是mipp间隙设置。默认的求解间隙对这类调度问题通常够了但如果你要输出“谁在哪个时段启停”这类结论必须把mipgap设得严格一些。否则求解器提前终止会留下次优解启停策略可能完全不是实际最优方案。我一般设0.001或更小。5.3 扩展方向这套代码框架其实可以延伸出很多实用变体。比如把单目标碳成本改成多种碳交易机制的对比分析把确定性风电预测扩展成基于场景的随机规划在储热罐之外加入电锅炉、电极锅炉等电转热设备比较不同灵活性改造路线的经济性。我个人最常做的一个扩展是在模型里增加“机组深度调峰状态”的逻辑。热电机组降到某个出力以下时煤耗率和排放强度都会急剧上升这时候继续压火电可能并不划算模型会在碳成本和弃风惩罚之间自动权衡。这个逻辑用一个阶梯约束就能实现但是能让结果更贴近真实调度员的决策行为。回到调度问题本身储热改造的本质就是买一个“时间平移”的能力。Matlab代码的价值在于把这个能力量化出来了让决策者看到每一MWth储热容量能换来多少风电上网空间、多少减排量、多少成本变化。按我上面这套流程来做就算从零开始两个工作日左右也能跑出一套完整可复用的低碳调度模型。