A Structure-preserving Model Reduction Algorithm for Dynamical Systems with Nonlinear Frequency Dependence

6th IFAC Symposium on System Structure and Control 6th IFAC Symposium on System Structure and Control 6th Symposium on June 22-24, 2016. Istanbul, TurkeyStructure 6th IFAC IFAC Symposium on System System Structure and and Control Control June 22-24, 2016. Istanbul, TurkeyStructure 6th IFAC Symposium on System and Control June 22-24, 2016. Istanbul, Turkey Available online at www.sciencedirect.com June 22-24, 2016. Istanbul, Turkey June 22-24, 2016. Istanbul, Turkey

ScienceDirect IFAC-PapersOnLine 49-9 (2016) 056–061

A Structure-preserving Model Reduction A Structure-preserving Model Reduction A Structure-preserving Model Reduction AAlgorithm Structure-preserving Model Reduction for Dynamical Systems with Algorithm for Dynamical Systems with Algorithm for Dynamical Systems with Algorithm for Dynamical Systems with Nonlinear Frequency Dependence Nonlinear Frequency Dependence Nonlinear Frequency Dependence Nonlinear Frequency Dependence

∗ ∗ ∗ Klajdi Sinani ∗ Serkan Gugercin ∗ Christopher Beattie ∗ Klajdi Sinani Gugercin Beattie ∗ ∗ ∗ ∗ Serkan ∗ Christopher ∗ Klajdi Sinani Serkan Gugercin Christopher Beattie Klajdi Sinani ∗ Serkan Gugercin ∗ Christopher Beattie ∗ Klajdi Sinani Serkan Gugercin Christopher Beattie ∗ ∗ Department of Mathematics, Virginia Tech, Blacksburg, VA, ∗ Department of Mathematics, Virginia Tech, Blacksburg, VA, ∗ Department of of Mathematics, Virginia Tech, Blacksburg, VA, USA (e-mail: {klajdi,gugercin,beattie}@vt.edu) Department Mathematics, Virginia ∗ 24061-0123, Department of Mathematics, Virginia Tech, Tech, Blacksburg, Blacksburg, VA, VA, 24061-0123, USA (e-mail: {klajdi,gugercin,beattie}@vt.edu) 24061-0123, USA (e-mail: {klajdi,gugercin,beattie}@vt.edu) {klajdi,gugercin,beattie}@vt.edu) 24061-0123, USA (e-mail: 24061-0123, USA (e-mail: {klajdi,gugercin,beattie}@vt.edu) Abstract: Very large-scale dynamical systems, even linear time-invariant systems, can present Abstract: Very Very large-scale large-scale dynamical dynamical systems, systems, even even linear linear time-invariant time-invariant systems, systems, can can present present Abstract: significant computational difficulties when used in numerical simulation. Model reduction is Abstract: Very large-scale dynamical systems, even linear time-invariant systems, can present Abstract: Very large-scale dynamical systems, even linear time-invariant systems, can present significant computational difficulties when used in numerical simulation. Model reduction is significant computational difficulties when used in numerical simulation. Model reduction is one response to this challenge but standard methods often are restricted to systems that significant computational difficulties when used in numerical simulation. Model reduction is significant computational difficulties when used in numerical simulation. Model reduction is one response to this challenge but standard methods often are restricted to systems that one response to this challenge but standard methods often are restricted to systems that are presented as standard first-order realizations; in the frequency domain such systems will one response this but methods are to one response to to this challenge challenge but standard standard methods often are restricted restricted to systems systems that are presented presented as standard standard first-order realizations; in the theoften frequency domain such such systems that will are as first-order realizations; in frequency domain systems will be linear in the frequency parameter. We consider here dynamical systems with a nonlinear are presented as standard first-order realizations; in the frequency domain such systems will are presented as standard first-order realizations; in the frequency domain such systems will be linear in the frequency parameter. We consider here dynamical systems with a nonlinear be linear in the frequency parameter. We consider here dynamical systems with a nonlinear frequency dependence; systems for which either aa standard first-order realization is unknown or be linear the parameter. We consider here systems with aa nonlinear be linear in in the frequency frequency parameter. We consider here dynamical dynamical systems with nonlinear frequency dependence; systems for which which either first-order realization is unknown or frequency dependence; systems for either aa standard standard first-order realization is unknown or inconvenient to obtain. Such systems may nonetheless have realizations that reflect important frequency dependence; systems for which either standard first-order realization is unknown or frequency dependence; systems for which either a standard first-order realization is unknown or inconvenient to obtain. Such systems may nonetheless have realizations that reflect important inconvenient to obtain. Such systems may nonetheless have realizations that reflect important structural features of the model and we may wish to retain this structure in any reduced model inconvenient to obtain. Such systems may nonetheless have realizations that reflect important inconvenient to obtain. Such systems may nonetheless have realizations that reflect important structural features of the model and we may wish to retain this structure in any reduced model structural features of ofInthe the model and we may wish wish to retain retain this structure structure inreduction any reduced reduced model used as surrogate. this work, we present a structure-preserving model algorithm structural features model and may to any model structural features ofIn model and we may wish to retain this this structure inreduction any reduced model used as as aa a surrogate. surrogate. Inthe this work, wewe present structure-preserving modelin reduction algorithm used this work, we present aaa structure-preserving model algorithm for systems having quite general nonlinear frequency dependence. We take advantage of recent used as a surrogate. In this work, we present structure-preserving model reduction algorithm used as a surrogate. In this work, we present a structure-preserving model reduction algorithm for systems having quite general nonlinear frequency dependence. We take advantage of recent for systems having quite general nonlinear frequency dependence. We take advantage ofrequire recent algorithms that produce high quality rational interpolants to transfer functions that only for systems having quite general nonlinear frequency dependence. We take recent for systemsthat having quitehigh general nonlinear frequency dependence. Wefunctions take advantage advantage ofrequire recent algorithms that produce high quality rational interpolants to transfer transfer functions that only onlyof require algorithms produce quality rational interpolants to that transfer function evaluation, thus allowing for nonstandard realizations that are nonlinear in the algorithms that produce high quality rational interpolants to transfer functions that only require algorithms that produce high quality rational interpolants to transfer functions that only require transfer function evaluation, thus allowing for nonstandard realizations that are nonlinear in the the transfer function evaluation, thusour allowing for nonstandard nonstandard realizations that are are that nonlinear in frequency parameter. However, final reduced model will have aa structure reflects the transfer function evaluation, thus allowing for realizations that nonlinear in transfer function evaluation, thus allowing for nonstandard realizations that are nonlinear in frequency parameter. However, our final reduced model will have structure that reflects the frequency parameter. However, our final reduced model will have a structure that reflects the structure of the original system, and indeed, may not have a rational transfer function. We frequency parameter. However, final reduced will aa structure that reflects frequency parameter. However, our final reduced model will have have structure thatfunction. reflects the the structure of of the original original system,our and indeed, maymodel not have have rational transfer function. We structure the system, and indeed, may not aaatranscendental rational transfer We illustrate our approach on a benchmark problem that offers a transfer function. structure of the original system, and indeed, may not have rational transfer function. We structure of the original system, and indeed, may not have a rational transfer function. We illustrate our approach on a benchmark problem that offers a transcendental transfer function. illustrate illustrate our our approach approach on on a a benchmark benchmark problem problem that that offers offers aaa transcendental transcendental transfer transfer function. function. illustrate our approach on a benchmark problem that offers © 2016, IFAC (International Federation of Automatic Control) Hostingtranscendental by Elsevier Ltd. transfer All rightsfunction. reserved. Keywords: dynamical systems, model reduction, structure-preservation, interpolation, H 2 Keywords: dynamical systems, model reduction, structure-preservation, interpolation, H 2 Keywords: dynamical systems, model reduction, structure-preservation, interpolation, H 2 approximation, nonlinear frequency dependency, generalized realization, Loewner framework Keywords: dynamical systems, model reduction, structure-preservation, interpolation, H 2 Keywords: dynamical systems, model reduction, structure-preservation, interpolation, H approximation, nonlinear frequency dependency, generalized realization, Loewner framework 2 approximation, nonlinear frequency dependency, generalized realization, Loewner framework approximation, nonlinear frequency dependency, generalized realization, Loewner framework approximation, nonlinear frequency dependency, generalized realization, Loewner framework 1. INTRODUCTION tems whose behavior is described via aa transfer function 1. INTRODUCTION tems whose behavior is described via transfer function 1. INTRODUCTION tems whose behavior is described via aamulti-input/multitransfer function 1. INTRODUCTION H(s). The extension of our approach to tems whose behavior is described via transfer function 1. INTRODUCTION tems whose behavior is described via a transfer function H(s). The extension of our approach to multi-input/multiH(s). The extension of our approach to multi-input/multioutput systems may be developed in an entirely analogous H(s). The extension of our approach to multi-input/multiDirect numerical simulation of dynamical systems plays H(s). The extension of our approach to multi-input/multioutput systems may be developed in an entirely analogous Direct numerical simulation of dynamical systems plays output systems developed in an entirely analogous and straightforward way, but here only focus on output systems may may be be developed inwe an will entirely analogous Direct numerical simulation of systems plays aDirect crucial role (and at times may be the only option) in numerical simulation of dynamical dynamical systems plays output systems may be developed in an entirely analogous and straightforward way, but here we will only focus on numerical simulation of dynamical systems plays aaDirect crucial role (and at times may be the only option) in and straightforward way, but here we will only focus on the SISO case to keep the presentation both simple and and straightforward way, but here we will only focus on crucial role (and at times may be the only option) in studying a great variety of complex physical phenomena a crucial role (and at times may be the only option) in and straightforward way, but here we will only focus on the SISO case to keep the presentation both simple and astudying crucial role (and at times may be the only option) in a great variety of complex physical phenomena the SISO case to keep the presentation both simple and concise. We assume that the transfer function H(s) of the the SISO case to keep the presentation both simple and studying a great variety of complex physical phenomena with applications ranging signal propagation in the studying great variety variety offrom complex physical phenomena SISOWe case to keep the both H(s) simpleof and concise. assume that thepresentation transfer function the studying aa great of complex phenomena with applications ranging from signalphysical propagation in the the concise. assume that transfer function underlying aa the generalized concise. We Wesystem assumehas that the transfer realization function H(s) H(s) of of the the with ranging from propagation in nervous system, to heat dissipation complex microwith applications applications ranging from signal signal in propagation in the the concise. We assume that the transfer function H(s) of the underlying system has generalized realization with applications ranging from signal propagation in the nervous system, to heat dissipation in complex microunderlying system has aa generalized T −1realization underlying system has generalized realization nervous system, to heat dissipation in complex microT K(s)−1realization electronic devices, vibration suppression in large wind nervous system, system, toto heat dissipation in complex complex microB(s) (1) H(s) underlying system has= generalized nervous to heat dissipation in microelectronic devices, to vibration suppression in large wind (1) H(s) =aC(s) C(s) T −1 T K(s) −1 B(s) electronic devices, to vibration suppression in large wind H(s) C(s) T K(s) turbines, and to timely prediction of storm surges before electronic devices, to vibration suppression in large wind H(s) = C(s) K(s) B(s) (1) n = n−1 B(s) n×n (1) electronic devices, to vibration suppression in large wind K(s) (1) H(s) C(s) turbines, and to timely prediction of storm surges before n, = n , B(s) n×n are where C : C → C B : C → C K : C → C turbines, and to timely prediction of storm surges before n n n×n where C :: C → C :: C → C :: C → C an advancing hurricane. However, the ever-present need turbines, and to timely prediction of storm surges before n ,, B n ,, K n×n are where C C → C B C → C K C → C are turbines, and to timely prediction of storm surges before n n n×n an advancing hurricane. However, the ever-present need analytic in the right half plane, and K(s) has full rank where C : C → C , B : C → C , K : C → C are an advancing hurricane. However, the ever-present need where C in : Cthe →right C , B : C → Cand , KK(s) : C has → Cfull are half plane, for model resolution leads to the inclusion of everan greater advancing hurricane. However, the ever-present need analytic analytic in the right half plane, and K(s) has full rank rank an advancing hurricane. However, the ever-present need for greater model resolution leads to the inclusion of everthroughout right analytic in right half plane, for greater model resolution leads to the inclusion of everanalytic in the the right half half plane. plane, and and K(s) K(s) has has full full rank rank throughout right half plane. greater detail at the modeling stage. The resulting largefor greater model resolution leads to the inclusion of everthroughout the right half plane. for greater model resolution leads to the inclusion oflargeever- throughout the greater detail at the modeling stage. The resulting half greater detail the stage. The resulting largethe right right half plane. plane. The state space dimension of the underlying dynamical scale dynamical even those that are ‘simple’ linear greater detail at atsystems, the modeling modeling stage. The resulting large- throughout The state space dimension of the underlying dynamical greater detail at the modeling stage. The resulting largescale dynamical systems, even those that are ‘simple’ linear The state space dimension of the underlying dynamical scale dynamical systems, even those that are ‘simple’ linear system typically is n and we will refer to this as the order The state space dimension of the underlying dynamical time-invariant systems, can present significant computascale dynamical systems, even those that are ‘simple’ linear The state space dimension of the underlying dynamical system typically is n and we will refer to this as the order scale dynamical systems, even those that are ‘simple’ linear time-invariant systems, can present significant computasystem typically is n and we will refer to this as the order time-invariant systems, can present significant computaof the system, even though H(s) may have a significantly system typically is n and we will refer to this as the order tional difficulties when used in numerical simulation due time-invariant systems, can present significant computasystem typically is n and we will refer to this as the order of the system, even though H(s) may have a significantly time-invariant systems, can present significant computational difficulties when used in numerical simulation due of the system, even though H(s) may have a significantly tional difficulties when used in numerical simulation due larger (even infinite) number of poles. We are interested of the system, even though H(s) may have a significantly to their sheer size; hence there is a persistent need to tional difficulties when used in numerical simulation due of the system, even though H(s) may have a significantly larger (even infinite) number of poles. We are interested tional difficulties when used in numerical simulation due to their sheer size; hence there is a persistent need to (even infinite) number of poles. We are interested to size; hence there is persistent need to in settings n could reach hundreds of thousands larger (even infinite) number of We interested approximate complex systems with to their their sheer sheersuch size;large hence there dynamical is a persistent need to larger larger (evenwhere infinite) number of poles. poles. We are are interested in settings where n could reach hundreds of to their sheer size; hence there is aa persistent need to approximate such large complex dynamical systems with in more. settings where nthat could reach hundreds of thousands thousands approximate such large complex dynamical systems with or We note the state-space representation in in settings where n could reach hundreds of thousands smaller, yet high accuracy, approximations. This is the approximate such large complex dynamical systems with in settings where n could reach hundreds of thousands or more. We note that the state-space representation in approximate such large complex dynamical systems with smaller, yet high accuracy, approximations. This is the or more. We note that the state-space representation in smaller, yet high accuracy, approximations. This is the (1) allows a much richer set of dynamical systems than or more. We note that the state-space representation in goal of model reduction: one constructs simpler (reduced smaller, yet high accuracy, approximations. This is the or more. We note that the state-space representation in (1) allows a much richer set of dynamical systems than smaller, yet high accuracy, approximations. This is the goal of model reduction: one constructs simpler (reduced (1) allows a much richer set of dynamical systems than goal of model reduction: one constructs simpler (reduced those represented by standard first-order realizations, i.e., (1) allows a much richer set of dynamical systems than order) models, which are much easier and faster to simgoal of model reduction: one constructs simpler (reduced (1) allows a much richer set of dynamical systems than those represented by standard first-order realizations, i.e., goal of model reduction: one constructs simpler (reduced order) models, which are much easier and faster to simn×n −1standard first-order T represented realizations, i.e., order) models, which much easier to n×n H(s) = C(sE − A)by B where E, A ∈R ,, and B, C ∈ those represented by standard first-order realizations, i.e., ulate requiring far fewer computational resources) order)(hence models, which are are much easier and and faster faster to simsim- those −1standard T those represented first-order realizations, H(s) = C(sE − A R and C ∈ order) models, which are much easier and faster to simulate (hence requiring far fewer computational resources) −1 n n×n −1 B H(s) = C(sE − A) A)by B where where E, E, A∈ ∈ prominent Rn×n and B, B, CTTTi.e., ∈ ulate (hence requiring far fewer computational resources) n×n ,, and n are all constant One example R H(s) = C(sE − A) B where E, A ∈ R B, C ∈ −1quantities. while retaining characteristics close to the original system. ulate (hence requiring far fewer computational resources) H(s) = C(sE − A) B where E, A ∈ R , and B, C ∈ n are all constant quantities. One prominent example R ulate (hence requiring far fewer computational resources) while retaining characteristics close to the original system. n are all constant quantities. One prominent example R while retaining characteristics close to the original system. nthe of form in (1) are systems with internal delays such are all constant quantities. One prominent example R These simpler reduced models can then serve as efficient while retaining characteristics close to the original system. are all constant quantities. One prominent example R of the form in (1) are systems with internal delays such while retaining characteristics close to the original system. These simpler reduced models can then serve as efficient −τ s −1 the form in (1) are systems with internal delays such These reduced models can serve as as H(s) = C(sE A B where ττ is the of the form in (1) are systems with internal delays such surrogates for the original, replacing them as components These simpler simpler reduced models can then then serve as efficient efficient of −τ s )−1 0 − 1e of the form in (1) − areA systems with internal delays such as H(s) = C(sE − A A where is the These simpler reduced models can then serve as efficient surrogates for the original, replacing them as components −τ −1 0 − 1e −τ s s) −1 B as H(s) = C(sE − A − A e ) B where τ is the surrogates for the original, replacing them as components 0 1 internal delay. Note that this transfer cannot as H(s) = C(sE − A − A e where τ is the −τ s )−1 B function in larger systems; facilitating rapid development of consurrogates for the original, replacing them as components 0 1 as H(s) = C(sE − A − A e ) B where τ is the internal delay. Note that this transfer function cannot surrogates for the original, replacing them as components in larger systems; facilitating rapid development of con0 1 internal delay. in Note that this thisfirst-order transfer function function cannot in larger systems; facilitating rapid development of conbe represented a standard form with finite internal delay. Note that transfer cannot trollers for real time applications; and enabling optimal in larger systems; facilitating rapid development of coninternal delay. Note that this transfer function cannot be represented in a standard first-order form with finite in larger systems; facilitating rapid development of controllers for real time applications; and enabling optimal represented in a standard first-order form with finite trollers for time applications; and dimensional choices E, A and A §4, we consider be represented in standard form with finite system uncertainty analysis. trollers design for real realand time applications; and enabling enabling optimal optimal be 0 ,, first-order 1 .. In be represented in aa for standard form with finite dimensional choices for E, A and A §4, we consider trollers for real time applications; and enabling optimal system design and uncertainty analysis. 0 , first-order 1 . In dimensional choices for E, A and A In §4, we consider system design and uncertainty analysis. 0 1 two models in the form of (1). For example, the model in dimensional choices for E, A , and A . In §4, we consider system design and uncertainty analysis. 0 1 dimensional choices for E, A , and A . In §4, we consider two models in the form of (1). For example, the model system design and uncertainty analysis. 0 1 In this paper, we will focus on model reduction of stable s in two models models in the the form form of (1). (1). For example, the model in In this paper, we will focus on model reduction of stable §4.1 has aa transfer function of the form H(s) = C((e − two in of For example, the model s in In this paper, we will focus on model reduction of stable two models in the form of (1). For example, the model in §4.1 has transfer function of the form H(s) = C((e s single-input/single-output linear dynamical sysIn this this paper, paper, we we will will focus focus(SISO) on model model reduction of stable stable s − §4.1 has aa transfer function of the form H(s) = C((e In on reduction of single-input/single-output (SISO) linear dynamical sys§4.1 has transfer function of the form H(s) = C((e − s − single-input/single-output single-input/single-output (SISO) (SISO) linear linear dynamical dynamical syssys- §4.1 has a transfer function of the form H(s) = C((e − single-input/single-output (SISO) linear dynamical sys-

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

2016 IFAC SSSC June 22-24, 2016. Istanbul, Turkey

Klajdi Sinani et al. / IFAC-PapersOnLine 49-9 (2016) 056–061

1)A2 + s2 A1 + I)−1 B where A2 and A1 are constant matrices, and C and B are constant vectors.

57

more details on model reduction of standard first-order systems, we refer the reader to Antoulas (2005); Beattie and Gugercin (2015); Baur et al. (2014); Antoulas et al. (2010) and the references there in. However, for generalized realization H(s) = C(s)K(s)−1 B(s) with nonlinear frequency dependence throughout the state-space quantities, many of the generic approaches do not apply except for special cases such as for second-order structure H(s) = C(s2 M + sK + G)−1 B; see, e.g., Bai and Su (2005); Su and Craig Jr (1991); Chahlaoui et al. (2005); Meyer and Srinivasan (1996); Reis and Stykel (2008). See also recent work of Breiten (2015) using frequency domain formulations of Gramians for integro-differential equations.

Our goal will be to generate, for some r  n, a reduced dynamical system using a Petrov-Galerkin projection: select Vr ∈ Rn×r and Wr ∈ Rn×r such that WTr K(s)Vr is nonsingular in the right-half plane (ideally). Then, the reduced transfer function is Hr (s) = Cr (s)Kr (s)−1 Br (s) (2) where Cr (s) = C(s)Vr ∈ Cr ; Br (s) = WTr B(s) ∈ Cr ; and Kr (s) = WTr K(s)Vr ∈ Cr×r . Consider the aforementioned delay example H(s) = C(sE − A0 − A1 e−τ s )−1 B. In this case, the reduced transfer function will have the form Hr (s) = CVr (WTr EVr − WTr A0 Vr − WTr A1 Vr e−τ s )WTr B. The associated reduced dynamical model preserves the original system structure yet the dynamics evolve in a much smaller state-space. We will choose Wr and Vr so as to enforce rational interpolation conditions; Hr (s) will interpolate H(s) at selected points in the complex plane.

Here we use interpolation. Beattie and Gugercin (2009) established the framework for interpolatory model reduction of systems with a generalized realization. The following result describes how to construct reduced-order interpolants via projection. Theorem 1. (Beattie and Gugercin, 2009, Theorem 3) Let H(s) = C(s)K(s)−1 B(s) be given as in (1). Suppose r distinct points, {si }ri=1 , are chosen in the right half plane. Define Vr ∈ Cn×r and Wr ∈ Cn×r as: Vr = [K(s1 )−1 B(s1 ), · · · , K(sr )−1 B(sr )] and   K(s1 )−1 C(s1 )T   .. WTr =  . .

We will assume that H(s) is an H2 -function where H2 denotes the set of complex functions H(s) that are analytic in the open right half plane {s = x + y ∈ C : x > 0} and  +∞ such that supx>0 −∞ |H(x + ıy)|2 dy < ∞. Note that H2 is a Hilbert space with the inner product  +∞ 1 H∗ (ıω)G(ıω)dω (3) G, HH2 = 2π −∞ and the corresponding norm   +∞  1 | H(ıω) |2 dω . (4) HH2 = H, HH2 = 2π −∞

K(sr )−1 C(sr )T

Define Kr (s) = WTr K(s)Vr and assume that Kr (si ) is nonsingular for i = 1, 2, 3, ..., r. Define further:

(5) Br (s) = WTr B(s), and Cr (s) = C(s)Vr . Then with Hr (s) = Cr (s)Kr (s)−1 Br (s) we have H(si ) = Hr (si ) and H (si ) = Hr (si ) (6) for i = 1, ..., r where H (si ) denotes the first derivative of H(s) evaluated at si and Hr (si ) denotes the first derivative of Hr (s) evaluated at si .

We will look for structure-preserving reduced order models that are high-fidelity approximations in the H2 norm.

The rest of the paper is organized as follows: In §2, we give the necessary background for model reduction that allows passage from generalized realizations as in (1) to structure-preserving reduced models as in (2). In §3, we introduce an algorithm for determining structurepreserving interpolatory reduced models that are highfidelity approximations with respect to the H2 norm. We present two numerical examples in §4.

Once the interpolation points have been selected, Theorem 1 describes precisely how one may construct a structurepreserving interpolant via a Petrov-Galerkin projection. 2.2 Loewner framework for Hermite Interpolation

2. INTERPOLATORY MODEL REDUCTION OF GENERALIZED REALIZATIONS

The projection-based framework of Theorem 1 requires access to internal dynamics to construct the model reduction bases Vr and Wr . However in some applications, one does not have access to a description of internal dynamics and only transfer function measurements are available. The Loewner framework of Mayo and Antoulas (2007) resolves this issue: It only requires evaluation of the transfer function (and its derivative H (s) if interpolation points are repeated). As long as the transfer function H(s) = C(s)K(s)−1 B(s) can be sampled, one can produce a rational reduced model that interpolates the original system. Given a transfer function H(s) and a set of initial points {si }ri=1 , define  − (H(si ) − H(sj )) if i = j (7) (Er )i,j := si − sj  −H (si ) if i = j

We review here the basic interpolatory framework for building reduced models for generalized realizations having the form in (1) that preserve structure as in (2). These tools form the foundation of our approach described in §3. 2.1 Interpolatory Projections for Generalized Realizations For linear dynamical systems with generic first-order realizations, i.e., H(s) = C(sE − A)−1 B, various model reduction methods exist to produced high-fidelity or in some cases optimal reduced models. Examples include Gramian based methods such as Balanced Truncation (Mullis et al. (1976); Moore (1981)) and Hankel Norm Approximation (Glover (1984)), or interpolatory methods such as Iterative Rational Krylov Algorithm (Gugercin et al. (2008)). For 57

2016 IFAC SSSC 58 June 22-24, 2016. Istanbul, Turkey

(Ar )i,j

Klajdi Sinani et al. / IFAC-PapersOnLine 49-9 (2016) 056–061

 − (si H(si ) − sj H(sj )) := si − sj  −[sH(s)] |s=si

if i = j

the first-order necessary conditions for optimality. For SISO systems that we consider here, the reduced model is guaranteed to be a local minimizer Flagg et al. (2012).

(8)

if i = j

(9) and Cr = BTr = [H(s1 ) · · · H(sr )] The matrix Er is known as the Loewner matrix and Ar is known as the shifted Loewner matrix. Then, the transfer function Hr (s) = Cr (sEr − Ar )−1 Br is a Hermite interpolant of H(s), i.e., it satisfies (6). Note that the generalized coprime projection framework of Theorem 1 constructs an interpolant that retains the structure of the original dynamics, yet it requires access to internal dynamics via explicit knowledge of a structured realization. On the other hand, the Loewner framework produces an interpolant using only transfer function measurements. The Loewner interpolant does not preserve structure and instead will produce a standard first-order realization. For some recent work on structure-preserving Loewner interpolants, see Schulze and Unger (2015); Duff et al. (2015).

Algorithm: TF-IRKA IRKA using transfer function evaluations Given: H(s); and initial interpolation points {s1 , . . . , sr } . Return: Hr (s) = Cr (sEr − Ar )−1 Br , a reduced rational approximation satisfying H2 -optimality conditions (10).

(1) repeat until convergence (a) Construct Er , Ar , Cr and Br as in (7), (8), and (9). (b) Compute the eigenvalues {λi }ri=1 of the reduced-matrix pencil λEr −Ar , i.e., the poles of Hr (s). (c) Update the interpolation points: σi ← −λi . (2) Construct Er , Ar , Cr and Br as in (7), (8), and (9).

3. A STRUCTURE-PRESERVING INTERPOLATION BASED ALGORITHM

2.3 Realization Independent Rational H2 Approximation

We have discussed two interpolation-based methods for model reduction of a dynamical system with generalized realization. On the one hand, the formulation in Theorem 1 uses projection and produces a structurepreserving interpolant retaining the same general realization structure as in the original. For example, if the original system has an internal delay, the interpolant has the same delay structure, yet evolves in a smaller state space. However, even though the reduced model is an exact interpolant, the resulting interpolation points do not necessarily posses optimality. On the other hand, Algorithm TF-IRKA produces an interpolatory reduced model that satisfies the first-order necessary conditions for H2 optimality yet the reduced-model is only a rational approximation and does not preserve the inherent physical structure. Preserving the internal structure of a dynamical system is crucial since the reduced model, then, has the same physical interpretation as the original. Moreover, structure preservation prevents mixing of state variables which are associated with distinct physical features of the system (such as displacement and velocity in the case of second-order dynamical systems). Also, in some cases, the underlying structure brings additional properties. For example, port-Hamiltonian systems are guaranteed to be stable and passive. If the model reduction process retains port-Hamiltonian structure, the reduced port-Hamiltonian system is automatically guaranteed to be stable and passive; see, e.g., Gugercin et al. (2012); Polyuga and van der Schaft (2011). Also, see Rapisarda and Antoulas (2015); Rapisarda and van der Schaft (2013) for structured model reduction of port-Hamiltonian systems in a Loewner framework. Therefore, preserving structure when approximating a large-scale system is highly desired. Although some necessary conditions are known, e.g., see Beattie and Benner (2014), a full characterization of optimality is as yet unknown, unlike the case for rational approximation (without structure preservation) where optimal interpolation conditions are given in (10), In this section, we offer a heuristic approach that appears to narrow this gap.

In §2.1 and §2.2, the interpolation points {si } are assumed to be given without a notion of optimality. In this section, we briefly review how to choose {si } for H2 optimal approximation. However, the optimality will require the reduced-order model to be a rational approximant Hr (s) = Cr (sEr − Ar )−1 Br as in §2.2 without structure preservation. The following theorem gives the interpolatory first-order necessary conditions for optimal rational approximation in the H2 norm. Theorem 2. (Meier and Luenberger (1967); Gugercin et al. (2008)) Let Hr (s) = Cr (sEr − Ar )−1 Br be the best rth order rational approximation of H(s) with respect to the H2 norm, i.e.,     r (s) H(s) − Hr (s)H2 = min H(s) − H  . H2  r (s))=r dim(H Then H(−λk ) = Hr (−λk ), and H (−λk ) = Hr (−λk ) (10) for k = 1, 2, ..., r where λk denotes the poles of the reduced system Hr (s). For a more detailed discussion of the H2 optimal approximation, see, e.g., Meier and Luenberger (1967); Gugercin et al. (2008); Beattie and Gugercin (2015); Wilson (1970). Theorem 2 reveals that an H2 optimal rational approximant Hr (s) needs to be a Hermite interpolant to H(s) at the mirror images of its own poles. Since the interpolation points depend on the reduced model to be computed, satisfying them requires an iterative process. Iterative Rational Krylov Algorithm (IRKA) Gugercin et al. (2008) precisely achieves this task. The original formulation of IRKA was projection-based and required H(s) to have the generic first order form. Recently Beattie and Gugercin (2012) incorporated the Loewner framework into IRKA, leading to an iterative algorithm for rational H2 optimal approximation for systems with the generalized realizations. Since the method uses only transfer function (TF) evaluations, it is called TF-IRKA. A brief sketch of this method is given in Algorithm TF-IRKA below. Upon convergence, TF-IRKA produces a rational approximation satisfying

To make the presentation clearer, we will illustrate the proposed method on the Hadeler model in §4. For this 58

2016 IFAC SSSC June 22-24, 2016. Istanbul, Turkey

Klajdi Sinani et al. / IFAC-PapersOnLine 49-9 (2016) 056–061

model, we have H(s) = C((es − 1)A2 + s2 A1 + I)−1 B. As in TF-IRKA, choose r initial interpolation points {si }ri=1 . Then, using Theorem 1, we construct the projection-based model reduction bases Vr and Wr and obtain the reduced model as in (5). In this case the interpolatory reduced model will have the form 

2 T T Hr (s) = CVr (es − 1)WT r A2 Vr + s Wr A1 Vr + Wr Vr

−1

59

The two numerical examples in the next section illustrate this. Algorithm: SPTF-IRKA Generalized Realizations and TF-IRKA

WT r B.

Given: H (s)=C(s)K(s)−1 B(s); and initial interpolation points {s1 , . . . , sr } . Return: Hr (s) = Cr (s)Kr (s)−1 Br (s), a reduced system preserving the structure of the original system.

(1) repeat until convergence (a) Compute Vr and Wr as shown above (b) Compute Kr (s) = WTr K(s)Vr , Br (s) = WTr B(s), and Cr (s) = C(s)Vr (c) Hr (s) = Cr (s)Kr (s)−1 Br (s) (d) Compute a rational H2 -approximation Grat (s) = Cr (sEr − Ar )−1 Br to Hr (s) via Algorithm TF-IRKA (e) Compute the eigenvalues {λi }ri=1 of the reduced-matrix pencil λEr −Ar , i.e., the poles of Grat (s) (f) Update the interpolation points: σi ← −λi (2) Compute Kr (s) = WTr K(s)Vr , Br (s) = WTr B(s), and Cr (s) = C(s)Vr

Due to the structure of Hr (s), we cannot continue the iteration described in Algorithm TF-IRKA. Note that in Steps 1.b and 1.c of Algorithm TF-IRKA, one computes the poles of Hr (s) and reflects them across the imaginary axis to produce the next set of interpolation points. When Hr (s) is an rth order rational map described via a canonical realization, the next set of interpolation points may be computed straighforwardly through resolution of an rth order matrix eigenvalue problem. However, when Hr (s) has a noncanonical structure inherited, say, from H(s), then proceeding from this point becomes problematic. For the Hadeler model we consider, note that Hr (s) still has infinitely many poles, notwithstanding that the reduced system evolves in a state-space of much lower dimension. We cannot compute all the poles of Hr (s) and even if we could, we don’t know which r of them to reflect for use as interpolation points in the next iteration.

4. NUMERICAL EXAMPLES In this section we apply SPTF-IRKA to two models with generalized realizations and compare the results with TFIRKA where possible.

Instead, we need an effective strategy to produce r “good” interpolation points based on Hr (s). We propose introducing a second-layer of model reduction by using TFIRKA to obtain a standard rth-order optimal rational interpolant to Hr (s). We denote this rational approximant by Grat (s). Note that TF-IRKA only uses transfer function evaluations and so, demands no particular structure in Hr (s). Moreover, since Hr (s) is itself a reduced model, the cost of this intermediate TF-IRKA step is negligible. Since Grat (s) is a rational map with degree r, its r poles may be computed and reflected to produce interpolation points for the next structured interpolant, Hr (s). We call this iteration Structure Preserving TF-IRKA (SPTF-IRKA).

4.1 Hadeler Model We consider the Hadeler model from the Nonlinear Eigenvalue Problems (NLEVP) package provided in Betcke et al. (2011). This model first was introduced as a nonlinear eigenvalue problem in Hadeler (1967). The eigenvalue problem satisfies a generalized form of overdamping condition that ensures the existence of a complete set of eigenvectors. K(s) is given by the problem and we pick the B and C vectors to obtain the transfer function H(s) = C((es − 1)A2 + s2 A1 + I)−1 B where A2 , A1 are symmetric 100×100 matrices, I is the identity matrix, and B ∈ R100 . We choose C = BT . For more details about this model, see Betcke et al. (2011); Hadeler (1967). We applied both TF-IRKA and SPTF-IRKA to the Hadeler model and reduced the order from n = 100 to r = 8, r = 10, and r = 12. The results in Table 1 show that SPTF-IRKA outperforms TF-IRKA drastically.

Since Hr (s) is a (structured) Hermite interpolant of H(s), say at r distinct points, {ˆ s1 , . . . , sˆr }, the error system H − Hr is H2 -orthogonal to a 2r-dimension subspace,   1 1 1 1 span ,..., , ,..., . s − sˆ1 s − sˆr (s − sˆ1 )2 (s − sˆr )2 Different choices of {ˆ s1 , . . . , sˆr } will induce different error systems, H(s) − Hr (s). One particularly successful choice of {ˆ s1 , . . . , sˆr } are the H2 -optimal interpolation points and upon convergence of SPTF-IRKA, Grat (s) will be a (locally) H2 -optimal rth order approximation to H(s). Furthermore, the related error system H(s) − Grat (s) will be H2 -orthogonal to the subspace described above with {ˆ s1 , . . . , sˆr } determined to be H2 -optimal interpolation points. SPTF-IRKA induces orthogonality of the structured error system, H(s) − Hr (s), to this same subspace.

Table 1. SPTF-IRKA vs. TF-IRKA: Relative H∞ errors for Hadeler model r 8 10 12

Note that the optimality of this or any other choice of subspace with respect to the class of similarly structured reduced models is still unknown. Nonetheless, the final structure-preserving reduced model may be expected to be a high-fidelity approximation of the original system.

SPTF-IRKA 4.9501 · 10−5 5.3172 · 10−6 4.6956 · 10−6

TF-IRKA 0.9578 0.9579 0.9580

The better performance of SPTF-IRKA is also illustrated by the amplitude Bode plots in Figures 1 and 2, and the amplitude Bode plot of the error model in Figure 3. 59

2016 IFAC SSSC 60 June 22-24, 2016. Istanbul, Turkey

Klajdi Sinani et al. / IFAC-PapersOnLine 49-9 (2016) 056–061

results shown in Tables 2 and 3 illustrate the success of SPTF-IRKA for this model as well. TF-IRKA leads to smaller error; indeed, for this example both reduced models are rational approximations and TF-IRKA yields locally optimal approximations without being restricted by structure preservation. For TF-IRKA with generic firstorder structure, the reduced model has order 16.

SPTF-IRKA Bode plots

10 2

Full model Reduced model 10 0

H(iw)

10 -2

10

-4

Table 2. SPTF-IRKA: Relative H2 error for beam model

10 -6

10

-8

10 -10 10 -2

10

-1

10

0

10

1

10

2

10

3

10

r 8

4

SPTF-IRKA 8.8186 · 10−4

Frequencies

Fig. 1. Hadeler model: SPTF-IRKA (r = 10)

Table 3. SPTF-IRKA: Relative H∞ error for beam model r 8

SPTF-IRKA 6.2263 · 10−4

TF-IRKA Bode plots

10 2

TF-IRKA 4.0076 · 10−5

TF-IRKA 8.2372 · 10−6

Full model Reduced model

The amplitude Bode plots of the full and reduced model obtained via SPTF-IRKA and TF-IRKA are given, respectively, in Figure 4, and Figure 5.

10 0

H(iw)

10 -2

10

-4

SPTF-IRKA Bode plots

10 8

Full model Reduced model

10 -6 10 6

10 -8

10

10 4

-10

10 -2

10 -1

10 0

10 1

10 2

10 3

10 2

10 4

Frequencies 10

Fig. 2. Hadeler model: TF-IRKA (r = 10)

0

10 -2

10

-4

10 -6 -4 10

10

-3

10

-2

10

-1

10

0

10

1

10

2

10

3

10

4

Error Bode plots

10 2

SPTF-IRKA TF-IRKA 10 0

10

Fig. 4. Amplitude Bode Plot for the beam model: SPTFIRKA (r = 8)

-2

10 -4

10 -6

10 -8

TF-IRKA Bode plots

10 8

Full model Reduced model

10 -10 10 6

10 -1

10 0

10 1

10 2

10 3

10 4

10 4

H(iω)

10 -12 10 -2

Fig. 3. Error Bode plots for model reduction of the Hadeler model with SPTF-IRKA and TF-IRKA (r=10)

10 2

10 0

10 -2

10 -4 10 -4

4.2 Beam Model

10 -3

10 -2

10 -1

10 0

10 1

10 2

10 3

10 4

Frequencies

We consider a second-order Cantilever beam model with transfer function H(s) = C(s2 M + sG + K)−1 B where M, G, K ∈ R200×200 and B = CT ∈ R200 . For more information about this system, see Wyatt (2012). We reduced the model order from n = 200 to r = 8. The

Fig. 5. Amplitude Bode Plot for the beam model: TFIRKA (r = 8)

60

2016 IFAC SSSC June 22-24, 2016. Istanbul, Turkey

Klajdi Sinani et al. / IFAC-PapersOnLine 49-9 (2016) 056–061

5. CONCLUSION

61

Glover, K. (1984). All optimal hankel-norm approximations of linear multivariable systems and their l∞ -error bounds. Internat. J. Control, 39(6), 1115–1193. Gugercin, S., Beattie, C., and Antoulas, A.C. (2008). H2 model reduction for large-scale linear dynamical systems. Siam J. Matrix Anal. Appl., 30(2), 609–638. Gugercin, S., Polyuga, R., Beattie, C., and Van der Schaft, A. (2012). Structure-preserving tangential interpolation for model reduction of port-hamiltonian systems. Automatica, 48(9), 1963–1974. Hadeler, K.P. (1967). Mehrparametrige und nichtlineare eigenwertaufgaben. Arch. Rational Mech. Anal., 27(4), 306–328. Mayo, A. and Antoulas, A. (2007). A framework for the solution of the generalized realization problem. Linear algebra and its applications, 425(2), 634–662. Meier, L. and Luenberger, D. (1967). Approximation of linear constant systems. IEE. Trans. Automat. Contr., 12, 585–588. Meyer, D.G. and Srinivasan, S. (1996). Balancing and model reduction for second-order form linear systems. Automatic Control, IEEE Transactions on, 41(11), 1632–1644. Moore, B. (1981). Principal component analysis in linear systems: Controllability, observability, and model reduction. IEEE Transactions on Automatic Control, 26(1), 17–32. Mullis, C.T., Roberts, R., et al. (1976). Synthesis of minimum roundoff noise fixed point digital filters. Circuits and Systems, IEEE Transactions on, 23(9), 551–562. Polyuga, R. and van der Schaft, A. (2011). Structure preserving moment matching for port-hamiltonian systems: Arnoldi and lanczos. Automatic Control, IEEE Transactions on, 56(6), 1458–1462. doi: 10.1109/TAC.2011.2128650. Rapisarda, P. and Antoulas, A. (2015). A duality perspective on loewner rational interpolation and state-space modelling of vector-exponential trajectories. Proceedings of IEEE-CDC. Rapisarda, P. and van der Schaft, A. (2013). Datadriven identification and reduced-order modeling for linear conservative port- and self-adjoint hamiltonian systems. Proceedings of IEEE-CDC. Reis, T. and Stykel, T. (2008). Balanced truncation model reduction of second-order systems. Mathematical and Computer Modelling of Dynamical Systems, 14(5), 391– 406. Schulze, P. and Unger, B. (2015). Data-driven interpolation of dynamical systems with delay. PreprintReihe des Instituts f¨ ur Mathematik Technische Universit¨ at Berlin, Preprint 17-2015. Su, T.J. and Craig Jr, R. (1991). Model reduction and control of flexible structures using krylov vectors. Journal of guidance, control, and dynamics, 14(2), 260– 267. Wilson, D. (1970). Optimum solution of model-reduction problem. Proc. IEE, 117(6), 1161–1165. Wyatt, S. (2012). Issues in interpolatory model reduction: Inexact solves, second- order systems and daes. PhD Dissertation, Virginia Tech.

In this paper, we considered the model reduction problem for dynamical systems with generalized realizations whose frequency dependence is nonlinear. We developed a bilevel model reduction strategy that combines TFIRKA with the projection-based generalized coprime reduction to construct structure-preserving interpolatory reduced-order models. During the iteration, the interpolation points are updated by an H2 -optimal rational approximation of an intermediate model. The two numerical examples illustrate the efficiency of the proposed method. REFERENCES Antoulas, A.C. (2005). Approximation of large-scale dynamical systems (Advances in Design and Control). Society for Industrial and Applied Mathematics, Philadelphia, PA, USA. Antoulas, A.C., Beattie, C.A., and Gugercin, S. (2010). Interpolatory model reduction of large-scale dynamical systems. Efficient Modeling and Control of Large-Scale Systems, 3–58. Bai, Z. and Su, Y. (2005). Dimension reduction of largescale second order dynamical systems via a second-order arnoldi method. SIAM J. Sci. Comput., 26(5), 1692– 1709. Baur, U., Benner, P., and Feng, L. (2014). Model order reduction for linear and nonlinear systems: a systemtheoretic perspective. Archives of Computational Methods in Engineering, 21(4), 331–358. Beattie, C. and Gugercin, S. (2009). Interpolatory projection methods for structure-preserving model reduction. Systems and Control Letters, 58, 225–232. Beattie, C. and Gugercin, S. (2012). Realizationindependent H2 -approximation. 51st IEEE Conference on Decision and Control, Maui, HI. Beattie, C.A. and Gugercin, S. (2015). Model reduction by rational interpolation. In P. Benner, A. Cohen, M. Ohlberger, and K. Willcox (eds.), To appear in Model Reduction and Approximation: Theory and Algorithms. Available as http://arxiv.org/abs/1409.2140. SIAM, Philadelphia. Beattie, C. and Benner, P. (2014). H2 -optimality conditions for structured dynamical systems. MPI Magdeburg Preprint MPIMD/14-18, October 2014. Betcke, T., Higham, N.J., Mehrmann, V., Shr¨ oder, C., and Tisseur, F. (2011). NLEVP: A collection of nonlinear eigenvalue problems. Manchester Institute for Mathematical Sciences. Breiten, T. (2015). Structure-preserving model reduction for integro-differential equations. Submitted. Chahlaoui, Y., Gallivan, K.A., Vandendorpe, A., and Van Dooren, P. (2005). Model reduction of secondorder systems. In Dimension Reduction of Large-Scale Systems, 149–172. Springer. Duff, I.P., Poussot-Vassal, C., and Seren, C. (2015). Realization independent single time-delay dynamical model interpolation and H2 -optimal approximation. Proceedings of the 54th IEEE Conference on Decision and Control, Atlanta, GA. Also available as arXiv:1504.06457v1. Flagg, G., Beattie, C., and Gugercin, S. (2012). Convergence of the Iterative Rational Krylov Algorithm. Systems & Control Letters, 61(6), 688–691. 61