1. 无功优化到底在优化什么做配电网研究的同学十有八九都绕不开无功优化这道坎。用IEEE33节点系统做基于二阶锥规划SOCP的配电网无功优化是论文和课程设计里出现频率最高的组合之一题目一报出来外行看到的是配电网优化算法内行看到的是潮流方程的非凸性怎么处理、补偿设备怎么选、目标函数怎么定这一连串实际问题。我自己从最早用粒子群硬搜到后面换成非线性规划再到现在用YALMIP把整个模型写成SOCP中间的坑踩了不少。这篇文章就按我实际跑通的思路把模型怎么建、参数怎么选、代码怎么写、结果怎么验从头到尾过一遍保证你能直接在Matlab里复现出可用的结果。1.1 配电网的无功电压问题为什么比输电网更棘手配电网和输电网最大的区别在于R/X比。输电网的线路电抗远大于电阻有功和无功基本可以解耦分析无功补偿通常只在变电站侧集中做配电网里电缆和架空线的电阻占比很大R/X比经常是1到3的量级意味着有功潮流会明显影响电压无功潮流不但影响电压还直接在线路上变成热损耗。更麻烦的是配电网末端负荷波动大、量测设备少很多台区连基本的自动化都没有电压越限的重灾区通常就发生在那些又长又细的末端支路上。无功优化要解决的本质上就是用最小的控制代价让系统电压回到合格区间同时把网损压下去。反应在物理上就是让无功尽量就地平衡——负荷在哪儿需要无功补偿设备就放在哪儿附近别让大量无功从变电站远端一路送过来。线路里跑的无功少了电流就小了线损自然降下来末端电压也会明显改善。这个道理说起来简单落到数学上却是一个高维、非线性的优化问题难点全在怎么又快又稳地找到最优解。1.2 目标、手段和约束先把框架搭清楚做优化第一步不是写代码而是先明确三件事目标函数、控制手段、约束条件。目标函数方面主流选项有四个网损最小、电压偏差最小、无功补偿经济性最优、多目标加权。IEEE33这类标准测试系统里最常用的就是网损最小因为它是无功优化的初心物理意义清楚、便于和其他文献对比。如果你想顺便改善电压可以把电压偏差作为一个惩罚项加进去但要注意权重设置这个我后面专门讲。控制手段方面模型里最方便处理的是连续无功补偿容量可以把它理解成逆变器型分布式电源的无功输出或者SVC/SVG这类连续无级补偿装置的输出值。如果要模拟传统的分组投切电容器就把补偿变量改成整数变量问题会从SOCP变成混合整数二阶锥规划MISOCPYALMIP加Gurobi或CPLEX能直接解只是求解时间会长一些。OLTC有载调压和无功优化联动在配网里也很常见但IEEE33基础版本里一般不建模变压器档位先把补偿设备做好再扩展也不迟。约束条件分四类。第一类是潮流平衡约束这是所有电力系统优化问题的宪法不能违背第二类是节点电压上下限国内工程上通常按0.931.07pu学术论文里更常用0.951.05pu第三类是支路电流上限给线路容量一个保护第四类是补偿设备自身的容量限制。约束不是越多越好尤其是IEEE33这种测试系统加太多约束反而会掩盖算法本身的对比价值新手容易犯的错就是什么约束都往里塞结果既难调参又看不出性能差异。1.3 IEEE33节点为什么它是配网优化的标准考场IEEE33节点系统来自Baran和Wu在1989年发表的论文原始用途就是配电网无功优化和网络重构研究。系统电压等级12.66kV33个节点、32条支路总负荷3715kW加2300kvar呈辐射状结构从变电站出发分出三条主馈线再加若干分支拓扑一看就很有代表性。这个系统最讨喜的地方在于规模不大做验证速度快但它有真实的电压越限问题——不装任何补偿设备时末端最低电压只有0.913pu左右全网有功损耗约202.7kW降损空间非常明显。正因为有这些病灶IEEE33成了检验配网优化算法的标准考场。你优化前后一对比电压抬升了多少、损耗降了多少数字摆在那里效果一目了然。再加上网上数据容易找新算法拿它先跑一遍几乎成了行规。做这个项目时我建议你把支路参数和负荷数据单独存成一个m文件或CSV别每次拷贝粘贴——后面调参、做对比实验会频繁用数据文件独立出来能省非常多事。2. 为什么偏偏是二阶锥规划2.1 数学上先看清楚潮流的killer非凸性配电网的潮流方程里藏着一个非常不友好的东西支路电流的平方等于视在功率的平方除以电压幅值的平方也就是 I² (P² Q²)/V²。这个关系里P、Q、V都是变量P²Q²除以V²是一个非凸约束直接把整个优化问题变成了非凸问题。非凸意味着什么意味着你没法保证找到的是全局最优解梯度类算法经常陷在局部最优里出不来初值选得不好跑几次得到的结果可能相差很大。这就是为什么以前大家做配网无功优化要么用内点法硬解NLP要么上粒子群、遗传算法这些启发式算法。内点法看着优雅但非凸NLP对初值极度敏感同一个题设换一个初值可能收敛到完全不同的局部解启发式算法虽然能全局搜索但计算量大、稳定性差而且没有任何最优性保证论文里用一用没问题工程上可不敢这么玩。2.2 DistFlow加锥松弛把硬骨头变软既然潮流方程非凸那就想办法让它变凸。配电网是辐射状结构这给了我们一个独门武器——DistFlow模型。DistFlow把每个节点的功率平衡和电压关系写成递归形式对每条支路(i,j)有v_j v_i - 2(r_ij·P_ij x_ij·Q_ij) (r_ij² x_ij²)·l_ij其中v_i是节点i电压幅值的平方l_ij是支路电流幅值的平方。这个方程本身是线性的真正非凸的是l_ij和v_i、P_ij、Q_ij之间的关系本来 l_ij·v_i P_ij² Q_ij² 是个等式约束现在把等号放松成不等号让 l_ij·v_i ≥ P_ij² Q_ij²再做一次等价变形就得到标准的二阶锥约束|| [2P_ij; 2Q_ij; v_i - l_ij] ||₂ ≤ v_i l_ij你可能想问把等式放松成不等式解还靠谱吗这就涉及SOCP最核心的工程直觉目标函数是最小化网损 sum(r_ij·l_ij)而电阻r_ij恒为正优化器会拼命把l_ij往小里压所以在最优解处这个锥松弛通常会被拉紧到恰好取等号也就是说放松前后的解是一样的。这就叫精确松弛。数学上Farivar、Low等人在2011到2013年的一系列工作证明了对于辐射状配电网在总负荷有上界、目标函数关于潮流单调递增等温和条件下SOCP松弛是精确的。实际做IEEE33优化时松弛间隙小到可以忽略。这一点很有意思——本来是个让人头疼的非凸问题靠着网络结构直接降维成凸问题求全局最优从一个奢望变成了多项式时间内轻松搞定的事。2.3 与粒子群、传统内点法的横向对比方法能否保证全局最优计算效率处理离散设备实现难度非线性规划内点法否易陷局部最优中等需特殊处理中粒子群/遗传算法否无最优性证明慢依赖种群参数天然支持较低二阶锥规划SOCP是多项式时间很快连续变量天然支持中混合整数二阶锥MISOCP是较慢支持分组电容器中高在IEEE33这种几十个变量的规模下SOCP的求解时间通常在秒级以内。我第一次用粒子群算同样的问题种群100、迭代200次跑了很久才出一个还不一定能复现的结果换成YALMIP加Cplex几秒钟就出全局最优。后来我也见过不少同学还在用遗传算法算IEEE33的无功优化不是说不能用但在需要反复对比、批量仿真的场景下SOCP的优势是压倒性的。3. 建模仿真前要先搞定的参数细节3.1 基准值怎么选一个小数点引发的求解器血案做IEEE33仿真最容易翻车的不是算法而是基准值。这个系统原始数据是12.66kV支路阻抗单位是欧姆负荷单位是kW/kvar不标幺直接建模方程里会出现0.0922和420这种数量级差几千倍的数求解器要么罢工要么给出精度极差的解。通常做法是取V_base12.66kV、S_base10MVA这样阻抗基准值Z_baseV_base²/S_base160.2756Ω支路阻抗标幺值落在0.001到0.02之间负荷标幺值落在0.004到0.06之间各变量的数量级基本一致数值条件非常好。有个细节值得提醒有些同学图省事取S_base100MVA算出来的支路阻抗就变成1e-4量级电压方程里2(rPxQ)这种项只有1e-5再叠加上二阶锥约束里的微小数字求解器的数值容差一挑动结果就容易出现假收敛。如果你用的是老版本求解器这种问题更明显。我直接给结论IEEE33就用10MVA基准值少走弯路。另外要注意单位换算。负荷表给的是kW和kvar标幺化要除以S_base*1000因为10MVA等于10000kVA不是直接除以10。好多同学第一步就在这错了后面电压约束再怎么调都无济于事检查半天最后发现负荷标幺值整体大了1000倍那种心情我太懂了。3.2 候选补偿节点怎么挑补偿点位选择直接决定优化效果。IEEE33里优先考虑的节点有三类一是末端节点因为末端电压最低、边际效益最大典型的就是18、22、33这三条馈线的末端二是重负荷节点比如24、25号节点负荷各420kW是整个系统最重的负荷群旁边放补偿效果立竿见影三是带重无功负荷的特殊节点比如30号节点带600kvar的无功负荷电压掉得特别明显也是热门候选。我比较推荐的最基础配置是[18; 22; 25; 30; 33]这五个候选点补偿容量上限分别设为300、300、600、600、300kvar。这个配置之所以这么定原因很实在三个末端点各给300kvar把最典型的低电压问题压住25和30号这两个重负荷点给到600kvar因为它们的负荷体量在那摆着给小了没效果给大了又会反向送无功。如果你想做得更深入可以把候选点扩到16、18、22、25、30、33甚至把每个候选点的容量上限也作为一个变量去做灵敏度分析看哪个点对降损贡献最大——这其实就是配电网无功规划里的选址定容问题了。3.3 目标函数和约束的取值范围建议目标函数我建议先只写网损最小obj sum(r .* Iij2)。这个目标对SOCP最友好因为电阻恒正能保证锥松弛被拉紧。如果你想兼顾电压可以用 obj sum(r.Iij2) lambdasum((Vi2-1).^2)。lambda的取值要小心IEEE33的网损标幺值在0.01到0.02量级而电压偏差平方和可能到0.01到0.05直接相加会让电压项喧宾夺主。我试下来lambda在0.1到1之间比较合适而且建议先把lambda0的纯网损结果跑出来再逐渐增大lambda观察网损和电压两个指标此消彼长的关系最后再定权重。电压约束方面IEEE33不加补偿时最低电压约0.913所以如果你把下限设成0.95那么优化前的工况在这个约束下就是不可行的。这不是bug反而说明问题有优化空间。做对比实验时基准工况无补偿的电压下限要放宽到0.88到0.90或者直接用解潮流的工具算出真实基态优化工况再用0.95的合格下限。这个细节很多新手没注意把基准工况和优化工况用同一套电压约束结果曲线一出来对比就失去了说服力。补偿容量的上限按节点负荷的无功分量估算最合理。经验值是取该节点Q_load的20%到50%或者参照实际工程里单台电容器的额定容量按30到100kvar的档位去凑。25号节点负荷重给到600kvar很合理30号节点本身带600kvar的无功负荷给600kvar也符合工程直觉末端节点给300kvar就够了。给太多会让电压抬过头、容性无功反向流动反而白白增加损耗。4. Matlab YALMIP 核心代码实现4.1 环境配置与求解器选型跑SOCPMatlab本身不算关键关键是配套的建模工具和求解器。我强烈建议用YALMIP做建模层它把变量声明、约束拼接、求解器调用统一起来比手写Cplex API或直接用CVX都省心。求解器方面优先Gurobi或Cplex它们都原生支持二阶锥约束且能解MISOCP对IEEE33这种规模几秒钟就完成求解。如果学校没有商业许可可以用免费开源的SeDuMi或SDPT3解连续SOCP完全没有问题只是速度稍慢、不能直接处理整数变量。这里特别说一个很多人问的点装了Gurobi还能不能装Cplex答案是完全可以8这两个工具在MATLAB路径里可以和平共处YALMIP会全部识别运行时用sdpsettings指定具体用哪个求解器就行。我的习惯是先用SeDuMi跑通逻辑再换Gurobi跑大规模和整数版本两者互不干扰。装好后可以用yalmiptest这个命令检查YALMIP到底识别了哪些求解器如果有显示No suitable solver多半是路径没加对。4.2 IEEE33完整数据和预处理先准备数据文件。我建议把基础数据独立存成一个函数文件case33_ieee.m返回基准电压、基准容量、支路表和负荷表四个量之后所有脚本都从它读取逻辑清晰、可复用。下面是完整数据摘自经典IEEE33节点算例可以直接保存使用。function [Vb, Sb, br, ld] case33_ieee % IEEE33节点配电网基础数据 % 输出: Vb-基准电压(kV), Sb-基准容量(MVA), br-支路表, ld-负荷表(kW/kvar) Vb 12.66; Sb 10; br [ 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 ]; ld [ 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 ]; end数据文件准备好之后在主脚本里调用然后做标幺化。这块代码并不长但每行都有讲究[Vb, Sb, br, ld] case33_ieee; nbus 33; nbranch size(br, 1); Zb Vb^2 / Sb; % 阻抗基准值单位欧姆 r br(:,3) / Zb; % 支路电阻标幺 x br(:,4) / Zb; % 支路电抗标幺 pL ld(:,2) / (Sb*1000); % 有功负荷标幺kW / 10000 qL ld(:,3) / (Sb*1000); % 无功负荷标幺kvar / 10000这里最容易错的就是最后两行。Sb单位是MVASb*1000才是kVA负荷的kW和kvar必须除以这个数才能得到正确的标幺值。我把这行写在最醒目的位置就是希望正在调试的你少走一次弯路。4.3 变量定义与约束构建代码逐段拆解变量定义直接用YALMIP的sdpvar。支路有功Pij、支路无功Qij、支路电流平方Iij2都是32维向量节点电压平方Vi2是33维向量无功补偿qC是33维向量其中只有候选节点才能真正出力。目标函数用网损加电压偏差惩罚项的形式后面你可以自己调lambda。% 决策变量 Pij sdpvar(nbranch, 1); % 支路有功 Qij sdpvar(nbranch, 1); % 支路无功 Iij2 sdpvar(nbranch, 1); % 支路电流平方 Vi2 sdpvar(nbus, 1); % 节点电压平方 qC sdpvar(nbus, 1); % 无功补偿容性为正 % 约束容器 cons []; % 根节点电压等于1 cons [cons, Vi2(1) 1]; % 电压幅值上下限 cons [cons, 0.95^2 Vi2 1.05^2];电压上下限这个地方优化工况用0.95到1.05这是学术论文最常用的区间。如果你做基准工况对比记得把下限放宽否则无补偿的基态在模型下不可行。接下来是最核心的部分支路电压方程和二阶锥约束。对每条支路先写电压降落关系再写SOC约束。这里的循环清晰直观每条支路两行约束不容易出错。% 支路方程 二阶锥约束 for k 1:nbranch i br(k,1); % 首端节点 j br(k,2); % 末端节点 % 电压平方关系 cons [cons, Vi2(j) Vi2(i) - 2*(r(k)*Pij(k) x(k)*Qij(k)) ... (r(k)^2 x(k)^2)*Iij2(k)]; % 二阶锥松弛约束 cons [cons, norm([2*Pij(k); 2*Qij(k); Vi2(i) - Iij2(k)], 2) ... Vi2(i) Iij2(k)]; enddistflow模型的这一行是支路潮流的核心它描述的物理过程是末端电压平方 首端电压平方 - 线路压降 线路损耗项。二阶锥约束这一行就是我前面讲的把非凸等式放松成凸不等式的关键一步YALMIP会把norm(x,2)y这种形式自动识别成二阶锥交给求解器时不必手动转换格式。节点功率平衡约束是另一个关键。对每个非根节点流入该节点的功率要等于流出功率加上净负荷。这里我用了find来定位每个节点的父支路和子支路循环清晰也方便你以后改成其他拓扑结构。% 节点功率平衡1号节点是平衡节点不写 for j 2:nbus par find(br(:,2) j); % 以j为末端的分支功率流入 chd find(br(:,1) j); % 以j为首端的分支功率流出 pbal 0; qbal 0; for kk par pbal pbal (Pij(kk) - r(kk)*Iij2(kk)); qbal qbal (Qij(kk) - x(kk)*Iij2(kk)); end for kk chd pbal pbal - Pij(kk); qbal qbal - Qij(kk); end cons [cons, pbal pL(j)]; % 有功平衡 cons [cons, qbal qC(j) qL(j)]; % 无功平衡 end无功平衡式的符号要特别注意qbal是支路送来的无功qC是本地补偿注入的容性无功为正两者之和等于负荷的无功消耗qL。如果把qC的符号搞反求解器要么给出明显不合理的补偿结果要么直接报不可行。补偿变量加容量约束也很简单。用零向量初始化上限只给候选节点赋值。这种写法比直接写五个单独约束清晰得多后面想调整候选节点也方便。% 无功补偿容量限制 capNodes [18; 22; 25; 30; 33]; qCapMax zeros(nbus, 1); qCapMax(capNodes) [0.03; 0.03; 0.06; 0.06; 0.03]; % 对应kvar数 cons [cons, 0 qC qCapMax];比如0.03的标幺值在10MVA基准下就是300kvar0.06就是600kvar。这个容量上限给的依据就是前面说的节点负荷体量和实际工程中电容器的常见容量档位。4.4 求解与结果导出目标函数和求解部分非常短。先把目标函数定义成网损加惩罚项然后调optimize函数。sdpsettings里指定求解器我一般用cplex或gurobi如果要用免费求解器就改solver为sedumi。% 目标函数网损最小 轻量电压偏差惩罚 lambda 0.5; obj sum(r .* Iij2) lambda * sum((Vi2 - 1).^2); % 求解 options sdpsettings(solver, cplex, verbose, 2); sol optimize(cons, obj, options); % 结果检查 if sol.problem 0 fprintf(求解成功\n); else yalmiperror(sol.problem); end求解完成之后把结果从标幺值还原成工程单位。网损是求和后乘以Sb再乘以1000得到kW补偿量乘以Sb再乘以1000得到kvar电压幅值直接对Vi2开方。% 结果换算 V_opt sqrt(value(Vi2)); loss_opt sum(r .* value(Iij2)) * Sb * 1000; % kW Qcomp_opt value(qC) * Sb * 1000; % kvar fprintf(优化后网损: %.2f kW\n, loss_opt); fprintf(最低电压: %.4f pu 节点%d\n, min(V_opt), find(V_opt min(V_opt))); for kk 1:length(capNodes) n capNodes(kk); fprintf(节点%2d 无功补偿: %.1f kvar\n, n, Qcomp_opt(n)); end到这里一个完整可跑的SOCP无功优化程序就成型了。整个脚本从数据准备到结果输出大约120行调试时间通常控制在几分钟内。如果你跑出来的结果跟常见文献对不上先回去查基准值和负荷单位这两个地方错一个后面全盘皆输。5. 典型结果分析与自检方法5.1 优化前后的电压曲线对比要拿电压曲线说明问题至少要有两条线一条是优化前的基态电压一条是优化后的电压。基态电压怎么来最简单的方法是不加补偿、把电压下限放宽用同一套SOCP模型再跑一次。% 基态求解无补偿 电压下限放宽 cons_base []; % ... 复用上面的约束但去掉补偿约束或令 qC0 % 电压下限改成 0.88^2 而不是 0.95^2 % 目标函数只用网损最小 % 求解后得到 V_base跑完之后画图比较figure; plot(1:nbus, V_base, -o, LineWidth, 1.5); hold on; plot(1:nbus, V_opt, -s, LineWidth, 1.5); grid on; xlabel(节点编号); ylabel(电压幅值 (pu)); legend(优化前, 优化后, Location, southwest); ylim([0.90, 1.06]);一个典型的可复现结果是优化前最低电压在18号节点附近约0.913pu优化后全网电压全部抬升到0.95以上最低电压点的位置和幅值都明显改善电压曲线整体变得更平缓。曲线图一出来优化效果非常直观。注意别把两条线画成同一种样式加粗、加标记评审和答辩时能省不少解释时间。5.2 网损下降率怎么算才严谨网损下降率是评审最认的硬指标。计算公式很简单降损率 (优化前网损 - 优化后网损) / 优化前网损 × 100%。难点在于优化前网损必须从同一个系统的同一种运行方式下算出来不能拿别的论文的数字顶替。我自己这次跑的典型结果是优化前约202.7kW优化后在175kW左右降损率大约13.6%。不同候选节点配置、不同容量上限会得到不同数字这个浮动是正常的你要关注的是自己的对比口径是否一致。关于降损率我想多说一句别只盯一个数字。你还要顺手看两个佐证指标全网总无功补偿量是多少电压偏差从多少降到多少。如果降损率很高但补偿容量大得离谱或者电压已经顶到上限那这个结果即使数学上最优工程上也不一定可行。IEEE33只是测试系统但做测试系统的严谨习惯会直接影响你以后做真实配网项目。5.3 松弛间隙检查你的解是不是真的稳SOCP最让人担心的问题就是松弛是否精确。虽然理论上有保证但数值实现和特定目标函数下还是有必要自己检查一遍。方法是求解后把所有最优值代回锥约束看不等号两侧的差距。% 松弛间隙检查 gap_soc zeros(nbranch, 1); for k 1:nbranch i br(k, 1); lhs norm([2*value(Pij(k)); 2*value(Qij(k)); value(Vi2(i)) - value(Iij2(k))], 2); rhs value(Vi2(i)) value(Iij2(k)); gap_soc(k) rhs - lhs; end fprintf(最大松弛间隙: %.3e\n, max(gap_soc));如果这个间隙在1e-6量级甚至更小说明锥约束在最优解处基本取等号你的解是可信的。如果你改了目标函数比如把电压偏差的权重调得很大锥约束可能不再拉紧这时求解器给的最优解就不一定对应真实物理潮流了。这个检查是我的必做步骤也是我判断一个SOCP无功优化代码是否靠谱的第一道门槛。如果条件允许还可以用电力系统潮流软件做交叉验证。Matpower新版本里带了case33bw这个IEEE33数据文件你可以跑一次runpf得到基态潮流再和SOCP模型在零补偿状态下的结果对比电压和网损能对上那基本就稳了。6. 常见问题与排查技巧实录6.1 求解器报错与兼容性问题YALMIP报No suitable solver for a class of problem是最常见的错误之一。这句话的意思是你的模型里包含了当前求解器不支持的约束类型最常见的情况是你用了不支持二阶锥的求解器或者装了求解器但没加进MATLAB路径。先运行yalmiptest看识别列表如果列表里有cplex或gurobi大概率是模型里混入了YALMIP无法自动转化的约束逐个检查约束类型即可。关于多个求解器共存的问题我在4.1已经说过Gurobi和Cplex可以同时安装YALMIP会自己识别。实际中还会遇到一种尴尬就是某个求解器版本太老对MISOCP支持不完整。我遇到过Gurobi老版本报二次约束不是凸的这种误导信息其实问题出在版本太旧。解决办法是升级求解器或者在sdpsettings里强制指定另一个求解器做交叉验证。6.2 不收敛或不可行的排查清单碰到不可行问题先别急着怀疑代码。按这个顺序排查第一步把电压下限从0.95放宽到0.88看问题是否变可行——如果放宽后就可行说明补偿容量不足或候选节点不够需要增大容量或增加候选点第二步去掉二阶锥约束以外的各类附加约束比如电压上限、支路电流上限逐条加回来找到是哪个约束导致不可行第三步检查目标函数和变量是否存在维度不匹配YALMIP的sdpvar乘法很容易在无意中制造非线性项第四步用sol.info和yalmiperror输出详细诊断信息很多不可行问题在诊断信息里会直接告诉你哪个约束冲突。求解器给了可行解但结果不合理比如补偿全为负、网损比基态还高优先检查符号约定。无功补偿的符号、功率平衡方程的正负号、支路功率流向的定义三个环节任何一处反了都会出这种诡异结果。IEEE33的对称性不强节点编号顺时针看是三条馈线但支路顺序不是严格按节点编号排列的用find函数定位分支时要特别留意。6.3 新手最容易踩的5个坑第一负荷单位没转对。kW和kvar要除以S_base*1000而不是除以S_base这个错误我见得太多了它在数值上表现为所有负荷标幺值大了1000倍电压约束形同虚设。第二基准工况和优化工况的约束不一致。跑基态时用了0.95的电压下限结果基态不可行就强行加补偿去救然后两个工况混在一起对比曲线完全没有说服力。正确的是基态放宽电压下限、不加补偿优化态加补偿、用合格下限两个工况的口径分开。第三只盯着网损下降率不检查松弛间隙。我在5.3讲过目标函数换成电压偏差主导后锥约束可能不紧。不紧的SOCP解看起来是凸问题最优解实际上对应的物理潮流是虚的。每跑一版都要带回代检查。第四补偿容量上限拍脑袋。容量上限最好依据节点无功负荷量级来设一个只有90kvar负荷的末端节点你给1200kvar的补偿求解器会算出补偿容量顶到上限但电压已经超限的荒唐结果。先做一次无功负荷分布统计再定容量。第五忽略YALMIP和求解器的版本匹配。老版本YALMIP对某些SOC约束的语法支持不完整新版本求解器又可能改变了默认容差。我建议固定一套经过验证的组合比如YALMIP R2021版配Gurobi 9.x然后所有对比实验都在同一套环境里跑免得被版本差异干扰结论。最后再分享一个我自己的习惯。做这个项目时我会把每一版结果都存成带时间戳的mat或Excel文件跑完一组实验先把电压曲线和网损数字记录到一张对比表里而不是只看命令行输出的最后三行。这样做的好处是论文返修或评审质疑结论时你能快速回溯每一个数字是在什么参数下跑出来的。你是要拿这个模型替自己说话的数据和结论的源头越清楚你说的话才越硬气。SOCP这个方向看起来是公式堆出来的但真正能让你在答辩和工程评审里站稳的是这些围绕模型的验证习惯和细节控制。 配电网无功优化实战：IEEE33节点二阶锥规划Matlab实现
1. 无功优化到底在优化什么做配电网研究的同学十有八九都绕不开无功优化这道坎。用IEEE33节点系统做基于二阶锥规划SOCP的配电网无功优化是论文和课程设计里出现频率最高的组合之一题目一报出来外行看到的是配电网优化算法内行看到的是潮流方程的非凸性怎么处理、补偿设备怎么选、目标函数怎么定这一连串实际问题。我自己从最早用粒子群硬搜到后面换成非线性规划再到现在用YALMIP把整个模型写成SOCP中间的坑踩了不少。这篇文章就按我实际跑通的思路把模型怎么建、参数怎么选、代码怎么写、结果怎么验从头到尾过一遍保证你能直接在Matlab里复现出可用的结果。1.1 配电网的无功电压问题为什么比输电网更棘手配电网和输电网最大的区别在于R/X比。输电网的线路电抗远大于电阻有功和无功基本可以解耦分析无功补偿通常只在变电站侧集中做配电网里电缆和架空线的电阻占比很大R/X比经常是1到3的量级意味着有功潮流会明显影响电压无功潮流不但影响电压还直接在线路上变成热损耗。更麻烦的是配电网末端负荷波动大、量测设备少很多台区连基本的自动化都没有电压越限的重灾区通常就发生在那些又长又细的末端支路上。无功优化要解决的本质上就是用最小的控制代价让系统电压回到合格区间同时把网损压下去。反应在物理上就是让无功尽量就地平衡——负荷在哪儿需要无功补偿设备就放在哪儿附近别让大量无功从变电站远端一路送过来。线路里跑的无功少了电流就小了线损自然降下来末端电压也会明显改善。这个道理说起来简单落到数学上却是一个高维、非线性的优化问题难点全在怎么又快又稳地找到最优解。1.2 目标、手段和约束先把框架搭清楚做优化第一步不是写代码而是先明确三件事目标函数、控制手段、约束条件。目标函数方面主流选项有四个网损最小、电压偏差最小、无功补偿经济性最优、多目标加权。IEEE33这类标准测试系统里最常用的就是网损最小因为它是无功优化的初心物理意义清楚、便于和其他文献对比。如果你想顺便改善电压可以把电压偏差作为一个惩罚项加进去但要注意权重设置这个我后面专门讲。控制手段方面模型里最方便处理的是连续无功补偿容量可以把它理解成逆变器型分布式电源的无功输出或者SVC/SVG这类连续无级补偿装置的输出值。如果要模拟传统的分组投切电容器就把补偿变量改成整数变量问题会从SOCP变成混合整数二阶锥规划MISOCPYALMIP加Gurobi或CPLEX能直接解只是求解时间会长一些。OLTC有载调压和无功优化联动在配网里也很常见但IEEE33基础版本里一般不建模变压器档位先把补偿设备做好再扩展也不迟。约束条件分四类。第一类是潮流平衡约束这是所有电力系统优化问题的宪法不能违背第二类是节点电压上下限国内工程上通常按0.931.07pu学术论文里更常用0.951.05pu第三类是支路电流上限给线路容量一个保护第四类是补偿设备自身的容量限制。约束不是越多越好尤其是IEEE33这种测试系统加太多约束反而会掩盖算法本身的对比价值新手容易犯的错就是什么约束都往里塞结果既难调参又看不出性能差异。1.3 IEEE33节点为什么它是配网优化的标准考场IEEE33节点系统来自Baran和Wu在1989年发表的论文原始用途就是配电网无功优化和网络重构研究。系统电压等级12.66kV33个节点、32条支路总负荷3715kW加2300kvar呈辐射状结构从变电站出发分出三条主馈线再加若干分支拓扑一看就很有代表性。这个系统最讨喜的地方在于规模不大做验证速度快但它有真实的电压越限问题——不装任何补偿设备时末端最低电压只有0.913pu左右全网有功损耗约202.7kW降损空间非常明显。正因为有这些病灶IEEE33成了检验配网优化算法的标准考场。你优化前后一对比电压抬升了多少、损耗降了多少数字摆在那里效果一目了然。再加上网上数据容易找新算法拿它先跑一遍几乎成了行规。做这个项目时我建议你把支路参数和负荷数据单独存成一个m文件或CSV别每次拷贝粘贴——后面调参、做对比实验会频繁用数据文件独立出来能省非常多事。2. 为什么偏偏是二阶锥规划2.1 数学上先看清楚潮流的killer非凸性配电网的潮流方程里藏着一个非常不友好的东西支路电流的平方等于视在功率的平方除以电压幅值的平方也就是 I² (P² Q²)/V²。这个关系里P、Q、V都是变量P²Q²除以V²是一个非凸约束直接把整个优化问题变成了非凸问题。非凸意味着什么意味着你没法保证找到的是全局最优解梯度类算法经常陷在局部最优里出不来初值选得不好跑几次得到的结果可能相差很大。这就是为什么以前大家做配网无功优化要么用内点法硬解NLP要么上粒子群、遗传算法这些启发式算法。内点法看着优雅但非凸NLP对初值极度敏感同一个题设换一个初值可能收敛到完全不同的局部解启发式算法虽然能全局搜索但计算量大、稳定性差而且没有任何最优性保证论文里用一用没问题工程上可不敢这么玩。2.2 DistFlow加锥松弛把硬骨头变软既然潮流方程非凸那就想办法让它变凸。配电网是辐射状结构这给了我们一个独门武器——DistFlow模型。DistFlow把每个节点的功率平衡和电压关系写成递归形式对每条支路(i,j)有v_j v_i - 2(r_ij·P_ij x_ij·Q_ij) (r_ij² x_ij²)·l_ij其中v_i是节点i电压幅值的平方l_ij是支路电流幅值的平方。这个方程本身是线性的真正非凸的是l_ij和v_i、P_ij、Q_ij之间的关系本来 l_ij·v_i P_ij² Q_ij² 是个等式约束现在把等号放松成不等号让 l_ij·v_i ≥ P_ij² Q_ij²再做一次等价变形就得到标准的二阶锥约束|| [2P_ij; 2Q_ij; v_i - l_ij] ||₂ ≤ v_i l_ij你可能想问把等式放松成不等式解还靠谱吗这就涉及SOCP最核心的工程直觉目标函数是最小化网损 sum(r_ij·l_ij)而电阻r_ij恒为正优化器会拼命把l_ij往小里压所以在最优解处这个锥松弛通常会被拉紧到恰好取等号也就是说放松前后的解是一样的。这就叫精确松弛。数学上Farivar、Low等人在2011到2013年的一系列工作证明了对于辐射状配电网在总负荷有上界、目标函数关于潮流单调递增等温和条件下SOCP松弛是精确的。实际做IEEE33优化时松弛间隙小到可以忽略。这一点很有意思——本来是个让人头疼的非凸问题靠着网络结构直接降维成凸问题求全局最优从一个奢望变成了多项式时间内轻松搞定的事。2.3 与粒子群、传统内点法的横向对比方法能否保证全局最优计算效率处理离散设备实现难度非线性规划内点法否易陷局部最优中等需特殊处理中粒子群/遗传算法否无最优性证明慢依赖种群参数天然支持较低二阶锥规划SOCP是多项式时间很快连续变量天然支持中混合整数二阶锥MISOCP是较慢支持分组电容器中高在IEEE33这种几十个变量的规模下SOCP的求解时间通常在秒级以内。我第一次用粒子群算同样的问题种群100、迭代200次跑了很久才出一个还不一定能复现的结果换成YALMIP加Cplex几秒钟就出全局最优。后来我也见过不少同学还在用遗传算法算IEEE33的无功优化不是说不能用但在需要反复对比、批量仿真的场景下SOCP的优势是压倒性的。3. 建模仿真前要先搞定的参数细节3.1 基准值怎么选一个小数点引发的求解器血案做IEEE33仿真最容易翻车的不是算法而是基准值。这个系统原始数据是12.66kV支路阻抗单位是欧姆负荷单位是kW/kvar不标幺直接建模方程里会出现0.0922和420这种数量级差几千倍的数求解器要么罢工要么给出精度极差的解。通常做法是取V_base12.66kV、S_base10MVA这样阻抗基准值Z_baseV_base²/S_base160.2756Ω支路阻抗标幺值落在0.001到0.02之间负荷标幺值落在0.004到0.06之间各变量的数量级基本一致数值条件非常好。有个细节值得提醒有些同学图省事取S_base100MVA算出来的支路阻抗就变成1e-4量级电压方程里2(rPxQ)这种项只有1e-5再叠加上二阶锥约束里的微小数字求解器的数值容差一挑动结果就容易出现假收敛。如果你用的是老版本求解器这种问题更明显。我直接给结论IEEE33就用10MVA基准值少走弯路。另外要注意单位换算。负荷表给的是kW和kvar标幺化要除以S_base*1000因为10MVA等于10000kVA不是直接除以10。好多同学第一步就在这错了后面电压约束再怎么调都无济于事检查半天最后发现负荷标幺值整体大了1000倍那种心情我太懂了。3.2 候选补偿节点怎么挑补偿点位选择直接决定优化效果。IEEE33里优先考虑的节点有三类一是末端节点因为末端电压最低、边际效益最大典型的就是18、22、33这三条馈线的末端二是重负荷节点比如24、25号节点负荷各420kW是整个系统最重的负荷群旁边放补偿效果立竿见影三是带重无功负荷的特殊节点比如30号节点带600kvar的无功负荷电压掉得特别明显也是热门候选。我比较推荐的最基础配置是[18; 22; 25; 30; 33]这五个候选点补偿容量上限分别设为300、300、600、600、300kvar。这个配置之所以这么定原因很实在三个末端点各给300kvar把最典型的低电压问题压住25和30号这两个重负荷点给到600kvar因为它们的负荷体量在那摆着给小了没效果给大了又会反向送无功。如果你想做得更深入可以把候选点扩到16、18、22、25、30、33甚至把每个候选点的容量上限也作为一个变量去做灵敏度分析看哪个点对降损贡献最大——这其实就是配电网无功规划里的选址定容问题了。3.3 目标函数和约束的取值范围建议目标函数我建议先只写网损最小obj sum(r .* Iij2)。这个目标对SOCP最友好因为电阻恒正能保证锥松弛被拉紧。如果你想兼顾电压可以用 obj sum(r.Iij2) lambdasum((Vi2-1).^2)。lambda的取值要小心IEEE33的网损标幺值在0.01到0.02量级而电压偏差平方和可能到0.01到0.05直接相加会让电压项喧宾夺主。我试下来lambda在0.1到1之间比较合适而且建议先把lambda0的纯网损结果跑出来再逐渐增大lambda观察网损和电压两个指标此消彼长的关系最后再定权重。电压约束方面IEEE33不加补偿时最低电压约0.913所以如果你把下限设成0.95那么优化前的工况在这个约束下就是不可行的。这不是bug反而说明问题有优化空间。做对比实验时基准工况无补偿的电压下限要放宽到0.88到0.90或者直接用解潮流的工具算出真实基态优化工况再用0.95的合格下限。这个细节很多新手没注意把基准工况和优化工况用同一套电压约束结果曲线一出来对比就失去了说服力。补偿容量的上限按节点负荷的无功分量估算最合理。经验值是取该节点Q_load的20%到50%或者参照实际工程里单台电容器的额定容量按30到100kvar的档位去凑。25号节点负荷重给到600kvar很合理30号节点本身带600kvar的无功负荷给600kvar也符合工程直觉末端节点给300kvar就够了。给太多会让电压抬过头、容性无功反向流动反而白白增加损耗。4. Matlab YALMIP 核心代码实现4.1 环境配置与求解器选型跑SOCPMatlab本身不算关键关键是配套的建模工具和求解器。我强烈建议用YALMIP做建模层它把变量声明、约束拼接、求解器调用统一起来比手写Cplex API或直接用CVX都省心。求解器方面优先Gurobi或Cplex它们都原生支持二阶锥约束且能解MISOCP对IEEE33这种规模几秒钟就完成求解。如果学校没有商业许可可以用免费开源的SeDuMi或SDPT3解连续SOCP完全没有问题只是速度稍慢、不能直接处理整数变量。这里特别说一个很多人问的点装了Gurobi还能不能装Cplex答案是完全可以8这两个工具在MATLAB路径里可以和平共处YALMIP会全部识别运行时用sdpsettings指定具体用哪个求解器就行。我的习惯是先用SeDuMi跑通逻辑再换Gurobi跑大规模和整数版本两者互不干扰。装好后可以用yalmiptest这个命令检查YALMIP到底识别了哪些求解器如果有显示No suitable solver多半是路径没加对。4.2 IEEE33完整数据和预处理先准备数据文件。我建议把基础数据独立存成一个函数文件case33_ieee.m返回基准电压、基准容量、支路表和负荷表四个量之后所有脚本都从它读取逻辑清晰、可复用。下面是完整数据摘自经典IEEE33节点算例可以直接保存使用。function [Vb, Sb, br, ld] case33_ieee % IEEE33节点配电网基础数据 % 输出: Vb-基准电压(kV), Sb-基准容量(MVA), br-支路表, ld-负荷表(kW/kvar) Vb 12.66; Sb 10; br [ 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 ]; ld [ 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 ]; end数据文件准备好之后在主脚本里调用然后做标幺化。这块代码并不长但每行都有讲究[Vb, Sb, br, ld] case33_ieee; nbus 33; nbranch size(br, 1); Zb Vb^2 / Sb; % 阻抗基准值单位欧姆 r br(:,3) / Zb; % 支路电阻标幺 x br(:,4) / Zb; % 支路电抗标幺 pL ld(:,2) / (Sb*1000); % 有功负荷标幺kW / 10000 qL ld(:,3) / (Sb*1000); % 无功负荷标幺kvar / 10000这里最容易错的就是最后两行。Sb单位是MVASb*1000才是kVA负荷的kW和kvar必须除以这个数才能得到正确的标幺值。我把这行写在最醒目的位置就是希望正在调试的你少走一次弯路。4.3 变量定义与约束构建代码逐段拆解变量定义直接用YALMIP的sdpvar。支路有功Pij、支路无功Qij、支路电流平方Iij2都是32维向量节点电压平方Vi2是33维向量无功补偿qC是33维向量其中只有候选节点才能真正出力。目标函数用网损加电压偏差惩罚项的形式后面你可以自己调lambda。% 决策变量 Pij sdpvar(nbranch, 1); % 支路有功 Qij sdpvar(nbranch, 1); % 支路无功 Iij2 sdpvar(nbranch, 1); % 支路电流平方 Vi2 sdpvar(nbus, 1); % 节点电压平方 qC sdpvar(nbus, 1); % 无功补偿容性为正 % 约束容器 cons []; % 根节点电压等于1 cons [cons, Vi2(1) 1]; % 电压幅值上下限 cons [cons, 0.95^2 Vi2 1.05^2];电压上下限这个地方优化工况用0.95到1.05这是学术论文最常用的区间。如果你做基准工况对比记得把下限放宽否则无补偿的基态在模型下不可行。接下来是最核心的部分支路电压方程和二阶锥约束。对每条支路先写电压降落关系再写SOC约束。这里的循环清晰直观每条支路两行约束不容易出错。% 支路方程 二阶锥约束 for k 1:nbranch i br(k,1); % 首端节点 j br(k,2); % 末端节点 % 电压平方关系 cons [cons, Vi2(j) Vi2(i) - 2*(r(k)*Pij(k) x(k)*Qij(k)) ... (r(k)^2 x(k)^2)*Iij2(k)]; % 二阶锥松弛约束 cons [cons, norm([2*Pij(k); 2*Qij(k); Vi2(i) - Iij2(k)], 2) ... Vi2(i) Iij2(k)]; enddistflow模型的这一行是支路潮流的核心它描述的物理过程是末端电压平方 首端电压平方 - 线路压降 线路损耗项。二阶锥约束这一行就是我前面讲的把非凸等式放松成凸不等式的关键一步YALMIP会把norm(x,2)y这种形式自动识别成二阶锥交给求解器时不必手动转换格式。节点功率平衡约束是另一个关键。对每个非根节点流入该节点的功率要等于流出功率加上净负荷。这里我用了find来定位每个节点的父支路和子支路循环清晰也方便你以后改成其他拓扑结构。% 节点功率平衡1号节点是平衡节点不写 for j 2:nbus par find(br(:,2) j); % 以j为末端的分支功率流入 chd find(br(:,1) j); % 以j为首端的分支功率流出 pbal 0; qbal 0; for kk par pbal pbal (Pij(kk) - r(kk)*Iij2(kk)); qbal qbal (Qij(kk) - x(kk)*Iij2(kk)); end for kk chd pbal pbal - Pij(kk); qbal qbal - Qij(kk); end cons [cons, pbal pL(j)]; % 有功平衡 cons [cons, qbal qC(j) qL(j)]; % 无功平衡 end无功平衡式的符号要特别注意qbal是支路送来的无功qC是本地补偿注入的容性无功为正两者之和等于负荷的无功消耗qL。如果把qC的符号搞反求解器要么给出明显不合理的补偿结果要么直接报不可行。补偿变量加容量约束也很简单。用零向量初始化上限只给候选节点赋值。这种写法比直接写五个单独约束清晰得多后面想调整候选节点也方便。% 无功补偿容量限制 capNodes [18; 22; 25; 30; 33]; qCapMax zeros(nbus, 1); qCapMax(capNodes) [0.03; 0.03; 0.06; 0.06; 0.03]; % 对应kvar数 cons [cons, 0 qC qCapMax];比如0.03的标幺值在10MVA基准下就是300kvar0.06就是600kvar。这个容量上限给的依据就是前面说的节点负荷体量和实际工程中电容器的常见容量档位。4.4 求解与结果导出目标函数和求解部分非常短。先把目标函数定义成网损加惩罚项然后调optimize函数。sdpsettings里指定求解器我一般用cplex或gurobi如果要用免费求解器就改solver为sedumi。% 目标函数网损最小 轻量电压偏差惩罚 lambda 0.5; obj sum(r .* Iij2) lambda * sum((Vi2 - 1).^2); % 求解 options sdpsettings(solver, cplex, verbose, 2); sol optimize(cons, obj, options); % 结果检查 if sol.problem 0 fprintf(求解成功\n); else yalmiperror(sol.problem); end求解完成之后把结果从标幺值还原成工程单位。网损是求和后乘以Sb再乘以1000得到kW补偿量乘以Sb再乘以1000得到kvar电压幅值直接对Vi2开方。% 结果换算 V_opt sqrt(value(Vi2)); loss_opt sum(r .* value(Iij2)) * Sb * 1000; % kW Qcomp_opt value(qC) * Sb * 1000; % kvar fprintf(优化后网损: %.2f kW\n, loss_opt); fprintf(最低电压: %.4f pu 节点%d\n, min(V_opt), find(V_opt min(V_opt))); for kk 1:length(capNodes) n capNodes(kk); fprintf(节点%2d 无功补偿: %.1f kvar\n, n, Qcomp_opt(n)); end到这里一个完整可跑的SOCP无功优化程序就成型了。整个脚本从数据准备到结果输出大约120行调试时间通常控制在几分钟内。如果你跑出来的结果跟常见文献对不上先回去查基准值和负荷单位这两个地方错一个后面全盘皆输。5. 典型结果分析与自检方法5.1 优化前后的电压曲线对比要拿电压曲线说明问题至少要有两条线一条是优化前的基态电压一条是优化后的电压。基态电压怎么来最简单的方法是不加补偿、把电压下限放宽用同一套SOCP模型再跑一次。% 基态求解无补偿 电压下限放宽 cons_base []; % ... 复用上面的约束但去掉补偿约束或令 qC0 % 电压下限改成 0.88^2 而不是 0.95^2 % 目标函数只用网损最小 % 求解后得到 V_base跑完之后画图比较figure; plot(1:nbus, V_base, -o, LineWidth, 1.5); hold on; plot(1:nbus, V_opt, -s, LineWidth, 1.5); grid on; xlabel(节点编号); ylabel(电压幅值 (pu)); legend(优化前, 优化后, Location, southwest); ylim([0.90, 1.06]);一个典型的可复现结果是优化前最低电压在18号节点附近约0.913pu优化后全网电压全部抬升到0.95以上最低电压点的位置和幅值都明显改善电压曲线整体变得更平缓。曲线图一出来优化效果非常直观。注意别把两条线画成同一种样式加粗、加标记评审和答辩时能省不少解释时间。5.2 网损下降率怎么算才严谨网损下降率是评审最认的硬指标。计算公式很简单降损率 (优化前网损 - 优化后网损) / 优化前网损 × 100%。难点在于优化前网损必须从同一个系统的同一种运行方式下算出来不能拿别的论文的数字顶替。我自己这次跑的典型结果是优化前约202.7kW优化后在175kW左右降损率大约13.6%。不同候选节点配置、不同容量上限会得到不同数字这个浮动是正常的你要关注的是自己的对比口径是否一致。关于降损率我想多说一句别只盯一个数字。你还要顺手看两个佐证指标全网总无功补偿量是多少电压偏差从多少降到多少。如果降损率很高但补偿容量大得离谱或者电压已经顶到上限那这个结果即使数学上最优工程上也不一定可行。IEEE33只是测试系统但做测试系统的严谨习惯会直接影响你以后做真实配网项目。5.3 松弛间隙检查你的解是不是真的稳SOCP最让人担心的问题就是松弛是否精确。虽然理论上有保证但数值实现和特定目标函数下还是有必要自己检查一遍。方法是求解后把所有最优值代回锥约束看不等号两侧的差距。% 松弛间隙检查 gap_soc zeros(nbranch, 1); for k 1:nbranch i br(k, 1); lhs norm([2*value(Pij(k)); 2*value(Qij(k)); value(Vi2(i)) - value(Iij2(k))], 2); rhs value(Vi2(i)) value(Iij2(k)); gap_soc(k) rhs - lhs; end fprintf(最大松弛间隙: %.3e\n, max(gap_soc));如果这个间隙在1e-6量级甚至更小说明锥约束在最优解处基本取等号你的解是可信的。如果你改了目标函数比如把电压偏差的权重调得很大锥约束可能不再拉紧这时求解器给的最优解就不一定对应真实物理潮流了。这个检查是我的必做步骤也是我判断一个SOCP无功优化代码是否靠谱的第一道门槛。如果条件允许还可以用电力系统潮流软件做交叉验证。Matpower新版本里带了case33bw这个IEEE33数据文件你可以跑一次runpf得到基态潮流再和SOCP模型在零补偿状态下的结果对比电压和网损能对上那基本就稳了。6. 常见问题与排查技巧实录6.1 求解器报错与兼容性问题YALMIP报No suitable solver for a class of problem是最常见的错误之一。这句话的意思是你的模型里包含了当前求解器不支持的约束类型最常见的情况是你用了不支持二阶锥的求解器或者装了求解器但没加进MATLAB路径。先运行yalmiptest看识别列表如果列表里有cplex或gurobi大概率是模型里混入了YALMIP无法自动转化的约束逐个检查约束类型即可。关于多个求解器共存的问题我在4.1已经说过Gurobi和Cplex可以同时安装YALMIP会自己识别。实际中还会遇到一种尴尬就是某个求解器版本太老对MISOCP支持不完整。我遇到过Gurobi老版本报二次约束不是凸的这种误导信息其实问题出在版本太旧。解决办法是升级求解器或者在sdpsettings里强制指定另一个求解器做交叉验证。6.2 不收敛或不可行的排查清单碰到不可行问题先别急着怀疑代码。按这个顺序排查第一步把电压下限从0.95放宽到0.88看问题是否变可行——如果放宽后就可行说明补偿容量不足或候选节点不够需要增大容量或增加候选点第二步去掉二阶锥约束以外的各类附加约束比如电压上限、支路电流上限逐条加回来找到是哪个约束导致不可行第三步检查目标函数和变量是否存在维度不匹配YALMIP的sdpvar乘法很容易在无意中制造非线性项第四步用sol.info和yalmiperror输出详细诊断信息很多不可行问题在诊断信息里会直接告诉你哪个约束冲突。求解器给了可行解但结果不合理比如补偿全为负、网损比基态还高优先检查符号约定。无功补偿的符号、功率平衡方程的正负号、支路功率流向的定义三个环节任何一处反了都会出这种诡异结果。IEEE33的对称性不强节点编号顺时针看是三条馈线但支路顺序不是严格按节点编号排列的用find函数定位分支时要特别留意。6.3 新手最容易踩的5个坑第一负荷单位没转对。kW和kvar要除以S_base*1000而不是除以S_base这个错误我见得太多了它在数值上表现为所有负荷标幺值大了1000倍电压约束形同虚设。第二基准工况和优化工况的约束不一致。跑基态时用了0.95的电压下限结果基态不可行就强行加补偿去救然后两个工况混在一起对比曲线完全没有说服力。正确的是基态放宽电压下限、不加补偿优化态加补偿、用合格下限两个工况的口径分开。第三只盯着网损下降率不检查松弛间隙。我在5.3讲过目标函数换成电压偏差主导后锥约束可能不紧。不紧的SOCP解看起来是凸问题最优解实际上对应的物理潮流是虚的。每跑一版都要带回代检查。第四补偿容量上限拍脑袋。容量上限最好依据节点无功负荷量级来设一个只有90kvar负荷的末端节点你给1200kvar的补偿求解器会算出补偿容量顶到上限但电压已经超限的荒唐结果。先做一次无功负荷分布统计再定容量。第五忽略YALMIP和求解器的版本匹配。老版本YALMIP对某些SOC约束的语法支持不完整新版本求解器又可能改变了默认容差。我建议固定一套经过验证的组合比如YALMIP R2021版配Gurobi 9.x然后所有对比实验都在同一套环境里跑免得被版本差异干扰结论。最后再分享一个我自己的习惯。做这个项目时我会把每一版结果都存成带时间戳的mat或Excel文件跑完一组实验先把电压曲线和网损数字记录到一张对比表里而不是只看命令行输出的最后三行。这样做的好处是论文返修或评审质疑结论时你能快速回溯每一个数字是在什么参数下跑出来的。你是要拿这个模型替自己说话的数据和结论的源头越清楚你说的话才越硬气。SOCP这个方向看起来是公式堆出来的但真正能让你在答辩和工程评审里站稳的是这些围绕模型的验证习惯和细节控制。