Nuclear Norm Subspace Identification Of Continuous Time State-Space Models

5th International Conference on Advances in Control and 5th International Conference on Advances in Control and Optimization of Dynamical Systems 5th International Conference on Advances in Control Available onlineand at www.sciencedirect.com Optimization of Dynamical Systems 5th International Conference on Advances in Control and February 18-22, 2018. Hyderabad, India Optimization of Dynamical Systems February 18-22, 2018. Hyderabad, India Optimization of Dynamical Systems February 18-22, 2018. Hyderabad, India February 18-22, 2018. Hyderabad, India

ScienceDirect

IFAC PapersOnLine 51-1 (2018) 530–535

Nuclear Norm Subspace Identification Of Nuclear Norm Subspace Identification Of Nuclear Norm Subspace Identification Of Continuous Time State-Space Models Nuclear Norm Subspace Identification Of Continuous Time State-Space Models Continuous Time State-Space Models Continuous Time State-Space Models ∗ ∗

Santhosh Santhosh Kumar Kumar Varanasi, Varanasi, ∗∗ Phanindra Phanindra Jampana Jampana ∗∗ Santhosh Kumar Varanasi, ∗∗ Phanindra Jampana ∗∗ Santhosh Kumar Varanasi, Phanindra Jampana ∗ ∗ ∗ Department of of Chemical Chemical Engineering, Engineering, Indian Indian Institute Institute of of Technology Technology ∗ Department Hyderabad, Kandi, Sangareddy, Telangana502285, India (e-mail: Department of Chemical Engineering, Indian Institute of Technology ∗ Hyderabad, Kandi, Sangareddy, Telangana- 502285, India (e-mail: Department of Chemical Engineering, Indian502285, Institute of Technology [email protected],[email protected]). Hyderabad, Kandi, Sangareddy, TelanganaIndia (e-mail: [email protected],[email protected]). Hyderabad, Kandi, Sangareddy, Telangana502285, India [email protected],[email protected]). (e-mail: [email protected],[email protected]). Abstract: Subspace identification Abstract: Subspace identification techniques techniques derive derive approximate approximate models models rather rather than than models models that are optimal with respect to a goodness of fit criterion. To obtain low rank models, aa Abstract: Subspace identification techniques derive approximate models rather than models that are optimal withidentification respect to atechniques goodness derive of fit criterion. To obtain low rank models, Abstract: Subspace approximate models rather than models nuclear norm that are optimal with respect to a goodness of fit criterion. To obtain low rank models, a minimization method for estimating the system matrices of linear time invariant nuclear minimization method for estimating thecriterion. system matrices of linear timemodels, invariant that arenorm optimal with respect to afor goodness of fit To obtain low is rank a continuous time state-space models in the presence of measurement noise proposed. In nuclear norm minimization method estimating the system matrices of linear time invariant continuous time state-spacemethod modelsforinestimating the presence of measurement noise is time proposed. In nuclear norm minimization the system matrices of linear invariant continuous time state-space models in the presence of measurement noise is proposed. In the proposed approach, Generalized Poisson Moment (GPMF) method is to the proposed approach, Generalized Poisson Moment Functional Functional (GPMF) method is used used In to continuous time state-space models in which the presence of in measurement noise is proposed. circumvent the time-derivative problem is inherent continuous time models. To make the proposed approach, Generalized Poisson Moment Functional (GPMF) method is used to circumvent theapproach, time-derivative problem whichMoment is inherent in continuous timemethod models.isTo make the proposed Generalized Poisson Functional (GPMF) used to the proposed algorithm consistent, instrumental variables (Hankel matrix of past inputs) are circumvent the time-derivative problem which is inherent in continuous time models. To make the proposed algorithm consistent, instrumental variablesin(Hankel matrix of models. past inputs) are circumvent the time-derivative problem which is inherent continuous time To make considered. The accuracy of method is the help numerical the proposed algorithm instrumental (Hankel with matrix pastof inputs) are considered. The accuracyconsistent, of the the proposed proposed methodvariables is demonstrated demonstrated theof help numerical the proposed algorithm consistent, instrumental (Hankel with matrix pastof are considered. The accuracy ofsystems. the proposed methodvariables is demonstrated with theofhelp ofinputs) numerical simulations on a variety of simulations on a accuracy variety ofofsystems. considered. The the proposed method is demonstrated with the help of numerical simulations on a variety of systems. simulations a variety ofFederation systems. of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. © 2018, IFACon(International Keywords: Keywords: System System Identification, Identification, Nuclear Nuclear Norm Norm Minimization, Minimization, Continuous Continuous Time, Time, State-Space State-Space Keywords: System Identification, Nuclear Norm Minimization, Continuous Time, State-Space Models, Generalized Poisson Moment Functionals, Semi-Definite Programming. Models, Generalized Poisson Moment Functionals, Semi-Definite Programming. Keywords: System Identification, Nuclear Norm Minimization, Continuous Time, Models, Generalized Poisson Moment Functionals, Semi-Definite Programming. State-Space Models, Generalized Poisson Moment Functionals, Semi-Definite Programming. 1. INTRODUCTION eralized Poisson Moment Functionals (GPMF) were used 1. eralized Poisson Moment (GPMF) were used 1. INTRODUCTION INTRODUCTION eralized Poisson Moment Functionals Functionals (GPMF) used previously in estimation of state-space modelswere [Merc` ere 1. INTRODUCTION previously in estimation of state-space models [Merc` eere eralized Poisson Moment Functionals (GPMF) were used previously of etstate-space [Merc` re et al. 2007].inInestimation [Haverkamp al. 1996], amodels Bi-linear transSubspace Identification (SID) methods for Linear Time et al. 2007]. [Haverkamp al. 1996], aamodels Bi-linear transpreviously inIn estimation of et state-space [Merc` ere s−a Subspace Identification (SID) methods for Linear Time et al. 2007]. In [Haverkamp et al. 1996], Bi-linear transformation w = has been applied prior to inversion into s−a Subspace (LTI) Identification (SID) methods for Linear Invariant state-space models in both discreteTime and formation s+a w = s−a has beenet applied priora to inversion into et al. 2007]. [Haverkamp al. 1996], Bi-linear transs+a Invariant models in discrete and Subspace (LTI) Identification (SID) methods for Linear s+a formation w In = s−a has beennoise applied prior to inversion into time domain for better attenuation. The identifiInvariant (LTI) state-space models in both both discrete and the continuous time state-space has received lot of attention in theTime past s−a s+a the time domain for better noise attenuation. The identififormation w = s+a has been applied prior to inversion into continuous time has received lot of attention in the past Invariantyears (LTI) state-space in both discrete and the cation algorithm behaves if the filtering parameter time domain for betterpoorly noise attenuation. The identificontinuous time hasOverschee receivedmodels lot of De attention in the Verpast twenty [Van and Moor 1994, twenty years [Van Overschee and De Moor 1994, Vercation algorithm behaves poorly if filtering parameter the time domain for better noise attenuation. The identificontinuous hasOverschee received lot of techniques attention the Verpast (a) is selected close to that of poles of system [Ohta 2011]. algorithm behaves poorly if the the filtering parameter twenty years [Van and De Moor in 1994, haegen andtime Verdult 2007]. These derive ap- cation (a) is selected close to that of poles of system [Ohta 2011]. cation algorithm behaves poorly if the filtering parameter haegen and Verdult 2007]. These techniques derive aptwenty years [Van Overschee and De Moor 1994, Ver(a) is selected close to that of poles of system [Ohta 2011]. haegen and Verdult 2007]. These techniques derive approximate models rather than models that are optimal As SID methods rely on SVD for low rank approximation, (a) is selected close to that of poles of system [Ohta 2011]. proximate models rather than models that are optimal haegen and Verdult 2007]. These techniques derive approximate models rather of than are optimal SID methods rely on SVD for low rank approximation, with respect to goodness fit models criterionthat [Verhaegen and As it is difficult to extend these existing methods to include As SID methods rely on SVD for low rank approximation, with respect to goodness of fit criterion [Verhaegen and proximate models rather than models that are optimal it isSID difficult to extend these existing methods to include with respect to goodness of fit criterion [Verhaegen and As Hansson 2016]. This is on methods rely on SVDproperties for low rank approximation, Hansson 2016]. This approximation approximation is based based on linear linear constraints on model errors, such is difficult extend these existing methods to include with respect to goodness of fit criterion [Verhaegen and it constraints onto model errors, properties such as as dissipativdissipativHansson 2016]. This approximation based onHankel linear transformations and factorizations of isstructured it is difficult to extend these existing methods to include ity and missing input or output data. The current work constraints on model errors, properties such as dissipativtransformations and factorizations of structured Hankel Hansson 2016]. This approximation is based on linear and missing inputerrors, or output data. such The as current work transformations and factorizations of structured matrices constructed from the input-output data. Hankel Unlike ity constraints on model properties dissipativfocuses on incorporating rank minimization into the goodity and missing input or output data. The current work matrices constructed from the input-output data. Unlike transformations and factorizations of structured Hankel focuses on incorporating rank minimization into the goodmatrices constructed from the[Ljung input-output the prediction error methods 1999] indata. whichUnlike non- ity and missing input or output data. The current work focuses rank minimization the goodness of on fit incorporating criteria and formulating a single into multi-criteria the prediction error methods [Ljung 1999] in which nonmatrices constructed from the input-output data. Unlike the prediction errorismethods [Ljung 1999] in which non- focuses ness of fit criteria and formulating a single multi-criteria linear optimization required, subspace methods estimate on incorporating rank minimization into the goodlinear optimization ismethods required,[Ljung subspace methods estimate optimization problem. ness of fit criteria and formulating a single multi-criteria the system prediction erroris 1999] in which nonoptimization problem. linear optimization required, subspace methods estimate the matrices with the help of projection ideas and ness of fit criteria and formulating a single multi-criteria optimization problem. the system matrices with the help of projection ideas and linear optimization is required, subspace methods estimate the system matrices thesingular help of projection ideas and optimization computational tools with such as value decomposition Starting from the problem. from the work work of of [Fazel [Fazel et et al. al. 2001], 2001], a a large large computational tools such value decomposition the system matrices with thesingular help of projection ideas and Starting computational tools such as as singular value decomposition Most of these these number of recent are made integrate the Starting from thedevelopments work of [Fazel et al. to 2001], a large (SVD) and RQ factorization. techniques (SVD) and RQ factorization. Most of techniques number of recent developments are made to integrate the computational tools such as singular value decomposition Starting from the work of [Fazel et al. 2001], a large (SVD) anddeveloped RQ factorization. Most of these techniques number low rank approximation step of SID with the goodness of of recent developments are made to integrate the have been for discrete time (DT) models. have been developed for discrete time (DT) models. low rank approximation step of SID with the goodness of (SVD) anddeveloped RQ factorization. these techniques number ofapproximation recent developments are made to integrate the have been for discreteMost time of (DT) models. fit in discrete time. As rank minimization is an NP-hard low rank step of SID with the goodness of fit inrank discrete time. As rank minimization is an NP-hard Identification of Continuous time (CT) models has gained have been developed for discrete time (DT) models. low approximation step of SID with the goodness of fit in discrete time. relaxation As rank minimization is anminimizaNP-hard Identification of Continuous time (CT) models has gained problem, a convex of nuclear norm Identificationinofrecent Continuous time (CT) models advantages has gained fit problem, a convex relaxation of nuclear norm minimizamomentum years because of several in discrete time. As rank minimization is an NP-hard momentum in recent years because of several advantages tion (NNM) is For discrete a convex relaxation nuclear time normstate-space minimizaIdentification Continuous time (CT) models has gained problem, tion (NNM) is considered. considered. Forof time momentum recent because of several advantages compared toinof that of years DT models [Rao and Unbehauen problem, a convex relaxation of discrete nuclear normstate-space minimizamodels [Verhaegen and Hansson 2016] proposed the tion (NNM) is considered. For discrete time state-space compared to that of DT models [Rao and Unbehauen momentum in recent years because of several advantages [Verhaegen and Hansson 2016] proposed the NuNucompared to that of ofDTCTmodels [Rao and Unbehauen 2006]. Identification state-space models can be models tion (NNM) is considered. For discrete time state-space clear Norm Subspace Identification (N2SID) [Verhaegen models [Verhaegen and Hansson 2016] proposed the Nu2006]. Identification of CT state-space models can be compared to that of DT models [Rao and Unbehauen clear Norm Subspace Identification (N2SID) [Verhaegen 2006]. Identification of categories CT state-space models can be broadly divided into two 1) Indirect methods, in models [Verhaegen and Hansson 2016] proposed the NuclearHansson Norm Subspace Identification (N2SID) [Verhaegen broadly divided into two categories 1) Indirect methods, in and 2016] this technique. 2006]. aIdentification of categories CT state-space models can be and Hansson 2016] using using this convex convex relaxation relaxation technique. broadly divided intoistwo 1) Indirect methods, in clear which DT model identified first and then converted to Norm Subspace Identification (N2SID) [Verhaegen which a DT model istwo identified first and then converted to The optimization solved using SemiHansson 2016] problem using thiswas convex relaxation technique. broadly categories in and The optimization solved using both both Semia divided DT model identified first1) and then converted to awhich CT model andinto 2)isDirect method, inIndirect which amethods, CT model and Hansson 2016] problem using(SDP) thiswas convex relaxation technique. Definite Programming as well as the Alternating The optimization problem was solved using both Semiaais CT model and 2) Direct method, in which a CT model which a DT model is identified first and then converted to Definite Programming (SDP) as well as the Alternating CT model and 2) Direct method, in which a CT model directly identified from the sampled input-output data. optimization wasassolved both Semidata. The Direction Method problem of Multipliers (ADMM) et al. Definite Programming (SDP) well asusing the [Fazel Alternating is directly identified from the sampled input-output a CT model and 2) Direct method, in which a CT model is directly identified from the sampled input-output data. Direction Method of (ADMM) [Fazel et The indirect method suffers from many limitations as (SDP) as well as of thethis Alternating The indirect method from suffers from many input-output limitations such such as Definite 2013, LiuProgramming et al. 2013]. Another advantage approach Direction Method of Multipliers Multipliers (ADMM) [Fazel et al. al. is directly identified the sampled data. The indirect and method from many limitations such as Direction 2013, Liu et al. 2013]. Another advantage of this approach [Unbehauen Raosuffers 1990] a) Difficult choice of sampling Method of Multipliers (ADMM) [Fazel et by al. is that the missing input-output data can be estimated 2013, Liu et al. 2013]. Another advantage of this approach [Unbehauen and Rao 1990] a) Difficult choice of sampling The indirect method from many limitations such as is that the missing input-output data can be estimated by [Unbehauen and Raosuffers 1990] a) Difficult choice of sampling period for stiff systems and b) Translating the zeros to CT 2013, Liu et al. 2013]. Another advantage of this approach formulating a matrix completion data problem et al. 2013]. is that the missing input-output can [Liu be estimated by period for stiff systems and the˚zeros to CT [Unbehauen and Rao 1990] a) Translating Difficult choice of sampling period for stiff systems and b) b) Translating aa matrix completion problem [Liu et al. 2013]. is that the missing input-output data can be estimated by domain is not easy compared to the polesthe [A str¨ omto etCT al. formulating ˚zeros All the aforementioned NNM methods are confined to DT formulating matrix completion problem [Liu et al. 2013]. period for stiffas systems and b) Translating the zeros to CT domain is not as easy compared to the poles [ A str¨ o m et al. All the aforementioned NNM methods are[Liu confined to DT a matrix completion problem et al. due 2013]. domain Astr¨ om et al. formulating 1984]. is not as easy compared to the poles [˚ models and are not easily extended to the CT case to All the aforementioned NNM methods are confined to DT ˚str¨ 1984]. domain is not as easy compared to the poles [A om et al. All models and are not easily extended to the case due to the aforementioned NNM methods are CT confined to DT 1984]. the difference between DT and CT Hankel matrices. models and are not easily extended to the CT case due to The 1984].major challenge in direct CT identification is the the difference between DTextended and CT Hankel matrices. models and are not easily to the CT case due The major challenge in direct CT identification is the difference between DT and CT Hankel matrices. to The majorofchallenge in direct identification is the the estimation time derivatives. ToCT compute the derivatives, main of work is the between DT andcurrent CT Hankel estimation of time derivatives. To compute the The major in direct CT identification is the The Thedifference main contribution contribution of the the current workmatrices. is to to develop develop estimation ofchallenge time derivatives. To compute the derivatives, derivatives, linear filters, integral methods and modulating functions a new algorithm that incorporates a nuclear norm minThe main contribution of the current work is to develop linear filters, integral methods and modulating functions estimation of time derivatives. compute thesuch derivatives, new algorithm that incorporates a nuclear norm minmain contribution of the current work is to develop linear filters, integral methods and modulating functions are used [Garnier et al. 2003]. To Linear filters as Gen- aThe a new algorithm that incorporates a nuclear norm minare used [Garnier et al. 2003]. Linear filters such as Genlinear filters, integral methods and modulating functions a new algorithm that incorporates a nuclear norm minare used [Garnier et al. 2003]. Linear filters such as Genare used [Garnier et al. 2003]. Linear filters such as Gen2405-8963 © 2018, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved.

Copyright © 2018 IFAC 562 Peer review© under responsibility of International Federation of Automatic Copyright 2018 IFAC 562 Control. Copyright © 2018 IFAC 562 10.1016/j.ifacol.2018.05.089 Copyright © 2018 IFAC 562

5th International Conference on Advances in Control and Optimization of Dynamical Systems Santhosh Kumar Varanasi et al. / IFAC PapersOnLine 51-1 (2018) 530–535 February 18-22, 2018. Hyderabad, India

imization method into the SID methods along with the goodness of fit constraint for identification of CT statespace models. As tuning the filtering parameter is easy in the GPMF method compared to that of the Bi-linear transformation based techniques, the former is used to estimate the derivative information. Hankel matrix of past inputs are considered as Instrumental variables to make the proposed algorithm consistent. The formulated NNM problem is then converted into an equivalent SDP and solved using CVX for MATLAB. Monte-Carlo simulations are performed to show that the proposed algorithm is unbiased. This paper is organized as follows. Section 2 gives the problem statement and explains subspace identification using GPMF method. In Section 3, the objective function in Nuclear norm minimization is defined and its equivalent SDP formulation is presented. Section 4 provides numerical results of the proposed method. Section 5 presents the comparison of proposed method with the N4SID method followed by conclusions in Section 6 2. CONTINUOUS TIME STATE-SPACE IDENTIFICATION WITH GPMF 2.1 Problem statement Consider a linear time invariant CT state-space system x(t) ˙ = Ax(t) + Bu(t) (1) y(t) = Cx(t) + Du(t) where u(t) ∈ Rm , y(t) ∈ Rl and x(t) ∈ Rn are input, output and state vector respectively. A ∈ Rn×n , B ∈ Rn×m , C ∈ Rl×n , D ∈ Rl×m are state, input, output and feed-forward matrices respectively. A noise corrupted measurement of y(t) is available as (2) z(tk ) = y(tk ) + (tk ), k = 1, 2, . . . , N. The objective of identification is to estimate A, B, C and D matrices from input-output data i.e., {(u(tk ), z(tk )); k = 1, 2, . . . , N } 2.2 Methodology

Yk,r

 y(tk+1 ) . . . y(tN )  y(t ˙ k+1 ) . . . y(t ˙ N)    = .. ..  ∈ R(r+1)l×(N −k) ..   . . . (r) (r) (r) y (tk ) y (tk+1 ) . . . y (tN ) 

y(tk ) y(t ˙ k) .. .

Xk = [xk xk+1 xk+2 . . . xN ] ∈ Rn×(N −k) Uk,r is identical to Yk,r in structure but with y(t) replaced by u(t) and  T Γ = C T (CA)T . . . (CAr )T ∈ Rl(r+1)×n   D 0 ... 0 D ... 0  CB ∈ Rl(r+1)×m(r+1) Hd =  . . . . ..   .. .. . . CAr−1 B CAr−2 B . . . D As the measured output data is corrupted by noise, time derivatives are not estimated directly. In the Generalized Poisson Moment Functional (GPMF) approach, time derivatives of the data are shifted onto some smooth filters by the process of convolution. The rth order GPMF transformation of a signal v(t) over [0, tk ] is given by the following convolution [Saha and Rao 1983] vfp (tk )

y(t) = Cx(t) + Du(t) y(t) ˙ = CAx(t) + CBu(t) + Du(t) ˙ 2 y¨(t) = CA x(t) + CABu(t) + CB u(t) ˙ + D¨ u(t) .. .

(3)

y (r) (t) = CAr x(t) + CAr−1 Bu(t) + · · · + Du(r) (t) For the multiple sampled time instants, the set of equations (Eq. 3) can be represented in matrix form, Yk,r = ΓXk + Hd Uk,r

(4)

where 563

=



tk 0

v(τ )gp (tk − τ )dτ = [gp  v](tk )

(5)

where gp (t) denotes the impulse response of the filter  r+1 β sp s+α , α ∈ R+ , β ∈ R+ are the Poisson filter constant and gain respectively (in most cases β = α is considered) and p = 0, 1, . . . , r. Thus the output equation is represented as Zfk,r = ΓXfk + Hd Ufk,r + Efk,r

Zfk,r

Consider Eq. 1 and assume that the input and output data is r times differentiable. Repeated substitution of x(t) ˙ into the output relation results into the following set of equations,

531



zf0 (tk ) zf0 (tk+1 )  zf1 (tk ) zf1 (tk+1 )  =  .. ..  . . (r) (r) zf (tk ) zf (tk+1 )

 . . . zf0 (tN ) . . . zf1 (tN )   ..  ..  . . (r) . . . zf (tN )

(6)

(7)

  Xfk = xfk xfk+1 xfk+2 . . . xfN

In Subspace Identification, the main idea is to use instrumental variables (IV) and orthogonal projection to remove the effect of future inputs and the noise from Eq. 6. A matrix of the past inputs as defined in Eq. 8 is considered as Instrumental variables [Merc`ere et al. 2007] to eliminate the effect of noise. 

u(t1 ) u(t2 )  u(t2 ) u(t3 ) Wp =  ..  ... . u(tk−1 ) u(tk ) where, q = N − k and r has to that r  n.

 ... u(tq ) . . . u(tq+1 )   (8) .. ..  . . . . . u(tq+k−1 ) be selected in such a way

5th International Conference on Advances in Control and 532 Optimization of Dynamical Systems Santhosh Kumar Varanasi et al. / IFAC PapersOnLine 51-1 (2018) 530–535 February 18-22, 2018. Hyderabad, India

It was proved in [Merc`ere et al. 2007] that Wp satisfies the following conditions, 1 f E WpT = 0 lim N →∞ N k,r and   1 f T Xk Wp = n. rank lim N →∞ N Therefore, in the asymptotic case, the output equation reduces to (9) Zfk,r WpT = ΓXfk WpT + Hd Ufk,r WpT ˆ = Uf W T and define projection onto the orthogLet U k,r p ˆ as Π ˆ ⊥ = I − onal complement of the row space of U U T ˆ ˆ T −1 ˆ T ˆ U (UU ) U . Applying this projection on Eq. 9, we obtain f Zfk,r Wp ΠU (10) ˆ ⊥ = ΓXk Wp ΠU ˆ⊥ It can be observed that the matrix G = Zfk,r Wp ΠU ˆ ⊥ is in the column space of matrix the Γ. The column spaces of G and Γ are equal when these two matrices have the same rank. Under mild conditions such as the system being ˆU ˆ T ) having full rank, it can observable and the matrix (U be shown that rank(G) = rank(Γ) = n (order of the statespace model) [Verhaegen and Verdult 2007]. An important step in subspace identification is the singular value decomposition (SVD) of matrix G, which is used to estimate the range of observability matrix (Γ). Once this range is estimated, Multi-variable Output Error StateSpace (MOESP) method [Verhaegen and Verdult 2007] is used to estimate the system matrices A, B, C, D and the initial states in DT case. Inspired from the same approach, in the current paper, a CT version of the MOESP method is proposed. The main steps in the CT MOESP (CTMOESP) method are given below.

is no such guarantee. The reliance on SVD for low rank approximation also makes it difficult to extend the existing subspace methods to the problem of missing input or output measurements, to incorporate prior knowledge (such as bounds on outputs, stability or dissipativity constraints) or to add regularization terms on model inputs and outputs. Being a convex relaxation of the rank minimization problem [Fazel et al. 2001], Nuclear norm minimization (NNM) provides an interesting alternative to tackle the aforementioned challenges. For a given linear map L(ˆ z ) and a goodness of fit an NNM problem can be formulated as a single multi-criteria optimization problem, min L(ˆ z )∗ subject to z − zˆ <  (11) ˆ Z

where, ||.||∗ denotes the nuclear norm. There exist two main methods to solve the above NNM problem 1) formulating it as an trace minimization [Fazel et al. 2001] and 2) using the Alternating Direction Method of Multipliers (ADMM) [Fazel et al. 2013]. In the current work, the formulation of trace minimization from NNM is considered. An equivalent formulation of the NNM problem as trace minimization is given below [Fazel et al. 2001]. 1 (trace(V1 ) + trace(V2 )) 2   z) V1 L(ˆ ≥0 subject to L(ˆ z ) V2   I z(tk ) − zˆ(tk ) ≥0 (z(tk ) − zˆ(tk ))T I minimize

(12)

Trace minimization can be also be solved directly using CVX package in MATLAB, 1 min ||z − zˆ||22 + λ||L(ˆ z )||∗ (13) zˆ 2 where λ is the regularization parameter.

(1) Compute [Un , Σn , Tn ] = svd(G) and let Cˆ be the first 1/2 l rows of κ = Un Σn . (2) for some s > n, κ(1 : (s − 1)l, :)Aˆ = κ(l + 1 : sl, :) is ˆ ˆ f Wp Π ˆ ⊥ used to obtain A. For the Nuclear Norm CT Identification, G = Z k,p U f ˆ ˆ (3) B, D are estimated along with the initial state x ˆ0 by is a linear map of zˆ and the problem is formulated as solving a least squares problem 1  k−1 minimize (trace(V1 ) + trace(V2 ))   2    ˆ D, ˆ x zf − Cˆ Aˆk x ˆfk − ufτ,r (B, ˆfk ) = argmin ˆ ˆ f Wp Π ˆ ⊥ V1 Z  k,r k,p U τ =0 ≥0 subject to (14) 2 ˆ f Wp Π ˆ ⊥    V2 Z k,p U  f k−τ −1   vec(B) − uk,r ⊗ Il vec(D) ⊗Cˆ Aˆ 2 I z(tk ) − zˆ(tk ) ≥0 (z(tk ) − zˆ(tk ))T I where, ˆ zfk,r represents the k th column of the matrix ˆ f and similar definition holds true for uf . vec(.) or Z k,r k,r represents the matrix formed by placing the columns of the matrix one below the other, Il denotes Identity 1 ˆ f Wp Π ˆ ⊥ ||∗ matrix of size l × l and ⊗ represents Kronecker min ||z − zˆ||22 + λ||Z (15) k,p U zˆ 2 product. In the current work, a range of λ values are considered 3. NUCLEAR NORM MINIMIZATION OF CT [Verhaegen and Hansson 2016] and the identification alSTATE-SPACE MODELS gorithm is performed at each λ and the λ at which the value of the objective is minimum is considered. Once zˆ The use of instrumental variables guarantees that the esti- is estimated from the algorithm, CT-MOESP method is mates are consistent i.e., range(Γ) is estimated correctly in used to identify the system matrices A, B, C and D and the limiting case (N → ∞). However for finite data, there the state sequence. 564

5th International Conference on Advances in Control and Optimization of Dynamical Systems Santhosh Kumar Varanasi et al. / IFAC PapersOnLine 51-1 (2018) 530–535 February 18-22, 2018. Hyderabad, India

3.1 Numerical Methodology For the purpose of numerical studies, a zero-order hold assumption is considered on the output. Hence Eq. 5 is written as

f = Sk,r



tk

s(τ )gr (tk −τ )dτ =

0

k  i=1

s(ti )gr (tk −ti )∆ti (16)

Using this, the matrix in Eq. 7 is simplified as Zfk,r = IfN,r Hk,N

(17)

where,

IfM,r



g0 (tN )∆tN g1 (tN )∆tN = ..  . gr (tN )∆tN

. . . g0 (tk )∆tk . . . g1 (tk )∆tk .. .. . . . . . gr (tk )∆tk

 . . . g0 (t1 )∆t1 . . . g1 (t1 )∆t1   .. ..  . . . . . gr (t1 )∆t1

533

systems were chosen since it is difficult to identify them with existing CT subspace techniques when the model order information is not available [Dai et al. 2016]. Even when the model order is known, the existing DT subspace methods and Indirect methods fail to identify the models accurately for System-1 [Merc`ere et al. 2007]. A pseudo random binary signal (PRBS) is chosen as input. The output signals from the systems are corrupted by additive white noise with signal to noise ratio (SNR, defined as 10 × log10 (Psignal /Pnoise ) where, Psignal is the power of the signal and Pnoise is the power of the noise; power of signal is defined as the root mean square value of given noiseless signal) of 20 dB. The input and the output signals are sampled at a constant rate of 0.05 seconds. The initial state vector in both the systems is considered as zero. 4.1 System-1 A 2-input and 1-output MISO model illustrated in [Merc`ere et al. 2007, Dai et al. 2016] is considered        x˙ 1 (t) 1 1 0 1 0 x1 (t) x˙ 2 (t) = −3 −2 −1 x2 (t) + 2 1 u(t) x˙ 3 (t) 1 2 −1 −2 −1 x3 (t)   x1 (t) y(t) = [1 0 0] x2 (t) x3 (t)

∈ Rr×(N −1)

 ... 0 0 z(t1 ) ... 0 z(t1 ) z(t2 )  ... z(t1 ) z(t2 ) z(t3 )    .. .. .. ..  . . . . Hk,M  . . . z(tN −k−1 ) z(tN −k ) z(tN −k+1 )   Monte-Carlo simulations are performed to evaluate the .. .. .. ..  performance of the proposed schemes. The magnitude . . . . . . . z(tN −2 ) z(tN −1 ) z(tN ) of the frequency response along with the output of the (N −1)×(N −k) estimated systems subjected to MCS are shown in Figs. 1 ∈R and 2 respectively. Thus the steps in overall algorithm can be summarized as follows Bode diagram 0 0 0  0  0 0   . ..  . . = . z(t ) z(t ) 2  1  . .  .. .. z(tk ) z(tk+1 )

(1) Divide the given data into two parts one for identification (identification data, z) and other for validation (validation data, zv ). (2) Select the range of λ as Λ = [λmin , λmax ]. (3) For each λi ∈ Λ, solve Eq. 14 or Eq. 15 to estimate zˆ using identification data. (4) Compute SVD of G and model order (ˆ n) is selected based on the singular values of G. This can be done manually or automatically [Ljung 2007]. (5) From SVD of G and the selected model order (ˆ n), Cˆ and Aˆ are estimated from step-1 and step-2 of MOESP method. ˆ D ˆ and state sequence are estimated from step-3 of (6) B, MOESP method. ˆ B, ˆ C, ˆ D} ˆ and (7) Use the estimated system matrices {A, validation data to compute the simulated output (ˆ zv ). N1 The cost function is evaluated as J(λi ) = k=1 ||zv − zˆv ||22 , where N1 is length of validation data. (8) The selected output of algorithm is the model that gives the minimal function value (J(λi ), λi ∈ Λ). 4. NUMERICAL EXPERIMENTS To illustrate the performance of the proposed method, two third order multi-variable systems are considered. These 565

from: input−1

from: input−2

10

10

0

Magnitude (dB)



0

−10 −10 −20 −20 −30 −30 −40

−40

−50

−60 −2 10

0

10

2

10

Estimated Actual Mean of Estimated

−50 −2 10

0

10

2

10

Frequency (rad/sec)

Fig. 1. Actual and estimated magnitude of frequency response of System-1 From the Figs 1 and 2, it can be concluded that the proposed method is able to identify the system accurately in the high frequency region and the mean of the estimated model is close to that of the actual system which shows that the estimates are likely unbiased.

5th International Conference on Advances in Control and 534 Optimization of Dynamical Systems Santhosh Kumar Varanasi et al. / IFAC PapersOnLine 51-1 (2018) 530–535 February 18-22, 2018. Hyderabad, India 30

1.5 Estimated Actual Mean of Estimated

Estimated Actual Mean of estimated

20

Output value

Output value

1

0.5

10

0

-10

0

-20 −0.5

-30 0

2

4

6

8

10

12

14

16

18

20

Time (in seconds) −1

0

2

4

6

8

10

12

14

16

18

20

Time (in seconds)

Fig. 2. Comparison of actual and estimated response of System-1 subject to PRBS input. 4.2 System-2 A 3-input and 1-output MISO model illustrated in [Dai et al. 2016] is considered as the other system for illustration. 





x˙ 1 (t) 0 0 10 x˙ 2 (t) = 0 −0.1 −10 x˙ 3 (t) −1 1 0   x1 (t) y(t) = [0 2 0] x2 (t) x3 (t)







 x1 (t) 111 x2 (t) + 2 1 3 u(t) x3 (t) 120

Fig. 4. Comparison of actual and estimated response of System-2 subject to PRBS input. the time but in a few cases, the difference between the actual and the estimated output responses is considerable. However, the mean of the output is close to the actual response suggesting that the model estimates are unbiased. 5. COMPARISON The proposed method is compared with the N4SID algorithm (Van Overschee and De Moor [1994]) using MATLAB System Identification toolbox. The poles of the identified systems with the proposed method and N4SID are considered for the purpose of comparison and the results are shown in Figs. 5 and 6 respectively.

Monte-Carlo simulations are performed to evaluate the performance of the proposed schemes and the results of the magnitude of frequnecy response and the output plots are shown in Figs. 3 and 4 respectively.

5

4

Estimated Actual

imaginary part of pole

3

2

1

0

−1

−2

−3

−4

−5 −0.06

−0.055

−0.05

−0.045

−0.04

−0.035

−0.03

−0.025

−0.02

real part of pole

Fig. 5. Actual and Estimated poles of System-2 with the proposed method

Fig. 3. Actual and estimated magnitude of frequency response of System-2 From Figs. 3 and 4, it can be concluded that the proposed method is able to identify the system accurately most of 566

From Fig. 5, it can be observed that the proposed method is able to provide good estimates of the poles with the real pole having a comparatively higher variance. As the estimated values are centered around the actual poles, the results seem to be unbiased. From Fig. 6, it can be concluded that with N4SID, the identification results are not accurate and biased.

5th International Conference on Advances in Control and Optimization of Dynamical SystemsSanthosh Kumar Varanasi et al. / IFAC PapersOnLine 51-1 (2018) 530–535 February 18-22, 2018. Hyderabad, India 6 Estimated Actual

imaginary part of pole

4

2

0

−2

−4

−6

−0.6

−0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

real part of pole

Fig. 6. Actual and Estimated poles of System-2 with N4SID 6. CONCLUSIONS In this paper, identification of CT state-space models is considered. A new algorithm which incorporates the goodness of fit into the traditional GPMF based subspace methods using Nuclear norm minimization (NNM) is proposed. The NNM is formulated as a SDP and solved using CVX for MATLAB. Two Numerical case studies are considered to illustrate the performance of the proposed method. In the current paper only measurement noise is considered. Extending the same problem to the case of state noise will be pursued in future. To reduce the computational complexity of the NNM problem, instead of SDP formulation, the ADMM algorithm can be used. REFERENCES ˚ Astr¨ om, K.J., Hagander, P., and Sternby, J. (1984). Zeros of sampled systems. Automatica, 20, 31–38. Dai, M., He, Y., and Yang, X. (2016). Continuous-time system identification with nuclear norm minimization and gpmf-based subspace method. IEEE/CAA Journal of Automatica Sinica, 3, 184–191. Fazel, M., Hindi, H., and Boyd, S.P. (2001). A rank minimization heuristic with application to minimum order system approximation. In Proceedings of the American Control Conference, 4734–4739. Fazel, M., Pong, T.K., Sun, D., and Tseng, P. (2013). Hankel matrix rank minimization with applications to system identification and realization. SIAM Journal on Matrix Analysis and Applications, 34, 946–977. Garnier, H., Mensler, M., and Richard, A. (2003). Continuous-time model identification from sampled data: implementation issues and performance evaluation. International journal of Control, 76, 1337–1357. Haverkamp, B., Chou, C., Verhaegen, M., and Johansson, R. (1996). Identification of continuous-time mimo state space models from sampled data, in the presence of process and measurement noise. In Proceedings of the 35th IEEE Conference on Decision and Control, 1539– 1544. Liu, Z., Hansson, A., and Vandenberghe, L. (2013). Nuclear norm system identification with missing inputs and outputs. Systems & Control Letters, 62, 605–612. 567

535

Ljung, L. (1999). System identification: Theory for the user. Prentice Hall, New Jersey. Ljung, L. (2007). System identification toolbox for use with matlab. version 7. Merc`ere, G., Ouvrard, R., Gilson, M., and Garnier, H. (2007). Subspace based methods for continuous-time model identification of mimo systems from filtered sampled data. In Proceedings of European Control Conference (ECC), 4057–4064. Ohta, Y. (2011). Stochastic system transformation using generalized orthonormal basis functions with applications to continuous-time system identification. Automatica, 47, 1001–1006. Rao, G. and Unbehauen, H. (2006). Identification of continuous-time systems. IEE Proceedings-Control theory and applications, 153, 185–220. Saha, D.C. and Rao, G.P. (1983). Identification of continuous dynamical systems: the Poisson moment functional (PMF) approach. Springer. Unbehauen, H. and Rao, G. (1990). Continuous-time approaches to system identificationa survey. Automatica, 26, 23–35. Van Overschee, P. and De Moor, B. (1994). N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica, 30, 75– 93. Verhaegen, M. and Hansson, A. (2016). N2SID: Nuclear norm subspace identification of innovation models. Automatica, 72, 57–63. Verhaegen, M. and Verdult, V. (2007). Filtering and system identification: a least squares approach. Cambridge university press.