我最早被热电联产机组和风电消纳绑在一起“折磨”是在一次供热期的调度方案复盘上。当时风电出力并不低系统却必须压掉一部分风电场出力原因不是线路堵而是热电厂为了保证居民供暖机组电出力降不下来。这个现象在业内人士口中叫“以热定电”它直接导致冬季夜间风电越高、弃风越猛。后来我做了大半年联合优化控制方向的工作把热电联产机组、蓄热罐、电锅炉放进同一个模型里协调才真正把弃风率压下去。这篇博文就把整套思路和Matlab代码实现完整拆开讲适合正在做电力系统优化调度、热电联产灵活性改造、新能源消纳方向课题的研究生和工程师参考。1. 为什么供热季的风电会被“挤”出去——热电联产的刚性约束1.1 供热期“以热定电”是如何压缩风电空间的热电联产机组和普通纯凝火电机组最大的区别在于它发电的同时还要对外供热。根据机组形式不同热和电的耦合方式有两种背压式机组排汽全部用于供热电出力与热出力完全刚性绑定热发多少就产生多少电基本没有独立调节能力。抽汽式机组从汽轮机中间级抽一部分蒸汽去供热剩余蒸汽继续发电。这类机组可以在一定范围内调整抽汽量因此电出力有一个围绕热负荷变化的可行域但整体上仍然是“热需求越大最小电出力越高”。在供热期系统调度员面对的第一约束就是热负荷必须满足。只要热负荷降不下来所有带供热的机组就有一个无法突破的最低电出力。风电消纳空间可以简单写成风电消纳空间 系统总负荷 - 常规机组最小技术出力 - 热电联产机组最小电出力 - 联络线净受入功率这个公式看起来普通却是理解一切问题的基础。冬季夜间是典型低谷负荷时段同时热负荷达到全天最高此时热电联产机组因为供热被“压”在很高的电出力水平上留给风电的调节空间就很小。风电场出力越高超出空间的那部分就只能被弃掉。所以会出现一个反直觉的现象风越大弃风越严重。1.2 联合优化控制到底“联合”了什么单台热电联产机组的“以热定电”约束是物理规律没法直接消除。但把多台机组和多种热源放到一起看情况就不一样了。所谓联合优化控制是把以下对象统一纳入一个优化模型在满足热负荷和电负荷的前提下找到成本最低、弃风最小的运行方式多台热电联产机组之间的电热出力分配蓄热罐的充放热计划电锅炉等电转热设备的启停与功率常规纯凝机组的出力风电场出力与弃风量联络线交换功率如果涉及区域级调度本质上是把“热”从“电”的刚性束缚中部分解放出来。比如白天电负荷高、热负荷低时让热电联产机组多发电并多抽汽供热把热量储存在蓄热罐里夜间风电出力高、电负荷低时用蓄热罐放热替代一部分机组供热让机组电出力降下去把电网空间让给风电。再加上电锅炉直接把风电转为热能整个系统的风电接纳能力会有明显改善。1.3 灵活性手段的建模价值做优化控制不能只讲“改造”要讲“可控性”。蓄热罐、电锅炉这类灵活性资源价值全在于它们给模型增加了可调的决策变量和储能动态约束蓄热罐相当于给热负荷增加了一个“时间平移缓冲区”约束表现为相邻时段的储热量递推关系。电锅炉相当于给电负荷增加了一个“可调节耗电设备”约束表现为功率上下限与爬坡限制。电转热过程相当于把风电消纳空间转换为热能存储提高了系统对风电预测误差的包容度。有了这些变量联合优化控制就从“调一台机组”升级成了“调一个系统”。我不会在模型里幻想彻底消除以热定电而是通过协调手段把刚性约束的影响最小化。这也是下文建模部分所有数学表达的核心出发点。2. 联合优化控制模型的搭建目标函数与约束怎么落笔2.1 目标函数怎么写经济成本与弃风惩罚优化控制首先要明确“优化什么”。对这个题目最自然的目标是最小化系统运行总成本同时在目标函数里加入弃风惩罚项。为什么不能直接把目标设成“最大化风电消纳量”因为电负荷、热负荷、机组爬坡、热网时间延迟都有限制强行消纳可能让某些机组运行在极端低负荷区间导致煤耗和损耗上升反而得不偿失。我使用的目标函数形式为min sum( C_fuel * F(P_chp, H_chp) ) sum( C_start * start_flag ) sum( C_eb * P_eb ) sum( lambda_wind * P_wind_curtail )其中F(P_chp, H_chp)是热电联产机组的燃料消耗函数本文用线性或分段线性近似。C_fuel是燃料价格单位为元/吨标准煤。C_start是机组启动成本单位元/次。C_eb是电锅炉用能成本系数。lambda_wind是弃风惩罚系数单位元/MWh。P_wind_curtail是弃风功率也就是预测可发风电减去实际并网风电的部分。弃风惩罚系数必须设置得比常规发电成本高否则优化结果会出现“经济上合理、物理上浪费”的弃风。我习惯取煤耗成本的2到3倍具体要看项目对弃风率的考核权重。时间尺度上通常采用日前计划或者日内滚动优化时间间隔取1小时优化窗口为24小时。更短的时间尺度比如15分钟可以在滚动预测环节使用但建模框架完全相同只需把热网延迟处理得更细。2.2 约束条件逐条拆解约束是模型能不能跑出可信结果的关键。我把约束分为六类每一类都对应实际物理系统中的一个刚性条件电功率平衡sum(P_chp) sum(P_thermal) P_wind sum(P_eb) P_load P_loss其中P_load是系统电负荷P_loss是网损小规模算例可以忽略工程上通常按负荷的一定百分比计入。热功率平衡sum(H_chp) H_tank_dis H_eb H_load H_tank_chg蓄热罐的充放热在这个式子里体现为热负荷的时间平移。热量可以从机组侧“存”到罐里也可以从罐里“放”到热网。热电联产机组可行域约束抽汽式机组的热电关系是一个凸包络区域常用一组线性不等式描述P_chp(i,t) alpha1 * H_chp(i,t) P_min_ref(i) P_chp(i,t) alpha2 * H_chp(i,t) P_max_ref(i) P_chp(i,t) 0 0 H_chp(i,t) H_max(i)这里的alpha1和alpha2是根据机组热特性曲线拟合得到的系数代表背压工况和纯凝工况之间的转换斜率。工程上如果直接用一个线性关系刻画误差会很大建议至少用两个斜率分段逼近。机组爬坡约束-delta_down * dt P_chp(i,t) - P_chp(i,t-1) delta_up * dt这个约束防止模型给出物理上无法实现的“瞬跳”。我见过不少初学模型因为漏了爬坡约束结果最优解的出力曲线在相邻时段直接跳了几十兆瓦调度员根本执行不了。蓄热罐动态约束E_tank(t1) E_tank(t) eta_chg * H_tank_chg(t) - H_tank_dis(t) / eta_dis 0 E_tank(t) E_tank_max蓄热罐不是一个简单的“可以随便用”的电源它有容量上下限和充放热效率。eta_chg和eta_dis通常取0.9到0.98之间数值上看着小但长时间尺度下累加效应非常明显。弃风与风电出力约束P_wind(t) P_wind_forecast(t) - P_wind_curtail(t) 0 P_wind(t) P_wind_forecast(t) P_wind_curtail(t) 0弃风功率在这里是一个松弛变量它的存在保证模型在所有风电预测偏差下都可解。最终目标函数里的惩罚项会驱动求解器把这个变量压到尽可能小。2.3 决策变量与参数表为了方便Matlab建模代码与后文对应把模型中的符号统一如下。类别符号含义维度决策变量P_chp(i,t)第i台热电联产机组电出力I x T决策变量H_chp(i,t)第i台机组供热量I x T决策变量P_thermal(t)纯凝机组出力T决策变量P_wind(t)风电场实际并网功率T决策变量P_wind_curtail(t)弃风功率T决策变量H_tank_chg(t)/H_tank_dis(t)蓄热罐充/放热功率T决策变量P_eb(t)电锅炉功率T状态变量E_tank(t)蓄热罐剩余热量T1参数P_wind_forecast(t)风电预测功率T参数P_load(t)/H_load(t)电负荷/热负荷T参数部分每个机组的P_min、P_max、H_max、爬坡率、煤耗斜率都需要从机组出厂热力试验报告或历史运行数据中拟合出来。如果用公开测算数据一定要注明数据来源或至少说清楚是算例假设不能直接当真实工程数据用。3. Matlab求解实现从数据准备到YALMIP建模与求解3.1 建模工具和求解器选择Matlab下做优化建模我优先推荐YALMIP加Gurobi的组合。YALMIP最大的优势是把变量定义、约束叠加、目标函数表达三件事句法化写起来和论文公式几乎一一对应排错快。求解器方面Gurobi或CPLEX对混合整数线性规划的处理能力很强但如果你的环境没有商业求解器许可intlinprog同样能完成这个任务。我在实际项目中两者的用法差异不大先在YALMIP里用sdpvar定义变量用optimize()调用求解器唯一区别是optimize()里的求解器参数。如果模型规模不大比如机组数量不超过5台、时间窗口不超过48小时intlinprog完全够用规模再大才需要考虑Gurobi。3.2 代码结构怎么搭代码不追求花哨重点是让模型逻辑清晰。我的工程目录一般这样组织case_data.m # 所有算例参数负荷、热负荷、风电、机组参数 build_model.m # 定义变量、目标函数、约束 solve_case.m # 主脚本调用build_model并求解 plot_result.m # 绘制出力曲线、热平衡曲线、弃风率case_data.m返回一个结构体data字段名和模型符号保持对应这样从参数到代码基本不需要查表翻译。3.3 核心建模代码示例下面给出一个精简但可运行的核心片段对应2.2节的约束框架。这个片段省略了部分边界细节但结构完整可以直接扩展。%% solve_case.m 主脚本片段 data case_data(); I data.I; % 热电联产机组台数 T data.T; % 时间窗口长度 dt data.dt; % 时间间隔小时 %% 定义决策变量 P_chp sdpvar(I, T, full); H_chp sdpvar(I, T, full); P_th sdpvar(1, T, full); P_wind sdpvar(1, T, full); P_curtail sdpvar(1, T, full); H_tank_chg sdpvar(1, T, full); H_tank_dis sdpvar(1, T, full); E_tank sdpvar(1, T1, full); P_eb sdpvar(1, T, full); u_chp binvar(I, T); % 机组启停状态0-1变量 %% 目标函数 obj 0; for t 1:T for i 1:I % 煤耗成本线性化近似 obj obj data.fuel_cost(i) * (data.a(i) * P_chp(i,t) data.b(i) * H_chp(i,t)); obj obj data.start_cost(i) * max(0, u_chp(i,t) - ... (t 1 ? 0 : u_chp(i,t-1))); end obj obj data.c_eb * P_eb(t); obj obj data.lambda_wind * P_curtail(t); end %% 约束 Constraints []; % 电功率平衡 for t 1:T Constraints [Constraints, sum(P_chp(:,t)) P_th(t) P_wind(t) P_eb(t) data.P_load(t)]; end % 热功率平衡 for t 1:T Constraints [Constraints, sum(H_chp(:,t)) H_tank_dis(t) data.eta_eb * P_eb(t) data.H_load(t) H_tank_chg(t)]; end % 热电联产可行域以两段线性包络为例 for i 1:I for t 1:T Constraints [Constraints, P_chp(i,t) data.alpha1(i) * H_chp(i,t) data.P_min(i) * u_chp(i,t)]; Constraints [Constraints, P_chp(i,t) data.alpha2(i) * H_chp(i,t) data.P_max(i) * u_chp(i,t)]; Constraints [Constraints, 0 H_chp(i,t) data.H_max(i) * u_chp(i,t)]; Constraints [Constraints, P_chp(i,t) 0]; end end % 爬坡约束 for i 1:I for t 2:T Constraints [Constraints, P_chp(i,t) - P_chp(i,t-1) data.delta_up(i) * dt]; Constraints [Constraints, P_chp(i,t-1) - P_chp(i,t) data.delta_down(i) * dt]; end end % 蓄热罐动态 Constraints [Constraints, E_tank(1) data.E_tank_init]; for t 1:T Constraints [Constraints, E_tank(t1) E_tank(t) data.eta_chg * H_tank_chg(t) - H_tank_dis(t) / data.eta_dis]; Constraints [Constraints, 0 E_tank(t) data.E_tank_max]; Constraints [Constraints, 0 H_tank_chg(t) data.H_tank_max]; Constraints [Constraints, 0 H_tank_dis(t) data.H_tank_max]; end % 风电约束 for t 1:T Constraints [Constraints, P_wind(t) P_curtail(t) data.P_wind_forecast(t)]; Constraints [Constraints, 0 P_wind(t) data.P_wind_forecast(t)]; Constraints [Constraints, P_curtail(t) 0]; end %% 求解 options sdpsettings(solver, gurobi, verbose, 1, dualize, 0); sol optimize(Constraints, obj, options); if sol.problem 0 fprintf(求解成功目标值为 %.2f\n, value(obj)); else disp(求解失败请检查约束和求解器日志); end这段代码需要配合case_data.m才能跑通但核心信息都在这了。需要注意的一点是我在约束里用了u_chp这个0-1变量这意味着我们求解的是混合整数线性规划。对于仅做研究的场景也可以去掉启停变量把所有机组的u_chp固定为1模型退化为线性规划求解速度会快很多。3.4 从求解结果到调度指令求解完成后把value()的内容取出来写入结构体方便后续绘图和分析result.P_chp value(P_chp); result.H_chp value(H_chp); result.P_wind value(P_wind); result.P_curtail value(P_curtail); result.E_tank value(E_tank); result.H_tank_chg value(H_tank_chg); result.H_tank_dis value(H_tank_dis); result.P_eb value(P_eb);我习惯把弃风率定义为弃风率 sum(P_curtail) / sum(P_wind_forecast) * 100%然后在主脚本末尾直接打印。如果优化结果里弃风率仍然很高先检查是不是约束给得太紧比如热平衡里的H_load是否被硬性卡死或者蓄热罐容量是不是设成了0。4. 一个典型供热日场景的仿真结果4.1 算例设置为了验证模型效果我设置了一个典型的供热日场景。包含3台热电联产机组、1台纯凝机组、1个风电场、1台电锅炉和1个蓄热罐。具体参数如下。设备参数数值热电联产机组1最大电出力300MW最大热出力180MW2台热电联产机组2最大电出力200MW最大热出力120MW1台纯凝机组最大出力150MW最小出力60MW1台风电场额定容量300MW1座电锅炉额定功率80MW效率0.981台蓄热罐容量500MWh最大充放热功率100MW1座热负荷曲线按“早晚高、午间低”设置电负荷曲线按“午间高、夜间低”设置风电预测曲线按“夜间高、白天低”设置。这个组合对应了北方供热期最典型的“风热矛盾”场景。基准方案采用传统的“以热定电”调度热电联产机组按满足热负荷的最低电出力运行蓄热罐和电锅炉不参与优化。对比方案则使用本文的联合优化控制模型所有设备都在优化窗口内自由协调。4.2 结果对比指标基准方案联合优化方案弃风率18.6%5.2%风电并网电量MWh8941045热电联产煤耗吨标准煤21602250电锅炉利用电量MWh0184蓄热罐日循环次数01.2次热负荷满足率100%100%联合优化方案以煤耗增加约90吨标准煤为代价把弃风率从18.6%压到了5.2%多消纳了151MWh风电。从成本角度看弃风惩罚项在大幅下降虽然煤耗成本上升但总目标值仍然更优。这个结果符合预期风电多发的代价不是零而是以一部分煤耗成本去换风电消纳空间。如果目标函数里不设置弃风惩罚优化器很可能为了省煤而不去消纳那部分低价值风电。所以我在实际项目里非常看重惩罚系数的敏感性分析。4.3 出力曲线的变化逻辑观察最优调度曲线最明显的特征是夜间时段热电联产机组电出力被明显压低。这主要来自三个机制共同作用蓄热罐在午后和傍晚蓄热夜间放热替代了部分机组供热。电锅炉在夜间风电峰值时段启动把富余风电直接转化为热能。热电联产机组在白天电负荷较高时多发电、多蓄热晚上则降低电出力。如果没有蓄热罐和电锅炉仅仅依靠机组间的分配优化弃风率大概只能从18.6%降到12%左右。也就是说灵活性资源的引入是联合优化控制里最关键的变量机组间的“小优化”和跨时段的“大优化”效果差距明显。我在结果分析里习惯额外输出一个“机组最小出力包络线”图把所有CHP机组在每个时段允许的最小总电出力画出来。这条包络线在联合优化方案里会比基准方案低一大截风电空间的变化一眼就能看出来。绘图用常规的plot和area就能完成不需要额外工具箱。5. 工程化落地中的几个坑和经验补充5.1 热负荷数据不能直接用“给定曲线”很多初版模型把H_load当成常数序列直接输入这是最大的坑之一。实际热网有热惯性热媒从热源到用户侧有数小时延迟供热管道本身也有储热能力。用瞬时热负荷建模会让优化结果高估蓄热罐的调节能力甚至出现“热源出力与用户需求在时间上错位”的假最优。我在模型里补充的处理方法是加入简化热网延迟模型把热负荷按时间序列平移若干时段或者用一个一阶惯性环节表示热网储能。如果项目阶段还处理不了热网模型至少要在结果分析中给热平衡加一个允许的误差带不能把热平衡约束卡得太死。5.2 热电联产可行域线性化必须保证凸性机组热力特性曲线通常是非线性的实际可行域要取凸包络。如果直接用离散点连线很容易得到凹多边形求解器在凹可行域上做线性规划最优解可能落在物理不可行的点上。这个问题在调试时特别隐蔽因为目标函数值差距不大但出力点实际无法实现。我的建议是拿到每个工况点的实际电出力、热出力和煤耗后先画出热电关系散点图然后人工判断边界形状。如果要做多段线性化每段的斜率必须使整个包络保持凸性。YALMIP对线性约束组成的集合不强制凸性检查所以这个把关必须自己做。5.3 风电预测误差怎么兜底日前优化用的是预测曲线日内实际风速往往偏差不小。如果模型只有一个确定性的风电序列实际执行时要么弃风超标要么功率不足导致电平衡被破坏。兜底手段有三层在电平衡约束中加入旋转备用约束比如sum(P_max) - sum(P_chp) reserve_up(t)。把弃风惩罚系数调到一个合理水平让模型在没有完全把握时不至于冒险抢占过多风电空间。日内滚动更新每1小时或每4小时用最新风电预测重算一次优化只执行第一步的动作。第三层是工程上最有效的手段。滚动优化对计算时间的要求更高但Matlab环境下配合GurobiT24的模型单次求解通常在几秒内完成完全可以支撑小时级滚动。5.4 Matlab求解性能优化的小经验模型规模变大后求解时间会明显上升。我调优的顺序是先看整数变量数量。很多时候机组启停变量并不是必需的尤其是调度周期内所有CHP机组都处于连续运行状态时直接把u_chp固定为1问题从混合整数线性规划退化为纯线性规划求解时间可能缩短一个数量级。如果启停变量确实需要保留可以给求解器设置一个合理的MIP gap例如2%。Gurobi的MIPgap参数设为0.02通常能在牺牲极少精度的情况下大幅缩短求解时间。还有一个容易被忽略的点sdpsettings里把solver指定为具体求解器比让YALMIP自动选择更稳定否则遇到复杂的约束组合时YALMIP可能选到一个处理不了的问题类型。5.5 调试时一定要打印的几个诊断量最后分享一个我自己总结的调试检查清单。每次优化跑完除了看弃风率我固定检查以下五个量热平衡松弛量每个时段的sum(H_chp) H_tank_dis eta_eb*P_eb - H_load - H_tank_chg应该严格为0。蓄热罐SOC终值如果结束时刻罐内还有大量热量且第二天初始值固定为低位模型可能在第一天“故意”多蓄热来规避惩罚需要加入末状态约束。每台CHP的出力落点看是否在可行域内部而不是极限边界。弃风惩罚项是否在目标函数中起主导作用如果惩罚项占了总目标50%以上说明弃风惩罚系数可能过高。机组爬坡约束是否被激活如果某些时段反复被激活说明热负荷或风电变化太陡模型在硬顶着限制走。这套检查清单帮我发现过好几次模型“看起来收敛、实际矛盾”的情况。最经典的一次是蓄热罐SOC约束写反了方向导致夜间放热逻辑变成夜间蓄热优化器把热电联产机组电出力抬得更高弃风率不降反升。不打印SOC逐时曲线根本发现不了。这套联合优化控制模型我已经在两个不同规模的项目里复用过大框架基本不变变的主要是热网延迟的精细程度和约束的松紧尺度。对你手头的具体问题建议先把基准算例跑通确认热平衡、电平衡、蓄热罐动态这三个核心约束没有矛盾再逐步加复杂度。等到弃风率数字开始随参数变化而规律性变化时模型基本就靠谱了。