A generalized structure-preserving doubling algorithm for
generalized discrete-time algebraic Riccati equations
E. K.-W. Chu
∗H.-Y. Fan
†W.-W. Lin
‡Abstract
In [8], an efficient structure-preserving doubling algorithm (SDA) was proposed for the solution of descrete-time algebraic Riccati equations (DAREs). In this paper, we generalize the SDA to the G-SDA, for the generalized DARE: ETXE = ATXA + (ATXB + CST)(R + BTXB)−1(BTXA + SCT) + CQCT. Swapping of matrix products and applications of the Shermann-Morrison-Woodburg formula circumvent any explicit inversions of (possibly ill-conditioned) R and E. Selected numerical examples illustrate that the G-SDA is efficient, out-performing other algorithms.
1
Introduction
Let matrices E, A ∈ Rn×n
with E being nonsingular, Q ∈ Rp×p with Q = QT > 0 (symmetric positive
definite or s.p.d.), B ∈ Rn×m
, S ∈ Rm×p
and C ∈ Rn×p with B, C possessing full column rank, and the
s.p.d. R ∈ Rm×m. Suppose further that the matrix Q − STR−1S is symmetric positive semidefinite (s.p.s.d.).
The generalized discrete-time algebraic Riccati equation (G-DARE) has the form
ETXE = ATXA + (ATXB + CST)(R + BTXB)−1(BTXA + SCT) + CQCT. (1) Equation (1) arises frequently in discrete-time optimal control problems and optimal filter problems [10, 11, 12, 14] for a given descriptor linear system:
Exk+1= Axk+ Buk, x0= x0
yk= CTxk
(2) with the control vectors {uk} chosen through
min uk J ≡1 2 ∞ X k=0 (yTkQyk+ uTkRuk+ ykTS Tu k+ uTkSyk). (3) Let C(Q − STR−1S)CT = bC0Cb0T ≥ 0 (4) be in full rank decomposition (FRD). The systems denoted by (E, A, B) and (E, A, bC0) are assumed to
be stabilizable (S) and detectable (D), respectively. Note that (E, A, B) is stabilizable if wTB = 0 and
wTA = λEwT imply |λ| < 1 or w = 0. The system (E, A, bC0) is detectable if (ET, AT, bC0T) is stabilizable.
The optimal feedback control {u∗k} for (2) and (3) are given by
u∗k = −(R + BTX+B)−1(BTX+A + SCT)xk, (5)
∗School of Mathematical Sciences, Building 28, Monash University, VIC 3800, Australia. [email protected] †Department of Mathematics, National Tsing Hua University, Hsinchu 30043, Taiwan. [email protected] ‡Department of Mathematics, National Tsing Hua University, Hsinchu 30043, Taiwan. [email protected]
where X+ ≥ 0 is the unique s.p.s.d. solution to (1). Furthermore, the closed-loop dynamics of the system
obtained with this control
Exk+1= (A + BK)xk=A − B(R + BTX+B)−1(BTX+A + SCT) xk (6)
is asymptotically stable, i.e., lim
k→∞xk= 0 (see, e.g., [14]).
It is well-known [14] that the s.p.s.d. solution of the G-DARE (1) can be obtained via the computation of the stable deflating subspace of the matrix pencil
A − λB = A 0 B −CQCT ET −CST SCT 0 R − λ E 0 0 0 AT 0 0 −BT 0 . (7)
If the columns of [In, ETX+, Z]T span the stable deflating subspace of A − λB, then X+ ≥ 0 solves the
G-DARE (1). Here In denotes the identity matrix of order n.
It is also well-known [14] that the s.p.s.d. solution X+of G-DARE can be solved by computing the stable
deflating subspace span
In
X+E
of the reduced matrix pencil of (7):
M0− λL0= A − BR−1SCT 0 −C(Q − STR−1S)CT ET − λ E BR−1BT 0 AT − CSTR−1BT . (8)
Furthermore, it is easily seen that the pencil M0− λL0 is equivalent to the symplectic pencil
M − λL = E−1A − E−1BR−1SCT 0 −C(Q − STR−1S)CT I − λ I E−1BR−1BTE−T 0 ATE−T − CSTR−1BTE−T . (9) If the column of I X∗
span the stable deflating subspace of M − λL, then E−TX∗E−1= X+ ≥ 0 solves
the G-DARE (1). Note that a 2n × 2n matrix pencil E − λF is symplectic if and only if E J ET = F J FT, where J =
0 In
−In 0
. The matrix pencil of the form in (9) is symplectic and is said to be a standard symplectic form (SSF), a stronger symplectic property. Being a SSF is the structure we try to preserve in the numerical algorithm (see [8] for more details of the SSF).
A well-known backward stable approach based on the reordering QZ-algorithm for computing the unique s.p.s.d. solution of DAREs has been proposed by [15]. The associated code, dare, has been developed in MATLAB control toolbox [13]. Unfortunately, QZ-like algorithms do not take into account of the symplectic structure, destroying it through the iterative process. Similarly, matrix disk function/inverse free methods [2, 3, 4, 5] have been developed for solving DAREs without preserving the symplectic structure. Recently, an efficient structure-preserving doubling algorithm (SDA) [8], based on the doubling algorithm [1, 9], has been proposed for solving DAREs, while preserving the SSF at each iterative step. G-DAREs can thus be solved by applying the SDA or other algorithms to the symplectic pencil (9). However, the symplectic form in (9) requires the explicit inversion of E and R, which may be ill-conditioned.
In this paper, based on the SDA algorithm [8], we develop a generalized structure-preserving doubling algorithm (G-SDA) for solving the G-DARE (1). Inversions of ill-conditioned matrices, such as E and R, are circumvented.
2
G-SDA and QR-SWAP algorithms for G-DAREs
In this section, we develop the G-SDA for solving the G-DARE (1). Let M = E−1A − E−1BR−1SCT 0 −H + CSTR−1SCT I , L = I E−1GE−T 0 ATE−T − CSTR−1BTE−T (10)
as given in (9), where H ≡ CQCT ≥ 0 and G ≡ BR−1BT ≥ 0. The pairs (E, A, B) and (E, A, bC
0) (with bC0
given in (4)) are assumed to be (S) and (D), respectively. Let b A0 ≡ E−1(A − BR−1SCT), (11) b G0 ≡ E−1BR−1BTE−T bB0Bb0≥ 0, (12) b H0 ≡ CQCT− CSTR−1SCT bC0Cb0T ≥ 0. (13) Here “ ” denotes the operation which expands a s.p.s.d. matrix into a FRD. The SDA [8] generates the sequences { bAk, bGk, bHk}: b Ak+1 = Abk(I + bGkHbk)−1Abk, (14) b Gk+1 = Gbk+ bAkBbk(I + bBTkCbkCbkTBbk)−1BbkTAbk bBk+1Bbk+1T ≥ 0, (15) b Hk+1 = Hbk+ bATkCbk(I + bCkTBbkBbTkCbk)−1CbkTAbk bCk+1Cbk+1T ≥ 0. (16) It was proven [8], under conditions (S) and (D), that the sequences { bAk, bGk, bHk} converge respectively to
zero and the s.p.s.d. solutions of the dual and prime DAREs corresponding to the symplectic pencil (9). Consequently, {E−THbkE−1} converges to the s.p.s.d. solution of the G-DARE.
In many applications, the matrices E or R in (10) and (11)–(13) are ill-conditioned, causing numerical instability. We shall modify the SDA (14)–(16) to the G-SDA for solving G-DAREs. The process utilizes on the basic swapping process [4, Lemma 1]:
Given ET ∈ Rr×r, FT ∈ Rq×r, let Q11 Q12 Q21 Q22 ET −FT = T 0 (17)
be the QR factorization of [ET, −FT]T, where ¯F = QT21∈ Rr×q, ¯E = QT 22∈ R
q×qand T is upper
triangular. Then
E−1F = ¯F ¯E−1, (18) Here the inverses of E or ¯E are merely notational — they do not have to be constructed explicitly unless required by the particular circumstances.
For the G-SDA, we first assume, without loss of generality, that only E is ill-conditioned. We shall later describe a pre-processing step when R is ill-conditioned. Let
A0 ≡ A − BR−1SCT, G0 ≡ BR−1BT B0B0T ≥ 0, H0 ≡ CQCT− CSTR−1SCT C0C0T ≥ 0. With Gk+1≡ Gk+ AkE−1Bk(I + BTkE −TC kCkTE −1BT k) −1BT kE −1AT k Bk+1Bk+1T ,
we can rewrite the sequences { bAk, bGk, bHk} of (14)–(16) in the SDA as
b Ak+1≡ E−1Ak+1 = E−1Ak(E + BkBkTE−THk)−1Ak, (19) b Gk+1≡ E−1Gk+1E−T = E−1Bk+1BTk+1E −T ≥ 0, (20) b Hk+1≡ Hk+1 = Hk+ AkTE−TCk(I + CkTE−1BkBkTE−TCk)−1CkTE−1Ak Ck+1Ck+1T ≥ 0. (21)
Applying the basic swapping process (17)–(18) we arrive at the G-SDA which generates the sequences {E−1A
k, E−1GkE−T, Hk} without explicitly inverting E:
E−1Ak+1 = E−1Ak( ¯Ekh) T(E( ¯Eh k) T+ B kBkTH¯k)−1Ak, (22) E−1Gk+1E−T = E−1Gk+ AkB¯k(( ¯Ekb) TE¯b k+ ¯B T kCkCkTB¯k)−1B¯kTA T k E −T = E−1Bk+1BTk+1E−T ≥ 0, (23) Hk+1 = Hk+ ATkC¯k( ¯Ekc( ¯Ekc)T + ¯CkTBkBkTC¯k)−1C¯kTAk= Ck+1Ck+1T ≥ 0, (24)
where swapping produces, for all k ≥ 0,
E−THk= ¯Hk( ¯Ekh)−T, E−1Bk= ¯Bk( ¯Ekb)−1, E−TCk= ¯Ck( ¯Ekc)−T. (25)
Notice that other inverses appear in (22)–(24) and their conditioning will be analyzed in Section 3.
From [8] under conditions (S) and (D), the sequences {E−1Ak, E−1GkE−T, Hk} constructed by (22)–(25)
converge to {0, E−1G∗E−T, H∗} in with E−TH∗E−1solving the G-DARE (1). By (6) the associated optimal
control matrix Ksis obtained by
Ks = −(R + BTE−TH∗E−1B)−1(BTE−TH∗E−1A + SCT)
= −( ¯E∗E¯TR + ¯K∗B)−1K¯∗A + ¯E( ¯ETE + ¯¯ BTH∗B)¯ −1E¯TSCT
where BTE−T = ¯E−TB¯T and ( ¯BTH
∗)E−1= ¯E∗−1K¯∗ are again computed by swapping.
We now describe the basic G-SDA algorithm for solving G-DARE:
G-SDA Algorithm
Input : E, A, B, C, Q, R, S, τ (a small tolerance); Note : E is ill-conditioned, R is well-conditioned; Output : the s.p.s.d. solution X for G-DARE.
Initialize : A ← A − BR−1SCT, G ← BR−1BT B0B0T ≥ 0, H ← CQCT − CSTR−1SCT C0C0T ≥ 0, B ← B0, C ← C0; Repeat : Swap as in (17)–(18): E−TH = ¯H ¯E1−T, E−1B = ¯B ¯E2−1, E−TC = ¯C ¯E3−T; Compute W1← E ¯E1T + BBTH,¯ W2← ¯E2TE¯2+ ¯BTCCTB,¯ W3← ¯E3E¯3T + ¯CTBBTC;¯
Solve for V1, V2, V3 from W1V1= A, W2V2= ¯BT, W3V3= ¯CT,
b A ← A ¯ET 1V1, bG ← G + A ¯BV2AT, bH ← H + ATCV¯ 3A; Compute the FRDs: bG bB bBT ≥ 0, b H bC bCT ≥ 0; If k bH − Hk ≤ τ k bHk, Then X ← E−THEb −1 Stop; Else A ← bA, G ← bG, H ← bH, B ← bB, C ← bC. Go to Repeat.
When R is ill-conditioned, we swap matrices as in (17)–(18):
BR−1 = ¯R−11 B,¯ R−1(SCT) = ¯CsTR¯−T2 . (26) Then the pencil M0− λL0 in (8) is equivalent to
¯ R1A ¯RT2 − ¯R1B ¯CsT 0 − ¯R2(CQCT) ¯RT2 + ¯CsR ¯CsT R¯2ETR¯T1 − λ ¯ R1E ¯RT2 BR ¯¯ BT 0 R¯2ATR¯1T− ¯CsBTR¯T1 (27) with the left and right transformations diag{ ¯R1, ¯R2} and diag{ ¯RT2, ¯RT1}. Compute the FRDs
¯
BR ¯BT B0BT0 ≥ 0 , R¯2(CQCT) ¯R2T− ¯CsR ¯CsT C0C0T ≥ 0
and let
A0≡ ¯R1A ¯RT2 − ¯R1B ¯CsT, E0≡ ¯R1(E ¯RT2) = ¯R1E2,
the matrix pencil in (27) corresponds to a G-DARE with (E, A, B, C, Q, R, S) = (E0, A0, B0, C0, In, Ir, 0).
The G-SDA can then be applied, with R = Irbeing perfectly conditioned. Let the corresponding s.p.s.d.
so-lution be E−T0 H∗E0−1. Transforming back using (26), we see that
X+= E−TR−12 H∗R−T2 E−1= E −T
2 H∗E2−1≥ 0
solves the original G-DARE (1) corresponding to M0− λL0in (8), and the associated optimal control matrix
Ks = −(R + BTE−T2 H∗E −T 2 B) −1(BTE−T 2 H∗E −1 2 A + SC T) = −( ¯E∗E¯2TR + ¯K∗B)−1K¯∗A + ¯E2( ¯E2TE¯2+ ¯BT2H∗B¯2)−1E¯2SCT ,
with the help of the swapping BTE−T 2 = ¯E
−T
2 B¯2T and ( ¯BT2H∗)E2−T = ¯E∗−1K¯∗.
Finally, the s.p.s.d. solution of the G-DARE can be solved by computing the stable deflating subspace of the matrix pencils in (7), (8) (when R is well-conditioned) or (27) by the generalized Schur algorithm. This is equivalent to applying the command dare [13] to the G-DARE (1). Similarly, the recently developed matrix disk function methods, based on the QR-SWAP process (17)–(18) [3, 4, 5], can also be applied. For comparison we briefly describe the QR-SWAP matrix disk function method:
QR-SWAP Algorithm [4]
Input : E, A, B, C, Q, R, S, τ (a small tolerance); Output : the s.p.s.d. solution X for G-DARE;
Initialize : T ← 0n, M ← M0, L ← L0, where M0 and L0 are computed by (8);
Repeat : Compute the QR-factorization: Q11 Q12 Q21 Q22 L −M = b T 0 ;
If k bT − T k ≤ τ k bT k, Then solves the least squared problem for X∗:
M(:, 1 : n) = M(:, n + 1 : 2n)X∗;
Set X ← X∗E−1,
Else Set L ← Q22L, M ← Q21M, T ← bT ;
Go to Repeat.
End of QR-SWAP Algorithm
By (6) the associated optimal control matrix is given by
Kq = (R + BTX∗E−1B)−1BTX∗E−1A = ( ¯E∗R + ¯BT∗B)−1B¯∗TA
3
Conditioning of inversions in G-SDA
The inversion of the matricesE( ¯Eh k) T + B kBkTH¯k, ( ¯Ekb) TE¯b k+ ¯B T kCkCkTB¯k and ¯Ekc( ¯E c k) T + ¯CT kBkBkTC¯k in (22)–(24) cannot be avoided in the G-SDA. We shall show that the condition numbers of these matrices are small. Let
Mk = E Gk −Hk ET (28) where the FRDs Gk = BkBkT ≥ 0 and Hk = CkCkT ≥ 0. It is easily checked by using the
Sherman-Morrison-Woodbury formula (SMWF) that Mk−1= (E + GkE−1Hk)−1 −E−1Bk(I + BkTE −TH kE−1BkT) −1BT kE −1 E−TCk(I + CkTE −1G kE−TCk)−1CkTE −1 (ET+ H kE−1Gk)−1 . (29)
An inspection of (19)–(21) reveals that the three inverses in (22)–(24) appear in Mk−1. From the properties of norms, the maximum singular value of Mk−1 is a upper bound of the maximum singular value of any of its submatrices. Consequently, we aim to bound kMk−1k from above with a moderate quantity, thus proving that the inversions in (22)–(24) are well-conditioned.
Since the stabilizability and the detectability of (E, Ak, Bk) and (E, Ak, Ck), respectively, are preserved
for all k [8], we shall try to measure the “quality” of stabilizability and detectability and subsequently estimate kMk−1k. We shall show that the relevant partial measures of stabilizability and detectability actually improves through the iteration process. Hereafter, k · k and k · kF denote the 2-norm and Frobenius norm,
respectively. For convenience, we drop the index k in (28) and (29).
Let (E, A) be a regular pencil and let ε > 0 be a small threshold. Suppose that E is nearly singular and the “large” eigenvalues of (E, A) have only linear divisor. Then there are orthogonal equivalence transformations such that the following equivalence relationship holds:
(E, A)eq.∼ E11 E12 E21 E22 , A11 A12 0 A22 , E22, A22∈ Rr×r, (30) in which
max {¯σ(E21), ¯σ(E22)} ≤ ε min{σ(E11), σ(A22)}. (31)
Here ¯σ(·) ≡ k · k and σ(·) denote, respectively, the maximal and minimal singular values of the given matrix, and r is the number of “large” eigenvalues.
Assume that r ≤ min{m, p}, where m and p are the rank of the control matrix B and the output matrix C, respectively, as in (1), and that λ((E11, A11)) ∩ λ((E22, A22)) = ∅ and
√
rεk(E12, A12)kF
δ2 <
1 4,
where λ ((·, ·)) denotes the spectrum for the given matrix pair, and δ = dif [(E11, A11), (E22, A22)] > 0 is the
difference between the matrix pairs (E11, A11) and (E22, A22) [16, pp. 307]. Then there are Pr, P`∈ Rr×(n−r)
with max {kPrk, kP`k} ≤ 2√rε δ satisfying I 0 P` I E11 E12 E21 E22 , A11 A12 0 A22 I 0 Pr I = " e E11 E12 0 Ee22 # , " e A11 A12 0 Ae22 #! , (32) where e E11≡ E11+ E12Pr, Ee22≡ E22+ P`E12, Ae11≡ A11+ A12Pr Ae22≡ A22+ P`A12. (33)
It is easily seen that the rows of [P`, I] and the columns of
− eE11−1E12
I
, respectively, form bases of the left and right invariant subspaces of the transformed pencil corresponding to the “large” eigenvalues.
Since (E, A, B) and (E, A, C) are respectively stabilizable and detectable, we have σB≡ σ YT 2 B = σ [P`, I] B1 B2 = σ( eB2) > 0 (34) and σC ≡ σ CTX2 = σ CT 1, C T 2 − eE11−1E12 I = σ( eC2T) > 0, (35) where B ≡ B1 B2 , CT ≡ CT 1, C T 2, Be2 ≡ P`B1 + B2, eC2T ≡ −C T
1Ee11E12+ C2T, and X2 and Y2 are
respectively the unitary bases of the right- and left-eigenvectors corresponding to the large eigenvalues of (E, A). Clearly σB and σC in (34) and (35), respectively, are partial measures of stabilizability and detectability.
Without loss of generality, we let
M = E11 E12 G11 G12 E21 E22 GT12 G22 −H11 −H12 E11T E21T −HT 12 −H22 E12T E22T , (36)
where G11= B1B1T ≥ 0, G12= B1B2T, G22= B2B2T ≥ 0. Multiplying M from the left by
L`≡ diag I 0 P` I , I 0 −ET 12Ee11−T I I PT r 0 I (37) and from the right by
Lr≡ diag I 0 Pr I I − eE11−1E12 0 I , I PT ` 0 I (38)
we obtain, with eE11and eE22 given in (32),
f M = L`M Lr= e E11 0 Ge11 Ge12 0 Ee22 GeT12 Ge22 − eH11 − eH12 Ee11T 0 − eH12T − eH22 0 EeT22 , (39) where " e G11 Ge12 e GT 12 Ge22 # = B1 P`B1+ B2 BT 1, B T 1P T ` + B T 2 = B 1 e B2 h B1T, eB2Ti, " e H11 He12 e HT 12 He22 # = C 1 −ET 12Ee11−TC1+ C2 h C1T, −C1TEe11−1E12+ C2T i = C 1 e C2 h C1T, eC2Ti.
In order to estimate kM−1k we first prove the following Lemmas. Lemma 3.1. Let E, G = GT ≥ 0, H = HT
(i) If E is nonsingular, G BBT ≥ 0 and H CCT ≥ 0, then E G −H ET −1 = (E + GE−TH)−1 −E−1G(ET + HE−1G)−1 E−TH(E + GE−TH)−1 (ET + HE−1G)−1 . (40) Furthermore, we have k(E + GE−TH)−1k ≤ σ(E)−1 1 + kGkkHk σ(E)2 and
kE−1G(ET + HE−TG)−1k ≤ σ(E)−2kGk, kE−TH(E + GE−TH)−1k ≤ σ(E)−2kHk. (ii) If G = GT > 0 and H = HT > 0, then
E G −H ET −1 = H−1ET(G + EH−1ET)−1 −(H + ETG−1E)−1 (G + EH−1ET)−1 G−1E(H + ETG−1E)−1 . (41) Furthermore, we have kH−1ET(G + EH−1ET)−1k ≤ kEk σ(H)σ(G), k(H + E TG−1E)−1k ≤ 1 σ(H) and k(G + EH−1ET)−1k ≤ 1 σ(G), kG −1E(H + ETG−1E)−1k ≤ kEk σ(G)σ(H). (42) Proof. (i) Equation (40) can easily be verified. As in (29), using the SMWF twice, we have
E−1G(ET + HE−1G)−1= E−1B(I + BTE−THE−1B)−1BTE−T. Consequently, with k(I + BTE−THE−1B)−1k ≤ 1, we have
kE−1G(ET+ HE−1G)−1k ≤ σ(E)−2kGk and similar
kE−TH(E + GE−TH)−1k ≤ σ(E)−2kHk. For the diagonal blocks in
E G −H ET
−1
on the right-hand-side of (40), by the SMWF, we have (E + GE−TH)−1 = (I + E−1GE−TH)−1E−1
= [I − E−1B(I + BTE−THE−1B)−1BTE−TH]E−1. Again with k(I + BTE−THE−1B)−1k ≤ 1, we obtain
k(E + GE−TH)−1k ≤ σ(E)−1 1 + kGkkHk σ(E)2 .
(ii) Equation (41) can easily be verified. Since G = GT > 0 and H = HT > 0, we have
H−1ET(G + EH−1ET)−1= H−1ETG−1/2(I + G−1/2EH−1ETG−1/2)−1G−1/2 and (H + ETG−1E)−1 = H−1/2(I + H−1/2ETG−1EH−1/2)−1H−1/2. Consequently, kH−1ET(G + EH−1ET)−1k ≤ kEk σ(H)σ(G), k(H + E TG−1E)−1k ≤ 1 σ(H). The inequalities in (42) can also be obtained in a similar fashion.
Lemma 3.2. Let Φ ∈ Rs×t (w.l.o.g. s ≤ t) with 0 ≤ σ1≤ · · · ≤ σs≡ ¯σ(Φ) singular values. Then it holds It 0 Φ Is = 2 + σ 2 s+ σs(4 + σ2s)1/2 2 1/2 ≤ ¯σ(Φ) + 1. (43) Proof. Applying the singular value decomposition of Φ and the definition ¯σ(Z) = λmax(ZTZ)1/2(43) follows
immediately.
Theorem 3.3. Let (E, A) be the matrix pair given by (30) satisfying (31). Assume that λ((E11, A11)) ∩ λ((E22, A22)) = ∅, √ rεk(E12, A12)kF δ2 < 1 4 and Pr, P`∈ Rr×(n−r), r ≤ min{m, p} with
max{kPrk, kP`k} ≤
2√rε δ satisfying (32) and (33). From (34), (35) we let
e11≡ σ( eE11) > 0, g22≡ σ2B= σ( eG22) > 0, h22≡ σ2C= σ( eH22) > 0, (44) e g11= ¯σ( eG11), eh11= ¯σ( eH11), eg12≡ ¯σ( eG12), eh12≡ ¯σ( eH12), ee12≡ ¯σ(E12), eε ≡ ¯σ( eE22), ηg≡eg11+ e g2 12 g22 , ηh≡ eh11+ eh212 h22 , η ≡ 1 + 1 e2 11 ηgηh, κ ≡ max η2+ η 2 h e2 11 e g2 12 e2 11g222 + η2+ η 2 g e2 11 ! e h2 12 e2 11h222 !1/2 , 1 g22h22 e h2 12εe 2 h2 22 +eg 2 12eε 2 g2 22 +ge212+ eh212 !1/2 .
If κε < 1, then the matrix of the form in (36) can be estimated bye kM−1k ≤ max 2η2 e2 11 +η 2 g+ η 2 h e4 11 !1/2 , 2εe 2+ g2 22+ h222 g2 22h222 1/2 · ω 1 − κεe, (45) where ω = max 2 √ rε δ + 1 , e e12 e11 + 1 2 · max ( e g12 g22 + 1 , eh12 h22 + 1 !)2 .
Proof. Let fM = L`M Lras in (39), where L`and Lrare given by (37) and (38), respectively. By assumption
(44) that eG22and eH22are s.p.d., we eliminate eG12 ( eG12T ) and eH12( eH12T), respectively, by using G22and H22
as pivot matrices. Note that G22 and H22do not possess full rank when r > m, p. Consequently, we have
e L`M efLr= fM0+ eE where f M0 ≡ e E11 0 G˘11 0 0 Ee22 0 Ge22 − ˘H11 0 Ee11T 0 0 − eH22 0 Ee22T , e E ≡ diag (" 0 − eG12Ge−122Ee22 − eH12THe22−1Ee22 0 # , " 0 − eET22He22−1He12T − eE22TGe−122GeT12 0 #) , e L` ≡ diag I − eG12Ge−122 0 I , I − eH12He22−1 0 I , (46) e Lr ≡ diag I 0 − eH22−1He12T I , I 0 − eG−122Ge12 I , (47)
and ˘G11≡ eG11− eGT12Ge−122Ge12, ˘H11≡ eH11− eH12THe22−1He12. This implies
M−1= LrLer( fM0+ eE)−1Le`L`= LrLer(I + fM0−1E)e−1Mf0−1Le`L`. (48)
Interchanging the second and third row- (and column-) blocks in fM0, we can show the similarity relationship
f M0∼ diag (" e E11 G˘11 − ˘H11 Ee11T # , " e E22 Ge22 − eH22 Ee22T #)
From Lemma 3.1, we can apply parts (i) and (ii), respectively, to the two diagonal blocks above. One can easily check that k fM0−1Ek ≤ κe eε < 1 and obtain
k fM0−1k ≤ max 2η2 e2 11 + η 2 g e4 11 + η 2 h e4 11 !1/2 , 2εe2 g2 22h222 + 1 g2 22 + 1 h2 22 1/2 . (49)
Furthermore, applying Lemma 3.2 to (37), (38), (46) and (47) we have ¯ σ(L`), ¯σ(Lr) ≤ max 2√rε δ + 1 , e e12 e11 + 1 (50) and ¯ σ( eL`), ¯σ( eLr) ≤ max ( e g12 g22 + 1 , eh12 h22 + 1 !) . (51)
Inequality (45) then follows from (48)–(51).
Remarks: (1) In Theorem 3.3 we give a upper bound for kM−1k in terms of e11, δ and the partial measures
of stabilizability and detectability g22 and h22. Let ε in (31) be chosen so that e11and δ are not too small,
and let the systems be reasonably stabilizable and detectable so that the partial measures are reasonably large. In such circumstances, the bound in (45) is reasonably moderate for kM−1k.
(2) The spectrum of (E, A) can spread out continuously in such a way that no small ε together with large e11and δ can be found. In such circumstance, a compromised r has to be chosen so as to optimize the
bounds in (45).
(3) We can show that the partial measures are not getting worse through the iteration with increasing k. From (22)–(24), the symmetric Gk and Hk are added to by low-rank updates through the iteration,
which increases their ranks as well as their minimum eigenvalues. Recall that the partial measures are really the minimum singular values of YT
2 B and CTX2, where Y2 and X2 are unitary bases for the eigenvectors
corresponding to the large eigenvalues. For various k, they span approximately the null spaces of matrices constructed from ET and E with components related to E
21 and E22 in (30) deleted. Let these null spaces
be span respectively by the unitary Nr and N`. With the help of the properties of norms and variational
properties of eigenvalues and singular values, for various values of k, we have g22 = Y2TBkBkTY2≈ NrTBkBkTNr≥ σ2(Bk) ≥ σ2(B0),
h22 = X2TCkCkTX2≈ N`TCkCkTN`≥ σ2(Ck) ≥ σ2(C0).
Thus the partial measures for stabilizability and detectability are improving, approximately, through the iteration. The G-SDA should perform well when the original system is reasonably stabilizable and detectable, in the sense that the starting values g22 and h22 for k = 0 are bounded reasonably far away from zero, with
a small enough ε and a large enough e11.
(4) The assumption in (30) that the large eigenvalues have only linear divisors can be removed. An arbitrarily small perturbation will perturb the matrix pencil (E, A) with defective eigenvalues to one with only linear divisors. The conclusions in the Section then still hold, effectively via a continuation argument.
Alternatively, the system (E, A) can be regularized via feedback [7] before the corresponding G-DARE is solved.
(5) We require that the number of large eigenvalues r ≤ min{m, p}. This only limits size of ε and the extent to which we can benefit from circumventing the inversion of E. This will not lead to a failure of the G-SDA as E and R, although possibly ill-conditioned, are assumed to be invertible.
4
Numerical experiments for G-DAREs
For the Tables in the following examples, data for various methods are lists in columns with obvious headings. The heading “dare” is for the dare command in MATLAB [13] applied to (1), and “G-SDA” stands for our G-SDA algorithm. The heading “QR-SWAP” stands for the QR-SWAP algorithm applied to the matrix pencil in (7). There is no iteration numbers to report for dare and an ‘∗’ in the Tables indicates a failure to obtain a solution. Besides, we also report the numbers of iterations (no. ite.) for the QR-SWAP and G-SDA algorithms, respectively, and the number of stable closed-loop eigenvalues (no. stab. ev.) in the examples. The symbol |λc
max| indicates the spectral radius of the closed-loop matrix pair (E, A + BK), i.e.,
|λc
max| = max{|λ| : λ ∈ λ(E, A + BK)}.
We use trid(a, b, c) to denote the tridiagonal matrix with the main-, sub- and super-diagonal elements being a, b and c, respectively. Also, we denote
Tn= 1 −1 −1 · · · −1 0 1 −1 · · · −1 .. . . .. 1 . .. ... .. . . .. . .. −1 0 · · · 0 1 ∈ Rn×n.
For the comparison of residuals computed by these three methods, we use the ‘normalized’ residual (NRes) formula proposed in [4]
NRes ≡ kA
TXA − E˜ TXE − A˜ TXB(R + B˜ TXB)˜ −1BTXA + Hk˜ kATXAk + kE˜ TXEk + kA˜ TXB(R + B˜ TXB)˜ −1BTXAk + kHk˜ , where ˜X is an approximate solution.
All computations were performed in MATLAB/version 6.5 on a PC of Intel Pentium-III processor at 866 MHz, with RAM of 768 MB, using IEEE double-precision floating-point arithmetic (ε ≈ 2.22 × 10−16). Example 1. Consider a linear descriptor system (E, A, B, C) with n = 6 and rank(B) = rank(C) = 3. The system matrices are
E = diag(1, 10−2, 10−4, 10−6, 10−8, 10−10), R = I, A = 4.0426 3.9258 2.6310 −2.1318 5.5853 −7.1839 3.5169 −0.0108 −1.7188 −8.5395 −5.2439 −0.2965 4.1518 5.7531 2.0055 4.6018 8.2394 5.7068 1.2700 −7.3705 −5.6308 3.8215 8.0503 2.2467 1.5915 0.6336 −2.9188 5.2129 0.1337 −6.8345 4.0271 −3.9175 −2.2047 2.2661 2.8700 0.1553 , BT = −0.4820 −0.4466 −0.8810 −0.8007 0.4766 −1.2284 1.2694 0.7538 −0.8847 −1.1809 0.5286 0.3069 −0.6425 1.2407 0.1126 0.7689 −0.8265 0.2993 , CT = 0.3285 −0.9312 1.0424 1.1712 −0.0214 0.6355 0.3685 0.6990 −0.3572 −0.5304 −1.7255 −1.3765 3.0559 −2.6376 −1.2290 −1.6608 0.0370 1.3068 .
In this case, one of the closed-loop eigenvalues achieved by QR-SWAP lies outside the unit circle, with modulus equals 24.255. The numerical results are given in Table 1.
G-SDA dare QR-SWAP NRes 1.71 × 10−16 ∗ 1.47 × 10−13
no. ite. 4 - 5
no. stab. ev. 6 ∗ 5 |λc
max| 3.48 × 10−1 ∗ 2.43 × 101
Table 1: Results for Example 1
Example 2. In this example, we consider a linear descriptor system (E, A, B, C) with E = Tn and R = In.
System matrices A, B, C are randomly generated with entries of A distributed normally in [−5, 5], and entries of B and C distributed normally in [−1, 1]. We set rank(B) = rank(C) = dn
2e (the nearest integer ≥ n 2) for
n = 5, 15, 25, 35, 45. Note that the matrix E becomes nearly singular for large values of n and its condition number varies from O(101) to O(1014). For n = 35, one of the closed-loop eigenvalues achieved by QR-SWAP
lies outside the unit circle, with modulus equals 3.1413. When n = 45, four of the closed-loop eigenvalues achieved by QR-SWAP lie outside the unit circle. The numerical results are reported in Table 2.
n cond(E) G-SDA dare QR-SWAP
NRes 9.13 × 10−17 6.19 × 10−16 1.70 × 10−14
5 2.9 × 101 no. ite. 6 - 8
no. stab. ev. 5 5 5
|λc
max| 6.32 × 10−1 6.32 × 10−1 6.32 × 10−1
NRes 2.25 × 10−16 7.90 × 10−16 8.12 × 10−12
15 9.5 × 104 no. ite. 5 - 6
no. stab. ev. 15 15 15
|λc
max| 1.76 × 10−1 1.76 × 10−1 1.76 × 10−1
NRes 1.04 × 10−16 ∗ 1.09 × 10−12
25 1.7 × 108 no. ite. 5 - 6
no. stab. ev. 25 ∗ 25
|λc
max| 1.72 × 10−1 ∗ 1.72 × 10−1
NRes 2.23 × 10−16 ∗ 5.08 × 10−13
35 2.4 × 1011 no. ite. 5 - 6
no. stab. ev. 35 ∗ 34
|λc
max| 2.51 × 10−1 ∗ 3.14 × 10 0
NRes 3.11 × 10−16 ∗ 1.18 × 10−13
45 3.3 × 1014 no. ite. 5 - 6
no. stab. ev. 45 ∗ 41
|λc
max| 5.05 × 10−1 ∗ 2.65 × 100
Example 3. Let E ∈ Rn×n be the Frank matrix E = n n − 1 n − 2 · · · 2 1 n − 1 n − 1 n − 2 · · · 2 1 0 n − 2 n − 2 · · · 2 1 0 0 n − 3 . .. ... ... .. . ... . .. . .. ... ... .. . ... . .. 2 1 0 0 · · · 0 1 1 , and let A = trid(20, −10, −10) ∈ Rn×n, R = In.
The control matrix B and output matrix C are randomly generated with entries distributed normally in [−1, 1]. We set rank(B) = rank(C) = dn2e. Note that the matrix E becomes nearly singular for increasing values of n and its condition number varies from O(1) to O(1014). The numerical results are reported in
Table 3 for n = 5, 8, 11, 13, 16. For n = 13, 16, some closed-loop eigenvalues achieved by QR-SWAP lies outside the unit circle, with moduli up to 2.7023 and 7.6843, respectively.
n cond(E) G-SDA dare QR-SWAP
NRes 4.14 × 10−17 4.70 × 10−16 2.71 × 10−14
5 6.5 × 102 no. ite. 6 - 7
no. stab. ev. 5 5 5
|λc
max| 3.62 × 10−1 3.62 × 10−1 3.62 × 10−1
NRes 3.90 × 10−16 2.28 × 10−15 8.70 × 10−10
8 2.8 × 105 no. ite. 6 - 7
no. stab. ev. 8 8 8
|λc
max| 4.79 × 10−1 4.78 × 10−1 4.79 × 10−1
NRes 9.81 × 10−17 ∗ 5.10 × 10−12
11 3.3 × 108 no. ite. 7 - 7
no. stab. ev. 11 ∗ 11
|λc
max| 5.68 × 10−1 ∗ 5.68 × 10−1
NRes 7.79 × 10−17 ∗ 6.24 × 10−11
13 5.9 × 1010 no. ite. 6 - 7
no. stab. ev. 13 ∗ 12
|λc
max| 4.46 × 10−1 ∗ 2.70 × 10 0
NRes 2.39 × 10−16 ∗ 2.26 × 10−14
16 2.3 × 1014 no. ite. 6 - 7
no. stab. ev. 16 ∗ 15
|λc
max| 7.36 × 10−1 ∗ 7.68 × 100
Table 3: Results for Example 3
Example 4. This example is modified from Example 15 in [6], which was presented originally in [15, Example 3]. Here we consider the G-DARE defined by
E = diag(1, 10−1, 10−2, · · · , 10−(n−1)) ∈ Rn×n, A = trid(0, 0, 1) ∈ Rn×n, and
If E = diag(e11, e22, · · · , enn), then it is easily seen that the stabilizing s.p.s.d. solution is
X = diag(x1, x2, · · · , xn),
where x1 = 1/e211 and xj = (xj−1+ 1)/e2jj, for j = 2, . . . , n. Note that the closed-loop eigenvalues are all
zero. The numerical results are given in Table 4.
n cond(E) G-SDA dare QR-SWAP
NRes 1.52 × 10−16 2.05 × 10−15 2.95 × 10−17
2 10 no. ite. 2 - 8
no. stab. ev. 2 2 2
|λc max| 0.00 × 100 1.52 × 10−33 0 NRes 2.32 × 10−16 2.14 × 10−5 7.75 × 10−13 4 103 no. ite. 3 - 7 no. stab. ev 4 4 4 |λc max| 0.00 × 10 0 2.16 × 10−7 5.51 × 10−7 NRes 8.15 × 10−17 ∗ 6.58 × 10−1 6 105 no. ite. 4 - 8
no. stab. ev. 6 ∗ 6
|λc
max| 0.00 × 100 ∗ 2.74 × 10−3
NRes 3.85 × 10−16 ∗ 9.99 × 10−1
8 107 no. ite. 4 - 11
no. stab. ev. 8 ∗ 6
|λc
max| 0.00 × 100 ∗ 1.27 × 105
NRes 1.95 × 10−16 ∗ 9.96 × 10−1
10 109 no. ite. 5 - 13
no. stab. ev. 10 ∗ 6
|λc
max| 0.00 × 100 ∗ 5.54 × 107
Table 4: Results for Example 4
Example 5. Here we consider a linear descriptor system (E, A, B, C) with E = Tn and R = TmTmT, where
m = rank(B). The matrices A, B, C are randomly generated with entries of A distributed normally in [−5, 5], and entries of B and C distributed normally in [−1, 1]. We set rank(B) = rank(C) = dn2e for n = 5, 15, 25, 35, 45. Note that the matrices E and R become nearly singular for increasing n and their condition numbers vary from O(101) to O(1015). For n = 35, one of the closed-loop eigenvalues achieved by
QR-SWAP lies outside the unit circle, with modulus equals 20.9885. When n = 45, four of the closed-loop eigenvalues achieved by QR-SWAP lie outside the unit circle. The numerical results are reported in Table 5.
5
Conclusions
We have develop the G-SDA algorithm which solves G-DAREs with conditioned R and E. Inversions of ill-conditioned matrices are circumvented via well selected swapping of matrix products and applications of the SMWF. Numerical results show that the G-SDA is competitive with QR-SWAP and dare, out-performing the other algorithms in the selected set of test examples. The advantage of our structure-preserving algorithm is evident from the absence of unstable closed-loop eigenvalues at the end of the iterative process, contrasting results from QR-SWAP. The corresponding solution from QR-SWAP may be accurate in the sense of a small residual, it is obviously useless in the sense of stabilizing the closed-loop system. The MATLAB command
n cond(E) cond(R) G-SDA dare QR-SWAP NRes 1.97 × 10−16 3.59 × 10−16 5.39 × 10−14 5 2.9 × 101 2.9 × 101 no. ite. 6 - 8
no. stab. ev. 5 5 5
|λc
max| 6.66 × 10−1 6.66 × 10−1 6.66 × 10−1
NRes 7.76 × 10−17 1.10 × 10−15 7.29 × 10−12 15 9.5 × 104 1.4 × 105 no. ite. 5 - 6
no. stab. ev. 15 15 15
|λc
max| 2.40 × 10−1 2.40 × 10−1 2.40 × 10−1
NRes 4.84 × 10−16 ∗ 2.74 × 10−12 25 1.7 × 108 4.2 × 108 no. ite. 6 - 6
no. stab. ev. 25 ∗ 25
|λc
max| 2.21 × 10−1 ∗ 2.21 × 10−1
NRes 2.25 × 10−16 ∗ 2.95 × 10−12 35 2.4 × 1011 8.6 × 1011 no. ite. 6 - 6
no. stab. ev. 35 ∗ 34
|λc
max| 3.12 × 10−1 ∗ 2.10 × 101
NRes 4.96 × 10−16 ∗ 4.35 × 10−14 45 3.3 × 1014 1.5 × 1015 no. ite. 6 - 6
no. stab. ev. 45 ∗ 41
|λc
max| 8.60 × 10−1 ∗ 4.11 × 10 1
Table 5: Results for Example 5
dare failed frequently for many ill-conditioned examples. Apart from having superior accuracy, convergence and structure-preserving properties, the operation count (per iteration) for G-SDA is a small fraction of those for the other algorithms, analogous to the superiority of the SDA for DAREs [8]. This efficiency is the consequence of the fact that the G-SDA operates in Rn×n while QR-SWAP and dare works with matrices
of higher dimensions.
References
[1] B. Anderson, Second-order convergent algorithms for the steady-state Riccati equation, Int. J. Control, 28 (1978), pp. 295–306.
[2] Z. Bai, J. Demmel, and M. Gu, An inverse free parallel spectral divide and conquer algorithm for nonsymmetric eigenproblems, Numer. Math., 76 (1997), pp. 279–308.
[3] P. Benner, Contributions to the numerical solutions of algebraic Riccati equations and related eigenvalue problems, PhD Dissertation, Fakult¨at f¨ur Mathematik, TU Chemnitz-Zwickau, Chemnitz, Germany, 1997.
[4] P. Benner and R. Byers, Evaluating products of matrix pencils and collapsing matrix products, Num. Lin. Alg. Appl., 8 (2001), pp. 357–380.
[5] P. Benner, R. Byers, R. Mayo, E. S. Quintana-Orti, and V. Hern´andez, Parallel algorithms for LQ optimal control of discrete-time periodic linear systems, J. Para. Distr. Comp., 62 (2002), pp. 306– 325.
[6] P. Benner, A. J. Laub, and V. Mehrmann, A collection of benchmark examples for the nu-merical solution of algebraic Riccati equations II: Discrete-time case, Tech. Rep. SPC 95 23, Fakult¨at f¨ur Mathematik, TU Chemnitz–Zwickau, 09107 Chemnitz, FRG, 1995. Available from http://www.tu-chemnitz.de/sfb393/spc95pr.html.
[7] A. Bunse-Gerstner, V. Mehrmann, and N. K. Nichols, Regularization of descriptor systems by derivative and proportional state feedback, SIAM J. Matrix Anal. Appl., 13 (1992), pp. 46–67.
[8] E. K.-W. Chu, H.-Y. Fan, C.-S. Wang, and W.-W. Lin, A structure-preserving doubling algorithm for periodic descrete-time algebraic Riccati equations, preprint 2002-18, NCTS, National Tsing Hua University, Hsinchu 300, Taiwan, 2003.
[9] M. Kimura, Convergence of the doubling algorithm for the discrete-time algebraic Riccati equation, Int. J. Syst. Sci., 19 (1988), pp. 701–711.
[10] A. J. Laub, A Schur method for solving algebraic Riccati equations, IEEE Trans. Auto. Control, 24 (1979), pp. 913–921.
[11] A. J. Laub, Algebraic aspects of generalized eigenvalue problems for solving Riccati equations, in Com-putational and Combinatorial Methods in Systems Theory, C. I. Byrnes and A. Lindquist, eds., Elsevier (North-Holland), 1986, pp. 213–227.
[12] A. J. Laub, Invariant subspace methods for the numerical solution of Riccati equations, in The Riccati Equation, S. Bittanti, A. J. Laub, and J. C. Willems, eds., Springer-Verlag, Berlin, 1991, pp. 163–196. [13] MathWorks, MATLAB user’s guide (for UNIX Workstations), The Math Works, Inc., 1992. [14] V. Mehrmann, The Autonomous Linear Quadratic Control Problem, Springer-Verlag, 1991. [15] T. Pappas, A. J. Laub, and N. R. Sandell, On the numerical solution of the discrete-time algebraic
Riccati equation, IEEE Trans. Auto. Control, 25 (1980), pp. 631–641.
