二阶锥松弛这四个字对刚接触配电网最优潮流的人而言往往是一道绕不过去的坎。你搜“IEEE33节点 最优潮流 yalmip cplex”能跳出来的结果不是只有一段含糊的原理描述就是一份把非凸模型直接硬塞给求解器然后报错的代码。这篇文章我把完整链路写清楚从配电网最优潮流为什么非凸到DistFlow模型如何通过二阶锥松弛变成凸优化问题再到MATLABYALMIPCPLEX环境下如何在IEEE33节点系统上一步步建模、求解和验证。写代码和跑算例的原始目的本来是为了服务一篇含分布式电源接入的配电网优化研究后来精简成了一个可以反复复现的标准算例。搞懂这一套你再去看其他带DG、储能或者重构的扩展模型基本就不慌了。1. 为什么配电网最优潮流必须先“松弛”一下1.1 非凸性根源二次等式和NP-hard配电网最优潮流OPF这个问题本质上是在满足潮流方程、电压约束、支路容量约束的前提下寻找使网损或运行成本最小的运行点。听起来不复杂但你只要把潮流方程展开就会发现麻烦来了。交流潮流方程中节点注入功率与节点电压之间是二次关系比如支路潮流里会出现电压平方项、电压乘积项、三角函数项。即使使用极坐标形式也仍然逃不掉电压幅值的二次项和相角的三角函数项。这些等式约束叠加起来使得整个可行域是一个非凸集合。非凸意味着什么意味着两个可行潮流解的线性组合很大概率不再是一个可行解。这就像一个山谷里的两个低点之间有一道陡壁凸优化工具只能在“从左到右单调上升或下降”的地形里找极值碰到陡壁就直接宣告失败。更麻烦的是通用最优潮流问题在数学上已经被证明是NP-hard的也就是说不存在一种多项式时间的算法能对任意规模的网络都保证找到全局最优解。所以就有了一个很朴素的思路既然非凸等式太硬我能不能把这个等式约束放成一个不等式先求一个“放宽松”的问题再看这个问题的解是否满足原来的等式这就是松弛。二阶锥松弛就是其中工程上最成功的一种做法。1.2 二阶锥松弛的几何直觉二阶锥松弛的核心操作是把原本必须取等号的约束改成允许取不等号。听起来像偷工减料但妙就妙在在特定的目标函数和网络结构下这个“放宽”并不会改变最终的最优解。可以这么理解原问题要求一个点在某条曲线上移动比如要求 voltage * current power松弛后允许这个点在曲线外侧的整个凸区域内移动。如果目标函数是朝着“让这个点尽量贴回曲线”的方向用力比如最小化网损那么优化器在追求最小目标的过程中会把点“推”到曲线边界上。这时候松弛问题的解和原问题的解重合了。这个“推回去”的过程并不是靠运气保证的。对于辐射状配电网在目标函数关于电流平方单调递增等若干条件下理论上有严格的精确性证明。这也是为什么近十几年的配电网优化文献里二阶锥松弛几乎成了标配工具。但有一点必须提醒不是所有问题都能保证松弛精确。如果目标函数有奇怪的形状或者电压下限约束太紧解有可能落在锥内部而不是锥面上。所以跑出结果后必须做后验校验这个在后面的章节我会专门讲。2. DistFlow模型与SOCP变换把非凸OPF改写成标准形2.1 DistFlow方程回顾在讨论配电网最优潮流时几乎绕不开DistFlow方程。这个模型专门针对辐射状配电网设计把潮流方程从节点功率注入形式改写成了沿着馈线逐段递推的形式。考虑一条支路首端节点为i末端节点为j支路电阻为r、电抗为x。令P_ij和Q_ij为从i流向j的有功和无功功率V_i为节点i的电压幅值I_ij为流过支路ij的电流幅值。DistFlow方程可以写成节点j的有功平衡P_ij - r * I_ij² 节点j下游所有出线有功之和 节点j的有功负荷节点j的无功平衡Q_ij - x * I_ij² 节点j下游所有出线无功之和 节点j的无功负荷电压降落关系V_j² V_i² - 2(r P_ij x Q_ij) (r² x²) I_ij²这套方程里最棘手的是最后那个等式关系它把P、Q、V、I全部耦合在一起。尤其是支路潮流本身还有一个必须满足的物理关系P_ij² Q_ij² V_i² I_ij²这个等式约束就是整个问题非凸的根源。只要它存在YALMIP就会告诉您“nonconvex”或者“无法识别”CPLEX也直接不接单。2.2 变量替换与旋转锥形式既然二次等式难处理那就做变量替换。引入两个新变量U_i V_i²代表节点i的电压幅值平方L_ij I_ij²代表支路ij的电流幅值平方。代入上面的物理关系式原来的等式就变成了P_ij² Q_ij² U_i L_ij右边出现了两个变量的乘积还是非凸。这时二阶锥松弛登场把等号松弛为大于等于号P_ij² Q_ij² ≤ U_i L_ij这个不等式在数学上描述的不是普通区域而是一个“旋转二阶锥”。它和标准的二阶锥是等价的。怎么等价做一下代数变形已知 (U L)² - (U - L)² 4 U L所以如果 U L ≥ P² Q²那么必然有(2P)² (2Q)² (U - L)² ≤ (U L)²也就是norm([2P; 2Q; U - L]) ≤ U L这才是标准二阶锥形式。为什么叫“二阶锥”因为这条不等式的几何图像是一个锥体锥体上每个截面的半径都正比于高度约束函数是范数属于二阶锥。在YALMIP里这条约束可以直接写成cone([2*P(i,j); 2*Q(i,j); U(i) - L(i,j)], U(i) L(i,j))少数情况下也可以直接写 U(i)*L(i,j) P(i,j)^2 Q(i,j)^2YALMIP如果识别到这是旋转锥结构会自动转换。但我在实际使用中遇到过老版本YALMIP把它当成双线性非凸约束的情况所以更推荐显式使用cone稳定且语义清晰。2.3 松弛精确性与后验校验松弛做完之后很多人最关心的就是把等号变成了大于等于号求出来的解还是原问题的解吗答案是不一定但在多数配电网算例中是。判断方法需要回到约束本身。我们定义“松弛间隙”为gap U_i L_ij - (P_ij² Q_ij²)如果松弛精确gap应该等于0也就是说原等式约束重新成立。如果gap明显大于0说明最优解落在了锥内部松弛不精确这时SOCP问题的解虽然满足所有约束但它不是原潮流方程的解物理上不可行。在实际工程中只要目标函数是网损最小或者与电流平方正相关且电压上下限没有把可行域压得过于极端gap往往都很小。这也是IEEE33节点这类标准算例能够反复跑通的原因。但如果你加入了很多DG且DG出力大到把电压抬到上限附近电压下限约束也可能夺走最优性这时候就需要认真检查gap必要时对目标函数加一个很小的惩罚项比如10^-6倍的ΣL迫使解回到锥面上。3. IEEE33节点算例从原始数据到YALMIP建模3.1 IEEE33节点系统参数与标幺化IEEE33节点系统是配电网研究里的“标准练习册”它由一个12.66kV中压配电网组成共33个节点、32条支路、5条联络开关。联络开关在正常运行状态下是打开的所以网络呈辐射状。在使用任何求解器之前第一件事是统一量纲。我的经验是千万千万不要直接在有名值下建模。12.66kV的电压平方量级是1.6e8而支路电流平方量级可能只有1e-7两者混在一个约束里CPLEX会直接报数值问题或者给出让你怀疑人生的结果。标准做法是标幺化。本算例建议取Sbase 10 MVAVbase 12.66 kV则阻抗基准Zbase Vbase² / Sbase 16.0256 Ω。所有支路电阻、电抗除以Zbase得到标幺值所有节点负荷有功和无功分别除以Sbase得到标幺值。这样处理后节点电压在1.0附近支路电流平方在0.05到0.3之间数值很均匀求解器跑起来又快又稳。下面我直接把IEEE33的完整数据写在代码里。支路数据格式为“首端节点、末端节点、电阻欧姆、电抗欧姆”负荷数据格式为“节点编号、有功kW、无功kvar”。3.2 在YALMIP中声明变量和书写锥约束建模时我习惯把支路变量定义成N×N的完整矩阵而不是一维向量。例如P是33×33矩阵P(i,j)表示从节点i流向节点j的有功功率。这样做的好处是节点功率平衡方程写起来非常直观不会有“变量序号和支路序号对不上”的烦恼。对于IEEE33这种规模完整矩阵完全够用如果将来做上千节点的网络我再建议改成稀疏向量形式。核心约束分为三块第一块二阶锥约束。对每条支路(i,j)添加Constraints [Constraints, cone([2*P(i,j); 2*Q(i,j); U(i) - L(i,j)], U(i) L(i,j))];第二块电压降落方程。同样对每条支路Constraints [Constraints, U(j) U(i) - 2*(r*P(i,j) x*Q(i,j)) (r^2 x^2)*L(i,j)];第三块节点功率平衡。对每个节点j根节点除外找到它的父节点i然后约束流入功率等于下游出线功率与节点负荷之和。这里要注意辐射状网络中每个节点只有一个父节点这是DistFlow能正确书写的前提。如果网络里有联络线闭合那就不能再简单这样写。3.3 目标函数与运行约束配电网最优潮流最经典的目标函数是网损最小。网损原本需要写成 r*(P²Q²)/U但在二阶锥模型里这个表达式已经被变量替换消化掉了。由于U_i L_ij ≥ P_ij² Q_ij²网损在目标函数中可以安全地线性化为loss Σ r_ij * L_ij这个“线性化红利”是二阶锥松弛非常漂亮的地方。原问题目标里全是二次项松弛后目标函数居然变成了线性函数整个问题变成“线性目标 二阶锥约束”在凸优化里属于非常容易求解的结构。电压约束方面本算例设置根节点电压为1.0 pu其余节点电压允许范围0.9~1.05 pu写成Constraints [Constraints, 0.9^2 U 1.05^2];注意这里约束的直接是电压平方U而不是电压本身。支路电流上限则根据实际线路载流量设置示例代码中我给了Ilimit 0.5对应基准电流约456A下的0.5pu也就是约228A这个值比较宽松不会影响正常结果。4. 完整代码、结果验证和DG扩展场景4.1 可直接运行的完整MATLAB程序以下程序是完整的复制到MATLAB中确保YALMIP和CPLEX可用直接运行即可得到IEEE33节点网损最小化的SOCP结果。%% IEEE33节点配电网最优潮流 —— 二阶锥松弛 (YALMIP CPLEX) % 基准值: Sbase 10 MVA, Vbase 12.66 kV clear; clc; % 支路数据: [首端节点 末端节点 电阻(ohm) 电抗(ohm)] branch [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; 3 4 0.3660 0.1864; 4 5 0.3811 0.1941; 5 6 0.8190 0.7070; 6 7 0.1872 0.6188; 7 8 0.7114 0.2351; 8 9 1.0300 0.7400; 9 10 1.0440 0.7400; 10 11 0.1966 0.0650; 11 12 0.3744 0.1238; 12 13 1.4680 1.1550; 13 14 0.5416 0.7129; 14 15 0.5910 0.5260; 15 16 0.7463 0.5450; 16 17 1.2890 1.7210; 17 18 0.7320 0.5740; 2 19 0.1640 0.1565; 19 20 1.5042 1.3554; 20 21 0.4095 0.4784; 21 22 0.7089 0.9373; 3 23 0.4512 0.3083; 23 24 0.8980 0.7091; 24 25 0.8960 0.7011; 6 26 0.2030 0.1034; 26 27 0.2842 0.1447; 27 28 1.0590 0.9337; 28 29 0.8042 0.7006; 29 30 0.5075 0.2585; 30 31 0.9744 0.9630; 31 32 0.3105 0.3619; 32 33 0.3410 0.5302]; % 负荷数据: [节点编号 有功(kW) 无功(kvar)] load_data [ 1 0 0; 2 100 60; 3 90 40; 4 120 80; 5 60 30; 6 60 20; 7 200 100; 8 200 100; 9 60 20; 10 60 20; 11 45 30; 12 60 35; 13 60 35; 14 120 80; 15 60 10; 16 60 20; 17 60 20; 18 90 40; 19 90 40; 20 90 40; 21 90 40; 22 90 40; 23 90 50; 24 420 200; 25 420 200; 26 60 25; 27 60 25; 28 60 20; 29 120 70; 30 200 600; 31 150 70; 32 210 100; 33 60 40]; % 标幺化 Sbase 10e6; Vbase 12.66e3; Zbase Vbase^2 / Sbase; nb size(branch, 1); N 33; Br branch(:,3) / Zbase; Bx branch(:,4) / Zbase; fb branch(:,1); tb branch(:,2); Pd zeros(N,1); Qd zeros(N,1); for k 1:size(load_data,1) idx load_data(k,1); Pd(idx) load_data(k,2)*1e3 / Sbase; Qd(idx) load_data(k,3)*1e3 / Sbase; end % 优化变量 U sdpvar(N,1); % 节点电压幅值平方 (pu) P sdpvar(N,N,full); % 支路有功 (pu) Q sdpvar(N,N,full); % 支路无功 (pu) L sdpvar(N,N,full); % 支路电流幅值平方 (pu) % 约束集合 Constraints []; % 支路二阶锥约束 电压降落方程 for k 1:nb i fb(k); j tb(k); r Br(k); x Bx(k); Constraints [Constraints, cone([2*P(i,j); 2*Q(i,j); U(i) - L(i,j)], U(i) L(i,j))]; Constraints [Constraints, U(j) U(i) - 2*(r*P(i,j) x*Q(i,j)) (r^2 x^2)*L(i,j)]; end % 节点功率平衡方程 for j 2:N parent_idx find(tb j); if isempty(parent_idx) continue; end i fb(parent_idx); kk parent_idx; % j节点下游所有出线之和 child_idx find(fb j); outP 0; outQ 0; for h 1:length(child_idx) c tb(child_idx(h)); outP outP P(j,c); outQ outQ Q(j,c); end Constraints [Constraints, P(i,j) - Br(kk)*L(i,j) outP Pd(j)]; Constraints [Constraints, Q(i,j) - Bx(kk)*L(i,j) outQ Qd(j)]; end % 根节点电压、电压上下限 Constraints [Constraints, U(1) 1.0]; Constraints [Constraints, 0.9^2 U 1.05^2]; % 支路电流上限 Ilimit 0.5; Constraints [Constraints, L Ilimit^2]; % 目标函数网损最小 loss 0; for k 1:nb i fb(k); j tb(k); loss loss Br(k)*L(i,j); end % 求解 options sdpsettings(solver, cplex, verbose, 1); options.cplex.display.function 1; sol optimize(Constraints, loss, options); if sol.problem ~ 0 disp(sol.info); end % 结果 Popt value(P); Qopt value(Q); Uopt value(U); Lopt value(L); fprintf(总网损: %.4f pu %.2f kW\n, value(loss), value(loss)*Sbase/1e3); Vopt sqrt(value(U)); for k 1:N fprintf(节点 %2d 电压: %.4f p.u.\n, k, Vopt(k)); end按我多次运行的经验这段代码在CPLEX求解器下的典型结果是总网损约0.0203pu换算成有名值大约是202~203kW最低电压出现在节点18附近约为0.913 pu。这个结果和用传统牛拉法做潮流得到的IEEE33系统网损约202.67kW几乎一致。数值上比较吻合说明SOCP模型和松弛都是有效的。4.2 松弛间隙验证与电压曲线分析结果跑出来只是一半另一半工作是验证松弛是否精确。很多人忽略这一步整篇论文只有一个目标函数值审稿人一问“你的二阶锥松弛精确吗”就露馅了。所以我建议每次求解完都把松弛间隙打印出来。松弛间隙的计算很简单对每条支路求gap zeros(nb,1); for k 1:nb i fb(k); j tb(k); gap(k) Uopt(i)*Lopt(i,j) - (Popt(i,j)^2 Qopt(i,j)^2); end max_gap max(abs(gap)); disp([最大松弛间隙: , num2str(max_gap)]);IEEE33这个算例max_gap通常能到10^-7甚至10^-8量级说明每条支路的P-Q-U-L都严格满足原等式松弛后的解就是原问题的可行解。电压曲线则建议画出来看。节点18和节点33附近的电压最低这个趋势在IEEE33里是常态因为线路末端负荷密集而供电距离远。如果发现问题电压低于0.9就要考虑加无功补偿或者DG支撑这也正好是很多扩展研究的切入点。4.3 扩展场景加入分布式电源出力优化当你理解了基础算例之后加入DG只是很小的改动。假设在节点7、节点13和节点27接入三台分布式电源每台出力可调目标是让运行成本最低。第一个改动是给这些节点增加注入功率变量P_g和Q_g有出力上下限。第二个改动是节点功率平衡方程的右侧加上“-P_g”项因为向电网注入功率相当于抵消负荷。第三个改动是目标函数从单纯网损最小改成min Σ c_g * P_g λ * loss其中c_g是DG发电成本系数λ是网损权重。改动后的节点平衡方程形如Constraints [Constraints, P(i,j) - Br(kk)*L(i,j) outP Pd(j) - P_g(j)]; Constraints [Constraints, Q(i,j) - Bx(kk)*L(i,j) outQ Qd(j) - Q_g(j)];DG接入之后一个常见现象是节点电压被抬高某些原本电压偏低的末端节点反而会抬升到更安全的水平。但也要注意如果DG出力上限设得太大电压上限约束会起作用松弛间隙可能开始变大。这时候建议把目标函数末尾加一个微小的惩罚项比如epsilon * sum(L(:))epsilon取10^-6左右帮助优化器把解压回锥面。5. 从调试到迁移几个容易让算例翻车的细节5.1 CPLEX解SOCP时的报错处理二阶锥规划在CPLEX里属于“能直接吃”的问题类型不需要额外设置但报错依然常见。我遇到过的主要有以下这几类。第一类约束里出现了两个变量直接相乘比如写了U(i)*L(i,j) P(i,j)^2 Q(i,j)^2而YALMIP版本没有自动识别为旋转锥。这时CPLEX会报“Nonconvex”或“QCP with non-convex constraints”。解决方式是显式改成cone形式让YALMIP不会产生歧义。第二类求解器报“Numerical issues”。这几乎都是由于量纲没有统一或者某个变量数量级差异过大。解决方式不是去调求解器容差而是回去把模型重新标幺化。标幺化做好了数值问题能消失一大半。第三类报“Infeasible problem”。多半是电压下限和电流上限设置得太紧。IEEE33基准算例本身是可行的如果无解先检查电压上下限范围是否给了过窄的区间再检查支路电流上限是否小到了不合理的程度。你可以把Ilimit从0.5改到0.3看结果趋势如果还是无解问题大概率在约束书写方向。5.2 量纲、缩放与数值稳定性量纲问题我再强调一次因为这是新人最容易踩的坑。假设你不做标幺化直接在有名值模型里写锥约束U是电压平方量级约1.6e8L是电流平方量级约10^-6P和Q是潮流功率量级约10^6。这四个量放在同一个范数表达式里两个大数加一个小数小数直接被数值误差淹没CPLEX内部的预处理阶段就可能判定问题无解或不可行。标幺化之后所有变量都落在0.01到1.1之间范数约束的各个分量是一个量级上的求解器内部计算非常舒服。另一个细节是如果将来算例扩展到含光伏逆变器光伏的容量可能是MW级别同样要除以Sbase不要只对负荷做标幺化而忘记对DG做标幺化。5.3 从IEEE33迁移到更大规模网络时要注意什么IEEE33的33×33变量矩阵在YALMIP里跑起来很轻松但如果换到IEEE123节点或者几百个节点的真实馈线再用full矩阵就会产生大量无效的自由变量模型规模虚胖求解时间迅速上升。迁移时我建议改成按支路索引的向量形式比如P定义成n_branch维向量搭配首端/末端数组来写约束。另外辐射状网络的“父节点-子节点”关系必须在建模前梳理清楚。DistFlow方程默认功率从根节点流向叶节点如果你拿到的网络数据里支路端点顺序是乱的一定要先根据拓扑定向重新整理。否则功率平衡方程里的正负号会乱掉结果自然也是错的。如果再进一步做配电网重构把联络开关闭合网络变成弱环网二阶锥模型依然可以用但功率流向不再天然固定这时需要额外引入环路变量和方向判断复杂度会明显上升。建议先把IEEE33的辐射状SOCP算例彻底跑明白再去动弱环和重构的扩展。我在实际项目里用这套模型跑了大量DG选址和储能配置的迭代优化二阶锥松弛作为内层子问题求解稳定性和速度都让我很满意。如果你按照这篇文章的操作顺序复现还是遇到“锥不紧”“求解失败”或者“电压结果和潮流对不上”的情况建议先从松弛间隙和电压上下限两个方向排查这两处是绝大多数问题的根源。也欢迎私下交流你复现时的具体报错信息。