配电网最优潮流的二阶锥松弛原理与MATLAB实现

关于配电网最优潮流这四个字凡是做过配电网研究的同学基本都绕不开。它做的事情可以概括为一句话在满足潮流方程、电压限值、支路容量这些物理约束的前提下找到一组控制量让网损、成本或者电压偏差这类运行指标最优。听起来就是标准优化问题可真的把交流潮流方程摆进去问题立刻变得麻烦——方程里有电压乘积项、功率二次项整个约束空间是非凸的传统算法很容易陷在局部最优里出不来算十次可能得到十个不同的“最优解”。二阶锥松弛SOCP relaxation这几年几乎成了解决配电网OPF的默认首选方法。它的思路很直接把原来潮流约束里唯一“惹事”的非凸等式替换成一个良结构的凸锥约束把整个问题转化成二阶锥规划。这样的模型可以被成熟求解器高效求解而且在辐射状配电网中松弛后的最优解往往恰好满足原等式约束也就是“松而不驰”结果完全可以作为真实最优解使用。这篇文章就以IEEE 33节点标准配电网算例为对象用Matlab YALMIP MOSEK做一套完整的二阶锥松弛最优潮流实现。代码环节会给出可以直接改着用的核心脚本并把变量定义、约束构建、结果校验都讲透。适合两类人看一是刚进配电网优化方向、想快速搭出第一个能跑模型的同学二是在工程里需要做配电网运行评估、想少踩坑的工程师。1. 为什么配电网最优潮流这么难算1.1 最优潮流的“难”不只在于规模输电网的最优潮流问题已经研究了几十年但配电网的OPF往往比输电网更麻烦。原因主要有三个第一配电网结构是辐射状的支路数等于节点数减一没有环网冗余任何一个支路潮流都承担着明确的供电任务对潮流方程的计算精度要求更高。第二配电网的电阻电抗比R/X通常比较大这意味着有功功率和无功功率的耦合比输电网强得多不能用输电网里常用的快速解耦假设来简化计算。第三近年来分布式电源大量接入光伏、风电、储能让配电网从“单向潮流”变成“双向潮流”电压分布不再是简单的从变电站往末端递减优化问题变得更加复杂。在这些特征下最优潮流的标准数学形式可以写成目标函数比如最小化网损 ( \sum_{k1}^{N_b} r_k l_k ) 或者最小化发电成本约束条件节点有功/无功平衡方程、支路潮流方程、节点电压上下限、支路容量限值、控制变量限值等。问题在于其中的支路潮流方程是二次等式约束。二次等式约束构成的可行域是一个弯曲的非凸曲面在这个曲面上做优化就像在连绵起伏的山脊上找最低点。经典内点法能不能找到全局最优很大程度上取决于初始点选在哪里启发式算法虽然能“撒网式”搜索但收敛速度慢而且同样不能保证全局最优。实际项目中如果你用通用非线性规划求解器去直接求解配电网OPF很容易遇到求解失败或者收敛到明显不合理的局部解的情况。1.2 传统方法为什么不顶用我经常看到有人一上来就用智能算法比如粒子群、遗传算法去解配电网最优潮流理由是“这类方法不用求导实现简单”。这话没错但它们有两个硬伤一是每次计算都要反复调用潮流计算作为适应度评估算一个中等规模33节点系统还能忍到了几百上千节点的实际馈线就非常吃力二是这些算法没有全局最优的理论保证参数稍微调不好结果可能比常规运行方式还差。另一种思路是直接非线性规划。把OPF原问题丢给内点法求解器能处理一定规模但还是那个问题非凸性决定了它只能保证局部最优。对于严格凸的目标函数和简单约束局部最优也许还凑合可一旦加上分布式电源的容量约束、支路容量约束、电压硬约束多局部最优的情况非常普遍。所以真正靠谱的路线是把非凸问题“凸化”。凸化有两类主流手段一类是半定松弛SDP另一类就是本文要展开的二阶锥松弛SOCP。SDP对网络结构的适应范围更广但计算量大SOCP虽然对网络结构有辐射状的要求但求解规模更大、速度更快。而配电网恰恰绝大多数是辐射状网络这就让SOCP有了天然的用武之地。2. 二阶锥松弛的核心原理先找非凸点再精确放掉它2.1 从DistFlow模型说起要理解二阶锥松弛最自然的切入点是支路潮流方程Branch Flow Model也叫DistFlow。它专门面向辐射状配电网把每条支路的潮流和节点电压用一组递推方程刻画。定义变量如下( v_i )节点i的电压幅值平方( l_k )支路k的电流幅值平方( P_k )支路k首端流出的有功功率( Q_k )支路k首端流出的无功功率( r_k, x_k )支路k的电阻、电抗。对从节点i流向节点j的支路kDistFlow写为[ v_j v_i - 2(r_k P_k x_k Q_k) (r_k^2 x_k^2) l_k ]这是电压降落方程。节点功率平衡方程则写成流入减去流出的形式。这几个式子本身不是问题它们都是线性方程真正导致非凸的是下面的关系[ l_k \frac{P_k^2 Q_k^2}{v_i} ]这个约束来自电流、功率和电压的固有物理关系。因为分母带有变量 ( v_i )等式两边都出现二次项所以它作为一个等式约束是非凸的。整体上这就像在模型里埋下了一颗“非凸的钉子”拔不掉整个问题就没法用凸优化求解器处理。2.2 松弛把等式放宽成二阶锥二阶锥松弛的做法很简单甚至可以说是“偷懒但讲道理”。把上面的非凸等式约束放宽成不等式[ l_k \ge \frac{P_k^2 Q_k^2}{v_i} ]或者写成更常见的形式[ P_k^2 Q_k^2 \le v_i l_k ]这个不等式可以进一步等价改写为标准二阶锥形式[ \left| \begin{bmatrix} 2P_k \ 2Q_k \ v_i - l_k \end{bmatrix} \right|_2 \le v_i l_k ]为什么要写成最后这种形式因为绝大多数商业求解器和高性能开源求解器对标准二阶锥约束都有专门的内点算法识别和求解效率远高于处理普通二次不等式。YALMIP在建模时能自动识别不少形式的二阶锥约束但你在代码里显式写出这种标准形式求解器处理起来会更直接。物理直觉上讲这个松弛相当于允许支路电流“不小于”功率实际要求的最小值。原问题要求电流与功率严格一致放宽之后电流可以多损耗一点点。但注意如果你优化的目标是网损或者发电成本电流增加只会增加损耗目标函数自然会驱使最优解回到锥的边界上。这就是SOCP松弛在配电网里“精确”的原因——不是数学上所有情况都精确而是在实际常见场景下最优解会自发地满足等式关系。2.3 松弛为什么往往就是精确的很多初学的人会担心松弛之后的最优解如果不在锥边界结果还能用吗这就需要理解“精确松弛”的成立条件。从文献研究的结论来看对于辐射状配电网如果目标函数是关于支路电流的单调递增函数网损就是典型例子同时在最优解处节点电压不越上限那么SOCP松弛就是紧的松弛后的最优解还原出的潮流与交流潮流的真实解一致。工程上的绝大多数算例包括各种标准测试系统都满足这些条件。所以跑出来的结果可以直接当交流最优潮流的最优解用。为了排除意外实操中建议每次求解完都顺手算一下“松弛间隙”[ \text{gap}_k v_i l_k - (P_k^2 Q_k^2) ]如果所有支路的gap都接近0说明解在锥边界上可以放心使用。如果某条支路的gap明显大于0就需要谨慎处理具体排查方法我会在第4节展开。3. Matlab YALMIP实现33节点算例从建模到求解3.1 环境准备与算例选择代码运行前提是三样东西Matlab、YALMIP工具箱、一个能解二阶锥规划的求解器。YALMIP安装不复杂从官网下载后把整个目录加入Matlab路径在命令行执行yalmiptest看到一堆测试通过即可。求解器方面最省心的选择是MOSEK或者Gurobi两者学术许可免费解SOCP都非常高效。手头没有商业求解器的话免费的SCS和ECOS也支持SOCP在小规模算例上足够用就是精度和速度略逊一筹。算例数据我推荐Matpower自带的case33bw也就是IEEE 33节点配电网系统。它的基准数据是基准电压12.66kV基准容量10MVA总共32条支路网络是标准的辐射状结构包含一个变电站节点和若干负荷节点。这个系统规模不大不小正好适合做教学验证和算法测试。mpc loadcase(case33bw); baseMVA 10; bus mpc.bus; branch mpc.branch; n size(bus, 1); br size(branch, 1); Vbase 12.66; % kV Zbase Vbase^2 / baseMVA; % 单位欧姆注意这里最关键的一步是标幺化。Matpower里的数据通常是有名值欧姆、kW/kvar而优化模型里如果用有名值数值可能跨好几个数量级求解器很容易出现数值稳定性问题。把电阻电抗除以Zbase、把负荷除以baseMVA之后所有变量都在0到1附近波动求解会顺畅很多。这是我实际跑配电网优化时最深的体会之一配电网优化90%的数值问题都是没做好标幺化。3.2 用BFS构建辐射状拓扑关系DistFlow方程是沿着支路方向递推的所以你必须知道每条支路的父节点是谁、子节点是谁。很多Matpower数据的支路表恰好就是按“从变电站往外”的方向排列的直接用fb、tb可以省事。但我不建议依赖这个约定更稳妥的办法是用一次广度优先搜索BFS来重建父子关系。children cell(1, n); parent_branch zeros(1, n); visited false(1, n); queue 1; visited(1) true; while ~isempty(queue) i queue(1); queue(1) []; for k 1:br if branch(k,1) i ~visited(branch(k,2)) j branch(k,2); visited(j) true; parent_branch(j) k; children{i}(end1) k; queue(end1) j; end end end这段代码从根节点出发按照支路连接关系遍历整个网络每个节点的父支路编号和子支路集合就都建好了。为什么要费这个事因为一旦你把算例从33节点换到更复杂的系统支路数据的排列顺序不一定还那么“规整”BFS能保证方向信息永远是对的。程序里写错一个方向潮流方程就反了得到的解基本是乱的。3.3 核心约束代码变量定义与约束拼装在YALMIP里把每个优化变量定义为sdpvar即可。对33节点系统变量规模很小全用全精度符号类型都没问题。V sdpvar(n, 1, full); % 节点电压幅值平方 l sdpvar(br, 1, full); % 支路电流幅值平方 P sdpvar(br, 1, full); % 支路有功 Q sdpvar(br, 1, full); % 支路无功 Pg sdpvar(1, 1, full); % 变电站有功注入 Qg sdpvar(1, 1, full); % 变电站无功注入然后逐条拼约束。先是节点功率平衡。对节点i变电站注入减负荷等于“从父支路流入的功率”减去“流向所有子支路的功率”。Constraints []; for i 1:n lhsP Pg * (i 1) - bus(i,3) / baseMVA; lhsQ Qg * (i 1) - bus(i,4) / baseMVA; rhsP 0; rhsQ 0; if parent_branch(i) 0 kp parent_branch(i); rhsP rhsP P(kp) - R(kp) * l(kp); rhsQ rhsQ Q(kp) - X(kp) * l(kp); end for k children{i} rhsP rhsP - P(k); rhsQ rhsQ - Q(k); end Constraints [Constraints, lhsP rhsP]; %#okAGROW Constraints [Constraints, lhsQ rhsQ]; %#okAGROW end注意父支路流入节点i的功率是支路首端功率减去线损。有功线损是 ( R_k l_k )无功线损是 ( X_k l_k )。这个细节很多人会漏漏掉的后果就是功率平衡算出来总有偏差网损对不上。然后是DistFlow电压降方程。对每条支路k节点i是父端节点j是子端for k 1:br i branch(k,1); j branch(k,2); Constraints [Constraints, ... V(j) V(i) - 2*(R(k)*P(k) X(k)*Q(k)) ... (R(k)^2 X(k)^2) * l(k)]; end再是二阶锥约束。直接使用cone命令把约束写成标准二阶锥形式for k 1:br i branch(k,1); Constraints [Constraints, ... cone([2*P(k); 2*Q(k); V(i)-l(k)], V(i)l(k))]; end如果你更喜欢可读性强的写法也可以直接写P(k)^2 Q(k)^2 V(i) * l(k)YALMIP会尝试自动识别并归入二阶锥类型。但我个人用下来显式cone写法对不同求解器的兼容性更好尤其当你从MOSEK切到SCS时不会出现模型类型识别失败的问题。最后加上电压上下限、支路容量约束Vmin 0.95^2; Vmax 1.05^2; Constraints [Constraints, Vmin V Vmax]; Constraints [Constraints, 0 l 0.1]; objective sum(R .* l); ops sdpsettings(solver, mosek, verbose, 2); sol optimize(Constraints, objective, ops);目标函数这里取最小化有功网损 ( \sum R_k l_k )。网损是支路电流平方乘以电阻本身是线性的所以整体模型是纯二阶锥规划不是混合整数问题求解速度极快。33节点模型在MOSEK下一个求解循环大概只需要一两秒。3.4 结果后处理与松弛紧性验证求解完成之后第一件事不是看网损数值而是验证松弛是不是紧的。这一步永远不要省。Vopt value(V); lopt value(l); Popt value(P); Qopt value(Q); loss_mw value(objective) * baseMVA; gap max(Vopt(branch(k,1)) .* lopt - Popt.^2 - Qopt.^2);gap如果小于 (1\times10^{-5}) 这个量级说明所有支路都处在锥边界上松弛是精确的网损结果可以直接作为交流最优潮流的解来用。如果gap跑到 (1\times10^{-2}) 甚至更大就说明某条支路的电流没有贴近最小允许值此时你需要回到第4节排查原因。电压恢复也很简单节点电压幅值就是sqrt(Vopt)每条支路的电流幅值就是sqrt(lopt)。因为求解的是电压幅值平方拿到的全是正值不需要考虑相位的问题。4. 常见问题与排查这些坑我基本都踩过4.1 松弛间隙大结果还能不能用松弛间隙偏大的情况不多但一旦出现就很让人头疼。根据我的经验主要原因通常有三个一是目标函数不是关于电流单调递增的。比如你优化的是电压偏差目标函数可能会故意让某些支路电流变大来把电压抬到更“居中”的位置这时候松弛不紧是正常的。解决思路是换一个单调目标或者对松弛间隙项加一个极小权重的惩罚引导最优解回到锥边界。二是支路容量约束在起作用。当一条支路满载或接近满载时电流的上限被约束夹住了目标函数没能力再把它压到最小间隙就会出现。这种情况如果出现在个别关键支路上可以适当放宽容量限值或者把容量约束改成软约束来观察模型变化。三是网络本身处于重载状态或者网络不是纯辐射状结构。环网结构下SOCP松弛的理论基础会变弱这时更保险的选择是换用SDP松弛或者干脆回到非线性规划配合多初值策略。4.2 标幺制和数值稳定性90%的报错根源说一个花了很久才想明白的经验配电网SOCP-OPF跑不通绝大多数时候不是模型写错了而是数值尺度出了问题。配电网有名值的特点是电压上万伏但标幺后接近1负荷几十千瓦但标幺后是小数值支路电阻可能只有零点几欧姆。如果把这些数直接丢进优化模型约束矩阵不同行的数量级差出好几个数量级求解器内部的数值分解会非常不稳定轻则警告重则直接报Numerical issues。标幺化之后电压、潮流、网损都在合理范围内求解器几乎不会抱怨。我自己的习惯是在每个模型里都先定义基准值然后所有数据都除以基准值再参与建模最后展示结果时再乘回去。这套流程虽然多写几行代码但能省掉后面所有的定位时间。4.3 求解器不识别SOCP怎么让模型更“标准”有些朋友遇到的问题是代码逻辑看着没问题但求解器一直报“无法处理次方/平方项”之类的错误。这通常是因为YALMIP没能把你的约束自动归类为二阶锥约束尤其是用了形如P^2 Q^2 V*l这种乘法二次约束时。解决办法很简单用cone命令显式写标准形式不要让工具去猜。另一个小技巧是在定义变量之初就指定sdpvar(n, 1, full)不要留无符号变量避免YALMIP为了处理符号问题把模型变成更复杂的类型。4.4 分布式电源接入的几个建模坑接入分布式电源后最常犯的错是直接把光伏逆变器当成标准PV节点建模。配电网里逆变器通常没有无功调节能力或者无功限值很小强行PV节点会让问题变成非凸。更稳妥的做法是把分布式电源都处理成PQ节点。具体来说如有功出力为 (P_g)无功出力范围 ([-Q_{\max}, Q_{\max}])那么在功率平衡方程左边加上这些注入量同时给无功加一个盒约束即可PG sdpvar(ngen, 1); QG sdpvar(ngen, 1); Constraints [Constraints, PG_min PG PG_max]; Constraints [Constraints, -Qmax QG Qmax];每个DG在哪条母线接入就在对应母线的功率平衡方程里加上相应的 (P_g) 和 (Q_g)。这样整个模型依然是SOCP可以参考复现。5. 还能往哪个方向扩展SOCP-OPF不只是“算个网损”二阶锥松弛配电网最优潮流这套框架最大的好处是“建模规范求解可规模化”所以稍微改改目标函数和约束就能适配很多实际课题。想把网损换成运行成本只要在目标函数里加入分布式电源的发电成本项。成本函数如果是分段线性的可以引入辅助变量写成线性目标如果是二次函数且凸可以直接放进SOCP模型。想做配电网重构把开关状态设为0/1整数变量目标函数和潮流约束都不动问题就变成混合整数二阶锥规划MISOCP求解器同样能处理。想接入智能软开关SOP或者储能系统也是一样。储能调度模型里的荷电状态递推方程是线性的充放电功率约束用两个非负变量表示就能避免非凸的双向耦合。SOP的有功/无功变量约束本质上是多边形容量约束可以拆成一组线性或锥约束加入模型。如果再往上走三相不平衡配电网的SOCP模型考虑三相之间的耦合约束复杂度会提升一个量级但思路仍然是一致的。很多实际工程里的配电网电压无功优化、馈线负载均衡问题最终都会落到这个框架上。最后分享一个我自己的经验跑这类模型最忌讳上来就搞大系统。花半天时间把33节点系统彻底跑透包括BFS拓扑、约束方向、标幺制、松弛间隙验证这些细节点都弄清楚再换到任何实际系统都是顺水推舟的事。后面即使遇到问题你也知道先检查哪个环节——像“支路方向反了”和“未标幺化”这两类问题每个初用SOCP的人都会遇到至少一次但只要你把检查和验证习惯建立起来它们基本只能浪费你十分钟而不是一个下午。