On the stability of projection-based linear reduced-order models: Descriptor vs non-descriptor forms

1.9 MB·English

by  David Amsallem

On the stability of projection-based linear reduced-order models: Descriptor
vs non-descriptor forms
D. Amsallema,∗, C. Farhatb
aDepartment of Aeronautics and Astronautics, Durand Building, 496 Lomita Mall, Stanford University, Stanford,
94305-4035, USA
bDepartment of Aeronautics and Astronautics, Department of Mechanical Engineering, and Institute for Computational &
Mathematical Engineering, Durand Building, 496 Lomita Mall, Stanford University, Stanford, 94305-4035, USA
2
1
0 Abstract
2
Two comprehensive approaches are considered for constructing projection-based reduced-order computa-
p
tional models for linear dynamical systems. The first one reduces the governing equations written in the
e
S descriptor form, using a Galerkin or Petrov-Galerkin projection onto a reduced-order basis or pair of them,
respectively. ThesebasescanbeconstructedbyanypreferredmethodsuchastheProperOrthogonalDecom-
5
position, Balanced Proper Orthogonal Decomposition, or Moment Matching method. The second approach
2
transforms first the governing equations into their non-descriptor form before applying the same projection-

] based model order reduction method. For several structural and coupled fluid-structure dynamical systems,
S
itisobservedthatthefirstreductionapproachleadstoreduced-ordermodelsthatexhibitsignificantlybetter
D
numerical stability and accuracy properties than those resulting from the second one. These observations
.
h are reinforced by theoretical results that anticipate and confirm the better stability properties obtained in
t general when reducing the descriptor rather than non-descriptor form of the equations governing a linear
a
m dynamical system.
[ Keywords: Balanced POD; Descriptor form; Galerkin projection; LTI dynamical system; Model reduction;

1 Moment matching; Non-descriptor form; Petrov-Galerkin projection; POD; Stability
v
4
9
6
1. Introduction
5
.
9 Linear Time-Invariant (LTI) systems are routinely used to model the dynamic behavior of electrical,
0
structural, fluid, thermal, and biological systems, to name only a few. Projection-based Model Order Re-
2
duction(MOR)methodsreducethedimensionalityofLTIsystemsbyprojectingthemontocarefullychosen
1
: Reduced Order Bases (ROBs). For many applications, MOR methods have demonstrated the ability to
v
preserve the accuracy of large-scale LTI systems [1, 2, 3, 4, 5, 6, 7, 8]. However, perserving their numerical
i
X stability properties has proved to be a much greater challenge [9, 10]. A few techniques have been success-
r fully proposed for stabilizing LTI Reduced-Order Models (ROMs) [11, 10]. However, these are not without
a limitations and do not necessarily lead to best-practices in MOR. On the other hand, this paper reports
on some observations that point to one simple and frequent source of numerical instability in MOR. It also
presents theoretical results that support these observations and lead to at least two best-practice guidelines
for performing MOR.
More specifically, this paper considers two different forms of LTI systems. The first one, known as the
descriptor form,correspondstotheoriginalformofanLTIsysteminwhichamatrixEthatisnotnecessarily
the identity pre-multiplies the state time-derivative. The second form, known as the non-descriptor form,
∗Correspondingauthor
Email addresses: [email protected](D.Amsallem),[email protected](C.Farhat)
URL: stanford.edu/~amsallem(D.Amsallem)
Preprint submitted to Elsevier September 26, 2012
results from the pre-multiplication of the descriptor form of the LTI system by E−1, when E is invertible.
Whereas both descriptor and non-descriptor forms of a given LTI system are mathematically equivalent,
their projections onto a given ROB generate two ROMs whose numerical properties are observed in this
paper to differ. Furthermore, it is also proved in this paper that Galerkin-based MOR methods applied to
the descriptor form of LTI systems deliver ROMs that are less prone to numerical instability than when
they are applied to the non-descriptor form of these systems. Therefore, two major contributions of this
paper are the advocacy of Galerkin-based projections over Petrov-Galerkin counterparts when both types of
projections are feasible for MOR, and their application to the descriptor rather than non-descriptor form of
LTI systems.
To this effect, the remainder of this paper is organized as follows. In Section 2, the linear time-invariant
equations of interest are introduced. In Section 3, the stability of LTI systems is defined for both descriptor
and non-descriptor forms. In Section 4, the reduction of both descriptor and non-descriptor forms of an LTI
system is discussed, together with the numerical stability of the resulting ROMs. In Section 5, theoretical
results associated with the MOR of structural dynamic and coupled fluid-structure systems are presented
and illustrated with numerical computations. Finally, conclusions are offered in Section 6.
2. LTI systems
In this work, first-order LTI systems of the following form are considered
dx
E (t)=Ax(t)+Bu(t)
dt (1)
y(t)=Cx(t)+Du(t),
where t denotes time, x ∈ Rn the vector of state variables, u ∈ Rp the vector of inputs and y ∈ Rq the
vector of outputs, and E ∈ Rn×n, A ∈ Rn×n, B ∈ Rn×p, C ∈ Rq×n and D ∈ Rq×p the time-invariant
matrices describing the behavior of the physical system of interest. Throughout this paper, E is assumed to
be non-singular. In Computational Fluid Dynamics (CFD), this matrix can be, for example, the diagonal
matrix of control volumes arising from the semi-discretization by the finite volume method of the governing
Euler or Navier-Stokes equations. In finite element methods for structural dynamics (in state space form)
and heat convection, it is typically the mass matrix.
Given an LTI system such as (1), the following complex-valued transfer function can be defined
H:s∈C(cid:55)−→C(sE−A)−1B+D∈Cq×p. (2)
When E (cid:54)= I , the LTI system (1) is said to be in descriptor form. Since E is assumed here to be
n
non-singular1, it can be rewritten as
dx
(t)=E−1Ax(t)+E−1Bu(t)
dt (3)
y(t)=Cx(t)+Du(t),
which is known as the non-descriptor form of the above LTI system. This form is popular, for example, in
CFD, as it allows to weight the residual in a given control volume by the inverse of this volume in order
to emphasize the regions of the computational domain associated with boundary layers where the mesh
resolution is typically finer. The transformation from the descriptor to the non-descriptor form is also often
performed to enable the application of a popular MOR method to an LTI system. In this case, a transfer
function associated with the system (3) is
H:s∈C(cid:55)−→C(s−E−1A)−1E−1B+D∈Cq×p. (4)
1ForthemoregeneralcasewhereEissingularandthereforethereisnoalternativetothedescriptorformofanLTIsystem,
modelreductionisstudiedin[12]and[13]
2
This transfer function is equal to its counterpart (2).
Throughouttheremainderofthispaper,bothsystems(1)and (3)areequippedwiththeinitialcondition
x(0)=x ∈Rn.
0
3. Stability
All High-Dimensional computational Models (HDM) of the form (1) or (3) considered in this paper are
assumed to be stable in the following sense.
Definition 1. The system (1) (respectively (3)) is said to be asympotically stable if lim x(t)=0, for
t→∞
every solution x(t) of Edx(t)=Ax(t) (cid:0)respectively dx(t)=E−1Ax(t)(cid:1).
dt dt
Remark 1. Since systems (1) and (3) are equivalent, asymptotic stability for one form is equivalent to
asymptotic stability for the other.
Next, a Lyapunov function is defined as follows.
Definition 2. V :Rn →[0,∞) is a Lyapunov function for system (1) (cid:0)respectively (3)(cid:1) if:
1. V(x)=0 if and only if x=0.
2. For every solution x(t) of Edx(t)=Ax(t) (cid:0)respectively dx(t)=E−1Ax(t)(cid:1), the following inequality
dt dt
holds
d
V(x(t))≤0. (5)
dt
Then, the following theorem holds [14].
(cid:0) (cid:1)
Theorem 1. If the system (1) respectively (3) has a Lyapunov function V, it is asymptotically stable.
Example 1. IfV(x)=xT (cid:0)ETPE(cid:1)x,wherePisasymmetricpositivedefinitematrix,theinequality(5)
can also be written in the case of the descriptor form as
ETPA+ATPE=−Q , (6)
D
where Q is a symmetric positive definite matrix.
D
Similarly, if V(x) = xTPx, where P is a symmetric positive definite matrix, the inequality (5) can also
be written in the case of the non-descriptor form as
PE−1A+ATE−TP=−Q , (7)
ND
where Q is also a symmetric positive definite matrix.
ND
4. Projection-based model order reduction
Whether applied to the descriptor (1) or non-descriptor (3) form of an LTI system, a projection-based
MOR method generates another LTI system of much smaller dimension k <<n. In general, such an MOR
method operates with two ROBs:
• A trial ROB V ∈Rn×k with full-column rank, introduced to approximate the state vector x(t) as
k
x(t)≈V x (t). (8)
k k
Theapproximatestatevectoristhenuniquelydefinedbythevectorofgeneralizedcoordinatesx ∈Rk.
k
SubstitutingtheaboveapproximationintotheoriginalLTIsystemyieldsanon-zeroresidualr(t)∈Rn.
• A test ROB W with full-column rank as well, introduced to limit the size of the residual r(t) by
k
constraining it to satisfy the orthogonality condition WTr(t)=0.
k
3
When W (cid:54)= V , the projection-based MOR method is a Petrov-Galerkin approximation method. When
k k
W =V , it becomes a Galerkin approximation method.
k k
The next two sections focus on the form of a projected LTI system originally written in descriptor or
non-descriptor form. Then, the impact of the form of the expression of the LTI system and the type of the
projection method on the numerical stability of the constructed ROM is discussed in the context of popular
methods for constructing ROBs.
4.1. Reduction of an LTI system written in descriptor form
Once a trial ROB V has been constructed, the state vector is approximated as in (8). Substituting this
k
approximation into the original system (1) gives
dx
EV k(t)=AV x (t)+Bu(t)+r(t)
k dt k k (9)
y(t)=CV x (t)+Du(t).
k k
Next, forcing the residual r(t) to be orthogonal to the test ROB W leads to the following ROM of the LTI
k
system
(cid:0)WTEV (cid:1)dxk(t)=(cid:0)WTAV (cid:1)x (t)+WTBu(t)
k k dt k k k k (10)
y(t)=CV x (t)+Du(t).
k k
ThecounterpartoftheaboveROMgeneratedbyaGalerkinprojectionisobtainedbysettingW =V ,
k k
which gives
(cid:0)VTEV (cid:1)dxk(t)=(cid:0)VTAV (cid:1)x (t)+VTBu(t)
k k dt k k k k (11)
y(t)=CV x (t)+Du(t).
k k
Ingeneral,thetrialandtestROBsarechosentobebi-orthogonalwithrespecttoawell-chosenoperator.
When operating on a descriptor form, V and W are often chosen so that WTEV =I . In this case, the
k k k k k
ROM of the LTI system can be written as
dx
k(t)=A x (t)+B u(t)
dt k k k (12)
y(t)=C x (t)+D u(t),
k k k
where
A =WTAV ∈Rk×k
k k k
B =WTB∈Rk×p
k k (13)
C =CV ∈Rq×k
k k
D =D∈Rq×p.
k
Denoting the approximated high-dimensional state by x(t) = V x (t), the following equivalent of the
(cid:98) k k
original high-dimensional system can be reconstructed
dx(cid:98)(t)=(cid:0)V WTA(cid:1)x(t)+(cid:0)V WTB(cid:1)u(t)
dt k k (cid:98) k k (14)
y(t)=Cx(t)+Du(t).
(cid:98)
This system will be particularly useful when studying the stability of the ROM.
4
4.2. Reduction of an LTI system written in non-descriptor form
For the LTI system expressed in non-descriptor form, the state approximation (8) results in
V dxk(t)=(cid:0)E−1AV (cid:1)x (t)+E−1Bu(t)+r(t)
k dt k k (15)
y(t)=CV x (t)+Du(t).
k k
Enforcing the orthogonality of the residual vector r(t) and the test ROB W yields
k
WTV dxk(t)=(cid:0)WTE−1AV (cid:1)x (t)+(cid:0)WTE−1B(cid:1)u(t)
k k dt k k k k (16)
y(t)=CV x (t)+Du(t).
k k
Here,itisnotedthatwhereasequations(15)areequivalenttotheirdescriptorformcounterparts(9),the
reduced LTI system (16) is not equivalent to its descriptor form counterpart (10). In the case of a Galerkin
projection, this system becomes
VTV dxk(t)=(cid:0)VTE−1AV (cid:1)x (t)+(cid:0)VTE−1B(cid:1)u(t)
k k dt k k k k (17)
y(t)=CV x (t)+Du(t).
k k
If the trial and test ROBs are chosen to be bi-orthogonal with respect to the identity matrix — that is
WTV =I , the resulting ROM becomes
k k k
dx
k(t)=A x (t)+B u(t)
dt k k k (18)
y(t)=C x (t)+D u(t),
k k k
where
A =WTE−1AV ∈Rk×k
k k k
B =WTE−1B∈Rk×p
k k (19)
C =CV ∈Rq×k
k k
D =D∈Rq×p.
k
As in the descriptor case, the following equivalent of the original high-dimensional system written in
non-descriptor form can be reconstructed
dx(cid:98)(t)=(cid:0)V WTE−1A(cid:1)x(t)+(cid:0)V WTE−1B(cid:1)u(t)
dt k k (cid:98) k k (20)
y(t)=Cx(t)+Du(t).
(cid:98)
4.3. Stability of the ROM
In general, a projection-based MOR method with arbitrary ROBs V and W does not preserve the
k k
stability of the LTI system to which it is applied. To illustrate this fact, consider the two LTI systems given
below
(cid:20) (cid:21) (cid:20) (cid:21) (cid:20) (cid:21)
1 0 dx 1 −3.5 dx 1 −3.5
(t)= x(t), and (t)= x(t). (21)
0 0.5 dt 0.6 −2 dt 0.6 −2
The first one is written in descriptor form, and the second one in non-descriptor form. Both systems are
dx
asymptoticallystable,buttheROM 1 =x resultingfromtheGalerkinprojectionofeitherofthemusing
dt 1
(cid:20) (cid:21)
1
the ROBs V =W = is not stable.
1 1 0
5
However, the following theoretical results show that when the HDM operators satisfy certain properties,
stability can be preserved after projection.
Theorem 2. IftheoperatorA+AT isstable—thatis,allofitseigenvalueshavenon-positiverealparts
— any ROM of the form (12) resulting from the Galerkin projection of an LTI system written in descriptor
form is asymptotically stable.
Proof. Let V(x)=xTEx, and let x denote the solution of the ROM problem (14) originating from the
(cid:98)
descriptor form of the LTI system of interest with u(t)=0 and W =V . Then,
k k
dV(x) dxT dx
(cid:98) = (cid:98) Ex+xTE (cid:98)
dt dt (cid:98) (cid:98) dt
= xT (cid:0)ATV VTE(cid:1)x+xT (cid:0)EV VTA(cid:1)x
(cid:98) k k (cid:98) (cid:98) k k (cid:98)
= xT (cid:0)ATV VTEV (cid:1)x +xT (cid:0)VTEV VTA(cid:1)x. (22)
(cid:98) k k k k k k k k (cid:98)
If VTEV =I and AT +A is stable, then
k k k
dV(x(cid:98)) =(cid:0)xTATV (cid:1)x +(cid:0)xTVTA(cid:1)x=(V x )T (AT +A)(V x )≤0. (23)
dt (cid:98) k k k k (cid:98) k k k k
Hence, V is a Lyapunov function for dx(cid:98) =(cid:0)V VTA(cid:1)x which, using Theorem 1, concludes the proof.
dt k k (cid:98)
Remark 1. For this choice of function V, there is no analog of the above theorem for the case of a
Petrov-Galerkin projection of an LTI system written in descriptor form, as the sign of
dV(x)
(cid:98) =xTVT(ATW VTE+EV WTA)V x (24)
dt k k k k k k k k
depends on the choice of the ROBs V and W . Moreover, a counterexample will be presented in Sec-
k k
tion 5.1.2 to show that a Petrov-Galerkin projection of an HDM written in descriptor form and satisfying
the properties formulated in Theorem 2 does not preserve stability.
For an LTI system expressed in the non-descriptor form, the following analog to Theorem 2 can be
established however.
Theorem 3. If the operator E−1A+ATE−T is stable, then any ROM of the form (18) resulting from
the Galerkin projection of an LTI system written in non-descriptor form is asymptotically stable.
Proof. The proof is identical to that given for the descriptor case using V(x)=xTx.
Theorem 2 and Theorem 3 highlight an advantage of the Galerkin projection method over its Petrov-
Galerkin alternative.
4.4. Sample methods for constructing a ROB
4.4.1. Proper orthogonal decomposition methods
The proper orthogonal decomposition method. TheProperOrthogonalDecomposition(POD)methodbased
on snapshots [15] computes a trial ROB V by compressing the information contained in solution snapshots
k
of the system of interest. For the case of LTI systems, these snapshots can be computed either in the time
or in the frequency domain. To simplify the notation, only the case of a single input (p = 1) and a single
output (q =1) is considered here. This case is easily generalizable to multiple inputs and outputs.
In the time domain, the snapshots are typically obtained by computing the dynamic response of an LTI
system for a given forcing input and collecting samples {x(t )}Nsnaps into a (snapshot) matrix
i i=1
(cid:2) (cid:3)
X=X(t ,··· ,t )= x(t ) ... x(t ) , (25)
1 Nsnaps 1 Nsnaps
where N denotes the number of snapshots.
snaps
In the frequency domain, each complex-valued snapshot [3] is computed as the solution of a frequency
response problem of the form
(jω E−A)x(ω )=B, (26)
i i
6
whereω denotesthesamplingcircularfrequencyofinterest. Next,thesesnapshotsarestoredinareal-valued
i
matrix as follows
(cid:2) (cid:0) (cid:1) (cid:0) (cid:1)(cid:3)
X=X(ω ,··· ,ω )= Re(x(ω )) ... Re x(ω ) Im(x(ω )) ... Im x(ω ) . (27)
1 Nsnaps 1 Nsnaps 1 Nsnaps
Whethercollectedinthetimeorfrequencydomain,thesnapshotsareindependentoftheforminwhichthe
LTIsystemiswrittenbecausebothdescriptorandnon-descriptorformsofsuchasystemaremathematically
equivalent. However, as it will be shown below, the trial ROBs V constructed using these snapshots differ.
k
In the case of an LTI system written in the descriptor form, when E is symmetric positive definite, the
POD method proceeds with computing the eigenvalue decomposition
XTEX=ΦdΛdΦdT, (28)
where XTEX is a small-size matrix. Next, the above decomposition is truncated to the matrices of first k
eigenvalues Λd and corresponding eigenvectors Φd, and Vd is constructed as follows
k k k
Vd =XΦdΛd−12 (29)
k k k
Alternatively, Vd can be constructed by first computing the Singular Value Decomposition (SVD) of the
k
matrix E12X, retaining the first k left singular vectors Yk, and defining Vkd =E−12Yk. Then, a ROM of the
form (12) can be built using a Galerkin projection based on Vd.
k
For an LTI system written in non-descriptor form, the following eigenvalue decomposition is performed
instead
XTX=ΦndΛndΦndT (30)
and Vnd is constructed as in (29)2. A ROM of the form given in (18) is then constructed using a Galerkin
k
projection based on Vnd.
k
In summary, because Φd (cid:54)= Φnd, Vd (cid:54)= Vnd and the two ROMs constructed by applying of the POD
k k k k
method to both descriptor and non-descriptor forms of the same LTI system using the same snapshots are
also different.
The balanced POD method. Information about the outputs of interest can be included in the construction
ofaROM.Forexample, theBalancedTruncation(BT)[1]methodconstructstrialandtestROBssuchthat
the ROM resulting from the Petrov-Galerkin projection is balanced in the sense that the most observable
and controllable states are retained. This method also preserves the stability of the underlying HDM.
Unfortunately for very large-scale systems, the computational cost associated with the BT method can be
prohibitive. For this reason, the Balanced Proper Orthogonal Decomposition (BPOD) method [4, 16] was
developed. This alternative method constructs at a reasonable CPU cost trial and test ROBs that lead
to an approximately balanced ROM system, but does not necessarily preserve stability (for example, see
Section 5.1.2). More specifically, BPOD relies on two sets of snapshots. The first one is known as the
set of primal snapshots. These are identical to the snapshots collected for the POD method. The second
set of snapshots is known as the set of dual snapshots. These are constructed using a dual LTI system as
summarized below.
For an LTI system written in descriptor form, the dual system is
dz
ET (t)=−ATz(t)−CTv(t)
dt (31)
w(t)=BTz(t)+DTv(t),
2ThetrialROBVnd canbeequivalentlyobtainedusingtheSVDofX
k
7
where z ∈ Rn is the dual state vector, v ∈ Rq the vector of inputs, and w ∈ Rp the vector of outputs.
Assume that the snapshots are collected in the frequency domain. In this case, each complex-valued dual
snapshot zd(ω ) associated with (31) is defined as the solution of the following frequency response problem
i
(−jω ET −AT)zd(ω )=CT. (32)
i i
As described above for the primal snapshots, the dual snapshots are also stored in a real-valued matrix
as follows
Zd =Zd(ω ,··· ,ω )=(cid:2)Re(cid:0)zd(ω )(cid:1) ... Re(cid:0)zd(ω )(cid:1) Im(cid:0)zd(ω )(cid:1) ... Im(cid:0)zd(ω )(cid:1)(cid:3).
1 Nsnaps 1 Nsnaps 1 Nsnaps
(33)
Next, the SVD of the small-size matrix ZdTEX is performed
ZdTEX=ΨdΣdΦdT, (34)
and the first k singular values Σd and associated left and right singular vectors Ψd and Φd are used to
k k k
construct
Wd =ZdΨdΣd−12, and Vd =XdΦdΣd−12. (35)
k k k k k k
Then, a ROM of the form (10) is obtained using a Petrov-Galerkin projection based on the ROBs Wd and
k
V .
k
Similarly, the dual system of an LTI system expressed in non-descriptor form is given by
dz
(t)=−ATE−Tz(t)−CTv(t)
dt (36)
w(t)=BTE−Tz(t)+DTv(t),
and the complex-valued dual snapshots znd(ω ) are defined as the solutions of the problems
i
(−jω I −ATE−T)znd(ω )=CT (37)
i n i
for multiple circular frequencies ω . Note that the comparison of Eq. (37) and Eq. (32) reveals that the dual
i
snapshots associated with both forms of a given LTI system are not identical, but related as follows
znd(ω )=ETzd(−ω ), i=1,··· ,N . (38)
i i snaps
The dual snapshots are stored in the matrix
Znd =Znd(ω ,··· ,ω )=(cid:2)Re(cid:0)znd(ω )(cid:1) ... Re(cid:0)znd(ω )(cid:1) Im(cid:0)znd(ω )(cid:1) ... Im(cid:0)znd(ω )(cid:1)(cid:3),
1 Nsnaps 1 Nsnaps 1 Nsnaps
(39)
and the following SVD is performed
ZndTX=ΨndΣndΦndT. (40)
TruncatingthefirstksingularvaluesΣnd andassociatedleftandrightsingularvectorsleadstothedefinition
k
of the matrices Ψnd and Φnd and to the construction of the test and trial ROBs as follows
k k
Wnd =ZndΨndΣnd−21, Vnd =XΦndΣnd−12. (41)
k k k k k k
Again, the Petrov-Galerkin projection based on the above test and trial ROBs leads to a ROM of the
form (16).
Focusing for simplicity on the case where the snapshots are computed in the frequency domain, the
following theorem shows that, unlike the POD method, the BPOD method constructs the same ROM
whether applied to the descriptor or non-descriptor form of the LTI system of interest.
Theorem 4. The ROMs obtained by the application of a Petrov-Galerkin projection based on trial and
test bases computed by the BPOD method to the descriptor and non-descriptor forms of an LTI system of
interest are identical if the snapshot frequency sampling {ω }Nsnaps is the same in both cases.
i i=1
Proof. See the Appendix.
8
4.4.2. Moment matching methods
The j-th moment of a transfer function H(s) at s=s is defined as
0
dH(j)
η (s )=(−1)j (s ), j =0, 1, 2, ··· (42)
j 0 ds 0
For the LTI systems considered in this paper, this moment can be expressed explicitly as
η (s )=C(s E−A)−1B+D
0 0 0
(43)
η (s )=j!C(s E−A)−(j+1)B, j ≥1.
j 0 0
It can be shown (for example, see [2, 14]) that the first k moments of an LTI system can be matched by a
projection-based MOR method provided that the trial basis V satisfies
k
range(V )=span(cid:8)(s E−A)−1B,··· ,(s E−A)−kB(cid:9). (44)
k 0 0
Since the transfer function (2) associated with the descriptor form of an LTI system is equal to its coun-
terpart (4) associated with the non-descriptor form of that system, the ranges of the trial bases associated
with both forms of the LTI system are identical. The right member of (44) describes a Krylov subspace
definedbythevector(s E−A)−1Bandthematrix(s E−A)−1. Algorithmsforconstructingthissubspace
0 0
can be found in [14, 17, 18], for example. In particular, the Arnoldi method constructs a trial ROB that
is orthogonal with respect to the identity. For an LTI system written in the descriptor form, a trial ROB
satisfyingVTEV =I cannextbeconstructedbyapost-processingstep. Then,aROMcanbebuiltusing
k k k
a Galerkin projection based on V .
k
5. Numerical stability: Case studies
5.1. Structural dynamics
5.1.1. LTI systems of interest
Linearized structural dynamics problems are typically modeled by second-order LTI systems of the form
d2z dz
M (t)+D (t)+Kz(t)=Fu(t)
dt2 dt
(45)
dz
y(t)=Gz(t)+C (t)+Hu(t),
dt
where z ∈ Rn denotes the vector of structural displacements, u and y the vectors of input and output
variables, respectively, M, D and K are the mass, damping, and stiffness matrices, respectively, and F, C,
GandHaretime-invariantmatricesdescribinginputsandoutputs. ThematricesM,DandKaretypically
Symmetric Positive Definite (SPD). Hence, throughout the remainder of this paper, they are assumed to
have this property.
Eq. (45) above can be recast in first-order form by introducing the state-vector
(cid:34)dz(cid:35)
x= dt (46)
z
containingboththedisplacementandvelocitydegreesoffreedom(dof). Thisleadstothefollowingfirst-order
LTI system written in descriptor form
(cid:20) (cid:21) (cid:20) (cid:21) (cid:20) (cid:21)
M 0 dx −D −K F
(t)= x(t)+ u(t)
0 J dt J 0 0 (47)
(cid:2) (cid:3)
y(t)= C G x(t)+Hu(t),
9
whereJ∈Rn×n isanarbitrarynon-singularmatrix. Theabovesystemcanalsobewritteninnon-descriptor
form as follows
dx (cid:20)−M−1D −M−1K(cid:21) (cid:20)M−1F(cid:21)
(t)= x(t)+ u(t)
dt I 0 0 (48)
(cid:2) (cid:3)
y(t)= C G x(t)+Hu(t).
In general, J is chosen to be the identity matrix. However, the following result suggests that a much
better choice can be made in the context of constructing a stable ROM.
Theorem 5. Ifinthedescriptorform(47)oftheLTIsystemofinterestJischosenasJ=K,thenevery
ROM constructed by the application of a Galerkin projection to the HDM (47) is asymptotically stable.
Proof. Since the matrix D is SPD,
(cid:20) (cid:21) (cid:20) (cid:21)T (cid:20) (cid:21)
−D −K −D −K −2D 0
+ = (49)
K 0 K 0 0 0
is stable. Therefore, the direct application of Theorem 2 concludes the proof of this theorem.
Remark 3. The Lyapunov function introduced in the proof of Theorem 2 corresponds to twice the sum
of the potential and kinetic energies
(cid:34)dz(cid:35)T (cid:20)M 0(cid:21)(cid:34)dz(cid:35) dzT dz
V(x)=xTEx= dt dt = M +zTKz. (50)
0 K dt dt
z z
Therefore, the decay in time of this Lyapunov function is related to the decay in time of the total energy of
the structural system due to damping.
Remark 4. Theorem 3 cannot be applied here because
(cid:20)−M−1D −M−1K(cid:21) (cid:20)−M−1D −M−1K(cid:21)T (cid:20)−M−1D−DTM−1 −M−1K(cid:21)
+ = (51)
I 0 I 0 −KM−1 0
(cid:20) (cid:21)
−2 −1
is not necessarily stable. For example if M = D = K = [1], the resulting matrix is unstable.
−1 0
Hence, Theorem 2 and this remark highlight an advantage of the descriptor form of first-order LTI systems
describing structural dynamics problems.
5.1.2. Analysis of a mass-damper-spring system
To illustrate the theoretical results presented above, a simple structural dynamics example is considered
first. It is of the academic type, and therefore offers the advantage of being easily reproducable by the
interested reader.
Thestructuraldynamicsystemassociatedwiththisexampleproblemhasaone-dimensionalchainstruc-
(cid:0) (cid:1)
ture Figure 1(a) with 20 blocks. Each of these, except the last one, contains four lumped masses, four
(cid:0)
dampers, and five springs. The last block contains one additional spring where an input u is applied see
(cid:1)
Figure 1(b) . Therefore, the corresponding one-dimensional first-order LTI system, referred to here as the
HDM, has n = 160 dof. The element properties of each block (cid:0){m }4 , {c }4 , {k }6 (cid:1), where k = 0
i i=1 i i=1 i i=1 5
except for the block where the input is applied, are given in Table 1. The input u is set to a unit step
horizontal displacement applied at t=0, and the response of the system is computed for t∈[0,100] s.
TwoMORmethodsaretestedusingthisexample: (a)aPetrov-GalerkinprojectionmethodusingBPOD
inthefrequencydomain,and(b)aGalerkinprojectionmethodusingMomentMatchingbyArnoldi’sscheme
(MMA). Each MOR method is applied to the reduction of both descriptor and non-descriptor forms of the
HDM and the stability of the resulting ROMs is assessed.
BPOD is applied here using 30 complex-valued primal and 30 complex-valued dual snapshots. These are
computed in the frequency domain, at evenly spaced frequencies in the interval ω ∈ [0,2π] rad/s. Then,
10