Application of Optimal Control and Fuzzy Theory for Dynamic Groundwater Remediation Design

DOI 10.1007/s11269-008-9293-1

Application of Optimal Control and Fuzzy Theory

for Dynamic Groundwater Remediation Design

Received: 1 May 2007 / Accepted: 9 June 2008 / Published online: 10 July 2008

© Springer Science + Business Media B.V. 2008

Abstract Obtaining optimal solutions for time-varying groundwater remediation design is a challenging task. A novel procedure first employs input/output data sets obtained by constrained differential dynamic programming (CDDP). Then the Adaptive-Network-Based Fuzzy Inference System (ANFIS), which is a fuzzy inference system (FIS) implemented in the adaptive network framework, is applied to acquire time-varying pumping rates. Results demonstrate that the FIS is an efficient way of groundwater remediation design.

Keywords CDDP· ANFIS · Remediation design · Ground water

1 Introduction

The pump-and-treat (P&T) method is one of the most common groundwater reme-diation methods. The method is primarily useful for decontaminating groundwater with highly soluble pollutants by pumping out contaminated groundwater and treat-ing the water. The feasibility of coupltreat-ing optimization techniques with groundwater flow and transporting simulation to design P&T systems has been extensively studied (Chang and Shoemaker1992; Culver and Shoemaker1993; McKinney and Lin1994; Wang and Zheng1998; Chang and Hsiao2002; Chu et al.2005; Chang et al.2007). Containing the contamination by removing contaminated ground water should be a dynamic process. Dynamic policies which allow changing pumping policies as the

H.-J. Chu· L.-C. Chang (

_B

Department of Civil Engineering, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu, Taiwan 30050, Republic of China e-mail: [email protected]

H.-J. Chu

contaminant plume moves, would expectedly be more cost-effective than steady state policies. Chang and Shoemaker (1992) employ an optimal control method, called the Successive Approximation Linear Quadratic Regulator (SALQR), to design a time-varying pumping system for contaminated aquifer remediation. Culver and Shoemaker (1992) find that varying policies are more cost-effective than time-invariant policies. Culver and Shoemaker (1993) extend the SALQR method by adding second derivatives governed by a quasi-Newtonian approach (QNDDP), accelerating algorithm convergence. Chang and Hsiao (2002) integrate Constrained Differential Dynamic Programming (CDDP) and Genetic Algorithm (GA) to op-timize total remediation cost, while attempting to resolve the planning problem of simultaneously considering fixed costs of well installation and operating costs of time-varying pumping rates. The CDDP used herein is a modification of SALQR. The CDDP algorithm has been shown efficient in solving time-varying problems. However, considering a more realistic problem and applying CDDP to large-scale field sites requires more computational time. The computational work for the CDDP model is proportional to the O(n3_{), where n is the total number of state variables} (Chang and Shoemaker1992).

Recent artificial intelligence (AI) aiding the process control system has been adopted widely in remediation design and water resources management (Rogers and Dowla 1994; Geng et al. 2001; Chang and Chang 2001; Hu et al.2003; Rao et al. 2003,2005; Spiliotis and Tsakiris 2007). Fuzzy logic and artificial neural network (ANN) are common tools in developing an intelligent control-engineering model. Such intelligent control is a technology resembling the human thinking process in decision making and learning. Geng et al. (2001) presents a fuzzy expert system for managing petroleum contaminated sites and uses AI techniques to construct a support tool for site remediation decision-making. Hu et al. (2003) develops a fuzzy control system controlling pumping rates based on the measured pollution level. They combine a fuzzy set theory and control technology for process control of a complex environmental system and apply the system to a real-world case study in western Canada. Rao et al. (2005) uses ANN and simulated annealing (SA) for planning groundwater development in coastal deltas and uses a trained ANN as the SEAWAT model to predict final groundwater concentration under variable pumping conditions. Jang (1993) combines ANN and fuzzy logic as ANFIS: Adaptive-Network-Based Fuzzy Inference System. The system is a neuro-adaptive learning technique which provides a method for the fuzzy modeling procedure to learn information about a data set. Chang and Chang (2001) apply the ANFIS and GA to water resources management and use the GA to search the optimal reservoir operation based on a given inflow series. The ANFIS is then built to create the fuzzy inference system (FIS), to determine optimal water release according to reservoir depth and inflow.

This study applies the ANFIS and CDDP approaches into remediation design. The CDDP finds optimal pumping rates and generates ANFIS training patterns. The ANFIS then trains the FIS to make decisions concerning how to pump out contaminants according to the groundwater head and concentration.

2 Formulation of the Management Model

The management model attempts to minimize the total cost of remediation, com-posed of pumping and treatment system operating costs. The problem can be formulated as min ui t, i∈I,t=1,...,N J= i∈I N  t=1  a1uit+ a2uit  Li_∗− hit+1  (1) subject to {xt+1} = T (xt, ut, t) , t = 1, 2, . . . , N (2) cN, j≤ cmax, j ∈  (3)  i∈I uit≤ utotal, t = 1, 2, . . . , N, i ∈ I (4) uimin≤ u i t≤ u i max, t = 1, 2, . . . , N, i ∈ I (5)

where I is a network design. The upper index i denotes a well in the network design (I). J(·) represents total cost of I; xt=

 ht: ct

T

∈ R(nh+nc)×1_{are the state continuous}

variables representing heads (ht) and concentrations (ct), nh and nc denote total number of hydraulic heads and concentrations, respectively; ut∈ Rm×1 represent the vector of control variables whose dimension depends on I, m is the number of control variables; T(xt, ut, t) represents the transition equation; Φ is the set of observation wells; N is the total number of time steps; a1 and a2 are factors used to convert treatment cost, and operating cost, respectively, into monetary values ($); L∗∈ Rm×1 denotes the distance from the ground surface to the lower datum of the well aquifer; ht+1denotes the hydraulic head for nodes at time step t + 1; utotal represents the maximum allowable total pumping rates from all extraction wells. The component in Eq. 1 is the operating cost, involving pumping and treatment costs. The transition equation, T, in Eq.2is solved with ISOQUAD (Pinder1978), a finite element of groundwater flow and the transport model for a confined two-dimensional aquifer. The numerical solution is obtained by applying the Galerkin finite element method for space derivative and an implicit finite difference scheme for time derivative. The transport model includes changes in head due to pumping and changes in contaminant concentration owing to advection, diffusion, dispersion, and linear equilibrium sorption. The set of observation wells (Φ) in practical applications is the group of sites that meet the water quality standard (Eq.3). The constraint in Eq.4specifies capacity constraint for treatment plants. Equation5specifies capacity constraints for each well. The groundwater remediation model defined by Eqs.1to 5is a time-varying optimization problem.

3 The Methodology

3.1 Constrained Differential Dynamic Programming (CDDP)

Nearly global optimization techniques, such as simulated annealing (Rao et al.2005), genetic algorithm (Huang and Mayer1997), or tabu search (Zheng and Wang1999), efficiently solve complex groundwater management problems. However, applying these techniques to solve time-varying policies dramatically increases computational resources required (Culver and Shoemaker 1997). The CDDP exceeds conven-tional dynamic programming (DP) and the algorithms solve time-varying policies efficiency. Accordingly, CDDP overcomes the “curse of dimensionality”, a serious limitation for conventional techniques. The CDDP is a successive approximation technique for solving optimal control problems, iteratively determining the optimal solution to the problem stated in (1) to (5). The objective function becomes a function of control and state variables with identical time index (t).The algorithm is employed in the backward and forward sweep to resolve the series of quadratic problems. State variables in the backward sweep are considered as unknown and op-timal control laws, which are unknown variables include derivatives of the objective function and the transition equation. The algorithm in the forward sweep computes updated policy using feedback function and performs ISOQUAD to determine new system status. Murray and Yakowitz (1979), Jones et al. (1987), and Chang and Shoemaker (1992) provide a detailed discussion of the CDDP algorithm.

3.2 Artificial Neural Network (ANN)

The ANN attempts simulation of the brain (Biological neural network). The network consists of an interconnected group of artificial neurons and processes information using a connectionist approach to computation. A number of neurons are arranged in an input layer, one or more hidden layers, and an output layer. An ANN has two neural processing phases. Learning is a process of adapting connection weights in an ANN to produce desired output. The recalling process attempts to retrieve information based on weights obtained from the learning process, and predict output data of the new example. The back-propagation learning algorithm is among the well-known algorithms used to train an ANN. This algorithm is applied to the multi-layer feed-forward network consisting of processing neurons.

3.3 Fuzzy Inference System (FIS)

Fuzzy set and fuzzy logic extend upon traditional Boolean logic and deal with imprecision in human experience. Zadeh (1965) introduced the fuzzy set theory for a mathematical description of imprecision and uncertainty.

A FIS is a mathematical model based on a fuzzy rule system. A fuzzy rule system is defined as the set of rules consisting of a set of input variables or premises in the form of fuzzy sets with membership function, and a set of consequence also in the form of a fuzzy set. Utilizing linguistic variables, fuzzy rules, and fuzzy reasoning provides a

tool to incorporate human expert experience. A fuzzy inference system is composed of five functional blocks.

1. A rule base containing a number of fuzzy if–then rules;

2. A database, which defines membership functions of the fuzzy set used in fuzzy rules;

3. A fuzzification interface, which transforms crisp inputs into degrees of match with linguistic values;

4. A decision-making unit, which performs inference operations on the rules; 5. A defuzzication interface, which transforms fuzzy results of the interface into

crisp values.

3.4 Adaptive Network-based Fuzzy Inference System (ANFIS)

This paper presents the architecture and learning procedure underlying ANFIS. ANFIS is a neuro-adaptive learning method for the fuzzy modeling procedure to learn information. ANFIS uses a hybrid learning algorithm to identify parameters of Sugeno-type fuzzy inference systems. The system applies a combination of the least-squares method and the backpropagation gradient descent method for training FIS membership function parameters to emulate a given training data set. By using a hybrid learning procedure, the proposed ANFIS constructs input–output mapping based on both human knowledge (in the form of fuzzy if–then rules) and stipulated input–output data pairs. The simplified ANFIS architecture is shown as Fig. 1,

∏ ∏ ∏ ∏ N N N N

including two inputs x, y and one output z by a first-order polynomial. Two rules are expressed as:

Rule 1: If x is A1, and y is B1, then z = f1= p1x+ q1y+ r1 (6)

Rule 2: If x is A1, and y is B2, then z = f2= p2x+ q2y+ r2 (7)

Rule 3: If x is A2, and y is B1, then z = f3= p3x+ q3y+ r3 (8)

Rule 4: If x is A2, and y is B2, then z = f4= p4x+ q4y+ r4 (9) Layer 1. Input nodes:

Every node in the layer performs a membership function

O1,i= μAi(x) , i = 1, 2 (10)

O1,i= μBi−2(y) , i = 3, 4 (11)

Where x and y are crisp inputs to node i and Ai, Bi are the linguistic labels char-acterized by appropriate membership functions. Usually the bell-shaped function is chosen as a membership function given by:

μAi = 1 1+x−ci ai  2bi andμBi−2= 1 1+x−ci ai  2bi

Where {ai, bi, ci} is the parameter set of the membership function in the premise part of fuzzy if then rule.

Layer 2. Rule nodes,:

After generating membership functions, layer links correspond to preconditions in fuzzy logic rule. The T-norm performs incoming signal multiplications to generate layer outputs. Therefore, layer outputs O2,kare products of corresponding degrees from layer 1.

O2,k= wk= μAi(x) μBj(y) , k = 1, . . . 4, i = 1, 2, j = 1, 2 (12) Layer 3. Normalized nodes, N:

The node output in the layer is the ratio of node output from the previous layer to the total. O3,i= wk= w i 4  i=1wk , i = 1, . . . .4 (13)

Layer 4. Consequest nodes:

In the layer, the node corresponds to the following function O4,i= wifi= wi(pix+

qiy+ ri), i = 1,. . . .,4

Where { pi, qi, ri} is the parameter set in the consequent part of the first-order Sugeno fuzzy model.

Layer 5. Output nodes:

Defuzzification inference transforms the model fuzzy results into crisp output, therefore, the sum of all incoming signals generates the decision crisp output.

O5,1= 4  i=1 wifi= 4  i=1 wifi 4  i=1wi (14)

Fig. 2 Flowchart of the model

for remediation design

N

3.5 Study Procedure

Numerical simulation for a large field-scale problem of groundwater remediation takes a great deal of computational time, as hundreds to thousands of simulation runs may be required during the search for optimal pumping strategy. Figure2further clarifies the following step-by-step description to solve the problem.

Step 1: Create a knowledge base

The CDDP is a successive approximation technique for solving optimal control problems, iteratively determining the optimal solution to the problem stated in (1) to (5) and creating a knowledge base of different training data patterns. The CDDP generates the data for determining optimal pumping rates for each well.

Step 2: Train FIS by ANFIS

The ANFIS applies a combination of the least squares method and the back prop-agation gradient descent method for training FIS membership function parameters to emulate a given training data set. The calculated output values are compared to example target value, and the weights are updated according to minimal difference between calculated values and the target on the next iteration. The ANFIS attempts to retrieve information on parameters determined from the training process. Step 3: Apply FIS to determine optimal pumping rates

After membership function parameters are adjusted, FIS determines pumping rates (ut) based on head (ht) and concentration (ct) at time step t.

41 40 38 36 37 39 42 (x,y)=(0,0) 16 9 2 1 8 15 3 10 17 30 23 22 29 24 31 y = 600 m h_b 4 11 18 5 12 19 25 32 26 33 21 14 7 6 13 20 35 28 27 34 60-90 ppm 120-180 ppm 79 no-flow boundary 59 57 58 51 44 43 50 45 52 65 72 64 71 66 73 x = 1200 m 86 85 78 87 80 61 60 46 53 47 54 67 74 68 75 no-flow boundary 62 63 49 56 48 55 70 77 69 76 hb 88 81 89 82 91 84 90 83 90-120 ppm 30-60 ppm 0-30 ppm =20m =20m

Table 1 Aquifer properties of

example application Coefficient Value

Hydraulic conductivity 4.31_{× 10}−4m/s Diffusion coefficient 1× 10−7m2_/s

Storage coefficient 0.001

Porosity 0.2

Sorption partitioning coefficient 0.245 cm3_/g

Media bulk density 2.12 g/cm3

Aquifer thickness 10 m

4 Numerical Results

This study presents solutions obtained for a hypothetical, isotropic confined aquifer with dimensions 600 m× 1,200 m to demonstrate algorithm performance described above. Figure3indicates the finite-element mesh which has 91 finite-element nodes. The north and south sides are no-flow boundaries. The west and the east constant-head boundary are 20 m. The initial peak concentration within the aquifer is approximately 150 mg/L, and the water quality goal at the end of five years must be less than 0.5 mg/L (cmax) at all observation wells. The planning horizon in the model is divided into twenty time steps over five years and each time step in the management model (t) is 91.25 days. Tables1and2list aquifer properties, the cost coefficient values and constraint values.

4.1 Case 1

This work generates the training data for ANFIS using the CDDP. Different initial concentration levels of pollution within the aquifer ranging from 120 mg/L to 180 mg/L create a knowledge base of 200 data patterns. Using a given input/output data set, the ANFIS constructs a FIS with tuned membership function parameters. Case1 discusses input selection based on the monitoring network alternative and alternative states. The alterative states (hydraulic head and concentration or only concentration) are selected as inputs of each well.

Table3lists the number of rules and parameters in ANFIS. For example, Case 1.1 chooses two monitoring wells with two states (hydraulic head and concentration) and four input pairs, as well as sixteen rules (24_{) and eight membership functions.} Twenty-four premise parameters (eight membership functions with three parameters in each) are added to the 80 consequent parameters (16(4 + 1)) for a total of 104 parameters. Fuzzy inference maps from a given input to an output using fuzzy logic after adjusting parameters using ANFIS. The following the statistical parameter ‘maximum relative

Table 2 Cost coefficient and

constraint values Coefficient Value

a1 $40, 000 m3 s· m · t a2 $1, 000  m3 s· t L_∗ 120 m utotal 2,000 L/s umax 120 L/s umin 0 L/s

Table 3 Case 1 setting and results The number o f T he number o f T he number o f T he number T he number T he number o f M aximum relative C PU time (s) observation w ells states inputs (observations) o f sets o f rules parameters e rror (%) Case 1.1 2 2 4 2 1 6 104 3.50 19.14 Case 1.2 1 2 2 2 4 24 3.77 1.45 Case 1.3 2 1 2 2 4 24 4.48 1.70 Case 1.4 1 1 1 2 2 10 6.80 1.31

26 28 30 32 34 36 Time step T o

tal pumping/injection rate(L/s)

FIS CDDP 26 28 30 32 34 36 Time step T o ta l pum pi ng/ injec tion r ate (L /s ) FIS CDDP 26 28 30 32 34 36 1 3 5 7 9 11 13 15 17 19 1 3 5 7 9 11 13 15 17 19 Time step T

otal pumping/injection rate(L/s)

26 28 30 32 34 36 1 3 5 7 9 11 13 15 17 19 1 3 5 7 9 11 13 15 17 19 Time step T

otal pumping/injection rate(L/s)

FIS CDDP FIS CDDP

(a)

(b)

(d)

(c)

Fig. 4 Optimal pumping rates at each period in Case 1. a Optimal pumping rates at each period in

Case 1.1. b Optimal pumping rates at each period in Case 1.2. c Optimal pumping rates at each period in Case 1.3. d Optimal pumping rates at each period in Case 1.4

error (MRE)’ quantifies error to show predictive accuracy. The optimal pumping rate (output) during t is estimated by FIS, and the target is the optimal pumping rate determined by CDDP. MRE= max _|T t− Ot| Tt  (15) where Tt: Target value at time step t; Ot: Output value at time step t.

Table 4 Case 2 setting and results

The number The number The number of Maximum relative CPU time(s) of sets of rules parameters error (%)

Case 2.1 2 2 10 6.80 1.31 Case 2.2 3 3 15 5.20 1.38 Case 2.3 4 4 20 4.63 1.48 Case 2.4 5 5 25 4.57 1.98 Case 2.5 6 6 30 4.47 2.25 Case 2.6 7 7 35 4.00 2.35 Case 2.7 8 8 40 3.93 2.47 Case 2.8 9 9 45 1.27 3.01

Fig. 5 Normalize CPU time

and MRE under different set numbers (Case 2) 0 0. 2 0. 4 0. 6 0. 8 1 0 2 4 6 8 10

The number of sets

Normalize MRE / CPU time

MRE CPU time

Figure 4shows optimal pumping rates determined by CDDP and estimated by FIS. Table3represents the MRE for all cases. Results show that a prediction is more accurate with more input pairs. Findings show that the number of inputs and rules significantly influence prediction accuracy.

Table3also represents CPU times for all cases. That is, CPU time increases with increased input numbers. This study finds that an increase in inputs decreases relative error and increases CPU time. The average computational time using CDDP is 31.3 s in Case 1. However, the average computational time using ANFIS reduces over 90% by Case 1.2.

4.2 Case 2

Case 2 explores the relationship between prediction accuracy and parameter num-bers in fuzzy inference. Case 2 is modified from Case 1.4 with more parameters used in fuzzy inference. In Table4, the number of sets ranges from two to nine and the number of parameters ranges from ten to 45. Fuzzy inference maps from a given input (a state) to an output using fuzzy logic. Figure5shows the relationship between normalized MRE/normalized CPU time and fuzzy set numbers. The figure shows that prediction error decreases with increasing numbers of fuzzy sets. However, CPU time increases with fuzzy set numbers. Figure 5also shows that normalized MRE and normalized CPU time is the same with five fuzzy sets. Therefore, this study assumes five fuzzy sets as adequate set numbers for balancing the trade-off between error and CPU time. Figure6shows optimal pumping rates computed by CDDP and FIS at each period for Case 2.4. Table3indicates that Case 1.4 has only

Fig. 6 Optimal pumping rates

at each period in Case 2.4

26 28 30 32 34 36 1 3 5 7 9 11 13 15 17 19 Time step T

otal pumping/injection rate(L/s)

one observation and thus has the largest maximum relative error (6.8%). Table 4 on the other hand, shows that the maximum relative error of Case 2.4 (4.57%) is less than that of Case 1.4 although they have the same basic setup. Compared with Figs. 4 and 6, the table indicates the same result. Therefore, increasing FIS parameters increases FIS performance.

5 Conclusions

This study extends artificial intelligence in the time-dependent remediation design. Searching for conventional optimal solutions of the nonlinear and time-variant system requires considerable computation. The current work uses the CDDP to generate solution pairs, and then uses ANFIS to adjust membership function pa-rameters to solve the problem. The FIS determines time-varying optimal pumping rates quickly and accurately after identifying membership function parameters. One advantage to using FIS is implicitly considering practical situation uncertainties. In this study, the FIS gains a controller obtained by CDDP based on different initial plume conditions. The FIS application for controlling pumping rates greatly simplifies complex dynamics, and relative error is small between CDDP and FIS. The cases also reveal that the number of observations and fuzzy sets significantly effect prediction performance. This paper serves as a reference to deal with future works.

References

Chang LC, Chang FJ (2001) Intelligent control for modelling of real-time reservoir operation. Hydrol Process 15:1621–1634

Chang LC, Hsiao CT (2002) Dynamic optimal groundwater remediation including fixed and opera-tion costs. Ground Water 40(5):481–490

Chang LC, Shoemaker CA (1992) Optimal time-varying pumping for groundwater remediation application of a constrained optimal control algorithm. Water Resour Res 28(12):3157–3173 Chang LC, Chu HJ, Hsiao CT (2007) Optimal planning of a dynamic pump–treat–inject groundwater

remediation system. J Hydrol 342(3–4):295–304

Chu HJ, Hsiao CT, Chang LC (2005) Optimal remediation design in groundwater systems by intel-ligent techniques. Knowl-Based Intelintel-ligent Information and Engineering Syst, Pt 2, Proceedings 3682:628–634

Culver TB, Shoemaker CA (1992) Dynamic optimal control for groundwater remediation with flexible management periods. Water Resour Res 28(3):629–641

Culver TB, Shoemaker CA (1993) Optimal control for groundwater remediation by differential dynamic programming with quasi-Newton approximations. Water Resour Res 29(4):823–831 Culver TB, Shoemaker CA (1997) Dynamic optimal ground-water reclamation with treatment

capi-tal costs. J Water Resour Plan Manage ASCE 123(1):23–29

Geng JQ, Chen Z, Chan CW, Huang GH (2001) An intelligent decision support system for manage-ment of petroleum-contaminated sites. Expert Syst Appl 20:251–260

Hu ZY, Huang GH, Chan CW (2003) A fuzzy process controller for in situ groundwater bioremedi-ation. Eng Appl Artif Intell 16:131–147

Huang C, Mayer AS (1997) Pump-and-treat optimization using well locations and pumping rates as decision variables. Water Resour Res 33(5):1001–1012

Jang J-SR (1993) ANFIS: adaptive-network-based fuzzy inference systems. IEEE Trans Syst Man Cybern 23(3):665–685

Jones L, Willis R, Yeh WW-G (1987) Optimal control of nonlinear groundwater hydraulics using differential dynamic programming. Water Resour Res 23(11):2097–2106

McKinney DC, Lin MD (1994) Genetic algorithm solution of groundwater-management models. Water Resour Res 30:1897–1906

Murray DM, Yakowitz SJ (1979) Constrained differential dynamic programming and its application to multireservoir control. Water Resour Res 15(5):1017–1027

Pinder GF (1978) Galerkin finite element models for aquifer simulation. Rep. 78-WR-5, Dept. of Civ. Eng., Princeton Univ., Princeton, NJ

Rao SVN, Thandaveswara BS, Bhallamudi SM et al. (2003) Optimal groundwater management in deltaic regions using simulated annealing and neural networks. Water Resour Manag 17(6): 409–428

Rao SVN, Bhallamudi SM, Thandaveswara BS, Sreenivasulu V (2005) Planning groundwater development in coastal deltas with paleo channels. Water Resour Manag 19:625–639

Rogers LL, Dowla FU (1994) Optimization of groundwater remediation using artificial neural networks with parallel solute transport modeling. Water Resour Res 30:457–481

Spiliotis M, Tsakiris G (2007) Minimum cost irrigation network design using interactive fuzzy integer programming. J Irrigat Drain Eng-ASCE 133(3):242–248

Wang M, Zheng C (1998) Groundwater management optimization using genetic algorithms and simulated annealing: formulation and comparison. J Am Water Resour Assoc 34(3):519–530 Zadeh LA (1965) Fuzzy sets. Inf Control 8(3):338–353

Zheng CM, Wang PP (1999) An integrated global and local optimization approach for remediation system design. Water Resour Res 35:137–148