Subspace-based system identification: weighting and pre-filtering of instruments

doi:10.1016/S0005-1098(01)00235-7

Automatica

Volume 38, Issue 3, March 2002, Pages 433-443

https://doi.org/10.1016/S0005-1098(01)00235-7 Get rights and content

Abstract

Subspace-based system identification is typically based on an estimate of the extended observability matrix. It is thus of great interest to investigate, and also optimize, the estimate of the observability matrix. Of special interest in this paper is the fact that the influence of certain weighting matrices is an unresolved issue in the literature on subspace identification. Here, an asymptotic analysis of the estimated observability matrix is presented. The main result of the analysis is that novel weighting matrices and pre-filters of instrumental variables are derived.

Introduction

Methods for Subspace-based State Space System Identification (4SID) have lately been suggested as alternatives to traditional system identification techniques, see e.g. Larimore (1983), Ljung (1999), VanOverschee and De Moor (1996), Verhaegen (1994) and Viberg (1995). The 4SID methods are attractive since they estimate a state-space realization directly from input–output data, without requiring canonical parameterizations and non-linear optimizations, which typically is the case for the Prediction Error Method (PEM), cf. Ljung (1999), Söderström and Stoica (1989). Subspace methods are often implemented using robust numerical tools such as the QR-factorization and the Singular Value Decomposition (SVD), which make them attractive from a numerical point of view.

Although 4SID methods have been studied for a while and several successful applications, also on real-world data, have been reported, there is still work to be done in terms of analysis and performance optimization. In the literature, several contributions in this direction have appeared, see e.g. Bauer and Jansson (2000), Jansson and Wahlberg (1996) and Viberg, Wahlberg, and Ottersten (1997). In Bauer and Jansson (2000) it was shown that the estimated transfer function, under a set of assumptions, asymptotically has a normal distribution assuming application of a MOESP-like (Multivariable Output-Error state SPace identification, cf. Verhaegen, 1994) 4SID algorithm. The expression for the variance of the estimated transfer function is however complicated and was not explicitly provided. In Viberg et al. (1997), an optimally weighted subspace fitting (WSF) procedure was proposed. The WSF approach is based on a parameterization of the left nullspace of the observability matrix. The optimality is in terms of asymptotically producing minimum-variance estimates of the parameters describing the left nullspace of the observability matrix, given an estimate of the observability matrix (i.e. the subspace estimate). In Viberg et al. (1997) it was not guaranteed that the applied subspace estimate was the best possible. Given this rather arbitrary subspace estimate, the parameters describing the left null-space of the observability matrix were however estimated in an optimal manner. An open issue here is then how to optimize the subspace estimate. In Viberg et al. (1997) also the asymptotic variance of the estimated poles was derived, assuming a “shift-invariant approach” for estimation of system poles. In Jansson and Wahlberg (1996), 4SID algorithms were derived using a reduced-rank weighted least squares approach, which produced natural choices of weighting matrices. Similar results were obtained in Katayama and Picci (1999), where a multi-stage least squares approach was studied. Although different approaches of analysis, it is interesting to note that the column weightings found in Jansson and Wahlberg (1996), Katayama and Picci (1999) and Verhaegen (1994) all correspond to the CVA weighting of Larimore (1983).

The available statistical results on 4SID methods are in conclusion rather complicated, and the influence of certain user-defined quantities is not fully understood. The objective of the work in the present paper is to carry out a statistical investigation of 4SID methods. We will only address parts of this difficult problem. The first step in 4SID methods is in general to estimate the observability matrix. The particular state-space realization that the estimated observability matrix belongs to is in general unknown, and only its range space is retrieved. The basic idea of the present paper is then to analyze the asymptotic variance of the estimated observability matrix, rather than to consider the more difficult problem of analyzing the variance of the estimated transfer function. From this analysis, the following issues are addressed:

The analysis of the present paper is based on a sub-optimal figure of merit that not necessarily is related to the variance of the estimated transfer function. To appreciate the obtained results on weightings and pre-filters, it is of course desirable to relate our findings to relevant system quantities. For this purpose we analyze available results on the asymptotic variance of the estimated poles, and we also relate these results to our findings. Several numerical examples illustrating the benefits of the new weighting and pre-filtering strategies are included in the paper.

Section snippets

Notations and general assumptions

Assume that the discrete-time linear time-invariant dynamical system under study can be described by the following n^th order state-space model with p outputs and m inputs collected in $y (t)$ and $u (t)$ , respectively: $x (t+1)= Ax (t)+ Bu (t)+ w (t), y (t)= Cx (t)+ Du (t)+ v (t).$ Here, $w (t)$ denotes the n-dimensional process noise, $v (t)$ denotes the p-dimensional measurement noise, and $x (t)$ denotes the n-dimensional state vector. The system matrices are of dimensions $A ∈ R^{n×n}, B ∈ R^{n×m}, C ∈ R^{p×n}$ , and $D ∈ R^{p×m}$ respectively. The

Analysis of $Γ ̂$

The subspace estimate $Γ ̂$ consists in general of components from both $R^{⊥} {Γ}$ (the orthogonal complement of $R {Γ}$ ) and $R {Γ}$ . On the other hand, only the part of $Γ ̂$ that lies in $R^{⊥} {Γ}$ actually is an error. With this in mind, we decompose $Γ ̂$ as $Γ ̂ = Π_{Γ} Γ ̂ + Π_{Γ}^{⊥} Γ ̂,$ where $Π_{Γ}^{⊥}$ denotes the orthogonal projection onto $R^{⊥} {Γ}$ , and $Π_{Γ}$ denotes the orthogonal projection onto $R {Γ}$ , i.e. $Π_{Γ} = Γ (Γ^{T} Γ)^{−1} Γ^{T}, Π_{Γ}^{⊥} = I − Π_{Γ} .$ Instead of analyzing the properties of the “subspace error” $Π_{Γ}^{⊥} Γ ̂$ , we have chosen to study the following

Results on weighting matrices

The purpose of this section is to derive weighting matrices, and to outline an algorithm for subspace identification. The starting point of the analysis lies in the asymptotic covariance matrix of ε_N.

Analysis of the variance of the pole estimate

The purpose of this section is to relate our findings to existing statistical results. For that purpose, recall the common “shift-invariant” approach for estimating $A$ , see e.g. Viberg et al. (1997). Define the selection matrices $J_{1} =[I_{(α−1)p} 0_{(α−1)p×p},],$ $J_{2} =[0_{(α−1)p×p} I_{(α−1)p}],$ where $I_{K}$ is the K×K identity matrix, and $0_{K×L}$ denotes an K×L matrix with zeros. The shift-invariant structure of $Γ$ implies that $J_{1} Γ A = J_{2} Γ .$ Given $Γ ̂$ , an estimate of $A$ can be found by solving (57) in the least-squares sense. In

Pre-filtering of $p (t)$

In the identification literature it has been noted that IV methods in general perform worse than the PEM. In order to make IV methods to perform as well as the PEM, a pre-filtering of the instruments is typically required, cf. Söderström and Stoica (1989, Chapter 8). Thus, considering IV-interpretations of 4SID methods, it is natural to investigate whether pre-filtering of the IV-vector $p (t)$ can improve the subspace estimate. For simplicity, only scalar “pre-filters” are considered.

The W4SID algorithm

In this section, previously obtained results are illustrated by means of a couple of numerical examples. In our first example, we consider identification of the following simple ARMAX SISO system: $(1−q^{−1} +0.5q^{−2})y(t)= (1−0.5q^{−2})u(t)+(1+q^{−1} +0.5q^{−2})e(t),$ where e(t) is white Gaussian noise. Three different identification algorithms are investigated:

Alg1:	The PO-MOESP weighting (i.e. $W_{c} = W ̂_{c}^{oa}$ ).
Alg2:	The W4SID algorithm.
Alg3:	PEM, implemented using the ARMAX routine in Matlab's System Identification

Conclusions

We have presented an analysis on the asymptotic error of the estimated observability matrix. Once a particular figure of merit was defined, useful results on how to choose weighting matrices were found. A particularly important result is that we found a novel weighting matrix, which in our numerical examples showed superior accuracy compared to PO-MOESP. From the numerical examples, we draw the following conclusions:

•
The W4SID algorithm performs favorable to MOESP.
•
Pre-filtering of the

Tony Gustafsson was born in Värnamo, Sweden, in 1969. He received the M.S. degree in electrical engineering from Chalmers University of Technology, Sweden, in 1994, and in 1999 he received the Ph.D. degree in signal processing from the same university. From 1999 to 2000, he was a postdoctoral researcher at University of California San Diego. Presently he is a research engineer at Switchcore corporation, Göteborg, Sweden.

References (14)

D. Bauer et al.
Analysis of the asymptotic properties of the MOESP type of subspace algorithms
Automatica
(2000)
M. Jansson et al.
A linear regression approach to state-space subspace system identification
Signal Processing
(1996)
M. Jansson et al.
On consistency of subspace methods for system identification
Automatica
(1998)
T. Katayama et al.
Realization of stochastic systems with exogenous inputs and subspace identification methods
Automatica
(1999)
M. Verhaegen
Identification of the deterministic part of MIMO state space models given in innovations form from input–output data
Automatica
(1994)
M. Viberg et al.
Analysis of state space system identification methods based on instrumental variables and subspace fitting
Automatica
(1997)
A. Albert
Regression and the Moore-Penrose pseudoinverse
(1972)

There are more references available in the full text version of this article.

Cited by (30)

Koopman operator dynamical models: Learning, analysis and control
2021, Annual Reviews in Control
The Koopman operator allows for handling nonlinear systems through a globally linear representation. In general, the operator is infinite-dimensional – necessitating finite approximations – for which there is no overarching framework. Although there are principled ways of learning such finite approximations, they are in many instances overlooked in favor of, often ill-posed and unstructured methods. Also, Koopman operator theory has long-standing connections to known system-theoretic and dynamical system notions that are not universally recognized. Given the former and latter realities, this work aims to bridge the gap between various concepts regarding both theory and tractable realizations. Firstly, we review data-driven representations (both unstructured and structured) for Koopman operator dynamical models, categorizing various existing methodologies and highlighting their differences. Furthermore, we provide concise insight into the paradigm’s relation to system-theoretic notions and analyze the prospect of using the paradigm for modeling control systems. Additionally, we outline the current challenges and comment on future perspectives.
Variance reduction in Covariance Based Realization Algorithm with application to closed-loop data
2020, Automatica
Citation Excerpt :
In addition, for MIMO state-space models, an optimal column weighting matrix was proposed for MOESP-like (Multivariable Output-Error state SPace identification) subspace methods in Gustafsson (2002) by using a similar idea from Stoica and Söderström (1983) and Viberg, Wahlberg, and Ottersten (1997). An optimal column weighting for CoBRA in an open-loop setting was then proposed theoretically in Hu and de Callafon (2017) based on the work in Gustafsson (2002) and Viberg et al. (1997). This paper analyzes a two stage CoBRA for closed-loop MIMO system identification, where the work from Hu and de Callafon (2017) is extended to the first stage.
The Covariance Based Realization Algorithm (CoBRA), one branch of subspace methods, enables the estimation of multivariable models using a large number of data points due to the use of finite size covariance matrices. In addition, the covariance pre-processing allows CoBRA to ignore any (high-order) noise dynamics and focus on the estimation of the low-order deterministic model. However, subspace methods and CoBRA in particular are not maximum-likelihood methods. In this paper, an in-depth study on the statistical behavior of the noise effects is conducted. An approach is provided to reduce the variance of an estimate obtained by CoBRA via the choice of optimal row and column weighting matrices. For closed-loop implementation of CoBRA, a two stage procedure is proposed with the estimate on an intermediate instrument. At the first stage, a least-length perturbation of a scalar accuracy function is used to obtain an analytic solution for the optimal weighting matrices. The resulting instruments produced by the first stage are used for the identification at the second stage to extract the plant model. Simulation examples are given to illustrate the efficiency of the proposed two stage CoBRA method on the basis of closed-loop data and compared with other subspace methods.
Filtered-covariance function-based subspace identification with bound effects integration
2014, IFAC Proceedings Volumes (IFAC-PapersOnline)
A novel subspace covariance-based identification method that achieves consistent system estimates in the presence of highly-colored noise without requiring the use of weighting matrices has been derived recently. However, the use of Fourier transform in the computation of correlation functions requires the need to include the estimation of initial and final conditions. In this article we propose a solution to get unbiased estimates regardless of initial or end effects. Moreover, to further improve the model order estimation a filter-based solution is considered. This results in a new covariance-based algorithm that uses frequency weights computed with limit conditions in mind. A simulation example illustrates the algorithm performance.
An extended closed-loop subspace identification method for error-in-variables systems
2012, Chinese Journal of Chemical Engineering
A closed-loop subspace identification method is proposed for industrial systems subject to noisy input-output observations, known as the error-in-variables (EIV) problem. Using the orthogonal projection approach to eliminate the noise influence, consistent estimation is guaranteed for the deterministic part of such a system. A strict proof is given for analyzing the rank condition for such orthogonal projection, in order to use the principal component analysis (PCA) based singular value decomposition (SVD) to derive the extended observability matrix and lower triangular Toeliptz matrix of the plant state-space model. In the result, the plant state matrices can be retrieved in a transparent manner from the above matrices. An illustrative example is shown to demonstrate the effectiveness and merits of the proposed subspace identification method.
A bias-eliminated subspace identification method for errors-in-variables systems
2012, IFAC Proceedings Volumes (IFAC-PapersOnline)
For model identification of industrial operating systems subject to noisy input-output observations, known as the error-in-variables (EIV) problem, a subspace identification method is proposed in this paper by developing an orthogonal projection approach to guarantee consistent estimation of the deterministic part of such a system. The rank condition for such orthogonal projection is analyzed in terms of the nominal state-space model structure. Using the principal component analysis (PCA), the extended observability matrix and low triangular block-Toeliptz matrix of the state-space model are analytically derived. Accordingly, the system state-space matrices can be retrieved in a transparent manner from the above matrices through linear algebra or an ordinary least-squares (LS) algorithm. A benchmark example used in the existing references is adopted to demonstrate the effectiveness and merit of the proposed subspace identification method.
A predictor form PARSIMonious algorithm for closed-loop subspace identification
2010, Journal of Process Control
In this paper, we present a novel subspace identification algorithm in which all non-causal terms are removed, and the specific Toeplitz structure of the Markov parameter’s matrices is fully exploited in the spirit of the so-called PARSIMonious algorithms. We use the state-space formulation in predictor form, and we show that consistent estimates of the Markov parameters are granted both for open-loop and closed-loop data. Furthermore, we propose to evaluate the system matrices $(B_{K} = B - KD, D, K)$ and the initial condition by a single Least Squares problem, which is well conditioned even for unstable systems. We present identification results for two multi-variable systems to show the main features of the proposed method, and to assess its performance against that achieved by other subspace methods. Furthermore, we compare the performance of Model Predictive Controllers, based on models identified from closed-loop data using the different subspace algorithms, in the control of the Wood–Berry distillation column. Results indicate that the proposed method is suitable for MPC design purposes, and compares favorably with the other subspace algorithms.

View all citing articles on Scopus

^☆: This paper was originally presented at the 12th IFAC Symposium on System Identification, 2000, Santa Barbara, USA. This paper was recommended for publication in revised form by Associate Editor Brett Ninness under the direction of Editor Torsten Söderström. This work was completed while visiting University of California San Diego. Support by the Swedish Foundation for International Cooperation in Research and Higher Education, and Telefonaktiebolaget LM Ericsson is gratefully acknowledged.

View full text

Subspace-based system identification: weighting and pre-filtering of instruments☆

Abstract

Introduction

Section snippets

Notations and general assumptions

Analysis of Γ̂

Results on weighting matrices

Analysis of the variance of the pole estimate

Pre-filtering of p(t)

The W4SID algorithm

Conclusions

Automatica

Signal Processing

Automatica

Automatica

Automatica

Automatica

Regression and the Moore-Penrose pseudoinverse

Analysis of $Γ ̂$

Pre-filtering of $p (t)$