风电最大化消纳：热电联产与储热联合优化控制Matlab实现

我去年接了一个调度优化的活，项目标题很长：风电最大化消纳的热电联产机组联合优化控制（Matlab代码实现）。说白了就是一件事——北方冬天供暖期，热电机组为了供热，夜间电出力压不下去，而风电在凌晨恰恰又是大发时段，电网总体就这么大消纳空间，火电把坑占满了，风电就只能白花花地弃掉。这篇博文想聊清楚：为什么会出现这种矛盾，如何通过热电机组与储热环节的联合优化控制把风电消纳空间腾出来，以及用Matlab从数学建模到代码落地需要走完哪些路。这篇东西适合正在做电力调度、新能源消纳、虚拟电厂方向的同学参考，也适合想用优化方法解决实际工程问题的初学者入门。

1. 问题拆解与总体思路

1.1 弃风矛盾的本质：藏在上游的“以热定电”约束

先说一个物理事实。热电联产机组（CHP）和普通纯凝火电机组不一样，它在发电的同时还要供采暖蒸汽，所以它的运行特性是“热、电耦合”的：为了保证供热，机组必须维持一定的锅炉蒸发量和汽轮机进汽量，而这部分蒸汽在抽汽供热之后，仍然会强制带动发电机转起来。

这就是业内常说的“以热定电”。

用数据感受一下。假设一台300MW的抽凝式机组，热出力与最小电出力的大致关系是：

[ P_{elec_min}(t) = 60 + 0.3 \times H_{chp}(t) \quad (MW) ]

如果夜间的热负荷比较高，比如180MW，那这台机组的最小电出力就被抬到了114MW。也就是说，哪怕电网完全不需要这么多电，它也得至少发这么多。如果系统里有两台这样的机组，总最小出力就是228MW左右（其中一台承担大部分热负荷，另一台在低负荷工况）。此时电网夜间电负荷如果是380MW，风电预测是300MW，火电下限已经占了228MW，留给风电的空间只有152MW，那必然弃掉148MW风，弃风率接近一半。

冬天为什么尤其严重？因为热负荷越高，机组电出力下限被抬得越高；而夜间电负荷又低，风电又大，三个因素叠在一起，系统的调峰能力根本不匹配。这也是“东风夜弃”的根源。

所以，想提高风电利用率，关键不是把风机本身改一改，而是要把火电机组尤其是热电联产机组的“最小出力”降下来。怎么降？最直接的办法不是硬性压机组负荷，那样可能破坏供热安全。

1.2 解锁思路：给热负荷加一个缓冲池

既然“以热定电”是因为“机组必须实时供应热负荷”，那如果我们能在机组和热网之间插一个储热装置——比如大型蓄热水罐，问题就变简单了。

夜间的运行逻辑可以改成这样：

1. 热负荷仍然要满足，但可以不完全由热电机组实时提供；
2. 让机组在白天热负荷高的时候多烧一点，把多余的热量存入储热罐；
3. 夜间热负荷不降的时候，由储热罐放热来替代一部分机组供热；
4. 机组热出力下降了，最小电出力随之下降，电网就腾出了空间接纳风电。

本质上，储热罐是在时间轴上把“热负荷曲线”整体搬了一个位，它是热力和电力之间的缓冲池。白天热负荷高、风电小的时候，机组多出力、多储热，相当于“以空间换时间”；夜间风电大的时候，机组少供热、少发电，风电就能补上来。

这个思路在学界叫“热电解耦”，在工程上叫“机组-储热联合运行”。它不需要新建大型电锅炉，也不需要改造汽轮机，一个蓄热水罐加一套优化调度策略就行，改造量和投资在电厂可接受范围内。这也是它成为主流解耦方案的原因。

所以，联合优化控制的核心不是“怎么调一台机组”，而是把热电机组、储热罐、风电甚至电负荷看成一个整体系统，用优化算法在24小时或96个时段里，统一决定每个时刻“谁出力、谁储热、谁弃风”。

2. 数学模型构建与约束解析

2.1 目标函数怎么设计

做优化第一件事，是回答“优化什么”。在这个场景里，最理想的目标是同时兼顾两件事：一是风电尽可能多上网，二是系统运行成本尽可能低。两者可以写成一个目标函数，用权重来体现优先级。

典型的写法是：

[ \min J = \sum_{t=1}^{T} \left[ \sum_{i} C_i(P_{chp,i}(t)) + \lambda \cdot \sum_{w} (P_{wf_fcst,w}(t) - P_{wind,w}(t)) \right] ]

第一项是热电机组的煤耗成本，(C_i(P_{chp,i})) 通常是一个关于电出力 (P_{chp}) 的二次函数，比如 (aP^2 + bP + c)。第二项是弃风惩罚，(\lambda) 是弃风惩罚系数，单位一般为元/MWh。

这里的关键是 (\lambda) 取多大。如果把弃风惩罚设成800元/MWh，而煤耗边际成本只有200-300元/MWh，优化器就会优先消纳风电，宁可让热电机组深度调峰，也不愿意让风机停机。如果你只是单纯想做“最大化消纳”，也可以直接把目标设为最小化弃风量，但那样可能会得到一个经济上不合理的解——比如为了多消纳10MW风电，让机组频繁大幅度调负荷，对设备寿命影响很大。

实际工程中更稳妥的做法是双目标：先设一个足够大的惩罚系数跑一版解，看弃风率和成本是否在可接受范围；再逐步降低惩罚系数，做灵敏度分析，找到成本和消纳的平衡点。在初版模型里，我习惯把 (\lambda) 设在500~1000元/MWh之间，跑通后再细调。

2.2 约束条件：把物理规则写进模型

模型不能只有目标函数，更核心的是约束条件。这个问题的约束大致可以分成六类，每一类背后都是一条物理或运行规则。

第一，电功率平衡约束。任意时刻，全网发电出力必须等于负荷：

[ \sum_i P_{chp,i}(t) + \sum_w P_{wind,w}(t) + P_{other}(t) = P_{load}(t) ]

这里的 (P_{other}) 代表其他固定出力，比如小水电或联络线功率。严格说还要考虑备用约束，但如果只是做日前计划，先按等式处理就行。

第二，热功率平衡约束。热电机组热出力加上储热罐放热，减去储热罐充热，必须等于热网热负荷：

[ \sum_i H_{chp,i}(t) + H_{sto_dis}(t) - H_{sto_ch}(t) = H_{load}(t) ]

这条约束是“以热定电”的直接数学表达。没有储热时，(H_{sto}) 两项都是0，机组必须实时跟踪热负荷；有了储热，热负荷的刚性就被打破了。

第三，CHP机组的热-电可行域约束。这是决定问题难度的关键。机组运行时，电出力的下限并不是固定值，而是随热出力上升而抬升：

[ P_{chp,i}(t) \ge P_{min0,i} + k_{he,i} \cdot H_{chp,i}(t) ]

同时还有上限和热出力上限：

[ P_{chp,i}(t) \le P_{max,i} ] [ 0 \le H_{chp,i}(t) \le H_{max,i} ]

要注意的是，真实的机组可行域并不是一条简单直线，通常是一个多边形，顶点可能包含背压、抽凝、纯凝等多种工况。工程上常用分段线性化来处理，把多边形拆成多个线性不等式，效果很好。

第四，爬坡约束。机组不能瞬时大幅度改变出力：

[ -\Delta_i \le P_{chp,i}(t+1) - P_{chp,i}(t) \le \Delta_i ]

(\Delta_i) 是机组每时段允许的最大爬坡速率，比如10MW/小时或20MW/小时。有些初版模型会忘掉这条，结果优化器给出的调度曲线在时间轴上剧烈振荡，到了现场根本执行不了。

第五，储热罐模型。首先定义SoC（State of Charge）表示罐内存储的热量，单位MWh。它的递推关系是：

[ SoC(t+1) = SoC(t) + \eta_{ch} \cdot H_{sto_ch}(t) - \frac{H_{sto_dis}(t)}{\eta_{dis}} ]

同时还有容量和功率限制：

[ 0 \le SoC(t) \le C_{sto} ] [ 0 \le H_{sto_ch}(t) \le R_{ch} ] [ 0 \le H_{sto_dis}(t) \le R_{dis} ]

这里 (\eta_{ch}) 和 (\eta_{dis}) 是充放热效率，一般在0.9~0.97之间。很多新手会忽略SoC的初值和终值，导致优化器“白嫖”热量——一开始就利用初始Soc把热负荷满足了，最后罐子被掏空，第二天无法运行。所以必须加约束：

[ SoC(T+1) = SoC(1) ]

这样保证一个调度周期内储热罐净蓄热为零。

第六，风电出力约束。风电场实际并网功率不能超过预测出力：

[ 0 \le P_{wind,w}(t) \le P_{wf_fcst,w}(t) ]

注意风电预测并不总是准确的，但在确定性优化模型里，我们暂时把预测值当作已知参数来用，不考虑预测误差的随机性。这也是这套方法能快速跑通的前提。

2.3 能直接求解吗：线性化与模型选型

上面的模型如果直接丢给求解器，里面涉及的二次成本函数、分段可行域、储热充放热效率，其实都是可以处理的，但具体要看你想用什么求解器。

最简单的版本是线性规划（LP）。把煤耗成本从二次函数分段线性化，机组可行域用多边形线性近似，储热效率用平均效率 (\eta) 统一表示，整个模型就是纯粹LP，Matlab自带的linprog就能求，速度快、调试容易，非常适合第一版验证。

如果想精度更高，把成本保留成二次项，那就是二次规划（QP），用quadprog也能解。Yalmip会自动识别模型类型，不用你自己判断该用哪个求解器。

真正麻烦的是储热同时充放的问题。如果约束里同时存在 (H_{sto_ch}) 和 (H_{sto_dis}) 两个非负变量，理论上优化器可能让它们同时为正——一边充、一边放。虽然这种操作在工程上很傻，但因为目标函数没有直接惩罚它，某些场景下会出现“无用功”。严格的办法是引入0-1变量做互斥约束：

[ H_{sto_ch}(t) \le M \cdot z(t) ] [ H_{sto_dis}(t) \le M \cdot (1 - z(t)) ]

加了二进制变量后，模型从LP变成MILP（混合整数线性规划），需要求解器支持，比如Gurobi、Cplex或者开源工具。不过我在实际项目中的经验是：只要效率 (\eta < 1)，同时充放会导致损耗，目标函数天然会避免它，所以初版不做互斥约束也经常没有异常。先跑简单的，再逐步加精度，是工程上最务实的路径。

3. Matlab实现架构与核心代码逻辑

3.1 环境与求解器选型

代码实现我推荐用Matlab + Yalmip的组合。Yalmip是一个建模层的工具箱，它最大的好处是不用自己拼A矩阵和b矩阵，而是像写数学公式一样去描述变量和约束。你定义一个sdpvar变量，然后用Constraints = [Constraints, ...]往里加约束，最后由Yalmip自动翻译成求解器需要的标准矩阵。

求解器方面，小规模线性模型直接用Matlab内置的linprog或quadprog就够；如果后续加了0-1变量或者规模变大，建议换成Gurobi或Cplex。Yalmip自带接口，只需要在代码里指定'solver','gurobi'即可。我这里演示用的版本是MATLAB R2024a，Yalmip从GitHub下载后放到路径里就能用，操作很简单。实际用R2026b也没有兼容性问题，Yalmip对版本依赖很低。

3.2 数据准备与程序结构

我习惯把所有输入参数放在一个结构体里管理，这样调试多个场景时很方便。比如：

• para.N_chp：热电机组台数
• para.N_wf：风电场个数
• para.T：调度时段数，取24表示1小时一个点，取96表示15分钟一个点
• para.P_load：1xT的电力负荷曲线
• para.H_load：1xT的热力负荷曲线
• para.P_wf_fcst：1xT的风电预测出力曲线
• para.C_sto：储热罐容量（MWh）
• para.R_ch、para.R_dis：最大充热/放热功率（MW）

程序结构建议分成三个文件：main_optimize.m负责加载数据和求解；build_model.m负责构建约束和变量；plot_results.m负责画图和统计。这样后续换场景、调参数，不动模型代码，只改数据文件就行。

3.3 关键代码实现

先定义决策变量。以下是核心片段：

T = 24; N_chp = 2; N_wf = 1; P_chp = sdpvar(N_chp, T, 'full'); % 热电机组电出力，MW H_chp = sdpvar(N_chp, T, 'full'); % 热电机组热出力，MW P_wind = sdpvar(N_wf, T, 'full'); % 风电并网功率，MW H_sto_ch = sdpvar(1, T, 'full'); % 储热罐充热功率，MW H_sto_dis = sdpvar(1, T, 'full'); % 储热罐放热功率，MW SoC = sdpvar(1, T+1, 'full'); % 储热罐蓄热量，MWh

然后构建约束。比较关键的几条是：

Constraints = []; % 电功率平衡 Constraints = [Constraints, sum(P_chp, 1) + sum(P_wind, 1) == para.P_load]; % 热功率平衡 Constraints = [Constraints, sum(H_chp, 1) + H_sto_dis - H_sto_ch == para.H_load]; % CHP机组热-电耦合：电出力下限随热出力上升 for i = 1:N_chp Constraints = [Constraints, P_chp(i,:) >= para.P_min0(i) + para.k_he(i) * H_chp(i,:)]; Constraints = [Constraints, P_chp(i,:) <= para.P_max(i)]; Constraints = [Constraints, 0 <= H_chp(i,:) <= para.H_max(i)]; end % 爬坡约束 for i = 1:N_chp for t = 1:T-1 Constraints = [Constraints, -para.ramp(i) <= P_chp(i,t+1) - P_chp(i,t) <= para.ramp(i)]; end end % 储热罐约束 SoC(1) = para.SoC0; Constraints = [Constraints, SoC(2:end) == SoC(1:end-1) + para.eta_ch * H_sto_ch - H_sto_dis / para.eta_dis]; Constraints = [Constraints, 0 <= SoC <= para.C_sto]; Constraints = [Constraints, 0 <= H_sto_ch <= para.R_ch]; Constraints = [Constraints, 0 <= H_sto_dis <= para.R_dis]; Constraints = [Constraints, SoC(end) == para.SoC0]; % 风电出力约束 Constraints = [Constraints, 0 <= P_wind <= para.P_wf_fcst];

目标函数我习惯先写弃风惩罚，再加煤耗成本：

Objective = 0; % 煤耗成本（简化线性形式，二次形式可自行扩展） for i = 1:N_chp Objective = Objective + sum(para.coal_a(i) * P_chp(i,:).^2 ... + para.coal_b(i) * P_chp(i,:) + para.coal_c(i)); end % 弃风惩罚 Objective = Objective + para.penalty * sum(sum(para.P_wf_fcst - P_wind));

求解设置也很简单：

ops = sdpsettings('solver', 'quadprog', 'verbose', 2); result = optimize(Constraints, Objective, ops); if result.problem == 0 P_chp_opt = value(P_chp); H_chp_opt = value(H_chp); P_wind_opt = value(P_wind); SoC_opt = value(SoC); disp('求解成功'); else disp(['求解失败: ', result.info]); end

跑通一次之后你会发现，Yalmip这套写法最大的优势是：代码的可读性和模型本身的数学表达几乎一一对应，出了问题直接看约束就能定位，而不需要去猜密集矩阵里哪个元素对应哪条物理规则。这一点在工程调试里非常重要。

3.4 结果输出与可视化

求解完之后，至少要画三张图：第一张是电力平衡图，包含电负荷、火电出力、风电并网功率的堆叠曲线；第二张是热力平衡图，展示热负荷由机组和储热罐共同承担的关系；第三张是储热罐SoC曲线，检查它是否在容量边界内、并回到初值。

画图代码不复杂，比如：

figure; stairs(1:T, para.P_load, 'k-', 'LineWidth', 1.5); hold on; stairs(1:T, sum(value(P_chp),1), 'r-', 'LineWidth', 1.2); stairs(1:T, value(P_wind), 'b--', 'LineWidth', 1.2); legend('电负荷', 'CHP电出力', '风电并网'); xlabel('小时'); ylabel('功率/MW'); grid on;

看曲线的时候最该关注三点：一是风电是否被压缩，也就是蓝线和负荷线的间隙是否过大；二是机组出力有没有频繁来回振荡，如果有，多半是爬坡约束没落地或惩罚系数设置不当导致的；三是储热罐SoC是否在合理区间，如果整天贴着上限或下限跑，说明储热容量选得不合适。

4. 典型场景仿真与关键参数调试

4.1 有储热与无储热的对比

为了验证联合优化到底有多大作用，我习惯做一组“基线对照实验”：场景A是没有储热罐，热电机组严格“以热定电”；场景B加上储热罐，做联合优化。其他数据全部保持一致。

下面是一组示例数据，不代表真实系统，但趋势和量级是典型的：

可以看到，储热罐加入之后，弃风率直接下降了一半左右。原因是夜间储热放热替代了部分机组供热，机组最小电出力被压低，风电并网空间变大。同时煤耗成本反而略微下降，因为白天热负荷高的时候机组带高负荷更接近经济工况，夜间压低负荷的损失被储热的转移效应补偿了。

当然这只是示例数值，真实收益取决于热负荷峰谷差、机组热-电耦合系数、储热容量和风资源波动情况。如果热负荷曲线很平缓，储热能发挥的空间就小；如果热负荷峰谷差大，储热的效果就显著。这也是为什么做项目时一定要拿实际数据跑，不能拍脑袋。

4.2 储热容量与惩罚系数的灵敏度

储热容量不是越多越好。我有一次做敏感性分析，把储热罐容量从0一直扫到400MWh，发现弃风率下降的边际效果呈明显的递减趋势。容量从100加到200MWh时，弃风率下降了4个百分点；但从300加到400MWh时，只下降了不到1个百分点。原因是储热罐容量一旦超过了夜间热负荷转移总量的需求，再大的罐子也没有多的热可以装。

所以选储热容量时不能只看罐子本身价格，要结合热负荷曲线、机组可调节范围、风电预测数据一起算。用一个简单方法：计算夜间风电大发时段的理论最大蓄热需求量，再乘以1.1~1.2的安全裕度，基本就是合理容量区间的上限。

惩罚系数 (\lambda) 的调整也讲究。我的经验是先把 (\lambda) 设得非常大（比如10000元/MWh），跑出“极限消纳解”，看成本和弃风率；再把 (\lambda) 降下来，比如设成300元/MWh，看系统在“经济优先”模式下怎么取舍。两个极端都跑一遍，中间那条权衡曲线就出来了。实际调度中，很多地方并不要求100%消纳，而是允许一个很小比例的弃风来换取火电运行的稳定性，这时候惩罚系数就相当于给调度员提供了一个“消纳优先级旋钮”。

5. 常见问题与排查技巧实录

5.1 求解不可行的排查路径

新手遇到infeasible problem的第一个反应往往是怀疑求解器坏了。其实绝大多数情况是约束之间自相矛盾。

我自己的排查顺序大致是：

1. 先检查电平衡和热平衡的等式是否严格成立。有时候输入曲线的采样时间不一致，比如电力负荷是15分钟一个点，热负荷是1小时一个点，直接相减就会导致数学上无解。
2. 再检查储热罐约束。如果SoC的初值和终值相等约束加上后，和充放热极限功率参数冲突，比如一个时段最多只能充50MW，但要求一天内把200MWh热量全部转存，那怎么都不可能满足。
3. 然后检查机组可行域约束。热-电耦合关系的斜率 (k_{he}) 如果设得太大，会让机组电出力下限被抬到超过上限，这属于参数错误。

还有一个很实用的调试技巧：在所有等式约束上临时加一个松弛变量，比如电平衡改成：

[ \sum P + P_{other} + P_{slack} = P_{load} ]

然后看最优解中 (P_{slack}) 在哪个时段非零，那个时段就是约束冲突最严重的地方。定位到具体时间段后，再回去检查对应时段的负荷和机组状态，问题基本一眼就能看出来。

5.2 数值尺度与求解器设置的坑

优化求解对数值尺度非常敏感。如果热负荷的单位是GJ/h、电出力的单位是MW、储热罐容量单位是kWh，三个量级差到百万倍，求解器内部的容差设置很容易出问题，最终结果出现“看似可行，实则错误”的情况。

我的建议是统一所有量纲：功率一律用MW，能量一律用MWh，成本一律用元或万元。如果数据原始单位是kW或GJ，先除以1000换算一下，再进模型。代码里可以在读取数据时就完成换算，并加一行注释说明，免得隔几天回来忘了单位。

另外，如果在Yalmip里选linprog却遇到检测不到求解器，通常是Yalmip把模型识别成了LP但系统路径里没有linprog的写权限。解决方法是：

ops = sdpsettings('solver', 'linprog', 'verbose', 2);

显式指定求解器，可以跳过Yalmip的自动探测。用Gurobi也有类似问题，安装后需要运行gurobi_setup或在sdpsettings里指定路径。

5.3 建模时容易忽略的细节

还有几个细节值得单独拿出来提醒。

第一，爬坡约束的维度方向。用diff(P_chp, 1, 2)时会得到T-1个差分值，如果直接和长度为T的向量比较，Yalmip会报维度不匹配。建议用双层for循环逐台机组、逐个时段构建，虽然代码长几行，但绝对不会因为维度问题翻车。

第二，储热罐的充放热效率别写反。SoC递推公式里，充热要乘效率，放热要除效率。如果写反了，优化器会发现“充1MW热量进去变成0.95MW，放出来却变成1.05MW”，这相当于无中生有，结果就是储热罐能量不守恒，解会偏离物理事实。

第三，不要急于引入非凸约束。第一版模型用线性约束跑通后，先验证趋势和结果是否符合物理直觉，再逐步加入更精细的非线性描述。我接手过不少项目，一上来就上MILP加非线性成本函数，结果模型又大又难调，最后反而是简化版本先解决了业务问题。工程上，简单模型跑出来的粗糙结果，远比复杂模型跑不出来的精确结果有价值。

我在实际项目中最大的体会是：风电最大化消纳这件事，真正的难点不在算法多高深，而在于把热和电两条能量流之间的矛盾用数学模型精确表达出来。一旦你理解了“以热定电”的机理，知道储热罐是解开这个死结的钥匙，后续的Matlab实现、代码调试都只是水到渠成的工作。做完这套优化控制之后，我再看到弃风率报表时，第一反应不再是感慨风机不行，而是去翻热负荷曲线和储热罐状态——因为问题的答案，很可能就藏在热网那边。