配电网最优潮流求解：二阶锥松弛与Matlab实现

从配电网调度说起。分布式光伏、储能、电动车充电桩大规模接入之后传统配电网最优潮流OPFOptimal Power Flow的计算压力一下子变大了。前阵子我在一个项目中需要用Matlab反复求解配电网的最优潮流目标很小就是在满足电压、容量约束的前提下让网损最小。但事情没有想象中顺利——用传统内点法去解初值稍微给得不好结果就跑到局部最优里去网络规模从33节点扩到123节点时求解耗时和收敛性都变得很难看。后来我把目光转向二阶锥松弛SOCPSecond-Order Cone Programming方法思路是把原问题中的非凸潮流等式做变量代换再松弛成二阶锥约束形成一个凸优化问题。凸问题有个很好的性质求解器能稳定地找到全局最优解而且用YALMIP、CVX这类建模工具在Matlab中实现起来非常顺手。这篇文章就把我踩过的坑、验证过的步骤和能直接复用的代码整理出来希望对做配电网优化调度、想做OPF研究或者正在被非线性求解器折磨的人有帮助。1. 配电网最优潮流为什么让人头疼非凸模型与求解瓶颈1.1 配电网OPF和输电网OPF差别到底在哪教科书上讲的OPF大多以输电网为背景网络结构是环网或网格状线路的R/X比值较小有功和无功的耦合度不高而且可以依托联络线做功率交换。传统求解方法像牛顿法、内点法在输电网场景下非常成熟因为网络拓扑相对规整、三相基本平衡建立的模型大多是连续可微的非线性方程组用迭代法能得到不错的结果。配电网则完全是另一副面孔。首先拓扑上大多是辐射状也就是树形结构支路之间有明显的“上游—下游”关系其次是线路参数配电网线路的电阻往往比电抗大得多R/X比值可能接近甚至超过1这让无功和有功的耦合非常强再加上单相用户多、三相不平衡、分布式电源接入后潮流方向可能反转整个问题的物理特征和输电网差别很大。用传统OPF模型去套配电网最直接的问题就是潮流方程里存在着电压幅值、电流幅值、有功功率、无功功率之间的大量非线性乘积项。这些项凑在一起造成可行域非凸求解器很容易被“卡”在局部最优。而配电网调度的现实需求是DG出力的随机波动、储能充放电切换、电动车充电负荷变化都要求更快更稳地算出调度策略不能容忍“这次收敛上次不收敛”的情况。1.2 非凸约束是怎么混进模型的来看一条典型的支路i→j电压方程可以写成V_j^2 V_i^2 - 2(r_ij·P_ij x_ij·Q_ij) (r_ij^2 x_ij^2)·I_ij^2支路功率方程则要满足P_ij - I_ij^2·r_ij ΣP_jk P_load_jQ_ij - I_ij^2·x_ij ΣQ_jk Q_load_j以及电流和功率的关系I_ij^2·V_i^2 P_ij^2 Q_ij^2注意最后这个式子I_ij^2 和 V_i^2 相乘等于 P_ij^2 Q_ij^2。这里不同变量之间是相乘关系还带着平方数学上这就是二次等式约束属于典型的非凸约束。非凸的含义用大白话说就是可行解集合的形状像个凹凸不平的山谷而不是一个光滑的碗。优化问题在这个山谷里找最低点很容易找到某个局部凹陷的“坑底”就停住了但那个坑底并不是整个区域的最低点。我在实际测试中感触很深的一件事是用FMINCON这类Matlab自带非线性优化器求解33节点配电网OPF从不同的初始点出发有时能得到网损约0.02 p.u.的解有时却在0.03 p.u.附近就收敛了结果差别很大。对于调度场景来说这样的不稳定性很难接受。所以在做研究和工程开发时我倾向于把模型整理成凸优化形式——这也正是二阶锥松弛的价值所在。2. 二阶锥松弛的核心思路从DistFlow到凸锥2.1 支路潮流方程DistFlow里的变量代换DistFlow模型是配电网最优潮流中非常经典的一种建模方式。它不直接对节点导纳矩阵做全维度的复数运算而是把每条支路的有功、无功、电流、电压之间的关系一条条写出来非常贴合配电网的辐射状拓扑结构。为了把非凸项消掉第一步是变量替换。让U_i V_i^2L_ij I_ij^2这样电压方程可以改写成U_j U_i - 2(r_ij·P_ij x_ij·Q_ij) (r_ij^2 x_ij^2)·L_ij功率方程变成P_ij - L_ij·r_ij ΣP_jk P_load_jQ_ij - L_ij·x_ij ΣQ_jk Q_load_j到这一步除了最后一个电流功率关系式其余约束都已经线性化了。变量U、L、P、Q之间不再出现电压与电流的直接乘积剩下的那双线性等式就是L_ij·U_i P_ij^2 Q_ij^2这个式子和原来的 I_ij^2·V_i^2 P_ij^2 Q_ij^2 等价但变量已经换成了更容易处理的U和L。不过它依然是一个非凸等式。下表能看清转换前后的变化原始形式变量代换后V_j^2 V_i^2 - 2(rPxQ) (r²x²)I²U_j U_i - 2(rPxQ) (r²x²)LP_ij - I²r 下游有功需求P_ij - Lr 下游有功需求I²V² P² Q²LU P² Q²仍非凸到这里整个问题就剩下这一颗“钉子”了。2.2 最后一个非凸等式如何变成锥约束等式约束 LU P² Q² 可以松弛成一个不等式L_ij·U_i ≥ P_ij² Q_ij²这个放宽的合理性在于原本必须在曲面上取值现在允许在曲面“上方”和“外部”取值可行域从一条曲线扩展成一个凸区域。数学上它等价于一个标准的二阶锥形式‖ [2P_ij; 2Q_ij; L_ij - U_i] ‖₂ ≤ L_ij U_i在YALMIP里可以直接这么写norm([2*P_ij; 2*Q_ij; L_ij - U_i]) L_ij U_i有些教材也称之为“旋转二阶锥”rotated second-order cone本质是一样的。这样处理后整个模型变成一个二阶锥规划问题属于凸优化范畴可以用Mosek、SDPT3、SeDuMi等专业求解器高效求解。需要说明的是松弛后的可行域比原问题更大所以松弛问题的最优目标值理论上不会比原问题更差。换句话说如果松弛后的解正好满足等式成立那么它一定也是原问题的全局最优解。这引出了下一个关键问题什么情况下我们能放心地用松弛解代替原问题解。2.3 什么情况下松弛是精确的辐射网结构是最大底气对辐射状配电网只要负荷有功和无功有合理上界且目标函数对电流是单调递增的典型如网损最小大量理论研究表明二阶锥松弛是精确的即在最优解处那些松弛成不等式的约束会自然取到等号LU P² Q²。直观理解也不难目标函数要最小化网损网损正比于 L_ij·r_ij求解器会倾向于把 L_ij 压小同时电压也不能太低U_i 太小时电压约束会卡住。在这两股力的拉扯下L_ij·U_i 会尽可能逼近 P_ij² Q_ij²所以不等式是“紧”的。但有两类情况要警惕一是网络存在环网比如联络开关闭合形成了闭环二是目标函数不再是网损最小这种“电流厌恶型”例如某些多目标优化里加入了电压偏差最小化。这些情况下松弛精度可能下降需要额外检查锥间隙具体方法我放在第5章详细讲。3. Matlab代码实现YALMIP与CVX两个建模路线3.1 环境准备优先推荐的工具链搭配在Matlab里做二阶锥优化最顺手的组合是 YALMIP Mosek。YALMIP是一个建模工具箱它的价值在于把优化模型和底层求解器解耦你不需要关心每种求解器的语法差异只需要用统一的建模语句描述变量、约束、目标然后告诉它调用哪个求解器即可。安装上没什么特别复杂的下载YALMIP解压后把文件夹加入Matlab路径再下载Mosek安装包安装后同样加入路径。检测是否装好可以用yalmiptest如果只想先跑通流程不想一开始就上商用求解器可以装SDPT3这是免费的YALMIP直接支持。但我的建议是如果你的项目会涉及上百节点以上的配电网尽量申请Mosek的学术授权速度和数值稳定性完全不是一个量级。CVX是另一个备选方案它的建模语法同样简单直观。不过CVX对二阶锥约束的表达稍微绕一点而且调试信息不如YALMIP直观。两条路线我在下面的代码里都给出来供你根据自己的习惯选。3.2 最小可复现案例3节点配电网SOCP-OPF完整代码为了让模型不依赖庞大的IEEE数据文件我先把一个3节点辐射网写完整。网络结构是节点1是松弛节点通过线路1-2连接节点2再通过线路2-3连接节点3。所有量都采用标幺值线路阻抗r0.01x0.01节点3接有功0.3和无功0.1的负荷节点2接有功0.1和无功0.05的负荷。% SOCP-OPF: 3节点辐射配电网YALMIP Mosek/SDPT3 clear; clc; % 支路数据: [首端, 末端, 电阻, 电抗] branch [ 1 2 0.01 0.01 2 3 0.01 0.01 ]; % 负荷数据: [节点, 有功, 无功] loads [ 2 0.10 0.05 3 0.30 0.10 ]; nB 3; % 节点数 nL size(branch, 1); % 支路数 % 决策变量 U sdpvar(nB, 1); % 电压幅值平方 L sdpvar(nL, 1); % 电流幅值平方 P sdpvar(nL, 1); % 支路有功 Q sdpvar(nL, 1); % 支路无功 % 约束集合 C []; % 松弛节点电压 C [C, U(1) 1.0]; % DistFlow等式约束 for k 1:nL i branch(k, 1); j branch(k, 2); r branch(k, 3); x branch(k, 4); % 末端负荷若非负荷节点则取0 Pd_j 0; Qd_j 0; idx find(loads(:,1) j); if ~isempty(idx) Pd_j loads(idx, 2); Qd_j loads(idx, 3); end C [C, U(j) U(i) - 2*(r*P(k) x*Q(k)) (r^2 x^2)*L(k)]; C [C, P(k) - L(k)*r Pd_j]; % 若节点j下游还有支路右侧需追加下游支路功率 C [C, Q(k) - L(k)*x Qd_j]; % 二阶锥松弛 C [C, norm([2*P(k); 2*Q(k); L(k) - U(j)]) L(k) U(j)]; end % 目标网损最小 objective sum(L .* branch(:,3)); % 求解 ops sdpsettings(solver, mosek, verbose, 0); optimize(C, objective, ops); % 输出 P_val value(P); Q_val value(Q); U_val value(U); L_val value(L); loss value(objective); fprintf(网损: %.6f p.u.\n, loss); fprintf(节点电压幅值: ); fprintf(%.4f , sqrt(U_val)); fprintf(\n); % 验证锥间隙 gap 0; for k 1:nL j branch(k, 2); gap_k L_val(k)*U_val(j) - (P_val(k)^2 Q_val(k)^2); gap max(gap, abs(gap_k)); end fprintf(最大锥间隙: %.2e\n, gap);这段代码的运行逻辑很简单YALMIP把锥约束直接写成norm形式求解器会识别为二阶锥问题。输出结果中我特意加了锥间隙的验证这是判断松弛是否精确的关键指标——如果间隙在1e-6量级说明这个松弛解就是原问题的最优解如果间隙到了1e-2量级那就要小心了。3.3 扩展到IEEE 33节点数据结构与建模思路3节点能跑通33节点只是数据规模变大建模逻辑完全一致。IEEE 33节点系统的常见Matlab数据格式是两张大矩阵bus矩阵第一列节点编号第三列有功负荷第四列无功负荷branch矩阵第一列首端节点第二列末端节点第三列电阻第四列电抗把上一节的3节点代码稍微改一下即可。核心在于将约束中“末端负荷”的部分改为循环累加下游功率% 对每条支路k找到末端节点j的所有下游支路集合 down find(branch(:,1) j); C [C, P(k) - L(k)*branch(k,3) sum(P(down)) Pd_j]; C [C, Q(k) - L(k)*branch(k,4) sum(Q(down)) Qd_j];注意这里的P和Q都是按“首端流向末端”的方向定义的有向变量YALMIP会自动把它们当成普通优化变量处理。迭代累加的顺序不需要手工排序因为这是等式约束求解器会一次性把整个方程组解出来。用IEEE 33节点系统测试配置为Mosek我这边通常0.05秒到0.2秒就能收敛而用传统非线性求解器如FMINCON去解同样的模型耗时会高出很多而且还有初值敏感的问题。SOCP在效率上的优势非常明显。4. 求解器选型与数值调优实测对比与工程建议4.1 主流求解器在配电网SOCP上的表现我先后试过Mosek、SDPT3、SeDuMi和Gurobi也稍微看过ECOS。粗略感受如下求解器许可证33节点耗时123节点感受数值稳定性备注Mosek商用/学术免费0.05~0.2s依然较快非常好我的首选SDPT3免费0.2~0.5s明显变慢中等适合中小规模SeDuMi免费0.3~1s较慢中等老旧不推荐Gurobi商用/学术免费0.1~0.3s较快好主攻LP/MIPSOCP略弱于MosekECOS开源依赖接口配置不稳定一般嵌入式场景为主上表的时间是典型量级具体数值会随机器配置和YALMIP版本浮动但它反映的趋势是真实的Mosek在处理配电网SOCP问题时无论求解速度还是数值鲁棒性都明显领先。SDPT3作为免费选项解决几十节点的小规模演示、学术验证完全够用但如果要做全天96时段的滚动优化我强烈建议上Mosek。4.2 让求解更稳的几个工程技巧第一所有电气量尽量用标幺值。配电网的线路电阻可能只有0.01 p.u.但如果用有名值电压可能是10kV电流几百安培功率几兆瓦变量之间量级差到6个数量级以上求解器的数值预处理会很难看。第二注意变量的初始化。虽然凸优化不像非线性优化那样对初值极其敏感但一个合理的初始点能显著减少内点法的迭代次数。YALMIP中可以用assign命令给变量赋初值assign(U, ones(nB,1)); assign(L, 0.01*ones(nL,1));第三适当调整求解器容差。Mosek里的相对间隙容差可以用参数控制ops sdpsettings(solver, mosek, mosek.MSK_DPAR_INTPNT_CO_TOL_REL_GAP, 1e-7);SDPT3则可以通过sdpsettings直接设置ops sdpsettings(solver, sdpt3, sdpt3.maxit, 1000);第四不要忽略YALMIP在建模层面的预处理。如果约束中存在大量元素级乘法或循环建议先把网络数据矢量化再用矩阵形式一次性构造约束避免几千行循环带来的建模开销。在Matlab中循环不是不能写但当节点数上万时建模时间可能比求解时间还长。这些技巧看着不起眼但在工程系统里一个“能跑通“的模型和一个“能稳定跑大量场景”的模型之间往往就差在这些细节上。5. 实操阶段的五个高频坑5.1 相角与参考节点最容易被忽略的物理约束在SOCP模型中经过变量代换后约束里已经看不到电压相角了。这不代表相角不重要。实际物理系统中必须要有一个参考节点提供相角基准通常就是松弛节点。如果你在建模时忘了固定松弛节点的电压幅值求解器可能会给出一个“电压全为0.8”的解目标函数倒是很小可物理上一塌糊涂。我在代码里用 U(1) 1.0 固定了参考节点电压幅值平方这个步骤不能省。如果你的网络中有多个电源节点也需要考虑是否需要把某些PV节点处理成给定电压幅值的形式。配电网中更多是PQ节点处理方式相对简单。5.2 锥松弛不精确的判别与处理前面提到过松弛并不保证在任何场景下都精确。判断是否精确的最直接方法是计算锥间隙gap L_ij·U_j - (P_ij² Q_ij²)在求解后遍历所有支路找到最大间隙。如果最大间隙超过1e-4具体阈值取决于你的精度需求说明这个松弛解可能不是原问题的可行解直接用可能带来误导性结果。处理不精确的情况常见思路是“惩罚法”在目标函数里给锥间隙加一个小惩罚项让求解器在优化的过程中主动把锥间隙压小。也可以用迭代收紧策略先解一次SOCP找出间隙最大的支路在原问题里加入针对这条支路的割平面或局部约束再迭代求解。5.3 电流、视在功率与容量约束的建模混淆配电网的线路容量约束原问题是 I_ij² ≤ I_max²写成新变量就是 L_ij ≤ I_max²这个非常直接。但很多人会在建模时把容量约束写成“流过支路的视在功率 S_ij ≤ S_max”然后又把 S_ij 近似成 P_ij这在R/X比较高、无功不可忽略的配电网中会产生明显误差。正确做法是如果要用功率来限制容量应该用 S_ij² P_ij² Q_ij²再配合电压约束一起考虑。但由于P和Q都是优化变量这个功率约束恰好可以复用锥约束的形式。实际工程中我建议优先使用电流约束因为它和导体发热直接相关物理意义更明确。5.4 数值病态当解出来是负电压电压平方U理论上必须非负。如果求解器返回的U中出现了明显的负值或者电压幅值sqrt(U)出现NaN首先检查是不是数据单位不一致。常见错误是把功率基值写成1MW但电阻和电抗却用了Ω/相的有名值导致电压方程里的单位混乱。其次检查约束中是否缺少U的下界。理论上二阶锥约束会限制U不能过小但在某些边界条件下数值求解器可能踩过界。稳妥的建模方式是对每个节点显式加入C [C, U 0.81]; % 电压不低于0.9 p.u. C [C, U 1.21]; % 电压不高于1.1 p.u.这样既符合调度规程也让数值求解更稳定。5.5 与非线性OPF对拍验证模型正确性的方法这是我个人强烈建议的步骤。任何SOCP-OPF代码写完先别急着扩展到大型网络先在IEEE 33节点或某个你熟悉的小网络上用原始非线性潮流模型比如直接用Matpower或FMINCON求解算一次再和SOCP结果对比网损、节点电压、支路潮流。为什么这样做因为SOCP模型在推导过程中做了变量代换如果不小心写错方程比如电压降落公式里的符号、下游功率累加的顺序整个结果会偏离很远。对拍是一种极好的调试手段“网损误差小于0.1%”这个标准基本能证明建模没有大问题。我在测试中发现SOCP结果和Matpower算出的交流OPF结果网损差异通常在0.01%以内电压幅值几乎重合。这给了我在实际项目中用它做日前调度的底气。6. 从基础模型走出去DG、储能与三相不平衡场景6.1 接入分布式光伏和风机时如何改模型分布式光伏和风机在配电网OPF中通常建模为负的有功负荷同时叠加无功出力范围。如果是定功率因数控制就把无功设为一个固定的比例如果是逆变器可调无功要加入容量约束P_dg² Q_dg² ≤ S_dg_max²这个圆形约束也是凸的可以直接作为一个二阶锥约束加进模型。关键是DG所在节点的功率平衡方程要从“只等于负荷”改成“负荷减去DG出力”P_ij - L_ij·r_ij 下游有功需求 - P_dg我在3节点的示例代码里没有写DG但扩展起来就多几个变量和约束逻辑完全一致。6.2 多时段储能运行与SOCP扩展储能的引入会把单时段OPF变成多时段问题。设储能节点s在时段t的充电功率为P_ch(t)放电功率为P_dis(t)则能量平衡约束为E(t1) E(t) η_ch·P_ch(t) - P_dis(t)/η_disE有上下限P_ch和P_dis也有上下限且两个变量不能同时非零。这个“不能同时充放电”的约束是非凸的但实际中只要目标函数合适求解器通常不会出现同时充放电的无意义解——因为那会增加网损、减少效率对目标函数不利。如果仍然担心可以通过设置充电费用或加入小惩罚项来规避。多时段模型会导致变量数量成倍增长但二阶锥规划的求解器对这类结构化问题的处理能力很强。我做96时段的33节点算例时Mosek求解时间仍在可接受范围内远快于非线性方法。6.3 三相不平衡配电网的前沿扩展三相不平衡是配电网的常态。扩展思路是把每个节点的电压、每条支路的功率都拆成A/B/C三相DistFlow方程逐相写支路间再加互感耦合项。变量规模大概变成单相模型的3到6倍但锥结构不变依然可以用SOCP求解。这项技术我目前还在研究中目前公开文献里已经有比较成熟的三相DistFlowSOCP框架。如果你要做这个方向建议先把单相模型的代码和验证跑熟再逐步扩展到三相否则调试难度会非常大。回到工程层面说点实在的。用二阶锥松弛做配电网最优潮流最大的价值不是让模型变得更“高级”而是让调度计算从“碰运气式收敛”变成“确定性收敛”。我在项目里最深刻的体会是凸优化让整个团队可以不再纠结求解器是否卡在局部最优而是把精力放在模型本身的业务约束上。建议初学的朋友一定先把3节点代码逐行跑通然后把锥间隙验证写进自己的工具箱最后再去追求规模和复杂度。这条路走顺了后面接DG、接储能、接三相不平衡都会顺畅很多。