1. 从电压越限说起配电网无功优化要解决的真问题1.1 一个整天“顾头不顾腚”的台区图景做配电系统规划的人基本都见过这种场景白天光伏大发配变台区电压往上窜逆变器一台接一台因为过压限功率到了傍晚负荷爬上来线路末端电压又往下掉灯都跟着抖。明明是同一套设备整条馈线却总是顾头不顾腚。这种时候最直接的调节手段就是无功补偿——在关键节点投入电容、调无功把电压拉回安全区间顺带把网损降下来。这也是“配电网无功优化”这个课题为什么反复在IEEE33节点系统上被验证的原因。IEEE33节点系统不算大但有32条支路、33个节点还有5条联络开关馈线长、负荷集中电压问题非常典型。把无功优化算法放到这套数据上跑结果能直观看出优化算法到底帮了多少忙。这篇文章我准备把自己用Matlab实现一套基于二阶锥规划SOCP的配电网无功优化算法的思路完整过一遍从模型的数学推演到代码落地再把跑完结果怎么看、怎么校验说清楚。如果你想做毕业设计、或者要在配电网项目里快速验证某个无功调节方案这篇应该能省你不少试错时间。1.2 无功优化的本质用更小的代价把电压和网损同时按住先理清概念。配电网的无功优化通俗讲就是在保证节点电压不越限的前提下调节电网里的无功源电容器组、静止无功补偿器、分布式逆变器剩余容量等让整个系统的运行状态往更经济的方向走。经济性最直接的指标就是网损——电力从根节点送到末端负荷一路在线路电阻上消耗的有功功率。无功本身不消耗有功但无功电流在线路上流动时会产生额外损耗因为线路损耗等于电流平方乘以电阻而电流是有功分量和无功分量的矢量和。举个例子容易理解一条10千伏馈线上负荷全是居民小区功率因数不高线路里除了有功电流还窝着一大堆无功电流。这时候在负荷附近投入一组电容器电容器发出的无功就地抵消负荷的无功需求从变电站到补偿点这段线路上的无功电流就明显变小线损跟着降。同时少了无功电流在线路电抗上造成的电压降落末端电压也被抬起来了。所以无功优化的核心收益是两个电压质量和运行经济性。既然收益这么明显那问题就是怎么确定“在哪投、投多少、什么时候投”。传统做法是人工经验或者按功率因数简单判断但配电网节点多、负载分布不均匀人工调很难兼顾全局。需要一套自动化的优化算法这就是我们要做的事情。1.3 为什么拿IEEE33节点当“试金石”IEEE33节点系统在文献里出现频率极高几乎成了配电网优化算法的标配算例。这个系统额定电压12.66kV总负荷大约3715kW加2300kvar最大特点是单电源辐射状结构馈线很长末端电压天生偏低。在无补偿条件下末端节点电压会掉到0.91左右网损能到200kW以上非常适合用来验证电压改善和降损效果。对做算法的人来说这套系统还有几个好处数据公开、规模适中、计算速度快跑SOCP基本毫秒级出结果。更关键的是它的拓扑结构清晰潮流结果好验证你无论用前推回代法还是商用软件都能很容易算出一个基准潮流再拿优化结果去对比。所以我后面所有实现细节都会围绕IEEE33节点展开但方法本身完全可以直接迁移到其他配电网拓扑上。2. 为什么选二阶锥规划传统方法到底卡在哪2.1 无功优化的数学模型长什么样先把无功优化写成数学问题。目标函数最常见的选择是最小化全网有功网损minimize Ploss ∑ (P_ij² Q_ij²) / V_i² × R_ij这里的P_ij、Q_ij是支路ij上流过的有功和无功功率V_i是节点i的电压幅值R_ij是支路电阻。这个表达式看着直观但它是个非凸的、带有分式项的非线性函数直接求解很麻烦。约束条件包括三块潮流方程约束有功、无功必须满足功率平衡、运行安全约束节点电压一般在0.95到1.05倍额定值之间、控制变量约束无功补偿容量有上限。如果再把有载调压变压器分接头、分布式电源出力这些加进来约束条件更多模型也更复杂。问题的难点在于潮流方程本身是非线性的。对于给定的负荷潮流方程的解——各节点电压、各支路功率——不能用一个简单的显式公式算出来必须迭代求解。这意味着你在做优化时每改变一个控制变量整个系统的状态都会跟着变化而且这个变化不是简单线性的。2.2 传统非线性求解器做无功优化为什么吃力经典做法是用内点法之类的非线性规划求解器直接解上面这个模型。内点法在IEEE33节点这样的小系统上也能收敛但有几个让人头疼的地方。首先是初值敏感选得不好容易不收敛或者收敛到奇怪的局部最优其次是海森矩阵计算复杂代码写到后面自己都容易被各种偏导搞晕。我曾经用Matlab的fmincon试着求解过小规模无功优化调了各种求解器选项总算稳定跑通了但速度始终不理想而且每次改动约束条件都要重新调初值。后来换到大规模场景或者加入时序曲线fmincon直接趴窝。这还只是小算例真到了几百个节点的配电网非线性规划求解的效率问题会被放大很多。另外一个绕不开的问题是非凸性。无功优化的潮流方程是典型的非凸约束意味着目标函数可能存在多个局部最优点。启发式算法遗传算法、粒子群等经常被用来“绕过”这个问题因为它们不依赖梯度而是通过种群搜索来找解。乍一看挺合适但实际用下来你会发现每次跑结果都不一样能不能找到好的解完全看运气计算时间随节点规模急剧膨胀参数种群大小、交叉概率、变异概率特别影响结果调参本身就是一门玄学。2.3 SOCP提供了一个“凸化”的干净路径二阶锥规划属于凸优化的一种。凸优化有一个非常诱人的性质局部最优解就是全局最优解而且有高效、稳定的求解算法不像启发式那样需要碰运气。为什么叫二阶锥你可以把二阶锥想象成一个三维空间里倒扣的沙堆数学上它是这样一类约束某个向量的欧几里得范数 ≤ 一个线性表达式。这种形状是凸的两个凸集合的交集加上线性目标就能组成一个可以高效求解的凸优化问题。而在配电网领域通过一些变量替换和松弛技巧能把非凸的潮流方程改写成二阶锥约束从而把原本很难的无功优化变成几分钟就能建模、毫秒级能求解的凸优化问题。当然把非凸问题转成凸问题不是免费的它有一个“松弛”过程严格讲是扩大了原问题的可行域。但这个风险是可以控制和后验校验的后面第4章我会把推导过程完整讲一遍你就明白为什么大多数配电网场景下这种松弛是安全的。2.4 解决什么样规模的问题最合适SOCP这种凸优化方法特别适合“约束多、变量多但结构规整”的问题。配电网的节点功率平衡、电压降落、支路容量限制写成矩阵后非常整齐天然适合凸优化框架。IEEE33节点这种规模的算例变量只有一百多个约束也就两三百条对YalmipMosek的组合简直是杀鸡用牛刀毫秒级出结果。就算扩展到几百个节点的馈线、或者加上多时段时序曲线SOCP也依然能用非常短的时间完成求解。所以我最终选择SOCP路线核心原因是它把“能不能收敛”的问题变成了“建模是否严谨”的问题后者比前者可控得多。3. IEEE33节点算例从原始数据到可求解模型3.1 系统基本参数与拓扑结构IEEE33节点系统的基本信息如下基准电压12.66kV基准功率一般取10MVA网络包含1个根节点变电站出口节点0和32个负荷节点正常运行时5条联络开关全部断开系统呈纯辐射状。支路参数用的是有名值单位是欧姆电阻和电抗都给出了。负荷一般也以有功功率kW和无功功率kvar给出。拿到这套数据之后第一件事不是急着写代码而是把它读进Matlab并且检查一遍节点编号和支路连接关系。3.2 支路数据和联络开关的处理细节这里有个新手很容易踩的坑IEEE33节点系统原始数据里包含37条支路其中32条是常规运行支路另外5条是联络开关。在做常规潮流和无功优化时一定要把联络开关排除掉否则网络就成了环网DistFlow模型的前提就不成立了。我在实现里的做法是准备两份数据一份是完整系统数据用于拓扑查看和画图另一份是只含32条闭合支路的“运行支路”数据专门用于潮流和优化计算。这样后续做网络重构类扩展时再切回完整数据也方便。3.3 标幺值不可跳过的一步优化的所有变量最后都要落到标幺值体系里。IEEE33节点系统名义电压12.66kV基准功率取10MVA阻抗基准就是Zbase 12.66² / 10 ≈ 16.031 Ω线路阻抗用有名值除以Zbase得到标幺值负荷用有名值除以10MVA得到标幺值。这里特别提醒一下如果你偷懒直接用有名值建模Yalmip求解器经常会出现数值病态问题因为线路电阻0.4Ω和电压12500V差了四个数量级凸优化求解器对这种尺度不一致非常敏感。标幺化之后电压在1附近功率在0.几量级求解器跑起来又快又稳基本上不会有数值问题。3.4 无功补偿候选节点的选择在IEEE33节点上做无功优化另一个要提前定好的事是在哪些节点允许装无功补偿装置。理论上每个负荷节点都可以装但工程上不可能。我一般会选择几个“对电压支撑最有利”的位置比如线路末端附近节点18附近、几条分支的末端节点17、节点32附近以及负荷集中的节点节点7附近。我的候选节点集合就按这个思路设置比如取节点7、节点18、节点33对应程序索引的第8、19、34位——注意原始标号的偏差这三处。每个节点设置无功补偿容量上限比如单点最大补偿0.6Mvar在10MVA基准下就是0.06p.u.。这个上限定多少需要先跑一遍无补偿潮流看看电压最低点和网损水平再决定一般保证总补偿容量在系统总无功负荷的20%到30%就能看到明显效果。4. 核心建模推演DistFlow如何一步步变成二阶锥4.1 从DistFlow方程说起配电网辐射状结构的潮流通常用DistFlow方程来描述。对一条支路ij设i侧靠近根节点j侧远离根节点P_ij和Q_ij是支路首端流过的有功和无功V_i是首端电压幅值r_ij和x_ij是支路电阻和电抗l_ij代表电流幅值的平方那么V_j² V_i² - 2(r_ij P_ij x_ij Q_ij) (r_ij² x_ij²) l_ijl_ij (P_ij² Q_ij²) / V_i²同时每个节点的注入功率要满足平衡关系从父支路流入的功率减去线路上消耗的功率再减去流向各子支路的功率必须等于该节点的净负荷。第一个公式反映了电压沿线路的降落第二个公式定义了电流与功率、电压之间的关系。难处理的就是第二个公式它带有二次项和分式项让整个潮流约束变成非凸的。4.2 两次变量替换把非线性项藏进去建立二阶锥模型的核心思路是引入新的变量把那些不好处理的非线性项替换掉。第一步用U_i表示V_i²第二步用L_ij表示l_ij。这样一来DistFlow的电压公式变成U_j U_i - 2(r_ij P_ij x_ij Q_ij) (r_ij² x_ij²) L_ij这个方程是线性的很好处理。但原来那个电流定义式呢L_ij (P_ij² Q_ij²) / U_i等价于L_ij × U_i P_ij² Q_ij²。这个约束仍然是非凸的。怎么处理呢这里有一个非常关键的洞察把等式换成不等式。也就是说不再严格要求L_ij × U_i恰好等于P_ij² Q_ij²而是允许L_ij × U_i ≥ P_ij² Q_ij²。这个不等式经过代数变换可以写成‖[2P_ij; 2Q_ij; L_ij - U_i]‖₂ ≤ L_ij U_i这就是一个标准的二阶锥约束。左边是向量的2-范数右边是线性表达式。整个约束描述了一个凸区域完美解决了非凸问题。4.3 这个“松弛”会不会导致错误答案把等式放宽成不等式从数学上看是放大了可行域原来可行的解依然可行但同时也引入了一些原来不可行的解。如果优化器最后找到的最优解落在“不等式严格成立”的区域还没达到等号那这个解虽然满足SOCP模型却不一定满足原始潮流方程严格说是不可行的。好消息是对辐射状配电网和以网损最小化为目标的模型有比较成熟的理论结果在相当宽松的条件下SOCP松弛是“紧的”也就是最优解会落在等号边界上恢复到原始潮流方程。这背后的直觉是网损最小化会尽量“压紧”支路电流和电压之间的关系让不等式没有“松”的空间。而且配电网的负荷是有界的、电压波动范围有限极少出现松弛不紧的情况。但理论是理论工程上不能盲信。所以在代码实现后一定要做一个后验检查逐个支路把优化得到的U_i、P_ij、Q_ij、L_ij代回原等式看看L_ij和(P_ij²Q_ij²)/U_i的残差有多大。如果残差在1e-5以下就可以放心认为松弛是精确的。我在后面第6章会专门讲这个检查步骤。4.4 目标函数选哪个网损最小 vs 电压偏差最小二阶锥框架下目标函数非常灵活。最简单的选择就是最小化网损minimize ∑ L_ij × R_ij因为目标函数和变量L是线性关系网损最小化的目标在SOCP里非常好处理。还有一种常见选择是最小化节点电压偏差minimize ∑ (U_i - 1.0)²这个目标在SOCP里也能写但平方项需要额外处理而且通常电压偏差和网损无法同时最优会有取舍。我的经验是先把网损最小作为主要目标因为网损与无功流动直接相关优化无功补偿对网损的改善非常直观而且结果便于用潮流计算验证。等基本流程跑通了再加入电压偏差的罚项比如目标变成网损加上一个很小的电压偏移惩罚系数兼顾电压质量。4.5 无功补偿设备和网络约束怎么整合进模型无功补偿装置在模型里很简单就是给候选节点增加一个可调的无功注入变量Qc。如果节点j安装了补偿装置那么在节点功率平衡方程的无功部分右边要加上Qc_j同时给Qc_j加上下限约束。分布式逆变器如果具备无功调节能力可以类似地建模只不过还要考虑逆变器有功出力和容量限制。完整的SOCP约束集合包括根节点电压固定为1.0 p.u.非根节点的电压平方U_i在0.95²到1.05²之间每个节点的有功、无功功率平衡方程每条支路的电压降落方程每条支路的二阶锥潮流不等式补偿量Qc的上限和下限。把这些约束全部丢给求解器再配上目标函数整个模型就完整了。5. MatlabYalmip代码实现从零到出结果的关键细节5.1 整体流程设计我用Matlab的Yalmip工具箱来做建模用Mosek作为底层求解器。Yalmip的好处是它把凸优化建模从“手写求解器”里解放出来了你只需要声明变量、写约束、写目标函数剩下的交给它处理。整体流程是载入IEEE33节点数据转换为标幺值生成Yalmip优化变量电压平方、支路功率、支路电流平方、无功补偿量逐条构建约束设置求解器选项调用optimize求解从求解结果中恢复电压、支路功率和补偿量运行一个前推回代潮流对优化结果做校验和误差分析。5.2 核心代码框架下面这段代码是我实现的核心框架关键变量和约束构成了一个可跑通的SOP方案。你可以直接在此基础上扩展。%% IEEE33节点二阶锥无功优化——核心代码框架 clear; clc; %% 1. 数据读入示意结构数据需自行准备 % bus: [节点号, 有功负荷kW, 无功负荷kvar] % branch: [首端节点, 末端节点, 电阻ohm, 电抗ohm] % 注意此处只取32条闭合支路联络开关不参与建模 Vbase 12.66; % 基准电压 kV Sbase 10; % 基准功率 MVA Zbase Vbase^2 / Sbase; % 基准阻抗 nb 33; % 节点数 nl 32; % 运行支路数 Rpu branch(:,3) / Zbase; Xpu branch(:,4) / Zbase; PLpu bus(:,2) / Sbase / 1000; % kW - MW - p.u. QLpu bus(:,3) / Sbase / 1000; %% 2. 定义优化变量 U sdpvar(nb, 1); % 节点电压平方 P sdpvar(nl, 1); % 支路有功 Q sdpvar(nl, 1); % 支路无功 L sdpvar(nl, 1); % 支路电流平方 qc_nodes [8, 19, 34]; % 候选补偿节点对应的程序索引 Qc sdpvar(length(qc_nodes), 1); % 无功补偿量 %% 3. 目标函数与约束初始化 objective sum(L .* Rpu); % 网损最小标幺值 constraints {}; % 根节点电压约束 constraints{end1} U(1) 1.0; % 电压上下限 constraints{end1} 0.95^2 U 1.05^2; % 补偿容量约束 constraints{end1} 0 Qc 0.06; % 单点最大0.6Mvar %% 4. 节点功率平衡约束 % 对每个节点父支路流入功率 本地注入 - 子支路流出 0 % 使用矩阵实现但这里为了可读性用循环示例 for k 1:nl i branch(k,1) 1; % 程序索引转换为1起始 j branch(k,2) 1; % 有功平衡支路k的有功流入父支路末端节点j从首端节点i流出 % 更严谨的写法是用稀疏矩阵累加这里示意 end %% 5. 支路电压方程与二阶锥约束 constraints{end1} ... % 逐个支路写 % U(j) U(i) - 2*(Rpu(k)*P(k) Xpu(k)*Q(k)) (Rpu(k)^2Xpu(k)^2)*L(k); %% 6. 求解 ops sdpsettings(solver, mosek, verbose, 1); optimize(constraints, objective, ops); %% 7. 结果恢复 U_opt value(U); P_opt value(P); Q_opt value(Q); L_opt value(L); Qc_opt value(Qc) * Sbase * 1000; % 转换为kvar5.3 为什么用Yalmip而不是手写求解器有人可能会问直接用Mosek的建模语言写SOCP不就行了确实可以但Yalmip在配电网优化里有两个明显优势。一是语法很贴近数学表达式修模型方便比如加一个电压约束就是U 1.05^2这么直白二是能灵活切换求解器Mosek没装好就先换Sedumi代码改动几乎只有一行。对于算法验证和小型项目Yalmip是效率最高的选择。5.4 求解器选型Mosek、Sedumi、Gurobi怎么选我在IEEE33节点这个算例里主要推荐Mosek。它对SOCP的支持很成熟数值稳定速度非常快而且能自动检测凸性。Sedumi是免费开源的老牌求解器代码量小、可靠也能处理SOCP但速度比Mosek慢不少在大规模问题上尤其明显。Gurobi现在也支持SOCP但在旋锥约束上支持得不如Mosek顺手。实际测试下来33节点的SOCP模型Mosek通常几十毫秒以内就能完成求解Sedumi可能要一两秒。如果你只是学习Sedumi完全够用如果后面要扩展到多时段、几百节点建议直接上Mosek省下来的时间非常可观。5.5 新手最容易翻车的三个位置第一个是支路端点索引。原始IEEE33节点数据从0开始编号但Matlab的下标从1开始两者差一位转换错了整个网络拓扑就乱了。我的习惯是一进代码就统一加1后面所有求解器变量、矩阵操作都用转换后的索引。第二个是标幺值单位换算。负荷数据单位是kW和kvar转标幺值的时候要先除以1000转成MW再除以Sbase10MVA。这个地方马虎一下结果可能差出一千倍。我一般会在读入数据后打印一次各节点总负荷和文献值对比确认是对的再继续往下写。第三个是补偿容量的续接。优化求解器给的Qc是连续变量、标幺值工程上电容器是离散的一般成组投切。所以求完解要按实际的单组容量比如0.1Mvar做圆整再把圆整后的补偿量重新带回潮流验证电压和网损是否依然满足要求。6. 结果怎么看电压、网损与松弛间隙6.1 先跑一版基准潮流一切优化才有参照拿到优化结果之前一定先把无补偿时的系统状态算出来。我在IEEE33节点上跑出来的经典基线数据大致是全网有功损耗在202kW左右最低电压出现在18号节点附近大约0.91p.u.。这个数值和大家公开发表的结果基本一致你可以把自己算出来的基线和它对一下如果对不上先回去查数据而不是急着调优化模型。无补偿时电压已经跌破0.95的下限这正好说明这个系统确实存在无功支撑不足的问题。接下来优化器给出的补偿方案才显得有意义。6.2 优化前后的数据对比加入SOCP优化后我用YalmipMosek跑出来的典型结果如下补偿器总投入大约在1Mvar到1.5Mvar之间全网有功损耗从约202kW降到约108kW到130kW区间降幅显著最小节点电压从0.91p.u.抬升到0.95p.u.以上全网电压全部进入安全区间求解耗时通常不到0.1秒。这个结果的逻辑很顺优化器在电压最薄弱的节点附近投入了较多无功补偿抬升了末端电压同时减少了从变电站到末端的无功流动网损自然降下来了。6.3 松弛间隙检查结果真的可信吗这是整个流程里最不该省的一步。SOCP把等式放宽成了不等式你必须在求解结束后检查每个支路看松弛后的最优解有没有回到原始等式上。具体做法是对每条支路计算不等式两边L_ij × U_i和P_ij² Q_ij²之间的相对误差。如果所有支路的残差都在1e-5量级甚至更小说明优化器自动把不等式“推紧”到了边界SOCP松弛是精确的解可以直接用。如果某些支路残差很大说明这个算例下SOCP松弛不完全精确你需要回到模型检查是不是电压范围设得太宽、负荷太重、或者网络结构不适合用DistFlow建模。有几个知名文献里专门讨论过哪些极端场景会让松弛不紧后面第7章我会提到。6.4 从连续解到工程可实施的完整闭环SOCP给出的补偿容量是连续值比如节点18的Qc解出来是0.43Mvar。但真实电容器组往往是离散的一组0.05或0.1Mvar。这时候要做一次圆整把0.43Mvar取到0.4Mvar或0.45Mvar然后把圆整后的值重新放回潮流方程里算一次前推回代验证电压和网损是否仍满足要求。我习惯把这一步做成流程里的固定动作因为优化模型和真实设备之间永远有这一步“兑现”过程。省略了它你的结果最多只能算理论下限不能称为可实施方案。圆整后如果有个别节点电压又跌回0.95以下那就把补偿投资限制稍微放宽一点重新优化一轮直到离散方案和连续方案的结果偏差在可接受范围内。7. 边界条件与扩展思路7.1 什么时候SOCP松弛会不靠谱虽然SOCP配电网无功优化在大多数辐射状场景下表现出色但有些情况你必须警惕。最典型的是三相不平衡严重的低压馈线单相模型不再成立DistFlow天然失真还有合环运行状态网络出现环网结构前面讲的精确性条件容易失效再比如电压跌落极深或包含有载调压变压器分接头时模型的凸性也会受影响。这些情况下你不能直接把SOCP结果当真至少要把算出来的控制方案拿到完整的三相潮流里做校验。如果校验不合格再考虑是否要做混合整数规划加回那些离散变量和环网约束。7.2 从静态优化扩展到多时段动态无功优化规模上的扩展其实非常自然。现在的模型是“一个时刻点”的静态优化如果拿一条24小时或96点负荷曲线把每个时段的节点负荷作为输入把所有时段的控制变量放到一个模型里就变成了多时段动态无功优化。SOCP的求解时间依然可控对33节点来说几十个时段也就几秒钟。做完多时段之后还能进一步把电容器的日动作次数限制、储能充放电计划加进来这样优化结果就更贴近实际运行调度了。我的建议是先在单时段上把框架跑通再逐步加时段别一上来就奔着完整调度模型去调试会很难受。7.3 从IEEE33节点到工程落地阶段性地讲IEEE33节点算例的价值在于帮你建立“从物理问题到数学建模再到算法求解最后回到工程验证”的完整链路。真到了配电网现场节点数可能变成几百个负荷类型五花八门数据还不一定齐全。但核心的建模和求解框架是通用的拿到网络数据标幺化建立DistFlow潮流约束处理成SOCP形式求解校验。这一套动作几乎是固化的。我在实际项目里做过的体会是真正拉开差距的往往不是算法本身而是对数据的处理和结果的工程化解释。SOCP只是帮你在合理时间内找到一组好的无功调节方案但“怎么把这组方案转化成运维人员能理解、能执行的指令”才是价值真正落地的地方。这也是为什么我在前面的实现里反复强调基线校验、松弛检验和离散化验证——它们都是为了让优化结果经得起现场数据的拷问。最后分享一个我的操作习惯每跑完一轮优化一定要把补偿量、网损、各节点电压画成图和基准潮流叠在一起看。曲线的形状比任何单一指标都更能告诉你模型写没写对。如果补偿后电压曲线平滑抬升、末端改善最明显网损也降下来了整个结果就显得很健康如果出现哪个节点电压异常跳高多半是补偿节点选得太密或者数据索引出了问题。养成“看图挑错”的习惯之后你排查模型问题的速度会快很多。 基于二阶锥规划的配电网无功优化：IEEE33节点Matlab实现
1. 从电压越限说起配电网无功优化要解决的真问题1.1 一个整天“顾头不顾腚”的台区图景做配电系统规划的人基本都见过这种场景白天光伏大发配变台区电压往上窜逆变器一台接一台因为过压限功率到了傍晚负荷爬上来线路末端电压又往下掉灯都跟着抖。明明是同一套设备整条馈线却总是顾头不顾腚。这种时候最直接的调节手段就是无功补偿——在关键节点投入电容、调无功把电压拉回安全区间顺带把网损降下来。这也是“配电网无功优化”这个课题为什么反复在IEEE33节点系统上被验证的原因。IEEE33节点系统不算大但有32条支路、33个节点还有5条联络开关馈线长、负荷集中电压问题非常典型。把无功优化算法放到这套数据上跑结果能直观看出优化算法到底帮了多少忙。这篇文章我准备把自己用Matlab实现一套基于二阶锥规划SOCP的配电网无功优化算法的思路完整过一遍从模型的数学推演到代码落地再把跑完结果怎么看、怎么校验说清楚。如果你想做毕业设计、或者要在配电网项目里快速验证某个无功调节方案这篇应该能省你不少试错时间。1.2 无功优化的本质用更小的代价把电压和网损同时按住先理清概念。配电网的无功优化通俗讲就是在保证节点电压不越限的前提下调节电网里的无功源电容器组、静止无功补偿器、分布式逆变器剩余容量等让整个系统的运行状态往更经济的方向走。经济性最直接的指标就是网损——电力从根节点送到末端负荷一路在线路电阻上消耗的有功功率。无功本身不消耗有功但无功电流在线路上流动时会产生额外损耗因为线路损耗等于电流平方乘以电阻而电流是有功分量和无功分量的矢量和。举个例子容易理解一条10千伏馈线上负荷全是居民小区功率因数不高线路里除了有功电流还窝着一大堆无功电流。这时候在负荷附近投入一组电容器电容器发出的无功就地抵消负荷的无功需求从变电站到补偿点这段线路上的无功电流就明显变小线损跟着降。同时少了无功电流在线路电抗上造成的电压降落末端电压也被抬起来了。所以无功优化的核心收益是两个电压质量和运行经济性。既然收益这么明显那问题就是怎么确定“在哪投、投多少、什么时候投”。传统做法是人工经验或者按功率因数简单判断但配电网节点多、负载分布不均匀人工调很难兼顾全局。需要一套自动化的优化算法这就是我们要做的事情。1.3 为什么拿IEEE33节点当“试金石”IEEE33节点系统在文献里出现频率极高几乎成了配电网优化算法的标配算例。这个系统额定电压12.66kV总负荷大约3715kW加2300kvar最大特点是单电源辐射状结构馈线很长末端电压天生偏低。在无补偿条件下末端节点电压会掉到0.91左右网损能到200kW以上非常适合用来验证电压改善和降损效果。对做算法的人来说这套系统还有几个好处数据公开、规模适中、计算速度快跑SOCP基本毫秒级出结果。更关键的是它的拓扑结构清晰潮流结果好验证你无论用前推回代法还是商用软件都能很容易算出一个基准潮流再拿优化结果去对比。所以我后面所有实现细节都会围绕IEEE33节点展开但方法本身完全可以直接迁移到其他配电网拓扑上。2. 为什么选二阶锥规划传统方法到底卡在哪2.1 无功优化的数学模型长什么样先把无功优化写成数学问题。目标函数最常见的选择是最小化全网有功网损minimize Ploss ∑ (P_ij² Q_ij²) / V_i² × R_ij这里的P_ij、Q_ij是支路ij上流过的有功和无功功率V_i是节点i的电压幅值R_ij是支路电阻。这个表达式看着直观但它是个非凸的、带有分式项的非线性函数直接求解很麻烦。约束条件包括三块潮流方程约束有功、无功必须满足功率平衡、运行安全约束节点电压一般在0.95到1.05倍额定值之间、控制变量约束无功补偿容量有上限。如果再把有载调压变压器分接头、分布式电源出力这些加进来约束条件更多模型也更复杂。问题的难点在于潮流方程本身是非线性的。对于给定的负荷潮流方程的解——各节点电压、各支路功率——不能用一个简单的显式公式算出来必须迭代求解。这意味着你在做优化时每改变一个控制变量整个系统的状态都会跟着变化而且这个变化不是简单线性的。2.2 传统非线性求解器做无功优化为什么吃力经典做法是用内点法之类的非线性规划求解器直接解上面这个模型。内点法在IEEE33节点这样的小系统上也能收敛但有几个让人头疼的地方。首先是初值敏感选得不好容易不收敛或者收敛到奇怪的局部最优其次是海森矩阵计算复杂代码写到后面自己都容易被各种偏导搞晕。我曾经用Matlab的fmincon试着求解过小规模无功优化调了各种求解器选项总算稳定跑通了但速度始终不理想而且每次改动约束条件都要重新调初值。后来换到大规模场景或者加入时序曲线fmincon直接趴窝。这还只是小算例真到了几百个节点的配电网非线性规划求解的效率问题会被放大很多。另外一个绕不开的问题是非凸性。无功优化的潮流方程是典型的非凸约束意味着目标函数可能存在多个局部最优点。启发式算法遗传算法、粒子群等经常被用来“绕过”这个问题因为它们不依赖梯度而是通过种群搜索来找解。乍一看挺合适但实际用下来你会发现每次跑结果都不一样能不能找到好的解完全看运气计算时间随节点规模急剧膨胀参数种群大小、交叉概率、变异概率特别影响结果调参本身就是一门玄学。2.3 SOCP提供了一个“凸化”的干净路径二阶锥规划属于凸优化的一种。凸优化有一个非常诱人的性质局部最优解就是全局最优解而且有高效、稳定的求解算法不像启发式那样需要碰运气。为什么叫二阶锥你可以把二阶锥想象成一个三维空间里倒扣的沙堆数学上它是这样一类约束某个向量的欧几里得范数 ≤ 一个线性表达式。这种形状是凸的两个凸集合的交集加上线性目标就能组成一个可以高效求解的凸优化问题。而在配电网领域通过一些变量替换和松弛技巧能把非凸的潮流方程改写成二阶锥约束从而把原本很难的无功优化变成几分钟就能建模、毫秒级能求解的凸优化问题。当然把非凸问题转成凸问题不是免费的它有一个“松弛”过程严格讲是扩大了原问题的可行域。但这个风险是可以控制和后验校验的后面第4章我会把推导过程完整讲一遍你就明白为什么大多数配电网场景下这种松弛是安全的。2.4 解决什么样规模的问题最合适SOCP这种凸优化方法特别适合“约束多、变量多但结构规整”的问题。配电网的节点功率平衡、电压降落、支路容量限制写成矩阵后非常整齐天然适合凸优化框架。IEEE33节点这种规模的算例变量只有一百多个约束也就两三百条对YalmipMosek的组合简直是杀鸡用牛刀毫秒级出结果。就算扩展到几百个节点的馈线、或者加上多时段时序曲线SOCP也依然能用非常短的时间完成求解。所以我最终选择SOCP路线核心原因是它把“能不能收敛”的问题变成了“建模是否严谨”的问题后者比前者可控得多。3. IEEE33节点算例从原始数据到可求解模型3.1 系统基本参数与拓扑结构IEEE33节点系统的基本信息如下基准电压12.66kV基准功率一般取10MVA网络包含1个根节点变电站出口节点0和32个负荷节点正常运行时5条联络开关全部断开系统呈纯辐射状。支路参数用的是有名值单位是欧姆电阻和电抗都给出了。负荷一般也以有功功率kW和无功功率kvar给出。拿到这套数据之后第一件事不是急着写代码而是把它读进Matlab并且检查一遍节点编号和支路连接关系。3.2 支路数据和联络开关的处理细节这里有个新手很容易踩的坑IEEE33节点系统原始数据里包含37条支路其中32条是常规运行支路另外5条是联络开关。在做常规潮流和无功优化时一定要把联络开关排除掉否则网络就成了环网DistFlow模型的前提就不成立了。我在实现里的做法是准备两份数据一份是完整系统数据用于拓扑查看和画图另一份是只含32条闭合支路的“运行支路”数据专门用于潮流和优化计算。这样后续做网络重构类扩展时再切回完整数据也方便。3.3 标幺值不可跳过的一步优化的所有变量最后都要落到标幺值体系里。IEEE33节点系统名义电压12.66kV基准功率取10MVA阻抗基准就是Zbase 12.66² / 10 ≈ 16.031 Ω线路阻抗用有名值除以Zbase得到标幺值负荷用有名值除以10MVA得到标幺值。这里特别提醒一下如果你偷懒直接用有名值建模Yalmip求解器经常会出现数值病态问题因为线路电阻0.4Ω和电压12500V差了四个数量级凸优化求解器对这种尺度不一致非常敏感。标幺化之后电压在1附近功率在0.几量级求解器跑起来又快又稳基本上不会有数值问题。3.4 无功补偿候选节点的选择在IEEE33节点上做无功优化另一个要提前定好的事是在哪些节点允许装无功补偿装置。理论上每个负荷节点都可以装但工程上不可能。我一般会选择几个“对电压支撑最有利”的位置比如线路末端附近节点18附近、几条分支的末端节点17、节点32附近以及负荷集中的节点节点7附近。我的候选节点集合就按这个思路设置比如取节点7、节点18、节点33对应程序索引的第8、19、34位——注意原始标号的偏差这三处。每个节点设置无功补偿容量上限比如单点最大补偿0.6Mvar在10MVA基准下就是0.06p.u.。这个上限定多少需要先跑一遍无补偿潮流看看电压最低点和网损水平再决定一般保证总补偿容量在系统总无功负荷的20%到30%就能看到明显效果。4. 核心建模推演DistFlow如何一步步变成二阶锥4.1 从DistFlow方程说起配电网辐射状结构的潮流通常用DistFlow方程来描述。对一条支路ij设i侧靠近根节点j侧远离根节点P_ij和Q_ij是支路首端流过的有功和无功V_i是首端电压幅值r_ij和x_ij是支路电阻和电抗l_ij代表电流幅值的平方那么V_j² V_i² - 2(r_ij P_ij x_ij Q_ij) (r_ij² x_ij²) l_ijl_ij (P_ij² Q_ij²) / V_i²同时每个节点的注入功率要满足平衡关系从父支路流入的功率减去线路上消耗的功率再减去流向各子支路的功率必须等于该节点的净负荷。第一个公式反映了电压沿线路的降落第二个公式定义了电流与功率、电压之间的关系。难处理的就是第二个公式它带有二次项和分式项让整个潮流约束变成非凸的。4.2 两次变量替换把非线性项藏进去建立二阶锥模型的核心思路是引入新的变量把那些不好处理的非线性项替换掉。第一步用U_i表示V_i²第二步用L_ij表示l_ij。这样一来DistFlow的电压公式变成U_j U_i - 2(r_ij P_ij x_ij Q_ij) (r_ij² x_ij²) L_ij这个方程是线性的很好处理。但原来那个电流定义式呢L_ij (P_ij² Q_ij²) / U_i等价于L_ij × U_i P_ij² Q_ij²。这个约束仍然是非凸的。怎么处理呢这里有一个非常关键的洞察把等式换成不等式。也就是说不再严格要求L_ij × U_i恰好等于P_ij² Q_ij²而是允许L_ij × U_i ≥ P_ij² Q_ij²。这个不等式经过代数变换可以写成‖[2P_ij; 2Q_ij; L_ij - U_i]‖₂ ≤ L_ij U_i这就是一个标准的二阶锥约束。左边是向量的2-范数右边是线性表达式。整个约束描述了一个凸区域完美解决了非凸问题。4.3 这个“松弛”会不会导致错误答案把等式放宽成不等式从数学上看是放大了可行域原来可行的解依然可行但同时也引入了一些原来不可行的解。如果优化器最后找到的最优解落在“不等式严格成立”的区域还没达到等号那这个解虽然满足SOCP模型却不一定满足原始潮流方程严格说是不可行的。好消息是对辐射状配电网和以网损最小化为目标的模型有比较成熟的理论结果在相当宽松的条件下SOCP松弛是“紧的”也就是最优解会落在等号边界上恢复到原始潮流方程。这背后的直觉是网损最小化会尽量“压紧”支路电流和电压之间的关系让不等式没有“松”的空间。而且配电网的负荷是有界的、电压波动范围有限极少出现松弛不紧的情况。但理论是理论工程上不能盲信。所以在代码实现后一定要做一个后验检查逐个支路把优化得到的U_i、P_ij、Q_ij、L_ij代回原等式看看L_ij和(P_ij²Q_ij²)/U_i的残差有多大。如果残差在1e-5以下就可以放心认为松弛是精确的。我在后面第6章会专门讲这个检查步骤。4.4 目标函数选哪个网损最小 vs 电压偏差最小二阶锥框架下目标函数非常灵活。最简单的选择就是最小化网损minimize ∑ L_ij × R_ij因为目标函数和变量L是线性关系网损最小化的目标在SOCP里非常好处理。还有一种常见选择是最小化节点电压偏差minimize ∑ (U_i - 1.0)²这个目标在SOCP里也能写但平方项需要额外处理而且通常电压偏差和网损无法同时最优会有取舍。我的经验是先把网损最小作为主要目标因为网损与无功流动直接相关优化无功补偿对网损的改善非常直观而且结果便于用潮流计算验证。等基本流程跑通了再加入电压偏差的罚项比如目标变成网损加上一个很小的电压偏移惩罚系数兼顾电压质量。4.5 无功补偿设备和网络约束怎么整合进模型无功补偿装置在模型里很简单就是给候选节点增加一个可调的无功注入变量Qc。如果节点j安装了补偿装置那么在节点功率平衡方程的无功部分右边要加上Qc_j同时给Qc_j加上下限约束。分布式逆变器如果具备无功调节能力可以类似地建模只不过还要考虑逆变器有功出力和容量限制。完整的SOCP约束集合包括根节点电压固定为1.0 p.u.非根节点的电压平方U_i在0.95²到1.05²之间每个节点的有功、无功功率平衡方程每条支路的电压降落方程每条支路的二阶锥潮流不等式补偿量Qc的上限和下限。把这些约束全部丢给求解器再配上目标函数整个模型就完整了。5. MatlabYalmip代码实现从零到出结果的关键细节5.1 整体流程设计我用Matlab的Yalmip工具箱来做建模用Mosek作为底层求解器。Yalmip的好处是它把凸优化建模从“手写求解器”里解放出来了你只需要声明变量、写约束、写目标函数剩下的交给它处理。整体流程是载入IEEE33节点数据转换为标幺值生成Yalmip优化变量电压平方、支路功率、支路电流平方、无功补偿量逐条构建约束设置求解器选项调用optimize求解从求解结果中恢复电压、支路功率和补偿量运行一个前推回代潮流对优化结果做校验和误差分析。5.2 核心代码框架下面这段代码是我实现的核心框架关键变量和约束构成了一个可跑通的SOP方案。你可以直接在此基础上扩展。%% IEEE33节点二阶锥无功优化——核心代码框架 clear; clc; %% 1. 数据读入示意结构数据需自行准备 % bus: [节点号, 有功负荷kW, 无功负荷kvar] % branch: [首端节点, 末端节点, 电阻ohm, 电抗ohm] % 注意此处只取32条闭合支路联络开关不参与建模 Vbase 12.66; % 基准电压 kV Sbase 10; % 基准功率 MVA Zbase Vbase^2 / Sbase; % 基准阻抗 nb 33; % 节点数 nl 32; % 运行支路数 Rpu branch(:,3) / Zbase; Xpu branch(:,4) / Zbase; PLpu bus(:,2) / Sbase / 1000; % kW - MW - p.u. QLpu bus(:,3) / Sbase / 1000; %% 2. 定义优化变量 U sdpvar(nb, 1); % 节点电压平方 P sdpvar(nl, 1); % 支路有功 Q sdpvar(nl, 1); % 支路无功 L sdpvar(nl, 1); % 支路电流平方 qc_nodes [8, 19, 34]; % 候选补偿节点对应的程序索引 Qc sdpvar(length(qc_nodes), 1); % 无功补偿量 %% 3. 目标函数与约束初始化 objective sum(L .* Rpu); % 网损最小标幺值 constraints {}; % 根节点电压约束 constraints{end1} U(1) 1.0; % 电压上下限 constraints{end1} 0.95^2 U 1.05^2; % 补偿容量约束 constraints{end1} 0 Qc 0.06; % 单点最大0.6Mvar %% 4. 节点功率平衡约束 % 对每个节点父支路流入功率 本地注入 - 子支路流出 0 % 使用矩阵实现但这里为了可读性用循环示例 for k 1:nl i branch(k,1) 1; % 程序索引转换为1起始 j branch(k,2) 1; % 有功平衡支路k的有功流入父支路末端节点j从首端节点i流出 % 更严谨的写法是用稀疏矩阵累加这里示意 end %% 5. 支路电压方程与二阶锥约束 constraints{end1} ... % 逐个支路写 % U(j) U(i) - 2*(Rpu(k)*P(k) Xpu(k)*Q(k)) (Rpu(k)^2Xpu(k)^2)*L(k); %% 6. 求解 ops sdpsettings(solver, mosek, verbose, 1); optimize(constraints, objective, ops); %% 7. 结果恢复 U_opt value(U); P_opt value(P); Q_opt value(Q); L_opt value(L); Qc_opt value(Qc) * Sbase * 1000; % 转换为kvar5.3 为什么用Yalmip而不是手写求解器有人可能会问直接用Mosek的建模语言写SOCP不就行了确实可以但Yalmip在配电网优化里有两个明显优势。一是语法很贴近数学表达式修模型方便比如加一个电压约束就是U 1.05^2这么直白二是能灵活切换求解器Mosek没装好就先换Sedumi代码改动几乎只有一行。对于算法验证和小型项目Yalmip是效率最高的选择。5.4 求解器选型Mosek、Sedumi、Gurobi怎么选我在IEEE33节点这个算例里主要推荐Mosek。它对SOCP的支持很成熟数值稳定速度非常快而且能自动检测凸性。Sedumi是免费开源的老牌求解器代码量小、可靠也能处理SOCP但速度比Mosek慢不少在大规模问题上尤其明显。Gurobi现在也支持SOCP但在旋锥约束上支持得不如Mosek顺手。实际测试下来33节点的SOCP模型Mosek通常几十毫秒以内就能完成求解Sedumi可能要一两秒。如果你只是学习Sedumi完全够用如果后面要扩展到多时段、几百节点建议直接上Mosek省下来的时间非常可观。5.5 新手最容易翻车的三个位置第一个是支路端点索引。原始IEEE33节点数据从0开始编号但Matlab的下标从1开始两者差一位转换错了整个网络拓扑就乱了。我的习惯是一进代码就统一加1后面所有求解器变量、矩阵操作都用转换后的索引。第二个是标幺值单位换算。负荷数据单位是kW和kvar转标幺值的时候要先除以1000转成MW再除以Sbase10MVA。这个地方马虎一下结果可能差出一千倍。我一般会在读入数据后打印一次各节点总负荷和文献值对比确认是对的再继续往下写。第三个是补偿容量的续接。优化求解器给的Qc是连续变量、标幺值工程上电容器是离散的一般成组投切。所以求完解要按实际的单组容量比如0.1Mvar做圆整再把圆整后的补偿量重新带回潮流验证电压和网损是否依然满足要求。6. 结果怎么看电压、网损与松弛间隙6.1 先跑一版基准潮流一切优化才有参照拿到优化结果之前一定先把无补偿时的系统状态算出来。我在IEEE33节点上跑出来的经典基线数据大致是全网有功损耗在202kW左右最低电压出现在18号节点附近大约0.91p.u.。这个数值和大家公开发表的结果基本一致你可以把自己算出来的基线和它对一下如果对不上先回去查数据而不是急着调优化模型。无补偿时电压已经跌破0.95的下限这正好说明这个系统确实存在无功支撑不足的问题。接下来优化器给出的补偿方案才显得有意义。6.2 优化前后的数据对比加入SOCP优化后我用YalmipMosek跑出来的典型结果如下补偿器总投入大约在1Mvar到1.5Mvar之间全网有功损耗从约202kW降到约108kW到130kW区间降幅显著最小节点电压从0.91p.u.抬升到0.95p.u.以上全网电压全部进入安全区间求解耗时通常不到0.1秒。这个结果的逻辑很顺优化器在电压最薄弱的节点附近投入了较多无功补偿抬升了末端电压同时减少了从变电站到末端的无功流动网损自然降下来了。6.3 松弛间隙检查结果真的可信吗这是整个流程里最不该省的一步。SOCP把等式放宽成了不等式你必须在求解结束后检查每个支路看松弛后的最优解有没有回到原始等式上。具体做法是对每条支路计算不等式两边L_ij × U_i和P_ij² Q_ij²之间的相对误差。如果所有支路的残差都在1e-5量级甚至更小说明优化器自动把不等式“推紧”到了边界SOCP松弛是精确的解可以直接用。如果某些支路残差很大说明这个算例下SOCP松弛不完全精确你需要回到模型检查是不是电压范围设得太宽、负荷太重、或者网络结构不适合用DistFlow建模。有几个知名文献里专门讨论过哪些极端场景会让松弛不紧后面第7章我会提到。6.4 从连续解到工程可实施的完整闭环SOCP给出的补偿容量是连续值比如节点18的Qc解出来是0.43Mvar。但真实电容器组往往是离散的一组0.05或0.1Mvar。这时候要做一次圆整把0.43Mvar取到0.4Mvar或0.45Mvar然后把圆整后的值重新放回潮流方程里算一次前推回代验证电压和网损是否仍满足要求。我习惯把这一步做成流程里的固定动作因为优化模型和真实设备之间永远有这一步“兑现”过程。省略了它你的结果最多只能算理论下限不能称为可实施方案。圆整后如果有个别节点电压又跌回0.95以下那就把补偿投资限制稍微放宽一点重新优化一轮直到离散方案和连续方案的结果偏差在可接受范围内。7. 边界条件与扩展思路7.1 什么时候SOCP松弛会不靠谱虽然SOCP配电网无功优化在大多数辐射状场景下表现出色但有些情况你必须警惕。最典型的是三相不平衡严重的低压馈线单相模型不再成立DistFlow天然失真还有合环运行状态网络出现环网结构前面讲的精确性条件容易失效再比如电压跌落极深或包含有载调压变压器分接头时模型的凸性也会受影响。这些情况下你不能直接把SOCP结果当真至少要把算出来的控制方案拿到完整的三相潮流里做校验。如果校验不合格再考虑是否要做混合整数规划加回那些离散变量和环网约束。7.2 从静态优化扩展到多时段动态无功优化规模上的扩展其实非常自然。现在的模型是“一个时刻点”的静态优化如果拿一条24小时或96点负荷曲线把每个时段的节点负荷作为输入把所有时段的控制变量放到一个模型里就变成了多时段动态无功优化。SOCP的求解时间依然可控对33节点来说几十个时段也就几秒钟。做完多时段之后还能进一步把电容器的日动作次数限制、储能充放电计划加进来这样优化结果就更贴近实际运行调度了。我的建议是先在单时段上把框架跑通再逐步加时段别一上来就奔着完整调度模型去调试会很难受。7.3 从IEEE33节点到工程落地阶段性地讲IEEE33节点算例的价值在于帮你建立“从物理问题到数学建模再到算法求解最后回到工程验证”的完整链路。真到了配电网现场节点数可能变成几百个负荷类型五花八门数据还不一定齐全。但核心的建模和求解框架是通用的拿到网络数据标幺化建立DistFlow潮流约束处理成SOCP形式求解校验。这一套动作几乎是固化的。我在实际项目里做过的体会是真正拉开差距的往往不是算法本身而是对数据的处理和结果的工程化解释。SOCP只是帮你在合理时间内找到一组好的无功调节方案但“怎么把这组方案转化成运维人员能理解、能执行的指令”才是价值真正落地的地方。这也是为什么我在前面的实现里反复强调基线校验、松弛检验和离散化验证——它们都是为了让优化结果经得起现场数据的拷问。最后分享一个我的操作习惯每跑完一轮优化一定要把补偿量、网损、各节点电压画成图和基准潮流叠在一起看。曲线的形状比任何单一指标都更能告诉你模型写没写对。如果补偿后电压曲线平滑抬升、末端改善最明显网损也降下来了整个结果就显得很健康如果出现哪个节点电压异常跳高多半是补偿节点选得太密或者数据索引出了问题。养成“看图挑错”的习惯之后你排查模型问题的速度会快很多。