Model Order Reduction of Stochastic Linear Systems by Moment Matching*

Proceedings of the 20th World The International Federation of Congress Automatic Control Proceedings of 20th Congress Proceedings of the the 20th World World The International Federation of Congress Automatic Control Toulouse, France, July 9-14, 2017 The International Federation of Automatic Control Proceedings of the 20th World Congress The International Federation of Automatic Control Toulouse, France, July 9-14, 2017 Available online at www.sciencedirect.com Toulouse, France,Federation July 9-14, 9-14, 2017 2017 The International of Automatic Control Toulouse, France, July Toulouse, France, July 9-14, 2017

ScienceDirect

IFAC PapersOnLine 50-1 (2017) 6332–6337 of Stochastic Model Order Reduction Model Order Reduction of Stochastic Model Order Reduction of Stochastic  Model Order Reduction of Linear Systems by Moment Matching ModelSystems Order Reduction of Stochastic Stochastic Linear by Moment Matching Linear Systems Systems by by Moment Moment Matching Matching  Linear LinearGiordano Systems by Moment Matching ∗ ∗∗ Scarciotti ∗ Andrew R. Teel ∗∗

Giordano Scarciotti Andrew R. Teel

∗ ∗∗ Giordano Giordano Scarciotti Scarciotti ∗∗ Andrew Andrew R. R. Teel Teel ∗∗ ∗∗ ∗ Giordano Scarciotti Andrew R. Teel ∗ Department of Electrical and Electronic Engineering, Imperial College Department of Electrical and Electronic Engineering, Imperial ∗ ∗London, Department of Electrical Electrical andUK Electronic Engineering, Imperial College College London SW7 2AZ, (e-mail: Engineering, [email protected]) Department of and Electronic Imperial College ∗London, London SW7 2AZ, UK (e-mail: [email protected]) ∗∗ Department ofand Electrical andUK Electronic Engineering, College London, London SW7 2AZ, (e-mail: [email protected]) Computer Engineering Department,Imperial University of ∗∗ Electrical London, London SW7 2AZ, UK (e-mail: [email protected]) Electrical and Computer Engineering Department, University of ∗∗ ∗∗ Electrical London, London SW7 2AZ, UK (e-mail: [email protected]) and Computer Engineering Department, University of California, Santa Barbara, CA 93106-9560 USA (e-mail: and Santa Computer Engineering Department, ∗∗ Electrical California, Barbara, CA 93106-9560 USA University (e-mail: of

Electrical and Santa Computer Engineering Department, California, Santa Barbara, CA 93106-9560 93106-9560 USA University (e-mail: of [email protected]) California, Barbara, CA USA (e-mail: [email protected]) California, Santa Barbara, CA 93106-9560 USA (e-mail: [email protected]) [email protected]) [email protected]) Abstract: In this paper we characterize the moments of stochastic linear systems by means of Abstract: In this paper we characterize the moments of stochastic linear systems by means of Abstract: this paper characterize the moments of linear by means the solutionIn a stochastic equation generalizes the classical Sylvester Abstract: Inof this paper we wematrix characterize the which moments of stochastic stochastic linear systems systems by equation. means of of the solution of a stochastic matrix equation which generalizes the classical Sylvester equation. Abstract: paper wematrix characterize the which moments of stochastic linear systems means of the solution a stochastic equation generalizes the classical Sylvester equation. The solutionInof ofthis the matrix equation is used to define the steady-state response ofby the system the solution of a stochastic matrix equation which generalizes the classical Sylvester equation. The solution of the matrix equation is used to define the steady-state response of the system the solution of a stochastic matrix equation which generalizes the classical Sylvester equation. The solution of the matrix equation is used to define the steady-state response of the system which is then exploited to define a family of stochastic reduced order models. In addition, the The solution of the matrix equation is used to define the steady-state response of the system which is then exploited to aa family of stochastic reduced order In the The of the matrix equation is used define the steady-state response the system whichsolution is of then exploited to define define family ofinto stochastic reduced order models. models. In ofaddition, addition, the notions stochastic reduced order model the mean and stochastic reduced order model which is then exploited to define a family of stochastic reduced order models. In addition, the notions of stochastic reduced order model in the mean and stochastic reduced order model which then exploited to define a family stochastic reduced order models. In addition, the notions of stochastic reduced order model in the mean and reduced order model in the isvariance are introduced. While theofdetermination of astochastic reduced order model based on notions of stochastic reduced order model in the mean and stochastic reduced order model in the variance are introduced. While the determination of aastochastic reduced order model based on notions of stochastic reduced order model in the mean and reduced order model in the variance are introduced. While the determination of reduced order model based on the stochastic notion of moment has high computational complexity, stochastic reduced order in the variancenotion are introduced. While the computational determination complexity, of a reduced order model based on the stochastic of moment has high stochastic reduced order in the variance are introduced. While the determination of a reduced order model based on the stochastic notion of moment has high computational complexity, stochastic reduced order models in the mean and variance can be determined more easily, yet they preserve some of the the stochastic notionand of moment has high computational complexity, stochastic reduced order models in the mean variance can be determined more easily, yet they preserve some of the the stochastic notionand of the moment has computational complexity, stochastic order models in the the mean and variance can be determined more easily, yet yet they preserve somefamilies of the stochastic properties of system to high be determined reduced. The differences between these reduced three models in mean variance can be more easily, they preserve some of the stochastic properties of the systemcan to be reduced. The differences between these three models in are the mean and variance be more easily, yet they preserve somefamilies of the stochastic properties of to be reduced. The differences between these families of models illustrated by system means of simulations. stochastic properties of the the system to numerical be determined reduced. The differences between these three three families of models are illustrated by means of numerical simulations. stochastic properties of the system to be reduced. The differences between these three families of models are illustrated by means of numerical simulations. of models are illustrated by means of numerical simulations. © 2017, IFAC (International Federation of numerical Automatic Control) Hosting by Elsevier Ltd. All rights reserved. of models are illustrated by means of Keywords: Model reduction; stochastic systems; simulations. stochastic modeling; steady-state. Keywords: Model reduction; stochastic systems; stochastic modeling; steady-state. Keywords: Keywords: Model Model reduction; reduction; stochastic stochastic systems; systems; stochastic stochastic modeling; modeling; steady-state. steady-state. Keywords: Model reduction; stochastic systems; stochastic steady-state. 1. INTRODUCTION active areamodeling; of research spawning new results on a multi1. INTRODUCTION active area of research spawning new results a 1. activeofarea area of research research spawning new results on onstochastic a multimultitude various applications suchnew as quantum 1. INTRODUCTION INTRODUCTION active of spawning results on a multitude of various applications such as quantum stochastic 1. also INTRODUCTION area research new results and onstochastic aNurdin multitude of various applications such as quantum systems, seeofe.g. Nurdinspawning (2014); Techakesari Stochastic systems, known as stochastic differential active tude of various applications such as quantum stochastic see e.g. Nurdin (2014); Techakesari and Nurdin Stochastic systems, systems, also also known known as as stochastic stochastic differential differential systems, tude of various applications such as quantum stochastic systems, see e.g. Nurdin (2014); Techakesari and Nurdin (2016), and system biology, see e.g. Bruna et al. (2014); Stochastic equations, are systems of differential equations in which see system e.g. Nurdin (2014); Techakesari and Nurdin Stochastic are systems, alsoofknown as stochastic (2016), and biology, see e.g. Bruna et al. (2014); equations, systems differential equationsdifferential in which which systems, seeAnderson e.g. Nurdin (2014); Techakesari and Nurdin (2016), and system biology, see e.g. Bruna et al. (2014); Stochastic systems, known as stochastic differential Sootla and (2014); Johnson et al. (2015). equations, are systems of differential equations in one (or more) of thealso terms is a stochastic process, i.e. a systems, (2016), and system biology, see e.g. Bruna et al. (2014); equations, are systems of differential equations in which Sootla and Anderson (2014); Johnson et al. (2015). one (or more) of the terms is a stochastic process, i.e. a (2016), and system biology, see e.g. Bruna et al. (2014); Sootla and Anderson (2014); Johnson et al. (2015). equations, are systems of differential equations in which In this paper we propose three families of stochastic one (or more) of the terms is a stochastic process, i.e. a time sequence representing the evolution of some variable Sootla and Anderson (2014); Johnson et al. (2015). one (or more) representing of the terms the is aevolution stochastic process, i.e. a In this paper we propose three families of stochastic time sequence of some variable Sootla and Anderson (2014); Johnson et al. (2015). In this paper we propose three families of stochastic one (or more) of the terms is a stochastic process, i.e. a reduced order models based on the moment matching time representing evolution of which is subject to random the variations. Stochastic systems In this paper we propose families of stochastic time sequence sequence representing evolution of some some variable variable order models basedthree on the moment matching which is subject subject to random random the variations. Stochastic systems reduced this paper we propose three families of stochastic reduced order models on the moment matching time sequence representing the evolution ofsystem some variable method. We define the based moments of linear stochastic syswhich is to variations. Stochastic systems are extensively studied in optimal control, biology In reduced order models based on the moment matching which is subject to random variations. Stochastic systems method. We define the moments of linear stochastic sysare extensively studied in optimal control, system biology reduced order models based on the moment matching method. We define the moments of linear stochastic syswhich is subject to random variations. Stochastic systems tems in terms of the solution of a stochastic Sylvester are extensively studied in optimal control, system biology and finance, see e.g. Karatzas and Shreve (1991); Yong and method. We define the moments of linear stochastic sysare extensively studied in optimal control, system biology tems in terms of the solution of aa linear stochastic Sylvester and finance, see e.g. Karatzas and Shreve (1991); Yong and method. We define the moments of stochastic systems in terms of the solution of stochastic Sylvester are extensively studied in optimal control, system biology equation. As a result the moments are processes. and finance, see e.g. Karatzas and Shreve (1991); Yong and Zhou (1999); Munk (2011). Some applications of stochastic tems in terms of the solution of a stochastic Sylvester and finance, e.g. (2011). Karatzas andapplications Shreve (1991); Yong and equation. As a result the moments are stochastic processes. Zhou (1999);see Munk Some of stochastic stochastic tems in terms of the solution of a Sylvester equation. As a result the moments are stochastic processes. and finance, see e.g. Karatzas and Shreve (1991); Yong and Exploiting this stochastic notion of moment we provide Zhou (1999); Munk (2011). Some applications of differential include theapplications modeling of stochastic uncertain equation. Asthis a result the moments are stochasticwe processes. Zhou (1999);equations Munk (2011). Some stochastic notion of moment provide differential equations include theapplications modeling of of stochastic uncertain Exploiting Asthis a of result the moments are stochastic processes. stochastic notion of moment provide Zhou (1999); Munk (2011). Some aExploiting first family stochastic models which havewe the same differential include the modeling uncertain systems, theequations filtering problem (which has been solved by equation. Exploiting this stochastic notion of moment we provide differential equations include the modeling of uncertain aa first family of stochastic models which have the same systems, the filtering problem (which has been solved by Exploiting this stochastic notion of moment we provide first family of stochastic models which have the same differential equations include the modeling of uncertain steady-state output response as the original system for systems, the filtering problem (which has been solved by the Kalman-Bucy filter, see Kalman (1960); Kalman and a first family of stochastic models which have the systems, the filtering problem (which(1960); has been solvedand by steady-state output response as the original systemsame for the Kalman-Bucy filter, see Kalman Kalman a first family of stochastic models which have the same steady-state output response as the original system for systems, the filtering problem (which has been solved by some selectedoutput input response signals. Since the resulting model the Kalman-Bucy filter, see Kalman (1960); Kalman and Bucy (1961)), the Dirichlet problem, the optimal stopping steady-state as the original system for the Kalman-Bucy filter, see Kalman (1960); Kalman and some selected input signals. Since the resulting model Bucy (1961)), the Dirichlet problem, the optimal stopping steady-state output response as the original system for some selected input signals. Since the resulting model the Kalman-Bucy filter, see Kalman (1960); Kalman and is computationally expensive, we propose two additional Bucy (1961)), the Dirichlet problem, the optimal stopping problem, the production planning problem, the stopping portfolio some selected input signals. Since the resulting model Bucy (1961)), the Dirichlet problem, the optimal computationally expensive, we propose two additional problem, the production production planning problem, the stopping portfolio is some selected input signals. Since the resulting model is computationally expensive, we propose two additional Bucy (1961)), the Dirichlet problem, the optimal families of stochastic reduced order models. The first famproblem, the planning problem, the portfolio management problem (all these problems are analyzed computationally expensive, we propose two additional problem, the production problem, are the analyzed portfolio is families of reduced order models. The first management problem (all (allplanning these problems problems computationally expensive, we propose additional families of stochastic stochastic reduced order models. The firstasfamfamproblem, the production problem, theand portfolio ily possesses the same steady-state outputtwo mean the management problem these are analyzed in the monographs Øksendal (2013) and Yong Zhou is families of stochastic reduced order models. The first fammanagement problem (allplanning these problems are analyzed ily possesses the same steady-state output mean as the in the monographs Øksendal (2013) and Yong and Zhou families ofbe stochastic reduced order models. The first family possesses the same steady-state output mean as the management problem (all these problems are and analyzed system to reduced, whereas the second family possesses in the monographs Øksendal (2013) and Yong Zhou (1999)). ily possesses the same steady-state output mean as the in the monographs Øksendal (2013) and Yong and Zhou system to be reduced, whereas the second family possesses (1999)). ily possesses the same steady-state output mean as the system to be reduced, whereas the second family possesses in the monographs Øksendal (2013) and Yong and Zhou the same steady-state output mean and variance. (1999)). The complexity of the algorithms solving these problems system to be reduced, whereas the second family possesses (1999)). same steady-state output mean and variance. The complexity of of the the algorithms algorithms solving solving these these problems problems the system to be reduced, whereas the second family possesses the same steady-state output mean and variance. (1999)). The rest of the paper is organized as follows. In Section 2 The complexity usually grows more than linearly with the dimension of same steady-state output meanasand variance. The complexity of the algorithms problems The rest of the paper is organized follows. In Section 2 usually grows more more than linearly solving with the thethese dimension of the the same steady-state output mean and variance. The rest of the paper is organized as follows. In Section 2 The complexity of the algorithms solving these problems we recall some properties of stochastic systems and model usually grows than linearly with dimension of the model, see Lindquist and Picci (2015). This motivates The rest of the paper is organized as follows. In Section usually grows more thanand linearly with theThis dimension of we recall some properties of stochastic systems and model2 the model, see Lindquist Picci (2015). motivates The rest of the properties paper is organized as follows. Inand Section 2 we recall some of systems model usually grows more than linearly with theThis dimension of reduction which are useful tostochastic streamline the presentation the model, Lindquist and Picci (2015). motivates need of see stochastic model reduction techniques that are we recall some properties of stochastic systems and model the model, see Lindquist and Picci (2015). This motivates which are useful to streamline the presentation the need need of see stochastic model reduction techniques that are are reduction we recall some properties of stochastic systems and model reduction which are useful to streamline the presentation model, Lindquist and Picci (2015). This motivates of the results of the paper. In Section 3 we provide a the of stochastic model reduction techniques that capable of reducing the order of the system maintaining at whichofare useful to In streamline the presentation the needofofreducing stochastic techniques that are of the results the paper. Section 33 we provide a capable themodel orderreduction of the the system system maintaining at reduction reduction which are useful to streamline the presentation of the results of the paper. In Section we provide a need of stochastic model reduction techniques that are stochastic notion of moment for stochastic linear systems. capable of reducing the order of maintaining at the same time some of the properties of the original systhe results of of the paper. for In stochastic Section 3linear we provide a capable oftime reducing the order of the system maintaining at of stochastic notion moment systems. the same some of the properties of the original sysof the results of the paper. In Section 3 we provide stochastic notion of moment for stochastic linear systems. capable of reducing the order of the system maintaining at In Section 4 we propose three families of stochastic rethe same time some of the properties of the original system. The literature of model reduction for linear stochastic stochastic notion of moment for stochastic linear systems. the same time someofofmodel the properties of the original sys- In Section 4 we propose three families of stochastic re-a tem. The literature reduction for linear stochastic stochastic notion of moment for stochastic linear systems. In Section 44 models we families stochastic rethe time theofproperties of the original sys- duced order that three have different of and tem. The literature of reduction for linear systems rich,some especially methods based on stochastic In Section we propose propose three families ofcomplexity stochastic retem.same Theis ofofmodel model for linear stochastic duced order models that have different complexity and systems isliterature rich, especially especially of reduction methods based based on stochastic stochastic In Section 4 we propose three families of stochastic reduced order models that have different complexity and tem. The literature of model reduction for linear properties. In Section 5 we compare these three families systems is rich, of methods on balancing and Hankel norm approximation, see e.g. Desai duced order models that have different complexity and systems is rich, especially of methods based on stochastic In Section 55 we compare these three families balancing and Hankel norm approximation, see e.g. Desai Desai properties. duced order models that have different complexity and properties. In Section we compare these three families systems is rich, especially of methods based on stochastic of models by means of numerical simulations. Finally, balancing and Hankel norm approximation, see e.g. and Pal (1982, 1984); Desai et al. (1985); Green (1988a,b); properties. In Section 5 we comparesimulations. these three families balancing and Hankel norm et approximation, see (1988a,b); e.g. Desai of by means of numerical Finally, and Pal (1982, (1982, 1984); Desai Desai al. (1985); (1985); Green Green In Section 5 concluding we compareremark. these three families of models models by means of numerical simulations. Finally, balancing and Hankel norm approximation, see (1988a,b); e.g. (1990, Desai properties. Section 6 contains someof and Pal 1984); et al. Harshavardhana et al.Desai (1984); Wang andGreen Safonov of models by means numerical simulations. Finally, and Pal (1982, 1984); et al. (1985); (1988a,b); Section 6 contains some concluding remark. Harshavardhana et al. (1984); Wang and Safonov (1990, of models by means ofconcluding numerical simulations. Finally, Section 66 contains some remark. and PalLindquist (1982, 1984); Desai et al. (1985); Green (1988a,b); Harshavardhana et (1984); Wang Safonov (1990, 1991); andal. Picci (1996); Xuand and Chen (2003). Section contains some concluding remark. Notation. We use standard notation. C denotes the Harshavardhana et al. (1984); Wang and Safonov (1990, <0 1991); Lindquist et andal.Picci Picci (1996); Xuand andSafonov Chen (2003). (2003). Notation. We use standard notation. C denotes the Section 6 contains some concluding remark. <0 Harshavardhana (1984); Wang (1990, 1991); Lindquist and (1996); Xu and Chen Model reduction of stochastic systems continues to be an Notation. We use standard notation. C denotes the set of complex numbers with negative real part; C de1991); Lindquist and Picci (1996); Xu and Chen (2003). <0 0 Notation. We use standard notation. C denotes <0 part; C0 the Model reduction of stochastic systems continues to(2003). be an an set of complex numbers with negative real de1991); Lindquist and Picci (1996); Xu and Chen Model reduction of stochastic systems continues to be Notation. We use standard notation. C denotes the set of complex numbers with negative real part; C denotes the set of complex numbers with zero real part. <0 Model reduction of stochastic systems continues to be an 0  This work set of complex numbers with negative real part; C de0part. was supported in part by Imperial College London under notes the set of complex numbers with zero real  Model reduction of stochastic to be an set of complex numbers with negative real part; C deThis work was supported in part bysystems Imperial continues College London under notes the set of complex numbers with zero real part. The symbol ι denotes the imaginary unit, I denotes the 0  notes the set of complex numbers with zero real part. theThis Junior Research Fellowship Scheme, by NSF grant no. ECCS work was supported in part by Imperial College London under The symbol ιι denotes the imaginary unit, II denotes the workResearch was supported in part by Imperial College London under theThis Junior Fellowship Scheme, by NSF grant no. ECCSnotes the set of complex numbers with zero real part. The symbol denotes the imaginary unit, denotes the  identity matrix and σ(A) denotes the spectrum of 1508757 andResearch by AFOSR grant no. AFOSR The symbol ι denotes the imaginary unit, I denotes the theThis Junior Research Fellowship Scheme, byFA9550-15-1-0155. NSF grant no. ECCSECCSwork was supported in part by Imperial College London under identity matrix and denotes spectrum of the the Junior Fellowship Scheme, NSF grant no. 1508757 and by AFOSR grant no. AFOSRbyFA9550-15-1-0155. The symbol ι denotes the imaginary unit, I denotes identity matrix and σ(A) σ(A) denotes the the spectrum of the 1508757 and by grant AFOSR the Junior Fellowship Scheme, NSF grant no. ECCSidentity matrix and σ(A) denotes the spectrum of the 1508757 andResearch by AFOSR AFOSR grant no. no. AFOSRbyFA9550-15-1-0155. FA9550-15-1-0155. identity matrix and σ(A) denotes the spectrum of the 1508757 and by AFOSR grant no. AFOSR FA9550-15-1-0155.

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

Proceedings of the 20th IFAC World Congress Giordano Scarciotti et al. / IFAC PapersOnLine 50-1 (2017) 6332–6337 Toulouse, France, July 9-14, 2017

matrix A ∈ Rn×n . The symbol ⊗ indicates the Kronecker product and ||A|| indicates the induced Euclidean matrix norm. The superscript  denotes the transposition operator. (Σ, A, P) indicates a probability space with a given set Σ, a σ-algebra A on Σ and a probability measure P on the measurable space (Σ, A). 2. PRELIMINARIES In this section we recall some definitions and notions which are instrumental for the remainder of the paper. We provide, in order, some properties of stochastic systems, then of moment matching and finally of the Kronecker product. Note that all the stochastic integrals in this paper are intended as Itˆ o integrals. Let Wt a standard Wiener process defined on a probability space (Σ, A, P). A stochastic process xt is a function of two variables such that for each t ∈ R, x(t, ·) is a random variable and for each w ∈ Ω, x(·, w) is called path of x. For ease of notation, we indicate the path of xt as x(t). Consider a linear, single-input, single-output, continuoustime, stochastic system described by the equations dx = (Ax + Bu)dt + (F x + Gu)dWt , n

y = Cx, n×n

(1) n×1

,B∈R , with x(t) ∈ R , y(t) ∈ R, u(t) ∈ R, A ∈ R F ∈ Rn×n , G ∈ Rn×1 , C ∈ R1×n and assume that AF = F A. The initial value problem associated to (1) has a unique solution described by  t Φ(t − τ )(B − F G)u(τ )dτ x(t) = Φ(t)x(0) + 0 t (2) Φ(t − τ )Gu(τ )dWτ , + 0

where, since AF = F A,

1 2 Φ(t) = e(A− 2 F )t+F Wt ,

(3)

see e.g. Gard (1988). We can  compute the expectation, or mean, m(t) = E[x(t)] = Σ xdP of the process xt . For system (1), the mean obeys the deterministic differential equation m ˙ = Am + Bu. (4) Similarly, the variance M (t) = E[x(t)x(t) ] is governed by the equation M˙ = AM + M A + F M F  + Bum + m(Bu) (5) +F m(Gu) + Gu(F m) + Gu(Gu) ,

6333

with ω(t) ∈ Rν , S ∈ Rν×ν , L ∈ R1×ν and L and ω(0) such that the triple (S, ω(0), L) is minimal 2 . It is useful to recall the notion of moment for deterministic systems. Definition 1. Let F = G = 0. The moments of system (1) at (S, L) are the elements of the matrix CΠ, where Π is the unique solution of the Sylvester equation AΠ + BL = ΠS. (7) Note that usually moments are defined in terms of the transfer function of a linear system. However, it has been shown that the moments can be defined equivalently as in Definition 1, or by means of a specific steady-state output response, see Astolfi (2010) and Gallivan et al. (2004, 2006). In fact, note that if σ(A) ⊂ C<0 (and F = G = 0), the interconnection of system (1) and (6) possesses a global invariant manifold described by M = {(x, ω) ∈ Rn+ν : x = Πω}. Hence, for all t ∈ R, y(t) = CΠω(t) + CeAt (x(0) − Πω(0)), (8) where CΠω(t) is the steady-state output response of system (1) driven by (6). These relations suggest that we may be able to define the moments of the stochastic system (1) in terms of its steady-state output response, without the need of defining the transfer function. The last notion that is useful to recall is the problem of determining the nearest Kronecker product approximation. Given a matrix Q ∈ RN ×N , the problem consists in determining the two matrices T1 ∈ Rn×n and T2 ∈ Rm×m , with N = nm, such that ||Q − T1 ⊗ T2 || is minimized. The solution of this problem is given in Genton (2007). After rearranging the elements of Q in a new matrix 2 2 called R(Q) ∈ Rn ×m (see Genton (2007) for details), we compute the singular value decomposition U  R(Q)V = diag(δ1 , . . . , δq ) of R(Q).  The solution to the  problem is given by vec(T1 ) = δ1 u1 and vec(T2 ) = δ1 v1 , where u1 and v1 are the first columns of U and V , respectively, and q = rank(R(Q)). 3. DEFINITION OF THE STOCHASTIC MOMENTS In this section we provide a stochastic notion of moment based on the steady-state response of system (1). Utilizing classical deterministic arguments, we expect that the steady-state of system (1) driven by (6) can be characterized as  t  t x(t) = Φ(t − τ )(B − F G)u(τ )dτ + Φ(t − τ )Gu(τ )dWτ −∞

= X ω(t),

−∞

see Gard (1988). Let A = I ⊗ A + A ⊗ I + F ⊗ F . If (1) is said asymptotically meanσ (A) ⊂ C<0 , then system  square stable. If σ A − 12 F 2 ⊂ C<0 , then system (1) is said asymptotically stable almost surely. For linear systems, asymptotic mean-square stability implies almost surely asymptotic stability. The converse is not true, see Damm (2004). S ∈ Rν×ν with Consider now a non-derogatory 1 matrix η characteristic polynomial p(s) = i=1 (s − si )ki , where η ν = i=1 ki , and the deterministic signal generator

(9) with  t  t X = Φ(t−τ )(B−F G)LeS(τ −t)dτ + Φ(t−τ )GLeS(τ −t)dWτ .

1

2 See Padoan et al. (2017) for detail on the significance of the minimality assumption.

ω˙ = Sω,

u = Lω,

(6)

A matrix is non-derogatory if its characteristic and minimal polynomials coincide.

−∞

−∞

(10) In the next lemma we show that if the function Φ(t) goes to zero asymptotically almost surely, ω(t) is bounded, and Wt is Wiener process, then the steady-state response of system (1) driven by the signal generator (6) is in fact characterized by (9) and (10). In the previous section we have recalled that for deterministic systems the steadystate can be described by means of the solution of a

6522

Proceedings of the 20th IFAC World Congress 6334 Giordano Scarciotti et al. / IFAC PapersOnLine 50-1 (2017) 6332–6337 Toulouse, France, July 9-14, 2017

Sylvester equation. We now show that a similar result holds for stochastic systems as well. Lemma 1. Consider the interconnection of system (1) and   the signal generator (6). Assume σ A − 12 F 2 ⊂ C<0 and σ(S) ⊂ C0 . Then the steady-state response of the output of such interconnection is y(t) = CX ω(t), almost surely, where X given in (10) solves the stochastic matrix equation (11) dX = (AX − X S + BL)dt + (F X + GL)dWt . We can now define the moments in the stochastic framework. Definition 2. Consider system (1) and the signal generator (6). The moments of system (1) at (S, L) are the elements of the matrix CX , where X is given in (10). Remark 1. The moments CX are functions of time and of the stochastic process Wt . Moments as functions of time have been introduced in Scarciotti and Astolfi (2016a) for discontinuous signal generators, in Scarciotti and Astolfi (2016b) for some classes of hybrid systems and in Scarciotti et al. (2017) for linear differential inclusions. Hence, the moments CX generalize both the classical timeinvariant moments (see e.g. Astolfi (2010); Scarciotti and Astolfi (2016c)) and the time-varying moments introduced in those papers. 4. REDUCED ORDER MODELS In this section we provide several families of reduced order models. We begin giving reduced order models which are based on the stochastic notion of moment. Then we show that these models are unsatisfactory from a computational point of view. As a consequence, we introduce other families of reduced order models, namely models in the mean and in the variance, that are computationally tractable. We begin providing the definition of model of system (1) at (S, L). Definition 3. Consider system (1) and the signal generator (6). The system described by the equations dξ = (Ar ξ + Br u)dt + (Fr ξ + Gr u)dWt , ψ = Cr ξ, (12) where ξ(t) ∈ Rν , ψ(t) ∈ R, Ar ∈ Rν×ν , Br ∈ Rν×1 , Fr ∈ Rν×ν , Gr ∈ Rν×1 , Cr ∈ R1×ν , is a stochastic model of system (1) at (S, L), if system (12) has the same moments of system (1) at (S, L). From this definition a result follows straightforwardly. Lemma 2. Consider the interconnection of system (1) and   the signal generator (6). Assume σ A − 21 F 2 ⊂ C<0 and σ(S) ⊂ C0 . Then system (12) is a stochastic model of system (1) at (S, L) if the steady-state solution P of the equation dP = (Ar P − PS + Br L)dt + (Fr P + Gr L)dWt , (13) is well-defined and such that (14) CX = Cr P, where X is the unique solution of (11). The family of models obtained from Lemma 2 has several free parameters that can be used to satisfy the two matrix equations (13) and (14). Hence, we provide the following result.

Proposition 1. Consider the interconnection of system (1)   and the signal generator (6). Assume σ A − 12 F 2 ⊂ C<0 and σ(S) ⊂ C0 . Then the system dξ = ((S − Br L)ξ + Br u)dt + (−Gr Lξ + Gr u)dWt , ψ = CX ξ,

(15) is a stochastic model L) for any Br   of system 1(1) at (S, and Gr such that σ S − Br L − 2 (Gr L)2 ⊂ C<0 .

The main drawback of the results based on the stochastic notion of moment is that they require solving equation (11) to determine X . Equation (11) is a stochastic matrix equation which consists of nν linear stochastic equations. Hence, while the models proposed in Proposition 1 possess the same moments or, equivalently, the same steady-state output response of system (1), they can be hardly considered reduced order models. Since this issue is intrinsic for all the models based on the process X , we look for an alternative way to characterize reduced order models which is computationally advantageous. To address this problem we introduce the definition of model in the mean. Definition 4. Consider system (1) and the signal generator (6). The system described by equation (12) is a stochastic model in the mean of system (1) at (S, L), if the mean of system (12) has the same moments at (S, L) of the mean of system (1). System (12) is a stochastic reduced order model in the mean of system (1) at (S, L) if ν < n. A family of reduced order models in the mean follows. Proposition 2. Consider the interconnection of system (1) and the signal generator (6). Assume σ (A) ⊂ C<0 and σ(S) ⊂ C0 . Then the system dξ = ((S − Br L)ξ + Br u)dt + (Fr ξ + Gr u)dWt , (16) ψ = CΠξ, is a stochastic model in the mean of system (1) at (S, L), for any Fr and Gr , and for any Br such that σ (S − Br L) ⊂ C<0 . Remark 2. The advantage of reduced order models in the mean is that we only need to determine CΠ, which can be computed with a plethora of efficient methods. For instance, suggestions on how to compute CΠ are given in Antoulas (2005) and Astolfi (2010). Alternatively the method proposed in Scarciotti and Astolfi (2017) allows computing the moments CΠ in less than O(ν 3 ). The main drawback of the models in the mean is that they do not preserve, in a systematic way, information regarding the matrices F and G of the system to be reduced. To preserve also this information we propose reduced order models in the variance. Definition 5. Consider system (1) and the signal generator (6). The system described by equation (12) is a stochastic model in the variance of system (1) at (S, L), if both the mean and the variance of system (12) have the same moments at (S, L) as the mean and the variance of system (1). System (12) is a stochastic reduced order model in the variance of system (1) at (S, L) if ν < n. Note that we have defined models in the variance in such a way that they preserve both mean and variance. To the end of determining reduced order models in the variance

6523

Proceedings of the 20th IFAC World Congress Giordano Scarciotti et al. / IFAC PapersOnLine 50-1 (2017) 6332–6337 Toulouse, France, July 9-14, 2017

we need a preliminary result, namely a description of the steady-state of equation (5). Lemma 3. Consider the interconnection of system (5) and the signal generator (6). Recall that A = I ⊗ A + A ⊗ I + F ⊗ F and assume σ (A) ⊂ C<0 and σ(S) ⊂ C0 . Then the steady-state response of such interconnection is M (t) = Kω(t)ω(t) K  , for all t ∈ R, where K is such that K = K ⊗ K is the unique solution of the augmented Sylvester equation AK + B = KS, (17) with S = I ⊗ S + S ⊗ I and B = BL ⊗ Π + Π ⊗ BL + GL ⊗ F Π + F Π ⊗ GL + GL ⊗ GL. We are now ready to give a family of reduced order models in the variance. Proposition 3. Consider the interconnection of system (1) and the signal generator (6). Assume σ (A) ⊂ C<0 and σ(S) ⊂ C0 . Let P be any invertible matrix such that CKP = CΠ, where Π is the unique solution of (7) and K is the unique solution of (17). Define Ar = P SP −1 − Br LP −1 for any Br such that σ(Ar ) ⊂ C<0 . Let Fr and Gr be such that Fr ⊗Fr + Gr L⊗Fr P + Fr P ⊗Gr L + Gr L⊗Gr L =

15 10 5 0 −5

5. EXAMPLE In this section we use a numerical example to show the differences between the three families of models that we have proposed. Most of the quantities in this example have been randomly generated in MATLAB. The command

0

1

2

3

4

5

6

7

8

9

10

6

7

8

9

10

t 10 8 6 4 2 0

0

1

2

3

4

5

t

Fig. 1. Top graph: time history of the output of system (1) (solid/blue) and time history of the output of the stochastic reduced order model (15) (dotted/red). The time is in seconds. Bottom graph: time history of the corresponding absolute error. 15 10

= I ⊗Ar + Ar ⊗I − I ⊗S − S ⊗I + Br L⊗P + P ⊗Br L, (18) and that σ (I ⊗ Ar + Ar ⊗ I + Fr ⊗ Fr ) ⊂ C<0 . Then the system dξ = (Ar ξ + Br u)dt + (Fr ξ + Gr u)dWt , (19) ψ = CKξ,

is a stochastic model in the variance of system (1) at (S, L). Remark 3. The main difficulty in the determination of the family of models (19) is to solve equation (18). However, since Gr is a free parameter, we can select it equal to zero. As a result, equation (18) becomes Fr ⊗Fr = I ⊗Ar +Ar ⊗I −I ⊗S−S⊗I +Br L⊗P +P ⊗Br L, from which we can determine Fr as the nearest Kronecker approximation of Fr ⊗ Fr . Another convenient solution is given by selecting Fr = Gr L and solving (18) with respect to Gr . Finally, note that Fr = −Gr L is not a solution. Remark 4. The determination of a reduced order model in the variance requires two nearest Kronecker approximations: the first is needed to compute K from K and the second to compute Fr from Fr . Note that the determination of the nearest Kronecker approximation introduces an error called separability approximation error, see Genton (2007). This error is zero if the rearranged matrix R(K) (or R(Fr )) has only one singular value different from zero. While the only way to influence the separability approximation error of R(K) is by selecting other matrices L and S, the separability approximation error of R(Fr ) can be also modified using the matrix Br .

6335

5 0 −5

0

5

10

15

10

15

t

10 8 6 4 2 0

0

5

t

Fig. 2. Top graph: time history of the output of system (1) (solid/blue), time history of the output of the first reduced order model in the mean (16) (dashed/black) and time history of the output of the second reduced order model in the mean (16) (dotted/red). The time is in seconds. Bottom graph: time history of the corresponding absolute errors. rng(’default’) has been used to set the random generator to the Mersenne Twister with seed zero. Thus, the numerical values and the figures can be reproduced. We have generated matrices A, B and C of order n = 50 with the function rss. The matrix A has been redefined subtracting the identity to itself to make its eigenvalues strictly negative. In addition, F = 0.05 · A, G = B, ν = 2, L = [1 0], all the initial states are defined as vectors of 1’s of appropriate dimension and the matrix S is taken with eigenvalues ±ι. Thus, the reduced order model has order 2. The system and the reduced order models are defined in MATLAB as sde objects and simulated with the method symByEuler. We have determined a stochastic reduced order model

6524

Proceedings of the 20th IFAC World Congress 6336 Giordano Scarciotti et al. / IFAC PapersOnLine 50-1 (2017) 6332–6337 Toulouse, France, July 9-14, 2017

with Fr = −Gr L can provide an approximation which is comparable with the more computationally expensive models in the variance.

15 10 5

6. CONCLUSION

0 −5

0

5

t

10

15

10

15

10 8 6 4 2 0

0

5

t

Fig. 3. Top graph: time history of the output of system (1) (solid/blue) and time history of the output of the reduced order model in the variance (19) (dotted/red). The time is in seconds. Bottom graph: time history of the corresponding absolute error. assigning the eigenvalues of S − Br L using Br . The two eigenvalues are selected as two of the eigenvalues of A. The matrix Gr is selected as Gr = 0.2·Br . Fig. 1 shows the time history of the output of system (1) (solid/blue line) and the time history of the output of the stochastic reduced order model (15) (dotted/red line). The bottom graph displays the corresponding absolute error. As expected the stochastic reduced order model has the same steady-state output response of system (1). We have determined two reduced order models in the mean. The first model maintains the same parameters of the stochastic reduced order model. Since Fr = −Gr L, the resulting model is deterministic. The second model instead is determined with Fr = −0.5 · Ar and, as a result, the model is stochastic. Fig. 2 shows the time history of the output of system (1) (solid/blue line), the time history of the output of the first reduced order model in the mean (16) (dashed/black line) and the time history of the output of the second reduced order model in the mean (16) (dotted/red line). The bottom graph displays the corresponding absolute errors. We note that the first model is a coarse approximation which does not maintain the stochastic nature of the original system. On the other hand, the second model shows features similar to the original model and the corresponding absolute error is smaller. Finally, we have determined a reduced order model in the variance. The eigenvalues of Ar are set to −1 and −2 using Br , whereas Gr is selected equal to zero. Fig. 3 shows the time history of the output of system (1) (solid/blue line) and the time history of the output of the reduced order model in the variance (19) (dotted/red line). The bottom graph displays the corresponding absolute error. We note that the reduced order model shows features similar to the original model. By direct comparison, this last absolute error is marginally smaller than the error associated to the second model in the mean. In summary, the example shows that with a “good” selection of the parameters, reduced order models in the mean

In this paper, exploiting the solution of a stochastic Sylvester equation, we have provided a stochastic characterization of the moments of stochastic linear systems. The moments, which are stochastic processes, are used to define a family of reduced order models which have the same steady-state output response for some specific classes of inputs. We have also introduced the notions of stochastic reduced order model in the mean and stochastic reduced order model in the variance. Further research includes the extension of the present framework to nonlinear stochastic systems and stochastic hybrid systems. REFERENCES Antoulas, A. (2005). Approximation of Large-Scale Dynamical Systems. SIAM Advances in Design and Control, Philadelphia, PA. Astolfi, A. (2010). Model reduction by moment matching for linear and nonlinear systems. IEEE Transactions on Automatic Control, 55(10), 2321–2336. Bruna, M., Chapman, S.J., and Smith, M.J. (2014). Model reduction for slow-fast stochastic systems with metastable behaviour. The Journal of Chemical Physics, 140(17), 174107. Damm, T. (2004). Rational Matrix Equations in Stochastic Control. Lecture Notes in Control and Information Sciences. Springer Berlin Heidelberg. Desai, U.B. and Pal, D. (1982). A realization approach to stochastic model reduction and balanced stochastic realizations. In Proceedings of the 21st IEEE Conference on Decision and Control, 1105–1112. Desai, U.B. and Pal, D. (1984). A transformation approach to stochastic model reduction. IEEE Transactions on Automatic Control, 29(12), 1097–1100. Desai, U.B., Pal, D., and Kirkpatrick, R.D. (1985). A realization approach to stochastic model reduction. International Journal of Control, 42(4), 821–838. Gallivan, K., Vandendorpe, A., and Van Dooren, P. (2004). Sylvester equations and projection-based model reduction. Journal of Computational and Applied Mathematics, 162(1), 213–229. Gallivan, K.A., Vandendorpe, A., and Van Dooren, P. (2006). Model reduction and the solution of Sylvester equations. In 17th International Symposium on Mathematical Theory of Networks and Systems, Kyoto, Japan. Gard, T.C. (1988). Introduction to Stochastic Differential Equations. Monographs and textbooks in pure and applied mathematics. M. Dekker. Genton, M.G. (2007). Separable approximations of spacetime covariance matrices. Environmetrics, 18(7), 681– 695. Green, M. (1988a). Balanced stochastic realizations. Linear Algebra and its Applications, 98, 211–247. Green, M. (1988b). A relative error bound for balanced stochastic truncation. IEEE Transactions on Automatic Control, 33(10), 961–965. Harshavardhana, P., Jonckheere, E., and Silverman, L. (1984). Stochastic balancing and approximation-

6525

Proceedings of the 20th IFAC World Congress Giordano Scarciotti et al. / IFAC PapersOnLine 50-1 (2017) 6332–6337 Toulouse, France, July 9-14, 2017

stability and minimality. IEEE Transactions on Automatic Control, 29(8), 744–746. Isidori, A. (1995). Nonlinear Control Systems. Communications and Control Engineering. Springer, Third edition. Johnson, T., Bartol, T., Sejnowski, T., and Mjolsness, E. (2015). Model reduction for stochastic CaMKII reaction kinetics in synapses by graph-constrained correlation dynamics. Physical Biology, 12(4), 045005. Kalman, R.E. (1960). A new approach to linear filtering and prediction problems. ASME Journal of Basic Engineering, 82(1), 35–45. Kalman, R.E. and Bucy, R.S. (1961). New results in linear filtering and prediction theory. ASME Journal of Basic Engineering, 95–108. Karatzas, I. and Shreve, S.E. (1991). Brownian Motion and Stochastic Calculus. Graduate Texts in Mathematics. Springer New York. Lindquist, A. and Picci, G. (1996). Canonical correlation analysis, approximate covariance extension, and identification of stationary time series. Automatica, 32(5), 709–733. Lindquist, A. and Picci, G. (2015). Linear Stochastic Systems: A Geometric Approach to Modeling, Estimation and Identification. Series in Contemporary Mathematics. Springer Berlin Heidelberg. Munk, C. (2011). Fixed Income Modelling. OUP Oxford. Nurdin, H.I. (2014). Structures and transformations for model reduction of linear quantum stochastic systems. IEEE Transactions on Automatic Control, 59(9), 2413– 2425. Øksendal, B. (2013). Stochastic Differential Equations: An Introduction with Applications. Universitext. Springer, sixth edition. Padoan, A., Scarciotti, G., and Astolfi, A. (2017). A geometric characterisation of the persistence of excitation condition for the solutions of autonomous systems. To appear on IEEE Transactions on Automatic Control. Pavlov, A., van de Wouw, N., and Nijmeijer, H. (2006). Uniform Output Regulation of Nonlinear Systems: A Convergent Dynamics Approach. Systems & Control: Foundations & Applications. Birkh¨ auser Boston. Scarciotti, G. and Astolfi, A. (2016a). Model reduction by matching the steady-state response of explicit signal generators. IEEE Transactions on Automatic Control, 61(7), 1995–2000. Scarciotti, G. and Astolfi, A. (2016b). Model reduction for hybrid systems with state-dependent jumps. In IFAC Symposium Nonlinear Control Systems, Monterey, CA, USA, 862–867. Scarciotti, G. and Astolfi, A. (2016c). Model reduction of neutral linear and nonlinear time-invariant time-delay systems with discrete and distributed delays. IEEE Transactions on Automatic Control, 61(6), 1438–1451. Scarciotti, G. and Astolfi, A. (2017). Data-driven model reduction by moment matching for linear and nonlinear systems. Automatica, 79, 340–351. Scarciotti, G., Teel, A.R., and Astolfi, A. (2017). Model reduction for linear differential inclusions: robustness and time-variance. In Proceedings of the 2017 American Control Conference, Seattle (to appear). Sootla, A. and Anderson, J. (2014). On projectionbased model reduction of biochemical networks part II:

6337

The stochastic case. In Proceedings of the 53rd IEEE Conference on Decision and Control, 3621–3626. Techakesari, O. and Nurdin, H.I. (2016). Tangential interpolatory projection for model reduction of linear quantum stochastic systems. To appear on IEEE Transactions on Automatic Control. Wang, W. and Safonov, M.G. (1990). A tighter relativeerror bound for balanced stochastic truncation. Systems & Control Letters, 14(4), 307–317. Wang, W. and Safonov, M.G. (1991). Relative-error bound for discrete balanced stochastic truncation. International Journal of Control, 54(3), 593–612. Xu, S. and Chen, T. (2003). H∞ model reduction in the stochastic framework. SIAM Journal on Control and Optimization, 42(4), 1293–1309. Yong, J. and Zhou, X.Y. (1999). Stochastic Controls: Hamiltonian Systems and HJB Equations. Stochastic Modelling and Applied Probability. Springer New York.

6526