做配电网潮流和电压分析的人大概都经历过这种时刻末端电压已经跌破0.95p.u.变压器分接头也调到了极限可远端节点的电压就是提不上去。这时候问题的根源往往不是分接头不够而是无功功率的布局不合理。我自己最开始做无功优化走的是智能算法路线结果每次运行数值都有波动后来在一次交流中接触到基于二阶锥规划SOCP的凸优化建模思路才意识到这才是能稳定复现、有全局最优理论保证的做法。这篇文章就要把“基于二阶锥规划的配电网无功优化算法”完整拆开用IEEE 33节点系统作为测试平台结合Matlab实现覆盖从数学建模、代码编写到结果验证的完整过程。电力系统方向的研究生、做配电网规划的工程师以及想快速复现优化算法的人都可以把这篇文章当成一份可以直接上手的参考。1. 无功优化到底难在哪为什么我从智能算法转向SOCP1.1 配电网无功优化的真实痛点在配电系统里我们关心三个指标电压质量、网络损耗、设备利用率。无功优化的直接目的就是通过调整无功补偿装置、变压器分接头、分布式电源的无功出力让这三个指标尽量达到理想状态。配电网和输电网不一样线路的R/X比值大不能忽略电阻的影响。这意味着有功功率从变电站传向负荷的过程中会明显拉动电压降落所以电压问题在这一层级特别突出。我遇到过的典型场景包括村镇线路末端重负荷时电压低到0.92p.u.或者大量分布式光伏并网后白天电压被抬到1.07p.u.以上。这两个方向的问题本质上都是无功功率分布失调。传统处理手段无非是几类灵敏度补偿法、智能优化算法、直接调用非线性内点法。灵敏度方法很快但属于近似方法处理不了多约束、强耦合的组合场景智能算法如遗传算法、粒子群算法在学术圈应用很广因为“跑出一个结果很容易”但它有明显的工程短板——初值随机导致结果不稳定参数敏感种群数、交叉概率、惯性权重稍微一变优化结果就跟着飘。更麻烦的是同一套参数多跑几次最优解也可能不同这在出报告、做方案比选时非常尴尬。非线性内点法能处理部分问题但非凸模型下对初值敏感容易收敛到局部最优点哪怕结果是可行的也不一定是最省电的方案。所以我倾向的做法是用凸优化里的二阶锥规划SOCP来建模。它保留精确潮流方程不是线性化近似同时有全局最优的理论保证。配上Matlab和YALMIP这一套成熟工具链工程上完全可用也需要不了多少运行时间。1.2 SOCP是怎么把非凸问题“熨平”的先简单说一下非凸问题为什么难办。配电网潮流方程是二次的电流等于功率除以电压功率又是电压和电流的乘积三个物理量互相纠缠导致可行域不是凸集。打个比方非凸问题的可行域像一片有很多凹陷的山谷里面有无数个低洼点传统算法从不同位置起步下山可能分别陷在不同的洼地没办法保证找到全局最低处。SOCP的思路是换一个角度看待问题。对于辐射状配电网潮流可以用一组称为DistFlow或Branch Flow的方程描述。经过变量替换不再直接用电压V和电流I而是用电压平方U和电流平方L绝大多数方程都会变成线性方程只剩一条支路上P²Q²和U、L的耦合关系是非凸的。再把这个耦合等式改写成不等式并转换成二阶锥约束整个优化问题就变成一个凸优化问题。凸问题的好处是可行域像一个光滑的碗没有隐藏的凹坑求解器可以稳定地找到全局最优解。YALMIP和主流求解器都能高效处理二阶锥约束即使IEEE 33节点这种规模的问题从建模到拿到结果通常也就秒级的事。1.3 松弛的精确性为什么敢把等式改成不等式这里有一个绕不开的问题我们把非凸等式L·U P²Q²松弛成了锥不等式L·U ≥ P²Q²松弛之后得到的最优解还能还原成原问题真正对应的潮流解吗理论上已经有成熟结论。只要网络拓扑是辐射状的也就是连通且无环并且目标函数是支路电流平方L的单调增函数那么最优点处锥约束一定是紧的也就是说不等号会恰好取到等号。网损最小化这个典型目标写出来是min Σ(r_ij·L_ij)网损确实随L单调递增。因此IEEE 33这种辐射状网络做网损最小化时SOCP的解可以精确还原物理可行的潮流。不过理论保证不等于实现不需要验证。我自己做项目时一直保留一个后处理步骤计算每条支路的gap L·U(fbus) - (P²Q²)检查最大gap是否接近数值零。如果gap明显非零说明模型里混入了破坏单调性的项比如电压偏差权重过大等这时候求出来的就不是想要的原潮流解。这个校验是很多复现SOCP的人容易忽略的第4章会给出具体的检查脚本。2. 从潮流方程到二阶锥松弛推导与建模要点2.1 IEEE 33节点测试系统的结构与特点IEEE 33节点系统是配电网优化领域的经典标准算例。拓扑上是一个辐射状网络节点1是变电站出线也就是根节点通过32条分支支路向其余节点供电系统电压等级是12.66kV基态总有功负荷约3.7MW无功负荷约2.3Mvar。它最典型的特征就是中后段节点电压偏低最远端节点在基态下普遍低于0.95p.u.非常适合用来验证无功优化算法对电压抬升和网损降低的效果。我拿到数据之后第一件事是把网络数据标准化成三个表节点表包含节点编号、节点类型、有功负荷、无功负荷支路表包含支路编号、首端节点、末端节点、支路电阻、支路电抗无功补偿候选点表写明哪些节点允许安装无功补偿设备容量上下限是什么。这里强烈建议把全部数据统一成标幺值系统。我习惯取基准功率Sb10MVA、基准电压Vb12.66kV于是阻抗基准Zb Vb²/Sb ≈ 16.02Ω。原始数据里的欧姆值除以16.02就是标幺阻抗功率除以10就是标幺功率。单位不统一是后续很多数值问题的根源这一步不能省。很多开源代码默认用有名值也能跑但在求解大规模多时段模型时数值稳定性会明显变差所以一开始就统一单位是最省心的。2.2 DistFlow方程推导从物理量到线性量辐射状网络中任意一条支路(i,j)上的功率传输都满足几个基本关系。记P_ij、Q_ij为从节点i流向节点j的有功和无功U_i为节点i电压幅值的平方L_ij为支路电流幅值的平方有功平衡流进节点j的有功功率扣除支路电阻上的损耗r_ij·L_ij后等于节点j所有下游支路的流出有功之和再加上该节点的净注入有功需求。无功平衡形式完全相同只是把电阻r换成电抗x。电压降落方程U_j U_i - 2(r_ij·P_ij x_ij·Q_ij) (r_ij² x_ij²)·L_ij。支路电流与功率的耦合关系L_ij·U_i P_ij² Q_ij²。前三个方程在(U, L, P, Q)这组新变量下全部是线性的。这一点极其重要它把潮流计算的复杂性集中到了第四个子方程上。第四个方程是二次等式属于非凸约束。我们把它松弛成不等式L_ij·U_i ≥ P_ij² Q_ij²。这个不等式可以等价改写为二阶锥形式也就是YALMIP里可以直接使用的cone约束|| (2P_ij, 2Q_ij, L_ij - U_i) ||₂ ≤ L_ij U_i。验证也很简单左边平方等于4P²4Q²(L-U)²右边是(LU)²要求左边小于等于右边并化简正好得到L·U ≥ P²Q²。这就是SOCP模型中锥约束的具体来源。2.3 目标函数与约束条件的设计对于配电网无功优化最常用的目标函数是系统总有功网损最小化min F Σ r_ij·L_ij。这是L的线性函数天然满足上一节说的单调性条件锥松弛紧度有理论保证。如果还想兼顾电压质量可以在目标函数里加一个电压偏差项比如min F Σ r_ij·L_ij λ·Σ|U_i - U_ref|。绝对值项可以通过引入非负辅助变量线性化问题依然是凸的。但注意λ的选择会影响锥松弛的紧度这点我之前踩过坑后面详细展开。约束条件按物理意义分成五类节点功率平衡约束对每个非平衡节点流入节点的功率扣除支路损耗后加上本地发电和补偿出力等于流出功率加上负荷。平衡节点电压固定U_root V_ref²。电压上下限U_min ≤ U_i ≤ U_max常见取0.95²到1.05²。支路电流上限L_ij ≤ I_max²。无功补偿容量范围Qc_min ≤ Qc_k ≤ Qc_max。每条支路的锥约束||(2P, 2Q, L-U_i)||₂ ≤ LU_i。到这一步优化模型的结构就已经完整了。剩下的事情是把纸上模型翻译成Matlab代码。3. IEEE33节点系统建模与Matlab代码实现3.1 用结构体组织网络数据写可复用的代码我一般不把网络数据硬编码进约束循环里而是先整理成结构体或者矩阵。Matlab下可以这样组织% 基准值定义 Sb 10; % 基准功率 10MVA Vb 12.66; % 基准电压 12.66kV Zb Vb^2 / Sb; % 阻抗基准约16.02欧姆 I_base 1000*Sb/(sqrt(3)*Vb); % 基准电流用于约束支路电流上限 % bus: 节点编号 类型 有功负荷kW 无功负荷kvar bus [ 1 1 0 0; 2 2 100 60; 3 2 90 40; ... ]; % branch: 支路编号 首端节点 末端节点 电阻欧姆 电抗欧姆 branch [ 1 1 2 0.0922 0.0470; 2 2 3 0.4930 0.2511; ... ];然后做标幺化处理r_pu branch(:,4) / Zb; x_pu branch(:,5) / Zb; Pload_pu bus(:,3) / 1000 / Sb; % kW - MW - p.u. Qload_pu bus(:,4) / 1000 / Sb;有一个细节容易被忽略IEEE 33节点的无功负荷数据在不同文献里有不同版本有些版本单位是kvar有些直接给了优化前的总无功还有些已经扣除了并联补偿。不同版本会直接影响基态网损和优化后的电压曲线所以建议第一次运行后先对照节点电压和总网损确认自己录入的数据与预期一致。这一步不核对后面算出来的结果再漂亮也很难对比验证。3.2 构建上下游支路关系在辐射状网络中每条支路有明确的方向。IEEE 33原始数据里支路的首端节点编号通常小于末端节点编号基本按从变电站向外扩展的顺序排列。但为保险起见最好先确认支路方向是朝负荷方向走的。在Matlab里我经常直接通过find找上游和下游支路索引这样功率平衡方程的物理含义最直观。fb branch(:,2); tb branch(:,3); % 找出每个节点的上游支路和下游支路集合 upstream_idx cell(nb,1); downstream_idx cell(nb,1); for j 2:nb upstream_idx{j} find(tb j); % 上游支路末端是j downstream_idx{j} find(fb j); % 下游支路首端是j end写循环加find只是为了可读性IEEE 33节点规模完全够用。如果以后扩展到几百上千节点的网络再用稀疏矩阵关联表优化。3.3 YALMIP变量定义与三类核心约束定义决策变量nb size(bus,1); nl size(branch,1); nc length(candidate_nodes); % 候选无功补偿节点数量 U sdpvar(nb,1); % 节点电压平方 L sdpvar(nl,1); % 支路电流平方 P sdpvar(nl,1); % 支路有功 Q sdpvar(nl,1); % 支路无功 Qc sdpvar(nc,1); % 无功补偿容量 Constraints []; % 电压上下限 Constraints [Constraints, U_min U U_max]; % 支路电流上限 Constraints [Constraints, 0 L Imax_pu^2]; % 无功补偿容量 Constraints [Constraints, Qc_min Qc Qc_max];然后是节点功率平衡约束。这里要特别注意支路损耗的处理。如果P_ij是支路在首端节点处的输送功率那么末端节点实际接收到的功率是P_ij - r_ij·L_ij。功率平衡方程应该写成for j 2:nb up upstream_idx{j}; dn downstream_idx{j}; % 接收上游功率(扣除损耗) 本地发电 - 负荷 流出到下游的功率 Constraints [Constraints, ... sum(P(up) - r_pu(up).*L(up)) Pgen(j) - Pload(j) ... - sum(P(dn)) 0]; end对应无功分支加上Q以及Qc补偿出力格式完全相同。电压降落方程也是逐支路写for k 1:nl f fb(k); t tb(k); Constraints [Constraints, ... U(t) U(f) - 2*(r_pu(k)*P(k) x_pu(k)*Q(k)) ... (r_pu(k)^2 x_pu(k)^2)*L(k)]; end锥约束是SOCP最核心的部分for k 1:nl f fb(k); Constraints [Constraints, ... cone([2*P(k); 2*Q(k); L(k)-U(f)], L(k)U(f))]; endYALMIP的cone函数语法是cone(v, y)表示y ≥ ||v||₂。因此必须确保第一个参数是差分向量第二个参数是求和项UL顺序不要写反。很多复现出错的人就是在这里把L-U和U-L写反导致锥的意义完全变了结果对不上潮流方程。最后补上根节点电压V0 1.0; % 根节点电压可按变压器出口分接头设定为1.0或1.05 Constraints [Constraints, U(1) V0^2];3.4 求解器选择与常见报错排查求解部分最省心的写法是用YALMIP配Gurobi或CPLEXops sdpsettings(solver,gurobi,verbose,1); sol optimize(Constraints, obj, ops); if sol.problem ~ 0 warning([求解失败: sol.info]); end U_opt value(U); L_opt value(L); P_opt value(P); Q_opt value(Q); Qc_opt value(Qc);求解器层面我的优先顺序是Gurobi、CPLEX、MOSEK、Sedumi、SDPT3。前三个商业求解器的数值稳定性更好尤其Gurobi处理二阶锥非常快IEEE 33节点这种规模通常在很短时间内就能收敛。Sedumi是社区经典求解器对付小规模SOCP没问题但模型复杂度上来后会慢一些。常见问题有三种无解(infeasible)优先检查电压上下限U_min和U_max是否矛盾以及无功补偿容量上限是否设得太小。求解器报numerical issue十有八九是单位混用。把阻抗、功率全部统一成标幺值再跑一遍问题就会消失。结果异常比如补偿点几乎不出力检查功率平衡约束的负荷符号是否搞反了先用基态潮流验算一遍。这里特别提醒在YALMIP里千万不要直接写P(k)^2 Q(k)^2 U(f)*L(k)。这个约束虽然数学上等价于锥约束但是以双线性约束形式给出的YALMIP不一定能识别其凸性可能无法调用合适的SOCP求解器运行时间和求解稳定性都会大受影响。用cone函数显式声明锥约束是我踩过坑之后总结出的固定习惯。4. 优化结果的三重验证电压、网损与锥紧度4.1 第一重节点电压曲线对比模型求解完成后第一个要看的指标是节点电压分布。我一般会把优化前基态潮流和优化后的电压画在同一条曲线图里figure; plot(1:nb, sqrt(U_base), o-, LineWidth, 1.5); hold on; plot(1:nb, sqrt(U_opt), s-, LineWidth, 1.5); yline(sqrt(U_min), --, 电压下界); yline(sqrt(U_max), --, 电压上界); xlabel(节点编号); ylabel(电压幅值 (p.u.)); legend(优化前基态,优化后,下界,上界); grid on;IEEE 33节点基态下最远端节点电压通常处于0.91到0.93p.u.之间具体数值取决于录入的无功负荷版本。加入无功补偿后末端电压会被抬升到0.95p.u.以上靠近中后段的关键节点也会回到安全区间。如果优化后电压曲线仍有节点低于下限不要急着改代码先检查无功补偿总容量是否不足。IEEE 33系统按典型数据配置需要0.3到0.4Mvar左右的总补偿容量就能把电压拉回正轨如果给的容量余量太小电压肯定提不上去。4.2 第二重网损降低效果网损是评价无功优化经济性的核心指标。计算式很简单loss_base_pu r_pu * value(L_base); loss_opt_pu r_pu * L_opt; loss_base_MW loss_base_pu * Sb; loss_opt_MW loss_opt_pu * Sb; reduction (loss_base_MW - loss_opt_MW) / loss_base_MW * 100;对于IEEE 33节点系统基态网损典型量级在0.2MW左右。做了无功优化后网损一般可以下降20%到40%具体取决于允许安装的无功容量和是否加入电压偏差权重。如果目标函数只写网损最小、补偿容量给得又宽裕网损会压得很低如果加入了电压偏差项电压分布会更平均一些但网损降幅会相应变小两者本质上是权衡关系。额外可以画一张“优化后各支路网损占比”的柱状图看看哪些线路被重点改善。通常会集中在距离变电站较远、负荷集中的馈线段这也从侧面验证了无功补偿落点选择的合理性。比如某个长线路末端电压低那么补偿点往往会落在它的下游节点附近而不是随便选一个靠近变电站的节点。4.3 第三重锥紧度校验决定你的解是否可信这一节是几乎所有复现SOCP的人都容易忽略的。很多教程和代码能算出数字却从不去验证这个数字对应的到底是不是原始潮流方程的解。锥紧度校验代码如下gap L_opt .* U_opt(fb) - (P_opt.^2 Q_opt.^2); max_gap max(abs(gap)); tol 1e-6; % 标幺制下容差 if max_gap tol disp(锥紧度校验通过: 优化解满足原始潮流方程); else warning(锥约束未紧: 最大gap%.2e, max_gap); [~, idx] max(gap); fprintf(最大gap对应支路: %d\n, idx); end我在项目里的经验是只要目标函数是纯网损最小并且根节点电压固定gap量级通常在1e-8以下说明松弛是精确的。但一旦引入电压偏差惩罚项并且权重λ取得过大就会出现gap上升到1e-3甚至更大的情况。原因是目标函数对U的单调性被削弱后锥约束不再被推挤到边界。这时求出来的解虽然满足模型的一切约束但它已经不能直接解读为物理潮流的可行点。解决思路是调整λ或者改用分段线性惩罚的形式并保持网损项在目标中的主导地位。校验这个环节能帮你避免很多“看起来对、实际不对”的结果也是有经验的工程师最容易追问的地方。5. 工程化扩展离散电容、分布式电源与多时段联调5.1 从连续无功到分组电容引入0-1整数变量前面的模型假设无功补偿容量可以连续调节这对应SVC、SVG或者逆变器无功这类设备。但实际配网里最常见的补偿设备是分组电容器容量按固定档位投入。此时决策变量从连续变量变成整数变量模型升级为混合整数二阶锥规划也就是MISOCP。YALMIP里的改动很简单% 假设每个补偿点最多有Ng组每组容量Qstep q_unit binvar(nc, Ng); % 0-1投切变量 Qc Qstep * sum(q_unit, 2); % 总补偿容量 Constraints [Constraints, Qc_min Qc Qc_max];求解器要换成支持混合整数的版本Gurobi、CPLEX都没问题。如果还在用Sedumi那它无法直接处理二元变量需要外部分支定界逻辑或者换MOSEK。引入整数变量后运行时间会从秒级上升到十几秒到几十秒对于33节点规模依然完全可接受。有一点需要额外注意整数变量会让SOCP的锥紧度不如连续模型那么完美。但只要目标仍然是网损最小大多数支路的gap仍能保持在可接受范围内校验环节照旧要做。不要因为用了MISOCP就跳过紧度检查。5.2 把分布式电源接入场景纳入模型现在配网里光伏渗透率越来越高白天光伏出力大、负荷小电压容易被抬升晚上负荷增大了光伏又不出力电压又偏低。很多地区的运行要求就是利用逆变器剩余容量提供无功支撑。模型上的改动是把每个光伏节点的有功出力Ppv作为给定参数由预测曲线事先确定无功出力Qpv作为优化变量同时加一条容量约束Q_pv² ≤ S_pv² - P_pv²。这条约束是二次的但它是标准的旋转二阶锥约束YALMIP里可以用cone改写后加入模型整个问题依然是SOCP。这样我们就能回答一个很实际的问题在给定的光伏出力曲线下各台逆变器应该发多少无功无功应该主要从哪几个节点发出才能在电压不越限的前提下让网损最小。5.3 多时段模型与设备动作次数约束单断面优化回答不了“分接头一天该动作几次”“电容器投切会不会太频繁”这类问题。做法是把24个时间断面的变量全部并列进一个模型每个断面沿用上述潮流约束断面之间通过共享设备状态变量建立联系。目标函数是24个断面的总网损之和另加动作次数约束例如Σ|q_unit(t) - q_unit(t-1)| ≤ MaxSwitchTimes。这类绝对值约束可以用辅助变量线性化。模型规模会成倍扩大但在现代求解器面前仍然可解只是建模时要注意矩阵维度对应关系不要写错。我在实际项目里通常先做单断面把候选补偿点筛选出来再进行24断面精细化校核两个步骤解耦处理效率比一次性直接建大规模模型高很多。做多时段模型最常碰到的坑是数值规模翻倍后求解器内存和迭代步数快速上升。这时候建议给每个线性约束分配独立的变量块适当放宽求解器的MIPGap参数。工程上3%以内的gap已经可以接受不需要一味追求绝对最优但必须保持优化结果与单目潮流校核的一致性。回到开头那个场景如果你也被配电网末端电压和网损问题困扰不妨把SOCP模型作为工具箱里的标准流程。我常用的组合是IEEE 33节点系统做验证、YALMIP加Gurobi做求解、锥紧度校验做兜底。这套方法跑通之后我最大感受是以前调智能算法参数可能要花大半天现在改约束、换场景基本是分钟级的事而且结果稳定、可复现、讲得出道理。希望这篇拆解能帮你少走我当年走过的弯路。