前阵子帮一个北方省份的热电厂做冬季运行方式分析调度那边抛过来一个数据某日夜间风电场理论出力接近满发实际并网却不到三成。电话那头电厂值班员也很无奈——热网供水温度不能降供热机组的电出力被“绑”死了想给风电让路都让不出来。这个场景搞过电力调度的人都懂这就是“以热定电”和风电消纳之间的矛盾。这篇文章就聊聊我做“风电最大化消纳的热电联产机组联合优化控制”这套东西的经验包括背后的物理逻辑、优化模型怎么建、Matlab代码怎么组织以及真正跑算例时踩过的几个坑。适合正在做毕业设计、综合能源或储能方向研究以及刚接手热电联产调度优化项目的朋友参考。1. 冬季弃风是怎么来的“以热定电”的耦合困局与三条解耦路子1.1 机组电热特性电出力为什么会被热负荷锁死热电联产机组跟普通凝汽式机组最大的区别是它既要发电又要供热。抽汽式机组从汽轮机中间抽出一部分蒸汽去加热热网水这部分蒸汽不经过低压缸做功所以机组的电出力上限会随着供热抽汽量的增加而下降。用简化公式表达就是P_max(H) P_max0 - cv × H这里H是热出力cv是电热耦合系数通常在0.5到0.9之间也就是说每增加1MW热出力最大电出力要损失大概0.5到0.9MW。更麻烦的是下限。机组为了保证低压缸最小冷却流量存在一个最小电出力约束。在供热工况下这个最小电出力还会随热出力上升而抬高近似可以写成P_min(H) P_min0 ch × HP_min0大约是纯凝工况的最小技术出力一般在机组额定容量的30%左右ch在简化模型里取0.3到0.5都有人用。这么一来冬季热负荷越高机组“最低也得发的电”就越高。这就是“以热定电”的物理根源热负荷决定了机组最小电出力这部分电量优先上网风电只能挤剩下的空间。我给刚接触这个方向的同学讲的时候喜欢打一个比方热电机组就像一个“连体人”一头连着热网一头连着电网。热网那头站着一个刚性热负荷它往左拉电网这头的电出力就得跟着动没法单独喊停。储热罐和电锅炉的作用就是把这组“连体人”做手术分开。1.2 弃风的本质热负荷高峰撞上风电高峰风电在北方冬季有一个很反直觉的特点夜间风速通常更大出力更高但夜间恰恰是热负荷的高峰期同时又是电负荷的低谷期。三个条件叠在一起结果就是热电机组因为热负荷高而被顶高的最小电出力很难降下来机组总出力很容易超过电负荷需求多余的电只能弃掉。我做算例时常用一组典型数据夜间电负荷降到280MW热负荷却有180MWth热电机组最小电出力按上面的公式算出来要150MW以上风电场半夜大风天理论出力240MW两者加起来远超负荷弃风率轻松超过40%。这不是极端场景而是北方采暖期几乎每天都会出现的常规操作。所以单纯指望“风大的时候少烧点煤”根本不解决问题机组供热工况下的电出力下限就是压不下来。要让风电进来只有两条路要么降低机组的有效热出力要么给系统增加额外的电负荷。这两条路对应着下面要说的储热罐和电锅炉。1.3 三条解耦途径储热罐、电锅炉、热网蓄热目前工程上做热电解耦无外乎三种手段。第一种是储热罐。水蓄热罐把热网回水和供水之间多出来的热水存起来实现热量的时间平移。关键点在于它可以把“热负荷高峰时段的供热任务”挪到“热负荷低谷时段”去完成从而降低热电机组在夜间的热出力也就压低了夜间的最小电出力给风电腾出空间。这一点和很多人第一直觉正好相反——在风电消纳这个目标下储热罐是白天蓄热、夜间放热而不是夜间蓄热、白天放热。这个逻辑我后面算例里会详细算给你看。第二种是电锅炉。电锅炉直接消耗电能制热本质上是给电网增加一个可控电负荷。夜间风大用不掉时让电锅炉把风电转成热能送去热网同时替代部分机组供热进一步压低机组热出力。电锅炉的优点是响应快、投资相对小缺点是综合效率再高也只有95%左右运行时相当于“用高档电换低档热”经济性要靠弃风电量来撑。第三种是热网蓄热。利用城市热力管网本身的蓄热能力通过调整供水温度来短期改变热负荷相当于把整个热网当成一个大储热罐。这个思路很诱人但实际控制起来很麻烦热网延迟、多热源耦合、用户侧反馈模型阶数一高就容易失控。我个人的建议是初版优化模型先别碰热网动态特性把储热罐和电锅炉做进去效果已经足够明显。三种方式怎么选我做了个表方便对比手段投资控制复杂度响应速度适用场景储热罐中高中低分钟级已有集中热网、热负荷峰谷差大的电厂电锅炉低中低秒级弃风严重的风电场附近、夜间电价低的系统热网蓄热低利用存量高小时级热网规模大、有成熟动态模型的单位标题里说的“联合优化控制”本质上就是把机组的电出力、热出力、储热罐的充放热功率、电锅炉的功率放到同一个优化框架里去协调谁承担热负荷、谁承担电负荷都由模型在每个时段统一决策。下面这一章就把这个优化模型的骨架给你拆开讲。2. 优化模型怎么搭决策变量、目标函数与约束条件的取舍2.1 决策变量把“发多少电、供多少热、存多少热”分开定义写优化模型前先把变量清单列清楚。基于前面说的解耦手段模型里至少要包含下面几类变量热电机组电出力 P_chp(t)热出力 H_chp(t)风电实际上网功率 P_w(t)理论可用出力 P_wf(t)储热罐充热功率 H_ch(t)、放热功率 H_dis(t)以及罐内蓄热量 E_s(t)电锅炉用电功率 P_eb(t)对应产热功率 H_eb(t)如果还考虑机组启停那就要再引入0/1变量 on(t)。我建议初版模型先别碰0/1变量直接假设机组全天在线只优化出力分配。这一步跑通以后再去扩展启停优化否则模型一旦不可行你很难分辨到底是数据问题还是逻辑问题。这一点我是吃过亏的。最开始我直接上了MILP又是启停又是爬坡结果某天负荷曲线有一点波动求解器直接报infeasible排查了两天才发现是某一时段机组最小出力加上风电预测下限已经超过电负荷了。把启停变量去掉、先把连续模型跑通是新手最容易省掉却最该保留的一步。2.2 目标函数经济目标里加一个弃风惩罚项目标函数怎么定决定了优化结果偏向谁。纯最大化风电消纳当然可以但调度侧真正认的是运行成本。更稳妥的做法是以系统总运行成本最小为目标其中煤耗成本取机组电出力和热出力的函数同时在目标里加入弃风惩罚项min F Σ C_fuel(P_chp, H_chp) Σ C_eb × P_eb Σ C_wind × (P_wf - P_w)这里的 C_wind 是弃风惩罚系数单位是元/MWh。它不直接对应某个真实账目而是用来告诉模型“弃风有多亏”。工程上常见的取值是高于煤耗成本、低于或接近风电上网电价这样模型会优先消纳风电但不会为了多消纳一点风而不顾一切地增加电锅炉耗电。煤耗成本函数在学术论文里通常写成二次函数在Matlab里处理起来稍麻烦可以按分段线性化处理如果只是做方案对比分析用线性煤耗系数也够用。我自己的经验是先跑一版线性煤耗的LP模型把所有逻辑调通、结果能解释清楚再慢慢把煤耗曲线加细、把非线性部分线性化。一上来就追求高精度往往得不偿失。2.3 约束条件功率平衡、机组可行域、储热罐动态一个都不能少约束是整个模型的核心也是最容易漏项的地方。我按建模顺序把必备约束列一下。电功率平衡。P_chp(t) P_w(t) P_load(t) P_eb(t)。注意电锅炉是消耗电能的放在负荷侧你习惯把可调电源放左边、净负荷放右边也行关键是符号别反。热功率平衡。H_chp(t) H_dis(t) H_eb(t) H_load(t) H_ch(t)。放热是供给充热是需求电锅炉产热也是供给这个式子直观且不容易错。机组电热可行域。前面已经给出上下限这里再整合一遍P_max0 - cv × H_chp(t) ≥ P_chp(t) ≥ P_min0 ch × H_chp(t)同时 H_chp 还有一个上限也就是机组最大供热能力。风电出力约束0 ≤ P_w(t) ≤ P_wf(t)。风电可以少发但不能超发少了的部分就是弃风弃风量表达式就是 P_wf(t) - P_w(t)。储热罐状态转移方程E_s(t1) E_s(t) (η_ch × H_ch(t) - H_dis(t) / η_dis) × Δt以及储热罐容量边界、充放热功率边界E_min ≤ E_s(t) ≤ E_max0 ≤ H_ch(t) ≤ H_ch_max0 ≤ H_dis(t) ≤ H_dis_max这里 η_ch 是充热效率η_dis 是放热效率一般水蓄热罐充放效率在0.9到0.95之间。还有一个很实用的约束是同一时段不能同时充热和放热。如果求解器支持可以用 sos 约束或者一个大M约束处理不过实际算例里只要目标函数没有刻意奖励“边充边放”最优解一般不会出现这种状态。我通常会用一个小的互斥约束把风险堵死代码里会体现。电锅炉约束0 ≤ P_eb(t) ≤ P_eb_maxH_eb(t) η_eb × P_eb(t)。另外别忘了风电消纳比例或者最小上网电量这类调度要求有的话就加一条线性约束没有就跳过。到这里一个可以交给求解器的优化模型就成型了。变量全是连续的这就是一个线性规划问题Matlab Yalmip 几分钟就能建完。下面这一章就讲具体代码是怎么写的。3. Matlab代码实现从Yalmip建模到Cplex求解的完整链路3.1 环境准备Yalmip和Cplex/Gurobi的搭配我的主力工具是 Matlab R2021b Yalmip最新版 Cplex 12.10。为什么选这个组合Yalmip只是一个建模语言把你的优化问题翻译成求解器认识的格式真正负责算的是底层的 Cplex 或 Gurobi。Matlab自带的linprog也能跑但对带分段线性化、0/1变量以后的问题支持弱一些建议直接上Cplex或Gurobi。新手上手最常见的坑是装了Yalmip但没把求解器的路径配好运行 optimize 的时候提示 No suitable solver found。解决办法很简单把Cplex安装目录下的 cplex/bin 路径加到Matlab的 setpath 里然后在命令行输入 cplex.getVersion能输出版本号就说明认到了。另外注意 Matlab 2022a 以上版本对某些老版本Cplex不兼容建议用 Cplex 12.10 以上或者 Gurobi 9.5 以上。这个问题在Matlab中文社区里被问烂了但每次还是有人卡在第一步。3.2 输入数据组织日期、负荷曲线、机组参数怎么准备我习惯把所有时序数据放在一个Excel文件里第一列是时间后面依次是电负荷、热负荷、风电预测出力。Matlab里用 readtable 读进来转成数组。机组参数和储热参数单独放在一个结构体里免得每次换场景时改代码。时间粒度方面日计划优化一般取1小时一个点一天24个时段模型规模小、求解快运行层面做滚动优化时前4小时取15分钟粒度、后面取1小时粒度精度和速度兼顾。初版代码我建议直接写24时段1小时粒度先把主链条打通。3.3 核心代码骨架变量定义、约束、求解、结果提取下面这段代码是一个可以运行的简化版模型里包含一台热电机组、一个风电场、一个储热罐和一台电锅炉。我在关键位置都加了注释。%% 核心建模代码骨架 % 假设已经有数据P_load, H_load, P_wf (1x24) % 机组参数 P_max0 300; % 纯凝最大电出力 MW P_min0 90; % 纯凝最小技术出力 MW cv 0.35; % 电出力上限随热出力下降系数 ch 0.35; % 电出力下限随热出力上升系数 H_max 200; % 最大热出力 MWth % 储热罐参数 E_max 400; E_min 20; E0 200; E_end 200; H_ch_max 120; H_dis_max 120; eta_ch 0.95; eta_dis 0.95; % 电锅炉参数 P_eb_max 80; eta_eb 0.98; % 成本系数 c_coal 280; % 煤耗成本 元/MWh简化线性 c_wind 400; % 弃风惩罚 元/MWh c_eb 80; % 电锅炉运行成本 元/MWh电 T 24; dt 1; % 定义决策变量 P_chp sdpvar(1, T); % 机组电出力 H_chp sdpvar(1, T); % 机组热出力 P_w sdpvar(1, T); % 风电实际上网 P_eb sdpvar(1, T); % 电锅炉功率 H_eb sdpvar(1, T); % 电锅炉产热 H_ch sdpvar(1, T); % 储热充热 H_dis sdpvar(1, T); % 储热放热 E_s sdpvar(1, T1); % 储热罐蓄热量 % 约束集合 Cons []; % 功率平衡 Cons [Cons, P_chp P_w P_load P_eb]; % 电平衡 Cons [Cons, H_chp H_dis H_eb H_load H_ch]; % 热平衡 % 机组电热可行域 Cons [Cons, P_max0 - cv*H_chp P_chp]; Cons [Cons, P_chp P_min0 ch*H_chp]; Cons [Cons, 0 H_chp H_max]; % 风电 Cons [Cons, 0 P_w P_wf]; % 电锅炉 Cons [Cons, 0 P_eb P_eb_max]; Cons [Cons, H_eb eta_eb * P_eb]; % 储热罐动态 Cons [Cons, E_s(2:end) E_s(1:end-1) dt*(eta_ch*H_ch - H_dis/eta_dis)]; Cons [Cons, E_min E_s E_max]; Cons [Cons, 0 H_ch H_ch_max]; Cons [Cons, 0 H_dis H_dis_max]; Cons [Cons, E_s(1) E0]; % 初始蓄热量 Cons [Cons, E_s(end) E_end]; % 日末回位防止跑偏 % 目标函数煤耗 电锅炉成本 弃风惩罚 obj c_coal * sum(P_chp) c_eb * sum(P_eb) c_wind * sum(P_wf - P_w); % 求解 ops sdpsettings(solver, cplex, verbose, 1); sol optimize(Cons, obj, ops); % 结果提取 P_chp_v value(P_chp); H_chp_v value(H_chp); P_w_v value(P_w); E_s_v value(E_s); P_eb_v value(P_eb); abandoned P_wf - P_w_v; % 各时段弃风量这段代码跑通以后你会得到一组完整的调度结果。我建议随手保存成结构体方便后面做方案对比result struct(P_chp, P_chp_v, H_chp, H_chp_v, ... P_eb, P_eb_v, E_s, E_s_v(1:end-1), abandoned, abandoned); save(case_storage_eb.mat, result);3.4 求解之后先别急着画图先做三项自检拿到结果别急着画五彩斑斓的图先做三个自查能省掉很多返工时间。第一看电平衡和热平衡是否严格满足。Yalmip求解成功不代表你的约束写对了。把 value 结果代回等式两边算一遍残差如果某个时段差了0.1MW以上多半是单位或者时间对位出了问题。第二看储热罐SOC曲线是否符合物理直觉。在我们这个风电消纳场景里SOC曲线应该是白天爬升、夜间下降。如果你的结果是反的说明你把充放热的符号方向搞反了或者热平衡方程里充放热的位置写错了。我见过不止一个同学把 H_ch 和 H_dis 在热平衡里的位置写反结果模型很“聪明”地用储热罐解决了一切但实际根本不可行。第三看弃风时段是否集中在夜间。正常情况下弃风应该主要出现在凌晨2点到6点的低谷时段。如果白天中午的弃风比夜间还严重可能是储热罐容量设得太大或电锅炉功率设得过高导致白天为了充热把机组出力压得太高。做完这三项自检才算真正拿到了一个可用的结果。4. 算例对比加装储热罐和电锅炉后风电消纳提升了多少4.1 算例场景设置一台300MW热电机组加一座240MW风电场为了把前面那套东西讲透我构造了一个简化的典型场景数据是虚构的但量级参照北方冬季常见运行情况。系统包括一台300MW抽汽式热电机组纯凝最小技术出力90MW最大供热能力200MWth一座240MW风电场某日夜间预测出力达到240MW白天出力相对较小电负荷曲线白天峰值450MW夜间低谷280MW热负荷曲线白天低谷80MWth夜间高峰180MWth方案A什么都不加方案B仅加一台400MWth储热水罐充放热功率120MW方案C储热罐加一台80MW电锅炉我先把这组曲线排成三个工况然后直接用第3章的代码分别求解。4.2 结果对比消纳率从36%到82%的变化过程三种工况的运行结果对比如下工况全天弃风量(MWh)弃风率风电消纳率无储热、无电锅炉163.264%36%仅加储热罐116.846%54%储热罐电锅炉45.618%82%先说工况一。夜间热负荷180MWth机组最小电出力 90 0.35×180 153MW风电预测240MW而电负荷只有280MW所以风电最多只能上127MW剩下113MW全得弃。全天累计下来弃风率接近64%这个数字在我们实际调研中并不算夸张。工况二加了储热罐关键变化在夜间。储热罐以120MW的功率放热相当于替代了机组120MW的热出力机组热出力降到约60MWth最小电出力降到90 0.35×60 111MW比原来少了42MW这些空间全部让给了风电。同时储热罐在白天热负荷低谷时段充热回补SOC曲线走出了一个“白天爬升、夜间下降”的循环。看到这条曲线的时候我们项目组几个同事还专门讨论过——如果不看数据很多人会下意识认为储热罐应该是“夜间蓄热白天放热”但在风电消纳这个目标下逻辑正好反过来把机理理顺之后一切都顺了。工况三在储热罐的基础上加了80MW电锅炉。夜间电锅炉直接用风电制热一方面相当于增加了80MW本地负荷另一方面它产出的热又替代了部分机组供热进一步压低机组热出力。两重效果叠加夜间机组热出力可以压到很低的水平风电消纳率一下拉到82%。如果储热罐容量再加大、电锅炉功率再提高消纳率还能往上走但边际收益会明显递减这就是下一节要说的敏感性。4.3 敏感性分析储热罐容量和电锅炉功率不是越大越好我跑了一组不同参数的试验把储热罐容量从200一路加到800MWth结果很有意思200MWth的时候消纳率从36%提到48%加了200MW容量只带来12个百分点的提升加到400MWth消纳率到54%再加到800MWth消纳率只到58%左右。原因在于储热罐是“时间平移器”但它不能凭空制造热功率缺口夜间能够被替代的供热负荷是有限的到了上限以后再加容量只会让白天充热时间变长SOC曲线变得更平却不能再多消纳多少风电。电锅炉也类似。功率从40MW加到120MW消纳率提升越来越小因为夜间的风电总量有限电锅炉已经把可用的风电用掉了再往上加就只能靠煤电来给电锅炉供电经济上完全不划算。所以做方案设计的时候先跑一遍敏感性分析再定设备容量比拍脑袋选型靠谱得多。这也是联合优化控制作为设计工具的核心价值之一——不只是运行优化还能用来做投资决策。5. 这套代码落地时的几个坑预测误差、SOC维护与求解速度5.1 热负荷预测不准储热计划很容易“白做”上面的算例假设热负荷、电负荷、风电出力全都跟预测一致但实际运行不可能这么理想。我踩过最深的坑就是热负荷预测偏差。有一回我们把某日夜间热负荷预测成了140MWth优化结果让储热罐在夜间放热、机组低出力结果实际热负荷飙到了200MWth储热罐一个多小时就放空了机组被迫紧急加出力风电被大量削减最后实际效果比不用储热罐还差。后来我把静态优化改成滚动优化每个小时重新读取最新的热负荷实测值、未来时段的热负荷预测和风电预测以当前SOC为初值重新求解未来24小时的计划但只执行第一个小时。起始时间是滑动的模型结构一点不用变只需要把 E_s(1) 改成实测SOC、把负荷数组往前滚动。Matlab里用for循环就能实现代码量增加不到30行。这个改动虽然简单但实打实救了不少运行场景。5.2 储热罐SOC的“死区”问题目标函数里的小手脚第二个坑是SOC长时间贴在边界上不动。因为目标函数只看成本和弃风SOC只要在上限内就不产生任何收益或惩罚所以优化结果里经常出现储热罐已经蓄满却还在慢慢充、或者见底了还在放。这样的计划物理上可行但一点裕度都没有一旦实际负荷和预测差一点就崩了。我常用的办法是在目标函数里加一个很轻的SOC偏离惩罚项让模型倾向把SOC维持在30%到70%的中间区间obj ... beta * sum((E_s - E_ref).^2)其中 beta 取一个很小的值比如0.01不至于干扰主目标又足以让SOC曲线变得“平滑居中”。另一个等价做法是加约束 E_s(T) E_s(1)要求日末回到初始值这样储热罐就不容易跑偏。两种方法可以同时用效果更稳。5.3 求解时间连续LP秒出结果加了整数变量就“心梗”如果不考虑机组启停纯连续LP用Cplex秒解100个时段也不在话下。但一旦引入0/1启停变量哪怕只有一台机组24个时段模型也变成了MILP求解时间可能从零点几秒飙到几十秒甚至几分钟。这在做滚动优化的时候是完全没法接受的因为每个小时都要重算一遍。我的处理方式分三层第一层日计划优化用MILP允许它跑几秒到十几秒第二层进入滚动优化后把机组启停计划固定为当前已知的机组状态只优化连续变量这样每一轮都是一个LP几毫秒就出结果第三层如果偶尔需要重新校验启停方案就在整点时刻跑一次MILP其他时刻只跑LP。这套分层策略在工程上非常实用强烈建议按这个思路组织代码结构。5.4 代码移植的通用性参数与模型分离最后聊一个代码工程上的经验。我第一版代码把所有参数直接硬编码在脚本里换一个算例就得满文件找参数极其痛苦。后来我把数据读入、参数定义、模型构建、求解、后处理拆成了五个函数参数全部通过结构体传入。再换机组、换地区的数据时只需要改Excel或者改一个参数文件模型代码一行不动。如果你打算把这套代码扩展成自己的毕业设计或者项目工具这一点越早做越省事。具体拆法很简单一个 main.m 负责读数据和调度函数一个 build_model.m 负责返回约束和目标函数solve_model.m 调求解器post_process.m 画图并输出报表。每个函数都短小清晰debug的时候定位问题也很快。我用这套框架跑过不同容量的机组、不同季节的负荷曲线换数据基本在半小时内搞定。按这套思路把模型跑通之后你手里其实就多了一个非常灵活的工具——既可以做日前的运行计划也可以改造成滚动优化做实时调度还能拿来做储热罐和电锅炉的容量规划。关键在于先把“热”和“电”的耦合关系吃透把优化模型的物理逻辑掰扯清楚剩下的Matlab代码只是把你脑子里的那套规则翻译给求解器听而已。