做综合能源系统调度有一段时间了电网侧的机组组合、经济调度写过不少但第一次摊上电-气-热三个网络耦合的优化调度时还是被结结实实地上了一课。独立的电网模型、独立的气网模型、独立的热网模型单独拎出来都能跑一旦把它们嵌进同一个优化框架变量数量、约束耦合和非线性程度完全不是一个量级。这篇文章就把我搭这套电-气-热综合能源系统耦合优化调度Matlab代码的完整过程拆开讲从数学建模到求解器选型从代码框架到实际跑算例踩过的坑一次性说清楚适合正在做综合能源系统、能源枢纽调度或者相关毕业设计的同学参考。这里先明确一下耦合到底耦合在哪。传统各干各的调度里电网和天然气网唯一的联系是燃气机组和天然气价格热网更是基本不参与。所谓耦合优化调度本质是把电力、天然气、热力这三个载体通过转换设备CHP热电联产机组、燃气锅炉、电锅炉、热泵、P2G电转气设备绑在同一个优化问题里让整体运行成本更低、可再生能源消纳空间更大。下面从建模到代码一步步说。1. 电-气-热耦合调度解决的问题三个网络的交叉点在哪1.1 单能调度的局限性为什么非要把三个网放一起算先看一个最简单的案例。一个区域里有风电、火电、燃气机组、CHP机组同时有热负荷。如果按照传统的电力调度方式只优化电功率热负荷就交给固定的CHP热出力或者锅炉来满足。问题就出在这里CHP机组电功率和热功率是强耦合的要满足热负荷就不得不维持一定的电出力哪怕这时候风电大发、电力市场电价很低。这就是所谓的以热定电约束是热电耦合调度里最典型的矛盾。加上天然气网之后情况更复杂。燃气机组和燃气锅炉都要从气网买气气网节点压力、管道流量和压缩机状态会影响天然气的输送而天然气价格又反过来决定燃气机组的发电成本。电、气、热三个网络通过设备和能源价格形成了完整的耦合链。单独做电力调度完全看不到这些约束的联动效应优化结果在工程上根本不可行。所以电-气-热耦合调度的第一个核心价值就是让以电定热和以热定电这种矛盾在一个统一的模型里被真正解决通过锅炉补热、电锅炉调峰、储气储热等手段找到全局最优的运行点。1.2 耦合设备是解耦的关键CHP、P2G、电锅炉的角色耦合优化调度之所以能联全靠中间的转换设备。这些设备的特性曲线就是不同能源载体之间的桥梁约束。最核心的几个CHP热电联产机组同时输出电和热电热出力关系一般用一个可行域来表示。电定热型和热定电型机组的可行域不同建模时可以直接用线性不等式组来近似这是后文会详细讲的。燃气锅炉与燃气轮机从气网取天然气转化为热或电是电力系统和天然气系统的直接交点。电锅炉与热泵用电力制热实现电-热耦合。尤其是可再生能源富余时段用电锅炉消纳弃风电给热网供热是综合能源系统调度的经典场景。P2G电转气Power-to-Gas用电解水制氢、再经甲烷化合成天然气注入气网。P2G让能量流动出现了从电到气的逆向路径这是电-气双向耦合的典型设备建模时通常描述为耗电产气、产出气量进入气网节点。理解了设备的耦合关系才谈得上建模。没有这些转换设备做桥梁三个网络就是三个并行的孤立优化问题耦合二字无从谈起。这部分是整个模型设计的起点也是代码里约束构建的核心逻辑。2. 各网络数学模型拆解电网、气网、热网如何写成约束2.1 电网模型直流潮流与节点功率平衡电力网络在日前调度中一般不做交流潮流太多人一开始就栽在这里。完整的交流潮流节点电压、无功、相角都是非线性再加上电网气网耦合求解器基本跑不动。工程上采用直流潮流近似忽略无功和电压幅值只用相角作为状态变量。对节点注入有功功率基准值有节点功率平衡节点注入有功 该节点连接的发电设备出力之和 - 电负荷线路潮流线路有功潮流近似等于两端相角差除以线路电抗即相关公式直接写表达式即可Yalmip里就是一行约束F_ij (theta_i - theta_j) / X_ij发电机约束出力上下限、爬坡速率参考节点平衡节点相角设为0直流潮流模型的精度在日前调度完全够用好处是全部线性不会给后续优化问题带来非线性负担。电网这部分是整个耦合模型里最容易写的绝大多数变量都是连续量如果没有机组启停的0-1变量甚至不涉及整数规划。2.2 天然气网模型管道流量方程与非线性的处理天然气网络的建模和电网差别很大。气网有气源、管道、压缩机、负荷节点燃气机组、燃气锅炉、居民用气状态变量是节点气压和管道流量。关键方程是管道稳态流方程Weymouth方程它把管道流量和两端气压联系起来。简化形式是管道流量平方等于两端气压平方之差乘以一个系数。这个方程是典型的二次非线性等式直接丢给求解器会被按非凸二次约束处理Cplex和Gurobi对这类约束的求解能力非常有限几乎无法在可接受时间内收敛。实际工程里常用的处理方式有三种分段线性化把气压平方或流量的二次关系在若干区间内用线性段逼近缺点是分段数越多越准但变量也越多。增量线性化和分段线性化类似引入连续分段变量精度可控是代码里最常见的实现方式。二阶锥松弛SOCP把等式松弛成不等式变成管道流量的平方小于等于两端气压平方之差乘以系数在某些条件下是严格可解的。但松弛解不满足原等式约束时需要通过迭代收紧来逼近代码复杂度明显上升。我在这套代码里采取的是分段线性化加增量变量的方式因为要保证模型是混合整数线性规划这样Cplex和Gurobi都能稳定求解不会出现非线性求解器的收敛问题。压缩机模型更麻烦压缩机消耗一部分天然气来增压这个耗气量和压缩比、流量耦合在一起通常简化为耗气量占流通量的一个固定比例或随压力比分段线性变化否则模型会变得异常复杂。2.3 热网模型热水网络与热平衡近似热网比气网还特殊。热力系统的主要状态是水温和流量完整的动态热网模型要考虑管道传输延迟这部分我在代码里做了简化——因为做的是日前优化调度时间尺度是小时级到日级热网稳态近似是合理的。热力平衡包括三个层面热负荷平衡热负荷由CHP供热、燃气锅炉供热、电锅炉供热共同满足节点回水/供水温度约束在稳态近似下简化为热功率的节点平衡热源供热能力各热源的供热上下限和爬坡约束实际需要重点注意的是热负荷曲线通常比电负荷曲线更平稳峰值集中在早晚两个时段这决定了耦合调度时热源设备的出力节奏。我在模型里直接用热功率平衡方程不引入温度变量因为温度变量会带来一组额外的等式约束对优化目标和耦合关系的影响很小但会让约束矩阵规模明显膨胀。这个取舍在工程上是划算的。2.4 时间尺度与调度场景为什么用日前调度做统一框架电、气、热三个网络的时间常数完全不同电力系统响应是秒级-分钟级气网是分钟级-小时级热网则是小时级。直接把它们放进同一个动态模型数值刚性问题会把你逼疯。我在代码里用的统一框架是稳态日前调度即把整个调度周期划分为24个时段每个时段内三个网络都按稳态处理。这是综合能源系统经济调度最通用的做法既能捕捉到日内的负荷波动和电价波动又不至于陷入动态模型的巨大计算量。3. 耦合优化调度数学规划建模目标函数、约束与求解器选型3.1 目标函数设计成本最小还是多目标耦合调度代码里最经典的目标函数是系统总运行成本最小。包括这几个部分电网购电成本从上级电网购电电价乘以购电量天然气购买成本从气源点购气量乘以气价机组启停成本如果有启停决策变量弃风弃光惩罚成本可选用来反映可再生能源消纳我一般把目标写成这个样子表述为公式说明 运行成本 电网购电费用 气源购气费用 启停成本 弃风惩罚有人会加碳排放成本做法是在目标函数里加入火电、燃气机组、锅炉的碳排放量乘以碳价。这样目标函数就从单目标变成含碳的多目标但代码结构不变只需要在成本矩阵对应位置加上排放系数。对于一篇讲耦合优化调度代码的博文建议先跑通基本成本模型再叠加上碳约束逻辑上更不容易出错。3.2 决策变量梳理哪些是连续、哪些是整数做变量定义之前一定要先把变量分类想清楚这是MatlabYalmip代码的第一道分水岭。我在这套代码里把变量分成三类电力侧变量各节点相角、线路潮流、发电机电出力、购电量、弃风量天然气侧变量节点气压平方、管道流量、气源注入量、压缩机耗气量热力侧变量各热源供热量、热负荷节点热功率、电锅炉耗电量整数变量主要是机组启停状态以及分段线性化引入的分段选择变量。注意一点分段线性化变量如果做成了二进制或者整数模型就升级为混合整数规划求解规模会明显上升。能用连续的分段变量就尽量用连续的只在必要时引入0-1变量。3.3 约束矩阵组装与求解器调用Yalmip是Matlab下的最佳选择Matlab环境下做优化调度绕不开Yalmip。Yalmip是一个建模工具箱好处是语法简单不需要手动组装约束矩阵直接把变量、约束、目标函数丢进去它帮你转换格式并调用底层求解器。搭配Cplex或Gurobi使用稳定性远高于Matlab自带的linprog和intlinprog尤其对大模型而言。求解器选择看模型规模纯线性规划几百个变量Matlab自带的linprog能凑合带整数变量、上千个变量建议直接用Gurobi免费学术许可很容易申请混合整数规划模型Cplex和Gurobi都能处理二者性能差距不大看自己机器上装了哪个值得一提的是Yalmip在2023版之后对Gurobi的接口支持很完整解算信息里能直接读出gap、迭代次数和各边界。强烈建议把求解器参数设置的代码单独写成一个函数方便切换测试。4. Matlab代码框架实战拆解从数据到求解的一气呵成4.1 代码整体结构与main函数流程我的代码框架遵循数据-模型-求解-后处理四个模块分离的原则。main函数逻辑很清晰load_data.m加载电网、气网、热网所有参数build_variables.m定义所有决策变量build_constraints.m构建全部约束build_objective.m组装目标函数solve_model.m调用求解器求解plot_results.m输出结果到图表模块化最大的好处是改参数不用动模型结构换算例不用改代码。我见过不少人把所有内容堆在单个脚本里50行的数据定义加200行的约束构建改一个节点参数要找半天这是攀爬综合能源调度这座山时最不需要的麻烦。4.2 数据加载结构体组织参数的关键设计数据处理是整套代码最容易乱的地方。我推荐用结构体或表table保存数据而不是散落的裸变量。以电网数据为例我定义的结构体里包含基准功率、发电机参数表节点、最大最小出力、成本系数、爬坡率、负荷时序矩阵、线路参数表首端、末端、电抗、容量、节点数、线路数。气网结构体类似节点数、气源节点、气源出力上下限与成本、管道参数表、负荷时序、压缩机参数、节点气压限值。热网结构体包括热负荷时序、各热源类型与成本系数、供热上下限、电锅炉效率。这里有个关键细节所有功率、流量、压力的单位必须统一。电网常用MW气网常用kCF/h或m³/h热网常用MWth压力单位可能是kPa或MPa。混用单位是模型不可行最隐蔽的原因之一。我的做法是全部统一到标幺制电网和热网都转为MW基准标幺气网流量也折算成等效功率标幺保证约束矩阵中所有系数在数量级上可比这一步能省掉后期大量调参时间。4.3 变量定义和一键构建约束Yalmip写法的核心操作用Yalmip定义变量非常直接如下所示%% 电网变量 P_g sdpvar(n_gen, T, full); % 发电机出力 theta sdpvar(n_bus, T, full); % 节点相角 P_line sdpvar(n_line, T, full); % 线路潮流 P_buy sdpvar(1, T, full); % 上级电网购电功率 P_curtail sdpvar(n_wind, T, full); % 弃风功率 %% 气网变量 Pi sdpvar(n_gas_bus, T, full); % 节点气压平方 F_pipe sdpvar(n_pipe, T, full); % 管道流量 F_source sdpvar(n_source, T, full); % 气源注入量 %% 热网变量 H_chp sdpvar(n_chp, T, full); % CHP供热功率 H_gb sdpvar(n_gb, T, full); % 燃气锅炉供热功率 H_eb sdpvar(n_eb, T, full); % 电锅炉供热功率约束定义就是一条条往里加。这里需要注意Yalmip里加约束用方括号和直接拼接而且Yalmip支持矩阵不等式一起加。构建约束时我习惯用变量名做前缀比如Constraints [Constraints; P_g 0]; Constraints [Constraints; P_g P_g_max];4.4 关键约束的代码实现示例电网直流潮流和节点平衡这是最基础的%% 节点功率平衡直流潮流 for t 1:T Constraints [Constraints; sum(P_g(:,t)) P_buy(:,t) - P_curtail(:,t) ... sum(P2G(:,t)) P_load(:,t) sum(P_eb(:,t))]; end气网管道流量的分段线性化我单独封装成函数输入是管道流量变量输出是线性化后的一组约束。这里给一个简单的增量线性化结构示意。首先把管道两端的压力平方差向量P_diff和流量F_pipe建立映射然后对每个分段区间引入非负分段连续变量并加SOS2约束或二进制选择约束。SOS2约束在Yalmip里可以直接用sos2命令配合Gurobi/Cplex求解时效率很高。为了避免代码无限膨胀这里不贴完整实现逻辑就是标准的增量线性化三段式。耦合约束部分最典型的是CHP电热耦合可行域。以常见的热电联产机组为例电出力和热出力的关系可以近似为一组线性不等式。实现时把每台CHP机组在每个时段的电出力、热出力联系起来%% CHP机组的电-热耦合约束 for t 1:T % 热出力上限随电出力变化 Constraints [Constraints; H_chp(:,t) c1 * P_chp(:,t)]; Constraints [Constraints; H_chp(:,t) c2 - c3 * P_chp(:,t)]; endP2G耦合约束更简单用固定效率将耗电量映射为注入气网的天然气量。注意P2G的耗电在电网侧是负荷产气在气网侧是气源两个网络通过效率系数连接%% P2G设备耗电转气 G_p2g eta_p2g * P_p2g; % 实际建模中G_p2g直接进入气网该节点流量平衡4.5 求解与结果后处理别忽略求解状态检查求解调用代码很短但很多人在这翻车。optimize函数会返回求解状态信息正常输出是求解成功。我建议在代码里强制检查求解状态不要直接提取结果绘图。如果不可行立刻定位冲突约束如果求解器卡住先看gap多大。ops sdpsettings(solver,gurobi,verbose,2, ... gurobi.MIPGap,0.001, gurobi.TimeLimit,600); result optimize(Constraints, Objective, ops); if result.problem 0 % 提取结果绘图 else disp(result.info); error(求解失败请检查模型可行性); end5. 典型算例演示用一个小型区域综合能源系统跑通全流程5.1 算例系统参数设计我通常用这个规模来验证代码一个6节点电力系统、一个6节点天然气系统、一个6节点热力系统时间范围24个时段。这样的规模足够展示耦合调度的全部特征又不会让求解时间拖太长。详细介绍如下电力系统3台常规发电机组火电、1台风电场、与上级电网相连负荷峰值为180MW左右接入气网的燃气机组1台、CHP机组1台。 天然气系统2个气源点、6条管道、1台压缩机负荷为1台燃气锅炉供热和1台燃气机组发电负荷高峰集中在傍晚。 热力系统热负荷峰值约90MW供热来源包括1台CHP机组、1台燃气锅炉、1台电锅炉。电锅炉用来在风电大发的夜间时段消纳风电。5.2 独立调度与耦合调度的结果对比跑完模型后我会做一个对照实验来验证耦合调度的价值先分别独立优化电网和气网热网只管热平衡不参与联动再运行整套耦合模型。典型结果大致如下独立调度下夜间风电大发电力系统被迫弃风电锅炉因为没有纳入调度框架没有参与消纳。耦合调度下夜间电锅炉启动将风电转化为供热同时CHP降低出力燃气锅炉跟着降出力整体燃气消耗降低。总体运行成本下降约6%-9%弃风率从约15%降至3%以下。举个具体时段的例子凌晨2点风电出力高电力负荷低电价为低谷。独立调度只能让部分风电上网剩余弃掉。耦合调度启动电锅炉吸收多余风电供热同时减少燃气锅炉的燃气消耗和CHP的热出力热负荷依然满足但天然气的用量下来了。这就是耦合调度最直接的经济效益代码跑出来一眼就能看懂适合写进报告或者论文结论。5.3 结果图表的输出要点结果可视化我一般分三张图第一张电力系统24小时出力分布堆叠图标出发电机、风电、购电、电锅炉和P2G耗电。第二张气网气源出力和管道关键节点压力分布图检验压力是否在允许范围内。第三张热网热源构成堆叠图展示CHP、燃气锅炉、电锅炉的供热量分配。绘图代码用Matlab自带的figureplotstackedArea就够不需要额外的包。图表的导出统一用exportgraphics输出高清版式方便直接贴到论文或项目报告中。6. 踩坑实录把这套代码跑通过程中最值得说的五件事6.1 求解不可行不一定是约束错先查单位我在这套代码上最开始遇到不可行大概花了整整两天时间排查。单看每个约束都合理单看每个变量都有限度但模型就是报infeasible。后来逐个放大限值测试才发现是气压的单位不一致天然气源的节点气压我用了kPa而管道流量方程的分段线性化系数用的是MPa两个量纲相差1000倍直接把线性化区间的指数撑爆了。所以再强调一次所有参数在进入Yalmip前先做标幺化把每个变量的取值范围控制在接近几百万数量级内这能规避掉大量数值问题。6.2 分段线性化的段数和精度如何权衡管道流量方程分段线性化是这套代码里唯一需要人为把握精度的地方。分段数太少线性逼近误差过大气压约束形同虚设分段数太多变量数量成倍膨胀求解时间指数上升。我的经验是对每个管道按流量范围分成8到12段配合SOS2约束已经能在精度和计算速度之间取得很好的平衡。如果你对结果精度要求很高可以做一个收敛性检验脚本逐步增加分段数看目标函数和关键变量变化不足0.1%时取当前分段数这个步骤能有效说服审稿人或者导师。6.3 冷启动问题为什么加一个初始可行点能快十倍带整数变量的混合整数规划对初始可行点非常敏感。最开始我让求解器从零开始搜索一个包含约4000个变量和约3000个约束的模型要跑七八分钟。后来我给所有机组启停状态设置了一个合理的初始解比如全天开机Gurobi的MIP gap收敛速度明显加快求解时间压缩到一分钟以内。在Matlab里完成这个设置的方式是构造一个初始解结构体再在sdpsettings里指定gurobi.StartNumber或初始化变量对应的值。要注意初始解只是加速搜索不能随便给必须满足所有约束否则求解器会忽略它。6.4 求解器许可证与版本兼容别忽视MatlabYalmip这套环境最隐蔽的坑是版本兼容。Yalmip对Gurobi接口在R2019b之后的版本基本都是兼容的但如果你用的是破解版或旧版本求解器就很容易出现参数不生效、解法器named接口报错等诡异情况。我的建议是优先使用官方提供的Gurobi优化器并申请学术许可安装路径不要有中文系统PATH变量里加上Gurobi的bin目录。Cplex早期版本对Yalmip的接口在部分Windows系统上偶发dll加载失败直接换Gurobi往往能省下半小时折腾时间。6.5 模型简化与扩展的方向最后说一点关于代码可扩展性的思考。耦合调度模型往更精细的方向走有两条路一条是增加储能设备电储能、储气、储热这需要在约束里新增储能设备的soc状态方程另一条是引入气网动态管存效应即管道末端的天然气流量不再等于首端流量而是与管存变化相关这会把稳态模型升级为准动态模型求解难度进一步加大。但核心代码框架不用推翻只需要在气网模型部分替换约束即可。这也是我坚持模块化拆分的另一个原因当前工程实践中很多项目需要反复调整模型复杂度从稳态到准动态、从单时段到多时段都是同一个框架下的局部改动。这套代码从入手到完全跑通我大概迭代了三轮第一轮是搭线性模型看结果合理性第二轮加整数变量和分段线性化第三轮优化求解参数和数据结构。如果你的最终目标也是实现一套可复用的综合能源系统耦合优化调度代码建议也按这个节奏来先把最简单的线性版本跑出图确认每个网络的潮流和节点平衡都正确再逐步加复杂性。这样每一步的错误都有限定范围排查起来清晰得多。