简介面向电气工程毕业设计及配电网优化研究者这份程序与论文配套资料聚焦主动配电网最优潮流与多源协调运行优化问题涵盖基于二阶锥规划的IEEE33节点潮流求解及综合能源低碳运行策略相关参考。资源提供MATLAB求解程序考虑风电、无功补偿、有载调压变压器、储能等多设备24小时多时段协调调度使用CPLEX商业求解器处理混合整数线性规划模型同时附带运行日志、配电网结构图、潮流计算说明PPT及可参考论文清单便于对照复现与结果分析。压缩包共25个文件以log运行日志、m源文件、txt说明、png结构图、pptx演示文稿为主整体仅507KB轻量易用。目前已有106人学习适合需要快速上手配电网最优潮流二阶锥规划建模的本科生或研究生。1. 二阶锥规划正在成为主动配电网最优潮流的主流通用解法主动配电网和传统配电网最大的区别在于分布式电源、储能和柔性负荷大规模接入之后电压、功率和网损之间的耦合不再能靠简化的潮流公式近似处理。最优潮流问题直接落到非线性非凸优化上直接求解既慢、又很难保证找到全局最优解。二阶锥规划SOCP通过变量替换和松弛把非凸潮流等式放进一个凸锥约束里使得问题能够被现代内点法稳定求解这也是《基于二阶锥规划的主动配电网最优潮流求解》和《主动配电网多源协同运行优化研究》两篇论文共用的核心方法。随程序附带的论文在知网按篇名检索即可找到算例和公式都对得上适合正在做配电网研究的学生、做电网规划的工程师以及想给调度系统加算法模块的底层开发人员。读这套程序时最有价值的不是抄一个模型而是理解松弛边界在什么时候是紧的、什么时候会失效。下面就从 DistFlow 模型开始把这一整套方案讲清楚。2. SOCP 模型为什么硬性非凸潮流会被松弛成凸锥2.1 DistFlow 方程里的非凸来源是电压与电流的乘积项配电网最优潮流问题里最常用的数学模型是 DistFlow。对一条从节点 i 指向节点 j 的支路 ij定义 P_ij 和 Q_ij 为支路有功和无功u_i 为节点 i 电压幅值的平方L_ij 为支路电流幅值的平方r_ij 和 x_ij 是支路阻抗。DistFlow 的电压方程写为$$ u_j u_i - 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}{u_i} $$第二个等式不是线性的也不是凸约束。这里的问题主要有两个一个是 u_i 在分母上属于分式非线性另一个是如果把等式两边同时乘 u_i得到 P_ij^2 Q_ij^2 u_i L_ij右边是两个变量的乘积在优化里是双线性非凸项。早期处理方式是直接把网损忽略掉变成直流潮流但配电网中线路电阻不小R/X 比较高忽略有功损耗带来的误差在电压约束和调度结果上都不可接受。把两条式子合起来看可以理解成一条支路的功率传输能力受制于首端电压平方和支路电流平方的乘积。这个耦合关系是问题非凸的根源。主动配电网里多个分布式电源接入之后功率流动方向变得双向非凸性对求解器的影响会被进一步放大。这就是为什么不能直接套用输电网最优潮流里常用的线性化或牛顿法方案。2.2 变量替换后二阶锥约束只放宽了“等式”而没有放宽功率平衡SOCP 做法的第一步是变量替换把节点电压模值的平方 u_i 和支路电流模值的平方 L_ij 直接作为优化变量。替换之后电压方程变成线性约束只剩 P_ij^2 Q_ij^2 u_i L_ij 这一个非凸等式需要处理。二阶锥松弛把等式放宽成不等式$$ \sqrt{(2P_{ij})^2 (2Q_{ij})^2 (u_i - L_{ij})^2} \le u_i L_{ij} $$这个约束等价于 u_i L_ij ≥ P_ij^2 Q_ij^2也就是把原等式放松成“电压电流平方乘积要足够大”。从物理上看它允许线路用略高于必要水平的电压或电流来传输同一段功率这是对实际物理约束的一种外逼近。在 CVXPY 里构造这个锥约束时我一般直接用 cp.SOC 函数避免手写二维范数带来维度错误。常见写法如下import cvxpy as cp # u[i] 为节点i电压平方变量L[k] 为支路k电流平方变量 # P[k]、Q[k] 为支路k首端流过功率 # 约束u[i] * L[k] P[k]^2 Q[k]^2 的SOC松弛形式 constraints.append( cp.SOC( u[from_node[k]] L[k], cp.hstack([2 * P[k], 2 * Q[k], u[from_node[k]] - L[k]]) ) )cp.SOC(t, x) 在 CVXPY 中表示 ||x||_2 t因此第一个参数是锥的半径第二个参数是锥内的向量。这里把 2P、2Q 和 u-L 组成向量半径取 uL展开之后正是前面那个二阶锥不等式。需要注意 from_node 的下标一定要取支路首端节点不要取成 j否则锥约束表达的是错误方向求解结果会完全失真。这一步的松弛并不改变节点功率平衡约束也没有改变电压方程所以最终得到的解仍然满足潮流模型里除“电流平方等于功率平方除以电压平方”之外的物理关系。判断 SOCP 求解结果能不能还原成真实潮流只需要看这个锥约束在最优解处是否取到了等号。这个检查点后面会专门讲方法。2.3 SOCP、SDP、LP 三种凸松弛的边界与成本选择 SOCP 而不是其他凸松弛是从精度和求解规模之间取的折中。把三种方案放在一起对比可以更清楚地看出这条技术路线的位置。松弛方法约束形式求解规模适用场景LP 线性化把潮流方程一次性线性化最小速度快规划近似、输电网粗算配电网误差偏大SOCP 二阶锥松弛保留电压电流二次关系只放松一个等式中等现代内点法高效主动配电网最优潮流、多时段运行优化SDP 半定松弛把整体潮流问题放入矩阵正半定锥变量多计算量大辐射状或弱环网的全局最优性验证SDP 在理论上比 SOCP 提供了一个更宽松且更接近原问题全局最优的边界但矩阵变量的维度会随着节点数迅速增长例如 IEEE 123 节点这样规模不大的配电网也可能让 SDP 求解变得非常慢。LP 线性化虽然快但在重负荷和低电压场景下容易低估网损导致调度结果在实际执行时电压越限。SOCP 在实际工程里还有一个额外优势它允许把储能、分布式电源、可调无功设备等运行约束全部保留在凸优化框架内。只要目标函数和额外约束本身是凸的整个问题依然是凸优化求解器给出的解就是全局最优解。这也是《主动配电网多源协同运行优化研究》里能够叠加多个设备约束的结构性基础。3. 用 Python CVXPY 复现《基于二阶锥规划的主动配电网最优潮流求解》的最小算例3.1 准备节点与支路数据把 DistFlow 写成可优化约束拿到这套程序时先不要急着改目标函数第一步应该是把数据接口理清楚。程序里的核心数据是 branch 数组和 load 数组branch 每一行保存支路首端节点、末端节点、电阻、电抗load 数组按节点保存有功负荷和无功负荷。下面这段代码给出了一个三节点辐射状配网的最小算例实际上换成 IEEE 33 节点数据只需要替换 branch、p_load、q_load 三个数组。import cvxpy as cp import numpy as np # 三节点示例节点0为根节点0-1-2 # branch: 首端, 末端, r(标幺), x(标幺) branch np.array([ [0, 1, 0.02, 0.01], [1, 2, 0.02, 0.01] ]) n_nodes 3 # 负荷节点1、2各 0.2MW 0.1Mvar p_load np.array([0.0, 0.2, 0.2]) q_load np.array([0.0, 0.1, 0.1]) # 分布式电源有功上限节点0不装DG pg_max np.array([0.0, 0.3, 0.3]) qg_max np.array([0.0, 0.1, 0.1]) # 优化变量定义 P cp.Variable(len(branch)) # 支路首端有功 Q cp.Variable(len(branch)) # 支路首端无功 L cp.Variable(len(branch)) # 支路电流平方 u cp.Variable(n_nodes) # 节点电压平方 pg cp.Variable(n_nodes) # DG有功注入 qg cp.Variable(n_nodes) # DG无功注入 constraints [u[0] 1.0] # 根节点电压设为1.0 p.u. (平方为1.0) for k, (i, j, r, x) in enumerate(branch): constraints.append( u[j] u[i] - 2 * (r * P[k] x * Q[k]) (r * r x * x) * L[k] ) constraints.append( cp.SOC(u[i] L[k], cp.hstack([2 * P[k], 2 * Q[k], u[i] - L[k]])) )这段代码把 DistFlow 的电压方程和二阶锥约束逐个支路添加到 constraints 列表。注意第 10 行使用了 u[j] u[i] -...它保证了线路末端电压平方等于首端电压平方减去线路压降是潮流模型的物理核心。后面的 SOC 约束用来替代电流平方的非凸等式是二阶锥规划的关键改动。节点功率平衡在这里隐含在 P、Q 与 pg、qg 的关系里需要按网络拓扑补充。对于上面这个简单串联网络约束写法为# 节点1注入 支路0末端功率 - 支路1首端功率 constraints.append(P[0] - P[1] pg[1] - p_load[1]) constraints.append(Q[0] - Q[1] qg[1] - q_load[1]) # 节点2注入 支路1末端功率 constraints.append(P[1] pg[2] - p_load[2]) constraints.append(Q[1] qg[2] - q_load[2])对于实际的大型配电网这种逐节点手写方式不可取。常见做法是用关联矩阵 incidence matrix 把支路功率和节点注入功率一次性表达左乘拓扑矩阵得到每个节点的流入减流出再与注入等式匹配。程序里如果已经有邻接表或拓扑矩阵建议直接用矩阵乘法写平衡约束避免在 33 节点里手写几十条等式。3.2 网损最小目标下 SOCP 的松弛紧性最优潮流里可以用不同目标函数。论文中常见的网损最小目标写为$$ \min \sum_{k}(r_k L_k) $$把线路电阻和电流平方变量相乘再求和是网损的标准表达。这个目标对 SOCP 松弛有很好的紧性系统为了降低网损会尽量降低支路电流平方同时又会尽量把电压水平稳定在合理区间所以最优解通常会逼到锥边界上。objective cp.Minimize(cp.sum(cp.multiply(branch[:, 2], L))) prob cp.Problem(objective, constraints) prob.solve(solvercp.ECOS) print(最优网损:, prob.value) print(节点电压平方:, u.value) print(DG出力:, pg.value)代码中的branch[:, 2]取的是 branch 数组第三列也就是每个支路的电阻 r。cp.multiply是 CVXPY 里元素级乘法不能直接用 NumPy 的*操作 Variable 对象。prob.value返回最优目标值u.value是电压平方结果实际电压幅值需要取平方根。求解完成后要立刻检查锥约束的紧性。一个简单做法是计算每条支路的 u_i * L_k - (P_k^2 Q_k^2)。如果这个值对全部支路都近似为 0说明二阶锥松弛是紧的SOCP 解就是原最优潮流问题的可行解如果某条支路的差值明显大于 0说明那条支路的电流或电压被人为放大直接拿去潮流计算会对不上。3.3 提取结果并用潮流校验判断解是否物理可行求解器给出的解需要经过潮流回代验证这一步不能省。原因是 SOCP 只保证在松弛空间里最优不保证在原来的物理可行域里严格成立。校验时我习惯把求解得到的 pg、qg 和根节点电压作为给定值重新用常规潮流计算各节点电压再和 SOCP 结果的 u.value 对比。# 简易校验把 SOCP 解中的注入功率代入 DistFlow 递推 v_sq np.ones(n_nodes) for k, (i, j, r, x) in enumerate(branch): p_ij P.value[k] q_ij Q.value[k] v_sq[j] v_sq[i] - 2 * (r * p_ij x * q_ij) \ (r * r x * x) * L.value[k] print(SOCP u:, u.value) print(回代 u:, v_sq) print(最大偏差:, np.max(np.abs(u.value - v_sq)))这段代码利用 DistFlow 的电压方程从根节点往下递推。如果最大偏差在 1e-4 以下说明 SOCP 解满足电压方程如果偏差达到 1e-2 甚至更大问题通常出在锥约束不紧或者是节点功率平衡约束写错。回代验证不能直接判断锥约束紧性但能快速抓住“电压方程写错”这类低级问题建议放在求解后的第一轮检查里。4. 把多源协同运行优化转到同一套二阶锥框架里4.1 多源协同的调度变量和约束拆分《主动配电网多源协同运行优化研究》不是只解决单时刻最优潮流而是把一组可控设备放到同一个优化框架里协调。多源协同的变量可以分为四类连续有功变量、连续无功变量、储能状态变量、离散投切变量。常见做法是把它们全部写成约束再和上一章的 SOCP 潮流约束共同求解。有功变量包括分布式光伏/风机的有功出力、储能充电功率、储能放电功率、上级电网购电功率。无功变量包括光伏逆变器无功、储能逆变器无功、并联电容器组输出和动态无功补偿装置的无功。储能状态变量主要是 SOC按时间递推$$ SOC_{t1} SOC_t - \frac{\eta_c P_{c,t}\Delta t}{E_{rated}} \frac{P_{d,t}\Delta t}{E_{rated}\eta_d} $$其中 eta_c 是充电效率eta_d 是放电效率P_c 和 P_d 是充放电功率E_rated 是储能额定容量。这个约束是线性的SOCP 框架可以直接容纳。真正带来计算负担的是“储能不能同时充电和放电”这样的互斥约束它需要引入 0-1 整数变量将问题变成混合整数二阶锥规划。配电网侧的离散设备主要包括有载调压变压器分接头和分组投切电容器。这类设备如果按连续变量处理得到的调度指令无法实际执行如果逐个离散档位枚举组合数会爆炸。更自然的做法是把每个分接头档位建模为整数变量并限制相邻时刻最大调节档数。4.2 储能、分布式电源和可调设备的典型参数在搭建多源协同优化程序之前最好把设备参数集中放在一个字典里方便换算成标幺值。下面是一组用于 10 kV 馈线的常见参数设备类型参数数值示例说明储能额定容量1 MWh决定 SOC 递推中的分母储能最大充/放电功率0.25 MW充放电功率上界储能效率0.92 / 0.92充电与放电分开填光伏逆变器额定功率0.5 MVA有功与无功功率耦合上限光伏逆变器无功能力±0.3 Mvar功率因数约 0.6 至 0.95电容器组单组无功0.1 Mvar离散投切整数倍出力上级电网购电功率上限5 MVA根节点注入受限这些参数在程序里通常已经做成可配置对象改参数时要注意单位一致。配电网程序最容易踩坑的地方是阻抗用欧姆、功率用千瓦、电压用千伏混在一起。我一般会把功率统一为 MW阻抗统一为标幺值根节点电压统一为 1.0 p.u.。4.3 从单时段 OPF 扩展到 24 小时多时段模型单时段最优潮流扩展到多时段运行优化时只需要给原本的节点变量、DG 变量和储能变量加一个时间维度。目标函数则改为一天内的购电费用、DG 运行成本和储能退化成本之和。下面给出一个多时段 SOCP 约束的构造片段其中只保留了储能部分的关键约束。import cvxpy as cp T 24 # 24小时 E_rated 1.0 # MWh dt 1.0 # 小时 eff_c 0.92 eff_d 0.92 soc cp.Variable(T 1) # 含初始时刻 p_ch cp.Variable(T, nonnegTrue) # 充电功率 p_dis cp.Variable(T, nonnegTrue) # 放电功率 b_ch cp.Variable(T, booleanTrue) # 1表示充电时段 b_dis cp.Variable(T, booleanTrue) # 1表示放电时段 constraints [soc[0] 0.5, soc[T] 0.5] # 初始与终态SOC各50% for t in range(T): constraints.append(soc[t1] soc[t] - p_ch[t]*dt/(E_rated*eff_c) p_dis[t]*dt/E_rated*eff_d) constraints.append(p_ch[t] 0.25 * b_ch[t]) constraints.append(p_dis[t] 0.25 * b_dis[t]) constraints.append(b_ch[t] b_dis[t] 1) constraints.append(soc[t1] 0.9) constraints.append(soc[t1] 0.1)这段代码里b_ch 和 b_dis 是布尔变量用来保证同一个时段不会同时充电和放电。注意 SOC 递推式的微分项带了效率系数在充电方向电能不能全部存入所以 SOC 增量要除以充电效率在放电方向电池放出的电量来自内部储存所以 SOC 减少量要乘以放电效率。如果把两个效率写反储能会在一天内凭空多出能量必须检查这个细节。将这段储能约束与前面章节中已经写好的 DistFlow 约束合并再给光伏、风机、逆变器无功加上上界就构成了一个完整的主动配电网多源协同运行优化模型。求解器面对的是一个混合整数二阶锥规划问题只有在求解器支持整数变量时才能直接处理这就进入了下一章要讲的求解器选择问题。5. 求解器选择与参数调优MOSEK、Gurobi、ECOS 的适用边界5.1 四种求解器的特点和适用场景同一套 CVXPY 模型交给不同求解器结果和速度差异可以非常大。纯连续变量的 SOCP 和带整数变量的 MISOCP适合的求解器并不完全相同。下表整理了四种常见选择求解器问题类型特点适用规模ECOSSOCP、LP稳妥、默认精度较好不支持整数连续变量小规模算例SCSSOCP、SDP内存占用小采用一阶方法大规模但精度要求不高的场景MOSEKSOCP、MISOCP内点法 分支定界数值最稳工商授权可用配网调度首选GurobiSOCP、MISOCP整数问题求解速度快参数丰富多时段多设备优化如果只是复现论文里的单时段最优潮流ECOS 就够了。如果要做 24 小时联合优化并且带储能、电容器组等整数变量ECOS 会直接报错因为模型里出现了布尔变量。这时可以用 MOSEK 或 Gurobi它们的 MIP 分支定界能力能处理中小规模 MISOCP。程序里如果看到 solver 参数可以配置建议先跑通 ECOS再切换到 MOSEK 对比目标值。5.2 在 CVXPY 里设置相对停止容忍度和 MIP gap多时段模型经常求解时间偏长默认参数不一定最优。以 MOSEK 为例通常需要把锥优化相对对偶间隙从默认值调小保证 SOCP 解的边界质量。# 用 MOSEK 求解纯SOCP收紧相对对偶间隙 prob.solve( solvercp.MOSEK, mosek_params{ MSK_DPAR_INTPNT_CO_TOL_REL_GAP: 1e-8, MSK_DPAR_INTPNT_CO_TOL_NEAR_REL: 1e-8 }, verboseTrue )第一个参数控制内部点法的相对对偶停止阈值越小结果越精确但耗时更长。第二个参数控制近可行解判定阈值配电网潮流约束本身数值尺度较接近设置过严反而会让求解器在后期做无意义迭代。一般先用默认值看目标值再把相对间隙从 1e-6 逐步收紧到 1e-8观察目标值是否还在变化。如果目标值稳定说明当前精度已经足够。使用 Gurobi 处理 MISOCP 时重点调两个参数MIPGap 和 TimeLimit。prob.solve(solvercp.GUROBI, MIPGap1e-4, TimeLimit300, verboseTrue)MIPGap 表示当前整数解与最优下界之间的相对间隙设置为 1e-4 意味着求解器接受 0.01% 的次优性。TimeLimit 给 300 秒上限非常常见配电网多时段模型的整数变量如果特别多不设时间限制可能跑几小时。建议在离线研究里先用 0% 的 gap 算出基准值再在工程部署中把 gap 适当放宽。5.3 常遇到的不可行、数值病态和松弛不成 tight 的处理思路不可行是配电网优化里最常见的问题。出现“infeasible”时先确认负荷和 DG 单位是否都是 MW/Mvar再检查根节点电压约束是否写成 u[0] 1.0 而不是 V[0] 1.0。还要注意节点功率平衡方向在 DistFlow 里支路功率流动方向由支路表序号决定如果拓扑数据里首末端方向不一致平衡等式就会差一个负号。数值病态通常来自量纲。线路电阻用欧姆表示时可能只有 0.5但如果把功率单位写成 kW负荷数值变成几百优化变量尺度相差三个数量级SOCP 的内点法就会在收敛精度上大幅退化。处理办法是统一转成标幺值取基准容量 10 MVA10 kV 配电网阻抗标幺乘以基准值分量负荷除以基准容量电压取 1.0 p.u.。这也是论文程序里常见的预处理方式。松弛不成 tight 的问题往往被忽视。如果目标函数不是网损也不是严格依赖电流平方的成本SOCP 松弛可能给出一个在锥内部的解。解决办法是在目标函数里加上一个小系数网损项例如 0.01 * sum(r_k * L_k)。这样既不会明显改变经济调度结果又能让锥约束在最优解处尽量贴近边界还原出物理可行的电压和功率分布。6. 验证 SOCP 松弛精度的两个实用技巧6.1 用锥边界间隙检查每条支路是否落在边界上拿到一组最优解之后最值得做的第一件事是逐条支路检查锥边界间隙。这项检查能直接说明 SOCP 解是否可信是写论文、改程序时最实用的诊断手段。# 计算每条支路的锥边界间隙 lhs u.value[branch[:, 0].astype(int)] * L.value # u_i * L_k rhs P.value**2 Q.value**2 # 归一化间隙方便观察量级 gap (lhs - rhs) / np.maximum(rhs, 1e-6) print(最大锥间隙:, np.max(gap)) print(最大间隙支路:, np.argmax(gap))如果 gap 总量接近 0说明 SOCP 松弛是紧的最优解可以直接当作真实潮流解使用如果 gap 达到 1e-3 甚至更大说明这条支路对应的锥约束没有触及边界。此时最优解的经济性值可能没问题但电压、电流值并不能真实还原到配电网场景。前文提过给目标函数增加网损项或者把该支路的首端电压固定到某一合理工作点通常可以恢复紧性。6.2 用对偶变量判断配电网瓶颈提升论文分析深度二阶锥规划自带一条其他启发式算法给不了的信息链路对偶变量。在 CVXPY 里求解完成后可以通过constraints[i].dual_value取出对应约束的影子价格。对节点功率平衡约束对偶值表示该节点单位功率增量带来的最优目标变化量也即节点的实时电价或阻塞信号对支路电压约束对偶值能定位电压越限最敏感的支路。多源协同运行优化里这部分信息可以用来解释“为什么储能要在某个时刻充电”这样的结果。比如在某个时段某节点储能放电约束的对偶值很高说明该时刻该节点附近的功率供应紧张储能的边际价值大反之如果对偶值接近 0说明储能动作对目标几乎没有影响可以优先调整其他设备。论文里如果需要做灵敏度分析和经济解释这是一个很容易出图、出表的小技巧。具体导出时注意 CVXPY 对约束取 dual_value 的顺序和 constraints 列表一致。在多时段模型中建议把每个时段的平衡约束单独保存到列表里而不是写在一个大矩阵里否则后面要对单一节点、单一时刻做分析会非常不方便。配合上一节提到的锥边界间隙检查这套流程基本覆盖了从“能出解”到“解可信”的全部关键环节。本文还有配套的精品资源点击获取
二阶锥规划在主动配电网最优潮流与多源协同中的应用
简介面向电气工程毕业设计及配电网优化研究者这份程序与论文配套资料聚焦主动配电网最优潮流与多源协调运行优化问题涵盖基于二阶锥规划的IEEE33节点潮流求解及综合能源低碳运行策略相关参考。资源提供MATLAB求解程序考虑风电、无功补偿、有载调压变压器、储能等多设备24小时多时段协调调度使用CPLEX商业求解器处理混合整数线性规划模型同时附带运行日志、配电网结构图、潮流计算说明PPT及可参考论文清单便于对照复现与结果分析。压缩包共25个文件以log运行日志、m源文件、txt说明、png结构图、pptx演示文稿为主整体仅507KB轻量易用。目前已有106人学习适合需要快速上手配电网最优潮流二阶锥规划建模的本科生或研究生。1. 二阶锥规划正在成为主动配电网最优潮流的主流通用解法主动配电网和传统配电网最大的区别在于分布式电源、储能和柔性负荷大规模接入之后电压、功率和网损之间的耦合不再能靠简化的潮流公式近似处理。最优潮流问题直接落到非线性非凸优化上直接求解既慢、又很难保证找到全局最优解。二阶锥规划SOCP通过变量替换和松弛把非凸潮流等式放进一个凸锥约束里使得问题能够被现代内点法稳定求解这也是《基于二阶锥规划的主动配电网最优潮流求解》和《主动配电网多源协同运行优化研究》两篇论文共用的核心方法。随程序附带的论文在知网按篇名检索即可找到算例和公式都对得上适合正在做配电网研究的学生、做电网规划的工程师以及想给调度系统加算法模块的底层开发人员。读这套程序时最有价值的不是抄一个模型而是理解松弛边界在什么时候是紧的、什么时候会失效。下面就从 DistFlow 模型开始把这一整套方案讲清楚。2. SOCP 模型为什么硬性非凸潮流会被松弛成凸锥2.1 DistFlow 方程里的非凸来源是电压与电流的乘积项配电网最优潮流问题里最常用的数学模型是 DistFlow。对一条从节点 i 指向节点 j 的支路 ij定义 P_ij 和 Q_ij 为支路有功和无功u_i 为节点 i 电压幅值的平方L_ij 为支路电流幅值的平方r_ij 和 x_ij 是支路阻抗。DistFlow 的电压方程写为$$ u_j u_i - 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}{u_i} $$第二个等式不是线性的也不是凸约束。这里的问题主要有两个一个是 u_i 在分母上属于分式非线性另一个是如果把等式两边同时乘 u_i得到 P_ij^2 Q_ij^2 u_i L_ij右边是两个变量的乘积在优化里是双线性非凸项。早期处理方式是直接把网损忽略掉变成直流潮流但配电网中线路电阻不小R/X 比较高忽略有功损耗带来的误差在电压约束和调度结果上都不可接受。把两条式子合起来看可以理解成一条支路的功率传输能力受制于首端电压平方和支路电流平方的乘积。这个耦合关系是问题非凸的根源。主动配电网里多个分布式电源接入之后功率流动方向变得双向非凸性对求解器的影响会被进一步放大。这就是为什么不能直接套用输电网最优潮流里常用的线性化或牛顿法方案。2.2 变量替换后二阶锥约束只放宽了“等式”而没有放宽功率平衡SOCP 做法的第一步是变量替换把节点电压模值的平方 u_i 和支路电流模值的平方 L_ij 直接作为优化变量。替换之后电压方程变成线性约束只剩 P_ij^2 Q_ij^2 u_i L_ij 这一个非凸等式需要处理。二阶锥松弛把等式放宽成不等式$$ \sqrt{(2P_{ij})^2 (2Q_{ij})^2 (u_i - L_{ij})^2} \le u_i L_{ij} $$这个约束等价于 u_i L_ij ≥ P_ij^2 Q_ij^2也就是把原等式放松成“电压电流平方乘积要足够大”。从物理上看它允许线路用略高于必要水平的电压或电流来传输同一段功率这是对实际物理约束的一种外逼近。在 CVXPY 里构造这个锥约束时我一般直接用 cp.SOC 函数避免手写二维范数带来维度错误。常见写法如下import cvxpy as cp # u[i] 为节点i电压平方变量L[k] 为支路k电流平方变量 # P[k]、Q[k] 为支路k首端流过功率 # 约束u[i] * L[k] P[k]^2 Q[k]^2 的SOC松弛形式 constraints.append( cp.SOC( u[from_node[k]] L[k], cp.hstack([2 * P[k], 2 * Q[k], u[from_node[k]] - L[k]]) ) )cp.SOC(t, x) 在 CVXPY 中表示 ||x||_2 t因此第一个参数是锥的半径第二个参数是锥内的向量。这里把 2P、2Q 和 u-L 组成向量半径取 uL展开之后正是前面那个二阶锥不等式。需要注意 from_node 的下标一定要取支路首端节点不要取成 j否则锥约束表达的是错误方向求解结果会完全失真。这一步的松弛并不改变节点功率平衡约束也没有改变电压方程所以最终得到的解仍然满足潮流模型里除“电流平方等于功率平方除以电压平方”之外的物理关系。判断 SOCP 求解结果能不能还原成真实潮流只需要看这个锥约束在最优解处是否取到了等号。这个检查点后面会专门讲方法。2.3 SOCP、SDP、LP 三种凸松弛的边界与成本选择 SOCP 而不是其他凸松弛是从精度和求解规模之间取的折中。把三种方案放在一起对比可以更清楚地看出这条技术路线的位置。松弛方法约束形式求解规模适用场景LP 线性化把潮流方程一次性线性化最小速度快规划近似、输电网粗算配电网误差偏大SOCP 二阶锥松弛保留电压电流二次关系只放松一个等式中等现代内点法高效主动配电网最优潮流、多时段运行优化SDP 半定松弛把整体潮流问题放入矩阵正半定锥变量多计算量大辐射状或弱环网的全局最优性验证SDP 在理论上比 SOCP 提供了一个更宽松且更接近原问题全局最优的边界但矩阵变量的维度会随着节点数迅速增长例如 IEEE 123 节点这样规模不大的配电网也可能让 SDP 求解变得非常慢。LP 线性化虽然快但在重负荷和低电压场景下容易低估网损导致调度结果在实际执行时电压越限。SOCP 在实际工程里还有一个额外优势它允许把储能、分布式电源、可调无功设备等运行约束全部保留在凸优化框架内。只要目标函数和额外约束本身是凸的整个问题依然是凸优化求解器给出的解就是全局最优解。这也是《主动配电网多源协同运行优化研究》里能够叠加多个设备约束的结构性基础。3. 用 Python CVXPY 复现《基于二阶锥规划的主动配电网最优潮流求解》的最小算例3.1 准备节点与支路数据把 DistFlow 写成可优化约束拿到这套程序时先不要急着改目标函数第一步应该是把数据接口理清楚。程序里的核心数据是 branch 数组和 load 数组branch 每一行保存支路首端节点、末端节点、电阻、电抗load 数组按节点保存有功负荷和无功负荷。下面这段代码给出了一个三节点辐射状配网的最小算例实际上换成 IEEE 33 节点数据只需要替换 branch、p_load、q_load 三个数组。import cvxpy as cp import numpy as np # 三节点示例节点0为根节点0-1-2 # branch: 首端, 末端, r(标幺), x(标幺) branch np.array([ [0, 1, 0.02, 0.01], [1, 2, 0.02, 0.01] ]) n_nodes 3 # 负荷节点1、2各 0.2MW 0.1Mvar p_load np.array([0.0, 0.2, 0.2]) q_load np.array([0.0, 0.1, 0.1]) # 分布式电源有功上限节点0不装DG pg_max np.array([0.0, 0.3, 0.3]) qg_max np.array([0.0, 0.1, 0.1]) # 优化变量定义 P cp.Variable(len(branch)) # 支路首端有功 Q cp.Variable(len(branch)) # 支路首端无功 L cp.Variable(len(branch)) # 支路电流平方 u cp.Variable(n_nodes) # 节点电压平方 pg cp.Variable(n_nodes) # DG有功注入 qg cp.Variable(n_nodes) # DG无功注入 constraints [u[0] 1.0] # 根节点电压设为1.0 p.u. (平方为1.0) for k, (i, j, r, x) in enumerate(branch): constraints.append( u[j] u[i] - 2 * (r * P[k] x * Q[k]) (r * r x * x) * L[k] ) constraints.append( cp.SOC(u[i] L[k], cp.hstack([2 * P[k], 2 * Q[k], u[i] - L[k]])) )这段代码把 DistFlow 的电压方程和二阶锥约束逐个支路添加到 constraints 列表。注意第 10 行使用了 u[j] u[i] -...它保证了线路末端电压平方等于首端电压平方减去线路压降是潮流模型的物理核心。后面的 SOC 约束用来替代电流平方的非凸等式是二阶锥规划的关键改动。节点功率平衡在这里隐含在 P、Q 与 pg、qg 的关系里需要按网络拓扑补充。对于上面这个简单串联网络约束写法为# 节点1注入 支路0末端功率 - 支路1首端功率 constraints.append(P[0] - P[1] pg[1] - p_load[1]) constraints.append(Q[0] - Q[1] qg[1] - q_load[1]) # 节点2注入 支路1末端功率 constraints.append(P[1] pg[2] - p_load[2]) constraints.append(Q[1] qg[2] - q_load[2])对于实际的大型配电网这种逐节点手写方式不可取。常见做法是用关联矩阵 incidence matrix 把支路功率和节点注入功率一次性表达左乘拓扑矩阵得到每个节点的流入减流出再与注入等式匹配。程序里如果已经有邻接表或拓扑矩阵建议直接用矩阵乘法写平衡约束避免在 33 节点里手写几十条等式。3.2 网损最小目标下 SOCP 的松弛紧性最优潮流里可以用不同目标函数。论文中常见的网损最小目标写为$$ \min \sum_{k}(r_k L_k) $$把线路电阻和电流平方变量相乘再求和是网损的标准表达。这个目标对 SOCP 松弛有很好的紧性系统为了降低网损会尽量降低支路电流平方同时又会尽量把电压水平稳定在合理区间所以最优解通常会逼到锥边界上。objective cp.Minimize(cp.sum(cp.multiply(branch[:, 2], L))) prob cp.Problem(objective, constraints) prob.solve(solvercp.ECOS) print(最优网损:, prob.value) print(节点电压平方:, u.value) print(DG出力:, pg.value)代码中的branch[:, 2]取的是 branch 数组第三列也就是每个支路的电阻 r。cp.multiply是 CVXPY 里元素级乘法不能直接用 NumPy 的*操作 Variable 对象。prob.value返回最优目标值u.value是电压平方结果实际电压幅值需要取平方根。求解完成后要立刻检查锥约束的紧性。一个简单做法是计算每条支路的 u_i * L_k - (P_k^2 Q_k^2)。如果这个值对全部支路都近似为 0说明二阶锥松弛是紧的SOCP 解就是原最优潮流问题的可行解如果某条支路的差值明显大于 0说明那条支路的电流或电压被人为放大直接拿去潮流计算会对不上。3.3 提取结果并用潮流校验判断解是否物理可行求解器给出的解需要经过潮流回代验证这一步不能省。原因是 SOCP 只保证在松弛空间里最优不保证在原来的物理可行域里严格成立。校验时我习惯把求解得到的 pg、qg 和根节点电压作为给定值重新用常规潮流计算各节点电压再和 SOCP 结果的 u.value 对比。# 简易校验把 SOCP 解中的注入功率代入 DistFlow 递推 v_sq np.ones(n_nodes) for k, (i, j, r, x) in enumerate(branch): p_ij P.value[k] q_ij Q.value[k] v_sq[j] v_sq[i] - 2 * (r * p_ij x * q_ij) \ (r * r x * x) * L.value[k] print(SOCP u:, u.value) print(回代 u:, v_sq) print(最大偏差:, np.max(np.abs(u.value - v_sq)))这段代码利用 DistFlow 的电压方程从根节点往下递推。如果最大偏差在 1e-4 以下说明 SOCP 解满足电压方程如果偏差达到 1e-2 甚至更大问题通常出在锥约束不紧或者是节点功率平衡约束写错。回代验证不能直接判断锥约束紧性但能快速抓住“电压方程写错”这类低级问题建议放在求解后的第一轮检查里。4. 把多源协同运行优化转到同一套二阶锥框架里4.1 多源协同的调度变量和约束拆分《主动配电网多源协同运行优化研究》不是只解决单时刻最优潮流而是把一组可控设备放到同一个优化框架里协调。多源协同的变量可以分为四类连续有功变量、连续无功变量、储能状态变量、离散投切变量。常见做法是把它们全部写成约束再和上一章的 SOCP 潮流约束共同求解。有功变量包括分布式光伏/风机的有功出力、储能充电功率、储能放电功率、上级电网购电功率。无功变量包括光伏逆变器无功、储能逆变器无功、并联电容器组输出和动态无功补偿装置的无功。储能状态变量主要是 SOC按时间递推$$ SOC_{t1} SOC_t - \frac{\eta_c P_{c,t}\Delta t}{E_{rated}} \frac{P_{d,t}\Delta t}{E_{rated}\eta_d} $$其中 eta_c 是充电效率eta_d 是放电效率P_c 和 P_d 是充放电功率E_rated 是储能额定容量。这个约束是线性的SOCP 框架可以直接容纳。真正带来计算负担的是“储能不能同时充电和放电”这样的互斥约束它需要引入 0-1 整数变量将问题变成混合整数二阶锥规划。配电网侧的离散设备主要包括有载调压变压器分接头和分组投切电容器。这类设备如果按连续变量处理得到的调度指令无法实际执行如果逐个离散档位枚举组合数会爆炸。更自然的做法是把每个分接头档位建模为整数变量并限制相邻时刻最大调节档数。4.2 储能、分布式电源和可调设备的典型参数在搭建多源协同优化程序之前最好把设备参数集中放在一个字典里方便换算成标幺值。下面是一组用于 10 kV 馈线的常见参数设备类型参数数值示例说明储能额定容量1 MWh决定 SOC 递推中的分母储能最大充/放电功率0.25 MW充放电功率上界储能效率0.92 / 0.92充电与放电分开填光伏逆变器额定功率0.5 MVA有功与无功功率耦合上限光伏逆变器无功能力±0.3 Mvar功率因数约 0.6 至 0.95电容器组单组无功0.1 Mvar离散投切整数倍出力上级电网购电功率上限5 MVA根节点注入受限这些参数在程序里通常已经做成可配置对象改参数时要注意单位一致。配电网程序最容易踩坑的地方是阻抗用欧姆、功率用千瓦、电压用千伏混在一起。我一般会把功率统一为 MW阻抗统一为标幺值根节点电压统一为 1.0 p.u.。4.3 从单时段 OPF 扩展到 24 小时多时段模型单时段最优潮流扩展到多时段运行优化时只需要给原本的节点变量、DG 变量和储能变量加一个时间维度。目标函数则改为一天内的购电费用、DG 运行成本和储能退化成本之和。下面给出一个多时段 SOCP 约束的构造片段其中只保留了储能部分的关键约束。import cvxpy as cp T 24 # 24小时 E_rated 1.0 # MWh dt 1.0 # 小时 eff_c 0.92 eff_d 0.92 soc cp.Variable(T 1) # 含初始时刻 p_ch cp.Variable(T, nonnegTrue) # 充电功率 p_dis cp.Variable(T, nonnegTrue) # 放电功率 b_ch cp.Variable(T, booleanTrue) # 1表示充电时段 b_dis cp.Variable(T, booleanTrue) # 1表示放电时段 constraints [soc[0] 0.5, soc[T] 0.5] # 初始与终态SOC各50% for t in range(T): constraints.append(soc[t1] soc[t] - p_ch[t]*dt/(E_rated*eff_c) p_dis[t]*dt/E_rated*eff_d) constraints.append(p_ch[t] 0.25 * b_ch[t]) constraints.append(p_dis[t] 0.25 * b_dis[t]) constraints.append(b_ch[t] b_dis[t] 1) constraints.append(soc[t1] 0.9) constraints.append(soc[t1] 0.1)这段代码里b_ch 和 b_dis 是布尔变量用来保证同一个时段不会同时充电和放电。注意 SOC 递推式的微分项带了效率系数在充电方向电能不能全部存入所以 SOC 增量要除以充电效率在放电方向电池放出的电量来自内部储存所以 SOC 减少量要乘以放电效率。如果把两个效率写反储能会在一天内凭空多出能量必须检查这个细节。将这段储能约束与前面章节中已经写好的 DistFlow 约束合并再给光伏、风机、逆变器无功加上上界就构成了一个完整的主动配电网多源协同运行优化模型。求解器面对的是一个混合整数二阶锥规划问题只有在求解器支持整数变量时才能直接处理这就进入了下一章要讲的求解器选择问题。5. 求解器选择与参数调优MOSEK、Gurobi、ECOS 的适用边界5.1 四种求解器的特点和适用场景同一套 CVXPY 模型交给不同求解器结果和速度差异可以非常大。纯连续变量的 SOCP 和带整数变量的 MISOCP适合的求解器并不完全相同。下表整理了四种常见选择求解器问题类型特点适用规模ECOSSOCP、LP稳妥、默认精度较好不支持整数连续变量小规模算例SCSSOCP、SDP内存占用小采用一阶方法大规模但精度要求不高的场景MOSEKSOCP、MISOCP内点法 分支定界数值最稳工商授权可用配网调度首选GurobiSOCP、MISOCP整数问题求解速度快参数丰富多时段多设备优化如果只是复现论文里的单时段最优潮流ECOS 就够了。如果要做 24 小时联合优化并且带储能、电容器组等整数变量ECOS 会直接报错因为模型里出现了布尔变量。这时可以用 MOSEK 或 Gurobi它们的 MIP 分支定界能力能处理中小规模 MISOCP。程序里如果看到 solver 参数可以配置建议先跑通 ECOS再切换到 MOSEK 对比目标值。5.2 在 CVXPY 里设置相对停止容忍度和 MIP gap多时段模型经常求解时间偏长默认参数不一定最优。以 MOSEK 为例通常需要把锥优化相对对偶间隙从默认值调小保证 SOCP 解的边界质量。# 用 MOSEK 求解纯SOCP收紧相对对偶间隙 prob.solve( solvercp.MOSEK, mosek_params{ MSK_DPAR_INTPNT_CO_TOL_REL_GAP: 1e-8, MSK_DPAR_INTPNT_CO_TOL_NEAR_REL: 1e-8 }, verboseTrue )第一个参数控制内部点法的相对对偶停止阈值越小结果越精确但耗时更长。第二个参数控制近可行解判定阈值配电网潮流约束本身数值尺度较接近设置过严反而会让求解器在后期做无意义迭代。一般先用默认值看目标值再把相对间隙从 1e-6 逐步收紧到 1e-8观察目标值是否还在变化。如果目标值稳定说明当前精度已经足够。使用 Gurobi 处理 MISOCP 时重点调两个参数MIPGap 和 TimeLimit。prob.solve(solvercp.GUROBI, MIPGap1e-4, TimeLimit300, verboseTrue)MIPGap 表示当前整数解与最优下界之间的相对间隙设置为 1e-4 意味着求解器接受 0.01% 的次优性。TimeLimit 给 300 秒上限非常常见配电网多时段模型的整数变量如果特别多不设时间限制可能跑几小时。建议在离线研究里先用 0% 的 gap 算出基准值再在工程部署中把 gap 适当放宽。5.3 常遇到的不可行、数值病态和松弛不成 tight 的处理思路不可行是配电网优化里最常见的问题。出现“infeasible”时先确认负荷和 DG 单位是否都是 MW/Mvar再检查根节点电压约束是否写成 u[0] 1.0 而不是 V[0] 1.0。还要注意节点功率平衡方向在 DistFlow 里支路功率流动方向由支路表序号决定如果拓扑数据里首末端方向不一致平衡等式就会差一个负号。数值病态通常来自量纲。线路电阻用欧姆表示时可能只有 0.5但如果把功率单位写成 kW负荷数值变成几百优化变量尺度相差三个数量级SOCP 的内点法就会在收敛精度上大幅退化。处理办法是统一转成标幺值取基准容量 10 MVA10 kV 配电网阻抗标幺乘以基准值分量负荷除以基准容量电压取 1.0 p.u.。这也是论文程序里常见的预处理方式。松弛不成 tight 的问题往往被忽视。如果目标函数不是网损也不是严格依赖电流平方的成本SOCP 松弛可能给出一个在锥内部的解。解决办法是在目标函数里加上一个小系数网损项例如 0.01 * sum(r_k * L_k)。这样既不会明显改变经济调度结果又能让锥约束在最优解处尽量贴近边界还原出物理可行的电压和功率分布。6. 验证 SOCP 松弛精度的两个实用技巧6.1 用锥边界间隙检查每条支路是否落在边界上拿到一组最优解之后最值得做的第一件事是逐条支路检查锥边界间隙。这项检查能直接说明 SOCP 解是否可信是写论文、改程序时最实用的诊断手段。# 计算每条支路的锥边界间隙 lhs u.value[branch[:, 0].astype(int)] * L.value # u_i * L_k rhs P.value**2 Q.value**2 # 归一化间隙方便观察量级 gap (lhs - rhs) / np.maximum(rhs, 1e-6) print(最大锥间隙:, np.max(gap)) print(最大间隙支路:, np.argmax(gap))如果 gap 总量接近 0说明 SOCP 松弛是紧的最优解可以直接当作真实潮流解使用如果 gap 达到 1e-3 甚至更大说明这条支路对应的锥约束没有触及边界。此时最优解的经济性值可能没问题但电压、电流值并不能真实还原到配电网场景。前文提过给目标函数增加网损项或者把该支路的首端电压固定到某一合理工作点通常可以恢复紧性。6.2 用对偶变量判断配电网瓶颈提升论文分析深度二阶锥规划自带一条其他启发式算法给不了的信息链路对偶变量。在 CVXPY 里求解完成后可以通过constraints[i].dual_value取出对应约束的影子价格。对节点功率平衡约束对偶值表示该节点单位功率增量带来的最优目标变化量也即节点的实时电价或阻塞信号对支路电压约束对偶值能定位电压越限最敏感的支路。多源协同运行优化里这部分信息可以用来解释“为什么储能要在某个时刻充电”这样的结果。比如在某个时段某节点储能放电约束的对偶值很高说明该时刻该节点附近的功率供应紧张储能的边际价值大反之如果对偶值接近 0说明储能动作对目标几乎没有影响可以优先调整其他设备。论文里如果需要做灵敏度分析和经济解释这是一个很容易出图、出表的小技巧。具体导出时注意 CVXPY 对约束取 dual_value 的顺序和 constraints 列表一致。在多时段模型中建议把每个时段的平衡约束单独保存到列表里而不是写在一个大矩阵里否则后面要对单一节点、单一时刻做分析会非常不方便。配合上一节提到的锥边界间隙检查这套流程基本覆盖了从“能出解”到“解可信”的全部关键环节。本文还有配套的精品资源点击获取
