DYNAMICS AND OPTIMIZATION OF A HUMAN MOTION PROBLEM
By
TUSHAR KANTI GHOSH
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE :.UNIk.RSITY OF FLORIDA IN PARTIAL
FULrILL~.I..'-I OF THE REQUIREMENT'S FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1974
DYNAMICS AND OPTIMIZATION OF A HUMAN MOTION PROBLEM
By
TUSHAR KANTI GHOSH
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1974
TO MY FATHER AND MY MOTHER
ACKNOWLEDGMENTS
It has been a very happy experience to work with Dr. W. H.
Boykin, Jr., during my stay at the University of Florida. I am
grateful to him for his help and guidance throughout this work, from
suggesting the topic of the dissertation to proofreading the manuscript.
I wish to express my deepest gratitude to him, both as an efficient and
enthusiastic research counselor and as a human being.
I wish to thank Dr. T. E. Bullock for the many discussions
I had with him about the theory and numerical methods of optimization.
These discussions provided me with understanding of many of the concepts
that were used in this work.
I chose Drs. L. E. Malvern, U. H. Kurzweg, and 0. A. Slotterbeck
to be on my supervisory committee as a way of paying tribute to them as
excellent teachers. It is a pleasure to thank them. Special gratitude
is expressed to Professor Malvern for going through the dissertation
thoroughly and making corrections.
I am thankful to Drs. T. M. Khalil and R. C. Anderson for
examining this dissertation when Dr. 0. A. Slotterbeck left the University.
I am thankful to Mr. Tom Boone for his interest in this work and
for volunteering his services as the test subject of the experiments.
Thanks are also due to the National Science Foundation which
provided financial support for most of this work.
Finally, it is a pleasure to thank my friend Roy K. Samras for
his help in the experiments and Mrs., Edna Larrick for typing the
dissertation.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS iii
LIST OF FIGURES vi
NOTATION vi i i
ABSTRACT ix
CHAPTER 1 INTRODUCTION 1
1.0. Why and What 1
1.1. Dynamics and Optimization in Human Motion 1
1.2. Previous Work 3
1.3. The Problem Statement 7
1.4. The Kip-Up Maneuver 8
1.5. Present Work . 8
CHAPTER 2 EXPERIMENTATION AND CONSTRUCTION OF THE
MATHEMATICAL MODEL 10
2.0. Introduction 10
2.1. Mathematical Model of the Kip-Up 10
2.2. The Equations of Motion 13
2.3. The Equations of Motion for the Experiment
and the Integration Scheme 21
2.4. Experimental Procedure 23-
2.5. Results and Discussion 26
2.6. Sources of Errors 31
2.6.1. Imperfections in the Model 31
2.6.2. Errors in Filming and Processing the Data 31
2.6.3. The Integration Scheme 32
CHAPTER 3 ANALYTIC DETERMINATION OF HIE MINIMUM-TIME
KIP-UP STRATEGY 34
3.0. Introduction 34
3.1. Mathematical Formulation of the Kip-Up Problem . 35
3.2. Bounds on the Controls 36
3.3. Torsional Springs in the Shoulder and
Hip Joints 38
3.4. Boundary Conditions 39
iv
TABLE OF CONTENTS (Continued)
Page
CHAPTER 3 (Continued)
3.5. The Necessary Conditions of Time
Optimal Control 40
3.6. The Solution Methods 47
3.7. A Quasilinearization Scheme for Solving
the Minimum-Time Problem 51
3.7.1. Derivation of the Modified
Quasilinearization Algorithm 52
3.7.2. Approximation of the Optimal Control
for the Kip-Up Problem 61
3.7.3. A Simple Example Problem for the
Method of Quasilinearization 69
3.7.4. The Results With the Kip-Up Problem ... 72
3.8. Steepest Descent Methods for Solving the
Minimum-Time Kip-Up Problem 79
3.8.1. Derivations for Formulation 1 81
3.8.2. Derivations for Formulation 2 93
3.8.3. Derivations for Formulation 3 106
3.8.4. The Integration Scheme for the
Steepest Descent Methods 108
3.8.5. Initial Guess of the Control Function . 109
3.9. Results of the Numerical Computations
and Comments 110
3.10.Comparison of the Minimum-Time Solution
With Experiment 127
CHAPTER 4 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK . 130
4.1. Conclusions 130
4.2. Recommendations for Future Work 134
APPENDIX A DETERMINATION OF THE INERTIA PARAMETERS OF THE
KIP-UP MODEL FROM THE-HANAVAN MODEL 138
APPENDIX B AN INVESTIGATION OF A STEEPEST DESCENT SCHEME
FOR FINDING OPTIMAL BANG-BANG CONTROL SOLUTION
FOR THE KIP-UP PROBLEM 143
LIST OF REFERENCES 149
BIOGRAPHICAL SKETCH 152
v
LIST OF FIGURES
Figure Page
1. Hanavan's Mathematical Model of a Human Being 12
2. Mathematical Model for Kip-Up 14
3. The Three-Link System 15
4. Sketch of Kip-Up Apparatus Configuration 24
5. Measured and Smoothed Film Data of Angles 0 and i|r
for the Kip-Up Motion 27
6. Measured and Computed Values of cp for Swinging Motion ... 29
7. Measured and Computed Values of cp for Kip-Up Motion .... 30
8. Unmodified Control Limit Functions 37
9. Approximation of Bang-Bang Control by Saturation Control 64
10. Graphs of Optimal and Nearly Optimal Solutions Obtained
via Quasilinearization for Simple Example 71
11. Modification of Corner Point Between Constrained and
Unconstrained Arcs After Changes in the Unconstrained
Arcs 90
12. Solution of Example Problem by the Method of Steepest
Descent 105
13. Initial Guess for the Control Functions 118
14. A Non-Optimal Control Which Acquires Boundary Conditions 119
15. Approximate Minimum Time Solution by Formulation 2
of the Method of Steepest Descent 122
16. Approximate Minimum Time Solution by Formulation 3
of the Method of Steepest Descent 124
vi
LIST OF FIGURES (Continued)
Figure Page
17. Angle Histories for Solution of Figure 16 125
18. Difference Between Measured Angle and Mathematical
Angle for Human Model Due to Deformation of Torso .... 129
19. Construction of Kip-Up Model from Hanavan Model 139
vii
NOTATION
Usage
Meaning
X
dx
; total derivative of the quantity x with respect to time t.
X
x is a column vector.
T .T
x (x)
Transpose of x; defined only when x is a vector or a matrix.
x,
y
3x
; partial derivative of x with respect to y.
oy
3x
x, ^
- y oy
Partial derivative of the column vector x with respect to the
scalar y. The result is a column vector whose i^ component
is the partial derivative of the i^a component of x with
respect to y-
9x
xv ¥
Partial derivative of the scalar x with respect to the column
vector y; also called the gradient of x with respect to y.
It is a row vector, whose i^h element is the partial deriva-
9x
tive of x with respect to the i^h element of y.
-V 3y
Partial derivative of a column vector x by a column vector y.
The result is a matrix whose (i,j)th element is the partial
derivative of the ith element of x with respect to the
element of y.
ii ill
T n 2
Norm of the column vector x. Defined as either x x or E K.x.
i=l 1 1
(to be specified which one) where x is a vector of dimension n
are given positive numbers, x^ is the i^h element of x.
viii
Abstract of Dissertation Presented to the
Graduate Council of the University of Florida in Partial
Fulfillment of the Requirements for the Degree of Doctor of Philosophy
DYNAMICS AND OPTIMIZATION OF A HUMAN MOTION PROBLEM
By
Tushar Kanti Ghosh
March, 1974
Chairman: Dr. William H. Boykin, Jr.
Major Department: Engineering Science, Mechanics,
and Aerospace Engineering
Questions about applicability of analytical mechanics and
usefulness of optimal control theory in determining optimal human
motions arise quite naturally, and especially, in the context of man's
increased activities in outer space and under water. So far very
little work has been done to answer these questions. In this disserta
tion investigations to answer these questions are presented.
A particular gymnastic maneuver, namely, the kip-up maneuver is
examined experimentally and theoretically. A mathematical model for
a human performer is constructed for this maneuver from the best per
sonalized inertia and joint centers model of a human being available
today. Experiments with the human performer and photographic data col
lection methods are discussed. Comparisons of the observed motion with
solutions of the mathematical model equations are presented. Discrepan
cies between the actual motion and the computed motion are explained in
terms of principles of mechanics and errors in measurements. Some
changes in the model are suggested.
IX
An approximate analytic solution of the kip-up maneuver performed
in minimum time is obtained for the model via optimal control theory.
Several numerical methods are used to determine the solution, which is
compared with the observed performance of the human subject. Difficul
ties in solving human motion problems by existing numerical, algorithms
are discussed in terms of fundamental sources of these difficulties.
Finally, recommendations for immediate'future work have been
made.
x
CHAPTER 1
INTRODUCTION
1.0. Why and What
Man's increased interest in the exploration of space and the
oceans was an impetus for a better understanding of the mechanics of
large motion maneuvers performed by human beings. Experience in space
walking and certain athletic events brought out the fact that human
intuition does not always give correct answers to questions on human
motion. For certain problems the solutions must be found by analytical
methods such as methods of analytical mechanics and optimal control.
Broadly speaking, this work deals with the application of the
principles of mechanics and optimal control theory in the analytical
determination of human motion descriptors.
1.1. Dynamics and Optimization in Human Motion
Dynamics provides the basic foundation of the analytic problem
while optimal control theory completes the formulation of the mathemat
ical problem and provides means to solve the problem. In any endeavor
of analytic determination of a human motion the first steps are construc
tion of a workable model having the same dynamics of the motion as that
of the human performer, and obtaining the equations of motion for the
human model. However, the principles of mechanics alone do not give
enough information for analytical determination of a desired maneuver
whenever there are sufficient degrees of freedom of movement. W'ithout
1
2
knowing what goes on in the human motor system, optimal control is
presently the only known analytical method which can provide the remain
ing necessary information.
After a workable dynamic model of the human performer has been
obtained, the remaining part of the analytic determination of a physical
maneuver is a problem of control of a dynamic system where the position
vectors and/or orientation of the various elements of the model and
their rate of change with respect to time (representing the "states" of
the "dynamic system" being considered) is to be determined. The state
components change from one set of (initial) values to another set of
(final) values at a later time. The "control variables," that is, the
independent variables whose suitable choice will bring the change are
torques of the voluntary muscle forces at the joints of the various
limbs. However, the problem formulation is not mathematically complete
with the above statements because there would, in general, be more than
one way of transferring a system from one state to another when such
transfers are possible. Constraints are required in the complete for
mulation. The concept of optimization of a certain physically meaning
ful quantity during the maneuver arises naturally at this point. A cost
functional to be maximized or minimized gives a basic and needed struc
ture to the scheme for exerting the control torques. It may be expected,
logically, that, unless given special orders, a human being selects its
own performance criterion for optimization while doing any physical
activity. Some of the physically meaningful quantities that may be
optimized during a physical activity are the total time to perform the
activity, and the total energy spent during the activity.
3
1.2. Previous Work
Very little work has been reported in the literature so far
where the optimization considerations have been used in the study of
human motion. Earlier work in the study of the mechanics of motion of
living beings was done primarily from the view of grossly explaining
certain maneuvers modeling the applicability of principles of rigid body
mechanics. Most of this work was done under either free fall or zero
gravity conditions. The models were made up of coupled rigid bodies
and conservation of angular momentum in the absence of external torques
was the most often used principle of mechanics.
The righting maneuver of a free falling cat in midair attracted
the attention of several authors in the early days of studies of living
*
objects. Marey's [1] photographs of a falling cat evoked discussions
in 1894 in the French Acadmie des Sciences on whether an initial angu
lar velocity was necessary in order to perform the righting maneuver.
Guyou [2] modeled the cat by two coupled rigid bodies and explained the
phenomenon with the aid of the angular momentum principle with the angu
lar momentum of the entire cat identically equal to zero. Later, more
photographic studies were made by Magnus [3] and McDonald [4,5,6].
McDonald made an extensive study of the falling cats with a high speed
(1500 frames/second) motion picture camera. His description of their
motion added many details to previous explanations. McDonald found no
Numbers in brackets denote reference numbers listed in the
List of References.
4
evidence for the simple motion of Magnus. In addition he studied the
eyes and the vestibular organ as motion sensors.
Amar [7] made one of the most complete of the early studies of
human motor activities in 1914. This study of the relative motion of
the head, limbs and major sections of the trunk was made with a view
to study the efficiency of human motion in connection with industrial
labor.
Fischer [8] considered the mechanics of a body made up of n links
and obtained equations of motion without introducing coordinates. He
made discussions of applications of his theory to models of the human
body, but did not give applications for the equations of motion he had
obtained.
In recent years most of the analytical studies of human motion
have been associated with human beings in free fall as applied to
astronauts maneuvering in space with and without external devices.
McDonald [9] made extensive experimental studies of human motions such
as springboard diving and the "cat-drop" maneuver. McCrank and Segar
[10] considered the human body to be composed of nine connected parts.
They developed a procedure for the numerical solution of their very
complex equation of motion. Although some numerical results were pre
sented, no general conclusions were drawn.
The most significant contribution to the application of rational
mechanics to problems in the reorientation of a human being without the
help of external torques was made by Smith and Kane [11]. Specifically,
they considered a man under free fall. In this paper the authors pointed
out that the number of the unknown functions exceeded the number of the
5
equations of motions that were obtained for the system and recognized
the need for optimization considerations. In order to get the adequate
number of equations, they introduced a cost functional to be optimized
which consisted of an integral over the total time interval of some
suitable functional of the undetermined generalized coordinates.
Optimization of this functional became a problem of calculus of varia
tions, which yielded the necessary number of additional equations (the
Euler-Lagrange equations) to solve the original problem completely.
The approach of Smith and Kane suffers from one major drawback
it ignores the internal forces of the system. The internal forces due
to muscle groups at the various joints of the body segments are mostly
voluntary, have upper bounds in their magnitudes and are responsible
for the partially independent movements of the various limbs. Because
they are internal forces, it is possible to eliminate them completely
in any equations of motion. These can be obtained, for example, if the
entire system is considered as a whole. However, such equations will
be limited in number (at most six, three from the consideration of
translation and three from the consideration of rotation) and will, in
general, be less than the number of the unknown functions. The intro
duction of a cost functional will yield via the calculus of variations
the proper number of equations. However, the maneuver obtained may be
beyond the physical ability of the individual. It is therefore essen
tial to recognize the role of all the internal voluntary forces that
come into play during a physical maneuver.
Ayoub [12] considered an optimal performance problem of the
human arm transferring a load from one point on a table to another point.
6
The motion considered was planar. Internal forces and constraints on
the stresses were considered. A two-link model was considered for the
arm and numerical solutions were obtained, using the methods of Linear
Programming, Geometric Programming, Dynamic Programming and simulation.
The performance criterion was a mathematical expression for the total
effort spent during the activity. The motion considered was simple
from the point of view of dynamics. Also, it required less than ten
points to describe the entire motion. This allowed consideration of
many physical constraints but the dynamics of the human body motion con
sidered was quite different than when large motion of the various
limbs are involved. A list of references of works on human performance
from the point of view of Industrial Engineering is given in Ayoub's
thesis.
The research in the area of free fall problems showed that
analytical solutions agreed qualitatively with observations in those
cases in which rational mechanics has been applied with care. Examples
of such problems are the cat-drop problem and the jack-knife-diver
maneuver (Smith and Kane). Attempts at analytical solutions of activ
ity problems which are not in the free fall category have not been com
pletely successful. Such problems are in athletic activities (such as
running, part of the pole vault maneuver and part of the high jump
maneuver) and in activities associated with working on earth.
Simultaneously with the study of the dynamics of motion, several
investigations were made for the determination of the inertia parameters
of human beings at various configurations. Knowledge of inertia param
eters is essential for performing any dynamic analysis. A list of the
7
research activities in this area is given in References 13, 14 and 15.
However, only Hanavan [14,15] has proposed a personalized model of a
human being. This inertia model has been used in the present investi
gation.
1.3. The Problem Statement
The present work belongs to a broader program whose objective
is to investigate the basic aspects of the applicability of rational
mechanics to the solutions of any human activity problem. The aspects
presently being considered are:
1. Construction of an appropriate presonalized model for the
individual's maneuver under consideration
2. Formulation of a well-posed mathematical problem for the
analytic description of the maneuver
3. Solution of the mathematical problem .by a suitable analytical
method
4. Comparison of the analytical solution with an actual motion
conducted in an experiment with a human subject
5. Determination of muscle activity and comparison with computed
muscle torque histories for the maneuver.
If a complete analysis based on the above steps results in the
correct motion as compared with the experiment, then the results can be
used for training purposes and design of man-machine systems, but, more
importantly, the results will establish the applicability of rational
mechanics to the solution of problems of human activity.
In the present investigation, a particlar gymnastic maneuver,
the kip-up, has been selected for analysis as outlined above. The methods
8
developed will, of course, be valid for other maneuvers but the ana-
lyti cal model will be personalized to the maneuver and individual.
Among all the common physical and athletic activities, the kip-up
maneuver was found to be particularly well suited for analysis. The
motion involves large motions, and is continuous and smooth. Also, it
is planar and needs relatively fewer generalized co-ordinates for its
complete description, since a correctly executed kip-up exhibits three
"rigid" links. At the same time it is not a trivial problem from the
point of view of our basic objective.
The physical quantity (the performance criterion) chosen for
optimization (minimization, in this case) is the total time to do the
maneuver.
1.4. The Kip-Up Maneuver
The kip-up maneuver is an exercise that a gymnast performs on
a horizontal bar. The gymnast starts from a position hanging vertically
down from the horizontal bar and rises to the top of the bar by swinging
his arms and legs in a proper sequence. During the maneuver the motion
is symmetric and does not involve bending of the elbows and the knees.
Normally, the grip on the bar is loose most of the time. For an inex
perienced person, the maneuver is not easy to perform.
1.5. Present Work
In Chapter 2, a mathematical model for the kip-up motion and
results of laboratory experiments to test the accuracy of the model
are presented. First, a mathematical model of a professional gymnast
for the kip-up motion has been constructed. The dynamic equations of
9
motion for the mathematical model were then obtained. Two sets of
equations of motion were derived: one for the purpose of verifying the
accuracy of the mathematical model from experiments and one for the
purpose of optimization of the kip-up maneuver. Next, the results of
laboratory experiments with the gymnast are presented. The gymnast was
told to perform symmetric maneuvers on the horizontal bar, including the
kip-up. The maneuvers were photographically recorded. Two of the many
records were selected, one of a simple swinging motion with relatively
small oscillations, and another of his quickest kip-up maneuver.
The angle measurements from these film records were then used in the
equations of motion to check the accuracy of the mathematical model.
An error analysis was then performed to explain the disagreement between
the experimental and the computed results.
In Chapter 3, an analytic solution of the minimum time kip-up
for the mathematical model was obtained by numerical computations.
First, the analytical problem of the determination of the kip-up in mini
mum time was stated in precise mathematical terms. This involved repre
senting the equations of motion in the state variable form, specifying
the boundary conditions on the state variables, establishing the bounds
on the control variables, and modeling the stiffness of the shoulder
and the hip joints at extreme arm and leg movements. A survey of the
necessary conditions for time optimality given by the optimal control
theory has been presented. Finally, several numerical schemes used for
solving the kip-up maneuver problem are presented and the results of the
numerical computations are discussed.
In Chapter 4, final conclusions of the present investigation
and recommendations for future work have been presented.
CHAPTER 2
EXPERIMENTATION AND CONSTRUCTION
OF THE MATHEMATICAL MODEL
2.0. Introduction
In this chapter, modeling of a human being for the kip-up
maneuver is considered. In Section 2.1, a mathematical model for a
professional gymnast is constructed, using the personalized model of
Hanavan [14]. In Section 2.2, equations of motion for the mathematical
model are derived for the purpose of using them in the analytic deter
mination of the subject's optimal kip-up motion. An equivalent system
of first-order differential equations are derived in Section 2.3 for
testing the inertia properties and structure of the model. In Sec
tion 2.4, the laboratory experiments performed are described. The
results of the experiments are discussed in Section 2.5. An analysis
of the comparisons between the observed and computed results are also
presented.
2.1. Mathematical Model of the Kip-Up
The equations of motion of a deformable body such as the human
body are usually partial differential equations. Presently, not enough
is known or measurable about the deformation of the human body under
voluntary motion to determine a partial differential equation model.
Also, such models are quite difficult to handle. However, for situations
10
11
where the deformation is small compared to the displacement of a body,
the deformable body may be considered as a rigid body in writing the
equations of motion. The rigid link assumption has been used widely
in modeling human beings. The personalized model of Hanavan [14] is
based on this assumption. The elements of the Hanavan model are shown
in Figure 1 and consist of fifteen simple homogeneous geometric solids.
This construction allows a large number of degrees of freedom for the
model and minimizes the deformation of the elements without undue com
plexity. As the number of degrees of freedom increases to the maximum,
it is likely that a mathematical model based on the Hanavan model
becomes more accurate. However, increases in the degrees of freedom
increase the model's complexity and make it difficult to analyze mathe
matically.
Observation of a correctly performed kip-up indicates that the
human performer might be modeled quite accurately as a system of only
three rigid links. The two arms form one link, the head-neck-torso
system forms another, and the two legs form the third link. The
shoulder and hip joints may be approximated as smooth hinges where
these links are joined together. Deformation of the link consisting
of the head, neck, and torso during certain periods of the maneuver
is detectable when observed via high speed filming. It was felt that
the effect of the deformation of the torso would be no more significant
than the standard deviation in the Hanavan inertia parameters because
this part of the orientation of the torso is determined by the line
joining the shoulder and the hip joint centers.
12
Figure 1
o
00
>>
Hanavan's Mathematical Model
of a Human Being.
13
The muscle forces acting at the hip and shoulder to cause link
motion have been replaced with their rigid body equivalent resultant
forces and couples. If the masses of the muscles causing the forces
are small compared to the other portions of the links, then the net
effects of the forces are the torques of the couples. The resultants
do not appear in the equations of motion.
The three-link kip-up model is shown in Figure 2. It is con
structed with elements of the Hanavan model. Twenty-five anthropometric
dimensions of the gymnast were taken and used in a computer program for
calculating first the inertia properties of the Hanavan model and then
those of the kip-up model. The determination of the inertia properties
of the kip-up model from those of the Hanavan model is presented in
Appendix A.
2.2. The Equations of Motion
The mathematical model for the kip-up motion is a three-link
system executing plane motion under gravity. The active forces in the
system are the pull of gravity acting on each of the elements of
the system and the two muscle torques. The muscle forces at the shoulder
acting at the joint between links 1 and 2 are replaced by the torque u^.
Likewise, u^ is the torque for the hip joint between links 2 and 3.
The system is suspended from a hinge at the upper end of link 1, repre
senting the fists gripping the horizontal bar which is free to rotate
on its spherical bearings. All joint hinges are assumed to be friction
less. The general three-link system for which the equations of motion
will be derived is shown with nomenclature in Figure 3. Equations of
motion of the system were obtained via Lagrange's equations as follows:
0 Center of Gravity of an Element
Figure 2. Mathematical Model for Kip-Up.
15
3
CGI
a
3
*1
I
= length of element 2 = distance between the hinges A and B
- distance between the center of gravity of element 1 and the hinge 0
= distance between the center of gravity of element 2 and the hinge A
= distance between the center of gravity of element 3 and the hinge B
, CG2, CG3 locations of centers of gravity of elements 1, 2, and 3,
respectively
1
= angle between O-CGl and OA
= angle between A-CG2 and AB
, m m = mass of elements 1, 2, and 3, respectively
Z o
I^, I = moments of inertia of elements 1, 2, and 3, respectively,
about axes perpendicular to the xz plane through their respective
centers of gravity
I =
g
moment of inertia of the horizontal bar at the hinge 0 about its
longitudinal axis
acceleration due to gravity
Figure 3. The Three-Link System
16
If we define
A = §^rl + h + m2r2 +I2+m2*l + I3+ m3r3 + m34 + m34 + Ir)
B = m2Vl
C = m3X1X2
D = m3Vl
E = m3r3^2
F = ^(m2r2+I2 + m3^2+ln3r3+I3)
J mgr3 + 13
M = m^g
N = (m2 + m3)
V = m3 V
W = ra2r2g
R = m3r3g
(2.2.1)
with Equations (2.2.1), we can express the Lagrangian of the system as:
2 "2
L = cp [A + B cos (6+3) + C cos 8 + D cos (6+ijf) + E cos i|f] + 9 [F+ E cos \jf]
1 '2
+ 6cp[2F+B cos (6+3) + C cos 6+D cos (6+i|0 + 2E cos i|f] + ^ t J
+ 8 i| [J + E cos ijr] + cp \[r [J + D cos (9+iJr) + E cos i|r] + M cos (cp+a)
+ N cos cp + V cos (cp+8) + W cos (cp+9+3) + R cos (cp+6+i)i). (2.2.2)
For the Hanavan model of
If we define and u2 as positive when they tend to increase
the angles 9 and i)r respectively, we can write the equations of motion
as:
17
d Bl Bl
_d_ BL
dt Be
3l
^0
= u.
d Bl 3l
dt "5F ~ U2
3^
(2.2.3)
(2.2.4)
(2.2.5)
Let us define
= 2(A + (B+C) cos 0 + D cos (0+ \|f) + E cos i|r)
a = 2F + (B+C) cos 9 + D cos (9+i|f) + 2E cos ijt
z
a^ = J + D cos (9+t) + E cos i)i
bl a2
b = 2(F + E cos f)
bg = J + E cos i|r
c3 = J >
d^ = cp9[2(B+C) sin 0 + 2D sin (0+\|/)] + 9\|i[2D sin (9+i|0 + 2E sin \|i]
+ 0i|r[2D sin (9+i|r) + 2E sin i|r] + 02[(B+C) sin 0 + D sin (0+\jf)]
2
+ i| [D sin (9+\J) + E sin \|r] (M+N) sin 9 (V+W) sin (9+0)
- R sin (cp+04-l)
*2
d = u + i|f[29 + 20 + i|r] E sin i|r 9 [(B+C) sin 0 + D sin (0+t|r)]
Z X
- (V+W) sin (9+9) R sin (9+0+ijf)
2 2
d = u -9 [D sin (+ \¡r) + E sin i|r] 9 E sin i|r -209 E sin ljf
o z
- R sin (9+0+i|i) .
(2.2.6)
18
With Equations (2.2.6) we obtain from Equations (2.2.2) (2.2.5) the
equations of motion
a P + = d1
V + 29 + c2* = d2
a3'P + b39 + C3* = d3 '
(2.2.7)
(2.2.8)
(2.2.9)
It will be helpful to express the equations of motion in
normal form in formulating optimal control problems for the system.
For that purpose we define the state variables as
x, =9 x =9 x_=e x. = e X = i|r x- = i|r
A =
(2.2.10)
(2.2.11)
A1 =
(2.2.12)
A2 =
(2.2.13)
(2.2.14)
19
Using these definitions, the equations of motion (2.2.7) (2.2.9) can
then be written in the normal form as
To write down the equations of motion in more convenient forms, we shall
further define the following quantities:
(2.2.15)
(2.2.16)
(2.2.17)
(2.2.18)
(2.2.19)
20
A -|T
r 3 l
A(X) = |_X2, x4, X6, -J
B(X) =
(_blC3+b3Cl)
("alC3_a3Cl)
(_aib3+a3bl)
(blC2 _b2Cl)
(_aiC2+a2Cl)
(aib2a2bl)
(2.2.20)
J
("blC3+b3Cl)
("alb3+ 3bl>
(2.2.21)
(blC2"Vl)
(_aiC2+ a2Cl}
(aib2 a2bl}
(2.22)
Using the definitions(2.2.16) (2.2.22), the equations of motion
(2.2.15) may now be expressed by any one of the following equations:
?1
, c.
.B j i.S ,-5rn. S_,-n,: V
v-
Up o this point, this analysis has been very gen
era": and three-dimensional witn "-id referente to a plate-
like two-dimensional problem. Equation >3} witn the aid
of (?) car, be applied to any three-oixensional elasticity
22
with 9 expressed in terms of p using Equation (2.3.1). Hamilton's
canonic equations are then given by
cp = -g (p,cp,t) (2.3.3)
and
P = (P,9,t). (2.3.4)
From Equations (2.3.2) (2.3.4) we obtain
.
^ p-0[2F+(B+C)cos 8+D cos (6+1|0+2E cos ifr] i|r[ J+D cos (8+i|f) + E cos i|f]
2[A + (B+C) cos 8 + D cos (6+i|) + E cos \|/]
(2. 3.5)
p = [ (M+N) sin cp + (V+W) sin (9+0) + R sin (9+0+i|r)]. (2.3.6)
If we wish to introduce the effects of the friction at the hinge 0,
the equation for p becomes
oh
p = - F
'5p
(2.3.7 )
where F^ is the generalized friction torque at the hinge 0.
The Integration Scheme
The integration scheme integrates Equations (2.3.5) and (2.3.7).
The initial condition of 9 is obtained from the measured values of 9.
The initial value of p is obtained from the definition of p given in
Equation (2.3.1). To compute the initial value of p, 9 has been com
puted for this starting point only. The 6 and i|r data are differenti-
ated numerically to generate 0 and \Jr values at every step. and
are computed at every step by using these values. The value of F^ is
not known a priori and it requires a separate experiment for its
determination. Integrations were done with F^ = 0 and F^= C sgn 9
23
with C determined experimentally. The difference between the two
cases was found to be insignificant even with C=1 ft-lb which was
well above the possible friction torque at the bearings.
2.4. Experimental Procedure
The test rig consisted of a horizontal bar and two motion
picture cameras. The horizontal bar was made of a short solid round
steel bar 1-3/16 inches in diameter and 58 inches long supported by two
very rigid vertical columns through a pair of self-aligning spherical
bearings. The bearings allowed free rotation of the bar with the arm
of the subject.
One movie camera was placed with its line of sight aligning
with the horizontal bar and about 30 feet away from the bar as shown
in Figure 4. Alignment of the camera's line of sight with the horizontal
bar would give the correct vertical projection of the two arms on the
film for determining the angle cp directly. The second camera was placed
in front of the horizontal bar at the same elevation as the bar. The
film taken in this camera showed whether a particular motion was sym
metric or not about the vertical plane of motion and also, the angle
between the two arms, which is required for computing the moment of
inertia of the arms. The film speed was determined from the flashes of
a strobe light regulated by a square wave generator.
Experiments
The experiments were done during two separate periods with the
same subject, a professional gymnast. An average of 15 experiments
of the subject's performance on the bar were recorded on film each day.
Figure 4. Sketch of Kip-Up Apparatus Configuration.
25
The subject was told to avoid bending his arms and legs and to maintain
symmetric motion. In the early days of experimentation he was told to
just swing on the bar by moving his stiff arms and legs relative to the
torso. These experiments were done with the idea of obtaining small
angle data for verifying the inertia properties of the Hanavan model.
The filming of the camera was done at speeds of 32 frames per second and
62.1 frames per second (as determined from stroboscopic measurements).
In later experiments the subject was told to perform the maneuver with
some specific objectives. He was told to perform what he thought would
be the kip-up (1) in minimum time, (2) with minimum expenditure of
energy, and (3) putting "least effort." Each of these maneuvers was
repeated several times. Between any two subsequent maneuvers, the
subject was given adequate rest periods to avoid fatigue. This exper
iment was conducted with the idea of making optimization studies as well
as testing the model.
Four white tapes were stuck to the subject, on the sides of
his upper and lower arms, sides of his torso, and on the sides of his
legs. These tapes were aligned between joint centers as suggested by
the Hanavan inertia model.
Processing the Data
The film speed was measured with the aid of a stroboscope by
running the camera with a developed film with the same speed setting.
After removing the lens the shutter was exposed to the stroboscope
flash. By arresting the shutter in the stroboscope light, the shutter
speed was obtained.
26
The films were run on an L-W Photo Optical Data Analyzer.
For the purpose of testing the inertia properties, two maneuvers were
selected, one from each day's filmings. The films were projected
plumb line perpendicularly on paper fixed to a vertical wall. As each
frame was projected, the white tapes fixed on the subject, now clearly
visible in the image, were marked out on the paper of the pad by means
of a pencil and a straight edge. In this way each frame was "trans
ferred" on separate sheets of paper. The angles 9 and i|r were then
measured from these traces. The angle 9 was measured with respect to
a vertical reference. The vertical reference was obtained from a sharp
window wall in the background.
Of the two sets of data processed, one was smoothed before
using it in the integration scheme. This was the one filmed at the
higher speed for the fast kip-up motion. A plot of the raw data and
those after preliminary smoothing for this set are shown in Figure 5.
2.5 Results and Discussion
The results of the integration of the equations of motion for
two of the data sets analyzed are shown in Figures 6 and 7. Figure 6
shows data for swinging motion, while Figure 7 is for the kip-up.
Also, in Figure 7 is given the computed results corresponding to
reinitiating the integration program with the measured data. Differ
ent starting points of integration were selected to eliminate the
errors that were generated before these points.
In Figure 6, curve 2 of the computed values for the unsmoothed
data agrees well with curve 1 of the measured values for a little more
Angles (Degree)
for the Kip-Up Motion. (Film Speed 62.1 frames/second.)
to
T h c* res o 1
t s o h a i
n v i
so fur a
re applicable
to
p: a
tes cons;stir.
of any
r. jiu2-
eo* la.
yens having ci
f f e r -
ent
properties an
; geomec
O
Cut -ro
n¡ here on, the
e n a 1 y
s i s
would depend
n t h s t
f o o
O pinte
c o n s i d l r e d v i
z. s 4 r.
gl e
layer cr mult
i p 1 a y e r
plates, thin
or moderately
thick
pi a
tes etc. and
also on
the
type Of
dispi a cement
rune-
t i on s V ?,rd V ,. Habi p [41 has demonstrated the apolica-
i j -
Dili rv of these results to a single layer plate.
In this dissertation, a thres-1ayered elate, popu
larly known as a sandwich plate will be considered.
San dvr i ch oate
A sandwich plate consisting of three layers is
shewn in figure (2). Tne face layers are much thinner
than the core. All layers are uniformly thick throughout-
The two faces are of the same thickness tt, the core thick
ness being t .
Tne present analysis is capable of treating aixed
boundary value problems. On the upper surface, the line
AB separates the two regions and over which
displacements and stresses respectively are prescribed.
On the lower surface, the line CD separates the correspond
ing regions. The line CD is located exactly below the line
AD. Also, at any two points locatec on the upoer and lower
surfaces of the plate and havine tne same >:, and x.
cp (Radian)
to
CD
(Radian)
(
Curve 1 Measured and Smoothed Values of cp
Curves 2,3,4- Computed Values of cp, Started
at Different Points
2.0
w
o
31
2.6. Sources of Errors
The following are considered to be the sources of errors
responsible for the disagreement between the computed and measured
values.
2.6.1. Imperfections in the Model
The human being for the motion studies was modeled as a system
of rigid bodies. The response of the system to be compared to that of
the rigid body model was a single generalized coordinate of the system.
The errors in modeling can be lumped into the overlapping categories of
(1) definition of the generalized coordinates of the rigid elements of
the system, (2) deformations of link lines during motions from prior
joint center measurements, and (3) significant variations during motion
in the inertia properties of the torso with respect to the fixed coor
dinate system. In several experiments the torso deformed with signif
icant movement of the shoulder joint centers. In these cases the con
stant inertia properties model is obviously incorrect. These variations
not only cause inaccuracies in the inertia parameters but also result in
errors in the link lines or additional errors in the mass center of
element 2. These errors are reflected in errors in the angles 0 and i|r,
which in turn cause a time varying "phase shift" in the computed angle cp-
2.6.2. Errors in Filming and Processing the Data
The primary source of errors in filming was considered to be
caused by inaccuracies in knowing precisely where the link lines were
in relation to the film plane. As mentioned in Section 2.4, care was
32
taken to minimize this error. That is, the cameras were aligned with
the horizontal bar so that the link line of the arms projected in the
plumb line plane of the film was very nearly the correct model refer
ence line. Since the link lines of the torso and legs were nearly
plumb line vertical, no corrections of the image data were required
for these links. Filming of static thin rods which were connected to
the bar at known angles to each other and the vertical plane produces
overall measurement errors of about 1 degree standard deviation.
Errors up to 3 degrees can be expected in data from films of the motion
experiments. These errors can cause the rates of the angles to be in
error by more than 30 percent. This is the primary cause of the error
in the amplitude of the angle cp computed from the dynamical equations.
(Second derivatives of the data can be in error by more than 100 percent.
This ruled out the use of equations such as Euler's or Lagrange's.)
2.6.3. The Integration Scheme
The integration scheme uses at any step the measured values of
0 and ijf and the computed values of 9 and i|f stored previously. Once a
difference between the measured (actual) and the calculated values of cp
has developed at a time t due to any of the sources of errors discussed
above, the system configuration determined by the calculated value of cp
and the measured (actual) values of 0 and ljr at the time t will be dif
ferent from that of the actual system at that time. This will cause
the model to have a different response after time t than that of the
actual system which has a different relative configuration. This
in turn will cause a further deviation between the actual and the
calculated motion that follows this instant of time. To reduce this
33
effect of propagation of error via the system equations, several
reinitializations of cp and p were done by restarting the integrations
at different points.
Various order differentiation and integraion schemes were
tested with insignificant differences in the results for smoothed
data. These were done with data from filming at F=62.1 frames per
second with a maximum integration step size of 2/F second. As men
tioned previously, the main cause of amplitude error was the errors in
the derivatives of the raw, unsmoothed data.
The results obtained from the experiments show that the model
for the kip-up motion constructed from the Hanavan model was reason
ably good considering its kinematical simplicity. In spite of imper
fections in the model, its dynamic behavior was quite similar to the
actual motion, so that this model could provide reasonable estimates
of optimal human performance via the theory of optimal processes and
numerical solution methods.
However, the results obtained from the experiments indicate
that the application of rational mechanics to the analysis and design
of man-machine systems could prove inadequate unless the model and the
data gathering techniques can be improved. This is especially true in
the design of high accuracy or low tolerance systems.
CHAPTER 3
ANALYTIC DETERMINATION OF
THE MINIMUM-TIME KIP-UP STRATEGY
3.0. Introduction
In this chapter the determination of an analytic solution of the
kip-up maneuver is presented. The problem of analytical determination
of the kip-up strategy in minimum time has been cast as a problem of
optimal control of dynamical systems. Before the techniques of the
optimal control theory may be applied to the problem, it is necessary
to state the physical problem in the language of mathematics and to
introduce the physical constraints that must also be considered for the
solution. Thus, the first four sections of this chapter have been
devoted to the formulation of the mathematical problem. In Section 3.5
a survey of the necessary conditions for optimality obtained from the
optimal control theory is presented. Since the problem under consider
ation cannot be solved in closed form, numerical methods were used to
obtain the solution. In Sections 3.6 3.9, the choice of the numerical
methods, their derivations and the results of the numerical computations
are discussed. In Section 3.10, results of the numerical computations
are compared with the actual motion.
34
35
3.1. ^Mathematical Formulation of the Kip-Up Problem
The problem is to determine the minimum time strategy for the
man model to kip-up without violating control constraints. These con
straints represent the maximum torques the man's muscles can exert for
any given configuration. Formulated mathematically, we have the follow
ing:
For the system equations
X = f(X,u) = A(X) + B(X)u (3.1.1)
and the boundary conditions
X(0) = X
o
(given)
(3.1.2)
where
and
given by
§(X(tf))= X(tf) Xf = 0 Xf = given
t = final time, to be determined
X(t) is the time-dependent state vector,
Xx(t) =qj(t), Xg(t) =cp(t), X3(t) = 9(t)
x4(t) = 0(t), x5(t) = \Kt), xg(t) I|r(t)
(3.1.3)
(3.1.4)
(3.1.5)
find a control
u(t) = [u1(t),u2(t)]
(3.1.6)
such that simultaneously Equations (3.1.1) (3.1.3) are satisfied,
t^ is minimized, and for all values of t, 0 ^ t £ t the inequalities
S*(X) ux(t) sJ(X)
s\a) S u2(t) S S2(X)
(3.1.7)
36
are satisfied. S1(X) are given functions of X and represent the bounds
on the control u.(t). The functions f(X,u), A(X) and B(X) were pre-
J
viously given in Section 2.2. § is the error in meeting the terminal
values of the state variables.
3.2. Bounds on the Controls
The control variable u^ is the muscle torque exerted at the
shoulder joint and u is that exerted at the hip. For the individual
being modeled, the functions u and u will have upper and lower limits
which are functions of the state X.
Samras [16] experimentally determined the maximum muscle torques
at the shoulder and hip joints for various limb angles at the joints.
This was done for the same subject modeled in the present study. These
measurements were made under static conditions and the mximums in
either flexion or extension were measured for the shoulder torque for
various values of 6 and the hip torque for various values of i|i. The ex
perimental bounds on the shoulder torque were then fitted by polynomials
in 9. The hip torque bounds were expressed in polynomials in \jf.
Even though each of these bounds might be expected to depend to
some degree on all four state variables X^, X^, X^, and X^, the bounds
on the shoulder torque u depend primarily on X and the bounds on the
1 O
hip torque u depend primarily on X The measurements of Samras do not
z o
include the rate dependence X and X Although the rate effect appears
4 6
to be measurable, it is a second-order effect and quite difficult to
obtain. The control limit functions are given in Figure 8. These func
tions are correct only for a certain range of values of the angles
Control Limit (ft-lb)
Figure 8. Unmodified Control Limit Functions
(Samras [16]).
38
X and X The values of S. and S0 can never be positive and those of
3 5 12
O O
and can never be negative. Whenever these sign conditions are
violated by extreme values of the states, S'? is set equal to zero. Also,
2
from extrapolated measurement data an upper limit has been set for S at
160.0 ft-lb and a lower limit has been set for at -100.0 ft-lb.
3.3. Torsional Springs in the Shoulder and Hip Joints
Our dynamical model and the control limit functions of the
shoulder and the hip do not account for the stiffness of the shoulder
and the hip joints at the extremities of shoulder and leg movements.
It has been observed that the shoulder joints produce a resistance to
raising the arm beyond an angle of 0 30. The hip joints resist move
ment for i|f > 120, or for i|r < -35. The effects of these "stops are
important and must be included in the model, since the film data showed
that these limits were reached. There are no data available for the
stiffness of these joint stops. It was observed that, although the
joints were not rigid, they were quite stiff. It was therefore decided
to use stiff torsional spring models at the model's shoulder and the
hip joints. These would be active when the stop angles were exceeded.
For the shoulder the spring is active for 9 ^ 0.5 radian. For the hip
joint the spring is active for ijj S -0.6 radian and i|r 2: 2.1 radians.
The springs have equal stiffnesses. One generates a 100 ft-lb torque
at the shoulder for a deflection of 0.1 radian. This corresponds to a
joint stop torque of the order of the maximum voluntary torque avail
able at the shoulder for the deflection of 0.1 radian. This gives a spring
constant of Kg = 1000 ft-lb/rad. The spring forces at the shoulders
39
would therefore be equal to -K (9-0.5) for 9 ^ 0.5 radian and those at
s
the hip joints would be -K (i|r-2.1) for 4 ^ 2.1 radians and -K (4+0.6)
s s
for 4 ^ -0.6 radian. These torques at the shoulder and the hip joints
were added to the voluntary control torques u^ and u^ when the stops
were activated.
3.4. Boundary Conditions
The boundary conditions for the kip-up maneuver were chosen
from the experimental data of Section 2.5. The initial values selected
correspond to motion which has already begun. This is beyond the
initial unsymmetrical motion which occurs on beginning the first swing.
This motion is difficult to model and is not important in this basic
research. The final values of the state variables represent the model
atop the horizontal bar still moving upward and just before body con
tact with the bar. (Once the torso contacts the bar, the model is no
longer valid.) The actual motion in the experiment terminated shortly
after this point when the gymnast used the impact of the horizontal bar
with his body to stop himself.
The initial and final values of the state variables for the
optimization problem are listed in Table 1.
40
TABLE 1
BOUNDARY CONDITIONS FOR A MINIMUM TIME KIP-UP MOTION
State Variables
Initial Value
Final
Value
V
X1 = 0.340
o
1
Xf =
-2. 84
2
2
X = -2.30
o
Xf -
-7.05
e
X3 = 0.305
x3 =
2.88
O
f
#
4
4
e
X -0.660
X-c =
0.163
o
f
X5 = -0.087
o
0.436
6
6
X = -1.20
o
Xf =
0.108
3.5. The Necessary Conditions for Time Optimal Control
In this chapter, we look into the necessary conditions for the
minimum time problem formulated in the previous chapter. The necessary
conditions for optimality of motion for the case when the constraints
on the control are not a function of the states are given in Reference
17. For the case where control constraints depend on the states, the nec
essary condition requires a modification in the adjoint equations.
These are obtained through a calculus of variations approach [18].
This approach is used in the following developments.
Writing the state equations of our system as
X = f(X,u) = A(X) + B(X)u
(3.5.1)
41
we can construct the cost function as
dt
(3.5.2)
where t
f
is free.
The Hamiltonian is then given by
H(X,u,X) = 1 + \Tf
= 1 + \TA(X) + \TB(X)u
= 1 + XTA(X) + XTb1(x)u1 + \TB2(X)U2
where X(t) is the time-dependent six-dimensional column vector of
adjoint variables. A, B, B1 B and f are the quantities as defined
by Equations (2.2.15) (2.2.22) in Section 2.2.
The minimum-time control policy u(t) will be given by the one
that minimizes the Hamiltonian (3.5.3), provided no singular arcs are
present. We note that in this case the Hamiltonian is a linear func
tion of the control u and therefore the minimum with respect to u occurs
only in the upper and the lower bounds of u if there is no singular
solution. Thus, we have, recalling the definitions of S^ in Section 3.2,
T
(1) If X B^(X) >0, u^ = the minimum allowable value of u^
= S*(X) (3.5.4a)
(3.5.3a)
(3.5.3b)
(3.5.3c)
(2) If \TB (X) < 0, U;L
(3) If XTB2(X) > 0, u2
the maximum allowable value of u^.
SJ(X)
(3.5,4b)
the minimum allowable value of ur
S2(
(3.5.5a)
(4) If \TB2(X) < 0, u2
the maximum allowable value of u
(3.5.5b)
42
(5) If \TB (X) = 0
T "
or X B (X) = 0 J
, u^ u^ = possible singular control.
(3.5.6)
u and u will be determined by investigating whether or not there is
-i- C*
a singular solution with respect to these variables.
The adjoint equations will be different for the portions of the
o
trajectories for u corresponding to constrained and unconstrained arcs.
The adjoint equations are, in general, given by
\T = -H, (3.5.7)
X
This yields
(a) When neither u nor u lie on a constraint
\T = -H, = -\TA, \TB u XTB u (3.5.8)
'X -X l,x 1 "2,x 2
(b) When any one or both u and u denoted by u. (i = 1 or 2), lie on
1 Z X
a constraint denoted by (j = 1 or 2) the right side of Equation (3.5.8)
of the adjoint variables has the additional term
- \tb. SJ
- -i i
X
and the equation can be written as
T_ j
X = X A, X B, U x Bn u Z 5. x B. S.
- -'X -1, 1 -2,x 2 i=1 -1 i,
(3.5.9)
where
6=0 if X B. = 0
l -l
6. = 1 if \ B. 4 0.
i -l
The boundary conditions on the state and the adjoint variables
are
43
X(0) given = X \(0) = free
- -o -
X(t ) = given = Xf \(tf) = free
and
H(tf) = (1 + \Tf)t ^ 0 .
(3.5.10)
(3.5.11)
The state and adjoint equations together with the control laws
and the boundary conditions written above form a two-point boundary
value problem (TPBVP) in the state and adjoint variables. If these
o
equations can be solved, the optimal control, u will be immediately
obtained from Equations (3.5.4a) (3.5.6).
Investigation of Singular Solutions
T o
It has been noted that if ¡V = 0, u^ cannot be determined
from the requirement that the Hamiltonian is to be minimized with
o T
respect to u The same is true for u when\ B = 0.
Since the treatment for u^ is the'same as for u^, we shall
investigate a singular control for only u^.
T
If the quantity \ B^ = 0 only for a single instant of time,
then the situation is not of much concern because the duration of the
interval is not finite and we can simply choose u = u^(t ) or u^(t+)
or 0, where u^(t ) = control at the instant preceding t, u^(t4) is the
instant exactly after t. The situation needs special attention when
T
\ B^ = 0 for a finite interval of time.
If t S t £ t is an interval for which u is singular, it is
1 Z a
clear that, for our system
^ t t^
0
for
(3.5.12)
44
and therefore,
d T .
Ht 5^
= 0
for
(3.5.13)
or
x\ + \TB = 0 for ti t t2
(3.5.14)
or, for the interval t ^ t ^ t the following results must hold:
X cj
Case 1
Only u is singular. u is nonsingular. Since u is not
1 ^ 2
singular, u is on a constraint boundary and is given by
2
u = S
2 2
j = 1 for the lower constraint
j = 2 for the upper constraint.
The adjoint equations are given by Equation (3.5.9)
T T T T i T i
l = -X A X B1 u XJ, S, X B S
X 1 } v 1 ~ 2 f 2 2 a y
(3.5.15)
X
§i = B X
'X
= §1, + ?1 ui + 52 ^
X
(3.5.16)
(3.5.17)
From Equations (3.5.14), (3.5.15), and (3.5.16), we obtain
4-h\4, jSi + ^Si, <* + s2 s2> = 0
X XX
It is to be observed that the necessary condition (3.5.17) is not
explicit in u^.
Case 2
Both u and u are singular. The value of u is no longer
X 2a 2 2
T j
and the term X B S in the X Equation (3.5.15) drops out in this
2 X
case. Accordingly, one obtains
-XT[A,x B1 A] XT[B2 B1 Bx B21 u2 = 0 (3.5.18)
X XX
45
Proceeding from the assumption that u is singular, one would also get,
£
for this case, when both u^ and u^ are singular
T
C*X
B_ A] \T[B1 B,
^X X '
(3.5.19)
From Equations (3.5.17), (3.5.18), and (3.5.19) we can see that
T T
only if both \ B^ and \ Bare zero simultaneously, is it possible to
find a singular solution by suitable choices of u^ and u^ from the
T T t
condition (3.5.13). If only one of \ B^ and \ B^, say \ B^, is zero,
the requirement (^ B,) = 0 does not yield an equation explicit in
dt -1
d T
u as observed in Equation (3.5.17). It is thus required that (\ B.)
1 i ' l
dt
= 0, which will be explicit in u during the interval t ^ t ^ t
J_ -L £
together with the requirement that the relation (3.5.17) is satisfied
T
at t r t^. These two conditions will ensure that \ B^ = 0 in the
interval t^ ^ t S t .
It is to be noted that singular control computed by the above
procedure has not been proved to be the minimizing control. Additional
necessary conditions analogous to the convexity conditions for singular
controls have been obtained by Tait [19] and Kelley, Kopp and Moyer [20]
for scalar control and by Robbins [21] and Goh [22] for vector control.
For the general case of vector control these conditions, summarized by
Jacobson [23], may be stated as on singular subarcs:
3
- Ldt
H, =0 if q is odd
(3.5.20)
and
:> 0 .
(3.5.21)
46
d2p
In these equations, * H, (X,\) is the lowest order time derivative
dt P -
of H, in which the control u appears explicitly, and q < 2p.
For a scalar control, Equation (3.5.20) is satisfied iden
tically.
Equations (3.5.20) and (3.5.21) also do not constitute suf
ficiency conditions for minimality. A complete set of sufficiency
conditions for singular arcs has not yet been established in the
literature of optimal control theory for a general nonlinear system.
We can see that there are quite severe restrictions on the
existence of singular arcs in the human motion problem. In the
numerical methods used in the present work to determine the optimal
solution, only in the method of quasilinearization is it necessary to
express the control (its optimal value) in terms of the state and
adjoint variables, while in the gradient methods where successive
improvements are made in the control variables, this is not so. In
the attempts with the quasilinearization method, singular solutions
were not considered in the construction of the two-point boundary value
problem in the state and adjoint variables. It was decided that if
a solution to the TPBVP was obtained by quasilinearization, singular
arcs would be looked for later. The gradient methods exhibit singular
arcs automatically if there are any. The additional necessary condi
tions for singular arcs should be checked when off-constraint arcs are
exhibited by the gradient method.
47
3.6. _The Solution Methods
The optimal control problem formulated in the preceding section
cannot be solved in closed form. Numerical methods must therefore be
used to find its solution. In the optimal control theory literature
several numerical methods have been proposed for solving the differ
ential equations and the optimality conditions that arise out of optimal
control problems such as the present one. None of these methods guaran
tees that a solution will be obtained readily, while some of the methods
do not guarantee that a solution may be obtained at all. The methods
are all iterative, necessitating the use of high-speed computers for
all nontrivial problems. A nominal guessed trajectory is improved
iteratively until the improved solution satisfactorily meets all the
necessary conditions.
Depending on whether the method requires finding the first or
both first and second derivatives of the system equations with respect
to the state and control variables, these methods are called First-Order
or Second-Order methods, respectively. This is so because they, in
effect, make first-order or second-order approximations of the system
equations with respect to the state and control variables. The first-
order methods, in general, have the property that they can start from
a poor guess and make fast improvements in the beginning. They need
fewer computations in each iteration. But their performance is not
good near the optimal solution where the convergence rate becomes very
poor. The second-order methods, on the other hand, need a good ini
tial guess to be able to start but have excellent convergence prop
erties near the optimal solution. Because the second-order methods
48
need computation of the second derivative of the system equations,
they need more computing time per iteration, which may be excessive
for some problems.
Apart from the first- and second-order methods mentioned above,
there is another class of methods which tries to combine the advantages
of both of these methods while eliminating the disadvantages of both.
The Conjugate Gradient Method, Parallel Tangent Method, and the Davidon-
Fletcher-Powell Method fall into this class. These methods work very
much like the first-order method except that, in the first-order
expansion, the coefficients of the first-order term, or the gradient
term, is modified by some transformations. These transformations are
generated from the modified gradient term of the previous iteration and
the gradient term of the current iteration. This has the effect of
using the information that is obtained from a second derivative.
It is not known which of the several methods used for solving
optimal control problems is good for a given problem and one may have
to try more than one method in order to obtain the solution. In the
published literature, most of the illustrations of these methods are
simple. In these simple problems control or state variable histories
do not have wide oscillations or the system equations themselves are
not complicated. This makes it truly difficult for someone without
previous experience to decide upon the merits of these methods.
There is no preference list, and it seems certain that there cannot
be one whereby a decision can be made as to which method should be
tried first so that a solution of a given problem will be obtained
most efficiently. In this respect, deciding upon a computing
49
method for a given problem is still an art and depends largely on the
previous experience of the individual trying to solve the problem.
In the attempts to solve the minimum time problem, the method
of quasilinearization was taken up first. This choice was based on
several factors. This is the only method where the two-point boundary
value problem obtained from the necessary conditions of optimality is
solved directly, and this feature was found very attractive. As a start
ing guess, this method requires the time histories of the state and
adjoint variables. Time histories of the state variables were available
from the experiments. (It was decided that if the method was success
ful for this guess, an arbitrary and less accurate initial guess would
be tried later.) When it converges, the method has a quadratic conver
gence rate. Also, in spite of its being a well-known method for solv
ing nonlinear two-point boundary value problems since it was first
introduced by Bellman and Kalaba [24] its applications in solving
optimal control problems have been very few. There was thus an added
incentive for using this methodto determine its usefulness in solv
ing fairly complicated optimal control problems.
Sylvester and Meyer [25] proposed, with demonstrations, an
efficient scheme for solving a nonlinear TPBVP using the method of
quasilinearization. This scheme was available in the IBM SHARE program
ABS QUASI and was used by Boykin and Sierakowski [26] who reported
excellent convergence properties of the scheme for some structural
optimization problems.
With this record of success, the program QUASI was taken up for
our problem. But with our problem several difficulties were encountered
50
from the very beginning. First, the bang-bang control law obtained
from the necessary conditions had to be replaced by a suitably steep
saturation type control law. Second, a slight modification in computa
tion scheme was necessary when it was found that the method was unable
to solve a simple example problem. The example problem could be solved
with these modifications. But, in spite of all these changes and sub
sequently, many attempts to generate a guess of the adjoint variables,
the method could not be made to work for the human motion problem.
Reasons for the difficulties encountered are discussed in detail in
Section 3.7.
During the attempts with quasilinearization, it was found that
computations of the second derivatives of the system equations were
taking an exorbitant amount of time and this was the deciding factor
for the next choice of a computing method. Also, the appearance of the
control function linearly in the Hamiltonian put restrictions on the use
of most of the other second-order methods.
The next attempts were based on the first-order steepest descent
method proposed by Bryson and Denham [27,28]. The most attractive
feature of this method is that the various steps involved in it render
themselves to clear physical understanding. This method directly reduces
the cost function in a systematic way and one obtains good insight into
the basic steps in the iterative computations and can make adjustments
to improve convergence and/or stability with relative ease. These
features of the method of steepest descent may more than offset the
advantages of other methods for some complicated problems. In the
attempts with this method, three different formulations of the minimum
51
time problem were tried. In the first formulation the computations were
not pursued beyond a certain point due to computational difficulties.
The solution was obtained by the second formulation and verified by the
third formulation. These attempts are discussed in Section 3.8.
3.7. A Quasilinearization Scheme for Solving
the Minimum-Time Problem
In Section 3.5 the adjoint equations and the optimal control
laws (Equations (3.5.9), (3.5.4a)-(3.5.6)) have been derived for the
minimum time kip-up problem. The system equations and the boundary
conditions on the state and the adjoint variables are given by Equa
tions (3.5.1) and (3.5.10), respectively. From these equations we can
readily see that if the control variables u^ and u^ appearing in the
system and adjoint equations are replaced by their optimal expressions
in terms of the state and the adjoint variables, one obtains a non
linear TPBVP in the state and adjoint variables. If these equations
are solved, the optimal state and adjoint variable trajectories will
be obtained and the optimal controls can be constructed by using the
state and adjoint variables and the optimal control laws.
In the TPBVP in the state and the adjoint variables, the final
time is not a given constant and is to be determined from the implicit
relation (3.5.11). This makes the problem one with a variable end
point.
The method of quasilinearization is formulated primarily for
a fixed-end-point TPBVP. In problems with variable end points, the
adjustment of the final time is usually done by a separate scheme, not
integral with the quasilinearization scheme. Long [29] proposed a
52
scheme for converting a variable end point problem into a fixed-end
point problem with the adjustment of the final time built into the
quasilinearization process. For the present system, however, this
scheme was not practicable because the boundary condition (3.5.11)
becomes too complicated to handle in this formulation.
It was decided that with a separate algorithm for adjusting
a
the finai time, described later in this section, the nonlinear TPBVP
with free final time would be converted to a sequence of nonlinear
TPBVP's with fixed final times. Each of these fixed final time problems
would then be solved by the modified quasilinearization algorithm until
the correct final time was obtained. The derivation of the modified
quasilinearization algorithm is described below.
3.7.1. Derivation of the Modified
Quasilinearization Algorithm
The fixed final time nonlinear TPBVP to be solved falls in the
general class of problems given by
dy
^=g(i:,t). (3.7.1)
With the boundary condition
Bj Z) + Br + c = 0 tf = given (3.7.2)
y, g, and c are of dimension n, B^ and B^ are matrices of dimension
(nxn). It is being assumed that the TPBVP has been defined for the
interval 0 £ t < t for some given t > 0.
In the state and adjoint equations, if the expressions for optimal
control in terms of the state and the adjoint variables are used for the
control variables, one obtains
53
X = f(X,u(X,\)) = F(X,\) (say)
(3.7.3a)
and
X = -H^(X,u(X,\),X) = G(X,X) (say). (3.7.3b)
In the formulations of the TPBVP given by Equations (3,7.1) and
(3.7.2), it may be seen that for the kip-up system, n=12,
X
~F
0
0
I
0
y =
_X
, g =
1
IO 1
1
i
0
II
u
CQ
0
0
and
-o
The 0 and I appearing in the matrices B. and B represent 6x6 order null
J Y
and unit matrices, respectively.
Let z(t) be an initial guess vector for y(t) which satisfies
the boundary conditions (3.7.2). If g(y,t) is approximated by its
Taylor series expansion about g(z,t), keeping only the first-order term,
one obtains
(y-z).
y=5
9g
g(y,t) = g(z,t) +
Let
W =
so that,
W. .
ij
(3.7.4)
or, W.. = partial derivative of the i
ij
th
. j element of y, evaluated at y = z.
. th
element of g with respect to the
54
With the above approximation of g(y,t), Equation (3.7.2) becomes
dy
= g(x,t) + W(z,t)(y-z)
at - - -
or,
de dz
= -rf + g(z,t) + W(z,t) e
dt dt - -
where,
e = y(t) z(t) = error in the guess z(t).
Rearranging the above equation, one obtains
de dz
W(z,t)s = + g(z,t). (3.7.5)
dt - dt -
Since z(t) is chosen to satisfy the boundary conditions,
Bj, z (0) + z(tf) + c = 0.
Subtracting this equation from Equation (3.7.2), one obtains the
boundary conditions on the error e(t) as
Bj £(0) + Br e(tf) = 0. (3.7.6)
Equations (3.7.5) and (3.7.6) form a linear TPBVP in e(t) which, when
solved, will give the values of the error between the guessed solution
z(t) and the actual solution y(t) based on the linearized expressions
of the right side of Equations (3.7.1) about the guessed solution z(t).
Because of using the linearized equation instead of the full nonlinear
"'n
equations, the values of e(t) obtained by solving Equations (3.7.5) and
(3.7.6) will not be the actual error between the guess z(t) and the
solution. However, a new guess of y(t) will be obtained from e(t) by
z (t) = z(t) + T) e(t) ,
0 < 7] < 1.
(3.7.7)
55
The algorithm of Sylvester and Meyer uses T| = 1 for all the
time, which is the usual quasilinearization algorithm. It was found,
while solving a simple example problem, that without the incorporation
of a multiplier 7] in the expression (3.7.7), i.e. using
z' (t) = z(t) + e(t)
the method was unstable. The convergence property of the scheme with
the incorporation of the multiplier 7] can be understood for small values
of 7] by comparison with the step-size adjustment procedure of the usual
steepest descent algorithms. The mathematical proof for the convergence
property follows the proof of Miele and Iyer [30] and is now given.
The integral squared norm of the error in the guessed solution
z(t) can be expressed by the integral
Similarly, the error in the solution z'(t) = z(t) + 7] e(t) is given by
t
f
If 7] is sufficiently small, one can write
g(z',t) = g(z,t) + g,z(z,t) Tie
where g, = g,
- z y
= W(z,t) (from Equation (3.7.4)).
y=z
Also, for all values of 7],
/
z
= z + 7] e .
56
From these results, one obtains, for small values of 7],
j' j = 27] J {z g(z,t)}T {e We}
Since e(t) satisfies the differential equation (3.7.5), this finally
yields:
t
j' J = 27] J || z g(z,t)||2dt
= a negative quantity.
Thus, for sufficiently small values of 7], the reduction in the cost is
assured.
In the quasilinearization algorithm z^(t) takes the role of z(t)
as the new guess of y(t) and the process is continued until the error in
satisfying the differential equations is reduced to an acceptable value.
The linear TPBVP of the error Equation's (3.7.5) and (3.7.6) is
solved as follows:
The time interval t = 0 to t = t^ is divided into m small inter
vals. This results in m+ 1 values of t at which the solution will be
computed. Equation (3.7.5) can be written in a finite central differ
ence scheme as
rl
/Z.+Z. t.+t.
[-1 -1+1 1 1+1]
e. + e. z. -z.
1 -i -i+l -i+l -i i
/z. + z. t. + t.\
f -i+l -l i+l i)
hi
' 2 2 /
2 h. +^'
i
^ 2 2 /
th
where h = t t The subscript i denotes values at the i
i i+1 i
station, i=1,2,...,ra+l. Rearranging and simplifying the above expres
sion one obtains
( + i h.w. )e. + (-1 + = h.W. )e = r. (3.7.8)
2 i i -i 2 i i -l+l -l
i = 1,2,... ,m
57
where
r. = z
-l -l+l
rui+ii i+i+^
-5i-hi 5 V 2 2 ')
/z. + Z. t. + t X
(3.7.9)
(3.7.10)
and
I unit matrix of dimension n xn .
The boundary conditions, Equation (3.7.2), reduce to
B. e + B e = 0
i -1 r -m+1
(3.7.11)
Equation (3.7.8) can be cast into the following convenient recursive
expression
e. + D. e. = s.
-l l-i+l -i
(3.7.12)
where
and
D. = (I + i h. W.)"1 (-1+ i- h. W.)
l 2 i l 2 i l
- 1 -1
s. = (I + h. W. ) r.
-l 2 i i -l
(3.7.13)
By repeated substitution, equation (3.7.12) yields the following rela
tionship between e_ and e
-1 -m+1
. m
m
e. = T + (-1) ( T D.) e
-1 i -m+1
l-l
where
m
i-1
i-1
T = s + E (-1) ( TT D.)s.
i=2 j=l J 1
(3.7.14)
(3.7.15)
or, multiplying by on both sides of Equation (3.7.14),
B. e = B. T + (-l)m B. TT D. e
JL -1 l l i -m+1
1=1
m
(3.7.16)
Equations (3.6.16) and (3.7.11) can be solved simultaneously for
and e to give
-m+1
58
-m+1
where
m
C = (-l)m B ( TT D.) + B .
I i x r
1=1
(3.7.17)
(3.7.18)
With e determined, e e ,. . e_ are determined in succession by
-m+1 -m -m-1 -1
using the recursive relation (3.7.12). Then, the new iterate is given
by
zil^1(t) = zi(t) + T) e(t).
The stopping condition of the algorithm is given by the fact that
e should be small. When they are small, it may be seen from Equation
(3.7.8) that the quantities r will also be small. From Equation (3.7.9)
it is also seen that small r correspond to satisfying the central differ
ence expression of the differential equations. That is, the finite dif
ference equation error must be small. This does not mean that the dif
ferential equation error is small unless the intervals for the differences
are sufficiently small.
The IBM SHARE program ABS QUASI is a program of the procedure
outlined above without the provision of the multiplier T) in Equation
(3.7.7), and therefore is for 7] = 1. The program was modified to intro
duce and to adjust 7] to get the desired convergence. The algorithm may
be described by the following steps:
1. Set up the matrices B. and B and guess a nominal trajectory
Ju r
z(t) that satisfies the boundary conditions (3.7.2).
Set ITER = 0
(ITER for Iteration)
59
2. Do the following for i = l,2,...,m.
a. Find r and W. as defined in Equations (3.7.9) and
-i i
(3.7.10). Find the largest element of r^, searching
between the elements of each r ^, for i = l,2,...,m.
Call it E2MAX. If E2MAX < specified maximum allowable
error, print out z and stop the computations.
b. Using Equation (3.7.13), find and s^.
3. Calculate T and C according to Equations (3.7.15) and (3.7.18).
Calculate the integral norm of the error
(here the norm is defined by the sum of squares of the elements
of the vector r.). Set J1 = J2.
-l
4. Find e using Equation (3.7.17).
-mf 1
5. Decide upon a value of T], Discussions on the choice of T] will
be presented in a later section.
6. Generate and store e e _,,e e. using Equation (3.7.12).
-m -m-l -m-2 -l
Generate the new guess z^, i = 1,2,. . ,nn-l, by doing
z. = z. + T) e..
-i -i -i
Set ITER = ITER + 1.
7. Find J2 and r i = l,2,...,m and find E2MAX as in step 2.
8. If J2 > Jl (unstable), go to step 11.
9. Stop if E2MAX < a prescribed value.
10. If J2 < Jl to to step 12 to continue to the next iteration.
11. If this step has been performed more than a specified number
of times in this situation, go to step 13. If not, store the
60
value of the current Tj and J2. Recover the values of
z of the previous iteration by doing
z = z T1 e .
-i -i -1
Reduce the value of T]. Generate the new by doing
z. = z. + T] e. V
-i -i -i
Go to step 7.
12. If this step has been performed more than a specified number
of times or if step 11 has been performed at least once in
this iteration, go to step 13. If not, store the value of J2
and T]. Recover the new z^ as in step 11, this time by increas
ing 7] to increase speed of convergence. Go to step 7.
13. Find out the minimum value of J2 that have been obtained in
steps 11 and 12. If this value is greater than Jl, stop com
putations and look for the cause of the instability. If this
value is less than Jl, recall the T] corresponding to this J2 and
regenerate the z^ for this case. Go to step 2.
In this program the value of T] was selected from the consider
ation that at the station m+1 the maximum value of
12
should be less than a prescribed value (j in parentheses represents the
th
element of the vector). This procedure limits the changes in the
j
existing trajectory by limiting the magnitude of the maximum fractional
change in the terminal values of the variables not specified at the
final time.
61
Whenever an iteration was found unstable, T] was reduced by
half. When there was an improvement, a linear extrapolation formula
was used to increase the value of T] so that the norm of the error J2
would decrease to a desired value. In such an attempt, however, T)
was not allowed to increase beyond a certain multiple of its existing
value.
3.7.2. Approximations of the Optimal Control
for the Kip-Up Problem
The method of quasilinearization disregards the question of
singular solutions in the present investigation. It was found in
Section 3.5 that there were many requirements for the existence of
a singular control arq and the necessary conditions for the existence
of singular controls are quite complicated. It was decided that, before
going into those cases extrema without singular arc would be looked for
first and if the computational method was successful, singular subarcs
v/ould be looked for later.
If the singular solutions are not considered, the optimal
o o
controls u^ and u^ become bang-bang and are given by Equations (3.5.4a)
through (3.5.5b).
It is convenient to rewrite these equations at this point
in the following way:
!f -X B (X)
< 0, u = S*(X)
= 0, ux = 0
> 0, u = sJ(X)
(3.7.19)
62
if -X B (X)
<0, u = Sg(X)
=0, u = 0
(3.7.20)
>0, u =
S2(
Let these control laws be expressed by the following expressions:
u = S1 sgn (-XT§1) (3.7.21)
and
where
and
i = 1 or 2
U2 = S2 Sgn (^B2)
sgn (x) = sign of x
= S1 when (-XTB.) < 0
1 -1
= S? when (-X^B.) > 0
i -x
(3.7.22)
(3.7.23)
(3.7.24)
With these expressions for the control laws, the state equations (3.5.1)
become
X = A(X) + B S1 sgn (-X11^) + §2 S2 sgn (-XTB2) (3.7.25)
= F(X,X) (cf. Equation (3.7.3a)).
In the present quasilinearization algorithm the derivatives of
F(X,X) with respect to both X and X in constructing the matrix
(Equation (3.7.4)) are needed. One can see that computation of the
derivative of F(X,X) with respect to X will occur only as a general
ized function with no numerical value for computation in the limit,
because X appears only in the argument of a sign function. Instead
63
of proceeding to the limit the following approximation of the bang-bang
control was made. For i = 1 and 2
o *
u sa u =
i i
rsJ(X) if AA (-\TB ) < S*(X)
AA.(-\TB.) if S2(X) > AA.(-\TB.) £ S1(X)
1--1 l- l -l l -
S2(X) if AA.(-\TB.) > S2(X)
_ 1 1 1 1
(3.7.26)
*
where AA^ and AA^ are two positive constants. This function of
T 12
AA^(-\ B^) is called the saturation function (sat) when and are
unity. The change made in the optimal control is shown graphically in
o *
Figure 9 for u^ and u^, i = 1 or 2.
o *
The controls u. and u. have been plotted against the function
i i
T
-\ B^ near a switch point in Figures 9a and 9b, respectively. It can
* o
be seen that the approximation u^ differs from the optimal control u^
*
only in the portion KL. By increasing the value of AA^, u^ can be made
o
closer to u..
i
With the above approximation of the bang-bang control by a
"saturation" control represented by Equation (3.7.26), one would
*
be able to find the derivative of the control u^ (and hence, of F) with
respect to X. The derivative will be zero when the control is on a
constraint (X) and will be nonzero on the arc KL which appears near
a switching time. In order to represent the saturation control (3.7.26)
by its linearization, it is necessary that at least one of the W ,
i = l,2,...,m, which contains the information of the first derivatives,
be computed on the arc KL. Otherwise, this vital portion of the control
would go unaccounted for in the linearized equations. This may happen
(b) Approximation of Optimal Control
(Saturation)
Figure 9. Approximation of Bang-Bang Control by Saturation Control.
o
65
if the arc KL is too steep for a given selection of integration stations.
In such a case, when none of the W is computed on the switching por
tions of the control variables like the arc KL,the TPBVP, Equations
(3.7.5) and (3.7.6),cannot be solved as explained below.
First, from a physical reasoning, it can be seen that the linear
ized state equations would get decoupled from the adjoint equations if
all of the W^, i = l,2,...,m, show zero derivatives of F with respect
to \. This means that the first six equations of (3.7.5) would get
decoupled from the last six. But the boundary conditions (3.7.6) is such
that they are for the first six variables of e(t) only. Therefore, this
results in a situation where there are six first-order equations with
twelve boundary conditions, a situation which, in general, does not
have a solution.
From the point of view of computations with the present algorithm
it can be shown that the matrix C in equation (3.7.17) cannot be used to
solve for e
-m+1
With the system equations defined as (Equations (3.7.3a) and
(3.7.3b))
X = F(X,X)
X = G(X,X)
one obtains
[W], =
th
where the subscript i represents the i
interval.
66
Let P., Q., E., DD., R. and M., j
J J
1,2,3, and 4, represent
6x6 matrices. Let
[W], =
P P "
1 2
P P
3 4J.
where P = F, = [0] for all points other than those lying on the
2 L 4J
arcs such as KL in Figure 9b. Then
1
+ II V'1
1 + f hi P1 I hi P2
2 hi P3
I + 77 h. P.
2 l 4
\ Q2
^3 Q4
It may be seen that Q = [0] if P = [0].
Z Z
From Equation (3.7.13),
Di = [" + \ hi 'l]'1 [-; + 5 hi wi]
r9! q2
rR.
R7
X
E_1
i
2
1
2
R
E .
J
L 3
4J
L 3
4J
(say),
where
R1 R2
R3 R4
r-i + \ h.w.
2xi
and
R2 = 2 hi p2 =
P9 = [0].
[0] if
67
It may again be seen that
E9 = QlR9 + Q2R4 = [0-* since = Q2 = [0] when P2 = [0],
Now, consider the product D. D. D. of any three succes-
i l+l i+2
m
sive D in the product T D.. Let this multiplication be expressed as
1 i=l 1
E11 E21
E31 E41
E12 E22
E32 E42
E13 E23
E33 E43
DD1 DD2
DD3 DD4
The expression for DD2 is
DD2 = E11 E12 E23 + E21 E32 E23 + E11 E22 E43 + E21 E42 E43
If the upper right elements E ^ = E?2 = E^ = [0], it may clearly be
seen that DD2 = [0].
In the product
m
TT D. =
i=l 1
I
M3 M4
one would therefore obtain
M = [0] if P = [0] for all i, i = l,2,...,m.
With M_ = [0], the expression for C (Equation (3.7.18))
D = (-l)m B. Tr D. + B
H l r
1=1
m
becomes, after using the values of B. and B ,
At ~!C
68
0
0
, B =
I
(f
r
I
0
0
0
"i 0
i 0
m
, x ni
c = (-1)
= (-D
_L M2_
1
o
j-1
Thus, C becomes singular if F,. = 0 for all i, i = l,2,...,m.
A
For a given step size or subdivision in the t-interval (0,t^)
the singularity of C sets an upper limit on the steepness of the arc KL
of Figure 9b whichcan be used for approximating the bang-bang control.
This steepness of the arc KL depends on AA^ (or AA^) for a given X(t)
and X(t). So, in selecting AA^ and AA we have to make sure that they
are not too large. When the solution is obtained, however, the slope of
KL does not depend on the choice of AA and AA if we decide to choose
1 z
AA1 = AA = AA, say. This is because AA will be "absorbed" within X(t) ,
acting as a scale factor on X. In fact, if we define X as AA*X, the state
and adjoint equations of our system do not change, and we could therefore
select AA^ = AA^ = 1. The slope of the arc KL in the solution depends
on the value of AA*X and not on AA alone. However, when the solution
has not yet been obtained, the best value of AA need not be 1.
If the final time, t^, selected for our problem is much larger
than the minimum time, for a bang-bang optimal solution, it can be
expected that the arc KL will have a relatively small slope. If the
final time is gradually reduced, this slope will increase until it
becomes so large that the matrix C becomes singular as explained above.
At that point the solution obtained for the smallest t^ will very
69
nearly be a bang-bang control and will represent the approximate solu
tion of the minimum time problem.
The computation of the optimal final time via the quasilinear
ization method was done according to the stopping condition outlined
above. A final time would be guessed, and the TPBVP (Equations (3.7.1)
and (3.7.2)) would be solved. The final time would then be reduced by
reducing the integration step size, and the TPBVP would again be solved
The process would be continued until the matrix C becomes singular and
the TPBVP could not be solved any further.
3.7.3. A Simple Example Problem for
the Method of Quasilinearization
A simple problem was first taken up to explore the various
features of the quasilinearization algorithm constructed above. The
system was defined as
X
X.
2
1
X
u
2
(3.7.27)
(3.7.28)
The constraints on u were
-1 £ u £ 1
(3.7.29)
The cost function to be minimized is the final time t^.
70
The adjoint equations for the system are
^ = 0
L = -x.
(3.7.30)
The optimal control u is
u = sgn (-X )
(3.7.31)
With the approximation of the optimal control
u = sat (-AA\ )
Ct
(3.7.32)
the state and adjoint equations become
X. = X
x =
sat (-AA X2)
X, 0
(3.7.33)
X0 = -X.
1
The boundary conditions are given by Equations (3.7.28). The analytical
solution of the optimal control is given by
t
f
2
X
1
-1
X
2
t 1
u = 1
Xx = | t2 ) for 0 S t < 1
*\
u = -1
X^ = -1 + 2t t 2 > for l
X = 2 t
2 J
The solution is shown graphically in Figure 10.
The problem was solved numerically by solving Equations (3.7.33)
and (3.7.28) for difieren' final times t by quasilinearization.
(C) (d)
1.0
0
-1.0
(e)
Curve 1 - X
Curve 2 X
Curve 3 -*
Curve 4 - \
Curve 5 - u
71
Figure 10. Graphs of Optimal and Nearly Optimal
Solutions Obtained via Quasilinearization
for Simple Example.
72
The solutions for t = 2.25, 2.05, and 2.005 were obtained and are shown
in Figures 10a, b, and c, respectively. With t = 2.00 the matrix C
became singular in iteration 10. The theoretical solution (t^ = 2.00)
is shown in Figure lOd. The results for t =2.005 is a reasonably good
approximation of the optimal solution.
For these problems, 25 time subdivision intervals were used .
The initial guess was deliberately poor, as shown in Figure lOe. The
solutions were obtained in 9 to 11 iterations from this guess in all
these cases.
The program originally written, according to Sylvester and
Meyer [25] did not have the provision for the amount of adjustment T)
in the iterations. Their method was found to be unstable for some of
the initial guesses.
3.7.4. The Results With the Kip-Up Problem
The kip-up problem was taken up after the method of quasilinear
ization was found successful in the case of the example problem. In the
kip-up problem, however, many difficulties were faced from the very
beginning and the problem could not finally be solved by this method.
A major difference between the human motion problem and the example
problem or the problem solved by Boykin and Sierakowski [26] is very
prominent. In the latter problems, the control variables switched only
once from one boundary to another in the entire trajectory, whereas,
in the human motion problem, there were many such switchings. This
made the human motion problem less amenable to iterative methods.
73
The program for the human motion problem was extremely lengthy,
taking more than eleven hundred statements and requiring the use of the
large core of the computer (IBM, 360-65). The subroutine (DIFEQ) which
generated the right side of the state and the adjoint equations and
their derivatives, turned out to be quite lengthy and required an
exorbitant amount of computing time. The quasilinearization program
called this subroutine at every station of the total interval and at
every iteration. As a result, the program required a tremendous amount
of computing timeabout 40 seconds per iteration. All computations
were done in double precision.
Several sources of difficulty were detected in the unsuccessful
attempts to solve the human motion problem by quasilinearization. The
central issue was the tremendous amount of computations with accuracy
required.
The total time interval was divided into 100 equal parts (seg
ments). As the first (and the only one) guess of the final time, these
segments were chosen to be 0.0216 second each. This made the final
time (2.16 seconds) equal to the time taken by the gymnast to do the
actual maneuver. With 100 segments, the program required the large
core of the computer. Even though a larger number of smaller segments
would have been preferred, this could not be done because it was
intended not to go beyond the large core. It was later found that
numerical integration of the nonlinear state equations needed time
steps at least as small as one-sixth of what was taken for quasilinear
ization. It was initially hoped that, since the quasilinearization
algorithm solves the linearized equations instead of the nonlinear
74
equations, the central difference solution would be stable for larger
integration step sizes. This was not the case.
It was proved theoretically that with perfect precision the
method was convergent for all initial guesses. However, it is well
known for the unmodified method of quasilinearization (7] = 1) that the
method is convergent only for certain initial guesses. It may therefore
be fair to expect that even with the modified method, the rate of con
vergence should depend on the initial guess and for some initial guesses,
this convergence may be extremely poor. For this reason, several initial
guesses were tried. The guess of the state variables was taken from the
experimental data which were shown to agree .veil with the computed motion
in Chapter 2. Different initial guesses for the adjoint variables were
tried. The first attempts simply used constant values for all the
adjoint variables. Convergence from this guess was nonexistent. Next,
attempts were made to generate the adjoint variables from forward inte
gration of the adjoint equation with the guessed values of the state
variables and a guessed value of adjoint variables at time t = 0. In
these cases the integrations were unstable with large numbers generated.
The method was not pursued further.
Lastly, the method suggested by Miele, Iyer and Well [31] was
tried for generating the initial guess of the adjoint variables. In
this method, an auxiliary optimization problem is constructed from
the original problem. It tries to make an optimal choice of the adjoint
variables such that the cumulative norm of the error in satisfying the
adjoint equations for a given state variable trajectory is minimized.
This is performed as follows.
75
Suppose u (X,X) is the optimal control. The state and adjoint
equations may then be written as
X = f (X,u (X, .)) = F(X,X)
X = -f^X,u*(X,X))X = G(X,\)
(3.7.34)
(3.7.35)
Suppose we have a guess of the state and the adjoint variables
given by X(t) and X(t). Since X(t) and X(t) do not satisfy Equations
(3.7.34) and (3.7.35), let
e = X + f,Tx>u*(X,X))X
so that e is the error in the adjoint equations (3.7.35). The above
equation can be rewritten as
T
X = f f,x(X,uX(X,X))X (3.7.36)
Now, consider the optimal control problem, where
a. X(t) is a given function of time t
b. X(t) is the state variable
c. e is the control variable
d. The system equation is given by Equation (3.7.36)
e. The cost function to be minimized is
1 r> T
j = I e e dt
2 2 J -
o
f. t^ is the final time (fixed) of the original problem
g. Boundary conditions are
X(0) and X(t ) are free.
76
For the optimal control problem posed above, we can construct
T
the Hamiltonian, using the Lagrange multipliers (6-dimensional
vector):
T
h2 = \ It + £*(£ ~ £,x X) (3.7.37)
The necessary conditions for optimization are
T
h9 ox
2,e "
and
e* = (H >
2 X
or, after performing the differentiations
H, = ^ + ?- 0T
2,f -* -
and
e = f, e
-* X -*
(3.7.38)
(3.7.39)
The boundary conditions are
f*(0) = f*(V = 2 (3.7.40)
because X(t) is free at t=0 and at t = tf.
Using Equation (3.7.38) in (3.7.39) and (3.7.40), one obtains
e = f, e (3.7.41)
~ ~ A.
and
e(0) = e(t ) 0 (3.7.42)
Clearly, Equations (3.7.36), (3.7.41) and (3.7.42) form a linear TPBVP
in X and e. This problem should be easier to solve, in theory, than the
original TPBVP of the state and adjoint variables and should give an
optimal choice for the multiples of X-
77
In several attempts this linear TPBVP could not be solved for
the guessed state trajectory by either the quasilinearization program or
by the IBM scientific subroutine package program DLBVP.
The prime reason for the failure of the method of quasilineariza
tion was finally found to be the numerical inaccuracies in the compu
tations which were dominant in spite of using the double precision
arithmetic. The problem was perpetuated and amplified by the large
numbers in the right side of the state and adjoint equations and in the
derivatives of these equations. The matrices to be inverted at the
various stages, one at every step of integration and another at the end
of the integration, were ill-conditioned for inversion. When a dif
ferent subroutine for inversion than what came with the QUASI program was
m
used, different numbers resulted. The entries of the matrices T and Tf D.
i=l i
were very large and resulted in large numbers for some entries of C.
This made the matrix C ill conditioned for inversion. Any error in
inverting the matrix C would be amplified in the values of e
This amplification was due to the multiplication of the inverse of C
by T, a matrix used in generating s If e was in error, all other
- -mfl -m+1
would be in error because they were generated from
If the inversion of C was accurate, then the first six entires
of computed from Equations (3.7.16,17,18) should be almost zero. But,
3
instead, large values of the order of 10 were obtained! This obviously
indicated that and hence all the were being computed inaccurately.
The quasilinearization program failed to solve the human motion
problem due to the above three reasons, and primarily due to the last
one, the numerical inaccuracies. This was felt to be rather difficult
78
to overcome since it was intimately related to the method used to
solve the linear TPBVP and the structure of the original TPBVP. So
far as the method was concerned, the key point was that the matrix C was
becoming ill conditioned. This occurred because the recursive rela
tionship between and ^ has been used to generate a relationship
m
between e and s which resulted in the product T D. with large
I IIH-1 jL=l ^
ra
entries. A look at the expression for C in terms of T D., B.,
i=l 1 i
and B (Equation (3.7.8)) would make it clear that with such a value
r
m
of T D C would automatically be ill conditioned.
i=l 1
The other standard methods of solving linear TPBVP (Equations (3.7.5),
(3.7.5)), for example, the transition matrix algorithm might have been
numerically more stable. However, other difficulties arose because these
methods require several forward integrations in one iteration. This
means calling the subroutine DIFEQ (the subroutine to generate the right
side of the differential equations and its derivatives) many more times.
This increases the computing time enormously.
With a step size small enough for the integration to be stable,
the storage requirements, computing time, and, therefore, the cost of
computing increases considerably. Even if these factors are absorbed,
it may still be necessary to try several initial guesses of the adjoint
variables to get the method to converge.
In view of the above problems, it was concluded that the
standard method of quasilinearization was not suitable for the human
motion problem and so should not be pursued further.
79 ,
3.8. Steepest Descent Methods for Solving
the Minimum-Time Kip-Up Problem
Three different formulations of first-order steepest descent
methods were used after the method of quasilinearization was unsuccess
ful in solving the minimum-time kip-up problem. These formulations
differ from each other in the construction of the cost functional,
handling of the terminal constraints, treatment of the control incre
ments and in the method of adjusting the final time. The basic features
of these three formulations are described below.
Formulation 1
a. The cost functional is the final time. The terminal errors
and the cost functional are reduced simultaneously.
b. The adjustment of the final time is done by extending or
truncating the final end of the trajectories.
c. The control functions take the form of a sequence of constrained
and unconstrained arcs. Improvements are made at the uncon
strained parts only, including the junctions of the constrained
and unconstrained arcs.
The method is based on the works of Bryson and Denham [27,28]
and Bryson and Ho [32].
Formulation 2
a. The cost functional is the sum of a scalar representing the
final time and a norm of the terminal error.
b. A change in the independent variable t is introduced by defin
ing the transformation
80
t = a t
where a is a constant, or
do
= 0
dT
so that a is treated as an additional state variable. The final
time t is directly proportional to a when is held fixed.
Long [29] used this transformation of the independent variable
to convert free end point TPBVP to fixed end point TPBVP for
solution by the method of quasilinearization.
The cost functional is reconstructed as
2 6 2
J = K O' + E K.S., K ,K. ,. . ,K > 0.
o.^xi o i 6
1=1
No terminal constraints were introduced in this formulation
since they (the §/s) are included in the cost functional,
c. The control functions are assumed to be free to change in any
direction while computing the gradient of the cost function.
That is, the control constraints were not considered when com
puting the gradient of the cost functional used to find a suit
able increment in the control and o'. The control constraints
were imposed in the next iteration during forward integration
of the state equations. When the computed control violated
a constraint in a subarc it was set equal to its limit, the
constraint function on the subarc. This approach for treating
constraints on the control has been used by Wong, Dressier and
Luenburger [33].
81
Formulation 3
This formulation consisted of the features (a) and (b) of
Formulation 2 and (c) of Formulation 1.
The derivations of the numerical algorithms for these formula
tions are now presented. The basic concepts on which these algorithms
are based are available in the literature [28,32]. The results are
derived in a manner suitable for analysis of the outcome of the numer
ical computations, and, thus, the derivations presented here are slightly
different from those found in the literature for any gradient method
approach to the computation of optimally controlled motion.
3.8.1. Derivations for Formulation 1
Suppose we have a continuous nominal control u^(t) and u^it)
and a nominal final time t These control histories have some parts
lying on the constraints S^(X), i = 1 or 2, and the remaining parts lie
away from the constraints. The parts lying on constraints will be
called "constrained arcs" and the parts lying off the constraints will
be called "unconstrained arcs." The interesections of constrained and
unconstrained arcs will be called "corner points." A nominal guess for
the control variables consists of specifying the corner points and
values of the control at unconstrained portions. On the constrained
arcs control variables are generated from constraint functions. An
initial choice of the control history and the final time will not, in
general, satisfy the boundary condition and will not do so in minimum
time either. One can improve upon the trajectories in the following way.
At first, we establish how a particular state variable X1 (the
th
i
component of the state vector X) changes at the final time due to
82
a small change in the control history and the final time. For that
purpose, let a cost functional be defined
J = X1(tf) (3.8.1)
Let \1(t) be an arbitrary time-varying vector of dimension six.
Since the system satisfies Equation (3.1.1), the final value of the
state variable will also be given by
*f T
j' = XX(t ) + f X1 (f-X) dt (3.8.2)
f do - -
If the control variables u (t) and u (t) change by a small
X j
amount 6u (t) and fiu (t) there will also be a small change in the
X i
state variable X(t), denoted by fix(t), throughout the trajectory.
It is clear that these changes in the control, denoted by fiu(t) and
in the state variables denoted by fix(t), will not be independent of
each other. Apart from changing the control history, an increment to
the final time t by a small amount 6t^ and small increment fiX(O) to the
initial state vector are also prescribed. The first-order change in J
due to the changes in the control and the final time is given by
Aj' = X1(tf) dtf + {Sx1 (\X)T 6x}t + (XX)T fiX(O)
tf tf
+ f {aV f, + (x1)} fix dt + f qV f, fiu dt.
(3.8.3)
The fiu is chosen in the following way: For the unconstrained parts,
fiu is completely free. The parts presently on the constraints will
remain on the constraints for the same periods of time as before.
83
For these portions, the change in control is given by the shift of the
constraint due to state changes according to the relation
6u (t) = S3 $X(t) (3.8.4)
1 "
Let c} denote all those portions of the trajectory of the
control u^(t), i = 1 or 2, which lie on any of the constraints s{ or
2 o
S.. Also, let C. denote all those portions of the control u. which
l i i
do not lie on a constraint. If the expressions for on the con
strained arcs given by Equations (3.8.4) are substituted in Equa
tion (3.8.3) and the integration of the last term of the right side of
o 1
Equation (3.8.2) is split into integrations over the intervals ,
o 1 , .
C and C one obtains
fij' = xx(t )dt + {ax1 (xV ax} + {(xV fix}
1 1 Z~Vf t=o
tf
+ Jo {(^1)T i'X+ ('X)} dt
+ f (\X)T f,u dt + f (xV #, S3 6X dt
o 1 j 1 X
(X1)T f, u dt + f (xV f, sJ 6x dt
" u2 2 1 2 X
If Xx(t) is computed such that
Xj(tf) = 1 for i = j
= 0 for i 4 j
where x{ = j element of X*
and
; i
X
l
X
1
(3.8.5)
(3.8.6)
(3.8.7)
84
where
6=0 on C?, 6-1 on
1 11 1
62 = 0 " 2' *2 = 1 n C2
and
u1 and u used in Equation (3.8.7) are computed only when
x,x 2X
the controls u^ and u^ are on a constraint and u^ and u^ are replaced by
the constraint expressions S^(X), one obtains
6j' = X1(t.) dt + (X1(0))T fix(0) + f (X1) f, 6U dt
f f - c ux 1
+ £'U S dt
= f1 (X (t_),u(t.))dt. + aX(0))T flx(0) + f (X1)1 f, flu.
- I I I o - u. 1
i.T
dt
+ lf0(X)T f,u 6u2 dt (3.8.8)
C2 2
where fi is the i^h element of the vector f(X,u).
Equation (3.8.8) is the desired expression for the change of
the state variable X1 at the final time t^ due to (1) a small arbitrary
increment §u(t) given to the control variable u(t) over the unconstrained
portions, (2) a small increase dt of the final time t^, and (3) an
arbitrary small change 6x(0) of the initial state vector X .
Similarly, one can find the expressions for the change in the
terminal value of any other state variable It can be seen that if
one constructs the (6 x6) matrix
85
R(t)
-[V
(t),
2
X (t),
X3(t)
4
X (t),
x5(t), x6(t)
]
so that R(t) satisfies
and
R(t^) = I (6x6 unit matrix)
R(t) =
6 f
l-u.
R(t)
(3.8.9)
(3.8.10)
where the meaning of the various terms in the parentheses of right
side of Equation (3.8.10) is the same as that in the same terms appear-
in Equation (3.8.7), then the change in the terminal value of the state
vector is given by
§X(t )=f dt + RT(0)SX(0) f RT f, Su dt + f RT f, §un dt .
f f ho u 1 u_o u 2
(3.8.11a)
If we choose that 6x(0) =0, we obtain
Sx(tf)
f dt + f RT f, 5u dt + f RT f, 6u dt.
- f J o -u 1 J o -u 2
C1 1 C2 2
(3.8.11b)
Following the method prescribed by Bryson and Ho [32], it will
now be attempted to make improvements in the terminal errors given by
§ = X(tf) Xf
and at the same time reduce the final time t^. Thus, since it is being
sought to minimize t^, one would maximize -dt^, or minimize dt^ with
respect to 6u_^, 6u^ subject to the constraint
6x(t ) = f(t )dt + f RT f, 6u dt + f RT f, §u dt
-f f f v o u 1 J o u 2
C, 1 C 2
(3.8.12)
86
with
6x(t ) = 6xp(t )
- f f
where OX (t^,) is a chosen decrement in the terminal error such that
§u maintains.the first-order approximations.
In this incremental minimization problem, the incremental cost
functional, dt^, to be minimized, and the constraints are linear in the
incremental control parameter §u. Such a problem does not have an
extremum. However, since these are linearized relations obtained from
a nonlinear system, the increments 6u §u and dt should be small for
1 ^ X
the first-order approximation to be valid. To limit the increments
6u 6u and dt the following quadratic penalty term to the incre-
X Zj X
mental cost dt^ is added:
I b dtf + \ So 6U1 Wl(t) dt + I Jo 6US W2(t) dt
(3.8.13)
where b is a positive scalar quantity and W (t) and W (t) are positive
X -J
scalar quantities specified as functions of time.
Adding these quantities to the cost functional and adjoining
the constraint relations (3.8.12) to the resulting expression by a
multiplier v (a six-dimensional vector), the following problem is
obtained.
Minimize wrt, &u $u and
X Z X
dt +
f 2
b dt^ + i f 6u7 W dt + i- f du^ W0dt + vT^f (t^)dt^
'f+2Jo VU1 Tl+2Jo 2 V
C1 C2
+ f Rrf, 6u dt+ f RTf, *u dt 6XP(tJ
C U1 1 C U2 2 f
}
(3.8.14)
87
If the derivatives of this functional with respect to 6u^, &u9, and
dt^ are set to zero, one obtains
dt = ^ [1 +vT f(t )] ,
f b f
(3.8.15)
and
s =
6u =
2
W.
iT
*2 2
Rv on C.
Rv on CL
P A
(3.8.16)
for an extremal.
Using Equations (3.8.15) and (3.8.16) in (3.8.12), one obtains
§XP = i fil + vTf3 -IMv (3.8.17)
b - M -
f
or
where
v = -
1 T
I. +- f f (tj
\jfiji b - f
n T -IT n T -IT
IM= R f, V,' f, Rdt+J R f, Wo R dt.
M 'io u 1 -u o -u 2 u
-i r
6*P
-1 T
(3.8.18)
C1 1 1 C2 2 "2
(3.8.19)
The value of 6xJ (t^) the desired change in the terminal values of the
state variables, may be chosen as a decrement in the terminal error S.
>XP(t ) = e. §.
if li
(3.8.20)
where
0 < e. £ 1
i
. th
6xP(t ) is the iL11 element of ^XP.
if
and
88
The steepest descent algorithm for Formulation 1 can now be
described as follows:
1. Guess a nominal control history u (t) and u (t) and a final
time t^.
2. Integrate the state equations forward with the nominal control
and the nominal final time with the initial values of the state
vector given by X(0) = X Store X (t). Compute and store
- -o -
Find the norm
INI
3.
4.
6.
7.
for some positive i = 1,...,6 (to be specified).
Save the controls u^ and u^ in another variable COLD, the corner
times N in the variable NOLD and the final time t^ in TFOLD.
Set R(t^) = I, the (6 X6) unit matrix. Integrate backward
T
Equations (3.8.10). At the same time compute and store R f,^
T ^
and R f, .It is not necessary to store R(t).
- u
2
Select the positive constant b and the positive quantities
W (t) and W (t). Unless there is some special reason, W (t)
i. Z JL
and W (t) may be taken as a positive constant equal to W (to be
z
specified). In such a case, the storage for W^(t) and W (t)
will not be needed. Select e^, i = 1,2,...,6 such that
o < e. ^ l.
Compute I^ using Equation
Compute v, §u^, Su^, anc* dt
and (3.8.16). Note that §u
unconstrained parts C and
(3.8.19).
from Equations (3.8.18), (3.8.15),
and u^ are defined for the
C only, and are to be computed only
on those parts.
89
Change the final time t^ to t + dt^. If the integration step
size is h, this means that the total interval is to be increased or
decreased at the final end (depending upon whether dt^ is positive or
negative) by dt^/h, rounded to the nearest integer number of integra
tion steps. In the case of a positive dt^, the controls u^(t) and
u (t) in the interval TFOLD ^ t ^ t (new) will be given by extrapola-
tion of the existing curves. If a control function is on a constraint
at the final time, TFOLD, it should remain on the same constraint in
this period and is to be generated during the forward integration in the
next steps.
8. Set the control u. = UOLD. + Su., i = 1 or 2, for the uncon-
l li
strained portions. Integrate the state equations forward with
the new control. The corner times N are to be generated again
during the forward integration as described graphically in
Figure 11.
Figure 11 shows the different cases which may arise during the
forward integration. It is seen that a new corner time is generated at
the point where the unconstrained u^ + 6u^ curve, extrapolated if neces
sary, meets the constraint. This part of the forward integration makes
the programming complicated with many logical program statements.
9. Find the errors in the boundary conditions § and find the norm
j| §^jj as in step 2. If this norm is less than that of the
previous iteration go to step 3 to continue the iterations.
If it is not reduced, then do the following:
a. If dt^ is too large, increase b, and go to step 9c below.
