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.

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

2. Problem Formulation

2.1. Problem Objectives
2.1.1. Real Power Loss

The objective is to minimize the real power loss in transmission lines expressed as where 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.

2.1.2. Voltage Deviation

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

2.1.3. Voltage Stability Index

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

2.2. Problem Constraints
2.2.1. Equality Constraints

These constraints represent load flow equation as where ; 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.

2.2.2. Inequality Constraints

The inequality constraints represent the system operating constraints: where 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;   is the bus voltage; is the apparent power flow at the th transmission line. The max and min represent the maximum and minimum values of the corresponding variables, respectively. , , , and are the numbers of generators, transformers, switchable VAR sources, and 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 position ) and its neighbors’ experiences (the best global position ). To extend PSO to handle multiobjectives, Coello uses the concept of Pareto dominance to decide the and . The of the particle is the position which is nondominated by its own past positions. To determine the best global position 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 , is the velocity of the particle , is the inertia weight, and are the acceleration constants, and are random numbers based on uniform probability distribution function in the range , 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.

3.2. Gaussian Distributed Random Numbers

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

3.3. 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 to 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 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 denotes chaotic variable at iteration . Although the above equation is deterministic, it exhibits chaotic dynamics when , the initial , and .

3.4. Time Variant Acceleration Coefficients

In addition to , parameters and influence the performance of PSO. Higher values of ensure larger deviation of the particle in the search space; the higher values of 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 is allowed to decrease linearly with iteration from to , and the value of is allowed to increase linearly with iteration from to 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 , , , and are random integers uniformly between 1 and the number of the population of the swarm and ; in other words, the indices are mutually different; , , , and are uniform distributed random numbers in the range .

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: where where is the number of objectives; , , and 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; is the of the th individual; is the variance of individuals which are neighbors of the th 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 size exceeds the maximum size limit , the archive truncation procedure is invoked, which can be described as follows.(1)Calculate DCDs of particles in the archive using (20).(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.

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 where  , , 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 . 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 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 . 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: where , , and are the maximal constraint violations of , , and 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: where is the length of the time window; is the maximal crowding distance at generation ; is the average of over generations; 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.

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 size and the maximum number of iterations are set at 50 and 100. The parameters of the stopping criteria and are set as 0.01 and 30. For the IEEE 118-bus system, , , , and . 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 , ; acceleration constants , ; chaotic control parameters , ;(2)MOPSO: ; ; the number of grids per each dimension ;(3)NSGA-II: crossover probability ; mutation probability , where is the number of decision variables; the distribution indices for crossover and mutation operators as and ;(4)MOEA/D: the neighborhood size ; the probability that parent solutions are selected from the neighborhood ; the parameters in DE operator and ; mutation probability ; the distribution index used in the polynomial mutation .

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 [] in discrete steps of 0.0125. The VAR compensators are installed at buses 10 and 24 within the interval [] 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 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 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.

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

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

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.

5.3. Multiobjective Performance Metrics Analysis

In multiobjective optimization processes, two goals are normally taken into account: convergence to the true Pareto-optimal set; 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 where   is the number of vectors in the reference Pareto-optimal 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. Consider where , 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 Pareto-optimal 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 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. 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.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This work was supported by the National High Technology Research and Development Program (2012AA050215) of the Ministry of Science and Technology of China.