热电联产机组联合优化调度：Matlab+YALMIP建模风电消纳与储热电锅炉算例-开发者社区

1. 冬季供暖季的弃风困局：热电联产机组到底卡在哪

每年供暖季一过，风电场的同事就开始盯着调度曲线叹气：白天风光还好，一到后半夜风速上来了，风电场却得压出力，甚至有整场停机的时候。而另一边，热电厂里锅炉烧得正旺，发电量被热负荷“绑架”着往下压也压不下来。这种供暖季弃风集中的现象，根源就在热电联产（CHP）机组的运行方式上。

所谓“以热定电”，意思是供热机组的电出力下限由热负荷决定——只要热网要那么多热量，机组就得维持对应的最小电出力。热负荷在冬天是刚需，尤其夜间居民采暖需求并不低，于是凌晨低谷时段电网调峰空间被热电机组堵死，风电再便宜也送不进来。这个问题在三北地区特别典型，也是新能源消纳研究里绕不开的老话题。

这篇内容把热电联产机组联合优化控制从建模到Matlab实现完整过了一遍：先讲清楚卡点是怎么形成的，再给出可求解的数学模型，然后是用YALMIP搭求解程序的完整流程，最后放一个含储热罐和电锅炉的算例对比，并附上我自己调试过程中踩过的坑。适合正在做新能源消纳、热电协同调度、灵活性改造方向研究或毕设的同学参考，代码思路可以直接复用到你自己的场景里。

1.1 弃风为什么集中出现在后半夜

很多人以为弃风是风资源预测不准的问题，实际上大部分弃风发生在电网“消纳能力”不足的时候。冬季后半夜有两个特点：一是负荷处于全天低谷，二是风电往往因为夜间风速大而出力高。两个因素叠加，系统需要压掉大量常规机组来给风电让路。

但热电联产机组让不开。供热管网有严格的热水温度要求，机组的热出力必须满足热负荷，而热出力一旦定下来，电出力就只能在很小的范围内调整。抽凝式机组还好一点，热、电区间是一个四边形区域，还能往下压一部分；背压式机组干脆就是热电线性耦合，发多少电由供多少热决定，基本没有调节空间。一个拥有多台背压机组的工业园区热电厂，在供暖期的调峰能力几乎为零，对电网来说就是一块铁板。

1.2 电网侧视角：调峰空间的挤压过程

把系统简化来看，某时刻电网的电功率平衡可以写成：

ΣP_CHP + P_wind + P_other = P_load

其中P_load是系统负荷，P_wind是实际接纳的风电功率。要最大化风电消纳，本质上是让常规机组（尤其是热电机组）尽量压低出力。但热电机组存在电出力下限P_MIN，这个下限随着热出力H升高而升高。于是风电实际可消纳的空间变成：

P_wind_available = P_load - P_other - ΣP_CHP_min(H)

热负荷越高，ΣP_CHP_min(H)越大，留给风电的空间就越小。当P_wind_available小于风电场预测出力时，多余的功率只能弃掉。这就是“以热定电”压掉风电空间的完整逻辑链路。明白了这一层，后面所有优化手段——加储热、加电锅炉、机组灵活性改造——本质上都是在想办法把这个P_CHP_min(H)往下拉，或者把热负荷与电出力之间的硬绑定关系解开。

2. 优化控制建模：把“以热定电”写成可求解的数学问题

搞清楚了物理机理，下一步就是把问题翻译成数学模型。建模的思路是：给定系统负荷、热负荷和风电预测曲线，在满足机组运行约束、管网热平衡约束和储热设备约束的前提下，决定各台机组每一时刻的出力，让目标函数最优。

这里的“联合优化”体现在两个维度：一是多台CHP机组之间的出力分配，二是CHP机组与储热罐、电锅炉、风电场之间的协同。不是单台机组自己调，而是整个供热-供电系统一起算。

2.1 目标函数怎么定：消纳最大化与经济性的博弈

目标函数有两种常见写法。第一种是纯物理量目标，直接最大化研究周期内的风电消纳量：

max Σ P_wind_t * Δt

这个写法简单直观，适用于学术研究和机理分析，但有个问题：完全不考虑机组运行成本，算出来的结果可能让某台机组长时间在很恶劣的工况点运行，经济上不现实。

第二种是经济性目标，把弃风作为惩罚项放进去：

min Σ (机组燃料成本 + 弃风惩罚成本 + 储热设备运行成本)

这里的弃风惩罚系数很关键，它代表弃风的社会成本或折算电价。系数设置多少没有统一标准，我一般先按当地风电上网电价和火电标杆电价的差值来估算，再在算例里做敏感性分析。目标函数的选择会直接影响求解结果，后面第5节会专门讲权重调参的教训。

2.2 机组运行可行域：抽凝式与背压式的关键约束

热电联产机组最核心的约束就是运行可行域。抽凝式机组的可行域通常用一个四边形（或者多边形）来描述，简化线性化后可以用以下不等式组表示：

P_min + c_v * H ≤ P ≤ P_max - c_d * H

其中P是电出力，H是热出力，c_v是低压缸最小凝汽流量对应的热电比斜率，c_d是最大进汽量下的热电比斜率。直观理解就是：热出力越大，电出力可调范围越窄，整体上移。

背压式机组更简单也更严格，电出力与热出力满足线性关系：

P = c_m * H

一台背压机组基本上就是一个“热电转换器”，没有独立调节能力。

在Matlab里，这些约束都是线性的，直接写成不等式矩阵交给求解器即可。需要注意P和H的量纲统一，我在第一次建模时吃过这个亏，后面会细说。

2.3 储热罐与电锅炉的模型化

储热罐的作用是“搬移”热量：在风电大发、热电机组需要压电出力的时段，从热网取热储起来；在热负荷高峰或风电低谷时，再放热补充供热，从而降低机组当前时刻必须提供的热出力，间接压低了电出力下限。

储热罐的状态方程：

S(t+1) = S(t) + (Q_ch(t) * η_ch - Q_dis(t) / η_dis) * Δt - S_loss(t)

约束包括：

0 ≤ S(t) ≤ S_max

0 ≤ Q_ch(t) ≤ Q_ch_max

0 ≤ Q_dis(t) ≤ Q_dis_max

S(0) = S(T)

最后一条是周期约束，保证储热罐在一整天调度结束后恢复到初始状态，否则昨天用掉的能量没有结算边界，长期下来算不出稳定运行方案。

电锅炉则是把电能转成热能：

Q_eb(t) = η_eb * P_eb(t)

电锅炉加入系统后，风电大发时段可以让电锅炉直接消耗多余风电制热，相当于给风电增加了“就地消纳”的出口。它和储热罐的配合逻辑是：电锅炉造热、储热罐存热、热电机组降出力，三者同步动作，把弃风空间腾出来。

2.4 功率平衡与热平衡约束

系统的电功率平衡和热功率平衡必须逐时段满足。电平衡：

Σ P_chp_i(t) + P_wind(t) - P_cur(t) + P_other(t) = P_load(t)

这里P_cur(t)是弃风功率，它在目标函数里被惩罚，所以优化会自动压低它。

热平衡：

Σ H_chp_i(t) + Q_dis(t) + Q_eb(t) = H_load(t) + Q_ch(t)

注意放热和电锅炉产热是热源，蓄热是热负荷，方向别搞反了。热网本身还有热损耗和传输延迟，严格模型里会用热网动态方程描述，但做中长期调度时通常简化为准稳态模型，损耗按比例系数处理即可。

3. Matlab代码落地：从公式到可跑通的求解程序

模型建完了，接下来是最耗时间的部分：把数学模型写成能跑的Matlab代码。这里我用的是YALMIP工具箱加外部求解器的方案，这是目前做这类优化调度比较顺手的一套组合。

3.1 工具链选型：为什么用YALMIP而不是手写求解代码

我见过不少同学自己写拉格朗日松弛或者内点法去解调度问题，精神可嘉但效率太低。这类热电联合优化本质上是一个线性规划（LP）或混合整数线性规划（MILP）问题，直接用成熟求解器又快又稳。YALMIP提供了一个建模层，让你用接近数学公式的语法描述优化问题，再由求解器处理计算，改约束和调权重都非常方便。

% 安装YALMIP后，只需要把路径加入Matlab即可 addpath(genpath('D:\tools\yalmip')); % 求解器推荐：CPLEX 或 Gurobi，学术许可可以申请

我在实际项目中用的是Gurobi 9.x版本，配合Matlab 2020a，稳定性不错。如果实验室没有商业求解器，Matlab自带的linprog和intlinprog也能解，只是大规模场景下速度会慢一些。

3.2 数据结构与参数准备

建模前先把所有系统参数整理成结构体，建议按“机组-热源-风场-负荷”四类分开管理。

%% 系统基础参数 T = 24; % 调度周期24小时，步长1h dt = 1; % 时间步长 load_elec = [560 545 ...]; % 电负荷曲线，长度24 load_heat = [420 405 ...]; % 热负荷曲线，长度24 wind_forecast = [180 200 ...]; % 风电预测出力，长度24 %% 热电机组参数（两台抽凝式） chp.chp_num = 2; chp.P_max = [180 180]; % 最大电出力 MW chp.P_min = [60 60]; % 最小电出力 MW chp.c_v = [0.55 0.55]; % 热电耦合系数（最小凝汽工况） chp.H_max = [220 220]; % 最大热出力 MW %% 储热罐参数 ts.S_max = 300; % 储热容量 MWh ts.Q_ch_max = 80; % 最大蓄热功率 MW ts.Q_dis_max = 80; % 最大放热功率 MW ts.eta_ch = 0.95; % 蓄热效率 ts.eta_dis = 0.95; % 放热效率 ts.S_init = 150; % 初始储热量 %% 电锅炉参数 eb.P_max = 50; % 最大功率 MW eb.eta = 0.98; % 电热转换效率

这里P_min其实是“纯凝工况下最小电出力”，有了机组供热后，实际电出力下限是P_min + c_v * H，这才是关键约束。

3.3 约束条件逐条写成代码

有了参数，约束条件的代码化就相对机械了。核心思路是用sdpvar声明决策变量，用竖线|累加约束，最后交给optimize求解。

% 声明决策变量 P_chp = sdpvar(T, chp.chp_num, 'full'); % 机组电出力 H_chp = sdpvar(T, chp.chp_num, 'full'); % 机组热出力 P_eb = sdpvar(T, 1, 'full'); % 电锅炉耗电功率 S_soc = sdpvar(T+1, 1, 'full'); % 储热罐状态 Q_ch = sdpvar(T, 1, 'full'); % 蓄热功率 Q_dis = sdpvar(T, 1, 'full'); % 放热功率 P_wind = sdpvar(T, 1, 'full'); % 实际风电消纳 P_cur = sdpvar(T, 1, 'full'); % 弃风功率

然后把每一类约束写清楚。机组可行域是最容易写错的地方，先展示抽凝式的写法：

Constraints = []; for t = 1:T for i = 1:chp.chp_num Constraints = [Constraints, ... P_chp(t,i) >= chp.P_min(i) + chp.c_v(i) * H_chp(t,i)]; Constraints = [Constraints, ... P_chp(t,i) <= chp.P_max(i)]; Constraints = [Constraints, ... H_chp(t,i) >= 0, H_chp(t,i) <= chp.H_max(i)]; end end

注意这里只写了左边界和下边界，实际上完整的抽凝机组可行域还应该包含右边界斜线（P ≤ P_max - c_d * H）。如果机组运行在部分负荷工况，右边界也需要带上，否则求解器可能给出“看起来合理但物理上达不到”的出力组合。

储热罐的状态转移和容量约束：

for t = 1:T Constraints = [Constraints, ... S_soc(t+1) == S_soc(t) + (Q_ch(t)*ts.eta_ch - Q_dis(t)/ts.eta_dis)*dt]; Constraints = [Constraints, ... 0 <= S_soc(t+1) <= ts.S_max]; Constraints = [Constraints, ... 0 <= Q_ch(t) <= ts.Q_ch_max]; Constraints = [Constraints, ... 0 <= Q_dis(t) <= ts.Q_dis_max]; end Constraints = [Constraints, S_soc(1) == ts.S_init]; Constraints = [Constraints, S_soc(T+1) == ts.S_init]; % 周期约束

电锅炉模型和功率平衡：

for t = 1:T % 电锅炉电热转换 Constraints = [Constraints, Q_eb(t) == eb.eta * P_eb(t)]; Constraints = [Constraints, 0 <= P_eb(t) <= eb.P_max]; % 电功率平衡：机组出力 + 风电 + 其他 = 负荷 + 电锅炉耗电 Constraints = [Constraints, ... sum(P_chp(t,:)) + P_wind(t) == load_elec(t) + P_eb(t)]; % 风电出力约束和弃风定义 Constraints = [Constraints, ... P_wind(t) + P_cur(t) == wind_forecast(t)]; Constraints = [Constraints, ... 0 <= P_wind(t) <= wind_forecast(t), P_cur(t) >= 0]; % 热功率平衡 Constraints = [Constraints, ... sum(H_chp(t,:)) + Q_dis(t) + Q_eb(t) == load_heat(t) + Q_ch(t)]; end

3.4 求解设置与结果可视化

目标函数选经济型写法，弃风惩罚成本设一个足够大的系数，让求解器优先消灭弃风：

% 弃风惩罚成本系数（元/MWh） penalty_cur = 800; % 机组煤耗成本近似用线性函数 cost_a = [0.28 0.28]; % 煤耗系数 吨/MWh cost_b = [10 10]; % 固定成本 Objective = 0; for t = 1:T for i = 1:chp.chp_num Objective = Objective + cost_a(i) * P_chp(t,i) * dt; end Objective = Objective + penalty_cur * P_cur(t) * dt; end ops = sdpsettings('solver', 'gurobi', 'verbose', 1); sol = optimize(Constraints, Objective, ops);

求解完以后，检查sol.info是否是成功标志，然后把结果画出来。我习惯画三张图：一是电功率平衡堆叠图（机组出力、风电、电锅炉耗电），二是热功率平衡图（机组供热、储热充放、电锅炉产热），三是储热罐SOC曲线和风电消纳率对比。这三张图基本可以把调度逻辑讲清楚，写报告和汇报都够用。

4. 算例对比：储热罐和电锅炉到底能多消纳多少风电

代码跑通了，接下来应该回答一个实际工程问题：加了储热罐和电锅炉之后，风电消纳率到底提高了多少？我用一个简化算例做了三组方案对比。

4.1 算例场景设置

算例假设一个区域系统包含两台抽凝式热电联产机组（参数见3.2节）、一座200 MW风电场、一个储热罐和一个50 MW电锅炉。冬季典型日负荷曲线取晚间较高、凌晨较低，风电预测为凌晨大风、白天平风。设计了三组方案对照：

方案A：热电联产机组独立运行，无储热、无电锅炉（基准场景）
方案B：加装储热罐，容量300 MWh，充放功率80 MW
方案C：储热罐+电锅炉联合运行

三组方案都用同一个目标函数和同一组负荷、风电数据，只改变可用设备集合。

4.2 典型日结果分析

先看方案A。凌晨1点到5点，热负荷约380 MW，两台机组按“以热定电”方式最低需要出力约120 MW左右，加上系统内其他常规电源，电网留给风电的空间只有90 MW上下，而风电预测出力接近180 MW，剩余全部弃掉。全天弃风率约28%，集中在凌晨时段。

方案B加入储热罐后，调度策略明显改变了：凌晨时段机组虽然仍然要供热，但一部分热负荷由储热罐放热承担，机组热出力可以降低，对应的电出力下限也跟着降了10到15 MW，腾出来的空间全部给了风电。储热罐SOC呈现出“夜间放热、白天蓄热”的典型模式。全天弃风率降到13%。

方案C在储热罐基础上加入电锅炉，效果又进了一步。凌晨风电大发时，电锅炉直接消耗约35 MW风电功率制热，这部分热量存入储热罐或者直接补入热网，机组供热压力进一步减小。全天弃风率压到4.5%，基本实现了风电的最大化消纳。

4.3 从结果反推：关键瓶颈在哪个约束

这个算例最值得玩味的是方案B和方案C的差异。方案B已经用储热罐搬移了大量热负荷，为什么弃风率还有13%？回头看约束条件：电锅炉没有启用，凌晨时段依赖储热罐放热降低机组热出力，但储热罐容量有限，300 MWh的储热量到凌晨4点就快放空了，后续时段机组被迫提高出力，风电空间再次被压缩。所以方案B的瓶颈是储热容量约束，而不是机组技术约束。

方案C里电锅炉相当于一个“可调热负荷”，直接把多余风电变成热能，相当于把终端用能结构从“必须烧煤产热”变成“风电多了就用电产热”，从源头上把多余风电消化掉。这个对比给工程实践的启示很直接：储热罐擅长“腾挪”热量改善时间分布，电锅炉擅长“吸收”瞬时多余风电，两者配合才是比较完整的解耦方案，只上其中一个都难以把消纳率做到极致。

5. 调试路上的坑：数值问题、求解器与参数敏感性

代码从能跑到结果可信，中间隔着不少坑。以下几条都是我自己实际踩过的，写出来供参考。

5.1 YALMIP安装与求解器对接的三个小坑

第一，YALMIP版本和Matlab版本不兼容会导致一堆奇怪的报错，最典型的是调用sdpvar时提示类未定义。建议从GitHub拉最新master分支，别用网上流传的老压缩包。第二，Gurobi和YALMIP对接后第一次optimize可能提示找不到求解器，原因往往是许可证环境变量没设置，检查GUROBI_HOME是否指向了正确目录。第三，sdpsettings里'solver','gurobi'的写法必须准确，大小写不对会静默退化到默认求解器，算出来的结果可能不对。我的习惯是先跑一个已知最优解的小例子验证工具链，再上正式模型。

5.2 量纲不一致导致的数值震荡

第一次完整建模时，我用了三个不同量纲的数据：电负荷单位取MW，热负荷单位取GJ/h，储热罐容量单位取MWh。单看每个约束都没问题，但放到同一个目标函数里加权求和时，量纲差异导致数值范围跨了好几个数量级，Gurobi迭代时一直报数值警告，结果也不稳定。后来统一把所有功率都用MW、能量用MWh、时间用小时，目标函数里各项的量纲就齐了。这个坑很基础，但很容易被忽略。

5.3 储热罐SOC初值：周期约束不写满的教训

有一版代码我图省事，只给了S(1)初值，没写S(T+1) == S(1)的周期约束。结果求出来的“最优解”里，储热罐从初始150 MWh一路蓄到接近满罐，全天结束时剩280 MWh。从数学模型看这个解没毛病，目标函数确实更小了，但物理上站不住脚——第一天结束时白赚了130 MWh的热量，第二天就没有这个“免费的午餐”了。加上周期约束后，储热罐才真正进入“一天一个循环”的合理运行状态。这个教训也适用于其它带储能设备调度的模型。

5.4 参数敏感性：弃风惩罚系数为什么不能乱定

惩罚系数penalty_cur从200调到800，系统的运行策略会有质的变化。系数偏低时，求解器优先保证机组运行经济性，弃一点风也无所谓；系数偏高时，求解器会不惜让机组在接近下限的低效工况运行来换风电。实际项目里，这个系数的合理取值应该反映真实的弃风补偿成本，比如我按风电上网电价0.35元/kWh、火电煤耗成本折算后，取了一个略高于两者差值的数据，结果是机组和风电各有让步，整体运行成本最低。而不是把系数拍脑袋设成一个很大的数——那样虽然“消纳率最美观”，但总成本可能反而更高。