- Email: [email protected]
Differential evolutionary algorithm for optimal reactive power dispatch
Electrical Power and Energy Systems 30 (2008) 435–441
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Differential evolutionary algorithm for optimal reactive power dispatch M. Varadarajan, K.S. Swarup * Department of Electrical Engineering, Indian Institute of Technology, Madras, Chennai 600 036, India
a r t i c l e
i n f o
Article history: Received 27 June 2007 Received in revised form 5 March 2008 Accepted 9 March 2008
Keywords: Evolutionary computation Differential evolution Optimal power flow Loss minimization Penalty function
a b s t r a c t This paper presents differential evolutionary algorithm for optimal dispatch for reactive power and voltage control in power system operation studies. The problem is formulated as a mixed integer, nonlinear optimization problem taking into account both continuous and discrete control variables. The optimal setting of control variables such as generator voltages, tap positions of tap changing transformers and the number of shunt reactive compensation devices to be switched for real power loss minimization in the transmission system are determined. In the proposed method, the inequality operational constraints were handled by ‘‘penalty parameterless” approach. This helps in avoiding the time-consuming trial and error process for fixing the penalty parameter and makes the process system independent. The algorithm was tested on standard IEEE 14,30,57 and 118-Bus systems and the results compared with conventional method. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Global optimization of a non-continuous, nonlinear function, arising from large-scale complex engineering problems, which may have a large number of local minima and maxima, is quite challenging. A number of deterministic approaches based on branch and bound and real algebraic geometry are found to be successful in solving these problems to some extend. Of late, stochastic and heuristic optimization techniques such as evolutionary algorithms (EA) have emerged as efficient tools for global optimization and have been applied to a number of engineering problems in diverse fields. For the secure and economic operation of largescale power systems, a variety of optimization problems have to be solved. The optimal power flow (OPF) problem, which was introduced in 1960s by Carpentier, [1] is an important and powerful tool for power system operation and planning. Reactive power optimization is a sub-problem of OPF calculation, which determines all the controllable variables, such as tap ratio of transformers, output of shunt capacitors/reactors, reactive power output of generators and static reactive power compensators etc., and minimizes transmission losses or other appropriate objective functions, while satisfying a given set of physical and operational constraints. Since transformer tap ratios and outputs of shunt capacitor/reactors have a discrete nature, while reactive power outputs generators, bus voltage magnitudes and angles are, on the other hand, continuous variables, the reactive power optimization problem is formulated as mixed-integer, nonlinear problem. * Corresponding author. Tel.: +91 44 2257 4440; fax: +91 44 2257 4401. E-mail address: [email protected] (K.S. Swarup). 0142-0615/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2008.03.003
A number of mathematical programming based techniques have been proposed to solve the OPF problem. For decades, conventional gradient-based optimization algorithms have been used for the solving optimal reactive power dispatch problem [2,3]. The gradient and Newton methods suffer from the difficulty in handling inequality constraints. Linear programming requires objective function and constraints have linear relationship, which may lead to loss of accuracy. Conventional methods are not efficient in handling problems with discrete variables. The combinatorial-search approaches, branch-and-bound and cutting plane algorithms, which are usually used to solve the mixed integer programming model, are ‘non-polynomial and all suffer from the problem of ‘‘curse of dimensionality” making them unsuitable for large-scale OPF problems. Wu and Ma [4] applied Evolutionary Programming (EP) for global optimization problems of large-scale power systems to achieve optimal reactive power dispatch and voltage control of power systems. Lai and Ma [5] showed that in optimization of non-continuous and non-smooth function, EP is much better than nonlinear programming and has applied it for reactive power planning. Lee et al. [6] solved the reactive power operational and investment-planning problem by using a Simple Genetic Algorithm (SGA) combined with the successive linear programming method. The Benders’ cut is constructed during the SGA procedure to enhance the robustness and reliability of the algorithm. Chebbo and Irving [7] proposed a linear programming based conventional approach for combined active and reactive power dispatch. Yoshida et al. [8] proposed a Particle Swarm Optimization (PSO) for reactive power and voltage control considering voltage security assessment. Zhao et al. [9] proposed a
436
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
solution to the reactive power dispatch problem with a PSO approach based on multi-agent systems. Differential evolutionary algorithm (DEA) is a technically simple; population based evolutionary algorithm (EA), which is highly efficient in constrained parameter optimization problems [10]. DEA employs a greedy selection process with implicit elitist features. It has demonstrated its robustness and effectiveness in a variety of applications, such as neural network learning and infinite impulse response filter design [11,12]. It presents no difficulty in solving mixed integer problems [13] and hence is highly suitable for reactive power optimization where the generator voltage is a real valued parameter while tap position and the number of shunt devices to be switched is integer parameters. DE differs from other EA’s in the mutation and recombination phase. Unlike stochastic techniques such as Genetic Algorithms (GA) and Evolutionary Strategies (ES), where perturbation occurs in accordance with a random quantity, DE uses weighted differences between solution vectors to perturb the population. Authors in [14–17] used differential evolution to solve problems in power systems. Method of constraint handling is extremely important in power system optimization problems. In all the previous works reported in literature, inequality constraints were handled by use of a penalty function approach, i.e., the constraint violation is multiplied by a penalty coefficient or parameter and added to the objective function. Deb [18] proposed a penalty parameterless scheme to overcome the difficulty of choosing penalty coefficients for GA based constrained optimization problems. It is important to realize that such a penalty parameterless strategy is only applicable to population based approach. This is because it requires the population to be divided into two sets: feasible and infeasible sets. The fitness function depends on the feasible and infeasible population members. Since in a point-by-point optimization approach, there is only one member in each iteration, such a penalty parameterless scheme cannot be applied. Although a penalty term is added to the objective function to penalize infeasible solutions, the method differs from the way the penalty term is defined in conventional methods and in earlier evolutionary algorithm implementations. In this paper an optimal reactive power dispatch using differential evolutionary algorithm with an efficient penalty parameterless scheme of constraint handling is employed. A performance comparison with conventional interior point technique is provided to highlight the efficiency of the differential evolutionary algorithm based optimization.
as capacitor and reactor banks, i.e., u ¼ ½P g ; V g ; t; Q c T . Of the control variable mentioned Pg and V g are continuous variables, while tap ratio, t, of tap changing transformers and reactive power output of compensation devices, Q c , are discrete in nature. Loss minimization is usually required when cost minimization is the main goal with generator active power generation as the control variable. When all control variables are utilized in a cost minimization, a subsequent loss minimization will not yield further improvements. Therefore in reactive power dispatch problem, such as loss minimization, active power generation of all generators, except slack generator, is fixed during the optimization procedure. 3. Problem formulation The solution of the optimal reactive power dispatch problem involves the optimization of the nonlinear objective function with nonlinear system constraints. 3.1. Objective function The objective function here is to minimize the active power loss (Ploss ) in the transmission system. Network losses either for the whole network or for certain sections are non-separable functions of dependent and independent variables. It is given as Ploss ¼
Nl X
g k ½ðtk V i Þ2 þ V 2j 2tk V i V j cos hij
where, Nl is the number of transmission lines; g k is conductance of branch k between buses i and j; t k the tap ratio of transformer k; V i is the voltage magnitude at bus i; hij the voltage angle difference between buses i and j. 3.2. Constraints The minimization of the above function is subjected to a number of equality and inequality constraints. The equality constraints are the power flow equations given by Pgi P di V i
Nb X
V j ðGij cos hij þ Bij sin hij Þ ¼ 0 for i ¼ 1; . . . ; N PV þ N PQ
j¼1
ð3Þ Q gi Q di þ Q ci V i
f ðx; uÞ gðx; uÞ ¼ 0
V j ðGij sin hij Bij cos hij Þ ¼ 0 fori ¼ 1; . . . ; N PQ ð4Þ
The optimal power flow (OPF) is a static, nonlinear optimization problem, which calculates a set of optimum variables from the network state, load data and system parameters. Optimal values are computed in order to achieve a certain goal such as generation cost or line transmission power loss minimization subject to equality and inequality constraints. The OPF problem can be presented as min
Nb X j¼1
2. Optimal power flow
s:t
ð2Þ
k¼1
ð1Þ
hðx; uÞ 6 0 where, f is the objective function that typically includes total generation cost, losses in transmission system etc. Generally, gðx; uÞ represents the load flow equations and hðx; uÞ represents transmission line limits and other security limits. The vector of dependent and control variables are denoted by x and u respectively. In general, the dependent vector includes bus voltage angles h, bus voltage magnitudes V L and generator reactive power Q g , i.e., x ¼ ½h; V L ; Q g T . The control variable vector consists of real power generation P g , generator terminal voltage V g , transformer tap ratio t and reactive power generation or absorption Q c of compensation devices such
where, Nb , N PV and N PQ are the number of buses, PV buses and PQ buses respectively; Gij , Bij are real and imaginary part of ði; jÞth element of bus admittance matrix; Pgi , Q gi are active and reactive power generation at bus i ; Pdi , Q di are active and reactive power demand at bus i ; Q ci the reactive power compensation source at bus i. The inequality constraints on security limits are given by max Pmin g;slack 6 P g;slack 6 P g;slack
V min 6 V Li 6 V max Li Li Q min gi
6 Q gi 6
jSl j 6 Smax l
Q max gi for
ð5Þ
for
i ¼ 1; . . . ; N PQ
ð6Þ
for
i ¼ 1; . . . ; N g
ð7Þ
l ¼ 1; . . . ; N l
ð8Þ
The inequality constraints on control variable limits are given by V min 6 V gi 6 V max gi gi
for
i ¼ 1; . . . ; N PV
ð9Þ
tmin 6 t k 6 t max k k
for
k ¼ 1; . . . ; N t
ð10Þ
i ¼ 1; . . . ; N c
ð11Þ
Q min 6 Q ci 6 Q max ci ci
for
where, N g , Nc and Nt are the number of generators, compensator devices and transformers; Sl the apparent power flow in line l ; Smax l
437
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
the maximum apparent power flow in line l; V gi the voltage magnitude at generator bus i; V Li the voltage magnitude at load bus i; min max Pg;slack the real power generation at slack bus; V min Li , V Li , Q gi , max min max min max min max min max Q gi , t k , t k , V gi , V gi , Q ci , Q ci , Pg;slack , and Pg;slack are minimum and maximum limits of the corresponding variables respectively. 4. Differential evolutionary algorithm
k 2 ½1; D
i 2 ½1; N p ; ð12Þ
Each variable k in a individual i in the generation G is initialized within its boundaries xkmin and xkmax . At each generation, two operators, namely mutation and crossover (recombination), are applied on each individual, thus producing the new population. Then, a selection phase takes place, where each individual of the new population is compared to the corresponding individual of the old population, and the best between them is selected as a member of the population in the next generation. In the following, the evolutionary operators are briefly described. 4.1. Mutation DE does not use a predefined probability density function to generate perturbing fluctuations. It relies upon the population itself to perturb the vector parameter. According to the mutation operator, for each individual, X Gi , i ¼ 1; . . . ; N p , at generation G, a Gþ1 Gþ1 mutant vector, V iGþ1 ¼ ðvGþ1 i1 ; vi2 ; . . . ; viD Þ is determined using the following equation V Gþ1 ¼ X Gr1 þ sðX Gr1 X Gr3 Þ i
X iGþ1 ¼
U Gþ1 i X Gi
if f ðU Gþ1 Þ < f ðX Gi Þ i
ð15Þ
otherwise
All solutions in the population have the same chance of being selected as parents independent of their fitness value. If the parent is still better, it is retained in the population, thus incorporating the feature of elitism. 4.4. Termination criteria
Differential Evolution Algorithm (DEA) is a simple populationbased, stochastic parallel search evolutionary algorithm for global optimization and is capable of handling non-differentiable, nonlinear and multi-modal objective functions. In DEA the population consists of real-valued vectors with dimension D that equals the number of design parameters. The size of the population is adjusted by the parameter N p . The initial population is uniformly distributed in the search space. xGi;k ¼ xkmin þ rand½0; 1 ðxkmax xkmin Þ
(
ð13Þ
where, s > 0 is a real parameter, called mutation constant, which controls the amplification of the difference between two individuals so as to avoid search stagnation and r 1 , r 2 and r 3 are mutually different integers, randomly selected from the set f1; 2; . . . ; N p g. 4.2. Crossover Following the mutation phase, the crossover (recombination) operator is applied on the population. For each mutant vector Gþ1 Gþ1 V Gþ1 , a trial vector, U iGþ1 ¼ ðui1 ; uGþ1 i i2 ; . . . ; uiD Þ is generated, with ( Gþ1 if ðrandj 6 CRÞ or ðj ¼ Irand Þ vj;i ð14Þ ¼ uGþ1 j;i xGj;i if ðrandj > CRÞ and ðj 6¼ Irand Þ where randj 2 ½0; 1 and Irand is chosen randomly from the interval ½1; . . . ; D once for each vector to ensure that at least one vector component originates from the mutated vector vi . CR is the DE control parameter that is called the crossover rate and is a user defined parameter within range [0,1]. Equation is applied for every vector component i 2 ½1; . . . ; N p , j 2 ½1; . . . ; D. 4.3. Selection To decide whether the vector U iGþ1 should be a member of the population comprising the next generation, it is compared to the corresponding vector X Gi .
The iterative procedure can be terminated when any of the following criteria is met, i.e., an acceptable solution has been reached, a state with no further improvement in solution is reached, control variables has converged to a stable state or a predefined number of iterations have been completed. In most of the cases, it is not easy to test whether the obtained solution is the most acceptable one. Here, in the present work, the iterations are terminated whenever the results remain constant for a fixed number of generations or maximum number of generations is reached, whichever occurs first. 5. Constraint handling 5.1. With penalty coefficients A common practice, for handling nonlinear constraints, is to use a penalty function approach because of its simplicity and ease of implementation. In most applications of EA to constrained optimization problems, this penalty function approach has been used. In this method for handling inequality constraints in minimization problems, the fitness function FðxÞ is defined as the sum of the objective function f ðxÞ and a penalty term which depends on the constraint violation hhðxÞi. FðxÞ ¼ f ðxÞ þ
n X
Rj hhj ðxÞi2
ð16Þ
j¼1
where, denotes the absolute value of the operand, if the operand is negative and returns a value zero, otherwise. The parameter Rj is the penalty parameter of the jth inequality constraint. The main difficulty arises here in choosing the appropriate penalty parameter R. The values of these parameters are chosen by trial and error method. To highlight the influence of the choice of penalty parameters on the objective function value, a simulation study is carried out on IEEE 14-Bus system [8]. The results obtained for three different cases are listed in Table 1. It can be observed that the real power loss value is greatly dependent on the choice of penalty parameters. 5.2. Without penalty coefficients In penalty parameter-less scheme, employed in this work, the composite fitness function for any f ðxÞ is given as follows f ðxÞ if x is feasible FðxÞ ¼ ð17Þ fmax þ CVðxÞ otherwise Here, fmax is the objective function value of the worst feasible solution in the population. In situation where none of the solutions in a
Table 1 Effect of penalty parameters on the solution Parameter
Case 1
Case 2
Case 3
R1 R2 R3 P loss (MW)
1000 100 100 13.42
2000 200 100 13.29
3000 300 100 13.32
438
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
population are feasible, fmax is not defined. Hence, such situations are handled by artificially inserting the base case solution into the population. CVðxÞ is the overall constraint violation of solution x. It is calculated as follows min CVðxÞ ¼ maxð0; P g;slack P max g;slack ; P g;slack P g;slack Þ
þ
N PQ X
min maxð0; V Li V max V Li Þ Li ; V Li
reactive power compensation capacitor/reactor banks at a bus is randomly generated from 0 to the number of existing equipment at the bus. This value is also modified in the search procedure, always limiting it between 0 and the number of existing compensator banks. 6.2. Algorithm
i¼1
þ
Ng X
min maxð0; Q gi Q max Q gi Þ þ gi ; Q gi
i¼1
Nl X
maxðjSl j Smax Þ l
l¼1
ð18Þ Since all feasible solution have zero constraint violation and all infeasible solutions are evaluated according to their constraint violations only, both the objective function value and constraint violation are not combined in any solution in the population. Thus there is no need to have any penalty coefficient R for this approach. The advantages of this scheme as compared to the usual penalty parameter based scheme are (i) The tedious process of choosing a suitable penalty coefficient R can be avoided, the inappropriate choice of which will affect the final solution and (ii) there is no need to evaluate the objective function value for individuals with constraint violation, which reduces the computation time. 6. Differential evolution for reactive power dispatch The implementation of differential evolutionary algorithm for reactive power dispatch consists mainly of identification objectives and handling of constraints. The control variables selected for reactive power optimization are the generator voltages except for slack generator; tap position of tap changing transformers and the number of shunt compensation banks to be switched. DEA is employed to find the best control variable setting starting from a randomly generated initial population of control variables. At the end of each generation, the best individuals, based on the fitness value, are stored. 6.1. Treatment of control variables In its basic form, differential evolutionary algorithm can handle only continuous variables. However, reactive power source installation and tap position of tap changing transformers are discrete variable or integer variable and generator terminal voltages are continuous variables in the reactive power dispatch problem. Here, differential evolutionary algorithm has been extended to handle mixed integer variables. For integer variables the value is rounded off to the nearest integer value of the variable. for continuous variables xi xi ¼ ð19Þ bxi c for integer variables The bxc function gives the nearest integer less than or equal to x. Different control variables in the optimal reactive power dispatch are treated as follows: initial generator terminal voltages are generated randomly between upper and lower limits of the voltage specification values. The value is then modified in the search procedure, within the specified limits. Transformer tap position is initially generated randomly between the minimum and maximum tap positions. The value is then modified in the search procedure among existing tap positions. Based on the tap position, the corresponding tap ratio is calculated as follows tk ¼ t min þ nk Dtk k
ð20Þ
where, nk is the number of tap positions and Dt k the step size of tap ratio. Using this tap ratio, the corresponding admittance of the transformer is determined for the load flow calculation. Initial number of
The details of the DE based optimization algorithm are as follows Step. 1 Generate an initial population randomly within the control variable bounds. Step. 2 For each individual in the population, run power flow algorithm such as Newton Raphson method, to find the operating points. Step. 3 Evaluate the fitness of the individuals according to Eqs. (17) and (18). Step. 4 Perform differentiation (mutation) and crossover as described in Sections 4.1 and 4.2 to create offspring from parents. Step. 5 Perform Selection as described in Section 4.3. between parent and offspring. While using the penalty parameterless method of constraint handling the following criteria are enforced while selecting the individuals for the next generation. Any feasible solution is preferred to any infeasible solution. Among two feasible solutions, the one having better objective function value is preferred. Among two infeasible solutions, the one having smaller constraint violation is preferred. Step. 6 Store the best individual of the current generation. Step. 7 Repeat steps 2 to 6 till the termination criteria is met (maximum number of generations).
7. Simulation results Results obtained by simulation of the algorithm using differential evolutionary algorithm, done in MATLAB 6.5 environment, are provided in this section. Simulation is carried out on IEEE 14 [8], IEEE 30, IEEE 57 and IEEE 118 [19] bus test systems. Floating point representation is used for encoding the control variables. There are several variants of DE. For this work, the DE scheme, which can be classified using notation DE=rand=1=bin strategy, is used. This strategy is the most popularly used one for solving practical problems [13]. The differential evolutionary computation parameters used for the simulation are summarized in Table 2. Number of individuals in a population for each test system is decided based on trial simulation run. The results, which follow, are the best solutions over 30 independent trials. 7.1. Case study of IEEE 14-Bus system In the modified IEEE 14-Bus system shown in Fig. 1, there are 14 buses, out of which 5 are generator buses. Bus 1 is the slack bus, 2, 3, 6 and 8 are taken as PV generator buses and the rest are PQ load buses. The network has 20 branches, 17 of which are transmission lines and 3 are tap- changing transformers. It is assumed that shunt compensation capacitor is available at buses 9 and 14 for voltage control. Totally, there are 9 control variables, which consists of 4 PV generator voltages, 3 tap changing transformers with 20 steps of 0.01 p.u. each and 2 shunt compensation capacitor banks with 3 steps of 0.06 p.u. each. The system data and initial operating conditions of the system are given in [8].
439
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441 Table 2 Simulation parameters used for differential evolution
1.12
DE control parameters
1.1 Step size ðsÞ ¼ 0:6 Crossover rate (CR) = 0.8
13 12
1.08
VG [p.u]
Population size = 30 Max. no. of generations = 200
VG8
1.06
VG2
VG6
1.04 VG3
14 11
1.02
10 6
1
8
9
1
7 5
0.98
4
0
10
20
30 40 Number of Generations
50
60
70
Fig. 3. Convergence of generator voltages for IEEE 14-Bus system.
2 3
25
Tap Position Number
20 Fig. 1. Network diagram of IEEE 14-Bus system.
15 10
Fig. 2 gives the performance of the optimization technique in terms of Ploss and Q loss . It can be observed that both P loss and Q loss get reduced over the evolutions and converge to a minimum value. The iterative procedure is terminated when there is no change in the result for 40 consecutive iterations or when 200 generations are reached, whichever occurs first. A good optimization results in convergence of all control variables to a steady state. Fig. 3 shows the variation of the continuous control variable, i.e., generator voltages V g , with respect to the number of generations. All four-generator voltages settle to a steady state before 30 generations. Fig. 4 shows the variation of the discrete control variables, tap position and capacitor bank switching. It can be observed that all discrete control variables also converge well before 30 generations. Table 3 gives the details of the control variables before and after optimization. The minimum and maximum limits of load bus volt-
T4–9
0 0
10
20
T5–6
T4–7
5 30
40
50
60
70
60
70
4 No. of Capacitor Banks Switched 3 Q
2
Q
C9
C14
1 0 0
10
20
30 40 Number of Generations
50
Fig. 4. Convergence of tap position and capacitor bank switching for IEEE 14-Bus system.
Table 3 Control variable setting and Ploss before and after optimization for IEEE 14-Bus system
Fig. 2. P loss and Q loss variation with generation for IEEE 14-Bus system.
Variable
Base Case
DE
PSO
IPM
V g2 V g3 V g6 V g8 T 47 T 49 T 56 Q c9 Q c14 P loss (MW)
1.0450 1.0100 1.0700 1.0900 0.9467 0.9524 0.9091 0.1800 0.1800 13.49
1.0449 1.0146 1.1000 1.1000 1.0600 1.0400 1.1000 0.1800 0.0600 13.239
1.0463 1.0165 1.1000 1.1000 0.9400 0.9300 0.9700 0.1800 0.0600 13.327
1.0449 1.0149 1.0971 1.0999 1.0238 1.0998 1.0550 0.1798 0.0739 13.246
ages are 0.95 and 1.1 p.u. respectively. From the base case value of 13.49 MW, the P loss was reduced to 13.239 MW with differential evolutionary computation. In order to evaluate the performance of differential evolutionary computation, the results were compared with popular Particle Swarm Optimization (PSO) [8] and conventional Interior Point Method (IPM), results of which are also presented in Table 3. OPF problem is modeled in AMPL [20] and
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
solved using the solver KINTRO [21] assuming all control variables to be continuous. To study the performance of the algorithm for different initial population, simulation is carried out with a constant population size of 30 on the test systems. Fig. 5 shows the effect of starting point on the final solution for five trial cases. It can be observed that even though the trajectory of convergence is different, all trials converge to the same final solution. As in the case of all random processes, it is useful to do repeated simulations and observe the solution statistics. For this particular case study, the simulations were repeated for 30 times from different starting points and the important statistical details are given in Table 4. It can be seen that the standard deviation of 30 trials is very small, guaranteeing a near best solution for any random trial.(See Table 4). It is essential to study the effect of population size on the optimization procedure. Even though a number of heuristic relations are available for finding the population size, there is no hard and fast rule, which can be universally adopted. In most cases, the population size is fixed by trial and error method. Fig. 6 shows the effect of population size on the solution, with all four trials having the same starting point. Except for a population size of 20, which has delayed convergence, all higher population sizes give more or less same speed of convergence. Hence, to achieve a compromise between fast convergence and reduced computation time, a population size of 30 was used for this study. 7.2. Case study of IEEE 30, 57 and 118-Bus systems In this section performance of differential evolutionary algorithm for optimal reactive power dispatch was evaluated on IEEE 30, 57 and 118 bus systems with simulation parameters given in Table 2. The network data for the test systems are given in [19]. As in the case of IEEE 14-Bus system, it was found that all of the control variables settled to a steady value when the termination
0.136 0.1355
population size 20 population size 30 population size 40 population size 50
0.135 Objective Value [p.u]
440
0.1345 0.134 0.1335 0.133 0.1325 0.132
0
10
20
30 40 Number of Generations
50
60
70
Fig. 6. Effect of population size on objective value.
Table 5 Ploss obtained by various methods for various test systems
Base Case DE PSO IPM a
IEEE 30-Bus
IEEE 57-Bus
IEEE 118-Bus
5.660 5.011 5.092a 5.101a
27.8637 25.0475 25.3047 26.0100
132.450 129.579 131.908a 132.110a
Taken from reference [22].
Table 6 Statistical details for IEEE 30,57 and 118-Bus systems Compared Item
IEEE 30-Bus
IEEE 57-Bus
IEEE 118-Bus
P loss – best (MW) P loss – worst (MW) P loss – average (MW) Standard deviation Average iteration Average CPU time (sec)
5.011 5.022 5.013 0.0026 66 13.647
25.0475 25.2016 25.1112 0.049 141 35.654
128.318 129.579 129.0817 0.345 193 42.1556
criterion is satisfied. It is assumed that tap changing transformers have 20 discrete steps of 0.01 p.u. each. The available reactive powers of capacitor banks are within the interval (0-30) MVAr in discrete steps of 1 MVAr. All bus voltages are required to be maintained within the range of 0.95-1.1 p.u. Table 5 shows the minimum P loss obtained by various methods namely DE, PSO and IPM. Differential evolutionary algorithm finds the minimum (optimal) value compared to other methods. To verify the robustness of the proposed methodology simulation is carried out for 30 trial runs with different initial population. The important statistical details are listed in the Table 6. To ensure a near optimum solution any random trial, the standard deviation should be very low. It can be seen that DEA is robust.
Fig. 5. Effect of starting point on objective value P loss .
8. Conclusion
Table 4 Statistical details for IEEE 14-Bus system Compared Item
DE
P loss – best (MW) P loss – worst (MW) P loss – average (MW) Standard deviation Average iteration Average CPU time (sec)
13.2396 13.2746 13.2506 1:6165 e2 62 8.172
A differential evolutionary computation based optimization method for reactive power dispatch was presented in this section. It has been observed that differential evolutionary algorithm is a simple but powerful tool for power system optimization problem with nonlinear objectives and constraints. Optimal control variables for real power loss minimization are determined using the differential evolutionary computation optimization strategy,
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
employing an efficient penalty parameterless method for constraint handling. Comparison with conventional techniques has verified the performance capability of differential evolutionary algorithm. References [1] Dommel HW, Tinney WF. Optimal power flow solutions. IEEE Trans Power App Syst 1968;PAS-87(10):1866–76. [2] Torres GL, Quintana VH. An interior-point method for nonlinear optimal power flow using rectangular coordinates. IEEE Trans Power Syst 1998;13(4):1211–8. [3] da Costa GRM, Costa CEU, de Souza AM. Comparative studies of optimization methods for the optimal power flow problem. Electric Power Syst Res 2000;56(3):249–54. [4] Wu QH, Ma JT. Power system optimal reactive power dispatch using evolutionary programming. IEEE Trans Power Syst 1995;10(3):1243–9. [5] Lai LL, Ma JT. Application of evolutionary programming to reactive power planning – comparison with nonlinear programming approach. IEEE Trans Power Syst 1997;12(1):198–206. [6] Lee KY, Bai X, Park YM. Optimization method for reactive power planning by using a modified simple genetic algorithm. IEEE Trans Power Syst 1995;10(4): 1843–50. [7] Chebbo AM, Irving MR. Combined active and reactive dispatch – Part 1: Problem formulation and solution algorithm. IEEE Proc Generat Transmission Distribut 1995;142(4):393–400. [8] Yoshida H, Fukuyama Y, Kawata K, Takayama S, Nakanishi Y. A particle swarm optimization for reactive power and voltage control considering voltage security assessment. IEEE Trans Power Syst 2001;15(4):1232–9. [9] Zhao B, Guo CX, Cao YJ. A multiagent-based particle swarm optimization approach for optimal reactive power dispatch. IEEE Trans Power Syst 2005;20(2):1070–8.
441
[10] Storn R, Price K. Differential evolution – a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical Report TR-95-012, ICSI; March 1995. [11] Storn R, Price K. Differential evolution A simple and efficient heuristic for global optimization over continuous spaces. J Global Optim 1997;11(4): 341–59. [12] Storn R. System design by constraint adaptation and differential evolution. IEEE Trans Evol Comput 1999;3(1):22–34. [13] Price K, Storn R, Lampinen J. Differential evolution: a practical approach to global optimization. Springer; 2005. [14] Chang CS, Lu LR, Wang F. Application of differential evolution in harmonic worst-case identification of mass rapid transit power supply system. Int J Syst Sci 2004;35(13–14):731–9. [15] Chiou JP, Chang CF, Su CT. Capacitor placement in large-scale distribution systems using variable scaling hybrid differential evolution. Int J Electric Power Energy Syst 2006;28(10):739–45. [16] Liang CH, Chung CY, Wong KP, Duan XZ. Parallel optimal reactive power flow based on cooperative co-evolutionary differential evolution and power system decomposition. IEEE Trans Power Syst 2007;22(1):249–57. [17] Chang CF, Chiou JP. Optimal feeder reconfiguration for loss reduction by robust searching hybrid differential evolution. Cybernetics and Systems 2007;38(4): 429–50. [18] Deb K. An efficient constraint handling method for genetic algorithms. Comput Meth Appl Mech Eng 2000;186(2–4):311–38. [19] System data (2006, December). Available from:. [20] Fourer R, Gay DM, Kernighan BW. AMPL: a modeling language for mathematical programming. Brooks Cole Publishing Company; 2002. [21] Waltz RA, Todd D. Plantenga TD. KNITRO 5.1 User Manual, Ziena Optimization Inc. [22] Vlachogiannis JG, Lee KY. A comparative study of particle swarm optimization for optimal steady state performance of power Systems. IEEE Trans Power Syst 2006;21(4):1718–28.
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Differential evolutionary algorithm for optimal reactive power dispatch M. Varadarajan, K.S. Swarup * Department of Electrical Engineering, Indian Institute of Technology, Madras, Chennai 600 036, India
a r t i c l e
i n f o
Article history: Received 27 June 2007 Received in revised form 5 March 2008 Accepted 9 March 2008
Keywords: Evolutionary computation Differential evolution Optimal power flow Loss minimization Penalty function
a b s t r a c t This paper presents differential evolutionary algorithm for optimal dispatch for reactive power and voltage control in power system operation studies. The problem is formulated as a mixed integer, nonlinear optimization problem taking into account both continuous and discrete control variables. The optimal setting of control variables such as generator voltages, tap positions of tap changing transformers and the number of shunt reactive compensation devices to be switched for real power loss minimization in the transmission system are determined. In the proposed method, the inequality operational constraints were handled by ‘‘penalty parameterless” approach. This helps in avoiding the time-consuming trial and error process for fixing the penalty parameter and makes the process system independent. The algorithm was tested on standard IEEE 14,30,57 and 118-Bus systems and the results compared with conventional method. Ó 2008 Elsevier Ltd. All rights reserved.
1. Introduction Global optimization of a non-continuous, nonlinear function, arising from large-scale complex engineering problems, which may have a large number of local minima and maxima, is quite challenging. A number of deterministic approaches based on branch and bound and real algebraic geometry are found to be successful in solving these problems to some extend. Of late, stochastic and heuristic optimization techniques such as evolutionary algorithms (EA) have emerged as efficient tools for global optimization and have been applied to a number of engineering problems in diverse fields. For the secure and economic operation of largescale power systems, a variety of optimization problems have to be solved. The optimal power flow (OPF) problem, which was introduced in 1960s by Carpentier, [1] is an important and powerful tool for power system operation and planning. Reactive power optimization is a sub-problem of OPF calculation, which determines all the controllable variables, such as tap ratio of transformers, output of shunt capacitors/reactors, reactive power output of generators and static reactive power compensators etc., and minimizes transmission losses or other appropriate objective functions, while satisfying a given set of physical and operational constraints. Since transformer tap ratios and outputs of shunt capacitor/reactors have a discrete nature, while reactive power outputs generators, bus voltage magnitudes and angles are, on the other hand, continuous variables, the reactive power optimization problem is formulated as mixed-integer, nonlinear problem. * Corresponding author. Tel.: +91 44 2257 4440; fax: +91 44 2257 4401. E-mail address: [email protected] (K.S. Swarup). 0142-0615/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2008.03.003
A number of mathematical programming based techniques have been proposed to solve the OPF problem. For decades, conventional gradient-based optimization algorithms have been used for the solving optimal reactive power dispatch problem [2,3]. The gradient and Newton methods suffer from the difficulty in handling inequality constraints. Linear programming requires objective function and constraints have linear relationship, which may lead to loss of accuracy. Conventional methods are not efficient in handling problems with discrete variables. The combinatorial-search approaches, branch-and-bound and cutting plane algorithms, which are usually used to solve the mixed integer programming model, are ‘non-polynomial and all suffer from the problem of ‘‘curse of dimensionality” making them unsuitable for large-scale OPF problems. Wu and Ma [4] applied Evolutionary Programming (EP) for global optimization problems of large-scale power systems to achieve optimal reactive power dispatch and voltage control of power systems. Lai and Ma [5] showed that in optimization of non-continuous and non-smooth function, EP is much better than nonlinear programming and has applied it for reactive power planning. Lee et al. [6] solved the reactive power operational and investment-planning problem by using a Simple Genetic Algorithm (SGA) combined with the successive linear programming method. The Benders’ cut is constructed during the SGA procedure to enhance the robustness and reliability of the algorithm. Chebbo and Irving [7] proposed a linear programming based conventional approach for combined active and reactive power dispatch. Yoshida et al. [8] proposed a Particle Swarm Optimization (PSO) for reactive power and voltage control considering voltage security assessment. Zhao et al. [9] proposed a
436
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
solution to the reactive power dispatch problem with a PSO approach based on multi-agent systems. Differential evolutionary algorithm (DEA) is a technically simple; population based evolutionary algorithm (EA), which is highly efficient in constrained parameter optimization problems [10]. DEA employs a greedy selection process with implicit elitist features. It has demonstrated its robustness and effectiveness in a variety of applications, such as neural network learning and infinite impulse response filter design [11,12]. It presents no difficulty in solving mixed integer problems [13] and hence is highly suitable for reactive power optimization where the generator voltage is a real valued parameter while tap position and the number of shunt devices to be switched is integer parameters. DE differs from other EA’s in the mutation and recombination phase. Unlike stochastic techniques such as Genetic Algorithms (GA) and Evolutionary Strategies (ES), where perturbation occurs in accordance with a random quantity, DE uses weighted differences between solution vectors to perturb the population. Authors in [14–17] used differential evolution to solve problems in power systems. Method of constraint handling is extremely important in power system optimization problems. In all the previous works reported in literature, inequality constraints were handled by use of a penalty function approach, i.e., the constraint violation is multiplied by a penalty coefficient or parameter and added to the objective function. Deb [18] proposed a penalty parameterless scheme to overcome the difficulty of choosing penalty coefficients for GA based constrained optimization problems. It is important to realize that such a penalty parameterless strategy is only applicable to population based approach. This is because it requires the population to be divided into two sets: feasible and infeasible sets. The fitness function depends on the feasible and infeasible population members. Since in a point-by-point optimization approach, there is only one member in each iteration, such a penalty parameterless scheme cannot be applied. Although a penalty term is added to the objective function to penalize infeasible solutions, the method differs from the way the penalty term is defined in conventional methods and in earlier evolutionary algorithm implementations. In this paper an optimal reactive power dispatch using differential evolutionary algorithm with an efficient penalty parameterless scheme of constraint handling is employed. A performance comparison with conventional interior point technique is provided to highlight the efficiency of the differential evolutionary algorithm based optimization.
as capacitor and reactor banks, i.e., u ¼ ½P g ; V g ; t; Q c T . Of the control variable mentioned Pg and V g are continuous variables, while tap ratio, t, of tap changing transformers and reactive power output of compensation devices, Q c , are discrete in nature. Loss minimization is usually required when cost minimization is the main goal with generator active power generation as the control variable. When all control variables are utilized in a cost minimization, a subsequent loss minimization will not yield further improvements. Therefore in reactive power dispatch problem, such as loss minimization, active power generation of all generators, except slack generator, is fixed during the optimization procedure. 3. Problem formulation The solution of the optimal reactive power dispatch problem involves the optimization of the nonlinear objective function with nonlinear system constraints. 3.1. Objective function The objective function here is to minimize the active power loss (Ploss ) in the transmission system. Network losses either for the whole network or for certain sections are non-separable functions of dependent and independent variables. It is given as Ploss ¼
Nl X
g k ½ðtk V i Þ2 þ V 2j 2tk V i V j cos hij
where, Nl is the number of transmission lines; g k is conductance of branch k between buses i and j; t k the tap ratio of transformer k; V i is the voltage magnitude at bus i; hij the voltage angle difference between buses i and j. 3.2. Constraints The minimization of the above function is subjected to a number of equality and inequality constraints. The equality constraints are the power flow equations given by Pgi P di V i
Nb X
V j ðGij cos hij þ Bij sin hij Þ ¼ 0 for i ¼ 1; . . . ; N PV þ N PQ
j¼1
ð3Þ Q gi Q di þ Q ci V i
f ðx; uÞ gðx; uÞ ¼ 0
V j ðGij sin hij Bij cos hij Þ ¼ 0 fori ¼ 1; . . . ; N PQ ð4Þ
The optimal power flow (OPF) is a static, nonlinear optimization problem, which calculates a set of optimum variables from the network state, load data and system parameters. Optimal values are computed in order to achieve a certain goal such as generation cost or line transmission power loss minimization subject to equality and inequality constraints. The OPF problem can be presented as min
Nb X j¼1
2. Optimal power flow
s:t
ð2Þ
k¼1
ð1Þ
hðx; uÞ 6 0 where, f is the objective function that typically includes total generation cost, losses in transmission system etc. Generally, gðx; uÞ represents the load flow equations and hðx; uÞ represents transmission line limits and other security limits. The vector of dependent and control variables are denoted by x and u respectively. In general, the dependent vector includes bus voltage angles h, bus voltage magnitudes V L and generator reactive power Q g , i.e., x ¼ ½h; V L ; Q g T . The control variable vector consists of real power generation P g , generator terminal voltage V g , transformer tap ratio t and reactive power generation or absorption Q c of compensation devices such
where, Nb , N PV and N PQ are the number of buses, PV buses and PQ buses respectively; Gij , Bij are real and imaginary part of ði; jÞth element of bus admittance matrix; Pgi , Q gi are active and reactive power generation at bus i ; Pdi , Q di are active and reactive power demand at bus i ; Q ci the reactive power compensation source at bus i. The inequality constraints on security limits are given by max Pmin g;slack 6 P g;slack 6 P g;slack
V min 6 V Li 6 V max Li Li Q min gi
6 Q gi 6
jSl j 6 Smax l
Q max gi for
ð5Þ
for
i ¼ 1; . . . ; N PQ
ð6Þ
for
i ¼ 1; . . . ; N g
ð7Þ
l ¼ 1; . . . ; N l
ð8Þ
The inequality constraints on control variable limits are given by V min 6 V gi 6 V max gi gi
for
i ¼ 1; . . . ; N PV
ð9Þ
tmin 6 t k 6 t max k k
for
k ¼ 1; . . . ; N t
ð10Þ
i ¼ 1; . . . ; N c
ð11Þ
Q min 6 Q ci 6 Q max ci ci
for
where, N g , Nc and Nt are the number of generators, compensator devices and transformers; Sl the apparent power flow in line l ; Smax l
437
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
the maximum apparent power flow in line l; V gi the voltage magnitude at generator bus i; V Li the voltage magnitude at load bus i; min max Pg;slack the real power generation at slack bus; V min Li , V Li , Q gi , max min max min max min max min max Q gi , t k , t k , V gi , V gi , Q ci , Q ci , Pg;slack , and Pg;slack are minimum and maximum limits of the corresponding variables respectively. 4. Differential evolutionary algorithm
k 2 ½1; D
i 2 ½1; N p ; ð12Þ
Each variable k in a individual i in the generation G is initialized within its boundaries xkmin and xkmax . At each generation, two operators, namely mutation and crossover (recombination), are applied on each individual, thus producing the new population. Then, a selection phase takes place, where each individual of the new population is compared to the corresponding individual of the old population, and the best between them is selected as a member of the population in the next generation. In the following, the evolutionary operators are briefly described. 4.1. Mutation DE does not use a predefined probability density function to generate perturbing fluctuations. It relies upon the population itself to perturb the vector parameter. According to the mutation operator, for each individual, X Gi , i ¼ 1; . . . ; N p , at generation G, a Gþ1 Gþ1 mutant vector, V iGþ1 ¼ ðvGþ1 i1 ; vi2 ; . . . ; viD Þ is determined using the following equation V Gþ1 ¼ X Gr1 þ sðX Gr1 X Gr3 Þ i
X iGþ1 ¼
U Gþ1 i X Gi
if f ðU Gþ1 Þ < f ðX Gi Þ i
ð15Þ
otherwise
All solutions in the population have the same chance of being selected as parents independent of their fitness value. If the parent is still better, it is retained in the population, thus incorporating the feature of elitism. 4.4. Termination criteria
Differential Evolution Algorithm (DEA) is a simple populationbased, stochastic parallel search evolutionary algorithm for global optimization and is capable of handling non-differentiable, nonlinear and multi-modal objective functions. In DEA the population consists of real-valued vectors with dimension D that equals the number of design parameters. The size of the population is adjusted by the parameter N p . The initial population is uniformly distributed in the search space. xGi;k ¼ xkmin þ rand½0; 1 ðxkmax xkmin Þ
(
ð13Þ
where, s > 0 is a real parameter, called mutation constant, which controls the amplification of the difference between two individuals so as to avoid search stagnation and r 1 , r 2 and r 3 are mutually different integers, randomly selected from the set f1; 2; . . . ; N p g. 4.2. Crossover Following the mutation phase, the crossover (recombination) operator is applied on the population. For each mutant vector Gþ1 Gþ1 V Gþ1 , a trial vector, U iGþ1 ¼ ðui1 ; uGþ1 i i2 ; . . . ; uiD Þ is generated, with ( Gþ1 if ðrandj 6 CRÞ or ðj ¼ Irand Þ vj;i ð14Þ ¼ uGþ1 j;i xGj;i if ðrandj > CRÞ and ðj 6¼ Irand Þ where randj 2 ½0; 1 and Irand is chosen randomly from the interval ½1; . . . ; D once for each vector to ensure that at least one vector component originates from the mutated vector vi . CR is the DE control parameter that is called the crossover rate and is a user defined parameter within range [0,1]. Equation is applied for every vector component i 2 ½1; . . . ; N p , j 2 ½1; . . . ; D. 4.3. Selection To decide whether the vector U iGþ1 should be a member of the population comprising the next generation, it is compared to the corresponding vector X Gi .
The iterative procedure can be terminated when any of the following criteria is met, i.e., an acceptable solution has been reached, a state with no further improvement in solution is reached, control variables has converged to a stable state or a predefined number of iterations have been completed. In most of the cases, it is not easy to test whether the obtained solution is the most acceptable one. Here, in the present work, the iterations are terminated whenever the results remain constant for a fixed number of generations or maximum number of generations is reached, whichever occurs first. 5. Constraint handling 5.1. With penalty coefficients A common practice, for handling nonlinear constraints, is to use a penalty function approach because of its simplicity and ease of implementation. In most applications of EA to constrained optimization problems, this penalty function approach has been used. In this method for handling inequality constraints in minimization problems, the fitness function FðxÞ is defined as the sum of the objective function f ðxÞ and a penalty term which depends on the constraint violation hhðxÞi. FðxÞ ¼ f ðxÞ þ
n X
Rj hhj ðxÞi2
ð16Þ
j¼1
where, denotes the absolute value of the operand, if the operand is negative and returns a value zero, otherwise. The parameter Rj is the penalty parameter of the jth inequality constraint. The main difficulty arises here in choosing the appropriate penalty parameter R. The values of these parameters are chosen by trial and error method. To highlight the influence of the choice of penalty parameters on the objective function value, a simulation study is carried out on IEEE 14-Bus system [8]. The results obtained for three different cases are listed in Table 1. It can be observed that the real power loss value is greatly dependent on the choice of penalty parameters. 5.2. Without penalty coefficients In penalty parameter-less scheme, employed in this work, the composite fitness function for any f ðxÞ is given as follows f ðxÞ if x is feasible FðxÞ ¼ ð17Þ fmax þ CVðxÞ otherwise Here, fmax is the objective function value of the worst feasible solution in the population. In situation where none of the solutions in a
Table 1 Effect of penalty parameters on the solution Parameter
Case 1
Case 2
Case 3
R1 R2 R3 P loss (MW)
1000 100 100 13.42
2000 200 100 13.29
3000 300 100 13.32
438
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
population are feasible, fmax is not defined. Hence, such situations are handled by artificially inserting the base case solution into the population. CVðxÞ is the overall constraint violation of solution x. It is calculated as follows min CVðxÞ ¼ maxð0; P g;slack P max g;slack ; P g;slack P g;slack Þ
þ
N PQ X
min maxð0; V Li V max V Li Þ Li ; V Li
reactive power compensation capacitor/reactor banks at a bus is randomly generated from 0 to the number of existing equipment at the bus. This value is also modified in the search procedure, always limiting it between 0 and the number of existing compensator banks. 6.2. Algorithm
i¼1
þ
Ng X
min maxð0; Q gi Q max Q gi Þ þ gi ; Q gi
i¼1
Nl X
maxðjSl j Smax Þ l
l¼1
ð18Þ Since all feasible solution have zero constraint violation and all infeasible solutions are evaluated according to their constraint violations only, both the objective function value and constraint violation are not combined in any solution in the population. Thus there is no need to have any penalty coefficient R for this approach. The advantages of this scheme as compared to the usual penalty parameter based scheme are (i) The tedious process of choosing a suitable penalty coefficient R can be avoided, the inappropriate choice of which will affect the final solution and (ii) there is no need to evaluate the objective function value for individuals with constraint violation, which reduces the computation time. 6. Differential evolution for reactive power dispatch The implementation of differential evolutionary algorithm for reactive power dispatch consists mainly of identification objectives and handling of constraints. The control variables selected for reactive power optimization are the generator voltages except for slack generator; tap position of tap changing transformers and the number of shunt compensation banks to be switched. DEA is employed to find the best control variable setting starting from a randomly generated initial population of control variables. At the end of each generation, the best individuals, based on the fitness value, are stored. 6.1. Treatment of control variables In its basic form, differential evolutionary algorithm can handle only continuous variables. However, reactive power source installation and tap position of tap changing transformers are discrete variable or integer variable and generator terminal voltages are continuous variables in the reactive power dispatch problem. Here, differential evolutionary algorithm has been extended to handle mixed integer variables. For integer variables the value is rounded off to the nearest integer value of the variable. for continuous variables xi xi ¼ ð19Þ bxi c for integer variables The bxc function gives the nearest integer less than or equal to x. Different control variables in the optimal reactive power dispatch are treated as follows: initial generator terminal voltages are generated randomly between upper and lower limits of the voltage specification values. The value is then modified in the search procedure, within the specified limits. Transformer tap position is initially generated randomly between the minimum and maximum tap positions. The value is then modified in the search procedure among existing tap positions. Based on the tap position, the corresponding tap ratio is calculated as follows tk ¼ t min þ nk Dtk k
ð20Þ
where, nk is the number of tap positions and Dt k the step size of tap ratio. Using this tap ratio, the corresponding admittance of the transformer is determined for the load flow calculation. Initial number of
The details of the DE based optimization algorithm are as follows Step. 1 Generate an initial population randomly within the control variable bounds. Step. 2 For each individual in the population, run power flow algorithm such as Newton Raphson method, to find the operating points. Step. 3 Evaluate the fitness of the individuals according to Eqs. (17) and (18). Step. 4 Perform differentiation (mutation) and crossover as described in Sections 4.1 and 4.2 to create offspring from parents. Step. 5 Perform Selection as described in Section 4.3. between parent and offspring. While using the penalty parameterless method of constraint handling the following criteria are enforced while selecting the individuals for the next generation. Any feasible solution is preferred to any infeasible solution. Among two feasible solutions, the one having better objective function value is preferred. Among two infeasible solutions, the one having smaller constraint violation is preferred. Step. 6 Store the best individual of the current generation. Step. 7 Repeat steps 2 to 6 till the termination criteria is met (maximum number of generations).
7. Simulation results Results obtained by simulation of the algorithm using differential evolutionary algorithm, done in MATLAB 6.5 environment, are provided in this section. Simulation is carried out on IEEE 14 [8], IEEE 30, IEEE 57 and IEEE 118 [19] bus test systems. Floating point representation is used for encoding the control variables. There are several variants of DE. For this work, the DE scheme, which can be classified using notation DE=rand=1=bin strategy, is used. This strategy is the most popularly used one for solving practical problems [13]. The differential evolutionary computation parameters used for the simulation are summarized in Table 2. Number of individuals in a population for each test system is decided based on trial simulation run. The results, which follow, are the best solutions over 30 independent trials. 7.1. Case study of IEEE 14-Bus system In the modified IEEE 14-Bus system shown in Fig. 1, there are 14 buses, out of which 5 are generator buses. Bus 1 is the slack bus, 2, 3, 6 and 8 are taken as PV generator buses and the rest are PQ load buses. The network has 20 branches, 17 of which are transmission lines and 3 are tap- changing transformers. It is assumed that shunt compensation capacitor is available at buses 9 and 14 for voltage control. Totally, there are 9 control variables, which consists of 4 PV generator voltages, 3 tap changing transformers with 20 steps of 0.01 p.u. each and 2 shunt compensation capacitor banks with 3 steps of 0.06 p.u. each. The system data and initial operating conditions of the system are given in [8].
439
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441 Table 2 Simulation parameters used for differential evolution
1.12
DE control parameters
1.1 Step size ðsÞ ¼ 0:6 Crossover rate (CR) = 0.8
13 12
1.08
VG [p.u]
Population size = 30 Max. no. of generations = 200
VG8
1.06
VG2
VG6
1.04 VG3
14 11
1.02
10 6
1
8
9
1
7 5
0.98
4
0
10
20
30 40 Number of Generations
50
60
70
Fig. 3. Convergence of generator voltages for IEEE 14-Bus system.
2 3
25
Tap Position Number
20 Fig. 1. Network diagram of IEEE 14-Bus system.
15 10
Fig. 2 gives the performance of the optimization technique in terms of Ploss and Q loss . It can be observed that both P loss and Q loss get reduced over the evolutions and converge to a minimum value. The iterative procedure is terminated when there is no change in the result for 40 consecutive iterations or when 200 generations are reached, whichever occurs first. A good optimization results in convergence of all control variables to a steady state. Fig. 3 shows the variation of the continuous control variable, i.e., generator voltages V g , with respect to the number of generations. All four-generator voltages settle to a steady state before 30 generations. Fig. 4 shows the variation of the discrete control variables, tap position and capacitor bank switching. It can be observed that all discrete control variables also converge well before 30 generations. Table 3 gives the details of the control variables before and after optimization. The minimum and maximum limits of load bus volt-
T4–9
0 0
10
20
T5–6
T4–7
5 30
40
50
60
70
60
70
4 No. of Capacitor Banks Switched 3 Q
2
Q
C9
C14
1 0 0
10
20
30 40 Number of Generations
50
Fig. 4. Convergence of tap position and capacitor bank switching for IEEE 14-Bus system.
Table 3 Control variable setting and Ploss before and after optimization for IEEE 14-Bus system
Fig. 2. P loss and Q loss variation with generation for IEEE 14-Bus system.
Variable
Base Case
DE
PSO
IPM
V g2 V g3 V g6 V g8 T 47 T 49 T 56 Q c9 Q c14 P loss (MW)
1.0450 1.0100 1.0700 1.0900 0.9467 0.9524 0.9091 0.1800 0.1800 13.49
1.0449 1.0146 1.1000 1.1000 1.0600 1.0400 1.1000 0.1800 0.0600 13.239
1.0463 1.0165 1.1000 1.1000 0.9400 0.9300 0.9700 0.1800 0.0600 13.327
1.0449 1.0149 1.0971 1.0999 1.0238 1.0998 1.0550 0.1798 0.0739 13.246
ages are 0.95 and 1.1 p.u. respectively. From the base case value of 13.49 MW, the P loss was reduced to 13.239 MW with differential evolutionary computation. In order to evaluate the performance of differential evolutionary computation, the results were compared with popular Particle Swarm Optimization (PSO) [8] and conventional Interior Point Method (IPM), results of which are also presented in Table 3. OPF problem is modeled in AMPL [20] and
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
solved using the solver KINTRO [21] assuming all control variables to be continuous. To study the performance of the algorithm for different initial population, simulation is carried out with a constant population size of 30 on the test systems. Fig. 5 shows the effect of starting point on the final solution for five trial cases. It can be observed that even though the trajectory of convergence is different, all trials converge to the same final solution. As in the case of all random processes, it is useful to do repeated simulations and observe the solution statistics. For this particular case study, the simulations were repeated for 30 times from different starting points and the important statistical details are given in Table 4. It can be seen that the standard deviation of 30 trials is very small, guaranteeing a near best solution for any random trial.(See Table 4). It is essential to study the effect of population size on the optimization procedure. Even though a number of heuristic relations are available for finding the population size, there is no hard and fast rule, which can be universally adopted. In most cases, the population size is fixed by trial and error method. Fig. 6 shows the effect of population size on the solution, with all four trials having the same starting point. Except for a population size of 20, which has delayed convergence, all higher population sizes give more or less same speed of convergence. Hence, to achieve a compromise between fast convergence and reduced computation time, a population size of 30 was used for this study. 7.2. Case study of IEEE 30, 57 and 118-Bus systems In this section performance of differential evolutionary algorithm for optimal reactive power dispatch was evaluated on IEEE 30, 57 and 118 bus systems with simulation parameters given in Table 2. The network data for the test systems are given in [19]. As in the case of IEEE 14-Bus system, it was found that all of the control variables settled to a steady value when the termination
0.136 0.1355
population size 20 population size 30 population size 40 population size 50
0.135 Objective Value [p.u]
440
0.1345 0.134 0.1335 0.133 0.1325 0.132
0
10
20
30 40 Number of Generations
50
60
70
Fig. 6. Effect of population size on objective value.
Table 5 Ploss obtained by various methods for various test systems
Base Case DE PSO IPM a
IEEE 30-Bus
IEEE 57-Bus
IEEE 118-Bus
5.660 5.011 5.092a 5.101a
27.8637 25.0475 25.3047 26.0100
132.450 129.579 131.908a 132.110a
Taken from reference [22].
Table 6 Statistical details for IEEE 30,57 and 118-Bus systems Compared Item
IEEE 30-Bus
IEEE 57-Bus
IEEE 118-Bus
P loss – best (MW) P loss – worst (MW) P loss – average (MW) Standard deviation Average iteration Average CPU time (sec)
5.011 5.022 5.013 0.0026 66 13.647
25.0475 25.2016 25.1112 0.049 141 35.654
128.318 129.579 129.0817 0.345 193 42.1556
criterion is satisfied. It is assumed that tap changing transformers have 20 discrete steps of 0.01 p.u. each. The available reactive powers of capacitor banks are within the interval (0-30) MVAr in discrete steps of 1 MVAr. All bus voltages are required to be maintained within the range of 0.95-1.1 p.u. Table 5 shows the minimum P loss obtained by various methods namely DE, PSO and IPM. Differential evolutionary algorithm finds the minimum (optimal) value compared to other methods. To verify the robustness of the proposed methodology simulation is carried out for 30 trial runs with different initial population. The important statistical details are listed in the Table 6. To ensure a near optimum solution any random trial, the standard deviation should be very low. It can be seen that DEA is robust.
Fig. 5. Effect of starting point on objective value P loss .
8. Conclusion
Table 4 Statistical details for IEEE 14-Bus system Compared Item
DE
P loss – best (MW) P loss – worst (MW) P loss – average (MW) Standard deviation Average iteration Average CPU time (sec)
13.2396 13.2746 13.2506 1:6165 e2 62 8.172
A differential evolutionary computation based optimization method for reactive power dispatch was presented in this section. It has been observed that differential evolutionary algorithm is a simple but powerful tool for power system optimization problem with nonlinear objectives and constraints. Optimal control variables for real power loss minimization are determined using the differential evolutionary computation optimization strategy,
M. Varadarajan, K.S. Swarup / Electrical Power and Energy Systems 30 (2008) 435–441
employing an efficient penalty parameterless method for constraint handling. Comparison with conventional techniques has verified the performance capability of differential evolutionary algorithm. References [1] Dommel HW, Tinney WF. Optimal power flow solutions. IEEE Trans Power App Syst 1968;PAS-87(10):1866–76. [2] Torres GL, Quintana VH. An interior-point method for nonlinear optimal power flow using rectangular coordinates. IEEE Trans Power Syst 1998;13(4):1211–8. [3] da Costa GRM, Costa CEU, de Souza AM. Comparative studies of optimization methods for the optimal power flow problem. Electric Power Syst Res 2000;56(3):249–54. [4] Wu QH, Ma JT. Power system optimal reactive power dispatch using evolutionary programming. IEEE Trans Power Syst 1995;10(3):1243–9. [5] Lai LL, Ma JT. Application of evolutionary programming to reactive power planning – comparison with nonlinear programming approach. IEEE Trans Power Syst 1997;12(1):198–206. [6] Lee KY, Bai X, Park YM. Optimization method for reactive power planning by using a modified simple genetic algorithm. IEEE Trans Power Syst 1995;10(4): 1843–50. [7] Chebbo AM, Irving MR. Combined active and reactive dispatch – Part 1: Problem formulation and solution algorithm. IEEE Proc Generat Transmission Distribut 1995;142(4):393–400. [8] Yoshida H, Fukuyama Y, Kawata K, Takayama S, Nakanishi Y. A particle swarm optimization for reactive power and voltage control considering voltage security assessment. IEEE Trans Power Syst 2001;15(4):1232–9. [9] Zhao B, Guo CX, Cao YJ. A multiagent-based particle swarm optimization approach for optimal reactive power dispatch. IEEE Trans Power Syst 2005;20(2):1070–8.
441
[10] Storn R, Price K. Differential evolution – a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical Report TR-95-012, ICSI; March 1995. [11] Storn R, Price K. Differential evolution A simple and efficient heuristic for global optimization over continuous spaces. J Global Optim 1997;11(4): 341–59. [12] Storn R. System design by constraint adaptation and differential evolution. IEEE Trans Evol Comput 1999;3(1):22–34. [13] Price K, Storn R, Lampinen J. Differential evolution: a practical approach to global optimization. Springer; 2005. [14] Chang CS, Lu LR, Wang F. Application of differential evolution in harmonic worst-case identification of mass rapid transit power supply system. Int J Syst Sci 2004;35(13–14):731–9. [15] Chiou JP, Chang CF, Su CT. Capacitor placement in large-scale distribution systems using variable scaling hybrid differential evolution. Int J Electric Power Energy Syst 2006;28(10):739–45. [16] Liang CH, Chung CY, Wong KP, Duan XZ. Parallel optimal reactive power flow based on cooperative co-evolutionary differential evolution and power system decomposition. IEEE Trans Power Syst 2007;22(1):249–57. [17] Chang CF, Chiou JP. Optimal feeder reconfiguration for loss reduction by robust searching hybrid differential evolution. Cybernetics and Systems 2007;38(4): 429–50. [18] Deb K. An efficient constraint handling method for genetic algorithms. Comput Meth Appl Mech Eng 2000;186(2–4):311–38. [19] System data (2006, December). Available from: