A GRADIENT COMPUTATIONAL TECHNIQUE -
FER A CLASS OF OPTTMAL CONTROL
PROBLEMS SUBJECT TO
ENEQUALTTY CONSTRAINTS

Timis fer the Degree cf Ph. D.
MICHEQAN STATE UNIVERSTTY
SHYH- JONG WANG
1959

 

‘LFYLT'Ib

LIBRARY '

Michigan State
University

 

This is to certify that the

thesis entitled

A GRADIENT COMPUTATIONAL TECHNIQUE
FOR A CLASS OF OPTIMAL CONTROL PROBLEMS
SUBJECT TO INEQUALITY CONSTRAINTS

presented by

Shyh Jong Wang

has been accepted towards fulfillment
of the requirements for

Ph. D. degree in E. E.

/”
‘ ’/ ‘ 6:" A1]
2 /; Major professor
Date 9/ 18169

0-169

m "Effie -w

1am? '
-. NS' '
5:, not must? we ‘

ABSTRACT

A GRADIENT COMPUTATIONAL TECHNIQUE
FOR A CLASS OF OPTIMAL CONTROL PROBLEMS SUBJECT TO
INEQUALITY CONSTRAINTS

BF

Shyh J ong Wang

A broad class of Optimal control problems with
isoperimetric constraints and instantaneous algebraic cons-
traints, the control problem of Bolza, are considered. An
important subclass of the general control problem of Bolza
which contains the bang-bang control problems and problems
with continuous control variables as well as discrete
control variables are also considered.

Local linearization and perturbation techniques
are used to obtain the computational algorithms for the
solutions of these problems. Iterative procedures for
these algorithms are given in detail. A sample program
for the general algorithm written in FORTRAN is presented
in APPENDIX A.

ISOperimetric inequality constraints are trans-
formed into equality constraints by introducing additional
control parameters, and a penalty function technique is

used for treating the instantaneous algebraic constraints.

Shyh J ong Wnag

‘Ihree numerical examples are given to illustrate
the application of the optimization process, and to demon-
strate the influences of the alternative choices of the
iterative parameters and their reSpective updating schemes.
The solutions of these examples are presented in curves

and in tabular forms, and the results are discussed.

A GRADIENT COMPUTATIONAL TECHNIQUE
FOR A CLASS OF OPTIMAL CONTROL PROBLEMS SUBJECT TO
INEQUALITY CONSTRAINTS

Shyh Jong Wang

-‘ m-

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 System Science

1969

1 .J

K.)

TO MY PARENTS

ii

ACKNOWLEDGEMENTS

The author wishes to express his sincere appre-
ciation to his adviser, Dr. Edgar C. Tacker, for his
constant support, motivation, and guidance throughout
this work. The author also wishes to thank the other
members of his Guidance Committee, Dr. Harry G. Hedges,
Dr. John B. Kreer, Dr. Richard C. Dubes, and
Dr. J. Sutherland Frame, for their guidance of his doctoral
program.and for their helpful suggestions during the final
preparation of this thesis. Thanks are also due to
Mr. Robert J. Greiner for his help in programming and for
many important suggestions.

The author gratefully acknowledges the support
provided to his work by the Division of Engineering
Research. Use of the Michigan State University computing
facilities was made possible through support, in part,
from the National Science Foundation.

Finally, the author wishes to thank his wife,
Nancy, for her patience, understanding, and encouragement

throughout his graduate work.

iii

TABLE OF CONTENTS

Page
ABSTRACT

ACKNOWIJEmmENTS O O O O O O O O O O O O O O O O O O 0 iii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . vii

LIST OF TABLES O O O O O O O O O O O O O O O O 0 O O 0 ix
LIST OF APPENDICES O C O O O O O O O O C O O 0 O O O O x
I INTRODUCTION . . . . . . . . . . . . . . . . . . . 1

II THE CLASSICAL VARIATIONAL THEORY AND

OPTIMAL CONTROL PROBLEMS . . . . . . . . . . . . . h

2.1 The First Variation of a Functional . . . . A

2.2 Extremum.of a Functional; A Necessary
Condition0.000000000000000 6

2.3 The Second Variation of a Functional;

A Necessary Condition Involving the

Second Variation . . . . . . . . . . . . . . 8
2.h. The Euler-Lagrange Equation . . . . . . . . 9
2.5 .A Generalization of Euler-Lagrange

Equation to the n-Dimensional Fixed

End-Point Problem. . . . . . . . . . . . . . 10
2.6 Weierstrass Necessary Condition (For

a Strong Extremum) . . . . . . . . . . . . . 13
2.7 Legendre's Necessary Condition . . . . . . . 15
2.8 Corner Conditions . . . . . . . . . . . . . 15

2.9 Jacobi's Necessary Condition; Conjugate
POints O O O O O O O O O O O O O O O O 0 O O 16

iv

2.10

2.11
2.12

2.13

2.111.
2.15

Necessary Conditions; Sufficient
Conditions . . . . . . . . . . . . . . .
Variable End-Point Problem . . . . . . .
A General Problem.of the Calculus of
Variations——the Problem of Bolza . . . .
The Multiplier Rule for the Problem

of Bolza . I . . . . . . . . . . . . . .
The Control Problem of Bolza . . . . . .
Summary . . . . . . . . . . . . . . . .

III A GRADIENT COMPUTATIONAL TECHNIQUE FOR THE
CONTROL PROBLEM OF BOLZA, ALGORITHM I . . . .

3.1

3-2
3-3

3.h

3-7

3.8

A Formulation of the Control Problem

of Bolza . . . . . . . . . . . . . . . .
A Simple Transformation . . . . . . . .
Derivation of Variational Formulas

for some Quantities of Interest . . . .
Constrained Effort and.Successive
Optimization,.Algorithm.I . . . . . . .
An Iterative Procedure for.Algorithm I .
Some Measures for the Validity of
Linear Approximation . . . . . . . . . .
The Instantaneous Algebraic Constraints
and the Penalty Function Technique . . .
Conclusion . . . . . . . . . . . . . . .

IV A SUBCLASS OF THE CONTROL PROBLEM OF BOLZA,
AMORITIM II 0 O O O O C O O O O O O O O O O 0

L1
ho2

ho3

Formulation of the Problem . . . . . . .
A Multiple-stage Formulation of the
Problem. . . . . . . . . . . . . . . . .
The Variations of the Multiple-stage
Problem.with.Fixed.Number of Stages . .

Page

17
18

21
23

26
33

36

37

52
S9

60

61
65

67

67

7O

73

Page

halt Successive Optimization Process for the

Problem with Fixed Number of Stages . . . . 76
ins A Computational Procedure for the
Problem with Fixed Number of Stages . . . . 78

Lt.6 Optimization Process for Adding a New
Stage and the Corresponding Computational

Procedure................. 80
14.7 A Computational Procedure for

AlgorithmII................ 85
vmmmICALmCAMPLEs................ 87

5.1 Brachistochrone Problem with Inequality
StageConstraint.............. 37

5.2 Orbit Transfer of a Solar Sail Ship . . . . 100

5.3 Low Thrust Trajectory Optimization
Problem.................. 110

VISUMMARYANDCONCLUSION.............. 121

REWCES O O O O O O O O O O O O O O O O O O O O O 125

vi

Figure

10
11

12

13

LIST OF FIGURES

The arcs x* and.x . . . . . . . . . . . .
Neighboring curves for variable endppoint
problem . . . . . . . . . . . . . . . . .
A nominal solution arc and its neighboring
arc with variable end-points . . . . . .
Brachistochrone problem . . . . . . . . .
Paths of brachistochrone problem - Case I
tf vs. N , brachistochrone - Case I . . .
tf vs. N , brachistochrone - Case II . .
tf vs. N , brachistochrone - Case IV . .
Control functions, brachistochrone - Case
Paths of brachistochrone problem - Case IV
Orbit transfer of a solar sail ship . . .
Trajectories of solar sail ship - Case II

Sail angle of solar sail problem - Case II

Orbit transfer of a low thrust rocket . .

vii

IV

Page

19

89
95
96
96
98
98
99
102
108

109
111

Figure Page

15 Trajectories of low thrust rocket Case I . . 116

16 Thrust angle of low thrust rocket Case I . . 117

17 Trajectories of low thrust rocket Case II . . 119

Case II . . 120

18 Thrust angle of low thrust rocket

viii

Table

LIST OF

Numerical results,
Numerical results,
Numerical results,

Numerical results,

TABLES

Page

solar sail Case I . . . . 106

solar sail Case II . . . . 10?

low thrust Case I . . . . 115

low thrust Case II . . . . 118

ix

LIST OF APPENDICES

Page

Appendix

A A Sample Program . . . . . . . . . . . . . . 129

I INTRODUCTION

The behavior of a physical system can usually be
described mathematically in terms of the parameters of the
system. One may call this representation the mathematical
model of the system. The behavior of a large class of such
physical systems can be represented or approximated by a
set of simultaneous ordinary differential equations of the
first order, and the scepe of this thesis is limited to
this class.

Once a mathematical model is given or identified,
one can, at least in principle, choose the parameters, or
the control variables of the system such that the perfor-
mance of the system is optimized in a sense which is
specified beforehand. The object or the functional which
is to be extremized subject to the mathematical constraints
of the system is called the performance index, or the cost
functional of the system.

A complete understanding and mastering of the
problem is not a simple matter; it has been the major work
of a branch of applied mathematics, the calculus of varia-
tions, for almost three hundred years.

The calculus of variations has been studied and
developed by many mathematicians, necessary and sufficient

conditions for many problems have been developed. However,

only relatively few problems can be solved using these
conditions directly.

The primary interest Of engineers and scientists
is perhaps the realization of a given problem, that is to
obtain a specific solution Of the problem. In most cases,
practical Optimization problems are sufficiently complicated
that solutions in compact form can not be obtained. There-
fore numerical techniques have been studied extensively
since high speed computers became available.

_ One of the most general problems in the calculus
Of variations is the sO-called problem of Bolza. This
problem was first formulated by Bolza in 1913 [BL-l], has
stimulated a great deal Of research and study since 1930,
and has become a very important problem in Optimal control
theory. Hestenes [HE-l, -2] formulated an equivalent
problem to the problem Of Bolza which is in a more desirable
form for the study of Optimal control theory, we shall call
this problem the control problem Of Bolza. In Chapter II
Of this thesis, we shall review briefly the work of the
classical variational theory and modern Optimal control
theory develOped in the past three hundred years. This
material will serve as background for the later develOpment
of this thesis.

The contributions of this work are in Chapter III
and Chapter IV. The contribution in Chapter III is the
extension Of the method Of gradients [KE-l, BR-l] to the

control problem Of Bolza with various constraints, and that
in Chapter IV is the formulation of an important subclass
of the control problem Of Bolza, and the derivation Of a
computational algorithm for the solution Of the problem.
Three numerical examples are given in Chapter V
to illustrate the application of the computational algorithm
for the general problem, and to show the effects Of the
iteration parameters on the speed of convergence.
In Chapter VI, we summarize this work and discuss
the possible extensions for further research.

A sample program is given in Appendix A.

II THE CLASSICAL VARIATIONAL THEORY AND
OPTIMAL CONTROL PROBLEMS

The variational calculus has been developing since
the late seventeenth century, having its beginning with
JOhann Bernoulli who posed the brachistochrone problem in
1696. This problem.was solved by the Bernoulli brothers,
Newton, de l'Hospital, and others. In 1697, Johann Bernoulli
solved another well known problem-—the problem Of geodesics,
later L. Euler and J. Lagrange solved the general problem.of
this type.

Besides these two problems there is another problem-—
the classical isOperimetric problem. The general method of
solving this problem was given.by L. Euler. These three
problems have had a great influence on the development Of the
variational calculus.

. For convenience as well as to introduce notation,
we shall now give a brief discussion of the classical varia-
tional theory, indicating the various mathematical arguments
employed and the results obtained. This will lead us to the

consideration of Optimal control problems.

2.1 The First Variation of a Functional

Let X be a normed linear space, whose elements are
real-valued functions defined on a closed interval I of R,
where R is the real line. ‘A functional J(-) is defined as

u

S

a mapping which maps X into R. We shall particularly be
interested in the functional

.151 ,

J(x) = ./t F(t,x,x)dt (2.1.1)
~ O
. A dx . . .

where x - d't' , and F is a real-valued function of its
arguments. The space X will be considered to be endowed
with the norm (H .H) defined by

Hxll = sup HUSH (2.1.2)
#561

where I = [to , 1, and x E X. The notion Of distance

t:L
between two elements x, y of the space X can be defined as

the norm of their difference (1.6., the distance induced by
11.“),
d(x,y) = iix - y“ (2.1.3)

Give a fixed function x(t) and its increment €(t)
of the space X, the corresponding increment of the functional

J(x) is a functional of g,

AJ(€) = J(x + 6) - J(x) (2.1.14)
Suppose we can write

AJ(€) = 1(5) + €ii€il (2.1.5)

where l(€) is a linear functional of §, and 3-DO as lléll-DO.

The functional J is differentiable if it has the above

 

property. 1(5) is called the first variation of J along x, and

 

we shall denote it by Jl(x,€). One can show that J'(x,€)
is unique (GE-1).

2.2 Extremmn of a Functional; A Necessary Condition

Let X1 be a normed linear Space of all real-valued

continuous functions on I having piecewise continuous
derivatives. It is clear that an element ole is also an
element of X. The elements of X.l will be called admissible

functions. It is convenient for the subsequent discussion

 

to define another norm, H '“l’

“x“l = “x“ + “x“ -"-’- sup lx(t)l + sup but“ (2.2.1)
{:61 t€I

Let x be an arbitrary but fixed element of X1. The strong
neighborhood of x corresponding to 6 > O, Ns(x,€), is the

set of all functions x in X such that “x - x" < 8, while

1
the weak neighborhood of 2': with e > o, Nw(x,8), is the set

of all functions x in x such that “x .. SEMl < 8. It is

1
clear that for a given 8 and 3%, the strong neighborhood
contains the weak neighborhood. We shall say that x = 32
yields a strong relative extremum of J if there is an 8 > 0
such that for all x 7 5’: and x E Nah-5,8), J(x) - J(x) has the
same sign. 0n the other hand, we say that x = 32 yields a

weak relative extremum of J if there exists a weak neigh-

 

borhood Nw(x,€), such that J(x) - J(x) has the same sign

for all x 7‘ 5': and x E Nw(x,8). Since “x - EH1 < e

7

implies llx - in s e, if 5% yields an extremum of J with
respect to all x with “x - i“ S 8, then 5: yields an extremum
of J with reapect to all x with “x - illl S 6. Hence it is
clear that if i furnishes a strong extremum of J it also
furnishes a weak extremum. A necessary condition for a weak
extremum is also a necessary condition for a strong extremum.
‘ It is important to note that a relative extremum is

determined by local properties of the functional in question and
it is always associated with a neighborhood of that which
yields the extremum. 0n the other hand, an absolute extremum
is determined by global properties of the functional. In the
subsequent discussions, only relative extrema are considered;
therefore, necessary and sufficient conditions for an
extremum are local although, for simplicity, we may not
mention it explicitly.

The following lemma is basic in the calculus of
variations and is known as the fundamental lemma [HE-LL].

Lemma 2.2.1 Let M(t), N(t) be piecewise continuous func-

tions on I. Then

t .
ftl (M(t)€(t) + N(t)€(t))dt = 0 (2.2.2)
0

holds for all 5 in X1 and €(to) = €(tl) = 0 if and only if

there exists a constant c such that

l/\
C‘-
M
d

t
N(t) = [to M(T)dT + c t0 1 (2.2.3)

Theorem 2.2.1 A Necessary Condition. If J is differen-
tiable, and if i furnishes an extremum.for J, then

J'(i,§) = 0 (2.2.10

for all g in.X1.
The proofs of this theorem.and.those which will be
stated in the subsequent sections can be found in.most of

the standard texts (see, for example, [BL-2, GE-l, HE-h1).

2.3 The Second Variation of a Functional; A Necessary
Condition Involving the Second Variation
we say that the functional J is twice differentiable

if its increment can be expressed as
Mtg) = 1(a) + q(&;) + ensug (2.3.1)

where 1(5) is the first variation of J, q(g) is a quadratic
functional, and e—po as Nan—>0. q(g) is called the
second variation of J along x. The second variation of J

is also unique and is denoted by J"(x,€).

Theorem 2.3.1 A Necessary Condition. Let J(x) be twice

 

differentiable, then x = i yielding a minimum.(maximum) of

J implies that
J"(i,€) 2 o (s 0) (2.3.2)

for all 5 in.X1.

2.11. The Euler-Lagrange Equation
Consider all the arcs x in X1 having fixed end
points x(to) = x0 , x(tl) = x1 . We ahall refer to these
arcs as the admissible arcs of fixed end-point. In this
section we are going to state the first necessary condition
of the following problem:
Among all the admissible arcs of fixed end-
point, find the one yielding a weak extremum
for the functional
1:1
J(x) = ft F(t,x,x)dt (2.u.1)
c
We shall assume that F is a known function of (t,x,5c) with
continuous first and second partial derivatives with respect
to all its arguments. This problem is known as the simplest

variational problem.

Theorem 2.1L.l If x = x“. furnishes the functional J(x) an
extremum in the class of all admissible arcs connecting its
end-points x0 and x1, then there exists a constant c such

that

t
F}; = f det + c (2.1+.2)
t
o
__ i:- A 61“
holds at every point of the are x - x , where Fx —- BE , and
Fe 9. é; .
x 5::

Equation (2.1+.2) is known as the integral form of

10
the Eul er-Lagrange equation

Ed? Pi = F (2014-03)

X

An admissible arc satisfying the Euler-Lagrange equation is
call an extremaloid.

Note that the extremum of J in Theorem 2.Lt.l is a
weak extremum. However, as mentioned earlier, a strong
extremum is also a weak extremum and any necessary condition
for a weak extremum is also a necessary condition for a
strong extremum. This implies that Theorem 2.L|..l gives a

necessary condition for both weak and strong extrema.

2.5 A Generalization of Euler-Lagrange Equation to the
n-Dimensional Fixed End-Point Problem

Let x = (JCJ',...,xn)T be an n-vector and each com-
ponent of x be a real-valued function on I. Let X1 be a
linear space of all n-dimensional functions x and each
component of x is continuous on I and has a piecewise con-
tinuous derivative. We shall call members of X1 admissible
arcs. Let F(t,x,x) be a continuous function having conti-
nuous first and second partial derivatives with respect to
all its arguments. The problem is stated as: find the
necessary conditions (Euler-Lagrange) for the functional

J to have a weak extremum furnished by an admissible arc

11

x‘= x* with fixed end-points x(to) = x0 and x(tl) = x1 .

J is defined as

t

1
J(x) = / F(t,x,x)dt (2.5.1)

1;0

Let g = (51,...,€n)T be an increment of an
admissible arc x. ‘We shall call s an admissible incre-
ment of an arc x if g is a.member of X1 and

€(to) = €(t1) = O. The correSponding increment AJ is

J(XHS) - J(X)

AME)

t

1 .
/ (F(t,x+g,5c+g) - F(t,x,x))dt (2.5.2)

t
o

The first variation of AJ can be obtained by expanding
(2.5.2) using Taylor's Theorem,

t
1 .
J'(x.€) = / mg. + Pinata (2.5.3)

to

where FI and Pi are the row vectors of partial derivatives

12

with respect to x and x reSpectively. We can write

%=(F1uu.Ffl (agm
X X

Pi: (F.1 eooee F.n) (20505)
X X

In (2.5.3) choosing all the components of {5. except the ith

one to be zero, we have

,t . .
J'(x,€) = Jtl (F 1&1 + F.i€i)dt i = 1,...,n (2.5.6)
o x x
Suppose x* extremizes J among the admissible arcs connecting
its end-points, by Theorem 2.2.1 the eXpression (2.5.6) must
be zero for x = x*, then by Lemma 2.2.1 there exists a

constant ci such that the equations
dt + Ci 1 = 1,...,n (205.7)

hold along x*. We thus have the following theorem:

Theorem 2.5.1 If the are x* extremizes J among all the

admissible arcs connecting its end-points, then there exist

:1 constants ci such that

F. ./t F 1 dt + c i = 1,...,n
X

hold along x*.

13

Equation (2.5.?) is the integral form of the Euler-

Lagrange equations

F = F i i = 1,...,n (20508)
I

2.6 Weierstrass Necessary Condition (For a Strong
Extremm)
Let x = x35 be a (strong) minimizing arc for J
among all admissible arcs with fixed end-points. The
arcs considered are n-dimensional. The Weierstrass

E-function is defined as

E(t,x,x,X) = F(t,x,X) - F(t,x,x) - Fi(t,x,x)(X-x)
(2.6.1)
In (2.6.1) (t,x,x) is a point on the arc x* at time t
and having derivative 5:. Similarly, (t,x,X) is a point
on an (neighboring) admissible are x (see Figure 1).
The following theorem is due to Weierstrass,

Elleorem 2.6.1 Weierstrass Necessary Condition. If
the arc in furnishes a strong minimum for J among the

class of admissible arcs connecting its end-points, then
E(t,x,:’c,x) 2 0 (2.6.2)

1'll-celds for every (t,x,x,X) with (t,x,5c) on 26:. and (t,x,X)
0-1:: an admissible arc and X 9‘ x.

 

11L

 

 

 

 

 

X1
X /
ii-
I
1
t2 t2+rd t2+d t1

Figure 1 The arcs x“-

and X

15

Remark 1: The condition (2.6.2) is trivially satisfied
ifX=2°c.

Remark 2: If x* (strongly) maximizes J the inequality
sign in.(2.6.2) should be reversed.

2.7 Legendre's Necessary Condition

Theorem 2.7.1 If the arc x* (weakly) minimizes the
functional J in the class of admissible arcs connecting
its end-points, then

aTFii(t,x,x)a 2 0 (2.7.1)

along x* for each non-zero vector a = (al,...,ap)T. Where

 

 

PF 0 O O O O F. . q
x1 1 1 n
Fee = .
F O O O O O F. .
xnxl n :1.4

This condition is a consequence of the weierstrass
necessary condition, although Legendre proved it very
differently and at a.much.earlier time.

2.8 Corner Conditions
The following is known as the Weierstrass-Erdmann

corner conditions,

16

meorem 2.8.1 Let x* be an admissible arc satisfying
the necessary conditions (2.5.7), the Euler-Lagrange
equations, then (I) Fi(t,x,x) is continuous along x*,
and (II) If x* also statisfies the Weierstrass necessary

condition then F - Fix is continuous along x*.

2.9 Jacobi's Necessary Condition; Conjugate Points

An admissible are x* without corners vmich is a
solution of (2.5.7) is called an extremal. Any solution
of (2.5.7) consists of a finite number of extremals. An
admissible are 1: along which the determinant of F3.“-c is

non-zero is said to be nonsingular.

An extremal which is also nonsingular is called
a nonsingular extremal. From the hypothesis on F and the
above theorem, a nonsingular extremal has a continuous
second derivative.

We say that a point (t2,x(t2)), to< taé t1 , on
the extremal x* is conjugate to (to,x(to)), if there is
an accessory extremal 5 such that €(to)= €(t2) = 0 but

€(t)¢0 on to< t<t2.

Theorem 2.9.1 Jacobi's Necessary Condition. If x65 is
a nonsingular minimizing are without corners for the

functional J (2:) defined by equation (2.5.1) in the class
of admissible arcs connecting its end-points, then there
does not exist any point (t2,x(t2)) with to< t2< t on x*

l
conjugate to the point (to,x(to)).

.17

2.10 Necessary Conditions; Sufficient Conditions
I. Necessary Conditions for a Weak Extremum: If the
admissible are x* on I provides the functional

t
J(x) = ftl F(t,x,x)dt (2.10.1)
O

a weak (relative) extremum among the class of admissible
arcs connecting its end-points x0 and x1, then
(Ia) The Euler-Lagrange equations are satisfied
by in (Theorem 2.5.1).
(lb) The Legendre's necessary condition holds
along x* (Theorem 2.7.1).
(Ic) The Jacobi's necessary condition holds
(Theorem 2.9.1).
II. Necessa_ry Conditions for StrongEbctremm: If x“. on I
provides a strong (relative) extremum for the functional
(2.10.1) among the class of admissible arcs connecting its
end-points, then the conditions (Ia), (Ib), (Ic) and (Ila)
are satisfied , where (IIa) is
(Ila) The Weierstrass necessary condition (Theorem
2.6.1).
III. Sufficient Conditions for a Weak Extremmn: If for
an admissible are x‘”, the following conditions hold simul-
taneously, then it.“ provides a weak (relative) minimum for
the functional (2.10.1):
(IIIa) The arc x* is an extremal, that is, x*

satisfies (Ia) and x“. has no corners.

18

(IIIb) The strengthened Legendre's necessary condi-
tion, that is, the inequality of (2.7.1) is
replaced by strict inequality.

(IIIc) The strengthened Jacobi's necessary condition,
that is, there is no point in (t°,tl] on x*
conjugate to (to,x(to)).

IV. Sufficient Conditions for a Strong Minimum: If the

admissible arc satisfies (IIIa), (IIIb), (IIIc) and (IVa)

then :6“ provides a strong (relative) minimum for the functional

(2.10.1) with fixed end-points , where (IVa) is

(IVa) There exists a neighborhood of the elements
(t,x*,x*) on x* such that the strict
inequality E(t,x,x,X) > 0 holds for all
(t,x,x,X) with (t,x,x) in the neighborhood,
(t,x,X) admissible, and X )6 5:.

For further discussion of sufficient conditions,
see [HE-)1, GE-l, BL-2].

2.11 Variable End-point Problem

Let x and y be admissible arcs defined on I = [to,t1]
and I' = [t5,t'] respectively, with end-points x0 = x(to),
x1 = x(tl), yo = y(t('>), yl = y(ti). Extending x and y to
some interval containing IUI'. The distance between the arcs

x and y is defined as

d(x " n 1 1
i=1

 

19

 

 

 

 

 

Figure 2 Neighboring curves for variable
end-point problem

20
where

d(xi.yi) = 83p lxi(t) - yi(t)| + 8:p lii(t) - §i(t)l

 

+ lxi(to) - vinyl + kiwi) - Sri(ti)l

(2.11.2)

Let x and y be neighboring arcs, let g be their

difference

€(t) = y(t) — x(t) (2.11.3)

and write the endppoints of y as

+
y (to + dto , xO dxo)

O

y1

(t1 + dt1 ’ "1'+ dxl)

The functional J(x) is defined as

t
J(x) =/t1 F(t,x,x)dt (2.11.h)
0

where F is a known function.having continuous partial
derivatives up to the second order with respect to all its
arguments.

The first variation of J is defined as the linear
part of the increment AJ(€) with respect to g, g, dto, dxO

and dxl, and is denoted by J'(x,€). We state the following

theorem,

21

Theorem 2.11.1 If the arc x”" minimizes the functional J
on the class of admissible arcs with variable end-points,
then x“. satisfies the necessary conditions stated in

Section 2.10, and the boundary condition

1:1

[(F - Fix)dt + Fidx] = 0 (2.11.5)

t
o

The boundary condition (2.11.5) is called the

transversality condition.

2.12 A General Problem of the Calculus of Variations—
the Problem of Bolza
The problem of Bolza can be stated as follows:
Find in the class of arcs x satisfying the equality

constraints
”8(t9x9i) = 0 S = 1,000’q < n (2.1201)
and the boundary constraints

9Y<t°.x(to).tf.X(tf)) = 0 Y = 1.--o.p S (2n+'1)

(2.12.2)

the one which minimizes the functional

22

t
J(x) = g(to,x(to),tf,x(tf))-+./tf L(t,x,x)dt (2.12.3)
o ,
where as, BY, g, and L are real-valued functions of their
arguments. The continuity properties of these functions will
be given in Section 2.13.
This form of the problem of Bolza was formulated by
Bliss [BL-3]. There are two other equally general problems,
the problem of Mayor, and the problem of Lagrange. These two
latter problems can be obtained from the problem of Bolza by
setting, reapectively, L E O and g .=_ 0 in (2.12.3). On the
other hand, one can obtain the problem of Bolza from the problem
of Mayer and that of Lagrange through suitable transformations
as introduced below: Consider the n+1 dimensional are x, an'
defined on I = [to’tf] satisfying the constraints

¢s(t,x,x) = O s = 1,...,q

in+1 - L(t,x,x) = 0 (2.12.4)
mnd

97(to.x(to),tf,x(tf)) = 0 Y = 1,...,p

x9*1(to) = 0 (2.12.5)

With this transformation, the functional J becomes

J(x) = g(to,x(to),tf,x(tf)) + xF*1(tf) , (2.12.6)

23

which defines the problem of Mayer.

If we consider the class of arcs x, x”1 on I
satisfying the constraints
3 - _
g (t,X,X) " O s = 1,...,q
inn“) = 0 (2.12.7)
and
eY(to,x(to),tf,x(tr)) = 0 Y = 1,...,p
Flog.) - FEET-3(to,x(to),tf,x(tf)) = 0 (2.12.8)
o
The functional J becomes
t
f . n+1
J(x) = / [L(t,x,x) + x (t)]dt (2.12.9)
t
o

This is the problem of Lagrange. Hence we have proved the

equivalence of the problems of Bolza, Mayer, and Lagrange.

2.13 The Multiplier Rule for the Problem of Bolza
Let x be an n-dimensional point which represents the

"state" of an arc x at time t E I, thus (t,x,x) is a

(2n+1)-tuple. Let R2nl-l be an open region (arcwise connected

open set) of elements (t,x,x), and let “s and L be continuous

on R2n+l and have continuous first partial derivatives.

Assume that the matrix ”3': has its maximum rank q on the

2h

subset R0 of R defined by the constraints (2.12.1).

21ml Zml
Note that ¢ is the vector notation for $3, 8 = 1,...,q.

Let x(to) = x , x(tf) = xf, and let S be an Open region

0 2n+2
of (2n+2)4tuples (to’xo’tf’xf) on which.6Y, and g are

continuous and have continuous first partial derivatives,

and the matrix (9 9 9 G ) has its maximum rank p
to xo tr x,f

Points in R

o
2nrl°

continuous and has piecewise continuous derivatives is

at each point of S are

2n+1’ $2n+2
An arc x which is

2n+2‘
admissible if they are also in.R

admissible if its elements (t,x,x) and (to,xo,tf,xf) are
admissible. With these assumptions we can state the
multiplier rule for the problem of Bolza formulated in the

last section [BL-2, HEph],

Theorem.2.13.l The Multiplier Rule. Let x* be an
admissible arc defined on [t°,tf], which.minimizes the
functional J, then there exist multipliers 1°, 18(t),

s = l,..,q«<n, and by, y = 1,...,p such that (1) 1° and h?
are constants and are not all zero; (2) 13(t) are continuous

on [to,tf] except possibly at those t corresponding to
corners of x* where the left and right limits exist, and 1°,

and 13(t) do not vanish at the same time on [to,tf]; (3) If

we define
F(t,x,$c,x) = x°L(t,x,:'c) + AT¢(t,x,5c) (2.13.1)

then

25

t
11.1 = / F.1dT + c1 i = 1,...,n (2.13.2)
x to x

¢s(t9x,;t) = O S = 1,...,q (201303)

hold along x* , where cl are constants.

(h) The conditions

t
f
[(F .. Fix)dt + Fidx] + A°dg + the = o (2.13-LL)
to
and
eY(to,xo,tf,xf) = O Y = 1,...,p (2.13.5)

hold at the end-points of x* for every choice of dto, dtf’

dxo and dx And furthermore (5) The Weierstrass condition

1..
E(t,x,5c,x,1) 2 0 (2.13.6)

hold for every (t,x,3c) on x* and for each “5:393” in 32ml

such that X 3‘ i, where

E(t,x,x,X,7\) = F(t,x,x,).) - F(t,x,i,1) - Fi(t,x,x,}.)(X - 5c)
(2.13-7)

It is easy to see, from (2.13.2), that along x“.
F3: is continuous and has piecewise continuous derivatives,

and

HE Pi = FX (201308)

26
‘JZhis is the Euler-Lagrange equation (in row-vector form).

2.111 The Control Problem of Bolza

We shall formulate a very important class of Optimal
control problems which have their origin in the classical
problem of Bolza. We shall call the new problem the control
problem of Bolza. Later in this section, we shall show that
these two problems of Bolza are equivalent.

Let the cost functional J be defined as

t

r
J = g(tf,x(tf),a) + / L(t,x,u,a)dt (2.1.11.1)

tn’o

where to is the initial time and t is the final or terminal

1'
time, which can either be fixed or variable; x is an n-vector,
(x1,...,xn)T, called the state variable (vector); u is an
m-vector, (u1,...,um)T, called the control variable (vector),
and a 18 an r-vector, (0-1,...,0-r)T, called the control
parameter (vector). We shall also use

x : t,x(t),u(t),a,t E [to,tf]

to represent the solution arc of the differential constraints
5: = f(t,x,u,a.) (2.11.1.2)

Equation (2.1h.2) may be considered as the mathematical model
of a system which is nonlinear in general. Let E be a region
in an (n+m+r|-l)-dimensiona1 Euclidean space whose elements

27

are the (n+m+r+1)-tup1es (t,x,u,c). We assume that the
real-valued functions f1, 1 = 1,...,n, are continuous and
have continuous first partial derivatives on E. Let Eo be

a. subset of E defined by the side constraints

¢8(t9x9u9a) s o s = 1,ooo,q' (2.11-L038.)

08(t,x'u,a) = O S = Q'+l’eeo,q (20m03b)

therefore the elements of EO are those in E that also satisfy
(2.11.)..3). Assume that each as is continuous and has con-
tinuous first partial derivatives on E, and the
qx(m+q)-matrix [flu I91] has its maximum rank q at each point
of E0, where {D = (¢1,...,¢q)T, ”u = 3% , and I is the qxq
identity matrix. Elements of Eo are called admissible
elements. We shall assume that the control functions u(t),

t E [to,tf] are piecewise continuous. A solution are x of
(2.11.1.2) with x(to) in a given initial manifold is called an
admissible arc if all its elements (t,x,u,a.) are admissible.

Consider the isoperimetric constraints

IY s O Y = 1,...,p. (Balk-(+3)
IY = 0 Y = p'+l,eee,p (20114-01413)
where
1‘1‘:
1* = 6Y(tr,x(tf),c) + / MY(t,x,u,c)dt (2.1LuS)
t

O

28

We shall use I, 9, and M to represent the vectors (11,...,Ip)T,
(91,...,9p)T, and (M1,...,MP)T resPectively. Assume that
97, and M‘{ are continuous and have continuous first partial
derivatives on E. Note that terminal constraints are Special
cases of isOperimetric constraints and hence are included in
(2.11-1.11).

Now we state the control problem of Bolza as follows:

Find in the class of admissible arcs
x : t,x(t),u(t),a. t G [to,tf]

satisfying the constraints (2.1L|..2), (2.11.1.3) and (2.11.1.11),
the one which minimizes the cost function J.

A problem, perhaps somewhat less general than the
one given here, was first formulated by Hestenes in 19119
[HE-1,2,3]. Since the new formulation and the classical
formulation of the problem of Bolza are equivalent, necessary
conditions for the new problem can be obtained from the
classical one through translation [HE-l, BE-l]. Necessary
conditions for this problem have been obtained directly by
Pontryagin, and the result is called the maximum principle,
see [PO-l]. The methods of McShane [MOS-l] and Pontryagin
were generalized by Hestenes, who obtained more general
results [HE-3], and the results were also obtained by Guinn
With weaker hypotheses [GU-l]. We shall state the first
necessary condition for the control problem of Bolza obtained

by Hestenes.

29

Setting
._ 0 T T

F(t,x,u,a,l,h) "' A L + X g + h M (2.11-L06)
and

H(t,x,u,c,p,l,h) = pr - F (2.1J+.7)
“where

T T
A = (llyoooglq), h = (ha-900091113): PT = (pl,""pn)
(2.1.11.9)

we have the following:
Theorem 2.11.91 Let x45 be a solution to the control problem

 

of Bolza, then there exist multipliers 1°, 13(t), pi(t), and
3hY such that (1) 1° 2 O, hl,...,hp are constants and for

v= 1,...,p', hYao if IY= o, hY = OifI <0 along x*,

Y
Furthermore 1° , hY , pi(t) do not vanish simultaneously
at any t E [to,tr]. (2) 11(t),...,1q(t) are continuous
except possibly at the discontinuities of u*(t) at which
the left and right limits exist. For s = 1,...,q',
as“) a o it ¢8= o and= o if as <0 along x”. (3)
P1(t),...,pn(t) are continuous and along the arc x* H is

continuous and

3LT = Hp {F = - Hx (Euler-Lagrange Equations) (2.111.10)

30
H=o dH=H (211111)
u 35 t ° '

hold along x*.
()1) The Weierstrass condition

'11: 41' 4!- *-
H(t,x (t),u,c ,p(t),o,h) s H(t,x*(t),u (t),a. ,p(t),o,h)
(2411.12)
holds for all t 6 [to,tf], and for all u with
(t,x*(t),u,e*) e E0 , that is, for all u satisfying the

(3 ans traints

s8(t,x*(t),u,a*) s o s = 1,...,q'
98(t,x*(t),u,c*) = 0 s = q'+1,...,q
(5) The transversality condition
T tr 1‘1'
(B + [-Hdt + p dx] - / (Hnda)dt a o
t t

O 0

holds on x* for every choice of dto, dtf, dxo, dxf and do.
To show the equivalence of the control problem of
Bolza and the classical problem of Bolza, consider the case

with equality constraints only, 1.9.,

¢8(t,x,u,a) = 0 3 = 1,...,q (2.1-(+013)

Y Y tr Y
I = e (tf,x(tf),a) + /t M (t,x,u,o.)dt y = 1,...,p

o
(2.111.111)

31

Let y1 = x1, 1 = 1,...,n; 53’” = M7, me(to) = o,

‘r = 1,...,p. And replace u3 by WW3, J= 1,...,m, and for
definiteness, set WP” = o, for j = 1,...,m. Substituting
these new variables into equations (2.1L)..1), (2.11.1.2).

(2.1mm) and (2.111.111), we have

t
J = 3(tf3y1,...,yn;a) + f: L(t:v1,....yn:§n+p+1.nuirmpmzfldt
° (2.1h.15)
51 - f1(t;y1,...,yn;yn+p+1,...,ymp+m;a) = 0 i = 1,...,n
(2.1h.16)
¢s(t;y1,...,yn;?fip+l,...,yn+p+m;a) = O s = 1,...,q
(2.1u.17)
if = emonluo),...,yfi(to);tf;y1(tf),...,y“(tf)) + WWI.)
Y = 1,...,p
(2.1h.18)

If we set y = (y1,...,yn+p+m)T, and

c1 = 52 - f1(t,y,y,o.) i = 1,...,n (2.1h.19a)
ng = W7 - MY(t,y,y,a) Y = 1,...,p (2.111.19b)
€2n+p+s == ¢8(t,y,y,a) s = 1,...,q (2.111.190)
and

aY = eY(tO,y(to),tf,y(tr)) + me(tf) y = 1,...,p (2.1h.20a)

32

op+3 = ym‘1(to) .1 = l,...,p+m (2.m.20b)

‘nlen, from equations (2.111.15) to (2.111.20), we have the
following,

és(t,y,§hya) = 0 3 = 1,...,mp‘l'q (ZeMeZl)

6*(to.y(to).tf.y(tf)) = 0

Y 1."..2mm (2.114.22)
tr
J = 8(tf.y(tf),a) + ft L(t,y,y,a)dt

O

(2.14.23)

We thus have showed that the control problem of Bolza can

be expressed in the same form as that of classical problem

of Bolza. Similarly, one can show the reverse. Note that

one may have to use the condition that 9):} has maximum rank
1 i

so that at least one set, x 1,...,x q, can be chosen from

1:1,.."1:n such that (2.12.1) can be written as

i i i i
X r g f r(t;x lgeaagx qzzlgeeepzm3ulgaeagum)

where m s n- q,and uj, j= 1,...,m are the rest of the

variables among x ,...,xn, and z‘1 satisfies the constraints

J: 1,000,171

1'

more x , r = 1,...,q are the state variables, and

113, J = 1,...,n, the control variables.

33

2.15 Summary

The variational calculus and the optimal control
theory is a very extensive subject which has proven to be
increasingly important to modern technology. In this chapter
we have briefly reviewed the history and the develOpment of
the subject, from the classical treatment of the so called
simplest problem of variational calculus to a reasonably
sophisticated modern control problem. For further generali-
zation of the subject, see [BE-2, GU-l, HE-h, PO-l, DR-l,
BR-2].

In the next chapter we shall develop a computational
method to obtain numerical solutions for the control problem

of Bolza.

III A GRADIENT COMPUTATIONAL TEDHNIQUE
FOR THE CONTROL PROBLEM OF
BOLZA, AIGORITHM I

In this chapter, we shall first restate the control
problem Of Bolza in a form which will be more convenient for
computational purposes. A technique. similar to that used by
Valentine [VA-l] is used to convert the inequality isoperimetric
constraints into equality constraints. The instantaneous
side constraints are not treated directly; instead, the
penalty function technique is used to treat the constraints
in a indirect fashion.

Variations of various quantities of the control
problem of Bolza are derived in terms of adjoint variables,
or influence functions, and the variations of control
function, control parameters, and initial time and state.

A technique due to Bryson and Denham [BR-1] is then employed
to obtain the Optimal variations of various quantities. A
computational procedure is given, and measures for prediction

errors in dJ and d1 are defined.

3.1 A Formulation of the Control Problem of Bolza
Let t, x, u, o. be the independent variable—time,
the n-state vector, m-control vector, and r-control parameter

respectively. Let E be a region in the Euclidean space with

311

35

elements (t,x,u,a). Let f1 (i = 1,...,n), L. No. MY (Y = lac-.13)
be defined on E, and assume that fi, L, 14°, M)r are of 01(3).
Futhemore, let g, 9°, BY (Y = 1,...,p) be functions of
(t,x,a), and of CI for each (t,x,a.) such that (t,x,u,o.) E E.
We shall assume that the control function uj(t) (j = 1,...,m)
is piecewise continuous on [to,tf], where to and tf are the
initial and terminal times respectively. We shall call u(t)
an admissible control function, if it is piecewise continuous,
and the corresponding (t,x,u,a) is in E for each t 6 [to,tf].

' Define the cost functional, the differential constraints,
and the isOperimetric constraints as those defined in Chapter II ,

t

f

J = g(tf,x(tf),a.) + ft L (t,x,u,a)dt (3.1.1)
0

5: = f(t’x,u,a) (3.102)

iY s o Y : 1,...,P' (301033)

iY = 0 Y = p'+l,eee’p (3010313)
t

"Y Y r Y

1 = o (tf,x(tf),a)+ ft M (t,x,u,a)dt (3.1.3c)

O

For computational convenience, we shall define a stOpping

functional,
tf
S(tf) = 9°(tf,x(tf),a) + / M°(t,x,u,c)dt (3.1.11)
t
O

A

where S = 0 defines the value of t and we assume that .

f,

36

s ,5 o for all t e [to,tr].
Let x be a solution arc of the system (3.1.2) which
can be written explicitly as,

x : t,x(t),u(t),c t 6 [to’tf]

We shall call elements Of E admissible elements, and we shall
call x an admissible are if x is a solution arc Of (3.1.2)

 

corresponding to an admissible control, and each element Of
x is in E.

The control problem of Bolza is to determine the
admissible control function u(t), the control parameter a,
and the corresponding admissible arc x, such that the cost
functional J is minimized and the constraints (3.1.3) are
satisfied.

Since the instantaneous constraints (2.1.11.3) are
treated indirectly, the admissible set E0 has been relaxed
to E. This is satisfactory especially when iterative

computation is concerned.

3.2 A Simple Transformation

It is convenient from the computational point of
view to transform the constraints (3.1.3) into equality
constraints at the elqaense Of introducing additional control
parameters.

Let B be the additional parameter vector with

components 81,...,(3P', and let Z be a p-vector whose

37

components are defined as ,

ZY = (BY)2 Y .-_-_ 1,...,p! (3.2.13)

z)r = o y p'+1,...,p (3.2.110)

If we set IY = 17 + 2?, then constraints (3.1.3) are equiva-
lent to

IY = 0 Y = 1,0009p (30202)

The control problem defined in Section 3.1 can be

restated in terms of the new constraints as follows:
Determine the admissible control u(t) , the
control parameters a, [3, and the corresponding
admissible are x such that the cost functional J
is minimized and the constraints (3.2.2) are
satisfied.

3.3 Derivation of Variational Formulas for Some Quantities
of Interest
Let x* be an admissible arc, we shall call it a
nominal solution are,

«a sit-«u- aa- s-zz- it- s «u-
x : t ,x (t ),u (t ),a. ,B t E [t°,tf]

Let s be an arbitrary small positive number, and let d(x,x*)

be a distance measure between the arcs x, and x* defined by

38

can?) = sup It-t*| + max Ixiwo) - x*1(t‘f,)l
t6 1w“ i=1, . . ,n

+ max Ixiuf) - x*1(t§)l
i=1,..,n

+ sup max Iuj(t) - u*J(t*)|
term?“ j=1,..,m

+ max Iak- a*k| + max IBY - (3*YI
k=1geepr F110'2p'
.. a- __ a- s .
where T - [to,tf], and T - [to , tr]. Let x be a neighboring
. «-
arc of x
x 3 t,X(t),U(t),a,B t E [toatf]

such that d(x,x*) < a. Let At, and Ax(t*) denote the

increments
At = t - t* (3.3.1)
Ax(t*) = x(t) - x*(t*) (3.3.2)

Let 5t and 5x(t*) be the variations corresponding to At and

Ax(t*) respectively, then for arbitrarily small a, we can write

At = 5t + o(e) (3.3.3)

me") = 6x(t*) + cos) (3.3.11)

where 0(8) satisfies

39

lim 1‘“ = 0 (3.3.5)
a +0

 

Let Ai(t*)be the increment
Ai(t*) = x(t*) .. x*(t*) = tan?) + 0(8) (3.3.6)

where 552(9)) is the corresponding variation. In a later
develOpment we shall need the relation between 6x(t*) and

6i(t*) which can be Obtained as

x(t) - x*(t*)

Ax(t*)

[x(t*) - x*(t*)1 + [x(t) - x(t*)]

5i(t*) + 93‘— At + 0(a)
dt*)
'1!-

t

'11-

oi(t*) + 95.. 5t + one) (3.3.7)

{-
dt l
17*

In the above derivation we used the facts that 0(8) + 0(8) = 0(8)
a-

and if; = if; 4» 0(8). Comparing (3.3.11) with (3.3.7) and

noting that both ox(t*) and oi(t*) + i*(t*)ot are principal

linear parts of Ax(t*) , we have

ox(t*) = oi(t*) + i*(t*)ot (3.3.8)
«3(- s- A dx*
where x (t ) - 7| . Now we shall derive the variational
_ dt
«-

t

L10

equation for 6x(t*). The state equations (3.1.2) corresponding

to x*, and x are
0-H- '31-
x (t ) = r(t*.x*(t*).u*(t*).a*) 6* 6 (ti. ti)
i(t) = f(t,x,u,a.) t e [to,tf]

Their difference is

 

 

Aim") = f<t.x<t).u(t>,a) - r(t*.x*(t*).u*(t*).a*)
_ 6f of «- af a
- (5%).“- MS "' (5-5).} AX“? ) + (5-5).} All“? )
+ (95 Ac + cm (3 3 9)
& * O 0
where
F l l -
ar 91.
SP axn
(is-i) = :
6r“ at?
3;”: each *
and (g—E) , a (g—E) and (3%) are defined similarly. The symbol
4!- 41- #-

(.)‘_meana that the quantity in the parentheses is evaluated

along the nominal are x*.

If the notions Of the variations 6t, 5x(t*), Ou(t*),
and 50. are introduced in (3.3.9) and the result compared with

Ai(t*) = oi(t*) + o(e) (3.3.10)

we have the variational equation:
- -::- a
Ox(t*) n (3%) 6t + (gé) ox(t*) + (gt—11') 6u(t ) + (35* do
a- * a-

(3.3011)
Let Q be the quantity Of interest, and be of the

general form
tf

Q = G(tr,x(tf),a.) + 2(5) + / F(t,x,u,a.)dt (3.3.12)

1so

Adjoining the differential constraint (3.1.2) to (3.3.12), and
forming the difference of Q along the arcs x and x* , we have

A0, = G(tf,x(t ,c) .. G(t:,x*(t:),a*) + 203) - MAX.)

1.)
tf

+ / [F(t.x(t).u(t).a) + 1T(t)(r(t.x(t).u(t).a)- i(t))1dt
t

0
t*

f I e‘ - 0".
- /, [F(t*.x*(t*).u*(t*).a*> + xT(t*)(r(t*.x*<t").

t ,

o .
u*(t*).a*) - i*<t*))1dt*

(3.3-13)
SO far we have treated the matter in a fairly general way,
from now on we shall focus our attention on our Specific

problem. Letting

cf
II

t: + Ato (3.3.1Aa)

t = t* + At (3.3.1Ab)

f f f

[12

'3!-

t = t for t e [t’f3 + At°,tr] (3.3.1Ac)

and extending the definitions of the arcs x and 25* to

(tint: + Atf]. (3.3.13) becomes,

AQ = G(tf+Atf.X(tffAtf).a) - G(tf.x*(tf).a*) + 2(a) - up")

1:
f
+ / [F'(t,x,u,o.) - F(t,x*,u*,o.*) + AT(f(t,x,u,0.)
t

O

t +At
f 1' T a
+ / [F(t,x,u,a) + A (f(t,x,u,a) - x)]dt
tf
t +At
O O T a
- f [F(t,x,u,a) - A (f(t,x,u,o.) - x)]dt
to (3.3.15)

Note that in (3.3.15) we have, for simplicity, drOpped the is
associated with t. Elcpanding (3.3.15), we have

AQ=.& “+99. At+aG A+a.ZA
(at)Mr I. (89m, x( 1.) ($1Mr a (21?)... a

t
f

+/ 5F A't-i- 6F A-t...a:: ,-
to {(5).} X( ) (35).)? u( ) (80).} a.

+ X A t + —- A t + -— A -A t (it
“'55).“ X( ) (311),, u( ) (5,1)“. <1 3d )1

+ (F)*,tf Atf - (F)*’to Ato + 0(8) (3.3.16)

113

In (3.3.16) we have used the relations

(F)|x,to At(o

(F)* t Ato + 0(a)
’ o

and

At + 0(a)

(Fth At = (F) I.

, f 1‘ *,tr

(3030173)

(IS-3.1713)

If we replace the differences by their correSponding variations

and combine the higher order terms relative to 8 with 0(a)

in (3.3.16), and compare the result with

AQ = 6Q + 0(8)

we have

50. = (353+ F), t étf + (g— i) anti.) - (F),,

f *Otf

t
f
+ / [(KT(§§)*+ (g— 312*) )6»: - 1 Todet

O
t

f
T of
+/t [x (619+ @— ) )6 udt
O
t
f Ta
'* bf (1 (§§)* + (% F)* )dt + (% G)* 15a
1‘0 (“tr

where 5i(t) satisfies the variational equation
-.‘. 6f - 8f - 8f
-ox=( ) 6x+( ) 8u+( ) 50.
E a E a E a-
with

5i(to) = 5x(to) .. (it), ,3 etc
’ 0

(3.3.18)

to OtO + (95%).}- OD

(3.3.19)

(303.208.)

(3.3.20b)

 

 

 

 

 

 

 

 

 

I 1 1 J J

 

+A
to to to t t

Figure 3 A nominal solution are and its neighboring
arc with variable end-points

115

6i(tr) = 5x(tf) - (22),,”f Ctr (3.3.20c)

In deriving (3.3.20), the relations (3.3.17) with f replacing
F are used. In this case, the variational equations (3.3.20)

and (3.3.11) are equivalent.
1:

f
Integrating / t A
O
to satisfy the following adjoint equation and the terminal

T15:"Edt by parts and choosing Mt)

condition
{T = - lag-jg) - (331E) (3.3.21a)
{- «n-
T _ as
1 (tr) - (an t (3.3.21b)
' f

and substituting these relations into (3.3.19). we have

.. ‘ T T. ,-
5Q - (G + F)*,tf6tf + (A 5X)*’to - (A x + F)-11-,t°0to
t1.
+ (§§)H61s+/ 11%;”) +(ag-F*)16u<it
t0
t‘f
+ 1/ (3M (1*) + 1%,- F*) hung; G) 16a (3.3.221
tau-,1:
o f
where
(é) = (86' + 8G i) . (66' + KT.
*’tf 3.5 E *,t 3'1? x)-::-,t

1‘
We shall call by the variation Of y, and by the total
variation of y. Therefore 6Q is the total variation of Q

L16

corresponding to the total variations OtO, 6t 6X(to).

f,
Ox(tr), 5a and OB, and the variations 6fi(t) and Oi(t). For
convenience we shall replace by by dy, and 65‘ by 5y, thus

(3.3.22) and (3.3.20) become, respectively,

= ' T T- az ..
dQ (G + F)* t dtf + (a dx)*,to - (A x + F)*’todt° + (5'15)ng

’ f
tf
T af OF
+ [to [A (55-h? + (5-5)*]8udt
t
+ 1/ f (9%) + (2%) )dt + (32 1a: (3.3.23)
to '31- «:1- *,tr
and
oi = (§§)*Ox + @331. + (§)*m (3.3.211a)
with
5x(to) = dx(to) - (31).} t <11;o (3.3.2Ab)
’ 0
anti.) = u(tf) - (a’c),,,tfdtf (3.3.211c)

The solution of the variational equation (3.3.21.1) is
not needed directly. Instead, the solution of the adjoint
equation (3.3.21) is required to evaluate the total variation
Of Q. A different approach to the derivation of the variation
of Q will show clearly that the adjoint variable Mt) defines
the variation of Q in toms of the variations or perturbations

of the control variable u(t), the control parameters a, p,

h?

and the initial and final times and states. For this reason,
x(t) is frequently called the influence function.

It is also important to point out that the adjoint
variable A(t) is a linear combination Of the column elements
of the transition matrix of the equivalent Mayer problem.
Therefore, one may call A(t) the modified transition

function [EV-l]. To show this, let X = (x0 xT)T , where x0
satisfies
= F(t,x,u,a) x°(to) = 0 (3.3-25)

The variational equation correSponding to the new state
variable X is

ox = (% a) OX.+ (——) ou +(% 2:) dc (3.3.26a)

* a

with

OX(to) = dX(to) - X(to)dto (3.3.26b)
where

— F

f = (303027)

f

the corresponding adjoint equation is

AT = — 11%;?) (3.3.28a)

1.1.8

A (3.3.28b)

The transition matrix <I>(t,t') of (3.3.26) is nonsingular and
satisfies the homogeneous linear differential equation (see,

for instance, Coddington and Levinson [CO-1]),
: _ of - _
ai-

Differentiaing the inverse of <I>(t,t'), we have

fisher) = - <I"l(t,t')<f>(t.t')é-l(t.t')
__= _ é-1(t,t,)(§_§)* (3.3.30)
Using the relations
9'1(t,t') = a(t1,t) (3.3.31a)
and
<I>T(t,t) = <I>(t,t) = 1 (3.3.3112)

(3.3.30) can be written as

AT tut = .. 51'; T eT(tv,t) aT(t,t) = I (3.3.32)
( 1 (31*

Therefore <I>T(t',t) is the transition matrix of the adjoint

equation (3.3.28). If we assign tf for t, and t for t',

#9

A(t) can be expressed as

Mt) = aT(t.tf)A(tf) (3.3.33)
01‘

ATM = AT(tf)<b(t.tr) (3.3.311)

Thus we have shown that A(t) (and thus 7((t)) is a linear
combination Of the column elements of <I>.

Now, we shall return to our main subject. In the
expression (3.3.23) for dQ, there is a term which involves
dtf explicitly. This term can be eliminated by using the
stOpping condition S(tf) = 0. Since 8 is defined as a

stopping function, 8 = 0 is always satisfied at t = t and

f,
f.
Q = S(tf), and drOpping the term involving (3 in (3.3-23):

therefore the total variation of S, as, is zero at t Let

we have

dS(tf) = (6° + M°) dt + (1 gdx) + M°)

a,t f f dto

{Wt'(s *’to

tr T of oM°
+ /' [15(55 ~+ (35 ) JOudt
t i:- e:-
O
t

f (360
+1/t (3% :1 +132 C’))<it+3-1<3a--—-o (3.3.35)

Since we have asmed that S )6 0, one can solve dtf from

(3.3.35). For convenience we set

50

(3),“: = (63° + M°)*'tf (3.3.36a)
(<31)Mr = (a + emf (3.3.36b)
Solving for dtf and substituting it into (3.3.23), we have
t t
__ f BHQS f OHQS
dQ — ft (33 )*6udt + [ft (3a- )*dt + YQJdo.
O O
H”) dB + K db (3 3 37)
3F ,, Q ' '

where

H05 = 135(t)f(t,x,u,a) + F(t,x,u,c) - (§)*,th°(t,x,u,a)
(3.3-38)

"os = AQ - 1.5%): t (3.3.39)
’ 1‘
Note that )‘S is the (nxl) influence function correSponding
to S, and (‘8 satisfies the adjoint equation (3.3.21a) and
the terminal condition (3.3.21b) with F replaced by Mo, and
G by 9° . AQ is the influence function corrOSponding to Q,
and AQ satisfies (3.3.21). If we assign J to Q, than F is
replaced by L, and G by g, and AJ is an n-vector, HJS is a
scalar; on the other hand, )‘I is an (nxp)-matrix function
which satisfies (3.3.21) with F replaced by M, and G by 9,

and consequently HIS is a p-vector-valued function . where

51

OH '

Q5 _-. AT 81' 6F .. 9'. 5M0 . e O
(3": )* Qs(gfi)* + (55).“. (8)*,tf(6—u )* (3 3 (4- )
KQ = I: -(HQS)*,t° )‘Qs(to) ] (3.301(1)
de = 1 at, de(to) 1 (3.3.112)

_ as 6) 39°
Y .. - — a a
Q (3;)*’tf (S)*,tf(a )*9tf (3 3 (L3)

where the partial derivatives with reapect to vector variables

are defined in the usual way. For example,

 

 

2!: = 6L 9.1.: , ,
(an * [31:1 0 a e a a e e e Bum] ‘3‘. (3 3 (414')
8M1 1121
36-1 Our
(391* - 3 (3.3.1151
9gp . . . . . . . . 213?
-6011 act’- *
F2131 0 . . . . o ‘1
o . .
az ' ' °
_ a e o . . 6)
(53).”. : 'Zép' (3 3 LL
. 0
6 O O O O O 6
L J *

 

 

52

3.1.1. Constrained Effort and Successive Optimization,
Algorithm I
In the previous section we derived the total varia-
tion of the generic quantity Q in terms of the influence
functions AQ and AS, and the variations of controls and
initial time and state. If we successively assign J and I

to Q, we have, reapectively

ea /tf (aHJS1 ff (BHJS)
= __ 5 dt + dt + Y do + K db
to du * u [ t E a." J] J
° (3.11.1)
ft (a1EIIS 1 /:f aHIS
d1 = 6 dt + dt + Y do
t 3— u [ (a )4? I]
O to
+(-B-) dB + KIdb (3.11.2)

Note that dJ and (11 are Obtained by adjoining the differential
constraint (3.1.2) to J and I and by choosing the adjoint
variables AJ, and 711 to satisfy the adjoint equations (3.3.21).
Thus far J and I are not related except that they are both
constrained by the state equation (3.1.2).

I In the optimization process we shall derive in this
section, dJ is minimized under the constraints (11 by choosing
511(t), (h, d5 and db Optimally in each iteration. Therefore,
we are not trying to minimize J with the constraints I directly.

In our previous derivations, perturbations and local

linearizations were used. The closeness between the predicted

53

improvement in J, dJ, and the actual improvement in J,

( Jnew

I, (11, and the actual correction (Inew - Iold) are heavily

- J 01d), and that between the predicted correction in

dependent upon the degree of validity Of these linearizations.
In order to avoid excessive error, constraint must be placed
on the perturbations 0u(t), d1, dB, and db. To this end, we
shall follow an approach which was first used by Bryson and
Denham [BR-1].
Let (d0) be a'small' positive number, and let U(t),
A, B, and W be mxm, rxr, (1+n)x(l+n), and p'xp' symmetric
positive definite weighting matrices, reSpectively. Set
1:
(do)2 = /’ f ouT(t)U(t)ou(t)dt + daTAda + dewdb + deBdb
t° (3.11.3)
thus, (dC)2 is the sum of the cumulated weighted squares of
the control variations Ou(t), d1, dB and db. Since (6.0) is
a measure of the effort in the (u,o. ,B,b)-Space, we shall

call (dC) the weighted control effort or the iterative step-

 

size in the (u,c,0,b)-Space-—the control effort Space.
3 In order to take account of the constraints (11 and
(dC)2, adjoin (3.11.2) and (3.1.1.3) to (3.11.1), to Obtain

t

f 5H cm
(H = f ( JS) - hT( IS - cOuTU oudt
to [ 5T." 3:- 5-1: )s'. 1
tr aHJS T aHIS T 213
+[ft ((3; ) -h(53 ))dt+YJ-hYI-cdo. Mda
'25 t!-

O

5L1-

g) + chTWJda + [KJ .. hTKI - cdeB]db
'11-

+ thI + c(dC)2 (3.11.11)

where h is a p-vector multiplier, and c is a scalar multiplier.
The variation of dJ corre3ponding to the perturbations
on Ou, d1, dB, and db is

t

f as as
6(dJ) = f ( JS - hT —13 - 2 o TU 62 dt
t13—51* (au1* quu
0
[tr ((aHJS T 8HIs T T 2
+ -h )dt+Y-hY -2cdaAd
I: to a )* ('55: )* J I J 0'
- [hT(g-§-) + 2chTW]d2(3 + [KJ - hTKI - 2cdeBJd2b
‘K'
T 2
+ 11 5(d1) + 05(dC) (3.8.5)

The stationarity of dJ due to perturbations in 011, d1, dB,
and db implies the vanishing of the variation of dJ, 6(dJ).
By applying a fundamental lemma in the calculus of variations,
see, for examle, [GE-l], the vanishing Of 5(dJ) implies the
vanishing of the coefficient of 6211 for all t E [to,tf], and
the vanishing Of the coefficients of (120., d2(3 and dab. Note
that 6.1 and (d0)2 are Specified quantities. Therefore, the
variations of <11 and mm2 are zero. By equating the

coefficients of various quantities to zero, we have, after

taking transposes

U-l

8H

8H

_ JS T I

511 "' E; [(a-fi' )* - (5";
tr
A-1 aH
da = 331/ ((EJSF-
1:0

d5 = - §-’1(5-9* h

_ 3'1 T T
db — 23 (KJ - Kih)

55

S1Th1
'31-

SH
(80'.

'2!-

In order to shmplify'the expressions,

notations:
= + '-
QII U11 AII + BII + wII
_ = T
PIJ ’ UIJ'+ AIJ + BIJ PJI
PJJ = UJJ + AJJ + BJJ
tf
aH 18H
_ IS IS T
UII -./f (5— )* U %( u )*dt -
0
t
f on
t a a
O
f an an
_ fiJs -1 JS T _
O
= - T = T
AII RlsA 1RIs AII
- -l T _ T
AIJ “ RISA RJs "AJI

T
‘ QII

I
C‘.

—IS)Th)dt-+'y§ - Yin]

(Pxp)

(pxl)

(scalar)

(3-h.6)

(3ch.7)

(Bohwe)

(3.h.9)

we define the following

(3.A.1o)

(Bohcll)

(3.A.12)

(3.h.13)

(3014'01-(4')

(Bah-15)

(3.h.16)

(Bah-017)

56

.A.JJ = RJSA'J'RES = AEJ (3.1+.18)
an = 193%? = 3% (34.19)
BU = KIB'J‘K}: = BE]: (3.LL.20)
BJJ = KJB'J-Kg = 33.1 (34.21)
_ oz -1 az T T
w — w = w . .22
II (35')* (330* 11 (3 LI- )
t
f 6HIs
RIS = ft (a )*dt + YI (PX?) (30’4023)
O
t
r GHJS
RJS = ft (go-L )*dt + YJ (lxr) (3.1+.2LL)

0

Substituting (34.6) to (3.h.9) and their transposes into
(3.L|..3), and using the notations (3.14.10) to (3.1;.21i), we have

T

(do)2 = fit?” - Pain - h PIJ + hTQIIh] (3.1+.25)

Similarly, substituting (34.6) to (34.9) into
(Boll-.2), and 113138 (Bah-010) to (Bel-Leah), “’3 have

Solving for h from (3.14.26), we have

b = q'I'JflPIJ - 2ch] (3.1M?)

57

Substituting (34.27) and its transpose into (34.25), we have

(<10)2 = i? [PJJ - Pi). Qfi PU] + dIT 42% 6.1 (34.28)

Solve for c using (34.28) and assuming that

2 T -1

we obtain

T -1 t
c = t 1 PJJ " PIJ Q11 PLI ’ (34.29)
2 (do): dIT Q52} (11

 

Then substituting c into (34.27), we obtain

T -1 it
P - P P
-1 JJ IJ Q11 IJ
h = Q P - (3:) T“ 6.1 (34.30)

 

Substituting (34.29) and (34.30) into (34.6) to (34.9).
we obtain the Optimal perturbations°

 

5H
._ -1 IS T -1
511(t) - U (13) [ (33 )* Q11 (11
an an (dC)2-dITQ dI $2
1' (___JS)T _ (___IS)T Q-J. P 111
an * au * II IJ PJJ-PEIJQII PIJ
(3.1+.31)

 

1
2 T -1 ‘2
(dC) — dI 2Q (11
_-lT-1+TT-1 II
d“ ' A [R15 Q11 ‘11 " (RJs'RIs Q11 PIJ)( T .-1

PJJ " PIJ QII PIJ
(Boll-62)

58

 

 

2 T -1 1a
(do) - dI Q (11
d" = W1‘§%)T QEI ‘11 " (“PIA T $1 ) ]
* PJJ ' PIJ Q11 PIJ
(3.u~33)
2 T -1 %
(do) - (11 Q (11
_ -1 T -1 T T -1 II
db " B KI QIIdIi(KJ'KI QII PIJ)( P Pf P ) ]
JJ ' IJ QII IJ
(Bah-3h)

Finally, substituting (34.31) to (3.4.311) into (34.1) we
obtain the predicted optimal dJ

35 32
dJ __ PIJ Q11 d1t((dc) .. dI Q11 ‘11) (PJJ "' PIJ QII PIJ)
(34.35)

Since the quantity
t *2
2 T -1 T _. -1
((d0) " ‘11 Q11 dI) (PJJ " PIJ Q11 PIJ)

is nonnegative, we have the following rule for choosing the
signs in (3.14.31) 1:0 (Lil-.35):
If J is to be minimized, choose the "-" sign
associated with the radical terms in
(34.31) to (34.35); Ir J is to be maxi-
mized, use the "+" sign instead.
Since we are discussing the minimization problem, we shall use
the "-" sign in our subsequent develOpment.

If we set d1 = O in (3.14.35), we have

1»:

_ T -1
dJ - - (dc)(PIJ - PIJ QII PIJ)

59

or
6.1! 1

66 = " (PIJ " PEEJ Q13 PIJré ‘3'“36)
Recalling that J is the cost functional which we want to
minimize and that (dC) is the weighted control effort, we
therefore, may interpret the right-hand side of (34.36) as
the negative gradient of J in the function space of control
effort.

By examining the Optimal perturbations (3.11,.31) to
(34.31;), one may recognize that each of the optimal pertur-
bations consists of two parts; one of them is for achieving
the error correction <11, and the other part is for the
minimization of M. The weighting between the two parts
depends on how large an isoperimetric error correction (11
was specified.

The predicted dJ in (3.14.35) also consists of two
parts; the term PEJ QE§ 6.1 is due to the effort made to
accomplish the correction of I, <11, and the rest of (3.1L.35)
is due to the effort made toward minimization of the varia-

tion of the cost functional dJ.

3.5 An Iterative Procedure for Algorithm I
The following iterative procedure is suggested:
(1) Initialization: Guess and store the nominal u*(t), (1*,
(3* and b*.
(ii) Forward integration: Integrate the state equation
(3.1.2) from to until 8 = 0, and store x*(t). Also
integrate L, M, and M°, add to g, e + z, and 6° at

(iii)

(iV)
(V)

(v1)

(vii)

(viii)
(ix)

(x)

(xi)

3.6

60

tr, respectively, to obtain J, I, and S.

Check the validity of linear approximation. Modify
the value of (d0). Stop the iteration if (10 S (101. .
Evaluate (3/é)-n-,tf , (i/é)*,tf , YJ , and Y1 .
Update the symmetric positive definite weighting
matrices U'1(t), A'l, B'l, W'l.

Backward integration: Integrate the adjoint equation

(3.3.21) from t to to to obtain lJ(t), AI(t), and

f
ls(t). Simultaneously integrate to obtain UII’ UIJ ,

UJJ , AII , AIJ , and AJJ , and store (BHJS/au)*’t

and (BHIS/Bu)*,t .

Evaluate BII , BIJ’ BJJ , and WII to obtain QII , PIJ
and PJJ and invert QII .

Cheese suitable dI and (dC) with (dC)2-dIT-QII d1 2 0.
Compute the optimal 5u(t), da, dp, db and the
predicted dJ.

Simultaneously with (ix) update u(t), a, (s, and b,

and store both new and old values of u(t), a, (i, and b.

Go back to (ii).

Some Measures for the Validity of Linear Approximation

As mentioned before, this optimization process

depends heavily upon linearization. Therefore, it is impor-

tant to check the validity of this approximation. In this

section we shall define some measures of the validity of

linear approximation in terms of the predictabilities of dJ

61

and d1.
Let eJ and eI be the relative prediction errors of

dJ and d1, respectively, and define them as

 

 

dJ - (J - J )
OJ = 23w Old (306.1)
and
HV(dI - (I - I ))H
_ new old
"' e 02
°1 HVdIn (3 6 )

where V is a diagonal weighting matrix, and the norm H'H is
defined as

(vu = sup Ivil (3.6.3)
i=1,..,p
_ 1 p T
where v - (v ... v ) is any bounded p-vector.

Let eJO and eIo be specified tolerances for eJ and
eI repectively, then eJ and eI can be used to define test
criteria to determine the step-size for the iterative process.

In practical optbmization problems, Judicious test
criteria must be used, otherwise great computational diffi-

culty may result.

3.7 The Instantaneous Algebraic Constraints and the
Penalty Function Technique
So far we have treated the control problem.of Bolza

with differential constraints and isoperimetric constraints

62

only, but we have not treated the problem with instantaneous
algebraic constraints given in (2.11.).3). For convenience,
we shall redefine the instantaneous algebraic constraints as,

¢8(t,x,u,a) ,

//
O

3 = 1,000’Q' (307018)

¢3(t.x.u.a) - q'+1.....q (3.7.1b)

I
O
to

II

more a is defined on E and is of Cl(E). And as can either
be a control variable constraint, or a state variable cons-
traint. Similar problems have been treated analytically by
many authors, see, for instance, Gamkrelidze [GA-1],

Berkovitz [BE-1,2], Hestenes [HE-3] and Guinn [GU-l]. Problems
with state variable constraints are more difficult than those
with control variable constraints in both analytical treatment
of the problem, and from the computational point of view.
Dreyfus [DR-l] considered the control problem of Mayer with
terminal equality constraints and Kth order state variable
inequality constraints, where both the inequality constraint
and the control variable are single scalar quantities. Denham
and Bryson [DE-l] obtained a steepest-ascent solution to this
problem by employing this direct approach to inequality cons-
traints. Denham [DE-2] discussed the case of two inequality
constraints with one or two control variables involved.
Although this method seems feasible for single constraints

and single control variable problems, it becomes more and

more involved for systems with increasing numbers of constraints

63

and control variables. A method known as the penalty function
technique was first introduced by Courant and Moser [COU-l]
for ordinary minimum problems. Kelley [KB-2] extended the
penalty function method to the control problem of Mayor by
introducing additional state variables. Okamura [OK-l]
applied the penalty function method to the control problem of
Lagrange with side constraints and established some mathe-
matical basis for this method.

There are a number of reasons why one would wish to
use the penalty function technique: (1) Due to the form of
the cost functional of our problem, no alternation is needed
in our algorithm, whether or not there are instantaneous
algebraic constraints. (2) No additional state variables
need to be added to our original system. (3) This method
is valid for any finite number of constraints and any finite
number of control variables.

Let PW) be a vector penalty function for the cons-
traints (3.7.1), the components of P(¢) can be defined as

PS(¢8) > o if as > o for s = 1,...,q' (3.7.2a)
P‘(¢‘) = o if as s o for s = 1,...,q' (3.7.2b)
P’(¢“) > o if as); o for s = q}+l,...,q (3.7.2c)
P°(¢’) = o if fl“ = o for s = q}+l,...,q (3.7.2d)

whore Paws) (s = 1,...,q) are real-valued continuous functions

61;

defined for all finite values of 33.

Let w be a q-vector, whose components wl’,...,wq
are real nonnegative numbers. We shall call w the weighting
vector for the vector penalty function PM). The weighted
penalty function is then defined as wTP(¢). Let JP be the
penalized cost functional defined as:

t;

. f T
JP = J + / w P(¢)dt (3.7.3)
to
or
t
f T
J,=s+/t[L+menfi 6J4)
0

If we define

L = L + pr(¢) (3.7.5)

P

then the penalized cost functional J can be written as

P
tr
J = g + f L dt (3.7.6)
P t P
0
Note that JP and J are of the same form. The quantity
tf T
/ WP(¢)dt
1so

is a measure of cumulated violation of the instantaneous
algebraic constraints.

By suitably choosing the weighting sequence {w} .
the corresponding sequence {JP} of the approximated control

65

problem will converge to J of the original problem under
some convergence hypotheses on the initial state, terminal
time, and control functions associated with the penalized
system [OK-l].

Therefore by introducing penalty terms to J, one
can apply Algorithm I to the general control problem of Bolza.
Further , it should be emphasized that the method applies
equally well to problems with control variable constraints,

or state variable constraints, or with mixed constraints.

3 . 8 Conclusion

In this chapter we have first formulated a compu-
tational version of the control problem of Bolza. Variations
of the various quantities of the problem are derived by using
the notion of variable region of integration. An algorithm
for the computation of the problem is then developed. A
penalty function technique is used to treat the instantaneous
algebraic constraints, and additional control parameters are
introduced to convert the inequality isoperimetric constraints
into equality constraints.

Before closing this chapter, we shall illustrate
some advantages of this formulation of the optimal control
problem over that of the earlier formulations, such as those
given by Bryson and Denham [BR-l], Kelley [KE-l], and
Vachino [VAC-l]. The present formulation is based on the

problem of Bolza, where the earlier ones are based on the

66

problem of Mayer. From the theoretical point of view, as

demonstrated in Chapter II, the problem of Bolza and that

of Mayer are equivalent. However, from the computational

point of view, this is not the case. For instance, if one

has a control problem of Mayer, the present computational

algorithm is readily applicable. One the other hand, if

one has a control problem of Bolza with isoperimetric

constraints and if the stopping function S also has an

integral part, then, in order to use the earlier formulations,

one will have to transform the problem by introducing (p+ 2)

additional state variables, and (16+ 2)‘2 additional adjoint

variables. Thus the number of equations needed to be integrated

in each iteration increases from nx(p+ 3)+p+ 2 to (n+p+2)x(p+3),

or a net increase of (p+ 2)2. This will increase the computing

time as well as the memory required to store the state history

and other quantities. The equations in the present

formulation of the problem are not significantly more complex.
In the next chapter we shall identify an important

subclass of problems treated in the present chapter; they

are so important that a separate chapter is devoted to

their study.

IV A SUBCLASS OF THE CONTROL PROBLEM
OF BOLZA, ALGORITHM II

The control problem we shall consider in this chapter
is a very important subclass of the control problem of Bolza.
Since this subclass contains a large number of problems of
practical interest, we shall treat this problem separately
from the general problem considered in Chapter III, and an
algorithm will be derived for the numerical solution of the
problem.

Bang-bang control problems and problems with both
bang-bang control variables as well as continuous control
variables with various constraints are contained in this
category.

Denham and Bryson [DE-1] considered a bang-bang
control problem with a single control variable, and an a
priori specified number of switches of the control variable.
Vachino [VAC-2] considered the case of multiple control
variables, and without the need of a priori information on
the number of switches. We shall consider more general
problems and we shall not assume an a priori knowledge of

the number of switches.

(4.1 Formulation of the Problem
me cost flmctional is defined as

67

68

t
f
J = g(tf,x(tf))+ ft L(t,x,u,v)dt 04.1.1)
0

The state equation (differential constraint) is

2.: = f(t,x,u,v) (Ll-0102)

l uml)T

where x = (x1 ... xn)T is the state vector, u = (u ...

m
is the continuous control vector, and v = (v1 ... v 2)T is

the discrete control vector. Assume u(t) is piecewise
continuous on [t°,tf], and v(t) 6 ( 1:1,sz , or

Vj(t) 6 {1%. kg} , j= 1,...,m2, wherlrle k
.. l 2 T
i "' (k1,...’ki)
is no loss of generality in assuming that k{ < kg for all

.1 J
1 and k2 are any

fixed real numbers, and k , i = 1,2. There

:1 = 1,...,m2 . The instantaneous algebraic constraints are

similarly defined

as S. O 8 = 1,...,Q' (LI-01033.)
a“ = o s = q'+1,...,q (h.i.3b)
where (is = ¢8(t,x,u,v). Define the isoperimetric constraints
II S o y = 1,...,p' (b.1421)
IY :.- 0 Y = phi-1,...,p (14-01-0143))
where
t
- 1'
IT = 97(tf,x(tf)) + ft MY(t,x,u,v)dt (14.1.!4c)

O

69

The stopping condition is determined from

S = 0 . (LL-1.5a)
which defines 'the terminal time tf’ where
tr
_ S(tf) = 9°(tf,x(tf)) +./; M°(t,x,u,v)dt (h.l.5b)
o

and it is assumed that S 7‘- O-for all t E [to,t The assump-

r3'
tions of continuity properties for the various functions are
essentially the same as those defined in Chapter III for the
general problem, and notions of admissible quantities are
similar to that defined in Chapter III.

As we did in Chapter III, we shall treat the instan-
taneous algebraic constraints by using the penalty function
technique, and we shall transform the inequality isoperimetric
constraints to their equivalent equality constraints. Let
B = ((31 ... Bp')T be an additional control parameter, and

let Z be a p-vector with components

- (pY)2 Y = 1,0e09p' (LI-01.68.)

N
..<
I

P'+1,o.o,p (hele6b)

zY = 0 Y
Set I = (11 ... 1P)T and let

IY=IY+ZY Y

1,...,P

rIhen

70

IV = 0 Y = 1.....p

are equivalent to the constraints defined in (14.14.).
We want to determine the controls u(t), v(t), (5,
and the corresponding state x(t) such that the cost functional

is minimized and the constraints are satisfied.

Li.2 A Multiple-stage Formulation of the Problem

The problem formulated in Section I4..l is a single-
stage problem. Due to the discrete nature of the control
variable v, it will be convenient to treat the problem as a
multiple-stage problem. Since the number of stages is deter-
mined by the number of discontinuities of v, and in general
the number of the discontinuities is not known for a given
problem, we, therefore, have to treat the problem as a
variable-multiple-stage problem. The following formulation
is similar to that used in [VAC-1].

Let u(t), t 6 [t°,tf] be a nominal continuous control
function, and v(t), t G [to,tf] be a nominal discrete control
variable. Since v (nominal) is discrete, both of its values
and switching times have to be specified. Let t3, ,
r3 = 1,...,N‘1 , be the nominal switching times forjv'j. For
convenience, assume to < ti and ti?I3 < tf for all
3 = 1,...,n.2 . Therefore to and t1. are not considered as

switching times. Let t(N) -- f t1,...,tN( be the ordered set

< < < = J
of switching times of v with t1 t2 ... tN and tr trj

71

for at least one 36 {1,...,m2( , and for some rj€{1"”’Nj} .
Note that N S N1 + ... + Nm . With these notations, we shall

2
reformulate the problem in Section L)..l as the following

multiple-stage problem: The cost functional

t
N+l r
' J = g(tIIH-l’xlhl‘th-l)) + r21 ft Lr(t’xr’ur’vr)dt
r-l

(h.2.1)

where tN+l is the terminal time t and xr, ur, vr are the

1"
state vector and the control vectors for the rth stage

respectively. The function L1‘ is defined as

Lr(t,x ,u ,v ) = L(t,x,u,v) , t E T , r = 1,...,N+1

r r r r
(I4..2.2)
where Tr = [tr-l’tr]'
The state equation for the rth stage is:

xr = fr(t,xr,ur,vr) , t e Tr , r = 1,...,N+l (4.2.3a)

23e009N+1 (LI-0203b)

x(t ) r

r r-l) = xr-l(t

r-l
where

’ fr(t,xr,ur,vr) = f(t,x,u,v) , t 6 Tr , r = 1,...,N+l
(LL-2J4)

The isoperimetric constraints are

IY = 0 Y = lgeee'p (LI-'205a)

72

t

Y Y Y N+1 r Y
I - e (tN+1,xN+l(tN+1)) + z «3)4- z. /" Mr(t,xr,ur,vr)dt
r=l tr-l
2 (h.2.5b)
ZY = (fiY) a Y = lav-onp' 3 ZY = O 1 Y = p'+l""’p
(u-Z-SC)
and the stopping condition is
S = O (h.2.6a)
o N+1 .tr 0
S“"N+l) = e (tN+l’x‘N+l(tN+l)) + Z ‘/ Mr(t’xr’ur’vr)dt
r=l t
r-l
(h.2.6b)
where
M;(t,xr,ur,vr) = MY(t,x,u,v) , t e Tr , r = 1,...,N+l
Y = 09-00:?
(II-2.7)

also assumes it o for all t e [to,t Note that v13. is

f].
constant in Tr for each ,1 and r, and furthermore, v3 may
be constant in more than one interval.
The problem is to determine the control variables

u(t), v(t), the control parameter (3, and the corresponding
state x(t) such that the cost functional J is minimized and
the constraints are satisfied.

I We shall divide this problem into two parts. One is
to consider the problem with a fixed number of stages, and
to obtain a computational algorithm for determining u(t), (5,

and the switching sequence of v(t). The other part is to

73

consider the problem of varying the number of stages,
especially, to add new stages if it is desired.

Ii.3 The Variations of the Multiple-stage Problem with
Fixed Number of Stages
The variations of quantities of interest can be
obtained in a similar way as those obtained in Chapter III.
We, therefore, will omit some of the details involved in
obtaining the variations for the present problem.
Let Q be the quantity of interest,

I:

' N+1 r
I) + 2(9) + 1:51 ft Fr(t’xr’ur’vr)dt
r-l

Q'= G(tN+1'xN+1(tN+l

(h.3.1)
where Q can be J, I, or S. The term 2((3) is dropped if we
designate J or S for Q.

Let x: be a nominal solution are of the system

(ll-0203). Then

r = l,...,N+l
(h.3.2)

is the explicit expression for the are x: . If we perturb the

nominal I; by perturbing u: , v: , (3*, to , "”tN+1’ and

x:(to), xN*+1(tN+l)' and adjoin the differential constraints

x: : t.x§(t).u§(t).v§<t).a* . t e T.- .

(LI.l.3) to Q, the corresponding total variation of Q is

7h

N+1 / r T arr aFr
dQ.= Z [A (5——) + (-—-) Jéu dt
r=l tr-l r “r * aur * r

+ [ifdxl - (@221 + slump,to + [G + FINJp’tzfduslf + (§.— f) as

N
+ 2 [(IT rxr+ Fr ) - (AT

r=l r+1 5W1 Fr+1H..,trdtr (u.3.3)

'vdlere we have used the fact that the state xr(t) is continuous,
and therefore, the total variatlons dxr(tr) and dxfi1(tr)
are equal, correct to first order for all r = 1,...,N. Also

oxr satisfies the variational equation

8f 61‘
= r 1‘ ‘ ::
(fin-511, + (5;)*Our , t E Tr , r 1, . . .,N+l

xr
(LI-030143)
and
5xr(tr-l) = dxr(tr_ 1)- xrt( M1)dt -l , r = l,...,N+l
(u.3.ub)
and Ar satisfies the adjoint equation and the boundary
conditions
of BF
'T _ T r r . =
Kr " - Kr(§i;)* - (ESE—1.).” , t 6 Pr , r 1,000’N+1
(h.3.5a)
T __ T _
Kr(tr) xfi1(tr ) 1' - l,eee,N (“.305b)
T _ aG (h.3.SCI
l (t ) — (3———-)
N+1tN+1 xN+l *,t

1'

Note that the only changes of vr considered here are the

75

r 0
Now we shall use the stOpping condition to eliminate

changes of switching times Tr-l and tr associated with v

t from dQ, Substitute S for Q, in (I4..3.3), and set A31, for

f
)(r in (h.3.3) to indicate that Ar is associated with s.
Furthemore, by the same reason employed in Chapter III,

setting dS = 0 , and solving for dtf , we have

t o
N'I'l 1' 1 of 5M
_. T r r
dt--Z/ —;— [6 (3—I+( )]5 at
f 1‘31 t (S) Sr ur * 311-; * LIT
r-l *,tr
1 I (IT )dt
' T [ Sldfi ' 81x1 Ml 1*,
( Tit- t
’ f
2T --1— [(IT + 11°) (AT M° )] dt
F1 (5‘).“- t er r Sr+l xr+1+ MrI-l *,tr
' f
(h.3.6)
where A satisfies the adjoint equation and boundary condition

Sr
04.3.5) with G = 9° and Fr = M: . and where

(é)*stf = [5° + MN+1]*,tr (LL.B.7)

Substituting dtf into dQ, we have

N+l tr HQSr o I:
dQ - :1 ft 1 (7:4,, ur dt + r—-1[HQ’ST - HQSfi1]*,trdt
r...

where

76

= T _ é o ...

HQ” T‘QSrfr + Fr (5)*,thr , r 1,...,N+l (u.3.9)
A = A - A (ET-)1. r = l N+l (LI- 3 10)
er Q1. Sr é*’tf 3000, e e

_ T
de = [dt dx$(t )] (u.3.13)
O O

and IQ}, satisfies the adjoint equation and the boundary
condition (14-0305).

I44 Successive Optimization Process for the Problem with
Fixed Number of Stages
Let Ur(t) (r = l,...,N+l), vr (r = 1,...,N). B and
W be (mlxml), (1x1), (1+n)x(l+n), and (p'xp') symmetric
positive definite weighting matrices, respectively, selected
in each iteration to provide suitable weighting and faster
convergence. Let (do) be the step-size in the control effort

space. For the same reason as that in Chapter III, set

1:

2 Nil 1' T .. 15
(dc) - 121/1; our(t)Ur(t)our(t)dt + rildtrvrdtr
r-l
T T
+ db Bdb + d(3 we (II4.1)

77

If we substitute J and IforQ in (t.3.3). we shall
obtain the eXpressions for dJ and d1 reSpectively. .Adjoin
the dI and (dc)2 constraints to dJ, the stationarity of the
constrained dJ corresponding to variations in Our , dtr , db
and db yields the following Optimal variations and the
predicted optimal value of dJ:

OH
_ -1 ISr T -1
6111,”) - Ur (t) (T—ur )* QII d1

 

2 T -1 35
_ ((aHJSr)T _ (aHISr)T Q-1 P (do) ' ‘11 QII ‘11
3u an II IJ T' -l
1' * 1‘ PJJ ' PIJ QII PIJ

t e T , r = 1,...,N+l

 

r
(halt-2)
_ -l T T -l
dtr vr (HISr ' HISrTl)* t QII dI ' ((HJSr I HJSr+l)* t
’ r ’ r
(dC)2- dIT q'l dI 3"
- (HT - HT ) Q‘1 P ) II
ISr ISr+1 * t II IJ P PT Q—l P
’ r JJ - IJ II IJ

r = 1,...,N (h h-B)

2 T -1
(d0) - dI QII dI

T .1
PJJ " PIJ QII PIJ

 

4.
d5 = WHSE): afi a: + PIJ( ) (tint)

db = 3'1

 

T .-
T -1 T T -1 “0’2" ‘11 QII ‘11 T
KI QII ‘11 ' (KJ " KI QII PIJ) P _ PT Q—l P
JJ IJ II IJ

(Lula-.5)

78

= T -l 2 T -1 a T -l a
dJ PIJ'QII d1 ' [(dc) ' d1 QII d1] [PJJ ' PIJ QII PIJ]

04.4.6)
Note that, except for dtr , the expressions for the Optimal
Our(t), dB, db and dJ are essentially the same as the
corresponding ones obtained in Chapter III. Except in the
expression for dtr , we can drOp the subscript r and let t
run from toto tf instead of from tr-l to tr .

We have used the following notations:

t

N+1 r 5H 8H
_ ISr -l ISr T

U - z /‘ ( ) U ( ) dt (h.u.7)

II r=1 tr-l Eur * r Eur *
v = T (H - H ) v‘1(HT - HT )

II r=1 ISr ISr*l T’tr r ISr ISr+l T’Tr

(huh-08) '
andUIJ , UJJ , VIJ , VJJ are defined similarly. And W ,

II
BII , BIJ , BJJ , QII , PIJ , and PJJ are precisely the same

as we defined in Chapter III.

Ii.5 A Computational Procedure for the Problem with Fixed
Number of Stages
For convenience, we shall name the following procedure
the Computational Procedure A:
(i) Initialization: Guess and store the nominal
t(N) = {t1,...,tN} , u:(t), v:(t), t e Tr , r = 1,...,N+l,

(3*, and b*.

(ii)

(iii)
IIV)
(V)

(VI)

(vii)

(viii)
(ix)

(X)

79

Forward integration: Integrate the state equation
(u.2.3) from to until 3 = o (i.e., until t = tf),
and store xT(t). Also integrate L, M, and.M°, and
add respectively to g, e + z, and 9° at tf to obtain
J, I, and 8.

Except for the first iteration, check to see if the
prediction errors are tolerable. Modify the value

of (do). StOp the iterative process if dc S d.Cf .

Evaluate (J/S)*,tf and (I/S)*,tf at tf .
Update the weighting matrices ‘U;l(t), V'l, B'l,

r
and‘W'l.
Backward integration: Backward integrate the adjoint
equation (h.3.5) fromtf to to obtainkJr , TIr ,
and TSr , and simultaneously evaluate the integrals
and sums to obtain'UII , UIJ , UJJ , VII , VIJ and
VJJ and store (aHJSr/aur)*,t , (aHISr/aur)*,t ,

(H - H and (H — H

JSr JSr+l)*,tr Isr ISr+lT*,tr '
Evaluate WII , BII , BIJ., and BJJ at to to obtain
Q11 , PIJ and PJJ . Invert QII .

Choose suitable dI and (dC).'

Compute the Optimal Our(t) (r = 1,...,N+l),

dtr (r = 1,...,N), ds, db and the predicted dJ.
Simultaneously with (ix) update ur(t), tr . B. b,

and

80

T

new r = ”new r-l ’ tnew r] = [t + dt

r-l r-l ’ tr+ dtr]

Store both new and old values.
(xi) Go back to (ii).
Note that, although this procedure will not add any
new stages to the system, it has the capability to decrease
the number of stages by moving consecutive switching times

to the same position.

I4.6 Optimization Process for Adding a New Stage and the
Corresponding Computational Procedure
So far we have only considered the problem with fixed
number of stages. In this section, we shall consider the
problem of adding new stages. There are two important

questions we have to answer in this section: first,
"Should we add a new stage?", and secondly "If the answer
is yes, where shall we add?".

For simplicity, we shall consider the case of adding
one stage at a time, and we shall not consider the isoperi-
metric constraints (I4.2.S). The generalization of this
simple case to a more general treatment is not difficult
and shall not be considered here.

The unperturbed system 04.2.3) with N + 1 stages is
called the nominal system. x“. is used to denote the nominal
solution are of (14.2.3), and ur(t) will be held unperturbed
for all r = l,...,N +1.

at t =

[t',t'4-dt'](: (t

The perturbed system.can be expressed

t!

81

Let a subarc (or stage) be added to the nominal are

x; = f; t 6 T1‘
xr(tr-l) — xr-l(tr-l)
where f; is defined as
f; = fr(t,xr,ur,vr) ,
r; = f;,
f; = fr(t,x£,ur,vr) ,
and f;, is
f;' = fr,(t,xr',ur',v
' f;, = fr,(t,xr,,ur,,v
f;, = fr'(t,x§',ur,,v

.As a first order effect, the partial variations Of fr

Let Q be the quantity of interest defined in

with duration dt' ,

r!

I
r!

r!

rg_l 9 t1”) 9

t e T ,

I‘

t E T

r!

t e Tr ,

)
)

)

where

r' E {1,...,N+1I
as
1,...,N+l

2,000,N+1

lgeee’r'“l

r'+l,...,N+l

[t t')

r'-l.’

[t',t'+dt']

t 6 (t'+dt',tr,]

(h.6.1)

(h.6.2a)

(h.6.2b)

(n.6-38)
(h.6.3b)

(h.6.3°)

(u.e.ua)
(u.6.ub)
(u.6.uc)

(h.3.1).

and Fr

‘with respect to the variation of vr due to the addition of a

82

subarc can be considered as impulses occurring at t' with

strengths
If}: -fr'TIt,dt' = [fr,(t,xr,,ur,,v1'fl)-fr,(t,:5‘,,ur',vr')]lt'dtt
(h-6.5)
and
[Fin Fr'JIt,dt' = [Fr'(t’ xr',ur',vr')- F’r'“: x I'"ur"vr')]lt,dt‘

(14.6.6)
where F1)" is defined by replacing f by F in (14.64).
After adjoining the system (1.4.2.3) to Q, the total

variation of Q is

N+1 [tr (oar) T arr ,
dQ=Z [ Ox+?\(( Ox-Ox)dt
r=l tr-l SEQ.” r r Fiji, r r]
+ [(3%)dt + (XNfil)d-XN+1]*’ tf - [Fldt]*’to
N
+ 2: [Fr -F 1 dt + [F dt],
1‘]. xvi-1 *,tr I' N+1 4-,tf
T
+ [Arl(fz"1 -frt) + (Fél-Fr|)]lttdt' (LI-.607)
t

r 0
Integrating / A: Oxrdt by parts and choosing I to
t

satisfy the r-l adjoint equation and the terminal
condition (14.3.5),

83

. = ' T To .
dQ [G + FN+1]*’tfdtf + [Aldxl _ (ilxl + Fl)dt]*,to

N
+ 2 (xT rr + Fr ) - (xT ) dt

F1 [ r+1f+r+1FrI-l 1s,tr r
+ [Ag (fb -rr,)-+ (F! Fr,)]| dt' (u.6.8)

1;!

If we substitute 8 and J for Q in (h.6.8), set dB to zero,

solve for dtf from d8 = O and then substitute dtf into

dJ, we have

N
dJ=de+Z[H -H ] dt
J r=1 JSr JSr+1 *,tr r
dt' (u.6.9)

+ [H33r' 'HJSr'1lt,

where HJSr , KJ and db are defined.in (u.3.9). (u.3.12),
and (u.3.13), respectively, with Q replaced by J. Hgsr.
is defined as

) M°' (“3.10)

1"
*,tr

I _. T 5'
HJSr' kJSrvfrt + Lrv “ (5

where AJSr , (J)*,tf , (3)-313%. are defined in 04.3.10),
(14.3.11), and (h..3.7) respectively.

It is understood that fr, , Lr, , and.M:' at t'
are evaluated on 316*, furthermore, xr'(t')' ur'(t') in I}, ,

Lé, , and Mg: are also on 2*, but v£,(t') is different from

vr,(t').

8h

In the iterative procedure we shall give later, the
process will be iterated a number of times to minimize dJ
with optimal choices in 6ur , dtr , dB, and db, and the
number of stages will be held fixed (with possible decreases)
during the iterations. Then the process will change its
phase to determine whether and how to add a new stage.

In order to make the process simple and workable,
we shall hold “r , tr , and b fixed in the process of adding
a new stage. Setting dtr and db to zero, 04.6.9) becomes

dJ = [H'

JSr' 'HJSrv3|t,dt'

= AHJSr,(tt)dtt (h.6.ll)

Choosing dt' positive, the minimization of dJ is equivalent
to the minimization of AHJSr,(t') . Note that if v is not

changed or if there is no stage to be added, then AH = O.

JSr'
Therefore any new value for v which will minimize dJ must make

AH JSr , negative .
Let At be a selected step-size for the evaluations
of AH (At may be chosen to be the same as the integration

JSr'
step-size or greater), and choose dt' the same size as the

integration step-size.
The following procedure is called the Computational

Procedure B:

(1) Start from r' = rl and t' = to up to r' = N+1, and

t'St —dt' 0

1‘
(ii) (a) At t' assign a value to v, call it v' .

(iii)

(iV)
(V)

1»?

(b)
(c)

(d)

(a)
(b)

85

Evaluate AH t' ) e

JSr'(

Compare AH (t!) with different values of

JSr'
vI'”(t') for fixed t' . Store the minimum
r,(t').
Compare the result of (c) with the Best—up-to-

AfiJSr,(t'), and the corresponding *7

date, store the smaller and call it the Best-
up-to-date, and denote it by AfiJ35,(t'). Also
store the corresponding 39(9).

Increase t' by At, call it t' .

Repeat (ii) for the new t' .

Repeat (iii) until t! reaches t - dt'

(a)
(b)

(e)

(d)

f 0
If AITIJsi-‘Jtd is zero, then make no change of v.

If Afi (t!) is negative, then add a new stage

JSi'"
to the system after the r'th stage with

TEE-H1 = “’5

Change the number N accordingly and again call

,, trt+1]=[t"t'+dt']

it N, and reorder the switching sequence as

t(N) = {t1’...,tr', tr'+1,ooo’tN+1}
Relabel the system equation and other functions

A Computational Procedure for Algorithm II

The overall computational procedure for Algorithm II

can be obtained by repeatedly executing Procedures A and B

and skipping Block (1) of Procedure A:

(i) Guess and store the nomimal t(N), u:(t), v':(t). 45*,

and b*.

86

(ii) Execute Procedure.A (skip (i) of A) for a given
number of times.
(iii) Execute Procedure B once.

(iv) Go back to (ii).

v NUMERICAL EXAMPLES

Three examples are given in this chapter to illus-
trate the application of the computational technique: (1)
Brachistochrone problem with inequality state constraint; (2)
Orbit transfer of a solar sail ship; (3) Low thrust trajec-

' tory optimization problem. In each of the examples, alter-
native choices of dC and automatic schemes for updating dC
and (11 are discussed.

A FORTRAN program for Algorithm I is presented in
APPENDIX A. The program is written in such a way that for
each specific prolbem one has only to rewrite the subroutines
corresponding to the given problem; however, the major part
of the main program will remain the same.

The computer used for the preparation of the
numerical data is a CDC-6500 digital computer.

5 .1 Brachistochrone Problem with inequality State Constraint

The classical brachistochrone problem was first
proposed by Johann Bernoulli in 1696. This problem can be
stated as follows: Find the path of a heavy particle falling,
under the influence of constant gravitational field, from a
given initial point to a final point on a specified vertical
line in minimum time. The initial point is not contained in
the vertical line.

87

88

We shall modify the classical problem so that the
path of the particle will be restricted in a given region.
The motion of the particle can be described by the differen-

tial equations:

5:1 = x3 cosu = fl (5.1.11!)
5:2 = x3 sinu = f2 (5.1.lb)
i3 = g sinu = f3 . (S.l.lc)

where x1 is the horizontal displacement (feet), x2 is the

vertical displacement (feet, directed downward), x3 is the
velocity of the particle (feet/second), u is the path angle in
radians , that is, the angle between the horizontal axis and
the path at time t. We shall call x = (x1 x2 x3)T the
state variable (vector), and u the control variable. Here g is
the gravitational acceleration (see Figure 14.). The constraint
is

¢=x-%~x-lSO (5.1.2)

This problem has been studied by many authors, for
instance, Dreyfus [DR-l],Denham[DE-2], Bryson et al. [BR-2],
and McGill [MC-1]. In order to compare results with other
authors, the following data are used and they were used in

[DR-l]:
to = 0

29(0) = 23(0) = o

89

 

 

' 11
x3\\ \<!(///r-CONSTRAINED PATH x

\x.
\\. ¢
‘\\

IA

3
\V
O

CONSTRAINT BOUNDARY
¢ = 0

Figure A Brachistochrone problem

 

 

 

 

90
x3(0) = l feet/second
x1(tf) = 6 feet

g = 32.2 feet/second/second
The cost functional J is the terminal thme
J = t - (5.1.3)

and the stopping functional is

sat) = x1(tf) - 6 (5.1.10

where S = 0 is the stopping condition which defines the
terminal time tf . Instead of minimizing J directly, we
shall minimize the penalized cost functional 3 defined by
stf
J'+ / wP(¢)dt
0

3
t

f
0

where w is a positive weighting coefficient, and P(.) is the
penalty function

P(m = (NZHM) (5.1.6)
H(.) is defined by
Hm) = 1 if {a > o (5.1.751)

H(¢) = 0 otherwise (5.1.710)

91

In order to apply Algorithm I, we identify the following

a“:f .x(tf)) = 131. (5.1.8)

L(t,x,u) = MMZHUD) (5.1.9)

9°(tf ,x(tf)) = x1(tf) - 6 (5.1.10)
and

M°(t,x,u) = 0 (5.1.11)

Applying the formulas derived in Chapter III formally, one
obtains the following adjoint equations:

' i} = w¢H(¢) A}(tf) = 0 (5.1.1220
i2 = - 2w¢H(¢) A2(t ) = 0 (5 1 12b)
J J r O C
i3 = - A1 cos u- A2 sinu A3“: ) = O (S l 12c)
J J J J r ’ °
' 1 _ 1 _.
AS - O AS(tf) "‘ 1 (5010133)
'2 __ 1 -
AS — o 13(tf) — o (5.1.13b)
lg = - A; cosu- A: sin 11 Ag(tf) = O (5.1.130)
and
_ .. ' +1.
"Js " >‘J )‘S(8'é_")*,tf (S'l'lb’)
_ T
HJS - 1J3: + L (5.1.15)

92

where
‘é’m, = 1 (5.1.16)
0 = 3
(3)-refit. (x cos u)*,tf (5.1.17)
r = (r1 r2 r3)T (5.1.18)

The constraint on the effort space, the predicted optimal
improvement in 3‘, and the optimal variation in u are,

respectively
t
2 r 2
(dc) =f0 U(t)(bu(t)) dt (5.1.19)
' _ 3: (5.1.20)
dJ - - (dC)(UJJ)
-1 8Has -35
511“?) = - U (tum—)(dCHUJJ) (5.1.21)
more
~tf -1 BHJS 2
UJJ = ./O U (t)(-—55-) (113 (5.1022)

Before one can obtain a numerical solution of the
problem, one has to face a number of difficult problems,
such as, how to choose the iterative step-size do, the
weighting matrix (scalar here) U, the weighting factor w
for the penalty function, and the integration step-size.
There are no straight forward answers to these questions:

they can only be found by guess and experiment. Since dc

93

is a measure of control effort, it has direct control on

the progress of the iterative process as well as the

linearization error. Therefore, the value of (10 must be
updated according to iteration number and error test. It

is possible in some iterations that excessive linearization

errors are observed. ‘Ihen one may have to go back to the

previous iteration, and therefore, both new and old control
functions have to be stored to provide suitable back-up

capability. The sample program presented in APPENDIX A

is fully automatic; there is no operator interface required.

In obtaining the following numerical results, the
same integration step-size (.001 second), and weighting
factor for penalty function (100) were used along with
different schemes in updating do and the alternative choices
of the weighting matrix U. An initial guess for the control
function u(t) = (7/6 (30°), t 6 [0,1] was used throughout
the computations.

Case I: (1) d0 = .1; (2) U'1 = 1; (3) In each iteration, a
linearization test was performed. If the prediction
error (due to linear approximation) exceeded the
prescribed tolerance, d0 was then reduced to 70%. (L1)
Subroutine SELCT was called by the main program once
in each iteration, and dC was then reduced to 70% if
1 $ N S 10, and no effect if- N > 10.

Figure 6 shows the terminal time tr as a function
of iteration number N. In Figure 6, N = 0 is the initial

91l-

iteration, 0 indicates that there was no violation of the
constraint boundary, and 0 indicates that the state constraint
was violated in that iteration. After 26 iterations (including
the initial iteration), the terminal time drOpped from..86767
second at N = 0 to .71.).589 second at N = 25. In Figure 5,
several of the descentpaths corresponding to N = 0, l, 2, 25

were plotted. The average computer time was 1.1.389 seconds/

iteration.
_1 t--.5t2 2
Case II: (1) d0 = .2; (2) U = a+ (b -a)(-—-‘-5-£-——) , 136 [0,1],
t2 = l, a = 2/3, b = 3/2; (3) Same as thgt in I(3)
except 95% instead of 70%; (h) Same as that in 1(h).

Withthese changes in (10, U’1

, and the updating
process, the terminal time decreased more rapidly than that

of Case I. As shown in Figure 7, t decreased from .86767

1'
second at N = 0 to .7L(J.(.78 second at N = 214.. The average
computer time was 11.14.23 seconds/iteration.
Case III: (1) d0 = .2; 0‘1 = 3511-? , t 6 [0,11, s = -1.5,
b '- 2; (3) Same as that in 1(3) except 80% instead of
70%; (1).) In subroutine SELCT, dC was reduced to 90% if
1 s N s 10, 85% if N > 25, and no effect otherwise.
After 25 iterations with (1.399 seconds/iteration
average execution time, the terminal time dropped from
.8676? second to .711289 second.
Case IV: (1) do = .2; (2) 11"1 = 1 ; (3) Same as that in

111(3); ((1) Same as that in 111(k).

95

 

 

 

H ommo ( 80.3099 ozoagoonnmwnomah .Ho mzpwm m oasmfim

 

_
_
m ....

pooh . Hun unoaoosammac HapnoNHAom

air—-

Ts

 

 

qeeJ ‘ 2x anemostdstp '[sonaeA

Terminal time t1. , second

Terminal time tf , second

.$

4%

.70

.$

.w

.5

.W

%

 

 

O
O
~r-
O
O
0 0

e 00.0000..o..oeooo
1 I 1 L I
T ' ' ' T

Figure 6 t1. vs. N , brachistochrone - Case I

 

 

o
o
‘ o

” 0 o o
e e e o°e°e°eoo°o°e
i 1 L 4 1
r* I T *T r

0 5 10 15 20 25 N

Figure 7 tr vs. N , brachistochrone - Case II

97

With this combination the terminal time drOpped to
.711221 second after 26 iterations (N = 25), .711207 after 37
iterations, and {Ill-199 second after 14.3 iterations (N = L12)
with an average computer time (1.1125 second/iteration (or
3.171 minutes for 113 iterations on CDC-6500 computer). This
result agrees with that obtained by Dreyfus [DR-l], in his
work, after 31 iterations, the terminal time was .71122 second,
and after a total of 50 iterations and L). corner modifications
the terminal time dropped to .71120 second with 10 minutes of
computer time (IBM-7090 computer).

Figure 8 shows the terminal time tf as a function
of iteration number N, N = 0,1,...,25, for Case IV. Several
control functions obtained during the iterative process
corresponding to the iteration numbers 0, l, 3, 25, L12 are
plotted in Figure 9, the corresponding paths of the particle
are presented in Figure 10. Here N = 0 represents the
initial iteration, N = 1, 3, 2S correspond to the intermediate
solutions, and N = [12 corresponds to the computed optimal
solution.

From the results presented in this section, we can
draw the following conclusions:

(l) The values of do and U and the updating scheme have
direct influence on the speed of convergence of the
Optimization process. As one can see from the results,
Case IV is superior to Case III, Case III is superior to

Case II, and Case II is superior to Case I.

Terminal time tf , second

, degrees

Control function u(t)

98

 

 

 

 

 

 

0
.85 ‘—
080 "" o
.1
O O
.75 P o
1 .‘ 0.. e°o°e°e°e°eoeoo
1. 1
.70 11 1r .L I 1
O 5 10 15 20 25 N
Figure 8 tr vs. N , brachistochrone - Case IV
901-
80 ‘>\
\
7o -1— \
25\\&112
60 ~N\ \\
\\\Q£ \\\\)
50 5- ' \
\\
(+0 "’ \\\1
\
\\3\ \\
30 w- 0
20 “‘1'" \\__.
\\
10 -w-
25\a112
0 1 1 1 1 1 1 1 \ 12 J_ 1
1 I l 1 r l I I l r

O .1 .2 .3 .h .5 .6 .7 .8 .9 1.0 t

1:1: ....-- o ”An'I-nn'I fluent-Inns- hann‘hishnnhnn‘nn .. (Inna TV

99

 

 

>H onwo 1 Scanoam Ononfioounwnomhn .Ho «Spam 0H. oazmfim

 

 

 

poo.“ . HM unoaoesannap Handouaaom

gee; ‘ 8x anemostdstp teens“

(2)

(3)

5.2

100

It seems not advisable to start with a large do in the
experimental stage of a given Optimization problem. In
general, a small (10 may cause a slow convergence;
however, too large a value of (if) may create computational
difficulty, such as poor predictability. Smaller (10 will
usually be associated with smoother tr vs. N plots (compare
Figure 6 and Figure 7).

It is inaccurate to say that the present iteration is
inferior to the proceeding one by an argument simply

> O. For instance,

based on the measurement Jn - Jo

ew ld
in Figure 6, tf(6) > tf(5), however, the 6th iteration is
not inferior to the 5th, since during the 5th iteration,
the state constraint was violated and that was not the
case for the 6th iteration. Furthermore, tf(17) > tf(l6),
in Figure 6, the 17th iteration is again not inferior to
the 16th although in both of these iterations the cons-
traint was violated. This is because the degree of
constraint violation incurred in the 17th iteration is
less than that in the 16th (the accumulation of the
weighted penalty for the 16th iteration was .02).).50, and

it was .00021 in the 17th iteration).

Orbit Transfer of a Solar Sail Ship
A solar sail ship develOps its propulsion energy

from solar radiation. The radiation pressure is inversely

proportional to the square of the distance between the sun

and the sail. By controlling the sail angle, the ship can

101

be driven to any destination in the solar system.

Assume that the orbits of the planets in the solar
system are circular and coplanar. Furthermore, only the
gravitational force due to the sun is considered. the

equations of motion of the ship are:

2
(v ) r
vr = 39 + ao(-r2)2 cos3u - go(-i;9-)2 (5.2.181)
, v v r
v¢ = - 2 —%:9'+ ao(;?)2 sin u 008211 (5.2.lb)
r = vr (5.2.10)

where u is the sail angle; vr and v¢ are the radial and

circumferential velocities of the ship, respectively; r

O,

r, and rp are the initial orbital, instantaneous, and final
orbital radii respectively; go is the solar gravitational

acceleration at the initial orbit; a0 is the acceleration of
the ship due to solar radiation pressure exerted normally on

the sail at r0; v is the orbital speed of the final orbit

(see Figure 11). F
We want to find the optimal sail angle and the
associated trajectory of the ship such that the ship is
driven from a given initial orbit to a specified final orbit
in shortest time.
Tsu [TS-l] studied the problem and found an approxi-
mate analytic solution by setting tr = 0 with constant

sail angle 0, The resulting trajectory of the ship was a

102

 

 

FINAL
ORBIT
1‘
p
INITIAL _.., v SAIL
ORBIT u
vr
SHIP
(1
SUN REFERENCE AXIS

Figure 11 Orbit transfer of a solar sail ship

103

logarithmic spiral. London [LO-l] obtained an exact solution
of the problem by assuming a logarithmic Spiral trajectory
and constant sail angle. Kelley [KE-l] obtained numerical
solutions for variable sail angle for a transfer from Earth
orbit to Mars orbit. We shall, in particular, consider the
case of transfering from Earth orbit to Venus orbit with
matching terminal orbital Speed.

Let vo be the initial orbital Speed, and VI) be the
terminal orbital Speed. One readily finds that

Vp = VOW-1171:)— (50202)

The following data were used for the Earth to Venus transfer

problem:
80 = .592 X 10'"2 m/sec/sec (at Earth orbit)

2

a = .2x10' m/sec/sec (at Earth orbit)

r 2 ”19.6 X 109 m
r = 108.2 x 109 111
v = 29.76 11:103 m/sec

VI) = 35 x 103 m/sec

For computational purposes,the system is normalized by setting

x1 x2 = lO-3v and x3 = 10‘9r. Note that x1 and x2

=-3
10v, (p

1‘

101.1.

are in km/sec, and x3 is in GM (1A9.6 GM==1 astronoumical unit).
With this notation, and the given data, (5.2.1) becomes

22
i1 = 10"6 £531 + .0hh76(x3)'2 008311 - ,132h9(x3)-2

x
= fl (502033)
1 2
x2 = - 21110-6 x x + .OI.(J.1.76(x3)'"2 sinu cosau
x
= f2 (5.203b)
i3 “ lO-6x1
= £3 (502030)
.. 1 2 3 T
where x - (x x x ) is the state vector, and u is the

control variable. The cost functional is defined as

where w = 1.15711. xlo'sis anormalization factor which

1
converts seconds into days, that is, J represents the

terminal time in days rather than seconds. The initial
states are x1(0) = O, x2 = 29.76 and x3(0) = 1119.6. The

terminal constraints are

1

I x2(tf) - 35 = 0 (5.2.5a)

I2

x3(tf) - 108.2 = 0 (5.2.5b)

and the stopping function was chosen to be

s(t x}(tf) (5.2.6)

f) =

105

S = 0 defines tf if the ship is inside the Earth orbit.

The adjoint equations and their terminal conditions

are
1'3 = .. 1§(§§)* 1§(tf) = [0 0 0] (5.2.7)
11$ = - 11(3— 13* flux.) =[g g; 3] (5.2.8)
1% = -fi—(f)*1g(tf)= [l 0 0] (5.2.9)

where f = (fl f2 f3)T , and (g— in.) is the nxn-matrix of
partial derivatives which can be obtained by differentiating
(5.2.3). The following were used in formulating the combined
multipliers AJS and AIS:

‘ _ l

(S)*,tf "' (f )*,tr (502010)
(2),”? = "1 (5.2.11)
((31)emf = (1'2),Mr (5.2.12a)
:2 ... 3 ‘ ,
(0 ) ,tr (1' )*,tf (5.2.l2b)

In the. following computations, a d0 of 200 was used,
the integration step-size was 21600 seconds or .25 day. The
initial guess for the control variable was «qr/1.1..

t “'05 Sat
Case I: (1) 0’1- = - (b- a) (T) , t e [0, 17280000]

106

a = .6666667, b = 1.5; (2) In the subroutine SELCT, (10

was reduced to 95%. and the error correction in I was

d11 -.05 11 if |11| s .7

= -.2 11 otherwise

d12 = -.1 12- if IIZI s 1.082
_ 2
- -.2 I otherwise

each time SELCT was called.

Table 1 shows the results for some iterations:

Table 1 Numerical results, solar sail--Case I

N tf(days) do I1 12

0 190.257 -- .66L132 4.5114112
1 183.708 200.00000 .87103 -1.382hh
5 177.872 162.90125 .60193 - .3015?
9 176.989 132.68A09 .33188 .A2920
20 178.796 75.17072 -.19132 .72elh
35 180.00A 31.96192 -.20613 .31793
50 180.228 16.1989A -.1255l .1077A
70 180.09h 5.80709 -.0507A .01969

From.this table one observes that the terminal time tf
decreases to a minimum at N = 9, then increases to 180.228
at N = 50, and then decreases to 180.09)... at N = 70.
However, the terminal constraint errors are much smaller

at N = 70 than those at N = 9.

107

Case II: (l) U'1 = 1; (2) In the main program, the value of
do was reduced to 90% if the prediction error exceeded
the tolerance; (3) Same as that in I(2).
Table 2 shows some of the results obtained during
the optimization process:

Table 2 Numerical results, solar sail--Case II

N tr(days) d0 11 I2

0 190.257 -~ .66432 -l.51AA2
1 183.861 200.00000 .86967 -1.AO991
5 178.270 162.90125 .59353 - .10070
8 177.695 139.667u6 .A7217 - .01108
20 178.57u. 75.17072 .126A5 .22976
35 179.092 25.h89h3 .03152 .0930?
50 179.7u6 3.33522 .0130? .02197
70 179.850 .1u536 .OOA76 .00273

Comparing these results with thoseobtained in Case I, we
conclude that Case II is more satisfactory than Case I, for,
N = 70, t

= 179.850 days in Case II, and t = 180.0911 days

if f
in Case I, furthermore I = (.OON76 .00273)T in Case II,
and I = (-.050711 .01969)T in Case I. In both cases the
average computer time was about 7.723 seconds/iteration.
Figure 12 shows the initial trajectory, the inter-
mediate trajectory (N = 8), and the computed Optimal
trajectory (N = 70). Their corresponding sail angles (the

control functions) were plotted in Figure 13 for Case II.

HH once 1 mfinm Adam asaom mo nowaoSOOnmmE NH casmwm

:6 02” end 0+2” ONH 00.... cm

108

b

or

ZDm

 

~1-

 

%

 

BHmmo BHmmo
szw>

       

zoEémS neon .. 1 1 1
zoEémaH new I. - 11..

109

HH ommo u aoflnoaa Haem amHom no OHwnm Hfimm ma oasmwm

when . 9 came

 

 

 

OON owH 00H OdH ONH OOH cm 00 0: ON 0
p _ b p _ p — _ _ _
_ a 4. d _ _ — d A A
Jrll
OTIII . K \\ / I
\ lllllhlll I (I 11r-
\\ // O
\ \ / /
/ on
\\ // Jl.
II (1
\ m LV41/
\\ 11
05
PI .1 \\

 

mm-

om:

ms:

or-

on...

mm-

0N1

sooafiep ‘ n crane I138

110

5.3 Low Thrust Trajectory Optimization Problem

We shall consider the orbit transfer of a low thrust
rocket. As in the solar sail problem, circular and ceplanar
orbits are assumed, and all but the sun's gravitational
field are neglected. Only Earth to Venus orbit transfer will
be studied numerically. Figure 111 shows the notation and the
configuration of the problem.

The motion of the rocket can be described by the
following set of differential equations:

° (v9) 1‘0 2

v1‘ = r + — sinu - g ( ) vr(0) == 0 (5.3.18)
. vrv T

v(p = - 2 “Jr + - cosu V¢(0) 7" V0 (SeBelb)
£- = vr r(O) = r0 (5.3.1c)

where r o = 1119.6 x 109 111 (average Earth orbital radius);

v0 = 29.76 x 103 m/sec (average Earth orbital Speed);

T = .56119336 newton (constant thrust); go is the same as that
given in the last section. The mass of the rocket and fuel is

m = me + mt (5.3.2)

where

mo = 679.602 kg (initial mass);

m = -l.012108 X10"5 kg/sec (fuel flow).

111

 

 

EARTH _
ORBIT
Vq) u vr

1‘

P
VENUS _‘
ORBIT 0 BK

‘P
SUN REFERENCE AXIS

Figure 111 Orbit transfer of a low thrust rocket

112

Let vp be the final orbital Speed and rp be the
final orbital radius, then, for Earth- Venus transfer,
vp == 35 )(103 m/sec and rp = 108.2 11109 m.

The problem is to determine the Optimal trajectory
and the optimal thrust angle of the rocket such that the
rocket is transferred from.Earth.orbit to Venus orbit with
matching final orbital Speed in shortest time. Similar
problems have been considered by Lindorfer and Mayer [LI-l] ,
Kelley [XE-2], Kelley et al. [KB-3], Moyer and Pinkham
[MO-l], Kenneth and McGill [KEN-l], McReynolds [MGR-1],
Sage [SA-l] and others.

If we use the same normalized variables as we did
in the last section, and substitute the given data into
(5.3.1), we obtain the normalized equations of motion:

 

. 2 2 -7
X1 = 10-6 LES)- + 8031 X10 T 81111]. _ .132u9(x3)-2

 

x 1-1.l(-9 x10' t (5,3.3a)
= fl
i2 = _ 2x10-6 x1x2+ filXIO-ZU cosu
x l-lJ-l-9 X 10 t
= r2 (5.3.3b)
i3 - 10"6 x1
= r3 (5.3.30)

x(0) = (0 29.76 1119.6)T

113

where x = (x1 x2 x3)T is the state variable (vector), and

u is the control variable. The cost functional is

J = wltf (5.3.11)
where wl = 1.157(1. 1110-5. The terminal constraints are

I1 = x2(tf) - 35 = 91 = 0 (5.3.521)

12 = X3(tf) - 108.2 = 92 = 0 (503-513)

the stopping functional is

- x1(t = 9° _ (5.3.6)

S(tf) -

f)

where S = 0 defines the terminal time tf if the rocket is
inside the Earth orbit.
The adjoint equations and their reSpective terminal

conditions are:

‘T _ Ta T _

xJ - - 1 Jr) :46) 13%,.) - [0 0 0] (5.3.7)
' 'T_ T6 T 0 1 0

AI — - 11(3—126) 1I(tf) = [0 0 1] (5.3.8)

.T - —

AS - nag—(51311,) -- [l 0 0] (5.3.9)

where f = (f1 1':2 1:3)T and (3%; is the partial derivative
'1"

matrix evaluated along the trajectory of the rocket. The

O . 01 02
quantities (Shut: , (3)41"tf , (e )*’tf , and (e )*’tf ,

11h

are the same as that given in the last section.

In the previous two examples, we used different
values of dC, and different updating schemes to obtain
different speeds of convergence. In this example we studied
two cases with everything identical except different initial
guesses of the control function 11. For both cases, do was
reduced to 70% if the prediction error exceeded 2, and the
present iteration would then be discarded, the program would
regenerate the trajectory of the immediate preceding iteration,
and then continued for further improvement. In the subroutine

SEICT, do was reduced to 97.5% for each 1 S N ff S 15, and no

effect otherwise. Where N6 is the effective iteration

ff
number which differs from N by 2 times the number of iterations

discarded preceding the Ne th effective iteration. The same

ff
scheme for choosing (11 was used as that in I(3) of Section 5.2.
Case I: u(t) = 2w/3, t E [0, 17280000].

Table 3 shows the numerical results for this case.
One Observes, from Table 3, that the terminal time tf increased
from 211.732 days to a maximum of 2118.993 days (at N= Nerf: 18),
however the constraint 12 was driven from 9.65560 to a much
more satisfactory value, 1.12011. The terminal time then

decreased to 202.370 days at N = 57 or Ne f = 53 with

1 = -.05333 and I2 = .030115.

f
I

Table 3 Numerical results, low thrust--Case I

115

N N eff t 1.( days) SC. I: if

0 0 211.732 -- -.073h9 9.65560
1 1 212.806 100.00000 -.06891 9.11231
3 3 215.677 95.00000 -.05A62 8.05579
5 5 219.369 90.25000 -.O3620 7.03335
7 7 223.707 85.73750 -.01318 6.037u1
9 9 228.5h2 81.A5062 -.01AA5 5.06257
18 18 2h8.993 69.83373 .1191u 1.12011
25 25 237.318 69.83373 .09829 .h823h
30 30 228.980 69.83373 .09199 .28691
no 10 21u.lh7 69.83373 .0869? .102A1
50 50 203.369 69.83373 .0311711 .ouoou
57 53 202.370 38.21853 -.05333 .030u5

Several trajectories are plotted in Figure 15, the

corresponding thrust angles are plotted in Figure 16. Here

N = 0 corresponds to the initial iteration, and Neff = 53

eff
corresponds to the computed Optimal trajectory.
Case II: u(t) = 511/6, t E [0, 17280000].

The numerical results for some iterations are given
in Table (4..
Figure 17 and Figure 18 Show some of the trajectories

and the correSponding thrust angle of the rocket reSpectively.

H 03.0 .. poxoo.“ pug» 33“ .Ho mowhopoonufia ma onswg
BHmmo

 

2.0 00H 00H 0+2” ONH OOH cm 00

Pp,_
__fi

BHmmo
mama

 

116

 

.HHmmo

8m 05 omm com 8H 8H oi oma ooH om 8 o: om o

.1111+1+\v\71 +1.

\ /
AUOFI \\\\\. Ax I
\ / ..
\ \\\ x /
\
\ \ /
7 /
Ho: ...
M... Q / \ \ /
/ \ \

\ / l

\ /

x
/ ll
\ z.
/
\ 1
mm?!\\ 1

H omao u poxooh punks» 3oH no mamas unsung 0H onswwm

mhdo . p oefie

 

 

 

OHH

ONH

on."

0...:

omH

00H

0:

00H

owfl

seeaflep ‘ '0 918m datum,

118

Table )1 Numerical results, low thrust--Case II

E 322$ tf(days) SE 1: 12
0 0 168.783 -_ 2.60378 .79175
1 1 170.559 100.00000 2.16572 .75791
3 3 171.602 95.00000 2.18211. .67727
5 5 179.361 90.25000 1.87816 .58379
7 7 185.083 85.73750 1.52581 .17010
9 9 192.192 81.15062 1.05155 .31809
11 11 196.051 77.37809 .61266 .23335
15 15 191.235 69.83373 .52016 .15881
20 20 192.706 69.83373 .39178 .09931
27 27 191.768 69.83373 .26816 .05286
30 28 191.785 18.88361 .25923 .01783
12 10 192.751 18.88361 .16138 .02926
52 50 193.157 18.88361 .0813? .03617

Comparing Table 3 and Table )1, we conclude that to
use u = Err/6 as an initial guess for the thrust angle is more
favorable than that of 2w/3. It is true in both cases that
longer travelling time (tr) pays the price for better satis-
faction of the terminal constraints.

In the computation of this example, the integration
step-size was 21600 seconds (or .25 day), and the average
computer time used was 8.788 seconds/iteration in Case I,

and 7.786 seconds/iteration in Case II.

HH omoo .. poxoon pagan Sod no 373:0...th Sn ofismwm

BHmmo BHmmo
mpzm> mama

 

5 8H. 8H 01: oma 03 oo

 

 

 

119

 

120

HH omen .. poxoon 3.9.39 33” .Ho pawns pug ma shaman
chop . 9 25m.

oom 8H 03 old oma oS om 8 o: ow

 

 

 

o
_ _ _ 1 _ _ _ _ _
1 4 _ A _ 1 _ 4 _
\\.\.l/
\ //
\
O, \ \\\|/ / .1
\ // /
\
\ \ / / ..
/
\ \\ / 11
\ \ /
.STI\\\ \ /
om \ /
fl \ / J
/ x. ,
/ \ / 1
/ \ .

[I

 

om...“

of”

own”

00H

ow."

00H

00H

oom

OHN

seesSep ‘ n 0181:! new

VI SUMMARY AND CONCLUSION

This thesis considers the gradient computational
technique for a class of optimal control problems, the control
problem of Bolza, with various constraints. The major objec-
tive is to derive computational algorithms and their respective
iterative procedures.

Some basic theorems and necessary and sufficient
conditions of the variational calculus, from the simplest
problem of the calculus of variations to a sophiscated control
problem of Bolza, are presented in Chapter II. . This
introductory treatment of the variational theory serves as
background material for the later development of the central
objective of this thesis.

The control problem of Bolza is redefined in Chapter
III as a computational version of the problem. Algorithm I
is then derived in some detail. Iterative formulas are given
in terms of the step-size dC in the control effort space, the
error correction d1 of the isoperimetric constraints I, the
adjoint variables, and the stored nominal solution are of the
system obtained in the preceding iteration. As discussed in
Chapter III, the present formulation of the optimal control
problem is more general and is superior to those formulated
earlier by other authors if computing time and memory require-

ment are considered.

121

122

An important subclass of the general control problem
of Bolza is defined in Chapter IV. Bang-bang control problems,
and problems with discrete control variables as well as conti-
nuous control variables belong to this category. Since it is
so important in application, a Special algorithm—Algorithm II
is derived for this special class of problems, although the
general algorithm—Algorithm I is also applicable. Iterative
procedures for both algorithms are given in detail.

In Chapter V, three numerical examples are given to
illustrate the application of the general computational
algorithm. Alternative choices of the iterative parameters
and their updating schemes are presented in detail. Solutions
of these examples are plotted, tabulated and discussed. A
sample program written in FORTRAN is given in APPENDIX A. The
major part of the main program is machine independent,
will not need to change for individual problems, however, the
user has to write his own subroutines to fit the specific
problem.

As a common character of all gradient techniques,
the convergence of a problem is relatively insensitive to the
initial guess of the problem. The convergence is fast in the
starting iterations, and the speed of convergence slows down
when the number of iterations increases.

Suitable choices of the iteration parameters (such
as the iteration step-size, the weighting matrices, the error

correction (11) and their updating schemes are the key factors

123

for obtaining a reasonable solution of a Specific problem
with a reasonable‘number of iterations. There seem to be no
general rules for choosing these parameters besides reason-
able guessing and conducting meaningful experiments.
Relatively larger integration step-size can be used during
the eaqaerimental stage of the problem so that computer time
can be saved for obtaining final solution. One should not
over emphasize the importance of the measurements of lineari-
zation errors. The suitability of linear approximation
directly affects the degree of predictabilities in dJ and d1;
however , it has no effects on the accuracy of the solution
arc of the system, and the latter can only be affected by the
integration step-size. It is important to note that each
time an iterative solution is rejected, it is necessary to
regenerate the nominal solution obtained in the preceding
iteration, hence one has to pay twice as much computer
time needed for one iteration before any attempt can be
made for further iteration. Therefore, one should not
reject any solution unless it is important to do so.
Throughout this work, a digital computer to carry
out the iterative process is assumed, however, the use of
a well equippedhybrid computer is obviously possible. The
problem of how convenient or inconvenient it will be and
the problem of how much cmnputation time can be saved

against the use of a digital machine are left Open.

121+

The following extensions or further research work

related to this study seem promising:

(1) Formulate a multiple-stage control problem of Bolza and
derive computational algorithms for the solution of the
problem.

(2) Search for new techniques for solving problems with
inequality constraints.

(3) Add further realism to the control problem of Bolza by

incorporating stochastic effects in the system model.

[BE-1]

[BE-2]

[BL-1]

[BL-2]

[BL-3]

[BLU-l]

[BR-1]

[BR-2]

[BRE-l]

[00-1]

[000-1]

REFERENCES

Berkovitz, L. D., "Variational Methods in Problems
of Control and Programming," J. Math. Anal. Appl.,
vol. 3, 1961.

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

Bliss, G. A., "The Evolution of Problems of the
Calculus of Variations," American Mathematical
Monthly, vol. L13, pp. 598-609, 1963.

Bliss, G. A., "Lectures on the Calculus of Varia-
tions," The University of Chicago Press, 1916.

Bliss, G. A., "The Problem of Bolza in the Calculus
of Variations," Annals of Mathematics, XXXIII,
pp. 261-2714., 1932. _

BLUM, E. K., "The Calculus of Variations, Functional
Analysis, and Optimal Control Problems ," in Topics
in Optimization, ed. by G. Leitmann, 1967.

Bryson, A. E. and Denham, W. F., "A Steepest-Ascent
Method for Solving Optimum Pro ramning Problems ,"
J. of Applied Mechanics, pp. 2 7-251, June 1962.

Bryson, A. E. Jr. and Denham, W. F. and Dreyfus, S. E.,
"Optimal Programming Problems with Inequality Con-
traints," I:Necessary Conditions for Extremal
Solutions,AIAA Journal, vol. 1, No. 11, pp. 25111-
2550, Nov. 1963.

Breakwell, J. V. and Speryer, J. L. and Bryson, A. E.,
"Optimization and Control of Nonlinear Systems

Using the Second Variation," SIAM J. Control,

861‘. A, V01. 1’ NO. 2, pp. 193-223, 19630

Coddington, E. A., Levinson, N., "Theory of Differ-
ential Equations," McGraw—Hill, 1955.

Courant, R., "Calculus of Variations and Supple-
mentary Notes and Exercises," 1915-1916, revised
and amended by J. Moser, New York University
Institute of Mathematical Sciences, New York,
1956-1957.

125

[DE-1]

[DEL2]

[DR-1]

[ED-1]

[EV-1]

[GA-1]

[GE-1]

[GU-1]

[HEkl]

[HE-2]

[HE-3]

[HE-1+]

126

Denham, W. F. and Bryson A. E. Jr., "Optimal
Programming Problems with Inequality Constraints ,"
II:Solution by Steepest-Ascent, AIAA Journal,

vol. 2, No. 1, pp. 25-31, January 1961.

Denham, W. F., "Steepest-Ascent Solution of Optimal
Programing Problems," Raytheon Rept., BR—2393
(April 1963)

Dreyfus, S. E., "Variational Problems with State
Variable Inequality Constraints," Rand Corp.
Paper P-2605 (July 1962)

Edwards, R. N. and Brown, H., "Ion Rockets for
Small Satellites," Preprint of ABS Controllable
Satellites Conference, MIT, Cambridge, Mass. ,

1959.

Eveleigh, Virgil, "Adaptive Control and Optimization
Techniques," McGraw-Hill, 1967.

Gamkrelidze, R. V., "Optimal Control Processes with
Restricted Phase Coordinates," Izvest. Akad. Nauk
S.S.S.R., Ser. Mat., vol. 214., 1960.

Gelfand, I. M. and Fomin, S. V., "Calculus of
Variations , " Prentice -Hall , 1963.

Guinn, T., "On First Order Necessary Conditions for
Variational and Optimal Control Problems,"
Ph.D. Dissertation, UCLA 1961.

Hestenes, M. R., "A General Problem in the Calculus
of Variations with Applications to Paths of Least
Time," Eund.Corporation RM-100, 1919; also Astia
Document No, AD112382, 1950.

Hestenes, M. R., "Variational Theory and Optimal
Control Theory, in Computing Methods in Optimi-
zation Problems, edited by Balakrishnan, A. V. and
Neustadt, L. w., 1961.

Hestenes, M. R., "0n Variational Theory and Optimal
Control Theory," J. SIAM Control, Ser. A, vol. 3,
NO. 1, 19650

Hestenes, M. R., "Calculus of Variations and
Opgzmal Control rIllusory," John Wiley 8: Sons, Inc.
19 .

[HE-5]

[HA-1]

[KEnl]

[K12]

IKE-3]

[KEN-1]

[LI-l]

[LO-1]

[LU-1]

[MC-l]

[MG-2]

[MGR-l]

127

Hestenes, M. R., "Numerical Method of Obtaining
Solutions of Fixed End Point Problems in the
Cigulus of Variations," Rand Corp. Rep. RM-102,
1 .

Kalaba, R., "0n Nonlinear Differential Equations,
The maximum Operation, and Monotone Convergence,"
J. of’Math. and.Mech., v61. 8, pp. 519-571. 1959.

Kelley, H. J., "Gradient Theory of Optimal Flight
Paths," ARS Journal 30, pp. 987-951;. Oct. 1960.

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

Kelley, H. J., Kopp, R. E., and Moyer, H. G., "A
Trajectory Optimization Technique Based Upon the
Theory of the Second Variation," Presented at

the AIAA Astrodynamics Conference, Yale University
Connecticut, August 1963.

Kenneth, P. and McGill, R., "Two-Point Boundary-
Value-Problem Techniques," Advances in Control
Systems, vol. 3, edited by Leondes, Academic Press,
New York, 1966.

Lindorfer, W. and Moyer, H. G., "An Application of
a Low-Thmst Trajectory Optimization Scheme to
Planar Earth-Mars Transfer," ARS Journal, Feb. 1962.

London, Howard 8., "Some Exact Solutions of the
Equations of Motion of A Solar Sail with Constant
Sail. Setting," ARS Journal, Feb. 1960.

Lurie, A. I., "Thrust Programming in a Central
Gravitational Field," G. Leitmann, ed., Topics in
Ogfiimization, Chapter 1]., Academic Press, New York,
1 7.

McGill, Robert, "Optimal Control, Inequality State
Constraints, and the Generalized Newton-Raphson
Algorithm," J. SIAM, Ser. A. Control, vol. 3, 1965.

McGill, R. and Kenneth, P., "Solution of Variational
Problems by Means of a Generalized Newton-Raphson
Operator," AIAA J ., .vol. 2, pp. 1761-1766, 1961.

McReynolds, Stephen R., "The Successive Sweep
Method and Dynamic Programing," J. MAA, vol. 19.
No. 3, Sept. 1967.

[mes-1]

[MO-1]

[OK-l]

[PO-1]

[SA-1]

[TS-l]

[VA-1]

[vac-1]

[VAC-2]

128

McShane, E. J., "On Multipliers for Lagrange
Problems," Amer. J. Math 61, 1939.

Moyer, H. Gardner and Pinkham, Gordon, "Several
Trajectory Optimization Techniques, Part II:
Application," Computing Methods in Optimization
Problems, edited by Balakrishnan and Neustadt,
Academic Press, Inc., New York, 1961.

Okamura, Kiyohisa, "Some Mathematical Theory of
the Penalty Method for Solving Optimum Control
Problems," J. SIAM Control, Ser. A, vol. 2,
No. 3, 1965.

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

Sage, Andrew P., "Optimum Systems Control,"
Chzgter 13, Prentice-Hall, Inc., New Jersey,
19 .

Tsu, T. G., "Interplanetary Travel by Solar
Sail," ARS Journal, June 1959.

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

Vachino, R. F., "A Generalized Steepest Descent
Algorithm for Multistage Optimization Processes,"
Ph.D. thesis, 1968, University of Michigan.

Vachino, R. F., "Steepest Descent with Inequality
Constraints on the Control Variable," J. SIAM
Control, vol. 11, No. l, 1966.

APPENDIX A

A SAMPLE PROGRAM

129

“flouz on u H o-vumv “Now. qua

33 «Gouz 6“ u n canvoxv .mom. wk~03
”HOKZ 9" u ~ oa~aoxv «New» Odma
33 OOON u E—JZ

n + <Q<QZ + a~4<a<QZ+mZOUZ.*mZOUZ + «N+mZOUZ.*h<PmZ u YUdmz

.m~0~. h<£a0l

33 kaOZoklwloamP~ZoQQZo<a<QZomZOUZoa0uZoFZOUZokdka auomv NP~G3
. FmJOZ

.muwz oflmbuz .QQZ o<a<QZ owZOUZ .GOMZ okZOUZ okdhmz auomu O<ma

33

wk~a3 do Odwa mamDm mmm oWZMP~ O mzudeOU OQOUMQ IU<M wJaidxm m—Ik Z—
mJ—k NIH Z— WOQOUWE mO OZulkZ omZO~F<a0wkZ~ OQ<3¥U<m no OZu¥U<mZ

momouma mo 02 OM3OJJ< Xdiui—JZ o.m0aouma WDO~>MEQ OZ. mom xmozuuhmJOZ
WZO~F<QMP~ m>~bummmw l0 OZnumMZ oWZO~P<GUF~ mo OZuEwF~Z

kaOU QOW~ mo ha<fl >P~J<DOMZ~ uO inauauz oflMhMZdada JOEFZOU m0 £~OW<G<QZ
mPZ~<aPmZOU U—flhwE~GMQOm~ l0 £~OflWZOUZ omZO~P<a0wFZu 0a<3EOk uO OZuaOuZ

NJm<~E<> JOQFZOU lo £~0uhZOUZ oMJm<~a<> MF<Pm to Z—Ouhdhmz
33

tint-......

.A~.W> ..a._> cauvniwk oaue< oa~.3 ca".
040m canv<hwm oadv3w2m o-.OJO< o.~.<IQJ< onuv3mZ< 6N. wUZwJ<>~DOw
anvfi> canouvn> .Auonvnzmh canauvd oéuonul
oa—VOJOm canvdkmm cauv3m2m o.~.0JO< o~—.<IQJ< caue3mz< ZOnmzwxno
m><aa< DmWDZD

agonumofifiav oaabuymou~.5_uv oganuumoauuunnv MUZMJ<>~DOM

~

U
#**#*U

UUUUUUUUU

a

a

33 OJOum oFmJOZ oklwz oflwhnz 02—42 .2 olPZ N

oho oOk oN3 o~3 elk o.N.Ok .00 .DJODX o3m25x ..NVOqux ..Nvlwzux
oauv3m10 o.~.N.Dm~U oa~o~.3 oa~VOJOD cau~3m23 oaflux ZOEZOU

acvox o.N.1—Q ..NoNunua o.N.~0 o~N~>QIMP 0.
Nv~k3 ..muhiafl oa¢uuw oaoqvmw oa0~o0mvx34 o.0~uam oaO—um ZO~WZNZ~O
anIU<aumZ o.N.IU<flumJ ZO~MZMZ~O
02504 .250< oOu om J(ZEUPXM

Ndmm N<mm Ndmm Ndmm Ndmm

03w: omm to Fame .0243 D m ow_mUIP O 11 20am

(NJOm kO EmeOEQ JOQkZOU MIF com MDO~ZIUMF J<ZO~F<FDOEOU FZM~O<do .

AFDQFDOunmadhohDuzuuNMQ<hobDQPDOoPDQZ~a #00 Z<QUOGQ

Sibanh flag-haunt“; a c Anilpal 4. 14-

n

—

(JL)UIJL)U

130

33
33

Rom wJ4awwkZ— NIP wk43J4>w 0044

FON .mON .mON

.~.um u a~v0m a.

.—

nmzouz. 1—

(x
(\I

.CX u “mum

flouZ .— u ~ «ON 00

k0
Foamw.~ + 0»
Oh

I “flukian
u .va200
".PZGQ OuN

u I mumz n uuwz CON

.4NJON to ZmeOGQ MIP

.0043d0u mzouh430m mh4bm NIP

mh4awwk2~

till-1...!

.0MUZ4IU mm OP meZ OmJ4 >42 000 ¥UOJ0 Zn MP~U3 024 04wa OP
mJJ4U MIh .¥m—0 20 Z4Ih UMIP4E MQ4F ZO OMEOPN mu 4k40 >d0h0w74ak

UIk k—

.PUJNW .0400 .m .02104 .—2304
to Pmma NIP MOZ4IU PWDZ Gum) th
024 DON WPZMZMP4Pm 2mw3kmm .S4EUOEQ W~Ik 024 QUQI UZ—kDOflmDm

33

33

33 N3 .u3

33

33

33

.mimeOGQ 024 mEMFDQZOU Fzmammm~0 aou F42P3 024 .wkua3
.2304 mwznhDOUMDm 024 £4QUOQQ W~IF
.hZ402MQm02— w2~IU4Z wfl4 .000

N**UQOFm n UQOFm
NttNUO u N00

a
n0 .~m

I Z n uFZ mm
.ufl an. un N0

N0 .fin..¢fl «NP I P. u—

aND + Pfiu3v\.

u + 2 u 2

#0 + F N F
MF~Q3 JJ4U

a u .~.~.D 40

m." + P0\.Ok 0 P. u 2
“an .Oemnum. k4tflOk mm
au~.m.O~mm. b42¢0u 00

~ .3m23 .ND .nD .NP .P «00.0. NP~¢3
~ .3wZD .ND .uD .NF .h “on.N. 04wU an

.000“
.m.o~

k0» h4ia0u m
k0. P42flou N

.1w .~w .k0 .Oh .UQOkW .NUO .m.n. whad3
N3 .u3..fiw .nw .F0 .Oh .UQOFm .NUO
«MZOUZ .~ u ~ .a~.~b3. “0.0. mhua3

“02002 .~ u m .a-~b3.

.N.N. 04w“

.N.N. 04wfl

a¥U4mZ .— u n ..uvmm. «0.0. mk—U3

A¥U4mZ .~ u u ..nvmm.

.N.Nv 04ma

«HOKZ .— u ~ ..uemm. «0.0. MF~G3

$18.31*...

U&)U(JL)U(JLJUL)U&)U

131

33
33

33

33

NIP NP4JDUJ4U ode .m0a43¥04m maOPUN> PZ.0§04 NNEIP NIP NP4GOUPZ.

.3 .4. P<ZP3 JJdU

A.NP.Z.uN0 N>.P.m00

Nm OJDOIN WNU.UP4Z NIP JJ4 e¥m~0 20 ONaOPm m. 3. .3 024 .3 .4

Z. WNNGN>Z. NIP ZUDPNU .3 024 .3 .4 NNU~RP4E 02.PIO.N3 NIP NP4OQD

. u PmJOZ

— + ukNZ n ukmz

OON. .CON. .004 ..N I P2~10\N. u—

0.0 .004 .0—0 .P2.50. k.

..040.X I 3NZ.X. I .0 10 ZEOZ £10k.23 0NPIO~N3 «~P3. NIP m. N

mazapzoo

. 2.21sz .... N
, mom .mon .0on ....>azm» n N. 1.
. ....onoix u ...:wzix. n ._._o .mm<*._.~»3 u .~.>azm»
...>azm. u .2110

won .mon ..on ..~.>azm. a .2110. 1.
.....o.mm<...._as u ...>a:w.

mzouz .. a a mom on

.o u »z_wo

.o u N

con .ooc .ooc .mzooz. a.

0JO~X I 3NZ.X IP~3 .0 to P2N2NNH04 UOK PWNP
00—. .00.. .400 .WN I «30\..0J01X I 3N25X. I 50..mm4. l.
0J05X I 3N27X IP~3 50 k0 PZNENNG04 acu PWNP

MON .004 .004 .PNJOZ. l—
. + GNP—Z I GNP—2

K.

«N00 .70 .XD4 .Ok .1 .uJI. .8012 .00 .N .PEQQ. UUQI JJ4U

.0P 024

.00 .3NZWX .3NZ.X NP4JDUJ4U 024 .OGNNZOZ «N.P2an PNN .0 u N ZNI3
.20.P.DZOU 02~QO0Pm NIP .m mNP4JDUJ4U 024 .20.P4Z&0&Z. >EOPUN14EP
0NE.NNO >24 WquaQ .¥m.0 NIP 20 X NNQOPN Cu .J4>CNP2~ N2.P

IU4N 001 N020 0NJJ4U N. .3NUO .0flQ50 .Pzaa .Z .lJIu .Q0 .0 0P. 01 0030

.P 024 N to ZO.PUZD& 4 04 am mNP4JDUJ4U .00 .N .P. u NZ.PDOQNDN

3N25X u 0405X
...3N2.X n ...040~X
mZOUZ .~ u . NON 00

.0e
00¢
can

won
000

NO”
hon

000
400

non

NON
NON
PON

UU

UUUU

U UL)U LJUI)

U(JL’UIJkiu

...1.0 u N
mzooz .. u . .mo 00

N u 13.....0

.x.1..>*.y...azw. + N u N

«can: .. u y «we on

.0.....a u N

mzooz .. u a mum oo

mzooz .. u . mmo on

N u .1...uzm.

.w.y.<..y....> + N u N

«1412 .. a x mum on

.o n N

41412 .. u w nmo oo

mzooz .. u . nNo oo

. .No .Nne .Nno .mzooz. 1.
.no .oco .060 .41412. 1.

334 .3.4 ...4 NPDQZOU

.¢*.3.4P00*.3...3*...4P00 + .3.....Q u .3.....Q
002 .. u 3 000 00
QQZ .. u a 000 00

N00 .ONO .ONO .QQZ. u.

132

..3 NP4JDUJ4U
33a .3.Q ...Q NP4JDUJ4U

33 .N00 .30 .X34 .02304 .2304 .lJI~ .¥U402 .00 .0 .PEGQ. 0UQI JJ4U
.20.P4£a0u2. >fl0PUN34flP OmanmNO >24

0P2.an 024 ¥0~0 NIP 20 3030 024 30—0 mNEOPm 02304 eJ4>GNP2~ N2.P

IU4N Rom N020 0NJJ4U 0. .3N00 .00Q30 .PIHQ .Z .kJI. .00 .0 .P. 02304 E030
.P 024 0 10 20.P023l 4 04 Q0 0NP4J3UJ4U .00 .0 .P. 2304 N2.P30a030

P0 I u .0.P2¢Q
0P I .N.P2¢Q
P0*..IkP2. + 0P u ...PZUQ

.0. .2304 4440
.0 2. PZ.0304 NIP com 020.P~0200 J4.P~Z. NIP PDQ

...0N " ...00
¥0402 .. I ~ .00 00

.3) 024 ..> .333 .3.3 ...3 044GONPZ.

wmo
0N0

MNO
NNO

.NO
.00
0N0

~00

UUU

UUUUU

UUUUU

133

33
33

33

33

33

33

02002 .. I . 000 00
0.0 GP 00
.C I 040

4PN00 .40 NP3Q200

Nt....0 I ....0

02002 .. I . 4.P 00
.N\N00.PQO0 I N

.00 .000 .fl.P .040. u.
N I N00 I 04a
...—0*...>Q2NP + N I N
02002 .. I . N.P 00
...)QZNP!...3.Q + 30 I 30
.0130

.0 I N

N I ...)QZNP
.3..0*.3.....n + N I N
02002 .. I 3 0.P 00

.0 I N
02002 .. I . ..P 00
000 0P 00

33n¢04d I I 30
.33Q\N00.PE00 I 04d
00h .ONP .ONP .02002. t.

30 0NP0.0N¢Q NP4430440

..0 .N00. P04N0 4440
N**00 I N00 2G3PNa ..0 024 00 PUN4N0

.......n\.. I ...—...Q

00F 0P 00

.I040002 .I040004 .N .02002 ...Q. >2.2 4440
.40 .N40 .OOP .. I 02002. I.
...3>*...>QINP + 330 I 330
40412 .. I . 000 00

N I ...)QENP

...3.4*.3.3> + N I N

4C4QZ .. I 3 0N0 00

.0 I N

4&4Q2 .. I . 0N0 00

N u ...3.Q

.3.3>*.3...EZNP + N I N
4E4QZ .. I 3 0N0 00

.00

((4

e:

0.P

N.P

..P
OHP

00b

ONF

00b

N40

.40
040
000

0N0
0N0

N00
FNO
0N0

UKJU

UIJL)

1311

N I ...)030P mmm
.3.3.0*...3..> + N I N «mm
02002 .. I 3 N00 00 000

000 .finm .flmm .WZOUZ. 0.
O<a*.n.3> I I N

40402 c. I . "HQ 00 .00

.00 .000 .000 .(UCQZ. In 000

40 NP30200

N + ...3020 I ...3020 mNm
...3020 I ...0400
.3.>0ZNP*.3...3 + N I N NNO
002 .. I 3 NNO 00
p .O I N
002 .. I . ”N0 00
...3.0#...4P00*.N I ...)020P .N0
002 o. I . .N0 00 0N0

P2N2NP4P0
N02N44>.30N N03 0P.N030 N0 .PZNJ4>.3GN N04 ...3020 024 ...4PNO

4PN00 NP30200

0N0 .000 .000 .002. 0.
N I ...3.0 0.0
.3.>0INP*.3.....0 + N I N N.0
02002 .. I 3 N.0 00
.0 I N
02002 .. I . 0.0 00
N*040 I 30 I 30
040*...3.0 + ....0 I ...)02NP ..0
02002 .. I . ..0 00 0.0

3.0 2. .040t3.0+.0.*...0.>2. NO0P0

.040.P000 I 040
N\040 I 040
...3.0*...>02NP I N I N 400
02002 .. I . 400 00
330 I N
N I ...>02NP 000
...3...0*.3.3.0 + N I N N00
02002 .. I 3 N00 00
.0 I N

UKJUIJ&)U L’UIJ

IJLIU

135

002 .~ I ~ GHQ“ 00 000”
OOOH .—~0_ .—~O— .002. 0— mmofi
.~VCJO< I «~03024 000"
40402 .~ I ~ 000— 00 FOO~
FOO~ .0N0~ .000“ .(04020 0n 000"
A~.OJO~X I «~0302~X 0C0"
02002 .~ I n 000— 00 coca

000" .000— .000— AwZOUZV 0—
OJO3X I 3023X N00—
00d OH 00 N00“

N00~ .0000 .0000 .UQObm I NUO. 0~
NUDt—0. I NUO 000*

O—N o—oo ._oo .ukz I 2. 0— 000
- ~+ZIZ

0F~03 JJ<U
.2 000000 OFZ~ 3030 024 .30~0 .3 .0JO3 .3023 .X P30F3O

N + ...:mzn u ...:wz: moo
.~.3m23 u ...CJOD
.1.>a:wh*.fi.~.3 + N n N coo
kzouz ._ u I «co oo
.o u N
.2002 ._ u . moo on
N u ...>a:m. moo
.1.fi.a..~.fi.3m.u + N u N «00
mzouz ._ I a N00 00 woo
woo .noo .noo .mzouz. u.
. o<c¢...3mww I I N
#2002 .~ I _ moo oo

0(00 4410 ~00
.2 000000 2000 3030 024 .30~0 .3 .OJOD .3023 .X #302—

u I 2 000
“#030 0F30ZOU

N + ...3m24 u .~.3mz< mnm
...3m24 n ...0404
.w.>nzw...1._.< + N u N «mm
(aqnz .. u 5 «nm on
.o u N
«aqaz .H u 0 mnm 00

U

UUU

136

U

2030.00 Q
0— I 0JI— 0

 

c OF 00
N— I 0JI— N
0
0 .N .0 ——XI—N~k2000*10 0—
——.>000 I ——.0.X34 —
.0 I ——.0—~X34
{—02 .— I — — 00
.0 I —0~P!00
.0 I —¢~b£00
—0~h200 I I
——~F200 I X
. 0 I 041—
—OON.0—.X34 .—OON.>000 .—OON.> .—00P£00 20—0202—0
0
N) 000 000\0 20— 2000 00—0—00! 0
0
33 —3000 .330—
33 .X34 .0k30 .k00 .0JI— .2—02 .>000 .> .FZ00. 0001 02—P30003m
mam
h—X0 JJ40 0000
00¢ 0F 00
33 ..000020— 00000 .Ih— .0.0—0 0
33 .043 002400J0b 01F .I0— .0.0—0 .0— 00000 2000—23 01% I—N N
33 .\ .0.0—0N .I —0 J43P04 0Ik .Ib— .0.0—0N —
33 .I —0 00k0—0000 0Ih .—0 2— 00000 20—h4N—0402—J 10¢. #42000 —0N—
33 —0 .—0000 .>0!0h .—0 ——0N— .0. 0P—03
33 ——.0J0—X I ——~302—X I ——e>020F O—N—
33 02002 .— I — o—N— 00
33 h2—30\N I —0000 OCN—
33 000— OF 00
33 —0.0—0 —
33 .043 002400J0h 01h I0— ..30 2— 00000 20—b4N—0402—J 100. #42000 —0——
33

um ._o__.n. m._a3 co.“
oo~.o. cm
0 n .mJoz
m».a3 44<u «.0.
.~.oqo: u ...:wza n~o_
.2002 .0 a . nqo. oo
o<wa 4440
a .. u z ._o— oo
.:.Jz.mm<_ u 5 —.o.
...OJom u .~.3mzm 0.0.

137

33

33

5.02 .0 u a on 00

MN 0# 00

..~.#~XD4 + #J00 + .~.00XD4.*I*MMMMMHM. + A~.~VXD4 I a->
#J00 + #J00 I #400

“N.0.XD4 + .~.0.XD4 I #J00

2—02 .~ I ~ 0N 00

com .om .bm .n t z. 0"
._.>000 I .~.¢+ZVXD4
.u.) I .~.Z.XD4
2.02 .d u _ cm 00
¢0N .com .mm .c a z. 0—
c .cm .0 ..m.hzaa. 0—
03N00;.330 .#$0Q .2—02 .0JI~ .>000 .> .x. Q#30 JJ40
.>awo .> .x. #00 JJ<U
g + z u 2
I + X I X
..~.>awO#boooouco. + ._.h.x34*mnmnmom. I
._.0.x34*#000—0#. + ~_.m.x34*m#n..tl + .~.~.x34.I ._.>
I . ._.>amo u .~.m.x3<
$.02 .— u u NM 00
._.Hzan u x
.>000 .> .x. #00 4440
_ I Z
*ttttttttt**#**#t*.tttttttttttttttttttt mmDJ<> 2301 01# 02.000

032~hzou

..~.>000 + .~.¢+2.x34.*1*m. + .—.Z.XD4 I .->

{—02 .u I _ No— 00

.>awo .> .x. bum JJ<U

I + X I X

.~.¢+Z.XD4*I + .~.Z.x34 I guy)

x_oz .— I _ «a. co
t.**¢***t*****ttttttt#tt*tt****¢#******tttttntttmtttttt** 230!

.->000 I .~.¢+Z~XD4

an.) I .~.Z.XD4

S~02 .u I a 0 00

4 .0 .Q a.0.#200~ 0—

a3N00 .310 .#i0n .2~02 .0JI— .>000 .> .X. 0#30 JJ40

N— .—— .uu .u | Z. 0—

~>000 .> .X. #00 4440

n .— I Z 0— 00

In+.—.#Z00 .IN+.a.#Z00 .I+.«~#200 ..~.#Z00 #4 >000 024 > 0#4DJ4>0

ON

#0

ON

EN
0N

0N

NN

—N

on

NO"

uOu

Nu
«d

138

.mv#200 .~¢.00 ..v00 ZO—0202~0

 

33 .3N00 .00050 .#200 .2 .041— .00 .m .#v 00 02~#DO0mDm
020
00N .00N..N~N ..Ivmm4tu. I ..N.#I0QIX00040 0~ G—N
. N—N .NuN .Q—N aaaN.#Z00IX.*I. 0— m—N
ZOD#00 NuN
NuN .M—N .NuN aa0.#200. 0—
33 .3N00 .370 .#200 .2~02 .04I_ .>000 .> .X. 0#30 4440

.>000 .> .Xv #00 4440 0~N

#400 I .4~#200 mum
onN .0—N .muN .#400 I ~¢~#!00. 0— .
00000 440040 0I# >40#4Z~X00004 mm .¢.#I0n

- ~.~.0~.XD4.004¥.~.0.XD4 + #400 I #400 OGN
Suoz .u I n 00N 00
.0 I #400

a~.0—0X34*#~0004#0. + #400 I an.) CON
#400 I .~.0~.XD4 I .~.0-X34
. XXX~.F.XD4 I a~.0.XD4 + u
-.0.X34 + .~.>000.*I*.n + .~.N.X34 I .~.¢.X34*.O.*0Nu. I #400
2~02 .u I.~ QON 00

«>000 .> .X» #00 4440

#400 I «n.0—0XD4 PCN
-.0u.XD4*00—00N0. I #400 I -0>
..~.0_XD4 + .~.0.X34 + u
.~.#.X34 I .~.00X34 + .~.0.XD4utltnnnnnn.u + .~.—.XD4 I #400
. 2~02 .— I — FON 00
I + X I X
anv>000 I «n.00XD4 MON
.-> I .~.¢.XD4
2—02 .u I a MON 00 ION

.~.Z.XD4 I a~.—IZ.XD4 MON
i~02 .n I n MON 00
0 .N I Z MON 00 CON
**#****t******#*tt 00#000000 00hU—0000 omuuuooz m02—2241 hadkm

nN 0# 00

.a~.0.XD4 + #400 + .~.m.XD4.tI*0#fl. I a—.—.XD4 I an.) 00
#400 + #400 + #400 I #400
.~.#.XD4 + ._.0.X34 I #400

139

33

N. .44 .- .00~0 0—

.44.4 .Xm .m.0~0N .0.0~0¢ .X0 .0~ .4”. #42000
#040~ .3000 .I .40000 ..0.~I~..~.m. .0301— .>40~ .0.0. 0#—03
+I— I X1.#040~

n .NN I 1 on 00

«In I ~

ON 0# 00

II~ I .§~#040—

0N .— I 3 F— 00

n+~ I n

0“ .ON .0~ X—N I a. 0—

*I— I .~.#04Q~

.Iu I .~N~#040~

I— I .5.#040~

, «4 .u I 1 4a 00
«4 I ~

4N .0N .0N an I N4. 0~

— I ~

NN .uN .uN an. 0~

0.~N + 00000000.13000 I ~

0# mON.#0*.—.3023 I 3000

Oh 00N.h0*.4~0 I 40000

a0.0\~—00.¢0N I .N#*aN.0 + N**.~.0.¥~00. I I
m. + ## ### ##N 000.*..00¢ 00*.>40-#4040 I #0 I 030I—
00 000 0N0. + #0 ¢F0 «no 000.*# I >40—

.000 N#— + 0!~#0 I 0I—#0

0.0.4 «wt—#QI#~ 0—

0#~03 4440

«0.0 I .0ox

«N.0 I .NuX

.«00 I auvx

0 I 00—

bu

mu

4—
IN
0N
NN
«N

.\ .nz. .xo .ux. .xo c
..I. .xo .nz. .xo .nx. .xr .mJoz<Im .xo .>uaw2mxo .xo. .m402<1m n
.xm .m3.o<azo .xn .4w> a.uI> .xn .4m> o<axr .xm .ax ><o Ir m
\.oo+xn .xr .on+In .x0 .01. .xm .onaIn .xb .oocxn .xm .JOapzoth .
.xo ...mmoxm .x.~ .mm4m<.a<> u»<.mrm. .xo. .us.. .:0. p<xaou
.n~.n. mh.a3

.ou mz.pa

m.... ...z. I.

m.. + »o\.o» u h. u z
040.m ..m402 .uuwz .am».z .2.42 .z .0»: w

..o .0» .m3 ..3 .uk ..m.o. .oo .OJowx .3mzwx ..~.040.x ..~.:mz.x .
....amwu ....m.3m.o ......3 ....0403 ....3m23 ..n.x zozzou

.... FOJQ. zo.mzm:.o

114.0

33
33
33
33

33
33

33
33
33
33
33
33
33
33
33
33
33
33
33
33

33
33
33

.0 I ...m

0...I. . 00

.0..0 Z0.0202.0
«0..2304 02.#300030

 

020

ZOD#00

...0 I 040.0

«0.0.0.I 30 2. 00000 0>.#4400 .I00.m.0.0.I 0003 00 I4.0

.\.0.0.0.I040N. .I0.0.0.0.I040.. .I0.0.0.0N
.I0403 .IF .0.0.0 .I 30 44D#04 .I0..
.0.0.0 .I50 00#0.000n I4. .0I# .IF .0. .I0#Z I0. #42000
30000.NNN00.«N.040.X....040.X.0401X.1030#.00010.0#Z .0F.0.0#.03 00
00 .00 .00 «.0000 I 10000. 0.
.00000\«1030#I00030..004I50000 00
00 0# O0
.0000 I.10000 00
00 .00 .00 «00050. 0.
0405XI3025XI3030#
FF .F .F «#0402. 0.
. «0.0.0N
.I3WZN. .Im .m.0.0 .I302.. .I0.0.0.0 .I30Zfi .IF .4.0.0 .I4#0I# .I0.
.0#40.I0.0.00 .I000I4.0.0.0.I0# .IM.0..I0Z .I0.0..IZI10.#42000 0.
302.X .3021X .40000 .>400# .0# .0002 .00#.Z «0..0. 0#.03
0F 00N.F0*«0# I #.*.4.00 I 40000 I 40000
«0.0.0N .0.0.04 .XN .0.F0 .X.. #42000 00
3000 .I .40000 ..0..I.....0. .0#40 «00.0.0#.03
.3N00.#000INN000
...00\«0.00 I .N.0#
...00\«N.00 I ...0#
«..00\.«4.00¢N3 + .3. I 00
N.00. I «0# I #.¢.0.00 I .0.0 I «N.30Z.X
0.00 I «0# I #.#«N.00 I .N.0 I «..30Z.X
N3*««0#I#.*«4.00 I «4.0. + .3*0# I 3025X
0000. + F0 4F0 ..0 000.## I 0>40
0000. + F0 4F0 ..0 000.¢0#I>400#
«040.0I...0.\«..0¥.0.#200 I # I 0#
OF 00N.F0*«..3023 I 3000
0F 00N.F0*«4.0 I 40000
«0.0\..00.40N I .N**«N.0 + N**«..0.*.00. I I
Z I 0#Z
.00. I «0.#200
0 .0 .F ««..0. 0.
F 0# 00
. I 00.
F .F .0. «...m. 0.

OF

FF

1&1

33

33

33

33

33

33

n##«n.X\
.m**3mou*mm 000. I 00 400.. + N**..M.X\.N.x.#.oo ooo.I
.n.X\.N.x*Noo ooo.
.....3023.mou.mm< DmOU
....3023.Z.m 32.0
0400 JJ40
. I Z
N .. .. «2. 0.
0.. + #O\«O# I k. I Z
.n.x0.NDO. ..N.x0..30. 00204<>.300

0.x0
N.X0

040.0 .#0J02 .0002 .00#.2 .E.JZ .Z .0#2 N

.ko .0# .N3 ..3 .uk ..N.oh .00 .OJOVX .3023x ..N.OJO.X ..N.302.x
....Dmdo ....N.Dm.0 ......) ....OJOD ....3023 ..n.x zozzou
.0..nm ..o..m zo.mzmz.o
XJ J<0a
.am .0 .h. Zfio< 02.#Docmbm
020
ZOD#00
.n.m\.N.mt.oo 000. I .c.nm
«..0#.oo 000. I .n.00
DzumthmDaIk +..n.m\.N.m*...thoo ooo.II .N.nm
N**.n.m\o¢ Nn.. I DmOU#hmDmI# + .n.m\N**.N.m*.oo 000. I ...am
>P.><ao hmbark m»<z.oaoou adJOQ .
N**..n.m\3wou.tor 440. I hmaaxh

 

024.040 .0JUZ4II«4.0

.#.ZD J4U.SOZO0#04 . I 200.04.. 20 .03.040II«0.0
0\£¥ .>#.UOJ0> J4.#20000ZDU0.UII.N.0

0\!¥ .>#.UOJ0> J4.040II«..0

.....zwza.mou.mm< u amou
“....3wza.z.m u Dz_m
o<ma 44<u

N .u .n .mp2 I z. u.

m.. + »0\.op I h. a 2

040.0 .#0JOZ .0002 .00#.Z .Z.JZ .2 .0#2 N

.po .0# .m: ..3 .m» ..~.o» .ow .0401x .3mzax ..~.oJo.x ..~.3wz_x
....meu .._.~.3m.o ......3 ....0403 ....3wza ..n.x zozzou
...nm ....m zo.mzmz_o
Km 5 .5... 3:30am?»
ozu
zaahwa
.u u .o..m
._ u .o.m
.— u .n.m

U

«(U U‘)L)U()U

11.3.2

33
33

33

33

 

«....30.0#«....3*«..3030 I I «F..00
«3030.#*«3.>Z.*«30.0. I I «0.0.0
«..N.30.0*«....3*«..N.30.0 I I «0..00
«4..00 I «0..00
«....30.0*«....3*«..N.30.0 I I .4..00
«....30.0*«....3#«....30.0 I I .0..00
«30.0.#*«3.>Z.*«30.0. I I «..0.0

030*««...0#00 I «0.0. + .30*««0..0*00 I «F.0. I «..3030
«3 #03 0 44.#040.*««0 400244.*«00. I «0 400244.. I «3050.

NDO*««...0*«N.0# I «4.0. + .30*««0..0*«N.0# I «N.0. I «..N.30.0
030*.«...0*«..0# I «0.0. + .30*««0..0¢«..0# I «..0. I «....30.0
«3 #03 0 44.#040.*««0 400244.*«0#. I «. 400244.. I «30.0.

N**.n.X\DmOU*.N**DZ.m*.n I ...tOF 440. I N30
NII..n.X\DmOU.IDZ.mImN 4n..I I .30

nmxmt....mI n.x0*.o..ml I.N..am
NNXII....0I N.XOI«o..ml I....nm
.nxm*.N..mI .qut....ml I«o..am
.x has 0 J<.ha<u.*.m <00£<J. I I «m <omz<4.o

x4 I anu*«m.m I n.x0*«b.m I I .0.00
.00 ooo.*~**«n.X\«N.x¢N3 I I x4

x4 I Nqut.m.m I N.X0*«F.m I I «0.00
.00 ooo.*«n.X\N3 I x4

.nx0*«o.m I .Nx0*«m.m I I «F.nm

«X #03 4 44.#040. I «x #03 0 44.#040.¥«3 400t44. I I «3 400244.0

anu*.¢.m I n.xu*.~.m I I .o.nm

nmxut.n.m I n.xu*...m I I .m.nm

memt.c.m I N.X0I.N.m I I .¢.nm

Nqut.n.m I N.xu*...m I I «n.am

.nxu*.o.m I .qu*.¢.m I I .N.am

.nxu*.m.m I .qu*.n.m I I ...nm

.x #03 z 4<.#0<u. I «x #03 0 44.#0<0.*«. 400244. I I «. 400244.0

.00 000. I .nx0

n**.m.X\N#*300U*DZ.0*N0 000. I N#*«n.X\«N.x*«..x*Noo 000. I nNXO
«0.X\«..X*NOO 000.II NNXO
N.X0II .Nx0

1&3

«2.0. 0300 I «....3
«2.0.0300 I «..0403
«2.4.0300 I «..3023
«2.0.0300 I «0.x
«Z.N.030U I «N.X
«2...030U I «..X

0 .0 .4 «2 I 2.42. 0.
..N.N «2.0.

n

 

 

 

 

33 040.0 .#0402 .0002 .00#.2 .2.42 .2 .0#2 N
.#0 .0# .N3 ..3 .0# .«N.0# .00 .0403x .3023x .«N.040.X .«N.302.X
....3010 .«..N.30.0 ......3 .«..0403 .«..3023 .«0.X 202200
33 «000N.0.030U \.\ 202200
0400 02.#300030
all
- 203#00
«3 .4. #42k3 02.#300030
020
33 203#00 4
03 0N00.*N00 I N00 0
33 0 .4 .4 «.I0002.0. 00
33 N.*«N.302.X I I «N..0 on
33 00 0# 00
33 ..*«N.302.X I I «N..0 0N
33 on .0N .0N «N00.. I ««N.302.X.004.0. ON
33 N.*«..302.X I I «...0 0.
33 ON 0... 00
33 00.*...302.X I I «...0 0.
33 0. .0. .0. «000F. I ««..302.X.004.0.
33 040.0 .#0402 .0002 .00#.2 .2.42 .2 .0#2 N
.#0 .0# .N3 ..3 .0# .«N.0# .00 .0403X .3023X .«N.040.X .«N.302.X
.«..3030 .«..N.30.0 .«....3 ....0403 .«..3023 .«0.X 202200
.N..0 20.0202.0
33 «.0.N00.#U400 02.#300030
020
203#00
0#.03 4440
«0.#200.«0..00.«0..0 20.0202.0
33 «3N00 .00050 .#200.2.0JI..00.0.#.

 

 

 

«..3030*«....3*«..3010 I
«3030.#*«3.>Z.*«3030. I

«0..00
«350.0

«..N.30.0#«....3*«..3010 I «0..00

 

 

 

 

m >Z.Z ¥I.¥.4
h >Z.z 2+¥ZI¥Z
0 >2.z Z..I¥ 00 00
m >Z.z ZII¥Z
4 >2.2 o..Io
n >2.z #202040 #000044 000 100400

N >Z.2 «..z....4.«..4 Zo.mzmz.o
. >2.2 .2.4.o.2.4.>2.z 02.#30003m
020

0 O# 00

2.42 I Z

2 .0.0. 0#.03

no._.00

. I 2

.00#003000 .I.. .0. ..000000 440044. 10.. #42000

2 «0.0. 0#.03

203#00

«..N.30.0 I .2.0.03mu

.....30.0 I «2.0.0300

«..3030 I «Z.F.0300

«....3 I «2.0.0300

...0403 I «2.0.0300

«..3023 I .Z.4.03mu

.n.x I .2.n.0300

«N.x I «Z.N.03mu

...x I .2...03mu

n .n .4 «Z I 2.42. 0.

. ..N.N .2.0.
33 040.0 .#m402 .0002 .00#.2 .2.42 .2 .0#2 N

.#o .0# .N3 ..3 .0# .«N.0# .00 .0403x .3023X ..N.o4o.x ..N.302.x

....3030 ....N.30.0 .«....3 .«..0403 ....3023 .«0.x 202200

33 «OOON.0.0300 \.\ 202200

0#.03 02.#300030

020

0 0# 00

2.42 I 2

2 «0.0. 0#.03

0 O# 00

. I Z

«00#003000 .I.. .0. ..000000 440044. I0.. #42000
2 «0.0. 0#.03

203#00
«2.0.0300 I «..N.30.0
«2.0.0300 I .....30.0

.Z.F.030U I ...3030

ILLS

mm
Nm
.m
00
04
00
F4
04
03
44
04
N4
.0
0c
on
mm
#0
on
mm
40
NM
Nfl
.0
ON
ON
ON
FN
0N
MN
4N
MN
NN
.N
ON
0.

0.

F.

0.

m.

4.

n.

N.

..

O.

>Z.Z
>Z.2
>Z.Z
>Z.Z
>Z.E
>Z.2
>Z.E
>Z.Z
>Z.Z
>Z.Z
>Z.Z
>Z.E
>Z.Z
>Z.E
>Z.Z
>Z.E
>Z.£
>Z.2
>Z.Z
>Z.E
>Z.£
>Z.2
>Z.Z
>Z.2
>Z.Z
>Z.E
>Z.Z
>Z.£
>2.2
>2.2
>2.2
>Z.2
>Z.£
>2.2
>Z.E
>2.Z
>Z.2
>Z.S
>Z.2
>2.2
>Z.Z
>Z.£
>Z.S
>Z.Z
>2.Z

.+¥2I¥.
z._u. 00 oo
x.a#<2 monoma
032.#zou
.<o_mu.\.x..<u.x..<
.+x2uy_
00.00.00 Ayn—.0.
z..«. 00 oo
203#00
o.ouo
m¢.oc.mc .<u.m.0.
.40.0 2. 002.<#zou
0. #202040 #o>.a 00 034<>. #o>.a 032.2 >0 223400 00.>.o
0401» ..5.<
..1.<u.¥1.<
.xw.<uu04ox
3+01I.fi
1+22uyw
z..u1 cc oo
..-..tzIaw
mn.m¢.mc .xu..0_
«2.2u.
0223400 00242u00#z.
QJOII ..1.<
..w.<u..y.<
w+¥u.¥u.1
..¥.<n-0402
2+.2I.¥
z..u. on 00
ZI¥u.x
mm.mn.0n «2-5.0.
«2.4I5
mzoa 002<Iuaw#z.
032.#zou
wu.y.2
.a.y.4
.3..<I<o.m
om.om.m. «..1..<.mm< I.<o.m.mm< .0.
.+N.I5.
2.2». cu oo
..-1.I2u~.
2.2u1 om oo
.xy.<u<u.m
2+22uy¥
¥I«¥.2

mm
00
04
04

MG

04

00

mm

on

0N

ON

0.
0.

114.6

00

3883

00
NO
.0
00
O0
00
F0

00
#0
00
N0
.0
00
OF
0h
FF
or
0b
05
0h
Nb
.h
or
00
00
#0
00
00
«0
00
N0
.0
00
00
00
#0
00

mm
#0

>Z.2

.>Z.Z

>Z.Z
>Z.Z
>Z.2
>Z.Z
>Z.E
>Z.2
>Z.Z
>Z.E
>Z.2
>Z.Z
>2.!
>Z.2
>Z.Z
>Z.Z
>Z.2
>2.Z
>Z.E
>Z.Z
>Z.2
>Z.2
>Z.Z
>Z.2
>2.2
>Z.2
>Z.2
>Z.2
>Z.Z
>Z.Z
>2.2
>Z.£
>Z.2
>2.2
>Z.2
>Z.2
>Z.2
>Z.2
>Z.2
>Z.£
>Z.Z
>Z.Z
>Z.Z

>Z.2
>Z.Z

203#00
oo. 0# 00
Q40II «.5.4
«.fi.4II«.¥.4
1+¥I.¥I.fi
«.¥.4I04OI
Z+.¥I.¥
Z..I. on. 00
ZI¥I.¥
mm..oc..oo. «¥I3.0.
.¥.2I5
040II «.fi.4
..1.4II«¥5.4
1+01I.1
«11.4IO4OI
3+03I¥fi
2..Ifi 0.. 00
«.I..*ZI03
«.I¥.*ZI05
00..om..om. .21..0.
«¥.4I.
mo..om..om. «2.0.
«.I¥.I¥
2I¥
00241000#2. 223400 024 300 442.0
032.#200
40.0\o..u.¥¥.4
440000.000 >0 #0>.0 0044000
40.0toIo
m#o>.n 00 #030000
032.#200
40.0\«1¥.4I«3¥.4
o#.0b.o# «¥I3.0.
2+5¥I1¥
2..I5 m# 00
20¥Ifi¥
#0>.0 >0 300 00.).0
032.#200
«W..4+«1¥.4*«¥..4I«3..4
1+.I3.I3¥
N0.00.No «2I3.0.
00.00.00 oxl..0.
z+fi.ufi.
2..I1 00 00
ZI.I1.

cu.

00.

00.

ON.
0..

00.

00.

00.

C0

0%
OF

NO
00

11L?

OD

.O¢00..OOOO.

0N00.
0N00.

>Z.Z

0N00.
ONMO.

.0000000F.l

ON ON

0000. 0N00.

0N00. 0N00.
0

..

.OOO.N

0000.
0N00.
OO0ONMNO.

0

.O

0

.0
.0

0000.
0N00.
0N00.

0¢0N00.0.
0.04.
O

N

.OOOO0NF.
.000.
m.
0N00.
0N00.
0N00.
0¢N0.000.
OFUON

4

.
O
.OON
0.
0N00.
0N0O.
0N00.

.00OOFFN.
O

020

0

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

1HIIHIflIHUIMHJH[l1fllfllfllflllllfllfllfllfllflllH