做配电网优化的人大概率都会在某个阶段被“无功优化”这个词拽进坑里。配电网的无功功率平衡直接关系到线损和末端电压质量而IEEE 33节点作为最经典的辐射状测试系统几乎成了各种算法验证的标配场地你用二阶锥规划SOCP在它上面做无功优化再用Matlab代码实现一遍基本上是这个领域的“入门仪式”。这篇不打算讲虚的主要把我自己在建模、编码、调参和结果校核过程中觉得最关键的东西整理出来DistFlow方程怎么变成锥约束、数据和基准值怎么处理、YALMIP配合求解器要注意什么以及怎么确认SOCP松弛是紧的、结果是可信的。无论你是刚进这个方向的研究生还是做配网项目需要快速搭一版优化方案的工程师按这条线走能省不少弯路。1. 为什么无功优化要用二阶锥规划三个绕不开的现实问题1.1 配电网无功调节的三类手段与各自的脾气配电网无功优化的核心任务是在满足电压、潮流、设备容量等约束的前提下通过调节各类无功源让系统网损最小、电压质量最好。这个任务之所以复杂是因为可调手段本身就五花八门脾气各不相同。第一类是变电站侧的无功调节设备比如有载调压变压器分接头OLTC和并联电容器组。OLTC通过改变变比来调整电压但分接头是离散档位动作一次有机械磨损通常一天内动作次数受限并联电容器组也是分组投切出力是阶梯式的投切过程伴随冲击电流不能太频繁。第二类是装设在馈线上的无功补偿装置比如SVC、SVG它们能连续、快速地输出无功调节性能好但单台容量有限成本也高。第三类是随着分布式电源普及越来越重要的一类——光伏逆变器和储能变流器的无功能力。逆变器本身可以工作在非单位功率因数状态发出或吸收无功调节粒度细、响应速度快几乎零边际成本但它会影响有功输出能力而且产权归属、调度方式更复杂。这三类设备的调节特性差异很大连续与离散并存快速与慢速并存全局与局部并存。做优化建模时最难受的就是离散变量和连续变量混在一起。OLTC档位和电容器投切组数是整数变量逆变器无功和SVC出力是连续变量这在数学上天然构成混合整数问题。叠加配电网的非线性潮流方程问题就变成混合整数非线性规划MINLP求解非常困难。这也是我从一开始就坚持用二阶锥规划的原因它能把非凸的潮流约束转化成凸约束以可接受的计算代价换来全局最优解的保证。1.2 传统求解方法在配电网上为什么吃力很多刚接触配电网无功优化的人第一反应是用粒子群、遗传算法这类启发式算法。好处是代码好写、概念好懂不需要处理复杂的导数信息但坏处也很明显没有全局最优性保证一组参数换一组参数跑出来的结果差异很大有时候优化半天出来的解还不如工程经验拍脑袋定下来的方案。还有一条路线是用内点法解非线性规划NLP把离散变量先松弛成连续变量求解后再取整或修补。问题在于配电网潮流方程是强非线性的目标函数又常常非凸内点法对初值非常敏感收敛到局部最优是常态。更麻烦的是NLP松弛加取整得到的解往往根本不满足设备档位约束需要额外设计启发式修补策略逻辑一下子变得很脆弱。还有直接调用商业求解器求解MINLP的尝试。MINLP全局求解器的效率和规模限制一直是瓶颈在几十个节点的配电网模型上还能凑合节点数一多、决策变量一多求解时间就非常难看了。所以从工程落地和学术研究的角度看把非凸问题转化成凸问题是更聪明的路。凸优化有一个非常关键的好处局部最优就是全局最优。只要能把配电网潮流约束合理改写成凸约束后续无论用哪个成熟求解器结果都可复现、可解释、可审查。1.3 二阶锥把“最难啃的约束”变成“好解的凸约束”二阶锥规划SOCP是凸优化的一种具体形式。它的标准约束长这样对一个向量要求它落在某个锥体区域内。可以把它直观理解成一个冰激凌筒筒壁的半径决定能装多少所有点必须落在筒壁以内。数学上写出来是sqrt(y1² y2²) ≤ x其中 x 是圆锥的纵轴y1、y2 是横截面上的两个坐标。这个约束看起来还是非线性的但它属于凸约束可以用专门的内点算法在多项式时间内高效求解而且数值稳定性好。配电网潮流问题里有大量二次项比如电压与电流的乘积、有功无功的平方和这些项让问题非凸。SOCP的妙处在于通过变量替换和不等式替换恰好可以把这些非凸二次约束转化为锥约束。后面我会详细推导这个过程。在实际使用中它和原始非线性问题的解非常接近——在网损最小化这类目标下凸松弛后求得的解几乎总是能回到原始问题的可行域。这就是SOCP能在配电网无功优化中站稳脚跟的原因。2. 建模推导从DistFlow方程到二阶锥约束的完整过程2.1 DistFlow支路潮流方程辐射状网络的“能量守恒”配电网通常是辐射状网络也就是一棵树每个节点有唯一的供电路径。针对这种拓扑结构支路潮流方程DistFlow是比完整交流潮流方程更精炼的模型。它最大的优点是不需要计算相角和复数运算直接把有功、无功、电压幅值的关系写出来。考虑一条支路 i - ji 是首端j 是末端。定义 r_ij、x_ij 为支路电阻、电抗P_ij、Q_ij 为从 i 流向 j 的有功、无功功率I_ij 为支路电流幅值V_i、V_j 为节点 i、j 的电压幅值。DistFlow方程由三组式子组成。第一组是节点功率平衡流过节点 j 的所有出线功率加上节点 j 的负荷与注入之差等于节点 j 从上游支路得到的净功率。第二组是支路电压降落关系末端电压的平方等于首端电压平方减去两条修正项一条对应线路上的功率损耗另一条对应电流产生的压降。第三组是支路电流定义电流平方与首端电压平方的乘积等于流过的视在功率平方。用 U_i 表示 V_i²用 l_ij 表示 I_ij²这三组式子可以写成P_j sum(k in N(j)) P_jk r_jk * l_jk P_load_j - P_gen_jQ_j sum(k in N(j)) Q_jk x_jk * l_jk Q_load_j - Q_gen_jU_j U_i - 2 * (r_ij * P_ij x_ij * Q_ij) (r_ij² x_ij²) * l_ijl_ij * U_i P_ij² Q_ij²我比较喜欢这个模型的原因是变量的物理意义非常清楚U 是电压幅值平方l 是电流幅值平方P、Q 是支路潮流的实部和虚部。比起面向输电网的完整交流潮流这个模型更适合配电网因为它的数据需求简单只需线路阻抗和节点负荷公式也简洁。2.2 非线性项是怎么被“松弛”进锥里的问题出在最后一组等式l_ij * U_i P_ij² Q_ij²。这是一个非凸约束因为左边是两个待求变量的乘积右边是平方和整个约束定义了一个双曲抛物面不是凸集。如果直接把这个约束丢给求解器问题就是MINLP或NLP难解。SOCP的思路是把等式放松成不等式l_ij * U_i ≥ P_ij² Q_ij²这个不等式的物理含义是流过支路的视在功率平方不超过电流平方乘以首端电压平方。从数学上说这个不等式定义了一个二次锥约束而它恰好是凸的。但直接这样写还不够求解器通常需要标准形式。通过变量变换可以把它化成真正的二阶锥标准形式。定义辅助向量要求它的2-范数满足特定关系|| [2 * P_ij; 2 * Q_ij; l_ij - U_i] || ≤ l_ij U_i这个范数不等式等价于 l_ij * U_i ≥ P_ij² Q_ij²。验证方法很简单两边平方展开就能看到它可以化为 (l_ij U_i)² ≥ (2P_ij)² (2Q_ij)² (l_ij - U_i)²化简后正是 4 * l_ij * U_i ≥ 4 * (P_ij² Q_ij²)。在YALMIP里这一条可以直接写成cone([2*P(k); 2*Q(k); l(k)-U(i)], l(k)U(i))。这里有个关键问题把等式放松成不等式凭什么不会导致错误解原因在目标函数。无功优化最常用的目标是最小化网损min sum(r_ij * l_ij)。网损正比于电流平方为了让目标最小求解器会倾向于让 l_ij 尽量小。如果某个电压水平下还能把电流再压小一点而不违反其它约束目标函数就会进一步下降。也就是说在网损最小这种“单调递减”的目标下最优解会自动把松弛间隙压到零不等式最终变成等式。这是所谓“松弛精确”exact relaxation的性质。当然这不是免费的午餐。松弛是否精确取决于目标函数形式、网络拓扑和负荷条件。如果把目标改成电压偏差最小化或者加入某些特殊设备约束松弛可能不再严格精确最后得到的结果对应到原始潮流上可能有偏差。所以在实际项目中我从来不会跳过求解后的校核单独计算每个支路的锥间隙也就是 l_ij * U_i 与 P_ij² Q_ij² 的差值如果最大间隙在1e-5量级以下才认为结果可信。2.3 完整数学模型目标函数与约束清单把上面这些拼在一起一个基于SOCP的配电网无功优化模型可以这样写。决策变量包括节点电压幅值平方 U_i、支路电流幅值平方 l_ij、支路有功 P_ij 和无功 Q_ij以及各类无功设备的出力 Qc。如果考虑OLTC和电容器还需要离散变量。目标函数一般写成三项加权min f w1 * sum(r_ij * l_ij) w2 * sum((U_i - U_ref)²) w3 * sum(Qc²)第一项是系统网损第二项是电压偏差第三项是设备动作成本或无功裕度惩罚。权重选择要结合工程需求网损通常权重最大电压偏差次之设备成本最小。纯学术验证场景下很多论文只保留网损项也可以但实际工程里很少只有一个目标。约束条件分为四类。第一类是DistFlow功率平衡约束和锥约束这是核心。第二类是节点电压约束用平方形式写就是 U_min ≤ U_i ≤ U_max标幺值下一般取 0.95² 和 1.05²。第三类是支路电流和容量约束l_ij ≤ l_ij_max对应导线载流量限制。第四类是无功源出力约束比如SVC输出有上下限电容器投切组数是整数区间逆变器无功受容量和功率因数限制。难点在于把控制变量和潮流变量耦合起来。比如节点 j 接一台SVC节点功率平衡方程里流入该点的无功要从 Q_load_j 变成 Q_load_j - Q_svc_j同时追加 0 ≤ Q_svc_j ≤ Q_max。如果是OLTC一般用虚拟节点加理想变压器模型会牵涉到对锥约束的修正代码复杂度明显上升。所以初期建议先做纯辐射状、不含OLTC的连续无功优化框架跑通后再逐步加入离散设备。3. IEEE33节点测试系统数据准备里最容易出错的三件事3.1 节点负荷与支路参数速查IEEE 33节点系统是目前配电网优化文献里出现频率最高的标准测试系统。它由一个12.66kV的平衡节点节点1供电包含32条普通支路和5条联络支路。基础状态下联络支路断开网络保持辐射状。系统总负荷为3715 kW加2300 kvar负荷作为常数功率处理这是无功优化最经典的假设。做Matlab实现时我强烈建议把支路参数和负荷数据单独整理成两个矩阵保存不要散落在代码里。支路矩阵每一行对应一条支路列分别是首节点、末节点、支路电阻欧姆、支路电抗欧姆负荷矩阵每一行是节点编号、有功负荷kW、无功负荷kvar。网上不同版本的数据格式不太统一有的用标幺值有的用有名值有的把联络支路顺序排在最前面。拿到任何一份数据第一步不是写代码而是先把数据校核一遍确认支路总数、节点总数、总负荷是否和标准值一致。比如我常用的部分数据长这样节点P负荷(kW)Q负荷(kvar)210060390408200100189040334530首端末端R(Ω)X(Ω)120.09220.0470230.49300.2511340.36600.186418190.93200.473032330.34100.2030先肉眼对一遍总负荷再开始建模能挡掉一大半的数据坑。3.2 标幺化用错基准值后面全废标幺化是无功优化编码里最常见的错误来源。它不产生算法层面的大错但会产生“结果全部跑飞”的后果而且很难排查。IEEE 33节点系统常用的基准值是基准功率 10 MVA基准电压 12.66 kV。阻抗基值是基准电压平方除以基准功率约等于 16.02 欧姆。线路电阻和电抗的有名值除以这个数得到标幺阻抗。负荷功率直接除以 10 MVA例如节点8负荷 200 kW 的话标幺值就是 0.02。为什么推荐用标幺值因为锥约束里的数值尺度如果不统一会出现量级差异巨大的问题。有名值下电压平方是几万伏特平方量级功率是几千瓦量级混在一个约束里求解器数值条件数很差收敛慢甚至发散。标幺化后所有变量都在0到1附近数值稳定性好很多。有些同学某个数据版本里用了有名值另一处又用了标幺值最后网损结果差好几个量级怎么调参都救不回来。3.3 节点-支路关联关系与网络拓扑的坑第三件容易出问题的事是拓扑表示。辐射状网络是一棵树节点1是根。我习惯把支路按首端到末端的顺序排列用两个数组 bus_i 和 bus_j 分别记录支路两端节点。所有DistFlow方程按支路遍历顺序从根附近往外走做功率平衡时可以方便地累加子节点的净注入。有一个隐蔽的坑如果支路顺序是任意给出的那么遍历到某条支路时它的首端电压 U_i 可能还没有被上游支路方程定义到“确定”状态。YALMIP 是直接把所有约束一次性交给求解器的不存在顺序执行问题但后续做敏感性分析时拓扑数组的顺序会影响可读性和调试效率。建议从根出发做一次BFS或DFS排序再按排序结果组织支路参数和变量索引。另外要注意平衡节点电压的处理。节点1是无穷大母线电压固定在1.0标幺值也就是 U_1 固定为1.0。很多入门代码忘了给根节点电压赋值导致解出来的电压整体偏移网损异常。初始化时给 U 都赋初值1.0然后单独固定 U(1)1.0这是个简单但非常管用的措施。4. MatlabYALMIP实现从变量定义到松弛校核的闭环4.1 工具箱与求解器选型参考Matlab 环境下做 SOCP最常用的建模工具箱是 YALMIP。它负责把变量、约束、目标翻译成求解器能识别的标准形式然后调度底层求解器。底层求解器选型直接影响可解性和速度。根据实际使用经验可以按下面这个表选求解器支持SOCP支持MISOCP授权情况适用场景MOSEK是是商业许可证大规模、稳定性要求高Gurobi是是商业许可证学术许可较容易申请工程与科研通用CPLEX是是商业许可证老牌部署ECOS是部分开源中小规模、教学SeDuMi是不支持开源纯SOCP科研SDPT3是不支持开源纯SOCP科研如果课题只需要纯连续变量 SOCP用开源的 ECOS 或 SeDuMi 就够了不必折腾商业许可证。如果要把 OLTC 和电容器整数变量加进去形成 MISOCP基本建议直接用 Gurobi 或 MOSEK求解效率和稳定性明显更好。我第一次做含离散变量的版本时用开源求解器硬跑等了一个小时没结果换 Gurobi 后几十秒就出来了完全不是一个量级。4.2 核心代码框架可复用下面给出一个我常用的代码框架覆盖从数据读取到结果校核的完整流程。为了篇幅可控省略完整数据数组只保留逻辑主线。% 1. 数据输入标幺值 baseMVA 10; baseKV 12.66; Zbase baseKV^2 / baseMVA; % branch: [bus_i, bus_j, R(ohm), X(ohm)] branch [...]; bus_load [...]; % [node, P(kW), Q(kvar)] Rpu branch(:,3) / Zbase; Xpu branch(:,4) / Zbase; Pload_pu bus_load(:,2) / (baseMVA*1000); Qload_pu bus_load(:,3) / (baseMVA*1000); % 2. 变量定义 nB 33; nL size(branch,1); U sdpvar(nB,1); % 节点电压幅值平方 l sdpvar(nL,1); % 支路电流幅值平方 P sdpvar(nL,1); % 支路有功 Q sdpvar(nL,1); % 支路无功 Qc sdpvar(nB,1); % 无功补偿装置出力 Pgen zeros(nB,1); % 分布式电源有功注入若不计可全置零 % 3. 约束 Constraints []; Constraints [Constraints, U(1) 1.0]; % 平衡节点电压 Constraints [Constraints, 0.95^2 U 1.05^2]; % 电压上下限 Constraints [Constraints, l 0]; % 无功源出力范围示例节点18、节点33各接一组补偿 Constraints [Constraints, 0 Qc(18) 0.3, 0 Qc(33) 0.3]; for k 1:nL i branch(k,1); j branch(k,2); % 功率平衡收集节点j的所有出线 Children find(branch(:,1) j); sumP 0; sumQ 0; for m 1:length(Children) km Children(m); sumP sumP P(km) Rpu(km)*l(km); sumQ sumQ Q(km) Xpu(km)*l(km); end Constraints [Constraints, P(k) sumP Pload_pu(j) - Pgen(j)]; Constraints [Constraints, Q(k) sumQ Qload_pu(j) - Qc(j)]; % 电压降落 Constraints [Constraints, U(j) U(i) - 2*(Rpu(k)*P(k) Xpu(k)*Q(k)) (Rpu(k)^2 Xpu(k)^2)*l(k)]; % 二阶锥松弛 Constraints [Constraints, cone([2*P(k); 2*Q(k); l(k)-U(i)], l(k)U(i))]; end % 4. 目标与求解 objective sum(Rpu .* l) 0.1*sum((U - 1).^2); ops sdpsettings(solver,gurobi,verbose,1); sol optimize(Constraints, objective, ops); % 5. 结果校核 if sol.problem 0 U_opt value(U); l_opt value(l); P_opt value(P); Q_opt value(Q); gap max(l_opt .* U_opt(branch(:,1)) - (P_opt.^2 Q_opt.^2)); fprintf(最大锥间隙: %.3e\n, gap); else disp(求解失败请检查数据或求解器配置); end这个框架有几个地方值得说明。第一节点功率平衡方程中由于网络辐射状每个节点最多一个父支路和多个子支路所以我用find(branch(:,1) j)收集节点 j 的所有出线在节点多时比手工排列清晰很多。第二平衡节点节点1的功率不需要用方程约束它的功率是自由变量用来平衡全网。第三目标函数里加了 0.1 倍的电压偏差项这是为了防止纯网损目标下出现电压偏高的问题也可以直接删掉看具体需求。4.3 结果校核锥间隙不检查等于白算很多教程讲完求解就结束但我会专门留一个校核环节。SOCP 松弛后得到的最优解对应到原始 DistFlow 方程可能有偏差。比如求解器为了数值便利可能让某个支路的锥间隙大到 1e-2这时候直接把电压和网损值拿出来写报告是不严谨的。校核分三步走。第一步计算最大锥间隙 gap即上面代码的最后一行理论上应趋近0我一般要求小于 1e-5 量级。第二步把求解得到的 P、Q、l 回代到原始 DistFlow 方程检查每个节点的功率不平衡量。第三步如果存在离散变量还需检验取整后偏差是否可接受必要时重新求解一次带离散变量修正的版本。如果锥间隙过大处理办法一般有两类。一是检查数据标幺化是否统一重点看阻抗和功率的量纲。二是调整目标函数权重比如把电压偏差项的权重大幅度调小因为偏差项的存在可能让松弛不精确。我遇到过的情况中绝大多数异常都出在这两个原因上。5. 从33节点走向实际工程结果解读与功能扩展5.1 怎么判断优化效果网损、电压、间隙三个指标拿到优化结果后第一步不是急着截图而是和不优化的基准场景对比。基准场景就是所有无功补偿为0、所有离散设备保持在初始状态下的潮流结果。对比时重点看三个指标网损下降百分比、最低电压改善幅度、最大锥间隙。网损指标的标准做法是先求基准潮流下的总网损再求优化后的总网损计算相对下降比例。在 IEEE 33 节点系统上如果允许在合适位置投入足够无功补偿网损下降 20% 以上很常见最低电压通常也会从 0.92 附近回升到 0.97 左右。如果优化后网损几乎没下降先怀疑补偿容量上限给得太小或者补偿位置落在低灵敏度节点上。场景总网损(kW)最低电压(pu)最大锥间隙基准无补偿202.70.9130连续无功优化151.20.9483.1e-7含离散设备优化158.60.9424.5e-7电压指标要注意看最薄弱节点的提升而不只是平均电压。配电网无功优化的一个主要收益就是改善末端电压画出优化前后的电压分布曲线一目了然。第三个指标是锥间隙前面反复强调过这里不再重复。再分享一个分析经验做灵敏度分析。固定网络结构把 SVC 容量从0逐步增大到上限记录每次优化后的网损和最低电压会看到一条先快速下降、后趋于饱和的曲线。这个饱和点其实就是系统最经济的无功补偿规模比你直接跑单点优化更有工程指导意义。5.2 高渗透分布式电源与多时段耦合的扩展思路IEEE 33节点只是第一步。真实配电网里分布式光伏和储能越来越多无功优化的建模边界也随之改变。光伏逆变器带来的主要变化是发电具有波动性出力的不确定性会影响无功调节效果。常用处理方式是把一个调度周期划分成多个时段每个时段的光伏出力用预测值在约束里给每个时段建一套潮流模型时段之间再通过储能电量、设备动作限制等耦合起来。这其实是从单时段 SOCP 变成多时段 MISOCP变量规模成倍增长但结构不变Gurobi 这类求解器仍能处理。储能接入后节点功率平衡里多了充放电功率变量同时多了储能电量的动态约束下一时段电量等于当前电量乘以自损耗系数加上或减去充放电功率乘以时间步长。储能作为“时间转移器”让无功优化具备了跨时段的能量调节能力这是单时段模型完全做不到的。还有一个扩展是 OLTC 和电容器组。把分接头档位建模为整数变量在目标里加入动作次数惩罚可以避免优化结果要求设备频繁动作。注意OLTC 在 DistFlow 框架里不是简单改一下电压约束就能处理通常需要增加一个理想变压器支路把变比变量引入电压降落方程数学上会复杂一些建议在基础版本完全跑通后再做。5.3 关于鲁棒性与三相不平衡的提醒最后说一个容易踩的坑IEEE 33节点是三相平衡模型所有方程都建立在单相正序假设上。真实低压配电网很多是三相四线制单相负荷比重高三相不平衡现象非常普遍。如果直接把单相 SOCP 模型套用到实际馈线上优化结果可能在现场严重失真尤其是零序和负序分量占比较大的台区。可行的做法有两种一是采用三相 DistFlow 模型把每个节点扩展成 A、B、C 三相支路参数用相间互阻抗矩阵约束数量直接乘三计算量增大但模型更真实二是在相平衡情况良好的馈线上继续用单相模型但把结果作为初步方案再通过三相潮流软件校核。我在实际项目里通常选第二种效率高工程上够用。鲁棒性方面如果光伏出力预测误差较大可以采用鲁棒优化或场景法。鲁棒优化的思路是把不确定性集合引入电压约束和网损目标求得一个在最恶劣情况下都可行的解场景法则是抽多个典型场景分别求解再聚合决策。这两种方法都会增加模型规模但工程意义很大。不过这就展开来说又是另一个大话题了。我自己做这个课题最深的感受不在于学会了某个求解器而在于想明白了一个道理建模时做一次凸松弛比调半年参数都管用。把非凸潮流写成 SOCP 之后整个问题性质变了求解变得稳定、高效、可解释YALMIP 配合 Gurobi 在 33 节点问题上的求解速度基本是秒级这让我有大量时间去打磨数据、做敏感性分析和扩展功能。如果你现在还在为配电网优化结果忽好忽坏发愁不妨先把数学模型清晰完整地列出来再一步步转成锥约束很多看似玄学的问题都会自动消失。代码框架和校核流程都在上面了按这个思路走应该能少踩几个我已经踩平的坑。
