Inverse analysis of cyclic constitutive models for unsaturated soil under consideration of oscillating functions

In order to assess the probability of foundation failure resulting from cyclic action on structures and to minimize the prediction error, various existing constitutive models considering cyclic loaded dry soils were extended to unsaturated soil conditions by the authors, thus requiring further calibration during application on existing slightly variable soil condition as well as the soil heterogeneities. The efficiency and effectiveness of these models is majorly influenced by the cyclic constitutive parameters and the soil suction. Little or no details exist in literature about the model based identification and the calibration of the constitutive parameters under cyclic loaded soils. This could be attributed to the difficulties and complexities of the inverse modeling of such complex phenomena. A wide variety of optimization strategies for the solution of the sum of least-squares problems as usually done in the field of model calibration exists, however the inverse analysis of the unsaturated soil response under oscillatory load functions has not been solved up to now. This paper gives insight into the model calibration challenges and also puts forward advanced optimization methods for the inverse modeling of cyclic loaded foundation response on unsaturated soils.


Introduction
The increased awareness in the impact of cyclic actions from natural and man-made sources on a structure during its lifespan necessitates understanding the complex soil structure interaction during these cyclic loads.The stresses induced by these cyclic loads lead to failure patterns such as flow-type failure (liquefaction), cyclic mobility or plastic strain accumulation.The failure type experienced depends on the initial stress state and degree of stress reversal as observed by [1].In order to design structures such that the effects of these stresses are minimized, various constitutive models based on different formulations such as -plasticity, elasto-plasticity, hypoplasticity and so on, to describe the soil structure interaction during complex cyclic loading on dry, saturated and unsaturated soils have been developed.The quality of these model responses is a function of the model-based parameter identification from experimental data.Various strategies based on different concpets exist in finding out the values of these parameters, however the complexities involved in modeling the soil response to cyclic loading limits the application of some of these strategies specifically with regards to the quality and computation effort required for the inverse analysis.
The last three decades heralded the development and extension of many constitutive models some of which includes the elasto-plastic strain hardening law model by [2] for monotonic loads which served as a foundation for many other models such as the model by [3] for cyclic loads and model by [4,5] for monotonic and cyclic loaded footings on unsaturated soil which takes into account the effect of matric suction on the bearing capacity of the soil as observed by [6].The models in [3] and [5] for cyclic loading is capable of simulating the foundation response under different loading conditions and soil states.However, the quality of the model's response is strongly influenced by certain constitutive parameters, thus in order to apply the models for the actualization of economical and efficient structural design, we have to bridge the gap between model simulation and experimental response.This is achieved by the calibration and optimization of the models' sensitive 'cyclic constitutive parameters' with experimental results with the aim of identifying parameters which are unique and robust with respect to change in experimental data.Thus, the implementation of an efficient, stable and robust algorithm is required for the solution of the complex ill-posed inverse problem for unsaturated soils under oscillatory load functions with minimum error and computational cost.
The importance of these constitutive models cannot be over-emphasized, especially when considering the ease with which it can be incorporated into finite element software for structural analysis, structural design and struc-tural health monitoring as well, therefore in order to yield reliable model responses an in-depth knowledge on the cyclic parameter behavior is required in addition to selecting an efficient optimization strategy.

Mathematical Model
The elasto-plastic strain hardening model proposed by [3] for response to cyclic loading on dry soil is an extension of the model proposed by [2] which was developed for monotonic loading.The plastic strain as a result of monotonic load is expressed in terms of generalized force and displacement as shown in Equation 1, where Λ according to [? ] is the plastic multiplier and it depends on the positive change in loading, the loading history, ρ c and the macroelement through the hardening modulus, H.
The generalized force Q consists of the vertical force V, horizontal force H and moment M, normalized with the maximum centric force the soil can resist V m .The generalized displacement q consists of the vertical displacement v, the horizontal displacement u, and the rotations θ normalized by V m .
In order to simulate cyclic loading, the model in [3] permits the formation of plastic strains inside the bounding surface.The magnitude of plastic strains induced is a function of the distance between the stress point and the image point I p on the bounding surface, the memory parameters (ρ c and ρ k ) and the cyclic constitutive parameters expressed in Φ in Equation 4. The parameters Λ(Q I ) and δg δQ (Q I ) are calculated on the image point and not on the current stress point as would be the case for monotonic loading The term Φ which is a function of δ and ρ k is a diagonal matrix (Equation 5) which acts as a weighting function during cyclic loading.Larger values of δ leads to small diagonal terms and also smaller plastic displacements [3] Plastic strain caused by cyclic loading is sensitive to the constitutive parameters ς i and κ i (where i = η, , ζ) in the Φ matrix.These parameters regulate the rate at which shakedown occurs in the model and also the the coupling between the horizontal loading and the corresponding vertical displacement at constant vertical force.The extension of these model to incorporate the multiphase soil behaviour was proposed by [5].It was observed that not only the intrinsic soil properties like maximum centric vertical force V m and initial slope of load deformation curve (R 0 ) are sensitive to matric suction but also the cyclic constitutive parameters.The proposed model expressed in Equation 6 is similar to Equation 4, except that the parameters are a function of 'ψ'

Calibration and optimization
The quality of any prognosis depends on the quality of the identified cyclic constitutive parameters.Besides accuracy, efficiency is also an important issue as dynamic non-linear models tend to be complex and time consuming, particularly during calibration.Therefore it was necessary to study the cyclic constitutive parameter dependency using contour lines before selecting an optimization strategy.A strong dependency between some parameters was observed and is expressed in Figure 1.The presence of many valleys (indicating minima) can be traced to the structure of the Φ matrix which has addition operations in each of the diagonal terms as seen in Equation 5. Thus different pairs of cyclic constitutive parameters in each diagonal term in the Φ matrix can yield similar results with regards to model response.Solving these obvious non-convex inverse problem by minimizing the conventional sum of squared error between the calculated vertical cyclic displacement, η [3] from the forward model and the experimental response in Equation 7 yields non-unique parameters which are a function of the initial value used to initiate the optimization process.The Nelder-Mead method ([10]), which is heuristic based, deterministic and known to be robust with respect to noisy cost function is applied in identifying these cyclic constitutive parameters.
where: -C f is the cost function to be minimized η mod is the vertical component of the generalized cyclic plastic displacement vector calculated by the model η exp is the vertical component of the generalized cyclic plastic displacement vector obtained from the experiment p is a vector of Λ I (ψ), Φ(δ, ρ k , ψ) and δg(ψ) δQ(ψ) I t i is the time step at each point n is the total number of data points In order to increase the probability of reaching the global minima and getting best quality results with the Nelder-Mead simplex algorithm, it was necessary to vary the starting point of the optimization until the best fit or lowest value of C f is attained.However, an attempt to obtain a more reliable lowest function value and unique parameters from the optimization, convexification of the problem (as seen from the multiple minimas in Figure 1) can be achieved by introducing a penalty term to the cost function (Equation 7)in line with Tikhonov regularization.This regularization/penalty process is influenced also by a priori knowledge and can be observed in Equation 8. where: -C f is the cost function to be minimized η mod is the vertical component of the generalized cyclic plastic displacement vector calculated by the model η exp is the vertical component of the generalized cyclic plastic displacement vector obtained from the experiment p is a vector of Λ I (ψ), Φ(δ, ρ k , ψ) and δg(ψ) δQ(ψ) I t i is the time step at each point  n is the total number of data points -Γ is a scaling factor for the penalty term -Ψ 0 is a vector of a priori cyclic constitutive parameters -Ψ mod is a vector of calculated cyclic constitutive parameters k is the total number of cyclic constitutive parameters The effect of this regularization can be observed in Figure 2. The value of Γ should be carefully chosen because it greatly influences the properties of the regularized solution.A large value of Γ results in over-smoothing and also large residuals, whereas very small Γ values (undersmooothing) may give a good fit but the solutions would be dominated by data errors.An optimal value for Γ is obtained from a plot of the residual norm and the solution norm (L-curve) according to [11].

Discussion of results
Model parameter identification was carried out by calibrating the model to experimental data obtained from [5].Before the model calibration, it was necessary to filter out measurement errors from the experimental data.The model parameter identification was carried out using the Nelder-Mead algorithm on saturated and unsaturated soils, with outcomes shown in Table 1 and Figure 3 to Figure 4, where a comparison between the experimental result is made with the optimized model response.
It is worth noting that as a result of the strong dependency between the cyclic constitutive parameters and the application of a deterministic local optimization strategy, the values used to initiate the optimization plays a major role in the computational time required and the quality of the output result obtained.Although values from literature (values published in [12]) were first used to initiate the optimization process, the quality of fit was not satisfactory, thus other randomly generated values were used as initial values until a best fit was attained.Furthermore, the outcome of the inverse modeling revealed a lower sensitivity to the parameters ς and κ .This is as a results of       the experimental data being dominated by vertical oscillatory loads as obtained in Table 2 where although different values were used to initiate the optimization algorithm, the cyclic constitutive parameters ς η and κ η related to η in the Φ matrix (Equation 5) remained constant while parameters related to and the computational cost varied.
Moreover, in some particular cases the solution to the inverse problem yields a good quality response, however the identified parameters are unrealistic as observed for the soil sample with ψ = 0 kPa and Γ = 0 in Table 1.An attempt to yield realistic parameters was made by implementing a modified Nelder-Mead algorithm to constrain the resulting optimized parameters.Although realistic values of parameters were obtained, a decrease in the quality of the model response was observed.An efficient and ro-

Conclusion
As observed in Figure 3 and Figure 4 and also corroborated in literature, the unsaturated soil yields lower plastic displacements when compared to other soil states, in addition to this and with respect identifying the model parameters it can be concluded that • Measurement error (noise) affects quality of fit during calibration and thus have to be filtered out before parameter identification • The initial guess is a key factor in determining the outcome (quality of parameters obtained) and efficiency of the the inverse modeling computation.
• The model exhibited a higher sensitivity towards parameters are related to the vertical component of the model.This is as a result of the experimental data being dominated by vertical oscillatory loads.
• Though the Nelder-Mead algorithm may converge to a local minima, the process of using random initial values and comparing the quality of fit shows the application of this algorithm to find the global minima with a high probability.
• Better quality fit was observed for unsaturated soil condition with minimal computational effort.This can be attributed to the fact that shakedown occurs earlier and as such the effect of the slight reduction in the accumulated plastic displacement during unloading by soil rebound which the model does not account for is minimized.
• With regularization, calibration yields more reliable results with minimal computational effort for these complex dynamic soil models as seen in Table 1.

Figure 3 .
Figure 3.Comparison between the optimized model response and experimental result for saturated soil (0[kPa]) Cyclic force-displacement plot (b) Accumulated plastic displacement plot

Figure 4 .
Figure 4. Comparison between the optimized model response and experimental result for unsaturated soil (1.7[kPa])

Table 1 .
Parameters obtained after model optimization with and without regularization

Table 2 .
Initial value and optimized constitutive parameters UNSATbust strategy to tackle this issue was to incorporate the regularization term Γ into the objective function as obtained in Equation 8 and explained in subsection 2.2.The process of convexification not only tackles the issue of unrealistic parameters but it also reduces the computational costs and yields unique and stable parameters irrespective of the values used to initiate the optimization.