Sima, V. (1)

Varga, A.

Verhaegen, M. (2) (1) Inst. for Informatics Bucharest

English

Institute of Transport Research:Publications

LAPACK-Based Software for Subspace IdentificationVasile Sima∗ Andras Varga† Michel Verhaegen‡June 28.19961 IntroductionComputer-aided control system design is usually based on linear time-invariant (LTI) statespace models. Surprisingly enough, multivariable state space system identification has onlyrecently become a topic of intense research. In contrast with classical input-output (I/O)identification approaches, the newly proposed Subspace Model Identification (SMI) techniquesessentially find a state sequence, or a column space approximation, and then determine thesystem matrices by solving some least squares problems. These techniques have promising ad-vantages over the classical ones. One advantage is that there is no need for parameterizations,which are notoriously difficult to choose or analyse, or could led to ill-conditioned problems.Another advantage is that robust numerical linear algebra techniques, like QR factorizationand singular value decomposition (SVD), can be largely used; this contrasts with the iterativeoptimization schemes required in the parametric model identification approach, documentede.g. in [1]. The attractiveness of SMI techniques is further increased by the small number ofparameters (essentially only one) to be selected for determining the model structure, withoutany restriction on model generality. See, for instance, [2]–[7] for a further discussion of theSMI features.MATLAB codes based on SMI algorithms have been recently developed, e.g. in [8]. Forefficiency, it is clearly useful to implement some algorithms in Fortran, using the state-of-the-art, public-domain linear algebra package LAPACK. This allows to exploit the potentialparallelism of many modern computer architectures.This paper briefly describes the basic approaches and LAPACK-based Fortran routines formultivariable system identification by subspace techniques. The developed software pack-age [9] includes implementations of two classes of SMI techniques to identify an LTI system.The first one is referred to as the State Intersection (SI) class of SMI techniques. The imple-mented method, described in [2], is known as the N4SID (Numerical algorithm for SubspaceState Space System IDentification) approach. The second one is referred to as the Multi-variable Output Error state SPace (MOESP) class of techniques. The implemented MOESPmethods are described in [3], [5] and [6]. Deterministic and combined deterministic-stochastic∗Research Institute for Informatics, Bd. Mares¸al Averescu, Nr. 8–10, 71316 Bucharest 1, Romania. e-mail:vsima@roearn.ici.ro†DLR-Oberpfaffenhofen, German Aerospace Research Establishment, Institute for Robotics and SystemDynamics, D-82234 Weßling, FRG. e-mail: Andreas.Varga@dlr.de.‡Delft University of Technology, Dept. of Electrical Engineering, P.O.Box 5031, NL-2600 GA Delft, TheNetherlands. e-mail: verhaege@dutentb.et.tudelft.nl1identification problems are dealt with. A state space model is computed from I/O data se-quences. It is possible to handle multiple I/O sequences, each one being collected by a possiblyindependent identification experiment. Sequential processing of large data sets is providedas an option. See [9], for further details. In addition to the main identification facilities,auxiliary routines are available to detect dead-times and shift the I/O signals appropriately,estimate initial states, identify systems operating in closed-loop, or simulate discrete-time LTIsystems.The theoretical algorithms and their MATLAB realizations have been largely reorganized. Forinstance, the N4SID algorithm makes use of the inverse of a part of the triangular factor ofan RQ factorization of a Hankel-like matrix constructed from I/O sequences. (This is neededbefore computing the SVD that gives the system order and then the quadruple of systemmatrices.) The MATLAB implementations solve several standard least squares problems toobtain the associated, so-called “oblique projection”. While the new, LAPACK-based imple-mentation only involves some orthogonal transformations for obtaining the various projectionsand residuals.One objective of the Fortran implementation has been to use the LAPACK routines as muchas possible. Preference has generally been given to the block variants and the Level III BLAScodes, which are responsible for efficient use of parallelism. Advantage was taken of theproblem structure, whenever possible. For instance, the calculations have been organizedsuch that singular value decompositions were required for triangular matrices only, and adedicated, very compact and efficient routine has been written for that purpose. Moreover,both in theory and in the MATLAB implementations, the computation of the system matricesinvolves a very large matrix. Fortunately, this matrix can be brought to a special block-triangular structure. A structure-exploiting QR row compression scheme has been devised,reducing the storage requirements and computational effort.Another related objective has been to allow performing the computations with as reducedmemory requirements as possible, but without sacrificing the speed, when enough memory isavailable. For instance, there is an option for sequential calculation of the triangular factor inQR factorization of the Hankel-like matrix. The user can specify that I/O data are enteredin batches (independent or not). Moreover, when there is not enough working memory forprocessing a given data batch, the codes automatically perform an inner sequential processingof that batch, to accommodate with the available memory space. No provisions for usingthe displacement rank techniques (based on the block-Hankel structure) are provided in thecurrent implementations.2 Basic ApproachesThe basic LTI state space models considered are described byxk+1 = Axk +Buk + wk,yk = Cxk +Duk + vk, (1)where xk is the n-dimensional state vector at time k, uk is the m-dimensional input (control)vector, yk is the `-dimensional output vector, {wk} and {vk} are state and output disturbanceor noise sequences, and A, B, C, and D are real matrices of appropriate dimensions. Thesystem order, n, and the quadruple of system matrices (A,B,C,D) are not known; the only2available information is given by an upper bound, s, on n, and by the input and output datasequences, {uk} and {yk}, k = 1, . . . t.2.1 The State Intersection ApproachThe main feature of the SI class of SMI techniques is the determination of the state sequenceof the LTI system to be identified, or of an observer to reconstruct its state sequence, viathe intersection of the row spaces of the Hankel-like matrices constructed from “past” and“future” I/O data. The basic idea was first described in [4]. An extension of this idea ledto the N4SID algorithm developed in [2]. This algorithm identifies LTI state space modelsin the so-called innovation form, where the additive disturbance wk in (1) is generated byan innovation model of the form wk = Kvk, {vk} is a zero-mean white noise sequence and{uk} is a deterministic input sequence (perfectly known to the user). Both the basic andextended variants of this class produce statistically consistent and efficient estimates whencertain assumptions hold.2.2 The MOESP ApproachThe main feature of this class of SMI techniques is the determination of an extended observ-ability matrix, Γs, of the deterministic part of the model (1). Note that the upper bound, s,on n, is used. The basic idea was first described in [5] and [6]. The simplest algorithm de-rived from this idea, called the ordinary MOESP scheme, allows to identify, in a statisticallyconsistent and efficient way, LTI systems that can be described by (1), with wk ≡ 0 and {vk}a zero-mean white noise sequence, independent of the input.Extensions based on using past input quantities and/or reconstructed state variables as in-strumental variables have been proposed, allowing to consistently identify a model in (1), forwk ≡ 0, and {vk} a zero-mean arbitrary stochastic disturbance, independent of the input.The price paid for this increased applicability is that the estimates are no longer efficient.These variants are referred to as the PI or RS schemes, when past inputs or reconstructedstate variables are used as instrumental variables, respectively.When the additive disturbance wk in (1) is generated by an innovation model of the formwk = Kvk, and vk is a zero-mean white noise independent of the input, it is again possibleto obtain both consistent and efficient estimates when, besides the past input, past outputquantities are also used as instrumental variables. The scheme derived from this extension ofordinary MOESP scheme is known as the PO scheme [3].For the deterministic-stochastic identification problem, it is possible to estimate the deter-ministic part by using the PO scheme, and the stochastic part, that is the noise covariancematrices, by using the SI approach [2]. Further, the corresponding Kalman gain K can becomputed to allow the design of state observers [3].The algorithms belonging to the MOESP class have the striking feature of being highlystreamlined, in the sense that the sequence of computations performed by these schemes isalmost independent from the type of problems to be analyzed. This has led to highly modularimplementations of the algorithms within this class.33 LAPACK-Based SoftwareThe collection of implemented subroutines forms a library called RASP-IDENT [9] whichis part of the RASP package [13]. All implementations were done according to the RASP-SLICOT mutual compatibility concept introduced in [14]. The RASP-IDENT library consistsof six driver routines wich are compatible with the RASP implementation and documentationstandards, five auxiliary routines for system identification, ten SLICOT compatible driverroutines, and nineteen SLICOT compatible computational and auxiliary routines.3.1 Driver RoutinesThe following driver routines are available for MOESP class of SMI:RPIMOE computes consistent and statistically efficient estimates of the matrices (A,B,C,D)in (1) when w ≡ 0 and v is zero-mean white noise independent of the input, using theordinary MOESP scheme.RPIMPI computes consistent estimates of the matrices (A,B,C,D) in (1) when w ≡ 0 andv is zero-mean, but of arbitrary statistical color, using the ordinary MOESP schemeextended with instrumental variables based on past input quantities.RPIMPO computes consistent and statistically efficient estimates of the matrices (A,B,C,D)in (1) when w is generated by an innovation model and v is zero-mean white noise in-dependent of the input, using the ordinary MOESP scheme extended with instrumentalvariables based on past input and past output quantities.RPIMRS computes consistent estimates of the matrices (A,B,C,D) in (1) when w ≡ 0 andv is zero-mean, but of arbitrary statistical color, using the ordinary MOESP schemeextended with instrumental variables based on reconstructed state variables.RPIMKG performs the computations in RPIMPO routine, but also estimates the noisecovariance matrices and the Kalman gain (to allow the design of a state observer as astate predictor), using the SI approach.The following driver subroutine is available for the SI class of SMI:RPIMN4 computes consistent and statistically efficient estimates of the matrices (A,B,C,D)in (1) when w is generated by an innovation model and v is a zero-mean white noisesequence, using the basic N4SID algorithm. This routine also estimates the noise co-variance matrices and the Kalman gain (to allow the design of a state observer as astate predictor).All driver routines have an option for sequential I/O data processing.3.2 Computational and Auxiliary RoutinesTo validate and extend the use of various subspace identification schemes, a number of com-putational and auxiliary routines have been implemented. The main ones are listed below:IMMPID estimates a finite set of Markov parameters from I/O data (no constraints oninput sequence to be white noise, or initial conditions be zero).4IMSHFT shifts the I/O data sequences to compensate for the dead-times in the input vectorcomponents.IMICID estimates the initial state x0 of a LTI system, given the matrix quadruple (A,B,C,D)and the input and output trajectories of the system.IMCLID computes the parameters of a LTI system operating in closed-loop with a LTIcontroller.RPIMIU computes the output trajectory of a LTI discrete-time system, given the inputtrajectory, the initial state and/or the state and output disturbances (or noise) trajec-tories.4 ConclusionsLAPACK-based Fortran software for multivariable system identification by subspace tech-niques have been briefly described. The underlying approaches seem very promising, and thesoftware components proved robust, flexible, and easy of use, even when the available memoryis quite limited. The software has been extensively tested; the results obtained using variousalgorithms (including their MATLAB versions) have been compared. New versions of thecodes, with improved modularity and efficiency, are under development. Implementation ofnew algorithms is also considered.References[1] L. Ljung, System Identification Toolbox. Natick, MA: The MathWorks, Inc., 1992.[2] P. Van Overschee, and B. De Moor, N4SID: Two subspace algorithms for theidentification of combined deterministic-stochastic systems, Automatica, vol. 30, no. 1,pp. 75–93, 1994.[3] M. Verhaegen, Identification of the deterministic part of MIMO state space modelsgiven in innovations form from input-output data, Automatica, vol. 30, no. 1, pp. 61–74,1994.[4] M. Moonen, B. De Moor, L. Vandenberghe, and J. Vandewalle, On- and off-line identification of linear state-space models, Int. J. Control, vol. 49, no. 1, pp. 219–232,1989.[5] M. Verhaegen, and P. Dewilde, Subspace Model Identification. Part 1: The output-error state- space model identification class of algorithms, Int. J. Control, vol. 56, no. 5,pp. 1187–1210, 1992.[6] M. Verhaegen, Subspace Model Identification. Part 3: Analysis of the ordinary output-error state-space model identification algorithm, Int. J. Control, vol. 58, no. 3, pp. 555–586, 1993.[7] P. Van Overschee, Subspace Identification : Theory – Implementation – Applications,Katholieke Universiteit Leuven, Ph. D. Thesis, 1995.5[8] M. Verhaegen, State Space Model Identification Toolbox, Delft University of Technol-ogy, Tech. Rep., Nov. 1993.[9] V. Sima, and A. Varga, RASP–IDENT Subspace Model Identification Programs, DLR– Oberpfaffenhofen, Tech. Rep. TR R888–94, Oct. 1994.[10] S. Van Huffel, and J. Vandewalle, The Total Least Squares Problem: Computa-tional Aspects and Analysis. Philadelphia, PA: SIAM, 1991.[11] M. Verhaegen, and A. Varga, Some Experience with the MOESP Class of SubspaceModel Identification Methods in Identifying the BO105 Helicopter, DLR – Oberpfaffen-hofen, Tech. Rep. TR R165–94, June 1994.[12] V. Sima, Algorithms and LAPACK-Based Software for Subspace Identification, Pro-ceedings of The IEEE International Symposium on Computer-Aided Control SystemDesign—CACSD’96, September 15–18, 1996, Dearborn, Michigan, U.S.A. (to appear).[13] G. Gru¨bel, and H.-D. Joos, RASP and RSYST - Two Complementary Program Li-braries for Concurrent Control Engineering, Prepr. 5th IFAC/IMACS Symp. CADCS’91,Swansea, UK, Pergamon Press, Oxford, pp. 101-106, 1991.[14] G. Gru¨bel, A. Varga, A. van den Boom, and A. J. Geurts, Towards a coordinateddevelopment of numerical CACSD software: the RASP/SLICOT compatibility concept,Prepr. of IEEE/IFAC Symp. CACSD’94, Tucson, Arizona, 19946

Algorithms and LAPACK-Based Software for Subspace Identiﬁcation,

Identiﬁcation of the deterministic part of MIMO state space models given in innovations form from input-output data,

N4SID: Two subspace algorithms for the identiﬁcation of combined deterministic-stochastic systems,

On- and oﬀ-line identiﬁcation of linear state-space models,

Overschee, Subspace Identiﬁcation : Theory – Implementation – Applications, Katholieke Universiteit Leuven,

RASP–IDENT Subspace Model Identiﬁcation Programs,

Some Experience with the MOESP Class of Subspace Model Identiﬁcation Methods

Subspace Model Identiﬁcation. Part 1: The outputerror state- space model identiﬁcation class of algorithms,

Subspace Model Identiﬁcation. Part 3: Analysis of the ordinary outputerror state-space model identiﬁcation algorithm,

System Identiﬁcation Toolbox.

The Total Least Squares Problem: Computational Aspects and Analysis.

Towards a coordinated development

LAPACK-Based Software for Subspace Identification

LAPACK-Based Software for Subspace Identification

Abstract

Similar works

Full text

Available Versions

Institute of Transport Research:Publications