# Computation of lower bounds
for the induced $\mathcal{L}_{2}$ norm of LPV systems

Tamas Peni, and Peter J. Seiler T. Peni is with the Systems and Control
Laboratory of Institute for Computer Science and Control, 1111 Budapest, Kende
u. 13-17., Hungary. E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>J.
Seiler is with Aerospace Engineering and Mechanics, University of Minnesota,
107 Akerman Hall, 110 Union St SE, Minneapolis, MN 55455-0153. E-mail:
<EMAIL_ADDRESS>
###### Abstract

Determining the induced $\mathcal{L}_{2}$ norm of a linear, parameter-varying
(LPV) system is an integral part of many analysis and robust control design
procedures. Most prior work has focused on efficiently computing upper bounds
for the induced $\mathcal{L}_{2}$ norm. The conditions for upper bounds are
typically based on scaled small-gain theorems with dynamic multipliers or
dissipation inequalities with parameter dependent Lyapunov functions. This
paper presents a complementary algorithm to compute lower bounds for the
induced $\mathcal{L}_{2}$ norm. The proposed approach computes a lower bound
on the gain by restricting the parameter trajectory to be a periodic signal.
This restriction enables the use of recent results for exact calculation of
the $\mathcal{L}_{2}$ norm for a periodic linear time varying system. The
proposed lower bound algorithm has two benefits. First, the lower bound
complements standard upper bound techniques. Specifically, a small gap between
the bounds indicates that further computation, e.g. upper bounds with more
complex Lyapunov functions, is unnecessary. Second, the lower bound algorithm
returns a ”bad” parameter trajectory for the LPV system that can be further
analyzed to provide insight into the system performance.

###### Index Terms:

Linear parameter varying systems, induced $\mathcal{L}_{2}$ norm, periodic,
linear time-varying systems

##  I Introduction

Determining the induced $\mathcal{L}_{2}$ norm of an LPV system is an integral
part of many analysis and robust control design procedures. In general, this
norm cannot be determined explicitly. Most prior work focuses on computing an
upper bound on the induced $\mathcal{L}_{2}$ norm. The method used to compute
the induced gain upper bound depends primarily on the structure of the LPV
system. One class of LPV systems involves state matrices that are rational
functions of the parameter. In this case the LPV system can be expressed as a
feedback interconnection (linear fractional transformation, LFT) of a known
linear time invariant (LTI) system and a parameter dependent block. The upper
bound for these LFT-type LPV systems is typically computed using scaled mall
gain theorems with multipliers and the full block S-procedure [1, 2, 3].
Another class of LPV systems involves state matrices with an arbitrary
dependence on the parameter. These systems, called ”Gridded” LPV systems, can
arise via linearization of nonlinear models on a grid of equilibrium operating
points [4, 5, 6, 7, 8]. The upper bound for Gridded LPV systems is computed
using the dissipation inequality evaluated over a finite set of parameter grid
points. In both cases the upper bound calculation reduces to a convex
optimization with linear matrix inequality (LMI) constraints.

This paper addresses the complementary problem, i.e. the approximation for the
lower bound of the induced $\mathcal{L}_{2}$ norm. Frozen point analysis is
one simple way to compute a lower bound. Specifically, the induced
$\mathcal{L}_{2}$ norm of the LPV system over all possible parameter
trajectories is clearly lower bounded by the restriction to constant (frozen)
parameter trajectories. A lower bound is thus obtained by computing the
maximum gain of the LTI systems obtained by evaluating the LPV system on a
grid of frozen parameter values. To our knowledge this is currently the only
method available to compute a lower bound on the induced $\mathcal{L}_{2}$
norm of an LPV system. Unfortunately, this approach produces very conservative
results in many cases as it neglects the variation of the scheduling
parameter. To compute a better estimate, the $\mathcal{L}_{2}$ norm has to be
evaluated over time varying parameter trajectories. This concept introduces a
complex optimization problem involving a maximization over both the allowable
parameter trajectories and the $\mathcal{L}_{2}$ inputs to the LPV system.

The proposed lower bound algorithm restricts the scheduling trajectories to
periodic signals. This restriction is useful because the LPV system evaluated
on a fixed periodic trajectory is simply a periodic linear time-varying (PLTV)
system. Moreover, the induced $\mathcal{L}_{2}$ norm of a PLTV system can be
exactly determined by using recently developed numerical methods, see e.g.
[9]. Thus the lower bound algorithm only requires the maximization over
periodic parameter trajectories and the maximization over $\mathcal{L}_{2}$
inputs is not required. Specifically, the lower bound algorithm maximizes a
cost function related to the induced $\mathcal{L}_{2}$ norm (described in
detail in Sec. III) over a finite dimensional subspace of periodic scheduling
trajectories. The output of the algorithm is the lower bound estimate and the
worst-case (bad) parameter trajectory. In addition, we construct a worst-case
input signal that approximately achieves the computed induced gain lower bound
for the LPV system. The construction of the worst-case input is based on the
PLTV results collected in [10] and [11]. In [10] the worst case input is
derived by using a special, frequency domain representation for the PLTV
system. In this paper the worst-case input construction and proofs are
provided using only time-domain operators. This approach, detailed in section
II-C, streamlines the construction of the worst-case input.

It is important to note that our algorithm does not assume any specific
structure on the LPV system. Thus it is developed for the ”Gridded” case. It
trivially applies for ”LFT” models as well although the additional rational
structure in LFT models may lead to faster lower bound algorithms than those
developed here. It will also be shown, by numerical examples, that the
proposed lower bound algorithm and the known upper bound methods complement
each other. Specifically, these methods, applied together, can yield a tight
bound for the induced $\mathcal{L}_{2}$ norm of an LPV system.

The paper is organised as follows: in the next section the properties of
periodic systems are summarized and the recent results related to the induced
$\mathcal{L}_{2}$ norm computation are reviewed. In Sec. II-C we provide an
algorithm to construct the worst-case input for PLTV systems. Our approach for
computing lower bounds on $\mathcal{L}_{2}$ norm of LPV systems is presented
in Sec. III and III-C. Sec. IV is devoted to the numerical simulations and
analysis. The conclusions are drawn in section V, where the future works
required to improve our results are also summarized.

The notations used in the paper are fairly standard. $\mathbb{R}$ and
$\mathbb{C}$ denote the fields of real and complex numbers, respectively.
$\mathbb{D}$ and $\partial\mathbb{D}$ are the unit disc and unit circle in the
complex plane $\mathbb{C}$. If $\mathbf{T}:E\rightarrow F$ is a bounded,
linear operator between Hilbert spaces $E$ and $F$, then the adjoint of
$\mathbf{T}$ is denoted by $\mathbf{T}^{*}$ and is defined as follows:
$\langle\mathbf{T}u,y\rangle_{F}=\langle u,\mathbf{T}^{*}y\rangle_{E}$,
$\forall u\in E$ and $\forall y\in F$. The space of square-integrable signals
$f:[0,\infty)\rightarrow\mathbb{R}^{n}$ is denoted by
$\mathcal{L}_{2}(\mathbb{R}^{n})$. The inner product and norm in
$\mathcal{L}_{2}(\mathbb{R}^{n})$ are defined as $\langle
f,g\rangle_{\mathcal{L}_{2}(\mathbb{R}^{n})}=\int_{0}^{\infty}f(t)^{T}g(t)dt$
and $\|f\|=\sqrt{\langle f,f\rangle_{\mathcal{L}_{2}(\mathbb{R}^{n})}}$,
respectively. If $\mathbf{G}$ is a bounded linear operator, such that
$\mathbf{G}:w\in\mathcal{L}_{2}(\mathbb{R}^{p})\mapsto
z\in\mathcal{L}_{2}(\mathbb{R}^{q})$ then its induced $\mathcal{L}_{2}$ norm
is defined as

$\displaystyle\|\mathbf{G}\|:=\sup_{0\neq
w\in\mathcal{L}_{2}(\mathbb{R}^{p})}\frac{\|z\|}{\|w\|}$

$\mathcal{L}_{2,[0,h)}(\mathbb{R}^{n})$ denotes the space of square-integrable
functions on the interval $[0,h)$ with inner product $\langle
f,g\rangle_{\mathcal{L}_{2,[0,h)}(\mathbb{R}^{n})}=\int_{0}^{h}f(t)^{T}g(t)dt$.
$\ell_{2}(E)$ denotes the square-summable sequences
$w=\\{w_{k}\\}_{k=0}^{\infty}$ in the Hilbert space $E$ with inner product
$\langle w,v\rangle_{\ell_{2}(E)}:=\sum_{k=0}^{\infty}\langle
w_{k},v_{k}\rangle_{E}$.

##  II Periodic systems

###  II-A Background

This section reviews known results on PLTV systems. Most results can be found
in [11]. In particular, we consider the linear time-varying system $G$ of the
form:

$\displaystyle
G:~{}~{}\\!\begin{array}[]{rcl}\dot{x}(t)&=&A(t)x(t)+B(t)w(t)\\\
z(t)&=&C(t)x(t)+D(t)w(t)\end{array}\\!$ (3)

that has the following properties:

  1. (i)

the system matrices $A(t):\mathbb{R}\rightarrow\mathbb{R}^{n\times n}$,
$B(t):\mathbb{R}\rightarrow\mathbb{R}^{n\times p}$,
$C(t):\mathbb{R}\rightarrow\mathbb{R}^{q\times n}$,
$D(t):\mathbb{R}\rightarrow\mathbb{R}^{q\times p}$ are bounded continuous
functions of time,

  2. (ii)

are $h$-periodic, i.e. $A(t+h)=A(t)$, $B(t+h)=B(t)$, $C(t+h)=C(t)$,
$D(t+h)=D(t)$, $\forall t$ and

  3. (iii)

the dynamics (3) are internally stable.

The state transition matrix $\Phi(t,\tau)$ associated with the autonomous
dynamics $\dot{x}(t)=A(t)x(t)$ is defined for all $(t,\tau)$ as the linear
mapping from $x(\tau)$ to $x(t)$, i.e. $x(t)=\Phi(t,\tau)x(\tau)$. It is easy
to check that $\Phi(t,\tau)$ is periodic with $h$, i.e.
$\Phi(t+h,\tau+h)=\Phi(t,\tau)$, and satisfies the matrix differential
equation $\frac{d}{dt}{\Phi}(t,\tau)=A(t)\Phi(t,\tau)$, $\Phi(\tau,\tau)=I$.
The state transition matrix over one period is called the monodromy matrix:
$\Psi(t):=\Phi(t+h,t)$. Clearly, the monodromy matrix is periodic with $h$,
i.e. $\Psi(t+h)=\Psi(t)$.

The system defined in (3) is stable and hence it defines a bounded operator
$\mathbf{G}:=w\in\mathcal{L}_{2}(\mathbb{R}^{p})\mapsto
z\in\mathcal{L}_{2}(\mathbb{R}^{q})$, such that

$\displaystyle
z(t)=\int_{0}^{t}C(t)\Phi(t,\tau)B(\tau)w(\tau)~{}d\tau+D(t)w(t)$ (4)

i.e. $z(t)$ is the response of (3) to the input signal $w(t)$ if $x(0)=0$. We
are interested in $\mathbf{G}$. It is shown in [12] that $\|\mathbf{G}\|$ is
equal to the norm of the ”lifted” operator
$\hat{\mathbf{G}}:\ell_{2}(\mathcal{L}_{2,[0,h)}(\mathbb{R}^{p}))\rightarrow\ell_{2}(\mathcal{L}_{2,[0,h)}(\mathbb{R}^{q}))$,
which, has a finite dimensional state-space realization in the following form
[13]:

$\displaystyle\xi_{k+1}$ $\displaystyle=\hat{A}\xi_{k}+\hat{B}\hat{w}_{k}$
$\displaystyle\hat{z}_{k}$ $\displaystyle=\hat{C}\xi_{k}+\hat{D}\hat{w}_{k}$
(5)

where $\xi_{k}\in\mathbb{R}^{n}$ and
$\hat{A}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$,
$\hat{B}:\mathcal{L}_{2,[0,h)}(\mathbb{R}^{p})\rightarrow\mathbb{R}^{n}$,
$\hat{C}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{q}$,
$\hat{D}:\mathcal{L}_{2,[0,h)}(\mathbb{R}^{p})\rightarrow\mathcal{L}_{2,[0,h)}(\mathbb{R}^{q})$
are linear operators defined as follows:

$\displaystyle\hat{A}\xi_{k}$ $\displaystyle:=\Psi(0)\xi_{k},$
$\displaystyle\hat{B}\hat{w}_{k}:=\int_{0}^{h}\Phi(h,\tau)B(\tau)\hat{w}_{k}(\tau)d\tau$
$\displaystyle\hat{C}\xi_{k}$ $\displaystyle:=C(t)\Phi(t,0)\xi_{k},$
$\displaystyle\hat{D}\hat{w}_{k}:=\int_{0}^{t}C(t)\Phi(t,\tau)\hat{w}_{k}(\tau)d\tau+D(t)\hat{w}_{k}(t)$
(6)

Although (5) is a finite dimensional system, it is still not suitable for
numerical computations, because its system matrices are operators. Therefore,
based on the same idea applied for LTI systems (see e.g. Chapter 21 in [14]),
a discrete-time, linear, time-invariant system $\underline{G}_{\gamma}$ is
introduced in [15],[9]. The state-space matrices of $\underline{G}_{\gamma}$
are computed as follows:

$\displaystyle\underline{A}_{\gamma}$
$\displaystyle:=\hat{A}+\hat{B}(\gamma^{2}I-\hat{D}^{*}\hat{D})^{-1}\hat{D}^{*}\hat{C}$
$\displaystyle\underline{B}_{\gamma}\underline{B}_{\gamma}^{*}$
$\displaystyle:=\gamma\hat{B}(\gamma^{2}I-\hat{D}^{*}\hat{D})^{-1}\hat{B}^{*}$
(7) $\displaystyle\underline{C}_{\gamma}^{*}\underline{C}_{\gamma}$
$\displaystyle:=\gamma\hat{C}^{*}(\gamma^{2}I-\hat{D}\hat{D}^{*})^{-1}\hat{C}$

More specifically, the matrices $\underline{B}_{\gamma}$ and
$\underline{C}_{\gamma}$ are defined to be full rank matrices that satisfy the
equalities that appear in (7). Note that
$\underline{A}_{\gamma},\underline{B}_{\gamma},\underline{C}_{\gamma}$ are now
real matrices $\underline{A}_{\gamma}\in\mathbb{R}^{n\times
n},\underline{B}_{\gamma}\in\mathbb{R}^{n\times
p},\underline{C}_{\gamma}\in\mathbb{R}^{q\times n}$. The following theorem,
taken from [9], proves that there is a strong relation between the induced
$\ell_{2}$-norm of $\underline{G}_{\gamma}$ and the norm of
$\hat{\mathbf{G}}$.

######  Theorem 1.

The following statements are equivalent:

  1. a)

$\text{eig}(\hat{A})\in\mathbb{D}\backslash\mathbb{\partial D}$ and
$\|\hat{\mathbf{G}}\|<\gamma$

  2. b)

$\text{eig}(\underline{A}_{\gamma})\in\mathbb{D}\backslash\mathbb{\partial D}$
and $\|\underline{G}_{\gamma}\|_{\infty}<1$

where $\|\underline{G}_{\gamma}\|_{\infty}$ is the standard
$\mathcal{H}_{\infty}$ norm of the LTI system $\underline{G}_{\gamma}$.

###  II-B Induced norm computation

Theorem 1 gives the base of the bisection algorithm, proposed in [9] to
compute the $\mathcal{L}_{2}$-norm of the periodic system (3). Therefore, an
efficient method is needed to determine the system matrices
$\underline{A}_{\gamma},\underline{B}_{\gamma},\underline{C}_{\gamma}$. For
this, consider the following differential equation associated with (3):

$\displaystyle\dot{e}(t)=H(t)e(t)$ (8)

where (suppressing the notation of time dependence)

$\displaystyle H_{11}$
$\displaystyle=-A^{T}-C^{T}D(\gamma^{2}I-D^{T}D)^{-1}B^{T}$ $\displaystyle
H_{12}$ $\displaystyle=-\gamma C^{T}(\gamma^{2}I-DD^{T})^{-1}C$ $\displaystyle
H_{21}$ $\displaystyle=\gamma B(\gamma^{2}I-D^{T}D)B^{T}$ (9) $\displaystyle
H_{22}$ $\displaystyle=A+B(\gamma^{2}I-D^{T}D)^{-1}D^{T}C$

Note that (8) is a Hamiltonian system ($JH(t)$ is symmetric with
$J=\left[\begin{smallmatrix}0&I\\\ -I&0\end{smallmatrix}\right]$). Let the
monodromy matrix of (8) be denoted by $\Psi_{H}(t)$ and let $Q=\Psi_{H}(0)$.
It was shown in [9], that

$\displaystyle\underline{B}_{\gamma}\underline{B}_{\gamma}^{T}$
$\displaystyle=Q_{21}Q_{11}^{-1}$
$\displaystyle\underline{C}_{\gamma}^{T}\underline{C}_{\gamma}$
$\displaystyle=-Q_{11}^{-1}Q_{12}$ $\displaystyle\underline{A}_{\gamma}$
$\displaystyle=Q_{22}-Q_{21}Q_{11}^{-1}Q_{12}$ (10)

where $Q=\left[\begin{smallmatrix}Q_{11}&Q_{12}\\\
Q_{21}&Q_{22}\end{smallmatrix}\right]$. Since $Q$ is symplectic [16], i.e.
$Q^{T}JQ=J$, thus $Q$ is uniquely determined by (10) in the following form:

$\displaystyle
Q=\left[\\!\begin{array}[]{cc}\underline{A}_{\gamma}^{-T}&-\underline{A}_{\gamma}^{-T}\underline{C}_{\gamma}^{T}\underline{C}_{\gamma}\\\
\underline{B}_{\gamma}\underline{B}_{\gamma}^{T}\underline{A}_{\gamma}^{-T}&\underline{A}_{\gamma}-\underline{B}_{\gamma}\underline{B}_{\gamma}^{T}\underline{A}_{\gamma}^{-T}\underline{C}_{\gamma}^{T}\underline{C}_{\gamma}\end{array}\\!\right]$
(13)

One possible method to determine the
$\underline{A}_{\gamma},\underline{B}_{\gamma},\underline{C}_{\gamma}$
matrices can be given as follows. Integrate first the Hamiltonian system (8)
on $[0,h)$ starting from the matrix initial condition $e(0)=I$. Then $e(h)=Q$.
Determine the system matrices from $Q$ by using (10). Since $H(t)$ is not
stable this approach is numerically not reliable. The method proposed in [9]
is based on the following relations:

$\displaystyle\underline{A}_{\gamma}=X(h),~{}~{}~{}\underline{C}_{\gamma}^{T}\underline{C}_{\gamma}=Z(0),~{}~{}~{}\underline{B}_{\gamma}\underline{B}_{\gamma}^{T}=Y(h)$

where $X(h),Z(0),Y(h)$ are point solutions of the differential Riccati
equations

$\displaystyle\dot{Z}$ $\displaystyle=-H_{22}^{T}Z-ZH_{22}-\gamma
ZB(\gamma^{2}I-D^{T}D)^{-1}B^{T}Z+\gamma C^{T}(\gamma^{2}I-DD^{T})^{-1}C$ (14)
$\displaystyle\dot{X}$ $\displaystyle=(H_{22}+\gamma
B(\gamma^{2}I-D^{T}D)^{-1}B^{T}Z)X$ (15) $\displaystyle\dot{Y}$
$\displaystyle=H_{22}Y+YH_{22}^{T}+\gamma
YC^{T}(\gamma^{2}I-DD^{T})^{-1}CY+\gamma B(\gamma^{2}I-D^{T}D)^{-1}B^{T}$ (16)

with boundary conditions $Z(h)=0,X(0)=I$ and $Y(0)=0$. Integrating these
Riccati equations is a numerically better conditioned problem than the direct
integration of (8) [9]. Moreover, (14) has a very useful property [9]:

######  Lemma 1.

Assume $\gamma^{2}I-D(t)^{T}D(t)>0$ for all $t\in[0,h)$. Then (14) has a
bounded solution over the interval $[0,h]$ if and only if
$\|\hat{D}\|<\gamma$.

The bisection algorithm proposed in [9] is based on iteratively solving
(14)-(16) and computing $\|\underline{G}_{\gamma}\|_{\infty}$ at different
$\gamma$ values that are tuned in a bisection loop. The output is a lower- and
upper bound pair $(\underline{\gamma},\overline{\gamma})$ satisfying
$\underline{\gamma}\leq\|\mathbf{G}\|\leq\overline{\gamma}$ such that
$\overline{\gamma}-\underline{\gamma}\leq\varepsilon$, where $\varepsilon$ is
a given tolerance.

######  Remark 1.

If we introduce $T=\left[\begin{smallmatrix}0&I\\\ \gamma
I&0\end{smallmatrix}\right]$ and apply the state transformation
$\tilde{e}(t)=T^{-1}e(t)$ in (8), we get the following transformed Hamiltonian
system:

$\displaystyle\dot{\tilde{e}}(t)=\tilde{H}\tilde{e}(t),~{}~{}\text{ with
}~{}~{}\tilde{H}(t)=T^{-1}H(t)T.$ (17)

This system is the same as that is used in [10]. Denote $\Psi_{\tilde{H}}(t)$
the monodromy matrix associated with (17) and let
$\tilde{Q}:=\Psi_{\tilde{H}}(0)$. Then, it can be shown that
$\tilde{Q}=T^{-1}QT$. We use (17) instead of (8) in the next section, because
(17) is more convenient for the forthcoming derivations.

###  II-C Construction of the worst-case input

We are also interested in constructing the worst-case input
$w^{\circ}\in\mathcal{L}_{2}$ that achieves the induced norm $\|\mathbf{G}\|$.
The signal achieving this gain is not in $\mathcal{L}_{2}$ and hence our
construction approximately achieves this gain. Specifically, for any
$\epsilon>0$ we construct an input $w^{\circ}\in\mathcal{L}_{2}$ such that
$z^{\circ}:=\mathbf{G}w^{\circ}$ satisfies
$\|z^{\circ}\|\geq\|\mathbf{G}\|\|w^{\circ}\|-\epsilon$. This is similar to
the LTI case where the worst-case input is a sinusoid and a truncated sinusoid
approximately achieves the system gain. The method used to construct the
worst-case input for the PLTV system is based on the proof of Lemma 2.6 in
[10]. In [10] the worst-case input is constructed using a special frequency-
domain representation of $\mathbf{G}$ defined over exponentially periodic
signals. In this section, an alternative construction and proof is provided
using only time-domain formulations. This alternative proof streamlines the
numerical construction of the worst-case input.

To this end, let the linear operator
$\mathbf{T}_{G}:\mathbb{R}^{n}\times\mathcal{L}_{2,[0,h)}(\mathbb{R}^{p})\rightarrow\mathbb{R}^{n}\times\mathcal{L}_{2,[0,h)}(\mathbb{R}^{q})$
be defined as follows:

$\displaystyle\mathbf{T}_{G}(x_{0},w)\rightarrow(x_{h},z):=(\hat{A}x_{0}+\hat{B}w,~{}\hat{C}x_{0}+\hat{D}w).$

$\mathbf{T}_{G}$ is equivalent to the (state-input) $\rightarrow$ (state-
output) map realized by (3) over the period $[0,h)$:

$\displaystyle\dot{x}(t)$ $\displaystyle=A(t)x(t)+B(t)w(t),~{}~{}x(0)=x_{0}$
$\displaystyle z(t)$ $\displaystyle=C(t)x(t)+D(t)w(t),~{}~{}x_{h}=x(h)$ (18)

If we introduce an inner product in
$\mathbb{R}^{n}\times\mathcal{L}_{2,[0,h)}(\mathbb{R}^{\cdot})$ as

$\displaystyle\langle(x,w),(y,v)\rangle=x^{*}y+\int_{0}^{h}w^{*}(t)v(t)dt,$

then we can define the adjoint operator $\mathbf{T}_{G}^{*}$ to satisfy the
equation

$\displaystyle\langle(\hat{x}_{h},\hat{z}),\mathbf{T}_{G}(x_{0},w)\rangle=\langle\mathbf{T}_{G}^{*}(\hat{x}_{h},\hat{z}),(x_{0},w)\rangle.$

The next lemma shows how the state-space realization of the adjoint operator
is related to the periodic system (3).

######  Lemma 2.

If $(\hat{x}_{0},\hat{w})=\mathbf{T}_{G}^{*}(\hat{x}_{h},\hat{z})$ then

$\displaystyle\dot{\hat{x}}(t)$
$\displaystyle=-A^{*}(t)\hat{x}(t)-C^{*}(t)\hat{z}(t),~{}~{}\hat{x}(h)=\hat{x}_{h}$
$\displaystyle\hat{w}(t)$
$\displaystyle=B^{*}(t)\hat{x}(t)+D^{*}(t)\hat{z}(t),~{}~{}\hat{x}_{0}=\hat{x}(0)$
(19)

###### Proof.

The proof can be found in the Appendix. ∎

The following lemma provides a useful interpretation for the Hamiltonian
dynamics (17).

######  Lemma 3.

Interconnect the dynamics of (18) and (19) with $\hat{z}(t):=z(t)$ and
$w(t):=\gamma^{-2}\hat{w}(t)$. The resulting autonomous dynamics has state
$\tilde{e}^{T}:=[x^{T}\,\hat{x}^{T}]$ with dynamics given by (suppressing
dependence on $t$):

$\displaystyle\dot{\tilde{e}}$ $\displaystyle=\tilde{H}\tilde{e}$ (20)
$\displaystyle w$
$\displaystyle=(\gamma^{2}I-D^{*}D)^{-1}\left[\\!\begin{array}[]{cc}D^{*}C&B^{*}\end{array}\\!\right]\tilde{e}$
(22) $\displaystyle z$
$\displaystyle=(\gamma^{2}I-DD^{*})^{-1}\left[\\!\begin{array}[]{cc}\gamma^{2}C&DB^{*}\end{array}\\!\right]\tilde{e}$
(24)

Moreover let $(\lambda,v)$ denote an eigenvalue/eigenvector of the monodromy
matrix $\tilde{Q}$ for (20). Partition $v^{*}:=[v_{1}^{*}\,v_{2}^{*}]$
conformably with the state $\tilde{e}^{T}:=[x^{T}\,\hat{x}^{T}]$. Then
$\mathbf{T}_{G}(v_{1},w)=(\lambda v_{1},z)$ and $\mathbf{T}_{G}^{*}(\lambda
v_{2},z)=(v_{2},\gamma^{2}w)$.

###### Proof.

The expression for the Hamiltonian dynamics as the feedback connection of (18)
and (19) is from [10]. The expressions for $\mathbf{T}_{G}$ and
$\mathbf{T}_{G}^{*}$ follow from the definitions of the forward and adjoint
operators. ∎

In what follows we show that worst-case input can be constructed from (22) if
the Hamiltonian system (20) is initialized as $\tilde{e}(0)=v$. For this we
need the following theorem, which is the direct application of Theorem 21.12
in [14] to the discrete-time system $\underline{G}_{\gamma}$:

######  Theorem 2.

The following statements are equivalent:

  1. (a)

$\|\underline{G}_{\gamma}\|_{\infty}<1$.

  2. (b)

$\tilde{Q}$ has no eigenvalues on the unit circle and
$\|\underline{C}_{\gamma}(I-\underline{A}_{\gamma})^{-1}\underline{B}_{\gamma}\|<1$

It follows form Theorems 2.1 and 2.2 that if $\gamma\leq\|\mathbf{G}\|$ then
$\tilde{Q}$ has a unit modulus eigenvalue $e^{j\omega h}$ for some $\omega$.
Let $v$ denote the corresponding eigenvector. Let $\tilde{e}$ denote the
solution of the Hamiltonian dynamics (20) and $(w,z)$ the corresponding
outputs with initial condition $e(0)=v$. Then Lemma 3 implies that
$\mathbf{G}$ maps the input $w$ and initial condition $v_{1}$ to the output
$z$. Moreover, it follows from Lemma 3 that:

$\displaystyle v_{2}^{*}v_{1}+\int_{0}^{h}z^{*}(t)z(t)$
$\displaystyle=\langle(e^{j\omega h}v_{2},z),\mathbf{T}_{G}(v_{1},w)\rangle$
(25) $\displaystyle=\langle\mathbf{T}_{G}^{*}(e^{j\omega
h}v_{2},z),(v_{1},w)\rangle$ (26)
$\displaystyle=v_{2}^{*}v_{1}+\gamma^{2}\int_{0}^{h}w^{*}(t)w(t)$ (27)

Thus $||z||_{[0,h)}=\gamma||w||_{[0,h)}$, i.e. the $\mathcal{L}_{2}$ gain of
$\mathbf{G}$ is equal to $\gamma$ on the interval $[0,h)$ provided the initial
condition of the system $\mathbf{G}$ is given by $x(0)=v_{1}$.

Note that the input/output pair $(w,z)$ for $t\in[0,h)$ is obtained by
integrating the Hamiltonian dynamics (20) starting from the initial condition
$e(0)=v$. The state of the Hamiltonian system after one period is given by
$e(h)=\tilde{Q}e(0)=e^{j\omega h}v$. Thus integrating the periodic Hamiltonian
dynamics forward over the next interval yields $w(t)=e^{jwh}w(t-h)$ and
$z(t)=e^{jwh}z(t-h)$ for $t\in[h,2h)$. Continuing to evolve the Hamiltonian
dynamics forward in time over subsequent periodic intervals yields an
input/output pair of $\mathbf{G}$ such that for any interval
$k\in\\{0,1,2,\ldots\\}$:

$\displaystyle w(t)$ $\displaystyle=e^{j\omega
kh}w(t-kh),\,\,\,t\in[kh,(k+1)h)$ (28) $\displaystyle z(t)$
$\displaystyle=e^{j\omega kh}z(t-kh)$ (29)

This input/output pair satisfies
$\|z\|_{[kh,(k+1)h)}=\gamma\|w\|_{[kh,(k+1)h)}$ over each interval
$[kh,(k+1)h)$. Thus $w$ is, loosely speaking, an input that achieves the gain
$\gamma$. There are three issues to be resolved to make this more precise.
First, the input (28) is persistent and hence is not in
$\mathcal{L}_{2,[0,\infty)}$. Second, the signals $(w,z)$ in (28-29) are an
input/output pair of $\mathbf{G}$ only if the periodic system starts with the
initial condition $x(0)=v_{1}$. However, the induced gain is defined assuming
$x(0)=0$. Third, the input can be complex if the eigenvector $v$ is complex.

The first two issues are resolved by noting that $(w,z)$ can be expressed
using the state transition matrix:

$\displaystyle z(t)=C(t)\Phi(t,0)v_{1}$
$\displaystyle+\int_{0}^{t}C(t)\Phi(t,\tau)B(\tau)w(\tau)\,d\tau+D(t)w(t)$

Define $w^{\bullet}$ as the truncation of $w$ after $K$ periodic intervals,
i.e. $w^{\bullet}(t):=w(t)$ for $t\in[0,Kh)$ and $w^{\bullet}(t):=0$
otherwise. Define $z^{\bullet}$ as the output of $\mathbf{G}$ driven by input
$w^{\bullet}$ but with initial condition $x(0)=0$. By causality,
$z^{\bullet}(t)=z(t)-s(t)$ for $t<Kh$ where $s(t):=C(t)\Phi(t,0)v_{1}$ is the
initial condition response. The norm of $z^{\bullet}$ can be bounded as:

$\displaystyle\|z^{\bullet}\|$
$\displaystyle\geq\|z\|_{[0,Kh)}-\|s\|_{[0,Kh)}$
$\displaystyle=\gamma\|w^{\bullet}\|_{[0,Kh)}-\|s\|_{[0,Kh)}$

The first line follows from the triangle inequality. The second line follows
from two facts. First, $\|z\|_{[kh,(k+1)h)}=\gamma\|w\|_{[kh,(k+1)h)}$ over
each interval as noted above. Second, $w^{\bullet}=w$ by construction for
$t\in[0,Kh)$. Next note that $\|s\|_{[0,Kh)}\leq\|s\|<\infty$ due to the
stability of $\mathbf{G}$. In addition,
$\|w^{\bullet}\|_{[0,Kh)}=K\|w^{\bullet}\|_{[0,h)}$ and hence
$\|w^{\bullet}\|_{[0,Kh)}\rightarrow\infty$ as $K\rightarrow\infty$. Thus it
is clear that for all $\epsilon>0$ there exists an integer $K$ such that the
input-output pair $(w^{\bullet},z^{\bullet})$ of $\mathbf{G}$ satisfies
$\frac{\|z^{\bullet}\|}{\|w^{\bullet}\|}\geq\gamma-\epsilon$. The input
$w^{\bullet}$ is a valid “worst-case” signal because it is in
$\mathcal{L}_{2}$ and the output has been generated with zero initial
conditions.

The remaining issue is the fact that $(w^{\bullet},z^{\bullet})$ may be
complex. The periodic system $\mathbf{G}$ is linear and the system matrices
are real. Thus neglecting the imaginary part of the complex valued input does
not change the gain, i.e. the ratio of the norms of the real input and the
corresponding real output remains the same. Therefore the real valued worst-
case input can be obtained by defining $w^{\circ}=Re(w^{\bullet})$. Now we can
summarize the complete algorithm.

######  Algorithm 1 (Worst-case input).

1: Let $\underline{\gamma}$ be given such that
$\underline{\gamma}\leq\|\mathbf{G}\|$. Set $\gamma:=\underline{\gamma}$.

2: Compute $Q$ from $A_{\gamma},B_{\gamma},C_{\gamma}$ using the equation
(13). Determine $\tilde{Q}$ by applying the similarity transformation $T$.

3: Compute the unit modulus eigenvalue $e^{j\omega h}$ of $\tilde{Q}$ and the
corresponding eigenvector
$v=\left[\\!\begin{array}[]{cc}\overline{x}^{*}&\overline{\hat{x}}^{*}\end{array}\\!\right]^{*}$.

4: Integrate the Hamiltonian dynamics (17) from $t=0$ to $t=h$ starting from
the initial condition $\tilde{e}(0)=v$.

5: Use (22) to compute $w(t)$ from the state trajectories $x(t),\hat{x}(t)$
over $[0,h)$.

6: Choose an integer $K\gg 1$ and define the (possibly complex valued) signal
$w^{\bullet}(t)$ by:  $\displaystyle w^{\bullet}(t)=\begin{cases}e^{j\omega
kh}w(t-kh)&\text{ if }t\in[kh,(k+1)h),\\\ &~{}~{}~{}k=0,1,\ldots,K-1\\\
0&\text{ if }t\geq Kh\end{cases}$

7: Let $w^{\circ}:=\text{Re}(w^{\bullet})$ and
$z^{\circ}:=\mathbf{G}w^{\circ}$.

##  III Computation of lower bounds for the induced $\mathcal{L}_{2}$ norm of
LPV systems

###  III-A Problem formulation

Linear parameter varying (LPV) systems are a class of systems whose state
space matrices depend on a time-varying parameter vector $\rho$. We assume
that $\rho:\mathbb{R}^{+}\rightarrow\mathbb{R}^{m}$ is a piecewise
continuously differentiable function of time and satisfies the known bounds:

$\displaystyle\underline{\rho_{i}}\leq\rho_{i}(t)\leq\overline{\rho}_{i},~{}~{}\underline{\mu_{i}}\leq\dot{\rho}_{i}(t)\leq\overline{\mu}_{i},~{}1\leq
i\leq m$ (30)

The set of allowable parameter vectors, denoted
$\mathcal{P}\subseteq\mathbb{R}^{m}$ consists of vectors $\rho$ that satisfy
the range bounds given in (30). The set of admissible trajectories, denoted
$\mathcal{A}$, consists of piecewise continuously differentiable trajectories
that satisfy both the rate and range bounds given also in (30).

The state-space matrices of an LPV system are continuous functions of the
parameter: $A:\mathcal{P}\rightarrow\mathbb{R}^{n\times n}$,
$B:\mathcal{P}\rightarrow\mathbb{R}^{n\times p}$,
$C:\mathcal{P}\rightarrow\mathbb{R}^{q\times n}$ and
$D:\mathcal{P}\rightarrow\mathbb{R}^{q\times p}$. An $n^{\textrm{th}}$ order
LPV system, $G$, is defined by the following state-space form:

$\begin{split}\dot{x}(t)&=A(\rho(t))x(t)+B(\rho(t))w(t)\\\
z(t)&=C(\rho(t))x(t)+D(\rho(t))w(t)\end{split}$ (31)

The performance of an LPV system $G$ can be specified in terms of its induced
$\mathcal{L}_{2}$ gain from input $w$ to output $z$. The induced
$\mathcal{L}_{2}$-norm is defined by

$\|G\|=\sup_{0\neq
w\in\mathcal{L}_{2}(\mathbb{R}^{p}),\rho(\cdot)\in\mathcal{A}}\frac{\|z\|}{\|w\|},$
(32)

The initial condition is assumed to be $x(0)=0$. The notation
$\rho(\cdot)\in\mathcal{A}$ refers to the entire (admissible) trajectory as a
function of time. The class of LPV system given above has an arbitrary
dependence on the parameter. For this class of systems there are known linear
matrix inequality (LMI) conditions to efficiently compute an upper bound on
the induced $\mathcal{L}_{2}$ gain [5], [6].

Now we address the complementary problem, i.e. our aim is to determine a lower
bound for the induced $\mathcal{L}_{2}$ norm by using the results of the
previous section.

###  III-B Lower bound for the induced $\mathcal{L}_{2}$ norm

The computation of the lower bound is based on restricting the scheduling
parameter trajectories to a finite-dimensional set of periodic signals. Let
$\rho:\mathbb{R}^{+}\times\mathbb{R}^{N}\rightarrow\mathbb{R}^{m}$ denote a
function that specifies a periodic scheduling trajectory for each value of
$c\in\mathbb{R}^{N}$. $\rho(\cdot,c)$ denotes the entire trajectory (as a
function of time) at the particular value $c$ and $\rho(t,c)$ denotes the
$m$-dimensional scheduling vector obtained by evaluating $\rho(\cdot,c)$ at
time instant $t$. The trajectory is assumed to be periodic, i.e. for each $c$
there is a period $h(c)$ such that $\rho(t+h(c),c)=\rho(t,c)$ $\forall t$. In
addition, we must ensure the trajectory is admissible in the sense that it
satisfies the range and rate bounds in (30). Let
$\mathcal{C}_{p}\subseteq\mathbb{R}^{N}$ denote the set of values that lead to
such admissible trajectories, i.e. $\rho(\cdot,c)\in\mathcal{A}$ for all
$c\in\mathcal{C}_{p}$. The corresponding set of periodic trajectories is
defined as

$\displaystyle\mathcal{A}_{p}:=\left\\{\rho(\cdot,c)\,|\,c\in\mathcal{C}_{p}\right\\}$
(33)

As a concrete example, $\rho$ can be specified as a linear combination of
periodic bases functions $\\{\phi_{k}\\}_{k=1}^{N}$, i.e.
$\rho(\cdot,c):=\sum_{k=1}^{N}c_{k}\phi_{k}(\cdot)$. Section III-C provides
alternative characterizations to specify periodic trajectories.

Define the lower bound $\gamma_{lb}$ on the induced $\mathcal{L}_{2}$ gain
(32) as

$\displaystyle\gamma_{lb}:=\sup_{\rho(\cdot,c)\in\mathcal{A}_{p}}\|\mathbf{G}_{\rho(\cdot,c)}\|=\sup_{c\in\mathcal{C}_{p}}\|\mathbf{G}_{\rho(\cdot,c)}\|$
(34)

where $\mathbf{G}_{\rho(\cdot,c)}$ denotes the periodic system (operator)
obtained by evaluating the LPV system $G$ along the periodic trajectory
specified by $\rho(\cdot,c)$. The algorithm described in [9] can thus be used
to evaluate the gain $\|\mathbf{G}_{\rho(\cdot,c)}\|$. It follows immediately
from $\mathcal{A}_{p}\subset\mathcal{A}$ that $\gamma_{lb}\leq\|G\|$. Equation
(34) defines a finite dimensional optimization problem, which is non-convex in
general. One further issue is that a single, accurate evaluation of
$\|\mathbf{G}_{\rho(\cdot,c)}\|$ requires many bisection steps and, as a
consequence the matrix differential equations (Equations (15), (16) and (14))
must be integrated many times for a single evaluation of the objective
function. Thus the evaluation of $\|\mathbf{G}_{\rho(\cdot,c)}\|$ is
computationally costly. A significant reduction in computation time can be
achieved by using Algorithm 2 described below. In this algorithm,
$\underline{G}_{\gamma,\rho(\cdot,c)}$ denotes the discrete-time system (7)
corresponding to the PLTV operator $\mathbf{G}_{\rho(\cdot,c)}$. Moreover,
define $\nu(c,\gamma):=\|\underline{G}_{\gamma,\rho(\cdot,c}\|$. With this
notation, $\|G_{\rho(\cdot,c)}\|<\gamma$ if and only if $\nu(c,\gamma)<1$. The
lower bound algorithm can now be stated.

######  Algorithm 2 (Lower bound on $\mathcal{L}_{2}$ norm-I).

1: Pick an initial parameter vector $c_{0}$. Let $c:=c_{0}$.

2: Compute the norm of $\mathbf{G}_{\rho(\cdot,c_{0})}$ by using the algorithm
in [9]. Take the upper bound $\overline{\gamma}$ from the bisection. It is
clear that $\nu(c_{0},\bar{\gamma})<1$ because $\bar{\gamma}$ is an upper
bound for $\|\mathbf{G}_{\rho(\cdot,c_{0})}\|$.

3: Starting from the initial value $c_{0}$, solve the nonlinear optimization
problem  $\displaystyle\sup_{c\in\mathcal{C}_{p}}\nu(c,\bar{\gamma})$ (35)
Three different outcomes are possible: (35) terminates at a (local) optimum
$c^{*}$ where $\nu(c^{*},\gamma)<1$ (case-a), $\nu(c^{*},\gamma)\geq 1$
(case-b), or the optimization stops because at some $c=c^{*}$ the solution
$Z(t)$ of (14) goes unbounded (case-c). Clearly, case-b and case-c mean that
$\bar{\gamma}$ is smaller than $\|\mathbf{G}_{\rho(\cdot,c^{*})}\|$ so the
optimization managed to improve the lower bound. Therefore, set $c_{0}:=c^{*}$
and go to step 2. In case-a the optimization above was unsuccessful in the
sense that the norm was not significantly improved. In this case, go to step
4.

4: Compute the norm of $\mathbf{G}_{\rho(\cdot,c^{*})}$ and take the lower
bound $\underline{\gamma}$ from the bisection. Let
$\gamma_{lb}:=\underline{\gamma}$ and stop.

The advantages of Algorithm 2 over the direct maximization of
$\|\mathbf{G}_{\rho(\cdot,c)}\|$ are the following: the bisection has to be
performed less times (only once before each optimization step (35) and then
once at the end of the procedure) and the computation of $\nu$ requires only a
single integration of the matrix differential equations, i.e. it can be
evaluated with significantly less computation than
$\|\mathbf{G}_{\rho(\cdot,c)}\|$. Testing the algorithm on numerical examples,
we found in most cases that only a few or even just one optimization step is
enough to get a good approximation for the lower bound. The problem is often
with the solution of the Riccati equation (14), which goes unbounded when the
norm of $\mathbf{G}_{\rho(\cdot,c)}$ becomes significantly larger than
$\bar{\gamma}$ at some $c=\hat{c}$. If this happens, a new $\bar{\gamma}$ has
to be computed by performing the time-consuming bisection algorithm. This
situation can be avoided if (35) is started with some
$\bar{\gamma}\gg\|\mathbf{G}_{\rho(\cdot,c_{0})}\|$. Clearly if $\bar{\gamma}$
is chosen to be an upper bound of the $\mathcal{L}_{2}$ norm of the LPV
system, then $\|\mathbf{G}_{\rho(\cdot,c)}\|\leq\|G\|\leq\bar{\gamma}$ for all
$c$, which implies that $Z(t)$ never goes unbounded. The modified, 1-step
algorithm is summarized as follows:

######  Algorithm 3 (Lower bound on $\mathcal{L}_{2}$ norm-II).

1: Compute an upper bound $\gamma_{ub}$ on the gain of the LPV system, i.e.
$\gamma_{ub}\geq\|G\|$. Such an upper bound can be determined by using
standard methods based on dissipativity relation (e.g. Bounded Real type LMI
conditions) [5], [6]. Let $\bar{\gamma}:=\gamma_{ub}$. Pick an initial
parameter vector $c_{0}$.

2: Solve the optimization problem (35).

3: Compute the norm of $\mathbf{G}_{\rho(\cdot,c^{*})}$ and take
$\underline{\gamma}$ from the outputs of the bisection. Let
$\gamma_{lb}:=\underline{\gamma}$ and stop.

This algorithm only requires the $\mathcal{L}_{2}$ nor of the periodic system
to be evaluated at the last step. This greatly reduces the number of required
integrations for the matrix differential equations associated with X, Y, and
Z. By performing either of the algorithms above we obtain a lower bound
$\gamma_{lb}$ for the induced $\mathcal{L}_{2}$ gain and the worst-case
scheduling trajectory $\rho(\cdot,c^{*})$, where the associated PLTV system
takes this norm. By using Algorithm 1 we can also construct a worst-case input
signal $w^{o}$ for $\mathbf{G}_{\rho(\cdot,c^{*})}$. This, together with
$\rho(\cdot,c^{*})$ gives the (worst-case input, worst-case scheduling
trajectory) pair, where the $\mathcal{L}_{2}$-gain of the LPV system is
$\gamma_{lb}$.

###  III-C Scheduling trajectories

Algorithms 2 and 3 optimize the $\mathcal{L}_{2}$ bound over the elements of
$\mathcal{A}_{p}$. Therefore it is important how this set is characterized,
i.e. how the periodic signals depend on the parameter vector $c$. The
structure of $\rho(\cdot,c)$ influences the convergence properties of (35) and
determines the final result. Therefore it has to be carefully chosen. There
are many ways to construct periodic signals. In this section we focus on
piecewise linear scheduling trajectories, while another construction method,
based on sinusoidal basis functions, is described in [17]. Both methods result
in a polytopic parameter set, i.e. $\mathcal{C}_{p}$ is defined by linear
inequality constraints. The linearity of the constraints simplifies the
optimization (35).

Let the scheduling trajectories be specified by restricting the time
derivative $\dot{\rho}$ to be piecewise constant. Specifically, consider the
one-dimensional case ($n_{\rho}=1$) where the time interval $[0,h]$ is
subdivided into $R$ intervals $[T_{j-1},T_{j})$ for $j=1,\ldots,R$ with
$T_{0}:=0$ and $T_{R}=h$. Let $c_{j}:=T_{j}-T_{j-1}$ denote the length of the
$j^{th}$ sub-interval. The parameter trajectory is then given as the periodic
function of period $h$ defined for $t\in[0,h]$ by:

$\displaystyle\rho(0,c)$ $\displaystyle:=\rho_{0}$ (36)
$\displaystyle\dot{\rho}(t,c)$ $\displaystyle:=r_{j}\mbox{ for
}t\in[T_{j-1},T_{j})$ (37)

where $r_{j}\in[\underline{\mu},\bar{\mu}]$ are user-selected values that
specify the rate on each interval. The vector
$c:=[\rho_{0},c_{1},\ldots,c_{R}]\in\mathbb{R}^{R+1}$ specifies the parameter
trajectory offset $\rho_{0}$ and the subinterval lengths. Equation (37) can be
integrated to give the explicit form of the trajectory as:

$\displaystyle\rho(t,c)$
$\displaystyle=\rho_{0}+\sum_{i=1}^{j-1}c_{i}r_{i}+(t-T_{j-1})r_{j}~{}~{}\mbox{
for }t\in[T_{j-1},T_{j})$ (38)

Note that the period $h=c_{1}+\cdots+c_{R}$ is a free parameter as it depends
on the interval lengths. Moreover, $\rho(t,c)$ is a linear function of the
offset $\rho_{0}$ and interval lengths $c_{j}$. Thus the magnitude constraints
in (30) can be transformed into $R+1$ linear constraints on $c$:

$\displaystyle\underline{\rho}\leq\rho_{0}+\sum_{i=1}^{j}c_{i}r_{i}\leq\bar{\rho},~{}~{}~{}~{}~{}j\in\\{0,\ldots,R\\}$
(39)

To ensure $\rho(0,c)=\rho(h,c)$, (39) has to be completed with an equality
condition $\sum_{i=1}^{R}c_{i}r_{i}=0$. Since all constraints are linear, thus
$\mathcal{C}_{p}$ is a polytope.

######  Remark 2.

The construction of the worst-case input requires additional considerations.
In particular, Algorithm 1 for constructing the worst-case input requires
integrating the unstable Hamiltonian system over $[0,h)$. This causes
numerical issues if the period length $h$ is ”too large”. Thus it is
reasonable to introduce an inequality constraint of the form $h\leq\bar{h}$ to
bound the period length when performing the optimization (24). This is simply
another linear constraint on the vector $c$.

##  IV Numerical examples

In this section three numerical examples are presented to demonstrate the
applicability of the proposed methods. To initialize our algorithms we need to
determine an upper bound $\gamma_{ub}$ for the induced $\mathcal{L}_{2}$ gain.
In all examples $\gamma_{ub}$ is computed by solving the following
optimization problem:

$\displaystyle\min_{V(x,\rho)}\gamma$ $\displaystyle
V(x,\rho)>0,\dot{V}(x,\rho,\dot{\rho})\leq\gamma^{2}w^{T}w-z^{T}z$ (40)

where the Lyapunov (storage) function $V(x,\rho)$ was chosen to be quadratic:
$V(x,\rho)=x^{T}P(\rho)x$, $P(\rho)=P(\rho,P_{0},P_{1},\ldots,P_{M})$, where
$P_{i}$-s denote the free (matrix) variables to be found. The infinite LMI
constraints obtained were transformed to a finite set by choosing a suitable
dense grid over the parameter domain $\mathcal{P}$ and only the inequalities
evaluated at the grid points are considered [4].

Our strategy to compute the lower bound on the $\mathcal{L}_{2}$ gain was the
following: first we performed Algorithm 3 and then we used the parameter
trajectory and lower bound we obtained as an initial value to run Algorithm 2.
Since in all numerical examples this second optimization cannot improve the
previous results, thus the numerical results we described in the forthcoming
sections are generated by Algorithm 3. The nonlinear optimization was
performed under MATLAB, by using the pattern search algorithm implemented in
the Global Optimization Toolbox.

###  IV-A LPV system with gain-scheduled PI controller

The first example, taken from [18] is a feedback interconnection of a first-
order LPV system with a gain-scheduled proportional-integral controller. The
state-space matrices of the closed-loop system are as follows

$\displaystyle
A(\rho):=\left[\\!\begin{array}[]{cc}-\frac{1}{\tau(\rho)}(1+K_{p}(\rho)K(\rho))&\frac{1}{\tau(\rho)}\\\
-K_{i}(\rho)K(\rho)&0\end{array}\\!\right],$ $\displaystyle
B(\rho):=\left[\\!\begin{array}[]{c}\frac{1}{\tau(\rho)}K_{p}(\rho)\\\
K_{i}(\rho)\end{array}\\!\right],$ (45) $\displaystyle
C(\rho):=[-K(\rho)~{}0],$ $\displaystyle D:=1$ (46)

where $\tau(\rho):=\sqrt{13.6-16.8\rho}$, $K(\rho):=\sqrt{4.8\rho-8.6}$ and

$\displaystyle
K_{p}(\rho)=\frac{2\xi_{cl}\omega_{cl}\tau(\rho)-1}{K(\rho)},~{}~{}K_{i}(\rho)=\frac{\omega_{cl}^{2}\tau(\rho)}{K(\rho)},~{}~{}\xi_{cl}=0.7,~{}~{}\omega_{cl}=0.25.$

The scheduling parameter is assumed to vary in the interval $[2,7]$ and
$\dot{\rho}\in[-1,~{}1]$. By performing the optimization (40) with

$\displaystyle
V(x,\rho)=x^{T}\left(P_{0}+\sum_{k=1}^{6}\rho^{k}P_{k}+\frac{1}{\rho}P_{7}+\frac{1}{\rho^{2}}P_{8}+\frac{1}{\rho^{3}}P_{9}\right)x$
$\displaystyle\Gamma:=\\{\rho_{1}=2,\ldots,\rho_{100}=7\\},~{}~{}\rho_{k+1}-\rho_{k}=5/99$

we got $\gamma_{ub}=2.964$ for the upper bound. Using the parameter values in
$\Gamma$ the frozen lower bound was also computed as
$\gamma_{lb,fr}=\max_{k}\|G_{\rho_{k}}\|$, $\rho_{k}\in\Gamma$, where
$G_{\rho_{k}}$ denotes the LTI system obtained by substituting
$\rho(t)=\rho_{k}$ for all $t$. The lower bound we obtained is
$\gamma_{lb,fr}=1.1066$. To compute the lower bound by Algorithm 3, we chose
the following pattern for the rate variation of the scheduling trajectory
$r=[\overline{\mu}~{}0~{}\underline{\mu}~{}0~{}\overline{\mu}~{}0]=[1~{}0~{}-1~{}0~{}1~{}0]$
and let the algorithm tune the period in the interval
$[h_{min}~{}h_{max}]=[12~{}20]$. The lower bound we obtained is
$\gamma_{lb,pwl}=2.84$, at $h=15.877$. The scheduling trajectory can be seen
in Fig. 1. This lower bound is is significantly larger than
$\gamma_{lb,\text{\cite[cite]{[\@@bibref{}{peni14_accinprep}{}{}]}}}$, where
$\gamma_{lb,\text{\cite[cite]{[\@@bibref{}{peni14_accinprep}{}{}]}}}$ is the
lower bound we obtained in [17] by using sinusoid scheduling trajectories.
Moreover, $\gamma_{lb,pwl}$ is very close to $\gamma_{ub}$: the difference is
only $\gamma_{ub}-\gamma_{lb,pwl}=0.13$. This means that we have very tight
bounds on the norm of $G$:
$2.84=\gamma_{lb,pwl}\leq\|G\|\leq\gamma_{ub}=2.964$.

Finally, by using Algorithm 1 we computed also the worst case input for both
scheduling trajectories. The input signals were constructed starting from the
unit modulus eigenvalue $-1$. The parameter $K$ in Algorithm 1 was chosen to
be $K=60$. The result can be seen in Fig. 1.

Figure 1: Worst case scheduling trajectories and worst-case input/output for
the LPV system in Section IV-A. The first two figures show the trajectory of
the worst-case scheduling parameter $\rho(t)$ and its time derivative
$\dot{\rho}(t)$. In the 3rd and 4th figures $w^{o}$ is depicted over a longer
and a shorter ($[0~{}8h]$) time intervals. In the 5th and 6th figures $z^{o}$
is displayed over the same time intervals.

###  IV-B Input and output scaled LTI system

The next example was constructed by taking two copies of the simple LTI system
$1/(s+1)$, scaling the input of the first and he output of the second and
computing the difference of the two outputs. The dynamics are given as
follows:

$\displaystyle\left[\\!\begin{array}[]{c}\dot{x}_{1}\\\
\dot{x}_{2}\end{array}\\!\right]$
$\displaystyle=\left[\\!\begin{array}[]{cc}-1&0\\\
0&-1\end{array}\\!\right]\left[\\!\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\\!\right]+\left[\\!\begin{array}[]{c}1\\\
\rho\end{array}\\!\right]w$ $\displaystyle z$
$\displaystyle=\left[\\!\begin{array}[]{cc}\rho&-1\end{array}\\!\right]\left[\\!\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\\!\right]$

We assume that $\rho\in[-1~{}~{}1]$ and
$\dot{\rho}\in[-\overline{\mu}~{}~{}\overline{\mu}]$. We are going to compute
the lower bound for different values of $\overline{\mu}$. It follows from the
structure of the system that if $\rho$ is constant then the difference between
the input and the output scaled systems is 0. This implies that
$\gamma_{lb,fr}=0$. Next, the upper bound on the $\mathcal{L}_{2}$ gain was
computed for different rate bounds. In all cases the storage function
$V(x,\rho)$ and the parameter grid $\Gamma$ was chosen as follows:

$\displaystyle V(x,\rho)$
$\displaystyle=x^{T}\left(P_{0}+\sum_{k=1}^{10}\rho^{k}P_{k}\right)x$
$\displaystyle\Gamma$
$\displaystyle=\\{\rho_{1}=-1,\ldots,\rho_{100}=1\\},~{}~{}\rho_{k+1}-\rho_{k}=2/99$

The upper bounds obtained are collected in Table I. To compute the lower bound
we used the following rate pattern to characterize the piecewise linear
scheduling trajectories:
$r=[\overline{\mu}~{}0~{}-\overline{\mu}~{}0~{}\overline{\mu}~{}0]$. The lower
bounds and the period of the worst-case scheduling trajectories obtained for
the different $\overline{\mu}$ values are collected in Table I. The range
bounds $h_{\min}$ and $h_{\max}$ used to constrain the period length in (35)
are also given in Table I. It can be seen that the upper and lower bounds are
very close to each other, which means that by using the upper and lower bound
algorithms we could precisely compute the norm of this LPV system. In two
particular cases, $\overline{\mu}=1$ and $\overline{\mu}=2$, we plotted also
the worst case scheduling trajectory and the worst-case input in Fig. 3. (The
worst case inputs were computed starting from the unit modulus eigenvalues
$[0.9723+0.2337i]$(in case $\overline{\mu}=1$) and $[0.9704+0.2415i]$ (in case
$\overline{\mu}=2$). $K$ was 60 in both cases.)

$\displaystyle\\!\begin{array}[]{c|c|c|c|c|c|c|c|c|c|c|c|c}\overline{\mu}&\gamma_{lb}&\gamma_{ub}&\gamma_{ub}-\gamma_{lb}&[h_{min}~{}h_{max}]&h\\\
\hline\cr 0.1&0.0979&0.1087&0.0108&[0.5~{}50]&49.9\\\
0.4&0.3309&0.3342&0.0033&[0.5~{}50]&44.2\\\
0.7&0.4783&0.4805&0.0022&[0.5~{}20]&19.78\\\
1&0.5645&0.5766&0.0121&[0.5~{}6]&6\\\ 1.3&0.6364&0.6435&0.0071&[0.5~{}6]&6\\\
1.6&0.6874&0.6924&0.0050&[0.5~{}6]&6\\\
2&0.7347&0.7403&0.0056&[0.5~{}6]&6\end{array}\\!$ TABLE I: Upper and lower
bounds on the $\mathcal{L}_{2}$ norm of the LPV system defined in Section IV-B
Figure 2: Worst case scheduling trajectories and worst-case inputs/outputs at
rate bounds $\overline{\mu}=1$ for the example in Section IV-B. The first two
subfigures show the trajectory of the worst-case scheduling parameter
$\rho(t)$ and its time derivative $\dot{\rho}(t)$. In the 3rd and 4th
subfigures $w^{o}$ is depicted over a longer and a shorter ($[0~{}8h]$) time
interval. In the 5th and 6th figures $z^{o}$ is displayed over the same time
intervals.  Figure 3: Worst case scheduling trajectories and worst-case
inputs/outputs at rate bounds $\overline{\mu}=2$ for the example in Section
IV-B. The first two subfigures show the trajectory of the worst-case
scheduling parameter $\rho(t)$ and its time derivative $\dot{\rho}(t)$. In the
3rd and 4th subfigures $w^{o}$ is depicted over a longer and a shorter
($[0~{}8h]$) time interval. In the 5th and 6th figures $z^{o}$ is displayed
over the same time intervals.  Figure 4: Upper and lower bounds on the induced
$\mathcal{L}_{2}$ norm computed for the example in Section IV-B.

###  IV-C Rotated LTI system

The next example was constructed by taking a stable LTI system and performing
a parameter dependent similarity transformation, which is actually a rotation,
on the $A$ matrix. The system can be given in state-space form as follows:

$\displaystyle\dot{x}$
$\displaystyle=R(\rho)^{T}\left[\\!\begin{array}[]{cc}-0.5&-0.4\\\
3&-0.5\end{array}\\!\right]R(\rho)x+\left[\\!\begin{array}[]{c}0.5\\\
0.5\end{array}\\!\right]w$ $\displaystyle z$ $\displaystyle=[1~{}-1]x+0.1w$

where $R(\rho)=\left[\footnotesize\begin{array}[]{cc}\cos(\rho)&\sin(\rho)\\\
-\sin(\rho)&\cos(\rho)\end{array}\right]$. The range and rate bounds for
$\rho$ were chosen to be: $\rho\in[\pi/4~{}~{}\pi/2]$ and
$\dot{\rho}\in[-0.1~{}~{}0.1]$. The upper bound on the induced
$\mathcal{L}_{2}$ norm was computed by using the following storage function
and parameter grid:

$\displaystyle
V(x,\rho)=x^{T}\left(P_{0}+\sum_{k=1}^{7}\rho_{1}^{k}P_{k}\right)x$
$\displaystyle\Gamma=\\{\rho_{1}=\pi/4,\ldots,\rho_{50}=\pi/2\\},$
$\displaystyle\rho_{k+1}-\rho_{k}=\pi/196,~{}~{}~{}k\in\\{1\ldots 50\\}$

We obtained $\gamma_{ub}=3.3$. The lower bound computed at frozen parameter
values over the grid above was $\gamma_{lb,fr}=2.696$. To compute the lower
bound we chose the following derivative pattern for the scheduling trajectory:
$r=0.1\cdot[1~{}0~{}-1~{}0~{}1~{}0]$. The algorithm was allowed to tune the
period length in the interval $[1~{}5]$. The lower bound we got is
$\gamma_{lb}=3.15$, ($h=3.11$), which is again significantly larger than
$\gamma_{lb,fr}$. The worst-case scheduling trajectories and the worst-case
input are plotted in Fig. 5. (The worst-case input was constructed by starting
from the unit modulus eigenvalue $-0.99997+0.00816i$. The parameter $K$ was
set to $60$.)

Figure 5: Worst case scheduling trajectories and worst-case input/output for
the LPV system in Section IV-C. The first two figures show the trajectory of
the worst-case scheduling parameter $\rho(t)$ and its time derivative
$\dot{\rho}(t)$. In the 3rd and 4th figures $w^{o}$ is depicted over a longer
and a shorter ($[0~{}8h]$) time intervals. In the 5th and 6th figures $z^{o}$
is displayed over the same time intervals.

###  IV-D A 2-parameter system

The next, more complex example was constructed from the examples in sections
IV-A and IV-B. We took two copies of the closed loop system (46) with
scheduling parameter $\rho(t):=\rho_{1}(t)$, then we scaled the input of the
first and the output of the second system by the same time-varying parameter
$\rho_{2}(t)$ and defined the output as the difference between the outputs the
two subsystems (like in Section IV-B). The state-space equations of the system
obtained can be written as

$\displaystyle\left[\\!\begin{array}[]{c}\dot{x}_{1}\\\
\dot{x}_{2}\end{array}\\!\right]$
$\displaystyle=\left[\\!\begin{array}[]{cc}A(\rho_{1})&0\\\
0&A(\rho_{1})\end{array}\\!\right]\left[\\!\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\\!\right]+\left[\\!\begin{array}[]{c}B(\rho_{1})\\\
\rho_{2}B(\rho_{1})\end{array}\\!\right]w$ $\displaystyle z$
$\displaystyle=\left[\\!\begin{array}[]{cc}\rho_{2}C(\rho_{1})&-C(\rho_{1})\end{array}\\!\right]\left[\\!\begin{array}[]{c}x_{1}\\\
x_{2}\end{array}\\!\right]$

The range and rate bounds were as follows: $2\leq\rho_{1}\leq 7$,
$-1\leq\rho_{2}\leq 1$ and $-1\leq\dot{\rho}_{1},\dot{\rho}_{2}\leq 1$. The
construction of the system implies that $\gamma_{lb,fr}=0$. The upper bound on
the induced $\mathcal{L}_{2}$ norm was computed by using the following storage
function and parameter grid:

$\displaystyle
V(x,\rho)=x^{T}\left(P_{0}+\sum_{k=1}^{3}\rho_{1}^{k}P_{k}+\sum_{\ell=1}^{5}\rho_{2}^{\ell}P_{\ell}+\frac{1}{\rho_{1}}P_{9}+\frac{1}{\rho_{2}}P_{10}\right)x$
$\displaystyle\Gamma=\\{\rho_{1,1}=2,\ldots,\rho_{1,30}=7\\}\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\{\rho_{2,1}=-1,\ldots,\rho_{2,10}=1\\},$
$\displaystyle\rho_{k+1,1}-\rho_{k,1}=2/29,~{}~{}\rho_{\ell+1,2}-\rho_{\ell,2}=2/9,$
$\displaystyle k\in\\{1\ldots 30\\},~{}~{}\ell=\\{1\ldots 10\\}$

We obtained $\gamma_{ub}=5.38$. To compute the lower bound we chose piecewise
linear scheduling trajectories with the following derivative patterns:
$r_{1}=[1~{}0~{}-1~{}0~{}1~{}0]$ and $r_{2}=[0~{}1~{}0~{}-1~{}0~{}1~{}0]$. The
lower bound we got is $\gamma_{lb}=5.047$, ($h=16.2142$), which is again very
close to the upper bound. The worst-case scheduling trajectories and the
worst-case input are plotted in Fig. 6. (The worst-case input was constructed
by starting from the unit modulus eigenvalue $-0.9997+0.0227i$. The parameter
$K$ was set to $60$.)

Figure 6: Worst case scheduling trajectories and worst-case input/output for
the LPV system in Section IV-D. The first two figures show the trajectory of
the worst-case scheduling parameter $\rho(t)$ and its time derivative
$\dot{\rho}(t)$. In the 3rd and 4th figures $w^{o}$ is depicted over a longer
and a shorter ($[0~{}8h]$) time intervals. In the 5th and 6th figures $z^{o}$
is displayed over the same time intervals.

##  V Conclusion

In the paper a numerical method is proposed for computing the lower bound of
the induced $\mathcal{L}_{2}$-gain of continuous-time, LPV systems. The
algorithm finds this bound – together with the worst-case parameter trajectory
– by using nonlinear optimization over periodic scheduling parameter
trajectories. Restricting the domain of parameter trajectories to periodic
signals enables to use the recent results for exact calculation of the
$\mathcal{L}_{2}$ norm for a periodic time varying system. It was shown that
the proposed algorithm can be reliably implemented by standard numerical tools
and provides precise approximation for the $\mathcal{L}_{2}$ bound.

##  Appendix A Proof of Lemma 2

To prove the theorem we need the following properties of the state transition
matrix: for all $t,\tau,\tau^{\prime}\in\mathbb{R}$

  1. (a)

$\Phi(t,\tau)\Phi(\tau,t)=I$, which implies $\Phi(t,\tau)=\Phi(\tau,t)^{-1}$

  2. (b)

$\frac{d}{dt}\Phi(\tau,t)=-\Phi(\tau,t)A(t)$

  3. (c)

$\frac{d}{dt}[\Phi(\tau,t)^{*}]=-A(t)^{*}\Phi(\tau,t)^{*}$.

Item (a) follows from the definition and from the regularity [11] of
$\Phi(t,\tau)$; item (b) is a consequence of item (a) and can be proved by
applying the chain rule of derivation to
$0=\frac{d}{dt}\Phi(t,\tau)\Phi(\tau,t)$ and finally, item (c) is the direct
consequence of item (b).

Then Lemma 2 can be proved directly from the definition of the adjoint
(Equation (19)) as follows:

$\displaystyle\langle(\hat{x}_{h},\hat{z}),\mathbf{T}_{G}(x_{0},w)\rangle=\hat{x}_{h}^{*}x(h)+\int_{0}^{h}\hat{z}(t)^{*}z(t)\,dt$
(47)

where $(x(h),z(t))$ are the final state and output of the PLTV system in (18).
The solution of (18) can be expressed in terms of the state transition matrix
as:

$\displaystyle x(h)$
$\displaystyle=\Phi(h,0)x_{0}+\int_{0}^{h}\Phi(h,\tau)B(\tau)w(\tau)\,d\tau$
(48) $\displaystyle z(t)$
$\displaystyle=C(t)\Phi(t,0)x_{0}+\int_{0}^{t}C(t)\Phi(t,\tau)B(\tau)w(\tau)\,d\tau$
(49)

Use these relations to substitute for $x(h)$ and $z(t)$ in Equation (47). Re-
arrange terms to obtain the following form:

$\displaystyle\langle(\hat{x}_{h},\hat{z}),\mathbf{T}_{G}(x_{0},w)\rangle=\hat{x}(0)^{*}x(h)+\int_{0}^{h}\left[B^{*}(\tau)\hat{x}(\tau)+D^{*}(\tau)\hat{z}(\tau)\right]^{*}w(\tau)\,d\tau$

where we have defined the signal

$\displaystyle\hat{x}(\tau):=\Phi(h,\tau)^{*}\hat{x}_{h}+\int_{\tau}^{h}\Phi(t,\tau)^{*}C(t)^{*}\hat{z}(t)\,dt$
(50)

Using property (c) above, we obtain

$\displaystyle\frac{d\hat{x}(\tau)}{d\tau}=-A^{*}(\tau)\hat{x}(\tau)-C^{*}(\tau)\hat{z}(\tau),~{}~{}\hat{x}(h)=\hat{x}_{h}$

If we also define
$\hat{w}(\tau)=B^{*}(\tau)\hat{x}(\tau)+D^{*}(\tau)\hat{z}(\tau),$ and
$\hat{x}_{0}=\hat{x}(0)$ then we can express the inner product
$\langle(\hat{x}_{h},\hat{z}),\mathbf{T}_{G}(x_{0},w)\rangle$ as

$\displaystyle\hat{x}_{0}^{*}x_{0}+\int_{0}^{h}\hat{w}(\tau)^{*}w(\tau)d\tau=\langle\mathbf{T}_{G}^{*}(\hat{x}_{h},\hat{z}),(x_{0},w)\rangle$

Thus the proof is complete.

## Acknowledgment

This work was supported by the National Science Foundation under Grant No.
NSF-CMMI-1254129 entitled “CAREER: Probabilistic Tools for High Reliability
Monitoring and Control of Wind Farms”. Any opinions, findings, and conclusions
or recommendations expressed in this material are those of the author and do
not necessarily reflect the views of the NSF.

The authors greatly acknowledge the help of Henrik Sandberg for making the
MATLAB code of the numerical example presented in [9] available for the
purpose of this research.

## References

  * [1] C. Scherer and S. Weiland, _Lecture Notes DISC Course on Linear Matrix Inequalities in Control_ , 1999. [Online]. Available: http://www.dcsc.tudelft.nl/ cscherer/lmi/notes99.pdf
  * [2] P. Apkarian and P. Gahinet, “A convex characterization of gain-scheduled $\mathcal{H}_{\infty}$ controllers,” _IEEE Transactions on Automatic Control_ , vol. 40, no. 5, pp. 853–864, 1995.
  * [3] C. W. Scherer, “LPV control and full block multipliers,” _Automatica_ , pp. 361–375, 2001.
  * [4] F. Wu, X. H. Yang, A. Packard, and G. Becker, “Induced $\mathcal{L}_{2}$ norm control for LPV systems with bounded parameter variation rates,” _International Journal of Robust and Nonlinear Control_ , vol. 6, no. 983-998, 1996.
  * [5] F. Wu, “A generalized LPV system analysis and control synthesis framework,” _International Journal of Control_ , vol. 74, no. 7, pp. 745–759, 2001.
  * [6] ——, “Control of linear parameter varying systems,” Ph.D. dissertation, University of California at Berkeley, 1995.
  * [7] W. J. Rugh and J. S. Shamma, “Research on gain scheduling,” _Automatica_ , vol. 36, pp. 1401–1425, 2000.
  * [8] J. S. Shamma, _Controls Handbook_. CRC Press, 1996, vol. 1, ch. 20.3. Linearization and gain-scheduling.
  * [9] M. Cantoni and H. Sandberg, “Computing the $\mathcal{L}_{2}$ gain for linear periodic continuous-time systems,” _Automatica_ , vol. 45, pp. 783–789, 2009\.
  * [10] P. Colaneri, “Continuous-time periodic systems in $\mathcal{H}_{2}$ and $\mathcal{H}_{\infty}$: Part i: Theoretical aspects,” _Kybernetica_ , vol. 36, pp. 211–242, 2000.
  * [11] S. Bittanti and P. Colaneri, _Periodic Systems_ , ser. Communications and Control Engineering. Springer, 2009.
  * [12] B. A. Bamieh and J. B. Pearson, “A general framework for linear periodic systems with applications to $\mathcal{H}_{\infty}$ sampled-data control,” _IEEE Transactions on Automatic Control_ , vol. 37, no. 4, pp. 418–435, 1992\.
  * [13] B. A. Bamieh, J. B. Pearson, B. A. Francis, and A. Tannenbaum, “A lifting technique for linear periodic systems with applications to sampled-data control,” _Systems & Control Letters_, vol. 17, pp. 79–88, 1991.
  * [14] K. Zhou, J. C. Doyle, and K. Glover, _Robust and optimal control_. Prentice-Hall, 1996.
  * [15] G. E. Dullerud and S. Lall, “Asynchronous hybrid systems with jumps - Analysis and synthesis methods,” _Systems & Control Letters_, vol. 37, no. 61-69, 1999.
  * [16] A. Neishtadt, _Lectures on Dynamical Systems, Part 1._ , Loughborough University, 2007. [Online]. Available: http://www.iki.rssi.ru/seminar/virtual/an_part1.pdf
  * [17] T. Peni and P. Seiler, “Computation of lower bounds for the induced $\mathcal{L}_{2}$ norm of parameter varying systems,” in _preparation for the American Control Conference (ACC)_ , 2015.
  * [18] H. Pfifer and P. Seiler, “Robustness analysis of linear parameter varying systems using integral quadratic constraints,” _Accepted to the International Journal of Robust and Nonlinear Control_ , 2014.