Identification of low-order models using Laguerre basis function expansions⁎

Proceedings,18th Proceedings,18th IFAC IFAC Symposium Symposium on on System System Identification Identification July Sweden Proceedings,18th IFAC System July 9-11, 9-11, 2018. 2018. Stockholm, Stockholm, Sweden on Proceedings,18th IFAC Symposium Symposium on System Identification Identification Available online at www.sciencedirect.com July Sweden Proceedings,18th IFAC Symposium July 9-11, 9-11, 2018. 2018. Stockholm, Stockholm, Sweden on System Identification July 9-11, 2018. Stockholm, Sweden

ScienceDirect

IFAC PapersOnLine 51-15 (2018) 72–77

Identification of low-order models using Identification Identification of of low-order low-order models models using using  Identification of low-order models using Laguerre basis function expansions Laguerre basis function expansions Identification of low-order models using  Laguerre Laguerre basis basis function function expansions expansions  Laguerre function expansions Mikael basis Manng˚ ard, Hannu T. Toivonen

Mikael Manng˚ ard, Hannu T. Toivonen Mikael ard, Mikael Manng˚ Manng˚ ard, Hannu Hannu T. T. Toivonen Toivonen Mikael Manng˚ ard, Hannu T. Toivonen Process Control Laboratory, Faculty of Science Process Control Laboratory, Faculty of Science and and Engineering, Engineering, ˚ ˚ Process Control Laboratory, Faculty of Science and Engineering, A bo Akademi University, FIN-20500 Turku ( A bo), Finland ˚ ˚ Process Control Laboratory, Faculty of Science and Engineering, A bo Akademi University, FIN-20500 Turku ( A bo), Finland ˚ ˚ A Akademi University, FIN-20500 Turku ( A bo), Finland Process Control Laboratory, Faculty of Science and Engineering, {mikael.manngard, hannu.toivonen}@abo.fi) ˚ ˚ Abo bo(e-mail: Akademi University, FIN-20500 Turku ( A bo), Finland (e-mail: {mikael.manngard, hannu.toivonen}@abo.fi) ˚ {mikael.manngard, hannu.toivonen}@abo.fi) Abo(e-mail: Akademi University, FIN-20500 Turku (˚ Abo), Finland (e-mail: {mikael.manngard, hannu.toivonen}@abo.fi) (e-mail: {mikael.manngard, hannu.toivonen}@abo.fi) Abstract: In In this this paper, paper, aa method method for for identification identification of of low-order low-order output-error output-error models models expressed expressed Abstract: Abstract: this method of low-order expressed in terms terms of of In Laguerre basisaa functions functions is identification presented. The The identification problem is ismodels formulated as aa Abstract: In this paper, paper, method for for identification of identification low-order output-error output-error models expressed in Laguerre basis is presented. problem formulated as in terms Laguerre basis is presented. The problem formulated as Abstract: thisoptimization paper, method for of identification low-order output-error expressed rank-constrained optimization problem in the the Laguerre Laguerre domain, which which can be beis solved efficiently in terms of of In Laguerre basisa functions functions is identification presented. The identification problem ismodels formulated as aa rank-constrained problem in domain, can solved efficiently rank-constrained optimization problem in the Laguerre domain, can be efficiently in terms of Laguerre basis functions is presented. The models identification problem as a using nuclear norm regularization. Since the identified models arewhich of low low order, minimal staterank-constrained optimization problem inthe the Laguerre domain, which can beissolved solved efficiently using nuclear norm regularization. Since identified are of order, aaformulated minimal stateusing nuclear norm regularization. Since the identified models are of low order, a minimal staterank-constrained optimization problem in the Laguerre domain, which can be solved efficiently space realization of the system can be obtained without the use of additional approximative using nuclear norm regularization. Since the identified models are of low order, a minimal statespace realization of the system can be obtained without the use of additional approximative space realization of the system be without the use additional approximative using nuclear norm Since the identified models order, a minimal statemodel-order reduction steps. space realization ofregularization. thesteps. system can can be obtained obtained without the are useofof oflow additional approximative model-order reduction model-order reduction steps. space realization of the system can be obtained without the use of additional approximative model-order reduction steps. © 2018, IFAC reduction (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. model-order steps. Keywords: System identification, Laguerre basis, basis, nuclear nuclear norm norm regularization regularization Keywords: System identification, Laguerre Keywords: System System identification, identification, Laguerre Laguerre basis, basis, nuclear nuclear norm norm regularization regularization Keywords: Keywords: System identification, Laguerre basis, nuclear norm regularization 1. INTRODUCTION INTRODUCTION reduce 1. reduce the the complexity complexity of of the the identified identified model. model. HankelHankel1. INTRODUCTION INTRODUCTION reduce the complexity of the identified Hankelnorm optimal model-order reduction has been 1. reduce the complexity of the identified model. Hankelnorm optimal model-order reduction has model. been presented presented norm optimal model-order reduction has been presented 1. INTRODUCTION reduce the complexity of the identified model. Hankelin (Glover, 1984), and Heuberger (1991) proposed an Identification of low-order models from input-output data norm optimal model-order reduction has been presented in (Glover, 1984), and Heuberger (1991) proposed an Identification of low-order models from input-output data in (Glover, 1984), and Heuberger (1991) proposed an norm optimal model-order reduction has been presented Identification of low-order models from input-output data iterative strategy for identification of orthonormal basis has been of general interest in the fields of systems en(Glover, 1984),forand Heuberger of (1991) proposedbasis an Identification of low-order models from input-output iterative strategy identification orthonormal has been of general interest in the fields of systemsdata en- in iterative strategy for identification of orthonormal in (Glover, 1984), and Heuberger (1991) proposed an has been of general interest in the fields of systems enfunction models, iterating between construction of basis Identification of low-order models from input-output data gineering and control. For example, when synthesizing strategy iterating for identification orthonormal has been of interest in the fields of synthesizing systems en- iterative function models, between of construction of basis gineering andgeneral control. For example, when function models, iterating between construction of basis iterative strategy for identification of orthonormal gineering and control. For example, when synthesizing functions based on prior knowledge of the system poles, has been of general interest in the fields of systems encontrollers from system models, the order of the controller function models, iterating between construction of basis gineering For example, when synthesizing controllersand fromcontrol. system models, the order of the controller functions based on prior knowledge of the system poles, functions based on prior of system models, iterating between construction of poles, basis controllers from system models, models, the order order of the the controller identification, and approximative model order reduction. gineering control. Forthe example, synthesizing is typically typicallyand determined by the order ofwhen the model. Thus, function functions based on prior knowledge knowledge of the the system poles, controllers from system the of controller identification, and approximative model order reduction. is determined by order of the model. Thus, identification, and approximative model order reduction. functions based on prior knowledge of the system poles, is typically determined by the order of the model. Thus, controllers from system models, the order of the controller for control purposes, it is often of interest to identify the and identification approximative ismodel orderpopular reduction. is determined byoften the order of thetomodel. Thus, fortypically control purposes, it is of interest identify the identification, another apSubspace system system another apidentification, and identification approximative ismodel orderpopular reduction. for control purposes, it is is often of interest tomodel. identify the Subspace is typically determined byoften the order of theto Thus, simplest model compatible with theinterest available information. for control purposes, it of identify the simplest model compatible with the available information. Subspace system identification is another popular approach for identification of low-order state-space models. Subspace system identification is another popular approach for identification of low-order state-space models. simplest model compatible with the available information. for control purposes, it is often of interest to identify the Standard approaches to identify low-order models have simplest compatible with thelow-order available models information. Standardmodel approaches to identify have Subspace proach for identification of low-order state-space models. system identification is another popular apHowever, since virtually all subspace methods rely on deproach for identification of low-order state-space models. However, since virtually all subspace methods rely on deStandard approaches to identify low-order models have simplest model compatible with the available information. either been beenapproaches to specify specify to theidentify model structure structure prior to idenidenStandard low-order prior models have riving either to the model to However, since virtually all subspace methods rely on deproach for identification of low-order state-space models. low-rank matrices from which the key subspaces can However, since virtually all subspace methods rely on can deriving low-rank matrices from which the key subspaces either been to specify the model structure prior to idenStandard approaches to identify low-order models have tification or to identify a high-order model and apply a either been model structure iden-a However, tification or to to specify identifythe a high-order modelprior and to apply riving low-rank matrices from which the key subspaces can since virtually all subspace methods rely on debe derived, and on singular-value decompositions to obriving low-rank matrices from which the key subspaces can be derived, and on singular-value decompositions to obtification or to identify a high-order model and apply a either been to specify the model structure prior to idenmodel order reduction strategy. Recent advances in sparse tification or reduction to identifystrategy. a high-order andinapply model order Recentmodel advances sparsea riving be derived, and on singular-value decompositions to oblow-rank matrices from which the key subspaces can tain low-rank approximations of the matrices (Verhaegen and on singular-value decompositions to obtainderived, low-rank approximations of the matrices (Verhaegen model order reduction strategy. Recent advances in sparsea be tification or reduction to identify a high-order model andin apply optimization have provided tools for developing developing methods model order strategy. Recent advances sparse optimization have provided tools for methods tain low-rank approximations of the matrices (Verhaegen be derived, and on singular-value decompositions to oband Dewilde, 1992; Van Overschee and De Moor, 1994), tain low-rank approximations of the matrices (Verhaegen and Dewilde, 1992; Van Overschee and De Moor, 1994), optimization have provided tools for developing methods model order reduction strategy. Recent advances in sparse that constrain constrainhave the provided model order order during the identification identification optimization toolsduring for developing methods tain that the model the and Dewilde, 1992; Van Overschee and De Moor, 1994), low-rank approximations of the matrices (Verhaegen output-error optimality cannot, in general, be guaranteed and Dewilde, 1992; Van Overschee and De Moor, 1994), optimality cannot, in general, be guaranteed that constrainhave the provided model order order during the identification identification optimization toolsduring for developing methods output-error process. that constrain the model the process. output-error optimality general, be guaranteed and Dewilde,and 1992; Vancannot, Overschee and Deand Moor, 1994), (Verhaegen and Hansson, 2014).in Fischer and Medvedev output-error optimality cannot, in Fischer general, be guaranteed (Verhaegen Hansson, 2014). Medvedev process. that constrain the model order during the identification process. (Verhaegen and Hansson, 2014). Fischer and Medvedev output-error optimality cannot, in general, be guaranteed (1998) have presented a Laguerre-domain subspace identiIn the the past past decades, decades, orthogonal orthogonal basis basis functions functions have have been been (Verhaegen and Hansson, 2014). Fischer subspace and Medvedev (1998) have presented a Laguerre-domain identiIn process. (1998) presented aa Laguerre-domain subspace identiand Hansson, 2014). Fischer and Medvedev In thefor past decades, orthogonal basis basis functions have have been fication method. Recently, the nuclear norm (Fazel et used for system identification. For example, example, King been and (Verhaegen (1998) have presented Laguerre-domain subspace identiIn the past decades, orthogonal functions ficationhave method. Recently, the nuclear norm (Fazel et al., al., used system identification. For King and fication method. Recently, the nuclear norm (Fazel et al., al., (1998) have presented a Laguerre-domain subspace identiused for system identification. For example, King and 2004; Recht et al., 2010) has been introduced to a subspace In the past decades, orthogonal basis functions have been Paraskevopoulos (1979) and Wahlberg (1991, 1994) have fication method. Recently, the nuclear norm (Fazel et used for system (1979) identification. For example, King have and 2004; Recht et al., 2010) has been introduced to a subspace Paraskevopoulos and Wahlberg (1991, 1994) 2004; Recht et al., 2010) has been introduced to a subspace fication method. Recently, the nuclear norm (Fazel et al., Paraskevopoulos (1979) and Wahlberg (1991, 1994) have identification framework to improve the low-rank approxiused for system identification. For example, King and the Laguerre and Kautz bases for identification of 2004; Recht et al., 2010) has been introduced to a subspace Paraskevopoulos (1979) and Wahlberg (1991, 1994) have identification framework to improve the low-rank approxiused the Laguerre and Kautz bases for identification of identification to the approxiRecht etframework al., 2010) been introduced to constraints a subspace used the Laguerre Laguerre and Kautz Kautz bases forVan identification of 2004; mations (Verhaegen and Hansson, 2014). Rank Paraskevopoulos (1979) and Wahlberg (1991, have single-input single-output systems; andfor Van den1994) Hof et et al. identification framework to improve improve the low-rank low-rank approxiused the and bases identification of mations (Verhaegen and has Hansson, 2014). Rank constraints single-input single-output systems; and den Hof al. mations (Verhaegen and Hansson, 2014). Rank constraints identification framework to improve the low-rank approxisingle-input single-output systems; and Van den Hof et al. have also been proposed by Grossmann et al. (2009b,a) used the Laguerre and Kautz bases for identification of (1995) and Ninness Gustafsson (1997) have presented and Hansson, 2014). Rank single-input single-output systems; and Vanhave den presented Hof et al. mations have also(Verhaegen been proposed by Grossmann et al.constraints (2009b,a) (1995) and Ninness and Gustafsson (1997) have also been proposed by Grossmann et al. (2009b,a) mations (Verhaegen and Hansson, 2014). Rank constraints (1995) and Ninness and Gustafsson (1997) have presented and Recht et al. (2010) for identification of finite-impulsesingle-input single-output systems; and Van den Hof et al. least-squares identification methods that use general orhave also been proposed by Grossmann et al. (2009b,a) (1995) and Ninness and Gustafsson havegeneral presented least-squares identification methods(1997) that use or- and Recht et al. (2010) for identification of finite-impulseand Recht et (2010) for identification finite-impulsealso(FIR) been proposed Grossmannof al. (2009b,a) least-squares identification methods that use use general or- have response approximate low-order state(1995) and bases Ninness Gustafsson (1997) have presented thonormal bases thatand generalize the natural, natural, Laguerre, and and Recht et al. al.models (2010) that forby identification ofet finite-impulseleast-squares identification methods that general orresponse (FIR) models that approximate low-order statethonormal that generalize the Laguerre, and response (FIR) models that approximate low-order and Recht et al. (2010) for identification of finite-impulsethonormal bases that generalize the natural, Laguerre, and space models. least-squares identification methods that use general orKauz bases. The generalized orthonormal bases allow prior (FIR) models that approximate low-order statestatethonormal generalize the natural, Laguerre, and response space models. Kauz bases.bases The that generalized orthonormal bases allow prior space models. (FIR) models that approximate low-order stateKauz bases.bases The generalized orthonormal bases allow prior prior thonormal that generalize theto Laguerre, and response information of the the pole locations tonatural, be incorporated incorporated into space models. Kauz bases. The generalized orthonormal bases allow information of pole locations be into In this paper, we generalize the procedure in Grossmann thismodels. paper, we generalize the procedure in Grossmann space information of the the pole locations locations to be beexpressed incorporated into In Kauz bases.identification. The generalized orthonormal bases allow prior the system system identification. Since models expressed in terms terms information of pole to incorporated into the Since models in In this paper, generalize procedure et al. (2009b,a) Recht et (2010) for identificaIn this paper, we we and generalize the procedure in Grossmann et al. (2009b,a) and Recht the et al. al. (2010) in forGrossmann identificathe system identification. Since models expressed in terms information of the pole locations to linear beexpressed incorporated into et al. (2009b,a) of orthonormal orthonormal basis functions are linear in parameters, parameters, the system identification. Since models in terms of basis functions are in and Recht et al. (2010) for In this paper, we generalize the procedure in Grossmann tion using the Laguerre basis. Since the Laguerre basis et al.using (2009b,a) and Recht etSince al. (2010) for identificaidentification the Laguerre basis. the Laguerre basis is is of orthonormal basis functions are linear in parameters, the system identification. Since models expressed in terms they are generally well suited for output-error identificaof orthonormal basis areoutput-error linear in parameters, they are generally wellfunctions suited for identifica- tion using the Laguerre basis. Since the Laguerre basis et al. (2009b,a) and Recht et al. (2010) for identificaable to capture exponentially decaying impulse responses tion using the Laguerre basis. Since the Laguerre basis is is able to capture exponentially decaying impulse responses they are generally well suited for output-error identificaof orthonormal basis functions are linear in parameters, tion (S¨ o derstr¨ o m and Stoica, 1988; Tomita et al., 1992) they are generally wellStoica, suited 1988; for output-error tion (S¨ oderstr¨ om and Tomita et identificaal., 1992) tion able to capture exponentially decaying impulse responses using the Laguerre basis. Since theof Laguerre basis is better than FIR models, the number model parameable to capture exponentially decaying impulse responses better than FIR models, the number of model parametion (S¨ derstr¨ mleast-squares and Stoica, 1988; Tomita formulations. et identificaal., 1992) 1992) they are generally wellStoica, suited type for output-error and allow allow linear least-squares type problem formulations. tion (S¨ ooderstr¨ oom and 1988; Tomita et al., and linear problem better FIR models, of model parameable to than capture impulse responses ters for identification is in smaller, thus, than FIRexponentially models, the the decaying number of model parameters needed needed for identification isnumber in general general smaller, thus, and allow linear least-squares type problem formulations. tion (S¨ o derstr¨ o m and Stoica, 1988; Tomita et al., 1992) However, least-squares identification strategies using or- better and allow linear least-squares type problem formulations. However, least-squares identification strategies using orters needed for identification is in general smaller, thus, better than FIR models, the number of model paramereducing the computational cost of solving the problem. ters needed for identificationcost is inofgeneral smaller, thus, reducing the computational solving the problem. However, least-squares identification strategies using orand allow linear least-squares type problem formulations. thonormal basis functions typically result in high-order However, identification strategies using or- ters thonormalleast-squares basis functions typically result in high-order reducing the computational cost of solving the problem. needed for identification is in general smaller, thus, Another advantage using the Laguerre basis compared to reducing the computational cost of solving the problem. Another advantage using the Laguerre basis compared to thonormal basis functions typically result in high-order However, least-squares identification strategies using ormodels and provide little information about the poles of thonormal functions typically result high-order models andbasis provide little information aboutinthe poles of reducing Another advantage using the Laguerre basis compared to the computational cost of solving the problem. the natural basis is that it allows some prior information Another advantage using the Laguerre basis compared to the natural basis is that it allows some prior information models and provide little information about the poles of thonormal basis functions typically result in high-order the underlying system. models and provide little information about the poles of the underlying system. the natural basis is that it allows some prior information Another advantage using the Laguerre basis compared to of the system poles to be included in the identification. thethe natural basis is that it included allows some prior information of system poles to be in the identification. the underlying system. models and provide little information about the poles of the underlying system. strategies have often been included the of the system poles to be included in the identification. natural basis is that it allows some prior information Furthermore, due to the specific structure of the Laguerre Model-order reduction of the system poles to be included in the identification. due to the specific structure of the Laguerre Model-order reduction the underlying system. strategies have often been included Furthermore, Furthermore, to specific structure the Laguerre thefunctions, systemdue poles to expressions be included in theof identification. Model-order reduction strategies have often often been been included basis functions, simple expressions relating the model paas an an additional additional step strategies in the the identification identification procedure to of Furthermore, due to the the specific structure ofthe themodel Laguerre Model-order reduction have included basis simple relating paas step in procedure to basis functions, simple expressions relating the model paFurthermore, due to the specific structure of the Laguerre as an additional step in the identification procedure to rameters, transform-domain impulse response, state-space Model-order reduction strategies have often been included  basis functions, simple expressions relating the model paas an additional step in the identification procedure to rameters, transform-domain impulse response, state-space M.M. gratefully gratefully acknowledge acknowledge the the financial financial support support from from the the  M.M. rameters, transform-domain impulse response, state-space basis functions, simple expressions relating the model parealizations, and model order can be derived. as an additional step in the identification procedure to  rameters, transform-domain impulse response, state-space realizations, and model order can be derived. Finnish Graduate School in Chemical Engineering (GSCE). M.M. gratefully acknowledge the financial support from the  Finnish in Chemical (GSCE). M.M. Graduate gratefullySchool acknowledge the Engineering financial support from the realizations, and model order can be derived. rameters, transform-domain impulse response, state-space  realizations, and model order can be derived. Finnish Graduate School in Chemical Chemical Engineering (GSCE). M.M. Graduate gratefullySchool acknowledge the Engineering financial support from the Finnish in (GSCE). realizations, and model order can be derived. Finnish Graduate School in Chemical Engineering (GSCE).

2405-8963 © IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 2018, 2018 IFAC 72 Copyright 2018 responsibility IFAC 72 Control. Peer review© under of International Federation of Automatic Copyright © 2018 IFAC 72 Copyright © 2018 IFAC 72 10.1016/j.ifacol.2018.09.093 Copyright © 2018 IFAC 72

2018 IFAC SYSID July 9-11, 2018. Stockholm, Sweden

Mikael Manngård et al. / IFAC PapersOnLine 51-15 (2018) 72–77

2. PROBLEM STATEMENT

73

such as the Laguerre or Kautz bases, is that a large number of basis functions is needed to approximate systems with exponentially decaying impulse responses. Thus, in this paper, we generalize the results presented by Grossmann et al. (2009b,a) to expansions using the Laguerre basis.

We consider the problem of identifying low-order models described by y(t) = G(q)u(t) + v(t), (1) where q is the forward shift operator, u(t) and y(t) are measured input and output signals at time instances t = 0, 1, 2, ... and {v(t)} is a sequence of noise. This model structure is referred to as the output-error model structure (Tomita et al., 1992).

The proposed method for solving Problem 1 can be outlined as follows: (1) Expand the system model in terms of a Laguerre basis (cf. (4)). (2) Estimate the model parameters {gk } subject to a rank constraint on the Laguerre domain Hankel matrix. (3) Find a minimal partial state-space realization in the Laguerre domain. (4) Compute a time-domain minimal state-space realization using the inverse Laguerre operator transform.

The relation between input and output signals of any linear, time-invariant (LTI) system G of finite order can be represented in a state-space form as x(t + 1) = Ax(t) + Bu(t) (2) y(t) = Cx(t) + Du(t), (3) where x(t) ∈ Rn is a state vector and (A, B, C, D) is the state-space realization. The system identification problem considered in this paper is defined as follows: Problem 1. Given a set of N input-output data pairs {u(t), y(t)}, identify a model G and a minimal statespace realization (A, B, C, D) of specified order n, which minimizes the output-error cost N −1  |y(t) − G(q)u(t)|2 .

In order to expand the system in terms of a Laguerre basis, the relation between the model parameters {gk }, the Laguerre-domain impulse-response coefficients, and the model order needs to be stated. Thus, in Section 3 we present key properties of the Laguerre basis and the signal and system transforms it induces, which are needed to pose a constraint on the model order in the Laguerre domain. In Section 4 we formulate Problem 1 as a rankconstrained minimization problem, which is relaxed using the nuclear norm to obtain a tractable convex optimization problem. The complete identification method proposed in this paper is presented in Section 4.

t=0

Minimization of the output-error cost in Problem 1 is in general non-linear in the parameters, resulting in nonconvex optimization problems. In subspace identification, e.g. (Verhaegen and Dewilde, 1992; Van Overschee and De Moor, 1994; Verhaegen and Hansson, 2014), one estimates the state sequence {x(k)} by performing a singular value decomposition, followed by a parameter estimation step. Another approach to the output-error identification problem has been to expand the system in terms of a complete orthonormal basis {Ψk (z)} in H2 as ∞  G(z) = gk Ψk (z). (4)

3. PROPERTIES OF THE LAGUERRE BASIS In this section, some properties of the Laguerre basis and of the transforms it induces are briefly reviewed. Following the convention of Nurges and Jaaksoo (1981), the discretetime domain Laguerre basis functions are defined as    k   t + k − j t+k−2j k+j k 2 α ψk (t) = 1 − α (−1) , k j j=0

(6) k ∈ {0, 1, 2, ...}, where |α| < 1, and the z-domain representations of the basis functions are given by (King and Paraskevopoulos, 1977) √  k z 1 − α2 1 − αz Ψk (z) = , k ∈ {0, 1, 2, ...}. (7) z−α z−α The parameter α corresponds to a real pole that occurs with multiplicity k + 1 in the basis functions. Observe that the natural basis {z −k } is a special case of the Laguerre basis corresponding to α = 0. Note also that the Laguerre basis is often defined as {1, z −1 Ψ0 (z), z −1 Ψ1 (z), ...}, thus avoiding a direct feed-through term in (4) (Heuberger et al., 2005). As the set of Laguerre basis functions {ψk (t)}k=0,1,2,... constitute a complete orthonormal basis in 2 , the elements of an 2 sequence {f (t)}t=0,1,2,... can be expressed in terms of a linear combination ∞  f (t) = fk ψk (t), (8)

k=0

By representing the system in terms of a basis function expansion, the output-error identification problem can be solved as a linear least-squares problem due to the linearin-parameter model structure. However, in practice a large number of terms may be required in the expansion for accurate representation of the underlying system, resulting in high-order models. Therefore, it is motivated to bound the model order during identification. Problem 1 was partially solved by Grossmann et al. (2009b,a) for the natural basis {Ψk (z)} = {z −k }. The key property they used to constrain the model order was that for a system expanded using the natural basis, the system order is given by the rank of the Hankel matrix   h1 h 2 h 3 · · · h2 h3 h4 · · ·  H= (5) h3 h4 h5 · · · .. .. .. . . . . . .

k=0

where {hk }k=1,2,.. is the impulse response of the system. This allows formulation of the identification problem as a rank-constrained optimization problem. A shortcoming of the natural basis compared to other orthonormal bases

with Laguerre expansion coefficients ∞  fk = f (t)ψk (t). t=0

73

(9)

2018 IFAC SYSID 74 July 9-11, 2018. Stockholm, Sweden

u(t)

uk

Mikael Manngård et al. / IFAC PapersOnLine 51-15 (2018) 72–77

G

˜ G

˜ k = C˜ A˜k−1 B, ˜ k = 1, 2, .... h (19) ˜ ˜ ˜ ˜ Furthermore, given a realization (A, B, C, D) in the Laguerre domain, a time-domain state-space realization (A, B, C, D) is obtained from Equations (14)-(17) by the reverse relations ˜ −1 (20) A = (A˜ + αI)(I + αA)  −1 ˜ 2 ˜ B = 1 − α (I + αA) B (21)  −1 2 ˜ ˜ (22) C = 1 − α C(I + αA) −1 ˜ ˜ ˜ ˜ (23) D = D − αC(I + αA) B.

y(t)

yk

Fig. 1. The Laguerre signal and system transform.

Thus, the Laguerre operator transform constitutes a one˜ B, ˜ C, ˜ D). ˜ to-one mapping between (A, B, C, D) and (A, Note that the one-to-one correspondence is due to the relatively simple structure of the Laguerre basis, and the result does not generalize as easily to other basis functions. ˜ B, ˜ C) ˜ is minTheorem 1. A Laguerre-domain system (A, imal if and only if the original time-domain system (A, B, C) is minimal.

The sequence of Laguerre coefficients {fk } is called the Laguerre signal transform of {f (t)}. This is a bijective mapping induced by the Laguerre basis (Heuberger, 1991). Thus, there also exists an inverse Laguerre signal transform, mapping signals yk in the Laguerre-domain to signals y(t) in 2 . Furthermore, given a system G ∈ H2 with input-output relation y(t) = G(q)u(t), by expanding the input u(t) and the output y(t) in terms of the Laguerre ˜ that describes the coefficients, there exists an operator G ˜ , input-output relation yk = G(λ)u k where λ is a shift operator in the Laguerre domain. Thus, the Laguerre basis also induces a Laguerre operator transform in H2 , ˜ mapping operators G in the time domain to operators G in the Laguerre domain, and an inverse Laguerre transform ˜ to G. The connection between the signals and mapping G operators is illustrated in Figure 1.

Proof. Follows from the one-to-one correspondence of the Laguerre transform. See e.g. Nurges and Jaaksoo (1981) or Heuberger (1991). From Equation (18) and Theorem 1 it directly follows that ˜ = rank (H). ˜ deg (G) = deg (G) (24) The following theorem establishes a connection between the model parameters gk in (4) and the Laguerre-domain ˜k. impulse response h Theorem 2. for a system expanded using Laguerre filters as in (4), the Laguerre-domain impulse-response coefficients are given by ˜0 = √ 1 ˜k = √ 1 h g0 , h (αgk−1 + gk ) , k ∈ N. 1 − α2 1 − α2

˜ is a linear and invariant causal system, To show that G consider a time-domain state-space system defined by Equations (2) and (3). In analogy with Nurges and Jaaksoo (1981), multiplying both sides with ψk (t) and summing over t = 0, 1, 2, ..., results in the following relation between the Laguerre coefficients: (10) x ¯k = Axk + Buk yk = Cxk + Duk , (11) ∞ ¯k + where x ¯k = t=0 x(t + 1)ψk (t). From the equality x α¯ xk+1 = αxk + xk+1 , and by defining a new state vector zk = (I − αA)xk − αBuk , the expansion coefficients yk can be expressed in a state-space form ˜ k + Bu ˜ k (12) zk = Az ˜ k + Du ˜ k, (13) yk = Cz

Proof. Consider the Laguerre-domain system yk = ˜ G(λ)u k . Let uk be an impulse at k = 0, i.e. u0 = 1 and uk = 0 for k = 0. Then, by the z-transform of (8) we have ∞  u(z) = uk Ψk (z) = Ψ0 (z), k=0

and the output response in the z-domain is given by ∞ ∞   gk Ψk (z)u(z) = gk Ψk (z)Ψ0 (z). y(z) =

where A˜ = (A − αI)(I − αA)−1  ˜ = 1 − α2 (I − αA)−1 B B  C˜ = 1 − α2 C(I − αA)−1 ˜ = D + αC(I − αA)−1 B. D

(14)

k=0

(15)

Since

k=0

k 1 − αz z−α  k+1 2 2 1 − αz z (1 − α ) = (z − α)(1 − αz) z − α  k+1  1 − αz 1 − αz + α(z − α) = z2 (z − α)(1 − αz) z−α  k+1  1 1 − αz α + =z z − α 1 − αz z−α α 1 =√ Ψk+1 (z) + √ Ψk (z) 2 1−α 1 − α2 and gk = 0 for k < 0, it follows that

(16)

z 2 (1 − α2 ) Ψk (z)Ψ0 (z) = (z − α)2

(17) Thus, the Laguerre operator transform of a system G described by (A, B, C, D) results in a causal linear system ˜ which is invariant with respect to the coefficient indices G, k = 0, 1, 2, .... From classical state-space theory, we have that for a ˜ B, ˜ C, ˜ D) ˜ with Laguerre-domain state-space realization (A, n×n ˜ , A∈R ˜ B, ˜ C, ˜ D) ˜ minimal ⇔ rank(H) ˜ = n, (A, (18) ˜ is the Hankel operator associated with the where H Laguerre impulse-response coefficients 74



2018 IFAC SYSID July 9-11, 2018. Stockholm, Sweden

Mikael Manngård et al. / IFAC PapersOnLine 51-15 (2018) 72–77

∞ ∞   αgk−1 g √ √ k y(z) = Ψk (z) + Ψk (z) 2 1−α 1 − α2 k=1 k=0

=

∞ 

k=0

√

where γ > 0 is a weight parameter, enforcing a tradeoff between the quality of fit and the model order. By solving (26) for various γ, Pareto-optimal solutions are obtained. Alternatively, a regularization criteria such as AIC or BIC can be incorporated in the identification by the SPARSEVA approach (Ha et al., 2015). To further promote sparsity of a solution, iterative reweighting (Mohan and Fazel, 2010; Ha et al., 2015) can be applied.

1 (αgk−1 + gk ) Ψk (z). 1 − α2

Thus, from the inverse z-transform and from Equation (8), ∞ ∞   1 √ y(t) = yk ψk (t) = (αgk−1 + gk ) ψk (t), 1 − α2 k=0 k=0

˜ k }k=0,1,2,...,m , such that Given a low-rank solution {h ˜ = n, a minimal partial state-space realization rank(H) ˜ B, ˜ C, ˜ D) ˜ of order n is for example obtained by applying (A, ˜k} the Ho-Kalman algorithm (Ho and Kalman, 1966) on {h as in De Hoog et al. (2002), or by the optimal Hankel norm partial realization method presented by Halikias et al. ˜ B, ˜ C, ˜ D), ˜ a minimal time(2010). Given a realization (A, domain state-space realization (A, B, C, D) of order n is given by Equations (20)-(23). The complete identification method proposed in this paper is summarized as follows:

for t = 0, 1, 2, .... Since the Laguerre-domain transfer function is invariant with respect to the coefficient indices k = 0, 1, 2, ..., the impulse response in the Laguerre˜ k = yk = √ 1 (αgk−1 + gk ) for domain is given by h 1−α2 all k ≥ 0, using the fact that gk = 0 for k < 0.  In Heuberger et al. (2005), an alternative proof for Theorem 2 is given, and the result is generalized for general orthonormal basis functions. 4. RANK-CONSTRAINED IDENTIFICATION

Algorithm: Proposed identification method. Define the Laguerre function pole |α| < 0, the number of Laguerre coefficients m, regularization parameter γ > 0, and the maximum model order s ≤ (m + 1)/2.

By Equation (24), Problem 1 can be formulated as a rank constrained optimization problem N −1 m    2 y(t) − gk Ψk (q)u(t) minimize g∈Rm+1 (25) t=0 k=0 subject to

1. Construct a Laguerre basis {Ψk (z)} as in (7). ˜ k }. 2. Solve optimization problem (26) to obtain {h 3. Compute a minimal partial state-space realization ˜ k }. ˜ B, ˜ C, ˜ D) ˜ based on {h (A, 4. Compute the minimal time-domain state-space realization (A, B, C, D) by Equations (20)-(23).

˜ s,m ) ≤ n, rank(H

˜ s,m is defined as where the truncated Hankel matrix H   ˜1 h ˜ m−s+1 ˜2 · · · h h h ˜ m−s+2  ˜2 h ˜3 · · · h  ˜ s,m =  H . . . .. . .. ..   .. . ˜ ˜ ˜ hs hs+1 · · · hm Here, the Laguerre-domain impulse-response coefficients ˜ k are given by Theorem 2, and s ≤ (m + 1)/2 can be h interpreted as an upper bound on the model order. Remark 1. In order to identify a system without a direct feedthrough, the minimization problem (25) can be solved m g subject to the linear constraint k=0 k ψk (0) = 0, where √ ψk (0) = 1 − α2 (−α)k (Nurges and Jaaksoo, 1981).

5. SIMULATION EXAMPLE In order to demonstrate how the proposed nuclear-norm regularized Laguerre-domain identification method performs, we give results for a simulated benchmark example that was originally introduced by Van den Hof and Ninness (2005). The proposed method is compared to MATLABs output-error routine OE, and to a standard method where a high-order Laguerre model is first identified by minimizing the least-squares cost, followed by a Hankel-norm model-order reduction step to find a low-order approximation of the system.

If a globally optimal solution to the optimization problem in (25) could be found, the proposed identification method would be output-error optimal in an 2 -norm sense. However, rank constrained optimization problems are in general NP-hard, and in many cases intractable. Thus, as a convex relaxation, the rank function is often relaxed with its convex envelope, the nuclear norm (Recht et al., 2010). Definition 1. The nuclear norm of a matrix X ∈ Rm×n , denoted by  · ∗ , is defined as

The simulated system has the structure as in (1), where G is a fifth-order system b1 z −1 + ... + b5 z −5 G(z) = , (27) 1 + a1 z −1 + ... + a5 z −5 with numerator coefficients b1 = 0.2530, b2 = −0.9724, b3 = 1.4283, b4 = −0.9493 and b5 = 0.2410, and denominator coefficients a1 = −4.15, a2 = 6.8831, a3 = −5.6871, a4 = 2.3333 and a5 = −0.3787. The system poles are located at 0.95 ± 0.20i, 0.85 ± 0.10i and 0.55 (rounded to two decimals). Colored noise v(k) = H(z)e(k) is added to the output, where {e(k)} is a zero-mean white-noise sequence with variance 0.0025 and the transfer function 1 − 1.38z −1 + 0.40z −2 H(z) = , (28) 1 − 1.90z −1 + 0.91z −2 which has a resonant peak at 0.1 rad/s. The input sequence {u(t)} is a zero-mean unit-variance white-noise sequence.

min{m,n}

X∗ :=



σi (X),

i=1

where σi (X) is the ith singular value of X. Introducing the nuclear norm as a convex relaxation in (25), we obtain the convex optimization problem N −1 m   ˜ s,m ∗ , (26) minimize γ |y(t)− gk Ψk (q)u(t)|2 +H m+1 g∈R

t=0

75

k=0

75

2018 IFAC SYSID 76 July 9-11, 2018. Stockholm, Sweden

Mikael Manngård et al. / IFAC PapersOnLine 51-15 (2018) 72–77

Magnitude (abs)

This results in a signal-to-noise ratio at the output of approximately 11.4 dB. To find an appropriate value for γ in (26), models using various weights γ ∈ [0.1, 100] were identified from a single data set of length N = 1200 using iterative reweighting. The root-mean-square error (RMSE) and the McMillan degree of the Pareto optimal models are presented in Figure 2. It is evident that for an appropriately selected γ > 2, the model order is the same as for the true model, and the RMSE equals the standard deviation of the noise. Hence it can be concluded that the proposed identification method can achieve output-error optimality.

10

10

RMSE

Magnitude (abs)

10

2

6

Model order

0.4 0.2 10 -1

10 0

Fig. 3. Bode amplitude plot of the true system (black), and of the 100 Laguerre models (gray) identified using the proposed identification method. Standard deviation error bars are indicated in red.

0.1

1

0.6

Frequency (rad/s)

0.15

0

0.8

0 10 -2

0.2

0.05 10 -1

1

5

1 0.8 0.6 0.4 0.2 0 10 -2

4 3

10 -1

10 0

Frequency (rad/s)

2 1 10 -1

10 0

10 1

Fig. 4. Bode amplitude plot of the true system (black), and of the 100 Laguerre models (gray) identified by the Hankel-norm reduction method. Standard deviation error bars are indicated in red.

10 2

γ

Magnitude (abs)

Fig. 2. Pareto-optimal soulutions for various γ > 0. The true standard deviation of {v(t)} is indicated with a dashed line. To evaluate the performance of the three different identification methods, a set of 100 different input-output data realizations of length N = 1200 were generated. The weighting parameter was set to γ = 5, and a truncated Laguerre basis consisting of m = 20 basis functions, with poles located at α = 0.7, was used in all experiments. Note that the pole selection results in a frequency weighting of the variance of the results, see e.g. Heuberger et al. (2005), Sec. 4.3.3. To promote solutions to be of low rank, 20 iterations of reweighing, as proposed by Mohan and Fazel (2010), where introduced to (26).

1 0.8 0.6 0.4 0.2 0 10 -2

10 -1

10 0

Frequency (rad/s)

Fig. 5. Bode amplitude plot of the true system (black), and of the 100 models (gray) identified using the OE routine in MATLAB. Standard deviation error bars are indicated in red

Out of the 100 different models identified with the proposed method from the different data sets, 96 had the correct model order n = 5. Frequency responses of the identified models are presented in Figures 3-5 along with standard deviation error bars. The Laguerre models appear to capture the resonant behavior at frequencies close to 0.2 rad/s better than the models obtained with the OE method (approximatively 50% of the OE models failed to capture the resonance peak). Note that the apparent high variance close to 0.2 rad/s results from only a few (four) models being overfitted to the data. The variance between different frequency responses, for frequencies higher than 0.3 rad/s, is smaller for models estimated with the proposed method than for the models estimated with the Hankel-norm reduction method. This behavior can be explained by the fact that the colored noise excites higher frequencies less than lower frequencies, whereas the

white-noise input signal excites all frequencies, and that the proposed identification method is (close to) outputerror optimal in an 2 -norm sense, whereas the Hankelnorm reduction procedure used to determine the models in Figure 4 treats all frequencies similarly. The accuracy of the identified models can be quantified ˆ H , where G denotes the true by the H2 norm G − G 2 ˆ model and G an identified model. The mean H2 error was 0.0190 ± 0.0062 for the proposed identification method, 0.0434 ± 0.0168 for the Hankel-norm reduction method, and 0.0340 ± 0.0787 for the OE method. The mean static gain for the models identified with the proposed method was 0.9949 ± 0.0423, 0.9925 ± 0.0454 for the Hankel-norm reduction method, and 1.0249±0.0499 for the OE method. 76

2018 IFAC SYSID July 9-11, 2018. Stockholm, Sweden

0

-1 -1

0 Re(z)

1

0

-1 -1

0 Re(z)

1

0

-1 -1

0 Re(z)

77

Halikias, G., Tsoulkas, V., Pantelous, A., and Milonidis, E. (2010). Hankel-norm approximation of FIR filters: a descriptor-systems based approach. International Journal of Control, 83(9), 1858–1867. Heuberger, P.S.C., Van den Hof, P.M., and Wahlberg, B. (2005). Modelling and identification with rational orthogonal basis functions. Springer. Heuberger, P.S.C. (1991). On approximate system identification with system based orthonormal functions. Ph.D. thesis, Delft University of Technology. Ho, B. and Kalman, R.E. (1966). Effective construction of linear state-variable models from input/output functions. Regelungstechnik, 14(1-12), 545–548. King, R. and Paraskevopoulos, P. (1977). Digital Laguerre filters. International Journal of Circuit Theory and Applications, 5(1), 81–91. King, R. and Paraskevopoulos, P. (1979). Parametric identification of discrete-time siso systems. International Journal of Control, 30(6), 1023–1029. Mohan, K. and Fazel, M. (2010). Reweighted nuclear norm minimization with application to system identification. In Proceedings of the American Control Conference (ACC), 2953–2959. Ninness, B. and Gustafsson, F. (1997). A unifying construction of orthonormal bases for system identification. IEEE Transactions on Automatic Control, 42(4), 515– 521. ¨ and Jaaksoo, U. ¨ (1981). Laguerre state equaNurges, U. tions of a multivariable discrete time system. IFAC Proceedings Volumes, 14(2), 1153–1158. Recht, B., Fazel, M., and Parrilo, P.A. (2010). Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM Review, 52(3), 471– 501. S¨oderstr¨om, T. and Stoica, P. (1988). System Identification. Prentice-Hall, Inc. Tomita, Y., Damen, A.A., and Van den Hof, P. (1992). Equation error versus output error methods. Ergonomics, 35(5-6), 551–564. Van den Hof, P. and Ninness, B. (2005). System identification with generalized orthonormal basis functions. In Modelling and Identification with Rational Orthogonal Basis Functions, 61–102. Springer. Van den Hof, P.M., Heuberger, P.S.C., and Bokor, J. (1995). System identification with generalized orthonormal basis functions. Automatica, 31(12), 1821–1834. Van Overschee, P. and De Moor, B. (1994). N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica, 30(1), 75–93. Verhaegen, M. and Dewilde, P. (1992). Subspace model identification part 1. the output-error state-space model identification class of algorithms. International Journal of Control, 56(5), 1187–1210. Verhaegen, M. and Hansson, A. (2014). Nuclear norm subspace identification (N2SID) for short data batches. IFAC Proceedings Volumes, 47(3), 9528–9533. Wahlberg, B. (1991). System identification using Laguerre models. IEEE Transactions on Automatic Control, 36(5), 551–562. Wahlberg, B. (1994). System identification using Kautz models. IEEE Transactions on Automatic Control, 39(6), 1276–1282.

1 Im(z)

1 Im(z)

Im(z)

1

Mikael Manngård et al. / IFAC PapersOnLine 51-15 (2018) 72–77

1

Fig. 6. True poles (cross) and pole estimates obtained with the proposed method (left), Hankel-norm reduction method (center), and the OE method (right). The estimated poles for all identified models are presented in Figure 6. We can conclude that the proposed identification method outperforms the other methods on this benchmark problem. Note that when using MATLABs OE routine and the Hankel-norm reduction method, the correct model structure was specified prior to identification, whereas using the proposed method, the correct model order was identified from the data. 6. CONCLUSIONS A framework for identification of low-order output-error models using the Laguerre basis has been presented. We have shown that properties of the signal and system transforms induced by the basis can be used to constrain the model order during identification. We have also shown that a minimal state-space realization of the same order as the identified Laguerre models can be obtained. Future work could be directed towards extending the presented identification method for systems represented in terms of general orthonormal basis functions. REFERENCES De Hoog, T.J., Szab´ o, Z., Heuberger, P.S.C., Van den Hof, P.M., and Bokor, J. (2002). Minimal partial realization from generalized orthonormal basis function expansions. Automatica, 38(4), 655–669. Fazel, M., Hindi, H., and Boyd, S. (2004). Rank minimization and applications in system theory. In Proceedings of the American Control Conference, volume 4, 3273–3278. Fischer, B.R. and Medvedev, A. (1998). Laguerre shift identification of a pressurized process. In American Control Conference, 1998. Proceedings of the 1998, volume 3, 1933–1937. IEEE. Glover, K. (1984). All optimal hankel-norm approximations of linear multivariable systems and their L∞ -error bounds. International Journal of Control, 39(6), 1115– 1193. Grossmann, C., Jones, C.N., and Morari, M. (2009a). System identification via nuclear norm regularization for simulated moving bed processes from incomplete data sets. In Proceedings of the Conference on Decision and Control (CDC), 4692–4697. IEEE. Grossmann, C., Jones, C.N., and Morari, M. (2009b). System identification with missing data via nuclear norm regularization. In Proceedings of the European Control Conference (ECC), 448–453. IEEE. Ha, H., Welsh, J.S., Blomberg, N., Rojas, C.R., and Wahlberg, B. (2015). Reweighted nuclear norm regularization: A sparseva approach. IFAC-PapersOnLine, 48(28), 1172–1177. 77