Iterative Identification for Multivariable Systems with Time-Delays Based on Basis Pursuit DeNoising and Auxiliary Model

This paper focuses on the joint estimation of parameters and time-delays of the multiple-input single-output output-error systems. Since the time-delays are unknown, an effective identification model with a high dimensional and sparse parameter vector is established based on overparameterization. Then, the identification problem is converted to a sparse optimization problem. Based on the basis pursuit de-noising criterion and the auxiliary model identification idea, an auxiliary model based basis pursuit de-noising iterative algorithm is presented. The parameters are estimated by solving a quadratic program, and the unavailable terms in the information vector are updated by the auxiliary model outputs iteratively. The time-delays are estimated according to the sparse structure of the parameter vector. The proposed method can obtain effective estimates of the parameters and time-delays from few sampled data. The simulation results illustrate the effectiveness of the proposed algorithm.


Introduction
This paper focuses on the identification of the multiple-input single-output output-error systems with unknown time-delays.The investigation is introduced from its background, the formulation of the problem, the literature survey, and the scope and the contributions in this section.

Background
System identification is a process of developing mathematical models of dynamical systems based on observed data [1].The dynamical system to be modeled can be a mechanical system, an industrial system or other types of system.Many identification methods have been proposed for the specific systems.For example, some identification techniques for mechanical systems can be found in [2][3][4].
In this paper, we study the identification of multivariable systems because most of the industrial processes can be modeled as multivariable systems.
Multivariable systems are a class of systems which often present significant interactions between the multiple-input and multiple-output channels.A wealth of literature is available on effective identification algorithms for such systems.For example, an expectation-maximization technique was employed for the estimation of the linear time-invariant (LTI) multiple-input multiple-output (MIMO) systems [5].A matchable-observable linear identification with zero-order oracle filter tuning method was developed for the LTI MIMO systems by comprising the problems of parameter estimation, filter tuning and model structure selection [6].The methods proposed in [5,6] are applicable for the MIMO systems without time-delays.In many industry processes, time-delays are usually unavoidable.Therefore, it is of great interest to study the identification of multivariable systems with unknown time-delays.
For the estimation of time-delays in multivariable systems, several methods have been proposed.In [7], a parametric model is first estimated by using generalized orthonormal basis filters.Then, the time-delay estimates are generated by analyzing the simulated step response of the resulting model.However, it is required that the model structure is available a priori.To address this problem, a non-parametric method was proposed in [8].The central idea is to break up a general m × n MIMO system into mn decoupled single-input single-output systems using partial coherence functions, and then to apply the frequency-domain method to estimate the time-delay of each subsystem.Since the time-delay estimation procedure was repeated for mn times, this method requires a large amount of computation and is limited for high order systems.In this paper, we attempt to find a simple and efficient method to address the joint estimation of time-delays and parameters of multivariable systems.

Formulation of the Problem of Interest for this Investigation
For simplicity, we consider the identification problem of multiple-input single-output output-error (MISO-OE) systems with unknown input time-delays.The main body of the MISO-OE model is the impulse transfer functions between the inputs and output.The identification challenge lies in that the system output is nonlinear in the polynomial parameters.
For the identification of the multivariable systems with time-delays, it is needed to form an appropriate identification model first.Taking into account the unknown time-delays, a high dimensional identification model is formed with a sparse parameter vector, which contains many zeros [9,10].For the identification of high dimensional models, traditional identification methods, such as the least squares method and the stochastic gradient method, require many observations as well as heavy computational burden.However, there are many situations in reality where only limited observations are available such as setpoint-operated processes, online estimation and linear time-variant system identification [11].It has proven that traditional identification is unreliable in these cases and the estimates often tend to be biased or have large variances [12].Therefore, it is necessary to develop new identification methods and reduce the identification cost.
The compressive sensing theory originally arose in the signal processing community, which enables the reconstruction of sparse or compressible signals under the undersampling condition [13].In view of the sparse characteristic of the identification model, the estimation of the sparse parameter vector can be addressed as a reconstruction problem of the sparse signal.

Literature Survey
Greedy algorithms and convex optimization approaches are most widely used reconstruction methods.Greedy algorithms are commonly applied for solving minimum l 0 norm problem with advantages of fast speed and easy implement [14].Convex optimization approaches are a class of global optimization methods, which address minimum l 1 norm problems with high stability and strong applicability [15].Basis pursuit, basis pursuit de-noising (BPDN) [16], the least absolute shrinkage and selection operator (LASSO) [17], and the least angle regression [18] are typical convex optimization approaches.These approaches can be converted to linear programming or quadratic programming form, and then the optimal solutions can be obtained [19].Recently, the reconstruction methods have been applied for system identification.For example, the orthogonal matching pursuit algorithm was applied to obtain the joint estimation of parameters and time-delays of the MISO finite impulse response systems and MISO controlled autoregressive systems [20,21].The compressive sampling matching pursuit algorithm was modified by employing the instrumental variable method to identify a class of closed-loop systems [22].The LASSO approach is investigated for computing efficient model structure of the overparameterized nonlinear autoregressive, moving average exogenous systems [23].In this paper, we aim to identify the MISO-OE systems with time-delays based on a convex optimization approach.
Compared with the systems studied in [20,21], the structure of the MISO-OE systems is more general and complex.Considering that the system output is nonlinear in the polynomial parameters, the MISO-OE system model can be transformed into a linear regression form by introducing intermediate variables.However, the intermediate variables are unmeasurable and make the identification difficult.Fortunately, the auxiliary model identification idea can solve this problem by replacing unmeasured variables with outputs of the auxiliary model which is constructed of measurable variables [24].Combining the auxiliary model idea with conventional methods, many new identification methods have been developed such as the auxiliary model based stochastic gradient algorithm [25], and the auxiliary model based least squares (AM-RLS) algorithm [26].It is indicated that the AM-RLS algorithm can obtain unbiased and consistent estimation of parameters under a persistent excitation condition.Inspired by the algorithms in [25,26], we attempt to combine the auxiliary model idea with the convex optimization approach to deal with the joint estimation of parameters and time-delays of MISO-OE systems.

Scope and Contribution of this Study
This investigation deals with the identification of the MISO-OE systems with unknown time-delays based on convex optimization and auxiliary model.The structure of the MISO-OE systems is different from the systems studied in [20,21], which means that different methods need to be applied to form the identification model.In this paper, the auxiliary model technique is employed to address the nonlinearity of parameters, and an effective identification model with a high dimensional and sparse parameter vector is derived due to the unknown time-delays.In addition, a different identification method is applied here.The methods proposed in [20,21] are modifications of a greedy algorithm.In this paper, a convex optimization approach, the BPDN approach, is modified for identification due to its robustness.By converting the BPDN approach to a quadratic programming form, the parameters are calculated from the optimal solutions, and the unmeasurable terms in the information vector are updated iteratively.The unknown time-delays can be read from the estimated parameter vector.The effectiveness of the proposed algorithm is tested by simulation examples.

Organization of the Paper
The rest of this paper is organized as follows.Section 2 introduces the MISO-OE systems with time-delays and describes the identification problem.Section 3 presents an auxiliary model-basis pursuit de-noising iterative (AM-BPDNI) algorithm.Section 4 provides a simulation example to show the effectiveness of the AM-BPDNI algorithm.Finally, some summaries are given in Section 5.

Problem Description
Consider an MISO-OE system where u i (t) and d i are the input and time-delay of the ith input channel, y(t) is the output, v(t) is a white noise vector with zero mean and variance σ 2 , and A i (z) and B i (z) are the time-invariant polynomials with constant coefficients in the unit backward shift operator Assume that the orders n ai and n bi are known, y(t) = 0, u i (t) = 0 and v(t) = 0 for t < 0.
To form the identification model, an intermediate variable is introduced [24], Since the time-delay d i of each input channel is unknown, an overparameterization method is applied by setting a maximum input regression length l which satisfies l ≥ max(d i + n bi ) [20].Then, x i (t) can be written in an impact form where The identification model of the system in Equation ( 1) can be formed as where It can be seen from Equations ( 3) and ( 6) that the parameter vector θ contains many zeros, therefore θ is a sparse vector and the system in Equation ( 4) is a sparse system.The sparsity level can be measured by K := r ∑ i=1 (n ai + n bi ), which denotes the number of non-zero elements in θ.The identification objective is to estimate the unknown parameters a ij , b ij as well as the time-delays d i from observations.

Identification Algorithm
If we have m observations from t = 1 to t = m, Equation ( 4) can be written in a stacked form where From Equations ( 2), ( 5) and (10), we can see that the information matrix Φ contains many unknown intermediate terms.Therefore, it is difficult to perform the identification directly.According to the auxiliary model identification idea [24,27], the information matrix Φ can be replaced with its estimate where φi,k (t Note that the unmeasurable terms x i (t − j) are replaced with their auxiliary model output estimates xi (t − j).Then, the parameter vector θ can be estimated by the auxiliary model based least squares iterative (AM-LSI) algorithm [28], where θk denotes the parameter vector estimate at the kth iteration.According to the least squares (LS) theory, the AM-LSI algorithm is efficient if it is satisfied that m n.However, from Equation ( 7), we can see that the dimension of the system in Equation ( 8) is high.Therefore, it would take a lot of time and efforts to obtain enough observations to meet the identification requirement.Moreover, Equation (14) shows that the AM-LSI algorithm requires computing the inverse matrix [ ΦT k Φk ] −1 at each iteration, which leads to a heavy computational burden.Furthermore, the sparse solution cannot be obtained [23], and the time-delays cannot be effectively estimated.Thus, the AM-LSI algorithm is infeasible for the high dimensional and sparse system identification.
Inspired by the compressive sensing theory, the identification of the sparse system in Equation ( 8) can be further expressed as an optimization problem where • 0 represents the l 0 norm, • 2 the l 2 norm, and ε the error tolerance, which is a priori chosen.However, it is a non-convex problem and is difficult to solve in practice.A commonly used alternative is to replace the l 0 norm with a relaxed convex l 1 norm [9], where • 1 is the l 1 norm.It is proved that the minimization of l 1 norm is equivalent to the minimization of l 0 norm when the Restricted Isometry Property (RIP) is satisfied [13].Then, the convex optimization problem in Equation ( 17) is solvable with the BPDN criterion [16].The BPDN approach can effectively reduce the interference of noise and has good robustness.Inspired by the auxiliary model and the iterative identification idea in the AM-LSI algorithm, we apply the BPDN criterion to modify the AM-LSI algorithm by replacing the LS step with the BPDN criterion.The parameters are estimated by the BPDN criterion and the unmeasurable terms are replaced by the outputs of the auxiliary model in each iteration.Let k = 1, 2, • • • be the iterative number.To get accurate reconstruction from Equation ( 17), the information matrix should satisfy some conditions, such as the RIP [29] or the exact recovery condition (ERC) [9].The ERC guarantee and the consistency properties for the identification of the controlled autoregressive models have been investigated in [12].To meet the ERC, we normalize the information matrix Φk defined in Equations ( 11)-( 13) by dividing the elements in each column by the l 2 norm of that column [30,31].Denote the element of the ith row and jth column in Φk by Φk,ij , and the normalized information matrix Φk,n is constructed by where Φk (j Similarly, the normalized parameter vector θ k,n can be defined as Note that the location of non-zeros in θ n is identical to that in θ.Accordingly, the constrained optimization problem in Equation ( 17) equals The problem in Equation ( 20) is closely related to the following unconstrained convex optimization problem where λ is a nonnegative parameter.Since the information matrix Φk is normalized, we can set λ to the value λ = σ 2 log(n) [16].
The key step of the BPDN approach is to express Equation ( 21) as a quadratic program (QP) problem [32].To begin with, two nonnegative vectors u n and v n are introduced to express θ n .Let θ nj , u nj and v nj be the jth element of the vectors θ n , u n and v n , respectively, where u nj ≥ 0 and v nj ≥ 0 for all j = 1, 2, • • • , n. Define u nj := (θ nj ) + , v nj := (−θ nj ) + , ( * ) + := max{0, * }.
Then, θ n can be rewritten as Accordingly, the l 1 regularization term θ n 1 can be expressed as where Note that all elements in z n are nonnegative.Similarly, the quadratic error term can be written as Equation ( 24) can be further written as From Equations ( 22) and ( 26), we have = min where Since 1 2 Y T Y is a constant, Equation ( 28) can be constructed in a standard QP framework, min For the inequality constrained QP problem in Equation ( 30), the common solution is the active set method [33].While for the sake of simplicity, the QP problem can be directly solved by calling the relevant function from the standard scientific software toolbox.For example, the MATLAB toolbox provides a function "quadprog".Then, the parameter vector θn can be obtained from the optimum solution ẑn , where ẑn (1 : n) represents a vector formed by the first n elements of ẑn , and ẑn (n + 1 : 2n) a vector formed by the last n elements of ẑn .Considering that the system in Equation ( 8) is contaminated with noise, the parameter estimation error can be large.To further reduce the estimation error, a small threshold TH= can be set to filter the elements close to zero in θn .Let θnj = 0 if θnj < and denote the filtered parameter vector as θn, .Then, the parameter vector estimate θk can be recovered according to Equation (19), The estimates of the intermediate variables x i,k (t) can be refreshed by θk as shown in Equations ( 15) and (16).
The unknown time-delay of each input channel can be estimated according to the location of zero-blocks and the number of zeros in θ .It can be seen from Equations ( 3) and ( 6) that there are 2r zero-blocks in θ.Denote the number of zeros in each zero-block by z i (i = 1, 2, • • • , 2r).Then, time-delays can be estimated by

Simulation Example
Example 1.Consider the following MISO-OE system with time-delays, The system in Equation ( 34) is a second order system with three inputs and one output.The inputs {u i (t)}, i = 1, 2, 3 are taken as random uncorrelated signal sequences with zero mean and unit variances, and {v(t)} as a white noise sequence with zero mean and variances σ 2 .Let the maximum input regression length be l = 50.Then, the parameter vector to be identified is θ 1 = [−0.1,0.7, 0 20 , 1.5, 0.9, 0 28 , 0.3, 0.5, 0 10 , 0.2, 1.8, 0 38 , −0.2, −0.4,0 30 , −0.1, 2, 0 18 ] T ∈ R 156 , where 0 i denotes the zero-block with i zeros.Note that the number of non-zero elements is Taking m = 130 and TH = 0.001, apply the AM-LSI algorithm and the AM-BPDNI algorithm to perform the identification, respectively.The parameter estimation errors δ := θ − θ / θ versus different noise levels are shown in Table 1.When σ 2 = 0.10 2 , the estimation errors δ versus the iterative number k are shown in Figure 1.It can be seen that the AM-BPDNI algorithm performs better than the AM-LSI algorithm and is insensitive to noise.Let the variance of {v(t)} be σ 2 = 0.10 2 .Apply the AM-BPDNI algorithm to obtain the estimated model of the system in Equation (34) with the first m = 130 data.Then, validate the estimated model by using m e = 200 samples from t = 131 to 330.The predicted output of the estimated model, the true output of the system and their errors are plotted in Figure 2. It is shown that the predicted outputs ŷ(t) are close to the true outputs y(t).Moreover, the average predicted output error is small and close to the standard deviation of the noise σ = 0.10.It follows that the estimated model can well capture the system dynamics.Let m = 130 and σ 2 = 0.10 2 .Using the AM-BPDNI algorithm to estimate the sparse parameter vector θ, the non-zero parameter estimates versus k are shown in Table 2 and Figure 3.   Obviously, the time-delay estimates are agreement with the true time-delays.
Example 2. Consider the following MISO-OE system with time-delays, Compared with the system in Equation (34), the system in Equation (36) has one more input.Thus, the number of parameters is increased.Let l = 50 and the true parameter vector is θ 2 = [−0.1,0.7, 0 20 , 1.5, 0.9, 0 28 , 0.3, 0.5, 0 10 , 0. The parameter estimation errors of the systems in Equations ( 34) and (36) versus k are shown in Figure 4.The running time of the proposed method for the systems is t 1 = 1.374637 s and t 2 = 2.398602 s.It can be concluded that the computational burden increases as the dimension of the parameter vector increases.
The simulation results show that for the MISO-OE model, the proposed AM-BPDNI algorithm can obtain efficient estimation of parameters from few observations (m < n) with good robustness.Moreover, the AM-BPDNI algorithm can effectively estimate the time-delays according to the sparse characteristic of the estimated parameter vector.However, as the number of the input channels increases, the computational burden of the proposed algorithm increases.

Conclusions
This paper proposes an AM-BPDNI algorithm for the identification of MISO-OE systems with unknown time-delays.Based on the BPDN criterion and the auxiliary model identification idea, the sparse parameters and multiple-input time-delays can be effectively and simultaneously estimated.The proposed algorithm requires few sampled data and is robust to noise.

Figure 4 .
Figure 4.The parameter estimation errors δ of different systems versus k.