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

, and


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, -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 -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.

Problem Formulation
where NE is the number of transmission lines;   is the conductance of the th line;   ∠  and   ∠  are the voltages at end bus  and bus  of the th line, respectively.

Voltage Deviation (VD).
The objective is to minimize the deviations in voltage magnitude at load buses that can be expressed as where  spec  is the prespecified reference value at load bus , which is usually set as 1.0 p.u., and NL is the number of load buses.

Voltage Stability Index (𝐿 max
).In this work, enhancing voltage stability can be achieved through minimizing the voltage stability indicator -index [1] values at every bus of the system and consequently the global power system indicator.This objective can be expressed as where  max is the global system indicator describing the stability of the complete system.The -index calculation for a power system is briefly discussed below.
For multinode system, By segregating the load buses (PQ) from generator buses (PV), (4) can be written as where   ,   and   ,   represent currents and voltages at generator buses and load buses.
Rearranging the above equation, we get where The -index of the th node is given by the expression Thus, a global system indicator describing the stability of the complete system is

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 position  best  ) and its neighbors' experiences (the best global position  best ).To extend PSO to handle multiobjectives, Coello uses the concept of Pareto dominance to decide the  best  and  best .The  best  of the  particle is the position which is nondominated by its own past positions.To determine the best global position  best 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]: where  is the index of iterations,   is the position of the particle , V  is the velocity of the particle ,  is the inertia weight,  1 and  2 are the acceleration constants,  1 and  2 are random numbers based on uniform probability distribution function in the range [0, 1], Rep  is the best global value of the particle  that is taken from the repository, and the index  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.

Gaussian Distributed Random Numbers.
In this paper, we employ Gaussian distribution sequences to generate the random numbers  1 and  2 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: where  1 and  2 are random numbers generated by Gaussian probability distribution with zero mean and unit variance.

Chaotic Inertia Weight.
In MOPSO, the parameter  controls the influence of the precious velocity on the present velocity.Suitable selection of  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  max to  min expressed as follows: 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: where cw  is a chaotic weight at iteration ,   is the weight factor from the IWA, and   is the chaotic parameter variable at iteration .Here,   is based on the Logistic map, whose formula is given by where  is a control parameter and  −1 denotes chaotic variable at iteration  − 1.Although the above equation is deterministic, it exhibits chaotic dynamics when  = 4, the initial  0 ∈ (0, 1), and  0 ̸ = {0, 0.25, 0.5, 0.75, 1.0}.

Time Variant Acceleration
Coefficients.In addition to , parameters  1 and  2 influence the performance of PSO.Higher values of  1 ensure larger deviation of the particle in the search space; the higher values of  2 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  1 is allowed to decrease linearly with iteration from  1, to  1, , and the value of  2 is allowed to increase linearly with iteration from  2, to  2, described as follows: 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: where  1 ,  2 ,  3 , and  4 are random integers uniformly between 1 and the number of the population of the swarm and

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: where , where  obj is the number of objectives;    ,   max , and   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 size  exceeds the maximum size limit   , the archive truncation procedure is invoked, which can be described as follows.
(2) Sort the nondominated solutions in   according to their DCDs in descending order.
(3) Remove the individual with the lowest DCD value from   .
(4) If  ≤   , stop the archive truncation; otherwise go to step 2 and continue.

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 where V   ,   , and    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 where Δ  and Δ   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(  ()).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.

Constraints Handling.
From Section 2, 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  best .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:

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: where  is the length of the time window;   is the maximal crowding distance at generation ;   is the average of   over  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 Figure 1.

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 size  pop and the maximum number of iterations Iter max are set at 50 and 100.The parameters of the stopping criteria  lim and  are set as 0.01 and 30.For the IEEE 118-bus system,  pop = 100, Iter max = 500,  lim = 0.02, and  = 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 weight  max = 0.7,  min = 0.4; acceleration constants  1, =  2, = 2.5,  1, =  2, = 0.5; chaotic control parameters  = 4,  0 = 0.48; (2) MOPSO:  = 0.729;  1 =  2 = 2.05; the number of grids per each dimension  grid = 50; (3) NSGA-II: crossover probability   = 0.9; mutation probability   = 1/, where  is the number of decision variables; the distribution indices for crossover and mutation operators as   = 20 and   = 20; (4) MOEA/D: the neighborhood size  size = 30; the probability that parent solutions are selected from the neighborhood   = 0.9; the parameters in DE operator CR = 1 and  = 0.5; mutation probability   = 1/; the distribution index used in the polynomial mutation   = 20.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.

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 Figure 2. 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   and VD minimization in IEEE 30-bus system.It is clear that the proposed HMOPSO can find better extreme solutions than other methods.
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.

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  3 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  max decreases and voltage stability is improved after the application of optimization and the results of the proposed HMOPSO are better than other two algorithms.

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.Table 5 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.

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 extreme solutions 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.

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 Paretooptimal 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 where  is the number of vectors in the reference Paretooptimal front and   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. where ,  obj is the number of objectives,  is the index of the solution which is marked as the current seed, and  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 Paretooptimal fronts.It calculates the volume covered by the members of the approximate Pareto-optimal  for problems where all objectives are to be minimized.Mathematically, for each solution  ∈ , a hypercube V  is constructed with a reference point  and the solution  is the diagonal corners of the hypercube.The reference point can simply be found by constructing a vector of worst objective function values.

Object
Initial GPAC [42] CA [42] APSODV [44] BBO [12] GSA [11] OGSA [ Thereafter, a union of all hypervolume (HV) is calculated as follows: 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 Table 8.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 9 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.

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.
min are the th objective of the th individual, of the maximum individual, and of the minimal individual, respectively, after sorting the population according to the th objective function values; CD  is the CD of the th individual; Var  is the variance of individuals which are neighbors of the th individual.
cos (  −   ) +   sin (  −   )] NB; NB is the number of buses;    and    are the generator real and reactive power, respectively;    and    are the load real and reactive power, respectively;   and   are the transfer conductance and susceptance between  and , respectively.  and    are voltage magnitude and reactive power output at the th generator bus, respectively;   is the tap setting of the th transformer;    is reactive power compensation of the th switched VAR source;  PQ  is the PQ bus voltage;    is the apparent power flow at the th Δ maxPQ  , Δ max   , and Δ max   are the maximal constraint violations of  PQ  ,    , and    achieved by any individual in the current population.
Initialize P best  for each particle in P 0 Determine G best  for each particle in P t Update the velocity and position for each particle in P t [40,41]pensators 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.05p.u.The detailed data are given in[40,41].

Table 1 :
Optimal settings of control variables for   and VD minimization in IEEE 30-bus system (Case 1).

Table 2 :
Comparison of best solutions of different algorithms for IEEE 30-bus system (Case 1).

Table 3 :
Optimal settings of control variables for   and  max minimization in IEEE 30-bus system (Case 2(a): heavy load).

Table 4 :
Optimal settings of control variables for   and  max minimization in IEEE 30-bus system (Case 2(b): contingency state).

Table 5 :
Optimal solutions of control variables for   , VD, and  max minimization in IEEE 30-bus system (Case 3).

Table 6 :
Optimal solutions for   , VD, and  max minimization in IEEE 118-bus system.

Table 7 :
Comparison of the best solutions of different algorithms for IEEE 118-bus system.

Table 8 :
Statistical results of performance metric for IEEE 30-bus and IEEE 118-bus systems.

Table 9 :
Computational time of different algorithm for IEEE 30-bus and 118-bus systems.