Reducing numerical work in non-linear parameter identification
Emoke Imre, Peter Berzi, Csaba Hegedus, Sandor Kovacs, Levente Kovacs
RReducing numerical work in non-linear parameter identification
Emoke Imre , Peter Berzi , Csaba Hegedus , Sandor Kovacs , Levente Kovacs HBM Systems research center, Obuda University, Hungary AIAM Doctoral School, Obuda University, Hungary Department of Numerical Analysis, Eotvos Lorand University, Budapest, Hungary * Corresponding author, e-mail: [email protected]
Abstract
The real-life merit functions have an unimaginable complexity of an M -dimensional topography, where M is the number of the parameters. It is shown that there is an underlying noise-free merit function, called follower merit function which can be constructed from simulated, noise-free data using the solution of a Least Squares minimization. The difference of the real-life merit and follower merit function is controlled by the norm of the noise-free error vector. Therefore, the non-important local minima arisen from the noise can be skipped during such minimization that apply adjustable, large step sizes, depending on the norm of the noise-free error vector. The non-linear minimization can be difficult even in the case of some strictly convex, non-negative merit functions if the geometry is highly asymmetric or if the global minimum is quasi-degenerated. The suggested hierarchical minimisation - based on the implicit function of a projection map - makes the minimisation in two steps. In the first step some parameters are eliminated with conditional minimisation, resulting in a kind of deepest section of the merit function, called clever section with respect to the remainder parameters. Then the clever section is minimized. The method can be used to handle a quasi-degenerate global minimum, in sensitivity analyses, reliability testing. Keywords non-linear minimization, implicit function theorem, sub-minimization, inverse problem, parameter identification, least-squares merit function
Let us assume that the model is given by a function u :R + x R M → R; ( t, p ) → u ( t, p ); where R + t is the time variable, D p is the parameter vector with M elements ; the D is the simply connected and compact parameter domain contained by R M . We can introduce a partial function v p ( t ) = u ( t, p ) and define its graph G :R M → R as F = G ( p ) = { R ( t,v p ( t ))}. The function v p ( t ) = u ( t, p ) is the so-called model solution function. In the direct problem - for a fixed value of p - the graph of the partial function v p ( t ) = u ( t, p ), expressed as F = G ( p ), is determined. In the inverse problem the goal is the determination of the parameter vector p for a given graph F, expressed as p = G -1 ( F ). In the discrete formulation, only some points of the graph of the partial function v p ( t ) are considered, denoted by v p ( t i ), i =1.. N . These values are compiled into a so-called data vector f . The measured data vector: == )( )( )( ....21 N tf tf tf m uf (1) where t i are the sampling times ( i = 1,2.. N ). Simulated data vector using a fixed value of the parameter vector: = ),( ),( ),( ....21 ***, pppf N tu tu tu (2) The function u ( p ): R M → R N , ( p ) → u ( p ) , called model response function, containing N elements, can be constructed from the copies of u ( t, p ) in such a way that times t ...t N - the elements of the so-called sampling time vector t - are substituted for the time variable t in u ( t , p ). = ),( ),( ),() ....21 pppu(p N tu tu tu (3) In the inverse problem, for a fixed value of f and t , the goal is the determination of the parameter vector p , from the equation h ( p ) = f - u ( p ) = 0, expressed symbolically as p = u -1 ( f ). === ),( ),( ),()( )( )()()( ....21....21 NN tu tu tutf tf tf (4) The solution of a system of non-linear of equations f - u ( p ) = 0 is needed, which generally has weak solution only, since generally N > M. The weak solution is the global minimizing vector p min of the least squares objective function called real-life merit function which is the norm squares of the error vector h ( t, p ) = u m ( t i ) - u ( t , p ): min!phphp == )()()( T F (5) The general linear inverse problem is to fit a set of data points to a model that is not just a linear combination of 1 and t (namely a + b t ), but rather a linear combination of any M specified functions of t . The general form of this kind of model is as follows: == i iiMk tfpppptftftfu )()()()()( p (6) In the general linear inverse problem, the model response function: == ),( ),( ),()()()( )()()( )()()()( ....21212....1 22221 11211 ppppu NMNkNN kk tu tu tuppptftftf tftftf tftftf (7)
The error vector and the direct formulation of the inverse problem: === ),( ),( ),()( )( )()()( ....21....21 NN tu tu tutf tf tf (8) === )( )( )()()()( )()()( )()()( ....21212....1 22221 11211 NMNkNN kk tf tf tfppptftftf tftftf tftftf fpA (9) The Gauss Normal Equations: fApAA TT = (10) The merit function: ( ) ( ) ( ) ( ) ( ) fpAfpA p,p u-fp u-fp −−= −== = T 21T )()()()()(
Ni ii tutfF (11)
Statement
In the general linear inverse problem, the second derivative of the merit function with respect to p is globally positive definite, being a constant Gram matrix. Proof
The derivative is a constant Gram matrix (see Eq 9):
AAp T = F (12) The parameter space is split into the direct sum of two subspaces:
RD RD, kMk = − ,][ ppp,pp (13) The partly linear inverse problem is to fit a set of data points to a model that is a linear combination of any M specified functions of t and p . The general form of this kind of model is == i iikk tfpppptftftfu ),()(),(),()( ,1,1 2,1 1,121 ppppp (14) In the partly linear inverse problem, the model response function: ppA ppp ppp pppp u )( ),(),(),( ),(),(),( ),(),(),()( ,1 2,1 1,12....1 22221 11211 = = kNkNN kk ppptftftf tftftf tftftf (15)
In the partly linear inverse problem, the error vector and the direct formulation of the conditional inverse problem: == ),( ),( ),()( )( )()( ....21....21 NN tu tu tutf tf tf (16) In the partly linear inverse problem, the Gauss Normal Equations of the conditional inverse problem for p , at fixed p : fpAppApA T212T2 )()()( = (17) In the partly linear inverse problem, the merit function can be written as follows: ( ) ( ) ( ) ( ) ( ) fppAfppA p,p u-fp u-fp −−= −== = )()( )()()()()( T 21T
Ni ii tutfF (18)
Statement
In the partly linear inverse problem, the second derivative of the merit function with respect to p is globally positive definit, being a constant Gram matrix at the fixed values of p . Proof
The derivative with respect to p is a constant Gram matrix in case of the fixed value of p (see Eq 15): )()(2'')( pApAp Tp = F (19) The multidimensional, non-linear minimization generally needs iterations at. The difficulty is generally more than this since (i) too many local minima of the real-life merit function may occur due to the noise, (ii) the geometry of the valley of the global minimum can be asymmetric and (iii) the minimum can be quasi-degenerated. The first problem - the unimaginable complexity of an N -dimensional topography [1] - is treated in the literature as follows. A trial and error process is used, when not significant decrease in the error is attained, the process may restart from a new starting vector. The non-linear minimisation methods, which include one-dimensional minimisation, are generally used together with some new starting vector determination methods (eg,. the simulated annealing method, which steps away from a local minimum using probability theory methods, see eg. [1] or by fitting a hypersurface locally [2]). The non-linear minimization methods without one-dimensional minimisation (eg., simplex method [1], or secant method see Section 6) are considered as effective tool in case of a few parameters [1], the reason of the better performance is related to the fact that these may skip some critical points, as shown here. The quasi-degenerated minimum problem is treated generally by regularization ([3]), the asymmetric shape of the vally is an unsolved problem. The singulariry of the Hesse matrix at a minimum p min gives information on local uniqueness. If the global minimum is degenerated, infinite many solutions may exists, some parameters are dependent. In case of linear inverse problems, the global uniqueness is met in case of local non-singularity. The solution of a real-life inverse problem p min is approximate with some error (see App. A). Using the linear model, the computation of the standard deviation of the parameters is possible, assuming normally distributed noise (residuals) [1]. The geometrical interpretation is related to the orthogonal projection of a confidence domain, which is defined by a level line of a kind of noise-free merit function (Fig. 1), with symmetry related to the individual parameters. The unsolved problem is the parameter error in case of non-normal noise. In case of non-linear inverse problems, the global uniqueness is not met in case of local non-singularity. The computation of the standard deviation is possible, similarly to the linear case, assuming normally distributed noise [1]. It has similar geometrical meaning: it is the orthogonal projection of the confidence domain, related to a noise-free, linearized merit function. Some unsolved, additional questions are the global uniqueness, and the parameter error in case of non-normally distributed noise and the non-symmetry with respect to the individual parameters.
Fig. 1
Confidence domain of 68.3%, the standard deviation in case of normally distributed noise (after [1]).
To overcome the foregoing difficulties, three new concepts are presented and discussed. The first concept is the Hierarchical Inverse Problem, which means split minimisation, such that different minimisation methods may be used in each part. The hierarchical inverse problem concept is mathematically explained by the implicit function theorem. The related map is the orthogonal projection of the graph of the merit function onto a plane < F , x >, and the implicit function y = g ( x ) is related to the fold of the image of the map. The second concept is related to the use of the section of the merit function along the fold, defined in the HIP concept. The so-called clever section can be used to estimate the global minimum without minimization if some independent information for the x min part of the global minimum is available. In this case y min = g ( x min ). The one-dimensional clever section - depending on a single parameter - is a deepest sensitivity section which can be used for sensitivity analyses and reliability testing. The third concept is the follower, noise-free merit function - the best possible noiseless approximation of the real-life merit function - which can equally be used for minimisation and reliability testing. The real-life merit function has too many local minima with nearly the same merit function value around the global minimiser p min . The global uniqueness of the solution can be interpreted by the follower here. The obstruction domain for minimisation is a parameter error domain, defined by a level line of the follower merit function here. The follower is not known a priori but can be determined afterwards using the global minimiser p min . It is shown that the difference of the two functions is dependent on the actual size of the noise-free residuals. Due to the similarity, the decreasing direction of the follower can be estimated during minimisation by the noisy (real-life) merit function on a part of the parameter domain - aiding minimisation - if the step size is large enough, eg. dependent on the actual size of the noise-free residuals. In the paper two examples are presented. The step size of the secant method is qualitatively analysed. The secant method is a non-linear minimization method, without one-dimensional minimisation ([6 to 8]), where the step size is determined by a linear equation in its original form. It is shown that this method is theoretically able to skip some local minima since its step size is dependent on the actual size of the noise-polluted residuals. The example also illustrates how theoretically the linearly dependent parameters can be eliminated in hierarchical minimisation in this special method. The second example shows an automatic global minimisation related to the injective solutions of linear PDE-s, giving reliability information as well. The hierarchical inverse problem concept was applied for a topography mapping with a coordinate map. The fact is used that the strict convexity of the follower merit functions related to injective solutions of system of coupled, linear partial differential equations may occur. Hyperplanes generated by coordinate a mesh points can bracket decreasing convex level surfaces of these. The injective solution was prepared for model fitting including a parametric function for the initial condition and by optionally adding a non-linearity in the form of an analytic function. Some model-variants were elaborated on the basis of a “maximal” model expressed with 8 parameters, 4 non-linearly dependent parameters („model solution H”). In its variants the number of the non-linearly dependent parameters was decreased up to 1 by using information concerning the analytical properties of the model (see eg., „model solution HC”). The number of the function value evaluation depended on the number of non-linearly dependent parameters exponentially, in case of a single non-linearly dependent parameters, it was less than 60, increasing by 1 order of magnitude with an additional non-linearly dependent parameter. The parameter space can be split into the direct sum of two subspaces:
RD RD,
JMJ = − ,][ ppp,pp (20) The first inverse problem is a conditional minimization: M RDF − == const, ,min!),( pppp (21) the solution of which is the relation p = a ( p ). If this inserted then an M - J dimensional minimal section of p the section of the original Least Squares merit function, denoted by F [ a ( p ),( p )] is resulted which depends only on p . The second inverse problem: RD F JM = − min!]),([ pppa (22) is the minimization of F ( a ( p ),( p )). The section F [ a ( p ),( p )] is called minimal or clever section of p The solution of the hierarchical inverse problem is trivially equivalent to the original one if the Hessian of the Least Squares merit function F with respect to p is globally positive definite, which is met if the function F is strictly convex, or in case of linear dependence in some parameters. The implicit function theorem is used for the interpretation, which is described in Appendix B. We change the notation so that subscripts can be omitted in the index. We use y and x instead of p and p in the direct sum decomposition p = ( x , y )= ( p , p ). Let us consider a strictly convex, non-negative, analytic Least Squares merit function F ( p ): D < R n + m → R with simply connected, compact domain D and with the arbitrary p = ( x , y ) arbitrary direct sum decomposition of D . The global minimum is ( a , b ), where F has the only critical point. Let us consider the orthogonal projection of the graph of F {( x, y ), F ( x , y )} onto the F - x coordinate plan in the total Euclidean space generated by the graph. The projection planes are x = constant coordinate planes. The critical values of the projection map - where the fold of the projection is found - are related to the F y ’ ( x, y ) = f ( x, y ) = . The implicit function theorem [4, 5] is applied in the case of the projection of the graph of a strictly convex, non-negative, analytic (Least Squares merit) function, to describe the implicit function related to the zero condition of the y -derivative F y = f ( x, y ) = as follows. Starting from a direct product decomposition p = ( x , y ), the goal is to construct the function g : R n → R m whose graph ( x , g ( x )) is precisely the set of all ( x , y ) such that F y ’ = f ( x, y ) = . Let us consider the critical point at the global minimum of the merit function ( a , b ) = ( a , ..., a n , b , ..., b m ) where F y ’ ( x, y ) = f ( x, y ) = , the zero is element of R m . If the matrix F y y ’’ ( x, y ) = [(∂ f i /∂ y j ) ( a , b )] has full rank, then there exists an open set U containing a , an open set V containing b , and a unique, analytic, injective function g : U → V such that == )(,),()(, UU ( ) Whenever we have the hypothesis that f is analytic inside U × V , then the same holds true for the implicit function g inside U. The implicit function is globally unique if F y y ’’ ( x, y ) = [(∂ f i /∂ y j ) ( a , b )] is globally positive definite. This condition is met - the Hessian is globally positive definite - in case of strict convexity for any direct product decomposition p = ( x , y ). Two consequences are as follows. The implicit functions related to various a direct product decomposition have the same lattice structure as the subspaces of the R M in the direct product decompositions. If the split is such that x’< x then the clever section F [ g ( x ), x] contains F [ g ( x’ ), x’] . It can be hypothetised that the real, analytic implicit function is injective at the single critical point (i.e. the global minimum of F ). Therefore, if the parameter vector x min is the part global minimiser vector, then the remainder part of the global minimiser vector can be computed as follows: y min = g ( x min ). The follower merit function F’ ( p ) is defined as the norm squares of the simulated error vector, computed by using the simulated, noise-free data, generated with p min : where u ( t i ,, p min ), the simulated, noise-free data and u ( t, p ) is the model solution, t i (i = 1.. N ) are sampling times, the notation is not related to derivation, it means noise-free. '( ) ( ) ( , ) ( , ) N Ni i ii i
F h u t u t = = = = − min p p p p (24) .1.2 The definition of the noise
The noise z i can be expressed from the following decomposition, after the determination of the global minimum of the merit function and the simulated data: ( ) ( , ) m i i i u t u t z = + min p (25) where p min minimiser, u ( t, p ) is the model solution. The level surface of the noise-free merit function related to the global minimum of the real-life merit function is used to define the parameter error domain, similarly to the statistical error domain (Fig. 1). The parameter error domain is defined as the follows, by a sublevel set of the follower merit function: )()(' min ppp
FcFDE == (26) In the parameter error domain the real-life merit function contains several local minima with about the same value. In the complement domain - far from the global minimum - is the following : )(' zh p (27) is met. This domain is defined as the so-called similarity domain, where the follower merit function can be approximated by the real-life merit function. The relation of the two merit functions
Statement
Using the foregoing, the two merit functions: the follower merit function F and the real-life merit function F has the following relationship:
21 1 ( ) ( ) '( ) 2 ( ) ( )
N Ni i i ii i
F h z F h z F = = = − = − + min p p p p p ( ) By using the Cauchy-Schwartz-Bunyakovsky inequality and the difference of the two functions is as follows: )()('2)(' )('2)(')( )('2)(')()('2)(' min 22 22min pphp phphp phphpphp FzF zzF zzFzF ++ =++ +−=+− ( ) )('2)(')()('2 zzFFzz +−+− phppph (30) According to Eq ( 30), the difference of the two function seems to be small for fixed, small noise, being controlled however by the norm of the noise-free error vector (denoted by a dash in this section).
5. The practical use of concepts and methods 5.1
Hierarchical inverse problem
The most important practical application of the split minimisation is the elimination of the linear part of parameter vector through sub-minimisation, being the related Hesse matrix globally positive definite. The linear sub-minimisation - a generalized linear inverse problem [1] - can be solved without iteration, eg., by the SVD or singular value decomposition (SVD) algorithm. The linearly dependent parameters can be eliminated by sub-minimisation before each function value evaluation in case of any non-linear minimization method.
The clever sections of convex merit functions
The implicit function g ( x ) can be used to estimate the global minimiser of F [ x , g ( x )] without minimization in case the value of x min is available from independent information. In this case y min = g ( x min ). This follows from the fact that the implicit function is injective at the single critical poinst. The one-dimensional clever sections of the follower merit function and the real-life merit functions can simultaneously be represented in the 2-dimensional plane, reflecting uniqueness and parameter error (Fig.2). Being the error domain non-symmetric, the orthogonal projection of its diameter onto parameter axis p i is reasonable to be used to characterize the parameter error. Minimisation with the follower
Statement
The follower and the real-life merit functions are pseudo-metrics. The metrics related to the real-life merit function F is introduced as follows: )()(),( baba FFd F −= (31) where D a,b are parameter vectors. In case b is the global minimum, the merit function value indicates the distance from it, being a pseudo-metrics. The definition of a metric ( ) ( )( ) ( ) ( ) ( ) ( ) ,,,)4,,,)3 ,0,)2,0,)1 ddddd dd += == (32) If the second condition is not met – which is the case of the norm of the error vector –, a pseudo- metric is defined. Consequences
According to Eq ( 30), the difference between the follower and the real-life merit functions is bounded by a term which is linearly related to the norm of the noise-free, ollower at a fixed noise vector. If the norm of the noise is relatively small, then theoretically (i) the follower merit function values can be approximated by real-life merit function values, (ii) decreasing directions can be assessed, (iii) the step size can be linked to the difference between the follower and the real-life merit functions. In the non-linear parameter identification methods not using line minimization, the step size can theoretically be selected to be large enough to skip the critical points due to the noise. In case of strictly convex follower merit function, some decreasing convex level surfaces (and a minimum) may approximately be bracketed on the similarity domain, by global coordinate grid, using real-life merit function values. The same global coordinate grid can be used - with a projection algorithm - to determine the one-dimensional clever sections. The grid dimension is determined by the number of the non-linearly dependent parameters, since the linearly dependent parameters can be eliminated. The real-life merit functions values are computed in the points of a coordinate grid in a nested cycle. The one-dimensional clever section point of each parameter x i can be determined by projections along each respective coordinate planes x i , j =constant (Fig. 3). Once the minimum of the real-life merit function is determined, the same algorithm can be used for the follower merit functions, too. In a local version, in the case of a general follower merit function, a local coordinate grid can be generated in to find a decreasing direction and a new starting vector. (a) (b) Fig. 2.
The one dimensional clever section of merit functions and the parameter error domain (a) in the total space (b) in the parameter space. See a parameter error value in App. A.
Fig. 3.
The determination of the approxinmate one-dimensional clever sections by elimination, projection and minimum search in a nested cycle system. A coordinate mesh is generated in the M-J dimensional space of the nonlinearly dependent parameters, J linearly dependent parameters are eliminated. The condition of the termination of the nested cycle is the visiting if all grid point.
The secant-type methods – where the variable step size is controlled by the real-life merit function F value (i.e., norm square of the noise-polluted error vector) – are shortly presented as such methods that are hypothetically able to tackle the effect of the noise and proceed quick into the vicinity of the global minimum In the secant method [6 to 8] n+ trial solutions for the parameter vector and the related residual or error vector ( p i , h ( p i )) ( i = 0,...,n ) are used to compute a new parameter vector: ( ) ( ) qAAphpphp + rp (33) where m is the number of sampling times, n is parameter number. A p = [ p - p , ..., p - p n ], A r = [ h ( p )- h ( p )-, ..., h ( p )- h ( p n )] and q T = [ q , ..., q n ] is a vector with n coefficients of secant method, a so-called an auxiliary variable . The q in the one-dimensional case is equal to h ( p )/ [ h ( p ) - h ( p )], similar is valid in multidimensional case for the elements of q . The absolute value of the multiplier h ( p ) is a pseudometric, expressing a distance from the global minimum. The step size is comparable with the controlling term of difference between the follower and the real-life merit function (see Eq (30)). Far from the global minimum, a large step from p is resulted and vice versa. a) (b) Fig. 4 . (a) Secant method in 1D. Note that if the function value at x (a) (b) Fig. 5 . Flow diagrams (a) modified secant method, the termination criterion is the sufficiently small change in the parameters; (b) the hypothetical completion with linear sub-minimisation (each function value evaluation is preceeded by a linear sub-minimisation). The question is that the linear sub-minimisation may amend the method being the secant method linear in each iteration cycle.
Wolfe ([5]) suggested to select n+ new trial solutions for the next iteration by dropping the trial solution for which the norm of the error vector is the largest. This procedure fails when the norm of the new residual vector is the largest one. In the amended method, after determining p n+1 , a second trial solution p n+2 is determined. With matrix equation: ( ) −= + ph ppqAAT0 0T rprP (34) where T p and T r are transformation diagonal matrices with elements: )(/)(, pp iniirii iiip hhtpp ppt ++ ++ =−−= (35) The distance of p n+1 and p n+2 is related to the decrease of the error and the merit function between p n+1 and p . The n +1 new trial solutions are selected from the new vectors as: ( ) inini +++ += p-pepp i1n (36) where e i is the i-th unit vector. The system of linear equations in the second rows of the matrix equations are solved by singular value decomposition (SVD). The new auxiliary variable q , is determined from the first line of the matrix equation (17). In the one-dimensional case, it is a multiple of q , being equal to h ( p )/ h ( p )* q [8]. A similar generalization is valid in the multi-dimensional case, described in [6 to 8] (Fig. 5). In the secant method, the step size (related to the q coefficient) depends on the norm of the residual vector. It can be noted that the norm square of the residual vector is a pseudo-metric, expressing a distance from the global minimum. The first consequence is that the step size depends on the distance from the global minimum, the step sizes will be larger farer from the global minimum and vice versa. The second consequence is that this term is found in Eq ( 30), too, which means that the effect of the noise could theoretically be tackled . The elimination of the linear part of parameter vector through sub-minimisation is indicated in Figure 5(b) technically, which can be added to any non-linear minimisation algorithm. However, being the secant method linear, this completion may not be necessary, the linear part of the parameter vector is probably automatically be solved in the first steps. able 1 Types and notation of one dimensional oedometric dissipation tests with constant boundary condition (Multistage) oedometric relaxation test MRT or ORT (Multistage) oedometric compression test MCT or OCT
Table 2
Types dissipation tests made with static penetrometers, modelled with cylindrical and spherical (ellipsoid) shaped domain. Measured variable dissipation test made by Pore water pressure dissipation test, sensor on the shaft and/or on the tip CPTu (static cone penetrometer) Total stress dissipation test, sensor on the shaft and/or on the tip CPT- piezo-lateral stress cell, CPT- DMT tip, CPT - CPT- qc Effective stress dissipation test, sensor on the shaft CPT- fs Table 3 Point-symmetric consolidation models [12] V o r boundary condition
1D (Oedometric models) 2D (Cylindrical pile models) 3D (Spherical pile models) no (uncoupled) Terzaghi (1923 ) Soderberg (1962) Torstensson (1975) v-v (coupled 1) Imre ( 1997-1999) Imre & Rózsa (1998) Imre & Rózsa (2002) v- (coupled 2) Biot (1941) Randolph at al (1979) [ Imre & Rózsa (2005) rr rr Fig. 6.
The displacement domain bounded by a (a) 0 dimensional sphere (oedometer model), (b) 1 dimensional sphere (cylindrical model), (c) 2 dimensional sphere (spherical model). Example 2 - dissipation test evaluation
The objective of the study was to elaborate automatic evaluaton methods for such dissipation tests of Geotechnics that can be evaluated with linear, point-symmetric, coupled consolidation models [9 to 16] (Tables 1 to 3). The follower merit function related to these injective, analytical consolidation model solutions, with fixed, convex initial condition and sufficiently long data series is strictly convex, and the related implicit function is injective, according to the experiences [11]. However, by the identification of the initial condition on a function class allowing non-convex analytic functions and by modelling of relaxation (time dependent constitutive law), this geometry may slightly change, but remains near-convex. The same method - applying the global coordinate grid projection algorithm - was applicable in every case (software is available for each case [17 to 22]). In case of the standard multistage oedometric compression tests (MCT) and its dual counterpart, the multistage oedometric relaxation tests (MRT), a short testing procedure (with shorter stages than 99% of the dissipation time except for the last stage) was suggested. Laboratory tests were made on identical soil samples with different plasticity [9, 10]. The model validation showed that (i) the fitting error decreased to half for those model-variants where the simultaneous creep/relaxation was considered, (ii) the relaxation qualitatively influenced the model response in a parameter dependent way, and (iii) the short stage data were successfully evaluated using the clever section of the single parameter c . In this paper a summary is given on the evaluation of the multistage oedometric relaxation tests MRT. It is shown, that by using the suggested technique, a class of inverse problems - being related to (near-)convex follower merit functions - becomes easily soluble without iteration or regularization, the reliability criteria of the solution can be be tested. The analytical solution of the linear, point-symmetric coupled consolidation model of the oedometric relaxation test for the so-called total stress and the so-called pore water pressure u (Appendix C, [8 to 11]): + e - = t T i-i1=i )( (37) e H iy - = yt,u T 2 i2-i1=i ]1)/([cos()( − (38) where T = ct / H is the time factor, c is the coefficient of consolidation, i are Fourier coefficients related to the initial condition, H is the half-width of a double-drained sample (length between the boundaries of the model). The linear, coupled consolidation model has some special features, which can be summarized as follows ([8 to 11]). eing the volume constant during a stage, the mean effective stress ( mean ) is equal to the final total stress ( ) in case of linear constitutive law. It is a model parameter, depending on the displacement load v and the H the sample height: )( '0 tHvE meanoed = (39) The time dependent constitutive law was included in an approximate form into the model in its most complex form, by adding a relaxation term to the total stress. The relaxation term includes a partial unloading effect. The equations used for model fitting were as follows. )()( t - + e - = t r T 2 i2-i1=i (40) The empirical relaxation term: t>t;t+t t+tsb-1 1(0) s= (t)
331 1r log (41) where s is the coefficient of relaxation, t is the pause of relaxation in case of partial unloading, b = log(( t + t )/ t ), t is time parameter. The initial condition was identified using the following parametric shape function in model H: y C +y B + y A = y)(0,u (42) where A , B and C are real valued initial condition parameters, with linear dependence. The related Fourier coefficients were included in the model. The ‘model H’ contained 8 parameters, 4 out them were with non-linear dependence (Table 4). The linearly dependent A , B and C parameters equally allowed convex, concave or partly convex and concave initial condition functions. Several model-versions were elaborated, differing in the number of the non-linearly dependent parameters, using additional information on the consolidation part-model and the relaxation part-model (described in Appendix C), in here only ‘model HC’ is mentioned where no relaxation was modelled. It can be noted that parameters σ and c have different meaning if relaxation is modelled or if not. If the relaxation is modelled, the compression curve constructed from the identified σ is time independent, this the so-called instant compression curve. However, if the relaxation is not modelled (ie., for model-version HC), the identified σ is depending on the stage duration, as it is dependent in the measurements (as it is well-known from the oedometer testing). The model-version H was fitted on measured data with the conjugate gradient method. However, the solution was not possible to be determined due to the large number of critical points. Therefore, a local grid system was used to determine some new starting vectors. The grids were 2-dimensional, considering the most important parameters ( c , initial condition). The conjugate gradient method was used with 5 to 8 parameters. The model-version H and HC were fitted on measured data with the projection algorithm including linear sub-minimisation [8 to 11], in the frame of an automatic global minimization (Fig. 3). The global coordinate grid was generated in the subspace of the non-linearly dependent parameters, the linearly dependent parameters were eliminated by sub-minimisation. The grid size was defined by using the half error domain diameter values of some preliminary simulations, regarding the clever section of noisy and the follower, noise-free merit functions. The one-dimensional clever sections were approximately determined by a projection algorithm, using real-life merit function values in the points of a coordinate grid. O nce the minimum of the real-life merit function was determined, the same algorithm was reused to determine the one-dimensional clever sections of the follower merit functions, too. These were used to determine the error domain and, to observe uniqueness under the effect of added non-linearities.
Three successive iteration steps, including the conjugate gradient starting and ending points and coarse grids centered the ending critical points determined by the conjugate gradient method are shown in Figure 7. The section along the fold of the projection within the local grids were represented, revealing the complex geometry of the merit function (Fig. 7). The success of the local grid as apparent but no reliability information was found for the solution and the computational work was excessive.
The projection algorithm was used with linear sub-minimisation according to Fig. 3. In case of model-version H a 4-dimensional coordinate grid was used, the four one-dimensional clever sections were determined by projection and a subsequent minimization for of the real-life and the follower merit functions. In case of model-version HC, 1-dimensional grid was used to determine the clever section of the real-life and the follower merit functions for parameter c. The clever clever section of parameter c determined for the two model-versions were similar, indicating two minima both models (Fig. 8), one minimum was related to a nice, physically admissible initial condition, another to a physically non-admissible one. The compression curve points ( σ ) were determined by model HC from each stage. In case of short stages, the good minimum was quasi-degenerated (Fig. 9). Evaluating the short stages with the simplest model HC, the implicit function g ( c ) was determined for each stage only. From function g ( c ), by using the c identified from the last, long stage, the global minimiser vector and parameters σ were selected. The latter gave compression curve points, being situated above the one determined by the long compression test in accordance with the expectations following from the time dependency of the constitutive law (Fig. 9).
Fig. 7.
The complexitiy of the topography, implicit functions at three successive steps, fine and coarse local grids. Conjugate gradient method (with arrows) ([13 to 16]).
Fig. 8.
Model H and HC, long stage, clever sections c. Note the non-unique solution for both models ([14]).
Fig. 9.
Model HC, short stages 1 to 11 and last, long stage 12. The clever sections c is quasi-degenerated for stages 1 to 11. The c identified from the long stage 12. was used for the evaluation of the short stages ([14]).
10 100 1000 10000 100000
TIME, t (s) ( t ) ( k P a ) Stage 6computedmeasured (a) Total stress.
10 100 1000 10000 100000
TIME, t (s) ' ( t, H ) ( k P a ) Stage 6 (b) Effective stress.
10 100 1000 10000 100000
TIME, t (s) u ( t, H ) ( k P a ) Stage 6 (c) Pore water pressure u(0,y) /u(0,H) H e i gh t o f s a m p l e , y ( c m ) _ (d) Identified initial condition in the various stage. Fig.10.
MRT, fitted and measured stresses, identified initial condition in the various stages.
Table 4.
Parameters of the most general oedometric model-version H Model Parameter Symbol Dependence Consolidation Initial condition A, B , C linear Coeff. of consolidation asymptotic total stress c σ nonlinear linear Relaxation coefficient of relaxation pause of relaxation delay time s t t nonlinear nonlinear nonlinear able 5. Parameters of the consolidation oedometric relaxation test model-version HC Model Parameter Symbol Dependence initial condition A, B , C linear Consolidation Coeff. of consolidation asymptotic total stress c σ nonlinear linear Table 6. Number of function value evaluations, oedometric relaxation test, model-version HC and H Model-version H HC Parameter ' [kPa] v [ mm ] shortrelaxation teststandardcompression test Fig. 11.
Compression curves measured by standard compression tests and identified from short stage MRT data, using model HC. The c identified from the last stage was used for the evaluation of the short stages with quasi-degenerate minimum (see Fig. 9). Discussion
The results can be summarized and discussed as follows.
The non-linear parameter identification related to a strictly convex, analytic Least Squares merit function can be split in terms of dimension. In the first part of the split inverse problem, a section of the merit function ( F [ g ( x ), x] , called clever section) is determined by sub-minimisation, then the clever section is minimised. Two different minimisation methods can be used at two parts of the direct sum decomposition of the parameter domain. The method can be applied repetitively. Geometrically, the fold of a map is determined, the map is the orthogonal projection of the graph of the merit function onto a coordinate plane < F , x >, where p = ( x , y ) is a direct sum decomposition). Analytically, the multi-dimensional implicit function theorem is applied for the y - derivative of a projection map related the graph of a merit function onto a vertical coordinate plane, generated by x , where p = ( x,y ) is a direct sum decomposition of the parameter space. The implicit function g ( x ) of the derivative of the projection map with respect to y can be determined by sub-minimization at every fixed x since the second derivative of the projection map with respect to y is positive definite. The implicit functions - for strictly convex merit functions, related to various direct sum decompositions - are globally unique and have the same lattice structure as the subspaces of the direct sum decomposition. The method can be used not only for strictly convex, analytic Least Squares merit function but also in the general case on condition that the Hessian of the merit function with respect to y is globally positive definite. The most important practical application is the elimination of the linear part of parameter vector through sub-minimisation, which can be added to any non-linear minimisation algorithm. The linearly dependent parameters can be eliminated by sub-minimisation before each function value evaluation.
An important practical application is the case of the quasi-degenerate minimum of convex merit functions. By using the known parameter vector value x min of the global minimum, y min = g ( x min ). The implicit function of convex merit functions is injective at the global minimum, according to the hypothesis and the experiences. The one-dimensional clever sections can approximately be determined by a projection algorithm suggested here. The one-dimensional clever section F [ g ( x i ), x i ] of the merit function related a parameter x i can be used to determine the sensitivity of the parameter and the diameter of the error domain, moreover, to test uniqueness. The real-life merit functions have an unimaginable complexity of an M -dimensional topography, where M is the number of the parameters. However, there is an underlying noise-free merit function, called follower merit function which can be constructed from simulated, noise-free data using the solution of a Least Squares minimization. The follower merit function is not known a priori but can be constructed once the global minimiser is known. In this way it can be used in reliability testing, can be used to define a parameter error domain and to define uniqueness criteria. The follower merit function may have much less critical points, can be even convex. Concerning convexity, the following comment can be made. There is a false opinion that no bracketing procedure can be used in multidimensional space [1]. However, there is an ideal exception which is the class of the is strictly convex minimisation which may occur in terms of noise-free merit functions. Hyperplanes generated y coordinate a mesh points can bracket decreasing convex level surfaces of these.
Examples
The most important practical application of the hierarchical minimisation is the elimination of the linear part of parameter vector through sub-minimisation, which can be added to any non-linear minimisation algorithm. The linearly dependent parameters can be eliminated by sub-minimisation before each function value evaluation. This idea was represented in case of the secant method in Figure 5(b). A note has to be added here. Being the secant method such a non-linear minimisation method that is inherently linear in each step, it may occur that a linear model part is automatically eliminated without the inclusion of a separate, linear algorithm for the elimination. Further research is needed on this. The secant method has two further advantages, relatively small difference between the follower, noise-free and the real-life, noisy merit functions is controlled by the value of the follower noise-free merit function, which is a pseudo-metric, (expressing the distance from the global minimum, where its value is zero).
In the secant method, the step size is linearly dependent on the merit function value. There are several consequences. The step size is comparable with the size of the irregularities of the real-life merit function with respect to the follower merit function on the similarity domain. If the function value at x is large, far from the global minimum, the stepsize is large and vice versa; the global minimum is approached effectively. the secant method In the case of strictly convex, follower merit functions, the bracketing of the level surfaces is theoretically possible on the M -dimensional topography by the hyperplanes of a properly selected coordinate grid (generated in the parameter space). In case of real-life inverse problems, only the follower, noise-free merit function can be strictly convex due to the noise and, the convexity can be used on the similarity domain. According to the experiences, this may be the case in case of the injective solutions of linear system of PDE-s, in case of fixed, convex initial condition and sufficiently long data series. In case of the injective solutions of linear, coupled point-symmetric consolidation models (Table 3), the follower merit function can be strictly convex. Some non-iterative, non-linear parameter identification methods were elaborated related to these injective solutions, by using the suggested projection algorithm and global coordinate grid. In the second example, an important practical application was illustrated on the treatment of the quasi-degenerate minimum in case of “short data series”. Conclusion
The following numerically-efficient strategies were suggested for non-linear inverse problem solution and related reliability testing. o The concept of the follower noise-free merit function (i.e. the best, noise-free approximation of the real-life merit-function) was introduced. Having less minima than the real-life merit-function, it can be used to skip the local minima of the real-life merit function. o The follower merit function is not known a priori but its decreasing direction can be assessed from the real-life merit-function on a part of the parameter domain if the step size is large enough. This can be realized by adjustable step size (e.g. being dependent on the actual fitting error). There are some „nice” models where the shape of the follower merit function is globally convex, e.g. the classical solutions of a system of linear partial differential equations. In this case a mesh can be used to bracket the vicinity of the global minimum. o The concept of a hierarchical inverse problem was introduced which can be used to decrease the parameter number in a non-linear minimization by eliminating the linear part of parameter vector through sub-minimisation which result in a deepest section of the Least Squares merit function. o The hierarchical inverse problem is equivalent to original one if the noise-free follower Least Squares merit function is convex. In this case the parameter error domain is connected and the error can be characterized by its diameter. The transformed form of the follower merit function where all parameters are eliminated except one are the best for determine the diameter of the error domain. The complementary methods were illustrated by the example of the modification of the shark secant method (with adjustable steps) and by the evaluation of a basic geotechnical test, the one dimensional oedometric compression test (with pre-determined coordinate mesh). Further research is suggested.
10. Suggestion for further research
Further research is suggested, among others • on the convexity of the follower merit functions related to injective solutions of system of coupled, linear partial differential equations; including the effect of non-convex initial conditions and an added non-linearity, on the injectivity of the implicit function; on the difference between real-life and convex follower merit functions, by using one-dimensional clever sections; • on the differences between the various (statistical and geometrical) parameter error domains in the light of the true error sizes; • on such step size of the non-linear, multidimensiona minimization methods that allows the use the follower merit function in minimisation; • on the application of the hierarchical minimisation of Least Squares merit function in case of system of non-linear equations and in the case of general non-linear optimization where the merit functions may have a slightly different type, for example it can not necessarily be expressed in the form of a Least Squares merit function defined in section 1; • on the secant – type non-linear parameter identification methods, the performance in case of linear models parts of general non-linear merit function and Least Squares merit functions. References [1] Press, W.H.; Flannery, B.P.; TeukoLeast Squares ky, S.A.; Wetterling, W.T. (1986): Numerical Recipes. Cambridge Univ. Press. [2]
Meier J, Rudolph, S,. Schanz T: Effective Algorithm for parameter back calculation - Geotechnical Applications published in: Bautechnik Volume 86 Issue S1, 86 - 97 [3]
Milnor, J.W. (1963) Morse theory. Princeton Universit Press. 153 p. [5]
Szőkefalvi-Nagy, Gy; Gehér, L.; Nagy, P. (1979) Differentialgeometry. Műszaki Könyvkiadó, Budapest. 256 p. [6]
Wolfe, P. (1959). The secant method for simultaneous nonlinear equations.
Comm. of ACM , 2, pp. 12-13. [7]
Berzi P, Imre E (1998). Parameter identification by the amended secant method. In: Delaunay, D; Jarny, Y; Woodbury, K A (szerk.) Proceedings of the 2nd International Conference on Inverse Problems in Engineering Fairfield (NJ), USA: American Society of Mechanical Engineers (ASME), (1998) submitted. [8]
Berzi P (2021) Root-finding of a nonlinear function with an updated secant method. Int.Conf. on opt. and Appl. submitted. [9]
Imre, E., Vijay P. Singh and Fityus S. (2013). The modelling of some point-symmetric tests 166-185.
Proc. of the 3rd Kézdi Conference. Budapest, Hungary , 2013.05.28.ISBN 978-963-313-081-0. [10]
Imre, E. (1995). Model discrimination for the oedometric relaxation test.
Proc. of 8-th Baltic Geotechnical Conference , 55-60. [11]
Imre, E. (1995). Model discrimination for conventional step-loaded oedometric test. Proc. of the Int. Symp. on Compression and Cons. of Clayey Soils, IS-Hiroshima'95, 525-530. [12]
Imre, E.; Rózsa, P.; Bates, L.; Fityus, S (2010) Evaluation of monotonic and nonmonotonic dissipation test results.
Computers and Geotechnics , 37, Issues 7-8, 885-904. DOI: 10.1016/j.compgeo.2010.07.008 [13]
Imre E, Schanz T, Hegedűs Cs (2013) Some thoughts in non-linear inverse problem solution. EURO:TUN . [14]
Imre, E (1998) Inverse problem solution with a geometrical method 1998 In: Delaunay, D; Jarny, Y; Woodbury, K A (szerk.) Proceedings of the 2nd International Conference on Inverse Problems in Engineering Fairfield (NJ), USA: American Society of Mechanical Engineers (ASME), (1998) pp. 331-338. , 8 p. [15]
Imre, E; Hegedüs Cs; Kovács S. (2018) Some Comments on the Non-Linear Model Fitting IEEE 18th International Symposium on Computational Intelligence and Informatics (CINTI 018) Budapest, Hungary [16]
Imre, E Trang, Q.P Fityus, I. Telekes, G (2011) A geometric parameter error estimation method. ComGeo II.405-414. [17]
Advanced Engineering Evaluation Software CPTu dissipation test, Imre-Rózsa model 2D [18]
Advanced Engineering Evaluation Software CPTu dissipation test, Randolph-Wroth model 2D [19]
Advanced Engineering Evaluation Software CPTu dissipation test, Imre-Rózsa model 3D [20]
Advanced Engineering Evaluation Software Compression test (MCT, modified Bjerrum model with and without the identification of the initial condition, AC, BC) models [21]
Advanced Engineering Evaluation Software Compression test (MCT), (modified) Terzaghi model(A,B) models [22]
Advanced Engineering Evaluation Software Oedometric relaxation test (MRT), Imre-models (H, HC, HCR, HCRT)
Appendix A Numerical xample for the introduced definitions at Least Squares inverse problem
Let us assume the following model with a single linearly dependent parameter a : xatx)u(t , = (43) where t is the time variable x is the space variable (fixed in here), a is the single parameter. Simulated data vector assuming === atx (i.e. fixed values for the sampling times, the space variable and parameter a ): = f (44) Noisy data vector and noise vector: ,2.0 ,2.0,1.0 ,2.188.7 1.2 −= = zf m (45) The direct formulation of the noisy inverse problem: m fAa = (46) The direct formulation of the noise-free problem: fAa = (47) where matrix A is as follows (with elements being equal to the derivative of the model with respect to a at fixed value of the sampling time vector and the space variable): = A (48) In detail, the direct formulation of the noise-free inverse problem: = a (49) In detail, the direct formulation of the noisy inverse problem = a (50) The solution can be determined by the Gauss Normal Equations: mTT fAaAA = ' (51) Solution with the Gauss Normal Equations for the noise-free problem: fAaAA TT = (52) In detail: === afAAA TT (53) Solution with the Gauss Normal Equation for the noisy problem: === afAAA mTT (54) Simulated data vector assuming === atx (i.e. fixed values for the sampling times, the space variable and parameter a ): = f (55) Real-life merit function: )92.18()48.7()1.2()( aaaaF −+−+−= (56) The simulated, noise-free merit function related to then true value of the parameter a : )918()48()2()(' aaaaF −+−+−= (57) The follower, simulated noise-free merit function, related to the solution of the noisy inverse problem, the identified value of parameter a : )9101016.18( )4044896.8()2,011224()(' a aaaF − +−+−= (58) The clever sections of the noisy and the follower merit function for parameter a are shown in Fig. 2. The true error (~0,011) is equal to about 1/5 of the width of the error domain (~0,056) or about 4/10 of the half width of the error domain called generalized standard deviation. a [-] F [ - ] Fig. 12. The clever sections of the noisy and the follower merit function for parameter a Appendix B The implicit function theorem
Let f : R n + m → R m be a continuously differentiable function. R n + m is the direct sum R n + R m , a point of this is ( x , y ) = ( x , ..., x n , y , ..., y m ). Starting from the given function f , the goal is to construct a function g : R n → R m whose graph ( x , g ( x )) is precisely the set of all ( x , y ) such that f ( x, y ) = . Fix a point ( a , b ) = ( a , ..., a n , b , ..., b m ) with f ( a , b ) = , where is the element of R m . If the matrix [(∂ f i /∂ y j ) ( a , b )] is invertible, then there exists an open set U containing a , an open set V containing b , and a unique continuously differentiable function g : U → V such that == )(,),()(, UU (1) Regularity
It can be proven that whenever we have the additional hypothesis that f is continuously differentiable up to k times inside U × V , then the same holds true for the implicit function g inside U and )}(,)}(,)( xgxxfxgxyfxxg −= − (2) imilarly, if f is analytic inside U × V , then the same holds for the implicit function g inside U . This generalization is called the analytic implicit function theorem. Appendix C Modelling details, example 2 C.1 The consolidation model
In this section the linear coupled consolidation models elaborated for the relaxation/compression test models (m = 1) are presented, denoted by ORT and OCT models here. System of differential equations
Two equations have been derived from the equilibrium condition and, from the continuity condition, assuming weightless soil and water. The Equilibrium Equation compiles the equilibrium condition, the effective stress equality, the geometrical and, the constitutive equations, as follows: oed (59) and, the Continuity Equation compiles the continuity equation, the Darcys law and the geometrical equation, as follows: v (60) where v is displacement, u is pore water pressure (neglecting the gravitational component of the hydraulic head), y and t are space and time co-ordinates respectively, Eoed is the oedometric modulus. The volumetric strain and the Laplacian operator, for dimension m = 1 are as follows: yv=ε , y=Δ (61) Boundary conditions
The displacement domain of the oedometric test is bounded by two zero dimensional circles, with radii r = 0 and r = H, respectively. For a double-drained sample, the upper-lower surface coincides with the outer boundary, respectively, the symmetry point coincides with the inner boundary. In the following four boundary conditions are presented . Three are common, one is different for the two models. (1) The (common) boundary condition Nr. 1 : y= (62) (2) The (common) boundary condition Nr. 2 : Hy= (63) (3) The (common) boundary condition Nr. 3 : Hy = (64) (4) Boundary condition Nr. 4 concerning the coupled 1 models: y= (65) (5) Boundary condition Nr. 5 concerning the compression test models: . y= (66) Analytical solution
Steady-state solution part The solution for the relaxation test model: − = Hyy vv p (67) Hvy E oedp )( −= (68) and, for the compression test model: ( ) .)( yH y Ev oedp −= (69) )( = y p (70) Transient solution part The transient part of the displacement v t for the oedometric relaxation test model: eyHka = y)(t,v tcH k-k=kt v222 πsin (71) where ak ( k = 1... ) are the Fourier coefficients of an odd initial displacement function and, for the compression test model: eyHkb = y)(t,v tcH4 k-k=kt v2 2 − )12(1 (72) where bk ( k = 1... ) are the Fourier coefficients of an even initial displacement function. The function u for the relaxation test model: eyHk = y)u(t, tcvH22 k2-=k − k1 (73) and, for the compression test model: e yH 2k = y)u(t, tcH4 k-1=k v2 2 − − )12(k )12(sinβ (74) where ;)12(β; k H 2kb= Hka= kkk −− k10 = − (75) The following two dimensionless arguments R and T can be derived for the compression and for the relaxation test model, respectively: HyY = or HyY = (76) HctT = or HctT = (77) It follows from the difference of the time factors that the model law is different for the two tests. It also follows that slower dissipation time can be expected for the OCT than for the ORT in case of identical initial conditions. C.2 Qualitative behaviour
Stress state in the sample
For the ORT model, the stress variation within the sample is partly rebound, partly compression as follows. The total stress is equal to the sum of a steady-state and a transient component. The former - the final total stress - depends on the displacement load v and the H the sample height: vE oed = (78) The latter is equal to the mean pore water pressure u mean : dy y) ,u(t H1tut H0meant == )()( (79) Effective stress decrease takes place at the sample top ( y = 0), and compression (increase in the effective stress ) takes place in the bottom ( H = 0) of the sample ),()(1),( ytutuEyt meanoedt −−= (80) For the OCT model, the total normal stress is constant with time. The effective normal stress increases with time, being its transient part equal to the pore water pressure. The consolidation model predicts basically compression: ),(1),( ytuEyt oedt = (81) C.3 The model parameters
The analytical solution of the linear, point-symmetric coupled consolidation model of the oedometric relaxation test for the so-called total stress and the so-called pore water pressure u (Appendix C, [8 to 11]): + e - = t T i-i1=i )( (82) e H iy - = yt,u T 2 i2-i1=i ]1)/([cos()( − (83) where T = ct / H is the time factor, c is the coefficient of consolidation, i are Fourier coefficients related to the initial condition, H is the half-width of a double-drained sample (length between the boundaries of the model). The time dependent constitutive law was included in an approximate form into the model in its most complex form, by adding a relaxation term to the total stress. The relaxation term includes a partial unloading effect. The equations used for model fitting were as follows. The equations used for model fitting were as follows. )()( t - + e - = t r T 2 i2-i1=i (84) The empirical relaxation term: t>t;t+t t+tsb-1 1(0) s= (t)
331 1r log (85) where s is the coefficient of relaxation, t is the pause of relaxation in case of partial unloading, b = log(( t + t )/ t ), t is time parameter. The analytical solution of model-version H contains 8 parameters, 4 out them are with non-linear dependence (Table 4). In model-version HCRT parameter t was specified, b was zero and, a minor change was made in the way of the identification of the relaxation parameter s as follows. The product of s and σ (0) denoted by sk (0) s= s k (86) was identified instead of s . The solution depended linearly on s k . As a result, the number of the non-linearly dependent parameters was decreased to 2 and, the computational work was less by more than two orders of magnitude than the one for model H. In model-version HCR value of t , t were specified, the modified relaxation parameter ( s k) was identified instead of s using the facts that term σ (0) can be expressed as the linear combination of other parameters: D + = (0) (87) where D is the mean initial pore water pressure. In model-version HCR only one among the three parameters of the relaxation part-model was identified ( s k ), value of t was set to zero. The relaxation branch of the model was left out in model-version HC (Table 5). The number of the non-linearly dependent parameters was 1 for both model-version HCR and HC. The computational work was less by more than three orders of magnitude than the one for model H (Table 6).was set to zero. The relaxation branch of the model was left out in model-version HC (Table 5). The number of the non-linearly dependent parameters was 1 for both model-version HCR and HC. The computational work was less by more than three orders of magnitude than the one for model H (Table 6).