做电力系统优化调度，尤其是冬季供暖期的运行分析时，八成以上的人都会碰到同一个痛点：热负荷一上来，热电联产机组就被“以热定电”摁得死死的，电网的调峰空间被压得很窄，风电场那边却经常是大风满发。这时候如果你手头只有常规火电的调度模型，算出来的弃风量往往会高得离谱，但又说不清楚问题到底出在哪。我去年做某区域电网冬季调峰方案时，就把抽凝式热电联产机组、背压式热电联产机组、蓄热罐和电锅炉放到同一个优化框架里，用Matlab结合Yalmip和Cplex求解器建模仿真，把风电消纳率从85%左右提到了96%以上。这篇文章就把这个项目的建模思路、约束怎么写、Matlab代码怎么组织、以及我踩过的坑完整梳理一遍。
这个题目看着是“热电联产机组联合优化控制”，本质上是一个带风电接入的多能源系统优化调度问题。适合两类人参考：一是做电力系统方向的研究生，需要复现论文里的算例；二是刚接触调度建模的工程师，想用代码快速验证蓄热罐和电锅炉对风电消纳的作用。下面我会把模型、代码、结果分析和排查经验分开讲。
1. 问题背景：为什么热电联产会和风电“抢空间”
1.1 “以热定电”是怎么束缚住调峰能力的
热电联产机组的特点是同时生产电和热。冬季供暖期，热负荷由热网需求决定，CHP机组为了满足热负荷，电出力不能随意降低。这就是行业里常说的“以热定电”。反映在数学模型上，抽凝式机组有一个电热耦合特性：电出力下限会随着热出力上升而抬高。简单写成约束就是：
P_chp >= P_min + cv * H_chp
其中 cv 是电热耦合系数。热出力 H_chp 越大，电出力下限 P_chp 就越高。比如一台额定最大电出力200MW的抽凝机组，纯凝模式下最低电出力是80MW，但如果带150MW的热负荷，电出力下限可能被抬到110MW甚至更高。
这里就会有一个明显的矛盾：深夜时段电负荷低、热负荷却不低（供暖需要），CHP机组为了供热必须维持较高电出力，风电却往往这时段风速大好发。系统里火电压不下去、风电上不来，弃风就这么产生了。可以说，热负荷一天不降，这个调峰压力就一天不消失。
1.2 解耦的两个经典方向：蓄热罐与电锅炉
既然矛盾根源是电和热被“绑定”，那就想办法解耦。工程上最常用的两个手段，一个是蓄热罐，一个是电锅炉。
蓄热罐的思路是“错峰供热”。风电大发、电负荷低谷时，让CHP机组降低电出力、少产热，缺的那部分热量由蓄热罐放热补上；等到电负荷升高、风电变小，再让CHP机组加大电出力、多产热，把多余热量存进蓄热罐。这样热负荷在时间维度上被“搬移”了，CHP电出力也因此获得了下调空间。
电锅炉的思路更直接，就是把风电转换成热。风电大发时，电锅炉投入运行，电功率变成热功率直接供给热负荷，替代一部分CHP的供热出力，CHP电出力也可以跟着降。它的优点是响应快、控制简单，缺点是电能转换为热能存在一定损耗，而且电锅炉本身并不储存能量。
这两个设备不是互斥关系。实际项目中我通常把它们一起配置：蓄热罐负责日内大范围的热量搬移，电锅炉负责短时深度压减CHP电出力。两者配合，风电消纳效果明显好过只上一种设备。
1.3 多机组联合优化的本质是“在约束里求最优”
你可能会问：为什么要“联合优化”？每个机组各调各的不行吗？还真不行。热电厂内一般有多台CHP机组，它们的电热耦合特性、容量、煤耗都不一样；再加上蓄热罐、电锅炉这些灵活性资源，系统的可调变量很多。如果靠人工经验去分配出力，不仅效率低，而且难以保证全局最优。
联合优化的本质，就是在一个统一的目标函数下，同时决策所有机组和设备的运行状态，满足电力平衡、热力平衡、设备出力上下限、爬坡速率、储热罐能量状态等约束条件。Matlab的优势在于，可以用Yalmip这类建模工具把数学模型直接“翻译”成代码，再交给Cplex、Gurobi这类商用求解器去算。你不需要手写求解算法，只需要把模型表达准确。
2. 数学模型设计：目标函数与约束条件怎么定
2.1 目标函数：运行成本里藏着弃风惩罚
优化调度模型的“指挥棒”是目标函数。很多论文里写的是“系统运行成本最小”，但其实要让风电最大化消纳，必须在目标里给弃风一个足够高的惩罚。否则求解器只盯着成本，该弃风时它就弃风了。
我采用的优化目标包含三部分：
- CHP机组的煤耗成本，近似为电出力和热出力的线性函数；
- 电锅炉的运行维护成本，很小，只起防止电锅炉乱开的作用；
- 弃风惩罚成本，单位惩罚系数设得足够大，保证优先消纳风电。
目标函数表达式可以写成：
min F = sum_t sum_i ( a_i * P_chp(i,t) + b_i * H_chp(i,t) ) + sum_t ( c_eb * P_eb(t) ) + sum_t ( w_wind * ( P_w_forecast(t) - P_w(t) ) )
其中 w_wind 是弃风惩罚系数，我取5000元/兆瓦时，远远高于煤耗成本，相当于给“少弃风”一个更高的优先级。电锅炉成本系数 c_eb 取50元/兆瓦时，主要是防止它无意义运行。
这里要注意，弃风量是用“预测可发功率减去实际上网功率”来刻画的。如果你直接用“风电出力变量”，模型会有投机空间，可能在不需要弃风时也压低风电出力。用预测值做差值之后，风电出力少的代价就是弃风惩罚，模型自然会尽量让风电多上。
2.2 电热耦合约束：抽凝机组可行域的线性化处理
CHP机组的电热耦合是模型里最容易写错的部分。抽凝式机组不是简单的一条线性关系，而是有一个可行域多边形。但为了工程计算方便，我做了简化处理，只保留最关键的下限耦合约束：
P_chp(i,t) >= P_chp_min(i) + cv(i) * H_chp(i,t)
同时热出力有上下限：
H_chp_min(i) <= H_chp(i,t) <= H_chp_max(i)
电出力也有上下限：
P_chp_min(i) <= P_chp(i,t) <= P_chp_max(i)
如果你查阅完整版学术模型，还会看到电出力上限随热出力下降的耦合约束、抽汽进汽量约束等。实际代码里把这些约束补上不难，但核心机制就是下限耦合这一条。我建议初次复现时，先跑通这个简化版本，确认结果合理，再逐步加严约束。
爬坡约束也得写进去：
-ramp_down(i) <= P_chp(i,t) - P_chp(i,t-1) <= ramp_up(i)
CHP机组的爬坡速率一般取额定容量的10%~20%，我取的是20MW/小时和15MW/小时。热出力变化速率在很多文献里默认足够快，可以不约束，但如果你研究的是背压式机组，热出力爬坡也需要纳入考虑。
2.3 蓄热罐与电锅炉约束：不能拍脑袋写
蓄热罐的核心约束是能量状态递推关系。我定义 H_ch(t) 是充热功率，H_dis(t) 是放热功率，S(t) 是储热量，dt取1小时，那么：
S(t+1) = S(t) + ( eta_ch * H_ch(t) - H_dis(t) / eta_dis ) * dt
充放热效率通常取0.9~0.95。蓄热罐容量、充放热功率上限、储热量上下限都要给约束。还有一个容易被忽略的约束是“周期始末储能相等”，也就是一天下来蓄热罐的净储热量不能变化太大，否则相当于从“地下凭空采热”，结果会失真。我一般写成 S(T) >= S(0) ，要求一天结束时的储热量不低于初始值。
比较麻烦的是充放热不能同时进行。如果模型里允许 H_ch 和 H_dis 同时大于0，最极端的情况下，求解器可能让蓄热罐一边充热一边放热，能量在罐内“空转”，虽然不违反约束，却造成资源浪费。解决办法是引入二进制变量把两者互斥：
0 <= H_ch(t) <= M * u(t) 0 <= H_dis(t) <= M * (1 - u(t))
这里的 u(t) 是0-1变量，u(t)=1表示充热，u(t)=0表示放热。引入二进制变量后，模型从线性规划LP变成混合整数线性规划MILP，求解时间会变长，但换来的是结果符合物理直觉。
电锅炉的约束就简单得多：
0 <= P_eb(t) <= P_eb_max
它消耗的电功率转化为热功率时，存在一个效率系数 eta_eb，我取0.95。转化关系体现在热功率平衡约束里：
eta_eb * P_eb(t) 进入热负荷侧。
电锅炉容量怎么选？我一般按风电预测出力的20%~30%配置。容量太小，压减CHP出力幅度有限；容量太大，投资回报不划算，而且电网也未必有那么多弃风电量给它用。
2.4 整体模型的性质与求解思路
把上述约束和目标函数组合起来，这是一个混合整数线性规划问题。决策变量包括：
- 每台CHP机组每个时段的电出力、热出力；
- 蓄热罐每个时段的充热/放热功率和储热量；
- 电锅炉每个时段的电功率；
- 风电场每个时段实际上网功率；
- 蓄热罐充放热互斥的二进制变量。
整体模型规模不大，24小时时段、2台机组、100多个变量，Cplex基本上秒解。就算你把机组数量扩展到10台、时段扩展到168小时（一周），求解速度也完全可以接受。我在实际项目中用Cplex求解168小时算例，时间一般控制在10秒以内。
3. Matlab实现要点：从数据到可运行代码
3.1 环境准备：Yalmip + 求解器的安装配置
这个项目我用的是Matlab R2022b、Yalmip工具箱和Cplex 12.10求解器。Yalmip是瑞典学者Lofberg开发的一个免费建模工具，它本身不算求解器，而是把优化模型“翻译”成求解器能认识的标准形式。比如你写了一个带二进制变量的优化问题，Yalmip会帮你自动转换成MILP格式，再交给Cplex去解。
安装上踩过的坑比较多，我简单说下要点。Yalmip不是MathWorks官方工具箱，需要去GitHub或者官网下载压缩包，解压后放到任意目录，然后在Matlab里用addpath(genpath('Yalmip目录'))加入路径。Cplex则需要先安装IBM官网的完整版，安装完成后，在你写的代码开头多写一行addpath('Cplex安装目录\cplex\matlab\x64_win64')。很多新手报“solver not found”错误，99%都是求解器路径没配好。
配置完成后，在命令行输入yalmiptest如果显示一堆测试通过，说明环境没问题。如果里面显示 Cplex 那一项是 failed，十有八九是Matlab版本和Cplex版本不兼容。比如Cplex 12.10不支持太老的Matlab版本，我建议还是用R2020b以上的版本。
3.2 数据准备与参数设定
求解之前，先把基础数据整理清楚。这里我列出一组典型冬季日的示例数据，方便读者对照理解。
| 时段 | 电负荷（MW） | 热负荷（MW） | 风电预测（MW） |
|---|---|---|---|
| 1 | 480 | 350 | 120 |
| 2 | 460 | 360 | 130 |
| 3 | 440 | 380 | 150 |
| 4 | 420 | 400 | 160 |
| 5 | 430 | 390 | 155 |
| 6 | 450 | 370 | 140 |
| ... | ... | ... | ... |
机组参数如下：
| 机组 | P_min (MW) | P_max (MW) | H_min (MW) | H_max (MW) | cv | 煤耗系数a |
|---|---|---|---|---|---|---|
| CHP1 | 80 | 200 | 0 | 150 | 0.2 | 0.25 |
| CHP2 | 60 | 150 | 0 | 80 | 0.25 | 0.30 |
蓄热罐容量设为300兆瓦时，最大充放热功率80MW，初始储热量150兆瓦时。电锅炉最大功率取60MW，效率0.95。
数据为什么要这么设？核心是要让“弃风矛盾”在基准场景下真实存在。从约束可以粗算：深夜时段热负荷400MW，如果由两台机组均摊热负荷（比如CHP1带300，CHP2带100），CHP1的电出力下限是80+0.2300=140MW，CHP2的电出力下限是60+0.25100=85MW，两台机组最小电出力合计225MW。此时电负荷只有420MW，可消纳风电空间只剩420-225=195MW，但风电预测有160MW，这还没问题。如果热负荷更高、风电更大，弃风就会出现。所以一定要先算一遍，确保你的数据真的能产生弃风场景，否则模型优化来优化去，结果只能是无用功。
3.3 核心约束的建模代码与逐段解释
下面给出核心建模代码。整个代码我按“定义变量 -> 写约束 -> 写目标 -> 求解 -> 画图”五段式组织。
%% 1. 基础数据 T = 24; P_load = [480 460 440 420 430 450 470 500 540 560 580 575 560 550 545 560 580 600 590 570 550 530 500 480]; H_load = [350 360 380 400 390 370 360 350 340 335 330 340 350 360 370 380 390 400 410 405 395 385 375 365]; P_w_forecast = [120 130 150 160 155 140 120 100 80 60 50 45 40 50 60 70 80 90 85 100 120 140 150 145]; n_gen = 2; P_chp_min = [80; 60]; P_chp_max = [200; 150]; H_chp_min = [0; 0]; H_chp_max = [150; 80]; cv = [0.2; 0.25]; a_fuel = [0.25; 0.30]; % 蓄热罐 S_max = 300; S_initial = 150; H_ch_max = 80; H_dis_max = 80; eta_ch = 0.95; eta_dis = 0.9; % 电锅炉 P_eb_max = 60; eta_eb = 0.95; % 风电惩罚系数 w_wind = 5000;%% 2. 定义优化变量 P_chp = sdpvar(n_gen, T); % CHP电出力 H_chp = sdpvar(n_gen, T); % CHP热出力 P_w = sdpvar(1, T); % 实际上网风电 P_eb = sdpvar(1, T); % 电锅炉电功率 H_ch = sdpvar(1, T); % 蓄热罐充热 H_dis = sdpvar(1, T); % 蓄热罐放热 S = sdpvar(1, T+1); % 蓄热罐储热量 u = binvar(1, T); % 充放热互斥标志这里sdpvar是Yalmip定义连续变量的函数，binvar定义0-1变量。变量维度我习惯写成分明矩阵：P_chp是2行24列，每一列是一个时段，每一行是不同机组。后面写约束时，按行循环会比矩阵运算更直观，也更容易检查错误。
%% 3. 约束条件 Constraints = []; % 3.1 CHP运行约束 for i = 1:n_gen Constraints = [Constraints, P_chp(i,:) >= P_chp_min(i) + cv(i) * H_chp(i,:)]; Constraints = [Constraints, P_chp(i,:) <= P_chp_max(i)]; Constraints = [Constraints, H_chp(i,:) >= H_chp_min(i)]; Constraints = [Constraints, H_chp(i,:) <= H_chp_max(i)]; end % 3.2 功率平衡 for t = 1:T Constraints = [Constraints, sum(P_chp(:,t)) + P_w(t) == P_load(t) + P_eb(t)]; Constraints = [Constraints, sum(H_chp(:,t)) + H_dis(t) - H_ch(t) + eta_eb * P_eb(t) == H_load(t)]; end这里物理含义是：电力平衡时，系统发电（CHP加风电）要等于电负荷加电锅炉消耗；热力平衡时，系统供热（CHP加蓄热罐放热加电锅炉产热）减去蓄热罐充热等于热负荷。我选择把蓄热罐充热写在等式左边，注意符号别搞反。充热是消耗热量，所以是负号；放热是补充热量，所以是正号。
% 3.3 风电约束 Constraints = [Constraints, 0 <= P_w <= P_w_forecast]; % 3.4 蓄热罐约束 Constraints = [Constraints, 0 <= H_ch <= H_ch_max]; Constraints = [Constraints, 0 <= H_dis <= H_dis_max]; Constraints = [Constraints, 0 <= S <= S_max]; Constraints = [Constraints, S(1) == S_initial]; for t = 1:T Constraints = [Constraints, S(t+1) == S(t) + eta_ch * H_ch(t) - H_dis(t) / eta_dis]; Constraints = [Constraints, H_ch(t) <= H_ch_max * u(t)]; Constraints = [Constraints, H_dis(t) <= H_dis_max * (1 - u(t))]; end Constraints = [Constraints, S(T+1) >= S_initial]; % 3.5 电锅炉约束 Constraints = [Constraints, 0 <= P_eb <= P_eb_max];注意蓄热罐能量递推式的两个效率：充热时是+ eta_ch * H_ch，放热时是- H_dis / eta_dis。也就是说充热有损耗、放热也有损耗。如果效率写反了，相当于给系统赠送能量，结果会低估弃风率。
%% 4. 目标函数 Objective = sum(sum(a_fuel .* P_chp)) + 50 * sum(P_eb) + w_wind * sum(P_w_forecast - P_w); %% 5. 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 1); optimize(Constraints, Objective, ops); %% 6. 结果提取 P_chp_opt = value(P_chp); H_chp_opt = value(H_chp); P_w_opt = value(P_w); P_eb_opt = value(P_eb); S_opt = value(S);value()是Yalmip提取求解结果的核心函数。所有变量计算完成后，把结果重新赋值给普通Matlab矩阵，后续画图和统计都用这些普通矩阵。
3.4 结果可视化：别光盯着一个数字
仿真做完后，我把结果画成了三张图。第一张是电功率平衡堆叠图，横轴是24小时，纵轴是功率，风电画在最下面，CHP各机组叠加起来，最上面是负荷线，电锅炉作为负负荷单独画出来。第二张是热功率平衡图，能看到CHP热出力、蓄热罐充放热、电锅炉产热是如何填补热负荷缺口的。第三张是蓄热罐储热量曲线，这张图最能验证模型逻辑：风电大发时段储热量应该下降（放热），风电小发时段储热量应该回升（充热）。
画图代码不复杂，但有个实用的小技巧，就是给关键时段加垂直线标记。比如凌晨2点到5点风电最大，可以在图上用xline画几条竖线，观察这个时段CHP电出力是不是真的压到了最低，蓄热罐是不是在放热。这种可视化检查比单纯看数据矩阵直观得多。
4. 仿真结果分析：怎么看懂优化结果
4.1 基准场景与优化场景的对比设置
为了说明方法的有效性，我设计了两个算例对比：
- 场景A（基准）：只有CHP机组，没有蓄热罐、没有电锅炉，CHP必须完全跟随热负荷；
- 场景B（优化）：加入蓄热罐和电锅炉，三者在统一优化框架下协调运行。
两种场景都采用同一组风电预测、电负荷和热负荷数据。风电预测曲线设定为凌晨大、白天小，模拟冬季典型大风日。
算完之后看三个指标：弃风电量、弃风率、系统运行成本。弃风电量直接反映消纳效果，运行成本反映经济性。单看弃风率会忽略成本变化，两个指标放一起才能说明方案是否可行。
4.2 结果对比：蓄热罐和电锅炉到底做了什么事
我在这组示例数据下得到的结果是：场景A弃风电量为86.4兆瓦时，弃风率约为3.9%；场景B弃风电量为12.7兆瓦时，弃风率降到了0.6%，风电消纳量接近预测满发水平。系统运行成本反而是场景B更低，因为CHP机组可以减少高煤耗的深度调节出力，电锅炉和蓄热罐的投入成本远低于因弃风造成的损失。
更细看曲线发现，场景A在凌晨4点到6点出现明显弃风，因为这几小时风电预测达到140MW以上，热负荷380~400MW，两台CHP机组电力平衡空间不足以消纳全部风电。场景B中，同样的时段，蓄热罐以80MW的功率放热，替代CHP供热，让CHP电出力比基准场景低了约15MW；电锅炉同时投入约50MW功率，直接吃掉风电，又给热负荷补了约47.5MW热量。两者合计，为风电腾出了将近65MW的消纳空间。
| 指标 | 场景A（基准） | 场景B（蓄热罐+电锅炉） |
|---|---|---|
| 弃风电量（MWh） | 86.4 | 12.7 |
| 弃风率 | 3.9% | 0.6% |
| CHP总煤耗成本（万元） | 132.6 | 128.4 |
| 蓄热罐放热时段 | 无 | 04:00-08:00 |
| 电锅炉投入时段 | 无 | 03:00-06:00 |
这个结果可以清楚说明：蓄热罐和电锅炉的联合配置，不只是“多装了设备”，而是真正改变了CHP机组的运行区间。热负荷被从弃风时段平移到非弃风时段，风电出力空间因此被释放。
4.3 灵敏度分析：蓄热罐容量和电锅炉功率怎么选
模型跑通后，我还做了一组灵敏度分析。蓄热罐容量从0逐步增加到400兆瓦时，其余参数不变，观察弃风率变化。结果是容量在0到150兆瓦时区间内，弃风率下降非常明显；超过200兆瓦时之后，弃风率接近0，曲线进入“平台期”，再加容量收益很小。
电锅炉也有类似的边际递减效应。在蓄热罐容量固定的情况下，电锅炉功率从0加到60MW，弃风率持续下降；超过80MW之后，即便电锅炉还有能力，系统也已无风可弃，再大的功率只会增加空转成本。
这个结论对工程项目很有指导意义：储能和电锅炉的容量配置要结合具体风电特性和热负荷曲线来算，不能拍脑袋选最大。我的经验是先用模型跑几组容量组合，画出“扩容边际收益曲线”，选拐点附近的值最经济。
5. 常见问题与排查技巧实录
5.1 求解器报错与路径配置
新手最容易卡在环境搭建上。常见的报错有No suitable solver found和Solver cplex not found。前者说明Yalmip路径没问题，但没找到任何求解器；后者说明指定的Cplex没有正确加载。排查思路很固定：
- 在命令行输入
which sdpvar，看Yalmip是否被正确加载； - 输入
which cplex，看Cplex的Matlab接口是否在路径里； - 运行
yalmiptest，看Cplex测试项是否通过。
我在一台新电脑上配置时，遇到过Cplex安装目录下有cplex.m文件但Matlab就是找不到的情况，后来发现是解压路径带了中文，Matlab对中文路径兼容不好。改成纯英文路径后一切正常。这个坑建议所有人记一下。
5.2 模型出现infeasible或者求解速度慢
模型报infeasible，先说结论：基本都是约束写矛盾了，不是求解器的问题。最多的情况是电力平衡和机组最小出力冲突。比如某个时段的电负荷特别低，风预测又很大，但你给CHP设了过高的最小电出力，导致等式左右永远不可能相等。快速排查方法：先把风电约束放开到0 <= P_w <= P_w_forecast，如果还不可行，就把风电固定为预测值的80%再算。逐层松绑，很快能定位是哪组约束出的问题。
求解速度慢则主要是二进制变量太多。蓄热罐充放热互斥用的是24个二进制变量，理论上MILP求解涉及分支定界，会比较慢。但实际上这个规模极小，Cplex秒解。如果速度很慢，看看是不是把爬坡约束或储能状态写成了非常复杂的非线性表达。Yalmip虽然支持非线性建模，但非线性问题求解难度大几个数量级。建议始终把模型保持在线性范畴内，遇到非线性项就做分段线性化。
5.3 结果里出现同时充放热，该怎么处理
我在初版模型里没有加互斥约束，结果确实看到了蓄热罐一边充热一边放热的“空转”现象。原因不是求解器傻，而是目标函数里充放热的成本太低，同时充放热虽然消耗能量，却可能帮助等式平衡找到更“便宜”的路径。解决办法就是之前写的，引入二进制变量强制互斥。
不过引入二进制变量之后，有时Cplex为了满足互斥约束，会让设备在很短的时段内来回切换，出现抖振。我的处理方法是给充放热功率变化加一个很小的惩罚项，比如eps * sum(abs(H_ch(t) - H_ch(t-1)))。因为线性模型里绝对值需要引入辅助变量，一般我不太用；更简单的办法是把最优解里出现抖振的时段看一遍，如果只是个别时段，直接人工调整即可。
5.4 风电消纳率上不去，检查这几个地方
如果模型跑完了，但是弃风率还是很高，先别急着加设备容量。我归纳了四个高频原因：
- 热平衡等式里蓄热罐充放热符号写反了，导致放热反而减少了系统可供热量；
- CHP电热耦合约束下限写成了
P_chp >= P_min - cv * H_chp，少了一个负号，耦合方向反了； - 风电预测曲线本身很差，大部分时段不超过20MW，那不管怎么调度都消纳不了多少；
- 电锅炉效率系数用在了热负荷侧而不是电功率侧，导致热平衡计算错误。
我建议每次修改约束后，先打印出约束表达式检查一遍：Constraints，Yalmip会显示完整约束内容。这个功能在调试时非常好用，能一眼看出符号问题。
5.5 关于Matlab版本和工具箱的一些提醒
如果你用的不是正版License，启动Matlab时可能会遇到许可证弹窗异常、工具箱加载失败的问题。这个项目其实完全不依赖Matlab自带的优化工具箱，只需要基础Matlab环境加Yalmip加Cplex即可。如果你的Matlab许可证里没有Optimization Toolbox，我的建议是不要依赖linprog这类官方函数，直接用Yalmip建模，它不检查你装没装官方优化工具箱。
版本方面，我在R2021a到R2023a上都跑过这段代码，均无问题。唯一需要留意的是Cplex版本的兼容性，老版本Cplex对新版Matlab支持往往滞后。如果你用的Matlab太新，可以考虑换用Gurobi求解器，Gurobi对Matlab版本支持和Yalmip集成度也很高。
6. 实操心得：一点扩充建议与个人体会
代码能跑通之后，下一步如果想往论文或实际工程方向扩，我个人建议按这个优先级来扩展。第一是加入更多种类机组的可行性域约束，比如把背压式CHP模型单独列出来；第二是加入热网传输延时和热损模型，这个对消纳计算的影响在真实项目中不容小觑；第三是把目标函数改成多目标，在弃风惩罚、煤耗成本、蓄热罐寿命损耗之间做权衡。
我在实际项目中还有一个心得：不要一开始就追求模型的复杂度，先把“纯CHP + 风电”这个最基础场景跑通，把弃风时段标出来，再逐步加入蓄热罐、电锅炉，每一步都记录弃风量的变化。这样做的好处是，一旦结果不合理，你能明确知道是哪一层设备引入后导致的问题，而不是在一堆复杂约束里找BUG。
最后再分享一个小技巧。在电力平衡等式里，风电实际上网功率变量P_w前面一定不要乘效率系数。风电上网就是上网，它不存在像电锅炉那样的转换效率。少了这一层认知，很容易在建模时给风电加一个0.95之类的效率，结果系统莫名多出损耗，弃风率虚高。这种细节错误不看数据根本发现不了，我在给别的项目做代码审查时就见过不止一次。
