LQ decomposition based subspace identification under deterministic type disturbance

ABSTRACT To overcome the influence from deterministic type load disturbance with unknown dynamics, a bias-eliminated subspace identification method is proposed for consistent estimation. By decomposing the output response into three parts, deterministic, disturbed and stochastic components, in terms of the linear superposition principle, an LQ decomposition approach is developed to eliminate the disturbance and noise effect for unbiased estimation of the deterministic system state. Subsequently, a shift-invariant approach is given to retrieve the state matrices. Consistent estimation on the state matrices is analyzed with a proof. Illustrative example of open-loop system identification is shown to demonstrate the effectiveness and merit of the proposed method.


Introduction
Subspace identification methods (SIMs) have been continuously developed in the last three decades, owing to uniform state-space description of multiple-input multiple-output (MIMO) systems for control system design such as model predictive control (MPC). There are a few well recognized subspace identification methods for practical application, e.g. the CVA approach (Larimore, 1990), the MOESP method (Verhaegen & Dewilde, 1992), the N4SID algorithm (Overschee & Moor, 1994), and the IVM algorithm (Viberg, 1995). It was clarified in the reference (Overschee & Moor, 1995) that some of these SIMs are equivalent to each other besides the use of different weights for data matrix analysis. The asymptotic properties of these SIMs were investigated in the references (Bauer, 2009;Bauer & Jansson, 2000;Chiuso & Picci, 2004). For closed-loop system identification, the plant input is correlated with the output measurement noise due to the feedback mechanism, which may reult in biased estimation if using the above SIMs based on open-loop identification tests. Closed-loop SIMs have therefore been explored to ensure consistent estimation in the recent years. The closed-loop SIMs can be roughly classified into three types: 1. The instrumental variable (IV) based SIMs, e.g. MOESP-variant (Chou & Verhaegen, 1997), N4SIDlike (Overschee & Moor, 1997), CSOPIM (Huang, Ding, & Qin, 2005), and ORT (Katayama, Kawauchi, & Picci, 2005), where the influence of noise was eliminated by using the CONTACT Tao Liu liurouter@ieee.org IVs constructed in terms of the input excitation and output observation data; 2. The predictor-based SIMs, e.g. SSARX (Jansson, 2003), WFA (Chiuso & Picci, 2005), and PBSID-opt (Chiuso, 2007), where a high order ARX model was used to estimate the Markov parameters; 3. The innovation pre-estimation based SIMs, e.g. PARSIM-E (Qin & Ljung, 2003), PARSIM-K (Pannocchia & Calosi, 2010), and OKID-rw (Phan, Horta, Juang, & Longman, 1995), where the innovation sequence was pre-estimated to avoid correlation between the input and noise for Markov parameter estimation. Note that many industrial processes and system operations suffer from deterministic type load disturbance which may occur randomly with varying magnitude or lasting time, such as the injection velocity response of a polymer injection moulding machine during the mould filling process (Liu, Zhou, Yang, & Gao, 2010). However, there are only a few papers that propose bias-eliminated SIMs method against certain types of load disturbance. To deal with a deterministic type disturbance with continuous-time dynamics, a bias-eliminated SIM was proposed (Liu, Huang, & Qin, 2015) which required a fixed occurrence time of load disturbance such that effective approximation along the time sequence could be obtained. Concerning periodic type disturbance, a difference operator was proposed to remove the disturbance effect for consistent estimation (Houtzager, van Wingerden, & Verhaegen, 2013).
In this paper, a bias-eliminated SIM is proposed to cope with load disturbance that may occur randomly with unknown dynamics but finally settle down to a steady value or zero as often encountered in engineering practice. The proposed SIM can be used for both open and closed-loop system identification in the presence of such disturbance. To identify the deterministic system response, the output response is decomposed into three parts, deterministic, disturbed and stochastic components, based on the linear superposition principle. Correspondingly, an LQ decomposition approach is developed to eliminate the influence from load disturbance and stochastic noise, such that unbiased estimation can be obtained for the deterministic system state. Consequently, a shift-invariant approach is given to retrieve the state matrices from the state estimation.
For clarity, the paper is organized as follows. In Section 2, the problem description is briefly stated. The proposed SIM for state estimation is presented in Section 3. In Section 4, consistent estimation is analyzed with a strict proof. Two illustrative examples are shown in Section 5. Finally, some conclusions are drawn in Section 6.
Throughout this paper, the following notations and operations will be used: m×n denotes a m × n real matrix space. 1 m×n ∈ m×n denotes a matrix with all entries of one. The identity (or zero) vector/ matrix with appropriate dimension is denoted by I (or 0), while I m (or 0 m×n ) means I m ∈ m×m (or 0 m×n ∈ m×n ). For any matrix P ∈ m×m of full rank, denote by P −1 the inverse of P, by P T the transpose of P, and by P 2 the matrix 2-norm; for P ∈ m×n of full row (or column) rank, P † denotes the Moore-Penrose pseudo-inverse of P. Denote byP an estimate of P. For Q 1 ∈ p×j , Q 2 ∈ q×j , and Q 3 ∈ r×j with appropriate dimensions, the orthogonal complement of the row space of Q 1 is denoted by Q ⊥ 1 .

Problem description
For an industrial process subject to load disturbance, consider the following state-space model of predictor form, S : where u(t) ∈ n u , x(t) ∈ n x , y(t) ∈ n y denote the input, process state, and output measurement, respectively; (A, B, C) are the state matrices with appropriate dimensions and w(t) ∈ n x indicates an unknown load disturbance which has determinsitic dynamics and settles down to a steady value including zero (due to feedback mechanism in a close-loop system) within a finite time, i.e. w(t) = c for t > t d 0 where t d 0 denotes a finit time length, as often encountered in engineering practice; v(t) ∈ n y denotes the output measurement noise which is assumed to a Gaussian white noise with zero mean and unknown variance, while assuming that w(t) and v(t) are uncorrelated, i.e. E[w(t)v T (t)] = 0. Assume that the system description in (1) is minimal in the sense that (A, B) is reachable and (A, C) is observable. It may be transformed into the Kalman innovation form, wherex(t) denotes the estimated state and K is the Kalman filter gain; The innovation e(t) is a zero-mean white noise which is independent of the input u(k) and output y(k) for k < t. For the convenience of analysis, in the rest of this paper x(t) is used instead ofx(t). DenoteĀ = A − KC, the innovation form in (2) can also be represented by the equivalent predictor form, where all the eigenvalues ofĀ lie inside the unit circle. Without loss of generality, the predictor form in (3) is studied herein for model identification, which can be applied for both an open-loop process and a closed-loop system. We decompose the system state into deterministic and disturbed components in terms of the linear superposition principle, where x u (t) denotes the state response arising from the input u(t), and x d (t) is the disturbance response. Correspondingly, the output response is decomposed into three parts, deterministic, disturbed and stochastic components, i.e.
where y s (t) denotes the output arising from measurment noise v(t), i.e. y s (t) = e(t). Therefore, the system description in (3) can be decomposed into two subsystems, namely deterministic (S u ) and disturbed (S d ) subsystems, as below S d : where ω(t) = B u(t) is assumed to facilitate analysis, with u(t) being regarded as the disturbance input that is unknown and time-varying, but finally recovers to a steady value including zero, i.e. u(t) = c. The identification task is therefore formulated as estimating the system matrices (A, B, C) from the deterministic subsystem (S u ) while eliminating the influence from the disturbed (S d ) subsystem and the innovation from measurement noise (e(t)).

Bias-eliminated subspace identification
Based on the above system description, a two-step SIM method is proposed for bias-eliminated model identification. The first step is estimation of the deterministic state and the second step is identification of the state matrices, which are detailed in the following subsections, respectively.

State estimation
Denote by p and f the past and future horizons. The past and future input stacked vectors are defined, respectively, by t), and e f (t).
Correspondingly, the past and future input block-Hankel matrices are denoted by Similar definitions are given for Y p , Y f , U p , U f , E p and E f . By iterating (3) using the above definitions on the 'past' and 'future' sequences, we obtain where the initial state is regarded as x(k − p). The extended controllability matrix is denoted by The extended observability matrix and the lower triangular Toeplitz matrices are, respectively, When p is sufficiently large, there existsĀ p → 0 as can be verified from (3). Correspondingly, substituting (12) into (13) yields Note that the last n y rows of Y f are in the form of It can be seen from (20) that In view of that [ L 1 ] f , [ L 2 ] f , H f , and G f contain all of the Markov parameters in L, H and G, model identification is therefore focused on (22) so as to avoid redundant estimation.
Based on the subsystem descriptions in (6) and (7), we decompose Y f (l) into three parts, the deterministic output denoted by Y u,f (l) , the disturbed output denoted by Y d,f (l) , and the stochastic output denoted by E f (l) , i.e.
It can be derived from (22) that Since the load disturbance response becomes a steady value after the time where Remark 1: For m = 1, it means that the disturbance is a constant type, which is indeed a special case of a deterministic type disturbance. Since the dynamics of load disturbance is often unknown in practice, it is suggested to take m as a sufficiently large value to ensure that the disturbance response has recovered to a steady value at the time step k − p + m − 1, for which the effectiveness can be verified by comparing the identification results based on taking different values of m.
Using (18) and (19), the estimates of L 1 and L 2 are retrieved fromθ 1 in (29) as followŝ Correspondingly, the product of the extended observability matrix and the deterministic state is estimated using (13) and (20) aŝ where the deterministic state is Performing a singular value decomposition (SVD) on (44), we obtain whereV 1 corresponds to the front n x eigenvalues ofˆ X u . Hence, the deterministic subsystem state is estimated asX Remark 2: Given the input and output observation data, there may be different state-space realizations, all of which are related to each other by a similarity transformation matrix T. Without loss of generality, it is suggested to take T = I for simplicity.

Identification of the state matrices
,J n y = [0 n y ×[(f −1)n y ] I n y ], and J n u = [I n u 0 n u ×[(f −1)n u ] ]. By iterating (5), (6) and (7) using the input and output observation data, we have Similar to (26), the disturbed state response can be written as where θ x d = [x d (t), . . . , x d (t + m − 1)], and Denote two short-hands, It follows from (47) and (48) that The least-squares (LS) estimates ofθ 2 andĈ are therefore obtained asθ Based on the estimates ofÂ,B andK fromθ 2 , the state matrix A can be estimated from the relationship ofĀ = A − KC asÂ

=Â +KĈ (56)
Remark 3: The disturbance effect on the deterministic state and output, θ x d and θ y d , can be similarly estimated as above, which may be used to estimate the dynamics of load disturbance in terms of the decomposed subsystem in (7), therefore facilitating model validation.

Consistent convergence analysis
The estimation error onθ 1 can be computed by (36) and (40) as Recall that the innovation e(t) is a zero-mean white noise which is independent of the input u(k) and output y(k) for k < t. Thus, the innovation process E f (l) is uncorrelated with the past inputs and outputs and future intputs.
Moreover, it follows from (31) that where R 11 and R 22 are nonsingular.
Hence, we have This implies that the orthogonal matrices Q T 1 and Q T 2 are uncorrelated with the future noise E f (l) .
It follows from (57) and (60) that The estimation error ofˆ L 1 andˆ L 2 can be derived from (42) and (43), It follows from (62) and (63) that whereθ 11 consists ofθ 1 with a multiple number of f , i.e. and Similarly, we obtain L 1 = θ 11 H 1 (71) where H 1 and H 2 are totally the same as the H 1 and H 2 in (64) and (65), and θ 11 consists of θ 1 with a multiple number of f , i.e.
Therefore, it can be seen from (73), (74) and (61) that The estimation error onˆ X u can be computed by using (44) as Clearly, it can be seen from (77) and (78) that Take the SVD of X u in the form of The estimation error ofX u can be computed by using (46) with the same T = I as It follows from (45) that By premultiplyingÛ T 1 and postmultiplying V ⊥ 1 to both sides of (83), and usingÛ T 1Û 1 = I and V T 1 V ⊥ 1 = 0, we obtain Also, by premultiplyingÛ T 1 and postmultiplying V ⊥ 1 to both sides of (83), and usingÛ T It follows from (84), (85) and (86) that It follows from (88) that X uy (k))V ⊥ 1 (89) By neglecting the second term of a higher order infinitesimal with respect to X uy (k) 2 in (89), we have The estimation error onX u can be computed by using (82) and (90) as According to (80), we have Using (47) and (48), we obtain The estimation errors onθ 2 andĈ can therefore be derived, respectively, as By omitting the second order infinitesimal with respect to X uy (k) 2 involved with the above computation, there follows In view of that e(t) is uncorrelated with the deterministic system state x u (t), there is Hence, it can be concluded that the estimates onÂ,B,Ĉ andK are surely consistent for N → ∞.

Illustration
Consider an injection moulding process studied in the references (Liu et al., 2015), where (t) denotes an inherent type load disturbance arising from the mould cavity pressure that affects the injection velocity, which was estimated in terms of a transfer function, , where a 1 = 0.15, by assuming the disturbance magnitude to be less than unity while the disturbance occurred at t 0 = 0 (Liu et al., 2015). In fact, the disturbance dynamics is time varying while having variable occurrence time from cycle to cycle. It is therefore assumed that , t 0 and a 1 take random values in these ranges for each test/cycle operation. The measurement noise, v(t), is assumed to be a zeromean random sequence with a variance of 0.2, in view of the fact that the maximal measurement error on the injection velocity is no larger than 1.0 (m/s). The excitation sequence of the valve opening is taken around the operational valve opening of 60% (i.e. u(t 0 ) = 0.6) as a PRBS signal with a standard deviation of σ u = 0.125, while each point of u(t) is executed for 5 sampling times to avoid wearing out the input valve so as to mimic practical implementation of an identification test. A typical scenario of the identification test is shown in Figure 1 under a load disturbance (t) with a 1 = −0.2, t 0 = 400, and δ(t) = 0.01. For model identification, 100 Monte Carlo (MC) tests are performed. By taking the data length of N = 5000, the proposed method using p = f = 10 and m = 1000 gives the estimation results on the eigenvalues of A listed in Table 1, where the result is shown by the mean value along with the standard deviation in parentheses, and the success ratio indicates the number of correctly identified model structures (of stable type) with respect to the total number of identification tests. The identification results obtained by using the N4SID (Overschee & Moor, 1994) and the bias-eliminated SIM (BESIM) (Liu et al., 2015) with p = f = 10 and q = 7 for approximation are also listed in Table 1 for comparison. Scattered plot of the estimated poles are shown in Figure 1.
It is seen that the proposed method gives good estimation accuracy with 100% success ratio, while the recently developed BESIM cannot provide consistent results owing to using the Maclaurin time series for approximation of the load disturbance response, which  requires t 0 = 0 for the occurrence of load disturbance. In contrast, the low success ratio resulted from the well known N4SID may confuse the determination of the true plant model structure in practice.

Conclusions
A bias-eliminated subspace identification method has been proposed to cope with deterministic type load disturbance with unknown but deterministic dynamics, which can be applied to both open-loop and closed-loop system identification. The output response is decomposed into three parts including deterministic, disturbed and stochastic components in the predictor form, based on the linear superposition principle. Correspondingly, an LQ decomposition approach is established to eliminate the influence from load disturbance and stochastic noise for estimating the deterministic system state. Meanwhile, the disturbance dynamics can also be estimated for reference. The state matrices are retrieved by using a shiftinvariant approach. The consistent estimation is clarified with a proof. An illustrative example has well demonstrated the effectiveness and advantage of the proposed method.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work is supported in part by the National Natural Sci-