Parameter identification in numerical modeling of tunneling using the Differential Evolution Genetic Algorithm (DEGA)

Tunnelling and Underground Space Technology 28 (2012) 109–123

Contents lists available at SciVerse ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Parameter identification in numerical modeling of tunneling using the Differential Evolution Genetic Algorithm (DEGA) Sotirios Vardakos a,⇑, Marte Gutierrez b, Caichu Xia c a

Parsons Brinckerhoff Geotechnical and Tunneling, New York City, NY, USA Division of Engineering, Colorado School of Mines, Golden, CO, USA c Department of Geotechnical Engineering, Tongji University, Shanghai, PR China b

a r t i c l e

i n f o

Article history: Received 2 October 2010 Received in revised form 28 September 2011 Accepted 6 October 2011 Available online 9 November 2011 Keywords: Back-analysis Differential evolution Genetic algorithm Finite difference

a b s t r a c t The paper investigates the potential use of the Differential Evolution Genetic Algorithm (DEGA) in the back-analysis of tunnel response in order to obtain improved estimates of model parameters by matching model prediction with monitored response. DEGA is implemented in the finite difference code FLAC using FLAC’s built-in language FISH. The performance of DEGA in back-analyzing tunnel response is demonstrated using two idealized tunnel geometries and synthetic monitoring data, and a real case involving a section of the Heshang Road Tunnel. The use of DEGA in back-analysis of tunnel response is analyzed in terms of the stability and efficiency of the procedure under highly non-linear problems, the sensitivity of the solution to the initial trial assumption and the sensitivity of the solution to the monitoring data. It is shown that DEGA exhibits excellent convergence and repeatability characteristics, and that the convergence is independent of the initial values of the parameters. Another strength of DEGA comes from its ability to handle highly nonlinear objective functions with potentially multiple optima resulting from the use of incomplete and relative monitoring data, and heterogeneous monitoring data from different sensors. The paper recommends DEGA as a very viable technique in the back analysis and computational modeling of tunnel response. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Computational models are powerful tools in the performancebased engineering analysis and design of geotechnical structures. They can provide detailed evaluation of response in terms of displacements, loads on structural support, failure and damage. The use of computational models can lead to more rational designs for geotechnical projects than provided by empirical design procedures. However, the key problem in using computational models is the paucity of reliable information on site geology, soil and rock properties and in situ stresses. It has been suggested that data gathering and input may be the ‘‘Achilles heel’’ of computational modeling (Brown et al., 2002). Geological, geophysical and geotechnical investigations are time consuming and expensive, and are carried out extensively only for very important projects. Even these are inevitably incomplete. Without adequate and reliable input data, predictions from computational models provide very little help in the design of geotechnical structures. Often, predicted response does not correspond to observed performance. Due to the difficulty in making accurate predictions of response, geotechnical design

⇑ Corresponding author. Tel.: +1 646 709 7555; fax: +1 312 782 1684. E-mail address: [email protected] (S. Vardakos). 0886-7798/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.tust.2011.10.003

tends to be conservative, although occasionally unconservative designs have resulted in failure. Field monitoring can provide the means to verify predicted response, and to check the validity of values of material properties and loads used in design. Direct measurements of field response from monitoring can provide faster and more economical means of determining model parameters and for improving the reliability of model predictions. The determination of parameters required in computational models using monitored data is referred to by different terms, including back-analysis, parameter or system identification, inverse analysis and model updating. The term back-analysis is used in this paper. The use of formalized backanalysis procedures to integrate field monitoring data in computational models is still limited in many geotechnical engineering applications. Even in cases where field data are available, back-calculation of model parameters from field data are often carried out manually and unsystematically. Back-analysis is particularly suited to underground constructions such as tunneling. Underground constructions provide unique situations where more information about ground conditions and response become available as the construction progresses. In contrast to forward modeling where response is predicted given the required input data (e.g., material parameters and loads), back-analysis involves the determination of the input parameters

110

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

given the response (i.e., from field measurements). Predicted response from computer simulation using an initial set of model parameters can be compared with observations. If predicted and observed response deviate, input data are iteratively revised. The process of adjusting model parameters so that the model objectively matches observed response during some historical period can be accomplished by back-analysis. Iterative (and sometimes manual trial-and-error) procedures are often used since direct inversion is only possible for very simple analytical models. The iterative process is essentially an optimization problem where the objective is to minimize (in a least squares sense) the differences between predicted and monitored response. A crucial component of back-analysis is an algorithm to find a set of input parameters that will minimize the difference between predicted and measured performance (e.g., in terms of deformations, stresses or tunnel support loads). Methods of back-analysis can be broadly classified as direct and gradient based optimization techniques. Direct optimization methods exploit the vector relation between two successive solutions, and perform linear combinations of sequential solutions to find the local optimum. In gradient-based methods, the infinitesimal change of the solution path and the corresponding gradients are used to control the solution processes. Usually such methods may offer higher solution precision at the expense of the time required to calculate the first or second order derivatives. Brief descriptions of the available back-analysis procedures are given below. The present paper focuses on the application of Differential Evolution Genetic Algorithm (DEGA) in modeling of tunneling using a commercially available computer code. DEGA is a heuristics-based global search algorithm belonging to a wide class of Genetic Algorithms that have been used in a wide range of optimization problems in a variety of applications. A description of DEGA and its implementation in the commercially available computer code FLAC (Itasca, 2005) are presented. The use of DEGA in back-analysis of tunnel response is analyzed in terms of the stability and efficiency of the code under highly non-linear circumstances, the sensitivity of the solution to the initial trial assumption and the sensitivity of the solution to the monitoring data. The application of the proposed back-analysis procedure is demonstrated using idealized tunnel geometries in conjunction with synthetic and model-generated monitoring data. Finally, the proposed procedure is applied to the combined back- and forward-modeling of the Heshang Highway Tunnel project in China. 2. Back-analysis procedures In back-analysis, an algorithm handles the function evaluations so that the predictions of a numerical simulation will ultimately match the measured performance. This involves a process of minimization of the differences between predictions and measurements expressed as error through an objective function. An example of an error function is:

e¼

m X ðui  umi Þ2

ð1Þ

i¼1

where e is the error, ui is the ith predicted value of performance, umi is the corresponding ith value of measured performance, and m is the number of observations. A normalized version of the error function is given as:

e¼

m X i¼1

(

 2 ) ui  umi wi umi

ð2Þ

where wi is the weighting factor that can be applied for each measurement. This factor can be related with the reliability and quality

of the various monitoring data. In Geotechnical Engineering, computational models describe highly nonlinear processes where model parameters are dependent on predicted response, and vice versa. Thus, Eqs. (1) and (2) are generally highly nonlinear functions of the unknown parameters that cannot be expressed analytically. An extensive review of various methods for back-analysis of tunneling and other geotechnical problems can be found in Vardakos (2007). Back-analysis procedures fall into two distinct categories: (1) direct and (2) gradient-based optimization techniques. Direct optimization includes methods such as the Univariate method, the Pattern Search, Powell’s algorithm and the ordinary Simplex method. Generally, direct optimization methods are easier to implement since they do not require gradient calculations of the error function. Most often, after a search direction has been calculated in the search space, a one-dimensional optimization technique has to be implemented so that the new vector is calculated properly along the calculated search direction. These methods can be employed for well-posed geotechnical problems, for example problems involving elastic materials with one to three parameters to be back analyzed. Their precision is also limited since they often have difficulty in closely tracking the optimal design vector when the search has reached the general area close to the optimum. Examples of gradient-based methods are the Steepest Descent, the Conjugate Gradient, the Newton, and the Quasi-Newton methods. The main advantage of gradient-based methods is the property of quadratic convergence which accelerates the solution progress. When constraints are required for the parameters, then either a penalty procedure must be superimposed on the optimization process, or a more elaborate method needs to be used. The gradient-based methods are powerful algorithms but present some challenges. The need to evaluate derivatives makes the use of gradient methods difficult in iterative numerical analysis. This could be dealt with when there are only a few unknowns (i.e., two or three parameters) involved, however, when the number of unknowns increases then the process becomes not only computationally expensive but sometimes impossible. Many direct and gradient methods have been applied successfully in back-analysis of ground deformations for various geotechnical problems. Cividini et al. (1981, 1983) give an insightful review on the use of back-analysis in Geotechnical Engineering, including principles and algorithms, with examples of both direct and inversion methods. They presented probabilistic-based backanalysis procedure which shares the concept presented by Eykhoff (1974) in parameter identification. The importance of the error involved in the measurements in the back analyses was also extensively discussed. Gioda and Sakurai (1987) also present a survey of back-analysis methods and principles with specific reference to tunneling problems. Their review and examples involved both deterministic and probabilistic approaches. A displacement-based back-analysis methodology can be found in Sakurai and Abe (1981) and Sakurai and Takeuchi (1983). This methodology yields the estimates on both the initial in situ stresses and Young’s modulus of the rock mass by assuming the rock as linearly elastic and isotropic. Contributions to the probabilistic-based methods in back-analysis for tunneling are also presented in Ledesma et al. (1996) and Gens et al. (1996), which describe a minimization procedure along with reliability estimates of the final calculated parameters in conjunction with the finite element method. The use of the boundary control method along with a local search algorithm to perform back-analysis is suggested in Swoboda et al. (1999). This method was later improved and promoted by Xiang et al. (2002, 2003). The numerical results revealed that this is a stable and fast-converging algorithm under certain circumstances. The limitations of local search algorithms in geotechnical parameter identification are discussed by Vardakos (2007). A

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

displacement-based back-analysis method formulated as a combination of a neural network, an evolutionary calculation, and numerical techniques was proposed by Feng et al. (2000). A model for the back-analysis of a shallow tunnel using the conjugate gradient method was presented by Chi et al. (2001). A method for displacement-based back analysis using a neural network and a genetic algorithm is described in Deng and Lee (2001). An interesting application of back-analysis of in situ stresses based on small flat jack measurements was performed by de Mello-Franco et al. (2002) in the case of a mine. Lecampion et al. (2002), performed back analysis of constitutive parameters for an elasto-viscoplastic constitutive model used in a lined tunnel based on displacement measurements. Back-analysis using neural networks (NN) was presented by Pichler et al. (2003). Their method utilizes an artificial neural network (ANN) which is trained to approximate the results of finite element simulations. A genetic algorithm (GA) uses the trained NN to provide an estimate of optimal model parameters. In Deng and Nguyen-Minh (2003), a back-analysis method based on minimization of error on the virtual work principle is presented. Feng and An (2004) suggested the integration of an evolutionary neural network and finite element analysis using a genetic algorithm for the problem of a soft rock replacement scheme for a large cavern excavated in alternating hard and soft rock strata. The method of neural networks in back-analysis has also been used by Chua and Goh (2005). Finno and Calvello (2005) performed backanalysis of braced excavations using a maximum likelihood type of objective function and local search optimization. Yang et al. (2006) presented a genetic algorithm approach to identify the value and orientation of boundary stress conditions in deep tunneling using tunnel wall displacements.

3. Back-analysis using the Differential Evolution Genetic Algorithm (DEGA) DEGA belongs to a family of global optimization procedures that aims to find a globally optimum solution (if there is one) to an optimization problem. Excellent reviews of global optimization methods can be found in Horst and Pardalos (1995) and Pardalos and Romeijn (2002). Some of these methods include dynamic programming, branch and bound methods, response surface methods, annealing methods, and genetic algorithms. Of those methods, the last two are very strong candidates for back-analysis in geotechnical engineering and both are based on emulation of physical processes. Recently, there has been a growing interest in the use of simulated annealing (SA) and genetic algorithm as global optimization procedures for geotechnical engineering applications. Annealing and genetic algorithm-based approaches are heuristic and stochastic in nature, thus, they do not behave in ‘‘greedy’’ fashion like gradient or pattern search direction-based methods. As will be shown below, implementation of heuristic methods can be advantageous in some geotechnical back-analysis problems especially when ‘‘a-priori’’ information may not be available or when it is unreliable. Genetic Algorithms (GA) have been used in geotechnical backanalysis by several researchers. An example is the work of Simpson and Priest (1993). According to Rao (1996), mixed continuous and discrete variables may be searched by GA during optimization in a discontinuous or non-convex search space. The development of GA is attributed to Holland (1992) and the core of the method rests on the Darwinian theory of the survival of the fittest. The main elements of a genetic algorithm can be summarized as follows: (1) The procedure starts from a population of trial vectors, instead of a single point. For an n-variable problem, the usual population size is 2n–4n. For a highly non-linear objective function, higher sizes may be required. (2) There is no gradient or pattern search

111

direction exploitation. (3) Each variable resembles a chromosome in genetics. (4) The objective function value is the equivalent of fitness in genetics. In minimization problems, a vector corresponding to a very low function value is a strong candidate and may survive in the future. (5) The process is based on repeated generations of new population vectors. The trial vectors of these populations are the result of randomized parental selection, crossover and mutation processes. They are essentially the offspring of previous trial vectors. A GA typically handles the optimization process using a binary system with string lengths to represent the design variables and subsequently the whole design vector. This may create programming difficulties and inefficiencies in problems of continuous variable optimization where floating point values are used. Evolution strategies deal with this issue by bypassing the requirement of binary system usage and also make the algorithm more user-friendly and transparent. A type of evolutionary strategy called the Differential Evolution (DE) was proposed by Price and Storn (1997) and Storn and Price (1997). As a heuristic and stochastic-based method, DE is based on the experimental behavior of a function. DE tests the fitness of successive values of a trial vectors instead of using a gradient to drive the direction of the optimization. The key to DE is on how it performs and organizes the determination of the trial vectors. With increasing computational power and advances in computational techniques, it is now more prudent to directly test for the fitness of trial vectors independently from pattern directions rather than from gradients whose existence may be questionable. Other areas of application of the DE include locating earthquake hypocenters (Ruzek and Kvasnicka, 2001), water resources optimization (Reed and Yamaguchi, 2004), and estimation of rock fracture sizes (Decker and Mauldon, 2006). The original concept presented herein was developed by Vardakos (2007) as part of a more extensive back-analysis toolbox for tunneling which involved the application of local optimization methods such as the multivariate Newton–Raphson method and global optimization methods such as the SA and the DEGA. This comparative study outlined the benefits of global search algorithms as in contrast to local search methods and presented a rapidly deployable toolbox which can be used with a widely available numerical code. Su et al. (2008) have presented an application of the DE algorithm coupled with a numerical code for identification of the constitutive model type and its parameters, using convergence measurements from tunneling. The current paper discusses in detail the principles, structure and implementation of the method, the efficient use of parametric constraints, and computational performance of the DEGA. Careful back-analyses studies are presented using both idealized tunnel problems and an actual case study from a construction project in China in order to demonstrate the efficiency, convergence characteristics and uniqueness of solutions obtained via the use of DEGA in tunneling problems. DEGA combines the concepts of Differential Evolution and Genetic Algorithms to handle optimization of non-binary valued nonlinear functions. Its implementation in parameter identification is outlined in Fig. 1. DEGA uses two arrays to store a population of NP, D-dimensional vectors of parameters that are being back-calculated (D = n = number of parameters). The parameter vector may consist of heterogeneous data from different types of monitoring (e.g., displacements, strain, loads on structures, pore pressures, etc.). The first array contains the primary values, which contains the present vector population, and the secondary array stores sequentially the products for the next generation. The algorithm starts by filling the primary array with NP vectors with randomly generated parameter values. The initial random generation should satisfy the constraints on the parameters. The primary array is also called as the ‘‘trial vector’’ since it contains NP vectors that will

112

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

Fig. 1. Flowchart for the Differential Evolution Genetic Algorithm (DEGA).

later be tried for fitness. Each of these individual randomly-generated vectors Xi is considered sequentially for genetic operations. For each of the chosen vectors, three other vectors XA, XB, XC are randomly chosen from the remaining vectors of the primary array. A mutation can then be performed by using the following relation to generate a new trial vector:

Xm 1 ¼ XA þ FðXB  XC Þ Xm 1

ð3Þ

where is the new mutant vector, and F is a scaling factor in the range 0 < F 6 1:2. According to Price and Storn (1997), the optimum value of F is in the range 0.4–1. One of the main issues of successful optimization strategies is the implementation of constraints. As is the case in geotechnical engineering problems, this is a required step in order to constrain the search procedure only

within a feasible and meaningful space. Often a-priori information from geotechnical investigation and laboratory testing programs will be available and can be used in order to get a first overview of the parameters involved. A modification of the main DEGA is applied in this study whereby the mutant vector Xm 1 is checked for constraint violation. Even though the initially generated XA, XB, XC vectors are feasible, their linear combination may be violating the constraints. If any parameter constraint is violated, then the sampled vectors XA, XB, XC are discarded and a new sampling is performed until a feasible vector Xm 1 is found. This intervention is seamlessly integrated in the algorithm and also keeps the true heuristic process of the problem intact as there is no change in the objective function since no additional penalty functions are used.

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

At this stage the crossover takes place. A random integer number randint(i) in the range [1, n] is generated. For each parameter j = 1, . . ., n, a random number randnum(j) is generated in the space [0, 1]. Then a new vector is created from the original Xi parent and the mutant vector using the crossover criterion:

( x0i;j

¼

xm i;j

if randnum 6 XR or randintðiÞ ¼ j

xi;j

if randnum > XR and randintðiÞ – j

113

constraints, the initial random filling of the primary array and definitions of the mutation, crossover and genetic comparison operations. Finally, the sixth block implements the DEGA iterative solution, which calls the previously stated functions. 4. Application of DEGA to back-analysis of tunnel response

ð4Þ

where XR is a crossover rate in the range [0, 1]. The crossover scheme essentially means that if randnum > XR, the new ith trial vector will receive the j parameter from the parent vector, otherwise the parameter will be obtained from the mutant vector Xm 1. In this way if XR = 1, then every trial vector will be obtained from the mutant vector, or if XR = 0, then all except for one parameter will be called from the parent trial vector. The new vector X0i;j is tested against Xi for fitness. Thus two function evaluations occur at this stage. For minimization problems, the vector corresponding to the lower value (fittest candidate) is entered in the secondary array. The same procedure is followed until all vectors of the original primary array are processed and an equal size secondary array has been formed. At this stage the secondary array values are transferred and update the primary array while the secondary array is purged. This marks the end of one generation. Obviously many generations are required for convergence. The above steps are repeated until a maximum number of generations is reached. When the algorithm converges to the global optimum based on the criterion, then all vectors of the primary array theoretically become equal. The flowchart for one generation of the method is shown in Fig. 1. The finite difference code FLAC (version 5.0) developed by Itasca (2005) is chosen for the implementation of DEGA. The proprietary built-in programming language FISH in FLAC, allows for the full implementation of DEGA in FLAC. As described in Vardakos (2007), this programming language can be successfully used for a wide variety of back-analysis implementations. The DEGA algorithm was implemented in FLAC with a six-block program structure. In typical cases, where the finite difference grid coordinates remain constant, only the input parameters will be changed before each FLAC solution. In this programming scheme, one FLAC solution (from stress initialization to end of construction) results in one objective function evaluation. The implementation of DEGA as an internal routine in FLAC, in comparison to an external implementation, has the important advantage of avoiding unnecessary and time consuming input/output and data exchange operations. The first block includes only the geometrical setup of the numerical model (no boundary conditions yet) while the second is used for the initialization of the genetic algorithm parameters and to predefine the sizes of the various needed arrays (tables to be filled in later during the DEGA process) via the use of function definitions. The third block includes the definitions for the locations and type of data to be monitored (which can be grid point or zone coordinates depending on the type of measurement), the values of the actual measurements, and the definition of the normalized error function as in Eq. (2). It is noted that all measurements from any number of stages can be defined at this block. The fourth block includes the actual run of the FLAC model via a function statement. Essentially this step is the objective function evaluation. All involved parameters are referred to inside this function with their generic parametric names (i.e., p1, p2, etc.) It is necessary to ensure at this stage that this function (i.e., the FLAC trial solution) be repeatable without leaving any unwanted old solution data in the memory of FLAC. All displacements, velocities, stress regime and plasticity states should be reset at the beginning of this function. The fifth block includes the definitions for the parametric

The applicability and performance of DEGA in the back-analysis of tunnel response is first investigated using numerical models of simple tunnel geometries with ‘‘synthetic monitoring data.’’ The synthetic monitoring data are nothing more than the model response at selected points in the model, which were generated with given model parameters. The objective of the back-analysis is, based on initial guesses of the values, to try to identify the model parameters that were used to generate the synthetic monitoring data. The use of DEGA is illustrated using: (1) a deep circular tunnel, and (2) a shallow circular tunnel. Both tunnels are assumed to be constructed in elastoplastic ground under anisotropic in situ stress condition. Plane strain conditions are assumed. The plane strain condition is still today the most often used numerical approximation. The main advantage is the shorter execution speed compared to three-dimensional models, which leads to timely solutions for a wide range of geotechnical scenarios of a particular problem. Nevertheless, DEGA can be also be equally implemented in 3D computational models as demonstrated by Su et al. (2008). 4.1. Example 1 – deep circular tunnel The first validation model is shown in Fig. 2 and consists of a circular tunnel with a radius R = 5.0 m. Due to symmetry, only a quarter of the model is discretized and the size of the model is assumed to extend 50 m radially from the center of the tunnel. Support is composed of a 20 cm thick shotcrete liner. Threedimensional loading effects and delayed ground-support interaction are simulated by letting the tunnel to relax (i.e., unload) to 50% of its initial in situ stress before the support is installed. This amount of relaxation is taken arbitrarily for this example and practically corresponds to support installed at some small distance from the tunnel face. The tunnel is excavated under anisotropic

Fig. 2. Deep circular tunnel model with monitoring locations.

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

Parameter

True optimum

Min

Max

Back-analysis

E (GPa)

6.0 0.25 14 6.0 0.8 30

5 0.22 12 5 0.6 26

7.5 0.3 15 8 1.5 35

5896 0.25 13.8 5.9 0.8 29.8

c0 (MPa) /0 (°)

Cohesion (MPa)

Frequency

Cohesion (MPa)

Frequency

Frequency

Frequency

Horizontal stress (MPa)

Poisson ratio

Poisson ratio

Vertical stress (MPa)

Horizontal stress (MPa)

Young's Modulus (MPa)

Frequency

Frequency

Young's Modulus (MPa)

Frequency

Frequency

Frequency

in situ stress conditions (i.e., rh < rv, where rh = horizontal in situ stress, rv = vertical in situ stress). As the tunnel is assumed to be deep enough, the variation of in situ stresses with depth along the vertical boundaries is neglected. The ground is assumed to behave elastoplastically following a Mohr–Coulomb failure criterion. The tunnel response was simulated using a known set of input model parameters as shown in Table 1, in order to create the ‘‘synthetic monitoring data.’’ These parameters where the elastic modulus E, the Poisson’s ratio m, the horizontal and vertical stresses rh and rv, and the strength parameters c0 and /0 . Relative displacements and lining loads were used as input for the back-analysis. To closely simulate a real case, the hypothesis is made that monitoring commences at the time of installation of the support liner (i.e., after 50% relaxation has occurred). Instrumentation is assumed to be three multipoint borehole extensometers. There is one 3-point extensometer at the wall, one 3-point at the 2 o’clock tunnel wall position, and one 6-point MPBX at the crown. The locations of the instrumentation are shown in Fig. 2. The axial loads developed in the shotcrete lining, were monitored at springline, crown and at 45° locations from the horizontal. The

Frequency

m ry (MPa) rx (MPa)

inclusion of lining loads allows the function to become more sensitive to changes of ground strength parameters and ground elastic modulus. All measurements were assumed to have the same weight factor of wi = 1. In the DEGA back-analysis, the number of unknowns is n = 6 and the population size was taken as NP = 60. The mutation scaling factor was taken as F = 0.5 and the crossover rate XR = 0.7. A total of 70 generations was executed during the back-analysis. Calculation time was clocked at about 5 min per generation in double precision using one core of a PC with a 2.67 GHz Intel Xeon X5550 Quad Core processors. Fig. 3 presents the results from one trial of the back-analysis problem using DEGA. The results are shown as frequency plots of the populations before and after the back-analysis for each of the six parameters. As can be observed, the solution shows convergence very close to the theoretical global optimum solution of the problem. At a higher number of generations, all parameters should converge to a more stable solution, but even at 70 generations the results are very close to the true value revealing very good performance of the back analysis procedure. In fact, during the back-analysis it was observed that useful results can be obtained as early as 20 generations, by examining the frequency distribution plots. The same analysis was repeated five more times, starting always from a different randomly generated population, to test the behavior of the code and check the repeatability of results. In all trials, the primary population array converged towards the true global optimum. Fig. 4 shows a comparison between the synthetic measured data and the predictions of the DEGA backanalysis. The final prediction was based on mean values of each parameter from the primary array population at 70 generations. As expected there is excellent agreement between the predictions and the synthetic monitoring data.

Frequency

Table 1 Back-analyzed parameters for deep the tunnel model.

Vertical stress (MPa)

Frequency

114

Friction angle (deg)

Friction angle (deg)

Fig. 3. Results from DEGA back-analysis for deep circular tunnel model. First and third column shows the initial random populations, and the second and fourth columns shows the population distribution after 70 generations.

115

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123 0

Extensometer anchorage (m)

1

MPBX 1

2 3 4 5 Monitoring

6

Back-analysis 7 0

1

2

3

4

5

6

Displacement (mm) Fig. 5. Shallow tunnel model with monitoring locations. 0

Extensometer anchorage (m)

MPBX 2 1

2

3

4

Monitoring Back-analysis

5 0

2

4

6

8

10

12

14

16

18

20

16

18

20

Displacement (mm)

Extensometer anchorage (m)

0

MPBX 3

1

2

Monitoring Back-analysis 3 0

2

4

6

8

10

12

14

Displacement (mm) Fig. 4. Comparison between measured extensometer displacements and backanalysis results for the deep circular tunnel model.

check the ability of the DEGA algorithm to obtain four parameters: the elastic modulus E, the in situ Ko ratio, and strength properties c0 and /0 by using two sets of deformation measurements. Lining stresses or strains are not employed in this back-analysis. The tunnel has a radius of 8.0 m and its crown is located at a depth of 16.0 m from the ground surface. Due to symmetry only half of the problem is modeled. Support was assumed to be installed after 40% relaxation of the initial in situ stress, and is composed of a 15 cm thick shotcrete. The ground is assumed to be homogenous and follows the Mohr– Coulomb failure criterion. The unit weight of the surrounding ground is assumed constant at 24 kN/m3 and the horizontal stresses vary linearly over depth with a Ko stress ratio which is assumed to be constant over the depth of the model. The Poisson’s ratio m was assumed to be constant for all trials and equal to 0.26. Horizontal displacements were not allowed for the two vertical sides and the bottom was fixed in the vertical direction. During back-analysis, it was assumed that the overburden displacements are monitored using a vertical 6-point extensometer above the crown and horizontal displacements are monitored with a 24 m long inclinometer installed at 10 m distance from the tunnel centerline. The monitored displacements along the extensometer AA0 were referenced to the MPBX anchor head at point A. Similarly all inclinometer displacements were relative to the bottom of the borehole. The monitoring is assumed to commence at a pre-construction stage, thus the majority of significant strains in the ground should be captured by the equipment. For this example only the tunnel relaxation before installation of the support is accounted for in the model. Synthetic monitoring data were generated in a forward model of the tunnel using the properties shown in Table 2. All measurements were assumed to have the same weight factor of wi = 1. In the DEGA back-analysis, the number of unknowns is n = 4 and the population size was taken as NP = 40. The mutation scaling factor was taken as F = 0.5 and the crossover rate XR = 0.7. A total of 70 generations were executed during the

4.2. Example 2 – shallow circular tunnel The second problem used to test the performance of DEGA is the back-analysis of the shallow tunnel problem shown in Fig. 5. This type of tunneling often occurs in urban environments where ground surface deformations should be kept to a minimum. In typical cases for shallow tunnels, monitoring most often starts before the excavation via inclinometers, multipoint extensometers installed from the ground surface, or with surface surveying. Shallow excavations through soft soils often suffer from ground or otherwise volume loss. The main purpose of this model is to

Table 2 Back-analyzed parameters for the shallow tunnel model. Parameter

True optimum

Min

Max

Back-analysis

E (kPa) Ko c0 (kPa) /0 (°) c0 (kN/m3)

30,000 0.5 8 25 24 – constant 0.26 – constant

20,000 0.3 2 20

45,000 0.6 20 35

30,015 0.5 7.99 25

m

116

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

Fig. 6. Results from DEGA back-analysis for shallow tunnel model. Left column shows the initial random populations, and the right column shows the converged values after 70 generations.

back-analysis. Calculation time was clocked at about 6 mins per generation in double precision using one core of a PC with a 2.67 Ghz Intel Xeon X5550 Quad Core processors. The results of the back analysis, shown in Figs. 6 and 7 demonstrate good convergence behavior of the back-analysis procedure. Fig. 6 shows the evolution of the frequency distributions as they become significantly narrower after 70 generations of the mutation of the parameters. All the components of the primary array have practically converged to the same optimal vector which is the global minimum solution after 70 iterations. Fig. 7 presents a comparison between the synthetic monitoring data and the predictions at the end of 70th generation. As in the previous example, the predictions are based on mean parameter values from the distribution of 70th generation of the primary array population. As can be seen, the predicted responses using the back-analyzed parameters are nearly identical to those of the synthetic monitoring data.

5. Application DEGA to the back-analysis of the Heshang Tunnel To illustrate its applicability to real problems, the DEGA-based back-analysis procedure was applied to the analysis of the response during construction of the Heshang Highway Tunnel. This tunnel is located in South East China at the Fujian province, about 5 km southeast of the Fuzhou City, and is a part of the transportation system between the local airport and Fuzhou city. The tunnel construction was completed on August 30, 2006. It is a twin tunnel project and approximately 450 m long. The geology along the tunnel is shown in Fig. 8, and a typical tunnel cross-section and construction stages are shown in Fig. 9. Each of the two tunnels is 11.5 m high and 15 m wide. The tunnel passes through highly weathered volcanic material. Towards the Northwest portal the tunnel penetrates through weathered tuff lava and residual loam. In the central section, the majority of the rock mass is weathered

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123 0

Extensometer anchorage (m)

2 MPBX 4 Monitoring 6

Back-analysis

8 10 12 14 16 0

2

4

6

8

10 12 14 16 18 20 22 24

Displacement (mm) 24 22

Inclinometer

20

Inclinometer height (m)

18 16 14 12 10 Tunnel side 8 6 Monitoring

4

Back-analysis

2 0 -20

-15

-10

-5

5

0

10

Displacement (mm) Fig. 7. Comparison between measurements and back-analysis results for shallow circular tunnel model.

117

tuff lava, while at the Southeast portal the rock mass quality drops again, with moderate and highly weathered tuff lava. Due to the poor quality of the rock mass, and also due to the lateral proximity of the two tunnels, a wide array of instrumentation methods was employed during the construction. Selected sections of the tunnel were fully instrumented using extensometer, internal convergence and surface displacements (see Fig. 10). The goal of the back-analysis is to determine model parameters using monitoring data from the earlier excavation stages and use the parameters for forward modeling of subsequent excavations. Thus, material parameters determined from the back-analysis of the left tunnel, which was constructed first, were then used in the prediction (or forward modeling) of the response of the right tunnel (i.e., experience gained from the behavior from one portal region is applied to predictions during excavation of the other portal). The combined back- and forward modeling was performed for all the stages of tunnel construction at Station K6 + 300. It is noted that based on the available information, no lateral differences in rock types were encountered in the left and right tunnels at Station K6 + 300. Thus, material parameters back-calculated from the left tunnel should be directly applicable for the right tunnel. The sequential excavation of the Heshang tunnel was designed in accordance with anticipated ground conditions, and different excavation sequences were used for the left and right tunnels. The multi-staged excavation of the left tunnel is a typical sequential type, with two side drifts, a top and a bottom core (Stages 1, 2, 3 and 5). The right tunnel is excavated by a top heading, two bench sections and followed by an invert (Stages 4, 6, 7 and 8). Due to the highly weathered type of the rock mass, a variety of primary support systems and ground improvement techniques was used at the Heshang tunnel. More specifically, the ground around the tunnel is pre-reinforced using pipe spilling composed of 50 mm diameter grouted steel pipes of 5.0 mm thickness and 5.0 m length. The circular interval of the pipes is 300 mm and the longitudinal spacing is 2.5 m. Spilling was also used for the pillar wall between the two tunnels. The primary support system consisted of 270 mm thick shotcrete. The shotcrete was reinforced with welded wire mesh and U-type steel beams were installed at every 0.5 m. The tunnel

Fig. 8. Geological profile along the length of Heshang Tunnel.

118

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

walls and roof were also reinforced with hollow rock bolts of 25 mm outside diameter and 4.0 m length. The bolts were installed at a 0.8 m  0.8 m grid. Similar to the example circular tunnels described above, the ground is modeled by a Mohr–Coulomb failure criterion. The prereinforcement was simulated using finite difference zones, instead of using discrete structural elements. The equivalent properties for the reinforced zone are calculated using a simple homogenization scheme, where the equivalent strength and deformability of the reinforced ground are estimated using weighted volume average contributions of the respective host rock and reinforcing material. For instance, the uniaxial compressive strength of the spilling region is estimated by summing the products of the strength of each participating material (steel pipe, grout, rock) by their corresponding cross sectional area and by dividing the sum with the total spilling section area. It was also assumed that the equivalent properties of the reinforced zone are known and remain constant during the back-analysis. This assumption is generally valid since the

Table 3 Back-analyzed parameters for the Heshang twin-tunnels. Parameter

Min

Max

Back-analysis

E (MPa)

100 0.3 16 0.02 24

1500 0.45 26 0.2 40

428 0.3 23.1 0.023 28

m c0 (kN/m3) c0 (MPa) /0 (°)

reinforced ground has more predictable properties than the native ground. The unknown properties that were back-analyzed are the rock mass elastic properties E and m, the cohesion c0 and friction angle /0 , and the average unit weight of the rock c. The back-analysis is performed using monitoring data obtained from Stages 2 and 4 in the excavation of the left tunnel. These include extensometer displacements measured at KO1 and KO2. It is noted that Stages

Fig. 9. Tunnel cross-sections and construction stages at Station K6 + 300 of Heshang Tunnel.

Fig. 10. Locations and details of extensometer measurements at Station K6 + 300 of Heshang Tunnel.

Fig. 11. FLAC model for the Heshang twin tunnel.

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

3 (left tunnel top core) and 4 (right tunnel top heading) had already been completed before the second set of measurements was obtained. A close-up of the finite difference grid used in region around the tunnels is shown in Fig. 11. During the finite difference grid generation, grid points corresponding to measurements points were relocated to conform to the actual monitoring locations. Roller boundaries are used for the left, bottom and right bound-

119

aries of the model. Initial stresses due to self weight of the materials are generated first, and then the excavation sequence begins as shown in Fig. 9. Due to the use of a two-dimensional model, estimates have to be made on how much tunnel deformations have occurred before support elements are installed. These deformations are based on the convergence measurements made during construction and the actual schedule of the support installation.

Fig. 12. Frequency plots from the population of the primary array at the end of the 1st and the 50th generations of the DEGA back-analysis.

120

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

Table 3 lists the estimated range of values for the model parameters involved in the back-analysis and the final back-analyzed values. The estimated ranges of parameters values are based on the available geotechnical information. As described earlier, there is no need to estimate a starting point, since the original population is randomly generated within its estimated range. In the backanalysis, the mutation scaling factor was taken as F = 0.5, the crossover rate was XR = 0.7, and the program was setup for 50 generations. Due to the complicated geometry and large number of grid elements, calculation time was clocked at about 3.5 h per generation in double precision using one core of a PC with a 2.67 GHz Intel Xeon X5550 Quad Core processors. The results of the back-analysis are shown in Figs. 12–14 and show interesting observations on the behavior of DEGA during mutations of the solution. Fig. 12 shows the distributions of the five back-analyzed parameters at the end of the first generation illustrating the significant variance in the histograms. As the solution progresses, the variance decreases and at the end of the 50th generation, the solution has nearly converged to the global optimum set of parameters. Fig. 13 shows a comparison between monitored displacements of the KO1 and KO2 extensometer and results from back-analysis. The comparisons are shown for the two construction stages described previously. As can be seen, good agreement between measured and back-analyzed results were obtained for both construction stages. There is very good agreement for the extensometer KO1 for both stages, and a reasonable agreement for extensometer KO2 with some small discrepancies. In any case, the discrepancies between simulated and measured values do not exceed 1 cm. However, given the unknown variation of the exact nature of the geomaterials around the tunnels and the inevitable

assumption of a homogenous ground around both tunnels, the analysis is considered successful for most engineering intents. Fig. 14 shows a comparison for the forward modeling part between monitored extensometer data and predicted response at the end of Stage 4 of the right tunnel. Due to excavation of the right tunnel, single top heading large displacements tend to develop around the right tunnel and all extensometers are now registering detectable displacements. There is relatively good agreement for extensometer KO5, but generally the back-analysis resulted in a slight overestimation of the displacements developed at KO3 and KO4, possibly due to local variations in the quality of the rock mass which were not explicitly modeled here. The differences between predicted and measured response are in the range of 1 cm in most cases. The relatively good agreement between predicted and monitored response indicates the validity of the model parameters obtained from the back-analysis of the tunnel response at the earlier stages of construction. As shown in the plasticity plot in Fig. 15, the close proximity of the twin tunnel to the slope and the ground surface results in significant ground strength mobilization around both excavations. This also verifies the requirements for ground pre-support and densely installed initial support systems in the Heshang tunnels. The FLAC forward-analysis model predicts a maximum displacement of approximately 5 cm which develops at the upper pillar side of the left tunnel and at the crown of the right tunnel as shown in Fig. 16 and Fig. 17. The predicted plastic behavior serves also to demonstrate that DEGA is fairly immune to being locked in an elastic behavior mode during the back-analysis process. The latter is a problem with many pattern search or gradient-based methods. For example, if the upper bound constraints of the strength parameters

7

7 KO1 - Stage 2

5 4 3 2

Monitoring Back-analysis

1

KO1 - Stage 4

6

Extensometer anchorage (m)

Extensometer anchorage (m)

6

5 4 3 2 Monitoring 1

Back-analysis

0

0 0

1

2

3

4

0

Displacement (cm)

4

6

11

11 Monitoring

10

10

Back-analysis

9

Extensometer anchorage (m)

Extensometer anchorage (m)

2

Displacement (cm)

8 7

KO2 - stage 2

6 5 4 3 2

9

KO2 - Stage 4

8 7 6 5 4 3

Monitoring

2 Back-analysis

1

1

0

0 0

1

2

Displacement (cm)

3

4

0

2

4

6

Displacement (cm)

Fig. 13. Comparisons of extensometer displacements at Station K6 + 300 from monitoring and back-analysis. Top: Extensometer KO1. Bottom: Extensometer KO2.

121

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123 10

10 9

Monitoring

8

Extensometer anchorage (m)

Extensometer anchorage (m)

9

Back-analysis

7 6 5

KO3 - Stage 2

4 3 2 1

Monitoring

8

Back-analysis

7 6 5 KO3 - Stage 4

4 3 2 1 0

0 0

1

2

3

0

Displacement (cm)

6

8

Back-analysis

9

7 6 KO4 - Stage 2

5 4 3 2

Monitoring

8

Back-analysis

7 6 5

.

9

Monitoring

Extensometer anchorage (m)

Extensometer anchorage (m)

4

10

10

KO4 - Stage 4

4 3 2 1

1

0

0 0

1

2

3

0

Displacement (cm) 10

10

9

9

8

Monitoring

7

Back-analysis

6 5 KO5 - Stage 2

4

2

4

6

Displacement (cm)

Extensometer anchorage (m)

Extensometer anchorage (m)

2

Displacement (cm)

3 2 1

8

Monitoring

7

Back-analysis

6 K05 - Stage 4

5 4 3 2 1

0

0 0

1

2

Displacement (cm)

3

0

2

4

Displacement (cm)

Fig. 14. Predicted extensometer displacements above the right tunnel from forward modeling.

Plasticity Indicator * at yield in shear or vol. X elastic, at yield in past o at yield in tension

Fig. 15. Predicted plastic regions around the tunnels from forward modeling.

6

122

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123 Y-displacement contours -4.00E-02 -3.00E-02 -2.00E-02 -1.00E-02 0.00E+00 1.00E-02 2.00E-02

Fig. 16. Contours of predicted vertical displacements (in m) from forward analysis.

Displacement vectors scaled to max = 2.000E-01 max vector = 5.040E-02

Fig. 17. Total displacement vectors (in m) around the tunnels predicted after back-analysis.

are high enough, and the lower bound constraint for the stiffness is low enough, the solution could easily deviate and restrain itself into an iterative search in an elastic-only domain of the search space, thus not being able to differentiate between iterative changes of the plastic parameters.

response. For the Heshang Road Tunnel, there is also satisfactory agreement between model and monitored response in both the back and forward modeling stages. The small discrepancies, which in all cases were less than 1 cm, were attributed to rock mass nonuniformities which were not accounted for in the modeling. The study revealed the following strengths that give DEGA strong potential in back-analyzing tunnel response:

6. Conclusions The paper investigated the potential use of the Differential Evolution Genetic Algorithm (DEGA) in the back-analysis of tunnel response. The goal of back-analysis is to obtain improved estimates of model parameters by matching model prediction with monitored response. DEGA is a variant of a class of heuristic type global optimization Genetic Algorithms. DEGA was implemented in the finite difference code FLAC using FLAC’s built-in language FISH. The performance of DEGA in back-analyzing tunnel response using FLAC was demonstrated using two idealized tunnel geometries and synthetic monitoring data. The objective of the back-analysis is, based on initial guesses of the values, to try to identify the model parameters that were used to generate the synthetic monitoring data. DEGA was also used in a real case involving a section of the Heshang Road Tunnel. DEGA was used to generate improved estimates of the model parameters from monitoring data on earlier excavation stages of the tunnel. The back-calculated model parameters were then used in forward modeling to predict the tunnel response on subsequent excavation stages. DEGA showed excellent convergence and repeatability characteristics. The frequency distribution of model parameters became narrower with iteration and subsequent parameter mutations. Parametric studies also showed that the convergence was independent of the initial values of the parameters. For the case of the two idealized tunnel cases, model predictions using the back-calculated parameters were nearly identical with the synthetic monitoring

(1) DEGA does not require the modification of the main computational model and does not require the computation of complicated quantities such as the first or second order gradients of the response with respect to the model parameters (i.e., as in the case of gradient-base methods). DEGA can be implemented simply as a master program that controls the execution of the computational model. Thus, DEGA can be readily implemented in any commercially available computer code. (2) DEGA can be used in conjunction with incomplete monitoring history as is commonly encountered in tunnel constructions where monitoring typically commences only after certain amounts of stress release and deformations have already occurred in the rock mass. The amount of pre-deformation can play a significant role in back-analysis as they can result in a very irregular objective function with potentially multiple local optima. In contrast to other methods, DEGA has the ability to avoid getting stuck in local optima. (3) DEGA can be easily used in conjunction with heterogeneous monitoring data consisting of measurements from different types of sensors. This is due to the fact that the objective function can be easily defined using any combination of load, stress, deformation and strain measurements, or in terms of relative instead of absolute measurements. The only requirement is that the back-analysis program can access the same set of data from the computational model. The use of heterogeneous and relative measurements can influence the shape of

S. Vardakos et al. / Tunnelling and Underground Space Technology 28 (2012) 109–123

the objective function and subsequently the convergence of the back-analysis. DEGA is capable of handling irregular and highly nonlinear objective functions. For highly complex numerical analysis problems, the DEGAbased back-analysis requires relatively long calculation times, however, for more typical applications the proposed back-analysis execution is considered very efficient. Even with more complex numerical applications, the DEGA application in FLAC showed that it is possible to obtain useful preliminary results from the early stages of the differential evolution, by closely examining the trends of the frequency distribution plots of the parameters at the end of a small number of generations. In all cases, the repeatability of results greatly supports its use. Furthermore, with newer and faster processing capabilities, calculation time is becoming less of an issue. The advent of parallel computing on computers with multiple processors is expected to speed up the back-calculation times. Better programming can also improve the performance of computer codes as was obtained in version 6.0 of FLAC, which is 35% faster than version 5.0. In summary, based on the results presented above, it is concluded that DEGA is a very viable technique in the back analysis and computational modeling of tunnel response. Acknowledgements The results presented in this paper are part of the AMADEUS (Adaptive real-time geologic Mapping, Analysis and Design of Underground Space) Project funded by the US National Science Foundation under grant number CMS 324889. This support is gratefully acknowledged. References Brown, E., Cundall, P., Desai, C. Hoek, E., Kaiser, P. 2002. The great debate in rock mechanics: data input – the Achilles Heel of rock mass constitutive modeling. In: Panel Discussion at NARMS-TAC 2002 Conference, Toronto, Canada, July 7– 10; 2002. Chi, S.-Y., Chern, J.-C., Lin, C.-C., 2001. Optimized back-analysis for tunnelinginduced ground movement using equivalent ground loss model. Tunneling and Underground Space Technology 16, 159–165. Chua, C.G., Goh, A.T.C., 2005. Estimating wall deflections in deep excavations using Bayesian neural networks. Tunnelling and Underground Space Technology 20, 400–409. Cividini, A., Jurina, G., Gioda, G., 1981. Some aspects of characterization problems in geomechanics. International Journal of Rock Mechanics and Mining Sciences 18, 487–503. Cividini, A., Maier, G., Nappi, A., 1983. Parameter estimation of a static geotechnical model using a Bayes’ rule approach. International Journal of Rock Mechanics and Mining Sciences 20, 215–226. de Mello-Franco, J.A., Armelin, J.L., Santiago, J.A.F., Telles, J.C.F., Mansur, W.J., 2002. Determination of the natural stress state in a Brazilian rock mass by backanalyzing excavation measurements: a case study. International Journal of Rock Mechanics and Mining Sciences 39, 1005–1032. Decker, J., Mauldon, M. 2006. Determining size and shape of fractures from trace data using a differential evolution algorithm. In: Proceedings of the 41st US Symposium on Rock Mechanics, Golden, CO, Electronic Proceedings. Deng, D., Nguyen-Minh, D., 2003. Identification of rock mass properties in elastoplasticity. Computers and Geotechnics 30, 27–40. Deng, J.H., Lee, C.F., 2001. Displacement back analysis for a steep slope at the Three Gorges Project site. International Journal of Rock Mechanics and Mining Sciences 38, 259–268. Eykhoff, P., 1974. System Identification – Parameter and State Estimation. John Wiley and Sons, Chichester, England.

123

Feng, X.-T., An, H., 2004. Hybrid intelligent method optimization of a soft rock replacement scheme for a large cavern excavated in alternate hard and soft rock strata. International Journal of Rock Mechanics and Mining Sciences 41, 655– 667. Feng, X.-T., Zhang, Z., Sheng, Q., 2000. Estimating mechanical rock mass parameters relating to the Three Gorges Project permanent shiplock using an intelligent displacement back analysis method. International Journal of Rock Mechanics 37, 1039–1054. Finno, R.J., Calvello, M., 2005. Supported excavations: the observational method and inverse modeling. Journal of Geotechnical and Geoenvironmental Engineering 131 (7), 826–836. Gens, A., Ledesma, A., Alonso, E.E., 1996. Estimation of parameters in geotechnical back-analysis – II. Application to a tunnel excavation problem. Computers and Geotechnics 18 (1), 29–46. Gioda, G., Sakurai, S., 1987. Back analysis procedures for the interpretation of field measurements in geomechanics. International Journal for Numerical and Analytical Methods in Geomechanics 11, 555–583. Holland, J.H., 1992. Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence. MIT Press, Cambridge, MA. Horst, R., Pardalos, P.M., 1995. Handbook of Global Optimization – Nonconvex Optimization and Its Applications. Kluwer Academic Publishers, Dordrecht, The Netherlands. Itasca, 2005. FLAC – Fast Lagrangian Analysis of Continua. Itasca Consulting Group, Minneapolis, MN. Lecampion, B., Constantinescu, A., Nguyen-Minh, D., 2002. Parameter identification for lined tunnels in a viscoplastic medium. International Journal for Numerical and Analytical Methods in Geomechanics 26, 1191–1211. Ledesma, A., Gens, A., Alonso, E.E., 1996. Estimation of parameters in geotechnical back analysis – I. Maximum likelihood approach. Computers and Geotechnics 18 (1), 1–27. Pardalos, P.M., Romeijn, H.E., 2002. Handbook of Global Optimization. Nonconvex Optimization and Its Applications, vol. 2. Kluwer Academic Publishers, Dordrecht, The Netherlands. Pichler, B., Lackner, R., Mang, H.A., 2003. Back analysis of model parameters in geotechnical engineering by means of soft computing. International Journal for numerical Methods in Engineering 57, 1943–1978. Price, K., Storn, R., 1997. Differential evolution – a simple evolution strategy for fast optimization. Dr. Dobb’s Journal 18, 24–78. Rao, S.S., 1996. Engineering Optimization, Theory and Practice, third ed. John Wiley and Sons, New York. Reed, P.M., Yamaguchi, S. 2004. Simplifying the parameterization of real-coded evolutionary algorithms. In: Proceedings of the World Water and Environmental Resources Congress, Salt Lake City, UT, pp. 1–9. Ruzek, B., Kvasnicka, M., 2001. Differential evolution algorithm in the earthquake hypocenter location. Pure and Applied Geophysics 158, 667–693. Sakurai, S., Abe, S., 1981. Direct strain evaluation technique in construction of underground opening. In: Proceedings of the 22nd US Symposium on Rock Mechanics. Cambridge, MA, pp. 278–282. Sakurai, S., Takeuchi, K., 1983. Back analysis of measured displacements of tunnels. Rock Mechanics and Rock Engineering 16 (3), 173–180. Simpson, A.R., Priest, S.D., 1993. Application of genetic algorithms to optimization problems in geotechnics. Computers and Geotechnics 15 (1), 1–19. Storn, R., Price, K., 1997. Differential evolution – a simple yet efficient heuristic for global optimization over continuous spaces. Journal of Global Optimization 11, 341–359. Su, G.-S., Zhang, X.-F., Chen, G.-Q., Fu, X.-Y., 2008. Identification of structure and parameters of rheological constitutive model for rocks using differential evolution algorithm. Journal of Central South University of Technology 15, 25–28. Swoboda, G., Ichikawa, Y., Dong, Q., Zaki, M., 1999. Back analysis of large geotechnical problems. International Journal of Numerical and Analytical Methods in Geomechanics 23, 1455–1472. Vardakos, S. 2007. Back-Analysis Methods for Optimal Tunnel Design. Ph.D. Thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA, January 2007, 180 p. Xiang, Z., Swobada, G., Cen, Z., 2002. Identification of damage parameters for jointed rock. Computers and Structures 80, 1429–1440. Xiang, Z., Swoboda, G., Cen, Z., 2003. Parameter identification and its application in tunneling. In: Beer, G. (Ed.), Numerical Simulation in Tunneling. SpringerVerlag, Vienna, Austria, pp. 161–200. Yang, C.X., Wu, Y.H., Hon, T., Chen, D.M. 2006. Technique for determination of boundary stress conditions in deep tunneling. In: Proceedings of the ISRM International Symposium 2006: Rock Mechanics in Underground Construction, Singapore, p. 148.