前阵子有位做配电网调度的同行找我一起看一个优化结果，程序跑了一个通宵，最后给出的无功补偿方案换个初始值就完全变了样。这个现象我见得太多了——配电网最优潮流（Optimal Power Flow）里的潮流方程本身就是强非线性约束，加上电压上下限、线路容量这些边界条件之后，整个可行域完全不是凸集，传统求解器很容易陷到局部最优里，解的质量全看运气。后来我帮他换了一套二阶锥松弛（SOCP）的写法，同一个优化目标在Matlab里用YALMIP建模，几十秒就能拿到全局最优解，结果拿去和潮流计算逐点核对也完全吻合。
这篇文章想把从问题建模、数学推导到代码实现、再到松弛精确性验证的完整流程拆开讲清楚。内容适配正在做配电网规划、分布式电源接入分析的研究生，也适合正在纠结无功电压控制切入点、自己搭优化模型做决策支持的工程师。我会尽量把“为什么要这样做”讲透，而不是丢一段能跑但看不懂的代码了事。
1. 配电网最优潮流：问题背景与“非凸”痛点
1.1 配电网最优潮流和输电网的最优潮流差在哪
输电网最优潮流大家比较熟悉，目标函数通常是发电成本最小化，约束是交流潮流方程、机组出力和线路限额。这套东西在输电网里虽然也是非线性的，但一来网络结构是环网、R/X比小，潮流方程的病态程度相对有限；二来输电网的调节手段集中在大型机组和高抗高容上，变量规整，传统内点法往往能顶得住。
配电网最优潮流就完全是另一回事了。配网是中低压、辐射状结构为主，线路电阻和电感的比值R/X往往接近甚至超过1，这意味着有功和无功的耦合非常强，不能像输电网那样做P-Q解耦近似。同时负荷分布在成百上千个节点上，三相不平衡、分布式光伏大量接入，调压手段分散在变压器分接头、并联电容器、静止无功补偿器以及逆变器无功能力上。配网又几乎没有高精度的实时量测，历史数据与在线数据的质量参差不齐。
这直接带来一个现实问题：把配电网最优潮流写成数学规划之后，约束里大量出现电压、电流、功率的乘积项和二次项，目标函数又经常同时包含网损、电压偏差、购电成本多个分量。这样一个非凸、大规模、强耦合的优化问题，传统数值方法很难稳定求解。
1.2 非凸可行域带来的两个麻烦：局部最优和求解器不适配
先想清楚“非凸”到底意味着什么。一个约束如果描述的是一个非凸集合，那么可行域可能像一片七拐八弯的地形，有多个局部最低点。这时候无论用牛顿法、序列二次规划还是内点法，本质上都是从某个初值出发沿梯度方向“往下走”，走到哪个谷底算哪个谷底。初值给得不好，收敛到全局最优的概率就很低。更麻烦的是，配电网潮流方程在重负荷情况下雅可比矩阵容易接近病态，迭代过程可能震荡甚至发散，很多人把这种现象误判成“这个系统无解”，其实只是算法本身不行。
另一个麻烦是求解器不适配。通用的非线性优化求解器可以处理非凸问题，但给不出全局最优性保证；基于全局搜索的算法（比如分支定界）又受限于大规模变量的组合爆炸，算到节点数上百就已经很难满足了。工程上需要的是一个既能在多项式时间求解、又能给出最优性保证的模型。
所以在配电网最优潮流这个方向上，近几年最主流的一条路就是“凸松弛”——把原本非凸的问题，通过某种技巧松弛成一个凸优化问题。凸问题的好处是，只要找到局部最优就是全局最优，并且有大量成熟的内点法求解器可以可靠求解。二阶锥松弛（SOCP）就是其中应用最广、效果最稳定的一种。
1.3 为什么要选二阶锥松弛：三条凸化路线对比
配电网最优潮流的凸化路线大致有三条：线性化近似、半定松弛（SDP）和二阶锥松弛（SOCP）。线性化近似最省事，把潮流方程里的高阶项直接忽略或者用一阶泰勒展开，得到的是线性约束，求解速度最快。但代价是精度随负荷加重而迅速恶化，尤其在配网R/X比较高、电压降落明显的时候，偏差大到没法工程使用。
SDP松弛精度最高，理论上对很多配电网情形都能做到精确松弛，但它把问题放到了矩阵变量空间，变量数量从向量级变成矩阵级，规模化之后内存和计算量增长极快，跑一个几百节点的三相不平衡系统就非常吃力。
SOCP居中：把最难的非凸约束松弛成一簇二阶锥约束，变量仍然是向量级，求解速度和内存消耗都比SDP友好得多，在辐射状配电网结构下，大量文献和工程实践都证明它的松弛结果是精确的。对比下来，SOCP是“精度—复杂度—工程实现”三点上最平衡的选择，也是YALMIP、CVX这类建模工具支持度最好的问题类型之一，Matlab生态里上手成本很低。
2. DistFlow支路潮流方程：SOCP松弛的前置基础
2.1 从功率平衡到支路潮流方程
要理解SOCP松弛在配电网里为什么能成立，得先熟悉DistFlow这套支路功率方程。很多近似潮流方法上来就是潮流计算或者牛顿迭代，但DistFlow的出发点不太一样：它对每条支路，把首端流出的有功、无功、电流平方和电压平方之间的关系写成递推式，顺着辐射状的拓扑关系一路推下去，最后耦合到根节点的电压参考值上。
以辐射状配电网里一条支路(i,j)为例，变量包括支路有功P_ij、无功Q_ij，节点i的电压幅值平方v_i，支路电流幅值平方l_ij。典型的一组DistFlow方程可以写成：
- v_j = v_i - 2(r_ij P_ij + x_ij Q_ij) + (r_ij^2 + x_ij^2) * l_ij
- l_ij = (P_ij^2 + Q_ij^2) / v_i
第一式描述的是电压降落关系，r_ij 和 x_ij 分别是支路的电阻和电抗，这个式子本身就是线路压降物理规律的精确表达。第二式把电流平方和功率、电压联系起来，是欧姆定律的直接推论。
再加上每个节点的功率平衡约束，即从上游支路进入的功率减去流向所有下游支路的功率，等于该节点净负荷（负荷减去发电），一个辐射状配电网的整体潮流关系就闭合了。这套方程的结构比传统“节点注入—导纳矩阵”的形式更适合做优化建模，因为它只用每条支路自身的变量就可以写全约束，不会出现全局稠密的雅可比矩阵。
2.2 唯一非凸项：电压、电流平方的恒等关系
如果把上面两个方程拆开看，第一式是线性的，节点功率平衡也是线性的。唯一的非线性来自第二式：
l_ij = (P_ij^2 + Q_ij^2) / v_i
这里v_i、l_ij、P_ij、Q_ij都是有界的连续变量。表面上只是一个除法，实际上它是一个双线性等式约束：l_ij * v_i = P_ij^2 + Q_ij^2。而左边是l_ij和v_i的乘积，两个变量乘在一起，这种等式定义出来的集合是非凸的。
为什么非凸？画个直观一点的图：把满足等式的点集想象成空间中一张弯折的曲面，曲面上任意两点连线不一定还落在曲面上，这就是非凸集合的典型特征。求解器在这样一个曲面约束下找最优，自然容易被局部最优坑住。
很多人在这一步会卡住，因为如果直接把第二式替换成近似式（比如把v_i近似看作常数），方程的精度就没了；如果保留等式，问题又非凸。二阶锥松弛的思路只改非常小的一步：把这个等式约束放宽成一个不等式方向。
2.3 松弛成二阶锥：数学形式与直观理解
二阶锥松弛的做法是允许l_ij比物理必需值更大，即把等式放宽为：
l_ij ≥ (P_ij^2 + Q_ij^2) / v_i
等价变形一下，两边同乘v_i（v_i在物理上有下界且大于0，所以这个操作安全）：
l_ij * v_i ≥ P_ij^2 + Q_ij^2
再把它整理成标准二阶锥形式：
|| [2P_ij; 2Q_ij; l_ij - v_i] || ≤ l_ij + v_i
这条约束限定的集合就是一个旋转二阶锥，它是凸的。因此在优化里，把非凸等式捅出一个“上方空间”后，可行域彻底变凸，问题可解性大幅提升。
直观理解也不难。l_ij是电流幅值平方，它越大代表线路损耗越大，物理上更“浪费”。松弛允许这个值不必精确等于欧姆定律给出的结果，而是可以偏大。偏大的电流平方在网损里是代价项，目标函数只要对网损单调增加，求解器为了把目标压到最小，会自动把每个l_ij往小里压，最后压到锥边界上，松弛约束就重新“绷紧”成等式。这也是为什么松弛在最小网损目标下通常不会丢精度。
3. Matlab+YALMIP搭建SOCP配电网最优潮流模型
3.1 环境准备：YALMIP安装与求解器选择
做SOCP最优潮流，我建议直接在Matlab里装YALMIP建模层，它是目前最省心的建模工具，能以接近数学表达式的方式直接写约束和目标函数。至于底层求解器，选择空间比较大：开源免费的可以用Sedumi、SDPT3，嵌入式速度可以看ECOS和SCS，商业求解器在内存和速度上更稳，但授权成本高，个人学习和学术验证用开源的那几个已经能覆盖绝大多数场景。
安装步骤不复杂：下载YALMIP源码，把文件夹加到Matlab路径，然后在命令行里用yalmiptest跑一下自检，能通过就说明环境基本就绪。最容易踩的坑有两个：一是Matlab版本太老导致YALMIP某些新函数不识别，尽量用近几年的版本；二是不要同时装多个求解器版本，YALMIP有时候会优先调到一个不兼容的旧版求解器，直接报错找不到函数，白白浪费时间。
3.2 输入数据组织：基准值、线路参数与负荷
做优化前第一件事不是写约束，而是把数据整理成标幺值。标幺化的基准通常取配电网额定线电压作为电压基准，取100 kVA或者1 MVA作为容量基准。线路电阻电抗都除以基准阻抗，负荷功率除以基准容量。标幺化这一步能在数值上把电压、功率、电流拉到接近1的数量级，避免出现10^4和10^-4这种跨度极大的数值，影响内点法的收敛精度。
数据组织建议按节点表、支路表、松弛节点三块来。节点表要包含节点编号、有功负荷、无功负荷；支路表要包含首端节点、末端节点、电阻标幺值、电抗标幺值、线路容量上限；松弛节点就是变电站母线侧，通常电压幅值固定取1.0 p.u.，作为全网的电压参考。如果后续要加分布式电源，再单独准备一份DG参数表，标明接入节点、有功上限、无功上下限。
我见过不少人跳过归一化直接拿有名值建模，结果YALMIP给出的警告信息非常奇怪，比如约束互相冲突、求解器误报不可行。排到最后基本都是数值尺度问题。所以这一节千万别省。
3.3 核心代码：变量、约束、目标函数与求解
下面这段是核心建模代码的主体框架。假设你已经把支路首端数组fb、末端数组tb、电阻数组rPU、电抗数组xPU、节点负荷Pd和Qd都定义好了，所有变量都是标幺值：
% 决策变量： % v(i) 表示节点i电压幅值的平方 % l(k) 表示支路k电流幅值的平方 % P(k) 和 Q(k) 表示支路k首端流出的有功、无功 v = sdpvar(nb, 1); l = sdpvar(nl, 1); P = sdpvar(nl, 1); Q = sdpvar(nl, 1); Constraints = []; % 电压上下限约束：0.95^2 <= v <= 1.05^2 Constraints = [Constraints, Vmin <= v <= Vmax]; for k = 1:nl i = fb(k); % 首端节点 j = tb(k); % 末端节点 % 电压降落方程 Constraints = [Constraints, ... v(j) == v(i) - 2*(rPU(k)*P(k) + xPU(k)*Q(k)) + ... (rPU(k)^2 + xPU(k)^2)*l(k)]; % 二阶锥松弛约束：l(k) >= (P(k)^2 + Q(k)^2)/v(i) % 等价写法 ||[2P; 2Q; l-v]|| <= l+v Constraints = [Constraints, ... norm([2*P(k); 2*Q(k); l(k) - v(i)]) <= l(k) + v(i)]; end注意这里的v是电压幅值的平方，不是幅值本身。很多初学者在写norm([...]) <= l + v的时候，会下意识把v当成幅值平方根，但DistFlow推导里v就是平方量，约束形式是完全一致的，不要混淆。节点功率平衡部分，可以用节点—支路关联矩阵来组织，也可以直接循环累加：
% 节点功率平衡：流入功率 - 流出功率 + 发电 = 负荷 % 这里简化写法：把发电变量留空，负荷直接由根节点供给 Pbal = zeros(nb, 1); Qbal = zeros(nb, 1); for k = 1:nl i = fb(k); j = tb(k); Pbal(i) = Pbal(i) - P(k); Pbal(j) = Pbal(j) + P(k) - rPU(k)*l(k); Qbal(i) = Qbal(i) - Q(k); Qbal(j) = Qbal(j) + Q(k) - xPU(k)*l(k); end % 除去松弛节点外，各节点注入功率应等于负荷 Constraints = [Constraints, Pbal(2:end) == Pd(2:end)]; Constraints = [Constraints, Qbal(2:end) == Qd(2:end)];目标函数先以网损最小为例：
objective = sum(rPU .* l); % 线路有功损耗总和 options = sdpsettings('verbose', 1, 'solver', 'sedumi'); sol = optimize(Constraints, objective, options); % 检查求解状态 if sol.problem == 0 disp('求解成功'); v_opt = value(v); l_opt = value(l); P_opt = value(P); Q_opt = value(Q); else disp(['求解失败: ', sol.info]); end这串代码跑通了，就意味着你已经把一个非凸的配电网最优潮流问题转换成了凸问题并求出全局最优解。从节点数10到500，计算时间通常在秒级到分钟级，具体取决于支路数和所配求解器。
3.4 结果提取与可视化：电压曲线、网损、迭代信息
求解完别急着收工，先把结果转成可检查的物理量。电压幅值就是sqrt(v_opt)，每个节点的实际电压可以直接画一条曲线，看有没有越限。支路潮流和电流平方也可以直接控制台打印。
我习惯把求解器返回的sol.solveroutput一起看一下，里面记录了迭代次数和对偶间隙等关键信息。如果sol.problem是0但电压曲线在某个节点出现明显凹陷，多半是负荷数据或者支路参数标幺化错了，回过头去对一遍数据表比改模型更快。可视化方面不用太花哨，用plot(sqrt(v_opt))画一条带节点编号的电压剖面图就够了，现场汇报的时候最直观。
4. 松弛精确性怎么验证：不能只看目标函数值
4.1 为什么松弛解可能是“假”的
前面说了，二阶锥松弛的本质是把等式放宽成了不等式。优化器找到一个可行的SOCP解，但这条解对应的l值可能真的比物理所需的l大，也就是松弛约束没有被“绷紧”。这种情况下你算出来的网损就可能大于真实网损，电压曲线也可能有偏差。问题在于，目标函数值看起来非常漂亮，可它对应的物理场景并不存在。
最典型的现象发生在目标函数不完全依赖网损的时候。比如目标函数里只有电压偏差项：
min Σ (v_i - 1)^2
这时候l不进入目标函数，松弛出来的l可以任意上浮而不受惩罚，解出来的电压虽然很整齐，但电流平方和电压平方的恒等关系早就被撑开了，不是真实解。所以SOCP松弛必须配一个“松弛精确性检查”，没有这一步，结果就只能算“半成品”。
4.2 验证方法：最大松弛间隙计算
松弛间隙的定义很简单，逐条支路检查二阶锥约束的等式还原程度。理想情况下松弛约束取等号，也就是：
|| [2P_k; 2Q_k; l_k - v_i] || == l_k + v_i
数值上不可能绝对相等，因此定义一个间隙指标：
gap_k = | ||[2P_k; 2Q_k; l_k - v_i]|| - (l_k + v_i) | / max(l_k + v_i, 1)
gap_k越接近0说明该支路松弛越紧。整个网络的松弛质量就看所有支路gap_k的最大值。我自己的经验阈值是：max_gap < 1e-4 基本可以认为松弛精确；如果到了1e-2以上，就必须处理了。
把这段检查写在求解之后，算是整个代码里性价比最高的一道防线：
gap = zeros(nl, 1); for k = 1:nl i = fb(k); gap(k) = abs(norm([2*P_opt(k); 2*Q_opt(k); l_opt(k) - v_opt(i)]) - ... (l_opt(k) + v_opt(i))) / max(l_opt(k) + v_opt(i), 1); end max_gap = max(gap);如果max_gap在1e-4以下，放心往下走；如果超了很多，先别急着调权，去检查原来的非凸问题是否真的可行、负荷是否接近极限运行点。有些时候松弛不精确恰恰说明原问题可行域已经被逼到很边缘了。
4.3 松弛不精确的典型场景与处理策略
处理不精确，工程上最常用的做法是给目标函数加一个很小的网损惩罚项，比如：
min Σ r_k * l_k + α * Σ l_k
其中α取一个相对很小的正数，例如1e-4左右。惩罚项会迫使l尽可能小，从而把松弛压回边界。代价是目标函数本身被扰动了一点，但通常影响可以忽略。
另一种思路是在模型里显式锁定几个关键支路的松弛等式，把原本的SOCP约束拆成两个方向的不等式（限制上限和下限），但这等于把问题重新推回非凸，不到万不得已不要用。多数情况下，只要目标函数与网损正向相关，或者加一个足够小的惩罚，SOCP松弛在辐射状配电网里几乎总是精确的。万一加惩罚后解仍在锥内部，那就要怀疑是不是模型本身错了，比如负荷方向或线路参数不对劲，排查优先级反而更高。
5. 实操拓展：分布式电源接入、无功补偿与三相扩展
5.1 给模型加分布式光伏：变量与约束写在哪
实际项目里很少有不带分布式光伏的配电网。给上面的SOCP模型加光伏，做法是给光伏接入节点新增加入变量，并把它写进节点功率平衡。比如节点i接入一个有功上限150 kW的光伏逆变器，同时具备无功调节能力：
pg = sdpvar(nb, 1); % 节点有功注入，发电机为正 qg = sdpvar(nb, 1); % 节点无功注入 % 光伏逆变器容量约束 Constraints = [Constraints, pg(i)^2 + qg(i)^2 <= S_pv(i)^2]; Constraints = [Constraints, 0 <= pg(i) <= P_pv_max(i)];容量球约束本身就是一个二阶锥，正好和SOCP框架无缝衔接。目标函数可以改成购电成本最小，也就是根节点的有功功率乘以电价，再加上光伏逆变器的无功调节成本。这样模型就从“纯网损”变成了一个更贴近实际的小型能量管理问题。
建议把多个DG接入节点的容量、限额都放在一个结构体里，循环添加约束，避免手动逐节点写死。代码跑通后再加一个“DG出力是否越限”的后处理，能快速定位是数据问题还是模型可行性问题。
5.2 储能与连续无功补偿设备的建模
储能本质上是在节点功率平衡方程里多两个符号相反的注入变量，一个充电功率，一个放电功率。如果不要求同时充放电，可以在优化里引入充电和放电两个非负变量，再用一个二元变量把它们隔开，但这会引入整数变量，变成混合整数二阶锥规划（MISOCP），计算量会上一个台阶。
如果仅仅是做日前计划，且储能容量比单时段充放电功率大得多，可以先忽略“同一天不能既充又放”的细节，用净功率变量来近似，模型保持纯SOCP，计算速度快，结果用于调度也够用。另外，配电网里常见的SVG连续无功补偿设备更好处理，直接增加一个无功注入变量，加上下限约束即可。
5.3 三相不平衡配电网的扩展思路与模型局限
单相模型的SOCP配电网最优潮流已经很成熟，但真实配电网绝大多数是三相不平衡的。要处理三相系统，最简单的办法是把每个节点展开成A、B、C三相三个节点，支路参数变成3x3的相阻抗矩阵，然后对每一相分别建立DistFlow方程和SOCP约束。原理上并没有新的难点，但变量数直接翻三倍，线路互联的耦合支路约束也会多出来，调试成本明显上升。
这里要提醒一句：三相SOCP松弛的精确性条件比单相更严格，目标函数如果完全不涉及网损，松弛出现不精确的概率更高。所以三相工程模型里更要贯彻前面的松弛间隙检查流程。
5.4 性能调优与工程落地中的坑
算例规模一大，求解时间就会成为瓶颈。一个实用经验是先用Sedumi这类免费求解器跑通小规模模型，确认逻辑无误之后，再切到商业求解器跑大规模算例，能明显看到求解速度提升。YALMIP切换求解器只需要改sdpsettings里的solver参数，建模代码完全不用动，非常方便。
工程落地中最影响成败的其实不是求解器，而是数据质量。配电网历史负荷数据缺失严重、拓扑档案更新不及时的情况很常见。有一种做法是直接拿SCADA量测数据校核完之后做优化，量测不到的地方用伪量测补，但这属于另一个话题。模型层面能把约束写对、松弛间隙验好，就已经解决了一大半问题。
最后提个实在建议：把松弛间隙检查写成求解后的固定输出，每次跑完都看一眼。我在实际项目里吃过亏，有次把目标函数换成电压偏差平方最小之后，解出来的电压曲线很漂亮，结果发现支路电流平方和电压平方的恒等关系明显不紧，如果直接拿这个结果去下指令，大概率会失真。后来养成了先查间隙的习惯，这个问题再没出现过。二阶锥松弛不是银弹，但它目前是配电网最优潮流里最值得优先尝试的凸化工具，值得长期留在自己的工具链里。
