配电网最优潮流OPF最近在配电自动化群里被聊得非常多。手里有分布式光伏、储能或者充电桩别管是做并网方案还是日内调度最后几乎都要落到同一个数学问题在网络允许的运行范围内找一组设备出力让网损、电压偏差或者购电成本尽可能小。回答这个问题的工具就是最优潮流计算。不过我接触过不少同行以为把Matlab里的fmincon或者粒子群拿过来把交流潮流约束当等式丢进去就能算。真算起来就发现问题了——模型非线性、非凸初值一换结果跟着变算半天还不一定找到全局最优。后来我在项目里换成了二阶锥松弛配电网模型直接写成凸优化Matlab调用YALMIP加免费求解器就能稳定求解计算规模大了也不慌最优解还能给出运行层面的明确物理含义。这篇就把整个思路、建模和Matlab实现完整写一遍重点讲清楚松弛在做什么、代码怎么落、坑在哪里。1. 配电网最优潮流到底难在哪1.1 算的不是传统潮流而是在网络里找最优运行点传统潮流计算是“给定负荷和发电求全网电压、功率分布”本质是解一组方程。但最优潮流不一样它多了一个“优化”的动作哪些分布式电源出力可调、无功设备补偿多少、能不能通过改变根节点电压改善末端电压、要不要调整联络开关开合状态这些都可以作为决策变量参与求解。所以它的数学结构是“在等式和不等式约束下求目标函数最小值”。这里面的等式约束不是别的正是潮流方程本身。很多人第一步就容易忽略潮流方程本来是非线性的直接嵌进优化问题之后整个可行域就变成了弯弯曲曲的非凸集合。非凸带来的后果很直接——你没法保证找到的解是全局最优可能一个看似合理的初值就把你带进局部极值里出不来。传统输电网做最优潮流早期常用内点法硬解非线性规划对初值、矩阵条件和计算规模都有讲究。配电网的网络规模虽然不如输电网大但分支多、R/X大、三相不平衡问题也常见反而更需要一套理论上更干净的方法。二阶锥松弛就是在“精确描述物理规律”和“保持凸优化可解性”之间找一个工程上最能接受的平衡点。1.2 很多人没想明白的一个点为什么潮流等式会让问题变凸性困难我先给个直观感受。一条支路从节点i送功率到节点j它满足的功率损耗关系本质上是Pij² Qij² lij * ui这里lij是电流幅值的平方ui是节点i电压幅值的平方。这个式子左边是平方和右边是两个非负变量的乘积。在二维平面里等号对应的是一条抛物面边界一旦优化过程中变量范围扩大满足等式的点集是条曲线不是凸集所以没法直接套凸优化理论。如果把这个等式“松一松”把它换成不等式Pij² Qij² ≤ lij * ui画出来就是一簇二阶锥也就成了凸集。这就是二阶锥松弛这个名字的由来。听起来好像只是把一个等号改成了小于等于但数学性质完全变了从非凸问题变成凸问题理论上可以保证收敛到全局最优Matlab里的商业求解器和开源工具也都有成熟接口。1.3 为什么在配电网里用二阶锥而不是更高阶的方法有人会问既然最优潮流非凸为什么不直接用半定松弛电网文献里半定松弛SDP确实很火输电网OPF经常用。但SDP要处理矩阵变量对配电网这种辐射状结构来说往往“杀鸡用牛刀”求解规模一大内存就吃不消。配电网和输电网一个最明显的区别是结构配电网绝大多数按辐射状单电源树状设计每条负荷支路只可能有唯一的功率供给路径。Baran-Wu提出的DistFlow方程把这种径向结构利用得非常充分方程里不出现相角变量只保留电压幅值平方、支路电流平方和功率量最后导出的松弛项正好是二阶锥。这种松弛对辐射状配电网足够精确松出来的解在绝大多数正常工况下能够回到原潮流等式的可行点上。计算代价比SDP小一个量级实用价值非常高这就是配电领域普遍采用SOCP的原因。2. 二阶锥松弛模型从DistFlow到可解凸问题2.1 DistFlow方程需要哪几类核心变量标准的DistFlow模型相对简单。假设配电网是辐射状的支路集合记为E支路e的起点和终点分别是i、j那么需要关心的变量大概分这么几类变量类型物理含义说明Pij支路有功功率从i流向j的有功方向按网络拓扑确定Qij支路无功功率从i流向j的无功l_ij支路电流幅值的平方相当于Iij²不是电流本身u_i节点电压幅值的平方相当于Vi²通常在0.95²到1.05²之间变化p_j / q_j节点净注入功率分布式电源出力减去负荷负荷方向为负在Matlab建模时为了数值稳定我建议全部采用标幺值。不管原始单位是千瓦、千乏还是千伏统一除以基准功率和基准电压平方后变量数量级都在0.01到1之间求解器数值条件会好非常多。有了这些变量DistFlow方程可以写成下面三组功率平衡方程∑(P_ij - r_ij·l_ij) - ∑P_jk p_j 0∑(Q_ij - x_ij·l_ij) - ∑Q_jk q_j 0电压降落方程u_j u_i - 2(r_ij·P_ij x_ij·Q_ij) (r_ij² x_ij²)·l_ij再加上每条支路的定义l_ij·u_i P_ij² Q_ij²前两组方程在对树状网络做变量累加时是线性的电压降落方程也是线性的真正让模型变难的是最后一组方程因为它把两个变量乘到了一起展开就是P、Q的二次函数等于u和l的乘积。这正好是我前面说的非凸来源。2.2 松弛那一步具体做了什么如果直接把l_ij·u_i P_ij² Q_ij²写成等式丢给求解器那模型还是非凸的。二阶锥松弛的思想非常直接把等式放宽为不等式l_ij·u_i ≥ P_ij² Q_ij²变成这个不等式之后通过简单的代数变换可以写成二阶锥的标准形式‖ [2P_ij; 2Q_ij; u_i - l_ij] ‖ ≤ u_i l_ij也就是把原来那个乘积等式等价表达成三维向量范数小于等于一个线性表达式。这个标准形式在YALMIP里可以直接用cone函数表达在CVX里也有现成写法不用手推复杂的二阶锥参数。可能有人担心放宽之后会不会把“正确解”放丢了工程实践中只要目标函数是网损最小或者有功购电成本最小这一类对l有单调递增惩罚的目标辐射状配电网上通常可以证明最优解会自然
配电网最优潮流的二阶锥松弛建模与Matlab实现
配电网最优潮流OPF最近在配电自动化群里被聊得非常多。手里有分布式光伏、储能或者充电桩别管是做并网方案还是日内调度最后几乎都要落到同一个数学问题在网络允许的运行范围内找一组设备出力让网损、电压偏差或者购电成本尽可能小。回答这个问题的工具就是最优潮流计算。不过我接触过不少同行以为把Matlab里的fmincon或者粒子群拿过来把交流潮流约束当等式丢进去就能算。真算起来就发现问题了——模型非线性、非凸初值一换结果跟着变算半天还不一定找到全局最优。后来我在项目里换成了二阶锥松弛配电网模型直接写成凸优化Matlab调用YALMIP加免费求解器就能稳定求解计算规模大了也不慌最优解还能给出运行层面的明确物理含义。这篇就把整个思路、建模和Matlab实现完整写一遍重点讲清楚松弛在做什么、代码怎么落、坑在哪里。1. 配电网最优潮流到底难在哪1.1 算的不是传统潮流而是在网络里找最优运行点传统潮流计算是“给定负荷和发电求全网电压、功率分布”本质是解一组方程。但最优潮流不一样它多了一个“优化”的动作哪些分布式电源出力可调、无功设备补偿多少、能不能通过改变根节点电压改善末端电压、要不要调整联络开关开合状态这些都可以作为决策变量参与求解。所以它的数学结构是“在等式和不等式约束下求目标函数最小值”。这里面的等式约束不是别的正是潮流方程本身。很多人第一步就容易忽略潮流方程本来是非线性的直接嵌进优化问题之后整个可行域就变成了弯弯曲曲的非凸集合。非凸带来的后果很直接——你没法保证找到的解是全局最优可能一个看似合理的初值就把你带进局部极值里出不来。传统输电网做最优潮流早期常用内点法硬解非线性规划对初值、矩阵条件和计算规模都有讲究。配电网的网络规模虽然不如输电网大但分支多、R/X大、三相不平衡问题也常见反而更需要一套理论上更干净的方法。二阶锥松弛就是在“精确描述物理规律”和“保持凸优化可解性”之间找一个工程上最能接受的平衡点。1.2 很多人没想明白的一个点为什么潮流等式会让问题变凸性困难我先给个直观感受。一条支路从节点i送功率到节点j它满足的功率损耗关系本质上是Pij² Qij² lij * ui这里lij是电流幅值的平方ui是节点i电压幅值的平方。这个式子左边是平方和右边是两个非负变量的乘积。在二维平面里等号对应的是一条抛物面边界一旦优化过程中变量范围扩大满足等式的点集是条曲线不是凸集所以没法直接套凸优化理论。如果把这个等式“松一松”把它换成不等式Pij² Qij² ≤ lij * ui画出来就是一簇二阶锥也就成了凸集。这就是二阶锥松弛这个名字的由来。听起来好像只是把一个等号改成了小于等于但数学性质完全变了从非凸问题变成凸问题理论上可以保证收敛到全局最优Matlab里的商业求解器和开源工具也都有成熟接口。1.3 为什么在配电网里用二阶锥而不是更高阶的方法有人会问既然最优潮流非凸为什么不直接用半定松弛电网文献里半定松弛SDP确实很火输电网OPF经常用。但SDP要处理矩阵变量对配电网这种辐射状结构来说往往“杀鸡用牛刀”求解规模一大内存就吃不消。配电网和输电网一个最明显的区别是结构配电网绝大多数按辐射状单电源树状设计每条负荷支路只可能有唯一的功率供给路径。Baran-Wu提出的DistFlow方程把这种径向结构利用得非常充分方程里不出现相角变量只保留电压幅值平方、支路电流平方和功率量最后导出的松弛项正好是二阶锥。这种松弛对辐射状配电网足够精确松出来的解在绝大多数正常工况下能够回到原潮流等式的可行点上。计算代价比SDP小一个量级实用价值非常高这就是配电领域普遍采用SOCP的原因。2. 二阶锥松弛模型从DistFlow到可解凸问题2.1 DistFlow方程需要哪几类核心变量标准的DistFlow模型相对简单。假设配电网是辐射状的支路集合记为E支路e的起点和终点分别是i、j那么需要关心的变量大概分这么几类变量类型物理含义说明Pij支路有功功率从i流向j的有功方向按网络拓扑确定Qij支路无功功率从i流向j的无功l_ij支路电流幅值的平方相当于Iij²不是电流本身u_i节点电压幅值的平方相当于Vi²通常在0.95²到1.05²之间变化p_j / q_j节点净注入功率分布式电源出力减去负荷负荷方向为负在Matlab建模时为了数值稳定我建议全部采用标幺值。不管原始单位是千瓦、千乏还是千伏统一除以基准功率和基准电压平方后变量数量级都在0.01到1之间求解器数值条件会好非常多。有了这些变量DistFlow方程可以写成下面三组功率平衡方程∑(P_ij - r_ij·l_ij) - ∑P_jk p_j 0∑(Q_ij - x_ij·l_ij) - ∑Q_jk q_j 0电压降落方程u_j u_i - 2(r_ij·P_ij x_ij·Q_ij) (r_ij² x_ij²)·l_ij再加上每条支路的定义l_ij·u_i P_ij² Q_ij²前两组方程在对树状网络做变量累加时是线性的电压降落方程也是线性的真正让模型变难的是最后一组方程因为它把两个变量乘到了一起展开就是P、Q的二次函数等于u和l的乘积。这正好是我前面说的非凸来源。2.2 松弛那一步具体做了什么如果直接把l_ij·u_i P_ij² Q_ij²写成等式丢给求解器那模型还是非凸的。二阶锥松弛的思想非常直接把等式放宽为不等式l_ij·u_i ≥ P_ij² Q_ij²变成这个不等式之后通过简单的代数变换可以写成二阶锥的标准形式‖ [2P_ij; 2Q_ij; u_i - l_ij] ‖ ≤ u_i l_ij也就是把原来那个乘积等式等价表达成三维向量范数小于等于一个线性表达式。这个标准形式在YALMIP里可以直接用cone函数表达在CVX里也有现成写法不用手推复杂的二阶锥参数。可能有人担心放宽之后会不会把“正确解”放丢了工程实践中只要目标函数是网损最小或者有功购电成本最小这一类对l有单调递增惩罚的目标辐射状配电网上通常可以证明最优解会自然