Subspace identification for industrial processes

Subspace identification for industrial processes

S.D.M. Borjas^I,1; C. Garcia^II,2

^IDCEN, Universidade Federal Rural do Semi-árido, UFERSA, Campus Mossoró, Av. Francisco Mota 570, Bairro Costa e Silva, 59625-900 Mossoró, RN, Brasil

^IIEscola Politécnica da Universidade de São Paulo - EPUSP, Av. Prof. Luciano Gualberto, trav. 3, nr. 158 - Butantã, São Paulo, SP, Brasil

ABSTRACT

Subspace identification has been a topic of research along the last years. Methods as MOESP and N4SID are well known and they use the LQ decomposition of certain matrices of input and output data. Based on these methods, it is introduced the MON4SID method, which uses the techniques MOESP and N4SID.

Palavras-chave: Subspace identification, state space models.

RESUMO

Identificação por subespaços tem sido um tema de pesquisa ao longo dos últimos anos. Métodos como MOESP e N4SID são bem conhecidos eles usam a decomposição LQ de certas matrizes de dados de entrada e saída. Com base nestes métodos, é apresentado o método MON4SID, que usa a técnica dos métodos MOESP e N4SID.

Palavras-chave: Identificação por subespaços, método MOESP, método N4SID, processos industriais.

1. Introduction

Nowadays, an engineer's work consists more and more of obtaining mathematical models of the studied processes [15]. A major part of the literature referring to system identification deals with how to find polynomial models as Prediction Error Methods (PEM) and Instrumental Variable Methods (IVM). In case of complex systems, the state space model appears as an alternative to PEM and IVM models [22]. For multivariable systems, these methods provide reliable state space models directly from input and output data. As systems of large dimensions are usually found in industry, the application of subspace identification algorithms in this field is very promising [3, 8, 10, 14, 16, 17, 24]. Each subspace identification method is different from the other ones in concept, interpretation and computational implementation. An identification generalization is presented in [18]. In [9] a comparative study among the three most commonly used algorithms, i.e., Canonical Variate Analysis (CVA) [13], Multivariable output-Error State sPace (MOESP) [21] and Numerical algorithms for Subspace State Space System IDentification (N4SID) [19] was made. For further details about these algorithms, the reader can consult [2, 6, 7, 12, 17, 21, 22, 23, 4].

In this work it is shown, in a summarized form, the operation of POMOESP [20] and N4SID [19] methods. Combining these methods, it is obtained the MON4SID algorithm, which estimates the extended observability matrix in the same way that POMOESP method; the state sequence is computed through the oblique projection as is done in N4SID method; from this sequence it is obtained the past and future states and finally a consistent estimate of the system matrices is obtained applying the least squares method.

Other form for integrating these methods can be found in [11], in which the author aimed at reducing the computational complexity, dividing MOESP and N4SID algorithms in modules, separating them in independent parts. Afterwards, he used the computationally fastest modules of each algorithm.

1.1. Subspace identification

Linear time invariant systems models, operating in discrete time, are dealt with in the dynamic subspace identification methods. These systems can be described by models in the innovation form

where the symbols represent the input u(k) ∈ ℜ^m, the output y(k) ∈ ℜ^l, the state x(k) ∈ ℜⁿ and the Kalman filter gain K . e(k) ∈ ℜ^l, is zero-mean Gaussian white noise and independent of past input and output data. A , B , C, and D are matrices with appropriate dimensions.

1.2. Subspace identification problem

The subspace identification problem: given u(k) and y(k) a set of input-output measurements aims to, determine the order n of the unknown system, the system matrices ( A , B , C , D ) up to within a similarity transformation and Kalman filter gain K [19].

1.3. Subspace matrix equation

Performing successive iterations in equation (1.1) and (1.2), one can derive the following matrix equation

where subscript f stands for the "future" and p for the "past". The matrices H_i^d and H_i^s are defined as

The past and future input block-Hankel matrices are defined as

where U_p , and U_f ∈ ℜ^mixN. The output and noise innovation block-Hankel matrices Y_p , e Y_f ∈ ℜ^lixN and E_p , and E_f ∈ ℜ^mixN, respectively, are defined in a similar way to (1.4). The states are defined as

The extended observability matrix Γ_i is given by

1.4. Orthogonal projection

The orthogonal projection of the row space of A_x into the row space of B_x is, [19],

where (.)† denotes the Moore-Penrose pseudo-inverse of the matrix (.) . The projection of the row space of A_x into the orthogonal complement of the row space of B_x is, [19],

1.5. Oblique projection

The oblique projection of the row space of G along the row space H into the row space of J is, [19],

Properties of the orthogonal and oblique projections

For a proof, see [19].

2. Identification Methods

2.1. MOESP identification Method

Inside the MOESP family, there is the POMOESP method, which solves the problem in section 1.3, by means of an approximation of the extended observability matrix Γ_i . Therefore, it is necessary to eliminate the last two terms in the right hand side of equation (1.3). That is done in two steps:

First, eliminating the term H_i^d U_f in (1.3), performing an orthogonal projection of equation (1.3) into the row space of U_f^┴ yields

By the property (1.11) it results in U_f/U_f^┴ = 0 . Equation (2.1) can be simplified to

Second, to eliminate the noises in (2.2), it is defined an instrumental variable Z = [ U_p^t Y_p^t ] ^t . Multiplication of (2.2) by Z yields

As it is assumed that the noise is uncorrelated with input and output past data [20], that means that E_f/U_f^┴ Z = 0 .Therefore, (2.3) is written as

In equation (2.4) X_f/U_f^┴ Z = _f is the estimate of the Kalman filter state. Equation (2.4) indicates that the column space of Γ_i can be calculated by the SVD decomposition of Y_f/U_f^┴ Z . For further details, see [19].

2.2. N4SID identification method

N4SID method solves the problem in section 1.3 by means of an approximation of past and future Kalman filter state sequence. From the Theorem 12 in [19],

where Θ_i = Y_f/U_f W_p which is achieved by performing an oblique projection of equation (1.3), along the row space U_f onto the row space of W_p , that is

It is easy to see that the last two terms of equation (2.6) are null, U_f/U_f W_p = 0 by the property of the oblique projection, equation (1.12); E_f/U_f W_p = 0 by the assumption that the noise is uncorrelated with input and output past data [19]. Thus, equation (2.6) can be simplified to

where

_i = X_f/U_f W_p and W_p = [ U_p^t Y_p^t ]^t . Then equation (2.7) is written as

Equation (2.8) indicates that the column space of Γ_i can be calculated by the SVD decomposition of Θ and Γ_i can be calculated as, [19],

Once are known Θ and Γ_i , it is easy to compute _i (calculated from (1.6)). Now _i+1 can be calculated as, [19],

where Θ_i+1 = Y_f^-/U_f^- W_p⁺ and Γ_i-1 denotes the matrix Γ_i without the last l rows. For further details, see [19] or [6].

2.3. MON4SID identification method

In this subsection, it is presented the MON4SID method. To solve the problem in section 1.3, it is used the POMOESP method to calculate the extended observability matrix Γ_i and the N4SID method is employed to calculate the matrices ( A , B , C , D ) through the least squares method.

Γ_i , given in (2.4), can be derived from a simple LQ factorization of a matrix constructed from the block-Hankel matrices U_f , U_p and Y_f , Y_p , in the form

and the orthogonal projection in the left side of (2.4) can be computed by matrix L₃₂ [20]. The SVD of L₃₂ can be given as

The order n of the system is equal to the number of non-zero singular values in S . The column space of U₁ approximates that of Γ_i in a consistent way [20], that is

The system (1.1)-(1.2) can be written as

In equation (2.14), suppose (ideally) that

_i+1 and _i are given, then the system matrices ( A, B, C, D ) could be computed through the least squares method. Therefore, the problem now is to find the state sequences. The oblique projection Θ_i given in equation (2.8) can be computed from (2.11) by

An estimate of the state sequence X is given by

The following matrices are defined:

_i = X (:,1 : N - 1) , _i+1 = X (:, 2 : N) . Thus, the system matrices can be estimated from equation (2.14). The Kalman gain K can be estimated from [19] or [20]

where P is the forward state covariance matrix, which can be determined as solution of the forward Riccati equation

We have then the following algorithm

Algorithm MON4SID

1) Compute the matrices U_f , U_p and Y_f , Y_p and the LQ factorization given in (2.11).

2) Compute the SVD of the matrix L₃₂ from (2.11).

3) Determine the system order by inspection of the singular values in S given in (2.12).

4) Determine Γ_i from equation (2.13) and the state sequence X from (2.16), determine X_i+1 and X_i . 5) Compute the matrices A, B, C, and D from equation (2.14).

3. Simulation

In this section it is employed the Shell benchmark [5], [25] to evaluate the performance of the MON4SID algorithm and to compare it with other existing identification algorithms (PEM, MOESP and N4SID). The Shell benchmark is a model of a two-input two-output distillation column. The inputs are overhead vapor flow (D) and reboiler duty (Q); the outputs are column pressure (P) and product impurity (X). The model of the process is given as follows

e_p(t) is the white noise which generates the disturbance in the pressure. The parameter Ns is used to set the noise level in the simulation. For this study three values are used: 0.2, 0.5 and 1 and they are called 20 % , 50 % and 100 % noise level simulation. The standard deviation of the noise is 1.231. The model of the impurity ( X ) is slightly nonlinear. It is given by

where e_x(t) is white noise with standard deviation 0.677. The Shell benchmark works around the nominal values D_nom = 20 ; Q_nom = 2500; P_nom = 2800; X_nom = 500. It was used a PRBS signal which is persistently exciting of any finite order, obtaining sufficiently informative open-loop experiments. It was collected 1000 samples, 800 samples were used for estimation and the remaining 200 samples were used for validation. The pre-treated signals were of the form F_x = (F - F_nom)/(std(F)) where F can be D, Q, P or X .

4. Model Order Estimation

There is an extensive literature on algorithms to estimate the model order of a linear system in state space. In the subspace algorithms the determination of the model order n is very subtle. Ideally, this information can be determined by the number of singular values different from zero of the orthogonal or the oblique projections of row spaces of the Hankel's block matrix of input and output data. Nevertheless, when the system data contains noise, this number is not easy to calculate. In Figure 1, it is shown the spectrum of SVD values of the process with different values of Ns. Another procedure is to select the value n, which minimizes the estimation errors, a technique used by algorithms such as PEM, which requires a higher computational effort. There is a statistical criterion that can help obtaining the model order: the Akaike Information Criterion (AIC), defined in [1],

where ρ is the number of samples used in the identification, σ²_error(n) is the variance of modeling error for a model of order n with P_n parameters. The index AIC (n) normally reaches a minimum for a certain number of parameters in the model. The application of the AIC criterion to the process with Ns = 0.2 is shown in Figure 2, where one can observe that the minimum value occurs for n = 6 .

5. Model Performance

When there are families of models, one wants to know which of them is the most appropriate for a certain objective, since there is not a model that represents the real system in all its aspects. To evaluate the quality of the model, it is necessary a criterion to measure the distance between the model and the real system. Performance indicators very used are mean relative square error (MRSE) and mean related variance (MVAF), which are defined as

where y is the real output and is the output estimated by the obtained model. The MRSE index is usually used in the literature and the MVAF index is used by the SMI Toolbox. Both performance indexes are used to evaluate the quality of the model obtained by each algorithm, as shown in Table 1 (for Ns = 0.2 and order n = 6 ), Table 2 (for Ns = 0.5 and order n = 6) and Table 3 (for Ns = 1 and order n = 7).

Analyzing the values of Table 1, the MON4SID model is the best one in terms of cross validation. It is verified that the processing time to obtain the model is smaller for MOESP and is much larger for PEM. The MON4SID model was chosen to identify the Shell benchmark process, due to the good balance between its performance and processing time. The PEM model gives better performance with n = 7.

Figure 3 shows the outputs generated by the process and by the identified model for Ns = 0.2. It can be observed that, for a given operating range, the identified model reproduces very well the main dynamic characteristics of the Shell benchmark process. Zero initial conditions were considered.

Analyzing the values in Tables 2 and 3, the MON4SID model is the best one in terms of cross validation. It is verified that the processing time to obtain the model is smaller for MOESP and is much larger for PEM.

6. Conclusions

In this work, it was presented the MON4SID algorithm, which uses LQ factorization in the same way as the MOESP method, which is used to compute the oblique and orthogonal projections; given these projections it is computed the state sequence and the extended observability matrix, respectively. The past and future state sequences are computed from the state sequences which have only one initial state, it does not happen in the N4SID method, since for each oblique projection ( Θ_i and Θ_i+1 ) different state sequences ( _i+1 and _i ) are computed, generating a problem of bias in the estimates. To evaluate this algorithm, three identification algorithms (MOESP, N4SID, PEM) were applied to the Shell benchmark process, to identify a MIMO model in discrete time state space and their results were compared. It is possible to observe that a linear model can provide a good description of a non-linear system within a certain operating range. The Akaike criterion provided the model order. The performance comparison was made according to cross validation for each algorithm, employing two different performance criteria. For this specific case, the MON4SID model was chosen, due to its performance and processing time. PEM was the slowest model out of the four tested ones. The obtained model is observable, controllable and asymptotically stable within a certain range of operation, and it can be applied in control and monitoring applications.

Recebido em 16 Setembro 2010

Aceito em 25 Maio 2011.

• [1] H. Akaike, Information theory and an extension of the maximum likelihood principle. In: Second International Symposium on Information Theory, Budapest, Hungary. (B.N. Petrov, F. Csaki, eds.), pp. 267-281, 1973.
• [2] T. Backx, "Identification of an Industrial Process: A Markov Parameter Approach", Ph.D. Thesis, Technical Univ. Eindhoven, The Netherlands, 1987.
• [3] S.D.M. Borjas, C. Garcia, Subspace identification using the integration of MOESP and N4SID methods applied to the Shell benchmark of a distillation column, artigo aceito 9th Brazilian Conference on Dynamics Control and their Applications DINCON 2010, Serra Negra/SP Brasil, pp. 57, 2010.
• [4] S.D.M. Borjas, C. Garcia, Modelagem de FCC usando métodos de identificação por predição de erro e por subespaços, IEEE América Latina, Revista virtual - na Internet, 2, No. 2 (2004), 108-113.
• [5] B. Cott, Introduction to the Process Identification, Workshop at the 1992 Canadian Chemical Engineering Conference, Journal of Process Control, 5, No. 2 (1995), 67-69.
• [6] K. De Cock, B. De Moor, Subspace identification methods, in Contribution to section 5.5, Control systems robotics and automation of EOLSS, UNESCO Encyclopedia of life support systems, (Unbehauen H.D.), 1 of 3, Eolss Publishers Co., Ltd., Oxford, UK, pp. 933-979, 2003.
• [7] B. De Moor, P. Van Overschee, W. Favoreel, Algorithms for subspace state space system identification - an overview, In Applied and computational control, signal and circuits, (B. Datta Ed.), Vol. 1, pp. 247-311. Birkhauser: Boston (Chapter 6), 1999.
• [8] W. Favoreel, B. De Moor, P. Van Overschee, Subspace state space system identification for industrial processes, Journal of Process Control, 10, No. 2-3 (2000), 149-155.
• [9] W. Favoreel, S. Van Huffel, B. De Moor, V. Sima, M. Verhaegen, Comparative study between three subspace identification algorithms, Niconet, 1998.
• [10] B. Haverkamp, Efficient implementation of subspace method identification algorithms, Niconet, 1999.
• [11] B. Haverkamp, M. Verhaegen, "SMI Toolbox: state space model identification software for multivariable dynamical systems", Vol. 1, Delft University of Technology, The Netherlands, 1997.
• [12] T. Katayama, "Subspace Methods for System Identification", Springer, London, 2005.
• [13] W. Larimore, Canonical variate analysis in identification, filtering and adaptive control, In "Proc. 29th Conference on Decision and Control", Hawai, USA, pp. 596-604, 1990.
• [14] W. Larimore, Automated multivariable system identification and industrial applications, In: "American Control Conference, ACC'99", San Diego, CA, Proceedings, Vol. 2, pp. 1148-1162, 1999.
• [15] L. Ljung, "System Identification Theory for the User", Prentice Hall Englewood Cliffs, NJ, 1999.
• [16] G. Mercere, L. Bako S. Lecouche, Propagator-based methods for recursive subspace model identification, Signal Processing, 88, No. 3 (2008), 468-491.
• [17] P. Roberto, G. Kurka, H. Cambraia, Application of a multivariable input-output subspace identification technique in structural analysis, Journal of Sound and Vibration, 312, No. 3 (2008), 461-47.
• [18] P. Van Overschee B. De Moor, A unifying theorem for three subspace system identification algorithms, Automatica, (Special Issue on Trends in System Identification), 31, No. 12, (1995), 1853-1864.
• [19] P. Van Overschee B. De Moor, "Subspace Identification for Linear Systems: Theory, Implementation, Applications", Dordrecht: Kluwer Academic Publishers, 1996.
• [20] M. Verhaegen, P. Dewilde, Subspace model identification. part i: the output-error state-space model identification class of algorithms, International Journal of Control, 56, No. 1 (1992) 1187-1210.
• [21] M. Verhaegen, Identification of the deterministic part of MIMO state space models given in innovation form from input-output data, Automatica (Special issue on Statistical Signal Processing and Control) , 30, No. 1, (1994) 61-74.
• [22] M. Viberg, Subspace methods in systems identification. In: 10th IFAC Symposium on System Identification, SYSID'94, Copenhagen, Denmark, Proceedings, Vol. 1, pp. 1-12,1994.
• [23] M. Viberg, Subspace-based methods for the identification of linear time-invariant system. Automatica, (Special Issue on Trends in System Identification), 31, No. 12 (1995), 1835-1851.
• [24] M. Viberg, Subspace-based state-space system identification, Circuits, Systems and Signal Processing, 21, No. 1 (2002), 23-37.
• [25] Y. Zhu, Multivariable process identification for MPC: the asymptotic method and its applications, Journal of Process Control, 8, No. 2, (1998) 101-115.

Publication Dates

• Publication in this collection
  05 Oct 2012
• Date of issue
  Dec 2011

History

• Received
  16 Sept 2010
• Accepted
  25 May 2011