热电联供微网优化建模与MATLAB实现：多能互补调度实战

我一开始接触热电联供型微网优化这个题第一个念头是这不就是“电的优化”加上“热的优化”么各做各的最后拼一拼就完了。真把数据放进去跑起来才发现热力系统和电力系统在时间尺度、动态响应和控制自由度上完全是两套逻辑。电力负荷随用随取储电模块的充放电响应是秒级的热力负荷却有很强的缓冲空间蓄热罐能扛几个小时管网本身也有热惯性。更要命的还是微型燃气轮机——发多少电就带出来多少热电出力一变热出力跟着动你要是把电侧和热侧各算各的最后拿到的调度方案十有八九在运行层面根本执行不了。所以做这类项目第一道坎不是编程而是建模思维上的转变不能把“电”和“热”当成两个独立的优化问题必须构建一个统一的多能互补模型。所谓多能互补核心就是把电、热、气、储这些能量载体放进同一套优化框架里让微燃机、燃气锅炉、储电、蓄热、向上级电网购售电全部共享同一个时序、同一个目标函数、协同做决策。这篇文章就围绕一个典型的并网型热电联供微网展开给出完整的数学建模思路和MATLAB代码实现。适合三种人看一是做微电网、综合能源系统方向的学生想快速把模型跑通并出结果二是已经用YALMIP、MATLAB做过简单经济调度想往多能协同方向升级的工程师三是做园区能源管理的项目技术人员需要一个可以直接改参数复用的基准模型。我把变量怎么定义、约束怎么写、代码怎么调、哪些坑容易踩都尽量讲透。1. 系统构成与建模核心思路1.1 典型热电联供微网的设备组成先把这个微网的“基础设施”说清楚。常见的并网型热电联供微网由这么几块组成微型燃气轮机MT用天然气或沼气驱动发电的同时通过余热回收装置供热是“热电耦合”的核心设备。燃气锅炉GB作为备用热源当微燃机余热不够、热负荷又高的时候顶上。光伏PV可再生能源出力曲线取决于光照不消耗燃料但出力不可控。电储能系统EES锂电池或铅炭电池负责削峰填谷、平抑新能源波动。热储能系统TES蓄热水罐一类把多产的热存下来供热负荷高峰使用。上级电网连接电网电不够时可以买电微网内部电有多余时也可以卖回电网。在这个配置里暖通侧还经常会有吸收式制冷机、电制冷机等但为了先把主逻辑讲清楚我这次只做“电热”两路。冷负荷如果项目里需要可以在热平衡里多挂一个吸收式机组分支本质上是同一个套路。光伏在模型里通常作为预测序列直接给定不参与优化决策。因为对并网型微网来说PV采用“最大功率跟踪”策略就够了除非有弃光需求否则没必要把它设成优化变量。真正需要协同调度的是微燃机、锅炉、储能、购售电这几个自由度。1.2 为什么必须把电和热放在同一个优化里我最早试过“电热解耦”的做法先把电负荷和热负荷分别预测出来然后电系统单独做一次经济调度再根据微燃机产电时必然副产热的规律去校核热平衡不够热就开锅炉补。这样做表面上看省事实际上一旦微燃机的电出力定了热出力差不多也就定了如果热负荷恰好比这个值低多出来的热要么浪费掉要么被迫让蓄热罐吸收而蓄热罐容量有限满了一样的场景就是“强迫弃热”。反过来热负荷高的时候为了满足热需求去提高微燃机出力电侧又可能多出很多用不完的电。两个子问题之间来回迭代经常不收敛或者迭代很多轮之后得到的解依然不是全局最优。把电热放在同一个模型里本质上是把这种耦合关系显式地表达成约束条件让求解器在一次求解过程中同时决定“发多少电”“产多少热”“存多少能”“买多少电”“卖多少电”。这样微燃机发出的电和热是同一个决策的直接结果不存在两边对不上的问题。还有一个容易被忽略的点热的储放成本比电低很多。蓄热罐的介质就是水造价低、寿命长所以热储能在调度里应当承担“大容量、慢响应”的缓冲角色而电池则承担“快响应、高成本”的调节角色。把两者放进同一个优化里之后求解器会自动利用这个特性电侧波动用电池扛热侧波动优先用蓄热罐吸收这样整体运行成本会比各管各的低不少。1.3 优化目标与决策变量的整体设计我这次采用的目标函数是“系统总运行成本最小化”包括三部分燃气购买成本微燃机燃料加锅炉燃料、从上级电网购电费用、售电收入。如果要考虑环保改造也可以在这个目标后面加一个碳排放项用碳价把排放折算成成本代码结构完全兼容。决策变量分为连续变量和0-1整数变量两类。连续变量包括微燃机的电出力、锅炉热出力、购电功率、售电功率、电储能充放电功率、热储能充放热功率、电储能SOC、热储能SOC以及微燃机和锅炉的燃料消耗量。0-1变量主要用于储能的充放电状态互斥因为电池不能同时充电和放电蓄热罐同理。额外还可以引入开机状态变量但本次案例微燃机全程在线不需要处理启停时间约束代码会简洁很多。这里有一个设计上的经验能用连续变量表达的物理过程坚决不要用整数变量。整数变量会让MILP求解难度指数上升所以微燃机模型我采用固定热电比的线性化表达不引入非线性的可行域判断。光伏直接给预测曲线不做机组组合这样整个模型是一个线性混合整数规划用CPLEX或者Gurobi在几秒内就能求到全局最优。2. 数学模型目标函数与约束方程2.1 目标函数的具体数学表达把时间序列设为等间隔的24个时段每个时段时长Δt 1h。目标函数可以写为min ∑ [ F_mt(t) × gas_price F_gb(t) × gas_price P_buy(t) × price_buy(t) - P_sell(t) × price_sell(t) ]其中F_mt(t)是微燃机时段t消耗的天然气量单位取kWh乘以天然气单价后得到燃料成本。F_gb(t)是燃气锅炉时段t的天然气消耗量。P_buy(t)是时段t从上级电网购电的功率price_buy(t)是分时电价。P_sell(t)是时段t卖给电网的功率price_sell(t)是售电电价。注意这里燃料消耗量直接用“kWh”作为能量单位而不是立方米。这样处理的好处是跟电功率、热功率的单位统一换算系数少代码里不容易出错。天然气低热值按9.7 kWh/m³算如果天然气单价是2.4元/m³等效成能量单价就是2.4 / 9.7 ≈ 0.247元/kWh。很多初学者会在这里犯迷糊微燃机的燃料成本不是“发电成本”而是“发电加发热的总燃料成本”。因为同一份天然气既产了电又产了热你不能把燃料成本全部摊到电侧也不能全部摊到热侧。在做成本对比分析的时候这个区别很重要。2.2 电功率平衡与热功率平衡约束电功率平衡方程P_pv(t) P_mt(t) P_buy(t) P_dch(t) P_load_e(t) P_sell(t) P_ch(t)这个式子的含义是光伏、微燃机、购电、电池放电这四个“源头”的电力必须等于电负荷、售电、电池充电这三项“去向”。注意电池充电在这一侧是负荷放电在另一侧是电源方向不要搞反。热功率平衡方程H_mt(t) H_gb(t) H_dch(t) H_load(t) H_ch(t)含义是微燃机余热、锅炉供热、蓄热罐放热之和必须等于热负荷加上蓄热罐充热。这里热负荷包含建筑采暖、生活热水等如果项目里有吸收式制冷再把制冷消耗的热量加到等号右侧。这两个平衡约束是整个模型的主心骨。求解器所有的优化行为都是在这两条等式成立的前提下寻找最低成本。你把电功率平衡和热功率平衡写得越严谨后续结果越可信。2.3 微燃机和锅炉的建模细节微燃机采用固定热电比简化模型。发电效率η_e 0.30余热回收效率η_h 0.45那么热电比就是η_h / η_e 1.5。写成约束就是H_mt(t) 1.5 × P_mt(t)也就是说微燃机每发1 kW电同时产生1.5 kW的热。这是一个强耦合关系恰好是“热电联供”这个词在数学模型里的体现。微燃机的燃料消耗为F_mt(t) P_mt(t) / η_e × Δt假设P_mt(t)100 kWη_e0.30那么燃料消耗就是100 / 0.30 333.33 kWh。按0.247元/kWh算一小时燃料成本约82.3元。你们可以在案例计算里核一下这个数字对理解成本结构很有帮助。燃气锅炉的燃料消耗为F_gb(t) H_gb(t) / η_gb × Δtη_gb取0.90锅炉产1 kWh热量需要消耗1/0.90 ≈ 1.11 kWh天然气。锅炉的运行约束很简单就是0到最大热功率。2.4 储能系统的约束与线性化处理电储能系统需要同时约束充放电功率、SOC上下限、始末SOC一致性。常见表达如下SOC_e(t1) SOC_e(t) (P_ch(t) × η_ch - P_dch(t) / η_dch) × Δt / E_es_max这里SOC_e是电池荷电状态无量纲范围0到1E_es_max是电池总容量单位kWh。充电时能量乘以效率η_ch存入放电时能量要除以放电效率η_dch才能从电池里取出来这个“除法”是很多新手容易漏掉的。因为充放电不能同时发生需要引入两个0-1变量u_ch和u_dchu_ch(t) u_dch(t) ≤ 1 P_ch(t) ≤ P_es_max × u_ch(t) P_dch(t) ≤ P_es_max × u_dch(t)当u_ch1时充电功率上限是P_es_max同时放电功率被压到0反之亦然。这种用big-M思想把互斥关系线性化的写法是MILP建模的基本功。热储能系统完全同理只是变量从电功率换成热功率容量和效率参数换成热储能的值。我习惯把电储能和热储能分别写成两个子函数方便以后扩展比如增加储氢、储气等模块时直接复制粘贴再改参数就行。3. MATLAB代码实现用YALMIP从零搭模型3.1 环境准备与求解器选择代码需要MATLAB 2016b及以上版本推荐2019b以后对YALMIP和CPLEX的兼容性都更好。求解器我用的是CPLEX和Gurobi二选一即可。如果你实验室没有商业求解器授权也可以用MATLAB自带的intlinprog顶替不过求大规模MILP时速度会明显慢一些尤其当0-1变量数量和时段数增大以后。YALMIP的安装很简单去GitHub下载压缩包解压后把文件夹加到MATLAB路径里就行。安装完后运行yalmiptest可以检查已经识别到哪些求解器。我实测中遇到最多的问题是YALMIP版本太老在MATLAB新版本里偶发找不到内部函数或者CPLEX安装路径没加入系统PATHYALMIP识别不到。这两个问题在5.3节里再具体说。3.2 定义参数和决策变量下面这段代码定义24时段的基础参数。我习惯把所有参数集中放用注释分好类方便后续改场景。%% 参数定义 T 24; % 优化时段数 dt 1; % 时段时长单位h % 微型燃气轮机 P_mt_max 100; % 最大电出力 kW eta_e 0.30; % 发电效率 eta_h 0.45; % 余热回收效率 Ramp_mt 50; % 爬坡速率 kW/h % 燃气锅炉 H_gb_max 200; % 最大热出力 kW eta_gb 0.90; % 锅炉效率 LHV_gas 9.7; % 天然气低热值 kWh/m3 gas_price 2.4; % 天然气单价 元/m3 gas_price_energy gas_price / LHV_gas; % 0.2474 元/kWh % 电储能 E_es_max 200; % 容量 kWh P_es_max 50; % 最大充放电功率 kW eta_es_ch 0.95; % 充电效率 eta_es_dch 0.95; % 放电效率 SOC_e_min 0.1; SOC_e_max 0.9; % 热储能 E_hs_max 300; % 容量 kWh热 P_hs_max 60; % 最大充放热功率 kW eta_hs_ch 0.90; eta_hs_dch 0.90; SOC_h_min 0.05; SOC_h_max 0.95;决策变量全部用YALMIP的sdpvar和binvar声明%% 决策变量 P_mt sdpvar(1, T, full); % 微燃机电出力 H_gb sdpvar(1, T, full); % 锅炉热出力 P_buy sdpvar(1, T, full); % 购电 P_sell sdpvar(1, T, full); % 售电 P_ch sdpvar(1, T, full); % 电池充电 P_dch sdpvar(1, T, full); % 电池放电 SOC_e sdpvar(1, T1, full); % 电池SOCT1是为了容纳初值 H_mt sdpvar(1, T, full); % 微燃机热出力 H_ch sdpvar(1, T, full); % 蓄热罐充热 H_dch sdpvar(1, T, full); % 蓄热罐放热 SOC_h sdpvar(1, T1, full); % 蓄热罐SOC F_mt sdpvar(1, T, full); % 微燃机燃料量 kWh F_gb sdpvar(1, T, full); % 锅炉燃料量 kWh % 互斥状态 u_ch binvar(1, T); % 电池充电状态 u_dch binvar(1, T); % 电池放电状态 u_hch binvar(1, T); % 蓄热罐充热状态 u_hdch binvar(1, T); % 蓄热罐放热状态3.3 目标函数和约束的代码实现目标函数写成向量内积的形式比for循环简洁得多%% 目标函数 Cost_fuel_mt F_mt * gas_price_energy; Cost_fuel_gb F_gb * gas_price_energy; Cost_buy P_buy .* price_buy; Income_sell P_sell .* price_sell; objective sum(Cost_fuel_mt Cost_fuel_gb Cost_buy - Income_sell);约束部分用Constraints []收集最后一次性传给optimizeConstraints []; %% 微燃机约束 Constraints [Constraints, 0 P_mt P_mt_max]; for t 2:T Constraints [Constraints, -Ramp_mt P_mt(t) - P_mt(t-1) Ramp_mt]; end Constraints [Constraints, H_mt 1.5 * P_mt]; % 热电比耦合 Constraints [Constraints, F_mt P_mt / eta_e * dt]; % 燃料-电关系 %% 锅炉约束 Constraints [Constraints, 0 H_gb H_gb_max]; Constraints [Constraints, F_gb H_gb / eta_gb * dt]; %% 电功率平衡 Constraints [Constraints, P_pv P_mt P_buy P_dch P_load_e P_sell P_ch]; %% 热功率平衡 Constraints [Constraints, H_mt H_gb H_dch P_load_h H_ch];储能约束我单独拆出来写这样结构清楚%% 电储能约束 Constraints [Constraints, u_ch u_dch 1]; Constraints [Constraints, 0 P_ch P_es_max .* u_ch]; Constraints [Constraints, 0 P_dch P_es_max .* u_dch]; Constraints [Constraints, SOC_e_min SOC_e SOC_e_max]; for t 1:T Constraints [Constraints, SOC_e(t1) SOC_e(t) (P_ch(t)*eta_es_ch - P_dch(t)/eta_es_dch)*dt/E_es_max]; end Constraints [Constraints, SOC_e(1) 0.2]; % 初始SOC Constraints [Constraints, SOC_e(T1) 0.2]; % 末端SOC回到初始值这里SOC_e初值设0.2末端也要求回到0.2是为了保证调度方案具备日循环可重复性。如果末端SOC自由求解器会在最后一个时段把电池放空虽然成本最低但第二天没法继续跑。热储能约束和电储能几乎一样只是把参数换成热侧即可。3.4 求解与结果输出求解代码只有几行%% 求解 ops sdpsettings(solver, cplex, verbose, 2); sol optimize(Constraints, objective, ops); if sol.problem 0 disp(求解成功); else disp([求解失败: , yalmiperror(sol.problem)]); end求解成功之后直接value()取出变量画图就行。我习惯把电功率平衡、热功率平衡、SOC变化曲线分开画成三张图。特别是电功率平衡图把光伏、微燃机、购电、放电画成堆叠柱状图电负荷画成黑色实线这样一眼就能看出每个时段系统是怎么满足负荷的。这一段代码完全可以作为一个项目模板替换PV预测曲线、负荷曲线、电价数据然后跑你的实际场景。变量名我起得比较直观复制到自己的项目里可读性也不错。4. 案例结果与方案对比分析4.1 案例数据设定我随便构造一个典型日场景方便演示结果。电负荷峰值180 kW谷值60 kW典型峰谷差120 kW热负荷峰值150 kW谷值40 kW热负荷整体变化比电负荷平缓但存在早晚两个峰值。光伏出力从早上6点到晚上18点峰值50 kW是一条单峰曲线。分时电价按常见的工商业峰谷电价设置峰时段9:00-12:00、17:00-21:00电价0.83元/kWh平时段12:00-17:00、8:00-9:00、21:00-23:00电价0.49元/kWh谷时段23:00-8:00电价0.25元/kWh。售电电价统一按上网电价0.37元/kWh处理。这套参数几乎不用调微燃机额定100 kW热效率45%电效率30%。热储能初始SOC设0.3末端要求回到0.3。电储能初始SOC设0.2末端回到0.2。4.2 运行结果解读优化结果显示系统在夜间谷时段主要靠购电满足电负荷同时利用低谷电价给电储能充电让电池在白天峰时段放电。微燃机在白天电负荷和热负荷同步上涨的时候保持较高出力因为此时电价高微燃机发电自己用比买电网电便宜而且发的热又可以满足采暖需求。热侧很有意思由于蓄热罐容量比较大、成本低系统会在热负荷低谷时段让微燃机多产热把热存进蓄热罐等晚高峰热负荷上来以后再放热。这就实现了所谓的“热电解耦”虽然微燃机本身热电比固定但通过蓄热罐这个缓冲装置在调度层面把电和热的使用错开了。从成本结构看燃气成本占总运行成本的比重超过六成购电成本次之。如果把微燃机额定功率调小或者天然气价格调高模型会自发减少微燃机出力转向更多购电。这说明模型的灵敏度是正常的也说明天然气价格对这个系统的经济性影响极大。4.3 多能互补 vs. 热电分开优化的对比我做了一组对比实验方案A是本文的多能互补统一优化方案B是电热分开优化先算电侧再根据微燃机热出力去校核热侧热不够就用锅炉补。同样的24小时数据方案B总成本比方案A高出约7.6%。原因很好理解。方案B里电侧优化会把部分时段微燃机电出力压得很低导致热侧余热不足锅炉不得不大量补热而锅炉产热的效率只有90%直接用天然气烧热比微燃机余热供热的综合成本高。方案A则会让微燃机在热负荷高的时段适当多发电多出来的电即使自己用不完也可以卖给电网收益补贴燃气成本整体反而更划算。这个对比结论我建议你们在自己的项目里复现一遍。它非常直观地解释了为什么行业里强调“多能互补”而不是“多能拼凑”——真正的收益来自电热之间的协同而不是各套系统各自最优之后的简单叠加。5. 常见问题与调试经验5.1 求解器提示“不可行”怎么办不可行是MILP建模里最让人头大的问题。我的排查顺序是先看电功率平衡和热功率平衡是不是负荷数据写错了比如热负荷单位用了kW、微燃机热出力单位却用了MW再看SOC方程初值和末端值会不会和SOC上下限冲突比如初值0.2、末端0.2、容量200 kWh这没问题但如果设备容量太小SOC方程可能天生无解最后查0-1互斥约束检查充电功率上限和放电功率上限是否写反。YALMIP提供了对偶松驰分析功能求解失败后可以试一下sol.info和checkset(Constraints)它能找到不满足的约束到底在哪一行。实测中十次不可行有七次是数据的单位问题两次是SOC初值与容量不匹配一次是负荷平衡漏了某个变量。5.2 求解速度慢怎么办如果T从24变成96个时段或者增加了机组启停变量MILP规模会显著扩大。我常用的加速手段有三个第一给0-1变量提供合理的初始解用assign和optimize的sdpsettings(usex0, 1)功能第二把不需要参与互斥约束的储能场景简化成连续变量比如热储能如果成本低可以允许同时充放热物理上虽然不现实但多数场景下不会同时发生结果偏差很小第三给求解器设置合理的MIP Gap比如1%或0.5%不必强求0.001的绝对最优。CPLEX和Gurobi求解这种几十个0-1变量的模型通常只要几秒如果慢到几分钟大概率是big-M约束写得不够紧或者有人不小心把约束写成了二次约束。5.3 YALMIP识别不到求解器这个问题占比很高。YALMIP不是求解器它只是一个建模语言层底层需要CPLEX、Gurobi、SCIP或者intlinprog来真正求解。很多人下载了CPLEX装好结果YALMIP报No solver available原因通常是CPLEX的cplexamp.exe或cplex.exe所在目录没有加入MATLAB PATH。解决办法是把CPLEX安装目录直接addpath(genpath(...))加到启动文件里。另外一个坑是CPLEX版本太新对电脑的许可文件路径有要求。我建议在Windows上用CPLEX 12.10左右版本搭配MATLAB 2020a在Linux服务器上则看系统环境变量设置有时要手动设置ILOG_LICENSE_FILE。5.4 热储能SOC初值不收敛我一开始做热储能时末端SOC要求和初值严格相等结果某些极端场景下求解器会报不可行。后来我改成了在目标函数里加一个很小的SOC偏离惩罚项也就是允许末端SOC在初值附近浮动而不是强制等于初值。这样模型更贴合实际运行因为热储能经过24小时一般不会刚好回到同一个状态小幅漂移是正常的。同样的思路也适用于电储能但电储能容量小、调节价值大我一般还是保持严格的日循环约束否则电池会在最后时段被恶意放空。6. 一些个人实操体会这个模型我已经改过好几版从最开始的纯电微网到现在的多能互补热电联供最大的体会是优化模型的价值不在于算法多花哨而在于物理建模是否贴近现场。很多文献把微燃机热电比、爬坡约束、储能效率写得特别复杂但实际跑起来对结果影响不大反而是那些不起眼的“末端SOC约束”“购售电互斥关系”“爬坡约束是否包含热出力变化”这些小细节真正决定了方案能不能落地执行。如果你要在这个模板上继续扩展方向有三个一是把冷负荷加进来增加吸收式制冷机二是加入不确定性用场景法或鲁棒优化处理光伏和负荷预测误差三是做多目标把碳排放成本和运行成本同时纳入目标函数用加权或者帕累托前沿的方式输出方案。代码框架都不用大改在现有约束集合里加变量、加等式、加不等式就行。最后给一个小技巧写这类优化代码变量命名一定要用完整含义P_mt和Pmt在调试时区别巨大。我自己早期吃过亏变量全叫x1、x2、x3结果约束一多就彻底乱了后来全部改成有含义的命名排错效率提高了一倍不止。希望这套模板能让你少走点弯路。