做电力系统优化调度尤其是冬季供暖期的运行分析时八成以上的人都会碰到同一个痛点热负荷一上来热电联产机组就被“以热定电”摁得死死的电网的调峰空间被压得很窄风电场那边却经常是大风满发。这时候如果你手头只有常规火电的调度模型算出来的弃风量往往会高得离谱但又说不清楚问题到底出在哪。我去年做某区域电网冬季调峰方案时就把抽凝式热电联产机组、背压式热电联产机组、蓄热罐和电锅炉放到同一个优化框架里用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(t1) 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风电预测MW148035012024603601303440380150442040016054303901556450370140............机组参数如下机组P_min (MW)P_max (MW)H_min (MW)H_max (MW)cv煤耗系数aCHP18020001500.20.25CHP2601500800.250.30蓄热罐容量设为300兆瓦时最大充放热功率80MW初始储热量150兆瓦时。电锅炉最大功率取60MW效率0.95。数据为什么要这么设核心是要让“弃风矛盾”在基准场景下真实存在。从约束可以粗算深夜时段热负荷400MW如果由两台机组均摊热负荷比如CHP1带300CHP2带100CHP1的电出力下限是800.2300140MWCHP2的电出力下限是600.2510085MW两台机组最小电出力合计225MW。此时电负荷只有420MW可消纳风电空间只剩420-225195MW但风电预测有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, T1); % 蓄热罐储热量 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(t1) 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(T1) 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蓄热罐电锅炉弃风电量MWh86.412.7弃风率3.9%0.6%CHP总煤耗成本万元132.6128.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那不管怎么调度都消纳不了多少电锅炉效率系数用在了热负荷侧而不是电功率侧导致热平衡计算错误。我建议每次修改约束后先打印出约束表达式检查一遍ConstraintsYalmip会显示完整约束内容。这个功能在调试时非常好用能一眼看出符号问题。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之类的效率结果系统莫名多出损耗弃风率虚高。这种细节错误不看数据根本发现不了我在给别的项目做代码审查时就见过不止一次。