开写之前我想先说明一点这篇文章不是那种“贴一段说明书式的代码然后让你自己琢磨“的凑数笔记。我会把二阶锥松弛在配电网最优潮流里的来龙去脉讲清楚再给出Matlab代码骨架最后把调通过程中容易踩的坑挨个拆开。你现在看到的这套东西是我自己从复现文献模型到跑通IEEE 33节点算例的真实路径希望能帮你少走弯路。为什么配电网最优潮流需要二阶锥松弛1. 配电网最优潮流为什么比输电网的难解1.1 配电网OPF和输电网OPF的“硬”不一样很多刚接触最优潮流OPF的人第一反应是“这不就是潮流加个优化目标嘛”。确实输电网OPF已经有了很成熟的处理方式内点法、原对偶方法二十年前就能解到几千节点规模。但配电网的情况完全两样。输电网X/R比高、电压等级高潮流方程可以近似成直流潮流处理模型差不多是线性的求解器很喜欢这种问题。配电网呢10kV、12.66kV这种等级线路电阻和电抗同量级甚至电阻更大直流潮流简化法在这里基本失效你必须回退到完整的交流潮流模型。交流潮流模型就不好惹了。主动配电网里的节点电压幅值、支路电流幅值、有功无功潮流之间全是非线性耦合关系电压降落的二次项、电流平方项、功率乘积项搅在一起整个可行域是一个非凸集合。目标函数再叠加上分布式电源出力、网损、电压偏差优化问题就变成了一个大规模非凸非线性规划NLP。这类问题对初值极度敏感陷入局部最优的概率不低而且没有一套通用的手段判断你得到的解到底是不是全局最优。1.2 传统NLP方法在配电网模型上的窘境我最初尝试直接用非线性规划的解法器去解配电网OPF用的还是经典的交流潮流方程并把所有等式约束原样写进去。这么做的问题是第一对初值要求高多点起算经常得到不同的结果你根本不知道哪一个是全局最优第二当配电网中有大量分布式电源、功率双向流动时雅可比矩阵条件数变差迭代步长很容易震荡甚至直接发散第三计算时间不可控。更麻烦的是配电网规划问题里常常还带着0-1变量比如网络重构中的开关状态、分布式电源的选址定容。一旦把非凸潮流约束和整数变量放在一起就成了一个非凸混合整数非线性规划MINLP这个问题连当前最前沿的商业求解器都不敢拍胸脯保证收敛。这也是为什么大家都想找一种把非凸约束转化成凸约束的办法——凸化之后哪怕带整数变量退化成MISOCPGurobi和Cplex也能稳定求解。二阶锥松弛Second-Order Cone ProgrammingSOCP就是这类凸化方法里应用最广的一个。1.3 二阶锥松弛的直观理解把“必须压在线上”变成“允许落在锥内”我学习SOCP时觉得最关键的思维转变是不要把它当成一个“近似技巧”而要理解它“松”掉了什么。原始潮流方程里有一个核心等式支路电流幅值的平方乘以首端电压幅值的平方等于该支路有功平方加无功平方。这个等式就是那条“必须压在线上”的曲线约束物理上强制功率、电压、电流三者完全自洽。二阶锥松弛做的事很简单把等式“”改成一个不等式“≥”也就是允许电流电压的组合落在锥形区域内部而不是必须落在锥面上。你可能会问这不是把约束放宽了吗得到的解还是原问题的解吗这就牵出了一个更微妙的性质在目标函数单调地拉大或压缩某个变量时最优解会自动回到锥面上让松弛后的最优解与原问题严格一致。这个性质才是SOCP能大规模落地的前提后面我会重点讲它成立的条件和失效的场景。2. 从DistFlow方程到锥约束的数学推导2.1 先约定记号每条支路的状态量要推导先得统一符号。考虑一个辐射状配电网节点集合为N支路集合为E根节点变电站出口编号为1。i、j支路(i,j)的首端和末端节点r_ij、x_ij支路电阻和电抗P_ij、Q_ij从节点i流向节点j的有功和无功潮流V_i节点i的电压幅值v_i V_i^2节点电压幅值的平方I_ij支路电流幅值l_ij I_ij^2支路电流幅值的平方P_inj,i、Q_inj,i节点i的有功、无功注入分布式电源或变电站注入减去负荷。这个记号和Baran和Wu在1989年提出的DistFlow模型保持一致配电网SOCP相关文献基本都沿用这套写法。符号统一的好处是你在读论文时不会因为流量方向定义不同而绕晕。2.2 DistFlow方程组的三个核心关系辐射状配电网中每条支路(i,j)满足三个方程。第一个是电压降落方程v_j v_i - 2(r_ij P_ij x_ij Q_ij) (r_ij^2 x_ij^2) l_ij这个方程的来源是欧姆定律和相量运算展开之后取幅值平方再把高次项整理成电流平方的形式。第二个是节点功率平衡方程说的是流入节点的功率加节点注入等于流出节点的功率加节点负荷sum_{k∈child(j)} P_jk P_inj,j - P_load,j (P_ij - r_ij l_ij)注意这里P_ij是上游支路流入j的功率扣掉线路损耗r_ij l_ij之后剩余部分才能参与本地平衡。第三个方程描述功率、电压和电流的关系l_ij (P_ij^2 Q_ij^2) / v_i问题就出在这第三个方程上。它是一个分数非线性等式而且关于v_i和l_ij是双线性的写成等价形式是l_ij v_i P_ij^2 Q_ij^2显然是非凸的。整个OPF模型之所以难解根源就在这个等式。2.3 把非凸等式改写成标准二阶锥约束SOCP的标准做法分两步。第一步做变量替换把v_i和l_ij直接设为优化变量而不是再写V_i^2和I_ij^2。这样电压降落方程和功率平衡方程中出现的都是线性项或简单二次项只有l_ij v_i P_ij^2 Q_ij^2这个关联等式是非凸的。第二步是松弛。把等式变成不等式l_ij v_i ≥ P_ij^2 Q_ij^2为什么这个不等式能改写成二阶锥约束把两边展开一下。定义向量x [2P_ij; 2Q_ij; l_ij - v_i]然后构造2-范数约束|| [2P_ij; 2Q_ij; l_ij - v_i] ||_2 ≤ l_ij v_i两边平方4P_ij^2 4Q_ij^2 (l_ij - v_i)^2 ≤ (l_ij v_i)^2整理后得到4(P_ij^2 Q_ij^2) ≤ 4l_ij v_i正好就是l_ij v_i ≥ P_ij^2 Q_ij^2。一个看似复杂的非线性等式就这样转换成了Yalmip里可以直接写给求解器的cone()约束。2.4 松弛为什么能在最优解处自动收紧理论上说松弛后的可行域比原问题大最优值不会比原问题差但可能比原问题“好得离谱”——因为原问题不可行的点松弛后可能可行。如果松弛后的最优解在锥的内部那它就不是原问题的可行解得到的结果没有物理意义。幸运的是在辐射状配电网中SOCP松弛经常是精确的也就是最优解会落在锥面上。这里头有个物理直觉。目标函数如果是网损最小就可以写成 sum(r_ij l_ij)由于电阻r_ij是正数目标函数会尽可能把每个l_ij往下压。而锥约束允许的范围是l_ij ≥ (P_ij^2 Q_ij^2)/v_i能压多低压多低。于是每个支路最后都会被压到锥的边界上松弛精确性天然成立。同理目标函数如果包含正系数倍的支路电流平方项也能达到类似效果。目前文献里比较公认的结论是当网络是辐射状、目标函数对电流/网损有单调递增趋势、并且不存在过紧的电流上限约束时SOCP松弛基本是精确的。但这不是免费的午餐后面我会讲怎么验证、失效了怎么办。3. Matlab代码实现从33节点数据到SOCP求解3.1 数据准备与基准值单位先说明一点我不建议依赖Matpower自带的case文件来跑SOCP-OPF即便Matpower里有配电网形式的数据格式依然围绕输电网习惯设计对你的模型构建没多大帮助。更直接的做法是手动把IEEE 33节点系统的数据整理成两个矩阵。第一个矩阵是bus_data两列分别是有功负荷PdkW和无功负荷Qdkvar。第二个矩阵是branch_data四列分别是首端节点、末端节点、线路电阻rohm、线路电抗xohm。IEEE 33节点标准数据在论文附录和很多公开资源里都能找到直接把那两张表粘贴到矩阵里即可。下面给一个片段格式% bus_data: [节点编号, Pd(kW), Qd(kvar)] % 1号节点是根节点负荷为0 bus_data [ 1, 0, 0; 2, 100, 60; 3, 90, 40; 4, 120, 80; 5, 60, 30; ... ]; % branch_data: [首端节点, 末端节点, r(ohm), x(ohm)] branch_data [ 1, 2, 0.0922, 0.0470; 2, 3, 0.4930, 0.2511; 3, 4, 0.3660, 0.1864; ... ];很多人在这第一步就翻车并不是拿不到数据而是单位搞错了。配电网你直接用Ω和kvar去建模数值尺度会差好几个量级Cplex和Gurobi的数值引擎会非常不稳定。所以必须标幺化。我习惯取基准电压Vbase 12.66 kV基准功率Sbase 10 MVA。这样基准阻抗Zbase Vbase^2 / Sbase 16.02 Ω负荷和注入也都除以Sbase换算成标幺值。3.2 Yalmip模型变量、目标函数、约束环境是Matlab配合Yalmip工具箱求解器可以选Cplex、Gurobi或者Mosek。Yalmip的作用是把优化问题建模成求解器能懂的格式尤其对SOCP约束它提供了直观的cone函数。下面是变量定义和约束构建的核心代码nb size(bus_data, 1); % 节点数 nl size(branch_data, 1); % 支路数 Vbase 12.66e3; % V Sbase 10e6; % VA Zbase Vbase^2 / Sbase; % ohm r_pu branch_data(:,3) / Zbase; x_pu branch_data(:,4) / Zbase; Pd_pu bus_data(:,2) / Sbase; Qd_pu bus_data(:,3) / Sbase; % 变量定义 v sdpvar(nb, 1); % 节点电压幅值平方 I2 sdpvar(nl, 1); % 支路电流幅值平方 P sdpvar(nl, 1); % 支路首端有功 Q sdpvar(nl, 1); % 支路首端无功 Pg sdpvar(nb, 1); % 节点有功注入 Qg sdpvar(nb, 1); % 节点无功注入 C []; % 支路约束电压降落方程 二阶锥松弛 for k 1:nl i branch_data(k,1); j branch_data(k,2); C [C, v(j) v(i) - 2*(r_pu(k)*P(k) x_pu(k)*Q(k)) ... (r_pu(k)^2 x_pu(k)^2) * I2(k)]; C [C, cone([2*P(k); 2*Q(k); I2(k) - v(i)], I2(k) v(i))]; end % 节点功率平衡 for j 1:nb inbound find(branch_data(:,2) j); outbound find(branch_data(:,1) j); C [C, -sum(P(inbound) - r_pu(inbound).*I2(inbound)) ... sum(P(outbound)) Pg(j) - Pd_pu(j)]; C [C, -sum(Q(inbound) - x_pu(inbound).*I2(inbound)) ... sum(Q(outbound)) Qg(j) - Qd_pu(j)]; end % 根节点电压固定变电站注入自由 C [C, v(1) 1.0]; C [C, Pg(1) 0]; C [C, Qg(1) 0]; % 无DG节点注入为0 C [C, Pg(2:end) 0]; C [C, Qg(2:end) 0]; % 电压上下限 0.95~1.05 p.u. C [C, 0.95^2 v 1.05^2]; % 支路电流上限按实际导线载流设定 Imax_pu 1.2; C [C, I2 Imax_pu^2]; % 目标函数网损最小 objective sum(r_pu .* I2);这里有几个细节我要特别强调。第一功率平衡方程里inbound变量是支路末端等于当前节点j的那些支路也就是从上级注入j的支路outbound则是从j流向子节点的支路。由于配电网是辐射状任意非根节点上方只有一条支路所以inbound多数时候只有一个元素但代码写成通用形式没有坏处。第二节点注入的正方向是“注入网络”所以根节点Pg(1)为正值表示变电站向网络输送功率如果未来考虑分布式光伏大发导致根节点倒送功率可以把约束改成Pg(1)的上下限。第三样例子中我把2号节点及以后所有节点的Pg设为0意思是不考虑分布式电源如果你要加DG就把该节点的Pg改为带上下限的变量而不是固定为0。3.3 求解器设置与结果恢复模型构建完后用下面这段代码调用求解器ops sdpsettings(solver, cplex, verbose, 2, cplex.mip.tolerances.mipgap, 1e-4); optimize(C, objective, ops); % 结果恢复 V_opt sqrt(value(v)); P_opt value(P); Q_opt value(Q); I2_opt value(I2); loss_opt sum(r_pu .* value(I2)) * Sbase; % 转成实际功率单位W % 电压分布图 figure; plot(V_opt, -o, LineWidth, 1.5); xlabel(节点编号); ylabel(电压幅值 (p.u.)); grid on;运行成功后你会看到33节点系统在网损最小化目标下末端节点电压比根节点低支路电流平方值也落在合理范围。用Cplex求解通常不到一秒就能收敛这比直接解NLP体感快得多。3.4 代码里最容易写错的三处第一Yalmip里的变量名尽量避免用字母l它和数字1长得太像调试时很难分清。上面代码统一用I2表示支路电流幅值平方就是不想给自己找麻烦。第二cone的写法。很多人第一次写SOCP约束时会把它写成norm([2P;2Q;I2-v]) I2v这在Yalmip里也支持但cone函数的签名是cone(x, y)表示norm(x)y别搞反了参数顺序。第一参数必须是向量第二参数必须是标量否则Yalmip会直接报错。第三功率平衡方向的符号。我在代码里明确用“注入网络为正方”的约定这是DistFlow文献的标准约定。如果你从网上找到一段代码它的符号可能用的是“负荷为正方”最后结果P和Q的值符号会整体相反但你如果不检查支路潮流的物理方向很可能意识不到这个问题。4. 松弛不精确怎么办验证与补救4.1 松弛是否精确的判断标准SOCP方法最大的隐患就是松弛不精确但你可以通过一个简单办法验证把SOCP求得的节点注入和根节点电压作为输入重新做一次精确的交流潮流计算得到网损然后和SOCP目标函数值对比。如果两者差值在0.1%以内说明松弛精确可以放心用如果差距过大说明最优解落在锥内部SOCP结果不可直接用。我通常在代码里加这么一段校验逻辑% flow_loss 是从潮流计算得到的网损 % loss_opt 是SOCP目标函数值 diff_loss abs(flow_loss - loss_opt); if diff_loss / flow_loss 1e-4 disp(SOCP松弛精确); else disp(SOCP松弛不精确需要检查约束和边界); end如果你不想自己写潮流计算可以设定好潮流数据格式后调用Matpower的runpf把SOCP得到的节点注入做一次潮流比较网损。个人经验是至少对每一篇论文的算例都做一次这样的校验不要上来就默认SOCP一定精确。4.2 哪些条件下SOCP会给出错误结果我复盘自己踩过的坑和文献里反复强调的案例SOCP松弛失效主要集中在几类场景。第一目标函数不“压”电流。比如目标只是最小化变电站购电成本而购电成本又只跟根节点有功有关跟支路电流没有直接单调关系松弛后的解可能停留在锥内部得到的网损偏低电压整体偏“美观”。这时候你就要小心了。第二电流上限约束太紧。如果某条支路的电流上限恰好卡在最优点上导致目标函数再怎么压也压不到锥面松弛就可能失效。辐条状网络中支路电流上限是SOCP精确性的边界条件。第三网络不是严格的辐射状。配电网实际运行时可能有联络开关闭合形成环网。在环网情况下SOCP松弛精确性不再有理论保证解出来的结果可能完全不符合潮流方程。这是很多人看SOCP文献时忽略的一个大前提。第四节点负荷过低。负荷接近零时潮流方程本身接近病态锥约束可能被轻易满足不需要贴着锥面走。4.3 可行的修复手段一旦发现松弛不精确有几种补救思路可以尝试按从简单到复杂的顺序说。第一种是加罚项。在目标函数里加一个很小的正数乘以所有支路I2之和也就是min sum(r_ij l_ij) ε*sum(l_ij)这样目标函数被迫继续“压”电流把解往锥面上推。ε一般取1e-4到1e-2之间但要注意罚项别大到扭曲原目标。第二种是收紧边界。检查一下是不是某些约束边界设置过宽比如支路电流上限给的太大给了锥约束“偷懒”的空间。适当收紧电流上限或者把节点电压上下限从0.95~1.05改到0.98~1.02往往能让松弛重新精确。第三种是交替迭代。解完SOCP后把得到的节点电压代入原始潮流方程算出精确的支路电流和网损然后把这些值作为新的初始点再解一次SOCP。这种序贯SOCP在工程中很常见虽然不是严格全局最优但实际效果稳定。第四种是切换到SDP或者更精确的凸松弛。当SOCP失效时半定松弛SDP在某些场景下更紧代价是计算量大不少。如果网络规模不大可以试试SDP作为对照。5. 调通代码后一定会遇到的几类实战问题5.1 求解器选型Cplex、Gurobi、Mosek、Sedumi怎么选理论上SOCP是凸优化任何支持二阶锥的商业求解器都能解。但我实际测试下来同样一个33节点网损最小问题不同求解器的求解时间和对Yalmip模型格式的容错度还是有差异。求解器求解速度大规模表现与Yalmip兼容性备注Cplex快优秀好配电网MISOCP场景常用Gurobi快优秀好参数多收敛行为稳定Mosek很快优秀好SOCP内点法非常专业Sedumi较慢一般可免费适合小规模教学我个人的建议是学习阶段用Cplex或Gurobi正式科研时再看看Mosek。配电网SOCP-OPF的问题规模一般不大33节点、123节点Cplex就足够用了如果做上百个时间断面的时序优化Mosek在SOCP内点法上的性能优势会更明显。5.2 量纲和数值尺度引起的玄学报错Yalmip本身能处理一部分量纲混乱问题但求解器内部对数值尺度是有偏好的。有次我把负荷数据直接写成kW注入写成标幺结果Cplex报infeasible我检查了很久才发现两条数据一个用kW一个用标幺。建议一开始就把所有电气量统一成标幺值不要嫌麻烦。根节点电压固定为1.0电压上下限写0.95^2~1.05^2支路电阻电抗都除以Zbase。这样整个模型的数值量级都在0.01到10之间对求解器非常友好。另一个常见问题是支路电流上限设置得过松。如果你不设置I2的上限模型可能会解出一个电流特别大的荒唐结果因为网损最小目标里电阻乘电流平方虽然会限制电流但某些边界场景下电压约束压力更大电流可能被顶到非物理上限。给每条支路一个合理的载流量上限既是物理约束也帮助数值稳定。5.3 求解时间长与迭代慢的原因对大系统SOCP虽然凸但求解时间不总是秒级。如果发现求解很慢第一步检查是不是Yalmip把SOCP误判成了非线性规划。Yalmip会检查约束类型但有时候前期建模不规范比如把原本应该是线性约束的电压降落方程写成了间接形式会导致求解器走NLP路径。你可以通过调用ops.solver查看实际选中的求解器如果显示fmincon之类大概率是被Yalmip当NLP处理了。第二步看模型是否有大量冗余变量。比如节点功率平衡用循环生成时每个循环都重新做一次find(branch_data(:,2)j)模型规模小时无所谓上千节点时积累的find耗时不可忽略。可以把inbound/outbound在循环外提前算好存成元胞数组再进循环。5.4 结果后处理与可视化经验SOCP直接解出来的变量是v和I2一个是电压幅值平方一个是电流幅值平方后处理时一定要记得开根号。我刚开始跑的时候直接拿value(v)去画电压分布图画出来全是0.9到1.1附近的数还觉得奇怪——后来才反应过来这是平方值不是电压本身。可视化除了画电压分布曲线我还会画一个“松弛间隙”的柱状图。对于每条支路计算I2_opt(k)*v_opt(i) - (P_opt(k)^2 Q_opt(k)^2)再除以I2_opt(k)*v_opt(i)。这个值表示该支路离锥面的相对距离所有支路的间隙都接近0说明SOCP松弛精确性在整个网络上都有保障而不是只看网损差值。6. 后续扩展方向与我的实操体会6.1 从静态OPF到网络重构MISOCP把SOCP-OPF跑通后很多人紧接着做的事就是网络重构通过控制联络开关的开合状态来降低网损或均衡负荷。重构问题的本质是每条支路要么在运行状态要么在停运状态相当于给支路引入0-1变量z_ij。支路停运时该支路的有功无功潮流必须为0电流平方也要为0支路运行时约束退回原来的DistFlow方程。用大M法可以把这些逻辑写进线性约束得到的模型就是混合整数二阶锥规划MISOCPCplex和Gurobi都能直接解。我个人建议先别急着上大算例先用33节点系统把开关状态变量和联络开关的表示练熟再扩展到大系统。重构问题对SOCP松弛精确性的要求更高因为目标函数里虽然还是网损但0-1变量的存在会把解推向不同的锥面组合容易出现松弛不精确的情况。6.2 从单一目标到多目标与时序优化配电网SOCP-OPF的另一个主流扩展是把单断面静态优化变成多断面时序优化。储能充放电会让节点注入变量带上时间维度DG出力会受光照和风速曲线影响调度目标也从单一时段的网损最小变成一整天的运行成本最小。此时模型规模会线性增长但SOCP的凸性保证了整个问题仍然是凸混合整数规划商业求解器依然有办法求解。我在做时序优化时的一个体会是时间断面数量不要一上来就设成24或96先把典型时段压缩成3到5个断面跑通逻辑再逐步增加断面数否则一旦模型有bug排查多个断面之间的错误会让人崩溃。6.3 规模化应用中的几条实在建议最后说点实战经验吧。第一无论你论文里写SOCP理论有多漂亮工程上一定要有潮流校验这一步。我在说的难听一点没有校验的SOCP结果等于没有兑现的支票。第二从IEEE 33节点开始复现是最快的路径这个系统的数据公开、规模小、全网辐射状非常适合理解模型和调通代码。第三每改一次模型就重新跑一遍校验和保存结果不要等到最后再回归测试否则出了错不好定位到底是模型问题还是数据问题。我在实际使用中发现SOCP在配电网最优潮流里的价值不在于“数学上完美”而在于它把原来不可控的非凸问题变成了可以反复调优的凸问题。你学会了这套建模方式后面无论接分布式电源规划、储能配置还是动态重构所有模型框架都是同一套东西只是往里面加变量和约束而已。如果你正卡在某个SOCP-OPF的复现细节上不妨从这一段代码开始把33节点算例跑通再对照这里的校验与排查思路大概率能把你从“代码不知道错在哪”的泥潭里拉出来。 二阶锥松弛在配电网最优潮流中的原理与Matlab实现
开写之前我想先说明一点这篇文章不是那种“贴一段说明书式的代码然后让你自己琢磨“的凑数笔记。我会把二阶锥松弛在配电网最优潮流里的来龙去脉讲清楚再给出Matlab代码骨架最后把调通过程中容易踩的坑挨个拆开。你现在看到的这套东西是我自己从复现文献模型到跑通IEEE 33节点算例的真实路径希望能帮你少走弯路。为什么配电网最优潮流需要二阶锥松弛1. 配电网最优潮流为什么比输电网的难解1.1 配电网OPF和输电网OPF的“硬”不一样很多刚接触最优潮流OPF的人第一反应是“这不就是潮流加个优化目标嘛”。确实输电网OPF已经有了很成熟的处理方式内点法、原对偶方法二十年前就能解到几千节点规模。但配电网的情况完全两样。输电网X/R比高、电压等级高潮流方程可以近似成直流潮流处理模型差不多是线性的求解器很喜欢这种问题。配电网呢10kV、12.66kV这种等级线路电阻和电抗同量级甚至电阻更大直流潮流简化法在这里基本失效你必须回退到完整的交流潮流模型。交流潮流模型就不好惹了。主动配电网里的节点电压幅值、支路电流幅值、有功无功潮流之间全是非线性耦合关系电压降落的二次项、电流平方项、功率乘积项搅在一起整个可行域是一个非凸集合。目标函数再叠加上分布式电源出力、网损、电压偏差优化问题就变成了一个大规模非凸非线性规划NLP。这类问题对初值极度敏感陷入局部最优的概率不低而且没有一套通用的手段判断你得到的解到底是不是全局最优。1.2 传统NLP方法在配电网模型上的窘境我最初尝试直接用非线性规划的解法器去解配电网OPF用的还是经典的交流潮流方程并把所有等式约束原样写进去。这么做的问题是第一对初值要求高多点起算经常得到不同的结果你根本不知道哪一个是全局最优第二当配电网中有大量分布式电源、功率双向流动时雅可比矩阵条件数变差迭代步长很容易震荡甚至直接发散第三计算时间不可控。更麻烦的是配电网规划问题里常常还带着0-1变量比如网络重构中的开关状态、分布式电源的选址定容。一旦把非凸潮流约束和整数变量放在一起就成了一个非凸混合整数非线性规划MINLP这个问题连当前最前沿的商业求解器都不敢拍胸脯保证收敛。这也是为什么大家都想找一种把非凸约束转化成凸约束的办法——凸化之后哪怕带整数变量退化成MISOCPGurobi和Cplex也能稳定求解。二阶锥松弛Second-Order Cone ProgrammingSOCP就是这类凸化方法里应用最广的一个。1.3 二阶锥松弛的直观理解把“必须压在线上”变成“允许落在锥内”我学习SOCP时觉得最关键的思维转变是不要把它当成一个“近似技巧”而要理解它“松”掉了什么。原始潮流方程里有一个核心等式支路电流幅值的平方乘以首端电压幅值的平方等于该支路有功平方加无功平方。这个等式就是那条“必须压在线上”的曲线约束物理上强制功率、电压、电流三者完全自洽。二阶锥松弛做的事很简单把等式“”改成一个不等式“≥”也就是允许电流电压的组合落在锥形区域内部而不是必须落在锥面上。你可能会问这不是把约束放宽了吗得到的解还是原问题的解吗这就牵出了一个更微妙的性质在目标函数单调地拉大或压缩某个变量时最优解会自动回到锥面上让松弛后的最优解与原问题严格一致。这个性质才是SOCP能大规模落地的前提后面我会重点讲它成立的条件和失效的场景。2. 从DistFlow方程到锥约束的数学推导2.1 先约定记号每条支路的状态量要推导先得统一符号。考虑一个辐射状配电网节点集合为N支路集合为E根节点变电站出口编号为1。i、j支路(i,j)的首端和末端节点r_ij、x_ij支路电阻和电抗P_ij、Q_ij从节点i流向节点j的有功和无功潮流V_i节点i的电压幅值v_i V_i^2节点电压幅值的平方I_ij支路电流幅值l_ij I_ij^2支路电流幅值的平方P_inj,i、Q_inj,i节点i的有功、无功注入分布式电源或变电站注入减去负荷。这个记号和Baran和Wu在1989年提出的DistFlow模型保持一致配电网SOCP相关文献基本都沿用这套写法。符号统一的好处是你在读论文时不会因为流量方向定义不同而绕晕。2.2 DistFlow方程组的三个核心关系辐射状配电网中每条支路(i,j)满足三个方程。第一个是电压降落方程v_j v_i - 2(r_ij P_ij x_ij Q_ij) (r_ij^2 x_ij^2) l_ij这个方程的来源是欧姆定律和相量运算展开之后取幅值平方再把高次项整理成电流平方的形式。第二个是节点功率平衡方程说的是流入节点的功率加节点注入等于流出节点的功率加节点负荷sum_{k∈child(j)} P_jk P_inj,j - P_load,j (P_ij - r_ij l_ij)注意这里P_ij是上游支路流入j的功率扣掉线路损耗r_ij l_ij之后剩余部分才能参与本地平衡。第三个方程描述功率、电压和电流的关系l_ij (P_ij^2 Q_ij^2) / v_i问题就出在这第三个方程上。它是一个分数非线性等式而且关于v_i和l_ij是双线性的写成等价形式是l_ij v_i P_ij^2 Q_ij^2显然是非凸的。整个OPF模型之所以难解根源就在这个等式。2.3 把非凸等式改写成标准二阶锥约束SOCP的标准做法分两步。第一步做变量替换把v_i和l_ij直接设为优化变量而不是再写V_i^2和I_ij^2。这样电压降落方程和功率平衡方程中出现的都是线性项或简单二次项只有l_ij v_i P_ij^2 Q_ij^2这个关联等式是非凸的。第二步是松弛。把等式变成不等式l_ij v_i ≥ P_ij^2 Q_ij^2为什么这个不等式能改写成二阶锥约束把两边展开一下。定义向量x [2P_ij; 2Q_ij; l_ij - v_i]然后构造2-范数约束|| [2P_ij; 2Q_ij; l_ij - v_i] ||_2 ≤ l_ij v_i两边平方4P_ij^2 4Q_ij^2 (l_ij - v_i)^2 ≤ (l_ij v_i)^2整理后得到4(P_ij^2 Q_ij^2) ≤ 4l_ij v_i正好就是l_ij v_i ≥ P_ij^2 Q_ij^2。一个看似复杂的非线性等式就这样转换成了Yalmip里可以直接写给求解器的cone()约束。2.4 松弛为什么能在最优解处自动收紧理论上说松弛后的可行域比原问题大最优值不会比原问题差但可能比原问题“好得离谱”——因为原问题不可行的点松弛后可能可行。如果松弛后的最优解在锥的内部那它就不是原问题的可行解得到的结果没有物理意义。幸运的是在辐射状配电网中SOCP松弛经常是精确的也就是最优解会落在锥面上。这里头有个物理直觉。目标函数如果是网损最小就可以写成 sum(r_ij l_ij)由于电阻r_ij是正数目标函数会尽可能把每个l_ij往下压。而锥约束允许的范围是l_ij ≥ (P_ij^2 Q_ij^2)/v_i能压多低压多低。于是每个支路最后都会被压到锥的边界上松弛精确性天然成立。同理目标函数如果包含正系数倍的支路电流平方项也能达到类似效果。目前文献里比较公认的结论是当网络是辐射状、目标函数对电流/网损有单调递增趋势、并且不存在过紧的电流上限约束时SOCP松弛基本是精确的。但这不是免费的午餐后面我会讲怎么验证、失效了怎么办。3. Matlab代码实现从33节点数据到SOCP求解3.1 数据准备与基准值单位先说明一点我不建议依赖Matpower自带的case文件来跑SOCP-OPF即便Matpower里有配电网形式的数据格式依然围绕输电网习惯设计对你的模型构建没多大帮助。更直接的做法是手动把IEEE 33节点系统的数据整理成两个矩阵。第一个矩阵是bus_data两列分别是有功负荷PdkW和无功负荷Qdkvar。第二个矩阵是branch_data四列分别是首端节点、末端节点、线路电阻rohm、线路电抗xohm。IEEE 33节点标准数据在论文附录和很多公开资源里都能找到直接把那两张表粘贴到矩阵里即可。下面给一个片段格式% bus_data: [节点编号, Pd(kW), Qd(kvar)] % 1号节点是根节点负荷为0 bus_data [ 1, 0, 0; 2, 100, 60; 3, 90, 40; 4, 120, 80; 5, 60, 30; ... ]; % branch_data: [首端节点, 末端节点, r(ohm), x(ohm)] branch_data [ 1, 2, 0.0922, 0.0470; 2, 3, 0.4930, 0.2511; 3, 4, 0.3660, 0.1864; ... ];很多人在这第一步就翻车并不是拿不到数据而是单位搞错了。配电网你直接用Ω和kvar去建模数值尺度会差好几个量级Cplex和Gurobi的数值引擎会非常不稳定。所以必须标幺化。我习惯取基准电压Vbase 12.66 kV基准功率Sbase 10 MVA。这样基准阻抗Zbase Vbase^2 / Sbase 16.02 Ω负荷和注入也都除以Sbase换算成标幺值。3.2 Yalmip模型变量、目标函数、约束环境是Matlab配合Yalmip工具箱求解器可以选Cplex、Gurobi或者Mosek。Yalmip的作用是把优化问题建模成求解器能懂的格式尤其对SOCP约束它提供了直观的cone函数。下面是变量定义和约束构建的核心代码nb size(bus_data, 1); % 节点数 nl size(branch_data, 1); % 支路数 Vbase 12.66e3; % V Sbase 10e6; % VA Zbase Vbase^2 / Sbase; % ohm r_pu branch_data(:,3) / Zbase; x_pu branch_data(:,4) / Zbase; Pd_pu bus_data(:,2) / Sbase; Qd_pu bus_data(:,3) / Sbase; % 变量定义 v sdpvar(nb, 1); % 节点电压幅值平方 I2 sdpvar(nl, 1); % 支路电流幅值平方 P sdpvar(nl, 1); % 支路首端有功 Q sdpvar(nl, 1); % 支路首端无功 Pg sdpvar(nb, 1); % 节点有功注入 Qg sdpvar(nb, 1); % 节点无功注入 C []; % 支路约束电压降落方程 二阶锥松弛 for k 1:nl i branch_data(k,1); j branch_data(k,2); C [C, v(j) v(i) - 2*(r_pu(k)*P(k) x_pu(k)*Q(k)) ... (r_pu(k)^2 x_pu(k)^2) * I2(k)]; C [C, cone([2*P(k); 2*Q(k); I2(k) - v(i)], I2(k) v(i))]; end % 节点功率平衡 for j 1:nb inbound find(branch_data(:,2) j); outbound find(branch_data(:,1) j); C [C, -sum(P(inbound) - r_pu(inbound).*I2(inbound)) ... sum(P(outbound)) Pg(j) - Pd_pu(j)]; C [C, -sum(Q(inbound) - x_pu(inbound).*I2(inbound)) ... sum(Q(outbound)) Qg(j) - Qd_pu(j)]; end % 根节点电压固定变电站注入自由 C [C, v(1) 1.0]; C [C, Pg(1) 0]; C [C, Qg(1) 0]; % 无DG节点注入为0 C [C, Pg(2:end) 0]; C [C, Qg(2:end) 0]; % 电压上下限 0.95~1.05 p.u. C [C, 0.95^2 v 1.05^2]; % 支路电流上限按实际导线载流设定 Imax_pu 1.2; C [C, I2 Imax_pu^2]; % 目标函数网损最小 objective sum(r_pu .* I2);这里有几个细节我要特别强调。第一功率平衡方程里inbound变量是支路末端等于当前节点j的那些支路也就是从上级注入j的支路outbound则是从j流向子节点的支路。由于配电网是辐射状任意非根节点上方只有一条支路所以inbound多数时候只有一个元素但代码写成通用形式没有坏处。第二节点注入的正方向是“注入网络”所以根节点Pg(1)为正值表示变电站向网络输送功率如果未来考虑分布式光伏大发导致根节点倒送功率可以把约束改成Pg(1)的上下限。第三样例子中我把2号节点及以后所有节点的Pg设为0意思是不考虑分布式电源如果你要加DG就把该节点的Pg改为带上下限的变量而不是固定为0。3.3 求解器设置与结果恢复模型构建完后用下面这段代码调用求解器ops sdpsettings(solver, cplex, verbose, 2, cplex.mip.tolerances.mipgap, 1e-4); optimize(C, objective, ops); % 结果恢复 V_opt sqrt(value(v)); P_opt value(P); Q_opt value(Q); I2_opt value(I2); loss_opt sum(r_pu .* value(I2)) * Sbase; % 转成实际功率单位W % 电压分布图 figure; plot(V_opt, -o, LineWidth, 1.5); xlabel(节点编号); ylabel(电压幅值 (p.u.)); grid on;运行成功后你会看到33节点系统在网损最小化目标下末端节点电压比根节点低支路电流平方值也落在合理范围。用Cplex求解通常不到一秒就能收敛这比直接解NLP体感快得多。3.4 代码里最容易写错的三处第一Yalmip里的变量名尽量避免用字母l它和数字1长得太像调试时很难分清。上面代码统一用I2表示支路电流幅值平方就是不想给自己找麻烦。第二cone的写法。很多人第一次写SOCP约束时会把它写成norm([2P;2Q;I2-v]) I2v这在Yalmip里也支持但cone函数的签名是cone(x, y)表示norm(x)y别搞反了参数顺序。第一参数必须是向量第二参数必须是标量否则Yalmip会直接报错。第三功率平衡方向的符号。我在代码里明确用“注入网络为正方”的约定这是DistFlow文献的标准约定。如果你从网上找到一段代码它的符号可能用的是“负荷为正方”最后结果P和Q的值符号会整体相反但你如果不检查支路潮流的物理方向很可能意识不到这个问题。4. 松弛不精确怎么办验证与补救4.1 松弛是否精确的判断标准SOCP方法最大的隐患就是松弛不精确但你可以通过一个简单办法验证把SOCP求得的节点注入和根节点电压作为输入重新做一次精确的交流潮流计算得到网损然后和SOCP目标函数值对比。如果两者差值在0.1%以内说明松弛精确可以放心用如果差距过大说明最优解落在锥内部SOCP结果不可直接用。我通常在代码里加这么一段校验逻辑% flow_loss 是从潮流计算得到的网损 % loss_opt 是SOCP目标函数值 diff_loss abs(flow_loss - loss_opt); if diff_loss / flow_loss 1e-4 disp(SOCP松弛精确); else disp(SOCP松弛不精确需要检查约束和边界); end如果你不想自己写潮流计算可以设定好潮流数据格式后调用Matpower的runpf把SOCP得到的节点注入做一次潮流比较网损。个人经验是至少对每一篇论文的算例都做一次这样的校验不要上来就默认SOCP一定精确。4.2 哪些条件下SOCP会给出错误结果我复盘自己踩过的坑和文献里反复强调的案例SOCP松弛失效主要集中在几类场景。第一目标函数不“压”电流。比如目标只是最小化变电站购电成本而购电成本又只跟根节点有功有关跟支路电流没有直接单调关系松弛后的解可能停留在锥内部得到的网损偏低电压整体偏“美观”。这时候你就要小心了。第二电流上限约束太紧。如果某条支路的电流上限恰好卡在最优点上导致目标函数再怎么压也压不到锥面松弛就可能失效。辐条状网络中支路电流上限是SOCP精确性的边界条件。第三网络不是严格的辐射状。配电网实际运行时可能有联络开关闭合形成环网。在环网情况下SOCP松弛精确性不再有理论保证解出来的结果可能完全不符合潮流方程。这是很多人看SOCP文献时忽略的一个大前提。第四节点负荷过低。负荷接近零时潮流方程本身接近病态锥约束可能被轻易满足不需要贴着锥面走。4.3 可行的修复手段一旦发现松弛不精确有几种补救思路可以尝试按从简单到复杂的顺序说。第一种是加罚项。在目标函数里加一个很小的正数乘以所有支路I2之和也就是min sum(r_ij l_ij) ε*sum(l_ij)这样目标函数被迫继续“压”电流把解往锥面上推。ε一般取1e-4到1e-2之间但要注意罚项别大到扭曲原目标。第二种是收紧边界。检查一下是不是某些约束边界设置过宽比如支路电流上限给的太大给了锥约束“偷懒”的空间。适当收紧电流上限或者把节点电压上下限从0.95~1.05改到0.98~1.02往往能让松弛重新精确。第三种是交替迭代。解完SOCP后把得到的节点电压代入原始潮流方程算出精确的支路电流和网损然后把这些值作为新的初始点再解一次SOCP。这种序贯SOCP在工程中很常见虽然不是严格全局最优但实际效果稳定。第四种是切换到SDP或者更精确的凸松弛。当SOCP失效时半定松弛SDP在某些场景下更紧代价是计算量大不少。如果网络规模不大可以试试SDP作为对照。5. 调通代码后一定会遇到的几类实战问题5.1 求解器选型Cplex、Gurobi、Mosek、Sedumi怎么选理论上SOCP是凸优化任何支持二阶锥的商业求解器都能解。但我实际测试下来同样一个33节点网损最小问题不同求解器的求解时间和对Yalmip模型格式的容错度还是有差异。求解器求解速度大规模表现与Yalmip兼容性备注Cplex快优秀好配电网MISOCP场景常用Gurobi快优秀好参数多收敛行为稳定Mosek很快优秀好SOCP内点法非常专业Sedumi较慢一般可免费适合小规模教学我个人的建议是学习阶段用Cplex或Gurobi正式科研时再看看Mosek。配电网SOCP-OPF的问题规模一般不大33节点、123节点Cplex就足够用了如果做上百个时间断面的时序优化Mosek在SOCP内点法上的性能优势会更明显。5.2 量纲和数值尺度引起的玄学报错Yalmip本身能处理一部分量纲混乱问题但求解器内部对数值尺度是有偏好的。有次我把负荷数据直接写成kW注入写成标幺结果Cplex报infeasible我检查了很久才发现两条数据一个用kW一个用标幺。建议一开始就把所有电气量统一成标幺值不要嫌麻烦。根节点电压固定为1.0电压上下限写0.95^2~1.05^2支路电阻电抗都除以Zbase。这样整个模型的数值量级都在0.01到10之间对求解器非常友好。另一个常见问题是支路电流上限设置得过松。如果你不设置I2的上限模型可能会解出一个电流特别大的荒唐结果因为网损最小目标里电阻乘电流平方虽然会限制电流但某些边界场景下电压约束压力更大电流可能被顶到非物理上限。给每条支路一个合理的载流量上限既是物理约束也帮助数值稳定。5.3 求解时间长与迭代慢的原因对大系统SOCP虽然凸但求解时间不总是秒级。如果发现求解很慢第一步检查是不是Yalmip把SOCP误判成了非线性规划。Yalmip会检查约束类型但有时候前期建模不规范比如把原本应该是线性约束的电压降落方程写成了间接形式会导致求解器走NLP路径。你可以通过调用ops.solver查看实际选中的求解器如果显示fmincon之类大概率是被Yalmip当NLP处理了。第二步看模型是否有大量冗余变量。比如节点功率平衡用循环生成时每个循环都重新做一次find(branch_data(:,2)j)模型规模小时无所谓上千节点时积累的find耗时不可忽略。可以把inbound/outbound在循环外提前算好存成元胞数组再进循环。5.4 结果后处理与可视化经验SOCP直接解出来的变量是v和I2一个是电压幅值平方一个是电流幅值平方后处理时一定要记得开根号。我刚开始跑的时候直接拿value(v)去画电压分布图画出来全是0.9到1.1附近的数还觉得奇怪——后来才反应过来这是平方值不是电压本身。可视化除了画电压分布曲线我还会画一个“松弛间隙”的柱状图。对于每条支路计算I2_opt(k)*v_opt(i) - (P_opt(k)^2 Q_opt(k)^2)再除以I2_opt(k)*v_opt(i)。这个值表示该支路离锥面的相对距离所有支路的间隙都接近0说明SOCP松弛精确性在整个网络上都有保障而不是只看网损差值。6. 后续扩展方向与我的实操体会6.1 从静态OPF到网络重构MISOCP把SOCP-OPF跑通后很多人紧接着做的事就是网络重构通过控制联络开关的开合状态来降低网损或均衡负荷。重构问题的本质是每条支路要么在运行状态要么在停运状态相当于给支路引入0-1变量z_ij。支路停运时该支路的有功无功潮流必须为0电流平方也要为0支路运行时约束退回原来的DistFlow方程。用大M法可以把这些逻辑写进线性约束得到的模型就是混合整数二阶锥规划MISOCPCplex和Gurobi都能直接解。我个人建议先别急着上大算例先用33节点系统把开关状态变量和联络开关的表示练熟再扩展到大系统。重构问题对SOCP松弛精确性的要求更高因为目标函数里虽然还是网损但0-1变量的存在会把解推向不同的锥面组合容易出现松弛不精确的情况。6.2 从单一目标到多目标与时序优化配电网SOCP-OPF的另一个主流扩展是把单断面静态优化变成多断面时序优化。储能充放电会让节点注入变量带上时间维度DG出力会受光照和风速曲线影响调度目标也从单一时段的网损最小变成一整天的运行成本最小。此时模型规模会线性增长但SOCP的凸性保证了整个问题仍然是凸混合整数规划商业求解器依然有办法求解。我在做时序优化时的一个体会是时间断面数量不要一上来就设成24或96先把典型时段压缩成3到5个断面跑通逻辑再逐步增加断面数否则一旦模型有bug排查多个断面之间的错误会让人崩溃。6.3 规模化应用中的几条实在建议最后说点实战经验吧。第一无论你论文里写SOCP理论有多漂亮工程上一定要有潮流校验这一步。我在说的难听一点没有校验的SOCP结果等于没有兑现的支票。第二从IEEE 33节点开始复现是最快的路径这个系统的数据公开、规模小、全网辐射状非常适合理解模型和调通代码。第三每改一次模型就重新跑一遍校验和保存结果不要等到最后再回归测试否则出了错不好定位到底是模型问题还是数据问题。我在实际使用中发现SOCP在配电网最优潮流里的价值不在于“数学上完美”而在于它把原来不可控的非凸问题变成了可以反复调优的凸问题。你学会了这套建模方式后面无论接分布式电源规划、储能配置还是动态重构所有模型框架都是同一套东西只是往里面加变量和约束而已。如果你正卡在某个SOCP-OPF的复现细节上不妨从这一段代码开始把33节点算例跑通再对照这里的校验与排查思路大概率能把你从“代码不知道错在哪”的泥潭里拉出来。