
简介基于二阶锥规划的IEEE33节点主动配电网最优潮流求解程序面向电力系统优化方向的科研人员与电气工程研究生可用于多源协同运行、配电网经济运行等场景的算例复现。程序基于MATLAB平台调用YALMIP和CPLEX求解器实现二阶锥规划针对IEEE33节点系统开展24小时多时段优化模型中综合了风电、CB、SVG、OLTC、ESS等常见配电网可控资源。压缩包内含2个文件以.m主程序为核心另附说明txt文档整体约4KB代码小巧但结构清晰注释被作者打磨得较为详尽具备骨灰级注释特点适合初学者对照算法流程逐步理解。目前已有37人浏览学习可作为课题研究或课程设计的直接参考。资源中额外提供了参考文献线索可帮助读者从方法论层面把握优化建模与求解思路。1. 为什么说二阶锥规划是主动配电网最优潮流的解药主动配电网里高比例分布式电源接入后光伏反送、电压越限、潮流倒送经常同时出现。调度要回答的不再是“算个潮流”而是“怎么调储能和逆变器让整个网在满足安全约束的前提下跑得最经济”这就是主动配电网最优潮流问题。问题的核心障碍在于传统交流最优潮流是非凸优化直接求解容易陷在局部最优里算出来的策略现场不敢用。二阶锥规划通过把非凸的潮流方程做凸松弛把最优潮流转成一个可全局求解的凸问题是目前配电网方向可落地性最强的一条路径。这篇笔记适合做配电网运行策略、规划研究和调度算法的工程师从模型推导讲到代码实现把能跑的代码和踩过的坑一起讲清楚。2. 从非凸到凸二阶锥规划求解最优潮流的模型演变2.1 交流最优潮流为什么让传统算法翻车在配电网里做优化最直接的想法是拿现成的交流最优潮流模型去套。套过就会知道那个模型是真的难解。节点注入功率方程长这样P_i V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij)电压幅值、相角、电导电纳全部耦合在一起还有三角函数项。这个方程本身不是凸函数由它构成的可行域也不是凸集。凸规划里“局部最优等于全局最优”的性质在这里不成立牛顿法、内点法从不同初值出发可能收敛到完全不同的结果初值给不好甚至直接发散。对于配电网这种动辄几百上千个节点的网络靠多初值试探来找全局解计算量完全失控。有人会退一步做直流潮流最优化把电压幅值全设成1忽略无功功率和网损。这在输电网里勉强能用因为输电网支路阻抗以电抗为主电压幅值接近1有功和无功近似解耦。配电网完全不同支路电阻和电抗在同一数量级R/X常常在0.5到2之间甚至更高。光伏和负荷波动带来的无功潮流对节点电压影响非常显著直流潮流模型算出来的电压分布和网损跟实际差出一大截拿去调度是危险的。2.2 DistFlow支路方程把潮流改写成“可凸”的样子放射状配电网有一条更顺手的建模路径DistFlow支路潮流模型也叫Branch Flow模型。这个模型不直接使用节点电压相量而是描述每条支路上的有功功率、无功功率、电流平方值和节点电压平方值之间的关系。对一个放射状配网每条支路(i,j)上有四个核心变量P_ij代表从节点i流向节点j的有功Q_ij代表对应的无功l_ij代表支路电流幅值的平方v_i代表节点i电压幅值的平方。DistFlow模型用四条方程描述稳态运行状态。第一条是节点有功平衡流入节点j的支路有功减去线路损耗再减去流向子支路的有功等于节点j的本地注入功率。第二条是无功平衡形式完全一致。第三条是电压降落方程末端节点电压平方等于首端电压平方减去支路电阻和电抗上的压降再加上电流引起的修正项。第四条是支路视在功率的等式定义P_ij的平方加Q_ij的平方等于首端电压平方乘以电流平方。前三条方程都是线性的只有第四条是二次等式非凸性全部集中在这条约束里。2.3 二阶锥松弛把等式放宽成锥约束二阶锥规划的关键操作就在这里把最后那条等式约束放宽成不等式。原本要求P_ij²加Q_ij²严格等于v_i乘以l_ij现在允许它小于等于。这样做物理含义很直观支路视在功率的平方不超过首端电压平方和电流平方的乘积相当于允许模型暂时低估线路电流。这个放宽后的约束可以直接改写成标准二阶锥形式‖ [2P_ij; 2Q_ij; v_i - l_ij] ‖₂ ≤ v_i l_ij这是一个旋转二阶锥约束可以用现成的凸优化求解器高效处理。但松弛不能白做要回答一个问题松弛之后的最优解还能不能映射回物理上真实可行的潮流答案是在放射状配电网且目标函数满足一定单调性条件时最优解会自动落到锥边界上也就是不等式取到等号松弛是紧的。实际操作中大多数以最小化网损、最小化购电成本为目标的最优潮流问题都能满足这个条件。遇到松弛不紧的情况通常需要检查目标函数或者是否加了一些把解往锥内部推的约束。下面这张表可以直观看到两个模型的差异对比项传统交流最优潮流二阶锥最优潮流优化变量V、θ、P、QvV²、P、Q、lI²潮流方程非线性三角函数方程线性方程二阶锥约束求解性质非凸局部最优凸全局最优适用网架任意拓扑放射状配电网结果校验无系统化校验手段可通过锥间隙定量检查松弛带来一个额外的好处求解器能给出全局最优解的下界这个下界对算法验证极其重要。用潮流计算器去核对SOCP给出的解时如果偏差大至少知道是模型本身的问题还是数值问题而不是在非凸空间里瞎猜。3. 用PythonCVXPY从零搭一个二阶锥最优潮流3.1 网络数据准备3节点放射网的最小案例理论讲完必须落到代码。我习惯用最小可运行案例验证建模思路再往真实网络扩展。这里用一个3节点放射状配电网节点0是变电站出口母线视为平衡节点电压固定在1.0标幺节点1和节点2带负荷节点1还带了一组分布式光伏。网络数据全部采用标幺值基准容量取1 MVA基准电压取10 kV。支路参数放在字典里每个支路包含首端节点、末端节点、电阻和电抗。负荷和光伏出力也存成字典键是节点编号值是有功和无功的标幺值。注意光伏出力在节点1是0.25负荷是0.50净注入是-0.25意味着这个节点仍然从电网取电但比原来少了。数据准备的关键是保持支路方向一致首端必须是靠近变电站的一端末端是远离变电站的一端。这样后面构建节点功率平衡时才不会把方向搞反。真实网络的配电网数据通常来自GIS导出的拓扑表第一步就要做方向规范化这是个容易踩坑的点。3.2 变量、约束、目标CVXPY建模三要素变量定义分两组。支路变量有三个字典P、Q、I键是支路编号值分别是支路首端有功、首端无功和电流幅值平方。节点变量有一个字典V键是节点编号值是电压幅值平方。为什么要用电压平方而不是电压幅值因为DistFlow方程本身就是关于v的线性方程用平方形式避免了开方操作保持了问题的凸性。约束构造围绕三个层次展开。第一层是变电站节点电压固定V[0]等于1.0。第二层是支路方程包括电压降落方程和二阶锥约束。第三层是节点功率平衡以本地净注入的形式写进等式约束。这里要特别注意支路变量P和Q不要设置非负约束分布式电源接入后配电网支路出现倒送功率是正常工况把变量限制为非负会把可行解直接截掉。目标函数选的是最小化网络总损耗。在DistFlow模型里每条支路的损耗就是电阻乘以电流平方对l_ij变量而言是线性项。目标函数线性约束条件是线性和二阶锥整个问题就是标准的SOCP。3.3 完整可运行代码3节点主动配电网最优潮流import cvxpy as cp import numpy as np # ---------- 网络数据标幺值基准容量1 MVA基准电压10 kV ---------- # 支路表: id - (首端节点, 末端节点, 电阻, 电抗) lines { 0: (0, 1, 0.010, 0.030), 1: (0, 2, 0.015, 0.040), } # 负荷: 节点 - (有功, 无功) loads { 1: (0.50, 0.20), 2: (0.30, 0.10), } # 光伏出力: 节点 - (有功, 无功) pv { 1: (0.25, 0.00), 2: (0.10, 0.00), } buses [0, 1, 2] # ---------- 优化变量 ---------- P {i: cp.Variable() for i in lines} # 支路首端有功 Q {i: cp.Variable() for i in lines} # 支路首端无功 I {i: cp.Variable(nonnegTrue) for i in lines} # 支路电流平方 V {b: cp.Variable() for b in buses} # 节点电压平方 cons [] # 变电站出口电压固定为1.0 cons.append(V[0] 1.0) # ---------- 支路方程 ---------- for i, (f, t, r, x) in lines.items(): # 电压降落方程 cons.append(V[t] V[f] - 2 * (r * P[i] x * Q[i]) (r * r x * x) * I[i]) # 二阶锥松弛: || [2P, 2Q, Vf-I] || Vf I cons.append(cp.norm(cp.hstack([2 * P[i], 2 * Q[i], V[f] - I[i]])) V[f] I[i]) # ---------- 节点功率平衡 ---------- parent_of {t: i for i, (f, t, r, x) in lines.items()} children_of {b: [] for b in buses} for i, (f, t, r, x) in lines.items(): children_of[f].append(i) for b in buses: if b 0: continue # 变电站节点不建平衡方程功率由目标函数决定 p_net 0.0 q_net 0.0 # 从父支路流入本节点的功率要扣除线路损耗 if b in parent_of: i parent_of[b] f, t, r, x lines[i] p_net P[i] - r * I[i] q_net Q[i] - x * I[i] # 流出到子支路的功率 for i in children_of[b]: f, t, r, x lines[i] p_net - P[i] q_net - Q[i] # 本地净注入 光伏 - 负荷 p_gen pv.get(b, (0, 0))[0] - loads.get(b, (0, 0))[0] q_gen pv.get(b, (0, 0))[1] - loads.get(b, (0, 0))[1] cons.append(p_net p_gen) cons.append(q_net q_gen) # ---------- 节点电压运行约束0.95 ~ 1.05 pu ---------- for b in buses: if b 0: continue cons.append(V[b] 0.95 ** 2) cons.append(V[b] 1.05 ** 2) # ---------- 目标最小化网络损耗 ---------- loss sum(r * I[i] for i, (f, t, r, x) in lines.items()) prob cp.Problem(cp.Minimize(loss), cons) prob.solve(solvercp.CLARABEL, verboseFalse) # ---------- 输出 ---------- print(f求解状态: {prob.status}) print(f最优网络损耗: {loss.value:.6f} pu ({loss.value * 1.0:.3f} MW)) print(f变电站出口有功: {sum(P[i].value for i in lines):.4f} pu) for b in buses: print(f节点 {b} 电压: {np.sqrt(V[b].value):.4f} pu) for i, (f, t, r, x) in lines.items(): print(f支路 {i}: P{P[i].value:.4f} pu, Q{Q[i].value:.4f} pu, fI{I[i].value:.4f} pu)这段代码有三个关键点要说明。第一V[f]在二阶锥约束里是支路首端电压平方不是末端方向写反会得到完全错误的结果。第二节点功率平衡里流入项扣掉了损耗r乘I这是DistFlow的标准写法漏掉损耗会让功率不守恒。第三光伏和负荷的符号约定必须统一代码里采用“注入为正”所有输出结果也按这个方向解读。运行这个3节点例子最优网损通常在0.002 pu量级也就是几十kW的水平节点电压会落在0.98到1.02之间。光伏出力大于本地负荷时支路有功出现负值代表功率倒送这是主动配电网的正常工况。3.4 参数怎么调从3节点扩展到真实配电网真实网络至少几百个节点数据结构要改成按列表组织三个基础参数要确认电压基准值、容量基准值和支路方向。我一般把电压基准设为配电网额定电压容量基准设为10 MVA或者1 MVA具体看量级。支路方向处理上配电网拓扑通常从变电站开始做广度优先搜索确定每个节点的父节点再统一生成lines表。求解器的选择也有讲究。CVXPY默认带的CLARABEL对SOCP问题非常稳小规模秒出结果。如果网络规模大ECOS和MOSEK都是好选择其中MOSEK对数值病态问题的鲁棒性最好工程上大量案例用MOSEK跑数万节点的配电网SOCP不会出问题。预算有限就选CLARABEL绝大多数场景足够了。4. 主动配电网的设备建模光伏、储能和可调无功怎么塞进SOCP4.1 光伏逆变器有功出力约束和无功补偿主动配电网和传统配电网最大的区别是分布式电源具备主动调节能力。光伏逆变器的有功出力不是固定值而是可以在一定范围内调节通常允许弃光。建模时节点注入的有功需要改写成一个变量而不是常量# 光伏有功出力变量: 0 到 1.2 倍额定装机之间可调 p_pv_min 0.0 p_pv_max 1.2 # 单位为pu假设装机1.2MVA P_pv cp.Variable(nonnegTrue) cons.append(P_pv p_pv_max) cons.append(P_pv p_pv_min) # 无功出力受逆变器视在容量约束 Q_pv cp.Variable() cons.append(cp.norm(cp.hstack([P_pv, Q_pv])) 1.2)无功容量约束实际上是一个圆盘约束用cp.norm写成二阶锥形式还能融入SOCP框架。这里要注意逆变器容量约束通常写成P_pv²加Q_pv²不超过S²这是凸约束可以直接加不需要额外近似。光伏出力的上限还要乘一个气象系数或者预测曲线真实场景里p_pv_max是一个时段序列不是常量。在静态单时段问题里先按当前时刻光照强度折算一个上限值就能满足大部分分析需要。4.2 储能系统跨时段耦合的SOC约束储能是主动配电网最有价值的调节资源但建模会引入跨时段耦合单时段最优潮流直接解不了。原因是储能电池的荷电状态SOC满足递推关系下一时刻的SOC等于当前SOC加上充电功率减去放电功率再乘效率系数。这个递推关系让每个时段的优化不能独立求解必须把多个时段放在同一个问题里联立求解这就成了动态最优潮流。在SOCP框架里储能模型也不复杂关键是把充放电功率拆成充电和放电两个非负变量并且加上两者不同时为正的约束。这个约束虽然是非线性的但因为目标函数通常会让储能尽量有效率地运行实践中可以省略互补约束只靠目标函数里的充放电成本项来避免同时充放电的荒谬结果。储能约束的SOCP写法# 储能变量 SOC cp.Variable(T 1, nonnegTrue) # T个时段多一个是初始时刻 P_ch cp.Variable(T, nonnegTrue) # 充电功率 P_dis cp.Variable(T, nonnegTrue) # 放电功率 Q_ess cp.Variable(T) # 无功出力 # SOC递推 cons.append(SOC[0] 0.5) # 初始SOC 50% for t in range(T): cons.append(SOC[t1] SOC[t] P_ch[t] * eta_ch * dt / E_cap - P_dis[t] / eta_dis * dt / E_cap) cons.append(SOC[t1] 0.1) cons.append(SOC[t1] 0.9) cons.append(P_ch[t] P_ess_max) cons.append(P_dis[t] P_ess_max) cons.append(cp.norm(cp.hstack([P_dis[t] - P_ch[t], Q_ess[t]])) S_ess_max)SOC递推公式里的dt和E_cap要换算成统一的标幺值体系不然数值会直接崩掉。举个例子假设储能容量2 MWh基准功率1 MW时间步长1小时dt除以E_cap就是1除以2等于0.5SOC方程的系数就是这个量级。如果基准值选得不合适SOC的数值在0到1之间变化但充放电功率的数值可能是几十两个数量级差太远求解器数值稳定性会很差。4.3 可调无功设备并联电容器和静止无功补偿器配电网里传统的无功调节设备主要有并联电容器组和有载调压变压器。电容器组是离散设备投切状态是整数变量放进SOCP会让问题变成混合整数二阶锥规划MISOCP。这倒不是不能解只是计算代价明显上升特别是网络规模大时求解时间会增加一到两个数量级。工程上的替代做法是先把电容器组当成连续变量求解得到最优解后再用启发式规则把连续值映射回最近的离散档位然后固定档位重新求解一次。这个两阶段法的误差通常在一个档位的容量范围内对工程调度来说够用。我一般会在代码里加一个后处理函数把连续的Q_cap值圆整到最近的投切档位。静止无功补偿器SVG是连续调节设备建模就简单得多直接作为一个无功注入变量加上容量上下限约束就好。一个配电网里SVG的调节速度和精度都比电容器组高得多适合应对光伏波动带来的快速电压变化在做动态最优潮流时优先用SVG做无功支撑。4.4 目标函数加惩罚项怎么调权重不翻车设备接入后目标函数往往不能只有网损。比如要最小化购电成本、弃光惩罚、储能老化成本等。多目标加权求和时权重的标定直接影响解的质量。网损以pu为单位通常数值很小0.001量级弃光惩罚如果直接写上百倍的系数等于强制保弃光为零这个方向不一定是对的。我常用的做法是先跑一次纯网损最优得到网损量级再跑一次纯弃光最优得到弃光量级然后用两个量级之比设定初始权重再按现场电价反推修正。这样权重至少在一个合理的数量级内不会被某个目标单方面支配。权重系数确定后目标函数仍然是线性的或者凸二次型问题性质不变求解器不需要换。5. 二阶锥规划求解最优潮流的五个高频坑现象、原因、解法5.1 求解器报infeasible标幺值系统不统一现象代码看起来没问题约束都对但求解器直接返回不可行甚至报“numerical issues”。这是最常见的翻车现场。原因网络数据里有的量是标幺值有的量是有名值混在一起用。比如支路电阻给了欧姆值电压基准却用的10 kV没有除以基准阻抗。在DistFlow模型里电压平方、功率、阻抗必须全部在同一个标幺值体系下混用会导致约束条件数量级混乱可行域被压成一个空集。解决写一个统一的数据转换函数把所有有名值转成标幺值后再进入优化模型。基准阻抗等于基准电压平方除以基准容量这个值一算出来整个网络参数就统一了。调试时打印每条支路的r和x确认范围在0.0001到0.1之间超过这个范围基本就是标幺值出问题了。5.2 二阶锥约束写错方向不等号方向颠倒现象求解器能算出结果但结果明显不合理比如电压接近下限时网损反而为零或者支路功率出现荒谬的数值。原因旋转锥约束有两种常见写法一个是‖[2P; 2Q; v_i-l_ij]‖≤v_il_ij另一个等价写法是P²加Q²小于等于v_i乘l_ij。有人会习惯性写成大于等于这个方向一错问题就从凸优化变成了非凸优化求解器给的所谓最优解毫无物理意义。解决写完之后做一次独立校验取任意一条支路把P、Q、v_i、l_ij的最优值代回到P²加Q²和v_i乘l_ij里打印两者的差值。正常松弛紧的情况下两者相差不超过1e-6如果差得很大或者方向反了立刻检查约束代码。5.3 结果和真实潮流对不上忽略松弛非紧现象SOCP最优解算出来电压分布和网损都不错但拿这个解去跑交流潮流计算器节点电压和支路功率全对不上误差超过5%。原因目标函数或者约束条件把最优解推到了锥的内部而不是边界。典型场景是最小化购电成本时如果节点电价设定异常求解器倾向于让某些支路电流虚高这时锥松弛不是紧的SOCP解在物理上是不可实现的。解决计算锥间隙。对每条支路计算(P_ij²加Q_ij²)除以(v_i乘l_ij)再减去1取所有支路的最大值。如果这个值大于1e-4说明松弛损失了精度。针对这种情况常用的修法是在目标函数里加一个很小的电流惩罚项比如0.001乘以所有支路的l_ij之和这样会鼓励求解器把电流压低让解回到锥边界上。5.4 储能SOC和功率数量级差距过大导致求解器震荡现象加了储能模型后求解器迭代次数暴涨有时几百次不收敛有时收敛了SOC曲线出现锯齿状跳变。原因储能容量基准和功率基准不一致。比如基准功率是1 MW储能容量是5 MWh时间步长是15分钟那么dt除以E_cap等于0.25除以5等于0.05而充放电功率的标幺值可能在0.8左右SOC数值在0到1之间三个量级完全不同导致线性方程组的条件数变得很差。解决统一容量和功率的基准。把储能容量也折算成标幺值或者直接把SOC递推方程的系数整体放大100倍让约束矩阵的条件数维持在合理范围。另一个更简单的做法是用秒作为时间单位确保dt和E_cap的量级差不超过100倍。5.5 电容器组离散档位处理不当现象把电容器组当成连续变量求解结果出来了但圆整到离散档位后重新跑潮流电压越限。原因圆整操作改变了无功注入量直接影响节点电压。如果原来的最优解里某个电容器的无功出力刚好在档位边界附近圆整到上一档或下一档都可能让电压越过上下限。解决不要直接圆整到最近档位而是优先选择能让电压留在安全区间的档位。具体做法是先固定其他设备把电容器组的每个离散档位依次代入原问题求解一次选目标函数最优且所有约束满足的档位。网络规模大时可以用灵敏度分析缩小候选档位集合一般不超过5个档位。这个方法是增加计算量的一般固定档位后再求解一次SOCP即可。6. 从单时段到动态最优潮流验证方法和一个进阶技巧单时段SOCP最优潮流跑通后下一步自然就是动态最优潮流把一天24个或者96个时段放进同一个问题里联立求解。储能SOC的递推约束把各时段耦合在一起模型不再是一个SOCP在独立运算而是一整个时间块的SOCP变量数量直接乘以时段数。求解器的性能在这里拉开差距同一批数据用CLARABEL可能几分钟换MOSEK只要几十秒体感差别很明显。我验证一个SOCP最优潮流结果的习惯动作有三步。第一步看锥间隙逐条支路检查松弛紧度数值大于1e-4就要加惩罚项修正。第二步做交流潮流回代把SOCP给出的节点有功无功注入值固定住用传统潮流计算器跑一次完整ACPF对比电压幅值和支路功率误差2%以内算正常。第三步做时序一致性检查动态场景里看储能SOC曲线的连续性SOC跳变意味着约束写错了或者数值有问题。一个值得尝试的进阶技巧是把SOCP当成“热启动器”用它的解作为非线性交流最优潮流的初值。具体做法是先解SOCP得到各DG的出力和储能策略然后把这个解作为内点法的初始点去解精确的ACOPF。这样ACOPF的迭代次数会大幅减少而且大概率能收敛到和SOCP接近的局部最优解。等于用凸问题画了一个靠谱的可行域入口再用精确模型做收尾。这个方法在多个配电网测试系统上都验证过比直接跑ACOPF稳定得多也比只用SOCP结果可信得多。做主动配电网调度这一年多下来最大的体会是SOCP最优潮流的价值不在理论研究而在于它给了一套可以反复使用的计算底座。模型从3节点扩展到实际馈线坑来来去去就那么几个把标幺值、松弛紧度、求解器选型这三件事管住大部分问题都能顺下来。希望这篇笔记能帮你少走一段弯路也希望你第一次跑通代码时能体会到“看电压曲线从越限被拉回来”的那种踏实感。本文还有配套的精品资源点击获取