ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

配电网最优潮流新解法:二阶锥松弛原理与Matlab实现

配电网最优潮流新解法:二阶锥松弛原理与Matlab实现 配电网最优潮流OPF这些年被分布式光伏、风机和储能大量接入推到了风口浪尖。真正上手做计算的人都有一个体会潮流约束和目标函数一凑到一起就是一个大规模非线性非凸优化问题传统方法要么算到怀疑人生要么卡在局部最优里出不来。二阶锥松弛SOCP就是专门来解决这个问题的通过变量替换和约束松弛把原本非凸的模型改写成凸优化问题交给成熟求解器一次性拿到全局最优解。这个思路在辐射状配电网里尤其好用实测下来求解速度快、数值稳定精度也完全够工程用。这篇文章不绕弯子直接从原理讲到Matlab代码落地我把每一步为什么这么做、有哪些坑都摊开说适合正在做配电网优化、分布式电源接入分析的硕博研究生和一线工程师。1. 为什么配电网最优潮流要让“非凸变凸”——SOCP方案的整体设计思路1.1 分布式电源接入后的两大痛点双向潮流与电压越限以前的配电网是纯粹的“被动网络”功率从变电站单向流动到末端负荷运行方式相对固定。那时候做潮流计算、可靠性分析用经典的前推回代法就够用了。但现在光伏、风电大规模接入中低压配电网情况完全变了分布式电源出力一高馈线末端的电压可能被顶到越限潮流方向也不再是从上到下单向流动而是可能出现局部反向。更麻烦的是配电网的R/X比远高于输电网有功和无功对电压的影响相互耦合输电网那套“PV节点调无功”的玩法在配电网里经常失灵。因此配电网的运行不能只靠固定的调控策略需要在考虑分布式电源出力的前提下对每个时段、每个节点的功率分配做优化计算这就是配电网最优潮流要解决的问题。它要回答的核心问题是在满足潮流方程、电压限值、支路容量和电源出力约束的前提下如何安排各节点注入功率才能让网损最小、电压质量最好或者分布式电源消纳最多。1.2 传统求解器的天花板局部最优与效率瓶颈最优潮流的传统建模方式是把潮流方程直接写进去得到一个非线性规划NLP问题然后用内点法或序列二次规划去求解。这类方法在小规模输电网里表现不错但放到配电网里就尴尬了。首先是非凸问题。潮流方程里的电压平方项、电流平方项、功率乘积项搅在一起可行域不是凸集。非线性求解器本质上是沿着梯度方向在局部搜索从不同初值出发可能收敛到完全不同的解。你无法确定算出来的结果是不是全局最优。对运行调度来说这很致命你说这个方案是“最优”的但实际可能只找到一个局部较优解损耗高个百分之几电压分布也更差。其次是效率问题。一旦系统规模变大或者要做多时段、多场景的联合优化非线性规划每次迭代都要重新计算海森矩阵和雅可比矩阵计算量成倍增长经常出现计算时间不可接受的情况。也有人尝试用粒子群、遗传算法这类启发式算法来绕开非线性但这类方法本质上靠随机搜索不保证最优性且每次适应度评估都要算一次完整潮流大规模场景下慢得离谱。1.3 二阶锥松弛的核心思想先替换变量再把等式“放松”成锥约束二阶锥松弛的思路很聪明它的核心可以拆成两步。第一步变量替换。把电压幅值的平方记作 (u_i)电流幅值的平方记作 (l_{ij})。这么一换潮流方程里最头疼的平方项大部分都变成线性的了只剩下电流平方和电压平方之间还拖着一个二次等式。第二步把那个二次等式“放松”成一个不等式。原来必须严格满足等式现在允许解落在更大的集合里。这个更大的集合恰好是一个凸锥数学上叫二阶锥。于是整个问题变成二阶锥规划这类问题有极其成熟的求解算法——内点法全局最优性有理论保证求解速度快对大规模问题也能稳定处理。你可以这样理解原问题就像在凹凸不平的山地上找最低点每一步只能靠脚底的感觉摸索凸松弛之后相当于把整个地形磨成一个光滑的碗状你从任何位置出发顺着坡度往下走最后一定会滚到碗底。这里还要说清楚一点这种“放松”并不是无原则的放水。在辐射状配电网的典型条件下最优解会被约束条件“推”到锥的边界上也就是说松弛后的最优解恰恰满足原来的等式。这在数学上叫精确松弛exact relaxation。所以放心用得到的解就是原问题的最优解。2. 从潮流方程到二阶锥核心数学细节逐层拆解2.1 三行方程吃透DistFlow支路潮流模型配电网最优潮流里我们很少用传统的节点导纳矩阵形式而是用支路潮流模型Branch Flow Model也叫DistFlow方程。原因很简单配电网是辐射状结构用“支路—节点”的父子关系来描述功率流动物理含义清晰也方便做凸松弛。对任意一条支路 (i \to j)定义(P_{ij})、(Q_{ij})支路首端从 (i) 流向 (j) 的有功、无功功率(U_i)、(U_j)节点 (i)、(j) 的电压幅值平方(I_{ij})支路电流幅值平方(r_{ij})、(x_{ij})支路电阻和电抗。DistFlow方程写成这样节点功率平衡。对任意节点 (j) [ \sum_{i \in \pi(j)} \left( P_{ij} - r_{ij} I_{ij} \right) P_j^{g} \sum_{k \in \delta(j)} P_{jk} P_j^{d} ] 其中 (\pi(j)) 是 (j) 的父节点集合(\delta(j)) 是子节点集合。这个式子的物理含义是从父支路流入本节点的功率减去支路本身消耗的功率等于流向所有子支路的功率之和加上本地负荷再减去本地分布式电源的注入。无功同理。电压降落方程 [ U_j U_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) I_{ij} ] 这本质上是欧姆定律和功率定义的结合描述了线路两端的电压幅值平方差。电流—功率—电压的耦合方程 [ I_{ij} U_i P_{ij}^2 Q_{ij}^2 ]前两条方程已经相当线性了真正让整个模型变成非凸的就是第三条这个二次等式。2.2 变量替换 (u_i) 和 (l_{ij})把二次项降到一等式为了把模型整得更好看我们直接引入新变量[ u_i U_i, \quad l_{ij} I_{ij} ]注意这里的 (u_i)、(l_{ij}) 本身就是“平方值”不是幅值。然后重写DistFlow方程节点功率平衡变成 [ \sum_{i \in \pi(j)} \left( P_{ij} - r_{ij} l_{ij} \right) P_j^{g} \sum_{k \in \delta(j)} P_{jk} P_j^{d} ] 全部是线性约束。电压降落方程变成 [ u_j u_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) l_{ij} ] 也全部是线性约束。唯一剩下的非线性约束 [ l_{ij} u_i P_{ij}^2 Q_{ij}^2 ] 这是一个二次等式就是它让整个问题非凸。所以现在所有“麻烦”都集中在最后一个等式的处理上。只要把这个等式搞定问题就彻底简化了。2.3 从不等式到标准二阶锥Schur补与cone函数推导二阶锥松弛的做法是把等式 (l_{ij} u_i P_{ij}^2 Q_{ij}^2) 直接放松成不等式[ l_{ij} u_i \ge P_{ij}^2 Q_{ij}^2 ]为什么能这么放松从物理上看(l_{ij} u_i) 比 (P_{ij}^2 Q_{ij}^2) 大意味着电压平方和电流平方的乘积大于功率平方和这在数学上是允许的只是把“精确相等”变成“至少这么大”。但光有这个不等式还不行求解器需要的是标准形式的二阶锥。我们把上式两边乘4并做一点代数变形[ (2P_{ij})^2 (2Q_{ij})^2 (l_{ij} - u_i)^2 \le (l_{ij} u_i)^2 ]验证一下右边展开是 (l_{ij}^2 2l_{ij}u_i u_i^2)左边展开是 (4P^2 4Q^2 l^2 -2l u_i u_i^2)两边同时消掉 (l^2u_i^2)就得到 (4P^24Q^2 \le 4l u_i)也就是 (P^2Q^2 \le l u_i)。写成二阶锥的标准形式[ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ l_{ij} - u_i \end{bmatrix} \right|2 \le l{ij} u_i ]这就是一个标准的二阶锥约束。在YALMIP里直接调用cone函数就能表达cone([2*P(k); 2*Q(k); l(k)-u(i)], l(k)u(i))到这里整个模型的凸化完成。第2.1节里那个非凸NLP问题被等价改写成三组线性约束加一组二阶锥约束的SOCP问题可以由Mosek、Gurobi、SeDuMi等求解器在多项式时间内高效求解并且全局最优性有保障。2.4 目标函数与约束的“凸性红线”松弛把可行性区域变得凸了但如果你在目标函数或约束里自己又引入非凸项那就前功尽弃了。这里总结一下常用的目标函数和约束哪些是安全的哪些会破坏凸性。目标/约束写法凸性说明网损最小(\min \sum r_{ij} l_{ij})线性凸最常用的目标直接用电压偏差最小(\min \sum (u_i - u_{ref})^2)二次凸可用但最好转为线性辅助变量电压偏差最小绝对偏差(\min \sum |V_i - 1|)绝对值非线性需引入辅助变量线性化DG出力最大(\min -\sum P_j^{g})线性凸常用于消纳分析DG无功约束(Q_j^{g2} \le S_j^2 - P_j^{g2})二阶锥凸注意是锥不是圆电压上限约束(u_j^{min} \le u_j \le u_j^{max})线性凸一般加0.9~1.1 pu对应平方有一点特别提醒如果你想表达“DG功率因数不小于0.95”可以用线性约束 (Q_j^g \le \tan(\arccos 0.95) P_j^g) 近似或者用二阶锥约束。但千万不要写 (P_j^{g2} Q_j^{g2} S_j^2) 这种圆等式那会立刻把问题打回非凸原形。如果不是特殊情况工程上宁可做保守处理也不要为了精确而引入非凸约束。3. Matlab实操从空白脚本到出结果3.1 环境准备YALMIP、Mosek/Gurobi/SeDuMiMatlab自带的optimoptions和fmincon可以求解非线性规划但对SOCP这种问题没有原生支持。我的建议是装YALMIP作为建模接口再配一个专业凸优化求解器。YALMIP一个Matlab下的优化建模工具箱语法非常简洁SOCP、SDP、LP都支持。从GitHub下载zip包解压后把整个文件夹加入Matlab路径addpath(genpath(D:\tools\yalmip-master)); savepath;Mosek商业求解器对学术用户免费授权SOCP求解速度业界顶级和YALMIP配合非常丝滑。Gurobi同样是顶级商业求解器学术免费对SOCP支持也很好。SeDuMi、SDPT3开源免费求解器适合小规模算例速度比Mosek/Gurobi慢一些但胜在免费不需要申请授权。安装完求解器后在Matlab里运行一次yalmiptest看到各个求解器状态为OK就说明环境没问题。版本上注意Mosek 10和Gurobi 10需要Matlab R2020a及以上版本我自己用Matlab R2022b配Mosek 10没有遇到兼容问题。3.2 数据准备IEEE 33节点算例的单位与格式IEEE 33节点系统是配电网优化研究最常用的标准算例33个节点、32条支路、5个联络开关基准电压12.66 kV总负荷约3715 kW 2300 kvar。做最优潮流时一般先把联络开关全部打开让它保持纯辐射状结构。数据格式上需要注意单位统一。我习惯用标幺值体系基准功率取10 MVA基准电压12.66 kV则阻抗基准为[ z_{base} \frac{U_{base}^2}{S_{base}} \frac{12.66^2}{10} 16.0276\ \Omega ]支路1-2的电阻是0.0922 Ω换算成标幺值就是0.0922 / 16.0276 0.00575。整个算例的负荷有功约0.3715 pu网损约0.02 pu也就是约200 kW。用这套基准时变量数值都在0.001到1之间求解器数值稳定性很好。完整数据可以从Matpower的case33bw.m中读取也可以网上搜索“IEEE 33节点配电网数据”下载Excel版。数据格式归纳一下支路表首端节点号末端节点号电阻Ω电抗Ω负荷表节点号有功kW无功kvar3.3 核心建模代码逐段拆解变量定义、约束构建、求解与后处理下面给出一个完整的YALMIP建模骨架。基于IEEE 33节点以网损最小为目标无分布式电源场景重点展示核心逻辑。clear; clc; % 1. 基准值与数据 baseMVA 10; % 基准容量 MVA basekV 12.66; % 基准电压 kV zbase basekV^2 / baseMVA; % 支路数据格式: [首端节点 末端节点 电阻Ω 电抗Ω] % 这里仅展示前8条示例完整数据请使用case33bw branch_mat [ 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 1.7114 1.2351 8 9 1.0300 0.7400 ]; % 负荷数据格式: [节点号 有功kW 无功kvar] load_mat [ 2 100 60 3 90 40 4 120 80 5 60 30 6 60 20
返回列表