考虑风电谐波影响与出力随机性的电力系统无功优化方法

Reactive Power Optimization Considering the Harmonic Effects and Output Uncertainty of Wind Generation

摘要

电容器组等无功补偿装置的优化调整不仅可以减小电力系统网损，还会对电力系统的谐波潮流与谐波网损产生影响.风电不仅存在出力不确定性，而且会产生谐波污染，影响谐波潮流与谐波网损.然而，传统无功优化方法没有同时考虑风电的谐波特性与出力不确定性对谐波潮流和谐波网损的影响，可能导致谐波超标问题，不利于电力系统降低网损.因此，本文首先构建计及风电谐波影响的电力系统无功随机优化模型.该模型通过场景法对风电出力随机性进行建模，目标函数同时考虑基波和谐波网损，约束条件包括基波和谐波潮流约束、电压谐波总畸变率约束等.然后，针对该谐波约束的无功随机优化模型，本文提出了融合驱动力可调粒子群算法和全连接型深度神经网络的高效求解方法 .最后，通过三个修改后的IEEE测试系统对所提模型与方法的有效性进行了验证.

Abstract

The optimal adjustment of reactive power compensation device such as capacitors can not only reduce the network losses， but also present influences on harmonic power flow and harmonic power losses. Wind generators can produce harmonic pollution and present impacts on harmonic power flow and harmonic power losses. However， the effects of the harmonic characteristic and output uncertainty of wind power on harmonic power flow and harmonic power losses are not considered simultaneously in the traditional reactive power optimization methods， which may result in the violation of harmonic standard and is adverse to network losses reduction. In this regard， this paper proposes the reactive power stochastic optimization model for power systems considering the impact of the wind power harmonics. In this model， the uncertainty of wind power is modeled by the scenario method， and the base-frequency network losses and the harmonic losses are considered in the objective function. Also， the constraints such as base-frequency power flow equations， the harmonic power flow equations and the total harmonic voltage distortion constraint are incorporated into the proposed model. After that， focusing on the proposed reactive power stochastic optimization model， the highly efficient method combining the adjustable driving force-based particle swarm optimization and a fully-connected deep neural network is proposed in this paper. Finally， the effectiveness of the proposed model and method is validated by three modified IEEE test systems.

Graphical abstract

关键词

无功功率 / 电力系统潮流 / 谐波分析 / 风电

无功潮流分布对电力系统经济运行具有重要的影响.通过电容器组等无功补偿装置的优化调整，可以改善无功电压分布、降低电力系统网损.由于电容器组等无功补偿装置影响电力系统的阻抗特性，而谐波注入电流对电力系统的影响大小又与电力系统阻抗特性密切相关，因此，电容器组等无功补偿装置的优化调整还会对电力系统的谐波潮流分布与谐波网损产生重要影响^［1］.在新能源电力系统中，诸如风力发电等新能源发电不仅存在出力不确定性，而且会对电力系统产生谐波污染，影响谐波潮流分布和谐波网损^［2］，进而影响电网的电能质量.因此，有必要在电力系统无功优化中同时考虑风力发电的谐波影响与出力不确定性的问题.

迄今为止，学术界已在电力系统无功优化领域开展了许多有益的研究工作.文献［3］提出了一种将无功调度模型转化为凸优化模型的新方法，并进一步针对该凸优化模型提出了基于半正定规划理论的快速求解算法.文献［4］针对含整数变量的无功优化调度问题，提出了一种融合紧圆锥松弛方法和取整方法的高效求解方法，该方法具有近似效果好、最优间隙小等优点.文献［5］构建了一种有载调压变压器新模型及其相应的无功优化半正定优化模型，并提出了一种基于弦分解方法的改进型分支定界算法用于求解前述半正定优化模型.文献［6］将无功优化调度问题松弛为一个混合整数二阶锥问题进行求解，它通过增加凸包络和切平面的方法改善松弛效果.此外，自启发式算法、机器学习等方法也在电力系统无功优化领域得到广泛应用.例如，文献［7］提出了基于控制参数分布范围自适应调整的改进型微分进化算法，以此来快速求解电力系统无功优化调度问题.文献［8］基于Q学习算法提出了一种面向电力系统无功优化调度的全分布式多智能体强化学习方法，该方法无需精确的系统模型且计算效率高.文献［9］考虑多区域电力系统无功优化的特点，提出了一种可克服维数限制、实现分布式协同无功优化的分层相关均衡Q学习算法.文献［10］充分考虑电力系统运行特点，提出了一种融合迁移强化学习和生物启发式优化算法的无功优化方法.文献［11］采用场景法考虑风电的不确定性，提出了基于自启发式算法的电力系统多目标随机无功优化方法.文献［12］针对交直流互联的电力系统，提出了基于动态进 化算法与鲁棒优化算法的多目标无功优化方法.文献［13］提出了基于i-NSGA-Ⅱ算法和风电预测值的电力系统多目标日前无功优化算法，以及考虑预测误差的小时级无功优化方案修正方法.

上述文献虽然在电力系统无功优化领域展开了深入研究，然而它们所提出的方法并没有同时考虑风电的谐波影响与出力不确定性，因而可能导致电力系统出现谐波超标问题，也不利于降低网损.因此，本文同时考虑风电的谐波特性与出力不确定性对谐波潮流与谐波网损的影响，构建计及谐波约束与谐波网损的电力系统无功随机优化模型.该模型采用场景法对风电出力随机性进行建模，目标函数涵盖了基波网损和谐波网损，约束条件包括基波潮流约束、谐波潮流约束、电压谐波总畸变率约束、单次谐波电压含有率约束等.针对所提谐波约束的无功随机优化模型，本文进一步提出了融合驱动力可调粒子群优化算法和全连接型深度神经网络的高效求解方法，该方法具有良好的求解效率.

1 元件谐波模型

1.1 负荷谐波模型

在电力系统中，不同类型的负荷具有不同的谐波特性.根据负荷的谐波特性，负荷可以分为线性负荷和非线性负荷，其中非线性负荷向电网注入谐波电流，产生谐波污染.本文通过阻抗模型对线性负荷进行建模，通过谐波电流源模型对非线性负荷建模，具体如下^［14］：

（1） 线性负荷

y_{L,h}=P_L-\mathrm{j}\frac{\mathrm{1}}{h}{Q}_L（1）

式中：P_L和Q_L分别表示基波频率下负荷的有功功率和无功功率，y_{L,h}表示负荷在次数为h的谐波频率下的等值导纳.

（2） 非线性负荷

{\tilde{I}}_L^{(\mathrm{1})}=\frac{{\left(-P_L-\mathrm{j}{Q}_L\right)}^{*}}{{\left({\tilde{V}}_L^{(\mathrm{1})}\right)}^{*}}（2）

{\tilde{I}}_L^h={\tilde{I}}_{\mathrm{L}}^{(\mathrm{1})}{\tilde{K}}^h={\tilde{I}}_L^{(\mathrm{1})}\left|K^h\right|{\mathrm{e}}^{\mathrm{j}{\phi }_L^h}（3）

式中：{\tilde{I}}_L^{(\mathrm{1})}表示基波频率下负荷的注入电流相量，{\tilde{V}}_L^{(\mathrm{1})}表示基波频率下负荷节点的电压相量，{\tilde{I}}_L^h表示次数为h的谐波频率下负荷的注入电流相量，{\tilde{K}}^h表示{\tilde{I}}_L^h与{\tilde{I}}_L^{(\mathrm{1})}之间的复倍数，\left|K^h\right|和{\phi }_L^h分别表示{\tilde{K}}^h的幅值和相角.

1.2 风电谐波模型

现有研究表明，风力发电可向电网注入谐波电流，产生谐波污染^［2］.对于谐波频率下的风电模型，本文采用谐波电流源对风电进行建模^［15］：

{\tilde{I}}_{\mathrm{W}}^{(\mathrm{1})}=\frac{{\left(P_W+\mathrm{j}{Q}_W\right)}^{*}}{{\left({\tilde{V}}_{\mathrm{W}}^{(\mathrm{1})}\right)}^{*}}（4）

{\tilde{I}}_{\mathrm{W}}^h={\tilde{I}}_{\mathrm{W}}^{(\mathrm{1})}{C}^h（5）

式中：{\tilde{I}}_{\mathrm{W}}^{(\mathrm{1})}表示基波频率下风电的注入电流相量，{\tilde{V}}_{\mathrm{W}}^{(\mathrm{1})}表示基波频率下风电节点的电压相量，{\tilde{I}}_{\mathrm{W}}^h表示次数为h的谐波频率下风电的注入电流相量，谐波电流系数C^h表示{\tilde{I}}_{\mathrm{W}}^h与{\tilde{I}}_{\mathrm{W}}^{(\mathrm{1})}之间的倍数.

1.3 串联和并联元件谐波模型

假定基波频率下，串联元件（譬如线路等值电路中的串联支路）的阻抗为z_{\mathrm{b}}=r_{\mathrm{b}}+\mathrm{j}{x}_{\mathrm{b}}，并联元件（譬如线路等值电路中的对地支路）的导纳为y_{\mathrm{s}\mathrm{h}}=g_{\mathrm{s}\mathrm{h}}+\mathrm{j}{b}_{\mathrm{s}\mathrm{h}}.在次数为h的谐波频率下，串联元件和并联元件的等值阻抗z_{\mathrm{b}}^h和等值导纳y_{\mathrm{s}\mathrm{h}}^h分别表示为

z_{\mathrm{b}}^h=r_{\mathrm{b}}+\mathrm{j}h{x}_{\mathrm{b}}（6）

y_{\mathrm{s}\mathrm{h}}^h=g_{\mathrm{s}\mathrm{h}}+\mathrm{j}h{b}_{\mathrm{s}\mathrm{h}}（7）

2 谐波约束的无功优化模型

2.1 目标函数

\mathrm{m}\mathrm{i}\mathrm{n}\sum _{ij\in {\Omega }_L}\left(∆P_{ij,\mathrm{1},s}^{loss}+\sum _{h\in {\Omega }_H}∆P_{ij,h,s}^{loss}\right)（8）

式中：∆P_{ij,\mathrm{1},s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}表示场景s中支路ij的基波有功损耗，∆P_{ij,h,s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}表示场景s中支路ij在次数为h的谐波频率下的有功损耗；{\Omega }_H表示所考虑的谐波次数集合；{\Omega }_L表示电力系统的支路集合.∆P_{ij,h,s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}可以按如下式子计算：

∆P_{ij,h,s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}={\left|{\tilde{V}}_{i,h,s}-{\tilde{V}}_{j,h,s}\right|}^{\mathrm{2}}\left|G_{ij}^h\right|（9）

式中：{\tilde{V}}_{i,h,s}和{\tilde{V}}_{j,h,s}分别表示谐波次数为h时场景s中节点i和节点j的电压相量，G_{ij}^h表示谐波次数为h的节点导纳矩阵中第i行第j列的元素的实部.令

{\tilde{V}}_{i,h,s}=V_{i,h,s,\mathrm{r}\mathrm{e}}+\mathrm{j}{V}_{i,h,s,\mathrm{i}\mathrm{m}}；

{\tilde{V}}_{j,h,s}=V_{j,h,s,\mathrm{r}\mathrm{e}}+\mathrm{j}{V}_{j,h,s,\mathrm{i}\mathrm{m}}（10）

因此进一步有

∆P_{ij,h,s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}={\left|\left(V_{i,h,s,\mathrm{r}\mathrm{e}}-V_{j,h,s,\mathrm{r}\mathrm{e}}\right)+j\left(V_{i,h,s,\mathrm{i}\mathrm{m}}-V_{j,h,s,\mathrm{i}\mathrm{m}}\right)\right|}^{\mathrm{2}}

\left|G_{ij}^h\right|=\left[{\left(V_{i,h,s,\mathrm{r}\mathrm{e}}-V_{j,h,s,\mathrm{r}\mathrm{e}}\right)}^{\mathrm{2}}+{\left(V_{i,h,s,\mathrm{i}\mathrm{m}}-V_{j,h,s,\mathrm{i}\mathrm{m}}\right)}^{\mathrm{2}}\right]\left|G_{ij}^h\right|（11）

对于基波网损{\sum }_{ij\in {\Omega }_L}∆P_{ij,\mathrm{1},s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}，它可以按照如下式子计算：

\sum _{ij\in {\Omega }_L}∆P_{ij,\mathrm{1},s}^{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}=\sum _{i\in {\Omega }_B}\left(P_{\mathrm{G},i,s}+P_{\mathrm{W},i,s}-P_{\mathrm{D},i,s}\right)（12）

式中：对于母线i而言，P_{\mathrm{G},i,s}表示基波频率下火电在场景s中的有功出力，P_{\mathrm{W},i,s}表示基波频率下风电在场景s中的有功出力，P_{\mathrm{D},i,s}表示基波频率下负荷在场景s中的有功功率，{\Omega }_{\mathrm{B}}表示母线集合.需要注意的是，对于任意非平衡节点的火电，它在各场景中的有功功率为相同的给定数值.

2.2 约束条件

2.2.1 火电约束

Q_{\mathrm{G},i,\mathrm{m}\mathrm{i}\mathrm{n}}\le Q_{\mathrm{G},i,s}\le Q_{\mathrm{G},i,\mathrm{m}\mathrm{a}\mathrm{x}}（13）

V_{\mathrm{G},i,\mathrm{m}\mathrm{i}\mathrm{n}}\le V_{\mathrm{G},i}\le V_{\mathrm{G},i,\mathrm{m}\mathrm{a}\mathrm{x}}（14）

式中：Q_{\mathrm{G},i,s}为基波频率下场景s中母线i处火电的无功功率，Q_{\mathrm{G},i,\mathrm{m}\mathrm{i}\mathrm{n}}和Q_{\mathrm{G},i,\mathrm{m}\mathrm{a}\mathrm{x}}分别为Q_{\mathrm{G},i,s}的下限和上限，V_{\mathrm{G},i}表示与火电相连的母线i的基波电压幅值，V_{\mathrm{G},i,\mathrm{m}\mathrm{i}\mathrm{n}}和V_{\mathrm{G},i,\mathrm{m}\mathrm{a}\mathrm{x}}分别为V_{\mathrm{G},i}的下限和上限.

2.2.2 并联电容约束

B_{\mathrm{C},i,\mathrm{m}\mathrm{i}\mathrm{n}}\le B_{\mathrm{C},i}\le B_{\mathrm{C},i,\mathrm{m}\mathrm{a}\mathrm{x}}（15）

式中：B_{\mathrm{C},i}表示投入节点i运行的并联电容器组的电纳，它是本文所提无功随机优化模型的离散型决策变量，B_{\mathrm{C},i,\mathrm{m}\mathrm{i}\mathrm{n}}和B_{\mathrm{C},i,\mathrm{m}\mathrm{a}\mathrm{x}}分别表示B_{\mathrm{C},i}的下限和上限.由于B_{\mathrm{C},i}为离散型决策变量，因此进一步令∆B_{\mathrm{C},i}表示投入节点i运行的并联电容器组电纳的调节步长.

2.2.3 风电约束

P_{\mathrm{W},i,s}\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{W},i}=Q_{\mathrm{W},i,s}（16）

{\beta }_{\mathrm{W},i,\mathrm{m}\mathrm{i}\mathrm{n}}\le {\beta }_{\mathrm{W},i}\le {\beta }_{\mathrm{W},i,\mathrm{m}\mathrm{a}\mathrm{x}}（17）

式中：对于母线i而言，Q_{\mathrm{W},i,s}表示场景s中风电的无功功率，{\beta }_{\mathrm{W},i}表示风电的功率因数角，{\beta }_{\mathrm{W},i,\mathrm{m}\mathrm{i}\mathrm{n}}和{\beta }_{\mathrm{W},i,\mathrm{m}\mathrm{a}\mathrm{x}}分别表示风电功率因数角的下限和上限.值得注意的是，在本文所提无功随机优化模型中，对于同一个风电功率因数角，它在任意场景中都是相同的数值.

2.2.4 基波潮流约束

P_{\mathrm{G},i,s}+P_{\mathrm{W},i,s}-P_{\mathrm{D},i,s}=V_{i,s}\sum _{j\in i}{V}_{j,s}\left(G_{ij}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{ij,s}+B_{ij}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ij,s}\right)（18）

Q_{\mathrm{G},i,s}+Q_{\mathrm{W},i,s}-Q_{\mathrm{D},i,s}=V_{i,s}\left[-V_{i,s}{B}_{\mathrm{C},i}+\sum _{j\in i}{V}_{j,s}\left(G_{ij}\mathrm{s}\mathrm{i}\mathrm{n}{\theta }_{ij,s}-B_{ij}\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{ij,s}\right)\right]（19）

式中：对于母线i而言，Q_{\mathrm{D},i,s}表示基波频率下负荷在场景s中的无功功率，V_{i,s}和V_{j,s}分别表示节点i和节点j在场景s中的基波电压幅值，{\theta }_{ij,s}表示场景s中节点i与节点j之间的基波电压相角差，G_{ij}与B_{ij}分别表示基波频率下不考虑并联电容器时节点导纳矩阵元素ij的实部和虚部.

2.2.5 谐波潮流约束

{\tilde{I}}_{i,h,s}=\sum _j{Y}_{ij}^h{\tilde{V}}_{j,h,s}（20）

式中：对于次数为h的谐波频率而言，{\tilde{I}}_{i,h,s}表示场景s中节点i的谐波注入电流相量，Y_{ij}^h表示节点导纳矩阵下标为ij的元素.由于式（20）为复数形式，因此将其转化为两个等效的实数方程：

I_{i,h,s,\mathrm{r}\mathrm{e}}=\sum _j{G}_{ij}^h{V}_{j,h,s,\mathrm{r}\mathrm{e}}-B_{ij}^h{V}_{j,h,s,\mathrm{i}\mathrm{m}}

I_{i,h,s,\mathrm{i}\mathrm{m}}=\sum _j{G}_{ij}^h{V}_{j,h,s,\mathrm{i}\mathrm{m}}+B_{ij}^h{V}_{j,h,s,\mathrm{r}\mathrm{e}}（21）

式中：B_{ij}^h表示Y_{ij}^h的虚部.

2.2.6 节点电压约束

V_{i,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{i}\mathrm{n}}\le V_{i,\mathrm{r}\mathrm{m}\mathrm{s},s}\le V_{i,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}（22）

V_{i,\mathrm{r}\mathrm{m}\mathrm{s},s}=\sqrt[]{{\left(V_{i,s}\right)}^{\mathrm{2}}+\sum _{h\in {\Omega }_H}{\left(V_{i,h,s}\right)}^{\mathrm{2}}}（23）

式中：对于次数为h的谐波频率而言，V_{i,h,s}表示场景s中母线i的谐波电压幅值，V_{i,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{i}\mathrm{n}}和V_{i,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}分别表示V_{i,\mathrm{r}\mathrm{m}\mathrm{s},s}的下限和上限.式（23）中的{\left(V_{i,h,s}\right)}^{\mathrm{2}}按如下式子计算：

{\left(V_{i,h,s}\right)}^{\mathrm{2}}={\left(V_{i,h,s,\mathrm{r}\mathrm{e}}\right)}^{\mathrm{2}}+{\left(V_{i,h,s,\mathrm{i}\mathrm{m}}\right)}^{\mathrm{2}}

2.2.7 支路电流约束

I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{i}\mathrm{n}}\le I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},s}\le I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}（24）

I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},s}=\sqrt[]{{\left(I_{ij,s}\right)}^{\mathrm{2}}+\sum _{h\in {\Omega }_H}{\left(I_{ij,h,s}\right)}^{\mathrm{2}}}（25）

式中：I_{ij,s}表示场景s中基波频率下支路ij的传输电流.对于次数为h的谐波频率而言，I_{ij,h,s}表示场景s中支路ij的传输电流.I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{i}\mathrm{n}}和I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}分别表示I_{ij,\mathrm{r}\mathrm{m}\mathrm{s},s}的下限和上限.

2.2.8 谐波约束

\tilde{I}\frac{V_{i,h,s}}{V_{i,s}}\le V_{h,\mathrm{L}\mathrm{I}\mathrm{M}}（26）

\sum _{h\in {\Omega }_H}{\left(V_{i,h,s}\right)}^{\mathrm{2}}\le {\left(V_{i,s}{V}_{T,\mathrm{L}\mathrm{I}\mathrm{M}}\right)}^{\mathrm{2}}（27）

式中：对于母线i而言，式（26）是关于单次谐波电压含有率的约束条件，式（27）是关于电压谐波总畸变率的约束条件，V_{h,\mathrm{L}\mathrm{I}\mathrm{M}}表示次数为h的谐波电压占有率的上限，V_{T,\mathrm{L}\mathrm{I}\mathrm{M}}表示电压谐波总畸变率的上限.

2.3 场景建模

令P_{\mathrm{W},i}={\xi }_{\mathrm{w},i}{P}_{\mathrm{W},i}^r，其中P_{\mathrm{W},i}^r表示风电的额定有功功率，{\xi }_{w,i}表示风电实际出力系数.受风速不确定性的影响，风电实际出力系数{\xi }_{\mathrm{w},i}的数值是不确定的.在文献［16］中，Beta分布被用于拟合风电出力预测误差，因此本文将Beta分布用于风电实际出力系数{\xi }_{\mathrm{w},i}建模，即认为{\xi }_{\mathrm{w},i}服从Beta分布，它的概率密度函数如下：

f({\xi }_{\mathrm{w},\mathrm{i}})=\frac{\mathrm{\Gamma }\left(A+B\right)}{\mathrm{\Gamma }\left(A\right)\mathrm{\Gamma }\left(B\right)}{\left({\xi }_{\mathrm{w},i}\right)}^{A-\mathrm{1}}{\left({\xi }_{\mathrm{w},i}\right)}^{B-\mathrm{1}}（28）

式中：\mathrm{\Gamma }\left(A\right)表示以A为自变量的\mathrm{\Gamma }函数，\mathrm{\Gamma }\left(A+B\right)表示以\left(A+B\right)为自变量的\mathrm{\Gamma }函数，\mathrm{\Gamma }\left(B\right)表示以B为自变量的\mathrm{\Gamma }函数.

给定风电实际出力系数的概率分布后，本文进一步采用文献［17］所提的场景生成方法构建风电实际出力系数集合，然后根据P_{\mathrm{W},i}={\xi }_{\mathrm{w},i}{P}_{\mathrm{W},i}^r得到风电有功出力场景.由于风电有功出力不确定性远大于负荷有功功率的不确定性^［18］，因此，为简化分析，本文仅仅考虑风电出力的不确定性.故在本文中，负荷在基波频率下的有功功率的数值不因场景变化而变化，负荷在基波频率下的无功功率的数值也不因场景变化而变化.

3 求解方法

本文所提谐波约束的无功随机优化模型考虑了离散控制的电容器组，是一个典型的非线性混合整数优化问题，求解难度高.基于此，本文针对所建立的谐波约束无功优化模型，提出融合驱动力可调粒子群优化算法（Adjustable driving force based particle swarm optimization，ADFPSO）^［19］与全连接型深度神经网络^［20］的两阶段求解方法.在该方法中，第一阶段采用ADFPSO算法来获取所提无功优化模型的整数型变量解.为加快ADFPSO算法的求解速度，本文将构建一个全连接型深度神经网络Net1来计算适应度.第二阶段根据前述整数型变量的解将本文所提模型转化为连续型优化问题M_cq，并通过联合使用深度神经网络与优化软件对优化问题M_cq 进行求解，得到所提无功优化模型的连续型决策变量数值.

3.1 融合ADFPSO算法与全连接型深度神经网络的第一阶段求解

3.1.1 ADFPSO算法

ADFPSO算法是一种融合粒子群算法与多驱动力思想的自启发式算法，它在复杂优化问题上表现出了良好的优化能力^［19］.ADFPSO算法融合了两种驱动力机制：高适应度和高新颖度，即通过适应度和新颖度来驱动粒子寻优搜索.

在ADFPSO算法中，粒子速度与位置的更新采用如下策略：

v_{ij}^{m+\mathrm{1}}=w{v}_{ij}^m+c_{\mathrm{1}}{r}_{\mathrm{1},j}\left(P_{i_{r\mathrm{1}}j}-x_{ij}^m\right)+c_{\mathrm{2}}{r}_{\mathrm{2},j}\left({\mathrm{\mathscr{H}}}_{i_{r\mathrm{2}}j}-x_{ij}^m\right)（29）

x_{\mathit{i},\mathit{m}+\mathrm{1}}=x_{\mathit{i},\mathit{m}}+v_{\mathit{i},\mathit{m}+\mathrm{1}}（30）

式中：x_{\mathit{i},\mathit{m}}表示第m次迭代时粒子i的位置；v_{\mathit{i},\mathit{m}+\mathrm{1}}表示第（m+1）次迭代时粒子i的速度；w为惯性权重系数，根据文献［19］，它是可以随着迭代次数而变化的.i_{r\mathrm{1}}和i_{r\mathrm{2}}都是分布于［1， M］的整数；P_{i_{r\mathrm{1}},j}与{\mathrm{\mathscr{H}}}_{i_{r\mathrm{2}}j}分别表示从高适应度个体极值库{\Omega }_F、高新颖度个体极值库{\Omega }_N随机选取的个体极值，M=⌈p\times N_p⌉，N_p表示种群规模，p的计算公式如下：

p=\mathrm{0.7}\frac{\mathrm{1}}{\mathrm{1}+{\mathrm{e}}^{\mathrm{0.001}\times \frac{F_{\mathrm{E}}-k_{\mathrm{m}\mathrm{a}\mathrm{x}}/\mathrm{2}}{D}}}+\mathrm{0.1}（31）

式中：F_{\mathrm{E}}表示适应度函数的计算次数，k_max表示适应度函数计算次数的上限，D表示粒子位置的维数.

对于粒子的新颖度，它的定义如下：

{\mathrm{𝒩}}_{i,m}=\frac{\mathrm{1}}{K}{\displaystyle \sum _{j=\mathrm{1}}^K}dist(x_{i,m}-{\mu }_{i,j,m})（32）

式中：{\mathrm{𝒩}}_{i,m}表示第m次迭代时粒子i的新颖度，K表示粒子i的邻居个数，{\mu }_{i,j,m}表示第m次迭代时粒子i的第j个邻居的位置，\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(x_{i,m}-{\mu }_{i,j,m})表示第m次迭代时粒子i与它的第j个邻居之间的欧式距离.

具体到本文所提的无功优化问题，对于ADFPSO算法的粒子i在第m次迭代时的位置x_{i,m}，它由本文所提模型的离散型决策变量组成.适应度函数与给定离散型决策变量数值下的目标函数最优值有关.在ADFPSO算法中，适应度函数的计算时间是影响计算速度的重要因素.为了加快适应度函数的计算速度，本文将采用全连接型深度神经网络计算个体所在位置的适应度值.

3.1.2 基于全连接型深度神经网络的适应度快速计算方法

如前所述，在本文所提方法中，ADFPSO算法的粒子适应度与给定离散型决策变量数值下的目标函数最优值有关.对于给定粒子，由于离散型决策变量的数值已知，此时本文所提无功优化模型变成连续型非线性优化模型M_cq .与原来的混合整数优化问题相比，连续型非线性优化模型M_cq 的求解难度低.但ADFPSO算法需要进行大量的适应度函数计算，因而对单次计算适应度函数的速度要求很高.如果直接采用诸如GAMS等优化软件求解非线性优化模型M_cq，那么ADFPSO算法的计算效率将难以得到保证.

全连接型深度神经网络在电力系统、通信等领域展现出了良好的性能^［21-23］.其中，全连接型深度神经网络已被用于提升电力系统最优潮流问题的求解速度^［22-23］.因此本文将构建全连接型深度神经网络Net1用于适应度值的近似计算，减少ADFPSO算法中适应度函数的单次计算时间，进而减少ADFPSO算法的总计算时间.本文将通过优化软件构建全连接型深度神经网络Net1的训练样本.需要注意的是，对于本文所提模型而言，给定离散型决策变量（电容器组电纳）的数值，非线性优化问题M_cq 可能出现无解的情况.由于在无解情况下无法通过目标函数的最优值给出样本标签，因此本文采用如下形式的样本标签（即ADFPSO算法的适应度函数）：

f\left(x\right)=\mathrm{O}\mathrm{b}\mathrm{j}+\mathrm{I}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{s}（33）

式中：\mathrm{I}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{s}表示模型不可行惩罚值，如果所提模型采用优化软件求解得到了不可行解，那么令\mathrm{I}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{s}等于给定的正数，否则\mathrm{I}\mathrm{n}\mathrm{f}\mathrm{e}\mathrm{a}\mathrm{s}=0；\mathrm{O}\mathrm{b}\mathrm{j}为优化软件返回的目标函数值.

对于全连接型深度神经网络Net1，采用如下向量作为它的输入：

{\mathit{X}}_{\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{1}}=\left[{\mathit{P}}_{\mathrm{L}},{\mathit{Q}}_{\mathrm{L}},{\mathit{A}}_{\mathrm{W}},{\mathit{B}}_{\mathrm{W}},{\mathit{B}}_{\mathrm{C}}\right]（34）

式中：{\mathit{P}}_{\mathrm{L}}和{\mathit{Q}}_{\mathrm{L}}分别表示有功负荷的行向量与无功负荷的行向量.{\mathit{A}}_{\mathrm{W}}\mathrm{和}{\mathit{B}}_{\mathrm{W}}分别表示由风电实际出力系数Beta分布的A参数与B参数构成的行向量；{\mathit{B}}_{\mathrm{C}}表示由电容器组电纳组成的行向量.需要注意的是，本文通过归一化的样本数据构建全连接型深度神经网络Net1，并将它直接用于第一阶段求解.

3.1.3 整数型决策变量求解流程

本文所提基于ADFPSO算法与深度神经网络的整数型决策变量求解流程如下：

步骤1：设置ADFPSO算法的参数初值、种群大小N_p 、适应度函数计算次数上限k_max.令高适应度个体极值库{\Omega }_{\mathrm{F}}=\mathrm{\varnothing }，高新颖度个体极值库{\Omega }_{\mathrm{N}}=\mathrm{\varnothing }.

步骤2：令适应度函数计算次数k=0，迭代次数m=0.

步骤3：随机初始化每个粒子的位置和速度.

步骤4：通过全连接型深度神经网络Net1计算每一个粒子的适应度值.

步骤5：确定个体极值：若m=1，则将每个粒子当前的位置设置为个体极值；若m>1，则对于每个粒子而言，根据它的当前适应度值、当前个体极值及适应度值，对它的个体极值进行更新.

步骤6：对所有粒子的当前个体极值根据适应度进行排序，然后选择排名前M个的个体极值进入高适应度个体极值库{\Omega }_{\mathrm{F}}.

步骤7：对所有粒子的当前个体极值根据新颖度进行排序，然后选择排名前M个的个体极值进入高新颖度个体极值库{\mathrm{\Omega }}_{\mathrm{N}}.

步骤8：根据式（29）更新粒子的速度，然后根据式（30）更新粒子的位置.

步骤9：判断粒子的位置是否超出允许范围.如果粒子位置的某个分量超出其允许范围，则将相应的粒子位置分量强制设置为其限值.例如，假设仅有一个节点含有电容器组，如果电容器组的电纳B_{C,i,m+\mathrm{1}}不在其限值范围内，则进行如下操作：

（35）

由于电容器组的电纳为离散型决策变量，因此进一步做如下处理：

B_{\mathrm{C},i,m+\mathrm{1}}=∆B_{\mathrm{C},i}\times \mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\left(B_{\mathrm{C},i,m+\mathrm{1}}/∆B_{\mathrm{C},i}\right)（36）

式中：函数\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\left(B_{\mathrm{C},i,m+\mathrm{1}}/∆B_{\mathrm{C},i}\right)表示对\left(B_{\mathrm{C},i,m+\mathrm{1}}/∆B_{\mathrm{C},i}\right)的数值进行取整.

步骤10：令k=k+N_p，m=m+1，并更新ADFPSO算法的其他参数.

步骤11：判断适应度函数的计算次数是否大于预先设定的最大值k_max.如果k>k_max，则结束迭代求解过程；否则，进入步骤4.

3.2 基于全连接型深度神经网络与优化软件的第二阶段求解

通过融合ADFPSO算法与全连接型深度神经网络，本文所提方法第一阶段可以获得离散型决策变量的解.但是，前述方法没有得到连续型决策变量的解.因此，本文所提方法的第二阶段将联合使用全连接型深度神经网络与优化软件，求解所提无功优化模型中连续型决策变量的数值.具体而言，首先根据第一阶段求解方法得到的离散型决策变量的数值，将所提无功优化模型转化为非线性连续型数学优化问题M_cq .然后，针对该连续型数学优化问题，通过全连接型深度神经网络与优化软件进行联合求解.

3.2.1 基于全连接型深度神经网络的连续型决策变量初值求解

本文将构建一个全连接型深度神经网络Net2，它的输出向量由所提无功随机优化模型中的连续型决策变量构成，它的输入向量{\mathit{X}}_{\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{2}}与全连接型深度神经网络Net1的输入向量{\mathit{X}}_{\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{1}}一致，即

{\mathit{X}}_{\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{2}}={\mathit{X}}_{\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{1}}（37）

在本文中，全连接型深度神经网络Net2的输出量将作为优化软件求解连续型数学优化问题M_cq 时决策变量的初值.

3.2.2 基于优化软件的连续型决策变量数值求解

针对连续型数学优化问题M_cq，本文利用优化软件GAMS中的非线性求解器对其进行求解，获取所提无功随机优化问题的连续型决策变量的数值.采用GAMS软件的非线性求解器进行求解时，连续型数学优化问题M_cq 决策变量的初值来自深度神经网络Net2的输出.

4 算例分析

为了验证所提谐波约束的无功随机优化模型及其求解方法的有效性，本文将采用修改后的IEEE 14节点系统、IEEE 30节点系统、IEEE 57节点系统进行仿真分析.上述三个IEEE系统的原始数据来自软件MATPOWER^［24］.本文将原始数据中PV节点的发电机全部考虑为火电.表1和表2分别给出了仿真分析采用的风电和电容器组参数.在本文所有仿真中，单次谐波电压含有率和电压谐波总畸变率的上限分别设为3%和5%.V_{i,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{a}\mathrm{x}}和V_{i,\mathrm{r}\mathrm{m}\mathrm{s},\mathrm{m}\mathrm{i}\mathrm{n}}的标幺值分别设为1.05 pu和0.95 pu，以下分析中电压均采用标幺制.通过借鉴文献［2］关于风电谐波影响的分析、文献［15］关于风电谐波次数的设置，本文假设风电谐波输出电流的次数为17和19.为了重点研究风电谐波电流的影响，本文在仿真分析中做了简化处理：假设所有的负荷为线性负荷.

4.1 风电谐波电流的影响

本小节将通过分析不同的风电谐波电流水平，研究所提无功优化模型与方法的有效性.本小节采用修改后的IEEE 14节点系统进行仿真，所采用的风电容量为40 MW（本文中的风电容量指风电的额定有功功率）.

为了分析考虑谐波约束的必要性，本小节将风电谐波电流系数C^h假设为0.02，并采用谐波电流系数C^h={\mathrm{10}}^{-\mathrm{9}}时的无功优化方案作为对比对象.当C^h={\mathrm{10}}^{-\mathrm{9}}时，风电注入电网的谐波电流已经非常小，因此可以近似认为：此种情况下通过本文方法得到的无功优化方案没有考虑风电的谐波电流.本文将此情况下得到的无功优化方案简称为无谐波方案.本节仿真分析中风电有功出力的概率参数为：A= 1.5，B= 1.5.上述参数下的仿真结果如图1~图4所示.

从图1可以看出，当C^h=\mathrm{0.02}时，若将无谐波方案作为无功优化方案，则节点电压幅值的最大值为1.051 0 pu，很明显这超出了它的上限（1.05 pu）.从图2可以看出，此时所有节点中17次谐波电压含有率的最大值为3.22%，这超出了单次谐波含有率的上限（3%）.这表明无谐波方案是不可行的.

如图3和图4所示，与上面现象形成鲜明对比的是，本文方法所得无功优化方案（即采用C^h=\mathrm{0.02}的条件得到的无功优化方案）可以保证节点电压幅值不超过1.05 pu，单次谐波电压含有率不超过3%，即能够保证节点电压和单次谐波电压含有率满足要求.因此，本文所提方法得到的无功优化方案是可行的.

表3给出了不同方案的决策变量数值，其中{\beta }_W表示风电的功率因数角，B_{{\mathrm{C}}_i}（i=1，2，3）表示第i个电容器组的电纳值，V_gi （i=1，2，3，4）表示第i个火电的节点电压幅值.从表中可以看出，考虑谐波方案和无谐波方案存在明显的不同.譬如，在考虑谐波方案中，B_{{\mathrm{C}}_{\mathrm{1}}}=0.06 pu，B_{{\mathrm{C}}_{\mathrm{3}}}=0.02 pu.而在无谐波方案中，B_{{\mathrm{C}}_{\mathrm{1}}}=0.1 pu，B_{{\mathrm{C}}_{\mathrm{3}}}=0.08 pu.因此，考虑谐波方案和无谐波方案的B_{{\mathrm{C}}_{\mathrm{1}}}和B_{{\mathrm{C}}_{\mathrm{3}}}数值明显不同，两个方案存在明显的差异性.这说明了本文所提方法考虑风电谐波电流的必要性.

4.2 不同风电容量下的综合网损分析

本小节通过修改后的IEEE 14节点系统分析风电容量对所提无功优化方法的影响.在本小节中，风电容量P_i^r表示为P_i^r=k_i{P}_{i,\mathrm{0}}^r，其中P_{i,\mathrm{0}}^r表示4.1小节所采用的风电容量，k_i表示风电容量系数.

本小节将风电容量系数k_i从0.85变化到1.05，变化间隔为0.05，所得到的仿真结果如图5~图6所示.在图5~图6中，负荷水平1中的节点负荷与4.1小节一致，负荷水平2中的节点负荷为4.1小节相应数值的95%.从图5-6的仿真结果可以看出：

1） 如果负荷水平、风电容量、参数A不变，那么在修改后的IEEE 14节点系统中，综合网损（谐波网损与基波网损之和）随着Beta分布参数B的增加而增加，即综合网损与Beta分布参数B之间呈现正相关关系.

2） 在Beta分布参数的不同组合下（参数组合1：A=1.5，B=4.5；参数组合2：A=1.5，B=3），如果修改后的IEEE 14节点系统负荷水平保持不变，那么它的综合网损随着风电容量系数的增加而不断减小.

4.3 计算性能分析

本小节将采用修改后的IEEE 14节点系统、IEEE 30节点系统、IEEE 57节点系统，分析所提融合ADFPSO算法与全连接型深度神经网络的模型求解方法的有效性.其中，本小节将修改后的IEEE 14节点系统与IEEE 57节点系统的负荷分别设为原始数据的105%和95%.在本小节中，首先采用的对比对象为直接通过GAMS软件求解得到的.仿真结果如图7~图12所示（图中标题在仿真系统的名称前省略了“修改后的”）.

从图7~图12可以看出：本文所提求解方法的计算速度显著优于通过GAMS软件直接求解这种方法.例如，在修改后的IEEE 14节点系统中，本文所提方法的计算时间减少率的最小值为18.68%；在修改后的IEEE 30节点系统中，本文所提方法的计算时间减少率的最小值为89.53%；在修改后的IEEE 57节点系统中，本文所提方法的计算时间减少率的最小值为94.76%.另外，从目标函数上看，无论是修改后的IEEE 14节点系统还是修改后的IEEE 30节点系统，本文所提方法的目标函数值增加率都小于0.5%.而对于修改后的IEEE 57节点系统而言，本文所提方法的目标函数值增加率的最大值为2.97%.因此，在上述修改后的IEEE系统中，本文所提方法的计算时间减少率的数值远大于目标函数值增加率的数值，即本文所提方法可以通过牺牲很小的最优性，取得模型求解速度的显著提升.

为了进一步分析本文方法的有效性，本文以修改后的IEEE 30节点系统为例，通过将所提方法中的ADFPSO算法替换为正余弦优化算法（Sine Cosine Algorithm，SCA）^［25］与经典PSO算法，分别得到作为对比的算法{\mathrm{𝒜}}_{\mathrm{1}}和算法{\mathrm{𝒜}}_{\mathrm{2}}.本小节采用两种概率参数组合（参数组合3：A=3，B=3；参数组合4：A=5，B=3）进行仿真分析，本文方法、算法{\mathrm{𝒜}}_{\mathrm{1}}和算法{\mathrm{𝒜}}_{\mathrm{2}}第一阶段求解的结果如图13~图14所示（图中结果为通过归一化的数据进行仿真得到的结果）.为了便于展示，图13~图14中ADFPSO、SCA、PSO分别指代本文所提方法、算法{\mathrm{𝒜}}_{\mathrm{1}}和算法{\mathrm{𝒜}}_{\mathrm{2}}.从这两张图可以看出，在风电出力的两种概率参数下，本文所采用的ADFPSO算法能够收敛到更好的适应度函数值，这说明本文采用的ADFPSO算法优于SCA算法和经典PSO算法.

5 结 论

本文计及风电谐波特性与出力不确定性的影响，建立了考虑谐波约束与谐波网损的电力系统无功随机优化模型，提出了融合数据驱动方法和ADFPSO算法的高效求解方法.本文通过修改后的IEEE系统对所提模型与方法的有效性进行了分析.仿真结果表明：

1）本文所提方法得到的无功优化方案可以避免节点电压谐波超标问题以及其他安全约束越限问题.

2）本文所提融合数据驱动方法和ADFPSO算法的模型求解方法可通过牺牲很小的最优性，带来模型求解效率的显著提升.

参考文献

[1] EAJAL A A， EL-HAWARY M E ．Optimal capacitor placement and sizing in unbalanced distribution systems with harmonics consideration using particle swarm optimization［J］．IEEE Transactions on Power Delivery，2010，25（3）：1734-1741．

[2] ZHAO Y Q， MILANOVIĆ J V ．Equivalent modelling of wind farms for probabilistic harmonic propagation studies［J］．IEEE Transactions on Power Delivery，2022，37（1）：603-611．

[3] DAVOODI E， BABAEI E， MOHAMMADI-IVATLOO B，et al ．A novel fast semidefinite programming-based approach for optimal reactive power dispatch［J］．IEEE Transactions on Industrial Informatics，2020，16（1）：288-298．

[4] BINGANE C， ANJOS M F， LE DIGABEL S ．Tight-and-cheap conic relaxation for the optimal reactive power dispatch problem［J］．IEEE Transactions on Power Systems，2019，34（6）：4684-4693．

[5] SANTIAGO G CONSTANTE F， LÓPEZ J C， RIDER M J ．Optimal reactive power dispatch with discrete controllers using a branch-and-bound algorithm：a semidefinite relaxation approach［J］．IEEE Transactions on Power Systems，2021，36（5）：4539-4550．

[6] KAYACıK S E， KOCUK B ．An MISOCP-based solution approach to the reactive optimal power flow problem［J］．IEEE Transactions on Power Systems，2021，36（1）：529-532．

[7] NIU M， XU N Z， DONG H N，et al ．Adaptive range composite differential evolution for fast optimal reactive power dispatch［J］．IEEE Access，2021，9：20117-20126．

[8] XU Y L， ZHANG W， LIU W X，et al ．Multiagent-based reinforcement learning for optimal reactive power dispatch［J］．IEEE Transactions on Systems，Man，and Cybernetics，Part C （Applications and Reviews），2012，42（6）：1742-1751．

[9] TAN M， HAN C J， ZHANG X S，et al ．Hierarchically correlated equilibrium Q-learning for multi-area decentralized collaborative reactive power optimization［J］．CSEE Journal of Power and Energy Systems，2016，2（3）：65-72．

[10] ZHANG X S， YU T， YANG B，et al ．Accelerating bio-inspired optimizer with transfer reinforcement learning for reactive power optimization［J］．Knowledge-Based Systems，2017，116：26-38．

[11] HASSAN M H， KAMEL S， EL-DABAH M A，et al ．Optimal reactive power dispatch with time-varying demand and renewable energy uncertainty using Rao-3 algorithm［J］．IEEE Access，2021，9：23264-23283．

[12] HUI Q， TENG Y， ZUO H，et al ．Reactive power multi-objective optimization for multi-terminal AC/DC interconnected power systems under wind power fluctuation［J］．CSEE Journal of Power and Energy Systems，2020，6（3）：630-637．

[13] LIU Y C， ĆETENOVIĆ D， LI H Y，et al ．An optimized multi-objective reactive power dispatch strategy based on improved genetic algorithm for wind power integrated systems［J］．International Journal of Electrical Power & Energy Systems，2022，136：107764．

[14] YOUSSEF K H ．Power quality constrained optimal management of unbalanced smart microgrids during scheduled multiple transitions between grid-connected and islanded modes［J］．IEEE Transactions on Smart Grid，2017，8（1）：457-464．

[15] PANDI V R， ZEINELDIN H H， XIAO W D ．Determining optimal location and size of distributed generation resources considering harmonic and protection coordination limits［J］．IEEE Transactions on Power Systems，2013，28（2）：1245-1254．

[16] LUO Y H， NIE Q B， YANG D S，et al ．Robust optimal operation of active distribution network based on minimum confidence interval of distributed energy beta distribution［J］．Journal of Modern Power Systems and Clean Energy，2021，9（2）：423-430．

[17] NIKNAM T， ZARE M， AGHAEI J ．Scenario-based multiobjective volt/var control in distribution networks including renewable energy sources［J］．IEEE Transactions on Power Delivery，2012，27（4）：2004-2019．

[18] HETZER J， YU D C， BHATTARAI K ．An economic dispatch model incorporating wind power［J］．IEEE Transactions on Energy Conversion，2008，23（2）：603-611．

[19] YU F， TONG L， XIA X W ．Adjustable driving force based particle swarm optimization algorithm［J］．Information Sciences：An International Journal，2022，609（C）：60-78．

[20] SZE V， CHEN Y H， YANG T J，et al ．Efficient processing of deep neural networks：a tutorial and survey［J］．Proceedings of the IEEE，2017，105（12）：2295-2329．

[21] DOBRE E I， LAMPE L ．Data-driven soft demapping for residual impairments channels［J］．IEEE Communications Letters，2022，26（11）：2611-2615．

[22] SINGH M K， KEKATOS V， GIANNAKIS G B ．Learning to solve the AC-OPF using sensitivity-informed deep neural networks［J］．IEEE Transactions on Power Systems，2022，37（4）：2833-2846．

[23] LOTFI A， PIRNIA M ．Constraint-guided deep neural network for solving optimal power flow［J］．Electric Power Systems Research，2022，211：108353．

[24] ZIMMERMAN R D， MURILLO-SÁNCHEZ C E， THOMAS R J ．MATPOWER：steady-state operations，planning，and analysis tools for power systems research and education［J］．IEEE Transactions on Power Systems，2011，26（1）：12-19．

[25] MIRJALILI S ．SCA：a sine cosine algorithm for solving optimization problems［J］．Knowledge-Based Systems，2016，96：120-133．

基金资助

国家自然科学基金资助项目（U22B200134）