EFFICIENT COMPUTATTONAL PROCEDURES
FOR OBTAINING OPTIMAL FEEDBACK
CONTROL OF DI’STRIBUTEB PARAMETER
SYSTEMS

Thesis for the Degree of Ph. D.
MICHIGAN STATE UNIVERSITY
TUMMALA RAMAMOHAN LAL

1970

L I B F; A R 3:"
Michigv 1 Star;-

‘ 9
b 1-111" * . .S: ty

nun-.-

  
 

THbs‘H

11y

* 1-”. v

This is to certify that the

thesis entitled

Efficient Computational Procedures for
Obtaining Optimal feedback Control of
Distributed Parameter Systems

presented by

Tummala)Ramamohan Lal

has been accepted towards fulfillment
of the requirements for

Ph. D. degree in Electrical Engineering & T
Systems Science

/&/«f 0. @«z/

i 5 Major professor [L

Bahama/El 3'0) M70

0.169

    
 

omomc av
all“ & SUNS'
BOOK BTHCERY IE

 
 
     

LlnnAnv amp:
swampy. Item _

._s_.—__

ABSTRACT
EFFICIENT COMPUTATIONAL PROCEDURES FOR

OBTAINING OPTIMAL FEEDBACK CONTROI.OF
DISTRIBUTED PARAMETER SYSTEMS

By

Tummala Ramamohan Lal

In this thesis techniques have been developed to
synthesize the sub-optimal feedback controls for a class of
distributed parameter systems. The original system, char-
acterized by partial differential equations is reduced to a
set of ordinary differential equations by means of a consis-
tent approximation along the spatial domain. The technique uses
no prior information about the optimal open-loop control. The
feedback parameters are obtained by solving a parameter optimiza-
tion problem.with differential constraints using a hybrid computer.

The difficulty of solving these problems on the hybrid
computer is the large set of differential equations that result
due to fine spatial discretization. The number of integrators
available on any analog computer is limited, so a decomposi-
tion principle is used to decompose a large set of differential
equations system into lower order independent subsystems. An
iterative method is used to obtain the solution. The convergence
theorems are stated and proved. With this treatment a larger
syStem (a finer Spatial discretization) can be treated than

otherwise would be feasible.

Tummala Ramamohan Lal

The second method uses the a priori information avail-
able about the optimal open-loop control to obtain the time-
varying feedback gains. The hybrid computer implementation
of this method is simple and straightforward. The timevarying
gains are obtained by sequentially solving the parameter
optimization problems on a smaller interval than given in the
problem. The number of parameters to be determined in the
parameter optimization problem is equal to the number of state
functions in the given problem. The method terminates when the

desired performance is obtained.

EFFICIENT COMPUTATIONAL PROCEDURES FOR
OBTAINING OPTIMAL FEEDBACK CONTROL OF
DISTRIBUTED PARAMETER SYSTEMS

BY

Tummala Ramamohan Lal

A THESIS
Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Electrical Engineering and Systems Science

1970

Dedicated to

Lord Sri Venkateswaraswami

ii

ACKNOWLEDGMENTS

The author wishes to express his gratitude to Dr.
John B. Kreer, for his supervision, constant encouragement
and his many helpful comments and suggestions. He also
expresses his sincere appreciation to Dr. Robert O. Barr
for his interest and constant guidance throughout the
period devoted to this thesis. The author also expresses
his sincere appreciation to Dr. J.S. Frame for giving many
helpful comments. Thanks are also due to Dr. G.L. Park and
Dr. R.C. Dubes for their encouragement and advice during the
writing of this thesis.

Finally the author owes a great debt of love and
gratitude to his wife, Droupathi, for her patience and under-
standing during his graduate study and her constant help

during the preparation of this thesis.

iii

Chapter

I

II

III

TABLE OF CONTENTS

INTRODUCTION ....................................
1.1 Introduction ...............................
1.2 Literature Survey Pertinent to the Study
of this Dissertation ............... . .......
1.3 Contribution of The Dissertation ...........
1.4 Problem Formulation ........ . ...............
1.4.1 Mathematical model ..................
1.4.2 Constraints .. .......................
1.4.3 Performance Indices .................
1.5 Controllability ............... .............
1.6 Observability ..... .........................
1.6.1 Definitions ................. . .......

ANALYTICAL SOLUTIONS FOR THE OPTIMAL ENDPOINT
CONTROL PROBLEMS ..................

NNNN
#UNH

Introduction ....... ..... .....

Necessary Conditions for Optimality ........

Example ....................

Application of Functional Analysis .........

EXISTING METHODS OF COMPUTING OPTIMUM CONTROL ...

3.1

Introduction ..... ........ ..

3.2 Existing Methods for Optimal Open- Loop
Control ................... ..
3.3.1 Sakawa' 8 method .......

3.4

3.3.2 Direct search on performance index ..

Existing Methods for Obtaining Optimal
Feedback Control ..........................
3.4.1 Seinfeld and Kumar' 3 method .........

3.4.2 Koivo and Kruh' 8 method

iv

18

18
18
24
26

29
29
29
31
34
36

37
39

Chapter Page

IV A DECOMPOSITION PRINCIPLE ....................... 42
4.1 Introduction.............. .................. 42
4.2 Definitions and Theorems ........ ...... ..... 43
4.2.1 Spectral radius of a matrix ......... 43
4.2.2 Spectral norm of a matrix ........... 43
4.2.3 Convergence of a matrix ..... ........ 45
4.2.4 Bounds for the spectral radius of a
matrix ........ . ..................... 47

4.2.5 Conditions for the existence of an
inverse of (I-M) when M is an

arbitrary matrix .................... 49

4.3 A De ecomposition Principle .................. 50

4.3.1 Algorithm ..... ..... ........ ......... 51

4. 3. 2 Convergence theorems ................ 56

4. 3. 3 Average rate of convergence ......... 63

4.4 Examples ......... . . ........ . .............. 64

4.4.1 Analytical example .................. 64

4.4 2 Computer results .................... 65

V DEVELOH’IENT OF ALGORITHM - I .............. . ..... 79
5.1 Introduction ............................... 79

5.2 Problem Formulation ........................ 80

5.3 Algorithm - I ...... ........ . ............... 80

5.4 Computer Results ....... .................... 88

5.5 Sensitivity Considerations ................. 99

5.5.1 Sensitivity coefficients ............ 99

5.5.2 Computational algorithm ............. 102

VI DEVELOPMENT OF AIGORITHM - II ........... . ....... 105
6.1 Introduction ............................... 105

6.2 Problem Formulation ........................ 105

6.3 Algorithm - II ............................. 106

6.4 Example ............ . ..... . ................. 109

6.5 Computer Results ........................... 119

VII CONCLUSIONS ................................ ..... 123
7.1 Conclusions ........ ........ ........... ..... 123

7.2 Possible Extensions .... .................... 124
REFERENCES ...................................... 126

Figure

1

.1

.5a

.5b

LIST OF FIGURES

Complete Null Controllability at Time t in
F'cl‘ ................................. ‘3 ...........

Complete Null 5-Controllability at to in P' C F .
Complete Controllability in F' ....................
Finite Tube of System Trajectories .................

Output Trajectory Corresponding to Output Trans-
formation M in the Transformed Space .............

Flow Chart of Algorithm ..... . ......................

Initial Guess, Exact Solution and Iterative Solutions
for us(t) in Example I ...........................

Convergence of the Algorithm for u2(t) in
Example I ........ ......... .........................

Initial Guesses for Figure 4.5(a) and Figure
4.5(b) .............................................

Convergence of the Algorithm for u2(t) in
Example II for Initial Guess A .....................

Convergence of the Method for u2(t) in Example
11 for Initial Guess B . ............................

Convergence of the Method for u4(t) in Example III

Initial Guess, Exact Solution, and the Iterative
Solution for u4(t) in Example III ................

Initial Guess, Exact Solution, and the Iterative
Solution for u2(t) in Example IV .................

Convergence of the Method for Three Partitions in
Example IV for a Sample Function u2(t) ............

Illustration of the Discussion in Section 5.3 ......

vi

Page

16
16
l6

17

17

53

67

68

70

71

72

74

75

77

78

84

6.2a

6.2b

6.2c

Flow Chart of the Algorithm ........................
Verification of Sakawa's Results ...................
Unconstrained Case .................................
Constrained Case ............. . .....................
Optimal Open Loop Control by Direct Search .........
Plot of 2 sinh pl/pl .............. . ...............
Plot of -2 sinh pz/p2 .. ............................
Plot of tan y 3 1/By ..............................
Plot of coth X ‘ -BX (for B > 0) ........ . ........

Hybrid Computer Set Up for the Example .............

vii

Page
89

95

110
114
114
114
116

121

CHAPTER I

INTRODUCTION

1.1 Introduction.

Recent contributions to the theory of optimal control
have been concerned primarily with Systems whose behaviour can
be described by ordinary differential equations. While many
physical systems have a Spatial energy distribution sufficiently
aggregated during the course of motion to be described by
ordinary differential equations, many others require formula-
tion by partial differential equations. As a result, the
development of optimal control theory for distributed systems
is of increasing interest from both theoretical and practical
points of view.

1.2 Literature Survey Pertinent to the Study of this
Dissertation

Research on optimal control of distributed parameter
systems was initiated by Butkovskii and Lerner (B-5), who
attempted to define certain types of control problems that
might arise. Butkovskii (B-4, B-6, B-7, B-8) subsequently con-
sidered the optimal control for a class of systems describable
by a set of non-linear integral equations, which can be derived
from linear partial differential equations. He derived a

maximum principle (in the sense of Pontryagin) embodying the

necessary conditions for optimality of such systems. However,
Butkovskii's result requires the explicit solution of system
equations, thus restricting the results to linear systems. In
addition, Butkovskii's maximum principle results in an optimal
control in the form of a solution to a non-linear integral equa-
tion involving multiple integrals. Such an integral equation

is not solvable in all cases.

The above deficiency was removed by Katz (K-2) who
formulated a general maximum principle which could be applied
to first order hyperbolic systems and parabolic systems as well
as lumped parameter systems and did not depend on the prior
representation of the system by integral equations. Egorov
(E-l, E-2) presented necessary conditions for second-order hyper-
bolic systems and parabolic systems. The most complete defini-
tion of the control problem was given by Wang (W-Z), and Wang
and'Tung GN-l), who introduced the concepts of controllability
and Observability and derived necessary conditions similar to
those of Katz (K-2) and Egorov (E-l) based on Dynamic program-
ming.

Several authors have recently considered necessary con-
ditions for specific systems. This is because the multiplicity
of possible control problems that can be conceived for dis-
tributed parameter systems is many orders of magnitude greater
than for lumped parameter Systems. Some of the reasons for

this are:

1. The boundary control has no analog in the lumped
parameter case,
2. The distributed and boundary controls are, in
general, functions of spatial variables as well
as of time,
3. There are many different ways of specifying the
admissible controls,
4. The state depends on Space as well as on time.
The different ways of Specifying "what an optimum con-
trol is" may be elaborated as follows. The fixed terminal
state problem for lumped parameter system means that the state
has to take on N Specified numerical values. For the dis-
tributed parameter case, the value of the state at the terminal
time can be Specified at every point in the Space, only at
certain points, or only in certain regions. Similarly there
are many possible cost functionals or performance criteria.
Finally there exists a very significant difference between
treating ordinary differential equations and partial differ-
ential equations. With ordinary differential equations, a
very nice uniformly applicable theory exists for treating an
nth order differential equation as an initial value problem.
With partial differential equations different classes of
equations, even of the same order, have very different char-
acteristics and must be treated differently in each case.
Brogan (B-l, B-2) extended Butkovskii's maximum prin-

ciple to systems with non-homogeneous boundary conditions.

Axelband (A-l) obtained the eigenfunction expansion for the
control of linear distributed systems. McCausland OM-l) used

a Fourier series representation of the temperature distribution
in a slab to select the input heating to bring the spatial dis-
tribution harmonics of the error distribution in a slab to zero.
Linear and non-linear programming Schemes were proposed by
Sakawa (S-l, S-2) to solve approximately an integral equation
resulting from Butkovskii's maximum principle.

The studies cited above have been based in general on
the linearity of the system or the ability to solve the system
equation analytically. Very little work is reported in the
area of non-linear distributed parameter systems. Denn (D-l)
studied a non-linear distributed control problem using varia-
tional methods and showed the linear system as a Special case.

AS in the lumped parameter systemsthe variational
calculus often yields the form of optimal control rather easily,
but the complete synthesis of optimal controls is a major prob-
lem. Seinfeld and Lapidus (S-7) applied direct search and
steepest ascent methods for solving a class of Systems described
by first order hyperbolic and parabolic equations. Wismer (W-3)
applied multilevel optimization techniques to a diffusion system
and stated that general convergence theorems are difficult to
prove. Sage and Chaudhuri (8-3) discussed the spatial and time
discretization schemes for approximately solving the problems

in distributed systemsby the known techniques of lumped systems.

Thus far we have considered only the case of obtaining
an open loop control law. Very little has been reported in the
synthesis of feedback controls for distributed parameter systems.
Seinfeld and Kumar (S-6) first obtained the Sub-optimal feed-
back controls for a class of distributed parameter systems.

Their method of obtaining the sub-optimal feedback controls

is based on the existence of the optimal open loop solution.

The feedback parameters are determined by choosing a criterion
that yields system performance that in some manner approximates
the optimal open loop behaviour. Koivo and Kruh (K-S) used the
same criterion for the design of feedback controller but deviated
from the above, by using a gradient technique in the parameter
Space to determine the optimal feedback parameters. Both used
discretization of the space variables for computational purposes.

The disadvantage in all the above cases is the complexity
of the computations because of the increased dimensionality in-
herent in these systems. The dimensionality is increased as
the Spatial discretization step becomes smaller. Either we can
discretize the necessary conditions or we can discretize the
original partial differential equations. Wang raises the
question of relative merit between these two types of dis-

cretizations.

1.3 Contribution of the Dissertation

 

In this dissertation, efficient computational procedures
for obtaining optimal feedback control of distributed parameter

systems are given. The first method uses no prior information

about the optimal open-loop control. A finite difference scheme
is used to approximate the infinite dimensional system by a
finite dimensional system. Then a multidimensional parameter
optimization technique is used to obtain the constant feedback
gains. (Chapter 5)

One of the difficulties in parameter optimization is
the large dimensionality of the approximate differential system.
The number of integrators available on any analog computer
installation is limited by the complexities involved in main-
tenance. So a decomposition principle, which divides the large
set of differential equations due to the above approximation,
into lower order independent subsystems is stated. The solution
in this case is obtained by an iterative technique. The con-
vergence theorems are proved and the theory is illustrated with
several examples. (Chapter 4)

In the second method, a priori information about the
optimal open loop control is used to obtain time varying gains
in contrast to the fixed gains. The implementation of this
method on the hybrid computer is straightforward. The time
varying gains are obtained by sequentially solving parameter
optimization problems with differential constraints. The
numberof parameters to be determined is equivalent to the number
of State functions in the distributed parameter system. The
method terminates when the desirable performance is obtained.

(Chapter 6)

1.4 Problem Formulation

 

The main prerequisites for the analytical design of

an optimum control system consists of:
i) establishing an adequate mathematical model of
the physical systems to be controlled,

ii) determining the constraints imposed by physical
limitations and design Specifications, and then
expressing them in terms of the pertinent physical
variables,

iii) selecting a realistic performance index.

1.4.1 Mathematical model

 

The dynamical behaviour of distributed parameter systems
can be described by a system of partial differential equations
or a set of non-linear integral equations, which result, in
general, from the solution of linear partial differential
equations. This thesis considers only distributed parameter
systems described by the partial differential equations of the

form:

3931(th

at = c<<2(x,t), m<x.c),x,t> (1.1)
Q(X.to) = 000:) . x e n (1.2)

where G has continuous first order derivatives with respect
to x and t and is twice continuously differentiable with

respect to the remaining arguments. In the above equation,

the following symbols are used:

Q(x,t) = Q(x1,x2,x3,...,xn,t), a p-dimensional State variable

m(x,t) = m(x1,x2,x3,...,xn,t), a q-dimensional control variable

u(xb,t) = u(t), an r-dimensional boundary control variable
independent of the Space variable x.

0 = a given finite (connected) region in Euclidean n-space;
and 0b, the boundary of 0.

8b = a linear Operator.

It can be seen from above that the state variable is not
only a function of time, but also function of spatial domain.
Thus the state of the dynamic system at any fixed time t can
be generally specified by a set of functions {Qi(x,t), i = 1,...p},
defined for all x E n. The set of all possible functions of
x defined on n, that Qi(x,t) can be any time t, will be
called the state component function space Pi, and the product
space .F = F1 X F2 X F3 x...x Pp will be called the sgggg
function Space. This definition is similar to the State space
in the case of lumped parameter systems.

The possible control variables can be placed conveniently
in two classes.

a) Distributed controls m(x,t), a q-dimensional con-

trol variable,

b) Boundary controls u(t) = (u1(t),...,ur(t)), where
m(x,t) and u(t) are piecewise continuous functions of their
arguments and are allowed to assume values from bounded convex

regions V and W reSpectively. Any control that belongs to

these convex regions is called an admissible control.

Finally, we will assume that all the problems con-
sidered in this study are "well-posed" and thus possess the
following properties.

a) The solutions to Equations l~1, 1—2 and 1-3 exist.

b) The solutions are uniquely determined.

c) The solution depends continuously on the initial

data. This says that small changes in the initial
data will cause correspondingly small changes in

the solution Q(x,t).

1.4.2 Constraints

 

In distributed parameter systems, the constraints may
be related to dynamic variables defined on certain subsets or
all of the Spatial domain a. They are essentially equality
and inequality constraints.

The class of equality constraints is of the form,

ZEx.t.Q(x.t).m(x.t)] = 0

where Z is a vector functional of its arguments defined on
certain subsets or all of a} where a. is the closure of 0.
Typical examples are:

1) Boundary conditions which represent certain inter-
actions between the dynamic system and its environment,

2) Physical quantities which are expressible as functionals

of the System dynamic variables that may be required to remain

invariant during the course of motion. A possible form of this

10

constraint is:
In F[x,Q(x,t),m(x,t)]dfl = constant

where F is a specified function of its arguments.

The class of inequality constraints is of the form:

8L S R[x,t,Q(x,t),m(x,t)] S gu

where and gu may be either functions of time, t, and/or

gt.
the Spatial variable x or constants. Typical examples are
bounded state functions of the form:

Max ‘Qi(x,t)‘ 5 Mi = constant,
x60

and bounded control variables of the form:

\m(x,t)‘ s gi(x) almost everywhere on a
or

‘u(t)‘ S M = constant.

1.4.3 Performance Indices

A generalized integral performance index for dis-
tributed parameter systems with fixed terminal time T can
be written in the form:
T
CI --= g £[Pl(t,x,Qd(x,t),Q(x,t),m(x,t))]d0dt (1.4)
For terminal control where the final time T is fixed, a per-

formance index can be defined in the form of a spatial integral:

CT =1; PO(Q(X.t).T,X)dO (1.5)

11

The problem of minimizing (maximizing) a performance
index in the form of Eq. (1.4) can be reduced to a terminal

control problem by defining

t
Q1(x st) = I£P1(T ,X ’Qd (X 31') 3Q (X,T) ,M(X,T))]d‘1’ (1'6)
0
T
and then Q1(x,T) = g PldT (1.7)

where Qd(x,t) is the desired state. Thus Eq. (1.4) is trans-

formed into the form of (1.5), that is

CI = J; Q1(x,T)dfl (1-8)

In other words CI represents the optimal transfer of the
initial Spatial distribution to a final desired distribution

in a specified time. We have seen that with T fixed, CI

could be transformed to CT’ and in case T is free, we seek
the first time when the state lies in some given e-neighborhood
of the desired state. If it is necessary and possible to choose
an admissible control such that the phase trajectory in function
Space exactly coincides with a desired distribution at T, then
the trajectory connecting the initial and desired states is

unique and hence optimal. This brings us to the concept of

controllability.
1.5 Controllability

In any control problem it is important to consider the
question ”Can any initial state of a given system be transferred

to any desired state in a finite period of time by admissible

12

control action?" We follow the definitions given by Wang and
Tung GW-l).

Let ¢(t,x,Q(x,to),to) be a solution of (1.1) with
Specified input functions and boundary conditions given in

(1.2) and (1.3). Then ¢ satisfies the following:

1) ¢Eto.x.Q(x.to).to] = Q(x,to)
ii) 3:5 = G£¢(x.t.Q(x.to),to).x,t,m(x.t)]

iii) TEtsxbsQ (Xsto):to] = u(tsxb)

The initial state of a distributed system: Q(x,to) is
said to be null controllable at time to, if there exist admissi-
ble controls (see Section 1.4.1) m(x,t) and u(t) that will
transfer Q(x,to) to the null state in a finite time T; that

is, the solution
¢[to +-T,x,Q(x,to),tO] = 0 almost everywhere in Q.

In general T depends upon both t and Q(x,to).

o,

 

The initial state is null 6-controllable at time to, if
H¢[to + T,x,Q (x,to),to]\\ s 6

where the norm is a Spatial norm and a typical Spatial norm is
MT -- x; $13. cmi

Obviously a null controllable state is also null 6-controllable.

However, the converse is not necessarily true.

13

In many systems, only the States belonging to F' (a
subset of the state function space F containing the null
state) are null controllable. This fact leads to the follow-
ing definitions:

A distributed parameter system is said to be completely

 

null controllable at time tO .13 F', F' C T, if there exist
admissible input functions which will transfer every state in
F' to the null state in finite time. (See Figure 1.1).

Similarly we can define complete null 6-controllability
.13 F'. Here, the null state must be an interior point of F'.
(See Figure 1.2).

By imposing the condition that the terminal state is
an arbitrary element in F', we have the stronger types of
controllability namely: Complete controllability in F' and
Complete 6-controllability in F'. (See Figure 1.3).

The notion of 6-controllability is useful when dealing
with approximate systems. For example the following result is
true.

If a convergent approximate system is completely con-
trollable, then the exact system is completely a-controllable.
This follows directly from the definition of the convergent
approximation, that is, for a given level of discretization,
the solutions of the approximate systems are within some
e-neighborhood of the exact solution.

Finally if the distributed system is asymptotically
stable about the null state for all initial states in F',

then the system is completely null 6-controllab1e in F'.

14

1.6 Observability

The notion of Observability of a dynamical system is
associated with processing of data obtained from observations
on the system. Thus the basic question is:

Given a mathematical model of a free dynamical system
(Control m(x,t) = O and u(t) = O) and the output transforma-
tion. ”L is it possible to determine the system State at any
time t by observing the output over a finite time interval,
(t, t+T), where T may depend on the system properties and

the output transformation Wfl

1.6.1 Definitions

 

Let eS(T) be a finite tube of system trajectories
(with no distributed or boundary control) ¢[t,x,Q(x,to),to]
defined on time interval (to, tO+T) and with Q(x,to) 6 F'(to)
the initial section of eS(T) (a subset of the state function
space) (See Figure 1.4). Let 90(T) be the tube of output
trajectories correSponding to a given continuous output trans-
formation. ”1 of all the trajectories in eS(T) (See Figure
1.5).

A distributed parameter system is said to be completely
observable in P'(to) at time to, if there exists a finite
time T and a one to one continuous mapping from 60(T) to
F'(to). If in addition to the above conditions, F'(to) = P,
then the system is said to be completely observable.

In contrast to the lumped systems, there are no precise

mathematical conditions to test the controllability and

15

Observability properties for a general class of distributed
parameter systems. The controllability is associated with
the ability of steering one system state to another in a
finite amount of time by means of certain admissible controls.
The lack of general methods to test this property justifies
the consideration of optimal end-point control problems. The
study of feedback control requires that the system be observ-
able, that is, it is possible to determine the system state
completely at any time from a finite amount of observed out-

put data.

16

F = F1 X F2 x F3 x...x EN

State Function Space

 

 

 

 

 

   

 

 

 

 

Null
State
F
Figure 1.1 Complete Null Controllability at Time to in F' c F
F
F
Figure 1.2 Complete Null 6-Control- Figure 1.3 Complete
lability at to in F' C P Controllability

in T'

17

Susan noeuowmcmua
mnu SH 2. cowumsuom
umcmue usauao ou wchSOQm
nouuoo mmHHOuommmuH usmuso m.H ousmwm mowSOuommmuH Emummm mo onus muwcwm ¢.H madman

 

 

1 mafia 1.meflu
3L
6:2 ‘I
o
95 m
mmfiuouoohmpe uamuno

o . ace .1
A uv.uz AHV o mmHHOUUSHSHH mumum .

 

./...T<\

 

CHAPTER II

ANALYTICAL SOLUTIONS FOR THE OPTIMAL ENDPOINT
CONTROL PROBLEMS

2.1 Introduction

In this chapter we obtain the necessary conditions and
analytical solutions for fixed time, free terminal state prob-
lems, where the differential constraints are either in the form
of linear partial differential equations or non-linear partial

differential equations with proper boundary conditions.

2.2 Necessary Conditions for Optimality

Before proceeding to obtain analytical solutions for
these problems, the necessary conditions are derived by using
the dynamic programming approach. This approach was used by
Wang and Tung (W-l) where integral constraints were considered.
Brogan (B-l) also applied the dynamic programming approach but
instead of integral constraints, he considered the differential
constraints in the form of linear partial differential equa-
tions. The following derivation closely follows the deriva-

tion of Brogan (B-2). Let us consider the cost functional

t
C gal PO(Q(x.tf).tf)dn +jt:l[; P1(Q(x.t).m(x.t).u(t).t)dndt (2.1)

subject to the linear partial differential equations with

side conditions

18

19

aQ/at = st(x.t> +D<x.c>m<x,t) (2.2)
Q(x,to) = QO(X) (2.3)
st<xb,t> = u(t) (2.4)

The optimal control problem is to find admissible
controls m(x,t) and/or u(t) so that the performance index
is optimized (minimized).

Let the minimum of (2.1) be denoted by

min C A 1'[(Q(x,to),tf - to) (25)
F67
where F denotes the general forcing function, and (7 is the
set of admissible controls. (Eq. 2.2 can be obtained in that
form by using the extended operator and thus converting the
non-homogeneous boundary conditions to homogeneous boundary
conditions, Brogan (B-2)).

Now (2.5) becomes

t:
_ . f
H(Q(x,to),tf - to) -;1;;{ Podn +J‘tog pldndt} (2.6)
= min { PO(Q(x.t).tf)dn
F67
tf t +3
+j‘ Pldodt +§t° {E Pldndt} (2.7)
+6 0
where
A
PC PO(Q(x.t).tf)
P A

1 - P1(Q (X,t) :Xat)

where F(x,t) denotes a general control variable.

20

Using Bellman's principle of optimality (B-9), if the
cost C is to be a minimum during the total period (to,Q(x,to))
to (tf,Q(x,tf)), then it is necessary that the cost incurred
during the shorter interval (to +e,Q(x,to + 3)) to
(tf,Q(x,tf)) be minimum also. The cost during this later
interval is equal to the sum of the first two integrals in (2.7)
so that

t: +3
1'1(Q(x,to),tf - to) = min {Ito gpldndt
F67 o

+ 1'1(Q(X:to + €)stf " t0 ' €)} (2'8)

The minimization in (2.8) is to be performed by optimizing
the first increment of the control F(x,t). After some manipula-
tions Brogan (B-2) has shown that the necessary condition for

optimality is,

2%(Q(X,t),T) = :2; {P1(Q(Xat),F(X,t),t) + (%g)t.g%}da (2.9)

In (2.9), T = tf - t has the meaning of time to go, or time
remaining to apply control to the system, and 6n/5Q is the
functional derivative (G-l). In view of the definition in (2.5)

for .n, (2.1) gives the initial condition for the differential

system in (2.9)
nEQ<x.tf),T = 0] =2]; PO(Q(x,tf>.tf)dn (2.10)

To simplify the notation, let Y(x,t) = éfl/SQ. The vector

Y(x,t) has N components, the same as Q(x,t). Now the

21

(N+l)th component can be added to Y and Q as is done in

the lumped parameter systems. Let

aQ/at
U(t,x) =
P1
(2.11)
Y(x,t)
F(x,t) =
1
Now (2.9) can be written as
all = min Pt(x,t)U(x,t)dfl . (2.12)

at Fed
When P and U are members of a Hilbert space, the inner

product notation can be used to define the Pre-Hamiltonian H

H(Q,P,F,t) 9;; PtU d0 = <P,U>a (2.13)
The Hamiltonian H0 is defined by
HO(Q,P,t) = min H(Q,P,F,t) (2.14)

F67
Thus Eq. (2.14) represents the minimum principle, which
obviously could be rewritten as maximum principle by a change
of Sign in the definition of\ P. Equation (2.12) can now be

written in the form of the Hamilton-Jacobi equation

311.. 0
3T H (Q,P,t) (2-15)

with the initial conditions given by (2.10). A pair of partial

differential equations analogous to Hamilton's cannonical equa-

tions can be found which are equivalent to the Hamilton-Jacobi

22

equation.

From the definition of Y(x,t),

3.- a__11 ......afl 2.16

at a (Q) m(at) < )
Making use of (2.15) and noting that T = tf - t, we obtain

anganelh all__

at aT at aT (2'17)
Therefore,

31:- 0

at 6H /éQ (2.18)

Directly from the definition of H, it is seen that

The above equation can be shown as follows:

For fixed x E G, Y(x,t), Q(x,t) are vector functions
of time. Now the first variation of H with respect to Y
can be obtained as follows.

For nth order cannonical equations, (2.13) gives
HQgPtUdfl= tEY':§%dn

Now,

H(Y + 6y) at

A“ +63!)t - 59d!)

H(Y)+£6Yt -§%dn

As “by“ a 0, we get for fixed x, 6H/6Y = aQ/at

23

Since x is any element of 0, we get (2.19). Kalman (K-7)

showed that if the solution to (2.15) is analytic, then

6H0/6Y (2.20)

6H/6Y
so that

6Ho/6Y (2.21)

30 /at

Thus the pair of nth order canonical equations, similar to

lumped parameter case, are

aQ/at = 6HO/6Y
(2.22)

aY/at =-an°/5Q
The initial conditions for Q are Q(x,to) = Qo(x). The

second set of conditions, for free terminal state problems, are

3P
Y(x,tf) 6Q (Q(x,tf),T 0) 5Q (2.23)

and if the terminal state is fixed then the other condition

would be,

Q(x,tf> = Qd(X) (2.24)

If the (n+l)th order canonical pair is desired, the (n+1)

components are obtained from (2.11) as

aQn+1 BYn+1

Since we are concerned with fixed time and free endpoint prob-

lems in this dissertation, we apply these necessary conditions

24

to a specific problem of this class.

2.3 Example

The problem considered here is to drive the temperature
distribution in a one dimensional conducting body from its
initial zero state to as near as possible to the desired dis-
tribution qd(x) at a fixed time t, by forcing the temperature
at one end of the body to have an optimal time history u(t).

The control u(t) is required to satisfy
|u(t)\ s 1 for all t (2.26)

The cost function is

1
2
c = £<qd(x> - q<c1,x)) dx (2.27)

The System equation is

aq/at = azq/axz (2.28)
with the initial and boundary conditions

q(x,to) = 0, q(o,t) = O, q(1,t) = u(t) (2.29)

Now applying the extended definition of the operator, Brogan
(B-2), the non-homogeneous boundary conditions are converted
to homogeneous boundary conditions to the system. Thus (2.28)
and (2.29) reduce to

aq/ac = an/ax2 + 6'(X - 1)u<c)

(2.30)
q(x,co> = o, q(o,t) = q(1,t) = o

25

where 6(x - L) is an impulse function.

The Pre-Hamiltonian as defined in (2.13) is

L
= in [5—‘1 + 6'0; - L)U(t)]d§

L 2
= £ gill [1.31365 _:§11 (6'(§ - L)U(t)')d§

2 d 6
(a—lemg - — (J

) _ u(t) (2.31)
ax d5 5Q ‘ént

ll
CDL‘—:o<~
as

The control uo(t) which minimizes (2.31) subject to the
constraint of (2.26) is, excluding the possibility of the

singular control,

uo(t) = sgn [5"an

a§ 6aka] (2.32)

Now let Y(x,t) Q (bu/Sq). Then the necessary conditions yield

aY/at = -6H°/6q (2.33)
subject to the condition at t1 that
Y(x,t) = aP/aq° = -2(qd<x) - q(x,t1>> (2.34)

Integration by parts within (2.31), so that q(x,t) appears

in undifferentiated form, facilitates finding that
2
6H°/6q = a Y/axz (2.35)
Substituting T = t - t into (2.33) we obtain,

1

aY/at = -aY/3T = -6H°/6q = -52Y/5x2

26
Therefore
2 2
aY/aT = a Y/ax (2.36)

Thus (2.34) and (2.36) form an initial value problem for the

diffusion equation. Thus the necessary conditions yield a

two point boundary value problem.which sometimes can be converted

to initial value problem. Because of the complexity of the

equations, the analytical solutions are very difficult to obtain,

and thus computational methods are used to obtain the solutions.
Similar reSults can be obtained for this class of prob-

lems by using functional analysis.

2.4 Application of Functional Analysis

 

Let the state space at a given time be denoted by H2

(Hilbert Space with L norm) and let the control variable

2

Space be H Then the solution to (2.2) with zero initial

1.
conditions and homogeneous boundary conditions represent a

mapping of elements from H into H2, and at time t can

1 1

be written as,

Q(x,tl) = Lt F (2.37)
1

The cost function to be minimized is

C = IEQdO‘) ‘ Q(xsthtEQdO‘) - Q(xstflda
O

\\Qd(X) - Q(x,tgnfiz

2

27

<Lt F,Lt F> +<Qd,Qd>H - 2<Qd’Lt F>H

1 1 H2 2 1 2
=<F,L* L F> +<Q ,Q > - 2<L*Q ,F>
:1 t1 H1 d d H2 :1 d H1
= <F,L* L F - 21* Q > +<Q ,Q > (2.38)
:1 t1 616111 d 6112

*
where L is an adjoint of L

t‘1 t1

Now the cost C will be minimized if the term depend-
ing on F is minimized, since Qd is fixed. If F0 is the
optimal control, then any other control F = F0 + 3F. will

satisfy

0 - 2 o 2
Mod - Ltlav + smug 2 HQd - Lt; HHZ
or

<F° + 95¢: L (F0 + 6F) - 2L* Qd>H 2 <F°,L* L F°-2L* Qd>H
1 t1 t1 1 t:1 t:1 "1 1
or
—' * o * 2 - * ._
6<F,2L L F - 21. Q > + e <F,L L F> 2 o (2.39)
t1 t1 t1 (1 t1 t1

But since

-*L
<F,Lt

i? — 1.: F - F2 o
1 >--<Lt ,Lt >-\\Lt1u 2

t1 1 1

Equation (2.39) requires that

23<F,L L F - L Q > 2 0 (2.40)
t1 t1 t1 d

for arbitrary e, and so

*L ° 241
LC F-LQd (.)

28

is the necessary condition for the optimality which is in the
form of an integral equation and this is the same as the solu-

tion of the two point boundary value problem expressed in the

integra 1 form.

CHAPTER III

EXISTING METHODS OF COMPUTING OPTIMAL CONTROL

3.1 Introduction
In this chapter, we will develop some of the existing

computational techniques for obtaining optimal open loop and

closed loop controls for distributed parameter systems. These

computational methods have been developed because of the dif-

ficulties encountered in solving these problems analytically

as shown in Chapter II.

3.2 ExistingMethods for Optimal Open Loop Control

As a consequence of the necessary conditions for
optimality discussed in Chapter II, two point boundary value
problems in terms of partial differential equations are obtained
This is similar to the lumped parameter case when Brogan's (B-2)

extended operator method is used to reduce the multiple boundary

value problems into initial value problems. Sage and Chaudhuri

(8-4) spatially discretized the necessary conditions and applied

the gradient and quasilinearization techniques available for

lumped parameter systems. The gradient method is based on

iteration on an assumed control trajectory to improve contin-

uous 1y the performance index. The quasilinearization technique

linearizes the state equations to generate a sequence of

29

3O

convergent approximations to the actual trajectory while re-
taining the boundary conditions (for example, the Newton-
Raphson method in function space). There are some other
methods called shooting methods where the boundary conditions
are iterated upon, while the actual State and adjoint equations
are retained. Seinfeld and Lapidus (S-7) developed two methods
called Direct Search technique and Steepest Descent method for
boundary value problems described by partial differential equa-
tions. The steepest descent method is an extension of Bryson's
method of steepest descent for optimal control problems in
lumped parameter systems. It is a gradient method based on
samll perturbations about a nominal trajectory. Sakawa (S-l)
converted the optimal control problem into a non-linear program-
ming problem. Khatri and Goodson (K-l) discussed approximate
methods of solving a class of optimal control problems using
calculus of variations. Their approximation consists of
harmonic truncation in the S-domain.

In the next section, the methods given by Sakawa (S-1)
and Seinfeld and Lapidus (S-7) are discussed in detail, since
some of the results obtained from these methods are utilized
in the computation of feedback control for example problems
in Chapter VI. AS in the case of lumped parameter systems,
each of the above methods have their advantages and dis-
advantages and no one method would serve as the best choice

for all the types of problems involved.

31

3 .3 .1 Sakawa's Method

Sakawa's method of obtaining the optimal control can

best be illustrated by an example. The process of one Sided

heating of a metal in a furnace is described by a diffusion

equation.

2
33 = 3.9. (3.1)
at 2
3x
with the boundary conditions
q(x.o) = 0
BS ... -
axIx=o a{q(o,t> v(t)} (3.2)
33‘ =
ax x=1
and the temperature v(t) is controlled by the fuel flow
u(t) and satisfies the following differential equation:
dv
r a; + v(t) = u(t) (3.3)
and
O s u(t) s 1 (3.4)
where r is the time constant of the furnace and u(t) is
normalized properly. The performance is:
1 * 2
Jiu<t>3 = j{q (x) - q(x,T)} dx (3.5)
o

where q*(x) is desired distribution and q(x,T) is the

actual distribution at time t = T. Equation (3.1) along with

the boundary conditions (3.2) can be converted to an integral

32

equation,

 

 

T
q(x,T) =3” sow: - t>u<t>dc (3.6)
o
where
2 2
8(x,t) = k cos ch-x) e~k t
cos k -- sin k
a
2 a cos (1-x)Bi (3'7)
+-2 k. E - .;r
._ 2- 2 .1. in
1-1 (k Bi)(a +' 2)cos Bi
1
where

1
k =7;- and Bi are the roots of 3 tan B = (1!

Thus the optimal control problem can be stated as follows:

Given (3.6) and the constraint 0 s u(t) s 1 on the

interval 0 s t s T, find u(t) such that the performance

index given in (3.5) is minimized.

Now the conversion of this problem into a non-linear

programming problem is given. After applying numerical integra-

tion formula to (3.5), the approximate performance index

DJ(u) is expressed as

n * 2
JEU] Z DJEu] = 2 Ci{q (Xi) - q(xi,T)} (3.8)
i=o

where C 'S are the weights assigned to the values of

integrand at the point xi. The values of x. and the weights

C1 are known for each integration formula. As an example, if

UnaSimpson's rule is used, the values of xi's and Ci are

given from standard tables as

33

xi = i/n (i = O,1,2,...,n)
Co = Cn = 1/3n
(3.9)
C1 = C3 = ..... = Cn-l = 4/3n
02 = c4 = ..... = cn_2 = 2/3n

where n is an even number.
Applying the same integration formula to (3.6), the

approximate value of q(xi,T) is given by,

n
q(x13T) : q(xiiT) g T jEO ng(xi,T - Tj)u(¢j) (3°10)
where
Tj = jT/n (j = O,1,...,n)
Putting

TC. x.,T - T. = a..
Jg( 1 J) 1]

u(Tj) = uj (3.11)

* *
q (xi) - qi

and substituting (3.10) into (3.8) yields,
n * n 2
DJ[u] ~ F[u] = .2 Ci(qi - Z a..u.) (3-12)

1=o j=o 1] J

The constraint in (3.4) is written as
O s uj s l (j = O,1,...,n) (3.13)

Consequently, the minimization problem of the functional in

Eq. (3.5) is approximately reduced to a minimization of the

function in (3.12) of n+1 variables uj's subject to the

34

constraints of (3.13). Thus the optimal control problem is

'reduced to a quadratic programming problem in this case and the

solution can be obtained with the known methods.

3.3.2 Direct Search on the Performance Index

Sakawa's algorithm essentially gives solutions for
linear problems that could be transformed to integral equations.
Seinfeld and Lapidus (8—7) extended the Direct Search method
of lumped parameter Systems. The ease of handling non-lin-
earities and control constraints as well as the success in
handling singular problems make the method attractive. Let us
discuss this method briefly, and note the advantages and dis-
advantages of this method over the others.

Assume that the interval (O,tf) is divided into L-

segments and (0,1) is divided into N-segments. We select

uk(t), k = l,2,...,q. The direct search algorithm can be out-

lined as follows:

1) Guess ui(t),u;(t),...,u:(t), the starting control

functions, where uk(t), k = l,2,...,q are the boundary controls.
2) The system equations are integrated over the given

domain with these starting control functions, to obtain the

value of the performance index Po.

3) Now allow u1(t) to vary and find u1(t) which

minimizes the performance index,

PEQ(X,t )stf] = FEQ(xatf) sqd(xstf)sulsu23°°°suq(t)]

In.other words, fix u2(t),u3(t),...,uq(t) at the assumed

....J

Y I‘-

 

 

35

values u;(t),u;(t),...,u:(t) and vary only u1(t) until
that value is found which minimizes the performance index.

We call the resulting index P1 with the control vector

1 _ 1 o o
u (t) — <u1<t>.u2<t)....,uq(t)>
At this point note that

P1 s P0 (3.14)

1 0
since at worst the control u1(t) = u1(t) is obtained.

4) Repeat the procedure for uk(t), k = 2,...,q.

5) Return to k = 1 and repeat the steps (1) through
(4) to obtain any improvement in the performance index so that

the consecutive values are within a prescribed error bound.

The direct search on the performance index offers the

following advantages for the distributed parameter systems.

1) Minimum storage capacity is required since only

the last Pj and the last control function

3 j '
u (t) = (ulcc),...,u;<t>)
has to be retained.

2) Control constraints are handled simply.

3) Knowledge of the variational formulation and the

two point'boundary value problem.is not required to use this

method.

4) Non-linear systems are handled in the same way as

linear systems .

36

The disadvantage of this method is the excessive amount
of cxxnputation time because of the large number of integrations
to be performed. As is the case in all the other methods, the
convergence to a global extremum.has not been proven in the
general case.

In contrast to the open-loop control, there are not
many studies on the computation of feedback controls for dis-

tributed parameter systems. In the next section, the existing

methods for obtaining the feedback control for terminal optimal

control problems is presented.

3.4 Existing Methods for Obtaining Optimal Feedback Control

The techniques developed in section 3.3 result in a
control which is a function of the independent variables and
the initial conditions, a so called open-loop control. From
an engineering point of view, it is desirable to have the
optimal control as a function of the state, and possibly time,
such a control is usually called a feedback control law.

Seinfeld and Kumar (S-6) first obtained the sub-optimal
feedback controls, for a class of distributed parameter systems.
Their method requires the existence of the optimal open-loop
control. The feedback parameters are chosen by minimizing a
system performance which in some manner approximates the
opthmfl.behavior. Koivo and Kruh (K-S) used the same criterion
flnrthe design of feedback controller but deviated from (S-6)

tithe actual design procedure. They used the gradient tech-

rfique in parameter Space. This method requires the transformation

of the syst

two methods
3.4

Cor
defined on
be denoted

Co

3.9
at

along with

where Q (
distribut
VECIOr’ é
ConditioI
8nd m(x

and v.

follow it
fUth 10‘

traject

37
of the system into corresponding integral equation form. These
two methods are discussed next.
3.4.1 Seinfeld and Kumar's method

Consider a parabolic or first order hyperbolic system
defined on a fixed spatial domain 0. Let the boundary of 0
be denoted by 0b.

Consider the system described by

39.:
at GEtsst(xat):Qx(xat)sQXX(xst)su(t)am(xst)] (3°15)

along with the boundary conditions

Q(x,to) Qo(x) x E Q , t E [O,tf] (3.16)

SbQ(xb,t) = u(xbst) X E Ob: t E [Oatf] (3°17)

where Q(x,t) is the p-dimensional state vector, m(x,t) the
distributed control, u(t) the boundary control an m-dimensional
vector, and (3.16) and (3.17) represent the initial and boundary
conditions respectively. In addition, we may constrain u(t)
and m(x,t) to assume values from bounded convex regions W
and V.

The open-loop optimal control problem is posed in the
following manner. Determine u(t) E‘W to minimize a scalar
functional of the state, desired state, and the control

trajectories.

c = (Undue) - Q(X.T)]T[Qd(X.T) - Q(x,T)]dn (3.18)

38

Let us assume that the open loop control is computed using any
one of the suitable techniques described in sections 3.2 and

3.3. Let us represent the open loop control (Optimal) as:

u*<t) = e<c,Q(x,o>) (3.19)

to stress the implicit dependence on the initial state.

The present problem is to determine the closed loop
control laws, denoted by uc(t), that yield system performance,
that in some manner approximates the optimal behavior. Thus
we require a criterion to compare the open loop and closed

loop system performance. One of the following criteria can

be used.
. tf 2
a) min J‘ \\e(t,Q(x,o)) - uC(t)H dt (3.20)
uc(t) 0
or
* 2
b) min & HQ <x.T) - QC(X.T)H an (3.21)
uc(t)

where Q*(x,T) is the optimal state trajectory obtained by
the application of the optimal open loop control law and
Qc(x,T) is the state resulting from the application of the
closed loop control law uc(t).

To carry out either of the above minimizations it is
necessary to assume functional forms for the feedback laws uc(t)
which include the adjustable parameters that can be determined
by the minimizations. Let us assume for convenience m = p = 1,

and

 

where kw
either of t

several was

 

feedback 1
Let
the sub-opt

mir
k
v
Expanding 1
min I
o

w

Differentj

We get

It is er
funct i0“

ditiOns

give“ Pr

39
uc(t) = kwhEg q(x,t)dn] (3.22)

where kw is the parameter to be determined by minimizing
either of the two equations (3.20) and (3.21). There are
several ways of choosing the functional relationship for the
feedback laws, and the form given in (3.22) is not unique.
Let us carry out the minimizations for k.W by using

the sub-optimality criterion (a) given in (3.20), i.e.
tf * * 2
min‘f {u (t) - k hq‘ q (x,t)dnj} dt
k 0 w 0

Expanding the terms under the integral sign, we desire
it

mini f{u’k2(t) + k:h2[z[’ q*(x,t)d0] - 2kwh[£q*(x,t)dfl]}dt (3.23)

k 0
w

Differentiating with respect to kw and equating to zero,
we get
t

J“ f{u*(t)h[l£ q*(x,t)dn]}dt

k = ° (3.24)

w tf 2 *

f h [gq (x.t)dn]dt

o

It is evident that the value of k.W obtained is an implicit
function of the initial condition, q(x,o), because of the
dependency of the open loop control laws on the initial con-

ditions of the system.

3.4.2 Koivo and Kruh's method

*
Let u (t) be the optimal open loop control for the

*
given problem, and the corresponding optimal state q (x,tf)

40

at t = tf. Then the feedback control is assumed in the form:
d
us(qsshst) = F[q8(x shot):hst] (3°25)
where
xd = col [x1,x2,...,xM]
h = col [hl,h2,...,hM]

and qs(x,t) denotes the solution of the given system when
the feedback control us(t) is applied, xd denotes the M
sensor locations, h represent the feedback constants of the
controller to be determined. The purpose of the design is to

obtain the feedback coefficients h, so as to minimize,
d *
DJ =£qu(x,x ,h,tf) - q (x,tf)‘d0 (3.26)

The method of obtaining the parameters is as follows:
Let us assume that the System is described by the
integral equation, which in the linear time invariant, constant

«coefficient case can be transformed into this form by the use

c>f Laplace transform techniques. Thus,

t
q(x,t) =j‘ g(X.t-T)U(r)d'r (3.27)
O .

Where g(x,t-"r) is the known characteristic of the system,

and let T = tf.

The first differential of the (3.26) with reSpect to
h can be written as,
M
m m
ADJEx,h; Ah] = 2 p Ah (3.28)

m=1

41

 

where
m d * aqs
p 1 sgn [qs(X.X ,h,'r> - q (xxm --; dx (3.29)
ah
d
qu T BFEQSOE shot) shat]
....a =J‘ 8(X:T"T) m d7 (3.30)
ah o ah
Requiring a constant step size
2 M m 2
“Ah” = 2 (Ah ) = constant (3.31)

m=
we have the following algorithm:

1) Compute the optimal trajectory q*(x,t) and the
correSponding performance index.

2) For each xm, approximate initially the value of

3) Compute the approximately optimum trajectory
qs(x,xd,h,t) and DJ from (3.26).

4) Compute pm, m = l,2,...,M

5) Change hm to hm + Ahm, m = l,2,...,M so as to
decrease DJ.

The Ahm used in reference (K-5) is

m
Ahm = - 2_UAEU___. m = l,2,...,M

M 3
< 2<pm>2>5

m=1

6) Repeat from (3) until the minimum of DJ is obtained.

In the next chapter, a decomposition principle which
decomposes a large differential system into smaller order in-
dependentlsubSystems is stated, and the convergence theorems

a re proved .

 

4.1 lntro<

All
equation 0:
approach".
finite diff

continuous

that the at

 

 

quantizati(
sTstems are
an the co:
approximat
SChEme has
In this ch
this thesi
composing
diSCret 12.-
systems.
complete :

thn in d

CHAPTER IV

A DECOMPOSITION PRINCIPLE

4.1 Introduction

A classical way of solving the partial differential
equation on analog computers is the so called "parallel
approach". The method replaces the space derivative by a
finite difference scheme while keeping the time derivative
continuous. Intrinsic to such an approach is the problem
that the amount of equipment required grows larger with finer
quantization of the space variable. Distributed parameter
systems are characterized by partial differential equations,
and the computational techniques require some kind of
approximation. The approximation by a finite difference
scheme has the disadvantage of demanding large computers.

In this chapter, a method which is an original contribution of
this thesis is proposed to circumvent this difficulty by de-
composing the large set of equations resulting from the space
discretization into a set of lower order independent sub-
systems. This requires an iterative technique to obtain the
complete solution. Of course the price paid for this reduc-

tion in dimensionality is increased computer time.

42

43

4.2 Definitions and Theorems

Before stating the central theorems of this chapter
with proofs, some of the necessary concepts are developed in

the following sections.
4.2.1 Spectral radius of a matrix

Let A = (aij) be an n X n complex matrix with eigen-

values xi’ 1 S i s n. Then

MA) = max “‘1‘ (4.1)
lsiSn

is the Spectral radius of the matrix A.

4.2.2 Spectral norm of a matrix A

Let A = (a,.) be an n X n complex matrix. Then

13
“AH = sup M (4.2)
“‘0 “X“

is the Spectral norm of the matrix A.

Theorem 1. If A and B are two n X n matrices, and a

is any scalar, then “AH > 0, unless A E 0,
(i) HaAH = \04 ' HAN
(ii) HA + B“ 5 “AH + NB”
um HA - 311 s uAu - nan (“'3’
(iv) “Ax“ s “A“ ”x“ for all vectors of x

where “x“ is an Euclidean norm.
Proof: The proof is given in Varge (V-l).

Corollary: For an arbitrary n X n complex matrix A,

44

HA“ 2 MA) (4.4)

Proof: If A is any eigenvalue of A, and x is any eigen-

vector associated with the eigenvalue 1. then Ax = 1x. Thus
N TX“ = HxxH = HAXH s HAMXH (45>
from.which we conclude
“A“ 2 ill for all eigenvalues of A,

which proves (4.4).

Theorem 2. Varga (V-l). If A = (aij) is an n X n complex

matrix, then

“An = [p (HA)?

*
where A is the conjugate transpose of A.

Corollary: If A is an n X n Hermitian matrix, then

\\A\\ = p (A) (4 . 6)

Moreover, if gm(x) is any real polynomial of degree m in

x, then,

Hgm(A)H = p(gm(A)) (4.7)
Proof: If A is Hermitian, then A = A*, and thus

“An?" = p<A*A> = p<A2> = fun

HAN = p (A)

Now since

(A) is a
8m
4 .2
Let
convergent

converges t
Theorem 3.
—-—~

convergent

 

 

 

 

Proof: For

| 'fi

11X [1 matr

normal for“

SAE

Rhere each

45

Now since gm(x) is a real polynomial in the variable x,

gm(A) is also Hermitian, and (4.7) is proven, because

\Tgmmuz = p<sm<A)-g;<A)> = p2(gm(A))-
Q.E.D.

4,2,3 gonvergence of a matrix A

Let A be an n X n complex matrix. Then A is
convergent (to zero) if the sequence of matrices A, A2, A3,...
converges to the null matrix, and is divergent otherwise.
Theorem 3. If A is an n X n complex matrix, then A is
convergent if and only if p(A) < 1.

.nggfz For the given matrix A, there exists a non-singular

n X n matrix S, which reduces the matrix A to its Jordan

normal form, i.e.

 

 

'.'J1
J O
2
-1 ~ ‘
SAS =A:—: O ‘\ (4.8)
\
J
L 5
where each of the nL X nL submatrices JL has the form
F 1
1. o
1
"t
J = l , l s s 4.9
t )‘t t r ( )

 

 

46

Since each submatrix JL is upper triangular, so is A.
Thus the set {XL}:'1 includes all the distinct eigenvalues
of the matrices A and A, which are similar matrices from

(4.8). By direct computation with (4.8), we get

 

 

 

 

 

 

pm
J1 '1
J: O
(.711)m = “ , m z 1 (4.10)
o
L H
The entries of the powers of the matrix JL are determined as
follows:
1 0 0“ 1 0 07
F1, F1,
0 f‘t 10 .0 0 i 10 0
2
J =00 10
L )1 O
' 1 O 1
O M.
— xL-J L XL“
2 '1
i- 2 10 0
AL KL

 

 

and in general if we define

“1 ___ (m) .
JL (dij (L)) 1 S 1, j 2 nc

then

47

df?)cc) = 0 j < i
“inf-j“ for i sj s min(n{’,m+i) (4.11)
0 m+i < j s n
t
where
(m) = -——E!——-—
k k!(m-k)!

Now if A is convergent, then by definition in 4.2.3 Am a 0
as m a a. But (A)m = SAmS-l. So it follows that Km‘q O as
m S.¢, Consequently each Tém) a O as m-4 a so that the
diagonal entries KL of JL must satisfy ‘XL‘ < 1 for all

1 s L s r. Clearly

p(A) = p(K> = max Ti,
ISLSr

T < 1

which proves the first part. On the other hand if

p(A) = p(A) < 1, then ‘XL‘ < l for all 1 s L s r. Then by
making direct use of (4.11) and the fact that 11%| < 1, it

follows that

. m _ . _ .
11m dij(L) — O for all 1 2 1, j 2 “L

Thus each J is convergent, and A is convergent. Finally,

L
Am = 3'1 Ems
This proves that the matrix A is convergent.

4.2.4 Bounds for the spectral radius of a matrix

It is generally difficult to determine precisely the

spectral radius of a given matrix. Nevertheless, upper bounds

48

can be easily found from the following theorem:
Theorem 4. Let A = (aij) be an arbitrary n X n complex
matrix, and let

n

A' E jgl laijl

j¥i

l s i s n

Then all the eigenvalues A of A lie in the union of the

disks,
|z - 811‘ 5 A1 1 s i s n (4.12)

Proof: Let A be any eigenvalue of the matrix A, and let x
be an eigenvector of A corresponding to 1- We normalize the
vector x so that its largest component in modulus is unity.

By definition,

n
(A - aii)xi = Z aijxj 1 S 1 S n
1‘1
j¥i
In particular, if ‘xr| = 1, then

n n
lA ' 8rr‘ 5 jEIlarj‘°lxj‘ s jgllarj‘ = Ar
j¢r j¥r

Thus, the eigenvalue 1 lies in the disk ‘2 - arrl s Ar'
But since 1 was an arbitrary eigenvalue of A, it follows
that all the eigenvalues of the matrix A lie in the union
of disks ‘2 - aii‘ 3 A1, 1 s i s n completing the proof.
Corollary. If A = (aij) is an arbitrary n X n complex

matrix and

49

n

u 5 max 2 ‘a. l (4.13)
lsisn j=1 lj

then p(A) S A.

Thus the maximum of the row Sums of the moduli of the entries

of the matrix A gives a Simple upper bound. Since A and

t .
A have the same eigenvalues,

n
”'5 max 2‘8

| (4.14)
lstn i=1 lj

then p(A) s u'.
4.2.5 Conditions for the existence of an inverse of

(I-M) when M is an arbitrapy matrix

Theorem 5. If M is an arbitrary complex matrix with p(M) < 1

then I-M is nonsingular, and
-l 2
(I-M) =I+M+M +... (4.15)

where the series on the right converges. Conversely, if the
series on the right converges, then p(M) < l.

‘nggf: First assume that p(M) < 1. If u is an eigenvalue
of M, then 1 - u is an associated eigenvalue of I-M, and,

as p(M) < 1, I-M is nonsingular. From the identity,

r+l

2 1-
1+u+u +...+nr=—-1-ETE

and substituting u = M, we get

2 +1
I - (I-M)(I +M+M +...+Mr) =11r

we have, upon premultiplying by (I-M)-1, that

50

1

(I-M)' - (I +M + M2 +...+ Mr) = (1-M)"1‘Mr+1

Thus,

-1 2

H(I-M) - (1 + M + M +...+ Mr)“ 5 \\(1-M)’1H-\\Mr+1\\

for all r 2 0. As M is convergent, it follows that

(Mr-*1“

a 0 as r a n. Thus the series in (4.15) converges
and is equal to (I-M)-1. Conversely, if the series converges,

let u be the eigenvalue of M, correSponding to an eigenvector

x, then

(I +M+M2 +...)x = (1 +p, +p,2 +...)x

Thus the convergence of the matrix series implies the con-
vergence of the series 1 +-u +~u2 +... for any eigenvalue

u of M. However, as is well known, for this series of complex
numbers to converge, it is necessary that ‘u‘ < l for all

eigenvalues of M, and thus p(M) < 1, completing the proof.

4.3 A Decomposition Principle

 

Many of the physical systems described by partial dif-
ferential equations involve at most three dimensions in the space
domain. To solve such problems in three space, a spatial dis-
cretization is used which yields a set of ordinary differential
equations. The number of ordinary differential equations in
this set increases rapidly as finer and finer Spatial discretiza-

tions are used.

51

One method of solving this set of differential equa-
tions is by the use of an analog computer. The difficulty is
that the number of integrators on any analog computer is limited
and the cost of the equipment increases markedly when additional
sophisticated integrators are added to the available facility.
The second method of solving this set of differential equations
is by the use of digital computers. The disadvantage in this
case is the large amount of memory required and the increase
of computation time with the number of equations. Hence a
hybrid computer solution obtains the advantages of both the
analog computer where the differential equations are solved
in parallel and the digital computer is used for the logical
and control functions.

Thus the necessity of an algorithm for obtaining at
least an approximate solution of the given system of equations
with smaller number of integrators and small amount of core
memory is established. In the next sections, an algorithm is
stated and proved, which increases the capabilities of the
analog and digital computers and thus obtains an approximate
solution to a larger set of equations with fewer integrators

and less core memory.

4.3.1 Algorithm

Many of the systems characterized by the partial dif-
ferential equations yield tridiagonal matrices when discretized.
If the discretized matrices are not in the tridiagonal form,

they can be reduced to this form without computing the eigenvalues

52

by using several existing methods, such as in Ellsworth (E-4).

Now, given the equation

x = A x , X(to) = x0 (4.16)

where A is a tridiagonal n X n real matrix, X is l X n
vector, and X0 is the initial condition vector. We define

a partition of X as follows:
= l
X (X1, X2, X3,...,XN)

where X has dimension nj. The decomposed problem thus also

i

partitions the matrix A into N2 blocks Aij such that,

+ A X +-A

= A21 x1 22 2 23

(4.17)

’hq-l = AN-l,N-2 xN-Z + AN-l,N-l hi4 1’ AN-1,N XN

in = AN,N-l 131-1 + AN,N xN

Then the iterative method of obtaining the solution of the
above equations is given in the flow chart in Figure 4.1.
This method is like the Gauss-Seidel method for the
solution of linear equations, since the values obtained for
the other variables are utilized immediately. The initial
vector P, which is the same for all the variables, is used
to reduce the storage in the digital computer. The Jacobi

method can also be used, but it increases the storage

 

 

Guess Xé?i(t) = P and compute Xé1)(t),

by integrating the differential system

 

 

correSponding to XN and store as Xé1)(t).

 

 

 

Guess x39; = Xégg =...= Xio) = P and

compute XN}i(t) by using Xéfi+1(t) and

the guessed value of xN-i-l(t) and store.
0 1 0
Compute 31:: =Hx1:_)i - fill“

l

I = I + 1‘]

- T
@‘—No ——I Is I=N 7]

 

 

 

 

Yes

A

J = J + 1 1

6

Figure 4.1 Flow Chart of the Algorithm

 

 

e

Compute and store XéJ) by using X;_1

<J)_(J-1)
XN

7.1T
' (J )

Compute and store XN_) by using XN_ 1+1

 

 

 

 

 

 

 

and (Jlli and compute
J J (J 1
suns).- x, W

_l__

Is I = N 7

l

<:)"‘NO Is sin desired value

 

 

 

 

 

 

 

 

“J

 

Yes

 

Print the values xl’°°°’le

Stop |

Figure 4.1 Flow Chart of the Algorithm (contd.)

 

 

 

55

required because the values of the variables corresponding to
two consecutive iterations are to be stored.

Thus loosely Speaking, the decomposition principle
starts with the last partition and sequentially proceeds to
the first and cycles from last to first to last etc., until
the error is less than or equal to the stipulated value.
tridiagonal matrix.

Assume that A is a 24 X 24

Let Z be divided into two partitions and let each partition

12 X 12 matrix, i.e.,

 

 

 

 

contain a

12 x 12 7

matrix Q

a
A = b
0 1. x ..

"’ matrix «H

where a and b are non-zero elements.

Then a hybrid computer that has twelve integrators,

24 X 24 matrix with the storage

could be used to solve this
If the matrix contains different

of only two functions.

numerical values, these can be adjusted by means of the digital

Now, if we suppose that each function is sampled

computer.
at 50 points in the interval of interest, then two functions

require only 200 words of core memory. In the next section,

the convergence theorem for this decomposition principle is

given.

 

56

4.3.2 Convergence theorems

First a theorem for N = 2 is proved and then it is
generalized.

Theorem 6. Given the equation
X = A X , X(0) = Xo on t E (0,tf)

Let A be partitioned into two partitions, i.e., N = 2 in

(4.17) and be rewritten as

A +
1 1X1 31x2

A
2 2x2+32x1

x.
II

(4.18)

x.
ll

while requiring that A1 and A2 are stable. Then a proper
choice of At such that p(AlAt) < 1, p(AzAt) < l,
l) guarantees the convergence of the algorithm in
section 4.3.1
2) p(AlAt),p(A2At) controls the rate of convergence.
.Egggf: Let (O,tf) be partitioned into M subintervals and
ti 6 (O,tf). Then (4.18) can be approximated at ti by using

the midpoint approximation

 

X (t ) - X (t. ) B
l i 1 1-1 = A_ ‘_l
or
A B A

-.- _1. __l _1

X1“? xl(ti-l) + [2 X1“? + 2 x2“? + 2 x1(t1-1)
B 2
+ ElX2(ti_1)]At + 9(At ) (4.19)

57

Similarly X2(ti) can be written as

A2 82 A2
x2021) = X2(ti_1) + [5—2 x2 (t) +— X1.(t) +--X2(ti_1)

B2

+ ‘2'" x1(ti-1>]At + 9(At2) (4.20)

Equations (4.19) and (4.20) could be written in the matrix forms

 

 

 

 

 

 

A B " ' T"
X1(ti) 52‘- 2—1- x1(tin I +2.19% B1Ag X (ti_1)-1
-= At + A_t_
X2“? B2 f; (3‘2“? L32 2 “:22? x2(ti-l)
— 2 .4 — L .1
2
+ 9(At )
or
X1“? .. A1 B1 A2 x1“? + "’11 “’12 x1("1-1)
~ 2
x2“? B2 A2 X2“? m21 22 x2(t1-1)
2
+ 8(At ) (4.21)
Substituting
xl(ti) = (1 A1 1[B1LX2(t1) + mllx 1021- 1) + m12x2(t1-1)]
+ 8(At2)
_ t -1 t
X2“? " “*2 2') [32 27x1“? +m21x1(t1-1) +m22x2(t1-1)3
+ 9(At2)
we get
' L' L
X1(ti) gar. A+B1(I-A22 ) 1132 2 0 X1(ti)
X2“? 2 0 Azfiza'Al £2EY131 2'2 X291)
X10314) 2
+1: + 6(At )
x2<ti_1)

 

58

 

 

where
" A5." “'1 AE AP. 1 AE .
(I A2 ) In21"“‘11 (1 A2 ) tn22"“‘12
N 9 (4.22)
- -AE. '1 A15. - At '1 fi
“1 A1 2 ) 2 m11“““21 (1 A1 2 ”7) 2 ““12“”22.

Now let At2 be sufficiently Small that all the terms multiplied
by At2 can be neglected. In the above the inverses (I-AlAt)-
and (I-A2At)-1 exists if and only if the matrices A1 and

A are such that p(AlAt) and p(AZAt) < l by theorem 5, in

2

section 4.2.5. The series then are

95 + 6(At2)

- £4.
(1 A12) I+A12

At. 2
22 +6(At)

(I -A22£)-1 I+A

and the approximation

- .AE '1 : .AE
(I A12) I-i-A12

- A£.'1 2 .AE
(1 A22) 1+A22

The neglected terms are of the order of At2/4. The approxima-

tions of the various terms in the matrix in (4.22) are as follows:

 

 

 

59

AEQ'IB .AE:,§ A

At. -
A12 +3“I A22 2

At
12

2

2

At _ At -1 At
A22 'H32(I A12) B1 2 “”35“”2

_ A2492 ..
(I A22) 2 m21+m11~m11 I+A1At/2

(4.23)
_ M -1 Lt =
(I A1 2 ) 2 m12 + m.22 m m22 I +A2 At/2

_ AS '1 A}. AP. = A3. AL

(1 A22) 2 m22+m12”m12+12 B12 +12
- A3 '1 A9. A3. = AE AS.

(1 A12) 2 m111”“21’*"“21"’12 B22 +I2

Hence (4.22) can be approximated as

At At Bl At
X1“? = A1 2 0 x1“? I+A1 2 2" M + I 2
AL AE AE AP.

x2051) 0 A2 2 X2“? B2 2 + I2 I + A2 2

Thus we can write for known values of X1(ti_1), X2(ti_1)

and At,
X(ti) = P X(ti) + C (4.24)

where G is a constant matrix.
Let us assume X* is the unique solution of the above equation.

Then writing

:11

* [ll
x-X=e

'where cm is the error at the mth iteration,

 

 

 

 

- *
em = P xm 1 - P x
- *
= F(xm 1 - x )
g P em-l
= Pmeo (4.25)
where
t 1
A2.
A1 2 0
p A
At.
0 A2 2 ,
b .J

 

 

Thus the error vectors tend to zero, if and only if the

spectral radius p(P) is less than one. This can be done

by proper choice of At. This also requires that the decomposed
systems be stable since any perturbations of the solution Should
not make the systems unstable.

The approximations made in (4.23) are accurate because B and

 

 

 

 

1
B2 are of the form
n- .1 F 1
0 0 0 0 O 0 0 O 0 *
0 and 0 respectively,
0 O O
* O
L .. L 4
where ‘* represents the non-zero element.
Q.E.D.

Now let there be N = 3 diagonal blocks Aii’ such that

then by similar

61

= A X +-A X

11 l 12 2

A X +-A X

21 1 22 2 + A X

23 3

A

X + A

32 2 X

33 3

computations mentioned in the above,

 

 

 

 

 

 

 

 

 

 

1 1
FX1(ti) F1111 A12 0 ‘1 rx1(ti)
=AE—
x2(t1) 2 A21 A22 A23 X2(t1) +
X (t.) O A A X (t )
b 3 1.1 b 32 33.1 b 3 1..
F ‘AE t
I + A11 2 A12 2 0 x1(ti-111
A£ AE AE
A21 2 I +"‘22 2 A23 2 x2(t -1)
AB AE
- 0 A32 2 I +"A33 2 x3(ti-1)~
2
o 0 ‘ t
and by suitable approx1mat1ons that neglect Az— terms,
= A-t- ' =
xj(ti) Ajj 2 xj(ti) + Gj’ j 1,2,3 (4.26)
where G is a constant matrix for a given Xj(ti_l), j = 1,2,3
and At.

Similarly for N partitions,

 

 

 

 

 

 

 

 

1 P 1 P- 1 r ‘
[x1(t1) A11 A12 0 0' x1(t1) x1(t'-1)
x2(t1) A21 A22 A23 0 0 X2(t1) x2(t1-1)
X3(t1) 0 A32 A33 A34 0 . . o X3(ti) X3(ti_1)
=AE

2 + M
(t) : \\ (c) (t )
51%! 1‘ L JZU‘J 1.; ‘38“ 1-1.1

62

where M is a known matrix of the form

 

 

fi
.AE AE . . .
F1 + All 2 A12 0 o
M Q A21 %5 1 + A22 gi A23 135 o o (4.27)
but d

Since At 15 fixed and X1(ti_1),X2(ti_1)... are known, the

matrix M is completely known. Hence we get,
At
X.(t. =A.. X. t. +G. 4.28
J 1) JJ 2 3(1) 3 ( )

where j = l,2,...,N

. . A; .
and thus the convergence 18 assured 1f p(Ajj 2 ) is less than
unity. Thus the following theorem is established.

Theorem 7. Given
X=AX ,X(0)=xo on t6 (0.tf)

where A is a time invariant tridiagonal matrix partitioned

such that

x1 = A11 x1 + A12 x2

X = A 1 X +'A X +'A

2 21 222 23x

3

(4.29)

° = A + A +
XN-l N-l,N-2 )31-2 N-1,N-1 xN-l Ala-LN xN
{(N = A x + A x

N,N-l N-l N,N N

A11,A22,...,ANN are

stable and if the Spectral radius of the matrices A11 At/Z

the algorithm.converges if the matrices

 

 

 

63

t t
Azz‘g- ,...,.ANN'%- is less than unity for a suitable choice

of At.

4.3.3 Average of rate of converggnce

 

From Equation (4.25) we have,

Using the matrix norm defined in section 4.2.2 and the vector

norm, we have [.

 

Hem“ s HPmH ' ne°n (4.30)

assuming “so“ is not a norm of a null vector then “Pm“

gives some idea about the rate of convergence. Thus if

“Pm“ < 1, then

m
R(Pm) = 4m [upmul/m = 153M (4.31)

m

is the average rate of convergence for m iterations. Eq.

(4.31) can be written as

m
l/m = e-R(P )

m
HP \\
Therefore,

ML 5 e-w‘“)
\\e°\\

Let a = “emu/“6°“ and [R(pm)]‘1 = um, then

N 1
omS; Mam

64

so that lfln is the measure of the number of iterations required

to reduce the norm.of the initial error vector by a factor e.
In the next section, some examples are presented illus-

trating the above theory. The examples include an analytical

example as well as computer examples.

4.4 Examples

 

4.4.1 Analytical example

Consider the 2 X 2 matrix equation

 

u1 -2 1 L11 u1(0) 0

o - ’ g
u2 1 2 u2 u2(0) O
for which the unique solution is u1(t) = u2(t) = O for all

t; and let the number of partitions be 2, i.e.,

-2 l
A =
l -2
' = _ = . 3
L11 2 u1 +u2 u1(0) 0 (4 3 )
' = -2 = .
u2 u2 +u1 u2(0) O (4 34)
Equations (4.33) and (4.34) can be rewritten as
'2t t 2T
u1(t) = e I e u2(T)dT (4.35)
0
-2t t 2t
u2(t> = e (f) e u1(T)dT (4.36)

assume u‘i(t) = constant = 1, then (4.35) and (4.36) gives

 

 

ugh.)

[1(1) (t)

Note that ui1)(t) is

computation, use,

u(l)

1ct>=

Then

u(2)
2

.9)

II

(t) =

(t) =

 

m
I
n:
n
oe—an (Dc—'1” o'L—wn
m
[Q
4
5
NA
H
V
A
4
V
a.
.1.

H (D
N N
4
A
H
I
(D
I
N
.1
v

I
N

 

always positive; so now to simplify

max ui1)(t) = 1/4
0< t<oo

 

 

t -2¢
-2t 2T.l =‘l l-e
e [e 4d'T 4(2)
0
-2c t 2T (2)
e i e u2 (T)dT
-2t -2t
1 l e t e
xix-4 ' 21"”

max ufz) (t) = 1L6-

t2 0

Thus it is obvious that,

uil+1>

u§i+l)(t)

and as i a a, the solution u1(t) = u2(t) = 0 is obtained.

(t) < uii)(t)

< uéi)(t)

4.4.2 Computer results

Example 1. Consider a matrix differential equation

X = A X ,

X(0) = 4.95

 

 

Then part

with u_(

l
The initi
same figu
tion afte
Convergem
gives the

defined a4

The OSCiIl

 

0f the mat

are giVEn

 

66

where,

01-21
L001-2J

 

 

Then partition A into two partitions as follows.

u -2 l u 0

.1 = 1 + u3 (4.37)
u2 l -2 112 1

a -2 1 u 1

3 .. 3 + u2 (4.38)
u4 l -2 u4 0

with ui(0) = 4.95, i = 1,...,4.

The initial guess of u2(t) is given in Figure 4.2. In the
same figure the exact solution u2(t) and the iterative solu-
.tion after iteration l and iteration 2 are given. Note the
convergence of the solution to the exact solution. Figure 4.3

gives the convergence of the algorithm. The error norm is

defined as,
“e(i)\\ = u2(i+1)(t) - u2(i)(t)‘ (4.39)

The oscillation of the error from .005 to .008 and back to
.005 is due to the limitation of the analog accuracy.

By the corollary to Theorem 4, the spectral radius
of the matrix given in (4.37) and the matrix given in (4.38)

are given as follows.

 

 

 

Iteration 2

r______ __._

  
 

 

 

 

 

 

I
l

1.01) '
: Initial Guess
I

0 ‘0..4 036 0.8 1.0—
time —o
1 second 44,

 

Figure 4.2 Initial Guess, Exact Solution and Iterative
Solutions for u2(t) in Example I.

“EU 68

 

 

100. Figure 4.3 Convergence of Algorithm for
u 2(t) in Example I.
i+1 1
\\e\\=n§ )< -u'§><t>u
\
104»
1.4»
0.14.
0.01..
Error Limited
0.008 by the Accuracy
of the Analog
.0050 Computer
0.005
0.001 ‘r = 3 3 ; = e
O 1 2 3 4 5 6 7

Iteration Number

69

0019 s3
and At = 1/50
p (A22) _<. 3

3 3
p(A11 At) "§0 ’ p(A22 At) 5 50

and hence the algorithm converges.

Example II. Consider a 4 x 4 matrix and divide it into two

 

partitions, i.e., N = 2. The matrix is given as follows.

 

 

F" '1
-2 l 0
1 -2 l 0 .
A = 0 1 _3 1 where X = A X , X(0) = 4.95
0 0 1 -3
L J
and thus,
111 -2 1 111 0 u1(0) 4.95
. = + “3 a
U2 1 -2 L12 1 u2(0) 4.95
u3 -3 l u3 l u3(0) 4.95
. = + u2 =
114 ’ 1 -3 L14 0 114(0) 4.95

Figure 4.4 shows the initial guesses given for example II.

Figure 4.5 shows the convergence of u with the guesses

2
given in Figure 4.4. The Spectral radii are estimated by

using the corollary of theorem 4 and At is chosen.

p(All) S 3
p(AZZ) s 4

_ 1 .2.
At - 50 ’ p(All At) S 50 3

4
9(A22 At) 5 50 ‘

so;

70

 

L120)

4.0 ..

3.04

2.0.7

1.0 d»

 

P
D

 

 

 

 

 

Figure 4.4 Initial Guesses for Figures 4.5a and 4.5b.

1 T0

time d

 

 

71

100. -;
10. -b
1. ..
0.1 q»
0.01 0

0-001 . r . : : : ‘r r

L 2 3 4 5 6 7

Iteration Number
Figure 4.5a Convergence of the Algorithm for u2(t) in Example II.

for the Initial Estimate A.

72

He“

Figure 4.5b Convergence of the Method for

100-‘ u2(t) in Example II with the

V

Initial Guess B.

10.4

1.0 l

0.014

 

0.001

 

 

NA?
UH}

B

Iteration Number

“I?
4.x:-

73

Example III. Given an 8 x 8 matrix of the following type.

 

The exact eigenvalues are known to be -4 sin2 n/18

,J

F-z 1 0 0 0 0 0 0

1 -2 1 0 0 0 0 0

0 1 -2 1 0 0 0 0

A = 0 0 1 -2 1 0 0 0
0 0 0 1 -2 1 0 0

0 0 0 0 1 -2 1 0

0 0 0 0 0 1 -2 1

0 0 0 0 0 0 1 -2
L a

 

 

and the matrix equation
X=AX ; x(0) =4.9S.

The number of partitions are two, i.e. N = 2. Thus
we have two independent 4 X 4 matrices. Figure 4.6 shows
the convergence of u4(t) and Figure 4.7 gives the initial
guess and the value of u4(t) after 5 iterations and the exact
u4(t). The spectral radii of these matrices are found in the
same way as in the previous examples. The value of At is

chosen as 1/50, and the interval of interest is (0,1).

Example IV. Finally a 6 X 6 matrix is considered and the

 

number of partitions are 3, i.e. N = 3. The partioned matrices

and the given matrix are given below.

F-z 1 0 0 0 0‘)

1-2 1 0 0 0

0 1 -2 0
A: 0 0 -2 1

0 0 0 1 -2 1

L0 0 O 0 1 -2..

 

 

ueu 7"

 

 

150.1
Figure 4.6 Convergence of the Method for
u4(t) in Example III with the
Given Initial Guess in Figure 4.7.
100.01%
10.0..
1.01>
0.14,
0.01.,
0.001 : ; : : ; r 1‘
0 1 2 3 4 5 6
Iteration Number

75

 

 

 

5.0 .1
”4“) W.
‘~¢L“-° After 5 Iterations
Exact
\o\
.1
4.0 4+
Initial Guess
3.0 _ . L
0 0 2 0.4 0' 6 0 8 1'0
time —.

 

 

 

Figure 4.7 Initial Guess, Exact Solution, and the Iterative
Solution for u4(t) in Example III.

"Ii-r!

 

and

 

Figure 4
initial g
are diff:
gives th1
Wesmc:

is choser

b—J

because ;

Which ena

 

76

and

L11 -2 1 u1 0 u u1(0) 4.95
02 = 1 -2 02 + 1 3 u2(0) g 4.95
03 -2 1 03 uz u3(0) 4.95
. = + , a

U4 1 -2 u4 u5 _ u4(0) 4.95
us -2 1 us 1 05(0) 4.95
. = + U4 , =

L16 1 -2 u6 0 u6(0) 4.95

Figure 4.8 shows a sample function u2(t) along with given
initial guess. The exact solution and the iterated solutions
are difficult to distinguish after 5 iterations. Figure 4.9
gives the convergence of u2(t) with the norm given in (4.39).
The spectral radii are same as given in example I, and At
is chosen to be 1/50. The interval of interest is (0,1).
Thus these examples show the usefulness of this method
because of the fast convergence exhibited by these examples.
The spectral radius can be obtained by using the theorems given

which enables one to choose At.

 

..w

3.01

 

 

F 0
lgur¢

 

77

 
  

  

 

u2(c>
K\\
0
Exact and Iterated
4.0‘_ Solution
3.0"
2.0..
1.01L
: : £ 2 c
0 2 0 4 0.6 0 8 1.0
time —0

Figure 4.8 Initial Guess, Exact Solution, and the Iterative
Solution for u2(t) in Example IV.

78

M

\\e\\ = 1x“) - x91

100.0 1» 2

10.0

 

 

 

 

1.0
0.1
0.01
an
0.001 4. é 4. i
0 1 2 3 4

Iteration Number

Figure 4.9 Convergence of the Method for Three Partitions
in Example IV for a Sample Function u2(t).

CHAPTER V

DEVELOPMENT OF AlGORITHM - I

5.1 Introduction

In this chapter, a description of Algorithm-I for

obtaining an optimal feedback control for a class of distrib-

uted parameter systems is given. This method is different

from the existing methods in the following way:

i) No a priori information of the existence of the

ii)

iii)

optimal open-loop control is necessary.

The disadvantage of computing optimal open-loop
control, whenever there is a change in the initial
distribution is removed.

The computational method for obtaining feedback
parameters is more efficient than the existing
methods in the sense that more accurate solutions
could be obtained. This is because of the extended
capabilities of obtaining solutions for larger
differential system using the decomposition

algorithm described in Chapter IV.

The problem is formulated in the following sections and

the algorithm is developed. The algorithm is illustrated by

an example with different constraints on the control.

79

80

5.2 Problem Formulation

Consider a distributed system characterized by a

vector partial differential equation
59-= c[c x Q(t x) Q (x C) Q (x t) u(t) m(x c)] (5 1)
at 3 3 3 3 X 3 3 xx 3 3 3 3 '

along with the boundary condition

Q(x,to) QO(X) x e n , t e [O,tf] (5 2)

SbQ(xb,t) = u(xb,t) xb 6 0b, t E [O,tf] (5.3)

where the symbols are explained as follows.

0 : a given finite (connected) region in Euclidean n-space
and ab, the boundary of a.

G : spatially varying differential operator on Q which
may include parameters which are linear functions of
Q,m,x or t.

Qo(x) : initial state vector, i.e., at t = 0, Q(x,t) = Qo(x)

u(t) : boundary control.

5.3 Algorithm-I

The algorithm considered here involves forming a semi-
discrete approximation of (5.1) through (5.3) by placing a grid
on the spatial domain.

The Spatial variables are discretized by defining a

vector,

xi = (11(Ax1), i2(3x2).--., ij(Ax )..-.. in(4xn))'

J

81

which in effect places the grid on the region 0. The prime
denotes the transpose. Here the elements of i = [i1,i2,i3,...,i ]

are integers defined by i, = O,1,...,Nj where

(Xj)max - (Xj)min
N, = (5.4)
J ij

Assuming that the operator G is at most second order in X,

 

it can be approximated as follows:
G(Q(xi,t).m(xi,t).xi.t) as

GiEQi(t).Q 1(t).Q (t). °"Qiiln(t)’mi(t)’u(t)’t] (5.5)

1:1 1112

where Ik = {i/i = 0 except for the kth element which equals to l}

i ranges over all the interior mesh points, and the functions
Gi are assumed to be real valued anc class C2. As an example,
consider a rectangular mesh in E2 and using the above notation,
i = (i1, i2) and consider the mesh point (1,1). Then
I1 = (1,0) and I2 = (0,1). Therefore the points that will be
considered are (1,1), (2,1), (0,1), (1,2), (1,0).

Thus following the above notation, the discretized

vector partial differential equation in (5.1) can be rewritten
along with (5.2) and (5.3) as follows.

39, _
at - GiEQi(t),Qiiil(t),Qii12(t).....Qiijnit).mi(t).t,u(t)]

along with the discretized versions of the boundary conditions,

82

Qi(t = 0) -- 010 xi 6 0 , c 2 o (5.7)
SbiQi = ui(t) xi 6 ab, t 2 0 (5.8)

Now, before stating the optimal feedback control problem, a
brief discussion to motivate the material in this chapter is
given.

In general, there are several methods of obtaining
open-loop control for distributed parameter systems. Some of
the methods are discussed in Chapter III. But in practice,
it is desirable to have a closed-loop control such as optimal
control as a function of state and possibly of time. Thus in
the case of distributed systems, feedback methods Similar to
lumped parameter systems can be discussed. Though it is dif-
ficult to obtain analytically the controller for a large class
of problems, assuming that it is possible, the implementation
of this control law is difficult. This difficulty arises
because of the infinite dimensional character of the state
vector which is a function of Spatial domain as well as time.
So some kind of approximation is necessary so that it is possible
to reconstruct the State function by a finite number of measure-
ments along the Spatial domain, while keeping the time con-
tinuous. Then a polynomial fit can be used to get the complete
state function. The coefficients of the approximating poly-
nomial will vary with time. Thus it can be visualized as a
black box containing a device which has as its inputs, the
values of the State measured at finite number of points along

the Spatial domain and as its output the coefficients of the

83

specified degree best fit polynomial in some given sense (as an
example Chebychev fit of nth degree). This is illustrated in

Figure 5.0, for fixed t,

Q(x,t) ... an(c>q“<t> + an_1(t>q“'1(t) +...+ a1(t)q(t)

+ a (t) (5.9)
0

Since the coefficients are time varying, the polynomial fit
is very difficult to perform. Another method could be to fit
the polynomial at each time. Since the coefficients are dif-
ferent at each time, this requires a large computation time at
each time interval and if the state variable is of higher
dimension than one, the polynomial fit is difficult to perform
partially because the theory of polynomial approximation is
not very well developed in higher dimensions.

In light of the above discussion, it is desirable to
obtain a feedback control in terms of measurements made at a
finite number of points in the spatial domain. Now let uc(t)

be denoted by the feedback control law which is written as,

new = KrF(Qr(t) ,erun (5.10)

1

"MD

r
1Jhere Kr is either zero or an unknown matrix and F is a
Stiitably chosen function of the State vector and the desired
stzate er(t) discretized by a finite difference Scheme as

discussed earlier. Thus, given

3Q1
3:- = Gi(Qi(t)3Qii11(t)9°-oan_iIn(t):t), 1 = 1,---3n

84

 

qlit) J

320:) a

 

 

 

v

qj(t)

 

 

qn(t)— fl

B LACK
BOX

 

4800;)

 

931(t)

 

+8200

 

 

 

a an(t)

Figure 5.1 Illustration for the Discussion in Section 5.3.

Q(x,t ) as an(t )qn(t ) + an_lqn-1(t ) +...+ a1(t)q(t)+ ao(t)

(t

fixed)

85

along with the boundary conditions

Qi(t=0)=Qio xien,t20
SbiQi(t) = ui(t) xi 6 ob, t 2 0
‘ui(t)| s 1

to find the optimal control of the form given in (5.10) such

that the following performance index is minimized.
T
C = gEQd(x’tf) ' Q(x,tf)] [Qd(xstf) ' Q(xstf)]d0

which by the same discretization scheme, becomes,
n

c as '5 = 210d<xi.tf> - Q(xi.tf>3T[Qd<xi.tf><2(xixfn (5.11)
i=1

xien

There are in general two ways of obtaining the feedback para-
meter Kr’ given in (5.10). One of the methods is to obtain
Kr such that (5.11) is minimized. The second method is to
obtain Kr such that
t
f * T *
f [u (t) - uc(t)] [u (t) - uc(t)]dt (5.12)
o
*
is minimized, where u (t) is the optimal open-loop control
obtained by one of the existing methods. The first method is
discussed in this chapter.
To Simplify the subsequent derivations, a linear dif-

fusion system will be considered. Thus given the system,

 

03
at

with the b1

53
3x1
hi
0" |

Q (X

The perform

Discretize

(5'4) 3 the

 

where

where Prime
are n x n
tridiagonal

u(t

Then inCOrp‘

86

2
3%.: h_%.’ a = [0,1]
5x

with the boundary conditions

g§'x=o = a{q(0.t) - u(t)} (5-13)
33. =
ax x=1 O

q(x,0) = 0 , and with the constraint,

0 s u(t) s 1

The performance index to be minimized is

1
2
J = [[qd<x> - q(x,cf>] dx (5.14)
0

Discretize (5.13) and (5.14) using the scheme explained in

(5.4), the discretized system can be written as,

x.
II

A x +~b u(t) X(0) = x (5 15)

0

where

X = [q1(t).q2(t),.-..qn(t>l'

where prime indicates the transpose. The matrices A and b
are n X n and n x 1 respectively. The matrix A is of the
tridiagonal form.
Now assume
n
u(t) = u (t) = 2 K q (t) =<K.X(.t)> (5-16)
c i=1 r r

Then incorporating (5.16) into (5.15), yields

x=Xx,xw)=x GAD

O

87

where A contains the unknown feedback parameters to be
determined. The solution of (5.17) with the given initial

condition is,

X(t) = eAt xo (5.18)

assuming that the eigenvalues of A in terms of the unknown

parameters are known. Expand (5.18) into a constituent matrix

expansion to get,

= , 3 ’00., 5.1
X(t) F(K1,K2,K ,Kr,t)XO f(K1,K Kr) ( 9)

3"' 2

Similar use of the discretization scheme for the performance
index in (5.14) gives,
“ 2
AJ = .2 [qd(xi) - Q(Xi,tf)] (5.20)

1=1

Substituting (5.19) into (5.20), we obtain

n
_ _ 2
AJ .. iE].\:qcl(xi) f(K1,K2,...,Kr)] (5.21)

Since uc(t) is constrained to be in the limits

0 s uc s 1 (5.21a)

Now (5.21) is minimized with a given constraint (5.21a), and
the desired distribution qd to obtain the parameters k1,k2,...,kr.
Thus the problem is reduced into a parameter optimization
problem in a parameter space of n-dimensions.

Since the matrix A in (5.15) contains unknown con-
stants k1,k2,...,kr, the eigenvalues are very difficult to

obtain in terms of k k ...,kr, for large matrices. An:

1’ 2’

88

alternative method is to implement this method on the hybrid
computer. Several automatic parameter optimization schemes
with differential constraints are discussed in the literature.
Some of these methods are discussed in Bekey and Karplus (B-lO).
In most of these methods, the differential equations
are simulated on the analog computer and the computation of
the gradient and adjustment of the parameters is done by the
digital computer. Since the number of integrators available
on an analog computer are limited, this method limits the order
of Spatial discretization. So the decomposition principle
discussed in Chapter IV is used to increase the capabilities
of the analog computer. This enables one to solve a higher
order differential equation than that is usually possible with
the available integrators. The flow chart of this algorithm

is given in Figure 5.2.

5.4 Computer Results

 

As an example, the following problem is considered.

Consider

- 133:3—3- (5.22)
at ax

where q(x,t) is the temperature distribution in the metal

in dependence on the Space coordinate x, (0 s x s l) and

time t (0 s t s T). The Space coordinate x is normalized
with reSpect to the thickness of the metal and t is normalized
so that the coefficients corresponding to the thermal dif-

fusivity is unity. The initial and boundary conditions are given by

89

 

Start

 

Approximate the partial
Differential Equation into ordinary

Differential System

 

 

X = Ax +'Bu

J

If the system matrix A

 

 

is not tridiagonal,
tridiagonalize it by

standard procedures

J _

Decide the number

 

 

 

of partitions N

i

Choose the form of the

 

 

control law

“C(t) = E krF (Qr3er)

 

 

 

 

Incorporate this in the

system equation

1

Approximate the

 

 

performance index J

Figure 5.2 Flow Chart for Algorithm - I

 

90

[Guess k
r

' :4

Apply decomposition

 

 

 

 

 

principle

Parameter

 

adjustment Subroutine

 

 

 

New k.
1

 

No_J Check to see if
the required performance

is obtained

 

 

 

Yes

 

Print optimum ki and J

 

 

Stop

Figure 5.2 Flow Chart for Algorithm - I

91

Q(X,0) = O

ggvfl=aMWJ)'VQH can
as. =

ax x=l

where a, the heat transfer coefficient assumed to be constant,
v(t), the temperature of the gas medium is controlled by the

fuel flow u(t) and they are related by

r §§'+-v(t) = u(t) (5.24)

where r is the time constant of the furnace, and u(t) is
normalized properly.

The problem here is to obtain u(t) (O s t s T) such that

1 * 2
Itucc>3 = jtq (x) - q(x,T)] dx (5.25)
O

is minimized.

Furthermore u(t) is constrained to

0 s u(t) s 1 (5.26)

The various constants in the above problem are:

a=10 OSXSI

*
r = 0.04 q (x) = 0.2

Now we discretize the system using the scheme mentioned in

section 5.2, to obtain

flu)" {10.25 0
q1(t), 0.5 -1.5
92(t) 0 1

d
32- = 100 (::)
1:190:31

 

 

 

92

0 0 .
1 0
-2 1 0
‘\\\::::::\\0
0 0 1 -2 1
0 0 0 1 -1

9

u(t) = ucu) = :1 qu<xr.t>

where k
r
number of points to be sample

points, then (5.27) becomes,

d/dt (X(t)) A

X(0)

II
>4

where, X(t) = q(t)/100, and A

FL0.25

0.5

A(k1,k2,...,k9) =

 

in

V(t) 1
q1(t)
q2(t)

 

99(tt

25

u(t)

 

 

LO.

is either a constant gain or zero, depending on the

d. Suppose we sample at all nine

(k1,k

q
k1 k2 k3 . . kg
-1.5 1 0 . 0
1 -2 1 . . . 0
<:::) 1 -2 1
0 1 -1.J

2,...,k9) X(t)

is given as follows.

 

Now to apply the decomposition principle, the matrix is

partitioned as follows.

(5.28)

(5.29)

(5.30)

(5.27)

 

93

where A A12, A21, and A22 are 5 X5 matrices. Also

11’
let
= l I

where prime denotes transpose. Then (5.28) can be written as

(5.31)

From the matrix given in (5.29), the decomposition principle
requires the Storage of all the functions in the vector Y2,
but a little modification using the superposition principle

avoids this difficulty.

Let

v(t) = v1(t) + v2(t)

where

dv t
-dé_l a v(t) +k1 q1(t) + k2 q2(t) +...+k9q9(t)

Then writing

d

EE-v1(t) = a v1(t) + k1q1(t) + k2q2(t) +...+ k494(t)
g: V2(t) = a v2(t) + k5q5(t) + ............ + k9q9(t)

two different partitions of the matrix in (5.29) are obtained

as follows.

 

 

 

 

 

-0.25 k1 k2 k
0.5 -1 5 1
Y1 = 0 l -2
0 0 -2
L.0 0 0 1
p
-2 l 0 0
l -2 l 0
0 l -2
A22 7 0 0 1 -2
0 O 0 l
Lk5 R6 7 k
0
0
and b = 0
0
Lo ..
and Y2 = A22 Y2 + b x

and thus storage of only three functions are necessary.

1 and Y2

Pvlml
ql(t)
1 q2(t)
q3(t)
Lq4(t).

vectors Y

 

 

The results are summarized below.

4

 

 

 

are as follows.

"95(t)0
96(t)
97(t)
q8(t)
99(t)

 

 

Lv2(t)4

 

(5.32)

(5.33)

The

(5.34)

Figure 5.3 represents

the verification of Sakawa's results using the hybrid computer.

The optimal control obtained by Sakawa for this problem is applied

and the resulting state is verified.

represent the State functions obtained with two feedback

Figure 5.4 and Figure 5.5

0.4

0.3

0.2

0.

l

  
   
     

 

AX

b.—.—.— .—'_O _ .—

95

Q (X ,T)

 

...T =0.2, u(t) = l

Desired Value

 

. —.t.|a'

Optimal Control
T = 0.2 sec.

=o.1 , 0$u(t)$1

Figure 5.3 Verification of Sakawa's Results.

q Om)

0.4

0.3

0.2

0.1

 

96

*
q (x,T)

 

Figure 5.4 Unconstrained Case
*
q (x,T) = 0.2

q(x,0) = 0.

and

0.4 u

0.3 q

0.1 _

 

97

~ I & II

Constrained Case

 

0
«Ir

0
hI“

0.2 0.3

0.4

5r

0.5

Figure 5.5 Contrained Case

0 s u
c

0 s u
c

$10

5 1

(I)
(II)

db

0.6

0.7

.L

0.8

9+-

T‘fir

98

parameters and the state in these cases is measured at x = 0.1

and x = 1.0 respectively. The different cases considered are,

i)

ii)

Unconstrained case, that is the control is not con-
strained. The time interval is (0, .4) and the result-
ing state is shown in Figure 5.4. The resulting per-
formance is 0.0105. The corresponding feedback co-

efficients are,

k1 = 55.16

k2 = 54.65

and the corresponding control is

uc(t) = k1(qd(0-1.t) - q(0.1.t))
+ k2(qd(1.0,t) - q(1.0,t))

Constrained case I, that is the control is constrained

 

to be within some prescribed limits. The time interval
of interest is (0,0.4) and the resulting state is

shown in Figure 5.5. The form of the control is,

uc(t) = k1(qd(0.l.t) - q(0-1.t))

+ k2(qd(1.0,t) - q(1.0,t))

and the performance obtained is 0.011. The correspond-

ing feedback coefficients are

k1 = 54.89

k2 = 64.03

The control is constrained to be

0 s uc(t) s 10

99

iii) Constrained Case II. The form of the feedback control

is as follows.

uC(t) = k1(qd(0.l,t) ' Q(O.1,t))

+ k2(qd(l.0,t) - q(1.0,t))

and the time interval of interest is (0,0.4) and hence

T = 0.4. The control is constrained as
O s u (t) s 1
C

The performance obtained is 0.011, and the correSpond-

ing feedback coefficients are,

k1 = 59.34

k2 = 52.47

The resulting state is shown in Figure 5.5.

5.5 Sensitivity Considerations

In the above, the feedback coefficients are assumed to
be constants. In general the coefficients are time varying.
So the performance index is sensitive with respect to the
initial conditions and final time for constant gains. This
is illustrated by the following example. A method of obtaining
these sensitivity coefficients is given, and the extension of

the method to the general case is Straightforward.

5.5.1 Sensitivity coefficients

 

Let us consider a one dimensional diffusion equation,

2
3% x: AA;- , n a: [0’1] (5.35)
ax

 

100

and the boundary conditions

fikw=amm¢>-wu)
(5.23)
33 = O
ax x=l
q(x,0) = 0
and v(t) is given by
21' ‘ 2
r dt + v(t) — u(t) (5. 4)
0 s u(t) s l
and it is required to minimize
1 *
2
I[u(t)] =3“ (q (x) - q(x,T)) dx (5.25)
o

Discretizing (5.22) through (5.26) using the scheme explained

in (5.14), the discretized system can be written as

X = A X + b u(t) , X(O) = C

Using (5.16) uC(t) = u(t) = <K, X(t)> then

x = A x + b<K,X> (5.35)

where <--> is the Scalar product.
For given ki’ the performance index is sensitive to both
the final time and the initial conditions. So,

1 * 2
I[u<t)] = ([9 (x) - q(x,T)] dx
0

or (5.36)
n

AJ1u<t>1 = z <q*<x,> - xi<r>>2
' 1

1:

“-

 

101

Now (5.35) gives the solution

x (T) = eBT c (5.37)

where B = (A +-bKt).
Substituting (5.37) into (5.36) and writing q(xi,T) = Xi(T)

n * 2

AJEC.T] = 2 (q (xi) - q(xi,T))

i=1
where AJ = AJ(C,T) to emphasize the dependence of AJ on
C and T. Now the sensitive coefficients of AJ with reSpect
to C and T are given by
k = 1,...,n

n
3A1 . * - as.
2 2 (q 9%) q(xi,T)) ack

aCk i=1 i = 1,...,n

n * aq.(x,,T) k = l,2,...,n
34% 2 2 (q <st - q(xia» —-¥-—1——
i=1

5T 5T i = l,2,...,n

Thus writing

*
e. = q (X,) ‘ Q(X,,T), i = 1:29°°°:n
1 1 1

then
E = (e ,e ,...,e )
1 2 n (5.38)
x = (qlaq2:°°°3qn)
we have
P W
2‘1
aCk
8X2
BAl=2E'3—X , a—X = —— (5.39)
30k 80k 60k ack
93
C
315.1

 

 

102

 

 

,
2’51
3T

30:13:23'33. 53.: :

T 5T ’ 3T
935.
LaT ..

Also from Eq. (5.37)
X(C,T) = eBT - 0

Therefore

a! = B'eBT-C = B eBT c = B X(T)

5T
r3151 3:;
acl acn
an BX“ .

 

 

80 (5.39) and (5.40) become,

3Al'= 2 E' ° eBT
3C

J n
§%-=2E °BX(T)

(5.40)

(5.41)

(5.42)

(5.43)

(5.44)

and represent the sensitivity of the performance functional

with respect to the initial conditions and the final time

5.5.2 Computational algorithm

T.

The sensitivity coefficients given in (5.43) and (5.44)

can be obtained by known eBT and X(T). But computation of

103

eBT by using the series

eBT = (AT)j/j!

or finding the constituent idempotents is quite cumbersome,
due to the necessity of finding the eigenvalues. However,

the following algorithm can be used. Now write (5.43),

 

 

 

 

q n
“.311 Faxl 951 f1
5C1 501 3C2 aCn
301 ax2
3C2 8C1
BA; 3 - 2 E'1
3C
311 3’31 3‘2 Big
C C C C
1.6 n.J b3 l 3 2 B nJ
where eBT 3%,
5C
Now X(T) = eBT-C
and C = (C1,C2,...,Cn) and hence we can write the
following
’ r
9:1 3’: 6X11 ‘
5C1 3C2 3Cn l
2
M) = 3’12 952 3’2 .
acl ac2 acn
BX“ ax“
7:" 35' On
La 1 n ..J L d

 

 

 

 

104

Now sequentially setting the initial condition

C1 = (c1,0,0,...,0)
2

C - (0,c2,0,...,0)
cn = (0,0, ...... ,c )

and integrating the system equations (5.35) for given values

1 2
to get X (T), X (T) ,..., where

of k1,k2,...
r33317 ‘35
ac1 ac2
X1(T) = . , X2(T) = . etc.
L5C1 .302,

 

 

 

 

Thus the sensitivity coefficients given in (5.43) are then

3%; = 2 E' (X1(T):X2(T),...,Xn(T))

Finally it can be seen that (5.43) and (5.44) for
sensitivity coefficients depend on the error at the final
time and hence conclude that if the state is reachable within
the specified time T, the coefficients obtained for this T
are approximately valid for all T. The approximation will

be good because the error vector E will be close to zero.

CHAPTER VI

DEVELOPMENT OF ALGORITHM - II

6.1 Introduction

 

In this chapter a description of Algorithm-II for
obtaining an optimal feedback-control for a class of dis-
tributed systems is given. This differs from Algorithm-I
where the feedback coefficients are assumed to be constant
with reSpect to time in contrast to time varying gains
obtained in Algorithm-II. In this algorithm, the optimal

open-loop control is assumed to be available.

6.2 Problem.Formulation

Consider a distributed system characterized by a

vector partial differential equation

3% = cit,x.Q(x,t) .Qx(x.t).Qxx(x.t),u(t),m(x,t>]

along with the boundary conditions

Q(x,to) =QO(X) x e 0. t €10,ch

SbQ(Xb,t) = u(xb,t) t e [O,tf], xb 6 ob (6.1)

where u(t) is the boundary control, Q(x,t), the vector of
state functions and the performance functional to be minimized

is

105

106
c 9 g10d<x¢f> - Q(x,tf>]T[Qd(x,tf) - Q(x,tfndn
with u(t) is constrained such that a s u(t) s b.

6.3 Algprithm - II

Loosely speaking, the algorithm involves measuring the
state vector at any finite number of points in the spatial
domain and then obtaining time varying feedback coefficients
such that the feedback control so obtained is close to the
optimal open loop control obtained in some given sense. Since
the closed-loop control obtained results in a degradation of the
optimum performance this control is called sub-optimal feedback
control.

Let Q(B,t) be the State vector measures at x = a,

X E n GiRn. Then let the feedback control be represented by

Uc(t) = F(Q (Bat) :Qd(eat):K(t)at)
where uc(t) is an r-dimensional vector,

X(t) is a r X n matrix with an off-diagonal term zero,
Q(B,t) is the state vector of n x 1.
For the following discussion, assume uc(t) in the following

form:
uc(t) = Ktcc>10(e.t> - Qd(B,t)] (6.2)

Then the time varying coefficients can be obtained by minimiz-

ing one of the following

0 Qgtodeef) - Q(x,tf)]t[qd<x.cf) - Q(x,tfndn (6.3)

107
A tf * t *
1 =3" [u (t) - uc(t)] [u (c) - uc(t)]dt (6.4)
o

where Qd(x,tf) is the desired distribution and u*(t) is

the optimal open loop control. Here the second type of func-
tional is used Since the aim of this algorithm is to obtain
the feedback controls utilizing the open-loop control.

Several methods exist for obtaining the open loop
control for different classes of problems. Two methods called
quadratic programming and direct search method are discussed
in section 3.3 of Chapter III. Having obtained the open-loop
'control by one of the above methods, the method of obtaining
feedback control is discussed.

First the system characterized by the partial dif-
ferential equation is discretized in the Spatial domain either
by truncating the higher order terms of the correSponding
integral equation using integral transform techniques or by a
finite difference method. For parabolic systems s-domain
approximation is very effective. The finite difference
approximation is discussed in section 5.3. Now to determine

the coefficients, we write

ij(t) = igl aijL(ti_1,ti) , j = 1,...,r (6.5)
where
1 t, s t s t.
L(ti_1,ti) = 1'1 1 (6.6)

0 otherwise

The feedback coefficients are approximated by piecewise con-

stant functions where N is the number of subdivisions of

108

the time interval (0,tf).

Hence

(L)

uéL) (t) L = 1,...,r

(t) = Ru“) 8

where (6.7)
emu) = tq“)(3.c)- q“)<a.t>]

Substituting (6.5) into (6.7),

N
6:0 (t) = Elana“) (t) ti_1 .< t s ti
1 (6.8)
0 otherwise

for all i = l,2,...,N; L = 1,...,r.
Substituting (6.8) into the given performance functional in (6.4)

t

= I f[u*(t) _ uc(t)]T[u*(t) - uc(t)]dt is written as
O
t N
"' f *(L) (L)
I= 2(u (t)- 20.10.1396) (t))2dt
g L=l i=1 1L 1 i
N t, r
= z t1 { z (u *(“m- e000»)2 }dt
i=1 1-1 i=1
N Iti r
= .3 5 {z 01*“) (t) - 01. 3(L)(t))2_}dt
1=l ti-l L=1 1L

Thus the problem is divided into N-independent subsystems

where the parameter optimization is performed and the number of
parameters are equivalent to the number of state functions, at
each Stage. The initial condition for the ith stage is the

final value of the state vector of the (i-l)th Stage.

109

6.4 Example

Consider a slab of material bounded by the planes x = 0
and x = 1 which is in contact with a heat transfer medium of
temperature u(t) at x = 0 and is perfectly insulated at
x = 1. The dimensionless slab temperature q(x,t) is governed

by the one-dimensional heat equation,

2
as =13-
at 3X2 (6.10)
q(X.0) = 0
33.160 = 6<q<o,t> - u(t))

(6.11)

as. = O
ax x=1

The optimal control problem consists of determining u(t),

0 s t s tf, tf specified, to minimize the integral average

deviation of the temperature at t = tf, from a desired dis-

tribution qd(x), namely, to minimize,
1 2
P = g [qd(x) - q(x,tf)] dx (6.12)

vvith constraint on u(t)

0 s u(t) s l (6.13)

Faere q (x) = 0.2, t = 0.2, a = 10.
d

f
The optimal open-loop control for the above problem is obtained

by using direct search on the performance index (see section

3.3.2), and is shown in Figure 6.1. The feedback-control is

110

 

 

 

 

 

 

*
u (t)
1.0
0.8 "
0.6 WI
0.4 "
f
0.2 “
i i- Q C
0 0.04 0.08 0.12 0.16 0.2 ,
time --9

Figure 6.1 Optimal Open Loop Control by Direct Search.

0 s u(t) s l

lll

obtained as follows. Applying the Laplace transform to (6.10)

subject to the initial condition yields,

2
L3”: sQ(x,s) (6.14)
ax

where Q(x,s) = £(q (x,t)). Similar transformation of the

boundary conditions yield,

PIX-=0 = 6{Q<o,s> - 6(3)} (6.15)
X S =
Max ‘x=1 0 (6.16)
where u(s) = £(u(t)).

The general solution of (6.14) is

Q(x,s) = C1(s)sinh ([3 x) + C2(s)cosh ([5 x) (6.17)

where C1(S) and C2(s) are arbitrary functions of 3. They
are determined such that the general solution (6.17) satisfies

the boundary conditions (6.15) and (6.16). Thus

mm = c203)

59* = les cosh [S x + 0/3 sinh ,/s x

ax

Therefore
§x=0 = Crfs = O’{Cz ' “(5” (6°18)
Bxix=1 = les cosh f8 + CZ/S sinh fs = 0 =

Cr/s = -CZ/s tanh(/s (6.19)

112

Equations (6.18) and (6.19) yield

-a u(s) sinh [s

 

 

 

C1(s) =/s sinh [S + a cosh fs (6’20)
= a u(s) cosh JS
C2(S) fs sinh fs +a cosh fs (6'21)
Thus (6.17) gives,
_ a cosh(1-xL/s
Q(x,s) - u(S) [S sinh ,/s +0! cosh [3 (6°22)
or
Elke). = COSMl'XVS (6.23)

 

u(s) gs; sinh [S + cosh [3
Now before proceeding further, the following lemmas are proved.
Lemma 1. The equation
cosh z +'Bz sinh z = 0 (6.24)

has only imaginary roots and if z = x + iy, the roots are

the solutions of the equation

y tan y = 1/B (6.25)
“2522:. Given

cosh z +'Bz sinh z = 0 (6.26)
implies cosh z = -Bz sinh z and now consider different

cases.
Case 1. We know 2 = x + iy and x # 0, y # 0, then (6.26)

can be written as

coth z = -Bz

113

If 2 = x + iy, then

sinh 2x - i Sin 2y
cosh 2x - cos 2y

 

coth (x + iy) = = -B (x + iy) (6.27)

Equating the real and imaginary parts on both sides,

 

 

1 sinh 2x

_._ = .2
x cosh 2x - cos 2y B (6 8)
1 sin 2y
— = .2
y cosh 2x - cos 2y B (6 9)

Equations (6.28) and (6.29) when equated yield,

-(1/x)(sinh 2x) = (1/y)(sin 2y)
i.e. letting p1 = 2x and p2 = 2y, we have
-(2 sinh P1)/p1 = (2 Sin pz)/p2

Both the left hand Side and right hand Side functions of (6.30)
are even functions of p1 and p2 respectively. They are
plotted in Figure 6.2 (a) and Figure 6.2 (b). It is clear that
there is no p1 and p2 to satisfy the above equation. Hence
there are no roots with x = 0 and y = 0.

Case 2. x = 0 and y ¢ 0, then (6.26) becomes

cosh iy +1i By sinh iy = O

or cosh iy +'i By(i sin y) = 0
or cos y = By sin y
or y tan y = l/B

which is exactly equation (6.25). The plot of tan y = l/By

is shown in Figure 6.2 (c).

114

 

 

 

 

 

 

 

 

. -2 sin
2 Slnh pl/p1 p2/p2
P2 "
2 0 -2.0
- 4L
I0 p1 0 p2
Figure 6.2a Plot of 2 sinh pl/p1 Figure 6.2b Plot of -2$in pz/p2
tan y
l/By
‘52 '31 I “
1 fl' _, fit- 1 .w .__
I 61 2 82 y
1/By

 

Figure 6.2c Plot of tan y = l/By

115

Case 3. x i 0, y 0, then (6.26) becomes

cosh x + Bx sinh x = 0

or cosh x = -Bx sinh x

and Since x # 0,

coth x -Bx

For positive B, there are no intersection points and hence
no roots, for x # 0 and y = 0. (See Figure 6.2 (d)).
Thus the only roots of (6.24) are imaginary and they

are the solutions of

y ta“ y = “B 0.3.1).

Lemma 2. The roots of cosh B2 = 0, B # 0, are completely

imaginary and are given by
y1 = i(21 - 1)/%§ , i = 1,2,3,...
Proof. Let 2 = x + iy, then
cosh B2 = 0 can be expanded in the following way,

0 = cosh B2 = cosh Bx cosh iBy + sinh Bx sinh iBy
or

cosh Bx cos By + i sinh Bx Sin By = 0
Equating the real and imaginary parts, we have the real part,
cosh Bx cos By = 0 implies cos By = 0.

Therefore y1 = i 1/B(2i - I)?- 1 -- 1,2,3,... (6.31)

and the imaginary part,

sinh Bx sin By = 0 implies from (6.31) that

116

 

 

   

-3x (B > O)
coth x

 

Figure 6.2d Plot of coth x = -BX (for B > 0).

117

sinh Bx = 0 implies Bx = 0 and hence x = 0 since

B # O by the hypothe81s. Q.E.D.

Now using the above lemmas the equation (6.23) can

be written at x = a, as

H (1 + s/y?)
Q(g,s) = i=1 1

 

 

, 0 S c S 1
u(S) co 2
n (1 + smi)
i=1
1 . . _
where yi =1: (1'0) (21 - l)n/2 , 1 - 1,2,...

5- are the roots of
1

8 tan 8 = a.

As can be seen from Figure 6.2 (c), the 51's increase very fast
and since their squares occur in the denominator, the infinite
product can be approximated by a finite product. Similarly by
choosing the point a as close as possible to the end point,

I

the roots yi s can be made very large, and thus approximating

the numerator by a finite product. Thus,

 

 

(1 + s—2>
Q(033) = y1 (6.32)
“(3) (1 +55) (1 + 137)
B1 B2

is quite a good approximation, i.e.

(1 + 35x1 +%>Q<e.s> = <1 + %>u<s>
B1 62 Y1

118

'Transform.this into time domain, and let v1(t) = Q(o,t)

 

 

V 9
212 + 212 (Bi + 6:) +‘V1(t) = u(t) +‘l§’&(t)
B152 B152 Y1
or

8282

. 2 2 2 .
V1 + (5% + 8%))“. + 6182 V1(t) = BIB; U(t) + 122 U(t)
y1

vahere dot denotes the differentiation with respect to time.

IJsing the following transformation

 

21 = v1(t)
2 2
. B152
Z2 - vl‘- 2 u(t),
y1

t:he equations given in (6.33) will be transformed into

2 o 22 b2

2
vahere a = 5152

m
I—‘
I
A
m
t—‘N
+
U)
NM
V

110w assuming,

u(t) = uc(t) = k(t)[Q(o,t) - Qd(o,t)]

N

k<t> = z ai(t)L(ti_1,ti) . [ti_1,ti] e [O,tf]
i=1

where

then

Thus (6.9) becomes

H!
II
IIMZ

i 1

Thus each “1

Thus the feedback gain k(t)

constant on sub-intervals and N

119

t, Stst,
1-1 1

otherwise

otherwise

f i <u*<t> - a, e<t>>2dc.
ti-l

can be obtained independent of the other a-

is approximated by a piecewise

can be chosen in an iterative

manner until the desired performance measure is obtained.

The next section gives the computer set up and the

results.

6.5 Computer Results

The sample value is taken at

model becomes,

v1 0

. = 2 2

v2 '5152
0
-52.98

where v(O) = 0.

x = 1. Hence the system
1 v1 0
+ MD
2 2 2 2
1 v1 0
+ u(t)
-21.3976 v2 52.98

120

The time interval of interest is (0, 0.2) and the analog
computer set up is shown in Figure 6.4.
We Start by dividing the interval (0, 0.2) into four

Sub-intervals as follows.

m

a1 belongs to t (0, 0.1)

(T\

02 belongs to t (0.1, 0.12)
a3 belongs to t E (0.12, 0.18)

belongs to t (0.18), 0.2)

m

0’4
Then the algorithm follows by setting the initial conditions on
the integrators 200, 201, and 241 to zero. Then the constant
a1 is obtained by one dimensional search by keeping the analog
computer in the repetitive mode. Having obtained the constant
the initial conditions are set up by using the digital computer
which are the final conditions for the first Stage. Then a
one dimensional search yields the second constant a2. Similarly
a3, etc. are obtained by following the above procedure. The

constants are as follows.

a =5.76 , t6 (0, 0.1)

1

62 = 13.57 , t e (0.1, 0.12)
a3 = 0.0 , c e (0 12, 0.18)
04 = 51.49 , r e (0.18, 0.2)

and i = 0.0872.
Then in the next step the interval is divided into

Six subintervals. The subintervals and the constants are,

121

A<HHUHQ

.mHQmem ecu u0w a: bum uou5a500 panama m.o muswflm

 

 

 

 

HOEHZQD

coma :5
Mmhbmzoo A<HHUHQ

 

 

 

 

 

 

 

 

_o H+ MNMOm h ,
_ - - Ham Hon
e... N

u . ODN

 

 

 

 

 

 

 

 

 

a1 = 5.52 ,

a2 = 7.09 ,
a3 = 10.38 ,
a4 = 16.52 ,
05 = 0.0 ,
a6 = -12.05 ,

and f = 0.0215.

122

t e (0, 0.05)

t E (0.05, 0.07)

n
m

(0.07, 0.1)

n
m

(0 1, 0 12)

c e (0 12, 0.18)

n
m

(0.18, 0.2)

Thus by increasing the number of intervals, the per-

formance index I can be made smaller and smaller, till it

satisfies the required performance.

CHAPTER VII

CONCLUSIONS

7.1 Conclusions

In this thesis, computational methods have been
developed to obtain the optimal feedback controls for a class
of distributed parameter systems. The class of problems con-
sidered in this these are "well posed" (see Section 1.4.1) and
possess the following properties.

a) Solutions to the system equations
3% = (m(x,t),m<x,t>.u<t>.x,t3

Q(x,to) =QO(X) ; x 6 n: Rn
SbQ(xb,t) = U(t) ; xb 6 (lb
with the given boundary conditions exist.

b) The solutions are uniquely determined.

c) The solution depends continuously on the initial
data. This says that small changes in the initial data will
cause correspondingly small changes in the solution, Q(x,t).

The algorithms for obtaining the sub-optimal feedback
control are discussed in Chapters V and VI. This is called sub-
optimal feedback control because the application of these feed-
back control laws result in a degradation of optimum performance.

In both the algorithms, the problem is reduced to a parameter

123

124

optimization problem with differential constraints. The a
priori information available about the optimal open loop con-
trol is used in the second method to obtain the time varying
feedback gains. The methods are easily implemented on the
hybrid computer. The hybrid system available in the Hybrid
Simulation and Control Laboratory (IBM 1800-AD-4) was used to
obtain solutions for the examples in the thesis.

The capability of obtaining a solution for a dif-
ferential system on an analog computer is limited by the
number of integrators available on a given facility. A decom-
position principle, which decomposes a large set of differ-
ential system equations into lower order independent Sub-
systems which are solved iteratively is described in Chapter
IV. (The convergence theorems are Stated and proved. With
this treatment, a larger system (a finer Spatial discretiza-
tion) can be considered which would not be feasible otherwise.

Thus a significant contribution is made in this thesis
in the area of distributed parameter systems by developing some
efficient computer algorithms for obtaining feedback-controls
and solving some of the problems encountered in the actual

implementation on the computer.

7.2 Possible Extensions

In this thesis linear systems were considered, but the
methods can be extended to non-linear Systems. These non-linear
problems must be well-posed. The verification of these conditions

for non-linear systems are very difficult. It may be possible

125

to apply linearization techniques about a nominal trajectory

in applying the above methods to some non-linear systems. These

yield approximate results. Another possible extension is to find
a class of problems where the state can be approximated by small

order polynomial fit so that the results available in the lumped

case could be applied.

This thesis emphasizes the fact that the results obtained
in the case of lumped parameter systems cannot be applied
directly for distributed parameter systems and thus new results
obtained in this thesis are necessary. In these lines the thesis
can be extended by changing the performance index such that the
number and location of the measuring instruments along the
spatial domain are optimized while penalizing the system for
using large number of sensors. For solving these systems de-
tailed comparison of the results, if possible, obtained by
using approximate techniques are desirable. A listing of the
best approximations for reducing several of the infinite
dimensional systems which are common to finite dimensional

systems will be very helpful.

REFEREN CES

B-3

B-4

B-5

B-8

B-9

REFERENCES

Axelband, Elliott, I., ”An Approximation Technique for

the Control of Linear Distributed Parameter Systems with
Bounded Inputs", IEEE Transactions of Automatic Control, I
Vol. Ac-ll, pp. 42-45, January 1966.

Athans, M. and Falb, P.L., "Optimal Control", McGraw-

 

 

Hill Book Co., New York, N.Y., 1966. 5
Brogan, W.L., "Dynamic Programming and a Distributed '
Maximum Principle", Proc. JACC, 1967.

Brogan, W.L., "Optimal Control Theory Applied to Systems
Described by Partial Differential Equations", Advances
in Control Systems, Vol. 6, 1968.

Berg, IRW. and McGregor, J.L., "Elementary Partial Dif-
ferential Equations", Holden-Day publications, San
Francisco, 1966.

Butkovskii, A.G., "Optimum Processes in Systems with
Distributed Parameters", Automation and Remote Control,
Vol. 21, pp. 13-21, 1961.

Butkovskii, A.G. and Larner, A.Y., "The Optimum Control
of Systems with Distributed Parameters", Automation and
Remote Control, Vol. 21, pp. 472-477, 1960.

Butkovskii, A.G., "The Maximum Principle for Optimum
Systems with Distributed Parameters", Automation and
Remote Control, Vol. 22, pp. 1156-1169, 1962.

Butkovskii, A.G., "The Broadened Principle of the Maxi-
mum for Optimal Control Problems", Automation and Remote
Control, Vol. 24, pp. 292-304, 1963.

Butkovskii, A.G., "Some Approximate Methods for Solving
Problems of Optimal Control of Distributed Parameter
Systems", Automation and Remote Control, Vol. 22, pp.
1429-1438, 1961.

Bellman, R., "Dynamic Programming", Princeton University
Press, Princeton, New Jersey, 1957.

126

E-2

E-3

E-4

K-Z

K-3

K-5

127

Bekey, G.A., Karplus, W., "Hybrid Computation", McGraw
Hill, 1969.

Denn, M.M., "Optimal Boundary Control for a Non-Linear
Distributed System”, Int. J. Control, pp. 167, 1966.

Egorov, L., "Optimal Control by Processes in Certain
Systems with Distributed Parameters", Automation and
Remote Control (English Translation), Vol. 25, p. 613,
1964.

Egorov, 1., "Optimal Processes in Systems Containing
Distributed Parameters", Automation and Remote Control,
p. 977, 1965.

Eveleigh, V.W., "Adaptive Control and Optimization
Techniques", McGraw-Hill Book Co., New York, N.Y., 1967.

Ellsworth, W.C., "Hybrid Computer Solution of Linear

State Models", Ph.D. Thesis, Department of Electrical
Engineering, Michigan State University, East Lansing,
Michigan, 1969.

Fernundo, L. Alvarado, R. Mukundan, "An Optimization
Problem in Distributed Parameter Systems", Int. J.
Control, Vol. 9, p. 665, 1969.

Gelfand, I.M. and Fomin, S.V., "Calculus of Variations",
Prentice Hall Incorporated, Englewood Cliffs, N.J., 1963.

Khatri and Goodson, "Optimal Control of Systems with
Distributed Parameters", JACC, p. 390, 1965.

Katz, S., "A General Minimum Principle for End Point
Control Problems", J. Electronics and Control, p. 189,
1964.

Kim and Erzberger, "On the Design of Optimum Distributed
Parameter Systems with Boundary Control Functions", IEEE
Transactions on Automatic Control, p. 22, 1967.

Kim and Gajwani, "A Variational Approach to Optimum.Dis-
tributed Systems", IEEE Transactions on Automatic Control,
p. 191, 1968.

Koivo, A.J. and Kruh, P., "On the Design of Approximately
Optimal Feedback Controllers for a Distributed Parameter
System", Int. J. Control, Vol. 10, p. 53, 1969.

Kreindler, B. and Athans, M., "Optimal Control with Piece-
wise Constant Gains", IEEE Transactions on Automatic
Control, August 1968.

K-7

R-l

S-2

S-3

s-4

S-5

S-7

S-8

128

Kalman, R.E., "The Theory of Optimal Control and Calculus
of Variations", RIAS Rept. 61-3, Res. Inst. for Advan.
Studied, Baltimore, Maryland, 1961.

Leitmann, "Topics in Optimization", Academic Press.

Lapidus, L. and Luus, R., "Optimal Control of Engineer-
ing Processes", Blaisdell Publishing Company, 1967.

McCausland, I.J., "On Optimum Control of Temperature
Distribution in a Solid", J. Electronics and Control,
Vol. 14, No. 6, pp. 655-68, 1963.

Russell, D.L., "Optimal Regulation of Linear Symmetric
Hyberbolic Systems with Finite Dimensional Controls",
SIAM J. on Control, Vol. 4, p. 276, 1966.

Sakawa, Y., "Solution of an Optimum Control Problem in
Distributed Parameter Systems", IEEE Transactions on
Automatic Control, Vol. AC-9, p. 420, 1964.

Sakawa, Y., "Optimal Control of Certain Types of Linear
Distributed Parameter Systems", IEEE Transactions on
Automatic Control, Vol. AC-11, p. 35, 1966.

Sage, A.P. and Chaudhuri, S.P., "Discretization Schemes
and the Optimum Control of Distributed Parameter Systems",
Proc. Asimilar Conference on Circuits and Systems,
Montery, California, p. 191, 1967.

Sage, A.P. and Chaudhuri, S.P., "Gradient and Quasi-
linearization Computational Techniques for Distributed
Systems", Int. J. on Control, p. 81, 1967.

Sage, A.P., "Optimum Systems Control", Prentice Hall
Inc., New Jersey, 1968.

Seinfeld, J.H. and Kumar, K.S.P., "Synthesis of Sub-
Optimal Feedback Controls for a Class of Distributed
Parameter Systems", Int. J. Control, Vol. 1, p. 417, 1968.

Seinfeld, J.H. and Lapidus, L., "Computational Aspects
of the Optimal Control of Distributed Parameter Systems",
Chemical Engineering Science, Vol. 23, p. 1461, 1968.

Seinfeld, J.H. and Lapidus, L., "Singular Solutions in
the Optimal Control of Lumped and Distributed Parameter
Systems", Chemical Engineering Science, Vol. 23, p. 1485,
1968.

Varga, R., "Matrix Iterative Analysis", Prentice Hall
Inc., Englewood Cliffs, New Jersey, 1962.

W-3

129

Wang, P.K.C. and Tung, F., "Optimum Control of Dis-

tributed Parameter Systems", J. Basic Engineering, Trans.
of ASME, Vol. 86D, p. 67, 1964.

Wang, P.K.C., "Control of Distributed Parameter Systems",
Advances in Control, Vol. 1, 1963.

Wismer, D.A., "An Efficient Computational Procedure for
the Optimization of a Class of Distributed Parameter
Systems", J. Basic Engineering, Trans. of ASME, p. 190,
June 1969.

Wendroff, B., "Theoretical Numerical Analysis", Academic
Press, 1966.

Wilde, D.L. and Beightler, C.S., "Foundations of Optimiza-
tion", Prentice Hall Inc., New Jersey, 1967

S

"'Tl't‘uifl‘lfltfll‘rfilluifl(H.1111511111111“