Reference Manual – MAINTENANCE - Airstream Pacific Northwest

13 Pages·2010·0.48 MB·English

HYDROLOGICALPROCESSES
Hydrol.Process.21,2897–2909(2007)
Publishedonline10January2007inWileyInterScience
(www.interscience.wiley.com)DOI:10.1002/hyp.6507
Multi-objective particle swarm optimization for generating
optimal trade-offs in reservoir operation
M. Janga Reddy and D. Nagesh Kumar*
DepartmentofCivilEngineering,IndianInstituteofScience,Bangalore-560012,India
Abstract:
A multi-objective particle swarm optimization (MOPSO) approach is presented for generating Pareto-optimal solutions for
reservoir operation problems. This method is developed by integrating Pareto dominance principles into particle swarm
optimization (PSO) algorithm. In addition, a variable size external repository and an efficient elitist-mutation (EM) operator
areintroduced.TheproposedEM-MOPSOapproachisfirsttestedforfewtestproblemstakenfromtheliteratureandevaluated
withstandardperformancemeasures.ItisfoundthattheEM-MOPSOyieldsefficientsolutionsintermsofgivingawidespread
of solutions with good convergence to true Pareto optimal solutions. On achieving good results for test cases, the approach
wasappliedtoacasestudyofmulti-objectivereservoiroperationproblem,namelytheBhadrareservoirsysteminIndia.The
solutionsofEM-MOPSOsyieldatrade-offcurve/surface,identifyingasetofalternativesthatdefineoptimalsolutionstothe
problem.Finally,tofacilitateeasyimplementationforthereservoiroperator,asimplebuteffectivedecision-makingapproach
was presented. The results obtained show that the proposed approach is a viable alternative to solve multi-objective water
resourcesandhydrologyproblems.Copyright2007John Wiley&Sons,Ltd.
KEYWORDS multi-objectiveoptimization;particleswarmoptimization;elitist-mutation;reservoiroperation;hydropower;
irrigation;water quality;Pareto optimalsolutions
Received2December2005;Accepted26June2006
INTRODUCTION approaches developed to deal with multiple objectives,
tradeoff methodologies have shown promise as effec-
Most of the water resources and hydrology problems are
tive means for considering non-commensurate objectives
characterized by multiple objectives and/or goals, which
that are to be subjectively compared in operation deter-
often conflict and compete with one another. Optimiza-
mination (Cohon and Marks, 1975). Therefore efficient
tion of multi-purpose reservoir systems involves solving
multi-objective problems. For example, for a reservoir generation of a set of alternatives for multiple objectives
system having hydropower and flood control as key pur- is very important with minimum computational require-
poses, the two major objectives can be maximization of ments.
the hydropower generation from the reservoir and mini-
mization of flood risk or flood damage. Obviously, these Scopefor multi-objective optimization using
two objectives are in conflict and compete with each meta-heuristic techniques
other. The higher the level of the reservoir, the more
Reservoir operation modeling has an exhaustive liter-
the hydropower generation possible because of the high
water head, yet less water storage will be available for aturepresenting variousoptimization techniques in order
flood control purposes and vice-versa. Clearly, one can to solve various kinds of problems (Yeh, 1985). How-
identify, within the active storage capacity of that reser- ever, in order to focus on the goal of this paper, a
voir,aParetooptimum regionwheretheenhancementof brief overview is given here. Most of the researchers
the first objective canbe achievedonly at the expense or on reservoir operation problems have tried conventional
degradationof the second, namelyflood control (Haimes techniques to generate tradeoffs among multiple objec-
etal., 1990). Also the units of these two objectives are tives. For example, Tauxe etal. (1979) have applied a
non-commensurable. The first objective, which maxi- multi-objective dynamic programming (DP) model for
mizes the hydroelectric power, is generally measured in analyzing a reservoir operation problem involving three
units of energy and not necessarily in monetary units, conflicting objectives. Thampapillai and Sinden (1979)
whereas the second objective can be measured in terms and Mohan and Raipure (1992) analyzed the tradeoffs
of acres of land, livestock, or human lives saved. If the for multiple objective planning through linear program-
objectives are non-commensurate, the classic methods ming (LP). To handle multiple objectives, many studies
of optimization cannot be applied easily. Of the several have used either the weighing approach or the constraint
method. The constraint method was used for generation
ofnoninferiorsetandtrade-offcurvesforreservoiroper-
*Correspondence to: D.NageshKumar,Department ofCivilEngineer-
ation problems (e.g. Croley and Rao, 1979; Yeh and
ing,IndianInstituteofScience,Bangalore-560012,India.
E-mail:[email protected] Becker, 1982; Liang etal., 1996 and Wang et al., 2005).
Copyright2007JohnWiley&Sons,Ltd.
2898 M.J.REDDYANDD.NAGESHKUMAR
The conventional optimization methods such as DP, PSO (MOPSO) is developed. To demonstrate the effi-
LP, and non-linear programming (NLP) are not suitable ciency of the proposed approach, results obtained are
to solve multi-objective optimization problems (MOOP), compared with NSGA-II, and evaluated with standard
becausethesemethodsuseapoint-by-pointapproach,and performance measures that are frequently used for per-
the outcome of these classical optimization methods is a formance evaluation of MOEAs.
single optimal solution. For example, the weighted sum The proposed approach for solving a multi-objective
method will convert the MOOP into a single objective decision problem in reservoir operation has a great
optimization.Byusingasinglepairoffixedweights,only potential for application, due to its attractive feature of
one point on the Paretofrontcanbe obtained. Therefore, generationoflargenumberofwellspreadParetooptimal
if one would like to obtain the global Pareto optimum, solutionsinasinglerun.Theotherapproachessuggested
all possible Pareto fronts must first be derived. This inthispaperfordecisionmaking,provideanopportunity
requires the algorithms to be executed iteratively, so as forthereservoiroperatortochoosethedesiredalternative
to ensure that every weight combination has been used. from a set of Pareto-optimal solutions.
Obviously, it is impractical to reiterate the algorithms
continuallytoexhaustalltheweightcombinations.Hence
PARTICLE SWARM OPTIMIZATION
the algorithms should have an ability to ‘learn’ from
previous performance to direct the proper selection of Swarm intelligence is a new area of research, from
weightsinfurtherevolutions.Alsoconventionalmethods which the PSO technique has been evolved through
may face problems, if the optimal solution lies on non- a simple simulation model of the movement of social
convex or disconnected regions of function space (Deb, groups such as birds and fish (Kennedy and Eberhart,
2001). 2001). The basis of this algorithm is that local inter-
Recently, meta-heuristic techniques such as evolu- actions motivate the group behaviour, and individual
tionary algorithms (EAs) and swarm intelligence tech- members of the group can profit from the discoveries
niques are becoming increasingly popular for solving and experiences of other members. Social behaviour is
optimization problems. In the recent past, evolution- modeled in PSO to guide a population of particles (so-
ary techniques have been successfully applied for single called swarm), moving towards the most promising area
objectiveoptimization(OliveiraandLoucks,1997;Ward- of the search space. The changes of the position of
law and Sharif, 1999; Raju and Nagesh Kumar, 2004) the particles within the search space are based on the
and MOOPs (Yapo etal., 1998; Vrugt etal., 2003; Khu social psychological tendency of individuals to emulate
and Madsen, 2005), due to their efficiency and ease in the success of other individuals. In PSO, each particle
handlingnon-linearandnon-convexrelationshipsofreal- represents a candidate solution. If the search space is
worldproblems.Thesetechniqueshavesomeadvantages D-dimensional, the ith individual (particle) of the popu-
over the classical optimization techniques (Deb, 2001). lation (swarm) can be represented by a D-dimensional
They use a population of solutions in each iteration and vector, Xi D(cid:1)xi1,xi2,...,xiD(cid:2)T. The velocity (position
offer a set of alternatives in a single run. They use ran- change) of this particle, can be represented by another
domizedinitializationandstochasticsearchintheiroper- D-dimensionalvector,Vi D(cid:1)vi1,vi2,...,viD(cid:2)T.Thebest
ation. Therefore, they can locate the search at any place previouslyvisitedpositionoftheith particleisdenotedas
over the entire search space and are able to overcome Pi D(cid:1)pi1,pi2,...,piD(cid:2)T. Defining g as the index of the
the problems of local optima. Thus population based globalguideoftheparticleintheswarm,andsuperscripts
stochasticsearchtechniquesaremoreappropriatetosolve denoting the iteration number, the swarm is manipulated
MOOPs. Achieving a well-spread and diverse Pareto according to the following two equations:
solution front is the primary goal of MOOP. Among the vnC1 D(cid:3)[w vn Cc rn(cid:1)pn (cid:1)xn(cid:2)/t
elitistmulti-objectiveEAs(MOEAs),strengthParetoEA id id 1 1 id id
(SPEA) (Zitzler and Thiele, 1999), Pareto-archived evo- Cc2r2n(cid:1)pngd(cid:1)xind(cid:2)/t] (cid:1)1(cid:2)
lutionary strategy (PAES) (Knowles and Corne, 2000)
xnC1 Dxn Ct vnC1 (cid:1)2(cid:2)
and non-dominated sorting genetic algorithm (NSGA-II) id id id
(Deb etal., 2002) have been successfully demonstrated where dD1,2,...,D; iD1,2,...,N; N is the size of
for solving MOOP. the swarm population; (cid:3) is a constriction factor which
Among the meta-heuristic techniques, until recently controls and constricts the velocity’s magnitude; w is
particle swarm optimization (PSO) was applied only to the inertial weight, which is often used as a parameter
single objective optimization tasks. However, the high to control exploration and exploitation in the search
speed of convergence of the PSO algorithm attracted space; c and c are positive constant parameters called
1 2
researcherstodevelopmulti-objective optimization algo- acceleration coefficients; r and r are random numbers,
1 2
rithms using PSO (Kennedy and Eberhart, 2001). Also, uniformlydistributedin[0,1];t isthetimestepusually
the PSO seems to have some advantages in terms of set as 1 and n is iteration number.
the better exploration and exploitation provided by local The successful application of PSO in many single
and global search capabilities of the algorithm. In the objectiveoptimizationproblemsreflectsitseffectiveness,
present study, a novel approach for multiple-objective anditseemstobeparticularlysuitableformultiobjective
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
MOPSOFORGENERATINGOPTIMALTRADE-OFFSINRESERVOIROPERATION 2899
optimization due to its efficiency in yielding better Elitist-mutated multi-objective particle swarm
quality solutions while requiring less computational time optimization
(Kennedy and Eberhart, 2001). The main difficulty in The main algorithm consists of initialization of pop-
extending PSO to multi-objective problems is to find the ulation, evaluation, and reiterating the search on swarm
best way of selecting the guides for each particle in the by combining PSO operators with Pareto-dominance cri-
swarm. The difficulty is noticeable, as there are no clear teria. In this process, the particles are first evaluated and
concepts of local and global bests that can be clearly checked for dominance relation among the swarm. The
identified,whendealingwithmanyobjectivesratherthan non-dominatedsolutionsfoundarestoredinanERP,and
asingleobjective.Recentlyafewproposalsonextensions areusedtoguidethesearchparticles.Itusesvariablesize
of PSO technique to multi-objective optimization have ERP, in order to improve the performance of the algo-
been reported (for example, Parsopoulos and Vrahatis, rithm to save computational time during optimization.
2002; Hu etal., 2003; Li, 2003; Coello et al., 2004). If the size of ERP exceeds the restricted limit, then it
In this paper, an efficient method is presented for is reduced by using the crowded comparison operator,
PSO to solve MOOPs. The approach uses Pareto dom- which gives the density measure of the existing parti-
inance criteria for selecting non-dominated solutions; cles in the function space. Also, an efficient EM strategy
an external repository (ERP) for storing best solutions is employed for maintaining diversity in the population
found (elitism); crowding distance operator for creating and for exploring the search space. The combination of
effective selection pressure among the swarm to reach these operators helps the algorithm to effectively prop-
true Pareto optimal fronts; and incorporates an effective agate the search towards true Pareto optimal fronts in
elitist-mutation (EM) strategy for effective exploration further generations.
of the search space. The proposed elitist-mutated multi-
objective particle swarm optimization (EM-MOPSO) EM-MOPSOalgorithm
algorithm isdiscussed indetail inthe following sections. ThedevelopedEM-MOPSOalgorithmcanbesumma-
rized in the following steps.
MULTI-OBJECTIVE PARTICLE SWARM Step 1.Initializepopulation.Setiterationcountert D0.
OPTIMIZATION 1.The current position of the i-th particle X is
i
initialized with random real numbers within
Brief concepts of multi-objective optimization are pre-
the specified decision variable range; each
sentedfirstandthentheproposedalgorithmisexplained.
particle velocity vector V is initialized with
i
uniformlydistributedrandomnumberin[0,1].
Multi-objective optimization and Pareto optimality
2.Evaluate each particle in the population. The
A general MOOP can be defined as: minimize a
personal best position P, is set to X.
i i
function f(cid:1)x(cid:2), subject to p inequality and q equality
Step 2.Identify particles that give non-dominated solu-
constraints.
tionsinthecurrentpopulationandstorethemin
an ERP.
min. f(cid:1)x(cid:2)Dff (cid:1)x(cid:2)f (cid:1)x(cid:2)...f (cid:1)x(cid:2)gT
1 2 m Step 3.t DtC1.
x 2D (cid:1)3(cid:2) Step 4.Repeat the loop (step through PSO operators):
1.Select randomly a global best P for the i-th
g
where x 2Rn, fi :Rn !R and particle from the ERP.
(cid:1) 2.Calculate the new velocity V, based on
x 2Rn : l (cid:2)x (cid:2)u, 8iD1,....,n i
i i
Equation (1),andthenewX byEquation(2).
DD g (cid:1)x(cid:2)½0, 8jD1,....,p (cid:1)4(cid:2) i
j
3.Repeat the loop for all the particles.
h (cid:1)x(cid:2)D0, 8k D1,....,q
k
Step 5.Evaluate each particle in the population.
where m is number of objectives; D is feasible search Step 6.Perform the Pareto dominance check for all the
space; x Dfx1x2...xngT is the set of n-dimensional particles:ifthecurrentlocalbestPiisdominated
decision variables (continuous, discrete or integer); R is by the new solution, then P is replaced by the
i
the set of real numbers; Rn is n-dimensional hyper-plane new solution.
or space; and l and u are lower and upper limits of i-th Step 7.Set ERP to a temporary repository, TempERP
i i
decision variable. and empty ERP.
TheMOOPshouldsimultaneously optimizethevector Step 8.Identify particles that give non-dominated solu-
function and produce Pareto optimal solutions. Pareto tions in current iteration and add them to
frontisasetofParetooptimal(non-dominated)solutions, TempERP.
beingconsideredoptimal,ifnoobjectivecanbeimproved Step 9.Find the non-dominated solutions in TempERP
without sacrificing at least one other objective. On the and store them in ERP. The size of ERP is
other hand, a solution xŁ is referred to as dominated by restricted to the desired set of non-dominated
another solution x, if and only if, x is equally good or solutions; if it exceeds, use the crowding dis-
betterthanxŁwithrespecttoallobjectives(Haimesetal., tance operator to select the desired ones. Empty
1990). the TempERP.
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
2900 M.J.REDDYANDD.NAGESHKUMAR
Step 10.Perform EM operation on specified number of 1.Getthenumberofnon-dominatedsolutionsintheERP
particles.
Step 11.Checkforterminationcriteria;ifthetermination lDjERPj
criterion is not satisfied, then go to step 3;
2.Initialize distance.
otherwiseoutputthenon-dominatedsolution set
from ERP.
For iD1 to l
Themainoperatorsusedinthisalgorithmareexplained ERP[i].distD0
below.
3.Compute the crowding distance of each solution.
For each objective m,
Variable size externalrepository
Sort using objective value.
The selection of the global best guide of the parti-
cle swarm is a crucial step in a multi-objective PSO ERPDsort (cid:1)ERP,m(cid:2)
algorithm. It affects both the convergence capability of
the algorithm as well as maintaining a good spread of Set the boundary points to a large value so that they
non-dominated solutions. As ERP stores non-dominated are always selected.
solutions found in the previous iteration, any one of the
solutions can be used as global guide. But we want to ERP[1].distDERP[l].distD1
ensure that the particles in the population move towards
For iD2 to (l-1)
the sparse regions of the non-dominated solutions and
speed up the convergence towards the true Pareto opti-
mal region. To performthese tasks, the global best guide
ERP[i].distDERP[i].distC(cid:1)ERP[iC1].m
of the particles is selected from a restricted variable size
ERP.ThisrestrictiononERPisdoneusingthecrowding (cid:1)ERP[i(cid:1)1].m(cid:2)/(cid:1)fmax(cid:1)fmin(cid:2)
m m
distance operator. This operator ensures that those non-
dominated solutions with the highest crowding distance Elitist-mutation operator
values are always preferred to be in the ERP. The other
To maintain diversity in the population and to explore
advantage of this variable size ERP is that it saves con-
the search space, a novel strategic mechanism called
siderablecomputationaltime duringoptimization. Asthe
EM is incorporated into the algorithm. This acts on a
ERP size increases, the computing requirement becomes
predefined number of particles. In the initial phase of
greaterforsortingandcrowdingvaluecalculations.Thus
thismechanism,ittriestoreplacetheinfeasiblesolutions
for effective exploration of the function space, the size
with the mutated least crowded particles of ERP and at
is initially set to 10% of maximum ERP, and then the
the later phase, it tries to exploit the searchspace around
value is increased in a stepwise manner, so that at the
the sparsely populated particles in ERP along the Pareto
start of 90% of maximum iterations, it reaches the max-
fronts. This is a special strategic mechanism, which
imum size of ERP. Selecting different guides for each
enhances the performance of MOPSO while extending
particle from a restricted repository allows the particles
from traditional PSO algorithm. Thus the EM operator
better exploration of the true Pareto optimal region. This
helpstouniformlydistributethenon-dominatedsolutions
kind of selection is novel and it effectively improves the
along the true Pareto optimal front. The pseudo-code of
performance of the algorithm.
the elitist mutation mechanism is given below.
Crowdingdistance assignment operator
1.Randomly select one of the objectives from m objec-
This operator is adopted from Deb et al. (2002). The tives. Sort the fitness function of particles in descend-
crowding distance value of a solution provides an esti- ing order and get the index number descending order
mate of the density of solutions surrounding that solu- sorted particles (DSP) for the respective particles.
tion. Crowding distance is calculated by first sorting the 2.Use crowdingdistance assignment operator andcalcu-
set of solutions in ascending objective function values. latethedensityofsolutionsintheERPandsortthemin
The crowding distance value of a particular solution is descending order of crowding value. Randomly select
the average distance of its two neighboring solutions. one of the least crowded solutions from the top 10%
The boundary solutions that have the lowest and high- of ERP as guide (g).
est objective function values are given infinite crowding 3.Perform EM on a predefined number of particles
distance values, so that they are always selected. This (NMmax).
process is done for each objective function. The final
crowding distance value of a solution is computed by Let R —be the size of repository; p - probability
p em
addingalltheindividualcrowdingdistancevaluesineach of elitist mutation; S - mutation scale used to preserve
m
objective function. For sorting, an efficient quick sorting diversity; rand - uniformly distributed random number
procedureisused.Thepseudo-codeofcrowdingdistance U(0,1); intRnd(a, b) - uniformly distributed integer
computation is given below. random number in the interval [a, b]; randn - Gaussian
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
MOPSOFORGENERATINGOPTIMALTRADE-OFFSINRESERVOIROPERATION 2901
random number N(0,1); and VR[i]- range of decision EM-MOPSO APPLICATION AND PERFORMANCE
variable i. EVALUATION
For iD1 to NM In order to demonstrate the efficiency of the proposed
max
EM-MOPSO, it is first tested for a few standard test
lDDSP[i]
problems taken from the MOEAs literature and its
gDintRnd(cid:1)1,0Ð1ðR (cid:2) performance is evaluated with results of NSGA-II.
p
For dD1 to dim
Testproblems
if (cid:1)rand<pem(cid:2) The four test functions considered to test the perfor-
X[l] [d]DERP[g][d]CS Ł VR[d]Ł randn mance of the proposed algorithm are given in Table I
m
(Deb, 2001). The first test problem (BNH) is a MOOP,
else with two objectives subject to two constraints. The sec-
ond test problem (KITA) is a MOOP, with maximiza-
X[l] [d]DERP[g][d]
tion of two objectives subject to three constraints. This
End For problem has non-convexity in its Pareto optimal region.
The third test problem (CONSTR) is a MOOP, with
End For
two objectives subject to two constraints. This problem
If the mutated value exceeds the bounds, then it is has the difficulty that a part of the unconstrained Pareto
limited to the upper or lower bound. The velocity vector optimal region is not feasible. Thus the resulting con-
of the particle remains unchanged during this EM step. strained Pareto optimal region is a concatenation of the
first constraint boundary and some part of unconstrained
Constraint handling Pareto optimal region. The fourth test problem (SRN) is
a MOOP, with two objectives subject to two constraints.
In order to handle the constrained optimization prob-
Here the constrained Pareto optimal set is a subset of
lems, this study adopts the constraint handling mecha-
the unconstrained Pareto-optimal set, which gives diffi-
nism proposed by Deb etal. (2002). This is a simple,
culty in finding the true Pareto optimal region for the
but very effective procedure reported in the literature. In
algorithm.
this approach, a solution i is said to be a constrained-
dominate solution j if any of the following conditions
Sensitivity ofEM-MOPSO parameters
hold good:
ThesensitivityanalysisofthePSOmodelisperformed
1.Solution i is feasible and solution j is not. with different combinations of each parameter. In this
2.Bothsolutionsiandjareinfeasible,butsolutionihas analysis, it is observed that by considering the proper
a smaller overall constraint violation. value for the constriction coefficient, the inertial weight
3.Both solutions i and j are feasible and solution i does not have much influence on the final result of the
dominates solution j. model (Nagesh Kumar and Janga Reddy, 2006). So in
this study the inertial weight (w) is fixed as 1. Also, it
By using all the above steps, the EM-MOPSO is foundthat the value of constriction coefficient(cid:3) equal
approachiscodedinuserfriendlymathematicalsoftware to 0Ð9 yields better results for the given model. After a
package MATLAB 6Ð5 and is run on PC/WindowsXP/ number of trials, it was found that constant parameters
256MB RAM/2GZ computer. The applicability and effi- c D1Ð0 and social parameter c D0Ð5 resulted in better
1 2
ciency of the proposed approach is demonstrated in the quality solutions. The same values are used for all the
following sections. test problems.
TableI.Test problemsusedinthestudy
Problem Variablebounds Objectivefunctions Constraints
BNH x 2[0,5] Minimize g (cid:1)x(cid:2)D(cid:1)x (cid:1)5(cid:2)2Cx2 (cid:2)25
1 1 1 2
x 2[0,3] f (x)D4x2C4x2 g (cid:1)x(cid:2)D(cid:1)x (cid:1)8(cid:2)2C(cid:1)x C3(cid:2)2½7Ð7
2 1 1 2 2 1 2
f (x)D(cid:1)x (cid:1)5(cid:2)2C(cid:1)x (cid:1)5(cid:2)2
2 1 2
KITA x 2[0,7] Maximize g (cid:1)x(cid:2)Dx /6Cx (cid:1)6Ð5(cid:2)0
i 1 1 2
iD1,...,3 f (cid:1)x(cid:2)D(cid:1)x2Cx g (cid:1)x(cid:2)D0Ð5x Cx (cid:1)7Ð5(cid:2)0
1 1 2 2 1 2
f (cid:1)x(cid:2)D0Ð5x Cx C1 g (cid:1)x(cid:2)D5x Cx (cid:1)30(cid:2)0
2 1 2 3 1 2
CONSTR x 2[0Ð1,1Ð0] Minimize g (cid:1)x(cid:2)Dx C9x ½6
1 1 2 1
x 2[0,5] f (x)Dx g (cid:1)x(cid:2)D(cid:1)x C9x ½1
2 1 1 2 2 1
f (x)D(cid:1)1Cx (cid:2)/x
2 2 1
SRN x 2[(cid:1)20,20] Minimize g (cid:1)x(cid:2)Dx2Cx2(cid:2)225
i 1 1 2
iD1,2 f (x)D(cid:1)x (cid:1)2(cid:2)2C(cid:1)x (cid:1)1(cid:2)2C2 g (cid:1)x(cid:2)Dx (cid:1)3x (cid:2)(cid:1)10
1 1 2 2 1 2
f (x)D9x (cid:1)(cid:1)x (cid:1)1(cid:2)2
2 1 2
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
2902 M.J.REDDYANDD.NAGESHKUMAR
Size of elitist-mutated particles (NM ). The number TableII.Resulting statistics by EM-MOPSO and NSGA-II for
max
ofparticlestobeelitistmutatedisselectedafterensuring test problems, considered in the study. In SC(A, B), A is
EM-MOPSOandBisNSGA-II.Boldnumbersindicatethebest
that the population does not lose control on search at the
performingalgorithm
cost of exploring for better non-dominated solutions. To
experiment with the size of elitist mutated particles, the Test Statistic Performancemetric
number of particles is varied as 5, 10, 15, 20, 25 and 30 problem
Set coverage Spacing
for the problems having a maximum population of 100.
metric(SC) metric(SP)
The best results are found for NM D20 and is kept
max
constantforallthetestproblemsconsideredinthisstudy. SC(A,B) SC(B,A) EM- NSGA-II
MOPSO
Probabilityofelitistmutation(p ).p isvariedfrom
em em Best 0Ð1400 0Ð1200 0Ð6357 0Ð6408
0 to 0Ð5 for sensitivity analysis. It is found that best
Worst 0Ð0900 0Ð0500 0Ð7559 0Ð8928
performance occurs at pem D0Ð2, and is kept constant BNH Mean 0.1111 0Ð0877 0.6941 0Ð7756
for all the test problems considered in the study. Variance 0Ð0003 0Ð0005 0Ð0015 0Ð0053
SD 0Ð0176 0Ð0233 0Ð0385 0Ð0727
Best 0Ð2900 0Ð2400 0Ð0374 0Ð0496
Elitist-mutation operator mutation scale (S ). The
m Worst 0Ð1200 0Ð1300 0Ð4254 0Ð5117
other parameter used in EM operation is mutation scale KITA Mean 0.2400 0Ð1811 0.1359 0Ð1464
(Sm). After various experiments, it is found that a value Variance 0Ð0015 0Ð0013 0Ð0196 0Ð0227
ofS intherangeof0Ð2to0Ð01givesgoodperformance. SD 0Ð0394 0Ð0355 0Ð1401 0Ð1507
m
Best 0Ð1600 0Ð1700 0Ð0374 0Ð0372
This is selected after ensuring that at the initial stage it
Worst 0Ð0700 0Ð1000 0Ð0431 0Ð0487
did not deteriorate the search while exploring the search
CONSTR Mean 0Ð1181 0.1344 0.0406 0Ð0437
space or stagnate the search at the end of iterations. Variance 0Ð0008 0Ð0004 0Ð0000 0Ð0000
SD 0Ð0281 0Ð0201 0Ð0017 0Ð0041
Simulation results Best 0Ð1400 0Ð1400 1Ð0768 1Ð3402
Worst 0Ð0400 0Ð0400 1Ð3929 1Ð7073
To run the EM-MOPSO algorithm, the following SRN Mean 0.0978 0Ð0944 1.2439 1Ð5860
parameters are used: size of populationD100; constant Variance 0Ð0014 0Ð0009 0Ð0114 0Ð0179
parametersc D1Ð0andc D0Ð5;inertialweightwD1; SD 0Ð0370 0Ð0305 0Ð1055 0Ð1337
1 2
constriction coefficient (cid:3) D0Ð9; size of ERPD100; the
size of elitist-mutated particles is set to 20, the value
Pareto optimal fronts than NSGA-II. With regard to the
of p was set to 0Ð2; and the value of S decreases
em m spacing metric (SP), as compared to NSGA-II, EM-
from 0Ð2 to 0Ð01 over the iterations. To run the NSGA-II
MOPSOgivessmallerSPvaluesforallthetestproblems
model, the initial population was set to 100, crossover
considered in the study. The smaller SP indicates that
probability to 0Ð9, and mutation probability to 1/n
the algorithm gives better distribution of solutions. Thus
(nisthenumberofrealvariables).Thedistributionindex
EM-MOPSO maintains the best distribution of solutions
values for real-coded crossover and mutation operators
for all the test problems. For illustration purposes, a
are set to 20 and 100 respectively (Deb etal., 2002).
sample result for each of the test problems considered is
Maximum number of iterations in a run is set to 250 for
shown in the plots in Figure 1. Thus the results obtained
both the algorithms. The same parameter settings were
clearlyshowthattheproposedmethoddoesnothaveany
used for all the problems. To evaluate the performance
difficulty in achieving a good spread of Pareto optimal
of the proposed EM-MOPSO algorithm, this study uses
solutions for constrained multi-objective optimization.
twoperformancemeasures,setcoveragemetric(SC)and
spacing metric (SP) (Deb, 2001). The details of these
performance metrics are presented in APPENDIX- A.
CASE STUDY
Table II shows the best, worst, mean, variance and
standard deviation (SD) values of the two performance To demonstrate the efficacy of the proposed approach,
metrics (SC and SP) obtained from 10 independent the Bhadra reservoir system in India is taken up as a
runs using EM-MOPSO and NSGA-II. The set coverage case study for developing optimal reservoir operation
metrics SC (A, B) and SC (B, A) give a measure of how policy. Figure 2 shows the location map of the Bhadra
many solutions of A are covered by B and vice versa. reservoir system. The Bhadra dam is located at latitude
Here, the value SC (cid:1)A, B(cid:2)D1 means that all solutions 13°420 N and longitude 75°3802000 E. The Bhadra reser-
in B are weakly dominated by A, while SC (cid:1)A, B(cid:2)D0 voirisamulti-purposeprojectprovidingfacilitiesforirri-
represents the situation when none of the solutions in gation,hydropowergenerationandmeetingwaterquality
B are weakly dominated by A. It can be seen that requirements downstream. The schematic diagram of the
with respect to the SC metric, the average performance reservoir system is given in Figure 3.
of EM-MOPSO is the best for test functions BNH, Mostoftheinflowsintothereservoirarereceiveddur-
KITA and SRN, whereas NSGA-II performs best for ing the monsoon season of 4 months. But the demands
the CONSTR problem. This metric shows the efficiency are distributed throughout the year. The reservoir pro-
of EM-MOPSO in achieving better convergence to true videswaterforirrigationof6367 haand87512 haunder
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
MOPSOFORGENERATINGOPTIMALTRADE-OFFSINRESERVOIROPERATION 2903
(a) BNH (b) KITA
50 8.6
EM-MOPSO EM-MOPSO
40 NSGA-II 8.4 NSGA-II
8.2
30
2 2 8.0
f f
20
7.8
10
7.6
0 7.4
0 30 60 90 120 150 −4.0 −2.0 0.0 2.0 4.0 6.0 8.0
f1 f1
(c) CONSTR (d) SRN
10
EM-MOPSO 0 EM-MOPSO
8 NSGA-II NSGA-II
−50
6
2 2 −100
f f
4 −150
2 −200
0 −250
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 0 50 100 150 200 250
f1 f1
Figure1.Non-dominatedsolutionsobtainedusingEM-MOPSOandNSGA-IIforthetestproblems(a)BNH,(b)KITA,(c)CONSTRand(d)SRN
left and right bank irrigation canals respectively. The Modelformulation
irrigated area spread over the districts of Chitradurga,
The objectives of the reservoir operation model
Shimoga, Chikmagalur, and Bellary in Karnataka state,
are: minimizing the irrigation deficits, maximizing the
comprises predominantly of red loamy soil, except in
hydropower generation and maximizing the satisfaction
someportionsoftherightbankcanalarea,whichconsists
level of water quality. These are conflicting and/or com-
of black cotton soil. Major crops grown in the com-
petitiveobjectives.Forexample,conflictmayariseduring
mand area are paddy, sugarcane, permanent garden, and
dry periods: for minimization of irrigation deficits, more
semidrycrops.Alsounderthisprojecttherearethreesets
wateristobereleasedtosatisfyirrigationdemands;while
of turbines, one set each on the left bank canal and the formaximizationofhydropowerproduction,higherlevel
right bank canal and the other set at the river bed level of storage in the reservoir is required to produce more
of the dam, generating hydropower. The operating head hydropower energy; and for maximizing the satisfaction
aboveriverbed,rangesfrom38Ð56 mto54Ð41 mfor the levelof water quality, steady releaseof water is required
right bank turbine (PH2), and from 36Ð88 m to 56Ð69 m to meet river water quality requirements downstream.
for left bank turbine (PH1) and bed turbines (PH3). The Thus, solving the allocation problems of this reservoir
mean tail water levels of right bank, left bank, and bed system is interesting from the multi-objective perspec-
turbines are at 32Ð736 m, 12Ð802 m and 6Ð706 m above tive. In order to simplify the water quality objective,
the bed level respectively. It can be noted that the water in this study it is assumed that if we discharge a pre-
released to left bank and right bank canals goes through specified amount of water into the downstream river, the
turbines only when the water is within the limits of tur- river water quality can be maintained. The water quality
bineoperatingrange,otherwiseitwillbereleaseddirectly demandsarechosenaftercarefullystudyingthehistorical
forirrigation.Waterqualityisalsoamajorconcerntothe data and previous studies on river water quality mainte-
reservoir authorities due to continuous development of nance.Tomaintain evendistribution of irrigationdeficits
industries in the downstream region. So the water qual- if any, the irrigation objective is taken as squared devi-
ityobjectiverequirescertainminimum waterlevelstobe ation of demand to release. The competing objectives of
maintained in the river downstream. the system are expressed as follows:
Salient features of the reservoir are given in Table III. Minimize sum of squared deviations for irrigation
Datapertainingtomonthlyinflowsandotherdetailswere annually:
collected from water resources development organiza-
(cid:2)12 (cid:2)12
tion (WRDO), Bangalore covering a period of 69 years
SQDVD (cid:1)D (cid:1)IR (cid:2)2C (cid:1)D (cid:1)IR (cid:2)2 (cid:1)5(cid:2)
(from 1930–1931 to 1998–1999). The monthly crop 1,t 1,t 2,t 2,t
tD1 tD1
waterrequirementswerecalculatedusingFoodandAgri-
cultural Organization (FAO) Penman-Monteith method where SQDV is the sum of squared deviations of
(Allen etal., 1998). irrigation demands from releases. D and D are the
1,t 2,t
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
2904 M.J.REDDYANDD.NAGESHKUMAR
Figure2.LocationmapoftheBhadraprojectcommandarea
Bhadra river TableIII.Salientfeaturesof Bhadrareservoirsystem
Description Quantity
Right turbine Left turbine
capacity = 13,200 kW capacity = 2,000 kW
Reservoir Grossstoragecapacity 2025 Mm3
PH2 PH1
Livestoragecapacity 1784 Mm3
Rcaigphatc ityb=a n7k1 mca3/nsal PcHa3paBceityd= t2u4rb,0in0e0 kWLcaepfta bciatyn k= 1ca0n mal3/s DAveeardagsteorAangneucaalpiancflitoyw 2284415MMmm33
Leftbank canal capacity 10 m3/s
Irrigated area = 87, 512 ha Irrigated area = 6, 367 ha
Annual demand = 1, 911 Mm3 Annual demand = 227 Mm3 Rightbankcanal capacity 71 m3/s
Leftbank turbinecapacity (PH1) 2000 Kw
Rightbankturbinecapacity (PH2) 13200 Kw
Riverbedturbinecapacity(PH3) 24000 kW
Figure3.SchematicdiagramoftheBhadrareservoirproject
Maximize annual hydropower production:
irrigation demands for the left bank canal and right bank
canal command areas respectively in period t in Mm3;
(cid:2)12
IR1,t and IR2,t are the irrigation releasesinto the left and PD p(cid:1)R H CR H CR H (cid:2) (cid:1)6(cid:2)
1,t 1,t 2,t 2,t 3,t 3,t
right bank canals respectively in period t in Mm3.
tD1
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
MOPSOFORGENERATINGOPTIMALTRADE-OFFSINRESERVOIROPERATION 2905
where P is the total energy produced in M kWh; p is Irrigation demands:
power production coefficient; R ,R and R are the
1,t 2,t 3,t
D (cid:2)IR (cid:2)D 8tD1,2,...,12(cid:1)16(cid:2)
releases to left bank, right bank and river bed turbines 1min,t 1,t 1max,t
respectively in period t in Mm3. H1,t,H2,t,H3,t are the D2min,t (cid:2)IR2,t (cid:2)D2max,t 8tD1,2,...,12(cid:1)17(cid:2)
net heads available to the left bank, right bank and bed
turbines respectively in meters during period t (here, where, D1min,t and D1max,t are minimum and maxi-
head is a nonlinear function of initial and final reservoir mum irrigation demands for left bank canal respectively;
storage). D2min,t and D2max,t are minimum and maximum irriga-
Maximize satisfaction level of river water quality: tion demands for right bank canal respectively in time
period t.
WQD min (cid:1)(cid:5)t(cid:2) (cid:1)7(cid:2) WaterQuality Requirements:
8tD1,2,...,12
R ½MDT 8tD1,2,...,12 (cid:1)18(cid:2)
where, (cid:5) is satisfaction level of water quality in period 3,t t
t
t and is given by,
where, MDT Dminimum release to meet downstream
 t
0 if R3,t (cid:2)QDmin,t water quality requirement in Mm3.
(cid:1)R (cid:1)QD (cid:2) Itcanbenotedthat,inthisstudy,underfavorablerange
(cid:5) D 3,t min,t if QD (cid:2)R (cid:2)QD
t (cid:1)QDmax,t(cid:1)QDmin,t(cid:2) min,t 3,t max,t of reservoir storage for power production, the releases
1 if R ½QD madethroughpowerturbinesalsoservetomeetirrigation
3,t max,t
(cid:1)8(cid:2)
demands of the left bank and right bank canals during
where, QD and QD are the minimum and maxi-
min,t max,t irrigation requirement periods. However, if the reservoir
mum water demands to maintain water-quality in period
storage is not within the limits of the turbine power
t in Mm3 for the river downstream of the dam.
production range, then the water is releasedonly to meet
The optimization is subject to the following con-
irrigationdemandsthroughsub-waystoirrigationcanals.
straints:
Since the water released for irrigation is restricted to its
Storage continuity:
demands, any excess water is not a penalizing problem
in the first objective function as given in Equation (5).
S DS CI (cid:1)(cid:1)R CR CR CE CO(cid:2)
tC1 t t 1,t 2,t 3,t t t
8tD1,2,...,12 (cid:1)9(cid:2)
RESERVOIR OPERATION MODEL APPLICATION
where S DActive reservoir storage at the beginning AND RESULTS
t
of period t in Mm3; I Dinflow into the reservoir during
t To apply the EM-MOPSO for reservoir operation model,
period t in Mm3; E Dthe evaporation losses during
t thefollowingparametersareselected.Theinitialpopula-
periodtinMm3 (here,E isanonlinearfunctionofinitial
t tionoftheEM-MOPSOissetto200;c andc aresetto
1 2
and final storages of period t); O Doverflow from the
t 1Ð0 and 0Ð5 respectively; w is set to 1Ð0; and (cid:3) is set to
reservoir in period t in Mm3;
0Ð9; the number of non-dominated solutions to be found
Storage limits:
is set to 200. For the EM step, the size of elitist-mutated
particlesissetto30,thevalueofp wassetto0Ð2;and
S (cid:2)S (cid:2)S 8tD1,2,...,12 (cid:1)10(cid:2) em
min t max the value of S decreases from 0Ð2 to 0Ð01 over the iter-
m
ations. Then EM-MOPSO is run for 500 iteration steps.
where S and S are the minimum and maximum
min max
ToruntheNSGA-IImodel,theinitialpopulationwasset
active storages of the reservoir in Mm3.
to 200, crossover probability to 0Ð9, and mutation prob-
Maximum powerproduction limits:
ability to 1/n (n is the number of real variables). The
pR H (cid:2)E 8tD1,2,...,12 (cid:1)11(cid:2) distribution index values for real-coded crossover and
1,t 1,t 1,max
mutation operators are set to 20 and 100 respectively.
pR H (cid:2)E 8tD1,2,...,12 (cid:1)12(cid:2)
2,t 2,t 2,max NSGA-II is also run for 500 generations. The average
pR H (cid:2)E 8tD1,2,...,12 (cid:1)13(cid:2) monthly inflow into the reservoir computed for each cal-
3,t 3,t 3,max
endarmonthoveraperiodof69 years(from1930–1931
where, E ,E , and E are the maximum to 1998–1999) is used as inflow data to implement the
1,max 2,max 3,max
amounts of power in M kWh, that can be produced above model.
(turbinecapacity)bytheleft,rightandbedlevelturbines
respectively. Two-objectivemodel
Canalcapacity limits: ToshowtheeffectivenessoftheproposedMOPSO,for
solvingthereservoiroperationproblem,firstitisapplied
IR (cid:2)C 8tD1,2,...,12 (cid:1)14(cid:2)
1,t 1,max to irrigation and hydropower as two objectives of the
IR (cid:2)C 8tD1,2,...,12 (cid:1)15(cid:2) model. The reservoir operation model consists of mini-
2,t 2,max
mizing irrigation deficits (Equation (5)) and maximizing
where,C andC arethemaximumcanalcarrying the hydropower (Equation (6)) subject to satisfying the
1,max 2,max
capacities of the left and right bank canals respectively. constraints from Equations (9) to (18).
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp
2906 M.J.REDDYANDD.NAGESHKUMAR
To check the performance, 10 independent runs were 235
carried out for the two-objective reservoir operation
225
model using both the algorithms. Table IV shows the
215
resulting statistics for both the EM-MOPSO and NSGA-
205
II models. It can be observed that with respect to set
195
coveragemetric,the averagevalueof SC(A,B)ishigher
2
f 185
thantheSC(B,A)value(hereAisEM-MOPSOandBis
NSGA-II). The metric SC(A, B) refers to the percentage 175
of solutions in B that are weakly dominated by solutions 165 EM-MOPSO
ofA.Thusinthiscase,EM-MOPSOisperformingbetter 155 NSGA-II
thantheNSGA-II.Similarly,forthespacingmetric,from 145
Table IV it can be observed that the mean value of SP 0 25000 50000 75000 100000 125000
metric for EM-MOPSO is lower than that for NSGA-II. f1
This indicates that best distribution of Paretosolutions is Figure4.Non-dominated solutions obtained using EM-MOPSO and
obtained in EM-MOPSO. For demonstration purposes a NSGA-IIfortwo-objectivereservoiroperationmodel.Thetwoobjectives
sample result corresponding to medianvalue of SC(A,B) aref1-squareddeviationofirrigation(cid:1)Mm3(cid:2)2,f2-hydropowerproduction,
MkWh
is shown in Figure 4. It can be seen that EM-MOPSO
is able to generate a set of well-distributed Pareto-
optimal solutions. In this case NSGA-II fails to yield
the extreme solutions, whereas our proposed method
1
is able to provide extreme solutions comfortably (the
minimumsquareddeviationofirrigationas0(cid:1)Mm3(cid:2)2and 0.8
maximum hydropower production as 227Ð845 MkWh).
0.6
Thus again these results demonstrate that the proposed 3
f
method is very useful for real life systems application. 0.4
0.2
Three-objectivemodel
0
Thereservoiroperationmodelforthreeobjectivescon- 250
sists of minimizing the irrigation deficits (Equation (5)), 200 2
1.5
maximizingthehydropower(Equation (6))andmaximiz-
150 1
ing the satisfaction level of water quality (Equations (7) 0.5 x 105
and (8)) subject to satisfying the constraints from f2 100 0 f1
Equations (9) to (18). Figure 5 shows the results of 200
Figure5.Non-dominatedsolutionsobtainedforthree-objectivereservoir
non-dominated solutions obtained using EM-MOPSO operationmodelusingEM-MOPSO.f1-sumofsquareddeficitsofirriga-
after 500 iterations. There are a number of alternatives tionreleases(cid:1)Mm3(cid:2)2;f2-hydropowerproduction,MkWh;f3-satisfaction
levelofwaterquality
that can be chosen at various satisfaction levels of the
multipleobjectives.Dependingonthecircumstancespre-
Decisionmaking
vailing under the reservoir system and by analyzing the
tradeoff between the multiple objectives, the reservoir In anyapplication, for final decision making, the deci-
operator can make an appropriate decision. The details sion maker might be interested in minimum possible
of decision making for application are presented in the numberofwellrepresentativesolutionsforfurtheranaly-
following section. sis. Soit is important that after obtaining many solutions
which are true Pareto Optimal with uniform spread and
wide coverage, we need to reduce the large set of solu-
TableIV.Resulting statistics by EM-MOPSO and NSGA-II for
tions to a few representative solutions. In order to do
the two-objective reservoir operation model. In SC(A,B), A is
EM-MOPSOandBisNSGA-II.Boldnumbersindicatethebest that, various clustering algorithms are available. A sim-
performingalgorithm pleclustering algorithm, whichreducesthe largenumber
offinalParetosolutions (N)toafewrepresentativesolu-
Statistic Performancemetric
tions (N), is described here.
Set coveragemetric Spacing metric
(SC) (SP) Clustering technique
First each solution in ERP is considered to reside in
SC(A,B) SC(B,A) EM-MOPSO NSGA-II
a separate cluster. Thus initially there are N clusters.
Best 0Ð2000 0Ð0400 206Ð1277 246Ð2255 Thereafter the cluster distances between all pairs of
Worst 0Ð8788 0Ð6350 294Ð1276 787Ð7765 clusters are calculated. Then the two clusters with the
Mean 0.5582 0Ð2878 258.2752 504Ð3212
minimum cluster distance are combined together to form
Variance 0Ð0435 0Ð0317 1419Ð9885 32583Ð4308
one big cluster. The procedure is repeated by calculating
SD 0Ð2085 0Ð1780 37Ð6827 180Ð5088
the cluster distances for all the pairs of clusters obtained
Copyright2007JohnWiley&Sons,Ltd. Hydrol.Process.21,2897–2909(2007)
DOI:10.1002/hyp