Phys. Chem. Emrh (B), Vol. 24, No. 1-2, pp. 139-144,
Pergamon
1999
0 1998 Elsevier Science Ltd All rights reserved 1464-1909/99/$--see front matter
PII: S 1464- 1909(98)00025-2
Optimization Models in Groundwater Management, Based on Linear and Mixed Integer Programming. An Application to a Greek Hydrogeological Basin A. A. Psilovikos Department Received
of Rural Engineering, IS January
1998;
Aristotle
accepted
University
54622, Thessaloniki,
Greece
6 July 1998
The paper compares two optimization Abstract. methods used in groundwater management, based on linear as well as on mixed integer programming (Lp and MIP). The solution obtained of the combined simulation - optimimtion models, consists of three steps (Psilovikos, 1996) : 1Simulation model - MODPLOW , 2Managementmodel MODMAN, 3.Optimization model - LINDO An hydrogeological basin in Northern Greece, was usedasacasestudyinthisproject.ThecollectedQtaare based on 26 managed wells, for a period of 12 months (Psilovikos et al, 1996). In order to reduce the amount of calculations, only 4 managed wells were selected for a period of 3 months. The results obtained from the solution with the two methods above, are : a) both models satisfy the same piezometric and balance constraints, b) the MIP model is more complicated and the feasible region of solutions is more constrained than the LP one, because a number of integer cons&aims is included, so the LP model can be considered as the relaxed version of the m One. 0 1998 Elsevier Science Ltd. All rights
of Thessaloniki,
(months) as a comparison between the two methods of mathematical programming.
2 MATHEMATICAL MODEL The mathematicalmodel consists of three steps, as described below : 2.1 Simulation Model : MODFLOW (MC Donald aad Harbagh, 1988). The governing equation is the basic three dimensional differential equation of groundwater movement, and given from the formula below (equation 1)
(1.1 where are values of hydraulic k-q Kw, Kzz conductivity along the x, y, z coordinate axes which is parallel to the- major axes of hydraulic assumedtobe conductivity [Lr ] h the potensiometric head [L] the volumettic flux per unit volume. Represents W sources and/or sinks of water [T’] the speci& storage of the porous material [L”] S, t timem. Afterspacediscretizationinx,y,zdirections,asshownin figure 1, the above equation 1 is transformed into a tinite difference form in three dimensions :
reserved.
1 Introduction The development of optimization models in groundwater management is examinedinthisprojec&basedonlinearas well as on mined - integer programming methods (LP and MIP). The solution of optimization models requims the combined use of simulation - optimii&on models that are linked together via the management model. The solution of thiskindofmodels,consistsofthreesteps: (MODFLOW, 3 - D finite 1. simulationmodel difTereuce method in groundwater movement). 2. Management Model (MODMAN, response matrix method, superp&tion in space and time). 3. Opbmization Problem (UNDO, linear and mixed integer pro-g method, objective-function and UXK&iints). AcatchmentareainNorthernGreecewasusedas casestudyinthisproject.~ecollecteddatawas~on26 managed wells, for a period of 12 months. This is a complicated case and an extensive response matrix is obtained from the management model. In order to simplify this case only 4 managed wells were selected for 3 periods
C%.j+t / zk
* bi.j-1.k.
- h.j,km)+
CCi-t,2,j,k
’ (hi-l.j.km
-$.j,km)cCCi+if2,j.k
‘~+l.j,km
-hi.j,km)+
- hi.j.km)’
’ &.j.k+lm
-h,.j,km)+
cvi.j.k-I,
q.jJ
‘4.j.k”
1
(h,j,k-I”
+Qi.j.k
=
%,j,k
CVi,j.k+l,Z
’
(b, j+l.km
CRi,j-1,l.k
-
h&km)+
hP~.j.t-~~~~ .Arj .Ac, .Av, ’ tm - tw
where :
(2)
the piaomemc LIP, b,;, h1amj2”, lx+,j;, hij*~“,~,,“,~~ levels in cells (ij-l,k), (ij+l,k), (i-l&k), (i+l,j,k), (i&k-l), (i,j,k+l), (i&k) respectively, at time step m [L] 139
A. A. Psilovikos:
140
Optimization
Models in Groundwater
the piezometric level in (ij,k) cell, at time step m-l [L] CR,.UU. CRJm, , C&?.j,k 9 C(h2J.k 3 cvu,k-IL? , cxj,k+,, Hydraulic conductance coefficients. For example CRstnl is the conductance in row i and layer k between nodes (ij-l,k) the conductance in column j and layer k and (ij,k) , cci-l&j,k behen nodes (i-l j,k) and (ij,k), etc. [L”T’] Coefficients corresponding on cell (ij,k) PUl..Q&j, p’T’], and &?‘I, respectively specific storage in cell (i.j.k) [L”] s (Uk V = Arj Ac, &k Volume of (ij,k) dell p3] time at end of time step m and m- 1 tm. b-1 respectively PI time derivative in piezometic level in cell t, - tnl-1 (i.j.k). obtained with backward difference approximation. This is a fully implicit scheme (backward difference form) and the slice successive overrelaxation method (SSOR) is use-dfor the solution of equation 2. 2.2 Management Model : MODMAN (MODflow MANagement, Gnxnwald, 1994) The large-scale problems in water management. have forced researchers to develop management models. Such a model is MODMAN and operates as a linkage from the simulation proc&re that takes place with MODFLOW model to the optimizationpmceduxthattakesplacewiththeLINDO model. The opimization pmcedme is the third step of the mathematical model. The MODMAN model is based on the response matrix method under the assumpdon of space supexposition for the steady state problems (equation3.) and both space and time superposition for unsteady state - transient problems (equation4.). So the p!yical @lem has transformedinto a mathematicalone and it ISexpressedas below : ?? Steady State (Space Supexposition) Hi =Ui +xv.Qj il ?? Unsteady State (Space (Psilovilcos, 1996) T
(3.) and
Time
Superposition).
N
k=l J=I
where : unmanagedhead at control point i at the end of the U, last managingperiodT. managedhead at control point i at the end of the last H,T F!PperiodT average drawdownin each i observationwell at the endoftheTpumpingperiodduetoauuitrateofpum@ngat the j managed well applied throughoutthe k pumping period (Colarullo, S. M. Heidari,T. MaddockIII, 19&1) Pumpingrateat well j duringthe k pumping period. QJk An example of 3 managed wells in 2 pumping periods (months), gives the response matrixbelow. where Ahh,T is the drawdownat pumping well i at the end of the last managingpe&dT:
(1)
Ally -
al.1
Ahh:” Lg’ @2’
Management (1) Q1.Z
(1) Q 2.1 (1)
at.2
(I)
0 a;‘;
a2,3
a3.1
a::
4:
a;:
a;;)
Ah?’
i2) a 2.1
a 2.1
Ah?'
a3.1
a3,2
(b
0
0 a;:’
(1)
a,‘;’
=
0 0
aI,)
(1)
h a 2,3
(b
0
0
Q;”
0
Q:”
0 a,‘;;
Qi” ’
(b
Q 2.1
Q2.2
a:f:
a3.1
a3,2
a3.3
(2) (2) (2) (1) 0) a3,3
0) _
Q:”
Q:" ,Q3
(2)
(5.) 2.3 Optimiiatioa
Model
: LINDO (S&rage
1990,
Winston 1996) 2.3.1 Linear Programming Method The last part of the mathematical model is the optimization procedure with the LlNDO model. The scope is the maximization of total pimping rates from all the wells while all the constraints must be satisfied. So the target is the optbizabon of a linear objective function of the form : Z=$“C:.Q; k--l j=l
(6.1
A number of linear constraints in piezometric level in specific control points in the aquifer, has to be satisfied based on the response matrix method (Fsilovikos, 1996). In this case the control points are the pumping wells :
Balance Cons&&s referred to the entire quantity of extracting water for each managing peri* that have to be equal to the water demand for irrigation.
&$Q:=d
(8.)
Finally, con&aints in minimum and maximum pumping rates are formulated : Oq++a.. (9.1 where i = 1,. ..,4 control point in piezometric level. j = 1,...,4 pumpingwell. k = 1,....3 managingperiod. maximummanageddrawdownat control point i hT at the end of the last irrigationperiodT. minimum allowable head at control point i at the Hi,tiT end of the last irrigationperiod T. The terms Cj” aIe Cost coefficients. In the present case the! areequaltoccl))becauSeweareinterestedonthequamitative management of water and not the total cost managementof punping. 2.3.2 Mixed Integer Programming. MANagement,Greenwald 1994)
(MODflO\\
This method is more or less the same with the linear programmingone. Moreoveran extra constraintto choose K active wells among J probablewell positions has to do with two types of integez constraints.
A. A. Psilovikos: Optimization Models in Groundwater Management ??
Well on/off constraints. and
??
Integer variable Summation constraints The on/off constraints are constructed with binary variables which are integer variables that can only have a value of 0 or 1. The on/off constmint for a well forces the binary variable to a value of 1 if the well is on. or 0 if the well is off. The integer variable summation constmint, based on the biuaty variables, enforces the limit on the number of active wells allowed. The form of the on/off constmint is : ??
. EXTBACTION Q, +M,I, 20
(Negative well rate)
INJECTION
(Positive well rate)
??
where rate at well j (Negative for pumping) 2 a large positive number with an absolute value greater than that of the largest well rate, and a binary variable acting as on/& switch for well j Ij The form of the integer variable summation constraints is
where K is the mtmber of requimd wells probable well sites J binary variables. If pumping take place from a well, I, the corresponded binary variable is equal to 1, if not it is equaltoo.
APPLICATION
The a@ication of the combined Simulation - Optimizatioa model. was carried out in the Eidomeni - Enones freatic aquifer of a 12 square kilometers basin in northern Greece. The Axios river flows through this area from north to south. The data used for the construction of the simulation model was precipitation, specific irrigation recharge. pumping rates in 26 wells, average head in the river, specific storage. hydraulic conductivitv, ground levels, well base level. The areaisdividedintosmallsquareareas(cellsof200m*200 m), where the mean depth of the aquifer is 18 m. The time increment of At was chosen to be equal to 1 month. The management period has taken 12 months (figure 2). The present example considers 4 managed wells for a period of 3 months. The target is the optimization of the aquifer based on linear and mixed integer programming, that maximization of pumping yield while hastodowiththe remaming within prescribed drawdown, balance, minim and maximum pumping rates. and moreover. integer constraints. The equations has the form below : 0 Objective function maxf(X,=&
,Q;.
(12)
k=l ,=I 0
Piezometxic constraints
H;,,,,” > 40, Hi,,
Balanceconstraints
$Q: =9300, k = 1,2,3 ??
(1%
Well on/off Constraints
Q,+I,.MrO, Q,+I,.Ms:O.
(16)
Integer variable summation constraints I, +I? +I, +I, $3,
??
I, =O.orl,
(17)
j = I,....4
Theresultsareshowninfigures 1,2and3atulalso in tables 1 and 2 at the end of the article. The comparison of the two methods is discussed below.
(11)
rl
3
.
Q,+I,.MzO (11)
j=l,..., J
z-+000 (14)
Qz+Iz.MrO. (10)
Q,-M,I,
iI,5K.
Minimum and maximum pmping rates O>Q,~-18OO,O~Q,~-25OO.O~Q, z-3500.02Q,
141
r 39, Hi,,, t 37, Hi.,,
2 36
(13)
4
RESULTS - CONCLUSSIONS
For the management of the present case. hvo different optimization models were consuuckd The first was based on LP and the second on MIP. The conclusions that can obtained from the comparison of the two models, are : The MIP model is more complicated and the feasible region of solution is more constmined than the feasible region of the LP model. This is due to the well on/off constraints and also to the integer variable summation constmint. so we can consider the LP solution to be the relaxed version of the MIP solution. The soIution of the MIP model is accomplished in 25 successive itenttive steps while the LP model solution is achieved in 4 steps. In the MIP model there is an integer variable summation constmint ((3 out of 4, wells and the pumping is applied inmmtber1.3and4wells.whilethenumbw2wellisnot operated (Table 2). The optimal sohnion in the LP model. is to pump at smaller rates from all wells (Table 1). Both the models satisfy the same balance cons&aims and the whole quantity of water being extracted is the same for both cases. equal to 27900 [m%Iay]. only the distribution of the pumping water is different (Table 1). The managed heads obtained from the two-optimization pro&ures are different to each other because of the different distribution of pumping recharge (Figure 3). Generally the LP model is better to be used when a uniform piezometric level and pumping recharge distribution is reqired Gn the other hand the MIP model is preferable. when from a field of wells we have to choose only K pumping wells from a set of J. The (J-K) wells are not going to be operated but the rates from the J operating wells will be quite increased So in this case the
A. A. Psilovikos: Optimization Models in Groundwater Management
142
piezometric constraints have to be more flexible otherwise the problemwill be infeasible. Further research will include optimization models applied in large scale catchmentareas based on LP and MIP models, and comparisonwith other methods of mathematical programming as Non-linear (NLP) and Quadratic programming(QP). 5 .
REFERENCES Bear. I., 1979. Hydraulics of Groundwater. HillBook
.
.
GIW
-
Colarullo. S. M. Heidari. T. Maddock III, 1984, Identification of an optimal groundwater management strategy in a contaminated aquifer, ~rarer ResowcesBulhn.
.
MC
Company,London.
vol. 20. no. 5, pp.747-76-O.
MC Donald M. G., and A.W. Harbaugh, 1988. A modular three - dimensional finite - difference groundwater flow model U.S.Geological Suryo~. Techmgves ofWater ResourcesInves~gations,Book6 ChapterAl. Winston W.,l994. Operations Research - Applications and ~gOritbRIS, Third Edition, Lhubury Press, Belmont,
Table 1.
Optimum PmmphS Rate8 fL?‘]
California
.
Gwnwak!, R M., 1994. MODflow MANagement : An Optimhatioa Module for MODPLOW, IGKVC, version
.
Aris. 19%. Optimum management in aquifer studies using the Linear Programming (LP) method. An application to Eidomeni - Enones ma. MScThew, Thessalonrkr.pp. IS5 Gorelick, S. M.. 1983. A Review of Distributed Parameter Groundwater Management Modelling Methods, UTR.vol.19. so. 2. pp. 305 - 319. Psilovikos A., C. Moutsopoulos, C. Tzimopoulos. S. Giamtopoulos,1996. Water mass balance estimation in the aquifer of Eidomeni - Evzanes using the Modflow model, Proc. r’ CCVII m Wafer Resources Management, Lorisa.
3.02.
.
.
F%ilovikos
I Linear ProgrammingSolution I2& I Period.5 / I 1” I 3ti I Wells 1 -GZU 0 1800 1800 2 -G~z,ss 1800 0 2500 3 - GE,_E2 3500 3500 1000 4 -Gzas 4000 4000 4000 Mixed - Integer Programming
Greece. pp. 210-222.
.
Psilovikos A.. and C. Tzimopoulos. 1998. Pumping cost analysis in groundwater management, using tbe MODMAN (MODflow MANagemeat) model. under publication Proc. IP Internatmnal Con/: On Computahona~Method3 In WaterResources.Crete. Greece. 1998
Table 2.
Rlnuy Vuiabks of the MIP aohtlon
A. A. Psilovikos:
I
I
Optimization
I
Models in Groundwater
Management
I
143
I
2500.00~
2000.007 ! \ , 1500.00_:
ooo.oo-
1
500.00~ ! I
1 I 1000.00 15Ob.00 2000.00 25Ob.00 3000.00 35Ob.00
144
A. A. Psilovikos: Optimization Models in Groundwater Management
41
40
39 E 1
‘*
g 1
ST 36
35
34 3
2 well Id
Flyre
3. Managed heads in well points obtaIned from the two models
4