Secure multiobjective real and reactive power allocation of thermal power units

Electrical Power and Energy Systems 30 (2008) 594–602

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Secure multiobjective real and reactive power allocation of thermal power units Lakhwinder Singh a,*, J.S. Dhillon b a b

Department of Electrical Engineering, Baba Banda Singh Bahadur Engineering College, Fatehgarh Sahib 140 407, Punjab, India Department of Electrical and Instrumentation Engineering, Sant Longowal Institute of Engineering and Technology, Longowal 148 106, District Sangrur, Punjab, India

a r t i c l e

i n f o

Article history: Received 21 November 2006 Received in revised form 22 June 2008 Accepted 15 August 2008

Keywords: Multiobjective optimization Non-inferior solutions Best weight pattern evaluation Generalized generation shift distribution factors Generalized Z-bus distribution factors

a b s t r a c t This paper presents best weight pattern evaluation approach to solve multiobjective load dispatch (MOLD) problem which determines the allocation of power demand among the committed generating units, to minimize a number of objectives. Operating cost, minimal impacts on environment, active power loss, are the objectives undertaken to be minimized subject to physical and technological constraints. MOLD problem is decomposed in two stages and is solved sequentially. In first stage, a optimization problem having multiple objectives which are function of only active power generation like operating cost, gaseous pollutant emissions, is solved to get optimal dispatch of active power generation, subject to meet the active power demand, generators’ capacity constraint and transmission active power line flow limits. In second stage, the system real power loss which is a function of reactive power generation is minimized, to get optimal reactive power generation, subject to meet reactive power demand, reactive power generators’ capacity constraint and transmission reactive power line flow limits, when active power generation is known in prior from first stage. The validity of the proposed method is demonstrated on 11-bus, 17-lines and 30-bus, 41-lines IEEE power systems. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Many real-life optimization problems have several conflicting objective functions that should be minimized or maximized simultaneously. In such multiobjective optimization problems, finding the best possible solution means trading off between the different objectives. Instead of a single optimal solution, a set of compromise solutions, so called Pareto optimal solutions have achieved, where none of the objective function value can be improved without impairing at least one of the others [1]. The ultimate goal in multiobjective optimization is to seek a most preferred solution among the set of non-inferior solutions. Since, multiobjective framework has been motivated by a number of factors such as its great capacity of modeling and simple computational implementations through existing mathematical programming packages [2]. The increasing energy demand and decreasing energy resources have necessitated the optimum use of available resources. Economic dispatch is the optimization scheme intended to find the generation outputs that minimize the total fuel cost [3]. In recent years, rigid environment regulations [4,5] forced the utility planners to consider emission control as an important objective and to modify their operating strategies to reduce pollution and atmospheric emissions of the thermal power plants. Coal produces particulate * Corresponding author. Tel.: +91 1763 233213; fax: +91 1763 232113. E-mail addresses: [email protected] (L. Singh), [email protected] (J.S. Dhillon). 0142-0615/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2008.08.003

matter such as ash and gaseous pollutants apart from heat. The harmful ecological effects caused by emission of particulates and gaseous pollutants like oxides of nitrogen (NOx), oxides of sulphur (SOx) and oxides of carbon (COx) from thermal plants, which are of greater concern to power utility and communities, can be reduced by adequate distribution of load between the plants of a power system. Combined economic emission dispatch (CEED) problem has been solved with various evolutionary programming techniques considering line flow constraints by Venkatesh et al. [6]. A novel approach has been applied to real-time economic emission power dispatch by Huang and Huang [7]. Abido [8] has developed Pareto based multiobjective evolutionary algorithms for solving a multiobjective nonlinear optimization problem. Specifically, nondominated sorting genetic algorithm (NSGA), niched Pareto genetic algorithm (NPGA), and strength Pareto evolutionary algorithm (SPEA) have been developed and successfully applied to an environmental/economic electric power dispatch problem. Fuzzy based mechanism is employed to extract the best compromise solution over the trade-off curve. Traditionally, gaseous pollutant emissions are treated as constraints or as weighted part of the objective function. Weights are usually associated with emissions through societal prices representing the harmful degree of such emissions on society. Allowed maximum levels of gaseous pollutant emissions are set by regulating organizations based on this policy, while solving such multiobjective optimization problems. It suggests that each pollutant should be treated on its own merit in assigning cost values usually referred to as evaluating environmental

595

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

externalities. So, the solution to this global problem still needs attention of the concerned agencies. Pollution control agencies restrict the amount of emissions of pollutants depending upon their relative harmfulness to human beings. System security and reliability is another aspect of power system to provide the secure and adequate power supply. Arya et al. [9] have undertaken the transmission line security constrained economic dispatch problem and have solved using Davidon–Fletcher–Powell’s optimization method. Line flow constraints were taken into account by exterior penalty function method. Therefore, a priority structure can be formed for the multiobjective problem. Ideally, the utilities would like to supply power to its customers with minimum environmental emission as well as minimum total fuel cost with security and reliability in operation of power plants [10]. The intent of the paper is to solve multiobjective load dispatch (MOLD) problem which is a multiple non-commensurable objective problem, minimizes operating cost, minimal impacts on environment, active power and reactive power losses subject to physical and technological constraints. The MOLD problem is decomposed in two stages and solved in sequential manner. In first stage, multiobjective P-optimization problem and in second stage, scalar Q-optimization problem is solved, respectively. Basically, Poptimization problem minimizes multiple objectives which are functions of only active power generation viz. operating cost, gaseous pollutant emissions to get optimal dispatch of active power generation, subject to meet the active power demand, generators’ capacity constraint and active power transmission line flow limits. In second stage, Q-optimization problem is solved that minimizes the system power loss which is a function of reactive power generation viz. active power loss, to get optimal reactive power generation, subject to meet reactive power demand, reactive power generators’ capacity constraint and reactive power transmission line flow limits, when active power generation is known in prior from first stage. The best weight pattern is evaluated by conventional statistical measures, which characterize the correlation coefficients matrix evolution. Two fundamental concepts of multiple criteria decision methodology are undertaken, which are the contrast intensity of the alternative’s performance in each single objective and the conflicting nature of the objectives with each other. The extraction and exploitation of these two features is beneficial in the decision making. The active and reactive power line flows are obtained with the help of generalized generation shift distribution factors (GGDF) and generalized Z-bus distribution factors (GZBDF), respectively. The validity of the proposed method is demonstrated on 11-bus, 17-lines and 30-bus, 41-lines IEEE power systems, comprising of three and six generators, respectively. The results of the proposed approach are compared with the results of Ref. [2] in which evolutionary optimization technique has been used to search the ‘preferred’ weightage pattern to get the ‘best’ optimal solution in non-inferior domain. The results are also compared with the results of Ref. [8], in which Pareto based multiobjective evolutionary algorithms has been applied to solve a multiobjective nonlinear optimization problem.

2.1. Generalized generation shift distribution factors (GGDF) Using the definition of the reactance matrix and the dc approximations, the generation shift distribution factors (GSDF) [11] are obtained as

Aðm; iÞ ¼

oPTm X pi  X qi ¼ ; oP Gi Xm

The use of sensitivity factors to compute line flows in system security and contingency analysis remain very popular. The excellent properties of simplicity, linearity, accuracy and rapidity in computation make it widely acceptable in real-time applications. The generalized generation shift distribution factors (GGDF) and generalized Z-bus distribution factors (GZBDF) are two sensitivity factors in conventional approach. Both reflects the generation shift to line flows, can be used to calculate real and complex line flows after generation scheduling.

ð1Þ

where PGi is the active power generation at ith bus. PTm is the real line flow on the mth line. Xpi and Xqi are the elements of the reactance-bus matrix, while Xm is the line reactance of mth line. A(m, i) is the GSDF factor of ith generator and mth line. The new active power line flows after shifting can be expressed in the incremental form

PTm ¼ P0Tm þ

Ng X

Aðm; iÞDP Gi ;

m ¼ 1; 2; . . . ; NL

ð2Þ

i¼1

where P 0Tm is a base case of the real power flow on mth line. NL is number of lines in the system. DPGi is the small change in real power generation at ith bus. Since all generation changes can be absorbed by reference bus, the total system generation remains unchanged. The generalized generation shift distribution factors (GGDF) can completely replace GSDF and calculates line flows in an integral form. The line flows are defined by Ng X

PTm ¼

Dðm; iÞPGi ;

m ¼ 1; 2; . . . ; NL

ð3Þ

i¼1

where

Dðm; iÞ ¼ Dðm; RÞ þ Aðm; iÞ P P0Tm  Ng i¼1 Aðm; iÞP Gi i–R Dðm; RÞ ¼ PNg i¼1 P Gi

ð3aÞ ð3bÞ

where R denotes the reference bus. D(m, i) is the GGDF for ith generator and mth line. 2.2. Generalized Z-bus distribution factors (GZBDF) On the basis of the Z-bus concept, the Z-bus distribution (ZBD) factor [12] can be expressed as

Hðm; iÞ ¼

  oSTm ðZ pi  Z qi Þ ¼  Zm oSi

ð4Þ

where Si is the complex power injection at ith bus. STm is the complex line flow on the mth line. Zpi and Zqi are the elements of the Zbus matrix, while Zm is the line impedance of mth line. * denotes the complex conjugate. H(m, i) is the ZBD factor of ith generator and mth line. Based on Eq. (4) the complex power flow on transmission lines can be expressed in incremental form as

STm ¼ S0Tm þ 2. Line flows using sensitivity factors

m ¼ 1; 2; . . . ; NL

Nb X

Hðm; iÞDSi ;

m ¼ 1; 2; . . . ; NL

ð5Þ

i¼1

where S0Tm is a base case of the complex power flow on mth line. DSi is the small change in complex power injection at ith bus. Since, reference bus compensates any change in generation on rest of buses of power system, so the total system generation is assumed to be unchanged. A sensitivity factor, GZBDF [11] has been expressed in the paper to update the ZBD factor. Instead of an incremental form, an integral generation form is defined as

STm ¼

Ng X i¼1

Gðm; iÞSGi ;

m ¼ 1; 2; . . . ; NL

ð6Þ

596

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

where

Gðm; iÞ ¼ Gðm; RÞ þ Hðm; iÞ PNg S0Tm  i¼1 Hðm; iÞSGi i–R Gðm; RÞ ¼ PNg i¼1 SGi

ð6aÞ ð6bÞ

where R denotes the reference bus. G(m, i) is the GZBDF for ith generator and mth line. Eq. (6) can be rewritten as

PTm þ jQ Tm ¼

Ng X ½GP ðm; iÞ þ jGQ ðm; iÞðPGi þ jQ Gi Þ;

m ¼ 1; 2; . . . ; NL

i¼1

ð7Þ where PTm is real part of complex power flow through mth line. QTm is reactive part of complex power flow through mth line. GP(m, i) and GQ(m, i) are real and imaginary parts of GZBDF, respectively. Separating real and imaginary parts of Eq. (7), active and reactive power flow on transmission lines is obtained as

PTm ¼

Ploss ¼ B00 þ

Ng X ½GP ðm; iÞPGi  GQ ðm; iÞQ Gi Þ;

m ¼ 1; 2; . . . ; NL

ð7aÞ

Ng X ½GP ðm; iÞQ Gi þ GQ ðm; iÞP Gi ;

m ¼ 1; 2; . . . ; NL

ð7bÞ

i¼1

3. Multiobjective optimization problem formulation The multiobjective optimization problem is decomposed into two stages and solved sequentially. In first stage, multiobjective P-optimization problem minimizes four important non-commensurable objectives of an electrical thermal power system simultaneously subject to systems’ constraints. In second stage, Qoptimization problem minimizes scalar active power loss subject to equality and inequality constraints, when active power generation is known in prior from first stage.

Ng X

Bi0 PGi þ

i¼1

Ng X Ng X i¼1

PGi Bij P Gj

ð9Þ

j¼1

with

B00 ¼

i¼1

Q Tm ¼

where ai, bi and ci are the cost coefficients of ith generator. d2i, e2i and f2i are the NOx emission coefficients of ith generator. d3i, e3i and f3i are the SOx emission coefficients of ith generator. d4i, e4i and f4i are the COx emission coefficients of ith generator. PGi is the active power output of ith generator and act as a control variable. PDi is active and P max power demand at ith bus. Pmin Gi Gi are the minimum and maximum values of active power output of ith unit, respectively. max Pmin Tm and P Tm are the minimum and maximum limits of active power flow on mth transmission line, respectively. PTm is calculated from Eq. (3). B-coefficients are evaluated by performing load flow analysis on the system using decoupled load flow (DLF) method. Ploss is active power loss in the transmission lines [16] and is given as

Nb X Nb X i¼1

PDi Bij PDj ;

Bi0 ¼ 

j¼1

Nb X ðBij þ Bji ÞPDj j¼1

Pi ¼ PGi  PDi ; i ¼ 1; 2; . . . ; Nb Rij cosðhi  hj Þ Bij ¼ ; i ¼ 1; 2; . . . ; Nb; j ¼ 1; 2; . . . ; Nb jV i jjV j j cos ui cos uj Q hi ¼ di  ui ; i ¼ 1; 2; ::; Nb and ui ¼ tan1 i Pi di is voltage angle at ith bus. Nb is number of buses. To generate the non-inferior solutions, the multiobjective Poptimization problem is converted into scalar P-optimization problem using weighting method and is stated as

Minimize

L X

wj F j ðPGi Þ

ð10aÞ

j¼1

Subject to

ðiÞ

L X

wj ¼ 1:0;

wj P 0:0

ð10bÞ

j¼1

ðiiÞ Eqs: ð8eÞ—ð8gÞ

3.1. Multiobjective P-optimization problem In multiobjective P-optimization problem, operating cost, NOx, SOx and COx gaseous pollutant emissions which are functions of active power generation of thermal power stations are minimized, subject to meet the active power system demand, generators’ active power capacity constraints along with the limits of transmission active power line flows, to arrive at optimal dispatch of active power generation. The active power line flows are obtained with the help of GGDF. Mathematically P-optimization problem is stated as given below:

Minimize operating cost : F 1 ðP Gi Þ ¼

Ng X ðai P2Gi þ bi PGi þ ci Þ $=h ð8aÞ

where wj are the levels of the weighting coefficients. L is the total number of objectives. Fj(PGi); j = 1, 2, . . . , 4 are the objective functions to be minimized over the set of admissible decision vector, PGi; i = 1, 2, . . . , Ng. This approach yields meaningful results to the decision maker when solved many times for different simulated values of wj; j = 1, 2, . . . , L in such a manner that the sum of all weights is equal to one. To find the solution, constrained problem is converted into an unconstrained problem. Equality and inequality constraints are clubbed with objective functions using Lagrangian multiplier and penalty method, respectively. The generalized augmented function is defined as

!

i¼1

Minimize NOx emission : F 2 ðPGi Þ ¼

Ng X

ðd2i P2Gi

þ e2i PGi þ f2i Þ kg=h

LðPGi ; kp Þ ¼

XL j¼1

wj F j  kp

i¼1

ð8bÞ Minimize SOx emission : F 3 ðPGi Þ ¼

i¼1

ð8cÞ Minimize COx emission : F 4 ðP Gi Þ ¼

Ng 1 X þ k UP r 1 i¼1 Gi

Ng X ðd3i P2Gi þ e3i PGi þ f3i Þ kg=h

Ng X

ðd4i P2Gi

þ e4i PGi þ f4i Þ kg=h

i¼1

XNg

P i¼1 Gi

Pmin Gi 6 Pmin Tm 6

Nb X

PDi þ Ploss i¼1 PGi 6 Pmax i ¼ 1; 2; . . . ; Ng Gi ; PTm 6 P max ; m ¼ 1; 2; . . . ; NL Tm ¼

U PGi ¼

ð8eÞ ð8fÞ ð8gÞ

!

PGi 

Nb X

i¼1

PDi  Ploss

i¼1

NL 1 X þ k VP r 2 m¼1 Tm

!

where

ð8dÞ Subject to

Ng X

V PTm ¼

8 2 > min > ; P Gi < Pmin > Gi < P Gi  PGi

max 0; P min > Gi 6 P Gi 6 P Gi >  > : max 2 max PGi  PGi ; P Gi > PGi 8 2 > min min > > < PTm  PTm ; PTm < PTm

0; > > 2 > : PTm  Pmax ; Tm

max Pmin Tm 6 P Tm 6 P Tm

PTm > Pmax Tm

ð11Þ

597

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

where kp is Lagrangian multiplier, rk1 and r k2 are penalty factors and should have small values. The approximate Newton–Raphson method [17] is applied to obtain the non-inferior solutions for simulated weight combinations, to achieve the necessary conditions. 3.2. Scalar Q-optimization problem In second stage, Q-optimization problem minimizes active power loss [13] of the system subject to reactive power demand equality equation, reactive power generators capacity constraints as well as reactive power transmission line flow limits, to find reactive power dispatch when active power generation is fixed at optimal values obtained from first stage. The reactive power line flows are obtained with the help of GZBDF. Mathematically Q-optimization problem is stated as:

Minimize PL ¼ f ðP Gi ; Q Gi Þ Ng X

Subject to

Q Gi ¼

i¼1

ð12Þ

Nb X

Q Di þ Q L

ð12aÞ

i¼1

max Q min Gi 6 Q Gi 6 Q Gi ; max Q min Tm 6 Q Tm 6 Q Tm ; min max P Tm 6 PTm 6 PTm ;

i ¼ 1; 2; . . . ; Ng

ð12bÞ

m ¼ 1; 2; . . . ; NL

ð12cÞ

m ¼ 1; 2; . . . ; NL

ð12dÞ

where QGi is the reactive power output of ith generator. QDi is reactive power demand at ith bus. QL is reactive power loss in the transand Q max are the minimum and maximum mission lines. Q min Gi Gi values of reactive power output of ith unit, respectively. Q min Tm and Q max Tm are the minimum and maximum limits of reactive power flow on mth transmission line, respectively. PTm and QTm are calculated from Eqs. (7a) and (7b), respectively. The active and reactive power transmission losses, PL and QL are given by following equations:

PL ¼

Nb X Nb X ½Aij ðPi Pj þ Q i Q j Þ þ Bij ðQ i Pj  Pi Q j Þ i¼1

QL ¼

ð13Þ

j¼1

Nb X Nb X ½C ij ðP i Pj þ Q i Q j Þ þ Dij ðQ i Pj  Pi Q j Þ i¼1

ð14Þ

j¼1

where Pi + jQi = (PGi  PDi) + j(QGi  QDi)

Rij cosðdi  dj Þ; jV i jjV j j X ij C ij ¼ cosðdi  dj Þ; jV i jjV j j

Aij ¼

Rij sinðdi  dj Þ jV i jjV j j X ij Dij ¼ sinðdi  dj Þ jV i jjV j j

Bij ¼

with Aij, Bij, Cij and Dij are loss coefficients, and are evaluated from line data by performing load flow analysis. di and dj are load angles at ith and jth buses, respectively. Vi and Vj are voltage magnitudes at ith and jth buses, respectively. Rij and Xij are the real and reactive components of impedance bus matrix, respectively. To solve the problem, constrained problem is converted into an unconstrained problem. Equality and inequality constraints are clubbed with objective function using Lagrangian multiplier and penalty method, respectively. The augmented objective function is defined as

LðQ Gi ; kq Þ ¼ PL  kq

Ng X i¼1

þ

Q Gi 

Nb X

! Q Di  Q L

i¼1

NL NL X 1 X V PTm þ V Q Tm k r 4 m¼1 m¼1

!

þ

1 r k3

Ng X

! U Q Gi

8 2 ðP  Pmin > Tm Þ ; > < Tm V PTm ¼ 0; > > : 2 ðPTm  Pmax Tm Þ ; 8 2 ðQ  Q min > Tm Þ ; > < Tm V Q Tm ¼ 0; > > : 2 ðQ Tm  Q max Tm Þ ;

PTm < Pmin Tm max Pmin Tm 6 P Tm 6 P Tm

PTm > Pmax Tm Q Tm < Q min Tm max Q min Tm 6 Q Tm 6 Q Tm

Q Tm > Q max Tm

where kq is Lagrangian multiplier, r k3 and r k4 are penalty factors. Newton–Raphson method is applied to obtain the optimal solution [16]. 4. Evaluation of best weight pattern The fuzzy sets are defined by equations called membership functions, which represent the goals of each objective function. The membership function represents the degree of achievement of the original objective function as a value between 0 and 1 with l(Fi) = 1 as completely satisfactory and l(Fi) = 0 as unsatisfactory. By taking account of the lower and upper values of each objective function together with the rate of increase of membership satisfaction, the decision maker must determine the membership function l(Fi) in a subjective manner [14]. Graphical representation of the membership function is given in Fig. 1 and mathematically membership function is defined as

8 F i  F Li > > ; F Li 6 F i 6 F Oi > > O L > F  F > i i > < lðF i Þ ¼ 0; F i < F Li and F i > F Ui > > > > FO  F > > i > ; F Oi 6 F i 6 F Ui : Ui F i  F Oi

ð16Þ

where F Li ¼ F Oi  10% of F Oi . F Li and F Ui are the lower and upper limits of the ith objective function in which the solution is expected and F Oi is the desired minimum objective function. For every objective Fi, a membership function l(Fi) is defined which maps the values of Fi within the interval [0, 1]. This transformation is based on the concept of the ideal point, F Oi . So, the value l(Fi) expresses the degree to which the alternative is close to the ideal value F Oi , which is the best performance of ith objective and far from the anti-ideal points, F Li or F Ui , which gives the worst performance of ith objective. F Oi is achieved by performing minimization of ith scalar objective with respect to constraints of the problem. F Ui is achieved by exploiting conflicting nature of the objectives, i.e. if one objective will have minimum value, other objectives will have maximum values. F Li is further taken as 10%

µ ( Fi )

1

i¼1

ð15Þ

where

U Q Gi

8 min 2 min > > < ðQ Gi  Q Gi Þ ; Q Gi < Q Gi min ¼ 0; Q Gi 6 Q Gi 6 Q max Gi > > : 2 ðQ Gi  Q max Q Gi > Q max Gi Þ ; Gi

0 Fi

L

Fi

O

Fi

U

Fig. 1. Membership function.

Fi

598

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

deterioration on the lower side of the objective from the ideal point F Oi . The initial matrix of evaluations is converted into a matrix of relative scores with membership function l(Fi) and standard deviation, which quantifies the contrast intensity of the ith objective of the generated non-inferior solutions [15]. So, the standard deviation of l(Fi) is a measure of the value of that objective to the decision making process. A symmetric matrix of correlation coefficients is constructed, with dimension L  L and a generic element G[l(Fi), l(Fj)] which is the linear correlation coefficient between the membership functions and is

G½lðF i Þ; lðF j Þ ¼

Cov½lðF i Þ; lðF j Þ

rlðF i Þ rlðFj Þ

;

i ¼ 1; 2; . . . ; L; j ¼ 1; 2; . . . ; L

1. 2. 3.

4. 5. 6. 7. 8. 9.

ð17Þ 10.

Input the system data Simulate weight patterns varying in a step size of 0.01, satisfying Eq. (10b) Find lower and upper values for each of the participating objectives DO Simulate next weight patterns Compute initial values of k0P and P0Gi , when Ploss = 0 Compute value of fuel cost, F prevs cost DO Run decoupled load flow (DLF) Compute initial values of real line flow, P 0Tm Compute loss coefficients B00, Bi0, Bij and Ploss by using Eq. (9) DO Apply approximate Newton–Raphson method [17] to solve Eq. (11) new ; Pnew Modify the values, knew P Gi ; P Tm new Compute new value  of fuel cost, F cost  new 0  WHILE kP  kP 6 Error

where l(Fi) and l(Fj) are the membership functions of ith and jth objectives, respectively. rlðF i Þ and rlðF j Þ are the standard deviations of membership functions, l(Fi) and l(Fj), respectively. Cov [l(Fi), l(Fj)] is covariance between l(Fi) and l(Fj). It has been observed that the more discordant the scores of the alternatives in objective l(Fi) and l(Fj), lower are the values of P G[l(Fi), l(Fj)]. In this sense, the sum Lj¼1 ½1  GflðF i Þ; lðF j Þg represents a measure of the conflict created by objective l(Fj) with respect to the decision situation defined by the rest of the objectives. Information contained in multiple criteria decision method problems is related to both contrast intensity and conflict of the decision objective. Hence, the amount of information Ai emitted by the ith objective can be determined by composing the measures which quantify the two notions through the following multiplicative aggregation formula

Table 1 Fuel cost ($/h) equations

L X Ai ¼ rlðF i Þ ½1  GflðF i Þ; lðF j Þg;

F 21 ¼ 64:83P 21  079:027P 1 þ 028:8249 F 22 ¼ 31:74P 22  136:061P 2 þ 324:1775 F 23 ¼ 67:32P 23  239:928P 3 þ 610:2535

i ¼ 1; 2; . . . ; L; j ¼ 1; 2; . . . ; L

11. 12.

F 11 ¼ 38:66P 21 þ 641:031P 1 þ 303:7780 F 12 ¼ 21:82P 22 þ 742:890P 2 þ 847:1484 F 13 ¼ 13:45P 23 þ 830:154P 3 þ 274:2241

Table 2 NOx emission (kg/h) equations

j¼1

ð18Þ According to the analysis, the higher the value Ai, larger is the amount of information transmitted by the corresponding objective and the higher is its relative importance for the decision making process. ‘Best’ objective weights are evaluated by normalizing these values to unity according to the following equation:

Ai W i ¼ PL

j¼1 Aj

;

i ¼ 1; 2; . . . ; L

ð19Þ

To get the best compromised solution, Eq. (11) is solved by assigning the best evaluated objective weights using Approximate Newton–Raphson method. 5. Solution procedure P-optimization and Q-optimization problems are solved sequentially in two stages. P-optimization problem is a multiobjective optimization problem, so non-inferior solutions are generated using weighting method for simulated weights. Best weight evaluation approach is applied to get optimal dispatch of active power generation. In Q-optimization problem, optimal reactive power dispatch is obtained, when active power generation is known in prior from the solution of P-optimization problem and hence active power loss is minimized. 5.1. Solution procedure for P-optimization Stepwise procedure to compute best compromising solution for multiobjective P-optimization by evaluating best weight pattern is given below:

Table 3 SOx emission (kg/h) equations F 31 ¼ 23:20P 21 þ 384:624P 1 þ 182:2605 F 32 ¼ 12:84P 22 þ 445:647P 2 þ 508:5207 F 33 ¼ 08:13P 23 þ 497:641P 3 þ 165:3433

Table 4 COx emission (kg/h) equations F 41 ¼ 1:40053P 21  2:9952210P 1 þ 3:824770 F 42 ¼ 1:05929P 22  0:9552794P 2 þ 1:342851 F 43 ¼ 1:06409P 23  1:2736420P 3 þ 1:819625

Table 5 Line data Line

From

To

R (p.u.)

X (p.u.)

B (p.u.)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

1 1 2 2 2 3 4 4 4 5 5 7 7 8 8 8 10

9 11 3 7 10 4 6 8 9 6 9 8 10 9 10 11 11

0.15 0.05 0.15 0.10 0.05 0.08 0.10 0.10 0.15 0.12 0.05 0.05 0.08 0.12 0.08 0.10 0.12

0.50 0.16 0.50 0.28 0.16 0.24 0.28 0.28 0.50 0.36 0.16 0.16 0.24 0.36 0.24 0.28 0.36

0.030 0.010 0.030 0.020 0.010 0.015 0.020 0.020 0.030 0.025 0.010 0.010 0.015 0.025 0.015 0.020 0.025

599

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602 Table 6 Load data Bus

P (MW)

Q (MVAR)

1 2 3 4 5 6 7 8 9 10 11

0.25 0.25 0.25 0.25 0.25 0.10 0.40 0.90 0.70 0.25 0.25

0.05 0.05 0.05 0.05 0.05 0.02 0.10 0.45 0.35 0.05 0.05

14. Construct a symmetric matrix of correlation coefficients using Eq. (17) 15. Compute amount of information Ai emitted by the ith objective using Eq. (18) 16. Compute optimal values of the normalized weights by using Eq. (19) 17. The best set of decision vector is found by solving the following problem

Minimize

" L X

# wi F i

i¼1

Subject to Eqs: ð8eÞ—ð8gÞ 18. Stop

Table 7 Lower and upper values of objectives: 11-bus system F L1 F L2 F L3 F L4

5.2. Solution procedure for Q-optimization FU 1 FU 2 FU 3 FU 4

¼ 4510:5360 $=h ¼ 581:8995 kg=h ¼ 2706:2890 kg=h ¼ 6:247609 kg=h

¼ 4743:3120 $=h ¼ 763:3317 kg=h ¼ 2845:2960 kg=h ¼ 8:931189 kg=h

Stepwise procedure to compute optimal value of the active power loss of the system is given below:

   prevs  6 Error or (slack bus convergence) WHILE F new cost  F cost WHILE (more non-inferior solutions required) 13. Compute membership functions of the obtained non-inferior solutions using Eq. (16)

1. Input the system data DO 2. Run decoupled load flow (DLF) 3. Compute loss coefficients Aij, Bij, Cij, Dij and real power loss, P prevs L

Table 8 Non-inferior solutions: 11-bus system Weights w1

Objectives and membership functions w2

w3

w4

Objectives

F1 $/h

F2 kg/h

F3 kg/h

F4 kg/h

Membership

l(F1)

l(F2)

l(F3)

l(F4)

0.01

0.01

0.49

0.49

Objectives Membership

4510.7620 0.9990288

751.1226 0.0672929

2706.4230 0.9990358

6.8944730 0.7589548

0.49

0.49

0.01

0.01

Objectives Membership

4553.0250 0.8174692

640.1637 0.6788651

2731.5740 0.8181031

6.7287940 0.8206928

0.01

0.49

0.49

0.01

Objectives Membership

4580.8500 0.6979328

617.8611 0.8017904

2748.1880 0.6985837

6.9605970 0.7343146

0.49

0.01

0.01

0.49

Objectives Membership

4510.6240 0.9996224

755.9305 0.0407934

2706.3460 0.9995908

6.9244520 0.7477836

0.02

0.02

0.48

0.48

Objectives Membership

4511.1640 0.9973024

743.3949 0.1098858

2706.6550 0.9973673

6.8480960 0.7762367

0.48

0.48

0.02

0.02

Objectives Membership

4552.4170 0.8200808

640.7937 0.6753927

2731.2110 0.8207148

6.7246930 0.8222212

0.02

0.48

0.48

0.02

Objectives Membership

4578.7110 0.7012263

619.2229 0.7942845

2746.9110 0.7077710

6.9404380 0.7418267

0.48

0.02

0.02

0.48

Objectives Membership

4510.7780 0.9989617

750.6294 0.0700114

2706.4320 0.9989725

6.8909460 0.7602692

0.03

0.03

0.47

0.47

Objectives Membership

4511.7470 0.9947978

736.3723 0.1485922

2706.9950 0.9949207

6.8089220 0.7908343

0.47

0.47

0.03

0.03

Objectives Membership

4551.7950 0.8227531

641.4443 0.6718070

2730.8400 0.8233826

6.7205630 0.8237601

0.03

0.47

0.47

0.03

Objectives Membership

4576.5840 0.7016259

620.6243 0.7865602

2745.6410 0.7169057

6.9207100 0.7491779

0.47

0.03

0.03

0.47

Objectives Membership

4511.0320 0.9978688

745.4919 0.0983279

2706.5780 0.9979223

6.8599070 0.7718353

0.04

0.04

0.46

0.46

Objectives Membership

4512.4680 0.9917017

729.9639 0.1839133

2707.4200 0.9918647

6.7757910 0.8031800

0.46

0.46

0.04

0.04

Objectives Membership

4551.1620 0.8254717

642.1166 0.6681014

2730.4620 0.8261032

6.7164130 0.8253065

0.04

0.46

0.46

0.04

Objectives Membership

4574.4700 0.7253406

622.0664 0.7786120

2744.3780 0.7259929

6.9014290 0.7563627

0.46

0.04

0.04

0.46

Objectives Membership

4511.3780 0.9963837

740.5120 0.1257753

2706.7800 0.9964681

6.8312290 0.7825218

600

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

8.

Table 9 Correlation coefficients: 11-bus system

9.

Case I 1.0 0.9774342 0.9999996 0.1274389

0.9774342 1.0 0.9772668 0.0849545

0.9999996 0.9772668 1.0 0.1282230

0.1274389 0.0849545 0.1282230 1.0

Case II 1.0 0.9769117 0.9999995 0.0690144

0.9769117 1.0 0.976738 0.1457117

0.9999995 0.9767380 1.0 0.0698244

0.0690144 0.1457117 0.0698244 1.0

10. 11.

6. Test system and results The validity of the proposed approach is illustrated on 11-bus, 17-lines [18] and 30-bus, 41-lines [8] IEEE power systems, comprising of three and six generators, respectively.

Compute initial values of k0q ; Q 0Gi and Q 0Tm DO Compute Hessian and Jacobian matrix values IF (kQGRADk 6 Error) GOTO 10 Apply Gauss elimination method to solve Eq. (15) IF (kDQi + Dkqk 6 Error) GOTO 10 ; Q new and Q new Modify the values of knew q Gi Tm

4. 5. 6. 7.

Compute power loss QL by  using Eq. reactive  (14) PNb  P WHILE  Ng i¼1 Q Di  Q L  6 Error i¼1 Q Gi  Calculate power loss, P new L   real by using Eq. (13) new prevs   6 Error WHILE PL  P L Write optimal values of PGi, QGi and PL Stop

6.1. Results on 11-bus, 17-lines power system The operating cost, NOx emission, SOx emission, COx emission characteristics equations are given in Tables 1–4. Line data and load data is depicted in Tables 5 and 6, respectively. Lower and upper values of the objectives, Fi; i = 1, 2, . . ., L are obtained by giving full weightage to one of the objectives and

Table 10 Optimal values of weights, objectives and their membership functions: 11-bus system Cost ($nh)

NOx emission (kg/h)

SOx emission (kg/h)

COx emission (kg/h)

Case I Weights

0.2394

0.4157

0.2393

0.1056

Proposed method considering security constraint

Objectives Membership

4557.7480 0.7971787

635.5349 0.7043776

2734.3930 0.7978229

6.7616050 0.8084664

Proposed method ignoring security constraint

Objectives Membership

4558.1810 0.7953181

635.2082 0.7061784

2734.6520 0.7959594

6.7656030 0.8069765

Case II Weights

0.2418

0.4076

0.2417

0.1089

Proposed method considering security constraint

Objectives Membership

4557.8600 0.8042988

637.1951 0.6972401

2734.4650 0.8049388

6.7536030 0.8127114

Proposed method ignoring security constraint

Objectives Membership

4562.9170 0.7824761

632.5594 0.7226158

2737.4710 0.7832175

6.8205170 0.7877668

Weights

0.3800

0.2070

0.2070

0.2060

Objectives Membership

4523.5910 0.9439173

685.7549 0.4275802

2714.0210 0.9443773

6.6359370 0.8552948

Ref. [2]

Table 11 Real and reactive line flows at optimal values considering security constraint and ignoring security constraint, respectively: case I Line (m)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 a

Sending bus

1 1 2 2 2 3 4 4 4 5 5 7 7 8 8 8 10

Ending bus

9 11 3 7 10 4 6 8 9 6 9 8 10 9 10 11 11

Overloaded lines in the system.

P-optimization

Q-optimization

Considering security constraint

Ignoring security constraint

Considering security constraint

Ignoring security constraint

PTm

PTm

QTm

QTm

0.01304 0.02131 0.04658 0.09664 0.05969 0.02736 0.05728 0.16968 0.07938 0.11511 0.08499 0.01170 0.07423 0.01965 0.07610 0.07433 0.00869

0.01390 0.02376a 0.04676 0.09771 0.06176a 0.02588 0.05712 0.16994 0.07966 0.11524 0.08514 0.01229 0.07416 0.01946 0.07642 0.07317 0.00999a

0.29280 0.37366 0.01383 0.09750 0.03783 0.18958 0.02620 0.07935 0.03527 0.04072 0.01276 0.03922 0.11019 0.09439 0.11381 0.44943 0.31874

a

0.29954 0.38985 0.01685a 0.09361 0.02765 0.18436 0.02612 0.08425 0.03394 0.03940 0.01445 0.03561 0.11160 0.09662 0.11277 0.46211 0.33028a

Maximum real and reactive power flow through transmission lines

P max Tm

Q max Tm

0.014 0.022 0.051 0.106 0.060 0.030 0.063 0.187 0.087 0.127 0.093 0.013 0.082 0.022 0.075 0.082 0.090

0.295 0.429 0.014 0.102 0.041 0.206 0.027 0.106 0.039 0.044 0.016 0.043 0.122 0.106 0.117 0.508 0.320

601

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

neglecting the others. In case the weight assigned to objective is 1.0, it means full weightage is given to the objective and when weightage is zero, the objective is neglected. Obtained lower and upper values of the objectives are shown in Table 7. P-optimization problem is solved to find the real power generations. Real power generations are fixed and reactive power generations are obtained by solving Q-optimization problem successively. In case, a slack bus has generator, real power generation is modified during Qoptimization whereas such effect is not realized when no generator is available on a slack bus. So, the following two different cases are undertaken to demonstrate results of the proposed method: Case I Slack bus having no generator. Case II Slack bus having generator. In multiobjective P-optimization problem, weights are varied with a step size of 0.01 in such a manner that Eq. (10b) is satisfied. 16 weight patterns are simulated to generate the non-inferior solutions. Achieved non-inferior solutions are shown in Table 8. A symmetric matrix of correlation coefficients is constructed, with dimension 4  4 and the linear correlation coefficients between the membership functions obtained from generated non-inferior solutions are given in Table 9. The dimensions of correlation coefficients matrix depends on the number of participating objectives.

Table 12 Generation schedule corresponding to optimal solutions after P-optimization: 11-bus system PG1 (p.u.)

PG2 (p.u.)

PG3 (p.u.)

Ploss

Case I Considering security constraint Ignoring security constraint

1.2137 1.2147

1.6791 1.6788

1.1691 1.1683

0.2120 0.2118

Case II Considering security constraint Ignoring security constraint

1.1898 1.2261

1.7026 1.6779

1.1729 1.1595

0.2154 0.2135

Ref. [2]

1.4882

1.6359

0.9337

0.2078

The best compromising values of operating cost, NOx emission, SOx emission, COx emission and their membership functions corresponding to the best evaluated weight pattern are shown in Table 10. Results of the proposed method are compared with the results of Ref. [2] in which weight patterns are searched. It is concluded that better overall membership function value is achieved in the proposed method. Table 11 depicts real and reactive line flows at the optimal solution achieved considering security constraint as well as ignoring security constraint, respectively. Comparison of results in terms of generation schedule corresponding to optimal solutions after P-optimization and Q-optimization are shown in Tables 12 and 13, respectively. The proposed method is implemented when: (i) security constraints are considered and (ii) the security constraints are ignored. The best generation schedule and voltage profile are given in Table 14 for both the cases undertaken. It has been observed that multiobjective P-optimization and Q-optimization problems are to be repeatedly solved till no further improvement is achieved when generator bus is taken as slack bus. So, time taken to achieve the result is comparatively more in case II as compared to case I. However, in case I, inferior solution is obtained as compared to case II. 6.2. Comparison of results: 30-bus, 41-lines power system The proposed method is implemented on standard IEEE 30-bus, 41-lines [8] power system, comprising of six generators to solve economic load dispatch (ELD), minimum emission dispatch and economic-emission load dispatch (EELD) problems to meet the systems’ demand within generators capacity constraints and line flow limits. Operating cost and NOx gaseous pollutant emission is undertaken for study. The obtained results are compared with the optimal results obtained by non-dominated sorted genetic algorithm (NSGA) [8], niched Pareto genetic algorithm (NPGA) [8] and strength Pareto evolutionary algorithm (SPEA) [8]. The comparison of results is given in Tables 15 and 16. The proposed method gives better results. Due to conflicting nature of objectives,

Table 13 Generation schedule corresponding to optimal solutions after Q-optimization: 11-bus system QG1 (p.u.)

QG2 (p.u.)

QG3 (p.u.)

PL

QL

Case I Considering security constraint Ignoring security constraint

0.3312 0.4964

0.1651 0.0086

0.1866 0.0997

0.2956 0.2849

0.6221 0.6943

Case II Considering security constraint Ignoring security constraint

0.5628 0.5871

0.5423 0.3693

0.3223 0.4403

0.2535 0.2490

0.1574 0.1323

Ref. [2]

0.5170

0.6359

0.2507

0.2078

0.1297

Table 14 ‘Best’ generation schedule and voltage profile with security and without security constraints Bus

1 2 3 4 5 6 7 8 9 10 11

Case I

Case II

Real and reactive demand of the system

Considering security constraint

Ignoring security constraint

Considering security constraint

Ignoring security constraint

V

d

V

d

V

d

V

d

PD

QD

1.0700 1.0880 1.0620 0.9619 0.8941 0.9267 0.9929 0.9587 0.8989 1.0252 1.0700

0.0775 0.1392 0.1695 0.0175 0.1459 0.0907 0.0037 0.0537 0.1229 0.0258 0.00

1.0700 1.0880 1.0620 0.9618 0.8941 0.9266 0.9928 0.9586 0.8988 1.0251 1.0700

0.0765 0.1412 0.1711 0.0165 0.1454 0.0898 0.0024 0.0528 0.1225 0.0271 0.00

1.0700 1.0880 1.0620 0.9473 0.8763 0.9103 0.9743 0.9318 0.8811 1.0045 1.0099

0.00 0.0658 0.0940 0.0936 0.2265 0.1691 0.0771 0.1286 0.2027 0.0457 0.0635

1.0700 1.0880 1.0620 0.9478 0.8766 0.9108 0.9746 0.9321 0.8813 1.0048 1.0101

0.00 0.0542 0.0818 0.1028 0.2333 0.1772 0.0862 0.1364 0.2086 0.0547 0.0675

0.25 0.25 0.25 0.25 0.25 0.10 0.40 0.90 0.70 0.25 0.25

0.05 0.05 0.05 0.05 0.05 0.02 0.10 0.45 0.35 0.05 0.05

602

L. Singh, J.S. Dhillon / Electrical Power and Energy Systems 30 (2008) 594–602

Table 15 Results for minimum cost and minimum emission optimized individually: 30-bus system Optimal values

PG1 PG2 PG3 PG4 PG5 PG6 Cost NOx emission

Proposed method

Ref. [8]

Ignoring line flow limits

Considering line flow limits

Ignoring line flow limits

Cost

Cost

Cost

Emission

Emission

Emission

Considering line flow limits Cost

Emission

0.1110 0.3038 0.5869 0.9738 0.5217 0.3518

0.3978 0.5197 0.5020 0.3674 0.5939 0.4778

0.1116 0.3045 0.5882 0.9751 0.5236 0.3460

0.4066 0.5249 0.4871 0.3673 0.5943 0.4782

0.1152 0.3055 0.5972 0.9809 0.5142 0.3542

0.4101 0.4631 0.5435 0.3895 0.5439 0.5150

0.1475 0.3340 0.7864 1.0096 0.1072 0.4806

0.4693 0.5223 0.6479 0.4734 0.1784 0.5761

603.7090 0.2106

646.5431 0.1890

608.9587 0.2151

652.7971 0.1935

607.78 0.2199

645.22 0.1942

618.50 0.2302

654.14 0.2016

Table 16 Best compromise solution for proposed method: 30-bus system Optimal values

Proposed method Ignoring line flow limits

PG1 PG2 PG3 PG4 PG5 PG6 Cost NOx emission

Ref. [8] Considering line flow limits

Ignoring line flow limits

Considering line flow limits

NSGA

NPGA

SPEA

NSGA

NPGA

SPEA

0.1314 0.3145 0.5780 0.9332 0.5320 0.3598

0.1319 0.3154 0.5787 0.9346 0.5340 0.3544

0.2935 0.3645 0.5833 0.6763 0.5383 0.4076

0.2976 0.3956 0.5673 0.6928 0.5201 0.3904

0.2752 0.3752 0.5796 0.6770 0.5283 0.4282

0.2712 0.3670 0.8099 0.7550 0.1357 0.5239

0.2998 0.4325 0.7342 0.6852 0.1560 0.5561

0.3052 0.4389 0.7163 0.6978 0.1552 0.5507

604.0377 0.2076

609.2942 0.2122

617.80 0.2002

617.79 0.2004

617.57 0.2001

625.71 0.2136

630.06 0.2079

629.59 0.2079

while solving EELD problem, the proposed method gives better cost but deteriorated emission when line flow limits are ignored or considered. 7. Conclusions In this paper, multiobjective load dispatch optimization problem is solved in two stages sequentially: multiobjective P-optimization problem and scalar Q-optimization problem. In multiobjective P-optimization problem, real power generation is optimized within the system constraints while minimizing multiple objectives simultaneously. The active power line flows are obtained with the help of generalized generation shift distribution factors (GGDF). The best objective weights are calculated by conventional statistical measures, which characterize the correlation coefficients matrix evolution. Objective weights derived in multiobjective P-optimization problem are found to embody the information which is transmitted from all the objectives participating in the multiple criteria problem. In scalar Q-optimization problem reactive power generation is optimized while active power loss is minimized subject to meet reactive power demand within reactive power generators capacity constraints along with reactive power transmission line flow limits, when active power generation is kept at the fixed optimal values obtained from P-optimization problem. The active and reactive power line flows are obtained with the help of generalized Z-bus distribution factors (GZBDF) considering combined effect of active and reactive power generation. The proposed method needs few non-inferior solutions and hence gives computational efficiency. The validity of the proposed method is demonstrated on 11-bus and 30-bus systems. The proposed method fully exploits the dependency of objective functions on participating control variables due to which the proposed method also avoids the oscillations in attaining the optimum point.

References [1] Miettinen K, Makela MM. Synchronous approach in interactive multiobjective optimization. Eur J Operat Res 2006;170:909–22. [2] Brar YS, Dhillon JS, Kothari DP. Multiobjective load dispatch by fuzzy logic searching weightage pattern. Electr Power Syst Res 2002;63:149–60. [3] Jayabarathi T, Jayaprakash K, Jeyakumar DN, Raghunathan T. Evolutionary programming techniques for different kinds of economic dispatch problems. Electr Power Syst Res 2005;73:169–76. [4] Dhillon JS, Kothari DP. The surrogate worth trade-off approach for multiobjective thermal power dispatch problem. Electr Power Syst Res 2000;56:103–10. [5] Chaaban FB, Mezher T, Ouwayjan M. Options for emissions reduction from power plants: an economic evaluation. Electr Power Energy Syst 2004;26(1):57–63. [6] Venkatesh P, Gnanadass R, Pandhy NP. Comparison and application of evolutionary programming techniques to combined economic emission dispatch with line flow constraints. IEEE Trans Power Syst 2003;18(2):688–97. [7] Huang CM, Huang YC. A novel approach to real-time economic emission power dispatch. IEEE Trans Power Syst 2003;18(1):288–94. [8] Abido MA. Multiobjective evolutionary algorithms for electric power dispatch problem. IEEE Trans Evolut Comput 2006;10(3):315–29. [9] Arya LD, Choube SC, Kothari DP. A non-decomposed approach for security constrained economic dispatch. IE(I) Journal-EL 1997;77:229–34. [10] Somasundaram P, Kuppusamy K. Application of evolutionary programming to security constrained economic dispatch. Electr Power Energy Syst 2005;27:343–51. [11] Ng WY. Generalised generation distribution factors for power system security evaluations. IEEE Trans Power Appl Syst 1981;100(3):1001–5. [12] Lin CE, Chen ST, Hunag CL. A two-step sensitivity approach for real-time line flow calculation. Electr Power Syst Res 1991;21:63–9. [13] Malange FCV, Alves DA, da Silva LCP, Castro CA, da Costa GRM. Real power losses reduction and loading margin improvement via continuation method. IEEE Trans Power Syst 2004;19(3):1690–2. [14] Niimura T, Nakahima T. Multiobjective trade-off analysis of deregulated electricity transactions. Electr Power Energy Syst 2003;25:179–85. [15] Diakoulaki D, Mavrotas G, Papayannakis L. Determining objective weights in multiple criteria problems: the critic method. Comput Operat Res 1995;22(7):763–70. [16] Wood AJ, Wollenberg BF. Power generation operation and control. New York: John Wiley; 2003. [17] Dhillon JS, Kothari DP. Power system optimization. New Delhi: PHI; 2004. [18] Singh L, Dhillon JS. Fuzzy satisfying multiobjective thermal power dispatch based on surrogate worth trade-off method. Electr Power Compon Syst 2008;36(1):93–108.