凌晨两点调度电话响起来10kV馈线跳闸重合失败故障点隔离后下游还有一片用户黑着。这时候最要紧的不是先查故障原因而是尽快给出一个“哪些开关合、哪些开关断”的供电恢复方案。这个决策背后就是配电网故障重构。我这两年复现和落地这类问题最顺手的组合是Matlab Yalmip搭模型把重构问题写成二阶锥规划SOCP来解稳定、快、还不容易掉进局部最优。这篇就把这套实现从头到尾掰开揉碎讲一遍问题怎么建模、二阶锥松弛为什么是首选、Yalmip代码骨架怎么搭、求解器怎么选以及我实际跑通之后踩过的几个坑。1. 故障重构到底在优化什么不只是“倒闸”1.1 问题边界故障隔离后的拓扑重构配电网和输电网最大的区别在于配电网基本都是辐射状结构闭环设计、开环运行。正常状态下网络是一棵树故障发生后保护装置把故障支路两侧的开关断开故障点被隔离但原本由这条支路往下送电的所有负荷会全部失电。这时候调度员能做的就是通过操作联络开关和分段开关把失电区域重新挂到别的馈线或者同一个变电站的其他出线上。故障重构问题的输入边界很明确给定一个故障前拓扑已知故障支路已经被强制隔离这一条支路在重构过程中不可恢复你需要找到一组新的开关状态使得整个网络的潮流满足运行约束同时让网损、开关操作次数、未恢复负荷量这些指标尽可能好。它本质上是一个带二进制变量的大规模混合整数优化问题——开关只有闭合和断开两种状态一旦加了潮流约束就成了混合整数非线性规划。但很多人第一次接触时容易把问题想窄了觉得“故障重构把联络开关合上就行”。实际上同时要考虑三件事第一合上联络开关会不会造成环网配电网闭环运行会导致保护配合混乱必须始终保持辐射状第二恢复后的电压是否越限远端节点很容易电压偏低第三线路容量能不能撑住负荷转移很多时候初始方案能恢复供电但某些支路过载照样不可行。1.2 目标函数与约束的完整画像一个完整的重构模型目标函数通常长这样min 网损 λ1 × 开关操作次数 λ2 × 未恢复负荷量网损所有闭合支路的电流平方乘以电阻之和公式上就是 Σ r_ij * L_ij。开关操作次数反映重构动作的代价。少动开关意味着倒闸快、风险低所以一般用重构后的开关状态与故障前状态做差的绝对值之和。未恢复负荷量如果系统容量足够这个量应该是0但如果发生严重故障或N-1不满足可能需要切除一部分负荷就要通过大权重强制优先恢复重要负荷。约束条件更琐碎但每一条都有实际物理含义约束类型表达式物理含义节点功率平衡Σ注入功率 - Σ流出功率 接入负荷每个节点的有功/无功必须守恒支路电压降支路两端电压满足欧姆定律线路压降不能随意设定二阶锥潮流约束v_i × L_ij ≥ P_ij² Q_ij²关联电压、电流和功率之间的非线性关系电压上下限V_min² ≤ v_i ≤ V_max²保证供电质量辐射状拓扑约束闭合支路数 节点数 - 1且全网连通避免环网与孤岛开关状态约束故障支路固定为0故障隔离不可恢复把这些东西全部写进模型才是一个可求解的完整问题。注意很多初学者写的模型只有潮流约束和开关变量忘了辐射状约束结果解出来是个环网还得靠人工后处理。我后面在Yalmip里实现的辐射状约束会结合虚拟潮流和支路数双重校验目前用下来比单独加一个“闭合支路数N-1”可靠得多。2. DistFlow的非凸困境与二阶锥松弛的价值2.1 DistFlow方程改写把潮流约束“翻译”成变量关系配电网潮流计算里用得最多的是DistFlow模型支路潮流模型。它的核心方程写出来并不复杂对每条支路 (i, j)功率平衡流入节点j的功率等于节点j的负荷加上从j流向所有子支路的功率和损耗电压降节点j的电压平方等于节点i的电压平方减去线路上压降带来的影响电流定义支路电流平方等于支路有功平方和无功平方之和除以节点电压平方如果直接把这些方程扔给求解器麻烦出在第三个方程它是电流、电压平方和功率之间的非线性隐式关系数学上是一个非凸的等式约束。非凸约束在优化里意味着多解、初值敏感、局部最优传统做法用遗传算法、粒子群这一类启发式算法去搜搜一次几十秒运气不好还搜不出可行解。有一个经典的改写技巧定义新变量把原来的物理量平方化。v_i U_i² l_ij I_ij²代入以后电压降方程从带平方根的复杂形式变成了关于 v、l、P、Q 的线性方程。这是DistFlow模型最漂亮的地方——把大部分非线性整理成了好处理的形态。但电流定义方程改写以后变成v_i × l_ij P_ij² Q_ij²左边是电压平方乘以电流平方右边是功率平方和。问题来了左边两个变量相乘这在优化里是双线性项依然是非凸的。这一步卡住了很多人。2.2 从双线性等式到二阶锥一步松弛的代价与边界关键转折点在于把上面的等式约束直接松成一个不等式。v_i × l_ij ≥ P_ij² Q_ij²这个不等式的几何意义是把一个抛物面非凸等式松弛成了一个二阶锥空间。二阶锥规划是凸优化里非常成熟的一类问题配合整数变量以后是MISOCP混合整数二阶锥规划商业求解器对付这种问题非常拿手。也许你会问松弛不等于原问题求出来的解能信吗这就要提到“精确松弛”的概念。在辐射状配电网中只要网络拓扑固定、负荷条件合理SOCP松弛在很多情况下是紧的也就是说松弛后的最优解恰好满足那个等号成立。你可以把它类比成“用一根直线去逼近一条曲线在绝大多数位置都贴得很紧”。当然不是所有情况都紧所以有经验的工程师在求解完成后都会主动检查松弛间隙不检查等于白做。这个坑我在第5章详细说。Yalmip里写这个锥约束并不需要手动推导出标准锥形式直接用norm写法即可% 标准二阶锥约束|| [2P; 2Q; v_i - l] || v_i l Constraints [Constraints, norm([2*P(k); 2*Q(k); v(fb(k))-l(k)], 2) v(fb(k)) l(k)];这一行的数学本质其实是在说“它们的取值被控制在一个圆锥曲面之内”而不要求严格落在这个曲面上。在配电网重构这个场景里只要求解结果把锥约束算到等号附近说明模型写对了。3. Yalmip建模实战一个33节点系统的完整代码骨架3.1 数据组织接线矩阵、支路参数与故障场景推荐用IEEE 33节点算例做测试平台它是配电网重构、故障恢复领域公认的标准算例。数据组织方式直接影响后面代码复杂度我习惯用四个数组fb每条支路的首端节点编号tb每条支路的末端节点编号r每条支路的电阻标幺值x每条支路的电抗标幺值另外还需要一个向量z0记录每条支路在故障前的闭合状态1表示闭合0表示断开。正常运行的33节点系统32条分段开关支路是闭合的5条联络开关支路是断开的。故障场景就用一个故障支路编号来模拟比如支路5-6发生永久性故障那么这条支路的z强制为0且在优化过程中不可操作。这里给一个数据格式的示例代码片段实际加载时直接从数据库或Excel读即可% 支路数据示例fb, tb, r(标幺), x(标幺), z0 branch_data [ 1, 2, 0.00922, 0.00470, 1; 2, 3, 0.00493, 0.00251, 1; % ... 中间略 ... 8, 21, 0.03420, 0.01690, 0; % 联络开关正常断开 ]; fb branch_data(:,1); tb branch_data(:,2); r branch_data(:,3); x branch_data(:,4); z0 branch_data(:,5);一定要全部转成标幺值再进模型。我见过不少直接把欧姆、安培、瓦特混在一起丢进Yalmip的结果要么求解严重病态要么收敛到离谱的答案。用标幺值之后节点电压在1.0附近支路电流也是0到几之间Gurobi和Cplex处理这种数值尺度最舒服。3.2 决策变量声明与约束块逐段拆解我们需要的决策变量包括z binvar(nbranch, 1); % 每条支路是否闭合 v sdpvar(nbus, 1); % 节点电压平方 l sdpvar(nbranch, 1); % 支路电流平方 P sdpvar(nbranch, 1); % 支路有功 Q sdpvar(nbranch, 1); % 支路无功 f sdpvar(nbranch, 1); % 虚拟潮流用于辐射状约束接着按约束块依次写。第一条是故障支路强制断开Constraints [Constraints, z(fault_branch) 0];第二条是节点功率平衡。对所有非电源节点所有流入该节点的支路功率减去流出该节点的支路功率等于该节点接入的负荷。负荷可以写成可切形式用一个二进制变量alpha表示负荷是否接入alpha binvar(nbus, 1); alpha(1) 1; % 电源节点始终接入 for n 2:nbus inflow_branches find(tb n); outflow_branches find(fb n); Constraints [Constraints, ... sum(P(inflow_branches)) - sum(P(outflow_branches)) alpha(n) * P_load(n)]; Constraints [Constraints, ... sum(Q(inflow_branches)) - sum(Q(outflow_branches)) alpha(n) * Q_load(n)]; end这里P_load(n)和Q_load(n)是节点n的有功和无功负荷标幺值。alpha为0意味着该节点负荷被切除功率平衡只要求注入等于流出相当于该节点变成一个纯转供节点。第三条是电压降约束。闭合支路必须满足DistFlow的电压降方程断开支路必须把这个约束松弛掉。这里用大M法处理for k 1:nbranch i fb(k); j tb(k); % 闭合时满足电压降方程 Constraints [Constraints, ... v(j) - v(i) 2*(r(k)*P(k) x(k)*Q(k)) - (r(k)^2 x(k)^2) * l(k) M_big * (1 - z(k))]; Constraints [Constraints, ... v(j) - v(i) 2*(r(k)*P(k) x(k)*Q(k)) - (r(k)^2 x(k)^2) * l(k) -M_big * (1 - z(k))]; % 二阶锥约束只对闭合支路有意义断开时功率为0锥自然满足 Constraints [Constraints, ... norm([2*P(k); 2*Q(k); v(i) - l(k)], 2) v(i) l(k)]; % 断开支路时支路功率和电流平方强制为0 Constraints [Constraints, -M_power * z(k) P(k) M_power * z(k)]; Constraints [Constraints, -M_power * z(k) Q(k) M_power * z(k)]; Constraints [Constraints, 0 l(k) M_power * z(k)]; end注意M_big和M_power不能随便取。M_big取1到10就够因为标幺电压本身在1附近压降绝对值远小于1M_power取系统最大视在功率的两倍以上即可但也不要大到1e6那种程度否则会拖慢MISOCP求解速度。电压上下限、变电站电压给定这些约束可以直接写V_min 0.95^2; V_max 1.05^2; Constraints [Constraints, V_min v V_max]; Constraints [Constraints, v(1) 1.0]; % 根节点电压恒定3.3 辐射状约束虚拟潮流避免环网辐射状约束是新手翻车率最高的地方。单独写一个“闭合支路数等于节点数减一”是必要条件但不是充分条件——一个“树孤岛”的网络也满足这个数量关系但孤岛内的节点根本没接上电源。所以必须在支路数约束之外再加连通性约束。我采用虚拟潮流方法让每个非电源节点向根节点发送1个单位的虚拟流只有闭合支路能传输虚拟流。数学上写成% 闭合支路数约束 Constraints [Constraints, sum(z) nbus - 1]; % 虚拟潮流上下界只有闭合支路能传流 Constraints [Constraints, -nbus * z f nbus * z]; % 每个非根节点发送1单位虚拟流 for n 2:nbus inflow_branches find(tb n); outflow_branches find(fb n); Constraints [Constraints, ... sum(f(inflow_branches)) - sum(f(outflow_branches)) 1]; end如果某个非根节点和根节点之间没有闭合路径它的虚拟流平衡方程就不可能满足整个问题不可行。反过来如果闭合支路数正好是N-1且每个节点都能向根节点传流那这个网络一定是一棵连通的生成树。这个写法的好处是不引入额外的大量整数变量只增加了一个连续向量f和一组线性约束对求解非常友好。3.4 目标函数与求解器设置目标函数把三部分加权组合起来。网损在标幺制下数值很小开关操作次数和未恢复负荷量要大得多所以权重设置上我建议手动归一化而不是拍脑袋给权重% 网损 Loss sum(r .* l); % 开关操作次数二进制变量直接平方等价于绝对值因为0/1的平方还是0/1 SwitchOps sum((z - z0).^2); % 未恢复负荷量 Unsupplied sum((1 - alpha) .* (P_load Q_load)); % 权重设置 lambda1 0.1; % 开关操作权重按基准值归一化 lambda2 100; % 失电惩罚必须远大于网损优先保证负荷恢复 Objective Loss lambda1 * SwitchOps lambda2 * Unsupplied; % 求解 ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; result optimize(Constraints, Objective, ops);求解结束后至少做三件事验证结果可靠性% 1. 检查求解状态 if result.problem 0 disp(求解成功); else disp(求解失败/不可行); end % 2. 提取重构后的拓扑 z_sol value(z); fprintf(闭合支路数: %d\n, sum(z_sol)); % 3. 检查二阶锥松弛间隙 % 理想情况下 v_i * l_ij 严格等于 P_ij^2 Q_ij^2 gap zeros(nbranch, 1); for k 1:nbranch if z_sol(k) 0.5 i fb(k); gap(k) value(v(i)) * value(l(k)) - (value(P(k))^2 value(Q(k))^2); end end fprintf(最大松弛间隙: %.6e\n, max(gap));如果最大松弛间隙在1e-5量级以下可以放心把结果当精确解用。如果间隙明显偏大说明SOCP松弛不紧需要额外处理这是第5章的坑之一。4. 求解器怎么选组合决定MISOCP的命运4.1 常见求解器在SOCP问题上的表现对比Yalmip只是一个建模语言真正求解靠后端的求解器。很多人以为装了Yalmip就完事了其实还需要一个能解二阶锥的求解器。我用过的组合大概有几种求解器是否支持整数SOCP学术License实测表现Gurobi支持提供最快推荐首选Cplex支持提供快和Gurobi接近Mosek支持提供凸优化极强稳健SCIP支持开源可用速度略慢Sedumi/SDPT3不支持整数开源只能解连续SOCPfmincon不支持内置别用容易陷局部最优对33节点这种规模的配电网重构Gurobi通常能在几十毫秒到一两秒内给出MISOCP的最优解。Cplex也差不多。如果模型规模到几百个节点Gurobi的并行加速优势会更明显。这里特别说明一下Yalmip Gurobi 13.0.3这个组合我特意在新版本上重新测过之前担心接口有变动实际跑下来完全兼容没有出现低级报错。如果只是做连续SOCP即固定拓扑后求最优潮流Sedumi这类开源内点法也能用速度可以接受。但一旦模型里有二进制开关变量Sedumi直接歇菜Yalmip会调用内置的bnb外循环那个速度在30个节点以上就非常感人。所以我的结论很简单别在求解器上省时间该上Gurobi就上Gurobi性能差距不是一级两级。4.2 Gurobi参数设置与性能调优实测Gurobi装上以后不是默认参数就能跑得最好我调了这么多次有几个参数值得单独设置ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; % MIP间隙设为0.01%避免过早停止 ops.gurobi.TimeLimit 300; % 防止极端情况卡死 ops.gurobi.NumericFocus 1; % 若遇到数值问题可以调整到更高 ops.gurobi.Presolve 1; % 默认开启即可MIPGap这个参数很关键。Gurobi默认的MIPGap是0.0001甚至更宽松但在配电网重构问题里如果目标函数中网损和失电惩罚量级差太大默认间隙可能导致提前停止在一个不是最优的方案上。我通常手动设成0.0001让它在33节点这种小规模问题上直接证明最优性。另一个细节Gurobi读取Yalmip传递过来的MISOCP模型时会把二阶锥约束自动转换为它内部的表达形式。你不需要在Gurobi里手动声明“这是一个二阶锥”Yalmip会处理好。但如果你的模型里既有二阶锥又有大M约束Gurobi的presolve阶段可能会因为大M取值不合理把问题变坏。所以还是那句话大M够用就行别贪大。实测下来33节点系统的故障重构模型大约有60~80个二进制变量、几百个连续变量和上千条约束Gurobi求解时间通常在0.2秒到2秒之间。如果遇到某个故障场景特别难解比如故障点靠近馈线末端可以考虑增加一个故障前可行解作为MIP start可以明显加速。5. 我实际跑通这套代码踩过的五个坑5.1 数值标幺与大M取值不当导致的“伪不可行”第一次跑通模型的时候我直接用有名值结果Gurobi报模型不可行。排查半天发现不是约束写错了而是数值问题——线路电阻是0.4欧姆量级潮流功率是千瓦到兆瓦量级电压是10kV量级数字差了好几个数量级求解器内部的容差容忍不了这种尺度差异。解决办法很简单全部转换成标幺值。基准功率取10MVA基准电压取12.66kV换算后所有支路的r和x都在0.0001到0.1这个区间节点电压在0.9到1.1之间功率在0.01到10之间。这个尺度下大M取值也有了合理参考电压相关的M取10功率相关的M取100就足够。另外一个大M取值特别大的典型后果是LP松弛非常弱分枝定界的下界差求解时间暴涨。我见过有人把M直接写成1e6结果同样的33节点案例求解时间从0.5秒涨到十几秒。所以别偷懒按物理量的合理范围去估值才是正确姿势。5.2 辐射状约束不能只看“支路数N-1”这个坑我印象太深了。早期版本只写了sum(z) nbus - 1解出来一个网络闭合支路数是对的但后来画图发现有一个孤岛三个节点自己连成一棵树跟根节点完全脱离。孤岛内部的支路数加上主网的支路数总数恰好是N-1但从供电角度看完全是废的。补上虚拟潮流约束之后这个问题才彻底解决。我做了一个小验证随机生成几千个故障场景逐一检查求解结果闭合支路数全部等于N-1且每个节点都能通过闭合支路到达根节点没有再出现孤岛。所以强烈建议任何配电网重构模型都必须把连通性校验做进模型里而不是靠事后检查。5.3 二阶锥松弛不紧时的处理思路SOCP松弛在大部分情况下是紧的但也不是铁律。我试过一个重载场景约束条件绷得很紧最后求出来的解某个关键支路的松弛间隙达到1e-2量级说明锥约束被“撑开了”。遇到这种情况一个有效手段是迭代收紧先正常求一次SOCP然后把锥约束右侧加上很小的惩罚项重新求解往复几次松弛间隙会逐步压缩到可接受范围。更简单的做法是检查是不是电压约束过紧或负荷模型太极限很多松弛不紧的案例本质是电压越限边界上存在互斥趋势调整一下电压上下限就能恢复。我个人的工程习惯是求解完成后统一计算每个闭合支路的max gap如果超过1e-4就把对应支路的负荷或拓扑单拎出来分析看是数据问题还是模型问题。这比盲目换算法靠谱很多。5.4 开关操作次数的建模陷阱开关操作次数最开始我写的是abs(z - z0)心想这多直观。放进Yalmip后它不仅引入了额外的二进制变量来分解绝对值还让模型多了不少Big-M约束求解变慢。后来改成(z - z0)^2因为z和z0都是0/1变量平方的展开结果正好等于绝对值。Yalmip识别这种形式之后不会引入额外整数变量模型更紧凑求解速度能提升30%以上。当然还有更精细的写法比如把开关操作次数放在一个单独的目标权重系数里用scalar化调参。33节点案例里我没做太复杂直接固定权重系数0.1效果已经很好。如果是大型网络可以考虑给每段开关定义操作代价用z z0 - 2*z*z0这类双线性形式但这样又会引入非凸通常改成线性约束加辅助变量的形式更稳妥。5.5 从33节点扩展到更大网络的加速经验在33节点上跑通之后我尝试把同一套模型搬到更大的118节点系统发现求解时间会明显上升。主要瓶颈在二进制变量的数量和虚拟潮流约束的规模。有几个操作实测有效第一把MIPGap从0.0001放宽到0.001求解时间可能从几十秒降到几秒而目标值差距只有0.1%以内工程上完全能接受。第二给Gurobi传递一个启发式初始解比如直接把故障后的“全联络开关闭合”状态作为start能显著缩小搜索空间。第三尽量把负荷数据整理成稀疏矩阵的形式减少遍历次数。对于几百个节点的大型配电网还有一种思路是先用连续SOCP松弛去掉整数变量快速得到拓扑方向的参考再用邻域搜索恢复整数解。但这个属于进阶玩法了基础的MISOCP框架在100节点以内已经足够高效。我自己的体会是Matlab Yalmip Gurobi这套组合最大的价值在于它把“从问题到代码”的距离压到了最短。你不用亲手写凸优化求解算法也不用费劲把模型转成某种求解器专用的算子格式只要把物理约束用数学表达式一列Yalmip就能帮你在后端生成并求解。对于配电网故障重构这类需要频繁调整目标函数和约束的实战问题开发效率比直接用C调Gurobi高了一个数量级。最后再提醒一点不管用什么求解器跑完一定要回验拓扑和潮流结果。配电网重构不是拿到一个优化目标值就算完事这张网络是要拿去给真实用户供电的每一步都要对物理世界负责。
配电网故障重构实战：Yalmip+二阶锥规划建模指南
凌晨两点调度电话响起来10kV馈线跳闸重合失败故障点隔离后下游还有一片用户黑着。这时候最要紧的不是先查故障原因而是尽快给出一个“哪些开关合、哪些开关断”的供电恢复方案。这个决策背后就是配电网故障重构。我这两年复现和落地这类问题最顺手的组合是Matlab Yalmip搭模型把重构问题写成二阶锥规划SOCP来解稳定、快、还不容易掉进局部最优。这篇就把这套实现从头到尾掰开揉碎讲一遍问题怎么建模、二阶锥松弛为什么是首选、Yalmip代码骨架怎么搭、求解器怎么选以及我实际跑通之后踩过的几个坑。1. 故障重构到底在优化什么不只是“倒闸”1.1 问题边界故障隔离后的拓扑重构配电网和输电网最大的区别在于配电网基本都是辐射状结构闭环设计、开环运行。正常状态下网络是一棵树故障发生后保护装置把故障支路两侧的开关断开故障点被隔离但原本由这条支路往下送电的所有负荷会全部失电。这时候调度员能做的就是通过操作联络开关和分段开关把失电区域重新挂到别的馈线或者同一个变电站的其他出线上。故障重构问题的输入边界很明确给定一个故障前拓扑已知故障支路已经被强制隔离这一条支路在重构过程中不可恢复你需要找到一组新的开关状态使得整个网络的潮流满足运行约束同时让网损、开关操作次数、未恢复负荷量这些指标尽可能好。它本质上是一个带二进制变量的大规模混合整数优化问题——开关只有闭合和断开两种状态一旦加了潮流约束就成了混合整数非线性规划。但很多人第一次接触时容易把问题想窄了觉得“故障重构把联络开关合上就行”。实际上同时要考虑三件事第一合上联络开关会不会造成环网配电网闭环运行会导致保护配合混乱必须始终保持辐射状第二恢复后的电压是否越限远端节点很容易电压偏低第三线路容量能不能撑住负荷转移很多时候初始方案能恢复供电但某些支路过载照样不可行。1.2 目标函数与约束的完整画像一个完整的重构模型目标函数通常长这样min 网损 λ1 × 开关操作次数 λ2 × 未恢复负荷量网损所有闭合支路的电流平方乘以电阻之和公式上就是 Σ r_ij * L_ij。开关操作次数反映重构动作的代价。少动开关意味着倒闸快、风险低所以一般用重构后的开关状态与故障前状态做差的绝对值之和。未恢复负荷量如果系统容量足够这个量应该是0但如果发生严重故障或N-1不满足可能需要切除一部分负荷就要通过大权重强制优先恢复重要负荷。约束条件更琐碎但每一条都有实际物理含义约束类型表达式物理含义节点功率平衡Σ注入功率 - Σ流出功率 接入负荷每个节点的有功/无功必须守恒支路电压降支路两端电压满足欧姆定律线路压降不能随意设定二阶锥潮流约束v_i × L_ij ≥ P_ij² Q_ij²关联电压、电流和功率之间的非线性关系电压上下限V_min² ≤ v_i ≤ V_max²保证供电质量辐射状拓扑约束闭合支路数 节点数 - 1且全网连通避免环网与孤岛开关状态约束故障支路固定为0故障隔离不可恢复把这些东西全部写进模型才是一个可求解的完整问题。注意很多初学者写的模型只有潮流约束和开关变量忘了辐射状约束结果解出来是个环网还得靠人工后处理。我后面在Yalmip里实现的辐射状约束会结合虚拟潮流和支路数双重校验目前用下来比单独加一个“闭合支路数N-1”可靠得多。2. DistFlow的非凸困境与二阶锥松弛的价值2.1 DistFlow方程改写把潮流约束“翻译”成变量关系配电网潮流计算里用得最多的是DistFlow模型支路潮流模型。它的核心方程写出来并不复杂对每条支路 (i, j)功率平衡流入节点j的功率等于节点j的负荷加上从j流向所有子支路的功率和损耗电压降节点j的电压平方等于节点i的电压平方减去线路上压降带来的影响电流定义支路电流平方等于支路有功平方和无功平方之和除以节点电压平方如果直接把这些方程扔给求解器麻烦出在第三个方程它是电流、电压平方和功率之间的非线性隐式关系数学上是一个非凸的等式约束。非凸约束在优化里意味着多解、初值敏感、局部最优传统做法用遗传算法、粒子群这一类启发式算法去搜搜一次几十秒运气不好还搜不出可行解。有一个经典的改写技巧定义新变量把原来的物理量平方化。v_i U_i² l_ij I_ij²代入以后电压降方程从带平方根的复杂形式变成了关于 v、l、P、Q 的线性方程。这是DistFlow模型最漂亮的地方——把大部分非线性整理成了好处理的形态。但电流定义方程改写以后变成v_i × l_ij P_ij² Q_ij²左边是电压平方乘以电流平方右边是功率平方和。问题来了左边两个变量相乘这在优化里是双线性项依然是非凸的。这一步卡住了很多人。2.2 从双线性等式到二阶锥一步松弛的代价与边界关键转折点在于把上面的等式约束直接松成一个不等式。v_i × l_ij ≥ P_ij² Q_ij²这个不等式的几何意义是把一个抛物面非凸等式松弛成了一个二阶锥空间。二阶锥规划是凸优化里非常成熟的一类问题配合整数变量以后是MISOCP混合整数二阶锥规划商业求解器对付这种问题非常拿手。也许你会问松弛不等于原问题求出来的解能信吗这就要提到“精确松弛”的概念。在辐射状配电网中只要网络拓扑固定、负荷条件合理SOCP松弛在很多情况下是紧的也就是说松弛后的最优解恰好满足那个等号成立。你可以把它类比成“用一根直线去逼近一条曲线在绝大多数位置都贴得很紧”。当然不是所有情况都紧所以有经验的工程师在求解完成后都会主动检查松弛间隙不检查等于白做。这个坑我在第5章详细说。Yalmip里写这个锥约束并不需要手动推导出标准锥形式直接用norm写法即可% 标准二阶锥约束|| [2P; 2Q; v_i - l] || v_i l Constraints [Constraints, norm([2*P(k); 2*Q(k); v(fb(k))-l(k)], 2) v(fb(k)) l(k)];这一行的数学本质其实是在说“它们的取值被控制在一个圆锥曲面之内”而不要求严格落在这个曲面上。在配电网重构这个场景里只要求解结果把锥约束算到等号附近说明模型写对了。3. Yalmip建模实战一个33节点系统的完整代码骨架3.1 数据组织接线矩阵、支路参数与故障场景推荐用IEEE 33节点算例做测试平台它是配电网重构、故障恢复领域公认的标准算例。数据组织方式直接影响后面代码复杂度我习惯用四个数组fb每条支路的首端节点编号tb每条支路的末端节点编号r每条支路的电阻标幺值x每条支路的电抗标幺值另外还需要一个向量z0记录每条支路在故障前的闭合状态1表示闭合0表示断开。正常运行的33节点系统32条分段开关支路是闭合的5条联络开关支路是断开的。故障场景就用一个故障支路编号来模拟比如支路5-6发生永久性故障那么这条支路的z强制为0且在优化过程中不可操作。这里给一个数据格式的示例代码片段实际加载时直接从数据库或Excel读即可% 支路数据示例fb, tb, r(标幺), x(标幺), z0 branch_data [ 1, 2, 0.00922, 0.00470, 1; 2, 3, 0.00493, 0.00251, 1; % ... 中间略 ... 8, 21, 0.03420, 0.01690, 0; % 联络开关正常断开 ]; fb branch_data(:,1); tb branch_data(:,2); r branch_data(:,3); x branch_data(:,4); z0 branch_data(:,5);一定要全部转成标幺值再进模型。我见过不少直接把欧姆、安培、瓦特混在一起丢进Yalmip的结果要么求解严重病态要么收敛到离谱的答案。用标幺值之后节点电压在1.0附近支路电流也是0到几之间Gurobi和Cplex处理这种数值尺度最舒服。3.2 决策变量声明与约束块逐段拆解我们需要的决策变量包括z binvar(nbranch, 1); % 每条支路是否闭合 v sdpvar(nbus, 1); % 节点电压平方 l sdpvar(nbranch, 1); % 支路电流平方 P sdpvar(nbranch, 1); % 支路有功 Q sdpvar(nbranch, 1); % 支路无功 f sdpvar(nbranch, 1); % 虚拟潮流用于辐射状约束接着按约束块依次写。第一条是故障支路强制断开Constraints [Constraints, z(fault_branch) 0];第二条是节点功率平衡。对所有非电源节点所有流入该节点的支路功率减去流出该节点的支路功率等于该节点接入的负荷。负荷可以写成可切形式用一个二进制变量alpha表示负荷是否接入alpha binvar(nbus, 1); alpha(1) 1; % 电源节点始终接入 for n 2:nbus inflow_branches find(tb n); outflow_branches find(fb n); Constraints [Constraints, ... sum(P(inflow_branches)) - sum(P(outflow_branches)) alpha(n) * P_load(n)]; Constraints [Constraints, ... sum(Q(inflow_branches)) - sum(Q(outflow_branches)) alpha(n) * Q_load(n)]; end这里P_load(n)和Q_load(n)是节点n的有功和无功负荷标幺值。alpha为0意味着该节点负荷被切除功率平衡只要求注入等于流出相当于该节点变成一个纯转供节点。第三条是电压降约束。闭合支路必须满足DistFlow的电压降方程断开支路必须把这个约束松弛掉。这里用大M法处理for k 1:nbranch i fb(k); j tb(k); % 闭合时满足电压降方程 Constraints [Constraints, ... v(j) - v(i) 2*(r(k)*P(k) x(k)*Q(k)) - (r(k)^2 x(k)^2) * l(k) M_big * (1 - z(k))]; Constraints [Constraints, ... v(j) - v(i) 2*(r(k)*P(k) x(k)*Q(k)) - (r(k)^2 x(k)^2) * l(k) -M_big * (1 - z(k))]; % 二阶锥约束只对闭合支路有意义断开时功率为0锥自然满足 Constraints [Constraints, ... norm([2*P(k); 2*Q(k); v(i) - l(k)], 2) v(i) l(k)]; % 断开支路时支路功率和电流平方强制为0 Constraints [Constraints, -M_power * z(k) P(k) M_power * z(k)]; Constraints [Constraints, -M_power * z(k) Q(k) M_power * z(k)]; Constraints [Constraints, 0 l(k) M_power * z(k)]; end注意M_big和M_power不能随便取。M_big取1到10就够因为标幺电压本身在1附近压降绝对值远小于1M_power取系统最大视在功率的两倍以上即可但也不要大到1e6那种程度否则会拖慢MISOCP求解速度。电压上下限、变电站电压给定这些约束可以直接写V_min 0.95^2; V_max 1.05^2; Constraints [Constraints, V_min v V_max]; Constraints [Constraints, v(1) 1.0]; % 根节点电压恒定3.3 辐射状约束虚拟潮流避免环网辐射状约束是新手翻车率最高的地方。单独写一个“闭合支路数等于节点数减一”是必要条件但不是充分条件——一个“树孤岛”的网络也满足这个数量关系但孤岛内的节点根本没接上电源。所以必须在支路数约束之外再加连通性约束。我采用虚拟潮流方法让每个非电源节点向根节点发送1个单位的虚拟流只有闭合支路能传输虚拟流。数学上写成% 闭合支路数约束 Constraints [Constraints, sum(z) nbus - 1]; % 虚拟潮流上下界只有闭合支路能传流 Constraints [Constraints, -nbus * z f nbus * z]; % 每个非根节点发送1单位虚拟流 for n 2:nbus inflow_branches find(tb n); outflow_branches find(fb n); Constraints [Constraints, ... sum(f(inflow_branches)) - sum(f(outflow_branches)) 1]; end如果某个非根节点和根节点之间没有闭合路径它的虚拟流平衡方程就不可能满足整个问题不可行。反过来如果闭合支路数正好是N-1且每个节点都能向根节点传流那这个网络一定是一棵连通的生成树。这个写法的好处是不引入额外的大量整数变量只增加了一个连续向量f和一组线性约束对求解非常友好。3.4 目标函数与求解器设置目标函数把三部分加权组合起来。网损在标幺制下数值很小开关操作次数和未恢复负荷量要大得多所以权重设置上我建议手动归一化而不是拍脑袋给权重% 网损 Loss sum(r .* l); % 开关操作次数二进制变量直接平方等价于绝对值因为0/1的平方还是0/1 SwitchOps sum((z - z0).^2); % 未恢复负荷量 Unsupplied sum((1 - alpha) .* (P_load Q_load)); % 权重设置 lambda1 0.1; % 开关操作权重按基准值归一化 lambda2 100; % 失电惩罚必须远大于网损优先保证负荷恢复 Objective Loss lambda1 * SwitchOps lambda2 * Unsupplied; % 求解 ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; result optimize(Constraints, Objective, ops);求解结束后至少做三件事验证结果可靠性% 1. 检查求解状态 if result.problem 0 disp(求解成功); else disp(求解失败/不可行); end % 2. 提取重构后的拓扑 z_sol value(z); fprintf(闭合支路数: %d\n, sum(z_sol)); % 3. 检查二阶锥松弛间隙 % 理想情况下 v_i * l_ij 严格等于 P_ij^2 Q_ij^2 gap zeros(nbranch, 1); for k 1:nbranch if z_sol(k) 0.5 i fb(k); gap(k) value(v(i)) * value(l(k)) - (value(P(k))^2 value(Q(k))^2); end end fprintf(最大松弛间隙: %.6e\n, max(gap));如果最大松弛间隙在1e-5量级以下可以放心把结果当精确解用。如果间隙明显偏大说明SOCP松弛不紧需要额外处理这是第5章的坑之一。4. 求解器怎么选组合决定MISOCP的命运4.1 常见求解器在SOCP问题上的表现对比Yalmip只是一个建模语言真正求解靠后端的求解器。很多人以为装了Yalmip就完事了其实还需要一个能解二阶锥的求解器。我用过的组合大概有几种求解器是否支持整数SOCP学术License实测表现Gurobi支持提供最快推荐首选Cplex支持提供快和Gurobi接近Mosek支持提供凸优化极强稳健SCIP支持开源可用速度略慢Sedumi/SDPT3不支持整数开源只能解连续SOCPfmincon不支持内置别用容易陷局部最优对33节点这种规模的配电网重构Gurobi通常能在几十毫秒到一两秒内给出MISOCP的最优解。Cplex也差不多。如果模型规模到几百个节点Gurobi的并行加速优势会更明显。这里特别说明一下Yalmip Gurobi 13.0.3这个组合我特意在新版本上重新测过之前担心接口有变动实际跑下来完全兼容没有出现低级报错。如果只是做连续SOCP即固定拓扑后求最优潮流Sedumi这类开源内点法也能用速度可以接受。但一旦模型里有二进制开关变量Sedumi直接歇菜Yalmip会调用内置的bnb外循环那个速度在30个节点以上就非常感人。所以我的结论很简单别在求解器上省时间该上Gurobi就上Gurobi性能差距不是一级两级。4.2 Gurobi参数设置与性能调优实测Gurobi装上以后不是默认参数就能跑得最好我调了这么多次有几个参数值得单独设置ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 0.0001; % MIP间隙设为0.01%避免过早停止 ops.gurobi.TimeLimit 300; % 防止极端情况卡死 ops.gurobi.NumericFocus 1; % 若遇到数值问题可以调整到更高 ops.gurobi.Presolve 1; % 默认开启即可MIPGap这个参数很关键。Gurobi默认的MIPGap是0.0001甚至更宽松但在配电网重构问题里如果目标函数中网损和失电惩罚量级差太大默认间隙可能导致提前停止在一个不是最优的方案上。我通常手动设成0.0001让它在33节点这种小规模问题上直接证明最优性。另一个细节Gurobi读取Yalmip传递过来的MISOCP模型时会把二阶锥约束自动转换为它内部的表达形式。你不需要在Gurobi里手动声明“这是一个二阶锥”Yalmip会处理好。但如果你的模型里既有二阶锥又有大M约束Gurobi的presolve阶段可能会因为大M取值不合理把问题变坏。所以还是那句话大M够用就行别贪大。实测下来33节点系统的故障重构模型大约有60~80个二进制变量、几百个连续变量和上千条约束Gurobi求解时间通常在0.2秒到2秒之间。如果遇到某个故障场景特别难解比如故障点靠近馈线末端可以考虑增加一个故障前可行解作为MIP start可以明显加速。5. 我实际跑通这套代码踩过的五个坑5.1 数值标幺与大M取值不当导致的“伪不可行”第一次跑通模型的时候我直接用有名值结果Gurobi报模型不可行。排查半天发现不是约束写错了而是数值问题——线路电阻是0.4欧姆量级潮流功率是千瓦到兆瓦量级电压是10kV量级数字差了好几个数量级求解器内部的容差容忍不了这种尺度差异。解决办法很简单全部转换成标幺值。基准功率取10MVA基准电压取12.66kV换算后所有支路的r和x都在0.0001到0.1这个区间节点电压在0.9到1.1之间功率在0.01到10之间。这个尺度下大M取值也有了合理参考电压相关的M取10功率相关的M取100就足够。另外一个大M取值特别大的典型后果是LP松弛非常弱分枝定界的下界差求解时间暴涨。我见过有人把M直接写成1e6结果同样的33节点案例求解时间从0.5秒涨到十几秒。所以别偷懒按物理量的合理范围去估值才是正确姿势。5.2 辐射状约束不能只看“支路数N-1”这个坑我印象太深了。早期版本只写了sum(z) nbus - 1解出来一个网络闭合支路数是对的但后来画图发现有一个孤岛三个节点自己连成一棵树跟根节点完全脱离。孤岛内部的支路数加上主网的支路数总数恰好是N-1但从供电角度看完全是废的。补上虚拟潮流约束之后这个问题才彻底解决。我做了一个小验证随机生成几千个故障场景逐一检查求解结果闭合支路数全部等于N-1且每个节点都能通过闭合支路到达根节点没有再出现孤岛。所以强烈建议任何配电网重构模型都必须把连通性校验做进模型里而不是靠事后检查。5.3 二阶锥松弛不紧时的处理思路SOCP松弛在大部分情况下是紧的但也不是铁律。我试过一个重载场景约束条件绷得很紧最后求出来的解某个关键支路的松弛间隙达到1e-2量级说明锥约束被“撑开了”。遇到这种情况一个有效手段是迭代收紧先正常求一次SOCP然后把锥约束右侧加上很小的惩罚项重新求解往复几次松弛间隙会逐步压缩到可接受范围。更简单的做法是检查是不是电压约束过紧或负荷模型太极限很多松弛不紧的案例本质是电压越限边界上存在互斥趋势调整一下电压上下限就能恢复。我个人的工程习惯是求解完成后统一计算每个闭合支路的max gap如果超过1e-4就把对应支路的负荷或拓扑单拎出来分析看是数据问题还是模型问题。这比盲目换算法靠谱很多。5.4 开关操作次数的建模陷阱开关操作次数最开始我写的是abs(z - z0)心想这多直观。放进Yalmip后它不仅引入了额外的二进制变量来分解绝对值还让模型多了不少Big-M约束求解变慢。后来改成(z - z0)^2因为z和z0都是0/1变量平方的展开结果正好等于绝对值。Yalmip识别这种形式之后不会引入额外整数变量模型更紧凑求解速度能提升30%以上。当然还有更精细的写法比如把开关操作次数放在一个单独的目标权重系数里用scalar化调参。33节点案例里我没做太复杂直接固定权重系数0.1效果已经很好。如果是大型网络可以考虑给每段开关定义操作代价用z z0 - 2*z*z0这类双线性形式但这样又会引入非凸通常改成线性约束加辅助变量的形式更稳妥。5.5 从33节点扩展到更大网络的加速经验在33节点上跑通之后我尝试把同一套模型搬到更大的118节点系统发现求解时间会明显上升。主要瓶颈在二进制变量的数量和虚拟潮流约束的规模。有几个操作实测有效第一把MIPGap从0.0001放宽到0.001求解时间可能从几十秒降到几秒而目标值差距只有0.1%以内工程上完全能接受。第二给Gurobi传递一个启发式初始解比如直接把故障后的“全联络开关闭合”状态作为start能显著缩小搜索空间。第三尽量把负荷数据整理成稀疏矩阵的形式减少遍历次数。对于几百个节点的大型配电网还有一种思路是先用连续SOCP松弛去掉整数变量快速得到拓扑方向的参考再用邻域搜索恢复整数解。但这个属于进阶玩法了基础的MISOCP框架在100节点以内已经足够高效。我自己的体会是Matlab Yalmip Gurobi这套组合最大的价值在于它把“从问题到代码”的距离压到了最短。你不用亲手写凸优化求解算法也不用费劲把模型转成某种求解器专用的算子格式只要把物理约束用数学表达式一列Yalmip就能帮你在后端生成并求解。对于配电网故障重构这类需要频繁调整目标函数和约束的实战问题开发效率比直接用C调Gurobi高了一个数量级。最后再提醒一点不管用什么求解器跑完一定要回验拓扑和潮流结果。配电网重构不是拿到一个优化目标值就算完事这张网络是要拿去给真实用户供电的每一步都要对物理世界负责。
找出最大离群值)
