拓冰网站定制 
最近总有做配电网方向的师弟问我同一个问题二阶锥松弛做最优潮流到底怎么落地。论文里满屏的SOCP、DistFlow、松弛紧性真到了自己上手写程序经常卡在IEEE33节点数据怎么组织、YALMIP里锥约束怎么写、CPLEX为什么一直报infeasible这些最基础的事情上。我干脆把自己常用的这套MATLAB YALMIP CPLEX算例从原理到代码完整整理了一遍从DistFlow方程推导、二阶锥转换开始到IEEE33节点数据准备、完整程序实现、结果校验和常见坑位排查全程走一遍。这篇内容适合正在做配电网优化、分布式电源接入、微电网调度的研究生和工程师直接对照复现看完你能真正理解SOCP每一步在干什么而不是只会复制粘贴代码。1. 为什么配电网最优潮流需要二阶锥松弛1.1 传统交流潮流模型为什么在配电网里不好使最优潮流Optimal Power FlowOPF的本质是在满足潮流方程、电压限值、设备容量等一系列约束的前提下最小化某个目标函数最常见的就是网络损耗。输电网里我们习惯用完整的交流潮流方程包含节点电压相角、支路导纳、cos和sin三角函数这些非线性项凑在一起形成一个高度非凸的优化问题。这种非凸问题用传统非线性规划方法比如内点法去求解最头疼的就是初值敏感。初值给得好能收敛到一个还不错的局部最优解初值给得差直接发散或者收敛到明显不合理的解。对配电网来说这个问题更严重因为配电网有几个和输电网不一样的特点一是R/X比值大线路电阻大有功和无功耦合严重潮流方程的非线性更强二是网络呈辐射状节点多、分支多但没有环网三是低压配电网里负荷波动大电压沿馈线下降明显约束常常是紧的。这意味着我们需要一种能够把配电网潮流方程转化成凸约束的办法让求解器能够稳定地找到全局最优解而不是靠调初值碰运气。二阶锥松弛Second-Order Cone Programming Relaxation就是目前配电网领域应用最成熟、效果最稳定的方案之一。1.2 二阶锥松弛如何把非凸约束变成凸约束二阶锥松弛的核心思想并不复杂把潮流方程中一个难以处理的非凸二次等式放宽成一个凸的二次锥不等式。听起来抽象我拆开说。配电网最优潮流里最麻烦的非凸项来源于线路电流和功率之间的关系。具体来说一条支路上的电流平方 (l_{ij})必须等于有功平方加无功平方再除以电压平方也就是[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]这个等式P和Q是变量v也是变量分母还带着变量数学上是一个非凸等式约束。求解器面对这种约束没办法保证全局最优。SOCP的做法是把等式“松弛”成不等式[ l_{ij} \ge \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]然后再把这个不等式等价变形为标准二阶锥形式。为什么可以这样松弛因为我们的目标函数是最小化网损而网损正比于电流平方 (l_{ij})目标函数会拼命把 (l_{ij}) 往下压。既然只允许 (l_{ij}) 大于等于右边这一项最优解就会让它正好取到等号。这种情况下松弛是“紧的”松弛前后的最优解一样。当然这个“目标函数会把l压下去”的直觉并不永远成立。如果目标函数不是网损最小而是购电成本最小、DG运行成本最小或者约束条件过于宽松也可能出现松弛不紧的情况。所以做SOCP算例最后一定要回来检查松弛间隙这一点我在第4部分会专门讲怎么做。2. 算例构建IEEE33节点系统与求解环境2.1 为什么大家都在用IEEE33节点做算例IEEE33节点系统是配电网领域最经典的测试算例出自Baran和Wu在1989年发表的配网重构论文后面几乎所有配电网研究——DG接入、无功优化、网络重构、储能调度——都会拿它作为基准系统。这个系统的参数很有代表性33个节点、32条支路基准电压12.66kV总负荷约3715kW加2300kVar网络呈辐射状单馈线结构。最妙的是它的30号节点带了一个600kVar的重无功负荷导致系统在无补偿状态下的末端电压明显偏低大概在0.904p.u.左右。这意味着它在无任何调节手段时就已经接近甚至越过电压下限非常适合用来测试电压调节类算法。用IEEE33做SOCP入门还有一层好处它的公开结果非常多。比如无DG、纯辐射状运行时的网损大约是202.7kW这个数字你可以拿来验证自己写的程序是否正确。一个SOCP模型跑出来如果网损和这个数字差太远那基本可以断定代码有问题而不是算法有问题。2.2 标幺化处理与基准值选择很多初学者拿到IEEE33的原始数据就直接往约束里塞结果CPLEX报出各种数值警告甚至infeasible。问题出在哪原始数据里支路电阻是0.0922Ω这种量级负荷是100kW这种量级导纳和功率之间差了好几个数量级。纯数值上看这是典型的病态问题求解器内部的容差设置很难同时照顾到所有约束。解决办法是标幺化。针对IEEE33通常取基准容量 (S_B 10) MVA基准电压 (V_B 12.66) kV那么基准阻抗[ Z_B \frac{V_B^2}{S_B} \frac{12.66^2}{10} \approx 16.03\ \Omega ]所有支路阻抗除以16.03就得到标幺阻抗。负荷功率同理1kW相当于 (1 / (10 \times 1000) 0.0001) p.u.。电压直接用标幺电压根节点设1.0p.u.。这样整个模型里的变量基本都在0.01到1这个区间求解器处理起来非常舒服。我见过太多人忽略这一步直接用有名值建模然后花大量时间在求解器参数上调来调去。坐标统一这件事做优化建模永远排在第一位。2.3 YALMIP与CPLEX的安装和验证步骤YALMIP是一个运行在MATLAB里的免费建模工具箱作用是把优化变量、约束、目标函数用很接近数学表达式的语法写出来然后自动转换成底层求解器能识别的格式。CPLEX则是IBM的商业求解器性能稳定对二阶锥规划支持得非常好高校一般能申请到学术版授权。安装步骤其实很简单从官方渠道注册安装MATLAB这里不展开拿到授权后正常安装就行。从YALMIP官网或GitHub仓库下载最新版YALMIP解压后把整个文件夹添加到MATLAB路径中。安装CPLEX安装完成后找到CPLEX安装目录下的cplex/matlab文件夹同样添加到MATLAB路径。在MATLAB里运行yalmiptest如果列表里能看到CPLEX的状态是可用就说明环境配置成功。这里有个容易踩的坑MATLAB路径里如果同时存在多个求解器的接口YALMIP默认会按自己的优先级选择求解器。所以当你明明装了CPLEX却一直显示在用别的求解器时记得在调用代码里用sdpsettings(solver,cplex)显式指定。我下面给出的程序里就是这么处理的。3. DistFlow方程到二阶锥约束从原理到代码3.1 DistFlow方程配电网潮流计算的基本框架配电网最优潮流里最常用的潮流模型是DistFlow它专门针对辐射状网络设计避开了完整交流潮流里那些复杂的三角函数。对于一条从节点i流向节点j的支路定义变量(P_{ij})支路有功功率(Q_{ij})支路无功功率(v_i V_i^2)节点电压幅值平方(l_{ij} I_{ij}^2)支路电流幅值平方DistFlow的三大方程如下有功平衡[ \sum_{k:(j,k)} P_{jk} P_{ij} - r_{ij}l_{ij} - P_{Lj} P_{Gj} ]无功平衡[ \sum_{k:(j,k)} Q_{jk} Q_{ij} - x_{ij}l_{ij} - Q_{Lj} Q_{Gj} ]电压降落方程[ v_i - v_j 2(r_{ij}P_{ij} x_{ij}Q_{ij}) - (r_{ij}^2 x_{ij}^2)l_{ij} ]再加上电流定义式[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]配电网的线路比较短对地导纳很小DistFlow忽略充电电容是合理的。这套方程组把一个复杂的交流潮流问题简化成了只包含实数变量的多项式方程为后面的凸松弛打下了基础。3.2 从非凸等式到二阶锥不等式关键一步在哪里上面方程组里前三个方程都是线性约束优化求解器很喜欢。真正麻烦的是最后一个电流定义式它把 (P)、(Q)、(v)、(l) 四个变量用非凸的二次等式绑在一起。SOCP的转换分两步。第一步把等式改成不等式[ l_{ij} \ge \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]第二步把这个不等式等价改写成标准二阶锥形式。具体来说通过对两边进行代数变形可以得到[ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ v_i - l_{ij} \end{bmatrix} \right|2 \le v_i l{ij} ]验证很简单把两边平方展开左边 (4P^2 4Q^2 (v-l)^2)右边 ((vl)^2 v^2 2vl l^2)。两边约掉公共项最终等价于 (P^2 Q^2 \le vl)也就是 (l \ge (P^2Q^2)/v)。这个变换妙就妙在二阶锥是一个凸集合。整个非凸最优潮流问题就变成了一个凸优化问题CPLEX这类求解器可以保证收敛到全局最优解不再依赖初值。3.3 用YALMIP表示SOCP时需要注意的细节在YALMIP里写二阶锥约束有两种常见写法% 方式一用cone函数 Constraints [Constraints, cone([2*P(k); 2*Q(k); v(i)-l(k)], v(i)l(k))]; % 方式二用norm Constraints [Constraints, norm([2*P(k); 2*Q(k); v(i)-l(k)]) v(i)l(k)];两种写法数学上等价但我建议优先用cone。原因是cone能显式告诉求解器这是一个二阶锥约束求解器可以走专门的SOCP算法而norm写法虽然YALMIP也能识别但在某些版本里有额外的转换开销也更容易触发求解器内部的预处理问题。还有一个细节DistFlow里的支路方向不是随便定的。IEEE33的辐射状结构可以看作从根节点1开始向下游分支我习惯按“起点靠近根节点、终点远离根节点”的方向组织支路数据。这样 (P_{ij}) 为正表示功率从根节点侧流向末端功率平衡方程写起来不会乱。4. 完整算例MATLAB YALMIP CPLEX程序与结果分析4.1 完整程序数据准备、建模、求解、后处理下面给出一个可以直接运行的完整程序。我在里面加了一个分布式电源场景节点18和节点33各接一个容量400kW、功率因数0.9的DG优化变量是DG的有功出力目标是网损最小。这个设置比纯无DG潮流更有“最优潮流”的味道能看出优化变量在起作用。首先给出数据部分%% IEEE33节点配电网SOCP最优潮流 % 基于YALMIP CPLEX实现目标函数网损最小 clear; clc; close all; %% 1. 基准值与数据输入 SB 10; % 基准容量 MVA VB 12.66; % 基准电压 kV ZB VB^2 / SB; % 基准阻抗 Ohm % 支路数据: [起点 终点 R(ohm) X(ohm)] branch [ 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 ]; % 负荷数据: [节点 有功(kW) 无功(kVar)]根节点1无负荷 loadData [ 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 ]; nb size(branch, 1); % 支路数 nn 33; % 节点数 % 标幺化 r_pu branch(:, 3) / ZB; x_pu branch(:, 4) / ZB; PL zeros(nn, 1); QL zeros(nn, 1); for t 1:size(loadData, 1) PL(loadData(t,1)) loadData(t,2) / (SB * 1000); % kW - p.u. QL(loadData(t,1)) loadData(t,3) / (SB * 1000); % kVar - p.u. end接下来是建模和求解部分%% 2. 定义优化变量 P sdpvar(nb, 1); % 支路有功 p.u. Q sdpvar(nb, 1); % 支路无功 p.u. l sdpvar(nb, 1); % 支路电流平方 p.u. v sdpvar(nn, 1); % 节点电压平方 p.u. % DG变量接在节点18和33 nodeDG [18; 33]; Pdg sdpvar(2, 1); % DG有功出力 p.u. pmax 400 / (SB * 1000); % 400kW pf 0.9; Qdg Pdg * tan(acos(pf)); % 恒功率因数控制 %% 3. 约束条件 C []; % 根节点电压 C [C, v(1) 1]; % DG出力限值 C [C, 0 Pdg pmax]; % 电压上下限配电网一般允许0.90~1.10 Vmin 0.90; Vmax 1.10; C [C, Vmin^2 v Vmax^2]; % 支路电压降落与二阶锥约束 for k 1:nb i branch(k, 1); j branch(k, 2); C [C, v(i) - v(j) 2*(r_pu(k)*P(k) x_pu(k)*Q(k)) - (r_pu(k)^2 x_pu(k)^2)*l(k)]; C [C, cone([2*P(k); 2*Q(k); v(i)-l(k)], v(i)l(k))]; end % 节点功率平衡 for t 2:nn inflowP 0; outflowP 0; inflowQ 0; outflowQ 0; for k 1:nb i branch(k,1); j branch(k,2); if j t inflowP inflowP P(k) - r_pu(k)*l(k); inflowQ inflowQ Q(k) - x_pu(k)*l(k); end if i t outflowP outflowP P(k); outflowQ outflowQ Q(k); end end gIdx find(nodeDG t); if isempty(gIdx) C [C, inflowP - outflowP PL(t)]; C [C, inflowQ - outflowQ QL(t)]; else C [C, inflowP - outflowP PL(t) - Pdg(gIdx)]; C [C, inflowQ - outflowQ QL(t) - Qdg(gIdx)]; end end %% 4. 目标函数网损最小 Objective sum(r_pu .* l); %% 5. 求解 ops sdpsettings(solver, cplex, verbose, 2); sol optimize(C, Objective, ops); if sol.problem ~ 0 error(求解失败%s, sol.info); end %% 6. 结果后处理 loss_kW value(Objective) * SB * 1000; V sqrt(value(v)); [minV, minIdx] min(V); fprintf(网损%.2f kW\n, loss_kW); fprintf(最低电压%.4f p.u.节点%d\n, minV, minIdx); % 根节点注入有功 pinj_pu 0; for k 1:nb if branch(k, 1) 1 pinj_pu pinj_pu value(P(k)); end end fprintf(根节点注入有功%.2f kW\n, pinj_pu * SB * 1000); % 功率平衡校验 Pdg_opt value(Pdg); balance_err (pinj_pu sum(Pdg_opt) - sum(PL) - value(Objective)) * SB * 1000; fprintf(功率平衡校验误差%.4f kW\n, balance_err); % 松弛间隙校验 max_gap 0; for k 1:nb i branch(k, 1); gap value(l(k)) - (value(P(k))^2 value(Q(k))^2) / value(v(i)); if gap max_gap max_gap gap; end end fprintf(最大对偶间隙%.2e\n, max_gap); % 电压剖面 figure; plot(1:nn, V, o-, LineWidth, 1.5); grid on; xlabel(节点编号); ylabel(电压幅值 (p.u.)); title(IEEE33节点优化后电压剖面);把这个脚本按顺序保存成一份MATLAB文件配置好YALMIP和CPLEX路径后直接运行即可。代码里的注释已经比较详细了下面说一下关键部分的设计意图。4.2 结果怎么判断网损、电压剖面与松弛间隙程序跑完后你会看到终端输出几组关键数据。判断程序是否正确我建议按下面顺序来。第一步看求解器返回状态。sol.problem 0表示求解成功其他值都代表有问题具体意思可以用sol.info查看。第二步看网损。如果你把DG变量删掉只跑纯无DG的SOCP潮流网损应该是大约202.7kW最低电压大约0.904p.u.最低电压节点在节点18附近。这两个数字是IEEE33系统的公开标准结果能对上说明你的模型和数据没有问题。第三步看功率平衡校验。程序里我做了全网功率平衡检查根节点注入有功加上DG有功总出力应该等于总负荷加上全网损耗。误差在kW量级的千分之一以下都算正常。这一步能一次性排查掉大部分建模错误。第四步看对偶间隙。SOCP松弛是否精确就看每个支路上计算出的 (l_{ij}) 与 (\frac{P_{ij}^2 Q_{ij}^2}{v_i}) 的差距。最大间隙在10的负5次方以下说明松弛紧结果可信。如果间隙很大说明松弛不紧那就要回头审视目标函数和约束条件了。加上DG优化后你会看到两个现象一是网损比无DG时下降二是系统最低电压被抬升。DG的无功出力和有功出力是绑定的功率因数0.9意味着它同时向系统注入无功这对电压支撑是有帮助的。不过注入无功过多也可能导致某些节点电压逼近上限所以DG出力不一定就会冲到400kW的满发上限这正是最优潮流的魅力——所有变量都在约束边界上自动寻找最优点。4.3 如何把这段代码快速扩展到自己的课题这套程序最大的价值在于扩展性。我自己写配电网SOCP相关程序都是在这个框架上改的。想加储能就在对应节点加一个新的变量 (P_s)把储能荷电状态SOC随时间变化的约束叠加上去目标函数里再加上充放电惩罚项。想加OLTC有载调压变压器就把变压器支路的电压比作为一个新变量配合变比范围和调节代价约束。想研究三相不平衡配电网那要把单相DistFlow换成三相DistFlow变量从标量变成3x1向量但SOCP的整体框架不变。换求解器也很简单。如果你手上有Gurobi或者Mosek的授权把solver参数从cplex改成gurobi或mosek就行YALMIP会自动做语法转换。这意味着你不需要因为换求解器而重写模型。5. 常见问题与排查技巧实录5.1 求解器报错类问题速查我把自己和周围人跑这些算例时最常遇到的报错整理了一下直接看表报错或现象可能原因解决办法No suitable solverYALMIP路径里没有找到CPLEX/Gurobi检查求解器是否加入MATLAB路径用yalmiptest验证License ErrorCPLEX授权未配置或到期检查授权环境变量重新激活学术版授权Infeasible problem电压约束设得太紧或模型约束写错先放宽电压上下限再逐条排查约束Numerical trouble/NaN变量量级不统一没做标幺化全部转成标幺值检查基准容量和基准电压求解很慢二阶锥约束写得太多或太碎用cone函数而非norm关闭verbose看耗时报错cone无法处理CPLEX版本过旧升级CPLEX到12.10以上版本这里面最坑的就是Infeasible。很多初学者一看到不可行就慌了其实排查思路很清晰先把所有约束全部注释掉只留根节点电压约束和最松的设备容量约束然后一条条加回去。哪一步加上去之后问题变成不可行哪一步就是出错的地方。用二分法排通常十分钟内能定位。5.2 结果异常与模型不可行的几类典型案例第一类网损数值对不上。IEEE33无DG标准网损是202.7kW左右你要是算出个500kW或者20kW先别怀疑算法查数据。重点检查支路阻抗有没有写反、负荷单位是不是kW写成了MW、节点编号有没有对错位。我甚至见过有人把33条负荷数据粘成32条末尾的负荷整体前移结果当然全错。第二类电压下限设置不当导致不可行。IEEE33在无DG时本身末端电压就只有0.904p.u.你要是把(V_{min})设成0.95求解器直接给你报infeasible。这不是代码问题是系统本身在重负荷下达不到这么高的电压下限。这时要么放宽下限到0.90要么加DG、无功补偿装置把电压抬起来后再收紧约束。第三类DG功率因数设置导致约束冲突。恒功率因数控制下DG的有功和无功是线性绑定关系。如果DG接入点的电压上限比较紧而这个节点负荷又很小DG发功率时就会把电压推高导致模型无解。实践中要么把DG容量调小要么改成无功可调的DG模型让DG既能发无功也能吸无功。5.3 松弛不精确怎么办最后说一个进阶问题对偶间隙不收敛到零。前面说过SOCP松弛紧不紧取决于目标函数是否在推动 (l) 往下压。当目标函数不是网损最小而是DG运行成本最小、或者某些节点注入功率有直接经济成本时最优解处的锥约束就有可能不紧。处理办法有三个思路。第一个思路最简单在目标函数里加一个很小的网损惩罚项比如 (\epsilon \sum r_{ij}l_{ij})(\epsilon)取1e-4量级既不影响原目标的主次又能把松弛“拉紧”。第二个思路是求解后检查哪些支路的对偶间隙超标然后对这几条支路单独加补偿。第三个思路是如果间隙始终大说明这个问题的目标函数本质上不鼓励松弛紧可以考虑更高精度的松弛比如SDP松弛或者用凸凹过程CCP做迭代修复。从我的经验来看配电网SOCP模型中90%以上的场景松弛都是紧的真正需要修复的情况很少。但检查这一步绝不能省因为这是你的模型结果可信度的最后一道保险。我自己做这套算例踩过最深的坑是单位不统一导致的一小时无意义排错。当时负荷数据里混用了kW和MWCPLEX怎么都给出一个离谱的网损结果。后来把所有量全部转成标幺值问题瞬间就干净了。所以每次拿到一个新的系统数据我第一件事永远是确认基准容量、基准电压然后才动手写约束。如果你准备拿这套程序去改自己的DG选址、储能调度课题建议先把无DG的版本跑通、跑出202.7kW这个标准结果再往上加设备。一次只加一个变量出问题的时候也容易定位是哪一步引入的。这个习惯能帮你省下大量调试时间。 二阶锥松弛在配电网最优潮流中的MATLAB+YALMIP+CPLEX实现
最近总有做配电网方向的师弟问我同一个问题二阶锥松弛做最优潮流到底怎么落地。论文里满屏的SOCP、DistFlow、松弛紧性真到了自己上手写程序经常卡在IEEE33节点数据怎么组织、YALMIP里锥约束怎么写、CPLEX为什么一直报infeasible这些最基础的事情上。我干脆把自己常用的这套MATLAB YALMIP CPLEX算例从原理到代码完整整理了一遍从DistFlow方程推导、二阶锥转换开始到IEEE33节点数据准备、完整程序实现、结果校验和常见坑位排查全程走一遍。这篇内容适合正在做配电网优化、分布式电源接入、微电网调度的研究生和工程师直接对照复现看完你能真正理解SOCP每一步在干什么而不是只会复制粘贴代码。1. 为什么配电网最优潮流需要二阶锥松弛1.1 传统交流潮流模型为什么在配电网里不好使最优潮流Optimal Power FlowOPF的本质是在满足潮流方程、电压限值、设备容量等一系列约束的前提下最小化某个目标函数最常见的就是网络损耗。输电网里我们习惯用完整的交流潮流方程包含节点电压相角、支路导纳、cos和sin三角函数这些非线性项凑在一起形成一个高度非凸的优化问题。这种非凸问题用传统非线性规划方法比如内点法去求解最头疼的就是初值敏感。初值给得好能收敛到一个还不错的局部最优解初值给得差直接发散或者收敛到明显不合理的解。对配电网来说这个问题更严重因为配电网有几个和输电网不一样的特点一是R/X比值大线路电阻大有功和无功耦合严重潮流方程的非线性更强二是网络呈辐射状节点多、分支多但没有环网三是低压配电网里负荷波动大电压沿馈线下降明显约束常常是紧的。这意味着我们需要一种能够把配电网潮流方程转化成凸约束的办法让求解器能够稳定地找到全局最优解而不是靠调初值碰运气。二阶锥松弛Second-Order Cone Programming Relaxation就是目前配电网领域应用最成熟、效果最稳定的方案之一。1.2 二阶锥松弛如何把非凸约束变成凸约束二阶锥松弛的核心思想并不复杂把潮流方程中一个难以处理的非凸二次等式放宽成一个凸的二次锥不等式。听起来抽象我拆开说。配电网最优潮流里最麻烦的非凸项来源于线路电流和功率之间的关系。具体来说一条支路上的电流平方 (l_{ij})必须等于有功平方加无功平方再除以电压平方也就是[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]这个等式P和Q是变量v也是变量分母还带着变量数学上是一个非凸等式约束。求解器面对这种约束没办法保证全局最优。SOCP的做法是把等式“松弛”成不等式[ l_{ij} \ge \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]然后再把这个不等式等价变形为标准二阶锥形式。为什么可以这样松弛因为我们的目标函数是最小化网损而网损正比于电流平方 (l_{ij})目标函数会拼命把 (l_{ij}) 往下压。既然只允许 (l_{ij}) 大于等于右边这一项最优解就会让它正好取到等号。这种情况下松弛是“紧的”松弛前后的最优解一样。当然这个“目标函数会把l压下去”的直觉并不永远成立。如果目标函数不是网损最小而是购电成本最小、DG运行成本最小或者约束条件过于宽松也可能出现松弛不紧的情况。所以做SOCP算例最后一定要回来检查松弛间隙这一点我在第4部分会专门讲怎么做。2. 算例构建IEEE33节点系统与求解环境2.1 为什么大家都在用IEEE33节点做算例IEEE33节点系统是配电网领域最经典的测试算例出自Baran和Wu在1989年发表的配网重构论文后面几乎所有配电网研究——DG接入、无功优化、网络重构、储能调度——都会拿它作为基准系统。这个系统的参数很有代表性33个节点、32条支路基准电压12.66kV总负荷约3715kW加2300kVar网络呈辐射状单馈线结构。最妙的是它的30号节点带了一个600kVar的重无功负荷导致系统在无补偿状态下的末端电压明显偏低大概在0.904p.u.左右。这意味着它在无任何调节手段时就已经接近甚至越过电压下限非常适合用来测试电压调节类算法。用IEEE33做SOCP入门还有一层好处它的公开结果非常多。比如无DG、纯辐射状运行时的网损大约是202.7kW这个数字你可以拿来验证自己写的程序是否正确。一个SOCP模型跑出来如果网损和这个数字差太远那基本可以断定代码有问题而不是算法有问题。2.2 标幺化处理与基准值选择很多初学者拿到IEEE33的原始数据就直接往约束里塞结果CPLEX报出各种数值警告甚至infeasible。问题出在哪原始数据里支路电阻是0.0922Ω这种量级负荷是100kW这种量级导纳和功率之间差了好几个数量级。纯数值上看这是典型的病态问题求解器内部的容差设置很难同时照顾到所有约束。解决办法是标幺化。针对IEEE33通常取基准容量 (S_B 10) MVA基准电压 (V_B 12.66) kV那么基准阻抗[ Z_B \frac{V_B^2}{S_B} \frac{12.66^2}{10} \approx 16.03\ \Omega ]所有支路阻抗除以16.03就得到标幺阻抗。负荷功率同理1kW相当于 (1 / (10 \times 1000) 0.0001) p.u.。电压直接用标幺电压根节点设1.0p.u.。这样整个模型里的变量基本都在0.01到1这个区间求解器处理起来非常舒服。我见过太多人忽略这一步直接用有名值建模然后花大量时间在求解器参数上调来调去。坐标统一这件事做优化建模永远排在第一位。2.3 YALMIP与CPLEX的安装和验证步骤YALMIP是一个运行在MATLAB里的免费建模工具箱作用是把优化变量、约束、目标函数用很接近数学表达式的语法写出来然后自动转换成底层求解器能识别的格式。CPLEX则是IBM的商业求解器性能稳定对二阶锥规划支持得非常好高校一般能申请到学术版授权。安装步骤其实很简单从官方渠道注册安装MATLAB这里不展开拿到授权后正常安装就行。从YALMIP官网或GitHub仓库下载最新版YALMIP解压后把整个文件夹添加到MATLAB路径中。安装CPLEX安装完成后找到CPLEX安装目录下的cplex/matlab文件夹同样添加到MATLAB路径。在MATLAB里运行yalmiptest如果列表里能看到CPLEX的状态是可用就说明环境配置成功。这里有个容易踩的坑MATLAB路径里如果同时存在多个求解器的接口YALMIP默认会按自己的优先级选择求解器。所以当你明明装了CPLEX却一直显示在用别的求解器时记得在调用代码里用sdpsettings(solver,cplex)显式指定。我下面给出的程序里就是这么处理的。3. DistFlow方程到二阶锥约束从原理到代码3.1 DistFlow方程配电网潮流计算的基本框架配电网最优潮流里最常用的潮流模型是DistFlow它专门针对辐射状网络设计避开了完整交流潮流里那些复杂的三角函数。对于一条从节点i流向节点j的支路定义变量(P_{ij})支路有功功率(Q_{ij})支路无功功率(v_i V_i^2)节点电压幅值平方(l_{ij} I_{ij}^2)支路电流幅值平方DistFlow的三大方程如下有功平衡[ \sum_{k:(j,k)} P_{jk} P_{ij} - r_{ij}l_{ij} - P_{Lj} P_{Gj} ]无功平衡[ \sum_{k:(j,k)} Q_{jk} Q_{ij} - x_{ij}l_{ij} - Q_{Lj} Q_{Gj} ]电压降落方程[ v_i - v_j 2(r_{ij}P_{ij} x_{ij}Q_{ij}) - (r_{ij}^2 x_{ij}^2)l_{ij} ]再加上电流定义式[ l_{ij} \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]配电网的线路比较短对地导纳很小DistFlow忽略充电电容是合理的。这套方程组把一个复杂的交流潮流问题简化成了只包含实数变量的多项式方程为后面的凸松弛打下了基础。3.2 从非凸等式到二阶锥不等式关键一步在哪里上面方程组里前三个方程都是线性约束优化求解器很喜欢。真正麻烦的是最后一个电流定义式它把 (P)、(Q)、(v)、(l) 四个变量用非凸的二次等式绑在一起。SOCP的转换分两步。第一步把等式改成不等式[ l_{ij} \ge \frac{P_{ij}^2 Q_{ij}^2}{v_i} ]第二步把这个不等式等价改写成标准二阶锥形式。具体来说通过对两边进行代数变形可以得到[ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ v_i - l_{ij} \end{bmatrix} \right|2 \le v_i l{ij} ]验证很简单把两边平方展开左边 (4P^2 4Q^2 (v-l)^2)右边 ((vl)^2 v^2 2vl l^2)。两边约掉公共项最终等价于 (P^2 Q^2 \le vl)也就是 (l \ge (P^2Q^2)/v)。这个变换妙就妙在二阶锥是一个凸集合。整个非凸最优潮流问题就变成了一个凸优化问题CPLEX这类求解器可以保证收敛到全局最优解不再依赖初值。3.3 用YALMIP表示SOCP时需要注意的细节在YALMIP里写二阶锥约束有两种常见写法% 方式一用cone函数 Constraints [Constraints, cone([2*P(k); 2*Q(k); v(i)-l(k)], v(i)l(k))]; % 方式二用norm Constraints [Constraints, norm([2*P(k); 2*Q(k); v(i)-l(k)]) v(i)l(k)];两种写法数学上等价但我建议优先用cone。原因是cone能显式告诉求解器这是一个二阶锥约束求解器可以走专门的SOCP算法而norm写法虽然YALMIP也能识别但在某些版本里有额外的转换开销也更容易触发求解器内部的预处理问题。还有一个细节DistFlow里的支路方向不是随便定的。IEEE33的辐射状结构可以看作从根节点1开始向下游分支我习惯按“起点靠近根节点、终点远离根节点”的方向组织支路数据。这样 (P_{ij}) 为正表示功率从根节点侧流向末端功率平衡方程写起来不会乱。4. 完整算例MATLAB YALMIP CPLEX程序与结果分析4.1 完整程序数据准备、建模、求解、后处理下面给出一个可以直接运行的完整程序。我在里面加了一个分布式电源场景节点18和节点33各接一个容量400kW、功率因数0.9的DG优化变量是DG的有功出力目标是网损最小。这个设置比纯无DG潮流更有“最优潮流”的味道能看出优化变量在起作用。首先给出数据部分%% IEEE33节点配电网SOCP最优潮流 % 基于YALMIP CPLEX实现目标函数网损最小 clear; clc; close all; %% 1. 基准值与数据输入 SB 10; % 基准容量 MVA VB 12.66; % 基准电压 kV ZB VB^2 / SB; % 基准阻抗 Ohm % 支路数据: [起点 终点 R(ohm) X(ohm)] branch [ 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 ]; % 负荷数据: [节点 有功(kW) 无功(kVar)]根节点1无负荷 loadData [ 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 ]; nb size(branch, 1); % 支路数 nn 33; % 节点数 % 标幺化 r_pu branch(:, 3) / ZB; x_pu branch(:, 4) / ZB; PL zeros(nn, 1); QL zeros(nn, 1); for t 1:size(loadData, 1) PL(loadData(t,1)) loadData(t,2) / (SB * 1000); % kW - p.u. QL(loadData(t,1)) loadData(t,3) / (SB * 1000); % kVar - p.u. end接下来是建模和求解部分%% 2. 定义优化变量 P sdpvar(nb, 1); % 支路有功 p.u. Q sdpvar(nb, 1); % 支路无功 p.u. l sdpvar(nb, 1); % 支路电流平方 p.u. v sdpvar(nn, 1); % 节点电压平方 p.u. % DG变量接在节点18和33 nodeDG [18; 33]; Pdg sdpvar(2, 1); % DG有功出力 p.u. pmax 400 / (SB * 1000); % 400kW pf 0.9; Qdg Pdg * tan(acos(pf)); % 恒功率因数控制 %% 3. 约束条件 C []; % 根节点电压 C [C, v(1) 1]; % DG出力限值 C [C, 0 Pdg pmax]; % 电压上下限配电网一般允许0.90~1.10 Vmin 0.90; Vmax 1.10; C [C, Vmin^2 v Vmax^2]; % 支路电压降落与二阶锥约束 for k 1:nb i branch(k, 1); j branch(k, 2); C [C, v(i) - v(j) 2*(r_pu(k)*P(k) x_pu(k)*Q(k)) - (r_pu(k)^2 x_pu(k)^2)*l(k)]; C [C, cone([2*P(k); 2*Q(k); v(i)-l(k)], v(i)l(k))]; end % 节点功率平衡 for t 2:nn inflowP 0; outflowP 0; inflowQ 0; outflowQ 0; for k 1:nb i branch(k,1); j branch(k,2); if j t inflowP inflowP P(k) - r_pu(k)*l(k); inflowQ inflowQ Q(k) - x_pu(k)*l(k); end if i t outflowP outflowP P(k); outflowQ outflowQ Q(k); end end gIdx find(nodeDG t); if isempty(gIdx) C [C, inflowP - outflowP PL(t)]; C [C, inflowQ - outflowQ QL(t)]; else C [C, inflowP - outflowP PL(t) - Pdg(gIdx)]; C [C, inflowQ - outflowQ QL(t) - Qdg(gIdx)]; end end %% 4. 目标函数网损最小 Objective sum(r_pu .* l); %% 5. 求解 ops sdpsettings(solver, cplex, verbose, 2); sol optimize(C, Objective, ops); if sol.problem ~ 0 error(求解失败%s, sol.info); end %% 6. 结果后处理 loss_kW value(Objective) * SB * 1000; V sqrt(value(v)); [minV, minIdx] min(V); fprintf(网损%.2f kW\n, loss_kW); fprintf(最低电压%.4f p.u.节点%d\n, minV, minIdx); % 根节点注入有功 pinj_pu 0; for k 1:nb if branch(k, 1) 1 pinj_pu pinj_pu value(P(k)); end end fprintf(根节点注入有功%.2f kW\n, pinj_pu * SB * 1000); % 功率平衡校验 Pdg_opt value(Pdg); balance_err (pinj_pu sum(Pdg_opt) - sum(PL) - value(Objective)) * SB * 1000; fprintf(功率平衡校验误差%.4f kW\n, balance_err); % 松弛间隙校验 max_gap 0; for k 1:nb i branch(k, 1); gap value(l(k)) - (value(P(k))^2 value(Q(k))^2) / value(v(i)); if gap max_gap max_gap gap; end end fprintf(最大对偶间隙%.2e\n, max_gap); % 电压剖面 figure; plot(1:nn, V, o-, LineWidth, 1.5); grid on; xlabel(节点编号); ylabel(电压幅值 (p.u.)); title(IEEE33节点优化后电压剖面);把这个脚本按顺序保存成一份MATLAB文件配置好YALMIP和CPLEX路径后直接运行即可。代码里的注释已经比较详细了下面说一下关键部分的设计意图。4.2 结果怎么判断网损、电压剖面与松弛间隙程序跑完后你会看到终端输出几组关键数据。判断程序是否正确我建议按下面顺序来。第一步看求解器返回状态。sol.problem 0表示求解成功其他值都代表有问题具体意思可以用sol.info查看。第二步看网损。如果你把DG变量删掉只跑纯无DG的SOCP潮流网损应该是大约202.7kW最低电压大约0.904p.u.最低电压节点在节点18附近。这两个数字是IEEE33系统的公开标准结果能对上说明你的模型和数据没有问题。第三步看功率平衡校验。程序里我做了全网功率平衡检查根节点注入有功加上DG有功总出力应该等于总负荷加上全网损耗。误差在kW量级的千分之一以下都算正常。这一步能一次性排查掉大部分建模错误。第四步看对偶间隙。SOCP松弛是否精确就看每个支路上计算出的 (l_{ij}) 与 (\frac{P_{ij}^2 Q_{ij}^2}{v_i}) 的差距。最大间隙在10的负5次方以下说明松弛紧结果可信。如果间隙很大说明松弛不紧那就要回头审视目标函数和约束条件了。加上DG优化后你会看到两个现象一是网损比无DG时下降二是系统最低电压被抬升。DG的无功出力和有功出力是绑定的功率因数0.9意味着它同时向系统注入无功这对电压支撑是有帮助的。不过注入无功过多也可能导致某些节点电压逼近上限所以DG出力不一定就会冲到400kW的满发上限这正是最优潮流的魅力——所有变量都在约束边界上自动寻找最优点。4.3 如何把这段代码快速扩展到自己的课题这套程序最大的价值在于扩展性。我自己写配电网SOCP相关程序都是在这个框架上改的。想加储能就在对应节点加一个新的变量 (P_s)把储能荷电状态SOC随时间变化的约束叠加上去目标函数里再加上充放电惩罚项。想加OLTC有载调压变压器就把变压器支路的电压比作为一个新变量配合变比范围和调节代价约束。想研究三相不平衡配电网那要把单相DistFlow换成三相DistFlow变量从标量变成3x1向量但SOCP的整体框架不变。换求解器也很简单。如果你手上有Gurobi或者Mosek的授权把solver参数从cplex改成gurobi或mosek就行YALMIP会自动做语法转换。这意味着你不需要因为换求解器而重写模型。5. 常见问题与排查技巧实录5.1 求解器报错类问题速查我把自己和周围人跑这些算例时最常遇到的报错整理了一下直接看表报错或现象可能原因解决办法No suitable solverYALMIP路径里没有找到CPLEX/Gurobi检查求解器是否加入MATLAB路径用yalmiptest验证License ErrorCPLEX授权未配置或到期检查授权环境变量重新激活学术版授权Infeasible problem电压约束设得太紧或模型约束写错先放宽电压上下限再逐条排查约束Numerical trouble/NaN变量量级不统一没做标幺化全部转成标幺值检查基准容量和基准电压求解很慢二阶锥约束写得太多或太碎用cone函数而非norm关闭verbose看耗时报错cone无法处理CPLEX版本过旧升级CPLEX到12.10以上版本这里面最坑的就是Infeasible。很多初学者一看到不可行就慌了其实排查思路很清晰先把所有约束全部注释掉只留根节点电压约束和最松的设备容量约束然后一条条加回去。哪一步加上去之后问题变成不可行哪一步就是出错的地方。用二分法排通常十分钟内能定位。5.2 结果异常与模型不可行的几类典型案例第一类网损数值对不上。IEEE33无DG标准网损是202.7kW左右你要是算出个500kW或者20kW先别怀疑算法查数据。重点检查支路阻抗有没有写反、负荷单位是不是kW写成了MW、节点编号有没有对错位。我甚至见过有人把33条负荷数据粘成32条末尾的负荷整体前移结果当然全错。第二类电压下限设置不当导致不可行。IEEE33在无DG时本身末端电压就只有0.904p.u.你要是把(V_{min})设成0.95求解器直接给你报infeasible。这不是代码问题是系统本身在重负荷下达不到这么高的电压下限。这时要么放宽下限到0.90要么加DG、无功补偿装置把电压抬起来后再收紧约束。第三类DG功率因数设置导致约束冲突。恒功率因数控制下DG的有功和无功是线性绑定关系。如果DG接入点的电压上限比较紧而这个节点负荷又很小DG发功率时就会把电压推高导致模型无解。实践中要么把DG容量调小要么改成无功可调的DG模型让DG既能发无功也能吸无功。5.3 松弛不精确怎么办最后说一个进阶问题对偶间隙不收敛到零。前面说过SOCP松弛紧不紧取决于目标函数是否在推动 (l) 往下压。当目标函数不是网损最小而是DG运行成本最小、或者某些节点注入功率有直接经济成本时最优解处的锥约束就有可能不紧。处理办法有三个思路。第一个思路最简单在目标函数里加一个很小的网损惩罚项比如 (\epsilon \sum r_{ij}l_{ij})(\epsilon)取1e-4量级既不影响原目标的主次又能把松弛“拉紧”。第二个思路是求解后检查哪些支路的对偶间隙超标然后对这几条支路单独加补偿。第三个思路是如果间隙始终大说明这个问题的目标函数本质上不鼓励松弛紧可以考虑更高精度的松弛比如SDP松弛或者用凸凹过程CCP做迭代修复。从我的经验来看配电网SOCP模型中90%以上的场景松弛都是紧的真正需要修复的情况很少。但检查这一步绝不能省因为这是你的模型结果可信度的最后一道保险。我自己做这套算例踩过最深的坑是单位不统一导致的一小时无意义排错。当时负荷数据里混用了kW和MWCPLEX怎么都给出一个离谱的网损结果。后来把所有量全部转成标幺值问题瞬间就干净了。所以每次拿到一个新的系统数据我第一件事永远是确认基准容量、基准电压然后才动手写约束。如果你准备拿这套程序去改自己的DG选址、储能调度课题建议先把无DG的版本跑通、跑出202.7kW这个标准结果再往上加设备。一次只加一个变量出问题的时候也容易定位是哪一步引入的。这个习惯能帮你省下大量调试时间。