- Email: [email protected]
Prediction of unconfined compressive strength of soft grounds using computational intelligence techniques: A comparative study
Computers and Geotechnics 33 (2006) 196–208 www.elsevier.com/locate/compgeo
Prediction of unconfined compressive strength of soft grounds using computational intelligence techniques: A comparative study B.S. Narendra a, P.V. Sivapullaiah a
a,*
, S. Suresh b, S.N. Omkar
b
Department of Civil Engineering, Indian Institute of Science, Bangalore 560 012, India Department of Aerospace Engineering, Indian Institute of Science, Bangalore, India
b
Received 14 December 2004; received in revised form 2 February 2006; accepted 13 March 2006 Available online 18 May 2006
Abstract Cement stabilization is one of the commonly used techniques to improve the strength of soft ground/clays, generally found along coastal and low land areas. The strength development in cement stabilization technique depends on the soil properties, cement content, curing period and environmental conditions. For optimal and effective utilization of cement, there is a need to develop a mathematical model relating the gain in strength in terms of the variables responsible. The existing empirical model in the literature assumes linear variation of normalized strength with the logarithm of curing period and hence, different empirical models are required for different conditions of the same soil. Also, the accuracy of strength prediction is unsatisfactory. Due to unknown functional relationships and nonlinearity in strength development, in this paper the computational intelligence techniques such as multilayer perceptron (MLP), radial basis function (RBF) and genetic programming (GP) are used to develop a mathematical model. To generate the mathematical model, an experimental study is conducted to obtain the strength of three inland soils namely, red earth (CL), brown earth (CH) and black cotton soil (CH) for different water contents, cement contents and curing periods. In order to generate a generic mathematical model using computational intelligence techniques, two saline soils (Ariake clay-3 and Ariake clay-4) and three inland soils are used. A detailed study of the relative performance of the computational intelligence techniques and the empirical model has been carried out. 2006 Elsevier Ltd. All rights reserved. Keywords: Soft ground; Saline soil; Cement stabilization; Empirical model; Multilayer perceptron; Radial basis function; Genetic programming; Unconfined compressive strength
1. Introduction Due to increased urbanization and industrialization, foundation soils with suitable geotechnical properties are not easily available. There is a growing demand for reclamation of marine and inland soft soils/grounds, which are characterized by high plasticity, higher fraction of fines and void ratio, with low strength and high compressibility [1]. In many cases, the natural water content of soft ground is high and at times even higher than the liquid limit of the soils. For such soils, deep mixing is a popular method for improving the strength [2]. In deep mixing, the cement grout or dry cement is injected into the natural soil at the *
Corresponding author. Tel.: +91 80 22932672; fax: +91 80 23600404. E-mail address: [email protected] (P.V. Sivapullaiah).
0266-352X/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2006.03.006
depth required and a blade is pushed into the ground to mix the soil and cement [3]. The extent of the strength improvement depends on the mineralogy, environmental conditions of the soft ground, curing period and the type and amount of cement used [1]. Porbaha [1], Tan et al. [3], Uddin et al. [4], Yin and Lai [5], Miura et al. [6], Suksun [7] and Horpibulsuk et al. [8,9] have carried out experimental studies to understand the strength improvement in soft ground using cement stabilization techniques. In the stabilization of soft ground, Miura et al. [6], Tan et al. [3] and Horpibulsuk et al. [9] have observed that the strength development depends on the initial clay water content and the cement content. For optimal and effective utilization of cement, a predictive model is required to relate the gain in strength of soil with time in terms of important variables involved. Tan
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
et al. [3] have established an empirical relationship to predict the strength development based on cement content, water content and curing period. The strength developed at specific cement and water contents and curing period is used as the reference compressive strength for a given soil and the strength developed under other conditions for the same soil is normalized using the reference strength. Their empirical model assumes linear variation of normalized strength and the model is different under different conditions of the same soil. The model also assumes that the ratio of the strength developed and the reference strength is constant for any soil if the reference conditions are the same. These assumptions will affect the accuracy of the strength prediction. However, it is mentioned in [3] that a careful statistical study needs to be undertaken to generalize the relationship to improve the performance of the model. One way to improve the performance of the empirical model is to consider the effects of other parameters. Miura et al. [6], Suksun [7] and Horpibulsuk et al. [9] proposed a parameter called the clay water cement ratio (wc/C), which is analogous to the water cement ratio used in concrete technology, which is the ratio of the initial water content of the clay (wc) to the cement content (C). Miura et al. [6] observed that the water cement ratio is the principal parameter for analysis of strength and deformation behavior of cement stabilized soft ground. Using the clay water cement ratio (wc/C), as the principal parameter Suksun [7] and Horpibulsuk et al. [9] developed the normalized empirical model to predict the strength of cement stabilized soft ground. The model requires the reference strength and assumes linear variation of the logarithm of the curing period. However, the empirical model is localized and different equations are considered for different soils depending on the liquidity index. The above limitations of the empirical models suggest the necessity to develop a generic mathematical model for strength prediction. Computational intelligence techniques are capable of approximating the nonlinear input–output relationship for a wide range of applications. Hence in this paper, computational intelligence techniques such as multilayer perceptron (MLP) network, radial basis function (RBF) network and genetic programming (GP) techniques are used to approximate the functional relationship between the various parameters and strength developed. Neural networks (NNs) have been used in geotechnical engineering for estimation of unsaturated shear strength [10], predicting settlement of shallow foundations [11–13], profiling stress history of clays [14], modeling of compaction curve [15], site characterization [16], triaxial compression behavior of sand and gravel [17], predicting the capacity of driven piles in cohesion less soils [18], predicting settlements during tunneling [19] and stress stain modeling of sands [20]. A collection of neural network applications in geotechnical engineering are presented in [21,22]. NNs gained prominence in the 1980s due to their ability to learn from past data, and their ability to generalize
197
knowledge obtained during learning. A variety of neural network architectures and associated learning algorithms are available in the literature for various types of applications. Of these, two neural network models namely, MLPs [23] and RBFs [24] are considered. MLPs are feed forward networks with a layered structure. In general, MLPs have an input layer, one or more hidden layers, and one output layer. Each of the nodes in a given layer are connected to the other nodes in the next layer. There are no connections within layers. In most of the applications a feed forward network with a single hidden layer is used with sigmoidal activation functions for the units. The back propagation (BP) algorithm is widely used for training the MLP networks. Although MLPs can be used with any number of hidden layers, it has been shown that only one hidden layer suffices to approximate any function with discontinuities to an arbitrary degree of precision, provided the activation function of the hidden units are nonlinear (the universal approximation theorem) [23]. The RBF network model also consist of three layers; input, hidden and output layers. The nodes within each layer are fully connected to the previous layer. The input variables are each assigned to nodes in the input layer, and pass directly to the hidden layer without weights. The hidden nodes or units contain the radial basis functions. The RBF is similar to the Gaussian density function, which is defined by parameters namely, a ‘center’ and a ‘width’ [24]. The Gaussian function gives the highest output when the incoming variables are closest to the center position and the output decreases monotonically as the distance from the center increases. The ‘width’ of RBFs controls the rate of decrease in the output. Even though RBFs are well suited for function approximation, for large sets of data the problem of ‘curse of dimensionality’ becomes predominant [24], i.e., the increase in number of features reduces the approximating ability of the RBF network. In NNs, the basic limitation is that the optimal configuration is not known a priori. Also, often it is not possible to express the relationship between the different parameters governing a phenomenon because the knowledge represented in the weights is generally opaque i.e., the weights do not reveal any information about the input–output relationship. One of the recent developments in computational intelligence techniques is genetic programming (GP) which has the ability to approximate the relationship between the inputs and outputs and express the relationship mathematically [25]. The generated mathematical expression may also provide an insight into the understanding of the problem. Hence in this paper, the MLP, RBF and GP techniques are used to develop a generic model for strength prediction. In order to generate such a model, three inland soils for which experimental studies are carried out, and two saline soils available in the literature [7] are considered. For these soils, empirical models are developed using the concept presented by Horpibulsuk et al. [9] and the performance is compared with the generic models developed using the computational intelligence technique discussed above.
198
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
2. Experimental procedure This section contains description of the materials used in the investigation and the experiments carried out for the determination of the unconfined compressive strength of cement stabilized brown earth, black cotton soil and red earth under different water and cement contents and curing periods. 2.1. Materials used Cement. Portland cement (43 Grade) was used as the stabilizer in the investigation. Clay water cement ratios of 5, 10, 15 and 20 were considered in the present study. The cement content required to prepare the specimen depends on the clay water cement ratio and initial clay water content. Soils. Three Indian soils namely, red earth, classified as CL (lean clay), brown earth and black cotton soil, both classified as CH (fat clay) are used in the study. Red earth is obtained from Indian Institute of Science campus, Bangalore, Karnataka state, India, Brown earth is obtained from a construction site on Airport road in Bangalore, Karnataka state, India and Black cotton soil is obtained from Bagalkot, Karnataka state, India. All the soils are collected by open excavation from a depth of 1 m from the natural ground. The geotechnical properties of the soils are given in Table 1. 2.2. Sample preparation The samples are prepared using split mould and are of size of 38 mm diameter and 76 mm length. Initially, the soil and de-ionized water are mixed to form a soil paste. Initial water content equal to liquid limit, 1.5 times the liquid limit and 2 times liquid limit are considered for specimen generation. For hydration process, the cement consumes water equal to 40% of its mass [26]. Hence, the cement paste is
2.3. Experimental results The specimens are cured for curing period (D) of 7, 14, 28 and 56 days. The unconfined compressive strength of specimens is determined using ASTM D2166-91 method [27]. Table 2 presents the results of unconfined compressive strength (kPa) of brown earth stabilized with cement for different initial clay water content (wc), clay water cement ratio (wc/C) and curing period. Similarly, the results for Table 2 Unconfined compressive strength (kPa) of brown earth stabilized with cement wc/C
Unconfined compressive strength (kPa) 7 Days
14 Days
28 Days
Water content, wc = liquid limit water content (%) 5 1182 1578 2552 10 569 585 656 15 150 225 232 20 15 49 69
56 Days 2699 736 284 98
Water content, wc = 1.5 · liquid limit water content (%) 5 438 644 846 10 279 419 497 15 161 171 193 20 58 62 74
932 595 238 93
Water content, wc = 2 · liquid 5 145 10 127 15 99 20 66
390 374 203 95
limit water content (%) 231 351 164 277 135 184 69 94
Table 3 Unconfined compressive strength (kPa) of black cotton soil stabilized with cement
Table 1 Geotechnical properties of soils Properties
prepared separately by mixing cement with de-ionized water equal to 0.4 times the mass of cement. Finally, both the soil paste and the cement paste are mixed thoroughly and put into the split mould. The specimens are cured in desiccators at 100% humidity, ensuring the availability of water for cementitious reactions.
Red earth
Brown soil
2.64
2.63
Specific gravity Atterberg limits (%) Liquid limit Plastic limit Plasticity index Shrinkage limit Free swell index
38 15 23 16 25
60 23 37 13 36
Standard proctor test Maximum dry density (kN/m3) Optimum moisture content (%)
16.4 18
14.3 28
Grain size distribution (%) Fine sand Silt Clay Classification
40 28 32 CL
20 34 46 CH
Black cotton soil 2.76 97 35 62 9 138 11.5 37 2 37 61 CH
wc/C
Unconfined compressive strength (kPa) 7 Days
14 Days
28 Days
56 Days
Water content, wc = liquid limit water content (%) 5 1842 2255 2636 10 515 630 781 15 120 207 250 20 82 117 155
3702 981 330 255
Water content, wc = 1.5 · liquid limit water content (%) 5 276 425 631 10 92 235 390 15 73 127 154 20 58 69 91
1136 500 223 119
Water content, wc = 2 · liquid limit water content (%) 5 122 284 343 10 79 132 212 15 77 99 135 20 51 75 89
548 334 166 105
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208 Table 4 Unconfined compressive strength (kPa) of red earth soil stabilized with cement wc/C
14 Days
28 Days
Water content, wc = liquid limit water content (%) 5 476 683 857 10 141 166 196
56 Days 1572 311
Water content, wc = 1.5 · liquid limit water content (%) 5 116 166 349 10 112 147 152
610 284
Water content, wc = 2 · liquid 5 88 10 76 15 61 20 34
552 310 196 95
limit water content (%) 158 323 112 161 79 128 60 85
black cotton soil and red earth are given in Tables 3 and 4, respectively. 3. Empirical formula The empirical model for brown earth, black cotton soil and red earth are developed using the concept presented in [9] are discussed below. 3.1. Development of empirical formula Suksun [7] and Horpibulsuk et al. [9] used Abram’s law as the basis for the model development. Abram’s law is extensively used in concrete technology and can be expressed as: A ; ð1Þ Bðw=CÞ where Q is the compressive strength of the concrete at a specific age, w/C is the water cement ratio by weight, and A and B are empirical constants. Using Abram’s law, a similar parameter, called the clay water cement ratio (wc/C) is derived and Abram’s law is modified as: Q¼
Q¼
A . Bðwc =CÞ
Qðwc =CÞ2 Qðwc =CÞ1
Unconfined compressive strength (kPa) 7 Days
Qðwc =CÞ1
ð2Þ
The modified power fit is used for the plot of the clay water cement ratio versus curing periods for different liquidity indices (LI). Using these plots, the values of A and B are calculated. The value of A varies over wide range and it depends on the type of clay, the liquidity index and curing period. However, the variation of the parameter B is only between 1.22 and 1.24, irrespective of the type of clay and the curing period. Hence, the value of parameter B is selected as 1.24. The parameter A has a specific value for a clay at a specified curing period (D) and liquidity index (LI) but does not change with clay water cement ratio. The parameter A is eliminated by taking the ratios of strength developed at different clay water cement ratios as follows:
Qðwc =CÞ2
¼
A=Bðwc =CÞ1 ¼ Bfðwc =CÞ1 ðwc =CÞ2 g ; A=Bðwc =CÞ2
¼ 1:24fðwc =CÞ1 ðwc =CÞ2 g ;
199
ð3Þ ð4Þ
where Qðwc =CÞ1 is the unconfined compressive strength at wc/C1 and Qðwc =CÞ2 is the unconfined compressive strength at wc/C2. Using the strength developed at the curing period (D = 28) as a reference point, the normalized strength for different types of soft ground, stabilized at different liquidity indices is calculated. The emprical relationship between the normalized strength and curing period for different liquidity indices is expressed as 0:038 þ 0:281 lnðDÞ if LI ¼ 1:0–2:5; QD ¼ ð5Þ Q28 0:126 þ 0:342 lnðDÞ if LI > 2:5; where QD is the strength to be estimated at D days of curing, and Q28 is the strength at the reference curing period (D = 28). Using the concept explained in [7], the empirical model is developed using Eqs. (4) and (5) 8 1:24fðwc =CÞ2 ðwc =CÞ1 g ð0:038 þ 0:281 lnðDÞÞ > > > < Qðwc =CÞ1;D if LI ¼ 1:0–2:5; ¼ Qðwc =CÞ2;28 > 1:24fðwc =CÞ2 ðwc =CÞ1 g ð0:126 þ 0:342 lnðDÞÞ > > : if LI > 2:5; ð6Þ where D is the curing period in days, Qðwc =CÞ1;D is the unconfined compressive strength at wc/C1 for the curing period = D, and Qðwc =CÞ2;28 is the unconfined compressive strength at wc/C2 for the reference curing period = 28 days. 3.2. Limitations of empirical formula Using the empirical model given in Eq. (6) and the reference strengths highlighted in Table 5, the unconfined compressive strength of black cotton soil has been predicted. The predicted and actual strength for black cotton soil for various water cement ratios and curing periods are listed in Table 5. The root mean square error deviation from actual strength is 226 kPa. The higher error may be due to the fact that certain parameters, such as sodium salts and pH, which play an important role in cement stabilization [28,29], are not considered. Also, the empirical model is based on a best fit of the normalized strength for different soils. However, the reference used for predicting the strength in Table 5 is different for each soil. It can be observed in Table 5 that two separate references have been used in the prediction of strength for liquidity indices of 1.0 and 1.8. Even though the liquidity index is less than 2.5, the use of a single reference strength in the equation is found to result in large errors. For example consider the 16th data sample in Table 5. By taking 781 kPa as the reference, the prediction is 106 kPa, whereas a reference of
200
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
Table 5 Actual and predicted unconfined compressive strength using empirical model for black cotton soil Serial number
Curing period (days)
Clay water cement ratio
Cement content (%)
Liquid limit (%)
Liquidity index
Water content (%)
Actual UCS (kPa)
Predicted UCS (kPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
7 7 7 7 14 14 14 14 28 28 28 28 56 56 56 56 7 7 7 7 14 14 14 14 28 28 28 28 56 56 56 56 7 7 7 7 14 14 14 14 28 28 28 28 56 56 56 56
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
19.40 9.70 6.47 4.85 19.40 9.70 6.47 4.85 19.40 9.70 6.47 4.85 19.40 9.70 6.47 4.85 29.10 14.55 9.70 7.28 29.10 14.55 9.70 7.28 29.10 14.55 9.70 7.28 29.10 14.55 9.70 7.28 38.80 19.40 12.93 9.70 38.80 19.40 12.93 9.70 38.80 19.40 12.93 9.70 38.80 19.40 12.93 9.70
97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55
97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 194 194 194 194 194 194 194 194 194 194 194 194 194 194 194 194
1842 515 120 82 2255 630 207 117 2636 781 (Ref) 250 155 3702 981 330 255 276 92 73 58 425 235 127 69 631 390 (Ref) 154 91 1136 500 223 119 122 79 77 51 284 132 99 75 343 212 (Ref) 135 89 548 334 166 105
1339 457 156 53 1785 609 208 71 2231 761 260 89 2677 913 311 106 669 228 78 27 891 304 104 35 1114 380 130 44 1337 456 156 53 335 114 39 13 483 165 56 19 630 215 73 25 777 265 90 31
390 kPa results in 53 kPa. This implies that the empirical model is not really generalized and more fine tuning the equation is required to improve the performance of each and every soil. 3.3. Parameter selection In order to obtain a generic model for unconfined compressive strength prediction, saline soils such as Ariake clay-3 and Ariake clay-4 presented in [7] are also considered. From the limitation of the empirical model, one can observe that the other parameters such as liquid limit
(LL), liquidity index (LI), pH and sodium ion concentration (Na+) play an important role in strength development. The liquid limit is one of the very important soil properties, which has been extensively correlated with the behaviour of soils. Liquidity index can indicate the interparticle cementation of the soil particles. The pH will affect the cementitious reaction to greater extent. It is well known that if the pH is greater than 12, the amorphous silica present in the clay gets dissolved and reacts with calcium hydroxide to form calcium–silicate–hydrate gel, which imparts strength to the stabilized soil [30]. Since, the salinity also affects the strength development the
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
sodium ion concentration present in soil is used to indicate the salinity in the soil. Hence, the liquid limit (LL), liquidity index (LI), water content (wc), cement content (C), curing period, (D), pH and sodium ion concentration (Na+) are used to develop the model. 4. Computational intelligence techniques In this section, the two popular neural network architectures, multilayer perceptron, radial basis function network, and genetic programming are described briefly. 4.1. Multilayer perceptron network Multilayer perception provides a general, nonlinear parameterized mapping between a set of inputs and a set of outputs. A MLP network with three layers and sigmoidal activation functions can approximate any smooth mapping, if sufficient number of neurons is included in the hidden layer [23]. The general architecture of MLP network is shown in Fig. 1. The MLP network is trained using supervised back propagation learning algorithm [23]. The functional relationship in the strength development is established through MLP network. The learning algorithm is designed to minimize a cost function J of the form
J¼
1XX i ðy k ^y ik Þ2 ; 2 i k
201
ð7Þ
where y ik and ^y ik are the desired output and the actual output vectors of the kth output of ith sample, respectively. Let the input and output vectors are denoted by u and ^y , respectively. The weight connection between the ith neuron in the input layer and the jth neuron in the hidden layer is represented as w1ij . Similarly, let the weight connection between the ith neuron in the hidden layer and the jth neuron in the output layer be w2ij . The output of the ith neuron in the hidden layer is described mathematically as ! X 1 hi ¼ f ðneti Þ ¼ f wij uj ð8Þ j
and the output of the ith neuron in the output layer is described mathematically as ! X 2 ^y i ¼ f ðneti Þ ¼ f ð9Þ wij hj ; j
where f(Æ) is the sigmoidal activation function. The interconnection weights are adjusted such that the cost function J can be minimized, and are given by
Fig. 1. Architecture of multilayer perceptron.
202
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
w1ij ¼ w1ij þ gd1i uj ; w2ij
¼
w2ij
þ
gd2i hj ;
ð10Þ
where dli is error term of jth neuron in lth layer and g is learning rate. The error terms are defined below: d2i ¼ ðy i ^y i Þf 0 ðnet2i Þ; ! X w2ij d2i f 0 ðnet1j Þ; d2j ¼
ð11Þ
i
where f 0 (Æ) is the first derivative with respective to their arguments. The training process is terminated when the cost function is sufficiently small. The details of the learning algorithm are given in [23]. 4.2. Radial basis function network The radial basis function network consists of three layers: an input layer, a kernel layer and an output layer as shown in Fig. 2. The neurons in input layer are assigned to kernel neurons in the hidden layer without any modification, i.e., the weight connection is unity. The kernel layer nonlinearly transforms the inputs to a higher dimension. The activation function in the kernel neurons is similar to the Gaussian function, which is defined by a ‘centre’ position and a ‘width’ parameter. The Gaussian function gives the highest output when the incoming variables are closest to the centre position and decreases monotonically as the distance from the centre increases. The width of the Gaussian function controls the rate of decreases. Let x be the input vector of dimension n. The output (hi) of the ith kernel neuron is given by , ! n X 2 hi ¼ exp ½xk xik r2i i ¼ 1; 2; . . . ; H ; ð12Þ k¼1
where H is the number of kernel neurons, xik is centre of the ith kernel neuron for the kth input variable and ri is the width of the ith kernel neuron. The output value of the ith output neuron is equal to the summation of the weighted outputs of the kernel neurons
Fig. 2. Architecture of radial basis function network.
and bias term of the output neuron, and is described mathematically as ^y i ¼
H X
wij hj þ wio ;
i ¼ 1; 2; . . . ; O;
ð13Þ
j¼1
where O is the number of output neurons and wio is bias value of the ith output neuron. The parameters of the RBF network are determined using a three step training scheme [24]. First, the kernel neurons centers are determined using the k-means clustering algorithm. Then, the width parameters are determined by a nearest-neighbor method. Finally, weights connecting kernel neurons and output neurons are calculated using multiple linear regression techniques. The details of the learning algorithm are given in [24]. 4.3. Genetic programming Genetic programming (GP) is a branch of genetic algorithms. The main difference between genetic programming and genetic algorithms is the representation of the solution. Genetic programming creates computer programs in the LISP or SCHEME computer languages as the solution, where as genetic algorithms create a string of numbers that represent the solution. GP has gained prominence in the 1990s and it has an ability to understand the underlying relationship between the input and outputs and express them mathematically. GP is provided with evaluation data, a set of primitives and fitness functions. The evaluation data describe the specific problem in terms of the desired inputs and outputs. They are used to generate the best computer program to describe the relationship between the input and output very well [25]. The following four steps are used to generate a mathematical relationship between the input/output data using: 1. Generate an initial population of random compositions of the functions and terminals of the problem (computer programs). 2. Execute each program in the population and assign it a fitness value according to how well it solves the problem. 3. Create a new population of computer programs. Copy the best existing programs Create new computer programs by crossover and mutation. 4. The best computer program that appeared in any generation, the best-so-far solution, is designated as the result of genetic programming. In GP, the best idea is to maximize the user defined fitness function. Thus, at the end of each GP run, a current population (computer program) with good fitness, and best fit individual is selected as solution to the given problem. The following subsections present a brief description on various components of GP.
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
4.3.1. Crossover operation The crossover operation is used to create new offspring computer programs from two parental computer programs selected based on the fitness value. The parental computer programs are typically of different sizes and shapes. The offspring computer programs are composed of subexpressions (subtress subprograms, subroutines and building blocks) from their parents, and typically of different sizes and shapes from their parents. The creation of offspring from the crossover operation is accomplished by swapping random fragments of the parents. Fig. 3 represents a typical crossover operation. 4.3.2. Mutation operation Mutation operation is shown in Fig. 4. Mutation operation protects the model against premature convergence and improves the non-local properties of the search. Mutation is implemented by randomly removing a subtree at a selected points and replacing it with randomly generated subtree. 4.3.3. Tournament selection Tournament selection is a process of selecting the computer program for genetic operations. Typically, the genetic programming implementation uses a model in which certain randomly selected computer programs in a subgroup compete and the fittest is selected. The tournament size defines the number of randomly selected computer programs for selection process.
203
Fig. 4. Mutation operation.
4.3.4. Fitness function GP is guided by fitness function to search for the most efficient computer program to solve a given problem. The negative of mean square error is selected as a simple measure of fitness in function approximation problem Fitness ¼
N 1 X 2 ðY i Y^ i Þ ; N i¼1
ð14Þ
where Y is the actual value and Y^ is the predicted value. 4.3.5. Functions and terminals The terminal and function sets are also important components of genetic programming. The terminal and function sets are the alphabets of the programs to be made. The terminal set consists of the variables and constants
Fig. 3. Crossover operation.
204
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
Table 6 List of parameters used in GPQUICK software Parameters
Values
Population size GP generations Crossover weight Mutation weight Mutation rate Tournament size Tree length
5000 1,000,000 70 20 60 6 25
of the programs and the function set consists of the functions such as addition (ADD), subtraction (SUB), multiplication (MUL), division (DIV), square root (SQ) and the natural logarithm (LN). 4.3.6. Termination criterion Koza [25] has shown that in GP the evolution is a neverending process, and hence a termination criterion is needed. The termination criterion for GP is generally based on the problem or is limited by the number of generations. In GP, a user-defined fitness function has to be maximized for his/her application. Thus at the end of a GP run, a current population of individuals (computer programs), and also the fittest individual (computer program) that appeared during the run is obtained. The fittest individual that has evolved for the given problem is its solution or desired mathematical model. GPQUICK software [31] has been used for this function approximation problem. In the GP runs, a maximum number of generations (1,000,000) or MSE is equal to 0.002 is used as the criterion to stop the program. The choice of parameters for GPQUICK software used in the runs is shown in Table 6. The parameters for genetic programming are selected using trial and error procedure. The parameter selection will affect convergence of GP [25].
computational intelligence techniques. Training error (TE) is the root mean sum of squared residuals (error) in the training data and generalization/testing error (GE) is the root mean squared error in the (validating) testing data vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N1 u1 X ^ i Þ2 ; TE ¼ t ð15Þ ðQ Q N 1 i¼1 i vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N2 u1 X ^ i Þ2 ; ð16Þ ðQ Q GE ¼ t N 2 i¼1 i ^ i the estiwhere Qi is the actual value of strength, and Q mated output of the computational intelligence technique for the ith input. N1 is the number of data points used in the training set and N2 is the number of data points used in the testing set. The total error (TOE) takes into account both the training and generalization error, i.e., sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N 1 TE2 þ N 2 GE2 TOE ¼ . ð17Þ ðN 1 þ N 2 Þ Since the TE, GE and TOE will not reveal substantial statistical information about the difference between the actual and the predicted value, we use statistical estimate like R2 as one of the performance measure. R2 is the square of correlation between the actual value and the predicted value or it is an indicator of how much the variation in the data is explained by the model [32]. It is calculated using the following equation: R2 ¼ 1
SS T ; SS R
ð18Þ
where SS R ¼
N X
2
ðQðtÞ QÞ ;
ð19Þ
N X 2 ^ ðQðtÞ QðtÞÞ ;
ð20Þ
t¼1
4.3.7. Function approximation using GP Koza [25] has applied GP for a function approximation problem (regression problem). In function approximation problem a single GP expression is evolved for each output. This GP expression evolved for each output is called as genetic programming approximator expression (GPAE). This GPAE approximates the functional relationship between the input–output data set. 5. Computational intelligence based unconfined compressive strength prediction In this section the procedures to develop a generic model to predict the unconfined compressive strength using MLP, RBF and GP are presented. The performances of these models with the empirical model are also compared. 5.1. Performance measures The training error, generalization error and total error are measured to study the qualitative performance of the
SS T ¼
t¼1
^ is the predicted where Q is the mean value of strength and Q strength. R2 can take any value between 0 and 1, with a value closer to 1 indicating a better fit [32]. For example, an R2 value of 0.85 means that the approximation explains 85% of the total variation in the data about the average. 5.2. Selection of training data To predict the unconfined compressive strength (kPa), the following eight input parameters are selected for the computational intelligence techniques, curing period (D), clay water cement ratio (wc/C), cement content (C) in percentage, liquid limit (LL) in percentage, liquidity index (LI), water content (wc) in percentage, pH and sodium ion concentration (Na+) in gram/litre. The input soil characteristic and the output strength are stored in a file for different soils. The 186 data are recorded for various conditions and these data are used for the model develop-
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
ment using computational intelligence techniques. Out of 186 data, 154 data are selected randomly for training and 32 data for testing the generalization capability of the computational intelligence techniques. 5.3. Model generation using neural networks The details of the training process adopted for neural networks and optimal network selection procedure are described below. 5.3.1. Multilayer perceptron The input to the network is eight soil characteristics and one bias. The output of the network is the strength developed in the soil. The network architecture is represented as N n1 ;n2 ;n3 , where n1 is the number of input nodes, n2 is the number of hidden nodes and n3 is the number of output nodes. A learning rate, g, of 0.0003 was used for training the MLP network. The network is trained for 75,000 epochs or training error less than 0.002. The optimal network configuration of MLP is N9,22,1. The optimal configuration of network to approximate the input–output relationship is found using the steps explained in [32]. The optimal configuration is tested with the same training data and the strength predicted and the actual strength are shown in Fig. 5. A few training points are selected for this purpose. From the figure, it can be observed that the predicted strength closely matches with the actual values at various sample points. 5.3.2. Radial basis function network Similar to the MLP network, the input to the radial basis function network are eight soil parameters and output of the network is strength. The network is represented as N n1 ;n2 ;n3 ;p , where n2 is number of radial basis function neu-
205
rons in the network and p is the overlap parameter. All eight inputs are used to calculate the centre and width of the RBF network using the k-means clustering algorithm and linear weight matrix between the RBF neurons and output neuron is calculated using regression procedure. The optimal RBF configuration is N8,97,1,95. Here also, the optimal configuration of network to approximate the input–output relationship is found using the steps explained in [32]. The trained RBF network is tested with same training data. The predicted and the actual strength are plotted in Fig. 5. From the above figure, it can be observed that the RBF network approximates the input–output relationship very well. 5.4. Model generation using genetic programming The same training data are used to generate the genetic programming expressions. The best GP program is evaluated until it reached fitness measure (0.002) or number of population is equal to 1,000,000. Four runs of GP with the training set were used and obtained the best computer program evolved at each run. At the end of each run, the best computer program evolved is in the form of LISP expression and this expression can be easily converted into a mathematical expression. Initially for the first run, the tree length of LISP expression is restricted to 1000. The resulting LISP expression is able to predict the strength development with considerable accuracy, but reducing the LISP expression to obtain mathematical equations is very difficult. The performance of the LISP expression is given in Table 7. Even though the performance of the above LISP expression is better, it is very difficult to study the effect of various parameters in strength development. Hence for the remaining three GP runs, the tree length was restricted to 25 such that the resulting LISP expression can be simplified to obtain a mathematical relationship.
4000 Actual GP RBF MLP
Unconfined compressive strength (kPa)
3500
3000
2500
2000
1500
1000
500
0 0
5
10
15
20
25
30
35
40
Data points Fig. 5. Comparison of trained data for different computational intelligence techniques with actual data.
206
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
Table 7 Performance measures for empirical and computational intelligence techniques based models Maximum error
Empirical model Multilayer perceptron Radial basis function GP (tree length = 25) GP (tree length = 1000)
Training
Testing
Total
Training
Testing
Total
Training
Testing
Total
– 221 413 505 251
– 178 420 572 324
1025 221 420 572 324
– 50 136 151 72
– 71 133 195 97
226 54 135 159 77
– 0.9948 0.9616 0.9527 0.9893
– 0.9914 0.9696 0.9345 0.9837
0.9080 0.9941 0.9632 0.9490 0.9881
Out of the three runs with restricted tree length, the best LISP expression obtained is given below: (ADD (MUL 35 B7) (SUB (DIV (MUL (MUL B1 B3) (DIV B4 B5)) (ADD (DIV (MUL B2 (DIV B2 B1)) B3) (ADD (DIV (MUL B7 (DIV (DIV (DIV B8 B5) B5) B5)) B3) (DIV B5 (DIV (MUL B7 (DIV 26 B3)) B3))))) (MUL (ADD (DIV 19 B8) (ADD (SUB (MUL B8 B5) (MUL B8 B3)) (DIV B3 (SUB 19 (MUL B8 B3))))) (MUL (MUL B8 B3) (DIV B2 B7))))) where B1 is the curing period (D) in days, B2 the clay water cement ratio (wc/C), B3 the cement content (C) in percentage, B4 the liquidity index (LI), B5 the water content (wc) in percentage, B6 the pH, and B7 is the sodium ion concentration (Na+) in gram/litre. The corresponding GPAE is after simplification is ðCÞðLLÞLnðDÞ i ðpHÞðNaþ Þ ðCÞ2 ðLIÞ þ ðCÞðLIÞ3 þ 26ð LI pHÞ ðwc =CÞðCÞðNaþ Þ 19 þ ðLIÞðNaþ Þ ðCÞðNaþ Þ ðpHÞ Naþ C þ . ð21Þ ð19 ðCÞðNaþ ÞÞ
Q ¼ 35ðpHÞ þ
R2 estimate
Root mean square error
h
ðwc =CÞ2 ðCÞLnðDÞ
The above expression is complex and is difficult to obtain any physical insight about the behavior of the soils in strength development. On the other hand, the mathemati-
cal expression is able to predict the strength development accurately over wide range of soils. The performance of the best LISP expression is given in Table 7. From Table 7, it is evident that the prediction based GP model developed using unrestricted tree length has better performance than the restricted tree length. The predicted strength and the actual strength are shown in Fig. 5. From the figure, it can be observed that the GP is able to approximate the underlying the relationship very well. 5.5. Performance of computational intelligence techniques models The quantitative performance of computational intelligence techniques for training data set is given in Table 7. From the table it can be observed that the TE for MLP is less than the RBF and GP. The same is reflected in the R2 estimate. 6. Prediction of unconfined compressive strength The MLP, RBF, GP and empirical models are used to predict the strength development in soft ground using cement stabilization technique and the performance of these models is compared. The predicted values of strength
3500 Actual GP RBF MLP Empirical
Unconfined compressive strength (kPa)
3000
2500
2000
1500
1000
500
0 0
5
10
15
20
25
30
Data points Fig. 6. Actual strength and predicted strength using empirical model and computational intelligence techniques.
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
for various experimental data are shown in Fig. 6. From the figure, one can observe that it is not possible to assess the accuracy of different predictive methods. From the performance measures listed in Table 7, it can be observed that for the unconfined compressive strength predicted using the empirical model, the maximum error equal to 1025 and the RMSE equal to 226. The errors are the higher when compared to the evolutionary computational intelligence techniques. Also the R2 estimate = 0.9080 is the lowest for empirical model. Hence, the prediction using computational intelligence techniques is more accurate than the conventional empirical model. Table 8 gives the details of some of the testing data points used for the comparative study of the efficiency of different techniques. In the computational intelligence techniques best prediction is obtained using MLP, with maximum error = 221, TOE = 54 and total R2 estimate = 0.9941. Even though the MLP and RBF perform better than the GP, the knowledge gained by the network weights is opa-
207
que. Hence, one cannot assess the role played by the different input parameters in the strength development. The mathematical expression generated using the GP has a potential to assess different parameters in strength development. The method has the advantage of continuous up-gradation when the new soil data are available, which improve the reliability of the model. 7. Conclusions A generic mathematical model is developed to predict the unconfined compressive strength of cement treated soils using multilayer perceptron, radial basis function network and genetic programming techniques for two saline soils in the literature and three inland soils. The quantitative performance of the computational intelligence techniques is measured and compared with the existing empirical model. The performance measures show that the MLP network has better performance followed by GP (tree length =
Table 8 Details of parameters of different soils used as testing data Liquidity index
Water content (%)
pH
Na+ (g/l)
UCS (kPa)
60 60 60 60 60 60 60
1 1 1 1.8 1.8 2.6 2.6
60 60 60 90 90 120 120
8.6 8.6 8.6 8.6 8.6 8.6 8.6
0.034 0.034 0.034 0.034 0.034 0.034 0.034
585 49 2552 171 846 69 277
9.7 19.4 4.9 29.1 7.3 19.4 12.9 38.8 12.9
97 97 97 97 97 97 97 97 97
1 1 1 1.8 1.8 2.6 2.6 2.6 2.6
97 97 97 145.5 145.5 194 194 194 194
8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8
0.159 0.159 0.159 0.159 0.159 0.159 0.159 0.159 0.159
630 2636 155 631 91 132 99 343 135
10 5 10 10 20 5
3.8 7.7 5.8 7.7 3.8 15.3
38 38 38 38 38 38
1 1 1.8 2.7 2.7 2.7
38 38 57 76 76 76
8.3 8.3 8.3 8.3 8.3 8.3
0.013 0.013 0.013 0.013 0.013 0.013
166 857 152 112 60 323
Ariake clay 3 28 7 28 7 14 28 14 14
7.1 7.1 8.7 8.7 10.7 10.7 12.3 9.3
15 15 15 15 15 15 15 20
106 106 106 106 106 106 106 106
1 1 1.5 1.5 2 2 2.5 2.5
106 106 130 130 160 160 185 185
8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5
1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6
1300 2737 1000 1867 589 857 714 1228
Ariake clay 4 7 7
12.2 10.1
10 15
122 122
1 1.5
122 152
8.8 8.8
1 1
661 1286
Curing period (days)
Clay water cement ratio
Cement content (%)
Brown earth 14 14 28 14 28 14 28
10 20 5 15 5 20 10
5.9 3.0 11.8 5.9 17.7 5.9 11.8
Black cotton soil 14 28 28 28 56 14 28 28 14
10 5 20 5 20 10 15 5 15
Red earth 28 28 14 14 28 56
Liquid limit (%)
208
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
1000), RBF, GP (tree length = 25) and empirical model. Even though the MLP and RBF perform very well, it is not able to express the relationship between the different parameters governing the strength development because the knowledge gained by the weights is often opaque. But the GP has the ability to learn the underlying input–output relationship and express them mathematically (GPAE). One can use these GP expressions to understand the influence of input parameters for strength development. References [1] Porbaha A. State of the art in deep mixing technology: Part I. Basic concepts and overview. Ground Improvement 1998;2(2):81–92. [2] Horpibulsuk S, Miura N, Koga H, Nagaraj TS. Analysis and strength development in deep mixing: a field study. Ground Improvement 2004;8(2):59–68. [3] Tan TS, Goh TL, Yong KY. Properties of singapore marine clays improved by cement mixing. Geotech Testing J ASTM 2002;25(4): 422–33. [4] Uddin K, Balasubramaniam AS, Bergado DT. Engineering behavior of cement treated Bangkok soft clay. Geotech Eng J 1997;28(1):89–119. [5] Yin JH, Lai CK. Strength and stiffness of Hong Kong marine deposit mixed with cement. Geotech Eng J 1998;29(1):29–44. [6] Miura N, Horpibussuk S, Nagaraj TS. Engineering behavior of cement stabilized clay at high water content. Soils Foundations 2001;41(5):33–45. [7] Suksun Horpibulsuk. Analysis and assessment of engineering behavior of cement stabilised clays. PhD Thesis, Saga University, Japan; 2001. [8] Horpibulsuk S, Bergado DT, Lorenzo GA. Compressibility of cementadmixed clays at high water content. Geotechnique 2004;54(2):151–4. [9] Horpibulsuk S, Miura N, Nagaraj TS. Assessment of strength development in cement admixed high water content clays with Abrams law as basis. Geotechnique 2003;53(2):439–44. [10] Lee SJ, Lee SR, Kim YS. An approach to estimate unsaturated shear strength using artificial neural network and hyperbolic formulation. Comput Geotech 2003;30(6):489–503. [11] Shahin MA, Maier HR, Jaksa MB. Predicting settlement of shallow foundations using neural networks. J Geotech Geoenviron Eng ASCE 2002;128(9):785–93. [12] Shahin MA, Maier HR, Jaksa MB. Data division for developing neural networks applied to geotechnical engineering. J Comput Civil Eng ASCE 2004;18(2):105–14. [13] Shahin MA, Maier HR, Jaksa MB. Settlement prediction of shallow foundations on granular soils using B-spline neurofuzzy models. Comput Geotech 2003;30(8):637–47.
[14] Kurup PU, Dudani NK. Neural networks for profiling stress history of clays from PCPT data. J Geotech Geoenviron Eng ASCE 2002;128(7):569–79. [15] Basheer IA. Empirical modeling of the compaction curve of cohesive soils. Canad Geotech J 2001;38(1):29–45. [16] Juang CH, Jiang T, Christopher RA. Three dimensional site characterization: neural network approach. Geotechnique 2001; 51(9):799–809. [17] Dayakar P, Rongda Z. Triaxial compression behavior of sand and gravel using artificial neural networks (ANN). Comput Geotech 1999;24(3):207–30. [18] Abu Kiefa MA. General regression neural networks for driven piles in cohesionless soils. J Geotech Geoenviron Eng ASCE 1998;124(12): 1177–85. [19] Jingsheng Shi, Ortigao JAR, Junle Bai. Modular neural networks for predicting settlements during tunneling. J Geotech Geoenviron Eng ASCE 1998;124(5):389–95. [20] Ellis GW, Yao C, Zhao R, Penumadu D. Stress strain modeling of sands using artificial neural networks. J Geotech Eng ASCE 1995;121(5):429–35. [21] Neural networks in geotechnical engineering, Session 68, TRB, 75th annual meeting; 1996. [22] Shahin MA, Maier HR, Jaksa MB. Artificial neural network applications in geotechnical engineering. Austr Geomech 2001; 36(1):49–62. [23] Haykins S. Neural networks – A comprehensive foundation. 2nd ed. Englewood Cliffs (NJ): Prentice Hall International Inc.; 1999. [24] Ranaweera DK, Hubele NF, Papalexopoulos AD. Application of radial basis function for short term load forecasting. IEE Proc Gener Trans 1995;142(1):45–50. [25] Koza JR. Genetic programming: On the programming of computers by means of natural selection. Cambridge (MA): MIT Press; 1992. [26] Neville AM. Properties of concrete. 4th ed.; 1996. [27] ASTM. Standard test method for unconfined compressive strength of cohesive soil. Designation D 2166-91. Annual book of ASTM Standards, vol. 4(8). American Society for Testing and Materials, Philadelphia; 1995. [28] Kezdi A. Stabilised earth roads. Amsterdam, Oxford, New York: Elsevier; 1979. [29] Miura N, Taesiri Y, Koga Y, Nishida K. Practice of improvement of Ariake clay by mixing admixtures. In: Proceedings of the international symposium on shallow sea and low land, Saga University, Saga, Japan; 1988. p. 59–68. [30] Davidson DT, Barnes HF. Improvement of lime stabilization of montmorillonite clay soils with chemical additives. Highway Res Board Bull 1960;262:33–50. [31] Singleton A. Genetic programming with C++. Byte 1994:171–6. [32] Suresh S, Omkar SN, Mani V, Guruprakash TN. Lift coefficient prediction at high angle of attack using recurrent neural networks. Aerospace Sci Technol 2003;7(8):595–602.
Prediction of unconfined compressive strength of soft grounds using computational intelligence techniques: A comparative study B.S. Narendra a, P.V. Sivapullaiah a
a,*
, S. Suresh b, S.N. Omkar
b
Department of Civil Engineering, Indian Institute of Science, Bangalore 560 012, India Department of Aerospace Engineering, Indian Institute of Science, Bangalore, India
b
Received 14 December 2004; received in revised form 2 February 2006; accepted 13 March 2006 Available online 18 May 2006
Abstract Cement stabilization is one of the commonly used techniques to improve the strength of soft ground/clays, generally found along coastal and low land areas. The strength development in cement stabilization technique depends on the soil properties, cement content, curing period and environmental conditions. For optimal and effective utilization of cement, there is a need to develop a mathematical model relating the gain in strength in terms of the variables responsible. The existing empirical model in the literature assumes linear variation of normalized strength with the logarithm of curing period and hence, different empirical models are required for different conditions of the same soil. Also, the accuracy of strength prediction is unsatisfactory. Due to unknown functional relationships and nonlinearity in strength development, in this paper the computational intelligence techniques such as multilayer perceptron (MLP), radial basis function (RBF) and genetic programming (GP) are used to develop a mathematical model. To generate the mathematical model, an experimental study is conducted to obtain the strength of three inland soils namely, red earth (CL), brown earth (CH) and black cotton soil (CH) for different water contents, cement contents and curing periods. In order to generate a generic mathematical model using computational intelligence techniques, two saline soils (Ariake clay-3 and Ariake clay-4) and three inland soils are used. A detailed study of the relative performance of the computational intelligence techniques and the empirical model has been carried out. 2006 Elsevier Ltd. All rights reserved. Keywords: Soft ground; Saline soil; Cement stabilization; Empirical model; Multilayer perceptron; Radial basis function; Genetic programming; Unconfined compressive strength
1. Introduction Due to increased urbanization and industrialization, foundation soils with suitable geotechnical properties are not easily available. There is a growing demand for reclamation of marine and inland soft soils/grounds, which are characterized by high plasticity, higher fraction of fines and void ratio, with low strength and high compressibility [1]. In many cases, the natural water content of soft ground is high and at times even higher than the liquid limit of the soils. For such soils, deep mixing is a popular method for improving the strength [2]. In deep mixing, the cement grout or dry cement is injected into the natural soil at the *
Corresponding author. Tel.: +91 80 22932672; fax: +91 80 23600404. E-mail address: [email protected] (P.V. Sivapullaiah).
0266-352X/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2006.03.006
depth required and a blade is pushed into the ground to mix the soil and cement [3]. The extent of the strength improvement depends on the mineralogy, environmental conditions of the soft ground, curing period and the type and amount of cement used [1]. Porbaha [1], Tan et al. [3], Uddin et al. [4], Yin and Lai [5], Miura et al. [6], Suksun [7] and Horpibulsuk et al. [8,9] have carried out experimental studies to understand the strength improvement in soft ground using cement stabilization techniques. In the stabilization of soft ground, Miura et al. [6], Tan et al. [3] and Horpibulsuk et al. [9] have observed that the strength development depends on the initial clay water content and the cement content. For optimal and effective utilization of cement, a predictive model is required to relate the gain in strength of soil with time in terms of important variables involved. Tan
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
et al. [3] have established an empirical relationship to predict the strength development based on cement content, water content and curing period. The strength developed at specific cement and water contents and curing period is used as the reference compressive strength for a given soil and the strength developed under other conditions for the same soil is normalized using the reference strength. Their empirical model assumes linear variation of normalized strength and the model is different under different conditions of the same soil. The model also assumes that the ratio of the strength developed and the reference strength is constant for any soil if the reference conditions are the same. These assumptions will affect the accuracy of the strength prediction. However, it is mentioned in [3] that a careful statistical study needs to be undertaken to generalize the relationship to improve the performance of the model. One way to improve the performance of the empirical model is to consider the effects of other parameters. Miura et al. [6], Suksun [7] and Horpibulsuk et al. [9] proposed a parameter called the clay water cement ratio (wc/C), which is analogous to the water cement ratio used in concrete technology, which is the ratio of the initial water content of the clay (wc) to the cement content (C). Miura et al. [6] observed that the water cement ratio is the principal parameter for analysis of strength and deformation behavior of cement stabilized soft ground. Using the clay water cement ratio (wc/C), as the principal parameter Suksun [7] and Horpibulsuk et al. [9] developed the normalized empirical model to predict the strength of cement stabilized soft ground. The model requires the reference strength and assumes linear variation of the logarithm of the curing period. However, the empirical model is localized and different equations are considered for different soils depending on the liquidity index. The above limitations of the empirical models suggest the necessity to develop a generic mathematical model for strength prediction. Computational intelligence techniques are capable of approximating the nonlinear input–output relationship for a wide range of applications. Hence in this paper, computational intelligence techniques such as multilayer perceptron (MLP) network, radial basis function (RBF) network and genetic programming (GP) techniques are used to approximate the functional relationship between the various parameters and strength developed. Neural networks (NNs) have been used in geotechnical engineering for estimation of unsaturated shear strength [10], predicting settlement of shallow foundations [11–13], profiling stress history of clays [14], modeling of compaction curve [15], site characterization [16], triaxial compression behavior of sand and gravel [17], predicting the capacity of driven piles in cohesion less soils [18], predicting settlements during tunneling [19] and stress stain modeling of sands [20]. A collection of neural network applications in geotechnical engineering are presented in [21,22]. NNs gained prominence in the 1980s due to their ability to learn from past data, and their ability to generalize
197
knowledge obtained during learning. A variety of neural network architectures and associated learning algorithms are available in the literature for various types of applications. Of these, two neural network models namely, MLPs [23] and RBFs [24] are considered. MLPs are feed forward networks with a layered structure. In general, MLPs have an input layer, one or more hidden layers, and one output layer. Each of the nodes in a given layer are connected to the other nodes in the next layer. There are no connections within layers. In most of the applications a feed forward network with a single hidden layer is used with sigmoidal activation functions for the units. The back propagation (BP) algorithm is widely used for training the MLP networks. Although MLPs can be used with any number of hidden layers, it has been shown that only one hidden layer suffices to approximate any function with discontinuities to an arbitrary degree of precision, provided the activation function of the hidden units are nonlinear (the universal approximation theorem) [23]. The RBF network model also consist of three layers; input, hidden and output layers. The nodes within each layer are fully connected to the previous layer. The input variables are each assigned to nodes in the input layer, and pass directly to the hidden layer without weights. The hidden nodes or units contain the radial basis functions. The RBF is similar to the Gaussian density function, which is defined by parameters namely, a ‘center’ and a ‘width’ [24]. The Gaussian function gives the highest output when the incoming variables are closest to the center position and the output decreases monotonically as the distance from the center increases. The ‘width’ of RBFs controls the rate of decrease in the output. Even though RBFs are well suited for function approximation, for large sets of data the problem of ‘curse of dimensionality’ becomes predominant [24], i.e., the increase in number of features reduces the approximating ability of the RBF network. In NNs, the basic limitation is that the optimal configuration is not known a priori. Also, often it is not possible to express the relationship between the different parameters governing a phenomenon because the knowledge represented in the weights is generally opaque i.e., the weights do not reveal any information about the input–output relationship. One of the recent developments in computational intelligence techniques is genetic programming (GP) which has the ability to approximate the relationship between the inputs and outputs and express the relationship mathematically [25]. The generated mathematical expression may also provide an insight into the understanding of the problem. Hence in this paper, the MLP, RBF and GP techniques are used to develop a generic model for strength prediction. In order to generate such a model, three inland soils for which experimental studies are carried out, and two saline soils available in the literature [7] are considered. For these soils, empirical models are developed using the concept presented by Horpibulsuk et al. [9] and the performance is compared with the generic models developed using the computational intelligence technique discussed above.
198
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
2. Experimental procedure This section contains description of the materials used in the investigation and the experiments carried out for the determination of the unconfined compressive strength of cement stabilized brown earth, black cotton soil and red earth under different water and cement contents and curing periods. 2.1. Materials used Cement. Portland cement (43 Grade) was used as the stabilizer in the investigation. Clay water cement ratios of 5, 10, 15 and 20 were considered in the present study. The cement content required to prepare the specimen depends on the clay water cement ratio and initial clay water content. Soils. Three Indian soils namely, red earth, classified as CL (lean clay), brown earth and black cotton soil, both classified as CH (fat clay) are used in the study. Red earth is obtained from Indian Institute of Science campus, Bangalore, Karnataka state, India, Brown earth is obtained from a construction site on Airport road in Bangalore, Karnataka state, India and Black cotton soil is obtained from Bagalkot, Karnataka state, India. All the soils are collected by open excavation from a depth of 1 m from the natural ground. The geotechnical properties of the soils are given in Table 1. 2.2. Sample preparation The samples are prepared using split mould and are of size of 38 mm diameter and 76 mm length. Initially, the soil and de-ionized water are mixed to form a soil paste. Initial water content equal to liquid limit, 1.5 times the liquid limit and 2 times liquid limit are considered for specimen generation. For hydration process, the cement consumes water equal to 40% of its mass [26]. Hence, the cement paste is
2.3. Experimental results The specimens are cured for curing period (D) of 7, 14, 28 and 56 days. The unconfined compressive strength of specimens is determined using ASTM D2166-91 method [27]. Table 2 presents the results of unconfined compressive strength (kPa) of brown earth stabilized with cement for different initial clay water content (wc), clay water cement ratio (wc/C) and curing period. Similarly, the results for Table 2 Unconfined compressive strength (kPa) of brown earth stabilized with cement wc/C
Unconfined compressive strength (kPa) 7 Days
14 Days
28 Days
Water content, wc = liquid limit water content (%) 5 1182 1578 2552 10 569 585 656 15 150 225 232 20 15 49 69
56 Days 2699 736 284 98
Water content, wc = 1.5 · liquid limit water content (%) 5 438 644 846 10 279 419 497 15 161 171 193 20 58 62 74
932 595 238 93
Water content, wc = 2 · liquid 5 145 10 127 15 99 20 66
390 374 203 95
limit water content (%) 231 351 164 277 135 184 69 94
Table 3 Unconfined compressive strength (kPa) of black cotton soil stabilized with cement
Table 1 Geotechnical properties of soils Properties
prepared separately by mixing cement with de-ionized water equal to 0.4 times the mass of cement. Finally, both the soil paste and the cement paste are mixed thoroughly and put into the split mould. The specimens are cured in desiccators at 100% humidity, ensuring the availability of water for cementitious reactions.
Red earth
Brown soil
2.64
2.63
Specific gravity Atterberg limits (%) Liquid limit Plastic limit Plasticity index Shrinkage limit Free swell index
38 15 23 16 25
60 23 37 13 36
Standard proctor test Maximum dry density (kN/m3) Optimum moisture content (%)
16.4 18
14.3 28
Grain size distribution (%) Fine sand Silt Clay Classification
40 28 32 CL
20 34 46 CH
Black cotton soil 2.76 97 35 62 9 138 11.5 37 2 37 61 CH
wc/C
Unconfined compressive strength (kPa) 7 Days
14 Days
28 Days
56 Days
Water content, wc = liquid limit water content (%) 5 1842 2255 2636 10 515 630 781 15 120 207 250 20 82 117 155
3702 981 330 255
Water content, wc = 1.5 · liquid limit water content (%) 5 276 425 631 10 92 235 390 15 73 127 154 20 58 69 91
1136 500 223 119
Water content, wc = 2 · liquid limit water content (%) 5 122 284 343 10 79 132 212 15 77 99 135 20 51 75 89
548 334 166 105
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208 Table 4 Unconfined compressive strength (kPa) of red earth soil stabilized with cement wc/C
14 Days
28 Days
Water content, wc = liquid limit water content (%) 5 476 683 857 10 141 166 196
56 Days 1572 311
Water content, wc = 1.5 · liquid limit water content (%) 5 116 166 349 10 112 147 152
610 284
Water content, wc = 2 · liquid 5 88 10 76 15 61 20 34
552 310 196 95
limit water content (%) 158 323 112 161 79 128 60 85
black cotton soil and red earth are given in Tables 3 and 4, respectively. 3. Empirical formula The empirical model for brown earth, black cotton soil and red earth are developed using the concept presented in [9] are discussed below. 3.1. Development of empirical formula Suksun [7] and Horpibulsuk et al. [9] used Abram’s law as the basis for the model development. Abram’s law is extensively used in concrete technology and can be expressed as: A ; ð1Þ Bðw=CÞ where Q is the compressive strength of the concrete at a specific age, w/C is the water cement ratio by weight, and A and B are empirical constants. Using Abram’s law, a similar parameter, called the clay water cement ratio (wc/C) is derived and Abram’s law is modified as: Q¼
Q¼
A . Bðwc =CÞ
Qðwc =CÞ2 Qðwc =CÞ1
Unconfined compressive strength (kPa) 7 Days
Qðwc =CÞ1
ð2Þ
The modified power fit is used for the plot of the clay water cement ratio versus curing periods for different liquidity indices (LI). Using these plots, the values of A and B are calculated. The value of A varies over wide range and it depends on the type of clay, the liquidity index and curing period. However, the variation of the parameter B is only between 1.22 and 1.24, irrespective of the type of clay and the curing period. Hence, the value of parameter B is selected as 1.24. The parameter A has a specific value for a clay at a specified curing period (D) and liquidity index (LI) but does not change with clay water cement ratio. The parameter A is eliminated by taking the ratios of strength developed at different clay water cement ratios as follows:
Qðwc =CÞ2
¼
A=Bðwc =CÞ1 ¼ Bfðwc =CÞ1 ðwc =CÞ2 g ; A=Bðwc =CÞ2
¼ 1:24fðwc =CÞ1 ðwc =CÞ2 g ;
199
ð3Þ ð4Þ
where Qðwc =CÞ1 is the unconfined compressive strength at wc/C1 and Qðwc =CÞ2 is the unconfined compressive strength at wc/C2. Using the strength developed at the curing period (D = 28) as a reference point, the normalized strength for different types of soft ground, stabilized at different liquidity indices is calculated. The emprical relationship between the normalized strength and curing period for different liquidity indices is expressed as 0:038 þ 0:281 lnðDÞ if LI ¼ 1:0–2:5; QD ¼ ð5Þ Q28 0:126 þ 0:342 lnðDÞ if LI > 2:5; where QD is the strength to be estimated at D days of curing, and Q28 is the strength at the reference curing period (D = 28). Using the concept explained in [7], the empirical model is developed using Eqs. (4) and (5) 8 1:24fðwc =CÞ2 ðwc =CÞ1 g ð0:038 þ 0:281 lnðDÞÞ > > > < Qðwc =CÞ1;D if LI ¼ 1:0–2:5; ¼ Qðwc =CÞ2;28 > 1:24fðwc =CÞ2 ðwc =CÞ1 g ð0:126 þ 0:342 lnðDÞÞ > > : if LI > 2:5; ð6Þ where D is the curing period in days, Qðwc =CÞ1;D is the unconfined compressive strength at wc/C1 for the curing period = D, and Qðwc =CÞ2;28 is the unconfined compressive strength at wc/C2 for the reference curing period = 28 days. 3.2. Limitations of empirical formula Using the empirical model given in Eq. (6) and the reference strengths highlighted in Table 5, the unconfined compressive strength of black cotton soil has been predicted. The predicted and actual strength for black cotton soil for various water cement ratios and curing periods are listed in Table 5. The root mean square error deviation from actual strength is 226 kPa. The higher error may be due to the fact that certain parameters, such as sodium salts and pH, which play an important role in cement stabilization [28,29], are not considered. Also, the empirical model is based on a best fit of the normalized strength for different soils. However, the reference used for predicting the strength in Table 5 is different for each soil. It can be observed in Table 5 that two separate references have been used in the prediction of strength for liquidity indices of 1.0 and 1.8. Even though the liquidity index is less than 2.5, the use of a single reference strength in the equation is found to result in large errors. For example consider the 16th data sample in Table 5. By taking 781 kPa as the reference, the prediction is 106 kPa, whereas a reference of
200
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
Table 5 Actual and predicted unconfined compressive strength using empirical model for black cotton soil Serial number
Curing period (days)
Clay water cement ratio
Cement content (%)
Liquid limit (%)
Liquidity index
Water content (%)
Actual UCS (kPa)
Predicted UCS (kPa)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
7 7 7 7 14 14 14 14 28 28 28 28 56 56 56 56 7 7 7 7 14 14 14 14 28 28 28 28 56 56 56 56 7 7 7 7 14 14 14 14 28 28 28 28 56 56 56 56
5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20 5 10 15 20
19.40 9.70 6.47 4.85 19.40 9.70 6.47 4.85 19.40 9.70 6.47 4.85 19.40 9.70 6.47 4.85 29.10 14.55 9.70 7.28 29.10 14.55 9.70 7.28 29.10 14.55 9.70 7.28 29.10 14.55 9.70 7.28 38.80 19.40 12.93 9.70 38.80 19.40 12.93 9.70 38.80 19.40 12.93 9.70 38.80 19.40 12.93 9.70
97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 1.78 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55 2.55
97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 97 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 145.5 194 194 194 194 194 194 194 194 194 194 194 194 194 194 194 194
1842 515 120 82 2255 630 207 117 2636 781 (Ref) 250 155 3702 981 330 255 276 92 73 58 425 235 127 69 631 390 (Ref) 154 91 1136 500 223 119 122 79 77 51 284 132 99 75 343 212 (Ref) 135 89 548 334 166 105
1339 457 156 53 1785 609 208 71 2231 761 260 89 2677 913 311 106 669 228 78 27 891 304 104 35 1114 380 130 44 1337 456 156 53 335 114 39 13 483 165 56 19 630 215 73 25 777 265 90 31
390 kPa results in 53 kPa. This implies that the empirical model is not really generalized and more fine tuning the equation is required to improve the performance of each and every soil. 3.3. Parameter selection In order to obtain a generic model for unconfined compressive strength prediction, saline soils such as Ariake clay-3 and Ariake clay-4 presented in [7] are also considered. From the limitation of the empirical model, one can observe that the other parameters such as liquid limit
(LL), liquidity index (LI), pH and sodium ion concentration (Na+) play an important role in strength development. The liquid limit is one of the very important soil properties, which has been extensively correlated with the behaviour of soils. Liquidity index can indicate the interparticle cementation of the soil particles. The pH will affect the cementitious reaction to greater extent. It is well known that if the pH is greater than 12, the amorphous silica present in the clay gets dissolved and reacts with calcium hydroxide to form calcium–silicate–hydrate gel, which imparts strength to the stabilized soil [30]. Since, the salinity also affects the strength development the
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
sodium ion concentration present in soil is used to indicate the salinity in the soil. Hence, the liquid limit (LL), liquidity index (LI), water content (wc), cement content (C), curing period, (D), pH and sodium ion concentration (Na+) are used to develop the model. 4. Computational intelligence techniques In this section, the two popular neural network architectures, multilayer perceptron, radial basis function network, and genetic programming are described briefly. 4.1. Multilayer perceptron network Multilayer perception provides a general, nonlinear parameterized mapping between a set of inputs and a set of outputs. A MLP network with three layers and sigmoidal activation functions can approximate any smooth mapping, if sufficient number of neurons is included in the hidden layer [23]. The general architecture of MLP network is shown in Fig. 1. The MLP network is trained using supervised back propagation learning algorithm [23]. The functional relationship in the strength development is established through MLP network. The learning algorithm is designed to minimize a cost function J of the form
J¼
1XX i ðy k ^y ik Þ2 ; 2 i k
201
ð7Þ
where y ik and ^y ik are the desired output and the actual output vectors of the kth output of ith sample, respectively. Let the input and output vectors are denoted by u and ^y , respectively. The weight connection between the ith neuron in the input layer and the jth neuron in the hidden layer is represented as w1ij . Similarly, let the weight connection between the ith neuron in the hidden layer and the jth neuron in the output layer be w2ij . The output of the ith neuron in the hidden layer is described mathematically as ! X 1 hi ¼ f ðneti Þ ¼ f wij uj ð8Þ j
and the output of the ith neuron in the output layer is described mathematically as ! X 2 ^y i ¼ f ðneti Þ ¼ f ð9Þ wij hj ; j
where f(Æ) is the sigmoidal activation function. The interconnection weights are adjusted such that the cost function J can be minimized, and are given by
Fig. 1. Architecture of multilayer perceptron.
202
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
w1ij ¼ w1ij þ gd1i uj ; w2ij
¼
w2ij
þ
gd2i hj ;
ð10Þ
where dli is error term of jth neuron in lth layer and g is learning rate. The error terms are defined below: d2i ¼ ðy i ^y i Þf 0 ðnet2i Þ; ! X w2ij d2i f 0 ðnet1j Þ; d2j ¼
ð11Þ
i
where f 0 (Æ) is the first derivative with respective to their arguments. The training process is terminated when the cost function is sufficiently small. The details of the learning algorithm are given in [23]. 4.2. Radial basis function network The radial basis function network consists of three layers: an input layer, a kernel layer and an output layer as shown in Fig. 2. The neurons in input layer are assigned to kernel neurons in the hidden layer without any modification, i.e., the weight connection is unity. The kernel layer nonlinearly transforms the inputs to a higher dimension. The activation function in the kernel neurons is similar to the Gaussian function, which is defined by a ‘centre’ position and a ‘width’ parameter. The Gaussian function gives the highest output when the incoming variables are closest to the centre position and decreases monotonically as the distance from the centre increases. The width of the Gaussian function controls the rate of decreases. Let x be the input vector of dimension n. The output (hi) of the ith kernel neuron is given by , ! n X 2 hi ¼ exp ½xk xik r2i i ¼ 1; 2; . . . ; H ; ð12Þ k¼1
where H is the number of kernel neurons, xik is centre of the ith kernel neuron for the kth input variable and ri is the width of the ith kernel neuron. The output value of the ith output neuron is equal to the summation of the weighted outputs of the kernel neurons
Fig. 2. Architecture of radial basis function network.
and bias term of the output neuron, and is described mathematically as ^y i ¼
H X
wij hj þ wio ;
i ¼ 1; 2; . . . ; O;
ð13Þ
j¼1
where O is the number of output neurons and wio is bias value of the ith output neuron. The parameters of the RBF network are determined using a three step training scheme [24]. First, the kernel neurons centers are determined using the k-means clustering algorithm. Then, the width parameters are determined by a nearest-neighbor method. Finally, weights connecting kernel neurons and output neurons are calculated using multiple linear regression techniques. The details of the learning algorithm are given in [24]. 4.3. Genetic programming Genetic programming (GP) is a branch of genetic algorithms. The main difference between genetic programming and genetic algorithms is the representation of the solution. Genetic programming creates computer programs in the LISP or SCHEME computer languages as the solution, where as genetic algorithms create a string of numbers that represent the solution. GP has gained prominence in the 1990s and it has an ability to understand the underlying relationship between the input and outputs and express them mathematically. GP is provided with evaluation data, a set of primitives and fitness functions. The evaluation data describe the specific problem in terms of the desired inputs and outputs. They are used to generate the best computer program to describe the relationship between the input and output very well [25]. The following four steps are used to generate a mathematical relationship between the input/output data using: 1. Generate an initial population of random compositions of the functions and terminals of the problem (computer programs). 2. Execute each program in the population and assign it a fitness value according to how well it solves the problem. 3. Create a new population of computer programs. Copy the best existing programs Create new computer programs by crossover and mutation. 4. The best computer program that appeared in any generation, the best-so-far solution, is designated as the result of genetic programming. In GP, the best idea is to maximize the user defined fitness function. Thus, at the end of each GP run, a current population (computer program) with good fitness, and best fit individual is selected as solution to the given problem. The following subsections present a brief description on various components of GP.
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
4.3.1. Crossover operation The crossover operation is used to create new offspring computer programs from two parental computer programs selected based on the fitness value. The parental computer programs are typically of different sizes and shapes. The offspring computer programs are composed of subexpressions (subtress subprograms, subroutines and building blocks) from their parents, and typically of different sizes and shapes from their parents. The creation of offspring from the crossover operation is accomplished by swapping random fragments of the parents. Fig. 3 represents a typical crossover operation. 4.3.2. Mutation operation Mutation operation is shown in Fig. 4. Mutation operation protects the model against premature convergence and improves the non-local properties of the search. Mutation is implemented by randomly removing a subtree at a selected points and replacing it with randomly generated subtree. 4.3.3. Tournament selection Tournament selection is a process of selecting the computer program for genetic operations. Typically, the genetic programming implementation uses a model in which certain randomly selected computer programs in a subgroup compete and the fittest is selected. The tournament size defines the number of randomly selected computer programs for selection process.
203
Fig. 4. Mutation operation.
4.3.4. Fitness function GP is guided by fitness function to search for the most efficient computer program to solve a given problem. The negative of mean square error is selected as a simple measure of fitness in function approximation problem Fitness ¼
N 1 X 2 ðY i Y^ i Þ ; N i¼1
ð14Þ
where Y is the actual value and Y^ is the predicted value. 4.3.5. Functions and terminals The terminal and function sets are also important components of genetic programming. The terminal and function sets are the alphabets of the programs to be made. The terminal set consists of the variables and constants
Fig. 3. Crossover operation.
204
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
Table 6 List of parameters used in GPQUICK software Parameters
Values
Population size GP generations Crossover weight Mutation weight Mutation rate Tournament size Tree length
5000 1,000,000 70 20 60 6 25
of the programs and the function set consists of the functions such as addition (ADD), subtraction (SUB), multiplication (MUL), division (DIV), square root (SQ) and the natural logarithm (LN). 4.3.6. Termination criterion Koza [25] has shown that in GP the evolution is a neverending process, and hence a termination criterion is needed. The termination criterion for GP is generally based on the problem or is limited by the number of generations. In GP, a user-defined fitness function has to be maximized for his/her application. Thus at the end of a GP run, a current population of individuals (computer programs), and also the fittest individual (computer program) that appeared during the run is obtained. The fittest individual that has evolved for the given problem is its solution or desired mathematical model. GPQUICK software [31] has been used for this function approximation problem. In the GP runs, a maximum number of generations (1,000,000) or MSE is equal to 0.002 is used as the criterion to stop the program. The choice of parameters for GPQUICK software used in the runs is shown in Table 6. The parameters for genetic programming are selected using trial and error procedure. The parameter selection will affect convergence of GP [25].
computational intelligence techniques. Training error (TE) is the root mean sum of squared residuals (error) in the training data and generalization/testing error (GE) is the root mean squared error in the (validating) testing data vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N1 u1 X ^ i Þ2 ; TE ¼ t ð15Þ ðQ Q N 1 i¼1 i vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N2 u1 X ^ i Þ2 ; ð16Þ ðQ Q GE ¼ t N 2 i¼1 i ^ i the estiwhere Qi is the actual value of strength, and Q mated output of the computational intelligence technique for the ith input. N1 is the number of data points used in the training set and N2 is the number of data points used in the testing set. The total error (TOE) takes into account both the training and generalization error, i.e., sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N 1 TE2 þ N 2 GE2 TOE ¼ . ð17Þ ðN 1 þ N 2 Þ Since the TE, GE and TOE will not reveal substantial statistical information about the difference between the actual and the predicted value, we use statistical estimate like R2 as one of the performance measure. R2 is the square of correlation between the actual value and the predicted value or it is an indicator of how much the variation in the data is explained by the model [32]. It is calculated using the following equation: R2 ¼ 1
SS T ; SS R
ð18Þ
where SS R ¼
N X
2
ðQðtÞ QÞ ;
ð19Þ
N X 2 ^ ðQðtÞ QðtÞÞ ;
ð20Þ
t¼1
4.3.7. Function approximation using GP Koza [25] has applied GP for a function approximation problem (regression problem). In function approximation problem a single GP expression is evolved for each output. This GP expression evolved for each output is called as genetic programming approximator expression (GPAE). This GPAE approximates the functional relationship between the input–output data set. 5. Computational intelligence based unconfined compressive strength prediction In this section the procedures to develop a generic model to predict the unconfined compressive strength using MLP, RBF and GP are presented. The performances of these models with the empirical model are also compared. 5.1. Performance measures The training error, generalization error and total error are measured to study the qualitative performance of the
SS T ¼
t¼1
^ is the predicted where Q is the mean value of strength and Q strength. R2 can take any value between 0 and 1, with a value closer to 1 indicating a better fit [32]. For example, an R2 value of 0.85 means that the approximation explains 85% of the total variation in the data about the average. 5.2. Selection of training data To predict the unconfined compressive strength (kPa), the following eight input parameters are selected for the computational intelligence techniques, curing period (D), clay water cement ratio (wc/C), cement content (C) in percentage, liquid limit (LL) in percentage, liquidity index (LI), water content (wc) in percentage, pH and sodium ion concentration (Na+) in gram/litre. The input soil characteristic and the output strength are stored in a file for different soils. The 186 data are recorded for various conditions and these data are used for the model develop-
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
ment using computational intelligence techniques. Out of 186 data, 154 data are selected randomly for training and 32 data for testing the generalization capability of the computational intelligence techniques. 5.3. Model generation using neural networks The details of the training process adopted for neural networks and optimal network selection procedure are described below. 5.3.1. Multilayer perceptron The input to the network is eight soil characteristics and one bias. The output of the network is the strength developed in the soil. The network architecture is represented as N n1 ;n2 ;n3 , where n1 is the number of input nodes, n2 is the number of hidden nodes and n3 is the number of output nodes. A learning rate, g, of 0.0003 was used for training the MLP network. The network is trained for 75,000 epochs or training error less than 0.002. The optimal network configuration of MLP is N9,22,1. The optimal configuration of network to approximate the input–output relationship is found using the steps explained in [32]. The optimal configuration is tested with the same training data and the strength predicted and the actual strength are shown in Fig. 5. A few training points are selected for this purpose. From the figure, it can be observed that the predicted strength closely matches with the actual values at various sample points. 5.3.2. Radial basis function network Similar to the MLP network, the input to the radial basis function network are eight soil parameters and output of the network is strength. The network is represented as N n1 ;n2 ;n3 ;p , where n2 is number of radial basis function neu-
205
rons in the network and p is the overlap parameter. All eight inputs are used to calculate the centre and width of the RBF network using the k-means clustering algorithm and linear weight matrix between the RBF neurons and output neuron is calculated using regression procedure. The optimal RBF configuration is N8,97,1,95. Here also, the optimal configuration of network to approximate the input–output relationship is found using the steps explained in [32]. The trained RBF network is tested with same training data. The predicted and the actual strength are plotted in Fig. 5. From the above figure, it can be observed that the RBF network approximates the input–output relationship very well. 5.4. Model generation using genetic programming The same training data are used to generate the genetic programming expressions. The best GP program is evaluated until it reached fitness measure (0.002) or number of population is equal to 1,000,000. Four runs of GP with the training set were used and obtained the best computer program evolved at each run. At the end of each run, the best computer program evolved is in the form of LISP expression and this expression can be easily converted into a mathematical expression. Initially for the first run, the tree length of LISP expression is restricted to 1000. The resulting LISP expression is able to predict the strength development with considerable accuracy, but reducing the LISP expression to obtain mathematical equations is very difficult. The performance of the LISP expression is given in Table 7. Even though the performance of the above LISP expression is better, it is very difficult to study the effect of various parameters in strength development. Hence for the remaining three GP runs, the tree length was restricted to 25 such that the resulting LISP expression can be simplified to obtain a mathematical relationship.
4000 Actual GP RBF MLP
Unconfined compressive strength (kPa)
3500
3000
2500
2000
1500
1000
500
0 0
5
10
15
20
25
30
35
40
Data points Fig. 5. Comparison of trained data for different computational intelligence techniques with actual data.
206
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
Table 7 Performance measures for empirical and computational intelligence techniques based models Maximum error
Empirical model Multilayer perceptron Radial basis function GP (tree length = 25) GP (tree length = 1000)
Training
Testing
Total
Training
Testing
Total
Training
Testing
Total
– 221 413 505 251
– 178 420 572 324
1025 221 420 572 324
– 50 136 151 72
– 71 133 195 97
226 54 135 159 77
– 0.9948 0.9616 0.9527 0.9893
– 0.9914 0.9696 0.9345 0.9837
0.9080 0.9941 0.9632 0.9490 0.9881
Out of the three runs with restricted tree length, the best LISP expression obtained is given below: (ADD (MUL 35 B7) (SUB (DIV (MUL (MUL B1 B3) (DIV B4 B5)) (ADD (DIV (MUL B2 (DIV B2 B1)) B3) (ADD (DIV (MUL B7 (DIV (DIV (DIV B8 B5) B5) B5)) B3) (DIV B5 (DIV (MUL B7 (DIV 26 B3)) B3))))) (MUL (ADD (DIV 19 B8) (ADD (SUB (MUL B8 B5) (MUL B8 B3)) (DIV B3 (SUB 19 (MUL B8 B3))))) (MUL (MUL B8 B3) (DIV B2 B7))))) where B1 is the curing period (D) in days, B2 the clay water cement ratio (wc/C), B3 the cement content (C) in percentage, B4 the liquidity index (LI), B5 the water content (wc) in percentage, B6 the pH, and B7 is the sodium ion concentration (Na+) in gram/litre. The corresponding GPAE is after simplification is ðCÞðLLÞLnðDÞ i ðpHÞðNaþ Þ ðCÞ2 ðLIÞ þ ðCÞðLIÞ3 þ 26ð LI pHÞ ðwc =CÞðCÞðNaþ Þ 19 þ ðLIÞðNaþ Þ ðCÞðNaþ Þ ðpHÞ Naþ C þ . ð21Þ ð19 ðCÞðNaþ ÞÞ
Q ¼ 35ðpHÞ þ
R2 estimate
Root mean square error
h
ðwc =CÞ2 ðCÞLnðDÞ
The above expression is complex and is difficult to obtain any physical insight about the behavior of the soils in strength development. On the other hand, the mathemati-
cal expression is able to predict the strength development accurately over wide range of soils. The performance of the best LISP expression is given in Table 7. From Table 7, it is evident that the prediction based GP model developed using unrestricted tree length has better performance than the restricted tree length. The predicted strength and the actual strength are shown in Fig. 5. From the figure, it can be observed that the GP is able to approximate the underlying the relationship very well. 5.5. Performance of computational intelligence techniques models The quantitative performance of computational intelligence techniques for training data set is given in Table 7. From the table it can be observed that the TE for MLP is less than the RBF and GP. The same is reflected in the R2 estimate. 6. Prediction of unconfined compressive strength The MLP, RBF, GP and empirical models are used to predict the strength development in soft ground using cement stabilization technique and the performance of these models is compared. The predicted values of strength
3500 Actual GP RBF MLP Empirical
Unconfined compressive strength (kPa)
3000
2500
2000
1500
1000
500
0 0
5
10
15
20
25
30
Data points Fig. 6. Actual strength and predicted strength using empirical model and computational intelligence techniques.
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
for various experimental data are shown in Fig. 6. From the figure, one can observe that it is not possible to assess the accuracy of different predictive methods. From the performance measures listed in Table 7, it can be observed that for the unconfined compressive strength predicted using the empirical model, the maximum error equal to 1025 and the RMSE equal to 226. The errors are the higher when compared to the evolutionary computational intelligence techniques. Also the R2 estimate = 0.9080 is the lowest for empirical model. Hence, the prediction using computational intelligence techniques is more accurate than the conventional empirical model. Table 8 gives the details of some of the testing data points used for the comparative study of the efficiency of different techniques. In the computational intelligence techniques best prediction is obtained using MLP, with maximum error = 221, TOE = 54 and total R2 estimate = 0.9941. Even though the MLP and RBF perform better than the GP, the knowledge gained by the network weights is opa-
207
que. Hence, one cannot assess the role played by the different input parameters in the strength development. The mathematical expression generated using the GP has a potential to assess different parameters in strength development. The method has the advantage of continuous up-gradation when the new soil data are available, which improve the reliability of the model. 7. Conclusions A generic mathematical model is developed to predict the unconfined compressive strength of cement treated soils using multilayer perceptron, radial basis function network and genetic programming techniques for two saline soils in the literature and three inland soils. The quantitative performance of the computational intelligence techniques is measured and compared with the existing empirical model. The performance measures show that the MLP network has better performance followed by GP (tree length =
Table 8 Details of parameters of different soils used as testing data Liquidity index
Water content (%)
pH
Na+ (g/l)
UCS (kPa)
60 60 60 60 60 60 60
1 1 1 1.8 1.8 2.6 2.6
60 60 60 90 90 120 120
8.6 8.6 8.6 8.6 8.6 8.6 8.6
0.034 0.034 0.034 0.034 0.034 0.034 0.034
585 49 2552 171 846 69 277
9.7 19.4 4.9 29.1 7.3 19.4 12.9 38.8 12.9
97 97 97 97 97 97 97 97 97
1 1 1 1.8 1.8 2.6 2.6 2.6 2.6
97 97 97 145.5 145.5 194 194 194 194
8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8 8.8
0.159 0.159 0.159 0.159 0.159 0.159 0.159 0.159 0.159
630 2636 155 631 91 132 99 343 135
10 5 10 10 20 5
3.8 7.7 5.8 7.7 3.8 15.3
38 38 38 38 38 38
1 1 1.8 2.7 2.7 2.7
38 38 57 76 76 76
8.3 8.3 8.3 8.3 8.3 8.3
0.013 0.013 0.013 0.013 0.013 0.013
166 857 152 112 60 323
Ariake clay 3 28 7 28 7 14 28 14 14
7.1 7.1 8.7 8.7 10.7 10.7 12.3 9.3
15 15 15 15 15 15 15 20
106 106 106 106 106 106 106 106
1 1 1.5 1.5 2 2 2.5 2.5
106 106 130 130 160 160 185 185
8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5
1.6 1.6 1.6 1.6 1.6 1.6 1.6 1.6
1300 2737 1000 1867 589 857 714 1228
Ariake clay 4 7 7
12.2 10.1
10 15
122 122
1 1.5
122 152
8.8 8.8
1 1
661 1286
Curing period (days)
Clay water cement ratio
Cement content (%)
Brown earth 14 14 28 14 28 14 28
10 20 5 15 5 20 10
5.9 3.0 11.8 5.9 17.7 5.9 11.8
Black cotton soil 14 28 28 28 56 14 28 28 14
10 5 20 5 20 10 15 5 15
Red earth 28 28 14 14 28 56
Liquid limit (%)
208
B.S. Narendra et al. / Computers and Geotechnics 33 (2006) 196–208
1000), RBF, GP (tree length = 25) and empirical model. Even though the MLP and RBF perform very well, it is not able to express the relationship between the different parameters governing the strength development because the knowledge gained by the weights is often opaque. But the GP has the ability to learn the underlying input–output relationship and express them mathematically (GPAE). One can use these GP expressions to understand the influence of input parameters for strength development. References [1] Porbaha A. State of the art in deep mixing technology: Part I. Basic concepts and overview. Ground Improvement 1998;2(2):81–92. [2] Horpibulsuk S, Miura N, Koga H, Nagaraj TS. Analysis and strength development in deep mixing: a field study. Ground Improvement 2004;8(2):59–68. [3] Tan TS, Goh TL, Yong KY. Properties of singapore marine clays improved by cement mixing. Geotech Testing J ASTM 2002;25(4): 422–33. [4] Uddin K, Balasubramaniam AS, Bergado DT. Engineering behavior of cement treated Bangkok soft clay. Geotech Eng J 1997;28(1):89–119. [5] Yin JH, Lai CK. Strength and stiffness of Hong Kong marine deposit mixed with cement. Geotech Eng J 1998;29(1):29–44. [6] Miura N, Horpibussuk S, Nagaraj TS. Engineering behavior of cement stabilized clay at high water content. Soils Foundations 2001;41(5):33–45. [7] Suksun Horpibulsuk. Analysis and assessment of engineering behavior of cement stabilised clays. PhD Thesis, Saga University, Japan; 2001. [8] Horpibulsuk S, Bergado DT, Lorenzo GA. Compressibility of cementadmixed clays at high water content. Geotechnique 2004;54(2):151–4. [9] Horpibulsuk S, Miura N, Nagaraj TS. Assessment of strength development in cement admixed high water content clays with Abrams law as basis. Geotechnique 2003;53(2):439–44. [10] Lee SJ, Lee SR, Kim YS. An approach to estimate unsaturated shear strength using artificial neural network and hyperbolic formulation. Comput Geotech 2003;30(6):489–503. [11] Shahin MA, Maier HR, Jaksa MB. Predicting settlement of shallow foundations using neural networks. J Geotech Geoenviron Eng ASCE 2002;128(9):785–93. [12] Shahin MA, Maier HR, Jaksa MB. Data division for developing neural networks applied to geotechnical engineering. J Comput Civil Eng ASCE 2004;18(2):105–14. [13] Shahin MA, Maier HR, Jaksa MB. Settlement prediction of shallow foundations on granular soils using B-spline neurofuzzy models. Comput Geotech 2003;30(8):637–47.
[14] Kurup PU, Dudani NK. Neural networks for profiling stress history of clays from PCPT data. J Geotech Geoenviron Eng ASCE 2002;128(7):569–79. [15] Basheer IA. Empirical modeling of the compaction curve of cohesive soils. Canad Geotech J 2001;38(1):29–45. [16] Juang CH, Jiang T, Christopher RA. Three dimensional site characterization: neural network approach. Geotechnique 2001; 51(9):799–809. [17] Dayakar P, Rongda Z. Triaxial compression behavior of sand and gravel using artificial neural networks (ANN). Comput Geotech 1999;24(3):207–30. [18] Abu Kiefa MA. General regression neural networks for driven piles in cohesionless soils. J Geotech Geoenviron Eng ASCE 1998;124(12): 1177–85. [19] Jingsheng Shi, Ortigao JAR, Junle Bai. Modular neural networks for predicting settlements during tunneling. J Geotech Geoenviron Eng ASCE 1998;124(5):389–95. [20] Ellis GW, Yao C, Zhao R, Penumadu D. Stress strain modeling of sands using artificial neural networks. J Geotech Eng ASCE 1995;121(5):429–35. [21] Neural networks in geotechnical engineering, Session 68, TRB, 75th annual meeting; 1996. [22] Shahin MA, Maier HR, Jaksa MB. Artificial neural network applications in geotechnical engineering. Austr Geomech 2001; 36(1):49–62. [23] Haykins S. Neural networks – A comprehensive foundation. 2nd ed. Englewood Cliffs (NJ): Prentice Hall International Inc.; 1999. [24] Ranaweera DK, Hubele NF, Papalexopoulos AD. Application of radial basis function for short term load forecasting. IEE Proc Gener Trans 1995;142(1):45–50. [25] Koza JR. Genetic programming: On the programming of computers by means of natural selection. Cambridge (MA): MIT Press; 1992. [26] Neville AM. Properties of concrete. 4th ed.; 1996. [27] ASTM. Standard test method for unconfined compressive strength of cohesive soil. Designation D 2166-91. Annual book of ASTM Standards, vol. 4(8). American Society for Testing and Materials, Philadelphia; 1995. [28] Kezdi A. Stabilised earth roads. Amsterdam, Oxford, New York: Elsevier; 1979. [29] Miura N, Taesiri Y, Koga Y, Nishida K. Practice of improvement of Ariake clay by mixing admixtures. In: Proceedings of the international symposium on shallow sea and low land, Saga University, Saga, Japan; 1988. p. 59–68. [30] Davidson DT, Barnes HF. Improvement of lime stabilization of montmorillonite clay soils with chemical additives. Highway Res Board Bull 1960;262:33–50. [31] Singleton A. Genetic programming with C++. Byte 1994:171–6. [32] Suresh S, Omkar SN, Mani V, Guruprakash TN. Lift coefficient prediction at high angle of attack using recurrent neural networks. Aerospace Sci Technol 2003;7(8):595–602.