简介基于二阶锥规划的主动配电网最优潮流求解MATLAB代码包面向电气工程及计算机、电子信息、数学等专业学生可直接用于课程设计、期末大作业与毕业设计。代码以参数化编程为核心参数可按需修改注释清晰内置案例数据可一键运行方便快速验证。压缩包共20个文件核心为两个.m源文件另配潮流计算PPT、IEEE33节点配电网结构图、参考论文CAJ/PDF、TXT说明以及运行过程日志等整体大小10.95MB。代码针对主动配电网中含Wind、CB、SVG、OLTC、ESS等多元素的多时段24小时最优潮流问题进行二阶锥建模基于IEEE33节点系统仿真涵盖数据导入、模型构建、求解器调用、结果分析等完整流程清晰展示了配电网优化运行的要点。已有130人学习下载适合希望系统掌握OPF模型原理、MATLAB实现及配电网优化运行实践的读者深入参考。1. 为什么主动配电网最优潮流绕不开二阶锥规划从非凸到凸的一次松弛做主动配电网最优潮流最怕的不是 DG 接入数量多而是模型建好之后求解器压根不收敛。二阶锥规划SOCP这几年几乎成了配网 OPF 工程落地的默认选型它把非凸的交流潮流方程松弛成凸锥约束用 matlab 加 YALMIP 建模几十行代码就能交给内点法求解器跑出全局最优解。这篇笔记把 DistFlow 二阶锥松弛从原理到代码实现完整拆开适合正在写配电网规划、做 DG 接入评估或者要给储能多时段优化搭模型的工程师——照着复现一遍比翻十篇论文都管用。2. 二阶锥松弛做在哪一步DistFlow 模型与那条核心锥约束2.1 潮流方程里的非凸项来自哪里配电网最优潮流最麻烦的从来不是约束数量而是潮流方程本身就带强非凸项。对一条连接节点 i 和 j 的支路如果忽略对地电容DistFlow 模型可以精确写成V_j V_i - 2 * (r * P x * Q) (r^2 x^2) * I这里 V 是节点电压幅值的平方标幺值P、Q 是从 i 流向 j 的支路有功和无功I 是支路电流幅值的平方。第二段关系是电流与功率之间的约束I (P^2 Q^2) / V_i问题出在最后这个分式项上P、Q、V_i 都是待求变量P^2、Q^2 除以 V 是一个非凸的分式二次约束。传统做法是把整个潮流模型扔给牛顿法或内点法做非线性求解但 DG 接入点多、负荷波动大时初值稍微给得不好就发散。我见过不少案例IEEE 33 节点系统只加两台分布式电源传统潮流计算在某些重负荷时刻就转不动了。这就是为什么最优潮流问题需要松弛技术——把非凸等式放成凸锥把难题交给内点法。2.2 二阶锥约束的本质把等式压成不等式二阶锥规划的基本形式是min c * x s.t. || A * x || b * x约束的几何意义是变量落在以原点为顶点、方向由 b 决定的锥体里。这类问题是凸优化内点法可以在多项式时间内稳定求解不存在局部最优和全局最优的纠缠。把上面的非凸等式I (P^2 Q^2) / V_i改写成等价形式I * V_i P^2 Q^2然后做一步数学变换它就能转成标准的二阶锥|| (2*P, 2*Q, I - V_i) ||_2 I V_i这个变换的妙处在于左侧是欧几里得范数右侧是线性项完全落在 SOCP 的模型框架里。在 YALMIP 中写起来也很直观一行代码的事% 支路 k 连接节点 i 和 j constraints [constraints, norm([2*P(k); 2*Q(k); I(k) - V(i)]) I(k) V(i)];需要强调的是这里的 V(i) 是电压幅值平方不是幅值本身。很多人第一次写这个约束时直接用电压幅值去套求解器当然不认。标幺值体系下V 的正常取值范围在 0.9 到 1.1 之间这个细节也是后续结果判读的关键。2.3 为什么辐射型配电网敢做这个松弛最优解会自己落回锥面把等式松弛成不等式之后理论上说求解出来的结果可能不满足原潮流方程。但实际工程里大家照用不误核心原因是对辐射状配电网在目标函数是网损最小化的场景下最优解会自动落回锥面上。也就是说松弛出来的不等式在最优解处取等号。直观理解是这样的网损是支路电阻乘以电流平方项的和目标函数在 I 上是单调递增的。求解器在最小化目标时会倾向于把 I 压到尽可能小而左侧的锥不等式就像一个下限弹力带——I 一旦小于(P^2Q^2)/V_i就会破坏约束。两边一挤最优解处锥约束自然取等号。这就是为什么网损最小化目标天然和 SOCP 松弛是绝配。反之如果目标函数改成最大化 DG 有功出力或者最大化某个节点的电压目标的推动方向就不一定是压紧锥约束此时松弛可能不精确。这个边界必须心里有数后面第 4 章会讲怎么在跑完结果后做核验。2.4 为什么不直接上 SDPSOCP 的性价比边界半定规划SDP同样能处理配电网最优潮流的凸松弛做法是把电压外积矩阵 W 作为决策变量约束 W 半正定。模型是最精确的但对规模极不友好33 节点系统的 W 矩阵是 34×34 的对称半定矩阵上百节点的系统就接近计算瓶颈。SOCP 只需要 V、P、Q、I 四个一维向量求解时间和内存基本随网络规模线性增长对 33、123、甚至有几百个节点的馈线都跑得动。工程上我的判断标准是三相平衡、辐射状网络SOCP 精度足够只有三相不平衡、低压合环这类场景才值得动用 SDP 或多相 SOCP 叠加模型。普通 DG 接入评估和网损优化二阶锥是性价比最高的入口。3. 用 matlab 把 ADN-OPF 建模成可求解的 SOCPYALMIP 代码逐段拆解3.1 输入数据先做标幺化拓扑按父子节点重排拿到任何一份配电网算例数据第一件事不是建模而是把原始数据整理成程序可用的格式。常见算例比如经典的 33 节点系统给出的是阻抗值欧姆、负荷有功千瓦、无功千乏以及每一条支路的首端节点和末端节点。这里有个工作习惯我坚持了很久所有电气量先标幺化再进模型。主动配电网的 DG 容量和负荷量级相差很大直接混着兆瓦、千瓦、千乏写YALMIP 里的数值尺度会差三到四个数量级求解器数值稳定性会很差。%% 参数标幺化基准值必须统一在一个系统里 S_base 10e6; % 基准功率典型配网算例取 10 MVA V_base 12.66e3; % 基准电压33 节点系统常见取 12.66 kV Z_base V_base^2 / S_base; % 基准阻抗 Pd Pd_kW * 1e3 / S_base; % 有功负荷单位从 kW 转到 p.u. Qd Qd_kvar * 1e3 / S_base; R R_ohm / Z_base; % 支路电阻 X X_ohm / Z_base; % 支路电抗 %% 拓扑整理parent/children 关系是 DistFlow 递归的基础 n 33; % 节点数 parent zeros(n, 1); children cell(n, 1); for k 1:m % m 为支路数 i Branch(k, 1); % 首端节点默认朝向根节点 j Branch(k, 2); % 末端节点 parent(j) i; children{i} [children{i}, j]; end这段代码的关键是 parent 数组它把每条支路的末端节点映射到唯一父节点从而把网络拓扑转换成一个递归关系。children 是 cell 数组记录了每个节点下有哪几条孩子支路——DistFlow 里节点功率平衡的递归求和要用到它。如果算例数据本身没按朝向根节点组织先做一次深度优先搜索调整支路顺序否则后面所有等式都会错位。3.2 决策变量与 DistFlow 等式约束的 YALMIP 表达决策变量一共四组节点电压平方 V、支路有功 P、支路无功 Q、支路电流平方 I。都用 sdpvar 声明然后写 DistFlow 等式和二阶锥松弛。%% 决策变量 V sdpvar(n, 1); % 节点电压幅值平方p.u. P sdpvar(m, 1); % 支路首端有功 Q sdpvar(m, 1); % 支路首端无功 I sdpvar(m, 1); % 支路电流幅值平方 %% 支路 DistFlow 等式 二阶锥松弛 constraints []; for k 1:m i Branch(k, 1); j Branch(k, 2); r R(k); x X(k); % 电压损耗等式这是精确的线性关系不解耦 constraints [constraints, ... V(j) V(i) - 2 * (r * P(k) x * Q(k)) (r^2 x^2) * I(k)]; % 二阶锥松弛替代非凸等式 I*V_i P^2Q^2 constraints [constraints, ... norm([2*P(k); 2*Q(k); I(k) - V(i)]) I(k) V(i)]; end %% 根节点电压变电站母线当作松弛节点 constraints [constraints, V(1) 1.0]; %% 节点电压上下限注意 V 是平方 constraints [constraints, V 0.95^2, V 1.05^2]; %% 节点功率平衡DG 注入 - 负荷 支路注入 - 子支路流出 for j 2:n sum_child_P 0; sum_child_Q 0; for c children{j} sum_child_P sum_child_P P(c) - R(c) * I(c); sum_child_Q sum_child_Q Q(c) - X(c) * I(c); end % 父支路首端功率减去子支路损耗就是从节点j流出的功率 constraints [constraints, ... Pg(j) - Pd(j) P(j) - sum_child_P]; constraints [constraints, ... Qg(j) - Qd(j) Q(j) - sum_child_Q]; end这段代码里的难点是功率平衡等式。P(j) 的隐含含义是“连接父节点与节点 j 的那条支路的首端功率”它的支路编号取末端节点编号 j。子支路流出的功率要减掉子支路上的电阻损耗这就是P(c) - R(c)*I(c)的含义。初次写 DistFlow 的人最容易在这里漏掉损耗项结果做出来的潮流分布系统性偏差网损对不上。DG 的出力变量 Pg 和 Qg 我建议直接定义成 n 维 sdpvar没有 DG 的节点用 0 固定住代码比动态拼索引简单得多。3.3 目标函数怎么选网损最小、电压偏差低、DG 消纳最大目标函数是 SOCP 模型里最影响松弛精确性的部分。网损最小化是最稳妥的选择因为目标在 I 上单调递增锥约束在最优解处最容易取紧。电压偏差最小化也不错但要做成平方和损失保证凸性。%% 目标一全网有功网损最小 objective_loss sum(R .* I); %% 目标二电压偏差最小用平方形式 objective_volt sum((V - 1).^2); %% 组合目标网损为主电压修正为辅 lambda 10; % 电压项权重工程经验值通常取 5~20 objective objective_loss lambda * objective_volt;权重 lambda 不是越大越好太大时会把电压硬拉到接近 1.0网损反而上升太小时电压越限约束起主导作用目标项形同虚设。我一般先跑一版纯网损目标看电压分布再决定要不要加这个修正项。对大多数 DG 接入评估场景纯网损目标已经够用电压问题基本能被硬约束兜住。3.4 求解器选择Cplex、Gurobi、Mosek 还是 SeDuMiYALMIP 本身不求解问题它负责任务调度真正的计算交给后端求解器。二阶锥问题可以交下的求解器选择直接影响你能跑到多大规模。下表是按工程经验整理的求解器对比许可证和数值稳定性都考虑进去求解器许可证大规模 SOCP 性能备注Cplex商业授权部分高校免费优秀同时支持 MI-SOCP配网混合整数最优潮流首选Gurobi商业授权部分高校免费优秀对锥约束的预处理很强求解速度快Mosek商业授权优秀凸优化领域专精数值稳定性最好SeDuMi开源免费一般小算例能跑几百维以上明显变慢SDPT3开源免费一般适合论文复现不适合反复调参调用方式统一用 sdpsettings 切换options sdpsettings(solver, cplex, verbose, 1); sol optimize(constraints, objective, options);我的配置习惯是32 位内存的机器跑 33 节点用 SeDuMi 也可以但 123 节点以上建议直接上 Cplex 或 Gurobi。另一个细节是提前设置求解器数值容差比用默认值更稳例如 Cplex 的cplex.optim.tolerance和 Gurobi 的gurobi.MIPGap。3.5 结果回读value、开方、画电压曲线求解完成后YALMIP 的解在变量对象里需要用 value 函数取出来。这里有一个老生常谈却年年有人犯的错V 是电压平方画电压曲线前必须开根号。sol optimize(constraints, objective, options); if sol.problem 0 P_opt value(P); Q_opt value(Q); V_opt_square value(V); % 电压平方单位 p.u. V_opt sqrt(V_opt_square); % 真实电压幅值标幺值 loss_total sum(R .* value(I)); % 有功网损p.u. loss_MW loss_total * S_base / 1e6; % 折回 MW figure; plot(1:n, V_opt, -o); grid on; xlabel(节点编号); ylabel(电压幅值 (p.u.)); ylim([0.93 1.07]); end到这里一套最小可复现的主动配电网最优潮流 SOCP 求解就闭环了。跑通这份代码之后下一步才轮到参数调优和松弛精确性验证。4. 三个必调参数与松弛精确性核验不让求解器变成黑匣子4.1 电压上下限范围越小锥边界越容易被目标压紧电压上下限是第一个必调参数因为它直接影响松弛能不能取紧。常见配网约束取 0.95 到 1.05 p.u.写成代码时要记得平方constraints [constraints, V 0.95^2, V 1.05^2];如果算例对电能质量要求更高比如低压台区要求 0.93 到 1.07V 的合法域会放宽。问题在于电压范围越宽锥松弛的“活动空间”越大目标函数需要足够的推力才能把锥约束压到边界。当电压允许范围太宽且 DG 出力很大时某些支路的锥可能取不到紧结果虽然满足松弛约束但回代原潮流方程会有明显偏差。所以调参逻辑是电压范围能收紧就收紧不要给求解器留多余的松弛空间。4.2 支路电流上限与 DG 视在功率上限不设限会让解失真第二个必调参数是支路电流上限。很多初版代码完全不写这个约束求解器会把某些重载支路的电流推到正常值的好几倍网损结果自然失真。配网支路电流上限可以按导线载流量折算典型架空线在标幺体系下取 1.5 到 2.0 p.u.。constraints [constraints, I 2.0^2]; % 电流上限平方DG 的逆变器约束也类似视在功率约束要写成二阶锥% 对每个装有 DG 的节点 for j 1:n if has_DG(j) constraints [constraints, ... norm([Pg(j); Qg(j)]) Sg_max(j)]; end end如果不加这个约束PQ 解耦的 DG 模型有时会出现无功功率超过逆变器容量的事故工况结果在工程上不可接受。加了之后Pn和Qn之间的关系被限制在一个圆内正好是 SOCP 天然支持的约束形式。4.3 惩罚系数 lambda当松弛不精确时的后悔药有时候网损目标也压不紧锥一个工程上常用的补救措施是在目标函数里加一项锥约束的惩罚。这个做法用到的量是锥约束两边之差写成代码是在目标里加罚项%% 附加惩罚项促进锥约束取紧 cone_violation 0; for k 1:m i Branch(k, 1); cone_violation cone_violation ... norm([2*P(k); 2*Q(k); I(k) - V(i)]) - (I(k) V(i)); end objective objective_loss 100 * cone_violation;罚系数取 10 到 100 之间比较常见。加罚的本质是牺牲一点目标精确度换取解更贴近原潮流方程。注意罚项本身是凸的不会破坏模型性质。我一般在松弛核验发现最大锥误差超过 1e-3 时才加这个罚项。4.4 松弛精确性核验跑完这一步才敢把结果拿去写报告跑完 SOCP 不等于拿到了可用的最优潮流结果。最后一步必须做核验——把解回代原始潮流方程看锥松弛到底松了多少。%% 松弛精确性核验 V_opt_square value(V); I_opt value(I); P_opt value(P); Q_opt value(Q); % 对每条支路计算锥误差 cone_gap zeros(m, 1); for k 1:m i Branch(k, 1); lhs I_opt(k) * V_opt_square(i); % 原始等式左边 rhs P_opt(k)^2 Q_opt(k)^2; % 原始等式右边 cone_gap(k) (lhs - rhs) / (rhs 1e-6); % 相对误差 end max_gap max(cone_gap); if max_gap 1e-4 disp(锥松弛精确结果工程可用); else disp([锥松弛最大相对误差: , num2str(max_gap), 建议收紧电压范围或加罚项]); end工程经验阈值是这样掌握的最大相对锥误差小于 1e-4结果可以直接用于报告1e-4 到 1e-3 之间结果可用但需要解释大于 1e-3建议回炉调参。这一步做下来求解器就不再是黑匣子所有关键结论都有了依据。5. 避坑手册从建模到出结果的 6 个日常翻车场景5.1 现象YALMIP 报 No suitable solver for the problem我第一次跑二阶锥模型时满以为装了 Cplex 就万事大吉结果 optimize 直接提示没有可用求解器。原因YALMIP 把这类约束识别成conic类型老版本 Cplex 接口对锥约束支持不完整或者求解器路径没加进 MATLAB。这也可能出现在用norm写法但 YALMIP 版本过老的场景——有些旧版本把三维norm识别成任意范数而不是自动转成二阶锥。解决先确认求解器带完整授权且接口文件在 MATLAB 路径下然后升级 YALMIP 到最新版。另外可以用check命令确认约束类型[constr_violation, conic_class] check(constraints)。遇到残留问题显式在 sdpsettings 里指定 Cplex 的求解类型为SOCP也能绕过误判。5.2 现象电压曲线在 DG 节点附近出现锯齿状跳动33 节点系统里加了几个 DG 之后电压曲线本来应该平滑上升结果在 DG 节点附近出现明显的锯齿反复调整参数也消不掉。原因大概率是支路功率平衡等式里漏了子支路损耗项。我在 3.2 节反复强调P(c) - R(c)*I(c)这里的R(c)*I(c)一漏节点功率平衡系统性偏大DG 节点的电压就会被抬高到不正常水平。解决逐条支路核对功率平衡等式同时把网损结果和潮流计算软件的结果做对比。如果总网损比交流潮流结果低 5% 以上基本就是损耗项写错了。5.3 现象网损算出来比不接 DG 还大DG 接入后网损应该下降结果算出来反而升高了不少。原因DG 的无功上限约束没加或者无功目标处理不当。逆变器型 DG 在最大化有功出力的同时如果无功上限设得太高可能会导致大量无功在馈线上长距离流动网损自然升高。另一种可能是目标函数里加入了电压偏差项且权重 lambda 调得太大求解器为了电压好看硬压 DG 无功网损被牺牲。解决先用纯网损目标跑一版确认 DG 消纳和网损的趋势再用带无功上限的完整模型跑一版。两版趋势一致才说明问题出在目标权重而不是模型错误。5.4 现象三节点只花 0.2 秒33 节点跑了 10 分钟还在转小算例秒出结果换 33 节点就卡死这是典型的大规模化翻车。原因变量虽然全是连续的 sdpvar但二阶锥约束的数量等于支路数。如果每个支路的锥约束都用三维norm写YALMIP 会为每个约束生成一个内部锥变量33 节点系统总共几十个锥约束理论上不会这么慢。慢在数值病态——最常见的是没有标幺化阻抗和功率的数值尺度差了太多。解决回到 3.1 节检查标幺化。另外把约束从I (P^2Q^2)/V这类隐式写法全部改成norm([2P;2Q;I-V]) IV的显式锥写法求解器预处理的效率会高很多。5.5 现象储能 SOC 约束加了之后模型直接不可行多时段模型里加入储能荷电状态约束后求解器报不可行怎么调都解不出来。原因储能模型里同时写了充电状态二进制变量和放电状态二进制变量把 SOC 更新等式强行做成了混合整数约束可行域在某些时段出现断点。另一个常见坑是充放电效率不等导致 SOC 更新等式两侧量纲不一致等式无解。解决对纯连续 SOCP 模型用双向功率变量即可SOC 更新写成SOC(t1) SOC(t) eta_ch * P_ch(t) - eta_dis * P_dis(t)但 P_ch 和 P_dis 不能同时为正的约束要靠额外二进制变量那就变成 MI-SOCP。如果不想引整数变量就用净功率变量加网损效率近似。先跑连续版本确认可行再加整数变量做精细调度。5.6 现象结果里电压全是 0.9~1.0 附近的平方数从 value(V) 里拿出来的电压数值在 0.95 附近怎么看都不像正常的 1.0 p.u. 电压曲线。原因V 变量是电压幅值平方。0.95 p.u. 的电压幅值平方就是 0.90250.98 p.u. 的平方是 0.9604。拿平方值直接画图当然所有点都偏小。解决输出前统一做sqrt(value(V))然后画在 0.93 到 1.07 的纵轴范围里。这个错基本每个人都会犯一次第一次跑通代码时建议写进自己的检查清单。6. 进阶多时段 SOCP 与最后那道松弛误差核验6.1 多时段把储能 SOC 和 OLTC 档位纳入同一个 SOCP静态最优潮流跑通之后下一步通常是多时段动态最优潮流。储能和可控负荷的时序约束让问题从单时段变成跨时段耦合但 SOCP 的凸性仍然能保持。变量扩展成V(:,t)、P(:,t)的二维 sdpvarDistFlow 等式对每个时段独立成立唯一跨时段耦合的是储能 SOC 更新T 24; % 24 个时段 SOC sdpvar(n_battery, T1); SOC(:,1) 0.2; % 初始荷电状态 20% for t 1:T constraints [constraints, ... SOC(:, t1) SOC(:, t) eta_ch * P_ch(:,t) - eta_dis * P_dis(:,t), ... SOC(:, t1) 0.1, SOC(:, t1) 0.9]; end只要不引入二进制充放电状态变量整个模型仍是 SOCPCplex 和 Gurobi 都能直接解。如果OLTC 有载调压变压器的离散档位必须精确建模那就升级成 MI-SOCP——这时候求解时间从秒级跳变到分钟级这是必要代价。如果档位很多建议先用连续逼近验证模型逻辑再切整数档位。6.2 最后核验拿回代误差当模型的体检报告多时段模型跑完核验工作比单时段更要仔细。我的习惯是把每个时段的锥误差做时间序列视图看它是不是只在尖峰负荷时段变差。如果误差在负荷尖峰时段集中变大说明电压范围设置偏宽或者 DG 容量约束没发挥足够的压紧作用。一个很实用的做法是把锥误差作为目标函数加进模型里再跑一遍比较加罚前后网损结果的差值。差值小于 0.5%说明模型本身足够紧差值大说明原问题的松弛有结构性偏差单纯调系数解决不了要回去检查电压范围或目标函数选择。我自己现在跑任何配电网 SOCP 模型都会把这种回代核验写成脚本末尾的固定动作——几行计算看起来不起眼却能在结果被质疑的时候用数据兜底。多花这三十秒希望帮你在写报告前就把所有隐患拦下来。本文还有配套的精品资源点击获取