# Online Identification of Multivariable Discrete Time Delay Systems Using a Recursive Least Square Algorithm

Authors: Saïda Bedoui; Majda Ltaïef; Kamel Abderrahim
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper addresses the problem of simultaneous identification of linear discrete time delay multivariable systems. This problem involves both the estimation of the time delays and the dynamic parameters matrices. In fact, we suggest a new formulation of this problem allowing defining the time delay and the dynamic parameters in the same estimated vector and building the corresponding observation vector. Then, we use this formulation to propose a new method to identify the time delays and the parameters of these systems using the least square approach. Convergence conditions and statistics properties of the proposed method are also developed.
Simulation results are presented to illustrate the performance of the proposed method. An application of the developed approach to compact disc player arm is also suggested in order to validate simulation results.


---

## Body

## 1. Introduction

Time delay system identification has received great attention in the last years since time delay is a physical phenomenon which arises in most control loops industrial systems [1, 2]. Several reasons cause the presence of time delay in control loops. In fact, it may be an inherent feature of the system such as processes of transport (mass, energy, and information), higher order processes, and accumulation of time lags in several systems that are connected in series. It may also be introduced by the devices of control loops, such as response times of sensors and actuators, computation time of control laws, and information transmission time in networks. This delay can be neglected if its value is too small for the system time constants. Otherwise, it cannot be neglected, and the dynamic representation of the system must be described by a time delay model. This model is, generally, constructed using the identification approach which allows building a mathematical model from input-output data.The identification of time delay systems is known to be a challenging identification problem because it involves both the estimation of dynamic parameters and time delay. Numerous methods have been proposed in the literature for the identification of time delay systems [3–15].Among these methods, the graphical approach has been the most popular since it represents the first method proposed in the literature for the identification of continuous time delay systems [16]. Moreover, it is frequently used to compute the parameters of PID controllers for industrial processes. It consists in determining the time delay and the dynamic parameters of the system from its step response. The main advantage of this approach lies in the simplicity of its implementation. However, it produces inaccurate results because it is very sensitive to measurement noises. Another popular approach is proposed in [9]. It is based on the approximation of the time delay by a rational transfer function using classical approximations such as Pade or Laguerre. This method can insure very satisfactory results in the case of linear systems with constant time delays and lower order. However, its performance degrades rapidly in the case of higher order systems or important time delays. Moreover, it raises the computational complexity because it increases the number of parameters to be estimated. The parametrisation approach can be considered as one of the interesting methods because it is based on theoretical concepts of discrete time systems [17]. It consists, firstly, in inserting a known time delay in the numerator of the discrete time model, secondly, in estimating the dynamic parameters of the system using a recursive algorithm, and, finally, in deducing the time delay from zero coefficients of the numerator. In practice, it is difficult or rather impossible to have zero coefficients from experimental data. Indeed, we must set a threshold, which is a delicate task, mainly in the case of a noisy output. The method developed in [3] consists, It consists, firstly, in using the recursive least square approach to identify the parameters assuming that the time delay is known, and secondly, in estimating the time delay, taking into account the results of the first step. The time delay may be identified either by maximizing the correlation function or by minimizing the quadratic error. This method assumes that the domain range of the time delay is a priori known. We also mention the method presented in [18]. It allows the identification of the time delay and the system parameters using the Levenberg-Marquaydt optimization approach to minimize the prediction error. An online identification algorithm for continuous-time single-input single-output (SISO) linear time delay systems with uncertain time invariant parameters is developed in [10]. It consists in constructing a sliding mode-based observer of an underlying system with uncertain parameters. This observer is then used to design an adaptive identifier of system parameters. A linear filtering method is introduced in [19] for simultaneous parameter and time delay estimation of transfer function models. This method estimates the time delay along other model using an iterative way through simple linear regression. Another method that identifies the time delay and the system parameters is presented in [20] which is based on the correlation technique. The method developed in [21] allows the identification of time delay and the parameters of a system operating in the presence of colored noise. This method is based on correlation analysis. The method developed in [12, 22] allows the identification of time delay and the parameters. It minimizes the error between the process output and the process predictive model output, and then the variable time delay parameter is identified. In our previous work, we have proposed two methods for the simultaneous identification of the time delay and dynamic parameters of monovariable time delay discrete systems. The first method is based on the least square approach [23, 24]. The second method consists in minimizing a quadratic criterion using the gradient approach [25].Most of these approaches deal with the problem of the identification of single-input single-output (SISO) time delay systems. However, the problem of multi-input multioutput (MIMO) time delay systems is one of the most difficult problems that represents an area of research where few efforts have been devoted in the past. The use of time delay approximation is extended to the MIMO case [26]. In fact, the authors deal with the problem of the identification of time delay processes using an overparameterization method. In [27], a method is developed for time delay estimation in the frequency domain of MIMO systems based on the combination of continuous wavelet transform (CWT) and cross-correlation. During the estimation, cross-correlation computations are carried out between the CWT coefficients of the input and the output data. The authors of [28] have proposed a simple method based on the combination of two well-known approaches: time delay estimation from impulse response and subspace identification.In this paper, we propose an alternative approach for the problem of simultaneous identification of linear discrete time delay multivariable systems. Indeed, we develop a new formulation of the problem allowing to define the time delay and the dynamic parameters in the same estimated vector and to build the corresponding observation vector. Then, we use this formulation to propose a new method to identify the time delays and the parameters of these systems using the least square approach. Convergence conditions and statistics properties of the proposed method are also developed. Simulation and experimental examples are presented to illustrate the effectiveness of the proposed methods and to compare their performance in terms of convergence and speed. Our approach presents several interesting properties which can be summarized as follows.(i)
The simultaneous identification of the time delays and parameters matrices is achieved by a new formulation of the parameters matrices.(ii)
No a priori knowledge of the time delay is required. In fact, most of the publications assume the knowledge of the time delay variation range or the initial condition.(iii)
The consistency of recursive least square methods has received much attention in the identification literature. In this paper, the proof of the consistency of the estimates is established.(iv)
It can be used to deal with control adaptive purposes.This paper is organized as follows. Section2 presents the model and its assumptions. In Section 3, we propose an extended least square algorithm for simultaneous online identification of unknown time delays and parameters of multivariable discrete time delay systems. Moreover, we develop the convergence properties of the estimates in order to show that the obtained estimates are unbiased. Simulation results and experimental test are provided in the last section.


## 2. Problem Statement

In this paper, we address the problem of identification of square linear multivariable delay system withp inputs and p  ARX model:
(1)A(q-1)Y(k)=B(q-1)U~(k)+V(k),
where Y(k)=[y1(k),…,yp(k)]T and U~(k)=[U1(k-d1),…,Up(k-dp)]T are the outputs and the delayed inputs of the system at time k and V(k)=[V1(k),…,Vp(k)]T is a vector of independent random variables sequences. Let {D=diag[d1,…,dp]/di∈ℕ,i=1,…,p} be the time delay diagonal matrix, also called the interactive matrix, and A(q-1) and B(q-1) two polynomial matrices in the unit backward shift operator q-1(i.e.,q-1yi(k)=yi(k-1),i=1,…,p), defined by
(2)A(q-1)=Ip+A1q-1+⋯+Anaq-na,dimAr=(p,p),r∈[1,na],B(q-1)=B1q-1+⋯+Bnbq-nb,dimBr=(p,p),r∈[1,nb].
The delayed inputs U~(k) can be expressed as
(3)U~(k)=ΩU(k),
where Ω is a diagonal matrix defined as
(4)Ω=(q-d10⋯00q-d2⋱⋮⋮⋱⋱00⋯0q-dp).
The following assumptions are made.(A1)
The two polynomial matricesA(q-1) and B(q-1) have no common left factor.(A2)
The ordersna and nb of the model are known.(A3)
The input sequences{U~=[U1(k-d1),…,Up(k-dp)]T} are independent of V(k), mutually independent and identically distributed with E[U(k)]=0 and E[U(k)U(k)T]=I, and are persistently exciting (PE).(A4)
The disturbanceV(k)=[V1(k),…,Vp(k)]T is sequences of independent, identically distributed random variables with zero mean and finite variance Σ={σ12,…,σp2}.(A5)
The inputs, the outputs, and the noises are causal; that is,U(k)=[0], Y(k)=[0], and V(k)=[0] for k⩽0.Problem Statement. The goal is to develop a recursive algorithm to estimate, simultaneously, the time delay matrix D and the matrices {Ai(k),Bi(k)} using the input/output measurement data {U(k),Y(k)}.In the following, we present three necessary definitions.Definition 1.
Operatorround(d) is defined by
(5)round(di)={int(di)+1ifdi-int(di)⩾0.5int(di)ifdi-int(di)<0.5,
where int(d) denotes the integer part of di,i=1,…,p.Definition 2.
Operatord~(·) is defined by
(6)d~i(·)=round(d^i(·)),i=1,…,p.Definition 3.
OperatorD~(·) is defined by
(7)D~(·)=(d~10⋯00d~2⋱⋮⋮⋱⋱00⋯0d~p).


## 3. The Proposed Approach

In this section, an extended least square algorithm for simultaneous online identification of time delays and parameter matrices is developed.Equation (1) can be rewritten as
(8)Y(k)=ΘTφ(k,D)+V(k),
where Θ is the parameter matrix and φ(k,D) is the observation vector defined as
(9)ΘT=[A1,A2,…,Ana,B1,B2,…,Bnb],φ(k,D)=[-YT(k-1),-YT(k-2),…,-YT(k-na),=..U~T(k-1),…,U~T(k-nb)]T.On the other hand, the estimated output is described by the following relation:(10)Y^(k)=Θ^Tφ(k,D~),
where Θ^ and D~=diag[d~1,…,d~p] represent, respectively, the estimated parameter matrix and the estimated time delay matrix.Let us consider the prediction errorΥ(k) given by
(11)Υ(k)=Y(k)-Θ^Tφ(k,D~).Since parameter matrixΘ^ does not contain the unknown time delays D~, then consequently it is not directly applicable to achieve our objective which is the simultaneous identification of the time delays and the parameter matrices of the multivariable discrete time delay systems (1).To overcome this problem, we suggest considering the time delay matrix in parameter matricesΘ to be estimated. Indeed, the new matrix, called generalized matrix, is given by
(12)ΘGT=[Θ,D].Moreover, we propose the use of the negative gradient of the error to obtain an appropriate observation vector which is given by(13)Φ(k,Θ^G)=-∂Υ∂Θ^G.
Then,
(14)Φ(k,Θ^G)=[φT(k,D~),-∂Υ∂D~]T,Φ(k,Θ^G)=[φT(k,D~),∂Y^∂D~]T,Φ(k,Θ^G)=[φT(k,D~),∂∂D^∑i=1nb‍UT(k-i)B^iΩ^]T.
The use of the approximation of Ln(q)≈1-q-1 (see the appendix) leads to
(15)Φ(k,Θ^G)=[φT(k,D~),-∑i=1nb‍ΔUT(k-i)B^iΩ^]T,
where ΔU(k)=U(k)-U(k-1).ReplacingφT(k,D~) by its expression, we obtain the generalized observation vector:
(16)Φ(k,Θ^G)=[-YT(k-1)⋮-YT(k-na)Ω^UT(k-1)⋮Ω^UT(k-nb)-∑i=1nb‍ΔUT(k-i)B^iΩ^]An estimationΘ^G of ΘG is denoted by the minimization of the following criterion:
(17)J(k,ΘG)=12∑i=0k‍Υ(i)2.
Then, the partial derivative of the criterion with respect to the generalized matrix is
(18)∂J∂Θ^G=∑i=0k‍∂Υ(i)∂Θ^GΥ(i)=-∑i=0k‍Φ(i,Θ^G)Υ(i).
So
(19)∂J∂Θ^G=-∑i=0k‍Φ(i,Θ^G)[Y(i)-Θ^Tφ(i,D~)].
Let us consider
(20)ψ=-∑j=1nb‍D~ΔUT(i-j)B^jΩ^.
Adding and subtracting from (19) the term ψ, given by (20), we have
(21)∂J∂Θ^G=-∑i=0k‍Φ(i,Θ^G)[Y(i)-ΦT(i,Θ^G)Θ^G+ψ]=-∑i=0k‍Φ(i,Θ^G)[Y(i)+ψ-ΦT(i,Θ^G)Θ^G].
Canceling the partial derivative of the criterion, we obtain
(22)Θ^G(k)=[∑i=0k‍Φ(i,Θ^G)ΦT(i,Θ^G)]-1∑i=0k‍Φ(i,Θ^G)×[Y(i)+ψ].
Let,
(23)R(k)=∑i=0k‍Φ(i,Θ^G)ΦT(i,Θ^G).
Based on assumption A3 which ensures that matrix R(k) is invertible [29], (22) can be rewritten as
(24)Θ^G(k)=R(k)-1∑i=0k‍Φ(i,Θ^G)[Y(i)+ψ].
It follows that
(25)Θ^G(k)=R(k)-1(∑i=0k-1‍Φ(i,Θ^G)[Y(i)+ψ]=R(k)-1bb+Φ(k,Θ^G)[Y(k)+ψ]∑i=0k-1‍).Using (22), we have
(26)Θ^G(k)=R(k)-1(R(k-1)Θ^G(k-1)=R(k.)-1+Φ(k,Θ^G)[Y(k)+ψ])=R(k)-1(R(k)Θ^G(k-1)=R(k)-1b-Φ(k,Θ^G)ΦT(k,Θ^G)Θ^G(k-1)=R(k)-1b+Φ(k,Θ^G)[Y(k)+ψ]).
So
(27)Θ^G(k)=Θ^G(k-1)+R(k)-1Φ(k,Θ^G)×(-ΦT(k,Θ^G)Θ^G(k-1)+[Y(k)+ψ]).
It follows from (27) that
(28)Θ^G(k)=Θ^G(k-1)+R(k)-1Φ(k,Θ^G)×(Y(k)-Θ^Tφ(k,D~)).
Thus,
(29)Θ^G(k)=Θ^G(k-1)+P(k)Φ(k,Θ^G)Υ(k),(30)P(k)=[∑i=0k‍Φ(i,Θ^G)ΦT(i,Θ^G)]-1=[R(k-1)+Φ(k,Θ^G)ΦT(k,Θ^G)]-1.Using the matrix inversion lemma given by [29],
(31)[B+CDT]-1=B-1-B-1CDTB-1[1+DTB-1C]-1.
Let B=R(k-1), C=Φ(k,ΘG), and D=Φ(k,ΘG)T, and then we have
(32)P(k)=P(k-1)-P(k-1)Φ(k,Θ^G)ΦT(k,Θ^G)P(k-1)1+ΦT(k,Θ^G)P(k-1)Φ(k,Θ^G).
The previous approach can be summarized by  Algorithm 1.Algorithm 1
Step
1. Data acquisition {U(k),Y(k)} and initialization:
setΘ^G=ΘG0 and P=βI where β is a scalar and
I is the identity matrix of size
(p,p(na+nb+nc+1)) and k=0,
Step
2. Increment k and construct the observation vector
φ(k,D~), the generalized observation vector
ΦT(k,Θ^G) using (9) and (16),
Step
3. Estimate Θ^G using the developed identification
method:
Υ(k)=Y(k)-Θ^T(k-1)φ(k,D~)
P(k)=P(k-1)-P(k-1)Φ(k,Θ^G)ΦT(k,Θ^G)P(k-1)1+ΦT(k,Θ^G)P(k-1)Φ(k,Θ^G)
Θ^G(k)=Θ^G(k-1)+P(k)Φ(k,Θ^G)Υ(k)
Step
4. Return to Step 2 until k=N where N is the
number of input/output data.


## 4. Convergence Properties

### 4.1. Consistency

For the system in (1), assume that (A3) and (A4) hold. Then for any c>1, the parameter estimation error, Θ^G→ΘG, associated with the LS algorithm in (29) and (30) satisfies
(33)∥Θ^G-ΘG∥2=O([Lnr(k)]cλmin[P-1(k)]),
where r(k) is the trace of the covariance matrix P-1(k) and λmin represents the minimum eigenvalues of P-1(k).Proof.
Define the parameter estimation error vector:(34)Θ-G:=Θ^G-ΘG.
Using (27) and (11) we have
(35)Θ-G(k)=Θ-G(k-1)+PΦ(k,ΘG)×[Θ^T(k-1)φ(k,D^)+P.c+V(k)-ΘT(k-1)φ(k,D)]=Θ-G(k-1)+PΦ(k,ΘG)[-ρ-(k)+V(k)],
where
(36)ρ-(k)=Θ^T(k-1)φ(k,D^)-ΘT(k-1)φ(k,D).
Equation (35) can be rewritten as (see the appendix)
(37)Θ-G(k)=Θ-G(k-1)+PΦ(k,ΘG)[-Y-(k)-ξ(k)+V(k)],
where
(38)ξ(k)=∑i=1nb‍Bi(Ω^-Ω)UT(k-i)+∑i=1nb‍(D^B^i-DBi)Ω^ΔUT(k-i),Y-(k)=Θ¯GT(k-1)Φ(k,Θ^G).
Let us define now a nonnegative definite function:(39)S(k):=Θ¯GT(k)P-1(k)Θ-G(k).
ReplacingΘ¯GT(k) by its expression, we obtain
(40)S(k)=[Θ-G(k-1)+P(k)Φ(k,Θ^G)bbb×(-Y-(k)-ξ(k)+V(k))]TP-1(k)×[Θ-G(k-1)+P(k)Φ(k,Θ^G)bbbbb×(-Y-(k)-ξ(k)+V(k))].
So
(41)S(k)=Θ¯GT(k-1)P-1(k)Θ-G(k-1)+2Θ¯GT(k-1)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]+ΦT(k,Θ^G)P(k)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]×[-Y-(k)-ξ(k)+V(k)]T.
SubstitutingP-1(k) by (30), we get
(42)S(k)=Θ¯GT(k-1)[P-1(k-1)bbbbbbbb.bbb.+Φ(k,Θ^G)ΦT(k,Θ^G)]×Θ-G(k-1)+2Θ¯GT(k-1)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]+ΦT(k,Θ^G)P(k)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]×[-Y-(k)-ξ(k)+V(k)]T.
Then,(43)S(k)=S(k-1)+Y-(k)Y-(k)T-2Y-(k)Y-(k)T-2Y-(k)ξT(k)+2Y-(k)VT(k)+ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)Y-(k)Y-(k)T-2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)Y-(k)×(-ξ(k)+V(k))T+ΦT(k,Θ^G)×P(k)Φ(k,Θ^G)(ξ(k)ξ(k)T+V(k)V(k)Tbbbbbbbbb..bbbbbb-2ξ(k)VT(k)).
It follows from (43) that
(44)S(k)=S(k-1)-[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)Y-(k)T+ΦT(k,Θ^G)×P(k)Φ(k,Θ^G)V(k)V(k)T+2[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)VT(k)-2×[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)ξT(k)+ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)ξ(k)ξ(k)T-2ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)ξ(k)VT(k).
We have
(45)-2[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)ξT(k)=[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)Y-(k)T+[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×ξ(k)ξ(k)T-[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×(Y-(k)+ξ(k))(Y-(k)+ξ(k))T,
and the use of the relation
(46)1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)=[1+ΦT(k,Θ^G)P(k-1)Φ(k,Θ^G)]-1⩾0
leads to
(47)S(k)⩽S(k-1)+ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)V(k)V(k)T+2×[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)VT(k)+ξ(k)ξ(k)T-2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)ξ(k)VT(k).
SinceY-(k), ΦT(k,Θ^G)P(k)Φ(k,Θ^G), ξ(k), and S(k-1) are uncorrelated with V(k), taking the conditional expectation on both sides of (47) and using (A3) and (A4), we obtain
(48)E[S(k)]⩽S(k-1)+2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)Σ+β.
Define
(49)W(k):=S(k)[Ln|P-1(k)|]c.
Since[Ln|P-1(k)|]c is nondecreasing, we have
(50)E[W(k)]⩽S(k-1)[Ln|P-1(k)|]c+2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)[Ln|P-1(k)|]cΣ+β[Ln|P-1(k)|]c⩽W(k-1)+2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)[Ln|P-1(k)|]cΣ+β[Ln|P-1(k)|]c.
Using this property∑i=1∞‍(ΦT(i,Θ^G)P(i)Φ(i,Θ^G)/[Ln|P-1(i)|]c)<∞ (the proof in the same way as the proof of Lemma 1 in [30]), we can see that the sum of the right-hand second term of equation (50) for k from 1 to ∞ is finite. Applying the martingale convergence theorem [30] to the previous inequality, we conclude that W(k) converges a.s. to a finite random variable, say W0; that is,
(51)W(k)=S(k)[Ln|P-1(k)|]c⟶W0<∞
is equivalent to
(52)S(k)=O([Ln|P-1(k)|]c).
From the definition of V(k), we obtain
(53)∥Θ-G∥⩽Θ-G(k)TP-1(k)Θ-G(k)λmin[P-1(k)]=S(k)λmin[P-1(k)].
Now, let us define the matrix trace(54)r(k)=tr[P-1(k)].
It follows that:
(55)|P-1(k)|⩽rn(k),r(k)⩽nλmax|P-1(k)|,Ln|P-1(k)|=O(Lnr(k))=O(Lnλmax|P-1(k)|).
We obtain, finally,(56)∥Θ^G-ΘG∥2=O([Ln|P-1(k)|]cλmin|P-1(k)|)=O([Lnr(k)]cλmin|P-1(k)|)=O([λmax|P-1(k)|]cλmin|P-1(k)|).


### 4.2. Lemma

For the estimate (22) with the assumption (A4), the following proprieties hold.
(P1)
Θ
^
G
is an unbiased estimate of
Θ
G
.
(P2)
The covariance matrix of
Θ
^
G
is given by
(57)
E
[
(
Θ
^
G
-
Θ
G
)
(
Θ
^
G
-
Θ
G
)
T
]
=
(
∑
i
=
0
k
‍
Φ
(
i
,
Θ
G
)
Φ
(
i
,
Θ
G
)
T
)
-
1
Σ
,
where
(58)
Σ
=
(
σ
1
2
0
⋯
0
0
σ
2
2
⋱
⋮
⋮
⋱
⋱
0
0
⋯
0
σ
p
2
)
.Proof.
If we replace (8) in (22), we have
(59)Θ^G=[∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)[ΘTφ(i,D)+V(i)+ψ].
Then,
(60)Θ^G=[∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)[ΘTφ(i,D)+ψ]+[∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)V(i).
So
(61)E[Θ^G]=ΘG+E([∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)V(i)).
Since V(i) is uncorrelated with the elements of Φ(i,ΘG) (13), then
(62)E[Θ^G]=ΘG
which proves (P1).
Consider the first-order Taylor series expansion around the real matrix ofΘG:
(63)∂J(k,Θ^G)∂Θ^G=∂J(k,ΘG)∂ΘG+∂2J(k,ΘG)∂ΘG2(Θ^G-ΘG).
Since∂J(k,Θ^G)/∂Θ^G=0, it derives from (63) that
(64)(Θ^G-ΘG)(Θ^G-ΘG)T=[∂2J(k,ΘG)∂ΘG2]-1∂J(k,ΘG)∂ΘG×[∂J(k,ΘG)∂ΘG]T[(∂2J(k,ΘG)∂ΘG2)-1]T.
The second partial derivative of the criterion with respect to the generalized matrix is(65)∂2J(k,ΘG)∂ΘG2=∑i=0k‍(Υ(i)∂2Υ(i)∂2ΘG-∂Υ(i)∂ΘGΦ(i,ΘG)).
So
(66)∂2J(k,ΘG)∂ΘG2=∑i=0k‍(Υ(i)∂2Υ(i)∂2ΘG+Φ(i,ΘG)ΦT(i,ΘG)).
The use of the small residual algorithms [31] leads to neglect the following term, then:
(67)∑i=0k‍Υ(i)∂2Υ(i)∂2ΘG⟶0.
Hence, an approach of ∂2J(k,ΘG)/∂ΘG2 is obtained:
(68)∂2J(k,ΘG)∂ΘG2≃∑i=0k‍(Φ(i,ΘG)ΦT(i,ΘG)).
Applying the mean value of (∂J(k,ΘG)/∂ΘG)(∂J(k,ΘG)/∂ΘGT), we get:
(69)E[∂J(k,ΘG)∂ΘG∂J(k,ΘG)∂ΘGT]=∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)TE(Υ(i)Υ(i)T).
So
(70)E[∂J(k,ΘG)∂ΘG∂J(k,ΘG)∂ΘGT]=∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)TΣ.
Then, we have
(71)E[(Θ^G-ΘG)(Θ^G-ΘG)T]=E[(∂2J∂ΘG2)-1∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)TΣ[(∂2J∂ΘG2)-1]T].
Finally, we obtain
(72)E[(Θ^G-ΘG)(Θ^G-ΘG)T]=(∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)T)-1Σ
which proves (P2).


## 4.1. Consistency

For the system in (1), assume that (A3) and (A4) hold. Then for any c>1, the parameter estimation error, Θ^G→ΘG, associated with the LS algorithm in (29) and (30) satisfies
(33)∥Θ^G-ΘG∥2=O([Lnr(k)]cλmin[P-1(k)]),
where r(k) is the trace of the covariance matrix P-1(k) and λmin represents the minimum eigenvalues of P-1(k).Proof.
Define the parameter estimation error vector:(34)Θ-G:=Θ^G-ΘG.
Using (27) and (11) we have
(35)Θ-G(k)=Θ-G(k-1)+PΦ(k,ΘG)×[Θ^T(k-1)φ(k,D^)+P.c+V(k)-ΘT(k-1)φ(k,D)]=Θ-G(k-1)+PΦ(k,ΘG)[-ρ-(k)+V(k)],
where
(36)ρ-(k)=Θ^T(k-1)φ(k,D^)-ΘT(k-1)φ(k,D).
Equation (35) can be rewritten as (see the appendix)
(37)Θ-G(k)=Θ-G(k-1)+PΦ(k,ΘG)[-Y-(k)-ξ(k)+V(k)],
where
(38)ξ(k)=∑i=1nb‍Bi(Ω^-Ω)UT(k-i)+∑i=1nb‍(D^B^i-DBi)Ω^ΔUT(k-i),Y-(k)=Θ¯GT(k-1)Φ(k,Θ^G).
Let us define now a nonnegative definite function:(39)S(k):=Θ¯GT(k)P-1(k)Θ-G(k).
ReplacingΘ¯GT(k) by its expression, we obtain
(40)S(k)=[Θ-G(k-1)+P(k)Φ(k,Θ^G)bbb×(-Y-(k)-ξ(k)+V(k))]TP-1(k)×[Θ-G(k-1)+P(k)Φ(k,Θ^G)bbbbb×(-Y-(k)-ξ(k)+V(k))].
So
(41)S(k)=Θ¯GT(k-1)P-1(k)Θ-G(k-1)+2Θ¯GT(k-1)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]+ΦT(k,Θ^G)P(k)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]×[-Y-(k)-ξ(k)+V(k)]T.
SubstitutingP-1(k) by (30), we get
(42)S(k)=Θ¯GT(k-1)[P-1(k-1)bbbbbbbb.bbb.+Φ(k,Θ^G)ΦT(k,Θ^G)]×Θ-G(k-1)+2Θ¯GT(k-1)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]+ΦT(k,Θ^G)P(k)Φ(k,Θ^G)×[-Y-(k)-ξ(k)+V(k)]×[-Y-(k)-ξ(k)+V(k)]T.
Then,(43)S(k)=S(k-1)+Y-(k)Y-(k)T-2Y-(k)Y-(k)T-2Y-(k)ξT(k)+2Y-(k)VT(k)+ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)Y-(k)Y-(k)T-2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)Y-(k)×(-ξ(k)+V(k))T+ΦT(k,Θ^G)×P(k)Φ(k,Θ^G)(ξ(k)ξ(k)T+V(k)V(k)Tbbbbbbbbb..bbbbbb-2ξ(k)VT(k)).
It follows from (43) that
(44)S(k)=S(k-1)-[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)Y-(k)T+ΦT(k,Θ^G)×P(k)Φ(k,Θ^G)V(k)V(k)T+2[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)VT(k)-2×[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)ξT(k)+ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)ξ(k)ξ(k)T-2ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)ξ(k)VT(k).
We have
(45)-2[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)ξT(k)=[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)Y-(k)T+[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×ξ(k)ξ(k)T-[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×(Y-(k)+ξ(k))(Y-(k)+ξ(k))T,
and the use of the relation
(46)1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)=[1+ΦT(k,Θ^G)P(k-1)Φ(k,Θ^G)]-1⩾0
leads to
(47)S(k)⩽S(k-1)+ΦT(k,Θ^G)P(k)Φ×(k,Θ^G)V(k)V(k)T+2×[1-ΦT(k,Θ^G)P(k)Φ(k,Θ^G)]×Y-(k)VT(k)+ξ(k)ξ(k)T-2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)ξ(k)VT(k).
SinceY-(k), ΦT(k,Θ^G)P(k)Φ(k,Θ^G), ξ(k), and S(k-1) are uncorrelated with V(k), taking the conditional expectation on both sides of (47) and using (A3) and (A4), we obtain
(48)E[S(k)]⩽S(k-1)+2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)Σ+β.
Define
(49)W(k):=S(k)[Ln|P-1(k)|]c.
Since[Ln|P-1(k)|]c is nondecreasing, we have
(50)E[W(k)]⩽S(k-1)[Ln|P-1(k)|]c+2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)[Ln|P-1(k)|]cΣ+β[Ln|P-1(k)|]c⩽W(k-1)+2ΦT(k,Θ^G)P(k)Φ(k,Θ^G)[Ln|P-1(k)|]cΣ+β[Ln|P-1(k)|]c.
Using this property∑i=1∞‍(ΦT(i,Θ^G)P(i)Φ(i,Θ^G)/[Ln|P-1(i)|]c)<∞ (the proof in the same way as the proof of Lemma 1 in [30]), we can see that the sum of the right-hand second term of equation (50) for k from 1 to ∞ is finite. Applying the martingale convergence theorem [30] to the previous inequality, we conclude that W(k) converges a.s. to a finite random variable, say W0; that is,
(51)W(k)=S(k)[Ln|P-1(k)|]c⟶W0<∞
is equivalent to
(52)S(k)=O([Ln|P-1(k)|]c).
From the definition of V(k), we obtain
(53)∥Θ-G∥⩽Θ-G(k)TP-1(k)Θ-G(k)λmin[P-1(k)]=S(k)λmin[P-1(k)].
Now, let us define the matrix trace(54)r(k)=tr[P-1(k)].
It follows that:
(55)|P-1(k)|⩽rn(k),r(k)⩽nλmax|P-1(k)|,Ln|P-1(k)|=O(Lnr(k))=O(Lnλmax|P-1(k)|).
We obtain, finally,(56)∥Θ^G-ΘG∥2=O([Ln|P-1(k)|]cλmin|P-1(k)|)=O([Lnr(k)]cλmin|P-1(k)|)=O([λmax|P-1(k)|]cλmin|P-1(k)|).


## 4.2. Lemma

For the estimate (22) with the assumption (A4), the following proprieties hold.
(P1)
Θ
^
G
is an unbiased estimate of
Θ
G
.
(P2)
The covariance matrix of
Θ
^
G
is given by
(57)
E
[
(
Θ
^
G
-
Θ
G
)
(
Θ
^
G
-
Θ
G
)
T
]
=
(
∑
i
=
0
k
‍
Φ
(
i
,
Θ
G
)
Φ
(
i
,
Θ
G
)
T
)
-
1
Σ
,
where
(58)
Σ
=
(
σ
1
2
0
⋯
0
0
σ
2
2
⋱
⋮
⋮
⋱
⋱
0
0
⋯
0
σ
p
2
)
.Proof.
If we replace (8) in (22), we have
(59)Θ^G=[∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)[ΘTφ(i,D)+V(i)+ψ].
Then,
(60)Θ^G=[∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)[ΘTφ(i,D)+ψ]+[∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)V(i).
So
(61)E[Θ^G]=ΘG+E([∑i=0k‍Φ(i,ΘG)ΦT(i,ΘG)]-1×∑i=0k‍Φ(i,ΘG)V(i)).
Since V(i) is uncorrelated with the elements of Φ(i,ΘG) (13), then
(62)E[Θ^G]=ΘG
which proves (P1).
Consider the first-order Taylor series expansion around the real matrix ofΘG:
(63)∂J(k,Θ^G)∂Θ^G=∂J(k,ΘG)∂ΘG+∂2J(k,ΘG)∂ΘG2(Θ^G-ΘG).
Since∂J(k,Θ^G)/∂Θ^G=0, it derives from (63) that
(64)(Θ^G-ΘG)(Θ^G-ΘG)T=[∂2J(k,ΘG)∂ΘG2]-1∂J(k,ΘG)∂ΘG×[∂J(k,ΘG)∂ΘG]T[(∂2J(k,ΘG)∂ΘG2)-1]T.
The second partial derivative of the criterion with respect to the generalized matrix is(65)∂2J(k,ΘG)∂ΘG2=∑i=0k‍(Υ(i)∂2Υ(i)∂2ΘG-∂Υ(i)∂ΘGΦ(i,ΘG)).
So
(66)∂2J(k,ΘG)∂ΘG2=∑i=0k‍(Υ(i)∂2Υ(i)∂2ΘG+Φ(i,ΘG)ΦT(i,ΘG)).
The use of the small residual algorithms [31] leads to neglect the following term, then:
(67)∑i=0k‍Υ(i)∂2Υ(i)∂2ΘG⟶0.
Hence, an approach of ∂2J(k,ΘG)/∂ΘG2 is obtained:
(68)∂2J(k,ΘG)∂ΘG2≃∑i=0k‍(Φ(i,ΘG)ΦT(i,ΘG)).
Applying the mean value of (∂J(k,ΘG)/∂ΘG)(∂J(k,ΘG)/∂ΘGT), we get:
(69)E[∂J(k,ΘG)∂ΘG∂J(k,ΘG)∂ΘGT]=∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)TE(Υ(i)Υ(i)T).
So
(70)E[∂J(k,ΘG)∂ΘG∂J(k,ΘG)∂ΘGT]=∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)TΣ.
Then, we have
(71)E[(Θ^G-ΘG)(Θ^G-ΘG)T]=E[(∂2J∂ΘG2)-1∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)TΣ[(∂2J∂ΘG2)-1]T].
Finally, we obtain
(72)E[(Θ^G-ΘG)(Θ^G-ΘG)T]=(∑i=0k‍Φ(i,ΘG)Φ(i,ΘG)T)-1Σ
which proves (P2).


## 5. Results

We now present a simulation example and an experimental validation to illustrate the performance of the proposed approach for the simultaneous identification of time delays and parameter matrices of square multivariable systems.


### 5.1. Simulation Example

The objective of the simulation is to compare the efficiency of the proposed method (DRLS) with that of the classic recursive least square approach (RLS) [29] which assumes that the delays are a priori known. In fact, we consider the following cases.Case 1.
The output is noise-free and the RLS method uses the true time delays.Case 2.
The output is noise free and the RLS method uses the misestimated time delays.Case 3.
The output is contaminated by additive noise and the RLS method uses the true time delays.We consider a square linear multivariable discrete time delay system with two inputs and two outputs described by the following equation [32]:
(73)A(q-1)Y(k)=B(q-1)U~+V(k),
where the delayed inputs and the outputs are defined, respectively, by U~(k)=(u1(k-d1)u2(k-d2)) and Y(k)=(y1(k)y2(k)).The two polynomials matricesA(q-1) and B(q-1) are given by
(74)A(q-1)=(1001)+(0.9048000.9048)q-1,B(q-1)=(0.095160.03807-0.02970.0475)q-1.The time delay matrixD is given by
(75)D=(3001).


#### 5.1.1. Case1

The proposed approach (DLSR) and the RLS algorithm are applied to estimate time delays and parameter matrices. The estimation starts with zero initial conditions. The obtained results are illustrated in Table1 and Figures 1, 2, and 3.Table 1
Simulation results: Case1.
True
Developed RLS
RLS
A
1
0.9048
0
0.9
-
0.0022
0.9048
-
0.0000
0
0.9048
-
0.0043
0.9021
-
0.0000
0.9048
B
1
0.09516
0.03807
0.0947
0.0378
0.0952
0.0381
-
0.02974
0.04758
-
0.0297
0.0475
-
0.0297
0.0476
D
3
0
3
0
Known
0
1
0
1Evolution of the true (-)  and the estimated  (- -)  delays: Case1.
(a)
(b)Evolution of the true and the estimated parameters:Ai: Case 1.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 1.
(a)
(b)
(c)
(d)Figures1–3 show the evolution of the estimated and the true parameter matrices.A validation of the obtained model is presented in Figures4 and 5 which show that the estimated outputs track fast and accurately the true outputs.Figure 4
The evolution of the true and the estimated outputsy1(k): Case 1.Figure 5
The evolution of the true and the estimated outputsy2(k): Case 1.


#### 5.1.2. Case2

We apply the proposed approach (DLSR) and the RLS algorithm to estimate time delays and parameter matrices. The RLS algorithm uses misestimated time delays(2001). The obtained results are illustrated in Table 2 and Figures 6 and 7.Table 2
Simulation results: Case2.
True
Developed RLS
RLS
A
1
0.9048
0
0.9
-
0.0022
0.5829
-
0.2376
0
0
.
9048
-
0.0043
0.9021
-
0.1288
-
0.1288
B
1
0.09516
0.03807
0.0947
0.0378
0.0692
0.0277
-
0.02974
0.04758
-
0.0297
0.0475
0.0277
0.0463
D
3
0
3
0
2
0
0
1
0
1
0
1Evolution of the true and the estimated parameters:Ai: Case 2.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 2.
(a)
(b)
(c)
(d)Figures6 and 7 show the evolution of the estimated and the true parameter matrices.


#### 5.1.3. Case3

The system's output is corrupted by additive zero mean white noisesV(k)=[v1(k),v2(k)] with variances σ1=0.0737, σ2=0.086.The result of the simulation is given in Table3.Table 3
Simulation results: Case3.
True
Developed RLS
RLS
A
1
0.9048
0
0.8949
-
0.0084
0.9032
-
0.0016
0
0
.
9048
-
0.0084
0.9
0.0003
0.951
B
1
0.09516
0.03807
0.0941
0.073
0.0952
0.0381
-
0.02974
0.04758
-
0.0303
0.0476
-
0.0297
0.0476
D
3
0
3
0
Known
0
1
0
1Figures8, 9, and 10 show the evolution of the estimated and the true parameter matrices.Evolution of the true  (-)  and the estimated  (- -)  delays: Case3.
(a)
(b)Evolution of the true and the estimated parameters:Ai: Case 3.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 3.
(a)
(b)
(c)
(d)A validation of the obtained model is presented in Figures11 and 12 which show that the estimated outputs track fast and accurately the true outputs.Figure 11
The evolution of the true and the estimated outputsy1(k): Case 3.Figure 12
The evolution of the true and the estimated outputsy2(k): Case 3.


#### 5.1.4. Observations

Based on Tables1–3 and Figures 1–12, we observe that (i)
the RLS method gives the better performance when the true time delays are used. However, it poorly performs for misestimated delays;(ii)
the proposed approach converges to the true delays with acceptable speed for the considered cases.


### 5.2. Experiment Example

The experimental data from a mechanical construction of a CD player arm is considered. The system has two inputs that are forces of the mechanical actuators (U) and two outputs that are related to the tracking accuracy of the arm (Y).The data set contains2048 sample points out of which 1000 were used for the identification procedure and the rest to validate the identified models. The input/output identification signals [33] are given in Figures 13 and 14.Evolution of input identification signals of CD arm.
(a)
(b)Evolution of output identification signals of CD arm.
(a)
(b)The system is described by the following equation:(76)A(q-1)Y(k)=B(q-1)U~(k)+V(k),
where Y(k) and U~(k) are the output and the delayed input of the system at time k, p=2 and the orders of A(q-1), B(q-1) are na=nb=2.The estimation starts with zero initial values for the parameter and the time delay matrices. Applying the proposed algorithm, we obtain(77)A^(q-1)=(1001)+(-0.9464-0.6600-0.1263-0.5278)q-1+(0.08090.5506-0.10660.2243)q-2,B^(q-1)=(-1.48391.02740.2900-1.1819)q-1+(1.2452-0.7638-0.18660.8853)q-2.The estimated time delay matrix is(78)D^=(0002).
Another data set is used for the validation test which is illustrated by Figures 15 and 16.Figure 15
The evolution of the true  (-)  and the estimated  (- -)  outputs.Figure 16
The evolution of the true  (-)  and the estimated  (- -)  outputs.We can see clearly that the estimated output tracks fast and accurately the true output.


## 5.1. Simulation Example

The objective of the simulation is to compare the efficiency of the proposed method (DRLS) with that of the classic recursive least square approach (RLS) [29] which assumes that the delays are a priori known. In fact, we consider the following cases.Case 1.
The output is noise-free and the RLS method uses the true time delays.Case 2.
The output is noise free and the RLS method uses the misestimated time delays.Case 3.
The output is contaminated by additive noise and the RLS method uses the true time delays.We consider a square linear multivariable discrete time delay system with two inputs and two outputs described by the following equation [32]:
(73)A(q-1)Y(k)=B(q-1)U~+V(k),
where the delayed inputs and the outputs are defined, respectively, by U~(k)=(u1(k-d1)u2(k-d2)) and Y(k)=(y1(k)y2(k)).The two polynomials matricesA(q-1) and B(q-1) are given by
(74)A(q-1)=(1001)+(0.9048000.9048)q-1,B(q-1)=(0.095160.03807-0.02970.0475)q-1.The time delay matrixD is given by
(75)D=(3001).


### 5.1.1. Case1

The proposed approach (DLSR) and the RLS algorithm are applied to estimate time delays and parameter matrices. The estimation starts with zero initial conditions. The obtained results are illustrated in Table1 and Figures 1, 2, and 3.Table 1
Simulation results: Case1.
True
Developed RLS
RLS
A
1
0.9048
0
0.9
-
0.0022
0.9048
-
0.0000
0
0.9048
-
0.0043
0.9021
-
0.0000
0.9048
B
1
0.09516
0.03807
0.0947
0.0378
0.0952
0.0381
-
0.02974
0.04758
-
0.0297
0.0475
-
0.0297
0.0476
D
3
0
3
0
Known
0
1
0
1Evolution of the true (-)  and the estimated  (- -)  delays: Case1.
(a)
(b)Evolution of the true and the estimated parameters:Ai: Case 1.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 1.
(a)
(b)
(c)
(d)Figures1–3 show the evolution of the estimated and the true parameter matrices.A validation of the obtained model is presented in Figures4 and 5 which show that the estimated outputs track fast and accurately the true outputs.Figure 4
The evolution of the true and the estimated outputsy1(k): Case 1.Figure 5
The evolution of the true and the estimated outputsy2(k): Case 1.


### 5.1.2. Case2

We apply the proposed approach (DLSR) and the RLS algorithm to estimate time delays and parameter matrices. The RLS algorithm uses misestimated time delays(2001). The obtained results are illustrated in Table 2 and Figures 6 and 7.Table 2
Simulation results: Case2.
True
Developed RLS
RLS
A
1
0.9048
0
0.9
-
0.0022
0.5829
-
0.2376
0
0
.
9048
-
0.0043
0.9021
-
0.1288
-
0.1288
B
1
0.09516
0.03807
0.0947
0.0378
0.0692
0.0277
-
0.02974
0.04758
-
0.0297
0.0475
0.0277
0.0463
D
3
0
3
0
2
0
0
1
0
1
0
1Evolution of the true and the estimated parameters:Ai: Case 2.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 2.
(a)
(b)
(c)
(d)Figures6 and 7 show the evolution of the estimated and the true parameter matrices.


### 5.1.3. Case3

The system's output is corrupted by additive zero mean white noisesV(k)=[v1(k),v2(k)] with variances σ1=0.0737, σ2=0.086.The result of the simulation is given in Table3.Table 3
Simulation results: Case3.
True
Developed RLS
RLS
A
1
0.9048
0
0.8949
-
0.0084
0.9032
-
0.0016
0
0
.
9048
-
0.0084
0.9
0.0003
0.951
B
1
0.09516
0.03807
0.0941
0.073
0.0952
0.0381
-
0.02974
0.04758
-
0.0303
0.0476
-
0.0297
0.0476
D
3
0
3
0
Known
0
1
0
1Figures8, 9, and 10 show the evolution of the estimated and the true parameter matrices.Evolution of the true  (-)  and the estimated  (- -)  delays: Case3.
(a)
(b)Evolution of the true and the estimated parameters:Ai: Case 3.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 3.
(a)
(b)
(c)
(d)A validation of the obtained model is presented in Figures11 and 12 which show that the estimated outputs track fast and accurately the true outputs.Figure 11
The evolution of the true and the estimated outputsy1(k): Case 3.Figure 12
The evolution of the true and the estimated outputsy2(k): Case 3.


### 5.1.4. Observations

Based on Tables1–3 and Figures 1–12, we observe that (i)
the RLS method gives the better performance when the true time delays are used. However, it poorly performs for misestimated delays;(ii)
the proposed approach converges to the true delays with acceptable speed for the considered cases.


## 5.1.1. Case1

The proposed approach (DLSR) and the RLS algorithm are applied to estimate time delays and parameter matrices. The estimation starts with zero initial conditions. The obtained results are illustrated in Table1 and Figures 1, 2, and 3.Table 1
Simulation results: Case1.
True
Developed RLS
RLS
A
1
0.9048
0
0.9
-
0.0022
0.9048
-
0.0000
0
0.9048
-
0.0043
0.9021
-
0.0000
0.9048
B
1
0.09516
0.03807
0.0947
0.0378
0.0952
0.0381
-
0.02974
0.04758
-
0.0297
0.0475
-
0.0297
0.0476
D
3
0
3
0
Known
0
1
0
1Evolution of the true (-)  and the estimated  (- -)  delays: Case1.
(a)
(b)Evolution of the true and the estimated parameters:Ai: Case 1.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 1.
(a)
(b)
(c)
(d)Figures1–3 show the evolution of the estimated and the true parameter matrices.A validation of the obtained model is presented in Figures4 and 5 which show that the estimated outputs track fast and accurately the true outputs.Figure 4
The evolution of the true and the estimated outputsy1(k): Case 1.Figure 5
The evolution of the true and the estimated outputsy2(k): Case 1.


## 5.1.2. Case2

We apply the proposed approach (DLSR) and the RLS algorithm to estimate time delays and parameter matrices. The RLS algorithm uses misestimated time delays(2001). The obtained results are illustrated in Table 2 and Figures 6 and 7.Table 2
Simulation results: Case2.
True
Developed RLS
RLS
A
1
0.9048
0
0.9
-
0.0022
0.5829
-
0.2376
0
0
.
9048
-
0.0043
0.9021
-
0.1288
-
0.1288
B
1
0.09516
0.03807
0.0947
0.0378
0.0692
0.0277
-
0.02974
0.04758
-
0.0297
0.0475
0.0277
0.0463
D
3
0
3
0
2
0
0
1
0
1
0
1Evolution of the true and the estimated parameters:Ai: Case 2.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 2.
(a)
(b)
(c)
(d)Figures6 and 7 show the evolution of the estimated and the true parameter matrices.


## 5.1.3. Case3

The system's output is corrupted by additive zero mean white noisesV(k)=[v1(k),v2(k)] with variances σ1=0.0737, σ2=0.086.The result of the simulation is given in Table3.Table 3
Simulation results: Case3.
True
Developed RLS
RLS
A
1
0.9048
0
0.8949
-
0.0084
0.9032
-
0.0016
0
0
.
9048
-
0.0084
0.9
0.0003
0.951
B
1
0.09516
0.03807
0.0941
0.073
0.0952
0.0381
-
0.02974
0.04758
-
0.0303
0.0476
-
0.0297
0.0476
D
3
0
3
0
Known
0
1
0
1Figures8, 9, and 10 show the evolution of the estimated and the true parameter matrices.Evolution of the true  (-)  and the estimated  (- -)  delays: Case3.
(a)
(b)Evolution of the true and the estimated parameters:Ai: Case 3.
(a)
(b)
(c)
(d)Evolution of the true and the estimated parameters:Bi: Case 3.
(a)
(b)
(c)
(d)A validation of the obtained model is presented in Figures11 and 12 which show that the estimated outputs track fast and accurately the true outputs.Figure 11
The evolution of the true and the estimated outputsy1(k): Case 3.Figure 12
The evolution of the true and the estimated outputsy2(k): Case 3.


## 5.1.4. Observations

Based on Tables1–3 and Figures 1–12, we observe that (i)
the RLS method gives the better performance when the true time delays are used. However, it poorly performs for misestimated delays;(ii)
the proposed approach converges to the true delays with acceptable speed for the considered cases.


## 5.2. Experiment Example

The experimental data from a mechanical construction of a CD player arm is considered. The system has two inputs that are forces of the mechanical actuators (U) and two outputs that are related to the tracking accuracy of the arm (Y).The data set contains2048 sample points out of which 1000 were used for the identification procedure and the rest to validate the identified models. The input/output identification signals [33] are given in Figures 13 and 14.Evolution of input identification signals of CD arm.
(a)
(b)Evolution of output identification signals of CD arm.
(a)
(b)The system is described by the following equation:(76)A(q-1)Y(k)=B(q-1)U~(k)+V(k),
where Y(k) and U~(k) are the output and the delayed input of the system at time k, p=2 and the orders of A(q-1), B(q-1) are na=nb=2.The estimation starts with zero initial values for the parameter and the time delay matrices. Applying the proposed algorithm, we obtain(77)A^(q-1)=(1001)+(-0.9464-0.6600-0.1263-0.5278)q-1+(0.08090.5506-0.10660.2243)q-2,B^(q-1)=(-1.48391.02740.2900-1.1819)q-1+(1.2452-0.7638-0.18660.8853)q-2.The estimated time delay matrix is(78)D^=(0002).
Another data set is used for the validation test which is illustrated by Figures 15 and 16.Figure 15
The evolution of the true  (-)  and the estimated  (- -)  outputs.Figure 16
The evolution of the true  (-)  and the estimated  (- -)  outputs.We can see clearly that the estimated output tracks fast and accurately the true output.


## 6. Conclusions

In this paper, we have addressed the problem of identification of linear discrete time delay multivariable systems. In fact, we have proposed a novel approach for the simultaneous identification of the unknown time delays and the parameter matrices of these systems. The proposed approach consists in constructing a linear-parameter formulation that is used to estimate the time delays and the polynomial matrices using recursive least square algorithm. The obtained estimates were shown to be unbiased, and an expression for their covariance matrix was given. Numerical simulation and experimental test are presented to demonstrate the performance of the proposed approach.



---

*Source: 658194-2013-09-16.xml*