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]).

Loading more pages...