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

• k

well's position

•
• i

FD grid point over X axis

•
• j

FD grid point over Y axis

•
• t

operation period

### Variables

• Qk

pumping rate of the well at position k

•
• hk

•
• bk

binary variables that indicates drilling and operation of well at position k

•
• Xi,j, Yi,j

non-negative slack variables used to avoid infeasibilities

•
• hi,j

hydraulic head at node (i,j) of FD grid

### Parameters

• N

total number of potential well locations

•
• a1

cost coefficient for drilling

•
• a2

cost coefficient for well and pump installation

•
• a3

cost coefficient for operation

•
• dk

well's depth at position k

•
• β

a constant

•
• Qd

total demanded pumping rate

•
• hmin

minimum allowed hydraulic head due to environmental issues

•
• Qmin

minimum pumping rate limit according to the available pumps

•
• Qmax

maximum pumping rate limit according to the available pumps

•
• dmax

maximum allowed well depth

•
• m

total number of wells required to be drilled and operated

•
• Kh

hydraulic conductivity of the aquifer

•
• R

surface recharge rate

•
• x, y

independent Cartesian coordinates

•
• xk, yk

Cartesian coordinates of the well at position k

•
• Δx, Δy

FD grid spacing

•
• δ

Dirac delta function

•
• α

penalty factor

•
• initial hydraulic head in the aquifer

•
• Sy

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.

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

Find the number and position of the wells to be drilled out of a set of known potential well locations, and their corresponding depth and operational policy represented by their pumping rate to minimize the following cost function (Tamer Ayvaz 2009):
1
subject to the following constructional and operational constraints.
Total discharge from the aquifer must exceed a given demand.
2
Hydraulic heads at the well locations must be greater than a specified lower bound.
3
The pumping rates from each operational well must be in the range of specified lower and upper bounds to account for the capacity of the available pumps in the market and to avoid dewatering problems.
4
The well depths should not be greater than an upper bound.
5
And finally, the total number of wells to be drilled must not exceed a specified number.
6
where the binary variable is defined as:
7
In the above equations, Qk is the pumping rate of the well at position k [L3 T−1], dk is the depth of the well at location k [L], hk is the hydraulic head at the well at location k [L], a1, a2, and a3 are the cost coefficients for well drilling, well and pump installation, and well operation, respectively, Qd is the water demand [L3 T−1], β is a constant, N is the total number of potential well locations, and bk is a binary variable indicating which well is to be drilled.
The hydraulic head hk 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):
8
with
where Kh is hydraulic conductivity [LT−2], R is the surface recharge rate [L T−1], x and y are the independent Cartesian coordinates [L], δ(x-xk,y-yk) is the Dirac delta function evaluated at (xk,yk), and Δ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

The first step in developing the proposed hybrid mixed-integer linear programming-linear programming-linear programming (MILP-LP-LP) methodology is to embed the differential Equation (8) as a new set of constraints into the problem formulation. For this, first the governing equation is converted to a linearized version of the Boussinesq equation for steady-state saturated unconfined flow by rewriting it in terms of h2 as its variable as follows:
9
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 (i,j) defined over the 2D aquifer domain:
10
11
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 bk, 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:
12
The problem under consideration can now be defined as follows.

Find the binary variables bk, the well depth dk, the pumping rates Qk, and the hydraulic head hk, 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.

#### Step 1

Assuming known values for the hydraulic heads hk and wells' depths dk, find the pumping rates Qk and the set of binary variables bk, k=1,…N, as the decision variables of the MILP sub-problem with the following objective function:
13
subject to:
14
15
16
17
18
where Xi,j and Yi,j are two sets of non-negative slack variables representing the possible violations of an arbitrary solution from the constraints (17) and (18), and α is a large enough penalty factor to ensure that any infeasible solution represented by non-zero values of Xi,j and Yi,j 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 Xi,j and Yi,j 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 Qk and a set of binary variables bk, k=1,…N, a second sub-problem is defined and solved for the hydraulic head hi,j as follows.

#### Step 2

Using known values of Qk and bk from the first step and values of dk and hi,j* assumed in the previous iteration, find the hydraulic head hi,j, i=1,…I, and j=1,…J as the solution of the following LP sub-problem:
19
subject to:
20
21
22
where is the initial hydraulic head in the aquifer. The solution hi,j obtained in the second step is then used to define and solve the third sub-problem.

#### Step 3

The third sub-problem is defined and solved for the wells' depths dk using the known values of Qk, bk from the first step, known values of hi,j from the second step, and assumed values of dk* from the previous iteration as follows:
23
subject to:
24
The solution dk obtained in the third step along with the Qk, bk and hi,j from previous steps is then used to iterate the three-step procedure until convergence is achieved.

### Transient flow condition

Considering the problem described previously for multiple periods of operation results in a transient condition, and introduces a new variable to the problem formulation, namely 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:
25
subject to:
26
27
28
29
30
31
32
33
where t indicates the period of the operation (any t superscript denotes the value for the tth period), T is the total operation periods, Δt is the time period duration, and Sy is the dimensionless aquifer specific yield coefficient.

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

The test example used to challenge the capability of the proposed method is the optimal operation of the utilization system of the 2D unconfined aquifer presented in Figure 1. The objective of both problems (i.e. steady-state and transient conditions) is to minimize the total pumping cost while satisfying the water demand, which is supplied from 10 potential well locations given in Figure 1.
Figure 1

Unconfined aquifer model for the first problem: (a) plan view, (b) elevation view (McKinney & Lin 1994).

Figure 1

Unconfined aquifer model for the first problem: (a) plan view, (b) elevation view (McKinney & Lin 1994).

### Problem one: steady-state flow condition

For the first problem the water demand is taken to be equal to 30,000 m3day−1 and it is assumed that all the pumping rates must operate in the range of 0 to 7,000 m3day−1. Plus, the hydraulic head values at well locations must be non-negative. The cost coefficients are set equal to a1 = 4,221 $m−1, a2 = 0.12$m−4, a3 = 0.03 $m−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. Table 1 Problem one: comparison of hybrid MILP-LP-LP and different solution methods (units: m3 day−1 pumping rates and$ for costs)

WellNLP*GA*GA**SA**SCE-UA***HS****Proposed MILP-LP-LP method
5,600 6,000 6,000 6,000 5,886 6,147 5,614
4,550 5,000 4,000 4,000 4,833 4,620
5,860 7,000 6,000 6,000 5,894 6,142 5,614
4,770
7,000 6,000 7,000 7,000 6,694 6,142 7,000
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
WellNLP*GA*GA**SA**SCE-UA***HS****Proposed MILP-LP-LP method
5,600 6,000 6,000 6,000 5,886 6,147 5,614
4,550 5,000 4,000 4,000 4,833 4,620
5,860 7,000 6,000 6,000 5,894 6,142 5,614
4,770
7,000 6,000 7,000 7,000 6,694 6,142 7,000
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
Table 2

Problem one: comparison of hybrid MILP-LP-LP and NLP solution in total solver iterations and CPU time

ModelTotal solver iterationsCPU time (s)
NLP 1,222,752 7,080
Proposed MILP-LP-LP 13,457 14
ModelTotal solver iterationsCPU time (s)
NLP 1,222,752 7,080
Proposed MILP-LP-LP 13,457 14

### Problem two: transient flow condition

The second problem considers a multi-period problem with four periods in a year (i.e. each period consists of 91.25 days) and the aquifer specifications demonstrated in Figure 2. The problem was solved using the same assumption as the first problem except for the following: the hydraulic conductivity of the aquifer media is set to 86.4 m day−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 Qd1 = 130,000 m3 day−1, Qd2 = 145,000 m3 day−1, Qd3 = 150,000 m3 day−1, and Qd4 = 130,000 m3 day−1. The lower and upper limits on the pumping rates are 0 and 30,000 m3 day−1, respectively. The cost coefficients of the problem formulations are a1=a2=a3 = 1 $m−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. Table 3 Problem two: comparison of hybrid MILP-LP-LP and different solution methods in every management period (units: m3 day−1 for pumping rates and$ for costs)

WellDDP+GA++SA++HS+++NLPProposed MILP-LP-LP method
30,000 28,000 30,000 29,476 30,000 30,000
30,000 28,000 30,000 29,458 30,000 30,000
21,924 28,000 17,000 18,600 24,736 22,934
21,924 4,000 17,000 20,597 30,000 22,967
7,494 12,000 16,000 10,764 7,410 7,855
7,494 12,000 11,000 12,801 7,501 7,985
5,582 14,000 7,000 2,452 2,812 4,128
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
28,000 30,000 29,919 30,000 30,000
28,000 30,000 29,721 30,000 30,000
10,000 21,000 25,125 30,000 28,265
20,000 21,000 19,792 30,000 28,303
8,000 10,000 11,482 9,040 9,423
14,000 12,000 13,757 9,153 9,572
16,000 12,000 7,865 3,421 4,719
20,000 9,000 7,339 3,385 4,719
Demand   144,000 145,000 145,000 144,999 145,001
28,000 30,000 29,795 30,000 30,000
30,000 30,000 29,599 30,000 30,000
12,000 25,000 24,188 30,000 30,000
28,000 18,000 27,917 30,000 30,000
6,000 12,000 12,160 10,924 9,991
8,000 15,000 12,947 11,058 10,152
26,000 10,000 3,299 4,026 4,930
12,000 10,000 10,096 3,990 4,925
Demand   150,000 150,000 150,001 149,998 149,998
28,000 30,000 29,876 30,000 30,000
28,000 30,000 29,976 30,000 30,000
20,000 22,000 20,200 24,736 22,976
22,000 12,000 22,185 24,759 22,967
4,000 5,000 10,902 7,410 7,855
8,000 13,000 10,506 7,501 7,985
6,000 6,000 160 2,812 4,128
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
WellDDP+GA++SA++HS+++NLPProposed MILP-LP-LP method
30,000 28,000 30,000 29,476 30,000 30,000
30,000 28,000 30,000 29,458 30,000 30,000
21,924 28,000 17,000 18,600 24,736 22,934
21,924 4,000 17,000 20,597 30,000 22,967
7,494 12,000 16,000 10,764 7,410 7,855
7,494 12,000 11,000 12,801 7,501 7,985
5,582 14,000 7,000 2,452 2,812 4,128
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
28,000 30,000 29,919 30,000 30,000
28,000 30,000 29,721 30,000 30,000
10,000 21,000 25,125 30,000 28,265
20,000 21,000 19,792 30,000 28,303
8,000 10,000 11,482 9,040 9,423
14,000 12,000 13,757 9,153 9,572
16,000 12,000 7,865 3,421 4,719
20,000 9,000 7,339 3,385 4,719
Demand   144,000 145,000 145,000 144,999 145,001
28,000 30,000 29,795 30,000 30,000
30,000 30,000 29,599 30,000 30,000
12,000 25,000 24,188 30,000 30,000
28,000 18,000 27,917 30,000 30,000
6,000 12,000 12,160 10,924 9,991
8,000 15,000 12,947 11,058 10,152
26,000 10,000 3,299 4,026 4,930
12,000 10,000 10,096 3,990 4,925
Demand   150,000 150,000 150,001 149,998 149,998
28,000 30,000 29,876 30,000 30,000
28,000 30,000 29,976 30,000 30,000
20,000 22,000 20,200 24,736 22,976
22,000 12,000 22,185 24,759 22,967
4,000 5,000 10,902 7,410 7,855
8,000 13,000 10,506 7,501 7,985
6,000 6,000 160 2,812 4,128
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
Figure 2

Plan view of unconfined aquifer model for the second problem (Wang & Zhang 1998).

Figure 2

Plan view of unconfined aquifer model for the second problem (Wang & Zhang 1998).

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.

Table 4

Problem two: comparison of hybrid MILP-LP-LP and NLP solution in total solver iterations and CPU time

ModelTotal solver iterationsCPU time (s)
NLP 41,356,543 175,600
Proposed MILP-LP-LP 54,390 248
ModelTotal solver iterationsCPU 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.

## REFERENCES

REFERENCES
Asher
M. J.
Croke
B. F. W.
Jakeman
A. J.
Peeters
L. J. M.
2015
.
Water Resour. Res.
51
,
5957
5973
.
Bostan
M.
Afshar
M. H.
M.
2015
.
Eng. Optimiz.
47
,
550
560
.
Gaur
S.
Chahar
B. R.
Graillot
D.
2011
.
J. Hydrol.
402
,
217
227
.
Jones
L.
Willis
R.
Yeh
W. W. G.
1987
.
Water Resour. Res.
23
,
2097
2106
.
Katsifarakis
K. L.
2008
.
Water Resour. Manag.
22
,
1089
1099
.
Katsifarakis
K. L.
Karpouzos
D. K.
1998
Minimization of pumping cost in zoned aquifer by means of genetic algorithms
. In:
Proceedings of the International Conference on Protection and Restoration of the Environment IV
, pp.
61
68
.
Katsifarakis
K. L.
Karpouzos
D. K.
Theodossiou
N.
1999
.
Eng. Anal. Bound. Elem.
23
,
555
565
.
M.
Afshar
M. H.
2013
.
Irrig. Drain.
62
,
120
128
.
McKinney
D. C.
Lin
M. D.
1994
.
Water Resour. Res.
30
,
1897
1906
.
Rosenwald
G. W.
Green
D. W.
1974
A method for determining the optimum location of wells in a reservoir using mixed-integer programming
.
Soc. Res. Eng. J.
14
,
44
54
.
Sidiropoulos
E.
Tolikas
P.
2004
.
Water Air Soil Poll.
4
,
227
239
.
Singh
A.
2014
.
Sci. Total Environ.
499
,
414
423
.
Tamer Ayvaz
M.
2009
.
32
,
916
924
.
Tamer Ayvaz
M.
Karahan
H.
2008
.
J. Hydrol.
357
,
76
92
.
Wang
M.
Zheng
C.
1998
.
J. Am. Water Resour. Assoc.
34
,
519
530
.
Wu
J.
Zhu
X.
2006
.
Frontiers of WWW Research and Development – APWeb 2006.
3841
,
986
995
.
Wu
J.
Zhu
X.
Liu
J.
1999
.
Sci. China Ser. E: Technol. Sci.
42
,
521
529
.