Elastic–plastic model identification for rock surrounding an underground excavation based on immunized genetic algorithm

To compute the stability of underground engineering, a constitutive model of surrounding rock must be identified. Many constitutive models for rock mass have been proposed. In this model identification study, a generalized constitutive law for an elastic–plastic constitutive model is applied. Using the generalized constitutive law, the problem of model identification is transformed to a problem of parameter identification, which is a typical and complicated optimization. To improve the efficiency of the traditional optimization method, an immunized genetic algorithm that is proposed by the author is applied in this study. In this new algorithm, the principle of artificial immune algorithm is combined with the genetic algorithm. Therefore, the entire computation efficiency of model identification will be improved. Using this new model identification method, a numerical example and an engineering example are used to verify the computing ability of the algorithm. The results show that this new model identification algorithm can significantly improve the computation efficiency and the computation effect.

very slowly. Different from parameter identification, in model identification, the constitutive model of surrounding rock is unknown. Based on the measurement information, such as displacement, the structure and parameters of the constitutive model can be identified. Because it is very complex to select the structure and parameters of the constitutive model at the same time, generally, the structure of the constitutive model is determined by the prior knowledge, such as engineering experience, and then only the parameters of the constitutive model are identified based on measurement information. Thus, the model identification can be simplified as similar as the parameter identification. However, the model identification is more complex than the parameter identification. In model identification, the model parameters and the mechanical parameters can be identified at the same time. And the structure of the constitutive model determined by the prior knowledge is generally more complex than the assumed simple model used in parameter identification. Moreover, different from the mechanical parameters, which have the clear physical and mechanical meaning, and can be determined by the tests easily, the model parameters, which only describe the constitutive model, can not be determined by the tests. In 1987, Gioda and Sakurai proposed that model identification based on displacement should be the main development for back analysis (Gioda and Sakurai 1987). In 1997, Sakurai demonstrated that the identification of a constitutive model is critical (Sakurai 1997). Therefore, researchers have investigated the identification of a constitutive model. Based on the displacement measurements of an underground roadway, Liu (2011) identified the visco-elastic constitutive model of rock mass based on the traditional nonlinear optimum technique. Wang et al. (2007) identified the geo-material constitutive model based on their constitutive model database and several identification algorithms. Yang and Wang (2009) presented a numerical model to identify the unknown equivalent constitutive model in the elastic layered rock mass of an underground opening by the Gauss-Newton technique.
Because model identification is a very complicated optimization problem (Gao and Liu 2009), whose objective function can be shown as in Fig. 1, traditional optimization techniques have some shortcomings. The relationship between objective function and the optimization method can be described as in Fig. 2 (Gao and Liu 2009). From Fig. 2, the global optimization is one suitable method to solve model identification problem. Therefore, some important contributions to global optimization algorithms have been achieved. Su et al. (2008) identified the structure and parameters of the rheological constitutive model of the surrounding rocks of the Jinping tunnels in China using a differential evolution algorithm. Feng et al. (2006) identified a visco-elastic model of surrounding rocks in the Goupitan hydroelectric power station in China using genetic programming Model Objective function Fig. 1 Optimization objective function of model identification for underground engineering and a modified particle swarm optimization algorithm. Meier et al. (2007) have performed model identification of surrounding rocks for tunneling by particle swarm optimization. Sha et al. (2011) presented three methods to identify the constitutive model of rocks based on a genetic algorithm, a back propagation neural network and genetic programming and compared all methods. Surajit and Wathugala (1996) calibrated a constitutive model using genetic algorithms. Gao (2007) proposed a method to identify the rheological constitutive model for rock mass of an underground roadway based on the fast-convergent genetic algorithm. Moreover, Gao et al. (2004) presented an identification method for the surrounding rock of an underground power house based on a new intelligent bionics algorithm. Because these global optimization algorithms comprise random search algorithms, their computational efficiency is very low. This low efficiency is problematic when the model identification is extremely complicated. In this study, a new method to identify the constitutive model based on an immunized genetic algorithm is proposed. The elastic-plastic constitutive model is the main model of geomaterials which has been comprehensively evaluated (Nakai 2012;Zheng et al. 2002). In the last years, many studies (Cui et al. 2015;Lee and Pietruszczak 2008;Spiezia et al. 2016;Zhang et al. 2012) have been carried out investigating the elasto-plastic behavior of rock surrounding underground excavations. And those researches have proved that the elastic-plastic constitutive model can describe the mechanical behaviour of rock surrounding underground excavation very well. Therefore, the elastic-plastic constitutive model is analyzed in this study.

Generalized constitutive law of elastic-plastic model for rock materials
The generalized constitutive law of surrounding rock in strain space is where D ijkl is the tensor of the elastic stiffness coefficient, g(σ ij , k) is the plastic potential function, f(σ ij , σ ij p , k) is the yield function of the stress space,  Fig. 2 Relationship between objective function and optimization method function of the strain space, σ ij p and ɛ ij p the plastic stress and strain, respectively, A is the hardening or softening parameter, whose form is as follows (Zheng et al. 2002;Zheng and Kong 2010): k is a plastic internal variable, h is a parameter, whose form is as follows (Zheng and Kong 2010): where w p is the plastic power, θ p is the plastic expansion, ε p is the equivalent plastic strain.
The switch function �·� describes the loading and unloading criterion, whose form is as follows (Zheng and Kong 2010): The relationship of the two yield functions is as follows: Equation (1) can be summarized by a matrix as where D p is a plastic matrix, D ep is the elastic-plastic matrix [D ep Therefore, the elastic-plastic constitutive model law can be described by a plastic matrix. According to the mechanics derivation (Koichi 2009), the plastic matrix [D p ] can be described as where Φ is the loading function, G is the plastic potential.
In this study, the surrounding rock is supposed as perfect elastic-plastic. Therefore, its elastic-plastic constitutive model can be simply described as follows (Nakai 2012): F can be described as (Zheng et al. 2002) where I 1 is the first invariant of stress, J 2 is the second invariant, J 3 is the third invariant of the stress deviator.
For the non-associative flow rule, the plastic potential G is different from the yield function F. And it is a very hard work to construct the function of G (Nakai 2012). Moreover, there is no generalized form for the function of G (Zheng et al. 2002). However, the elastic-plastic constitutive model with the associative flow rule can describe the main mechanical behaviours of geo-meterials well (Zheng and Kong 2010). Thus, in this study, simply, G = F. Therefore, only two terms are to be determined in this elastic-plastic constitutive law: the stress yield function F and the elastic matrix [D e ]. The problem of identification for the elastic-plastic constitutive model can be transformed to the problem of identification for the stress yield function F and the elastic parameters E and μ.
The generalized form of the stress yield function is (Zheng et al. 2002) where is a function to describe the shape of yield function (Liu and Zheng 1997;Zheng and Kong 2010), θ is the Lode angle, α, β and K are unknown parameters.
Because the Eq. (10) is the generalized form of the yield functions, with the different parameters, mainly α and K, the different yield functions can be obtained (Zheng et al. 2002). These yield functions include almost all widely used yield functions, such as Mohr-Coulomb (M-C) yield criterion, Mises yield criterion, two improved Mises yield criterions (exterior angle circle and interior angle circle), Drucker-Prager (D-P) yield criterion, Tresca yield criterion, three Zienkiewicz-Pande (Z-P) yield criterions (elliptic curve, hyperbolic curve and parabolic curve) and twin-shear yield criterion, etc. The detailed specific relationship between different parameters and different yield functions can be found in reference (Zheng et al. 2002). The curves of the Eq. (10) for some generally used yield functions are shown in Fig. 3.
Using this generalized constitutive law for an elastic-plastic model, the problem of model identification can be transformed to a parameter identification.
To simply analyze without a loss of generalization, we assume that the parameter β = 0; thus, only the two parameters α and K need to be identified. In addition to the elastic parameters E and μ, four model parameters need to be identified: E, μ, α and K.

Model identification by immunized genetic algorithm
To improve the main operations of genetic algorithm, such as creation of initial population, reproduction operation, crossover operation, mutation operation and selection operation, one new immunized genetic algorithm has been proposed by author, whose flow chart is as shown in Fig. 4. The details of the immunized genetic algorithm can be found in reference (Gao 2012).
Based on new immunized genetic algorithm, the model identification method can be proposed, whose flow chart is as shown in Fig. 5.
The detailed procedure of the model identification is described in the following.

Parameter initialization
Before the algorithm is performed, all parameters must be initialized. These parameters include the number of individuals in the initial population, the maximum number of

Individual expression
In this study, an individual is expressed as X = (E, μ, α, K). An individual corresponds to an elastic-plastic constitutive model.

Creation of the initial population
Initial population is created by SRCM (Gao and Yin 2011).

Fitness function
In this study, the objective function is defined as where e(i) is the error between the displacement computations and the displacement measurements and Num is the number of monitor points. The fitness function is as follows:

Optimization process
The optimization operation is as same as in immunized genetic algorithm (Gao 2012).

Termination condition
In this study, the termination criterion specifies that the difference of the maximum fitness value and the average fitness value is less than 10e-5. To avoid infinite iteration, a maximum number of evolutionary generations is also specified.

Numerical example
A numerical testing example is employed to test the proposed algorithm. The example is an underground tunnel with a radius of 3 m. The details of this example can be found in reference (Gao 2016).
For computation, the mechanical property of the surrounding rock is assumed as follows: The surrounding rock is an elastic-plastic material. The associated parameters are Young's Modulus E = 2100 MPa, Poisson's ratio μ = 0.2, cohesion C = 1.1 MPa and friction angle ϕ = 30°. Because the yield function is the Drucker-Prager theory, the definitions of α and K (Zheng et al. 2002) are as follows, According to the definition of α and K, we can obtain the values of α and K from the values of C and ϕ, which are as follows: α = 0.165872, K = 0.875726 MPa.
Because the ranges of the parameters only affect the computing efficiency and have minimal influence on the computing results, the ranges for the four parameters are as follows: For comparison study, the traditional genetic algorithm (GA), fast genetic algorithm (FGA) (Gao 2007) and the immunized genetic algorithm (IGA) are all applied for this example.
Based on testing and experience, the parameters for three algorithms are summarized in Table 1.
To compare the computation effect, the results of GA, FGA and IGA are summarized in Table 2.  Table 2, the results of the IGA are superior to the results of the GA and FGA; therefore, the computation effect in this study is reasonable.

As shown in
To simply evaluate the computational efficiency of the GA, FGA and IGA, the number of objective function evaluations during the search for the optimum, which is denoted by NOF and represents the computation time required by the optimization algorithm, is applied. The NOFs of the GA, FGA and IGA methods are given in Fig. 6.
As shown in Fig. 6, the NOF of IGA is much less than those of other algorithms, and then the computational speed of IGA in this study is very faster than those of other algorithms in the literature. In other words, the computational efficiency of the IGA in this study is the best and is superior to those of other algorithms.
The results of these studies conclude that the IGA method can be used to obtain a suitable model with higher accuracy and less effort.

Engineering example
The Huainan coal mine is located north of Huainan city, of Anhui Province in China. As an old mining, the entire mine has been integrated into the stage of deep mining. The different geological environments range from shallow rock roadway to the deep rock roadway; as a result, the mechanical characteristics are very complicated. To analyze the failure mechanism of the surrounding rock for deep rock roadway, a constitutive model of the surrounding rock based on the displacement measurement results must be identified. In this study, two main types of surrounding rock are identified: types III and II (Liu et al. 2010). The surrounding rock of the −780 C S 13 floor haulage roadway in the Xieyi mine is characterized as the moderate poor type, which is associated with type III. The depth of this roadway is 800 m. Its main lithology consists of fine sandstone and limestone, and its integrity is satisfactory. The roadway cross-section is shown in Fig. 7. Its width is 3.8 m, the height of the side wall is 1.3 m, the height of the crown is 1.9 m and the depth of the floor arch is 0.5 m.
To analyze the stability of the roadway, some site monitoring studies have been performed, including surface displacement measurements and deep displacement measurements.
For the surface displacement measurements, the convergence between both sides and the convergence between the roof and the floor were completely measured. The layout of the monitoring points is illustrated in Fig. 7. The layout of the monitoring points for the deep displacement measurements is also illustrated in Fig. 7. The depth of each monitoring hole was 15 m. Six monitoring points were installed along each monitoring hole, with distances of 1, 3, 5, 7, 10 and 15 m, which are designated point 1, point 2, point 3, point 4, point 5 and point 6. Assuming the point at a depth of 15 m is stationary, the relative displacement of the remaining points can be obtained.
The surface convergence displacement measurement and deep multi-point displacement are shown in Figs. 8 and 9.
According to the field monitoring of the hydraulic fracturing technique, the vertical initial stress is 13.2 MPa and the horizontal initial stress is 19.5 MPa. The FEM model is shown in Fig. 10. The parameters of model identification method are as follows: The number of individuals is 150, and the maximum number of evolutionary generations is 500.
The identified model parameters are as follows: In previous studies (Liu et al. 2010), based on some laboratory experiments and engineering experience, the suggestion values for mechanical parameters of typical surrounding rock were given to conduct the support design. Because the suggestion values can consider the complicated geological environment for surrounding rock, they are very suitable to be used to test the identified parameters in this study. The suggestion values of mechanical parameters for surrounding rock of type III are as follows (Liu et al. 2010). E is between 6 and 15 GPa, μ is almost 0.3, C is between 1.0 and 1.5 MPa and ϕ is between 30° and 45°.
Thus, the identified values for parameters E and μ are agree with the suggestion values well. However, the two model parameters α and K can not be obtained directly from the tests, and they can be determined from the parameters c and φ based on the yield function. To obtain the model parameters, the yield function of the surrounding rock must be constructed, but it is a very hard work. Because this study is only to analyze the stability of surrounding rock, the hard work to construct the yield function of the surrounding   Fig. 11. For comparison, the stable displacement of each monitoring point is utilized.
As shown in Fig. 11, the computation displacements correspond with the measured displacements. Therefore, the identified model for the surrounding rock of the −780 C S 13 floor haulage roadway in the Xieyi mine can adequately describe the real rock performance.
The surrounding rock of the −720 rock level crosscut in the Xieqiao mine is the moderate type, which corresponds to type II. The depth of this roadway is 720 m. Its main lithology consists of fine sandstone, fine siltstone and medium fine sandstone, and its integrity is satisfactory. The cross-section of the roadway shown in Fig. 12. Its width is 4.5 m, the height of the side wall is 1.5 m and the height of the crown is 2.25 m.
The layout of the monitoring points is illustrated in Fig. 12. The layout of the monitoring points for the deep displacement measurements is also illustrated in Fig. 12. The depth of each monitoring hole was 15 m. Six monitoring points were installed along each monitoring hole, with distances of 1, 3, 5, 7, 10 and 15 m, which are designated point 1, point 2, point 3, point 4, point 5 and point 6. Assuming that the point at a depth of 15 m is stationary, the relative displacement of the remaining points can be obtained.
The surface convergence displacement measurement and deep multi-point displacements are shown in Figs. 13 and 14.
According to the field monitoring results, the vertical initial stress is 20.1 MPa and the horizontal initial stress is 15.6 MPa. The FEM model is shown in Fig. 15.
Because the ranges of the parameters only affect the computing efficiency and have minimal influence on the computing results, the ranges for the four model parameters are as follows: 1 GPa ≤ E ≤ 30 GPa, 0.15 ≤ µ ≤ 0.35, 0.0 ≤ α ≤ 0.85, 0.0 < K ≤ 4.5 MPa. The parameters of model identification method are as follows: The number of individuals is 150, and the maximum number of evolutionary generations is 500.
The identified model parameters are as follows: As the same studies for surrounding rock of type III, the two mechanical parameters E and μ are verified by the suggestion values of mechanical parameters in previous E = 14.6 GPa, µ = 0.26, α = 0.1871, K = 0.8841 MPa. The suggestion values of mechanical parameters for surrounding rock of type II are as follows (Liu et al. 2010).
E is between 10 and 25 GPa, μ is almost 0.25, C is between 1.2 and 2.0 MPa and φ is between 40° and 55°.
Thus, the identified values for parameters E and μ are agree with the suggestion values well. The model parameters α and K are verified by the comparison of measured deep As shown in Fig. 16, the computation displacements coincide with the measured displacements. Therefore, the identified model for the surrounding rock can adequately describe the real rock performance.
In this study, the surrounding rock is assumed as one equivalent homogeneous material. Therefore, the computed displacements at both sides are equal. In fact, the measurement displacements at both sides are different. In other words, the inhomogeneous displacement as shown in Figs. 11 and 16 can not be described in this study.
The engineering applications prove that the new model identification method proposed in this study can be used to identify a suitable constitutive model for the surrounding rock of the underground engineering with reasonable efficiency using only surface displacement measurements.

Conclusion
The identification of a rock constitutive model is very important. From the theory analysis of the elastic-plastic constitutive law of rock materials, the generalized law of constitutive model is presented. Using this generalized law, the problem of identification of a constitutive model can be transformed to parameter identification. Therefore, the new immunized genetic algorithm proposed by the author is applied to an elastic-plastic model identification problem. Using a numerical example and an engineering example, this new method is verified. The results indicate that the proposed method can be used to identify a suitable constitutive model with high accuracy and efficiency.
Authors' contributions WG carried out the new algorithm studies and engineering applications, participated in the sequence alignment and drafted the manuscript. DC carried out the constitutive law studies. XW conceived of the study, and participated in its design and helped to draft the manuscript. All authors read and approved the final manuscript.

Acknowledgements
The financial supports from The Fundamental Research Funds for the Central Universities under Grant Nos. 2014B17814, 2014B07014, 2016B10214 and B15020060 are all gratefully acknowledged.