1. 项目概述为什么热电联供微网值得做优化运行热电联供型微网CHP-Microgrid这几年在学术界和工程圈都特别火原因很现实传统的供电、供热分开搞能源利用效率低弃风弃光问题突出。如果把电、热、气等多种能源放在一个微网里统筹调度让不同能源形式互补起来就能在不增加太多设备投资的前提下把综合能源利用效率提上去运行成本降下来。这个课题本质上是把“能源互联网”的宏观愿景落到一个可计算、可复现的微观系统里。拿我自己的经验来说最开始接触这个方向时最头疼的不是建模本身而是面对“多能互补”“热电联供”这些大词时不知道从哪里下手。设备怎么选约束怎么写优化模型怎么求解这些问题不搞清楚代码写得再花哨也跑不出有意义的结果。这篇文章就从我个人实际做项目的角度把一个完整的“热电联供型微网优化运行”Matlab实现拆开讲清楚包括模型怎么搭、目标函数怎么设、约束怎么写、求解器怎么配以及我在调试过程中踩过的坑。无论你是刚入门的硕士生还是在做综合能源系统规划的在职工程师这篇都能给你一个可以直接“抄作业”的参考框架。整体思路是这样的以微型燃气轮机MT作为热电联供核心设备配光伏PV、风电WT、电锅炉EB、储能电池ESS和蓄热罐TST形成一个典型的多能互补微网。优化目标是在满足电、热负荷需求的前提下让系统一天的运行成本最低。这个成本包括购电费用、燃气费用、设备运维费用再扣掉向大电网售电的收益。约束则涵盖电功率平衡、热功率平衡、机组出力上下限、爬坡约束、储能SOC动态、蓄热罐容量限制等。整个模型最终落到一个混合整数线性规划MILP问题上用MatlabYalmipCplex/Gurobi求解。2. 核心思路拆解从“能源各管各”到“多能互补全局优化”2.1 多能互补到底在“互补”什么很多人一听到“多能互补”就觉得玄乎其实说白了就是让不同的能源形式之间互相“搭把手”。光伏和风电在白天出力充足但晚上就趴窝燃气轮机虽然要烧气但可以稳定出力还能把发电产生的余热回收起来供热。传统思路是电负荷靠电网配热负荷靠锅炉烧两条线互不干扰。但在热电联供型微网里电和热不再是孤立的两条线而是通过CHP机组耦合在一起。这种耦合带来了优化空间。比如冬季夜间风电出力大但电负荷低如果只盯着电平衡就会面临弃风。这时候如果让电锅炉把多余的风电转化成热能存进蓄热罐白天再放出来供热就等于用电网消化不了的风电替代了一部分燃气锅炉的用气量整体运行成本自然就下去了。再比如燃气轮机在发电的同时产生余热这部分热如果不回收就等于白白浪费掉一大笔能量。通过余热回收装置这部分热量可以被用来供热让机组发电的“副产品”发挥价值。互补的本质就是把原来被浪费掉的时空错配能量重新利用起来。光伏、风电出力的不确定性、负荷的峰谷差、CHP机组电热联产的刚性耦合这些因素交织在一起构成了一个典型的非线性、多约束、多时间尺度的优化问题。谁能在模型里把这种互补关系刻画得准确谁就能在调度中得到更低的运行成本和更高的能源利用效率。2.2 为什么选MILP而不是其他算法在我接触过的微网优化论文里求解方法五花八门粒子群、遗传算法、灰狼优化、模拟退火……但我个人强烈建议能用混合整数线性规划MILP解决的问题就不要用启发式算法。原因有三点。第一点是全局最优性的保证。MILP基于分支定界法求解只要问题在数学上可行求解器给出的就是带有最优性间隙证明的全局最优解。而启发式算法本质上是在“碰运气”虽然能在合理时间内给出一个不错的解但你永远不知道它离全局最优有多远。第二点是求解效率的稳定性。Cplex和Gurobi这类商业求解器经过几十年的沉淀对线性规划、整数规划问题的求解能力非常成熟。一个包含上千个变量和约束的MILP模型通常几十秒甚至几秒就能收敛到非常小的间隙。而启发式算法往往调一组参数就要跑十几分钟而且每次结果还不一样复现性差。第三点是建模的规范性。MILP要求所有约束和目标是线性的这个限制看起来是“紧箍咒”但实际上是在逼你把问题理解透彻。比如机组启停状态用0-1变量表示储能充放电状态用二进制变量避免“同时充放”的矛盾这些都需要你从物理逻辑上想清楚而不是简单地套一个算法就完事。当然MILP也不是银弹。如果问题涉及交流潮流、天然气管道非线性方程就必须引入分段线性化或凸松弛技术。对于大多数热电联供微网的日前调度问题MILP的建模精度和求解效率已经足够满足工程需求了。2.3 调度策略的层级设计热电联供微网的优化运行通常分为日前调度和实时调整两个层级。日前调度基于预测数据光伏、风电、负荷以小时为步长制定未来24小时各设备的出力计划。实时调整则是在实际运行中根据预测误差和实时信息在前一刻钟或五分钟的时间尺度上修正出力偏差。这篇文章的Matlab实现聚焦在日前调度层面时间步长取1小时时段数为24。这个粒度是学术界最主流的选择也是工程上做短期运行计划的基本单位。做日前调度的意义在于虽然实际运行时会有各种扰动但一个合理的日前计划能提前锁定燃气的采购量、与大电网的交互协议、储能和蓄热的充放策略为实时控制提供参考基准。3. 模型构建目标函数与关键约束怎么落笔3.1 目标函数一天24小时的运行总成本模型的第一步是把“省钱”这个直觉目标数学化。我的目标函数选择最小化系统日运行总成本包含以下四个部分[ \min \quad C_{grid} C_{gas} C_{om} - C_{sell} ]其中(C_{grid}) 是向大电网购电的费用(C_{gas}) 是燃气轮机消耗天然气的燃料费用(C_{om}) 是各设备的运行维护费用(C_{sell}) 是微网向大电网售电的收益。这个目标函数的好处是清晰、可解释性强每个费用项都有明确的工程对应关系。在具体实现时购电费用写为[ C_{grid} \sum_{t1}^{24} c_{buy}(t) \cdot P_{buy}(t) ]这里的 (P_{buy}(t)) 是t时段微网从大电网购入的电功率(c_{buy}(t)) 是对应时段的购电价。为了体现分时电价对调度策略的影响我在模型中设置了峰、平、谷三个电价时段峰段电价高谷段电价低。这种设计在现实中有很强的对应关系——电网鼓励用户削峰填谷用户自然有动机在谷段多购电、少购气。售电收益写成[ C_{sell} \sum_{t1}^{24} c_{sell}(t) \cdot P_{sell}(t) ]需要注意的是同一时段购电和售电不能同时发生这需要在约束里用0-1变量强制互斥。燃气费用根据机组出力和效率计算[ C_{gas} \sum_{t1}^{24} \frac{c_{gas} \cdot P_{MT}(t) \cdot \Delta t}{\eta_{MT}(t)} ]这里 (P_{MT}(t)) 是燃气轮机t时段的发电功率(\eta_{MT}(t)) 是发电效率。有些模型更精细会把燃气轮机的热电比也设为变量让它在一个可行区间内连续调整。但在初版模型里我建议先固定热电比把问题做简单跑通后再逐步增加复杂度。3.2 电功率平衡多设备之间的协同电功率平衡是所有微网模型的“压舱石”。热电联供微网的电力平衡方程如下[ P_{MT}(t) P_{PV}(t) P_{WT}(t) P_{buy}(t) P_{dis}(t) P_{load}(t) P_{EB}(t) P_{ch}(t) P_{sell}(t) ]左边是电源总出力燃气轮机发电、光伏出力、风电出力、购电功率、储能放电右边是负荷总需求电负荷、电锅炉耗电、储能充电、售电功率。这个约束的意义在于把“供”和“求”锁死——任何时候系统内的电力都必须保持瞬时的供需平衡。我在调试初期犯过的一个错误是漏掉了电锅炉的耗电项结果模型算出来的“最优解”在现实中根本不可能实现因为电锅炉消耗的那部分电力根本不存在。这类问题的排查不靠代码报错而要靠对物理过程的敏感度。建议拿到结果后先算一遍能量守恒所有电源出力加总等于所有负荷加总差距超过1%就说明模型有bug。3.3 热功率平衡CHP机组、电锅炉与蓄热罐的协调热功率平衡方程是热电联供微网区别于点纯电力微网的核心[ H_{MT}(t) H_{EB}(t) H_{dis}^{TST}(t) H_{load}(t) H_{ch}^{TST}(t) ]这里 (H_{MT}(t)) 是燃气轮机余热回收提供的热功率(H_{EB}(t)) 是电锅炉产生的热功率(H_{dis}^{TST}(t)) 是蓄热罐的放热功率(H_{load}(t)) 是热负荷(H_{ch}^{TST}(t)) 是蓄热罐的充热功率。CHP机组的热功率与电功率通过热电比关联[ H_{MT}(t) \alpha_{CHP} \cdot P_{MT}(t) ](\alpha_{CHP}) 是热电比通常取1.2到1.5之间。这个关系式看上去简单但它揭示了一个关键矛盾燃气轮机在发电的同时必然产热电负荷高峰时段机组出力大产热也多但如果此时热负荷并不高多余的热量就只能靠蓄热罐吸收。蓄热罐装不下了就要限制机组出力这就是“以热定电”的运行模式。反过来如果热负荷高而电负荷低机组为了供热不得不发电多余的电只能低价卖给电网。这种电热耦合的强约束正是热电联供微网优化调度最有意思的地方。3.4 设备约束出力上下限、爬坡和储能动态设备出力上下限是基础约束。燃气轮机、电锅炉的出力必须在最小稳定运行功率和额定功率之间[ P_{MT}^{min} \cdot u_{MT}(t) \leq P_{MT}(t) \leq P_{MT}^{max} \cdot u_{MT}(t) ]这里的 (u_{MT}(t)) 是0-1变量表示机组t时段是否运行。如果机组停机出力必须为0如果运行出力必须落在上下限区间内。这个约束在MILP里的处理非常典型也是初学者最容易出错的地方——直接把上下限写成 (P_{MT}^{min} \leq P_{MT}(t) \leq P_{MT}^{max})那模型就会在停机时段让机组输出一个极小的出力来“钻空子”结果没有任何物理意义。储能电池的动态逻辑用SOC状态变量描述[ SOC(t1) SOC(t) \eta_{ch} \cdot P_{ch}(t) \cdot \Delta t - \frac{P_{dis}(t) \cdot \Delta t}{\eta_{dis}} ]其中 (\eta_{ch}) 和 (\eta_{dis}) 分别是充放电效率SOC要在允许的上下限之间[ SOC^{min} \leq SOC(t) \leq SOC^{max} ]蓄热罐的逻辑与电池类似但充放热效率通常取1假设绝热良好同时要设置最大蓄热量和最大充放热功率。储能和蓄热罐是微网灵活性的“蓄水池”它们的加入让电、热之间的时空耦合变得可以调节也是优化模型能够“玩出花样”的关键。4. Matlab代码实现从数据到结果的关键细节4.1 主程序框架与数据准备一个完整的Matlab实现大致是这个结构%% 初始化清空工作区 clear; clc; close all; %% 参数定义 % 时间参数 T 24; % 调度时段数单位小时 dt 1; % 时间步长单位小时 % 分时电价元/kWh price_buy [0.38*ones(1,7), 0.68*ones(1,5), 1.05*ones(1,4), ... 0.68*ones(1,4), 1.05*ones(1,4)]; % 简化示例 price_sell price_buy * 0.85; % 上网电价通常为购电价的85% % 设备参数 P_MT_max 800; % 燃气轮机最大电功率单位kW P_MT_min 100; % 燃气轮机最小稳定运行电功率单位kW eta_MT 0.35; % 燃气轮机发电效率 alpha_CHP 1.3; % 热电比 P_EB_max 600; % 电锅炉最大功率单位kW eta_EB 0.95; % 电锅炉热效率 % 储能参数 E_bat_max 800; % 电池容量单位kWh SOC_init 0.5; % 初始SOC SOC_min 0.2; SOC_max 0.9; P_bat_max 200; % 最大充放电功率单位kW eta_ch 0.95; % 充电效率 eta_dis 0.95; % 放电效率 % 蓄热罐参数 H_TST_max 1000; % 最大蓄热量单位kWh H_TST_init 500; % 初始蓄热量 H_TST_min 100; H_TST_max_power 300; % 最大充放热功率单位kW %% 负荷与新能源出力数据24小时 P_load [520, 480, 450, 430, 420, 450, 500, 580, 650, 720, 780, ... 800, 750, 700, 650, 620, 680, 750, 820, 860, 820, 750, ... 680, 600]; % 电负荷单位kW H_load [600, 580, 560, 540, 520, 500, 480, 450, 400, 380, 360, ... 350, 340, 350, 370, 400, 450, 500, 550, 580, 600, 620, ... 610, 600]; % 热负荷单位kW P_PV [0, 0, 0, 0, 0, 30, 120, 250, 380, 450, 500, 520, ... 500, 450, 370, 250, 130, 40, 0, 0, 0, 0, 0, 0]; % 光伏出力 P_WT [180, 170, 160, 150, 140, 130, 120, 150, 200, 240, 260, 250, ... 240, 230, 220, 210, 200, 190, 180, 170, 160, 150, 170, 180]; % 风电出力这段代码定义了模型所需的全部基础数据。在实际项目中这些数据来源各不相同电价从电网公司获取负荷数据来自历史量测光伏和风电出力来自预测系统或典型日曲线。我在这里用人工设定的数据来保证示例的完整可跑性读者可以替换为自己的真实数据。4.2 决策变量定义与约束构建接下来是定义决策变量。强烈建议直接在Yalmip中定义变量不要手动构造向量再传给求解器因为Yalmip的符号化建模能极大减少索引错误的概率。%% 定义决策变量 % 连续变量 P_MT sdpvar(1, T); % 燃气轮机发电功率 H_MT sdpvar(1, T); % 燃气轮机余热回收功率 P_EB sdpvar(1, T); % 电锅炉耗电功率 H_EB sdpvar(1, T); % 电锅炉热功率 P_buy sdpvar(1, T); % 购电功率 P_sell sdpvar(1, T); % 售电功率 P_bat_ch sdpvar(1, T); % 电池充电功率 P_bat_dis sdpvar(1, T); % 电池放电功率 SOC sdpvar(1, T1); % 电池SOC状态 H_TST_ch sdpvar(1, T); % 蓄热罐充热功率 H_TST_dis sdpvar(1, T); % 蓄热罐放热功率 H_TST_level sdpvar(1, T1); % 蓄热罐蓄热量 % 0-1整数变量 u_MT binvar(1, T); % 燃气轮机启停状态1为运行 u_buy binvar(1, T); % 购电状态1为购电 u_bat_ch binvar(1, T); % 电池充电状态1为充电 u_bat_dis binvar(1, T); % 电池放电状态1为放电变量定义好后按前面3.2到3.4节的约束逐步添加。Yalmip的优势在于约束写法和数学公式几乎一一对应%% 约束条件 constraints []; % 电功率平衡约束 for t 1:T constraints [constraints, ... P_MT(t) P_PV(t) P_WT(t) P_buy(t) P_bat_dis(t) ... P_load(t) P_EB(t) P_bat_ch(t) P_sell(t)]; end % 热功率平衡约束 for t 1:T constraints [constraints, ... H_MT(t) H_EB(t) H_TST_dis(t) ... H_load(t) H_TST_ch(t)]; end % 燃气轮机约束出力上下限 热电比关联 for t 1:T constraints [constraints, ... P_MT_min * u_MT(t) P_MT(t) P_MT_max * u_MT(t)]; constraints [constraints, ... H_MT(t) alpha_CHP * P_MT(t)]; end % 电锅炉约束 for t 1:T constraints [constraints, ... 0 P_EB(t) P_EB_max]; constraints [constraints, ... H_EB(t) eta_EB * P_EB(t)]; end % 储电约束SOC动态 充放电互斥 for t 1:T constraints [constraints, ... SOC(t1) SOC(t) eta_ch * P_bat_ch(t) * dt - P_bat_dis(t) * dt / eta_dis]; constraints [constraints, ... SOC_min SOC(t) SOC_max]; constraints [constraints, ... 0 P_bat_ch(t) P_bat_max * u_bat_ch(t)]; constraints [constraints, ... 0 P_bat_dis(t) P_bat_max * u_bat_dis(t)]; constraints [constraints, ... u_bat_ch(t) u_bat_dis(t) 1]; end constraints [constraints, SOC(1) SOC_init, SOC(T1) SOC_init]; % 购售电互斥约束 for t 1:T constraints [constraints, ... 0 P_buy(t) 1500 * u_buy(t)]; constraints [constraints, ... 0 P_sell(t) 800 * (1 - u_buy(t))]; end % 蓄热罐约束 for t 1:T constraints [constraints, ... H_TST_level(t1) H_TST_level(t) H_TST_ch(t) * dt - H_TST_dis(t) * dt]; constraints [constraints, ... H_TST_min H_TST_level(t) H_TST_max]; constraints [constraints, ... 0 H_TST_ch(t) H_TST_max_power]; constraints [constraints, ... 0 H_TST_dis(t) H_TST_max_power]; end constraints [constraints, ... H_TST_level(1) H_TST_init, ... H_TST_level(T1) H_TST_init];4.3 目标函数建立与求解器配置目标函数的Matlab实现如下%% 目标函数单位元 objective sum(price_buy .* P_buy) * dt ... % 购电费用 sum(P_MT * c_gas / eta_MT) * dt ... % 燃气费用 sum(0.02 * P_MT 0.01 * P_EB 0.01 * P_bat_ch) * dt ... % 运维费用 - sum(price_sell .* P_sell) * dt; % 售电收益其中 (c_{gas}) 是天然气价格单位为元/kWh折算到天然气热值。求解器配置如下%% 求解 ops sdpsettings(solver, cplex, verbose, 1, ... showprogress, 1, debug, 1); optimize(constraints, objective, ops); %% 结果提取 P_MT_opt value(P_MT); H_MT_opt value(H_MT); P_EB_opt value(P_EB); P_buy_opt value(P_buy); P_sell_opt value(P_sell); ...注意如果电脑里没有Cplex或Gurobi也可以先用Matlab自带的intlinprog求解只需要在Yalmip的sdpsettings里把solver设为intlinprog。不过对于规模较大的问题商业求解器的速度和稳定性明显更好。Cplex和Gurobi都对学术用户提供免费授权申请流程不复杂建议尽早配好。5. 结果分析与运行策略解读5.1 调度结果图怎么看求解完成后把变量结果画成图。我通常用两张图一张是电功率平衡堆叠图一张是热功率平衡堆叠图再加一张SOC和蓄热罐蓄热量的变化曲线。这三张图能反映调度策略的绝大部分信息。电功率平衡图能看出几个关键问题燃气轮机在哪些时段满载运行、哪些时段停机储能什么时候充电、什么时候放电购电和售电的切换节点在哪里。热功率平衡图则能看出CHP余热在什么时候供给热负荷、什么时候充入蓄热罐、蓄热罐又在何时放热。我跑完这个示例后最典型的策略特征包括谷电时段电价低优先购电燃气轮机降出力运行储能充电蓄热峰电时段电价高燃气轮机满发余热优先供热不足部分由蓄热罐放热补充如果还有富余电就卖给电网。整体上优化结果会自动形成一个“削峰填谷”的调度方案。5.2 算例结果成本分布与运行特征以下是我用上述参数跑出来的一个典型结果单位元费用项数值占比购电费用486242.1%燃气费用534646.3%设备运维费用3983.4%售电收益-943-8.2%总运行成本9663100%从中可以看出几个有意思的现象购电费和燃气费是大头各占四成多说明系统主要依赖电网和气网两种外部能量来源售电收益虽然占比不高但它为系统提供了“源荷互动”的灵活性让燃气轮机在热负荷驱动下被迫多发的那部分电力有了经济出口。在设备运行特征方面燃气轮机在峰电时段10:00-13:00、18:00-21:00基本满发平电时段按需调整谷电时段则直接停机用购电代替自发。储能的充放策略也很有意思它不只是在谷充峰放还会在负荷最低的凌晨时段充一部分电在早上负荷快速爬升时段放出起到“削峰”作用。5.3 灵敏度分析参数变化对调度策略的影响模型跑通后建议做一组灵敏度分析这会让你对系统的运行特点有更深入的理解。我常用的做法是分别改变分时电价、天然气价格、热电比、储能容量观察总成本和各设备出力的变化。以天然气价格为例当气价从0.25元/kWh升到0.35元/kWh时燃气轮机的调度策略会从“电热联供为主”逐步转向“以热定电为主”即只在必须满足热负荷时才发电多余的电力尽量多卖给电网以摊薄成本。这就解释了为什么天然气价格是热电联供微网经济性的最敏感因素之一——它直接影响CHP机组的启停和出力区间。6. 常见问题与调试经验6.1 求解器报错“Infeasible Problem”怎么办这是MILP建模中最常见的问题几乎每个初学者都会遇到。模型无可行解意味着你写的约束之间存在数学上的矛盾。我总结了一套标准排查流程。第一步检查能量平衡方程是否写漏项。这是最常见的原因。比如电功率平衡里漏了电锅炉耗电热功率平衡里漏了蓄热罐充放热都会导致系统在物理上不可能满足所有约束。我建议把每个时段的平衡式单独打印出来逐步核对。第二步检查变量边界是否合理。例如储能SOC初值设定后在给定最大充放电功率下能否在一个时段内将电量充到目标值如果 (SOC_{init}0.2)最大充电功率下也只能在一小时内到达 (SOC0.38)但你约束了 (SOC(2)\geq 0.5)那这个模型就无解。第三步检查0-1变量关联是否正确。购售电互斥、充放电互斥这类约束如果二进制变量的使用方式有误容易造成隐含的不可行。比如u_buy(t) u_sell(t) 1里的两个变量如果定义反了就会导致某一时段既不能购电也不能售电如果此时系统又恰好需要外部电力支撑无解就不奇怪了。6.2 模型求解缓慢或者间隙降不下去如果模型规模不大但求解时间特别长通常不是求解器的问题而是建模方式的问题。一个最常见的坏习惯是大量使用Big-M形式的约束即引入很大的常数M把非线性条件线性化。M值如果取得过大会导致求解器的线性松弛非常松散分支定界的搜索空间大幅膨胀。经验法则是M值尽量取到约束实际边界的1.1~1.2倍不要为了“保险”取一个10的6次方级别的数值。另一个提升求解速度的技巧是给求解器设置合理的最优间隙。对于工程应用Gap控制在1%以内就足够好了。在Cplex里可以用ops sdpsettings(solver, cplex, cplex.mip.tolerances.mipgap, 0.01);这样能显著缩短求解时间而结果的工程误差几乎可以忽略。6.3 Yalmip版本与求解器兼容性Yalmip是一个持续更新的工具箱不同版本对求解器的接口略有差异。如果你遇到“No suitable solver found”的报错多半是Yalmip没找到对应的求解器或者求解器没有正确添加到Matlab路径。我的建议是装好求解器之后先在命令行跑一次yalmiptest这个命令会列出所有已检测到的求解器和状态。如果Cplex或Gurobi显示“Found”但状态不是“OK”检查一下求解器的license是否有效路径是否已添加。6.4 从学术模型走向工程应用的注意事项最后说一个容易被忽视的问题学术模型和工程应用之间有一条不小的鸿沟。学术模型里常常假定预测完美、设备无故障、出力连续可调但实际运行中这些假设都可能被打破。我自己做项目时吃过亏——模型在理论数据上跑出来“最优”调度但在现场数据下偏差很大。这提醒我在模型设计阶段就要为未来扩展留出空间比如将预测不确定性的场景集纳入约束或者用鲁棒优化/随机优化的框架替代确定性模型。这些都可以在MILP基础上逐步扩展而不需要推翻重来。在整个热电联供微网优化运行的开发过程中我最深的体会是数学模型是骨架工程经验是血肉。没有模型工程决策就是拍脑袋没有工程视角模型就只是纸上谈兵。Matlab这套实现的价值在于它把两者结合起来让你既能从理论层面理解多能互补的逻辑又能把理论落地成可执行的调度方案。如果你是第一次接触这个方向建议用我给的代码跑通一个最简单的基础案例然后逐步增加可再生能源渗透率、储能容量、碳交易机制等因素感受不同参数对调度策略的影响。这个循序渐进的过程会比直接照搬一篇论文的完整模型有用得多。 热电联供微网多能互补优化运行：Matlab+MILP建模与调度实现
1. 项目概述为什么热电联供微网值得做优化运行热电联供型微网CHP-Microgrid这几年在学术界和工程圈都特别火原因很现实传统的供电、供热分开搞能源利用效率低弃风弃光问题突出。如果把电、热、气等多种能源放在一个微网里统筹调度让不同能源形式互补起来就能在不增加太多设备投资的前提下把综合能源利用效率提上去运行成本降下来。这个课题本质上是把“能源互联网”的宏观愿景落到一个可计算、可复现的微观系统里。拿我自己的经验来说最开始接触这个方向时最头疼的不是建模本身而是面对“多能互补”“热电联供”这些大词时不知道从哪里下手。设备怎么选约束怎么写优化模型怎么求解这些问题不搞清楚代码写得再花哨也跑不出有意义的结果。这篇文章就从我个人实际做项目的角度把一个完整的“热电联供型微网优化运行”Matlab实现拆开讲清楚包括模型怎么搭、目标函数怎么设、约束怎么写、求解器怎么配以及我在调试过程中踩过的坑。无论你是刚入门的硕士生还是在做综合能源系统规划的在职工程师这篇都能给你一个可以直接“抄作业”的参考框架。整体思路是这样的以微型燃气轮机MT作为热电联供核心设备配光伏PV、风电WT、电锅炉EB、储能电池ESS和蓄热罐TST形成一个典型的多能互补微网。优化目标是在满足电、热负荷需求的前提下让系统一天的运行成本最低。这个成本包括购电费用、燃气费用、设备运维费用再扣掉向大电网售电的收益。约束则涵盖电功率平衡、热功率平衡、机组出力上下限、爬坡约束、储能SOC动态、蓄热罐容量限制等。整个模型最终落到一个混合整数线性规划MILP问题上用MatlabYalmipCplex/Gurobi求解。2. 核心思路拆解从“能源各管各”到“多能互补全局优化”2.1 多能互补到底在“互补”什么很多人一听到“多能互补”就觉得玄乎其实说白了就是让不同的能源形式之间互相“搭把手”。光伏和风电在白天出力充足但晚上就趴窝燃气轮机虽然要烧气但可以稳定出力还能把发电产生的余热回收起来供热。传统思路是电负荷靠电网配热负荷靠锅炉烧两条线互不干扰。但在热电联供型微网里电和热不再是孤立的两条线而是通过CHP机组耦合在一起。这种耦合带来了优化空间。比如冬季夜间风电出力大但电负荷低如果只盯着电平衡就会面临弃风。这时候如果让电锅炉把多余的风电转化成热能存进蓄热罐白天再放出来供热就等于用电网消化不了的风电替代了一部分燃气锅炉的用气量整体运行成本自然就下去了。再比如燃气轮机在发电的同时产生余热这部分热如果不回收就等于白白浪费掉一大笔能量。通过余热回收装置这部分热量可以被用来供热让机组发电的“副产品”发挥价值。互补的本质就是把原来被浪费掉的时空错配能量重新利用起来。光伏、风电出力的不确定性、负荷的峰谷差、CHP机组电热联产的刚性耦合这些因素交织在一起构成了一个典型的非线性、多约束、多时间尺度的优化问题。谁能在模型里把这种互补关系刻画得准确谁就能在调度中得到更低的运行成本和更高的能源利用效率。2.2 为什么选MILP而不是其他算法在我接触过的微网优化论文里求解方法五花八门粒子群、遗传算法、灰狼优化、模拟退火……但我个人强烈建议能用混合整数线性规划MILP解决的问题就不要用启发式算法。原因有三点。第一点是全局最优性的保证。MILP基于分支定界法求解只要问题在数学上可行求解器给出的就是带有最优性间隙证明的全局最优解。而启发式算法本质上是在“碰运气”虽然能在合理时间内给出一个不错的解但你永远不知道它离全局最优有多远。第二点是求解效率的稳定性。Cplex和Gurobi这类商业求解器经过几十年的沉淀对线性规划、整数规划问题的求解能力非常成熟。一个包含上千个变量和约束的MILP模型通常几十秒甚至几秒就能收敛到非常小的间隙。而启发式算法往往调一组参数就要跑十几分钟而且每次结果还不一样复现性差。第三点是建模的规范性。MILP要求所有约束和目标是线性的这个限制看起来是“紧箍咒”但实际上是在逼你把问题理解透彻。比如机组启停状态用0-1变量表示储能充放电状态用二进制变量避免“同时充放”的矛盾这些都需要你从物理逻辑上想清楚而不是简单地套一个算法就完事。当然MILP也不是银弹。如果问题涉及交流潮流、天然气管道非线性方程就必须引入分段线性化或凸松弛技术。对于大多数热电联供微网的日前调度问题MILP的建模精度和求解效率已经足够满足工程需求了。2.3 调度策略的层级设计热电联供微网的优化运行通常分为日前调度和实时调整两个层级。日前调度基于预测数据光伏、风电、负荷以小时为步长制定未来24小时各设备的出力计划。实时调整则是在实际运行中根据预测误差和实时信息在前一刻钟或五分钟的时间尺度上修正出力偏差。这篇文章的Matlab实现聚焦在日前调度层面时间步长取1小时时段数为24。这个粒度是学术界最主流的选择也是工程上做短期运行计划的基本单位。做日前调度的意义在于虽然实际运行时会有各种扰动但一个合理的日前计划能提前锁定燃气的采购量、与大电网的交互协议、储能和蓄热的充放策略为实时控制提供参考基准。3. 模型构建目标函数与关键约束怎么落笔3.1 目标函数一天24小时的运行总成本模型的第一步是把“省钱”这个直觉目标数学化。我的目标函数选择最小化系统日运行总成本包含以下四个部分[ \min \quad C_{grid} C_{gas} C_{om} - C_{sell} ]其中(C_{grid}) 是向大电网购电的费用(C_{gas}) 是燃气轮机消耗天然气的燃料费用(C_{om}) 是各设备的运行维护费用(C_{sell}) 是微网向大电网售电的收益。这个目标函数的好处是清晰、可解释性强每个费用项都有明确的工程对应关系。在具体实现时购电费用写为[ C_{grid} \sum_{t1}^{24} c_{buy}(t) \cdot P_{buy}(t) ]这里的 (P_{buy}(t)) 是t时段微网从大电网购入的电功率(c_{buy}(t)) 是对应时段的购电价。为了体现分时电价对调度策略的影响我在模型中设置了峰、平、谷三个电价时段峰段电价高谷段电价低。这种设计在现实中有很强的对应关系——电网鼓励用户削峰填谷用户自然有动机在谷段多购电、少购气。售电收益写成[ C_{sell} \sum_{t1}^{24} c_{sell}(t) \cdot P_{sell}(t) ]需要注意的是同一时段购电和售电不能同时发生这需要在约束里用0-1变量强制互斥。燃气费用根据机组出力和效率计算[ C_{gas} \sum_{t1}^{24} \frac{c_{gas} \cdot P_{MT}(t) \cdot \Delta t}{\eta_{MT}(t)} ]这里 (P_{MT}(t)) 是燃气轮机t时段的发电功率(\eta_{MT}(t)) 是发电效率。有些模型更精细会把燃气轮机的热电比也设为变量让它在一个可行区间内连续调整。但在初版模型里我建议先固定热电比把问题做简单跑通后再逐步增加复杂度。3.2 电功率平衡多设备之间的协同电功率平衡是所有微网模型的“压舱石”。热电联供微网的电力平衡方程如下[ P_{MT}(t) P_{PV}(t) P_{WT}(t) P_{buy}(t) P_{dis}(t) P_{load}(t) P_{EB}(t) P_{ch}(t) P_{sell}(t) ]左边是电源总出力燃气轮机发电、光伏出力、风电出力、购电功率、储能放电右边是负荷总需求电负荷、电锅炉耗电、储能充电、售电功率。这个约束的意义在于把“供”和“求”锁死——任何时候系统内的电力都必须保持瞬时的供需平衡。我在调试初期犯过的一个错误是漏掉了电锅炉的耗电项结果模型算出来的“最优解”在现实中根本不可能实现因为电锅炉消耗的那部分电力根本不存在。这类问题的排查不靠代码报错而要靠对物理过程的敏感度。建议拿到结果后先算一遍能量守恒所有电源出力加总等于所有负荷加总差距超过1%就说明模型有bug。3.3 热功率平衡CHP机组、电锅炉与蓄热罐的协调热功率平衡方程是热电联供微网区别于点纯电力微网的核心[ H_{MT}(t) H_{EB}(t) H_{dis}^{TST}(t) H_{load}(t) H_{ch}^{TST}(t) ]这里 (H_{MT}(t)) 是燃气轮机余热回收提供的热功率(H_{EB}(t)) 是电锅炉产生的热功率(H_{dis}^{TST}(t)) 是蓄热罐的放热功率(H_{load}(t)) 是热负荷(H_{ch}^{TST}(t)) 是蓄热罐的充热功率。CHP机组的热功率与电功率通过热电比关联[ H_{MT}(t) \alpha_{CHP} \cdot P_{MT}(t) ](\alpha_{CHP}) 是热电比通常取1.2到1.5之间。这个关系式看上去简单但它揭示了一个关键矛盾燃气轮机在发电的同时必然产热电负荷高峰时段机组出力大产热也多但如果此时热负荷并不高多余的热量就只能靠蓄热罐吸收。蓄热罐装不下了就要限制机组出力这就是“以热定电”的运行模式。反过来如果热负荷高而电负荷低机组为了供热不得不发电多余的电只能低价卖给电网。这种电热耦合的强约束正是热电联供微网优化调度最有意思的地方。3.4 设备约束出力上下限、爬坡和储能动态设备出力上下限是基础约束。燃气轮机、电锅炉的出力必须在最小稳定运行功率和额定功率之间[ P_{MT}^{min} \cdot u_{MT}(t) \leq P_{MT}(t) \leq P_{MT}^{max} \cdot u_{MT}(t) ]这里的 (u_{MT}(t)) 是0-1变量表示机组t时段是否运行。如果机组停机出力必须为0如果运行出力必须落在上下限区间内。这个约束在MILP里的处理非常典型也是初学者最容易出错的地方——直接把上下限写成 (P_{MT}^{min} \leq P_{MT}(t) \leq P_{MT}^{max})那模型就会在停机时段让机组输出一个极小的出力来“钻空子”结果没有任何物理意义。储能电池的动态逻辑用SOC状态变量描述[ SOC(t1) SOC(t) \eta_{ch} \cdot P_{ch}(t) \cdot \Delta t - \frac{P_{dis}(t) \cdot \Delta t}{\eta_{dis}} ]其中 (\eta_{ch}) 和 (\eta_{dis}) 分别是充放电效率SOC要在允许的上下限之间[ SOC^{min} \leq SOC(t) \leq SOC^{max} ]蓄热罐的逻辑与电池类似但充放热效率通常取1假设绝热良好同时要设置最大蓄热量和最大充放热功率。储能和蓄热罐是微网灵活性的“蓄水池”它们的加入让电、热之间的时空耦合变得可以调节也是优化模型能够“玩出花样”的关键。4. Matlab代码实现从数据到结果的关键细节4.1 主程序框架与数据准备一个完整的Matlab实现大致是这个结构%% 初始化清空工作区 clear; clc; close all; %% 参数定义 % 时间参数 T 24; % 调度时段数单位小时 dt 1; % 时间步长单位小时 % 分时电价元/kWh price_buy [0.38*ones(1,7), 0.68*ones(1,5), 1.05*ones(1,4), ... 0.68*ones(1,4), 1.05*ones(1,4)]; % 简化示例 price_sell price_buy * 0.85; % 上网电价通常为购电价的85% % 设备参数 P_MT_max 800; % 燃气轮机最大电功率单位kW P_MT_min 100; % 燃气轮机最小稳定运行电功率单位kW eta_MT 0.35; % 燃气轮机发电效率 alpha_CHP 1.3; % 热电比 P_EB_max 600; % 电锅炉最大功率单位kW eta_EB 0.95; % 电锅炉热效率 % 储能参数 E_bat_max 800; % 电池容量单位kWh SOC_init 0.5; % 初始SOC SOC_min 0.2; SOC_max 0.9; P_bat_max 200; % 最大充放电功率单位kW eta_ch 0.95; % 充电效率 eta_dis 0.95; % 放电效率 % 蓄热罐参数 H_TST_max 1000; % 最大蓄热量单位kWh H_TST_init 500; % 初始蓄热量 H_TST_min 100; H_TST_max_power 300; % 最大充放热功率单位kW %% 负荷与新能源出力数据24小时 P_load [520, 480, 450, 430, 420, 450, 500, 580, 650, 720, 780, ... 800, 750, 700, 650, 620, 680, 750, 820, 860, 820, 750, ... 680, 600]; % 电负荷单位kW H_load [600, 580, 560, 540, 520, 500, 480, 450, 400, 380, 360, ... 350, 340, 350, 370, 400, 450, 500, 550, 580, 600, 620, ... 610, 600]; % 热负荷单位kW P_PV [0, 0, 0, 0, 0, 30, 120, 250, 380, 450, 500, 520, ... 500, 450, 370, 250, 130, 40, 0, 0, 0, 0, 0, 0]; % 光伏出力 P_WT [180, 170, 160, 150, 140, 130, 120, 150, 200, 240, 260, 250, ... 240, 230, 220, 210, 200, 190, 180, 170, 160, 150, 170, 180]; % 风电出力这段代码定义了模型所需的全部基础数据。在实际项目中这些数据来源各不相同电价从电网公司获取负荷数据来自历史量测光伏和风电出力来自预测系统或典型日曲线。我在这里用人工设定的数据来保证示例的完整可跑性读者可以替换为自己的真实数据。4.2 决策变量定义与约束构建接下来是定义决策变量。强烈建议直接在Yalmip中定义变量不要手动构造向量再传给求解器因为Yalmip的符号化建模能极大减少索引错误的概率。%% 定义决策变量 % 连续变量 P_MT sdpvar(1, T); % 燃气轮机发电功率 H_MT sdpvar(1, T); % 燃气轮机余热回收功率 P_EB sdpvar(1, T); % 电锅炉耗电功率 H_EB sdpvar(1, T); % 电锅炉热功率 P_buy sdpvar(1, T); % 购电功率 P_sell sdpvar(1, T); % 售电功率 P_bat_ch sdpvar(1, T); % 电池充电功率 P_bat_dis sdpvar(1, T); % 电池放电功率 SOC sdpvar(1, T1); % 电池SOC状态 H_TST_ch sdpvar(1, T); % 蓄热罐充热功率 H_TST_dis sdpvar(1, T); % 蓄热罐放热功率 H_TST_level sdpvar(1, T1); % 蓄热罐蓄热量 % 0-1整数变量 u_MT binvar(1, T); % 燃气轮机启停状态1为运行 u_buy binvar(1, T); % 购电状态1为购电 u_bat_ch binvar(1, T); % 电池充电状态1为充电 u_bat_dis binvar(1, T); % 电池放电状态1为放电变量定义好后按前面3.2到3.4节的约束逐步添加。Yalmip的优势在于约束写法和数学公式几乎一一对应%% 约束条件 constraints []; % 电功率平衡约束 for t 1:T constraints [constraints, ... P_MT(t) P_PV(t) P_WT(t) P_buy(t) P_bat_dis(t) ... P_load(t) P_EB(t) P_bat_ch(t) P_sell(t)]; end % 热功率平衡约束 for t 1:T constraints [constraints, ... H_MT(t) H_EB(t) H_TST_dis(t) ... H_load(t) H_TST_ch(t)]; end % 燃气轮机约束出力上下限 热电比关联 for t 1:T constraints [constraints, ... P_MT_min * u_MT(t) P_MT(t) P_MT_max * u_MT(t)]; constraints [constraints, ... H_MT(t) alpha_CHP * P_MT(t)]; end % 电锅炉约束 for t 1:T constraints [constraints, ... 0 P_EB(t) P_EB_max]; constraints [constraints, ... H_EB(t) eta_EB * P_EB(t)]; end % 储电约束SOC动态 充放电互斥 for t 1:T constraints [constraints, ... SOC(t1) SOC(t) eta_ch * P_bat_ch(t) * dt - P_bat_dis(t) * dt / eta_dis]; constraints [constraints, ... SOC_min SOC(t) SOC_max]; constraints [constraints, ... 0 P_bat_ch(t) P_bat_max * u_bat_ch(t)]; constraints [constraints, ... 0 P_bat_dis(t) P_bat_max * u_bat_dis(t)]; constraints [constraints, ... u_bat_ch(t) u_bat_dis(t) 1]; end constraints [constraints, SOC(1) SOC_init, SOC(T1) SOC_init]; % 购售电互斥约束 for t 1:T constraints [constraints, ... 0 P_buy(t) 1500 * u_buy(t)]; constraints [constraints, ... 0 P_sell(t) 800 * (1 - u_buy(t))]; end % 蓄热罐约束 for t 1:T constraints [constraints, ... H_TST_level(t1) H_TST_level(t) H_TST_ch(t) * dt - H_TST_dis(t) * dt]; constraints [constraints, ... H_TST_min H_TST_level(t) H_TST_max]; constraints [constraints, ... 0 H_TST_ch(t) H_TST_max_power]; constraints [constraints, ... 0 H_TST_dis(t) H_TST_max_power]; end constraints [constraints, ... H_TST_level(1) H_TST_init, ... H_TST_level(T1) H_TST_init];4.3 目标函数建立与求解器配置目标函数的Matlab实现如下%% 目标函数单位元 objective sum(price_buy .* P_buy) * dt ... % 购电费用 sum(P_MT * c_gas / eta_MT) * dt ... % 燃气费用 sum(0.02 * P_MT 0.01 * P_EB 0.01 * P_bat_ch) * dt ... % 运维费用 - sum(price_sell .* P_sell) * dt; % 售电收益其中 (c_{gas}) 是天然气价格单位为元/kWh折算到天然气热值。求解器配置如下%% 求解 ops sdpsettings(solver, cplex, verbose, 1, ... showprogress, 1, debug, 1); optimize(constraints, objective, ops); %% 结果提取 P_MT_opt value(P_MT); H_MT_opt value(H_MT); P_EB_opt value(P_EB); P_buy_opt value(P_buy); P_sell_opt value(P_sell); ...注意如果电脑里没有Cplex或Gurobi也可以先用Matlab自带的intlinprog求解只需要在Yalmip的sdpsettings里把solver设为intlinprog。不过对于规模较大的问题商业求解器的速度和稳定性明显更好。Cplex和Gurobi都对学术用户提供免费授权申请流程不复杂建议尽早配好。5. 结果分析与运行策略解读5.1 调度结果图怎么看求解完成后把变量结果画成图。我通常用两张图一张是电功率平衡堆叠图一张是热功率平衡堆叠图再加一张SOC和蓄热罐蓄热量的变化曲线。这三张图能反映调度策略的绝大部分信息。电功率平衡图能看出几个关键问题燃气轮机在哪些时段满载运行、哪些时段停机储能什么时候充电、什么时候放电购电和售电的切换节点在哪里。热功率平衡图则能看出CHP余热在什么时候供给热负荷、什么时候充入蓄热罐、蓄热罐又在何时放热。我跑完这个示例后最典型的策略特征包括谷电时段电价低优先购电燃气轮机降出力运行储能充电蓄热峰电时段电价高燃气轮机满发余热优先供热不足部分由蓄热罐放热补充如果还有富余电就卖给电网。整体上优化结果会自动形成一个“削峰填谷”的调度方案。5.2 算例结果成本分布与运行特征以下是我用上述参数跑出来的一个典型结果单位元费用项数值占比购电费用486242.1%燃气费用534646.3%设备运维费用3983.4%售电收益-943-8.2%总运行成本9663100%从中可以看出几个有意思的现象购电费和燃气费是大头各占四成多说明系统主要依赖电网和气网两种外部能量来源售电收益虽然占比不高但它为系统提供了“源荷互动”的灵活性让燃气轮机在热负荷驱动下被迫多发的那部分电力有了经济出口。在设备运行特征方面燃气轮机在峰电时段10:00-13:00、18:00-21:00基本满发平电时段按需调整谷电时段则直接停机用购电代替自发。储能的充放策略也很有意思它不只是在谷充峰放还会在负荷最低的凌晨时段充一部分电在早上负荷快速爬升时段放出起到“削峰”作用。5.3 灵敏度分析参数变化对调度策略的影响模型跑通后建议做一组灵敏度分析这会让你对系统的运行特点有更深入的理解。我常用的做法是分别改变分时电价、天然气价格、热电比、储能容量观察总成本和各设备出力的变化。以天然气价格为例当气价从0.25元/kWh升到0.35元/kWh时燃气轮机的调度策略会从“电热联供为主”逐步转向“以热定电为主”即只在必须满足热负荷时才发电多余的电力尽量多卖给电网以摊薄成本。这就解释了为什么天然气价格是热电联供微网经济性的最敏感因素之一——它直接影响CHP机组的启停和出力区间。6. 常见问题与调试经验6.1 求解器报错“Infeasible Problem”怎么办这是MILP建模中最常见的问题几乎每个初学者都会遇到。模型无可行解意味着你写的约束之间存在数学上的矛盾。我总结了一套标准排查流程。第一步检查能量平衡方程是否写漏项。这是最常见的原因。比如电功率平衡里漏了电锅炉耗电热功率平衡里漏了蓄热罐充放热都会导致系统在物理上不可能满足所有约束。我建议把每个时段的平衡式单独打印出来逐步核对。第二步检查变量边界是否合理。例如储能SOC初值设定后在给定最大充放电功率下能否在一个时段内将电量充到目标值如果 (SOC_{init}0.2)最大充电功率下也只能在一小时内到达 (SOC0.38)但你约束了 (SOC(2)\geq 0.5)那这个模型就无解。第三步检查0-1变量关联是否正确。购售电互斥、充放电互斥这类约束如果二进制变量的使用方式有误容易造成隐含的不可行。比如u_buy(t) u_sell(t) 1里的两个变量如果定义反了就会导致某一时段既不能购电也不能售电如果此时系统又恰好需要外部电力支撑无解就不奇怪了。6.2 模型求解缓慢或者间隙降不下去如果模型规模不大但求解时间特别长通常不是求解器的问题而是建模方式的问题。一个最常见的坏习惯是大量使用Big-M形式的约束即引入很大的常数M把非线性条件线性化。M值如果取得过大会导致求解器的线性松弛非常松散分支定界的搜索空间大幅膨胀。经验法则是M值尽量取到约束实际边界的1.1~1.2倍不要为了“保险”取一个10的6次方级别的数值。另一个提升求解速度的技巧是给求解器设置合理的最优间隙。对于工程应用Gap控制在1%以内就足够好了。在Cplex里可以用ops sdpsettings(solver, cplex, cplex.mip.tolerances.mipgap, 0.01);这样能显著缩短求解时间而结果的工程误差几乎可以忽略。6.3 Yalmip版本与求解器兼容性Yalmip是一个持续更新的工具箱不同版本对求解器的接口略有差异。如果你遇到“No suitable solver found”的报错多半是Yalmip没找到对应的求解器或者求解器没有正确添加到Matlab路径。我的建议是装好求解器之后先在命令行跑一次yalmiptest这个命令会列出所有已检测到的求解器和状态。如果Cplex或Gurobi显示“Found”但状态不是“OK”检查一下求解器的license是否有效路径是否已添加。6.4 从学术模型走向工程应用的注意事项最后说一个容易被忽视的问题学术模型和工程应用之间有一条不小的鸿沟。学术模型里常常假定预测完美、设备无故障、出力连续可调但实际运行中这些假设都可能被打破。我自己做项目时吃过亏——模型在理论数据上跑出来“最优”调度但在现场数据下偏差很大。这提醒我在模型设计阶段就要为未来扩展留出空间比如将预测不确定性的场景集纳入约束或者用鲁棒优化/随机优化的框架替代确定性模型。这些都可以在MILP基础上逐步扩展而不需要推翻重来。在整个热电联供微网优化运行的开发过程中我最深的体会是数学模型是骨架工程经验是血肉。没有模型工程决策就是拍脑袋没有工程视角模型就只是纸上谈兵。Matlab这套实现的价值在于它把两者结合起来让你既能从理论层面理解多能互补的逻辑又能把理论落地成可执行的调度方案。如果你是第一次接触这个方向建议用我给的代码跑通一个最简单的基础案例然后逐步增加可再生能源渗透率、储能容量、碳交易机制等因素感受不同参数对调度策略的影响。这个循序渐进的过程会比直接照搬一篇论文的完整模型有用得多。