NASA Contnletor Report 178033

lCASE REPORT NO. 85-58

leASE

NASA-CR-178033

19860009522

FACTORIZATION AND REDUCTION METHODS

FOR OPTiMAL CONTROL OF DISTRIBUTED PARAMETER SYSTEMS

John A. Burns .-.. :".,....... . !- (" :-; .: ',,! .

..' :

Robert K. Powers

- r-- -,-,' .1' r" J..::..~ ... ~ll 1~,-~_.-,

Contract Nos. NASI-1707D, NASI-18ID7

November 1985

INSTITUTE FOR CO~WUTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

National Aeronautics and

Space Administration

Langley Research cent.

Hampton, Virginia 23665

FACTORIZATION AND REDUCTION METHODS

FOR OPTIMAL CONTROL OF DISTRIBUTED PARAMETER SYSTEMS

John A. Burns

Virginia Polytechnic Institute and State University

Robert K. Powers**

Institute for Computer Applications in Science and Engineering

ABSTRACT

A Chandrasekhar-type factorization method is applied to the linearquadratic optimal control problem for distributed parameter systems. An

aeroelastic control problem is used as a model example to demonstrate that if

computationally efficient algorithms, such as those of Chandrasekhar-type, are

combined with the special structure often available to a particular problem,

then an abstract approximation theory developed for distributed parameter

control theory becomes a viable method of solution. A numerical scheme based

on averaging approximations is applied to hereditary control problems.

Numerical examples are given.

Research was supported in part by NSF Grant No. ECS-8109245, AFOSR Grant No.

AFOSR-85-0287. Partial support was provided under NASA Contract Nos. NASI17070 and NASI-18107 while the author was a visitor at ICASE, NASA Langley

Research Center, Hampton, VA 23665-5225.

Research was supported in part by the National Aeronautics and Space

Administration under NASA Contract Nos. NASI-17070 and NASI-18107 while the

author was in residence at ICASE, NASA Langley Research Center, Hampton, VA

23665-5225 and in part by NSF Grant No. ECS-8109245 while in residence at

Virginia Polytechnic Institute and State University, Blacksburg, VA in

connection with his Ph.D. Thesis.

i

I. INTRODUCTION

The problem of developing computational algorithms for distributed

parameter control systems has been the subject of a large number of recent

papers. These articles run the gamut from very abstract papers dealing with

general approximation theory for infinite dimensional systems to explicit

numerical algorithms derived expressly for a particular type of system and/or

application. In any investigation of approximation schemes for such systems,

we believe that it is important to keep in mind the ultimate goal of the

approximation. For example, an approximation scheme that produces excellent

results if used for parameter estimation might be inappropriate for computing

feedback gain operators. The point at which approximations are introduced

into the analysis is something that also varies with the problem and with the

particular approach used to analyze the algorithm.

In this paper we concentrate on the linear quadratic optimal control

problem for certain distributed parameter systems. We employ the

approximation theory developed by Gibson [25] to formulate and analyze fast

computational algorithms for approximating optimal feedback gain operators.

These methods are based on factorization schemes of Chandrasekhar type.

A primary objective of this paper is to show that if one combines a

computationally efficient algorithm such as a Chandrasekhar type method with

the special structure of often available in "real problems," then many

practical problems can be attacked using distributed parameter control theory

and sound computational techniques. As noted in Casti and Ljung [15], this

has been considered a major stumbling block between theoretical results for

distributed parameter systems and their application to practical problems.
-2-

Our approach is based on an approximation theory specifically developed

for approximating infinite dimensional control systems. Although the paper

does contain some theoretical results on the existence, uniqueness, and

smoothness of strong solutions to Chandrasekhar integral equations, we feel

that the most significant aspects of this paper lie in the numerical results

and applications. However, such theoretical results (especially smoothness)

play an important role in the analysis of numerical schemes for direct

integration of these equations (see Sorine [42]). These results have also

been used to establish differentiability of strong solutions to operator

Riccati differential equations (see [27], [44]).

Previous authors have studied infinite dimensional Chandrasekhar equations

in connection with quadratic control problems [8], [15], [16], [44], and many

of these authors have suggested that such algorithms should lead to very

efficient computational techniques. However, except for a few examples [12],

[15], [16], [40], very few numerical results have appeared in the open

literature that substantiate these claims. Therefore, we have included a

number of numerical examples to illustrate the computational aspects of these

algorithms.

In Section 2 we present an aeroelastic control problem and briefly outline

the derivation of the model. The control problem is then formulated as an

infinite dimensional linear quadratic control problem. We introduce this

model problem in order to provide an example to motivate our work and to check

our computational algorithms. Section 3 is devoted to the development of the

Chandrasekhar equations and computational algorithms for a general distributed

parameter control problem. This section contains the statement of the major

theoretical results. Proofs of these results are given in the Appendix.
-3-

In Section 4 we apply the general ideas developed in Section 3 to control

problems governed by retarded functional differential equations. We

restricted our attention to these problems for two reasons. First, except for

some minor extensions, Gibson's approximation results [25] can be directly

applied to obtain convergence of the Chandrasekhar algorithm. Secondly, we

have considerable numerical experience with this class of problems, which

allows us to compare the efficiency of a number of different numerical

algorithms. We begin Section 4 by summarizing the basic results (see [2],

[3], [25]) concerning the averaging approximation scheme for control problems

governed by retarded functional differential equations. We also consider a

"reduced averaging scheme" that takes advantage of the special structure that

occurs in many problems and which often leads to considerable computational

savings. Finally, we combine the Chandrasekhar algorithm with the reduced

approximation scheme to produce a convergent and computationally efficient

algorithm for approximating gain operators. The infinite time problem is also

discussed.

In Section 5 we present a number of numerical examples to illustrate the

computational aspects of the theoretical results. These examples illustrate

the potential applicability of these factorization methods for very largescale control problems. As noted above proofs will be given in the Appendix.

The notation used in this paper is fairly standard. Given two real

Hilbert spaces X and Y, I(X,Y) shall denote the space of all bounded

linear operators B:X . Y with the usual operator norm IIBII. The

innerproduct and the norm of the space X shall be denoted by <'''>X and

If-fIx, respectively; subscripts will be dropped if it is clear from the

context which space is intended. The adjoint of a closed linear operator
-4-

A:D(A)Cx + Y shall be denoted by A*. The resolvent set of A is denoted

by p(A). Finally, for an interval (a,b) of the real line, L2(a,b;X)

denotes the set of L2-integrable functions Whose values lie in the Hilbert

space X. If the Hilbert space is understood from the context, L2(a,b) shall

be used.

II. AN AEROELASTIC CONTROL PROBLEM

In order to provide the reader with some concrete example and to motivate

our work, we present a brief description of a problem that we shall use to

test the computational algorithms. A more detailed derivation of the model

may be found in [I0].

Consider the problem of controlling the aeroelastlc structure (i.e., the

typical section) shown in Figure 2.1. The airfoil is placed in a flow field

with undisturbed stream velocity U and allowed to plunge and pitch in the

_low.

i ! i I I I i I I 1.1 I1.11 I I I I I_

' h

IJ

Figure 2.I
-5-

Let h(t) denote the plunge and _(t) the pitch of the airfoil at time

t. The equations of motion can be written in the form

(2.1) M y(t) + K y(t) = F(t), S S

where y(t) = (h(t), _(t)) T and F contains the aerodynamic and applied

loads on the airfoil. In particular

= [L(t)+ u(t)1

(2.2) F(t) L Ms(t) j,

where L(t) and M (t) are the aerodynamic loads corresponding to the total

wing lift per unit depth and total moment about the 1/4 chord per unit depth,

respectively. For the airfoil considered here, it follows that (see [9],

[I0], [46])

(2.3) M (t) = _pb3[h(t)/2 + 3b_(t)/8 + U_(t)]

and

(2.4) L(t) = _pb2[h(t) + b_(t) + U_(t)] + (2_pUb)D(t),

where p is the density of air, b is the semlchord length and D(t) is the

"Duhamel integral." In particular, D(t) is given by

t

(2.5) D(t) = f _(U (t - s))Q(s)ds O

where _(Ut/b) is the Wagner function (see [9]) and

t ao

(2.6) Q(t) = f [h(s) + b_(s) + U_(s)]ds.

0
-6-

In order to obtain a state space model that is suitable for control, it is

necessary to provide a useful representation for D(t). Clearly, the key to

this problem is the Wagner function 4. It is possible to show that (see

[11])

(2.7) D(t) = Q(t) - g(t)

where g(t) is the output to a hereditary system (i.e., a functional

differential equation) with input _(t). In particular, it follows that

(2.8) _(Ut/b) = 1 - W(t)

and

(2.9) W(t) = CeAt B

is the weighting pattern for a hereditary control system. Almost all

approaches to modeling aeroelastic structures can be reduced to some scheme

for approximating W(t). For example, R. T. Jones (see [30]) used the twoterm exponential function

-81t -B2t

(2.10) W(t) = =i e + =2e

with _i = 0.165, _2 = 0.335, 81 = 0.0455 and 82 = 0.3. If W(t) is

substituted into (2.5), then D(t) is approximated by

(2.11) D(t) = Q(t) - g(t) = Q(t) _ [=i _i (t) + =2 _2 (t)]'

where _l(t), _2(t) satisfy the second-order ordinary differential equation
-7-

d-_ = + Q(t).

[x2(t) j -8 2 [_2(t)J

In particular, g(t) is the output of the system (2.12) with input Q(t) and

the Jones approximation is equivalent to approximating the weighting pattern

W(t) by

i:-8012

(2.13) W(t) = [=1 =2 ]e "

As indicated above, one can show that the model should include hereditary

terms (see [ii]). We shall consider the simplest possible model of this

nature. This shall be accomplished by replacing (2.12) with the delaydifferential equation

(2.14) xd(t) = gl xd(t) + g2 xd(t - r) + Q(t)

and approximate D(t) by

(2.15) Dd(t) = Q(t) - gd(t) = Q(t) - c I xd(t)

with system parameters gl, q2, Cl' and delay r > 0 (normally these

parameters must be identified using a parameter estimation scheme). The

problem of estimating the time delay r has been considered in other papers

[4], [I0]. The delay is clearly dependent on the chord length b and the
-8-

undisturbed stream velocity U. We have found that for problems considered

below a reasonable asymptotic "first estimate" of r is r = kb where k =

U

I0. In this model, xd(t) represents a generalized aerodynamic "lag state."

To complete the model we augment equations (2.1) - (2.4) with equation (2.14)

and replace D(t) with the approximation Dd(t) defined by (2.15). Let

x(t) denote the flve-dimensional vector

(2.16) x(t) = (_(t), _(t), h(t), _(t), xd(t)) T.

The basic model becomes

(2.17) Ex(t) = F0 x(t) + F1 x(t - r) + Gu(t)

with initial data

(2.18) x(0) = qE_; xd(s) = ¢(s)€L2(-r,0;_

and output

(2.19) y(t) = Cx(t) €I_P.

The control problem we consider is to find u :[0,T] + _ that minimizes

T

(2.20) J(u) = f [fly(s)1, 2 + Ru2(s)]ds

0

where R > 0 and y(t) is the output to the delay-differential system (2.17)

- (2.19).
-9-

It is important to note that only the aerodynamic lag state xd(t)

involves a time delay (i.e., the first four columns of F I are zero). This

observation will be important when we consider the computational algorithms in

Section 5 below. The matrices E, FO, F I, G, and C are given in the

technical report [I0]. Parameters needed to complete the modeling process can

be obtained by applying the parameter identification methods presented in [I0]

to experimental wind tunnel data. A specific example of this process is given

in [I0].

There are a number of state space formulations for this problem. The

approach we take here is a slight variation of the "standard" formulation

given in references [9], [24], [46]. Let A 0 = E-I F0, A I = E-I FI,

-I

B = E G and note that A I has the form

where A d1 is a 5 x I (i.e., column) vector. We choose H = _ x L2(-r,0;l_)

as our state space and let (n,_(')) = (_I' _2' n3' n4' n5' _(')) denote a

typical element in H. Define the linear operator A on H

(2.22) D(A) = {(n, _('))l_(')€wl'2(-r,0;IR), _(0) = _5}

by

d _(-r), $(.)) (2.23) A(_, _(.)) = (A0 _ + A I •

Moreover, we define B:R+ H and V:H + _P by

(2.24) Bu = (Bu,8)
-I0-

and

(2.25) V(n,_(.)) = Cq,

respectively. The delay-differential system (2.17) - (2.19) can be realized

(see [I0]) as the system

(2.26) z(t) = Az(t) + Bu(t)

(2.27) z(0) = z0 = (n,_(-))

(2.28) y(t) = Vz(t).

Moreover, the optimal control problem is equivalent to finding u :[0,T] .

that minimizes

T

(2.29) J(u) = f [lly(s)ll2+ Ru2(s)]ds,

0

where y(t) is the output to (2.26) - (2.28).

The problem defined by (2.26) - (2.29) will be used to test the numerical

schemes described in Section 4. However, it also serves to motivate the

theoretical developments presented below.

III. THE CL_NDRASEKIL_RALGORITHM

In this section our attention is focused on the time-invariant infinite

dimensional linear regulator problem. For a study of control problems

governed by a general evolution equation we refer the interested reader to

Curtain and Pritchard [18], Gibson [26], and Datko [19]. The evolution
-IIprocesses usually arise from control of systems governed by partial

differential equations (see Lions [35] and Lukes and Russell [36]), or

functional differential equations (see Delfour and Mitter [23] and Manitius

[38]). The approach we follow is similar to the development of the optimal

control problem in Gibson [25]. The importance and potential usefulness of

Chandrasekhar equations have been known for some time. There are a number of

papers that discuss the application of these equations to finite dimensional

linear quadratic optimal control problems (i.e., see [13], [14], [31],

[22]). However, the derivation of these equations relied upon being able to

twice differentiate the solution to a matrix Riccati differential equation.

In an infinite dimensional setting there are a number of questions concerning

the existence of these derivations (in a strong sense) that limit the

usefulness of this approach for infinite dimensional systems. The derivation

of the infinite dimensional equations in [45] by Tung is purely formal. The

differentiations used in his derivation were not justified. In [8] and [15]

the authors used a Lions-type setting and derived a set of Chandrasekhar

differential equations satisfied in a distributional sense. We intend to give

an alternate derivation which will lead to a set of Chandrasekhar integral

equations that have unique strong solutions. Moreover, under fairly weak

assumptions it can be shown that the solution to these integral equations is

strongly differentiable (a result that can be applied to establish smoothness

properties of solutions to Riccati operator differential equations [27]).

In [44], Sorine established that the gain operator satisfies a set of

Chandrasekhar equations if the underlying semigroup of the system is

analytic. For systems governed by hyperbolic PDE's or differential-delay

equations (such as the aeroelastic system above) the associated semigroups are
-12-

not analytic, and hence Sorlne's results do not apply. In this section we

give a set of Chandrasekhar equations in an integral form for a large class of

optimal control problems governed by a general distributed parameter model.

The results are established by an approximation technique and leads to

computationally feasible methods. Because the proofs are technical and

provide little insight into the structure of the problem, we have placed them

in an appendix.

Let H and U be real Hilbert spaces. Throughout this section

T(t):H + H, t _ 0, will denote a strongly continuous C0-semigrou p of bounded

linear operators with infinitesimal generator A. We shall always assume that

the following basic hypothesis holds.

H O) i) The operators H:U . H, Q:H . H and R:U + U are continuous

linear operators.

ii) The operators Q and R are self-adjoint and non-negative and

R satisfies UR_ > m > O.

ill) The operator Q can be factored into the form Q = V* V where

V:H + A is a bounded linear operator and A is a Hilbert space.

The infinite dimensional linear quadratic (LQ) optimal control problem

is to find a u CL2(0,tf;U) which minimizes

tf

(3.1) J(z,u) = f [<Qz(s),z(s)> + <Ru(s),u(s)>]ds

0

where z(t) is defined by
-13-

t

(3.2) z(t) = T(t)z + f ST(t - s) u(s)do

0

for 0 _ t _ tf and zcH. Although our discussion is mainly restricted to

the control problem on finite intervals, we will have occasion to discuss its

relationship to the infinite interval control problem. In this case the cost

functional (3.1) becomes

e_

(3.3) J (z,u) = f [<Qz(s),z(s)> + <Ru(s),u(s)>]ds.

0

Under the assumptions given above on Q and R, it is known that (see

[18], [25] for example) there exists a unique u CL2(0,tf;U) , (respectively

L2(O,_;U) for each _ > O) which minimizes J(z,u) (respectively

J (z0,u)). The following characterization of this optimal control may be

found in Gibson [25]. This characterization is discussed in detail in order

to develop the notation we shall need for the derivation of the Chandrasekhar

equations.

For 0 _ s _ tf define the Hilbert spaces Hs and U s by

Hs = L2(s'tf;H) and Us = L2(s'tf;U)' respectively. Let the operators

TsCL(H, Hs), TsCL(Hs,Hs), and FsCL(H,H) be defined by

(Ts z)(t) = T(t - s)z, zcH,

t

(Ts €)(t)= f T(t-n)€(n)dn,0isit_< if, €€Hs, S

and

_ (Ts _)(tf),
-14-

respectively. Straightforward calculations imply that T is given by S

tf

(7"s _)(t) = [ T (q - t)_(q)dn, _EH s

t

and that F has the form

S

(F z)(t) = T (tf - t)z, zEH.

These representations will be useful in the proofs of Lemmas 3.1 and 3.2 in

the Appendix.

For the optimal control problem on the interval [s,tf] with 0 _ s _ tf,

define the cost functional

tf

J(s,z(s),u) = + f (<Qz(_),z(n)> + <Ru(_),u(n)>)dn,

s

where z(t) is defined by (3.2). In [25] Gibson has shown that the optimal

control may be expressed in the form

(3.4) u (t) = - (s), 0 <__s <__t,

^ ^_

for almost all t in [0,tf] where Rs and Bs are defined by

^ _

(3.5) Rs = R + B Ts QT s

and

(3.6) Bs = B Ts QT s,
-15-

respectively. Note that RsEL(Us'Us) and B*_L(H,Us).s Since T(t)z =

T(t - s)z*(s) it follows directly from (3.2) and (3.4) that the optimal

trajectory z*(t) has the representation

(3.7) z (t) = T(t - s)z (s) - f T(t - _)B( (n)z (s)d_,

S

for 0 J s _ t _ tf. If S(t,s)EL(H,H) is defined by

t

(3.8) S(t,s)z : T(t - s)z - f T(t - _)B I Bs)(n)zdn,

s

then (3.7) becomes

z (t) = S(t,s)z, 0 _ s _< t _< tf.

In [25], Gibson shows that S(t,s) is the bounded perturbation of T(t)

by -BR -I B* H(t) where H(t) is defined by

tf ,

(3.9) H(t)z = f T (_ - t)QS(n,t)zdn, zcH,

t

Furthermore, he shows that the optimal control has the representation

* -i * * *

u (t) = -R B H(t)z (t) = -K(t)z (t).

It has also been shown (see [18] and [25]) that _(t) is the unique selfadjoint solution to the two Riccati-type integral equations

tf

(3.10) _(t)z = f T*(s - t)[Q - _(s)BR -1B H(s)]T(s - t)zds

t
-16-

and

tf , ,

(3.11) K(t)z = f S (s,t)[q + _(s)BR -I B _(s)]S(s,t)zds.

t

A formal differentiation of either of these two equations yields a Riccati

differential equation analogous to the case for the finite dimensional linear

regulator problem. Such formal differentiations will be avoided in our

derivation of the infinite dimensional Chandrasekhar equations.

As in the finite dimensional case, under appropriate stabilizability and

detectability requirements, the solution of the infinite time problem may be

viewed as the limit of the solution to the finite time problem as tf + _.

Our only interest in this problem is to show that the Chandrasekhar equations

may also lead to a computationally feasible method for computing solutions to

these problems. Numerical examples will be given in Section 5 to illustrate

this idea and to discuss some of its limitations.

Since our main objective is the development of computational algorithms

for numerically calculating the feedback gain operator for the optimal control

problem (3.1) - (3.2), we shall approach the development of the Chandrasekhar

equations via approximation theory. This theory is based on Gibson's work

[25], [26], and our presentation makes use of a slight extension of his

results. Again, we want to emphasize that these theoretical extensions are

minor and that the real advantage of our approach lies in the computational

savings that come from employing the problem structure and Chandrasekhar

methods. However, this extension will enable us to show that the gain

operator K(t) for the system (3.1) - (3.2) satisfies, in a strong sense, a

set of Chandrasekhar integral equations. Moreover, examples will be given to

show that approximation of the gain through the Chandrasekhar equations can be
-17-

computationally more efficient than the usual method of solving a Riccatl

equation.

{TN(t)}t>O be a sequence of Co-semigroups on H with corresponding

Let

infinitesimal generators _N}. We also assume that BN_L(U,H) and QNEL(H,H)

define sequences of operators with each QN self-adjoint and non-negatlve.

The following hypothesis is needed.

H 1 ) i) There exist constants tO, FI, r2 such that

lITN(t)II< r0

for all N and tE[O,tf] and

"sN"! q, ,,QNH ! r2

for all N.

ii) The semigroups TN(t) and [TN(t)] * converge strongly to T(t)

and T*(t), respectively and the convergence is uniform in t on

[O,tf].

iii) The operators BN, [BN]* and QN converge strongly to B, B*

and Q, respectively.

It is important to note that condition H I - (ii) requires convergence of

the semigroups and their adjoints. This is a critical assumption if one is

interested in applying these results to non-self-adjoint problems.

The Nth approximate LQ optimal control problem is to find [uN]*(t)

which minimizes
-18-

tf

(3.12) jN(zN' uN) = f [<QN zN(s),zN(s)> + <RuN(s),uN(s)>]ds

0

where zN(t) is defined by

N t

(3.13) z (t) = _(t)z(0) + f TN(t - s)_ uN(s)ds

0

for 0 < t < t . Under the above assumptions, Gibson [25] has shown that the

optimal control [uN] * * converges in L2(0,tf;U) to the optimal control u

for the original problem. In particular, he proved that the Riccati operators

HN(t) converge strongly to H(t), and uniformly in t on compact t-intervals

This yields strong convergence of the gains KN(t) m R-I BN HN(t) to

K(t) E R-I B_(t). If the control space is finite dimensional, the convergence

of KN(t) is uniform in norm.

In order to establish analogous results for Chandrasekhar equations, one

needs the following technical lemmas, the proofs of which appear in the

Appendix.

LEMI_ 3.1: If H0 and HI are satisfied_ then IIBN[RN]-I(t)II is -- _ S L S J --

uniformly bounded in N, s, and t for 0 d s d t d tf.

Gibson [25] established that conditions H0 and H I are sufficient to

ensure the strong convergence of sN(t,s) to S(t,s). We make use of a

slightly more general result.
-19-

LEMMA 3.2: If H 0 an__dH I hold, then for each zcH, 0 < s < t < tf

sN(t,s)z . S(t,s)z

and

[sN(t,S)]* z + S (t,s)z,

where the convergence is uniform in s an___dt. In addition, sN(t,s) and

[sN(t,s)] * are uniformly bounded in N, s, t f____r0 d s d t _ tf.

In order to derive an infinite dimensional version of the Chandrasekhar

equations, it is necessary to make an additional assumption on the

approximating sequences. This condition essentially requires that each of the

approximating LQ control problems (3.12) - (3.13) have optimal gain operators

that satisfy a form of the Chandrasekhar differential equations.

H2) i) There exists a sequence of approximating LQ optimal control

problems (3.12)- (3.13)satisfying H I.

ii) The operators QN can be factoredinto QN = [vN]* VN with

VN + V and [vN]* + V* strongly.

iii) The optimal gain operator KN(t) for the problem (3.12) - (3.13)

is stronglycontinuouslydifferentiableand satisfies

(3.14) KN(t)z = -R-I BN*[LN(t)] * LN(t), KN(tf) = 0

(3.15) _N(t)z = -LN(t)_ N - B N KN(t)]z, LN(tf) = VN
-20-

for zEH, where A N is the infintesimal generator of TN(t) and

LN(t) = VN sN(tf,t).

Condition H2) - (ill) is always satisfied if each AN is a bounded

linear operator. In this case the semigroups TN(t) are differentiable, and

the derivation of (3.14) - (3.15) proceeds as in the finite dimensional

case. It should also be noted that conditions H I) - (iii) and H 2) - (li)

are independent in that strong convergence of QN to Q does not imply H2)

- (iii) and conversely that H2) - (iii) is not sufficient to imply strong

convergence of QN to Q.

We may now state the fundamental existence result.

THEOREM 3.3: Suppose that conditions H0, HI and _ hold. If K(t)

denotes the optimal gain operator for problem (3.1) - (3.2)_ then K(t)

satisfies the system of Chandrasekhar integral equations

tf

(3.16) KCt)z = f R-I B* L (s)LCs)zds

t

tf

(3.17) L(t)z = VT(tf - t)z - f L(s)BK(s)T(s - t)zds

t

for all zcH and0 ! t ! tf, and L(t)z = VS(tf,t)z. Moreover, the

approximate g@in operators KN(t) converge strongly to K(t). If U is

finite dimensional, then

(3.19) lim lIKN(t) - K(t)ll = 0

N+_

and the convergence is uniform in t o__n [0,tf].
-21-

Theorem 3.3 does not imply that K(t) and L(t) are the only strongly

continuous solutions to (3.16) - (3.17). The proof that these equations have

unique strongly continuous solutions is non-trivlal, but essential if we hope

to use these equations as a basis for computational algorithms.

THEOREM 3.4: Assume that H0, H 1 and H 2 hold. Then K(t) an___dL(t) =

VS(tf_t) are the unique strongly continuous solutions to the Chandrasekhar

equations (3.16) - (3.17).

Some comments concerning the assumptions H 0, HI and H 2 are in

order. Although these conditions may at first glance seem severe, they are in

fact the properties one would llke to have in a scheme that is to be used to

numerically approximate (3.1) - (3.2). In this case the convergence

properties are assured and the Chandrasekhar integral equations may be

approximated directly in order to obtain the gain operator. The averaging

approximation scheme which is discussed in Section 4 satisfies these

conditions when applied to linear regulator problems governed by delaydifferential equations. For distributed parameter systems, Lukes and Russell

[36] state conditions which also satisfy these hypotheses. Their

approximations are essentially eigenfunction expansions and apply to problems

whose infinitesimal generators have normal extensions (i.e., heat and wave

equations with the appropriate boundary conditions).

While our major concern is to illustrate that approximation of the gain

operator through the Chandrasekhar equations may have computational advantages

over approximation of the Riccati operator, it should be pointed out that the

theory developed is not restricted to finite dimensional approximations. This
-22-

led Kazufumi Ito (see [27]) to note that if condition H0 was satisfied, then

the Yosida approximations could be used to define approximating sequences that

satisfy H 1 and H2. In particular, let A N = NA(NI - A) -1 (for N

sufficiently large so that NCp(A)) and define BN = B, QN = Q, and VN = V

for each N. It can be shown that these approximating sequences satisfy H 1

and H 2 and hence one has the following theorem on the existence and

uniqueness of Chandrasekhar integral equations.

THEOREM 3.5: If H0 is satisfied, then K(t) and L(t) = VS(tf,t) ar__._ee

the unique strongly continuous solutions to the Chandrasekhar equations (3.16)

- (3.17).

Although the Yosida approximates provide a set of approximating sequences

that satisfy H I and H2, in general these approximations are not very useful

for numerical schemes. For example, the A N operators are infinite

dimensional and not easily constructed. It is for this reason that we have

stated the conditions that are needed to not only obtain existence and

uniqueness but to also guarantee the convergence of the gain operators.

We turn now to the application of these ideas to a particular class of

distributed parameter systems. In particular, we shall concentrate on control

systems governed by functional differential equations.
-23-

IV HEREDITARYCONTROLPROBLEMS

In this section we consider the control of a linear hereditary system of

the form

- Lxt + B0 u(t), 0 _< t _< tf

(4.1)

x(O)= n, x0 = €,

where n E I_n, ¢€L2(-r,0;l_n), uEL2(0,tl;_ m) for each tl < _, B€i(_m, _n),

and xt = x(t + 8) for -r J O J 0. The linear operator L has the form

0

Lx t = _ A i x(t - h i ) + f D(O)x(t + O)dO, i=0 -r

where each AiEi(_n,_n) , 0 = h0 < h I <...< h v = r, and DCL2(-r,0;i(_n,l_n)).

Let Q be a real, symmetric, nonnegatlve n × n matrix, and R a real,

symmetric, positive m x m matrix. The optimal control problem is to find

u _(t)EL2(0,tf;l_m) _ _ which minimizes

tf

(4.2) J(u) f [<Qx(s),x(s)> - <Ru(s),u(s)>]ds

0

where x(t) is the solution to (4.1) corresponding to u(t).

It is well-known that (see [3] for example) the system (4.1) has a unique

solutionfor each u(t)EL2(0,tf;l_)and (n,_)E_× L2(-r,0;Ka)_ Z. In

addition,the solutiondepends continuouslyon _, $, and u. Under the above

stated conditionson Q and R there existsa unique u EL2(0,tf;_m)

which minimizes(4.2) subject to (4.1).
-24-

In recent years a standard technique for obtaining the solution to (4.1)-

(4.2) has been to reformulate the system as an evolution equation on the

Hilbert space Z. The solution to problem (4.1) - (4.2) is then approximated

via the abstract formulation. Taking this approach, we will see that the

dynamical system falls into the framework described in the previous sections

and the optimal gain operator is thus characterized via the Chandrasekhar

equations (3.16) - (3.17).

For (n,_(.))EZ,define the operator T(t)Ei(H,H) by

(4.3) T(t)(_,#(.)) = (x(t),xt(.))

where x(t) is the solution to the homogeneous version of (4.1), i.e.,

u(t) m 0. It is well-known (see [3]) that T(t) is a C0-semigrou p on Z

with infinitesimal generator A characterized by

(4.4) _(A) = {(_,_(.))l!(.)Ewl'2(-r,0;]Rn), n = _(0)}

and

(4.5) =

for (_,_(-))E_(A). Moreover, the delay differential equation (4.1) is

equivalent to the abstract evolution equation (AEE)

t

(4.6) z(t) = T(t)(n,#(.)) + f T(t - s)Bu(s)ds

0

where Bu = (B0 u,0). In particular, the following result may be found in

[4].
-25-

THEOREM 4.1: Let (n,_)EZ be given. If x(t;u) is the solution of

(4.1) for uEL2(0,tf;_m) , then z(t) defined by (4.6) satisfies z(t) =

(x(t;u),xt(.;u)),t _ 0.

Define the operator Q€i(Z,Z) by

where 0 denotes the zero operator on the appropriate spaces. The control

problem (4.1) - (4.2) is equivalent to finding u (-)€L2(0,tf;_Im) which

minimizes

tf

(4.7) ] = f [<Qz(s),z(s)> + <Ru(s),u(s)>]ds

0

where z(t) is defined by (4.6). Since the problem defined by (4.6) - (4.7)

falls into the framework above, the optimal control has the representation

(4.8) u*(t) = -R-I B* E(t)z*(t) = -K(t)z*(t)

where E(t) satisfies the Riccati integral equations (3.10) and (3.11) and

z*(t) is the optimal trajectory. Due to the special structure of the state

space Z = _ x L2(-r,O;_n), K(t) may be expressed as the matrix

of operators

(4.9) _(t) = , 0 _ t ! tf

t_10(t) _ll(t)J
-26-

where there exists a real, nonnegative, symmetric n x n matrix H00(t), and

^ A

_00(t)n = _00(t)n, _10(t) is a bounded linear operator from L2(-r,O;lln) ,

and can be represented by

A

[_10(t)_](s) = Nl0(t,s)n

A

where Hl0(t,s) is an n x n matrix function, _01(t) is the bounded linear

operator from L2(-r,0;_n ) into _n given by

0

(4.10) H01 (tIC = f _I0T (t;s)¢(s)ds, ¢€L2(-r,0;Rn )

-r

and ill(t) is a real, non-negative, self-adjoint operator on L2(-r,0;_n ).

Further properties of these operators may be found in [25], [33], and [35].

Combining (4.8) - (4.10), it follows that the optimal control may be written

as

, _R_I 0

-r

for 0 J t J tf. Therefore, the gain operator has the form

0

(4.12) K(t)(_,_(.)) = K0(t)n + f Kl(t,s)¢(s)ds

-r

where K0(t) = R-I BO I H00(t) is an m x n matrix and

Kl(t,s ) = R-I BT0 _10(t,s) T is an m x n matrix valued function. For each

t _ 0 the function s . Kl(t,s) is called the functional gain.

It is important to point out that the Riccati operator maps an infinite

dimensional space to itself while the gain operator maps an infinite
-27-

dimensional space to a finite dimensional space. Note that the gain is the

actual operator needed to produce the optimal control. Thus it should be

advantageous to approximate K(t) instead of _(t) since it involves

approximations in one spatial variable only. This difference can be accounted

for by examining equations (4.9) and (4.10). The infinite dimensionality in

^

two spatial variables of K(t) is due to the operators Kll(t) and _lO(t).

A

However, the operator Kll(t) does not appear in the representation (4.11),

A ^

and the information in Kl0(t) is contained in the operator E01(t). The

Chandrasekhar equations exploit this reduced dimensionality, and for the

particular approximation scheme below it will be seen that this reduction can

be quite significant.

We turn now to a particular approximation scheme (the so-called AVE scheme

[3], [4], [25]) for the simple delay differential equation

x(t) = _ x(t) + A1 x(t - r) + B0 u(t)

(4.13)

x(0) = _, x(s) = _(s), -r < s < 0.

This scheme has been extensively studied by Banks and Burns [3] and Gibson

[25]. The treatment of a more general case which includes multiple discrete

delays and a continuously delayed term may be found in [3].

For any positive integer N, partition [-r,0] into subintervals

[t_, tj_l] N for j = 1,2,-.-,N, where tjN = _ jrN" Let X_ denote the J

N tN N

characteristic function on [tj, j_l ) for j = 2,3,...,N, and X1 the

characteristic function on [-r/N,0]. Define the finite dimensional

subspaces zNcz by

N

ZN = {(B,_)EZ I D E]_n, _ = _ vj X_ v. E ]Rn}

j=l 'J
-28-

and the projections PN:z + ZN by

N

N N

pN(+)=(+o +i×j) j=l

whe re

tN

i-I

N N N

40 = n, _i = Y f _(s)ds, for j = 1,2,-..,N.

tN.

J

Let the operators AN:z + ZN be defined by

AN(n'#) = (A0 _0 + AI _N' _ j=l - _j)X ).

If TN(t) denotes the semigroup generated by A N , then the approximation to

the abstract formulation of (4.13) is

t

zN(t) = TN(t)pN(n,_) + f TN(t - o)(B 0 u(a),0)do. 0

Since AN is bounded for each N, this is equivalent to the initial value

problem

zN(t) = AN zN(t) + (B0 u(t),O)

(4.14)

zN(o) = eN(n,¢)

Note that the operator AN is reduced by the Hilbert space ZN for

each N. Since (B0u(t),O)EzN for every N, (4.14) is a differential

equation in a finite dimensional space. Upon choosing the appropriate basis,

equation (4.14) has the representation (see [3])
-29-

wN(t) = AN wN(t) + BN u(t)

(4.15)

w N (0) = col I+NO,_I,...,@N N N) ,

where AN is the n(N + i) x n(N + I) matrix

A 0 0 • • • 0 A 1

N I i N__ •

r r

e

• •

e

(4.16) 0 • • • ,

e

e

• • • 0

0 .... 0 --NI NI

r r

I is the n x n identity matrix, and BN is the n(N + I) x m matrix

BN = col(Bo,O,...,O).

The Nth approximate control problem becomes: find the

N*

[u ] EL2(O,tf;_) which minimizes

tf

(4.17) jN(wN(0),uN) = f I<QN wN(s),wN(s)> + <RuN(s),uN(s)>)ds,

0

where wN(t) is the solution to (4.15). The matrix QN is given by

where 0 denotes the appropriate n × n zero matrix.
-30-

In [3] the authors proved that the AVE scheme satisfies a Trotter-Kato

type approximation theorem. This established the strong convergence of the

semigroups TN(t) to T(t) defined by (4.3). Recall that one of the

requirements of _ is that the adjoint semigroups [TN(t)] * converge

strongly to T*(t). This result was obtained by Gibson [25] and is the key

fact used to prove Lemma 3.2. The special forms of BN and QN immediately

satisfy Ho,HI, and H2; thus the solutions to the Chandrasekhar equations

(3.14) - (3.15) of the AVE approximation converge to the solutions of the

Chandrasekhar integral equations (3.16) - (3.17) in the strong sense.

Before demonstrating that approximation of the gain operator via the

Chandrasekhar equations has significant computational advantages over

approximation of the gain operator via the Riccati equation, we introduce an

additional reduction technique which takes advantage of a special structure

that frequently occurs in hereditary systems. For these special systems, this

technique will further reduce the number of differential equations that it is

necessary to solve in order to approximate the gain operator, K(t).

In applications of hereditary control problems, the delay does not always

appear in each component of the state (see Examples 5.1, 5.2, and 5.3). In

this type of system the initial data actually lie in a "smaller" state

space. The AVE scheme as described projects each of the n components of the

history portion of (x(t),x t) into an N-dimensional space. If only q < n

of the components are delayed, it should be possible to project the history

portion into a space of dimension q × N instead of a space of dimension

n × N. The following discussion develops this idea.

The reduction as discussed here was originally studied in Cliff and Burns

[17] for spline approximations in the context of parameter identification, and
-31-

then for AVE in [I0]. This reduction is a special case of the F-reduction

discussed in [21] and [22].

For the remainder of this section it will be assumed that the matrix A 1

defined in equation (4.13) has the form

A I =

A 1

where All is p x q, AI2 is q × q, and p + q = n. We shall write x(t) =

col(w(t),y(t))€_ p x Rq and not distinguish between columns and rows for ease

of notation. Also, we shall not distinguish between the space Z defined

earlier in this section and _P x _q × L2(-r,0;]RP) x L2(-r,0;]Rq).

Let ZR denote the "reduced space" _P × _q x L2(-r,0;_q). Define the

projection operator P:Z + ZR by

P(nl,n2,#l,_ 2) = (nl,n2,_2),

and the injection I:ZR + Z by

^

I(nl,n2,¢2) = (nl,n2,nl,¢2),

A

where n I is the function in L2(-r,0;_P) with constant value BI" For

A and T(t) as defined by (4.3) - (4.5) define the reduced operator

_:D(_)CZ R . ZR by

(4.18) AR = PAl,

with domain
-32-

(4.19) D(AR) = {(nl'n2'i2)EZR I _2Ewl'2 _2 = _2 (0)}

In [I0], Burns and Cliff stated but did not prove the following theorem. We

include the proof for completeness.

THEOREH _.2: The operator AR is the generator of a C0-semigrou p

{TR(t)}t>0 on ZR, and TR(t ) =PT(t)I.

Proof: Let IR be the identity operator on ZR. It is easy to show

that P and I are bounded and it follows that PT(t)I is strongly

continuous in t. Also, since p is a left inverse for I, we have that

PT(0)I = IR. In order to show that PT(t)I satisfies the semigroup property,

first note that the solution to (4.13) (with u = 0) is independent of the

initial choice _I" Hence, for any _I' _ICL2 (-r'0;]Rp)' the equality

(4.20) T(t)(nl,n2,_l,# 2) = T(t)(_l,n2,_l,_2) , t _ 0

is obtained• Let (w(t),y(t),wt,Yt) denote the solution to (4.13)

corresponding to (nl,n2,_l,_2) and u _ 0 (recall with our notation x(t) =

col(w(t),y(t))). For tI and t2 _ 0,

A

pT(tl)IPT(t2)I(n I,n2,_2 ) = PT(tl)IPT(t2)(n l,n2,n 1,12 )

(4•21) = pT(tl) l(w(t2),y(t2),yt2 )

= pT(tl)(W(t2),Y(t2),w(t2),Yt2).
-33-

In view of (4.20), (4.21) becomes

pT(tl)_PT(t2)I(_l,n2,_ 2) = PT(tl)(W(t2),Y(t2),wt2,Yt2 )

= P(w(t I + t2),Y(t I + t2),Wtl+t 2'ytl+t2 )

= (w(t I + t2),Y(t I + t2),Ytl+t2).

A straightforward calculation yields

pT(t I + t2)_(Bl,n2,_ 2) = (w(t I + t2),Y(t I + t2),Ytl+t2 ),

and the semigroup property is satisfied. Using standard agruments one can

show that the domain of the infinitesimal generator of PT(t)I is given by

(4.19). For each (nl,_2(0),_2)€_(A R) it follows that

PT(t)l - IR

lim t (_I'_2(0)'_2) = (A0(w(0)'y(0)) + Al(W(-r)'y(-r))'Yt)- +

t+_

It is easily verified that

AR(nI,_2(0),_2) = (A0(w(0),y(0)) + Al(W(-r),y(-r)),yt)

for (_I,@2(O),_2)€D(AR); therefore AR is the infinitesimal generator of

pT(t)[ and rR(t) = PT(t)I.
-34-

The "reduced" abstract evolution equation becomes

t

(4.22) ZR(t) = TR(t)(nl,n2,12) + f TR(t - s)_ u(s)ds 0

where B R u = (B0 u,0) (here B0 uE_ n should be written as an element of

I_p x Rq). Theorem 4.2 and the action of the operators P and I yield the

equivalence of (4.6) to (4.22), and hence ZR(t) = (w(t),y(t),y t) with x(t)

= col(w(t),y(t)). The cost functional for the control problem becomes (4.7)

with the zero operators in Q adjusted to the appropriate spaces. As a

result of these observations, we conclude that solving the optimal control

problem associated with (4.20) is equivalent to solving the optimal control

problem of the original system (4.14).

Define the reduced AVE operators A N by A N = PA N I Computations R R "

similar to those in the proof of Theorem 4.2 reveal that the semigroup T_(t)

generated by A_ satisfies T_(t) = pTN(t)_ and furthermore, it is easily

verified that A NR and T_(t) are the same operators that result if the AVE

scheme is applied directly to (4.21). For (nl,n2,_2)EZR, it follows that

T_(t)(_l,n2,_2) =pTN(t)I(_I,n2,_2 )

= p(TN(t)(nl,_2,_l,_2)).

As commented earlier, Banks and Burns [3] proved that TN(t)z . T(t)z

^

uniformly in t on compact intervals for each zEZ; thus TN(t)(_I,_2,_I,#2 )

^

. T(t)(nl,n2,nl,_2). Since P is bounded, T_(t) + TR(t) strongly, and

uniformly in t on compact intervals. These results are summarized in the

following theorem.
-35-

THEOREM 4.3: Let z_(t) denote the solution of the reduced AVE scheme

approximation to (4.22) and suppose that ZR(t) = (w(t),y(t),y t) where x(t)

= col(w(t),y(t)) is the solution to (4.1) with initial data (_i,n2,_2).

Then z_(t) + ZR(t) i___nZR, and the convergence is uniform in t on compact

intervals.

If the AVE scheme is applied to the reduced equation (4.20) and the

appropriate basis chosen, then the system of equations

N wNR(t) + N u(t)

wN(0) = col(nl,_2,(_2)l,...,(_2)N) N N

is obtained. The (n + qN) x (n + qN) matrix _ has the form

All

A0 0 • • • 0

_= AI2

0 Nz -Nz 0 r r

• • •

• • • •

• • 0

0 • • • 0 --NI - --NI

r r
-36-

where I is the q × q identity matrix, 8 is the q × p zero matrix, and

0 represents the zero matrix of appropriate size. The reduced scheme has

effectively discarded the zero columns of A1 that appear in (4.16).

The strength of the Chandrasekhar algorithm lies in a comparison between

the number of differential equations in its approximation to the number of

differential equations in the approximation of the Riccati operator. We now

demonstrate that the Chandrasekhar algorithm drastically reduces the number of

equations that it is necessary to solve in order to approximate the gain

operator. Moreover, if the structure of the hereditary system is such that

the additional reduction technique discussed above applies, then further

computational reductions occur.

Let the rank of the matrix Q in (4.2) be P0" For the AVE scheme, the

approximations QN also have rank P0 and it is possible to obtain a

factorization QN = [_]T VN where VN is a P0 x n(N + i) matrix. Thus

the set of Chandrasekhar equations associated with (4.15) and (4.17) contains

C(N) _ (m + p0)[n(N + i)]

equations. The Riccati differential equation for (4.15), (4.17) has the form

(see Gibson [25])

[_N(t)] = _QN _ [AN]T pN(t ) _ pN(t)AN

(4.23) + pN(t)BN[R]-I[BN]T pN(t )

pN(tf) = O.
-37-

Taking into account the symmetry of the the system, (4.23) contains

n(N + l)[n(N + i) + i] R(N)

2

differential equations. Similar counts for the reduced system yields

CR(N) = (m + p0)(n + qN)

and

RR:_N t = (n + qN)(n 2+ qN + i)

respectively.

The special form of QN yields P0 _ n and thus P0 is independent of

N. Tables (4.1) and (4.2) give a comparison of the sizes of the different

systems for various values of the parameters. The parameters in Table (4.2)

are for the wind tunnel problem of Example 5.2.

Table 4.1: n = 2, q = i, m = i, P0 = 2.

N R(N) C(N) RR(N) CR(N)

16 595 102 171 54

32 2775 198 595 102

64 8515 390 2211 198
-38-

Table 4.2: n = 3, q = I, m = i, P0 = i.

N R(N) C(N) RR(N ) CR(N)

16 1326 102 190 38

32 4550 198 630 70

64 19,110 390 2278 134

Based on an equations count only, for sufficiently large N (which is

small for most problems) the Chandrasekhar equations should offer substantial

savings in time necessary to compute the gain. A comparison of R(N) to

C(N) shows a reduction in the number of equations from 0(N 2) to 0(N). If

in addition the F-reduction technique applies, then the number of equations

is further reduced by a factor approximately equal to q/n.

A few remarks concerning generalizations of the problem discussed above

are worthwhile. The inclusion of multiple discrete delays in equation (4.13)

does not reduce the effectiveness of the Chandrasekhar equations. These

additions affect the size and form of AN given by (4.16) (see [3]). In

particular, for 2 delays AN becomes an n[mN + I] x n[mN + i] matrix, and

the equation counts above are correct if C(N) is replaced by C(2N), R(N)

by R(2N), etc. The presence of a continuous delay term does not affect these

equation counts; it only affects the form of AN . A numerical example

involving two delays is given in Section V.

Another generalization would be to include a term of the form <Gx(tf),

x(tf)> in the cost functional (4.2), which penalizes the final state. The

addition of this term is motivated by the fact that in the finite dimensional
-39-

case the inclusion of such a term has virtually no effect on the derivation of

the Chandrasekhar equations other than to change the final values of K(t)

and L(t) (see [14], [31]). In this case, the Chandrasekhar equations for

the AVE approximation scheme may be written

_N(t ) = R-I[BN]T[LN(t)]T cLN(t)

KN(tf) = R-I[BN] T GN

_N(t) = -LN(t)[A N - BN KN(t)]

N

LN(tf) = B 1

where

N and C are matrices satisfying 0 being the appropriate zero matrix and B 1

We point out that in the case where only discrete delays appear, the special

structures of AN , BN, GN, and QN allow us to find C independent of N. In

this case the matrices BN1 have a constant rank P0 with PO < (k + l)n

where k is the number of discrete delays. In the case of a continuous

delay, it follows that for each N, P0 = P0(N) and C = CN are functions

of N. However, it still follows that Po(N) _ (k + l)n for all N.
-40-

No claim is made concerning the convergence of these equations to an

infinite dimensional form, but we do point out that the equations appear to

converge numerically (see Examples 5.1 and 5.5), and it is interesting to note

that an infinite dimensional version for this case has been derived in [27].

The use of the Chandrasekhar equations is not restricted to the AVE

scheme. They may be applied to other schemes as well. Approximations using

splines have become increasingly popular ([4], [6], and [7]). The spline

scheme introduced by Banks and Kappel [6] does not satisfy the the strong

convergence criteria on the adjoint semigroups that our work requires.

However, a modification of this scheme by Kappel and Salamon [30] does. In

fact, Kappel and Salamon have now applied a Chandrasekhar algorithm to their

numerical examples and obtained satisfactory results. Another scheme for

which the Chandrasekhar equations are applicable is the Legendre-tau

approximation scheme introduced by Ito and Teglas in [28] and later applied to

hereditary control systems in [29]. This scheme satisfies the strong

convergence criteria on the adjoint semigroups and retains the low rank

condition on QN which makes the implementation of the Chandrasekhar

equations effective.

V. NUMERICAL EXAMPLES

In this section we present several examples which illustrate the numerical

efficiency of the Chandrasekhar equations and the F-reduction technique when

applied to delay-differential systems. These examples are presented to show

the efficacy of the methods and are not intended to be a complete numerical

test. The "applications" in this section are twofold. First, the two
-41-

computational reduction techniques are applied to the optimal control problem

on a finite interval. Secondly, since the infinite interval problem is the

"limiting solution" of the finite interval problem under appropriate

conditions, the Chandrasekhar approximations are integrated backward in order

,to obtain a "steady-state" value. This steady-state value is then used in the

forward integration of the state equations, and the optimal control is thus

approximated.

Each of the examples presented here was numerically solved on several

different computers. In order to make a consistent comparison between

examples, all of the computer run times recorded in this section will refer to

the CPU time required on an IBM 3081 (located at VPI & SU). However, the

plotted data given here may be from different machines, specifically, either a

VAX 11/750 or the Cyber VPS-32 (Cyber 205 with enhanced memory located at NASA

Langley Research Center) running in scalar mode. A standard fourth-order

Runge-Kutta scheme was implemented to solve all differential equations, and

unless stated otherwise, a fixed stepsize of h = .01 was chosen. For the

finite interval optimal control problems, the Chandrasekhar equations were

integrated backward from tf and the gain values stored in increments of

10*h. In the forward integration of the state equations, linear interpolation

was used to obtain the intermediate gain values.

Example 5.1: This problem has been numerically solved in [3] and [5] and

an analytical solution given in [3]. The optimal control problem is to

minimize

2

J(u) = 10x_(2) + f u2(s)ds 0

where
-42-

x2(t) j x2(t)j - [x2(t 1

Xl(S - i for -I _< s _< 0.

[x2(s)j

The equation is the vector formulation of a harmonic oscillator with delayed

damping given by

y(t) + y(t - I) + y(t) = u(t).

The F-reduction technique combined with the Chandrasekhar algorithm was used

to compute the optimal controls for N = 8, 16, 48, and i00. The total CPU

times required for these computations are given in Table 5.1. Of particular

note are the results for N = 48.

Table 5.1

N 4 8 16 48 I00

CPU(sec) .62 .94 1.57 4.18 8.3

In [3], the authors state that the CPU time required to approximate the

optimal control by solving the Riccati equations for the AVE scheme with N =

48 was 3700 seconds. These computations were performed on an IBM 370/158

system at VPI & SU in 1974. The IBM 3081 on which the present results were

obtained is approximately 6 - 7 times "faster" than the IBM 370/158 (this is a

rough estimate given by the computer center at VPI & SU). The CPU time in [3]
-43-

also included an additional computation to obtain j48 which the present

results do not. The total CPU time for our computations was only 4.18

seconds, and even after taking into account the differences between the two

computers, the reductions are still quite substantial. The reduction in CPU

time may be credited to the reduction in the number of differential equations

solved. In [3], the authors obtained the gain via a Riccati equation and thus

had to solve 4851 differential equations in addition to 98 state equations.

Using the techniques discussed in this paper, it is necessary to solve only

150 differential equations to obtain the gain and 50 additional equations to

obtain the optimal trajectories.

Figures 5.1 - 5.3 illustrate the convergence of the approximation

IN,t)]* to the optimal control u*(t). Note that in Figure 5.3 the

approximation [ul00(t)] * follows the true solution closely until a time

value of approximately i. and then appears to lose accuracy. Implementing a

smaller stepsize in the Runge-Kutta scheme corrected this.

Also included here are graphs of the approximations to the functional

gain a . Kl(t,a), of the gain operator K(t). Recall that the averaging

scheme produces piecewise constant approximations to functions in L2(-r,0) ;

thus the elements O and A in Figures 5.4 - 5.11 represent these constant

values on the appropriate subintervals [ai_l,ai] of [-r,0]. In order to

obtain illustrative clarity, these points were connected to form a smooth

curve for N = i00. Pictured are approximations of the functional gains

K_(t,a), for N = 8, 32, I00, at the values of time t = 0, .25, .5, .75, I.,

1.25, 1.5, and 1.75. These figures show not only the convergence in N of

the functional gains, but also the evolution in time of these approximations.

Note that the approximate functional gains calculated at t = 1.25, 1.5, and
-44-

1.75 do not clearly indicate the discontinuities that actually appear in the

functional gain due to the presence of a terminal cost term in the functional

to be minimized (see Delfour [20]). This phenomenon will be discussed in more

detail in Example 5.5.

Finally, graphs of the resulting state approximations for N = 8 and N =

I00 which show rapid convergence of the states are given in Figures 5.12 and

5.13.

Example 5.2: This next example is an application of the Chandrasekhar

equations and of the F-reduction technique to an infinite time optimal

control problem. This illustrates that the Chandrasekhar equations may be a

viable method of computing the constant gain associated with the infinite time

problem. The motivation for studying this particular example is not only that

the problem has the special structure which fits our framework, but that it is

a practical problem as well. The system is a model for fine tuning the mach

number in a cryogenic windtunnel (National Transonic Facility) constructed by

NASA at its Langley Research Center in Hampton, Virginia. This model has been

studied in [37], [39], and the LQR problem numerically approximated in [7],

[33].

The controller of the system is an actuator attached to a wind guide vane

and finely tunes the mach number by changing the angle of the vane. The state

consists of the variation of the mach number, the variation in the guide vane

angle velocity, and the variation of the guide vane angle. The equation has

the form
-45-

x(t) = -2bm - x(t) + 0 x(t - .33) + u(t),

1 0

x(s) = coi(-.I,0.,8.547), -.33 ! s ! 0,

2

where the parameters i/a, m , k, and b have the values 1.964, 36., -.0117,

and .8, respectively. The quantity to be minimized is

J(u) = f (xT(s)Qx(s) + u2(s))ds,

0

where Q = diag(104,0,0).

In [7], the authors were comparing the AVE approximation scheme to a

spline approximation scheme. They applied a Newton iteration scheme to the

Riccati matrix equations that resulted when the AVE scheme was used (for N =

2, 4, and 8). However, the Newton iteration scheme did not converge to the

solution of the Riccati equation for N = 8.

In our approach, we solved the finite time optimal control problem on an

interval [0,tf] of sufficient length so that the gain KN(t) satisfied

N

KN(t) . K0 as t . 0+, (i.e., we integrated backward to a steady-state

K_). This method resulted in a convergent scheme for N = 8, and in

solution

fact, for larger values of N as well.

48(0,a) Figure 5.14 illustrates the convergence of the functional gain K 1

to a steady-state value for the values of tf equal to 2, 3, 5, and I0. Note

the rapid convergence of the gain to a steady-state value and that the values

for tf = 5 and I0 are almost indistinguishable.
-46-

The states and optimal control computed by using the steady-state gain,

K48(0) are given in Figures 5.15 - 5.18. As can be seen, the response of the

system is very good, and the states are driven to zero rapidly. These results

agree well (pictorially) with those found in [7].

The total computer time required to compute the gain, optimal control, and

states for this problem (N = 48) was 7.76 seconds. This computation used

both the Chandrasekhar equations and the additional reduction technique

outlined earlier. The same problem was also solved on the computer using only

the Chandrasekhar equations and not the additional reduction. In this case

the total CPU time was 21.55 seconds, a significant increase.

Example 5.3: In this example the Chandrasekhar algorithm and F-reduction

techniques are applied to the hereditary model of the two-dimensional airfoil

discussed in Section 2. The resulting delay model has the form (see equations

(2.16) and (2.17))

x(t) = A0 x(t) + A I x(t - s) + Bu(t)

(5.1)

x(s) = ¢(s) -r < s < 0

with output

(5.2) y(t) = Cx(t).

Here, equation (2.17) has been multiplied by the inverse of E. The matrices

in (5.1) - (5.2) are
-47-

--7 .3865 -33 .8517 -4479 .50 -5711 .27 3 .0804--

.9378 -12.4685 568.73 -7068.77 -.3911

A 0 = 1.0 0.0 0.0 0.0 0.0

0.0 1.0 0.0 0.0 0.0

-2.7443 1404.29 -1664.26 -40701.67 -253.856 m

AI(5,5) = -47.00, Al(i,j) = 0, (i,j) € (5,5)

B = [-82.333 191.286 0.0 0.0 864.533] T

C = diag[_50.00 _50.0 _i0.0 _I0.0 i]

i(s) = coi[-.80 .50 .055 .029 50.0], s€[-r,0],

and the time delay is r = .05. The cost functional that we desire to

minimize is

tf

J(u) = f (yT(t)y(t) + Ru2(t))dt

0

where the final time is tf = .25 sec., and R = I0.

Because of the small delay, the stepsize h used in the Runge-Kutta

scheme was reduced to .001. The computations for this particular example

exhibited in Figures 5.19 - 5.24 were performed on a VAX 11/750 computer using

single precision arithmetic.

The parameters used to construct the above matrices are a slight variation

of those found in [i0]. In particular, the non-dimensionalized distance, x ,
-48-

from the airfoil center of mass to the elastic axis (e.a.) is .5 which places

the e.a. at the quarter-chord point. In addition, the free-stream velocity

U was chosen to be 1500 in/sec. This particular value was selected since it

produced unbounded oscillations in our open loop system (a few numerical tests

determined a "flutter speed" of approximately 1375 in/sec for this model).

This unstable plant models the airfoil operating under flutter conditions,

i.e., the pitch and plunge modes couple in a manner that allows unstable

oscillations to occur. Since it is desirable to eliminate flutter, the

controller should be able to damp effectively the pitching and plunging

motions.

The reduced AVE scheme with N = 16 was implemented in order to

approximate (5.1) - (5.2). The time-varying gain values produced by the

Chandrasekhar equations were used to obtain the optimal control, ul6(t).

Graphs plotting the closed loop responses (dotted graph) of the LQR design are

compared with the responses of the open loop (solid graph) system in Figures

5.19 - 5.23. As can be seen in the figures, the pitching and plunging

motions (x3 and x4 respectively) and their velocities (xI and x2) are

quickly driven to zero in the closed loop response, illustrating the

effectiveness of an LQR approach for our model. The optimal control which

produced these responses is given in Figure 5.24.

The computational advantages of the reduction techniques discussed in this

paper are very apparent in this example. For N = 16, in order to solve for

the gain via the Riccati equation, one must solve 3,655 differential equations

(this includes taking advantage of the symmetry of the equation); moreover,

there are 85 state equations to solve. In contrast, if the two techniques

discussed here are employed, it is necessary to solve only 126 equations to
-49-

obtain the gain and 21 equations to obtain the state. If N is further

increased to 48, then the number of equations becomes 30,135 Riccati and 245

state equations compared to 318 Chandrasekhar and 53 state equations. A

similar problem (i.e., same system with different aerodynamic parameters) was

solved on an IBM 3081. For N = 16, the total CPU time required by the

methods here was 8.47 seconds. When N was increased to 48, the CPU time

increased to 20.13 seconds.

Example 5.4: This example is a variation on Example 5.3. The parameters

x and U defined previously are set equal to .297 and 1325.,

respectively. This value of x moves the e.a. from the quarter chord point

(note also that from the figures U is below the "flutter speed"). The

resulting matrices for equations (5.1) - (5.2) may be found in [12]; all other

values are as in Example 5.3. The purpose of this example is to compare the

state response obtained by using the steady-state value of the gain to the

state response obtained by using the time-varying gain. In applications, use

of a constant gain is often preferable since the gain values do not have to be

continually updated at each instant of time.

The Chandrasekhar equations were integrated on the interval [0,.25], and

the value of KN(t) at t = 0 was used as the steady-state value. Figures

5.25 - 5.29 compare the open loop response to the closed loop response on

[0,.5] obtained by integrating the state equations using the constant gain.

The closed loop responses produced by the tlme-varying gain yielded results

that appeared graphically the same as the responses produced by the constant

gain. There were, however, slight numerical differences. Graphs of the

absolute values of these differences for _(t) and _(t) are given in
-50-

Figures 5.30 and 5.31. These graphs are representative of the results for all

of the states. Figure 5.32 is included to show how the kernel portion of the

48

functional gain K 1 (t,_) evolves in time. In [i], the authors treated a

similar problem. However, their results showed high frequency oscillations _

near t = .25, which may be due to numerical instabilties.

Example 5.5: This last example is a one-dimensional equation with two

delays considered in [20]. The optimal control problem is to minimize

2

J = i0x2(2) + f u2(s)ds

0

where x(t) satisfies

x(t) = x(t - .5) + x(t - i.) + u(t)

x(s) - I, -i < s < 0.

This problem illustrates the type of discontinuities that may occur in the

kernel of the gain, and how they are affected by the approximation scheme.

Since the problem has delays of i. and .5, Kl(t,a) will have a fixed

discontinuity at a = -.5 for each t. However, because of the presence of a

final cost term, Kl(t,a) will also have discontinuities at (see Delfour

[20]) aI = -i. + tf - t and a2 = -.5 + tf - t (for values of t which

100(t,a) at t = 1.9, place a I or a2 in [-i,0]). The approximations K 1

1.4, I., and 0. are given in Figures 5.33 - 5.36. Note that Figure 5.33

should show discontinuities at a = -.4 and a = -.9, and Figure 5.34 should

show a discontinuity at a = -.4. The smoothing that occurs is due to the
-51-

averaging of values which the Runge-Kutta scheme uses to produce each new

iterate, and is not a characteristic of the Chandrasekhar algorithm. Though

not exhibited here, use of a lower order Euler scheme (which does not average

values for the next iterate) shows definite jumps at the points of

discontinuity, but sacrifices overall accuracy. The effect of the smoothing

of the kernel is minimized, however, since the action of the kernel is through

an integral term. Consequently, accurate values for the optimal trajectories

and optimal control are still obtained.
-52-

APPENDIX

PROOF OF _ 3.1

Let RN(t) s and ^ *(t) be the Nth approximation analogues of (3.5)

^N

and (3.7). Note that RN(t) is self-adjoint and satisfies fIR(t)ll > M for S S --

each N, s, and t. It follows that If[R ]-l(t)li _ uniformly in N, s,

and t. Also, it follows that IT_.)(t), (T_.)(t), and (F_.)are uniformly

bounded in N, s, and t. Combining these observations with the uniform

bounds on BN and QN, it follows that llB_(t)i; is uniformly bounded in N,

s, and t. Therefore, ;fIR -1 *(t)ll is also uniformly bounded. Recalling

that [R ]-l(t) is self-adjoint and that the norm of a bounded operator

equals the norm of its adjoint, we have established the lemma.

PROOF OF LEMMA 3.2

We first observe that strong, pointwise convergence of TN(t) and

[TN(t)] * along with an application of the dominated convergence theorem

implies that (T_.)(t), (T_*.)(t), and (F_.) converge strongly and pointwise

to (Ts.)(t) , IT_')(t), and (Fs.), respectively. In addition, (T_.)(t),

IT_*.)(t), and IF_*')(t) converge strongly and pointwise to (Ts.)(t),

Ir_.)(t), and (F$.)(t), respectively. Since BN, and QN converge strongly,

it follows that RN(t)s . _s(t) and B_*(t) . B_(t) strongly and pointwise

for 0 _ s _ t _ tf. Taking the adjoint of (3.6) and using similar reasoning,

it can be shown that B_(t) + is(t) strongly and pointwise for

^ --I

0 < s < t < tf. The identity R_ (t)- R-l(t) RN-I(R s ^ )R] = - RN l(t) now -- -- -- S S S

implies that RN-I s (t) + R-l(t)s strongly and pointwise. Therefore, it now
-53-

-i

follows that BN R_ (t) . BR-I(t) S strongly and pointwise for

0 _ s _ t _ tf. As indicated previously, the convergence of sN(t,s) to

S(t,s) was shown by Gibson [25]. Using equation (3.8) we see that

sN_(t,s) may be expressed as

(A.I) [sN(t,s)] * z = [TN(t - s)]* z - ft Bs[RsN]-I(_)[ ^BN]*[TN(t - n) zd_.

S

Combining the assumptions on TN*(t) and B N* with Lemma 3.1, the integrand

of (A.I) is seen to be uniformly bounded in N, s, and t. As a result of the

comments above, the integrand also converges pointwise to

^ R-I(n)B* T*( Bs s z - n)z, and an application of the dominated convergence

theorem yields sN*(t,s) . S*(t,s) strongly and pointwise. The uniform

boundedness of sN*(t,s) (and hence sN(t,s)) follows directly from (A.I).

PROOF OF THEOREM 3.3

An integration of (3.14) results in the equation

tf

(A.2) KN(t)z = f R-I_N]*[LN(_)]* LN(_)zdB.

t

Equation (3.15) may be rewritten as

LN(t)z = -LN(t)A N z + LN(t)_ KN(t)z

and a variation of parameters formula yields
-54-

tf

(A.3) LN(t)z = VN TN(tf - t)z - f L(_)B N KN(_)TN(_ - t)zd_.

t

Using the assumption that LN(t) = V N sN(tf,t), (A.2) and (A.3) become

tf

N* VN

(A.4) KN(t)z = f R-I BN* S (tf,n)V N* sN(tf,n)zdn,

t

and

tf

(A.5) VN sN(tf,t)z = VN TN(tf - t)z - f VN sN(tf,n)B N KN(n)TN(n - t)zdn,

t

respectively. The dominated convergence theorem, Lemma 3.2, and Gibson's

convergence results (see [25]) for KN(t) imply that the limit as N + = of

(A.4) exists and

tf

(A.6) K(t)z = f R-I B S (tf,n)V VS(tf,n)zdn.

t

Since KN(t) and TN(t) are each strongly continuous in t and converge

to K(t) and T(t) in the strong sense, they are uniformly bounded in N

and t on compact t-intervals. Therefore, another application of the

dominated convergence theorem and Lemma 3.2 imply that the limit as N . = of

(A.5) exists and

tf

(A.7) VS(tf,t)z = VT(tf - t)z - _ VS(tf,_K(_)T(_ - t)zdn.

t

Defining L(t) = VS(tf,t) we have shown that the pair K(t), L(t) satisfy

the Chandrasekhar equations (3.16) - (3.17) and KN(t)z + K(t)z for all

zEH. If U is finite dimensional, then it follows (see Gibson [25]) that

KN(t) . K(t) in the uniform operator topology.
-55-

PROOF OF THEOREM 3.4:

Let (K(t),L(t)) be another solution pair for (3.16) - (3.17). It

follows from (3.16) that

tf ^_ ^

ll(K(t) - K(t))zll < f fiR-I /3(L (n)L(_) - L (_)L(n))zlld_

t

tf

_< fi-RI B*Ill ilL^*(_)L(rl) ^ -L *(rl)L(n)ll flzlld_,

t

which implies that

tf

IlK(t)- K(t)il< nR-I 8*tif I^lL*(_)L(n) ^ -L *(r_)L(_)Ildn.

t

Adding and subtracting the appropriate term yields the inequality

tf ^ ^_ ^ _

11K(t) - K(t)ll < fI-IR /_*llf (ilL (n)L(_) - L (_)L(rl)ll+ IlL(_)_,(B)

t

- L (_)L(_)it)dR

tf

< UR-I B* llf (iIL(n)llII^*L (n) -L * (n)ll+ IlL* (n)l] ll^ L(n)

t

- L(_)ll)d_.

Since the norm of an operator equals the norm of its adjoint, it follows that

tf

liK(t) - K(t)ll _ fI-IR 13*Ill (IIL(B)II+ tIL(n)ll)llL(n) - L(n)lld_.

t
-56-

The principle of uniform boundedness implies the existence of a constant

M 1 _ 0 such that

A tf ^

(A.8) IlK(t) - K(t)ll ! M 1 f IlL(n) - L(n)lldn.

t

Similarly, equation (3.17) yields the inequality

A tf A ^

ilL(t) - L(t)ll ! f IIL(n)BK(n) - L(n)BK(q)II ;IT(q- t)lldn

t

tf

_< f llr(n- till Itl_.(nlBK(n)- L(nl_K(n) + L(nlBK(n)

t

- L(n)BK(n)lldn

tf

_< f ;IT(n- t)ll (IIL(n)B;IIlK(n)- K(n)ll

t

+ fiLCh) - L(n)II llBK(n)ll)dn.

Again, the principle of uniform boundedness can be applied to obtain the

existence of M 2 > 0 such that

tf

(A.9) ll[(t) - L(t)ll _< M 2 _ (liE(n) - L(n)ll + IlK(n) - K(n)ll)dn.

t

Substituting (A.8) into (A.9) yields

^ tf A t ^

(A.10) liE(t) - E(t)ll _< M2 f (liE(n) - E(n)ll + M I f liE(s) - E(s)llds)dn.

t n
-57-

Since t _ _ in (A.10), the lower limit may be extended to obtain

tf tf tf ^

llL(t) - L(t)ll _ M2 f IlL(n) - L(_)fldn + M I M2 f f IrE(s)- L(s)lldsd_.

t t t

Applying Fubini's theorem and interchanging the order of integration, we have

tf

ilL(t) - L(t)il d (M2 + M 1M2(t f - t)) f lIL(s) - L(s)llds,

t

and hence

tf ^

(A.II) llL(t) - L(t)ll d (M2 + M I M 2 tf) f llL(s)- L(s)11ds.

t

In order to justify the existence of each of the integrals appearing above, we

note that each operator is strongly continuous on [0,tf], and hence the norms

of these operators are lower semi-continuous on [0,tf] (see Kato [34]).

Lower semi-continuity implies measurability, and since the norms are uniformly

bounded, they are integrable.

Gronwall's inequality holds for integrable functions almost everywhere

(Reid [41]), hence (A. II) and Gronwall's inequality imply that

iiL(t) - L(t)tl = 0 almost everywhere in [0,tf]. Since llL(t)z - L(t)zil

continuous for each zEH and equal to zero almost everywhere, it follows that

L(t) = L(t) everywhere on [0,tf]. Equation (A.8) immediately yields

^

K(t) = K(t) for tC[0,tf], and the proof is complete.
-58-

REFERENCES

[I] A. V. Balakrishnan and J. W. Edwards, Calculation of the transient

motion of elastic airfoils forced by control surface motion and

gusts, NASA TM-81351, NASA Dryden Flight Research Center, Edwards,

CA, 1980.

[2] H. T. Banks and J. A. Burns, An abstract framework for approximate

solutions to optimal control probelms governed by hereditary

systems, Proc. International Conference of Differential Equations,

H. A. Antosiewicz, ed., Academic Press, New York, 1975, pp. 10-25.

[3] H. T. Banks and J. A. Burns, Hereditary control problem: numerical

methods based on averaging approximations, SIAM J. Control Optim.,

16 (1978), pp. 169-208.

[4] H. T. Banks, J. A. Burns, and E. M. Cliff, Parameter estimation and

identification for systems with delays, SIAM J. Control Optim., 19

(1981), pp. 791-828.

[5] H. T. Banks, J. A. Burns, E. M. Cliff, and P. R. Thrift, Numerical

solutions of hereditary control problems via an approximation

technique, LCDS Technical Report 75-6, Brown University,

Providence, RI, October 1975.
-59-

[6] H. T. Banks and F. Kappel, Spline approximations for functional

differentlal equations, J. Differential Equations, 34 (1979), pp.

496-522.

[7] H. T. Banks, I. G. Rosen, and K. Ito, A spline based technique for

computing Riccati operators and feedback controls in regulator

problems for delay equations, ICASE Report No. 82-31, NASA Langley

Research Center, Hampton, VA, 1982.

[8] J. S. Baras and D. Lainiotis, Chandresekhar algorithms for linear timevarying distributed systems, J. Information Sci., 17 (1979), pp.

153-167.

[9] R. L. Bisplinghoff, H. Ashley, and R. L. Halfman, Aeroelasticity,

Addison-Wesley, Reading, MA, 1955.

[I0] J. A. Burns and E. M. Cliff, Hereditary models for airfoils in unsteady

aerodynamics, numerical approximation and parameter estimation,

AFWAL-TR-81-3173, Flight Dynamics Laboratory, Wright-Patterson Air

Force Base, OH, 1982.

[ii] J. A. Burns, E. M. Cliff, and T. L. Herdman, A state-space model for an

aeroelastic system, Proc. 22nd IEEE Conference on Decision and

Control, Vol. 3, San Antonio, TX, 1983, pp. 1074-1077.
-60-

[12] J. A. Burns, K. Ito, and R. K. Powers, Chandresekhar equations and

computational algorithms for distributed parameter systems, Proc.

23rd IEEE Conference on Decision and Control, Vol. I, Las Vegas,

NV, 1984, pp. 262-267.

[13] J. L. Casti, Dynamical Systems and Their Applications: Linear Theory,

Academic Press, New York, 1977.

[14] J. L. Casti, Matric Riccait equations, dimensionality reduction, and

generalized X and Y funtions, Utilitas Math., 6 (1974), pp. 95-110.

[15] J. L. Casti and L. Ljung, Some new analytic and computational results

for operator Riccati equations, SIAM J. Control Optim., 13 (1975),

pp. 817-826.

[16] J. L. Casti and L. Ljung, Reduciton of the operator Riccati equation,

Proc. International Symposium on Control Theory, Numerical Methods,

and Computer Systems Modelling, Paris, A. Bensoussan, J. L. Lions,

eds., 1974, pp. 384-396.

[17] E. M. Cliff and J. A Burns, Reduced approximations in parameter

identification of hereditary systems, Lecture Notes in Control and

Information Sciences, A. V. Balakrishnan and M. Thomm, eds.,

Springer-Verlag, New York, 1982.
-61-

[18] R. Curtain and A. J. Pritchard, The infinite dimensional Riccati

equations for systems defined by evolution operators, SlAM J.

Control Optim., 14 (1976), pp. 951-983.

[19] R. Datko, A linear control problem in abstract Hilbert space, J.

Differential Equations, 9 (1971), pp. 346-359.

[20] M. C. Delfour, The linear quadratic control problem for hereditary

differential systems: theory and numerical solution, Appl. Math.

Optim., 3 (1977), pp. 101-162.

[21] M. C. Delfour, E. B. Lee, and A. Manitius, F-reductlon of the operator

Riccati equation for hereditary differential systems, Automatica -

J. IFAC., 14 (1978), pp. 385-395.

[22] M. C. Delfour and A. Manltius, The structural operator F and its role in

the theory of retarded systems I, Centre de Recherches

Mathematiques, CRM-765, University de Montreal, March 1978.

[23] M. C. Delfour and S. K. Mitter, Controllability, observability and

optimal feedback control of afflne hereditary differential systems,

SlAM J. Control Optlm., i0 (1972), pp. 298-328.

[24] J. W. Edwards, Unsteady aerodynamic modeling and active aeroelastic

control, SUDAAR 504, Stanford University Report, February 1977.
-62-

[25] J. S. Gibson, Linear-quadratic optimal control of hereditary

differential systems: infinite dimensional Riccati equations and

numerical approximations, SlAM J. Control Optim°, 21 (1983), pp.

95-139.

[26] J. S. Gibson, The Riccati integral equations for optimal control

problems on Hilbert spaces, SlAM J. Control Optim., 17 (1979), pp.

537-565.

[27] K. Ito and R. K. Powers, Chandrasekhar equations for infinite

dimensional systems, NASA CR-172525, ICASE Report No. 84-67, NASA

Langley Research Center, Hampton, VA, June 1984.

[28] K. Ito and R. Teglas, Legendre-tau approximation for functional

differential equations, NASA CR-172155, ICASE Report No. 83-17,

NASA Langley Research Center, Hampton, VA, June 1983.

[29] K. Ito and R. Teglas, Legendre-tau approximation for functional

differential equations, part II: the linear quadratic optimal

control problem, NASA CR-172397, ICASE Report No. 84-31, NASA

Langley Research Center, Hampton, VA, July 1984.

[30] R. T. Jones, Operational treatment of the nonuniform llft theory to

airplane dynamics, NACA TN-667, 1938.
-63-

[31] T. Kailath, Some Chandrasekhar-type algorithms for quadratic regulators,

Proc. IEEE Conference on Decision and Control, New Orleans, LA,

1972, pp. 219-223.

[32] T. Kailath, Some new algorithms for recursive linear estimation in

constant linear systems, IEEE Trans. Inform. Theory, IT-19 (1973),

pp. 750-760.

[33] F. Kappel and D. Salamon, Spline approximations for retarded systems and

the Riccati equation, MRC Technical Report No. 2680, University of

Wisconsin, Madison, WI, Aprll 1984.

[34] T. Kato, Perturbation Thoery for Linear Operators, Springer-Verlag, New

York, 1966.

[35] J. L. Lions, Optimal Control of Systems Governed by Partial Differential

Equations, Springer-Verlag, New York, 1971.

[36] D. L. Lukes and D. L. Russell, The quadratic criterion for distributed

systems, SlAM J. Control Optim., 7 (1969), pp. 101-121.

[37] A. Manitius, Feedback controllers for a wind tunnel model involving a

delay: analytical design and numerical solution, IEEE trans.

Automat. Control, 29 (1984), pp. 1058-1068.
-64-

[38] A. Manitius, Optimal control hereditary systems, Control Theory and

Topics in Functional Analysis Vol. III, International Atomic Energy

Agency, Vienna, 1976.

[39] A. Manitius and H. Tran, Computation of closed loop eigenvalues

associated with the optimal regulator problem for functional

differential equations, Proc. Conference on Decision and Control,

San Antonio, TX, 1983.

[40] R. Powers, An application of Chandrasekhar-type algorithms for optimal

control of delay-differential systems, Proc. 15th Southeastern

Symposium on System Theory, Huntsville, AL, 1983.

[41] W. T. Reid, Properties of solutions of an infinite system of ordinary

linear differential equations of the first-order with auxiliary

boundary conditions, Trans. Amer. Math. So¢., 32 (1930), pp. 284-

318.

[42] M. Sorine, Schema d'approximation pour des _quations du type de

Chandrasekhar intervenant en contr$1e optimal, C. R. Acad. Sci.

Paris S_r. I Math., t. 284 (1977), pp. 61-64

[43] M. Sorine, Sur les _quatlons de Chandrasekhar associ_es au probl_me

de controle d'un syst_me parabolique: cas du controle et de

l'observation distribu_s, C. R. Acad. Sci. Paris S_r. I Math., t.

285 (1977), pp. 863-865.
-65-

[44] M. Sorlne, Sur le semi-groupe nonlin_aire associ_ _ l'equation de

Riccati, Rapport INRIA, No. 167, 1982. ,_

[45] F. Tung, Identification of Aircraft Parameters in Turbulence with NonRational Spectral Density, Ph.D. Thesis, University of California,

LA, 1979.

[46] D. L. York, Analysis of Flutter and Flutter Suppression via an Energy

Method, M.S. Thesis, Aerospace and Ocean Engineering Department,

Virginia Polytechnic Institute and State University, Blacksburg,

VA, June 1980.
-66-

2.50

1.6 2.0

6.

o

.8 1.2

u

o N =8

6. N =16

.50

o

2.00

1.00

1.50

U(T)

T

FIGURE 5.1
-67-

2.50

.8 1.2 1.6 2.0

-u

C:. N =48

.4

.50

o

LOO

1.50

2.00

U(T)

T

FIGURE 5.2
-68-

2.50

2.00

1.50

U (T)

--U

1.00 z_N = I00

.50

l I I L____L I

0 .4 .8 1.2 1.6 2.0

T

FIGURE 5.3
-69-

-.20 -

-.25

__,,

-o30

<> N=8

"_ ,_ N=32

K (O,A) - o _ -- N = 100

-.35 - "'_'_

,,\

-\

-.40- o o __//

_

-.45 i I I I I I I I I I

-I. -.8 -.6 -.4 -.2 0

A

FIGURE 5.4
-70-

-.35 -

K (.25, A) - <> N= 8 o

z_ N=32 ,,"

-'_T -. : loo o ,_//

o

°\ o _;/

-.65 - _N,_..o o ,, ,,_

A Z_

_

-°75 . I I_ I I I .I I 1 I I

-1. -.8 -.6 -.4 -.2 0

A

FIGURE5.5
-71-

-.9

K L5, A)

-1.1

o

o N = 8

Ll N =32

- N = 100

-1.3

-.8 -.6 -.4 -.2 a

-1.5 "-_...Io.-_-'--_~_--,-_----,-_---a __.a...-_...I..-_-'-_~

-1.

A

FIGURE 5.6
-72-

m

-.30-

- o#y z_

-.60 - <>N = 8 o ,,/_/

A

/

nN=32 ,,/

K L75,A) - --N = I00 o ,,/

& /

-.90 - o A_/i /

/

-- 0 _/

oY

oY

-1.20 - _> ",-5"/"

A _

-1.50 I I I I I I I I I !

-I. - 8 -.6 -.4 -.2 0

A

FIGURE 5.7
-73-

-.4 - o

z_/

oY

-.8 - o _/

/

K (I.A) - ,,Y

' o ,,/

-1.2- "Y

o "/

__ ,,>/

-1.6 - "/ z_N = 32

,,,_"/ --N = I00

-2.0 I I I I I I I I I I

-I. -.8 -.6 -.4 -.2 0

A

FIGURE5.8
-74-

-.50 -

o

-1.00 - ,,,_

K (1.25,A) - ,,/

z_/ /

-1.50- o _/

_- _'/ ON=8

-2.00 / --AN=32 N = 100

-2.50 l I i I I I I i I /

-1. -.8 -.6 -.4 -.2 0

A

FIGURE 5.9
-75-

0 - _,, .,,-. _ <> A 0

A

z,.<)

-.4 A o

A A

-.8

K (1.5,A)

-1.2

¢N=8

z_N 32

--N I00

-1.6-

-2.0 I I I I I I I I. I I

-i. -.8 -.6 -.4 -.2 0

A

FIGURE 5.I0
-76-

0

O

K (1.75A),

-.6 ON=8

_N=32

--N= 100

-1.0 1 I I I I I I I I 1

-1. -.8 -.6 -.4 -.2 0

A

FIGURE 5.II
-77-

10.0

7.6

m

5.2 -

X1 (T) -

2.8 -

- z_N=8

--N

.4 --

-2.0 I I I I I I I I ] I

0 .4 .8 1.2 1.6 2.0

T

FIGDq_5.1 E 2
-78-

o

-2

-4

X2 (T) /

-6

A

fJ.N=8

-8 - N = 100

.4 .8 1.2 1.6 2.0

-10 L-_...J.-_----L-_--J._---I__"'"'--_...J.-_......l-_---L-_----I._-----l

o

T

FIGURE 5.13
-79-

.23 tf = 10

=3

.21

K (0,A1

.19

N=48 .17 -

.15 I I I I I I I I I I

-.33 -.264 -.198 -.132 -.066 0

A

FIGURB 5.14
-80-

__ f

-.04

-.08

XI (T)

-.12 -

N=48

-.16 -

-.20 - I I I I I 1 I i I 1

0 2 4 6 8 10

T

FIGURE 5.15
-81-

6

-2

-10

X2 (T)

-18

N=48

-26

-34 I I I I I I I I I I

0 2 4 6 8 10

T

FIGURE 5.16
-82-

10.0

7.0

4.0

X3(T)

1.0

-2.0

N=48

-5.0 I I I I I I I I 1

0 2 4 6 8 10

T

FIGURE 5.17
-83-

0

-1.0

-2.0

u (T)

-3.0

-4.0

N= 48

-5.0 I I I I . ! I I I I I

0 2 4 6 8 10

T

FIGURE 5.18
-84-

8,-------------------------,

6

4

2

Xl (\) a

-2

-4

N= 16

0.05 O. 1 O. lS 0.2 0.25

-f!>L..-. --"- ...L- ---I --.J... ---J.

a

T

FIGURE 5.19
-85-

N = 16

2-

\

-2 /

/

-3 --./

!

-4 / I I

0 0.05 O.I O.15 0.2 0.25

T

FIGURE5.20
-86-

O.12

o. 1 N 16

0.136

0.06

0.04

\

0.02

X3 (T) ',,.

o ".,.

-0.02

-.0,04

-0.06

-0.08

-o.1 I I

o o,o5 o.! o.is o.2 o.2s

T

FIGURE5.21
-87-

N= 16

..........'..'......... o '., _ _........•........_ - _......... _ .

0.03

0.01

0.02

-0.02

-0.01

X4 (T)

0.05 O. 1 O. 15 0.2 0.25

-0.03 1..- --J, --L. -L.. -1- ---I

o

T

FIGURE 5.22·
-88-

6O

50--

_ N = 16 •

40--

X5 (T)

10

0

/

t

-I0 /

-m I I I

0 0,05 o,_ o._5 0.2 0.25

T

FIGURE5.23
-89-

2

0 _ .....

-2 - N = 16

-4 -

-8 -

U (T)

-10 -

-12 -

-14. --

-16 --

__ I I I I

o 0.05 o.1 o.,5 0.2 o._3

T

FIGURE5.24
•

5

3

•..

..

..

-90-

..

1

Xl <T) 0

-1

-3 ..

...

• ..

•

• ..

•

...

t •

..

• • ...

•

.. ...

•

.. .. • • ....

+ Unforced

Closed loop

N =48

...

• •

.. .....

..... .. .

-5

o .10 .20

T

FIGURE 5.25

.30 .40 .50
-91-

5

3

• _, 4,

X2 (T) 0

tl- _ t11?,

-1 .

+Unforced

-.3 - Closedloop

N =48

-.5 I I ! I I I I I

0 .10 .20 .30 .40 .50

T

FIGURE 5.26
-92-

.06 -

.04 .

01 , " "*"

X3 (T) 0 I

4

m • • '

• !

+ Unforced

_.04 - • .. - Closedloop

" N=48

-.05 "'1 ! I I [ I I I I I

0 .10 .20 .30 .40 .50

T

FIGURE5.27
-93-

.04

.02

.01

.. .. ..

. .

+........ ." .. +. .+ ......

....... .

•• ••

.. .. .. " .. +.++...t

+ Unforced

- Closed loop

N= 48

•

.... .. . • to

..

•

..

. .. ••t

.. . •

..

.... ... ..

•

•

to ..

••

•

"

..

+•

•

..

..

..

•

"

..

•

o

-.01

-.02

X4 (T)

-.04

o .10 .20 .30 .40 .50

T

FIGURE 5.28
-94-

55

39

23

X5 (T)

4_ o

t

-9 + Unforced

- Closedloop

t

" N=48

-25 ! I 1 ! I I I I I

0 .10 .20 .30 .40 .50

T

FIGURE 5.29
-95-

-3 10x 10

m

X1 -

m

o _-_ I 1 ! I_ ! I_

O. .05 .10 .15 .20 ,,25

T

FIGLq_E5.30
-96-

-5 50x I0

40-

3O

X3 _

20 -

10 -

o_ [ l J l 1 l l

O. .05 .lO .15 .20 .25

T

FIGURE 5,31
-97-

2.0

_0 0°°170

1.0 o o

o0000000

0

0

C]

0 o ooo_oooooooooooooogg

O lO000000000r1000000

K IT,A) o

-1.0

0

-T=O

z_ T = .15

n T=.2

0 T = .225

-2.0 N = 48

-3.0 I I ! I I ! t I

-.05 -.04 -.03 -.02 -.01 0

A

FIGURE5,32
-98-

5°0

4.0

3.0

K (1.9, A) -

2.0 --

N= I00

Io0 --

0 I L I ! l I I I J

-1. -.8 -.6 -.4 -.2 0

A

FIGURE 5.33
-99-

3.0 --

K (1.4,A) -

2.0-

N- 100

1.0-

O, I i I I I I I I i 1

. -.8 -.6 -.4 -.2 0

A

FIGURE 5.34
-I00-

3,0

K (I.0A), -

2.0-

1.0 - N = I00

0 1 ! 1 I I I I I 1 I

-1. -.8 -.6 -.4 -.2 0

A

Fibre 5.35
-101-

.5.0 --

4.0

3.0

K (0,A)

2.0

1.0 - N = I00

0 I ! I I I I I I I I

-1. -.8 -.6 -.4 -.2 0

A

Figure 5.36
1 ReportNo. NASA CR-178033 2. Government Acc_4ion No. 3. R_ipient's CmUlog No.

ICASE Report No. 85-58

4. T,tleand Subtitle 5. Report Date

November 1985

FACTORIZATION AND REDUCTION METHODS 6, PerformingOrganizatioCnode

FOR OPTIMAL CONTROL OF OISTRIBUTED PARAMETER SYSTEMS

7. Author(s) 8. Performing Organ;zation Repo4"tNo.

John A. Burns and Robert K. Powers 85-58

10. Work Unit No.

9. PerformingOrganization Name and Address

Institute for Computer Applications in Science

and Engineering 11. Contract or Grant No.

Mail Stop 132C, NASA Langley Research Center NASI-17070; NASI-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

Washington, D.C. 20546 14. SponsoringAgencyCode

505-31-83-01

15. Supplementary Notes

Langley Technical Monitor: Submitted to Mathematica Applicada E

J. C. South Jr. Computacional

Final Report

16. Abstract

A Chandrasekhar-type factorization method is applied to the linearquadratic optimal control problem for distributed parameter systems. An

aeroelastic control problem is used as a model example to demonstrate that if

computationally efficient algorithms, such as those of Chandrasekhar-type, are

combined with the special structure often available to a particular problem,

then an abstract approximation theory developed for distributed parameter

control theory becomes a viable method of solution. A numerical scheme based

on averaging approximations is applied to hereditary control problems.

Numerical examples are given.

"17. Key Words ISuggested by Author(s)) 18. DistributiOn Statement

hereditary control 59 - Mathematics and Computer Science

LQR (General)

reduced computations

Unclassified - Unlimited

19. Security Classif.(of this report) 20. Security Classif. (of this page) 21. No. of Pages 22. Price

Unclassified Unclassified 103 A06

Forsalebythe NationalTechnicalInformationService,Springfield,Virginia 22161