# Application of Hybrid MOPSO Algorithm to Optimal Reactive Power Dispatch Problem Considering Voltage Stability

Authors: Yujiao Zeng; Yanguang Sun
Journal: Journal of Electrical and Computer Engineering
 (2014)

---

## Abstract

This study presents a novel hybrid multiobjective particle swarm optimization (HMOPSO) algorithm to solve the optimal reactive power dispatch (ORPD) problem. This problem is formulated as a challenging nonlinear constrained multiobjective optimization problem considering three objectives, that is, power losses minimization, voltage profile improvement, and voltage stability enhancement simultaneously. In order to attain better convergence and diversity, this work presents the use of combing the classical MOPSO with Gaussian probability distribution, chaotic sequences, dynamic crowding distance, and self-adaptive mutation operator. Moreover, multiple effective strategies, such as mixed-variable handling approach, constraint handling technique, and stopping criteria, are employed. The effectiveness of the proposed algorithm for solving the ORPD problem is validated on the standard IEEE 30-bus and IEEE 118-bus systems under nominal and contingency states. The obtained results are compared with classical MOPSO, nondominated sorting genetic algorithm (NSGA-II), multiobjective evolutionary algorithm based on decomposition (MOEA/D), and other methods recently reported in the literature from the point of view of Pareto fronts, extreme, solutions and multiobjective performance metrics. The numerical results demonstrate the superiority of the proposed HMOPSO in solving the ORPD problem while strictly satisfying all the constraints.


---

## Body

## 1. Introduction

The optimal reactive power problem (ORPD) has attracted great attention in the past decades because it can greatly improve economy and security of power system. Generally, the aim of the ORPD is to minimize the network real power loss and improve voltage profile by regulating generator bus voltages, switching on/off static VAR compensators, and changing transformer tap-settings, while satisfying various constraints.Voltage stability is a major aspect of power system security analysis. It is well known that voltage instability and collapse have led to major system failures. Due to the continuous growth in the demand of electricity with unmatched generation and transmission capacity expansion, power systems are operated nearer against the voltage stability limit than before. Voltage instability is emerging as a new challenge to power system planning and operation. Hence, there is a need to consider voltage stability enhancement as one of the objectives in the ORPD. Voltage stability analysis often requires examination of a wide range of system conditions and a large number of contingency scenarios. For such application, many analysis methods of system voltage stability determination have been developed on static analysis techniques based on the power flow model, since these techniques are simple, fast, and convenient to use. In this work, a fast indicator of voltage stability,L-index, presented by Kessel and Glavitsch [1], is chosen as the objective for the voltage stability enhancement.In the literature, a large number of optimization techniques have been proposed to solve the ORPD problem. Many conventional methods such as nonlinear programming (NLP) [2], quadratic programming [3], linear programming (LP) [4], and interior point algorithms [5] have been adopted to solve the problem. However, these techniques are failed in handling nonconvexities and nonsmoothness and susceptible to be trapped in local minima. In order to overcome these limitations, a wide variety of the heuristic methods have been proposed to solve the ORPD problem, such as genetic algorithms (GA) [6], evolutionary algorithms (EP) [7], particle Swarm optimization (PSO) [8], differential evolution (DE) [9], seeker optimization algorithm (SOA) [10], gravitational search algorithm (GSA) [11], and biogeography based optimization (BBO) [12]. But these heuristic methods only result in a single optimal solution. In order to provide a means to assess trade-off between two conflicting objectives, the ORPD problem is formulated as a multiobjective problem, which is solved as a single objective problem by weighted sum of objectives. In order to obtain desired Pareto-optimal solutions, these methods require multiple runs and tend to find weakly nondominated solutions.In the last few years, the use of evolutionary algorithms for multiobjective optimization has significantly grown. These multiobjective evolutionary algorithms (MOEAs) can find multiple Pareto-optimal solutions in one single run. Nondominated sorting genetic algorithm II (NSGA-II) [13], strength Pareto evolutionary algorithm 2 (SPEA2) [14], Pareto archive evolution strategy II (PAES-II) [15], and multiobjective evolutionary algorithm based on decomposition (MOEA/D) [16] are the representatives of the state of the art in this area.Particle Swarm optimization (PSO) [17] is proposed by Kennedy and Eberhart, which is bioinspired metaheuristic that simulates the movement of a flock of birds or fish which seek food. Since this technique is relatively simple to implement, it has received widespread attention. Also, it has been extended to deal with multiple objectives by Coello Coello et al. [18] and other researchers [19, 20].Recently, some of these multiobjective optimization algorithms have been proposed to the MORPD problem. Strength Pareto evolutionary algorithm (SPEA) [21], multiple evolutionary algorithms with adaptive selection strategy (MEAASS) [22], and modified DE (MDE) [23] have been applied to solve a multiobjective ORPD model minimizing network real power loss and voltage deviation of the system, simultaneously. In [24, 25], voltage stability index and real power loss were considered as two objective functions and NSGA-II was utilized as an optimization technique. A multiobjective teaching learning algorithm based on decomposition (MOTLA/D) [26] was also utilized to solve the same biobjective problem. In [27], three objectives including real transmission losses, voltage deviation, and L-index were considered to formulate multiobjective ORPD model and an opposition-based self-adaptive modified gravitational search algorithm (OSAMGSA) was proposed as a solution technique. In [28, 29], MOPSO was utilized to solve the multiobjective optimal power flow problem by considering generation cost, transmissions loss, voltage deviations, and other functions as objectives.In this work, a new hybrid MOPSO algorithm (HMOPSO) is proposed for solving the multiobjective ORPD problem, in which three competing objectives such as real power loss minimization, voltage profile improvement, and voltage stability enhancement are optimized simultaneously in a single run. MOPSO is one of the modern heuristic algorithms that has already been applied successfully to solve complicated problems. With few parameters, MOPSO is simple to implement and efficient and robust compared to other algorithms. However, the classical MOPSO often converges to local optimal solutions so as to be premature convergence. Moreover, it suffers the problem of maintaining good diversity and uniform spacing among the solutions in the Pareto front. Therefore, some modifications have been incorporated into the classical MOPSO to improve its performance. In order to enhance overall stochastic search capability of algorithm and prevent premature convergence, the proposed approach employs Gaussian probability distribution to generate random numbers, chaotic sequences to tune the inertia weight, and a time variant structure for adjusting acceleration coefficients. In addition, a self-adaptive mutation operator is adopted to enhance the diversity of the population of MOPSO, which can avoid getting trapped in local optima. In order to obtain Pareto front with good diversity, dynamic crowding distance is adopted. The performance of the proposed approach is tested and evaluated on the standard IEEE 30-bus and IEEE 118-bus systems in both nominal and contingency situations. To validate the performance of the proposed algorithm, three performance metrics, namely, generational distance (GD) [30], minimal spacing (MSP) [31], and hypervolume (HV) [32], are used to compare the performance of the proposed HMOPSO with MOPSO, NSGA-II, MOEA/D, and other methods recently reported in the literature. The numerical results show the superiority of HMOPSO for handling the ORPD problem.


## 2. Problem Formulation

### 2.1. Problem Objectives

#### 2.1.1. Real Power Loss(PL)

The objective is to minimize the real power loss in transmission lines expressed as(1)min⁡PL=min⁡∑k=1NEGk[Vi2+Vj2-2ViVjcos⁡(δi-δj)],
where NE is the number of transmission lines; Gk is the conductance of the kth line; Vi∠δi and Vj∠δj are the voltages at end bus i and bus j of the kth line, respectively.


#### 2.1.2. Voltage Deviation(VD)

The objective is to minimize the deviations in voltage magnitude at load buses that can be expressed as(2)min⁡⁡VD=min⁡∑i=1NL|Vi-Vispec|,
where Vispec is the prespecified reference value at load bus i, which is usually set as 1.0 p.u., and NL is the number of load buses.


#### 2.1.3. Voltage Stability Index(Lmax⁡)

In this work, enhancing voltage stability can be achieved through minimizing the voltage stability indicatorL-index [1] values at every bus of the system and consequently the global power system indicator. This objective can be expressed as
(3)min⁡F3=min⁡Lmax⁡,
where Lmax⁡ is the global system indicator describing the stability of the complete system. The L-index calculation for a power system is briefly discussed below.For multinode system,(4)Ibus=Ybus×Vbus.By segregating the load buses (PQ) from generator buses (PV), (4) can be written as
(5)[IGIL]=[YGGYGLYLGYLL][VGVL],
where IG, IL  and VG, VL represent currents and voltages at generator buses and load buses.Rearranging the above equation, we get(6)[VLIG]=[ZLLFLGKLGYGG][ILVG],
where
(7)FLG=-[YLL]-1[YLG].TheL-index of the jth node is given by the expression
(8)Lj=|1-∑i=1NGFjiVGiVLj|.Thus, a global system indicator describing the stability of the complete system is(9)Lmax⁡=max⁡(Lj),j=1,2,…,NL.


### 2.2. Problem Constraints

#### 2.2.1. Equality Constraints

These constraints represent load flow equation as(10)PGi-PDi-Vi∑j=1NBVj[Gijcos⁡(δi-δj)+Bijsin(δi-δj)]=0,QGi-QDi-Vi∑j=1NBVj[Gijsin(δi-δj)+Bijcos⁡(δi-δj)]=0,
where i=1,2,…,NB; NB is the number of buses; PGi  and QGi are the generator real and reactive power, respectively; PDi and QDi are the load real and reactive power, respectively; Gij and Bij are the transfer conductance and susceptance between i and j, respectively.


#### 2.2.2. Inequality Constraints

The inequality constraints represent the system operating constraints:(11)VGimin⁡≤VGi≤VGimax⁡,i=1,2,…,NG,Timin⁡≤Ti≤Timax⁡,i=1,2,…,NT,QCimin⁡≤QCi≤QCimax⁡,i=1,2,…,NC,VPQimin⁡≤VPQi≤VPQimax⁡,i=1,2,…,NPQ,QGimin⁡≤QGi≤QGimax⁡,i=1,2,…,NG,SLk≤SLkmax⁡,k=1,2,…,NE,
where VGi and QGi are voltage magnitude and reactive power output at the ith generator bus, respectively; Ti is the tap setting of the ith transformer; QCi  is reactive power compensation of the ith switched VAR source; VPQi  is the PQ bus voltage; SLk is the apparent power flow at the kth transmission line. The max and min represent the maximum and minimum values of the corresponding variables, respectively. NG, NT, NC, and NPQ are the numbers of generators, transformers, switchable VAR sources, and PQ buses, respectively.


## 2.1. Problem Objectives

### 2.1.1. Real Power Loss(PL)

The objective is to minimize the real power loss in transmission lines expressed as(1)min⁡PL=min⁡∑k=1NEGk[Vi2+Vj2-2ViVjcos⁡(δi-δj)],
where NE is the number of transmission lines; Gk is the conductance of the kth line; Vi∠δi and Vj∠δj are the voltages at end bus i and bus j of the kth line, respectively.


### 2.1.2. Voltage Deviation(VD)

The objective is to minimize the deviations in voltage magnitude at load buses that can be expressed as(2)min⁡⁡VD=min⁡∑i=1NL|Vi-Vispec|,
where Vispec is the prespecified reference value at load bus i, which is usually set as 1.0 p.u., and NL is the number of load buses.


### 2.1.3. Voltage Stability Index(Lmax⁡)

In this work, enhancing voltage stability can be achieved through minimizing the voltage stability indicatorL-index [1] values at every bus of the system and consequently the global power system indicator. This objective can be expressed as
(3)min⁡F3=min⁡Lmax⁡,
where Lmax⁡ is the global system indicator describing the stability of the complete system. The L-index calculation for a power system is briefly discussed below.For multinode system,(4)Ibus=Ybus×Vbus.By segregating the load buses (PQ) from generator buses (PV), (4) can be written as
(5)[IGIL]=[YGGYGLYLGYLL][VGVL],
where IG, IL  and VG, VL represent currents and voltages at generator buses and load buses.Rearranging the above equation, we get(6)[VLIG]=[ZLLFLGKLGYGG][ILVG],
where
(7)FLG=-[YLL]-1[YLG].TheL-index of the jth node is given by the expression
(8)Lj=|1-∑i=1NGFjiVGiVLj|.Thus, a global system indicator describing the stability of the complete system is(9)Lmax⁡=max⁡(Lj),j=1,2,…,NL.


## 2.1.1. Real Power Loss(PL)

The objective is to minimize the real power loss in transmission lines expressed as(1)min⁡PL=min⁡∑k=1NEGk[Vi2+Vj2-2ViVjcos⁡(δi-δj)],
where NE is the number of transmission lines; Gk is the conductance of the kth line; Vi∠δi and Vj∠δj are the voltages at end bus i and bus j of the kth line, respectively.


## 2.1.2. Voltage Deviation(VD)

The objective is to minimize the deviations in voltage magnitude at load buses that can be expressed as(2)min⁡⁡VD=min⁡∑i=1NL|Vi-Vispec|,
where Vispec is the prespecified reference value at load bus i, which is usually set as 1.0 p.u., and NL is the number of load buses.


## 2.1.3. Voltage Stability Index(Lmax⁡)

In this work, enhancing voltage stability can be achieved through minimizing the voltage stability indicatorL-index [1] values at every bus of the system and consequently the global power system indicator. This objective can be expressed as
(3)min⁡F3=min⁡Lmax⁡,
where Lmax⁡ is the global system indicator describing the stability of the complete system. The L-index calculation for a power system is briefly discussed below.For multinode system,(4)Ibus=Ybus×Vbus.By segregating the load buses (PQ) from generator buses (PV), (4) can be written as
(5)[IGIL]=[YGGYGLYLGYLL][VGVL],
where IG, IL  and VG, VL represent currents and voltages at generator buses and load buses.Rearranging the above equation, we get(6)[VLIG]=[ZLLFLGKLGYGG][ILVG],
where
(7)FLG=-[YLL]-1[YLG].TheL-index of the jth node is given by the expression
(8)Lj=|1-∑i=1NGFjiVGiVLj|.Thus, a global system indicator describing the stability of the complete system is(9)Lmax⁡=max⁡(Lj),j=1,2,…,NL.


## 2.2. Problem Constraints

### 2.2.1. Equality Constraints

These constraints represent load flow equation as(10)PGi-PDi-Vi∑j=1NBVj[Gijcos⁡(δi-δj)+Bijsin(δi-δj)]=0,QGi-QDi-Vi∑j=1NBVj[Gijsin(δi-δj)+Bijcos⁡(δi-δj)]=0,
where i=1,2,…,NB; NB is the number of buses; PGi  and QGi are the generator real and reactive power, respectively; PDi and QDi are the load real and reactive power, respectively; Gij and Bij are the transfer conductance and susceptance between i and j, respectively.


### 2.2.2. Inequality Constraints

The inequality constraints represent the system operating constraints:(11)VGimin⁡≤VGi≤VGimax⁡,i=1,2,…,NG,Timin⁡≤Ti≤Timax⁡,i=1,2,…,NT,QCimin⁡≤QCi≤QCimax⁡,i=1,2,…,NC,VPQimin⁡≤VPQi≤VPQimax⁡,i=1,2,…,NPQ,QGimin⁡≤QGi≤QGimax⁡,i=1,2,…,NG,SLk≤SLkmax⁡,k=1,2,…,NE,
where VGi and QGi are voltage magnitude and reactive power output at the ith generator bus, respectively; Ti is the tap setting of the ith transformer; QCi  is reactive power compensation of the ith switched VAR source; VPQi  is the PQ bus voltage; SLk is the apparent power flow at the kth transmission line. The max and min represent the maximum and minimum values of the corresponding variables, respectively. NG, NT, NC, and NPQ are the numbers of generators, transformers, switchable VAR sources, and PQ buses, respectively.


## 2.2.1. Equality Constraints

These constraints represent load flow equation as(10)PGi-PDi-Vi∑j=1NBVj[Gijcos⁡(δi-δj)+Bijsin(δi-δj)]=0,QGi-QDi-Vi∑j=1NBVj[Gijsin(δi-δj)+Bijcos⁡(δi-δj)]=0,
where i=1,2,…,NB; NB is the number of buses; PGi  and QGi are the generator real and reactive power, respectively; PDi and QDi are the load real and reactive power, respectively; Gij and Bij are the transfer conductance and susceptance between i and j, respectively.


## 2.2.2. Inequality Constraints

The inequality constraints represent the system operating constraints:(11)VGimin⁡≤VGi≤VGimax⁡,i=1,2,…,NG,Timin⁡≤Ti≤Timax⁡,i=1,2,…,NT,QCimin⁡≤QCi≤QCimax⁡,i=1,2,…,NC,VPQimin⁡≤VPQi≤VPQimax⁡,i=1,2,…,NPQ,QGimin⁡≤QGi≤QGimax⁡,i=1,2,…,NG,SLk≤SLkmax⁡,k=1,2,…,NE,
where VGi and QGi are voltage magnitude and reactive power output at the ith generator bus, respectively; Ti is the tap setting of the ith transformer; QCi  is reactive power compensation of the ith switched VAR source; VPQi  is the PQ bus voltage; SLk is the apparent power flow at the kth transmission line. The max and min represent the maximum and minimum values of the corresponding variables, respectively. NG, NT, NC, and NPQ are the numbers of generators, transformers, switchable VAR sources, and PQ buses, respectively.


## 3. Proposed Approach

### 3.1. Brief View of Multiobjective Particle Swarm

Up to now, there have been several proposals to extend PSO to handle multiobjective optimization problems. Here, we choose one of popular MOPSO algorithms developed by Coello Coello et al. [18] to briefly explain how to extend PSO to the multiobjective optimization (MOO) problem.In PSO, the group is a community composed of all the particles, which fly around in a multidimensional search space. During each iteration, each particle adjusts its position according to its past experiences (the best local positionPbesti) and its neighbors’ experiences (the best global position Gbest). To extend PSO to handle multiobjectives, Coello uses the concept of Pareto dominance to decide the Pbesti and Gbest. The Pbesti of the i particle is the position which is nondominated by its own past positions. To determine the best global position Gbest for each particle in the population, Coello adopts an external (or secondary) repository to store the nondominated vectors found along the search process. Thus, the adjustment for each particle during each flight iteration can be expressed as follows [18]:
(12)vi(t+1)=w·vi(t)+c1·r1·[Pbesti(t)-xi(t)]+c2·r2·[Repj(t)-xi(t)],(13)xi(t+1)=xi(t)+vi(t+1),
where t is the index of iterations, xi is the position of the particle i, vi is the velocity of the particle i, w is the inertia weight, c1 and c2 are the acceleration constants, r1 and r2 are random numbers based on uniform probability distribution function in the range [0,1], Repj is the best global value of the particle i that is taken from the repository, and the index j is determined through two levels using roulette-wheel selection method and random selection method.Though MOPSO has been applied to the MOO problem with impressive success, it suffers from premature convergence and bad diversity of the Pareto front. In this paper, MOPSO is improved by incorporating Gaussian probability distribution, chaotic descending inertia weight, time variant acceleration coefficients, self-adaptive mutation operator, and dynamic crowding distance, which are explained below in detail.


### 3.2. Gaussian Distributed Random Numbers

In this paper, we employ Gaussian distribution sequences to generate the random numbersr1 and r2 instead of uniform probability distribution in the original MOPSO. This approach can allow particles to move away from the current point and escape from local minima. Thus, the velocity update equation (12) is modified as follows:
(14)vi(t+1)=w·vi(t)+c1·G1·[Pbesti(t)-xi(t)]+c2·G2·[Repj(t)-xi(t)],
where G1 and G2 are random numbers generated by Gaussian probability distribution with zero mean and unit variance.


### 3.3. Chaotic Inertia Weight

In MOPSO, the parameterw controls the influence of the precious velocity on the present velocity. Suitable selection of w can provide a balance between global exploration and local exploitation and results in a lower number of iterations to find the optimal solution. Thus, Shi and Eberhart [33] proposed a time variant inertia weight (IWA), linearly descending with iteration from wmax⁡ to wmin⁡ expressed as follows:
(15)wt=wmax⁡+(wmin⁡-wmax⁡)·ttmax⁡.This strategy can enhance the convergence velocity of MOPSO, but it pays out cost: MOPSO usually gets into the local optimum when it solves the question of more apices’ function. In order to overcome these draws, this paper implements chaotic inertial weight approach (CIWA) defined as follows:(16)cwt=wt×γt,
where cwt is a chaotic weight at iteration t, wt is the weight factor from the IWA, and γt is the chaotic parameter variable at iteration t. Here, γt is based on the Logistic map, whose formula is given by
(17)γt=u×γt-1×(1-γt-1),
where u is a control parameter and γt-1 denotes chaotic variable at iteration t-1. Although the above equation is deterministic, it exhibits chaotic dynamics when u=4, the initial γ0∈(0,1), and γ0≠{0,0.25,0.5,0.75,1.0}.


### 3.4. Time Variant Acceleration Coefficients

In addition tow, parameters c1 and c2 influence the performance of PSO. Higher values of c1 ensure larger deviation of the particle in the search space; the higher values of c2 signify the convergence to the present global best. To incorporate better compromise between the exploration and exploitation of the search space, time variant acceleration coefficients introduced in [34] are employed. The value of c1 is allowed to decrease linearly with iteration from c1,i to c1,f, and the value of c2 is allowed to increase linearly with iteration from c2,i to c2,f described as follows:
(18)c1,t=c1,i+(c1,f-c1,i)·ttmax⁡,c2,t=c2,i+(c2,f-c2,i)·ttmax⁡.


### 3.5. Self-Adaptive Mutation

In order to diversify the population of MOPSO and avoid being trapped in local optimal solutions, a self-adaptive mutation strategy is newly introduced to the algorithm. The mechanism mixes several different mutation functions; each mutation function has different searching capability. This feature causes searching of the entire search space in different directions to be improved. Consequently, the diversity of the swarm increases intensely and the probability of being trapped in local minima decreases drastically; that will insure the global optimal value being found by the proposed HMOPSO. In this mutation technique, four different mutation functions are put forward to generate different mutant vectors described as follows:(19)xmutant,1(t)=xr1(t)+Y1×(Gbest(t)-xr2(t))+(1-Y1)*(xr3(t)-xr4(t)),xmutant,2(t)=xr1(t)+Y2×(xr2(t)-xr3(t))+(1-Y2)×(xr4(t)-xr1(t)),xmutant,3(t)=Gbest(t)+Y3×(xr1(t)-xr2(t)),xmutant,4(t)=Gbest(t)+0.5×Y4×Gbest(t),
where r1, r2, r3, and r4 are random integers uniformly between 1 and the number of the population of the swarm and r1≠r2≠r3≠r4; in other words, the indices are mutually different; Y1, Y2, Y3, and Y4 are uniform distributed random numbers in the range [0,1].


### 3.6. Dynamic Crowding Distance

In order to improve the diversity of the external archive, this paper employs a diversity maintenance strategy (DMS) based on dynamic crowding distance (DCD) [35], a revised version of crowding distance (CD) proposed in NSGA-II. The DCD of the individual is given as follows:
(20)DCDi=CDilog⁡(1/Vari),
where
(21)CDi=1Nobj∑m=1Nobj|Ii+1m-Ii-1m|,Vari=1Nobj∑m=1Nobj(|Ii+1m-Ii-1m|-CDi)2,
where Nobj is the number of objectives; Iim, Imax⁡m, and Imin⁡m are the mth objective of the ith individual, of the maximum individual, and of the minimal individual, respectively, after sorting the population according to the mth objective function values; CDi is the CD of the ith individual; Var⁡i is the variance of individuals which are neighbors of the ith individual.The application of DCD to maintain the diversity of the external archive is explained below. At each iteration, the archive gets updated by incorporating nondominated solutions from the current population into archive. If the archive sizeM exceeds the maximum size limit Na, the archive truncation procedure is invoked, which can be described as follows.(1)
Calculate DCDs of particles in the archiveAt using (20).(2)
Sort the nondominated solutions inAt according to their DCDs in descending order.(3)
Remove the individual with the lowest DCD value fromAt.(4)
IfM≤Na, stop the archive truncation; otherwise go to step 2 and continue.


## 3.1. Brief View of Multiobjective Particle Swarm

Up to now, there have been several proposals to extend PSO to handle multiobjective optimization problems. Here, we choose one of popular MOPSO algorithms developed by Coello Coello et al. [18] to briefly explain how to extend PSO to the multiobjective optimization (MOO) problem.In PSO, the group is a community composed of all the particles, which fly around in a multidimensional search space. During each iteration, each particle adjusts its position according to its past experiences (the best local positionPbesti) and its neighbors’ experiences (the best global position Gbest). To extend PSO to handle multiobjectives, Coello uses the concept of Pareto dominance to decide the Pbesti and Gbest. The Pbesti of the i particle is the position which is nondominated by its own past positions. To determine the best global position Gbest for each particle in the population, Coello adopts an external (or secondary) repository to store the nondominated vectors found along the search process. Thus, the adjustment for each particle during each flight iteration can be expressed as follows [18]:
(12)vi(t+1)=w·vi(t)+c1·r1·[Pbesti(t)-xi(t)]+c2·r2·[Repj(t)-xi(t)],(13)xi(t+1)=xi(t)+vi(t+1),
where t is the index of iterations, xi is the position of the particle i, vi is the velocity of the particle i, w is the inertia weight, c1 and c2 are the acceleration constants, r1 and r2 are random numbers based on uniform probability distribution function in the range [0,1], Repj is the best global value of the particle i that is taken from the repository, and the index j is determined through two levels using roulette-wheel selection method and random selection method.Though MOPSO has been applied to the MOO problem with impressive success, it suffers from premature convergence and bad diversity of the Pareto front. In this paper, MOPSO is improved by incorporating Gaussian probability distribution, chaotic descending inertia weight, time variant acceleration coefficients, self-adaptive mutation operator, and dynamic crowding distance, which are explained below in detail.


## 3.2. Gaussian Distributed Random Numbers

In this paper, we employ Gaussian distribution sequences to generate the random numbersr1 and r2 instead of uniform probability distribution in the original MOPSO. This approach can allow particles to move away from the current point and escape from local minima. Thus, the velocity update equation (12) is modified as follows:
(14)vi(t+1)=w·vi(t)+c1·G1·[Pbesti(t)-xi(t)]+c2·G2·[Repj(t)-xi(t)],
where G1 and G2 are random numbers generated by Gaussian probability distribution with zero mean and unit variance.


## 3.3. Chaotic Inertia Weight

In MOPSO, the parameterw controls the influence of the precious velocity on the present velocity. Suitable selection of w can provide a balance between global exploration and local exploitation and results in a lower number of iterations to find the optimal solution. Thus, Shi and Eberhart [33] proposed a time variant inertia weight (IWA), linearly descending with iteration from wmax⁡ to wmin⁡ expressed as follows:
(15)wt=wmax⁡+(wmin⁡-wmax⁡)·ttmax⁡.This strategy can enhance the convergence velocity of MOPSO, but it pays out cost: MOPSO usually gets into the local optimum when it solves the question of more apices’ function. In order to overcome these draws, this paper implements chaotic inertial weight approach (CIWA) defined as follows:(16)cwt=wt×γt,
where cwt is a chaotic weight at iteration t, wt is the weight factor from the IWA, and γt is the chaotic parameter variable at iteration t. Here, γt is based on the Logistic map, whose formula is given by
(17)γt=u×γt-1×(1-γt-1),
where u is a control parameter and γt-1 denotes chaotic variable at iteration t-1. Although the above equation is deterministic, it exhibits chaotic dynamics when u=4, the initial γ0∈(0,1), and γ0≠{0,0.25,0.5,0.75,1.0}.


## 3.4. Time Variant Acceleration Coefficients

In addition tow, parameters c1 and c2 influence the performance of PSO. Higher values of c1 ensure larger deviation of the particle in the search space; the higher values of c2 signify the convergence to the present global best. To incorporate better compromise between the exploration and exploitation of the search space, time variant acceleration coefficients introduced in [34] are employed. The value of c1 is allowed to decrease linearly with iteration from c1,i to c1,f, and the value of c2 is allowed to increase linearly with iteration from c2,i to c2,f described as follows:
(18)c1,t=c1,i+(c1,f-c1,i)·ttmax⁡,c2,t=c2,i+(c2,f-c2,i)·ttmax⁡.


## 3.5. Self-Adaptive Mutation

In order to diversify the population of MOPSO and avoid being trapped in local optimal solutions, a self-adaptive mutation strategy is newly introduced to the algorithm. The mechanism mixes several different mutation functions; each mutation function has different searching capability. This feature causes searching of the entire search space in different directions to be improved. Consequently, the diversity of the swarm increases intensely and the probability of being trapped in local minima decreases drastically; that will insure the global optimal value being found by the proposed HMOPSO. In this mutation technique, four different mutation functions are put forward to generate different mutant vectors described as follows:(19)xmutant,1(t)=xr1(t)+Y1×(Gbest(t)-xr2(t))+(1-Y1)*(xr3(t)-xr4(t)),xmutant,2(t)=xr1(t)+Y2×(xr2(t)-xr3(t))+(1-Y2)×(xr4(t)-xr1(t)),xmutant,3(t)=Gbest(t)+Y3×(xr1(t)-xr2(t)),xmutant,4(t)=Gbest(t)+0.5×Y4×Gbest(t),
where r1, r2, r3, and r4 are random integers uniformly between 1 and the number of the population of the swarm and r1≠r2≠r3≠r4; in other words, the indices are mutually different; Y1, Y2, Y3, and Y4 are uniform distributed random numbers in the range [0,1].


## 3.6. Dynamic Crowding Distance

In order to improve the diversity of the external archive, this paper employs a diversity maintenance strategy (DMS) based on dynamic crowding distance (DCD) [35], a revised version of crowding distance (CD) proposed in NSGA-II. The DCD of the individual is given as follows:
(20)DCDi=CDilog⁡(1/Vari),
where
(21)CDi=1Nobj∑m=1Nobj|Ii+1m-Ii-1m|,Vari=1Nobj∑m=1Nobj(|Ii+1m-Ii-1m|-CDi)2,
where Nobj is the number of objectives; Iim, Imax⁡m, and Imin⁡m are the mth objective of the ith individual, of the maximum individual, and of the minimal individual, respectively, after sorting the population according to the mth objective function values; CDi is the CD of the ith individual; Var⁡i is the variance of individuals which are neighbors of the ith individual.The application of DCD to maintain the diversity of the external archive is explained below. At each iteration, the archive gets updated by incorporating nondominated solutions from the current population into archive. If the archive sizeM exceeds the maximum size limit Na, the archive truncation procedure is invoked, which can be described as follows.(1)
Calculate DCDs of particles in the archiveAt using (20).(2)
Sort the nondominated solutions inAt according to their DCDs in descending order.(3)
Remove the individual with the lowest DCD value fromAt.(4)
IfM≤Na, stop the archive truncation; otherwise go to step 2 and continue.


## 4. Approach Implementation

### 4.1. Mixed-Variable Handling Method

In the decision variables of the ORPD problem, generator voltages are continuous variables, which can be running at any real number within the limit boundary; transformer taps and reactive compensations are discrete variables, which can only be given a value from a fixed discrete values set. Here, we use real numbers to represent continuous variables and the integer to other variables. The vector of new control variables can be written as(22)X*T=[vG1⋯vGNG,t1⋯tNT,qc1⋯qcNC],
where  vGi, ti, and  qCi are the generator bus voltage, the position of transformer tap, and the step of reactive compensation capacitor, respectively. Then, the relationship between practical parameters and new control variables can be written as
(23)VGi=vGi,Ti=1.0+ti·Δt,QCi=qCi·ΔqCi,
where Δti and ΔqCi express the ratio of transformer tap and the step capability of reactive compensation capacitor.Because the basic forms of MOPSO can only handle continuous variables, we employ a mixed-variable handling method [36]. In the initialization process, all variables are randomly generated within their upper and lower bounds at first. Then, the integer part of the variable value is picked to be its value, denoted as int⁡(xi(t)). In the iterative process, the new position for the integer variables will be selected in the neighbor integer values of the former position according to the velocity of the particle. If the velocity is more than zero, the new position is forward to plus 1. Otherwise, the position is back to minus 1.


### 4.2. Constraints Handling

From Section2, we know that the ORPD problem is subjected to several equality and inequality constraints. The equality constraints are satisfied by running Newton-Raphson method; the inequality constraints on the control variables can be automatically satisfied by setting the boundary of the search space; only the inequality constraints on the state variables need to be considered. In order to handle these constraints without introducing heavy work of turning penalty factors into the HMOPSO algorithm, this paper adopts an efficient constraint handling method [37], which uses the simple feasibility-based rules to update Pbest. The rules are described as follows.(1)
The feasible solution is better than the infeasible solution.(2)
For two feasible solutions, the one who has the better objective function value is better.(3)
For two infeasible solutions, the one who has the smaller constraint violations is better.According to the above criteria, objective and constraint violation information are considered separately. The update process of every iteration generation can be completed by comparing using the objective function value and the sum of all the constraints violations. In this paper, the sum of constraint violations is calculated as follows:(24)FV=∑i∈NPQ|VPQi-VPQilim⁡|ΔVPQimax⁡+∑i∈NG|QGi-QGilim⁡|ΔQGimax⁡+∑k∈NE|SLk-SLkmax⁡|ΔSLkmax⁡, where
(25)VPQilim⁡={VPQimax⁡;VPQi>VPQimax⁡VPQimin⁡;VPQi<VPQimin⁡QGjlim⁡={QGjmax⁡;QGj>Qjmax⁡QGjmin⁡;QGj<QGjmin⁡.ΔVPQimax⁡, ΔQGjmax⁡, and ΔSLkmax⁡ are the maximal constraint violations of VPQi, QGj, and SLk achieved by any individual in the current population.


### 4.3. Stopping Criteria

The most commonly used stopping criterion in MOEAs is a priori fixed number of generations. Unfortunately, it may involve a waste of computational resources, as the optimal algorithm goes on, being applied after a point where iterations get no improvement over the current solution. Thus, this paper adopts another efficient stopping criterion [38]; the standard deviation of maximum crowding distance is expressed as follows:
(26)σL=1L∑l=1L(dl-dL¯)2<δlim⁡,
where L is the length of the time window; dl is the maximal crowding distance at generation l; dL¯ is the average of dl over L generations; δlim⁡ is the threshold of the stopping criterion. In this paper, all the multiobjective algorithms will be terminated if the above criterion is satisfied or the number of iterations reaches the maximum allowable number.In addition, the follow chart of the proposed approach for solving ORPD problem is shown in Figure1.Figure 1
Computational flow of the proposed HMOPSO for ORPD problem.


## 4.1. Mixed-Variable Handling Method

In the decision variables of the ORPD problem, generator voltages are continuous variables, which can be running at any real number within the limit boundary; transformer taps and reactive compensations are discrete variables, which can only be given a value from a fixed discrete values set. Here, we use real numbers to represent continuous variables and the integer to other variables. The vector of new control variables can be written as(22)X*T=[vG1⋯vGNG,t1⋯tNT,qc1⋯qcNC],
where  vGi, ti, and  qCi are the generator bus voltage, the position of transformer tap, and the step of reactive compensation capacitor, respectively. Then, the relationship between practical parameters and new control variables can be written as
(23)VGi=vGi,Ti=1.0+ti·Δt,QCi=qCi·ΔqCi,
where Δti and ΔqCi express the ratio of transformer tap and the step capability of reactive compensation capacitor.Because the basic forms of MOPSO can only handle continuous variables, we employ a mixed-variable handling method [36]. In the initialization process, all variables are randomly generated within their upper and lower bounds at first. Then, the integer part of the variable value is picked to be its value, denoted as int⁡(xi(t)). In the iterative process, the new position for the integer variables will be selected in the neighbor integer values of the former position according to the velocity of the particle. If the velocity is more than zero, the new position is forward to plus 1. Otherwise, the position is back to minus 1.


## 4.2. Constraints Handling

From Section2, we know that the ORPD problem is subjected to several equality and inequality constraints. The equality constraints are satisfied by running Newton-Raphson method; the inequality constraints on the control variables can be automatically satisfied by setting the boundary of the search space; only the inequality constraints on the state variables need to be considered. In order to handle these constraints without introducing heavy work of turning penalty factors into the HMOPSO algorithm, this paper adopts an efficient constraint handling method [37], which uses the simple feasibility-based rules to update Pbest. The rules are described as follows.(1)
The feasible solution is better than the infeasible solution.(2)
For two feasible solutions, the one who has the better objective function value is better.(3)
For two infeasible solutions, the one who has the smaller constraint violations is better.According to the above criteria, objective and constraint violation information are considered separately. The update process of every iteration generation can be completed by comparing using the objective function value and the sum of all the constraints violations. In this paper, the sum of constraint violations is calculated as follows:(24)FV=∑i∈NPQ|VPQi-VPQilim⁡|ΔVPQimax⁡+∑i∈NG|QGi-QGilim⁡|ΔQGimax⁡+∑k∈NE|SLk-SLkmax⁡|ΔSLkmax⁡, where
(25)VPQilim⁡={VPQimax⁡;VPQi>VPQimax⁡VPQimin⁡;VPQi<VPQimin⁡QGjlim⁡={QGjmax⁡;QGj>Qjmax⁡QGjmin⁡;QGj<QGjmin⁡.ΔVPQimax⁡, ΔQGjmax⁡, and ΔSLkmax⁡ are the maximal constraint violations of VPQi, QGj, and SLk achieved by any individual in the current population.


## 4.3. Stopping Criteria

The most commonly used stopping criterion in MOEAs is a priori fixed number of generations. Unfortunately, it may involve a waste of computational resources, as the optimal algorithm goes on, being applied after a point where iterations get no improvement over the current solution. Thus, this paper adopts another efficient stopping criterion [38]; the standard deviation of maximum crowding distance is expressed as follows:
(26)σL=1L∑l=1L(dl-dL¯)2<δlim⁡,
where L is the length of the time window; dl is the maximal crowding distance at generation l; dL¯ is the average of dl over L generations; δlim⁡ is the threshold of the stopping criterion. In this paper, all the multiobjective algorithms will be terminated if the above criterion is satisfied or the number of iterations reaches the maximum allowable number.In addition, the follow chart of the proposed approach for solving ORPD problem is shown in Figure1.Figure 1
Computational flow of the proposed HMOPSO for ORPD problem.


## 5. Results and Discussion

In this section, the effectiveness and efficiency of the proposed HMOPSO algorithm for solving the ORPD problem are tested on the standard IEEE 30-bus and 118-bus power systems and the results are compared with existing popular algorithms, MOPSO, NSGA-II, MOEA/D, and other previous methods. All the techniques and simulations developed in this study are implemented on 1.83 GHz PC using MATLAB language. The load flow is run using MATPOWER 4.1 software [39] with necessary alterations in the coding.After conducting a series of experiments with different parameter settings, the optimal parameters are selected as follows. For the IEEE 30-bus system, the population sizeNpop and the maximum number of iterations Itermax⁡ are set at 50 and 100. The parameters of the stopping criteria δlim⁡ and L are set as 0.01 and 30. For the IEEE 118-bus system, Npop=100, Itermax⁡=500, δlim⁡=0.02, and L=40. The maximum size of the Pareto-optimal set is selected as 50 solutions for the two test systems. The other parameters are given below:(1)
HMOPSO: inertia weightwmax⁡=0.7, wmin⁡=0.4; acceleration constants c1,i=c2,f=2.5, c1,f=c2,i=0.5; chaotic control parameters u=4, γ0=0.48;(2)
MOPSO:w=0.729; c1=c2=2.05; the number of grids per each dimension Ngrid=50;(3)
NSGA-II: crossover probabilitypc=0.9; mutation probability pm=1/n, where n is the number of decision variables; the distribution indices for crossover and mutation operators as ηc=20 and ηm=20;(4)
MOEA/D: the neighborhood sizeTsize=30; the probability that parent solutions are selected from the neighborhood pu=0.9; the parameters in DE operator CR=1 and F=0.5; mutation probability pm=1/n; the distribution index used in the polynomial mutation ηm=20.


### 5.1. IEEE 30-Bus Power System

The standard IEEE 30-bus test system consists of 6 generators, 4 transformers, and 2 VAR compensators. Thus the number of optimized variables is 12 in this problem. Four branches (6, 9), (6, 10), (4, 12), and (27, 28) are under load tap setting transformer branches within the range [0.9,1.1] in discrete steps of 0.0125. The VAR compensators are installed at buses 10 and 24 within the interval [0,5] MVAr in discrete steps of 1 MVAr. Bus 1 is the slack bus; the buses 2, 5, 8, 11, and 13 are taken as PV-buses; the others are PQ-buses. All PV bus voltages are required to be maintained within the range of 0.95–1.1 p.u. and for PQ bus 0.9–1.05 p.u. The detailed data are given in [40, 41].To demonstrate the effectiveness of the proposed algorithm, three different cases have been considered as follows:Case 1: minimization of power loss and voltage deviation.Case 2: minimization of power loss and voltage stability index.Case 3: minimization of power loss, voltage deviation, and voltage stability index.Owing to the randomness of the proposed HMOPSO, MOPSO, NSGA-II, and MOEA/D, each case has 30 independent runs for different algorithms.


#### 5.1.1. Minimization of Power Loss and Voltage Deviation

In this case, two objectives are considered, that is, power loss and voltage deviation. The best Pareto fronts obtained out of 30 runs using HMOPSO, MOPSO, NSGA-II, and MOEA/D are shown in Figure2. It can be seen that the HMOPSO can obtain better convergence and spread solution over the decision space. Table 1 summarizes the best extreme solutions obtained by using different algorithms for PL and VD minimization in IEEE 30-bus system. It is clear that the proposed HMOPSO can find better extreme solutions than other methods.Table 1
Optimal settings of control variables forPL and VD minimization in IEEE 30-bus system (Case  1).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestPL
Best VD
BestPL
Best VD
BestPL
Best VD
V
G
1
1.05
1.0708
1.0183
1.0658
1.0086
1.0800
1.0112
1.0715
1.0120
V
G
2
1.04
1.0630
1.0084
1.0561
1.0008
1.0703
1.0041
1.0642
1.0022
V
G
5
1.01
1.0415
1.0030
1.0346
1.0000
1.0482
1.0025
1.0401
0.9998
V
G
8
1.01
1.0405
0.9995
1.0351
1.0013
1.0494
1.0060
1.0409
1.0000
V
G
13
1.05
1.0744
1.0300
1.0417
0.9933
1.1000
1.0261
1.0773
1.0982
V
G
1
1.05
1.0620
1.0418
1.0694
1.0650
1.0995
1.0500
1.0598
1.0376
T
6
–
9
1.078
1.075
1
1.025
1
1.0125
0.9875
1.0375
1.05
T
6
–
10
1.069
0.9
0.9125
0.9
0.9125
0.9125
0.9875
0.95
0.9875
T
4
–
12
1.032
1
1
0.9875
1.025
1.0125
1
0.975
0.9875
T
27
-
28
1.068
0.9625
0.925
0.95
0.925
0.95
0.925
0.9625
0.925
Q
10
0.0
0.05
0.01
0.01
0.04
0.1
0.02
0.02
0.01
Q
24
0.0
0
0.05
0.02
0
0.1
0.03
0.01
0.05
P
L (MW)
5.8327
4.9849
5.5968
5.0465
5.6749
4.9595
5.6553
4.8536
5.8650
VD (pu)
1.1582
0.7196
0.1222
0.7248
0.1274
0.7452
0.1237
0.8842
0.0913Figure 2
Best Pareto fronts obtained using HMOPSO, MOPSO, NSGA-II, and MOEA/D forPL and VD minimization of IEEE 30-bus system.For a fair and complete comparison, the extreme solutions obtained by the proposed approach are compared with SPEA [21], enhanced general passive congregation (GPAC) [42], coordinated aggregation (CA) [42], BBO [12], GSA [27], OSMGSA [27], DE [23], and modified DE (MDE) [23], in which the constraints, control variables, and initial settings of the test systems are the same as our works. By comparing the obtained results in Table 2, the proposed algorithm can converge to lower values of real power loss and voltage deviations with respect to other reported algorithms.Table 2
Comparison of best solutions of different algorithms for IEEE 30-bus system (Case  1).
Objective
Initial
SPEA [21]
GPAC [42]
CA [42]
BBO [12]
GSA [27]
OSAMGSA [27]
DE [23]
MDE [23]
HMOPSO
P
L (MW)
5.8327
5.1167
5.0923
5.0921
4.9650
5.0924
5.0713
4.8750
4.8736
4.8536
VD (pu)
1.1582
0.1442
0.1274
0.1225
0.15499
0.1133
0.1126
0.0927
0.0910
0.0913


#### 5.1.2. Minimization of Power Loss and Voltage Stability Index

To validate the effectiveness of the proposed HMOPSO algorithm for enhancing the voltage security under stress condition and contingency states, two different conditions are considered in this case. Case  2(a) is heavy load condition. The load on each bus is uniformly increased to 150% of the base load condition to analyze the voltage stability level of the system under severe conditions; Case  2(b) is contingency state and the single line outrage (27–30) is simulated. Tables3 and 4 give the optimal control variables obtained by HMOPSO, MOSPO, NSGA-II, and MOEA/D for the above two cases, respectively. From the results, it is clear that the value of Lmax⁡ decreases and voltage stability is improved after the application of optimization and the results of the proposed HMOPSO are better than other two algorithms.Table 3
Optimal settings of control variables forPL and Lmax⁡ minimization in IEEE 30-bus system (Case  2(a): heavy load).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
V
G
1
1.05
1.1000
1.0965
1.0978
1.0947
1.1000
1.0959
1.0923
1.0714
V
G
2
1.04
1.0807
1.0844
1.0865
1.0882
1.0871
1.0771
1.0871
1.0545
V
G
5
1.01
1.0496
1.0393
1.0547
1.0522
1.0400
1.0241
1.0400
1.0142
V
G
8
1.01
1.0522
1.0547
1.0346
1.0313
1.0502
1.0525
1.0502
1.0103
V
G
13
1.05
1.0962
1.1000
1.0864
1.0897
1.0831
1.0626
1.0831
1.0881
V
G
1
1.05
1.0918
1.0784
1.0323
1.0344
1.0560
1.0815
1.0560
1.0678
T
6
–
9
1.078
1.0125
1
0.9875
0.9875
1.0125
1
1.0125
0.9875
T
6
–
10
1.069
0.95
0.9625
0.9
0.9
0.925
0.9
0.925
0.925
T
4
–
12
1.032
1.025
0.9375
0.9375
0.9
0.9875
1
0.9875
0.9375
T
27
-
28
1.068
0.95
0.9
0.925
0.925
0.925
0.9
0.925
0.925
Q
10
0.0
0.03
0.03
0
0.01
0.02
0.02
0
0.02
Q
24
0.0
0.03
0.03
0.04
0.03
0.01
0.03
0.01
0.01
P
L (MW)
13.9971
12.1380
12.5101
12.8269
13.0480
12.1136
12.7399
12.0156
12.9451
L
max
⁡
0.2410
0.1837
0.1715
0.1856
0.1824
0.1858
0.1739
0.1883
0.1636Table 4
Optimal settings of control variables forPL and Lmax⁡ minimization in IEEE 30-bus system (Case  2(b): contingency state).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
V
G
1
1.05
1.0844
1.0837
1.0848
1.0803
1.0888
1.0754
1.0772
1.0512
V
G
2
1.04
1.0735
1.0729
1.0723
1.0689
1.0734
1.0668
1.0733
1.0451
V
G
5
1.01
1.0537
1.0533
1.0525
1.0456
1.0504
1.0297
1.0624
1.0175
V
G
8
1.01
1.0558
1.0566
1.0525
1.0475
1.0506
1.0413
1.0555
1.0136
V
G
13
1.05
1.0963
1.0962
1.0688
1.0674
1.0762
1.0771
1.0862
1.0747
V
G
1
1.05
1.0793
1.0788
1.0805
1.0758
1.0649
1.0656
1.0646
1.0531
T
6
–
9
1.078
1.0125
1.0125
0.9875
0.9875
0.9875
1.0375
1.0125
0.9875
T
6
–
10
1.069
0.9
0.9
0.9125
0.925
0.9
0.9375
0.9875
0.9125
T
4
–
12
1.032
0.9875
0.9875
1
1
1
1.05
1.025
1.05
T
27
-
28
1.068
0.95
0.95
0.95
0.9375
0.95
0.925
0.95
0.925
Q
10
0.0
0.01
0
0.01
0.02
0.02
0.01
0.02
0.02
Q
24
0.0
0.04
0.3
0.02
0.02
0.01
0.05
0.01
0.03
P
L (MW)
6.5179
5.4436
5.5493
5.5094
5.5548
5.4504
5.6932
5.3076
5.9774
L
max
⁡
0.2956
0.2162
0.2061
0.2177
0.2168
0.2197
0.1980
0.2386
0. 1842


#### 5.1.3. Minimization of Power Loss, Voltage Deviation, and Voltage Stability Index

In this case, three objectives, that is, power loss, voltage deviation, and voltage stability index are optimized simultaneously. Table5 gives the best optimal settings obtained results using different algorithms. It can be observed from the results that the proposed HMOPSO outperforms other algorithms in terms of the values of the objective functions.Table 5
Optimal solutions of control variables forPL, VD, and Lmax⁡ minimization in IEEE 30-bus system (Case  3).
Variable
Initial
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
V
G
1
1.05
1.0688
1.0256
1.0594
1.0828
1.0176
1.0706
1.0886
1.0425
1.0748
V
G
2
1.04
1.0570
1.0246
1.0053
1.0727
1.0044
1.0642
1.0813
1.0316
1.0655
V
G
5
1.01
1.0338
1.0248
1.0780
1.0478
0.9998
1.0349
1.0547
1.0143
1.0452
V
G
8
1.01
1.0340
1.0057
1.0518
1.0480
0.9999
1.0362
1.0591
1.0132
1.0461
V
G
13
1.05
1.0877
1.0543
1.0831
1.0938
1.0938
1.0627
1.0986
1.0664
1.0456
V
G
1
1.05
1.0538
1.0287
1.0503
1.0951
1.0359
1.0597
1.0714
1.0238
1.0542
T
6
-
9
1.078
0.9875
1
1.025
1.025
1.025
0.9875
1.025
1.05
1
T
6
–
10
1.069
0.9375
0.9
0.9625
0.9
1
0.9625
0.9375
0.9125
0.9
T
4
–
12
1.032
0.9625
0.95
0.9375
1.0125
0.9875
1.0125
1.0125
0.975
0.9375
T
27
-
28
1.068
0.95
0.925
0.95
0.975
0.925
0.9
0.9625
0.925
0.925
Q
10
0.0
0.03
0.02
0
0
0
0.02
0.02
0.01
0.04
Q
24
0.0
0.01
0.03
0.02
0.08
0.06
0.01
0.03
0.05
0.04
P
L (MW)
5.8327
5.2064
5.8076
5.8909
4.9803
5.6340
5.3730
4.8638
5.8846
5.4867
VD (pu)
1.1582
0.6114
0.1268
0.6963
0.7145
0.1137
0.876
0.9625
0.1098
0.9886
L
max
⁡
0.1728
0.1394
0.1437
0.1315
0.1385
0.1472
0.1301
0.1367
0.1567
0.1203


### 5.2. IEEE 118-Bus Power System

To validate the possibility of the proposed HMOPSO to the large-scale power system, the IEEE 118-bus system is adopted. This system consists of 54 generators, 9 transformers, and 12 capacitor banks. Thus, the dimension of control variables in this case is 75. For more information about the system, one can refer to [43]. The step size and interval of the decision variables are same as in IEEE 30-bus system.The best extremesolutions obtained by different algorithms for IEEE 118-bus system are given in Table 6. Due to the space limitation, the final optimal settings of 75 decision variables are not given in Table 6. For more validation, the obtained results are compared with GPAC [42], CA [42], and PSO with differentially perturbed velocity hybrid algorithm with adaptive acceleration coefficient (APSODV) [44], BBO [12], GSA [11], and OGS [45], which are summarized in Table 7. The comparison results show that the proposed HMOPSO can yield better solutions than other reported approaches.Table 6
Optimal solutions forPL, VD, and Lmax⁡ minimization in IEEE 118-bus system.
Object
Initial
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
P
L (MW)
132.8629
124.1573
131.781
128.5268
121.458
134.647
131.332
120.0564
138.8756
132.4865
VD (pu)
1.4393
1.3478
0.5841
0.9013
1.255
0.5231
0.968
1.5536
0.4237
1.2524
L
max
⁡
0.0694
0.0569
0.0571
0.0552
0.0571
0.0586
0.0534
0.0582
0.0598
0.0512Table 7
Comparison of the best solutions of different algorithms for IEEE 118-bus system.
Object
Initial
GPAC [42]
CA [42]
APSODV [44]
BBO [12]
GSA [11]
OGSA [45]
HMOPSO
P
L (MW)
132.8629
131.9083
131.8638
—
128.97
122.6762
126.99
120.0564
VD (pu)
1.4393
1.2849
1.2756
0.4666
0.5020
0.5091
0.3666
0.4237
L
max
⁡
0.0694
—
—
0.0538
0.0523
0.0600
0.0512


### 5.3. Multiobjective Performance Metrics Analysis

In multiobjective optimization processes, two goals are normally taken into account:(1) convergence to the true Pareto-optimal set; (2) maintenance of diversity in solutions of the Pareto-optimal set. In order to obtain an accurate quantitative comparison of the performance of the proposed approach, this paper implements three different performance metrics: generational distance (GD) [30], minimal spacing (MSP) [31], and hypervolume (HV) [32]. A brief introduction of these metrics is given as follows.Generational Distance (GD). GD is used to evaluate the closeness of the nondominated set obtained by an algorithm to the reference Pareto-optimal front. Consider
(27)GD=∑i=1ndi2n,
where  n is the number of vectors in the reference Pareto-optimal front and di is the Euclidean distance (measured in objective space) between each of these and the nearest member of the Pareto optimal set.Minimal Spacing (MSP). MSP is used to evaluate how evenly the nondominated solutions are distributed in the objective. Consider
(28)MSP=1n-1∑i=1n(di-d¯)2,
where di=mink∑m=1Nobj|fmi-fmk|, Nobj is the number of objectives, i is the index of the solution which is marked as the current seed, and k is the index of the solution which is still unmarked at the current iteration.Hypervolume (HV). HV provides the combined qualitative information about closeness and diversity in obtained Pareto-optimal fronts. It calculates the volume covered by the members of the approximate Pareto-optimal Q for problems where all objectives are to be minimized. Mathematically, for each solution i∈Q, a hypercube vi is constructed with a reference point W and the solution i is the diagonal corners of the hypercube. The reference point can simply be found by constructing a vector of worst objective function values. Thereafter, a union of all hypervolume (HV) is calculated as follows:
(29)HV=volume(⋃i=1|Q|vi).In order to evaluate the performance metrics of these algorithms, each algorithm is executed 100 times for two test systems. The statistic comparison results of performance metrics for two test systems are shown in Table8. It is clear that the proposed algorithm can obtain better Pareto front with respect to other algorithms, because the average and standard deviation performances of GD, MSP, and HV in HMOPSO are the best.Table 8
Statistical results of performance metric for IEEE 30-bus and IEEE 118-bus systems.
Performance measures
IEEE 30-bus system
IEEE 118-bus system
NSGA-II
MOAE/D
MOPSO
HMOPSO
NSGA-II
MOAE/D
MOPSO
HMOPSO
GD
Mean
0.113627
0.015378
0.015862
0.012735
0.330574
0.022916
0.024811
0.018853
Std. dev.
0.037141
0.00414
0.003135
0.002283
0.064578
0.003867
0.005095
0.003456
MSP
Mean
0.263546
0.32752
0.289632
0.186752
0.287658
0.41543
0.325624
0.253628
Std. dev.
0.040234
0.01564
0.02466
0.011338
0.058462
0.003133
0.038236
0.031247
HV
Mean
0.185746
0.22676
0.211651
0.24836
0.132385
0.153652
0.165727
0.201435
Std. dev.
0.035541
0.031572
0.041701
0.025813
0.042356
0.055814
0.063217
0.040472Table9 shows the statistic results of computational time over 100 independent runs on the two test power systems. It is quite evident that the run time of HMOPSO is computationally more efficient than other techniques.Table 9
Computational time of different algorithm for IEEE 30-bus and 118-bus systems.
Time
IEEE 30-bus system
IEEE 118-bus system
NSGA-II
MOAE/D
MOPSO
HMOPSO
NSGA-II
MOAE/D
MOPSO
HMOPSO
Best (s)
168.42
54.34
40.16
23.52
405.78
167.59
142.65
110.17
Mean (s)
195.76
65.29
48.35
35.47
451.69
191.58
174.34
138.88
Worst (s)
213.84
78.67
55.37
41.78
510.36
215.28
191.28
165.54
Std. dev.
7.48
4.13
3.18
3.46
23.57
11.38
16.55
13.76


## 5.1. IEEE 30-Bus Power System

The standard IEEE 30-bus test system consists of 6 generators, 4 transformers, and 2 VAR compensators. Thus the number of optimized variables is 12 in this problem. Four branches (6, 9), (6, 10), (4, 12), and (27, 28) are under load tap setting transformer branches within the range [0.9,1.1] in discrete steps of 0.0125. The VAR compensators are installed at buses 10 and 24 within the interval [0,5] MVAr in discrete steps of 1 MVAr. Bus 1 is the slack bus; the buses 2, 5, 8, 11, and 13 are taken as PV-buses; the others are PQ-buses. All PV bus voltages are required to be maintained within the range of 0.95–1.1 p.u. and for PQ bus 0.9–1.05 p.u. The detailed data are given in [40, 41].To demonstrate the effectiveness of the proposed algorithm, three different cases have been considered as follows:Case 1: minimization of power loss and voltage deviation.Case 2: minimization of power loss and voltage stability index.Case 3: minimization of power loss, voltage deviation, and voltage stability index.Owing to the randomness of the proposed HMOPSO, MOPSO, NSGA-II, and MOEA/D, each case has 30 independent runs for different algorithms.


### 5.1.1. Minimization of Power Loss and Voltage Deviation

In this case, two objectives are considered, that is, power loss and voltage deviation. The best Pareto fronts obtained out of 30 runs using HMOPSO, MOPSO, NSGA-II, and MOEA/D are shown in Figure2. It can be seen that the HMOPSO can obtain better convergence and spread solution over the decision space. Table 1 summarizes the best extreme solutions obtained by using different algorithms for PL and VD minimization in IEEE 30-bus system. It is clear that the proposed HMOPSO can find better extreme solutions than other methods.Table 1
Optimal settings of control variables forPL and VD minimization in IEEE 30-bus system (Case  1).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestPL
Best VD
BestPL
Best VD
BestPL
Best VD
V
G
1
1.05
1.0708
1.0183
1.0658
1.0086
1.0800
1.0112
1.0715
1.0120
V
G
2
1.04
1.0630
1.0084
1.0561
1.0008
1.0703
1.0041
1.0642
1.0022
V
G
5
1.01
1.0415
1.0030
1.0346
1.0000
1.0482
1.0025
1.0401
0.9998
V
G
8
1.01
1.0405
0.9995
1.0351
1.0013
1.0494
1.0060
1.0409
1.0000
V
G
13
1.05
1.0744
1.0300
1.0417
0.9933
1.1000
1.0261
1.0773
1.0982
V
G
1
1.05
1.0620
1.0418
1.0694
1.0650
1.0995
1.0500
1.0598
1.0376
T
6
–
9
1.078
1.075
1
1.025
1
1.0125
0.9875
1.0375
1.05
T
6
–
10
1.069
0.9
0.9125
0.9
0.9125
0.9125
0.9875
0.95
0.9875
T
4
–
12
1.032
1
1
0.9875
1.025
1.0125
1
0.975
0.9875
T
27
-
28
1.068
0.9625
0.925
0.95
0.925
0.95
0.925
0.9625
0.925
Q
10
0.0
0.05
0.01
0.01
0.04
0.1
0.02
0.02
0.01
Q
24
0.0
0
0.05
0.02
0
0.1
0.03
0.01
0.05
P
L (MW)
5.8327
4.9849
5.5968
5.0465
5.6749
4.9595
5.6553
4.8536
5.8650
VD (pu)
1.1582
0.7196
0.1222
0.7248
0.1274
0.7452
0.1237
0.8842
0.0913Figure 2
Best Pareto fronts obtained using HMOPSO, MOPSO, NSGA-II, and MOEA/D forPL and VD minimization of IEEE 30-bus system.For a fair and complete comparison, the extreme solutions obtained by the proposed approach are compared with SPEA [21], enhanced general passive congregation (GPAC) [42], coordinated aggregation (CA) [42], BBO [12], GSA [27], OSMGSA [27], DE [23], and modified DE (MDE) [23], in which the constraints, control variables, and initial settings of the test systems are the same as our works. By comparing the obtained results in Table 2, the proposed algorithm can converge to lower values of real power loss and voltage deviations with respect to other reported algorithms.Table 2
Comparison of best solutions of different algorithms for IEEE 30-bus system (Case  1).
Objective
Initial
SPEA [21]
GPAC [42]
CA [42]
BBO [12]
GSA [27]
OSAMGSA [27]
DE [23]
MDE [23]
HMOPSO
P
L (MW)
5.8327
5.1167
5.0923
5.0921
4.9650
5.0924
5.0713
4.8750
4.8736
4.8536
VD (pu)
1.1582
0.1442
0.1274
0.1225
0.15499
0.1133
0.1126
0.0927
0.0910
0.0913


### 5.1.2. Minimization of Power Loss and Voltage Stability Index

To validate the effectiveness of the proposed HMOPSO algorithm for enhancing the voltage security under stress condition and contingency states, two different conditions are considered in this case. Case  2(a) is heavy load condition. The load on each bus is uniformly increased to 150% of the base load condition to analyze the voltage stability level of the system under severe conditions; Case  2(b) is contingency state and the single line outrage (27–30) is simulated. Tables3 and 4 give the optimal control variables obtained by HMOPSO, MOSPO, NSGA-II, and MOEA/D for the above two cases, respectively. From the results, it is clear that the value of Lmax⁡ decreases and voltage stability is improved after the application of optimization and the results of the proposed HMOPSO are better than other two algorithms.Table 3
Optimal settings of control variables forPL and Lmax⁡ minimization in IEEE 30-bus system (Case  2(a): heavy load).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
V
G
1
1.05
1.1000
1.0965
1.0978
1.0947
1.1000
1.0959
1.0923
1.0714
V
G
2
1.04
1.0807
1.0844
1.0865
1.0882
1.0871
1.0771
1.0871
1.0545
V
G
5
1.01
1.0496
1.0393
1.0547
1.0522
1.0400
1.0241
1.0400
1.0142
V
G
8
1.01
1.0522
1.0547
1.0346
1.0313
1.0502
1.0525
1.0502
1.0103
V
G
13
1.05
1.0962
1.1000
1.0864
1.0897
1.0831
1.0626
1.0831
1.0881
V
G
1
1.05
1.0918
1.0784
1.0323
1.0344
1.0560
1.0815
1.0560
1.0678
T
6
–
9
1.078
1.0125
1
0.9875
0.9875
1.0125
1
1.0125
0.9875
T
6
–
10
1.069
0.95
0.9625
0.9
0.9
0.925
0.9
0.925
0.925
T
4
–
12
1.032
1.025
0.9375
0.9375
0.9
0.9875
1
0.9875
0.9375
T
27
-
28
1.068
0.95
0.9
0.925
0.925
0.925
0.9
0.925
0.925
Q
10
0.0
0.03
0.03
0
0.01
0.02
0.02
0
0.02
Q
24
0.0
0.03
0.03
0.04
0.03
0.01
0.03
0.01
0.01
P
L (MW)
13.9971
12.1380
12.5101
12.8269
13.0480
12.1136
12.7399
12.0156
12.9451
L
max
⁡
0.2410
0.1837
0.1715
0.1856
0.1824
0.1858
0.1739
0.1883
0.1636Table 4
Optimal settings of control variables forPL and Lmax⁡ minimization in IEEE 30-bus system (Case  2(b): contingency state).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
V
G
1
1.05
1.0844
1.0837
1.0848
1.0803
1.0888
1.0754
1.0772
1.0512
V
G
2
1.04
1.0735
1.0729
1.0723
1.0689
1.0734
1.0668
1.0733
1.0451
V
G
5
1.01
1.0537
1.0533
1.0525
1.0456
1.0504
1.0297
1.0624
1.0175
V
G
8
1.01
1.0558
1.0566
1.0525
1.0475
1.0506
1.0413
1.0555
1.0136
V
G
13
1.05
1.0963
1.0962
1.0688
1.0674
1.0762
1.0771
1.0862
1.0747
V
G
1
1.05
1.0793
1.0788
1.0805
1.0758
1.0649
1.0656
1.0646
1.0531
T
6
–
9
1.078
1.0125
1.0125
0.9875
0.9875
0.9875
1.0375
1.0125
0.9875
T
6
–
10
1.069
0.9
0.9
0.9125
0.925
0.9
0.9375
0.9875
0.9125
T
4
–
12
1.032
0.9875
0.9875
1
1
1
1.05
1.025
1.05
T
27
-
28
1.068
0.95
0.95
0.95
0.9375
0.95
0.925
0.95
0.925
Q
10
0.0
0.01
0
0.01
0.02
0.02
0.01
0.02
0.02
Q
24
0.0
0.04
0.3
0.02
0.02
0.01
0.05
0.01
0.03
P
L (MW)
6.5179
5.4436
5.5493
5.5094
5.5548
5.4504
5.6932
5.3076
5.9774
L
max
⁡
0.2956
0.2162
0.2061
0.2177
0.2168
0.2197
0.1980
0.2386
0. 1842


### 5.1.3. Minimization of Power Loss, Voltage Deviation, and Voltage Stability Index

In this case, three objectives, that is, power loss, voltage deviation, and voltage stability index are optimized simultaneously. Table5 gives the best optimal settings obtained results using different algorithms. It can be observed from the results that the proposed HMOPSO outperforms other algorithms in terms of the values of the objective functions.Table 5
Optimal solutions of control variables forPL, VD, and Lmax⁡ minimization in IEEE 30-bus system (Case  3).
Variable
Initial
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
V
G
1
1.05
1.0688
1.0256
1.0594
1.0828
1.0176
1.0706
1.0886
1.0425
1.0748
V
G
2
1.04
1.0570
1.0246
1.0053
1.0727
1.0044
1.0642
1.0813
1.0316
1.0655
V
G
5
1.01
1.0338
1.0248
1.0780
1.0478
0.9998
1.0349
1.0547
1.0143
1.0452
V
G
8
1.01
1.0340
1.0057
1.0518
1.0480
0.9999
1.0362
1.0591
1.0132
1.0461
V
G
13
1.05
1.0877
1.0543
1.0831
1.0938
1.0938
1.0627
1.0986
1.0664
1.0456
V
G
1
1.05
1.0538
1.0287
1.0503
1.0951
1.0359
1.0597
1.0714
1.0238
1.0542
T
6
-
9
1.078
0.9875
1
1.025
1.025
1.025
0.9875
1.025
1.05
1
T
6
–
10
1.069
0.9375
0.9
0.9625
0.9
1
0.9625
0.9375
0.9125
0.9
T
4
–
12
1.032
0.9625
0.95
0.9375
1.0125
0.9875
1.0125
1.0125
0.975
0.9375
T
27
-
28
1.068
0.95
0.925
0.95
0.975
0.925
0.9
0.9625
0.925
0.925
Q
10
0.0
0.03
0.02
0
0
0
0.02
0.02
0.01
0.04
Q
24
0.0
0.01
0.03
0.02
0.08
0.06
0.01
0.03
0.05
0.04
P
L (MW)
5.8327
5.2064
5.8076
5.8909
4.9803
5.6340
5.3730
4.8638
5.8846
5.4867
VD (pu)
1.1582
0.6114
0.1268
0.6963
0.7145
0.1137
0.876
0.9625
0.1098
0.9886
L
max
⁡
0.1728
0.1394
0.1437
0.1315
0.1385
0.1472
0.1301
0.1367
0.1567
0.1203


## 5.1.1. Minimization of Power Loss and Voltage Deviation

In this case, two objectives are considered, that is, power loss and voltage deviation. The best Pareto fronts obtained out of 30 runs using HMOPSO, MOPSO, NSGA-II, and MOEA/D are shown in Figure2. It can be seen that the HMOPSO can obtain better convergence and spread solution over the decision space. Table 1 summarizes the best extreme solutions obtained by using different algorithms for PL and VD minimization in IEEE 30-bus system. It is clear that the proposed HMOPSO can find better extreme solutions than other methods.Table 1
Optimal settings of control variables forPL and VD minimization in IEEE 30-bus system (Case  1).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestPL
Best VD
BestPL
Best VD
BestPL
Best VD
V
G
1
1.05
1.0708
1.0183
1.0658
1.0086
1.0800
1.0112
1.0715
1.0120
V
G
2
1.04
1.0630
1.0084
1.0561
1.0008
1.0703
1.0041
1.0642
1.0022
V
G
5
1.01
1.0415
1.0030
1.0346
1.0000
1.0482
1.0025
1.0401
0.9998
V
G
8
1.01
1.0405
0.9995
1.0351
1.0013
1.0494
1.0060
1.0409
1.0000
V
G
13
1.05
1.0744
1.0300
1.0417
0.9933
1.1000
1.0261
1.0773
1.0982
V
G
1
1.05
1.0620
1.0418
1.0694
1.0650
1.0995
1.0500
1.0598
1.0376
T
6
–
9
1.078
1.075
1
1.025
1
1.0125
0.9875
1.0375
1.05
T
6
–
10
1.069
0.9
0.9125
0.9
0.9125
0.9125
0.9875
0.95
0.9875
T
4
–
12
1.032
1
1
0.9875
1.025
1.0125
1
0.975
0.9875
T
27
-
28
1.068
0.9625
0.925
0.95
0.925
0.95
0.925
0.9625
0.925
Q
10
0.0
0.05
0.01
0.01
0.04
0.1
0.02
0.02
0.01
Q
24
0.0
0
0.05
0.02
0
0.1
0.03
0.01
0.05
P
L (MW)
5.8327
4.9849
5.5968
5.0465
5.6749
4.9595
5.6553
4.8536
5.8650
VD (pu)
1.1582
0.7196
0.1222
0.7248
0.1274
0.7452
0.1237
0.8842
0.0913Figure 2
Best Pareto fronts obtained using HMOPSO, MOPSO, NSGA-II, and MOEA/D forPL and VD minimization of IEEE 30-bus system.For a fair and complete comparison, the extreme solutions obtained by the proposed approach are compared with SPEA [21], enhanced general passive congregation (GPAC) [42], coordinated aggregation (CA) [42], BBO [12], GSA [27], OSMGSA [27], DE [23], and modified DE (MDE) [23], in which the constraints, control variables, and initial settings of the test systems are the same as our works. By comparing the obtained results in Table 2, the proposed algorithm can converge to lower values of real power loss and voltage deviations with respect to other reported algorithms.Table 2
Comparison of best solutions of different algorithms for IEEE 30-bus system (Case  1).
Objective
Initial
SPEA [21]
GPAC [42]
CA [42]
BBO [12]
GSA [27]
OSAMGSA [27]
DE [23]
MDE [23]
HMOPSO
P
L (MW)
5.8327
5.1167
5.0923
5.0921
4.9650
5.0924
5.0713
4.8750
4.8736
4.8536
VD (pu)
1.1582
0.1442
0.1274
0.1225
0.15499
0.1133
0.1126
0.0927
0.0910
0.0913


## 5.1.2. Minimization of Power Loss and Voltage Stability Index

To validate the effectiveness of the proposed HMOPSO algorithm for enhancing the voltage security under stress condition and contingency states, two different conditions are considered in this case. Case  2(a) is heavy load condition. The load on each bus is uniformly increased to 150% of the base load condition to analyze the voltage stability level of the system under severe conditions; Case  2(b) is contingency state and the single line outrage (27–30) is simulated. Tables3 and 4 give the optimal control variables obtained by HMOPSO, MOSPO, NSGA-II, and MOEA/D for the above two cases, respectively. From the results, it is clear that the value of Lmax⁡ decreases and voltage stability is improved after the application of optimization and the results of the proposed HMOPSO are better than other two algorithms.Table 3
Optimal settings of control variables forPL and Lmax⁡ minimization in IEEE 30-bus system (Case  2(a): heavy load).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
V
G
1
1.05
1.1000
1.0965
1.0978
1.0947
1.1000
1.0959
1.0923
1.0714
V
G
2
1.04
1.0807
1.0844
1.0865
1.0882
1.0871
1.0771
1.0871
1.0545
V
G
5
1.01
1.0496
1.0393
1.0547
1.0522
1.0400
1.0241
1.0400
1.0142
V
G
8
1.01
1.0522
1.0547
1.0346
1.0313
1.0502
1.0525
1.0502
1.0103
V
G
13
1.05
1.0962
1.1000
1.0864
1.0897
1.0831
1.0626
1.0831
1.0881
V
G
1
1.05
1.0918
1.0784
1.0323
1.0344
1.0560
1.0815
1.0560
1.0678
T
6
–
9
1.078
1.0125
1
0.9875
0.9875
1.0125
1
1.0125
0.9875
T
6
–
10
1.069
0.95
0.9625
0.9
0.9
0.925
0.9
0.925
0.925
T
4
–
12
1.032
1.025
0.9375
0.9375
0.9
0.9875
1
0.9875
0.9375
T
27
-
28
1.068
0.95
0.9
0.925
0.925
0.925
0.9
0.925
0.925
Q
10
0.0
0.03
0.03
0
0.01
0.02
0.02
0
0.02
Q
24
0.0
0.03
0.03
0.04
0.03
0.01
0.03
0.01
0.01
P
L (MW)
13.9971
12.1380
12.5101
12.8269
13.0480
12.1136
12.7399
12.0156
12.9451
L
max
⁡
0.2410
0.1837
0.1715
0.1856
0.1824
0.1858
0.1739
0.1883
0.1636Table 4
Optimal settings of control variables forPL and Lmax⁡ minimization in IEEE 30-bus system (Case  2(b): contingency state).
Variable
Initial
MOPSO
NSGA-II
MOEA/D
HMOPSO
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
BestPL
BestLmax⁡
V
G
1
1.05
1.0844
1.0837
1.0848
1.0803
1.0888
1.0754
1.0772
1.0512
V
G
2
1.04
1.0735
1.0729
1.0723
1.0689
1.0734
1.0668
1.0733
1.0451
V
G
5
1.01
1.0537
1.0533
1.0525
1.0456
1.0504
1.0297
1.0624
1.0175
V
G
8
1.01
1.0558
1.0566
1.0525
1.0475
1.0506
1.0413
1.0555
1.0136
V
G
13
1.05
1.0963
1.0962
1.0688
1.0674
1.0762
1.0771
1.0862
1.0747
V
G
1
1.05
1.0793
1.0788
1.0805
1.0758
1.0649
1.0656
1.0646
1.0531
T
6
–
9
1.078
1.0125
1.0125
0.9875
0.9875
0.9875
1.0375
1.0125
0.9875
T
6
–
10
1.069
0.9
0.9
0.9125
0.925
0.9
0.9375
0.9875
0.9125
T
4
–
12
1.032
0.9875
0.9875
1
1
1
1.05
1.025
1.05
T
27
-
28
1.068
0.95
0.95
0.95
0.9375
0.95
0.925
0.95
0.925
Q
10
0.0
0.01
0
0.01
0.02
0.02
0.01
0.02
0.02
Q
24
0.0
0.04
0.3
0.02
0.02
0.01
0.05
0.01
0.03
P
L (MW)
6.5179
5.4436
5.5493
5.5094
5.5548
5.4504
5.6932
5.3076
5.9774
L
max
⁡
0.2956
0.2162
0.2061
0.2177
0.2168
0.2197
0.1980
0.2386
0. 1842


## 5.1.3. Minimization of Power Loss, Voltage Deviation, and Voltage Stability Index

In this case, three objectives, that is, power loss, voltage deviation, and voltage stability index are optimized simultaneously. Table5 gives the best optimal settings obtained results using different algorithms. It can be observed from the results that the proposed HMOPSO outperforms other algorithms in terms of the values of the objective functions.Table 5
Optimal solutions of control variables forPL, VD, and Lmax⁡ minimization in IEEE 30-bus system (Case  3).
Variable
Initial
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
V
G
1
1.05
1.0688
1.0256
1.0594
1.0828
1.0176
1.0706
1.0886
1.0425
1.0748
V
G
2
1.04
1.0570
1.0246
1.0053
1.0727
1.0044
1.0642
1.0813
1.0316
1.0655
V
G
5
1.01
1.0338
1.0248
1.0780
1.0478
0.9998
1.0349
1.0547
1.0143
1.0452
V
G
8
1.01
1.0340
1.0057
1.0518
1.0480
0.9999
1.0362
1.0591
1.0132
1.0461
V
G
13
1.05
1.0877
1.0543
1.0831
1.0938
1.0938
1.0627
1.0986
1.0664
1.0456
V
G
1
1.05
1.0538
1.0287
1.0503
1.0951
1.0359
1.0597
1.0714
1.0238
1.0542
T
6
-
9
1.078
0.9875
1
1.025
1.025
1.025
0.9875
1.025
1.05
1
T
6
–
10
1.069
0.9375
0.9
0.9625
0.9
1
0.9625
0.9375
0.9125
0.9
T
4
–
12
1.032
0.9625
0.95
0.9375
1.0125
0.9875
1.0125
1.0125
0.975
0.9375
T
27
-
28
1.068
0.95
0.925
0.95
0.975
0.925
0.9
0.9625
0.925
0.925
Q
10
0.0
0.03
0.02
0
0
0
0.02
0.02
0.01
0.04
Q
24
0.0
0.01
0.03
0.02
0.08
0.06
0.01
0.03
0.05
0.04
P
L (MW)
5.8327
5.2064
5.8076
5.8909
4.9803
5.6340
5.3730
4.8638
5.8846
5.4867
VD (pu)
1.1582
0.6114
0.1268
0.6963
0.7145
0.1137
0.876
0.9625
0.1098
0.9886
L
max
⁡
0.1728
0.1394
0.1437
0.1315
0.1385
0.1472
0.1301
0.1367
0.1567
0.1203


## 5.2. IEEE 118-Bus Power System

To validate the possibility of the proposed HMOPSO to the large-scale power system, the IEEE 118-bus system is adopted. This system consists of 54 generators, 9 transformers, and 12 capacitor banks. Thus, the dimension of control variables in this case is 75. For more information about the system, one can refer to [43]. The step size and interval of the decision variables are same as in IEEE 30-bus system.The best extremesolutions obtained by different algorithms for IEEE 118-bus system are given in Table 6. Due to the space limitation, the final optimal settings of 75 decision variables are not given in Table 6. For more validation, the obtained results are compared with GPAC [42], CA [42], and PSO with differentially perturbed velocity hybrid algorithm with adaptive acceleration coefficient (APSODV) [44], BBO [12], GSA [11], and OGS [45], which are summarized in Table 7. The comparison results show that the proposed HMOPSO can yield better solutions than other reported approaches.Table 6
Optimal solutions forPL, VD, and Lmax⁡ minimization in IEEE 118-bus system.
Object
Initial
NSGA-II
MOEA/D
HMOPSO
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
BestPL
Best VD
BestLmax⁡
P
L (MW)
132.8629
124.1573
131.781
128.5268
121.458
134.647
131.332
120.0564
138.8756
132.4865
VD (pu)
1.4393
1.3478
0.5841
0.9013
1.255
0.5231
0.968
1.5536
0.4237
1.2524
L
max
⁡
0.0694
0.0569
0.0571
0.0552
0.0571
0.0586
0.0534
0.0582
0.0598
0.0512Table 7
Comparison of the best solutions of different algorithms for IEEE 118-bus system.
Object
Initial
GPAC [42]
CA [42]
APSODV [44]
BBO [12]
GSA [11]
OGSA [45]
HMOPSO
P
L (MW)
132.8629
131.9083
131.8638
—
128.97
122.6762
126.99
120.0564
VD (pu)
1.4393
1.2849
1.2756
0.4666
0.5020
0.5091
0.3666
0.4237
L
max
⁡
0.0694
—
—
0.0538
0.0523
0.0600
0.0512


## 5.3. Multiobjective Performance Metrics Analysis

In multiobjective optimization processes, two goals are normally taken into account:(1) convergence to the true Pareto-optimal set; (2) maintenance of diversity in solutions of the Pareto-optimal set. In order to obtain an accurate quantitative comparison of the performance of the proposed approach, this paper implements three different performance metrics: generational distance (GD) [30], minimal spacing (MSP) [31], and hypervolume (HV) [32]. A brief introduction of these metrics is given as follows.Generational Distance (GD). GD is used to evaluate the closeness of the nondominated set obtained by an algorithm to the reference Pareto-optimal front. Consider
(27)GD=∑i=1ndi2n,
where  n is the number of vectors in the reference Pareto-optimal front and di is the Euclidean distance (measured in objective space) between each of these and the nearest member of the Pareto optimal set.Minimal Spacing (MSP). MSP is used to evaluate how evenly the nondominated solutions are distributed in the objective. Consider
(28)MSP=1n-1∑i=1n(di-d¯)2,
where di=mink∑m=1Nobj|fmi-fmk|, Nobj is the number of objectives, i is the index of the solution which is marked as the current seed, and k is the index of the solution which is still unmarked at the current iteration.Hypervolume (HV). HV provides the combined qualitative information about closeness and diversity in obtained Pareto-optimal fronts. It calculates the volume covered by the members of the approximate Pareto-optimal Q for problems where all objectives are to be minimized. Mathematically, for each solution i∈Q, a hypercube vi is constructed with a reference point W and the solution i is the diagonal corners of the hypercube. The reference point can simply be found by constructing a vector of worst objective function values. Thereafter, a union of all hypervolume (HV) is calculated as follows:
(29)HV=volume(⋃i=1|Q|vi).In order to evaluate the performance metrics of these algorithms, each algorithm is executed 100 times for two test systems. The statistic comparison results of performance metrics for two test systems are shown in Table8. It is clear that the proposed algorithm can obtain better Pareto front with respect to other algorithms, because the average and standard deviation performances of GD, MSP, and HV in HMOPSO are the best.Table 8
Statistical results of performance metric for IEEE 30-bus and IEEE 118-bus systems.
Performance measures
IEEE 30-bus system
IEEE 118-bus system
NSGA-II
MOAE/D
MOPSO
HMOPSO
NSGA-II
MOAE/D
MOPSO
HMOPSO
GD
Mean
0.113627
0.015378
0.015862
0.012735
0.330574
0.022916
0.024811
0.018853
Std. dev.
0.037141
0.00414
0.003135
0.002283
0.064578
0.003867
0.005095
0.003456
MSP
Mean
0.263546
0.32752
0.289632
0.186752
0.287658
0.41543
0.325624
0.253628
Std. dev.
0.040234
0.01564
0.02466
0.011338
0.058462
0.003133
0.038236
0.031247
HV
Mean
0.185746
0.22676
0.211651
0.24836
0.132385
0.153652
0.165727
0.201435
Std. dev.
0.035541
0.031572
0.041701
0.025813
0.042356
0.055814
0.063217
0.040472Table9 shows the statistic results of computational time over 100 independent runs on the two test power systems. It is quite evident that the run time of HMOPSO is computationally more efficient than other techniques.Table 9
Computational time of different algorithm for IEEE 30-bus and 118-bus systems.
Time
IEEE 30-bus system
IEEE 118-bus system
NSGA-II
MOAE/D
MOPSO
HMOPSO
NSGA-II
MOAE/D
MOPSO
HMOPSO
Best (s)
168.42
54.34
40.16
23.52
405.78
167.59
142.65
110.17
Mean (s)
195.76
65.29
48.35
35.47
451.69
191.58
174.34
138.88
Worst (s)
213.84
78.67
55.37
41.78
510.36
215.28
191.28
165.54
Std. dev.
7.48
4.13
3.18
3.46
23.57
11.38
16.55
13.76


## 6. Conclusions

In this paper, a multiobjective ORPD problem with three conflicting objectives, such as real power loss minimization, bus voltage profile improvement, and voltage stability enhancement, is considered. To improve convergence and maintain good diversity, a new HMOPSO algorithm is proposed by incorporating Gaussian probability distribution, chaotic descending inertia weight, time variant acceleration coefficients, self-adaptive mutation operator, and dynamic crowding distance into the classical MOPSO. The proposed approach has been successfully applied to solve the ORPD problem on IEEE 30-bus and IEEE 118-bus systems under both normal and contingency states. The performance of the proposed approach is compared with existing popular MOEAs including NSGA-II and MOEA/D and other methods recently reported in the literature with respect to various multiobjective performance measures. The comparison results demonstrate the superiority of HMOPSO in terms of solution quality and computational efficiency and confirm its potential to solve the ORPD problem in the practical power system operation.



---

*Source: 124136-2014-06-18.xml*