Multi-objective global optimization for hydrologic models

The development of automated (computer-based) calibration methods has focused mainly on the selection of a single-objective measure of the distance between the model-simulated output and the data and the selection of an automatic optimization algorithm to search for the parameter values which minimize that distance. However, practical experience with model calibration suggests that no single-objective function is adequate to measure the ways in which the model fails to match the important characteristics of the observed data. Given that some of the latest hydrologic models simulate several of the watershed output fluxes (e.g. water, energy, chemical constituents, etc.), there is a need for effective and efficient multi-objective calibration procedures capable of exploiting all of the useful information about the physical system contained in the measurement data time series. The MOCOM-UA algorithm, an effective and efficient methodology for solving the multiple-objective global optimization problem, is presented in this paper. The method is an extension of the successful SCE-UA single-objective global optimization algorithm. The features and capabilities of MOCOM-UA are illustrated by means of a simple hydrologic model calibration study.


Introduction and scope
To calibrate a hydrologic model, the hydrologist must specify values for its "parameters" in such a way that the model's behavior closely matches that of the real system it represents.In some cases, the appropriate values for a model parameter can be determined through direct measurements conducted on the real system.However, in a great many situations, the model parameters are conceptual representations of abstract watershed characteristics and must be determined through a trial-and-error process which adjusts the parameter values so that the model response matches the historical input-output data.
Because of the time-consuming nature of manual trial-and-error model calibration, there has been a 0022-1694/98/$19.00© 1998 Elsevier Science B.V. All rights reserved PI1 S0022-1694(97)  great deal of research into the development of automated (computer-based) calibration methods (see, e.g., Gupta and Sorooshian, 1994).These efforts have focused mainly on the selection of a singleobjective measure of the distance between the model-simulated output and the data and the selection of an automatic optimization algorithm to search for the parameter values which minimize that distance.Objective functions that have been shown to work well in practice include the Mean Squared-Error Estimator (MSE) and the Heteroscedastic Maximum Likelihood Estimator (HMLE) criterion (Sorooshian and Dracup, 1980).Research into optimization methods has led to the use of population-evolutionbased search strategies (e.g.Brazil and Krajewski, 1987;Brazil, 1988;Wang, 1991;Duan et al., 1992Duan et al., , 1993;;Sorooshian et al., 1993;among others).In this regard, the Shuffled Complex Evolution (SCE-UA) global optimization algorithm has been found to be consistent, effective, and efficient in locating the globally optimal model parameters of a hydrologic model with respect to a given objective function (Duan et al., 1992(Duan et al., , 1993;;Sorooshian et al., 1993;Luce and Cundy, 1994;Gan and Biftu, 1996;Tanakamaru, 1995;Tanakamaru and Burges, 1996, among others).
Practical experience with model calibration suggests that any single-objective function, no matter how carefully chosen, may not adequately measure the ways in which the model fails to match the important characteristics of the observed data.This is reflected in the fact that the U.S. National Weather Service typically uses as many as 10 different objective functions to measure the goodness-of-fit of the Sacramento Soil Moisture Accounting model (SAC-SMA) during a multi-stage semi-automated calibration procedure (Brazil, 1988).Further, many of the latest hydrologic models simulate several of the watershed output fluxes (e.g.water, energy, chemical constituents, etc.) for which measurement data are available, and all these data must be properly utilized to ensure proper model calibration (e.g.Beven and Kirkby, 1979;De Grosbois et al., 1988;Kuczera, 1982Kuczera, 1983;Woolhiser et al., 1990;Yan and Haan, 1991 a, b).In particular, watershed hydrochemical models may simulate several physical and chemical properties of the hydrograph, in addition to the rate of flow (Wolford and Bales, 1996), while land-surface hydrology models designed for coupling with General Circulation Models typically simulate several energy and water fluxes and state variables including latent heat, sensible heat, temperature, runoff, and soil moisture at various depths (Dickinson et al., 1993).Because many of these models employ distributed representations of the watershed, the state variables and output fluxes may also be simulated and measured at numerous locations.
Procedures for the proper calibration of complex hydrologic models must effectively and efficiently utilize the various measurement data time series that provide useful information about the physical system.In Gupta et al. (1997), the advantages of a multipleobjective representation of the model calibration problem were discussed, and this representation was shown to be applicable and desirable, even in the case of the classical single-output-flux hydrological model.In this paper, an effective and efficient methodology for solving the multiple-objective global optimization problem is presented, and its features and capabilities are demonstrated by means of a simple example.

Formulation
The multi-objective hydrologic model calibration problem can be stated as the optimization problem: where fl(O) .... fro(0) are the m non-commensurable objective functions to be simultaneously minimized with respect to the parameters 0 of the model.For example, each f/(0) might correspond to least-squares matching of a particular model output flux to the available data on that flux.The multi-objective formulation has its foundations in 18 th century economics and philosophical discussion of welfare theories and competitive equilibrium (Edgeworth, 1881;Arrow, 1968;Pareto, 1971).It was introduced into mathematics by Cantor (1895Cantor ( , 1897) ) and has since evolved into a well-established mathematical discipline (Kuhn and Tucker, 1951;Hurwicz, 1958) that has been applied to several fields, including game theory (e.g. von Neumann, 1928;Borel, 1953), production theory (e.g.Koopmans, 1951), engineering (e.g.Stadler, 1984a, b;Zadeh, 1963), and quite extensively in water resources (e.g.Marglin, 1967).A detailed historical retrospect on multiobjective optimization can be found in Stadler (1986).

Pareto optimality
The overriding characteristic of a multi-objective problem is that the solution will not, in general, be unique.In fact, it is common to have several solutions with the property that moving from one solution to another results in the improvement of one objective function while causing a deterioration in the value of at least one other objective function. .Because the individual objectives are non-commensurate (each function is its own entity), it is not possible to determine objectively among such solutions which one is the best, and such solutions are called Pareto, nondominated, non-inferior, or efficient solutions.The formulation presented above divides the feasible parameter space into two parts: the part consisting of all Pareto solutions 0p is called the Pareto set, and the part consisting of the remaining space (non-Pareto solutions 0d) is called the dominated set.By definition, any member 0p of the Pareto set has the properties that: 1. F(0p) is strictly less than F(Oa) for all members 0d not contained in the Pareto set; i.e. {f~(0p) <fi(0d), for i = 1 ..... m} 2. It is not possible to find a 0p within the Pareto set such that F(Op) is strictly less than F(0p).
Put simply, the first of these statements says that the feasible parameter space can be partitioned into "good" solutions (Pareto solutions) and "bad" solutions.The second says that, in the absence of additional information, it is not possible to distinguish any one of the "good" (Pareto) solutions as being objectively better than any of the other "good" solutions (i.e.there is no uniquely "best" solution).For a formal definition of Pareto optimality, please refer to Vincent and Grantham ( 1981).
While it may be relatively simple to state the multiobjective optimization problem as in Eq. (1), solving it to identify the Pareto set is not easy and has been the subject of much research.While explicit analytic solutions to the Pareto set can be derived for simple problems (Vincent and Grantham, 1981), this is typically not possible for many practical problems and, in general, only an approximation of the Pareto set by a finite number of solutions can be obtained.Several methods for the generation of solutions that sample the Pareto set have been proposed, two of the more popular ones being the weighing method (Zadeh, 1963;Goicoechea et al., 1982) and the econstraint method (Marglin, 1967;Cohon and Marks, 1975;Szidarovszky et al., 1986).Some methods attempt to generate only the Pareto subset of interest to a decision maker.Such methods are usually interactive; e.g. the tradeoff development method (Trade) of Goicoechea et al. (1976) and the Surrogate Worth Tradeoff (SWT) method of Haimes et al. (1975).Other methods require choosing a solution among several alternatives; e.g.compromise programming and goal programming.Descriptions of these methods can be found in Harboe (1992) and Goicoechea et al. (1982), among others.

Limitations of classical solution methods
The major strength of the classical "generating" methods that attempt to sample the entire Pareto set is that they have the potential to provide "global" solutions to the multi-objective optimization problem (just as in single-objective optimization, it is possible for local Pareto solutions to coexist with global Pareto solutions).However, a major weakness of these techniques arises in their strategy of converting the multiobjective optimization problem into a large number of single-objective optimization problems.For problems with a large number of objectives, the computational burden grows rapidly (at an exponential rate).The task of estimating the entire Pareto set becomes inefficient, cumbersome, and time consuming.Further, such methods find it difficult to generate solutions at the boundary of the Pareto set.In contrast, the strength of interactive methods is that they generate only a relatively small number of efficient solutions by moving from one feasible solution to the next using subjective articulation of preference among the different objectives.However, such methods may fail to explore the entire feasible space and therefore not arrive at a "global" solution (the final solution may depend on the initial point used in the procedure).

Discussion
The classical strategies used to solve the multiobjective optimization problem are, in principle, quite similar to the strategies used in single-objective optimization.The generating technique is similar to "global random search", and the interactive approach is similar to "local deterministic search".Therefore, it is not surprising that these multi-objective techniques have strengths and weaknesses which are similar to those of their single-objective counterparts (for a discussion of single-objective approaches, see Duan et al., 1992).Among single-objective approaches, a recent development is the strategy of "population evolution", which draws on the strengths of both random and deterministic methods to achieve, in a synergistic fashion, both global search ability and rapid convergence.Examples of population evolution strategies are the Genetic Algorithm (Holland, 1975) and the SCE-UA algorithm (Duan et al., 1992).
Because population-based strategies achieve their performance by the simultaneous evolution of numerous potential solutions towards the region of the global optimum, they are ideally suited to solving for the Pareto solution set of a multi-objective problem.In principle, multi-objective versions of such strategies would be able to search for an approximate representation of the Pareto set in a single optimization run.With this aim in mind, we have developed an algorithm entitled the Multi-Objective Complex Evolution (MOCOM-UA) method, based on extensions of the successful SCE-UA single-objective method previously developed by our research group.

The multi-objective complex evolution (MOCOM-UA) global optimization method
The MOCOM-UA method is a general-purpose global multi-objective optimization algorithm designed to be effective and efficient for a broad class of problems.The MOCOM-UA strategy combines the strengths "controlled random search" (Price, 1987) with the "competitive evolution" (Holland, 1975), Pareto ranking (Goldberg, 1989), and a newly developed strategy of multi-objective downhill simplex search.Although the complex shuffling strategy implemented in the SCE-UA method has some inherent strengths (improved global search and better computational efficiency) that are desirable in a population-based method, it has not yet been implemented into the algorithm.The essence of the MOCOM-UA method is described below; detailed descriptions and explanations appear in Yapo (1996).
The MOCOM-UA strategy is illustrated in Figs. 2  and 3. We begin with an initial sample of s points distributed randomly throughout the n-dimensional feasible parameter space U°(0) that represents the initial parameter uncertainty.In the absence of prior information about the location of the Pareto optimum, a uniform sampling distribution is used.For each point, the multi-objective vector F(O) is computed, and the population is ranked and sorted using a Pareto-ranking procedure suggested by Goldberg (1989) (see discussion in next subsection).Simplexes of n + 1 points are then selected from the sample according to a robust rank-based selection method (Whitley, 1989).Each simplex is evolved in an improvement direction using a multi-objective  Like the SCE-UA method, MOCOM-UA treats the global search as a process of natural evolution.The s sampled points constitute a population.Each member of the population is a potential parent with the ability to participate in a process of reproduction.A simplex selected from the population is like a pair of parents, except that a simplex contains more than two members (except in the trivial case that n = 1).To ensure that the evolution process is competitive, we require that the probability that "better" parents contribute to the generation of offspring is higher than that of "worse" parents.The use of a triangular probability distribution for parent selection ensures this competitiveness.The multi-objective simplex strategy is applied to each simplex to generate the off-spring.The information contained in the simplex is used to direct the evolution in an improvement direction.Each new offspring replaces the worst point of the current subcomplex.This ensures that each parent gets at least one chance to contribute to the reproduction process before being replaced or discarded; thus, none of the information in the sample is ignored.

Pareto ranking
Because in a multi-objective problem several objectives are to be considered simultaneously, ordered ranking of the population by conventional scalar sorting is not possible, and the concept of inferiority-superiority in multi-objective theory is used instead.The special sorting used in MOCOM-UA is called "Pareto ranking" (Goldberg, 1989).The procedure begins by identifying all the non-dominated individuals in the population and assigning them the rank "one".These individuals are set aside (temporarily removed from the population), and the nondominated individuals identified in the remaining population are assigned the rank "two".The process is repeated until every point has been assigned a rank.This sorting procedure essentially assigns equivalent rankings to all points that lie on the same Pareto frontier.In a population ofs individuals, the assigned ranks are 1,2,3 .... ,Rmax, where Rmax ~-~ S; the smallestranked individuals are closest to the Pareto optimal set, while the largest-ranked individuals are the furthest from it.
Points marked for evolution

Rank-based selection
Rank-based selection utilizes an individual's rank to determine its likelihood of producing offspring.In the SCE-UA algorithm, this is achieved by first sorting the individuals from smallest to largest function value and assigning to them a corresponding rank.Based on the individual ranks, a trapezoidal probability mass function is assigned so that points with lower rank have a greater probability of selection.In a similar manner, MOCOM-UA utilizes a rank-based selection procedure that favors the selection of points closest to the Pareto set (smaller Pareto ranking).Let ri be the computed rank of point i.Assignment of selection probability is done by fitting the following probability distribution to the population: s Pi=(Rmax-ri+ l)/ ~(Rmax-rj+ l) J The denominator in Eq. ( 2) ensures that {Pi} i=l to s is a probability mass function (pmf) (i.e.Z~iPi = 1).It has been shown that rank-based selection is a more robust scheme than function-based selection in multiobjective applications (Whitley, 1989).

Multi-objective complex evolution
Once all individuals have been ranked, evolution proceeds by selecting from the population a number of simplexes equal to the number of points NRm,x which have largest rank (rank equal to {Rmax}).Each simplex consists of n + 1 points, where one point is selected from the worst NRmax points, and the remaining n points are selected at random, with replacement, from the P-NRm~ remaining points according to the pmf in Eq. (2).Fig. 4 illustrates the procedure of simplex selection for a simple twodimensional problem with two objectives to be minimized (n = 2) and a population of seven points.
The numbers in parentheses next to each point indicate its Pareto rank; points A, B, and C are of rank one, points D and E are of rank two, and points F and G are of rank three (Rma x = 3).Because there are two points with worst rank, we form two simplexes for evolution.Each simplex contains one of the worstranked points and two (n) of the better-ranked points.
Although the simplexes shown in this illustration have no points in common, in practice, the simplexes can share better-ranked points.This strategy for simplex generation allows the worst points to be identified by a global ranking of the population (instead of a local ranking of the simplex) and determines global improvement directions.Because the simplexes are evolved independently of each other, the evolution procedure is well suited for parallel processing with its advantages of significantly reduced computer execution time.
When all simplexes have been evolved exactly once, the new points replace all worst points in the original population, thus creating a new population also called the next generation.All points that are ranked better than Rma x in the previous population are passed unchanged onto the next generation.The new population is then re-sorted by way of Pareto ranking.The steps of Pareto ranking, simplex generation, and simplex evolution constitute a main loop of the MOCOM algorithm; iteration through these steps evolves the population in the direction of the global Pareto optimal set.The algorithm terminates, in a natural manner, when all points in the population have rank equal to one (i.e.there exist no inferior points in the population), which indicates that no further global directions of improvement can be found by this strategy.

Multi-objective downhill simplex evolution (MOSIM)
The evolution procedure used to improve the worst point is a multi-objective extension of the downhill simplex method (Nelder and Mead, 1965) and is therefore named MOSIM.The MOSIM method generates a new point that replaces the point in the simplex which has worst global rank.The new point is generated by the application of two operations, namely reflection and contraction, to the simplex as follows: where Sne w is the new point, Sw is the worst ranked point in the simplex, and Sg is the location of the centroid of the n better ranked points in the simplex.When 3' = 2, a reflection step is obtained and, when 3' = 0.5, a contraction step is obtained.
The rule for choosing either the reflection (Srer) or contraction (Soon) point is based on the concept of multi-objective dominance.The reflection point Srer is accepted if (and only if) it is non-dominated with respect to the n points that were used to compute the centroid.However, if the reflection point is found to be dominated, it is rejected, and the contraction point Sco, is accepted instead.This process causes each simplex to evolve exactly once and produce a single offspring.
An illustration of the MOSIM procedure for a twoobjective (FI,F2) problem with two parameters is shown in Fig. 5.The objectives to be minimized are quadratic functions centered at (0,0) and (1,0) and represented by the contour lines in Fig. 5

Simple illustration of MOCOM-UA evolution
The progression of the MOCOM-UA algorithm is illustrated in Fig. 6 for the simple two-dimensional two-objective (two quadratic functions) problem of Fig. 5.In Fig. 6(a) and (b), an initial population of 30 randomly selected individuals in the feasible space is shown.The four black dots indicate non-dominated points, while the 26 circles denote dominated points.After one loop (Fig. 6(c) and (d)), there are seven nondominated individuals (black dots).With successive loops, the Pareto front approximated by the black dots "close in" on the true Pareto set in both parameter and objective space.The final points are all mutually non-dominated and provide a close approximation of the Pareto region (Fig. 6(e)) and the tradeoff curve (Fig. 6(t)).

The model and data
The application of the MOCOM-UA algorithm to hydrologic modeling is illustrated by using it to calibrate the Sacramento Soil Moisture Accounting model (SAC-SMA) of the National Weather Service River Forecasting System (NWSRFS) using historical data from the Leaf River watershed (1950km 2) located north of Collins, Mississippi.A reliable 40water-year data set that represents a variety of hydrological conditions and phenomena is available for this watershed.The SAC-SMA model is a conceptual rainfall-runoff model with 13 parameters to be determined by the process of calibration.Although approximately eight consecutive relatively wet water years of daily data are required to ensure reliability of the parameter estimates (Yapo et al., 1996), the illustrative study presented here uses only a one year calibration period consisting of the wettest water-year on record (October 1, 1979to September 30, 1980), i.e. the year with the largest mean annual flow.Because the SAC-SMA model and the Leaf River data have been discussed extensively in previous work (e.g.Burnash et al., 1973;Peck, 1976;Kitanidis and Bras, 1980a, b;Brazil and Hudlow, 1981;Brazil, 1988;Sorooshian and Gupta, 1985;Sorooshian et al., 1982Sorooshian et al., , 1983Sorooshian et al., , 1993;;Duan et al., 1993Duan et al., , 1994;;and Yapo et al., 1996), the details of these will not be described here.

Objective functions
Automatic calibration procedures typically require the selection of a single-objective function for optimization in order to identify the best set of model parameters.This implies proper choice of an objective function.For hydrological model calibration, the objective function can be chosen to match some assumption regarding the distribution of the errors present in the observed data.For instance, the DRMS criterion: is the unbiased, minimum variance estimator, and it is the Maximum Likelihood Estimator under the assumption that measurement errors, approximated by et = q~im-qtbS, are normally distributed with zero mean and constant variance a 2. Similarly, the HMLE criterion: _ n sire obs 2 1 Xwt(X) [q, (O)-qt ] minimize HMLE(0, ~.) = n t= 1 (5) is the Maximum Likelihood Estimator under the assumption that the measurement errors are normally distributed with zero mean but having heteroscedastic (non-stationary) variance a 2 proportional to the observed flows (Sorooshian and Dracup, 1980).The weights are computed from the relationship wt(~,)=f 2~x-l), where f~ is the expected value of obs qt .Sorooshian and Dracup (1980) suggested using ft =q~im; however, it has been our experience that, in ~ t, obs practice, the use orjt=qt results in a more stable estimator.
In previous work (Yapo et al., 1996), we have shown that these two objective functions tend to provide different parameter estimates which result in different simulated hydrographs.The DRMS criterion  tends to emphasize minimization of peak flow error, while the the HMLE tends to provide more consistent performance across all flow ranges.In this study, we calibrated the SAC-SMA model using these two objective functions simultaneously in a multiobjective context and solving for the Pareto solution set using the MOCOM-UA algorithm.The multiobjective problem was defined as: minimize F=(DRMS, HMLE) (6) (0, x)

Calibration results
A population size of 500 was used to estimate the Pareto solution set with respect to the 13 model parameters.The algorithm required 68 890 function evaluations to converge to a solution.The results of this study are presented in Figs.7 and 8.A linear-linear plot of the Pareto solutions in the DRMS versus HMLE objective function space is given in Fig. 7(a).Notice that the MOCOM-UA algorithm has generated a uniform set of solutions spanning the Pareto space.The tradeoff between the two objectives is clearly illustrated: within the Pareto region, an improvement in one objective is obtained only at the expense of the other.For clarity, the 10 best DRMS and HMLE solutions are highlighted with different colored shading.A plot of the corresponding parameter sets is illustrated in Fig. 7(b).It is encouraging that the Pareto solutions tend to cluster closely in the parameter space.However, certain parameters (notably PCTIM, LZTWM, LZPK, LZSK, and PFREE) tend to be different for the two objectives; these parameters play a primary role in determining  the shape of the hydrograph during recession (between storm) periods.A log scale plot of the observed hydrographs (circles) and the hydrograph uncertainty (shaded area) associated with the Pareto solution set is displayed in Fig. 8.The larger relative uncertainty found during low flow and recession periods is consistent with the relative uncertainty in the parameters mentioned above.Notice also that the hydrograph uncertainty ranges do not bracket the observed flows during the low flow and recession periods, while the matching for medium and high flow events is very good.Clearly, the SAC-SMA model, as calibrated here, is better able to reproduce storm period flows than recession period flows, and these results may point towards aspects of the model and/or the calibration procedure that need to be further refined.We stress the illustrative nature of this study; clearly, issues such as the proper choice of objective functions for use in the multi-objective context need to be addressed before any conclusions regarding the model itself can be derived.

Sensitivity of the results to population s&e
The only MOCOM-UA algorithm parameter that must be specified by the user is s: the population size.To study the sensitivity of the calibration results to population size, we conducted a series of optimization runs in which all the conditions were identical, except that the population size was varied as s = 50, 100, 200, 500, 1000, and 2000.The final results, shown as a linear-linear plot in the objective space (DRMS versus HMLE) in Fig. 9, indicate that: 1.The Pareto front moves closer to the origin, indicating improvement in the results as the population size increases.2. The improvement is significant as s increases from 50 to 1000 points; however, the benefit of going from 1000 to 2000 can be considered to be marginal, given the order of magnitude increase in number of function evaluations (from 265 124 to 2 819 128).3. The populations seem to span the Pareto front evenly.4. The result obtained with 500 points approximates the best solution (2000 points) over most of the Pareto frontier, except at the lowest DRMS values.This supports the findings of Sorooshian and Dracup (1980) and Sorooshian and Gupta (1983) that the response surface associated with HMLE is smoother and better behaved than that of DRMS.
These results suggest that a population size of only 500 to 1000 is required to guarantee a reliable estimate of the Pareto front.This is encouraging because the number of required function evaluations increases dramatically (by a factor of 10) from a relatively reasonable 265 124 for s = 1000 to 2 819 128 for s = 2000.

Conclusions
Research into the development of automated (computer-based) calibration methods has focused mainly on the selection of a single-objective measure of the distance between the model-simulated output and the data and the selection of an automatic optimization algorithm to search for the parameter values which minimize that distance.However, practical experience with model calibration suggests that no single-objective function is adequate to measure the ways in which the model fails to match the important characteristics of the observed data.Given that many of the latest hydrologic models simulate several of the watershed output fluxes (e.g.water, energy, chemical constituents, etc.), there is a need for effective and efficient multi-objective calibration procedures which are capable of exploiting all of the useful information about the physical system contained in the measurement data time series.
In this paper, we have presented an effective and efficient methodology for solving the multipleobjective global optimization problem and demonstrated its features and capabilities by a simple example involving calibration of a conceptual rainfallrunoff model using two objectives.The MOCOM-UA multi-objective optimization method is relatively simple to implement.However, further research into a number of theoretical and experimental issues related to its application in hydrologic model calibration are ongoing.These include: (1) the proper manner for selecting the set of objective functions, and (2) the sensitivity of the results to number of objective functions and amount of data.In collaboration with colleagues, the MOCOM-UA approach is currently being applied to some of the more sophisticated physically based hydrologic models such as soil-vegetation transfer schemes and hydrochemical watershed models.The results of these studies will be reported in due course.We welcome dialogue on these and other ideas related to hydrologic model calibration.The code for the MOCOM multi-objective optimization algorithm is available from the second author by request (send email to hoshin@hwr.arizona.edu).

FigF1Fig. 1 .
Fig. 1.Illustration of Pareto optimality.this for a simple problem where two objectives 0cJ2) are to be minimized with respect to two parameters (Xl,X2).In Fig. l(a), point A indicates the parameters (X~l,X A) that minimize function fl, while point B B B indicates the parameters (xl, x2 ) that minimize function f2.The solution to the multi-objective problem consists of all the points falling on the line joining A and 13.These points have the characteristic that moving along the line from points A to B results in successively larger values of fl and successively smaller values off2 (see Fig. l(b)).Because the individual objectives are non-commensurate (each function is its own entity), it is not possible to determine objectively among such solutions which one is the best, and such solutions are called Pareto, nondominated, non-inferior, or efficient solutions.The formulation presented above divides the feasible parameter space into two parts: the part consisting of all Pareto solutions 0p is called the Pareto set, and the part consisting of the remaining space (non-Pareto solutions 0d) is called the dominated set.By definition, any member 0p of the Pareto set has the properties that:

Fig. 4 .
Fig. 4. Illustration of simplex selection in a two-objective problem.

Fig. 7 .
Fig. 7. Pareto solution obtained by application of MOCOM-UA to calibration of the SAC-SMA model: (a) objective space, and (b) parameter space.
reviewers, Dr Richard Ibbitt and Dr Qingyun Duan, for their perceptive and positive comments which helped us improve the paper.Partial financial support for this research was provided by the National Science Financial assistance provided to Dr Yapo by The University of Arizona Graduate College is gratefully acknowledged.