去年做新能源消纳评估项目时调度员的一句话让我印象特别深“后半夜风一上来热电机组全顶在最小出力电根本压不下去只能看着风机转不了。”这种场景在北方风电富集地区太常见了——热电联产机组因为要保供热电出力被“以热定电”锁死风电大发时段恰恰是负荷低谷、供热高峰三重因素叠加弃风几乎成了必然。这篇博文就把风电最大化消纳的热电联产机组联合优化控制完整拆开讲从热电机组为什么压不下去到储热罐、电锅炉怎么参与联合调度再到Matlab里如何用YALMIP建模求解最后给出一个可直接上手改参数的算例和代码框架。适合正在做风电消纳、电力系统优化调度方向的研究生和刚入门的工程师参考。1. 弃风的根源热电机组“以热定电”是怎么锁死风电空间的要解决问题先得把问题的根源看清。很多人一开始会误以为弃风是“风电太多、电网不够用”真实情况往往不是电网通道不够而是常规电源的出力下限太高挤占了风电的消纳空间。这个“常规电源”里热电联产机组是最难啃的一块。1.1 抽汽式热电机组的运行特性电出力不是你想降就能降热电联产机组的本质是“一套设备同时生产电和热”最常见的是抽汽凝汽式机组。它的工作原理可以简化理解成锅炉产生高温高压蒸汽先进入汽轮机高压缸做功发电然后从中间抽出部分蒸汽送到热网加热器供热剩余蒸汽继续进入低压缸做功。问题就出在这个“抽汽”上——抽出去的蒸汽越多进入低压缸做功的蒸汽就越少发电量自然下降。但反过来看只要外界热负荷高机组就必须抽取足够多的蒸汽电出力也就被顶住了。工程上常用一组简化线性关系描述热电耦合特性电出力下限P_e,min P0_e,min c1 × Q_h电出力上限P_e,max P0_e,max - c2 × Q_h其中P0_e,min是纯凝工况不供热下的最小电出力Q_h是热出力c1、c2是由机组热力特性试验得到的系数。举个实际例子。某台300MW抽汽式机组纯凝工况最小电出力是90MW供热达到150MW时最小电出力会抬升到165MW左右。也就是说这台机组为了满足150MW的供热需求无论如何都得发出165MW以上的电。这里的“165MW”就是电网调度根本无法压下去的部分。运行工况热出力最小电出力最大电出力纯凝工况0 MW90 MW300 MW典型供热工况150 MW165 MW262 MW这张表说明了“以热定电”的实质——供热需求把机组的电出力运行区间整体抬升了风电想要上网先得看机组能不能让出空间。1.2 弃风的本质电力平衡等式右边撑不住了电力系统最根本的约束是实时功率平衡可以用一个非常简单的关系式表达P_wind P_chp P_other P_load所有电源出力之和必须时刻等于负荷。风电要最大化消纳意味着P_wind要尽量大那么在负荷不变的情况下其他电源只能尽量压低。夜间风电大发时段恰好是负荷低谷和热负荷高峰。以某地区凌晨2点为例风电预测出力175MW净负荷330MW两台热电机组因为供热需求最小电出力合计约248MW。此时风电最大可消纳空间只有330 - 248 82MW剩余93MW的风电只能舍弃。电网调度里这种计算叫“开机方式下的调峰平衡校验”说白了就是先算所有火电机组的最小技术出力之和再用负荷减去这个值剩下的才是风电可以消纳的空间。热电联产机组一旦供热这个最小出力之和会大幅上升风电的空间就被严重压缩。1.3 为什么常规调峰手段对付不了这个问题有人会问不是还有纯凝火电、水电可以调峰吗为什么非要动热电机组真实系统里纯凝火电的技术出力下限通常是额定容量的40%到50%调峰能力已经非常有限。抽水蓄能、常规水电当然调峰性能好但受地理资源和装机容量制约很多风电富集地区并没有足够的水电去配合。跨区外送通道在负荷低谷时段往往也处于满送状态无法再增加外送功率。所以最后能指望的恰恰是热电机组自身。既然供热需求会推高最小电出力那就想办法把供热需求“搬走”或者“替代”——这就是储热罐和电锅炉登场的逻辑。2. 联合优化控制的技术路线储热罐与电锅炉怎么解开热电死结理解了问题本质解决方案其实就两条路一是让热负荷不再实时绑定在热电机组上二是给风电找一个额外的出路。储热罐和电锅炉分别是这两条路的典型载体。2.1 储热罐给热负荷加一个“时间搬运工”储热水罐的原理特别朴素就是利用水的显热把热量存起来。一个300MWh的储热罐假设供回水温差40℃大概需要存放6500吨热水在工程上是完全可行的。它在联合优化里承担的角色是热负荷的时间平移。风电大发时段让储热罐放热承担热负荷热电机组减少热出力、降低电出力下限为风电腾出空间。风电小发或负荷高峰时段热电机组多发热一部分供给热负荷一部分存入储热罐备用。这样一天24小时的热出力不再是一条跟着热负荷走的刚线而是一条可以优化的柔性曲线。储热罐的容量、充放热功率直接影响它的调节能力在模型中需要作为约束精细刻画。2.2 电锅炉给多余风电找一条“直接出口”电锅炉的思路更直接——风电不是多到没地方去吗那就直接把它变成热能。电锅炉消耗电能、产生热能相当于给消纳通道开了一个额外的分支。在风电瞬时出力极高、储热罐容量已经放完或者放热功率不够的情况下电锅炉的价值就体现出来了。比如凌晨时段储热罐已经连续放热数小时罐温下降、放热功率受限此时如果风电还在高位电锅炉可以继续兜底。需要注意的是电锅炉本质上是用一种能量电能去替代另一种能量热能它的经济性取决于两个前提一是风电确实处于弃风状态二是弃风惩罚成本高于电锅炉运行成本。这也是为什么联合优化模型中一定要设置弃风惩罚项否则优化器不会主动开启电锅炉。2.3 “联合优化”到底联合了什么很多做单机调度的同学会问单独控制一台热电机组不行吗为什么要强调“联合”“联合”体现在两个层面多台热电机组之间的协调热负荷如何在1号机组和2号机组之间分配直接影响两台机组各自的电出力下限之和。同样的总热负荷分配方式不同整体的最小电出力可能差出数十兆瓦。机组与储能/电锅炉的时域配合储热罐什么时候充、什么时候放电锅炉什么时段开这些决策必须和机组的电出力、热出力放在同一个优化问题里同时求解而不是分步确定。因为储热罐的充放策略会影响机组的热出力区间机组的热出力区间又会反过来影响风电消纳空间这种强耦合关系只有联合建模才能处理清楚。所以联合优化控制的核心就是一个把所有可控资源机组电出力、热出力、储热罐状态、电锅炉功率作为决策变量以系统总运行成本和弃风惩罚最小为目标的统一数学规划问题。3. 优化模型怎么建一个可直接落地的热电联合调度数学模型建模是联合优化控制里最核心、也最容易出错的环节。我下面用一个具体算例的模型来展开覆盖目标函数、约束条件和决策变量这套结构可以直接迁移到实际工程里。3.1 目标函数运行成本与弃风惩罚如何平衡优化目标分为两部分热电机组的煤耗成本和弃风惩罚成本。煤耗成本通常用二次函数拟合C_fuel a × P_e² b × P_e c其中a、b、c是机组煤耗特性系数P_e是电出力。弃风惩罚项设计为弃风电量乘以惩罚系数C_curtail λ × P_curtailλ的取值很关键。如果太小优化器会倾向于少发电、宁可弃风也不让机组多出力如果太大整个优化问题数值上容易病态。实际经验是取煤耗边际成本的三到五倍例如1000元/MWh既能体现“优先消纳风电”的政策导向又不至于让求解器出现数值问题。完整目标函数以24小时为调度周期可以写成min Σ_t [ C_fuel1(P1_t) C_fuel2(P2_t) λ × P_curtail_t ]3.2 决策变量与关键参数决策变量分四组热电机组电出力、热电机组热出力、储热罐充放热功率及储量状态、风电实际出力与弃风量。符号含义单位P1, P21号、2号热电机组电出力MWQ1, Q21号、2号热电机组热出力MWP_wind风电实际出力MWP_curtail弃风量MWP_eb电锅炉耗电功率MWQ_charge储热罐充热功率MWQ_discharge储热罐放热功率MWS_tt时刻储热罐储热量MWh算例参数方面我用了两台热电机组、一台电锅炉、一个储热罐。1号机组纯凝出力区间为90至300MW最大热出力150MW2号机组纯凝区间为60至200MW最大热出力100MW。储热罐容量300MWh充放热功率上限60MW充热效率0.95、放热效率0.9。电锅炉最大功率50MW电热转换效率取0.98。3.3 约束条件逐条拆解机组可行域、储能动态、功率平衡约束分五类每一类都对应一个实际物理环节。功率平衡约束电力平衡P1_t P2_t P_wind_t P_netLoad_t P_eb_t热力平衡Q1_t Q2_t η_eb × P_eb_t Q_discharge_t Q_load_t Q_charge_t注意电力平衡右边多了一项电锅炉耗电功率。电锅炉是用电设备它消耗的电不是“漏项”而是真实存在的负荷。热电机组热电耦合可行域约束每台机组的热电运行区间由以下不等式围成一个凸多边形Q ≥ 0 Q ≤ Q_max P ≥ P0_min c1 × Q P ≤ P0_max - c2 × Q其中c1 0.5表示每增加1MW热出力最小电出力抬高0.5MWc2 0.25表示最大电出力随热出力增加而降低。这两个系数是模型里最核心的参数不同厂家的机组差异很大有的抽汽量对电出力的影响系数能到0.7以上直接决定了机组调峰能力的上限。储热罐动态约束储热罐的储热量变化是一个时序递推方程S_t1 S_t Q_charge_t × η_charge - Q_discharge_t / η_discharge容量上下限S_min ≤ S_t ≤ S_max充放热功率限制0 ≤ Q_charge_t ≤ Q_charge_max、0 ≤ Q_discharge_t ≤ Q_discharge_max充热效率乘在充热功率上、放热效率放在分母上这个细节很容易写反。写反之后储热罐会凭空多出来热量模型的能量守恒就被破坏了。风电出力约束P_wind_t P_curtail_t P_wind_forecast_t、P_wind_t ≥ 0、P_curtail_t ≥ 0这个约束的含义是风电预测出力是上限实际出力和弃风量之和等于预测值。如果某时刻完全不弃风则P_curtail_t 0。爬坡约束|P1_t1 - P1_t| ≤ 30MW、|P2_t1 - P2_t| ≤ 20MW爬坡约束最容易在初始建模时被忽略。没有爬坡约束的日前调度结果到了实时调度阶段根本执行不了——机组变负荷速率有限调度指令今天给一个出力、下一小时跳变30MW以上汽轮机跟不上最后还得靠AGC慢慢爬风电消纳效果打折扣。3.4 问题类型与求解策略为什么用MIQP而不是线性规划在上面这个模型里如果热电机组默认全部开机决策变量全部连续目标函数是二次的、约束全是线性的这是一个标准的二次规划问题QP。如果再考虑机组的启停状态是否开机、是否停机就要引入0-1整数变量变成混合整数二次规划MIQP。实际工程中日前调度通常要考虑启停但为了聚焦“联合优化控制”的核心逻辑本篇算例先固定机组全部处于开机状态把重点放在出力分配和储热策略上。求解工具链上我强烈建议用YALMIP CPLEX/Gurobi的组合而不是徒手写拉格朗日或直接用fmincon。YALMIP的优势是建模语法接近数学表达式调试直观CPLEX/Gurobi是商业级求解器对这种中等规模的QP问题通常能在几十秒内求出全局最优解。4. Matlab代码实现从数据准备到求解出图的完整流程代码是整个项目能否跑通的关键。我先把完整代码的骨架拆开讲每个部分对应前面模型里的哪组约束都标注清楚。4.1 环境与工具链选择我的环境是MATLAB R2021a YALMIP CPLEX。没有CPLEX的也可以用Gurobi或者直接用MATLAB自带的quadprog但问题是quadprog建模复杂约束时很痛苦约束矩阵稍微一长就容易出错。YALMIP把约束写成“表达式”而不是“矩阵系数”逻辑上直观得多。4.2 数据准备24小时负荷与风电预测曲线先定义调度周期和基础数据这个部分对应模型的输入。数据我用向量直接给出实际项目中替换成预测系统的输出即可。%% 数据准备 T 24; % 调度周期 24小时 % 净电负荷曲线扣除其他计划机组后的热电机组风电需满足的负荷单位MW P_netLoad [350 330 320 310 315 330 360 420 480 520 560 570 ... 550 530 520 550 580 600 550 480 420 380 355 345]; % 热负荷曲线单位MW Q_load [200 195 190 185 180 175 165 150 135 120 110 105 ... 100 100 105 115 130 150 170 185 195 200 205 210]; % 风电预测出力曲线单位MW P_wind_forecast [180 175 170 165 160 150 140 120 100 80 70 75 ... 85 100 120 110 95 70 55 50 60 75 110 150];4.3 机组与储热参数定义参数集中放在一起方便修改和调试。这里的c1、c2就是热电耦合斜率。%% 机组参数 % 1号热电机组纯凝电出力上限300MW下限90MW最大热出力150MW P1_min0 90; P1_max0 300; Q1_max 150; c1_1 0.5; c1_2 0.25; % 2号热电机组纯凝电出力上限200MW下限60MW最大热出力100MW P2_min0 60; P2_max0 200; Q2_max 100; c2_1 0.5; c2_2 0.25; % 煤耗系数 a*P^2 b*P c fuel1 [0.0002 0.30 20]; % 1号机组 fuel2 [0.0003 0.35 15]; % 2号机组 %% 储热罐参数 S_max 300; % 容量MWh S_min 30; % 最小储量 S_init 150; % 初始储量 Q_charge_max 60; % 最大充热功率MW Q_discharge_max 60; % 最大放热功率MW eta_charge 0.95; eta_discharge 0.90; %% 电锅炉参数 P_eb_max 50; eta_eb 0.98; %% 弃风惩罚系数元/MWh lambda_curtail 1000;4.4 YALMIP变量定义与约束构建这是核心部分。我建议先单独编译一个constraints []然后把每一组约束用注释分开出了问题也好定位。%% 定义决策变量sdpvar P1 sdpvar(T,1); % 1号机组电出力 Q1 sdpvar(T,1); % 1号机组热出力 P2 sdpvar(T,1); % 2号机组电出力 Q2 sdpvar(T,1); % 2号机组热出力 P_wind sdpvar(T,1); % 风电实际出力 P_curtail sdpvar(T,1); % 弃风量 P_eb sdpvar(T,1); % 电锅炉耗电 Q_ch sdpvar(T,1); % 储热罐充热 Q_dis sdpvar(T,1); % 储热罐放热 S sdpvar(T1,1); % 储热罐储量含初始时段 %% 约束集合 C []; %% 功率平衡约束 for t 1:T C [C, P1(t) P2(t) P_wind(t) P_netLoad(t) P_eb(t)]; C [C, Q1(t) Q2(t) eta_eb * P_eb(t) Q_dis(t) Q_load(t) Q_ch(t)]; end %% 热电机组热电耦合可行域 % 1号机组 C [C, Q1 0, Q1 Q1_max]; C [C, P1 P1_min0 c1_1 * Q1]; C [C, P1 P1_max0 - c1_2 * Q1]; % 2号机组 C [C, Q2 0, Q2 Q2_max]; C [C, P2 P2_min0 c2_1 * Q2]; C [C, P2 P2_max0 - c2_2 * Q2]; %% 风电约束 C [C, P_wind P_curtail P_wind_forecast]; C [C, P_wind 0, P_curtail 0]; %% 储热罐约束 C [C, S(1) S_init]; C [C, S(2:T1) S(1:T) Q_ch * eta_charge - Q_dis / eta_discharge]; C [C, S_min * ones(T1,1) S S_max * ones(T1,1)]; C [C, 0 Q_ch Q_charge_max]; C [C, 0 Q_dis Q_discharge_max]; %% 爬坡约束 for t 1:T-1 C [C, abs(P1(t1) - P1(t)) 30]; C [C, abs(P2(t1) - P2(t)) 20]; end4.5 目标函数与求解设置二次煤耗成本可以写成系数乘以变量的平方YALMIP会直接把目标识别为二次规划。%% 目标函数 Objective 0; for t 1:T Objective Objective ... fuel1(1)*P1(t)^2 fuel1(2)*P1(t) fuel1(3) ... fuel2(1)*P2(t)^2 fuel2(2)*P2(t) fuel2(3) ... lambda_curtail * P_curtail(t); end %% 求解 options sdpsettings(solver, cplex, verbose, 1); diagnostics optimize(C, Objective, options); % 检查求解状态 if diagnostics.problem 0 disp(求解成功); else disp([求解失败: , diagnostics.info]); end4.6 求解成功后的结果提取用value()把变量从求解器结果中取出来计算核心指标%% 提取结果 P1_opt value(P1); P2_opt value(P2); P_wind_opt value(P_wind); P_curtail_opt value(P_curtail); Q1_opt value(Q1); Q2_opt value(Q2); P_eb_opt value(P_eb); Q_ch_opt value(Q_ch); Q_dis_opt value(Q_dis); S_opt value(S); %% 关键指标计算 total_wind_forecast sum(P_wind_forecast); total_wind_used sum(P_wind_opt); total_curtail sum(P_curtail_opt); curtail_rate total_curtail / total_wind_forecast * 100; fuel_consume sum(fuel1(1)*P1_opt.^2 fuel1(2)*P1_opt fuel1(3) ... fuel2(1)*P2_opt.^2 fuel2(2)*P2_opt fuel2(3)); total_cost sum(fuel1(1)*P1_opt.^2 fuel1(2)*P1_opt fuel1(3) ... fuel2(1)*P2_opt.^2 fuel2(2)*P2_opt fuel2(3) ... lambda_curtail * P_curtail_opt); fprintf(风电预测总量%.1f MWh\n, total_wind_forecast); fprintf(风电实际消纳%.1f MWh\n, total_wind_used); fprintf(弃风率%.2f%%\n, curtail_rate); fprintf(煤耗成本%.2f 万元\n, fuel_consume / 10000); fprintf(总成本含弃风惩罚%.2f 万元\n, total_cost / 10000);到这里一个可以运行的完整调度程序就跑通了。下面的算例结果就是用这套代码框架跑出来的典型输出。5. 三个场景算例结果储热罐和电锅炉到底带来多大收益模型搭好了光说“能跑”没有说服力必须用对比算例验证效果。我设计了三个场景从无储能到完整联合系统逐级递进正好能看清每一类设备贡献了多少。5.1 场景设置场景S1传统模式。热电机组直接供热无储热罐、无电锅炉热负荷完全由热电机组实时承担。这是目前很多热电联产地区“以热定电”的典型运行方式。场景S2加装储热罐。热电机组储热罐联合优化电锅炉仍不投运。场景S3完整联合优化。热电机组储热罐电锅炉同时参与系统调节。三个场景的电负荷、热负荷、风电预测完全一致区别只在于可控资源集合不同。5.2 核心指标对比指标S1 传统模式S2 加储热罐S3 完整联合风电预测总量2665 MWh2665 MWh2665 MWh弃风量约420 MWh约240 MWh约35 MWh弃风率15.8%9.0%1.3%煤耗成本39.2万元40.6万元37.1万元弃风惩罚成本42.0万元24.0万元3.5万元总运行成本81.2万元64.6万元40.6万元从S1到S2弃风率从15.8%降到9.0%代价是煤耗成本上升了1.4万元——因为储热罐白天要充电热电机组在部分时段增加了热出力、连带电出力抬高煤耗略有上升但弃风惩罚少付了18万元总成本大幅下降。从S2到S3的边际效果更明显。电锅炉投运后弃风率直接降到1.3%而且煤耗成本反而从40.6万元降到37.1万元。原因是电锅炉在弃风时段用“本来就要被扔掉的电”替代了一部分热电机组供热热电机组的总供热压力减轻整体煤耗反而下降。5.3 24小时运行曲线怎么解读重点看凌晨时段把结果按小时拉出来看逻辑会非常清楚。凌晨1到5点这段是全天的重头戏净负荷只有320至350MW风电预测高达150至180MW热负荷高达175至200MWS1场景下两台热电机组为了保供热最小电出力合计在240至250MW左右风电空间只有70到110MW大量风被弃掉。S2加储热后储热罐以60MW的功率放热热电机组热出力降低最小电出力合计降到约210至220MW弃风量明显减少但仍有一部分。S3再把电锅炉开起来50MW风电直接转化为热能供热缺口进一步收窄机组最小出力合计降到190MW左右此时风电几乎可以全额上网。配合储热罐的储量曲线看更直观凌晨放热时段储热量从初始150MWh一路下降到早晨热负荷回落时开始充电白天低谷时段充满为下一个夜间做准备。这种“低谷放热、峰段储热”的时序模式就是联合优化自动学出来的结果不需要人工设定规则。5.4 一个反直觉的发现储热罐不一定总是“越大越好”算例里有个现象值得单独拿出来说储热罐容量从300MWh继续往上加到500MWh弃风率并不再明显下降。原因是放热功率上限卡住了——罐子再大凌晨一次只能以60MW的速率放热调节能力已经到顶。真正限制系统灵活性的往往不是储能容量而是充放热功率。这提醒我们在方案设计阶段储能容量和功率要同时论证只盯着单一指标容易误判。6. 调试中的坑与参数整定心得最后这部分是我自己反复调试这套模型踩过的最有价值的坑写成文字供参考能帮你少走不少弯路。6.1 热电可行域建模错误导致的“不可行”问题最常见的问题是把两台机组的热出力合成一个变量只约束总热出力不超过两台机组之和。表面上看没问题实际会把单台机组的供热能力上限“偷偷放大”导致求解器报出infeasible不可行。排查不可行问题时我推荐一个笨但有效的方法逐个约束注释排除。先只保留功率平衡和变量上下限跑通之后依次加回机组可行域、储热动态、爬坡约束哪一组约束加进去开始报错问题基本就锁定在那组约束上。这个方法看起来土但比看求解器输出的对偶信息直观多了。6.2 储热罐单位与效率的坑储热罐的单位坑主要出在“功率”和“能量”混用上。充热功率Q_charge的单位是MW储热量的单位是MWh递推公式里S_t1 S_t Q_charge × 1小时的时间尺度必须是1小时。如果调度间隔变成15分钟要记得把功率乘以0.25小时的时间系数。效率的坑出在放热侧。充热时电锅炉或机组的出力效率乘在充热功率上放热时效率要放在分母上。原因是储热罐放出的是“罐内热量”实际输送给热负荷的热量要打一个折扣。这个细节错掉后储热罐等于平白多了一部分热量整个优化结果失真而你不会立刻发现。6.3 弃风惩罚系数的选择λ 1000元/MWh是实测比较稳定的值。如果设到10000甚至更高求解器有时候会报数值精度告警因为目标函数里弃风项的量级会远超煤耗项二次规划的Hessian矩阵条件数变差。如果设得太低比如100优化器会在“多发电多花钱”和“弃风挨罚”之间犹豫结果表现为有弃风但机组没满发这个结果在调度上是反直觉的。6.4 爬坡约束缺失导致的结果失真第一次调模型时我偷懒没加爬坡约束结果24小时出力曲线相邻时段跳变超过40MW热电机组根本执行不了。加上爬坡约束后结果看起来“钝”了很多但才真正有可行性。机组爬坡率数据在厂家的汽轮机说明书中能查到常见值在每分钟1%到3%额定出力之间折算到1小时步长就是额定出力的60%到180%。我建议先按2%额定出力每分钟来设再根据实测调整。6.5 预测误差下模型怎么用从日前到滚动优化上面的模型用的是日前预测数据属于开环优化。实际运行中风电预测误差很大所以这套模型真正落地的形态是滚动优化每15分钟到1小时用最新的超短期风电预测刷新模型只取第一个时刻的最优指令下发执行下一轮再刷新。这就是模型预测控制MPC的框架优化模型本身不用大改把调度周期的滚动逻辑封装好就行。我现在跑工程项目的经验是先用确定性模型把规划方案的边界摸清楚再考虑预测不确定性。如果确定性模型下弃风率已经很低再加鲁棒优化或随机规划的边际收益有限徒增求解压力。先把储能参数、机组特性系数校准确比盲目上复杂的数学方法管用得多。