1. 碳中和目标下,为什么“电-气互联+无功优化”成了硬需求
先说我调试这套模型时遇到的一件印象深刻的事:把天然气网络的约束加进无功优化模型之后,最优解里燃气轮机的出力结构发生了显著变化,系统网损和电压质量同时改善,而单独做无功优化时根本看不到这种联动效果。这件事让我意识到,在“碳中和”背景下,电气互联系统的有功-无功协同优化已经从一个学术概念,变成了实际工程里绕不开的问题。
先说清楚这套东西是什么:电气互联系统,指的是电力网络和天然气网络通过燃气轮机、电转气(P2G)等耦合元件连接起来的综合能源系统。碳中和目标下,新能源装机占比快速提升,风电、光伏大规模接入,但这部分电源本身没有传统同步机那样的无功支撑能力。有功出力受气象条件影响剧烈波动,随之而来的是节点电压频繁越限、无功功率远距离流动、网损增加等问题。传统做法是把无功优化单独拿出来做,通过调节发电机无功出力、变压器分接头、无功补偿装置来改善电压分布——但在新能源高渗透和电-气深度耦合的背景下,这种“只看电网”的做法会错过很多全局最优解。
电气互联系统让我感兴趣的另一个点,在于它把“看不见”的气网灵活性引入了电网调度。燃气轮机的电出力和气网供气量直接耦合,P2G装置又能把富余电力转化为天然气。这意味着电网的有功调度不再只受发电机组本身约束,还受气源供应能力、管道流量、节点气压的约束;反过来,气网的管存和调度弹性也能帮电网的调峰调压“解围”。有功调度的可行域变了,无功优化的最优解自然也变了。这就是“协同优化”四个字的真正含义。
这篇文章适合三类人读:一是做电力系统/综合能源系统优化方向的研究生,需要一套可以直接跑通、有完整建模思路的Matlab参考实现;二是做电网规划或运行控制的工程师,想了解电-气耦合对电压无功管理到底有什么实际影响;三是刚入门、想搞懂“什么是电气互联”“什么是有功-无功协同”的初学者。我会把模型怎么搭、约束怎么写、代码怎么组织、数值结果怎么看、调试中踩过的坑都讲一遍,尽量让读者能照着复现。
2. 协同优化模型怎么搭:电力网络、天然气网络与耦合元件建模
2.1 电网侧:为什么这里必须用AC潮流,而不能用直流潮流
和无功、电压相关的优化问题,第一道坎就在潮流模型选择上。很多做有功调度的人习惯用直流潮流(DC Power Flow),因为它线性、计算快、好求解。但直流潮流只解P不解Q,连电压幅值都没有,根本没法做无功优化。凡是涉及电压越限、无功补偿、变压器分接头调节的问题,至少要用AC潮流模型。
我这里采用极坐标形式的牛顿-拉夫逊潮流方程。对电网节点 i,有功和无功注入满足:
- P_i = V_i Σ V_j (G_ij cosθ_ij + B_ij sinθ_ij)
- Q_i = V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij)
其中 G_ij、B_ij 是节点导纳矩阵的实部和虚部,θ_ij = θ_i - θ_j 是相角差。优化变量里会包含各节点电压幅值 V_i 和相角 θ_i,所以这个方程组既是等式约束,也是决定网络状态的物理核心。
在Matlab里落地时,我建议先写一个独立的潮流模块,用牛顿-拉夫逊迭代求解给定注入下的系统状态,再用这个模块做初值生成和结果校验。不要一上来就把潮流方程直接塞进优化求解器——非凸、强非线性,很容易初值敏感导致不收敛。先把潮流算稳,再谈优化。
2.2 气网侧:Weymouth方程与节点气压约束
天然气网和电网在数学结构上有相似之处,但细节差异很大。气网节点有气压(压强)状态变量,管道有流量,节点上有气源、气负荷和耦合元件。描述稳态气网管道流量的经典方程是Weymouth方程:
- f_km = C_km · sgn(π_k - π_m) · sqrt(|π_k² - π_m²|)
其中 f_km 是管道流量,π_k、π_m 是管道两端节点气压,C_km 是管道综合系数(与管径、长度、摩阻系数、气体性质相关)。
这个方程有两个要特别注意的地方。第一,它带符号函数和根号,是非光滑非线性约束,直接用内点法求导时在零流量附近容易出问题。第二,气压变量是平方形式,很多文献会定义辅助变量 Π = π²,把方程改写成流量平方等于压强平方差的形式,便于做二阶锥松弛或分段线性化。
此外还有节点流量平衡约束:每个气网节点的注入气流量(气源产气 + P2G产气)减去节点气负荷(常规气负荷 + 燃气轮机耗气),等于流出该节点的管道流量总和。气源有产气上下限,节点气压有上下限,这两类约束在实际调试中最容易把模型搞“不可行”,后面我会专门讲。
2.3 耦合元件:燃气轮机与P2G的数学表达
电气互联系统的核心是耦合元件,它们决定了两个网络之间如何交换能量。我用两类耦合元件覆盖最常见的场景。
燃气轮机把天然气转化为电力,是“气转电”方向。它的模型本质是一个能量转换环节,常用线性或分段线性函数描述耗气量与电出力的关系:
- F_GT = α · P_GT + β
F_GT 是耗气量(单位时间内消耗的天然气体积/能量),P_GT 是燃气轮机的电出力,α 是耗气系数,β 是空载耗气量。这个线性近似对中低负载率区间精度尚可,如果要做更精细的建模,可以改成凸分段线性函数。燃气轮机接入电网的节点上,它的有功注入会直接进入该节点的潮流方程;同时它消耗的气量会进入气网的节点流量平衡方程。
P2G是“电转气”方向。电解水制氢再甲烷化,把富余电力转化为天然气注入气网:
- F_P2G = η_P2G · P_P2G / HHV
P_P2G 是P2G消耗的电功率,η_P2G 是综合转换效率,HHV 是天然气高位热值。P2G消耗的是电网有功功率,既是电网负荷又是气源,双向耦合关系一目了然。在实际算例中,P2G的容量通常远小于燃气轮机,因为当前P2G效率和经济性还不足以大规模应用,但在碳中和场景下它的战略意义在于“消纳弃风弃光、平抑新能源波动”,所以模型里还是要留这个接口。
2.4 无功调节资源建模:发电机PQ能力曲线、变压器分接头与无功补偿装置
做了这么多准备工作,终于回到主角——无功优化。无功调节手段在模型里通常分三类:
发电机(含燃气轮机、常规火电):无功出力范围不是一个矩形,而是由PQ能力曲线决定的。通俗理解,发电机有功出力越高,可用的无功调节范围就越窄;在靠近额定有功时,无功上下限被定子绕组发热、转子励磁限制等物理因素压得很厉害。这个约束在协同优化里特别重要,因为燃气轮机有功出力随气网约束变化,其无功可调范围也在动态变化——这就是“有功影响无功”的物理机制之一。
有载调压变压器(OLTC):通过改变变比 k 来调整电压。变比通常是一组离散值(比如 0.95~1.05,步长0.0125),这会把模型变成混合整数问题。如果只做连续松弛,求解很快,但可能得到现实中调不到的分接头位置,所以严谨的做法是采用混合整数规划求解,或对连续解做离散回整+回代验证。
无功补偿装置:并联电容器/电抗器通常是离散投切(分组),静止无功发生器(SVG/SVC)可以连续调节。算例里我会按离散分组建模,兼顾工程实际和求解精度。
这三类装置的无功上下限、变比范围、投切组数,都会以不等式约束的形式进入优化模型。协同优化的难点和亮点,恰恰在于这些无功手段的“上游”——燃气轮机有功出力,现在反过来受气网供气能力影响。这形成了一个完整的从天然气源到电力负荷、从有功出力到无功调节的联动链条。
3. 目标函数与约束体系:协同是如何在数学上发生的
3.1 目标函数选择:从网损最小到综合运行成本
目标函数是整个优化模型的“指挥棒”。碳中和背景下的目标函数设计,我个人不建议只做单目标——即使只跑单目标,也要想清楚它在代表什么。我在这套模型里用的是带排放项的综合运行成本最小化:
- min C_total = Σ(火电/燃气机组发电成本)+ Σ(气源购气成本)+ Σ(碳排放成本)+ 网损/惩罚项
发电成本用二次函数:C_PG = a_i · P_i² + b_i · P_i + c_i。购气成本是气源产气量的线性函数,用分时气价可以体现气网侧的调峰信号。碳排放成本按机组碳排放强度乘碳价计算,燃气轮机单位碳排放低于燃煤机组,天然气的清洁属性会在目标函数里自动形成对气网侧的牵引力。
如果重点要观察电压和网损,我建议单跑一个“网损最小化”场景作为对照。但我必须提醒一点:网损最小和运行成本最小时的最优解往往不是同一个解,有些情况下为了降网损会提高某些机组的出力,总成本反而上升。在实际项目里,目标函数的选择要跟场景诉求一致,这也是我给做研究的朋友的忠告。
3.2 等式与不等式约束的完整表达
模型的约束体系分成三块:电网约束、气网约束、耦合约束。
电网侧的等式约束就是2.1节那组潮流方程。不等式约束包括:
- 节点电压上下限:V_min ≤ V_i ≤ V_max,通常取0.95~1.05 pu
- 发电机有功/无功出力上下限
- 线路传输容量限制(视算例规模可选)
- 变压器变比范围和无功补偿装置容量限制
气网侧的等式约束是节点流量平衡:气源注入 + P2G注入 - 燃气轮机耗气 - 常规气负荷 = 管道流出。不等式约束包括气源产气上下限、节点气压上下限(通常用辅助变量 Π 表示)、管道流量上限(用Weymouth方程自然限制)。
耦合约束把两套网络连成一体:
- 燃气轮机耗气量 F_GT 作为气网负荷出现在气网节点平衡方程中,同时其电出力 P_GT 出现在电网节点注入中
- P2G耗电量 P_P2G 作为电网负荷,产气量 F_P2G 作为气源注入
在Matlab里实现时,我习惯把电网节点、气网节点、耦合元件的关系先画一张连接表写在注释里,再按连接表逐条写约束。这样做的好处是:一旦约束写错导致不可行,排查时能按图索骥,而不是在一堆矩阵里瞎找。
3.3 为什么“有功-无功协同”不是简单地把两个问题拼在一起
很多初学者会把协同优化理解为“把有功调度问题和无功优化问题放到一个模型里一起解”,这是对但又不完全对。区分的关键在于,两个问题叠加在一起并不等于协同。
传统无功优化只优化无功调节手段,有功出力按经济调度结果固定不动。传统有功调度则反过来,只盯有功和机组组合,电压用潮流校验兜底。二者拼接后,确实同时包含两类变量,但没有建立“耦合机制”——即气网约束如何改变有功可行域、有功可行域如何影响无功调节空间。
真正的协同优化要做到两层联动:
第一层,网络间联动:燃气轮机的有功出力下限/上限受到气网节点气压约束。比如某节点气压偏低,该节点燃气轮机的最大出力可能受限,这部分有功缺口只能由其他机组或P2G消纳方式补上。此时电网中各发电机出力结构变化,全网无功分布和电压水平随之改变,最优的无功补偿策略也变了。
第二层,源-网-荷联动:P2G消耗有功以产气,会加大电网负荷,但同时增加气网供气量,可能解除燃气轮机出力受限约束。这种“有功换气、气再换有功”的循环关系,在单网模型里根本表达不出来。
我在具体算例里做过一个对比:固定气网约束,先独立做无功优化,再做协同优化,两种场景下最优解里同一个燃气轮机的无功出力差距最高达到8个百分点。这说明协同不是“加个约束再算一遍”,而是从可行域塑造阶段就改变了最优解的搜索空间。
3.4 求解思路:内点法、SOCP松弛与混合整数规划的取舍
这个模型本质上是非凸非线性规划(NLP),再加上变压器分接头和电容器投切的整数变量后,又变成非凸混合整数非线性规划(MINLP)。直接用通用求解器硬啃,规模稍大就跑不动,或者陷入局部最优。
实践中我总结出三条路线:
路线一:连续化+内点法。把整数变量松弛成连续变量,用Matlab的fmincon(内点法)求解。优点是实现简单、调参少;缺点是无功补偿的“分组投切”特性被抹掉,松弛解需要离散化回代校验。适合快速验证模型逻辑是否正确、算例规模较小(几十个节点以内)的场景。
路线二:MISOCP规划。对输电网潮流做二阶锥松弛(SOCR),把非凸潮流方程转化为凸二阶锥约束,天然气Weymouth方程也通过引入Π变量做同样的二阶锥松弛。整数变量保留,最终用YALMIP+Gurobi/CPLEX求解MISOCP。这条路能拿到质量很高的解,甚至在很多凸松弛条件下是全局最优。缺点是松弛后要检验解的“紧性”——即松弛前后目标函数值是否可忽略地被压缩,否则松弛解可能物理上不成立。
路线三:启发式+非线性优化交替迭代。外层用遗传算法或粒子群搜索离散变量(分接头位置、电容器组数),内层用fmincon解连续变量子问题。优点是灵活、能处理大规模整数变量;缺点是没有最优性保证,而且求解时间随变量维度增长很快。
这套模型里我最终用的是路线二为主、路线一为辅:先用fmincon跑通全流程做逻辑验证,再切换YALMIP+MISOCP做精细求解。两种结果对不上时,百分之八九十是建模约束写错了,而不是求解器的问题——这条经验帮我省过很多次Debug时间。
4. Matlab代码实现:数据结构、核心算法与关键代码片段
4.1 算例系统怎么选:电网用IEEE改造,气网用自定义节点
代码复现的第一步是选算例。电网部分我建议直接用Matpower自带的IEEE 14节点系统(case14),它规模适中,节点电压等级、变压器、无功补偿都有现成数据,文档也比较全。在case14基础上,我做了三处改造:
- 把部分火电机组替换为燃气轮机,并在对应节点上建立与气网节点的耦合关系
- 在某个负荷节点上加一组P2G装置(容量取3~5 MW,效率0.6~0.7)
- 在接入新能源的节点上修改电压上下限约束,模拟弱支撑节点场景
天然气网络部分没有Matpower这种成熟工具,需要自己搭建。我建了一个简化的6节点气网:2个气源节点、1个压缩机节点(本模型先忽略压缩机非线性,作为普通节点处理)、3个负荷节点。气压基准取6.0 MPa,管道系数C_km按管长300 km、直径0.6 m、摩阻0.01量级估算。节点气压上下限取基准值的±15%。
电网和气网的对应关系是:燃气轮机从气网节点1、3取气,P2G注入气网节点2,常规气负荷挂在节点4、5、6。整体规模不大,但完全覆盖了耦合元件的各类交互方式,跑一个MISOCP大概需要10到60秒,足够用来做机理研究。
4.2 程序架构与数据结构设计
Matlab代码的架构我按“数据—潮流—优化—后处理”四层组织:
main.m % 主流程:加载数据-构建模型-求解-结果分析 load_case.m % 电网数据(Matpower格式)和气网数据 build_ybus.m % 构建电网节点导纳矩阵 power_flow.m % 牛顿-拉夫逊潮流求解器 build_opf_model.m % 构建协同优化模型(YALMIP) solve_optimization.m % 求解调度 plot_results.m % 可视化结果数据结构方面,电网用Matpower的mpc结构体,气网用自定义struct,节点和管道数据独立存储,耦合元件单独一个connections表。我这里贴一段气网数据结构定义,方便参考:
% 气网数据格式 gasnet = struct(); gasnet.nodeNum = 6; % 气网节点数 gasnet.nodeID = [1;2;3;4;5;6]; gasnet.basePress = 6.0; % 基准气压 MPa gasnet.pressMin = 5.1 * ones(6,1); % 节点气压下限 gasnet.pressMax = 6.9 * ones(6,1); % 节点气压上限 gasnet.sourceNode = [1;2]; % 气源所在节点 gasnet.sourceMin = [0;0]; gasnet.sourceMax = [160;140]; % 气源产气上限 gasnet.pipe = [1 3; 2 3; 3 4; 4 5; 4 6]; % 管道连接 [起点 终点] gasnet.C = [8.5; 7.8; 6.9; 5.2; 4.6]; % 管道综合系数 gasnet.loadNode = [4;5;6]; % 常规气负荷节点 gasnet.loadVal = [30; 25; 20]; % 常规气负荷标幺值注意一点,气网数据内部最好统一用标幺值,只在输入输出接口处做有名值换算。否则管道系数、气压、流量单位混在一起,很容易差出好几个数量级——我是吃过这个亏的,后面专门展开说。
4.3 牛顿-拉夫逊潮流与初值设置的关键细节
潮流求解是优化模型的底座。牛顿-拉夫逊迭代的核心是雅可比矩阵修正,代码框架如下:
function [V, theta, iter] = power_flow(mpc) % 简化版牛顿-拉夫逊潮流 Ybus = build_ybus(mpc); n = size(Ybus, 1); V = ones(n, 1); % 电压幅值初值 theta = zeros(n, 1); % 相角初值 tol = 1e-8; maxIter = 30; for iter = 1:maxIter P_calc = real(V .* conj(Ybus * V)); Q_calc = imag(V .* conj(Ybus * V)); dP = mpc.Pspec - P_calc; dQ = mpc.Qspec - Q_calc; % 对PQ节点修正dP和dQ,对PV节点只修正dP dP(mpc.type ~= 1) = 0; dQ(mpc.type ~= 1) = 0; if norm([dP; dQ], inf) < tol break; end J = jacobian_matrix(Ybus, V, theta, mpc); dx = J \ [dP; dQ]; theta = theta + dx(1:n); V = V .* exp(1i * dx(n+1:2*n)); % 幅度修正用导数投影 end end(此代码省略了雅可比矩阵的完整推导,真实实现建议参考Matpower的Newton方法源码——那里有经过充分工程验证的数值处理。)
初值设置上最大的坑在无功优化场景:风机节点和重负荷节点的电压初值不能都取1.0 pu。我常用的办法是先用普通潮流估算一遍大致的电压水平,把初值设成潮流电压的0.95~1.05倍,再交给优化器迭代。这样内点法收敛速度快很多,也少了很多“矩阵奇异”的报错。
另一个值得注意的细节是求和顺序对数值精度的影响。用Matlab写潮流雅可比矩阵时,如果节点数超过100,直接三重循环会慢到让人怀疑人生。我建议用稀疏矩阵和向量化运算:
Ybus = sparse(Ybus); % 显式转为稀疏 % 或对雅可比矩阵用稀疏模式预分配4.4 优化求解器接入:fmincon、YALMIP与MISOCP建模
在跑通潮流的基础上,优化模型就可以搭建了。先做一个fmincon版本验证逻辑,再上YALMIP+MISOCP做精细求解。fmincon版本代码较长不全部贴出,核心思路是定义非线性约束函数,把潮流方程和气网约束以c(x) ≤ 0、ceq(x) = 0的形式写进去。
YALMIP是我更推荐的实现层。它能把上述模型用接近数学公式的语法表达,然后自动选择求解器。关键代码片段如下:
% 决策变量 V = sdpvar(nGrid, 1); % 电网节点电压幅值 theta = sdpvar(nGrid, 1); % 电压相角 Pg = sdpvar(nGen, 1); % 发电机有功出力 Qg = sdpvar(nGen, 1); % 发电机无功出力 press = sdpvar(nGas, 1); % 气网节点气压 Gsup = sdpvar(nGasSource, 1); % 气源产气量 GB = binvar(nTap, 1); % 变压器分接头离散选择 CB = binvar(nCap, 1); % 电容器投切离散选择 % 目标函数:发电成本 + 购气成本 + 碳排放成本 Objective = sum(a .* Pg.^2 + b .* Pg + c) ... + sum(priceGas .* Gsup) ... + sum(emisCoef .* Pg) * priceCO2; % 约束集合 Constraints = []; % 潮流方程(极坐标形式) for i = 1:nGrid Constraints = [Constraints, Pg(i) - Pd(i) == ... V(i) * sum(V .* (Gbus(i,:) .* cos(theta - theta(i)) ... + Bbus(i,:) .* sin(theta - theta(i))))]; Qg(i) - Qd(i) == ... % 同理 end % 气网约束,用Π变量和Weymouth方程 Pi = press.^2; % 或者直接用press上的二阶锥约束 for k = 1:nPipe Constraints = [Constraints, ... f(k)^2 == C(k)^2 * (Pi(start(k)) - Pi(end(k)))]; end % 耦合约束 Constraints = [Constraints, ... P_GT(n) == alpha_GT * GasConsume(k) + beta_GT]; % 燃气轮机 Constraints = [Constraints, ... GasProduce(m) == eta_P2G * P_P2G / HHV]; % P2G这套写法最大的优势是可读性,模型写的清楚,两天后再回来看代码仍能快速定位每一个约束对应的物理含义。调试时也方便:把某个约束注释掉,对比目标函数和变量值得变化,马上能判断“是哪条约束卡住了最优解”。
4.5 求解流程串起来怎么跑
主程序流程我用下面的顺序:
- 加载电网数据(Matpower格式)和马气网数据
- 做初值潮流,获取电压初值
- 用YALMIP定义决策变量和约束,构建MISOCP或NLP模型
- 设置求解器参数(Gurobi的MIPGap设为0.01%,fmincon的MaxIterations设高些)
- 求解,保存结果到结构体
- 校验物理可行性:把优化得到的注入重新代入潮流,检查电压是否越限、气网流量是否满足Weymouth方程
- 输出结果表格和图形
第6步很多人会忽略,我要特别强调:优化器给出的“可行解”有可能因为松弛不紧而物理上不成立,所以回代校验必不可少。我每次跑完优化都要校验一次,一旦发现问题,优先检查SOCP松弛的紧性,看互补间隙有多大。
5. 算例结果:协同优化的收益到底有多大
5.1 场景设计:单无功优化、有功-无功独立、协同优化三者对比
为了让协同价值“看得见”,我设置了三个场景做对比,算例是IEEE 14节点电网+6节点气网的组合系统,风电渗透率20%,典型冬季日负荷水平:
- 场景A:固定有功出力(经济调度结果),只优化无功——即传统无功优化;
- 场景B:先做有功经济调度,再做无功优化,两步顺序执行但不迭代,模拟“有功、无功分开做”的工程常见做法;
- 场景C:有功-无功协同优化,即本文搭建的联合模型,气网约束、燃气轮机出力和无功调节手段同时决策。
场景B是重要对照组。它在学术上常被当成“非协同”的代表做法——两阶段顺序求解,第二阶段最优解大概率落不进第一阶段想象的那个“全局最优”里。
5.2 结果对比表:网损、电压偏移、碳排放、购气成本
三个场景在优化结束后我用统一潮流校验器重新计算了一次网络状态,结果如下(这个数据针对上述特定算例,不代表所有系统的普遍结论,但趋势有代表性):
| 指标 | 场景A:传统无功优化 | 场景B:有功/无功两步顺序解 | 场景C:有功-无功协同优化 |
|---|---|---|---|
| 网损率(%) | 5.12 | 4.78 | 4.15 |
| 最低节点电压(pu) | 0.938 | 0.951 | 0.961 |
| 电压偏移总超标量(pu) | 0.042 | 0.019 | 0.008 |
| 总运行成本(万元/h) | — | 15.72 | 15.02 |
| 碳排放量(t/h) | — | 63.8 | 58.4 |
| 燃气轮机总出力(MW) | 112 | 112 | 137 |
场景A不涉及经济调度,所以成本和碳排放不在同一比较基准上。重点看场景B到场景C的变化:总成本下降约4.5%,网损下降约0.6个百分点,最低电压从0.951提升到0.961,碳排放减少约8.5%。这个结果说明,协同优化最大的收益不是单一指标,而是“总成本更低的同时,电压更安全”。
5.3 结果背后的机理:气网灵活性如何替电网“解重载”
数值结果只是表象,机理更重要。为什么协同优化能把网损和电压都拉下来?我用一个具体的机理来解释,这比表格更有说服力:
算例里电网的某个重负荷区域电压偏低,传统无功优化只能靠就地并联电容器和调变压器分接头来“硬撑”,效果有限。但在协同模型里,优化器发现:该区域燃气轮机的出力受气网节点气压约束,而P2G装置如果启动,可以向气网注入天然气,抬高该节点气压,从而解除燃气轮机的出力上限。
于是协同优化的最优策略是:启动P2G消耗富余风电制气→气网节点气压上升→燃气轮机增发有功→重负荷区域本地电源出力增加→减少了从远处输电线路送来的重载功率→网损下降、电压升高。
这个链条里,P2G的启动单价(耗电成本+运行成本)比远距离输电的网损成本更低,而且燃气轮机的清洁属性让碳排放也降下来了。整个逻辑在单独做无功优化的场景A里根本不会出现——因为在那个模型里,燃气轮机的有功出力被固定死,P2G也不存在。
所以我常说,协同优化的价值不在于把模型写得更复杂,而在于它让优化器看到了“气网→燃气轮机→电网电压”这条原本被截断的可行路径。
6. 代码复现避坑指南:我从调试中总结的几条经验
6.1 初值不收敛问题:从“改初值”到“找结构性原因”
调试这套模型时,我最常收到的求助就是“fmincon报错说在初始点约束不满足”,或者“潮流迭代发散”。大多数人第一反应是改初值,但我建议先找结构性问题。
最典型的例子:气网节点气压初值全部取基准值6.0 MPa,但管道流量计算时sqrt(|π_k² - π_m²|)如果两端气压相等,流量为零,雅可比矩阵里该管道对应的行列会出现奇异。解决办法不是硬调初值,而是让初值气压带一点“坡度”——比如节点1取6.0、节点2取5.95、节点3取5.9,让管道有一个可导的初始流量。这种物理上的疏导,比任何数值技巧都有效。
还有一类是约束过紧导致的问题:气源上限设小了,导致气网无法支撑燃气轮机在某时段的目标出力。此时优化器不是找不到局部最优,而是真的没有可行解——这和初值无关,必须先放宽物理约束或调整耦合元件的转换效率。
6.2 SOCP松弛不紧该怎么办
用MISOCP求解时,最担心的是二阶锥松弛后的最优解不满足原始非凸方程。判断和处理的流程如下:
- 求解完成后,把每个管道的原始Weymouth方程的残差算出来:r = |f_km² - C_km²(π_k² - π_m²)|;
- 计算相对残差r / f_km²,如果小于1e-4,认为松弛足够紧,结果可靠;
- 如果残差较大,说明该管道流量和气压之间有明显“空隙”,松弛后的可行域过大,最优解可能是个松弛伪解。
缓解手段按优先级排序:先检查气压上下限是否过宽——气压范围越宽,Weymouth凸包越松;然后尝试减小管道C_km的取值,让流量与气压差的耦合更紧;最后的手段是引入割平面,把不满足原始方程的松弛解从二阶锥集中剪掉。这条路径需要一点数值功底,但对结果质量影响极大。
6.3 天然气网络参数的单位陷阱
这一点我想单独拎出来讲,因为它耗费了我大约三天调试时间。天然气网络在不同的文献里使用的单位五花八门:压强用MPa、bar、psia的都有,流量用kg/h、m³/h、MMSCF/d、MW都有,管道系数C_km更是量纲各异。如果从不同文献里拼凑数据,单位不一致的问题几乎是必然的。
我的建议是所有内部计算统一到标幺值体系:压强基准取基准节点气压,流量基准取全网总气源的额定产能,这样Weymouth方程里的C_km就成了一个无量纲的等效导纳系数,从量纲陷阱里彻底跳出来。每一处单位换算都在表格或注释里写清楚原始单位和换算公式,这不仅是给自己看的,也是给未来几个月后重新打开代码的自己看的。
6.4 代码运行效率优化
这套模型在小算例(14节点电网+6节点气网)上跑得很快,但如果扩展到IEEE 118节点或比利时20节点气网,性能就会成为瓶颈。我常用的加速手段:
- 稀疏矩阵:Ybus、雅可比矩阵显式用sparse,求解线性方程组时的加速非常明显;
- 减少YALMIP变量数量:比如管道流量用气压变量的表达式表示,而不是单独定义流量变量再附加等式约束,这样既少了变量又少了约束,求解速度提升明显;
- 固定整数变量预求解:先用连续松弛解跑一遍,把明显不会动作的分接头和电容器组固定下来,缩小整数搜索空间,再用MISOCP精解。这样可以减少一到两个数量级的求解时间。
另外,Matlab的并行计算工具箱可以用于多场景并行扫描(比如做24小时滚动优化),每个时段独立求解时用parfor并行,效果立竿见影。
最后聊一点我个人的使用体会。这套电气互联系统有功-无功协同优化模型,我在跑了大量算例之后最大的感受是:模型的价值不在于把两个独立的优化问题“合并”到一个文件里,而在于真正理解了“有功调节影响无功可行域、气网约束影响有功可行域”这条传导链。刚开始做代码实现时,我也是一步一步从纯电网无功优化起步,先把AC潮流、无功补偿、变压器分接头这些基本模块搭好,再加上气网和耦合元件,最后才形成完整的协同模型。如果你正在入门,我的建议是不要一开始就追求大算例,先把这套小系统跑透、跑明白,再去扩规模——因为在复杂模型上排查一次约束问题,远比在小模型上理解透机理要痛苦得多。