Robust formulas for optimal H∞ controllers [original]

Technische Universit¨at Berlin

Institut f¨ur Mathematik

Robust Formulas for H∞Controllers

Peter Benner, Ralph Byers, Philip Losse,

Volker Mehrmann, and Hongguo Xu

Preprint 20-2010

Preprint-Reihe des Instituts f¨ur Mathematik

Technische Universit¨at Berlin

http://www.math.tu-berlin.de/preprints

Report 20-2010 December 2010

Robust formulas for optimal H∞controllers

Peter Benner a,1, Ralph Byers b, Philip Losse c,1, Volker Mehrmann c,1,

Hongguo Xu d,2

aFakult¨at f¨ur Mathematik, TU Chemnitz, D-09107 Chemnitz, FRG.

bDeceased, last address: Department of Mathematics, University of Kansas, Lawrence, Kansas, USA; correspondence to

Hongguo Xu.

cInstitut f¨ur Mathematik MA 4-5, TU Berlin, Straße des 17. Juni 136, D-10623 Berlin, FRG.

dDepartment of Mathematics, University of Kansas, Lawrence, Kansas, USA.

Abstract

We present formulas for the construction of optimal H∞controllers that can be implemented in a numerically robust way. We

base the formulas on the γ-iteration developed in [6]. The controller formulas proposed here avoid the solution of algebraic

Riccati equations with their problematic matrix inverses and matrix products. They are also applicable in the neighborhood

of the optimal γ, where the classical formulas may call for the inverse of singular or ill-conditioned matrices. The advantages

of the new formulas are demonstrated by several numerical examples.

Key words: H∞control; controller design; optimal controller; CS decomposition; Lagrangian subspaces; even pencil.

1 Introduction

The optimal infinite-horizon output (or measurement)

feedback H∞control problem is one of the central tasks

in robust control, see, e.g., [13,18,21,22], but the devel-

opment of robust numerical methods for the H∞con-

trol is unusually difficult [20]. The classic γ-iteration of-

ten used in optimal H∞control computations encoun-

ters several finite precision arithmetic hazards that of-

ten limit its accuracy as a numerical method. A new nu-

merical method for the γ-iteration suggested in [6] has

significantly better robustness in the presence of round-

Email addresses: [email protected]

(Peter Benner), [email protected] (Philip Losse),

[email protected] (Volker Mehrmann),

[email protected] (Hongguo Xu).

1These authors were partially supported by Deutsche

Forschungsgemeinschaft, Research Grant Me 790/16-1,

Be 2174/6-1.

2This author was partially supported by National Science

Foundation grant 0314427, and the University of Kansas

General Research Fund allocation # 2301717.

ing errors. Based on this approach, this paper proposes

a numerical method for the implementation of the asso-

ciated optimal controllers. Note that another variant of

controller formulas based on the the γ-iteration from [6]

is suggested in [15]. Our approach differs in the deriva-

tion and form of the controller formulas. Moreover, we

have implemented our formulas and we will present nu-

merical results obtained with these formulas.

Consider the linear control system

˙x=Ax +B1w+B2u, x(t0) = x0,

z=C1x+D11w+D12u, (1)

y=C2x+D21w+D22u,

where A∈Rn×n,Bi∈Rn×mi,Ci∈Rpi×n, and Dij ∈

Rpi×mjfor i, j = 1,2. (By Rn×kwe denote the set of

real n×kmatrices.) As usual, see [13,22], we assume

p1≥m2and m1≥p2. In this system, x(t)∈Rnis

the state vector, u(t)∈Rm2is the control input vec-

tor, and w(t)∈Rm1is an exogenous input that may in-

clude noise, linearization errors and unmodeled dynam-

Preprint submitted to Automatica 4 December 2010

ics. The vector y(t)∈Rp2contains measured outputs,

while z(t)∈Rp1is a regulated output or error.

The optimal H∞control problem: Determine a dy-

namic controller

ˆx=ˆ

Aˆx+ˆ

B1y+ˆ

B2ˆu,

u=ˆ

C1ˆx+ˆ

D11y+ˆ

D12 ˆu,

ˆy=ˆ

C2ˆx+ˆ

D21y,

(2)

with ˆ

A∈RN×N,ˆ

B1∈RN×p2,

B2∈RN×q1,ˆ

C1∈Rm2×N,ˆ

C2∈Rq2×N,ˆ

D11 ∈Rm2×p2,

D12 ∈Rm2×q1,ˆ

D21 ∈Rq2×p2such that the closed-loop

system resulting from the combined system of (1) and

(2),

(1) is internally stable, i.e., the solution of the system

with w≡0 is asymptotically stable, and

(2) the closed-loop transfer function Tzw from wto z

is minimized in the H∞norm.

For a matrix valued rational function F(s) that is ana-

lytic in the open right-half plane, the H∞norm is given

by kFk∞= supω∈Rσmax[F(ıω)], where σmax[F(ıω)] de-

notes the maximal singular value of the matrix F(ıω).

If F(s) is the transfer function of a control system with

noise or disturbance inputs, then kFk∞is a measure of

the worst case influence of the disturbances on the out-

put. The solution of this problem is usually approached

via the modified optimal H∞control problem:

The modified optimal H∞control problem: Let Γ

be the set of numbers γ > 0 for which there exists an in-

ternally stabilizing dynamic controller (2) such that the

closed loop transfer function Tzw satisfies γ > kTzwk∞.

Determine γmo = inf Γ.

Because there may be no dynamic controller that leads

to a transfer function that actually achieves H∞norm

equal to γmo, in general, one must use a controller whose

transfer function has larger H∞norm, i.e., an inter-

nally stabilizing dynamic controller such that the closed

loop transfer function satisfies kTzwk∞< γ for some

γ > γmo. Such a controller is usually called a suboptimal

controller. The γ-iteration is the iterative root finding

process of determining an approximation to γmo. Classi-

cal numerical methods for determining γmo are based on

the solution of Riccati equations or Lagrangian invari-

ant subspaces, see [11,13,14,18,22] and are implemented

in software packages like MATLAB R

or SLICOT, [7–9].

A more robust method for carrying out the γ-iteration

has recently been proposed in [6].

Once a sufficiently accurate approximation to γmo is de-

termined, a suboptimal controller can be constructed

using the mathematically correct, but numerically haz-

ardous formulas suggested in [14,22] which we recall in

Section 2, or by the more robust formulas that we present

in Section 3. We will demonstrate the quality of the new

formulas with several numerical examples in Section 4

and give some final remarks in Section 5.

2 Preliminaries

The formulas for designing optimal controllers are quite

technical and only hold under some suitable assump-

tions. In this section we review the classical formulas

and the assumptions under which H∞norm calculations

typically operate.

A typical set of assumptions for the solution of the mod-

ified optimal H∞control is as follows [13,14,18,22]:

A1. The pair (A, B2) is stabilizable and the pair

(A, C2) is detectable, i.e., rank[A−λI, B2] =

rank[AT−λI, CT

2] = nfor all λ∈Cwith Re λ≥0.

A2. D22 = 0 and both D12 and D21 have full rank.

A3. The matrix hA−ıωI

D12 ihas full column rank for

all real ω.

A4. The matrix hA−ıωI

D21 ihas full row rank for all

real ω.

Remark 2.1 The requirement that D22 = 0 (Assump-

tion A2) is for convenience. It is not a fundamental re-

striction, since systems that have a direct link from input

to output, i.e., for which D22 6= 0, can be synthesized by

first studying the problem without this term, see [22].

Following the notation in [22], we introduce the following

two symmetric matrices formed from the matrices Dij

and a parameter γ∈R,

RH(γ) := "DT

12 #hD11 D12 i−"γ2Im10

0 0 #,

RJ(γ) := "D11

D21 #hDT

11 DT

21 i−"γ2Ip10

0 0 #.

(3)

These matrices play an essential role in the theory of op-

timal H∞control problems, see [14,22], and the classical

numerical methods require both RH(γ) and RJ(γ) to

be nonsingular. Under Assumption A2, there exist only

a finite number of nonnegative values γfor which (at

least) one of the matrices RH(γ) and RJ(γ) is singular.

Let ˆγbe the largest γvalue for which this is the case. If

D11 = 0, then ˆγ= 0; otherwise, ˆγis typically positive.

Note that by definition, γmo >ˆγ.

Let

D12 =U12 "0

Σ12 #VT

12, D21 =V21 h0 Σ21 iUT

21,(4)

be (slightly permuted) singular value decompositions

(see [12]) of D12 and D21 with real orthogonal matrices

U12,U21,V12,V21 and nonnegative diagonal matrices

Σ12, Σ21. The diagonal entries of Σ12 and Σ21 are the

singular values of D12 and D21, respectively. Then define

D11,¯

D12, and ¯

D21 in terms of D11,D12,D21 and (4) by

"D11 D12

D21 0#=

"U12 0

0V21Σ21 #"¯

D11 ¯

D12

D21 0#"UT

21 0

0 Σ12VT

12 #.

(5)

Note that by assumption D22 = 0 and the described

transformation does not change this. It follows from (4)

that ¯

D21 = [0, Ip2] and ¯

D12 =h0

Im2i. This induces a

finer partition of ¯

D11 so that

"¯

D11 ¯

D12

D21 0#=





D1D20

D3D4Im2

0Ip20







.(6)

Now, under Assumption A2, for ˆγas defined above, we

have

ˆγ= max σmax hD1D2i, σmax "D1

D3#!,

where σmax(M) denotes the maximal singular value of

the matrix M.

The classical approach for the computation of γmo, see,

e.g., [21,22], employs the solution of algebraic Riccati

equations (AREs). Consider the Hamiltonian matrices

H(γ) = "H1(γ)H2(γ)

H3(γ)−H1(γ)T#

="A0

−CT

1C1−AT#(7)

−B1

−CT

1D11

−CT

1D12 R−1

H(γ)DT

11C1

12C1

2,

J(γ) = "J1(γ)J2(γ)

J3(γ)−J1(γ)T#

="AT0

−B1BT

1−A#(8)

−CT

−B1DT

21 R−1

J(γ)D11BT

D21BT

C2,

and the associated γ-dependent AREs

H1(γ)XH(γ) + XH(γ)H1(γ)T

+XH(γ)H2(γ)XH(γ)−H3(γ)=0,(9)

and

J1(γ)XJ(γ) + XJ(γ)J1(γ)T

+XJ(γ)J2(γ)XJ(γ)−J3(γ)=0.(10)

Classically, one computes the unique symmetric positive

semidefinite (stabilizing) symmetric solutions XH(γ)

and XJ(γ) of (9), (10), respectively, or what is more

numerically stable, invariant subspaces of the associ-

ated Hamiltonian matrices, see [22, Ch. 16–17]. The

latter approach determines symmetric matrices XH, XJ

matrices such that

H(γ)"In

XH#="In

XH#TH, J(γ)"In

XJ#="In

XJ#TJ,

for some n×nmatrices THand TJ, respectively, with

all their eigenvalues in the open left half complex plane.

Remark 2.2 The columns of the matrices hIn

XHiand

hIn

XJiform unique Lagrangian invariant subspaces. Un-

der some further assumptions, see [10,19], such unique

Lagrangian invariant subspaces still exist, even when

eigenvalues are on the imaginary axis.

In terms of XH, XJ(we leave off the dependency on γin

the following) and the original data we then define the

matrices

F=−R−1

H "DT

12 #C1+"BT

2#XH!=: "F1

F2#,

(11a)

L=−B1hDT

11 DT

21 i+XJhCT

1CT

2iR−1

J(11b)

=: hL1L2i,

Z=In−γ−2XJXH−1,(11c)

where RHand RJare defined in (3). Once γmo, the

optimal value of γ, has been determined, then for

all γ > γmo, one has (see [6,22]) that RH(γ), RJ(γ)

are nonsingular; the matrices XH, and XJexist; and

γ2> ρ(XJXH), where ρ(XJXH) is the spectral radius

of XJXH. Therefore, for every γ > γmo, the matrices

F, L, Z are well defined.

Then, for a given number γ≥γmo, a suboptimal con-

troller (2) is usually constructed by using the following

formulas ([22, Theorem 17.1]).

(a)ˆ

D11 =−V12Σ−1

12 ·

D3DT

1γ2I−D1DT

1−1D2+D4Σ−1

21 VT

21,

(b)ˆ

D12 ˆ

12 =V12Σ−1

12 ·

Im2−D3γ2I−DT

1D1−1DT

3Σ−1

12 VT

12,

(c)ˆ

21 ˆ

D21 =V21Σ−1

21 ·

Ip2−DT

2γ2I−D1DT

1−1D2Σ−1

21 VT

21,

(d)ˆ

B2=Z(B2+L1D12)ˆ

D12,

(e)ˆ

B1=Z[(B2+L1D12)ˆ

D11 −L2],

(f)ˆ

C2=−ˆ

D21(C2+D21F1),

(g)ˆ

C1=F2−ˆ

D11(C2+D21F1),

(h)ˆ

A=A+hB1B2iF−ˆ

B1(C2+D21F1).

(12)

The basis for the robust method derived in [6] to com-

pute γmo is to avoid all inversions, matrix products and

sums, as well as solutions to AREs, that are used in

the classical γ-iteration. For this, the eigenvalue prob-

lems for the Hamiltonian matrices H(γ) and J(γ) are

first replaced by generalized eigenvalue problems for the

following two even (skew-symmetric/symmetric) matrix

pencils that only contain original data from (1):

λN −MH(γ) :=







0In000

−In0 0 0 0

0 0 0 0 0







−







0−AT0 0 −CT

−A0B1B20

0BT

1γ2Im10DT

0BT

20 0 DT

−C10D11 D12 Ip1







λN −MJ(γ) :=







0In000

−In0 0 0 0

0 0 0 0 0







−







0−A0 0 −B1

−AT0CT

1CT

0C1γ2Ip10D11

0C20 0 D21

−BT

10DT

11 DT

21 Im1







Denote by ˆγIthe largest γvalue for which at least one

of these pencils has a purely imaginary eigenvalue. It

has been shown in [6] that if the assumptions A1–A4 are

satisfied, then for all γ > ˆγI, the pencils λN−MH(γ) and

λN −MJ(γ) each have a unique n-dimensional deflating

subspace corresponding to the eigenvalues in the open

left half plane and for γ→ˆγIthere still exists a unique

n-dimensional deflating subspace corresponding to the

eigenvalues in the closed left half plane.

If the columns of the matrices

QH=





n QH,1

n QH,2

m1QH,3

m2QH,4

p1QH,5







, QJ=





n QJ,1

n QJ,2

p1QJ,3

p2QJ,4

m1QJ,5







span these unique deflating subspaces, then the columns

of the submatrices

"QH,1

QH,2#,"QJ,1

QJ,2#

span the desired Lagrangian invariant subspaces of the

Hamiltonian matrices H(γ) in (7) and J(γ) in (8), re-

spectively, i.e., QT

H,1QH,2=QT

H,2QH,1and QT

J,1QJ,2=

J,2QJ,1, see [6,10,19]. Furthermore, the symmetric pos-

itive semidefinite (stabilizing) solutions of the AREs (9)

and (10), if they exist, can be expressed as

XH=QH,2Q−1

H,1, XJ=QJ,2Q−1

J,1.

In order to avoid explicitly forming XH, XJ, in [6] the

following technique was introduced. Let the columns of

"XH,1

XH,2#,"XJ,1

XJ,2#(13)

form orthonormal bases of the Lagrangian invariant sub-

spaces hQH,1

QH,2i,hQJ,1

QJ,2i, respectively. These may be de-

termined by QR factorization or the modified Gram-

Schmidt process [12]. Note again that we have dropped

the explicit γdependency in XH,i, XJ,i,i= 1,2.

Introduce the symmetric matrix

Y(γ) := "γXT

H,2XH,1XT

H,2XJ,2

J,2XH,2γXT

J,2XJ,1#

="XT

H,20

0XT

J,2#"γXH,1XJ,2

XH,2γXJ,1#.

If γ≤ˆγI, then the pencils λN −MH(γ) and λN −MJ(γ)

may or may not have a unique n-dimensional deflating

subspace corresponding to the eigenvalues in the closed

left-half plane. If for a particular γ≤ˆγI, such unique

deflating subspaces do not exist, then Y(γ) is not defined.

The optimal value γmo is then determined in [6] and [22,

Theorem 16.16] by detecting a rank change in the ma-

trix Y(γ). For this and for the formulas in Section 3 we

need the following well-known lemma on the CS decom-

position.

Lemma 2.3 [17] If X1, X2∈Rn,n and the columns of

hX1

X2iform an orthonormal basis of a Lagrangian sub-

space, i.e., XT

1X2=XT

2X1, then there exist orthogonal

matrices U∈Rn×nand V∈Rn×nsuch that UTX1V=

Cand UTX2V=Sare both diagonal and C2+S2=I.

If for a given value of γwe apply Lemma 2.3 to hXH,1

XH,2i,

hXJ,1

XJ,2i, then we get

HXH,1VH=CH

=: 



rHkHn−tH

rH0 0 0

kH0 ΣH0

n−tH0 0 I

,(14a)

HXH,2VH=SH

=: 



rHkHn−tH

rHI0 0

kH0 ∆H0

n−tH0 0 0 

,(14b)

JXJ,1VJ=CJ

=: 



rJkJn−tJ

rJ0 0 0

kJ0 ΣJ0

n−tJ0 0 I

,(14c)

JXJ,2VJ=SJ

=: 



rJkJn−tJ

rJI0 0

kJ0 ∆J0

n−tJ0 0 0 

,(14d)

where kH+rH=tH,kJ+rJ=tJ, ΣH, ∆H, ΣJand

∆Jare diagonal, nonsingular and satisfy Σ2

H+ ∆2

H=I

and Σ2

J+ ∆2

J=I.

The following theorem then is the basis of the variant

γ-iteration introduced in [6].

Theorem 2.4 [6] For all γ > γmo, the matrix Y(γ)is

positive semidefinite and rank Y(γ) = kH+kJis con-

stant. For all ˆγ < γ < γmo, either Y(γ)is not defined,

or rank Y(γ)< kH+kJ, or Y(γ)is not positive semidef-

inite.

The problem of finding γmo thus reduces to the prob-

lem of finding the largest value of γ(≥ˆγ) at which Y(γ)

changes rank or fails to exist. This can be done, see [6],

by applying a one-variable root-finding procedure to the

eigenvalues of Y(γ). In this way a good approximation

of γmo can be determined without explicitely forming

the Hamiltonian matrices H(γ), J(γ), while still using

structure preserving methods [5]. Once γmo has been de-

termined, it remains to construct a corresponding con-

troller.

In order to avoid the hazards of solving ill-conditioned

systems of equations, we will transform the formulas for

the optimal controllers (12). For this, we will make fre-

quent use of the following refactorization, that can be

computed without forming explicit inverses or matrix

products, [1–4].

Proposition 2.5 Given a matrix product MP −1with

M∈Rη×µand P∈Rµ×µ, there exist P ∈ Rη×ηand

M ∈ Rη×µsuch that

MP−1=P−1M,(15)

where the rows of [P,M]span the η-dimensional left null

space of h−M

Pi.

To construct a refactorization as in (15), observe that

MP−1=P−1Mif and only if PM=MPor, equiva-

lently,

[P,M]"−M

P#= 0.

So, each basis of the left nullspace gives rise to such a

refactorization. A convenient way to calculate such a left

null space and refactorization that was used in [1] (and

in this paper) is to use a QR factorization of h−M

Pi, i.e.

"Q11 Q12

Q21 Q22 #"R

0#="−M

with R∈Rµ×µ,Q11 ∈Rη×µ,Q21 ∈Rµ×µ,Q12 ∈Rη×η,

and Q22 ∈Rµ×η. A suitable left null space is the row

space of [QT

12, QT

22] — one may use P=QT

12,M=QT

22.

3 Formulas for the (Sub)optimal Controller

This section discusses the construction of (sub)optimal

H∞-controllers as in (2) from γmo or an approximation

γ > γmo to it. The new approach is to reorganize the

controller formulas (12) that use the computed data in

(11a)–(11c) into a descriptor system form by avoiding

numerical hazards where possible.

Using (4) and (5), we introduce

WH="U21 0

0V12Σ−1

12 #, WJ="U12 0

0V21Σ−1

21 #,(16)

and

B1=B1U21,¯

B2=B2V12Σ−1

12 ,

C1=UT

12C1,¯

C2= Σ−1

21 VT

21C2.(17)

Note that in these formulas the inverses of the diago-

nal matrices Σ12 and Σ21 occur. These inverses can be

formed in a numerically stable way. Ill-conditioning of

these matrices, however, indicates that the resulting ro-

bust controller may be itself sensitive to small perturba-

tions.

Using (5), we define ¯

RHand ¯

RJby

RH=WT

HRHWH="¯

12 #h¯

D11 ¯

D12 i−"γ2Im10

0 0 #,

RJ=WT

JRJWJ="¯

D11

D21 #h¯

11 ¯

21 i−"γ2Ip10

0 0 #,

and let hXH,1

XH,2i,hXJ,1

XJ,2ibe the orthonormal matrices

defined in (13). Using (14a)–(14d) in Lemma 2.3, the

unique symmetric positive semidefinite (stabilizing)

ARE solutions XH, XJcan be expressed as ([22, Theo-

rem 17.1])

XH=XH,2X−1

H,1=UHSHC−1

HUT

H=UHC−1

HSHUT

XJ=XJ,2X−1

J,1=UJSJC−1

JUT

J=UJC−1

JSJUT

Here, again, ill-conditioning of the diagonal matrices

CH, CJmay indicate high sensitivity of the controller.

See [6] for more details.

With the quantities defined in (5) and (17) we introduce

FM=−"¯

12 #¯

C1UHCH+"¯

2#UHSH,

LM=−CJUT

J¯

B1h¯

11 ¯

21 i+SJUT

Jh¯

1¯

2i.

The matrices F,L, and Zin (11a)–(11c) can then be

expressed as

F=WHˆ

F, where ˆ

F=¯

R−1

HFMC−1

HUT

H=: 

m1ˆ

m2ˆ

F2,

(18a)

L=ˆ

LWT

J,where ˆ

L=UJC−1

JLM¯

R−1

J=: h

p1p2

L1ˆ

L2i,

(18b)

Z=UHCHˆ

Z−1CJUT

J,where

Z=CJUT

JUHCH−γ−2SJUT

JUHSH.(18c)

Then the coefficient matrices ˆ

A,ˆ

B1,ˆ

B2,ˆ

C1,ˆ

C2in the con-

troller (2) can be expressed as

B2=Z(¯

B2+ˆ

L1¯

D12)Σ12VT

12 ˆ

D12,

B1=Z[( ¯

B2+ˆ

L1¯

D12)(Σ12VT

12 ˆ

D11V21Σ21)−ˆ

L2]Σ−1

21 VT

21,

C2=−ˆ

D21V21Σ21(¯

C2+¯

D21 ˆ

F1),

C1=V12Σ−1

12 [ˆ

F2−(Σ12VT

12 ˆ

D11V21Σ21)( ¯

C2+¯

D21 ˆ

F1)],

A=A+h¯

B1¯

B2iˆ

F−ˆ

B1V21Σ21(¯

C2+¯

D21 ˆ

F1).

Note that in these formulas still some unwanted inverses

arise which we like to avoid. To do this, using Propo-

sition 2.5, we can refactor the products ¯

R−1

HFMand

LM¯

R−1

Jas

R−1

HFM=FM¯

R−1

H, LM¯

R−1

J=¯

R−1

JLM.

Using the same notation as in [22], we partition LMand

FMas

LM= [

p1−m2m2p2

LM,11,∞LM,12,∞LM,2,∞],

FM="

m1−p2FM,11,∞

p2FM,12,∞

m2FM,2,∞#.

Then, in (18a)–(18c) we obtain new factors

F1="FM,11,∞

FM,12,∞#¯

R−1

HC−1

HUT

F2=FM,2,∞¯

R−1

HC−1

HUT

L1=UJC−1

J¯

R−1

JhLM,11,∞LM,12,∞i,

L2=UJC−1

J¯

R−1

JLM,2,∞.

Let

D1=˜

U1"Σ10

0 0 #˜

be the singular value decomposition of D1in (6). Define

D1,¯

D2,¯

D3,¯

D4, ∆13, ∆23, ∆31 and ∆32 in terms of ˜

U1,

V1,D1,D2,D3and D4in (6) by

"¯

D1¯

D3¯

D4#="˜

0I#"D1D2

D3D4#"˜

V10

0I#

=





Σ10 ∆13

0 0 ∆23

∆31 ∆32 D4







If d1= rank D1, then Σ1∈Rd1,d1, ∆13 ∈Rd1,p2, ∆23 ∈

Rp1−m2−d1,p2, ∆31 ∈Rm2,d1and ∆32 ∈Rm2,m1−p2−d1.

Using these factorizations, we can rewrite

D11 := −(D3DT

1γ2I−D1DT

1−1D2+D4)

=−(¯

D3¯

1(γ2I−¯

D1¯

1)−1¯

D2+D4)

and define ˜

D12 and ˜

D21 via the Cholesky factorizations

D12 ˜

12 =Im2−D3(γ2I−DT

1D1)−1DT

=Im2−¯

D3(γ2I−¯

1¯

D1)−1¯

21 ˜

D21 =Ip2−DT

2(γ2I−D1DT

1)−1D2

=Ip2−¯

2(γ2I−¯

D1¯

1)−1¯

D2.

By these factorizations and using M+to denote the

Moore-Penrose pseudoinverse of M, we obtain

(a)˜

D11 =−∆31Σ1(γ2I−Σ2

1)+∆13 +D4,

(b)˜

D12 ˜

12 =Im2−∆31(γ2I−Σ2

1)+∆T

31 −γ−2∆32∆T

32,

(c)˜

21 ˜

D21 =Ip2−∆T

13(γ2I−Σ2

1)+∆13 −γ−2∆T

23∆23.

(19)

Note that it is not necessary to use Cholesky factoriza-

tions, mathematically, any factorizations satisfying (19)

(b)–(c) can be used for ˜

D12 and ˜

D21.

In (19) we have replaced the inverses on the diagonal

matrix γ2I−Σ2

1by Moore-Penrose generalized inverses.

This allows to use the formulas even when for some γ

value the matrix becomes singular.

By using (19) and the forms ¯

D12 =h0

Im2iand ¯

D21 =

[0, Ip2], we then have the reformulation of (11c) as

Z=UHCHˆ

Z−1CJUT

where ˆ

Z=CJUT

JUHCH−γ−2SJUT

JUHSH, and (12)

becomes

(a)ˆ

D11 =V12Σ−1

12 ˜

D11Σ−1

21 VT

21,where

D11 = ∆31Σ1(γ2I−Σ2

1)+∆13 +D4,

(b)ˆ

D12 =V12Σ−1

12 ˜

D12,where

D12 ˜

12 =Im2−∆31(γ2I−Σ2

1)+∆T

31 −γ−2∆32∆T

32,

(c)ˆ

D21 =˜

D21Σ−1

21 VT

21,where

21 ˜

D21 =Ip2−∆T

13(γ2I−Σ2

1)+∆13 −γ−2∆T

23∆23,

(d)ˆ

B2=Z(¯

B2+ˆ

L1¯

D12)˜

D12 =UHCHˆ

Z−1¯

R−1

J·

(¯

RJCJUT

J¯

B2+LM,12,∞)˜

D12,

(e)ˆ

B1=Z[( ¯

B2+ˆ

L1¯

D12)˜

D11 −ˆ

L2]Σ−1

21 VT

=UHCHˆ

Z−1¯

R−1

J[( ¯

RJCJUT

J¯

+LM,12,∞)˜

D11 − LM,2,∞]Σ−1

21 VT

21,

(f)ˆ

C2=−˜

D21(¯

C2+¯

D21 ˆ

F1)

=−˜

D21(¯

C2UHCH¯

RH+FM,12,∞

)¯

R−1

HC−1

HUT

(g)ˆ

C1=V12Σ−1

12 [ˆ

F2−˜

D11(¯

C2+¯

D21 ˆ

F1)] = V12Σ−1

12 ·

[FM,2,∞−˜

D11(¯

C2UHCH¯

RH+FM,12,∞)]·

R−1

HC−1

HUT

(h)ˆ

A=A+h¯

B1¯

B2iˆ

F−(ˆ

B1V21Σ21)( ¯

C2+¯

D21 ˆ

F1)

=nAUHCH¯

RH+h¯

B1¯

B2iFM−(ˆ

B1V21Σ21)·

(¯

C2UHCH¯

RH+FM,12,∞)}¯

R−1

HC−1

HUT

(20)

Theorem 17.1 in [22] states that the optimal controllers

of the H∞control problem are obtained by composing

ˆx=ˆ

Aˆx+ˆ

B1y+ˆ

B2ˆu,

u=ˆ

C1ˆx+ˆ

D11y+ˆ

D12 ˆu, (21)

ˆy=ˆ

C2ˆx+ˆ

D21y,

with a system

˜x=˜

A˜x+˜

Bˆy,

ˆu=˜

C˜x+˜

Dˆy,

whose transfer function Q(s) has H∞norm less than

γ. In any application, some such controllers are likely

to be more robust, or less expensive, or more elegant,

or otherwise more desirable than others. A convenient

choice is Q(s) = 0, the “central controller” discussed in

[22, Ch. 16–17] and [14].

This suggests the following procedure to avoid the re-

maining explicit matrix inverses. Refactor CHˆ

Z−1¯

R−1

CH(¯

RJˆ

Z)−1=E−1CH,

as described in Proposition 2.5. Set ˇx=¯

R−1

HC−1

HUT

Hˆx

and multiply the first equation in (21) by CH¯

RJˆ

ZC−1

HUT

to obtain the descriptor system

E˙

ˇx=ˇ

Aˇx+ˇ

B1y+ˇ

B2ˆu,

u=ˇ

C1ˇx+ˆ

D11y+ˆ

D12 ˆu, (22)

ˆy=ˇ

C2ˇx+ˆ

D21y,

where ˆ

D11,ˆ

D12, and ˆ

D21 are as in (20), and

(a)ˇ

E=ECH¯

RH,

(b)ˇ

B1=CH[¯

RJCJUT

J¯

B2+LM,12,∞˜

D11 −LM,2,∞]·

Σ−1

21 VT

21,

(c)ˇ

B2=CH¯

RJCJUT

J¯

B2+LM,12,∞˜

D12,

(d)ˇ

C1=V12Σ−1

12 ·

[FM,2,∞−˜

D11(¯

C2UHCH¯

RH+FM,12,∞)],

(e)ˇ

C2=−˜

D21(¯

C2UHCH¯

RH+FM,12,∞),

(f)ˇ

A=EUT

HAUHCH¯

RH+h¯

B1¯

B2iFM

−(ˇ

B1V21Σ21)( ¯

C2UHCH¯

RH+FM,12,∞).

(23)

A descriptor system expression for the central controller

as suggested in [22, Ch. 16–17] and [14] is then of the

form

E˙

ˇx=ˇ

Aˇx+ˇ

B1y

u=ˇ

C1ˇx+ˇ

D11y, (24)

where ˇ

D11 =ˆ

D11.

In this section we have presented new formulas for op-

timal H∞conrollers that avoid unnecessary numerical

hazards. In the next section we illustrate the quality of

the new formulas via several numerical examples.

4 Numerical examples

This section demonstrates the robustness of (22) in form-

ing the suboptimal control for γclose to γmo. Note that

in the presented examples for the computed values of

γmo typically the solutions to the associated AREs do

not exist, and as a consequence, the classical controller

formulas cannot be applied. This is also the reason, why

we do not present comparisons with other approaches.

Unfortunately, we were also not successful in compar-

ing our results for the examples presented here with the

formulas presented in [15], since our implementation of

these formulas did not produce correct results.

Example 4.1 This is Example 6.1 from [6] which first

appeared in [22, p. 461]:







−a0 1 −2 1

0−100 0 0 0

0 0 0 −2a a

0 0 0 0 1

0 0 0 3 2







, B1=













, B2=







−90







C1=



1 0 0 0 0

0 1 0 0 0



, D11 =





, D12 =





,

C2=h0 0 1 −2 1 i, D21 =h1i, D22 =h0i.

In this example, γmo is independent of the choice

of a. Using a= 1 and the variant γ-iteration de-

scribed in [6], our experimental program determined

γmo = 7.853923684022. With γ= 7.8541, our MAT-

LAB implementation returns the central suboptimal

controller (24) as







0.1843 −0.1737 −0.0414 −0.1018 −0.0341

−0.1460 0.1434 0.4586 −0.0903 −0.0900

0.1083 −0.0368 0.0465 0.0126 0.0453

0.2221 0.4605 −0.1595 0.0046 0.0315

0.2766 −0.0515 0.1694 0.1360 0.2079













−3.7597 9.6660 0.0414 −3.6630 8.3078

5.7528 −14.9391 −0.4586 5.8541 −12.5123

−2.3836 5.3517 −0.0465 −1.9928 4.5766

9.7209 −26.6163 0.1595 10.2078 −22.4336

−6.4265 16.2923 −0.1694 −6.5070 13.7647







B1=h−0.2857 0.2052 −0.1335 0.0756 0.5193 iT

C1=h0.0937 0.2055 −0.0000 0.1356 −0.7939 i,

D11 = [0].

The MATLAB robust control toolbox function

normhinf [9] reports that the closed loop H∞norm is

7.8540366769208774.

The MATLAB function hinfsyn computes γ= 7.8545

and constructs the controller accordingly.

Example 4.2 This is Example 3.1 and Example 6.2

from [6]. Consider the system







A B1B2

C1D11 D12

C2D21 0













−1 0 0 0 1

0−1 0 0 1

1 0 1

200

0 1 0 1

1 1 0 1 0







The variant γ-iteration described in [6] determined

γmo = ˆγ=.5000000000000 which agrees with the

theoretical value to thirteen significant digits. With

γ= 0.50001, our implementation returns the central

suboptimal controller (24) as

E=



−0.448419067323971 0.986059138740681

−0.047683732968496 0.379074901512008



×10−4,

A=



−0.097230000237421 0.087478934743184

−0.180766934371941 0.162712715362113



,

B1=



−0.193452124650871

−0.359503744256900



,

C1=h−0.2514183866414180.226354558822012 i,

D11 = [−0.5].

The function normhinf [9] reports that closed loop H∞

norm is 0.500009995.

The MATLAB function hinfsyn computes γ= 0.50625

and constructs the controller accordingly.

Example 4.3 This is Example 6.3 in [6]. Consider the

system in Example 4.2 with the (2,2) element of B1set

to one, i.e.,







A B1B2

C1D11 D12

C2D21 0













−1 0 0 0 1

0−1 0 1 1

1 0 1

20 0

0 1 0 1

1 1 0 1 0







The variant γ-iteration described in [6] reports γmo =

0.8062257748299. With γ= 0.80623, our implementa-

tion returns the central suboptimal controller (24) as

E=



−0.116338530296271 −0.308766348091874

0.187484278723965 −0.137335303455170



,

A=



0.188596106770156 0.745011728081016

−0.103219543410531 0.479206276745429



,

B1=



0.066440383622532

−0.376057188147526



,

C1=h−0.005120459018602 −0.213142135493424 i,

D11 = [−0.5].

The function normhinf [9] reports that closed loop H∞

norm is 0.80622598.

The MATLAB function hinfsyn computes γ= 0.80878

and constructs the controller accordingly.

Example 4.4 This is Example 6.4 in [6]. Let







A B1B2

C1D11 D12

C2D21 0













2 0 0 1 −1

0−1 0 1 −2

1 0 3 0 0

0 1 0 −1 1

4−2 0 1 0







In this example, ˆγ=γmo = 3. With γ= 3.0001, our im-

plementation returns the central suboptimal controller

(24) as

E=



0.205571073196335 0.034158098463962

−0.402143396055301 −0.066598384223367



,

A=



−0.616633677960107 −0.102994497996504

1.210370083410015 0.174898262642715



,

B1=



−0.235063370281313

0.459795474769660



,

C1=h1.752972656836075 0.265994049023116 i,

D11 = [1].

The function normhinf [9] reports that the closed loop

H∞norm is 3.00000006.

The MATLAB function hinfsyn computes γ= 3.00645

and constructs the controller accordingly.

Our last example is a standard textbook example.

Example 4.5 We consider a four-disk control system

taken from [22]. Due to space limitations, we refer to [22,

Example 19.4] for the data matrices. In [22] the optimal

H∞norm is given as γopt = 1.1272 and a controller is

computed for γ= 1.2. The MATLAB function hinfsyn

computes γ= 1.1292 and constructs the controller ac-

cordingly. Our experimental code computes γ= 1.1267

and is able to construct a controller for that value.

5 Conclusion

In this paper we have introduced new formulas for com-

puting optimal H∞controllers, that avoid all explicit

matrix inverses (except for some diagonal matrices, i.e.,

row or column scalings) and many potential numerical

difficulties that classical formulas may face. The formu-

las are based on the formulation used for the variant

version of the γ-iteration presented in [6]. Since this γ-

iteration has recently been extended to descriptor sys-

tems [16] we expect that these formulas can also be ex-

tended in a similar way. We have demonstrated the nu-

merical properties of the new formulas with several nu-

merical examples.