Task III.A Development of standard software for linear, time-invariant state space model identification

A toolbox for linear systems identification has been developed and included in the SLICOT package. The development of this toolbox was the subject of Task III.A of the European Community BRITE-EURAM III Thematic Networks Programme NICONET (contract number BRRT-CT97-5040). The main functionalities of the toolbox include: identification of linear, time-invariant discrete-time multivariable state space systems (A, B, C, D);  identification of state and output (cross-)covariance matrices for such systems;  estimation of the associated Kalman gain matrix;  estimation of the initial state of discrete-time linear systems. In the development of the toolbox, computational reliability  and numerical efficiency  were our main concerns. The coding of all subroutines uses extensively the linear algebra standard package LAPACK and is done according to the SLICOT implementation and documentation standards.  All developed subroutines  have been extensively  tested  on various test examples and are fully documented.  Their performance has also been evaluated in several industrial applications:  a glass furnace, an industrial evaporator, a CD player arm, etc. For MATLAB and Scilab users' convenience, mex- and m-functions have been developed so that the major functionality of SLICOT routines is available within these user-friendly environments. The linear system identification toolbox is described in detail in the SLICOT Working Note report 2000-4, together with  the above-mentioned, and other industrial applications.

Description of the Linear System Identification Toolbox

The toolbox subroutines

The Fortran routines for TASK III.A  are aimed to estimate state space models and covariance matrices of linear discrete-time multivariable systems, using the available input-output data sequences. Reliability, efficiency, and ability to solve large, industrial identification problems received a special consideration in the toolbox development process. Flexibility of usage is enhanced by providing various options: two algorithmic subspace-based approaches (MOESP - Multivariable Output-Error state SPace identification, and N4SID - Numerical algorithms for Subspace State Space System IDentification) and their combination, standard or fast techniques for data compression (including abilities to exploit the block-Hankel structure), multiple (possibly connected) data batches processing, availability of both, fully documented drivers and computational routines, the use of structure exploiting algorithms and dedicated linear algebra tools, etc. In addition, optionally, the quality of the intermediate results can be assessed by inspecting the associated condition numbers. A list of the implemented Fortran routines with links to the associated .html documentation is given in the following table:
 
IB01AD Input-output data preprocessing and finding the system order 
IB01BD Estimating the system matrices, covariances, and Kalman gain
IB01CD Estimating the initial state and the system matrices B and D
IB01MD *Upper triangular factor in the QR factorization of the concatenated input and output block-Hankel matrices
IB01ND *Singular value decomposition giving the system order
IB01OD *Estimating the system order
IB01OY *User's confirmation of the system order
IB01PD *Estimating the system matrices and covariances
IB01PX *Estimating the matrices B and D of a system using Kronecker products
IB01PY *Estimating the matrices B and D of a system exploiting the special structure
IB01QD *Estimating the initial state and the matrices B and D of a system, using A, C, and the input-output data sequences
IB01RD *Estimating the initial state of a system, using (A, B, C, D), and the input-output data sequences

The documentation of routines are also accessible from the SLICOT Library main index (in case of the drivers, or user-callable routines), or from the SLICOT Supporting Routines index (in case of the auxiliary routines, marked with * in the list). The SLICOT Supporting Routines index is also accessible from the main Library index.

The SLICOT library is available from the SLICOT homepage www.slicot.de. All routines in the above list are included in the files slicot.tar.gz (Unix version) and slicotPC.zip (Windows version), downloadable from the SLICOT homepage.
 

Mex- and m-functions

The mex- and m-functions implemented within TASK III.A  are aimed to provide convenient, efficient and reliable solutions to the system identification problems,  for linear discrete-time multivariable systems, in the user-friendly MATLAB environment. They are built on the SLICOT Library routines implemented within TASK III.A. Executable SLICOT mexfiles are provided for MATLAB running under WINDOWS 95/98/NT/ME/2000. Note that the Fortran source texts for the mexfiles have been written using the Fortran 90 dynamic allocation feature.

The following table contains the list of mexfiles with their function:

order Input-output data preprocessing, possibly sequentially, and finding an estimate of the system order
sident System matrices, Kalman predictor gain, and covariance matrices estimation, using MOESP, N4SID, or their combination
findBD Estimating the initial state and/or the matrices B and D, using A, C, and the input and output trajectories

These three mexfiles provide interfaces to the main user-callable identification routines IB01AD, IB01BD, and IB01CD, respectively. All functionality available in the SLICOT identification routines has been implemented. For users' convenience, several easy-to-use interface m-functions have been additionally included in the identification toolbox, explicitly addressing some of supported features. Moreover, whenever possible, these m-functions allow to work with system objects defined in the Control Toolbox of MATLAB. The following table contains the list of implemented m-functions:
 

findR Input-output data preprocessing, using Cholesky or (fast) QR factorization and MOESP or N4SID identification techniques, and estimating the system order
findABCD System matrices and Kalman gain estimation, using MOESP, N4SID, or their combination
findAC Estimating the matrices A and C, using MOESP or N4SID
findBDK Estimating the matrices B, D, and Kalman gain K (given A and C), using MOESP or N4SID
findx0BD Estimating the initial state and/or the matrices B and D, given the matrices A, C, and a set of input-output data
inistate Estimating the initial state, given the system matrices, and a set of input-output data
slmoesp System matrices and the Kalman gain estimation, using MOESP technique
sln4sid System matrices and the Kalman gain estimation, using N4SID technique
slmoen4 System matrices and the Kalman gain estimation, using combined MOESP and N4SID techniques: A and C found via MOESP, and B and D, via N4SID
slmoesm System matrices, the Kalman gain, and initial state estimation, using combined MOESP and system simulation techniques

The last four m-functions are method-oriented, while the first five allow to flexibly identify various system and covariance matrices. The method-oriented m-functions also enable to efficiently estimate models of various orders. The executable mexfiles for MATLAB running under WINDOWS 95/98/NT/ME/2000, together with the source code, the above m-files and several MATLAB test scripts are available as ident_mex.zip from the SLICOT homepage

Similar interfaces have been realized to integrate the SLICOT system identification software into Scilab. The implemented Scilab functions are mostly compatible with the MATLAB interface and, in particular, the names of the mex-files and m-files are the same.

A demonstration package for a PC machine is available as slident_demo.zip from the SLCIOT homepage. After decompressing that file in the current directory, the demo can be directly tried by typing slidemo in a MATLAB 6.1 session (from that directory).

Industrial case studies

Several industrial examples have been considered for testing the SLICOT system identification software: a glass furnace, an industrial evaporator, a CD player arm, etc. These, as well as other industrial examples are described in the following NICONET Report:
A. van den Boom, T. Backx, Y. Zhu. Benchmarks for Identification, NICONET Report 1999-19, July 2000.
Several industrial examples are used as applications in another NICONET Report,

V. Sima. SLICOT Linear Systems Identification Toolbox, SLICOT Working Note 2000-4, July 2000.

Glass furnace model. This data set from the DAISY collection contains data for a glass furnace (Philips); there are 3 inputs, 6 outputs, and 1247 samples. The number of block rows s was chosen 10 and the model order was set to 5. The estimated output follows well the original output. Click here to see the good agreement obtained for the fifth output and its estimate. Standard MOESP approach (the slmoesp m-function) has been used to obtain these results.

Industrial evaporator. The input-output data set for a four-stage evaporator is also included in the DAISY collection; there are 3 inputs (feed flow to the first evaporator stage, vapor flow to the first evaporator stage, and cooling water flow), 3 outputs (dry matter content, flow of the outcoming product, and temperature of the outcoming product), and 6305 samples. The number of block rows s was chosen 10 and the model order was set to 4. The temperature of the outcoming product, and its filtered estimate, obtained using the identified system matrices and Kalman gain, are displayed here. The filtered output follows well the original output. The m-function slmoesm has been used.

CD-player arm. This data set, also from the DAISY collection, contains data for the mechanical construction of a CD-player arm; the system has 2 inputs (the forces of the mechanical actuator), and 2 outputs (tracking accuracy of the arm); 2048 input-output samples are available. The inputs are highly colored. The number of block rows s was chosen 15 and the model order was set to 8. Click here to see the good agreement obtained for the second output and its filtered estimate, obtained using the identified system matrices and Kalman gain. The m-function slmoesm has been used.

Test examples

Extensive testing of the implemented system identification software has been performed using several benchmark problems. Some examples used for tests are described in Section 5 of the report V. Sima. "SLICOT Linear Systems Identification Toolbox", SLICOT Working Note 2000-4, July 2000. Included is here, as an example, a mechanical flutter application (1 input, 1 output); 1024 input-output samples are available. The number of block rows s was chosen 20 and the model order was set to 6. The input trajectory is plotted here. This application is particularly difficult for the standard MOESP approach, because ill-conditioned matrices appear during the computation of B and D. However, there is no difficulty in solving the identification problem using the other methods, implemented in the m-files sln4sid, slmoen4, or slmoesm. Click here to see the good agreement between the output, its estimate, and the filtered output estimate, obtained using the identified system matrices, without and with the Kalman gain, respectively. The m-function slmoesm has been used. The filtered output is practically indistinguishable from the output.

Accuracy and efficiency comparisons of the SLICOT identification software and the available subspace-based techniques for 9 applications, including those discussed above, are presented in the report V. Sima. "SLICOT Linear Systems Identification Toolbox", SLICOT Working Note 2000-4, July 2000. Click here to see such a timing comparison when the SLICOT routines used the Cholesky factorization algorithm for the input-output data compression. This algorithm worked for all, but the application #4. In the figure, SLmoesp and SLn4sid denote the SLICOT implementations of the standard MOESP and N4SID approaches, OMOESP denotes a previous MATLAB code using Total Least Squares for finding B and D, and MOESP denotes a new MATLAB code which computes the B and D matrices by solving a large least squares problem built using the matrices A and C, and the input and output trajectories of the system. When successful, the Cholesky factorization algorithm significantly increased the efficiency, while preserving the same accuracy as for the QR factorization algorithm. The speed-up factors vary between 10 and 20 when comparing to the SLICOT QR factorization algorithm, and between 15 and 40 when comparing to MATLAB codes.

Available Reports

1.  W. Favoreel, V. Sima, S. Van Huffel, M. Verhaegen, B. De Moor "Subspace Model Identification of Linear Systems in SLICOT", SLICOT Working Note 1998-6, Sept. 1998.
2.  B. Haverkamp "Efficient Implementation of Subspace Method Identification", Niconet Report 1999-3, Jan. 1999.
3. A. van den Boom, T. Backx, Y. Zhu "Benchmarks for Identification", NICONET Report 1999-19, July 2000.
4.  V. Sima "SLICOT Linear Systems Identification Toolbox", SLICOT Working Note 2000-4, July 2000.

Vasile Sima,  July 28, 2000
Updated: February 22, 2001;  April 18, 2002;  January 18, 2004