Environ. Sci. Techno/. 1995, 29, 1145-1155
Optimal Field-Scale Groundwater Remediation Using Neural Networks and theGenetic Algorithm LEAH L. ROGERS,+ FARID U. DOWLA,t AND VIRGINIA M. JOHNSON*,* Earth Sciences Division and Environmental Restoration Division, L-206, Lawrence Livermore National Laboratory, Livermore, California 94551
The escalating costs of environmental cleanup together with the conflicting concerns of various stakeholders motivate the search for improved management methodologies to reduce costs. The primary computational obstacle to groundwater management using formal optimization techniques is the flow and transport numerical model. Presented here is an approach to nonlinear optimization that harnesses the strengths of artificial neural networks and the genetic algorithm to search for more cost-effective remediation strategies within practical time frames. Results drawn from a major western Superfund site demonstrate the potential for saving tens of millions of dollars.
Introduction The enormous expense associated with the nation's environmentalproblems is now receiving considerable public attention (1-3). While much of this debate focuses on the need to change environmentallaws to better reflect actual risks ( 1 , 3 ) ,such changes may not be realized for years. In the meantime, there are techniques drawn from the fields of contaminant transport modeling and optimizationthat can be used to help reduce cleanup costs. Contaminant transport models have a vital role to play in environmental remediationplanning. Whatever their current limitations, models are the primary tools planners can employ to forecast the likely impact of their proposed cleanup strategies over time. If it can be granted then that models will continueto be developed and used for decision-making, questions arise as to how best to exploit the effort that goes into the construction and calibration of a model for a particular site. Combiningthe methodologies of contaminant transport modeling and formal optimization allows planners to evaluate extremely large numbers of treatment alternatives. These exhaustive searches can yield more efficient and costeffective cleanup designs than are found through informal design strategies alone (4). In the groundwater management area, both linear and nonlinear optimization methods have been applied to various problems in capture and cleanup (5-10). In general, hydraulic control of a contaminant plume usually requires only linear optimization methods. But when contaminant concentrations are incorporated into the problem, nonlinear methods are often required. The nonlinearity stems from the multiplication of unknown contaminant concentrations by groundwater velocities and dispersion coefficients, elements which in turn are dependent on the decision variables (e.g., well locations and pumping rates) being optimized. The computer hardware and software required to implement nonlinear optimization techniques on groundwater problems involving more than a handful of wells quickly tax the resources of most organizations, mainly because of the time required by the model to evaluate each potential treatment alternative. In this paper, we present an approach to the nonlinear optimization of groundwater remediation that is designed to make field-scaleapplications practical using conventional computer equipment. First, artificial neural networks (A"s) are trained to predict selected outcomesof a groundwater contaminanttransport model. Then, a genetic algorithm (GA) searches through millions of possible treatment alternatives, rapidly evaluating the effectiveness of each alternative with predictions generated by the A " s . Traditional approaches to optimization rely on search proceduresthat use the model itself to generate these predictions. In the current approach, A " s take the place of this time-consumingmodeling step duringthe search, greatly reducingthe computationaltime needed to complete the optimization. The approach is demonstrated on a contaminant transport model and * Corresponding author; e-mail address:
[email protected]; Fax: (510) 422-3118. + Earth Sciences Division. Environmental Restoration Division.
*
0013-936X'95/0929-1145$09.00/0
@ 1995 American Chemical Society
VOL. 29, NO. 5 , 1995 I ENVIRONMENTAL SCIENCE &TECHNOLOGY rn 1145
B
Number of wells
I
Regulatory constraint Proportion of contaminant removed
I
Cost reduction
a
S
13
Met
Met
Mat
0.97
0.98
1.05
74%
69%
associated data drawn from an actual groundwater remediation problem.
Remediation Example Figure la shows the distribution of total volatile organic compounds (VOCs) in the groundwater at a well-charac1146. ENVIRONMENTAL SCIENCE
&TECHNOLOGV/VOL.29. NO.
5.1995
66%
terized Superfund site. The primary components are TCE and PCE from solvents originally dumped onto the ground and which haveleachedintothe groundwater over aperiod of 50 years. The site is also contaminated by gasoline constituents and chromium. This problem considers only the remediation of VOCs. A pump-and-treat strategy was
optimized using a numerical groundwaterflow and transport code (GFTC) as the basis for forecasting the effectiveness of different pumping patterns. SUTRA (11),a hybrid finite-element/finite-differencemodel developed by and available through the U.S. Geological Survey, was used to solve the governing equations for confined areal groundwater flow and areal solute transport. The irregular numerical grid contained 1650 nodes to describe the site shown in Figure l a and portions of the surrounding basin, with the greatest density of nodes in the area of highest contamination. Although it supports only 2 D simulation, SUTRA permits the aquifer to be defined by variable thickness values. Local modifications to the public domain source code permit testing of drawdown. Heterogeneous values for the key feature of the subsurface medium, hydraulic conductivity, can be represented at each node. Referring to Figure la, the flow boundaries were noflow fault zones to the northeast and southeast, with flux boundaries directly east and more distantly downgradient to the west. Severalinnuent and effluent streams,municipal and irrigation wells, and constant annual recharge from rainfall were included. Variable hydraulic conductivitywas represented by four values corresponding to four main hydrogeologic zones, with a typical value being 3 mlday. A steady-state flow model was originally calibrated to the piezometric surface, which has had fairly constant topology over time. With no remediation, the low-concentration plume on the left side of the figure migrates westerly along the arroyo, while the major plume moves mainly westward but with a marked southwesterly component as well. Detailed descriptions of the site are contained in Thorpe et al. (12) and Tompson et al. (13). Engineering designs to contain and clean up the contamination at the site have gone through many stages as more field data have been collected and more sophisticated tools for modeling and prediction have been employed. The 23 extraction and 5 injection well locations shown in the figure were chosen at one point as a design that could contain the contaminantplume within its current boundaries and clean the site up to regulatory standards within 50 years. The placement of extractors served both contaminant removal and hydraulic control purposes. Injectors served both to aid contaminant removal and to prevent dewatering. Since the cleanup goal was only attainable by assuming that natural degradation processes will remove about half the contamination and there was reason to question the applicability of that assumption to this site, the emphasis for this example was shifted from cleanup defined in regulatory terms to removal of mass. The management question forming the basis for the optimization was as follows: What is the lowest-cost subset of these 28 locations
that will contain the plume behind its current boundaries and extract as much solute mass as the f u l l design over the course of 50 years? While the particular numerical model used in this example required 2-3 h on a SUN SPARCstation I1 to assess one pumping pattern, the trained A " s evaluated roughly 2000 patterns per second. As a result, the GA was able to evaluate over 4 million of the 268 million possible patterns in 36 h on the same computer. Three top-ranked patterns resulting from this search are shown in Figure lb. All three met the regulatory constraint of containing the contamination behind its original boundaries and extracted approximately as much mass as the full 28-well pattern.
However,because they employed so few wells and pumped substantially less water to the surface, the installation, maintenance, and treatment costs associated with these patterns over 50 years were a small fraction of the original plan. If the cost of pumping all 28 locations for 50 years is $155 million, the patterns in Figure l b yield savings of $102-$114.2 million. These patterns certainly might have been found by expert judgment, trial-and-error, or some other informal design strategy. But the systematic evaluation of millions of alternatives that is made possible by the A " s increases the likelihood of finding these minimalist patterns.
Methodology The majority of ANN optimization applications use A " s to distinguish "good"or optimal patterns from nonoptimal ones (14,15). In the current approach, A " s are trained to predict information, such as whether or not a pattern will contain the plume or the volume of solute mass it will extract, that would normally be generated by the GFTC. The nets learn to associate pumping patterns with model predictions from training examples created by the running of the GFTC on a representative sample of pumping patterns. A test set is used to assess the A " ' s ability to generalize its knowledge to new pumping patterns. Once trained to acceptablelevels of accuracy on the test patterns, the A " s stand in for the original GFTC, predicting the effectiveness of new pumping patterns generated by the GA. It is the GA that is responsible for generating increasingly successful patterns and which eventually determines the optimal set of patterns. A flow chart summarizingthe elements of the methodology is shown in Figure 2. Establishing the GETC Ihowiedge Base. SUTRA evaluated 325 pumping patterns. The patterns representedboth a random sample from the domain of possible subsets of the 28 locations and patterns thought by a hydrogeologist to reflect important boundary conditions and good pumping strategies. In this way, the A " s can take advantage of professional expertise when it is available. Results of these simulations formed the knowledge base from which 275 training and 50 testing examples for the A " s were drawn. Another important feature of this approach is that timeconsuming simulations are mined for the maximum amount of information they can yield. Various A " s can be trained and used in separate searches having different objectives without repeating this initial knowledge base creation step. In the current remediation example, the management question was originally formulated in terms of the time it would take to clean up the site to regulatory standards. When the decision was made to shift to mass extraction, it was possible to train new A " s for this prediction task without repeating the SUTRA runs because the needed mass extraction informationhad been saved in the knowledge base. TrainingtheANNs. ANNs were trained to predict model outcomes relevant to each of the three objectives in the management question: (1) whether or not a pumping pattern held the contamination at or below 10 ppb at a line of monitoringwells throughout the 50-year duration of the cleanup; (2) the total solute mass (in kg) extracted by a pattern over 50 years; and (3)the cost of implementingthe pattern. Cost, which was calculated by a function separate from the GFTC itself, included calculations of facilities, piping, operating, and water treatment. Capital costs VOL. 29, NO. 5, 1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
1147
I
I
Generate training patterns
1
I
Ground water flow and transport code (GFTC)
I
Run training f low-and-transport simulations in perallel c I W
Train ANNs to predict objectives Of remediation (Le., outcome of flow and
~~~1
cost objective
.
Artifkiai neural network (ANN)
regulatory objective
inltlaiize population of pumping patterns with fitness measure as a sum of ail three obiectives
Evaluate population fitness for cost objective
Evaluate population fitness for regulatory objective
associated with injection wells were higher than those of extraction wells ($75 000 vs $63 000 per well). The cost of treating groundwater pumped to the surfacewas averaged over several treatment methods to $0.25/L. The type of ANN used in the methodology is a multilayer perceptron that learns by the back-propagation algorithm (16). A simplified diagram of the network used to predict plume containment is shown in Figure 3. The top or input layer uses one node to represent each of the 28 prospective well locations. Each input node, which can be thought of as an independent variable, is valued at 1 or 0, depending on its on or off status in the pumping pattern. This simple binary scheme was employed because for the current problem the pump rate at each location was fixed at the maximum appropriate level (anywherefrom 10 to 50 gpm)for the location. Had pump rate been a design variable being optimized, the input nodes could have taken on continuous values. A 29th node, not shown in the figure, was used to supplythe networkwith additionalinformation 1148. ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 29, NO. 5 , 1 9 9 5
Evaluate population fitness for mass objective
Genetic algorithm (GA)
about the pattern as a whole. The bottom or output layer represents the dependent variable that the network is learning to predict, in this case containment expressed as a dichotomy. For mass and cost, the output node was a continuous variable. The basic task of the ANN is to come up with a set of weights, one for each connecting line in the diagram, that will predict the dependent variable from the input variables as accuratelyas possible. This is exactly the same problem as is solved by conventional linear or nonlinear regression techniques. However, regression techniques fit data to some prespecified model. ANNs empirically develop a model, expressed by the weights, that fits the data. As a result, they can be used to predict complex, poorlyunderstood relationships. These complexitiesare captured in the middle or hidden layer. Unlike the input and output layers whose sue is fured by the problem, the number of nodes in the hidden layer is set during the course of training to whatever size leads to the most accurate predictions.
Input layer (28 fixed locatlons either ON or OFF)
Hidden layer
Output layer (plume contained? yln) FIGURE 3. Simplified diagram of an artificial neural network (ANN).This net is used to predict whether or not a particular pumping pattern will contain the contamination. In the actual representation,every node in the input is connected to every node in the hidden layer; however, for visual simplicity only the connections of the first and last input nodes are shown in the diagram. A 29th input node, also not shown, is used to indicate each pattern's proportion of nonpumping wells, a summary statistic that was shown to enhance the ANNs' predictive . . accuracy.
The mechanics of how a back-propagation ANN learns from examples are beyond the scope of this discussion. Several good introductory texts are available to the interested reader (I7-19].What is most important to remember is that the training of an ANN is a process of adjusting the weights until the outcome predicted by the ANN matches (within some specified tolerance) the outcome predicted by the GFTC on the training examples. This final set of weights is then applied to examples from the testing set to assess the network's predictive accuracy. Obtaining an ANN with an acceptable level of predictive accuracy requires the systematic variation of several network training parameters. Table 1 presents selected accuracy results for all three ANNs. Accuracyon the training set was defined by the product-moment correlation (for the cost and mass variables) or the proportion-of-agreement (forthe regulatoryvariable)between the ANN and the GFTC predictions on the 275 examples used for training. Testing set accuracy is the same correlation/agreementon 50 new examples. In the current problem, accuracy was affectedprimarily by architecture and length of training. The architecture specification in the table lists the number of nodes in the input, hidden, and output layers, with variations in the hidden layer being the component that is manipulated to improve accuracy. A highly linear variable such as cost shows no accuracy differencesbetween the 29-7-1 and 291-1 architectures. The most nonlinear variable, the regulatory constraint, is sensitive to the number of hidden nodes and shows best performance at the 29-7-1 level. The mass variablefalls in between these two extremes, requiring only four hidden nodes to capture the relationship. The function evaluations column in the table reflects how long the networks were allowed to refine their weights on the training examples. In general, a high degree of accuracy on the training set is a prerequisite for good accuracy on the testing set. However, the more time a
TABLE 1
ANN Performancea on 275 Training and 50 Testing Examples architecture
no. of function evaluations
training sat accuracy
testing sat accuracy
ANN No. 1: Regulatory Constraint 29-7-1 29-7-1 29-1-1 29-1-1
1000
5000 1000 5000
0.999 0.990 0.959 0.959
0.982" 0.982 0.809 0.821
ANN No. 2 Total Mass of Solute Extracted 29-7-1 29-7-1 29-4-1 29-4-1 29-1-1 29-1-1 29-7-1 29-1-1
1000 0.998 5000 0.999 1000 0.992 5000 0.998 1000 0.965 5000 0.966 ANN No. 3 Cost of Remediation 1000 0.999 1000 0.998
0.952 0.865 0.965* 0.925 0.952 0.952 0.998 0.999"
a Performance values for ANN no. 1 are proportions of agreement; for ANNs nos. 2 and 3, they are correlation coefficients. An asterisk ("1 indicates the particular network incorporated into the GA.
network spends refining its weights to reflect the idiosyncracies of the training set, the greater the likelihoodthat it is simply memorizingthe input-output pairings rather than capturing critical features needed to generalize its knowledge to the unseen patterns in the test set. This time, it was the mass variable that was most vulnerable to this overfitting. Network training, in summary, is a largely experimental task that continues until acceptable levels of accuracy on the testing set examples are achieved. [For interested readers: The specificvariant of back-propagation used in the current example was the conjugate gradient method (23),an algorithm designed to speed convergence and avoid VOL. 29, NO. 5,1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY m 1149
entrapment in local minima. The learning rate was held constant throughout training at 0.1. A sigmoid function with a slope of 1 and a range of 2 served as the activation function.] In all cases,the particular networkvariation with the best testing set accuracy (0.982 for regulatory,0.965 for mass, and 0.999 for cost) was selected for inclusion in the search phase of the optimization. They are identifiedwith an asterisk (*) in the table. Searchingfor OptimalPatterns. Once trained,the three A " s were linked to a GA. Similar to the A " s , the GA represented a pumping pattern as a string of 28 bits, each bit indicating whether or not the corresponding well location was on or off. The goal of the GA was to find patterns that maximized a fitness function defined as the sum of three equally-weighted components: meeting the regulatoryconstraint (yes = 1,no = 01, the amount ofsolute mass extracted over 50 years (normalized to the highest mass extractionvalue found in the originalknowledge base and scaled to 0.1-0.9), and 1-the cost (normalized to the most costly pattern in the knowledgebase and scaled,again, to 0.1-0.9). A score of 2.8 would be attained if a pattern met the regulatory constraint, extracted as much mass as the highest extractionvalue found in the original knowledge base, and yet magically cost nothing. The GA provides several mechanisms for handling constraints (20). The mechanism adopted here of simply including it as one element in an overall fitness function is analyzed below in the section on Multiple-Objective Search. Because it uses mechanisms derived from concepts of natural selection, the GA operates on populations of patterns rather than a single pattern at a time. The initial population in this example consisted of the 239 pumping patterns (Le., strings) from the original knowledge base that met the regulatory constraint. This population was cycled through 20 000 generations of the GA using simple versions of the three basic GA operators: reproduction, crossover, and mutation (21, 22). The reproduction operator ranks strings according to their total fitness scores, as defined above. The higher a string's fitness, the greater the likelihood that it will be reproduced in the next generation. A roulette wheel selection method was used to copy stringsinto the survivor's pool, maintaining a constant population size of 239. The crossover operator determines how surviving strings will pair for mating. The randomly-selected crossover site between two mated stringsis the point at which both strings are split into two parts. The first part of string no. 1 is joined to the second part of string no. 2 (and vice versa) to form two new strings. The mutation operator serves to keep the characteristics of highly fit strings from completely dominating the mating pool, forcing out less promising but potentially valuable information. In this simple GA, mutation was defined as the random switching of bits in a string and was used at a rate of 1 mutation per 1000 bits. An example of reproduction, mating, and mutation is given in Figure 4. Although the regulatoryconstraint was simply included in the fitness function for this problem, some constraints can also be handled when the crossover and mutation operators create new patterns. If, for example, a limit of 20 wells was placed on the total number of pumps that could be turned on at one time, new strings which violated that constraint could either be discarded or repaired according to some rule. 1150
ENVIRONMENTAL SCIENCE &TECHNOLOGY / VOL. 29, NO. 5,1995
In the GA, search continues until some arbitrarycriterion is met. In this example, the algorithm was allowed to run for 20 000 generations. Anytime a pattern having a total fitness score above 2.45 appeared in the population, it was saved for later analysis. The highest score ever achieved was 2.486 and was found at the 11 246th generation. The patterns in Figure l b were 3 out of over 100 with similar levels of fitness. Since the A " s serve as substitutes for the GFTC, they could be harnessed to any one of several search techniques. The GA was selected for use in this example because of its flexibility and robustness for field applications. First, by searching from a population of pumping patterns and using probabilistic rather than deterministic rules to generate new patterns, the GA reduces the odds that the search will terminate prematurely in what are called "local minima". Second, the researcher has wide latitude in the formulation of the GA's fitness function, being free of assumptions regarding continuity and differentiability that gradientbased search procedures impose (24). These properties are shared by the similarly robust technique of simulated annealing (25). Verification of Optimal Patterns. Since the fitness information driving the search is being estimated by the A " s rather than calculated directly by the GFTC, it is necessary to verlfy each pattern being seriouslyconsidered for implementation. Results of SUTRA runs on the top 10 patterns confirmed that all met the regulatory constraint. The correlation coefficient between the ANN and GFTC predictions of total extracted mass was 0.97.
Analysis of Multiple-Objective Search Since the GA imposes no restrictionson the fitnessfunction it is trying to maximize, it becomes all the more important that the formulation chosen for a particular problem be closely scrutinized. In the current example, the three components of fitness all have some degree of correlation with each other, a situation that could have considerable bearing on the outcomes of the search. Relationships among the three components drawn from the original 325 model runs are as follows: the regulatory constraint measure was moderately correlated with amount of mass extracted ( r = 0.6231, the regulatory constraint vs cost correlationwas 0.352, and the mass extractionvscost correlation was 0.449. In addition, the choice to handle the regulatory constraint as simply another element of fitness rather than an absolute constraint deserves further analysis. For these reasons, the search phase was repeated six times to test the sensitivityof the optimization results to changes in the weights assigned to the components of fitness. Table 2 shows the 25 highest scoring patterns found in the original search, where each component was given equal weight. The three patterns identified by asterisks are the examples shown in Figure lb. The winning 8-well pattern having appeared by the 11 246th generation, the followup searches were ended at 12 000 rather than 20 000 generations. All other search parameters, including the seed used to initialize the random number generator, were held constant. Results from the two searches testing the influence of the regulatory constraint component, by first assigning a 0.0 and then a 2.0 weight to that variable, produced 25 top-scoring patterns identicul to those in Table 2. In other words, the effect of the moderate correlation between meeting the regulatoryconstraint and extracting mass was
A Current C m t i n n No. PunpingPattern
Total Fitness
w1 w75
1111111110011000011111111100 0000011100000000100111111111 X l l O 1 1 00 0 1 00 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 t239 1 1 0 1 1 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1
1.943
Probability of Selection
1.245 1.921 1.564
0.024 0.015 0.024 0.019
110
5
B Matine Pa01 After Reamduction
#1 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 1 0 0 #110 1 1 0 0 0 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0
parent1#1
11111 11110011000011111111100
parent2 #110 1 1 0 0 0 10011111011111111111000
child 1 child 2
1111110011111011111111111000 1100011110011000011111111100
D NewGeneration No. Pumping Pattern
#2
1111110011111011111111111000 1100011110011000011111111100
It1 #2
1110110011111011111111111000 1100011110011000011111111100
#1
FIGURE 4. Example of reproduction, mating, and mutation in a GA population. (a) Four strings from a population of 239 are shown. Total fitness was a linear combination of (1) meeting the regulatory constraint (1 = yes, 0 = no); (2) the normalizedmass of solute that was removed over the 50 years; and (3) one minus the normalized value of the cost of remediation over 50 years. Probability of selection is each string's fitness divided by the total fitness for the entire population. (b) Two of the example patterns with higher probabilities make it into the mating pool. No. 1 randomly draws no. 110 as a mate. The randomly selected crossover point is after position 5. (c) The parent strings are split at the crossover point and recombined to form two children. (d) The children become members of the next generation. (e) A position's status "mutates" from on to off (or vice versa) at an average rate of 1 switch per loo0 positions.
such that any pattern with good mass-extraction performance characteristics was also capable of containing the plume. Neither doublingnor zeroing out the contribution made any difference to the final outcomes of the searches. The mass and cost components,on the other hand, were affected by the weighting schemes. Table 3 shows the results of zeroing out mass performance. This corresponds to the management question: What is the cheapest way to contain the plume? As the top five patterns show, containment does require fewer wells than mass extraction but not dramatically so. The correlation between containment and mass extraction acts to keep the mass-extraction scores of the top patterns from being totally random, but they are markedly lower than those in Table 2. A final comment concerns the role of injection wells. In Table 2, only two of the top 25 patterns contain an injection well. In Table 3,12 of the top patterns contain an injector. [Only one of the five injection wells, no. 26, ever appeared in any ofthe top-rankedpatterns.] This suggests that re-injection,
at least in the scheme shown in Figure la, reduces the efficiency of mass extraction but does not interfere with containment. Doubling the weight assigned to mass extraction (see Table 4) has the expected effect of giving preference to the better mass-extraction patterns over the cheaper patterns. Note too that with the balance thrown to mass extraction over cost, the 8-well and 9-wellpatterns from Figure l b are lost. The most interestingresultscome about by manipulating the cost component. Table 5 might be called the moneyis-no-object search, since zero weight is assigned to cost. Predictably,the mass-extraction component dominatesthe search, and the winning patterns all have high numbers of extractor locations. The values of 0.909,0.908,etc. for mass indicate that these patterns extract slightly more mass than the best mass-extractionvalue in the original knowledge base. The appearance of one particular injection well, no. 26, in the top patterns in this table adds some qualification VOL. 29, NO. 5,1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
1151
TABLE 2
25 Highest Scoring Pstreras under the Equal Weights Weighting Scheme: Fitness = l*Regulatory 1"Mass 1"Cost
+
+
normalizedC total fitness
pumping pattern"
ext wells
inj wells
rag flagb
mass
cost
2.486* 2.479 2.476 2.475 2.474" 2.472 2.472 2.471 2.471 2.470 2.470 2.469 2.468 2.468 2.467 2.467 2.467 2.466 2.466 2.465* 2.465 2.464 2.464 2.464 2.464
00111000011001000110000 00000 1111101101100100011000000000 1111100001100110011000000100 1111101101100010111000000000 1111101101100110011000000000 0111101001100111101100000000 1111101101100000010000000000 1111101101100010001100000000 11111011011001100010000 00000 01111001011000110011000 00000 01 111010010000100110000 00100 11111010111001010110000 00000 01111000010000110110000 00000 1111101000100010011100000000 0111101101100100011100000000 0111101111100110001000000000 1111101101110000011000000000 11111001111001100110000 00000 1111101111100010001000000000 01111000010001000111000 00000 1111100111110000011000000000 0111101101100010111100000000 0111101101110010011000000000 1101100111110000011000000000 11111000011100000011000 00000
8 12 11 13 13 13 10 12 12 11 9 13 9 11 12 12 12 13 12 9 12 13 12 11 10
0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.797 0.834 0.829 0.838 0.841 0.834 0.801 0.829 0.794 0.816 0.796 0.840 0.780 0.838 0.865 0.803 0.819 0.847 0.798 0.815 0.831 0.870 0.819 0.812 0.791
0.689 0.645 0.648 0.637 0.633 0.637 0.671 0.643 0.678 0.654 0.674 0.629 0.687 0.630 0.602 0.664 0.648 0.619 0.668 0.649 0.635 0.594 0.644 0.652 0.673
0
0 0 0 0
Positions 1-23 of the pattern string reflect the on/off status of the extraction wells. The five remaining positions indicate the status of the injection wells. The number (1)indicates that the regulatory constraint was met. Both mass and 1-cost have been scaled t o the range 0.1 -0.9. In both cases, goodness is indicated by a higher number. So, the cost figure may be thought of as a cheapness index. An asterisk ( * ) indicates a pattern shown in Figure I b . a
TABLE 3
25 Highest Scoria4 Patterns under tbe Cost of Containment 1"Regulatory 0"Mass 1"Cost
+
+
Ow Weighting Scheme:
Fitness = normalized
total fitness
pumping pattern
ext wells
inj walls
reg flag
(mass)'
cost
1.739 1.723 1.720 1.720 1.715 1.712 1.709 1.709 1.709 1.707 1.705 1.704 1.704 1.703 1.703 1.702 1.702 1.702 1.702 1.701 1.701 1.700 1.700 1.700 1.700
0110000000110000011000000100 0101101000110011001000000000 1111000100000011001000000100 11111011000000110000000 00000 00110010100100100010000 01000 11111010001100010010000 00000 0111001001110011001000000000 0111101010000010001000000100 1011101010000010001000000100 0111101100010011001000000000 1101101100100110001000000000 0111101100010110001000000000 11110010001100110010000 10000 0111101000010111101000000000 0111101100110010001000000000 11110000001101110010000 00100 11110010011100110010000 00000 1111001010110100001000000000 11111010001100110010000 00000 01101010111001110010000 10000 1011100000010111101000000100 01110010011000100010000 01100 1110101101010001001010000000 1111101010000010001000000100 1111101100000010001000001000
6 9 8 9 7 10 10 8 8 10 10 10 10 11 10 10 11 10 11 11 10 8 11 9 9
1 0 1 0 1 0 0 1 1 0 0 0 1 0 0 1 0 0 0 1 1 2 0 1 1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
(0.488) (0.665) (0.673) (0.695) (0.302) (0.703) (0.720) (0.739) (0.672) (0.706) (0.697) (0.706) (0.620) (0.721) (0.726) (0.700) (0.718) (0.662) (0.725) (0.652) (0.678) (0.653) (0.584) (0.732) (0.586)
0.739 0.723 0.720 0.720 0.715 0.712 0.709 0.709 0.709 0.707 0.705 0.704 0.704 0.703 0.703 0.702 0.702 0.702 0.702 0.701 0.701 0.700 0.700 0.700 0.700
a
Mass scores for each pattern are given in parentheses, even though they were not included in the total fitness measure.
to the earlier speculation concerningthe effect of injection of mass extraction. Apparently, when cost is not a factor, the injector at location no. 26 actually facilitates mass 1152
ENVIRONMENTAL SCIENCE &TECHNOLOGY / VOL. 29, NO. 5,1995
extraction. When cost efficiency is folded into the management formulation, as is the case in Table 2, no injectors facilitate mass extraction.
TABLE 4
25 Highest Scoring Patterns under the Emphasis on Mass Extraction Weiyhtiny Scheme: Fitness = l*Regulatory 2*Mass 1"Cost
+
+
normalized total fitnessa
pumping pattern
ext walls
inj wells
3.336 3.335 3.333 3.330 3.327 3.327 3.327 3.326 3.325 3.324 3.323 3.321 3.321 3.320 3.319 3.317 3.316 3.31 5 3.315 3.315* 3.315 3.315 3.314 3.314 3.313
1111101011100010011100000000 01111011011000101111000 00000 0111101101100100011100000000 11111011011001000111000 00000 1011101011100110011100000000 11111001011000101111000 00000 11111011011001101111000 00000 11111011111001000111000 00000 0111101111100110011100000000 111 110111110001001 11000 00000 01111011011000101111000 00100 0111101111100110111100000000 11111011111001100111000 00000 1111100011110100111100000000 1111100011100010111100000100 1111101101100010111100000100 0111101111100000011100000100 1011101011100110011100000100 1111100001100011011100000100 1111101101100110011000000000 1111101101100111011100000000 1111101111100000011100000100 0111101101100011011100000100 01111011111000110111000 00000 11111001111001100110000 00000
13 13 12 13 13 13 15 14 14 14 13 15 15 14 13 14 12 13 12 13 15 13 13 14 13
0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 1 1 1 0 0 1 1 0 0
reg flag
mass
cost
1
1.749 1.741 1.731 1.737 1.745 1.732 1.762 1.762 1.765 1.759 1.745 1.774 1.773 1.746 1.747 1.751 1.738 1.749 1.715 1.682 1.751 1.749 1.737 1.750 1.694
0.587 0.594 0.602 0.593 0.582 0.595 0.564 0.564 0.560 0.565 0.578 0.547 0.547 0.574 0.572 0.566 0.578 0.566 0.600 0.633 0.564 0.566 0.577 0.564 0.619
1
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
An asterisk ( * ) indicates a pattern shown in Figure lb.
TABLE 5
25 Highest Scoring Patterns under the Money Is No Object Weighting Scheme: Fitness = 1"Regulatory l*Mass 0"Cost
+
+
normalized total fitness
pumping pattern
ext wells
inj walls
reg flag
mass
(cost)'
1.go9 1.go9 1.908 1.908 1.908 1.907 1.907 1.907 1.907 1.907 1.907 1.907 1.907 1.906 1.906 1.906 1.905 1.905 1.905 1.905 1.905 1.905 1.905 1.905 1.905
11111011111011101111101 00000 11111011111111101111101 00100 11111011111011101111111 00100 11111011111111101111101 00000 11111011111111101111111 00100 1111101111101100111110100100 1111101111101110111111100000 11111011111011111111101 00000 11111011111011111111101 00100 1111101111101111111111100100 1111101111111110111111100000 11111011111111111111101 00000 11111011111111111111101 00100 1111101111101111111111100000 11111011111111001111101 00100 11111011111111111111111 00000 10111011111111101111101 00100 1111100111101110111110100100 11111011111001101111101 00100 11111011111011001111111 00000 11111011111011011111101 00100 11111011111011100111101 00100 11111011111101101111101 00000 11111011111111001111111 00000 11111011111111011111101 00000
19 20 20 20 21 18 20 20 20 21 21 21 21 21 19 22 19 18 18 19 19 18 19 20 20
0 1 1 0 1 1 0 0 1 1 0 0 1 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.909 0.909 0.908 0.908 0.908 0.907 0.907 0.907 0.907 0.907 0.907 0.907 0.907 0.906 0.906 0.906 0.905 0.905 0.905 0.905 0.905 0.905 0.905 0.905 0.905
(0.366) (0.338) (0.297) (0.355) (0.284) (0.366) (0.312) (0.355) (0.339) (0.284) (0.299) (0.341) (0.325) (0.299) (0.355) (0.287) (0.351) (0.364) (0.407) (0.325) (0.355) (0.362) (0.409) (0.315) (0.358)
1
0 1 1 1 0 1 1
0 0 0
Cost scores for each pattern are given in parentheses, even though they were not included in the total fitness measure.
Finally, Table 6 examines the case where the cost component is given twice its original weight. Predictably, preference is given to the less expensive patterns: and the 8-well pattern from Figure l b re-emerges as the winner. The 9- and 13-well patterns from Figure lb, having been
chosen from the original list in Table 2 primarily for their mass-extraction performance, do not appear in Table 6. Each followup search required 22 h to complete. In that span of time, the original GFTC could have been called only 8.8 times. VOL. 29, NO. 5. 1995 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
1153
TABLE 6
25 Higkest Scming PlrttsRII d e r the Money Doesn't Grow on Trees Weighting Scheme: Fitness = 1"Regulatory 1'Mass 2"Cost
+
+
normalized total fitness'
pumping pattern
ext wells
inj wells
reg flag
mass
cost
3.176* 3.157 3.155 3.149 3.145 3.144 3.144 3.144 3.143 3.142 3.142 3.140 3.139 3.138 3.137 3.137 3.135 3.134 3.134 3.133 3.133 3.133 3.132 3.131 3.131
00111000011001000110000 00000 01111010100000100010000 00100 01111000010000110110000 00000 11111011011001100010000 00000 11111011011000110010000 00000 01111000011000100010100 00000 01111010010000100110000 00100 11111011001000110010000 00000 11111011011000000100000 00000 0111101010100001001000000100 11111001011000111010000 00000 11111000010001110010000 00100 11111010011100110010000 00000 01110010011100110010000 00000 10111001111100001010000 00000 11111000011100000011000 00000 1111101011110100001000000000 11111011000000110000000 00000 1111101111100010001000000000 01111011001000100100000 00000 0111101101100101001000000100 11111010100000100010000 00100 11111000010000110110000 00000 01111011001100100010000 00000 0111101111100110001000000000
8 8 9 12 12 9 9 11 10 9 12 10 12 10 11 10 12 9 12 9 11 9 10 10 12
0 1 0 0 0 0 1 0 0 1 0 1 0 0 0
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0.797 0.739 0.780 0.794 0.783 0.774 0.796 0.753 0.801 0.747 0.778 0.745 0.764 0.720 0.775 0.791 0.787 0.695 0.798 0.765 0.778 0.732 0.777 0.726 0.803
1.379 1.418 1.375 1.355 1.362 1.370 1.348 1.392 1.343 1.395 1.364 1.395 1.375 1.418 1.362 1.346 1.348 1.440 1.336 1.367 1.354 1.401 1.355 1.405 1.327
a
0 0 0 0 0 1 1 0 0 0
A n asterisk ("1 indicates a pattern shown in Figure l b .
Discussion Comparisons with Other Search Techniques. It was mentioned earlier that, while SUTRA required an average of 2.5 h to evaluate a single pumping pattern, the trained A " s could evaluate patterns at a rate of 2000 patternsls. These rates explain how the GA could evaluate 4 million patterns in 36 h on a conventional Unixworkstation. They do not, however, reflect the entire process from the creation of the knowledge base through the search. To do that, it is more instructive to compare the total number of GFTC runs required by different approaches. Earlier work (26, 27) on a 20-well, single-objectiveremediation example compared the G A - A " approach described here with a nonlinear optimization procedure, NPSOL, that has been used in several groundwater optimization studies. 245 SUTRA runs were required to create the initial knowledge base; NPSOL converged after 206 SUTRA calls. NPSOL's single optimal pattern was almost identical to one of the top 10 patterns generated by the GA. So, it would appear that NPSOL solves the problem more efficiently were it not for the following considerations: (1) In the G A - A " method, the model runs are completely independent of each other and are run in parallel on different machines; some degree of sequential dependency in the calls to the model is inherent in all other approaches except random search. (2) Gradient-based nonlinear optimization methods do not guarantee that the resulting solution is the "global optimum". As a result, the search must be repeated from different starting points to increase the likelihood of finding the best overall solution. Since each repetition of the search is likely to require a similar number of calls to the model, the one-search advantage in number of model runs shown by NPSOL is quickly lost. (3) In field situations, engineers and planners change their minds. So, searches have to be repeated not only to verify 1154
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 29, NO. 5,1995
the solutions but also to reflect changing management formulations. Within certain bounds to be discussed below, different management formulations can be accommodated by the G A - A " approach without recreating the original knowledge base. Caveats on ANN Usage. The most critical step in the G A - A " methodology is the creationof the knowledge base. The A " s can only generate meaningful predictions over the problem space defined by the initial modeling runs. If a major change in the definition of the problem is required, such as increasing the pool of prospective locations or the timeframe over which predictions are being made, the knowledge base must be augmented or recreated. Arelated caution concerns the degree of discrimination the A " s can display. If the original sampling of model runs is too sparse over a region in the problem space that later turns out, in the course of the search phase, to hold promising patterns, the A " s will be unable to distinguish subtle differences among the patterns in that region. The patterns will all be seen as comparable, and no further refinements will be made. The best way to avoid this problem is to develop the knowledge base iteratively, beginning with a crude sampling of the problem domain and adding more coverage to particular regions as they show promise of containing effective patterns. As described here, optimization is a value-adding procedure, exploitingwhatever level of complexityhas been put into the construction of a contaminant transport model for a given site. As the models improve, so does the quality of the results obtained through the optimization. However, since the advantages of the G A - A " approach rest on the abilityo f A " s to function as rapid yet accurate substitutes for a model, the impact of increasing model scale and complexity is an important issue to address. Current applications of the G A - A " approach are only available for 2D, single-phase groundwater simulations. However,
ANNs have been successfullyused to predict complexphysical phenomena such as seismic patterns associated with earthquakes vs underground explosions (28)and sun spot activity (29).Furthermore,theoretical proofs exist to show that a multi-layer perceptron with one hidden layer can approximateany continuous function and that at most two hidden layers are required to approximate any particular set of functions (30). So, complexity by itself is not necessarily an obstacle to the successful application of A " s . Theory aside, there are two additional factors that may explain how a complex, nonlinear model can be estimated by A " s , as in the current problem. First, optimizationis a site-specific activity. Most of the initial conditions are held constant, while only the impact of the design factors (number,location, and function ofwells, for example)need to be estimated by the A " s . Second, the data being predicted are generated by a deterministic model. The usual noise stemming from unreliability of measurements and the high variability often associated with natural phenomena are not encountered here. Current Directions. Work is currently underway to apply the methodology to a 225-location version of the 2D design problem described here and to 3D multiphase models. The major challenge that has been encountered so far involves completing the initial modeling runs in a reasonble amount of time. As the number of design factors increases, so does the size of the initial knowledge base. Also, the increasing scale and complexity of the models increase the amount of time required to complete a single run. These problems are currently being addressed by exploitation of the networks of workstations available at Lawrence Livermore National Laboratory and elsewhere. A simple job-scheduling routine, using no additional resources beyond those routinely available through standard network services, distributes model runs for simultaneous execution among 50-75 workstations. Since the runs are independent of each other, no special messagepassing software such as PVM (31)need to be installed. This system is so simple to implement that the only real barrier to a muchwider utilization of computersthroughout the organization involves public relations. Other approaches that will be investigated in the near future are sampling schemes drawn from the experimental design literature to reduce the number of runs needed to cover the problem space and up-scaling techniques for reducing the resolution (and, as a result, time) needed by the models.
Conclusions G A - A " management methodology is a robust and flexible tool to aid decision making in groundwater remediation applications. It begins and ends with the classical basis of groundwater transport prediction, the contaminant transport code. However, harnessingthe abilities of the genetic algorithm and artificial neural networks allows several million pumping patterns to be examined with the same computational burden that traditional steepest descent optimization algorithms use to evaluate several hundred pumpingpatterns. This renders heretofore intractable fieldscale applications practical.
Acknowledgments We are grateful for permissionto use the back-propagation code of Erik Johansson and Sean Lehman of Lawrence Livermore National Laboratory (LLNL). The paper benefited from critical review by Richard Knapp and Kenneth Jackson of LLNL. This work was made possible by support
from the LLNL Environmental Restoration Division under the auspices of the U.S. Department of Energy under Contract W-7405-Eng-48.
Literature Cited (1) Schneider, K. N.Y. Times 1993, 142 (Mar 21), 1. (2) Richards, B. Wall Sf.I. 1992, 220 (Oct261, 6. (3) Abelsen, P. H. Science 1993, 259, 159. (4) Gorelick, S. M. Water Resour. Res. 1983, 19, 305. (5) Gorelick, S. M.; Freeze, R. A.; Donohue, D.; Keely, J. F. Groundwatercontamination: optimal capture and containment; Lewis Publishers: Boca Raton, FL, 1993; pp 1-13, 101-211. (6) Wagner, B. J. Recent advances in simulation-optimization groundwater management modeling; US. Geological Survey: Menlo Park, CA, 1994. See also Rev. Geophys. (in press). (7) Gailey, R. M.; Gorelick, S . M. Groundwater 1993, 31, 107. (8) Karatzas, G. P.; Pinder, G. F. Water Resour. Res. 1993, 29,3371. (9) Marryott, R. A.; Dougherty, D. E.; Stollar, R. L. WaterResour. Res. 1993, 29, 847. (10) McKinney, D. C.; Gates, G. B.; Lin, M. D. In Numerical Methods in WaterResources; Peters, A. Ed.; Kluwer Academic Publishers: New York, 1994; p 845. (11) Voss, C. I. A finite-element simulation model for saturatedunsaturated, fluid-density-dependent groundwater flow with energy transport or chemically-reactive single-species solute transport;U.S. GeologicalSurvey,Water Resources Investigations Report 84-4369; US. Geological Survey: Denver, CO, 1984. (12) Thorpe, R. K.,Ishenvood, W. F., Dresen, M. D.,Webster-Scholten, C. P., Eds. CERCLA Remedial Investigations Report for the LLNL Livermore Site; Lawrence Livermore National Laboratory: Livermore, CA, 1990; UCAR-10299. (13) Tompson, A. F. B.; Nichols, E. M.; McKereghan, P. R.; Small, M. C. Summary of preliminary ground water simulations in the Livermore Regional Modeling Study: CFESTfinite element code; Lawrence Livermore National Laboratory: Livermore, CA, 1991; UCRL-AR-107049. (14) Hopfield, J. 7.; Tank, D. W. Biol. Cybernetics 1985, 52, 141. (15) Hegde, S. U.;Sweet, J. L.; Levy, W. B. In Proceedings of the International Conference on Neural Networks; IEEE: Piscataway, NJ, 1988; pp 291-298. (16) Rumelhart, D. E.; McClelland, J. L. Parallel Distributed Processing: Explorations in the Microstructure of Cognition; MIT Press: Cambridge, MA,1986; Vol. 1, pp 318-362. (17) Dayhoff, J. E. Neural Network Architectures: An Introduction; Van Nostrand Reinhold New York, 1990. (18) Maren, A. J.; Harston, C. T.; Pap, R. M. Handbook of Neural Computing Applications; Academic Press: New York, 1990. (19) Wasserman, P. D. Neural Computing: Theoryand Practice;Van Nostrand Reinhold: New York, 1989. (20) Michalewicz, 2.;Janikow, C. E. In Proceedings of the Fourth International Conference on Genetic Algorithms; Belew, R. K., Booker, L. B., Eds.; Morgan Kaufman: San Mateo, CA, 1991;pp 151-157. (21) Goldberg, D. E. Geneticalgorithms in search, optimization, and machine learning; Addison-Wesley: Reading, MA,1989. (22) Michalewicz,2. Geneticalgorithms + datastructure= evoluation programs, 2nd ed.; Springer-Verlag: Berlin, 1994. (23) Johansson, E. M.; Dowla, F. U.; Goodman, D. M. Int. J. Neural Syst. 1992, 2, 291. (24) Gill, P. E.; Murray, W.; Wright, M. H. Practical Optimization; Academic Press: London, 1981. (25) Dougherty, D. E.; Marryott, R. A. Water Resour. Res. 1991, 27, 2493. (26) Rogers, L. L.; Dowla, F. U. Water Resour. Res. 1994, 30, 457. (27) Rogers, L. L. Ph.D. Dissertation, Stanford University, 1992. (28) Dowla, F. U.; Taylor, S. R.; Anderson, R. W. Bull. Seismol. SOC. Am. 1990, 80 (5), 1346. (29) Weigend, A. S.;Rumelhart, D. E.; Huberman, B. A. In Proceedings of the 1990 ConnectionistModels Summer School; Touretzky, D. S . , Ed.; Morgan Kaufman: San Mateo, CA, 1990; p 105. (30) Hertz, J.; Krogh, A.; Palmer, R. G. Introduction to the theory of neural computation;Addison-Wesley: Redwood City, CA, 1991; pp 141-147. (31) Geist, A. PVM: parallel virtual machine: a userk guide and tutorial for networked parallel computing; MIT Press: Cambridge, MA,1994.
Received f o r review April 21, 1994. Revised manuscript received December 29, 1994. Accepted January 17, 1995.@ E339402462 @
Abstract published in Advance ACS Abstracts, February 15,1995.
VOL. 29, NO. 5 . 1 9 9 5 / E N V I R O N M E N T A L SCIENCE & TECHNOLOGY
1115
