1. BT-based - We first compute the reduced-order models using balanced truncation and then transform the resulting system into an eigenvector basis such that ˆA is diagonal. Then, we choose the mirror images with respect to the vertical lineLα :={s∈C:(s) =α}, of the poles of this reduced system as starting interpolation points. The columns of matrix Bˆ∗ and rows of matrix ˆC∗,transformed accordingly to the matrix ˆA, are selected as right and left tangential directions, respectively.
2. Random-Complex (R-C) - The interpolation data are all complex and generated randomly. The requirement for interpolation points is that they should be elements of C+α and should not lie too far away from the vertical line Lα.
3. Random-Real (R-R) - Interpolation points are randomly chosen real points lying in C+α, tangential directions are purely real and also gene-rated randomly with values varying between 0 and 1.
4. Uniform-Real-Dir1 (U-R-D1) - The left and right tangential directions are taken to be the same as in the set above (R-R), while the interpolation points are uniformly distributed real numbers varying between 2α+ 0.1 and 2α+ 100.
5. Uniform-Real-Dir2 (U-R-D2) - The interpolation points stay the same as in (U-R-D1) but now we generate a second (different) random set of real tangential directions.
The first set of interpolation data is chosen in order to investigate the possible H2,α-improvement of results obtained by balanced truncation. From our exten-sive earlier experiments we learned that eigenvalues of reduced systems received by balanced truncation or satisfying the first order necessaryH2,α-norm opti-mality conditions are often very good approximation of eigenvalues of the original system lying closest to the vertical line Lα. Therefore the complex points (R-C) are chosen randomly in C+α with the constraint that their real part is not too much bigger than α. In order to examine the dependence of the method on the choice of starting points we fix, for sets R-R and U-R-D1, the same starting set of interpolation directions but different starting points.
Analogously, for U-R-D1 and U-R-D2, we fix the starting set of points and let the directions differ from each other so that the dependence of the technique on the set of starting directions could be investigated.
By using the MIRIAm algorithm, we aimed to construct the reduced-order systems satisfying, for α varying from −0.024 to 1.6 and different starting data, the first order necessary H2,α-optimality conditions. In most cases, these systems turn out to be asymptotically stable, even forα > 0. Thus, we were able to measure the quality of approximation by computing the relative
4.5. NUMERICAL RESULTS 81
2 4 6 8 10 12 14 16 18 20
10−5 10−4 10−3 10−2 10−1
Dimension of the reduced order model
Relative H2−norm of the error system
Relative H
2−norm of the error vs. dimension of the reduced order model
α=−0.0240 α=−0.0160 α=−0.0080 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 α=1.1 α=1.2 α=1.3 α=1.4 α=1.5 α=1.6
(a) Set R-C (H2-norm)
2 4 6 8 10 12 14 16 18 20
10−4 10−3 10−2 10−1
Dimension of the reduced order model
Relative H2,α−−norm of the error system
Relative H
2,α−norm of the error vs. dimension of the reduced order model
α=−0.0240 α=−0.0160 α=−0.0080 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 α=1.1 α=1.2 α=1.3 α=1.4 α=1.5 α=1.6
(b) Set R-C (H2,α-norm)
2 4 6 8 10 12 14 16 18 20
10−5 10−4 10−3 10−2 10−1
Dimension of the reduced order model
Relative H2−norm of the error system
Relative H
2−norm of the error vs. dimension of the reduced order model
α=−0.0240 α=−0.0160 α=−0.0080 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 α=1.1 α=1.2 α=1.3 α=1.4 α=1.5 α=1.6
(c) Set R-R (H2-norm)
2 4 6 8 10 12 14 16 18 20
10−5 10−4 10−3 10−2
Dimension of the reduced order model
Relative H2,α−norm of the error system
Relative H2,α−norm of the error vs. dimension of the reduced order model
α=−0.0240 α=−0.0160 α=−0.0080 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 α=1.1 α=1.2 α=1.3 α=1.4 α=1.5 α=1.6
(d) Set R-R (H2,α-norm)
2 4 6 8 10 12 14 16 18 20
10−5 10−4 10−3 10−2 10−1
Dimension of the reduced order model
Relative H2−norm of the error system
Relative H
2−norm of the error vs. dimension of the reduced order model
α=−0.0240 α=−0.0160 α=−0.0080 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 α=1.1 α=1.2 α=1.3 α=1.4 α=1.5 α=1.6
(e) Set BT-based (H2-norm)
2 4 6 8 10 12 14 16 18 20
10−5 10−4 10−3 10−2 10−1
Dimension of the reduced order model
Relative H2,α−norm of the error system
Relative H
2,α−norm of the error vs. dimension of the reduced order model
α=−0.0240 α=−0.0160 α=−0.0080 α=0 α=0.1 α=0.2 α=0.3 α=0.4 α=0.5 α=0.6 α=0.7 α=0.8 α=0.9 α=1 α=1.1 α=1.2 α=1.3 α=1.4 α=1.5 α=1.6
(f) Set BT-based (H2,α-norm)
Figure 4.4: The relative H2- and H2,α-norm of the error system between original system and its reduced-order approximant obtained by MIRIAm with different starting vectors vs. the dimension of the reduced-order model for α varying from −0.024 to 1.6.
8 10 12 14 16 18 20 0
1 2 3 4 5 6 7 8 9x 10−5
Dimension of the reduced order model Relative H2−norm of the error system
Relative H
2−norm of the error vs. dimension of the reduced order model
BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(a)α=−0.024 (H2-norm)
8 10 12 14 16 18 20
0 1 2 3 4 5 6 7 8 9x 10−5
Dimension of the reduced order model Relative H2,α−norm of the error system
Relative H2,α−norm of the error vs. dimension of the reduced order model
BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(b)α=−0.024 (H2,α-norm)
8 10 12 14 16 18 20
0 1 2 3 4 5 6 7 8 9x 10−5
Dimension of the reduced order model Relative H2−norm of the error system
Relative H
2−norm of the error vs. dimension of the reduced order model BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(c)α= 0 (H2-norm andH2,α-norm)
8 10 12 14 16 18 20
0 1 2 3 4 5 6 7 8 9x 10−5
The dimension of the reduced order model Relative H2−norm of the error system
Relative H
2−norm of the error vs. dimension of the reduced order model BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(d)α= 0.1 (H2-norm)
8 10 12 14 16 18 20
0 0.2 0.4 0.6 0.8 1 1.2x 10−4
Dimension of the reduced order model Relative H2,α−norm of the error system
Relative H2,α−norm of the error vs. dimension of the reduced order model
BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(e)α= 0.1 (H2,α-norm)
8 10 12 14 16 18 20
0 0.2 0.4 0.6 0.8 1 1.2x 10−4
Dimension of the reduced order model Relative H2−norm of the error system
Relative H2−norm of the error vs. dimension of the reduced order model
BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(f)α= 1.0 (H2-norm)
8 10 12 14 16 18 20
0 1 x 10−4
Dimension of the reduced order model Relative H2,α−norm of the error system
Relative H2,α−norm of the error vs. dimension of the reduced order model
BT−based R−C R−R U−R−D1 U−R−D2 BT SPA
(g)α= 1.0 (H2,α-norm)
Figure 4.5: The relative H2- and H2,α-norm of the error system between the original system and its reduced-order approximant obtained by MIRIAm, ba-lanced truncation and singular perturbation approximation vs. the dimension of the reduced-order model for different starting data
4.5. NUMERICAL RESULTS 83 H2- and H2,α-norms of the error between original system and reduced-order one, satisfying the first order necessary H2,α-optimality conditions. We per-formed several experiments with even and odd dimensions of the reduced-order model and we observed that in the latter case the number of iterations needed for MIRIAm algorithm to converge was very large. It could be due to the fact that the considered systems were real and therefore their eigenvalues occurred in complex conjugate pairs. After performing the reduction we obtained a reduced-order model of odd dimension, such that all eigenvalues except one appeared in complex conjugate pairs, and the real one was jumping between several different values until it finally converged to a certain number.
Model order reduction of a large-scale real system to a lower but odd dimension can be not very beneficial. For instance, we have reduced the order of a system using balanced truncation. This method is known to be a good approximation technique in bothH2- andH∞-norms, but in its standard form it is expensive.
We observed that the reduction to an odd dimension was of course possible but the comparison of the norms showed that for real systems, in most cases, the improvement of the approximation was very small compared to preceding even dimension. We conclude that, it is better to take a slightly larger but even number as the dimension of the reduced-order model. Therefore, in this thesis we consider only reduced models of even dimension. In what follows, we discuss the results obtained from these experiments.
In Figure 4.4 we present a comparison of the relative H2- and H2,α-norms of the error systems. As starting data we choose: Random-Complex data (Figure 4.4a and 4.4b), Random-Real data (Figure 4.4c and 4.4d) and BT-ba-sed data (Figure 4.4e and 4.4f). We plotted these results versus even dimen-sions of the reduced-order model varying from 2 to 20. The figure illustrates that for small variations of α, deviations of the relative H2,α- andH2-norm of the errors are also small. The difference is hardly distinguishable for reduced-order models of sufficiently large dimension. However, we could also observe that the starting points as well as tangential directions influenced the quality of approximation. Two upcoming questions are how to decide when the dimen-sion is sufficiently large and what does it mean that α do not differ too much.
In case of gramian-based model reduction methods the dimension could be deduced from information about Hankel singular values. However, in case of interpolation-based methods both aforementioned problems are still not solved and further research is needed. In our extensive experiments we observed that the correct choice ofαstrongly depends on the considered system. For the CD player model, the eigenvalues with the largest and with the smallest real part are λmax = −0.024344 + 2.4343i and λmin = −800.9−40037i, respectively.
In this case, we chose α to vary from −0.024 to 1.6. For almost all starting data the resulting H2,α-optimal reduced systems were both α-bounded and asymptotically stable. It is not the case, when dynamical systems with only
very small eigenvalues are considered. Then one needs to be very careful with the selection ofα.
In all experiments performed with asymptotically stable systems, mostly with medium-size examples coming from the Niconet benchmark collection [1], we also constructed reduced-order models using balanced truncation and singular perturbation approximation. For these systems we computed theH2-norm as well as theH2,α-norm. In Figure 4.5, we present a comparison of these results with those obtained by MIRIAm method with several different sets of starting data. These experiments were done for even dimensions of the reduced-order model varying from 2 to 20 and for several different values of parameter α.
Forα ∈ {−0.024,0,0.1,1} the results are plotted in Figure 4.5. As we already explained, the quality of MIRIAm method strongly depends on the choice of starting data. In general, we get results very close to the ones obtained by balanced truncation, although with inappropriate selection of starting data we run into local minima, which do not give a good approximation. If sufficiently large dimensions of the reduced-order model were chosen, we always obtained a very good approximation of the original input-output behavior. The best results, usually after very few iterations, were obtained when ”BT-based” in-terpolation data were taken. In the H2-norm sense, taking random, real or complex, choices of data usually lead to results close to the ones obtained by balanced truncation.
Although the reduced-order models constructed by means of the MIRIAm method are usually very good approximation of the original large-scale model, more sophisticated numerical methods for constructing the projection matrices as well as for computing the eigenvalues and eigenvectors have to be imple-mented. Further research has to be done in order to prove and improve the convergence of the MIRIAm method as well as to ensure the preservation of α-boundedness of the original model in its reduced-order approximant.
The dependence of the quality of approximation on the choice of α and the starting data has to be investigated in more details. Some intuitive selections for starting points have been explained earlier, but the choice of tangential directions is still an open problem. It is important since they both influence the convergence of MIRIAm algorithm as well as the quality of approximation.
Chapter 5
h 2 ,α -optimal model reduction
This chapter is devoted to the study of a model reduction technique, which aims to minimize the h2,α-norm of the error between α-bounded large-scale model and its α-bounded reduced-order approximant. The results presented here are a counterpart of the results established in Chapter 4 for discrete-time LTI systems. To a large extend they are the same. However, for the sake of clarity of the presentation, it is necessary to write them separately.
The chapter is organized as follows. After introducing the concept of α-bounded-ness, we derive an explicit formula for the h2,α-norm of theα-bounded system.
Next, the minimization problem is formulated and new tangential interpolation-based first order necessary h2,α-optimality conditions are established. In Sub-section 5.2.3, an analogue of MIRIAm algorithm for discrete-time systems is presented. In the last section, the potential of the new method is illustrated on the shallow water model. For the special case, when α= 1, we will prove that the gramian-based Wilson conditions [24] and Hyland and Bernstein conditions [64] are equivalent to the h2,α-optimality conditions introduced in this thesis.
Then we will argue, how these gramian-based conditions can be generalized to the case ofα-bounded systems.
5.1 h
2,α-Norm
In this thesis, we consider both stable and unstable LTI MIMO discrete-time systems. However, the standard h2-norm, commonly used in control theory, is not defined for unstable ones. Therefore, we introduce the concept of α-boundedness for discrete systems.
Definition 5.1.1 (α-boundedness). A discrete-time LTI system of the form (3.1.6) is called α-bounded, if all eigenvalues λ1, . . . , λN of the state matrix A
85
lie inside a disc around the origin with radius α, i.e., for all k = 1, . . . , N, λk∈D−α, where D−α :={s∈C:|s|< α, α >0}.
For all α-bounded systems, the h2,α-norm of the transfer function exists and is given by (2.3.13). Below, we derive an explicit formula for this norm.
Theorem 5.1.1. Let H(s) be a transfer function of the α-bounded discrete-time SISO LTI system (3.1.6), whose polesλ1, . . . , λN are pairwise distinct and belong toD−α, where D−α:={s∈C:|s|< α, α >0}. Let φk= Res (H(s), s=λk) be the corresponding residues. Then the h2,α-norm of H is given by
H2h2,α =
XN k=1
φ∗k λ∗kH
α2 λ∗k
. (5.1.1)
Proof. Lets=eiθ, θ ∈[0,2π]. From equation (2.2.2), we see that H(αs) =
XN k=1
φk αs−λk =
XN k=1
φk
α
s− λαk.
This means thatλk, for all k = 1, . . . , N is a pole of H(s) if and only if λk
α is a pole of H(αs). Note that for all k = 1, . . . , N,
φk= Res (H(s), s=λk) ⇔ φk
α = Res
H(αs), s= λk α
.
From the assumption that λk ∈ D−α, we conclude that λαk ∈ D−, for all k = 1, . . . , N, where D− :={s∈C:|s|<1}. One can also show that, for s=eiθ, θ∈[0,2π],
[H(αs)]∗ =
" N X
k=1
φk αs−λk
#∗
=
XN k=1
φ∗k
α s −λ∗k.
Now, the fact that λ∗k ∈ D−α allows us to deduce that the poles of [H(αs)]∗ belong toD+,where D+ :={s∈C:|s| ≥1}. Consequently, the product [H(αs)]∗H(αs) has inD− the same poles as H(αs), namely λα1, . . . ,λαN. By substituting eiθ →s in the expression (2.3.13) of the h2,α-norm, we obtain that
H2h2,α= 1 2π
Z2π
0
H(αeiθ)∗H(αeiθ)dθ = 1 2πi
Z
γ
1
s[H(αs)]∗H(αs)ds, whereγ(t) :=eit, t∈[0,2π],is a parametrization of the unit circle.
5.1. h2,α-NORM 87 Applying Cauchy’s Residue Theorem (Theorem 4.2.1) to the above formula yields
H2h2,α=
XN k=1
Res
1
s[H(αs)]∗H(αs), s= λk α
=
XN k=1
Res
H(αs), s= λk α
lim
s→λkα
1
s[H(αs)]∗
=
XN k=1
φk λk
H
α2 λ∗k
∗
. The last equality follows from the fact that
lim
s→λkα
1
s[H(αs)]∗= lim
s→λkα
1 s
XN l=1
φ∗l
α s −λ∗l
!
= α λk
XN l=1
φ∗l
α2 λk −λ∗l
= α λk
XN
l=1
φl
α2 λ∗k −λl
∗
= α λk
H
α2 λ∗k
∗
.
Since Hh2,α is a real number, we have H2h2,α =
XN k=1
φk λk
H
α2 λ∗k
∗
=
XN k=1
φ∗k λ∗kH
α2 λ∗k
and the proof is complete.
We are now ready to generalize Theorem 5.1.1 to the MIMO case.
Theorem 5.1.2. Let H(s) be a transfer function of the α-bounded discrete-time MIMO LTI system (3.1.6). Suppose that system matrices B and C are partitioned as in (2.2.3) andAis given in the diagonal form: A=diag(λ1,. . . ,λN) where λi =λj for i=j andλk∈D−α for allk = 1, . . . , N. Then the h2,α-norm of H is given by
H2h2,α = trace
( 1 λ∗k
XN k=1
H
α2 λ∗k
b∗kc∗k
)
. (5.1.2)
Proof. The h2,α-norm, (2.3.13), of a MIMO system can be rewritten in the
following form:
H2h2,α= 1 2π
Z2π
0
trace[H(αeiθ)]∗H(αeiθ)dθ
= 1 2π
Z2π
0
Xm l=1
Xp q=1
[Hql(αeiθ)]∗Hql(αeiθ)dθ
=
Xm l=1
Xp q=1
1 2π
Z2π
0
Hql(αeiθ)∗Hql(αeiθ)dθ
=
Xm l=1
Xp q=1
Hql2h2,α.
Here,Hqldenotes the (q, l)-th component of the transfer functionH,see (2.2.4).
After inserting (5.1.1) into the above formula, we conclude that H2h2,α =
Xm l=1
Xp q=1
XN k=1
(φqlk)∗ λ∗k Hql
α2 λ∗k
,
where residuesφqlk are defined as in (2.2.6), namelyφqlk =cqkbkl. Now the proof of (5.1.2) is a matter of straightforward computation.
Corollary 5.1.1. Let Σ = (A, B, C) and Σ = ( ˆˆ A,B,ˆ C)ˆ be state-space repre-sentations of the original and reduced systems (3.1.6) and (3.1.7), respectively.
Assume that the state matricesA and Aˆ are diagonal and their poles belong to D−α. LetH andHˆ be the corresponding transfer functions. Then theh2,α-norm of the error systemΣe = Σ−Σ, denoted byˆ J, can be represented by
J =H−Hˆ2h2,α = trace
8<
: PN k=1
1 λ∗k
H
α2 λ∗k
−Hˆ
α2 λ∗k
b∗kc∗k + Pn
k=1 ˆ1 λ∗k
Hˆ
α2 λˆ∗k
−H
α2 ˆλ∗k
ˆb∗kˆc∗k
9=
;,
(5.1.3)
where λ∗k, λˆ∗k are the poles of Σ and Σ, respectively.ˆ
Proof. The system matrices of the error system Σe are given by Ae = A 0
0 Aˆ
!
, Be = B Bˆ
!
, Ce =C −Cˆ (5.1.4) with transfer functionHe=H−H. Asˆ Ae is also diagonal,λ1,. . . ,λN,λˆ1,. . . ,λˆn are the poles of Σe. The h2,α-norm of the error system is an obvious conse-quence of Theorem 5.1.2.