This paper proposes a hybrid mixed-integer linear programming-linear programming-linear programming (MILP-LP-LP) methodology for simultaneous optimal design and operation of unconfined groundwater utilization systems. The proposed model is an extension to the earlier LP-LP model introduced by the authors for the optimal design and operation of confined groundwater utilization systems. The proposed model can be used to minimize the total cost, including the well drilling, pump installation cost and energy cost, of utilizing a two-dimensional unconfined aquifer under both steady-state and transient flow conditions. The solution of the problem is defined by the well numbers and location, well drilling depth and the corresponding pumping rates, satisfying a downstream demand, lower/upper bound on the pumping rates, and lower/upper bound on the water level drawdown in the wells. A discretized version of the differential equation governing the aquifer is first embedded into the model formulation as additional constraints. The resulting mixed-integer non-linear programming problem is then decomposed into three sub-problems with different sets of decision variables, namely, hydraulic head, well numbers and locations and their pumping rates, and well depths. Starting with a set of random values for all decision variables, the three sub-problems are solved assuming fixed values for the other variables. This process is iterated until convergence is achieved.

## NOMENCLATURE

### Indexes

### Variables

### Parameters

*N*total number of potential well locations

*a*_{1}cost coefficient for drilling

*a*_{2}cost coefficient for well and pump installation

*a*_{3}cost coefficient for operation

*d*_{k}well's depth at position

*k**β*a constant

*Q*_{d}total demanded pumping rate

*h*_{min}minimum allowed hydraulic head due to environmental issues

*Q*_{min}minimum pumping rate limit according to the available pumps

*Q*_{max}maximum pumping rate limit according to the available pumps

*d*_{max}maximum allowed well depth

*m*total number of wells required to be drilled and operated

*K*_{h}hydraulic conductivity of the aquifer

*R*surface recharge rate

*x, y*independent Cartesian coordinates

*x*_{k}, y_{k}Cartesian coordinates of the well at position

*k**Δx, Δy*FD grid spacing

*δ*Dirac delta function

*α*penalty factor

initial hydraulic head in the aquifer

*S*_{y}specific yield coefficient of the aquifer

*T*total operation period

*Δt*duration of each time step

## INTRODUCTION

Groundwater is not only an important component of water resources, but also a reliable source of fresh water for a variety of purposes including domestic, industrial and irrigation uses that is known to mostly remain intact from pollution. Nowadays, with increasing population and living standards, there is a growing need for the utilization of groundwater resources. However, due to some anthropogenic causes such as unplanned urbanization and industrialization, the quantity of groundwater resources continues to decrease. Therefore, decision makers in the field should obtain proper design and management strategies to optimally utilize the groundwater resources and prevent adverse losses.

The groundwater utilization problems can be theoretically categorized into three main classes, namely design problems, operational problems, and design and operation problems. Design problems are defined as finding the number, location and depth of each well so that the total cost of the utilization system's construction including well drilling and installation is minimized' assuming a fixed and pre-defined operation policy. In operational problems, often referred to as groundwater management problems, the optimal operation of an existing utilization system represented by the pumping discharge from each well is sought so that the total operational cost is minimized over the service life of the project. In design and operation problems, as the most general form of groundwater utilization problems, both design parameters and operational policies are determined simultaneously to minimize the total cost of the project including both the construction and operational cost of the utilization system. Due to the strong coupling between the design and operation part of the problem and the fact that the operational cost constitutes the major part of any groundwater utilization systems' costs, the design problems are considered less practical and useful compared to the two other problems. Problems in this field are solved using different methods which are discussed below in two categories: operational problems and design and operational problems. In the field of operational problems, Katsifarakis & Karpouzos (1998) applied a simulation-optimization model based on a genetic algorithm (GA) as the optimizer and a boundary element method (BEM) as the simulator to minimize the pumping cost (i.e. pumping energy) of a pre-defined set of wells in a confined aquifer under steady-state flow condition. Katsifarakis *et al.* (1999) extended the GA-BEM model to solve three different problems, namely (a) determining the transmissivity of a zoned aquifer, (b) minimizing the pumping cost of a pre-defined number of fixed wells in a confined aquifer under steady-state flow, and (c) hydrodynamic control of a contaminant plume by means of pumping/injection wells. Katsifarakis (2008) proposed an analytical approach seeking the minimum pumping cost in an infinite and semi-infinite homogeneous confined aquifer for a pre-defined number of wells using the method of images. In the field of operational and design problems, Rosenwald & Green (1974) employed a transient response matrix and a branch and bound method to solve a mixed-integer programming problem to obtain the best locations for a specified number of wells. A harmony search (HS) algorithm was used by Tamer Ayvaz (2009) to obtain the solutions to three different problems of (a) total pumping rate maximization for a pre-defined number of wells, (b) minimization of the total utilization cost of a pre-defined number of wells, including well drilling, and capital and operating costs to satisfy the given demand, and (c) minimization of the total utilization cost of a pre-defined number of wells satisfying the given demand for multiple operation periods. Sidiropoulos & Tolikas (2004) used a combination of GA and BEM to determine the optimal well configurations and corresponding pumping rates to minimize the pumping cost from a confined aquifer to provide a given downstream demand. Tamer Ayvaz & Karahan (2008) implemented a moving sub-domain approach to determine well locations to be used for groundwater model calibration. The pumping rates were then obtained using a coupled GA-finite difference method (FDM) optimization-simulation model.

Groundwater utilization problems are typically solved by researchers using a simulation-optimization approach, in which an optimization model is coupled with a groundwater flow/transport model to establish the best utilization policy for the problem in hand (Gaur *et al.* 2011). Most real-world groundwater utilization problems are mixed-integer non-linear non-convex problems for which the use of gradient-based optimization algorithms may fail or result in unexpected solutions (McKinney & Lin 1994). These algorithms usually require a good initial guess to produce an optimal one. They also rely on local gradients of the objective function and/or the constraints to determine the search direction, and thus may converge to local optimal solutions. Therefore, use of heuristic algorithms is usually preferred due to their ability to find solutions without requiring an initial guess or any information on the gradients of the objective function or the constraints (Tamer Ayvaz 2009). Modern heuristic methods are shown to be effective in avoiding local optimums and dealing with mixed variables, but they are also known to incur very high computational costs due to their evolutionary nature. Furthermore, these models have some free parameters that should be tuned for the best performance of the algorithm before the main calculation, contributing to the high computational cost incurred by these methods. More importantly, modern heuristic search methods basically are unconstrained optimization solvers requiring the use of a penalty method to handle the constrained optimization problems introducing another free parameter into the method. Moreover, the application of modern heuristic search methods to problems with constraints cast in the form of a differential equation, such as the one considered here, requires a simulation/optimization approach which in turn increases the computational inefficiency of these methods as a result of calling the simulation model for every single function evaluation (Bostan *et al.* 2015). Asher *et al.* (2015) and Singh (2014) further reviewed groundwater modelling techniques and approaches.

In this paper, the hybrid LP-LP approach of Khadem & Afshar (2013), originally proposed for optimal operation of confined aquifer management systems, is extended for the simultaneous optimal design and operation of unconfined aquifer utilization systems. The problem to be addressed here is to optimally determine the design parameters, such as well numbers, depths and positions, and the operational parameters, namely the pumping rates of the operational wells, to minimize the total cost, including the well drilling cost, pump installation cost, and energy cost of utilizing a two-dimensional (2D) unconfined aquifer under both steady-state and transient flow conditions. For this, a discretized version of the differential equation governing the aquifer is first embedded into the model formulation as additional constraints. The resulting highly constrained mixed-integer non-linear programming (MINLP) problem is then decomposed into three sub-problems with different sets of decision variables, namely hydraulic head, well numbers, locations and corresponding pumping rates, and well depths. Starting with a set of random values for all decision variables, the three sub-problems are solved in turn using proper methods assuming fixed values for the other variables. This process is iterated until convergence is achieved. The efficiency and effectiveness of the proposed method is tested on two examples from the literature and the results are presented and compared with those of other methods. Results indicate the efficiency and robustness of the proposed method for design and operation of real-world unconfined aquifer utilization systems under steady state and transient flow conditions.

## PROBLEM FORMULATION

As economic aspects of utilization problems play an important role for decision makers, the objective of this utilization model is to minimize the total utilization cost of a 2D unconfined aquifer including well drilling, well installation, and pumping costs. The problem is defined as determining the numbers, position and the corresponding depth of the wells to be drilled out of a set of potential locations, and their operational policy. This is defined by their pumping rate over the operational horizon, such that the total cost is minimized while the operational and construction constraints are satisfied under steady-state/transient conditions.

### Steady-state flow condition

The problem of design and operation of a 2D unconfined aquifer under steady-state condition can be formulated as follows.

*Q*is the pumping rate of the well at position

_{k}*k*[

*L*],

^{3}T^{−1}*d*is the depth of the well at location

_{k}*k*[

*L*],

*h*is the hydraulic head at the well at location

_{k}*k*[

*L*],

*a*,

_{1}*a*, and

_{2}*a*are the cost coefficients for well drilling, well and pump installation, and well operation, respectively,

_{3}*Q*is the water demand [

_{d}*L*],

^{3}T^{−1}*β*is a constant,

*N*is the total number of potential well locations, and

*b*is a binary variable indicating which well is to be drilled.

_{k}*h*as the state variable of the problem is dependent on the pumping discharge via the equation governing the steady-state saturated flow in unconfined aquifers, which is a non-linear, second-order partial differential equation known as the Boussinesq equation (Wu & Zhu 2006): with where

_{k}*K*is hydraulic conductivity [

_{h}*LT*],

^{−2}*R*is the surface recharge rate [

*L T*,

^{−1}]*x*and

*y*are the independent Cartesian coordinates [

*L*],

*δ(x-x*is the Dirac delta function evaluated at (

_{k},y-y_{k})*x*), and

_{k},y_{k}*Δx*and

*Δy*are the grid spacing used to discretize the governing equation [

*L*]. The optimization problem defined by Equations (1)–(8) is a MINLP problem which can theoretically be solved using any of the conventional or modern optimization algorithms, with corresponding deficiencies and merits. While conventional gradient based methods are known to be efficient in solving convex non-linear problems, they often fail to find the global optimum for non-convex problems. A common way out of this difficulty is to solve the problem using different initial guesses at the expense of higher computational cost. The inefficiency of the gradient-based methods is resonated when they are used to solve MINLP problems, such as the one considered in this paper, in which a branch and bound method is often used to deal with the integer variables. Modern heuristic methods are shown to be effective in avoiding local optimums while dealing with mixed variables, but they also suffer from high computational cost. This is due to the fact that solution of such a problem using modern heuristic algorithms, such as GA, requires the use of a simulation-optimization approach in which an optimization method is combined with a simulation tool. In this approach, the optimization module creates a trial solution to the problem while the simulation module evaluates the performance of each solution. In the simulation module, the hydraulic head is found for each set of decision variables (i.e. well numbers, positions, depth and discharge) by solving Equation (8) using any of the available numerical methods such as FDM, finite element method (FEM) or BEM. This requires that for every trial solution created by the optimization algorithm, a system of non-linear algebraic equations resulting from the discretization of Equation (8) is solved for the hydraulic head, requiring high computational effort in particular when large-scale problems are considered.

#### Proposed hybrid mixed-integer linear programming-linear programming-linear programming methodology

*h*as its variable as follows: Next, the governing Equation (9) is discretized using a central differencing for the second order derivatives of Equation (9), leading to the following system of equations on the finite differences (FD) grid points (

^{2}*i,j*) defined over the 2D aquifer domain: Equations (10) and (11) are now considered as a set of constraints on the hydraulic head and the well discharges. A set of binary variables

*b*,

_{k}*k*

*=*

*1, 2,…N*corresponding to each of the wells is now defined and used to reformulate the constraints (4) and (7) in the following form: The problem under consideration can now be defined as follows.

Find the binary variables *b _{k}*, the well depth

*d*, the pumping rates

_{k}*Q*and the hydraulic head

_{k},*h*as the decision variables of the problem to minimize the objective function of Equation (1) subject to Equations (2), (3), (5), (6), (7), (10), (11), and (12). The MINLP problem defined by these equations is reformulated in the form of a hybrid MILP-LP-LP problem and solved in a three-step manner to exploit the efficiency and effectiveness of the LP solvers.

_{k},#### Step 1

*h*and wells' depths

_{k}*d*, find the pumping rates

_{k}*Q*and the set of binary variables

_{k}*b*

_{k}, k*=*

*1,…N*, as the decision variables of the MILP sub-problem with the following objective function: subject to: where

*X*and

_{i,j}*Y*are two sets of non-negative slack variables representing the possible violations of an arbitrary solution from the constraints (17) and (18), and

_{i,j}*α*is a large enough penalty factor to ensure that any infeasible solution represented by non-zero values of

*X*and

_{i,j}*Y*has a total cost greater than any feasible solution. The problem defined by Equations (13)–(18), is clearly a MILP problem which can be efficiently solved using branch and bound and simplex methods. The slack variables

_{i,j}*X*and

_{i,j}*Y*are introduced in the MILP stage to avoid failure of the simplex method at the early stages of the computation where a feasible solution to the corresponding MILP problem could not be found. The first step problem so defined will ensure a feasible final solution (Khadem & Afshar 2013). Having solved the above problem for the pumping rates

_{i,j}*Q*and a set of binary variables

_{k}*b*

_{k}, k*=*

*1,…N*, a second sub-problem is defined and solved for the hydraulic head

*h*as follows.

_{i,j}#### Step 2

*Q*and

_{k}*b*from the first step and values of

_{k}*d*and

_{k}*h*assumed in the previous iteration, find the hydraulic head

_{i,j}^{*}*h*

_{i,j}, i*=*

*1,…I,*and

*j*

*=*

*1,…J*as the solution of the following LP sub-problem: subject to: where is the initial hydraulic head in the aquifer. The solution

*h*obtained in the second step is then used to define and solve the third sub-problem.

_{i,j}#### Step 3

*d*using the known values of

_{k}*Q*,

_{k}*b*from the first step, known values of

_{k}*h*from the second step, and assumed values of

_{i,j}*d*from the previous iteration as follows: subject to: The solution

_{k}^{*}*d*obtained in the third step along with the

_{k}*Q*,

_{k}*b*and

_{k}*h*from previous steps is then used to iterate the three-step procedure until convergence is achieved.

_{i,j}### Transient flow condition

*t*. The problem defined by Equations (1), (2), (3), (5), (6), (7), (10), (11), and (12) can now be mathematically redefined with the following cost function: subject to: where

*t*indicates the period of the operation (any

*t*superscript denotes the value for the

*t*th period),

*T*is the total operation periods,

*Δt*is the time period duration, and

*S*is the dimensionless aquifer specific yield coefficient.

_{y}The time domain for this problem is discretized into several time sub-domains each with the steady-state flow condition explained previously. Therefore, the transient problem is treated as multiple linearized problems for the steady-state condition and solved by the proposed three-step MILP-LP-LP methodology.

## TEST EXAMPLE

### Problem one: steady-state flow condition

For the first problem the water demand is taken to be equal to 30,000 m^{3}day^{−1} and it is assumed that all the pumping rates must operate in the range of 0 to 7,000 m^{3}day^{−1}. Plus, the hydraulic head values at well locations must be non-negative. The cost coefficients are set equal to *a _{1}* = 4,221 $m

^{−1},

*a*= 0.12 $m

_{2}^{−4}

*, a*= 0.03 $m

_{3}^{−4}, and

*β*= 0.299, which are adopted from McKinney & Lin (1994). This problem was previously solved using non-linear programming (NLP) and GA by McKinney & Lin (1994), GA and simulated annealing (SA) by Wang & Zheng (1998), shuffled complex evolution-University of Arizona (SCE-UA) by Wu

*et al.*(1999), and HS by Tamer Ayvaz (2009). In this paper, this problem is solved using the proposed MILP-LP-LP approach.

Table 1 summarizes the results of the hybrid MILP-LP-LP model as well as the results of other studies. As can be seen, not only are the pumping rates calculated by the hybrid MILP-LP-LP model rather better than those obtained by the other solution methods, but also the total supply is in good agreement with the demand (which doesn't occur for the NLP model). It is clearly seen that the hybrid MILP-LP-LP model obtains a slower objective function value (106,107$ compared to that of NLP, GA, SA, SCE-UA, and HS with 106,990$, 107,134$, 106,992$, 106,896$, 106,896$, respectively). This fact may be due to the evolutionary nature of applied heuristic methods, which have to be coupled with a simulation package (consider the non-linear nature of the simulation process of unconfined aquifers) and tune many of the model parameters, all leading to attaining a near optimal solution. In addition, NLP models, as mentioned before, are not efficient for non-convex problems, are based on local gradients, and require good initial guesses to find the optimal solutions. Table 2 compares the NLP model and the proposed MILP-LP-LP model in iterations used to achieve the optimal solution and the corresponding CPU time, showing the great superiority of the proposed method. Results indicate a significant reduction in computation cost in comparison to those of the NLP method while obtaining a better solution (considering the cost minimization objective function, the value of the objective function obtained by the proposed method is smaller than that of NLP), emphasizing the efficiency and effectiveness of the proposed hybrid MILP-LP-LP methodology. It should be noted that the problem was solved 10 times using different initial random guesses expressing the robustness of the method.

Well | NLP* | GA* | GA** | SA** | SCE-UA*** | HS**** | Proposed MILP-LP-LP method |
---|---|---|---|---|---|---|---|

1 | 5,600 | 6,000 | 6,000 | 6,000 | 5,886 | 6,147 | 5,614 |

2 | 4,550 | 5,000 | 4,000 | 4,000 | 4,833 | 4,620 | 0 |

3 | 5,860 | 7,000 | 6,000 | 6,000 | 5,894 | 6,142 | 5,614 |

4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

5 | 0 | 0 | 0 | 0 | 0 | 0 | 4,770 |

6 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

7 | 7,000 | 6,000 | 7,000 | 7,000 | 6,694 | 6,142 | 7,000 |

8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

10 | 7,000 | 6,000 | 7,000 | 7,000 | 6,693 | 6,543 | 7,000 |

Total pumping | 30,010 | 30,000 | 30,000 | 30,000 | 30,000 | 30,000 | 30,000 |

Total costs | 106,990 | 107,134 | 106,992 | 106,992 | 106,896 | 106,891 | 106,107 |

Well | NLP* | GA* | GA** | SA** | SCE-UA*** | HS**** | Proposed MILP-LP-LP method |
---|---|---|---|---|---|---|---|

1 | 5,600 | 6,000 | 6,000 | 6,000 | 5,886 | 6,147 | 5,614 |

2 | 4,550 | 5,000 | 4,000 | 4,000 | 4,833 | 4,620 | 0 |

3 | 5,860 | 7,000 | 6,000 | 6,000 | 5,894 | 6,142 | 5,614 |

4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

5 | 0 | 0 | 0 | 0 | 0 | 0 | 4,770 |

6 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

7 | 7,000 | 6,000 | 7,000 | 7,000 | 6,694 | 6,142 | 7,000 |

8 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |

10 | 7,000 | 6,000 | 7,000 | 7,000 | 6,693 | 6,543 | 7,000 |

Total pumping | 30,010 | 30,000 | 30,000 | 30,000 | 30,000 | 30,000 | 30,000 |

Total costs | 106,990 | 107,134 | 106,992 | 106,992 | 106,896 | 106,891 | 106,107 |

Model | Total solver iterations | CPU time (s) |
---|---|---|

NLP | 1,222,752 | 7,080 |

Proposed MILP-LP-LP | 13,457 | 14 |

Model | Total solver iterations | CPU time (s) |
---|---|---|

NLP | 1,222,752 | 7,080 |

Proposed MILP-LP-LP | 13,457 | 14 |

### Problem two: transient flow condition

^{−1}with specific yield equal to 0.1. The distance between the ground surface and the aquifer bottom is 150 m everywhere. The initial hydraulic head is assumed to be 100 m at all grid points. Downstream demand for period one to four is

*Q*= 130,000 m

_{d}^{1}^{3}day

^{−1},

*Q*= 145,000 m

_{d}^{2}^{3}day

^{−1},

*Q*= 150,000 m

_{d}^{3}^{3}day

^{−1}, and

*Q*= 130,000 m

_{d}^{4}^{3}day

^{−1}. The lower and upper limits on the pumping rates are 0 and 30,000 m

^{3}day

^{−1}, respectively. The cost coefficients of the problem formulations are

*a*

_{1}*=*

*a*

_{2}*=*

*a*= 1 $ m

_{3}^{−4}. This problem was also previously solved using different methods such as differential dynamic programming (DDP) by Jones

*et al.*(1987), GA and SA by Wang & Zheng (1998), and HS by Tamer Ayvaz (2009). Note that the non-linear problem defined by Equations (1)–(7) was also solved by the authors for comparison purposes, the results of which are shown in Table 3.

Well | DDP^{+} | GA^{++} | SA^{++} | HS^{+++} | NLP | Proposed MILP-LP-LP method |
---|---|---|---|---|---|---|

1 | 30,000 | 28,000 | 30,000 | 29,476 | 30,000 | 30,000 |

2 | 30,000 | 28,000 | 30,000 | 29,458 | 30,000 | 30,000 |

3 | 21,924 | 28,000 | 17,000 | 18,600 | 24,736 | 22,934 |

4 | 21,924 | 4,000 | 17,000 | 20,597 | 30,000 | 22,967 |

5 | 7,494 | 12,000 | 16,000 | 10,764 | 7,410 | 7,855 |

6 | 7,494 | 12,000 | 11,000 | 12,801 | 7,501 | 7,985 |

7 | 5,582 | 14,000 | 7,000 | 2,452 | 2,812 | 4,128 |

8 | 5,582 | 4,000 | 2,000 | 5,853 | 2,781 | 4,129 |

Demand | 130,000 | 130,000 | 130,000 | 130,001 | 135,240 | 129,998 |

1 | 28,000 | 30,000 | 29,919 | 30,000 | 30,000 | |

2 | 28,000 | 30,000 | 29,721 | 30,000 | 30,000 | |

3 | 10,000 | 21,000 | 25,125 | 30,000 | 28,265 | |

4 | 20,000 | 21,000 | 19,792 | 30,000 | 28,303 | |

5 | 8,000 | 10,000 | 11,482 | 9,040 | 9,423 | |

6 | 14,000 | 12,000 | 13,757 | 9,153 | 9,572 | |

7 | 16,000 | 12,000 | 7,865 | 3,421 | 4,719 | |

8 | 20,000 | 9,000 | 7,339 | 3,385 | 4,719 | |

Demand | 144,000 | 145,000 | 145,000 | 144,999 | 145,001 | |

1 | 28,000 | 30,000 | 29,795 | 30,000 | 30,000 | |

2 | 30,000 | 30,000 | 29,599 | 30,000 | 30,000 | |

3 | 12,000 | 25,000 | 24,188 | 30,000 | 30,000 | |

4 | 28,000 | 18,000 | 27,917 | 30,000 | 30,000 | |

5 | 6,000 | 12,000 | 12,160 | 10,924 | 9,991 | |

6 | 8,000 | 15,000 | 12,947 | 11,058 | 10,152 | |

7 | 26,000 | 10,000 | 3,299 | 4,026 | 4,930 | |

8 | 12,000 | 10,000 | 10,096 | 3,990 | 4,925 | |

Demand | 150,000 | 150,000 | 150,001 | 149,998 | 149,998 | |

1 | 28,000 | 30,000 | 29,876 | 30,000 | 30,000 | |

2 | 28,000 | 30,000 | 29,976 | 30,000 | 30,000 | |

3 | 20,000 | 22,000 | 20,200 | 24,736 | 22,976 | |

4 | 22,000 | 12,000 | 22,185 | 24,759 | 22,967 | |

5 | 4,000 | 5,000 | 10,902 | 7,410 | 7,855 | |

6 | 8,000 | 13,000 | 10,506 | 7,501 | 7,985 | |

7 | 6,000 | 6,000 | 160 | 2,812 | 4,128 | |

8 | 14,000 | 12,000 | 6,196 | 2,781 | 4,129 | |

Demand | 130,000 | 130,000 | 130,001 | 129,999 | 130,040 | |

Total costs | 28,693,336 | 29,779,432 | 29,572,110 | 29,540,860 | 29,564,447 | 28,992,860 |

Well | DDP^{+} | GA^{++} | SA^{++} | HS^{+++} | NLP | Proposed MILP-LP-LP method |
---|---|---|---|---|---|---|

1 | 30,000 | 28,000 | 30,000 | 29,476 | 30,000 | 30,000 |

2 | 30,000 | 28,000 | 30,000 | 29,458 | 30,000 | 30,000 |

3 | 21,924 | 28,000 | 17,000 | 18,600 | 24,736 | 22,934 |

4 | 21,924 | 4,000 | 17,000 | 20,597 | 30,000 | 22,967 |

5 | 7,494 | 12,000 | 16,000 | 10,764 | 7,410 | 7,855 |

6 | 7,494 | 12,000 | 11,000 | 12,801 | 7,501 | 7,985 |

7 | 5,582 | 14,000 | 7,000 | 2,452 | 2,812 | 4,128 |

8 | 5,582 | 4,000 | 2,000 | 5,853 | 2,781 | 4,129 |

Demand | 130,000 | 130,000 | 130,000 | 130,001 | 135,240 | 129,998 |

1 | 28,000 | 30,000 | 29,919 | 30,000 | 30,000 | |

2 | 28,000 | 30,000 | 29,721 | 30,000 | 30,000 | |

3 | 10,000 | 21,000 | 25,125 | 30,000 | 28,265 | |

4 | 20,000 | 21,000 | 19,792 | 30,000 | 28,303 | |

5 | 8,000 | 10,000 | 11,482 | 9,040 | 9,423 | |

6 | 14,000 | 12,000 | 13,757 | 9,153 | 9,572 | |

7 | 16,000 | 12,000 | 7,865 | 3,421 | 4,719 | |

8 | 20,000 | 9,000 | 7,339 | 3,385 | 4,719 | |

Demand | 144,000 | 145,000 | 145,000 | 144,999 | 145,001 | |

1 | 28,000 | 30,000 | 29,795 | 30,000 | 30,000 | |

2 | 30,000 | 30,000 | 29,599 | 30,000 | 30,000 | |

3 | 12,000 | 25,000 | 24,188 | 30,000 | 30,000 | |

4 | 28,000 | 18,000 | 27,917 | 30,000 | 30,000 | |

5 | 6,000 | 12,000 | 12,160 | 10,924 | 9,991 | |

6 | 8,000 | 15,000 | 12,947 | 11,058 | 10,152 | |

7 | 26,000 | 10,000 | 3,299 | 4,026 | 4,930 | |

8 | 12,000 | 10,000 | 10,096 | 3,990 | 4,925 | |

Demand | 150,000 | 150,000 | 150,001 | 149,998 | 149,998 | |

1 | 28,000 | 30,000 | 29,876 | 30,000 | 30,000 | |

2 | 28,000 | 30,000 | 29,976 | 30,000 | 30,000 | |

3 | 20,000 | 22,000 | 20,200 | 24,736 | 22,976 | |

4 | 22,000 | 12,000 | 22,185 | 24,759 | 22,967 | |

5 | 4,000 | 5,000 | 10,902 | 7,410 | 7,855 | |

6 | 8,000 | 13,000 | 10,506 | 7,501 | 7,985 | |

7 | 6,000 | 6,000 | 160 | 2,812 | 4,128 | |

8 | 14,000 | 12,000 | 6,196 | 2,781 | 4,129 | |

Demand | 130,000 | 130,000 | 130,001 | 129,999 | 130,040 | |

Total costs | 28,693,336 | 29,779,432 | 29,572,110 | 29,540,860 | 29,564,447 | 28,992,860 |

The results shown in Table 3 clearly indicate that the solution obtained using the proposed MILP-LP-LP method has a sound agreement with water demand for every period while attaining a better cost solution. Yet, the DDP model reflects lower costs. Computational burden is a concern for researchers, especially for large-scale real-life management problems. Among all available models with similar accuracy and abilities, the faster the model the more applicable it is and the wider use it gets. Modern heuristic search methods are known to be time-consuming as a result of the iterative use of simulation models for evaluating every single function. In contrast, conventional mathematical methods often bear fast convergence rates while being reported to have their own defects. To prove the claim, the second problem was solved using the NLP model, which shows a reduction of approximately 99% in computational time can be achieved by implementing the proposed MILP-LP-LP model indicated in Table 4.

Model | Total solver iterations | CPU time (s) |
---|---|---|

NLP | 41,356,543 | 175,600 |

Proposed MILP-LP-LP | 54,390 | 248 |

Model | Total solver iterations | CPU time (s) |
---|---|---|

NLP | 41,356,543 | 175,600 |

Proposed MILP-LP-LP | 54,390 | 248 |

As can be seen from the results, the proposed model optimally solved the transient condition problem without any remarkable constraint violations, with a significant reduction in computational time expressing an efficient and effective method. It should again be noted that the problem was solved 10 times using different initial random guesses, emphasizing the robustness of the method.

## DISCUSSION AND CONCLUSION

A hybrid MILP-LP-LP method is proposed in this paper for simultaneous optimal design and operation of the groundwater utilization systems which play a significant role in life nowadays. For this, a discretized version of the differential equation governing the aquifer is first embedded into the model formulation as additional constraints. The resulting MINLP is then decomposed into three sub-problems with different sets of decision variables. Having started with random values for any of the decision variables, the three sub-problems are solved iteratively until convergence is achieved. In order to test the performance of the proposed management model, a management problem from the literature is solved for both steady-state and transient flow conditions and the results are presented and compared with those obtained using other methods. Solutions to the problem indicate that hybrid MILP-LP-LP gives a rather better objective function value than other methods without any violation in constraints. Plus, the proposed method requires less iterations and CPU time. The results demonstrate the efficiency, effectiveness, and robustness of the proposed hybrid MILP-LP-LP methodology.