Model reduction for LPV systems based on approximate
modal decomposition
T. Luspay, T. Péni, I. Gőzse, Z. Szabó, B. Vanek
Systems and Control Laboratory

arXiv:1609.06948v1 [cs.SY] 22 Sep 2016

of Institute for Computer Science and Control,
1111 Budapest, Kende u. 13-17., Hungary
E-mail: peni.tamas@sztaki.mta.hu

Abstract
The paper presents a novel model order reduction technique for large-scale linear parameter varying
(LPV) systems. The approach is based on decoupling the original dynamics into smaller dimensional
LPV subsystems that can be independently reduced by parameter varying reduction methods. The
decomposition starts with the construction of a modal transformation that separates the modal subsystems. Hierarchical clustering is applied then to collect the dynamically similar modal subsystems
into larger groups. The subsystems formed from the groups are then independently reduced. This
approach substantially differs from most of the previously proposed LPV model reduction techniques,
since it performs manipulations on the LPV model and not on a set of linear time-invariant (LTI)
models defined at fixed scheduling parameter values. Therefore the model interpolation, which is the
most challenging part of most reduction techniques, is avoided. The applicability of the developed
algorithm is thoroughly investigated and demonstrated by numerical case studies.
Keywords: linear parameter-varying systems, model order reduction, balanced realization,
modal transformation, clustering

1

Introduction

The LPV models have proven to be useful in system analysis and control design, because they are
able to represent a wide class of nonlinear systems while preserving the advantageous properties
of the linear time invariant (LTI) structure [26], [7], [27]. Using LPV models the analysis and
control synthesis tasks are generally casted to convex optimization problems involving linear matrix
inequality (LMI) constraints [42]. As long as the complexity of the system is low (the dimension of
the state is smaller than 20-30, the number of scheduling variables is at most 2 or 3) these problems
can be efficiently solved by off-the-shelf semidefinite solvers. On the other hand, the modeling of
complex systems (e.g. flexible structures [24]) results in high dimensional LPV systems (even with
100-1000 states), the dimension of which has to be reduced in order to make them numerically
tractable. It is important to emphasize that our aim is to reduce the dimension of the state vector
and not of the scheduling variable. The latter problem is fundamentally different and is addressed
e.g. in papers [16], [30].

1

Model order reduction for LTI systems is a well-studied topic, see e.g. reference [3]. The same
problem for LPV systems was first addressed in [40], [41], where the concept of balanced realization
based model truncation [23] was extended to parameter varying systems. The approach is based
on the balancing state transformation computed from the parameter-varying controllability and
observability Gramians. The computation of the Gramians involves LMI optimization, which suffers
from the same computational limitations as the LPV analysis and synthesis problems. Therefore,
this model reduction technique can only be applied to systems of moderate complexity. In order
to avoid these limitations an approximate balanced truncation method is proposed in [35]. The
approach is based on reformulating the balanced truncation as a Petrov-Galerkin (oblique) projection
[3] that involves two parameter-dependent transformation matrices. These transformations depend
again on the controllability and observability Gramians, but instead of computing them by LMI
optimization, [35] proposes to directly approximate the transformation matrices by using the local
Gramians associated with the LTI systems obtained at frozen scheduling parameter values. Since
the algorithm uses only QR factorization and singular value decomposition (SVD) it is numerically
attractive. On the other hand, it works only for stable systems and the approximation rises some
technical questions. Finding answers to these questions is part of the ongoing research.
In the literature some other model reduction methods can also be found for LPV systems, see
e.g. [28], [36], [1]. Although these techniques are different, but share a common feature; they are
based on frozen parameter models, i.e. they independently reduce the LTI models obtained at fixed
scheduling parameter values and then seek a suitable interpolation algorithm to construct an LPV
system from the reduced LTI model set. The latter step is very difficult in general [9], because the
independently reduced (transformed and projected) local, LTI systems have to be transformed into
a consistent state-space representation.
Regarding parameter-dependent systems it is important to mention the family of parametric
model reduction methods, see e.g. [6], [12], [25], [2]. Though these approaches show similar characteristics to the LPV model reduction the problem they address is fundamentally different. The
parametric model reduction starts from a parameterized set of large-scale LTI systems. The systems
are reduced and a model database is constructed from the obtained reduced order models. If a parameter value is given, the corresponding reduced order model is constructed from the stored model
set by a suitably chosen interpolation algorithm. It is important to emphasize that the parameters
in this framework are considered to be constant in time. This is significantly different from the
problem studied in this paper, since here the parameters change in time and thus the system to be
reduced is time-varying.
The model reduction method proposed in this paper returns to the LPV reduction technique
presented in [40], but instead of constructing a numerically tractable approximation for the balanced truncation, it decouples the large-scale system into smaller dimensional LPV subsystems that
can be independently reduced by the original algorithms of [40]. The decomposition starts with the
construction of the parameter-varying modal transformation that separates the parameter varying
modes of the system. Hierarchical clustering is applied then to group the modes into larger LPV
subsystems, which can then be independently reduced. Since the model reduction is applied on
LPV subsystems and not on frozen LTI models, the difficulties of system interpolation are avoided.
Furthermore, the reduction of unstable systems can be naturally integrated in the proposed methodology.
Although the (approximate) modal decomposition is part of our algorithm too, it is important
to emphasize that the concept is fundamentally different from that of the modal truncation methods
elaborated for LTI systems [3], [38]. Modal truncation aims at reducing the large-scale model
by keeping only the dominant modal subsystems. This is not efficient in general, because the pole
location does not necessarily indicate the contribution of the corresponding state to the input-output

2

behavior of the overall system [3]. This is a reason why we propose to cluster the modal subsystems
into larger groups and use a model reduction algorithm to reduce the larger subsystems obtained.
The paper is organized as follows. In the following section the model reduction problem is
formulated. Section 3 is devoted to the construction of the parameter-varying modal transformation.
In Section 4 the hierarchical clustering and the balanced reduction methods are discussed. The
numerical case studies are presented in Section 5. At the end of the paper the main results are
summarized and the most important conclusions are drawn.

2

Problem formulation

A continuous-time LPV system can be given in state-space form as follows
G(ρ) :

ẋ(t) = A(ρ(t))x(t) + B(ρ(t))u(t)
y(t) = C(ρ(t))x(t) + D(ρ(t))u(t),

(1)

where x : R+ → Rnx , u : R+ → Rnu and y : R+ → Rny , are respectively the state, input
and output. The matrix functions A : Rnρ → Rnx ×nx , B : Rnρ → Rnx ×nu , C : Rnρ → Rny ×nx ,
D : Rnρ → Rny ×nu are assumed to be continuous functions of the scheduling parameter ρ : R+ → R.
We moreover assume that both ρ and ρ̇ are bounded: ρ(t) ∈ Ω := [ρmin , ρmax ] and |ρ̇(t)| ≤ δ for
all t1 . (By this definition we restrict ourselves to LPV models having only 1 scheduling parameter.
The case of vector valued ρ will be commented briefly in the conclusion.)
The LPV system (1) is often given by a pair (G, Γ), where G is a set of LTI models obtained by
evaluating (1) at the finite set of scheduling parameter values stored in Γ. Formally,
Γ = {ρ1 = ρmin , ρ2 , . . . , ρN = ρmax }, ρ1 < ρ2 < . . . < ρN and
i
o
n
h
Ak = A(ρk ), Bk = B(ρk ),
k Bk
, k = 1 . . . N.
G = Gk Gk = A
Ck Dk , Ck = C(ρk ), Dk = D(ρk )

(2)

Between two consecutive grid points linear interpolation is assumed. This grid-based representation
is typical in practice, because the LPV model is often generated by trimming a nonlinear system at
different operating points [21].
 nx , such that the
The aim of the LPV model reduction is to find Gred (ρ, ρ̇) of order nred
x
difference between the input-output behavior of the reduced and the full order models is as small as
possible. Beside time-domain simulations, several metrics can also be used to measure the similarity
of the two models, see e.g. Section 5. These metrics give quantitative information on the applicability
of the reduced order model for model-based control design.

3

Approximate modal decomposition for parameter-varying
systems

For LTI systems the concept of modal decomposition is theoretically sound. Let the LTI system
be given by its state-space matrices A, B, C, D and let A be diagonalizable. Then there exists
a similarity state transformation T̄ such that the matrix Ā = T̄ −1 AT̄ is block diagonal and the
1 For

parameter-dependent variables we use the following notational convention: a(ρ) denotes the variable a as a
function of the parameter ρ, while a(ρ(t)) or a(ρk ) denotes the value of the variable at the specific parameter value
ρ(t) or ρk .

3

transformed system can be written in the

Ā1 0

 0 Ā2
ẋ = 

..

.
y=



C̄1

C̄2

form



B̄1

..
 B̄2

.
x + 
 ..

..
 .
.
0 
B̄m
0 Ām

. . . C̄m x + Du.




 u,


(3)

where each block of the structured matrix Ā corresponds to one dynamical mode of the system,
defined by the eigenvalues of A as follows: let λ1 , . . . , λm be an arbitrary sequence of the eigenvalues
such that every eigenvalue is repeated as many times as its multiplicity, but the complex conjugate
pairs are represented by only one member of the pair. Then
 

R (λi ) I (λi )

if I(λi ) 6= 0,
−I (λi ) R (λi )
Āi =
(4)

[λi ]
if I(λi ) = 0
where R(λi ) and I(λi ) denote the real and imaginay part of the complex eigenvalue λi . The similarity
transformation T̄ can be constructed from the eigenvectors of A as follows: let v1 , . . . , vm be the
set of distinct eigenvectors corresponding to the eigenvalue sequence above. Again, the complex
conjugate eigenvector pairs are represented by only one member of the pair. Then
(
R(vi )
if I(vi ) = 0
T̄ = [ν1 , . . . , νm ], where νi =
(5)
[R(vi ) I(vi )] otherwise
Note that the modal transformation is not unique, because it can be constructed from any eigenvectorset spanning the same eigen-space.
Our aim is to extend this concept to parameter-varying systems and construct a parameterdependent state transformation T̄ (ρ) that is differentiable in ρ and decouples (at least approximately)
the LPV system into modal subsystems. For this, we start from the grid-based representation (2)
and compute first the eigen-decomposition of the Ak matrices. Then an algorithm is constructed for
shaping the eigenvectors at each grid point such that the generated modal transformation sequence
T̄1 , . . . T̄N becomes suitable for smooth interpolation in ρ.

3.1

Assumptions

In order to proceed, we make the following assumptions on the full order LPV model:
(A) A(ρ) is diagonalizable and has a differentiable eigenvalue decomposition, i.e. there exist a
diagonal Λ(ρ) and invertible V (ρ) matrices such that both are differentiable in ρ and A(ρ) =
V (ρ)Λ(ρ)V (ρ)−1 . The diagonal entries λi (ρ) of Λ(ρ) are the parameter-dependent eigenvalues,
while the columns vi (ρ) of V (ρ) give the parameter-dependent eigenvectors.
(B) The multiplicity of each λi (ρ) is constant over the parameter domain.
(C) The Ak matrices do not have eigenvalues on the imaginary axis.
Assumption (A) implies first that each Ak matrix is diagonalizable, so the eigen-decomposition Λk =
Vk−1 Ak Vk exists in every grid point. Second, the analyticity of V (ρ) makes it possible to construct
4

a differentiable and invertible parameter-varying state transformation from the eigenvectors Vk .
Matrices satisfying assumption (A) have been subject of intensive research for decades, see e.g. [17],
[18], [29]. These papers provide rigorous mathematical conditions for the existence of an analytic
eigenvalue/eigenvector system. However, these conditions are rather theoretical and as explained
next, they can be relaxed in the most engineering applications. Assumption (B) is rather technical,
it guarantees that the dimension of the eigenspace associated with each λi (ρk ) is constant for all
ρk ∈ Γ. This makes it easier to transform the eigenvectors into a smooth sequence, from which a
differentiable state transformation can be constructed. Assumption (C) is necessary, because the
hyperbolic metric introduced in Section 3.2 is not defined for eigenvalues on the imaginary axis.
Note that assumption (C) does not exclude the case, when an eigenvector crosses the imaginary axis
at some ρ̄ value. If the parameter grid is chosen such that ρ̄ ∈
/ Γ, the proposed algorithms remain
applicable.
Though the assumptions above seem restrictive, they can be easily satisfied in the most real
applications. This is because the LPV formalism is only an approximation of the underlying physical
system, hence it is possible to apply minor numerical modifications in the systems description, as
long as the overall input/output behavior remains almost uneffected. Therefore, small perturbations
on the A(ρ) matrix is allowed as long as they do not significantly influence the input/output map
realized by the dynamical system. In real applications this is generally sufficient to satisfy the
assumptions above. For further details on these possible numerical corrections see Remark 1 and
the numerical example in Section 5.1.

3.2

Eigen-decomposition

If assumption (A) holds, then each Ak in (2) is diagonalizable, hence the following eigen-decomposition
can be computed for each k:
Ak → Vk−1 Ak Vk =: Λk .

(6)

Here Λk = blockdiag(λk,1 , . . . , λk,nx ) collects the eigenvalues and Vk = [vk,1 , . . . , vk,nx ] stores the
normalized eigenvectors of Ak . Note that the ordering of the eigenvalues may vary over the grid
points, so the sequence λ1,i , . . . , λN,i may not correspond to λi (ρ1 ), . . . , λi (ρN ), where λi (ρ) denotes
the i-th parameter-dependent eigenvalue of A(ρ). In order to ensure the consistency of the Λk
matrices, the ordering of the eigenvalues has to be modified in each grid point. The right ordering
can be found if the eigenvalues at every grid point are correctly paired with the eigenvalues at the
succeeding grid point. To find this pairing, two ingredients are needed: first, a distance metric
to compare the eigenvalues and second, an algorithm to find the pairing. Moreover, the algorithm
should work efficiently for large number of eigenvalues as well. In the next section a suitably distance
metric is proposed, which compares the eigenvalues based on dynamic similarity. Then the pairing
problem is reformulated as finding a perfect matching in a bipartite graph. This is beneficial, because
the latter problem can be solved efficiently in polynomial time.
Hyperbolic distance metric.
Our goal is to connect the series of local eigenvalues in such way that the resulting continuous trajectories correspond to the parameter-varying eigenvalues of the A(ρ) matrix. In case of an LTI system
each eigenvalue represents a dynamical mode, i.e. a subsystem with particular dynamical characteristics. Extending this concept to the parameter varying case, we intend to pair two eigenvalues if
they belong to the same parameter varying subsystem. For this purpose a metric is adopted, which
measures dynamic similarity between the eigenvalues. Our approach is based on the fundamental
5

results of system theory, which relate the location of the dominant poles to the transient response.
This property is more pronounced in a discrete time setup, therefore all the stable eigenvalues are
mapped into the unit circle for further discussion. Each z 0 ∈ C discrete-time eigenvalue defines a
first-order SISO transfer function G0 (z) with a single pole located at z 0 :
G0 (z) =

a0 (z)
,
z0 − z

(7)

where a0 (z) is an arbitrary polynomial in z. The collection of all possible transfer functions, generated
by a single eigenvalue z 0 is denoted by G[z 0 ]. It is clear, that G[z 0 ] ⊂ H 2 (D), where H 2 (D) denotes
the Hardy space of complex-valued functions2 [37]. This is an infinite dimensional vector space.
Among its possible orthonormal bases the discrete Laguerre system has received distinct attention
in system and control theory [33]. The Laguerre system is generated by any complex number z 00 as
follows:
p

m
1 − |z 00 |2
z − z 00
00
(m = 0, 1, 2 . . . ),
(8)
Φm (z) =
1 − (z 00 )∗ z 1 − (z 00 )∗ z
Using the Laguerre system as a basis, each transfer function in G[z 0 ] can be represented through the
linear combination of a Φm (z), i.e.
G0 (z) =

∞
X

lm Φ00m (z)

(9)

m=0

where lm are the corresponding linear coefficients. It follows from the results in [15] and [32] that
under the given conditions, the ratio of two consecutive coefficients in (9) is always constant and
depends only on the generator elements z 00 and the z 0 . That is: lm+1 /lm = r for all m = {1, 2, . . .}
and r is called the convergence factor, because it characterizes the convergence of the Laguerre series
(9). Small coefficient implies that only a few elements are dominant in (9), i.e. few basis functions are
enough to capture the dynamic behaviour of G0 (z). The metric r is thus characterizes the similarity
between the transfer functions G0 (z) generated by z 0 and the basis elements parameterized by z 00 .
Therefore, it can be considered as a possible measure of the dynamic similarity between z 0 and z 00 .
Furthermore, it can be proved (see e.g. [15]) that r can be computed by the following simple formula
r = h(z 0 , z 00 ) =

z 0 − z 00
,
1 − z 0∗ z 00

(10)

which is also known as the pseudo-hyperbolic metric between z 0 and z 00 [5], [34].
Being defined on the unit disk, this metric can only be applied on stable, discrete eigenvalues. In order to compare continuous eigenvalues, they have to be ”discretized” first by the formula
λd = exp(λc Ts ), where Ts is a sufficiently small sampling time and λc , λd denote the continuous
eigenvalue and its discrete counterpart, respectively. Formula (10) can be applied to discrete, unstable eigenvalues as well provided that they are transformed into the unit disk by the mapping
f (z) = 1/z ∗ . (Note that f (z) reflects the unstable eigenvalue across the unit circle, which guarantees that the distance between a stable and a (transformed) unstable eigenvalue is small if they
are close to each other in the sense that they can be considered as two consecutive points of the
same eigenvalue trajectory. This property is important when mixed stability eigenvalues have to
2 which

f ∈

are holomorphic on the unit disk and have finite 2-norm, i.e. kf k22 :=

H 2 (D)

6

1
2π

R 2π
0

|f (reiθ )|2 dθ < ∞ for all

be identified.) For notational convenience, the hyperbolic distance of two stable/unstable continuous/discrete eigenvalues will also be denoted by h(·, ·) and the transformations necessary to use (10)
are assumed to be performed beforehand.
Note also that (10) is 1 if either z 0 or z 00 is on the unit disk. Therefore the metric has no use
for continuous eigenvalues on the imaginary axis. This is the reason why assumption (C) has been
introduced in Section 3.1.
Returning to our original problem, the hyperbolic distance provides a possible method to compare the eigenvalues λk,i and λk+1,j based on dynamic similarity. In order to take the associated
eigenvectors vk,i and vk+1,j into consideration as well, the distance metric h(λk,i , λk+1,j ) can be
∗
∗
weighted by 1 − |vk,i
vk+1,j |, where |vk,i
vk+1,j | is the modal assurance criterion (MAC) [10] that
is a suitable measure for the directional similarity of the normalized eigenvectors. Denoting the
weighted distance metric by hw (·, ·) the final formula we use to compare two eigenvalues (and their
eigenvectors) can be given as follows:
∗
hw (λk,i , λk+1,j ) := h(λk,i , λk+1,j ) · (1 − |vk,i
vk+1,j |)

(11)

Perfect matching in complete bipartite graphs.
Let Lk and Lk+1 denote the ordered sets containing the eigenvalues at grid points k and k + 1,
respectively, i.e. Lk = {λk,1 , . . . , λk,nx }, Lk+1 = {λk+1,1 , . . . , λk+1,nx }. In order to find the pairing
between these values, a graph-theoretic reformulation is given here. For this, let the weighted,
complete, bipartite graph [8] Bk = (Nk , Ek , Ck ) be defined with vertices Nk = Lk ∪ Lk+1 , edges
Ek = {eij | eij := (λk,i , λk+1,j ), ∀(i, j) pairs} and edge-costs Ck = {cij | cij = hw (λk,i , λk+1,j )}.
So Bk is defined such that its edges connect every element in Lk with all elements in Lk+1 and
the cost of an edge characterizes the dynamical similarity between the two eigenvalues on the edge.
Note that, every λk,i has exactly one pair in Lk+1 and every λk+1,j is a pair of exactly one vertex
in Lk . In terms of the graph-theoretic setup, this is a pairing between Lk and Lk+1 , which is a
perfect matching ([19], [8]) in Bk , characterized by a set of nx independent edges in Ek . The cost
of a perfect matching is defined by the sum of the costs for the corresponding set of independent
edges. Accordingly, finding the right pairing between the eigenvalues of Lk and Lk+1 is formulated
as finding the minimum cost perfect matching in Bk [19]. The obtained matching problem can be
efficiently solved in polynomial time by using the Hungarian Method (Kuhn-Munkres Algorithm)
[8], offering a numerically attractive solution for the eigenvalue pairing problem. Considering the
ordering of the eigenvalues at ρ1 as a reference, the pairing problem can be solved successively for
k = 1, 2, . . . , N − 1. As a result the Lk sets (and consequently the Λk matrices) will be consistently
ordered. This ordering has to be applied to the columns of Vk as well, in order that the eigenvectors
become consistent with the reordered eigenvalues.

3.3

Continuity of the eigenvectors and the Procrustes-problem

The next step is to follow the idea of modal decomposition and construct (by using (5)) a parameter
varying state transformation from the eigenvectors Vk , k = 1 . . . N . For this, a sequence of local
modal transformations T̄1 , . . . T̄N is generated from the eigenvectors and then these matrices are
interpolated along the scheduling parameter. In order to facilitate the interpolation, the eigenvectors
have to be shaped first such that their entries form a smooth function along the parameter grid. Of
course, the shaping has to preserve the eigenspace spanned by the eigenvectors.
To this end, the first step is to identify the multiplicity of the eigenvalues in order that the
eigenvectors associated with a same eigenvalue can be grouped and handled together. For this, let

7

τi denote the i-th eigenvalue trajectory, i.e. τi = (λ1,i , . . . , λN,i ). By using (10), we can introduce a
distance metric between two eigenvalue trajectories as follows3 :


H(τi , τj ) = min max h(λk,i , λk,j ) , max h(λk,i , λ∗k,j ) .
(12)
k

k

This metric does not distinguish the complex pairs, i.e. the distance of τi from τj and τj∗ are the same.
This is important to render the complex conjugate eigenvalue sequences together. By computing
the distance between every (τi , τj ) pair, the multiple eigenvalue trajectories can be identified. Two
eigenvalue trajectories are considered to belong to a single repeated eigenvalue, if the distance
between them is smaller than a given threshold. (This threshold-based decision is necessary here,
because after the numerical manipulations (e.g. the eigen-decomposition) performed on the LPV
model we might not expect that the eigenvalue sequences corresponding to a multiple eigenvalue of
A(ρ) will be perfectly equal in each grid point and give 0 distance. On the other hand, if the original
system contains distinct eigenvalues that are very close to each other, then from numerical point of
view, it is generally better to handle them as a single, repeated eigenvalue.)
By assigning to every eigenvalue the eigenvectors spanning its eigenspace, an ordered set {(λk,1 , Vk,1 ), . . . , (λk,n , Vk,n )}, n ≤ nx can be defined at each grid point. Here Vk,i ∈ Cnx ×di collects
the eigenvectors associated with λk,i . Due to assumption (B) the dimension di of the eigen-space is
constant and the same at every grid point, for any i. This property will be exploited in the algorithm
below.
The next step is to transform the eigenvector sequence V1,i , . . . , VN,i for all i. This can be done by
a right multiplication of Vk,i with an invertible matrix Qk,i , which changes the eigenvectors, however
leaves the eigenspace intact. To obtain the required smooth interpolation, Qk,i should transform the
respective eigenvectors at consecutive grid points as close as possible. This condition is formulated
as a complex, unconstrained Procrustes problem [13], [12], [2] as follows:
Q̄k+1,i := arg min kVk,i − Vk+1,i Qk+1,i kF , ∀k, i
Qk+1,i

(13)

where k goes from 1 to N − 1, i ∈ {1, . . . , n}, Qk+1,i ∈ Cdi ×di , with di denoting the dimension of
the eigenspace associated with the i-th eigenvalue. The solution for (13) can be analytically given
in the following closed form:


†


R (Vk+1,i ) −I (Vk+1,i )
R (Vk,i )
vec(Q̄k+1,i ) = I ⊗
vec
(14)
I (Vk+1,i ) R (Vk+1,i )
I (Vk,i )
where vec(Q̄k+1,i ) is a vector formed formed from Q̄k+1,i by stacking its columns below each other.
Accordingly, V̄k+1,i = Vk+1,i Q̄k+1,i is the appropriately rotated eigenvector, consistent with Vk,i .
The described Procrustes problem hence can be solved successively for k = 1, 2, . . . N , for each
eigenvalue of the system. It should be noted, that the Procrustes iteration can be started from any
other grid point as well. The iteration has to be performed then in both directions: for the larger
and for the smaller indices. (For choosing the starting grid point it is a reasonable approach to
analyze the numerical conditioning of the eigenvector matrices and choose a grid point, where the
associated Vk matrix is well-conditioned. This improves the numerical stability and reliability of the
further computations.)
Remark 1 If assumption (B) does not hold the Procrustes algorithm cannot be applied because the
dimension of the eigenspaces changes at some grid points. If assumption (A) is not satisfied the
3 Note

that, (12) can also be defined by the weighted hw (·, ·) metric as well.

8

Procrustes algorithm fails to generate smooth eigenvector sequence. However, if either problem is
detected it is possible in general to correct the eigenvalues and eigenvectors at the specific grid points
such that the input-output behavior of the overall system does not change significantly, but after
the corrections, the Procrustes algorithm can run and returns a smooth solution. To illustrate the
idea, consider the case when a repeated real eigenvector trajectory switches at certain grid points to a
complex pair. In real applications this is a common source of discontinuity, which is caused in general
by the numerical inaccuracies of the eigen-decomposition. In order to resolve this issue, the complex
eigenvalues can be substituted by their real part and the corresponding eigenvector pairs can be simply
replaced by the eigenvectors associated with the closest real eigenvalue. This, of course, modifies the
A(ρ) matrix, but as long as this modification has small impact on the dynamical behavior, it is an
allowed correction to eliminate the discontinuity. The case study in Section 5.1 gives an example for
the application of this idea in practice.

3.4

Approximate modal transformation

The sequence V̄1 , . . . V̄N of the shaped eigenvector matrices (with V̄k = [V̄k,1 . . . V̄k,n ]) can be used to
construct the local modal transformations T̄1 , . . . , T̄N . Due to the corrections above, the entries of T̄k
can be smoothly interpolated (e.g. by splines) over the parameter domain. Therefore, a parameter
dependent, differentiable transformation T̄ (ρ) can be obtained. Defining a new state vector x̄ such
that T̄ (ρ)x̄ = x, the original LPV system (1) transforms into


∂ T̄ (ρ)
ρ̇ x̄ + T̄ −1 (ρ)B(ρ)u
x̄˙ = T̄ −1 (ρ)A(ρ)T̄ (ρ) − T̄ −1 (ρ)
∂ρ
(15)
= (Ā(ρ) + Ē(ρ, ρ̇))x̄ + B̄(ρ)u
y = C(ρ)T̄ (ρ)x̄ + D(ρ)u = C̄(ρ)x̄ + D(ρ)u.
where Ā(ρ) = T̄ −1 (ρ)A(ρ)T̄ (ρ) similarly to (3) is block diagonal with 1 or 2 dimensional blocks.
Ē(ρ, ρ̇) contains the ρ̇-dependent term in (15) and the difference T̄ −1 (ρ)A(ρ)T̄ (ρ) − Ā(ρ). The latter
is zero only at the grid points. The term Ē(ρ, ρ̇) represents the main challenge in the application of
the transformation, as noticed in the LPV literature. The grid-based representation of (15) can be
given as (Ḡ, Γ × Ω), where
h
i
o
n
Ēk,s B̄k
Āk = Ā(ρk ), B̄k = B̄(ρk ),
1...N
,
, k=
(16)
Ḡ = Ḡk,s Ḡk,s = Āk +
ρk ∈ Γ
C̄ = C̄(ρ ), D = D(ρ )
C̄
D
k

k

k

k

k

k

and
Ēk,s = Ē(ρk , νs ) = −T̄ −1 (ρk )

∂ T̄ (ρ)
∂ρ

νs ,

with

s ∈ {1, 2}, νs ∈ {−δ, δ}

(17)

ρ=ρk

Note that the difference term T̄ −1 (ρ)A(ρ)T̄ (ρ) − Ā(ρ) is automatically neglected, because it is 0
for every ρk ∈ Γ. It is clear, that without the Ē(ρ, ρ̇) (or Ēk,s ) term, the transformed system is
fully decoupled and is similar to the modal form (3). At this point the transformed model has to
be analyzed in order to decide if the ρ̇-dependent term Ēk,s is kept or neglected. If the original
and the modal system (i.e. (16)-(17) without Ēk,s ) do not differ significantly in terms of inputoutput behavior, the Ēk,s term can be dropped out and the computations can be proceeded with
the decoupled structure. Otherwise, it has to be kept and taken into account in the final, model
reduction step discussed in Section 4.2.

9

4

Clustering and balanced reduction

The algorithm we have developed so far is able to decouple (at least approximately) the LPV system
into a set of independent parameter-varying modal subsystems. In order the reduce the dimension
of the overall model, the idea is to group the modal subsystems of similar dynamical properties into
clusters, in order to construct larger dimensional subsystems that can be efficiently reduced. The
use of dynamic similarity is important, because we can expect significant dimension reduction only
if the subsystems are built up from similar dynamical components.
In the forthcoming two subsections, we show first how the hierarchical agglomerative clustering
(HAC) methodology [20], [14] can be adapted to solve our specific clustering problem. Then, based
on [40], the main steps of the parameter-varying balanced reduction algorithm are recalled.

4.1

Hierarchical clustering

In this section the hierarchical agglomerative clustering methodology is adopted to group the similar
modal subsystems together. The HAC framework is a bottom-up clustering approach, where each
data object is treated as a singleton cluster and successively merged until a single cluster is obtained.
In the context of the model reduction, the eigenvalue trajectories τi , i = 1, . . . , nx are considered
as the individual data objects. For measuring similarities between these trajectories, formula (12) is
used. As this metric ensures merging of complex pairs into one cluster, the parameter-varying modes
– (τi , τi∗ ) pairs – take place at the lowest level of the HAC. For the comparison of two clusters the
complete link clustering is used, i.e.: the similarity of two clusters is determined by the similarity
of their most dissimilar members. Formally, if Cm and Cn are two clusters, then the corresponding
merging criterion is:


L(Cm , Cn ) = max(H(τi , τj )), τi ∈ Cm , τj ∈ Cn .
(18)
i,j

Consequently, in the HAC framework, at each algorithmic step those two clusters are merged together
for which the (18) value is the smallest. Note that (18) is non-local, i.e. the entire structure of the
clustering can influence merge decisions. This results in a preference for compact clusters with small
diameters (i.e. the most similar dynamics are grouped together) and causes sensitivity to outliers
(i.e. uncommon dynamical components). The merging is repeated until all the objects have been
grouped into a single cluster. This can be done in O(n2x log nx ) steps [20]. The result of the HAC
is visualized by a dendrogram, which is a tree diagram illustrating how the data objects are merged
into larger clusters until the one single cluster is reached. The final cluster structure is obtained by
cutting the dendrogram at a user-defined level of similarity. The careful choice of this threshold is
important, because it determines the number and size of the clusters generated. In the decomposition
of dynamical systems the number of clusters is mainly limited by the size of the largest cluster, since
the cluster size gives the dimension of the underlying dynamical system, which cannot be arbitrarily
large due to the numerical limitations of the balanced reduction algorithm to be applied in the next
step.
Assume M clusters have been generated. By rearranging the state x̄ of (15) according to these
clusters the system structure presented in Fig. 1 is obtained. If the new state vector is denoted by
x̃ and Π is the permutation matrix mapping x̃ to x̄, i.e. x̄ = Πx̃, then
Ã(ρ) = ΠT Ā(ρ)Π,
B̃(ρ) = ΠT B̄(ρ)
T
Ẽ(ρ, ρ̇) = Π Ē(ρ, ρ̇)Π
C̃(ρ) = C̄(ρ)Π,
D̃(ρ) = D(ρ)
10

(19)

Figure 1:
The structure of the system matrices after clustering:
Ã(ρ)
=
blockdiag(Ã(1) (ρ), . . . , Ã(M ) (ρ)), Ẽ(ρ, ρ̇) is decomposed such that Ẽ(ρ, ρ̇) = Ẽ1 (ρ, ρ̇) + Ẽ2 (ρ, ρ̇) such
that the structure of Ẽ1 (ρ, ρ̇) is aligned with the structure of Ã(ρ). The white areas denote the zero
entries.
or by using the grid-based representation it can be given by the pair (G̃, Γ × Ω), where
h
n
i
o
. . . N, s = 1, 2
k = Ã(ρk ), B̃k = B̃(ρk ),
G̃ = G̃k,s G̃k,s = Ãk + Ẽk,s B̃k , Ã
, kρk=∈1Γ,
νs ∈ {−δ, δ}
C̃k = C̃(ρk ), Dk = D(ρk )
C̃k
Dk

(20)

and Ẽk,s =
Ẽ(ρ, ρ̇) the dynamics are fully decoupled into M subsystems
k , νs ). Without
h Ẽ(ρ
i
Ã(`) (ρ) B̃ (`) (ρ)
(`)
G̃ (ρ) =
, ` = 1, . . . , M . The dimension of each subsystem equals to the size
(`)
(`)
C̃

(ρ) D

(ρ)

of the corresponding cluster. The cluster size is limited by the numerical complexity of the balanced
reduction algorithm, which will be performed on each subsystems next. Taking this numerical limitation into account, the maximum cluster size should not be larger than 30-40. For more details,
see the discussion in the next subsection. If we have decided earlier to neglect the coupling term
Ẽ(ρ, ρ̇), then these subsystems can be handled separately. Otherwise, the coupling term has to be
taken into consideration. Since a decoupled structure is needed to continue our algorithm, this is
only partially possible. For this, let Ẽ(ρ, ρ̇) be expressed as a sum of Ẽ1 (ρ, ρ̇) and Ẽ2 (ρ, ρ̇) such
that the structure of Ẽ1 (ρ, ρ̇) is aligned with the structure of Ã (see Fig. 1). Assuming that the
effect of Ẽ2 (ρ, ρ̇) on the input-output behavior is negligible, we can proceed the computations with
(`)
the subsystems defined by the system matrices Ã(`) (ρ) + Ẽ1 (ρ, ρ̇), B̃ (`) (ρ), C̃ (`) (ρ) and D(`) (ρ).
These systems are also independent, so they can be separately reduced. The neglected term can be
treated as a modeling uncertainty during the analysis or control synthesis procedures performed on
the reduced order model. (To further minimize the approximation error it is a reasonable idea to
complete the distance metric (12) with an additional term penalizing large entries in Ẽ2 (ρ, ρ̇).)

11

4.2

Balanced reduction

Balanced reduction is a fundamental approach for the model reduction of linear (time invariant
and varying, as well as parameter-dependent) systems [41], [40]. The key concept is the balanced
realization which reveals the controllability and observability properties of the system. The similarity
transformation T̂ (ρ), which transforms a general LPV system (1) into balanced form is obtained
from the observability Xo (ρ) and controllability Xc (ρ) Gramians4 . If the LPV system is given in
a state-space
is a priori fixed (e.g. in the form Xo (ρ) =
Pnb form and the structure of the Gramians
Pnb
Xo,0 + i=1
Xo,i fi (ρ)Xo,i and Xc (ρ) = Xc,0 + i=1
Xc,i gi (ρ)Xc,i where fi (ρ), gi (ρ) are fixed basis
functions and Xo,i , Xc,i , i = 1, . . . nb are free variables), then Xo (ρ) and Xc (ρ) can be obtained as
a result of the following optimization problem [40]:
X
min
trace Xo (ρk )Xc (ρk )
Xo,i ,Xc,i ,i=1...nb

dXo (ρ)
dρ

ρ=ρk

dXc (ρ)
−
dρ

ρ=ρk

k

νs + A(ρk )T Xo (ρk ) + Xo (ρk )A(ρk ) + C(ρk )T C(ρk ) ≺ 0
(21)
νs + A(ρk )Xc (ρk ) + Xc (ρk )A(ρk )T + B(ρk )B(ρk )T ≺ 0

Xo (ρk )  0, Xc (ρk )  0, ∀k ∈ {1, . . . , N }, ρk ∈ Γ and s ∈ {1, 2}, νs ∈ {−δ, δ}
This is a nonconvex optimization problem, but if either Xo (ρ) or Xc (ρ) is fixed, then the cost function becomes linear in the remaining variables, hence the problem reduces to a linear optimization
problem with Linear Matrix Inequality (LMI) constraints. As suggested in the literature by alternately fixing Xo (ρ) and Xc (ρ) a numerically tractable iterative algorithm is obtained, where an
initial Xo (ρ) (or Xc (ρ)) can be calculated from the point-wise (at fixed ρk values) time-invariant
solutions.
Although the above modifications ease the computational burden, the approach still suffers from
the curse of dimensionality. The number of decision variables involved, and hence the numerical
complexity grows with the state-dimension of the system. If the number of states goes over 30-40,
the LMI optimization problem becomes intractable by the off-the-shelf semidefinite solvers. The
decomposition of the system into smaller, independent subsystems, offers a remedy for this problem.
Starting from the decoupled structure (20) the next step is to reduce the M subsystems separately
by the parameter varying balanced reduction method. After computing the controllability and
(`)
(`)
observability Gramians Xo (ρ) and Xc (ρ) for every subsystem, ` = 1, . . . M , the parameter-varying
balanced transformations are constructed in the following way. First the unique Cholesky factors
(`)
(`)
Ro (ρ) and Rc (ρ) (both are analytic in ρ) are determined by the algorithm given in the proof of
Theorem 7.8.1 in [40]:

T
Xo(`) (ρ) = Ro(`) (ρ) · Ro(`) (ρ), such that Ro(`) (ρ) is upper triangular
(22)

T
(`)
(`)
(`)
Xc (ρ) = Rc (ρ) Rc (ρ) , such that Rc (ρ) is lower triangular
(`)

(`)

Then a parameter-varying singular value decomposition is performed on the product Ro (ρ)Rc (ρ):

T
U (`) (ρ)S (`) (ρ) V (`) (ρ)
= Ro(`) (ρ)Rc(`) (ρ)
(23)
4 The presented balanced reduction algorithm can only be applied to quadratically stable LPV systems. The
extension of the method to unstable systems is well documented in [40] and [41].

12

(`)

where S (`) (ρ) is the diagonal matrix of the parameter-dependent singular values σj (ρ) and S (`) (ρ),
U (`) (ρ) and V (`) (ρ) are all analytic in ρ. The analyticity is ensured by computing (23) by the Analytic
Singular Value Decomposition (ASVD) algorithm proposed in [22]. The balancing transformation
T̂ (`) (ρ) for the `-th subsystem can be computed as

− 21
.
T̂ (`) (ρ) = Rc(`) (ρ)V (`) (ρ) S (`) (ρ)

(24)

The balancing transformations can be applied on the subsystems and the states corresponding to
small singular values can be eliminated. Since T̂ (`) (ρ) transformations are parameter-dependent, the
reduced systems explicitly depend on ρ̇ as well [40]. Having reduced the subsystems individually,
the reduced dynamics are finally joined together to obtain the low dimensional approximation of
(2).

5

Case studies

5.1

A benchmark example

As a first example a 80 dimensional, 2-input-2-output LPV system has been generated for the
numerical evaluation of the developed algorithm. The procedure used to generate the model consists
of 4 steps. First, the parameter-variation of each eigenvalue is defined over the parameter domain
Ω = [0, 1]. The λi (ρ) functions obtained are then used to construct the parameter-varying blockdiagonal A0 (ρ) ∈ R80×80 matrix. This is then completed with a randomly generated constant
input B0 ∈ R80×2 and output C0 ∈ R2×80 mappings, as well as a constant direct feedthrough term
D0 ∈ R2×2 . To make the problem more realistic and industrially relevant, the generated modal
system is transformed in the second step by a parameter-varying matrix T (ρ). T (ρ) is constructed
by randomly generating an invertible matrix and making its certain (randomly selected) blocks
parameter-dependent. The resulting LPV system can be given as follows:

ẋ = T (ρ)−1 A0 (ρ)T (ρ)x + T (ρ)−1 B0 u
G0 (ρ) :
(25)
y = C0 T (ρ)x + D0 u
The third step is evaluating the system above over an equidistant grid Γ0 containing N0 parameter
values. Then a d-th degree polynomial is fitted to each entry of the system matrices in order to
deform the eigenvalue trajectories and thus make the model more realistic. Finally, in the fourth
step, a grid-based representation is generated by evaluating the model over the grid Γ. This way, the
LPV system is free from any special structure and can be considered as a generic parameter-varying
model. To obtain the particular model used in this section the length N0 of the initial grid, the
length N of the final grid and the degree d of the interpolating polynomials have been chosen to be
60, 100 and 14, respectively. Furthermore, to evaluate and challenge the algorithm, a wide range of
dynamical behaviors have been covered including:
• parameter varying real- and complex conjugate eigenvalues,
• parameter independent dynamics,
• higher order complex and real eigenvalues (with algebraic multiplicity),
• integrators,
• mixed stability eigenvalues (that cross the imaginary axis),
13

Figure 2: Pole migration of the LPV benchmark system. The value of the scheduling parameter is
indicated by colors according to the color-bar on the right.
• complex conjugate - real transitions.
These modes can be recognized on the eigenvalue trajectories depicted in Figure 2. The numerical
testing of the algorithm can be performed according to the following steps:
• Eigen-decomposition. Using standard numerical methods, the eigenvalues and the corresponding eigenvectors of the Ak matrices are computed at each grid point.
• Identification of integrators. An initial ordering is carried out, based on the absolute value
of the eigenvalues in order to locate and label integrators for further computations. This is a
rather technical step, necessary to satisfy assumption (C). The eigenvalues corresponding to
the integrators are simply skipped from the reduction procedure, and will be added back only
at the end, when the reduced order LPV model is assembled.
• Finding the eigenvalue trajectories. The samples of the pointwise modes over the parameter domain are connected by the Hungarian Algorithm, based on their weighted hyperbolic
distance. Unstable poles are projected into the unit circle during the construction of the distance matrix, as discussed in Section 3.2. Continuous eigenvalue trajectories are restored as a
function of the scheduling parameter. Integrators are excluded from the pairing and matched
individually.
• Detection of multiple eigenvalue trajectories. Eigenvalue trajectories which are in close proximity for every value of the scheduling parameter are grouped together and handled as muliple
ones. In a same manner, using the hyperbolic metric between trajectories, irregular behaviours
(i.e. complex-real transitions) can be detected and labelled. For example, in this particular
example there are two eigenvalues very close to each other, accompanied by occasional complexreal transitions. This represents discontinuity in the eigenspace (assumptions (A) and (B) are
14

violated), hence a numerical correction (see Remark 1) is needed. The absolute variance for
the corresponding poles over the scheduling parameter domain was found to be in the range
of 10−7 clearly indicating a constant real pole with multiplicity of 2. Therefore we can replace
the values with their average, and the complex eigenvectors can be replaced by the closest real
eigenvectors. These steps succesfully resolve the mathematical discontinuity without sensible
change in the input-output behaviour.
• Eigenvector smoothing. The complex Procrustes problem is solved for the matched eigenvectors, where the subspaces of the previously labelled poles are treated accordingly. The smoothness of the modal transformations are compared numerically, by comparing the ρ̇ dependent
terms. For this purpose, a cubic spline interpolation is applied for the eigenvectors obtained
from the eigen-decomposition, and for the one obtained through the Procrustes problem. These
continuous functions are then evaluated between grid points. The results are summarized in
Table 1, where the large entries indicate discontinuity in ∂ T̄ (ρ)/∂ρ, which implies that the state
transformation T̄ (ρ) is probably not differentiable. The results illustrate the effectiveness of
the proposed complex Procrustes smoothing.

Before Procrustes ∂T∂ρ(ρ)
After Procrustes ∂T∂ρ(ρ)

maxρ max {·}

meanρ max {·}

maxρ k·k2

meanρ k·k2

100.52

73.6

312.28

246.28

4.57

1.07

12.6

2.85

Table 1: Numerical comparison of the maximal elements and the matrix norms of modal transformations before- and after smoothing, over the parameter domain. The numerical values were
obtained by applying the functions labelling the columns to the matrices labelling the rows, i.e. the
1,1 entry is obtained by computing the maximal entry of ∂ T̄ (ρ)/∂ρ as a function of the parameter
and maximizing it over ρ.
• Modal form. Applying transformation T̄ (ρ), the parameter-varying modal form is obtained.
At this point a numerical test has to be carried out to investigate the effect of the ρ̇-dependent
term. It can be checked via time domain simulations that in this particular example this term
can be neglected without significant change in the input/output response. (see also Table 1).
• Stable-unstable decomposition. By using the continuous eigenvalue trajectories and the parametervarying modal form, unstable (as well as mixed stability) dynamics can be separated and removed from the system. In the underlying example, five states can be separated (see Figure 2)
and thus the model reduction is performed for the remaining 75 (stable) states. In this example, this is a reasonable step, as the number of unstable/mixed stability eigenvalue trajectories
is small. If this was not the case, these eigenvalues had to be preserved, which would result
in unstable subsystems after clustering. In that case the variant of the balanced reduction
algorithm extended to unstable systems [40] would have to be applied.
• Clustering. Hiearchical clustering of the matched eigenvalue trajectories are performed next,
aiming to reveal dynamical redundancies of the system. Figure 3 shows the corresponding
dendrogram plot, which is used for representing the arrangement of the clusters. In Figure 3
the height of each line equals to the distance between the two data objects (either eigenvalue
trajectory or cluster of trajectories, computed by (18) ) below. Based on the dendrogram,
different number of clusters can be generated and it is the task of the user to decide. To make
this decision the following trade-off has to be taken into consideration: large number of clusters
15

corresponds to smaller sized groups, in which state elimination is generally more conservative
(consider the extreme case, where each cluster contains a single mode). On the other hand,
small number of clusters imply higher dimensional groups with increasing numerical burden
(on the other extreme: the entire system considered as a single cluster). Taking this into
account, the given system are subdivided into five clusters (see Figure 3), with the following
dimensions: 26 (red), 34 (purple), 7 (blue), 6 (green) and 2 (black). Note that the last cluster
cannot be reduced, since it contains a single complex-conjugate mode. According to Figure 3,
a dynamic dissimilarity can be observed, implying the preservation of this mode in question.

Figure 3: Dendrogram for clustering in the benchmark example. The horizontal axis is labelled by
the indices of the data objects (eigenvalue trajectories) while the vertical axis shows the similarity
distance computed by using (18). The dashed line indicates the threshold where the dendrogram
was cut to obtain the 5 clusters.
In addition, the cophenetic correlation coefficient is computed, which is often used for characterizing a dendrogram: how faithfully does it represent the similarities among data objects
[20]. The magnitude of the cophenetic correlation should be close to 1 for the case of a good
description. The computed value was 0.83, which shows that the selected hyperbolic metric is
indeed a good indicator for comparing different dynamical behaviours.
• Computing Gramians. The parameter-dependent controllability and observability Gramians
are computed for each individual cluster. By fixing the structures of the Gramians as Xc (ρ) =
Xc,0 + ρXc,1 and Xo (ρ) = Xo,0 + ρXo,1 the iterative optimization (Section 4.2) is carried
out for each sub-system by using the MOSEK optimization tool [4]. This corresponds to
(351 + 595 + 28 + 21) × 4 decision variables in four separate optimization problems, which is
a significant decrease compared to the 5550 × 4 variables involved in the single LMI problem
for the entire 75 dimensional system.

16

• Analytic Singular Value Decomposition. In order to determine the number of states with
the largest contributions to the input-output behaviour of each subsystem, the ASVDs are
computed individually according to (23). The number of most significant singular values
followed by the dimension of the corresponding cluster are 5/26, 5/34, 3/7 and 2/6.
• Balanced Transformation. The parameter-varying balancing transformations and their inverses
can be computed for each subsystem, using the singular value decomposition (23). Applying
these transformations, balanced forms are obtained, resulting in a ρ, ρ̇ dependent system. The
balanced models are then truncated according to the most significant singular values, computed
in the previous step. That is, the 75 dimensional stable part is reduced to 17 dimension.
• Reduced model construction. Finally, the individually reduced subsystems are joined together
and the five dimensional unstable dynamics are added back. Hence the final, reduced-order
LPV model has 22 states (compared to the original 80) and given in a grid-based fashion,
depending on ρ and ρ̇.
• Evaluation of the reduced model. In order to validate the reduced model, various numerical
properties have been investigated. One of the main motivation of model reduction is to obtain
a smaller dimensional representation, which is suitable for controller design. From this point
of view, similartity of closed-loop behaviours between the high-order and low-order models
is of capital importance. In order to take the feedback control objective into account, the
widely-used ν-gap metric is evaluated, defined between LTI systems, G1 (jω) and G2 (jω), as
[11], [39]:
δν (G1 (jω), G2 (jω)) =
1

1

(I + G2 (jω)G∗2 (jω))− 2 (G1 (jω) − G2 (jω))(I + G∗1 (jω)G1 (jω))− 2

∞

.

(26)

Essentially, if δν (G1 , G2 ) ≤ β, then a controller, which stabilizes G2 also stabilizes G1 , with a
stability margin of β [39]. For identical systems δν (G1 , G2 ) = 0, otherwise it is 0 ≤ δν ≤ 1.
For LPV systems (26) can be interpreted in different ways. One possibility is to compare
point-wise LTI systems of the LPV dynamics by taking:
max δν (Gk (jω), G̃k (jω))
ω

(27)

at each grid point k ∈ [1 100]. Figure 4 shows the evaluation of the above expression for
models interpolated between the original grid points. The maximal value is 0.12, implying a
satisfactory similarity between the full and reduced order models. The second representation
of (26) is illustrated on Figure 5, where the interpolated systems are compared at every ωi
frequency, i.e.:
max δν (G( ρ)(jωi ), G̃(ρ)(jωi )).
ρ

(28)

Figure 5 also provides an insight on the frequency-domain properties of the LPV model reduction.
Finally, the proposed algorithm is compared with a local reduction approach. At each grid
point a balanced transformation based model reduction is performed for the corresponding
high dimensional LTI model. The local Hankel singular values imply lower dimensional models,
within the range of 2 − 10, for the sake of consistency the dimension has been fixed for 10
at each grid-point. This fact indicates some conservativeness of the proposed methodology,
17

Figure 4: ν-gap metric of the full (80) and reduced (22) order model over the scheduling parameter
domain

Figure 5: Frequency-wise maximum ν-gap metric of the full (80) and reduced (22) order model over
the scheduling parameter domain

18

Figure 6: Pole migration map for the reduced order model (left) and for the collection of locally
reduced systems (right). The value of the scheduling parameter is indicated by colors according to
the color-bar on the right.
which, in general preserves more states as a consequence of the modal decomposition. The
pole migration maps for the two approaches are compared in Figure 6. It can be observed,
that the resulting parameter varying model has a smooth pole map, which makes the LPV
model generation straightforward and less challenging. On the other hand, the pole map for
the set of locally reduced models shows large variation (see the right plot in Figure 6). This
implies the need of a more refined interpolation technique to successfully recover time domain
behaviour of the original plant, which also makes control design virtually infeasible.

5.2

B-1 aircraft

Our second example is the publicly available nonlinear aeroelastic model of the Rockwell B-1 aircraft.
This model is widely known and investigated due to the aeroelastic effects observed at subsonic
speeds. To obtain an LPV system for model reduction the simulation model from [31] is used. By
linearizing the nonlinear model at the flight altitude of 15000f t with different Mach numbers between
0.6 − 0.75 grid-based LPV dynamics with 20 state variables are obtained. The 20 states consist of
the 10 dimensional rigid body dynamics (bank angle φ, pitch angle θ, yaw angle ψ, roll, pitch and
yaw rates of the center of gravity: p, q, r, and x − y − z axis velocities in the body coordinate frame:
U , V , W and altitude h) and the additional 10 states of the five flexible modes. Eight measured
outputs are considered: flight path angle γ, the roll, pitch and yaw rates of the center of gravity, and
the lateral and vertical accelerations of the cockpit and center of gravity. The model has 9 inputs:
throttle T , right and left (symmetric) horizontal tails: δHR , δHL , upper and lower split rudder
surfaces: δRU , δRL , wing upper-surface spoilers: δSR , δSL and canard control vanes: δCV R , δCV L .
For these actuators a 13 dimensional LTI dynamics is given and used through the augmentation of

19

Figure 7: Pitch rate response to a positive one degree deflection of the symmetric horizontal tail.
the LPV model. Accordingly, the final model used for testing the model reduction algorithm is 33
dimensional.
The reduced order model can be obtained by performing the following steps. After obtaining the
parameter-varying modal form, three mixed stability modes are removed. The resulting 30 states
can be grouped into 3 clusters with dimensions of 23, 5 and 2. Based on the parameter-varying
Gramians, the ASVD plots indicated 7, 5 and 2 significant dimensions (i.e. only the largest block
had been reduced). Hence, the resulting reduced-order model has a 14 dimensional stable and 3
dimensional unstable part.
To investigate the numerical properties of the reduced order model, time-domain simulations are
carried out first. In these simulations the flight speed is set to vary as:
V (t) = 0.675 + 0.05 sin(0.2t),

(29)

m
which corresponds to a maximal acceleration of ≈ 3.3 sec
2 . Figure 7 compares the pitch rate responses
to a positive-one-degree commanded symmetric-horizontal-tail deflection of the full order model
(with 33 states) and the reduced order one (17 states). It can be depicted that the two models
give almost identical responses. Figure 8 shows the roll rate responses, where the horizontal-tail
deflections are set asymmetrically δHL = −1◦ , δHR = 1◦ . Again, only a small discrepency can be
observed between the two model.
As noted previously, the LPV balanced transformation and reduction suffers from the curse of
dimensionality, i.e. it runs into numerical problems for higher dimensional systems. Accordingly, it
cannot be applied directly to the benchmark example in Section 5.1. At the same time, for the case
of the 33 dimensional model it is possible to solve the underlying LMI problems and consequently
obtain a reduced description. Accordingly, after a stable-unstable decomposition the parametervarying singular values have been computed for the stable part by using the parameter-varying
Gramians and the corresponding balancing transformation [40]. The singular-value plots are given
in Figure 9 implying the elimination of 17 states, which corresponds to a 16 dimensional reduced

20

Figure 8: Roll rate response to a one degree asymmetric horizontal tail deflection.
system. This shows and verifies that the proposed approach indeed give a good estimation for the
order of the reduced system.
Finally, the frequency-domain properties are compared in Figures 10 and 11. The Bode amplitude
and phase diagrams of the i) full (blue), ii) reduced with the proposed approach (red) and iii) reduced
by the method of [40] (black), models have been compared from the upper rudder surface δRU to
the roll rate r at two different flight speed (0.6M and the interpolated point 0.7312M ). Note that,
according to the applied balanced residualization, model mismatch appears in the higher frequency
domain. Still, the frequency domain approximation is satisfactory in the operation domain.

6

Conclusion

A novel model order reduction algorithm, based on approximate system decomposition, has been
developed for LPV systems. The fundamental element of the algorithm is a parameter varying state
transformation, which transforms (at least approximately) the LPV system into modal form. For
this purpose the local transformation matrices have been appropriately modified and interpolated.
The obtained structure is then subjected to a dynamic similarity analysis by using a hierarchical
clustering method. This step reveals dynamically redundant modes of the system, which are grouped
together to form parameter-varying subsystems. Finally, these smaller dimensional subsystems are
separately reduced by parameter varying balanced reduction algorithm. The effectiveness of the
methodology is illustrated by in-depth numerical studies.
The present framework requires some assumptions on the LPV system to be reduced. It is
shown that most of these assumptions can be relaxed if slight modifications are allowed on the
system matrices. The most stringent constraint is the one limiting the dimension of the scheduling
variable to 1. Though the main concept of the proposed model reduction method can be extended
to systems having more than one scheduling parameters, the numerical details of the algorithm
have to be elaborated. The most critical point here is the separation of the modal subsystems, i.e.
the generalization of the Hungarian algorithm and the iterative Procrustes smoothing. But once
21

Figure 9: Singular Values for the stable part of the B1 aircraft model. The value of the scheduling
parameter is indicated by colors according to the color-bar on the right.

Figure 10: Bode plots from δRU to r of the full model (33 states) at 0.6 Mach compared to the two
reduced order models. The 17-state model has been obtained by the proposed algorithm, while the
balanced reduction algorithm from [40] has generated the 16-state model.

22

Figure 11: Bode plots from δRU to r of the full model (33 states) at 0.7312 Mach compared to the
two reduced order models. The 17-state model has been obtained by the proposed algorithm, while
the balanced reduction algorithm from [40] has generated the 16-state model.
this is done the procedure can be continued since the clustering and the balanced reduction do not
depend on the dimension of the scheduling parameter. The extension of the modal decomposition
to multi-parameter LPV systems is one important direction for the future research.
The well-known numerical problems related to eigen-decomposition is certainly one area which
also needs further considerations, especially for very large dimensional and ill-conditioned problems. Furthermore, currently we are seeking more rigoruous methods for the decomposition of LPV
systems. Guaranteed error bounds of the method will also be part of our future research.

Acknowledgement
The research leading to these results is part of the FLEXOP project. This project has received
funding from the European Unions Horizon 2020 research and innovation programme under grant
agreement No 636307. The authors would like to thank Peter Seiler (University of Minnesota)
and Arnar Hjartarson (formerly MUSYN Inc.) for providing insight related to the B-1 simulator
model, originally developed by David Schmidt. The authors also would like to thank prof. Gary
Balas (University of Minnesota) for initiating the research on control oriented modeling of flexible
aircrafts, which gave motivation for working on large-scale LPV systems.

References
[1] F. Adegas, I. Sonderby, M. Hansen, and J. Stoustrup. Reduced-order LPV model of flexible
wind turbines from high fidelity aeroelastic codes. In IEEE Conference on Control Applications,
pages 424–429, 2013.

23

[2] D. Amsallem and Ch. Farhat. An online method for interpolating linear parametric reducedorder models. SIAM Journal on Scientific Computing, 33(5):2169–2198, 2011.
[3] A. Antoulas. Approximation of Large-Scale Dynamical Systems. Advances in Design and Control. SIAM, 2005.
[4] MOSEK ApS. The MOSEK optimization toolbox for MATLAB manual, Version 7.1 (Revision
28) edition, 2015.
[5] A.F. Beardon and D. Minda. The hyperbolic metric and geometric function theory. In Proceedings of the International Workshop on Quasiconformal Mappings and their Applications,
2000.
[6] P. Benner, S. Grundel, and N. Hornung. Parametric model reduction with small H2 -error using
radial basis functions. Advances in Computational Mathematics, 41:1231–1253, 2015.
[7] J. Bokor and G. Balas. Detection filter design for LPV systems–a geometric approach. Automatica, 40:511–518, 2004.
[8] R. Burkard, M. Dell’Amico, and S. Martello. Assignment Problems. SIAM, 2009.
[9] J. de Caigny, R. Pintelon, J. F. Camino, and J. Swevers. Interpolated Modeling of LPV Systems.
IEEE Transactions on Control Systems Technology, 22(6):2232–2246, 2014.
[10] D. J. Ewins. Modal Testing: Theory and Practice. John Wiley & Sons, 1984.
[11] T. T. Georgiou. On the computation of the gap metric. Systems Control Letters, 11:253–257,
1988.
[12] M. Geuss, H. Panzer, and B. Lohmann. On parametric model order reduction by matrix
interpolation. In European Control Conference, pages 3433–3438, 2013.
[13] G. H. Golub and Ch. F. Van Loan. Matrix computations. Johns Hopkins University Press,
2013.
[14] T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning. Data Mining,
Inference, and Prediction. Springer, 2009.
[15] P. S. C. Heuberger and P. M. J. Van den Hof. Modelling and Identification with Rational
Orthogonal Basis Functions. Springer, Dordrecht, 2006.
[16] A. Kwiatkowski and H. Werner. PCA-Based Parameter Set Mappings for LPV Models With
Fewer Parameters and Less Overbounding. IEEE Transactions on Control Systems Technology,
16(4):781–788, 2008.
[17] P. Lancaster. On eigenvalues of matrices dependent on a parameter. Numerische Mathematik,
6:377–387, 1964.
[18] P. Lancaster, A. S. Markus, and F. Zhou. Perturbation Theory for Analytic Matrix Functions:
The Semisimple Case. SIAM Journal on Matrix Analysis and Applications, 25(3):606–626, 2003.
[19] L. Lovasz and M. D. Plummer. Matching Theory, volume 29 of Annals of Discrete Mathematics.
North-Holland Mathematics Studies, 1986.

24

[20] Ch. D. Manning and P. Raghavan abd H. Schütze. Introduction to Information Retrieval.
Cambridge University Press, 2008.
[21] A. Marcos and G. Balas. Development of linear-parameter-varying models for aircraft. Journal
of Guidance, Control and Dynamics, 27(2):218–228, 2004.
[22] V. Mehrmann and W. Rath. Numerical methods for the computation of Analyic Singular Value
Decompositions. Electronic Transactions on Numerical Analysis, 1:72–88, 1993.
[23] B. C. Moore. Principle Component Aanalysis in Linear Systems: Controllability, Observability,
and Model Reduction. IEEE Transactions on Automatic Control, 26(1):17–32, 1981.
[24] C. P. Moreno, P. J. Seiler, and G. J. Balas. Model Reduction for Aeroservoelastic Systems.
Journal of Aircraft, 51(1):280–290, 2014.
[25] H. Panzer, J. Mohring, R. Eid, and B. Lohmann. Parametric Model Reduction by Matrix
Interpolation. Automatisierungstechnik, 58:475–484, 2010.
[26] T. Peni and P. Seiler. Computation of lower bounds for the induced L2 -norm of LPV systems.
International Journal of Robust and Nonlinear Control, 26(4):646–661, 2016.
[27] H. Pfifer and P. Seiler. Robustness analysis of linear parameter varying systems using integral
quadratic constraints. International Journal of Robust and Nonlinear Control, 25(15):2843–
2864, 2015.
[28] C. Poussot-Vassal and C. Roos. Generation of a reduced-order LPV/LFT model from a set of
large-scale MIMO LTI flexible aircraft models. Control Engineering Practice, 20:919–930, 2012.
[29] J. Quian, D. Chu, and R. C. E. Tan. Analyticity of Semisimple Eigenvalues and Corresponding
Eigenvectors of Matrix-Valued Functions. SIAM Journal on Matrix Analysis and Applications,
36(4):1542–1566, 2015.
[30] S. Z. Rizvi, J. Mohammadpour, R. Toth, and N. Meskin. A Kernel-Based PCA Approach to
Model Reduction of Linear Parameter-Varying Systems. IEEE Transactions on Control Systems
Technology, 24(5):1883–1891, 2016.
[31] D. Schmidt. A nonlinear simulink simulation of a large flexible aircraft. Technical report, 2013.
[32] A. Soumelidis, J. Bokor, and F. Schipp. Realizing system poles identification on the unit disc
based on Laguerre representations and hyperbolic metrics. In Mediterranean Conference on
Control and Automation (MED), pages 1208–1213, 2013.
[33] A. Soumelidis, F. Schipp, and J. Bokor. Pole Structure Estimation from Laguerre Representations Using Hyperbolic Metrics on the Unit Disc. In Conference on Decision and Control
(CDC) and European Control Conference, pages 2136–2141, 2011.
[34] Z. Szabó and József Bokor. Non-Euclidean Geometries in Modeling and Control. Széchenyi
University Press, 2015.
[35] J. Theis, P. Seiler, and H. Werner. Model Order Reduction by Parameter-Varying Oblique
Projection. American Control Conference, 2016.
[36] J. Theis, B. Takarics, H. Pfifer, G. J. Balas, and H. Werner. Modal matching for LPV model
reduction of aeroservoelastic vehicles. In AIAA Science and Technology Forum, 2015.
25

[37] R. Toth. Modeling and Identification of Linear Parameter-Varying Systems. Springer-Verlag,
2010.
[38] A. Varga. Enhanced Modal Approach for Model Reduction. Mathematical Modelling of Systems,
1:91–105, 1995.
[39] G. Vinnicombe. Measuring Robustness of Feedback Systems. PhD thesis, Department of Engineering, University of Cambridge, 1993.
[40] G. D. Wood. Control of Parameter-Dependent Mechanical Systems. PhD thesis, University of
Cambridge, 1995.
[41] G. D. Wood, P. J. Goddard, and K. Glover. Approximation of linear parameter-varying systems.
In Conference on Decision and Control (CDC), pages 406–411, 1996.
[42] Fen Wu. Control of Linear Parameter Varying Systems. PhD thesis, University of California
at Berkeley, 1995.

26

