基于知识共享的算法在配电网多目标电压无功协调控制优化中的应用

版权说明

    本博客翻译自会议论文Jinjin Ding, Xunting Wang, Bin Xu, Mingxing Zhu, Wei Liu*. Multi-Objective Voltage and Reactive Power Coordinated Control Strategy for Distribution Networks Utilizing Gaining-Sharing Knowledge Based Algorithm [C]// 2023 IEEE 6th International Electrical and Energy Conference, Hefei, China: IEEE.
    (*指代通讯作者)

摘要

    现有的启发式优化算法在处理多目标电压和无功功率协调控制优化问题时存在容易停滞于局部最优、Pareto前沿多样性差等缺点。为了解决上述不足，本文采用了一种新开发的基于共享知识的算法(GSK)来处理多目标电压和无功功率协调控制优化问题。然后，以系统损耗最小、平均电压偏差最小和弃电率最小为优化目标，以改进的IEEE 33节点配电系统作为算例系统。将灰狼优化(GWO)和平衡优化器(EO)作为比较来验证GSK在收敛性和多样性方面的改进。结果表明，GSK可得到更多样化的配电网多目标电压和无功功率协调控制方案，进而可以更好地应用于配电网多目标电压和无功功率协调控制优化的实际场景。

关键词

    多目标电压和无功功率协调控制优化、配电网、分布式光伏(PV)、基于知识共享的算法(GSK)

I. 引言

    为达到“碳达峰、碳中和”目标，推动分布式光伏进入配电网是一项重要举措。对发展以可再生能源为主体的新型电力系统，实现“双碳”目标和乡村振兴战略具有重要意义。然而，大量分布式光伏接入配电网，造成有功潮流波动，进而引发电能损耗增加、末端节点过压、功率倒送、三相不平衡等一系列电能质量问题，对配电系统的安全、稳定和经济运行造成不利影响。因此，配电网调度领域亟需关注分布式光伏接入的多目标电压和无功功率协调控制优化问题。
    在分布式光伏接入配电网的情况下，一般需要综合系统损耗最小、平均电压偏差最小、分布式光伏弃光率最小等多个优化目标来制定配电网多目标电压和无功功率协调控制方案。然而，各种优化目标之间的权衡导致多目标电压和无功功率协调控制优化模型的优化结果不唯一，这也增加了解决问题的难度和复杂性，因此不可避免地要讨论如何解决多目标电压和无功功率协调控制优化上的问题。目前，应用于多目标电压和无功功率协调控制优化的现有算法可分为两大类：直接搜索算法和启发式优化算法。在启发式优化算法中，随机生成一定规模的初始种群，再计算种群中每个个体的适应度值。然后根据启发式优化算法的更新规则，向当前最优适应度方向迭代更新种群，直至收敛到全局最优解。与直接搜索算法相比，启发式优化算法具有以下优点：迭代规则简单，不依赖于初值，不需要目标函数的连续性和可导性以及约束[5] ，更适合求解复杂的问题具有多个峰值和非线性的优化问题。此外，结合Pareto优势理论的启发式算法也被广泛应用于高维度、多目标的复杂应用场景。
    经过几十年的研究，以遗传算法(GA)[6] 、进化算法(EA)[7] 、粒子群算法(PSO)[8] 、免疫算法(IA)[9] 、蚁群算法(ACO)[10] 、灰狼优化算法(GWO)[11] 、樽海鞘群优化算法(SSO)[12]和平衡优化器算法(EO)[13]为代表的诸多启发式算法均已应用于求解多目标电压和无功功率协调优化问题。然而，上述算法由于种群更新和外部存储库维护机制的弱点，收敛速度慢，容易停滞在局部最优，得到的Pareto前沿多样性不佳。 基于知识共享的算法(GSK)是一种新颖的启发式算法，于2020年首次提出[14] ，其灵感来自于在人类的整个生命周期中获取和共享知识的过程[15]。A. W. Mohamed等人在[16]中验证了GSK在用于解决大量基准优化问题时能够呈现出惊人的性能。在这种情况下，为了解决以上问题，GSK采用Pareto支配理论来解决多目标电压和无功功率协调控制优化问题，以改进的IEEE 33节点配电网作为算例系统，并将有功功率损耗、电压偏差和分布式光伏弃光率作为目标函数。结果表明，GWO和EO，GSK得到更为优质的多目标电压和无功功率协调控制优化方案有功功率损耗更小、电压偏差更低、太阳能利用率更高，展现出了更优质的求解性能。
    本文的主要工作可以总结如下：以有功功率损耗、电压偏差和分布式光伏弃光率作为为目标函数，并利用GSK求解改进的IEEE 33节点系统中的多目标电压和无功功率协调控制优化。算例分析的结果证实，与GWO和EO相比，GSK在处理多目标电压和无功功率协调控制优化问题时可以展现更为优秀的优化性能。
    本文的其余部分结构如下：第II节建立了配电网多目标电压和无功功率协调控制优化问题的数学模型；第III部分描述了GSK的基本机制；第IV节描述了利用GSK求解配电网多目标电压和无功功率协调控制优化问题的流程；第V节介绍了仿真环境并展示了算例分析的结果；第VI节对本文的主要工作进行总结。

II. 配电网多目标电压和无功功率协调控制优化的数学模型

A. 控制变量

    对于考虑分布式光伏接入的配电系统，其功率潮流和电压优化控制的方式主要包括调节有载调压变压器抽头、无功补偿装置的补偿容量和分布式光伏设备的有功和无功出力等，因此以各时段有载调压变压器变比$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}$分别为电力系统内的发电机节点总数、有载调压变压器可调抽头总数、并联电容器总组数、分布式电源节点总数。

B. 目标函数

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$之间的电压相位差。

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$处的电压标幺值。

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}$包含的元素个数。

C. 约束条件

1) 等式约束

    配电系统需满足如式(5)所示的配电系统功率潮流平衡约束，该约束是配电网功率潮流和电压优化模型等式约束的主要组成部分
{ P 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 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_{i}(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_{i}(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} ⎩ ⎨ ⎧​Pi​(t)Qi​(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_i(t)$和$Q_i(t)$分别为第$t$个时段内流入节点$i$的总有功功率和总无功功率，计算方法如式(6)所示
{ P i ( t ) = P G i ( t ) − P L i ( t ) Q i ( t ) = Q G i ( t ) − Q L i ( t ) + Q C i ( t ) i ∈ { 1 , 2 , … , N B } \begin{equation} \begin{cases} P_{i}(t) &= P_{Gi}(t)-P_{Li}(t) \\ Q_{i}(t) &= Q_{Gi}(t)-Q_{Li}(t)+Q_{Ci}(t) \end{cases} i∈ \{1,2,…,N_B\} \end{equation} {Pi​(t)Qi​(t)​=PGi​(t)−PLi​(t)=QGi​(t)−QLi​(t)+QCi​(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$处无功补偿设备补偿的无功功率。

2) 不等式约束

    为保障配电系统的安全稳定运行，配电网各节点处的电压、各节点处的总进线功率因数、各支路的电流应满足如式(7)和式(8)所示的不等式约束
{ V i , m i n ≤ V i ( t ) ≤ V i , m a x P i ( t ) P i ( t ) 2 + Q 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_{i}(t)}{\sqrt{P_{i}(t)^2 + Q_{i}(t)^2}} ≥ PF_{min} \end{cases} i∈ \{1,2,…,N_B\} \end{equation} ⎩ ⎨ ⎧​Vi,min​≤Vi​(t)≤Vi,max​Pi​(t)2+Qi​(t)2 ​Pi​(t)​≥PFmin​​i∈{1,2,…,NB​}​​
$$\begin{array}{c}\begin{array}{cc}{I}_{(i,j)}(t)\le I_{(i,j)max}(t)& (i,j)\in L\end{array}\end{array}$$
其中：$I_{(i,j)}(t)$为第$t$个时段内支路$(i,j)$的电流。

III. GSK

A. 基本原理

    GSK是一种2020年提出的受人类行为启发的新型启发式算法，它模拟人类在整个生命周期中获取和分享知识的行为。GSK将人类在整个生命周期中获取和分享知识的过程分为以下两个阶段：
    1) 初级知识获取和分享阶段，即人一生中的早期和中期阶段。在这个阶段，人们更多地通过小型网络（如家庭、邻居、亲戚等）获取知识，而不是通过大型网络（如工作、社交网络、朋友等）。这个阶段的人虽然思想上尚不成熟，但会尝试分享自己的想法。
    2) 高级知识获取与分享阶段，即人一生中的中后期。这个阶段的人通常通过大型网络（如工作、社交网络、朋友等）获取知识。例如，这个阶段的人通常更加追求学习过程中的成效，相信成功者的观点，从而避免失败。

B. 数学描述

    假设第$i$个个体的适应度 x i = [ x i , 1 x i , 2 ⋯ x i , N D ] , i ∈ { 1 , 2 , … , N P } \boldsymbol{x_i} = \left[\begin{matrix} x_{i,1} & x_{i,2} & \cdots & x_{i,N_D}\\\end{matrix}\right],i∈\{1,2,…,N_P\} xi​=[xi,1​​xi,2​​⋯​xi,ND​​​],i∈{1,2,…,NP​}值为$f_i$ ，其中$N_D$是个体${\mathbf{x}}_{\mathbf{i}}$的维数， $N_P$为种群数，任一个体在初级知识获取和分享阶段和高级知识获取与分享阶段中的变化情况如图1所示[14] 。
    根据图1可知，人类个体初级维度和高级维度的数量随时光流逝而变化，GSK则将初级维度和高级维度数量的变化分别量化成式(9)和式(10)所示的形式。
$$\begin{array}{c}{N}_{D,junior}=N_D{\left(1-\frac{k}{K}\right)}^{\alpha }\end{array}$$
$$\begin{array}{c}{N}_{D,senior}=N_D-N_{D,junior}\end{array}$$
其中，$N_{D,junior}$为初级维度的数量， $N_{D,senior}$为高级维度的数量； $k$表示当前迭代次数， $K$表示最大迭代次数； $\alpha (\alpha >0)$为知识学习率。

图1 在初级知识获取和分享阶段和高级知识获取与分享阶段中的任一人类个体 Figure 1: An arbitrary individual in the junior and the senior knowledge gaining and sharing phase

C. 基本流程

    假设使用GSK来处理最小问题，则GSK的基本流程描述如下：
    在初级知识获取和共享阶段，种群中的每个个体更新如下：
    1) 按照适应度值从小到大的顺序对种群中的每个个体进行排序。排序结果为${\mathbf{x}}_{\mathbf{worst}},\dots ,{\mathbf{x}}_{\mathbf{i}-1},{\mathbf{x}}_{\mathbf{i}},{\mathbf{x}}_{\mathbf{i}+1},\dots ,{\mathbf{x}}_{\mathbf{best}}$。
    2) 更新种群中的每个个体。对于个体${\mathbf{x}}_{\mathbf{i}}$，若服从单位均匀分布的随机数$r$小于知识比${\lambda }_r$，则按照式(11)更新种群中的个体
$$\begin{array}{c}{\mathbf{x}}_{\mathbf{i}}^{\mathbf{new}}=\{\begin{array}{l}\begin{array}{cc}{\mathbf{x}}_{\mathbf{i}}+{\lambda }_f\left[{\mathbf{x}}_{\mathbf{i}-1}-{\mathbf{x}}_{\mathbf{i}+1}+{\mathbf{x}}_{\mathbf{r}}-{\mathbf{x}}_{\mathbf{i}}\right]& f({\mathbf{x}}_{\mathbf{i}})>f({\mathbf{x}}_{\mathbf{r}})\\ {\mathbf{x}}_{\mathbf{i}}+{\lambda }_f\left[{\mathbf{x}}_{\mathbf{i}-1}-{\mathbf{x}}_{\mathbf{i}+1}-{\mathbf{x}}_{\mathbf{r}}+{\mathbf{x}}_{\mathbf{i}}\right]& f({\mathbf{x}}_{\mathbf{i}})\le f({\mathbf{x}}_{\mathbf{r}})\end{array}\end{array}\end{array}$$
其中，${\mathbf{x}}_{\mathbf{i}}^{\mathbf{new}}$为新一代个体， ${\mathbf{x}}_{\mathbf{r}}$为从当前代个体中随机选择的个体； ${\lambda }_f$为知识因子；${\lambda }_f>0$且${\lambda }_r\in [0,1]$。否则，即如果$r\ge \lambda r$时，${\mathbf{x}}_{\mathbf{i}}^{\mathbf{new}}={\mathbf{x}}_{\mathbf{i}}$。
    之后，在高级知识获取和共享阶段，种群中的每个个体都按以下方式更新：
    1) 按照适应度值从小到大的顺序对种群中的每个个体进行排序。然后，将排序好的个体分为三类，即最佳个体、一般个体和最差个体。最佳个体所占比例为$p$，最差个体所占比例为$p$，一般个体所占比例为$1-2p$，$p$一般取0.1。
    2) 更新种群中的每个个体。对于个体${\mathbf{x}}_{\mathbf{i}}$，若服从单位均匀分布的随机数$r$小于知识比${\lambda }_r$，则按照式(12)更新种群中的个体
$$\begin{array}{c}{\mathbf{x}}_{\mathbf{i}}^{\mathbf{new}}=\{\begin{array}{l}\begin{array}{cc}{\mathbf{x}}_{\mathbf{i}}+{\lambda }_f\left[{\mathbf{x}}_{\mathbf{rbest}}-{\mathbf{x}}_{\mathbf{rworst}}+{\mathbf{x}}_{\mathbf{rm}}-{\mathbf{x}}_{\mathbf{i}}\right]& f({\mathbf{x}}_{\mathbf{i}})>f({\mathbf{x}}_{\mathbf{r}})\\ {\mathbf{x}}_{\mathbf{i}}+{\lambda }_f\left[{\mathbf{x}}_{\mathbf{rbest}}-{\mathbf{x}}_{\mathbf{rworst}}-{\mathbf{x}}_{\mathbf{rm}}+{\mathbf{x}}_{\mathbf{i}}\right]& f({\mathbf{x}}_{\mathbf{i}})\le f({\mathbf{x}}_{\mathbf{r}})\end{array}\end{array}\end{array}$$
其中，${\mathbf{x}}_{\mathbf{rbest}}$为从最佳个体群中随机选择的个体， ${\mathbf{x}}_{\mathbf{rworst}}$为从最差个体群中随机选择的个体， ${\mathbf{x}}_{\mathbf{rm}}$为从一般个体群中随机选择的个体。否则，即当$r\ge {\lambda }_r$时，${\mathbf{x}}_{\mathbf{i}}^{\mathbf{new}}={\mathbf{x}}_{\mathbf{i}}$。

IV. 利用GSK求解配电网多目标电压和无功功率协调控制优化问题

    利用GSK处理配电网多目标电压和无功功率协调控制优化问题的流程如下：
    步骤1) 检索并处理配电网各监测点的运行数据和各电力设备的铭牌参数。
    步骤2) 根据步骤1)中的运行数据和标称参数，建立配电网配电网多目标电压和无功功率协调控制优化问题的数学模型。
    步骤3) 在GSK中设置种群规模和最大迭代次数等初始参数，初始化种群后计算初始种群中每个个体的目标函数值，并判断是否满足约束条件。基于Pareto支配理论，对所求非支配解集中的每个个体采用[17]所提的快速非支配排序法进行非支配解排序。
    步骤4) 将每个个体分为初级和高级维度，根据式(11)更新初级获得和共享维度，同时利用式(12)更新高级获得和共享维度，再计算新生代中每个个体的目标函数值，判断是否可行。
    步骤5) 判断迭代是否满足终止条件。若不满足终止条件，则转到步骤3)继续迭代；若满足终止条件，停止迭代并得出优化方案。
    据此，利用GSK解决配电网多目标电压和无功功率协调控制优化问题的流程图如图2所示。

图2 基于GSK的配电网多目标电压和无功功率协调控制优化流程图 Figure 2: Flowchart of solving MOVRPOC for the distribution network utilizing GSK

V. 算例分析

A. 仿真环境

    OpenDSS是一种用于配电系统的开源潮流仿真软件，支持各种辐射状、网状配电系统的基波潮流分析、谐波潮流分析、故障分析等相频域分析。OpenDSS也提供了对应MATLAB的接口，这使得MATLAB可以调用OpenDSS计算配电系统的功率潮流，用以解决配电网的OPF问题。
    在本次仿真中，测试系统的CPU为Intel六核i7-9750H 2.60GHz，RAM大小为16GB，开发环境为MATLAB R2021A，OpenDSS版本为9.2.0.1。
    为了验证GSK对Pareto前沿收敛性和多样性的改善效果，本文采用改进的IEEE 33节点配电系统作为算例系统，其拓扑结构如图3所示，系统参数见文献[18]。

图3 改进IEEE 33节点配电系统的拓扑结构 Figure 3: Topology of revised IEEE 33-bus distribution system

B. 优化效果

    本文以$F_{loss}$，$F_{\Delta V}$和$F_{PV}$作为目标函数，将新型的多目标多目标GSK (MOGSK)与现有的多目标GWO (MOGWO)[9]、多目标EO (MOEO)[13]进行比较，验证MOGSK在优化配电网多目标电压和无功功率协调控制问题时能够得到多样性更为优良的Pareto前沿。以1h为时间间隔，以配置参数如表I所示的动态无功功率补偿器(DRPC)作为并联补偿器，进行基于改进IEEE 33节点配电系统的日前功率潮流与电压优化。为公平起见，本文将以上启发式算法的种群大小和最大迭代次数均设置为200。上述算法得到的Pareto前沿如图4所示，MOGSK所得与目标空间中的原点归一化欧氏距离最近的非支配解的电压特性如图5所示。

表I 算例系统的配置参数 TABLE I: Configuration parameters of three selected benchmark networks

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

图4 多目标优化场景下以IEEE 33节点系统为算例系统时的Pareto前沿 Figure 4: Pareto fronts under multi-objective optimization scenarios for revised IEEE 33-bus distribution system

图5 MOGSK所得与目标空间中的原点归一化欧氏距离最近的非支配解的电压特性 Fig. 5: Voltage profile of the solution obtained by MOGSK with the nearest normalized Euclidean distance to the origin point in the objective space

    由于多样性对估计多目标启发式算法性能的重要性，应对MOGWO、MOEO 和MOGSK算法所得的分布性进行定量表征。本文采用中提出的间距(Spacing, SP)指标定量评价所获Pareto前沿的均匀性，$SP$越小则Pareto前沿在目标空间中分布越均匀，最理想的情况即为所得Pareto前沿中的所有非支配解在目标空间中等距均匀分布，此时$SP=0$。

表II 所得Pareto前沿的间距指标 TABLE II: Spacing values of the obtained Pareto fronts

启发式算法	改进IEEE 33节点配电系统
MOGWO	0.0184
MOEO	0.0176
MOGSK	0.0140

    结合图4 、表II和图5所示的Pareto前沿可以看出：
    1) MOGWO得到的Pareto前沿分布较为均匀，但在目标空间的整体分布范围不够广泛。
    2) MOEO在目标空间中得到的Pareto前沿分布范围太窄，过于集中在实际的中间解内，空间分布不均匀。
    3) 改进后的MOGSK得到的Pareto前沿可以广泛均匀地分布在目标空间中。
    4) 各节点电压满足配电网多目标电压无功协调控制优化问题中的节点电压约束。
    总之，新型启发式算法GSK所得配电网多目标电压无功协调控制优化方案多样性优良，能够满足多目标优化场景的要求。

VI. 结论

    本文将改进IEEE 33节点配电系统作为算例系统，以有功功率损耗、电压偏差和分布式光伏弃光率作为目标函数，采用新型启发式算法GSK求解多目标电压和无功功率协调控制优化问题，并以现有的GWO和EO作为为对照，验证了GSK在避免局部最优和Pareto前沿分布多样性方面的改进效果。
    仿真结果表明，GSK得到的Pareto前沿可以更广泛、更均匀地分布在目标空间中。可以看出，改进后的GSK得到的优化结果在收敛性和多样性方面表现出优异的性能，可以为决策者提供更好的多目标电压和无功功率协调控制优化方案，从而更好地满足实际配电网的多样化优化需求。

参考文献

[1] Anco S. Blazev, “Energy Generation and the Environment in the 21st Century,” in Solar Technologies for the 21st Century. New York, NY, USA: River Publishers, 2013, pp. 635-686.
[2] R. A. Walling, R. Saint, R. C. Dugan, J. Burke and L. A. Kojovic, “Summary of Distributed Resources Impact on Power Delivery Systems,” IEEE Transactions on Power Delivery, vol. 23, no. 3, pp. 1636-1644, Jul. 2008, doi: 10.1109/TPWRD.2007.909115.
[3] J. D. Watson et al., “Impact of solar photovoltaics on the low-voltage distribution network in New Zealand,” IET Generation, Transmission & Distribution, vol. 10, no. 1, pp. 1-9, Jan. 2016, doi: 10.1049/IET-GTD.2014.1076.
[4] D. Pappalardo, V. Calderaro and V. Galdi, “Microgrids-Based Approach for Voltage Control in Distribution Systems by an Efficient Sensitivity Analysis Method,” IEEE Systems Journal, vol. 16, no. 3, pp. 4210-4221, Sept. 2022, doi: 10.1109/JSYST.2022.3156798.
[5] S. Mirjalili, S. M. Mirjalili, and A. Lewis, “Grey wolf optimizer,” Advances in engineering software, vol. 69, pp. 46-61, Mar. 2014, doi: 10.1016/J.ADVENGSOFT.2013.12.007.
[6] S. S. Reddy, P. R. Bijwe and A. R. Abhyankar, “Optimal posturing in day-ahead market clearing for uncertainties considering anticipated active-time adjustment costs,” IEEE Systems Journal, vol. 9, no. 1, pp. 177-190, Mar. 2015, doi: 10.1109/JSYST.2013.2265664.
[7] R. H. A. Zubo, G. Mokryani, R. Abd-Alhameed and T. Butler, “A Stochastic Market-based Approach for Operation of Active Distribution Networks”, 2019 IEEE Power & Energy Society General Meeting (PESGM), 2019, pp. 1-5, doi: 10.1109/PESGM40551.2019.8973573.
[8] C. Cai, J. Chen, M. Xi, Y. Tao and Z. Deng, “Multi-Objective Planning of Distributed Photovoltaic Power Generation Based on Multi-Attribute Decision Making Theory,” IEEE Access, vol. 8, pp. 223021-223029, 2020, doi: 10.1109/ACCESS.2020.3042010.
[9] M. Amroune, “Wind integrated optimal power flow considering power losses, voltage deviation, and emission using equilibrium optimization algorithm”, Energy, Ecology and Environment, vol. 7, pp. 369-392, May 2022, doi: 10.1007/S40974-022-00249-2.
[10] A. Q. Abdulrasool and L. T. Al-Bahrani, “Multi-objective Constrained Optimal Power Flow Based on Enhanced Ant Colony System Algorithm”, 2021 12th International Symposium on Advanced Topics in Electrical Engineering (ATEE), 2021, pp. 1-5, doi: 10.1109/ATEE52255.2021.9425082.
[11] S. H. Hemati, G. H. Sheisi, and 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, vol. 45, no. 3, pp. 921-944, Apr. 2021, doi: 10.1007/S40998-021-00419-8.
[12] S. A. El-Sattar, S. Kamel, M. Ebeed, and F. Jurado, “An improved version of salp swarm algorithm for solving optimal power flow problem,” Soft Computing, vol. 25, pp. 4027-4052, Jan. 2021, doi: 10.1007/s00500-020-05431-4.
[13] A. Mostafa, M. Ebeed, S. Kamel and M. A. Abdel-Moamen, “Optimal Power Flow Solution Using Levy Spiral Flight Equilibrium Optimizer With Incorporating CUPFC”, IEEE Access, vol. 9, pp. 69985-69998, 2021, doi: 10.1109/ACCESS.2021.3078115.
[14] A. W. Mohamed, A. A. Hadi and A. K. Mohamed, “Gaining-sharing knowledge based algorithm for solving optimization problems: a novel nature-inspired algorithm”, International Journal of Machine Learning and Cybernetics, vol. 11, pp. 1501–1529, Jul. 2020, doi: 10.1007/S13042-019-01053-X.
[15] H. A. A. Nomer, A. W. Mohamed and A. H. Yousef, “GSK-RL: Adaptive Gaining-sharing Knowledge algorithm using Reinforcement Learning”, 2021 3rd Novel Intelligent and Leading Emerging Sciences Conference (NILES), 2021, pp. 169-174, doi: 10.1109/NILES53778.2021.9600551.
[16] A. W. Mohamed, H. F. Abutarboush, A. A. Hadi and A. K. Mohamed, “Gaining-Sharing Knowledge Based Algorithm With Adaptive Parameters for Engineering Optimization”, IEEE Access, vol. 9, pp. 65934-65946, 2021, doi: 10.1109/ACCESS.2021.3076091.
[17] K. Deb, A. Pratap, S. Agarwal and T. Meyarivan, “A fast and elitist multi-objective genetic algorithm: NSGA-II,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182-197, Apr. 2002, doi: 10.1109/4235.996017.
[18] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” IEEE Transactions on Power Delivery, vol. 4, no. 2, pp. 1401-1407, Apr. 1989, doi: 10.1109/61.25627.

附录

附录1 GSK程序代码

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Gaining-Sharing Knowledge Based Algorithm for Solving Optimization
%%Problems: A Novel Nature-Inspired Algorithm
%% Authors: Ali Wagdy Mohamed, Anas A. Hadi , Ali Khater Mohamed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Best_score, Best_pos, Convergence_curve]=GSK(pop_size,Max_iter,lb,ub,dim,fobj,fopt)
    max_nfes = 10000 * dim;
    %Initialize the positions of search agents
    popold=initialization(pop_size,dim,ub,lb);

    bsf_error_val=[];
    run_funcvals = [];
    G_Max=fix(max_nfes/pop_size);

    %% Initialize the main population
    pop = popold; % the old population becomes the current population

    nfes = 0;
    Best_score = inf;
    Convergence_curve = inf;

    for i = 1 : pop_size
        fitness(i, :) = fobj(pop(i,:),fopt);
        nfes = nfes + 1;
        %%      if nfes > max_nfes; exit(1); end
        if nfes > max_nfes; break; end
        if fitness(i, :) < Best_score
            Best_score = fitness(i, :);
            Best_pos = pop(i, :);            
        end
        run_funcvals = [run_funcvals;Best_score];
    end
%     Convergence_curve(1) = Best_score;

    %%%%%%%%%%%%%%%%%%%%%%%% Parameter settings%%%%%%%%%%
    KF=0.5;% Knowledge Factor
    KR=0.9;%Knowledge Ratio
    K=10*ones(pop_size,1);%Knowledge Rate

    iter=0;
    lu = [lb; ub];
    %% main loop
    while nfes < max_nfes && iter < Max_iter
        iter=iter+1;
        D_Gained_Shared_Junior=ceil((dim)*(1-iter/G_Max).^K);
        D_Gained_Shared_Senior=dim-D_Gained_Shared_Junior;
        pop = popold; % the old population becomes the current population

        [valBest, indBest] = sort(fitness, 'ascend');
        [Rg1, Rg2, Rg3] = Gained_Shared_Junior_R1R2R3(indBest);

        [R1, R2, R3] = Gained_Shared_Senior_R1R2R3(indBest);
        R01=1:pop_size;
        Gained_Shared_Junior=zeros(pop_size, dim);
        ind1=fitness(R01)>fitness(Rg3);

        if(sum(ind1)>0)
            Gained_Shared_Junior (ind1,:)= pop(ind1,:) + KF*ones(sum(ind1), dim) .* (pop(Rg1(ind1),:) - pop(Rg2(ind1),:)+pop(Rg3(ind1), :)-pop(ind1,:)) ;
        end
        ind1=~ind1;
        if(sum(ind1)>0)
            Gained_Shared_Junior(ind1,:) = pop(ind1,:) + KF*ones(sum(ind1), dim) .* (pop(Rg1(ind1),:) - pop(Rg2(ind1),:)+pop(ind1,:)-pop(Rg3(ind1), :)) ;
        end
        R0=1:pop_size;
        Gained_Shared_Senior=zeros(pop_size, dim);
        ind=fitness(R0)>fitness(R2);
        if(sum(ind)>0)
            Gained_Shared_Senior(ind,:) = pop(ind,:) + KF*ones(sum(ind), dim) .* (pop(R1(ind),:) - pop(ind,:) + pop(R2(ind),:) - pop(R3(ind), :)) ;
        end
        ind=~ind;
        if(sum(ind)>0)
            Gained_Shared_Senior(ind,:) = pop(ind,:) + KF*ones(sum(ind), dim) .* (pop(R1(ind),:) - pop(R2(ind),:) + pop(ind,:) - pop(R3(ind), :)) ;
        end
        Gained_Shared_Junior = boundConstraint(Gained_Shared_Junior, pop, lu);
        Gained_Shared_Senior = boundConstraint(Gained_Shared_Senior, pop, lu);

        D_Gained_Shared_Junior_mask=rand(pop_size, dim)<=(D_Gained_Shared_Junior(:, ones(1, dim))./dim);
        D_Gained_Shared_Senior_mask=~D_Gained_Shared_Junior_mask;

        D_Gained_Shared_Junior_rand_mask=rand(pop_size, dim)<=KR*ones(pop_size, dim);
        D_Gained_Shared_Junior_mask=and(D_Gained_Shared_Junior_mask,D_Gained_Shared_Junior_rand_mask);

        D_Gained_Shared_Senior_rand_mask=rand(pop_size, dim)<=KR*ones(pop_size, dim);
        D_Gained_Shared_Senior_mask=and(D_Gained_Shared_Senior_mask,D_Gained_Shared_Senior_rand_mask);
        ui=pop;

        ui(D_Gained_Shared_Junior_mask) = Gained_Shared_Junior(D_Gained_Shared_Junior_mask);
        ui(D_Gained_Shared_Senior_mask) = Gained_Shared_Senior(D_Gained_Shared_Senior_mask);

    %     children_fitness = feval(fhd, ui');
    %     children_fitness = children_fitness';

        for i = 1 : pop_size
            children_fitness(i, :) = fobj(ui(i,:),fopt);
            nfes = nfes + 1;
            if nfes > max_nfes; break; end
            if children_fitness(i) < Best_score
                Best_score = children_fitness(i);
                Best_pos = ui(i, :);
            end
        end

        [fitness, Child_is_better_index] = min([fitness, children_fitness], [], 2);

        popold = pop;
        popold(Child_is_better_index == 2, :) = ui(Child_is_better_index == 2, :);

        % Display Results
        str1 = strcat('Best estimates: Best_pos=');
        disp(str1);
        str2 = strcat(num2str(Best_pos));
        str2 = strcat(str2,'  iteration='); str2=strcat(str2,num2str(iter));
        disp(str2);
%         if rem(iter, 10) == 1 && rem(iter, 100) ~= 11
%             fprintf('After the %d st iteration, best-so-far value = %f\n', iter, Best_score);
%         elseif rem(iter, 10) == 2 && rem(iter, 100) ~= 12
%             fprintf('After the %d nd iteration, best-so-far value = %f\n', iter, Best_score);
%         elseif rem(iter, 10) == 3 && rem(iter, 100) ~= 13
%             fprintf('After the %d rd iteration, best-so-far value = %f\n', iter, Best_score);
%         else
%             fprintf('After the %d th iteration, best-so-far value = %f\n', iter, Best_score);
%         end
%         Convergence_curve(iter+1) = Best_score;
        Convergence_curve(iter) = Best_score;
    end % end while loop


%     fprintf('\nAfter %d iterations, best-so-far value = %f\n', iter, Best_score);
end

附录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;

附录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