Linear time-periodic dynamical systems: An H2 analysis and a model reduction framework

TL;DR

This paper introduces a new model reduction framework for linear time-periodic systems using an extended H2 analysis, enabling efficient simulation of complex systems in fluid dynamics, electronics, and mechanics.

Contribution

It generalizes H2 optimal model reduction methods from LTI to LTP systems, providing an algorithm that preserves periodic structure and includes an a posteriori error bound.

Findings

Successfully applied to numerical examples demonstrating effectiveness.

Retains periodic structure in reduced models.

Provides an a posteriori error bound for model accuracy.

Abstract

Linear time-periodic (LTP) dynamical systems frequently appear in the modeling of phenomena related to fluid dynamics, electronic circuits, and structural mechanics via linearization centered around known periodic orbits of nonlinear models. Such LTP systems can reach orders that make repeated simulation or other necessary analysis prohibitive, motivating the need for model reduction. We develop here an algorithmic framework for constructing reduced models that retains the linear time-periodic structure of the original LTP system. Our approach generalizes optimal approaches that have been established previously for linear time-invariant (LTI) model reduction problems. We employ an extension of the usual H2 Hardy space defined for the LTI setting to time-periodic systems and within this broader framework develop an a posteriori error bound expressible in terms of related LTI systems. Optimization of this bound motivates our algorithm. We illustrate the success of our method on two numerical examples.

Equations

$$\displaystyle{\mathcal{G}}:\left\{\begin{aligned} \dot{{\mathbf{x}}}(t)&={\mathbf{A}}(t){\mathbf{x}}(t)+{\mathbf{b}}(t)u(t)\\ y(t)&={\mathbf{c}}^{T}(t){\mathbf{x}}(t)\end{aligned}\right.\left(\mbox{abbreviated }{\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{A}}(t)&{\mathbf{b}}(t)\\ \hline\cr{\mathbf{c}}^{T}(t)&\end{array}\right]\right),$$

$${\mathsf{L}}_2:=\left\{u:{\int }_0^{\infty }|u(t){|}^{2}dt<\infty \right\}.$$

$$\displaystyle\widetilde{\mathcal{G}}:\left\{\begin{aligned} \dot{\widetilde{\mathbf{x}}}(t)&=\widetilde{{\mathbf{A}}}(t)\widetilde{{\mathbf{x}}}(t)+\widetilde{{\mathbf{b}}}(t)u(t)\\ \widetilde{y}(t)&=\widetilde{{\mathbf{c}}}^{T}(t)\widetilde{{\mathbf{x}}}(t)\end{aligned}\right.\left(\mbox{abbreviated }\tilde{\mathcal{G}}=\left[\begin{array}[]{c|c}\tilde{\mathbf{A}}(t)&\tilde{\mathbf{b}}(t)\\ \hline\cr\tilde{\mathbf{c}}^{T}(t)&\mathbf{0}\end{array}\right]\right),$$

$$G:\{\begin{array}{rl}\dot{\mathbf{x}}(t)& =\mathbf{A}\mathbf{x}(t)+\mathbf{Bu}(t)\\ \mathbf{y}(t)& =\mathbf{C}\mathbf{x}(t)\end{array}$$

$$\mathbf{G}(s)=\mathbf{C}[s\mathbf{I}-\mathbf{A}{]}^{-1}\mathbf{B}$$

$$\tilde{G}:\{\begin{array}{rl}\dot{\tilde{\mathbf{x}}}(t)& =\tilde{\mathbf{A}}\tilde{\mathbf{x}}(t)+\tilde{\mathbf{B}}\mathbf{u}(t)\\ \tilde{\mathbf{y}}(t)& =\tilde{\mathbf{C}}\tilde{\mathbf{x}}(t)\end{array}$$

$$\tilde{\mathbf{G}}(s)=\tilde{\mathbf{C}}(s\tilde{\mathbf{I}}-\tilde{\mathbf{A}}{)}^{-1}\tilde{\mathbf{B}}$$

$${\mathbf{W}}^T\left(\mathbf{V}\dot{\tilde{\mathbf{x}}}(t)-\mathbf{A}\mathbf{V}\tilde{\mathbf{x}}(t)-\mathbf{B}\mathbf{u}(t)\right)=0,\tilde{\mathbf{y}}(t)=\mathbf{CV}\tilde{\mathbf{x}}(t).$$

$$\tilde{\mathbf{A}}={\mathbf{W}}^T\mathbf{A}\mathbf{V},\tilde{\mathbf{B}}={\mathbf{W}}^T\mathbf{B},\text{mbox}and\tilde{\mathbf{C}}=\mathbf{C}\mathbf{V}$$

$$⟨\mathbf{G},\mathbf{H}{⟩}_{{\mathsf{H}}_2}=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\text{mbox}trace\left(\overline{\mathbf{G}(\dot{ıı}\omega )}{\mathbf{H}}^T(\dot{ıı}\omega )\right)d\omega =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\text{mbox}trace\left(\mathbf{G}(-\dot{ıı}\omega ){\mathbf{H}}^T(\dot{ıı}\omega )\right)d\omega$$

$${\Vert \mathbf{G}\Vert }_{{\mathsf{H}}_2}={\left(\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\Vert \mathbf{G}(\dot{ıı}\omega ){\Vert }_F^{2}d\omega \right)}^{1/2}={\left({\int }_0^{\infty }\text{mbox}trace({\mathbf{g}}^{\ast }(t)\mathbf{g}(t))dt\right)}^{1/2}$$

$${\Vert \mathbf{y}\Vert }_{{\mathsf{L}}_{\infty }}\le {\Vert \mathbf{G}\Vert }_{{\mathsf{H}}_2}{\Vert \mathbf{u}\Vert }_{{\mathsf{L}}_2}.$$

$${⟨\mathbf{G},\mathbf{H}⟩}_{{\mathsf{H}}_2}=\sum _{k=1}^n\text{mbox}res[\text{mbox}Tr\left(\overline{\mathbf{G}}(-s)\mathbf{H}(s{)}^T\right),{\mu }_k].$$

$$⟨\mathbf{G},\mathbf{H}{⟩}_{{\mathsf{H}}_2}=\sum _{k=1}^n{\mathsf{c}}_k^T\overline{\mathbf{G}}(-{\mu }_k){\mathsf{b}}_k\text{mbox}and{\Vert \mathbf{H}\Vert }_{{\mathsf{H}}_2}={\left(\sum _{k=1}^n{\mathsf{c}}_k^T\overline{\mathbf{H}}(-{\mu }_k){\mathsf{b}}_k\right)}^{1/2}.$$

$$y(t)$$

$$g(t+T,\tau +T)=g(t,\tau ),\text{mbox}forallt\ge \tau .$$

$$u(t)=\sin (\omega t)\stackrel{y=\mathcal{G}u}{⟹}y(t)\to |G(\dot{ıı}\omega )|\sin [\omega t+\text{mbox}arg(G(\dot{ıı}\omega ))].$$

$$g(t,t-r)=\sum _{k=-\infty }^{\infty }{g}_k(r)e^{\dot{ıı}k{\omega }_{0}t},\text{mbox}where{g}_k(r)=\frac{1}{T}{\int }_0^T{e}^{-\dot{ıı}k{\omega }_{0}t}g(t,t-r)dt.$$

$$u(t)=\sin (\omega t)\stackrel{y=\mathcal{G}u}{⟹}y(t)\to \sum _{k\in \mathbb{Z}}|g_k(\dot{ıı}\omega )|\sin [(\omega +k{\omega }_0)t+\text{mbox}arg(g_k(\dot{ıı}\omega ))].$$

$$\dot{\mathbf{x}}=\mathbf{A}(t)\mathbf{x},\mathbf{A}(t)=\mathbf{A}(t+T)$$

$$\mathbf{X}(t)=\mathbf{P}(t)e^{\mathbf{Q}t},\text{mbox}where\mathbf{P}(t)\in {\mathbb{R}}^{n\times n}\text{mbox}with\mathbf{P}(t+T)=\mathbf{P}(t),$$

$$\mathcal{G}:\{\begin{array}{rl}\dot{\mathbf{x}}(t)& =\mathbf{A}(t)\mathbf{x}(t)+\mathbf{b}(t)u(t)\\ y(t)& ={\mathbf{c}}^T(t)\mathbf{x}(t)\end{array}⟹\mathcal{G}:\{\begin{array}{rl}\dot{\mathbf{z}}(t)& =\mathbf{Q}\mathbf{z}(t)+{\mathbf{P}}^{-1}(t)\mathbf{b}(t)u(t)\\ y(t)& ={\mathbf{c}}^T(t)\mathbf{P}(t)\mathbf{z}(t)\end{array}$$

$$⟨\mathcal{G},\mathcal{H}{⟩}_{{\mathsf{H}}_2}$$

$$⟨\mathcal{G},\mathcal{H}{⟩}_{{\mathsf{H}}_2}=\sum _{k\in \mathbb{Z}}⟨{\hat{g}}_k,{\hat{h}}_k{⟩}_{{\mathsf{H}}_2}$$

$${\Vert \mathcal{G}\Vert }_{{\mathsf{H}}_2}^2=\sum _{k\in \mathbb{Z}}{\Vert {\hat{g}}_k\Vert }_{{\mathsf{H}}_2}^2.$$

$$\mathbf{c}(t)=\sum _{k\in \mathbb{Z}}{\mathbf{c}}_k{e}^{-\dot{ıı}k{\omega }_{0}t}\text{mbox}and\mathbf{b}(t)=\sum _{k\in \mathbb{Z}}{\mathbf{b}}_k{e}^{-\dot{ıı}k{\omega }_{0}t}$$

$${\Vert {\hat{g}}_k\Vert }_{{\mathsf{H}}_2}$$

$$=\Vert [s\mathbf{I}-\mathbf{Q}{]}^{-1}{\Vert }_{{\mathsf{H}}_2}\sum _{\ell \in \mathbb{Z}}\Vert {\mathbf{c}}_{k-\ell }{\Vert }_2\Vert {\mathbf{b}}_{\ell }{\Vert }_2,$$

$$\displaystyle\left\{\sum_{\ell\in{\mathbb{Z}}}\|{\mathbf{c}}_{k-\ell}\|_{2}\|{\mathbf{b}}_{\ell}\|_{2}\right\}_{k=-\infty}^{\infty}=\Big{\{}\|{\mathbf{c}}_{k}\|_{2}\Big{\}}_{k=-\infty}^{\infty}*\Big{\{}\|{\mathbf{b}}_{k}\|_{2}\Big{\}}_{k=-\infty}^{\infty}\in\ell_{1}.$$

$$⟨G,H{⟩}_{{\mathsf{H}}_2}$$

$$\text{mbox}and\Vert G{\Vert }_{{\mathsf{H}}_2}^2$$

Full text

We develop here an algorithmic framework for constructing reduced models that retains the linear time-periodic structure of the original LTP system. Our approach generalizes optimal approaches that have been established previously for linear time-invariant (LTI) model reduction problems. We employ an extension of the usual ${{\mathsf{H}}_{2}}$ Hardy space defined for the LTI setting to time-periodic systems and within this broader framework develop an a posteriori error bound expressible in terms of related LTI systems. Optimization of this bound motivates our algorithm. We illustrate the success of our method on two numerical examples.

{classcode}

37M05, 37M99, 93A15, 65K99, 93C05

1 Introduction

We consider linear continuous-time-periodic (LTP) single-input/single-output (SISO) dynamical systems realized through a state-space representation,

[TABLE]

where ${\mathbf{A}}(t)\in{\mathbb{R}}^{n\times n}$ and ${\mathbf{b}}(t),{\mathbf{c}}(t)\in{\mathbb{R}}^{n}$ are $T$-periodic, with ${\mathbf{A}}(t)={\mathbf{A}}(t+T)$, ${\mathbf{b}}(t)={\mathbf{b}}(t+T)$ and ${\mathbf{c}}(t)={\mathbf{c}}(t+T)$ for some fixed $T>0$ and we assume that ${\mathbf{x}}(0)=0$. Let $\mathsf{L}_{2}$ denote the vector space of real-valued square integrable functions (having “bounded energy”).

[TABLE]

We assume that the system ${\mathcal{G}}$ in (3) is causal and is a bounded, linear mapping from $\mathsf{L}_{2}$-inputs to $\mathsf{L}_{2}$-outputs, ${\mathcal{G}}:\mathsf{L}_{2}\to\mathsf{L}_{2}$.

For any given order $r\ll n$, our goal is to find a reduced-order model,

[TABLE]

where $\widetilde{{\mathbf{x}}}(0)=0$ and with $\widetilde{{\mathbf{A}}}(t)\in{\mathbb{R}}^{r\times r}$ and $\widetilde{{\mathbf{b}}}(t),\widetilde{{\mathbf{c}}}(t)\in{\mathbb{R}}^{r}$ also $T$-periodic and chosen in such a way so that $\widetilde{y}(t)\approx y(t)$ over a wide class of inputs $u(t)$.

Starting with a Floquet transformation of the LTP system (3), our approach to this problem involves conversion of the LTP model reduction problem into a closely related linear time-invariant (LTI) multiple-input/multiple-output (MIMO) model reduction problem, connecting the model reduction error in the original LTP setting directly to a related LTI model reduction error. In particular, we are able to provide a posteriori bounds for the error between the full- and reduced-order LTP systems directly in terms of the error in an LTI MIMO model reduction counterpart. This in turn invites the use of the Iterative Rational Krylov Algorithm (IRKA) [1] as a way of minimizing the error bound we establish. We find that this approach generates superior performance over other applicable and commonly used model reduction techniques for the examples considered in this paper.

2 Background

2.1 Previous work on the analysis and reduction of LTP systems

A significant body of literature has developed focussing on model reduction of linear time-invariant (LTI) systems [2, 3, 4]. Analogous developments for continuous-time LTP systems have been far more limited and usually appear as a special case of the more general problem class of model reduction for general linear time-varying systems where computationally effective strategies for large scale problems have not yet emerged (see, e.g., [5, 6, 7]). Available strategies for this problem have focussed on extensions of balanced truncation to general linear time-varying systems. The computational challenges are formidable even for modest order and are typified by the need to solve two large-scale Lyapunov differential inequalities [5]. Notably, the recent work [7] is a step towards addressing this computational challenge. By way of contrast, it is worth noting that for the case of linear discrete-time periodic systems, strategies built upon balanced truncation have been somewhat more successful, producing strategies that still are computationally demanding but that could remain tractable for problems with modest order (see, e.g., [8, 9, 10, 11, 12]). Throughout this work, our focus is on continuous-time LTP systems.

Aside from model reduction, there is a significant body of literature that has developed on LTP systems focussing on a variety of other important topics including control, spectral analysis, and harmonic response, among others. Wereley and Hall [13, 14] develop a Fourier analysis of state-space systems with periodic matrices via the harmonic balance method. They represent the LTP system as generalization of a frequency response operator called the Harmonic Transfer Function which can be viewed as an operator on a family of functions called exponentially modulated periodic inputs. The operator is infinite-dimensional and consists of a countable number of time-invariant LTI components.

Sandberg et al. [15, 16] expand the impulse response of LTP systems via Fourier series expansion. The frequency response of the system is understood as a countable collection of infinite-dimensional LTI transfer functions called the Floquet-Fourier representation, making analytic expressions of frequency responses considerably easier to represent .

Zhou and Hagiwara [17, 18, 19, 20, 21] address dynamical system norm computation using truncations of the Harmonic Transfer Function. Convergence of the truncations is proven. Additionally, the authors discuss the spectral characteristics of LTP systems as input/output operators.

2.2 Linear Time-Invariant Dynamical Systems

We review briefly the basic features of linear time-invariant dynamical systems that are relevant to our approach and that will serve to establish our notation. In that context, we consider MIMO linear dynamical systems given in the following state-space form:

[TABLE]

where \mathbf{u}(t)\in{\mathbb{R}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{u}}}}, \mathbf{y}(t)\in{\mathbb{R}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{y}}}}, and $\mathbf{x}(t)\in{\mathbb{R}}^{n}$ are respectively, the inputs, outputs, and states of ${{G}}$ so that $\mathbf{A}\in{\mathbb{R}}^{n\times n}$, \mathbf{B}\in{\mathbb{R}}^{n\times{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{u}}}} and \mathbf{C}\in{\mathbb{R}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{y}}}\times n}. Let $\hat{\mathbf{Y}}(s)$ and $\hat{\mathbf{U}}(s)$ denote the Laplace transforms of $\mathbf{y}(t)$ and $\mathbf{u}(t)$, respectively. Then by taking the Laplace transform of (7), we obtain the transfer function $\mathbf{G}(s)$ of the system ${{G}}$, where $\hat{\mathbf{Y}}(s)=\mathbf{G}(s)\hat{\mathbf{U}}(s)$ and

[TABLE]

is the transfer function of the associated dynamical system ${{G}}$ in (7). $\mathbf{G}(s)$ is a matrix-valued function \mathbf{G}:{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{{\mathbb{C}}\to{\mathbb{C}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{u}}}}\times{\mathbb{C}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{y}}}}}} with entries consisting of proper rational functions with poles at the eigenvalues of ${\mathbf{A}}$.

2.2.1 Projection-Based Model Reduction

When the state-space dimension of ${{G}}$ is large, it may be useful to replace the original model (7) with a reduced model that has much smaller state space dimension yet still is able to replicate significant features of the original input/output dynamics. Toward this end, we seek a reduced model

[TABLE]

where $\widetilde{{\mathbf{A}}}\in{\mathbb{R}}^{r\times r}$, \widetilde{{\mathbf{B}}}\in{\mathbb{R}}^{r\times{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{u}}}}, and \widetilde{{\mathbf{C}}}\in{\mathbb{R}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{y}}}\times r} with $r\ll n$, and such that $\widetilde{\mathbf{y}}(t)\approx{\mathbf{y}}(t)$ for a large class of inputs $\mathbf{u}(t)$. In order to make the quality of this approximation largely independent of input, it becomes necessary that the transfer function of (9) given by

[TABLE]

should approximate $\mathbf{G}(s)$ well, to the extent that the error $\mathbf{G}(s)-\widetilde{\mathbf{G}}(s)$ is small with respect to natural metrics that we discuss below.

The most common way to obtain the reduced-order models as in (9) is via projection: Construct two matrices $\mathbf{V}\in{\mathbb{R}}^{n\times r}$ and $\mathbf{W}\in{\mathbb{R}}^{n\times r}$ such that $\mathbf{W}^{T}\mathbf{V}=\mathbf{I}_{r}$. Then, approximate the full-order state $\mathbf{x}(t)$ by $\mathbf{V}\widetilde{\mathbf{x}}(t)$, and enforce the Petrov-Galerkin condition, forcing orthogonality of the residual dynamics to the range of $\mathbf{W}$:

[TABLE]

Then the reduced model state-space matrices become

[TABLE]

Obviously the choice of ${\mathbf{V}}$ and ${\mathbf{W}}$ will determine the accuracy of the reduced model approximation.

2.2.2 The ${\mathsf{H}}_{2}$ Inner-product Space

Let ${\mathsf{H}}_{2}$ denote the set of n_{y}\times{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{u}}} matrix-valued functions ${\mathbf{H}}(s)$ analytic in the open right half-plane such that \sup_{x>0}\int_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{-\infty}}}^{\infty}\|{\mathbf{H}}(x+{\dot{\imath\!\imath}}y)\|_{F}^{2}\,dy<\infty. The ${\mathsf{H}}_{2}$ space is a Hilbert space endowed with the inner product

[TABLE]

with the associated norm

[TABLE]

where \bm{g}(t)\in{\mathbb{R}}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{y}}}\times{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{n_{u}}}} is the matrix-valued impulse response of the MIMO dynamical system, ${{G}}$. The ${{\mathsf{H}}_{2}}$ norm is of particular interest because it bounds the dynamical system output for bounded-energy inputs as,

[TABLE]

The next result, used frequently throughout this paper, provides an alternative representation of the ${\mathsf{H}}_{2}$ inner product using the residue calculus.

Theorem 2.1**.**

[1, Lemma 2.4, p.615]** Suppose that ${\mathbf{G}}(s)$ and ${\mathbf{H}}(s)$ are stable (poles contained in the open left halfplane) and suppose that ${\mathbf{H}}(s)$ has poles at $\mu_{1},\,\mu_{2},\,\ldots\,\mu_{n}$. Then

[TABLE]

In particular, if ${\mathbf{H}}(s)$ has only simple or semi-simple poles at $\mu_{1},\,\mu_{2},\,\ldots\,\mu_{n}$, then ${\mathbf{H}}(s)=\sum_{i=1}^{n}\frac{1}{s-\mu_{i}}\mathbf{\mathsf{c}}_{i}\mathbf{\mathsf{b}}_{i}^{T}$, where $\mathbf{\mathsf{c}}_{i}\mathbf{\mathsf{b}}_{i}^{T}$ is the residue of ${\mathbf{H}}(s)$ at $s=\mu_{i}$, and

[TABLE]

2.2.3 Optimal ${{\mathsf{H}}_{2}}$ Approximation

We seek reduced models for LTP systems that are accurate with respect to a metric that is equivalent to the ${\mathsf{H}}_{2}$ norm for LTI systems (which are special cases of LTP systems). Indeed, for the case of LTI systems, we are able to do this optimally, finding a degree-$r$ LTI reduced model with transfer function ${\mathbf{H}}_{r}(s)$ that minimizes (at least locally) the ${\mathsf{H}}_{2}$ error norm $\|{\mathbf{H}}-\widetilde{{\mathbf{H}}}_{r}\|_{{\mathsf{H}}_{2}}$ over all degree-$r$ reduced models with transfer functions $\widetilde{{\mathbf{H}}}(s)$. The optimal LTI ${\mathsf{H}}_{2}$ approximant satisfies certain Hermite tangential interpolation conditions; for details, see, e.g., [22, 1]. Significant for the present work is the availability of a numerically effective algorithm with which one may construct (local) ${\mathsf{H}}_{2}$-optimal approximants, namely the Iterative Rational Krylov Algorithm (IRKA) of [1]. IRKA is an iterative algorithm; it typically converges quite rapidly though its speed of convergence may slow as the number of inputs and outputs grows. Beattie and Gugercin in [23] developed a modified version of IRKA, in order to address this slowing of convergence and improved the performance of IRKA significantly for MIMO problems; see [23] for more details. Convergence of IRKA is guaranteed for special cases [24] even though there are known, but rare, cases where convergence may fail [1, 24]. Upon convergence, the resulting reduced model is guaranteed to be a local $\mathcal{H}_{2}$-minimizer. In [25], Beattie and Gugercin developed a trust-region framework for IRKA, which is globally convergent to a local $\mathcal{H}_{2}$-minimizer. The added guarantees this approach provides seem unnecessary for most problems at hand and the original formulation of IRKA has been successfully applied to large-scale problems in various application settings in order to provide (locally) optimal reduced models; see, e.g., [26, 27, 1]. We will employ IRKA as a critical step in the model reduction framework for LTP systems that we propose here.

2.3 Linear Time-Periodic Dynamical Systems

We view the time-periodic system (3) as a causal linear map ${\mathcal{G}}:{\mathsf{L}_{2}}\to{\mathsf{L}_{2}}$, that maps an input signal $u\in{\mathsf{L}_{2}}$ to an output signal $y\in{\mathsf{L}_{2}}$. This can be expressed via a convolution integral with an generalized impulse response $g(t,\tau)$ which, for $t\geq\tau$, gives the response of the system (3) to an impulsive input at $t=\tau$. Thus,

[TABLE]

Unlike LTI systems, the impulse response of a time-varying system depends on the time of the impulse $\tau$ and not simply on the time elapsed since the impulse was applied, $t-\tau$. Since the state-space parameters are $T$-periodic, so too are the impulse responses. That is,

[TABLE]

In this setting, causality implies that $g(t,\tau)=0$ for $t<\tau$.

2.3.1 Frequency Coupling

An important characteristic of stable LTI systems is that their steady state response to a single-frequency sinusoid is another sinusoid of the same frequency. That is, if $G(s)$ is the transfer function of an LTI system ${\mathcal{G}}$, then

[TABLE]

Note that we write “$y(t)\to$” to describe the steady-state behavior of $y(t)$ after the transient response has decayed.

By way of contrast, LTP systems produce a countable number of harmonics of the input frequency. Let $r=t-\tau$. Since $g(t,t-r)$ is $T$-periodic in $t$, the impulse response can be expanded into a Fourier series in $r$,

[TABLE]

We call $\{g_{k}\}$ the subsystems of ${\mathcal{G}}$. Now, let $u(t)=\sin(\omega t)$ and $\displaystyle\omega_{0}=\frac{2\pi}{T}$ be the fundamental frequency of the LTP system ${\mathcal{G}}$. Then,

[TABLE]

For a single frequency input, an LTP system will have output frequencies $\omega_{k}=\omega+k\omega_{0}$, where $k\in\mathbb{Z}$; see Sandberg [15].

2.3.2 Floquet Transformations

The Floquet transformation is a time-dependent change-of-variable that transforms the time-periodic differential equation

[TABLE]

to an equivalent system of differential equations with constant coefficients, $\dot{\mathbf{z}}={\mathbf{Q}}\mathbf{z}(t)$ via a periodic, time-dependent change of variables, $\mathbf{P}(t)\mathbf{z}(t)={\mathbf{x}}(t)$; see, e.g., [28] for details. In particular, we have

Theorem 2.2** (Floquet’s Theorem).**

Let ${\mathbf{A}}(t)\in\mathbb{R}^{n\times n}$ be continuous and periodic with period $T$ for $-\infty<t<\infty$ Then any fundamental matrix $\mathbf{X}(t)$ of the differential equation $\dot{{\mathbf{x}}}={\mathbf{A}}(t){\mathbf{x}}$ has a representation of the form

[TABLE]

and ${\mathbf{Q}}\in\mathbb{R}^{n\times n}$ is a constant matrix.

The constant matrix, ${\mathbf{Q}}$ is defined so that ${\mathbf{Q}}=\frac{1}{T}\mathsf{log}(\mathbf{M}(T))$ where the matrix $\mathbf{M}(T)=\mathbf{X}(T)\mathbf{X}(0)^{-1}$ is the monodromy matrix associated with (16). Note that ${\mathbf{Q}}$ is stable precisely when $\mathbf{M}(T)$ is a contraction. The $T$-periodic matrix ${\mathbf{P}}(t)$ is defined as $\mathbf{P}(t)=\mathbf{X}(t)\mathbf{X}(0)^{-1}e^{-{\mathbf{Q}}t}$ for $t\geq 0$.

Given an LTP system, ${\mathcal{G}}$, a Floquet transformation using the periodic change of variable $\mathbf{z}(t)={\mathbf{P}}^{-1}(t){\mathbf{x}}(t)$ from Theorem 2.2, can be performed so that

[TABLE]

The Floquet transformation can be costly to determine for large-scale LTP systems. Methods have been developed recently to make the transformation computationally tractable for modest order; see, e.g., [29, 30] for a discussion. The main computational cost in most cases will be the calculation of the monodromy matrix, which necessitates the solution of $n$ independent initial value problems associated with (16). The remaining tasks scale with a complexity of $n^{3}$, so this can be feasible for modest orders of $n$. In this work, we will not consider the practical difficulties associated with constructing this transformation of the original LTP system; we will suppose that our LTP system is presented with a realization having the form (17) with constant ${\mathbf{Q}}$ and $T$-periodic input and output maps.

3 An ${{\mathsf{H}}_{2}}$ Analysis for LTP systems

We develop here for LTP systems the metrics and associated analysis that extend certain notions of ${{\mathsf{H}}_{2}}$-approximation that are well-established for LTI systems.

3.1 The Periodic ${{\mathsf{H}}_{2}}$ Inner Product

Let ${\mathcal{G}}$ and ${\mathcal{H}}$ be two LTP systems with fundamental frequency $\displaystyle\omega_{0}=\frac{2\pi}{T}$. The ${{\mathsf{H}}_{2}}$ inner product $\langle{\mathcal{G}},{\mathcal{H}}\rangle_{{\mathsf{H}}_{2}}$ is defined in terms of their impulse responses $g(t,\tau)$ and $h(t,\tau)$,

[TABLE]

The next results express this inner product in terms of the kernel subsystems; this can be found as Corollary 1 in Sandberg et al. [16]. Our formulation is framed in the frequency domain whereas the original formulation of Sandberg et al. [16] places it in the time domain.

Theorem 3.1**.**

Let ${\mathcal{G}}$ and ${\mathcal{H}}$ be two LTP systems with subsystems $g_{k}(t)$ and $h_{k}(t)$, respectively. Then,

[TABLE]

where $\hat{g}_{k}(s)$ and $\hat{h}_{k}(s)$ denote the Laplace transformations of $g_{k}(t)$ and $h_{k}(t)$, respectively, and $\langle\hat{g}_{k},\hat{h}_{k}\rangle_{{{\mathsf{H}}_{2}}}$ is the regular ${{\mathsf{H}}_{2}}$-inner product for LTI systems defined in (12).

It follows immediately that

[TABLE]

This expression for the ${{\mathsf{H}}_{2}}$ norm will allow us to develop a simple test for the boundedness of the ${{\mathsf{H}}_{2}}$ norm of an LTP system ${\mathcal{G}}$. We adapt the following (standard) definition to our context:

Definition 3.2**.**

For $\alpha>0$, $\mathsf{Lip}(\alpha)$ denotes the set of continuous periodic functions $f(t)$ (say, with period T and fundamental frequency $\omega_{0}=\frac{2\pi}{T}$) such that for some finite $M>0$, $|f(t_{1})-f(t_{2})|\leq M|t_{1}-t_{2}|^{\alpha}$ uniformly for all $t_{1},t_{2}\in[0,T]$.

Note that $\mathsf{Lip}(\alpha)$ with $\alpha>1$ consists only of constant functions and $\mathsf{Lip}(1)$ consists of Lipschitz continuous periodic functions which in turn will be contained within $\mathsf{Lip}(\alpha)$ for any $\alpha\in(0,1)$.

Our main result related to the ${{\mathsf{H}}_{2}}$ analysis of LTP systems makes use of the following two lemmata:

Lemma 3.3**.**

*([31, Theorem 1, p.176]).

Let $\ell_{1}$ and $\ell_{2}$ denote the usual normed spaces of absolutely summable and square summable (scalar) sequences, respectively. If {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\{a_{k}\}}}\in\ell_{1} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\{x_{k}\}}}\in\ell_{1} then the convolution $y_{k}=\sum_{\ell}a_{k-\ell}x_{\ell}$ is absolutely convergent; {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\{y_{k}\}}}\in\ell_{1}. If {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\{a_{k}\}}}\in\ell_{1} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\{x_{k}\}}}\in\ell_{2} then the convolution $y_{k}=\sum_{\ell}a_{k-\ell}x_{\ell}$ is an $\ell_{2}$ sequence; {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\{y_{k}\}}}\in\ell_{2}.*

Lemma 3.4**.**

Let the vector-valued function ${\mathbf{x}}(t)\in{\mathbb{R}}^{n}$ be $T$-periodic with components in $\mathsf{Lip}(\alpha)$ with $\alpha>\frac{1}{2}$. Then

the Fourier expansion ${\mathbf{x}}(t)=\sum_{k\in{\mathbb{Z}}}{\mathbf{x}}_{k}e^{jk\omega_{0}t}$ is absolutely convergent and $\{\|{\mathbf{x}}_{k}\|_{2}\}\in\ell_{1}$ (Bernstein’s Theorem, **[32, Theorem VI.3.1]**). 2. 2.

the Fourier coefficients $\{{\mathbf{x}}_{k}\}$ satisfy $\sqrt{k}|{\mathbf{x}}_{k}|\rightarrow 0$ (componentwise) as $k\rightarrow\infty$ (**[32, II.4]**)

The following theorem can be understood as a special case of Lemma 1 of Sandberg et al. [15], which pertains to a more general class of impulse response functions. The proof provided here is instead a concise frequency domain-based proof for dynamical systems with state space representations.

Theorem 3.5**.**

Given system {\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{A}}(t)&{\mathbf{b}}(t)\\ \hline\cr\mathbf{c}^{T}(t)&\end{array}\right] with ${\mathbf{A}}(t)$, ${\mathbf{b}}(t)$ and $\mathbf{c}(t)$ having components that are T-periodic and in $\mathsf{Lip}(\alpha)$ with $\alpha>\frac{1}{2}$. If the Floquet-transformed state space matrix ${\mathbf{Q}}$ is Hurwitz, then $\left\|{\mathcal{G}}\right\|_{{{\mathsf{H}}_{2}}}<\infty$.

Proof.

Without loss of generality, assume ${\mathcal{G}}$ is in the Floquet form. Let

[TABLE]

denote the Fourier expansions of ${\mathbf{c}}(t)$ and ${\mathbf{b}}(t)$, respectively. Then from Sandberg (2006) [15] we know $\hat{g}_{k}(s)=\sum_{\ell\in{\mathbb{Z}}}{\mathbf{c}}_{k-\ell}^{T}[s_{\ell}{\mathbf{I}}-{\mathbf{Q}}]^{-1}{\mathbf{b}}_{\ell}$ with $s_{\ell}=s+{\dot{\imath\!\imath}}\ell\omega_{0}$. Additionally, $\left\|{\mathcal{G}}\right\|_{{{\mathsf{H}}_{2}}}^{2}=\displaystyle\sum_{k\in{\mathbb{Z}}}\left\|\hat{g}_{k}\right\|_{{{\mathsf{H}}_{2}}}^{2}$ from (20). First, show that $\|\hat{g}_{k}\|_{{{\mathsf{H}}_{2}}}<\infty$:

[TABLE]

where in the last step we used $\|[s_{\ell}{\mathbf{I}}-{\mathbf{Q}}]^{-1}\|_{{{\mathsf{H}}_{2}}}=\|[s{\mathbf{I}}-{\mathbf{Q}}]^{-1}\|_{{{\mathsf{H}}_{2}}}$. One can immediately observe that

[TABLE]

where “ * ” denotes the convolution operator. Note that $\{\|{\mathbf{c}}_{k}\|_{2}\},\{\|{\mathbf{b}}_{k}\|_{2}\}\in\ell_{1}$ by Lemma 3.4. Therefore, it follows from Lemma 3.3 that the sequence $\left\{\sum_{\ell\in{\mathbb{Z}}}\|{\mathbf{c}}_{k-\ell}\|_{2}\|{\mathbf{b}}_{\ell}\|_{2}\right\}$ is an $\ell_{1}$ sequence indexed in $k$. Since ${\mathbf{Q}}$ is Hurwitz, then $\|[s{\mathbf{I}}-{\mathbf{Q}}]^{-1}\|_{{{\mathsf{H}}_{2}}}<\infty$. Thus, $\{\|\hat{g}_{k}\|_{{{\mathsf{H}}_{2}}}\}$ forms an $\ell_{1}$ sequence indexed in $k$. Since $\ell_{1}\subset\ell_{2}$, $\|{\mathcal{G}}\|_{{{\mathsf{H}}_{2}}}^{2}=\displaystyle\sum_{k\in{\mathbb{Z}}}\|\hat{g}_{k}\|_{{{\mathsf{H}}_{2}}}^{2}<\infty$. ∎

3.2 Pole-Residue Representation of ${{\mathsf{H}}_{2}}$ Inner Product

An important result from LTI system theory is the representation of ${{\mathsf{H}}_{2}}$ inner products and ${{\mathsf{H}}_{2}}$ norms in terms of poles and residues as shown in (14). For SISO LTI systems, $G$ and $H$, that are stable, and assuming simple poles, $\{\lambda_{i}\}_{i=1}^{n}$ for ${\mathbf{G}}(s)$ and $\{\mu_{i}\}_{i=1}^{n}$ for for ${\mathbf{H}}(s)$, these inner product/norm representations simplify to:

[TABLE]

Even though this formulation of the inner product is rarely used for computation, it has been used in deriving optimality conditions for ${{\mathsf{H}}_{2}}$ approximation, see [1, 22]. In what follows, we extend this result to SISO LTP systems.

In the LTP setting, the meromorphic functions in question no longer have a finite number of poles. Although the pole residue formulation of the inner product provided above extends naturally in the way one would hope, the proof of this extension is considerably more difficult as the next theorem shows.

Theorem 3.6**.**

Let {\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{Q}}&\widetilde{{\mathbf{b}}}(t)\\ \hline\cr\widetilde{{\mathbf{c}}}^{T}(t)&\end{array}\right] and {\mathcal{H}}=\left[\begin{array}[]{c|c}{\mathbf{Q}}&{\mathbf{b}}(t)\\ \hline\cr{\mathbf{c}}^{T}(t)&\end{array}\right] denote two LTP systems sharing state space matrix ${\mathbf{Q}}$ with $\widetilde{{\mathbf{b}}}(t)$, $\widetilde{{\mathbf{c}}}^{T}(t)$, ${\mathbf{b}}(t)$, ${\mathbf{c}}^{T}(t)$ having components that are $T$-periodic and in $\mathsf{Lip}(\alpha)$ with $\alpha>\frac{1}{2}$. Let $\hat{g}_{k},\hat{h}_{k}\in{{\mathsf{H}}_{2}}$ be the $k$-th subsystems of ${\mathcal{G}}$ and ${\mathcal{H}}$ respectively where

[TABLE]

Assume that the eigenvalues $\{\lambda_{j}\}_{j=1}^{n}$ of ${\mathbf{Q}}$ are in the open left-half plane and that the fundamental frequency $\omega_{0}$ bounds the largest imaginary component of the spectrum of ${\mathbf{Q}}$. That is,

[TABLE]

If ${\mathbf{Q}}$ has simple poles, then $\hat{g}_{k}$ and $\hat{h}_{k}$ have poles at {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\lambda_{j}^{(\ell)}=\lambda_{j}({\mathbf{Q}})-{\dot{\imath\!\imath}}\ell\omega_{0},\ \ell\in{\mathbb{Z}}}} and the inner product, $\langle\hat{g}_{k},\hat{h}_{k}\rangle_{{{\mathsf{H}}_{2}}}$ can be written as

[TABLE]

The ${{\mathsf{H}}_{2}}$-norm of $\hat{g}_{k}(s)$ is given by

[TABLE]

Proof.

Fix $k$. We note that residue calculus immediately provides the pole residue formulation for subsystems with a finite number of poles. However, the subsystems $\hat{g}_{k}$ and $\hat{h}_{k}$ we consider have a countable number of poles that tesselate along the imaginary axis, $\lambda_{j}^{(\ell)}=\lambda_{j}({\mathbf{Q}})-{\dot{\imath\!\imath}}\ell\omega_{0}$ for $\ \ell\in{\mathbb{Z}}$.

First, consider the case where ${\mathbf{Q}}$ is scalar: $n=1$, $\mathbf{Q}=\lambda$ with $\textrm{Re}\{\lambda\}<0$. Without loss of generality we assume $\omega_{0}=1$. Then $\lambda^{(\ell)}=\lambda-\ell{\dot{\imath\!\imath}}$. Define $\phi_{\ell}=c_{k-\ell}b_{\ell}$ and $\widetilde{\phi}_{\ell}=\widetilde{c}_{k-\ell}\widetilde{b}_{\ell}$. Thus, we obtain

[TABLE]

with $\{\phi_{\ell}\},\{\widetilde{\phi}_{\ell}\}\in\ell_{1}$.

Consider the rectangular contour $\Gamma_{M}$ defined by the vertices $(0,M+\frac{1}{2})$, $(0,-M-\frac{1}{2})$, $(-M,-M-\frac{1}{2})$, $(-M,M+\frac{1}{2})$, for integer $M>|\lambda|$, so that $\Gamma_{M}$ contains $2M+1$ poles. Define $\Gamma_{M}^{(2)}$, $\Gamma_{M}^{(2)}$, $\Gamma_{M}^{(3)}$, $\gamma_{M}$ to represent the the top, left, bottom, and right contours of the rectangle respectively (see Figure 1). Finally, write $\Gamma_{M}=\Gamma_{M}^{(1)}+\Gamma_{M}^{(2)}+\Gamma_{M}^{(3)}+\gamma_{M}$ to represent the entire contour.

From residue calculus we have $\displaystyle\int_{\Gamma_{M}}\hat{g}_{k}(-s)\hat{h}_{k}(s)ds=\sum_{\ell=-M}^{M}\hat{g}_{k}(-\lambda^{(\ell)}){\textsf{res}}[\hat{h}_{k}(s),\lambda^{(\ell)}]$. We first show that $\int_{\Gamma_{M}^{(2)}}\hat{g}_{k}(-s)\hat{h}_{k}(s)ds\to 0$ as $M\to\infty$. Observe that $\Gamma_{M}^{(2)}=\{z=-M+{\dot{\imath\!\imath}}y:-M-\frac{1}{2}\leq y\leq M+\frac{1}{2}\}$. Then,

[TABLE]

We know that if $f(s)$ is continuous and bounded on the finite contour $\Gamma$, then

[TABLE]

Therefore, using (25), we can bound the contour integral over $\Gamma_{M}^{(2)}$ with $\textrm{length}(\Gamma_{M}^{(2)})=2M+1$,

[TABLE]

Next, we show that $\int_{\Gamma_{M}^{(1)}}\hat{g}_{k}(-s)\hat{h}_{k}(s)ds\to 0$ as $M\to\infty$:

[TABLE]

Then for $s\in\Gamma_{M}^{(1)}$,

[TABLE]

Note that $\mbox{length}(\Gamma_{M}^{(1)})=M.$ The inequality (25) yields

[TABLE]

We need to show that the quantity on the right approaches [math] as $M\to\infty$. To accomplish this, we perform a change of variable on the first summation,

[TABLE]

Then we can pass the limit,

[TABLE]

since M|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widetilde{\phi}_{\hat{\ell}-M}}}|=M|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widetilde{c}_{k-(\hat{\ell}-M)}\widetilde{b}_{\hat{\ell}-M}}}|\to 0 by Lemma 3.4.

We can employ a similar argument for $\Gamma_{N}^{(3)}$ to show that

[TABLE]

Then we obtain

[TABLE]

To consider the case where ${\mathbf{Q}}\in{\mathbb{R}}^{n\times n}$ with $n\geq 2$, we need to guarantee that the poles of $\hat{g}_{k}$ and $\hat{h}_{k}$, i.e., $\lambda_{j}^{(\ell)}=\lambda_{j}({\mathbf{Q}})-{\dot{\imath\!\imath}}\ell\omega_{0}$ for $\ell\in{\mathbb{Z}},$ will not coincide and consequently remain simple. The assumption that $\max_{j}|\textrm{Im}\{\lambda_{j}\}|<\omega_{0}$ is a sufficient condition to keep the poles simple and hence the case for finite dimension $n\geq 2$ is a natural generalization,

[TABLE]

The expression for the ${{\mathsf{H}}_{2}}$-norm follows immediately. ∎

Now we are in a position to give a result analogous to (21), but for LTP systems, representing the ${{\mathsf{H}}_{2}}$ inner product and norm using poles and residues.

Theorem 3.7**.**

Let ${\mathcal{G}}$ and ${\mathcal{H}}$ be two LTP systems as in Theorem 3.6 sharing the same state-matrix ${\mathbf{Q}}$ with $\widetilde{{\mathbf{b}}}(t)$, $\widetilde{{\mathbf{c}}}^{T}(t)$, ${\mathbf{b}}(t)$, ${\mathbf{c}}^{T}(t)$ having components that are $T$-periodic and in $\mathsf{Lip}(\alpha)$ with $\alpha>\frac{1}{2}$. Let $\hat{g}_{k},\hat{h}_{k}\in{{\mathsf{H}}_{2}}$ be the $k$th subsystems of ${\mathcal{G}}$ and ${\mathcal{H}}$ respectively. Assume $|Im\{\lambda_{j}({\mathbf{Q}})\}|<\omega_{0}$ where $\omega_{0}$ denotes the fundamental frequency of ${\mathcal{G}}$ and ${\mathcal{H}}$. Then,

[TABLE]

and

[TABLE]

where $\lambda_{j}^{(\ell)}=-\lambda_{j}+\ell{\dot{\imath\!\imath}}$ is as defined in Theorem 3.6.

Proof.

The result follows from combining Theorem 3.6 and Theorem 3.1. More specifically, combining (23) with (19) yields (26), and combining (24) with (20) yields (27). ∎

Theorem 3.7 extends the formulae (21)-(22) to the LTP setting. However, note that the inner product formula in the LTP setting assumes that the systems share the same ${\mathbf{Q}}$ matrix. Even though formulae are more complicated including two infinite sums, it is analogous to the LTI case in the sense it is a weighted sum of the residues where the weight is the subsystem transfer function evaluated at the mirror images of the poles. In the LTI case, this formulation has lead to the interpolatory ${{\mathsf{H}}_{2}}$ optimality conditions for model reduction [1]. This is a potential research direction to pursue for the LTP setting.

4 An ${{\mathsf{H}}_{2}}$-based Model Reduction Framework for LTP systems

Given the full-order LTP system ${\mathcal{G}}$ as in (3) with state-space dimension $n$, we seek a low-order LTP approximant $\widetilde{\mathcal{G}}$ as in (6) with state-space dimension $r$, where $r\ll n$. The mumerical implementation of the proposed method presented below assumes the availability of a Floquet-transformed equivalent system. Since we know such a transformation exists for every LTP system, without loss of generality, we assume ${\mathbf{A}}(t)={\mathbf{Q}}$ in this section. We note that in some prominent applications, e.g., in the analysis of mechanical systems with moving loads as considered in [33], the full-order system has already a constant state-mapping matrix ${\mathbf{A}}(t)={\mathbf{Q}}$. In other applications, e.g., nonlinear circuit modeling and analysis, one can compute the Floquet-transformation from the sampled simulation data as we did in our numerical example in Section 5.3.

Our approach uses a Petrov-Galerkin projection framework to perform reduction. Given the LTP system ${\mathcal{G}}$ as in (3) with now ${\mathbf{A}}(t)={\mathbf{Q}}$, we will construct two matrices ${\mathbf{V}},{\mathbf{W}}\in\mathbb{R}^{n\times r}$ with ${\mathbf{W}}^{T}{\mathbf{V}}=\mathbf{I}_{r}$ such that the reduced LTP system in (6) is given by $\widetilde{\mathbf{A}}(t)={\mathbf{W}}^{T}{\mathbf{Q}}{\mathbf{V}}$, $\widetilde{\mathbf{b}}(t)={\mathbf{W}}^{T}{\mathbf{b}}(t)$ and $\widetilde{\mathbf{c}}(t)={\mathbf{c}}(t){\mathbf{V}}$. Our proposed approach converts the LTP model reduction problem into an equivalent MIMO LTI problem, utilizing the optimal, numerically efficient $\mathsf{H}_{2}$ model reduction techniques for LTI systems to construct projection subspaces, ${\mathbf{V}}$ and ${\mathbf{W}}$ for the LTP problem. We are then able to provide an error bound for the approximation error in the original LTP setting.

Stykel and Vasilyev [33] introduced a model reduction scheme for linear time-varying (LTV) systems that was developed independently of the earlier thesis [34] upon which the present work is based. Stykel and Vasilyev considered a special class of LTV systems describing mechanical systems with moving loads; creating a time-dependent linear dynamical system equivalent in structure to a ${\mathbf{Q}}$-${\mathbf{b}}(t)$-${\mathbf{c}}(t)$ realization considered in Theorem 3.6 but for arbitrarily time-varying ${\mathbf{b}}(t)$ and ${\mathbf{c}}(t)$. The approach we pursue here originates with the presumed periodicity of ${\mathbf{b}}(t)$ and ${\mathbf{c}}(t)$, leading to a substantially different approach from that of [33]. We are able to take advantage of this presumed periodicity both analytically and computationally. Indeed, the extended-${{\mathsf{H}}_{2}}$ space of LTP systems introduced above allows for familiar performance guarantees entirely analogous to those in the LTI setting.

4.1 Connection of LTP Systems to LTI MIMO Systems

If ${\mathbf{b}}(t)$ and ${\mathbf{c}}(t)$ have finite order Fourier series expansions,

[TABLE]

then the LTP system has finite number of nontrivial subsystems, each of which have only finite spectra. This system has $2N+1$ inputs formed by the Fourier coefficients of ${\mathbf{b}}(t)$ and $2N+1$ outputs formed by the Fourier coefficients of ${\mathbf{c}}(t)$ resulting in $4N+1$ nontrivial subsystems. The next result connects the ${{\mathsf{H}}_{2}}$ norm of such an LTP system to the ${{\mathsf{H}}_{2}}$ norm of an LTI system.

Theorem 4.1**.**

Let \displaystyle{\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{Q}}&{\mathbf{b}}(t)\\ \hline\cr{\mathbf{c}}^{T}(t)&\mathbf{0}\end{array}\right] be an LTP system where ${\mathbf{b}}(t)$ and ${\mathbf{c}}(t)$ have the finite Fourier expansions as in (28). We define an associated LTI MIMO system

[TABLE]

Then the ${{\mathsf{H}}_{2}}$ norm of the LTP dynamical system ${\mathcal{G}}$ can be bounded in terms of the ${{\mathsf{H}}_{2}}$ norm of the LTI dynamical system ${{H}}$, namely

[TABLE]

Proof.

The above system, ${\mathcal{G}}$, has subsystems $\hat{g}_{k}(s)$ that can be written

[TABLE]

Each of these terms are finite-dimensional LTI systems. Using equation (20),

[TABLE]

Now we use that \|a_{1}+\ldots+a_{K}\|^{2}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{K}}(\|a_{1}\|^{2}+\ldots\|a_{K}\|^{2}) (application of the triangle inequality and the property, $2ab\leq a^{2}+b^{2}$ for $a,b\in{\mathbb{R}}$) to conclude

[TABLE]

Then we have

[TABLE]

where we recognize the last term as the ${{\mathsf{H}}_{2}}$ norm of the MIMO system,

[TABLE]

Thus, $\left\|{\mathcal{G}}\right\|^{2}_{{{\mathsf{H}}_{2}}}\leq(2N+1)\left\|{{H}}\right\|_{{{\mathsf{H}}_{2}}}^{2}$. ∎

This result leads to the model order reduction algorithm for LTP systems outlined in the next section.

4.2 A Model Reduction Method for a Special Case of LTP Systems

Inspired by Theorem 4.1, here we introduce a model-reduction scheme that works for a special case of LTP systems with time-invariant ${\mathbf{A}}(t)={\mathbf{Q}}$. We begin by introducing a performance guarantee of reduced-order LTP systems constructed from projection matrices ${\mathbf{V}}$ and ${\mathbf{W}}$. By Theorem 4.1, we immediately obtain the following corollary.

Corollary 4.2**.**

Given an LTP system of the form {\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{Q}}&{\mathbf{b}}(t)\\ \hline\cr{\mathbf{c}}^{T}(t)&\mathbf{0}\end{array}\right] where ${\mathbf{b}}(t)$ and ${\mathbf{c}}^{T}(t)$ have finite Fourier expansions as in (28). Let ${\mathbf{V}},{\mathbf{W}}\in{\mathbb{R}}^{n\times r}$ be projection matrices such that ${\mathbf{W}}^{T}{\mathbf{V}}={\mathbf{I}}$ and define the reduced LTP system

[TABLE]

Then

[TABLE]

where ${{H}}$ and $\widetilde{{H}}$ are the LTI MIMO counterparts, i.e.,

[TABLE]

Corollary 4.2 connects the model reduction error resulting from the analogous LTI MIMO problem to the original LTP model reduction problem. This gives a lot of flexibility since the tools for LTI model reduction are well established and one can employ various model reduction methods, such as Balanced Truncation [35, 36], optimal Hankel Norm Approximation [37], or Iterative Rational Krylov Algorithm (IRKA) [1]. However, since the error bound is given in terms of the ${{\mathsf{H}}_{2}}$ error and IRKA produces locally optimal approximations in the ${{\mathsf{H}}_{2}}$ sense, IRKA is a natural candidate for producing $\widetilde{{H}}$ that minimizes $\|{{H}}-\widetilde{{H}}\|_{{{\mathsf{H}}_{2}}}$.

Corollary 4.2 analyzed the case for finite Fourier expansion. The error analysis for the case with infinitely many Fourier coefficients is given next.

Corollary 4.3**.**

Given an LTP system of the form {\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{Q}}&{\mathbf{b}}(t)\\ \hline\cr{\mathbf{c}}^{T}(t)&\mathbf{0}\end{array}\right] with Fourier expansions ${\mathbf{b}}(t)=\sum_{-\infty}^{\infty}{\mathbf{b}}_{k}e^{{\dot{\imath\!\imath}}k\omega_{0}t}$ and ${\mathbf{c}}(t)=\sum_{-\infty}^{\infty}{\mathbf{c}}_{k}e^{{\dot{\imath\!\imath}}k\omega_{0}t}$, let ${\mathcal{G}}_{[N]}$ denote the Fourier truncated LTP system

[TABLE]

and the system ${{H}}_{[N]}$ is the associated LTI MIMO system for ${\mathcal{G}}_{[N]}$,

[TABLE]

*Moreover, let ${\mathbf{V}},{\mathbf{W}}\in{\mathbb{R}}^{n\times r}$ be projection matrices such that ${\mathbf{W}}^{T}{\mathbf{V}}={\mathbf{I}}$ and define the reduced system

---

[TABLE]

Then the dynamical system error $\|{\mathcal{G}}-\widetilde{\mathcal{G}}_{[N]}\|_{{{\mathsf{H}}_{2}}}$ is bounded by

[TABLE]

where $\widetilde{{H}}_{[N]}$ is the Petrov-Galerkin approximation of ${{H}}_{[N]}$, i.e.,

[TABLE]

A sketch of the proposed model reduction algorithm is given in Algorithm 1.

5 Numerical Results

We demonstrate the proposed model reduction scheme on three examples: (1) a heat equation with a moving point source, (2) a nonlinear transmission line LTP model, and (3) a constructed example, the structural model of component 1r (Russian service module) of the International Space Station. As most of the linear-time varying model reduction methods are computationally challenging for the problems we would like to consider, we compare our method to the Linear Time Varying Balanced Truncation method of [7] and Proper Orthogonal Decomposition (POD) [38], which remains computationally tractable even for large-scale problems since it only requires a time-domain simulation. Steih and Urban [39] introduce a space-time reduced basis method for time-periodic parametric partial differential equations. However, since we only consider nonparametric LTP systems here, for our purposes, it is enough to consider regular POD.

The proposed method as described in Algorithm 1 applies IRKA to a system with $2N+1$ inputs and outputs. As mentioned briefly in Section 2.2.3, the convergence of IRKA may slow down as the number of inputs and outputs increases. Thus, for modest $N$, one might expect that IRKA in Step 3 of Algorithm 1 may be slow to converge. The residue correction step introduced by Beattie and Gugercin in [23] has largely resolved this issue and improved the MIMO behavior of IRKA significantly; see [23]. For the problems studied here, we have found that the (regular) IRKA algorithm converged after a modest number of iterations even without the residue correction methodology of [23].

5.1 Computation of the ${{\mathsf{H}}_{2}}$ norm for error comparisons

Even though Algorithm 1 does not require computing the ${{\mathsf{H}}_{2}}$ norm at any point, in order to provide a detailed comparison to the reader between the reduced models resulting from Algorithm 1 and POD, we present the resulting ${{\mathsf{H}}_{2}}$ error norms. However, computing the ${{\mathsf{H}}_{2}}$ error of LTP systems is a nontrivial exercise and we discuss a practical implementation based on Zhou and Hagiwara, [19].

To compute the ${{\mathsf{H}}_{2}}$ norm of an LTP system {\mathcal{G}}=\left[\begin{array}[]{c|c}{\mathbf{Q}}&{\mathbf{b}}(t)\\ \hline\cr{\mathbf{c}}^{T}(t)&\end{array}\right], construct a new system, \mathfrak{G}_{N}=\left[\begin{array}[]{c|c}\mathfrak{A}_{N}&\mathfrak{B}_{N}\\ \hline\cr\mathfrak{C}_{N}&\end{array}\right] from a finite number of Floquet-Fourier coefficients $\{{\mathbf{b}}_{i}\}_{-N}^{N}$ and $\{{\mathbf{c}}_{i}\}_{-N}^{N}$ where

[TABLE]

Then we can approximate the LTP ${{\mathsf{H}}_{2}}$ norm $\|{\mathcal{G}}\|_{{{\mathsf{H}}_{2}}}$ with arbitrary accuracy as we keep more of the Floquet-Fourier coefficients.

Theorem 5.1** (Zhou & Hagiwara, [19]).**

If \mathfrak{G}_{N}=\left[\begin{array}[]{c|c}\mathfrak{A}_{N}&\mathfrak{B}_{N}\\ \hline\cr\mathfrak{C}_{N}&\end{array}\right] where $\mathfrak{A}_{N}$, $\mathfrak{B}_{N}$ and $\mathfrak{C}_{N}$ are defined according to (34), then

[TABLE]

where $\mathcal{V}_{N}$ and $\mathcal{W}_{N}$ are, respectively, the solutions of the finite-dimensional Lyapunov equations

[TABLE]

We note that in practice the Lyapunov equations (35) can be enormous in size as $\mathfrak{A}_{N}\in{\mathbb{R}}^{Nn\times Nn}$ where $N$ is the number of Fourier coefficients retained and $n$ is the order of the system. However, due to the block diagonal structure of $\mathfrak{A}_{N}$, the dominant cost of the Lyapunov solver is the Schur decomposition of ${\mathbf{Q}}$, which can be computed once and reused for ${\mathbf{Q}}-{\dot{\imath\!\imath}}k\omega_{0}{\mathbf{I}}$. Then the back substitution in the Bartels-Stewart algorithm can be done on each block individually, avoiding the cost of an $Nn\times Nn$ Schur decomposition. For description of the Bartels-Stewart algorithm, see Sorensen and Zhou [40]. We found this to be a critical step in making the ${{\mathsf{H}}_{2}}$ norm computations feasible for illustrating the model reduction errors.

In the three numerical examples below, to illustrate the formula (33), we plot the $\mathcal{H}_{2}$ error and its upper bound. However, as the discussion above illustrates, computing the exact Fourier truncation error $\|{\mathcal{G}}-{\mathcal{G}}_{[N]}\|_{{{\mathsf{H}}_{2}}}$ is not numerically feasible if ${\mathcal{G}}$ has an infinite Fourier expansion. Indeed, in the numerical example of Section 5.3, we do not even have full access to ${\mathcal{G}}$. Therefore, in Sections 5.2 and 5.3, we keep all the $N$ Fourier coefficient that we are able to obtain via a numerical simulation and treat this system approximately as ${\mathcal{G}}$, i.e., ${\mathcal{G}}\approx{\mathcal{G}}_{[N]}$. This means that in Sections 5.2 and 5.3, the $\mathcal{H}_{2}$ norm plots show the $\mathcal{H}_{2}$ error $\left\|{\mathcal{G}}_{[N]}-\widetilde{\mathcal{G}}_{[N]}\right\|_{{{\mathsf{H}}_{2}}}$ together with its upper bound $\sqrt{2N+1}\left\|{{H}}_{[N]}-\widetilde{{{H}}}_{[N]}\right\|_{{{\mathsf{H}}_{2}}}$. On the other, the numerical example of Section 5.4 is constructed such that the LTP system has only $N=5$ Fourier coefficients. Therefore, in that example we have ${\mathcal{G}}={\mathcal{G}}_{[N]}$ and the norm computations are exact.

5.2 1D Heat Model

The following model is a modified version of Example 1 in [7] of the 1D heat equation,

[TABLE]

with a moving point source where $\delta(t)$ is the Dirac delta function and $\xi(t)$ and $u(t)$ denote the heat source position and thermal flux respectively.

The PDE is discretized via a finite difference discretization with $2502$ equidistant grid points. Due to the Dirichlet boundary condition the discretized system has $n=2500$ degrees of freedom. Simulations are discretized in time via Backwards Euler with final time $T=100$ [sec] and $\Delta t=1$ [sec]. To make the source term periodic, we choose $\xi(t)=0.5+0.4\sin(8\pi t/T)$.

For this model we include a comparison to POD and the Linear Time-Varying Balanced Truncation (LTV BT) method introduced by Lang et al. [7] This approach solves the two Lyapunov equations at every time step and produces different ROM trial and test spaces (${\mathbf{V}}(t_{i})$ and ${\mathbf{W}}(t_{i})$) for every time step as well. While LTV BT produces very accurate reduced models, solving many Lyapunov equations is expensive. The authors alleviate this issue with an iterative process that warm-starts the Lyapunov solver at each time step with the solution from the previous time step; however, the computational cost of LTV BT is still considerably larger than our approach as is demonstrated in the numerical results that follow. It is important to note that our approach exploits periodicity of the state space representation, whereas LTV BT does not and indeed can handle a broader class of time-varying problems than we consider here.

We train POD with the constant input function $u(t)=1$. Both POD and H2-Based reduced models use $r=14$ for dimension of the reduced system. On the other hand, since the model reduction bases vary at every time step, LTV BT uses a varying reduced order throughout the simulation, $5\leq r\leq 9$, depending on the time step. The LTV BT simulations take considerably longer to run ($1355$ [sec]) whereas POD and the proposed $\mathcal{H}_{2}$-Based approach, i.e., Algorithm 1, each took less than five seconds. In Figure 2, we plot the output $y(t)$ (the top plot) and the output errors (the bottom plot) due to the three reduced models obtained via POD, LTV BT, and the proposed method (labeled as “H2-Based”). The input $u(t)$ for these simulations is the same input function that was used to train POD. The first observation is that the input/output based approaches outperform POD as illustrated by the error plot in Figure 2. For this example, LTV BT and the H2-based proposed approach perform similarly. However, one might expect LTV BT to outperform the proposed approach in general since it uses time varying model reduction bases ${\mathbf{V}}(t_{i})$ and ${\mathbf{W}}(t_{i})$ at every time step; as opposed to the proposed approach where the model reduction bases are fixed. Therefore, it is encouraging that the proposed method is able to mimic the accuracy of LTV BT for this example. Since LTV BT has time varying bases ${\mathbf{V}}(t_{i})$ and ${\mathbf{W}}(t_{i})$, the resulting reduced model does not allow $\mathcal{H}_{2}$ norm computations. Therefore, we compute the $\mathcal{H}_{2}$ error norms only for POD and the H2-based method. Results in Figure 3 show the ${{\mathsf{H}}_{2}}$ error as the dimension of the ROM increases. To make the ${{\mathsf{H}}_{2}}$ computation affordable, we have chosen a FOM with $n=100$ degrees of freedom. Once again the proposed method significantly outperforms POD.

5.3 Nonlinear Transmission Line

The following model is the “Test Network 2” from Noda et al. [41]. It represents an RLC circuit with nonlinear inductors to model the effect of saturable transformers. The network is split into 10 sections, each consisting of 6 states. Therefore, the full model has 60 states representing the current and voltages through each section.

In Figure 4, $e=E_{m}\cos(\omega_{0}t)$ where $\omega_{0}=2\pi 60$ [rad/sec] and $E_{m}=\sqrt{2/3}\times V_{p.u.}\times 500$ [kV]. The circuit is supplied with an “overvoltage” by a factor of $V_{p.u.}=1.25$ to induce a noticeable saturation in the nonlinear inductors, $L_{N}$. The current through $L_{N}$ is $i_{N}=\alpha\psi+\beta\psi^{7}$ where $\psi$ is the magnetic flux through the inductor and $\alpha$ and $\beta$ depend on the section. For a table of the resistor, capacitor, and inductor values, see Figure 8 in Noda et al. [41]. See Figure 4 for the nonlinear circuit schematic.

The steady-state solution to the problem is time-periodic. Noda et al. [41] are concerned with the transient deviations from the steady-state solution from perturbations of the input $e\mapsto e+u$. The voltage $e$ and current $i_{0}$ are the input and output of the dynamical system. As both quantities are degrees of freedom in the dynamical system, the input matrix ${\mathbf{b}}\in{\mathbb{R}}^{60}$ and output matrix ${\mathbf{c}}\in{\mathbb{R}}^{60}$ do not vary with time. Therefore, a linearization of the nonlinear dynamical system results in a time-periodic system of the form (3) where ${\mathbf{A}}(t),{\mathbf{b}},{\mathbf{c}}$ are $T$-periodic with $T=1/60$ [sec] and ${\mathbf{x}}(t)\in{\mathbb{R}}^{60}$ represents the perturbation in the state from the steady-state solution and $u(t)\in{\mathbb{R}}$ represents the perturbation from $e$.

Unlike the previous 1D Heat model, the full LTP models in this example and the next are stiff, requiring many time-samples to resolve the numerical simulation, even with the implicit backward Euler scheme. This results in significantly more Lyapunov equations to solve increasing the cost of LTV BT even further compared to the previous example. Therefore, for these last two examples, we only provide comparison to POD. However, even though LTV BT might be computationally more intensive, we still accept it to provide very accurate reduced models due to the time-varying model reduction bases. To run POD, we generated simulations of the system when the control was the Heaviside function, $u(t)=H(10^{-3}-t)$ [kV]. After generating the numerical simulation of the state ${\mathbf{x}}(t)$, we collect the simulation into a snapshot matrix, ${\mathbf{X}}=[{\mathbf{x}}_{t_{1}},\ldots,{\mathbf{x}}_{t_{m}}]$ and perform an SVD to construct the subspace ${\mathbf{V}}\in{\mathbb{R}}^{n\times r}$. The result of the reduced-order model generated according to (6) with ${\mathbf{W}}={\mathbf{V}}$ from above is labeled as POD in Figure 5.

For the model reduction via Algorithm 1, we begin by performing a Floquet transformation on the periodic ${\mathbf{A}}(t)$. For this example, we have the snapshots ${\mathbf{A}}(t_{i})$ for 256 uniformly spaced temporal points in the period, $t_{i}\in[0,1/60]$ [sec]. Therefore, the Floquet-Fourier transformation resulted in a state-space matrices of the form: ${\mathbf{Q}}\in{\mathbb{R}}^{60\times 60}$, ${\mathbf{B}},\mathbf{C}\in{\mathbb{R}}^{60\times 256}$. We preserve all 256 Fourier modes of the input and output and used IRKA to construct the Petrov-Galerkin subspaces ${\mathbf{W}}$ and ${\mathbf{V}}$.

Figure 5 shows the outputs (the top plot) and the output errors (the bottom plot) in time domain simulations due to two reduced models obtained via POD and the proposed method (labeled as “H2-Based”), each with reduced order $r=10$. We also perform a similar comparison for a variety of reduced order sizes in Figure 6, measuring the ${{\mathsf{H}}_{2}}$ error of the LTP systems. Note that the POD-based reduced model becomes unstable for $r\geq 12$; hence, the ${{\mathsf{H}}_{2}}$ error is infinite in these cases and are not plotted.

As Figures 5 and 6 illustrate the proposed approach demonstrates superior performance to POD, in terms of stability and dynamical system error measured in the ${{\mathsf{H}}_{2}}$ norm for LTP systems. In terms of both the time domain simulation and the ${{\mathsf{H}}_{2}}$ error norm, the error due to the proposed method is almost two orders of magnitude smaller than that of POD. We note that this superior performance is obtained in the best case scenario for POD since the model reduction error is measured for the same input which was used to train POD. Still the ${{\mathsf{H}}_{2}}$-based proposed method produces a significantly better reduced model. We note also that the error bound predicts the true error behavior well.

5.4 Structural Model of Component 1r (Russian service module) of the International Space Station

We consider the 1r Russian service module that has $270$ states, $3$ inputs and $3$ outputs. We make this model LTP by placing modulators on inputs $2$ and $3$. We do the same for outputs $2$ and $3$. Define \mathbf{B}=\left[\begin{array}[]{ccc}{\mathbf{b}}_{0}&{\mathbf{b}}_{1}&{\mathbf{b}}_{2}\end{array}\right] to be the original input vector and $\mathbf{C}=[{\mathbf{c}}_{0},{\mathbf{c}}_{1},{\mathbf{c}}_{2}]^{T}$ to be the original output vector. Then we construct a SISO system by feeding the input $u(t)$ into two modulators with local oscillator frequencies of $\omega_{0}$ and $2\omega_{0}$ respectively,

[TABLE]

To see this abstracted into a diagram see Figure 7.

The system ${\mathcal{G}}$ is already in the Floquet-Fourier form of dimension $n=270$ and $N=2$ since the Fourier expansions of ${\mathbf{b}}(t)$ and ${\mathbf{c}}(t)$ have 5 nontrivial terms. We tried to reduce the model by POD, using a test function of $u(t)=\sin(19.2875t)$ but each reduced model was unstable for reduced orders greater than or equal to 5. The frequency $\omega=19.2875$ [rad/sec] was chosen after determining that it excited many of the system harmonics. Since the resulting reduced systems were unstable, we omitted the POD simulation results from these comparisons.

For $r=30$, we run simulations of the full model and reduced model obtained via Algorithm 1. As shown in Figure 8, the reduced LTP model is almost indistinguishable from the full LTP model.

As in the previous example, we construct reduced order models of dimensions $r=4$ to $r=80$ using Algorithm 1 and compare their respective error system ${{\mathsf{H}}_{2}}$ error norm $\|{\mathcal{G}}-\widetilde{\mathcal{G}}\|_{{{\mathsf{H}}_{2}}}$ in Figure 9. The error bound accurately predicts the true error.

6 Conclusions

We develop a model reduction scheme for LTP systems by converting the model reduction problem into an analogous LTI problem and employing existing model reduction techniques to the new problem. Numerical results demonstrate the success of the proposed method. Moreover, we extend the analysis of the certain notions of ${{\mathsf{H}}_{2}}$-approximation established for LTI systems to LTI systems.

In practice, current approaches for computing Floquet transformations do not scale well to large-scale systems, so this aspect may remain a bottleneck for effective model reduction for general LTP systems. In a variety of circumstances, however, this step is not difficult or problematic. We do not consider this aspect of the problem in the present work and the algorithm we propose here truncates the Floquet-Fourier coefficients, keeping $2N+1$ centered coefficients. The error introduced by this step becomes arbitrarily small as we increase the number of coefficients conserved, $N\to\infty$. Note that for any given $N$, the coefficients that are kept need not be an optimal choice and investigating what may constitute a better selection of coefficients would benefit any further refinement of this approach.

7 Acknowledgements

We thank to Dr. Taku Noda and Dr. Jens Saak for generously providing us with the MATLAB code used to generate the results in their paper [41] and [7], respectively. The work of C. Magruder was supported in part by the ExxonMobil Ken Kennedy Institute 2015/2016 High Performance Computing Graduate Fellowship. The work of S. Gugercin was supported in part by NSF through Grants DMS-1217156 and DMS-1522616, and by the Alexander von Humboldt Foundation. The work of C. Beattie was supported in part by NSF through Grant DMS-1217156 and by the Einstein Foundation - Berlin.