# Identifiability for a Class of Discretized Linear Partial Differential Algebraic Equations

Authors: Begoña Cantó; Carmen Coll; Elena Sánchez
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

This paper presents the use of an iteration method to solve the identifiability problem for a class of discretized linear partial differential algebraic equations. This technique consists in replacing the partial derivatives in the PDAE by differences and analyzing the difference algebraic equations obtained. For that, the theory of discrete singular systems, which involves Drazin inverse matrix, is used. This technique can also be applied to other differential equations in mathematical physics.


---

## Body

## 1. Introduction

Mathematical models are used extensively to simulate the behavior of experimental and practical situations. These models can be constructed in different ways and differ greatly in complexity but they have in common that they can be used to predict the behavior of the initial system. Usually, the model is constructed using some of the properties of a system by tinkering with adjustable parameters and sometimes the equations of the models involved unknown parameters. In this case, it is interesting to obtain these parameters to make the model accurate. This property is known as identifiability problem. Given a parameterized state-space model, structural identifiability is concerned with whether the unknown parameters within the model can be identified uniquely from the experiment considered. Thus structural identifiability analysis is an important step in the modeling process, and it is necessary for theoretical prerequisites to experiment design and system identification, that is, to estimate the unknown parameters of the model using experimental data. A wide variety of these models is often defined by differential algebraic equations (DAEs) or partial differential algebraic equations (PDAEs). These equations relate internal variables and their derivatives with inputs to the system and appear in several scientific fields. Recently, several researches have been devoted to PDAEs, and some theoretical and numerical results are obtained [1–4].In this paper, we limit ourselves to solve the identifiability problem for second order partial differential equations which has added additional algebraic equations. The relations can be described by the general explicit formulation (see [1, 5]) using a matrix equation in the form(1.1)E∂z∂t+F(p)∂2z∂x2+G(p)z=f(t,x),
where f(t,x)∈ℝn is a smooth function, E≠0,F(p)≠0,G(p) are matrices in ℝn×n, all of which can be singular, and p=(α,β)∈ℝ2 is a unknown parameter vector. The solution of (1.1), denoted by z=z(t,x), depends on the parameter vector p, the spatial variable x∈Ω and the time value t∈(0,T).For system (1.1) initial values and boundary conditions may depend on the parameter vector to be estimated. In general, if an arbitrary initial condition is imposed, it may be inconsistent with the PDAE; for that, we assume that initial values may be decomposed into the data which can be prescribed arbitrarily and the consistent data. The boundary values are in similar form.The outline of this paper is as follows. Problem description and general formulations are discussed in Section2, also it contains the numerical procedure to discretize the partial differential algebraic equation. Section 3 contains the solution to the identifiability problem of a structural discrete system and the application of the obtained results to solve the initial identifiability problem. For that, an explicit expression of the Drazin inverse of a structured matrix with parameters is given and the structure of the Markov parameters of the singular system is obtained. Finally, in last the section, an algorithm and an illustrative example are shown.


## 2. Problem Description and Discretization

In this section we consider the second-order partial differential algebraic equations (PDAEs) given by(2.1)E∂z∂t+F(p)∂2z∂x2+G(p)z=O,
where z=(z1z2), z1∈ℝn1, z2∈ℝn2 and E=diag(I,O), F(p)=αI, α>0, and G(p)=βI+E, β>4α. The differential algebraic aspect of the system is due to singularity of E.The basic idea to solve the identifiability problem is the use of an iteration method based on the numerical procedure that is useful to discretize partial differential algebraic equations. This technique consists in replacing the partial derivatives in the PDAE by differences. For that, we discretize our domain, both time and space can be discretized by an uniform grid with gridparameterΔx, where Δx is the distance between two neighboured nodes of the grid, and we discretize the interval (0,T) by a one-dimensional grid with stepsize Δt.In this case we use the forward-looking difference operator in time and the approximation is(2.2)∂z∂t≈zk+1,i-zk,iΔt,
and the second-order derivative is approximated by the second-order central difference operator in space given by (2.3)∂2z∂x2≈zk,i-1-2zk,i+zk,i+1(Δx)2.To obtain the discretized system we takeΔx=Δt=1. The discretized problem can by written as (2.4)zk+1,i1=pzk,i1+q(zk,i-11+zk,i+11),0=pzk,i2+q(zk,i-12+zk,i+12),i=1,…,N-1,
where p=2α-β and q=-α.Moreover, we need to add some equations to show the transition among the variables. These new equations are linearly independent of the other equations and have the form(2.5)zk+1,01=pzk,01+qzk,11,zk+1,02=pzk,02+qzk,12+uk,zk+1,N1=pzk,N1+q(zk,N-11+zk,N2),zk+1,N2=pzk,N2+q(zk,N-12+zk,N1).To obtain a structure which simplifies the solution we write the full system of equations as follows:(2.6)Ex(k+1)=A(p)x(k)+Bu(k),
where the state vector x(k)=(z1(k)Tz2(k)T)T with z1(k)=(zk,j1)j=0,…,N and z2(k)=(zk,N-j2)j=0,…,N, the control vector u(k)=uk∈ℝ, and the coefficient matrices are given by(2.7)E=(IlOOO),A(p)=(pq0⋯0qpq⋯00qp⋯0⋮⋱⋯⋱q00⋯qp),B=(0⋮01),
with l=(N+1)n1. These matrices have a fixed structure with the parameter vector p=(p,q) belonging to a subset 𝒫⊆ℝ2, which suggests the relation among variables in the process. It is well-known (see [6]) that the system has solution if there exists an scalar λ∈ℂ such that det[λE-A(p)]≠0 (regularity condition). The matrix A(p) is a symmetric, tridiagonal, and regular matrix. This follows from the eigenvalues of A(p) given by(2.8)λk=p+2qcos⁡(kπn+1),k=1,2,…,n,
and using the expressions of p and q, since it is easy to check that(2.9)λk<-2α(1+cos⁡(kπn+1))<0.
This establishes the regularity of A(p) and this fact guaranties the regularity of the system (2.6).


## 3. The Identification Problem

The system (2.6) is a singular discrete-time system which, from the regularity of A(p), we can rewrite as follows: (3.1)E¯(p)x(k+1)=x(k)+B¯(p)u(k),
where E¯(p)=A-1(p)E and B¯(p)=A-1(p)B. This system is denoted by S(p). Since there exists a scalar λ∈ℂ such that det[λE¯-I]≠0, the explicit solution of the system is given by (3.2)x(k)=E¯D(p)k+1E¯(p)x(0)+∑i=0k-1E¯D(p)k-iB¯(p)u(i)-(I-E¯D(p)E¯(p))∑j=0ν-1E¯(p)jB¯(p)u(k+j),
where E¯D(p) denotes the Drazin inverse of matrix E¯(p), ν is the index of E¯(p), and x(0) is an admissible initial condition, x(0)∈X0, with (3.3)X0=Im⁡[E¯D(p)E¯(p),H0(p),…,Hν-1(p)],
where Hj(p)=(I-E¯D(p)E¯(p))E¯(p)j, j=0,1,…,ν-1.The input-output behavior of the systemS(p) is obtained from the Markov parameters associated which the system. These parameters are defined as(3.4)V(j,p)=E¯D(p)jB¯(p)j≥1.


### 3.1. Identifiability of Structural Singular Systems

To get a good formulation of the model we need to know if the unknown parameters can be determined uniquely from the obtained experiment data. This is known as identifiability problem.Our system (2.6) has a fixed structure, and we want to solve the identifiability problem when this structure holds on. The identifiability helps us test the unique relationship between parameter sets and model response and guarantees that the parameters can be estimated under ideal conditions.A model is globally structurally identifiable if and only if there exists a unique input-output behavior for every parameter set [3]. That is, the system S(p) is globally structurally identifiable if and only if, for almost any two candidates parameter vector values p,q∈𝒫, io(p)=io(q) imply p=q, where io(·) denotes the input-output behavior of the system S(p). From now on, we use the concept globally identifiable as globally structurally identifiable.For linear models there are many well established techniques to analyze structural identifiability; see, for example, [4, 7, 8] and the references therein. The main contribution of this part consists in using the special structure of the Markov parameters associated with S(p) to the identifiability analysis. For that, we need to know the structure of E¯D(p).By definition ofE¯(p) to obtain E¯D(p) we use the inverse of the tridiagonal matrix. A concise expression of A-1(p)=(𝔞ij)i,j is obtained as follows:(3.5)aij=(-1)i+jqαdi-1dn-jdn,
where α=|j-i| and di=pdi-1-q2di-2, for i=2,…,n, with d0=1 and d1=p. Note that di=Mi(p), where Mi(p) is the principal minor of order i of the matrix A(p) and the symmetry property of matrix A(p) is transferred to A-1(p). For more information on the inverse of a tridiagonal matrix see [9, 10].Using the above inverse, an explicit expression of the Drazin inverse ofE¯(p) is given in the following result.Proposition 3.1.
Consider the matrices given in (2.7) and construct E¯(p)=A-1(p)E. Then, the Drazin inverse matrix E¯D(p)=(𝔢ij)i,j is given by
(3.6)eij={pi=1,…,l-1,j=i,q|j-i|=1,0|j-i|>1,qi=l,j=l-1,dn-l+1dn-lj=l,0j≠{l-1,l},(-1)i-l+1qi-l+1dn-idn-li=l+1,…,n,j=l-1,q-1ei,l-1el,lj=l,0j≠{l-1,l}.Proof.
Write the matrixA(p) in a block form
(3.7)A(p)=[A1A2A2TA4].
where A1 and A4 are tridiagonal matrices and A2 is given by
(3.8)A2=[00⋯0⋮⋮⋯⋮00⋯0q0⋯0].
MatrixA-1(p) can be obtained from the Schur complement of matrix A4, S=A1-A2A4-1A2T. Then, we obtain
(3.9)A-1(p)=[S-1*-A4-1A2TS-1*].
From the above expression and by definition of E¯(p)=A-1(p)E,
(3.10)E¯(p)=[S-1O-A4-1A2TS-1O].
Using the expression to obtain the Drazin inverse given in [11] we obtain
(3.11)E¯D(p)=[SO-A4-1A2TS-1S2O]=[SO-A4-1A2TSO].
Using the expression given in (3.5) to obtain A4-1 and by simple mathematical calculations expression (3.6) is obtained.
Now, the next step is to check that this matrixE¯D(p) verifies the properties of the Drazin inverse matrix. First, we need to know the index of E¯(p). From the symmetry of the inverse matrix obtained in (3.5) we obtain rank(E¯(p))=rank(E¯(p)2), then ind(E¯(p))=1.
By technical calculations we check that this matrix verifies the Drazin properties (i)E¯D(p)E¯(p)E¯D(p)=E¯D(p), (ii) E¯D(p)E¯(p)=E¯(p)E¯D(p), and (iii) E¯2(p)E¯D(p)=E¯(p), to prove the proposition.It is important to focus our attention on the structure obtained in (3.6) since it is useful to analyze the input-output behaviour in terms of the vector parameter p. This fact allows to study directly the identification problem.Proposition 3.2.
Consider the systemS(p). Then, the first Markov parameter V(1,p)=(v(1,p)i)i is given by
(3.12)v(1,p)i={0i=1,…,l-1,qal-1,n+dn-l+1dn-lal,ni=l,(-1)i-l+1qi-ldn-idn-lv(1,p)li=l+1,…,n.Proof.
First we prove for alli=1,…,n that
(3.13)q(ai-1,n+ai+1,n)+pai,n=0.
Using the expression given in (3.5) we have that
(3.14)ai-1,n+ai+1,n=(-1)n+i-1qn-i-1dn(q2di-2+di),
and since di=pdi-1-q2di-2 we have
(3.15)q(ai-1,n+ai+1,n)+pai,n=(-1)n+i-1qn-idn(q2di-2+di-pdi-1)=0.
Expression (3.12) is easily seen using the above relation and the Drazin inverse of matrix E¯D(p) given in (3.6).Using Proposition3.2 we give the structure of the rest of Markov parameters.Proposition 3.3.
Consider the systemS(p). Then, the Markov parameters V(j,p)=(v(j,p)i)i, j>1, are given by
(3.16)v(j,p)i={0i=1,…,l-j,qv(j-1,p)i+1+pv(j-1,p)i+qv(j-1,p)i-1i=l-j+1,…,l-1,qv(j-1,p)l-1+dn-l+1dn-lv(j-1,p)li=l,(-1)i-l+1qi-ldn-idn-lv(j,p)li=l+1,…,n.Next, we solve the identifiability problem in the next result.Proposition 3.4.
The systemS(p) is identifiable.Proof.
We consider two structured systemsS(p1) and S(p2) with p1,p2∈𝒫, p1=(p1,q1), and p2=(p2,q2) such that they have the same input-output behavior (io), that is, V(j,p1)=V(j,p2), j≥1.
From expression (3.12), since v(1,p)n=-(q/p)v(1,p)n-1, and taking i=n, we have that
(3.17)V(1,p1)=V(1,p2)implies-q1p1=-q2p2.
On the other hand, from expression (3.16), since v(2,p)l-1=qv(1,p)l,
(3.18)v(2,p1)l-1=v(2,p2)l-1impliesq1=q2.
Hence, p1=p2 and the system is identifiable.


## 3.1. Identifiability of Structural Singular Systems

To get a good formulation of the model we need to know if the unknown parameters can be determined uniquely from the obtained experiment data. This is known as identifiability problem.Our system (2.6) has a fixed structure, and we want to solve the identifiability problem when this structure holds on. The identifiability helps us test the unique relationship between parameter sets and model response and guarantees that the parameters can be estimated under ideal conditions.A model is globally structurally identifiable if and only if there exists a unique input-output behavior for every parameter set [3]. That is, the system S(p) is globally structurally identifiable if and only if, for almost any two candidates parameter vector values p,q∈𝒫, io(p)=io(q) imply p=q, where io(·) denotes the input-output behavior of the system S(p). From now on, we use the concept globally identifiable as globally structurally identifiable.For linear models there are many well established techniques to analyze structural identifiability; see, for example, [4, 7, 8] and the references therein. The main contribution of this part consists in using the special structure of the Markov parameters associated with S(p) to the identifiability analysis. For that, we need to know the structure of E¯D(p).By definition ofE¯(p) to obtain E¯D(p) we use the inverse of the tridiagonal matrix. A concise expression of A-1(p)=(𝔞ij)i,j is obtained as follows:(3.5)aij=(-1)i+jqαdi-1dn-jdn,
where α=|j-i| and di=pdi-1-q2di-2, for i=2,…,n, with d0=1 and d1=p. Note that di=Mi(p), where Mi(p) is the principal minor of order i of the matrix A(p) and the symmetry property of matrix A(p) is transferred to A-1(p). For more information on the inverse of a tridiagonal matrix see [9, 10].Using the above inverse, an explicit expression of the Drazin inverse ofE¯(p) is given in the following result.Proposition 3.1.
Consider the matrices given in (2.7) and construct E¯(p)=A-1(p)E. Then, the Drazin inverse matrix E¯D(p)=(𝔢ij)i,j is given by
(3.6)eij={pi=1,…,l-1,j=i,q|j-i|=1,0|j-i|>1,qi=l,j=l-1,dn-l+1dn-lj=l,0j≠{l-1,l},(-1)i-l+1qi-l+1dn-idn-li=l+1,…,n,j=l-1,q-1ei,l-1el,lj=l,0j≠{l-1,l}.Proof.
Write the matrixA(p) in a block form
(3.7)A(p)=[A1A2A2TA4].
where A1 and A4 are tridiagonal matrices and A2 is given by
(3.8)A2=[00⋯0⋮⋮⋯⋮00⋯0q0⋯0].
MatrixA-1(p) can be obtained from the Schur complement of matrix A4, S=A1-A2A4-1A2T. Then, we obtain
(3.9)A-1(p)=[S-1*-A4-1A2TS-1*].
From the above expression and by definition of E¯(p)=A-1(p)E,
(3.10)E¯(p)=[S-1O-A4-1A2TS-1O].
Using the expression to obtain the Drazin inverse given in [11] we obtain
(3.11)E¯D(p)=[SO-A4-1A2TS-1S2O]=[SO-A4-1A2TSO].
Using the expression given in (3.5) to obtain A4-1 and by simple mathematical calculations expression (3.6) is obtained.
Now, the next step is to check that this matrixE¯D(p) verifies the properties of the Drazin inverse matrix. First, we need to know the index of E¯(p). From the symmetry of the inverse matrix obtained in (3.5) we obtain rank(E¯(p))=rank(E¯(p)2), then ind(E¯(p))=1.
By technical calculations we check that this matrix verifies the Drazin properties (i)E¯D(p)E¯(p)E¯D(p)=E¯D(p), (ii) E¯D(p)E¯(p)=E¯(p)E¯D(p), and (iii) E¯2(p)E¯D(p)=E¯(p), to prove the proposition.It is important to focus our attention on the structure obtained in (3.6) since it is useful to analyze the input-output behaviour in terms of the vector parameter p. This fact allows to study directly the identification problem.Proposition 3.2.
Consider the systemS(p). Then, the first Markov parameter V(1,p)=(v(1,p)i)i is given by
(3.12)v(1,p)i={0i=1,…,l-1,qal-1,n+dn-l+1dn-lal,ni=l,(-1)i-l+1qi-ldn-idn-lv(1,p)li=l+1,…,n.Proof.
First we prove for alli=1,…,n that
(3.13)q(ai-1,n+ai+1,n)+pai,n=0.
Using the expression given in (3.5) we have that
(3.14)ai-1,n+ai+1,n=(-1)n+i-1qn-i-1dn(q2di-2+di),
and since di=pdi-1-q2di-2 we have
(3.15)q(ai-1,n+ai+1,n)+pai,n=(-1)n+i-1qn-idn(q2di-2+di-pdi-1)=0.
Expression (3.12) is easily seen using the above relation and the Drazin inverse of matrix E¯D(p) given in (3.6).Using Proposition3.2 we give the structure of the rest of Markov parameters.Proposition 3.3.
Consider the systemS(p). Then, the Markov parameters V(j,p)=(v(j,p)i)i, j>1, are given by
(3.16)v(j,p)i={0i=1,…,l-j,qv(j-1,p)i+1+pv(j-1,p)i+qv(j-1,p)i-1i=l-j+1,…,l-1,qv(j-1,p)l-1+dn-l+1dn-lv(j-1,p)li=l,(-1)i-l+1qi-ldn-idn-lv(j,p)li=l+1,…,n.Next, we solve the identifiability problem in the next result.Proposition 3.4.
The systemS(p) is identifiable.Proof.
We consider two structured systemsS(p1) and S(p2) with p1,p2∈𝒫, p1=(p1,q1), and p2=(p2,q2) such that they have the same input-output behavior (io), that is, V(j,p1)=V(j,p2), j≥1.
From expression (3.12), since v(1,p)n=-(q/p)v(1,p)n-1, and taking i=n, we have that
(3.17)V(1,p1)=V(1,p2)implies-q1p1=-q2p2.
On the other hand, from expression (3.16), since v(2,p)l-1=qv(1,p)l,
(3.18)v(2,p1)l-1=v(2,p2)l-1impliesq1=q2.
Hence, p1=p2 and the system is identifiable.


## 4. Algorithm and Example

In this section we introduce an algorithm to obtain the parameters of a PDAE given by (2.1), using the Markov parameters technique.Step 1.
Introduce the size of the state vector,n. Introduce the matrices {V(j)=(v(j)i)i=1,…,n,j=1,…,n} that determine the known input-output behavior of the process.Step 2.
Fori=1,…,n, if v(1)i≠0, then rename l=i.Step 3.
Constructq=v(2)l-1/v(1)l.Step 4.
Constructp=-qv(1)n-1/v(1)n.Step 5.
Introduced0=1 and d1=p.Step 6.
Fori=2,…,n, define di=pdi-1-q2di-2.Step 7.
Fori=l-1,l, define 𝔞i,n=(-1)i+nq|n-i|(di-1/dn).Fori=1,…,n, apply the following.Step 8.
Prove thatV(1) verifies
(4.1)v(1)i={0i=1,…,l-1,qal-1,n+dn-l+1dn-lal,ni=l,(-1)i-l+1qi-ldn-idn-lv(1)li=l+1,…,n.
If the result is false, then stop the process. This means that the given input-output behavior is not suitable to a system of type (2.1).Step 9.
Prove thatV(j), j=2,…,n, verifies
(4.2)v(j)i={0i=1,…,l-j,qv(j-1)i+1+pv(j-1)i+qv(j-1)i-1i=l-j+1,…,l-1,qv(j-1)l-1+dn-l+1dn-lv(j-1)li=l,(-1)i-l+1qi-ldn-idn-lv(j-1)li=l+1,…,n.
If the result is false, then stop the process. This means that the given input-output behavior is not suitable to a system of type (2.1).Step 10.
Obtainα=-q and β=2α-p.Step 11.
Obtain the matrices of the system (2.6).Step 12.
Obtain the matrices of the system (2.1).This algorithm has been implemented in MATLAB. To illustrate this algorithm, we present the following academic example.Example 4.
Consider the PDAE given in (2.1), with n1=n2=1. The discrete system associated with this PDAE is given by (2.6) with A(p)∈ℝ6. The input-output behavior of a process is given by
(4.3)V(1)=(000-43-89-169),V(2)=(0049-1427-89-169),V(3)=(0-4272081-85243-2881-5681),V(4)=(481-26243151729-9554374-170729-340729),V(5)=(32729-22621875633938-11127223-4393016-4391508),V(6)=(2746561-51364526976360-108310072-107510474-10755237).
Then we apply the above algorithm.First, we introduce the size of the state vector,n=6. Introduce the matrices {V(j)=(v(j)i)i=1,…,n,j=1,…,n} that determine the known input-output behavior of the process, that is,(4.4)v(:,1)=[0;0;0;-43;-89;-169],v(:,2)=[0;0;49;-1427;-89;-169],v(:,3)=[0;-427;2081;-85243;-2881;-5681],v(:,4)=[481;-26243;151729;-9554374;-170729;-340729],v(:,5)=[32729;-2262187;5633938;-11127223;-4393016;-4391508],v(:,6)=[2746561;-5136452;6976360;-108310072;-107510474;-10755237],V1=v(:,1),V2=v(:,2),V3=v(:,3),V4=v(:,4),V5=v(:,5),V6=v(:,6).Applying MATLAB algorithm we obtainl=4, q=-1/3, and p=1/6. The comparison is ok; a=1/3; b=1/2.The discrete model (2.6) is given by(4.5)E=diag⁡(I4,O2),A=(16-130000-1316-130000-1316-130000-1316-130000-1316-130000-1316),B=(000001).


## 5. Conclusions

In this paper we show some results to identify and determine parameters in a mathematical model with a particular structure. To obtain these results we study a discrete model associated with a partial differential algebraic equation. This model has been obtained using finite difference method from the partial differential equation. Some results about the Markov parameters of the discrete system are given, and the identifiability problem is solved in the discrete case. By a similar way, the given results can be applied to different models with a similar structure to different areas.In this case, the results are used to know if a partial differential algebraic equation with unknown parameters is identifiable and to obtain the parameter values. Finally, an algorithm, implemented in MATLAB, is given to obtain the parameters.



---

*Source: 510519-2011-04-26.xml*