NASA ContractorReport178081

l_t.,SE REPORT NO. 86-15 NAS 19860015705 A-CR-178081

._J

ICASE

NUMERICAL APPROXIMATION FOR THE INFINITE-DIMENSIONAL

DISCRETE-TIME OPTIMAL LINEAR-QUADRATIC REGULATOR PROBLEM

J. S. Gibson

and

I. G. Rosen

Contract Nos. NASI-17070 and NASI-18107

March 1986

INSTITUTE FOR CO_UTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

LIB, ARYCOPY

MAY2 o !986 National Aeronautics and

SpaceAdministration LAI,IGLEY RESEARCH CENTF't_

Langley Research Gerlter LI_,_ARYN, ASA

Hampton, V_rginia 23665 HAMR]ON,VIRGI..AN'

NUMERICALAPPROXIMATIOFORN THEINFINITE-DIMENSIONADILSCRETE-TIME

OPTIMALLINEAR-QUADRATICREGULATORPROBLEM$

J.S. Gibson*

Department of Mechanical

Aerospace and Nuclear Engineering

University of California, Los Angeles

Los Angeles, CA 90024

and

**

I.G. Rosen

Department of Mathematics

University of Southern California

Los Angeles, CA 90089

ABSTRACT

An abstract approximation framework is developed for the finite and

infinite time horizon discrete-time linear-quadratlc regulator problem for

systems whose state dynamics are described by a linear semigroup of operators

on an infinite dimensional Hilbert space. The schemes included in the

framework yield finite dimensional approximations to the linear state feedback

gains which determine the optimal control law. Convergence arguments are

given. Examples involving hereditary and parabolic systems and the vibration

of a flexible beam are considered. Spline-based finite element schemes for

these classes of problems, together with numerical results, are presented and

discussed.

(*)This research was supported in part by the Air Force Office of Scientific

Research under contract No. AFOSR-84-0309.

(**)This research was supported in part by the Air Force Office of Scientific

Research under contract No. AFOSR-84-0393.

(?)Part of this research was carried out while the authors were visitin_

scientists at the Institute for Computer Applications in Science and

Engineering (ICASE), NASA, Langley Research Center, Hampton, VA, 23665 which

is operated under NASA contracts No. NASI-17070 and No. NASI-18107.

i

I. Introduction

Recent advances in micro-processor technology have led to increased

interest in digital or discrete-tlme control systems. In addition, because

many current application areas involve complex systems which are most

appropriately modelled using functional and/or partial differential equations,

it has become important to study digital control techniques in the context of

infinite dimensional or distributed systems.

A great deal of attention has been given to the contlnuous-time infinitedimensional llnear-quadratlc regulator problem. The general theory and

characterization of the linear state feedback form of the optimal control are

discussed in [5], [6], [8], [9], [21] and [22], while its application to

hereditary, parabolic and hyperbolic systems with emphasis on approximation is

treated in [2], [3], [7], [I0], [II], [14] and [17] to mention just some of

the work that has been done.

On the other hand, relatively little can be found in the literature

concerning the corresponding dlscrete-tlme problem. The major contributions

in this area can be found in the papers by Lee, Chow and Barr [20] and Zabczyk

[28]. In these studies the Riccatl difference equations that characterize the

linear feedback form of the optimal control for the finite time problem are

given and limiting properties as the length of the time horizon tends to infinity

are discussed. However, the issue of approximation is not considered.

In the present paper, we develop numerical approximation schemes that

yield finite dimensional approximations to the feedback gain operators which

determine the discrete-time optimal control law. We consider control systems

whose dynamics can be described in terms of a linear semigroup of operators on

an infinite dimensional Hilbert space. The basis of our approach is the

construction of a sequence of finite dimensional (presumably finite element
-2-

based) state approximations which in turn leads to a sequence of finite

dimensional discrete-tlme llnear-quadratic regulator problems each of which

can be solved using standard techniques.

Under appropriate assumptions on the nature of the original problem and

the convergence of the state approximation, we are able to prove that the

approximating optimal controls and feedback gains converge to the true optimal

control sequences and feedback laws for the original infinite dimensional

system. Depending upon the convergence properties of the state approximation,

we are able to establish strong or uniform norm convergence of the

i approximating gain operators and the corresponding weak or strong convergence

of the approximating feedback kernals which are used in the implementation of

the optimal control. We treat both the finite and inflnlte-tlme horizon

problems.

We have tested our schemes on a wide variety of examples. This paper

includes numerical results for problems with state dynamics given by

hereditary and parabolic (heat/dlffuslon) differential equations and a hybrid

system of partial and ordinary differential equations for the vibration of an

Euler-Bernoulll beam connected to a rigid body and a lumped mass. We

implemented and tested the methods on an IBM Personal Computer.

We give a brief outline of the remainder of the paper. In section 2 we

breifly outline previous results concerning the characterization of the

optimal control and feedback gains for both the finite and infinite time

horizon dlscrete-tlme regulator problem for distributed systems. The Riccatl

difference and algebraic equations whose solutions determine the optimal

feedback control law are discussed. In section 3 we develop the abstract

approximation framework and convergence arguments. Section 4 contains a

discussion of particular schemes for the classes of problems mentioned above
-3-

together with the results of our numerical studies. Some concluding remarks

are given in Section 5.

We employ standard notation throughout. For an interval (a,b), we denote

by Hk(a,b) the usual Sobolev spaces of real-valued functions defined on (a,b)

whose (k-l)st derivatives are absolutely continuous and whose kth derivatives

are L 2. The standard Sobolev inner product on Hk(a,b) is denoted

by <'">k" For X and Y normed linear spaces we denote by I(X,Y) the space of

bounded linear operators from X into Y. When Y = X, we use the shorthand

notation i(X).

2. The Optimal Control Problem

2.1 Optimal Control on a Finite Interval

Let Z and U be Hilbert spaces with inner products <'">Z and <'">U'

respectively, with U finite dimensional. For {H, <.,.>H } a Hilbert space, let

£2(t 0, x(t)}_t0

tf;H) denote the usual Hilbert space of sequences x = { _ with

x(t) _ H together with the inner product

tf

(2.1) <x,Y>£2 = _ <x(t), y(t)> H.

t=t 0

The discrete-time linear quadratic regulator problem on the finite time

interval [t0,tf] is

(PI) Choose u € £2(t0,tf;U ) to minimize the quadratic performance index
-4-

J(G;to,tf,z(t0),u) =

(2.2) tf-I

[<Qz(t), z(t)> Z + <Ru(t),u(t)>u] + <Gz(tf),z(tf)> Z

t=t 0

subject to the discrete-time control system

z(t+l) = Tz(t) + Bu(t), t _ to

(2.3)

z(t0) € Z,

where T and B are bounded linear operators from Z into Z and U into Z,

respectively, Q and G are bounded, nonnegative self-adjoint operators on Z,

and R is a positive definite self-adjoint operator on U.

Of primary concern to us will be applications where (2.3) is the sampled

form of the continuous-time control system

(2.4) z(s) =nz(s) + Bu(s)

where n is the infinitesimal generator of a C0-semigrou p of bounded linear

operators T(s), s _ O, on Z, and B is a possibly unbounded linear operator

from U into Z. In this case we have

(2.5) T = T(T) and B = f0 T(s)Bds,

where T is the sampling interval. If, as in our subsequent example

discussed in Section 4.1 where u is a boundary control in a heat equation, B

is unbounded (more precisely, B maps U not into Z but into some larger space),

then the integral in (2.5) is not interpreted literally.
-5-

The solution to Problem (PI) has been given for infinite dimensional

control systems in [20],[28], and the equations representing the solution have

the same form as in the finite dimensional case. We will give now the version

of the solution that is most useful for our purposes.

For given z(t0), J(G;t0,tf,z(t0),u) is a bounded linear-quadratic

functional on %2(t0,tf;U) with coercive quadratic part. Therefore,

for each

z(t0) , there exists a unique optimal control sequence in £2(t0,tf;U). Also,

the minimum value of the _erformance index is a quadratic functional of z(t0) ,

so that there exists a unique nonnegatlve, self-adjolnt H(t 0) E i(Z) such

that

(2.6) J, = min J(G;t0,tf,z(t0),u) = <]l(t0)z(t0),z(t0)> z.

Application of the principle of dynamic optimality establishes that the

optimal control has the feedback form

(2.7) u,(t) = -F(t)z,(t), to • t • tf-I

where

(2.8) F(t) = R(t)-IB*_(t+I)T,

^

(2.9) R(t) = R + B H(t+I)B

and H(t) satisfies the Riccati difference equation

(2.10) H(t) = T [H(t+l) - H(t+I)BR(t)-IB*E(t+I)]T + Q, t • tf-l,
-6-

with the final condition

(2.11) n(tf) = G.

The optimal trajectory z, is given by

(2.12) z,(t+l) = S(t)z,(t), t ) to

where

(2.13) S(t) = T-BF(t).

2.2 Control on the Infinite Interval

Here, tf = = and G = 0. To simplify notation, we will write

J(t0,_,z(t0),u) instead of J(0;t0,=,z(t0),u).

Definition 2.1. A control sequence u _ £2(0,=;U) is an admissible control

for the initial condition z if J(0,=,z,u) < = .

The discrete-time linear-quadratic regulator problem on the infinite

interval is

(P2) Choose an admissible control u, to minimize J(0,=,z,u), if an

admissible control exists for the initial condition z.

That a unique optimal control u, exists whenever at least one admissible

control exists follows from the fact that the quadratic part of J(0,=,z,u) is

coercive on a subspace of £2(0,_;U). See the discussion following Definition

4.1 of [9].
-7-

Definition 2.2. A bounded linear operator _ on Z is a solution to the

Riccati algebraic equation if

* * -I *

(2.14) _ = T [_-HB(R+B _B) B H]T + O.

The following theorem summarizes results from Zabczyk [28].

Theorem 2.3. The following are equivalent:

(i) There exists an admissible control for each z s Z;

(ll) for each z s Z, sup _ql(t)z,z>Z < _, where H(t) is the Riccatl

t<tf

operator in (2.10) and H(tf) = 0 for fixed tf;

(lii) as t . -_, H(t) converges strongly to a nonnegatlve self-adjoint

solution to the Riccatl algebraic equation;

(iv) there exists a nonnegative self-adjolnt solution to the Riccati

algebraic equation.

For uniqueness of the solution to the Riccatl algebraic equation and

characterization of the optimal control, Zabczyk treated two cases: when Q is

coercive, and when the spectral radius of T is less than 1 (i.e., the openloop system is uniformly exponentially stable). Since neither is the case in

the example we discuss in Section 4.2 and other applications in which we are

interested, we will need the following hypothesis and theorem.

Rypothesls 2.4. The operators T, B and Q are such that, if z(0) s Z and u is

an admissible control for z(0), then

(2.15) lim Iz(t)Iz = 0.

t+_
-8-

Theorem 2.5. When Hypothesis 2.4 holds, there exists at most one nonnegative

self-adjoint solution to the Riccati algebraic equation. If such a solution

exists, then there exists a unique solution to problem (P2) for each

initial condition z(0) _ Z, the minimum value of the performance index is

(2.16) J* = u acmlssiDie'mln .,i J(0,=,z(0),u) = <_z(0),z(0)>Z,

the optimal control has the feedback form

(2.17) u,(t) = -Fz,(t), t ) 0,

where

(2.18) F = R-IB*_ T,

(2.19) R = R+B E B

and the optimal trajectory z,(t) satisfies

(2.20) z,(t+l) = Sz,(t), t ) 0,

with

(2.21) S = T-BF.

Proof. Let E be such a solution and note that, for any finite tf, E is a

constant solution to (2.10) and (2.11) with G = E. Then the corresponding
-9-

F(t) and R(t) defined by (2.8) and (2.9) are the constant operators in (2.18)

and (2.19). For z(0) g Z, define _(0) = z(0),

(2.22) z(t+l) = (T-BF)z(t), t • O,

and

(2.23) _(t) = -F_(t), t • 0.

Now suppose that u is an admissible control for z(0) and that z(t) is the

corresponding solution to (2.3). For tf > 0, the preceding results about the

solution to Problem (PI) with G = _ imply

J(_;0,tf,z(0),u) _ J(0;0,tf,z(0),u) + _qlz(tf),z(tf)> Z

(2.24)

J(0;0,_,z(0),u) + _Iz(tf),z(tf)> Z.

Also,

J(H;O,tf,z(0),u) = _qlz(0),z(0)> Z

(2.25)

= J(O;O,tf,z(O),u) + _qlz(tf),z(tf)> Z.

Since z(tf) + 0 as tf + _, (2.24) shows that u is both admissible and

optimal for Problem (P2). Since _(tf) . 0 as tf + _, (2.25) shows (2.16).

As we see now, (2.16) must hold for any nonnegative self-adjoint solution of

the Riccatl algebraic equation; therefore, such a solution is unique.

Remark 2.6. When Hypothesis 2.4 does not hold, the Riccati algebraic equation

may have more than one nonnegative self-adjoint solution. In this case, the

minimal such solution -- there will be one -- gives the solution to Problem
-I0-

(P2) as in Theorem 2.5. Throughout this paper, we assume that Hypothesis 2.4

holds.

Lemma 2.7. Suppose that Q > m for some positive constant m, and set

n

* tQTt Cn = Z (T) , for n = 1,2, ... . Then ICnZlz is bounded in n for each

t=0

z _ Z if and only if Cn converges in norm to the operator

(2.26) C = l (T*)tQT t

t=O

and

(2•27) ITtl< (JCllm)(1- mllCl)t, t = 1,2, .e• •

Proof, Since Cn is an increasing sequence of bounded self-adjoint linear

operators, Cn converges strongly to some bounded self-adjoint C if and only if

<CnZ,Z> Z is bounded in n for each z, if and only if ICnZlz is bounded in n for

each z. This is a standard result• The proof of the Lemma is then a standard

exercise using the Lyapunov functional <Cz(t), z(t)> Z for the homogeneous part

of (2.3).

Corollary 2.8. If Q _ m > 0 and the Riccati algebraic equation has a

nonnegative self-adjoint solution H, then the spectral radius of the operator

S in (2.21) is less than I, and

(2.28) JstJ_ (InJ/m)(1-mlnJ)t, t = 1,2,•ee •

Proof. This follows from Lemma 2.7 and
-II-

*)t[Q + F RF]S t. (2.29) <]Iz,z>Z = 7 (S t=O

For Q coercive, Zabczyk proved a stronger result than part (iii) of

Theorem 2.3: if a nonnegative self-adjoint solution to the Riccati algebraic

equation exists, then l_(t) - H I . 0 geometrically fast as t . -_ (Also, see

[13]). We will need such a result, along with an explicit convergence rate,

for the approximation theory in Section 3.2. Since Zabczyk's proof does not

yield an explicit convergence rate, we give the following.

Theorem 2.9. Suppose that there exists a nonnegative self-adjoint solution

to (2.14) and that

(2.30) Istl ( Mr t, t = 1,2, ... ,

where M and r are positive constants with r < 1 and S is the optimal closedloop operator in Theorem 2.5. If H(.) is the operator in (2.10) with tf = 0

and

(2.31) H(O) • _,

then

(2.32) <Hz,z> Z < <]l(-t)z,z> z < <_z,z> Z + (Mrt)21H(O)l, t = 1,2, ....

Proof. For to a negative integer, let u 0 be the optimal control sequence for

the finite-time Problem _I) on the interval [t0,0] with initial condition
-12-

z(t0) _ Z, with z0 the corresponding optimal trajectory. Also, let u, be the

optimal control sequence on the infinite interval for Problem (P2) with

initial condition z(t0) , with z, the corresponding optimal trajectory.

Since H is a constant solution to (2.10) for the final condition G =H, we

have

<]Iz(to),Z(to)> Z = J(H;to,O,z(to),U,)

(2.33) (J(O;to,O,z(to),Uo) + <Hzo(O),zo(O)> Z

J(O;to,O,z(to),Uo) + <H(O)zo(O),zo(O)> Z

= _l(t0)z(t0),z(t0)> Z.

-t

On the other hand (note that z,(t0) = S z(t0)),

<_(to)Z(to),Z(t0)> Z

(2.34) _ J(0;0'-t0'z(t0),u,) + _ql(O)z,(-t0),z,(-t0)> Z

J(0;O,_,z(t0),u,) + lql(0)z,(t0),z,(to)> Z

-to

_Iz(t0),z(t0)> Z + IH(0)I(IS llZ(to)Iz )2
-13-

3. Approximation Theory

3.1 The finite time interval problem

In this section we develop a general approximation framework for the

finite time interval problem (PI) and describe associated convergence results.

For each N = 1,2, ... , let ZN c Z be a finite dimensional subspace of Z

and let PN: Z . ZN denote the orthogonal projection of Z onto ZN with respect

to the <'">Z inner product. We require the following hypotheses.

Hypothesis 3.1 There exist operators TN: ZN + ZN, BN: U + ZN

QN: ZN . ZN and GN: ZN . ZN which satisfy

TNP N . T strongly,

TN PN . T strongly,

BN . B strongly,

QNPN + Q strongly,

GNP N . G strongly,

as N . = with TN and BN bounded and QN and GN bounded, self-adjoint and

nonnegative.

Hypothesis 3.2 The spaces ZN are approximating subspaces in the sense that

the projections PN satisfy PN . I strongly on Z as N . =.

We note that since U has been assumed to be finite dimensional,

Hypothesis 3.1 above necessarily implies that BN . B and BNP N . B in the

uniform norm topology on I(U,Z) and I(Z,U) respectively.
-14-

We definea sequenceof approximatingdlscrete-tlmelinearquadratic

regulatorproblemson the finitetimeinterval[t0,tf] as follows:

(P1N)

u_ £2(t0,tf-l;U) which minimizes Find

tf-I

(3.1) JN(GN;t0,tf,z(t0),u)= [ [<QNZN(t),ZN(t)>z +

t=t0

<Ru(t),u(t)> U] + <GNZN(tf), ZN(tf)> z

subject to

(3.2) ZN(t+l) = TNZN(t) + BNU(t), t > to

ZN(t 0) = PNZ(t0).

The results stated in Section2.1 concerningthe existenceand uniquenessof

solutionsto Problem (PI) apply to the Problems (PIN) as well. Indeed, there

exists a unique solutionu,N € £2(t0,tf-l;U) to Problem (PIN)which is given

in feedback form by

(3.3) u,(t) N = -FN(t ) z_(t) , to _ t tf-I

where

^ -I *

(3.4) FN(t ) = RN(t ) BN EN(t+I)T N
-15-

with

^ ,

(3.5) RN(t) = R + BN HN(t+I)B N

tf

and the operators {_N(t)}t=t0 on ZN satisfying the Riccati difference

equation

* ^ -I *

(3.6) HN(t) = TN[gN(t+I) - _N(t+I)B_N (t) B_N(t+I)]TN + QN

with terminal condition

(3.7) HN(t f) = GN.

N

The optimal trajectory z, is given by

N

z,(t+1) = SN(t)z_(t), t _ to, (3.8)

N

z, (to ) = eNz(to )

where

(3.9) SN(t) = TN - BNFN(t) , t • to •

tf

The operators {EN(t)}t=to are bounded, self-adjoint and nonnegative. The

minimum value of the performance index (3.1) is given by

= N N

(3.10) jN JN(GN;to,tf,Z(to),U N) = _qlN(to)Z,(to) ' z,(to)>z

The fundamental convergence result is given in the following theorem.
-16-

N

Theorem 3.3 Let u, and u, be the unique solutions to problems (PIN) and (PI),

N

respectively, with z, and z, the corresponding optimal trajectories

generated by (3.8) and (2.12). Let JN' HN and FN and J, H and F be given by

(3.1), (3.6) and (3.4) and (2.2), (2.10) and (2.8). Then, if Hypotheses 3.1

and 3.2 hold, we have

N

(i) N+°_limlu, - u, 1£2 = 0,

N

(ii) N+oolim Iz.- z,I 2 = 0,

(iii) lim IJN - J*l = O,

N+_

(iv) lira IHN(t)PNZ -n(t)Zlz = 0, z _ z, to < t _ tf N+_o

and

(v) lim IFN(t)P N - F(t)l = 0, to < t < tf -I. N+_o

Proof

HN(t) being nonnegative implies that l&(t)l > IRl

We first note that

consequently that IRN(t)-II _ IRl-I. It follows therefore that for

and

u € U.

(3.11) I(RN(t)-I - ^R(t)-l)ulu

= l_(t)-l(R(t) - _(t))R(t)-lulu
-17-

^ ^ i

IRI-II(R(t) - _(t))R(t)- Ulu.

The above estimate together with (2.9), (2.11), (3.5), (3.7) and Hypothesis

3.1 imply that

(3.12) R N(tf_l)-I . ^R(tf-l) -I

as N . = strongly on U. Since U is finite dimensional the convergence in

(3.13) is in fact uniform. It then follows immediately from (2.8), (3.4) and

Hypothesis 3.1 that

(3.13) FN(t f - I) PN . F(tf - I),

uniformly as N + =, and from (2.10) and (3.6) that

(3.14) HN(t f - I)PN + H(tf-l)

strongly on Z as N . =. A simple induction yields (iv) and (v) from which

(i), (ii) and (iii) then follow trivially.

Remark It will , on occassion, be the case that in constructing a particular

approximation scheme TNP N . T strongly but TNP N + T only weakly (see, for

example, [3]). However, by using the fact that

(3.15) (T_N(t+I)) = HN(t+I)T N
-18-

implies that T_N(t+I ) . T H (t+l) weakly if HN(t+I) . H(t+l) weakly, we

conclude that Theorem 3.3 continues to hold under these somewhat weaker

hypotheses with the strong convergence in (iv) replaced by weak and the

uniform convergence in (v) replaced by strong.

Under certain additional hypotheses it can be shown that the operators

H(t), to _ t _ tf given by (2.10), (2.11) are trace class (see [15]) and that

(3.16) lim U_N(t)P N - H(t)D 1 = 0, to _ t _ tf, N+_

where N.Xl1 denotes the trace norm, the strongest of all common operator

norms. We require the following lemmas.

Lemma 3.4 If {ai}i= 1 is an absolutely summable sequence of real numbers then

there exist sequences {bi}i= 1 and {ci} =I such that limb i = 0, {ci} =I is i+_

absolutely summable and ai --bic i.

Proof

Let

(3.17) e = y lail 1=1

and for j = 0,1,2,... define nonnegative integers nj as follows. Let nO = 0

and let nj denote the first index for which

n°

3 1

I latl >= - i=l 3

j = 1,2, .... Set
-19-

1

(3.19) bi = _ , ci = Jai, i = nj_l+l,...,nj, j = 1,2,...

Then bic i = ai, i = 1,2,..., limb i = 0 and

i+_

n

GO OO j OO

(3.20) I Ici[ = I J I [ak] _ _ + I 7f< i _-

i=l j--1 k=nj_1+l j=l 3

Lemma 3.5 If L is a self-adjoint trace class operator on a separable Hilbert

space H, then L can be written as LIL2 where L 1 is compact and L2 is trace

class.

Proof

Let {li}_= 1 denote the eigenvalues of L repeated according to

multiplicity and let {@i}_=l denote the corresponding eigenvectors. Then

{li}_= I is a sequence of real numbers, each of finite multiplicity, and

oo

(3.21) I Ill I =]L] 1 <" •

i=l

Applying the oprevious o lemma there exist sequences {Bi i=l and {vi} =I with

lim _i = 0, Z {vi[< _ and I i = Bivi . Defining L 1 and L 2 by

i+= i=l

(3.22) LI@ = [ Bi<_,_i>H#i, # € H

i=l

and

GO

(3.23) L2_ = [ vi<_'@i>_i' # g H

i=l
-20-

respectively, the lemma immediately follows.

Lemma 3.6 Let {SN}N= 1 be a sequence of bounded linear operators on a

seperable Hilbert space H which converges strongly to a bounded linear

operator S. Let {LN}_= I be a sequence of trace class operators on H which

converges in trace norm to an operator L. If L can be written as L = L1L 2

with L I compact and L 2 trace class then the sequence {SNLN}N= I converges in

trace norm to SL.

Proof

The result follows immediately from

(3.24) iiSNLN - SLIt1 _ gSN(L N - L);I1 + II(SN - S)LIL2111

[SNI'ILN - LI;1 + i(SN - S)LINIIL2i,I.

Theorem 3.7 If Q and G are trace class operators then the operators

tf

{H(t)}t=t0 given by (2.10) and (2.11) are trace class. Moreover, if

Hypotheses 3.1 and 3.2 hold and QNPN . Q and GNP N . G in trace norm as

N . = then we have

(3.25) lim IIEN(t)PN - H(t)ll1 = 0, to _ t _ tf. N+_

Proof

That the operators H(t), to ( t _ tf are trace class is an immediate

consequence of the hypotheses of the theorem, (2.10), (2.11) and the fact that
-21-

the trace class operators form a two sided ideal of [(Z), the space of bounded

linear operators on Z (see [15]).

The trace norm convergence stated in (3.25) will follow once we have

shown that

(3.26) llm IIHN(t+I)PN - H(t+l)llI = 0

N+_

implies

(i) lim IIT_N(t+I)TNP N - T H(t+l)Tfl I = 0

N+_

and

* ^ -I *

(ii) llm liTN HN(t+I) B_N(t) B_N(t+I)TNPN -

N+_

T _ (t+I)BR(t)-IB*_ (t+l)Tlli 0.

To argue (1) we first note that Hypothesis 3.1 and Lemmas 3.5 and 3.6

imply

(3.27) lim lITN EN(t+I)PN - T H(t+l)ll1 = 0.

N+_

Taking adjoints we obtain

(3.28) llm II_N(t+I)TNPN -E(t+I)TU 1 = 0.

N+_

Another application of the previous two lemmas yields

(3.29) lim liTN _N(t+I)TNPN- T H(t+I)TIII< lim ITNiflgN(t+l)TNP N - g(t+l)TflI

N+= N+_
-22-

* 1

+ lira I(TNPN - T ) _ (t+l)lllH2(t+l)T1ll N+_

where _(t+l) = Hl(t+l)_2(t+l) is the factorization of H(t+l) described in

Lemma 3.6.

The verification of (ii) is analogous and the theorem is proven.

We note that if Hypotheses 3.1 and 3.2 hold and if the operators Q and G

are trace class with QN and GN defined by

(3.30) QN = PNQ

and

(3.31) GN = PNG,

then Lemmas 3.5 and 3.6 imply that the trace norm convergence hypotheses in

Theorem 3.7 hold. The significance of this observation will become apparent

when examples are discussed in Section 4. Indeed, it is frequently the case

in practice that Q and G have finite rank (and consequently are trace class)

and the operators QN and GN are defined as in (3.30), and (3.31).

3.2 Approximation on the Infinite Interval

Problem (p2N) is Problem (P2) for the control system in (3.2) and the

performance index

(3.32) JN(0,_,ZN(0),u) = _ [<QNZN(t),zN(t)> Z + <Ru(t),u(t)>u]" t=O
-23-

Hypothesis 3.8. For each N, there exists exactly one nonnegative self-adjolnt

solution to the Riccati algebraic equation

* * -1 *

(3.33) HN = TN[HN - H_N(R + B_N) B_N]TN + QN"

By Theorem 2.3, this implies that

(3.34) lim [HN - HN(t) [ = 0

t+ -_o

for each N, since dlm(Z N) < _ •

As in Theorem 2.5, we write

(3.35) FN = RN I BN HNTN,

(3.36) RN = R + B_N,

and

(3.37) SN = TN - BNFN.

From here on, H will be the nonnegative self-adjoint solution to the infinite

dimensional Riccati algebraic equation (2.14) -- when it exists -- F will be

the corresponding feedback operator in (2.18) and S will be the corresponding

closed-loop operator in (2.21).
-24-

Theorem 3.9. If HNP N converges strongly to some bounded linear operator H,

then H is a nonnegatlve self-adjoint solution to (2.14), FNPN converges in

norm to F and SNP N converges strongly to S.

Proof. This follows from Hypotheses 3.1 and 3.2, (3.33) and (3.35) - (3.37),

and the fact that the control space U has fixed finite dimension.

Theorem 3.10. Suppose that there exist positiveconstantsM and r,

independentof N, with r < I, such that

(3.38) HN ( M, N = 1,2, ... ,

and

(3.39) ISNI ( Mr t, t = 1,2, ... , N = 1,2, ....

Then a nonnegatlve self-adjolnt solution H to (2.14) exists, and as

N+_,

(3.40) HNP N . H strongly.

If there exists a positive m, independent of N, such that

(3.41) QN _ m, N = 1,2, ...,

then (3.38) implies the existence of an r less that one and independent of N

for which (3.39) holds.
-25-

Proof. For each N, let nN(') satisfy (3.6) with tf = 0 and nN(0) = MI, where

I denotes the identity operator on ZN. From (2.32),

(3.42) InN - nN(-t) ] + 0 as t + _,

uniformly in N. Now, for z E Z, write

(3.43) <(n N - nN,)Z,Z> Z = <(nN - nN(-t))z,z> Z + <(nN(-t) - nN,(-t))z,z> Z

+ <(_N,(-t) - nN,)Z,Z> Z.

For E > 0 choose t > 0 such that I(HN - nN(-t))Zlz < _ and

l(nN, - _N,(-t))zl < _. Then, for N and N' large enough,

l(nN(-t ) - HN,(-t))Zlz < €. This shows that HNZ is a Cauchy sequence in Z

for each z. Therefore, n N converges strongly to a nonnegative self-adjoint

solution to (2.14).

An important application of this theorem is when the approximating openloop operators TN have an exponential decay rate independent of N, Q is

coercive and QN = PNQIz N" In this case, the zero control gives an upper

bound, independent of N, on HN. Such is the case in the example discussed in

Section 4.1 and in applications to flexible structures with no rigid-body

modes and coercive structural damping.

Theorem 3.11. Suppose that HNP N converges strongly to n, QNPN converges in

trace norm to Q (hence Q is trace class), and (3.39) holds for positive M and

r independent of N with r less than one. Then nNP N converges in trace norm

to _.
-26-

Since IISNQNSNIII( 'ISNI211QNII,1 the series in (2.26) converges in

Proof.

trace norm, uniformly in N. The current result follows then from Lemmas 3.5

and 3.6.

Note that QNPN converges in trace norm to Q if Q is trace class and

Q = PNQPNIz NTheorem 3.12 if IHNI is bounded in N, then a nonnegative self-adjolnt

solution H to (2.14) exists, HNP N converges weakly to E, and FNP N and SNPN

converge strongly to F and S, respectively.

Proof. According to [II Theorem 6], HNP N converges weakly to some

nonnegatlve self-adjoint bounded H. It follows from (3.33) and Hypotheses 3.1

and 3.2 that H satisfies (2.14) and that FN and SN converge as indicated.

Note that Theorem 3.12 holds if SNPN converges strongly but SNPN

converges only weakly.

3.3 Implementation of the Approxlmatlon Schemes

In constructing the approximating operators TN, BN,QN and GN a standard

Galerkln approach is often taken; that is, TN = PNT, BN = PNB,

QN = PNQ and GN = PNG. We note however that explicit representations for the

operators T and B are frequently not available. In particular, this can occur

when the dlscrete-time system (2.3) arises from the sampling of an infinite

dimensional contlnuous-tlme system of the form (2.4). In this case it is the

operators A and B which are approximated by a sequence of finite dimensional
-27-

operators A N and BN on ZN, from which an approximation to the semigroup

[ T(s) : s _ 0} is obtained as TN(S) = exp (ANS) , s _ 0. The operators TN

and BN are then TN= T N (T) and BN = f_ TN(S)BNdS , respectively. The strong

convergence TNP N . T and BN . B is then usually argued using an appropriate

formulation of the Trotter-Kato theorem, a well known semigroup approximation

result (see [15] [23]).

The expressions given by (3.3) - (3.7) are operator equations and

although they are finite dimensional, they are not appropriate for

computations. To make use of our approximation framework, we must first

determine equivalent matrix formulations. Toward this end we assume, without

loss of generality, that U = Rm with the standard basis and inner product and

_ i_KN

let i_N}i= 1 be a basis for ZN. Define the _ × KN Gram matrix M N by

i j

(3.44) [MN]ij = <#N,_N>Z .

For an operator A we denote its matrix representation with respect to the

bases defined above by [A]. Similarly, for an element z m Z or u _ U, we let

its vector representation be given by [z] or [u] respectively. Standard

calculations yield

n

(3.45) [TN] = MN I[TN]TM N

and

(3.46) [BN] = [BN]TMN .

Defining
-28-

A

(3.47) _N(t) = _[_N(t)],

A

(3.48) QN = MN [QN]'

and

A

(3.49) GN = MN[GN]

we obtain

N N

(3.50) [u,(t)] = -[FN(t)][z,(t)], to • t • tf -I,

A

(3.51) [FN(t)] = [RN(t)] -I[BN]T_N(t+I)[TN] ,

(3.52) [RN(t)] = [R] + [BN]THN(t+I)[BN ],

(3.53) _N(t)- [TN]T(_N(t+I ) -

A A

_N(t+I)[BN ][RN(t)]-I[BN]T_N(t+I))[TN ] + QN' to • t • tf-l,

A A

(3.54) HN(tf) = G N.

A

Note that since QN and GN are self-adjoint and nonnegative so too are QN

A

and GN. Equations (3.50) - (3.54) are therefore in the form of the standard

ones obtained for the feedback law for a discrete-time linear-quadratic

KN

regulator problem in R . Consequently they can be solved using conventional

techniques. The minimum value of the performance index is given by
-29-

(3.55) jN N T_ N , = [z, (to)] N(t0) [z,(t 0) ]•

Analogously, for the infinite time horizon problem, (3.33), (3.35) and

(3.36) yield

N N

(3.56) [u,(t)] = -[FN][Z,(t)], t • to,

(3.57) [FN] = [RN]-I[BN]T_N[TN ]

(3.58) [RN ] = [R] + [BN]T_N[BN ]

A

where HN is the solution to the matrix algebraic Riccati equation

(3.59) HN = [TN]T(_N - _N[RN]-I[BN]T_N)[TN ] + QN'

with QN given by (3.48). The minimum value of the performance index is given

by

* N ]T_ N

(3.60) JN = [z*(t0) N [z,(to)].
-30-

4. Examples and Numerical Results

In this section we describe the application of the general a_proximation

framework developed above to a variety of examples. In addition to

theoretical considerations, in each of the examples below, we discuss some

numerical results for an inflnite-time horizon problem of the form given in

Problem (P2). All numerical studies were performed on an IBM Personal

Computer. The machine we used was equipped with an Intel 8086 math coprocessor chip and 640K bytes of random access memory (of which less than 384K

was required).

Matrix exponentials were computed from elgenvalue-eigenvector

decompositions obtained using the QR algorithm. The matrix Riccati equations

(3.59) were solved using a Schur-vector decomposition of the Hamiltonlan

matrix (see [18][24]). It should be noted that if the elgenvalue pairs of the

Hamiltonlan matrix for a continuous-time linear-quadratic regulator problem

are asymptotic to _+Y(n) as n . _, then the eigenvalue pairs of the

Hamiltonlan matrix for the corresponding discrete-time problem will be

±y(n)T

asymptotic to e as n . _. Consequently, for all but very small T,

conditioning problems arise more quickly than in the continuous-time case when

the approximating matrix algebraic Riecati equations are solved.

4.1. The Heat Equation with Boundary Input

In this example we consider the scalar parabolic system with boundary

control given by

(4.1) _s w(s,x) = _x a(x) _x w(s,x), s > 0, x _ (0,I)
-31-

(4.2) w(s,0) = 0, w(s,l) = v(s)

(4.3) w(0,x) = _(x)

with a g HI(0,1), a(x) > 0, x _ [0,I], _ g H0(0,1) = L2(0,1) and v g L2(0,=).

To formulate the discrete-time state equation for this system we let T

denote the sampling interval and consider only piecewise constant controls v

given by

(4.4) v(s) = u(t) s g [tT, (t+l)T),

t = 0, 1,2, .... We choose as our state space Z the Sobolev space

H0(0,1) with the usual inner product

(4.5) <€'€>Z = <€'€>0= .f_,(e)_(e)de.

The state z(t) _ Z is

(4.6) z(t) = lim _w(s,.), t = 1,2, ...

s+tT

(4.7) z(0) = #.

For t _ {0,1,2...}, we define y(s) g Z by

(4.8) y(s) = w(s,') - _0u(t), s € (tT, (t+l)T)
-32-

(4.9) y(tT) = z(t) - _0u(t),

where _0 £ Z is given by _0(x) = x, x _ [0,I]. A straight forward

calculation reveals that y(s) = y(s,.) satisfies

(4.10) y(s) = DaDy(s) + a'u(t), s g (tT, (t+l)T)

(4.11) y(s)10 = 0, y(s)[1 = 0, s g (tT,(t+l)T)

(4..12) y(tT) = z(t) - _0u(t),

where D denotes the differentiation operator on HI(0,1).

Let A: dom(A) c Z + Z be given by

H2(0, I) Dora(A)

n H_(0,1) = {0 v H0(0,1): € H2(0,1),

€(o) : €(1) : o}

(4.13) A_ = DaD_ .

The operator A is densely defined and self-adjoint. It satisfies

2

(4.14) <Az,z>Z < -_IZ[z , z _ Dom(A)

for some m > 0 and has compact resolvent. Also, A is the infinitesimal

generator of an analytic semigroup of contractions {T(s) : s _ O} on Z which,

-_0s

in light of (4.14), satisfies IT(s)[ _ e , s _ O. It follows therefore,

that
-33-

(4.15) y(s) = T(s-tT)y(tT) + fS T(s-o)a'dau(t), s _ [tT (t+l)T). tT

The continuity of y, (4.6), (4.8) and (4.9) imply

(4.16) z(t) = y(tT) + _0u(t)

and

(4.17) z(t+l) = y((t+l)T) + _0u(t),

and hence that

(4.18) z(t+l) = y((t+l)T) + _0u(t)

= T(T)(z(t) - _0u(t)) + -r(t+l)_tTT((t+l)T-o)a'dou(t) + _0u(t).

I

Defining the operators T _ /(Z) and B E /(R ,Z) by

(4.19) Tz = T(T)z, z _ Z

and

T T(o)a'do]u, u € R I

(4.20) Bu = [(I - T(T))_ 0 + f0

we obtain

(4.21) z(t+l) = Tz(t) + Bu(t), t = 0,1,2,..
(4.22) z(O) = </> •

-34-

We take the performance index to be '

(4.23)

with qo,go > 0 and r > O.

Applying the theory developed in Section 2.1, we have, for the finite

time interval problem, that the optimal control is given by

(4.24) t = 0, 1,2,••• ,tf

-1,

where for each t, F(t) is the continuous linear functional on Z given by (2.8)

o - (2.11). It follows that F(t) has a representation f(t,·) € H (0,1) and

that

(4.25) 1 F(t)1/I = JOf(t,' )1/1(8 )d8 ,

For the infinite time interval problem (tf = =, ~O = 0), it is

immediately clear that Hypothesis 2.4 is satisfied. It is also clear that

(4.14) imp1es that for each z(O) €

0H (0,1), u(t) = 0, t = 0,1,2, ••• is an

admissible control, and hence that there exists a unique nonnegative se1fadjoint solution of the Riccati algebraic equation (2.14). From (2.17) -

(2.19) we obtain
-35-

(4.26) u,(t) = -Fz,(t), t = 0,1,2,...

where F is a continuous linear functional on Z and

(4.27) F_ = f_f(0)_(0)d0, _ _ H0(0,1),

with f _ H0(0,1).

We define an approximation scheme using a standard Ritz-Galerkin

approach. We note that the operator -nis coercive and that it can be written

as -A= L L where the operator L and its adjoint L are given by

(4.28) L_ = al/2D_, @ _ H_(0,1),

and

(4.29) L*_ = _Dal/2_ , _ E HI(0 ,I),

1

respectively. We define the space V = H0(0,1) together with the inner

product

(4.30) <_'_>V = <L_'L_>z' _,_ € V.

We note that V is the energy space associated with the operator - A,

v = Dom((- n)I/2), and it is the closure of Dom(A) with respect to the energy

norm I-IV , which satisfies
-36-

(4.31) 1€12 = <L€,14> z = <-A_,¢>Z, _ _ Dom(A).

For each N = 2,3,... let A N denote the uniform partition of the interval

_ j_N-I

1 2 N-I I}. Let {eN_j= I denote the usual linear [0,i] given by {0, _ , N ''''' N '

B-splines on [0,i] corresponding to the partition AN and which satisfy

e_(0) = e (I) = 0, j = 1,2,... N-I. The e are given by

(4.32) e3 (8) = N((_- N 0) 0 e:: [_ , N ]

0 elsewhere

j = 1,2,...,N-I and are depicted in the figure below.

0 t I

o 1 1 N N N

Figure4.1

r j_N-I

Let ZN = spanteN2j= I . We note that ZNc V and dim ZN = N-I, N = 2,3, ....

Define PN:Z . ZN to be the orthogonal projection of Z onto ZN with respect to

the < , >Z inner product and PN : V + ZN to be the orthogonal projection of V

onto ZN with respect to the < _ >V inner product.

We define the operator A N : ZN . ZN as the inverse of the operator given
-37-

by

(4.33) A-I P A-I

N = N IZN"

The invertibility of A of course follows from the coercivity of - A (see

(4.14)). On the other hand, a straight forward calculation yields

(4.34) < ANI ZN, ZN> v = IzNIZ2 ' ZN _ ZN"

Consequently the operator AN 1 is invertible and the operator A N is well

defined. It is also self-adjoint. Indeed

(4.35) <_zN,YN> Z = - <zN,YN> v , zN,y N _ ZN.

From (4.35) it also follows that

(4.36) <ANZN,ZN>Z < -_IZNlZ2 , ZN _ ZN"

It can be concluded therefore, that AN is the infinitesimal generator of a

semigroup of bounded linear operators on ZN, { TN(S) : s _ O} with

TN(S) = exp(ANS) , s _ 0 and satisfying

(4.37) .N_ o.I!Cs_Iz _ e , s ) 0.

Elementary properties of spline functions (see [27]) imply PN . I

A-I as N . _, strongly on V. Furthermore, compact and
-38-

(4.38) IPN A-Iz - n-lzlz _ IPNA-Iz -A-Izlv = I(PN - I) n-lzlv

imply that PN A-I + A-I as N + _ in the uniform operator topology on i(Z).

We have therefore, that

(4.39) IANIPN - A-IIz + 0

as N + _.

From (4.37) and (4.39) we conclude (see [4], [15])

(4.40) TN(S)PNZ + T(s)z

and

(4.41) TN(S)PNZ + T (s)z

as N . = for each z _ Z, uniformly in s for s in compact subintervals.

With TN = TN(T)' QN = q0PN ' GN = g0PN and

(4.42) BN = (I - TN)PN_ + f_ TN(O)PNa'do ,

(4.40), (4.41) and elementary approximation properties of spline functions

guarantee that Hypotheses 3.1 and 3.2 hold and hence that the convergence

results for the finite time interval problem given in Theorem 3.3 apply.

For the infinite time interval Problem, (4.37) implies that Hypothesis

3.8 is satisfied. Moreover, if q0 > 0 the conditions given in the statement

of Theorem 3.10 are satisfied (see the Remark following the proof of Theorem 3.10)

and consequently the convergence results for the infinite time horizon problem
-39-

given in Theorem 3.9 hold.

Define the RN-I vector valued function EN on [0,I] by

(4.43) EN(8) T (e (8), e (8),..., e

8 _ [0,I] and the (N - I) x (N - I) Gram matrix associated with the basis

elements {e_} N-1 j=l

(4.44) MN = <EN, E_> Z •

A straight forward calulation yields

I

2 1

0 0

1 1 2 1

(4.45) MN = _ 60 _ _

!

6

A

3

!o6

1 2

0 0

6 3

Let the (N - I) x (N - i) matrix HN be given by

(4.46) HN = <EN, A_> Z •

Then

[HN]Ij = <eNi' AN eNJ >z

(4.47)

=- <e i , eJ>

N NV
-40-

=- <Le i, LeJ> Z

= - <aDe i, DeJ> Z .

In the case a(x) = a, x E [0,I], a constant, we have

m

-2 I 0 0

I -2 I

(4.48) HN = aN 0

\

o

1 -2 1

0 , 0 1 -2

The matrix representation [ AN ] for the operatornN is given by

-I

(4.49) [nN] = M N HN'

while for the operators TN, QN and GN we have

(4.50) [TN] = exp(_N]T),

(4.51) [QN] = q01N

and

(4.52) [GN] = g01N
-41-

where IN denotes the (N - i) x (N - i) identity matrix. If we define

_0N' a_ g RN-I by

(4.53) _0N = <EN'_0>Z

and

(4.54) _ = <EN,a'> Z,

j = 1,2,...,N-I we obtain

T

(4.55) [BN] = (IN - [TN])_I_0N + f0 exp([AN]°)MNI_ do

= (I N -[TN])MT.IlOoN + [AN]-I([TN ] - IN)MITIla_ .

When the finite dimensional approximating gain matrices [FN(t)], t =

0,1,2,...tf-I for the finite time interval problem have been computed using

(3.51)-(3.54), an approximation fN(t,') to the feedback kernal f(t,.) can be

obtained from

(4.56) fN(t,e) = EN(e)TMNI[FN(t)] T,

t = O,l,2,...tf-l, e g [0,I]. We have

(4.57) fN(t,-) + f(t,')
-42-

in H0(0,1) as N + _ for each t = 0,1,2,...tf-l.

Similarly, for the infinite time interval problem, an approximation fN to

f is given by

(4.58) fN(e) = EN(0)T_I[FN]T,

e _ [0,i] where the matrix [FN] is computed from (3.57) - (3.59). We have

(4.59) fN . f

in H0(0,1) as N + _ .

We demonstrate the feasibility of our schemes on an infinite time

interval problem of the form discussed above. Taking q0 = 1.0, r = 1.0,

a(x) = a = 1.0, x _ [0,I] and T = .01 we obtained the approximating feedback

kernals given in Table 4.2 and shown in Figure 4.3 below.
-43-

8 f3(8) f5(8) f9(8) fll(8) f13(8)

0.00 .0000 .0000 .0000 .0000 .0000

0.05 .0250 .0228 .0224 .0222 .0221

0. I0 .0500 .0456 .0449 .0445 .0443

0.15 .0750 .0684 .0673 .0668 .0664

0.20 .I000 .0911 .0898 .0891 .0887

0.25 .1250 .1147 .1122 .1115 .1109

0.30 .1500 .1382 .1346 .1338 .1331

0.35 .1727 .1617 .1571 .1561 .1554

0.40 .1908 .1853 .1798 .1787 .1778

0.45 .2088 .2129 .2028 .2013 .2005

0.50 .2269 .2405 .2286 .2261 .2245

0.55 .2450 .2681 .2543 .2513 .2496

0.60 .2630 .2958 .2819 .2793 .2769

0.65 .2811 .2956 .3097 .3056 .3035

0.70 .2584 .2954 .3246 .3275 .3275

0.75 .2153 .2953 .3331 .3352 .3385

0.80 .1723 .2951 .3169 .3260 .3263

0.85 .1292 .2213 .2697 .2163 .2170

0.90 .0861 .1475 .2097 .2163 .2170

0.95 .0431 .0738 .1049 .1125 .1175

1.00 .0000 .0000 .0000 .0000 0000

Table 4.2
-44-

I I I I I I I I I I I I I I I i I i J I I I I I I I I I I

._- ..<:&.7,.-:.._ N = 5

•e...i ." "% _,'l

"_-' -- ...Zh".... : '_,' 9

II { I ? 'II ' . Ii. I/

-- .,.$: ,',/

_."...... ./" ,,,,.,,.'_,

',+I-

_ ..Y, , 6 /--_ 13

,.../. ,. ,,;t,,/

..."_,.:.._ ". ,_

e-_,'r'_ -- "J_:!" ',. ,ll,

- ....,.Y ,_

P".-_"O -- .-_:'/"::'_ "_ " "'I!_'_I

/ :I

":":':: --/v/ / "!I

'XI-- Ill I I I I I I I I I 1

'X' 0"0 20"0 40"0 '--<0"0 :--:0"0 100"0

_.01

Figure 4.3
-45-

4.2. Hereditary or Time Delay Systems

In this example we consider linear hereditary control systems of the form

(4.60) x(s) = LXs + B0v(s) s > 0

where x(s) 8 Rn, xs E L2((-r,0);Rn ) for some r > 0, v 8 L2((0,Sl);Rm) for all

s 1 < _ and B0 _ L(Rm,Rn). The function x s represents the history on the

interval [s-r,s]; that is x (6) = x(s+e), e _ [-r,0]. The operator L is s

assumed to be of the form

_)

(4.61) L_ = [ Ai_(-r i) + f0_rA(O)_(e)de i=0

where A i € i(Rn), i = 0,1,2,...,v, A(.) _ L2((-r,0) ; L(RN)) and

0 = r 0 < r I < ... < rv = r. The initial data is given by

(4.62) x(0) = n, x0 = 4,

where n E R n and _ € L2((-r,0); Rn).

Once again, we formulate the dlscrete-tlme control problem by lettln£ T

denote the sampling interval and considering plecewlse constant control inputs

of the form

(4.63) v(s) = u(t) s £ [tT, (t+l)T),

t = 0,I,2,..., where u(t) € Rm for each t. The state space is

(4.64) Z = Rn x L2((-r,0);R n)
-46-

with the inner product

(4.65) <(q'_)' (_'_>Z = <_'_>R n + <_'_>L 2"

The state of the dlscrete-time control system is

(4.66) z(t) = (x(tr;n,_,v), xt_(_,_,v)) , t = 0,1,2,...,

where x(s;n,_,v) denotes the solution at time s to the system (4.60) - (4.62)

and x (n,_,v) its history on [s-r,s]. Then s

T(o) Bdou(t), t = 0,1,2, (4.67) z(t+l) = T(_)z(t) + f0 "'''

(4.68) z(0) = (n,_)

where {T(s) : s > 0} is the CO semigroup of bounded linear operators on Z

with infinitesimal generator A: Dom(n) c Z . Z given by

Dom (A) = {(_,_) _ Z : _ _ Hl((-r,O); Rn), _ = _(0)}

(4.69) A(_(O), _) = (L_, D_),

R m and the operator B: . Z is defined by

(4.70) Bu = (Bou, 0), u E Rm.

T T(o) Bdo, we obtain Letting T = T(T) and B = f0
-47-

z(t+l) = Tz(t) + Bu(t), t = 0,1,2,...

(4.71)

z(0)=

Let Q0' GO € i(R n) be symmetric and nonnegatlve and let R _ i(R m) be symmetric

and positive definite. We consider the dlscrete-time linear-quadratic

regulator problem with state given by (4.71) and performance index

tf-I

(4.72) J(G;O,tf,n,_,u) = _ {<Qz(t),z(t)>z+

t=O

<Ru(t),u(t)>Rm} + <Gz(tf),z(tf)> Z

where the operators Q : Z . Z and G : Z . Z are given by

(4.73) Q(_,_) = (Q0_,0)

and

(4.74) G(_,_) = (G0$,0)

respectively. For the infinite time problem we of course have

tf = _ and GO = 0.

For the finite time problem with the operators, T,B,Q and G as defined

above, we apply the theory developed in Section 2.1 and conclude that the

optimal control is given by

(4.75) u,(t) = -F(t)z,(t), t = O,l,2,...,tf-l,
-48-

where the operators F(t) g L(Z,R n) are given by (2.8) - (2.11). The operator

F(t) can be represented by a matrix of operators, [F0(t), Fl(t)] with

F0(t) _ L(Rn,R m) and Fl(t) € i(L2((-r,0); Rn); Rm). It follows from (4.70)

therefore that

(4.76) u,(t) = -f0(t)x,(tT) -fO r fl(t,8)(x,)tT(8)dS, t = 0,1,2,...,tf-l,

where f0(t) is an m x n matrix and fl(t,-) is a square integrable m x n

matrix valued function on (-r,0).

For the infinite time problem we assume that our original hereditary

system and QO are such that there exists an admissible control for each

z(O) = (n,#) _ Z and that Hypothesis 2.4 is satisfied. Then Theorems 2.3 and

2.5 imply that there exists a unique nonnegative self-adjoint solution H to

the Riccati algebraic equation (2.14) with the optimal control u, given by

(4.77) u,(t) = -Fz,(t), t = 0,1,2,...,

where F _ L(Z,R m) is given by (2.18) - (2.19). The feedback gain F can be

represented by a matrix of operators [F0,F1] with F0 _ L(Rn,R m) and

F 1 € i(L2((-r,0);Rn),Rm). We have

(4.78) u,(t) = -f0x,(tT) - f0-rfl(8)(x,)tr(8)d8' t = 0,I'2,"'''

where fO is an m x n matrix and fl is a square integrable

m x n matrix valued function on (-r,0).

Numerical methods for the approximate solution of the continuous time

linear quadratic regulator problem for hereditary systems in closed-loop form
-49-

have been studied extensively (see [7],[16],[19],[25],[26]). Most closely

related to the approximation framework which we have developed here for the

discrete-tlme problem are the treatments for the continuous-time problem in

[2],[II],[14] and [17]. The first approximation scheme applied to the

continuous-time linear-quadratic control problem for hereditary systems was

the AVE scheme used in [I],[11],[25] and [26] which approximates the history

by piecewise constant functions and derivatives with finite differences. In

[I], Banks and Burns cast AVE in its modern semigroup approximation form,

proved strong convergence of the approximating open-loop semigroups and used

the scheme to compute the open-loop optimal control. Gibson [II] used the

convergence results in [I] to prove strong convergence of the solutions to the

approximating continuous-time Riccati equations and uniform norm convergence

of the approximating feedback control laws. However, the rate of convergence

as it was observed in numerical studies was relatively slow. Spline based

schemes for continuous- time systems were developed for the finite time

interval problem in [17] and for the infinite time interval problem in [2].

With regard to the minimization of the performance index these schemes

represented an improvement over the AVE scheme. However, they appeared to

yield only weak convergence of the solutions to the approximating Riccati

equations and strong operator convergence of the approximating feedback gains.

Recently, a new spline-based state approximation for hereditary systems

has been proposed in [14]. When applied to the continuous time control

problem, this new method performs at the level of the earlier spline-based

schemes and yields strong convergence of the approximating feedback gains. We

have chosen this scheme to describe and implement here for the discrete-time

problem.
-50-

To simplifythe presentation,we considersystemsof the form (4.60)

having only a single discretedelay (v = I) and no distributeddelay term

(A(e)_ 0). A detaileddescriptionof the scheme in completegeneralitycan

be found in [14].

For each N = 1,2,...define

(4.79) e0 = (In,0) and ^jeN = (0,eJ), j = 1,2,...N+I,

where In denotes the n × n identitymatrix and the e_ are the

n x n matrix valued piecewiselinearelementson [-r,0)given by

N (8 + r -r 7 _)In e _ [ _ , 0)

e_(e) =

0 elsewhere

r r r

• Nr(8 + j g )In 8 s [-j g , -(j-l) g)

(4.80) e (e) = -Nr -- (8 + (j-2) _ )I n 8 s [-(j-l) _r , -(j-2). _)r

0 elsewhere

j = 2,3,4,...N
-51"

I-N N J r

(0 + (N-l) )I n 0 _ [-r,'(n-1) _)

N+I

eN (O) =

0 elsewhere

e

Let

N+I

(4.81) ZN = {z _ Z : z = _ ej e_ o_._ Rn}. j=0 ' 3

Note that dim ZN --n(N + 2). We shall refer to the collection {ej}j=0 N+I as a N+I

"basis" for ZN and a vector = _ x Rn as being a "coordinate vector" for an

j=O

element in ZN. Defining

(4.82) ET = (eN, ^0 eN,..- ^i ^N+I_ e )

we have

N+I

^T

(4.83) ZN = {z € Z : z = E_,_ = € x Rn}.

j=0

r^j_N+l

Let MN denote the Gram matrix corresponding to the basis teNlj= 0 . A

straight forward computation yields

r

(4.84) M N = diag (In, _ mN ®I n )

where the (N+I) x (N+I) matrix mN is given by
-52-

B

I I

o o

1 2 1

_ _ 0

(4.85) mN = 6 3 6

0

0

1 2 1

0

6 3 6

1 1

0 0

6 3

L..

and ® denotes the Kronecker product.

Let PN denote the orthogonal projection of Z on to ZN with respect to the

inner product (4.65). It follows that

(4.86) PN(_,_) = (_,pN_)

• j_N+I

where PN is the orthogonal projection of L2((-r,0);Rn ) on to spanleN_j= l with

respect to the usual L 2 inner product. We have

(4.87) PN(_,_) = ETcN(_,¢)

N+I

where _N(_,_) € x Rn is given by j=O

(4.88) CN(_,_) = M_I col (_'¢N'1 _N"''_N 2 N+I )

with

(4.89) *J = fO r eJ(e)_(O)d@.
-53-

We shall set TN = exp(ANT) where A N : ZN . ZN is an appropriately

defined finite dimensional approximation to the operator A given by (4.69).

Noting that ZN _ Dom(A), we motivate the definition of AN by first formally

extending the operator A to an operator defined on ZN.

For zN = (_N,_N) € ZN define

(4.90) ANzN = (A0_N+AICN(-r)' D+_N + _(_N - lim-_N(8))) 8+0

where _ is the Dirac delta impulse concentrated at zero and D% denotes the

right hand derivative of _. For each N = 1,2... let A N : ZN . ZN be given by

+

(4.91) ANZ N = (A0_N+AI_N(-r), pN D _N ) + _N(_N - lim__N(O))

8+0

where

AT

(4.92) 6N = ENYN

with

(4.93) YN = MNI col (0, lim_ eN1(8), ..- , lim_ eN+Ii_)). .^

O.0 O+0

To compute [AN], the matrix representation for the operator AN, we let

[zN] denote the coordinate vector representation for an element zN g ZN.

Then from

(4.94) [ANZN] = [AN][ZN]
-54-

and

(4.95) [ANZN] = _1 aN(ANZN)--MN I- _N,ANZN> Z = MNI- _N_NET>z[zN]

we obtain

I

(4.96) [AN] = MNI HN

where

^ ^T

(4.97) HN = <EN,ANEN>z .

^

Using the definitions of A N and EN a straight forward calculation (see

[14]) yields

B

M

A0 0 0 A 1

I

n

(4.98) HN =

0

bN® I n

_ m

where the (N+I) x (N+I) matrix hN is given by
-55-

I i

I 1

-_ -_ o o

/

1 1 0 / (4.99) _ = _ 0 - _ _

I

o

o

1 1

o - o m I

2 2

1 1

0 0 --2 2 "

i I

The matrix representation for the operator TN = exp (ANt) can then be

computed from

(4.100) [TN] = exp([AN]_).

We define the operators BN : Rm . ZN, QN : ZN . ZN and GN : ZN . ZN by

(4.101) BN = PN B

(4.102) QN = PNQ

and

(4.103) GN = PNG

from which we obtain
-56-

(4.104) [BN] = col (Bo, 0,...,0) m Rn(N+2)xm

m

Q0 o o

0 o

(4. 105) [QN] = I _ _ Rn(N+2)xn(N+2)

0 0

m m

and

m

GO 0 0

0 0

(4. 106) [GN] = I _ € Rn(N+2)xn(N+2)

o 0

D

Finally, defining

(4.107) BN = fO exp (_s) BNdS

we have

(4.108) [BN] = f_ exp (_N]S)[_]ds.

Once the matrix representations for the approximating feedback gains have

been computed, [FN(t)] , t = 0,I,2,... tf-I from (3.51) - (3.54) for the finite

time interval problem and [FN] from (3.57) - (3.59) (assuming, for the moment,

that solutions to (3.34) exist) for the infinite time interval problem,
-57-

approximations for f0, fl(t,.), t = 01,2,...tf-1 and f0,fl(.) can be

computed from

(4.109) ((fN(t)) 0 , (f (t,')) T) = [FN(t)] T, t = 0,1,2,...,tf-l,

and

(4.110) ((fN) =

respectively.

For the approximation scheme defined above, it is shown in [14] that

PN . I strongly on Z. Using a Trotter-Kato like result it is also shown that

(4.111) exp _NS)P N + T(s)

and

(4. 112) exp(ANS)P N . T (s)

strongly on Z and uniformly in s for s in compact intervals. Hypothesis 3.1

is a simple consequence of these results. The present scheme, therefore,

satisfies all of the hypotheses of Theorem 3.3 and we may conclude that the

convergence results for the finite time interval problem given in the

statement of the theorem hold. In particular, we have

(4.113) f_(t) + f0(t)
-58-

in Rmxn and

(4.114) f_(t,.) + fl(t,.)

in L2((-r,0);R mxn) for each t = 0,1,2,... tf-l.

With the operators Q and G given by (4.73) and (4.74) and the operators

QN and GN defined as in (4.102) and (4.103) it is clear that the hypotheses

given in the statement of Theorem 3.7 are satisfied. We have therefore that

tf

for the present example the operators {H(t)}t= 0 are trace class and

(4.115) lim II_N(t)PN _(t)ll1 = 0, t = 0,1,2,...,tf. N+_

For the infinite time problem and the approximation scheme discussed

here, the situation with regard to convergence is much the same as it is for

the continuous time problem (see [14]). We are unable to demonstrate the

existence of an M and an r < 1 for which (3.38) and (3.39) hold. In fact, our

numerical studies point to the conclusion that condition (3.39) is violated by

the present scheme. We observe the existence of a sequence of closed-loop

eigenvalues of the approximating discrete-time control problems (p2N) which

tend toward the unit circle as N + =.

On the other hand, upon solving the approximating problems it is also

apparent that I_N[ remains bounded in N. Consequently we may apply Theorem

3.12 to conclude that a solution _ to (2.14) does in fact exist,

_NPN + H weakly and FNP N + F strongly as N + = . We have therefore that

0 f0 (4.116) fN +
-59-

in R m×n and

1 fl

(4.117) fN +

weakly in L2((-r,0);R mxn) as N + _.

We applied the scheme to the infinite-time problem with state

(4.118) y(s) + y(s-l) = v(s).

Transforming (4.118) to an equivalent first order system we obtain a system of

the form given in (4.60), (4.61) with n = 2, r = I, m = I,

A0 = 0 ' A1 = - ' B0 = , A = 0 and x(s) = .

y(s)J

Taking Q0 = I_ Ol I the performance index takes the form

(4.119) J(0;0,_,y(0),y(0),y0,Y0,U ) = _ y(tT) 2 + _ (tT)2+ Ru(t) 2

t=O

where T is the length of the sampling interval. The optimal feedback control

is given by

u,(t) = _[f0]lY(t T) _[f0]2_(t T)

(4.120)

- f01{[fl(0)]lY(tT+0 ) + [fl(0)]2Y(tT+0)}d0

where [fO]i and [fl(8)]i,i = 1,2 are the ith componentsof the 1 x 2

matrices fO and fl(8).
-60-

Since by taking the initial conditions

(4. 121) y(0) = 0, ;(0) = 0, yo(8) = 0, -i _ e _ 0

we have y(s) = 0, s > 0 regardless of how ;0(8), -1 _ 8 < 0 is chosen, it

follows that [fl(e)] 2 = 0, -I _ e _ 0. Indeed, the optimal control

corresponding to the initial conditions (4.121) with Y0 arbitrary is u(t) =

0, t = 0,1,2, .... Furthermore, the nature of the approximation scheme is

such that we must have [f_(e)] 2 = 0, -I _ e _ 0, N > I.

Setting T = .01 we obtained the results given in Tables 4.4 and 4.5 and

Figure 4.8 below when R = .05. With R = 1.0, the results given in Tables 4.6

and 4.7 and Figure 4.9 were obtained. As the cost of control increases the

effect that the optimal control for the infinite dimensional problem has on

higher modes decreases. Consequently, the finite dimensional approximations

are more effective and convergence is more rapid.
-61-

N 2 4 8 I0

[f_]l 4.5483 4.5452 4.5451 4.5451

[f$]2 5.2954 5.2948 5.2948 5.2948

Table 4.4

°

.00 .1277 .0942 .0878 .0872

-.05 .0948 .0690 .0659 .0658

-.I0 .0619 .0437 .0439 .0445

-.15 .0289 .0185 .0163 .0154

-.20 -.0040 -.0068 -.0171 -.0137

-.25 -.0369 -.0321 -.0506 -.0556

-.30 -.0698 -.1046 -.1054 -.0976

-.35 -.1027 -.1772 -.1603 -.1650

-.40 -.1357 -.2497 -.2375 -.2324

-.45 -.1686 -.3223 -.3371 -.3361

-.50 -.2015 -.3949 -.4367 -.4399

-.55 -.5891 -.6156 -.6090 -.6054

-.60 -.9767 -.8363 -.7813 -.7709

-.65 -1.3643 -1.0570 -1.0200 -1.0283

-.70 -1.7519 -1.2777 "1.3251 -1.2858

-.75 -2.1395 -1.4984 -1.6301 -1.6899

-.80 -2.5271 -2.1836 -2.1603 -2.0941

-.85 -2.9147 -2.8688 -2.6905 -2.7232

-.90 -3.3023 -3.5540 -3.4158 -3.3523

-.95 -3.6899 -4.2391 -4.3361 -4.3269

-I.00 -4.0775 -4.9243 -5.2565 -5.3017

Table 4.5
-62-

N 2 4 8 10

[fO] I 1.4050 1.4054 1.4054 1.4054

[f0N]2 1.9477 1.9479 1.9479 1.9479

Table 4.6

8 [f_(8)] I [f_(8)] I [f_(8)] I [flO1 (e)]1

.00 -.2813 -.2831 -.2835 -.2835

-.05 -.3433 -.3402 -.3383 -.3379

-.I0 -.4052 -.3972 -.3931 -.3923

-.15 -.4052 -.4543 -.4522 -.4534

-.20 -.5292 -.5114 -.5156 -.5146

-.25 -.5911 -.5684 -.5790 -.5798

-.30 -.6531 -.6439 -.6478 -.6450

-.35 -.7151 -.7195 -.7165 -.7179

-.40 -.7770 -.7950 -.7902 -.7907

-.45 -.8390 -.8705 -.8689 -.8685

-.50 -.9010 -.9460 -.9476 -.9462

-.55 -1.0044 -1.0347 -1.0328 -1.0324

-.60 -1.1078 -1.1233 -1.1181 -1.1185

-.65 -1.2112 -1.2120 -1.2089 -1.2103

-.70 -1.3146 -1.3007 -1.3053 -1.3021

-.75 -1.4180 -1.3894 -1.4016 -1.4030

-.80 -1.5214 -1.5013 -1.5057 -1.5040

-.85 -1.6248 -1.6132 -1.6098 -1.6112

-.90 -1.7282 -1.7251 -1.7198 -1.7185

-.95 -1.8316 -1.8370 -1.8357 -1.8352

-I.00 -1.9350 -1.9489 -1.9516 -1.9519

Table 4.7
-63-

0.0

2

4

8

N = 10

-----...---------.: . ..----;:::.-'"='=.---

.-z.....R-~:::-· .._·- --...,.""-

./~:~..~_.:_:~-~~~:~..;,;~~ ......_~_·......-:_;·__-----N =

,"

lSI

( ..

'r'". ( .. :+: (r·.

I

...-j-

-:5:1·( ..

I

'.J)

.......

t··)

I

0:-

( ..·

V

I

C::'

v·

1.0

I

:0:.01.

Figure 4.8

-:::a3. (1

tS,1

_i

I

........ (0

:"': V·

I

"J)

(I)·

I

V

(,~·

-<

( ..

..J:'·

....... I

-=::'

I:::;)·( ..

I

Figure 4.9
-64-

4.3 Control of a Flexible Structure

We consider an Euler-Bernoulll beam cantilevered to a rigid hub which is

free to rotate about its fixed center, point O. Also, a point mass m I is

attached to the other end of the beam. The control is a torque u applied to

the hub, and all motion is in the plane. See Figure 4.10 and Table 4.11.

Figure 4.10
-65-

r = hub radius I0 in

£ = beam length I00 in

I0 = hub moment of inertia about axis

perpendicular to page through 0 100 slug in2

mb = beam mass per unit length .01 slug/in

m I = tip mass I slug

EI = product of elastic modulus and

second moment of cross section for beam 13,333 slg in3/sec 2

fundamental frequency of undamped structure .9672 rad/sec

Table 4.11

The angle 8 represents the rotation of the hub (the rigid-body mode),

w(t,n) is the elastic deflection of the beam from the rigid-body position,

and wl(t) is the displacement of m I from the rigid-body position. For

technical reasons, we do not yet impose the condition Wl(t ) = w(t,_).

The control problem is to stabilize rigid-body motions and linear (small)

transverse elastic vibrations about the state 0 = 0 and w = 0 . Our linear

model assumes not only that the elastic deflection of the beam is linear but

also that the axial inertial force produced by the rigid-body angular velocity

has negligible effect on the bending stiffness of the beam. The rigid-body

angle need not be small.

For this example, it is a straight forward exercise to derive the coupled

ordinary and partial differential equtions of motion in 0, w and w I. However,

rather than writing these equations explicitly, it is easier and more useful
-66-

for our purposes to derive an abstract second order evolution equation for the

structure. To do this, we define the generalized displacement vector

(4.122) x = (8,W,Wl) _ H = R 1 x L2(0,£ ) x R I.

The kinetic energy in the system is then

(4.123) Kinetic Energy = i/2 <M0x,x> H

where M0 is the unique bounded self-adjolnt linear mass operator M 0 on H such

that

A A _ A A A

(4.124) <M0x,x> = 1088 + mb<w + _08,w + _08>L2 + ml(w I + _0(£)8)(wi + _0(£)8),

where _0 € L2(0,£ ) is given by _0(_) = r + n. It is easy to show that M 0 is

also coercive. The elastic strain energy is

(4.125) Strain Energy = i/2 a(x,x)

with

(4.126) a(x,x) = E1 <D2w, D2W>L2 .

We make a(',') into an inner product by setting

A A

(4. 127) <x,x> v = a(x,x) + 88
-67-

and define the straln-energy space

(4.128) V = {x = (O,{,{(_)): { s H2(0,£), {(0) = D{(0) = 0}.

The last term in (4.127) is necessary for the V-inner product because there is

no strain energy associated with the rotation of the hub.

We define the stiffness operator A0 by

(4.129) Dom(A 0) = {x = (e,{,{(%)) s v: { _ H4(0,£), D2{(%) = 0}

and

0 0 0 1

(4. 130) A0 = 0 E1 D4 0 .

0 -El D3 0

This operator is self-adjoint with compact resolvent and all positive

eigenvalues except the one zero eigenvalue corresponding to the rlgid-body

mode. Note that V is the domain of the _square root of AO.

With these mass and stiffness operators, we can write the equations of

motion as

(4.131) M0x(s) + c0A0x(s) + A0x(s) = B0u(s) , s _ 0,

where cO is a positive constant and the term coAo_ represents viscoelastic

damping in the beam. The input operator is
-68-

(4.132) B0 = (I,0,0).

A A _ A

Letting Z = V x H with inner product <(v,h),(v,h)> Z = <v,v> v + <M0h,h>H,

the first order form of this system is given by

(4.133) z(s)= Az(s)+ Bu(s), s ; 0,

where z = (x,x) g Z and A is the unique extension of the operator

(4.134) _ I 0 I 1 o

= , Dom(A) = Dom(A 0) x Dom(A 0),

that generates a Co-semigroup on the space Z. Of course, B is

E°]

(4.135) B = .

Mo1Bo

See [I0] and [12]. The hub-beam-tip mass structure here is discussed in more

detail in [12], along with the continuous-time problem.

The dlscrete-time control system for sampling interval T is

(4.136) z(t+l) = Tz(t) + Bu(t), t = 0,1,2,...,

where

T

(4.137) T =T(T), B = fO T(s)B ds

and {T(s): s _ O} is the semigroup generated by the A in (4.133).
-69-

As in the previous examples, we will solve a discrete-time optimal

control problem on the infinite interval. In the performance index, we take

the state weighting operator Q to be the identity on Z. This means that

<Qz,z> Z is twice the total energy in the structure plus the square of the

rigid-body rotation. Since there is one input, the control weighting R is a

scalar. The optimal control has the feedback form

(4.138) u,(t) = - <f,x(t)> v - <M0g,x(t)> H

where x(t) has the form (4.122) and

(4.139) f = (fl ,f2,f3 ) s V

(4.140) g = (gl ,g2 ,g3 ) _ H .

Our approximation of the structure is based on a finite element

approximation of the beam which uses Hermite cubic splines as basis functions

([27]). We define the sequence of spaces VN = span {e_}j_IN = with

1

(4.141) eN = (I,0,0),

(4.142) e = (0,@N,@N(£)), j J j = 2,3,... JN'

where the

_'s are the cubic splines. Each VN is a subspace of V, and our

approximation scheme is a Ritz-Galerkin approximation obtained by projecting

(4.131) onto VN. See [12] for details. Writing
-70-

JN

(4.143) XN(S) = _ [XN(S)]je_, j=1

we have

(4.144) MN [_(s)] + Co_ [XN(S)] + _[XN(S)] = BoNU(S )

to solve for the vector [XN(S)] of time-dependent coefficients [XN(S)] j. The

mass matrix MN and the stiffness matrix KN are given by

(4.145) [MN]ij = <MOeiN, jeN>H, [KN]ij = <el, eJ>v

and the input matrix is

(4.146) BON = [i 0 0 ... 0]T.

With zN = (XN,XN) _ VN x VN, (4.144) is the matrix representation of an

evolution equation

(4.147) ZN(S) = ANZN(S) + BNU(S )

where A N and B N approximate Aand B. It is shown in [12] that, as N

increases, the semigroup {TN(S): s _ O} generated by A N converges strongly to

the semigroup {T (s) : s ) 0} and that the adjoint semigroup {T_(s): S ) 0}

converges strongly as well. Since B N is the Z-projection of B onto VN x VN,

it converges strongly to B .

For the approximating discrete-time control systems, we replace z(t), T,
-71-

B, T(') and B in (4.136) and (4.137) with ZN,TN,BN, TN(.) and BN,

respectively. For each N, the solution to the infinlte-time optimal control

problem is based on the Nth Riceatl operator equation (3.33). As in the

A

previous examples, we solve the Riccati matrix equation (3.59) for _N' which

is related to [_N ] (the matrix representation of the operator _N ) as in

Section 3.3, except here we have

A

(4. 148) nN = WN[RN] ,

where

and % is the stiffness matrix with I added to the first element. Since Q =

I in the infinite dimensional problem, QN is the identity on VN x VN and it

follows from (3.48) that the matrix QN for (3.59) is W N.

The optimal feedback control for the N th problem is given by (3.56) with

the matrices in (3.57) and (3.58), and it has the equivalent representation

(4.150) u_(t) = - <fN,XN(t)>v - <M0gN,XN(t)> H

where

1 2 3

(4.151) fN -- (fN'fN'fN) s V,

1 2 3

(4.152) gN = (gN'gN'gN) g H
-72-

as in (4.139) and (4.140). From (3.56), (4.143) and (4.150), it follows that

(4.153) = - [FNIT ,

gN E

where ETN = (e_ ,eN,...,e 2 JNN ).

For the sampling interval T = .01, the damping coefficient cO = .001 and

the control weighting R = 1, Tables 4.12-4.15 give the values of the

i i

corresponding scalar and functional gains, fN' gN' i = 1,2,3 for various

values of N. The values of the functional gains D2"2 2 _N and gN along the length

of the beam also are plotted in Figures 4.16 and 4.17. We tabulated and

22

plotted D fN because this is what appears in the V inner product in (4.150)

and also to show the H 2 convergence. We note that the form of the V inner

3

product given by (4.127) is such that fN does not appear in the feedback law.
-73-

N 3 4 5 7

1 .9991 .9992 9990 9992 fN " "

1

gN .1030 .1040 .1043 .1044

Table 4.12

N 3 4 5 7

f_ .1750 .1769 .1774 .1777

3

__ -18.1231 -18.3385 -18.3902 -18.4158

Table 4.13
-74-

2 2

D2f3(n)xl 2 05 D2f4(n)xl 2 05 D2f_(n)x10 5 D f7(n)x105

0.0 6.0266 7.0718 7.6281 8.1441

5.0 4.8991 5.4251 5.6142 5.6817

I0.0 3.7717 3.7784 3.6003 3.2193

15.0 2.6442 2.1317 1.5865 1.1727

20.0 1.5167 .4850 .2048 .7757

25.0 .3893 .1353 .4057 .3786

30.0 -.7382 .4896 .6066 .4481

35.0 1.0908 .8438 .8075 .8303

40.0 1.2767 1.1980 1.3790 1.2126

45.0 1.4626 1.5523 1.5713 1.6121

50.0 1.6484 2.0984 1.7636 1.8049

55.0 1.8343 1.9858 1.9559 1.9976

60.0 2.0202 1.8731 2.0543 1.9613

65.0 2.2061 1.7605 1.8179 1.8003

70.0 1.6333 1.6478 1.5815 1.6392

75.0 1.3545 1.3288 1.3450 1.3391

80.0 1.0757 1.0562 1.0303 1.0505

85.0 .7970 .7837 .7667 .7620

90.0 .5182 .5112 .5031 .4921

95.0 .2394 .2387 .2395 .2400

100.0 -.0394 -.0339 -.0242 -.0122

I

Table 4.14
-75-

2 2 2 2

n g3(n) g4(n) g5(n) g7(n)

0.0 .0000 .0000 .0000 .0000

5.0 -.2348 -.2576 -.2673 -.2745

I0.0 -.8710 -.9445 -.9727 -.9903

15.0 -1.8063 -1.9320 -1.9715 -1.9854

20.0 -2.9383 -3.0915 -3.1190 -3.1171

25.0 -4.1649 -4.2941 -4.2940 -4.3031

30.0 -5.3837 -5.4478 -5.4540 -5.4733

35.0 -6.4982 -6.5450 -6.5720 -6.5842

40.0 -7.5155 -7.5840 -7.6210 -7.6289

45.0 -8.4647 -8.5629 -8.5918 -8.6032

50.0 -9.3609 -9.4797 -9.5044 -9.5179

55.0 -10.2188 -10.3493 -10.3761 -10.3910

60.0 -11.0535 -11.1963 -11.2240 -11.2407

65.0 -11.8797 -12.0341 -12.0666 -12.0838

70.0 -12.7146 -12.8759 -12.9158 -12.9327

75.0 -13.5719 -13.7350 -13.7790 -13.7988

80.0 -14.4502 -14.6201 -14.6632 -14.6861

85.0 -15.3469 -15.5278 -15.5720 -15.5946

90.0 -16.2595 -16.4533 -16.5000 -16.5232

95.0 -17.1858 -17.3917 -17.4414 -17.4659

I00.0 -18.1231 -18.3385 -18.3902 -18.4158

Table 4.15
1:-;':1

(---I·

((I

........ I::;)

;~': -0::1--·

-oJ)

1:::;:-

-.J)·

-::r

1:S:-

CO

(--~

I:S:I

I::::~I·,-i

c:;:,

((I

I

-76-

:0:1

Figure 4.16

1::';::1

I~:~I·

,-I ((I

;~': V:O·

I

..J)

roo.·

I

-::r

.......

,-I

I~--~

ll")·

.......

I

I:~I

0-··

.......

I

--~

0.0 20.0 4(1. (1 .~.O. 0

Figure 4.17

::;:121. (1 100.0
-77-

5. Concluding Remarks

We have presented an approximation theory for numerical solution of the

discrete-time optimal linear regulator problem in Hilbert space, on both

finite and infinite tlme intervals. The motivation for thls theory comes from

optimal control problems for systems involving diffusion equations, hereditary

differential equations and distributed models of flexible structures. We have

demonstrated the application of the theory to examples from all three areas

The solution to the infinite dimensional optimal control problem is based

on an infinite dimensional Riccati operator equation -- a difference equation

in the finlte-time problem and an algebraic equation in the inflnite-tlme

problem. We have shown that the solution to the infinite dimensional problem

can be approximated by the solutions to a sequence of finite dimensional

problems each of which involves a finite dimensional Riccati matrix equation

to be solved numerically. The finite dimensional problems are just the

corresponding optimal control problems for finite element approximations to

the infinite dimensional control system. For the infinlte-time problem, the

finite dimensional Riccati equations usually are solved vla eigenspace

decomposition of the Hamiltonian matrix.

In both continuous and discrete-tlme optimal regulator problems for

distributed systems, the numerical solution often involves solution of large

Riccati matrix equations. As we observed at the beginning of Section 4, the

asymptotic relationship between the eigenvalues of a contlnuous-time

Hamiltonian system and the eigenvalues of the corresponding discrete-time

Hamlltonlan system is exponential. This means that the approximating finite

dimensional discrete-time Rlccati equations for a given distributed system

invariably are not as well conditioned as the corresponding contlnuous-time

Riccati equations. Nonetheless, as our examples should illustrate, the
-78-

numerical solution of such problems is well within the reach of current

computing. To emphasize this, we obtained all of the numerical results in

this paper on an IBM Personal Computer (not an XT or AT) with 640K of random

access memory and an Intel 8087 math coprocessor chip. The largest Riccati

matrix equation that we solved here was a 30 x 30 steady state equation for

the hub-beam-tip mass example. This solution takes 15 to 20 minutes on the

PC. We have solved much larger Riccati equations easily on larger mainframe

computers.
-79-

REFERENCES

[I] H.T. Banks and J.A. Burns, Hereditary control problems: Numerical

methods based on averaging approximations, SlAM J. Control and Opt.

16 (1978), 169-208.

[2] H.T. Banks and K. Kunlsch, The linear regulator problem for parabolic

systems, SlAM J. Control and Opt. 22 (1984), 684-698.

[3] H.T. Banks, I.G. Rosen and K. Ito, A spllne based technique for

computing Riccatl operators and feedback controls in regulator

problems for delay equations, SlAM J. Scl. Stat. Comput., 5 (1984),

830-855.

[4] A. Belleni-Morante, Applied Semlgroups and Evolution Equatlons_

Clarendon Press, Oxford, 1979.

[5] R.F. Curtain, The inflnlte-dlmenslonal Riccatl equations with

application to afflne hereditary differential systems, SlAM J.

Control and Opt., 13 (1975), 1130-1143.

[6] R.F. Curtain and A.J. Prltchard, The infinite dimensional Riccati

equation for systems defined by evolution operators, SlAM J. Control

and Opt., 14 (1976), 951-983.

[7] M.C. Delfour, The linear quadratic optimal control problem for

hereditary differential systems: Theory and numerical solution,

Appl. Math. and Optlm., 3 (1977), 101-162.

[8] M.C. Delfour, C. McCalla and S.K. Mitter, Stability and the infinitetime quadratic cost problem for linear hereditary differential

systems, SlAM J. Control and Opt. 13 (1975), 48-88.

[9] J.S. Gibson, The Riccatl integral equations for optimal control

problems in Hilbert spaces, SlAM J. Control and Opt. 17 (1979), 537-

565.

[I0] J.S. Gibson, An analysis of optimal modal regulation: convergence

and stability, SlAM J. Control and Opt, 19 (1981), 686-707.

[II] J.S. Gibson, Linear quadratic optimal control of hereditary

differential systems: infinite dimensional Riccatl equations and

numerical approximations, SlAM J. Control and Opt. 21 (1983), 95-139.

[12] J.S. Gibson and A. Adamian, Approximation theory for optimal LQG

control of flexible structures, Report, DeDartment of Mechanical,

Aerospace and Nuclear Engineering, University of California, Los

Angeles, Los Angeles, CA (1986).

[13] W. Green and T.D. Morley, Operator means, fixed points and the norm

convergence of monotone approximates, Mathematica Scandinavlca, to

appear, (1986).
-80-

[14] F. Kappel and D. Salamon, Spline Approximation for retarded systems

and the Riccati Equation, SIAM J. Control and Opt., to appear,

(1986).

[15] T. Kato, PerturbationTheory for Linear Operators,Springer-VerlaK,

New York, 1966.

[16] N.N. Kravsovskii,The approximationof a problem of analytic

design of controls in a system with time-lag, J. Appl. Math. and

Mech., 28 (1964), 876-885.

[17] K. Kunisch, Approximation schemes for the linear quadratic optimal

control problem associated with delay equations, SIAM J. Control and

Opt. 20 (1982), 506-540.

[18] A.J. Laub, A Shur method for solving algebraic Riccati equations,

IEEE Trans. Automat. Contr., AC-24 (1979), 914-921.

[19] E.B. Lee and A. Manitius, Computational approaches to synthesis of

feedback controllers for multivariable systems with delays, Proc.

1974 IEEE Conference on Decison and Control, 1974, 791-792.

[20] K.Y. Lee, S. Chow and R.O. Barr, On the control of discrete-time

distributed parameter systems, SIAM J. Control and Opt., I0 (1972),

361-376.

[21] J.L. Lions, Optimal Control of Systems Governed by

Partial Differential Equations, Springer-Verlag, New York, 1971.

[22] D.L. Lukes and D.L. Russell, The quadratic criterion for distributed

systems, SIAM J. Control and Opt., 7 (1969), 101-121.

[23] A. Pazy, Semigroups of Linear Operators and Applications to Partial

Differential Equations, Springer-Verlag, New York 1983.

[24] T. Pappas, A.J. Laub and N.R. Sandell, Jr., On the numerical solution

of the discrete-time algebraic Riccati equation, IEEE Trans. Automat.

Contr., AC-25 (1980), 631-641.

[25] D.W. Ross, Controller design for time lag systems via a quadratic

criterion, IEEE Trans. Automat. Contr., AC-16 (1971), 644-672.

[26] D.W. Ross and I. Flugge-Lotz, An optimal control problem for systems

with differential-difference equation dynamics, SlAM J. Control and

Opt. 7 (1969), 609-623.

[27] M.H. Schultz, Spline Analysis, Prentice Hall, Englewood Cliffs, N.J.

1973.

[28] J. Zabczyk, Remarks on the control of discrete-time distributed

parameter systems, SIAM J. Control and Opt. 12 (1974), 721-735.

1. Report No. NASA'CR-178081 \2. Government Accession No. 3. Recipient's Catalog No.

ICASE Report No. 86-15

4. Title and Subtitle 5. Report Date

NUMERICAL APPROXIMATION FOR THE INFINITE-DIMENSIONAL March 1986

DISCRETE-TIME OPTIMAL LINEAR-QUADRATIC REGULATOR 6. Performing Organization Code

P1UHH.RM

7. Author(s) 8. Performing Organization Report No.

J. S. Gibson and I. G. Rosen 86-15

9. Performing Organization Name and Address 10. Work Unit No.

Institute for Computer Applications in Science

and Engineering 11. Contract or Grant No.

Mail Stop 132C, NASA Langley Research Center NASl-17070, NASl-18107

Hampton, VA 23665-5225 13. Type of Report and Period Covered

12. Sponsoring Agency Name and Address Contractor Report

National Aeronautics and Space Administration 14. Sponsoring Agency Code

Washington, D.C. 20546 505-31-83-01

15. Supplementary Notes

Langley Technical Monitor: Submitted to SIAM J. Cont. Opt.

J. C. South

Final Report

16. Abstract

An abstract approximation framework is developed for the finite and

infinite time horizon discrete-time linear-quadratic regulator problem for

systems whose state dynamics are described by a linear semigroup of operators

on an infinite dimensional Hilbert space. The schemes included in the

framework yield finite dimensional approximations to the linear state feedback

gains which determine the optimal control law. Convergence arguments are

given. Examples involving hereditary and parabolic systems and the vibration

of a flexible beam are considered. Spline-based finite element schemes for

these classes of problems, together with numerical results, are presented and

discussed.

17. Key Words (Suggested by Authors(s)) 18. Distribution Statement

linear quadratic regulator, 64 - Numerical Analysis

discrete time, distributed system, 67 - Theoretical Mathematics

approximation

Unclassified - unlimited

19. Security Classif.(of this report) 120. Security Classif.(of this page) 21. No. of Pages 122. Price

Unclassified Unclassified 82 A05

For sale by the National Technical Information Service, Springfield, Virginia 22161

NASA Langley Form 63 (June 1985)

NASA-Langley, 1986