1. 项目概述与整体思路这几年双碳目标提出来之后电力系统调度的问题越来越有意思。以前做经济调度就是让机组怎么发得便宜怎么来现在不行了还得看碳排放还得考虑新能源消纳。而火电机组作为电网的压舱石光靠改造锅炉、升级汽轮机还不够储热改造是近期特别火的一条技术路线。我接手这个项目的初衷很简单课题组需要一套能够同时体现“储热改造”和“低碳经济调度”耦合关系的Matlab代码用来做仿真验证和论文支撑。标题里这几个字拆开看都不陌生但真要在一个模型里把它们串起来会遇到不少实际问题比如储热罐的数学模型怎么建才不冗余、碳排放成本怎么折算成目标函数里的一个项、混合整数线性规划MILP能不能稳定求解。这篇文章就把我踩过的坑、试过的解法、最后跑通的代码框架完整梳理一遍。适合什么人看如果你是电力系统方向的研究生或者在做火电灵活性改造的工程师手里有Matlab基础想快速搭建一个“含储热火电的低碳经济调度”模型这篇文章可以直接给你一套能跑的代码骨架。如果你只是对概念感兴趣也可以看前面的建模思路代码部分跳过不影响理解整体逻辑。2. 储热改造与低碳经济调度的模型构建2.1 火电机组储热改造的物理模型先说储热改造到底改的是什么。传统火电机组是“热电联产”模式供热季机组出力受热负荷约束电负荷调节范围很窄也就是常说的“以热定电”。储热改造就是在机组的热力系统里加一个蓄热罐用电低谷时把多余的热量储起来用电高峰时再释放出来供热。这样一来机组的电出力就和热出力解耦了调峰能力大幅提升。在数学模型里储热罐可以简化成一个带输入输出和容量上限的能量存储装置。我用的是一阶动态模型储热功率H_ch(t)表示t时段从机组抽热存入罐内的功率放热功率H_dis(t)表示t时段从罐内释放出来供给热负荷的功率储热罐蓄热量S(t)满足S(t) S(t-1) H_ch(t) * η_ch - H_dis(t) / η_dis容量约束S_min ≤ S(t) ≤ S_max这里η_ch和η_dis分别是储热和放热效率一般取0.9左右。要注意的是储热罐不能同时储热和放热所以需要引入一个0-1整数变量来互斥。我一开始没加这个约束跑出来的结果出现了同时储热和放热的无效工况虽然目标函数值更“好看”了但物理上完全不可行。这个问题在后面的代码实现里我会专门说明怎么处理。改造后的火电机组电出力约束也要重新写。传统热电联产机组的电出力范围是一个关于热出力的函数加了储热罐之后机组对外供热的需求不再直接等于机组热出力而是等于热出力减去储热功率再加上放热功率。这样机组的电出力区间就变得更宽了调峰能力就出来了。2.2 碳排放成本与低碳经济调度的耦合经济调度模型的目标函数传统上就是煤耗成本、启停成本、购电成本这一套。低碳调度需要在这些成本之上再加入碳排放相关的成本项。目前主流做法有两类一类是设置一个碳排放基准额度超出部分需要购买碳配额另一类是直接用碳税每吨碳一个固定价格乘以碳排放量。我项目里用的第一类因为国内碳市场目前的运作逻辑就是配额制论文里也更常见。机组碳排放量怎么算最简单的方法是排放强度法即E_i(t) e_i * P_i(t) * Δt其中e_i是第i台机组的碳排放强度t/MWhP_i(t)是该机组在t时段的电出力。更精细一点的做法是用二次函数拟合碳排放量和出力的关系类似煤耗函数的处理。我建议如果数据允许尽量用二次函数因为碳排放强度和负荷率不是线性关系低负荷率下单位电耗会明显上升线性模型会低估低负荷工况的碳排放。这时调度目标函数就变成了min F Σ(燃料成本 启停成本 碳交易成本) Σ(弃风惩罚)注意这里要有一个权衡的思想。加了碳成本之后调度结果会倾向于让高排放的机组少发让清洁能源多发。但并不是碳排放越少越“经济”因为少发火电可能导致系统运行约束被破坏或者需要付出更大的启停成本。这就是“低碳经济调度”这个名词里两个维度互相博弈的地方。用Matlab做这个模型本质上就是在优化一个多目标加权之后的单目标问题。2.3 约束条件全清单模型跑得稳不稳、结果合不合理一半看目标函数另一半看约束条件写没写全。列一下我最后定稿的约束条件清单供参考系统电功率平衡约束所有机组出力加上风电功率等于负荷需求系统热功率平衡约束机组对外供热量加上储热罐放热量等于热负荷需求机组出力上下限约束带有储热后的电出力范围动态变化机组爬坡约束相邻时段出力变化不能超过爬坡速率限制最小启停时间约束这是MILP模型的难点需要引入状态变量旋转备用约束系统要有一定比例的备用容量应对预测误差风电出力约束实际出力不超过预测值储热罐容量与储放热功率约束这里我最想强调最小启停时间约束。很多初学者会把机组启停简单建模成二进制变量然后忘了最小运行/停机时间约束。结果出来的调度方案是机组频繁启停一分钟开一分钟关现场完全没法执行。我用的处理方式是引入两组状态变量加两组不等式约束标准的MILP建模方法后面代码里我会贴出来。3. 求解环境搭建与算法选型3.1 为什么选MatlabYALMIP这个项目本质上是一个混合整数线性规划问题。有人会问直接写C调CPLEX不就行了或者用Python的PyomoGurobi为什么非要用Matlab我的回答是Matlab在这个场景里的优势不在求解速度而在建模速度和调试体验。YALMIP这个工具箱允许你直接用符号变量表达约束和目标函数代码写起来几乎和数学公式一一对应特别适合我们这种“先有数学模型再转代码”的工作流。你不需要手动把约束转化成矩阵形式YALMIP内部会自动帮你整理成求解器需要的标准形式。而且对课题组和学生党来说Matlab的License一般学校都买了不需要额外申请。Python生态虽然也很好但有的服务器上没有配置好Python环境Matlab的跨平台部署反而更省心。3.2 求解器选型与参数设置YALMIP只是一个建模层真正求解还需要底层的求解器。我优先推荐CPLEX和Gurobi两者对MILP问题的求解性能在业界是最顶尖的。如果手头没有商业求解器也可以用开源的SCIP。我实测下来小规模算例比如10台机组、24时段用SCIP也能跑但到了几百个整数变量的规模CPLEX和Gurobi的求解时间优势就非常明显了。在Matlab里的调用方式很简单ops sdpsettings(solver, cplex, verbose, 2, showprogress, 1); ops.cplex.mip.tolerances.mipgap 0.001;MIP Gap这个参数值得说道说道。默认值是1e-4理论上更精确但求解时间会显著拉长。做仿真对比实验的时候没必要追求绝对最优解把MIP Gap设为0.001甚至0.01求解速度快很多结果在工程上几乎没有差别。我一般用0.001兼顾精度和速度。3.3 为什么不用启发式算法硬算还有人和我说这种调度问题为什么不直接用粒子群算法、遗传算法Matlab自带ga函数写起来不是更省事这个问题我用了很多年才想明白。对于小规模问题启发式算法的确能跑但是有几个致命问题第一结果每次跑都不一样复现性差第二整数变量和约束多了之后启发式算法的收敛性没有理论保障第三也是最重要的论文审稿人看到你用遗传算法解MILP大概率会问你“为什么不用精确求解器”。所以只要问题能建模成MILP而且规模在求解器能力范围内我建议优先用精确算法。储热改造调度问题的变量规模通常在几千到几万个量级对CPLEX这种级别的求解器来说完全在射程之内。4. 核心代码实现与关键模块解析4.1 主程序架构设计代码架构我采用了“数据加载-模型构建-求解-结果处理”四段式。没有把功能堆在一个脚本里而是分成了几个函数模块这样后续改机组参数、改调度时段、换目标函数权重都方便。%% 主程序 main.m clear; clc; close all; %% 1. 加载系统数据 load_data; % 定义机组参数、负荷曲线、风电出力、热负荷曲线 %% 2. 构建优化模型 [model, x] build_model(data); % 返回YALMIP变量和约束集 %% 3. 求解 ops sdpsettings(solver, cplex, verbose, 1); result optimize(model.Constraints, model.Objective, ops); %% 4. 结果分析与可视化 if result.problem 0 output_analysis(x, data); else disp([求解失败: , result.info]); end这种模块化的好处有两个一是每个函数可以单独调试出错了不会拖累整个脚本二是方便做敏感性分析比如想比较“储热改造前”和“储热改造后”的调度结果只要改load_data里的参数再跑一遍就行了。4.2 目标函数代码实现目标函数的代码实现并不复杂但有几个隐藏的坑。先说燃料成本。为了线性化我把煤耗函数做了分段线性化处理。具体说就是把出力区间切成几段每段用一条直线近似原来的二次曲线。分段越多精度越高但变量数也会增加。我实际测试中分成三段精度就足够了。%% 目标函数构建片段 Objective 0; for t 1:T for i 1:N % 燃料成本分段线性化 Objective Objective sum(FC_seg{i}(:, t)); % 启停成本 Objective Objective SU_cost(i) * u_on(i, t) SD_cost(i) * u_off(i, t); % 碳交易成本 Objective Objective carbon_price * (E_emit(i, t) - E_allowance(i, t)); end % 弃风惩罚 Objective Objective curtail_penalty * (wind_forecast(t) - wind_actual(t)); end写这段代码时我最开始犯了一个错误就是忘记了碳交易的成本是“超出配额的部分”才需要付费所以我先算了一个总排放量再减去总配额然后乘以碳价。这个写法在总层面没问题但如果你想做逐机组的碳排放成本分析就得在每台机组层面分别计算不然结果输出的时候看不出每台机组分担了多少碳成本。4.3 储热罐建模与互斥约束储热罐的建模是这里面的核心难点原因是储能模型天然带时间耦合。某一时段的蓄热量会影响后续所有时段的可用容量所以不能把每个时段单独拿出来看。%% 储热罐约束组 S sdpvar(1, T1); % 蓄热量多一个维度方便写初值 H_ch sdpvar(1, T); % 储热功率 H_dis sdpvar(1, T); % 放热功率 b_ch binvar(1, T); % 储热状态1表示储热 b_dis binvar(1, T); % 放热状态1表示放热 % 容量约束 S_min S(t) S_max; % 动态约束 S(t1) S(t) eta_ch * H_ch(t) - H_dis(t) / eta_dis; % 储放热功率上限 H_ch(t) b_ch(t) * CH_max; H_dis(t) b_dis(t) * DIS_max; % 互斥约束 b_ch(t) b_dis(t) 1;这里我想特别强调互斥约束也就是最后一行。如果没有这个约束求解器有可能让同一时段既有储热又有放热相当于储能罐在“空转”对系统没有任何实际帮助但会因为效率损耗而白白增加成本所以目标函数自然会排斥这种情况。等等那实际跑的时候会不会出现同时储放热我发现智能算法比如GA确实会出现因为它的搜索方向不严格朝向可行域内部。但是用CPLEX这类精确求解器加上互斥约束之后就没见过这种情况了因为约束已经彻底锁死了可行性空间。当然S_min、S_max、CH_max这些参数必须用实际数据初始化我见过有人忘了初始化S的值导致整个约束组失效求解结果完全莫名其妙。4.4 机组最小启停时间约束的实现最小启停时间是MILP建模里的经典问题。我用了Yu等人的经典建模方法通过三个状态变量来表达u(i,t)机组i在t时段的运行状态1运行0停机y(i,t)机组i在t时段是否启动z(i,t)机组i在t时段是否停机约束表达如下%% 最小启停时间约束简化形式 for i 1:N for t 2:T % 逻辑关系 u(i, t) - u(i, t-1) y(i, t) - z(i, t); end for t 1:T-MinUp(i)1 sum(u(i, t:tMinUp(i)-1)) MinUp(i) * y(i, t); end for t 1:T-MinDown(i)1 sum(1 - u(i, t:tMinDown(i)-1)) MinDown(i) * z(i, t); end end这段代码原理上不难理解如果机组在t时刻启动了那么它至少要保持MinUp个时段的运行状态若在t时刻停机至少要保持MinDown个时段停运。这样写还有个好处——自动把启动和停机动作解耦逻辑清晰也便于后期加入启动燃料、停机维护等细节。我在实际调试里发现这里的难点往往不是代码而是初始时段的处理。仿真从t1开始但机组的初始状态是开着还是关着连续运行了多少小时也得反映到约束里。如果你忽略了机组的初始状态前几个时段的调度结果可能完全不合理比如一台初始状态为运行的机组第一时段就被错误地停机了。处理办法就是在t小于最小启停时间时用初始状态对约束进行修正。具体代码有点啰嗦但原理就是多写几个“前时段”的辅助约束。4.5 数据准备与参数设置示例数据这块我给你一套能直接跑通的小系统参数10台机组、24时段。火电机组参数我参考了几个公开算例修改了一些数值让它更像实际机组。参数数值范围单位机组数量10台单机容量100 ~ 600MW煤耗系数a0.01 ~ 0.05t/MW²h煤耗系数b20 ~ 50t/MWh煤耗系数c100 ~ 500t/h碳排放强度0.5 ~ 1.0t/MWh碳配额总量系统总排放的85%t碳价30 ~ 80元/t储热罐容量100 ~ 300MWh储/放热功率上限30 ~ 100MW负荷数据、风电预测数据、热负荷数据我建议你根据自己的研究对象合理设定如果是做论文尽量贴近实际电网公开数据。有一个细节值得注意热负荷曲线和电负荷曲线的峰谷特性不一样。北方冬季的热负荷白天和晚上差别不大有时候甚至夜里更高这和电负荷的“晚高峰”并不同步。储热罐的价值恰恰就在于利用这种不同步性——夜间电负荷低但热负荷高机组被迫维持较高出力来供热储热罐可以在白天热负荷低的时候多存热夜间放热从而把机组的电出力压低实现深度调峰。这个逻辑是模型能否体现出“储热改造价值”的关键后面结果分析时你会看得非常明显。5. 求解结果分析与效果验证5.1 改造前后的调度结果对比我把储热改造前后的结果放一起对比最直观的差别体现在两个方面一是弃风率二是系统总碳排放。在引入储热罐之前夜间风电大发时段火电因为供热约束不能压得太低风电只能被舍弃。加了储热罐之后夜间让机组降低电出力、同时把多余热量存进储热罐风电消纳空间一下子就大了。我跑的那个10机系统改造前弃风率约在11%改造后降到了3%左右效果可以说是立竿见影。碳排放方面因为增加了储热罐的充放热损耗机组的发电量并没有减少因为它还是要满足同样的电力负荷只是调整了出力时机。那碳排放为什么还会下降原因是夜间风电出力被充分消纳后部分原本由火电承担的电力被风电替代火电总发电量下降碳排放自然就降下来了。具体数值上我的算例从改造前的日均排放8200吨降到了7700吨左右降幅约6%。这个结果很有意义它说明储热改造降低碳排放的机制不是“让火电更高效”而是“给风电让路”。所以分析结论时一定要强调这个因果关系很多文章写的模糊最后读者看了也不知道储热罐到底起了什么作用。5.2 目标函数各项成本拆解我习惯把结果输出成表格方便直接粘贴到论文里。一般我会输出以下几个关键指标指标改造前改造后变化总运行成本万元328.5312.7-4.8%燃料成本万元265.3252.1-5.0%启停成本万元9.27.8-15.2%碳交易成本万元32.528.6-12.0%弃风惩罚万元21.56.2-71.2%弃风率%11.23.1-8.1%从表中能看出成本下降的主要贡献来自弃风惩罚的大幅减少。这也从经济性角度验证了储热改造的可行性虽然储热罐本身有投资成本和运行维护成本但在调度层面它通过减少弃风带来的系统运行效益已经相当可观。要注意的是这里没有考虑储热改造的一次性投资成本如果要做完整的项目可行性分析还要加上投资回收期等内容。5.3 储热罐运行状态可视化我画了三张图辅助分析第一张是电负荷、火电出力、风电出力、弃风功率的时序曲线第二张是储热罐蓄热量和储放热功率随时间的变化曲线第三张是各机组在24个时段的启停状态甘特图。看储热罐的蓄热量曲线最有意思。你会发现它的变化规律基本是白天热负荷低的时候蓄热量逐渐上升傍晚开始放热夜间维持在一个较高水平凌晨风电大发时又开始蓄热。这个模式和电负荷曲线并不同步而是和“电负荷与热负荷的差值”相关。理解了这个逻辑你就能预测不同负荷场景下储热罐的运行模式也能解释为什么某些文献里会出现储热罐一天多次充放的现象。绘制这些图的Matlab代码不复杂核心就是用stairs或者plot配合legend。我比较喜欢用stairs画蓄热量的变化因为储能系统的蓄热量本质上是分段常值变化的楼梯图比平滑曲线更真实。6. 调试记录与常见问题排查6.1 求解器报错信息与排查思路这个项目里我遇到的报错也算经典整理成一个速查表下次你遇到了直接对照排查就行。报错信息可能原因解决办法Warning: Solver not applicableYALMIP认为问题类型与求解器不匹配检查是否引入了非线性项如两个变量相乘Infeasible problem约束过紧或数据矛盾先去掉部分约束测试可行性逐步加回定位冲突Out of memory变量规模过大减少时段数或机组数检查是否有冗余变量NaN in solution数据中有未初始化变量检查数据加载函数打印关键变量检查是否为NaNMIP gap not converging求解时间不够或问题规模大设置最大求解时间适当放松MIP GapIndex exceeds matrix dimensions循环变量越界检查tMinUp-1这类索引是否超出T最常见的还是Infeasible problem。我排查的经验是先看功率平衡约束和机组出力上下限是否一致。比如我把机组的总最小出力设置得比系统最低负荷还高而系统又没有储能和外来输电通道那这个模型天然就无解。这种问题不是代码写错了而是数据组合本身不可行。调试方法就是先把所有约束注释掉一个一个加回来加上哪个开始无解问题就在哪里。6.2 数值尺度问题一个容易被忽视的大坑还有一个经验想分享给你就是数值尺度问题。最开始我的目标函数里燃料成本动辄几百万而储热罐的蓄热量只有几百两者相差好几个数量级。CPLEX这类求解器在内部处理时要靠容差来判断约束是否满足如果尺度差太远会出现两个问题一是收敛速度极慢二是约束被求解器“忽略”掉。解决办法有两个。第一个是统一单位把所有成本都换算成万元所有功率都换算成百MW这样量级就接近了。第二个是给YALMIP设置数值容差但我觉得不如直接改单位来得干净。我当时把成本项除以10000功率项除以100之后求解时间直接下降了一半以上这个优化非常值得做。6.3 调试技巧小规模测试先行最后建议你拿到一个新模型别一上来就跑完整24时段、10台机组的规模。先做一个3台机组、6时段的玩具模型跑通了再逐渐扩大规模。这样做的好处有三个第一小规模模型的解你可以手算验证能确认模型逻辑正确第二求解速度快改一版代码几分钟就能看到结果第三小规模时储热罐的运行逻辑更容易可视化观察方便你理解调度结果。我自己写代码的习惯是先把所有约束都放到一个大数组里然后加一些disp输出每个约束的数量确认约束数量和你预期是否一致。比如10台机组、24时段最小启停约束理论上应该有10*(24-MinUp1)条左右如果数量对不上一定是有循环边界写错了。这个检查技巧看起来土但确实帮我抓到过好几次隐藏bug。7. 扩展方向与个人心得代码跑通之后很多人会问接下来还能做什么。我的建议是往三个方向扩展第一加入不确定性。风电预测是有误差的把这种误差建模成随机场景就变成了随机规划或者鲁棒优化问题。调度结果会更保守但也更贴合实际运行场景。第二把储能系统换成电化学储能或者抽水蓄能模型框架几乎不用动只需要修改储能设备的约束参数就行。这样一套代码就能覆盖多种储能技术对比写论文时特别有用。第三考虑电网网架约束。现在的模型是单节点模型没考虑线路潮流约束。如果想研究储热改造对线路阻塞的影响需要在模型中引入直流潮流方程那变量和约束都要增加一个维度。就我个人体会这个项目对我最大的启发是储能技术的价值必须在系统层面才能完全体现出来。你要是单独看储热罐本身它只是一个装热水的罐子但放到电力系统调度模型里它就成了连接热力系统和电力系统的关键节点是火电灵活性的放大器。建模的时候不要只盯着储能设备本身的公式要多想想它在整个系统里的角色——是削峰填谷、是促进消纳、还是提供备用不同定位下模型的约束表达会有微妙差别。最后再补充一句Matlab代码写得好不好很大程度上取决于你对YALMIP的熟练程度。建议花一个下午把YALMIP自带的几个例子跑熟理解sdpvar、binvar、optimize这三个核心函数的基本用法后面的项目就顺了。希望这份笔记能帮你少走一些弯路。