Modification of the Quasilinearization Method for the Inverse Problem

We propose a new modification of Bellman’s quasilinearization method such that at any iteration step, it works with an approximate solution of the original nonlinear system and with new approximation of parameters α(k+1) which are close enough to the previous ones. As an output, this approach provides a construction of a convergent sequence of parameters where the limit is the best approximation of parameters of a given system. We apply this method to a mathematical model describing BSP-kinetics in the human liver.


Introduction
For solving the inverse problems, in particular, for identification of systems with known structure, the quasilinearization method (QM) is a standard tool.Designed by Bellman et al. [1], this method was later applied to different kinds of identification problems (cf.[2] or [3] for references).We were interested in application of QM to solve the parameter identification problem for the BSP-kinetics in the human liver [4][5][6][7].One of the possible descriptions of this kinetics can be given by the nonlinear system of ordinary differential equations where X(t), Y (t), Z(t) mean the amount of BSP in the blood, in the membranes of hepatic cells, inside the cells at the time t, respectively, and α = (c 1 ,c 2 ,c 3 ,K 1 ,K 2 ) is a vector of unknown positive parameters [6].Suppose a "single injection" in which the amount I (mg) of BSP is injected into the blood at once.This leads to the initial conditions In order to uniquely determine the unknown positive parameters α = (K 1 ,K 2 ,c 1 ,c 2 ,c 3 ) , we have to know at least two different data sets.From practical point of view, we can obtain data describing the decreasing level of BSP in the blood (Table 1.1) and in Table 1.2, they are presenting the measurements of BSP in the bile.These data were obtained through medical experiments by Hrnčíř [6].
The first data set corresponds to the function X(t).The second one corresponds to the function V (t) = I − X(t) − Y (t) − Z(t) describing the level of BSP in the bile.
However, the standard approach like in [2,3], or recent [8,9] does not provide the reasonable outputs corresponding to the nature of parameters, especially if we solve an identification problem for nonlinear system of ordinary differential equations.(We can obtain negative values of determinated parameters, see Section 5.) Therefore we propose a modification of the quasilinearization method (MQM).The algorithm of the modified QM consists of the steps displayed below.Let us briefly introduce the MQM (see Section 3 for details).
The classical approach used by Bellman (see [2,3]) is similar to Algorithm 1.1 with the exception of Step 3 (which requires the computation of the solution of the given differential equation in every step of the algorithm) and with the exception of Steps 6 and 7.In the existing sources, like [2,3,8,9], only the linearized differential equation given in Step 4 is used only.This makes things easier from the viewpoint of computation and works properly especially for linear systems of differential equations.The development of computing devices since the eighties of the last century and the software (like the package Mathematica) allow to do the computations fast even if the given differential equation is solved approximately in every step of determining a better approximation of the values of parameters.The problem is that the solution of the differential equation (1.1) for the certain value of the parameter can be far from the solution of this equation linearized around the fixed solution x (k) from Step 3.This obstacle is removed by Steps 6 and 7 especially in the case of nonlinear differential equations.In this way, the final value of the Lenka Čelechovská-Kozáková 3 Step 1.Consider a nonlinear autonomous initial problem ẋ = f (x,α), x(0) = c, where x ∈ R n , α ∈ R N , and f : R n+N → R n is a continuous function.This problem is equivalent to the Cauchy problem ẋ = g(x), x(0) = c, where Step 2. Choose the initial approximation α (1) , the tolerance ε > 0, and put k = 1.
Step 9.If S(x (k+1) ) > S(x (k) ), then go back to Step 2 and start the algorithm with a better choice α (1) .
Algorithm 1.1 parameters is reached (according to the criteria for stopping the computatuion given in Steps 8-9).The organization of this paper is as follows.In Section 2 we give a basic notations and definitions.In Section 3 we describe the modification of quasilinearization method in detail, and in Section 4 we give the convergence theorem.Section 5 includes the numerical results.

Notations and definitions
Let R m be a vector space with the scalar product )

.2)
Let A = (a i j ), i, j = 1,...,m, be an m × m matrix.Then the matrix norm is given by The matrix A is called positive definite if there is a constant K > 0 such that (u,Au) ≥ K u 2 (2.4) Let M be an m × m symmetric matrix of the form where m, and E is the m × m identity matrix.Then the matrix M is positive definite. Proof.Denote We can write the matrix M kk in the form Lenka Čelechovská-Kozáková 5 where e i = (0,...,0,1,0,...,0) is the k-dimensional vector with 1 on the ith position, i = 1,...,k.The minor detM kk of the matrix M can be evaluated as follows: det det e 1 ,...,e l−1 ,Γ l ,e l+1 ,...,e k + where Q j are the matrices with at least two columns Γ r ,Γ s .These k-dimensional vectors Γ r ,Γ s are not linearly independent since for all i = 1,...,k.Therefore, and the matrix M is positive definite by Sylvester criterion [10, page 248].
is positive definite too.
The proof is clear. (2.12) Let the matrix M d have the form (2.11).Then is a scalar product on L2 m [0,T] too.The proof follows easily by Lemma 2.2.

T] be the normed space of continuous m-dimensional vector functions with the norm
(2.16) where the norm h is defined by (2.14).
Proof.We can write (2.18) From this inequality, the assertion of Lemma 2.5 follows.
Let D ⊂ R m be a convex set.The function S : D → R is called a strictly convex function if there is a constant χ > 0 such that for every u,v ∈ D and for every α ∈ [0,1], the inequality is satisfied.The constant χ is called the constant of the strict convexity of the function S on the set D.
Lemma 2.6.Let D ⊂ R m be a convex closed set.Let S(u) have the form

.21)
where A is a positive definite m × m matrix, b ∈ R m , and c ∈ R. Then S is a strictly convex function.
The proof is clear.

Modification of the quasilinearization method
Let Q ⊂ R n be a closed convex set of the variables x = (x 1 ,...,x n ) and let D ⊂ R n be a closed convex set of the parameters α = (α 1 ,...,α N ) .Let f : Q × D → R n have continuous bounded partial derivatives up to the second order.Consider a nonlinear autonomous system of ordinary differential equations with the initial condition In order to avoid considering two different types of vectors, we will suppose that the vector α satisfies the differential equation with the initial condition where β = (β 1 ,...,β N ) .Define a new vector x by and a vector c (corresponding to the initial condition) by The vector x(t) satisfies the nonlinear differential equation where g(x) = ( f (x,α),0,...,0

N
) , with the initial condition The aim is to find the unknown parameters α such that the solution of the initial problem (3.1) fits in some sense with a given tolerance ε > 0 to the measured data or to the continuous function which approximates these data, respectively.
Assume that the approximating fuction r(t) = (r 1 (t),...,r n (t)) corresponding to the measured data is given and let e(t) be an approximating function appropriate to a certain linear combination of the components of the solution of (3.1) which is again measured during the experiment (in our case, r(t) ≈ (X(t),Y (t),Z(t)) , e(t) ≈ V (t)).In this context, let us point out that in practice, the values of r(t) and e(t) are measured in discrete instants of time, {t 1 ,...,t L } and {s 1 ,...,s M }, L,M ∈ N, and the functions r(t), e(t) have to be produced from given measured values.The procedure how to do this is in fact a matter of taste and intuition.It seems to be reasonable to get the functions r(t) and e(t) using spline interpolation.Our motivation is the Cauchy problem given by (1.1), (1.2) described in the intoduction.
The weighted deviation, Γ : C n [0,T] → R, of a given function z(t) ∈ C n [0,T] from the approximating functions r(t) and e(t) can be expressed, in sense of the least-square method, in the form where γ,γ l are given real weighting constants (in our case, γ = X(0) = I and γ l = −1 for l = 1,2,3).
Proof.Matrix A k+1 is the Gramm matrix which is real and symmetric.Since the vectors h ( j,k+1) (t) are linearly independent, we have detA k+1 = 0. Let λ j , j = 1,...,N, be the eigenvalue of the matrix A k+1 and let u ( j) be the corresponding eigenvector, u ( j) = 0. Then λ j ∈ R and 0 < u ( j) ,u ( j) = u ( j) A k+1 u ( j) = u ( j) λ j u This inequality implies that all eigenvalues are positive.There are orthogonal matrix O k+1 and diagonal matrix D k+1 = diag(λ 1 ,...,λ N ) so that In the next lemma, we give a set and its property in which we look for the minimum of the function (3.24).Lemma 3.4.Let S k+1 (β) have the form (3.24).Denote V k := S(x (k) ), where x (k) is a solution of (3.1) Then M αk is a convex set for all k = 1,2,....
The last inequality holds since A is positive definite.
The necessary conditions for determining the local extreme on the set M αk are given by the equations Let us denote the solution of (3.36) by β * = (β * 1 ,...,β * N ) .Since the matrix A k+1 is positive definite by Lemma 3.3 and the function Ψ k+1 (β) is the strictly convex function by Lemma 2.6, β * is the unique point of minimum (see [11, page 186]).Put In this way, we obtain new initial condition for the solution x (k+1) (t) of (3.6).Computing this solution, we get the solution x (k+1) of the equivalent system (3.1) for α = α (k+1) .Determine the deviation (3.9).If the inequality (3.10), that is, holds and the distance between α (k) and α (k+1) is small, that is, for a given ζ k small, then we can repeat the whole process of enumeration until the condition where ε > 0 is a given tolerance, is satisfied.If the inequality (3.10) is fulfilled, but Lenka Čelechovská-Kozáková 13 we have to modify the value of the parameter α (k+1) .The modification is based on the following lemma.
Lemma 3.5.Let M αk have the form (3.34) (cf.Lemma 3.4).Then for arbitrary ζ k > 0, there is a parameter α (k+1) ∈ M αk such that Proof.Let β * ∈ M αk be an argument of minima of Ψ k+1 (β).Since M αk is a convex set, we can look for the parameter α (k+1) in the form where a ∈ (0,1).The object is to find a proper value a such that the vector α (k+1) has to satisfy the inequality (3.43).We would like to have Hence, we have to choose a such that a ≤ ζ k / β * − α (k) .
We are able to shift the parameter α (k+1) to α (k) such that the distance between α (k+1) and α (k) is arbitrarily small, in particular less than a given tolerance ζ k .
If S(x (k+1) ) > S(x (k) ) (the value of deviation has increased), we have to stop the whole process of computation and to start with a better choice of the initial approximation α (1) .
If S(x (k+1) ) = S(x (k) ) holds, we get the required values of parameters α = α (k) and the algorithm cannot produce better parameter values (for a given α (1) ) and we are finished.

Convergence of the method
We did not manage to formulate the sufficient conditions for convergence of the sequence k=1 generated by the modified quasilinearization method (MQM) for arbitrary initial approximation α (1) .Nevertheless, the method, if it is successful, constructs a convergent sequence of parameters {α (k) } ∞ k=1 .We can choose a sequence {ζ k } ∞ k=1 such that it is decreasing, liminf ζ k = 0, and in addition k=1 be a sequence generated by MQM.Then {α (k) } ∞ k=1 is a Cauchy sequence.Proof.The sum ∞ k=1 ζ k is a convergent sum which consist of positive real numbers, therefore for every ε > 0, there is k From the facts above, it follows that for every ε > 0, there is natural number k 0 such that for every natural number p and for every k ≥ k 0 , the inequality is true.This means that the sequence {α (k) } ∞ k=1 is a Cauchy sequence.Corollary 4.2.The sequence {α (k) } ∞ k=1 has a limit α (∞) .The ideal situation is a construction of the sequence α (k) → α ( * ) such that S(x ( * ) ) = 0, where x ( * ) is a solution of (3.1) for α = α ( * ) .From practical point of view, this ideal situation is very rare, consequently we take up with a sequence for which the condition (3.41) is satisfied.Using MQM, we receive the best possible approximation α (∞) depending on an initial choice α (1) .
Our modification was proved on the simple linear mathematical model of the human liver published in [7].The advantage of the system describing the simple mathematical model is a knowledge of the exact analytic solution.Modification of the quasilinearization method applied to this simple linear model provides identical results as classical Bellman's quasilinearization method for the inverse problem.

. 47 )
Proof.The proposition follows from the continuous dependence of the solution x(t) of (3.6) on the initial conditions[12, page 94].

Table 1 .
1.The amount of BSP in the blood.

Table 1 .
2. The amount of BSP in the bile.