xcai@mail.utexas.edu

lasdon@mail.utexas.edu

Nonconvex nonlinear programming (NLP) problems arise frequently in water resources management, e.g., reservoir operations, groundwater remediation, and integrated water quantity and quality management. Such problems are usually large and sparse. Existing software for global optimization cannot cope with problems of this size, while current local sparse NLP solvers, e.g. MINOS (Murtagh and Saunders, 1987), or CONOPT (Drud, 1994) cannot guarantee a global solution. In this paper, we apply the Generalized Benders Decomposition (GBD) algorithm to two large nonconvex water resources models involving reservoir operations and water allocation in a river basin, using an approximation to the GBD cuts proposed by Floudas et. al.(1989) and Floudas (1995). To insure feasibility of the GBD subproblem, we relax its constraints by introducing elastic slack variables, penalizing these slacks in the objective function. This approach leads to solutions with excellent objective values in run times much less than the GAMS nlp solvers MINOS and CONOPT2. Using these solutions as initial points for MINOS or CONOPT2 often leads to further improvements.

For many water resources management models, nonlinear programming (NLP) offers a general mathematical formulation handling a non-separable objective function and nonlinear constraints. However, in many cases the NLP’s are nonconvex, but local solvers are applied. Methods for obtaining global solutions to nonconvex mathematical programming problems have not appeared frequently in the water resources literature, even though these problems are often unavoidable in water resources management modeling. In recent years, genetic algorithms (GAs) have been proposed as a promising method to solve nonconvex NLP problems in water resources systems planning (McKinney and Lin, 1995; Ritzel et al., 1994). Watkins and McKinney (1997) used Generalized Benders Decomposition (GBD) and Outer Approximation (OA) to solve a mixed-integer nonlinear (MINLP) water resources optimization model with fixed costs.

In this paper, we apply the GBD technique to two large-scale, nonconvex NLP problems in water resources management, in an approximate way suggested in Floudas et. al.(1989) and Floudas (1995). To demonstrate the modified approach, we solve two large NLP models in water resources management. One is a reservoir system operation model, considering both water supply and power generation. The other is a river basin water allocation model, including salt control. In Section 2, we describe these two classes of models. In Sections 3 and 4, we apply Generalized Benders Decomposition (GBD) to our models, relaxing some of the sub-problem constraints to insure feasibility. We show that the true GBD cuts are partially separable piecewise linear functions, so the GBD relaxed master program is a reverse convex program. Our computations use linear approximations to the true GBD cuts, as suggested in Floudas et al. (1989) and Floudas (1995). Our GAMS implementation is described in section 5. The results of Section 6 show that, for the problem instances we consider, GBD using these cuts lead to objective values only slightly worse than those achieved by MINOS and CONOPT2 (Drud, 1994), in much less computing time.

In this section we formulate two typical water resources management models, one a reservoir operation model and the other a river basin water allocation and salinity control model.

2.1.1 Sets and Indices

2.1.2 Data

2.2.1 Objective Function

There are several objectives to consider in the management of a river basin.

• Maximize the ratio of water delivered to that demanded at demand sites

which insures that no demand is oversatisfied.

• Maximize the smallest water deficits among all demand sites and time periods

• Maximize the sum of the total amount of hydropower generated in the basin and the minimum hydropower generation over all periods

• Minimize salt concentration in the system

where flows out of demand sites represent return flow from demand sites, including surface and sub-surface drainage.

• Hydroelectric power generation

which expresses hydropower generated as being proportional to the flow through the turbine times the difference between average surface elevation and tail water elevation tw(n).

• Reservoir head - volume relationship:

• Salt balance at river, canal, or drainage nodes

• Salt balance at reservoir nodes

• Salt balance at groundwater nodes

• Salt balance at demand nodes

(2.18)

• Variable bounds

All variables are nonnegative and some have specific lower and upper bounds to represent some physical or policy limits. For example, flows through river nodes have a lower bound for some instream water uses such as flow augmentation and recreation, and an upper bound for flood control; the ratio of water supply to water demanded cannot be lower than a specified value in some months; salt concentrations in some locations have upper bounds.

For this model, the objective function is a weighted combination of objectives (2.1), (2.3), and (2.7):

where the non-negative weights w_i reflect the relative importance of the individual objectives. The constraints are the water and salinity balance constraints (2.8) – (2.12) and (2.15) – (2.18), omitting the power generation constraints (2.13)-(2.14).

For the reservoir operation model, the objective function is a weighted combination of objectives (2.1), (2.3), and (2.5):

The constraints are the water balance constraints (2.8)-(2.12) and the power generation constraints (2.13) – (2.14), omitting the salinity balance constraints (2.15)-(2.18).

Generalized Benders Decomposition (GBD), developed in Geoffrion (1972), was originally derived for problems of the form:

s.t.

(3.1)

where g is the vector of coupling constraint functions. The vector y contains complicating variables, in the sense that the problem is assumed to be considerably easier to solve when y is fixed. At iteration r, the procedure involves successive solutions of the subproblem:

subject to (3.2)

whose optimal objective value is v(y_r-1), and the relaxed master program

RMP(r): maximize y₀

subject to

and

whose solution is . The function is defined as:

and is an optimal multiplier vector for .

In the following, we assume that has an optimal solution and an optimal multiplier vector for each . This is guaranteed for our problems by relaxing some constraints of P(y) by introducing "elastic" slack variables where necessary. Under these assumptions, the steps of the GBD algorithm are:

2. X is nonempty and convex.

3. It must be possible to solve each RMP(r) globally.

   In order for GBD to be computationally attractive, the following additional assumption must be satisfied:

4. It must be possible to obtain a closed form expression (or one that can be evaluated rapidly) for the GBD cut functions for each fixed . This is "property P" as defined in (Geoffrion, 1972)

   In our problems, assumptions (1) and (2) are satisfied because is bilinear, X and Y are polyhedral, and f is linear. Some of the coupling constraints (3.1) are equalities, but these define a convex region in x space for fixed y because g is bilinear. We show in the next section that assumption (4) is also satisfied, but (3) is not, because each RMP is a reverse convex program.

   4. Applying GBD to Water Resources Management Models

   All the constraints described in Section 2.2 are linear except for power generation (2.13) and salinity balance (2.15)-(2.18) which are bilinear, i.e. all nonlinear terms are products of two variables. In the salinity balance constraints, we find products of flows Q(n, n₁, t) and storages S(n,t) with salt concentrations C(n, t), while the hydropower generation relations involve products of reservoir outflows Q(n, n₁, t) and reservoir heads H(n, t). Hence if one member of the set of variables ((Q,S), C) is fixed, the salinity balance constraints become linear, and if one member of the set ((Q,S), H) is fixed, the power generation constraints are linear. This structure motivates the application of GBD, and the selection of either the (Q,S) or C variables as complicating variables y in the salinity control model, and either (Q,S) or H as complicating variables in the reservoir operation model. As pointed out in (Floudas, 1989), these selections apply to any problem with bilinear objective and/or constraint functions, and theorem 1 below applies to any such problem.

   4.1 Salinity Control Model

   It is natural to choose the salt concentrations C as the complicating variables y, since there are fewer of them. Surprisingly, this is not the best choice. The computational results of Section 6 show that the approximate GBD cuts described by Floudas (1989 and 1995), using this choice of y, lead to solutions of quality comparable to those achieved by the GAMS solvers MINOS or CONOPT2, using the same starting points, in 15 or fewer major iterations. Choosing C as y, GBD using these approximate cuts did not converge.

   Let x be the vector of salt concentration variables C(n, t) and y be the vector of all other salinity model variables (Q, RI, S, G). Then most of the salinity model constraints, i. e., the water balance constraints (2.8) – (2.12), and the constraints defining the objective (2.1)-(2.4) involve only y, and these are linear. Letting Y represent the polyhedron defined by these constraints, they can be written as:

   (4.1)

   The salt balance constraints (2.15) – (2.18) involve both x and y, so these are the coupling constraints g(x, y) of the problem (3.1). In these constraints, each concentration, x_j, multiplies one or more flow or storage variables (components of y). Hence we can write the vector of salt balance constraints as

   (4.2)

   where b is a vector of constants, a_j(sy_j) is the column of coefficients associated with x_j in the subproblem P(y), and sy_j is the subvector of y upon which a_j depends. If has nonzero entries as shown in Table 2 below:

   Row	Coefficient
   (n,t)
   	Q(n,n1,t)
   (n, t+1)	S(n,t)

   Table 2 – Nonzero coefficients of C(n,t)

   The variable S(n,t) does not appear if n is not a storage node. No other column a_j involves these variables, because only C(n,t) multiplies flows out of node (n,t) and S(n,t). Hence a_j(sy_j) is a linear function of sy_j, and the subvectors sy_j partition y.

   It is desirable to insure that the subproblem constraints, (4.2) and , are feasible for all fixed yÎ Y, in order to avoid having to find an extreme ray of the GBD subproblem when they are infeasible. This is typically done by introducing two nonnegative vectors of "elastic" or "deviation" variables dp and dn, rewriting (4.2) as

   (4.3)

   Then the objective becomes

   (4.4)

   where M is a sufficiently large positive constant. This form for the penalty term is exact so, if the problem is feasible, there is a positive threshold value for M (the max absolute optimal multiplier for (4.2)) such that any value above that threshold yields optimal solutions with all elastic variables equal to zero and the remaining variables optimal for the original problem. To complete the model statement, the only constraints involving the salt concentrations C alone are simple upper and lower bounds, written as

   (4.5)

   Applying the GBD algorithm as stated in Section 3 to this form of the model, the subproblem P(y) is to maximize (4.4) over x, dp, and dn for fixed y subject to (4.3) and (4.5). P(y) is a linear program. The Lagrangian function for P(y) is

   (4.6)

   where u is a vector of Lagrange multipliers for (4.3), introduced with a minus sign to give them the same signs as the simplex multipliers of P(y). The GBD cut function is

   . (4.7)

   Collecting terms in y, x, dp, and dn yields:

   . (4.8)

   The coefficient of x_j above is the reduced cost or dual slack corresponding to x_j, defined as

   .

   Since L is linear in x, dp, and dn for fixed y and u, the max in (4.7) may be applied separately to each term of (4.8), so

   

   (4.9)

   By complementary slackness, the maxima over dp and dn in (4.9) are zero. For any fixed y and u, the max over each x_j occurs either at lx_j or ux_j so

   . (4.10)

   Thus we can characterize as follows:

   Theorem 1. is a convex piecewise linear function of y, and is separable in the subvectors sy_j for any fixed u.

   Proof

   In (4.10), the max of the pair of linear functions is convex and piecewise linear. Denoting this max by pl_j(sy_j ,u) we have

   (4.11)

   Since is a sum of convex functions, it is convex. Since the sy_j partition y, is separable in these sub-vectors.

   --------------

   Our previous discussion of the column a_j showed that sy_j has dimension equal to the number of successors of the node associated with x_j. This dimension is almost always small (e.g., <10), independent of network size. Hence the piecewise linear functions pl_j in (4.11) are easy to evaluate.

   The master program RMP(r) for the water allocation model (or any problem with bilinear nonlinearities) is a nonconvex nonlinear program, known as a "reverse convex" program (Horst, 1992). As such, it may have distinct local optima. Its constraints include the inequalities

   

   Since is convex but, in general, not linear, is a convex nonlinear function of The inequalities in RMP(r) thus have the "wrong" direction to define a convex set. Each inequality defines the complement of a convex set. There is no known polynomially bounded algorithm for solving the reverse convex RMP(r) globally, and GBD cannot guarantee a global solution of the original problem without such a guarantee for RMP(r).

   Assume that each upper bound in ux is greater than the corresponding component of lx. Then the max in (4.10) is attained at lx_j if ds_j is nonpositive, and at ux_j if ds_j is nonnegative. By restricting the signs of each ds_j we force the max to occur at a particular bound, and these restrictions define polyhedra over which L* is linear. Let xe^k be the kth extreme point of X and define

   (4.12)

   Then, for all y in P_k(u), L* is the linear function

   (4.13)

   If x has nx components, there are at most 2^nx linear pieces . In (Visweswaren and Floudas, 1996), a branch and bound algorithm for finding a global solution to RMP(r) is devised based on this piecewise linear structure.

   Because our problems are so large, we have opted for a simpler but approximate approach. Floudas ((1989) and (1995)) has suggested that the functions L* be approximated as follows: let and let be primal and dual solutions for P. Then use

   L*(y, = L (4. 14)

   This holds for problems where f and g are separable in x and y. However, for nonseparable problems, while (4.14) holds at , it does not hold in general at other y, since need not maximize L for y different from . In general, the linear function on the right hand side of (4.14) does not coincide with any piece of L*, unless the reduced costs ds_j satisfy special conditions like those in (4.12). For example, if is strictly between its bounds, then . However, as y varies while u remains fixed, is no longer a maximizer unless ds_j remains zero. Hence these approximate cuts need not provide valid upper bounds. Nonetheless, we have obtained good solutions using them, as we discuss in section 6. However, our results include instances where the optimal RMP value is lower than the objective value of a locally optimal solution found by another solver.

   4.2 Reservoir operation model

   The previous results also hold for the reservoir operation model. We choose y = H and x = (Q, RI, G, S, P). Using (4.14) as the GBD cut leads to feasible solutions of high quality, as we describe in Section 6. With y = H, the coupling constraints are the hydropower generation constraints (2.13) and the reservoir head-volume relationship (2.14), and these, plus the water balance constraints (2.8)-(2.12) are the constraints of the GBD subproblem. They are linear for fixed y. As with the salt balance constraints of the salinity control model, the constraints (2.13) are relaxed by introducing elastic variables, insuring that the subproblem is feasible for all y in Y. The set Y is determined by the upper and lower bounds on y. The GBD relaxed master program involves only the GBD cuts in the approximate form (4.14), and the bounds on y. Both master and subproblem are linear programs.

5. Implementation

This algorithm has been implemented using the algebraic modeling language GAMS (Brooke et. al., 1988), using GAMS version 2.50. Both subproblem and relaxed master program (RMP) models are defined, and a loop statement is used to drive the GBD algorithm. This loop contains two SOLVE statements, one for the master and subproblem, both using the OSL LP solver. The GBD cut is created by a GAMS statement which involves the optimal values of the primal and dual variables from the previous subproblem solution, and the new cut is indexed by the loop index. The GBD termination criterion is that of step 5 of Section 3, with tol = 1.0E-3.

Our GAMS program exploits the fact that the GBD subproblem structure remains the same at each iteration, but with different parameters (values of y), and that only one equation is added to the relaxed master program (RMP). All solutions of these linear programs after the first take advantage of GAMS' automatic warm starting capabilities. The optimal basis from each LP is saved and used as a starting basis for the next SOLVE. This restarting facility saves significant computing time (Brooke et. al., 1996).

As to initial values for y, the reservoir model instances solved here have either 12 or 24 complicating variables, the heads of one or two reservoirs in each period, and it is easy to provide good initial guesses for these. The water allocation model has many more sets of complicating variables, e.g. all flows, storages, etc. Some instances solved in section 6 have 1499 complicating variables. Initial values for y are best chosen by first solving the relaxed master program with no GBD cuts, i.e. choosing the flow and storage variables to optimize the portion of the model in which only they appear. Section 6 also shows results where initial y’s are selected by assigning identical "ballpark" values to each variable within a set of related variables, e.g., flows into and out of river nodes, from rivers to canals, etc.

Each model discussed in Section 2 is first expressed directly, and either MINOS or CONOPT2 (the new version of CONOPT available with GAMS 2.50) is used as a nonlinear solver for both models. Then the models are implemented in the GBD form described above. The results for the two models, solved by either MINOS or CONOPT2 and by GBD, are discussed in the following section.

The instances of the reservoir system operation model used in our computational experiments are described in McKinney and Cai (1997). These models have been used to analyze the tradeoff between upstream power generation and downstream irrigation in the Syrdarya River basin in Central Asia, and the model data is specific to this region (McKinney and Cai, 1997). This river basin network includes 13 reservoirs, 30 river reaches, 6 aquifers, 6 demand sites (irrigation, municipal, and industrial water demand), and 5 hydroelectric power generation stations.

Among the 5 hydropower stations, the linear head/storage relations (2.14) were only defined for the Toktogul (14.5 km³ of active storage volume) which is the major reservoir for flow control and power generation in the basin and the Andijan reservoirs !.64 km³ of active storage volume). The reservoir surfaces of the other hydropower stations are maintained at constant elevation, and thus their heads are treated as constants. The time periods are months, with horizons ranging from 12 to 60 months. The model with 12 time periods consists of 895 equations, 1151 variables, and 2673 nonzero Jacobian elements. Considering the monthly average heads of the Toktogul Reservoir for one year as complicating variables, with the Andijan elevations taken as constant, there are 12 complicating variables and 24 coupling constraints, 12 power generation constraints (2.13), and 12 head/volume equations (2.14). If the Andijan heads are allowed to vary, the number of complicating variables increases to 24, and the number of coupling constraints becomes 48.

Figure 1 shows the lower bound, upper bound, and penalty cost vs. iterations for the model with 12 complicating variables. In this case, the penalty cost goes to zero, the upper bound always decreases, and the lower bound fluctuates until shortly before the penalty reaches zero, increasing thereafter until it reaches the upper bound. Figure 2 shows the upper and lower bounds vs. iterations for two sets of initial values for the Toktogul Reservoir heads in all time periods. In the first scenario, the heads in all 12 periods were 836 m, which corresponds to the dead storage of the reservoir. In the second, the initial values are all 880 m, which corresponds to the full storage of the reservoir. In both scenarios, the gap between the 2 bounds is substantial until around iteration 13, when it shrinks sharply as the lower bounds improve, and then converges slowly to zero at iteration 37. The upper bounds provided by the relaxed master programs are nearly optimal after about 10 iterations. Solving this model directly using MINOS with the same starting points leads to the same solution.

Figure 3 shows the lower bound, upper bound, and penalty cost vs. iterations for the model with 24 complicating variables. Here convergence is much slower, requiring about 155 iterations, with computation time increasing accordingly. The penalty weight, M, also strongly affects GBD convergence and computing time. Figure 4 shows that, for weights of 1 and 10, the penalty cost converges to zero, while for M = 0.1 it converges to 1.6. Figure 5 shows the effect of the penalty weight on GBD convergence. The run with M = 10 converges fastest, while GBD does not converge for M = 0.1.

In these figures, MINOS is used as the LP solver. If MINOS is replaced by OSL, the number of iterations does not change significantly (for the model with 24 complicating variables, the number of iterations is 155 for MINOS, and 137 for OSL), but the total computing time is reduced to about 25% of the MINOS time.

The river basin water allocation model used in this paper was developed by McKinney et al. (1997) for the Karshi region of the Amudarya River basin of Uzbekistan. Its mathematical structure is described in Section 2. It is used to support water allocation decisions with multiple goals including satisfying water demand, maintaining river flow for ecological protection, balancing water use rights among demand sites, and controlling salinity. The network used to model this basin includes 6 reservoirs, 6 aquifers, 8 river reaches, 2 canals, 7 agricultural drainage water collectors, and 10 demand sites (irrigation, municipal and industrial water demand). There are 12 monthly time periods. When salt concentrations are included in the objective, this model instance has 1567 constraints, 2039 structural variables (not including the deviation variables), and 9129 nonzero Jacobian elements, of which 4773 (about 52%) are nonconstant. Hence this model is highly nonlinear, due to the presence of 540 nonlinear salt balance constraints. There are 1499 complicating variables y, 540 salt concentrations, and 540 coupling constraints.

The storage of each reservoir and aquifer is divided into a constant part, called dead storage, and a variable part called live storage. Only the live storage can be accessed, so only live storage variables appear in water balance equations. However, salt balance equations for reservoirs and aquifers involve both types of storage, so they contain terms where this constant storage volume is multiplied by variable salt concentrations. Hence these salt balance equations contain linear terms involving x alone. The objective term also introduces a linear term in x into the Lagrangian. We conjectured that such terms would degrade the performance of GBD, because they would force the optimal flows and storages to be further from their starting values, those which optimize the portion of the model not including salt concentrations. However, the salt concentrations of these aquifers and reservoirs do not change much over time, so deleting the dead storage terms has little impact on the model’s solution. In addition, we can delete or include the salt concentration terms in the objective, which has a greater effect. Hence, we can define 4 cases as shown in table 1, representing increasing distance from the "best" situation of case 1. In each case, we first solve the problem using GBD and then, using the final GBD solution as an initial point, solve it using MINOS and (separately) CONOPT2. In all GBD runs, the initial y values are computed by solving the relaxed master program with no cuts, the "optimal flow" solution. Despite the excellent initial points, MINOS and CONOPT2 were unable to improve on the first four significant figures of the final GBD values in cases 1 and 2, so this value is at least locally optimal to within the rather tight default tolerances of CONOPT2. There are slight improvements to cases 3 and 4, confirming our hypothesis that these cases are less favorable to GBD. Still, the largest improvement, by CONOPT2 in case 4, is only 0.14% better than the GBD value. The CONOPT2 final objective value for case 4 of 4.981 is is not the best we have obtained for this problem instance-CONOPT2 attains 4.9838 using the starting point of case 4-2 in table 2. It is unclear if these values represent 2 distinct local optima or if the objective is just very flat near these points.

As discussed above, the most effective way to choose initial y’s is to solve the relaxed master program with no GBD cuts. Figure 6 shows the behavior of the lower bound, upper bound and the penalty term versus iteration count for the model of case 4 of table 1 when the initial y is chosen in this way. This model, with 1499 complicating variables, converges in 13 iterations, much less than the 155 iterations required for the reservoir operation model, which has only 24 complicating variables. We believe that this faster convergence is due to the "tighter" relaxed master program of the water allocation models. The RMP of the reservoir operation model has only one set of constraints, the GBD cuts, with a new cut added at each iteration. In early iterations, when the RMP has few constraints, the proposals provided by this RMP are far from optimal. For the water allocation model, the RMP includes the flow balance equations, which constrain this RMP more tightly than the RMP of the reservoir operation model. The tighter RMP provides proposals that are of higher quality, resulting in faster convergence.

Table 2 shows the final objective values, GBD iterations, and run times of GBD, MINOS and CONOPT2 when they use the same initial point when solving case four (the most realistic case) in table 1. MINOS and CONOPT2 use all default tolerances and options. Run times shown are the "resource usage" times reported on the GAMS listing file for the solution phase. The computer used is a pentium II 300 mhz PC. Four different initial points are used. The best initial point, case 4-1, chooses initial y’s as the "optimal flow" solution described above, and initial x’s by solving the GBD subproblem. The other three cases use different "ballpark" initial y values as discussed in section 5, again choosing x by solving the GBD subproblem. GBD produces objective values slightly worse than the other 2 solvers in all cases, but the largest difference, in case 4-3, is only 0.74%. GBD is much faster than the other two solvers, often by more than a factor of 20 over MINOS, and by factors of 3 to 9 over CONOPT2. GBD time is affected very little by the starting point. In the column "MINOS, Third Run", we show the final objective value resulting from 3 successive applications of MINOS, using 3 consecutive GAMS SOLVE statements. SOLVE 1 uses the same initial solution as GBD, and SOLVES 2 and 3 use the result of the previous SOLVE as starting points. This easy-to-implement strategy improves the MINOS objective value achieved with one SOLVE by 1.56% in case 4-3. Using a single SOLVE, CONOPT2 achieves slightly better objective values than the other solvers in two of the four cases, but fails to find a feasible solution in case 4-2. However, CONOPT2 achieves the best objective value of all using three successive SOLVES in this case.

Using a relaxed formulation of the GBD subproblems and a sufficiently large penalty weight, GBD has performed well in solving the large, nonconvex, bilinear problems studied here. For the water allocation models, GBD solution quality is comparable to that of MINOS or CONOPT2, and GBD is considerably faster. Further, GBD can be used to advantage in conjunction with these or any other local solvers, by using the final GBD solution as an initial point for the local solver. In our experiments with the water allocation model, MINOS and CONOPT2 were able to improve this solution to a small degree. Computation time for GBD solution of the reservoir operation models increases as the number of complicating variables increases from 12 to 24. However, experience with the water allocation models indicates that GBD can handle roughly 1500 complicating variables in situations where the relaxed master program is tightly constrained by the constraints defining Y. Hence analysts considering using GBD to solve nonlinear problems should consider carefully which variables are designated as complicating, and consider reversing the "natural" assignment which tries to minimize the number of these variables.

Brooke A., D. Kendrick, and A. Meeraus. 1996. GAMS: A User’s Guide. GAMS Development Corporation, Washington, DC.

Drud, A. S. 1994. CONOPT - A Large Scale GRG Code. ORSA J. on Computing 6, 207-216.

Floudas, C. A., A. Aggarwal, A. R. Ciric. 1989. Global Optimum Search for Nonconvex NLP and MINLP Problems. Computers Chem. Eng. 13, 1117-1132.

Floudas, C. A. 1995. Nonlinear and Mixed-Integer Optimization. Oxford Univ. Press, New York.

Geoffrion, A. M. 1972. Generalized Benders Decomposition. J. Optimization Theory Applic. 10, 237-259.

Horst, R., and Hoang Tuy. 1996. Global Optimization-Deterministic Approaches. Springer, New York.

McKinney D. C., and M. D. Lin. 1994. Genetic Algorithm Solution of Groundwater Management Models. Water Resources Res. 30, 1897-1906.

McKinney D. C., Karimov, A. Kh., and Cai, X. 1997. Model development: Aral Sea Regional Allocation Model for the Amu Darya River. Project Report, Environmental Policy and Technology Project, US Agency for International Development, Central Asian Regional Office, Almaty, Kazakstan.

McKinney D. C., and Cai, X. 1997. Multi-Objective Water Resources Allocation Model for the Toktogul Reservoir. Technical Report, Environmental Policy and Technology Project, US Agency for International Development, Central Asian Regional Office, Almaty, Kazakstan.

Murtagh, B. A. and M. A. Saunders. 1982. A Projected Lagrangian Algorithm and its Implementation for Sparse Nonlinear Constraints. Mathematical Programming Study 16, 84-117.

Ritzel, B., J., J. W. Eheart, and S. Ranjithan. 1994. Using Genetic Algorithms to Solve a Multiple Objective Groundwater Pollution Contaminant Problem. Water Resour. Res. 30, 1589-1603.

Viswesaran, V. and C. A. Floudas. 1996. New Formulations and Branching Strategies for the GOP Algorithm. Global Optimization in Chemical Engineering, I. E. Grossman, ed., 75-100.

Watkins, D. W. and D. C. McKinney. 1997. Decomposition Methods for Water Resources Optimization Models with Fixed Costs. Advances in Water Res. to appear 1997.

Figure 1. The lower bound (LBD), upper bound (UBD) and penalty for the reservoir operation model with 12 complicating variables. The final solution converges to 18.05 at iteration 35.

Figure 2. Upper and lower bounds for various initial values of the complicating variables (reservoir surface elevation in each month). All scenarios converge to the same final value, 18.05.

Figure 3. Lower bound (LBD), upper bound (UBD) and penalty for the reservoir operation model with 24 complicating variables. The final solution converges to 18.77 at iteration 155.

Figure 4. Penalty vs. iteration for the reservoir operation model with 24 complicating variables under various penalty weights.

Figure 5. Upper bound vs. iteration for the reservoir operation model with 24 complicating variables using various penalty weights.

Figure 6. GBD Lower bound (LBD) and upper bound (UBD) for the problem in case 4 of table 1.

Table 1. GBD and MINOS final objective values for four water allocation model cases

Table 2. Performance of GBD, MINOS, and CONOPT using 4 different initial points on the problem defined as

case 4 of Table 1.

Case 4-1: initial y values obtained by solving the GBD master program with no cuts.

Cases 4-2, 4-3, 4-4 use 3 sets of "ballpark" initial guesses for y.

All programs were run on the same PC-300, Pentium-II. The computational time here is defined as the "resource usage" in the GAMS output file.