加过班的都知道热电联供微网这个方向听起来像是把几个设备包进一个优化模型里跑一跑的事但真上手做Matlab代码实现的时候几乎所有人都会在同一个地方卡住——不是不会写约束而是搞不明白电热耦合到底在调度里扮演什么角色。我最早做综合能源系统项目时也是这个状态第一版代码跑出来CHP机组在电价低谷时段疯狂满发蓄热罐纹丝不动夜间热负荷明明很高燃气锅炉却在加班。看到那个结果我才反应过来这个问题的核心从来不是会不会建模而是你有没有把多能互补的物理逻辑写进模型里。这篇博文就围绕基于多能互补的热电联供型微网优化运行展开。我会把CHP机组怎么建模、蓄热罐和电储能如何解耦电热约束、目标函数和约束集怎么设计、MatlabYALMIP这条代码链路怎么搭以及我实测中踩过的坑全部摊开来讲。适合刚接触微网优化方向的研究生也适合做园区综合能源设计的工程师拿来当参考。1. 电热耦合热电联供微网优化真正难的起点1.1 为什么以热定电是绕不开的核心约束先把场景搭起来。一个典型的热电联供型微网里通常包含这么几类设备CHP机组燃气轮机或内燃机加余热回收装置它同时产电和产热燃气锅炉只产热光伏和风电只产电但出力不可控电储能和蓄热罐分别负责储存电和热再加上与上级电网的购售电交互。光看设备清单你会觉得这不就是个混合整数线性规划问题嘛把所有设备建模丢给求解器就完事了。但真正的难点藏在CHP机组那一栏。CHP机组最特殊的地方在于它的电出力和热出力不是两个独立变量。很多教材和简化论文会用热电比来描述这种关系——热出力等于热电比乘以电出力。这个近似在窄负载范围内勉强能用但真实机组的可行域是一个多边形区域而不是一条直线。原因不复杂燃气轮机在不同负载率下余热回收系统的效率会变抽汽式汽轮机的供热抽汽量会直接影响凝汽发电的效率和出力上限。所以工程上更可靠的建模方式是可行域法——用一组线性不等式把电出力与热出力的所有可行组合圈出来。有了这个耦合关系以热定电就成了一个非常现实的调度难题。电负荷要求实时平衡热负荷也有自己的曲线两者往往不同步。假设某个冬夜热负荷很高CHP机组被迫高发有功电出力跟着一起抬高可此时的电网电价正处于低谷更经济的做法本来是少发电、多买电但热负荷硬逼着CHP必须开着这就是以热定电约束在起作用——明明买电更便宜却买不了因为CHP得产热。反过来白天热负荷低的时候CHP又得压低出力可能连电负荷都带不动只能高价从电网买电。这个耦合关系就是热电联供微网优化区别于普通电力调度的根本所在。普通微网调度里发电机想开就开、想停就停而在热电联供微网里每一个让CHP出力多少的决策都同时绑定着电力和热力两本账。1.2 多能互补的物理实质用时间维度上的自由度去破解耦合既然电热天生耦合那多能互补到底补的是什么我的理解是它补的不是能量总量而是调度自由度。第一种互补是设备层面的。当CHP机组被热负荷绑定时电储能可以削峰填谷光伏风电补充电力缺口当CHP被电负荷绑定时——比如电价高峰需要多发电——蓄热罐可以填补热缺口让CHP不必完全跟着热负荷走。储能设备在这里的本质作用是把电热同时平衡这一个硬约束松弛成两个可以被时间平移的软平衡。第二种互补是时间层面的。热负荷本身有惯性蓄热罐可以把后半夜的多余热量储存起来供白天热负荷高峰使用电储能则可以把午间光伏大发时的电量存下来留到晚高峰放电。两种储能介质的充放时段不同配合起来相当于给调度问题增加了多个时间维度的决策变量。第三种互补是价格信号层面的。分时电价本身就是一种很自然的经济引导电价高时多让CHP发电电价低时多买电、让蓄热罐蓄热。优化算法做的事情本质上就是在这些价格信号和设备自由度之间找最优平衡点。所以这一章最关键的一句话是多能互补建模不是把设备清单堆进约束里就行而是要把每一种设备提供的自由度建模清楚。设备越多约束矩阵越大解空间也越大最终调度结果的好坏完全取决于你有没有把这些自由度正确地写进模型。2. 数学模型搭建从设备特性到完整优化方程2.1 设备建模的工程取舍线性化是最稳的起点很多论文会把设备模型写得很细CHP机组用非线性效率曲线、燃气锅炉用分段线性模型、储能考虑充放电效率随SOC变化。但说实话对一个24时段、时间间隔1小时的日前调度问题来说线性化模型完全够用而且线性模型在求解稳定性上远优于非线性模型。原因有两点一是日前调度的精度要求还没高到需要逐点非线性效率曲线的程度二是求解器对MILP混合整数线性规划的处理能力远比对MINLP成熟。真要做高精度仿真那是另外一套时域仿真工具的事不是优化调度该干的活。我的建议是CHP用可行域线性不等式组建模燃气锅炉用恒定效率模型电储能和蓄热罐用一阶递推方程加充放互斥约束。下面逐个说关键建模方式。CHP可行域用一组不等式描述电出力P_chp和热出力H_chp的可行组合。以典型的小型燃气轮机CHP为例可行域可以写成H_chp_min H_chp H_chp_max; P_chp_min P_chp P_chp_max; P_chp a1*H_chp b1; P_chp a2*H_chp b2;其中a1、a2来自机组性能曲线的包络线拟合。注意环境温度和大气压对出力上限的影响可以折算进边界系数里修正。燃气锅炉相对简单H_boiler eta_boiler * F_boiler; 0 F_boiler F_boiler_max;这里F_boiler是燃料输入功率单位要统一成kWeta_boiler一般取0.85到0.92。千万别把天然气体积单位m3/h直接当成kW塞进公式我在后面会专门讲这个坑。电储能用一阶SOC递推方程SOC(t1) SOC(t) eta_ch * P_ch(t) * dt / E_cap - P_disch(t) * dt / (eta_disch * E_cap);注意这个公式里充电效率乘在充电项上放电效率除在放电项上很多人第一次写会搞反。同时要加充放互斥约束推荐用两个0-1变量线性化不要写P_ch * P_disch 0这种非线性约束。蓄热罐的建模和电储能完全同构只是把电量换成热量E_tank(t1) E_tank(t) H_ch_tank(t) * dt - H_disch_tank(t) * dt; E_tank_min E_tank E_tank_max;2.2 目标函数成本最小到底在算哪些账目标函数一般包括四块购电成本、燃料成本、运维成本以及可选的弃风弃光惩罚项。数学形式如下min sum_t ( C_buy(t)*P_buy(t) - C_sell(t)*P_sell(t) C_gas*(F_chp(t) F_boiler(t)) C_om_chp*P_chp(t) C_om_boiler*P_boiler(t) C_om_ess*(P_ch(t)P_disch(t)) )这里有两个容易忽略的细节。第一上网电价C_sell通常显著低于购电价C_buy所以模型大多数情况下会倾向于自产自用而不是卖电这符合实际运营逻辑。第二CHP的燃料消耗F_chp是电和热的共同输入严谨的建模要把它作为独立变量用线性关系约束住像F_chp c1P_chp c2H_chp c3这样。如果图省事直接写成P_chp/eta_e等于忽略了热电联产中热效率对燃料成本的贡献算出来的成本会虚高调度结果也会偏向锅炉而非CHP。2.3 约束条件全景一张表说清楚完整的约束集合可以用下面这张表梳理每个约束的含义和方程式一目了然约束类型数学形式说明电功率平衡P_chp P_pv P_wind P_disch P_buy - P_ch - P_sell P_load每个时段成立热功率平衡H_chp H_boiler H_disch_tank - H_ch_tank H_load每个时段成立CHP可行域一组线性不等式描述电热耦合关系锅炉上下限0 F_boiler F_boiler_max燃料输入受限电储能SOC递推SOC(t1)SOC(t)eta_chP_chdt/E_cap - P_dischdt/(eta_dischE_cap)时段间耦合电储能充放互斥P_ch Mu_ch, P_disch M(1-u_ch)防止同时充放蓄热罐递推E_tank(t1)E_tank(t)H_ch_tankdt-H_disch_tankdt时段间耦合购售电互斥P_buy Mu_buy, P_sell M(1-u_buy)同一时段不同时买卖旋转备用(可选)P_chp_max - P_chp R(t)增强可靠性特别说明一点热功率平衡用的是功率而不是热量因为蓄热罐的蓄放热本质是热功率随时间平移展开到24时段就是一条时间序列用功率平衡加罐内容量约束调度结果直接就是每个时段的功率指令工程上可以直接下发给现场控制系统。3. Matlab代码实现YALMIP建模链路与求解器选型3.1 为什么选YALMIPCPLEX/Gurobi这套组合这个项目我用的Matlab版本搭配YALMIP加Gurobi。选择理由很实际YALMIP把sdpvar变量、约束和目标函数封装得非常干净写模型的手感和数学公式几乎一一对应调试时能直接对照公式找问题CPLEX和Gurobi的MILP求解速度比Matlab自带的intlinprog快一个数量级——24时段含二进制变量的模型商业求解器秒级解完intlinprog经常要几十秒甚至更久。另外之后如果想从确定性调度换成鲁棒优化或随机优化YALMIP的框架不需要推翻重写只需增加场景集或不确定性集合。如果你暂时没有商业求解器的license用intlinprog也能跑通小规模算例但注意要在sdpsettings里指定solver,intlinprog。我实测下来对24时段、30个左右变量的模型intlinprog尚可接受一旦加入蓄热罐和储能多个时段的递推耦合变量性能和稳定性差距就很明显了。经费允许的话建议直接装一个Gurobi学术版免费且安装省心。3.2 代码框架从数据初始化到结果输出完整代码的结构大致是数据初始化、变量定义、约束组装、目标函数、求解、结果输出。我贴一个核心骨架细节参数用注释说明%% 1. 数据初始化 T 24; dt 1; % 时段数和时间间隔(小时) % P_load, H_load, P_pv, P_wind, C_buy, C_sell 均为 1xT 向量 % CHP可行域系数 a1,b1,a2,b2; 电储能参数 eta_ch, eta_disch, E_cap %% 2. 定义决策变量 P_chp sdpvar(1, T); % CHP电出力 H_chp sdpvar(1, T); % CHP热出力 F_chp sdpvar(1, T); % CHP燃料输入功率 P_boiler sdpvar(1, T); % 锅炉燃料输入 H_boiler sdpvar(1, T); % 锅炉热出力 P_buy sdpvar(1, T); % 购电功率 P_sell sdpvar(1, T); % 售电功率 P_ch sdpvar(1, T); % 电储能充电 P_disch sdpvar(1, T); % 电储能放电 SOC sdpvar(1, T1); % 电储能SOC注意T1个点 H_ch_tank sdpvar(1, T); % 蓄热罐蓄热功率 H_disch_tank sdpvar(1, T); % 蓄热罐放热功率 E_tank sdpvar(1, T1); % 蓄热罐储热量注意T1个点 u_ch binvar(1, T); % 充电0-1变量 u_disch binvar(1, T); % 放电0-1变量 u_tch binvar(1, T); % 蓄热0-1变量 u_tdisch binvar(1, T); % 放热0-1变量 %% 3. 约束组装 Constraints []; % 电功率平衡 Constraints [Constraints, P_chp P_pv P_wind P_disch P_buy - P_ch - P_sell P_load]; % 热功率平衡 Constraints [Constraints, H_chp H_boiler H_disch_tank - H_ch_tank H_load]; % CHP可行域与燃料关系 Constraints [Constraints, H_chp_min H_chp H_chp_max]; Constraints [Constraints, P_chp_min P_chp P_chp_max]; Constraints [Constraints, P_chp a1*H_chp b1]; Constraints [Constraints, P_chp a2*H_chp b2]; Constraints [Constraints, F_chp c1*P_chp c2*H_chp c3]; % 锅炉 Constraints [Constraints, 0 P_boiler P_boiler_max]; Constraints [Constraints, H_boiler eta_boiler * P_boiler]; % 电储能递推向量切片形式避免写循环 Constraints [Constraints, SOC(2:T1) SOC(1:T) eta_ch*P_ch*dt/E_cap - P_disch*dt/(eta_disch*E_cap)]; Constraints [Constraints, P_ch M_ch*u_ch, P_disch M_ch*(1-u_ch)]; Constraints [Constraints, SOC_min SOC SOC_max, SOC(1) SOC_init, SOC(T1) SOC_init]; % 蓄热罐递推 Constraints [Constraints, E_tank(2:T1) E_tank(1:T) H_ch_tank*dt - H_disch_tank*dt]; Constraints [Constraints, H_ch_tank M_tank*u_tch, H_disch_tank M_tank*(1-u_tdisch)]; Constraints [Constraints, E_tank_min E_tank E_tank_max, E_tank(1) E_tank_init, E_tank(T1) E_tank_init]; % 购售电互斥 Constraints [Constraints, P_buy M_buy*u_buy, P_sell M_buy*(1-u_buy)]; %% 4. 目标函数 Objective sum(C_buy.*P_buy - C_sell.*P_sell ... C_gas*(F_chp P_boiler) ... C_om_chp*P_chp C_om_boiler*P_boiler); %% 5. 求解 ops sdpsettings(solver,gurobi,verbose,1); optimize(Constraints, Objective, ops); %% 6. 结果输出 P_chp_opt value(P_chp); P_buy_opt value(P_buy); SOC_opt value(SOC); % 其余变量同理绘图和制表这里省略这段代码里的M是大数用于线性化互斥约束。注意M不能取太大否则会造成数值病态经验值是取该支路最大可能出力的1.2倍左右。3.3 几个绕不开的实现细节第一燃料变量的处理。我在3.2里用了F_chp作为独立变量而不是直接写P_chp/eta_e。这是为了准确反映CHP的热电联产效益——燃料同时转化为电和热两个出力共享一份燃料输入线性关系F c1P c2H c3是工程上最常用的简化。我见过不少代码直接写目标函数里燃料成本等于P_chp除以电效率结果整个模型偏向锅炉供热CHP的热效率优势完全没体现出来。第二首末状态约束。调度周期首末SOC或罐内储热量要给定相等值否则优化结果会出现首时段白嫖能量的假象——比如初始SOC设0模型会故意在最后一个小时把储能放空把成本压力转移到周期外。我通常是让SOC(1)和SOC(T1)都等于0.5倍的容量。第三时间索引对齐。SOC和E_tank都有T1个点而P_ch、P_disch、H_ch_tank等只有T个点。递推方程里如果用循环索引从2到T1如果用向量切片SOC(2:T1) 对应 SOC(1:T)长度必须一致。我写过一版把SOC(1:24)和SOC(2:25)直接对齐的代码跑出来全是NaN查了半天才发现是索引错位。4. 三种典型场景的调度结果与解读4.1 场景设计与参数设置为了验证模型我设计了三个递进式的场景场景ACHP燃气锅炉电储能无蓄热罐场景B在A的基础上加入蓄热罐场景C在B的基础上加入高比例光伏和风电分时电价取三段式峰段1.2元/kWh8:00-11:0018:00-21:00平段0.7元/kWh谷段0.35元/kWh23:00-次日7:00。CHP容量2MW锅炉1.5MW电储能1MWh/0.5MW蓄热罐2MWh冬季日热负荷峰值1.2MW电负荷峰值1.8MW。天然气价2.5元/m3机组发电效率0.35热电比取1.2。4.2 结果对比蓄热罐和可再生能源的贡献跑完之后我汇总了三个场景的关键指标指标场景A场景B场景C总运行成本(元)152301289010120购电量(kWh)3100980620CHP发电量(kWh)83001180010500弃风弃光电量(kWh)00850蓄热罐日循环次数-2.11.4场景A的结果非常有代表性因为没有蓄热罐CHP机组被迫完全跟着热负荷走。夜间热负荷高哪怕电价处于谷段CHP也得满发发出来的电用不完只能低价卖给电网白天热负荷降下来CHP跟着降出力可这时候电负荷正是高峰缺口全靠高价购电补。这就是典型的谷段卖电、峰段买电倒挂局面成本自然最高。场景B加蓄热罐之后效果立竿见影。蓄热罐相当于给CHP的热尾巴装了一个缓冲池夜间热负荷高时CHP多发的热能存在罐里电出力不用压得过狠白天需要高负荷出力的时候罐里的热释放出来补热缺口CHP的电出力就被解绑了。CHP发电量提高了42%购电量下降了约68%这就是多能互补在设备层面发挥作用的直接证据。场景C加入高比例风光后问题变成了弃风弃光。850kWh的弃电量说明光伏大发的中午时段连储能都充不满这时候如果缺乏上网激励政策模型选择弃掉是最经济的选择。这也反过来验证了一个道理可再生消纳率不是越高越好它要和经济目标放在同一个天平上权衡。4.3 结果里值得留意的规律多跑几组参数之后我总结出一个规律蓄热罐容量不是越大越好而是要和CHP的热灵活性窗口匹配。所谓热灵活性窗口就是CHP在夜间能输出的最大热出力乘以夜间时长再减去夜间实际热负荷。罐子容量超过这个理论上限它就永远装不满属于投资浪费低于这个值电热耦合约束又会局部复发调度重新被以热定电绑架。工程上一般取理论值的0.8到1.2倍做经济性校核。5. 实测排坑跑模型时最容易翻车的五个地方5.1 求解器路径和license是第一个坑YALMIP装好之后默认是找不到Gurobi的。假如你把Gurobi装好了却在优化时报错说solver未找到多半是YALMIP的路径没配好。更隐蔽的是Windows环境下Gurobi需要设置环境变量GRB_LICENSE_FILE否则明明装好了却报license缺失。我建议装完YALMIP先跑一次yalmiptest看输出里gurobi和cplex的检测状态是不是OK这一步能省掉后面所有为什么求解器没生效的排查时间。5.2 模型不可行时先砍耦合再砍平衡Problem infeasible大概是这个方向出现频率最高的报错。我的排查顺序是固定的先把蓄热罐和电储能的SOC递推约束注释掉用独立变量替代。如果可解问题出在时间耦合约束上。再把电热平衡约束里的可再生能源出力项先置0。如果可解问题出在风光预测与可用容量的矛盾上。最后检查二进制变量对应的M值是否过小。如果M小于该支路的最大可能出力约束会人为制造不可行。这套从强耦合到弱耦合逐步松绑的方法比盯着报错信息瞎猜有效得多。5.3 量纲不统一是结果诡异的最大来源模型里同时出现kW、kWh、元、MW、m3/h这些单位任何一个换算错结果都会合理得离谱。我踩过一次大坑天然气热值约9.7kWh/m3锅炉效率0.9我把燃料的m3/h直接当成了kW乘以效率结果锅炉出力凭空小了约9倍模型拼命让CHP多发电去供热整个调度策略完全变形。最后排查的办法是画出每个设备的能量流从输入端到输出端每一步都换算到kW再进模型。你可以在初始化区加一段统一的单位换算注释别让单位问题散落在各个约束里。5.4 互斥约束的两种写法按场景选储能充放互斥和购售电互斥理论上用0-1变量最严谨我在3.2的代码里也是这么写的。但实际做滚动优化或灵敏度分析时二进制变量会显著拖慢求解速度而且M值选取不当还会引起数值奇异。这种情况我会改用惩罚大M法目标函数里给P_ch和P_disch同时非零的情况加一个很大的惩罚系数效果和互斥约束接近但求解稳定性好很多。简单说单次日前调度用0-1变量滚动优化用惩罚法。5.5 别忘了给变量加非负约束这听起来像废话但sdpvar默认对变量没有符号限制。我第一版代码里锅炉燃料和SOC都出现过负值目标函数甚至算出了负成本——因为你允许了负购电相当于凭空卖电赚钱。所有物理功率变量P_chp、H_chp、P_boiler、P_buy、P_ch、P_disch等都必须手动追加0约束。这一步在代码里看起来多余但漏掉一个就会让整个优化结果失真且极难发现。6. 扩展方向与我的最终一点体会这套代码再往后扩展性价比比较高的方向有三个一是把日前确定性调度升级成日前日内滚动两阶段日内用实测数据更新预测能显著抵消风光预测误差的影响二是把目标函数从纯成本最小扩展到考虑热舒适度的需求响应让可调温控负荷和蓄热罐联动三是多微网互联时电热耦合会跨网络传递单微网优化要升级为分布式协同优化YALMIP里的多目标或bilevel结构可以承担一部分工作。最后说句掏心窝的话做热电联供微网优化最大的坑往往不是代码写不出来而是模型假设与实际设备特性的错位。我见过不少论文里的模型蓄热罐效率取0.95、锅炉效率取0.9、CHP可行域画成完美的矩形——跑出来的结果确实漂亮但拿到工程现场根本没法用。真正能落地的模型是把每一个参数的合理性问一遍这个数哪来的再开始写代码。参数错了求解器再快也是白搭。试过的都知道这活儿急不得慢工出细活。
