简介这份资源面向能源系统优化、综合能源调度方向的研究生与科研人员围绕电-气-热综合能源系统的耦合建模与优化调度展开适合具备一定MATLAB基础、希望复现多能互补调度案例的读者。压缩包共28个文件约5.13MB以20个m脚本为核心配套xls、xlsx与txt数据说明、docx文档及md说明文件覆盖模型构建、耦合机制、优化求解与仿真分析等环节。内容涉及电力、燃气、热力三个子系统的数学建模能量转换与供需平衡的耦合关系以及以运行成本、碳排放为目标的优化调度实现并借助MATLAB优化工具箱与SIMULINK完成算法验证与动态仿真。已有151人学习可帮助读者快速理解多能源系统联合调度的建模思路、约束设置与求解流程为课题复现与二次开发提供参考。1. 电-气-热综合能源系统耦合调度从三张网到一套求解器园区能源站里最头疼的事往往不是某台设备效率不够而是电、气、热三张网各算各的。电侧按分时电价排产气侧按管网压力调阀热侧按回水温度启停锅炉三套逻辑互不通气结果就是电低谷时气价高、热负荷上来时电又在峰值。电-气-热综合能源系统耦合调度要解决的正是把这三类能流放进同一个优化模型里让电转气、燃气轮机、余热锅炉、热泵这些耦合设备在时间维度上协同出力。它适合做园区级能源管理、微网调度、多能互补项目的人也适合想把单一能流优化扩展到多能流的研究者。核心难点不在建模本身而在耦合矩阵的构造和求解器对非线性约束的处理。2. 耦合调度的数学骨架从能流平衡到矩阵组装2.1 三网耦合的物理关系怎么写成约束电-气-热综合能源系统的耦合本质是能量转换设备把一种能流变成另一种。燃气轮机消耗天然气产生电和热电转气设备消耗电产生天然气热泵消耗电产生热余热锅炉回收燃气轮机排烟产生热。每条转换关系对应一组等式约束写进优化模型时通常用效率系数和转换比来表达。以燃气轮机为例输入天然气流量为 $F_{gt}$输出电功率 $P_{gt}$ 和热功率 $H_{gt}$关系可以写成$$P_{gt} \eta_{gt}^e \cdot F_{gt} \cdot LHV$$ $$H_{gt} \eta_{gt}^h \cdot F_{gt} \cdot LHV$$其中 $LHV$ 是天然气低位热值$\eta_{gt}^e$ 和 $\eta_{gt}^h$ 分别是发电效率和余热回收效率。这两个等式把气网节点和电网节点、热网节点绑在一起是耦合矩阵里最核心的非零块。电转气的方向相反消耗电功率 $P_{p2g}$ 产生天然气流量 $F_{p2g}$$$F_{p2g} \eta_{p2g} \cdot P_{p2g} / LHV$$热泵则连接电和热$$H_{hp} COP_{hp} \cdot P_{hp}$$这些等式在建模时不能只写符号必须明确每个变量的单位和量纲。我见过不少模型因为把 kW 和 MW 混用求解结果偏差三个数量级排查半天才发现是单位问题。2.2 目标函数里电价、气价、热价怎么加权目标函数通常是运行成本最小化包含购电成本、购气成本、设备运维成本和弃风弃光惩罚。写法如下# 目标函数总运行成本最小化 # 变量说明 # P_buy[t] : t时段购电功率 (kW) # F_buy[t] : t时段购气流量 (m3/h) # price_e[t] : t时段电价 (元/kWh) # price_g[t] : t时段气价 (元/m3) # C_om[i][t] : 设备i在t时段的运维成本 (元) # C_cur[t] : t时段弃风弃光惩罚 (元) def objective(model): cost_e sum(price_e[t] * P_buy[t] * dt for t in model.T) cost_g sum(price_g[t] * F_buy[t] * dt for t in model.T) cost_om sum(C_om[i][t] for i in model.devices for t in model.T) cost_cur sum(lambda_cur * P_cur[t] * dt for t in model.T) return cost_e cost_g cost_om cost_cur这里的关键参数是dt即时间步长。园区调度常用 15 分钟或 1 小时步长越小精度越高但变量数成倍增长。lambda_cur是弃风弃光惩罚系数设得太小会导致优化结果大量弃风设得太大又可能让求解器找不到可行解。我一般先设 0.5 元/kWh 试算再根据实际弃风率调整。2.3 用 Python 搭一个最小可跑的三网耦合模型下面是一个简化版的三网耦合调度模型用 Pyomo 建模调用开源求解器 CBC 求解。模型包含电平衡、气平衡、热平衡和三个耦合设备。import pyomo.environ as pyo # 时间步长1小时共24个时段 T range(24) dt 1.0 model pyo.ConcreteModel() model.T pyo.Set(initializeT) # 决策变量 model.P_gt pyo.Var(model.T, domainpyo.NonNegativeReals) # 燃气轮机发电 model.H_gt pyo.Var(model.T, domainpyo.NonNegativeReals) # 燃气轮机产热 model.F_gt pyo.Var(model.T, domainpyo.NonNegativeReals) # 燃气轮机耗气 model.P_p2g pyo.Var(model.T, domainpyo.NonNegativeReals) # 电转气耗电 model.F_p2g pyo.Var(model.T, domainpyo.NonNegativeReals) # 电转气产气 model.H_hp pyo.Var(model.T, domainpyo.NonNegativeReals) # 热泵产热 model.P_hp pyo.Var(model.T, domainpyo.NonNegativeReals) # 热泵耗电 model.P_buy pyo.Var(model.T, domainpyo.NonNegativeReals) # 购电 model.F_buy pyo.Var(model.T, domainpyo.NonNegativeReals) # 购气 # 参数示例值实际项目需按设备铭牌填写 eta_gt_e 0.35 # 燃气轮机发电效率 eta_gt_h 0.45 # 余热回收效率 eta_p2g 0.60 # 电转气效率 COP_hp 3.5 # 热泵能效比 LHV 9.7 # 天然气低位热值 kWh/m3 # 负荷曲线示例24小时电、气、热负荷 P_load [500,480,460,450,470,520,600,700,750,780,800,820, 810,790,770,760,780,820,850,830,780,700,620,560] F_load [80,75,70,68,72,85,100,120,130,135,140,142, 138,132,128,126,130,138,145,140,130,115,100,90] H_load [200,190,180,175,180,210,260,320,360,380,400,410, 405,395,385,380,390,410,430,420,390,350,300,250] # 电平衡约束 def power_balance(m, t): return m.P_gt[t] m.P_buy[t] m.P_p2g[t] * 0 m.P_hp[t] * 0 \ P_load[t] m.P_p2g[t] m.P_hp[t] model.power_balance pyo.Constraint(model.T, rulepower_balance) # 气平衡约束 def gas_balance(m, t): return m.F_gt[t] m.F_buy[t] m.F_p2g[t] F_load[t] model.gas_balance pyo.Constraint(model.T, rulegas_balance) # 热平衡约束 def heat_balance(m, t): return m.H_gt[t] m.H_hp[t] H_load[t] model.heat_balance pyo.Constraint(model.T, ruleheat_balance) # 耦合设备等式约束 def gt_power(m, t): return m.P_gt[t] eta_gt_e * m.F_gt[t] * LHV model.gt_power pyo.Constraint(model.T, rulegt_power) def gt_heat(m, t): return m.H_gt[t] eta_gt_h * m.F_gt[t] * LHV model.gt_heat pyo.Constraint(model.T, rulegt_heat) def p2g_gas(m, t): return m.F_p2g[t] eta_p2g * m.P_p2g[t] / LHV model.p2g_gas pyo.Constraint(model.T, rulep2g_gas) def hp_heat(m, t): return m.H_hp[t] COP_hp * m.P_hp[t] model.hp_heat pyo.Constraint(model.T, rulehp_heat) # 目标函数购电购气成本 price_e [0.35,0.35,0.35,0.35,0.35,0.35,0.60,0.60,0.60,0.90,0.90,0.90, 0.90,0.90,0.90,0.60,0.60,0.60,0.90,0.90,0.90,0.60,0.35,0.35] price_g [2.5]*24 def obj(m): return sum(price_e[t]*m.P_buy[t]*dt price_g[t]*m.F_buy[t]*dt for t in m.T) model.obj pyo.Objective(ruleobj, sensepyo.minimize) # 求解 solver pyo.SolverFactory(cbc) results solver.solve(model, teeTrue) print(f总运行成本: {pyo.value(model.obj):.2f} 元)这段代码的逻辑是先定义 24 个时段的决策变量再写三网平衡约束和耦合设备等式最后用分时电价和气价构成目标函数。求解器 CBC 是开源混合整数规划求解器适合中小规模问题。如果模型里加入机组启停状态变量就变成混合整数线性规划CBC 也能处理但规模大了会慢。参数方面eta_gt_e和eta_gt_h要按实际燃气轮机型号填常见小型燃气轮机发电效率在 0.25 到 0.35 之间余热回收效率在 0.4 到 0.5 之间。COP_hp热泵能效比在 3.0 到 4.5 之间低温工况取低值。LHV天然气低位热值一般取 9.7 kWh/m³ 左右不同气源有差异。注意这个模型没有考虑管网动态和储能实际项目里热网惯性往往不能忽略后面章节会讲怎么处理。3. 求解器选型与参数调优为什么你的模型跑不动3.1 线性、非线性、混合整数三类求解器怎么选电-气-热耦合调度模型按复杂度分三档。第一档是纯线性模型耦合设备效率设为常数所有约束都是线性的用 CBC、GLPK、HiGHS 这类开源求解器就能跑变量数几千以内秒级出结果。第二档是混合整数线性模型加入机组启停、储能充放状态等 0-1 变量CBC 和 HiGHS 也能处理但变量数上万后求解时间可能到分钟级。第三档是非线性模型比如考虑燃气轮机效率随负荷率变化、热泵 COP 随环境温度变化这时需要 IPOPT、SCIP 这类非线性求解器或者用分段线性化把非线性转成线性。我一般先用线性模型跑通确认约束逻辑没问题再逐步加入非线性因素。直接上非线性模型容易在调试阶段就卡住因为非线性求解器对初值敏感约束写错一点就报不可行排查起来很痛苦。3.2 时间步长、预测误差与鲁棒性之间的取舍时间步长决定了模型的精度和规模。1 小时步长适合日前调度15 分钟步长适合日内滚动5 分钟步长适合实时控制。步长减半变量数翻倍求解时间大约增加 2 到 3 倍。园区级调度一般用 1 小时做日前计划15 分钟做日内修正。预测误差是另一个关键参数。电负荷预测误差通常在 3% 到 5%热负荷预测误差在 5% 到 8%光伏出力预测误差在 10% 到 20%。如果模型里用确定性预测值实际运行时会出现功率缺额。常见做法是加旋转备用约束或者用鲁棒优化把预测误差范围写进约束。# 鲁棒约束示例电负荷预测误差 ±5% epsilon_e 0.05 def power_balance_robust(m, t): return m.P_gt[t] m.P_buy[t] P_load[t] * (1 epsilon_e) model.power_balance_robust pyo.Constraint(model.T, rulepower_balance_robust)这个约束保证即使负荷比预测高 5%系统也有足够出力。代价是运行成本会上升因为备用容量平时也在消耗成本。epsilon_e设多大取决于预测系统的实际精度设太大经济性差设太小风险高。3.3 用滚动时域把日前调度和实时修正串起来日前调度给出 24 小时计划但实际运行中负荷和新能源出力会偏离预测。滚动时域的做法是每 15 分钟用最新预测数据重新求解未来 4 小时的调度计划只执行第一个时段的指令下一个周期再滚动。# 滚动时域调度框架 horizon 4 # 预测窗口4小时 step 0.25 # 滚动步长15分钟 total_steps int(24 / step) for k in range(total_steps): # 获取最新预测数据 P_load_pred get_load_forecast(k, horizon) F_load_pred get_gas_forecast(k, horizon) H_load_pred get_heat_forecast(k, horizon) # 构建并求解滚动窗口内的优化模型 model build_model(P_load_pred, F_load_pred, H_load_pred) result solver.solve(model) # 只执行第一个时段的指令 execute_command(result, k)这个框架的关键是get_load_forecast返回的预测数据要随 k 更新。实际项目中预测系统可能每 15 分钟更新一次滚动调度就跟着更新。如果预测系统更新慢滚动步长可以放宽到 30 分钟或 1 小时。注意滚动时域会增加计算频次如果单次求解超过滚动步长就会导致指令堆积。我一般会限制单次求解时间超时就用上一周期的解做退化处理。4. 避坑与排查耦合调度模型最常见的五类翻车4.1 求解器报 infeasible 但约束看起来没问题现象模型写完后求解器返回 infeasible逐条检查约束都合理。原因最常见的是耦合设备容量约束和负荷需求矛盾。比如热负荷 400 kW但燃气轮机最大产热 300 kW热泵最大产热 80 kW加起来不够模型自然不可行。另一个原因是储能荷电状态约束写反了导致能量无法平衡。解决先去掉目标函数只保留约束用求解器的不可行诊断功能找冲突约束。Pyomo 可以调用pyo.SolverFactory(cbc).solve(model, teeTrue)查看日志或者用model.write(model.lp)导出 LP 文件手动检查。我一般会先把负荷乘以 0.8 试算如果可行就说明是容量不够再逐项加回负荷找临界点。4.2 结果里购电和购气同时很大耦合设备却不出力现象优化结果里购电和购气都接近上限但燃气轮机和电转气出力接近零。原因目标函数里耦合设备的运维成本设得太高或者耦合设备效率设得太低导致优化器认为直接购能更便宜。另一个可能是耦合设备约束里变量方向写反了比如把P_gt eta * F_gt * LHV写成了F_gt eta * P_gt * LHV导致量纲不匹配。解决检查耦合设备等式约束的变量方向确认效率系数乘在正确的一侧。然后检查运维成本参数燃气轮机运维成本一般在 0.05 到 0.15 元/kWh电转气在 0.1 到 0.3 元/kWh热泵在 0.02 到 0.08 元/kWh。如果设得比购能价差还大耦合设备就不会出力。4.3 热网惯性没考虑导致热平衡约束频繁越限现象日前调度计划里热平衡满足但实际运行时热网温度波动大热平衡约束频繁越限。原因模型里把热网当成瞬时平衡忽略了管道和建筑的热惯性。实际热网从热源到用户有十几分钟到几小时的延迟热源出力变化后用户侧不会立刻响应。解决在热平衡约束里加入热惯性项用一阶惯性环节描述# 热惯性约束热网温度动态 tau_h 0.5 # 热惯性时间常数小时 def heat_inertia(m, t): if t 0: return pyo.Constraint.Skip return m.H_net[t] m.H_net[t-1] (m.H_gt[t] m.H_hp[t] - H_load[t]) * dt / tau_htau_h根据实际管网长度和建筑热容确定一般园区级热网在 0.3 到 1.0 小时之间。加入惯性后模型从代数方程变成差分方程变量维度增加但更贴近实际。4.4 分时电价时段划分和实际账单对不上现象优化结果显示在电价低谷时段大量购电但实际电费账单并没有明显下降。原因分时电价时段划分和实际执行的不一致。很多地区的峰谷时段不是简单的早中晚三段而是有尖峰、高峰、平段、低谷四段甚至五段而且季节不同时段也不同。解决从实际电费账单里反推分时电价曲线或者从当地能源管理部门获取最新分时电价文件。把 24 小时的电价数组按实际时段填写不要用理想化的三段式。我一般会保留最近三个月的电费账单逐月核对电价曲线。4.5 模型跑通了但结果不收敛到全局最优现象求解器返回成功但结果明显不是最优比如明明有更便宜的购能组合却选了贵的。原因非线性模型陷入局部最优或者混合整数模型的 gap 设得太大。CBC 默认 gap 是 0.01%但有些版本默认 1%导致求解器在离最优解还有 1% 时就停止。解决检查求解器的 gap 设置把mipgap设到 0.001 或更小。如果是非线性模型尝试多组初值求解取最好的结果。Pyomo 里可以设置solver.options[mipgap] 0.001。另外如果模型里有大量 0-1 变量可以加割平面或启发式初值加速收敛。5. 从能流平衡到碳流追踪耦合调度还能怎么用碳流追踪是这两年多能流调度里比较实用的扩展。电-气-热耦合系统里每个设备的碳排放因子不同电网购电的碳因子随电源结构变化天然气碳因子相对固定热泵如果消耗绿电则碳因子接近零。把碳流写成约束或目标函数就能在调度时兼顾经济性和低碳性。具体做法是在目标函数里加入碳成本项# 碳成本项 def carbon_cost(model): carbon_e sum(cf_e[t] * P_buy[t] * dt for t in model.T) # 购电间接碳排 carbon_g sum(cf_g * F_buy[t] * dt for t in model.T) # 购气直接碳排 carbon_total carbon_e carbon_g return carbon_total * carbon_pricecf_e[t]是 t 时段电网购电的碳因子单位 kgCO2/kWh随电网电源结构变化。cf_g是天然气碳因子约 0.2 kgCO2/kWh。carbon_price是碳价单位元/kgCO2。把碳成本加到运行成本里优化器就会在电价低但碳因子高的时段减少购电转向低碳设备。验证碳流追踪结果是否合理可以看两个指标一是总碳排是否比纯经济调度低二是碳排降低的幅度是否和碳价设置匹配。如果碳价设了 0.1 元/kgCO2 但碳排几乎没降说明碳成本在总成本里占比太小需要提高碳价或检查碳因子是否填错。我自己的习惯是每次改完模型参数后先跑一遍纯经济调度做基准再跑一遍加碳成本的调度对比两者的购能结构和碳排。如果碳排降了但成本涨了超过 15%就要重新审视碳价是否设得过高。这个对比过程比单看一个结果更能说明模型是否在按预期工作。希望帮到你。本文还有配套的精品资源点击获取