# On an Identification Problem on the Determination of the Parameters of the Dynamic System

Authors: Dossan Baigereyev; Nevazi Ismailov; Yusif Gasimov; Atif Namazov
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

An inverse problem is considered for the determination of the parameters, involved in the right-hand side of the system of nonlinear ordinary differential equations by given initial and final conditions. The solution of the problem is reduced to the minimization of the quadratic functional, which indeed is a deviation of the value of the solution from the given values at the end points. Using the quasilinearization method a calculation method is proposed to the solution of the considered problem. The application of this method is demonstrated on the example of the determination  of the hydraulic resistance in the tubes.


---

## Body

## 1. Introduction

As is known, different classes of inverse and identification problems play an important role in the solition of many applied problems from physics, hydrodynamics, and industry [1–5]. There exist various methods to solve these problems, as well as optimization methods [6–10]. One of the principle steps of these methods is a choice of suitable functional. Since many applied problems are described by the nonlinear systems, the choice of such functional and further solution of the corresponding optimization problem is problematic [2]. These difficulties may be avoid, for example, by the iterational quasilinearization method, convergence of which is in detail studied in [11, 12].In the present work a multidimensional identification problem is considered for the determination of the parameters involving the right-hand side of the system of nonlinear differential equations by given initial and final conditions. Solution of the problem is reduced to the optimization problem, in which the functional under minimization is constructed as a quadratic deviation of the solition of the system from the given data at the end points. Since solution of the problem in the stated nonlinear formulation presents certain difficulties [13, 14], the considered problem is reduced to the linear case with respect to phase coordinates and vector of parameters by the help of quasilinearization method. A quadratic functional is constructed and an expression for its gradient is derived. Using Gram-Schmidt orthogonalization method a calculation algorithm is proposed, which allows one to define sought parameters. This algorithm is applied to the example, describing the flow in the pipes.


## 2. Problem Statement

Let the movement of the object be described by the system of differential equations:(1)y˙t=fyt,α,where y is n-dimensional phase vector, f is n-dimensional differentiable function continuous in the interval (0,T), and α is m-dimensional constant vector to be found.Let the following initial conditions be given:(2)yi0=y0i,i=1,N¯,where Nn, m are given natural numbers and t is an independent variable y0i, i=1,N¯ given n-dimensional vector. The problem consists in the finding of the vector α, by which the solution of the Cauchy problem (1)-(2) in the point T satisfies the given condition(3)yiT=yTi,i=1,N¯.Such problems are often met in applications [1, 3, 6, 15], when initial data (2) are given and final ones are statistically measured. In these cases it is required to find the vector α such that the solution of the problem by initial data (2) is maximally close to the measured data at the end points. As an example the problem in oil-gas production may be shown when it needs to define the coefficient of the hydraulic resistance.


## 3. Solution Method

Since the functionf(y,α) is nonlinear, to solve the problem (1)–(3) it is expedient to use any numerical method, as well as quasilinearization method [11, 13]. So in the first step we linearize the problem (1)–(3). For this purpose some nominal trajectory y0t and parameter α0 are chosen and it is assumed that kth iteration is already hold. If we linearize (1) with respect to these data we obtain (4)y˙k=∂fyk-1,αk-1∂yyk+∂fyk-1,αk-1∂ααk+fyk-1,αk-1-∂fyk-1,αk-1∂yyk-1-∂fyk-1,αk-1∂ααk-1.After integration of the linear differential equation (4) with condition (2) we get the representation [16] (5)ykt=Φk-1t,t0·ykt0+Φ1k-1t,t0·αk+Φ2k-1t,t0,where Φk-1t,t0 is a fundamental matrix of the system of homogeneous equations(6)y˙kt=∂fyk-1,αk-1∂yy=yk-1α=αk-1·ykt0and the matrices Φ1k-1t,t0, Φ2k-1t,t0 are defined as in [16](7)Φ1k-1t,t0=∫t0tΦt,τ∂fyk-1,αk-1∂αdτ,Φ2k-1t,t0=∫t0tΦt,τfyk-1,αk-1-∂fyk-1,αk-1∂yyk-1-∂fyk-1,αk-1∂ααk-1dτ.To provide that the solutionyi(T) of the linearized differential equation (4) with initial data (2) coincides with the values of the measurements yTi=yTi in the point T we construct the following quadratic functional in the kth iteration:(8)Ik=12∑i=1NyikT,α-yTikTAikyikT,α-yTik,where the sign T stands for transpore, Ak is a constant symmetre m×n-dimensional weight matrix, that is chosen in each iteration, considering the specifics of the concrete problem; yTik is n×1-dimensional vector of observation; yik is n×1-dimensional vector defined by (5). Then the solution of the stated problem is reduced to the problem: Find a constant vector α, by which the solution of (1) with conditions (2) minimizes the functional (8).Various algorithms exist for the minimization of the functional (8). However in the solution of the concrete problem, as well as problem arising in the oil production, these algorithms met some difficulties [14] (e.g., to reach the necessary accuracy and speed of convergence). Therefore in [6] the use of Gram-Schmidt orthogonalization method is proposed.


## 4. Algorithm for the Minimization of the Functional (8)

Here we consider the minimization of the functional (8) by the help of the relation (4) with conditions (2). Putting yk(T) from (5) into (8) we get (9)Ik=12∑i=1NΦik-1T,t0yikt0+Φ1ik-1T,t0·α+Φ2ik-1T,t0-yTiTAikΦik-1T,t0yikt0+Φ1ik-1T,t0·α+Φ2ik-1T,t0-yTi=∑j=116Ikj.Considering the symmetricity of the matrix Ak the relation (9) may be written as (10)Ik1≝12∑i=1NΦik-1T,t0yikt0TAikΦik-1T,t0yikt0,Ik2≝12∑i=1NΦik-1T,t0yikt0TAikΦ1ik-1T,t0·α,Ik3≝12∑i=1NΦik-1T,t0yikt0TAikΦ2ik-1T,t0,Ik4≝-12∑i=1NΦik-1T,t0yikt0TAikyTik,Ik5≝12∑i=1NΦ1ik-1T,t0·αTAikΦik-1T,t0yikt0,Ik6≝12∑i=1NΦ1ik-1T,t0·αTAikΦ1ik-1T,t0·α,Ik7≝12∑i=1NΦ1ik-1T,t0·αTAikΦ2ik-1T,t0,Ik8≝-12∑i=1NΦ1ik-1T,t0·αTAikyTik,Ik9≝12∑i=1NΦ2ik-1T,t0TAikΦik-1T,t0yikt0,Ik10≝12∑i=1NΦ2ik-1T,t0TAikΦ1ik-1T,t0·α,Ik11≝12∑i=1NΦ2ik-1T,t0TAikΦ2ik-1T,t0,Ik12≝-12∑i=1NΦ2ik-1T,t0TAikyTik,Ik13≝-12∑i=1NyTikTAikΦik-1T,t0yikt0,Ik14≝-12∑i=1NyTikTAikΦ1ik-1T,t0·α,Ik15≝-12∑i=1NyTikTAikΦ2ik-1T,t0,Ik16≝-12∑i=1NyTikTAikyTikand gradient of the functional (9) has a form(11)∂Ik∂α=∑j=116∂Ikj∂α.Since the terms  Ik1, Ik3, Ik4, Ik9, Ik11, Ik12, Ik13, Ik15, Ik16 do not depend on the parameter α, we have(12)∂Ik1∂α=∂Ik3∂α=∂Ik4∂α=∂Ik9∂α=∂Ik11∂α=∂Ik12∂α=∂Ik13∂α=∂Ik15∂α=∂Ik16∂α=0.Based on the formulas(13)∂xTa∂x=∂aTx∂x=a,∂xTBx∂x=B+BTxfrom [17, 18], for the gradients of Ik2, Ik5, Ik6, Ik7, Ik8, Ik10, Ik14, we get the formulas(14)Ik2∂α=Φik-1T,t0yikt0TAikΦ1ik-1T,t0T=Φ1ik-1TT,t0AikTΦik-1T,t0yikt0,Ik5∂α=αTΦ1ik-1TT,t0AikΦik-1T,t0yikt0∂α=Φ1ik-1TT,t0AikΦik-1T,t0yikt0,Ik6∂α=αTΦ1ik-1TT,t0AikΦ1ik-1T,t0α∂α=Φ1ik-1TT,t0AikΦ1ik-1T,t0+Φ1ik-1TT,t0AikTΦ1ik-1T,t0α=2Φ1ik-1TT,t0AikΦ1ik-1T,t0α,Ik7∂α=αTΦ1ik-1TT,t0AikΦ2ik-1T,t0∂α=Φ1ik-1TT,t0AikΦ2ik-1T,t0,Ik8∂α=-Φ1ik-1T,t0αTAikyTik∂α=-αTΦ1ik-1TT,t0AikyTik∂α=-Φ1ik-1TT,t0AikyTik,Ik10∂α=Φ1ik-1TT,t0AikTΦ2ik-1T,t0,Ik14∂α=-yTikTAikΦ1ik-1T,t0α∂α=-Φ1ik-1TT,t0AikTyTik.Finally if we consider these results, then the gradient of the functional (9) will be defined by the formula (15)∂Ik∂α=12∑i=1NIk2∂α+Ik5∂α+Ik6∂α+Ik7∂α+Ik8∂α+Ik10∂α+Ik14∂α=12∑i=1NΦ1ik-1TT,t0AikTΦik-1T,t0yikt0+Φ1ik-1TT,t0AikΦik-1T,t0yikt0+2Φ1ik-1TT,t0AikΦ1ik-1T,t0α+Φ1ik-1TT,t0AikΦ2ik-1T,t0-Φ1ik-1TT,t0AikyTik+Φ1ik-1TT,t0AikTΦ2ik-1T,t0-Φ1ik-1TT,t0AikTyTik=12∑i=1N2Φ1ik-1TT,t0AikTΦik-1T,t0yikt0+2Φ1ik-1TT,t0AikΦ1ik-1T,t0α+2Φ1ik-1TT,t0AikΦ2ik-1T,t0-2Φ1ik-1TT,t0AikyTik=0.Then, for the gradient of the functionalIk relatively to the parameter α we get the expression (16)∂Ik∂α=∑i=1NΦ1ik-1TT,t0AikΦ1ik-1T,t0α-∑i=1NΦ1ik-1TT,t0AikyTik-Φ1ik-1TT,t0AikTΦik-1T,t0yikt0-Φ1ik-1TT,t0AikΦ2ik-1T,t0.Taking equal to zero the expression (16) we get(17)∑i=1NΦ1ik-1TT,t0AikΦ1ik-1T,t0α=∑i=1NΦ1ik-1TT,t0AikyTik-Φ1ik-1TT,t0AikTΦik-1T,t0yikt0-Φ1ik-1TT,t0AikΦ2ik-1T,t0.Solution of (17) with respect to α gives(18)α=∑i=1NΦ1ik-1TT,t0AikΦ1ik-1T,t0-1∑i=1NΦ1ik-1TT,t0AikyTik-Φ1ik-1TT,t0AikTΦik-1T,t0yikt0-Φ1ik-1TT,t0AikΦ2ik-1T,t0,where it is assumed that ∑i=1NΦ1ik-1TT,t0AikΦ1ik-1T,t0-1 exists.Value of the parameterα, defined by (18) is a solution of the multiparameter optimization problem for the functional (9) that gives minimum to the cost functional.Considering the above the following algorithm may be proposed to the solution of the identification problem (1), (2).Algorithm 1.
(
1
) Construct the function fx from (1), initial and final data y0iand yTii=1,N from (2) and (3) correspondingly.
(
2
) Calculate the derivatives ∂fyk-1,αk-1/∂y, ∂fyk-1,αk-1/∂α taking as initial approaches yiand αi.
(
3
) Calculate the fundamental matrix Φk-1t,t0 from (6); reconstruct Φ1k-1t,t0, Φ2k-1t,t0 from (7) and functional Ik from (8).
(
4
) Solving the system of algebraic equations (14) relatively α find the value of the m-dimensional vector αk in the kth iteration.
(
5
) Check the condition (19)∂Ik∂α<ε,where ε is a given small enough number. If the condition (16) is satisfied the process stops; otherwise go to Step (2). The convergence of this algorithm may be proved similar to [13, 14].Now we discuss the realization of this algorithm.


## 5. Calculational Algorithm

As one can see, in the realization of the above algorithm the main step is a calculation of the fundamental matrixΦk-1t,t0 and the matrices Φ1k-1t,t0, Φ2k-1t,t0. Note that, as is mentioned in [16], construction of these matrices is an enough difficult procedure. So for simplicity we try (4) and find the corresponding derivatives by using the Euler method. Really ∂fyt,α/∂y and ∂fyt,α/∂α everywhere are replaced by f~y=fy+δ,α-fy,α/δ and f~~α=fy,α+δ1-fy,α/δ1 correspondingly, where δ and δ1 are small enough numbers.To calculate the fundamental matrixΦk-1t,t0 and the matrices Φ1k-1t,t0, Φ2k-1t,t0 it is proper to replace (4) by the following discrete one:(20)ykt2N=∏i1=2N-1jE+δf~~yk-1ti1yktN+1+Φ1k-1t,t0α+Φ2k-1t,t0,where(21)Φik-1t2N,t0=∏i1=2N-1jE+δf~~yk-1ti,Φ1ik-1t2N,t0=∑j=N+22N-1∏i1=2N-1jE+δ·f~~yk-1tiδf~~αj-1k-1+δ·f~~α2N-1k-1,Φ2ik-1t2N,t0=∑j=N+22N-1∏i1=2N-1jE+δ·f~~yk-1tiδ·f~yk-1tj-1,αk-1-f~~yk-1tj-1yk-1tj-1-f~~αk-1αk-1+δ·f~yk-1t2N-1,αk-1-f~~yk-1t2N-1yk-1t2N-1-f~~αk-1αk-1.E is unit matrix of proper dimension.Then from (20) we get that Φk-1t,t0 is a fundamental matrix for the system of homogeneous equations:(22)ykti+1=E+δf~yy=yk-1α=αk-1·ykti.Therefore, similar to the nongradient methods we propose an algorithm, based on the orthogonalization of gradient directions using the Gram-Schmidt procedure.Step 1.
Using Gram-Schmidt orthogonalization, vectorsωij, i=j,j+1,…,n, j=2,3,…,n2 are calculated and the set of vectors(23)∂fλc1∂λc1,∂fλc2∂λc2,…,∂fλci∂λci,ωi+1j,…,ωnjis found which form an orthogonal basis in Rn.
If we apply the orthogonalization algorithm, then a linearly independent systema1,a2,…,ak should be formed, that is, orthogonal system b1,b2,…,bk and each vector bi should be linearly expressed through a1,a2,…,ai. Here ai and bi are upper triangular matrices. Thus it is possible to ensure that the systems bi were orthonormal, where the diagonal elements of the transition matrix are positive; by these conditions the system bi and the transition matrix are uniquely determined.
The algorithm considersb1=a1, if the vectors, b1,b2,…,bi-1, are constructed. Then (24)bi=ai-∑j=1i-1ai,bjbi,bjbj,where  is the sign of the scalar product of vectors.Step 2.
For the orthogonalization of the gradient directions we computeνij in the form (25)νij=fλc+δωij-fλc-δωij2δ,i=j,j+1,…,n.Here δ>0 is any small parameter.Step 3.
The orthogonal gradient directions are chosen in the form(26)lj=∑i=jnνijωij.Replacing ∇f(λc(k)) by lj in (17) the nongradient iterative minimization procedure will be(27)λck+1=λck-χ∗klk,where χ∗(k) is a scalar, which is determined by golden section method.Now we apply the above proposed technique to the example of 15 production by gas-lift method.Example 2.
It is known that nonstable motion of gas in tubes and gas liquid mixture (GLM) in vertical tubes, that is, in the lift pipe of the gas-lift well with constant across profile, is described by the following system of linear partial differential equations of hyperbolic type (see Figure1):(28)-∂P∂x=∂ρωc∂t+2αρωc,-∂P∂t=c2∂ρωc∂x,where ρ=P(x,t), ωc=ωc(x,t), is an additional pressure on its stationary value and averaged over across section speed of motion of GLM; t,x: time and coordinate; c: speed of sound in gas and GLM; ρ: gas, oil, and GLM in correspondence with coordinates; 2α=g/ωc+λcωc/2D; g, λc: free fall accel …; D: interval effective diameter of the tube [19].
The partial differential equation of gas and GLM motion by are averaging over timet may be reduced to the following ordinary differential equation [20](29)Q˙=2αλcρFQ2c2ρ2F2-Q2,Q0=u,where c≫ωc and all quantities are assumed constant Q=ρωcF and F is area of across section of the pump-compressor tubes and is constant relative to axes.Figure 1It is assumed that the passing from the end of tube through the layer to the beginning of the lift (x=l) is described by the following difference equations:(30)Ql+0=γQl-0+γ1Ql-0Q-,γ1Ql-0=-δ3Ql-0-δ22+δ1,where γ and δ1, δ2, δ3 are constant numbers to be found. For the sake of simplicity we suppose that the parameters γ, δ1, δ2, δ3 are known and it is required to reconstruct λc involved in (19) due to α(λc).Then some nominal trajectoryQ0x and parameter α0 are chosen assuming that kth iteration is already held. Let us linearize (29) among these data(31)Q˙kx=AQk-1,αk-1·Qkx+BQk-1,αk-1αk+CQk-1,αk-1,where (32)∂fyk-1,αk-1∂y=AQk-1,αk-1=4αk-1c2ρ3F3Qk-1c2ρ2F2-Qk-122,∂fyk-1,αk-1∂α=BQk-1,αk-1=2αk-1ρFQk-1c2ρ2F2-Qk-12,fyk-1,αk-1-∂fyk-1,αk-1∂yyk-1-∂fyk-1,αk-1∂ααk-1=CQk-1,αk-1=4αk-1c2ρ3F3Qk-12+4αk-1ρFQk-14c2ρ2F2-Qk-122.Note that by the help of the relations (20) and (21) the matrices Φk-1x,0, Φ1k-1x,0, Φ2k-1x,0 are calculated as follows: (33)Φik-1x,0=∏i=2N-1jE+AQk-1xi1,αk-1h,Φ1ik-1x2N,0=∑j=N+22N-1∏i=2N-1jE+AQk-1xi1,αk-1hBQk-1xj-1,αk-1h+BQk-1x2N-1,αk-1h,Φ2ik-1x2N,0=∑j=N+22N-1∏i=2N-1jE+AQk-1xi1,αk-1hCQk-1xj-1,αk-1h+CQk-1x2N-1,αk-1h,where h is small enough number.Let some statistical data be given.Let us assume that some observation points foryil+0 and yi2l are given (see Table 1).Table 1
y
i
l
+
0
5.5698
5.5732
5.5761
5.5810
5.5848
5.5852
5.5824
y
i
2
l
4.4242
4.4248
4.4254
4.4262
4.4266
4.4263
4.4251We give in Table2 the values of λ obtained by using MATLAB by given input parameters.Table 2
δ
0.05
0.1
0.5
0.54
0.55
0.6
λ
0.1966
0.2141
0.2295
0.2298
0.2299
0.2302As we see from Table2 by δ=0.55λ gets value 0.2299, with error estimation 10−3.Here is the dependence ofλ or δ (see Figure 2).Figure 2The above algorithm reaches given accuracy after 4 iterations and givesλ=0.2299.Note that the inequality(34)λi-λ4≤10-15holds for any i>4 that shows the stability of the proposed quasilinearization algorithm.



---

*Source: 570475-2015-11-01.xml*