SOLUTION OF A LINEAR REGULATOR
PROBLEM WITH QUADRATIC
PERFORMANCE INDEX AND STATE
VARIABLE CONSTRAINTS

Thesis for the Degree of Ph. D.
MICHIGAN STATE UNIVERSITY
HENRY STANLEY MIKA
1970

' mama? L

Michigan State
University .

       
    

This is to certify that the

thesis entitled
SOLUTION OF A LINEAR REGULATOR PROBLEM
WITH QUADRATIC PERFORMANCE INDEX AND
STATE VARIABLE CONSTRAINTS

presented by

Henry Stanley Mike

has been accepted towards fulfillment
of the requirements for 3

Ph D . 5.5.
___degree m______

 

U v Major, professor

N

Date ‘ ., ‘4' 0

0-169

 
  
 

   
    

y DINSING BY

'IIOAG & SIIIIS'
509K BINDERY INC.
1'... ' +2 tows

 

 

ABSTRACT
SOLUTION OF A LINEAR REGULATOR PROBLEM

WITH QUADRATIC PERFORMANCE INDEX AND
STATE VARIABLE CONSTRAINTS

By

Henry Stanley Mika

A direct method has been developed to obtain the optimal
linear constant control law for the linear, time invariant
regulator problem with quadratic performance index and state
variable constraints. For a scalar control and a diagonal state
weighting matrix Q, it is Shown that the constant gain matrix
K can be obtained without solving the n(n+l)/2 matrix Riccati
equations. Instead, it is demonstrated that a simple algebraic
algorithm defines all the elements of the K matrix given the
n elements forming the bottom row which are simply and linearly
related to the n coefficients of the closed loop characteristic
equation.

A computational algorithm is develOped which consists
of two basic parts: a digital program that locates the minimum
quadratic cost as a function of n elements of the K matrix;
and an hybrid program which checks and, if necessary, adjusts
these elements to meet the state variable constraints. This
adjustment is implemented by transforming the problem into an

n-parameter Optimization problem and using sensitivity techniques.

Henry Stanley Mike

The minimum quadratic cost search yields a global
minimum and, if the corresponding K matrix yields a feedback
control such that all state variable constraints are satisfied,
then this is the desired solution. If adjustment of the K
matrix is required to meet state variable constraints, then the
hybrid program determines an upper bound for the desired solution.
Thus the desired solution is bounded by the global minimum and
a known upper bound. Further search within this restricted region
is continued until the desired solution is obtained. The final
result yields the complete feedback loop engineering design in-
formation for the optimal control.

The method is illustrated with a solution of a third
order system which represents the automatic control of longitudinal

spacing between two moving vehicles with acceleration constraints.

SOLUTION OF A LINEAR REGULATOR PROBLEM
WITH QUADRATIC PERFORMANCE INDEX AND
STATE VARIABLE CONSTRAINTS
By

Henry Stanley Mika

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

(D

— bid:

V'

“".‘
/-- (d I.

.’

.a'
/.J

PLEASE NOTE:

Some pages have indistinct
print. Filmed as received.

UNIVERSITY MICROPILMS.

To my wife

Margaret

ii

ACKNOWLEDGEMENTS

In looking back at the course of events leading to the
publication of this thesis, the author finds himself indebted to
Dr. John B. Kreer, Chairman of the Guidance Committee and thesis
advisor. The successful completion of this work is in large part
the result of many stimulating discussions with him where his
insight and guidance proved invaluable. The author will always
recall this period with some nostalgia and a great sense of
gratitude.

The author wishes particularily to acknowledge the con-
structive connents and suggestions of Dr. Gerald L. Park and
Dr. Robert O. Barr both in the preparation of this thesis and
during the many occasions for personal discussion. Special
thanks are due Dr. Alex C. Bac0poulus for his interest in this
work and for his probing questions which the author found both
challenging and inspiring. Finally, appreciation is expressed
to Dr. Herman E. Koenig who, deSpite the pressure of many
activities, made himself available whenever needed.

The author wishes to express his gratitude to the Division
of Engineering Research, Michigan State University for the
financial support of this research.

It has been said that no man stands alone. Never has

this been more true than in this instance. Many times during

iii

the period of study and research the author's wife, Margaret,

has had to make personal sacrifices and share many near-traumatic
experiences. Her understanding patience and encouragement have
made this possible, and the author can only repay with his love

and gratitude.

iv

Chapter
I.

II.

III.

IV.

TABLE OF CONTENTS

ABSTRACT

DEDICATION

ACKNOWLEDGEMENTS

LIST OF TABLES

LIST OF FIGURES

LIST OF APPENDICES

INTRODUCTION

REVIEW OF EXISTING CONTRIBUTIONS BASIC TO
SOLUTION OF PROBLEM

Constraints

Controllability and Observability
Analytic Solution for Regulator Problem
The Inverse Problem

Sensitivity Analysis

DEVELOPMENT OF SOLUTION

Canonical Transformation

Gain Matrix

Autonomous Form

Minimum Cost Algorithm

Search Procedure

Constraint Function

3.6.1 Geometric Interpretation

3.6.2 C-Function Iterative Procedure
3.6.3 Discontinuity in p-Space

COMPUTER IMPLEMENTATION

4.1

4.2

Minimum Cost Algorithm Program
4.1.1 Incrementation

4.1.2 Computation

C-Function Algorithm Program
4.2.1 Criteria Checks

4.2.2 Analog Computation

V

Page

ii
iii
vii
viii

ix

11
13
15
23

27

27
29
32
36
41
47
50
55
59

63

63
71
71
73
73
77

VI.

omputation
stimation

4.2.3 C C
4.2.4 5 E
EXAMPLES

5.1 Example 1
5.2 Example 2

SUMMARY AND CONCLUSIONS

REFERENCES

83
83

86

89
96

100

103

LIST OF TAB LES

Table Page
I. Search for Absolute Minimum 90
II. Results of C-Function Search 93
III. C-Function Behavior under Peak Switching 98

vii

3-3

3-4

3-5

4-l(a)
4-1(b)
4-l(c)
4-1(d)
4-1(e)
4-1(r)
4-2(a)
4-2(b)
4-2(c)
4-3(a)
4-3(b>
4-3(c)

5-1

5-3

5-4

LIST OF FIGURES

Minimum Cost Search Procedure

hi Function

C-Function in Two-Space

Algorithm to Obtain the Minimum Cost

K Matrix

Peak 1 Dominant

Peak 2 Dominant

Minimum Cost Program

Minimum Cost

Minimum Cost

Minimum Cost

Minimum Cost

Minimum Cost

Program
Program
Program
Program

Program

C-Function Program

(con't.)
(con't.)
(con't.)
(con't.)

(con't.)

C-Function Program (con't.)

C-Function Program (con't.)

Solution of

x = Gx

Sensitivity Analysis

Analog Logic Control

Vehicle Longitudinal Spacing Problem
C-Function Iteration - Example 1

Step-size History of C-Function Solution -

Example 1

Solution Stability under Peak Switching -

Example 2

viii

Page
42
51

52

61
62
62
65
66
67
68
69
7O
74
75
76
78
79
80
87

92

97

LIST OF APPENDICES
Appendix Page
A. Sample Minimum Cost Algorithm Program 107

B. C-Function Algorithm Program 110

ix

CHAPTER I

INTRODUCTION

The last ten or fifteen years have produced significant
advances in the field of optimal control. However, deSpite the
current advanced status of the theory, a very limited number of
engineering applications have been made. For example, one of
the simplest and most useful optimization criteria is the quadratic
performance index which leads directly to determination of stable
closed loop systems, eliminating much of the "cut and try" effort
associated with conventional servomechanism design. In the case
of the linear system with the control u(t) unconstrained, the
general solution is well known. However, in contrast to the
general solution, if the engineering solution is defined as the
solution which yields specific design values for a stable feed-
back control and satisfies all state variable constraints, then
to the author's knowledge there exists no direct method to obtain
such a solution. In the method described in this thesis, the
general solution is used as a necessary condition to obtain the
engineering solution for a linear time-invariant system.

In general mathematical form, known as the Problem of
Bolza, the problem can be stated as follows:

Given the dynamic system

.2 = f[x(t) ,u(t) ,t]

and some performance index

tr

J = ¢(tf) +.It Lix(t).u(t),t]dt.
O

*
find the control u (t), t0 3 t S t which minimizes (or

f)
maximizes) the functional J. When the terminal "cost" m(tf)
is not a consideration, i.e., when

tr

J = It L[x(t),u(t),t]dt
0

it is known as the Problem of Lagrange.

To obtain u*(t), in general, it is necessary to solve
a two-point boundary value problem. In a few instances, it has
been possible to obtain analytical solutions; but in most cases,
particularily where constraints are imposed on the control u(t)
or on the State x(t), the only practical method of solution is
with the use of numerical techniques and a digital computer.

The choice of the form of the performance index J affects
the complexity of the solution to the two-point boundary value
problem as well as the complexity of the physical implementation
of the desired control u*(t). The choice may be obvious and
leave no room for flexibility as, for example, the engineering
specification that the system perform its function in minimum
time. However, if some flexibility of choice is permitted,
motivated perhaps by cost of physical implementation, then the
form of J is made on the basis of judgement, experience,
physical insight and mathematical tractability.

The complicated and sophisticated hardware often needed

for the implementation of an optimal controller may very well

be the reason why few actual applications of optimal controls
currently exist. One has only to consider the complicated feed-
back elements required for a minimum-time optimal control for a
third order system. The switching surface is represented by an
extremely complex mathematical function so that the control or
feedback elements may represent the major cost of implementing
the entire system.

In the case of linear systems
x = Ax + Bu,
the quadratic performance index
tf T T
J = g ft [x Qx + u Ru]dt
o

is an example of a functional which possesses many desirable

engineering features. It is mathematically convenient and leads

to a linear control law. In addition, it possesses the intuitively

appealing feature of "least-square-error" minimization used in

many areas of scientific and engineering analysis and application.

Specifically, it penalizes the system for large excursions from

some desired state trajectory while minimizing the control ”energy”.
There are other, less obvious, advantages to the use of

a quadratic performance index which mitigate against the use of

other performance criteria. For example, a smooth control u*(t)

is assured as contrasted to a "bang-bang" type, and the physical

implementation of the control u*(t) can be relatively simple

and inexpensive.

There is a large class of problems where a smooth control
is of paramount importance. For example, in the field of trans-
portation, including automobiles, aircraft, trains, etc., where
human comfort and performance are important considerations, it
is obvious that a "bang-bang" control of vehicle acceleration
would be highly undesirable.

The general solution for the optimal control of a linear
system with quadratic performance index and u(t) not constrained

is well known and is given by
* -l
u (t) = -R BTK(t)x(t).

The matrix K determines the feedback coefficients and is related
to the A, B, Q and R matrices through the nonlinear matrix
Riccati equation. Since the K matrix contains the feedback gain
information, it controls the behavior of the resulting closed
loop system. However, since this matrix depends on the weighting
matrices Q and R, it follows that the system trajectories in
state Space are a function of Q and R. Thus, when the designer
chooses a specific performance index, he also determines the
system behavior in state Space.

The relationship between the Q and R matrices and
the system trajectories is not readily apparent, and the follow-
ing trial and error design procedure is typical:

I) assume initial Q and R matrices;

2) solve the matrix Riccati equation;

3) simulate the system;

4) if the system performance is not satisfactory,

repeat the process.

This procedure is far from satisfactory because it requires
repeated solution of the matrix Riccati equation and there is
no systematic scheme for correcting Q and R.

The solution of the matrix Riccati equation may in itself
be a very difficult computational problem. It requires the solu-
tion of n(n+l)/2 simultaneous nonlinear equations, either in
the algebraic or in the differential form. In general, the matrix
Riccati is a first order differential equation which, when solved
backward in time, yields a limiting value for K which is a
constant matrix equivalent to the solution obtained for the
algebraic or steady state matrix Riccati equation. Satisfactory
computational methods for systems of moderately high order do
not appear to be available.

The recognition of the weighting matrices Q and R
as "design parameters" has led to the postulation of the inverse
problem: given K, what are Q and R? On the surface, this
seems to be a trivial engineering problem -- since knowing K
implies knowing the desired feedback gains which are the primary
engineering design results. This has been a major shortcoming
of the published results in this area: dependence on a-priori
knowledge of either K or closed loop eigenvalues. In addition,
solution of the matrix Riccati equation is necessary, sometimes
within an iterative loop. In summary, a practical engineering
solution has not been available.

In the following chapters a solution is presented to the
problem:

Given the linear, time-invariant system with scalar control

x=Ax+Bu
and a quadratic performance index

J=%

[xTQx + ruzjdt ,

OL—a8

find the constant matrix K which minimizes J, defines an
optimal control u*(t) such that no state variables Xi(t)
exceed their constraints, and does not depend on the solution
of the matrix Riccati equation.

The result is in the form of a computational procedure
utilizing both a digital and a hybrid configuration, and the
discussion is divided as follows:

Chapter II reviews the existing theoretical background with
Specific emphasis on the inverse problem;

Chapter III develops the theoretical foundation for the algorithm
used to obtain the solution;

Chapter IV describes in detail the computer implementation of
the algorithm;

Chapter V shows typical results for a third order system based
on the problem of automatic longitudinal Spacing between two
vehicles;

Chapter VI summarizes the results and Suggests future extensions
to this work.

All of Chapters III and IV are devoted to a detailed dis-
cussion of the original contributions; and any material which is
not original is explicitly referenced. Specifically, I) a Minimum

Cost Algorithm is developed which locates the arbitrarily small

neighborhood of the global quadratic cost minimum under the con-
straints that the Q and K matrices be positive definite and
the control is asymptotically stable; 2) a Constraint Function
Algorithm is developed which establishes an upper bound on the
minimum quadratic cost such that all state variable constraints
are satisfied and significantly reduces the search region which
contains the K matrix corresponding to this minimum; 3) a simple
K matrix algorithm is deve10ped, an integral part of the Minimum
Cost and Constraint Function Algorithms, which eliminates the

need to solve the n(n+l)/2 Riccati equations.

CHAPTER II

REVIEW OF EXISTING CONTRIBUTIONS
BASIC TO SOLUTION OF THE PROBLEM

The optimal control solution to the linear regulator
problem with quadratic performance index has been known since
1960 when Kalman published his results. In the past seven or
eight years recognition has been given the inverse problem, i.e.,
how to find the weighting matrices in the quadratic performance
index given an optimal control law. This chapter reviews these
results and other background material pertinent to the new
results to be discussed in later chapters. Included in this
chapter are brief discussions on topics of controllability and
Observability, state and control variable constraints and
sensitivity analysis.

In general, a linear time-invariant dynamic system may

be represented by the state equation

x.
ll

Ax + Bu (1)

and an output equation

Cx +-Du (2)

‘<
II

where x is the nXl state vector, u is the le, m s n, con-
trol vector and y is the er output vector. A, B, C and D

are constant and conformable matrices.

J

[C .
The optimal control u (t), to s t 5 t1, 18 one that

minimizes the quadratic performance index

t1

J = % Ito [XTQX +-uTRu]dt (3)
where Q is a state-variable weighting matrix which is either
positive definite or positive semidefinite; and R is the con-
trol variable weighting matrix which is positive definite.

The search for this minimum, if it exists, may be
implemented in several ways, for example, by the Maximum
Principle of Pontryagin [FUN-l] which is an efficient generaliza-

tion of the necessary conditions for a local minimum.

2.1 Constraints.

The solution of an optimal control problem inherently
involves constraints. For example, the solution of equation (3)
depends on the constraints imposed by equation (1). As in

ordinary calculus, equality constraints of the form
f[x(t),u(t),t] - x = o

are readily adjoined to J by a Lagrange multiplier or co-state
vector I, so that the problem is one of minimizing
t1 T T T -
J = 15 It {X Qx + u RU “I" )\ [Axi-Bu-x]}dt .
o

A more complex problem results when inequality constraints

are imposed on the control vector

U.(max) 2 ui(t) 2 Ui(min); i = 1,2,.....,m .

10

The usual procedure is to change the inequality constraint into
an equality constraint as first suggested by Valentine [VAL-1]

to form a function

2
Vi(ui) = [ui<max> - u,(c>][ui<t> - ui(min>] - g,<t> = o

and again adjoining with the Lagrange multiplier Y so that
t m
J = % jtl {xTQx + uTRu + XT[Ax+Bu-x] - .2 Vivi(ui)}dt .
0 1=l
The minus sign precedes the last term because Vi(ui) is
negative when constraints are violated and this must reflect
an added cost.
In the case of state variable inequality constraints,
augmentation of the state vector is implied. Several approaches

[REL-1], [BER—l] of which the following [McG-l] is representative

have been successfully applied.

L

. - _ = 2
Define xn+1 — fn+1 i=1[hi(x,t)] H(hi)

2
where [hi(x,t)] = [xi(max) — Xi(t)][xi(t) - xi(min)]
{0 if hi(x,t) 20
11011) =

IKi if hi(x,t) < O
K.

L s n is the set defining the constrained state variables

I
O

Xn+l(to) -
Xn+l(tl) =

Using again the method of Lagrange multipliers results in

11

t m
J = g I 1 {xTQx+uTRu + xT(Ax+Bu-x) - Z Y,v,(u )
t ._ 1 1 1
o 1—1
- . d o
+ xn+l(fn+l Xn+l)} t (4)

Thus, it is possible to formulate a minimization problem
with constraints, and computer techniques using the gradient
method have been deve10ped that converge to a solution [SAG-l],
[WAN-l]. It will be shown in the development of this investiga-
tion that in the case of the quadratic performance index, an
optimal solution with control and state variable constraints is
obtained without the need to solve the complex two-point boundary

value problem implied by equation (4).

2.2 Controllability and Observability.

The concepts of controllability and Observability were
first introduced by Kalman [KAL-l] who showed that they appear
as necessary and sometimes sufficient conditions for existence
of solutions in control problems, particularily those involving
multiple inputs and multiple outputs. One of the best "state of
the art" papers is that of Kreindler and Sarachik [KRE-l] and is
the source for the following definitions. Discussion will be
limited to linear time invariant systems defined by equations
(1) and (2).

Definition 1: A plant is said to be completely state-controllable
if for any tO each initial state x(to) can be transferred to
any final state x(tf) in a finite time tf 2 t.

It can be shown that a linear time-invariant plant is

completely state-controllable if and only if the nan matrix

12

2 n-l
B

E=[BIABIABI ------- IA 3
I I I I

satisfies the relation
rank E = n .

One important application of complete State-controlla-
bility is that it can be shown [WON-l] that if the system (1)
is completely controllable, then there exists a nonsingular

matrix H such that

 

 

 

 

X = O O
. 1
Lial -a2 -a3 .......... -an d
r n
0
_ O
B = .
l
. J
y = Ay + Bu .

That is, the system (1) can always be transferred into the phase

variable canonical form.

Definition 2: An unforced plant is said to be completely observable

on [to’tf] if for given t0 and tf, every state x(to) in X

can be determined from the knowledge of y(t) on [t0,tf]. If

the above is true for every t0 and some finite tf > tO then

the plant is said simply to be completely observable. If the

above is true for every t0 and every tf 2 to, the plant is said

13

to be totally observable.

As in the case of controllability, there exists a
necessary and sufficient condition for total Observability.
A linear time-invariant plant is totally observable if and only

if the anm matrix

. . T -1
F = [GTE ATCT: ------- ' (A )n CT

satisfies the relation

rank F

ll
:1

In the Specific problem to be discussed hereafter, it will
be assumed that C is the unit matrix, and, therefore (1) - (2)

is a totally observable system.

2.3 Analytic Solution for Regulator Problem.
Let x = Ax + Bu (1)

represent the state equation for a regulator, i.e., a control
u(t) is desired such that the state vector x(t) is returned
to zero.

T

T
Let J = % [x Qx + u Rujdt . (5)

0918

If all the matrices are constant, Q and R positive definite,
u(t) unconstrained, and the system completely controllable,

then it is well known [ATH-l] that a unique optimal control

* I o I
u (t) eXists and IS given by

u*(t) = -R-1BTKx(t). (6)

The constant positive definite Symmetric matrix K satisfies

14

the algebraic matrix Riccati equation

KA+ATK-KBR-LBTK+Q=O. (7)

Since equation (7) is a system of quadratic algebraic
equations, it does not, in general, have a unique solution.
Kalman [KAL-2], however, has shown that if there exists a matrix

H such that

then total Observability of the pair [A,H] aSSures the existence
of a unique positive definite solution. In addition, this condition
relaxes the requirement on Q in that it need be only positive
semidefinite.

The solution of the matrix Riccati equation may be a
formidable task since it requires the solution of a system of
n(n+1)/2 simultaneous nonlinear equations. Computational
stability and computer capacity [BLA-l] set an upper bound on
the order of the system which can be handled. A recent technique
[KLE-l] indicates that the algebraic matrix Riccati solution can
be obtained for n = 10. However, it may be more efficient
computationally to solve the matrix Riccati equation in dif-
ferential form by using established Runge-Kutta [ATH-Z] pro-
cedures. This follows from the fact that the algebraic Riccati

equation represents the steady state solution of

R(t) = -K(t)A - ATK(t) + K(t)BR-1BTK(t) - Q (8)

with

15

lim K(tO) = K.

t -~-co
0
An algorithm developed in Chapter III for a scalar con-
trol u(t) eliminates the need to solve the matrix Riccati
equation (7) or (8) to obtain the K matrix, and the computa-

tional limitation on the order of the system is relaxed.

2.4 The Inverse Problem.

It has been recognized for sometime that the Q and R
weighting matrices in the quadratic performance index are really
engineering design parameters. Since no generality is sacrificed
by aSSuming R as the unit matrix, the design problem is reduced
to finding an acceptable Q matrix. This is the so-called inverse
problem.

The general inverse problem was first posed by Kalman

[KAL-Z] in the following fashion:
t1
Let J = lim f0 L[x(t,to),u(t)]dt
tldm

then, given a completely-controllable constant linear plant and
constant control law, determine all loss functions L of the
quadratic form such that the control law minimizes J.

Kalman then proceeded to derive the necessary and suf-
ficient conditions for the solution of this problem in terms of
a frequenCy domain characterization. Let

T . . . .

Q = H H be a non-negative definite matrix

K be a symmetric positive definite matrix satisfying

the steady state Riccati equation

16

-KA - ATK = HTH - KBBTK
R = [l], i.e, a scalar control

(SI - A).1

n(S)
Y(S)

IsI - AI
6 = A - BBTK

YG(S) = IsI - GI .

Theorem (Kalman):

Consider a completely controllable plant and the associated
variational problem with a quadratic performance index such that
the pair [A,H] is totally observable. Let KB be a fixed
control law. Then a necessary and sufficient condition for KB
to be an optimal control law is that KB be a stable control

law and that the condition

I1 + BTKm(im)BI2 = l + HHm(iw)BH2 (9)

hold for all real (1).

This is an important fundamental result which establishes
a relationship between the time domain of the state variable
approach and the frequency domain of conventional feedback theory.
In fact, the term 1 + BTKm(im)B represents the return difference
in feedback theory.
Theorem (Kalman):
A control law KB is optimal if and only if
a) IsI - GI satisfies the Routh-Hurwitz conditions;
b) m?) = IwcumIz - New?

This theorem implies that it is possible to characterize

optimality in terms of closed loop and open loop characteristic

l7

equations. Unfortunately, the condition b) is not known in
general form for application to any nth order system.

These results are fundamental but are not easily applied.
Equation (9) implies that an H matrix, if it exists, can be
found if K is known. Kalman proposes a solution using spectral
factorization, but this implies the non-uniqueness of H and,
therefore, non-uniqueness of Q. Furthermore, the actual engineer-
ing design solution is the K matrix which determines the feed-
back coefficients for an optimal control system; and, thus, to
obtain H it is pre-supposed that the design information is
already known. For these reasons there appears to be little or
no evidence that these results have been applied to engineering
problems.

Three years after the publication of these results,

' procedure

Athans, et. a1. [ATH-3] were proposing a "cut and try'
based on the work of Bryson and Ho to obtain acceptable Q and

R matrices. Initial values were chosen on the basis of

-l _ 2
qii — max xii(t)
rTI = max u?.<t>
JJ JJ

that is, the diagonal entries of the Q and R matrices were
an inverse function of the maximum magnitude of the variable to
be weighted. Again, an a-priori knowledge of system behavior
is implied. If the initial "guess" was not satisfactory then
adjustments were made in a "cut and try" fashion. Of course,
each time Q or R is adjusted a new K matrix must be

determined from the matrix Riccati equation. An essentially

18

"cut and try" procedure similar to the above is still common
practice in engineering work.

Rekasius and Hsia [REK-l] attempted a more formal
approach which was motivated by the Problem of Letov [LET-1].
In this problem the control u(t) is assumed a scalar with the

constraint Iu(t)I S 1 and

(xTQx + ru2)dt .

(4
II
NY"
OL‘1 8

It was found that under the condition, a < O,

QBTK = BTKEA - (BBTK)/r] (10)
the control
u*(x) = ~(BTKx)/r ; IuI s 1
u(x) = 7'c
Sgn u (x) ; IuI 2 l

is optimal.
ASSuming a phase variable canonical form for the A

and B matrices, they defined

F 1
kl/kn

kZ/kn

-(KB)/r = kn . (11)

 

 

Substituting A, B and (11) into equation (10) results in

(n-l) equations

19

5‘
w
w

 

 

 

 

 

 

1 n-l l _
k k ' an R + a1 ” 0
n n n
k k
2 n-1 n-1 _
kn kn an kn + a2 ’ O
kn-l 2 . kn-l _ kn-2 + a = 0
kn ’ an kn kn n-l

For an optimal solution to exist, these equations must be
independent and have at least one set of real solutions such that
u(x) is a stable control law. Once the ratios ki/k = vi,

n

i = l,2,.....,n-l, are known, substitution in equation (ll)'yields

r

 

 

L. A
where krm is the bottom-row, last column element of the K
matrix, and once known the Optimal control law follows as well
as the kni’ i = l,2,...,n-l, elements of the K matrix. AS in
Kalman's work, Rekasius and Hsia assumed a kn as well as the
rest of the elements that determine the entire K matrix. They
then solved for the control law as well as the Q matrix by
substitution in the steady state Riccati equation.

Tyler and Tuteur [TYL-l] approached the inverse problem

by assuming Q to be a diagonal matrix, R the unit matrix and

20

obtaining a relationship in the frequency domain which was a
function of Q. The procedure involves trial and error adjust-
ments of the Q elements until satisfactory reSponse is achieved.
Once Q is established in this fashion, the K matrix is obtained

by solving the matrix Riccati equation.

ASSuming x Ax +-Bu

y Cx

m
J = % I (xTCTQCx +-uTu)dt
o

and using the canonical equations
r .I r- Wr‘ fi P‘ ‘N

8 A -BBT x X
0000 = 00000000000000. 000 —F ...

i -CTQ 0 -AT

I A
L J L J L. 4 L .J

.the characteristic equation of the optimal system is defined by

 

 

 

 

 

 

 

 

n
IsI - A + BBTKI = O = H (s - a ) . (12)

However, it is well known that the 2n eigenvalues of the

FC matrix consist of the n eigenvalues equivalent to the
roots in equation (12) plus n eigenvalues which are the mirror
image about the imaginary axis of the s-plane. Therefore,

IsI - FCI = .

(S - oi) (5 + (xi)

":35
H

Pre-multiplying FC yields the following

21

(SI - A)”1 o (51 - A) BBT
-(CTQC)(SI - A).1 I CTQC (ST + AT)
I (SI - A)- BBT
0 -CTQC(sI - A)-1BBT + (SI + AT)

and

IsI - FCI = IsI - AI I(sI +-AT) - CtQC(sI - A)-1BBTI = 0 .

By further manipulation
2(n-l)
I31 - AI I31 + ATI + (-1)n z kiNi(s)Ni(-s) = 0
i=1
where Ni(s) is a polynomial in 8, each ki consists either of
a qii element or a product of qii elements.

This rather formidable relationship of (2n)th degree
containing 2n terms must be satisfied by adjusting the qii
elements via root locus techniques until the desired system
response is obtained. This leads to a solution for the entire
Q matrix which is then used in the matrix Riccati equation to
obtain the K matrix. The difficulty of using root locus
techniques with a high order system is alone a challenging task.

A technique proposed by Chen and Shen [CHE-l] appears
to overcome some of the difficulties associated with the Tyler
and Tuteur method. Given a quadratic performance index and
specified closed loop eigenvalues, a sensitivity relationship

between the eigenvalues and the Q matrix is used to solve

iteratively for the desired Q matrix.

22

Letting G = A - BB K

it follows that

where *1 and ui are the distinct eigenvalues and eigenvectors,
respectively. This leads to the Jacobi sensitivity formula

[Van-l]

dAi = <vi,dcui> (13)

where V1 is the reciprocal basis of ui, i.e., [MGR-l] if

_ I I I
U [Ul:u2:ooooo.:un]

_ T '1 _ I l I
then V [U ] [VIIVZ ....1 vn] .
Since KA +'ATK - KBBTK = -CTQC
then dKA + ATdK - dKBBTK - KBBTdK = -CTdQC
T T
or dKG +'G dK + C dQC = 0 (14)
and d6 = -BBTdK (15)

By combining equations (13), (14) and (15) it can be shown that

T
dxi = -tr(uiviBB dK)

(16)
«Ire + GTdK = -CTdQC

Assuming Q is diagonal and given a set of desired eigen-
values, equation (16) and the steady state matrix Riccati equation
are used within an iterative loop to obtain a Q matrix.

It is clear that this method suffers from two disadvantages:

1) a desired set of eigenvalues is not always available to the

23

designer; and 2) an iterative scheme that involves solving the
nfltrix Riccati equation, particularily when the order of the
system is moderately large, will be either extremely time-con-
suming or computationally unstable.

In summary, although the classic inverse problem can be
solved, the engineering designer requires a much better tool to

do his job.

2.5 Sensitivity Analysis.

In the proposed procedure to obtain a solution to the
linear regulator problem with a quadratic performance index and
constrained state variables, the problem is essentially reduced

to optimizing a n-parameter system using partial derivatives of

ax,
the form -——£ . It is common practice [TOM-1], [PER-l] to define

BP-

J
this type of derivative as the sensitivity of state variable
xi to the parameter pj.
Following the usual method of development, let V be
the an matrix whose elements are
axis)

Vij(t) = SST——

Differentiating with reSpect to time

ax.(t)
to.) = —-d ———l
1] dt apj

where x. = f,(§,t).
1 1
Before further expansion it is first necessary to show that

matters are considerably simplified if the following assumptions

are made:

Then,

Since

and

Since

and

24

2
a xi(t)
1) -—-—-—- is continuous
BtBPj

2) the system is time invariant, i.e., p is not an

explicit function of time.

2 2 2
~ axi(t) _ a xi(t) a x(t) a a x(t)

 

 

 

 

 

 

 

 

 

 

:2 -— + ~ .—-
t . t . . t t
BPJ . 8 BPJ BPBPJ a BPj
. . 2
dxi(t) = a_— axi(t) +Iaxifll) SE = a Xi(t)
aP. dt ap. at 53 dt ap.at
J J J
2
a xi(t)
——————- is continuous, then
BtBP-
J
2 2
a Xi(t) a Xi(t)
t . = , t
5 SP] BPJB
Q (t) = d__ dxi(t) = axi(t)
ij , dt .
BPJ BPJ
x = Gx
ii(t) = gi(i.5.t)
= a__
V1] (t) apj [81(Xspst)]
axl<t> apj 5x2 an axn (t) av:

n

n
r2

 

Isa: ‘J ax (t) 581
z + -—-
r=1 ax (t) apj apj

 

[as: W‘] ( ) + 531
t _
_1 OX (t) r3 an

. th
In matrlx form for an n order system

88.
+ISPL
j

(17)

25

 

 

 

 

 

 

 

 

r a r' Ag 53 n r' 1
6 v v -——l-—— ‘-l—-— v v v
11 12 ... 1n ax1(t) ... an(t) 11 12 °°° 1n
. . 58 as
v v ... v -—Jl——- -Jl—-— 'v v ... v
L 111 n2 nnJ L ax1(t) an(t)J L'nl n2 nnJ
r- n
881 581
apl 0000 apn
+ . . (18)
58“ Bgn

Since the matrix C can always be transformed into the phase

variable form

g1(t) = X2(t)
g2(t) = x3(t)

gn_1(t> = xn<t>

gn(t) = -p1x1(t) - P2X2(t) - ........ - pnxn(t) (19)

substituting into equation (18) and then decomposing into

 

 

 

j = l,2,...,n first order differential equations yields
r . -n r a r ‘1 1
vljno o 1 0 ..“H..... o vljno Io
v2j(t) . 0 1 ...... .... 0 v2j(t) .
0 0 0 0 1 000000 0 0 0
. = - . . 0 . . - . (20)
O Q 1 0 0

 

 

 

 

 

or V. = GV. + Bx.(t)
J J J

26

where Vj is the jth column of the an matrix V, and

 

 

Since Vj(to) = O for normal physical systems, [ROE-1]

Gt t
i 1 -GT
= d 22
Vj(ti) € Ito e ij(T) T ( )

Equation (22) will be basic to some of the development

in Chapter III.

CHAPTER III

DEVELOPMENT OF A SOLUTION

In this chapter the new developments pertinent to the solution

of the linear regulator problem with state variable constraints

are presented. A simple algebraic algorithm is developed for a
large class of problems to obtain the gain matrix K, associated
with the matrix Riccati equation, from an n-parameter optimization
solution. Thus it is not necessary to solve the usual n(n+l)/2
simultaneous nonlinear equations. The parameter optimization
solution is based on a quadratic cost minimization algorithm

that follows from the K matrix algorithm and a special constraint
function developed for this purpose. Both of these algorithms

are described in detail.

3.1 Canonical Transformation.
Given the linear dynamic system
x=Ax+Bu (1)

where x is the nXl state vector, u is the mxl control
vector, m s n, A is a constant an matrix and B is a constant
nXm matrix, then for a physically meaningful problem an originally
unstable System can always be stabilized by a properly chosen
controller [KAL-Z]. For the optimal control designed to minimize

the quadratic performance index

27

 

28

to
J = gflxTQX + uTRu‘Idt (2)

o
where the constant an matrix Q and the constant me
matrix R are positive definite, it has been shown [ATH-l] that

the unique optimal control exists and is defined by

u* = -R-]'BTKX (3)

where K is a unique positive definite symmetric nxn matrix

satisfying

T -T
KA+AK-KBRlBK+Q=O (4)

Lemma 1: If the symmetric matrix K is unique then the an
matrix Q is unique.
Proof: Q = -KA - ATK + KBR-lBTK

Since for a given system and performance index A, B, R are

unique, it follows that Q is unique.

Q.E.D.
It is known that [KAL-Z] any completely controllable

system with a scalar control u(t)
y = Fy + Hu

can always be transformed into the phase variable canonical form

 

 

F" P' S P 1
21? 0 l O ...... .. O 21 O
22 O O l O 22
0 0 0 O 0 O 0
. = . . . . . + . u (5)
z -a -a -a ... -a z 1
n o 1 2 n-l n
L 4 k .d L A

 

 

 

 

29

Quite often the phase variable form follows directly

from the original differential equation for the system

(H) (n-l) (1) (0)

x +an_1x +.........+a1x + aox =bu (6)

By suitable scaling of the state vector 2, i.e.,
z. = x,/b
1 1

the form in (5) can be transformed into the more convenient form

")c ‘ F0 1 0 ... 0 x ‘ r()1

 

 

 

 

 

 

l 1
x2 . O 1 x2 .
. O . . .
= . . . + . u (7)
. l .
I.XnJ L-ao -a1 ...... "an- an LI)J

which will be used throughout the sequel. In matrix form (7)

becomes

x = Ax + Bu

3.2 Gain Matrix.

One of the major objectives of this research was to
eliminate the need for the usual iterative procedure [ATH-2],
IKLE-l] to solve the n(n+l)/2 equations represented by the
algebraic matrix Riccati equation (4) in order to obtain the gain
matrix K. It will be shown that this can be done for the system
(7), assuming a quadratic performance index structure.

First, it will be assumed that R is the unit matrix.

This is commonly done and does not affect the generality of the

30

results. In effect, the matrix Q is redefined to be the matrix
Q1 such that it "absorbs" the difference between any other R

and R = I. From (4)

T -1 T
Q = -KA - A K + KBR B K

and it follows [KAL-Z] that, for a completely controllable system

and Q positive definite,

J(t ) = a
O

FTC—‘38

IXTQx + uTRujdt = % XT<to)KX(to) (8)
0

where K is a unique solution. Since A, B and K are unique,
by Lemma 1 it is possible to generate a unique matrix Q1 with

R = I Such that

Q1 = -KA - ATK + KBBTK (9)

J(t0) = ,1, I [XTle + uTu'Idt = 1g XT(tO)KX(tO) (8a)
0

Since K is the same in both cases, the quadratic performance
index in (8a) is equivalent to the one in (8).

The second aSSumption is that the Q matrix be diagonal.
In practice this represents the largest class of problems Since
it is quite often very difficult to ascribe weighting coefficients
to state variable cross products.
K Matrix Algorithm:
Given the completely controllable system (7), R the unit matrix
and Q a diagonal (with unknown entries) matrix, then the complete
symmetric K matrix may be obtained from the n entries in the

th .
n row of the K matrix.

31

Proof is by expansion of

T T
Q = -KA - A K + KBB K

Taking advantage of the fact that kij = kji’ it can be readily

established that

qu = Zackn’l + bzki’l (10)
_ 2 2
qii — -Zki-l,i + 281-1 n,1 + b kn,i , 1 = 2,3, .,n (11)
2
qij = -kl,j-l + aokn,j + aj-lkn,l + b kn,1kn,j
j = 2,3,.....,n (12)
q1j - -ki,j-1 i-l,j i-lkn,j j-l n,1 b2kn,1 n,j
1 = 2,3, ..... ,n
j = 1,2, ..... ,n (13)

Since Q is diagonal, all qij = 0 when 1 ¢ j. Then,

beginning with equation (12)

k a k k

= + a k + b
l,j-1 o n,j j-l n,1 kn,ln,j

j = 2,3, ..... ,n (14)

showing that all first row entries in the K matrix can be

obtained if the k , entries are known, 1 = 1,2, ..... ,n.
n,1

Continuing with equation (13), Since all the ki 1 j
' s
are known from the previous row computation, a relation for
. . is obtained
13.]-

= -k + a k + a +
i,j-l i-l,j i-l n,j j-lkn,i b kn,1 n,j

i = 2,3,....,n

(15)

L_.I
II
H
v
N
u
o
o
u
:3

32

The originally unknown qii elements are then computed
directly from the knowledge of the K matrix obtained via equa-
tions (10) and (11).

For example, with n = 3

_ 2
k11 _ aok32 + a1k31 + b k311‘32

2
R12 ‘ a01‘33 + a2k31 + b k311‘33

_ 2
k22 ‘ 'k31 + a11‘33 I a2k32 + b k321‘33

, th
Together W1th the three elements of the n row of the K matrix,

these define the complete matrix since K is symmetric. And

2
C111 2‘301‘31 + b k31
= -2k + 2a k + bzkz
q22 12 1 32 32
— 2k + 2 k + bzkz
q33 " ' 32 a2 33 33

It is obvious that if kn i’ i = l,2,...,n, are available,
3

the entire K matrix is readily determined with Simple algebra;
and solution of the n(n+l)/2 simultaneous nonlinear equations

becomes unnecessary. The method of obtaining the kn 1 forms
,

the bulk of this research effort.

3.3 Autonomous Form.
A -
Since u = -R 1BTKX

for R = l, i.e., a scalar control,

3%
U — "b(knlxl + knzxz + 0.000. + knnxn)

*
Thus, it is clear that the optimal control u is a function

33

t
of only the n h row K matrix elements. Substituting in (1)

 

 

 

 

x = (A - BBTK)x = Gx (16)
where
I. 1
O l 0 O
O l
O .
G= . . . . (17)
l
2 2 2
- +b - +b ......... - +b
(a0 knl) (a1 kn2) (an-1 knn)
L- d
L t r a + bzk
8 pi i-l n1
rho 1 0 0 0 -I
. O 1 0
. O l
G = . . O (18)
0 0 1
-p 'P -p oooooooooooooo "p
L 1 2 3 “J
The problem can now be posed as, given the autonomous
system

x = Gx
find the parameters p, so that the quadratic performance index
1
cost

_ T
J(to) - t X (tO)KX(tO)

is minimized. In this sense the problem becomes one of parameter

optimization; and since the parameters are simply and linearly

34

related to the kni elements, once the optimum parameters pi
have been found, the entire K matrix is known by use of the K
matrix algorithm.

Constraints on State variables often exist and preclude
a simple choice of K which minimizes the cost J. The follow-
ing theorem is new proof that the cost J is a Lyapunov function
and leads to a general cost minimization criterion which can be
used to minimize the cost in any subinterval of [to,m].
Lemma 2: Given the completely controllable linear time invariant

system
x=Ax+Bu
and the cost functional

J = 15 IEXTQX 'l' UTRU:Idt
0

where u is not constrained and K is a symmetric positive
definite matrix satisfying

T -l T
-KA - A K + KBR B K - Q = O
Q positive semi-definite and R positive definite, then the

matrix KB(R-1)TBTK is positive definite.

Proof:
uTRu is positive definite by hypothesis
u* = -R-1BTKX is the optimal control
Substituting,
u*TRu* = -xTKB(R-1)TR - R—lBTKX

T

xT KB(R-1)TB K x

35

-l T T
which implies that the matrix KB(R ) B K is positive definite.
Q.E.D.

Corollary: If R = l, and K is the an matrix and

 

 

then the matrix

r")
B
L..l

II

[R k.]

ij in nj
is positive definite.
Proof is by direct substitution of R and B into the result
of Lemma 2.
Theorem I:

Given the completely controllable linear time invariant system

Ax + Bu

)4.
II

“1 T
-R B KX

and u

Q and R positive definite and K the Symmetric positive
definite matrix satisfying the steady state Riccati equation

-1T
KA+ATK-KBR BK+Q=O,

then xTKx defines a Lyapunov function.
Proof:
T . . . . . .
x Kx IS pos1t1ve def1n1te by hypothes1s

xTKx = 0 if and only if x = O

x = Ax - BR-LBTKX = A - BR-LBTKX x = Gx

36

T T T T
—E XTKX = XTKX +-x KX = (GX) KX + X KGX = X B K + KG X

Substituting in the Riccati equation, and since

A = G + BR-lBTK

K G + BR-IBTK + c + BR—lBTK TK - KBR-IBTK + Q = 0
KC + GTK + KB(R-1)TBTK + Q = 0
then KG + GTK = -Q - KB(R-1)TBTK

Since Q is positive definite by hypothesis and KB(R-1)TBTK

is positive definite by Lemma 2, then KG + GTK is negative
T T T

definite. Since gg'x Kx = x KG + G K x then it follows that
d

EE'XTKX is negative definite. Since the System is asymptotically
stable by hypothesis, then xTKx is a Lyapunov function.
Q.E.D.
It follows that given the constant matrix K and t1,

t s t s m, the cost J1 for the time interval [t ,t1] is
0

J1 = %[xT(tO)KX(tO) - xT(tl)Kx(t1)]

3.4 Minimum Cost Algorithm.

The objective is to find the minimum value of the performance
index J Subject to the following restrictions:
1) the system x = Gx is asymptotically stable
2) the weighting matrices Q and R are positive
definite
3) the gain matrix K is positive definite

4) the constraints on the state variables are satisfied.

37

It is convenient to let t = m, then

1
J = ‘5 Jo me)
O O

For convenience of notation the argument of the state vector

will be dropped,so that
J=1§XKX (19)
It can be readily shown that by expanding (19)

2 2

2
= L + 0000 + 0.0 +
J ZIkllxl + k22x2 knnxn + 2(k + + k x x

12X1X2 ln 1 n

k23x2x3 + ....... + kn-1,an-lxn)]

where the kij are the elements of the K matrix.
In more compact form
n n
J = g 121 jil kijxixj (20)

It has been shown earlier that given the completely con-

trollable system
X = Ax + Bu

u a scalar control, J a quadratic cost functional with Q
diagonal and positive definite and R = [1], then the positive
definite gain matrix which defines the Optimal control can be
determined from equations (14) and (15). That is, all elements
of the K matrix are uniquely determined once the bottom row

elements are known. Therefore, J can be found if the vector

38

r kn1q
n2

KN = . (21)

 

 

g nn J
is known.
. 2 .
Since p, = a, + b k , ; 1 = 1,2,.....,n
1 1-1 n1
a region in p-space of the system x = Gx can be readily
established which assures asymptotic stability by the use of

the Routh Criterion, i.e., forming the array from the char-

acteristic equation

n n-l
+ + ....... + =
1 put + p21 pl 0
1 pn_1 pn_3 .... .....
pn pn_2 pn_4 .........

 

Since a zero entry in the first column is not permitted,
there will be, in general, n constraints that must be satisfied
such that each element in the first column, starting with the
second row, must be positive.

The Routh Criterion is used to set bounds on the minimum

values of k ,, i = 1,2,....,n. For example, assuming a third
n1

order system

39

3 2
A +P3A +P2A+P1=O

 

1 p2
p3 p1
p p -p
2 3 1 0
p3
p1
the - a + bzk > O
“ p3 2 33
2
k33 > —aZ/b
= a - bzk > 0
pl 0 31
2
k31 > ‘ao/b
p2P3 > p1
+ bzk )(a + bzk + bzk
(a1 32 2 33) > (a0 31)
2
a + b k
k > 1/b2I- O 31 - a
32 + bzk 1
I82 33

The strict inequalities define open sets so that one must in
practice adjust the relationships by choosing €i > O, arbitrarily
close to zero, such that

, _ 2
m1n k31 aO/b + 31
2

 

+ b k a
2 a
min k = 1/b ° 31 - —l-+ e
32 + bzk b2 2
32 32
' k = a /b2 +
ml“ 33 2 63

This is defined as the e-minimal constraint set and may lead to
an e-optimal solution, i.e., there exists a neighborhood, however

small, so that if 5(a) is the e-Optimal parameter solution,

40

then there exists a 6 = HE - 5(a)“ > 0 such that

J(5) < J(5(€))-

Quite often the designer must obtain the minimum cost
optimal solution within a given range of eigenvalues as determined
by the desired reSponse of the closed loop system. In this case,
the maximum and minimum values of the kni can be found by the
following procedure. Given the eigenvalues *1 for an

asymptotically stable system then [TUR-l]

pH = -()\1 + A2 + ...... + In)

= +... + ... ...
Pn_1 (AIAZ Alhn A2A3 + + AZAn +

+ I‘n-lxn)

_ n
p1 - (-1) (ilkz ...... In)

2

and k . = (p. - )/b

a
n1 1 i-l
Criteria for the positive definiteness of the Q matrix
are readily established in terms of the KN vector, since
q..>0g ; i=1,2,.....,n

q. =0 ; i,j=1,2,.....,n

then equations (10), (ll), (12) and (13) yield a set of equations

(22) which is defined as the Qii Criterion.

41

2 2

= 2a + k + b >
q11 0 n1 kn1 O

‘ 2k + 2 k + bzk2 > 0
C122 ‘ 12 a1 n2 n2

2 2
= -2 + 2 + b ‘ O

qn-l,n-l n-2,n-l an-2kn,n-l kn,n-l /

_ 2 2
q — -2k + 2a k + b k > O (22)

nn n-l,n n-l n,n n,n

The Qii Criterion has the effect of a filter which
passes only the K matrices that satisfy the constraint that
the Q matrix be positive definite. However, it is used in
the Minimum Cost Algorithm as a means of generating a succession
of KN vectors in a systematic fashion so that the arbitrarily
small neighborhood of the absolute minimum of J is located.
Computationally, it is convenient to start with the last equa-
tion in (22), i.e., start with an initial kn n and find a

3

satisfactory kn Having obtained k the relationship

. ’
,n-l n,n-l

f is used to solve for a satisfactor k .
or qn-1,n-1 y n-2,n-l
But k is in general a function of the vector K , there-
n-2,n-l N
fore, kn n-2 is adjusted until qfi_1’n_1 > O is satisfied.
The process is continued with qn 2 n 2, etc. This procedure is
- , -

flow-charted in Figure 3-1.

3.5 Search Procedure.

The solution for an absolute minimum or e-minimum of J
is achieved by the adjustments on the KN vector, as shown in

Figure 3-1, within the region defined by the restrictions and

 

42

 

[initial KN I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Adjust
k
nn
. 0
Adjust kn2
' i
0 ..
‘L VINO
‘ Yes
NO Adjust
9 I a I? ._—‘ I
q22 > 0 . kn1 kn1 < m1n kn1
’Yes '
Adjust k i

 

n
i Routh I NO I

Criterion I

IYes

Compute Pl; Print
J I I

Adjust

 

 

 

 

 

I!

 

 

 

DD

Figure 3-1

 

Minimum Cost Search Procedure

43

constraints.

One iterative scheme would be a gradient search using

r . r . a
kn1(1+1)I kn1(1)
kn2(i+1) kn2(i)
. = - 471- (23>
HV JH

 

 

 

 

Lknn(i+1)J L knn(i)J

where the constant B determines the step size and

hi__
5knl

Bi...
9kn2

(24)

VJ

ai__
Laknn

 

 

II
DJ

For example, assuming n

J = 1/2[k x2 +-k x2 - k X2 + 2(k )<)( 4-k. x x +~k x x )1

11 l 22 2 33 3 11 1 2 13 l 3 23 2 3

Using equations (14) and (15)

- 1/2[ k + k + bzk k 2 + k + k + k +

J ' (a0 32 a1 31 31 32>x1 (‘ 31 a1 33 a2 32

2 2 2
b + + k + k +

k32k33)X2 + k33x3 I 2k31X1x3 2k23X2X3 2("Io 33 a2 31

2
b k31k33)"1x2]

53 _ 2 2 _ 2 2
Ok31 - 1/2I_(a1 + b k32)x1 x2] +-x1x3 + (a2 + b k33)x1x2

x(to)

unique vector KN.

 

44

It is useful to know that there are initial conditions

such that if the gradient

VJ = O

(20)
2
VJ = b MKN - F
2 2
r- 2 -
x1 2x1X2 (X1X3 x2) .....
2
0 X2 2X2x3 0000000
2 2
2x1x2 x2 C) x3 .........
I
b 0000000000 00 00000000000000 0 00
that is, M is a symmetric an

the diagonal.

 

Ffleia) 1

2

f (ai,x)

~

 

exists,

it defines a

Using equations (10) through (15) in equation

(25)
2("1"n-1'X2Xn-2) I
2(XZXn--1.x3Xn-2)
2(X3xn-l-X4Xn-2)

2(

Xn-3Xn-1-Xn-2Xn-2)

 

matrix with zero entries on

45

If VJ = 0, then equation (25) becomes

and if M is nonsingular then the vector K.N can be uniquely
obtained for a given system from the initial conditions. Under
these conditions if KN is not in the region of interest then
there exist no extrema in the region of search and the minimum
of J should be on the boundary of the search region.
Unfortunately, if KN lies within the region of interest
nothing can be said in general about this extremal, if it exists,
since sufficiency conditions for a minimum or a maximum cannot

be satisfied. To Show this

 

 

 

6% a—J—R— oooooooooo a—‘l—T
aknl 9kn1a n2 6knla nn
2
a_J
2 2
ii—___ .....OOOOOOOOOOOQOOOOO. LL
2
O nnaknl 8k 3
‘ nn
Since the m11 element of the M matrix is always zero, it is

clear that M can never be positive or negative definite; and,
therefore, the sufficiency condition for a maximum or a minimum
cannot be satisfied. In fact, it is possible for KN to be a

saddle point.

46

One disadvantage with the gradient method is that the
computer program becomes quite involved for the general case
which includes the possibility of M being Singular. On the
other hand, because of the Simplicity of the function to be
minimized and because the region of search is easily defined,
it has been possible to prepare a simple digital computer program
based on the exhaustive search technique that is extremely fast
and very satisfactory. For example, the neighborhood of the
absolute minimum of a third order system was obtained with
sufficient accuracy in about one minute on an IBM 1800 computer,
most of the time being consumed by the printer. Adequate resolu-
tion was obtained after four successive contractions of the search
region.

The procedure is summarized as follows:

1. Determine the maximum and minimum kni’ i = l,2,...,n, and
use as input to the computer program;

2. Choose 01’ 0 < 01 < l, to increment k i’ i.e.,
n

Ak . = a.[max k , - min k ,1
n1 1 n1 n1

It has been found that a. = 0.1 gives satisfactory reSults and
1

has been incorporated in the Computer program;

max k

3. Starting with the n-tuple [max kn n2""°’ max k

l’ nn:I

use the scheme in Figure 3-1 until the whole region is scanned;
4. Determine the minimum J obtained from Step 3 and form a
neighborhood about the correSponding KN in KN-space that in-
cludes the absolute minimum of J. This contracted region now
defines new maxima and minima for the kni;

5. Repeat Steps 3 and 4 until the desired accuracy in locating

47

the minimum of J is achieved;

6. If the value obtained in Step 5 is such that all state
variable constraints are satisfied, then this is the desired
solution to the minimization problem.

In general, the above procedure does not yield a solution
where the state variable constraints are satisfied. To implement
the search for this solution a special constraint function has been
developed which forms the basis of a hybrid computer program that

complements the Minimum Cost Algorithm program.

3.6 Constraint Function.

Any type of constraint function is acceptable providing
it can be formulated in mathematical terms and at the same time
reflects a measure of "goodness" with reSpect to a physically
meaningful criterion. For example, Kalman [KAL-3] discusses
the Lyapunov function of a stable system as useful in evaluating
transient reSponse. Thus, given a constraint in terms of
"Lyapunov distance" from the origin in State space as a function
of time, an upper bound can be found and invoked as criterion
of goodness.

If V(x,t) is some Lyapunov function, then

_ V(x,t)

V(x,t) — V(x,t)

V(x,t) s cV(x,t)

where c is a Suitably chosen constant. Then by a well-known

Lemma [BEL-l]

V(x,t) g epr-c(t - t0)] V(x0,to)

It has been shown that the matrices P and H exist for a

48

time invariant system such that P and H are symmetric and
positive definite and

T
x PX

V(x)

T
-x HX

V(X)
Kalman continues, showing that if xmin is the minimum eigen-

value of the matrix HP-1, then

C = xmin
Thus, the Lyapunov distance is bounded by a function of the
longest time constant. Such a bound is not new in the field
of differential equations ICES-1].
It is well known that the characteristic equation for

the G matrix in equation (16) is

n n
..... .00 + = O
i + pnx + + p21 P1

As shown previously,

= - + + ....... + =
pn (11 12 In) 0
= + 000 + 000000
I’n-1 (I1A2 HA3 + + >‘1)‘n + "2I3 + kn-lxn)
= (~1)“( 1 1 >
pl )(1 2 ..... n

and, therefore, there exists a direct relationship between the
eigenvalues of G and the parameters pi. However, it is quite
possible for a reasonable change in pi to have a relatively small
effect on the minimum eigenvalue; and there could be regions in

p-5pace such that the criterion based on xmin may be too

49

conservative, i.e., relatively insensitive.

In contrast to the minimum eigenvalue type of constraint,
a constraint based on maximum allowable State variable excursion
is usually a well-known value. The designer is well aware of the
limitations such as maximum acceleration, stress, etc. It is

this type of constraint that has been chosen:
xi(max) 2 Xi(t) 2 xi(min); all t (26)

The constraint region defined by equation (26) is still
not in a completely satisfactory form since it only provides
upper and lower bounds on the constrained state variables. While
it reduces the region of search, it still represents an uncountable
number of solutions. The search for the best solution within
the constraint region is made possible by changing equation (26)
from an inequality constraint to an equality constraint. To
implement this need, a constraint function C is defined which
serves as a test function that senses the deviation and the
direction of the deviation from the constraint bounds. In its
structure the function is similar to the penalty functions used
in adjoining to the performance index J in conventional Optimal
control design involving constraints.

Define
hi = [xi(max) - xi(ti)]Ixi(ti) - xi(min)]
(27)

I:i a Ixi(ti)I 2 xi(t) ; ti,t 6 [1:03]

Then it is clear that

50
h, 2 0 if x (max) 2 x,(t,) 2 x,(min)
1 1 1 1 1
h,<0 if x,(t,)>x_(max) or x,(t,)<x,(min)
1 1 1 1 1 1 1

Let H be a function of hi SUCh that

H(h.) = (28)
1

if at least one state variable exceeds its bounds

Let H be function of hi such that

H(hi)=1;i€LSn (29)

if all the state variables lie at or inside their reSpective
constraints and L defines a subspace of En which contains
the constrained state variables.
A constraint function G is formed such that
C = X [x,(max) - x.(t.)][x.(t.) - x.(min)] H(h.) (30)
16L 1 1 1 1 1 1 1
which is used to test the solution obtained with the Minimum
Cost Algorithm for constraint violation; and if the constraints
are violated, it is used to obtain a new estimate for the KN

vector such that all constraints are satisfied.
3.6.1 Geometric Interpretation.

Assume that x1 and x, are two constrained state
J

variables. Expanding equation (27),

:3"
II

-Xi(ti) + [xi(max) + xi(min)]xi(ti) - xi(max)xi(min)

2
hj -Xj(tj) + [Xj(max) + xj(min)]xj(tj) - xj(max)xj(min)

51

Each of the h-functions represents a parabola as shown in Figure

3-2 with its maximum value occurring at

xi(max) +-xi(min)

xiep = 2 (31)

 

Since the C-Function represents the sum of two parabolas
in x. X x. X R space, it can be readily visualized as shown
in Figure 3-3. The shaded region represents the surface of
allowable solutions, C 2 0. This surface meets the plane C = O
at exactly four points a,b,c and d which define the four minima.
Below the plane C = O, the surface C < 0 is a paraboloid trun-

cated at the top by the circle a,b,c,d in the plane C = O.

hi‘I

hi<max> - — —

 

+xi(ti)

 

[X 1 (min) \X i (max)

Figure 3-2

h, Function
1

 

52

 

 

 

Figure 3-3

C-Function in Two-Space

> Xi(ti)

53

Theorem II:

Given the constraint function

G = iELEXimaX) ' Xi(ti.)][xi(ti) ' Xi(min)]H(hi)

xi(max) 2 xi(ti) 2 xi(min)

where L is the number of constrained state variables, L < n,

 

dQ—-——— and a C , i,j E L, exist and are continuous,
9Xi(ti) axj(tj>axi(ti)
then C has at most 2% minima, all of which occur at C = 0.

Proof:

Clearly, if xi(max) 2 Xi(ti) 2 xi(min) then

hi = [xi(max) - Xi(ti)][xi(ti) - xi(min)] 2 0

but h1 = 0 if and only if either xi(ti) = xi(max) or
Xi(ti) = xi(m1n). If Xi(ti) > xi(max) or xi(ti) < xi(m1n)
then hi < O and these Xi(ti) are not allowable.
Since all H(hi) = l, i E L,

C = 2 h.

iCL 1

then C = 0 if either each h1 = 0 or some hi < 0-
But hi < 0 implies xi(max) 2 Xi(ti) 2 xi(min) does not hold.
Therefore, C = 0 if and only if each h1 = 0.

Then the number of ways that C can be zero is the number of

ways that 2 bi = 0 holds.
iEL
Since each hi can be zero in two ways, then 2 hi can be zero
iEL
in at most 2 ways.

Since C < O is not allowable, then the zeros of C correSpond

to the minima of C when all X'(ti)’ i E L, attain one of their
1

respective bounds. Therefore, C can have at most 2% minima

when xi(ti) — xi(max) or Xi(ti) xi(m1n).

To show no other minima of C exist inside the constraint region:

the necessary condition for an extremum is that

 

ag—TE—T = O: i E If,
5X1 1

A ah,

Q— : ______ = -2 + + '

ax.(t.) AX.(C.) xi(ti) in(max) xi(m1n)]

1 1 1 1
X (max) + X (min)

. , 1 -
implies Xi(ti) = 2 at extremum, 1 e., each

x,(t ) corresponding to an extremum is unique since x,(max)
1 1 1

and x (min) are unique. These extrema then correspond to a
1

 

 

 

 

unique maximum for C since
2
ag— = -2 < 0
axi(tg)
32C 52C
= = 0, i #1
axi(ti)oxj(tj) 5xj(tj)axi(ti)
i.e., the {XL matrix
r- 1 r-
2 2
B—9-— ........ 57C 2
2
x (ti) Xi(ti)xj(tj)
2 0
0 = - 2
. O 2
2 2
a C a C
><.(t.)X.(t.) ' 2
j j 1 1 x-(t.)
J J
L

 

 

 

 

 

55

is negative definite. Therefore, all the minima of C occur

at C = 0 when constraints are invoked.

3.6.2 C-Function Iterative Procedure.

In the event that the parameter solution Obtained by the
Minimum Cost Algorithm is such that at least one of the state
variable constraints is violated, then it is necessary to adjust
the 5 vector; and the C-Function provides a satisfactory means.
Given C < O, the objective is to develop a convergent iterative
procedure so that successive estimates of the parameter vector 5
yield a solution for C = 0. Clearly, this solution yields a cost
Jc which is an upper bound for the desired solution.

Since the iteration is on p, it is necessary to transform
the C-Function into p-space. The gradient of C in x-Space is

defined by the LX]. vector

VG: . aisje'f/ (32)

_ac__
ij(tj)

L J

 

 

while in p-space it is defined as the nxl vector

56

r2211
891

v c = . (33)

a9.
L apn J

 

 

Slnce Xi(ti) = fi[p1,p2,.....,pn,xi(to)] 1f the part1al

 

 

 

 

derivatives '—S————— exist, then
i
V' a
8P1 Cpl 3Xi(ti)
I
V C =
P
2.1.1:). 2*: “1) 50
t
Lapn apn ‘4 L ij( J.)_J
= SVXC (34)

In the computation Scheme it was assumed that these
partial derivatives exist and successive estimates of B were
obtained using equation (34). However, in the event that these
derivatives do not exist -- and this can happen, it will be

shown -- special measures have been incorporated in the program

to allow the iteration to continue to a solution. A more detailed

discussion is presented in the following section of this chapter.

It is common practice to define the partial derivative
Axi(ti)

8P as the sensitivity of state variable xi
m

to the parameter

57

pm at t = ti' Thus, the S matrix in equation (34) may be

defined as the sensitivity matrix of the system. Following the

development on sensitivity analysis in Chapter II,

I" 9X1(ti) 7
39m

ax2(ti)
5P

m
Gt, t.

= -e 1I l e-GTBXm(‘T)d‘T (3S)

3x3(ti)

 

 

L Opm .4

In the hybrid computer implementation of the parameter
search, equation (35) is solved repeatedly on the analog computer
at the various times ti’ tj, etc. until the entire S matrix is
determined.

An iterative algorithm for the C-Function search must in
the least indicate the direction in which the p vector must be
changed for a closer estimate to the solution. The gradient
VPC obviously has this characteristic. This leads to the steepest

ascent or descent scheme [PER-l], [REL-l], [BEK-l]

V C(i)
'13(i+1) = 13(1) + 01 . (maximization) (36)
‘R—Invpcun
V C(i)
§(i+l) = 3(1) - a “V C(i) (minimization) (37)
p

where NV C(i)” is the Euclidean norm of the gradient vector
p .

58

at the ith iteration and a is a constant that determines the step
size and is usually obtained empirically. The determination of

q is usually a part of a "cut and try" loop within the iterative
procedure since it depends on the dynamics of the particular
system, i.e., the rate at which the gradient may be changing.

It is obvious that an iterative procedure that not only
senses the direction but also provides an estimate of the required
step size is highly desirable. This indeed is the basic appeal
of the Newton-Raphson [PER-l], [BEK-l] method.

C(i)v C(i)
fi<i+1) = 3(1) - P (38)
vacumz

 

It has further advantage of quadratic convergence. However,
care must be exercised in case vpC(i) a 0.

Both the gradient and a modified version of the Newton-
Raphson methods are used in the C-Function algorithm to take
advantage of both procedures. The Newton-Raphson method in

equation (38) is modified to include a Step size constant v

C(i)vpC(i)

vac<i>u2

 

5(i+1) = 5(1) - Y (39)

which has been found to be very helpful in
1) stabilizing any effects due to discontinuities induced by ti;
2) overcoming the inherent accuracy limitations, particularily
near C = 0, of an analog computer when measuring the S
matrix elements.
In summary, the entire computational procedure is outlined

in Figure 3-4. The initial region defined by max KN and min KN

59

is searched for the global quadratic cost minimum, JG, by the
Minimum Cost Algorithm. The K matrix corresponding to an
arbitrarily close estimate of JG is checked for state variable
constraint violations. If all constraints are satisfied, then
this is the desired solution and the search is ended. If one

or more state variable constraints are violated then the C-Function
Algorithm is used to Obtain the C = 0 solution, thus defining

a K matrix such that all constraints are satisfied and a
correSponding quadratic cost JC which is clearly an upper

bound for the desired solution. Using the same search region,
increment KN in the fashion described in Section 3.5, comparing

the cost J(i) at each point with J If J(i) is greater

C'
than JC then continue search; but if J(i) is less than J ,
check state variable constraints. If constraints are violated
continue search; but if constraints are satisfied then J(i)
becomes the new (lower) upper bound. After the entire region
is scanned, the entire process may be repeated for a closer
estimate to the desired solution by choosing the smaller region
obtained from the results Of the previous search.

The detailed computational description of the Minimum

Cost and C-Function Algorithms is covered in Chapter IV.

3.6.3 Discontinuity in p-Space.

The parameter search is based on measurements made at
ti’ 1 E L, tO a ti S T, when xi(t) attains its maximum
magnitude as shown in Figure 3-5. If Xi(ti) exceeds the con-

straint an adjustment is made in B. It is not uncommon that in

60

adjusting B to decrease the amplitude of peak 1, peak 2 may

become the dominant one as shown in Figure 3—6, resulting in a

sudden jump in the value of ti and a discontinuity in V C-

Safeguards are incorporated in the computer program so that
in this case:

1) possible oscillation of the Solution between peak 1 and peak
2 is quickly ”damped” out by the step size constant y in
equation (39). The programming details will be described in
Chapter IV, and an example of an actual case will be shown in
Chapter V.

2) the peak closest to its constraint is chosen automatically

as the solution after a finite number of iterations.

61

 

Input

max KN
min KN

 

 

 

 

 

 

 

Minimum Cost]

 

 

 

 

 

 

Algorithm
1, i A
Constraints
Satisfied Yes D I j 2 No Solution

 

 

 

 

 

 

 

Yes

L

 

 

‘3

L4

II

{—4
A

H.
v

b...

 

 

1 No
C-Function
1 Algorithm,

J
c

‘T;

Increment

 

 

 

 

 

 

 

 

V c_.

 

 

 

 

 

 

 

 

 

 

 

Figure 3-4
Algorithm to Obtain the Minimum

Cost K Matrix

X.(t)
1

xi(max)

x,(min)
1

Kim

xi(max)

Xi (min)

62

 

 

 

a...\\\\
t \
1
T >
i
’- --(-- -——---—-——-—- —~rv--—-
’ /
1
Figure 3-5
Peak 1 Dominant
2
p. _. __ .. __ __ __
?’

 

 

Figure 3-6

Peak 2 Dominant

CHAPTER IV

COMPUTER IMPLEMENTATION

The application of the concepts presented in the pre-
vious chapter generated a digital computer program for the
Minimum Cost Algorithm and an hybrid program for the C-Function
Algorithm. The computational details as well as the various

procedures are described in this chapter.

4.1 Minimum Cost Algorithm Program.

This digital computer program locates the arbitrarily

small neighborhood of the absolute minimum of

J = 1/2 xT(tO)KX(tO)

by adjusting the KN vector within the constraints imposed by

the Routh and the Qii Criteria. Basically, an exhaustive search

procedure is used, but the program has been written so that a

gradient search can be readily implemented if desired. As men-

tioned in Chapter III, the exhaustive search approach is practic-

able because of the simplicity of the function to be minimized

and because the region of search can usually be identified readily.
In reSpect to region identification, the initial search

region is the hypercube in n-Space, each ”side" of which is equal

to

[max kni - min kni] ; i = 1,2,.....,n

63

64

This defines a vector

max KN = [max kn , max k ,........., max knn]

1 n2

which is chosen such that if the vector

corresponds to the absolute minimum of J, then

7‘:
k . Smax k . ; i = 1,2,......,n
1'11 n1

A starting value for max KN can be readily obtained by either
invoking the relationships between the eigenvalues and the kni
in the case a range of acceptable eigenvalues is available to
the designer, or using the C-Function Algorithm and a "large
enough” KN guess.

The method of choosing a min KN vector

min KN = [min knl’ min kn ,......, min k n]
n

2
was discussed in Chapter III.

One of the advantages to the exhaustive search procedure
is that successive contractions of the search region can be
readily identified and easily implemented within the computer
program.

In the interest of clarity, a computer program for a third
order system will be described. Extension to any nth order system
follows directly and presents no computational difficulties. The
flow chart for this program is shown in Figures 4-l(a), (b), (c),

(d), (e) and (f), and the listing is contained in Appendix A.

 

 

IInput

Ak Ak32, Ak3 3|

 

“ i

max k33

1
l

 

“Law I———o
ii; 59

 

 

 

 

 

 

(Routh) J < 3

 

K Matrix I

 

 

 

J ’1 Storel
1 i No
V'J qll > 0 '3] ’1 Stop

 

 

 

 

 

 

Figure 4-1(a)

Minimum Cost Program

66

 

 

No

Yes

 

 

 

 

 

 

 

(k32

 

- min k32) S O ? I

 

Yes

 

 

33 33

 

 

 

 

 

I(k33 - min k33) s o ?

 

No

 

Figure 4-1(b)

Minimum Cost Program (con't.)

Yes
‘ Stop

67

 

 

 

an22>ogl Yes
[

 

 

 

 

 

 

 

 

 

 

k - Ak
k12 31 31
No
L (k31- min k31) s O ?

 

 

 

W168
”1

No
- ' 9
@- [(k32 m1n R32) S 0 .

Yes

 

 

 

 

 

Figure 4-1(c)

Minimum Cost Program (con't.)

68

V

Routh I No
Criterion l
9

 

 

 

 

 

 

 

Yes
i I Yes
J(i+1)< J(i) ?[ a®
No

 

 

 

 

 

 

l

W 3 HVJH

 

 

 

 

 

 

k31 -

1% 3
(E 1

Figure 4-1(d)

 

 

Minimum Cost Program (con't.)

 

69
Yes

1 Wk - min k31) S O ?

[k31 = min k3ll +IH

K Matrix

ML

 

 

 

 

 

 

 

 

 

1

Yes
@4 J(i+l) < J(i)
?

 

 

 

 

 

Figure 4-1(e)

Minimum Cost Program (con't.)

7O

 

 

 

 

 

 

 

 

 

 

 

 

Yes

 

 

 

 

 

No
- ' 9
33 min k33) s O . ————+®

Yes

 

 

 

 

StOp I

 

Figure 4-1(f)

Minimum Cost Program (con't.)

71

4.1.1 Incrementation.

The exhaustive search procedure requires incrementing
each element k i of the KN vector. In this program incrementa-
n

tion proceeds on the basis of
Akni = aiEmax kni - m1n kni]

Although ai may be any value Such that 0 < 01 < l, satisfactory
results are obtained by setting ai = 0.1. Since the Qii Criterion
defines an open set, not all the values in the hypercube defined

by max KN and min KN will be attained; but it is of interest

to note that a lower bound is automatically obtained from
k . > k . .. - . max k , - min k .
n1 n1(Q11) 0/1[ n1 n1]

where k '«2'i) is the minimum value of k , which satisfies
n1 1 n1

the Qii and Routh Criteria.

4.1.2 Computation.

The computation is initiated by starting with the n-tuple
[max knl’ max kn2,......, max knn]. A slight modification has
been incorporated in this program which expedites the search:

the search is begun with the (n-1)-tuple [max k32, max k33]

 

and an initial value‘for max k31 based on the Routh Criterion
(a + bzk )(a + b2k ) - a
_ 1 32 2 33 0
max k - - 0.0001
31 2
b
where = 0.0001 is used to assure a strict inequality in the

61

Routh Criterion.

72

An initial K matrix is computed using the algorithm
described in Chapter III and used to determine an initial cost

J. The Qii Criterion is then invoked

q11 = 2a0k31 + b k31 > 0
= -2k + 2a k + b2k2 > o
q33 32 2 33 33
= 2k + 2a k + bzkz > o
q22 12 1 32 32
— 2 k + k + bzk k ) + 2a k + bzkz > o
‘ ‘ (a0 33 a2 31 31 33 1 32 32

incrementing k and k in that order until the

k33’ 32’ 31

criterion is satisfied. Next, the Routh Criterion is used to
check for stability and the kni adjusted as necessary.

If both criteria are satisfied, vJ and “VJ“ are
computed. At this point the new estimate for KN may be
computed using either the gradient or the exhaustive search methods.
In this program the choice was made to use a modified exhaustive
search procedure where k32 and k33 are incremented with

Ak32 and Ak33, respectively, but R31 is incremented using

BL
. . Bk31
k31(1+1) = k31(1) - B “63W

This procedure further expedites the search; and for the example
used, ”VJ“ # O in the region of interest.
The new cost J(i+l) is compared to J(i), and if the

former is smaller, KN(i+l) and J(i+l) are stored, k32 in-

< min k . This is

cremented and the process repeated until k 32

32

followed by incrementing k33, etc., until the whole region is

scanned.

73

If J(i+l) is greater than or equal to J(i), then the

step size 3 is decreased
B<i+1) = 1/2 3(1)

until a lower cost is achieved or until a fixed number (20) of
adjustments in B are exceeded.

After the entire region is scanned, a new contracted
region containing within it the minimum J is identified, and
the whole program repeated until the desired resolution in

identifying minimum J is obtained.

4.2 C-Function Algorithm Program.

The hybrid computer program configuration consists of
an IBM 1800 digital computer and an Applied Dynamics AD 4 analog
computer. The IBM 1800 is used primarily to control the logic
and for algebraic computations while the AD 4 is used for analog
signal sensing and solving the differential equations. It is
clear that numerical techniques, Such as the Runge-Kutta, could
be substituted with a resulting all-digital configuration,
sacrificing computation speed for accuracy. The Fortran program
is written for a general nth order system with L state variable
constraints. The flow charts in Figures 4-2(a), (b) and (c) show
the hybrid system and the program listing is contained in Appendix

B.

4.2.1 Criteria Checks.

Every new estimate of the 5 vector is checked to aSSure

that the system x = Gx is asymptotically stable and the resulting

74

 

Input

2 Routh I No 9-
Criterion 31 Stop

Yes

 

 

 

 

 

 

 

 

t

K Matrix I No
Pos. Def. I Stop
9

 

 

 

 

 

I Max Ixi(ti)I I

 

|_._____._.1

 

Compute
C

 

Figure 4-2(a)

C-Function Program

75

 

 

 

 

No
Newton- I
Raphson
'-' I
Yes

 

 

 

 

 

 

 

 

 

 

 

 

Yes -
C(i+l) < C(i)
?
No
7' No ; GradientI_
y 2 0.004 ? Search -—.®
Yes

 

 

y(i+l) = 1/2 y(i) I

 

 

 

 

 

 

 

Figure 4-2(b)

C-Function Program (con't.)

76

 

 

 

Yes I C(i+l) < C(i) ? I

 

No

 

Yes '

 

 

 

1
d(i+1) = 1/2 01(1) I‘—_ cy 2 0.004 ? I

 

 

No

 

 

 

 

1 _ I 1..

Print
K Matrix

 

 

 

 

vpc ; HVPCH

 

 

 

 

é Print
...—1:...
I Stop

 

Figure 4-2(c)

C-Function Program (con't.)

77

K matrix is positive definite. The program determines the Routh
array and checks each entry in the first column whether it is
positive. The Hayes [HAY-l] algorithm is incorporated within
the program to check the positive definiteness of the K matrix.
Since violations of either criterion will be due usually to some
pathological condition, the program stops in case of a violation

and prints the violation.

4.2.2 Analog Computation

The analog computation, shown by the dashed-line box in
Figure 4-2(a), has three basic operations requiring 2n+1 integrators:
n integrators forthe solution of x = Gx; n integrators for the
sensitivity analysis; and one integrator for signal timing. The
analog connection diagram is shown for a third order system with
a single state variable constraint in Figures 4-3(a), (b) and (c).
It can, of course, be readily generalized to any nth order system
with L state variable constraints.

The system solution operation is shown in Figure 4~3(a).

The primary function is to solve the differential equation

X = Gx ; x(to) = x0

1, P2 and P3 are set with each

The servo-set potentiometers P
successive iteration of the 5 vector so that the final parameter

solution is represented by these potentiometer settings. The

x and x are transferred to

analog signals representing x1, 2 3

the digital computer via analog trunk lines TA 16, TA 15 and

TA 10, reSpectively.

 

 

 

 

 

X.

IC

 

78

 

X ___
3 +1 TA 10

 

 

 

 

 

 

 

 

 

 

 

 

 

/"\J

Figure 4-3(a)

Solution of x

=GX

gé4h1TA ISI

 

 

x 1
F4a-ITA 16I

 

 

 

 

 

 

 

 

79

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(9)

Figure 4-3(b)

Sensitivity Analysis

 

 

 

 

 

 

5X3
AP].
TA 14
6X2
P.
a J
TA 13
1’“ —-i_31
8P

3,-—

 

 

80

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

....2 --—~---' --

J

 

 

 

 

TD 10:

1

 

~~~---- TD 20 2

 

 

b

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11

SJ

Figure 4-3(c)

 

 

 

 

 
   

Analog Logic Control

 

 

 

 

 

 

 

 

 

 

TD 25 I

81

The sensitivity analysis required to obtain the 8 matrix,
equation (34) in Chapter III, is performed by the three integrators
in Figure 4-3(b). This configuration is essentially a duplicate
of the system solution configuration with the exception that the

system is forced, i.e.,

 

 

 

 

 

 

r 6X1 1
BPJ
° Gt, t, GT
. =-e lIle BX.(“F)dT
t J
. o
3’31
LapJ
r0 7 r 5x (t ) 1
l o
. 5p.
B = . ; J
0 =6
L11 -
axn(to)
_ 591 J
Ax1 ax2 ax3
The analog signals corresponding to -——— , _—— and -———
BPj BPj BPj

are transferred to the digital computer via analog trunk lines
TS 12, TA 13, and TA 14, reSpectively. An additional amplifier
was found useful when small signal levels were encountered,
particularily for the constrained variable x3. Compensation for
this "scaling" was made in the digital program.

The signal sensing and timing functions were performed
by comparators C1, C2 and C3, an OR gate and an integrator shown

in Figure 4-3(c). The basic function is to determine the maximum

82

amplitude of each constrained state variable in the time interval
[tO,T] and the corresponding time ti. The final time T is pre-
selected to be greater than the largest time constant in the system.

The first three peaks of each constrained state variable
are observed by sensing the derivative of the variable as it goes
through zero. In Figure 4-3(c), this function is performed by
comparator CI for x3. The procedure is as follows: switch 8
is set so that with the given initial conditions the first peak
of x3 is sensed and its value transferred to the digital computer
via analog trunk line TA 10; simultaneously, the time at which
this occurs is observed by monitoring analog trunk line TA 11;
these values are stored and switch S is changed to the opposite
polarity signal and the procedure repeated twice, each time storing
x3 and ti. The digital computer program then evaluates these
signals and chooses the pair [x3,ti] corresponding to the largest
magnitude of x3 which is then used in computing the C-Function
and in the sensitivity analysis. Ordinarily, switching of switch
8 would be under the control of the digital program, but equip-
ment limitations necessitated manual control.

Comparator C2 controls the total integrator running time
which is pre-selected by potentiometer T. Both Cl and C2 are
enabled simultaneously by a logic signal transferred from the
digital computer via digital trunk line TD 21.

Potentiometer ti is servo-set and comparator C3 then
controls the integration interval of the sensitivity circuit.

All three comparators and a logic signal from the digital computer

via digital trunk line TD 27 feed into an OR gate so that any one

83

of these signals places all integrators into the ”Hold" mode.
Simultaneously, the digital program is signalled of a "Hold"
condition by digital trunk line TD 10. Comparator C3 is enabled
by digital trunk line TD 20 from the digital computer. All three
comparators have a common initialization, E-CMPK, controlled by
digital trunk line TD 22.
Digital trunk lines TD 23, TD 24, and TD 25 switch the
forcing inputs required for the sensitivity analysis. Digital
trunk line TD 26 controls the ”operate" mode of all integrators,
and a logic ”0" on both TD 26 and TD 27 is used for the "IC" mode.
Summarizing, for each iteration on p the system x = Gx,
Figure 4-3(a), is run once for each of the L constrained variables
to obtain the maximum amplitudes and the corresponding ti's.
Then for each t1 the sensitivity computation, Figure 4-3(b),
is run n times. This correSponds to L(l+n) total analog runs.
In the case of the third order system with a single constrained
state variable, a total of four analog runs are required for each

iteration.

4.2.3 C Computation.

This digital computation is straightforward, using
C = '2 [Xi(max) - Xi(ti)][xi(ti) - xi(min)]H(hi) .
161
4.2.4 5 Estimation.
Both the Newton-Raphson and the gradient iteration schemes

are incorporated in the digital program. Initially, the search

for a parameter solution utilizes the modified Newton-Raphson

111.11“ Hilfihllllinli

 

84

algorithm
V C(i)

sum = 5(1) - v 2
iivpC<i>H

The step size coefficient y was found very useful in preventing
oscillation, particularily near C = 0 because of the accuracy
limitations imposed by the analog computer. In the program,
C(i+l) is compared to C(i) and if it is smaller, the Newton-
Raphson scheme with v = l is continued; but if C(i+l) is
larger, the smaller C(i) is stored while v is adjusted by
successive halving until either C(i+l) becomes smaller than
C(i) or v < 0.004. The latter limit was determined by experi-
ment to be consistent with the accuracy limitations of the analog
computer. At all times the smallest value of C and the correspond-
ing 5 vector are stored.

To insure safeguards in the event VPC « 0, the search
automatically switches into the gradient procedure

v C(i)

i5(i+1) = 13(1) 1 a 11vp0<i>i

every time v becomes less than 0.004. The step size coefficient
a is determined by successive halving as in the case of v.
When 0 becomes less than 0.004, the K matrix is computed for
parameter vector 3 corresponding to the solution closest to
C = 0. The vector KN obtained from the bottom row of the K
matrix yields the desired feedback coefficients.

In the event that VPC a 0 faster than C a 0, then it

is possible that the procedure will not continue to the solution

85

at C = 0. In this case, after a finite number of iterations,
the program automatically switches to the gradient method; and
the extremum is located in terms of the pe vector. The user
can then choose a new starting 5 vector sufficiently far away
from Be and in the direction of C = 0 and continue with a

new search using the C-Function Algorithm.

CHAPTER V

EXAMPLES

A problem of current interest is that of automatic spacing
of vehicles in a single lane, i.e., longitudinal control of a
vehicle [ATE-3], [PEN-l], [MIK-l]. Assuming a given desired
spacing between two vehicles, the objective is to design an
optimal control that would accelerate or decelerate the follow-
ing vehicle to maintain the desired headway distance and velocity
relative to the lead vehicle. The vehicles chosen for this example
are assumed to be conventional automobiles.

Blackwell [BLA-l] has shown that the longitudional dynamics
of a typical automobile may be quite accurately approximated by

a simple transfer function of the form

V(s} _ k/T
e(s) — s + a (1)

where s is the uSual Laplace complex variable
V is the vehicle velocity

e is the throttle position

1

k = __—--k‘
2 3

k1 is the ratio of rear wheel traction force to throttle
position

k2 is the ratio of rear wheel traction force to vehicle

Speed

86

87

 

k3 is the ratio of road load to vehicle Speed
M
. e
T = time constant 3 k +——-
2 3
a = l/T

Me = effective mass of the vehicle

Actually, equation (1) represents the small-signal transfer

function since k k k and Me are not constants over the

1’ 2’ 3
whole Operating range of the vehicle. For example, at 40 mph k
may have a value of 240 sec.1 and at 60 mph k may be reduced to

-1 . .
165 sec . For this example nominal values were chosen

 

 

 

 

 

 

 

 

 

 

 

k = 200 sec-1
T = 20 sec
10
d = —-—————— 2
an V(S) s + 0.05 ( )
I
... H i
I-r D -—-——~h1~t___ I. ——q’!
: yl ' yz :31
y2 y1
y2 y1
y2 y1
Figure 5-1

Vehicle Longitudinol Spacing Problem

Figure 5-1 is a schematic representation of the problem

where y, y, and y are the position, velocity and acceleration,

.liN’im H. .... .1..—

88

reSpectively, of each vehicle. The distance HD iS the desired
Spacing between the two vehicles. There is justification based
on traffic studies [TRE-l] that drivers tend to choose a headway

corresponding to the relation

HD = C0 + 2Ty1 (3)

where cO is a constant and 1 is the reaction time of the
driver. The actual choice of HD for an automatic system would
maximize traffic flow consistent with the driver ”psychological"
comfort.

Assuming that the lead vehicle is moving at a constant
velocity, equation (2) is transformed into the regulator form

in state Space by defining the state variables

x1 = y1"y2 ‘11)- L

 

 

 

 

 

 

° =- =- -be
X3 y2 ax3
'1 'fi r
H .1 F0 1 0 r3: 0 1
1 1
:22 = 0 0 1 x2 + 0 <3 (4)
x 0 0 -0.05 x ~10
L. 34 L- JL34 L- J

 

 

Equation (4) is in the form compatible with the algorithms described
in Chapters III and IV such that constraints on the vehicle accelera—

tion can be invoked:

10.0 = x3(max) 2 x3(t3) 2 x3(min) = -10.0

89

The autonomous form of equation (4) becomes

 

 

 

 

 

 

”.1 ' a r ‘1
X1 0 1 0 x1
x2 = 0 0 1 x2
X3 I'Pl 'P2 “P3 x3

L a 4L 4

where p1 = 100k

31
p2 = 100k32
p3 = 0.05 + 100k33

5.1 Example 1.

Assuming initial conditions

 

 

 

 

r- ‘1
x (t ) F450 1
1 o
x2(t0) = -60
x3(t0) 0
L A L J
minimum values of k3l’ k32, and k33 determined by the Routh

Criterion

0.0001

min k

31

min k 0.0006

32

min k33 = -0.0004

and maximum values based on a preliminary exploration using the

Minimum Cost Algorithm and the C-Function Algorithm

max k31 = 0.0226
max k32 = 0.0715
max k33 = 0.0488

90

Table I summarizes the results of the search for an absolute
minimum of J using the Minimum Cost Algorithm. It is to be
noted that after the fourth contraction of the search region,
the neighborhood of the absolute minimum was located, with the

following lower bounds:

k 2 0.0001

 

31
k32 > 0.0029 - 0.1(0.0029 - 0.0006) = 0.00267
k33 > 0.0083 - 0.1(0.0083 + 0.0004) = 0.00743
and parameter values
p1 = 0.01
= 7
p2 0.-9
p3 = 0.88
TABLE I
Search for Absolute Minimum
Step k31 R32 k33 J
Initial 0.0226 0.0715 0.0488 13,956.9
1 0.0001 0.0170 0.0192 72.1
2 0.0001 0.0088 0.0133 26.9
3 0.0001 0.0046 0.0105 10.6
4 0.0001 0.0029 0.0083 5.2
5 0.0001 0.0029 0.0083 5.2

The parameter values were then used as the input to the

C-Function program to check for constraint violations. Since

91

max Ix3(t3)I = 10.75

as determined by the analog run, the constraint on acceleration
was violated and the search with the C-Function Algorithm pro-
vided new parameter solutions as shown in Figure 5-2 and Table
II. Accuracy limitations of the analog computer resulted in the

closest Solution

max Ix3(t3)I = 9.75

 

 

E__82££E£
r' ‘w
0.000039976 0.00012159 0.00013827
0.00012159 0.0037216 0.0028910
0.00013827 0.0028910 0.0082936 J
L, .
J = 7.46
9 Matrix
r 1
0.0000019 0 0
0 0.00059 0
L 0 0 0.0019 J

 

 

92

(t3)

15

 

10~_:N_.______.__

 

-10 .—— —_— .__ “a. ... ... ... .__ .__. ___ ... ... ... __ __

-15

 

ITERATION

Figure 5-2

F-Wnnprinn Thornrinn — Fvnmnlp l

93

TABLE II

Results of C-Function Search

Iteration

 

1 2 3
p 0.01 0.0136 0.0138
1 39
62 0.29 0.2891 0.2891 I I
p3 0.88 0.8794 0.8794 I
1 . 5 1 .25 . 5 I
x3(t3) 0 7 0 9 7 1
1
t3 2.81 2.40 2.41 i I
is?
c -0.156 -0.051 0.049
AX3
-4.75 -4.03 -4.83
591
5X3
——— 1.08 0.78 0.23
892
8x3
——- 8.03 -0.08 -0.08
0p3
49— 9.84 8.25 9.41
8P1
;§_ -2.31 -1.59 -0 44
AP,
59— -1.73 0.0154 -0.0146
5P3

The actual Steps required within the C-Function Algorithm
to obtain the results in Figure 5-2 are graphed in Figure 5-3.
The circled points represent the successive acceptable estimates
while the rest of the points demonstrate adjustments required in

the y and a step size coefficients to obtain rapid damping

94

 

ea

H mHQmem u coflusaom coHuocsmuo mo mouUmH: muwmuamum

m-m muswae
NH 3 m o q

u

'4
.4
-4
‘4

AMHm

 

HZMHQ<mw ZOmmm<muzOHBmz

 

 

,Q.01

N.Ou

N.O+

Eé

o.¢+

uoinoung-Q

 

J". . .
LL " (I'm-fl

95

of the oscillation as C a 0.
Since a solution at C = 0 does not necessarily imply

the final solution, a neighborhood was chosen about the point

k31 = 0.00013827
k32 = 0.0028910
k33 = 0.0082936

such that

I!

max K [0.00013827, 0.0147, 0.0192]

[0.0001, 0.0006, -0.004]

min KN

With this input to the Minimum Cost Algorithm 3 new minimum was

located at

k31 = 0.0001

k32 = 0.0020

k33 = 0.0133
J = 2.5

This KN vector yielded the following parameter values:

p1 = 0.01
= 0.20

p2

p3 = 1.38

where were used as input to the C-Function Algorithm to check

constraints. The results were as follows:

96

2
max Ix3(t3)I = 5.25 ft/sec

 

 

t3 = 3.38 sec
K Matrix
r- ‘
0.00002 0.000138 0.0001
0.000138 0.002622 0.0020
I-0.0001 0.0020 0.0133 J
9.. 1131.213
1
0.000001 0 0
0 0.000124 0
O 0 0.01512j
L

 

 

Since the constraints are satisfied, this is the desired solu-

tion within the limits of resolution in the search region.

5.2 Example 2.

The possibility of discontinuity in p-Space was dis-
cussed in Chapter III. This example is typical of the results
obtained when the largest amplitude of x3(t) Shifted from the
first peak to the second. Figure 5-4 demonstrates the behavior

of the C-Function and x3(t3) under the following conditions:

5.0 2 x3(t) a -5.0

 

 

 

 

r " r fl
x1(to) 450
x2(to) = -60
x3(t0) 0

J L .

r p1(0) r- 9.996

[32(0) = 18.26

 

L p3(0) 12.68

+2

-10

-12

97

 

 

 

 

 

 

 

r '1
1
x3(max)
. __. ._ __ _. _. ... __ —- 2+5
1. ... (I)
x3(min) r
1.. -- - — — -— -- — — -— —- -— —- — — — “‘1 ‘J
A. , -1O
’ ’0". V‘— "V 1
p “-15
C
p .1-20
» . -25
1
L. -1-30
1' .1 -35
$ ITERATION
1 a J 1 L L
0 1 2 3 4 5 6

Figure 5-4

Solution Stability under Peak Switching - Example 2

98

for the system shown in equation (4). The initial values and the
first iteration values were extremely large and are not shown in

the figure. There was no oscillation of the C-Function as the
dominant peaks changed at the fifth iteration, and final convergence,
within the limits of analog accuracy, was obtained on the sixth
iteration. Table III Summarizes the values at each step in the

Computation.

TABLE III
C-Function Behavior Under Peak Switching

Iteration

 

2 3 4 5 6
p1 3.577 2.997 2.713 2.5913 2.2616
p2 17.60 17.663 17.694 17.708 17.805
p3 12.18 12.218 12.23 12.229 12.235
x3(t3) -35.75 -18.25 -9.75 6.75 5.75
t3 0.18 0.17 0.16 2.56 1.83
c -11.25 -3.08 -0.70 -0.21 -0.08
,‘jX
-—-§ -2.97 -2.93 -2.93 0.43 0.18
891
dX3

0.32 0.32 0.32 -0.13 -0.13
151132
7}}(3
.——— 1.97 0.92 0.48 -0.08 -0.08
5p3
;g_ -21.27 -10.68 -5.70 -0.57 -0.20
1“ p 1
99- 2.32 1.19 0.63 0.17 0.14
3P2
LE— 1.41 0.34 0.09 0.01 0.009

r

0.051148

0.35153

0.022616

 

L

r
0.051148

0

L O

 

99

t__me2£££
0.25153
2.82691

0.17805

J = 23,241.5

91.889525
0
2.46712

0

1
0.022616

0.17805

 

0.12185 .1

 

1.14083 .J

CHAPTER VI

SUMMARY AND CONCLUSIONS

An algorithm has been developed resulting in a computa—
tional procedure which determines the optimal scalar control for
a linear, time-invariant regulator with state variable constraints
and a quadratic performance index. The basic contribution which
makes possible a practicable engineering solution to this problem
is the development of a simple algebraic algorithm for the deter-
mination of the gain matrix K. It has been shown that the K-
matrix can be obtained from an n-parameter Optimization of the
closed loop system. The need to solve the matrix Riccati equation
is eliminated, and, therefore, there are no computational stability
difficulties which impose a practical limit on the order of the
System.

The K-matrix algorithm also provides a criterion for the
positive definiteness of a diagonal Q weighting matrix which is
used as part of the search procedure to obtain a minimum cost
solution. An efficient digital computer program allows only the
positive definite constant K matrices which represent stable
control laws and define positive definite Q matrices to be con-

sidered in the minimization of the quadratic performance index

J = l/2xT(t0)Kx(tO).

100

101

The result of this minimization is then evaluated for
state variable constraint violations with the aid of a constraint
function specifically developed for this purpose. This constraint
function is an integral part of an iterative procedure that locates
the K matrix nearest the minimum cost solution and which also
satisfies the constraints. The new value of K is then the basis
for the next estimate in the cost minimization program. Successive
iterations produce the desired K matrix.

An example of the application of this algorithm to a third
order system is presented, showing the results obtained in defining
the optimal feedback coefficients for vehicle Spacing control with
acceleration constraints.

Since the K matrix contains all the necessary engineering
design information, the procedure described in this thesis should
have greater application than any of the currently available methods.

There are obvious extensions to this work which should be
considered. It would be desirable to investigate the possibility
of developing a correSponding K matrix algorithm for the general
case of a vector control u(t). Additional research needs to be
done on the sensitivity of the cost function J to state variable
initial conditions. For example, it would be desirable to obtain
a solution for minimum J such that J has minimum sensitivity
to the initial conditions. Another possible extension of these
results is to the work of Kleinman, Fortmann, Athans [KLE-Z]
where a time varying matrix K(t), representing a finite time
cost interval, is approximated by a series of piecewise-constant

K matrices. Finally, it would be of practical interest to

102

investigate the computational limits imposed by the errors inherent
in an analog computer. For example, in the course of this re-
search it was found that potentiometer setting tolerances imposed
a limit on the accuracy with which the constraint function could

be obtained in the vicinity of C = 0. Results of such work would
show the conditions under which it would be advisable to consider

an all-digital configuration.

REFERENCES

[ATH-l]

[ATH—z]

[ATH-B]

[EEK-1]

[BEL-1]

[BER-1]

[BLA-l]

[BLA-z]

LCES-l]

[CHE-1]

REFERENCES

Athans, M., and Falb, P., "Optimal Control: An
Introduction to the Theory and Its Applications",
McGraw-Hill Book Company, 1966.

Athans, M., and Levine, W.S., "On the Numerical
Solution of the Matrix Riccati Differential Equation
using a Runge-Kutta Scheme", Report ESL-R-276,
Electronic Systems Laboratory, Massachusetts Institute
of Technology, Cambridge, Mass., July, 1966.

Athans, M., Levine, W.S., Levis, A.H., "On the
Optimal and Suboptimal Position and Velocity Control
of a String of High-Speed Moving Trains", Report
ESL-R-291, Electronic Systems Laboratory, M.I.T.,
Cambridge, Mass., November, 1966.

Bekey, G.A., Karplus, W.J., "Hybrid Computation",
John Wiley and Sons, Inc., New York, 1968.

Bellman, R.E., "Stability Theory of Differential
Equations", McGraw-Hill Book Company, Inc.,
New York, 1953.

Berkovitz, L.D., "On Control Problems with Bounded
State Variables", J. Math. Anal. Appl., vol. 5, 1962.

Blackburn, T.R., "Solution of the Algebraic Matrix

Riccati Equation via Newton-Raphson Iteration",
pp. 940-945, JACC, 1968.

Blackwell, L.M., "A Study of the Adaptivility of the
Automobile to Automatic Control”, Report EES-276A-3,
Communications and Control Systems Laboratory, Ohio
State University, Columbus, Ohio, February, 1967.

Cesari, L., ”Asymptotic Behavior and Stability
Problems in Ordinary Differential Equations", Berlin,
Springer, 1963.

Chen, R.T.N., Shen, D.W.C., "Sensitivity Analysis
and Design of Multivariate Regulators Using a

Quadratic Performance Criterion", pp. 229-237,
JACC, 1968.

103

\rHN—l]

[HAY-1]

[KAIrl]

[KAL-Z]

[KAL-3]

[KEL-l]

[KLE-l]

[KLE-z]

{KOE-l]

iKRE-I]

[LET-1]

[McG-l]

104

Fenton, R.E., Cosgriff, R.L., ()lson, K.W.,
Blackwell, L.M., "One Approach to Highway
Automation", Proceedings of the IEEE, Vol. 56,
No. 4, pp. 556-566, April, 1968.

Hayes, R.M., Jr., "An Algorithm for the Determina-
tion of the Definiteness of a Real Square Matrix",
IEEE Transactions on Automatic Control, vol. 13,
no. 2, pp. 219-220, April, 1968.

Kalman, R.E., "Mathematical Description of Linear
Dynamical Systems", J. SIAM on Control, ser. A, vol.
1, pp. 152-192, 1963.

Kalman, R.E., "When is a Linear Control System
Optima17", pp. 1-15, JACC, 1963.

Kalman, R.E., Bertram, J.E., "Control System
Analysis and Design Via the Second Method of
Lyapunov", J. Basic Eng., v01. 82, pp. 371-393, 1960.

Kelley, H.J., "Method of Gradients", G. Leitman,
ed., Optimization Techniques, Chapter 6, Academic
Press, New York, 1962.

Kleinman, D.L., "On the Iterative Technique for
Riccati Equation Computations", IEEE Transactions
on Automatic Control, vol. 13, no. 1, p. 114,
February, 1968.

Kleinman, D.L., Fortmann, T., Athans, M., "On the
Design of Linear Systems with Piecewise-Constant
Feedback Gains", IEEE Transactions on Automatic
Control, vol. 13, no. 4, pp. 354-361, August, 1968.

Koenig, R.E., Tokad, Y., Kesavan, H.K., "Analysis
of Discrete Physical Systems", McCraw-Hill Book
Company, New York, 1967.

Kreindler, E., Sarachik, P.E., "On the Concepts of
Controllability and Observability of Linear Systems",
IEEE Transactions on Automatic Control, vol. 9, no. 2,
pp. 129-136, April, 1964.

Letov, A.M., "Analytic Controller Design II",
Automation and Remote Control, vol. 21, pp. 561-568,
May, 1960.

McGill, R., "Optimal Control, Inequality State
Constraints and the Generalized Newton-Raphson
Algorithm", J. SIAM on Control, vol. 3, pp. 291-298,
1965.

[MlK-l]

[MGR-1]

[IER-l]

[PON-l]

[REK-l]

{SAG-1]

[TOM-1]

[TRE-l]

[TUR-l}

[TYL-l]

[VAL-1]

LVan-l]

105

Mika, H.S., "Optimal Control for Automatic Vehicle
Longitudinal Guidance", Report 1206-1, Ford Motor
Company, September, 1969.

Morgan, B.S., Jr., "Sensitivity Analysis and
Synthesis of Multivariable Systems", IEEE
Transactions on Automatic Control, vol. II, no. 3,
pp. 506-512, July, 1966.

Perkins, W.R., Cruz, J.B., Jr., ”Engineering of
Dynamic Systems", John Wiley and Sons, Inc., New
York, 1969.

Pontryagin, L.S., Boltyanskii, V.G., Gamkrelidze, R.V.,
Mishchenko, E.F., ”The Mathematical Theory of Optimal
Processes", Interscience Publishers, New York, 1962.

Rekasins, S.V., Hsia, T.C., "On an Inverse Problem
in Optimal Control", IEEE Transactions on Automatic
Control, vol. 9, no. 4, pp. 370-375, October, 1964.

Sage, A.P., "Optimum Systems Control", Prentice-
Hall, Inc., 1968.

Tomovic, R., ”Sensitivity Analysis of Dynamic
Systems", McGraw-Hill Book Company, Inc., New York,
1963.

Treiterer, J., "Improvement of Traffic Flow and
Safety by Longitudional Control", TranSportation
Research, vol. 1, pp. 231-251, Pergamon Press,
Great Britain, 1967.

Turnbull, H.W., "Theory of Equations”, Interscience
Publishers, New York, 1957.

Tyler, J.S., Tuteur, F.B., "The Use of a Quadratic
Performance Index to Design Multivariable Control

Systems", IEEE Transactions on Automatic Control,

vol. II, no. 1, pp. 84-92, January, 1966.

Valentine, F.A., "The Problem of Lagrange with
Differential Inequalities as Added Side Conditions",
Contributions to the Calculus of Variations 1933-37,
University of Chicago Press.

Van Ness, J.B., Boyle, J.M., Imad, F.P., "Sensitivities
of Large Multiple-Loop Control Systems”, IEEE

Transactions on Automatic Control, vol. 10, no. 3,
pp. 308-315, July, 1965.

[WAN - 1] :

LWON-l] :

106

Wang, S.J., "A Gradient Computational Technique for
a Class of Optimal Control Problems Subject to
Inequality Constraints", Ph.D. Dissertation,
Michigan State University, 1969.

Wonham, W.M., Johnson, C.D., "Optimal Bang-Bang
Control with Quadratic Performance Index",
pp. 101-112, JACC, 1963.

APPENDIX A

SAMPLE MINIMUM COST AUSORITHM PROGRAM

TASK l

...... .-

// JHF’
// fin?
“LIST

800 TSX V3M494

*- - s o- ...--.....“ -....- . _ . .H . 7 .... 7.

 

 

- - .. ...- .. .-..»._*.—.-._...~._ - ».. -. 0‘- . -.. .-..-uw——..—..—- ... . ,_- -..-

.. ...—.-- .-.. ..- 1,1,, ... .. .<W - -m—v. ..-...h. ....

2‘". IN K
ILL

22mm :F-IP‘RHCF- F78 PRUGRMI

3:11 11‘ r 1
*lHSb

5H0

Ilwfiv It iEi%?RS

(CAHU,1443 PRINItR,UISK1

REAL 1’~(5931
“17f:l\bltm1'A(319>“1131

1’51111129 511011..91'19\l111191ml];2 1,1".
FUVMh111295F10031 &

1<r:/‘111/9 5011K13‘ 9.119K(.'9/19K139-1

.Libtz,:01)x0111.xu(2),xo(:)

FHHrn113F10.4)
Nt£Ui295061RK3mx,RKZMX,RKimx9RK3MN,RK2MP,RK1tL
1%Ieflalitfi310.41

“:1 1- -. 1 H --,-_-

J21

“C=1.0 -- *‘* '

nKl=O.l*(RKlMX-RK1NF1

eK2=u.l*iRK2flx-Rk2tt1

31K32iiol‘i‘ (RKBI‘X-RK' ... :111

I.(L.’~931=RK3MX

1(29212‘11K2MX

'(5911‘11111214'11’ """ 2«1* K139711*(A(31+(ti’i"i‘21"i'31.(5'9211‘[‘~(211/iw-‘-’----}‘-‘G.LmUl
“(391)-12K11vX1(19719.12

1--’911=1-’\K11'~'1XK‘ ’

1119113A111*v9a1+A121‘K1J9l1+1V**21WK1)921 '(3911

191121—Ai11«K(3951+Ais1'Kia911’rii 21’ K(-9111 (19))

(2421:“K(39 l1+/-‘(21‘3 K11H931+A131K15921+(H* */’1 1’1’9/1‘11313951

111123111911’3HX0111«‘21+K12921’v‘(X[l(21*’"21+K(.‘1931*(Aft?( 1* c1+r.t

] ( 1:;911*WX11(11 thib1-+K(359 21* X11(21 X1 (3 $1'*K(l.921 ’bi 11‘fJJ(21 1

1108-

VN=0.5*VN
""“" 1 t11=2.°f)*f\1“11*‘K1'3'9‘T‘1‘¥1‘B**2‘1*K1391'1**2"“”"
IF1U117397397
1.3 'W1‘iITs'51395031’"‘ ' """ ' ' “' “" "”"
508» F1JF‘8'1AT113H CHANGE RK1MX1

(”7211 TH - 80-2.2-0---” . .-.. .--. -.- ., _ ...--
7 U3=“2. 0%‘K139 21+2.U*A151*K15931+1R ~~~~~ 1*K15931 """
"W'11‘11-131f59‘7596 _. H ' ‘_-n
75 K13921=K13921—DK2

“-"- IF(K(3,2)*RK2MN)SB,7}7“'”' - -“"~ ”~*

«3 K(395)=K(3931~DK3
~ IF1K13931-RK3WN1809 40*, 0*"‘ ~ ~
5 n2=-2. K11 2)+2. o A(2) K(5,2)+(H«22)»K(5 2)r*'

'I~'(H”2)7n,74, 3“- --»-—- ---- -
74 <1391)=K(3, 11-DK1

11:1”13911"R‘K1M'K’111‘391‘7' 1.7.... , ...-2.--.-. .-. .
1’1‘’1121"(1131(13951‘1/5141 K139111+1R’*"*‘21*K13,1)?;--l-;(3,3)

1711 T11 ’ -- —-- - - -. -._.-. .-

10 L139212K13921“UK2
“' IF(K139ZV“RK2MN105954954“ ' ' ”
H P3=A(b)+(B**21*K(3931
'IF19317697699 “
7b MRITE1395091
5119' 141113154 AT113‘H CHANUF K K,31*’il\'1'
(1C1 '11} 811
“9'P2=A121+18**21*K13971 "'
18 P1=A'(l1+1B**21*K(3911
“‘ ”R3=(P2*P3=?11fP3
.(R3154954910
10 1F(P1177977982
77 WHITE1395101
‘3 1.‘fr‘13“t.)1~tF-*A‘T‘1 13H"'CHAN (B E; 141% MN )“
CU T11 80
‘ 82'1r(L~ ~215985, 85 " “ ‘ -
b=1’(1911"(X0111’h12)+4<(2971*(xh(21 """ «2)+V13 1 (X1 (31*W2'1+2.W1
1(1<1391.1“Xf1111 “X[}(.$1+K13 921 ‘”X1(2i1 ’X1113141K119 1 M1 (11 >U(2 )1
VN=0.b’VN
'“‘lb PAR1=1A121+18**21v-K139?11«‘TX[(11**Z1-XU(21 """ M2+7 U XU111 XH€31
1+7.0 A131*X0(1)=- X112)

35 PA22=(A(11+(B**2)*K1391))*(XU(1)**2)+(A(5)+(H*$2)*K(39511*(ku(21**
12)+2.U=#XU(2) =XO(3)
PflR5“1A121+(BH“21 K(3,2))*(XU(2) ...... 2)+x0(3) “’+/.u*(n(l)+(hmw2)*
1K139l))»XU(1) XU(2)
SUFLzszTf9AR1**2+PAR2**2+p223**?r»

30 CURlZC’kPARl/SDEL
"'CHR2=C*PfiR2/SHEL
CKK33C*PAR3/S11LL
1F1fiBS1CQR'11‘“A 1‘511‘159111
36 IF(AHS(CUR21-AHS(F(39211
"'1‘7 If:1Afi5'1CU1t3'1"/\1’>S1K139311
;2 CT”.D*C
1 T1: 38' ‘
5 K1.911:K13911‘1NPH1l/Q1'tL
1? 1K1-.9 11 -4<11111.'1 {1’ 9 {(19 {11
'1'1-411 911-1<K11‘H
71.1 1 11911“A111"K159)1+A121K13911'1113‘1'5'2131K15'921 113911
n11921‘A111‘K139”1+A151*‘K13911+1H*’21’K1“911 K159i)
1 921““K15911*L1121‘1K‘151931‘1’A1 K19921+1 31/) 1159/ 1:1'1'\1:'27"J1

"169

vm=u(1,1)*1x0111**2)+x(2, 2)“1’1(/) ------ z>+x(5.l) (x013)**z)+2.u*

‘~-11K133115x0111*13131+«1321*xw11 )zxnt3)+511 Wfl1113XU12>>
\‘mr-U . b =i=v N
“1((VN-VUTfir“79;79“‘“'"'““““'“" *'-~*'"

79 l’(391)=K(3,1)+C PARI/SDtL
—2 CUJJ 5*c ---~‘ ~
J=J+l
"“ “II-"(J-“2013y3'; 100' ‘
81 L=L+l
' 171] T1! 1 "" -.--..
5‘5 V1h111. (315()2)K‘371)9K(39&)9F(-QJ)9V|"
502 1:11in Ana-150:4 F2111") " ‘
11 11(A:)5(VN-VO)~1.OO()1)12912941
13+ I—l
N=M+l
1|:(7‘~1-10)15,15761
(>1 WHIFEKD, 507)
5177 Fun-1 MUM“ N)
1;; 11(3,Z)=K(3,2)-Ui<2
'C=1.('J ' " ‘
(’15? J—

r -.

I:
I1;\i/(_,,d)-RK‘{MI\‘)"5,- 11,54
“33 ll(_19.)):i\(3 3)-1K3
hf) 11(“(59 3"”RK3I’N) '~’),’1(‘940'
1."!11 1.'] 1(39503)
>305 FWVNI‘ T(7H ‘J LUMP)
:1' TI‘ 12
'1) CAN" EXIT

r .
;- INF 1 \

APPENDIX B

C -FUNCTION ALGORITHM PROGRAM

HIM IfEbHY iflrflhtfi'
TASK 1500 TSX V3M4.4

-.- .. . - . ..-. ...”k- . ..-. ..~: . ..-: .-v- .. -7- . ..1 . 1..-... ~‘ 4 . , . ... ....

// JUN
// HIP? HIV}.
3LI§T AL!
¢HnerwCISS PRUGRAN’
vnmr ~HQH INTEGERS
’TIIJCS (CZ/WI. I91443 ‘PRIIITFI’WI‘IISKT
1. I11:IIJN RI690I91’KI59‘7I9SI595‘I9W(5'91)yII\(177‘IvI’LIIRIt’IVI’II’I7AI3I7
17H“)9A1'1AXI5I9XIIIVI‘3I9I‘TIcI9XI5I9TI5I9FIEI9I.(-)19~IPI5I9PI‘I’I‘II9HIAIJI9
fiT‘III,HPP(5IQPCHV(berQ]P(HI
C'\I-I_1":YI'.NT
SAIL 11Y1
WEN? (291IIOI'I‘J9TF9IAIII9 1:19 IQI9I7
W95" (Z9200I(P(II9I=19I\ I
3“}“/"II (292013‘I’(XU(1191=19E‘)9 SCALI’
Mil-‘11.(29400IIXIW/1XIII91=19I\I9(I\1 II"(II91=19II
111’.) 1‘4111f+"L\I'I12‘9“7F1f‘.3
gwu HHHmhIIBFlflg4I
11f“) F11E'1I"l/3T‘IIIIF'R‘bZI
NH 1 I=l9N
II IPIIIIIyzyl
1 :n'i‘IIIIJIIITT
,’\1‘~9::1\‘
I'1\Io‘,"‘.I‘I/2‘.O
kat‘lm.
1'“. 2f 121-1 1'.

1:]_.(:
: I.“
‘2'"21
”10011.0
’2 I. II‘I-NIJ)3,4,3

l2 \,I=I.‘+lI/2

:1'

 

. ..-—.p.~-.—-—.—-...»--.- .1fi-1-..<- u.- “.7

GU TL) 10

- "4*";1—11/'2+'1*~“5"*”"5' .-- - _ _ W . _
10 R(19 11:1. 0
1m zor 1:20, 21 ~ - ~ .-_,_.._-_~_._.._.~ _ .- 3...- 53.,

201 CALL HYUO(I:O) ,
“”'CALL HYDOf22911W
CALL HYDU(2690)

CALL HYDUIY79UI” ' ' “*‘55‘

CALL HYDLYIbOI

CALL HYUUIZZan‘ A I IV“—
PAUSE;

RI79TI=‘PINI"

I’d-‘31

1311 51:29:] ' " ' ' ‘ ~ 5 ‘-

I<I<:'\‘I-’I’\
‘RII9II3-PIKKI
b K2K+?
“19:? _
IF Iu-lUIll 12911
13 I- 1’1” '
I1

IHI A I: 29 J
I/.I/="I—K
13 :(Zle:‘PIKKI
C) 3354‘2 ‘
I‘I11r-‘1‘xi-l

‘IFIRLZ911I9199197
7 III} «”1 1:39I‘I
MU H K:1,J
H(I9KI=IRII-l9lI*K(l-29K+lI-RI1-2911*R(1-19K+1II/R(I-l9l)
IE1H1I91II929929H * ' "““”
h CUMTINUE
’HIN+191)=RLN-I921
IFIRIN-19ZIIUB9939H
:11) 14' 151751"
AKIM9II=-(A(I)+P(I)I/(H**?)
NKII9NI=RKINHII
IIII IL) 1:19NN
IFII‘1718917318
17 Mn 15 K229N
RKII9K-1)=A(I)*RK(W,K)+A(K)*RKII9NT+RKII9NIIRKIN9K)*(b**2)
1m z-’1’(1-1--l9II=1;'.K(I9I\'-11
8H Th 15'
II”: IIii I") K339“

9...:
1\
+

‘(1‘K'11=-RKII-1HQLHHI1*RK1N9KI+K1KVWUH19' +5K111w) KKL'HU*1

1(C 21
IN :KLX- I9I)tRK(I9F-1I
12 CUHTINUE
‘ DU 34 1:19N
DU :4 J~19W
A. S119J):1-:1«119J1
1m :13 r-..—.1,r\1
Sr:1;‘=§5(1-1911)
1)“ .3/1- I-1V19HI
1 (“55131191I)-Ah§(§1AX1)?4 ?49 ZR
25 :511911
L: l
24 CH :111LI
I (L- 120926927
1’1 233 1:149“!

TC(I9L)=S(I9M)

3119m1as119t1*"“" '**“““
25 S(I9L)=TC(I9L)
~ ~nu 29 ltM9N”“ 1.- -~ ~-*~

TIEIL91I=SIM9II
"?(fi9l)=StL9TI' "" ""
2v 8(L9I)=TR(L9I)
26 NH 30‘12M'N“”“""
in ‘(n9l)~S(M9I)/APS(91A/‘
1111M +1 '“““""
PH 31 I =MM9N
~v-nn 31 J: M9 I "‘
31 SII9J)= >(l9J)-E(I9w) 3(M9L)*S(M9J)
” 'HH’BZ I: M9NN " "'” " '
Ill=l+l
I111 32 J=III"9i‘-‘1
32 3(19J)=E(J9I)
113(Sham-’7)3fi49'~)49?3?3
f; CMPTIWUL
9% r2 1:1,m
i~('11E(!’(I) 1-111.H),329;fi5923
3 '1’ I1I1I39anII
9111;; I-zlw ".TIIQI
‘" PT1112091*(ARS(P(I)I-LO.O)"
0H 1H 82
22 PTtI)=U.1*AHS(P(I))
15.2 {,1}! 71M”?
VPI11:(398fH11‘PT
PM 12L K219N
121‘ Jr<1<)=11100¢1.0i#)111<)
(LIL PUTSIl9Zlb9JVIIII
1nLL‘PU15(19221 JP(I)J
LALL PHI5(19 ZUE9JP(2‘))
CALL ruTStl92179 JP12 )1
CALL PUTS(192069JP(3))
CALL PUTSI19 2259 JP(3))
UH an 1:20921
40 CALL'HYUUII9U)
1,1411 11‘1’11U(279O)
CfLL MYDUt229l)
CALI. I‘1YUUI269III
** "CL!1 HYDLY16fl)
CALL HYUOI2290}
CPII'HYUUV21917
LHII HYUUI2691I
IL I‘"{UI(109I'\I
1(1)41941942
44 um 00 J=19M
X(J)=U.0
IF(XHAXIJ))20292039202‘
203 IEIXMIN(J))2029609202
"2fi2 HPITLt395001J "
WHITEI391004I
1(Hw+ Lnutnrf124rith;wm1(71 TR1(1'|U XIHITIJII
”ANSI -
CAIL HYMK
CALL .Y11R(1(929~1(1)91A(]))
CALL IYUMK
IE(1I-TA(1))439439411

‘

-."
H
.9 (\

 

WM:

411 D0 412 I= 293

”*“CALLMIflmWfi?TTTT ””““”””'”"“‘” ”“

CALL HYDO(21 0)

""“ WPITE (3;TUUO‘)“ *“'”“""‘°‘"""”“‘"‘"“"‘ " "" “”"”‘ "
1000 FDRMAT(20H CHANGE C10 PDLARITY)
PAUSE""" " M‘ ” '“’

CALL HYDO (27,0)

‘”“”CALL'HYULYTEU7””“ '”" "”*'"” “ '" “'0'

CALL HYUO(2791)
'CALL"HYUOT?T?T) "‘””“ “"” ' "”“
CALL HYUU(2730)
‘413“CALL"HYDITTOVK)'
IF(K)4139413, 414
414 CALL WK' ““"'""‘ " "‘
CALL HYAIR(10, 2, AM(I),TA(I))
‘CALL HYUMK‘ “
412 CHHTINUE
‘ ""' DH 415 f: 273 "" ""
Ir(ABS(AM(1))-AH5(AM(I))141I9416,416
416WXTJL=AMTTW”“"“' ’“'
T(J)=TA(1)
'GD‘TU ‘4‘1’5""‘""‘
41/ X(J)=AM(I)
AM(1)=AM(11“
T(J)=TA(I)
TA(1)=TA(17'
4 n t) Cuis': MU};
'”‘ITL(3, 9001X1J)

‘17.) 1"(l‘1‘(X(J) )--‘:\HS(XI(J) ) )C)()9‘fn9(()

4“, (1)—X01117

an CHITILUE

44 DD 51 Jilfm“'

1F(XHAX(J))46,47,46
47 IF(XWIN(U?)46,48,4$
48 F(J)=0.0
(3|) TU“5‘1 ”"“' .
4b H(J )=(XMAX(J)-X(J))<(X(J)-Vx XMIK(J))
‘ 1F1HtJ))49}50, 50 " “

49 V(J)=100

., (AH TH 51..--.-4..‘.7.._,...

50 F(J)=0.0

HI'CHNTINUE"‘

Sili'~nl‘=(}.()
DH biw'd 1 N"

52 *W i=SUHF+F(J)

J21

101 11(SUDF)53,54,53
‘ 53‘1F(ITJ))Y?,74,7?

G4 ln~'>’(J)-XH(J))‘:b bt»,55

‘16"(.1)=H.Cn

74 L» M? IK 1,N

5713(J,K11)=0.0”

J2J+1 .
1:(Lbfi*)1()1.1()1,ll)3

Ht) 11"(X'wAX(J))10411051104

105 15(XHIW(J))104,50,10£

104 I?(.H:-']..H

7; ”“1T:(39900)T(J) ‘““

585 11(LKK-3)7159715,102 ——————J

-—.- --_ —..._.._... . . . - ...--.

’/1_5 ilLL: 1(jf).()3‘1'( J )

CALLfiPUTS+172397KK3-- mm ,, -U,,m, .CWw

DU 130 1:20721
130“CALL“fiYDOfffiTr" "" ”
CALL HYDOTZZol)
“ TX) 58 M=Z3VZ5 “” ’ " ’ ’”"
SH CALL HYDU(M90)
“’II=J¥11‘ ““"”“'“ "" ‘”‘“"‘” ‘ w *"
131} 5‘9 1"1323925
' "CALL ’HYD‘UI’Z‘Z? 1’7" ‘ " "' "
CALL HYUUH‘dyl)
'CALL-HYmflt26.0J
CALL HYHU(27,0)
”(JH L HYULYWffl11W"‘“
CAL1_ HYLKJ(22,())
CALL HYUUT2071) ' " ' ' '"‘
CidJ_ LYLKJ12691)
AZ‘CALL‘HYUIC103K1“‘
IF(K)62962963

6'3 JJ=W~22 W.-- ‘ ‘ “' ° ' “ .-..-.
CALL hYMK _
‘ CALL HY’ATRT'T'I‘Vi 13( J, JJ) )' ' T' ' “

(:ALL. HYlHflK
"fiRITEf3,9ODTS(J,JJ)
‘39 CALI. HYUU (M,0)
“J:J+r*‘*““““ ~+ '~~
IF(J—N)101,101,103
103 SHMJI=SUM32” “"‘ ‘
SHLJ2=SUMJ
FHMJ:0:O"“”“‘W“
LN 54 J=I,N
7*; 3081'=“r—i"("3‘)*FTJ‘)“ ‘ “ ’ " ' '“’
04 SHHJZSJMJ+CUST
uRITL(3¢900)SUMJ
2H3 IG(LKK~2)708,YOH,YUZ
708'IF(nHS(SUMUZI-AH§(SUNJT)7P1,702,702 "'
701 fiTSzSUHJZ
' 01‘. 7'03 It’r'm
705 :3T9(1)=PPP(I)
IFtKKK—27711;707;553"
711 DK=0.b*DK
’rlo inn 704“”I=Iyfi!'
704 P(I)=LSTP(I)+DK*PCHR(I)
”KKL22 "-LM_ ...... A . - “w
(m MI 405
n? I?(AHSLFSTSTHAHS(SUMJ11709,709{710‘ “““'
n KKL=1
- (an ’ Ls~7t;z' ' ' "‘
}WW3:H<:0.5*DK
1|:("K‘000041732,712,713“C
12 kw 714 1:1,N
71A V(I)=LSTPTL)
KKKzB
00 TU «05“
702 IF(KKK-2)121,160,122
1A0 IF(AHSLSUMUY3AHSTL5TS)T151.709,709
lfinl L‘YIS=\AJmJ
D0 162 ItTyN
1&2 r?flW’(I)=P(II

115

CH TM 181
1F1fiHS1SUMJ1-KPS11S1S1112491?:51?3
Ki-LKZQ
1S15=3HMJ'
HM 105 I=I,N
1K3 LS1P1I1=P111

an 1U 715

lco’
‘5

’-J

1:) ;=c .v~c -- v -a----~~
11"(L“(1.1101114(1‘5,]()(9121)
12 mm 1:7 1:1,m <5

1C11'1I11C5 PC”P11)
1:" (FEIIJ)I?8,L(1:11/:"j
P1112LSTP1I1~PCUL111
an TM 127
F'% P111=t ‘TP1I1+PCH’1I1
1“? CHUTIMUE "
an Tu 4H5
21 SK71:C.O
1n n; K21,N
“PAR:0.0”"
1‘1"; I'm J:
78 1F1K—?110H,109,INW
1M4 3(J,K1TH.1*S1J,K)
lUM PAUJL1XKAX1J1+VPI11J1-?.05X1J115 F1‘15 S111,K 1
kn frnrrSPAK+DARJ
1K1=SPAR
: I<'11 1 , ,111111 119711 1 K 1

*4

.‘j K“ I" S": fxr’ AR :1: :'- p

'7 “ i=1 LL+PAR2
“H 121“ 1:1,Ix
111=PHP1I)
”PP 1111=P(1)

;:> {F U;KK—21802,802,Z1
3112 J’1' 11 11— $1H7J5W3AF311 1/53H15L
1931 11.1 (:1 '1'
6“ 1K1HK*1. 077199029802
71 11 (I1 C*PAP(I1/SukT1Slfil1
KKI M '— " "
Ir1CUIJ1bQ,GU?1
“7 P11 =11I1+PCUP1I1
MH115
1? V111- P1I1-P1UP111
HUB ‘I11r13y90111PCUP111
47‘}! (_'1‘~"TIM1F
z+I‘-;\ 1 1: I 1'1- 1 :3 ,1; 1111 ) ( 1) 1 I 1 , I ==1 , 1,)
CI:QB.U”P111-6. 0:P121
t111115111‘191},

\
I
11

1:217: 1r :1
113., 11-11'1I"123,5(3(111C
an” uuuup1Istzn.41
147.11 1S1:
1311 11' 77
4H5 1'11HJ1151S1“M“‘1311J11019h19I“
CI '“11111 9001ESTS
CH IU ll
71) 11"1'111-1.‘ 1 3,90015111’K1
71 V111L13960011X11111319W1

4r Z) Y

426

1315‘:

.‘"3

116

“v :1;- (:3,(.-(')«"1)( 1 ( I 3‘,

wRITt(3,bUn)K(RK(I

I:?a:1.u

KKVll‘ ‘ ~

X1<§<i3=().(}

MW 425 I=11N‘ ' “

”'15-‘19“!

VALHt=HK(I;J)*XH(i)*7U(JT*SCALF**2

Xi"X}‘?.\"<i1\XL)+\/I:\LUE

KKVH:V.B*XKXH

LHITr(B,1001)

FHerT(ZéHrCHfiWG¥
(

1:1,?)
1J)’J:1-9?~5)91319H)

Hi 425‘)

CH“ f-i“/'A,“L f". Tt‘

f’lkllJL
H (Lyn 1:2(1,gl
(“@le r1\"i}ij(1 ,1")

LxLL MYHUt22,l)
{ALL MYLMJ(2(\,U}
CflLL HYUU(27,H)

LAIL :YULY(UU)
flit IHWN3(??,“)
ALI i‘YUU(21y1)

F—
r" .

HYUU‘QFJ 1 1)

HYHi(lUyK)

“)427,4?7,A2x

,L HleK
l‘HYAIR110,1,‘(3
HYAIR(15,2,X(/
HYUMK

KAT-=0.“

HM A2v 1:1,N

’4 ('3? J31,”

rIWflL=HK(I,J)*X(I)*X(J)*§CALF**2

‘A’S’l 7' '1” :XKXTFF Il‘xs/‘l-

XKXT=U.U*XKXT

kl. VT J=>{:’\XI':-7\r\>( '1‘

U“Iit(3,1002)

.‘
f- -
I

I

))
),X(

“J n (.3 r—d a“. f“ (‘3 r)
. l .“
r-\

3))

m7:
I

‘a>-.':.-. “Hm CHE/{J}

a": I ‘1 t ( 5,500 ) KKXIVI ,mix‘: ,(Lt'S TJ
xmilm(:,1c03)

gm HA1(20H CHAT3$ L1H VNAHLE TU TRZI)

3») Al i Ski

:*Lf1'ft:(ii, 8(‘(3) (SJ (1 ) , l=—] ,f“)
FHRHAT(5E2O.4)

GU Th 30'"'
HHITt(3,900)H(2,L)
FQPVATtEZG.4)

CH TJ no
wHITE(3,1%O)I,P(I.1)
PUP-511W ( i Z ,10X,'F}’(;. (a)
CH 80 '
WRITE(?,?5O)H,”(W+1,1)
unPrHVr(I2,]rn<,fi?fl.re)
(We TH 10
Pk77f(3,350)9,§(“9“)
g-‘Hm‘. {.1 T ( ‘i X 1 5 ("EX 7 ‘9 7’ . ’4.)
C:LL IXIT

l'. ’2!
l 1‘

T:‘
zU