Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
Volume 2013, Article ID 258721, 11 pages
Research Article

Application of Taguchi Method and Genetic Algorithm for Calibration of Soil Constitutive Models

1Civil and Environmental Engineering Faculty, Tarbiat Modares University, Tehran 1411713116, Iran
2Civil Engineering Faculty, Islamic Azad University of Tehran, Central Tehran Branch, Tehran, Iran
3Department of Civil Engineering, University of Tabriz, Tabriz 51666-14766, Iran

Received 30 April 2013; Revised 27 August 2013; Accepted 18 September 2013

Academic Editor: Pengcheng Fu

Copyright © 2013 M. Yazdani et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A special inverse analysis method is established in order to calibrate soil constitutive models. Taguchi method as a systematic sensitivity analysis is conducted to determine the real values of mechanical parameters. This technique was applied on the hardening soil (as an elastoplastic constitutive model) which is calibrated using the results from pressuremeter test performed on “Le Rheu” clayey sand. Meanwhile, a genetic algorithm (GA) as a well-known optimization technique is used to fit the computed numerical results and observed data of the soil model. This study indicates that the Taguchi method can reasonably calibrate the soil parameters with minimum number of numerical analyses in comparison with GA which needs plenty of analyses. In addition, the contribution of each parameter on mechanical behavior of soil during the test can be determined through the Taguchi method.

1. Introduction

One of the most important aspects of geotechnical problems is to adopt a suitable constitutive model for each material. Then, one or more appropriate experimental and/or field tests should be conducted to find the mechanical parameters of each constitutive model. When a set of parameters used in a model is selected so that it creates the most precise coincidence with the soil behavior, then the constitutive model is said to be calibrated. Generally, there are different methods ranging from simple to advanced for calibration of soil constitutive models. Simple conventional calibration techniques typically use stress and strain levels at certain states in which a material undergoes during specific types of laboratory tests. Sometimes, this method of calibration fails to capture the overall behavior of a material, that is, behavior at every point in stress-strain path [1]. For example, in a direct shear test, sets of normal and shear stresses in failure condition are used to find the peak values of internal friction angle and cohesion. However, often the other features of soil behavior such as the variation of shear stress versus shear displacement are not considered. Therefore, there is a vital need to fill this gap and find a much more comprehensive way to calibrate the constitutive models for soils. The best-proposed method to satisfy this requirement is the inverse analysis technique, which is based on mathematical solutions to find the best match between stress-strain curves.

Many researchers adopted inverse analysis method with different modifications. Cekerevac et al. [2] proposed an inverse calibration approach in which quasi-Newton and stochastic methods were used as optimization tools. They employed this method to calibrate Hujeux constitutive model for the results of isotropically consolidated drained triaxial compression tests. Quasi-Newton and stochastic methods were used to search for local and global minimums, respectively [2]. Calvelloet applied the inverse analysis techniques to calibrate hardening soil (HS) constitutive model for Chicago glacial clays. They used the results of triaxial compression tests along with the displacement profile recorded from inclinometer readings in a supported excavation in glacial clays [3]. In these researches, classical optimization tools were used. These methods are based on the derivatives of the objective function. However, such optimization techniques may lead to computational difficulties during the calculation of error function derivatives [46].

In this research, a new systematic search technique is proposed on the basis of genetic algorithm (GA) [712] and Taguchi method [1316]. GA as a well-known metaheuristic algorithm can be utilized to calibrate any soil constitutive model by means of the results obtained from any laboratory test and/or in situ experiment [2]. GA needs only an objective function rather than its derivatives. In this way, the shortcomings of classical methods will be eliminated as a result. However, in order to decrease the computational time, sensitivity analyses are required to select only the dominant parameters when the input parameters affecting the mechanical behavior are numerous. In this study, sensitivity analyses are carried out systematically using the well-known Taguchi method. This method which is conventionally used for the design of laboratory experiments can be treated as a modern technique in geotechnical application.

Genichi Taguchi, who first introduced this method during the late 1940s, utilized the conventional statistical tools in a simplified form by identifying a set of stringent guidelines for experiment layout and the analysis of results [13]. He made an applicable method for design and analysis of factorial experiments which is mainly used in quality engineering. This method, well known for its industrial applications to identify sensitive parameters for a given target, has fewer applications in geotechnics, particularly on material property identification [13].

In this paper, the results of pressuremeter tests [17, 18] which are performed on clayey sand in “Le Rheu” site located in France have been adopted for calibration of soil constitutive model [19]. The proposed method for inverse calibration is expressed using a special example which entails the curve obtained by pressuremeter test in a particular depth [20]. The constitutive law of this soil is assumed to be HS model due to the behavior that is exhibited during laboratory results. Thereafter, the inverse calibration is repeated with the reduced number of input parameters, obtained from the Taguchi method.

2. Specifications of the Soil in “Le Rheu” Site

The site is located in the west part of France, in a region called “Le Rheu.” The soil of this site contains reddish sand for tens of meters. Several in situ and laboratory tests have been performed on this soil to identify its mechanical and engineering characteristics. The main reason for selection of this site in current research is the uniformity of the soil type in different depths and the existence of water table at very low levels. These conditions reduce the complexity of modeling process and let all efforts be concentrated on the mathematical solution for inverse calibration.

The results of pressuremeter tests are available at three points of B4, P1, and P2 (Figure 1) in various depths of 2 m, 3 m, 4 m, and 5 m [19]. However, in this study, only the curve related to point B4 at the depth of 2 m was selected. Figure 1 illustrates the results of tests at point B4 in the form of curve.

Figure 1: Pressuremeter curves at a depth of 2 m after being modified by lift-off method.

3. Hardening Soil Model

The hardening soil model is an advanced model for simulating the behavior of both soft and stiff soils. When subjected to primary deviatoric loading, the soil shows a decreasing stiffness and simultaneously irreversible plastic strains developing. In the special case of a drained triaxial test, the observed relationship between the axial strain and the deviatoric stress can be well approximated by a hyperbola function, as (Figure 2):

Figure 2: Hyperbolic stress-strain relation in primary loading for a standard drained triaxial test.

In contrast to an elastic-perfectly plastic model, the yield surface of a hardening plasticity model is not fixed in principal stress space, but it can expand due to plastic straining. Distinction can be made between two main types of hardening, namely, shear hardening and compression hardening. Shear hardening is used to model irreversible strains due to primary deviatoric loading. Compression hardening is used to model irreversible plastic strains due to primary compression in oedometer loading and isotropic loading. Both types of hardening are contained in the present model. Table 1 presents input parameters of the HS model.

Table 1: HS input parameters [3].

4. Numerical Modeling of Menard Pressuremeter Test

The first step of numerical modeling is generating the geometry. A mass of soil should be considered in which a borehole is dug to model the pressuremeter test. Then the pressure is induced to the boundary of the soil element adjacent to the middle cell of the probe. Because of the symmetric geometry and loading, only a half of the geometry is modeled (Figure 3). Three regions are identified in Figure 3, as follows:

Figure 3: Geometry of Menard pressuremeter test model.

Region 1 shows soil mass around the borehole in which the stress-field induced by loading can be assumed negligible.

Region 2 is the area exposed to direct impact of induced pressure, so it needs a finer mesh. This area consists of a 20 cm × 20 cm square, where the height implies the height of pressuremeter probe middle cell. Loading occurs on the inner boundary of this square (left side of the square in Figure 3).

Region 3 stands for the borehole which will be eliminated at the first phase of analysis. Illustrated dimensions in Figure 3 are as follows:: distance between probe center and ground surface (i.e., depth of experiment) = 2 m;: probe height = 1.5 m;: probe center distance to the bottom of the borehole (i.e., half height of the probe) = 75 cm;: depth of the probe = 1.5 m;: 50 times of the borehole diameter (50 × 0.06 m = 3 m).

The boundary conditions and loading position are defined in Figure 4. The pressure, which is made by expansion of the middle cell of the probe, will be induced in analysis phase.

Figure 4: Boundary conditions and loading position.

Mesh is generated in the next step as shown in Figure 5. Since the displacements and stresses produced in region 2 are very important, a finer mesh is considered for this part. According to high value of height-width ratio of region 3, a refined mesh is needed in this region. To increase the precision of calculations, 15-node triangular elements have been used.

Figure 5: Mesh generation for pressuremeter model.

For the current pressuremeter test modeling, analysis phases have been defined as follows:

Phase 1. Borehole is excavated and the stresses due to the excavation are calculated. Calculations in this phase are in the plastic zone of the soil.

Phase 2. Pressuremeter apparatus is planted in the desired depth of the borehole (2 meters in this case) and the experiment starts by inducing a 100 kPa pressure. In this phase, displacements of the previous phase, due to the borehole excavation, are set to zero.

Phases 3 to 46. In subsequent phases, pressure increases gradually. In this experiment a 100 kPa incensement is considered for each step. Therefore, in phase 3, we have  kPa and in phase 4,  kPa, and so forth until phase 46 which it is  kPa. It should be mentioned that from phase 2 on, calculations are updated according to the mesh type and produced large displacements and they may not necessarily continue to phase 46. The final step depends on the time of failure.

5. Inverse Analysis for Calibration of Soil Constitutive Models

In inverse analysis, a given model is calibrated by iteratively changing input values until the simulated output values match the observed data [3]. The basic form of inverse analysis technique can be categorized as a trial and error approach (Figure 6). When the number of input parameters is too large, this method may be inefficient or impractical. Therefore, to avoid this troublesome effort, providing a systematic approach seems to be necessary. In the following section, an optimization tool is introduced in order to systematically minimize the difference between numerical and experimental results.

Figure 6: General inverse analysis diagram for calibration of soil constitutive models.
5.1. Systematic Inverse Analysis Method

The given constitutive model is calibrated by a repetitive procedure in systematic inverse analysis. In this cycle, input parameters of the constitutive model are changed until the results of numerical simulation match the experimental responses. In this research, the results of Menard pressuremeter tests have been considered as the soil response used for calibration of HS model. A set of input parameters for soil constitutive model which leads to the coincidence of in situ pressuremeter curve and model pressuremeter simulation curve is desired. There is an extreme need for a quantity, which shows the degree of coincidence between the two mentioned curves in order to solve the problem. This quantity which is error function is generally defined as “area between the two curves,” as This concept is illustrated in Figure 7.

Figure 7: Concept of error function.

In this paper, an error function with the following form is used, as where represents the summation of its subsequent term ( discrete values) and is the number of used experimentally obtained data in the process.

Therefore, the calibration is changed into a familiar optimization problem in which finding a feasible set of soil’s model parameters leads to the least value for error function. Soil constitutive model parameters are those 6 parameters previously introduced in Table 1. Since there is a need to change the level of each parameter without any limitations, the parameter is eliminated from the inverse analysis procedure. As there is not the possibility for to be changed freely, this parameter should be removed from the cycle and the default value for this parameter will be accepted (). Thus, the number of input parameters reduces to 5.

Now, this idealized problem is ready to be solved. The optimization tool used in this research is GA. There are many computer programs written for GA, but none is able to communicate with PLAXIS. To solve this problem, instead of using available programs for GA, a code is written for GA by Visual Basic (VB), which has the ability to interface with the PLAXIS, a useful finite element program which can perform the analysis according to predefined stages. Therefore, this code can change the value of each parameter in that optimization process and obtain the objective function. Figure 8 presents the algorithm with more details.

Figure 8: The algorithm of the written code for systematic inverse analysis.

The best set of parameters obtained by this method is gained after 496 cycles as shown in Table 2. Figure 9 illustrates a very good coincidence between the in situ and simulation curves. Inverse analysis algorithms allow simultaneous calibration of multiple input parameters [3]. On the other hand, the required time for inverse analysis intensively increases by increasing the number of parameters. However, the computational time can be reduced to a large extent by removing some unimportant parameters. A sensitivity analysis attains the degree of importance of each parameter [21, 22]. In this paper, “Taguchi method” is used to fulfill this aim.

Table 2: The best set of parameters obtained by systematic inverse analysis via GA optimization tool.
Figure 9: In situ pressuremeter curve in comparison with the best simulation curve, pressuremeter simulation curve obtained by inverse analysis (performed after sensitivity analysis), and pressuremeter simulation curve obtained by substituting field attained-parameters (direct method).
5.2. Sensitivity Analysis by Taguchi Method

Taguchi method is conventionally an approach for sensitivity analysis method, by changing a selected factor in different levels, while the other factors are kept constant. Then, the same process repeats exactly for each of the remaining factors. In Taguchi method, all factors are changed simultaneously according to predefined tables called “orthogonal arrays.” Choosing the appropriate orthogonal array for a given problem is called “experiment design.” The first step to perform a systematic sensitivity analysis is to define experiment design. In order to generate design experiments (i.e., finding the suitable orthogonal array), “degrees of freedom” is needed, which is obtained as follows: In this study, for each of the 5 factors, 4 levels are considered. Therefore, the degree of freedom for each factor equals 3 (, = number of levels for factor ). Interactions’ degree of freedom will be zero since no interaction is considered. By substituting the mentioned values into (4), will be 15.

The smallest orthogonal array with the degree of freedom greater than (or equal to) the experiment degree of freedom should be found in this step. Degree of freedom for L16 array is 15: , so L16 array can be obtained (Table 3). But L16 contains only 2-level factors, while an orthogonal array with 4-level factors is needed. Therefore, using the rule of converting 2-level columns into 4-level columns, M16 orthogonal array is achieved (Table 4). Variation interval of each factor is divided into 4 equal divisions as mentioned before, thus factor levels will be as Table 5. Factors can be assigned to columns of orthogonal array M16, now. Here, as interaction between factors has not been taken into account, the factors will arbitrarily be assigned to any desired column of M16.

Table 3: L16 orthogonal array.
Table 4: Modified L16 orthogonal array (M16).
Table 5: Considered levels for each factor.

Final plan of experiments is shown in Table 6. In this table, each row stands for an experiment, so the pressuremeter finite element model should be run 16 times, according to the conditions of the orthogonal array M16. The results obtained after running these experiments are shown in the last column of Table 6.

Table 6: Final plan of experiments for this project.

Obtained data of Table 6 are analyzed according to Taguchi ANOVA table (Analysis of variance). Results are shown in Table 7. The last column of Table 7 shows the contribution percent of each parameter. Contribution percent shows the sensitivity degree of numerical model response with respect to each parameter variations. As it can be seen in this table, parameter has the most, and parameter has the least degree of importance (sensitivity degree).

Table 7: Results of ANOVA table.

The parameter with the degree of importance less than 10% of the most significant factor will be assigned a constant value and removed from the inverse analysis process. As a result, parameter has a very small degree of importance (4.8%). This value is less than 10% of the importance degree of the most significant parameter (here with contribution percent of 53.2%). Therefore, a constant value is assigned to (here, ). Now, inverse analysis can be performed with the 4 remaining parameters. The result of this analysis is shown in Table 8. Figure 9 illustrates the simulated pressuremeter curve obtained from numerical analysis based on Table 8 parameters in comparison with the in situ pressuremeter curve.

Table 8: The best set of soil constitutive model parameters obtained by systematic inverse analysis after performing sensitivity analysis and removing the unimportant parameter.
5.3. Comparing to Results of Direct Calibration

The HS constitutive model for “Le Rheu” soil has been calibrated directly [5]. In situ and laboratory tests were utilized to assess the value of HS soil model parameters as shown in Table 9. For example, vane shear test is used to estimate the values of and . After substituting obtained parameters into the simulated pressuremeter test model and running the numerical model, stress-volumetric strain curve is attained. This curve is shown in company with in situ pressuremeter curve in Figure 9. As it can be seen, the two curves have a similar trend and they are nearly parallel but there is no close coincidence.

Table 9: The set of parameters obtained by direct method (experimental method).

6. Discussion

The main purpose of the paper is to introduce a systematic approach to derive mechanical parameters of a typical soil constitutive model based on an available test data. For adopted example of “Le Rheu” clayey sand, three different test data related to three points of B4, P1, and P2 were available at depth of 2 m. The proposed approach can be applied to the curve of each point independently, and then the corresponding mechanical parameters may be averaged to represent the mean values of the soil mechanical parameters of “Le Rheu” clayey sand at depth of 2 m. According to Figure 1, the parameters obtained from the points P1 and P2 should be similar but they might be different from the point of B4. Therefore, there was no specific reason for selection of point B4, since the target was a presentation of the method.

Taguchi method was originally proposed to design experiments. However, in this paper it was adopted to derive the mechanical parameters of a soil through systematic inverse analyses. On the other hand, GA is an optimization technique which was utilized here to obtain the optimum parameters fitting to an available soil test data. Though, the above two methods are different tools in engineering and scientific practice, in this paper they were utilized for a single specific application, that is, the calibration of a soil constitutive model. Accordingly, the comparison achieved in the paper between Taguchi and GA methods is only attributed to the precision of the results and the number of analyses needed in each method. In addition, giving the relative significance of each mechanical parameter in soil constitutive model is another ability of the method based on Taguchi approach.

The results of obtained parameters (Table 8) and its corresponding Figure 9 have been obtained only through the Taguchi method without any need to GA. In the first cycle of Taguchi method, 5 soil parameters were selected (5 factors). However, since one of those parameters observed to have little significance respect to the others, it was decided to assign it a constant value and run the second cycle of Taguchi with only 4 parameters. For this study, in each cycle of Taguchi method 16 analyses have been carried out based on orthogonal arrays of L16 (or M16). Table 4 presents the level of every factor (parameter) for each of 16 analyses. The values of factors in each of the 16 tests (analyses) were presented in Table 6. Thus, there is no need for GA in this approach.

Taguchi method is a systematic approach for designing experiments which investigates how different parameters affect the mean and variance of a process performance characteristic. However, it is very important to determine the most important parameters (factors) governing the process since the total number of parameters involving the process might be high. In addition, the variation range of each parameter should be introduced as much as limited in order to define minimum number of levels. These considerations may need some experiences and, without such information, the method may not be effective and useful. Having a large number of parameters (factors) with a wide range of variation for each parameter tends to select the orthogonal arrays with numerous tests. This will be time consuming and expensive from computational costs point of view.

Regarding the ability of the method to be applied on other tests or constitutive behaviors, it can be useful to say that we have already utilized the method in order to extract the Mohr-Coulomb perfect plastic parameters of soil from the results of pile load tests [23]. In another research, the HS constitutive parameters of rock masses in site of “Siah-bisheh” were estimated from the monitoring results of powerhouse cavern [24].

7. Conclusions

In this research, a systematic inverse analysis approach is introduced for calibration of soil constitutive models. The capability of this method has been shown in the case of calibrating HS constitutive model for “Le Rheu” soil in pressuremeter stress path. The benefits of using this method are being able to be used for many laboratory or field tests, and also constitutive models, giving the whole parameters simultaneously, automatic procedure of calibration with least interpretation, and considering overall soil behavior (i.e., behavior at every point in stress-strain path).

The Taguchi method is a useful tool for parametric analysis which can be beneficial in geotechnical engineering due to its relatively high precision and low time consumption. Furthermore, the significance of the parameters can be evaluated quantitatively using the Taguchi method. In the current research, it was exhibited that the parameters of soil cohesion and internal friction angle have the most influence on the hardening soil elastoplastic constitutive model and the dilatancy angle has the least influence. This conclusion is probably valid only for clayey sand located in Le Rheu site. For granular soils with large size grains such as gravels in which the dilatancy angle is large, it is possibly expected to observe more contribution of dilatancy.

As illustrated in Tables 2, 8, and 9, based on the error function values and the calculation time, it is obvious that the Taguchi method is faster than both direct method and the single GA and more precise than the direct method. The results obtained from the Taguchi method are close to the GA, but with less computational time. As shown in Table 10, an error function of 21.8 was achieved with 496 analyses of the GA method. The direct method gave an error function of 75.3 with a very low precision. However, an error function of 28 was concluded with mere 16 analyses, using the Taguchi method. Hence, it is obvious that the Taguchi method is a cheap and fast method to gain acceptable results.

Table 10: Comparison of used methods.


  1. S. Pal, G. W. Wathugala, and S. Kundu, “Calibration of a constitutive model using genetic algorithms,” Computers and Geotechnics, vol. 19, no. 4, pp. 325–348, 1996. View at Google Scholar
  2. C. Cekerevac, S. Girardin, G. Klubertanz, and L. Laloui, “Calibration of an elasto-plastic constitutive model by a constrained optimisation procedure,” Computers and Geotechnics, vol. 33, no. 8, pp. 432–443, 2006. View at Publisher · View at Google Scholar · View at Scopus
  3. M. Calvello and R. J. Finno, “Selecting parameters to optimize in model calibration by inverse analysis,” Computers and Geotechnics, vol. 31, no. 5, pp. 411–425, 2004. View at Publisher · View at Google Scholar · View at Scopus
  4. A. Kaveh and S. Talatahari, “Particle swarm optimizer, ant colony strategy and harmony search scheme hybridized for optimization of truss structures,” Computers and Structures, vol. 87, no. 5-6, pp. 267–283, 2009. View at Publisher · View at Google Scholar · View at Scopus
  5. A. Kaveh and S. Talatahari, “A novel heuristic optimization method: charged system search,” Acta Mechanica, vol. 213, no. 3-4, pp. 267–289, 2010. View at Google Scholar
  6. A. Kaveh and S. Talatahari, “Optimal design of skeletal structures via the charged system search algorithm,” Structural and Multidisciplinary Optimization, vol. 41, no. 6, pp. 893–911, 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. D. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Boston, Mass, USA, 1989.
  8. M. D. Vose, The Simple Genetic Algorithm: Foundations and Theory, MIT, Cambridge, Mass, USA, 1999. View at MathSciNet
  9. M. Gen and R. Cheng, Genetic Algorithms and Engineering Design, John Wiley & Sons, New York, NY, USA, 1997.
  10. H. Demuth and M. Beale, Genetic Algorithm and Direct Search Toolbox User’s Guide for Use with MATLAB, version 1, MathWorks, 1st edition, 2004.
  11. W. M. Spear, Adapting Crossover in a Genetic Algorithm, Naval Research Laboratory, Washington, DC, USA, 2003.
  12. S. Levasseur, Y. Malécot, M. Boulon, and E. Flavigny, “Soil parameter identification using a genetic algorithm,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 32, no. 2, pp. 189–213, 2008. View at Publisher · View at Google Scholar · View at Scopus
  13. R. K. Roy, A Primer on the Taguchi Method, Van Nostrand Rein hold, New York, NY, USA, 1990.
  14. P. J. Ross, Taguchi Techniques for Quality Engineering: Loss Function, Orthogonal Experiments, Parameter and Tolerance Design, McGraw-Hill, New York, NY, USA edition, 1996.
  15. K. Balaraman, S. Mukherjee, A. Chawla, and R. Malhotra, Inverse Finite Element Characterization of Soft Tissues Using Impact Experiments and Taguchi Method, SAE International, 2005.
  16. R. K. . Roy, Design of Experiments Using the Taguchi Approach: 16 Steps to Product and Process Improvement, Wiley, New York, NY, USA, 2001.
  17. B. G. Clarke, Pressuremeters in Geotechnical Design, Blackie Academic and Professional, an Imprint of Chapman and Hall, 1st edition, 1995.
  18. F. Baguelin, J. F. Jezequel, and D. H. Shields, The Pressuremeter and Foundation Engineering, Trans Tech, Clausthal, Germany, 1978.
  19. Simulation of the Menard Pressuremeter Test Using the F.E.M Method, Ecole Centrale, Paris, France, 1997.
  20. “Standard test method for pressuremeter testing in soils,” ASTM Standard D4719-87, 1988.
  21. Z. Yang and A. Elgamal, “Application of unconstrained optimization and sensitivity analysis to calibration of a soil constitutive model,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 27, no. 15, pp. 1277–1297, 2003. View at Publisher · View at Google Scholar · View at Scopus
  22. N. Dien, F. Jean-Marie, and M. Arezou, “Identification of parameters and validation of Hujeux model for a coarse soil,” Laboratory of Soil, Structure and Material Mechanics CNRS UMR 8579, Ecole Centrale, Paris, France.
  23. M. Yazdani and P. Ghahremani, “Determination of bearing capacity for piles using systematic numerical analyses of pile load testing based on Taguchi method,” Civil Engineering Sharif Journal, vol. 28-2, no. 4, pp. 39–51, 2013. View at Google Scholar
  24. A. Majnouni, Estimation of rock mass parameters in Siah-Bisheh powerhouse cavern using Taguchi method [M.S. thesis], Tarbiat Modares University, Tehran, Iran, 2012.