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/MWh），P_i(t)是该机组在t时段的电出力。更精细一点的做法是用二次函数拟合碳排放量和出力的关系，类似煤耗函数的处理。我建议如果数据允许，尽量用二次函数，因为碳排放强度和负荷率不是线性关系，低负荷率下单位电耗会明显上升，线性模型会低估低负荷工况的碳排放。
这时调度目标函数就变成了：
min F = Σ(燃料成本 + 启停成本 + 碳交易成本) + Σ(弃风惩罚)
注意这里要有一个权衡的思想。加了碳成本之后，调度结果会倾向于让高排放的机组少发，让清洁能源多发。但并不是碳排放越少越“经济”，因为少发火电可能导致系统运行约束被破坏，或者需要付出更大的启停成本。这就是“低碳经济调度”这个名词里两个维度互相博弈的地方。用Matlab做这个模型，本质上就是在优化一个多目标加权之后的单目标问题。
2.3 约束条件全清单
模型跑得稳不稳、结果合不合理，一半看目标函数，另一半看约束条件写没写全。列一下我最后定稿的约束条件清单，供参考：
- 系统电功率平衡约束：所有机组出力加上风电功率，等于负荷需求
- 系统热功率平衡约束：机组对外供热量加上储热罐放热量，等于热负荷需求
- 机组出力上下限约束：带有储热后的电出力范围动态变化
- 机组爬坡约束：相邻时段出力变化不能超过爬坡速率限制
- 最小启停时间约束：这是MILP模型的难点，需要引入状态变量
- 旋转备用约束：系统要有一定比例的备用容量应对预测误差
- 风电出力约束：实际出力不超过预测值
- 储热罐容量与储放热功率约束
这里我最想强调最小启停时间约束。很多初学者会把机组启停简单建模成二进制变量，然后忘了最小运行/停机时间约束。结果出来的调度方案是机组频繁启停，一分钟开一分钟关，现场完全没法执行。我用的处理方式是引入两组状态变量加两组不等式约束，标准的MILP建模方法，后面代码里我会贴出来。
3. 求解环境搭建与算法选型
3.1 为什么选Matlab+YALMIP
这个项目本质上是一个混合整数线性规划问题。有人会问：直接写C++调CPLEX不就行了？或者用Python的Pyomo+Gurobi？为什么非要用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, T+1); % 蓄热量，多一个维度方便写初值 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(t+1) = 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:t+MinUp(i)-1)) >= MinUp(i) * y(i, t); end for t = 1:T-MinDown(i)+1 sum(1 - u(i, t:t+MinDown(i)-1)) >= MinDown(i) * z(i, t); end end这段代码原理上不难理解：如果机组在t时刻启动了，那么它至少要保持MinUp个时段的运行状态；若在t时刻停机，至少要保持MinDown个时段停运。这样写还有个好处——自动把启动和停机动作解耦，逻辑清晰，也便于后期加入启动燃料、停机维护等细节。
我在实际调试里发现，这里的难点往往不是代码，而是初始时段的处理。仿真从t=1开始，但机组的初始状态（是开着还是关着，连续运行了多少小时）也得反映到约束里。如果你忽略了机组的初始状态，前几个时段的调度结果可能完全不合理，比如一台初始状态为运行的机组，第一时段就被错误地停机了。处理办法就是在t小于最小启停时间时，用初始状态对约束进行修正。具体代码有点啰嗦，但原理就是多写几个“前时段”的辅助约束。
4.5 数据准备与参数设置示例
数据这块我给你一套能直接跑通的小系统参数，10台机组、24时段。火电机组参数我参考了几个公开算例，修改了一些数值让它更像实际机组。
| 参数 | 数值范围 | 单位 |
|---|---|---|
| 机组数量 | 10 | 台 |
| 单机容量 | 100 ~ 600 | MW |
| 煤耗系数a | 0.01 ~ 0.05 | t/MW²h |
| 煤耗系数b | 20 ~ 50 | t/MWh |
| 煤耗系数c | 100 ~ 500 | t/h |
| 碳排放强度 | 0.5 ~ 1.0 | t/MWh |
| 碳配额总量 | 系统总排放的85% | t |
| 碳价 | 30 ~ 80 | 元/t |
| 储热罐容量 | 100 ~ 300 | MWh |
| 储/放热功率上限 | 30 ~ 100 | MW |
负荷数据、风电预测数据、热负荷数据我建议你根据自己的研究对象合理设定，如果是做论文，尽量贴近实际电网公开数据。
有一个细节值得注意：热负荷曲线和电负荷曲线的峰谷特性不一样。北方冬季的热负荷白天和晚上差别不大，有时候甚至夜里更高，这和电负荷的“晚高峰”并不同步。储热罐的价值恰恰就在于利用这种不同步性——夜间电负荷低但热负荷高，机组被迫维持较高出力来供热，储热罐可以在白天热负荷低的时候多存热，夜间放热，从而把机组的电出力压低，实现深度调峰。这个逻辑是模型能否体现出“储热改造价值”的关键，后面结果分析时你会看得非常明显。
5. 求解结果分析与效果验证
5.1 改造前后的调度结果对比
我把储热改造前后的结果放一起对比，最直观的差别体现在两个方面：一是弃风率，二是系统总碳排放。
在引入储热罐之前，夜间风电大发时段，火电因为供热约束不能压得太低，风电只能被舍弃。加了储热罐之后，夜间让机组降低电出力、同时把多余热量存进储热罐，风电消纳空间一下子就大了。我跑的那个10机系统，改造前弃风率约在11%，改造后降到了3%左右，效果可以说是立竿见影。
碳排放方面，因为增加了储热罐的充放热损耗，机组的发电量并没有减少（因为它还是要满足同样的电力负荷），只是调整了出力时机。那碳排放为什么还会下降？原因是夜间风电出力被充分消纳后，部分原本由火电承担的电力被风电替代，火电总发电量下降，碳排放自然就降下来了。具体数值上，我的算例从改造前的日均排放8200吨降到了7700吨左右，降幅约6%。
这个结果很有意义，它说明储热改造降低碳排放的机制不是“让火电更高效”，而是“给风电让路”。所以分析结论时一定要强调这个因果关系，很多文章写的模糊，最后读者看了也不知道储热罐到底起了什么作用。
5.2 目标函数各项成本拆解
我习惯把结果输出成表格，方便直接粘贴到论文里。一般我会输出以下几个关键指标：
| 指标 | 改造前 | 改造后 | 变化 |
|---|---|---|---|
| 总运行成本（万元） | 328.5 | 312.7 | -4.8% |
| 燃料成本（万元） | 265.3 | 252.1 | -5.0% |
| 启停成本（万元） | 9.2 | 7.8 | -15.2% |
| 碳交易成本（万元） | 32.5 | 28.6 | -12.0% |
| 弃风惩罚（万元） | 21.5 | 6.2 | -71.2% |
| 弃风率（%） | 11.2 | 3.1 | -8.1% |
从表中能看出，成本下降的主要贡献来自弃风惩罚的大幅减少。这也从经济性角度验证了储热改造的可行性：虽然储热罐本身有投资成本和运行维护成本，但在调度层面，它通过减少弃风带来的系统运行效益已经相当可观。要注意的是，这里没有考虑储热改造的一次性投资成本，如果要做完整的项目可行性分析，还要加上投资回收期等内容。
5.3 储热罐运行状态可视化
我画了三张图辅助分析：第一张是电负荷、火电出力、风电出力、弃风功率的时序曲线；第二张是储热罐蓄热量和储放热功率随时间的变化曲线；第三张是各机组在24个时段的启停状态甘特图。
看储热罐的蓄热量曲线最有意思。你会发现它的变化规律基本是：白天热负荷低的时候蓄热量逐渐上升，傍晚开始放热，夜间维持在一个较高水平，凌晨风电大发时又开始蓄热。这个模式和电负荷曲线并不同步，而是和“电负荷与热负荷的差值”相关。理解了这个逻辑，你就能预测不同负荷场景下储热罐的运行模式，也能解释为什么某些文献里会出现储热罐一天多次充放的现象。
绘制这些图的Matlab代码不复杂，核心就是用stairs或者plot配合legend。我比较喜欢用stairs画蓄热量的变化，因为储能系统的蓄热量本质上是分段常值变化的，楼梯图比平滑曲线更真实。
6. 调试记录与常见问题排查
6.1 求解器报错信息与排查思路
这个项目里我遇到的报错也算经典，整理成一个速查表，下次你遇到了直接对照排查就行。
| 报错信息 | 可能原因 | 解决办法 |
|---|---|---|
Warning: Solver not applicable | YALMIP认为问题类型与求解器不匹配 | 检查是否引入了非线性项，如两个变量相乘 |
Infeasible problem | 约束过紧或数据矛盾 | 先去掉部分约束测试可行性，逐步加回定位冲突 |
Out of memory | 变量规模过大 | 减少时段数或机组数，检查是否有冗余变量 |
NaN in solution | 数据中有未初始化变量 | 检查数据加载函数，打印关键变量检查是否为NaN |
MIP gap not converging | 求解时间不够或问题规模大 | 设置最大求解时间，适当放松MIP Gap |
Index exceeds matrix dimensions | 循环变量越界 | 检查t+MinUp-1这类索引是否超出T |
最常见的还是Infeasible problem。我排查的经验是：先看功率平衡约束和机组出力上下限是否一致。比如我把机组的总最小出力设置得比系统最低负荷还高，而系统又没有储能和外来输电通道，那这个模型天然就无解。这种问题不是代码写错了，而是数据组合本身不可行。调试方法就是先把所有约束注释掉，一个一个加回来，加上哪个开始无解，问题就在哪里。
6.2 数值尺度问题：一个容易被忽视的大坑
还有一个经验想分享给你，就是数值尺度问题。最开始我的目标函数里，燃料成本动辄几百万，而储热罐的蓄热量只有几百，两者相差好几个数量级。CPLEX这类求解器在内部处理时要靠容差来判断约束是否满足，如果尺度差太远，会出现两个问题：一是收敛速度极慢，二是约束被求解器“忽略”掉。
解决办法有两个。第一个是统一单位，把所有成本都换算成万元，所有功率都换算成百MW，这样量级就接近了。第二个是给YALMIP设置数值容差，但我觉得不如直接改单位来得干净。我当时把成本项除以10000，功率项除以100之后，求解时间直接下降了一半以上，这个优化非常值得做。
6.3 调试技巧：小规模测试先行
最后建议你，拿到一个新模型，别一上来就跑完整24时段、10台机组的规模。先做一个3台机组、6时段的玩具模型，跑通了再逐渐扩大规模。这样做的好处有三个：第一，小规模模型的解你可以手算验证，能确认模型逻辑正确；第二，求解速度快，改一版代码几分钟就能看到结果；第三，小规模时储热罐的运行逻辑更容易可视化观察，方便你理解调度结果。
我自己写代码的习惯是，先把所有约束都放到一个大数组里，然后加一些disp输出每个约束的数量，确认约束数量和你预期是否一致。比如10台机组、24时段，最小启停约束理论上应该有10*(24-MinUp+1)条左右，如果数量对不上，一定是有循环边界写错了。这个检查技巧看起来土，但确实帮我抓到过好几次隐藏bug。
7. 扩展方向与个人心得
代码跑通之后，很多人会问接下来还能做什么。我的建议是往三个方向扩展：
第一，加入不确定性。风电预测是有误差的，把这种误差建模成随机场景，就变成了随机规划或者鲁棒优化问题。调度结果会更保守，但也更贴合实际运行场景。
第二，把储能系统换成电化学储能或者抽水蓄能，模型框架几乎不用动，只需要修改储能设备的约束参数就行。这样一套代码就能覆盖多种储能技术对比，写论文时特别有用。
第三，考虑电网网架约束。现在的模型是单节点模型，没考虑线路潮流约束。如果想研究储热改造对线路阻塞的影响，需要在模型中引入直流潮流方程，那变量和约束都要增加一个维度。
就我个人体会，这个项目对我最大的启发是：储能技术的价值必须在系统层面才能完全体现出来。你要是单独看储热罐本身，它只是一个装热水的罐子；但放到电力系统调度模型里，它就成了连接热力系统和电力系统的关键节点，是火电灵活性的放大器。建模的时候，不要只盯着储能设备本身的公式，要多想想它在整个系统里的角色——是削峰填谷、是促进消纳、还是提供备用，不同定位下模型的约束表达会有微妙差别。
最后再补充一句：Matlab代码写得好不好，很大程度上取决于你对YALMIP的熟练程度。建议花一个下午把YALMIP自带的几个例子跑熟，理解sdpvar、binvar、optimize这三个核心函数的基本用法，后面的项目就顺了。希望这份笔记能帮你少走一些弯路。
