A CS Recovery Algorithm for Model and Time Delay Identification of MISO-FIR Systems

This paper considers identifying the multiple input single output finite impulse response (MISO-FIR) systems with unknown time delays and orders. Generally, parameters, orders and time delays of an MISO system are separately identified from different algorithms. In this paper, we aim to perform the model identification and time delay estimation simultaneously from a limited number of observations. For an MISO-FIR system with many inputs and unknown input time delays, the corresponding identification model contains a large number of parameters, requiring a great number of observations for identification and leading to a heavy computational burden. Inspired by the compressed sensing (CS) recovery theory, a threshold orthogonal matching pursuit algorithm (TH-OMP) is presented to simultaneously identify the parameters, the orders and the time delays of the MISO-FIR systems. The proposed algorithm requires only a small number of sampled data compared to the conventional identification methods, such as the least squares method. The effectiveness of the proposed algorithm is verified by simulation results.


Introduction
In many industry processes, time delay is an unavoidable behavior in their dynamical models.For example, in a chemical industry, analyzers need a certain amount of time to process an analysis, which is used as a part of the control loop [1].Analyzing and designing a controller for a system with time delays depend on an effective dynamical system model.System identification is an effective way to build mathematical models of dynamical systems from observed input-output data.Conventional identification methods, such as the least squares methods [2], the stochastic gradient methods [3] and the maximum likelihood estimation methods [4], usually require a large number of sampled data and that the model structures, including time delays, are a priori known [5,6].However, this is not the case in many situations in practice for which the time delays and systems orders have to be estimated together with the parameters.Most of the process industries, such as distillation columns, heat exchangers and reactors, can be modeled as multiple input multiple output (MIMO) or multiple input single output (MISO) systems with long time delays [1].For a multiple input single output finite impulse response (MISO-FIR) system with many inputs and unknown input time delays, the parameterized model contains a large number of parameters, requiring a great number of observations and a heavy computational burden for identification.Furthermore, the measuring cost is an issue.Especially for chemical processes, it will take a long time to get enough data because of the long analysis time.Therefore, it is a challenging work to find new methods to identify the time delays and parameters of multivariable systems with less computational burden using a small number of sampled data.
System models of practical interest, especially the control-oriented models, often require low order, which implies that the systems contain only a few non-zero parameters.Because the time delays and orders are unknown, we can parameterize such an MISO-FIR system by a high-dimensional parameter vector, with only a few non-zero elements.These systems can be termed as sparse systems.On the other hand, a sparse MISO-FIR system contains a long impulse response, but only a few terms are non-zeros.The sparse systems are widely discussed in many fields [7][8][9][10][11].For such sparse systems, the identification objective is to estimate the sparse parameter vector and to read the time delays and orders from the structure of the estimated parameter vector.
In this paper, we consider the identification of the sparse MISO-FIR systems by using the compressed sensing (CS) recovery theory [12][13][14][15].CS has been widely applied in signal processing and communication systems, which enables the recovery of an unknown vector from a set of measurements under the assumption that the signal is sparse and certain conditions on the measurement matrices [16][17][18][19].The identification of a sparse system can be viewed as a CS recovery problem.The sparse recovery problem can be solved by the l 1 -norm convex relaxation.However, the l 1 regularization schemes are always complex, requiring heavy computational burdens [20][21][22].The greedy methods have speed advantages and are easily implemented.In this literature, a block orthogonal matching pursuit (BOMP) algorithm has been applied to identify the determined linear time invariant (LTI) and linear time variant (LTV) auto-regressive with external (ARX) models without disturbances [9].In this paper, we will consider identifying the parameters, orders and time delays of the MISO-FIR systems from a small number of noisy measured data.
Briefly, the structure of this paper is as follows.Section 2 gives the model description of MISO-FIR systems and describes the identification problem.Section 3 presents a threshold orthogonal matching pursuit (TH-OMP) algorithm based on the compressed sensing theory.Section 4 provides a simulation example to show the effectiveness of the proposed algorithm.Finally, some concluding remarks are given in Section 5.

Problem Description
Consider an MISO-FIR system with time delays: where u i (t) is the i-th input, y(t) the output, d i the time delay of the i-th input channel and v(t) a white noise with zero mean and variance σ 2 , B i (z), i = 1, 2, • • • , r are the polynomials in the unit backward shift operator z −1 (i.e., z −1 y(t) = y(t − 1)) and: Assume that y(t) = 0, u(t) = 0 and v(t) = 0 for t 0 and the orders n bi , the time delays d i , as well as the parameters b ij are unknown.Thus, the identification goal is to estimate the unknown parameters b ij , orders n bi and the time delays d i from the observations.Define the information vector: where l is the input maximum regression length.l must be chosen to capture all of the time delays and the degrees of the polynomials B i (z).Thus, l max(d i + n bi ).The dimension of ϕ(t) is N = lr.The corresponding parameter vector θ is: Then, Equation ( 1) can be written as: Because of the unknown orders and time delays, the parameter vector θ contains many zeros, but only a few non-zeros.According to the CS theory, the parameter vector θ can be viewed as a sparse signal, and the number of non-zero entries K = r i=1 n bi is the sparsity level of θ [8].Let θ 0 denote the number of non-zero entries in θ.Then, the identification problem can be described as: where θ is the estimate of θ and ε > 0 is the error tolerance.
For m observations, define the stacked matrix and vectors as: Then, we have: If there are enough measurements, i.e., m N , according to the least squares principle, we can get the least squares estimate of θ, However, from Equation (2), it is obvious that N is a large number, especially for systems with a large number of inputs.Therefore, it will take a lot of time and effort to get enough measurements.In order to improve the identification efficiency, in this paper, we aim to identify the parameter vector θ using less observations (K < m < N ) based on the CS recovery theory.

Identification Algorithm
The CS theory, which was introduced by Candès, Romberg and Tao [12,13] and Donoho [14], has been widely used in signal processing.It indicates that an unknown signal vector can be recovered from a small set of measurements under the assumptions that the signal vector is sparse and some conditions on the measurement matrix.The identification problem can be viewed as a CS recovery problem in that the parameter vector θ is sparse.
From Equation (4), we can see that the output vector y m can be written as a linear combination of the columns of Φ m plus the noise vector v m , i.e., where θ i is the i-th element of θ and φ i is the i-th column of Φ m .Because there are only K non-zero parameters in θ, the main idea is to find the K non-zero items at the right-hand side of Equation (5).
In this paper, we propose a TH-OMP algorithm to do the identification.First, we give some notations.Let k = 1, 2, • • • be the iterative number and λ k be the solution support of the k-th iterative number; Λ k is a set composed of λ i , i = 1, 2, • • • , k; r k denotes the residual of the k-th iterative number; Φ Λ k is the sub-information matrix composed by the k columns of Φ m indexed by Λ k .θk is the parameter estimation of the k-th iterative number.The TH-OMP algorithm is initialized as r 0 = y m , Λ 0 = ∅ and θΛ 0 = 0.At the k-th iteration, minimizing the criteria function: with respect to θ i yields: Plugging it back into ε(i), we have: From the above equation, we can see that the quest for the smallest error is equivalent to the quest for the largest inner product between the residual r k−1 and the normalized column vectors of Φ m .Thus, the k-th solution support can be obtained by: Update the support set Λ k and the sub-information matrix Φ Λ k by: The next step is to find the k-th provisional parameter estimate from the sub-information matrix Φ Λ k and measurement vector y m .Minimizing the cost function: leads to the k-th parameter estimate: Then, the k-th residual can be computed by: If the system is noise free, then the residual will be zero after K iterations, and the non-zero parameters, as well as their positions can be exactly determined.This is the so-called orthogonal matching pursuit (OMP) algorithm.However, for systems disturbed by noises, the OMP algorithm cannot give the sparsest solution, because the residual r K is non-zero after K iterations, and the iteration continues.In order to solve this problem, an appropriate small threshold can be set to filter the parameter estimate θΛ k .If | θΛ | < , let θΛ = 0, where θΛ is the i-th element of θΛ k .The choice of will be discussed in the simulation part.Denote θΛ k as the parameter estimate after the filtering.Then, the residual can be computed by: if r k < ε, then stop the iteration; otherwise, go to apply another iteration.If the iteration stops at k, we can recover the parameter estimate θ ∈ R N from θΛ k ∈ R k based on the support set Λ k .
It is worth noting that taking the derivative of J(θ Λ k ) with respect to θΛ k gives: The above equation shows that the residual r k is orthogonal to the sub-information matrix Φ Λ k .Therefore, it does not require computing the inner products of r k and the columns in Φ Λ k during the (k + 1)-th iteration.Modify Equation ( 7) as: Then, the computational burden can be reduced.
From the recovered parameter estimation θ, the orders and time delays of each input channel can be estimated.The orders n bi can be determined by the number of elements of each non-zero block.Then, the time delay estimates can be computed from the orders n bi and the input maximum regression length l.It can be seen from Equation (2) that θ contains (r + 1) zero blocks.Denote the number of zeros of each zero-block as n i , i = 1, 2, • • • , r + 1.Then, the time delay estimates can be computed by: The proposed TH-OMP algorithm can be implemented as the following steps.
3. Update the support set Λ k and the sub-information matrix Φ Λ k using Equation (8). 4. Compute the parameter estimate θΛ k by using Equation (10). 5. Using the threshold to filter θΛ k to get θΛ k .6. Compute the residual using Equation (12).If r k < ε, stop the iteration; otherwise, go back to Step 2 to apply another iteration.7. Recover the parameter estimate θ ∈ R N from θΛ k ∈ R k based on the support set Λ k .8. Read the orders n bi from θ ∈ R N and estimate the time delays using Equation (15).

Simulation Example
Consider the steam water heat exchanger system in Figure 1, where the steam flow u 1 and the water flow u 2 are two manipulated variables to control the water temperature y.The discrete input-output representation of the process is expressed as the following equation.
Let l = 50; then the parameter vector to be identified is: θ = [0 42 , 0.95, 0.72, 0.40, 0 25 , −0.63, −0.47, −0.25, 0 27 ] T ∈ R N It can be seen from the above equation that the dimension of the parameter vector θ is N = 100, but only six of the parameters are non-zeros.
In the identification, the inputs {u i (t)}, i = 1, 2 are taken as independent persistent excitation signal sequences with zero mean and unit variances and {v(t)} as white noise sequences with zero mean and variance σ 2 = 0.10 2 .Let = 0.05, applying the proposed TH-OMP algorithm and the least squares (LS) algorithm to estimate the system with m measurements.The parameter estimation errors δ := ( θ − θ )/ θ versus m are shown in Figure 2.
From Figure 2, we can see that under a limited dataset (m < 100), the parameter estimation error of the TH-OMP algorithm approaches zero, and the estimation accuracy is much higher than that of the LS algorithm.
Let m = 80; using the TH-OMP algorithm to estimate θ, the parameter estimation errors δ versus iteration k are shown in Figure 3; the non-zero parameter estimates versus iteration k are shown in Table 1 and Figure 4, where the dash dot lines are true parameters and the solid lines are parameter estimates.
If we choose = 0.01, then we can get the estimated parameter vector as: It is obvious that the parameter vector θ has only three zero blocks in this example.However, we can see from the above equation that there exist three undesirable parameter estimates θ5 , θ31 and θ87 .Moreover, the three parameter estimates are much smaller than other parameters.Therefore, according to the structure of θ, the parameters θ5 , θ31 and θ87 can be set to zeros.Then, the estimation result is the same as the case when = 0.05.This implies that even if a too small is chosen, we can still obtain the effective parameter estimates based on the model structures.
From Equations ( 15) and ( 16), we can obtain the orders n b 1 = n b 2 = 3 and get the time delay estimates as: From Figures 2-4, Table 1 and Equations ( 16) and ( 17), we can conclude that the TH-OMP algorithm converges faster than the LS algorithm, can achieve a high estimation accuracy by K iterations with finite measurements (m < N ) and can effectively estimate the time delay of each channel.

Conclusions
This paper presents a TH-OMP algorithm for MISO-FIR systems with unknown orders and time delays.The parameter estimates, orders and time delays can be simultaneously estimated from a small number of observation data.The proposed method can be simply implemented and can reduce the measuring cost, as well as improve the identification efficiency.The simulation results are given to demonstrate the performance of the proposed method.

Figure 4 .
Figure 4.The parameters estimates versus iteration k.