Application of Genetic Algorithms for the Estimation of Hydraulic Conductivity

The objective methods of parameter calibration are usually based on inverse analysis. The general assumption here is that the estimation should lead to the maximum possible imitation of the experiment by the calibrated model. With the use of inverse methods, the obtained parameters values would, in general, have rather ‘effective’ than ‘physical’ meaning. This is because of the maximal correspondence between the experiment and the simulation result would cause, as the not modelled phenomena would also manifest in the obtained parameters values. This could be a drawback for theoretical analysis, but for practical (engineering) modelling this is actually an advantage, as it allows us to include in the results aspects that are not formally included in the model.

Inverse calibration techniques are based on some search algorithms. Typical implementation would search through the calibrated parameters space to find a set of parameters leading to maximal concurrency with experimental data. The correspondence between both must be expressed as a search criteria and the most common is the least square error between experimental data and simulation results. The typical estimation procedure is presented in Figure 1.

Figure 1

Our experiments are based on a model of multiphase medium, assuming that the solid has hydraulically connected pores or microcracks, which allow for the flow of liquid and gas filtration. A mathematical model of the creep of a porous medium defined as a multicomponent diphase body was first introduced by Biot [13,14]. In a very general way the model has been analyzed on the basis of asymptotic homogenization theory for the case when the diphase medium’s pores are filled with a poorly compressible liquid. For this purpose [3,5,12,21] used the methods for periodic structures while [20,23] used statistical methods. For the case when the medium’s pores are filled with a gas, a mathematical model based on the asymptotic homogenization method was presented by [4]. One of the fundamental assumptions of porous medium mechanics is the porous medium continuity postulate. It is a considerable simplification when this theory is applied to soils. Despite this assumption, models describing the creep of the Biot body represent well the actual processes of deformation of cohesive soils (argillaceous and clayey soils and mudclays). The processes are described also taking into account temperature changes by [26]. In [9] PhD thesis and in [8,10] one can find consolidation equations based on the Biot model with the Kelvin-Voight rheological skeleton, showing that besides volumetric compressibility and shape elasticity, the viscosity of the soil skeleton also plays a significant role.

The classical constitutive relations for a poroelastic body in isothermal problems were given by [14]. The effective parameters of the Biot model (Biot’s constants) were described for the first time in 1957 [15] as A, N, Q and R . The original form of the constitutive relations for a medium consisting of a porous elastic skeleton and a liquid was written as follows:

σ<!-- s -->ij=2Nε<!-- e -->ij+Aε<!-- e -->δ<!-- d -->ij+Qθ<!-- ? -->δ<!-- d -->ijσ<!-- s -->−<!-- - -->σ<!-- s -->a=Qε<!-- e -->+Rθ<!-- ? -->$$\begin{array}{} \displaystyle \sigma_{ij}=2N\varepsilon_{ij}+A\varepsilon\delta_{ij}+Q\theta\delta_{ij}\\ \sigma-\sigma_{a}=Q\varepsilon+R\theta \end{array}$$(1)

where A, N, Q and R are Biot’s constants interpreted in the study of [15].

N is equivalent to shear modulus μ, and A are is equivalent to Lame coefficient λ; R is the modulus of volume elasticity of the liquid; Q is a coupling coefficient stemming from the mutual interaction between the solid phase and the liquid phase, σ_ij is a tensor of stress in the skeleton, related to the total RVE surface, defined as a fuzzy stress tensor, σ is the stress in the liquid, related to (similarly as the stress in the skeleton) the total surface area of the RVE cross section, also defined as fuzzy stress, ε_ij is a skeleton deformation tensor, ε = ε_ii, a skeleton dilatation and θ is a liquid dilatation. Moreover, denoting the pressure in the liquid as p and porosity as f₀, the fuzzy pore stress in the liquid is expressed by the formula σ = −f₀p

Consolidation equations were derived by [14] and have the following form:

N∇<!-- ? -->2ui+(M+N)ε<!-- e -->,i=−<!-- - -->HRσ<!-- s -->,ikf2∇<!-- ? -->2σ<!-- s -->=1Rσ<!-- s -->˙<!-- ? -->−<!-- - -->HRε<!-- e -->˙<!-- ? -->$$\begin{array}{} \displaystyle \left\{\begin{array}{} N\nabla^{2}u_{i}+(M+N)\varepsilon,_{i}=-\frac{H}{R}\sigma,_{i}\\ \frac{k}{f^{2}}\nabla^{2}\sigma=\frac{1}{R}{\dot\sigma}-\frac{H}{R}{\dot\varepsilon} \end{array}\right. \end{array}$$(2)

where: kf2=k,$\begin{array}{} \frac{k}{f^{2}}=k, \end{array}$ with k being the Darcy’s filtration coefficient and f the skeleton’s porosity.

In his postdoctoral thesis [27] put forward an idea of using the genetic computation technique to solve the selected geotechnical engineering problems (the stability of scarps and slopes, the interpretation of the results of triaxial compression tests, etc.). He demonstrated that the state-of-the-art numerical techniques can restore the full functionality of some solutions that have existed in geotechnics for decades, but which have not been applied in practice because of, for example, too high computing costs. Biot’s model calibration based on inverse analysis, carried out by us using the genetic algorithm, has proved the belief that genetic algorithms constitute a highly effective calibration tool for the automatic (based on simple rules and simultaneous for all the parameters) determination of effective parameters, also for the more complex problems of mechanics. The results presented here can be used as the basis for further research on the application of genetic algorithms to the analysis of the effective parameters of the Biot model with the viscosity of the soil skeleton taken into account.

Genetic algorithms (GAs) can be defined as a numerical technique inspired by natural processes of information processing [17]. They are classified as resistant optimization methods, or methods showing similar efficacy for a wide range of issues [19]. This flexibility derives from the fact that the optimization is performed based only on the value of the objective function. The algorithm does not have any more information of the model or of a modeled process [27]. This is a very important property of GA, as it means that for different models geometries, processes, parameters sets and experiments, they would have similar effectiveness and can be used without modification.

The methodology is, in fact, the inverse estimation. During the estimation process, the algorithm generates sets of parameters, which are used for process simulation. Based on results obtained, the estimator least-squared error is applied for valuation of the obtained solutions (Figure 1). Based on error values, more sets are generated and the process is continued until satisfying results are obtained.

GAs have proved to be very successful in many scientific and engineering areas [2,25]. A great theoretical introduction for the calibration with GAs was given by [18]. The complicated model calibration was used by [1,22,28] and many others.

The terminology connected with genetic algorithms includes, beside typical optimization terms, also terms derived from genetics. This is linked with the way genetic algorithms work, which, as mentioned above, imitates natural evolutionary processes.

GAs operate on variables represented by code combinations called chromosomes. Such combinations consist of a set of characteristics (or characters) called genes. The characters belong to an (usually binary) alphabet (selected for the particular application) whose possible states are called alleles. The entire genetic information about a given element is called a genotype.

The data structures (code combinations) on which GA operates are converted into values of the each variable. Then the variables are inserted into the calibrated model and used for simulation and then evaluation by means of the objective function calculation (least square error estimator). The set of parameter values and the evaluation constitute a phenotype. A set of phenotypes processed at a given stage of optimization is called a population. Genetic operators are used when creating next populations. Due to the application of the genetic operators, a newly created generation has a different gene pool (and thus parameters values) than the preceding generation.

Figure 2 shows a block diagram illustrating how genetic algorithms work. Seven principal steps can be distinguished in the diagram [24].

Figure 2

The first step in the operation of GAs is the selection of an initial population. In our case, the initialization was done through the random selection of the values of consecutive genes.

The adaptation of the chromosomes in the population is evaluated on the basis of the adaptation function values of all the individuals in the current generation. This process takes place in several stages. First, using an appropriate decoding method, the chromosomes are converted into values of variables. Then the variables are used to calculate the objective function for all the individuals in a given population.

The stop condition is defined on the basis of the optimization requirements. In the simplest case, the performance of a set number of iterations, the testing of a proper number of points of space, or time out can be the condition. In our case, the algorithm was stopped when there was no significant improvement in the results.

The selection of chromosomes is conducted in order to isolate chromosomes for a parental pool, from which a new population (a set of points) is created using genetic operators. We used the roulette method with the scaling function. Regardless of the method used, selection is always made solely on the basis of the adaptation function values for the individuals.

The chromosomes of the individuals selected for a parental pool are subjected to the operation of genetic operators, as a result of which a population of descendants is obtained (Figure 3). We used two classic genetic operators, i.e. the crossover operator and the mutation operator. The crossover operator operates simultaneously on two chromosomes selected from a parental pool, causing the exchange of fragments of their chromosomes. Mutation is performed on one gene, with an assigned probability, causing a change in its state.

Figure 3

As a result of the operation of the genetic operators, a new population is formed from the individuals. This population becomes a current population and objective function and adaptation function values are calculated for it. Then the stop condition is checked and if the computations need to be continued, another cycle is started. When the stop condition is reached, the result is given as the output. The result is a set of values of the decision variables, calculated for the chromosome with the highest (when a maximum is sought) value of the objective function.

Turoszow clays, obtained from the mine dump soil KWB Turow, due to its structural complexity, were used in the presented research. Because of its structure, this material is extremely interesting for research, and problematic from the perspective of engineering practice. The presence of water associated with the boundary surface of clay particles causes that the mechanical processes occurring in these soils are much more complex and should take into account the influence of clay particles’ contact through bound water. Therefore, defining individual parameters with the greatest accuracy is a key task in the context of the modeling of engineering processes of this soil.

Basic soil parameters are:

–

liquid limit: w_L =37.4%,

The consolidation chamber used in the studies allowed the application of vertical load P and measurement of displacement, filtration flow through the sample and the pore pressure in the middle of height of the sample. The scheme of the measurement station and the picture of consolidation chamber are shown in Figure 4. The load applied to the sample placed in consolidation chamber (1) was produced by the regulator (2) and transmitted through a flexible tube (8). Similarly, pressures applied to the bottom and top plates are generated with the controllers (3) and (4) and transferred to the chamber through the tubes (9) and (10) respectively. The measurement of pore pressure inside the sample was measured with a pressure sensor (11) connected via a tube (12) with the interior chamber. The measurement of the progress of the consolidation was made by the displacement sensor (13) and also, for inspection, by the clock (14). Signals from the measuring devices were transferred by conductors (15), (16) and (17) to the transducer (5) where from the conditioned signal was routed with the cable (18) to PC (6). The computer allowed simultaneous automatic recording of measurement data and visualization of the experiment progress.

Figure 4

The presented consolidation chamber was used for the oedometer test. The pressures applied were: 12.5, 25, 50 and 100 kPa. Initial sample height h₀ = 20 mm and diameter d =75 mm was used. The results of the experiment are presented in Figure 5.

Figure 5

CV=KMγ<!-- ? -->w$$\begin{array}{} \displaystyle C_{V}=\frac{KM}{\gamma_{w}} \end{array}$$(3)

CV=0.196h2t50$$\begin{array}{} \displaystyle C_{V}=0.196\frac{h^2}{t_{50}} \end{array}$$(4)

CV=0.848h2t90$$\begin{array}{} \displaystyle C_{V}=0.848\frac{h^2}{t_{90}} \end{array}$$(5)

Based on the results, the bulk modulus was calculated: M = 1662, 9KPa. Based on consolidation coefficient equation (Eq.3), using Cassagrande (Eq.4, Fig.5) and Taylor (Eq. 6, Fig. 7) methods, hydraulic conductivities were calculated: t₅₀ =193,5min; = 524,6min. The calculated k values are 1.189 × l0^-8 m/s and 1.1898×l0^-8 m/s for Cassagrande and Taylor methods respectively.

Figure 6

Figure 7

Oedometer test results were used as calibration data in the investigations. The calibration itself was carried out using a program implementing the classical genetic algorithm, designed by the authors, and the FlexPDE software for running FEM simulation needed for an inverse analysis. Thanks to this software combination, calibration can be performed exactly on the model, which is later used to solve actual engineering problems. In addition, owing to the FlexPDE software capabilities a complex experiment design involving the gradual loading of an oedometer sample can be used. A mathematical model was implemented in the FlexPDE program.

The algorithm was implemented in the form of an independent program written in the C# language in the .NET environment. The program was written to make it possible to calibrate any mathematical model with the help of an external program supplying simulation results.

The scheme of C# software operation in connection with the FlexPDE program is presented in Fig. 8. At the beginning of its operation, the GA reads laboratory measurement results from a file. During each cycle, GA generates parameter values and inserts it to basic pde file. Then it starts FlexPDE and instructs it to run simulation from this file and save its results. Next, GA reads simulation results and compares it with the data read earlier. The comparison is based on least square error estimator (objective function). The operation is stopped when completion criteria is met.

Figure 8

The calibration was performed with the use of a specially developed software implementing genetic algorithm and FlexPDE. The experimental data from four-stage oedometer test (as described in chapter 6) were used for the calibration. The mathematical model was a third Biot model with a rheological Kelvin–Voight skeleton. The subject of calibration was only one parameter . To check the reproducibility of the method, 20 independent estimations were performed.

The results of the calibration are presented in Table 1. The values of hydraulic conductivities obtained varies, but the relatively small standard deviation σ_k = 4.38 compared to average value k_av = 64.75 proves good reproducibility. The values obtained with inverse optimization technique are one order smaller than values obtained with Terzaghi model.

Table 1

Parameter estimation results.