写配电自动化项目的调试笔记，我一般都从“凭什么相信这个优化结果”开始讲。因为很多朋友拿到二阶锥规划（SOCP）这个词，第一反应是“高大上”，但落到配电网无功优化里，真正要解决的是这么一个问题：辐射状网的非线性潮流方程，怎么放进一个能保证全局最优的凸优化框架里。这半年我在33节点标准测试系统上完整实现了一遍，从DistFlow建模、锥松弛、YALMIP代码、求解器选型到工程化扩展，踩了不少坑，也沉淀下一套可以直接复用的流程。这篇就按我实际做项目的顺序，把关键环节拆开讲清楚，适合配电方向的研究生、刚接触凸优化的工程师，也适合想把算法落地到分布式电源场景的朋友。1. 无功优化到底在优化什么：从目标函数到潮流约束的数学骨架1.1 网损最小不等于“省电费”无功优化最常见的目标函数是网损最小，但只盯着网损会吃大亏。真实配电网里，分布式光伏中午大发时末端电压可能冲到1.07 pu，夜间轻载时长线路末端电压可能掉到0.93 pu，这些电压越限问题往往比网损几个百分点更致命。所以我在做模型时，通常把目标函数写成网损项加电压偏差惩罚：$$ \min \quad \sum_{ij \in E} r_{ij} l_{ij} + \lambda \sum_{i \in N} (V_i - 1)^2 $$其中 $r_{ij}$ 是支路电阻，$l_{ij}$ 是支路电流平方，第一项就是各支路有功损耗之和。$\lambda$ 是电压偏差惩罚系数，取多少很敏感：太小，优化结果可能把电压推向边界；太大，又会稀释网损优化的权重。我的经验是先不加惩罚项，只把节点电压作为硬约束，比如 $0.95^2 \le V_i \le 1.05^2$，跑一遍看有没有可行解；如果边界场景下无解或最优解贴着边界，再引入电压软约束，并给 $\lambda$ 一个递增的试探序列。电网运行考核里还有一层容易被忽略的诉求：无功倒送。现在很多配电线路没有无功表计，但关口表能看到功率因数异常偏低，这本质上是变电站向线路输送了大量无功。做目标函数时可以把变电站关口无功也纳入惩罚，或者用约束硬性限制注入无功上限。总之，纯学术demo可以只写网损最小，工程落地版的模型一定要加电压和关口无功的“安全护栏”。1.2 DistFlow潮流约束与失去凸性的那一刻配电网是辐射状拓扑，用支路潮流方程（DistFlow）建模非常自然。对每条支路 $i \to j$，要满足四组物理约束：$$ P_{ij} - \sum_{k \in \text{child}(j)} P_{jk} - r_{ij} l_{ij} = P_{j}^{\text{load}} - P_{j}^{\text{gen}} $$$$ Q_{ij} - \sum_{k \in \text{child}(j)} Q_{jk} - x_{ij} l_{ij} = Q_{j}^{\text{load}} - Q_{j}^{\text{gen}} $$$$ V_j = V_i - 2(r_{ij} P_{ij} + x_{ij} Q_{ij}) + (r_{ij}^2 + x_{ij}^2) l_{ij} $$$$ l_{ij} V_i = P_{ij}^2 + Q_{ij}^2 $$第一行是节点有功平衡：父支路送过来的有功，扣掉线上损耗，分配给各子支路以及本地负荷。第二行是无功平衡，含义一样，只是把电阻损耗替换成电抗压降带来的无功损耗。第三行是电压降落关系，第四行是电流平方乘以电压平方等于功率平方。问题就出在第四行的等式上。它把支路电流、首端电压和传输功率用二次形式强行捆在一起，这个等式在当前变量空间里是非凸的。直观理解的话，优化问题的可行域变成了一条“拧巴”的曲面，直接丢给内点法求解，大概率收敛到局部最优，而且你根本不知道当前解离全局最优有多远。这也是为什么很多启发式算法（粒子群、遗传算法）在配网无功优化里虽然能“找到解”，却很难给出严格最优性证明。1.3 决策变量、补偿设备边界与等式梳理我习惯把决策变量分成四组：支路有功 $P$、支路无功 $Q$、支路电流平方 $l$、节点电压平方 $V$，外加每个节点上的无功补偿量 $Qc$。注意这里电压和电流都处理成“平方”变量，不是一时兴起，而是为二阶锥松弛铺路——如果直接使用原始电压 $U$ 和原始电流 $I$，约束还是非凸的，松弛起来非常别扭。补偿设备的边界必须认真对待。连续可调的SVG可以直接加上下限，比如 $-0.3 \le Qc_i \le 0.3$（标幺值）；分组投切电容器则要引入整数变量，后续会变成混合整数二阶锥规划（MISOCP）。实际项目里还有一个隐藏坑：台账上的设备容量不等于实际可用容量，比如某节点电容柜一个支路坏了，可用无功上限其实只有额定值的60%。模型里宁可把边界设保守一点，也不要盲目相信台账，否则优化结果会在执行时被打脸。2. 二阶锥松弛的关键一步：把非凸潮流等式换成凸锥约束2.1 旋转锥与标准二阶锥的换算前文提到的非凸等式：$$ l_{ij} V_i = P_{ij}^2 + Q_{ij}^2 $$这一步的思路很直接：把等号放宽成不等号：$$ l_{ij} V_i \ge P_{ij}^2 + Q_{ij}^2 $$可行域从一条非凸曲线变成一个“碗状”凸锥。这看起来像作弊，但数学上是有边界条件的：当网络是辐射状、目标函数在变量上单调递增、且最优时该不等式取到等号，那么松弛前后的最优解完全一致。这个性质叫“紧性”（exactness）。为了喂给二阶锥求解器，还要把上面这个旋转锥换算成标准二阶锥形式。两边同时乘4，再补一个 $V_i - l_{ij}$ 的平方项：$$ 4 P_{ij}^2 + 4 Q_{ij}^2 + (V_i - l_{ij})^2 \le (V_i + l_{ij})^2 $$写成向量形式就是：$$ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ V_i - l_{ij} \end{bmatrix} \right|2 \le V_i + l{ij} $$在YALMIP里，这个形式直接用cone([2*P(k); 2*Q(k); V(i)-L(k)], V(i)+L(k))一句话写进去，非常清爽。为什么这一步能保证全局最优？因为凸优化问题的任何局部最优解都是全局最优，这是凸分析的基本结论。配电网支路潮流约束一旦变成二阶锥，整个可行域变凸，目标函数又是线性的，求解器自然可以给你一个带全局最优性保证的解。这在工程评估里特别重要：领导问“你这套方案比起粒子群好在哪”的时候，我可以直接回答“粒子群给你一个看起来不错的解，SOCP给你一个证明过的最优解”。2.2 目标函数保持线性，模型整体变成SOCP松弛之后回头看整个模型，会发现一个幸运之处：目标函数里的网损 $\sum r_{ij} l_{ij}$ 关于新变量 $l$ 是线性的；电压上下限约束 $0.95^2 \le V_i \le 1.05^2$ 是线性的；节点功率平衡方程全部是线性等式组合；调度变量边界也是线性不等式。唯一非线性的地方就是那个二阶锥约束。所以整个问题最终是一个标准SOCP：线性目标 + 线性约束 + 二阶锥约束。这种结构下的求解理论和算法都非常成熟，YALMIP + MOSEK/Gurobi/SeDuMi 都能处理。常见的困惑是“为什么不用SDP松弛”？半定规划（SDP）在配电网重构等领域确实也用，但SOCP的约束是逐条支路独立锥形，规模更小、求解更快；对辐射状网的支路潮流模型来说，SOCP已经足够，不需要更广义的SDP，尤其在节点数上千的场景里，SOCP的数值稳定性优势很明显。2.3 Matlab+YALMIP代码骨架怎么搭搭建代码骨架时，我强烈建议把“拓扑组织”和“建模”分开。否则约束里一会儿用下标、一会儿用矩阵，很容易出错。下面是一段核心骨架：%% 定义变量 P = sdpvar(nb,1); % 支路有功 Q = sdpvar(nb,1); % 支路无功 L = sdpvar(nb,1); % 支路电流平方 V = sdpvar(nnode,1); % 节点电压平方 Qc = sdpvar(nnode,1); % 节点无功补偿注入 %% 拓扑准备：先构建 parent、children % children{j,1} = 节点 j 的所有子支路编号 % parent_br(j) = 节点 j 的父支路编号，根节点为0 % 该结构需要从 from/to 数组通过 BFS 或拓扑排序生成 Con = []; Con = [Con, V(1) == 1.0]; % 根节点电压 1.0 pu Con = [Con, 0.95^2 <= V <= 1.05^2]; % 电压上下限 Con = [Con, Qc_min <= Qc <= Qc_max]; % 补偿设备边界 %% 节点功率平衡（按节点遍历，最直观） for j = 1:nnode child_br = children{j}; if j == 1 || parent_br(j) == 0 inflow_p = 0; inflow_q = 0; else fk = parent_br(j); inflow_p = P(fk) - mpc.R(fk) * L(fk); inflow_q = Q(fk) - mpc.X(fk) * L(fk); end outflow_p = sum(P(child_br)); outflow_q = sum(Q(child_br)); Con = [Con, inflow_p - outflow_p == mpc.Pgen(j) - mpc.Pd(j)]; Con = [Con, inflow_q - outflow_q == mpc.Qgen(j) - mpc.Qd(j) + Qc(j)]; end %% 电压降落方程 for j = 2:nnode fk = parent_br(j); i = mpc.from(fk); Con = [Con, V(j) == V(i) - 2*(mpc.R(fk)*P(fk) + mpc.X(fk)*Q(fk)) ... + (mpc.R(fk)^2 + mpc.X(fk)^2) * L(fk)]; end %% 二阶锥约束 for k = 1:nb i = mpc.from(k); Con = [Con, cone([2*P(k); 2*Q(k); V(i)-L(k)], V(i)+L(k))]; end %% 目标函数与求解 obj = sum(mpc.R .* L); ops = sdpsettings('solver','mosek','verbose',1); sol = optimize(Con, obj, ops);这段代码可以直接铺开来跑小算例。里面有几个细节：根节点用V(1) == 1.0，因为我们变量本身就是电压平方，所以根节点电压平方等于1.0；功率平衡里，inflow_p是父支路送到本节点的净有功，outflow_p是本节点送往所有子支路的总有功，右边是发电减负荷；补偿无功Qc放在无功平衡等式的右边，正值代表向节点注入无功。3. IEEE33节点算例搭建与Matlab代码实现细节3.1 数据准备：标幺化、拓扑脚本和从哪里修改IEEE33节点测试系统是配电网算法验证的“标准题”：33个节点、32条支路、基准电压12.66kV、总负荷3715kW+2300kvar、根节点连接变电站。这里最容易被新手坑的就是标幺化。如果直接把有名值塞进去，电压、电流、功率三个物理量量级差好几档，求解器数值稳定性和迭代速度都会受影响。我用的是 10MVA 基准功率，12.66kV 基准电压，总负荷折算下来大约 0.3715 + j0.23 pu，支路电阻电抗都在 $10^{-4} \sim 10^{-2}$ 这个量级，数值上很舒服。数据脚本我建议单独放一个case33.m，返回一个结构体：mpc.nnode = 33; mpc.nb = 32; mpc.from = [...]; % 32x1 首端节点编号 mpc.to = [...]; % 32x1 末端节点编号 mpc.R = [...]; % 32x1 支路电阻(pu) mpc.X = [...]; % 32x1 支路电抗(pu) mpc.Pd = zeros(33,1); % 节点有功负荷(pu) mpc.Qd = zeros(33,1); % 节点无功负荷(pu) mpc.Pgen = zeros(33,1); % 节点有功电源注入(pu)，根节点1通常是变电站 mpc.Qgen = zeros(33,1);支路数据顺序很重要，建议直接按尾端编号递增排列，这样用拓扑排序生成parent、children结构时不容易乱。负荷数据要放在末端节点，例如支路1-2上面的负荷放在节点2，而不是支路上。3.2 主程序流程：建模-求解-取结果主程序我把流程拆成固定的七步，日常调参基本就是在这个骨架上改参数、改约束：步骤做什么关键变量/函数1读入网络数据mpc = case33;2生成拓扑关系children,parent_br3定义优化变量sdpvar4写约束根电压、电压限值、设备边界、平衡方程、电压方程、锥约束5写目标函数obj = sum(mpc.R .* L);6调用求解器optimize(Con, obj, ops);7提取并检查结果value(P),value(L),value(V),value(Qc)第七步特别容易被忽略。很多同学跑完只看一个目标函数就完事，但value(V)的逐节点曲线、value(Qc)的配置方案才是工程上最关心的。我会在后处理里计算三个量：网损现值：loss = sum(mpc.R .* value(L));最低电压：Vmin = min(sqrt(value(V)));松弛间隙：gap = max(value(V(mpc.from(k))) .* value(L) - value(P).^2 - value(Q).^2);其中gap是判断二阶锥松弛“紧不紧”的核心证据，后面专门讲。3.3 求解器选择与结果初勘求解器选择这块，我按不同场景给出建议：求解器类型适用场景注意事项SeDuMi免费SOCP小算例验证、学生练习数值容差相对弱，模型病态时容易不收敛SDPT3免费SOCP/SDP中等算例对稀疏大模型速度一般MOSEK商业SOCP生产级应用、大规模网络数值稳健，能很好地处理大规模SOCPGurobi商业SOCP/MISOCP带整数变量的场景做混合整数二阶锥很顺手在我自己的项目里，仿真阶段首选MOSEK，因为它的容差控制和警告信息更友好；如果模型里面有电容器组投切这种0/1变量，我就切到Gurobi，省得外部循环烷法。跑通之后第一次看结果，别急着画图。先看三件事：sol.problem是不是0，松弛间隙gap是不是小于1e-6，电压曲线上有没有突变点。如果这三关过了，再往下做场景对比和参数敏感性分析。4. 松弛紧性、求解器与收敛性：调优过程中的真实坑4.1 三组仿真场景的对比结果我在33节点系统上做了三组场景，分别对应工程里最常见的三种状态：场景配置网损最低电压松弛间隙1无任何无功优化约202kW（0.0202pu）约0.903pu无2连续无功补偿优化约150kW>0.95pu<1e-73补偿 + 分布式电源无功联合优化约130kW0.97~1.03<1e-8场景1的网损在202kW左右，这是IEEE33标准算例的经典参考值。场景2我在末端电压薄弱的节点附近增加了可连续调节的无功注入，网损降了约四分之一，最低电压提到0.95以上。场景3进一步把分布式电源的无功能力纳入优化，损失更低，更重要的是白天光伏大发和夜间轻载两组工况下电压都稳住了。这些具体数值依赖于基准值和设备边界，不必迷信，但趋势是稳定的：单纯调无功对网损的改善通常能做到15%~30%，而电压合格率往往从“临界”到“稳定”，收益比网损更可观。4.2 松弛间隙检查：判断“凸化是否安全”SOCP结果出来以后，第一件事永远是看松弛间隙。所谓“紧”，可以简单理解成不等式被“顶满”到了原等式的状态；如果gap超过了1e-4，说明模型在某处被松弛得太厉害，数学上可能是“安全”的凸解，物理上却已经不忠实于真实的潮流方程。这个结果拿去工程汇报是要出问题的。我遇到过的gap变大的原因主要有三个：一是目标函数被改坏了。比如原本应该最小化网损，结果想当然加了一个电压偏差的二范数惩罚，比如1e3*sum((V-1).^2)，这个强二次项在某些边界场景下会干扰松弛紧性。解决办法是尽量用线性惩罚或者干脆不加，电压交给硬约束。二是根节点电压范围设置得过松。如果你把根节点电压也放进优化变量，还给了0.9~1.1的宽范围，那么整个系统都有了“钻空子”的空间，求解器完全可能在某个支路上用偏大的电压偏差换取更低的损耗，但实际电网根节点电压是受变电站控制的。稳妥做法是根节点直接V(1)==1.0，把这个自由度锁死。三是设备边界设得离谱。曾有一版代码我把某个补偿点上限写成2pu（对应20Mvar），本来只是测试用的宽松值，结果gap直接到1e-2。原因很直白：补偿上限宽到离谱之后，SOCP在可行域里找到一条“非物理但凸可行”的路，让某条中段支路变相“吸收”了大量虚拟无功。把设备上限改回真实的0.05pu后，gap立刻降到1e-8。这个案例后来被我写进了团队的调试手册：设备边界必须按现场台账的保守值来，越宽越危险。4.3 建模阶段容易踩的四个错误我在带团队做类似算法时，发现有些错误反复出现，值得单独列出来：第一个错误，是在YALMIP里直接写二次非凸约束，比如P(k)^2 + Q(k)^2 <= V(i)*L(k)。YALMIP能接住这个写法，但会把它当成一般非线性约束，交给内点法去硬解，既慢又不保证全局最优。正确做法是一步到位写成cone([2*P(k); 2*Q(k); V(i)-L(k)], V(i)+L(k))，把问题明确标记成二阶锥。第二个错误，是节点功率平衡漏掉支路损耗。很多示例代码会直接写成“流入=流出+负荷”，却忘了父支路送进来的功率在到达本节点时已经扣掉了r*l的损耗。漏掉这一项，优化出来的网损会比真实值低，结果对不上实际。第三个错误，是补偿设备重复配置。一次排查时我发现gap居高不下，翻数据才发现同一个节点在“新增无功节点”和“原有电容器节点”两张表里各自出现了一次，相当于给同一点塞了两组补偿变量。这类问题在真实配网拓扑数据里尤其常见，节点编号一旦合并过就会踩坑。第四个错误，是初值设置欠佳。SOCP本质上迭代求解，给一个物理合理的初值能显著减少迭代次数。我现在习惯在求解前用前推回推潮流先算一组静态工作点，然后用assign(P, value0)把初值塞给变量；即使懒得做潮流，也可以用一段全零或全一初值给求解器一个更稳的起点。4.4 性能与数值稳定性小经验33节点秒级可解，但这个模型扩展到几百上千节点时会遇到性能瓶颈。一个经验是：等节点规模大了以后，别再用一堆for循环逐条加约束，YALMIP对大规模模型的解析时间比求解时间还长。我后来把功率平衡和电压方程全部改写成稀疏矩阵形式，用A*[P;Q;L] == b一次性构建约束，速度提升非常明显。另外，ops里建议加两样东西：sdpsettings('saveyalmipmodel',1)和sdpsettings('debug',1)。前者方便排查“模型是不是非凸”，后者在YALMIP报错时会输出更详细的约束分析，基本能定位到具体是哪一行约束出问题。有一次报错说“quadratic constraint nonconvex”，我用这两个设置立刻追到了漏写锥的那一行。5. 从33节点到工程场景：扩展思路与落地思考5.1 分布式电源接入后的时序无功优化33节点只是实验室里的“验证题”，真正工程现场是分布式光伏、储能、可调负荷同时存在的动态场景。光伏出力一天内波动非常大，如果只做单断面优化，会出现设备频繁动作的问题：上午9点投一组电容器，下午2点又要求退出，一天折腾十几次，现场开关设备根本受不了。处理思路一般有三种：一是多时段联合优化，把24小时分成96个15分钟断面，目标函数里加入设备动作惩罚，让求解器自己权衡“多投一组电容省下的网损”和“动作一次带来的开关损耗”；二是滚动优化，每15分钟滚动求解一次，但把前一时段的结果作为下一时段的初始可行点，减少突变；三是给动作次数加硬约束，比如一天之内同一补偿点的动作次数不超过5次。值得多说一句，逆变器无功能力不是固定的，有功出力高时无功上限会收缩，约束关系是 $P^2 + Q^2 \le S_n^2$。这个约束本身就是一个旋转锥，二阶锥规划正好能直接表达，所以在SOCP框架里做光伏无功优化反而比传统内点法更顺手。5.2 离散设备、三相不平衡与更大规模网络再往工程走，会遇到三个更现实的问题。第一个是离散设备：分组电容器、有载调压变压器分接头都是整数变量。把0/1变量加进SOCP，模型就升级成MISOCP。如果整数变量数量不多，直接用Gurobi可以一次求到最优；如果规模太大，整数变量上百个，业界常用两层算法：外层用启发式枚举离散组合，内层还是SOCP，兼顾速度和全局搜索能力。第二个是三相不平衡。园区型配电网、农村电网很多都存在严重三相不平衡，单相模型会掩盖电压偏相问题。扩展思路是把每一相都作为独立网络建模，支路潮流从一维变量变成三维向量，约束数量翻三倍，但锥结构仍然保留。工程上更折中的做法是先做单相等效优化，再用最终配置去三相潮流里校验一遍不平衡度，这样改动小，实战中也够用。第三个是网络规模。到了几千节点，Matlab + YALMIP的建模开销会很大，此时可以考虑把建好的模型导出成MPS或LP格式，直接交给求解器批量算；或者直接切换到Python生态，用成熟的建模工具重写一版，后续对接配网系统也更方便。5.3 我从这套流程里沉淀下来的实操习惯在多次迭代和排错之后，我现在每次上新网络数据都固定执行一套校验动作，这些习惯帮我避免了很多“看起来能跑、结果不能用”的情况。第一，拓扑校验。每个非根节点有且仅有一个父节点，全图无环。不要以为数据表里写着“配电网”就真是辐射状——很多历史台账里混着联络开关闭合状态的网络，闭合后就是环网，SOCP算出来照样有解，但不是你要的辐射状结果。第二，设备容量校验。进模型前把电容柜投运状态、损坏支路、额定容量重新核对一遍，取保守可用值作为边界。这个在项目验收阶段特别重要，因为优化方案最终是要被执行机构操作的，参数写错，执行机构就乱动。第三，结果报告固定化。每次算完都输出一张“四加一”报告：网损、最低/最高电压、松弛间隙、设备执行表。评审会上，这张表比任何算法解释都更有说服力。第四，也是我最想提醒的：二阶锥规划不是万能钥匙，但它确实是目前辐射状配电网无功优化里“理论保证”和“工程效率”权衡得最好的工具。只要松弛间隙检查到位、设备边界贴合实际，这个框架从研究到落地是走得通的。希望这篇整理能帮你少走几个弯路，尤其是把那个容易被忽略的cone约束一次写对。
写配电自动化项目的调试笔记，我一般都从“凭什么相信这个优化结果”开始讲。因为很多朋友拿到二阶锥规划（SOCP）这个词，第一反应是“高大上”，但落到配电网无功优化里，真正要解决的是这么一个问题：辐射状网的非线性潮流方程，怎么放进一个能保证全局最优的凸优化框架里。这半年我在33节点标准测试系统上完整实现了一遍，从DistFlow建模、锥松弛、YALMIP代码、求解器选型到工程化扩展，踩了不少坑，也沉淀下一套可以直接复用的流程。这篇就按我实际做项目的顺序，把关键环节拆开讲清楚，适合配电方向的研究生、刚接触凸优化的工程师，也适合想把算法落地到分布式电源场景的朋友。1. 无功优化到底在优化什么：从目标函数到潮流约束的数学骨架1.1 网损最小不等于“省电费”无功优化最常见的目标函数是网损最小，但只盯着网损会吃大亏。真实配电网里，分布式光伏中午大发时末端电压可能冲到1.07 pu，夜间轻载时长线路末端电压可能掉到0.93 pu，这些电压越限问题往往比网损几个百分点更致命。所以我在做模型时，通常把目标函数写成网损项加电压偏差惩罚：$$ \min \quad \sum_{ij \in E} r_{ij} l_{ij} + \lambda \sum_{i \in N} (V_i - 1)^2 $$其中 $r_{ij}$ 是支路电阻，$l_{ij}$ 是支路电流平方，第一项就是各支路有功损耗之和。$\lambda$ 是电压偏差惩罚系数，取多少很敏感：太小，优化结果可能把电压推向边界；太大，又会稀释网损优化的权重。我的经验是先不加惩罚项，只把节点电压作为硬约束，比如 $0.95^2 \le V_i \le 1.05^2$，跑一遍看有没有可行解；如果边界场景下无解或最优解贴着边界，再引入电压软约束，并给 $\lambda$ 一个递增的试探序列。电网运行考核里还有一层容易被忽略的诉求：无功倒送。现在很多配电线路没有无功表计，但关口表能看到功率因数异常偏低，这本质上是变电站向线路输送了大量无功。做目标函数时可以把变电站关口无功也纳入惩罚，或者用约束硬性限制注入无功上限。总之，纯学术demo可以只写网损最小，工程落地版的模型一定要加电压和关口无功的“安全护栏”。1.2 DistFlow潮流约束与失去凸性的那一刻配电网是辐射状拓扑，用支路潮流方程（DistFlow）建模非常自然。对每条支路 $i \to j$，要满足四组物理约束：$$ P_{ij} - \sum_{k \in \text{child}(j)} P_{jk} - r_{ij} l_{ij} = P_{j}^{\text{load}} - P_{j}^{\text{gen}} $$$$ Q_{ij} - \sum_{k \in \text{child}(j)} Q_{jk} - x_{ij} l_{ij} = Q_{j}^{\text{load}} - Q_{j}^{\text{gen}} $$$$ V_j = V_i - 2(r_{ij} P_{ij} + x_{ij} Q_{ij}) + (r_{ij}^2 + x_{ij}^2) l_{ij} $$$$ l_{ij} V_i = P_{ij}^2 + Q_{ij}^2 $$第一行是节点有功平衡：父支路送过来的有功，扣掉线上损耗，分配给各子支路以及本地负荷。第二行是无功平衡，含义一样，只是把电阻损耗替换成电抗压降带来的无功损耗。第三行是电压降落关系，第四行是电流平方乘以电压平方等于功率平方。问题就出在第四行的等式上。它把支路电流、首端电压和传输功率用二次形式强行捆在一起，这个等式在当前变量空间里是非凸的。直观理解的话，优化问题的可行域变成了一条“拧巴”的曲面，直接丢给内点法求解，大概率收敛到局部最优，而且你根本不知道当前解离全局最优有多远。这也是为什么很多启发式算法（粒子群、遗传算法）在配网无功优化里虽然能“找到解”，却很难给出严格最优性证明。1.3 决策变量、补偿设备边界与等式梳理我习惯把决策变量分成四组：支路有功 $P$、支路无功 $Q$、支路电流平方 $l$、节点电压平方 $V$，外加每个节点上的无功补偿量 $Qc$。注意这里电压和电流都处理成“平方”变量，不是一时兴起，而是为二阶锥松弛铺路——如果直接使用原始电压 $U$ 和原始电流 $I$，约束还是非凸的，松弛起来非常别扭。补偿设备的边界必须认真对待。连续可调的SVG可以直接加上下限，比如 $-0.3 \le Qc_i \le 0.3$（标幺值）；分组投切电容器则要引入整数变量，后续会变成混合整数二阶锥规划（MISOCP）。实际项目里还有一个隐藏坑：台账上的设备容量不等于实际可用容量，比如某节点电容柜一个支路坏了，可用无功上限其实只有额定值的60%。模型里宁可把边界设保守一点，也不要盲目相信台账，否则优化结果会在执行时被打脸。2. 二阶锥松弛的关键一步：把非凸潮流等式换成凸锥约束2.1 旋转锥与标准二阶锥的换算前文提到的非凸等式：$$ l_{ij} V_i = P_{ij}^2 + Q_{ij}^2 $$这一步的思路很直接：把等号放宽成不等号：$$ l_{ij} V_i \ge P_{ij}^2 + Q_{ij}^2 $$可行域从一条非凸曲线变成一个“碗状”凸锥。这看起来像作弊，但数学上是有边界条件的：当网络是辐射状、目标函数在变量上单调递增、且最优时该不等式取到等号，那么松弛前后的最优解完全一致。这个性质叫“紧性”（exactness）。为了喂给二阶锥求解器，还要把上面这个旋转锥换算成标准二阶锥形式。两边同时乘4，再补一个 $V_i - l_{ij}$ 的平方项：$$ 4 P_{ij}^2 + 4 Q_{ij}^2 + (V_i - l_{ij})^2 \le (V_i + l_{ij})^2 $$写成向量形式就是：$$ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ V_i - l_{ij} \end{bmatrix} \right|2 \le V_i + l{ij} $$在YALMIP里，这个形式直接用cone([2*P(k); 2*Q(k); V(i)-L(k)], V(i)+L(k))一句话写进去，非常清爽。为什么这一步能保证全局最优？因为凸优化问题的任何局部最优解都是全局最优，这是凸分析的基本结论。配电网支路潮流约束一旦变成二阶锥，整个可行域变凸，目标函数又是线性的，求解器自然可以给你一个带全局最优性保证的解。这在工程评估里特别重要：领导问“你这套方案比起粒子群好在哪”的时候，我可以直接回答“粒子群给你一个看起来不错的解，SOCP给你一个证明过的最优解”。2.2 目标函数保持线性，模型整体变成SOCP松弛之后回头看整个模型，会发现一个幸运之处：目标函数里的网损 $\sum r_{ij} l_{ij}$ 关于新变量 $l$ 是线性的；电压上下限约束 $0.95^2 \le V_i \le 1.05^2$ 是线性的；节点功率平衡方程全部是线性等式组合；调度变量边界也是线性不等式。唯一非线性的地方就是那个二阶锥约束。所以整个问题最终是一个标准SOCP：线性目标 + 线性约束 + 二阶锥约束。这种结构下的求解理论和算法都非常成熟，YALMIP + MOSEK/Gurobi/SeDuMi 都能处理。常见的困惑是“为什么不用SDP松弛”？半定规划（SDP）在配电网重构等领域确实也用，但SOCP的约束是逐条支路独立锥形，规模更小、求解更快；对辐射状网的支路潮流模型来说，SOCP已经足够，不需要更广义的SDP，尤其在节点数上千的场景里，SOCP的数值稳定性优势很明显。2.3 Matlab+YALMIP代码骨架怎么搭搭建代码骨架时，我强烈建议把“拓扑组织”和“建模”分开。否则约束里一会儿用下标、一会儿用矩阵，很容易出错。下面是一段核心骨架：%% 定义变量 P = sdpvar(nb,1); % 支路有功 Q = sdpvar(nb,1); % 支路无功 L = sdpvar(nb,1); % 支路电流平方 V = sdpvar(nnode,1); % 节点电压平方 Qc = sdpvar(nnode,1); % 节点无功补偿注入 %% 拓扑准备：先构建 parent、children % children{j,1} = 节点 j 的所有子支路编号 % parent_br(j) = 节点 j 的父支路编号，根节点为0 % 该结构需要从 from/to 数组通过 BFS 或拓扑排序生成 Con = []; Con = [Con, V(1) == 1.0]; % 根节点电压 1.0 pu Con = [Con, 0.95^2 <= V <= 1.05^2]; % 电压上下限 Con = [Con, Qc_min <= Qc <= Qc_max]; % 补偿设备边界 %% 节点功率平衡（按节点遍历，最直观） for j = 1:nnode child_br = children{j}; if j == 1 || parent_br(j) == 0 inflow_p = 0; inflow_q = 0; else fk = parent_br(j); inflow_p = P(fk) - mpc.R(fk) * L(fk); inflow_q = Q(fk) - mpc.X(fk) * L(fk); end outflow_p = sum(P(child_br)); outflow_q = sum(Q(child_br)); Con = [Con, inflow_p - outflow_p == mpc.Pgen(j) - mpc.Pd(j)]; Con = [Con, inflow_q - outflow_q == mpc.Qgen(j) - mpc.Qd(j) + Qc(j)]; end %% 电压降落方程 for j = 2:nnode fk = parent_br(j); i = mpc.from(fk); Con = [Con, V(j) == V(i) - 2*(mpc.R(fk)*P(fk) + mpc.X(fk)*Q(fk)) ... + (mpc.R(fk)^2 + mpc.X(fk)^2) * L(fk)]; end %% 二阶锥约束 for k = 1:nb i = mpc.from(k); Con = [Con, cone([2*P(k); 2*Q(k); V(i)-L(k)], V(i)+L(k))]; end %% 目标函数与求解 obj = sum(mpc.R .* L); ops = sdpsettings('solver','mosek','verbose',1); sol = optimize(Con, obj, ops);这段代码可以直接铺开来跑小算例。里面有几个细节：根节点用V(1) == 1.0，因为我们变量本身就是电压平方，所以根节点电压平方等于1.0；功率平衡里，inflow_p是父支路送到本节点的净有功，outflow_p是本节点送往所有子支路的总有功，右边是发电减负荷；补偿无功Qc放在无功平衡等式的右边，正值代表向节点注入无功。3. IEEE33节点算例搭建与Matlab代码实现细节3.1 数据准备：标幺化、拓扑脚本和从哪里修改IEEE33节点测试系统是配电网算法验证的“标准题”：33个节点、32条支路、基准电压12.66kV、总负荷3715kW+2300kvar、根节点连接变电站。这里最容易被新手坑的就是标幺化。如果直接把有名值塞进去，电压、电流、功率三个物理量量级差好几档，求解器数值稳定性和迭代速度都会受影响。我用的是 10MVA 基准功率，12.66kV 基准电压，总负荷折算下来大约 0.3715 + j0.23 pu，支路电阻电抗都在 $10^{-4} \sim 10^{-2}$ 这个量级，数值上很舒服。数据脚本我建议单独放一个case33.m，返回一个结构体：mpc.nnode = 33; mpc.nb = 32; mpc.from = [...]; % 32x1 首端节点编号 mpc.to = [...]; % 32x1 末端节点编号 mpc.R = [...]; % 32x1 支路电阻(pu) mpc.X = [...]; % 32x1 支路电抗(pu) mpc.Pd = zeros(33,1); % 节点有功负荷(pu) mpc.Qd = zeros(33,1); % 节点无功负荷(pu) mpc.Pgen = zeros(33,1); % 节点有功电源注入(pu)，根节点1通常是变电站 mpc.Qgen = zeros(33,1);支路数据顺序很重要，建议直接按尾端编号递增排列，这样用拓扑排序生成parent、children结构时不容易乱。负荷数据要放在末端节点，例如支路1-2上面的负荷放在节点2，而不是支路上。3.2 主程序流程：建模-求解-取结果主程序我把流程拆成固定的七步，日常调参基本就是在这个骨架上改参数、改约束：步骤做什么关键变量/函数1读入网络数据mpc = case33;2生成拓扑关系children,parent_br3定义优化变量sdpvar4写约束根电压、电压限值、设备边界、平衡方程、电压方程、锥约束5写目标函数obj = sum(mpc.R .* L);6调用求解器optimize(Con, obj, ops);7提取并检查结果value(P),value(L),value(V),value(Qc)第七步特别容易被忽略。很多同学跑完只看一个目标函数就完事，但value(V)的逐节点曲线、value(Qc)的配置方案才是工程上最关心的。我会在后处理里计算三个量：网损现值：loss = sum(mpc.R .* value(L));最低电压：Vmin = min(sqrt(value(V)));松弛间隙：gap = max(value(V(mpc.from(k))) .* value(L) - value(P).^2 - value(Q).^2);其中gap是判断二阶锥松弛“紧不紧”的核心证据，后面专门讲。3.3 求解器选择与结果初勘求解器选择这块，我按不同场景给出建议：求解器类型适用场景注意事项SeDuMi免费SOCP小算例验证、学生练习数值容差相对弱，模型病态时容易不收敛SDPT3免费SOCP/SDP中等算例对稀疏大模型速度一般MOSEK商业SOCP生产级应用、大规模网络数值稳健，能很好地处理大规模SOCPGurobi商业SOCP/MISOCP带整数变量的场景做混合整数二阶锥很顺手在我自己的项目里，仿真阶段首选MOSEK，因为它的容差控制和警告信息更友好；如果模型里面有电容器组投切这种0/1变量，我就切到Gurobi，省得外部循环烷法。跑通之后第一次看结果，别急着画图。先看三件事：sol.problem是不是0，松弛间隙gap是不是小于1e-6，电压曲线上有没有突变点。如果这三关过了，再往下做场景对比和参数敏感性分析。4. 松弛紧性、求解器与收敛性：调优过程中的真实坑4.1 三组仿真场景的对比结果我在33节点系统上做了三组场景，分别对应工程里最常见的三种状态：场景配置网损最低电压松弛间隙1无任何无功优化约202kW（0.0202pu）约0.903pu无2连续无功补偿优化约150kW>0.95pu<1e-73补偿 + 分布式电源无功联合优化约130kW0.97~1.03<1e-8场景1的网损在202kW左右，这是IEEE33标准算例的经典参考值。场景2我在末端电压薄弱的节点附近增加了可连续调节的无功注入，网损降了约四分之一，最低电压提到0.95以上。场景3进一步把分布式电源的无功能力纳入优化，损失更低，更重要的是白天光伏大发和夜间轻载两组工况下电压都稳住了。这些具体数值依赖于基准值和设备边界，不必迷信，但趋势是稳定的：单纯调无功对网损的改善通常能做到15%~30%，而电压合格率往往从“临界”到“稳定”，收益比网损更可观。4.2 松弛间隙检查：判断“凸化是否安全”SOCP结果出来以后，第一件事永远是看松弛间隙。所谓“紧”，可以简单理解成不等式被“顶满”到了原等式的状态；如果gap超过了1e-4，说明模型在某处被松弛得太厉害，数学上可能是“安全”的凸解，物理上却已经不忠实于真实的潮流方程。这个结果拿去工程汇报是要出问题的。我遇到过的gap变大的原因主要有三个：一是目标函数被改坏了。比如原本应该最小化网损，结果想当然加了一个电压偏差的二范数惩罚，比如1e3*sum((V-1).^2)，这个强二次项在某些边界场景下会干扰松弛紧性。解决办法是尽量用线性惩罚或者干脆不加，电压交给硬约束。二是根节点电压范围设置得过松。如果你把根节点电压也放进优化变量，还给了0.9~1.1的宽范围，那么整个系统都有了“钻空子”的空间，求解器完全可能在某个支路上用偏大的电压偏差换取更低的损耗，但实际电网根节点电压是受变电站控制的。稳妥做法是根节点直接V(1)==1.0，把这个自由度锁死。三是设备边界设得离谱。曾有一版代码我把某个补偿点上限写成2pu（对应20Mvar），本来只是测试用的宽松值，结果gap直接到1e-2。原因很直白：补偿上限宽到离谱之后，SOCP在可行域里找到一条“非物理但凸可行”的路，让某条中段支路变相“吸收”了大量虚拟无功。把设备上限改回真实的0.05pu后，gap立刻降到1e-8。这个案例后来被我写进了团队的调试手册：设备边界必须按现场台账的保守值来，越宽越危险。4.3 建模阶段容易踩的四个错误我在带团队做类似算法时，发现有些错误反复出现，值得单独列出来：第一个错误，是在YALMIP里直接写二次非凸约束，比如P(k)^2 + Q(k)^2 <= V(i)*L(k)。YALMIP能接住这个写法，但会把它当成一般非线性约束，交给内点法去硬解，既慢又不保证全局最优。正确做法是一步到位写成cone([2*P(k); 2*Q(k); V(i)-L(k)], V(i)+L(k))，把问题明确标记成二阶锥。第二个错误，是节点功率平衡漏掉支路损耗。很多示例代码会直接写成“流入=流出+负荷”，却忘了父支路送进来的功率在到达本节点时已经扣掉了r*l的损耗。漏掉这一项，优化出来的网损会比真实值低，结果对不上实际。第三个错误，是补偿设备重复配置。一次排查时我发现gap居高不下，翻数据才发现同一个节点在“新增无功节点”和“原有电容器节点”两张表里各自出现了一次，相当于给同一点塞了两组补偿变量。这类问题在真实配网拓扑数据里尤其常见，节点编号一旦合并过就会踩坑。第四个错误，是初值设置欠佳。SOCP本质上迭代求解，给一个物理合理的初值能显著减少迭代次数。我现在习惯在求解前用前推回推潮流先算一组静态工作点，然后用assign(P, value0)把初值塞给变量；即使懒得做潮流，也可以用一段全零或全一初值给求解器一个更稳的起点。4.4 性能与数值稳定性小经验33节点秒级可解，但这个模型扩展到几百上千节点时会遇到性能瓶颈。一个经验是：等节点规模大了以后，别再用一堆for循环逐条加约束，YALMIP对大规模模型的解析时间比求解时间还长。我后来把功率平衡和电压方程全部改写成稀疏矩阵形式，用A*[P;Q;L] == b一次性构建约束，速度提升非常明显。另外，ops里建议加两样东西：sdpsettings('saveyalmipmodel',1)和sdpsettings('debug',1)。前者方便排查“模型是不是非凸”，后者在YALMIP报错时会输出更详细的约束分析，基本能定位到具体是哪一行约束出问题。有一次报错说“quadratic constraint nonconvex”，我用这两个设置立刻追到了漏写锥的那一行。5. 从33节点到工程场景：扩展思路与落地思考5.1 分布式电源接入后的时序无功优化33节点只是实验室里的“验证题”，真正工程现场是分布式光伏、储能、可调负荷同时存在的动态场景。光伏出力一天内波动非常大，如果只做单断面优化，会出现设备频繁动作的问题：上午9点投一组电容器，下午2点又要求退出，一天折腾十几次，现场开关设备根本受不了。处理思路一般有三种：一是多时段联合优化，把24小时分成96个15分钟断面，目标函数里加入设备动作惩罚，让求解器自己权衡“多投一组电容省下的网损”和“动作一次带来的开关损耗”；二是滚动优化，每15分钟滚动求解一次，但把前一时段的结果作为下一时段的初始可行点，减少突变；三是给动作次数加硬约束，比如一天之内同一补偿点的动作次数不超过5次。值得多说一句，逆变器无功能力不是固定的，有功出力高时无功上限会收缩，约束关系是 $P^2 + Q^2 \le S_n^2$。这个约束本身就是一个旋转锥，二阶锥规划正好能直接表达，所以在SOCP框架里做光伏无功优化反而比传统内点法更顺手。5.2 离散设备、三相不平衡与更大规模网络再往工程走，会遇到三个更现实的问题。第一个是离散设备：分组电容器、有载调压变压器分接头都是整数变量。把0/1变量加进SOCP，模型就升级成MISOCP。如果整数变量数量不多，直接用Gurobi可以一次求到最优；如果规模太大，整数变量上百个，业界常用两层算法：外层用启发式枚举离散组合，内层还是SOCP，兼顾速度和全局搜索能力。第二个是三相不平衡。园区型配电网、农村电网很多都存在严重三相不平衡，单相模型会掩盖电压偏相问题。扩展思路是把每一相都作为独立网络建模，支路潮流从一维变量变成三维向量，约束数量翻三倍，但锥结构仍然保留。工程上更折中的做法是先做单相等效优化，再用最终配置去三相潮流里校验一遍不平衡度，这样改动小，实战中也够用。第三个是网络规模。到了几千节点，Matlab + YALMIP的建模开销会很大，此时可以考虑把建好的模型导出成MPS或LP格式，直接交给求解器批量算；或者直接切换到Python生态，用成熟的建模工具重写一版，后续对接配网系统也更方便。5.3 我从这套流程里沉淀下来的实操习惯在多次迭代和排错之后，我现在每次上新网络数据都固定执行一套校验动作，这些习惯帮我避免了很多“看起来能跑、结果不能用”的情况。第一，拓扑校验。每个非根节点有且仅有一个父节点，全图无环。不要以为数据表里写着“配电网”就真是辐射状——很多历史台账里混着联络开关闭合状态的网络，闭合后就是环网，SOCP算出来照样有解，但不是你要的辐射状结果。第二，设备容量校验。进模型前把电容柜投运状态、损坏支路、额定容量重新核对一遍，取保守可用值作为边界。这个在项目验收阶段特别重要，因为优化方案最终是要被执行机构操作的，参数写错，执行机构就乱动。第三，结果报告固定化。每次算完都输出一张“四加一”报告：网损、最低/最高电压、松弛间隙、设备执行表。评审会上，这张表比任何算法解释都更有说服力。第四，也是我最想提醒的：二阶锥规划不是万能钥匙，但它确实是目前辐射状配电网无功优化里“理论保证”和“工程效率”权衡得最好的工具。只要松弛间隙检查到位、设备边界贴合实际，这个框架从研究到落地是走得通的。希望这篇整理能帮你少走几个弯路，尤其是把那个容易被忽略的cone约束一次写对。