1. 项目定位为什么值得做一套电气热综合能源SOCP程序在综合能源系统优化这个方向上用Matlab做电气热综合能源二阶锥优化模型已经是目前比较主流的工程做法。我这里要说的这套程序选的是39节点电网、6节点天然气网再加热网规模算不上大但足够覆盖电-气-热三网耦合时最核心的建模问题非线性怎么处理、多能流怎么平衡、求解器怎么收敛。整段代码用YALMIP搭二阶锥松弛交给MOSEK或Gurobi这类商业求解器求解最后得到一天24个时段各个发电机的出力、气源供气量、热源供热功率以及系统总运行成本。这套东西适合谁看一是刚接触综合能源优化、想把电热气耦合模型跑通的研究生二是做园区或区域综合能源规划需要快速搭建可行性评估模型的工程师三是本来会线性优化、想搞懂二阶锥到底怎么落地到代码里的朋友。它不是一个园区级别的精细化物理仿真程序而是一个面向日前调度和规划评估的优化模型追求的是在可接受的精度内让大规模多能流问题能被全局求解器稳定解出来。1.1 三网耦合到底在耦合什么先明确一下模型范围。电网是简化的39节点系统参考IEEE New England 39节点拓扑包含发电机、负荷、输电线路重点考虑网损、节点电压和支路潮流。气网是6节点系统包含气源、管网和负荷重点考虑天然气流量、节点压力以及加压站对压力水平的支撑。热网我没有做完整的水力-热力动态模型而是按“热功率网络”来处理包含热源、热负荷、供回水温度和循环流量适合日前尺度的热力平衡分析。三网之间的耦合关系是整套模型的灵魂。最常见的耦合设备有三类燃气轮机一头消耗天然气一头发电属于电网和气网的耦合热电联产机组同时输出电和热是电、气、热三个系统共同交汇的地方电锅炉或热泵消耗电能产生热能把电网和热网联系在一起。子系统规模核心状态变量关键约束电网39节点发电机出力、支路潮流、节点电压平方潮流平衡、电压上下限、支路容量气网6节点气源流量、管存流量、节点压力平方气源容量、节点压力、管道流量约束热网热源/热负荷节点热源供热量、管道热功率热功率平衡、热源出力上下限把这三个网络放到同一个优化问题里约束矩阵会出现明显的块结构块与块之间靠耦合设备连接。这种结构越早用YALMIP或同类建模工具做好变量分区后面排查问题就越轻松。1.2 为什么偏偏选SOCP而不是LP或NLP很多第一次做综合能源优化的同学会问既然有线性规划也有通用的非线性求解器为什么还要绕一圈做二阶锥我实际对比过答案很直接线性模型算得快但把电网潮流、气网流量这些关键非线性全部拍平了算出来的调度方案往往偏乐观放到真实物理系统里根本不满足节点电压或管道压力约束。通用非线性模型比如直接在Matlab里用fmincon求解潮流方程理论上是精确的但因为是非凸问题初值稍微给得不对就会掉进局部最优解而且气网压力方程和热网温度方程连在一起之后收敛性会变得非常差。二阶锥规划的好处在于它把非凸的等式约束先松弛成凸的二阶锥不等式再用现代内点法求解得到的是全局最优解而且对大规模模型的求解效率远高于传统非线性求解器。代价是结果满足的是松弛后的不等式不一定是原始等式所以需要在求解完成之后做一个“紧性检验”。这个环节我会在后面的章节专门展开。用个生活化的类比一条不规则的河道你想知道船能不能通过直接精确测量每一处弯曲很费劲二阶锥的做法是把河道外面加一根笔直的护栏让所有可能的航线都落在护栏内部然后在这个更大的安全区域里找最省油的路。大部分情况下最优航线刚好贴在护栏上那就是真实航线如果贴不住说明护栏本身加宽得太多需要调整。2. 建模细节电网、气网、热网三个模块的取舍2.1 电网39节点系统的DistFlow与二阶锥松弛电网部分我没有用完整的交流潮流而是采用分支潮流模型简称DistFlow。这个模型在配电网和输电网的凸优化里非常常见。对每条支路从节点 i 到节点 j定义有功潮流 P_ij、无功潮流 Q_ij、线路电流平方 L_ij 和节点电压平方 V_i。约束可以写成这样节点有功平衡注入的有功等于流出支路有功加上线路电阻损耗。节点无功平衡同理需要计入线路电抗损耗。电压降落关系V_j V_i - 2(r_ij * P_ij x_ij * Q_ij) (r_ij^2 x_ij^2) * L_ij。这里面有一个非凸等式P_ij^2 Q_ij^2 V_i * L_ij。直接丢给Gurobi它不认这种等式因为可行域不是一个凸区域。二阶锥松弛的思路就是把这个等式放宽成 V_i * L_ij P_ij^2 Q_ij^2然后把 V_i 和 L_ij 的关系改写成标准二阶锥形式for k 1:n_branch for t 1:T i bus_fr(k); j bus_to(k); Constraints [Constraints, ... cone([2*Pij(k,t); 2*Qij(k,t); Vi(i,t)-Lij(k,t)], ... Vi(i,t)Lij(k,t))]; end end这个写法在YALMIP里很直接。cone的第一个参数是向量第二个参数是标量含义是第一个参数的二范数不超过第二个参数。有两点需要特别提醒第一Lij是电流平方变量不要习惯性写成电流第二松弛后要检查等式两侧的误差如果Lij * Vi明显大于P^2Q^2说明松弛不紧电压幅值算出来会失真。一般情况下目标函数里只要包含网损项优化过程会自动把电流压到尽量低松弛自然就是紧的。2.2 气网6节点系统与管道流量方程处理气网我选了6个节点规模不大但足够展示天然气网建模的典型套路。天然气系统里节点压力是一个重要状态量管道内流量的平方和管道两端压力平方差近似成正比。如果直接用压力作为变量会出现开根号和非线性乘积很难处理。工程上通行的办法是把“节点压力平方”当作一个新变量例如用 pi_i 表示压力平方。这样管道流量约束可以写成 f_km^2 C_km^2 * (pi_k - pi_m)其中 C_km 是管道参数f_km 是从节点 k 流向节点 m 的流量。这里有个方向性问题。如果管道流量方向已知比如从气源到负荷那么直接保证 pi_k pi_m上面这个不等式就是凸的。如果方向未知需要引入0/1变量判断流向问题就从SOCP变成了MISOCP求解难度立刻上一个台阶。我在程序里采取的是折中方案先根据气源位置和负荷位置固定每条管道的流向实测下来求解速度很快结果也符合物理直觉。只有在阀门或双向管道较多时才需要上整数变量。MATLAB里对应的约束片段可以写成Gpipe sdpvar(n_pipe, T); % 管道流量 Pi sdpvar(n_gas_node, T); % 节点压力平方 Constraints [Constraints, Pi 0]; for k 1:n_pipe from gas_fr(k); to gas_to(k); Constraints [Constraints, ... Gpipe(k,:).^2 Cpipe(k)^2 * (Pi(from,:) - Pi(to,:))]; end注意这是一个二次约束不是标准的二阶锥锥体约束但YALMIP会把它识别成二次锥MOSEK和Gurobi都能处理。使用Gpipe.^2这种写法在Matlab里是逐元素平方配合YALMIP生成的是一个二次锥约束的矩阵形式。2.3 热网热功率平衡与温度水力的简化热网部分我的处理比电网和气网都要粗。完整的热网动态模型需要管道的传输延迟、散热损失、供回水温度变化、水泵压力等多个维度这些要全部做进优化模型里会让二阶锥问题瞬间变成复杂的混合整数非线性问题。所以程序里采用“热功率网络”的抽象每个热负荷节点需要一定热功率热源节点通过供热管道把热功率送到负荷节点管道热损按固定比例近似。热功率本身由流量、比热容和供回水温差决定如果保持供水温度和回水温度不变热功率就和流量成正比。另一种更常见的近似是直接把热功率当作“广义功率流”类似电网有功功率这样热网就是一个线性网络只需要满足热源出力上下限、热负荷需求和管道容量约束。这样的处理在日前调度层面已经足够准因为温度动态的时间常数通常是几十分钟到几小时日前优化只需要保证总热量的供需平衡和热源调节范围。在程序中热网的约束大概是这样每个节点的热平衡流入热功率加上热源出力等于热负荷。热源出力上下限热电联产机组的热出力和电出力还要满足联合可行域。管道容量约束每条管道传递的热功率不超过上限。这样处理的好处显而易见模型线性没有额外非线性项不会破坏整个SOCP问题的凸性。坏处是无法评估热网末端用户的实际供水温度如果后续需要做台区或楼宇级精确调控需要在这个基础上扩展热力学动态约束。3. Matlab程序整体结构与求解流程3.1 主程序框架和数据结构整套程序我按模块拆开写避免在一个脚本里堆两千行代码。主程序负责读取数据、连接各子网模型、调用求解器、输出结果。数据结构上用结构体保存电网、气网、热网参数这样后续扩展时不需要频繁改动主脚本。建议的主程序流程是这样加载系统参数电网线路表、气网管道表、热网负荷表。定义优化变量发电机有功无功、支路潮流、电压平方、气源流量、节点压力平方、热源功率。生成各子系统的约束电网平衡、气网平衡、热网平衡。加入耦合设备约束燃气轮机、热电联产、电锅炉。构建目标函数发电成本、购气成本、供热成本之和。调用optimize求解。结果后处理提取变量值、检验约束紧性、绘制调度曲线。这种模块化写法还有一个好处调试的时候能单独把气网模块拿出来固定电网或热网的变量先看气网能不能收敛。我一直坚持一个原则多能流系统程序最大的风险不是算法复杂而是子系统模型本身的错误被耦合关系掩盖住。模块化就是用来对抗这种风险的。3.2 YALMIP建模关键代码下面这段是程序里比较核心的建模片段只保留关键部分。T 24; n_gen 10; n_branch 46; n_gas_source 2; n_heat_source 3; Pg sdpvar(n_gen, T); % 发电机有功 Qg sdpvar(n_gen, T); % 发电机无功 Pij sdpvar(n_branch, T); % 支路有功 Qij sdpvar(n_branch, T); % 支路无功 Vi sdpvar(39, T); % 节点电压平方 Lij sdpvar(n_branch, T); % 支路电流平方 Fg sdpvar(n_gas_source, T); % 气源供气量 Gpipe sdpvar(n_pipe, T); % 气网管道流量 Pi sdpvar(6, T); % 气网节点压力平方 Hh sdpvar(n_heat_source, T); % 热源热功率 Constraints []; % 电网二阶锥约束前面已给出 % 气网管道流量约束前面已给出 % 热电联产耦合约束 Hh(1,:) 0.8 * Pg(1,:); Hh(1,:) 1.2 * Pg(1,:); % 目标函数 Cost sum(sum(c_gen .* Pg)) sum(sum(c_gas .* Fg)) sum(sum(c_heat .* Hh)); ops sdpsettings(solver, mosek, verbose, 2); optimize(Constraints, Cost, ops);注意代码里的c_gen、c_gas、c_heat是单位成本矩阵维度要和对应的变量一致。很多同学在YALMIP里最容易报错的原因就是维度对不上矩阵加法都会直接崩掉。建议每构造一个约束就用size或者length检查一下变量的维度。3.3 求解器配置与数据尺度问题二阶锥优化模型对求解器有要求不是所有求解器都能直接处理。YALMIP自带的sedumi和sdp3虽然能用但速度慢且容易内存爆炸。我测试下来MOSEK 的SOCP引擎最稳定Gurobi在加入整数变量之后更有优势。如果你只有学术版授权Gurobi和MOSEK通常都能申请Cplex也可以但界面和授权相对老一些。求解器设置上有几个参数值得关注。对于MOSEKmosek.tol_rel_gap可以控制最优间隙默认值通常就够用不需要调太紧。对于Gurobi如果模型是纯连续SOCP我会把gurobi.Threads设成实际物理核数如果加了机组启停0/1变量变成MISOCP则把gurobi.MIPGap设为1e-4或更松。目标函数量级相差太大时比如总费用几百万但电压约束数量级只有1.0极易出现收敛判定问题。解决办法很粗暴对变量统一标幺或归一化。功率基准取系统容量100MW压力基准取10bar温度基准取100摄氏度把这些量纲调到1附近求解器能省掉大量无效迭代。4. 我踩过的坑收敛、精度与求解时间4.1 常见现象对照表整理一下我在调试这套39节点电热气模型时遇到的典型问题以及对应的解决思路。现象可能原因处理方式求解报infeasible子系统约束方向写反或耦合设备可行域过小分模块单测先验证电网、气网、热网各自可行SOCP收敛但结果明显偏离物理二阶锥松弛不紧电流平方变量虚高增加网损项到目标函数检查紧性误差求解时间突然暴涨加入了整数变量或目标函数尺度差太大用连续变量试算归一化数据检查是否误用整数变量气网流量全为负管道压力方向设置错误固定气源到负荷方向或引入大M判断方向热网热功率不平衡热负荷节点对应关系写错画出节点连接表逐点核对热源和负荷编号YALMIP报错“Unable to assign”变量维度不对或矩阵拼接错误用length和size逐行检查不要一次性拼接大矩阵这个表我在实际调模型的时候基本贴在旁边每遇到一个症状就对照一遍比重新翻公式高效得多。4.2 排查顺序与初值技巧遇到模型不收敛我的调试顺序固定如下。先断开所有耦合设备也就是把燃气轮机、热电联产、电锅炉各自拆开分别让电网独立优化、气网独立优化、热网独立优化。这一步能暴露出最基础的建模错误。比如电网不收敛多半是DistFlow里电压更新公式的正负号有问题气网不收敛多半是压力平方差方向不对。子系统都能跑通后再把热电联产和燃气轮机一个个接回去。接一个解一次不要一口气全接上。我印象最深的一次电网和气网分别都很好但接上燃气轮机后气网管道流量全部越限查了很久才发现燃气轮机耗气量公式里的效率系数写反了。单独看每部分都不显眼合在一起就变成了全局不可行。如果模型本身可行只是求解器告诉你不收敛还可以做“热启动”。YALMIP不直接支持热启动所有求解器但MOSEK支持assign机制。先把线性化版本模型解一遍用结果作为SOCP变量初始值可以明显减少迭代次数。还有一种土办法把目标函数里的网损项或运行成本项先放大10倍让优化方向更明确等找到可行解后再恢复原始系数。这个方法虽然看似粗暴实测对难收敛模型很管用。4.3 松弛紧性检验要不要做要而且一定要做。二阶锥优化模型天然是松弛模型求出的未必是原始非凸问题的全局最优解。如果最优解处所有被松弛的锥约束都取等号那么松弛是紧的结果等效于原始问题。如果不紧最优解可能在物理上不存在调度方案拿过去根本没法执行。我的检验方法很简单。拿电网的支路约束来看计算每个时段每条支路的松弛间隙gap Lij_value .* Vi_from_value - (Pij_value.^2 Qij_value.^2); max_gap max(gap(:));如果max_gap相对于支路潮流平方的量级小于1e-4可以认为松弛足够紧。如果偏大说明目标函数里缺少必要的“压紧”机制。最简单的调整办法是把网损项放进目标函数或者对电流平方变量加一个小惩罚系数例如0.001 * sum(Lij(:))。多加一个线性惩罚项不会破坏凸性却能把松弛误差压下来。气网的管道约束也类似要检验Gpipe^2是否贴近C^2*(Pi_from - Pi_to)不过气网允许少量松弛因为压力本身有调节裕度。4.4 求解时间优化心得从纯连续SOCP升级到含机组启停的MISOCP求解时间通常不是线性增加而是爆炸式增加。我这套39节点模型纯连续版本大概几十秒到一两分钟就能解完一旦加入10台机组一天的启停状态240个0/1变量直接变成几十分钟甚至跑不完。解决办法有三个方向。第一把一天24个时段聚合成几个典型时段先做一个粗粒度模型验证方案。第二对机组启停变量启用对称破缺约束比如相同机组按成本排序强制递增或递减能有效减少分支定界树的搜索空间。第三把气网和热网约束在初始版本中放松成线性化约束先求得整数解再固定整数变量回头精算SOCP连续解。这个“整数-连续”两步法在实际工程里非常好用。5. 程序还能怎么改扩展方向与个人心得5.1 从单目标到碳交易和需求响应这套程序的目标函数目前是发电成本、购气成本和供热成本之和。如果需要做碳电协同只需要增加碳排放约束或者碳成本项。天然气的碳排放系数可以直接乘到气源供气量上电网侧则根据各机组燃料类型确定排放系数。加了碳价之后优化结果会出现一个很明显的现象低效煤电机组的出力下降热电联产机组和气电的出力上升热网侧会更倾向于利用余热。这个逻辑对于做园区碳达峰方案特别有价值。需求响应也可以很容易加进来。电网侧增加可削减负荷变量气网侧增加可中断气负荷热网侧增加可调热负荷。每个可调负荷都有调节成本目标函数里增加对应项约束里加入总削减量上限和单时段削减量上限就能模拟综合灵活性的调度空间。二阶锥框架完全兼容这样的线性扩展不需要改变求解器类型。5.2 从日前优化到滚动调度我后来把这套程序从“一次性求解24小时”改成了滚动调度也就是模型预测控制思路。每个控制周期前用负荷预测更新未来4小时或8小时数据只执行当前时段的决策下一个周期重新优化。这样改有个好处能规避气网和热网动态时间常数带来的预测偏差。但要注意每次求解都需要保留上一周期的状态变量作为当前周期约束比如热网的热源出力不能瞬间跳变需要增加爬坡约束。这在实际中比单纯拉长优化周期更贴近运行需求。5.3 一点避坑心得如果你准备复现类似程序我的建议是千万不要试图一次性把39节点电网、6节点气网和复杂热网全部搭完再求解。先把39节点电网跑通一天96时段再去接气网和热网。我一开始三张网同时上出了问题都不知道是潮流越限还是气压掉压。先把每个网的模型单独验证好再耦合调起来会快得多。另外YALMIP的报错信息有时候比较抽象尤其是变量维度不匹配时提示信息不会直接告诉你哪一行出错。我的土办法是每构造一组约束后立刻运行一次check(Constraints)让YALMIP报告最大约束违例量逐步缩小问题范围。做综合能源优化一半时间在建模一半时间在排错把排错流程理顺了开发效率至少提升一倍。