Distributed-parameter optimal control via mathematical programming

Distributed-Parameter Optimal Control Via Mathematical ProgrammingT by

S.G.TZAFESTAS

Control Systems Laboratory,

University of Patras, Patras, Greece

ABSTRACT: In this paper a variety of optimal control (OC) problems for distributedparameter (DP) systems are approached using mathematical programming (MP). First, the principal DP models in current use are given, a variety of DP objective functions is provided, and the OC problems based on them are formulated. Second, these models and objective functions are converted in algebraic form, as required by MP, and the solution procedure of the OC problems via MP is outlined. Third, a representative set of nonlinear programming results applied to DP systems is presented, and finally, a number of application examples is given.

I. Introduction Mathematical programming (MP) techniques have found increasing use in many areas of modern life. One such area is optimal control (OC) which embraces technological, economic and operational research/management applications. Following the initial work dealing with MP and OC (l), a large amount of research has been perfomed in this area; it is reviewed in (2-4). An OC problem of a finite-dimensional dynamic system is, in fact, equivalent to a MP problem in infinite-dimensional space. Therefore, the main research effort has concentrated on the reduction of this problem to an approximate finitedimensional problem amenable to conventional MP theory. A dynamic system under control may be given in continuous-time form or in discrete-time form. The latter is treated directly, while the former is transformed into the discrete-time form with specified or unspecified sampling intervals. Work on stochastic systems using MP techniques includes stochastic compensator/controller design, identification, and filtering (2-3). Some salient features of the treatment of system design and control problems by MP [especially by linear programming (LP), and quadratic programming (QP)] are: (a) The solution can be determined easily for both fixed and free endpoints. (b) Equality and/or inequality constraints on state and/or control variables can be easily accommodated and treated. (c) The states of the system can be specified, not only at the terminal points, but also at other intermediate points. (d) The computational effort is generally reduced. ? This is a revised and amplified version of a paper presented Simulation, Sorrento, September 1979.

at IMACS Congress on

S. G. Tzafestas Optimal control and system design problems for distributed-parameter systems (DPS) are difficult to treat, since a DPS has an infinite-dimensional state-space. One has therefore to approximate the system by one in finitedimensional space, and also the objective function by a finite-dimensional (algebraic) function. This calls for time and space discretization, which must be accurate while not leading to excessive numbers of variables. Clearly, the advantages of applying MP techniques carefully to OC problems are expected to be greater in the case of DPS. Our purpose in this paper is to present a comprehensive treatment of OC problems for DPS via MP, including existing results (5-12). The dynamic programming approach to DP control, which is included by some authors in MP, will not be considered here [see (13-14)]. Other techniques which are used in DP control are: the maximum principle (l&16), variational calculus (17), and functional analysis (18). Two survey papers encompassing most of the work on DPS are (19) and (20).

II. Distributed-Parameter

Models

Models for both linear models are the differential 11.1. Linear

differential

and nonlinear DPS will be described. The basic equation model and the integral equation model.

model

A unique differential equation model, covering all possible DPS encountered in practice, does not exist. This is particularly true in the nonlinear case. The most popular linear model is the state-space one, i.e. Xt(x, t) =%X(x, PX(x, t)=O

t) + U(x, t) XEaD;

x E D;

X(x, 0) = X,(x),

tzo,

x E D;

I

(1)

where Xt(x, t) =8X(x, t)/dt; D is a spatial domain with boundary *D; X(x, t), x E 0, is the state vector function; X,(x), x E D, is the initial condition; (3, /3} is a well-posed pair of spatial operators; and U(x, t), x E 0, is the vector control function over D. Sometimes the system also involves boundary control Ub(x, t), x E I’D. The boundary condition can be written as t) = Ub(x, t),

pX(x,

x E 8D.

W-4

Clearly, it is possible to have just U(x, t), or just Ub(x, t), or both. This model can, with care, be made to accommodate heating-diffusion-like systems, wavelike systems, heat exchangers, ‘etc. (5-20). 11.2. Linear When both form (18)

integral equation U(x, t) and

model

L&(x, t) exist, the integral

equation

model

takes the

t X(x,

t) =

X0(x,

t) +

K(x,

H0

D

is0

6D

f

+

400

K,(x,

r;

5, ~)U(&T)

t;

&,,

7) U,(S,,,

d5

d7

7)

d&

dT

(2) Journal of The Frankfin Institute Pergamon Press Ltd.

Distributed-Parameter

Optimal Control Via Mathematical

with X%,

J

t) =

K(x, t; 5, 4,)X&)

Programming

d5,

D

where the matrix kernels K(x, r; 5, 7) and &(x, either analytically, or experimentally by measuring testing inputs at different spatial points (21,22). 11.3. Lumped-parameter

t; 4,~) can be determined the system response due to

and pointwise control models

When the system is subject to lumped-parameter L&(x, t) in (1) and (2) take the forms U(x, t) = U”(x) U(t),

XED;

U(x, t) and

control,

Q(x, t)= U;(x)U,(t),

XEIYD,

(3)

where U*(x) and U:(x) are known spatial functions, and U(t) and I.&(t) are the lumped-parameter controls. An important example, treated in (5,8, lo), is a special case of the composite distributedand lumped-parameter model studied in (16). The model in this case is given by Xt(x, t)=X,,(x, X,(0,

t),O
t)=c-u[X(O,

&t)+

x(x,0)=x,(x),

t)-V(t)];

OCxSl;

t)=O;

X,(1,

(4

v(t) = Q(f).

One finds the integral

model

X(x, t)=X’(x,

t)+

J

‘g(x, t-T)U,(T)dT,

(5)

0

where

x0(x,

t) and g(x, t) are determined

x”(x, s) = X,(x, s) + a

by Laplace

[(sf sinh si)A(O,

s) + A(l,

inversion

of

s)] cash (1- x)sf

si sinh &[s; sinh si+ cy cash s”]

cash (l-x)&

G(x’ ‘) = (1 + ys)[s; sinh s! + (Ycash sf] ’

A(0, s) = &(O, s)- ajz;(O, s), A(l,

Xh(x, s) = aZo(x, s)/dx,

s) = c&,(1,

’ (6)

s),

and &(x, s) is determined from X,,(x). For example, if X,(x) = Z,, (a constant), then X0(x, s) = Eo/s, xh(x, s) = 0, x0(0, s) = Eo/s, and a(8,/s)

x0(x,

t) = 28, f

cash (1 - X)S~

0, (x)e-“‘,

(7)

7ri = p:,

,=l

@(x) = where

pi are the roots

Vol. 309, No. 6, June 1980 Printed in Northern Ireland

ices (lvx)Pil p:{( 1+ a)/Pf + l/(11}cos pi ’

J

of p tan p = (Y.

401

S. G. Tzafestas When the system is subject to pointwise control (e.g. a nuclear control rods located at the points xk E D), we have U(x, t) = z

U,(t)s(x

reactor

with

-x!J,

k=l

ub(x,

where

t)

=

6(a) is the unit impulse x(x,

f k=md+l

Uk(t)a(x

-

(8)

xk), 1

function.

The model

(2) reduces

to

x0(x, t) + 'K(x, t; T)U(T) d7,

t) =

(9)

where K(x, t; T)=[K(x,

The model by inserting

t;

x1,

T), . . . , K(x,

t; x,, T)];

(9) also covers the case of lumped-parameter (3) into (2), or by (5).

11.4. Nonlinear

control,

as is shown

models

For the nonlinear case, there models: the model (14,17) Xt(x, t) = N,(X, Nb(X,x,

appear

x, t) + U(x, t),

t)=O,xeD;

and the hyperbolic X(x,

model

to be two main

differential

x ED;

X(x,0)=X,(x);

t20,

XED;

equation

>

(10)

(10,23)

0 = FAX, 6 X, Xx, Xx,, U(x,

01, x E(0, 1);

X(x, 0) = X,(x), x E (O,l); Fbo[x, t, Xx, UbO(f)] = 0,

x = 0;

&J-x,

x = I;

6 Xx, &l(f)1

= 0,

(11)

t 20, I

where Nd and Nb are well-posed nonlinear spatial operators, Fd, Fbo, and Fbl are sufficiently smooth conventional functions, and U(x, t), Ubo(f), and Ubl(f) are the control inputs. The nonlinear integral DP model is (15) X(x, t) = X%,

t) + lot ID GCx, t; 5, 7, W5, ~)l d5 dr f +

II0

G,Cx,

t;

Q,

7,

&(SI,,

T)l

d$

(114

dT,

lYD

where x0(x, t) is the response component due to the initial condition, and Gb are conventional nonlinear functions of their arguments.

402

and

G

Journal of The Franklin Institute Pergamon Press ud.

Distributed-Parameter 11.5. Stochastic

Optimal

Control

Via Mathematical

Programming

models?

The principal stochastic models are obtained by inserting terms in the models described by (l), (2), (5), (9), (lo), (ll), for example, the stochastic counterpart of model (1) is (24) X,(x, t) = YX(x, /3X(x, t) = 0, Here W(x, t) is a zero-mean covariance matrix

stochastic and (lla).

t) + U(x, t) + W(x, t); (12)

x E IYD.

1

DP Gaussian

and white

(in time)

process

ELW(x,t) W’(5, ~)l= Qb,5,tP%r - 71, where E[.] Similarly,

is the expectation operator, the stochastic counterparts X,(x, t) = N(X,

with (13)

and I?(.) is the unit impulse function. of models (10) and (11) are

x, t) + U(x, t) + W(x, t),

t)=0

&(X,x,

input Thus,

(14) 1

and Xt(x, t) = Fd(. . .)+ W(x, t), I&(.

I&(.

. .)+ W,,(t)=O,

(1.3

. .)+ W,,(t)=%

I

where W(x, t) is as above, and WbO(t) and W, 1(1) are boundary zero-mean white stochastic inputs (23). Finally, the stochas.tic counterpart of model (2), when U,(&, 7) = 0 for all & E IYD, is (21)

+JJf K(x, t;6,T)w(t, 7)d5 dT> (16) K(x,

0

t;

5,

7) u(5,

7)

d7

D

where it is assumed that there is no boundary noise. All DP models presented above are continous-time models. DP models can be derived by starting from the continuous-time shown in Section IV.l. III. Distributed-Parameter The usual

objective

111.1. Energy-like

Objective Functions and Control Problems

(cost) functions

objective

This has the quadratic J= CT{ j 0

Discrete-time ones, as is

types:

function

form

x=(x,

are of the following

(T is the final time)

t)Q(x,

t)T?(x, t) dx+

U.‘.(t)R(r)U(t)}

dt

(ITa)

D

t The stochastic models outlined here are formal, since direct use is made of the white noise process (21,24). However, there is no difficulty in treating by MP, rigorous stochastic models which involve Ito integration and Hilbert-space-valued Wiener processes (67-69).

Vol. 309, No. 6, June 1980 Printed in Northern Ireland

403

S. G. Tzafestas for the lumped-parameter and point-wise control cases, where 2(x, t) = X(x, t) Xd(x, t) is the difference between the actual state trajectory X(x, t) and the nominal or desired one Xd(x, t), x E D, t E [0, T]; Q(x, t) is a semipositive, and R(t) is a positive definite matrix function. Sometimes one has an objective function involving only the state penalty term, or only the control penalty term which is known as the control effort term. In more general situations, the objective function involves a control-rate penalty term JOToT(t)R,(t)ti(f) dt and/or a penalty term of the spatial state derivative, i.e. a term Jo’JD zx(x, t)Q,(x, t)kx(x, t) dx dt, where zx(x, t) = X,(x, t) - X,,(x). These constrain the variations of o(t) and/or X,(x, t) (22,28) and so minimize the time and/or the spatial ripple of the response. terminal objective function

111.2. Quadratic J=

JD

[X(x, T) - X,(x)ITQ~(x)[X(x,

T> - X,(x)]

dx,

(17b)

where Q, is a given semipositive matrix. This implies that the final state X(x, T) must be as near as possible to the desired one Xd(x), x E D. 111.3. Terminal state distance objective function J=

2 IXi(x, T)-X&(x)(

5D i where i extends over all components

dx,

(17c)

of X(x, t), or (17c’)

J = max max [X,(x, T) - X,(x)\. I

111.4. Quadratic

xsD

terminal state and control energy objective function T

J=

[X(x, T)-X,(x)ITQ,[X(x,

T)-X,(x)]dx+

UT(t)R(t)U(t)dt. (174

111.5. Fuel-like objective function J=

i0

where i extends over all components

TV IUi(t)l dt,

(17e)

i

of U(t).

111.6. Control amplitude objective function

(17f) 111.7. Terminal state distance and fuel-like objective function FIX,(x,

404

TJ-Xh(x)/dx+/oTF

Iq(t)IdC

Distributed-Parameter where i extends over components of U(r). 111.8. Minimum-time

the

Optimal Control

of X(x, t), and

components

objective

Via Mathematical Programming i extends

over

the

function .l=

J0

Tdr=T.

O’W

In some cases one has simultaneously to minimize a group of objective functions, e.g. fuel-like and terminal distance. This is known as mulri-objective (multi-criteria) control problem. In the stochastic case, one can use the expected values of the above objective functions, for example .T=E [ lT { 1 kT(x, r)CXx, r)%x,

t) dx + UT(t)R(t)

U(t)] dt],

D

or J=E[s,

(x)3'Qf[Xb, ~Mk4ld*], -

[X(x, T)-X,

or

1

T

J=E

[I 0

dr = E[T], etc.

Now using the models of Section II and the above objective functions, we can formulate a repertory of optimal control problems for DPS of the type: “Given a DP model and an objective function, choose the control functions so as to minimize this objective function, while satisfying certain additional control and/or state constraints of the amplitude or energy type, etc.” For example, one of the usual control constraints is I W)l

s 1

or

06

U(r)sl.

(18)

IV. Solution of DP Control Problems by Linear and Quadratic Programming To formulate a DP control problem as a MP problem, one has to reduce the system model, the objective function, and all constraints to the form required by MP theory. If the model, objective function, and constraints are linear, one must employ LP theory, i.e. the simplex algorithm. If the model and the constraints are linear, and the objective function is of the energy type (i.e. quadratic), QP theory should be used. Finally, if at least one of the constraints, or the system, or the objective function is nonlinear, NLP has to be employed. IV. 1. Reduction of DP models to mathematical programming form We shall now outline reduced to the algebraic

the methods form Av=b

required

by which

linear

DP

models

can

be (19)

by LP or QP

Vol. 309, No. 6, June 1980 Printed in Northern Ireland

405

S. G. Tzafestas The main methods for a differential equation model are: (i) Finite difference method. (ii) Assumed-mode method (Galerkin, etc.). (iii) Method of lines. (iv) Laplace and Z transform or digital integrator operators. If the integral equation model is available, one can use: (v) Numerical integration formulae. (vi) Assumed-mode expansion (Walsh-Galerkin). (vii) Time discretization. IV. l.(i) Finite difference method. Using finite difference approximations for the spatial derivatives, we obtain a conventional finite-dimensional state-space model of the type Y(t)=

X(t)=AX(t)+BU(t), where A, B, and C are block circulant matrix of index p has the form

&,f=

matrices

m.

m,...m,_,

mp-l .

mp..+-2 .

where

A is said to be block

A,

are circulant

(20)

of the same index.

A circulant

I * .I

matrices

.

.

m2 . . . m,

ml

A matrix

CX(t),

circulant

of index

if it has the form

p. For example,

the DP model

X1 u1 1 [*I .

X,(x, t) = Xx,(x, t) + U(x, t), 0
x is an angular

coordinate

-2 [

/+L

a2

1 1 . 10

-2 .

such that x and x + 27~ are the same, leads to

0

...

0

1

1

...

0

0

.

. ..l-2

)

x=

x2.

%

)

u=

u2.

u,

with X,(t) = X(ai, t) and a = 25~/p. The model (20) can be converted in algebraic form by time discretization (4), or by using Legendre polynomials (3), or by Walsh function expansion (22), or through its transition function @(t) = exp (At). IV.l.(ii) Assumed-mode method. A general presentation of this method may be found in (2526). We will outline the Walsh-Galerkin method (27) here. Consider the linear DP model (l), and the family of assumed modes {qi(x)} which satisfy the homogeneous boundary conditions. Then, consider the Walsh

406

Journal of The Franklin Institute Pergam0n Press Ltd.

Distributed-Parameter

Optimal Control Via Mathematical

Programming

which constitute a complete orfunctions a(t) = [Q(t), al(t), . . . , @m-l(r)]’ thonormal family of binary functions (22, 27, 29, 30). We start by expanding Xt(x, t) = 8X(x, t)/& and U(x, t), thus:

x*(x, t, =

Tii(r)Ti(x),

i=l

U(x, t) =

2 i=l

(214

1

J

mi(r)qi(x),

with

c:@(f), m,(t) = VT@(t),

i,(t) = where obtain

c, and vi are m-dimensional

Walsh-series

vectorial

@lb) coefficients.

Then we

T

X,(x, t) df + X,(x)

X(x, t) =

=

T (c~P+c~=)@(r)lIqx),

(21c)

i=l

where c,(O) = cpT@(t) are the coefficients of the expansion X,(x) = Ci &(O)?,(x), and P is the Walsh-matrix integration operator of order m given by

Li 1 I --

1

p, = __i!?~i___~___??_!? zy& L,, :

Introducing residual

the approximate

PI=;.

(22)

%2

expansions

(21a-c)

(1) yields

the

c,(t) are chosen so that jD ET(x, t)E(x, t) dx is a minimum, = 0, i = 1,2, . . . , IV,,. This gives

or

E(x, t) = X<(x, t) - 2X(x, The coefficients ID E’(x, t)qi(x)

)

if? fiji I=1 from which one obtains

Vol. 309. No. 6, June 1980 Printed in Northern Ireland

-,El

Sij$

-

into

model

t) - U(x, t), x ED.

jg

fijU,(t) =

0,

(23)

S. G. Tzafestas where A = F-‘G,

fij

= [

F= [&I,

‘I’T(X)Wi(X)

d-X,

G = [gii],

gii =

u = I’m,

dx.

{cFP~((X)}~*~(X)

Obviously, if 9, (x) are orthonormal, then F = I. If qj (x) are orthonormalized eigenfunctions of {Y, p}, then G = diag [h,, . . . , AND], where Ai are the associated eigenvalues. Since the qi(x) terms satisfy the homogeneous boundary condition there is no residual over IYD. Using the Walsh expansions in (23) gives (p-- G‘+)c = Gc”+ Fy, from which we obtain (25a)

r = (E- G*)-9, F=I@F, Here, the symbol A@B B; it is defined by

denotes

y = (F- G*)-l&O,

C?=I@G,

&=P=@G.

the Kronecker

product

of the matrices

A and

(25c)

Clearly, an algebraic equation of the type (25a) is obtained if we apply the Walsh expansion to any state-space model such as (20) etc. IV.l.(iii) Method of lines. In this method, the DPS is converted to an equivalent lumped-parameter system, described by ordinary differential equation along one of its dimensional paths of integration or lines, in the space-time domain. These equations can be readily written in the state-space form such as (20) or (24a). Among the work dealing with the method of lines are (31, 32), while (33) describes the application of this method to DPS identification. IV.l.(iv) Laplace and z transformation method. Using the Laplace transformation method, we obtain the tranfer matrix of a multidimensional DP system. Either the differential equation or the integral equation model (5, 8) can be used as the starting point for this method. Applying the 2 transform, or using digital integrator operators (34, 35), yields the pulse transfer matrix of the system, which then provides the discrete-time model (difference equations) of

Distributed-Parameter

Optimal

Control

Via Mathematical

Programming

the system. Let G(x, z) be the transfer matrix of the DPS from the lumpedparameter control input to the DP response X(x, z), where z is the variable of the 2 transform. Then X(x, z) = G(x, z) U(z). In general,

G(x, z) is a matrix

with elements

Gij(x,2) = F

(26)

which can be put in the form (8)

Kij,l(x)z-l 1

_

a,,

*-I

(27)

7

Y.1

where

1 extends

over a non-negative Xij,Jx, kAt) =

where

At is the sampling

period.

integer

set. Now define

Kij,l(x)z-'

(28)

uj(kAt),

1- aij,l(x)z-’

This gives the discrete-time

Xii,r(~, kAt) = aij,,(x)Xii,,(x, k - lAt)+Ki,,l(x)uj(k which,

if applied

at a set of discrete

spatial

Xj,,(xq, kAt)=aij,,(x,)x,j,,(x,, Now using

points

model - lht)

x4, gives

k-IAt)+Kij,,(x,)~(k-IAt).

(26), (27), and (28), the ith element Xi(x, 2, = C C i 1

(29)

of X(x, z) is found

Xij,lCXT

to be (30)

z)Y

where j extends over all columns of G(x, z). Equations (29) and (30) provide the required algebraic form of the DPS in hand. The key-point of this method is the derivation of the transfer-matrix element (27). An example of such a derivation is given in (8) for the heat-conduction system studied in (5) by the Laplace transform method. In this example (see Section 11.3.) we have

T!;y:_: , Ki(X)

Gb, z) = j.

=

Ai(x)(l-ai)

where ai = exp (-niAt)>

no = K2,

rri = p;,

K2 = l/y,

(31)

’

7-ri

1

pi

are

the

roots

of

/3 tan p = a, and K’cos(l-x)K A,(x)=

Ai

K cos K--sin CY

’

K

2K2 cos [( 1 - x)pi]

=

(K’-P:)[++-j The corresponding

transfer

function

cos pi ’

(in the Laplace

G(x, s) =

i=l,2,....

izo 2.

domain)

is

(314

I

Vol. 309, No. 6, June 1980 F’rinte.5in Northern Ireland

409

S. G. Tzafestas Now we turn lumped-parameter

our attention to the integral equation control model (5) or (9), i.e. X(x, t) =

x0(x, t) +

I0

model.

‘K(x, t;T)U(T)

Consider

dr.

the

(32)

IV. l.(v) Numerical integration. This is the most straightforward approach. Any numerical integration formula can be used in (32) to determine a discrete model with the required accuracy. In general, the resulting algebraic model has the form X(x,, 1,) = T i

c,K(xi,

tk,

Tj)u(Tj)+X”(Xi,

(33)

fk),

j=l

where tk = kAr; T= N,At, At = AT; and cj, j = 1,2,. . . , are coefficients depending on which numerical integration formula is employed. For example, if the Simpson formula is used, then c,, = ck = 1/3Nr, . . = ck_ 1= 413 NT,

c*=c3=.

(34)

. . . = ck = 2/3N,.

cz=c4=

IV.l.(ui) Walsh-Galerkin expansion. If the objective function involves LJ(t> only [i.e. it does not involve the control rate U(t)] expand U(t) in Walsh series. If o(t) is also involved, start by expanding ri(r) (27). Thus, for example, on setting U,(t) = v:@(t), or in matrix form

U(t) =

L -. 1 @T(f) @>‘(t)

(32) takes X(x, t) =

2,

=@(t)v,

(35)

K(x, t; T)@(T) dT,

(36)

. aT(t)

0

the integral

0

the algebraic

form

X0(x, t>+T(x, t)v,

Iyx, t) =

where u is a column vector composed by the column Walsh coefficient vectors u,, i = 1,2, . . . , m. IV.l.(uii) Time discretization. Applying the model (32) in discrete time, with initial time kAT and final time (k + l)Ar, yields X(x, k + 1) =

fD

K(x, k + 1,5, k)X(&

k) d5+

K(x, k + 1,

T)

dT U(k), (364

which, if discretized numerical formula),

in space gives

X(X,, k + 1) = f

as well (computing

CjK(Xi,

the spatial

k + 1, Xj)X(Xj, k)+B(xi,

integral

k + l)U(k)

by some

(36b)

j=l

for xi, xj ED

410

with

i, j, = 1,2, . . . , ND.

lou,,,al of The

Franklin Institute Pergamon Press Ltd.

Distributed-Parameter Then,

using

(36b)

Optimal Control Via Mathematical

at all available

points

xi, i = 1,2, . . . , N,, k=O,l,...,

X k+l=A(k+l,k)Xk+B(k+l)U(k),

Programming we obtain

NT-l,

(37)

where

1 ,

‘c,K(x,, A(k+l,

k + 1, xi). . . c&W,,

k + 1, xi.+,)

k)=

(374 .c,K(x,~,

Finally, desired

B(k+l)=

k + 1, xi). . . qv,Wx~D, k + 1, XND)

the discrete-time model (37) at k = 0, 1, . . . , NT-i model XNT=AXo+Bv,

applying algebraic

yields the (38)

where A = A(N,, B=[A(N,,

..A@, 1MlL AU%-, NT -1)

NT- 1). . . A(2, l)A(l, NT--l).

01,

. . . A(35 2)B(2),

. . . , AN-, N-PWT-J,

BW,)I,

1I (384

.

IV.2.

Reduction of objective functions to mathematical programming form

To reduce the objective function of a DP control problem to the algebraic form required by MP, one must apply a method compatible with that used for the reduction of the system model to the desired algebraic form. The two general methods are: (i) Space-time discretization. (ii) Assumed mode expansion. IV.2.(i) Space-time discretization. By this method integrals are converted into sums, with coefficients depending on the numerical integration formula used. The various objective functions listed in Section III take the following forms. Energy-like objective function [see Eq. (17a)]: L = T k=l

CI,k[,z

G,JT(Xi,,)Q(Xj~

k)X(+,

k)+

UT(k)R(k)U(k)],

(3%

integration coefficients, where cr,k and cZ,, are the time and space numerical respectively. This form is appropriate for utilization with the discrete model (33). If the discrete model (29)/(30), or (38) used, the control summation vol. 309, NO. 6, June 1980 Printed in Northern Ireland

411

S. G. Tzafestas be C$3 cI.k+l U=(k)R(k)U(k) and U(l), . . . , U(N, - 1). The objective function (39a) can be written

must

J=

T

Cl&C&

include

the

in the compact

sequence form

+ UT&

(39b)

k=l

where

I% and v are defined

as in (37a) and (38a), and 0

‘%Q(xl,k) ok =

(39c) 1. 'c,,,Qh,, k)

0

L Quadratic

U(O),

terminal objective function [See Eq. (17b)]:

where z(xj,

T)

=

x(xj,

T)-xd(xj),

and

Terminal state distance objective function J=

?

%,j

C

j=l

(XCxj3

[See Eq. (17c)]: T)-xd,Cxj)(

(414

T)-X,,(xj)(.

(41b)

I

or [see Eq. (17c’)l: J=maxmax[Xi(xj, I i

Quadratic terminal state and control energy objective function [See Eq. (17d)]: J = 2

czj%“(xj,

T)Q$(xj,

T) + F

j=l _

=

c,,,UT(k)R(k)U(k)

k=l

x&,Q,X,

_

+ v=Rv,

(42)

where JirN,, of, v, and I? are defined as in (39b) and (40b). Fuel-like objective function [See Eq. 17e)]: J= k$I

Cl,k

T

(43)

IV,(k)I*

Control amplitude objective function [See Eq. (17f)]: J=maxxIU,(k)J, k i 412

k=l,2

,...,

A$.

(44) Jounml of The

Franklin fnstitute Pergallwn Press Ltd.

Distributed-Parameter

Optimal Control Via Mathematical

Terminal state distance and fuel-like

Programming

objective function [See Eq. (17g)]:

Minimum-time objective function [See Eq. (17h)]. This cost can be treated directly by the Pontryagin principle. In order to apply MP, an objective function must be constructed, which imposes a penalty when the state response has not reached its desired terminal value, the severity of the penalty increasing with increasing time. Such an objective function is

2 10k-‘Nd2’i

J=

j-

i=l

k=l

IXi(x, k)-X,,(x)1

dx

(46)

D

or

“c’

lok-‘N,O)

J=

i

In discretized

space,

% 1 c,jIX(x,,

j-

i=l

k=l

the spatial k)-Xdi(xj)]

{X,(XTk) - X,$(x)I’ dx.

integrals and

form

must be replaced T c,j{X(xj,

by the sums

k)-Xdi(xj)>‘.

j=l

j=l

Another

(47)

D

can be constructed

using

the concept

of the

attainability

set

(l&36). The absolute values involved in the objective functions (41) and (43-46) can be converted into linear forms, and so all these objective functions can be minimized by LP (i.e. the simplex algorithm). Two methods for such conversion are: the difference variables method and the bounding variables method (2-5). The difference variables method. In this method we set Xi(xj, T) - X,>(x,) = X;(xj, T) - X:(x,, T), subject

(48a)

to XI(xj, T) 20

and

Xy(xi, T) a0

for all

i, j,

(48b)

and similarly,

U,(k) = U;(k)subject

U’,‘(k),

(48~)

to

U:(k)>0

and

U:(k)*O,

where X’, X”, U’, and u” are new variables which replace and U. If Xi(xj, T)- X,)(xj) 2 0, we set X’/(x,, T) = 0. If Xi(xj, T)-X,L(xj)
(4W the original

ones X

413

S. G. Tzafestas Similarly: If U,(k)~O, we set U’i(k)=O. If Ui(k)
in the cost functions

IXi(xi, T)-Xdi(xj)\

= Xl(x,, T)+X:‘(x,>

can be replaced

by (49a)

T)

and IU,(k)l=

U;(k)+

U’;(k).

(49b)

Obviously, one now has double the number of variables, though no additional constraint is introduced. The bounding variables method. In this method we use the new variables zi(xj, T) and V,(k), which satisfy the following inequality constraints: IXi(Xj, T)-X,i(~j)J~~i(Xj,

(504

T)

and

(5Ob)

I Ui(k)I s Vi(k), or, equivalently, T),

Xi(Xj, T)-XdG(~i)2-Zi(~i,

T)aXi(xj,

Ei(Xj,

T)-X,,(X~),

(5 14

and U,(k)a-Vi(k), Thus,

the objective

function

Vi(k)2

(45), for example,

(5 lb)

becomes

ND J = 1 c,,~ c Ei(xj, T)+ k=l 2 j=l

U,(k).

Cl,k;

(52)

V,(k),

I

which is linear with respect to the new variables Bi and Vi. Obviously, we again have double the number of variables but now we also have two sets of additional constraints. Hence, the first method is usually preferred. IV.2.(ii) Assumed-mode expansion (Walsh-Galerkin). In this case we use the expansions U(r) = a,(t) v U(X, t) = ~

[see (35)],

dint

Vi = aT(X,

t)V,

i=l

X(x,

t) =

r

T’,(x)@=(t)[P’c,

+

CP]

=

aT(x,

t)[P*TC

+

CO],

1

i=l

where

aT(X,t) = [*,(x)+=(t), . . . >*ND(X)@T(f)l, vT=[vT, VT,. . . v&,1, )

P” = diag [P, <, . : . ,,P].

414

I

(54)

Distributed-Parameter

Optimal Control Via Mathematical

Programming

Also xd(x)

= ?

&f,Ti(x)

i=l

=

z

cT;,@(t)qi(x)

=

z

qi(xPT(T)cd,

i=l =

aT(x,

t)c,

=

aT(x,

T)c,,

where

with

cf=[5d,, 0, f . . 901. Thus,

the objective

function

(17a),

for example,

Wb) becomes

T J=

I0 [I

= +

(P*=c + EO)=a(x, t)Q(x, t)a=(x, t)(P*=c + Z”) dx D

1

dt

v’@=(t)R(t)@(t)u

(P*=c + ~“)TQo(P*T~ + Z”) + v=R,u,

where

(564

T

Qo= R, = EO=

a(x, t)Q(x, t)a=(x, t) dx dt,

H0 10

D

=@=(r)R(t)@(t)

dt,

(56b)

CO-C,.

1

The objective function (56a) should be used in conjuction (25a). Alternatively, one may directly use the Walsh-Galerkin the objective function (17a). IV.3.

with the model model (36) in

Treatment of constraints

Suppose

we have the constraints X’(x, t)aO,

[see (48b) and (48d)]

X”(X, t)s0,

If we use the discretization

method,

X’(X~, k)~O,

LJ’(t)aO,

these

X”(Xi, k)~O,

for all i = 1,2, . . . , ND and If we use the Walsh-Galerkin

constraints U’(k)~O,

U”(t)z=O. are replaced U”(k)~O

(57) by (58)

k = 1,2, . . , I$. method, we obtain

X’(x, t) =

cYT(X,t)[P*=C’+C0’]30,

X”(X, t) = aT(x, t)[P*TC”+ CO”]3 0, (59)

U(t) =Gqt) u’z=O, u”(t)=@(t) Vol. 309, No. 6, June 1980 FVinted in Northern Ireland

u”S0.

i

415

S. G. Tzafestas In general, X(x, t) and U(t) in the constraints must be replaced by their Walsh-Galerkin expansions. For example, (48a) gives the following constraint for the Walsh coefficient vectors: c-c, which,

if substituted

in (25a),

= C’-cCr’,

(60)

gives

c’-c”=T(U-v”)+y-c,. Using the above equations one can treat the absolute values of the variables. IV.4.

Treatment

of objective

functions

(61) all objective

involving

We shall illustrate the treatment of objective control with a stochastic system in the integral The objective function is

functions

which involve

the rate of control functions involving the rate of representation (16).

J=E[/OTL(l)dl],

L(t) =; J {x=(x, OQob, D

o-vx,t)+xT(x,

(624

OQ,(x,

w,(x, t)

+ U=(x, t)Ro(x, t)U(x, t)+ CJT(x, t)R,(x, t)Q(x, The pointwise

control

t)} dx.

(62b)

case

2 tJ,(t)t3(x-x,)

U(x, t)=

(63)

k=l

is assumed.

Then ri,(t)=~=(t)uk,ir(t)=~(t)u, U(0) = Q,(t) u”,

U(t) = aqt)(Pu •t UO), where a(t) Introducing

= diag [QT( t), a,‘(t), . . . , Ca’( t)]. these relations into the model (16) gives X,(x, t) = x(x,

t) + A,(x,

X(x, t) = X0(x, t) + A,(x, where A(x, t), A,(x, is the operator

t) 2)+ a(x, t),

N(0,JJt 0

(65)

t) 2, + b(x, t) + N(0, r) W(x, t),

t), a(x, t), and b(x, t) have obvious

t) =

Using

(64)

K(x,

t

definitions,

and NO,

; 5, 7)L.l d5 d7.

t)

(66)

D

(63) and (64) in (62b) yields

J

X=(X,

t)Qo(x,

t)X(x,

D

+ U=(O)R,(t)@(t)Pu

t) dx +;

X:(x,

t)Ql(x,

t)Xt(x, t) dx

1

+;

u=fi(t)v. JournalofTheFranklin Institute

416

Pergam0n Press Ltd.

Disttibuted-Parameter where

Optimal Control Via Mathematical

&(x1, 0.0 , R,(t)= R,(d = 0 kAxm,, 1 [ t) R(t)=

%(x1,

P=
0

0

0

Thus,

Programming

the objective

function

(62a) finally

takes

&(x,,,

the algebraic

t) 3.

form [see (65)]

J=uTSv+%+q

(67)

where T

S=iE[I (s(A:Q,A,+A:QIA,)dx+fi(t) C=E [/pi/ 1

IV.5.

dt

I

,

;(X”+b)TQoA,+(~+~)QIA~dx]+UT(0)Ro(~)@(~)P]dl],

0

(T=-E 2

I

D

= {(x”

+

b)‘Q,(X” + b) + (X, + a,‘Q,(x

+ cx)+ WTNTQONW}

1

dx d t .

H OD i

Formulation of some distributed-parameter

system design problems

IV. 5 .(i) Problem (i) . It is desired to design a DPS which is physically realizable with its impulse response following a desired output Xd(x, t), x E D, t E [0, T], as closely as possible. Let G(x, s) be the transfer function of the system which (without loss of generality) can be assumed to have the form

G(x, s)= with n sufficiently to minimize the function is

,zo z

large (5, 8) [see (31a)]. Then G(x, s) is to be designed so as error 2:(x, t) = X(x, t) - Xd(x, t). An appropriate objective

xi k=l

E

D,

j=l

where

2(X, tk) =

i

Ai

eem~'k-Xd(xj, tk).

i=O

Clearly, J has to be minimized with respect to the Ai terms and/or the q terms. Constraints may be placed on the poles q to limit the rise-time and the plant bandwidth. Constraints on Ai can be incorporated to take into account stability, realizability, and time-domain performance to specified inputs. IV.5. (ii) Problem (ii). It is desired to design a digital series compensator G,(z) so that its response to a specific input U(z) follows a desired output Xd(x, z) Vol. 309, No. 6, June 198G Printed in Northern Ireland

417

S. G. Tzafestas as closely

as possible.

This is similar

to the previous

problem.

Here

we have

X(x, z) = G(x, z)G,(z) U(z),

X(x, 2) =X(x, z)-Xd(X, z). One

must

choose

G,(z)

so as to minimize

k=l

the error

function

j=l

or J= The system

pulse

transfer

2

2

k=l

i=I

function

xT(xi,

kAht)A(x,

kht).

is of the form

G(x, z) = i

K(x)H,(z),

i=O

where Hi(z) = z-‘/(1 An example

- a,z-‘).

is [see (31)] Ki(X) =

4(x)(1 - ai)

with

ai = exp (-rfiAt)_

=i

One can also treat closed-loop configurations with forward transfer G(x, z)G,(z), and the case of stochastic inputs utilizing expected objective functions (e.g. mean-square error, etc.). The case in which the input signal U(z) has to be selected so as to minimize .I can also be treated; in this instance one puts U(z) = u,+ uiz-l + . . . + upzmp, where the ui terms are constant parameters to be selected. IV.6.

Solution procedures by mathematical programming

Having reduced the system equations and the objective functions in algebraic form, one can apply linear programming (LP), or quadratic programming (QP), or nonlinear programming (NLP) as the case may be. Thus, for example, the objective functions (41a) and (43)-(46) can be treated by LP if the DPS is given in the form (29)/(30), (33), or (37), and all other constraints are linear. Under the same conditions, the objective functions (39b), (40a), (42), and (47) can be minimized by QP which is reduced to LP(4). An analytical solution is possible in these circumstances if no inequality constraints are imposed. For example, the optimal value of u for the objective function (67) is the solution of _$=2Su+ZT=0,

i.e. uopt= -;

Similarly, in the unconstrained case the objective given by (25a) can be minimized analytically.

418

S_‘CT.

function

(56a) with c being

Journal of The Franklin Institute Perganml Press Ltd.

.

Distributed-Parameter

Optimal Control Via Mathematical

Programming

It is useful to remark here that, when discretizing the time of a system, the time intervals may not be equal and, more importantly, may not be known a priori; the time intervals are then part of the variables to be selected by the optimization procedure. This is best illustrated by using the discretized integral model (33), i.e.

cjK(xi, fk, Tj)u(Tj)

X(X,, tk) = T i

+p(Xi, ‘k)>

j=l

where Ik= i

At,

and

7i = $, Ati

i=l

are unknown constraints

and

;

must

be

At, s T,,,

specified.

and

Usually,

one

applies

the

additional

A& 2 Tmin (say lop3 set).

i=l

In all cases this problem leads to a NLP problem. In a similar way, one may include the discretized positions xi among the variables to be selected.

V. Review

of Nonlinear

Programming

Methods Applied

to DPS

Studies of nonlinear DP control may be grouped according to the method employed: 1. Gradient methods (10, 18, 37-40). 2. Conjugate gradient methods (41, 42). methods (43, 44). 3. Quasilinearization 4. Penalty methods (45-48). 5. Second-order methods (14, 17, 48, 49). 6. Functional analysis methods (16, 18, 48, 50, 51). In (10) an approximately optimal feedback DP controller is designed which is based on the optimal open-loop controller found by the simplex method of LP. The feedback controller is determined by a gradient procedure. In (18) a steepest-descent algorithm is given for a terminal DP control problem, and also a computational scheme is developed for time-optimal control of DPS with multi-norm constraints. In (37-39) the gradient technique is utilized in various nonlinear DP control problems, and in (40) a steepest-descent technique is developed for second-order nonlinear DPS. In (41, 42) the conjugate gradient technique is applied to linear DPS. In (43) the quasilinearization technique of Bellman and Kalaba is used for solving nonlinear DP TPBVB, and in (44) a direct search of objective function, combined with quasilinearization, is applied. In (45, 46) some results concerning the application of the penalty method to DPS are given, in (47) the conjugate gradient method is used for penalized convex programming in Hilbert space, and in (48) a second-order Newton-like Vol. 309, No. 6, Juna 1980 Printed in Northern Ireland

419

S. G. Tzafestas algorithm is provided for general systems in Hilbert space, in which the inequality constraints are treated by the penalty method. In (14) a secondorder optimization algorithm for nonlinear DPS is given which is based on differential dynamic programming (iterative solution of Hamilton-JacobiBellman equation). In (17) a second-order algorithm is developed for nonlinear DPS, based on a successive application of LQP, and in (49) a second-order method is put forward for nonlinear DPS with nonlinear boundary conditions, based on a second-order variational expansion of the objective function. In this method, which is an extension of Mitter’s technique for lumped-parameter systems (52), the solution of the DP adjoint equations is not required. A unified presentation of the methods in (14) and (49) may be found in (53). In (16) the Newton-like algorithm is applied to a class of nonlinear DP control problems with the aid of Green’s identity, in (50) a treatment of DP control problems using function space techniques is given, and in (51) a modified steepest-descent method in Hilbert space is devised, which possesses better convergence properties than conventional steepest-descent methods. Below the techniques developed in (10, 16, 18, 37, 41,46, 48) will be briefly reviewed. V.l.

Application

of gradient methods

V.l. (i) We start with the application made in (lo), where a specific DP optimal controller (SOC), i.e. a feedback controller whose input-output relation is specified in advance within some adjustable parameters, is designed. The parameters of this SOC are chosen in an optimal way by a gradient algorithm so that the output of the controller generates a “best approximation to the optimal open-loop controller”. The output measurements are assumed continuous in time, but discrete in space at a finite number of points. The system considered is the same as that used in (5), and the optimal open-loop control is computed using the simplex method by LP. Consider a SOC of the form

UcWc,h, t) = HX(xd, h, 0, k tl,

(68)

where X,(x, t) = XC(xd, h, t) is the closed-loop response, h is the parameter vector to be determined, and x, is the column vector of measurement points at which the sensors are located. The parameter vector h = [h,, hZ, . . . , h,lT is chosen so as to minimize J=

IIIXc(Xd,

Uc, h, T) - -X,,,tb, T) I dx

(69)

0

which expresses the accumulated difference Xopt(x, T) is the optimal open-loop response. Clearly, the first differential of J is 6J = pTAh

pT = [PI,

in X(x, t) over

~2,

. . . ,

piA

D = [0, 11. Here

(704

Distributed-Parameter

Optimal Control Via Mathematical

Programming

where pm = aJ/ah, =

ay, -=

T

g(x,

T-r)

Io

ak

Then, with Ah of constant respect to h yields

(7Ob)

I0 $dr m

(704

[see (5)].

step size ((Ah\\ = (Ah=Ah)i,

the minimization

of J with

(71)

Ah = -(]lAhlllllp]])p. Hence,

the updating

algorithm

is

h new= hold+ (Ah), which

improves

the value

of J, i.e. makes

(72)

&I negative,

6J= p=Ah = -PT;“;h”= P

thus:

-]]p]].]]Ah]]
The steps for computing U,(X,, h, t) then are: (a) Compute the optimal open-loop response and make initial guesses h and the sensor position vector xd. (b) Using h and xd, compute X, and J. (c) Update h using (72), and then repeat steps (b) and (c) until the minimum of J is obtained to a desired accuracy. The value of h corresponding to this final .I is hopt. The desired SOC is found by using h,,, in (68). The same procedure can also be used for choosing optimal measurement locations, i.e. for choosing both h and xd. V.1. (ii) We now consider the steepest-descent algorithm given in (18). The system is assumed to have the integral form (2), which is written as X(x, t) = X0(x, t) + 2(x, where

f 6p(x, [I0 JI, t) =

u(x,

The problem

t, = ["d(x,

t) U(x, t),

K(x, t; 5, 7X.1d5 dt,

K,(x, t; &, 7X.1d& dt , I

t), u,,(x, t)lT.

is to minimize J = 11X(x, T) - Xd(x)ll+

where

FZ? is the Hilbert

space

in which X(x, T) and X,(x) belong, subject space in which U(x, t) belongs. Minimizing the constraint, is equivalent to minimizing

IIuJI&GE’, where H is the Hilbert with respect

to U, without l(U) =(A(%

Vol.309,No.6,June1980 PrintedinNorthemIreland

VU,

&r-2(2*(%

T)&(x),

to J

U>,, 421

S. G. Tzafestas where gd(x)= Xd(x)-x0(x, T), L&*(x,T) is the A(x, T) = 2*(x, T)Y(x, T). The solution is obtained &91/19U=A(x, In the steepest-descent U, and seek a correction

T)U(x, t)-2*(x,

adjoint of by solving

9(x,

T),

and

T)&(x)=O.

method, we choose an initial AU’ E H such that

guess U1 of the optimal

a1( U1 + y1A U’)/ 13~1

(73)

is maximized for y1 = 0. Then we choose y1 such that I(U’+ ylAU1) is minimized. If this value of -yl is designated yypt, the next approximation is U” = U’ + yypfAU1, and so on. Here,

condition

(73) gives T)U’-2*(x,

AUApt=A(x,

of I( U1 + ylAU’),

and the minimization

YlOpt=

Hence,

with respect

-bU~,h.x/(Ab,

the computational

procedure

= Nx,

to y1 gives

T)AU:pt,AU:p,)~.

is defined

Cl:,:’ = Uzpt+ yT’AUgpt, Aup,,

T)X,(x),

T)U:,,-2*(x,

yy* = -[iA U:,,llfIl(Nx,

by the relations

VA,, given, T)%(x), T)A UEpt, A UEpt>H,

where p determines the iteration number. In order that U satisfies the constraint UE %, where the above algorithm is modified to E;:(p) Avp,,(p) $%4

= U,“,&)

+ yiP%)A

UZp&),

= (A + ccW:.,&)-2*(x, = - IlAU:,&)II”l((A

Ou = { UE H:ljUl~~ E’},

UA,, given,

T)&(x), + CLQAU$,,

AU,P,,)W

p+l = U$l(O). Otherwise p must be chosen If U:,:‘(O) E %, one must choose U,,, to be pp, such that Uf$:(pp) E 94 i.e. IIU~,:1(p,)Ij&=E2 and l-J:,:‘= Ui,:‘(pp). In practice, the iteration is stopped when the ratio 77= IlAG,,

- s*%

is less than a specified value qO. V.l. (iii) The results of (37) concern a2x

-g=F,

ak-lx

422

a nonlinear

IIH hyperbolic

DPS of the type

X,s,+;,u,x,t> ,

X(x, 0) = &o(x), = CYl, axk-l1l?D

IIJll~*X,

E&J!_)=p&(x),

1

ak-2x axk-2 4D

=(Yz,

ax -

. . . , ax

1

I

aD=ak-l~

(744

I

Journal of The Franklin Institute Pergamon Press Lid.

Distributed-Parameter with objective

Optimal Control Via Mathematical

function

X,g,?$U,x,t)dxdr+[b Since the system form as

ax,/at

dk-‘&/dxk-‘]~D and the objective

= X2(x, t),

= (Yl,

function

is expressed

of the problem

in state-space

X,(x, 0) = X00(x), X2(x, 0) =X,,(x),

. . . , ~x&x],,

X,(x, t) The Hamiltonian

O(X,~cix]~. (74b)

in t, Eqs. (74a) are first written

is second-order

ax*/ar = Fd(. . .),

= (Y&l,

in terms

and

of the state variables

X2(x, t).

is

H = L + hT(X, t)XJx, where

Programming

A and p are DP Lagrange

t) + /_LT(X,f)Fd(. . .),

multipliers.

Thus,

we obtain

(H-hTX,-&d)dxdf+

J=

The system and objective function are discretized in space, and a model of the lumped-parameter type dX/dt = F(X, U, t), X(0) =X0 [see (20)] is obtained. The state vector is then augmented by a component X0(t), where dX,/dt = L, X,(O) = 0. Now J = X,(T). The necessary condition of optimality is 6J = 0, where SJ is the first variation of J. Here,

where A;(t) = h=(t)[aF/aU]. The greatest

6J is obtained

when au(t)

where usual,

R(t)

is a positive

definite

G.,(t)

= -R(t)/+,(t), matrix

= U,,,(0

function.

- WW”(0.

(75) Thus,

the algorithm

is, as

(76)

V.l. (iv) Lastly, we consider the second-order (Newton-like) gradient algorithm of (16). By using standard variational theory, the canonical equations Vol. 309, No. 6, June 1980 F’rint.4 in Northern Ireland

423

S. G. Tzafestas below

which

form a TPBVP

are derived:

X = NdW, X,, U,, x, 0, X(x, 0) = X,(x), Nb (X, X,, U,, x, t) = 0, % = fd Wd, Ud, 0,

% =fbK, --A, =

x

E

x E 0,

IYD,

(774

X,(O) = xfo,

Ub,t), X,(O) = XbO,

(NJ*,,

(N&A(x,

1

a J/ax, x E D,

A (x, T) = t) = 0,

x

E

IYD,

The difficulty in solving these equations is that X(x, t), X,(t), and X,(t) are given at t = 0, whereas A (x, t), Ad(f), and Ah(t) are given at t = T The algorithm, which is actually a hill-climbing on A(+ 0), ~~(0), and pb(0), involves the following steps: (a> Make a guess for ydx) = y(x, t>L, = co1 [A(x, t), kd(t), pb(t)ltcO. (b) Using X(x, 0), X,(O), X,(O), and y,,(x), integrate the canonical equations in forward time from t = 0 up to t = T, so as to compute the error Ayf(x)= [y(x, t) - Y~(x)]~=~, x E D, where yf(x) is the true final value of y(x, t). (c) Integrate the sensitivity equations for

R(x, Xl, 0 =

2,

Q(x,x~,

adx, 0 aYoh) ’

t)= -

where S(x, t) = co1 [X(x, t), Xd(t), Xb(t)] is the state vector. equations are found by differentiating the canonical equations Y,,(X), x E D. (d) Improve

(784

x,x,~D

The sensitivity with respect to

y,,(x) by using c

J

Yo(X)new = YO(XL,d - 8

D

0-k

XI)A~~(XI) dxl

G’8b)

for some 0
424

Journal of The

Franklin Institute Pergamon Press Ltd.

Distributed-Parameter

Optimal Control Via Mathematical

Programming

V.2. Conjugate gradient method The conjugate gradient method is a second-generation gradient method having quadratic convergence properties. It may be summarized as follows: (a) Guess the initial control function Up for p = 0. (b) Integrate the state equations forwards from t = 0 up to t = T, and the co-state equations backwards from t = T to t = 0. (c) Compute the gradient gp = g(Up) of the objective function, and the conjugate gradient parameter Op, where W = (g”, gP)/(gP_l, with 0’ = 0. (d) Compute the improve Up by

direction

RP of search,

UP+ePRP,

u p+l=

(794

gp-‘)

where

p=O, 1,2,

RP = -gp + OPRPel, and ..>

U’9b)

where sp is found by a one-dimensional minimization along the direction search. The algorithm is terminated when a specified accuracy criterion satisfied. In (41) the inner

product

(. , .) means

notation

UT(x, t) Wx, t) dx dt,

(U(x, 0, V(x, t)>D = j T j 0

of is

XED,

D

T (“(x,,

t),

The DP model

v(x,,

t)>SD

=

considered

0

H

UT(xb, t)V(x,,

pX(x,, t) +&L&(x, t) = 0, function

L = The gradient

(X’Q,X+

U,TR,U,) dx +

of J is

g(U)=[-~~-~~;where

the co-state

variable

t),

XE D,

xb E IYD,

T) dx , 1

(81)

(XTQbX + U,TR, U,) dx,.

(82) I[--___________---1) (RdUd +B,Th) XE D = (RbUb+B;A),xEfiD

h(x, t) is described

by

4(x, t) = -Z*A(x, t) - Qd(x, t)X(x, t), A (x, T) = Q,XCx, T). Vol. 309,No. 6,June 1980 PrintedinNorthem Ireland

(80)

is

XT(x, T)QfX(x,

J=;{joTLdt+jD

x,, EII~D.

is

Xt(x, t)=,SPX(x, t)+BdUd(x,

and the objective

t) dx, dt,

6D

(83)

425

S. G. Tzafestas V.3. Quasilinearization

method

Quasilinearization is a technique with quadratic covergence and monotonicity properties, by which a nonlinear multi-point BVP is converted into a more readily solvable linear BVP. It is actually a modification of the NewtonRaphson method (54, 55). In (37, 43) the DPS is discretized by finite differences and transformed into the lumped-parameter form, which is treated by the quasilinearization technique as follows: Consider the system k(t) = F[X(t)], X(0) = X,,. Let X”(t) be the initial guess of the solution. Linearization of this dynamic equation about XP(t) gives kp+l=

F(XP)+D{F(XP)}[Xp+l(t)-XP(f)],

where D = [Dii] = [IYFJI~X~]is the Jacobian tal matrix of

and BP”(t)

is the particular gP+l=

then

the complete

vector

solution

to (84) is given

xp+l(t) In a TPBVP

of optimal

Cap+‘(O)= I,

solution

F(XP) + D{F(XP)}(BP+l

If Qp+l(t) is the fundamen-

matrix.

6)“+‘(t) = D{F(XP)jW+‘(t),

(84)

of - SP),

B”+‘(O) = 0,

by

=~~+l(f)X(0)+E~+l(f).

control, X(t)=

+_

[

X(t) is given

I[ =

(85)

by

State vector _______~__‘_______

AdJoInt

vector

1 )

with X(0) and A(T) given. The initial value A(0) required for the adjoint vector can be found using (85). An initial guess X”(t) of the solution is provided by the method of gradients. The computational solution through quasilinearization converges very quickly. V.4. Penalty methods We shall describe two such methods, namely, a method based on (16) and (48), and the method presented in (46). V.4. (i) Consider the TPBVP (77a, b) with the additional inequality constraint G(X, X,, X,)
426

coefficients,

X,, X,)H,

and F{.} is a vector-valued

penalty

Journal of The Franklin Institute Pergamon Pres Ltd.

Distributed-Parameter where

fi (.) are continuous

Now define

Optimal Control Via Mathematical

non-decreasing

functionals,

Qi = 1

if

G,(X, X,, X,) > 0,

Q = 0

if

G,(X, X,, X,)
the augmented

vector

Programming

and

functions

where @ is the vector, whose components are the final conditions written so that the ith component ai = 0 when the ith condition is the vector with components

of A, /A~, etc., is satisfied; s

[ X.X&Xb 1,

Zi = Ci Qi max {Gi(X, X,, X,) - Ei}

maximum penetrations where Ei, i= 1,2,. . . , q, are the permissible forbidden region, and C > 1 are constants which are chosen according importance of Zi. In the present case the updated Eq. (78b) becomes yf;”

= yg-

&

J[$p(x,

x1,

in the to the

t)]-'~Pdxl,

D

where C&x, x1, t) = X%x, t)/+?(x,, t). If the inequality constraints also involve the control variable, then one must first convert the problem to one of terminal control by the technique of increasing the state vector [see section V.l.(iii)]. The computation of the penalty of sensicoefficient CL vector along ho, A,,, hb,, makes possible reduction tivities; reduction of computation time; and determination of the maximum trajectory penetration in the forbidden region. V.4. (ii) The work in (46) concerns a nonlinear discrete-time DPS of the type (ll), i.e. X(x, k + I) = X(x, k)+F,,{X(x, with X(x, 0) = X,(x)

k), X,(x,

and X(x, N) = X,(x);

k), Xx,(x, k), U(x, k))

the objective

function

is

Jk[Wx, k), Wx, k)l dx. The problem is to choose U(x, k) E 92, k = 0, 1, . . . , N- 1, so as to the solution is found by the ~-technique which in (56) is applied to time DPS. The s-problem here is in fact to choose both U(x, k) E Q and [where % and Z are assumed to have special properties (46)] so as

minimize J: continuousX(x, k) E 2 to minimize

N-l

J,=J+

c k=OE

and satisfy

all initial

Vol. 309, No. 6, June 1980 Print.4 in Northern Ireland

i J /X(x, k+l)-X(x,

k)-F,#dx

D

and final conditions.

427

S. G. Tzafestas Define h(s) as h(s) = min J, aJE(X”(x, 0, E), v(x, 0, E)), where p and Uc’ are the optimal values and the minimization is performed over all X(x, k, E) E Z and U(x, k, E) E Ou with X(x, 0, E) =X,(x) and X(x, N, E) = X,(x). Now h(e)< inf J(X(x, .), U(x, .)), U(x, .) E ‘II. If e1
P(E)=~~~ j-IIX% k+1, k=O D

E) -

x0(x, k, E) - Nk (x”, g,

x:,, u”)112 dx.

If &i<&2, then .4,(X% . , ~~1, uO(x, . , Ed)) ~.J,,W”(x, . , 4, vb, . , -d), where &(X0(x, . , E), U”(x, . , E)) = J(X’(x, . , E), U”(x, . , E)) + (l/e)p(e). Furthermore, there exists a sequence {ei}, with si >si+i for all i, which tends to zero as i-+m, so lim h(Ei) = lim

im

J(?‘?(X,

. , Ei),

v(X,

. , Ei))

i-xa

=

inf J(X(x,

.), U(x, .), U(x, .)) E %.

Mild sufficient conditions for the existence of a solution to the c-problem are also given in (46), and a sub-optimal computational scheme (with the reduced dimensionality requirements of the e-problem) is given for the subclass of nonlinear DPS with

Eik(X, Xx, Xx,, U) = Nk (X(x, k)) + h(x) u(x, k). This scheme involves the following steps: (a) Construct the appropriately restricted subspaces %’ and %’ of Z and W (b) For each value of E, determine a sub-optimal solution x0(x, k, E), in aP’ and Q’, respectively. U”(x, k, E), k = 0, 1,2,. . . , N- 1, contained (c) Determine SV(x, k, E) using x0(x,

k+l,

E)=~(x,

k,

&)+Nk(*(x,

k, s))+&(x)[I?(x,

k,

k, E)].

E)+~@(x,

. E)) and compare the values of (d) Compute J(p(x, . E), u(x, . , E) + 6u(x, J corresponding to the various values of E. There exists some value E = E’ for which the corresponding objective function is a minimum. Then x0(x, k, E’) and v(x, k, so) + H?‘(x, k, E’) provide the sub-optimal solution. Some features of this scheme are as follows: (A) It is different from general finite difference schemes where the optimal control is found for the resulting discretized system. (B) No discrete approximations to the system equations are required, but only solutions in the restricted subspaces 85’ and %’ in which the search is made for the solution of the c-problem for each fixed E. (C) Each s-solution, thus obtained, is perturbed so as to satisfy the system model, and the corresponding values of the objective function are compared. (D) Finally, the perturbed solutions are extended in a small neighbourhood of the restricted subspaces to give a satisfactory sub-optimal solution of the original problem.

428

Journal of The

Franklin Institute Pergam0n press Ltd.

Distributed-Parameter VI. Application VI.1. Example

Optimal Control Via Mathematical

Programming

Examples 1. Linear programming

with integral model (5)

and parameters The system modelled by (4), with zero initial condition The model (4) was a! = 10, y = 0.04, X,(x) = 0.2, o==x c 1, was considered. transformed into the integral form (5) as shown in Section 11.3, and then it was discretized in space and time using Simpson’s composite formula with N= 20 [see (33), (34)]. The terminal state distance objective function (17~) with 0~ U(t) G 1 was considered and discretized using the same formula with ND = 20 [see (41a)]. A computer code for simplex was employed. The results are shown in Figs. 1 and 2. One can see that the temperature response distribution obtained with the optimal control is acceptable when T= 0.2, and coincides almost exactly with the desired X,(x) = 0.2, Osx < 1, when T = 0.4.

u @I4 1.0. -

Q5.

T=Q4

n Q24

6 FIG.

01

FIG.

2. Temperature

Vol. 309, No. 6. June 1980 Printed in Northern Ireland

1. Optimal

I

0

response.

0.2

control

function

Q32

t

for T = 0.4.

I

0.4

T--0.4

,

0.6

0.0

1

(a) Constant control: U(t) = 1. (Curves control. (Curves II).

X

I). (b) Optimal

429

S. G. Tzafestas VI.2. Example

2. Quadratic programming with pulse transfer function model (8)

The same system as in Example 1 was considered. First its pulse transfer function model was derived as shown in Section IV.l. [see (26)-(31)]. The discrete quadratic objective function l+N

J=

c c X(Xi, ky,~=X-Xd,

k=I

O
jcJ

was assumed for minimization. The space and Ax = 0.1 and At = 0.1. The simulation time for than 10 set on a CDC 6400 computer. Better results might be expected for a larger increasing T further towards the transient time performance (Fig. 3). For a time varying system, change of the parameters.

time increments taken were 30 sampling periods was less final time t = NTAf; actually, does not improve the system NT will depend on the rate of

X (0.5,kAt) t

0 FIG.

3.

Optimal temperature

4

8

12

16

2o

kAt

response. Curve (A): iVT= 0; curve (B): NT = 2; curve (C): NT=9.

In the present example the optimal control sequence was always less than one, and hence the analytic unconstrained solution (see Section IV.6) gave the same result as the simplex algorithm for QP. V1.3. Example

3. Multi-criteria

linear programming

control (9)

The same system as above was considered in the integral representation. The problem of simultaneously minimizing the final distance and the fuel consumption objective functions, i.e.

was solved by applying the space-time discretization technique with NT = ND = 20, Ax = l/20, and At = l/100. Table I shows the 39 non-dominated extreme solutions. The individual minimum of J, and of J2 correspond to solutions 1 and 39, respectively. Solution 1 coincides with the result of Example 1.

430

Journal of The Franklin Institute Pergamon Press Ltd.

Distributed-Parameter

Optimal Control Via Mathematical

Programming

TABLE I

1 2 3 4 5

6 I 8 9 10 11 12 13 14 15 16 17 18 19 20

0.0372662 0.0374765 0.0375273 0.0378839 0.0403418 0.0419588 0.0427904 0.0483905 0.0488581 0.0518858 0.0526976 0.0594691 0.0609402 0.0650078 0.0656572 0.0673165 0.0725064 0.0726738 0.0754747 0.0781380

VI.4. Example

JI

JZ

Jl

0.2558615 0.2520715 0.2516105 0.2501709 0.2430849 0.2387172 0.2365142 0.2219092 0.2207082 0.2130694 0.2110556 0.1945298 0.1911297 0.1819015 0.1804548 0.1768045 0.1661544 0.1658211 0.1603864 0.1552476

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

J2

0.0816508 0.0873010 0.0900488 0.0917818 0.0927104 0.0969318 0.1015699 0.1024379 0.1055101 0.1077262 0.1233609 0.1192414 0.1216673 0.1238002 0.1246673 0.1255254 0.1549703 0.1699552 0.1999995

0.1486729 0.1384733 0.1336514 0.1306560 0.1291062 0.1223247 0.1151897 0.1139017 0.1094980 0.1064104 0.0858845 0.0910728 0.0879959 0.0853398 0.0843067 0.0833334 0.0500001 0.0333334 0.0

4. Design of specific optimal controller (10)

Again, the same system as used in Example 1 was considered, with T= 0.3, X,(x) = 0.3, X,(x) = 0.5, 0 G x s 1, NT = 20. The terminal objective function by LP as in Example 1, and the J = I: [X(x, T) -X,(x)1 d x was minimized open-loop optimum control sequence {U,,(k)} and the system response were obtained. A SOC with only one measurement point xd =x1 = 1, i.e. at the boundary, was used, namely U,(X,,

h, t) = h’ Xc(xl, h’, t).

The accumulated difference objective function (69) was to h1 by using the gradient technique of Section V.l; Jp”’ = 0.0577 with h& = 1.23. The dependence of J, on Next a SOC with two measurement points x1 = 0 and namely U, = hlX,(x’,

minimized with respect the minimum value is h’ is shown in Fig. 4. x2 = 1 was considered,

xd, h, t) + h2X,(x2, xd, h, t).

In this case, the optimal solution is hA,,= -9.90, h;,,= 12.78, Jyp’= 0.0317. Clearly, possible constraints on U,, e.g. saturation, could be accommodated without difficulty. vol. 309,No. 6,June 1980 F’rinted in Northern

Ireland

431

S. G. Tzafestas

FIG. 4. Dependence of J1 on h’.

Lastly, .I, was minimized not only with respect to h’ and h2 but also with respect to x1 and x2 (this is an optimal selection problem of measurement points) by an incremental procedure. That is, for some initial point (x’, x*)l, .I1 is computed. Then, one parameter only is changed. If J1 decreases, another move in the same direction is made, until the minimum of J1 in this direction is reached. The same procedure is then repeated in the direction of the other variable. The result in the present case is 59’ = 0.0246 at xApt = 0.35, xzpt = 1.0 with h&= -12.06 and h&,,= 15.9. A significant reduction of .I;“’ over the case with fixed x1 and x2 has resulted. The accuracy of this discrete sensor location method can be improved by interpolation, or by increasing the number of gridpoints. Possible constraints on the controller parameters and on the measurement points can be accommodated. VI.5. Example

5. Nonlinear system with gradient and quasilinearization

techni-

que (37) The nonlinear

DPS considered

a2x Tg=-*x(x,

is given by

t)‘+p

X(x, 0) = xxx,, ax(x, t) ax

=-

1x=0

$+wx,t),

O
$(x7o]t=o=x,‘(x),

1

= ax(x, t) at x=L ’ o

with

432

Journal of The Franklin Institute Pergamon Press Ltd.

Distributed-Parameter

Optimal Control Via Mathematical

Programming

2.0 1.6 1.2

FIG. 5. (a) Optimal

control

function.

(b) Optimal

response

trajectory.

and Xi(x)= Initially, solution Here

l+(~x,

Xi(x)=

lx

with

(Y= 1,

5=2.

the solution was found by the gradient method. This was then used as an initial guess in the quasilinearization

L = 1.5,

p = 2.0,

5 = 1.0,

T= 0.4,

q = 0.5,

p = 1.0,

In the finite difference

approximate technique.

R = 1.0.

scheme

d2X(x, t) =xi+l(t)-2Xi(t)+xi-l(t)

ax2

@xj2

for discretizing the system equations, the spatial increment used was Ax = 0.50. The fourth-order Runge-Kutta method was employed for integrating the resulting ordinary differential equations. As an initial guess of the control function U(x, t), the optimal control corresponding to the linear model [i.e. without the non-linear term -(X(x, t)‘] was used. The optimal control and optimal state trajectories are shown in Fig. 5. Values of the objective function at the various iterations are given in Table II. TABLE

Values of .I al Iteration 0 1 2 3 4 5

Vol. 309, No. 6, June 1980 Printed in Northern Ireland

J 1.53864 1.53666 1.53429 1.53222 1.53043 1.52893

II

various iterations Iteration 6 7 8 9 10 11

J 1.52771 1.52677 1.52609 1.52568 1.52553 1.52553

433

S. G. Tzafestas VI.6. Example

6. Conjugate gradient control of vibrating string (41)

The system considered

here is described by the equations

a2x a2x -g=ax”+

o
U(x, t),

1,

X(x, 0) = sin 7rx, X(0, t) = X(1, t) = 0, [ax/at],=,

= 0,

X(x, T) = X,(x) = 0,

with the objective function being U2(x,

t) dx dt.

The terminal constraint X(x, 7’) = 0 was taken into account by adding to J the penalty term P(X) = s,,/~‘X~(X, T) dx, where the penalty parameter minimize

E, was selected

arbitrarily.

Thus, one had to

J’= J+P(x) without the terminal constraint. The adjoint equations here were found to be

a% a2A -=A (4 at2 ax2'

n =0,

A(0, t) = h(1, t) = 0, WaOl,=T

= ~EAx,

T),

and the gradient of J’ is g = U(x, t)+A(x,

t).

The results with increasing E, are shown in Table III. TABLE III Results for increasing ualues of E,

Penalty parameter 2 5 50 100 500 1000

434

Objective function: J=J’-P(X)

Penetration in the forbidden region: P(X)

0.2495 0.5501 1.0938 1.1454 1.1894 1.1957

0.14859 0.05241 0.0010421 0.00027283 0.000011353 0.0000028465

Distributed-Parameter

Optimal Control Via Mathematical

Programming

VIZ. Conclusion Mathematical programming [60-651 provides useful tools for dealing with optimization and control problems of DPS. In this paper, we have attempted to give a unified and comprehensive presentation of the solution of a variety of DP control and system-design problems, which cover the majority of practical situations. Naturally, the application of MP techniques (LP, QP, NLP) leads to open-loop optimal controllers, and these controllers may be used as a basis for designing feedback controllers as well [as described, for example, in Section V.I. (i)]. In addition to the DP models studied here, there are many others: for example, DPS with moving boundaries (57), or DPS with time delays (58), etc. All these models can be treated by many of the techniques employed in this paper. It is useful to point out that MP can also be used to design optimal DP filters of the Kalman type (59), including the sensor position selection, which may be used either on their own or in conjunction with the optimal controllers when the system is stochastic or the state of the system under control is not accessible to direct measurement. Finally, it may be noted that one can also treat DP problems by generalized programming (66) which concerns the case where the matrix A [see (19)] in a LP problem involves some elements which are unknown and must be chosen together with the other control variables.

References (1)L. Zadeh

and B. Whalen, “On optimal control Control, Vol. AC-7, pp. 45-46, July 1962.

and LP”,

IRE

Trans,

autoin.

of mathematical programming techniques in optimal (2) D. Tabak, “Applications control-A survey”, IEEE Trans. autom. Control, Vol. AC-15, pp. 688-690, Dec. 1970. (3) K. Fegley et al., “Stochastic and deterministic design and control via linear and quadratic programming”, IEEE Trans. autom. Control, Vol. AC-16, pp. 759766, 1971. (4) D. Tabak and B. Kuo, “Optimal Control by Mathematical Programming”, Prentice-Hall, Englewood Cliffs, N.J., 1971. (5) Y. Sakawa, “Solution of an optimal control problem in a DPS”, IEEE Trans. autom. Control, Vol. AC-9, pp. 420-426, 1964. (6) V. Lorchirachoonkul, “Optimal sampled-data control of DPS”, Ph.D. Dissertation, Montana State University, 1967. (7) V. Lorchirachoonkul and D. Pierre, “Optimal control of multivariable DPS through LP”, Proc. 1967 Joint Automatic Control Conference, pp. 702-710, 1967. (8) M. Sheirah and M. Hamza, “Optimal control of DPS”, Int. J. Control, Vol. 19, pp. 89 l-902, 1974. (9) M. Sakawa, “Solution of multicriteria control problems in certain types of linear DPS by a multicriteria simplex method”, J. Math. Analys. Applic., Vol. 64, pp. 181-188, 1978. (10)A. Koivo and P. Kruh, “On the design of approximaely optimal feed-back controllers for a DPS”, Znt. J. Control, Vol. 10, pp. 53-63, 1969. (11)E. Barnes, “Computing optimal controls in systems with distributed parameters”, IFAC Symp. on the Control of DPS, Banff. Canada, 1971. (12) J. Rosen, “Iterative graphical solution of boundary value problems using LP”, in Vol. 309, NO. 6, June 1980 Printed in Northern Ireland

435

S. G. Tzafestas “Techniques of Optimization”, (Edited by A. Balakrishnan), Academic Press, New York, pp. 183-194, 1972. (13) P. Wang and F. Tung, “Optimum control of DPS”, ASME Trans., J. bas. Engng, Vol. 86D, pp. 67-79, 1964. (14) S. Tzafestas and J. Nightingale, “Differential-dynamic-programming aDDrOaCh to optimal nonlinear DP control systems”, Proc. IEE, Vol. 116, pp. 1078-1084, 1969. (15)A. Butkovskii, “The maximum principle for optimum systems with distributed parameters”, Automn remote Control, Vol. 21, pp. 472-477, 1961. (16) S. Tzafestas, “Final-value control of nonlinear composite distributedand lumpedparameter systems”, J. Franklin Inst., Vol. 290, pp. 439-451, 1970. (17) S. Tzafestas, “Optimal DP control using classical variational theory”, Int. J. Control, Vol. 12, pp. 593-608, 1970. (18) S. Tzafestas, “DP control in function space”, J. Franklin Inst., Vol. 295, pp. 317-342, 1973. (19) A. Robinson, “A survey of optimal control of DPS”, Automatica, Vol. 7, pp. 371-388, 1971. (20) W. Ray, “Some recent applications of DPS theory-A survey”, Automatica, Vol. 14, pp. 281-287, 1978. (21) S. Aidarous and M. Ghonaimy, “A direct method for optimization of stochastic DS”, Automatica, Vol. 11, pp. 203-207, 1975. (22) S. Tzafestas and N. Chrysochoides, “Nuclear reactor control using Walsh function variational synthesis”, Nucl. Sci. Engng, Vol. 62, pp. 763-770, 1977. (23) T. Yu and J. Seinfeld, “Suboptimal control of stochastic DPS”, A.I.Ch.E.JI, Vol. 19, p. 389, 1973. (24) S. Tzafestas and J. M. Nightingale, “Optimal control of a class of linear stochastic DPS”, Proc. IEE, Vol. 115, pp. 1213-1220, 1968. “Analog/hybrid solution of PDE in the nuclear industry”, (25) R. Vichnevetsky, Simulation, Vol. 11 pp. 269-281, 1968. (26) R. Vichnevetsky, “Use of functional approximation methods in the computer solution of initial value partial differential equation problems”, IEEE Trans. Comput., Vol. C-18, pp. 499-512, 1969. (27) S. Tzafestas, “Design of distributed-parameter optimal controllers and filters via Walsh-Galerkin expansions”, Proc. IFAC Symp. on Control of DPS, Warwick University, June 1977, (Edited by S. Banks and A. Pritchard), Pergamon Press, Oxford. (28) G. McGlothin, “Optimal control of DPS with penalties on spatial derivatives of the state”, Int. J. Control, Vol. 24, pp. 445-466, 1976. (29) C. Chen and C. Hsiao, “State-space approach to Walsh series solution of linear systems”, Int. J. syst. Sci., Vol. 6, pp. 833-858, 1975. (30) S. Tzafestas, “Walsh series approach to lumped and distributed system indentification”, J. Franklin Inst., Vol. 305, pp. 199-220, 1978. (31) R. Vichnevetsky and R. Shyam, “Frequency analysis of methods of lines for hyperbolic systems”, NAM 50, Dept of Computer Science, Rutgers University, N.J., 1972. (32) R. Vichnevetsky and Y. Shieh, “On the numerical methods of lines for onedimensional water quality equations”, DCS Rept No. 20, Dept of Computer Science, Rutgers University, N.J., 1972. (33) S. Tzafestas, “Parameter estimation in DP dynamic models”, Inst. &em. Engng Symp., Series No. 35, pp. 5.43-5.50, Inst. Chem. Eng., London, 1972.

436

Distributed-Parameter

Optimal

Control

Via Mathematical

Programming

(34) S. Tzafestas,

(35) (36) (37) (38) (39) (40) (41) (42) (43) (44) (45) (46) (47)

(48) (49) (50) (51) (52) (53)

(54) (55) (56) (57)

“Advanced techniques of digital-model design and real-time simulation for direct digital control”, Int. AICA Symp. on Simulation of Complex Systems, Tokyo, August 197 1. J. Rosko, “Digital Simulation of Physical Systems”, Addison-Wesley, Reading, Mass., 1972. A. Fath, “Approximation to the time-optimal control of linear state-constrained systems,” Proc. 1968 JACC, Ann Arbor., Michigan, pp. 962-969, 1968. S. Chaudhuri, “Optimal control computational techniques for a class of nonlinear DPS”, Int. J. Control, Vol. 15, pp. 419-432, 1972. H. Lee and D. Shen, “Optimal control of a class of DPS using Gradient methods”, Proc. IEE, Vol. 116, pp. 1237-1244, 1969. J. Bansal and K. Chang, “Optimal control of dispersion type DPS with time-delay Znt. J. Control, Vol. 16, pp. 481-500, 1972. in interacting TPB conditions”, M. Denn, “Optimal boundary control for a nonlinear DPS”, Int. J. Control, Vol. 4, pp. 167-178, 1966. optimization of DPS by the conjugate D. Cornick and A. Michel, “Numerical gradient”, IEEE Trans. autom. Control, Vol. AC-17, pp. 358-362, 1972. approach to a DP optimal control D. Ball and J. R. Hewit, “An alternative problem”, ht. J. Syst. Sci., Vol. 5, pp. 309-316, 1974. A. Sage and S. Chaudhuri, “Gradient and quasilinearization computational techniques for DPS”, Int. J. Control, Vol. 6, pp. 81-98, 1967. J. Seinfeld and L. Lapidus, “Computational aspects of the optimal control of DPS”, Chem. Engng Sci., Vol. 23, pp. 1461-1483, 1968. H. Sasai, “A note on the penalty method for DPOC problems”, SIAM J. Control, Vol. 10, pp. 730-736, 1972. M. Tarng and L. Nardizzi, “A computational technique for the control of systems described by PDE“, Int. J. Control, Vol. 17, pp. 1189-1199, 1973. P. Kenneth et al., “Penalization techniques for optimal control problems governed by PDE”, in “Computational Methods in Optimization Problems”, Edited by L. Zadeh et al., Vol. 2, pp. 177-186, Academic Press, New York, 1969. D. Efthymiatos and S. Tzafestas, “A generalized approach to feedback systems in Hilbert space”, Znt. J. syst. Sci., Vol. 4, pp. 87-95, 1973. G. Zone and K. Chang, “A successive approximation method for nonlinear DP control systems”, Znt. J. Control, Vol. 15, pp. 225-272, 1972. E. Axelband, “Optimal control of linear DPS”, in Advances in Control Systems (Edited by C. Leondes), Vol. 7, pp. 257-310, Academic Press, New York, 1969. K. Chang, “On the numerical computation of a class of DPS”, IEEE Trans. autom. Control, Vol. AC-15, pp. 514-516, 1970. S. Mitter, “Successive approximation methods for the solution of optimal control problems”, Automatica, Vol. 3, p. 133, 1966. K. Chang, “Second-order computational methods for DP optimal control problems”, in “DPS Identification-Estimation and Control”, Edited by W. Ray and D. Lainiotis pp. 321-384, Marcel Dekker, New York, 1978. R. Bellman and R. Kalaba, “Quasilinearization and Boundary Value Problems”, Elsevier, New York, 1964. R. Kalaba, “Nonlinear Differential Equations and Nonlinear Mechanics”, Academic Press, New York, 1964. S. DeJulio, Ph.D. Dissertation, University of California, Los Angeles 1968. W. Ray and J. Seinfeld: “Filtering in DPS with moving boundaries”, Automatica, Vol. 11, p. 509, 1975.

Vol. 309, No. 6, June 1980 Printed in Northern Ireland

437

S. G. Tzafesras “Optimum control of DPS with time delays”, IEEE Trans. autom. Vol. AC-9, pp. 13-22, 1964. S. Tzafestas, and “DP state estimation”, in “DPS Identification-Estimation Control”, (Edited by W. Ray and D. Lainiotis), pp. 135-208, Marcel Dekker, New York, 1978. G. Dantzing, “Linear Programming and Extensions”, Princeton University Press, Princeton, N.J., 1963. G. Zoutendijk, “Method of Feasible Direction”, Elsevier, New York, 1960. J. Rosen, “The gradient projection method for NLP. Part 1, Linear constraints”, SIAM J. Control, Vol. 8, pp. 181-217, 1960. “Part II. Nonlinear constraints”, ibid, Vol. 9, pp. 514-532, 1961. A Fucco and G. McCormick, “Nonlinear Programming-Sequential Unconstrained Minimization Techniques”, Wiley, New York, 1968. G. Hadley, “Nonlinear and Dynamic Programming”, Addison-Wesley, Reading, Mass. 1964. J. Pearson, “On variable metric methods of minimization”, Comput. J., Vol. 11 pp. 171-178, 1969. “Generalized programming solution of optimal problems,” in G. Jizmagian, “Computational Methods in Optimization Problems”, (Edited by L. Zadeh et al.), Vol. 2, pp. 165-176, Academic Press, New York, 1969. P. Falb, “Infinite dimensional filtering”, Inf. Control, Vol. 11, p. 102, 1967. dimensional filtering”, SIAM J. Control, Vol. 13, pp. R. F. Curtain, “Infinite 89-104, 1975. smoothing in Hilbert spaces”, Inf. Control, Vol. 34, S. Omatu et al., “Fixed-point pp. 324-338, 1977.

(58) P. Wang,

Control

(59)

(60) (61) (62)

(63) (64) (65) (66)

(67) (68) (69)

438