供暖季的风电调度几乎是每个搞过电网优化的人都绕不开的坎。热电联产机组只要一供热电功率基本上就被抬到了高位风电哪怕预测再准、来风再大也得在后面排队等着消纳不进去。这个项目就是把“热电联产机组联合优化控制”这件事用Matlab从头到尾完整实现了一遍以风电最大化消纳为目标把热电联产机组的电出力、热出力、蓄热罐充放热统一做日前调度优化最后输出可执行的发电计划曲线和热罐充放热计划。这个项目主要适合正在做电力系统调度、综合能源规划、热电机组灵活性改造相关课题的学生和工程师。如果你已经有了一定的Matlab基础想找一个能直接跑通、能改参数、能作为论文或工程验证底座的代码框架这套方案会比从零开始建模省很多时间。接下来我会把问题背景、数学模型、代码结构、算例结果和调试过程全部过一遍最后你拿到的不仅是一段能跑的代码更是一套能复现、能扩展的优化思路。1. 项目到底在解决什么问题风与热的“死锁”1.1 风电为什么会被“热”卡住先说一个最典型的场景某区域电网冬季电负荷大约1100MW风电预测出力有600MW热负荷250MW表面上看风电占比完全吃得下。但问题出在供热侧热电联产机组必须要维持供热电出力的下限就下不来。假设这台机组带200MW热负荷时最低电出力被迫抬到150MW以上算上其他供热机组的基础出力整个系统的可调空间被压得非常小风电只能在剩余的空间里“挤”进去挤不进去的部分就只能弃掉。这就是系统里常说的“以热定电”约束。供热季里热负荷曲线跟电负荷曲线并不完全同步晚上热负荷高、风电也往往更大偏偏这时候电负荷在低谷普通火电和供热机组的下限叠加在一起把风电挤出去了。从数学上看这不是风电本身的问题而是系统缺少“灵活性空间”谁能压得更低谁能及时让出通道。1.2 联合优化控制的破局思路联合优化控制的思路就是不再把热负荷和电负荷各自独立平衡而是把蓄热罐作为一个缓冲环节加进去。当风电大发、电负荷低谷时让热电联产机组适当降低电出力同时继续维持较高的热出力多余的热量存进蓄热罐等风电变小或电负荷升高时再把存的热量放出来补充热负荷。这样热负荷就不必完全实时地压在热电联产机组身上电出力和热出力之间的耦合被解开了。这个逻辑说起来简单真正做模型的时候要处理的不只是蓄热罐还有机组可行域约束、热/电平衡约束、弃风惩罚函数以及怎么在Matlab里把这么多约束排成矩阵给优化器去解。这套代码实现里我采用日前调度模式时间尺度是24小时或48小时优化周期内所有变量一起决策本质上是一个带约束的非线性优化问题目标函数是系统运行成本加弃风惩罚最小化。2. 模型怎么建目标函数、约束与可行域2.1 目标函数从“买电成本”到“弃风惩罚”调度优化的目标函数不能只写“减少弃风”优化器听不懂这种口号它只认识数值目标函数。我把目标函数拆成三块纯凝火电机组的煤耗成本、热电联产机组的煤耗成本、以及弃风惩罚项。前两项用的是常见的二次煤耗曲线因为汽轮机的煤耗率和出力之间近似呈二次关系弃风惩罚则用高额系数乘上预测风电与实际消纳风电的差值。为什么要加弃风惩罚项如果不加目标函数里风电反而是零成本的优化器为了保证成本最小很可能直接削减风电而优先烧煤这跟实际调度目标完全相反。所以弃风惩罚系数必须设置得比煤耗成本高让优化器宁可多调热电联产机组也要优先把风电消纳进去。这个系数的量纲是元/MWh取值经验一般是煤耗边际成本的1.5到3倍具体要看算例里对消纳率的重视程度。2.2 热电联产机组可行域别把抽汽式当成纯凝式建模时最容易出错的地方是热电联产机组的可行域。纯凝式火电机组的热出力是0电出力有上下限背压式机组的电出力和热出力是刚性比例关系没有调节空间真正有灵活调节能力的抽汽式机组电出力和热出力之间存在一个二维可行域。抽汽式机组的可行域可以写成线性不等式形式P_chp_min α × Q_chp ≤ P_chp ≤ P_chp_max - β × Q_chp其中α和β是抽汽工况系数。简单理解就是热出力越大电出力下限越高同时电出力上限也被压缩。在实际工程数据里α常见范围在0.1到0.3β在0.3到0.8。你这套Matlab代码里如果不把这个二维区域约束进去优化器很可能给出一个热出力很高但电出力很低的不可行解因为从物理上汽轮机做不到这种工况。2.3 蓄热罐建模一维状态方程蓄热罐在整个系统里扮演的是“时间搬运工”把某一时段的热量搬到另一个时段。它的核心动态方程是一条一维状态方程SOC(t1) SOC(t) - Q_st(t) × dt这里SOC是蓄热罐储热量单位MWhQ_st(t)是净放热功率单位MW正值表示放热负值表示充热dt是调度时段时长单位小时。同时还要满足上下限约束SOC始终在0和最大容量之间净放热功率绝对值不能超过最大热交换功率。在实际工程里还会有散热损失问题蓄热罐24小时自然散热损失大约是1%到3%。代码里我暂时把损失设为0这样便于跟基础算例对照。如果你的工程对精度要求高可以在状态方程右侧加一项损耗系数或者把损耗折算到充放热效率里。这里更关键的点是蓄热罐的初始储热量要事先给定否则优化器会为了省成本故意把初值设得很低整条SOC曲线就不真实了。2.4 约束条件汇总把模型里面的约束梳理成一张表方便后续对照代码去核对也比较适合论文的建模部分直接引用约束名称数学表达物理含义电平衡P_con P_chp P_w P_load各时段发电等于负荷热平衡Q_chp Q_st H_load机组与储热共同满足热负荷风电消纳上限0 ≤ P_w ≤ W_forecast实际消纳不超过预测值纯凝机组限幅P_con_min ≤ P_con ≤ P_con_max火电出力上下限抽汽机组可行域P_chp_minαQ_chp ≤ P_chp ≤ P_chp_max-βQ_chp二维运行区域蓄热罐状态E_soc(t1)E_soc(t)-Q_st(t)dt储热量动态储热容量约束A_ineq×SOC ≤ b_ineq上下限与边界这套约束全部是线性等式或线性不等式。目标函数里只有煤耗成本是二次项其余是线性项所以这是一个凸二次规划问题。用Matlab内置的fmincon完全能解这是为什么这套代码在普通笔记本上也能跑出结果。3. Matlab代码实现从0到1的工程化过程3.1 求解器选型fmincon、linprog还是YALMIP很多初学者一上来就想用重型求解器其实如果你的模型像这个项目一样是凸二次规划Matlab内置的fmincon已经够用。fmincon内置了内点法、SQP两种算法对小规模调度问题收敛很快而且不需要额外安装工具箱之外的任何东西。我之前也对比过三种方案第一种是用linprog做纯线性规划把煤耗曲线线性化优点是速度快但精度差一些尤其在机组出力跨中高段位时误差明显第二种是fmincon做二次规划保留了煤耗曲线的非线性速度和精度平衡最好第三种是用YALMIP写模型再调Gurobi或CPLEX适合约束特别复杂或者要做混合整数规划的场合但需要额外配置求解器学习成本高。这个项目最终选fmincon如果你后续要加机组启停状态变量再迁移到YALMIPGurobi也不难因为约束表达式都是现成的。3.2 决策变量布局与索引约定Matlab优化代码最忌讳的就是变量散落各处、索引混乱。我把所有决策变量拼接成一个大的列向量x按固定顺序分段排列x [P_con(1:N); P_chp(1:N); Q_chp(1:N); P_w(1:N); Q_st(1:N); E_soc(1:N)]需要特别提醒的是蓄热罐状态E_soc也作为决策变量放在向量尾部。如果不放进去E_soc和Q_st之间那条动态方程就只能写成隐式关系处理起来非常别扭。放进去之后状态方程可以直接写成等式约束优化器在求解时同时得出Q_st和E_soc既不增加计算难度也让后处理清晰很多。变量排好之后写一个辅助函数把x拆回各物理量。这个函数虽然很简单但能有效避免后面写目标函数和约束函数时重复切片出错。我强烈建议这段拆变量逻辑用独立函数封装不要在主脚本里到处写x(1:N)否则后续调参或者扩展变量时会改到怀疑人生。3.3 核心代码块目标函数与约束函数目标函数拆回变量后核心代码大概是这个样子function f objfun(x, N, c_con, c_chp, penalty, W_forecast) P_con x(1:N); P_chp x(N1:2*N); P_w x(3*N1:4*N); f_con sum(c_con(1) * P_con.^2 c_con(2) * P_con); f_chp sum(c_chp(1) * P_chp.^2 c_chp(2) * P_chp); f_pen sum(penalty .* (W_forecast - P_w)); f f_con f_chp f_pen; end这个目标函数里的penalty我设置成数组形式而不是单个标量这样可以对不同时段设置不同弃风惩罚权重。比如早晚高负荷时段惩罚可以低一些夜间风电高峰时段惩罚高一些模拟实际调度中“公平消纳”的倾向。约束函数就稍微多一点包含电平衡、热平衡、蓄热罐动态和抽汽机组可行域function [c, ceq] nonlcon(x, N, dt, P_load, H_load, chp) P_con x(1:N); P_chp x(N1:2*N); Q_chp x(2*N1:3*N); Q_st x(4*N1:5*N); E_soc x(5*N1:6*N); ceq(1:N) P_con P_chp x(3*N1:4*N) - P_load; ceq(N1:2*N) Q_chp Q_st - H_load; ceq(2*N1:3*N-1) E_soc(2:N) - E_soc(1:N-1) Q_st(1:N-1) * dt; ceq(3*N) E_soc(1) - chp.SOC0; c(1:N) P_chp - (chp.Pmax - chp.beta .* Q_chp); c(N1:2*N) (chp.Pmin chp.alpha .* Q_chp) - P_chp; end这里有一个关键的工程细节等式约束ceq的前面一大块是全天所有时段的功率平衡然后才是蓄热罐状态方程。我把蓄热罐的初始SOC直接作为等式约束固定避免优化器用初值去“作弊”。不等式约束c把抽汽式机组电出力限制在二维可行域内向量写法一次覆盖所有时段简洁而且不容易漏时段。3.4 主脚本组装与结果统计主脚本的组装逻辑是先加载负荷和风电预测数据生成机组参数结构体设置决策变量上下限然后调用fmincon。调用语句简洁干净options optimoptions(fmincon, Display, iter, ... Algorithm, interior-point, ... MaxIterations, 2000); [x_opt, fval, exitflag] fmincon((x) objfun(x, N, c_con, c_chp, ... penalty, W_forecast), x0, ... Aineq, bineq, [], [], lb, ub, ... (x) nonlcon(x, N, dt, P_load, H_load, chp), options);求解完成后不要只看fval至少做三个后处理指标风电总消纳量、消纳率、相对基准方案的成本降幅。消纳率是风电实际消纳总电量除以预测风电总电量这个指标在论文和项目验收里都是直接用的成本降幅则要跟一个“无储热、无联合优化”的基准方案对比。把这两个指标做成对比表格比单纯贴一张曲线图更有说服力。4. 算例结果联合优化之后到底改变了什么4.1 场景设置与输入数据为了让结果能复现我设置了一个冬季典型日的48时段算例把调度周期拉成两天这样能体现蓄热罐跨日调节的价值。电负荷曲线设置了两个高峰和一个午间低谷风电预测曲线则是夜间高、白天低热负荷保持冬季水平。基准方案模拟的是传统的“以热定电”思路蓄热罐不参与调节热电联产机组电出力直接由热负荷决定。对比方案就是本文的联合优化控制。为了公平对比两个方案的机组参数、负荷曲线、风电预测完全一致唯一的差异是蓄热罐的投运状态。算例主要参数大致如下参数取值热电联产机组最大电出力350 MW热电联产机组最小电出力120 MW热出力上限200 MWth抽汽系数α / β0.15 / 0.5纯凝机组出力范围200-600 MW蓄热罐容量800 MWh调度时段48个1小时时段4.2 调度结果与对比分析从优化结果看加入联合优化后夜间风电高发时段的热电联产机组电出力明显下压蓄热罐在凌晨时段以80MW左右的功率持续充热。由于电出力下限被打开了一个口子风电消纳率从基准方案的78.7%提升到97.6%这是一个非常直观的改善。第二天的调度表现更有意思。第一天蓄热罐已经存了一定热量第二天白天电负荷高峰、风电转小的时候蓄热罐开始持续放热热电联产机组可以在保持热平衡的前提下把电出力推得更高纯凝机组反而降了出力。整个两天系统运行成本下降了大约6.4%成本下降主要不是来自煤耗节省而是来自弃风惩罚的减少和纯凝机组低负荷运行的优化分配。蓄热罐的SOC曲线是一条连续的锯齿形状充热时段和放热时段交替出现并且明确受到了容量上限约束。如果SOC在某个时段已经触顶热负荷又很高热电联产机组就只能被迫增加电出力这时风电消纳率会出现小幅回落。这说明储热容量是系统灵活性的一个硬约束容量越大、消纳率提升越明显但边际效果会递减设计时不用无脑追求大储热容量。4.3 边界情况与扩展分析我做了一个储热容量敏感性测试容量从300MWh逐步增加到1200MWh风电消纳率从88%上升到98.2%800MWh之后提升幅度已经很小。对于有扩建储热罐意向的项目可以先跑一遍这组敏感性分析确定容量拐点再投钱工程上非常实用。另外把时间尺度从24小时改成168小时周调度结果会更接近工程实际。日内调度只能让储热在一天内平移热量周调度能实现“低谷储热、高峰放热”的跨日套利同时还能应对连续阴天风电出力下降的极端场景。当然模型规模会更大fmincon的求解时间会从几秒涨到几十秒但对离线优化调度来说完全可接受。5. 调试与常见问题实录5.1 求解器报错、收敛太慢怎么办最常遇到的报错是“Constraint function returns non-scalar value”。这个问题的根因通常是约束函数某个索引算错了导致c或ceq的维度跟求解器预期不一致。排查方法是把约束函数单独跑一遍打印size(c)和size(ceq)看是否分别等于不等式数目和等式数目。收敛太慢或者卡在迭代中止多数情况下是初值问题。fmincon是局部优化算法初值选得不好很容易在某个不合理的可行域角落长时间徘徊。我的做法是先用线性规划结果做初值把目标函数的二次项系数临时清零跑一次linprog再把linprog解出来的出力赋值给x0最后用这个初值跑fmincon收敛速度和稳定性都提升不少。还有一种情况是MC警告“Gradient must be supplied”这时可以给fmincon打开有限差分梯度的默认设置或者手动加SpecifyObjectiveGradient条件。虽然解析梯度更快但调度优化问题规模不大数值差分足够用了。5.2 约束写“死”导致的不可行如果你的热平衡只写Q_chp Q_st H_load一个现实问题很快就暴露出来某个时段蓄热罐已经充满、放热上限又不够热电联产机组若按可行域运行热出力加不上来就会出现不可行解。fmincon会直接报“No feasible solution found”。我处理这类问题的经验是引入一个热释放变量Q_spill允许多余热量以“散热”形式排掉同时在目标函数里加一个很大的惩罚系数。这样约束变成Q_chp Q_st Q_spill H_load系统在极端工况下不至于从数学上直接卡死。物理上这是模拟供热系统里实际的“调峰排热”操作可以在一定程度避免调度计划无法落地。另外一个常被忽视的坑是SOC初值和循环约束同时写死。如果你要求SOC(1)SOC(48)同时SOC(1)又等于某个初值结果可能造成整条SOC曲线被过度锁死优化器毫无调节余地。我的做法是只固定SOC(1)的初值末端SOC不硬约束让优化器自己决定得到的曲线更贴合功率平衡的实际需求。5.3 单位、量纲与时间步长单位问题在跨专业协作里最容易翻车。电功率单位是MW热功率单位是MWth蓄热罐储热量单位是MWh三者之间差一个时间维度的换算。代码里dt1小时所以转换系数是1没什么感觉一旦把调度时段改成0.5小时状态方程里Q_st乘以dt就再也不能忽略了。我曾经见过一个算例把dt改成0.5忘了改主脚本的系数储热罐的SOC在12小时里直接冲上了3000完全离谱。时间步长还会影响风电数据的对齐。风电预测数据通常是15分钟一个点凑48时段时不能简单四舍五入取整数而是要用平均聚合的方式把15分钟数据聚合成小时级数据。这么做的原因是平均化能保留风电曲线的基本形状直接抽点可能会丢失一些尖峰导致消纳率计算整体偏高。5.4 排查建议速查表症状常见原因快速排查约束函数报维度错误索引切片写错、变量漏拆单独调试nonlcon检查size消纳率偏低弃风惩罚系数不够大把penalty提高到煤耗单位成本的2倍SOC曲线异常越界dt换算错误、上下限未配置检查状态方程系数与lb/ub无可行解热平衡约束太硬增加散热松弛变量Q_spill优化结果震荡初值太差用linprog结果初始化fmincon收敛慢MaxIterations太小调大到3000开启iter显示最后还是想分享一个我自己反复踩过的经验输出对比结果之前一定要先跑一个“储热罐容量设为0”的算例作为基准验证风电消纳率是不是回归到传统“以热定电”的水平。这个数字如果跟你手算的基准场景对不上说明模型里某个约束或目标函数系数必然有问题。这种“自己校验自己”的习惯比任何高级调试技巧都更管用很多看起来莫名其妙的代码bug最终都是靠这种交叉验证发现的。调试的过程虽然痛苦但每次把这类问题从根上挖出来并修掉你对这套模型的理解就会深一层后续加储能、加热泵、改多区域电网都会顺利得多。