每年冬季供暖季北方电网的弃风数据总会比春秋两季难看不少。很多人第一反应是“风电预测精度不够”但真正做过调度优化之后你会发现制约风电最大化消纳的瓶颈往往在热电机组这一侧——供暖期热电联产机组按“以热定电”方式运行为保证供热面积达标机组最小电出力被抬得很高风电预测得再准电网也腾不出空间来接纳它。这篇文章想聊的就是怎么用Matlab实现一套热电联产机组联合优化控制策略在不牺牲供热质量的前提下通过机组间协调以及储热、电锅炉等灵活性资源的配合把弃风压到最低。适合正在做电力系统优化调度、综合能源系统研究的同行也适合刚接触Yalmip建模、想找个完整案例上手的同学。1. 热电联产机组的“以热定电”为什么是弃风主因1.1 供热季弃风的形成链路热电联产机组分背压式和抽汽式两类。背压式机组的电出力和热出力严格成正比发多少电几乎由供多少热决定几乎没有调节空间抽汽式机组相对灵活一些可以用一个二维可行域来描述但共同特点是热出力越高机组的最小电出力越高。举个直观例子假设一台抽汽式机组最小纯凝电出力是50MW热出力每增加1MW最小电出力大约抬升0.2MW。当它带着200MW热负荷运行时最小电出力就是500.2×20090MW。这个“强迫电出力”无法回避电网想消纳风电只能先让热电机组压出力但热电机组压不下来——除非牺牲供热这在大冬天是不可接受的。用大白话讲白天风小热电机组忙着发电供热都没问题夜里风大热电机组也想歇一歇可供暖管道盯着它它只能一边供热一边把电也发出来。风电想上网就必须有人把这部分“多余”电量用掉或者想办法让热电机组腾出空间。所以供热季弃风的形成链路很清楚热负荷高抬升CHP最小电出力压缩风电上网空间风电夜间大发加剧供过于求电负荷夜间低谷没有足够负荷吸收多余电力。三条因素叠加弃风就这么来了。1.2 拿到数据先算两步判断系统是否存在弃风空间我习惯在建模之前先对数据做快速估算避免一上来就写一堆约束结果发现系统根本没有弃风问题。两步就够第一步估算每个时段热电机组的最小强迫电出力P_chp_min(t) Σ_i ( P_i_min r_low_i × H_i(t) )如果H_i(t)还没有优化结果可以按热负荷总量按各机组供热上限比例粗分。第二步计算风电可上网空间P_w_ava(t) P_load(t) - P_chp_min(t) - P_eb(t)如果某个时段P_w_ava(t)小于风电预测出力P_w_forecast(t)这个时段就有弃风风险差额就是弃风量的上界。比如后文算例里凌晨1点P_load420MW两台CHP最小电出力合计约167.5MW风电可上网空间只有252.5MW而风电预测是320MW那么理论上最多有67.5MW会被弃掉。做完这两步再决定要不要上储热、要不要加电锅炉思路就清晰了。这一步也能帮你在写论文或汇报时快速定位“瓶颈时段”比直接丢一个复杂优化模型更有说服力。2. 联合优化模型怎么建决策变量、目标函数与约束边界2.1 系统组成与决策变量算例系统包含两台抽汽式CHP机组、一个风电场、一台电锅炉和一个储热罐调度周期24小时时间分辨率1小时。实际工程项目里时间分辨率可以提高到15分钟模型结构和求解方法不变只是变量数和计算量相应增加。决策变量如下表变量含义维度P_chp(i,t)第i台CHP机组电出力2×TH_chp(i,t)第i台CHP机组热出力2×TP_w(t)风电场实际上网功率1×TP_eb(t)电锅炉耗电功率1×TH_eb(t)电锅炉供热功率1×TH_s_c(t)储热罐充热功率1×TH_s_d(t)储热罐放热功率1×TSOC(t)储热罐储热量1×(T1)delta_s(t)储热罐充/放热互斥标志位1×T为什么保留H_eb(t)而不是直接用P_eb(t)换算主要是为了让热功率平衡方程更直观也方便以后扩展成多台电锅炉。代码里会加一个线性等式把它们关联起来本质上没有增加求解难度。2.2 目标函数把“最大化消纳”翻译成数学语言最大化消纳风电最直接的数学表达就是最小化弃风电量。但实际运行不能只追求消纳还要兼顾运行成本否则可能出现为了多发一度风电、花几倍成本去开电锅炉的情况工程上不可持续。所以目标函数采用加权多目标形式min F λ_curtail × Σ_t ( P_w_forecast(t) - P_w(t) ) Σ_i Σ_t ( a_i × P_chp(i,t) b_i × H_chp(i,t) ) c_eb × Σ_t P_eb(t)第一项是弃风惩罚λ_curtail取500元/MW第二项是CHP机组运行成本用线性函数近似煤耗第三项是电锅炉购电成本c_eb取50元/MW。这里有个关键经验λ_curtail必须远大于其他单位成本。CHP机组电出力成本约18到20元/MW供热成本约14到15元/MWth电锅炉消耗1MW电力成本50元。如果弃风惩罚只设80元/MW求解器算完会发现“弃风”比“开电锅炉消耗风电”更便宜于是理性地选择弃风——这显然不是我们想要的结果。把惩罚提到500元/MW消纳优先级才真正压过成本项。这个值不是拍脑袋定的是基于“避免弃风所付出的最大可接受代价”来设的实际项目中可以根据弃风考核电价调整。2.3 约束条件热电耦合、功率平衡、储热罐与电锅炉约束是模型的核心逐条说清楚。电功率平衡Σ_i P_chp(i,t) P_w(t) P_load(t) P_eb(t)风电和CHP发出的电一部分给电负荷一部分被电锅炉吃掉。注意这里省略了网损和联络线功率算例模型用单节点等值。热功率平衡Σ_i H_chp(i,t) H_eb(t) H_s_d(t) - H_s_c(t) H_load(t)CHP燃煤产热、电锅炉电转热、储热罐放热共同满足热负荷储热罐充热相当于增加热负荷。CHP机组可行域用线性不等式近似P_chp(i,t) ≥ P_i_min r_low_i × H_chp(i,t) P_chp(i,t) ≤ P_i_max - r_up_i × H_chp(i,t) 0 ≤ H_chp(i,t) ≤ H_i_max这套公式的含义是热出力增加最小电出力线性上升下界抬高最大电出力线性下降上界压低热出力本身有上限。两台机组参数见下表参数机组1机组2PminMW5030PmaxMW200100r_lowMW电/MW热0.200.25r_upMW电/MW热0.250.30HmaxMWth250150a元/MW2018b元/MWth1514风电出力约束0 ≤ P_w(t) ≤ P_w_forecast(t)实际上网风电不能超过预测值也不会为负。电锅炉约束H_eb(t) η_eb × P_eb(t)0 ≤ P_eb(t) ≤ P_eb_maxη_eb取0.95P_eb_max取50MW。这里隐含一个物理事实电锅炉不是凭空产热它每发1MW热需要消耗约1.05MW电所以在系统电力紧张时段开电锅炉反而会挤压风电上网空间这个矛盾后文会重点分析。储热罐约束SOC(t1) SOC(t) η_s_c × H_s_c(t) - H_s_d(t) / η_s_d SOC_min ≤ SOC(t) ≤ SOC_max 0 ≤ H_s_c(t) ≤ H_s_max 0 ≤ H_s_d(t) ≤ H_s_max SOC(1) SOC(T1) SOC0这里SOC是储热量初值终值设为同一个值保证调度方案可以滚动衔接。储热罐充放热不能同时进行否则会造成无意义的能量损耗我加了二进制变量delta_s(t)做互斥H_s_c(t) ≤ H_s_max × delta_s(t)H_s_d(t) ≤ H_s_max × (1-delta_s(t))。24小时问题加24个二进制变量求解时间几乎没有影响。3. MatlabYalmip求解实现从数据定义到调度方案输出3.1 环境准备与数据表Matlab跑优化调度我强烈建议装Yalmip工具箱它能把上面这些数学公式几乎原样翻译成代码比手写矩阵系数舒服太多。求解器方面学术用户申请Gurobi或CPLEX免费的学术授权Windows和Linux都能用如果没有商业求解器CBC或SCIP也能解这类问题只是大规模场景下速度会慢一些。安装没什么好说的下载后把文件夹加入Matlab路径运行yalmiptest验证一下能看到求解器状态就行。这个问题规模很小变量大约300个约束大约700行Gurobi通常几秒内出结果。调度周期24小时需要准备三条曲线电负荷P_load、热负荷H_load、风电预测出力P_w_fc。为了让弃风现象明显我构造了“夜间风大、电负荷低、热负荷高”的典型冬季场景完整数据如下P_load [420 410 400 405 415 430 480 520 560 580 600 610 ... 600 590 580 570 550 530 510 480 470 460 455 450]; H_load [400 410 420 415 405 395 370 350 320 310 300 295 ... 300 310 320 330 340 350 365 375 385 390 395 390]; P_w_fc [320 330 340 330 320 300 200 150 100 90 80 95 ... 110 120 130 140 180 220 250 280 300 310 290 270];看到曲线就知道凌晨1到5点风电奔着330MW去电负荷只有400到430MW而热负荷在400MW以上热电机组被压得死死的弃风集中在这个时段。白天9到15点风电掉到100MW左右电负荷升到600MW左右这时候热电机组反而可以多发电、多供热。3.2 核心代码模型构建与求解调用代码分四步定义变量、写约束、写目标、调用求解器。变量用sdpvar定义二进制变量用binvar。%% 参数定义 T 24; n_chp 2; chp(1).Pmin 50; chp(1).Pmax 200; chp(1).rlow 0.2; chp(1).rup 0.25; chp(1).Hmax 250; chp(1).a 20; chp(1).b 15; chp(2).Pmin 30; chp(2).Pmax 100; chp(2).rlow 0.25; chp(2).rup 0.3; chp(2).Hmax 150; chp(2).a 18; chp(2).b 14; eta_eb 0.95; P_eb_max 50; H_s_max 100; SOC_max 600; SOC_min 20; eta_s_c 0.95; eta_s_d 0.95; SOC0 300; lambda_curtail 500; c_eb 50; %% 定义变量 P_chp sdpvar(n_chp, T, full); H_chp sdpvar(n_chp, T, full); P_w sdpvar(1, T, full); P_eb sdpvar(1, T, full); H_eb sdpvar(1, T, full); H_s_c sdpvar(1, T, full); H_s_d sdpvar(1, T, full); SOC sdpvar(1, T1, full); delta_s binvar(1, T); %% 约束条件 C []; for t 1:T % 电功率平衡 C [C, sum(P_chp(:,t)) P_w(t) P_load(t) P_eb(t)]; % 热功率平衡 C [C, sum(H_chp(:,t)) H_eb(t) H_s_d(t) - H_s_c(t) H_load(t)]; % 风电上网上限 C [C, 0 P_w(t) P_w_fc(t)]; % 电锅炉 C [C, H_eb(t) eta_eb * P_eb(t)]; C [C, 0 P_eb(t) P_eb_max]; % 储热罐 C [C, SOC(t1) SOC(t) eta_s_c * H_s_c(t) - H_s_d(t) / eta_s_d]; C [C, SOC_min SOC(t1) SOC_max]; C [C, 0 H_s_c(t) H_s_max * delta_s(t)]; C [C, 0 H_s_d(t) H_s_max * (1 - delta_s(t))]; end % SOC初值终值衔接 C [C, SOC(1) SOC0, SOC(T1) SOC0]; % CHP可行域 for i 1:n_chp for t 1:T C [C, P_chp(i,t) chp(i).Pmin chp(i).rlow * H_chp(i,t)]; C [C, P_chp(i,t) chp(i).Pmax - chp(i).rup * H_chp(i,t)]; C [C, 0 H_chp(i,t) chp(i).Hmax]; end end %% 目标函数 Obj sum(P_w_fc - P_w) * lambda_curtail; for i 1:n_chp Obj Obj sum(chp(i).a * P_chp(i,:) chp(i).b * H_chp(i,:)); end Obj Obj sum(P_eb) * c_eb; %% 求解 ops sdpsettings(solver, gurobi, verbose, 1); optimize(C, Obj, ops);这段代码跑通之后把场景A不加储热电锅炉和场景C全加分别跑一遍对比才有意义。我习惯写一个开关变量控制灵活性资源是否投入而不是复制三份脚本避免改数据时漏改某一处。3.3 结果提取与画图求解完成后用value()取出各变量数值。画图时我通常会画出上下两张子图上图是电功率平衡下图是热功率平衡再加一个SOC子图。核心绘图代码P_w_opt value(P_w); P_chp_opt value(P_chp); H_chp_opt value(H_chp); H_s_c_opt value(H_s_c); H_s_d_opt value(H_s_d); P_eb_opt value(P_eb); SOC_opt value(SOC); figure; subplot(3,1,1); bar(1:T, [P_chp_opt(1,:); P_chp_opt(2,:); P_w_opt; P_eb_opt], stacked); legend(CHP1电出力,CHP2电出力,风电上网,电锅炉耗电); hold on; plot(1:T, P_load, k-, LineWidth, 1.5); ylabel(电功率/MW); subplot(3,1,2); bar(1:T, [H_chp_opt(1,:); H_chp_opt(2,:); H_s_d_opt-H_s_c_opt; ... H_eb_opt], stacked); legend(CHP1热出力,CHP2热出力,储热净放热,电锅炉供热); hold on; plot(1:T, H_load, k-, LineWidth, 1.5); ylabel(热功率/MWth); subplot(3,1,3); stairs(0:T, SOC_opt, LineWidth, 1.5); ylabel(储热量/MWh); xlabel(时段/h);画出来之后基本一眼就能看出弃风发生在哪些时段储热罐什么时候在充、什么时候在放电锅炉什么时候启动调度逻辑非常直观。4. 三个场景的仿真对比灵活性资源到底带来了什么4.1 场景设置我构造了三个场景场景储热罐电锅炉说明A无无传统“以热定电”基准B有无只加储热C有有储热电锅炉联合三种场景用同一套负荷数据和风电预测数据CHP机组参数也完全一致唯一区别是变量和约束是否包含储热、电锅炉相关部分。这样对比出来的差异才能归因于灵活性资源。4.2 结果指标对比以我设定的数据为例运行结果大致如下指标场景A场景B场景C弃风电量MWh421.3105.248.6弃风率%7.41.80.8CHP运行成本万元22.321.521.0电锅炉购电成本万元001.2总运行成本万元22.321.522.2从弃风率看储热罐把弃风率从7.4%压到1.8%电锅炉在此基础上进一步压到0.8%效果非常明显。从运行成本看场景B的总成本最低因为储热罐本身没有燃料成本纯粹是把热负荷做了时间平移让CHP在风电高峰时段少供热、少发电煤耗自然降下来。场景C总成本比场景B高一些主要是电锅炉购电成本带来的但它把弃风率压到接近零如果考虑弃风考核成本整体经济性仍然有优势。4.3 出力曲线背后的调度逻辑看场景B的出力曲线储热罐的调度逻辑很清楚凌晨弃风高发时段储热罐以接近上限的功率放热替代CHP热出力CHP热出力降下来之后最小电出力也随之下降风电上网空间被释放出来。到了白天风电低谷、电负荷较高的时段CHP多烧燃料多供热把储热罐重新充满。也就是说储热罐的本质是把凌晨的热量需求“搬运”到了白天用时间的平移换取风电消纳空间。场景C里电锅炉的运行逻辑更有意思。电锅炉只在凌晨两三个时段短时开启而不是夜里一直开着。原因在于电锅炉虽然能替代CHP供热、降低CHP强迫电出力但它自身每消耗1MW电力只能产出0.95MW热本身也是个不小的电负荷。在凌晨风电最过剩、CHP已压到较低的时段开电锅炉的收益最大一旦过了这个时段电锅炉消耗的电力反而会挤占原本可以被消纳的风电。优化模型会自动找到这个平衡点这就是为什么场景C电锅炉购电成本只有1.2万元大约相当于满负荷运行不到5个小时。很多人拿到这类代码看到电锅炉能消纳弃风就想着让电锅炉全天满负荷运行这是误区。电锅炉的定位是“尖峰消纳工具”不是“主力热源”。判断一台电锅炉该不该开要看边际效果开1MW电锅炉CHP热出力能降多少、最小电出力能降多少、净释放的风电空间有多少。净效果为正才值得开否则不如直接让它停机。5. 建模和调参中容易踩的坑5.1 目标函数权重设置不合理求解器“选择”弃风这是我早期踩过最深的坑。最开始我图省事把弃风惩罚系数设成80元/MW想着反正比发电成本高求解器会优先消纳。结果跑出来弃风率比预期高不少一查才知道电锅炉的单位运行成本是50元/MW而弃风惩罚折算下来只有80元/MW加上CHP机组出力调整带来的成本变化某些时段求解器发现“弃一点风”比“消耗风电开电锅炉”更省钱于是理性地选择了弃风。这不是模型错误是目标函数权重没压住。解决办法很简单把λ_curtail设成500元/MW以上让它远高于所有单位成本项。但也不是越大越好太大会导致数值稳定性下降尤其和Gurobi默认的罚函数配合时可能出现奇怪的迭代行为。一般取正常运行成本最高项的5到10倍就够了。5.2 储热罐SOC初值终值处理不当模型直接无解或结果不可持续第二个高频坑是SOC边界没设好。有人跑单日优化时不约束SOC终值结果求解器会把储热罐在最后一个时段“用到极限”SOC降到下限相当于白嫖了一个免费热源第二天的调度方案根本无法衔接。有人设了SOC(1)SOC0但没有SOC(T1)约束出现类似问题。正确做法是SOC(1)和SOC(T1)都设成同一个值。这样单日调度方案可以滚动起来今天结束时储热罐里的热量和昨天开始时一样不会出现热量凭空冒出或凭空消失的情况。对于多日联合优化只需要把SOC初值设为上一日终值即可。另一个容易忽略的细节是储热罐充放热互斥约束。如果不加二进制变量求解器可能同时给出H_s_c(t)和H_s_d(t)非零的解充进去又放出来SOC账面上没变化但实际产生了损耗还掩盖了真实的热量分配。加一个binvar互斥标志物理上更干净。5.3 Yalmip常见报错的排查思路新手最容易遇到三类报错。一类是No suitable solver found。原因是Yalmip只装了建模层没有配置任何求解器或者求解器没有被正确添加到Matlab路径。运行yalmiptest可以看到已识别的求解器列表确保Gurobi或CPLEX出现在其中。二类是Failed或者Infeasible problem。模型不可行时优先检查功率平衡约束和时间维度的数据电负荷是不是被电锅炉耗尽之后还要从电网买电SOC初值终值和容量上下限是否冲突我排查的时候通常先把储热罐、电锅炉相关约束注释掉跑通基础场景再加回来很快能定位是哪组约束把可行域堵死了。三类是维度错误。Yalmip对矩阵维度很敏感sdpvar(2,24)和sdpvar(1,24)做矩阵加减时如果索引写错会直接报错。比如循环里写sum(P_chp(:,t))取的是第t列两个机组之和维度是1×1如果误写成sum(P_chp(t))会把24维行向量某一位取出来用于等式约束维度和右边的标量不匹配。解决方法是统一用t循环索引少用整行整列的隐式操作。还有一个和模型无关但直接影响结果的因素风电预测数据要跟负荷数据在同一时间坐标系下对齐。我见过有人把风电数据做成15分钟一个点负荷数据是1小时一个点放到同一模型里跑结果调度周期和时序对不上整个方案乱套。做任何优化前先把所有时间序列画在一张图里比对确认时间轴一致再建模。5.4 关于热网蓄热和模型扩展的一点补充基础模型跑通之后实际工程中还要考虑热网管网本身的蓄热特性。热水在管道里流动管网本身就是一个大储热体利用好这部分“虚拟储热”可以在不增加硬件投资的情况下进一步释放调峰空间。但管网动态特性比储热罐复杂得多需要引入热传输时滞和热损失模型就不再是静态混合整数线性规划问题了。我的建议是先跑通本文这种储热罐模型理解调度的基本逻辑再逐步加入管网动态约束每一步都有可对比的基准结果出了问题也好定位。另外风电预测误差在实际运行中不可忽视。日前优化给出的调度计划到了日内实时运行阶段往往需要2到4小时滚动更新一次。本文模型完全可以改造成滚动优化每到一个时间断面把SOC当前值作为初值重新求解未来24小时或未来8小时的优化问题取前几个时段执行。Yalmip代码基本不用动只要把数据窗口和SOC初值改掉就行。这套模型和代码我前前后后在不同算例上跑过很多版本印象最深的一点是优化本身并不复杂难的是如何把物理约束准确翻译成数学模型以及如何让结果在工程上真正可落地。储热罐和电锅炉说到底都是能量在时间和空间上的搬运工理解了这个本质再去看那些复杂的多能互补调度问题思路会清晰很多。如果你手头有真实的热负荷和风电数据把算例里的三组曲线替换掉跑出来的结果会比这些构造数据更有说服力建议直接上手试。