☰
基于二阶锥松弛的主动配电网故障重构建模与MATLAB实现
2026/10/14 7:11:37 网站建设 项目流程

配电网和输电网的故障处理有一个特别大的区别:输电网可以靠多重冗余路径和就地保护完成自动恢复,而配电网一旦上游故障隔离,下游往往直接失压,能不能快速恢复供电,很大程度上取决于你能否在短时间内拿出一套"把非故障区转移出去"的可行方案。这件事在行业里的正式名字叫配电网故障重构。我第一次被这个问题彻底“卡住”,是在一个接入了多台分布式光伏和储能的小型试验馈线上:故障发生在干线上游,按老经验把联络开关合上倒是简单,可合上去之后电压越限、部分线路过载,反而比不复电更麻烦。后来我采用基于二阶锥松弛的主动配电网故障重构方法,在MATLAB里建模求解,几分钟内就得到了既有拓扑切换又有分布式电源出力调节的联合恢复方案。这篇文章想把从“看懂公式”到“跑通代码”再到“能向别人解释结果”的完整过程写出来,适合正在做配电网重构研究的同学,也适合配电自动化工程师快速理解这类方法的落地路径。

1. 故障重构到底做的是什么?先厘清问题边界

1.1 重构的优化本质:一份开关组合与运行方式的联合决策

配电网故障重构的完整含义是:故障发生后,通过保护装置隔离故障点,系统内部分分段开关断开、部分联络开关闭合,把非故障失电区域重新挂接到其他馈线,使全网恢复辐射状供电。重构策略不只是回答“哪些开关合、哪些开关断”,还要回答“合上之后电压合格吗、线路会不会过载、分布式电源出力怎么调”。

这个问题的数学形式是一个优化决策模型,目标函数通常是多个指标的加权和:减供负荷量最小(最好为0)、网络损耗最小、开关操作次数最少、节点电压偏差最小。决策变量分两类:一类是表示支路投切状态的0-1变量,另一类是连续的运行变量,包括支路有功/无功潮流、节点电压幅值、分布式电源的有功/无功出力。约束条件由潮流方程、节点电压上下限、支路容量、辐射状连通结构、DG可调范围共同构成。

这里有个很容易被忽略的点:配电网重构和输电网重构的诉求不一样。输电网重构关心功角稳定、断面潮流这些系统级问题;配电网本来就是辐射状或弱环网运行,所以必须强制约束重构后的网络仍然是辐射状,否则闭环会改变故障电流水平、影响保护配合。这也是为什么重构模型里“辐射状约束”是不可缺少的一环,后面建模章节会专门讲。

1.2 主动配电网带来的变量:重构从“拓扑切换”变成“拓扑+运行协同”

传统配电网故障重构可以近似认为负荷是固定不变的,决策重心就是开关组合。但主动配电网里加入了光伏、风电、储能、无功补偿装置,情况完全不同。

举一个很典型的场景:某条馈线故障后,把失电负荷转移到相邻馈线,如果转移过去的负荷对应区域里有大量分布式光伏,倒送功率可能过大,或者某个节点电压越过上限。反过来,如果DG容量不足,又可能撑不住电压,末端电压偏低。此时“怎么调整DG出力”和“怎么切开关”必须放在同一个优化问题里联合求解,不能像传统做法那样先定拓扑再调运行。

正因如此,主动配电网故障重构的求解难度比传统版本高一个量级:拓扑变量是离散的0-1变量,潮流方程是非线性的,两者耦合在一起,直接构成混合整数非线性规划(MINLP)。这类问题在数学上通常是NP-hard的,节点数一多,穷举开关组合根本不现实。33节点系统的联络开关组合已经有数万种,工程上几十上百节点的系统更没法直接枚举。这也是为什么一定要引入二阶锥松弛——把非凸问题转成凸问题,才能用通用求解器在可接受时间内算出高质量解。

1.3 为什么不能直接枚举开关组合,再逐个算潮流?

有人会问:既然开关组合是有限的,为什么不能全列出来然后挑一个最好的?关键在于“检验一个开关组合好不好”这个步骤不便宜——每组开关都要做一次潮流计算,而潮流本身在弱环网和有DG情况下需要迭代求解;更麻烦的是,一旦DG出力可变,即使开关组合固定,运行方式仍然连续可变,最优运行点还得再解一个优化问题。所以“枚举-验算”在规模稍大时根本不现实。

比较现实的路径是:把开关状态建模成0-1变量,把潮流关系建模成约束,把运行指标写成目标函数,然后交给成熟求解器搜索。这个思路在学术界和工业界已经是主流。二阶锥松弛的价值就在这个框架里体现:它把非凸的潮流非线性约束转换为凸锥约束,使问题最终成为一个混合整数二阶锥规划(MISOCP),商业求解器对这类问题通常有稳定的全局求解能力。

2. 从MINLP到SOCP:二阶锥松弛的建模推导与精确性

2.1 DistFlow潮流方程:配电网优化里最常见的简化模型

在配电网重构建模中,最常用的潮流模型是DistFlow分支潮流方程。它假设网络是辐射状的,只考虑每条支路的功率流和电压降,避开完整牛顿-拉夫逊迭代,形式简洁,适合嵌入优化问题。

考虑一条支路(i, j),令 P_ij、Q_ij 分别为从节点i流向节点j的有功和无功功率,令 U_i = V_i²(电压幅值平方的变量代换)、l_ij = I_ij²(电流幅值平方),DistFlow方程可以写为:

  • 节点j的有功功率平衡:P_j = P_ij - r_ij * l_ij - Σ P_jk
  • 节点j的无功功率平衡:Q_j = Q_ij - x_ij * l_ij - Σ Q_jk
  • 支路电压降方程:U_j = U_i - 2(r_ij * P_ij + x_ij * Q_ij) + (r_ij² + x_ij²) * l_ij
  • 电流幅值等式:l_ij = (P_ij² + Q_ij²) / U_i

前三个方程在变量代换后是线性的,最后一个等式是非线性的。也就是说,DistFlow本身已经把“电压类”和“电流类”二次量都替换成了新变量,剩下一个电流平方的等式还带着非凸性。接下来的思路就非常清晰了:把非线性等式松弛成不等式,再进一步改写为凸锥。

2.2 松弛与锥变换:用一个“≥”换来整个问题的凸性

将最后一个等式 l_ij = (P_ij² + Q_ij²) / U_i 松弛为不等式:

l_ij ≥ (P_ij² + Q_ij²) / U_i

这个不等式等价于:

(2P_ij)² + (2Q_ij)² + (l_ij - U_i)² ≤ (l_ij + U_i)²

写成标准二阶锥形式:

|| (2P_ij, 2Q_ij, l_ij - U_i) || ≤ l_ij + U_i

在YALMIP里可以直接用cone函数写,也可以用norm(vals) <= y的形式。这个约束的几何意义是:把一个三维向量限制在一个二阶锥内部,这个集合是凸的。原本非凸的MINLP潮流部分,就此变成凸约束,整体问题变为MISOCP。

为什么可以放心地把等式“松弛”成不等式?关键在于目标函数。如果目标函数里有让 l_ij 尽可能小的方向(比如最小化网络损耗),最优解会有动力把 l_ij 压到不等式下界上,从而“自动”恢复等式成立。这种情况在文献里叫做精确松弛。当然,精确性不是无条件成立的,具体条件后面专门聊。

2.3 辐射状约束与开关逻辑怎么进模型

辐射状约束是重构建模中最容易写错、写漏的一块。最常用的写法是:

  1. 边数约束:闭合支路数等于节点数减1,即 Σ z_ij = N - 1。
  2. 连通性约束:单靠边数约束不能保证连通,可能出现多个互不相连的环路。常见做法是引入虚拟“商品流”,要求每条闭合支路能向子节点传播单位流量,源点提供N-1单位虚拟流量。这样既保证连通性,也能保证无环。

实际建模时,我通常用“边数约束+虚拟流约束”。虚拟流变量记为 f_ij,要求:

  • |f_ij| ≤ M * z_ij,M取节点数,保证断开支路上虚拟流为0;
  • 每个节点的净流出为1(源节点特殊处理)。

开关逻辑约束的典型处理是大M法。对每条支路,如果开关断开(z_ij=0),支路潮流必须为0;如果闭合,潮流在容量范围内自由变化:

  • -bigM * z_ij ≤ P_ij ≤ bigM * z_ij
  • -bigM * z_ij ≤ Q_ij ≤ bigM * z_ij

电压降方程也需要用大M“接管”:支路断开时,U_i 和 U_j 之间不要求满足DistFlow方程;闭合时才要求成立。具体做法是把电压降约束拆成两部分,一部分是闭合时的DistFlow电压降,另一部分是带 z_ij 的松弛大M。这部分是代码里最容易出bug的地方,后面会专门讲。

2.4 松弛的精确性:SOCP解什么时候就是原问题的解?

SOCP松弛的精确性不是免费赠送的。辐射状网络、正常运行电压范围、目标函数单调性这些条件,都会影响“松弛紧不紧”。研究文献里最常见的结论是:在辐射状网络、目标函数关于电流单调递增的条件下,二阶锥松弛通常精确;但在电压幅值被约束到系统下限附近、或者网络存在多重连通可能时,可能出现松弛后的最优解对应一个不可实现的运行点。

从工程实操角度,我一般会在求解完成后做一个后验校验:把SOCP得到的最优 P_ij、Q_ij、U_i 代回原始等式 l_ij = (P_ij² + Q_ij²) / U_i,检查每条闭合支路的紧性残差。如果某个支路的残差超过1e-4量级,说明松弛可能不紧,需要检查目标函数设置,或考虑加入电压幅值惩罚等修正手段。这一步非常重要,很多论文里“算例验证松弛紧性”做的就是这个事。

3. MATLAB+YALMIP建模代码:把数学模型落成可求解脚本

3.1 算例数据准备:先把系统参数组织成矩阵

我使用的是33节点测试系统,含32条支路和5条联络开关,初始状态是正常辐射状运行。为了演示故障场景,我设置了节点7与节点8之间发生永久性故障,故障隔离后节点8到18这一大段负荷失电。

先说数据组织的思路。你需要准备四类数据:支路表、负荷表、DG表、运行参数。支路表里包含首端节点、末端节点、电阻、电抗、初始开合状态;负荷表里包含各节点有功和无功负荷;DG表里包含DG接入节点、有功上限、无功上下限;运行参数包含电压上下限、支路容量。这些数据在MATLAB里建议全部用结构体保存,引用方便,排错也容易。

3.2 变量定义与二阶锥约束的YALMIP写法

YALMIP是MATLAB下的优化建模工具箱,最大的好处是可以把优化模型写得接近数学公式,后端再切换不同求解器。对于MISOCP,它需要调用支持混合整数锥规划的求解器,比如Gurobi、CPLEX、MOSEK。

变量定义部分:

% 系统规模 nNode = 33; % 节点数量 nBranch = 37; % 支路数量(含联络开关) % 决策变量 z = binvar(nBranch, 1); % 支路开关状态,1闭合,0打开 P = sdpvar(nBranch, 1); % 支路有功潮流 Q = sdpvar(nBranch, 1); % 支路无功潮流 U = sdpvar(nNode, 1); % 节点电压幅值平方 l = sdpvar(nBranch, 1); % 支路电流幅值平方 f = sdpvar(nBranch, 1); % 虚拟流量,用于辐射状约束 Pg = sdpvar(nDG, 1); % DG有功出力 Qg = sdpvar(nDG, 1); % DG无功出力

构建DistFlow约束时,需要对照方程逐条写。以节点功率平衡为例,方向符号非常容易搞混。我习惯的统一约定是:支路功率从首端节点流向末端节点,节点功率平衡写成“流出 = 负荷 - 注入”。

Constraints = []; for i = 1:nNode % 流出节点的有功功率 outflowP = 0; inflowP = 0; for k = 1:nBranch if branchFrom(k) == i outflowP = outflowP + P(k); end if branchTo(k) == i inflowP = inflowP + P(k); end end % 节点净流出 = 负荷 - DG注入(这里需要根据负荷和DG数据组织好) netLoadP = loadP(i) - dgP(i); Constraints = [Constraints, outflowP - inflowP == netLoadP]; end

这段代码写成循环是为了直观,实际工程中可以向量化。有一点要提醒:节点功率平衡里的正负号,必须和你定义的潮流方向一致。如果首末端搞反,模型会直接不可行,而且排查起来很费劲。

二阶锥约束可以直接写:

for k = 1:nBranch vals = [P(k); Q(k); l(k) - U(branchFrom(k))]; Constraints = [Constraints, cone(vals, l(k) + U(branchFrom(k)))]; end

我实测下来,cone(vals, y)在多个求解器上都能正确识别为二阶锥约束;直接写norm(vals) <= y也等价,但部分求解器内部转换效率会有细微差异。支路多的时候,用循环逐条添加Constraints即可。

辐射状约束里,虚拟流部分可以写成:

% 每条支路的虚拟流只在开关闭合时存在 Constraints = [Constraints, -nNode * z(k) <= f(k) <= nNode * z(k)]; % 节点平衡:每个节点的净流出为1,源点特殊 for i = 1:nNode outF = 0; inF = 0; for k = 1:nBranch if branchFrom(k) == i outF = outF + f(k); end if branchTo(k) == i inF = inF + f(k); end end if i == sourceNode Constraints = [Constraints, outF - inF == nNode - 1]; else Constraints = [Constraints, outF - inF == 1]; end end

加上闭合支路总数约束:

Constraints = [Constraints, sum(z) == nNode - 1];

到这里,辐射状约束就完整了。注意sum(z) == nNode - 1这个约束千万别漏,否则即使虚拟流第一眼看不出问题,也可能解出多个环。

3.3 完整求解流程与结果提取

求解调用非常简洁:

options = sdpsettings('solver', 'gurobi', 'verbose', 2); optimize(Constraints, Objective, options); % 结果提取 P_opt = value(P); Q_opt = value(Q); U_opt = value(U); l_opt = value(l); z_opt = round(value(z)); % 0-1变量有时返回0.999/0.001,就近取整

这里要特别提示:value(z)得到的结果在数值上可能出现0.99999或0.00001,不要直接拿去展示拓扑,先round再校验。很多初学代码里的奇怪拓扑,就是没取整导致的。

3.4 求解器选择:不是每个求解器都能吃MISOCP

一个容易踩的坑是:MATLAB自带的intlinprog不支持二阶锥约束,fmincon也不支持0-1变量与锥约束的混合问题。跑MISOCP,你需要Gurobi、CPLEX、MOSEK这类商业求解器,或者至少支持二阶锥的开源求解器。YALMIP的好处是后端可以随意切换,换求解器只改一个sdpsettings选项。

我把常用求解器和特点整理成表格:

求解器MISOCP支持许可证我的实测感受
Gurobi支持商用,有学术授权收敛最稳,MIP搜索能力强
CPLEX支持商用,有学术授权同样稳定,部分锥问题略逊
MOSEK支持商用,有学术授权连续锥问题非常强
SCIP部分支持开源能算,速度有差距

我的建议很简单:有学术授权就优先Gurobi;没有商业授权,先在小算例上跑SCIP验证模型正确性,再找有授权的求解器资源来算大规模场景。

4. 求解过程中最容易翻车的四个环节:数值、参数与收敛性

4.1 大M取值:多少算大,多少算小?

大M法看起来简单,取值却直接影响求解稳定性。取1e6甚至1e8时,约束矩阵条件数变差,求解器可能把本该有明显数值差异的约束混为一谈,出现数值病态;取太小又可能把可行域错误剪掉,甚至直接得出不可行的结论。

我的经验是分两类取:

  • 支路潮流大M:取该支路容量上限的3-5倍,或者全系统最大负荷的2倍。对33节点系统,100这个量级足够。
  • 电压降方程大M:取(Vmax² - Vmin²)的上限即可,不需要再放大。

我习惯的做法是:先用稍大的M跑一遍,检查最优解里哪些约束的潮流值接近bigM边界;如果某个约束被打到边界,说明M取得太低,或者该约束建模本身有问题,需要修正。

4.2 MIPGap与TimeLimit:重构问题的实用平衡点

MISOCP本质上是MIP问题,Gurobi在求解过程中会不断更新可行解和目标上下界。默认的MIPGap是1e-4,这在某些重构场景下可能导致求解时间过长。

实操中我会先把TimeLimit设成120秒,MIPGap设成0.01,也就是1%,先拿一个工程可用的可行解;如果后续想把结果做得更精细,再把MIPGap收紧到1e-3或1e-4,同时放宽时限。逻辑很简单:配电网故障重构在实际调度场景里对时间敏感,一个1%以内的次优拓扑,往往比等十几分钟才拿到的全局最优更有工程价值。尤其是实时故障处理时,“快”本身就是硬约束。

4.3 infeasible的排查链路:先注释哪一条约束?

模型不可行时,求解器会返回Infeasible,但怎么定位写错的那条约束才是难点。

我的排查路径是固定的:

  1. 先把开关约束和辐射状约束全部注释掉,只保留原拓扑下,也就是z_ij固定为1的连续潮流约束,算一次连续SOCP。如果这个都不可行,说明DistFlow方程或节点功率平衡有明显问题,先检查正负号和负荷数据。
  2. 还原开关约束,注释掉辐射状约束,看求解结果。如果不可行,问题极可能出在大M约束或潮流方程与开关状态的耦合上。
  3. 最后还原辐射状约束。如果这时才不可行,再检查虚拟流约束和边数约束的组合是否有冲突。

这个方法本质是一个“约束二分查找”。我靠这个思路多次在几分钟内定位到不可行原因,比对着报错信息瞎猜效率高很多。

4.4 锥约束不紧的后验校验

前面提到过“后验校验”,这里给出具体代码:

% 计算每条闭合支路的紧性残差 U_from = U_opt(branchFrom); residual = l_opt - (P_opt.^2 + Q_opt.^2) ./ U_from; % 只看闭合支路 active = find(z_opt > 0.5); maxAbsResidual = max(abs(residual(active)));

如果maxAbsResidual在1e-6量级或更小,可以认为松弛精确;如果达到1e-3或更大,就要小心。一个常见原因是:目标函数里损耗项系数为零或过小,或者电压幅值约束把解推到了锥边界之外,导致锥松弛不紧。处理办法是在目标函数里加一个小权重乘以所有支路的l_ij之和,或者适当放宽电压约束边界,让问题有空间回到精确解。

5. 可视化:从矩阵结果到能直接汇报的图谱

5.1 拓扑切换对比图是最重要的第一张图

光有数据和表格,别人很难直观理解重构的价值。我强烈建议第一步先把拓扑画出来。

MATLAB里可以用graph对象快速画网络。把闭合支路和打开支路分别提取,用实线、虚线和颜色区分,再用plot配合Layout选项自动分布节点,效果比手动坐标好很多。重构前后两个状态并排放在同一张figure里,视觉对比非常清楚。

G = graph(fromNode(active), toNode(active)); p = plot(G, 'Layout', 'force', 'LineWidth', 2); highlight(p, G.Edges.EndNodes, 'NodeColor', 'r');

这里我习惯把故障隔离支路单独标红,把重构后新合上的联络开关用绿色虚线标出,一张图就能说完整个调度方案。

5.2 电压剖面对比:调度最关心的指标

电压是配电运行人员最关心的指标之一。重构前失电区域没有电压,恢复后各节点电压可能偏低或偏高。我会画一个横轴为节点编号、纵轴为电压幅值的折线图,把重构前、重构后、电压上限、电压下限四条线放一起。

用U_opt开方得到电压幅值,注意标出越限点。这张图在汇报时几乎是必放的一张,尤其当你想说明“这个重构方案确实解决了电压问题”的时候。

5.3 支路功率热力图与开关动作表

除了拓扑和电压,我还喜欢把支路功率用线宽和颜色表达在图上。闭合支路的线宽正比于支路视在功率,颜色从绿色到红色渐变表示负载率。这样一眼就能看出负载转供后的热点路段,帮助运行人员判断安全裕度。

开关动作表可以用表格输出,列字段:支路编号、首端节点、末端节点、初始状态、最终状态、是否动作。配电网调度人员更关心“动了哪几台开关”,而不是所有支路细节。输出时把动作的开关排在最前面,用醒目颜色标出即可。

5.4 图形导出与报告整理

我习惯在画图结束后统一导出300 dpi的PNG,方便嵌入报告:

exportgraphics(gcf, 'reconfiguration_result.png', 'Resolution', 300);

如果你的MATLAB版本较老,可以用print('-dpng', '-r300', ...)。导出前把图例、字体大小调好,避免后续反复重导。

6. 往工程落地走:扩展方向与我的实操建议

6.1 从单一故障扩展到N-1预想故障集

实际系统的重构不只是针对一个已经发生的故障,还经常要应对“预想故障集”。为了让运行方式更有弹性,离线阶段要把多个典型故障场景分别求解,得到一套“开关预案表”。

这个扩展在建模上几乎没有改动,把故障场景参数封装成一个函数,循环调用即可。真正要花心思的是:不同场景下的DG出力场景怎么设定、故障隔离范围怎么自动识别。这部分工程化工作比重往往比建模本身还大,我在实际推进中深有体会。

6.2 考虑DG不确定性的鲁棒重构思路

当DG渗透率进一步提高,固定出力假设会让重构结果在实际运行中吃不消。一个自然的扩展是把DG出力当成不确定参数,构造盒式鲁棒重构模型,或者用场景法描述多组可能的出力曲线来做随机重构。

此时SOCP基础模型仍然保留,只是在不确定约束处理上多出一层对偶变换或场景复制。代码量会大不少,但只要把基础SOCP重构理解扎实,扩展起来并不难。

6.3 给初学者的学习路径与操作习惯

如果让我带人,路径会是这样:

  • 第一步:在33节点系统上先别管开关变量,把DistFlow和SOCP连续模型跑通,和潮流计算的电压结果比对,确保物理正确。
  • 第二步:固定一条故障支路,加开关变量和大M约束,跑MISOCP,目标只设为“负荷恢复+辐射状”,不要一上来就整多目标。
  • 第三步:逐步加损耗目标、电压偏移目标、DG出力调节。
  • 第四步:再做可视化。

这个路径能避开大量弯路。重点在于第一步一定不要跳:很多同学一上来就建全模型,结果不可行时,分不清是潮流方程错了还是开关约束错了。两步并一步走,排查成本成倍上升。

最后再分享一个实操习惯:跑这类MISOCP时,每改一次数据或约束,我都会记录当次的求解时间、目标值和关键约束残差,汇总到一个表格里。时间久了,你能很清楚地总结出哪类参数组合最稳定、哪类故障最难解。这种积累,比任何教程都有用。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询