基于改进蜉蝣算法的一种新型无功优化补偿方法及其应用

彭泽森，舒 恺，高飞翎，余萃卓

（1.福州大学电气工程与自动化学院，福建福州 350108；2.宁波市电力设计院有限公司，浙江宁波 315000）

0 引言

电力系统的无功优化是一个在电网运行中时刻存在的问题，如何将系统的无功功率分布调节到最优、系统运行效益达到最高，是无功优化的主要任务。无功优化属使电网的电压水平合格并使整个系统的运行经济效益对于混合整数规划带约束问题，既包含连续变量也包含离散变量，求解时需要满足系统的潮流约束要求以及各变量的操作和安全运行要求，这就导致无功优化变得十分复杂。为了解决无功优化算法问题，陆续有学者提出了线性规划法（如灵敏度分析法[1]、内点法[2]等）、非线性规划法（简化梯度法[3]、牛顿法[4]等）、动态规划法等应用到无功优化领域，也取得的较好的效果。

随着人工智能的发展，群体智能搜索算法发展迅速。群体智能搜索算法通过模拟自然界各类动物的社会行为[5]，经过多次的迭代，从而找到问题的最优解。群体搜索算法容易陷入局部最优，通过对基础算法进行改良，能够使算法性能变优，使算法能够更好地应用于各类实际问题中。无功优化领域自20世纪以来，引进了各种启发式智能寻优方法。如文献[6]通过算例将遗传算法（Genetic Algorithm，GA）引入到无功优化领域，验证了人工智能类算法相较于传统方法求解无功优化问题的优势，同时指出GA 求解速度较慢的缺点。为了改进基础智能算法的性能使智能算法更适用于无功优化问题，文献[7]提出了在萤火虫算法（Firefly Algorithm，FA）中引入精英保留策略和新的排序机制；而文献[8]则综合粒子群算法（Particle Swarm Optimization，PSO）的随机搜索性和禁忌搜索算法（Tabu Search Algorithm，TSA）的短期记忆性提出了禁忌搜索-粒子群算法，两种思路均对基础算法的性能起到了一定的提升效果。但文献[6-8]中所提及的方法，都是将自然界中生物的日常行动以及生物繁衍两类生物现象分离开或者是将两种智能算法简单结合，算法思路与自然规律仍存在一定差距，从而影响了算法性能。

蜉蝣算法（Mayfly Algorithm，MA）是一种新型的群体智能搜索算法[9]，在2020 年由Konstantinos Zervoudakis 提出。MA 是基于自然界中蜉蝣的社会行为，特别是交配行为，而提出的一种群体“启发式”寻优算法。MA 比其他算法更加符自然规律，具有寻优速度快、搜索能力强等优点，是粒子群算法、遗传算法的结合成果。这种算法提出以后，其在各学科应用中倍受关注。

本文主要研究MA 在电力系统无功优化中的应用。首先针对MA 的算法流程进行了阐述，将其与电力系统无功优化模型相结合，对蜉蝣寻优算法做了相应的优化和改进，提出一种新型智能搜索无功优化补偿方法，并利用标准案例进行了有效验证。

1 无功优化模型

电力系统的无功优化是一个混合整数规划带约束问题，控制变量有发电机机端电压、支路变压器分接头位置、静止电容器投入组数等，状态变量有节点电压大小、相位等[10-11]。以下将从目标函数、变量约束条件两个方面描述其数学模型。

1.1 目标函数

对电力系统无功优化的目标函数进行选取时，应当结合系统本身以及系统的技术指标和经济效益[12]。常以有功功率损耗、发电成本、无功补偿设备动作次数等作为目标，寻求最优调节策略。降低有功网损能够减小电网运行成本并提高电力系统的电能质量[13]，因此本文以有功功率损耗最小为优化的目标。假设系统中有L 个节点，N 条支路，则目标函数为：

式中：Vs(s=1,2,…,L)为节点s的电压幅值；Vz(z=1,2,…,L)为节点z的电压幅值；Gsz为节点s和节点z的互导；θsz=θs-θz为节点s和节点z的相角差值。

1.2 控制变量约束条件

无论系统运行在何种状态下，都应满足系统有功、无功潮流约束方程（即无功优化的等式约束条件），而各设备和系统状态量在实际运行时则需要满足操作要求和运行要求（即无功优化的不等式约束条件）[14]。

1.2.1 等式约束条件

尽管系统的拓扑结构复杂多变，在对其进行问题求解时，都应时刻保证整个系统的功率平衡，即满足潮流方程：

式中：PGs为节点s的有功出力；PLs为节点s的有功负荷；QGs为节点s的无功出力；QLs为节点s的无功负荷；QCs为补偿设备在节点s投入的无功大小；Bsz为节点s，z之间的电纳。

1.2.2 不等式约束条件

无功优化策略中输出的控制变量大小必须满足设备的动作要求[15]，同时不能使系统状态量越限。

可直接调节变量的不等式约束：

式中：NG为系统可调节的发电机数；NT为系统中可调变压器的个数；NQ为系统安装了无功补偿设备的节点总和；UGsmin，UGsmax分别为系统发电机端电压的最小值和最大值；Tsmin，Tsmax分别为变压器分抽头位置的最小位置和最大位置；QCsmin，QCsmax分别为无功最小补偿量和最大补偿量。

不可直接调节变量约束条件：

式中：Usmin，Usmax分别为负荷节点电压变动的下限值与上限值；NL为负荷节点数。

2 蜉蝣智能优化算法

2.1 基础蜉蝣算法

MA 是通过观察蜉蝣的社会行为而提出一种新的群体搜索智能算法[16]。蜉蝣的社会行为主要包括飞行、交配以及死亡。MA 将蜉蝣的社会行为抽象出来，利用蜉蝣的飞行位置作为问题的解，通过不断地修改蜉蝣的飞行位置并在每次迭代结束后对个体进行筛选，只留下确定数量的蜉蝣个体，从而找到问题的最优解。为了模拟蜉蝣的交配行为，MA中将个体分为雄性个体与雌性个体。雄性个体的移动行为受到自身历史最优以及全局历史最优影响，而雌性个体的移动取决于其雄性配偶。当雌性个体的适应度差于配偶的适应度时，则雌性向配偶位置移动，反之则自由飞行[17]。同时，雌性个体与匹配的雄性配偶在每次迭代中会产生两个后代，从而增大算法的搜索范围。算法的具体实现方法如下：

假设有一个D维问题，第i个蜉蝣的飞行位置xi=[xi1，xi2，xi3，...，xiD]，其对应的速度为vi=[vi1，vi2，vi3，...，viD]。MA 的具体实现方法如下完成。

2.1.1 雄性个体的移动

雄性蜉蝣喜欢成群聚集，其位置的移动主要依据：自身所到过最优位置以及全局最优位置来进行。将自身适应度与全局最优适应度进行比较，当自身适应度优于全局最优适应度时，则按照自身所到过最优位置、全局最优位置与自身位置的关系进行调整；若自身适应度差于全局最优适应度，则依据婚礼舞蹈系数以及惯性权重进行调整[18]。蜉蝣移动速度的变化可表述如下：

因此，雄性个体的位置更新公式为：

2.1.2 雌性个体的移动

与雄性个体不同，雌性蜉蝣不会成群聚集，飞行位置的移动主要取决于与其配对的雄性蜉蝣。当雌性蜉蝣被雄性蜉蝣吸引时，即适应度差于雄性蜉蝣时，会追随雄性蜉蝣的运动轨迹，否则雌性蜉蝣将在空间中随机飞行。雌性蜉蝣的移动速度的变化可描述如下：

2.1.3 雄雌交配行为

蜉蝣的社会行为中还包括交配行为，其交配过程可简单描述为：先按照适应度分别对雌性蜉蝣和雄性蜉蝣进行排序，相同适应度等级的雌性蜉蝣与雄性蜉蝣匹配并进行交配[19]，繁衍出两个子代，其子代的产生方式如下：

式中：offspring1为生成的第1 个子代；offspring2为生成的第2 个子代；A为随机数，取值范围为（0，1）；xm，ym分别为雄性个体与雌性个体。

该步骤使MA 的搜索范围扩大化，并在子代产生后通过筛选操作加速MA 收敛至最优解处。

同时，为了避免蜉蝣在飞行时采用过大的飞行速度而飞出可行域空间的情况，应对其飞行速度加以限制，即设定一个最大飞行速度Vmax，则飞行速度的更新规则如下：

式中：Vmax=rand·(xmax-xmin)，rand 为（0，1）的随机数。

2.1.4 子代变异行为

与其他智能算法类似，MA 在对问题进行寻优求解时，会出现“早熟”行为，即收敛于局部最优解，从而导致最终结果不是全局最优[20]。因此，本文对子代进行了变异操作，并使婚礼舞蹈系数与随机飞行系数随迭代次数减小，具体过程如下：

为了进行子代的变异操作，通过将高斯变异加入到MA 中，对子代的单个维度进行变异。由于子代变异的数量不能过多，故设变异率mu，种群大小为Pop。以Pop为基准，变异的子代数量为[Pop·mu]。变异行为描述如下：

式中：offspringp为生成的第p个子代；σ为正态分布的标准差；N（0，1）是平均值0，方差为1 的标准正态分布。

为了增强寻优过程的全局搜索能力，使婚礼舞蹈系数d与随机飞行系数fl随着迭代次数的增加而减小，描述如下：

式中：dn，d0，ddamp分别为婚礼舞蹈系数的第n次迭代数值、初始值以及衰减系数；fln，fl0，fldamp分别为随机飞行系数的第n次迭代数值、初始值以及衰减系数。

2.2 蜉蝣算法的改进与优化

通过分析已有文献对MA 的研究[21]，发现基础MA 在对问题进行求解时，收敛速度过慢，且在寻优时容易陷入局部最优解。因此，本文提出将Levy 飞行与随机惯性权重系数引入到基础MA 中，对常规MA 进行优化，形成一种改进的MA，这种新型算法能有效避免寻优时陷入局部最优解，增强其全局搜索能力。

2.2.1 Levy飞行

Levy 飞行属于一种连续时间的随机游动，由较长时间的短步长和较短时间的长步长组成，其大概率在小范围内游走，而小概率在大范围内游走，具有平稳、独立的步长。通过观察发现，自然界中许多生物的觅食轨迹以及飞行轨迹，都符合Levy 飞行的特征[22]。

Levy（λ）表示服从当前迭代次数t 的随机分布，其概率分布式如下：

式中：λ为幂次数；μ为服从正态分布的随机数。

Levy 分布是一个带有重尾的概率分布，利用计算机模拟Levy 飞行时，通常采用Mantegna 提出的模拟Levy 飞行路径的计算式：

式中：s为Levy 飞行的步长；η为服从正态分布的随机数；参数γ为常数，一般取值为1.5。

参数μ、η为服从式（17）所示的正态分布的随机数。

式（17）中，方差δμ，δη定义如下：

式中：β取值范围为[0,2]。

结合Levy 飞行摆脱局部最优解的能力以及MA 本身的特征，将Levy 飞行引入到MA 后雌雄个体的位置更新公式如下：

式中：Levy（λ）为跳跃步长，服从Levy 分布的随机搜索向量；⊗为矢量运算。

雄性蜉蝣喜欢成群聚集，因此当本次迭代适应度优于个体最优时，采用式（19）进行位置更新，增大解的搜索范围。而雌性蜉蝣的位置移动受雄性蜉蝣的影响，故当雌性蜉蝣适应度优于雄性蜉蝣时，采用式（20）进行更新。采取两种不同的位置更新规则对雌雄蜉蝣进行位置更新，能将Levy 飞行更好地应用到MA 中。

2.2.2 随机惯性权重系数

惯性权重系数g的取值大小，决定了迭代步长受历史位置的影响程度。在求解初期，需要增大解的搜索范围；而到了求解后期，则需要对解变化范围进行缩小，从而逐渐收敛至最优解处。为了使算法在不同阶段能够获得合适的惯性权重系数，采用随机惯性权重系数取代算法原来算法固定的惯性权重系数[23]，更新公式如下：

式中：gmin，gmax分别为惯性权重系数的最小值与最大值；rand（）为（0，1）内均匀分布的随机数；randn（）为正态分布的随机数；σ为方差，表示g与其数学期望的偏差程度，有利于算法向期望方向进化。

采用式（21）作为惯性权重系数，有利于增强算法的种群的多样性和搜索能力。在搜索初期，惯性权重系数可能取到较大的值，从而使算法的前期搜索范围增大；到了接近收敛的末期，惯性权重系数可能取到较小的值，使得解逐渐收敛于全局最优。

3 改进MA新型无功优化方法

结合电力系统无功优化问题的数学模型以及蜉蝣算法的实现过程，将蜉蝣的位置对应于无功优化模型中控制变量的大小，包括发电机节点电压UGi、有载调压变压器变比Ti、无功补偿投入组数QCi，可调节变量数则为蜉蝣个体的维数D，第i个蜉蝣Xi的位置表示为Xi=[xi1，xi2…，xiD]=[UG，Ti，QCi][24]。

利用蜉蝣算法进行电力系统无功优化的具体实现流程如下:

步骤1）读取系统原始数据，包括节点信息和支路信息，初始化基本参数，包括种群大小、迭代次数等。

步骤2）依据控制变量取值范围随机生成初始种群，将种群划分为雌性和雄性两类，并计算初始潮流和初始网损，从而获得各种群的适应度值，设置初始速度为0。

步骤3）根据式（7）和式（20）分别更新雌蜉蝣的速度与位置。同时将更新后的位置代入潮流计算程序中获得新适应度值，新适应度值优于个体最优则更新个体最优。

步骤4）根据式（5）和式（19）分别更新雄蜉蝣的速度与位置。通过潮流计算程序获得雄蜉蝣的新适应度，若优于个体最优则更新个体最优，同时与全局最优比较，若优于全局最优则更新全局最优。

步骤5）分别对雌雄个体按照适应度由优到差进行排列，同等级的雌雄个体按照式（9）和式（10）进行交配操作以产生子代蜉蝣。

步骤6）根据式（12）对部分子代蜉蝣进行变异操作。

步骤7）对空间中所有的蜉蝣个体按照适应度值进行排序，筛选出与初始种群数量相同的优秀个体。

步骤8）判断迭代计算退出条件是否达到，未达到则返回步骤3，达到则不再计算并输出最优解。

4 算例验证

为了验证改进的MA（Professional Mayfly Algorithm，ProMA）的可行性，选择标准IEEE30 节点为测试系统，分别采用GA，PSO，TSA，FA，MA 以及本文提出的改进的MA 算法对系统状态进行优化，优化时通过控制发电机端电压、可调变压器变比以及并联静止电容器来达到降低网损的目的。

IEEE30 系统由30 个节点共同组成，拓扑中存在6 个发电机节点（节点1、2、5、8、11、13），其余节点为PV 节点（即负荷节点）、4 条支路包含可调变压器（支路4—12，6—9，6—10，27—28）及存在两个无功补偿节点10，24（通过并联电容器调节），选择节点1 为平衡节点。因此，选定的控制变量为发电机电压（UG1，UG2，UG3，UG4，UG5）、变压器变比（T6-9，T6-10，T4-12，T27-28）、补偿设备容量（Qc10,Qc24）。控制变量上下限及调节档位见表1。其中，系统的基准容量为100 MVA。

表1 控制变量约束条件Table 1 Control variable constraints p.u.

各算法迭代时适应度进化曲线如图1 所示，最终的寻优结果如表2、表3 所示。

表2 IEEE30各算法无功优化结果对比Table 2 Comparison of reactive power optimizations with different algorithms in IEEE 30-node system p.u.

图1 适应度进化曲线对比图Fig.1 Comparison of fitness evolution curves

从图1 可看出，GA，PSO，TSA，FA，MA 以及ProMA 均能达到降低网损的目的，在经过60 次的迭代计算后，各算法均能收敛至各自的最优解处。分析各算法的适应度曲线可发现，ProMA 的收敛速度最快，在迭代20 次后便收敛于最优解处且最终适应度最小，而GA、PSO、TSA、FA、MA 分别在第40次、第50 次、第60 次、第41 次、第30 次才收敛。除ProMA 外，其他几种算法均出现了不同程度的“早熟”现象，即收敛至局部最优，其中FA 最为严重。改进前的MA 在迭代过程中出现了两次“早熟”现象，分别出现在第8 次迭代后以及第16 次迭代后，严重影响收敛速度，最终寻优结果也略差于PSO。在将Levy 飞行以及随机惯性权重系数加入MA 后，ProMA 能跳脱出局部最优解，从而沿着全局最优的方向继续进化。ProMA 引入了新的位置更新策略，增强了寻优过程中种群的多样性以及跳脱局部最优的能力，使MA 算法的全局搜索能力得到提升，算法性能优于GA、PSO、TSA、FA、MA 算法。

分析表2 中6 种算法寻优后的控制变量大小可知，GA，PSO，TSA 和MA 的控制策略中，两台并联电容器组的投入组数总和分别为10 组、9 组、9 组、9 组，而ProMA 投入的电容器组数总和为5 组，明显降低了电容组的投入数量，从而降低了系统无功优化成本。相比于其他算法提供的优化策略经济效益高，充分利用了发电机组来对系统进行优化，在保证系统运行效率最高的前提下，不过多外加补偿设备投入运行，同时也满足了系统运行的经济性。

对比表3 中各算法的总迭代时间可知，ProMA运行时间最短，而MA 次之，说明了本文提出的ProMA 具备优良的求解效率。分析表3 中各算法的降损效益可知，各算法得出的优化策略均可使系统的有功网损下降，其中ProMA 相较于其他算法降损效益最高，证明了ProMA 优良的求解性能并验证了MA 在电力系统无功优化中的可行性。

表3 IEEE 30 节点系统优化后的有功网损Table 3 Active network loss after optimization in IEEE 30-node system

5 结论

本文提出将MA 引入到电力系统无功优化中，利用Levy 飞行以及随机惯性权重系数对基础MA进行优化改进，使MA 更加适应于无功优化问题中。通过对比PSO、MA、ProMA 在IEEE30 节点系统的优化结果，证明了ProMA 在无功优化问题的求解效率高、跳脱局部最优能力强以及优化策略中无功补偿设备投入量少等方面都要具有明显的优势，是一种新型有效的无功优化方法。该方法的提出为以后MA 在电力系统调压降损问题的求解上提供了新思路以及理论基础。