An Improved MOEA/D Algorithm for the Solution of the Multi-Objective Optimal Power Flow Problem

Abstract

The optimal power flow (OPF) is an important tool for the secure and economic operation of the power system. It attracts many researchers to pay close attention. Many algorithms are used to solve the OPF problem. The decomposition-based multi-objective algorithm (MOEA/D) is one of them. However, the effectiveness of the algorithm decreases as the size of the power system increases. Therefore, an improved MOEA/D (IMOEA/D) is proposed in this paper to solve the OPF problem. The main goal of IMOEA/D is to speed up the convergence of the algorithm and increase species diversity. To achieve this goal, three improvement strategies are introduced. Firstly, the competition strategy between the barnacle optimization algorithm and differential evolution algorithm is adopted to overcome the reduced species diversity. Secondly, an adaptive mutation strategy is employed to enhance species diversity at the latter stage of iteration. Finally, the selective candidate with similarity selection is used to balance the exploration and exploitation capabilities of the proposed algorithm. Simulation experiments are performed on IEEE 30-bus and IEEE 57-bus test systems. The obtained results show that the above three measures can effectively improve the diversity of the population, and also demonstrate the competitiveness and effectiveness of the proposed algorithm for the OPF problem.

1. Introduction

The OPF optimizes specified objective functions by adjusting control variables and is subject to the operating constraints [1]. In power systems, the OPF problem is a complex and constrained optimization problem that frequently requires the optimization of many objectives. At the same time, different objective functions in a multi-objective OPF problem often contradict each other, which might make classical optimization approaches, such as traditional interior point methods [2], Newton methods [3], and programming methods [4], difficult to solve. Unlike classical optimization methods, which can only find one optimal solution in one trial, multi-objective evolutionary algorithms (MOEAs) can obtain a set of multiple non-dominated Pareto solutions in one trial and have demonstrated great performances on the multi-objective OPF problem.

Depending on different methods, MOEAs can be classified into Pareto-based approaches and decomposition-based approaches. The NSGA-II [5] and SPEA-II [6] algorithms are classical Pareto-based algorithms. However, the performance of these algorithms is unsatisfactory when applied to the multi-objective OPF problem [7]. Because as the number of objectives increases, the population size of non-dominant individuals also increases. Compared with Pareto-based approaches, decomposition-based approaches provide a new idea for solving the multi-objective OPF problem [8,9,10], in which MOEA/D [11] is the most prominent and has attracted great interest from scholars in recent years.

In MOEA/D, a multi-objective optimization problem is decomposed into a set of scalar optimization sub-problems, which are solved at the same time. The traditional MOEA/D and its variants have been widely used in the multi-objective OPF problem and have been proven to have a good performance. Zhu et al. adopted the MOEA/D algorithm to solve the environmental economic dispatch problem [12]. Zhu et al. also used MOEA/D to optimize the cost and emissions of a wind thermal power generation system and obtained better solution quality and computational efficiency [13]. Li and Fang proposed an improved MOEA/D-M2M framework to solve the environmental/economic dynamic dispatching problem [14]. Biswas et al. combined MOEA/D with the superiority of feasible solutions to deal with constraints in an OPF problem [15]. Zhang et al. proposed a multi-stage dynamic resource allocation strategy MOEA/D to solve the OPF problem [16]. In the above-mentioned literature, the improvement of the update method for the individual is relatively limited. Therefore, it is necessary to carry out deeper research on the individual update approach of MOEA/D to solve the multi-objective OPF problem effectively.

In the standard MOEA/D algorithm, new individuals are generated by the differential evolution algorithm (DE). The way of updating individuals is single, which may lead to slow convergence and insufficient species diversity. In this paper, the population diversity of the MOEA/D algorithm is systematically studied, and three improvements are proposed: Firstly, we use the barnacle optimization algorithm (BMO) to compete with DE to improve population diversity. Secondly, an adaptive mutation strategy is employed instead of a fixed mutation rate. Finally, the use of similarity selection balances the exploration and exploitation capabilities of the algorithm.

The main contributions of this paper are as follows:

(1) The competitive evolution strategy, adaptive mutation strategy, and selective candidate with similarity selection are introduced to improve the performance of MOEA/D.

(2) An improved MOEA/D is proposed to better solve the OPF problem.

(3) IEEE 30 and 57-bus systems are employed to verify the efficiency of the proposed IMOEA/D. The simulation results demonstrate that the convergence and diversity of the IMOEA/D are enhanced when dealing with the OPF problem.

The rest of the paper is organized as follows: In Section 2, we first briefly introduce the mathematical model of the multi-objective OPF problem. Then, the framework of the MOEA/D and related concepts are described in Section 3. In Section 4, we describe the improved strategy and IMOEA/D in detail. In Section 5, we explain the results and data analysis of the proposed IMOEA/D for the IEEE 30-bus system and the IEEE 57-bus system. Section 6 is the conclusion.

2. OPF Mathematical Model

When solving a multi-objective OPF problem, it is necessary not only to deal with multiple conflicting objective functions but also to satisfy the basic constraints. Therefore, the general formulation of the multi-objective optimization problem is described as follows:

$$\begin{array}{c}min\left[f_1\left(x\right)f_2\left(x\right)\cdots f_M\left(x\right)\right]\hfill \\ \mathrm{s}.t.\hfill \\ \left\{\begin{array}{c}{g}_i\left(x\right)⩾0i=\mathrm{1,2},\cdots ,G\\ h_j\left(x\right)=0j=\mathrm{1,2},\cdots ,H\end{array}\right.\hfill \end{array}$$

(1)

where $M$ is the number of objectives, $f_k\left(x\right)$ is the $k$-th objective function, $k$ = 1,2, …$M.$ $g_i\left(x\right)$ is the $i$-th inequality constraint, $h_j\left(x\right)$ is the $j$-th inequality constraint, $G$ and $H$ are the numbers of inequality and equality constraints, respectively. The objectives and constraints of the multi-OPF problem are presented next.

2.1. Objectives

This paper selects four common objective functions, which are fuel cost, emission, real power loss, and voltage deviation [17,18].

Cost: Fuel cost is the most fundamental objective function in the power system. The relationship between cost and power is shown in Equation (2).

$$\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}=\sum _{i=1}^{NG}{a}_i+b_i{P}_{Gi}+c_i{P}_{Gi}^2$$

(2)

where $NG$ is the number of generators, $a_i$, $b_i$, and $c_i$ are the cost coefficients of the $i$-th generator, and $P_{Gi}$ is the active power of the $i$-th generator.

Emission: Harmful gases such as SO_x, NO_x, and CO_x are generated when using fuel to generate electricity. With the improvement of people’s awareness of environmental protection, power plants must reduce the emission of harmful gases while ensuring the quality of power generation. The relationship between emission and real power is given by Equation (3).

$$\mathrm{Emission}=\sum _{i=1}^{NG}\left[\left({\alpha }_i+{\beta }_i{P}_{Gi}+{\gamma }_i{P}_{Gi}^2\right)+{\omega }_i{e}^{\left({\mu }_i{P}_{Gi}\right)}\right]$$

(3)

where ${\alpha }_i$, ${\beta }_i$, ${\gamma }_i$, ${\omega }_i$, and ${\mu }_i$ are the emission coefficients.

Power loss: Power loss is inevitable during power transmission due to the resistance of power transmission lines. The power loss is computed by Equation (4).

$$\mathrm{Loss}=\sum _{q=1}^{nl}{G}_{q\left(ij\right)}\left[V_i^2+V_j^2-2V_i{V}_j\cos \left({\delta }_{ij}\right)\right]$$

(4)

where ${\delta }_{ij}={\delta }_i-{\delta }_j$ is the voltage angle difference between buses $i$ and $j$. $V_i$ and $V_j$ are the voltage magnitude of the $i$-th bus and the $j$-th bus, respectively (either generator or load). $G_{q\left(ij\right)}$ is the transfer conductance of branch $q$ connecting buses $i$ and $j$, and $nl$ is the number of transmission lines.

Voltage deviation: The voltage deviation is an important indicator to measure the safety of a power system. The indicator is defined as the cumulative sum of deviations of voltages at all load buses in the network from the desired value. The aggregate voltage deviation is calculated by Equation (5).

$$\mathrm{V}\mathrm{D}=\left(\sum _{p=1}^{NL}\left|V_{Lp}-1\right|\right)$$

(5)

where $V_{Lp}$ is the voltage magnitude of the $p$-th load bus.

2.2. Constraints

Equality constraints: The equality constraints are the power balance equations, presented in Equations (6) and (7).

$$P_{Gi}-P_{Di}-V_i\sum _{j=1}^{NB}{V}_j\left[G_{ij}\cos \left({\delta }_{ij}\right)+B_{ij}\sin \left({\delta }_{ij}\right)\right]=0\forall i\in NB$$

(6)

$$Q_{Gi}-Q_{Di}-V_i\sum _{j=1}^{NB}{V}_j\left[G_{ij}\sin \left({\delta }_{ij}\right)-B_{ij}\cos \left({\delta }_{ij}\right)\right]=0\forall i\in NB$$

(7)

where $NB$ is the total number of buses, and $P_D$ and $Q_D$ are active and reactive power demands, respectively. $G_{ij}$ is the transfer conductance and $B_{ij}$ is the susceptance between buses $i$ and $j$.

Generator constraints:

$$V_{Gi}^{min}\le V_{Gi}\le V_{Gi}^{max}\forall i\in NG$$

(8)

$$P_{Gi}^{min}\le P_{Gi}\le P_{Gi}^{max}\forall i\in NG$$

(9)

$$T_j^{min}\le T_j\le T_j^{max}\forall j\in NT$$

(10)

where $T_j$ is the tap of the transformer at the $j$-th branch, and $NT$ is the number of transformers.

Shunt compensator constraints:

$$Q_{Ck}^{min}\le Q_{Ck}\le Q_{Ck}^{max}\forall k\in NC$$

(11)

where $Q_{Ck}$ is the shunt compensation at the k-th bus, and $NC$ is the number of shunt VAR compensators.

Security constraints:

$$V_{Lp}^{min}\le V_{Lp}\le V_{Lp}^{max}\forall p\in NL$$

(12)

$$S_{lq}\le S_{lq}^{max}\forall q\in nl$$

(13)

where $S_{lq}$ is the loading of the $q$-th line. $V_{Lp}$ is the voltage magnitude of the $p$-th load bus. $NL$ and $nl$ are the numbers of load buses and transmission lines, respectively.

3. MOEA/D Algorithm and Related Concepts

3.1. Pareto Optimality

The model of a multi-objective problem is shown in Equation (1), where $\Omega$ represents the decision space, and the $R^m$ represents the target space. Let $u,v$∈$R^m$. If $u_i\ge v_i$ for any i and $u_j>v_j$ for at least one index $j$ $(i,j\in \{1,2,...,m\left\}\right),$ then we say that $u$ dominates $v$. If there is no point $F\left(y\right)$ that can dominate the point $F\left(x\right)$ in the decision space, then $x$ is the Pareto optimal, and $F\left(x\right)$ is called the Pareto optimal vector. In other words, the improvement of the Pareto optimal point on a certain objective function will cause the degradation of at least one other objective function. The set of all Pareto optimal solutions is called the Pareto set, and the set of all optimal vectors is called the Pareto front [19].

3.2. Tchebycheff Aggregation Function

In a multi-objective problem, when the number of objectives is more than three, the sorting of the quality of the optimization results will be more complicated. MOEA/D uses the Tchebycheff aggregation function to convert multiple target values into a single target value, making it easier to judge the quality of the results [10].

$$g^{te}\left(x\mid {\lambda }^m,z^{*}\right)=\underset{1\le i\le M}{max}\left\{{\lambda }_i^m\left|\frac{f_i\left(x\right)-z_i^{*}}{z_i^{\mathrm{nad}}-z_i^{*}}\right|\right\}$$

(14)

where ${\lambda }^m$ is the weight vector. $m$ is the number of the objectives, and $m\in \{1,2,...,N\}$. $f_i\left(x\right)$ is the $i$-th objective function value. $z_i^{\mathrm{nad}}$ is the maximum value of the $i$-th objective function., $z_i^{*}$ is the minimum value of the $i$-th objective function.

3.3. Introduction of the MOEA/D Algorithm

As mentioned above, MOEA/D uses the decomposition method to divide the overall multi-objective problem into $N$ individuals. Each individual is optimized simultaneously in the iterative process, and the aggregation function is used to normalize them to obtain a single objective value to reflect the quality of the optimization results. The MOEA/D algorithm process is as follows [10]:

Step1: Initialization

(1) Initialize the population variable {$x^1,x^2,\dots {,x}^N$} and calculate the target value { $F\left(x^1\right),F\left(x^2\right),\dots ,F\left(x^N\right)$ }, where $N$ is the population size, $F\left(x^i\right)=\left\{f_1\left(x^i\right),f_2\left(x^i\right),\dots ,f_M\left(x^i\right)\right\}(i\in \{1,2,\dots ,N\left\}\right)$, $M$ is the number of objectives.

(2) Initialize the weight vectors{${\lambda }^1{,\lambda }^2,\dots ,{\lambda }^N$} and determine the neighborhood of each individual.

(3) Let $z^i$ be the optimal value of $F\left(x^i\right)$ in the current generation, and the ideal point is $z^{*}\in \left\{z^1,z^2,\dots ,z^N\right\}.$

Step2: Update

For each individual in the population, use the following three steps to update:

(1) Generation of new individuals: Two solutions are randomly selected from the neighborhood of the individual $x$, and the DE algorithm is used to perform crossover and mutation operations on the individual to generate a new solution $y$. Calculate the target value $F\left(y\right)$ and aggregation function value $g\left(y\right)$.

(2) Update individual: If $g\left(y\right)>g\left(x\right)$, replace $x$ with $y$ and update the neighborhood. Otherwise, no update will be performed.

(3) Update ideal point: If $F_i\left(y\right)<z^i$, then $z^i=F_i\left(y\right)$.

Step3: Stop until the maximum number of iterations is reached. Otherwise, go back to step 2.

4. Improved MOEA/D Algorithm (IMOEA/D)

4.1. Algorithm Improvements

4.1.1. The Evolutionary Strategy of Competition between BMO and DE Operator

BMO is a primitive heuristic algorithm whose principle comes from the mating behavior of barnacles in nature [20]. Barnacles are hermaphroditic creatures that produce offspring in two ways: self-mating and normal mating. The choice of the barnacle mating method depends on the distance between paternal and maternal generations. As marine organisms, sperm and eggs can form fertilized eggs in water even if there is a certain distance between the father and the mother. However, when the distance between the father and the mother exceeds the $pl$ value, the barnacles will produce offspring by self-mating. The iterative steps are represented as follows:

$$\begin{array}{c}\mathit{I}\mathit{f}distance<pl\hfill \\ x_i^{N_{-}new}=p{x}_{barnacl{e}_{-}d}^N+q{x}_{barnacl{e}_{-}m}^N\hfill \\ \mathit{e}\mathit{l}\mathit{s}\mathit{e}\hfill \\ x_i^{n_{-}new}=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}()\times x_{\mathit{barnacle}\_m}^n\hfill \\ \mathit{E}\mathit{n}\mathit{d}\hfill \end{array}$$

where $p$ is a random number in the range 0 to 1, and $q=1-p$. $x_{barnacl{e}_{-}d}^n$ and $x_{barnacl{e}_{-}m}^n$ are the selected variables for paternal and maternal barnacles, respectively. It can be understood that $p$ and $q$ represent the percentage of the characteristics of the father and mother in the next generation, respectively. Thus, the offspring inherit the behavior of their father and mother according to a random number probability between 0 and 1. It is worth emphasizing that the value of $pl$ plays an important role in determining the process of exploration and exploitation. In this paper, $pl$ is set to 9.

The competitive evolution strategy: At the beginning of the algorithm, the DE algorithm is used to perform the iteration and record the number of updated individuals after each round of updates. When the number of updated individuals is reduced, it means the evolution speed slows down. At this time, the BMO algorithm replaces the DE algorithm to perform the iteration. When using BMO for iterative operations, if the number of updated individuals is reduced after a round of updates, the DE algorithm replaces the BMO algorithm to perform the iteration again. The two algorithms are used alternately to update the population through this mechanism, which improves the diversity of species.

4.1.2. Adaptive Mutation

In the process of population evolution, species diversity changes from high to low. The mutation rate is adjusted according to this rule. Figure 1a shows that the mutation rate changes exponentially with the number of iterations [21]. In the early stage of evolution, species diversity was better and the mutation rate was higher. In the process of evolution, the mutation rate increases continuously with the decrease of species diversity and finally tends to be stable. The probability of mutation is computed by Equation (15).

$$\mathrm{Muta}\_\mathrm{rate}=0.1\ast \left(-e^{-\frac{gen}{\gamma }}+1\right)$$

(15)

where $gen$ is the number of generations, and $\gamma$ is the scaling factor. Here, $\gamma$ is set to 500.

4.1.3. Selective Candidate with Similarity Selection (SCSS)

To prevent the algorithm from falling into local optima, SCSS is utilized to balance exploration and exploitation [22]. The rule is based on DE. It executes the DE algorithm twice to generate two candidate solutions—$u_1^i$ and $u_2^i$. Select the candidate solutions according to the following rules to complete the replacement of the current solution.

$$\begin{array}{c}\mathit{I}\mathit{f}rand()\times 2\times GD>\frac{rank\left(i\right)}{N}\hfill \\ \mathrm{Select}\mathrm{the}\mathrm{closest}\mathrm{candidate}\mathrm{from}\{u_1^i,u_2^i\}\mathrm{for}\mathrm{the}\mathrm{current}\mathrm{solution}{x}^i\\ \mathit{e}\mathit{l}\mathit{s}\mathit{e}\hfill \\ \mathrm{Select}\mathrm{the}\mathrm{farthest}\mathrm{candidate}\mathrm{from}\{u_1^i,u_2^i\}\mathrm{for}\mathrm{the}\mathrm{current}\mathrm{solution}{x}^i\\ \mathit{E}\mathit{n}\mathit{d}\hfill \end{array}$$

where $rank\left(i\right)$ is the solution quality ranking, and $N$ is the population size. High-quality solutions tend to be replaced by closer solutions, while the inferior solution prefers the farther one. $GD$ is a fixed value. The value of $GD$ is the key to the selection rule, which controls the balance between exploitation and exploration. The larger $GD$ is, the more current solutions select closer candidates, and the algorithm becomes more exploitative. On the contrary, the smaller $GD$ value is, the more current solutions select further candidate solutions, and the exploration capability of the algorithm is enhanced. After several attempts, we found a proper value for $GD$. When $GD$ is taken to be 0.5, the algorithm has a more balanced exploitation and exploration capability.

4.2. Superiority of Feasible Solutions

This paper adopts the superiority of feasible solutions to deal with many constraints in the OPF problem [23]. The total constraint violation is calculated by Equation (16).

$$\epsilon \left(x\right)=\frac{\sum _{k=1}^{N_C}{w}_k\left(U_k\left(x\right)\right)}{\sum _{k=1}^{N_C}{w}_k}$$

(16)

where $U_k\left(x\right)$ is the violation value of solution $x$ for all constraints. $w_k=1/U_{k,max}$, $U_{k,max}$ is the maximum value of $U_k\left(x\right)$, and $N_C$ is the number of constraints. The quality of the solution is judged by the degree of a constraint violation. The smaller the value of $\epsilon \left(x\right)$, the better the quality of the solution. $\epsilon \left(x\right)=0$ means no constraint violation, where $x$ is a feasible solution.

Given that $x_i$ and $x_j$ are any two solutions in the population, the quality of the solution $x_i$ is better than $x_j$ in the following cases:

(1) $\epsilon \left(x_i\right)<\epsilon \left(x_j\right)$

(2) $\epsilon \left(x_i\right)=\epsilon \left(x_j\right)$ $andg\left(x_i\right)<g\left(x_j\right)$

When a new individual is created, compare the quality between the new solution and the original solution according to the above rules and judge whether to update.

4.3. IMOEA/D

The IMOEA/D algorithm is based on the framework of the MOEA/D algorithm. The implementation steps of the IMOEA/D algorithm are as follows:

Step1: Initialization

Initialize $N$ individuals; each individual includes the following parameters:

(1) The solution is $x_i$, and its target value is $F\left(x_i\right)=\{f_1\left(x_i\right),f_2\left(x_i\right),\dots ,f_M\left(x_i\right)\}$,$(i\in \{1,2,\dots ,N\left\}\right)$; M is the number of the objectives.

(2) The weight vector is ${\lambda }^i=\{{\lambda }_1^i,{\lambda }_2^i,\dots ,{\lambda }_M^i\}$.

(3) The neighborhood is $B\left(i\right)$, $T$ presents the size of the neighborhood. $T$ is set to 30.

(4) The minimum value of the $i$-th objective function is $z_i^{*}$, the maximum value of the $i$-th objective function is $z_i^{\mathrm{nad}}$.

(5) The utility value for each individual is ${\pi }_i$.

(6) The total constraint violation is $\epsilon \left(x_i\right)$.

Step2: Selection of individuals and update of solutions

In each generation, about 20% of the total individuals are selected using the 10-tournament selection method [24]. The current solutions associated with the selected indices of individuals are picked up in the generation. The algorithm performs the following operations for each selected solution $x_i$ [25]:

(1) Calculate the mutation rate based on the current iteration number

(2) According to the evolutionary strategy, two individuals $x_j$ and $x_k$ are randomly selected from the neighborhood $B\left(i\right)$, and the BMO or DE operator is selected for cross-evolution based on the population evolution speed, resulting in a new individual $y_1$.

(3) Perform step (2) again to generate another candidate solution $y_2$ and select one solution from $\{y_1,y_2\}$, according to the similarity selection. The final selected new solution is named $y$.

(4) For the individual $x$, if $y$ produces fewer constraint violations than $x$, replace $x$ with y. If the constraint violation values for $x$ and $y$ are the same, compare the target fitness values of the related individuals, and the solution with the smaller target value will be retained. This step ensures that the new solutions will move toward the feasible region, and the improvement of the solution is completed when the violation constraint of the new individual gradually decreases and finally becomes “0”.

Step3: Update of individual utility value

The utility values of all individuals will be updated every 50 generations [24]. The individual utility value is updated according to Equations (17) and (18).

$${\pi }^m=\left\{\begin{array}{c}1if{\Delta }^m>0.001\hfill \\ \left(0.95+0.05\frac{{\Delta }^m}{0.001}\right){\pi }^{m}otherwise\end{array}\right.$$

(17)

$${\Delta }^m=\frac{g^{tch}\left(x_{t-\Delta t}^i\mid {\lambda }^i,z^{*}\right)-g^{tch}\left(x_t^i\mid {\lambda }^i,z^{*}\right)}{g^{tch}\left(x_{t-\Delta t}^i\mid \lambda ,z^{*}\right)}$$

(18)

where $t$ is the current generation, $g^{tch}$ is the Tchebycheff aggregation function, $x_{t-\Delta t}^i$ is the solution of individual $i$ at the generation of $t-\Delta t$ and $x_t^i$ is the one at the current generation $t.$ ${\Delta }^m$ represents the degree of change in the objective function.

Step4: Stop after reaching the maximum number of iterations. Otherwise, go back to step 2.

Table 1. Parameter settings for the IMOEA/D algorithm.

Symbol	Parameter	Value
$N$	Population size	200 for 2 objectives’ problem
		300 for 3 objectives’ problem
		455 for 4 objectives’ problem
$T$	Neighborhood size	0.1 × $N$
$n_r$	Max number of members replaced with better new solutions for each individual at every generation	0.01 × $N$
$Max$_$veal$	Maximum number of fitness evaluations	100,000
$CR$	Crossover rate	0.7
$pl$	Distance between paternal and maternal generations	9
$GD$	Greedy degree	0.5
${\Delta }^{req}$	Frequency of utility function update	50 generations

shows the parameters involved in the IMOEA/D algorithm when dealing with the OPF problem. Some of these parameters, such as population size, crossover rate, and the maximum number of fitness evaluations, are different from the conventional MOEA/D. Different population sizes are set for the cases with different target numbers. In this way, a uniform Pareto front can be obtained, and the computation time can be shortened. Greedy degree ($GD$), $pl$, and crossover rate ($CR$) values are chosen after a few trials. The optimization effect of the algorithm is related to the neighborhood, so it is necessary to maintain a high probability of neighbor updates. The algorithm is implemented by using MATLAB software, and a computer equipped with Intel Core i5 CPU @2.7 GHz and 4 GB RAM is used to run the simulations.

5. Results

To illustrate the effectiveness of IMOEA/D, in this paper, IEEE 30-bus and IEEE 57-bus systems are employed as the test cases in solving the multi-objective OPF problem. The optimization results are compared with other multi-objective algorithms. Except for MOEA/D-SF, the optimization results of other multi-objective optimization algorithms are directly quoted. The main parameter settings of MOEA/D-SF are the same as IMOEA/D.

5.1. Comparison of Multi-Objective Algorithm Indicator

The IGD metric [26] is the average of the minimum distance of all individuals in the actual Pareto front to the ideal Pareto front, which reflects the degree of convergence of the algorithm. The IGD definition is shown in Equation (19).

$$\mathit{IGD}\left(S,P^{*}\right)=\frac{\sum _{x\in P^{*}}\mathit{dist}\left(x,S\right)}{\left|P^{*}\right|}$$

(19)

where $P^{*}$ is the ideal Pareto front and $S$ is the actual Pareto front. $\mathit{dist}\left(x,S\right)$ is the Euclidean distance between individual $x$ on the ideal Pareto front $P^{*}$ to the nearest individual on $S$. $\left|P^{*}\right|$ normalizes the IGD value. The smaller the IGD value, the better the degree of convergence and diversity of $S.$ HV metric [27] is the volume of the target space enclosed by the set of non-dominated solutions of the algorithm and the reference points. The HV metric is calculated by Equation (20).

$$HV\left(S\right)=\mathrm{V}\mathrm{O}\mathrm{L}\left({\cup }_{x\in s}\left[f_1\left(x\right),r_1^{*}\right]\times \dots \left[f_m\left(x\right),r_m^{*}\right]\right)$$

(20)

where $r^{*}$ is a set of reference points in the target space, $r^{*}=(r_1^{*},r_2^{*},\dots ,r_m^{*})$, $r^{*}$ is dominated by all solutions in $S$. $\mathrm{V}\mathrm{O}\mathrm{L}\left(\right)$ indicates the Lebesgue measure [28]. The higher the HV value, the more $S$ is approximated to the entire Pareto front.

Table 2. Comparison of HV for study case 1-case 10.

Test cases	IMOEA/D			MOEA/D-SF
	Mean	MaxHV	Std	Mean	MaxHV	Std
Case 1	0.8274	0.8288	0.0011	0.672	0.7053	0.0245
Case 2	0.7857	0.8686	0.0909	0.7848	0.8956	0.1055
Case 3	0.6996	0.7013	0.0011	0.6986	0.7002	0.001
Case 4	0.8262	0.8306	0.0028	0.8212	0.8244	0.0023
Case 5	0.7792	0.7856	0.0035	0.7646	0.7681	0.0032
Case 6	0.6773	0.686	0.0043	0.6706	0.6858	0.0091
Case 7	0.8389	0.8509	0.0084	0.8258	0.8332	0.005
Case 8	0.8667	0.9475	0.0555	0.8484	0.9327	0.0647
Case 9	0.8479	0.8619	0.0075	0.8201	0.8568	0.0268
Case 10	0.7443	0.7581	0.0086	0.7436	0.7582	0.0133

Bold indicates the best values.

lists the mean, maximum, and standard deviations of HV over 10 runs for each case. Likewise,

Table 3. Comparison of IGD for study case 1-case 10.

Test cases	IMOEA/D			MOEA/D-SF
	Mean	MinIGD	Std	Mean	MinIGD	Std
Case 1	0.2221	0.2067	0.014	0.98	0.448	0.3667
Case 2	0.1	0.0095	0.0878	0.0756	0.0144	0.0999
Case 3	0.2889	0.2446	0.0366	0.3246	0.265	0.0365
Case 4	0.3618	0.2993	0.05	0.3762	0.2844	0.067
Case 5	0.4906	0.4134	0.0733	0.5281	0.4435	0.0756
Case 6	0.4691	0.3609	0.0647	0.6148	0.5196	0.0957
Case 7	7.2109	5.2896	2.0687	8.976	1.8343	3.8577
Case 8	8.624	6.5318	1.5779	15.0639	2.4823	20.5443
Case 9	7.4198	5.8653	1.3495	8.2303	5.9423	1.8555
Case 10	8.624	6.5318	1.5779	11.2508	8.914	1.579

Bold indicates the best values.

lists the minimum and average IGD values in 10 runs for each case, along with the standard deviation. It can be seen from Table 3, except in case 2 and case 10, that the maximum HV value of IMOEA/D is slightly lower than that of MOEA/D-SF. The HV value obtained by IMOEA/D is significantly higher than that of MOEA/D-SF. It can be seen from

Table 4. Comparison of best compromise solutions for study case 1.

Methods	Cost	Emission
IMOEA/D	828.674	0.2512
MOEA/D [29]	833.72	0.2438
MOEA/D-SF [15]	829.515	0.259
MOPSO [29]	833.86	0.2483
NSGA-II [29]	833.59	0.2449
MOMICA [30]	865.066	0.2221
ESDE [31]	833.474	0.254
ESDE-MC [31]	830.718	0.2483
ISPEA [6]	865.95	0.2234
GBICA [32]	830.852	0.2488
MGBICA [32]	830.851	0.2484

Bold indicates the best values.

that the average IGD of IMOEA/D is significantly lower than that of MOEA/D-SF in case 1 and cases 3–10. In cases 7 and 8, IMOEA/D achieves better minimum IGD values. In addition, the standard deviation of HV and IGD of IMOEA/D is better than that of MOEA/D-SF in most cases.

5.2. IEEE 30-Bus System: Detailed Results

Table 4,

Table 5. Comparison of best compromise solutions for study case 2.

Method	Cost	VD
IMOEA/D	802.374	0.1391
MOEA/D	799.99	0.354
MOEA/D-SF	802.406	0.1362
MOPSO	800.03	0.4422
NSGA-II	800.06	0.4486
MOMICA	804.96	0.0952

Bold indicates the best values.

,

Table 6. Comparison of best compromise solutions for study case 3.

Methods	Cost	Emission	Loss
IMOEA/D	878.652	0.2164	4.1257
MOEA/D	902.54	0.2107	3.4594
MOEA/D-SF	881.02	0.2164	4.4144
MOPSO	891.48	0.2144	3.9557
NSGA-II	903.79	0.2103	3.7917

Bold indicates the best values.

,

Table 7. Comparison of best compromise solutions for study case 4.

Methods	Cost	Emission	VD
IMOEA/D	836.650	0.2462	0.1153
MOEA/D	850.28	0.2332	0.1155
MOEA/D-SF	842.446	0.2406	0.1093
MOPSO	846.93	0.2386	0.2188
NSGA-II	825.86	0.2648	0.1421

Bold indicates the best values.

,

Table 8. Comparison of best compromise solutions for study case 5.

Methods	Cost	VD	Loss
IMOEA/D	835.8	0.1244	5.8176
MOEA/D	831.81	0.1355	5.9926
MOEA/D-SF	836.711	0.1258	5.8233
MOPSO	827.82	0.1588	6.5929
NSGA-II	843.14	0.1931	6.4917
B-MMOFPA [33]	843.18	0.1745	5.7886

Bold indicates the best values.

and

Table 9. Comparison of best compromise solutions for study case 6.

Methods	Cost	Emission	VD	Loss
IMOEA/D	869.260	0.2251	0.1214	4.7533
MOEA/D-SF	883.322	0.2187	0.1322	4.4527
MOMICA	830.188	0.2523	0.2978	5.585

Bold indicates the best values.

compare the best compromise solutions of IMOEA/D and other algorithms in different cases. In case 1, both cost and emission are optimized simultaneously. The best compromise solutions for IMOEA/D are 826.674 $/h and 0.2512 t/h. MSFLA has the best cost value but also gives the highest emission value of all algorithms. MOMICA provides the lowest emission value, but the cost is also significantly higher than other algorithms. In case 2, cost and voltage deviation are optimized at the same time. Although the cost values of MOEA/D, MOPSO, and NSGA-II are lower than IMOEA/D, the best compromise solutions provided by these three algorithms are infeasible. In 3-objective optimizations of cases 3, 4, and 5, according to the data in the table, IMOEA/D obtains better compromise solutions. Among the best compromise solutions of IMOEA/D, the cost objective is the lowest in case 3, while the voltage deviation is the least in cases 4 and 5. In the 4-objective case 6, IMOEA/D achieves the best value of cost and voltage deviation.

Figure 2a compares the Pareto fronts (PFs) of IMOEA/D and MOEA/D with the penalty factor approach (MOEA/D-penalty) in case 1. All PFs here represent the best PF with maximum HV value among all trial runs in a case study. In Figure 2a, the Pareto front of IMOEA/D is similar to that of MOEA/D-penalty. However, it is worthwhile that some solutions of MOEA/D-penalty are infeasible. Figure 2b shows that the Pareto solutions of MOEA/D-SF dominate the majority of the non-dominated solutions given by MOEA/D-penalty in case 3. In general, the Pareto front distribution of IMOEA/D is relatively uniform in most cases. Figure 3 compares the Pareto fronts of IMOEA/D and MOEA/D-penalty in cases 2, 5, and 6, respectively. As shown in Figure 3, IMOEA/D achieves PFs with better convergence and more uniform distribution, which confirms the effectiveness of the improved algorithm.

5.3. IEEE 57-Bus System: Detailed Results

Table 10. Comparison of best compromise solutions for study case 7.

Method	Cost	Emission
IMOEA/D	42,297.43	1.298
ISPEA	42,444.55	1.2904
NSGA-II	43,567.77	1.2979
ESDE	42,863.32	1.2662
ESDE-MC	42,857.49	1.2191
GBICA	42,138.37	1.3941
MGBICA	42,369.07	1.2940
MOMICA	41,886.8	1.4784
MOICA [30]	41,919.71	1.6010
MODFA [34]	43,174.57	1.2679

Bold indicates the best values.

,

Table 11. Comparison of best compromise solutions for study case 8.

Method	Cost	Emission
IMOEA/D	41,794.49	0.6068
MDE [35]	41,843	0.5962

Bold indicates the best values.

,

Table 12. Comparison of best compromise solutions for study case 9.

Method	Cost	Loss	VD
IMOEA/D	42,295.469	12.1544	0.6462
MDE	42,070	12.4024	0.6933
CMICA4 [36]	41,781.73	13.9936	0.8127
NSGA-II	41,930.94	21.5325	2.6699
MOPSO	41,901.36	16.8022	2.0059

Bold indicates the best values.

and

Table 13. Comparison of best compromise solutions for study case 10.

Method	Cost	Emission	VD	Loss
IMOEA/D	42,609.03	1.382	0.668	11.936
MOMICA	41,983.06	1.496	0.797	13.697
MOICA	41,998.57	1.7605	0.8748	13.335

Bold indicates the best values.

compare the best compromise solutions for various multi-objective algorithms. In case 7, MOMICA achieves the lowest cost but with higher emission. ESDE-MC [26] achieves the lowest emission level with a higher cost than some other comparable algorithms. While the cost value of IMOEA/D is better than most of the comparison algorithms, the emission value is lower than GBICA, MOICA, and MOMICA. In case 8, IMOEA/D achieves better cost. In case 9, the best values of voltage deviation and loss are given by IMOEA/D. In 4-objective optimization case 10, although the cost value of IMOEA/D is the highest, the other three objective function values of emission, voltage deviation and power loss are the lowest.

Figure 4 shows the comparison of PFs for IMOEA/D and MOEA/D-penalty in cases 7 and 9. Similar to case 1, the best PF of the two algorithms almost coincides with case 7. Considering constraint violations, IMOEA/D performs slightly better. It can be seen from Figure 4 that the diversity and the convergence ability of IMOEA/D are stronger than MOEA/D-penalty. In case 9 of the 3-objective optimization, IMOEA/D achieves a set of Pareto solutions with more regular distribution, which further demonstrates the effectiveness of the algorithm improvement.

6. Conclusions

The OPF problem is a complex nonlinear optimization problem [35,36]. In this paper, an improved MOEA/D (IMOEA/D) algorithm is proposed for solving the multi-objective OPF problem. Despite the applicability of the MOEA/D in solving complicated problems, its performance is degraded when the dimension size of the OPF’s test system is increased. In this respect, the evolution strategy of MOEA/D is improved by introducing three improvement strategies: BMO and DE competitive evolution, adaptive mutation strategy, and similarity selection. The main purpose of IMOEA/D is to improve convergence speed and maintain population diversity. The effectiveness of the proposed IMOEA/D algorithm was experimentally evaluated using standard IEEE 30-bus and IEEE 57-bus test systems to optimize multi-objective functions of the OPF under the system constraints. At the same time, the optimization results are also compared with classical multi-objective optimization algorithms. The comparison of the simulation results shows that the IMOEA/D algorithm can provide better solutions than other comparative algorithms in solving the multi-objective OPF problem.