基于改进灰狼平衡算法的配电网功率潮流和电压优化

Abstract

  Nowadays, existing heuristic optimisation algorithms have disadvantages that they are easily trapped into the local optimum and appear a poor diversity performance of the Pareto front when handling multi-objective power flow and voltage optimisation problems. In order to overcome the aforementioned disadvantages, an improved grey wolf equilibrium optimiser (GWEO) is created in this study by integrating the search mechanism of grey wolf optimisation (GWO) and equilibrium optimiser (EO). In the improved GWEO, the disturbance component of EO is mixed into the search mechanism of GWO, which can further expand the search range of the algorithm and effectively boost the convergence and the diversity of the search results. Next, GWEO is utilised to multi-objective power flow and voltage optimisation for the distribution network with the penetration of distributed photovoltaics (PVs) where the minimal system losses, the minimal average voltage deviation, the minimal curtailment rate are optimisation objectives, and IEEE 12-bus, 33-bus, and 118-bus distribution system are employed as the benchmark networks. The results reveal that power flow and voltage optimisation solutions with higher convergence and diversity can be obtained by the improved GWEO, which can be better applied to the actual power flow and voltage optimisation scenarios for the distribution network.

Keywords

  Power flow and voltage optimisation; distributed photovoltaic; grey wolf equilibrium optimiser (GWEO); multi-objective optimisation

摘要

    针对现有群体智能优化算法在处理多目标功率潮流和电压优化问题时出现的易陷入局部最优、Pareto前沿分布性能不佳等问题，本文结合灰色狼群算法(GWO)和平衡优化器(EO)算法的搜索机制，开发出改进的灰狼平衡算法(GWEO)。该算法在GWO搜索机制的基础上加入了EO的扰动机制，进一步扩大了算法的搜索范围，有效提升了算法搜索结果的收敛性和分布性。然后以IEEE12、33、118节点配电系统为算例，以系统损耗最小、平均电压偏差最小和分布式光伏弃光率最小为优化目标，将GWEO应用到考虑分布式接入场景的配电网多目标功率潮流和电压优化中。优化结果表明，改进的GWEO可为决策者提供更优质、更多样的功率潮流和电压优化方案，因此更能满足实际配电网功率潮流和电压优化的场景需求。

关键词

    功率潮流和电压优化；分布式光伏；灰狼平衡算法(GWEO)；多目标优化

0 引言

    为实现“碳达峰、碳中和”的战略目标，推动分布式光伏接入配电网成为发展以可再生新能源为主体的新型电力系统、实现“双碳”目标和乡村振兴战略的重要举措。然而，分布式光伏的大量接入配电网将会引起配电网内部的有功功率波动，进而引起电网损耗增加、有功功率波动、末端电压偏高、首端功率潮流倒送、三相不平衡等一系列电能质量问题[1]，对配电系统的安全、稳定、经济运行带来了全新的挑战。因此，对考虑分布式光伏接入的配电网进行功率潮流和电压优化已成为保障配电网电能质量和安全稳定运行的迫切需求。

    在考虑分布式光伏接入配电网的工况下，为保证配电系统运行的安全性、稳定性和经济性，一般需要综合系统损耗最小、平均电压偏差最小、弃光率最低等多个优化目标制定配电网功率潮流和电压优化方案，而不同优化目标之间的相互冲突使得多目标功率潮流和电压优化模型的求解结果不唯一，这也增大了问题求解的难度和复杂度[5]。因此，在实际的功率潮流和电压优化领域的研究中不可避免地需要讨论模型的求解问题，而优化算法的选取则是提高求解结果准确性和求解速度的关键。目前，应用于功率潮流和电压优化领域的算法可分为以单个体目标函数值计算和迭代求解为主要思想的传统优化算法和以群体适应度评估和迭代更新为主要思想的智能优化算法两大类[2]。其中，智能优化算法的核心思想是随机生成一定规模的初始种群并计算种群的适应度，然后按照一定的更新规则朝着当前适应度最优的方向迭代更新，最终收敛于全局最优解。相比于传统优化算法，智能优化算法具有迭代规则简单、不依赖初始值、无需要求目标函数和约束条件连续可微等方面的优势，更适用于对多峰值、非线性的复杂优化问题求解，而结合多目标Pareto支配理论的智能优化算法也已广泛应用于高维多目标的复杂应用场景[4]。

    经过国内外学者数十年的研究，现已开发出遗传算法(Genetic Algorithm，GA)、进化算法(Evolutionary Algorithm，EA)、粒子群算法(Particle Swarm Optimisation，PSO)、免疫算法(Immune Algorithm，IA)、蚁群算法(Ant Colony Optimisation，ACO)、烟花爆炸算法(Fireworks Explosion Algorithm，FEA)、灰色狼群算法(Grey Wolf Optimisation，GWO)[8]、樽海鞘群算法(Salp Swarm Algorithm，SSA)、平衡优化器算法(Equilibrium Optimiser，EO)[12]等诸多类型的群体智能优化算法，但在配电网多目标功率潮流和电压优化的应用场景中，由于群体更新方法和外部档案维护机制的局限性，上述算法普遍存在收敛速度慢、易陷入局部最优、Pareto前沿分布性能不佳的问题。针对以上问题，本文结合灰色狼群算法和平衡优化器算法的搜索机制，开发出改进的灰狼平衡算法(Grey Wolf Equilibrium Optimiser，GWEO)并将其应用到考虑分布式接入场景的配电网多目标功率潮流和电压优化中，最终以三种典型的配电系统为仿真算例验证了该算法的搜索结果具有良好的收敛性和分布性。

1 配电网功率潮流和电压优化模型

1.1 控制变量

    对于考虑分布式光伏接入的配电系统，其功率潮流和电压优化控制的方式主要包括调节有载调压变压器抽头、无功补偿装置的补偿容量和分布式光伏设备的有功和无功出力等，因此以各时段有载调压变压器变比$K_{T,i}(t)$、并联补偿装置的无功补偿量$Q_{C,i}(t)$、分布式光伏设备的有功出力$P_{DG,i}(t)$和无功出力$Q_{DG,i}(t)$作为优化模型中的控制变量，并对控制变量设置如式(1)所示的上下限约束
$$\begin{array}{c}\{\begin{array}{l}{K}_{T,i,min}(t)\le K_{T,i}(t)\le K_{T,i,max}(t)\\ Q_{C,i,min}(t)\le Q_{C,i}(t)\le Q_{C,i,max}(t)\\ P_{DG,i,min}(t)\le P_{DG,i,i}(t)\le P_{DG,i,max}(t)\\ Q_{DG,i,min}(t)\le Q_{DG,i}(t)\le Q_{DG,i,max}(t)\end{array}\end{array}$$
其中：$N_G$、$N_T$、$N_C$、$N_{DG}$分别为电力系统内的发电机节点总数、有载调压变压器可调抽头总数、并联电容器总组数、分布式电源节点总数。

1.2 目标函数

1.2.1 系统损耗

    为保证发出的功率能够在配电网内被有效利用，提高配电网运行的经济性，应按照如式(2)所示的形式将全时段内配电网平均系统损耗最小作为其中一个优化目标
$$\begin{array}{c}\begin{array}{rl}\text{min}{F}_{loss}& =\frac{1}{N_t}\sum _{t=1}^{N_T}{P}_{loss}(t)\\ & =\frac{1}{N_t}\sum _{t=1}^{N_T}\left[\sum _{(i,j)\in L}^{N_L}{G}_{ij}(t)\left(V_i^2(t)+V_j^2(t)-2V_i(t)V_j(t)\cos {\theta }_{ij}(t)\right)\right]\end{array}\end{array}$$
其中：$N_L$分别为电力系统内的电力线路总数，$N_t$为模型的时段总数，$P_{loss}(t)$为第$t$个时段内电网的总损耗，$L$为电力系统中所有电力线路的集合，$G_{ij}(t)$和$B_{ij}(t)$分别为第$t$个时段内节点导纳矩阵中节点$i$和节点$j$之间的电导和电纳，$V_i(t)$和${\theta }_{ij}(t)$为第$t$个时段内节点i处的电压有效值和节点$i$和节点$j$之间的电压相位差。

1.2.2 电压偏差

    为保证配电系统的安全稳定运行，应保证配电网内各节点的电压均能维持在额定值附近，因此需按照如式(3)所示的形式将配电网内各节点的平均电压偏差最小作为优化目标
$$\begin{array}{c}\begin{array}{rl}\text{min}{F}_{\Delta V}& =\frac{1}{N_t}\sum _{t=1}^{N_T}\frac{1}{N_B}\sum _{i=1}^{N_B}\frac{\left|V_i(t)-V_{i,base}\right|}{V_{i,base}}\times 100\%\\ & =\frac{1}{N_t}\sum _{t=1}^{N_T}\frac{1}{N_B}\sum _{i=1}^{N_B}\left|V_i^{\ast }(t)-1\right|\times 100\%\end{array}\end{array}$$
其中：$V_{i,base}$为节点$i$处的基准电压，$V_i^{\ast }(t)$为节点$i$处的电压标幺值。

1.2.3 分布式光伏利用率

    为提高分布式光伏利用率，使得分布式光伏设备能够尽可能多地向电网输送清洁的电能，应将分布式光伏平均弃光率最小作为其中的一个优化目标，该优化目标的计算方法如式(4)所示
$$\begin{array}{c}\text{min}{F}_{PV}=\frac{1}{N_{PV}}\sum _{t=1}^{N_T}\frac{1}{\left|T_{PV}\right|}\sum _{i=1}^{N_B}\frac{P_{PV,i}(t)}{P_{PV,i,max}}\times 100\%\end{array}$$
其中：$P_{PV,i}(t)$为第$t$个时段内第$i$个分布式光伏设备实际发出的有功功率，$P_{PV,i,max}(t)$为第$t$个时段内第$i$个分布式光伏设备能够向配电网发出的最大有功功率，$T_{PV}$为由所有分布式光伏设备最大有功出力大于$0$的时段组成的集合，$\left|T_{PV}\right|$为集合$T_{PV}$包含的元素个数。

1.3 约束条件

1.3.1 等式约束

    配电系统需满足如式(5)所示的配电系统功率潮流平衡约束，该约束是配电网功率潮流和电压优化模型等式约束的主要组成部分
{ P G i ( t ) − P L i ( t ) = V i ( t ) ∑ j = 1 , j ≠ i N B V j ( t ) ( G i j ( t ) cos ⁡ θ i j ( t ) + B i j ( t ) sin ⁡ θ i j ( t ) ) Q G i ( t ) − Q L i ( t ) + Q C i ( t ) = V i ( t ) ∑ j = 1 , j ≠ i N B V j ( t ) ( G i j ( t ) sin ⁡ θ i j ( t ) − B i j ( t ) cos ⁡ θ i j ( t ) ) i ∈ { 1 , 2 , … , N B } \begin{equation} \begin{cases} P_{Gi}(t)-P_{Li}(t) &= V_i(t) \displaystyle\sum_{j=1,j≠i}^{N_B}V_j(t){\left(G_{ij}(t)\cos{θ_{ij}(t)}+B_{ij}(t)\sin{θ_{ij}(t)}\right)} \\ Q_{Gi}(t)-Q_{Li}(t)+Q_{Ci}(t) &= V_i(t) \displaystyle\sum_{j=1,j≠i}^{N_B}V_j(t){\left(G_{ij}(t)\sin{θ_{ij}(t)}-B_{ij}(t)\cos{θ_{ij}(t)}\right)} \end{cases} i∈ \{1,2,…,N_B\} \end{equation} ⎩ ⎨ ⎧​PGi​(t)−PLi​(t)QGi​(t)−QLi​(t)+QCi​(t)​=Vi​(t)j=1,j=i∑NB​​Vj​(t)(Gij​(t)cosθij​(t)+Bij​(t)sinθij​(t))=Vi​(t)j=1,j=i∑NB​​Vj​(t)(Gij​(t)sinθij​(t)−Bij​(t)cosθij​(t))​i∈{1,2,…,NB​}​​
其中：$P_{Gi}(t)$和$Q_{Gi}(t)$分别为第$t$个时段内节点$i$处发电设备发出的总有功功率和总无功功率，$P_{Li}(t)$和$Q_{Li}(t)$分别为第$t$个时段内节点$i$处负荷消耗的总有功功率和总无功功率，$Q_{Ci}(t)$为第$t$个时段内节点$i$处无功补偿设备补偿的无功功率。

1.3.2 不等式约束

    为保障配电系统的安全稳定运行，配电网各节点处的电压、各节点处的总进线功率因数、各支路的电流应满足如式(6)和式(7)所示的不等式约束
{ V i , m i n ≤ V i ( t ) ≤ V i , m a x P G i ( t ) − P L i ( t ) [ P G i ( t ) − P L i ( t ) ] 2 + [ Q G i ( t ) − Q L i ( t ) + Q C i ( t ) ] 2 ≥ P F m i n i ∈ { 1 , 2 , … , N B } \begin{equation} \begin{cases} V_{i,min} ≤ V_i(t) ≤ V_{i,max} \\ \dfrac{P_{Gi}(t)-P_{Li}(t)}{\sqrt{\left[P_{Gi}(t)-P_{Li}(t)\right]^2 + \left[Q_{Gi}(t)-Q_{Li}(t)+Q_{Ci}(t)\right]^2}} ≥ PF_{min} \end{cases} i∈ \{1,2,…,N_B\} \end{equation} ⎩ ⎨ ⎧​Vi,min​≤Vi​(t)≤Vi,max​[PGi​(t)−PLi​(t)]2+[QGi​(t)−QLi​(t)+QCi​(t)]2 ​PGi​(t)−PLi​(t)​≥PFmin​​i∈{1,2,…,NB​}​​
$$\begin{array}{c}{I}_{(i,j)}(t)\le I_{(i,j)max}(t)\end{array}$$
其中：$I_{(i,j)}(t)$为第$t$个时段内支路$(i,j)$的电流。

    综上，配电网功率潮流和电压优化模型可抽象为如式(8)所示的形式
min ⁡ y = f ( x ) = [ f 1 ( x ) , . . . , f m ( x ) ] s.t. { x ‾ ≤ x ≤ x ‾ g i ( x ) ≤ 0 , i ∈ { 1 , 2 , … , p } h j ( x ) ≤ 0 , j ∈ { 1 , 2 , … , q } \begin{equation} \begin{matrix} \min & \boldsymbol{y} = \boldsymbol{f}(\boldsymbol{x}) = \left[\boldsymbol{f_1}(\boldsymbol{x}),...,\boldsymbol{f_m}(\boldsymbol{x})\right]\\ \text{s.t.} & \begin{cases} \underline{\boldsymbol{x}} ≤ \boldsymbol{x} ≤ \overline{\boldsymbol{x}} \\ g_i(\boldsymbol{x}) ≤ 0, i∈ \{1,2,…,p\} \\ h_j(\boldsymbol{x}) ≤ 0, j∈ \{1,2,…,q\} \end{cases} \end{matrix} \end{equation} mins.t.​y=f(x)=[f1​(x),...,fm​(x)]⎩ ⎨ ⎧​x​≤x≤xgi​(x)≤0,i∈{1,2,…,p}hj​(x)≤0,j∈{1,2,…,q}​​​​
其中：$\mathbf{y}$为由$m$个优化目标组成的目标向量，$\mathbf{x}$为维度为$n$的控制变量，$\underline{\mathbf{x}}$和$\overline{\mathbf{x}}$分别为控制变量的下界和上界，$p$为所有不等式约束的个数，$q$为所有等式约束的个数。

2 灰狼平衡算法

2.1 灰色狼群算法

    灰色狼群算法(Grey Wolf Optimisation，GWO)[8]是一种模拟自然界中狼群等级制度和觅食机制的群体智能优化算法，于2014年首次提出，其更新方法如式(9)所示
$$\begin{array}{c}\{\begin{array}{ll}{\mathbf{D}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)& ={\mathbf{C}}_{1,2,3}\circ {\mathbf{X}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)-\mathbf{X}(k)\\ {\mathbf{X}}_{1,2,3}(k)& ={\mathbf{X}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)-{\mathbf{A}}_{1,2,3}(k)\circ {\mathbf{D}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)\\ \mathbf{X}(k+1)& =\frac{{\mathbf{X}}_1(k)+{\mathbf{X}}_2(k)+{\mathbf{X}}_3(k)}{3}\end{array}\end{array}$$

图1 灰狼攻击猎物和搜寻猎物 Figure 1: Attacking prey versus searching for prey

其中：$k$为算法的当前迭代次数，$K$为算法的最大迭代次数，$\circ$为元素群运算的积；$\mathbf{X}(k)$、$\mathbf{X}(k+1)$分别为第$k$代和第$k+1$代群体中任一灰狼个体对应的$n$维控制变量，${\mathbf{X}}_{\mathbf{\alpha }}(k)$、${\mathbf{X}}_{\mathbf{\beta }}(k)$、${\mathbf{X}}_{\mathbf{\delta }}(k)$分别为$\alpha$个体、$\beta$个体、$\delta$个体对应的$n$维控制变量，$\alpha$个体、$\beta$个体、$\delta$个体为第$k$次迭代后筛选出的适应度排名前三位的个体；$\mathbf{C}$为$n$维随机向量，服从区间$[-2,2]$上的均匀分布，主要用于体现灰狼群体搜索猎物的随机性；$\mathbf{A}(k)$为$n$维随机向量，服从区间$[-a,a]$上的均匀分布；如图1所示，当$|\mathbf{A}(k)|<1$时灰狼个体靠近猎物，$|\mathbf{A}(k)|>1$时灰狼个体远离猎物并搜寻其它猎物，该变量可扩大狼群的搜索范围，避免算法陷入局部最优，在多目标优化搜索的场景中还能提升Pareto前沿分布范围的广度；参数$a$按照如式(10)所示的规则计算，其数值随着迭代次数的增加由$2$线性减小到$0$。
$$\begin{array}{c}a=2\left(1-\frac{k}{K}\right)\end{array}$$

2.2 平衡优化器算法

    平衡优化器算法(Equilibrium Optimiser，EO)[12]是根据流体在控制体积容器内的流动过程而开发的一种新型群体智能优化算法，其更新方法如式(11)所示
X e q ( k ) ∈ { X e q 1 ( k ) , X e q 2 ( k ) , X e q 3 ( k ) , X e q 4 ( k ) } X ( k + 1 ) = X e q ( k ) − [ X ( k ) − X e q ( k ) ] ∘ F + G λ ∘ ( 1 − F ) \begin{equation} \begin{split} \boldsymbol{X_{eq}}(k) &∈ \left\{\boldsymbol{X_{eq1}}(k), \boldsymbol{X_{eq2}}(k), \boldsymbol{X_{eq3}}(k), \boldsymbol{X_{eq4}}(k)\right\} \\ \boldsymbol{X}(k+1) &= \boldsymbol{X_{eq}}(k)-\left[\boldsymbol{X}(k)-\boldsymbol{X_{eq}}(k)\right]∘\boldsymbol{F}+\dfrac{\boldsymbol{G}}{\boldsymbol{λ}}∘\left(1-\boldsymbol{F}\right) \end{split} \end{equation} Xeq​(k)X(k+1)​∈{Xeq1​(k),Xeq2​(k),Xeq3​(k),Xeq4​(k)}=Xeq​(k)−[X(k)−Xeq​(k)]∘F+λG​∘(1−F)​​​
其中：${\mathbf{X}}_{\mathbf{eq}1}(k)$、${\mathbf{X}}_{\mathbf{eq}2}(k)$、${\mathbf{X}}_{\mathbf{eq}3}(k)$、${\mathbf{X}}_{\mathbf{eq}4}(k)$为第$k$次迭代后筛选出的适应度排名前四位个体对应的n维控制变量，${\mathbf{X}}_{\mathbf{eq},\mathbf{ave}}(k)$为上述四者的平均值，${\mathbf{X}}_{\mathbf{eq}}(k)$为从以上五个个体中进一步随机筛选得出的精英个体对应的$n$维控制变量；$\mathbf{\lambda }$为服从单位均匀分布的$n$维随机向量；$\mathbf{F}$为$n$维随机向量，满足如式(12)所示的关系；$\mathbf{G}$为$n$维随机向量，满足如式(13)所示的关系
$$\begin{array}{c}\mathbf{F}=2\text{sign}(\mathbf{r}-0.5)\circ ({\text{e}}^{-l\mathbf{\lambda }}-1)\end{array}$$
$$\begin{array}{c}\begin{array}{rl}\mathbf{G}& ={\mathbf{G}}_0\circ \mathbf{F}\\ {\mathbf{G}}_0& =\mathbf{GCP}-{\mathbf{A}}_{1,2,3}(k)\circ {\mathbf{D}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)\\ \mathbf{GCP}& =0.5r_1(r_2\ge 0.5)\text{ones}(1,n)\end{array}\end{array}$$

其中：$\text{sign}$为符号函数，$\mathbf{r}$为服从单位均匀分布的$n$维随机向量，$r_1$和$r_2$均为服从单位均匀分布的随机数，参数$l$满足如式(14)所示的关系，$\text{ones}(1,n)$为$n$维全$1$行向量。
$$\begin{array}{c}l={\left(1-\frac{k}{K}\right)}^{\frac{k}{K}}\end{array}$$
    EO引入了随机向量$\mathbf{F}$和$\mathbf{G}$，增大了算法跳出局部最优的概率，在一定程度上可提高群体的多样性，从而保证算法搜索结果更加逼近于真实的最优值。

2.3 改进的灰狼平衡算法

2.3.1 群体更新方法改进

    根据上述对GWO和EO更新机制的对比，不难发现：GWO中随机向量$\mathbf{A}$内每个元素的取值范围为$[-2,2]$，而EO中随机向量F内每个元素的取值范围为$[-1.2642,1.2642]$，因此EO的搜索范围相对GWO更为狭窄，这将对EO算法搜索结果的收敛性和多目标Pareto前沿的分布性造成不利影响；GWO缺少随机的扰动因子，这将增加算法陷入局部最优的概率，进而使得GWO的收敛性弱于EO。为解决上述两种优化算法在搜索过程中存在的问题，本文结合如式(9)和式(11)所示的GWO和EO的更新机制，开发出改进的灰狼平衡算法(Grey Wolf Equilibrium Optimiser，GWEO)，该算法采用如式(15)所示的方法进行群体更新
$$\begin{array}{c}\{\begin{array}{ll}{\mathbf{D}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)& =\left|{\mathbf{C}}_{1,2,3}\circ {\mathbf{X}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)-\mathbf{X}(k)\right|\\ {\mathbf{X}}_{1,2,3}(k)& ={\mathbf{X}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)-{\mathbf{A}}_{1,2,3}(k)\circ {\mathbf{D}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}(k)+\frac{{\mathbf{G}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}}{{\mathbf{\lambda }}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}}\circ \left(1-{\mathbf{F}}_{\mathbf{\alpha },\mathbf{\beta },\mathbf{\delta }}\right)\\ \mathbf{X}(k+1)& =\frac{{\mathbf{X}}_1(k)+{\mathbf{X}}_2(k)+{\mathbf{X}}_3(k)}{3}\end{array}\end{array}$$
    与现有的GWO和EO相比，改进的GWEO结合了EO在收敛性上的优势和GWO在多目标优化场景中分布性能上的优势，能够为决策者提供更为优质和多样的决策方案。

2.3.2 参数计算方法改进

    通过如图1所示的灰狼群体狩猎过程和如式(15)所示的群体更新公式，可知参数$a$的大小将会改变灰狼群体远离猎物的概率，进而改变灰狼群体搜索的范围，从而影响迭代结果的收敛性。

    设灰狼远离猎物为概率事件$E$，迭代过程恰好进行到第$k$次为概率事件$I_k$，则在算法的整个迭代搜索的过程中事件$E$发生的概率$p(E)$满足
$$\begin{array}{c}\begin{array}{rl}p(E)& =\sum _{k=0}^{K-1}p(E\cap I_k)\\ & =\sum _{k=0}^{K-1}p(E|I_k)p(I_k)\\ & =\frac{1}{K}\sum _{t=1}^{N_T}p(E|I_k)\end{array}\end{array}$$
其中：概率事件$E\cap I_k$为概率事件$E$与概率事件$I_k$的交事件，用以表示迭代过程恰好进行到第k次且灰狼远离猎物；$p(E|I_k)$表示迭代过程进行到第$k$次时灰狼远离猎物的概率。当$A$服从在$[-a,a]$上的均匀分布时，$p(E|I_k)$满足如式(17)所示的关系
$$\begin{array}{c}p(E|I_k)=\{\begin{array}{l}\begin{array}{cc}\frac{a-1}{a}& a\ge 1\\ 0& a<1\end{array}\end{array}\end{array}$$

    在群体智能优化算法的搜索过程中，迭代前期的个体搜索范围应相对较大，以防止群体智能算法陷入早熟收敛；迭代后期的个体搜索范围应相对较小，以使得群体中的每个个体都能够收敛到搜寻的最优结果。

    对于改进的GWEO算法，参数$a$的取值与算法的搜索范围呈正相关，因此参数a应随着迭代次数的增加而减小。此外，参数$a$不再随迭代次数线性递减，而是以如式(18)和图2所示的规律呈非线性减小的趋势
$$\begin{array}{c}\begin{array}{rl}a& =2l\\ l& ={\left(1-\frac{k}{K}\right)}^{\frac{\mu k}{K}}\end{array}\end{array}$$

其中：$k$为算法的当前迭代次数，$K$为算法的最大迭代次数，$\mu$为比例因子。

图2 参数a随迭代次数的变化趋势 Figure 2: Variations of a with the number of iterations

    设优化算法的最大迭代次数$K=200$，则当参数$a$以如式(13)所示的形式计算时，灰狼远离猎物的概率随迭代次数的变化趋势如图3所示。

图3 灰狼远离猎物的概率随迭代次数的变化趋势 Figure 3: Variations of the probability of each grey wolf searching for prey with the number of iterations

    根据如式(16)所示的方法，可计算出参数$a$随迭代次数的增加线性递减以及按如式(18)所示的形式非线性递减时灰狼远离猎物的概率，其计算结果如表1所示。

表1 灰狼远离猎物的概率 Table 1: Probability of each grey wolf searching for prey

参数a的计算方法		p(E)
依照式(10)线性递减		0.1547
依照式(18)非线性递减	μ = 1	0.2479
	μ = 2	0.1861
	μ = 3	0.1560
	μ = 4	0.1372
	μ = 5	0.1240

    结合图3和表1，本文选取参数$\mu =2$用于参数$a$的更新计算，其主要目的为：

1. 增强算法前期的全局搜索能力：此时灰狼群体远离猎物的概率较线性递减时更高，算法前期陷入局部最优的概率有所降低，从而有利于提升算法搜索结果的收敛性；
2. 实现前期全局搜索和后期局部搜索的统筹兼顾：当参数$\mu =2$时，GWEO算法的搜索过程恰好以迭代过程进行到迭代中期时为界，这与参数$a$随迭代次数的增加线性递减时的情况一致；当迭代次数$k<\frac{K}{2}$时，参数$a$的取值在$1$到$2$之间，此时灰狼群体能够以一定的概率远离猎物，可防止算法在迭代前期陷入局部收敛；当迭代次数$k>\frac{K}{2}$时，参数$a$的取值在$1$以下，此时灰狼群体全部向猎物靠近，可保证算法在迭代后期具有良好的局部搜索能力。
3. 若均为不可行个体，如果它们的相对可行度相等，则比较它们的约束违约度，违反约束程度小的个体占优；否则，相对可行度大的个体占优。

3 基于Pareto支配理论的多目标优化

3.1 多目标Pareto支配理论

    设集合$X_f$为$n$维决策空间中所有可行解组成的集合，定义Pareto支配关系为：对于${\mathbf{x}}_1,{\mathbf{x}}_2\in X_f$，当且仅当对 ∀ i ∈ { 1 , 2 , … , m } ∀i ∈ \{1,2,…,m\} ∀i∈{1,2,…,m}都有$f_i({\mathbf{x}}_1)\le f_i({\mathbf{x}}_2)$，且 ∃ j ∈ { 1 , 2 , … , m } ∃j ∈ \{1,2,…,m\} ∃j∈{1,2,…,m}使得$f_i({\mathbf{x}}_1)<f_i({\mathbf{x}}_2)$。此时称作${\mathbf{x}}_1$支配${\mathbf{x}}_2$，记为${\mathbf{x}}_1\prec {\mathbf{x}}_2$。

    进一步地，定义Pareto最优解为：对于${\mathbf{x}}^{\ast }\in X_f$，当且仅当$\neg \exists \mathbf{x}\in X_f$，使得$\mathbf{x}\prec {\mathbf{x}}^{\ast }$时，则称${\mathbf{x}}^{\ast }$为Pareto最优解，亦称非劣解或非支配解。更进一步地，将由所有Pareto最优解组成的集合称为Pareto最优解集，亦可称为非劣解集或非支配解集；然后，将Pareto最优解集向目标空间的映射定义为Pareto前沿。

3.2 基于栅格的精英个体筛选

    对于多目标优化问题，由于不同优化目标之间相互冲突、相互影响，基于Pareto支配机制的多目标群体智能优化算法在每次迭代过程结束后一般会得到多个非支配个体，因此只需在非支配解集中筛选出精英个体用于群体更新即可。而在每次迭代更新的过程中，应尽量避免从非支配解集中筛选密度过高的个体，以防止过多的个体分布在目标空间中的某个狭小的区域内，从而降低多目标群体智能优化算法搜索结果的分布性能。

    据此，本文采用文献[6]中基于自适应栅格的排序方法在外部档案中筛选出三个个体作为GWEO算法中的$\alpha$个体、$\beta$个体、$\delta$个体，再以精英个体的位置为基准更新群体中每个个体，从而使搜索结果能够以有限的迭代次数逼近于真实的Pareto前沿。此外，为避免算法搜索结果陷入全局最优，在每次迭代过程中应尽量不重复地筛选出多个个体作为精英个体。

3.3 外部档案的更新和多样性维护

    在多目标群体智能优化算法中，一般需要将算法迭代过程中得到的所有非支配个体保存在专门的存储库（即外部档案）中。由于计算速度和存储容量的限制，当外部档案内非支配解的数目达到预设的最大规模时，必须对外部档案进行维护，以保证外部档案的规模始终不超过预设值。

    据此，本文采用如下所示的更新机制更新外部档案：在每次迭代过程完成后，利用Pareto支配理论比较新一代个体与外部档案中的已有个体，将不受已有个体支配的所有新一代个体存入外部档案，并从外部档案中剔除所有受支配的已有个体；当外部档案数据容量溢出时，则计算外部档案中每个个体的个体密度，然后将外部档案中个体密度最小的个体丢弃，直至外部档案的规模等于预设的最大值。

    为提高多目标优化算法搜索结果在目标空间上分布的广泛性和均匀性，本文采用文献[4]中基于拥挤距离的计算方法计算外部档案中每个个体的个体密度，将极端个体的拥挤距离设为无穷大，其它个体的拥挤距离等于两个相邻解的函数值的绝对归一化差。当外部档案规模溢出时，个体密度越大的个体越容易被从外部档案中剔除，以防止Pareto前沿过于集中地分布在目标空间中的某些特定区域。

4 配电网功率潮流和电压优化问题的求解

4.1 控制变量的编码

    配电网功率潮流和电压优化问题中，应对各时段有载调压变压器变比、并联补偿装置的无功补偿量、分布式光伏设备的有功和无功出力进行优化，并采用如式(19)所示的编码方式将其转化为式(8)中的控制变量$\mathbf{x}$
$$\begin{array}{c}\mathbf{x}=\left[\begin{array}{ccccc}\mathbf{x}(1)& \cdots & \mathbf{x}(t)& \cdots & \mathbf{x}(N_T)\end{array}\right]\end{array}$$

其中：$\mathbf{x}(t)$对应第$t$个时段下的功率潮流和电压优化策略，满足
$$\begin{array}{c}\mathbf{x}(t)=\left[\begin{array}{cccc}{\mathbf{K}}_{\mathbf{T}}(t)& {\mathbf{Q}}_{\mathbf{C}}(t)& {\mathbf{P}}_{\mathbf{DG}}(t)& {\mathbf{Q}}_{\mathbf{DG}}(t)\end{array}\right]\end{array}$$
其中：
$$\begin{array}{c}\begin{array}{rl}{\mathbf{K}}_{\mathbf{T}}(t)& =\left[\begin{array}{ccccc}{\mathbf{K}}_{\mathbf{T},1}(t)& \cdots & {\mathbf{K}}_{\mathbf{T},\mathbf{i}}(t)& \cdots & {\mathbf{K}}_{\mathbf{T},{\mathbf{N}}_{\mathbf{T}}}(t)\end{array}\right]\\ {\mathbf{Q}}_{\mathbf{C}}(t)& =\left[\begin{array}{ccccc}{\mathbf{Q}}_{\mathbf{C},1}(t)& \cdots & {\mathbf{Q}}_{\mathbf{C},\mathbf{i}}(t)& \cdots & {\mathbf{Q}}_{\mathbf{C},{\mathbf{N}}_{\mathbf{T}}}(t)\end{array}\right]\\ {\mathbf{P}}_{\mathbf{DG}}(t)& =\left[\begin{array}{ccccc}{\mathbf{P}}_{\mathbf{DG},1}(t)& \cdots & {\mathbf{P}}_{\mathbf{DG},\mathbf{i}}(t)& \cdots & {\mathbf{P}}_{\mathbf{DG},{\mathbf{N}}_{\mathbf{T}}}(t)\end{array}\right]\\ {\mathbf{Q}}_{\mathbf{DG}}(t)& =\left[\begin{array}{ccccc}{\mathbf{Q}}_{\mathbf{DG},1}(t)& \cdots & {\mathbf{Q}}_{\mathbf{DG},\mathbf{i}}(t)& \cdots & {\mathbf{Q}}_{\mathbf{DG},{\mathbf{N}}_{\mathbf{T}}}(t)\end{array}\right]\end{array}\end{array}$$

4.2 目标函数和约束条件的处理

    对于第1章中配电网多目标功率潮流和电压优化模型中的待优化目标函数，只需以如式(22)所示的方式将其转化为式(8)中的目标向量$\mathbf{y}$即可。
$$\begin{array}{c}\mathbf{y}=\left[\begin{array}{ccc}{F}_{loss}& F_{\Delta V}& F_{PV}\end{array}\right]\end{array}$$

    功率潮流和电压优化模型的约束条件包括等式约束和不等式约束两个部分。式(5)所示的配电系统潮流平衡方程约束是优化模型等式约束的重要组成部分，只需借助专业的电力系统潮流计算工具进行配电系统的潮流计算即可得出满足潮流平衡方程约束的配电网功率潮流和电压优化策略；而对于如式(6)-(7)所示的不等式约束，本文选用文献[10]提出的可行性比较法处理约束条件，其基本原则为，对于群体中的两个待比较个体：

1. 若均为可行个体，则目标函数值小的个体占支配地位；
2. 若存在一个可行个体和一个不可行个体，则可行个体占支配地位；
3. 若均为不可行个体，如果它们的相对可行度相等，则比较它们的约束违约度，违反约束程度小的个体占优；否则，相对可行度大的个体占优。

4.3 求解流程

    综上所述，利用GWEO处理配电网多目标功率潮流和电压优化问题的求解流程为：

1. 采集并处理配电网内各监测点的运行数据和各电力设备的铭牌参数；
2. 根据运行数据和铭牌参数对配网进行系统建模，建立多目标功率潮流和电压优化控制模型；
3. 设置功率潮流和电压优化算法的初始参数并进行群体初始化，计算初始群体中每个个体的目标函数值并判断其是否满足约束条件，根据多目标Pareto支配原则从初始群体中选取非支配个体存放到外部档案；
4. 结合GWO和EO的选择机制，从外部档案中选择精英个体，根据所提出的群体更新方法得到新一代群体，然后计算新一代群体中每个个体的目标函数值，并判断其是否满足约束条件；
5. 依照Pareto支配原则从新一代群体中筛选出不受支配的个体并更新外部档案，通过多样性维护机制保证外部档案的规模始终不超过预设的最大值；
6. 判断是否满足迭代终止条件，若不满足条件则转到步骤3)继续进行，若满足条件则停止迭代并从外部档案中获得最终的功率潮流和电压优化控制方案。

    据此，采用改进的GWEO求解配电网多目标功率潮流和电压优化问题的具体流程如图4所示。

图4 基于GWEO的配电网多目标功率潮流和电压优化流程图 Figure 4: Flowchart for power flow and voltage optimisation of the distribution network based on GWEO

5 算例验证

5.1 仿真环境

    OpenDSS软件是一款开源的配电系统潮流仿真软件，支持对各类辐射状和网状配电系统进行基波潮流分析、谐波潮流分析、故障分析等相频域分析。此外，OpenDSS仿真软件还提供与MATLAB相对应的接口，因此可在MATLAB平台中调用OpenDSS仿真软件求解配电网功率潮流，从而实现配电系统功率潮流与电压优化。

    本次算例分析中，测试系统的CPU为Intel六核i7-9750H 2.60GHz，RAM大小为16GB，开发环境为MATLAB R2021a，潮流计算软件OpenDSS版本为9.2.0.1。

5.2 算例系统

    为体现改进的GWEO算法在收敛性和Pareto前沿分布性方面的提升效果，本文选用三种典型的辐射状配电系统——改进的IEEE12、33、118节点系统作为算例系统进行验证，其拓扑结构如图5所示，运行参数分别见文献[19]、[20]、[21]。

a) IEEE 12节点系统 a) IEEE 12-bus distribution system

b) IEEE 33节点系统 b) IEEE 33-bus distribution system

c) IEEE 118节点系统 c) IEEE 118-bus distribution system 图5 所选算例系统的拓扑结构 Figure 5: Topologies of selected benchmark networks

    依照文献[23]的要求，本文首先在上述三种算例系统中安装分布式光伏设备和无功补偿装置，其中分布式光伏设备和无功补偿设备的安装位置如图5所示，分布式光伏设备和无功补偿装置的额定容量、优化模型的节点电压约束和功率因数约束如表2所示，分布式光伏设备的典型日有功出力曲线和负荷设备的日功率曲线如图6所示。在完成对分布式光伏设备和无功补偿装置的选址定容后，以1h为时间间隔进行日前功率潮流与电压优化。

表2 三种算例系统的配置参数 Table 2: Configuration parameters of three selected benchmark networks

配置参数	IEEE12节点配电系统	IEEE33节点配电系统	IEEE118节点配电系统
分布式光伏设备的额定容量	100kW	300kW	500kW
无功补偿设备的最大容量	60kvar	200kvar	1002kvar
系统电压等级	10kV
节点电压约束	0.93p.u.~1.07p.u.
功率因数约束	≥ 0.95

图6 日功率波动曲线 Figure 6: Daily power fluctuation curves

5.3 仿真结果

5.3.1 单目标场景下的收敛性

    分别以配电网平均损耗最小和各节点平均电压偏差最小为优化目标，利用现有的PSO、GWO[8]、SSA、EO[12]以及改进的GWEO对IEEE12、33、118节点配电系统进行功率潮流和电压优化。为公平起见，本文设置上述算法的种群规模和最大迭代次数均为200，最终可得到如表3-表5所示的优化结果。

表3 以IEEE12节点系统为算例系统时的优化目标值 Table 3: Optimal objective values when IEEE 12-bus distribution system is employed as the benchmark network

优化算法	配电网平均损耗(kW)	各节点平均电压偏差(%)
未优化	0.9323	2.7647
PSO	0.4981	2.8512
GWO	0.4741	2.7625
SSA	0.4678	2.7979
EO	0.4491	2.7501
GWEO	0.4491	2.7505

表4 以IEEE33节点系统为算例系统时的优化目标值 Table 4: Optimal objective values when IEEE 33-bus distribution system is employed as the benchmark network

优化算法	配电网平均损耗(kW)	各节点平均电压偏差(%)
未优化	1.6856	2.3894
PSO	1.3940	2.3811
GWO	1.3708	2.2630
SSA	1.3652	2.3425
EO	1.3408	2.2336
GWEO	1.3459	2.2402

表5 以IEEE118节点系统为算例系统时的优化目标值 Table 5: Optimal objective values when IEEE 118-bus distribution system is employed as the benchmark network

优化算法	配电网平均损耗(kW)	各节点平均电压偏差(%)
未优化	4.6584	2.7382
PSO	2.5441	2.9178
GWO	2.6002	2.7254
SSA	2.4829	2.9014
EO	2.3939	2.6773
GWEO	2.3782	2.6861

    根据如表3-表5所示的结果可知，利用改进的GWEO得到的最终优化结果明显优于现有的PSO、SSA和GWO，且接近于EO，因此GWEO在单目标优化场景中体现了良好的收敛性能。

5.3.2 多目标场景下Pareto前沿的分布性

    以配电网平均损耗最小、各节点平均电压偏差最小、分布式光伏设备弃光率最低为优化目标，利用文中现有的多目标灰色狼群算法(Multi-Objective Grey Wolf Optimisation, MOGWO)[9]和多目标平衡优化器算法(Multi-Objective Equilibrium Optimiser, MOEO)[13]与改进的多目标灰狼平衡算法(Multi-Objective Grey Wolf Equilibrium Optimiser, MOGWEO)对IEEE12、33、118节点配电系统进行多目标功率潮流和电压优化。为公平起见，本文设置上述算法的种群规模、最大迭代次数、外部档案的最大规模均为200，并设置栅格数均为10，最终可得到如图7-图9所示的Pareto前沿。

    结合如图7-图9所示的Pareto前沿可知，当进行多目标场景下的优化时，利用MOGWO搜索得到的Pareto前沿分布较为均匀，但在目标空间上的整体分布范围不够广泛；利用MOEO得到的Pareto前沿在目标空间中的分布范围过于狭小，且集中在中间解附近，空间分布不均匀；利用改进MOGWEO得到的Pareto前沿能够在目标空间中广泛且均匀地分布。

图7 以IEEE12节点系统为算例系统时的Pareto前沿 Figure 7: Pareto fronts when IEEE 12-bus distribution system is employed as the benchmark network

图8 以IEEE33节点系统为算例系统时的Pareto前沿 Figure 8: Pareto fronts when IEEE 33-bus distribution system is employed as the benchmark network

图9 以IEEE118节点系统为算例系统时的Pareto前沿 Figure 9: Pareto fronts when IEEE 118-bus distribution system is employed as the benchmark network

6 结论

    本文首先开发了改进的GWEO，该算法在GWO群体更新公式的基础上加入了EO的扰动机制，指定参数$\mu =2$，并设置参数$a$以式(18)所示的形式非线性递减。再以IEEE 12节点、33节点、118节点配电系统为算例，将现有的PSO、GWO、SSA、EO以及改进的GWEO应用到考虑分布式接入的配电网功率潮流和电压优化中。最终对比各种算法的优化结果以体现改进的GWEO在收敛性和Pareto前沿分布性方面的提升作用。

    仿真结果表明，利用GWEO得到的优化结果的收敛性明显优于多数现有的优化算法；而在进行多目标优化时，利用GWEO得到的Pareto前沿能够更为广泛且均匀地分布在目标空间中。因此，利用改进的GWEO得到的优化结果具有优良的收敛性和分布性，可为决策者提供更为优质的功率潮流和电压优化方案，因此更能满足实际配电网功率潮流和电压优化的场景需求。

参考文献

[1] 朱永强，朱凌志．分布式电源接入电网的电能质量[M]. 北京：中国电力出版社，2015．
[2] 谢成旺．多目标群体智能优化算法[M]．北京：北京理工大学出版社，2020．
[3] 温正，孙华克．MATLAB智能算法[M]．北京：清华大学出版社，2017．
[4] 丁蕊．复杂高维多目标优化方法[M]．北京：电子工业出版社，2020．
[5] Deb K．，Pratap A．，Agarwal S．A Fast and Elitist Multiobjective Genetic Algorithm：NSGA-Ⅱ [J]．IEEE Transactions on Evolutionary Computation，2002，6(2)：182-197．
[6] C．A．C．Coello，G．T．Pulido，M．S．Lechuga．Handling multiple objectives with particle swarm optimization [J]．IEEE Transactions on Evolutionary Computation，2004，8(3)：256-279．
[7] Q. Zhang，H．Li．MOEA/D：A Multiobjective Evolutionary Algorithm Based on Decomposition [J]．IEEE Transactions on Evolutionary Computation，2007，11(6)：712-731．
[8] S．Mirjalili，S．M．Mirjalili，A. Lewis．Grey Wolf Optimizer [J]．Advances in Engineering Software，2014，69：46–61．
[9] S．Mirjalili，S．Saremia，S．M．Mirjalili，L．S．Coelhod．Multi-objective grey wolf optimizer: A novel algorithm for multi-criterion optimization [J]．Expert Systems With Applications，2016，47：106-119．
[10] Deb K．An efficient constraint handling method for genetic algorithms [J]．Computer Methods in Applied Mechanics and Engineering，2000，186：311-338．
[11] 李智勇，黄滔，陈少淼，等. 约束优化进化算法综述[J]．软件学报，2017，28(6)：1529-1546．
[12] A．Faramarzi，M．Heidarinejad，B．Stephens，et al．Equilibrium optimizer：A novel optimization algorithm [J]．Knowledge-Based Systems，2020(191)：112-132．
[13] A．B．Mohamed，R．Mohamed，S．Mirjalili，et al．MOEO-EED：A multi-objective equilibrium optimizer with exploration-exploitation dominance strategy [J]．Knowledge-based systems，2021，214：1-26．
[14] M．Farivar，S．H．Low．Branch flow model：relaxations and convexification—Part Ⅰ [J]．IEEE Transactions on Power Systems，2013，28(3)：2554-2564．
[15] M．Farivar，S．H．Low．Branch flow model: relaxations and convexification—Part Ⅱ [J]．IEEE Transactions on Power Systems，2013，28(3)：2565-2572．
[16] 吴小刚，刘宗歧，田立亭，等．基于改进多目标粒子群算法的配电网储能选址定容[J]．电网技术，2014，38(12)：3405-3411．
[17] S．H．Hemati，G．H．Sheisi，S．Karimi．Allocation-based optimal reactive power dispatch considering polynomial load model using improved grey wolf optimizer．Iranian Journal of Science and Technology，Transactions of Electrical Engineering，2021，45(3)，921-944．
[18] 杨蕾，李胜男，黄伟，等．基于平衡优化器的含高比例风光新能源电网无功优化[J]．电力系统及其自动化学报，2021，33(4)：32-39．
[19] D．Das，H．S．Nagi，D．P．Kothari．Novel method for solving radial distribution networks [J]．IEE Proceedings – Generation，Transmission & Distribution，141(4)，291-298，1994．
[20] M．E．Baran，F．F．Wu．Network reconfiguration in distribution systems for loss reduction and load balancing [J]．IEEE Transactions on Power Delivery，1989，4(2)：1401-1407．
[21] D．Zhang，Z．Fu，L．Zhang．An improved TS algorithm for loss-minimum reconfiguration in large-scale distribution systems [J]．Electric power systems research，2007，77(5-6)：685-694．
[22] H．R．E．H．Bouchekara，Y．Latreche，K．Naidu，et al．Comprehensive Review of Radial Distribution Test Systems for Power System Distribution Education and Research [J]．Resource-Efficient Technologies，2019，3：1-12．
[23] 国家电网公司．Q/GDW 1480-2015 分布式电源接入电网技术规定[S]．北京：中国电力科学研究院，2015．

附录

附录1 灰狼平衡算法程序代码

function [Alpha_score,Alpha_pos,Convergence_curve]=GWEO(SearchAgents_no,Max_iter,lb,ub,dim,fobj,fopt)
    
    % GWEO parameters
    V = 1;
    a1 = 2;
    a2 = 1;
    GP = 0.5;

    % initialize alpha, beta, and delta_pos
    Alpha_pos=zeros(1,dim);
    Alpha_score=inf; %change this to -inf for maximization problems

    Beta_pos=zeros(1,dim);
    Beta_score=inf; %change this to -inf for maximization problems

    Delta_pos=zeros(1,dim);
    Delta_score=inf; %change this to -inf for maximization problems

    %Initialize the positions of search agents
    Positions=initialization(SearchAgents_no,dim,ub,lb);

    Convergence_curve=zeros(1,Max_iter);

    l=1;% Loop counter

    % Main loop
    while l<=Max_iter
        for i=1:size(Positions,1)  

           % Return back the search agents that go beyond the boundaries of the search space
            Flag4ub=Positions(i,:)>ub;
            Flag4lb=Positions(i,:)<lb;
            Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;               

            % Calculate objective function for each search agent
            fitness=fobj(Positions(i,:),fopt);

            % Update Alpha, Beta, and Delta
            if fitness<Alpha_score 
                Alpha_score=fitness; % Update alpha
                Alpha_pos=Positions(i,:);
            end

            if fitness>Alpha_score && fitness<Beta_score 
                Beta_score=fitness; % Update beta
                Beta_pos=Positions(i,:);
            end

            if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score 
                Delta_score=fitness; % Update delta
                Delta_pos=Positions(i,:);
            end
        end

        t = (1-l/Max_iter)^(2*a2*l/Max_iter);                             % Eq (9)
        a = 2.*t;

        % Update the Position of search agents including omegas
        for i=1:size(Positions,1)
            % For ALPHA
            Ceq = Alpha_pos;
            % Eq.(3.4) in the paper
            c = 2.*rand(1, dim);
            % Eq.(3.1) in the paper
            D = abs(c.*Ceq-Positions(i,:));
            % Eq.(3.3) in the paper
            A = 2.*a.*rand()-a;
            lambda = rand(1, dim);                                 % lambda in Eq(11)
            r = rand(1, dim);                                      % r in Eq(11)
            F = a1*sign(r-0.5).*(exp(-lambda.*t)-1);                % Eq(11)
            r1 = rand(); r2 = rand();                               % r1 and r2 in Eq(15)
            GCP = 0.5*r1*ones(1, dim)*(r2 >= GP);                  % Eq(15)
            G0 = GCP.*(Ceq-lambda.*Positions(i,:));        % Eq(14)
            G = G0.*F;                                              % Eq(13)
            Position_Alpha = Ceq + ((Positions(i,:)-Ceq).*F...
                + (G./lambda*V).*(1-F)) - A.*abs(D);                            % Eq(16) for ALPHA
            %         Position_Alpha = Ceq + (Positions(i,:)-Ceq).*F...
            %             - A.*abs(D);                            % Eq(16) for ALPHA
            %         Position_Alpha = Ceq - A.*abs(D) + (G./lambda*V).*(1-F); % Eq(16) for ALPHA
            
            % For BETA
            Ceq = Beta_pos;
            % Eq.(3.4) in the paper
            c = 2.*rand(1, dim);
            % Eq.(3.1) in the paper
            D = abs(c.*Ceq-Positions(i,:));
            % Eq.(3.3) in the paper
            A = 2.*a.*rand()-a;
            lambda = rand(1, dim);                                 % lambda in Eq(11)
            r = rand(1, dim);                                      % r in Eq(11)
            F = a1*sign(r-0.5).*(exp(-lambda.*t)-1);                % Eq(11)
            r1 = rand(); r2 = rand();                               % r1 and r2 in Eq(15)
            GCP = 0.5*r1*ones(1, dim)*(r2 >= GP);                  % Eq(15)
            G0 = GCP.*(Ceq-lambda.*Positions(i,:));        % Eq(14)
            G = G0.*F;                                              % Eq(13)
            Position_Beta = Ceq + ((Positions(i,:)-Ceq).*F...
                + (G./lambda*V).*(1-F)) - A.*abs(D);                            % Eq(16) for BETA
            %         Position_Beta = Ceq + (Positions(i,:)-Ceq).*F...
            %             - A.*abs(D);                            % Eq(16) for BETA
            %         Position_Beta = Ceq - A.*abs(D) + (G./lambda*V).*(1-F); % Eq(16) for ALPHA
            
            % For DELTA
            Ceq = Delta_pos;
            % Eq.(3.4) in the paper
            c = 2.*rand(1, dim);
            % Eq.(3.1) in the paper
            D = abs(c.*Ceq-Positions(i,:));
            % Eq.(3.3) in the paper
            A = 2.*a.*rand()-a;
            lambda = rand(1, dim);                                 % lambda in Eq(11)
            r = rand(1, dim);                                      % r in Eq(11)
            F = a1*sign(r-0.5).*(exp(-lambda.*t)-1);                % Eq(11)
            r1 = rand(); r2 = rand();                               % r1 and r2 in Eq(15)
            GCP = 0.5*r1*ones(1, dim)*(r2 >= GP);                  % Eq(15)
            G0 = GCP.*(Ceq-lambda.*Positions(i,:));        % Eq(14)
            G = G0.*F;                                              % Eq(13)
            Position_Delta = Ceq + ((Positions(i,:)-Ceq).*F...
                + (G./lambda*V).*(1-F)) - A.*abs(D);                            % Eq(16) for DELTA
            %         Position_Delta = Ceq + (Positions(i,:)-Ceq).*F...
            %             - A.*abs(D);                            % Eq(16) for DELTA
            %         Position_Delta = Ceq - A.*abs(D) + (G./lambda*V).*(1-F); % Eq(16) for ALPHA
            
            % Global
            Positions(i, :) = (Position_Alpha + Position_Beta + Position_Delta)/3;
%             for j=1:size(Positions,2)     
% 
%                 r1=rand(); % r1 is a random number in [0,1]
%                 r2=rand(); % r2 is a random number in [0,1]
% 
%                 A1=2*a*r1-a; % Equation (3.3)
%                 C1=2*r2; % Equation (3.4)
%                 
%                 D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); % Equation (3.5)-part 1
%                 X1=Alpha_pos(j)-A1*D_alpha+(G./lambda*V).*(1-F); % Equation (3.6)-part 1
% 
%                 r1=rand();
%                 r2=rand();
% 
%                 A2=2*a*r1-a; % Equation (3.3)
%                 C2=2*r2; % Equation (3.4)
% 
%                 D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); % Equation (3.5)-part 2
%                 X2=Beta_pos(j)-A2*D_beta+(G./lambda*V).*(1-F); % Equation (3.6)-part 2       
% 
%                 r1=rand();
%                 r2=rand(); 
% 
%                 A3=2*a*r1-a; % Equation (3.3)
%                 C3=2*r2; % Equation (3.4)
% 
%                 D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); % Equation (3.5)-part 3
%                 X3=Delta_pos(j)-A3*D_delta+(G./lambda*V).*(1-F); % Equation (3.5)-part 3             
% 
%                 Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7)
% 
%             end
        end   
        Convergence_curve(l)=Alpha_score;
        % Display Results
        str1 = strcat('Best estimates: Alpha_pos=');
        disp(str1);
        str2 = strcat(num2str(Alpha_pos));
        str2 = strcat(str2,'  iteration='); str2=strcat(str2,num2str(l));
        disp(str2);
        l=l+1;
    end
end

附录2 基于MATPOWER的IEEE标准节点系统仿真建模代码

附录2.1 IEEE 12节点系统仿真建模代码

function mpc = case12da
%CASE12DA  Power flow data for 12 bus distribution system from Das, et al
%   Please see CASEFORMAT for details on the case file format.
%
%   Data from ...
%       D. Das, H.S. Nagi, D.P. Kothari, "Novel method for solving radial
%       distribution networks", IEE Proc. C, Vol. 141, No. 4, pp. 291-298, 1994.

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 1;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [ %% (Pd and Qd are specified in kW & kVAr here, converted to MW & MVAr below)
	1	3	0	0	0	0	1	1.05	0	10	1	1.1	0.9;
	2	1	60	60	0	0	1	1       0	10	1	1.1	0.9;
	3	1	40	30	0	0	1	1       0	10	1	1.1	0.9;
	4	1	55	55	0	0	1	1       0	10	1	1.1	0.9;
	5	1	30	30	0	0	1	1       0	10	1	1.1	0.9;
	6	1	20	15	0	0	1	1       0	10	1	1.1	0.9;
	7	1	55	55	0	0	1	1       0	10	1	1.1	0.9;
	8	1	45	45	0	0	1	1       0	10	1	1.1	0.9;
	9	1	40	40	0	0	1	1       0	10	1	1.1	0.9;
	10	1	35	30	0	0	1	1       0	10	1	1.1	0.9;
	11	1	40	30	0	0	1	1       0	10	1	1.1	0.9;
	12	1	15	15	0	0	1	1       0	10	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
	1	0	0	10	-10	1.05	100	1	10	0	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [  %% (r and x specified in ohms here, converted to p.u. below)
	1	2	1.093	0.455	0	0	0	0	0	0	1	-360	360;
	2	3	1.184	0.494	0	0	0	0	0	0	1	-360	360;
	3	4	2.095	0.873	0	0	0	0	0	0	1	-360	360;
	4	5	3.188	1.329	0	0	0	0	0	0	1	-360	360;
	5	6	1.093	0.455	0	0	0	0	0	0	1	-360	360;
	6	7	1.002	0.417	0	0	0	0	0	0	1	-360	360;
	7	8	4.403	1.215	0	0	0	0	0	0	1	-360	360;
	8	9	5.642	1.597	0	0	0	0	0	0	1	-360	360;
	9	10	2.89	0.818	0	0	0	0	0	0	1	-360	360;
	10	11	1.514	0.428	0	0	0	0	0	0	1	-360	360;
	11	12	1.238	0.351	0	0	0	0	0	0	1	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0	20	0;
];


%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3;      %% in Volts
Sbase = mpc.baseMVA * 1e6;              %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);

%% convert loads from kW to MW
mpc.bus(:, [PD, QD]) = mpc.bus(:, [PD, QD]) / 1e3;

附录2.2 IEEE 33节点系统仿真建模代码

function mpc = case33
mpc.version = '2';
mpc.baseMVA = 100;
%%bus data
%   bus_i type    Pd  Qd  Gs  Bs  area    Vm  Va  baseKV  zone    Vmax    Vmin
mpc.bus = [
    1	3	0	    0	    0	0	1	1.05	0	12.66	1	1.1	0.9;
%     1	3	0	    0	    0	0	1.05	1	0	12.66	1	1.05	0.95;
	2	1	0.100   0.060	0	0	1	1	0	12.66	1	1.05	0.95;
	3	1	0.090	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
	4	1	0.120	0.080	0	0	1	1	0	12.66	1	1.05	0.95;
	5	1	0.060	0.030	0	0	1	1	0	12.66	1	1.05	0.95;
	6	1	0.060	0.020	0	0	1	1	0	12.66	1	1.05	0.95;
	7	1	0.200   0.100	0	0	1	1	0	12.66	1	1.05	0.95;
	8	1   0.200   0.100   0	0	1	1	0	12.66	1	1.05	0.95;
	9	1	0.060	0.020	0	0	1	1	0	12.66	1	1.05	0.95;
	10	1	0.060	0.020	0	0	1	1	0	12.66	1	1.05	0.95;
	11	1	0.045	0.030	0	0	1	1	0	12.66	1	1.05	0.95;
	12	1	0.060	0.035	0	0	1	1	0	12.66	1	1.05	0.95;
	13	1	0.060	0.035	0	0	1	1	0	12.66	1	1.05	0.95;
	14	1	0.120	0.080	0	0	1	1	0	12.66	1	1.05	0.95;
	15	1	0.060	0.010	0	0	1	1	0	12.66	1	1.05	0.95;
	16	1	0.060	0.020	0	0	1	1	0	12.66	1	1.05	0.95;
	17	1	0.060	0.020	0	0	1	1	0	12.66	1	1.05	0.95;
	18	1	0.090	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
	19	1	0.090	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
	20	1	0.090	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
	21	1	0.090	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
	22	1	0.090	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
	23	1	0.090	0.050	0	0	1	1	0	12.66	1	1.05	0.95;
	24	1	0.420	0.200	0	0	1	1	0	12.66	1	1.05	0.95;
	25	1	0.420	0.200	0	0	1	1	0	12.66	1	1.05	0.95;
	26	1	0.060	0.025	0	0	1	1	0	12.66	1	1.05	0.95;
	27	1	0.060	0.025	0	0	1	1	0	12.66	1	1.05	0.95;
	28	1	0.060	0.020	0	0	1	1	0	12.66	1	1.05	0.95;
	29	1	0.120	0.070	0	0	1	1	0	12.66	1	1.05	0.95;
	30	1	0.200	0.100	0	0	1	1	0	12.66	1	1.05	0.95;
    31	1	0.150	0.070	0	0	1	1	0	12.66	1	1.05	0.95;
    32	1	0.210	0.100	0	0	1	1	0	12.66	1	1.05	0.95;
    33	1	0.060	0.040	0	0	1	1	0	12.66	1	1.05	0.95;
];
%%generator data
%   bus    Pg  Qg  Qmax    Qmin    Vg  mBase   status  Pmax    Pmin
mpc.gen = [
    1    0       0      60 -60    1.05   100   1   100   0  
        
];
%%branch data
%   fbus   tbus    r   x   b   rateA   rateB   rateC   ratio   angle
%   status
mpc.branch = [  %% (r and x specified in ohms here, converted to p.u. below)
    1   2	0.0922	0.0470	0	130	130	130 0	0	1;
	2	3	0.4930	0.2512	0	130	130	130	0	0	1;
	3	4	0.3661	0.1864	0	130	130	130	0	0	1;
	4	5	0.3811	0.1941	0	130	130	130	0	0	1;
	5	6	0.8190	0.7070	0	130	130	130	0	0	1;
	6	7	0.1872	0.6188	0	130	130	130	0	0	1;
	7	8	0.7115	0.2351	0	130	130	130	0	0	1;
	8	9	1.0299	0.7400	0	130	130	130	0	0	1;
	9	10	1.0440	0.7400	0	130	130	130	0	0	1;
	10	11	0.1967	0.0651	0	130	130	130	0	0	1;
	11	12	0.3744	0.1298	0	65	65	65	0	0	1;
	12	13	1.4680	1.1549	0	32	32	32	0	0	1;
	13	14	0.5416	0.7129	0	65	65	65	0	0	1;
	14	15	0.5909	0.5260	0	65	65	65	0	0	1;
	15	16	0.7462	0.5449	0	65	65	65	0	0	1;
	16	17	1.2889	1.7210	0	65	65	65	0	0	1;
	17	18	0.7320	0.5739	0	32	32	32	0	0	1;
	2	19	0.1640	0.1565	0	32	32	32	0	0	1;
	19	20	1.5042	1.3555	0	32	32	32	0	0	1;
	20	21	0.4095	0.4784	0	16	16	16	0	0	1;
	21	22	0.7089	0.9373	0	16	16	16	0	0	1;
	3	23	0.4512	0.3084	0	16	16	16	0	0	1;
	23	24	0.8980	0.7091	0	16	16	16	0	0	1;
	24	25	0.8959	0.7071	0	32	32	32	0	0	1;
	6	26	0.2031	0.1034	0	32	32	32	0	0	1;
	26	27	0.2842	0.1447	0	32	32	32	0	0	1;
	27	28	1.0589	0.9338	0	32	32	32	0	0	1;
	28	29	0.8043	0.7006	0	32	32	32	0	0	1;
	29  30	0.5074	0.2585	0	32	32	32	0	0	1;
	30	31	0.9745	0.9629	0	16	16	16	0	0	1;
	31	32	0.3105	0.3619	0	16	16	16	0	0	1;
	32	33	0.3411	0.5302	0	16	16	16	0	0	1;
];

%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3;      %% in Volts
Sbase = mpc.baseMVA * 1e6;              %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);
return;

附录2.3 IEEE 118节点系统仿真建模代码

function mpc = case118zh
%CASE118ZH  Power flow data for 118 bus distribution system from Zhang, et al
%   Please see CASEFORMAT for details on the case file format.
%
%   Data from ...
%       Zhang D, Fu Z, Zhang L (2007) An improved TS algorithm for
%       loss-minimum reconfiguration in large-scale distribution systems.
%       77:685-694. doi: 10.1016/j.epsr.2006.06.005
%       URL: https://doi.org/10.1016/j.epsr.2006.06.005

%% MATPOWER Case Format : Version 2
mpc.version = '2';

%%-----  Power Flow Data  -----%%
%% system MVA base
mpc.baseMVA = 10;

%% bus data
%	bus_i	type	Pd	Qd	Gs	Bs	area	Vm	Va	baseKV	zone	Vmax	Vmin
mpc.bus = [ %% (Pd and Qd are specified in kW & kVAr here, converted to MW & MVAr below)
% 	1	3	0	0	0	0	1	1	0	11	1	1.1	0.9;
    1	3	0	0	0	0	1	1.05	0	11	1	1.1	0.9;
	2	1	133.84	101.14	0	0	1	1	0	11	1	1.1	0.9;
	3	1	16.214	11.292	0	0	1	1	0	11	1	1.1	0.9;
	4	1	34.315	21.845	0	0	1	1	0	11	1	1.1	0.9;
	5	1	73.016	63.602	0	0	1	1	0	11	1	1.1	0.9;
	6	1	144.2	68.604	0	0	1	1	0	11	1	1.1	0.9;
	7	1	104.47	61.725	0	0	1	1	0	11	1	1.1	0.9;
	8	1	28.547	11.503	0	0	1	1	0	11	1	1.1	0.9;
	9	1	87.56	51.073	0	0	1	1	0	11	1	1.1	0.9;
	10	1	198.2	106.77	0	0	1	1	0	11	1	1.1	0.9;
	11	1	146.8	75.995	0	0	1	1	0	11	1	1.1	0.9;
	12	1	26.04	18.687	0	0	1	1	0	11	1	1.1	0.9;
	13	1	52.1	23.22	0	0	1	1	0	11	1	1.1	0.9;
	14	1	141.9	117.5	0	0	1	1	0	11	1	1.1	0.9;
	15	1	21.87	28.79	0	0	1	1	0	11	1	1.1	0.9;
	16	1	33.37	26.45	0	0	1	1	0	11	1	1.1	0.9;
	17	1	32.43	25.23	0	0	1	1	0	11	1	1.1	0.9;
	18	1	20.234	11.906	0	0	1	1	0	11	1	1.1	0.9;
	19	1	156.94	78.523	0	0	1	1	0	11	1	1.1	0.9;
	20	1	546.29	351.4	0	0	1	1	0	11	1	1.1	0.9;
	21	1	180.31	164.2	0	0	1	1	0	11	1	1.1	0.9;
	22	1	93.167	54.594	0	0	1	1	0	11	1	1.1	0.9;
	23	1	85.18	39.65	0	0	1	1	0	11	1	1.1	0.9;
	24	1	168.1	95.178	0	0	1	1	0	11	1	1.1	0.9;
	25	1	125.11	150.22	0	0	1	1	0	11	1	1.1	0.9;
	26	1	16.03	24.62	0	0	1	1	0	11	1	1.1	0.9;
	27	1	26.03	24.62	0	0	1	1	0	11	1	1.1	0.9;
	28	1	594.56	522.62	0	0	1	1	0	11	1	1.1	0.9;
	29	1	120.62	59.117	0	0	1	1	0	11	1	1.1	0.9;
	30	1	102.38	99.554	0	0	1	1	0	11	1	1.1	0.9;
	31	1	513.4	318.5	0	0	1	1	0	11	1	1.1	0.9;
	32	1	475.25	456.14	0	0	1	1	0	11	1	1.1	0.9;
	33	1	151.43	136.79	0	0	1	1	0	11	1	1.1	0.9;
	34	1	205.38	83.302	0	0	1	1	0	11	1	1.1	0.9;
	35	1	131.6	93.082	0	0	1	1	0	11	1	1.1	0.9;
	36	1	448.4	369.79	0	0	1	1	0	11	1	1.1	0.9;
	37	1	440.52	321.64	0	0	1	1	0	11	1	1.1	0.9;
	38	1	112.54	55.134	0	0	1	1	0	11	1	1.1	0.9;
	39	1	53.963	38.998	0	0	1	1	0	11	1	1.1	0.9;
	40	1	393.05	342.6	0	0	1	1	0	11	1	1.1	0.9;
	41	1	326.74	278.56	0	0	1	1	0	11	1	1.1	0.9;
	42	1	536.26	240.24	0	0	1	1	0	11	1	1.1	0.9;
	43	1	76.247	66.562	0	0	1	1	0	11	1	1.1	0.9;
	44	1	53.52	39.76	0	0	1	1	0	11	1	1.1	0.9;
	45	1	40.328	31.964	0	0	1	1	0	11	1	1.1	0.9;
	46	1	39.653	20.758	0	0	1	1	0	11	1	1.1	0.9;
	47	1	66.195	42.361	0	0	1	1	0	11	1	1.1	0.9;
	48	1	73.904	51.653	0	0	1	1	0	11	1	1.1	0.9;
	49	1	114.77	57.965	0	0	1	1	0	11	1	1.1	0.9;
	50	1	918.37	1205.1	0	0	1	1	0	11	1	1.1	0.9;
	51	1	210.3	146.66	0	0	1	1	0	11	1	1.1	0.9;
	52	1	66.68	56.608	0	0	1	1	0	11	1	1.1	0.9;
	53	1	42.207	40.184	0	0	1	1	0	11	1	1.1	0.9;
	54	1	433.74	283.41	0	0	1	1	0	11	1	1.1	0.9;
	55	1	62.1	26.86	0	0	1	1	0	11	1	1.1	0.9;
	56	1	92.46	88.38	0	0	1	1	0	11	1	1.1	0.9;
	57	1	85.188	55.436	0	0	1	1	0	11	1	1.1	0.9;
	58	1	345.3	332.4	0	0	1	1	0	11	1	1.1	0.9;
	59	1	22.5	16.83	0	0	1	1	0	11	1	1.1	0.9;
	60	1	80.551	49.156	0	0	1	1	0	11	1	1.1	0.9;
	61	1	95.86	90.758	0	0	1	1	0	11	1	1.1	0.9;
	62	1	62.92	47.7	0	0	1	1	0	11	1	1.1	0.9;
	63	1	478.8	463.74	0	0	1	1	0	11	1	1.1	0.9;
	64	1	120.94	52.006	0	0	1	1	0	11	1	1.1	0.9;
	65	1	139.11	100.34	0	0	1	1	0	11	1	1.1	0.9;
	66	1	391.78	193.5	0	0	1	1	0	11	1	1.1	0.9;
	67	1	27.741	26.713	0	0	1	1	0	11	1	1.1	0.9;
	68	1	52.814	25.257	0	0	1	1	0	11	1	1.1	0.9;
	69	1	66.89	38.713	0	0	1	1	0	11	1	1.1	0.9;
	70	1	467.5	395.14	0	0	1	1	0	11	1	1.1	0.9;
	71	1	594.85	239.74	0	0	1	1	0	11	1	1.1	0.9;
	72	1	132.5	84.363	0	0	1	1	0	11	1	1.1	0.9;
	73	1	52.699	22.482	0	0	1	1	0	11	1	1.1	0.9;
	74	1	869.79	614.775	0	0	1	1	0	11	1	1.1	0.9;
	75	1	31.349	29.817	0	0	1	1	0	11	1	1.1	0.9;
	76	1	192.39	122.43	0	0	1	1	0	11	1	1.1	0.9;
	77	1	65.75	45.37	0	0	1	1	0	11	1	1.1	0.9;
	78	1	238.15	223.22	0	0	1	1	0	11	1	1.1	0.9;
	79	1	294.55	162.47	0	0	1	1	0	11	1	1.1	0.9;
	80	1	485.57	437.92	0	0	1	1	0	11	1	1.1	0.9;
	81	1	243.53	183.03	0	0	1	1	0	11	1	1.1	0.9;
	82	1	243.53	183.03	0	0	1	1	0	11	1	1.1	0.9;
	83	1	134.25	119.29	0	0	1	1	0	11	1	1.1	0.9;
	84	1	22.71	27.96	0	0	1	1	0	11	1	1.1	0.9;
	85	1	49.513	26.515	0	0	1	1	0	11	1	1.1	0.9;
	86	1	383.78	257.16	0	0	1	1	0	11	1	1.1	0.9;
	87	1	49.64	20.6	0	0	1	1	0	11	1	1.1	0.9;
	88	1	22.473	11.806	0	0	1	1	0	11	1	1.1	0.9;
	89	1	62.93	42.96	0	0	1	1	0	11	1	1.1	0.9;
	90	1	30.67	34.93	0	0	1	1	0	11	1	1.1	0.9;
	91	1	62.53	66.79	0	0	1	1	0	11	1	1.1	0.9;
	92	1	114.57	81.748	0	0	1	1	0	11	1	1.1	0.9;
	93	1	81.292	66.526	0	0	1	1	0	11	1	1.1	0.9;
	94	1	31.733	15.96	0	0	1	1	0	11	1	1.1	0.9;
	95	1	33.32	60.48	0	0	1	1	0	11	1	1.1	0.9;
	96	1	531.28	224.85	0	0	1	1	0	11	1	1.1	0.9;
	97	1	507.03	367.42	0	0	1	1	0	11	1	1.1	0.9;
	98	1	26.39	11.7	0	0	1	1	0	11	1	1.1	0.9;
	99	1	45.99	30.392	0	0	1	1	0	11	1	1.1	0.9;
	100	1	100.66	47.572	0	0	1	1	0	11	1	1.1	0.9;
	101	1	456.48	350.3	0	0	1	1	0	11	1	1.1	0.9;
	102	1	522.56	449.29	0	0	1	1	0	11	1	1.1	0.9;
	103	1	408.43	168.46	0	0	1	1	0	11	1	1.1	0.9;
	104	1	141.48	134.25	0	0	1	1	0	11	1	1.1	0.9;
	105	1	104.43	66.024	0	0	1	1	0	11	1	1.1	0.9;
	106	1	96.793	83.647	0	0	1	1	0	11	1	1.1	0.9;
	107	1	493.92	419.34	0	0	1	1	0	11	1	1.1	0.9;
	108	1	225.38	135.88	0	0	1	1	0	11	1	1.1	0.9;
	109	1	509.21	387.21	0	0	1	1	0	11	1	1.1	0.9;
	110	1	188.5	173.46	0	0	1	1	0	11	1	1.1	0.9;
	111	1	918.03	898.55	0	0	1	1	0	11	1	1.1	0.9;
	112	1	305.08	215.37	0	0	1	1	0	11	1	1.1	0.9;
	113	1	54.38	40.97	0	0	1	1	0	11	1	1.1	0.9;
	114	1	211.14	192.9	0	0	1	1	0	11	1	1.1	0.9;
	115	1	67.009	53.336	0	0	1	1	0	11	1	1.1	0.9;
	116	1	162.07	90.321	0	0	1	1	0	11	1	1.1	0.9;
	117	1	48.785	29.156	0	0	1	1	0	11	1	1.1	0.9;
	118	1	33.9	18.98	0	0	1	1	0	11	1	1.1	0.9;
];

%% generator data
%	bus	Pg	Qg	Qmax	Qmin	Vg	mBase	status	Pmax	Pmin	Pc1	Pc2	Qc1min	Qc1max	Qc2min	Qc2max	ramp_agc	ramp_10	ramp_30	ramp_q	apf
mpc.gen = [
% 	1	0	0	10	-10	1	100	1	10	0	0	0	0	0	0	0	0	0	0	0	0;
    1	0	0	10	-10	1.05	100	1	10	0	0	0	0	0	0	0	0	0	0	0	0;
];

%% branch data
%	fbus	tbus	r	x	b	rateA	rateB	rateC	ratio	angle	status	angmin	angmax
mpc.branch = [  %% (r and x specified in ohms here, converted to p.u. below)
	1	2	0.036	0.01296	0	0	0	0	0	0	1	-360	360;
	2	3	0.033	0.01188	0	0	0	0	0	0	1	-360	360;
	2	4	0.045	0.0162	0	0	0	0	0	0	1	-360	360;
	4	5	0.015	0.054	0	0	0	0	0	0	1	-360	360;
	5	6	0.015	0.054	0	0	0	0	0	0	1	-360	360;
	6	7	0.015	0.0125	0	0	0	0	0	0	1	-360	360;
	7	8	0.018	0.014	0	0	0	0	0	0	1	-360	360;
	8	9	0.021	0.063	0	0	0	0	0	0	1	-360	360;
	2	10	0.166	0.1344	0	0	0	0	0	0	1	-360	360;
	10	11	0.112	0.0789	0	0	0	0	0	0	1	-360	360;
	11	12	0.187	0.313	0	0	0	0	0	0	1	-360	360;
	12	13	0.142	0.1512	0	0	0	0	0	0	1	-360	360;
	13	14	0.18	0.118	0	0	0	0	0	0	1	-360	360;
	14	15	0.15	0.045	0	0	0	0	0	0	1	-360	360;
	15	16	0.16	0.18	0	0	0	0	0	0	1	-360	360;
	16	17	0.157	0.171	0	0	0	0	0	0	1	-360	360;
	11	18	0.218	0.285	0	0	0	0	0	0	1	-360	360;
	18	19	0.118	0.185	0	0	0	0	0	0	1	-360	360;
	19	20	0.16	0.196	0	0	0	0	0	0	1	-360	360;
	20	21	0.12	0.189	0	0	0	0	0	0	1	-360	360;
	21	22	0.12	0.0789	0	0	0	0	0	0	1	-360	360;
	22	23	1.41	0.723	0	0	0	0	0	0	1	-360	360;
	23	24	0.293	0.1348	0	0	0	0	0	0	1	-360	360;
	24	25	0.133	0.104	0	0	0	0	0	0	1	-360	360;
	25	26	0.178	0.134	0	0	0	0	0	0	1	-360	360;
	26	27	0.178	0.134	0	0	0	0	0	0	1	-360	360;
	4	28	0.015	0.0296	0	0	0	0	0	0	1	-360	360;
	28	29	0.012	0.0276	0	0	0	0	0	0	1	-360	360;
	29	30	0.12	0.2766	0	0	0	0	0	0	1	-360	360;
	30	31	0.21	0.243	0	0	0	0	0	0	1	-360	360;
	31	32	0.12	0.054	0	0	0	0	0	0	1	-360	360;
	32	33	0.178	0.234	0	0	0	0	0	0	1	-360	360;
	33	34	0.178	0.234	0	0	0	0	0	0	1	-360	360;
	34	35	0.154	0.162	0	0	0	0	0	0	1	-360	360;
	30	36	0.187	0.261	0	0	0	0	0	0	1	-360	360;
	36	37	0.133	0.099	0	0	0	0	0	0	1	-360	360;
	29	38	0.33	0.194	0	0	0	0	0	0	1	-360	360;
	38	39	0.31	0.194	0	0	0	0	0	0	1	-360	360;
	39	40	0.13	0.194	0	0	0	0	0	0	1	-360	360;
	40	41	0.28	0.15	0	0	0	0	0	0	1	-360	360;
	41	42	1.18	0.85	0	0	0	0	0	0	1	-360	360;
	42	43	0.42	0.2436	0	0	0	0	0	0	1	-360	360;
	43	44	0.27	0.0972	0	0	0	0	0	0	1	-360	360;
	44	45	0.339	0.1221	0	0	0	0	0	0	1	-360	360;
	45	46	0.27	0.1779	0	0	0	0	0	0	1	-360	360;
	35	47	0.21	0.1383	0	0	0	0	0	0	1	-360	360;
	47	48	0.12	0.0789	0	0	0	0	0	0	1	-360	360;
	48	49	0.15	0.0987	0	0	0	0	0	0	1	-360	360;
	49	50	0.15	0.0987	0	0	0	0	0	0	1	-360	360;
	50	51	0.24	0.1581	0	0	0	0	0	0	1	-360	360;
	51	52	0.12	0.0789	0	0	0	0	0	0	1	-360	360;
	52	53	0.405	0.1458	0	0	0	0	0	0	1	-360	360;
	53	54	0.405	0.1458	0	0	0	0	0	0	1	-360	360;
	29	55	0.391	0.141	0	0	0	0	0	0	1	-360	360;
	55	56	0.406	0.1461	0	0	0	0	0	0	1	-360	360;
	56	57	0.406	0.1461	0	0	0	0	0	0	1	-360	360;
	57	58	0.706	0.5461	0	0	0	0	0	0	1	-360	360;
	58	59	0.338	0.1218	0	0	0	0	0	0	1	-360	360;
	59	60	0.338	0.1218	0	0	0	0	0	0	1	-360	360;
	60	61	0.207	0.0747	0	0	0	0	0	0	1	-360	360;
	61	62	0.247	0.8922	0	0	0	0	0	0	1	-360	360;
	1	63	0.028	0.0418	0	0	0	0	0	0	1	-360	360;
	63	64	0.117	0.2016	0	0	0	0	0	0	1	-360	360;
	64	65	0.255	0.0918	0	0	0	0	0	0	1	-360	360;
	65	66	0.21	0.0759	0	0	0	0	0	0	1	-360	360;
	66	67	0.383	0.138	0	0	0	0	0	0	1	-360	360;
	67	68	0.504	0.3303	0	0	0	0	0	0	1	-360	360;
	68	69	0.406	0.1461	0	0	0	0	0	0	1	-360	360;
	69	70	0.962	0.761	0	0	0	0	0	0	1	-360	360;
	70	71	0.165	0.06	0	0	0	0	0	0	1	-360	360;
	71	72	0.303	0.1092	0	0	0	0	0	0	1	-360	360;
	72	73	0.303	0.1092	0	0	0	0	0	0	1	-360	360;
	73	74	0.206	0.144	0	0	0	0	0	0	1	-360	360;
	74	75	0.233	0.084	0	0	0	0	0	0	1	-360	360;
	75	76	0.591	0.1773	0	0	0	0	0	0	1	-360	360;
	76	77	0.126	0.0453	0	0	0	0	0	0	1	-360	360;
	64	78	0.559	0.3687	0	0	0	0	0	0	1	-360	360;
	78	79	0.186	0.1227	0	0	0	0	0	0	1	-360	360;
	79	80	0.186	0.1227	0	0	0	0	0	0	1	-360	360;
	80	81	0.26	0.139	0	0	0	0	0	0	1	-360	360;
	81	82	0.154	0.148	0	0	0	0	0	0	1	-360	360;
	82	83	0.23	0.128	0	0	0	0	0	0	1	-360	360;
	83	84	0.252	0.106	0	0	0	0	0	0	1	-360	360;
	84	85	0.18	0.148	0	0	0	0	0	0	1	-360	360;
	79	86	0.16	0.182	0	0	0	0	0	0	1	-360	360;
	86	87	0.2	0.23	0	0	0	0	0	0	1	-360	360;
	87	88	0.16	0.393	0	0	0	0	0	0	1	-360	360;
	65	89	0.669	0.2412	0	0	0	0	0	0	1	-360	360;
	89	90	0.266	0.1227	0	0	0	0	0	0	1	-360	360;
	90	91	0.266	0.1227	0	0	0	0	0	0	1	-360	360;
	91	92	0.266	0.1227	0	0	0	0	0	0	1	-360	360;
	92	93	0.266	0.1227	0	0	0	0	0	0	1	-360	360;
	93	94	0.233	0.115	0	0	0	0	0	0	1	-360	360;
	94	95	0.496	0.138	0	0	0	0	0	0	1	-360	360;
	91	96	0.196	0.18	0	0	0	0	0	0	1	-360	360;
	96	97	0.196	0.18	0	0	0	0	0	0	1	-360	360;
	97	98	0.1866	0.122	0	0	0	0	0	0	1	-360	360;
	98	99	0.0746	0.318	0	0	0	0	0	0	1	-360	360;
	1	100	0.0625	0.0265	0	0	0	0	0	0	1	-360	360;
	100	101	0.1501	0.234	0	0	0	0	0	0	1	-360	360;
	101	102	0.1347	0.0888	0	0	0	0	0	0	1	-360	360;
	102	103	0.2307	0.1203	0	0	0	0	0	0	1	-360	360;
	103	104	0.447	0.1608	0	0	0	0	0	0	1	-360	360;
	104	105	0.1632	0.0588	0	0	0	0	0	0	1	-360	360;
	105	106	0.33	0.099	0	0	0	0	0	0	1	-360	360;
	106	107	0.156	0.0561	0	0	0	0	0	0	1	-360	360;
	107	108	0.3819	0.1374	0	0	0	0	0	0	1	-360	360;
	108	109	0.1626	0.0585	0	0	0	0	0	0	1	-360	360;
	109	110	0.3819	0.1374	0	0	0	0	0	0	1	-360	360;
	110	111	0.2445	0.0879	0	0	0	0	0	0	1	-360	360;
	110	112	0.2088	0.0753	0	0	0	0	0	0	1	-360	360;
	112	113	0.2301	0.0828	0	0	0	0	0	0	1	-360	360;
	100	114	0.6102	0.2196	0	0	0	0	0	0	1	-360	360;
	114	115	0.1866	0.127	0	0	0	0	0	0	1	-360	360;
	115	116	0.3732	0.246	0	0	0	0	0	0	1	-360	360;
	116	117	0.405	0.367	0	0	0	0	0	0	1	-360	360;
	117	118	0.489	0.438	0	0	0	0	0	0	1	-360	360;
	46	27	0.5258	0.2925	0	0	0	0	0	0	0	-360	360;
	17	27	0.5258	0.2916	0	0	0	0	0	0	0	-360	360;
	8	24	0.4272	0.1539	0	0	0	0	0	0	0	-360	360;
	54	43	0.48	0.1728	0	0	0	0	0	0	0	-360	360;
	62	49	0.36	0.1296	0	0	0	0	0	0	0	-360	360;
	37	62	0.57	0.572	0	0	0	0	0	0	0	-360	360;
	9	40	0.53	0.3348	0	0	0	0	0	0	0	-360	360;
	58	96	0.3957	0.1425	0	0	0	0	0	0	0	-360	360;
	73	91	0.68	0.648	0	0	0	0	0	0	0	-360	360;
	88	75	0.4062	0.1464	0	0	0	0	0	0	0	-360	360;
	99	77	0.4626	0.1674	0	0	0	0	0	0	0	-360	360;
	108	83	0.651	0.234	0	0	0	0	0	0	0	-360	360;
	105	86	0.8125	0.2925	0	0	0	0	0	0	0	-360	360;
	110	118	0.7089	0.2553	0	0	0	0	0	0	0	-360	360;
	25	35	0.5	0.5	0	0	0	0	0	0	0	-360	360;
];

%%-----  OPF Data  -----%%
%% generator cost data
%	1	startup	shutdown	n	x1	y1	...	xn	yn
%	2	startup	shutdown	n	c(n-1)	...	c0
mpc.gencost = [
	2	0	0	3	0	20	0;
];


%% convert branch impedances from Ohms to p.u.
[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
Vbase = mpc.bus(1, BASE_KV) * 1e3;      %% in Volts
Sbase = mpc.baseMVA * 1e6;              %% in VA
mpc.branch(:, [BR_R BR_X]) = mpc.branch(:, [BR_R BR_X]) / (Vbase^2 / Sbase);

%% convert loads from kW to MW
mpc.bus(:, [PD, QD]) = mpc.bus(:, [PD, QD]) / 1e3;

附录3 完整版代码

附录3.1 不同优化算法的无功优化效果对比

    1. 几种单目标优化算法的无功优化效果对比
    2. 基于改进灰狼平衡算法的配电网功率潮流和电压优化

附录3.2 无功补偿策略探系列

    1. 无功补偿策略探（一）：基于二进制电容器分组的无功补偿优化策略（上）
    2. 无功补偿策略探（二）：基于二进制电容器分组的无功补偿优化策略（下）
    3. 无功补偿策略探（三）：基于任意电容器分组比例的无功补偿优化策略
    4. 无功补偿策略探（四）：基于任意分组比例的电容器组和电抗器组的无功补偿优化策略
    5. 无功补偿策略探（五）：电容器组和电抗器组混合安装的无功补偿优化策略

附录4 相关已发表文献

    1. 考虑分布式光伏接入的配电网多目标潮流电压优化方法（发明专利）
    2. Multiobjective Optimal Power Flow for Distribution Networks Utilizing a Novel Heuristic Algorithm—Grey Wolf Equilibrium Optimizer（IEEE Systems Journal期刊论文）
    3. Multi-Objective Voltage and Reactive Power Coordinated Control Strategy for Distribution Networks Utilizing Gaining-Sharing Knowledge Based Algorithm（CIEEC会议论文）

附录5 博主常用联系方式

    E-mail：williamliu19980525@outlook.com
    QQ：2532213077
    微信号：wxid_lrgrjky8zplf22