# Model reduction for linear parameter-varying
systems through parameter projection

Sil Schouten, Daming Lou and Siep Weiland 1 ∗ Accepted by $58^{\text{th}}$
IEEE Conference on Decision and Control.1 Department of Electrical
Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands.
Emails: {S.Schouten, D.Lou<EMAIL_ADDRESS>
###### Abstract

For affine linear parameter-varying (LPV) systems, this paper develops two
parameter reduction methods for reducing the dimension of the parameter space.
The first method achieves the complexity reduction by transforming the affine
LPV system into a parameter-ordered form and establishing an affine upper
bound of the system Gramians, which is extended to time-varying rate-bounded
parameters. The second method is based on considering the sensitivity function
of the transfer function and time evolution equations. Both methods are
applied to an academic example and a thermal model. Simulation results
together with some analysis are given.

##  I Introduction

In the past decade, the class of linear parameter-varying (LPV) systems has
been developed and established as a reliable and very efficient model class
for characterizing nonlinear systems, representing parametric uncertainty and
gain scheduling purposes. Many successful applications ranging from very-
large-scale integration (VLSI) to aircraft designs[1, 2, 3] have been based on
implementations of the LPV framework. The inherent complex nature of physical
systems often results in high dimensional models with large dimensional state
spaces and large dimensional parameter spaces. Typically, the dimension of the
parameter space grows as the complexity of the system increases. In practice,
it is often necessary to evaluate system performance over substantial ranges
of parameter values. We find this theme in problems where design parameters
need to be tuned, calibration problems, geometrical optimisation in circuits
[4] and MEMS devices [5] and robustness analyses in control systems. To have a
reasonable computational complexity in terms of synthesis and simulation,
model order reduction for parametrised systems (pMOR) is often required.

In the aforementioned applications, a high dimensional state space model is
often derived from a high resolution spatial discretization of partial
differential equations (PDEs). Typically, if high precision is required, this
process results in many first order ordinary differential equations (ODEs)
approximating the solutions of the PDE. Moreover, the physical constraints and
design parameters of such high-fidelity models lead to large dimensional
parameter and state spaces. A number of relevant approximation problems can
then be phrased as follows:

  1. 1.

the _state reduction problem_ involves the reduction of the dimension of the
state space, while preserving accuracy and the physical meaning of the
parameters.

  2. 2.

the _state and parameter reduction problem_ involves the simultaneous
reduction of the dimension of the state space and the dimension of the
parameter space.

  3. 3.

the _parameter reduction problem_ involves the reduction of the dimension of
the parameter space only.

In most pMOR work, the primary goal has been to solve the _state reduction
problem_. A number of such methods have been proposed, mainly focusing on
sampling techniques [6, 7, 8]. On this topic, two influential survey papers
[7, 1] have appeared. Methods based on moment matching of the parametrised
transfer function are of particular interest, but tend to become complex with
the number of matched moments. For example, in [7], the dimension of the
reduced model increases exponentially as the number of parameter interpolation
points and moments increase. One solution to this problem is to interpolate
both the state space and the parameter state within the predefined variation
range [9]. Even though it is not always stated, simple sampling schemes imply
static dependence on the parameter. Techniques which explicitly deal with
time-varying parameters are more involved [10]. All these works are confined
to state reduction without considering the problem to reduce the number of
parameters. Approaches which consider parameter reduction either require
typical trajectories of the parameter [11], lack interpretation or are limited
in application [12]. The development of more general parameter reduction
techniques can significantly improve the efficiency of simulations, often
without loss of generality, as was shown in [13]. This leads to the second
challenging problem: the _state and parameter reduction problem_. In [13] a
two-step approach is introduced. First, parameter reduction is employed to
find a low-dimensional parameter space. The second step amounts to
constructing a state reduction via moment-matching. However, the reduced rank
regression method used in the first step only quantifies the relation between
the parameters and the outputs, which is limited by the type of input
excitation used. Besides, the system dynamics are not taken into account over
ranges of parameters. Given the fact that parameter spaces are usually
determined by the physics and the design constraints, it is relevant to
explore the correlation between the parameter space and system theoretical
properties such as reachability and observability of parametrised systems.

Our intention, therefore, is to focus on the _parameter reduction problem_.
Firstly, we exploit the relationship between the parameter space and system
Gramians. A projection-based method is proposed as the means to reduce the
dimension of the parameter space using Hankel-norm approximation. Secondly, we
give an analysis of the sensitivity of the evolution equations to parameter
changes. In doing so, some definitions of system norms for LPV systems are
introduced, aimed at characterizing approximation errors in a consistent
manner. The methods in this paper are developed for time-invariant parameters
with an extension of the Gramian based approach to the time-varying rate-
bounded case.

The paper is organised as follows. In Section II we provide a brief
introduction to LPV systems and introduce the system norms that will be used
in this paper. After a formalisation of the problem, an argument for affine
Gramian reduction is presented. Next, the cross-correlation among the
parameter space in both frequency domain and time domain is developed. Section
IV provides a sensitivity analysis of the parameters on the transfer function.
In Section V, the results are illustrated in an academic example and in a real
application of a thermal model which consists of several interconnected
components. Conclusions and future work are presented in Section VI.

##  II Preliminaries and notation

Consider a system $\Sigma(\theta)$ defined by the following LPV state-space
representation

$\displaystyle\dot{x}(t)$ $\displaystyle=A(\theta)x(t)+B(\theta)u(t)$ (1a)
$\displaystyle y(t)$ $\displaystyle=C(\theta)x(t)+D(\theta)u(t)$ (1b)

where $x\in\mathbb{R}^{n}$ is the state variable, $u\in\mathbb{R}^{m}$ and
$y\in\mathbb{R}^{q}$ denote the input and output, respectively. Furthermore,
$\theta$ represents the parameter that is assumed to reside in the parameter
space $\Theta\subseteq\mathbb{R}^{\ell}$. A distinction is made between
systems where $\theta$ is time-variant and time-invariant. In this work, we
consider $\dot{\theta}=0$ unless stated otherwise. The state space matrices
are assumed to have affine dependence on $\theta$, i.e., the system matrix
$A(\theta)$ is given by

$A(\theta)=A_{0}+A_{1}\theta_{1}+\cdots+A_{\ell}\theta_{\ell}.$ (2)

To exploit the parameter dependence of the system, we rewrite the system
matrices as linear fractional representations, i.e., $A(\theta)$ is
represented as

$A(\theta)=\bar{\theta}_{n}A^{\theta}$ (3)

where

$\bar{\theta}_{n}=[1,\theta_{1},\cdots,\theta_{\ell}]\otimes
I_{n}\in\mathbb{R}^{n\times n(\ell+1)}$ (4)

and

$A^{\theta}=\begin{bmatrix}A_{0}\\\ \vdots\\\
A_{\ell}\end{bmatrix}\in\mathbb{R}^{n(\ell+1)\times n}$ (5)

is the matrix of column stacked matrices. All superscripted matrices represent
stacked matrices obtained in this manner. Here $\otimes$ denotes tensor
multiplication. Similar expressions apply to $B({\theta})$, $C({\theta})$ and
$D({\theta})$. Furthermore, the parameter space $\Theta$ is assumed to be the
convex hull of $k$ generating vectors $w_{j}$, $j=1,\ldots,k$, so that
$\Theta=\text{Co}\left\\{w_{1},\dots w_{k}\right\\}$ represents the parameter
space. Without loss of generality, the system is assumed to be scaled such
that $\theta_{i}\in[0,1],i=1,\dots,\ell$. It is also assumed that the LPV
system is quadratically stable, observable and reachable for all
$\theta\in\Theta$. Using the above notation the continuous time affine LPV
system can be rewritten as

$\displaystyle\dot{x}(t)$
$\displaystyle=\bar{\theta}_{n}A^{\theta}x(t)+\bar{\theta}_{n}B^{\theta}u(t)$
(6a) $\displaystyle y(t)$
$\displaystyle=\bar{\theta}_{n}C^{\theta}x(t)+\bar{\theta}_{n}D^{\theta}u(t).$
(6b)

The parameter projection methods proposed in this paper are defined as linear
transformations $T\in\mathbb{R}^{(\ell+1)\times(\ell+1)}$ of the original
parameter space. Specifically, let the transformed parameter $\tilde{\theta}$
be defined as

$\displaystyle\tilde{\theta}:=[1,\theta_{1},\cdots,\theta_{\ell}]\cdot T$ (7)

where $T$ has full rank and satisfies the orthonormality property

$\displaystyle TT^{T}=T^{T}T=I_{(\ell+1)}.$ (8)

Any such $T$ defines the state space matrices of a transformed system
according to

$\displaystyle{A}(\tilde{\theta})$ $\displaystyle=(\tilde{\theta}\otimes
I_{n})(T^{T}\otimes I_{n})A^{\theta}$ (9a) $\displaystyle{B}(\tilde{\theta})$
$\displaystyle=(\tilde{\theta}\otimes I_{n})(T^{T}\otimes I_{n})B^{\theta}$
(9b) $\displaystyle{C}(\tilde{\theta})$ $\displaystyle=(\tilde{\theta}\otimes
I_{n})(T^{T}\otimes I_{n})C^{\theta}$ (9c) $\displaystyle{D}(\tilde{\theta})$
$\displaystyle=(\tilde{\theta}\otimes I_{n})(T^{T}\otimes I_{n})D^{\theta}$
(9d)

where $(T\otimes I_{n})$ denotes the block-diagonal matrix
$\text{diag}(T,\ldots,T)$. Then it is easily seen that, by (8), the condition

$\displaystyle(\tilde{\theta}\otimes I_{n})(T^{T}\otimes I_{n})$
$\displaystyle=([1,\theta_{1},\cdots,\theta_{\ell}]\cdot T\otimes
I_{n})(T^{T}\otimes I_{n})$
$\displaystyle=([1,\theta_{1},\cdots,\theta_{\ell}]\otimes
I_{n})(TT^{T}\otimes I_{n})$
$\displaystyle=[1,\theta_{1},\cdots,\theta_{\ell}]\otimes I_{n}$ (10)

guarantees system equivalence in the sense that
${A}(\tilde{\theta})=A(\theta$), ${B}(\tilde{\theta})=B(\theta)$,
${C}(\tilde{\theta})=C(\theta$), ${D}(\tilde{\theta})=D(\theta)$.

After introducing transformation of the parameter space, a reduction of the
parameter space is achieved by choosing $T_{r}\in\mathbb{R}^{(\ell+1)\times
n_{r}}$ with rank $n_{r}<\ell$ and by setting

$\tilde{\theta}_{r}:=[1,\theta_{1},\cdots,\theta_{\ell}]\cdot T_{r}$ (11)

as the lower dimensional parameter vector. In this case, $T_{r}$ no longer
satisfies (8) but defines an orthonormal projection whenever
$T_{r}^{T}T_{r}=I_{n_{r}}$. By applying the parameter projection $T_{r}$
defined above, the transformation (9) results in a system
${\Sigma}(\tilde{\theta}_{r})$ with $n_{r}<\ell$ parameters. The state space
matrices of the reduced system ${\Sigma}(\tilde{\theta}_{r})$ are

$\displaystyle{A}(\tilde{\theta}_{r})$
$\displaystyle=(\tilde{\theta}_{r}\otimes I_{n})(T_{r}^{T}\otimes
I_{n})A^{\theta}$ (12a) $\displaystyle{B}(\tilde{\theta}_{r})$
$\displaystyle=(\tilde{\theta}_{r}\otimes I_{n})(T_{r}^{T}\otimes
I_{n})B^{\theta}$ (12b) $\displaystyle{C}(\tilde{\theta}_{r})$
$\displaystyle=(\tilde{\theta}_{r}\otimes I_{n})(T_{r}^{T}\otimes
I_{n})C^{\theta}$ (12c) $\displaystyle{D}(\tilde{\theta}_{r})$
$\displaystyle=(\tilde{\theta}_{r}\otimes I_{n})(T_{r}^{T}\otimes
I_{n})D^{\theta}.$ (12d)

Noticing that the reduced system matrices have the same rank as the matrices
of the original system. To evaluate the performance of the reduced system, an
error system is defined as the LPV system with input $u$ and output $y-y_{r}$
with $y$ and $y_{r}$ the outputs of (1) and (12). See Fig. 1. This system is
compactly denoted as
$\Sigma_{e}(\theta)=\Sigma(\theta)-{\Sigma}(\tilde{\theta}_{r})$ where
$\tilde{\theta}_{r}$ is defined in (11). As such, $\Sigma_{e}(\theta)$ is
viewed as an error system in the parameter $\theta$ only. Note that the error
system is affinely dependent on the parameters.

$1/s$${\Sigma}(\theta)$${\Sigma}(\tilde{\theta}_{r})$$1/s$$+$$-$$x_{r}$$\dot{x}_{r}$$\dot{x}$$x$$u$$y_{r}$$y$$e$$\Sigma_{e}(\theta)$
Figure 1: Interconnection of the reduced parameter error system, indicated by
the blue area.

In the LTI setting, the error system can be evaluated by many well established
and computable norms. Among these, the $H_{\infty}$-, $H_{2}$\- and Hankel-
norm are commonly used in model reduction. In this work, we introduce a
composite error measure on the LPV system that consists of a norm over the
parameter space and a system norm

$\left\lVert\Sigma_{e}(\theta)\right\rVert_{p_{\infty,H}}:=\underset{\theta\in\Theta}{\max}\left\lVert\Sigma_{e}(\theta)\right\rVert_{H},$
(13)

which evaluates the maximal Hankel-norm of the system $\Sigma(\theta)$ when
ranging over the feasible parameter space. Based on what we have discussed so
far, we give the problem formulation of parameter reduction for LPV systems.

Problem: Given $\Sigma(\theta)$ such as (1) and $n_{r}<\ell$, find
$T_{r}\in\mathbb{R}^{(\ell+1)\times n_{r}}$ and $\Sigma(\tilde{\theta}_{r})$
such that

$\left\lVert\Sigma_{e}(\theta)\right\rVert_{p_{\infty,H}}=\underset{\theta\in\Theta}{\max}\left\lVert\Sigma(\theta)-{\Sigma}(\tilde{\theta}_{r})\right\rVert_{H}$
(14)

is minimal.

##  III Hankel-norm reduction

It is well known that the Hankel norm of a stable LTI system can be expressed
in terms of reachability and obervability Gramians [14]. We establish a
similar result for LPV systems first. The Hankel operator associated with a
stable LTI system $\Sigma$ is defined by

$\displaystyle\mathcal{H}:\mathcal{L}_{2}(\mathbb{R}_{-}^{m})\longmapsto\mathcal{L}_{2}(\mathbb{R}_{+}^{q}),u_{-}\longmapsto
y_{+}$
$\displaystyle\text{where}\quad\mathcal{H}(u_{-})(t)=\int_{-\infty}^{0}H(t-\tau)u(\tau)d\tau,t\in\mathbb{R}_{+},$
(15)

it maps the past inputs $u_{-}$ into future outputs $y_{+}$. Here, $H$ is the
impulse response of $\Sigma$. The $\mathcal{\ell}_{2}$-induced norm of
$\mathcal{H}$ is defined as

$\displaystyle\left\lVert\mathcal{H}\right\rVert_{\mathcal{L}_{2}-\text{ind}}=\sup_{\left\lVert
u_{-}\right\rVert_{2}}\frac{\left\lVert y_{+}\right\rVert_{2}}{\left\lVert
u_{-}\right\rVert_{2}}=\left\lVert\Sigma\right\rVert_{H}.$ (16)

The quantity $||\Sigma||_{H}$ is the Hankel norm of the system $\Sigma$ and
equals the spectral norm $||\Sigma||_{H}=\sigma_{\text{max}}(\mathcal{H})$.
The Hankel singular values of the system $\Sigma$ are defined as the singular
values of the Hankel operator $\mathcal{H}$ associated with $\Sigma$.

###### Lemma 1

Given a reachable, observable and stable LTI system $\Sigma$ of dimension $n$,
the Hankel singular values are equal to the absolute value of the eigenvalues
of the product of $\mathcal{PQ}^{T}$

$\displaystyle\sigma_{i}(\Sigma)=\sqrt{\lambda_{i}(\mathcal{PQ}^{T})},i=1,...,n$
(17)

where $\mathcal{P}$ and $\mathcal{Q}$ are the controability Gramian and
observability Gramian of $\Sigma$.

For a time-invariant LPV system, the Hankel operator is defined as

$\displaystyle\mathcal{H}_{\theta}:\mathcal{L}_{2}(\mathbb{R}^{m}_{-}\times\mathbb{P})\longmapsto\mathcal{L}_{2}(\mathbb{R}^{q}_{+}),\
u_{-},\theta\longmapsto y_{+}$ where (18)
$\displaystyle\quad\mathcal{H}_{\theta}(\theta,u_{-})(t)=\int_{-\infty}^{0}H(\theta,t-\tau)u(\tau)d\tau,t\in\mathbb{R}_{+},\theta\in\Theta,$

it maps the past inputs and parameters into future outputs. It follows that
the $\ell_{2}$-induced norm of $\mathcal{H}_{\theta}$ is now parameter
dependent. Since $\Sigma(\theta)$ is assumed to be stable
$\forall\theta\in\Theta$, the Hankel norm of $\Sigma(\theta)$ can be expressed
as

$\displaystyle\left\lVert\Sigma(\theta)\right\rVert_{H}=\sigma_{\text{max}}(\mathcal{H}_{\theta})=\sqrt{\lambda_{\text{max}}\left(\mathcal{P}(\theta)\mathcal{Q}^{T}(\theta)\right)},$
(19)

where $\mathcal{P}(\theta),\mathcal{Q}(\theta)$ are the reachability and
observability gramians, defined as the unique solutions of the parametrized
Lyapunov equations

$\displaystyle
A(\theta)\mathcal{P}(\theta)+\mathcal{P}(\theta)A^{T}(\theta)+B(\theta)B^{T}(\theta)$
$\displaystyle=0$ (20) $\displaystyle
A(\theta)^{T}\mathcal{Q}(\theta)+\mathcal{Q}(\theta)A(\theta)+C(\theta)^{T}C(\theta)$
$\displaystyle=0.$ (21)

Finding an exact solution to the Lyapunov equation for the whole parameter
space is not trivial and often intractable. In the literature, a static
Gramian [15] is proposed which often suffices, but leads to conservative
solutions. For parameter dimension reductions, the parameter dependent
Gramians are necessary as they express changes in the system due to parameter
variations.

The following result shows that a relaxation of (20) to an inequality
naturally leads to an upper bound on $\mathcal{P}(\theta)$ and
$\mathcal{Q}(\theta)$ for all $\theta$. Subsequently, an upper bound of the
Hankel norm (19) is derived.

###### Theorem 1

Consider an LPV system (1) that is stable, reachable and observable for all
$\theta\in\Theta$. There exist unique solutions
$\mathcal{P}(\theta)=\mathcal{P}^{T}(\theta)\succeq
0,\mathcal{Q}(\theta)=\mathcal{Q}^{T}(\theta)\succeq 0$ which satisfies (20)
and (21). Furthermore, there exists $P(\theta)=P^{T}(\theta)\succeq 0$ and
$Q(\theta)=Q^{T}(\theta)\succeq 0$ which satisfy

$\displaystyle
A(\theta)P(\theta)+P(\theta)A^{T}(\theta)+B(\theta)B^{T}(\theta)$
$\displaystyle\preceq 0.$ (22) $\displaystyle
A(\theta)^{T}{Q}(\theta)+{Q}(\theta)A(\theta)+C(\theta)^{T}C(\theta)$
$\displaystyle\preceq 0.$ (23)

for all $\theta\in\Theta$. All solutions $P(\theta)$ an $Q(\theta)$ upper
bound $\mathcal{P}(\theta)$ and $\mathcal{Q}(\theta)$ in the sense that
$P(\theta)\succeq\mathcal{P}(\theta)\succeq 0$ and
$Q(\theta)\succeq\mathcal{Q}(\theta)\succeq 0$. Moreover, the Hankel norm of
$\Sigma(\theta)$ is upper bounded by

$\displaystyle\lVert\Sigma(\theta)\rVert_{H}\\!=\sigma_{\text{max}}(\mathcal{H}_{\theta})$
$\displaystyle=\\!\sqrt{\\!\lambda_{\text{max}}\left(\mathcal{P}(\theta)\mathcal{Q}^{T}(\theta)\right)}\\!$
$\displaystyle\leq\\!\sqrt{\\!\lambda_{\text{max}}\left(P(\theta)Q^{T}(\theta)\right)}.$
(24)

The proof is given in Appendix A. Note that solutions of (22) and (23) are not
unique. To measure the importance of the parameter $\theta\in\Theta$ to the
system Gramians, a function which is affinely dependent on $\theta$ is of
particular interest to find an upper bound of $\mathcal{P}(\theta)$ and
$\mathcal{Q}(\theta)$ for all $\theta\in\Theta$.

###### Theorem 2

Consider an LPV system defined by (1) that is stable, reachable and observable
for all $\theta\in\Theta$. Suppose that
$\Theta=\text{Co}\\{w_{1},\dots,w_{k}\\}$. Then, for $\mathcal{P}(\theta)$
which satisfies (20) there always exists an affine function
$f:\Theta\rightarrow\mathbb{R}^{n\times n}$ such that

$f(\theta)\succeq\mathcal{P}(\theta)\text{ for all }\theta\in\Theta.$ (25)

Here $f(\theta)=\mathbf{\bar{\theta}}_{n}P^{\theta}$ and $P^{\theta}$ denote
stacked matrices $P(\theta)$ which satisfy (22).

_Theorem 2_ can be proved by making use of the convexity of affine functions
$f(\theta)$ with convex hull of $\theta\in\Theta$. The upper bound affine
function for observability Gramian is omitted for clarity.

###  III-A Parameter reduction using Hankel-norm approximation

The main ingredients of parameter reduction for LPV systems are to construct
an upper bounded Hankel norm defined in _Theorem 1_ which satisfies the affine
form in _Theorem 2_, and then to find a reduced model with the least error for
the given reduced parameter space. By choosing affine upper bounding Gramians
of $\Sigma(\theta)$ and using the transformation introduced previously, the
following equivalence can be obtained

$P(\theta)Q^{T}(\theta)=\bar{\theta}_{n}P^{\theta}Q^{\theta^{T}}\bar{\theta}_{n}^{T}.$
(26)

Following the same arguments (9) as for
$P(\bar{\theta})=([1,\theta_{1},\cdots,\theta_{\ell}]\cdot T\otimes
I_{n})(T^{T}\otimes I_{n})P^{\theta}$, pick up an appropriate
$T_{r}\in\mathbb{R}^{(\ell+1)\times n_{r}}$ with $n_{r}<\ell$, use the
relation in (II) and the equivalence is found

$\small{P(\tilde{\theta}_{r})Q^{T}(\tilde{\theta}_{r})=\bar{\theta}_{n}(T_{r}T_{r}^{T}\otimes
I_{n})P^{\theta}Q^{\theta^{T}}(T_{r}T_{r}^{T}\otimes
I_{n})\bar{\theta}_{n}^{T}}.$ (27)

Once the upper bounded Hankel norm of the original model and the reduced model
has been defined and transformed, the problem stated in _Section II_ can be
reformulated as

$\displaystyle\underset{\tiny{T_{r}}}{\text{min}}$
$\displaystyle\underset{\theta\in\Theta}{\text{max}}\lVert
P(\theta)Q^{T}(\theta)-P(\tilde{\theta}_{r})Q^{T}(\tilde{\theta}_{r})\rVert_{H}$
(28) s.t $\displaystyle
P(\theta)Q^{T}(\theta)=\bar{\theta}_{n}P^{\theta}Q^{\theta^{T}}\bar{\theta}_{n}^{T},$
$\displaystyle
P^{\theta^{T}}\mathbf{\bar{\theta}}^{T}_{n}=\mathbf{\bar{\theta}}_{n}P^{\theta}\succeq
0,\quad Q^{\theta^{T}}\bar{\theta}_{n}^{T}=\bar{\theta}_{n}Q^{\theta}\succeq
0$ $\displaystyle T_{r}\in\mathbb{R}^{(\ell+1)\times n_{r}}.$

Remark that the explicit expression of
$P(\tilde{\theta}_{r})Q^{T}(\tilde{\theta}_{r})$ is given in and it is $T_{r}$
dependent. For solving the above $\min-\max$ optimisation problem, the first
step is to find an upper bound $P(\theta)$ that satisfies _Theorem 1_ and
_Theorem 2_. This is equivalent to finding a solution of the following
inequalities

$\displaystyle\mathbf{\bar{\theta}}_{n}\left[{A}^{\theta}P^{\theta^{T}}+P^{\theta}{A^{\theta}}^{T}+{B^{\theta}}{B^{\theta}}^{T}\right]\mathbf{\bar{\theta}}^{T}_{n}\preceq
0$ , (29a) $\displaystyle
P^{\theta^{T}}\mathbf{\bar{\theta}}^{T}_{n}=\mathbf{\bar{\theta}}_{n}P^{\theta}\succeq
0$ . (29b)

The second step is to construct a $T_{r}$ which minimises the error between
the original system and the reduced one with the maximal $\theta$. This
discussion is summarised in _Algorithm 1_. Indeed, the ’optimal’ solution of
(28) is not guaranteed since it is not convex problem.

Algorithm 1 : Parameter reduction for affine LPV systems

1:Construct parameter affine dependent upper bounds
$\bar{\theta}_{n}P^{\theta},\bar{\theta}_{n}Q^{\theta}$ that satisfy _Theorem
1_ and _Theorem 2_.

2:Obtain $v_{n_{r}}$ by solving (28) for given $n_{r}$,

3:Orthonormalise the $v_{n_{r}}$ such that $v_{n_{r}}^{T}v_{n_{r}}=I_{n_{r}}$.

4:Assign $T_{n_{r}}\leftarrow v_{n_{r}}$.

5:Return $T_{n_{r}}$

The optimum of the quadratic function (29) of $\theta$ can be evaluated on the
vertices of $\Theta$. A formal proof of this observation is given in [16].
Finding a solution to this Lyapunov inequality is done using linear matrix
inequalities (LMIs). Every vertex of the parameter space results in two LMIs
(29). In the case of interval bounding constraints on each entry of the
parameter vector, this results in a total of $2^{(\ell+1)}$ LMIs. Like many
problems in the LPV setting, this approach becomes intractable with growing
dimensions of the parameter space. The difference with other problems is that
there is no tuning involved in this process. Therefore there is a one time
cost of solving the LMIs.

The previously presented method assumed static parameter dependence
$\dot{\theta}=0$. In the case where the parameter is dynamic but is rate
limited $\dot{\underline{\theta}}_{i}<\dot{\theta}_{i}<\dot{\bar{\theta}}_{i}$
for every $\theta_{i}\subseteq{\theta}$, the Lyapunov inequality may be
adjusted according to (LABEL:eq:timevariantlyapunov). The resulting LMIs are
again quadratic in $\theta$, now with a constant offset. Therefore it can
still be evaluated on the vertices of the parameter space. The formal proof is
again found in [16].

$\displaystyle\dot{\mathcal{P}}(\theta)+{A}(\theta)\mathcal{P}(\theta)+\mathcal{P}(\theta)A(\theta)^{T}+{B}(\theta){B}(\theta)^{T}\preceq
0$ (30a)
$\displaystyle{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\small
AffineGramian}\Downarrow\text{\small Affine Gramian}$ $\displaystyle
P^{\theta}\dot{\theta}+\mathbf{\bar{\theta}}_{n}\left[{A}^{\theta}P^{\theta^{T}}+P^{\theta}{A^{\theta}}^{T}+{B^{\theta}}{B^{\theta}}^{T}\right]\mathbf{\bar{\theta}}^{T}_{n}\preceq
0$ (30b)
$\displaystyle{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\small
Boundedparametervelocity}\Downarrow\text{\small Bounded parameter velocity}$

$L(\theta,\dot{\underline{\theta}})=P^{\theta}\mathbf{\dot{\underline{\theta}}}_{n}+\mathbf{\bar{\theta}}_{n}\left[{A}^{\theta}P^{\theta^{T}}+P^{\theta}{A^{\theta}}^{T}+{B^{\theta}}{B^{\theta}}^{T}\right]\mathbf{\bar{\theta}}^{T}_{n}\preceq
0$

(30c)

$L(\theta,\dot{\overline{\theta}})=P^{\theta}\mathbf{\dot{\bar{\theta}}}_{n}+\mathbf{\bar{\theta}}_{n}\left[{A}^{\theta}P^{\theta^{T}}+P^{\theta}{A^{\theta}}^{T}+{B^{\theta}}{B^{\theta}}^{T}\right]\mathbf{\bar{\theta}}^{T}_{n}\preceq
0$

(30d)

Extending the approach to parameters which are rate bounded increases the
number of LMIs to $3^{\ell+1}$. By imposing an additional constraint
$P_{i}^{\theta}\succ 0$ (We abuse the notation $P_{i}^{\theta}$ for $i$th
stacked matrix of $P^{\theta}$ which is referred to the notation in (5).)
results in $L(\theta,\dot{\underline{\theta}})\prec
L(\theta,\dot{\overline{\theta}})$ and thus (30c) can be dropped. This reduces
the number of LMIs to $2^{\ell+1}+\ell$. The previously presented methods to
obtain an affine reachability Gramian can be applied for finding an affine
observability Gramian $Q(\theta)$ as well.

In LTI approximation the loss function is related to the energy of the error
system since the Hankel singular values are a measure of energy in each state
respectively. For the proposed method, the interpretation is more nuanced as
the affine Gramians constitute an upper bound on the actual Gramians. If this
upper bound is tight, the loss function is a good indication of the error. The
tightness of this upper bound is heavily dependent on the system, and
determining the tightness is not trivial. Inspired by Adamjan-Arov-Krein (AAK)
theorem, we provide a relative $p_{\infty,\infty}$-error which asses the
approximation error

$\displaystyle\frac{\lVert{\Sigma}(\theta)-\tilde{\Sigma}_{r}(\theta_{r})\rVert_{p_{\infty,\infty}}}{\lVert{\Sigma}(\theta)\rVert_{p_{\infty,\infty}}}.$
(31)

Here $||\cdot||p_{\infty,\infty}$ is defined as

$\|\Sigma_{e}(\theta)\|_{p_{\infty,\infty}}:=\underset{\theta\in\Theta}{\max}\|\Sigma_{e}(\theta)\|_{H_{\infty}}$
(32)

which evaluates the $H_{\infty}$-norm of the system $\Sigma(\theta)$ when
varying over the feasible parameter space.

##  IV Sensitivity analysis

Another method for parameter reduction is sensitivity analysis, which is
similar to principal component analysis [17, 18]. This method evaluates how
sensitive the outputs are to small changes in the parameter values. It can be
viewed as principal component analysis without the requirement of needing a
typical parameter trajectory. Evaluation can be done in either the time or the
frequency domain.

###  IV-A Frequency domain

Since it is assumed that $\dot{\theta}=0$ the notion of transfer functions is
still applicable. The transfer function of an affine LPV system is given below

$\displaystyle
H(\theta,q)=D_{\theta}+C_{\theta}q(I-A_{\theta}q)^{-1}B_{\theta}.$ (33)

The _transfer function sensitivity_ is defined as the Jacobian of the transfer
function with respect to $\theta$, as shown in the following

$\displaystyle{{}_{\theta}}H(\theta,q)$
$\displaystyle=\frac{dD_{\theta}}{d\theta}+\hbox{\pagecolor{blue!10}\large$\frac{dC_{\theta}}{d\theta}q(I-A_{\theta}q)^{-1}B_{\theta}$}$
(34)
$\displaystyle+\hbox{\pagecolor{green!10}\large$C_{\theta}q\frac{d(I-A_{\theta}q)^{-1}}{d\theta}B_{\theta}$}+\hbox{\pagecolor{red!10}\large$C_{\theta}q(I-A_{\theta}q)^{-1}\frac{dB_{\theta}}{d\theta}$}$

The $i^{th}$ element of the Jacobian can be represented as the system given in
(42), where the colours correspond to separate elements of the Jacobian as in
(34). Note that the resulting system again has affine dependency on the
parameters.

(42)

From the Jacobian an ordering of the parameters can be determined from the
$p_{\infty,\infty}$-norm of each respective element. However, the transfer
function sensitivity does not take into account the cross-correlations between
parameters. Multiplying the transfer function sensitivity by its complex
conjugate will result in a transfer function sensitivity covariance matrix
(TSCM). The elements of this matrix are defined as

$\displaystyle\Pi_{ij}:=\left\lVert{{{}^{i}_{\theta}}H(\theta,q)}^{*}{{{}_{\theta}^{j}}H(\theta,q)}\right\rVert_{p_{\infty,\infty}}.$
(43)

This product of systems gives element $i,j$ of the TSCM, in the notation of
(42), and is equal to

$\displaystyle\Pi_{i,j}=\left\lVert\left[\begin{array}[]{cc|c}\mathcal{A}^{*}_{\theta,i}&\mathcal{C}^{*}_{\theta,i}\mathcal{C}_{\theta,j}&\mathcal{C}^{*}_{\theta,i}D_{j}\\\
0&\mathcal{A}_{\theta,j}&\mathcal{B}_{\theta,j}\\\
\hline\cr\mathcal{B}^{*}_{\theta,i}&D^{*}_{i}\mathcal{C}_{\theta,j}&D_{i}^{*}D_{j}\end{array}\right]\right\rVert_{p_{\infty,\infty}}$
(47)

All elements of this system have affine dependency on $\theta$, except for
$\mathcal{C}^{*}_{\theta,i}\mathcal{C}_{\theta,j}$ which in general does not
possess such a property. Still the $p_{\infty,\infty}$-norm can be evaluated
over a finite set of points. This can be shown by extending the parameter
space with all quadratic elements of $\theta$. The system will then have
affine dependence on the extended parameter space and thus can be evaluated
over an extended convex hull. If the extended convex hull is not properly
chosen it may lead to conservatism.

The resulting TSCM is a symmetric $n_{\theta}\times n_{\theta}$ matrix. By
taking the singular value decomposition (SVD) of this matrix, $\Pi=TST^{*}$,
the parameter transformation matrix is found. This transformation orders the
parameter directions in terms of transfer sensitivity covariance and is
orthonormal, constituting a valid transformation as defined in (7).

###  IV-B Time domain

Unlike principal component analysis, the proposed time-domain approach does
not require a typical trajectory of the parameters. Instead the solution is
written as time dependent equation. In discrete time the output evolution
equation is given as

(48)

The sensitivity function at time $k$ is simply calculated as the Jacobian of
this equation towards $\theta$. For affine parameter dependency, (49) can be
calculated in terms of the stacked matrices $\bullet^{\theta}$ and state space
matrices $\bullet_{\theta}$.

$\displaystyle{{}_{\theta}}y(k)=$
$\displaystyle\frac{dC_{\theta}A_{\theta}}{d\theta}x_{0}+\frac{dD_{\theta}}{d\theta}u(k)+$
(49)
$\displaystyle\left[\sum_{i=1}^{k}\frac{dC_{\theta}A_{\theta}^{i-1}B_{\theta}}{d\theta}u(k-i)\right]$

In discrete time the sensitivity can be written as a row vector multiplied by
a column vector stacking all inputs. Stacking all outputs,
${{}_{\theta}}y_{0},\ {{}_{\theta}}y_{1}\ \dots\ {{}_{\theta}}y_{k_{\max}}$,
into a column vector gives a matrix multiplied by the input vector
${{}^{i}_{\theta}}Y={{}^{i}}M(\theta)U$. In the particular case where
$A_{\theta}=A_{0}$ and either $C_{\theta}=C_{0}$ or $B_{\theta}=B_{0}$, the
matrix ${{}^{i}}M(\theta)$ is parameter independent. To take cross correlation
into account ${{}^{i}_{\theta}}Y$ is pre multiplied by
${{}^{j}_{\theta}}Y^{*}$, being equal to
${{}^{j}_{\theta}}Y^{*}{{}^{i}_{\theta}}Y=U^{*}{{}^{j}}M^{*}{{}^{i}}MU$. The
largest singular values of ${{}^{j}}M^{*}{{}^{i}}M$ result in the peak gain
between the sensitivity functions of the $i^{th}$ and $j^{th}$ parameter. The
collection matrix of all the maximum singular values $\bar{\sigma}_{ij}$ gives
an $n_{p}\times n_{p}$ matrix, the sensitivity covariance matrix (SCM).

$\displaystyle\bar{S}=\begin{bmatrix}\bar{\sigma}_{11}&&\bar{\sigma}_{12}&&\dots&&\bar{\sigma}_{1n_{\theta}}\\\
\bar{\sigma}_{21}&&\bar{\sigma}_{22}&&\dots&&\bar{\sigma}_{2n_{\theta}}\\\
\vdots&&\vdots&&\ddots&\vdots\\\
\bar{\sigma}_{n_{\theta}1}&&\bar{\sigma}_{n_{\theta}2}&&\dots&&\bar{\sigma}_{n_{\theta}n_{\theta}}\end{bmatrix}$
(50)

From the SCM a similar approach is taken to derive the transformation matrix
as with the TSCM. One note on this approach is the choice of $k_{\max}$ which
should be chosen at least as large as the largest time constant of the system.
Choosing $k_{\max}$ small will exclude dynamics and may lead to bad
approximations.

##  V Results

To illustrate these methods of parameter reduction through sensitivity
analysis and Hankel-norm approximation, two examples systems are used. The
first system is an illustrated LPV system and the second a real thermal
system.

###  V-A Illustrative example

Consider an affine LPV system

$\displaystyle\dot{x}$
$\displaystyle=A_{0}x+B_{0}u+\sum_{i=1}^{5}\left(A_{i}\theta_{i}x+B_{i}\theta_{i}u\right)$
(51a) $\displaystyle y$
$\displaystyle=C_{0}x+D_{0}u+\sum_{i=1}^{5}\left(C_{i}\theta_{i}x+D_{i}\theta_{i}u\right),$
(51b)

where $\theta_{i}\in[0,\ 1],u\in\mathbb{R}^{2}$, $y\in\mathbb{R}^{2}$ and
$x\in\mathbb{R}^{45}$. To ensure stability all $A_{i}\prec 0$. For this system
Hankel-norm and TSCM approximations are determined. As a comparison subsystem
Hankel-norm approximation is used where reduction is based singular values of
actual Hankel norm associated with the remaining parameter space. The
$n_{p}-n_{r}$ parameters having the smallest ’subsystem’ Hankel-norm are
removed. In Fig. 2 the output evolution which consists of the transient
response and the steady-state response is shown for an input $u(t)=$ constant
$\neq 0$ for $t=[0,25]$ and $u(t)=0$ otherwise. To illustrate the accuracy of
the reduced models, the errors are shown in Fig. 3. For $t=[0,25]$, both
figures show a tendency that more parameters preserved, the less error is
between the reduced and the original model. For $t>25[s]$, all reduced models
with $n_{r}=1,...,4$ convergence to the same point. It will change for a
different parameter realisation and therefore Fig. 4 shows the relative
$p_{\infty},\infty$ error of the reductions to be non-increasing with growing
parameter order for both methods. It also shows that for this example, the
sensitivity analysis outperforms the Gramian based method.

Figure 2: System evolution of the randomised LPV system, at a random parameter
realisation, for different reduction orders and methods. Figure 3: Output
error of the randomised LPV system, at a random parameter realisation, for
different reduction orders and methods. Figure 4: Relative
$p_{\infty,\infty}$-error of the randomised LPV system for different parameter
orders and approximation methods.

###  V-B Thermal simulation

The second system is a thermal simulation consisting of five coupled metal
blocks, all having parameter dependent heat capacity. Using COMSOL the system
is generated with the following structure.

$\displaystyle\dot{x}$
$\displaystyle=A_{0}x+B_{0}u+\sum_{i=1}^{5}\left(A_{i}\theta_{i}x+B_{i}\theta_{i}u\right)$
(52a) $\displaystyle y$ $\displaystyle=C_{0}x$ (52b)

Where $u\in\mathbb{R}^{2}$, $y\in\mathbb{R}^{2}$ and $x\in\mathbb{R}^{45}$.
The inputs and outputs represent heat power in $[W]$ and temperature $[K]$
respectively. In Fig. 5 an illustration of the system is shown.

Figure 5: Thermal model with five different material blocks. Orange and pink
are heat inputs. Purple and blue are outputs.

To evaluate the error of the reduced model, a simulation is performed with
constant power from $u(t)=[50;45][W]$ for $t\in[0\ 250]$ and $u(t)=[0;0]$
afterwards. For different realisations of the parameter values, the error
between the reduced model and full order model is shown in Fig. 6. Clearly the
non-parametric model $n_{r}=0$ performs the worst showing that the nominal
model is not a good approximation of the system. The best performing model is
the $n_{r}=4$ model as the error is almost zero. This result shows the
performance with respect to a specific parameter space. Therefore, a global
error bound is used to infer conclusions on performance. In Fig. 7 the
relative $p_{\infty,\infty}$ error is plotted for different reduction orders.
This figure illustrates that reducing the system using sensitivity analysis is
comparable to reduction in subsystem using Hankel-norm. It is also clear that
approximation in the $p_{\infty,H}$ norm using transformation optimisation
yields improved results for $n_{r}<4$. For approximation of $n_{r}=4$ the
Hankel-norm optimisation method performs worse in $p_{\infty},\infty$ error.
This is due to a combination of issues, the non-convex optimisation (28) and
the error introduced in finding the upper bound of the affine Gramian (25).

Figure 6: Simulation error of the thermal system for different parameter
orders at a randomly selected parameter in the parameter space. Figure 7:
Relative $p_{\infty},\infty$ error of the thermal system for different
parameter orders and different reduction methods.

For both systems computation of the Gramians takes in the order of
approximately $1000$ seconds(with a dual core PC and using YALMIP [19]).
Because the optimisation in (28) has not shown to be convex, it is ran using
different initial condition to find a close approximation of the Hankel
singular values. Due to this non-convexity, the resulting transformation
matrix does not guarantee the global optimal approximation.

##  VI Conclusions and future work

In this paper, two methods for parameter reduction have been given for time-
invariant LPV systems. The first method approximates the system using the
Hankel-norm over the parameter space. The second method uses principal
component analysis on the covariance matrix. A system norm analysis of the
performance of these methods has been presented in the paper together with
simulation results. It can be concluded that both methods, though
substantially different, are computationally feasible and provide good
approximations over the parameter space.

Two reduction methods have been presented which focus on reducing the
parameter space. The proposed methods are not limited to be only used for
simplifying the complexity of parameters, but also can be integrated with
state reduction problem as state and parameter reduction problem.

## References

  * [1] P. Benner, S. Gugercin, and K. E. Willcox, “A survey of projection-based model reduction methods for parametric dynamical systems,” SIAM Review, vol. 57, no. 4, pp. 483–531, 2015.
  * [2] A. Marcos and G. J. Balas, “Development of Linear Parameter Varying Models for Aircraft,” Journal of Guidance, Control and Dynamics, vol. 27, no. 2, pp. 218–228, 2004.
  * [3] M. Dettori and C. W. Scherer, “LPV design for a CD player : an experimental evalua- tion of performance,” Selected Topics in Signals, Systems and Control, vol. 12, no. September, 2001.
  * [4] L. Daniel, O. C. Siong, and L. S. Chay, “A Multiparameter Moment Matching Model Reduction Approach for Generating Geometrically Parameterized Interconnect Performance Models,” Tcad, vol. 23, no. 5, pp. 1–15, 2004\.
  * [5] U. Baur, P. Benner, A. Greiner, J. G. Korvink, J. Lienemann, and C. Moosmann, “Parameter preserving model order reduction for MEMS applications,” Mathematical and Computer Modelling of Dynamical Systems, vol. 17, no. 4, pp. 297–317, 2011.
  * [6] H. Panzer, J. Mohring, R. Eid, and B. Lohmann, “Parametrische ordnungsreduktion mittels matrixinterpolation,” Automatisierungstechnik, vol. 58, no. 8, pp. 475–484, 2010.
  * [7] P. Benner and L. H. Feng, “A robust algorithm for parametric model order reduction based on implicit moment matching,” in Reduced Order Methods for Modeling and Computational Reduction (A. Quarteroni and G. Rozza, eds.), ch. 6, pp. 159–186, Springer International Publishing Switzerland, 9 ed., 2014\.
  * [8] T. Bui-Thanh, K. E. Willcox, and O. Ghattas, “Model reduction for large-scale systems with high-dimensional parametric input space,” SIAM Journal on Scientific Computing, vol. 30, no. 6, pp. 3270–3288, 2008.
  * [9] D. Lou and S. Weiland, “Parametric model order reduction for large-scale and complex thermal systems,” in 2018 European Control Conference (ECC), pp. 2593–2598, IEEE, 2018.
  * [10] G. D. Wood, Control of Parameter-Dependent Mechanical Systems.  Phd, Cambridge University, 1995.
  * [11] S. Z. Rizvi, J. Mohammadpour, R. Tóth, and N. Meskin, “A kernel-based PCA approach to model reduction of linear parameter-varying systems,” IEEE Transactions on Control Systems Technology, vol. 24, no. 5, pp. 1883–1891, 2016.
  * [12] C. Sun and J. Hahn, “Parameter reduction for stable dynamical systems based on Hankel singular values and sensitivity analysis,” Chemical Engineering Science, vol. 61, no. 16, pp. 5393–5403, 2006.
  * [13] Z. Feng and P. Li, “Performance-Oriented Statistical Parameter Reduction of Parameterized Systems via Reduced Rank Regression *,” in 2006 IEEE/ACM International Conference on Computer Aided Design, pp. 868–875, 2006.
  * [14] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems.  SIAM, 2005.
  * [15] C. W. Scherer and S. Weiland, “Linear Matrix Inequalities in Control,” in Linear Matrix Inequalities in Control (C. Scherer and S. Weiland, eds.), ch. 9, pp. 257–280, 2017.
  * [16] P. B. Cox, S. Weiland, and T. Roland, “Affine Parameter-Dependent Lyapunov Functions for LPV Systems with Affine Dependence,” IEEE Transactions on Automatic Control, pp. 1–8, 2018.
  * [17] A. Kwiatkowski and H. Werner, “LPV control of a 2-DOF robot using parameter reduction,” Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference, CDC-ECC ’05, vol. 2005, pp. 3369–3374, 2005.
  * [18] A. Kwiatkowski and H. Werner, “Parameter reduction for LPV systems via pricipal components analysis,” in 16th IFAC World Congress, pp. 596–601, IFAC, 2005.
  * [19] J. Löfberg, “YALMIP : A toolbox for modeling and optimization in MATLAB,” in Proceedings of the CACSD Conference, (Taipei), pp. 284–289, 2004.
  * [20] R. Bhatia, Matrix analysis, vol. 169.  Springer Science & Business Media, 2013.

##  Appendix A Proof of _Theorem 1_

This proof consists of two parts. In the first part, we prove the Lyapunov
inequality as an upper bound to the Lyapunov equality. The second part we show
that upper bound on the Hankel-Singular value from the upper bounded Gramian
which is derived from the first part.

###  A-A Proof of the Lyapunov inequality as an upper bound to the Lyapunov
equality

Consider a stable LTI system with stable matrix $A$. Lyapunov equation theorem
shows that there is a unique solution $P\succ 0$ to

$\displaystyle A\mathcal{P}+\mathcal{P}A^{T}+Q=0$ (53)

where $Q=BB^{T}\succ 0$. The same system admits a solution space
$P\in\tilde{\mathcal{P}}$, being all solutions that satisfy

$\displaystyle AP+PA^{T}+Q\preceq 0$ (54)

The equality can be rewritten and substituted into the inequality to give

$\displaystyle A(P-\mathcal{P})+(P-\mathcal{P})A^{T}\preceq 0$ (55)

From Lyapunov equation theorem and stable $A$, we have $P-\mathcal{P}\succeq
0$ as the solution of the above inequality. Extending to stable, observable
and reachable time-invariant LPV systems, it suffices to show that for every
$\theta\in\Theta$ a solution can be found to the Lyapunov (in)equality
concluding the proof.  $\square$

###  A-B Proof of the upper bound on the Hankel-singular values from the upper
bounding Gramians

From [20] we take two properties of eigenvalue arithmetic of symmetric
matrices $A,B\in\mathbb{R}^{n_{x}\times n_{x}}$,

$\displaystyle\lambda_{i}(A+B)\geq 0$ $\displaystyle\text{ if
}\lambda_{i}(A)\geq 0\text{ and }\lambda_{i}(B)\geq 0,\ $ (56)
$\displaystyle\lambda_{i}(AB)\geq 0$ $\displaystyle\text{ if
}\lambda_{i}(A)\geq 0\text{ and }\lambda_{i}(B)\geq 0,$

for $i=1,...,n$. For clarity we drop $\theta$, and define a function $F$
below:

$\displaystyle F$ $\displaystyle:=PQ-\mathcal{PQ},$
$\displaystyle=(P-\mathcal{P})(Q-\mathcal{Q})+\mathcal{P}(Q-\mathcal{Q})+(P-\mathcal{P})\mathcal{Q}.$
(57)

Given $P(\theta)\succeq\mathcal{P}(\theta)\succ 0$ and
$Q(\theta)\succeq\mathcal{Q}(\theta)\succ 0$, applying properties (56) and
_Theorem 1_ to (57) yields $\lambda_{i}(F)\geq 0$, $i=1,\dots,n$. With the
properties of Hankel matrix $PQ=Q^{\frac{1}{2}}PQ^{\frac{1}{2}}$, for the
symmetric $F=(Q^{1/2}PQ^{1/2}-\mathcal{Q}^{1/2}P\mathcal{Q}^{1/2})=F^{T}$ it
is proven. Therefore, $PQ\succeq\mathcal{PQ}$ is concluded. Next consider the
eigenvector $x_{1}$ associated to the largest eigenvalue of $\mathcal{PQ}$.
Then the following holds

$x^{T}_{1}\mathcal{PQ}x_{1}\leq
x^{T}_{1}(PQ)x_{1}\Rightarrow\lambda_{1}(\mathcal{PQ})\leq\frac{x^{T}_{1}(PQ)x_{1}}{||x_{1}||_{2}^{2}}.$
(58)

By the definition of eigenvalue decomposition, we have

$\frac{x^{T}_{1}(PQ)x_{1}}{||x_{1}||_{2}^{2}}\leq\sup_{y_{1}}\frac{y^{T}_{1}(PQ)y_{1}}{||y_{1}||_{2}^{2}}=\lambda_{\text{max}}(PQ),$
(59)

here $y_{1}$ is the eigenvector associated to the largest eigenvalue of $PQ$.
Thus, the $\lambda_{\text{max}}(\mathcal{PQ})\leq\lambda_{\text{max}}(PQ)$ is
proved.  $\square$