is
{yfztltl‘i
. .55.... ii!

JI...)).KI! P
$5.11!!) .3
all... 5...?- If.»

E) 53...!) I?!PI..R:I,¢I

Rik: .gén‘li‘i
(1:)!-

 

!! {P1101 ;.
Sent!!! 5‘ '05:...55.
.21.,tlrvrflfilrl

 

4532(1)
1 :I \(p t?
((

 

03...! 11,1}. xr:

. 11.5.; tr: pvt:
5.
. .i 4
5.2.1 2.; t ,M. 135:}... ..
4 <5. .15.

51...)? £g_u\ti€§‘¥=4£
. E?!§E:.!.IE3 ‘r‘frrl‘l
.- (Iv)! ca‘!:!§.l<$..\.l
magi-14.25... 1.55:2...
It}!,(llkl(<!yt(u\lttl¢('ti(li :-

. girls-v)! I...

 

. :érlllkneltrlscc $5553.91?
. .‘(rlgitirxitllli ll.

Ectlrir‘n; .
(fi‘.(rv/l§(fir:r(ra.
5.2.; 1%.... 5:35:53
5275?}:

 

5.33."? 2.. «a

 

r§,tfir..$§.
. . . S.¢‘¢rl.¢3x. I «(11121Ipitré..rr¢(¢-i\tiié}.(:i \iiéiktlffilslrisf‘lll

.. 1 2... .prII :1. .3 39!}: .. :1 ‘52). isgilsg‘isrioiil‘lfli

. : ... x 5.13.5! .44 .31....)(5I‘tu¢olilloltlllril

/ , AI(.€I....II..\.51

. s... ‘ . ., : . . ‘ .. ‘ .

1 1 ‘ ‘

/7 IT." l k... "I." ' ‘
/”"C / ~> J

“C ESTES“

1|quliiiflllfliilfl'Iililiiilliiiiilliiifl

0790 3226

This is to certify that the

dissertation entitled ‘

OPTIMAL CONTROL IN
CONSTRAINED MULTI-BODY MECHANICAL SYSTEMS

presented by

Moon Suk Suh

has been accepted towards fulfillment
of the requirements for

 

Ph. D. degreein Mechanical Engineering

 

 

MAL/AM '

Major professor I

Date November 8, 1989

 

MSU is an Affirmative Action/Equal Opportunity Institution 042771

 

 

.r
LIBRARY
Michigan State
University

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

JU
”.12 a. 20
, 5 j 3' it

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
c:\cIrc\datedm.pm3-p.t

OPTIMAL CONTROL IN
CONSTRAINED MULTI-BODY MECHANICAL SYSTEMS

By
MoonSuk Suh

A DISSERTATION

submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

1989

b054QA7

ABSTRACT

OPTIMAL CONTROL
IN CONSTRAINED MULTI-BODY MECHANICAL SYSTEMS

By

Moon Suk Suh

Multi-body mechanical systems modeled in non-minimal coordinates have vari-
ous advantages over systems in minimal sets of generalized coordinates. For example
a system in non—minimal coordinates is simpler to model and important system
measures(descriptor variables) are directly available. Unfortunately the many advan—
tages of the nonminimal forms are nearly overcome by disadvantages due to the
singularities of the system. Especially the standard forms of optimal control are not
convenient for the systems represented in non-minimal coordinates (descriptor equa-
tions).

The goal of this dissertation is to derive and test new and flexible forms of the
necessary conditions of optimal control appropriate for multi-body systems in non-
minimal coordinate sets. In order to derive the necessary conditions for optimal con-
trol, the problem is posed directly on an implicit system represented in non-minimal
coordinates sets without reducing the second order form to the first order form. The
approach follows the standard method of forming an augmented cost functional
wherein the system is adjoined to a nonnegative by introducing adjoint variables. The
second order implicit form necessary conditions for optimal control problms on
multi-body mechanical systems models of non-minimal or descriptor form are

presented.

ACKNOWLEDGEMENTS

I would especially like to express my deep gratitude to my adviser Dr. Joseph
Whitesell, not only for his significant contributions to this dissertation, but also his

friendly advice and steady encouragement.

I would also like to thank the members of the Guidance Committee, Dr. Alejan-
dro Diaz, Dr. R. Lal Tummala and Dr. Byron Drachman for their thoughtful com-
ments, inspiration and continued interest during the course of this research.

Special thanks is due to Dr. Hassan Khalil for the informative and helpful
suggestions concerning singular perturbation theory.

This thesis is dedicated to my late father Chung Kook Suh and my mother In
Keun Suh who supported me all the way to where I am today.

And, of course, a great debt is owed to my lovely wife Hwa Sook for her con-
tinuous assistence and patience.

I would like to thank my cute two sons, Chul Won and Jang Won and my

brother and sister.

Finally, I deeply appreciate the support and consideration of my elder cousin,

Min Suk Suh.

iii

TABLE OF CONTENTS

 

 

 

 

 

 

 

 

 

LIST OF TABLES vi
LIST OF FIGURES vii
Chapter 1: INTRODUCTION 1
1.1 Motivation 1

1.2 Literature survey 2

1.3 Organization of the dissertation 5
Chapter 2: TECHNICAL PRELIMINARIES 7
2.1 Multibody dynamics on constrained Lagrangian systems ................ 7

2.1.1 The choice of generalized coordinates ............................. 7

2.1.2 Cartesian coordinates 11

2.1.3 Constraint equations 13

 

2.1.4 The formulation of the equations of motion

 

 

 

on constrained planar systems 15

2.1.5 Generalized coordinate partitioning ................................. 17

2.2 Optimal control general 20
2.2.1 Typical optimal control problem 20

2.2.2 Optimal control on descriptor system .............................. 24

2.2.3 Mathematical concepts for optimal control ...................... 26

Chapter 3: THE SECOND ORDER IMPLICIT FORM NECESSARY
CONDITIONS ON OPTIMAL CONTROL IN

CONSTRAINED MULTI-BODY MECHANICAL SYSTEMS ........

3.

p—A

Introduction

 

3.2 The second order implicit form necessary conditions

 

on optimal control in descriptor systems

Chapter 4: NUMERICAL VHEW

 

4.1 Introduction

 

 

4.2 The Gear stiff integration method

4.3 Formulation of corrector for system dynamic equations and

adjoint equations

 

Chapter 5: EXAMPLES

 

 

5.1 Two constrained sliding masses problem

5.2 Simple pendulum problem

 

5.3 Simple two degree of freedom manipulator ..................................

Chapter 6: SUMMARY AND CONCLUSIONS

 

LIST OF REFERENCES

 

32

32

32

48

48

49

51

56

56

66

78

88

9O

LIST OF TABLES

 

‘/
Table 5.1 The necessary conditions for J = if 142d: 70
'0

Table 5.2 The necessary conditions for J = %f( u2 + 13+ 1% )dr ......................... 74

vi

LIST OF FIGURES

Figure 2.1 The two degree of freedom manipulator model .............................
Figure 2.2 The Cartesian generalized coordinates for the i—th body .................

Figure 2.3 Revolute joint coordinates

 

 

Figure 3.1 Relation between Sq/ and 8q(tf)

Figure 5.1 Two constrained sliding masses

 

 

Figure 5.2 The optimal control for case (1a)

Figure 5.3 The performance measure for case (1a)

 

Figure 5.4 The optimal trajectory for case (1a)

 

Figure 5.5 The Lagrange multiplier for case (1a)

 

Figure 5.6 The optimal control for case (1b)

 

 

Figure 5.7 The performance measure for case (1b)

Figure 5.8 The optimal trajectory for case (1b)

 

Figure 5.9 The Lagrange multiplier for case (1b)

 

Figure 5.10 Comparison between non-minimal approach and exact solution

 

Figure 5.11 Planar simple pendulum

 

Figure 5.12 The optimal control for case (2a)

 

 

Figure 5.13 The optimal trajectory for case (2a)

Figure 5.14 X-direction reaction force for case (23)

 

Figure 5.15 Y-direction reaction force for case (2a)

 

 

Figure 5.16 The performance measure for case (2a)

vii

39

57

59

59

60

6O

62

63

63

64

65

66

71

71

72

72

73

Figure 5.17 The optimal control for case (2b)

 

 

Figure 5.18 The optimal trajectory for case (2b)

Figure 5.19 X-direction reaction force for case (2b)

 

Figure 5.20 Y-direction reaction force for case (2b)

 

Figure 5.21 The performance measure for case (2b)

 

Figure 5.22 The optimal controls for manipulator model .................................
Figure 5.23 The optimal trajectories for manipulator model .............................
Figure 5.24 The first reaction forces for manipulator model ............................
Figure 5.25 The second joint reaction forces for manipulator model .................

Figure 5.26 The performance measure for manipulator model ..........................

viii

76

76

77

77

78

85

85

86

86

87

CHAPTER 1

INTRODUCTION

1.1. Motivation

Pontryagin Principle [12] brought a significant advance in the general theory of
optimal control problem. The theorem provides a set of local (necessary) conditions
which an optimal control must satisfy in order to minimize a given performance cri-
terion. The result says that if a control function minimizes the objective functional,
its values at each time must also minimize the Hamiltonian. This point of view gives
broad guarantees for existence of solutions, but it is non-trivial to specialize it to
specific cases since the optimal control problem has been treated in a quite general
form. For example, the standard forms of optimal control are not convenient for the
second order implicit form differential-algebraic systems which can result when
multi-body mechanical systems are represented by equations involving non-minimal

coordinate sets.

In this dissertation a non-minimal set of Cartesian coordinates together with a
Lagrange multiplier set is chosen to construct mechanical system dynamic equations.
In order to derive the necessary conditions for optimal control we directly apply the
optimal control problem to the implicit system of non-minimal coordinates sets
without reducing the second order form to the first order form. There are several
motivations for this.

0 A system in non-minimal coordinates is simpler to model

0 In constructing performance measures it is convenient to have access to a large

set of system (descriptor[25,26]) variables such as local body positions, velocities

and reaction forces. For example, we might wish to construct a performance
measure which depends on elements of the Lagrange multiplier set in order to
minimize reaction forces which can contribute to component wear or or other

undesirable effects.

0 While in some cases we might be able to transform from a set of Cartesian coor-
dinates to some minimal set of Lagrangian generalized coordinates, it may not be
convenient or possible to do this in all cases. In general, it has not been shown

how to find minimal Lagrangian models over finite time intervals.

0 Deriving information from minimal Lagrangian models may be difficult since a

complex set of transformations which may become singular may be needed.

0 There are advantages in assembling and computing the necessary conditions of
optimality if the problem is posed directly on the second order implicit system.
For example if the necessary conditions for standard state space form are used in
the first order explicit system, the necessary conditions not only require matrix
inversions and products for each value of state vectors, but also they involve
computing the Jacobian of the matrix—vector products, which is difficult to pro-
gram and inefficient to compute. In addition the adjoint equations will have quite
different forms as compared to the dynamic equations thus making automatic

assembly difficult.

Unfortunately the many advantages of the non-minimal forms are nearly over-
come by their disadvantages. Computations on them are difficult. When the descriptor

system is applied to the control problem, we may encounter fast transients.

1.2. Literature Survey

A conventional approach for analyzing multi-body mechanical systems is to use
Lagrangian methods to formulate the system equations. Readers can to refer several
standard texts [1,2,3] for fundamentals of dynamics required for mechanical system
dynamic analysis. An important issue is coordinate choice. There are many dynamic
analysis programs which employ Lagrangian coordinates whose coordinates consist of
a minimal set of independent generalized coordinates[49]. Since the equations of
motion are written in terms of the minimal set of coordinates, a small number of
governing equations of motion is obtained, which may be integrated by standard
numerical integration techniques. However the reduced set of equations is highly
nonlinear and complex. In addition some important variables may not appear in the
method. For example, the system associated with Lagrangian coordinates does not
directly evaluate the Lagrange multipliers which are associated with the forces of
constraint or reaction forces. Typically, the computation of reaction forces is done

after the motion of the system is evaluated.

An alternative choice of coordinates is a non-minimal set of Cartesian coordi-
nates. This approach was used in computer programs [4,5,6]. The Cartesian coordi-
nates typically lead to large sparse sets of mixed differential-algebraic equations. For
example, multi-body systems are often represented by mixed differential-algebraic
equations involving non-minimal coordinate sets such as Eq. (1.1) below [44].

M<q>ii — <I>qTx = f<q.<'1.u) (1'1)

<I>(q) = 0
where M is a n x n matrix, the vector q(t) is an n dimensional state vector and u(t)
is an k dimensional control vector. The vector 7» with dimension m denotes the well
known Lagrange multiplier set. In particular, in system (1.1), the term <1)qu may be
interpreted as the reaction forces which result when m independent algebraic con-
straint function (I) are imposed on the system. It is obvious that m is zero and that the

system consists of just differential equations when the dynamic system is

unconstrained. We note that if system (1.1) is represented in the first order form
Ei(t) = a(z,u,t)

where z = [ qT, (1T ,KT ,AT ]T, the matrix E is, in general, singular. Such systems are
referred as singular systems, differential-algebraic systems or descriptor systems.

The models such as the system (1.1) of Cartesian coordinates are simple to
create and important variables such as Lagrange multipliers appear explicitly as
shown in system (1.1). In addition the Jacobian matrix is simple to form. While it is
easy to formulate large systems of mixed differential-algebraic equations of motion
using the Cartesian coordinate approach, solving the mixed differential-algebraic
equations is difficult [12]. Systems of the mixed differential-algebraic equations can
not normally be solved by conventional numerical methods that are used to solve
ordinary differential equations [9]. Nikravesh [7] reviews three integration algorithms
that are used to solve mixed differential-algebraic equations. These are the direct
integration method, the constraint stabilization method [8], and the generalized coor-
dinate partitioning method [5]. But there are often numerical difficulties with these
methods [11,25]. Gear [11,12] has shown that his stiff integration algorithm[15,l6]

can be used to solve only certain kinds of descriptor systems.

Descriptor systems have recently attracted the attention of many authors.
Knowledge and understanding for descriptor systems can be obtained from references
[28,30,31]. Some theoretical aspects of the topic have been given in [25,26], which
discuss singular systems with the aid of the matrix generalized inverse called the
Drazin inverse. Pandolfi[34] studied the controllability and stabilization for the sys-
tem using the Drazin inverse. Important research was done by Cobb[1983] who
proved that the optimal control for linear descriptor systems can be found by solving
a reduced order Riccatti equation through a sequence of coordinate transformations

without restrictions on initial conditions. He noted that fast transients may appear in

the optimal control unless special case is taken. Lovass [29] has considered the prob-
lem of unconstrained optimal control of continuous-time singular systems. Recently,
Verghese [28] extended the results of Lovass [29] to the case where the final time tf
is not fixed. Wu [32] derived the necessary conditions, called the generalized max-
imum principle, for optimal control of descriptor systems in a more general form by
making use of a method of the calculus of variations for a fixed final state with a
fixed time. Barman [50] performed the design sensitivity analysis and optimization on
multi-body systems. In [33], the problem of free final time optimal control of

descriptor systems has been considered.

In this dissertation, the optimal control problem is posed on directly on the
second order equations of system (1.1). The approach follows the standard method of
forming an augmented cost functional wherein the system is adjoined to a nonnega-
tive performance measure by introducing adjoint variables [35]. After demonstrating
the regularity of the constraint set(i.e. system (1.1)), we use the standard Lagrange
multiplier theorem of constrained optimization to guarantee the existence of adjoint

variables for the problem.

1.3. Organization of the Dissertation

The dissertation is divided into six chapters.

In chapter 2, there is a short review of the multi-body mechanics. The chapter
includes general concepts of the generalized coordinates, Cartesian coordinates, the
formulation of Lagrange’s equation on constrained systems, and the method of gen-
eralized coordinate partitioning. Even if the field of multi-body mechanics is not new,
the chapter will help to understand the overall problem in the dissertation. It is gen-
erally recognized that many of the problems of optimal control theory are approached

through the calculus of variations. The chapter also reviews the conventional optimal

control problems by the calculus of variations. Some disadvantages which occur
when the conventional optimal control approach is applied to the dynamical system
are stated. Important mathematical concepts for developing the necessary conditions

for optimal control follow.

The chapter 3 is the most important contribution of this dissertation. After dis-
cussing the regularity of the descriptor system, the second order implicit form of
necessary conditions will be introduced by posing the optimal control problem
directly on the second order descriptor system. In particular the performance meas-
ure minimized can include system descriptor variable such as dynamic Lagrange mul-
tipliers.

The chapter 4 presents a numerical method for the system dynamical equations
and the adjoint equations of optimal control which both are the mixed-differential
algebraic. The formulation of corrector for the system dynamic equations and the

adjoint equations are described in the chapter.

Chapter 5 contains several examples which are a constrained spring-mass prob—
lem, pendulum problem and simple two degree of freedom manipulator model. The
computer simulation results confirms the validity of the necessary conditions which
are obtained although it is noted that fast terminal transient appear in some of the

results.

Finally in chapter 6, summary and conclusions of the dissertation are given, and

suggestions for future research are presented.

CHAPTER 2

TECHNICAL PRELIMINARIES

In this chapter technical background material is presented.

2.1. Multibody Dynamics on Constrained Lagrangian Systems

A multibody mechanical system is a collection of bodies (rigid links or flexible
links) whose relative motion is limited by joints or constraints. Dynamic analysis of
a mechanical system involves determining the time history of position, velocity, and
acceleration of one or more of the bodies due to the action of external and internal
forces. The equations of motion are differential or mixed differential-algebraic equa-
tions. This section represents a review of the equations and principles of multibody
mechanics. For simplicity we will describe only the planar case, however the results

of the this thesis are not limited to this case.

2.1.1. The Choice of Generalized Coordinates

Any set of parameters qj that completely specify the configuration of a system
are called generalized coordinates. A proper set of generalized coordinates whose
number equals the number of degrees of freedom of the system and which are not
restricted by the constraints will be called a set of independent generalized coordi-
nates. A set of n independent generalized coordinates can be considered to be n
components of a vector in an n-dimensional configuration space (or C-space) in

which each point represents a possible configuration of the system. The choice of a

 

 

set of generalized coordinates for description of a system is not unique. The coordi—
mates for defining a mechanical system can be chosen from two large categories.
One is a minimal set of Lagrangian generalized coordinates. The other is the non-
minimal Cartesian coordinate set.

The Cartesian approach was introduced in computer programs in the later
1970’s and early 1980’s[4,5,6] while computer implimentation of the Lagrangian
approach appeared in the early 1970’s[45,46,47]. The simple two degree of freedom
manipulator model is shown in Figure 2.1 to illustrate the two systems.

The symbols m1 and m2 are the masses of link 1 and link 2 respectively. The sym—
bols 11 and 12 are the moments of inertia about the centers of each mass and the g is
gravity. The length of the links is defined as 211, 212 for simplicity. Since this model
has two degree of freedom, the minimal number of coordinates is two. If the angular
orientations 61, 92 of the link are chosen as the Lagrangian generalized coordinates as

shown in the Figure 2.1, the two second order differential equations can be written as

[ m111f+ 4m21% + 2m21112cos(02 — 91) + J1 19,
+ [ m21§+ 2m2l112cos(02 — 9,) + 12 192
— 2m211l2(Ol + 92x92 — Ol)sin(02 - 9,)
+ mlgllcosOI + m2g(lzcos(92 + lecosel) = 0
and
[ "121% + 2m21112cos(62 — 91) + 12 191 + ("121% + J,)92

+ 2m21112025in(02 — 9,) + nglzcosez = 0

 

 

 

Figure 2.1 The two degree of freedom manipulator model
In order to formulate the model in terms of Cartesian coordinates the angular orienta-
tions 61 and 02, the coordinates x1 and y] of the center of mass of link 1, the coordi-
nates x2 and y2 of the center of mass of link 2 respectively are chosen as Cartesian
coordinates. The equations of the motion become the following six differential equa-

tions f and four algebraic constraint equations d).

 

That is,

mlic'l + Al — A3
mri’i + 71-2 “ A4 + ”‘18
1161+1110~1+ A3)Sin91-(7v2 + 3.4) c0361 ] - u1
”12552 + As
m7y'2 + 7L4 + ng
1262 + 12( MsinOZ — Mcost) — uz

In.)
II

II
O

 

 

and

‘X1 ‘l' 1100861
—y1 + llsinOl
(I) = x1 + 1100591 — x2 + 1200592 =

yl '1' 11811191 — y2 + leIIIGZ

Though the advantage of Lagrangian coordinates is that a small number of
governing equations of motion is obtained, which may be integrated by standard
numerical integration techniques, there are disadvantages. The internal representation
of the equations of motion of the approach is very complex and some important vari-
ables such as Lagrange multipliers which are related to generalized reaction forces in
the system do not appear in the system equations. In addition linearized equations are
difficult to form. A major disadvantage is that a minimal set Lagrangian generalized
coordinate can not always be found in a practical sense. Furtherinore the necessary
coordinates for a given system may change as the system moves. On the other hand
the models in Cartesian coordinates are simple to obtain, introducing special effects
are easy, important variables such as Lagrange multipliers appear explicitly in the
system equations and linearized equations are simple to form. The major disadvan-
tage is the difficulty of the integration because Cartesian coordinates leads to a mixed

differential algebraic equations.

 

 

 

11

2.1.2. Cartesian Coordinates

For simplicity we will illustrate the Cartesian approach for planar systems. In
order to specify the configuration of a planar mechanical system, consider the xy
coordinate system shown in Figure 2.2 as a global reference frame. As a local refer-

ence frame, we define a body-fixed coordinate system Em,- embedded in body 1'.

111

 

 

 

Figure 2.2 The Cartesian generalized coordinates for the i—th body

Usually a local reference frame in each link of the mechanism or machine is
located at the center of mass. Since a local reference frame 12m,- moves with its part,
the origin of the local reference frame can be located by specifing in the global coor-

dinates Rl- = [x,y],-T. and angle 6,- of rotation of this system relative to the global

12

coordinate system. This angle is considered positive if the rotation from positive x
axis to positive é,- axis is counterclockwise.

If a point P,- on body i can be located from the origin of the I‘M-t],- axes by the
vector rf , the coordinates of the point P,- with respect to the body-fixed coordinate
gin,- are £5) and 11fJ . The body-fixed components of vector rf can be written as
r? = [amp]? . Since P,- is a fixed point on body 1' , g? and 1]? are constants, and
therefore r’{’ is a constant vector regardless of the motion of the body. The global
xy components of vector rf should be expressed as a function of the rotation angle 6,-
because the elements of rf vary when body i rotates. It is clear that the components

of rf are not necessarily constant, since body i may be in motion.
The relation between the local and global coordinates of point Pi is
R]? = R, + rf’ (2.1)
where

rf’ = T-r’f’ (2.2)

I

Here T; is the rotational transformation matrix for body i. T,- is expressed by Eq.

(2.3)

__ cosO —sin0
TIT [sine cos9]i (2‘3)

Equation (2.2) provides the relationship between the local and global components of

vector rf
Equation (2.1) in expanded form can be written as
BY = k] + [cosO —sin9] [12],)
i i srnO cost) in i
or

xf = x,- + éfcosei — nfsinOi

 

13
y? = yi + firsinei + nfcosei

Equation (2.1) uniquely defines the position of any point P,- on body i , in terms of
its Em,- body-fixed coordinates r’f and variables x,- ,y,- , and O,- that define the position
and orientation of the 13m,- local reference frame. Thus, xi, yi, and 9,- serve as Carte-

sian generalized coordinates for body 1' , denoted by the vector

[x]
q,-= y
6i

For a planar mechanical system with n rigid bodies, the generalized coordinate vector
is the 3n vector

q = [ (ring... q}; IT

= [ x1,}’1961,x2,)’2,92, ----- a memen ]T

where q denotes the vector of generalized coordinates for the entire system.
2.1.3. Constraint Equations

A mechanical system often consists of links or bodies that are interconnected by
one or more kinematic joints. In the Cartesian coordinate formulation the joints are

represented by algebraic constraints between the bodies that make up the system.

In the planar system there are two principal types of joint which are the revolute
and translational joints. Figure 2.3 shows schematic representation of a revolute joint
connecting to bodies i and j with body-fixed coordinate systems Em,- and éjnj ,
respectively. The origins of these local reference frames are located in the global

reference frame xy by vectors R,- and R-, respecu'vely.

 

 

14

 

 

Y A
711' rp
rl-j
onj
.............. 6...... nj

(xi yi) l

E'

é, rfl' J

............................................. 9 1

R1. (xi y)

R.
I

X

Figure 2.3 Revolute joint coordinates

Let point P9- on body 1’ be located by a body-fixed vector rij and point P}; on
body j be located by a body-fixed vector rJ-i . Since points Pi]- and Pji are connected
by a vector rP , the loop closure relation of the position vectors gives

R,+r,j+r”-rfi—Rj=0 (2.4)
Since points Pa and Pi; can be considered to be coincident rotational joint between
bodies i and j ,

H’=0

15

Thus the constraint equation (2.4) for the revolute joint yields
R‘- + rij — r17- Rj = 0 (2.5)
Let (xbyi) and (xj,yj) be the global coordinates of body-fixed coordinates gm,- and éjnj

respectively, the constraint equations for a revolute joint between bodies i and j are

obtained from Eq. (2.5) as
(I)x = x,- + fifcosei — nfsinei — xj - éfcosej + nfsinej = 0 (2.6a)

_ P - P P - P _
(by — yi + §i51n6i+ T1,- cosOi — yj — Q 31an — Tl; cost — 0 (2.6b)

2.1.4. The Formulation of the Equations of Motion on Constrained Planar Sys-

terns

In the Lagrangian formulation [1,2,3] the general form of the equations of

motion can be written by z

1 [fl] _ .31. = Q, (2.7)

where T = kinetic energy of the system
q,- = generalized coordinates
('1,- = generalized velocities
Q = generalized forces with the particular coordinate q,- being considered.

The Q,- are a summation of conservative and nonconservative applied

force effects.

The applied forces may be arbitrary functions of displacements and

velocities.

16

Examples of applied forces would include gravity, springs, and viscous

dampers etc.

Equation (2.7) is valid for constrained and unconstrained mechanical dynamic
systems. If a convenient set of independent generalized coordinates is found, the
problems are considered as unconstrained systems and the equations of motion

reduce to the form of only differential equations such as :
M(q) ii = f(q.<'r) (2.8)

where M(q) is the system matrix and f(q,q) is the generalized force term. It is not
easy(or possible) however to find the independent generalized coordinates in all
cases. Therefore if we retain Cartesian generalized coordinates, the problem leads to

constrained systems.

When constraint equations (1) are present in a system, the coordinates involved
in the system are not independent and the terms in the Lagrange equations include
Lagrange multiplier terms

The Lagrange’s equations of motion with constraints for mechanical systems

can be written by :

m q).
£[3_T]_fl_2 8—1:; (2.9)

where 7t is a unknown Lagrange multiplier and m is the number of independent con-

straints equations.

Equations (2.9) have the form of mixed differential-algebraic equations such as :

[Mona — daft = mm] (2.10)

(D(q,t) = 0

17

2.1.5. Generalized Coordinate Partitioning

The state of a system that consists of n particles and which is subject to m con-
straints that connect some of the 3n Cartesian coordinates is completely specified by
3n—m generalized cordinates. But sometimes it is desirable to replace the Cartesian
coordinates with a minimal set of Lagrangian coordinates. For example, to facilitate
numerical integration with standard methods, one method to reduce a set of Cartesian
coordinates to a minimal generalized coordinate set is to use coordinate partition-
ing[5].

Define the vector of generalized coordinates of body i by q; = [xbybOi ]T and the
composite vector of the whole system generalized coordinates by q = [q{,q§,...,qZ]T .

Assume algebraic constraint equations
<I>(q.t) = rem). <I>2(q.r),....<1>m<q,r> 17 (2.11)

are independent where 3n > m . If the m algebraic constraint equations are indepen-

dent, the Jacobian matrix for Equation (2.11)
3(1)
(bq = [a—qlmx3n

has full row rank. The implicit function theorem says that if 3n—m generalized coor-
dinates v and m generalized coordinates W which satisfy the constraint equation
(2.11) are specified and (I),v is nonsingular, m generalized coordinates w can be
solved as a function of 3n—m generalized coordinates v and t. The 3n—m generalized
coordinates v are called independent generalized coordinates and the remaining m

generalized coordinates w are dependent generalized coordinates.

With this definition, the vector q is partitioned as

q=IWT.VT]T

 

 

18

where w = b(v,t) by implicit function theorem

Differentiation of equation (2.11) with respect to time gives the velocity relation for

m algebraic constraint equations
(qu = —<l>, (2.12)
Since equation (2.11) can be expressed as
<D(W,v,t) = 0
The Jacobian matrix Q, can be partitioned as
<1)q = [ (Dw, <1)v ]
where (1)“, is m x m nonsingular matrix and (I)v is an m X (3n—m) matrix
The velocity relation of Eq. (2.12) can be written in partitioned form as
(wa + on = —(I>, (2.13)
Since the coefficient matrix wa of w in Eq. (2.13) is nonsingular, the dependent
velocity w can be expressed as
iv = -<r>,,-1<r>,& — (aw-1o, (2.14)
which is a function of independent coordinate position and velocity.
Similarily, the acceleration vector 2] is partitioned as
21 = [ WT, VT ]T
Differentiation of equation (2.13) with respect to time gives the acceleration relation
in partitioned form as
(wa + (by? = -’yww — 7,9 — 2<I>mw — 2¢,,9 - (1),,
where

we,» a awe + car

19
Since the (1),, is nonsingular, the dependent acceleration at can be expressed as
66 = —<r>,,-1[ cm + ywv'v + 7,6 + 2<I>,,,€v + 2am + on] (2.15)

which may be reduced in terms of independent coordinates v,v,'v after substituting

equation (2.14)

Equation (2.10) may be partitioned as

Mwwaw + MW — (DJ). = mm,» (2.16a)
Mvwea + MWV — oft = n(v,v,z) (2.16b)

where

MWW MWV
M = [MW MVV ]
Since (1),, is nonsingular , Eq. (2.16a) may be solved for it such as
it = [owT1‘1[ M‘Ww + MWVV - f"(v,v,t) 1 (2.17)
After all It may be expressed as a function of the independent coordinates v.93 and

input forces.

Using the equations (2.15), (2.16b) and (2.17) we can obtain a set of differential
equations in only the independent generalized coordinates v. Since the a: of Eq.

(2.15) may be written in terms of v,v,°v, we can express :
a = -<1>,,-1<1>,v — (aw-1n (2.18)
where
n(v,v,t) a 7th + 7,9 + 24>,» + 2(1),,v + (1),,
Substituting Eq. (2.17) and Eq. (2.18) into the Eq. (2.16b), we may reduce from the

system with an non-minimal set of coordinates to it with an minimal set of coordi-

nates as follow :

 

J

20
(MW + SMWWST — MWSTW = (Mvw — SMW)<1>:,1n + 1*(v,t,t) - wa(v,v,t) (2.19)

where

MM)“ = s

2.2. Optimal Control General

Optimal control has long been of concern to the applied scientists and engineers
and a broad, general theory based on variational techniques resulted. But earlier
advances in the study of important variational problems was severely limited by the
immense amount of computation required to obtain numerical solutions. The emer-
gence of the high speed digital computer and its appreciation by engineers as a tool
for carrying out complex numerical integrations provided an incentive to continue the

study of nonlinear systems and variational techniques.

Excellent bibliographies on both theoretical and practical aspects of the optimal
control problem can be referred in the survey papers [42,48] and in several text-
books[23,43]. The reference [38] will give good help for studying the calculus of

variation method.

2.2.1. Typical Optimal Control Problems

A typical optimal control problem begins with the following statements :
(i) State Differential Equations

99) = a(z(t),u(t),t) (2.20)

(ii)

21

where z(t) is an state vector with n dimension and u(t) is an control vector with

m dimension and a is a vector function and t is time.
Initial Conditions

The initial state at given initial time to is denoted by z(t0) and the final state at

final time If is denoted by z(tf). The final time tfmay be specified or not.

(iii) Control Constraints

If there are constraints in u, the control history 11 which satisfies the control con-

straints is called an admissible control which is denoted by u e U

(iv) Targets(or final state constraints)

(V)

The final state z(tf) may be required to lie in some set S c R”. S could be a
point, moving point, curve, plane, moving set, etc.
Performance Index ( or cost functional )

‘I

J(U) = h<z<rf).zf) + J gee). um.» dr
1‘9

By h and g, the optimal control problems can be classified as the following

categories. :

0 minimum time problem when h = 0, g = 1, tffree, target given

In
0 minimum fuel problem when h = 0, g = fuel rate, often g = Earl-lull with If fixed,

I=l

or free.
. tracking problem when h = u 20) — r(tf) ”2 g = u z(:) — r(t) ”2 where r(t) =
desired state vector for IDS: Stf

The optimal control problem is to find an admissible control u*(t) on [10, if]

which drives ( with (i), (ii), (iii) and (iv) satisfied ) the system from z(t0) to z(tf) and

among all such controls minimizes J :

 

 

22

that is,

J(u*) s J(u)

where u*(t) = optimal control, to S t Stf

z*(t) = optimal trajectory
f = J(u*) = optimal cost

To solve the above optimal control problems, there are some issues which we

have to consider such that :

a) Does any control it exist for which the system hits the target(means that (i), (ii),

b

V

(iii) and (iv) are satisfied ) ? Here the controllability of (i) is important. The If
"free or fixed" is important for existence. Function a and sets U, S affect this
question but h,g,r does not. Stability of (i) or of equilibrium solutions of (i) is
not so important for the existence question, especially when fixed If case.

Properly posed ploblems would seem to have an infinite number of controls u(t)
for which (i),(ii),(iii) and (iv) can be satisfied. The issue then becomes finding
an ( means not necessarily unique ) optimal control u*(t). This leads the compu-
tational problem. Of course, we also wish to compute z‘(t) and 1". Here it is

noticed that 1" could be unique even when z*(t) and u*(t) are not unique.

c) Unfortunately, for the general problems with n and m bigger than 1 or 2, the

computation problem is so difficult. As we know [24] , for some linear time-
invariant problems of low order, we can compute explicit formula for an optimal
control and the resulting optimal trajectory and performance index. Sometimes
u” can be expressed in the form of a feedback control law (i.e. u*(z*(t) ) which
is a function of the current state by constructing switching surfaces in the state

space even though it is very difficult if n > 2.

23

d) Often, we must be content with expressing necessary conditions for an optimal
control, which may tell us the forms of bang—bang in minimum-time problems

and of bang-off-bang in minimum-fuel problem of u*

Considering the above issues and using variational approach, the necessary con-

ditions of the typical optimal control problems can be derived as following :

The Hamiltonian, defined as :
H(z,u,p,t) = g(z,u,t) + pT(t) a(z,u,t)
where p(t) is the adjoint (or costate) vector with n dimension.

0 State Equations

2(t) =

8H(z,u,p,t) _
_—8p —- a(z,u,t)

0 Adjoint Equations

. _ BH(z,u,p,t) _ 8g(z,u,t) _ Taa(z,u,t)
p(t) - Oz _ 82 p Oz
0 Stationarity

_ 8H(z,u,p,t) _ 8g(z,u,t) Taa(z,u,t)
0 _ au _ Bu + p Bu

0 Transversality for free final state

3h(z(t ),t)
N9) = -—3zf-—[-

Here, it is noticed that many mechanical systems are often represented in the second
order form ( Eq. (2.8) or Eq. (2.10) ). The equations with the second order form

might be reduced to the equations with the first order form such that :

Mia) = fi(z,u,t)

24

where z = [ qT, if ]T

In order to obtain the explicit form of the differential equation such as Eq. (2.20), if

matrix M is nonsingular, the matrix M should inverted such that :
in) = M‘la(z,u,:) (2.21)

Since the standard state equation (2.21) includes a matrix inverse, the necessary con-

ditions for the standard optimal control become :

0 State Equations

 

z(t) = W = M'la(z,u,t) (2.22)
o Adjoint Equations
139) = 3mm”) = -—3———3 ‘2‘”) — thiIW‘atzmm (2.23)
Oz dz dz
0 Stationarity
0 = w = M + pT{-a—[M'la(z,u,t)]} (2.24)
an Bu Bu

0 Transversality for free final state

ah I ,t

Oz (2.25)

The above equations (2.22), (2.23) and (2.24) not only require matrix inversions and
products for each value of 2, but also it involves computing the J acobian of the
matrix-vector products, which is extremely difficult to program and inefficient to

compute.

2.2.2. Optimal Control on Descriptor System

 

25
Consider linear time-invariant systems of the form
Ez(t) = A2 + Bu (2.26)

where E is nxn matrix, which may be singular, and B is nxm matrix, with a specified
initial condition z(t0)=zo.
when E is singular in Eq. (2.26), the system is called a generalized state-space

system(or singular system,or descriptor system [26,27]).

Note that a state variable system such as Eq. (2.20) may be viewed as a special case
of Eq. (2.26) with E=I.

The descriptor theory employs Drazin inverses and Matrix pencil theory [26,27].
For example, solvability rules for descriptor equations depend on the invertability of
the matrix pencil. The solution of a linear time invariant set of descriptor equations
can be stated in terms of matrix exponential involving the Drazin inverse of the
matrix pencil. In equation (2.26), The index of a descriptor system is based on the
order of nilpotency of a matrix related to the system. Solvability results for systems
of index greater than 1 are not straightforward in the time varying linear case and in
the nonlinear case. Gear [12] proved that the constrained Lagrangian equation set is
index three or greater ( greater if the algebraic constraints become dependent). These
equations are notorious for their numerical behavior. For example, incorrect fast
transients may arise after the integration step size changes [9]. In addition descriptor
equations in general are sensitive to initial and/or terminal conditions and exhibit fast

initial and/or final transients[Sincovec].

Now, consider the multi-body dynamic systems which are often represented by
equations involving non-minimal coordinate sets such as Eq. (2.10)

that is,

26

Eq. (2.10) can be reduced to the first order form of Eq. (2.27):
Ez(t) = a(z,u,t) (2.27)
where z = [qT, if, 1T, ATV
Obviously, the matrix E is singular since Lagrange Multiplier It has no accera-
tion term in Eq. (2.10). Therefore, dynamic state equation (2.27) can be regarded as
a descriptor system. Here, there is no information for Lagrange Multiplier at initial
time to (i.e. 100)) because Mto) depends on known q(t0) and £100) and unknown
u(t0). Cobb [27] considered the case that some initial conditions on the descriptor
systems are not specified. All previous works for descriptor system treated the sys-
tem as only first order form. When we attack the descriptor systems by the conven-
tional optimal control method which is discussed in previous section, to get the expli-
cit form of state equations such as Eq. (2.20), matrix E should be inverted. The step
causes the trouble because E is singular. In order to solve the problem, we can rely
on the method of differential equation constraints[24]. The method not only forms
very complicated mixed—differential-algebraic equations, but also its matrix size is
increased because of conversion from second order to first order form and because of
control vectors are added to the state vectors. That makes the automatic generation of

the adjoint equations using the dynamic state equations too difficult.

2.2.3. Mathematical Concepts for Optimal Control

In this section important definitions and theorems for developing the necessary
conditions for optimal control are introduced. The reader is referred to the original
sources [35,37] for detailed inforrnations such as definitions and proofs of theorems.

Variable quantities called functionals perform an important role in many prob-

lems arising in analysis, mechanics, geometry, etc. We define a functional is a rule of

 

 

27

correspondence which assigns a definite (real) number to each function belonging to
some class. Thus, one might say that a functional is a function of function, where the
independent variable is itself a function. The concept of continuity plays an important
role for functionals, just as it does for ordinary functions considered in classical
analysis. In order to formulate this concept for functionals, we have to present a con-
cept of closeness for elements in a functional space whose elements are functions.
This can be conveniently done by introducing the concept of norm of a function.
Reference [39] shall be useful to obtain the basic mathematical concepts such as vec-

tors and normed linear spaces.
Definition 2.2-l : Let a functional J be a map from a normed linear space Q into a
normed linear space E1 which is of dimension one. The function J(q) is said to be
continuous at the point qo e Q if for every 8 > 0, there is a 5 > 0 such that
I J(q) — J(qo) I < 8 provided that H q — qo H < 5

In following lemma 2.2-2, 2.2-3 and 2.2-4, the normed linear space D is defined

on closed interval [to, if] and is continuous and has continuous partial derivatives.

Lemma 2.2-2 : If 0t(t) is continuous in [10, tf], and if
’r

j 0t(t)8q(t) dt = 0

to
for every function 8q(t) e D[t0, If] such that 8q(t0) = 8q(tf) = 0, then a(t) = O for all t
in [[0, If]
Lemma 2.2-3 : If (1(1) is continuous in [10, tf], and if

’r

j «064(2) dt = 0

to

28

for every function 8q(t) e D[t0, tf] such that 8q(t0) = 8q(tf) = 0, then 0t(t) = c for all t

in [to, II] where c is a constant.
Lemma 2.2-4 : If or(t) and B(t) are continuous in [t0, :1], and if

‘1

Jr amaqm + Bear 1 dt = o

‘o

for every function 8q(t) e D[t0, If] such that 8q(t0) = 8q(tf) = 0, then [5(t) is
differentiable and [3(t) = a(t) for all t in [to, :1]. Lemma 2.2-2 is called the fundamen-
tal lemma of the calculus of variations. Proofs of those lemmas are shown in Ref.
[37]

We now introduce the concept of the differential on normed linear spaces. There
are two important definitions of differentiability. One is the Gateaux differential
(weak differential) and the other is the Fréchet differential (strong differential). The
Gateaux differential is defined on a map from a vector space X into a normed linear
space Q. Since the definition of Gateaux differential does not require any norm on
X, the existence of the Gateawc differential is a rather weak requirement and the pro-
perties of the differential are not easily related to continuity. The definition of
Fréchet differential is more satisfactory when X is normed. The definition of Fréchet
differential, which shall be discussed here, generalizes the definition of
differentiability at a point of a real-valued function of one real variable from ele-
mentary calculus. The relation between the Gateaux differential and the F réchet

differential can be found in [37,40].

Definition 2.2-5 : Let Q, W be two normed linear spaces and T: Q —-> W a map,
q, h 6 Q. If there is a d1‘(q,h) which is linear and continuous with respect to h such
that

lim ll T(q + h) — T(_q) -— dT(q, h) II = 0
Ilhll—)0 Ilhll

7

 

 

29

then T is said to be F réchet differentiable at q and dT(q,h) is said to be the F réchet

differential of T at q with increment h
Another way of expressing the existence of the dI‘(q,h) is to say :
there exists a linear map L : Q —> W such that

T(q + h) — T(q) = Lh + n(q,h)||h||

where n(q,h) e W, lln(q,h)ll -—>0 for Ilhll —) 0. Then Lh is called the Fréchet
diflerential of T at q with increment h. It is denoted by dT(q,h). In general L
depends on q, by definition,

dr<q. h) = Lh = T'tqm
The linear map T(q) can be regarded as a function of q and is called Fréchet deriva-

tive of T at q. It is not difficult to show that if such an Fréchet differential exists then

it is unique.

Now we shall apply the concept of the Fréchet differential to the problems of

calculus variations.
Definition 2.2-6 : Let a functional J(q) be defined from any normed linear space Q
into an one dimensional normed linear space E1. Let

M(q.5q) = M + 5C1) — M)

be the increment of the functional J, corresponding to the increment Sq. If there

exists a linear map 1’ (q) such that
Al<q.8q) = J’<q>5q + n<q,8q)IISqII
where n(q,8q) —> O for Il8q|| —> 0,

then the functional J(q) is said to be Fréchet differentiable of J(q) at q with incre-
ment Sq. The linear functional J'(q)8q is denoted by 5.](q,8q). It is called the

Fréchet difierential (or variation) of J(q) at q with increment Sq. The linear map

30

J'(q) is called the Fréchet derivative of the functional J at q. Sometimes J, is called

the gradient of J with respect to q and is denoted by %(or Jq).

Definition 2.2-7 : Let J be a real-valued functional defined on a subset Q of a
normed space Q. The functional J(q) has a relative(or local) minimum at q“ if there
is an E > 0 such that J(q') S J(q) for Ilq — q‘ll < a, any q e 9

On the other hand, the functional J(q) has a absolute(or global) minimum at q‘
if J(q‘) S J(q) for any q e Q. In above definition, if “J(q‘) S J(q)” is replaced by
“J(q') < J(q)” and “any q e Q” is replaced by “any q e Q, q at q'”, the relative
minimum is a strict relative minimum and the global minimum is a strict global
minimum respectively. The maximum J can be defined by minimum of negative J.
Maximum and minimum are called extremum.
Theorem 2.2-8 : A necessary condition for the Fréchet differentiable functional J(q)

It . . . . . III .
to have an extremum at q rs that us varratron SJ must vanrsh on q : that 18,

SJ(q*, Sq) = 0 for all admissible Sq

Definition 2.2-9 : Let T be a continuously Fréchet difierentiable transformation from
an open set D in a Banach space Q into a Banach space W. If (10 in D is such that
T(qo) maps Q onto W, the point qo is said to be a regular point of the transformation
T.

Theorem 2.2-10 : If the continuously Fréchet differentiable functional J has a local
extremum under the constraint G(q) = 0 at the regular point q", then there exists an

adjoint vector p which satisfies such that :
J'(q*> + p"G'(q*> = 0
Definition 2.2-ll : Let E, A be a linear time invariant n><n matrix in complex Cnx"

and to be a real scalar, the n—dimensional vector 20 in C" is said to be a consistent

initial vector [26,27] passing through 2000) associated with Eq. (2.28) if Eq. (2.28)

31

possesses at least one solution. The Eq. (2.28) is said to be solvable(or tractable)
[26,27] at the point to if the initial value problem (2.28) has a unique solution for

each consistent initial vector zo.
Ei = Az + f, z(t0) = 20 (2.28)

Theorem 2.2-12 : Let E, A be a linear time invariant an matrix in complex 0“”.
The Eq. (2.29) is said to be solvable [26,27] if and only if there is a scalar s such

that (sE — A)‘1 exists.
Ez = A2 (2.29)

We note that for time varying descriptor systems and nonlinear descriptor sys—

tems characterizations of solvability are not so easy [26,27] as in Theorem 2.2-12.

CHAPTER 3

THE SECOND ORDER IMPLICIT FORM NECESSARY CONDITIONS
ON OPTIMAL CONTROL PROBLEM
IN MULTI-BODY MECHANICAL SYSTEMS

3.1. Introduction

The purpose of this chapter is to present a set of necessary conditions with
second order form for the solution of the optimal control problem using the calculus
variation approach. The necessary conditions are obtained by posing the optimal con-
trol problem directly on a second order descriptor system. Especially the performance
measure to be minimized may include system descriptor variables such as dynamic
Lagrange multipliers. After discussing the regularity of the descriptor system, the

necessary conditions are derived.

3.2. The Second Order Implicit Form Necessary Conditions on Optimal Control

in Descriptor Systems

Now consider a performance measure of the form :
‘r

J(ua» = h(‘l(’f):‘1(‘f)Jf) +j { gtqtt).c’1<r).x<r).u<r).r) } d: (3.1)
to

subject to the dynamic system constraint of Eq. (3.2) :

M(qfii — daft = f<q,d.u>

<I>(q) = 0 (3.2)

32

33

where M(q) is a n x n matrix, the vector q(t) is an n dimensional state vector and
u(t) is an k dimensional control vector. (1) is m dimensional independent constraint
function(m < n). And the A. with an m dimension is well known as the dynamic
Lagrange multipliers.

Here equation (3.2) may be represented in general second order implicit form of

equation (3.3) :
G(ii(t),ti(t).q(t).l(t).u(t),t) = 0 (3.3)

Thus G is a mapping of Q x A x U into Q X A where q(t) e Q, Mr) E A, u(t) e U.
The symbol h denotes a penalty function designed to ensure arrival of the system at
the desired final position. We note that the performance measure depends explicitly
on the Lagrange multipliers set as well as the state and controls.

The optimal control problem is then that of finding functions qt, ('1‘, 1‘, u‘
minimizing J while satisfying the constraint Eq. (3.3). In order to develop the neces-
sary conditions for optimal control with fixed to, q(t0), (100) and free tf, q(tf),q(tf) the

following assumptions are initiated.

Assumptions
i) For any given control u(t) e U and an initial state q(t0),q(t0), the dynamic con-
straint equation (3.3) has a unique solution q(t), t e [t0,tf].

ii) The functions h, g and G have continuous partial derivatives with respect to

their arguments.
iii) The constraint G has a regular point.

iv) J and G are continuously Fréchet differentiable and J has a local extremum at

the regular point q‘ under G = 0.

34

v) The system is controllable.

vi) Admissible state and control are not bounded.

To show the regularity (Assumption iii) ) of the constraint set G, it is necessary

to show that Fréchet derivative of G maps Q X A X U onto Q X A. that is,

G’:QxAxU0L’0QxA

In other words if the linear equation SG = c(t) where c 6 Q X A has at least one
solution for each c(t) , it is onto. Here SC is Fréchet differential of G. If we take
Fréchet differential of G at some point qo, lo, uo and reduce it to first order form,

then the following linear equation (3.4) can be obtained :
ESi(t) = ASz(t) + BSu(t) (3.4)
where z = [qT,qT,kT]T.
Thus if equation
ESz(t) — ASz(t) — BSu(t) = c(t) (3.5)

is solvable for each C(t), the regularity condition is satisfied. It is easy to show that
(3.5) is solvable via the coordinate partitioning method of section 2.1.5 if the con-

straints are independent.

Assumption iv) is needed to apply the standard Lagrange multiplier theorem

[2.2-10] of constrained optimization theory.

Controllability of descriptor systems is described in Campbell[4], where control-
lability is defined in terms of the equation’s R-reachable set. Briefly, a system is R-
controllable if it can be driven from a specified point in the reachable set to another

specified point in the reachable set. Here we can interpret the R-reachable set as

35

equivalent to the Lagrangian configuration space(C-space) associated with system G.

Now to develop the second order implicit form of necessary conditions for
minimizing J of Eq. (3.1) subject to system equation (3.3) for fixed to, q(t0),q(to) and
free If, q(tf), (1(9). we follow a standard variational approach[24].

Since h is a differentiable function, we can write

‘I

h<q(r,),<1<tf).tf) = fl §;[h(q(t).£1<r>.r) 1 } dr + htqoxioro ) (3.6)
‘0

so that the performance measure becomes :

’/

J(um) = j{ g<q<r>.<1(z).1<r).u<r>,t) + $144611) 1 } at + h<qe<1eto>
’0

Since ('10 ,qo , to are assumed to be fixed, the minimization does not affect the

h(qo,q0,t0) term, so we need consider only the functional

’/

M9» = j{g<q(t).d(t),t(r).u(t).9 + %[h(q(t).<i(t).t) 1} dt (3.7)
to

using the chain rule of differentiation Eq. (3.7) becomes :

‘I

T
J(u(r)) = j{ g(q(t).<°1(t).l(t).u(t).t) + [g—Z(q(t).<i(t).t)l ('10)
‘o

T
+ 18—i(q<t),<‘1<r>.t>1 tic) + §£(q(,),¢(,),,) } dr
aq at

Due to assumption iii) it is possible to apply the Lagrange multiplier theorem to our
system. Thus by the multiplier theorem [2.2-10] there exists an n+m dimensional

vector-valued function p(t) which is called the adjoint vector. Therefore the

 

 

 

36

augmented functional is defined by ;

‘I
J.(u(9) = j{ g.(ii(9,<'1(9.q(9.u(9,x(9.p(9.9 } d: (3.8)
to

where

T
g,,(ti(9,£1(:),q(:),u(z),x(:),p(9,9 =--' g(q(t),(°|(t),2.(t),u(t),t) + [gs-(«minim ('10)
+ {—(q(t).q(t).t)]T ii(t) + —(q(t).q(t) t)
+ pT(9G(ii(9.<1(9,q(9.u(9.x(9.9 (3.9)

Suppose that u‘(t) is an optimal control and that q'(t) is the corresponding optimal
trajectory. The increment of the functional Ja corresponding to the increment Su at u*

is defined by
Ala(u*,Su) = Ja(u‘ + 5n) - 1001*)

Since the J is assumed to be Fréchet differentiable and the Fréchet differential Sla(or

variation) of Ja(u’) is a linear functional,

For free If and free final state, the SI“ at u' is given by :

 

‘/
at a a at .1: «t a: a: :1:
Wu > = ]{ t 3‘: (q (94 (9.1 (9.7» (9.u (9p (991T 2399)
to
880 t 0* 00* t * * T 9
+ [SE-(q (On (On (9.1 (9,1. (9.1) (99] Wt)

*(9,<°1*(9.ii‘(9.x*(9.u*(9.p"(9.917821(9

       

+:[g ,f*(q (9419,4191(9u**(9.p (9 91T 8M9

37

a 0 U 00‘ t t *
+ [£44194 (9.4 (9.1(9.u (9.1) (9.9]T 8u(9

age

ap(4'(9.('1"'(9.ii"(9.2t“(9.u‘(9.p"(9.9]7‘ 8p(9 } dt

+[

 

+ [ga(q‘(tf).ci’(If).ii'(tf).l‘(tf).u‘(tf).p'(tf).tf)l69 (3.10)

where

849=u9—4m
599), Sii(t), 8M1) and Sp(t) are defined like similarly. The superscript T denotes tran—
spose of a matrix.

Integrate by parts :

‘I

a a e .t “at g. c a: .
jt «fin (9.4 (9.4 (9.9» (9.u (9p (9.91T84(9 d:

to
aga t e t «t y t y T
= 1 9—9“ (9.4 (9.4 (9).?» (9.11 (9.1: (9.91 64(9)

38..
84

 

- [ (4‘(to).4‘(:o).4‘(taxeoiu‘ootp‘(to).zo) 1751.90)

 

’1
a t 00* j y t
— ]{ i g." (4*(9.4 (9.4 (9.1(99 (9.1: (9.9 }T 849) d: (3.11)
to dt 8‘1

and

 

‘/

a a n at «t g g a: u
j [ 3‘2 (4 (9.4 (9.4 (9.19)» (9.p (9.91" 849) d:
to

a a n c It on t 41 a: g: o
= [ [Biz-(4 (9.4 (9.4 (9.19)» (99 (991T 2349) 1:2,

38

             

dd, 9: (9.4 (9.4 (91 (9. u *‘(9,p (x), 917 Sq(t) 1:,

‘r
d2 g
+£“dfia

Since q(t0), 490) are specified, Sq(t0)=Sq(to)=0, in equation (3.11) and (3.12) terms

*(9.4“(9.4'(9.x*(9,u’(9.p*(9.9 1T 549) } dr (3.12)

     

which include q(to), and 490) vanish. If we substitute the results of Eq. (3.11) and

Eq. (3.12) with 6q(t0)=64(:0)=0 into Eq. (3.10), the 41,,(u‘) becomes :

38a

al.(u*) = { [ a4 (4‘(:,).4“(:,).419).x’(9).u‘(9).p‘(9).rf)1

 

d 380 t 0‘ 00* t t . T
— [1374“ (9.4 (9.4 (9).). (94 (9)9 (9.9) 1} 84(9)
380 t 0* 00* t t t T 0
+ { a—am ((9.4 (9M (9).?» (tf).u (99.1) ((9.9)) 54(9)

+ {g.(4’(9).419.479.119.1149.11199) 159

‘1
380 . .9 «t a: t t
+ 9 [37m (t),q (9.4 (9.1(9.u (9p (99

     

 

   

a

9g
+ {—(4' (9.4 (9.4" (9 x (9 u *"‘(9.p (9 9]" 5M9

aga a .4 «t g .1 g. T
+ [Fa-(4 (9.4 (9.4 (9.7» (9.u (9.1: (9.9] 8u(9

39

a a t J «t t t r t
+ ['51—«1 (0.4 (0.4 (0.1 (0.11 (0.!) (OJ)? 5P0) } d! (3.13)

From Figure (3.1) we define that :
8Qf 5- 84(9) + 4’(tf)8tf 84fE 84(9) + $0,959

Since 8q(tf) depends on Sq, and 8:}; we write

 

 

 

 

 

 

 

 

 

84(9) 5 84,— 4“(r,)5r, 84(9) 2 84,— 4‘(rf>5r, (3.14)
Q“) A A A
4(tf) " \
CI“) 1 Squ)
a:
400) ‘1 (9
to t; W?

Figure 3.1 Relation between Sqf and 5q(tf)

By the Theorem [2.2-8], 8.10 must vanish at u‘ for the Fréchet differentiable func-

tional J(q) to have an extremum.

 

40

Thus substitute 8q(tf) and 831(9) of Eq. (3.13) into Eq. (3.14), and 8/0 becomes :

alam“) = 0
aga a: .t . "a. 4. 4. 4.
= { [871(4 (9.4 (9.4 (9.9. (9).u (9).p (9.91
d ago 4 .t ”t g 9 .9 T
— [Ea—am (9.4 (9.4 (9.1 (9.4 (9).p (9.91} 84;
ago 4. .t “at: . .9 . T .
+ { 321"“ (9.4 (9.4 (9.1(9.u (94» (99} (4
+ { g..(4*(9.4"(9).4’(9.1*(9).u*(9.p*(9.9
aga 0 ct «t a. a; t
— [371m (9.4 (9.4 (9).1(9).u (9).p (9.9)
d 88 .t at x: at t T .9:
— Ea—(q‘ (9.4 (9.4 (9), (9).u(9).p(9.91 4 (9)

— t—(4‘ (9) 4 (94 ".(9 a: (9) u (9). p “.(9) 9? 4(9) }89

‘1
a a .9 . at: «It a: It It
+ {fl {—a‘i—(q (9.4 (9.4 (9.9 mm (9.9 (99

     

_d8
-d-; a: ".(94 (94 (9X (9a (94)" (9 9
7d8 ea in ”(9.4 (9.4“ (91" (9.4 """(9.p (9.91" 84(9

aga t .1: «t a: a: at T
+ [Ex—(4 (9.4 (9.4 (9.1(9.u (9P (99] 8M9

8 a t . t “at g .9 *
+ [-%(4 (9.4 (9.4 (9.10)» (9.p (9.9]T 6u(9

41

a a t 0 t 00. t t t
+ [7:44 (9.4 (9.4 (91(9): (9:) (9.91T 8p<9 1 d: (3.15)

From Eq. (3.8)

g..(4(9.4(9.4(9.1(9.u(9.p(9.9 E g(4(9.4(9.1(9.u(9.9 + l-gfi-(q(9.4(9.91T4(9
+ [—.(4(t).4(t). t)1T 40) + —(q(t).4(t) t)
+ pT(9G(4(9.4(9.4(9.M9.u(9.9

Rewrite the coefficients of Sqf, 84):, 8t], 8q(t), 8u(t), 51,0), 5p(t) respectively:

To get the coefficient of 8q(t)

     

# out a. g g g 0 t «It 2 at 0 It 0 It
‘(9.4 (9.4 (9.1(94! (9.p (9.9= gfi-(qaz (9.4 (9.4 (9.9+ [Sq—’34 (9.4 (9.914 (9

82h
aqa

 

+[ (1"(4 (9.4 (9. 914 (9 +8 —%(4 (9.4 (9. 9
+ [p 137%) G(4* (9.4 (9.4 (9 7L (9 u (9. 91T (3.16)

(9.4 (9.4 (9 X (9. u (9 p (.991

       

ddtagf:

=d El —é‘-(4* (9.4 *.".*(9x(9u (99H;ll

 

:qu (9.4 (9.914 (9
2
+ %E(4‘(t).4t(t).t) + [%(q'(9,4'(9,9]4*(t)

+ aa——(4' (9.4 (99+ [Him—$4 (9.4 (9.4 (97» (t)u (091” (3.17)

=7}%@ (9.4 (t), 9. (t), u “,(9 t) + —t[——(4* (9.4" (9. 914 *(9)

32h
+—[——(4 (9.4 (9.91+ —[[-.—-(4 (9.4 (9.914 (9}

42

2 1 0 t t a i ‘0‘ $ t
+ 1(fi(4‘(9.4 (9.9) + [p Tog—(gm (9.4 (9.4 (9.1(9.u (9.9]T

 

dt aaa:
+ tWe);13—9(4‘(9.4‘(9.4‘(9.x‘(9.u‘(9.91 1’ (3.18)
1‘ aq
‘12 [ag,( *(9 “(9 ”(91*(9 *(9 *(991
dtz Ba q ’q ’q ’ ’u ’p ’
2 2
= fili—ZWYOhYOJfl + ff;(fit)?(4‘(9.4‘(9.4*(9.x'(r).u*(9.9]T
Since
d2 ah d d ah d 82h . a2}... a2}:
__ = _ __..—_. = — ———— -— — 3.19
a? aq dt d: 3;, atlaqaq“ 3:12“ aaat] ( )

d2 Tac -11 Tac

(122“) 34] dt[dt(p 74)]
_£_ .TaG 713G
—dt[p 321 +9 dz 3&1

..Tac .1 d ac .Td ac 9.12 ac
= __ + ___.__ + __ + _.__.
p a4 p(drabi) pdzaii (1:2 a4

«Tao .Td ac 9.12 ac
= _ +2 __ + __
" aa " dz a4 " .112 aa

therefore

 

''2 800 q ]

_1a2h... .. 19:9,...
— (1,118 a. (4 (9.4 (9.914 (9}+ dtl [aqz (4 (9.4 (9.914 (9}

2 I! "I! y .11: "It * 1:
+ (jg—3.” (4‘(9.4 (9.9} + { p «fa—9m (9.4 (9.4 (9.1(9.u (9.9)T
: a a 9.;

daG

+ 2( 449737(4‘(9.4‘(9.4*(z).1*(9.u*(9.9 )T

43

”1 [31407289701 (1)"! (09d (1)9x.(t)’u*(t)al) }T

and from Eq. (3.18)

2
113-44094 (9 9)- - [—(4 (9.4 (9. 914‘ (9 +1aq a h

 

q'(4 (9.4' (9. 914 (9

82h t 0 ‘
— t , ,t
a at(4 ()4 (9)
After cancelling out all common terms, finally the coefficient of 8q(t) becomes :

13%(4‘(9.4‘(9.x‘(9.u‘(9.9 1T+‘ p Tag—(9 G(q‘ (9.4 ‘(9.4 “.‘(9x (9.”"u (9 9

-[——g-(4 (9.4 '(9.X'.'(9u (9 9 1T [13 ”(938— (4 (9.4 (94 '.'.'(971(9u (9 9 1T

- { p' Tt()-g- :6 (4 (94' (9.4 '.'(91 (9.u' (9 9 + 4' T(-—9 :(4 ".(94' (9.4 '.'(97» am" (9. 9
+ 2:3 ”(9112(4‘ (9.4 (9. 4 ".(91 (9 u *(9 9
+ P TI()— 3:: 839 (4 (0.4 (0.4 (1).?» (t) 11".(t) t) } (3.20)

If we rewrite the Eq. (3.20), it becomes :

[-—.(4 (9.4‘ (94 (9349449914549

+{12 1%.‘3-(4 (9.4‘ (9.4‘ (91. (9 u ‘.(9 91 — ——.-(4‘ (9.4‘ (9.4 ‘(9.x*(9.u*(9.9}Tb‘(9

+{ [Era—(4' (9.4' (9.4 '.(9 71 (9 u '(9. 91 - [1339(4' (94' (9.4 49.1 (9. u '(9. 91

+ '44(4‘(:).4‘(9.4‘(9.x‘(9.u‘(9.9}Tp‘(9

+1%:-(4'(9.4'(t).1'(9.u'(t).t)1 ~1—3fi-(4 (9.4' (t), 1* (t) u ‘(r), 1)] (3.21)

The coefficient of 8u(t) is

a a :tr -* n! a: t It
341 (4 (9.4 (9.4 (9.1 (9.u (9p (99

 

= 3%(4‘(9.429.149.4799+ [13—3(4‘(9.4‘(9.4‘(9.x‘(9.u‘(9.91Tp‘(9 (3.22)

The coefficient of 81 becomes :

a a II t 1' 00* t t 1
3i (4 (9.4 (9.4 (9.1(9.u (9p (99

 

= 3%(4‘(9,4‘(9.x‘(9.u‘(99 + 118%(4‘(9.4‘(94(9).“(9.u‘(9.91Tp*(9 (3.23)

There are still the terms outside the integral in Eq. (3.15).
Since

a 2 2 2
1% =§1+fl+—ah.4+fl"+—a.h MEG—V
aq Bq 34 aqaq 8612 aqa: aq

a
li _ 1% + i[pTa_G]T

dt 34 dt 35, dz 34

_a2h. a_211_.. 82h

 

aqai;q aaz aaa: [a4] p [dz a4] p

The coefficient of Sqf becomes :

11(4‘(9).4‘(9),2*(9.u*(9).9) + 1979,4799
aq aq

8G 1 4* “a: at: at T at
+ [33-(4 (9.4 (9.4 (9).?» (9).u (9.9)] p (9

— 118%(4‘(9.4‘(9).4‘(9.x’(9).u‘(9.91Ti)‘(9

45

- {-537.44 (99.4" (9.4' (9)1 (9). u (9). 91p p(9 (3.24)

Similarly, the coefficient of 64f becomes :
444494499 + 1%?(4"(9.4‘(9.4‘(9>.1‘(99441191441) (325)

To get the coefficient of Stf , in Eq. (3.15) we can get the following equations.(For

convenience, we drop * and ti)

Using Eq. (3.17) and Eq. (3.19),

 

 

 

 

since
a 2 2 2
—g.—" =§@+——ah.4+-a-'1+M4+——a.h 41139—171)
aq aq aqaq 34 8612 aqat Bq
1330-124 3&2
draii -dt[a4+(a4)p]
= a2”4+1'14Jr—a'h +(a—G—)’b+11(1‘3-)’1p
aqa4 3:2 848: aa d: a4
38.. dag”
{ga-t 277814-271 ..1’4 1
=1g+(-g—2’c’1)+-2—(2)’4+%—h+p’c1
_% a__h_+_a_c’p_%’._ dacT
[ aq (7(1):) ( ) ”(d—77MM
_ .14., T"
'8q+ '§a— (5)54]

Coefficients of Stf will be obtained by :
g(4‘(9.4‘(9).1‘(9.u‘(9.9) + p"(9G(4"(9.4‘(9.4‘(993(9).u*(9.9)

- { 32-2(4‘(9.4‘(9.x‘(9.u‘(9.9) + 132949.449,»

46

3G d 8G a: .4 ..4 a: a: 4
+15 — 27-53-144 (4),4 (4)4 (4)31 (4).u (4),4)p (4)

— 12—:(4”(4).4"(4).4“(4w(4).u’(4),4)1’i)*(4) }T 4*(4)

— { 1%?-(4’(4),4*(4).4‘(4),x‘(4),u‘(4),4)1’p*(4) }T4‘(4) (3.26)

Since the 8],, must be zero to obtain the necessary conditions for a extremum,
the necessary conditions of the second order implicit form which determine the
optimal control u’ and state vectors q‘, ('1* and 74' can be expressed as follows :

1) State Equations

Coefficients of 5p become zero because the constraint G should be

satisfied.
G(4“(r).4‘(r).4‘(r).1*(r).u*(t).r) = 0 (3.27)

11) Adjoint Equations

For convenience we will drop the * and time t.

ac 4.. (1 ac ac T. d2 ac d ac ac T
a d 3g = 2
+ iaq _ —dt( ad) 0 (3.8a)
and
a ac T :
’g'ax +(—a).) p () (3.28b)

1H) Stationarity

Coefficients of 6u are 0 because 8n is arbitrary

2a 1, [32915, = 0 (3.29)
Bu an

47

IV) Transversality
If the final state is free and final time is fixed, coefficients of Sqf and 564

will be zero and 82f is also zero.

:34 if}. a_G_1_3£T_4_Gr- =0 3.30
{aa+aq+[(aa) (dream (aa)p}‘=" ( a)
And
{—8}? + (_ac)4p] =0 (3.30b)
8q 34 ”I

As a simple example, we can consider the following problem such that :

Minimize the performance measure of the form :
1 x’ 1
J = 3[ 23(4) + 12(4) ] +j 303’ + x2) d:
‘0

subject to Mjc'+ Ci + Kx — u = 0
where M, C, K are constant and x, u are scalars.

The second order form necessary conditions become :

M55 + C5: + Kx - u = 0 (State)
Mfi — Cp' + Kp + x = 0 (Adjoint)
u — p =0 (Stationarity)

x(tf) + Cp(tf) — Mp'(tf) = 0 and x'(tf) + Mp(tf) = 0 (Transversality)

It is shown that the above results finally match to those of conventional approaches
which are the standard state space optimal control[10,43], the method of differential

equation constraints[24] and Wu’s method[32].

CHAPTER 4

NUMERICAL VIEW

4.1. Introduction

The equation of motion arising from rigid body dynamics are described by a set
of differential equation or a set of mixed differential-algebraic which is difficult to
solve analytically. Therefore a numerical integration procedure must be used to solve
the dynamic equations of motion. In addition the optimal control problems by varia-
tional approach necessarily lead to a nonlinear two-point boundary problem which is
extremely diffiult to solve. Since finding solutions to nonlinear two-point boundary
value problems are, in many cases, not necessarily a trivial extension of finding solu—
tions to initial value problem, the difficulty of solving nonlinear two-point boundary
value problems (which is usually beyond analytic tractability) make us resort to

numerical methods [20,21,22] for the solution of such problems.

There are many different types of integration methods, each with its own advan-
tages and disadvantages. Integrators used in mechanism programs are multistep
integrators. These integrators give high accuracy for the amount of numerical work
required. Within the class of multistep integrators there are those that work well on

stiff equations and those that do not.

In this chapter the second section deals with the Gear stiff implicit integration
method. The third section shows how the Jacobian is formed when the integration
method is applied to the system dynamic equations and adjoint equations of the

optimal control.

48

49

4.2. The Gear Stiff Integration Method

Any system of differential equations that has widely seperated eigenvalues, at
least, locally, is called a stiff system. In other words a system is stiff when the actual
solution is changing slowly but there is a solution nearby which could change
rapidly. The possiblity of a rapidly changing solution indicates that the system has an
eigenvalue of large magnitude. The fact that the solution is changing slowly indicates
that this eigenvalue is not affecting the solution. It has been shown in references
[15,16] that neither Runge- Kutta nor Adams-Bashforth nor Adams-Moulton algo-
rithms are suitable for the solutions of such systems, for stability reasons[15,16]. In
1969, Gear developed stiffly stable [15,16,17,l8] multistep algorithms, which are
better suited for the solution of stiff systems. Later he showed [19] that the method
can be used for the DAB although numerical difficulties have been reported for

descriptor systems at high index [12].

The general form of the k-th order Gear implicit algorithm can be given by the

following equation (4.1)

k—I
qn+l = 2 Win—i + hBanH (4~1)

i=0
where h is the time step, qn+1, q", . . . are the values of q at time instants tn+lt tn, . . .
and a0, a1, . . . , (1k_1, [30 are the (k+1) coefficients known as Gear coefficients for

this multistep algorithm[,]. Since the k-th order Gear’s algorithm is an implicit mul-
tistep algorithm, it is necessary to solve an implicit equation (4.1) in each time. The
k-th order Gear’s algorithm for solving equation (4.1) can be recast into the general

form at t = In“, i.e.

2( in”, qua, In“ > = 0 (4.2)

50

Applying the Newton-Raphson algorithm to equation (4.2), we obtain

2g.(M)Aq(m) + 22.0") Ado") : _g(m) (4.3)
3q at]

where
Aqon) = q(m+1) __ q(m)
A510?!) = awn _ (ion)

The superscript m is the iteration number. The time step counter n has been dropped

k-l
for srmplrcrty of notation. Since .2 ajqnfi- remains invariant for Newton s iteration

at the (n+1)st time step, one obtains Eq.(4.4) from Eq. (4.1)

1
Aq<m> = — Aq<'"> (4.4)
Bah
Hence, Eq.(4.3) becomes Eq.(4.5)
22"” ., 4.22"” W) = .530» (4 5)
Bq Bob ad

Equation (4.5) is the Newton-Raphson corrector for implementing Gear’s algo-
rithms. In Eq. (4.5) the coefficient of Aq is called Jacobian. As a consequence of the
Mean Value Theorem, the Jacobian is not usually updated for every iteration to
reduce the computational cost and a sufficiently accurate solution will be attained.
The only penalty for not updating the Jacobian may be a requirement for more itera-
tion steps in the corrector procedure. This is more than compensated for by not hav-

ing to re-invert the Jacobian at each step.

The stiff integration method is a predictor-corrector method which has nice
features such as self-starting, variable time step and variable order for automatically
choosing the optimal time step and order. Besides being stiff stability, the fact the

Gear stiff method works for the mixed-differential-algebraic equations is important

51

for the types of equations that arise in many constraint mechanical dynamic system
and network analysis. One of the most important aspects of this method is the spar-
sity of the Jacobian matrix. Taking advantage of the sparsity of the Jacobian matrix
can save both execution time and storage space while processing of the decomposi-

tion and back substitution efficient1y[4,14].

4.3. Formulation of Corrector for System Dynamic Equations and Adjoint Equa-

tions

The corrector formula for the system dynamic equations will be employed for a

combination of Eq. (2.10) and Eq. (4.2)

i.e. At (n+1) time step Eq. (2.10) can be expressed by the implicit form as :
.. . .. .. ad) T .
f(q,n+1aQn+1aQn+l,xn+lvtn+l) = M(qM1)Qn+l _ [3;]Mler-l — f((ln+1vqn+1) = (14-63)
¢(qn+1) = 0 (46b)

If we combine Eq. (4.6) and Eq. (4.2), Eq. (4.7) can be obtained at the (n+1) time
step :
Retake
G(q,é,ii,l,t) = 261,51) = 0 (4.7)
<D(q,t)

The subscript n+1 has been dropped for convenient notation in Eq. (4.7). Since Eq.
(4.7) can be solved by Newton-Raphson iteration, at the m-th iteration the Eq. (4.7)
becomes Eq. (4.8).
i.e.

. ~ (m) m * (m) ~ (m)
1"") + if- Aq<m> + g!— Ado") + £4, Air“) + % AW) = 0 (4.8a)
q q

( )
gm) + 33 MW + 22 Aqua) _ 0
q
(1)0") + gym Aq‘”) - 0
q

where Aq‘m) = q‘m‘q) — q‘m) and Ado"), Ali“), A10") are defined similarly.

Since
(m )
An<m>_§(i(m Amt), 2i” _ I
q aq q aq (— 7-30)
Eq. (4.8) can be rearranged into Eq. (4.9)
31qu (”0+ [_ + _1__ 3? Jo") A510") + a_i(’")M(m) = _;<m>
aq hBO aq a).
.. (m) (In)
33. A (m) + 22. A‘(m) : _A(M)
a q q 351 q g
acb ("0 ( )
.___ "1:=._qp0n)
aq Aq
Eq. (4.9) can be written in matrix form of Eq. (4.10)
' .. .. . (m)
2: (a. . _1__a_t_ E
aq Bq hBO 3Q a}. Aq (m) f (m)
33 3g 0 Ar. =- g
q q -Ai <r>
fl 0 0
La“

 

 

(4.8b)

(4.8c)

(4.9a)

(4.9b)

(4.9c)

(4.10)

In a similar manner, we can obtain the corrector for the adjoint equations for

optimal control discussed in chapter 3. In order to minimize the reaction forces we

might choose a performance measure such as :

‘1

J = [ g(u(t),i(t) )dt
‘0

(4.11)

53

subject to the system dynamic constraint of Eq. (2.10).

If we apply the second order implicit necessary conditions of Eq. (3.28) to Eq. (2.10
) and Eq. (4.11), the adjoint equations become 2

- .. f . f
f = Mqu + {—8. lqu— [($51)q + a—]qu + <I>£pl (4.12a)
aq Bq
~ 3
(I) = <1>quq — 5i- (4.12b)

where pT = [ qu, pr ]T. pq and p; are the adjoint vectors corresponding to the sys-
tem vector q and 1 respectively. Here, since the adjoint equations are integrated
backwards in time, the negative time step is adopted in the k—th order Gear implicit

algorithm of Eq. (4.1) as :
~ . k—l .
g( pn+l9 Pmi’ tn+l ) = —Pn+1 + 2 aan—i — hB0pn+1 = 0 (412(3)
11:0

By Newton-Raphson, the corrector of adjoint equations becomes in matrix form :

a_? if. _ _1_a_‘f a_‘f (”0
afiq apq {1150 apq 3P1. qu (m) ~f ("0
321 fig 0 M», = _ ~g (4.13)
apq apq ..
~ Apr ‘1’
ac»
— 0 0
apq

It is noted that the relations between each block element of Jacobian matrices for the

dynamical system equations and the adjoint equations can be written as :

_a_i
Bq
a_i
851
ii.
Bii
91
a).

A

_g_
Bq

A

_8_
3f:
E
3q

54

a_'r
apq

_3_~fT=
(8‘3)

)T=

(

q

8f
_(E)_
-_3"g_ =
apq
22.
35,,
£1
apq —

> (4.14)

 

A

Since the elements of these Jacobian matrices are the same except for the matrix

transpose and the sign, the results make their automatic generation easier.

In order to start the integration step forward or backward we need initial or final

values. But there are no information for ii(t0), 1(t0),i5q(tf) and p1(tf). Because the

M10) depend on q(t0), q(t0) and unknown u(t0) they can not be specified in advance.

Similarly the adjoint vectors pl(tf) corresponding to the unknown Lagrange multiplier

are not known. At the starting point a direct integration method[7] can be used to

find the initial ii, A and the final liq, p1. If we simply form the first and second

derivatives of (I) and ED with respect to time and rearrange respectively, Eq. (4.15)

can be obtained :

 

55

T u .
M ‘Dq [q] _ [f(q,q,u,t)]
(Dq 0 mo ‘7‘ t=10 — Will) t=to

(4.15)
t=l/

(bq 0 V(q,q,q,Pq,Pq,1) (:1!

p}. Mf—

where
vth) = —(<I>q<i)q€r - 2<I>q<i
T: — [83—21qu + [(431), + 37:]qu
v<q,4,ii,pq,bq,i> = [<¢,pq>q4]qd + <¢qpq>qii + <I>qi5

. . . . d2 a
+ [((quq)qQ]pqpq + ((quq)qq '- E(5§T)

Since the coefficient matrix of Eq. (4.15) is nonsingular[7], we can find the initial
and final values.

CHAPTER 5

EXAMPLES

In this chapter, we give three examples to illustrate the applicability of the
theory developed in previous chapters. The examples show the second order form of
necessary conditions for the different performance measures and the results of the
computer simmulations. The first example consists of two constrained sliding masses.
The second example is based on a simple pendulum problem. Finally simple two

degree of freedom manipulator model will be followed.

5.1. Two Constrained Sliding Masses Problem

The problem to be analyzed is shown in Figure 5.1. It consists of the two
masses which are each allowed one degree of freedom of translation along the
horizontal(x)-axis. A constraint requiring that the two masses have equal displace-

ments is also imposed. The governing equations in Lagrange form for this system

are :
mpr'l + klxl + )V = u (5.1a)
mzjc'z + k2x2 — 7L = 0 (5.1b)
x1 —x2=0 (5.1c)

where m1, m2 2 two masses respectively.
k1, k2 = stiffness of the springs connected to the respective masses.

7» = Lagrange multiplier corresponding to the equation of constraint.

56

57

 

 

 

 

 

Figure 5.1 Two constrained sliding masses
Case (la) : Minimizing the control effort only

Consider minimizing the performance measure :
1 t’
J = — u2 dt
2 l,

Using the second order form of necessary conditions the adjoint equations become :
mlfil + klpl + p1 = O (5.23.)
Wizfiz + k2p2 — pl = 0 (5.2b)

p1 _ p2 = o (5.20)

58

The Stationarity condition is given by

u = p1 (5.3)

And the transversality condition is :
prof) = 122(9) = 0 (5.4a)
251(9) = p'2(zf) = 0 (5.4b)

The exact solutions of this system can be obtained easily, that is, the optimal control

* . . . . . . .
u (t) is zero at given time interval and the optimal trajectories are :
t t 1 - .
x1(t) = x20) = x100) cos wt + ——x1(t0) srn wt (5.5a)
w

where

_ (k1 + k2)
_ (m1 + m2)

The corresponding Lagrange multiplier will be ;

mk—mk . m *

___12 “um—2 u(t) (55b)
1

m1+m2 m1+m2

W) =

To obtain numerical results the following parameter values are used. They are :
m1=10 k1=20 m2=5 k2=10
’0 = 0 tf=1 x100) =1 x100) = 2
initial guess u0 = 10

Figure (5.2), (5.3), (5.4) and (5.5) show that the numerical results are exactly the

same as the analytical solutions.

 

 

 

 

 

 

0.1 51”““— I l r l I“—‘"“‘!“ I T I T
0.1 0‘1...” ---_m___._---._--..._-.... ..._..._.._2
1 _

l
0.05-1 —«

l
4 4

i
0.004 ------------------------------------------------------------------ _.
J: ———- non—minimal “
4.05%“. em?“ 4 , . r -—.—«-—-—,——-———.—_—4

0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.2 The optimal control for case (1a)

 

 

 

60. T T" “l—' ”r "_—T'“ ~‘T”""‘!“ " 'Th—m’jm' ' 1" l" 'ml ' I T" I" l
o exact
. o 4> non—minimoL
40. —4 ~
4 4
o i
l
20. ”‘1 ..i

 

No. of Iterations

Figure 5.3 The performance measure for case (la)

X1

60

 

 

 

 

2.0—i' 1 l ‘1 l 1 Ti _ r l l
1.5— T
1 0M —
0.5— ~
“ — non—minimal “
e ct
O-OJV'" . X0) . l r F ' l r
0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.4 The optimal trajectory for case (1a)

 

0.10",» . I . , ,_,- To?

0.05~ ..

0.00

 

l —— non—minimal
_0.05_i":'—;r exaft | i r r l
0.0 0.2 0.4 0.6 0.8 1.0

 

 

r r "—4

 

Time ( sec )

Figure 5.5 The Lagrange multiplier for case (1a)

61

Case (1b) : Minimizing the control effort and the reaction force

Consider minimizing the performance measure :

J = (u2 + if) d:

i
2

5“—-~:.‘

The necessary conditions are same to Eq. (5.2a), (5.2b), (5.3) and (5.4). Since the
performance measure includes the control and the Lagrange multiplier, instead of Eq.

(5.2c) the adjoint algebraic equation becomes :
pl-p2+)\.=0 (5.6)

. . . . 1‘1 k2 .
For srmplrcrty, if — = — = wz, then from the system dynamic Eq. (5.1) the
m1 m2

Lagrange multiplier can be expressed by :

M!) = [514(1) (57)
where
_ _’Z’2__
B - m1 + m2

From the Eq. (5.6) we obtain such that :
[320) = (1 + B) 510)
The adjoint equations are reduced to the form of :
[m1 + m;(1 + [3)] 1510) + [k1 + k2(1 + [3)] 1010‘) = 0 (5.8)

By the transversality condition Eq. (5.4), Eq. (5.8) implies that p1(t) is zero in the
given time interval. Therefore, the optimal control u*(t) and the corresponding

Lagrange multi lier WU) are both zero. The 0 timal trajectories are expressed by :
P P

1 . .
xi(t) = x30) = x100) cos wt + —x1(t0) srn wt
w

62

The figure (5.6), (5.7), (5.8) and (5.9) show the numerical results with the

parameters :
m1=10 k1=20 m2=5 k2=10
t0 = 0 tf=1 x100) =1 i100) = 2
initial guess uo = 10
Figures indicate that results of non-minimal approach by second order form necessary

conditions are close enough to exact solutions. We note that a fast terminal transient

has appeared in u(t) and Mt).

 

 

 

 

 

 

0.15 l I l I l I r I I
0.10-
V.
:3 0.05— —
0.00
‘ -— non—minimal ‘
—0 05 exact
I l I l I I I I I T
0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.6 The optimal control for case (lb)

X1

63

 

 

 

 

 

 

 

 

60. . , . I . I . I , I ,
4 o exact
‘ o non—minimal_
40.4 —
O
20.4 _
.1 0 _
O
O
O. T j l .r % Q fi Q 9| {i
1 3 5 7 9 1 1 13
No. of Iterations
Figure 5.7 The performance measure for case (lb)
2.0 ‘l I T l l l 5 l I
1_5_. _
1.01 _
~ 4
0.5— _
‘ — non—minimal ‘
---- t
0.0 r a“? . f 4 , . r f
0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.8 The optimal trajectory for case (lb)

 

 

 

 

 

 

 

0.10 i I l I I I l I l
—I —
0.05— —
,. .1
0.00
— non-minimal
t
—0.05 r exo? 1 I r I F If r
0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.9 The Lagrange multiplier for case (lb)

Figure 5.10 shows the result of following parameters :
m1=0.5 k1=1 m2=0.5 k2: 10
t0 = 0 If: 1 x100) = 1 .2100) = 0.
initial guess u0 = 10
By using symbolic mathematics software (mathematica) we could find exact solution.

The results shows that our approach is close enough to the exact solution except at

final time. We also note that a fast terminal transient has appeared in u(t) and Mt).

 

65

=ou£8 888 was 585% REEEéoc 5053 somEQEoU ofim Ram—m

 

 

pomxm

 

 

 

05:

04 md 0.0 4.0 Nd 0.0
h _ _ _ _ _ _ _ OBVI
.lO.Ntl

10.0

 

 

 

at:
O._. md ed #20 Nd 0.0
_ _ — — F n _ _ L _ _ O.Pl
X 11
Ed:
10.0
Imd
IO.—

 

soeocaq,‘ Eewcwsucoz in.—

66

5.2. Simple Pendulum Problem

 

 

 

Y
11
0 R0
0 :
X
r
R1
T1
1
m,I 0
R0 - r - R1 — 0
§

Figure 5.11 Planar simple pendulum

Consider a simple pendulum shown in Figure 5.11. It is attached to ground by a
revolute joint 0. The figure shows the global reference frame xy and the local refer-

ence frame CT]. The pendulum body has a mass m and inertia I through the center of

mass.

 

67

At any point in time, the body must be in equilibrium with its inertial (D’Alembert)
forces and the forces of constraint holding it to ground. Three generalized
coordinates(i.e. q = [x y 0]T) associated with the pendulum body will be used to for—
mulate the problem. The dynamic equations of motion in the Lagrangian form of
these coordinates are :
ma? + 7&1
i: my + M + mg = 0 (5.9)
m10 + Illsine — szcose — u

where u is a control torque to rotate the pendulum body in plane motion. 1122 are
the rectangular components of the reaction forces required to maintain joint contact,
as exerted on the pendulum from ground. The length of pendulum is 2] and g is

gravity. we assume that the center of mass is in the origin of the pendulum.

Since the constraint requires that the pendulum pin coincide with the ground pin,

from Eq. (2.6), the constraint equations for the revolute joint become :
—x + 10050
(D = -y + lsin0 = 0 (5'10)

The system dynamic equations (5.9) and (5.10) can be expressed in the form of :

A

f

G=¢

=0 (5.11)

 

 

If a corrector formula is employed for integrating Eq. (5.11), then
i
G = g = 0 (5.12)
(I)

where g is a corrector formula.

 

 

 

 

 

68

The Jacobian matrix for this pendulum problem can be written as :

0009—00—10
h
oooofloo—i
h
1
00e00——lslc
h
"710010000
J: _1 (5.13)
0—001000
h
00;h1-00100

 

 

where e = l( llcose + Msine ) and s 5 sin , c 5 cos. Note that Eq. (5.13) results

of when the first order Euler backward integration method is used.

The system dynamic equation (5.9) consists of an non-minimal set of coordi-
nates which can be transformed into a minimal set of coordinates by proper partition-
ing of coordinates q into w and v. Since this system is one degree of freedom, the
dimension of the independent coordinate v is 1. If we choose 0 as an independent
coordinate v and [Jay]T as dependent coordinate W, by the method discussed in sec-
tion 2.1.5, a differential equation in only the independent generalized coordinate 0

can be obtained :

é: — glkcose + flu (5.14)
1 ml2

The Lagrange multipliers can be expressed by :

1.1 = mlc050 02 + 0.75( 3131—9 u —mgc050 sin0) (5.15a)

k2 = mlsin002 + 0.75( mgcosze — %cos@ it ) — mg (5.15b)

 

 

 

 

 

69

The tables (5.1) and (5.2) show the necessary conditions for two different per-
formance measures. The term "non-minimal" shown in the top row of the tables
means that we treat the system constraint equations as the implicit descriptor form of
Eq. (5.9) with a non-minimal set of coordinates by Cartesian approach and the
second order form of necessary conditions discussed in chapter 3 are applied to for-
mulate the conditions. The term "minimal" shown in the tables indicates that the sys-
tem constraint equations are posed on the explicit form of Eq. (5.14) with a minimal
set of coordinates by using the generalized coordinate partitioning method and that
the necessary conditions of the conventional standard state space form are obtained.
The character "A", “S" and "T" in the tables indicate the adjoint equations, the sta—

tionarity and the transversality conditions respectively.
Case (2a) : Minimizing the Control Torque only
Consider minimizing the performance measure :
1 t’
J = — u2 dt
2 1.
Note that Eq. (5.14) should be reduced to the first order form of Eq. (5.14c) and Eq.

(5.14d) in order to obtain the standard space form of necessary conditions for the

"mimimal" case.

 

x, = x2 (5.140)
22 = — flcosxl + 0'75 u (5.14d)
1 ml2

where x1 = 0, x2 = 0

The necessary conditions are shown in the Table 5.1. It looks as if that the minimal

case is simpler than non—minimal. But we have to recognize that in general to

70

transform a non-minimal representation into a minimal representation and to reduce
the second order form of the equations to the standard space form. For Figures 5.12,

5.13, 5.14, 5.15 and 5.16 the following parameters are used. :
m = 10 1= 1 900) = —0.5236 radian (300) = 1

t0 = 0 If: 1 initial guess u0 = 10

 

 

Non-minimal Minimal

"1131 " P4 = 0

ma ‘ P5 = 0 p‘1 = 4.7531- sin 9 02
A 1153 + l (llcose + Msin0)p3 — lsin0p4 + lcoseps = 0 p3 = -p1

[)1 + Isin0 p3 = 0

p2 — Icost) p3 = 0

 

 

 

 

 

S u - = 0 u + — = 0
p3 4m]2 p2
T [71(9) = [12(9) = P3(ff) = 0 P1(If) = 0
1510f) = PfiUf) = FEW) = 0 MW) = 0

 

 

'/
Table 5.1 The necessary conditions for J = % juzdt
10

 

6 (rad)

71

 

 

 

 

 

 

 

 

 

 

0.08 I I I I I I’ T I I
0.04-« -
4'- ‘“ _ T“ - -.
0.00
.I .
-0.O4-4 2
4 —~ non—minimal '
-0.08 ":fimlnIanal r T r l r I r
0.0 0.2 0.4 0.6 0.8 1.0
Time ( sec )
Figure 5.12 The optimal control for case (2a)

000 I I I I I l I l I l
—1.0~ _
—2.01 -

J _
~~ non—minimal
—3.0 ----r mmlrmalf I T r T r r
0.0 0.2 0.4- 0.6 0.8 1.0

Time ( sec )

Figure 5.13 The optimal trajectory for case (2a)

 

80.1——4-»—T-———Ta, . , . , .

 

72

 

 

 

 

  

 

 

 

 

 

‘4 ~— non—minimal ‘ \I
~80. minImal r r T r T f r
0.0 0.2 0.4 0.6 0.8 1.0
Time ( sec )
Figure 5.14 X—direction reaction force for case (2a)
0-“F““"T"‘ TI 1 l . r 1 I r "1
4 4
-50.— —«
4 .
—100.~I _
4 t‘
—-150.~ ~
‘ —~ non—minimal ......... ’ ‘
—2oo. ""4 mini?“ . r r r r r .4
0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.15 Y—direction reaction force for case (2a)

73

 

r 1
minimal I
4» non—minimaljt

O

 

 

O
O
0. T i e "+--“+-“f——¢-~+-~O~-+
5 9 13

No. of Iterations

Figure 5.16 The performance measure for case (2a)
Case (2b) : Minimizing the control and the reaction forces with h term

Consider minimizing the performance measure :
‘/
J=h+%I(u2+x%+7v§)dt (5.16)
to
where
h = 401 00)) — xiv)? + 00,) — My)? + (00)) — 6.0))? 1
+ gm 05(1)) - 40(th + 00)) — 1.0)? + (6(9) — 6.0))? 1

Here at and B are adjusting values which can weight the relative importance of the
deviation of each states from their desired values. The necessary conditions are same

to Table 5.2 except the transversalty conditions.

74

 

Non-minimal

Minimal

 

mp1 ‘P4=0

mpz—P5=0

A Ip'3 + l (Alcose + )vzsin0)p3
— lsin0p4 + lcos0p5 = 0
p1+lsin9p3 + x1: 0

pz—lcosep3+)vz=0

2 2
151 = %usin0 — 15%—cos0sin0

+ m2g102cos0 - Eff-sin 0 p2

152 = —p1 — 2m21203 + 2m2g10sin0

 

 

 

 

 

 

S u-p3=0 l+—9—u+ 3 +3—mg-cose=0
( 1612) 4m12 pz 161
T Piaf) = [92(9) = 193(tf) = 0 plaf) = 0
1010]) = [52(9) = DEW) = 0 P2(‘f) = 0

 

'/
Table 5.2 The necessary conditions for J = % I ( u2 + 7.1+ 765‘ )dt
‘0

For non-minimal approach, the transversality conditions become :

mpl(tf) + l3[x'(tf) — 224(9)] = 0

mpzof) + 1300)) — 11(9)] = 0

IP3(l‘f) + Bl9(tf) - 0d(tf)] .-. 0

mti1(tf) + oc[x(tf) “ X40)” = 0

”1152(9) + aUOf) — yd(tf)] = 0

1153(0) + a190,) - 94(9)] = 0

 

75

For the minimal approach, the transversality conditions become :
121(9) = 0490) — 6.0)] 1220/) = 1116(4) — 6.10)]

For 0,, = —1.57 04(tf) = 0. , Figure 18 shows that the optimal trajectory pursue the
desired final position. Same parameters of case (2a) were used for Figures 17-21.
Those Figures also confirm the applicability of the theory developed in previous

chapters.

9 (rad)

76

804.44.

 

 

 

 

 

 

 

 

 

1 1 _ 1_ 1 I ‘ ‘1' “ 1 ““1‘ ‘ ~—.~ ‘"
,__ non—minimal
* -‘-- minimal ‘
60.~ ~
.1
40.—
1
I
20 H II I, ‘t
i i
04* MT“ r r r m r 1 J1
0.0 0.2 0.4 0.6 0.8 1.0
Time ( sec )
Figure 5.17 The optimal control for case (2b)
0.0—~'—"Vil— 'T'I‘"'T""I’"" 'TIT""" ' I'T‘T’” 'I‘NT“ '1""""T'T‘I'-"’TI‘"' 7"
—0.5—4 --
J Ts
—1.04
-1 .5-1
1 minimal ‘
~— non—minimal
—2.0 1 f 'r I T l r
0.0 0.2 0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.18 The optimal trajectory for case (2b)

N
r<

O. T'"""' 'TTT'T‘T" ““I‘T—T" I ' fl ' "‘7'“ ‘T"r""“‘-I—“ ' "Tr “" " T“_“'i

—4-0.

—-120.

 

l
i

l ~--- minimal
_ 80. 4;“ nonr—minimal f , r Y r fl
0.0 0.2 0.4 0.6 0.8 1.0

r...“ .r,._.._r-._.

 

77

'i““‘ T _I‘T—”T—T~'”I" "‘T‘r’ 'TT’ 7"" 1

 

 

Time ( sec )

Figure 5.19 X-direction reaction force for case (2b)

l

i

—-l

 

0.0

     

---' minimal

M— non—minimal ..................
r Y

T
0.2

—<

 

T r r r r "‘
0.4 0.6 0.8 1.0

Time ( sec )

Figure 5.20 Y-direction reaction force for case (2b)

M

7R

 

 

 

8000. l l I I l I l
. ’.o....’.‘.0.9.o.o.ooooooooooooooooooooooooon
6000.— -
m 4000.— ~
2000.— _
i o non—minimal '

0 minimal
0. T l I I I l I
O 12 24 36 48

 

No. of Iterations

Figure 5.21 The performance measure for case (2b)

5.3. Simple Two Degree of Freedom Manipulator

As shown in chapter 2.1.1, Figure 2.1 is the two degree of freedom manipulator
model. The symbols m1, m2 denote the masses of link 1 and link 2. The 11, J2
respectively are the moments of inertia and g is gravity. The length of each link is

defined as 211, 212.

Consider minimizing the performance measure :

t

J: (u%+ug)dz

'—-.\

i
2

n

0

The mixed differential algebraic dynamic equations can be written by equations

(5.17a) and (5.17b).

 

79

As one follows procedures discussed in chapter 4, the following equations can be

obtained for integration.

0 State Equations

where

7711;] + x1 — l3
m15’.1+ 7V2 _ 14 + ”'18

"hi2 + 7&3
mjz + 14 + ng
1252 + 12( 7L3sin92 — Mcost) -— uz

 

—x1 + llcosel
—'y1 + [13in91
(p =

y1 + llsinel — yz + lzsinez

—l O 1 0

0 —1 0 1 -1 0
—11s1 llcl —llsl [lo] 0 -1 llcl

0 0 —1 o ‘I’q = 1 0

0 0 0 —1 0 1

0 0 —1252 [2C2

s] E sin91 s; E sinez c1 E 00391 c2 5 c0592

1161 + 11[ (k1 + l3)Sin91'—(M + 9&4) C0891 ] - u]

 

x1 + llcosel — x2 + lzcosez =

O 0 0

0 (5.17a)
(5.17b)
0 0 0
O —1252
—1 12C2

 

m10 0 0 0 00

0 mIO 0 o 00

[Egg _ 0 0 h 0 0 00
86 ' o 0 02M 0 00

o 0 0 0 sz 0

0 0 0 0 0.50
000000—1 0
000000 0 —1

a_GT _ 008100 0 —115111C1
[aq] ‘ 00 00()o 0 0
000000 0 0

00 00()Q 0 0

where

el = [1[ (1.1 + X3)COSGI + (l2 + X4)Sin91 ]

oooooo
loooooo.

 

1 0
0 1
“[131 11C]
(5.18)

’1252 1202

e2 = 12[ 1300502 + Msinez ]

31 E sinGl 52 E sinez c1 5 c0301 c2 E 00502
1 0 [151 0 0 0 0 0 0 0
{—3ng _ 0 1 —llc100 0 0000 (519)
al —1 0 [1S] 1 0 [2S2 O 0 O 0 I
0 —1 —llc1 0 1 —12c2 0 0 0 0
[§§f=[0000] 62m

81
o Adjoint Equations
”11131 ‘P7 +199 = 0
"1152 " P8 +1310 = 0
le'3 + 11[ (K1 + k3)cosel + (702 + l4)sin01 ]p3
+ 11[ -sin91p7 + cosGlpg — sin01p9 + coselplo ] = 0

..

"1W4 _ P9 = 0
"12.55 - P10 = 0
12136 + 12[ (l3cosez + K4sin02)p6 — sin02p9 + cos02p10 ] =0
[)1 + llsin01p3 = 0
p2 — llcoselp3 = 0

—p1 + llsin01p3 + p4 + Izsinezpé = 0

—p2 — llcoselp3 + p5 — 12c0502p6 = 0

(5.21a)

(5.21b)

(5.210)
(5.21d)
(5.21e)
(5.211)
(5.21 g)
(5.2111)
(5.211)

(5.21j)

82

0 To get initial and final values

T [q] _ [Maw]
(:10 —1 H0 — “(19(1) t=t°

 

M <I>q
(Dq 0
MT <1)
(I)q 0
F 0
_mlg
“1
0
[f] -ng
Y = ”2
Ilciéi
[isléi
llclf)? + lzczég
11310? + lzszégj

T 00
q [911] _
1:1, pl m,

l

 

 

where
V] = llcoselélp'3
v2 = llcos0101p3
v3 = llsin010¥p3
v7 = 12c059202p'6
v8 = 120059202126

v9 = 1251119263126

v4 = llsin0101p3
v5 = llsin0101p3
v5 : llcoseléim
v10 = [25in0202p'6
v11 = lzsinezézpé

v12 = [zoosezézyé

Rq’q!xapq9bqat)
V(q,qiaquafiq9x) [2:9

0
0

0
0
lzl 1302 + A432 1P6
—2v1—v2+v3
—2v4—v5—v6
—2v1—v2+v3—2v7—v8+v9

_ —2v4—v5—v6—2v 1 O—vl 1—v 1 2

[1[ (A1 + l3)C1 + (12 + x4031 ]p3

 

.l

(5.22)

83

0 Corrector for State

 

 

T. . . .‘(m)
2f. (91+;fl _fi
_g_ 2%. 0 Ad :—g
a“ 3“ -AA <0
Bq 4
"000000-
000000
fii 00e1000
aq 000000
000000
00000e2

 

 

where

e] = 11[ (X1 + 7&3)C0391 + (M + K4)Sin01 ]

 

 

lmlooo 00'
0m100 00
31310011000
(331' fifi‘Zooomoo
0 000m20
0 000 OJ2

a? 3" —1a*
-5x=<1>§ 35:71 i=1

.—1 1 .

g = TqIH-l + 7111;: + C[n+1

(5.23)

62 = [z[ 13(30892 + 91.4sin02 ]

 

 

84

o Corrector for adjoint

_a? a_"f_;_a“f_ a_“f"")
aliq abq 0B0 35¢] 3px qu (m) .f (m)
523- _8g_ 0 qu =—'g (5.24)
apq apq -
~ Apl (I)
M)
— 0 0
apq
000000
000000
a? 0061000
a—pq‘000000
000000
00000122
m100000
Om10000
8f 13? 10011000
ap_hBOa'—,5q)"7000m200
0000m20
0000012

 

 

8px q apq h abq
__a&) = (I) ~ — __1 + l + '
apq q g h pct»: h pq. pqn

where

p1 + llsin01p3

P2 " l100391P3
—p1 + llsin01p3 + p4 + lzsin02p6
—p2 - llcoselp3 + p5 — [2C0362p6

85

 

 

 

 

 

 

0010 I I I l l [ l I I
l _
0.05— -
1 _
0.00— —
__.... “2
— u
”0.05 l I 1 l I l l l I I
0.0 0.2 0.4 0.6 0.8 1.0
Time ( sec )
Figure 5.22 The optimal controls for manipulator model
1. . 1 - 1 . 1 1 1

6 (rad)

 

92

 

 

 

0.0

1
0.2

r l
014 0.6

Time ( sec )

Figure 5.23 The optimal trajectories for manipulator model

86

 

200.

-200.

 

 

\
I
\‘ I

 

 

—400. I

r 1
0.4 0.6

Time ( sec )

Figure 5.24 The first joint reaction forces for manipulator model

 

r< —100.—~

 

 

l'

 

 

—300. I

1 1
0.4 0.6

Time ( sec )

r
0.8

Figure 5.25 The second joint reaction forces for manipulator model

87

 

 

 

 

120. l I T l I r I I I
1» 7
90.— —
60w —
9
30A —
_ . _
O . 1 1 9, 1 fi : a a
1 3 5 7 9 11

No. of Iterations

Figure 5.26 The performance measure for manipulator model

CHAPTER 6

SUMMARY AND CONCLUSIONS

We have introduced more flexible forms of necessary conditions for optimal

control problems on multibody mechanical systems models of non-minimal or

descriptor form. After demonstrating the regularity of the optimization problem,

necessary conditions of descriptor form were derived using the standard Lagrange

multiplier theorem of constrained optimization. We can state several advantages to

the approach taken.

(a)

(b)

(c)

A minimal coordinate model is not needed in order to pose optimal control

problems.

The performance measure can be posed on a broader and larger coordinate set.
We have emphasized this by choosing a performance measure which depends on

a Lagrange multiplier in the example problems.

It was shown that the system equations and the adjoint equations are of similar
form. This is also observed in Barman[50]. This results in simplications in for-
mulating the necessary conditions since similar methods for formulating the both
the system equations and the adjoint equations can be used. Computational
advantages also result. The approach made effective application of optimal con-

trol methods to non-minimal coordinate systems.

As we have seen in the numerical results in chapter 5, the goal of this disserata-

tion was achieved satisfactorily. In other words the variational approach was properly

applied to the descriptor systems. But results have pointed out difficulties either

numerical or theoretical. Singular perturbation methods [51] may be the solution. An

important question is how to set up a proper formulation of the multi-body equation

88

89

for a singular perturbation analysis of fast transient and their minimization. This

issue is discussed in [52]. See also reference [27].

 

LIST OF REFERENCES

LIST OF REFERENCES

[1] Greenwood, D.T., "Principles of Dynamics," Prentice-Hall, Englewood Cliffs,
New Jersey, 1965.

[2] Kane, T.R., "Dynamics," Holt, Rinehart, and Winston, N.Y., 1968.

[3] Goldstein, H., "Classical Mechanics," Addison-Wesley, Reading, Massachusetts,
1980.

[4] Orlandea, N., Chace, M. a., and Calahan, D. A., "A Sparsity-Oriented Approach
to Dynamic Analysis and Design of Mechanical Systems, Part I and II," ASME J.
Engr. Indust., Vol. 99, 773-781, Aug. 1977.

[5] Wehage, R. A., and Haug, E. J., "Generalized Coordinate Partitioning of Dimen-
sion Reduction in Analysis of Constrained Dynamic Systems," ASME J. Mech.
Design, Vol. 104, 247-255, Jan. 1982.

[6] Nikravesh, PE, and Chung, I. S.,"Application of Euler Parameters to the
Dynamic Reduction in Analysis of Constrained Dynamic Systems,” ASME J.
Mech. Design, Vol. 14, 785-791, Oct. 1982.

[7] Nikravesh, P. E., "Some Methods for Dynamic Analysis of Constrained Mechani-
cal Systems: A Survey," Computer Aided Analysis and Optimization of Mechani-
cal System Dynamic Systems, E. J. Haug, ed., Springer-Verlag, Heidelberg, West
Germany, 1984.

[8] Baumgarte, J., "Stabilization of Constraints and Integrals of Motion in Dynamic
Systems," Computer Methods in Applied Mechanics in Engineering, Vol. 1- 16,
1972.

[9] Petzold, L., "Differential- -A1gebraic Equations are Not ODE’s " SIAM J. of
Scientific and Statistical Computing, Vol. 3,No.3, 1982

[10] Sage A. P., "Optimum Systems Control", Prentice-Hall, Englewood Cliffs, New
Jersey, 1977.

[11] Gear, C. W., and Petzold, L. R., "ODE Methods for the Solution of
Differential/Algebraic Systems," SIAM J. Number Analysis, Vol. 21, NO. 4,
716-728, 1984.

[12] Gear. C. W., "Differential-Algebraic Equations," Computer Aided Analysis and
Optimization of Mechanical System Dynamics, ed. E. J. Haug, Springer Verlag,
Heidelberg, 1984.

[13] Haug, E. J., Wehage, R. A., and Barman, N. C., "Design Sensitivity Analysis of
Planar Mechanisms and Machine Dynamics," ASME Journal of Mechanical

90

91

Design, Vol. 103, No. 3, July 1981.

[14] Hachtel, G. D., Brayton, R. K., and Gustavson, F. G., "The Sparse Tableau
Approach to Network Analysis and Design," IEEE Transactions on Circuit
Theory, Vol. CT-l8, No. 1, Jan. 1971.

[15] Gear, C.W., "Numerical Initial Value Problems in Ordinary Differential Equa-
tions," Prentice-Hall, Inc., 1971.

[16] Chua, L.O., and Lin, P.M., "Computer-Aided Analysis of Electronic Circuits Cir-
cuits", Prentice-Hall, Inc., Englewood Cliffs, NJ, 1975.

[17] Grenin, y. ., "A New Approach to the Synthesis of Stiffly Stable Linear Multistep
Formulas", IEEE Trans. Circuit Theory, 352- 360, July 1973.

[18] Gear, C. W., "The Automatic Integration of Stiff Ordinary Differential Equa-
tions", Information Processing 68, A.J.H. Morrell, Ed., North Holland Publishing
Company, Amsterdam, The Netherlands, 187-193, 1969.

[19] Gear, C. W., "Simultaneous Numerical Solution of Differential Algebraic Equa-
tions", Trans. Circuit Theory, Vol. CT-18, No. 1,89- 95, Jan. 1971.

[20] Keller, H.B., "Numerical Methods for Two-Point Boundary Problems," Blaisdell,
1968.

[21] Lambert, J.D., "Computational Methods in Ordinary Differential equations," John
Wiley & Sons, 1973.

[22] Osborne, M.R., "On shooting methods for boundary value problem,” J. Math.
Anal. Appl., Vol. 27, 417-433, 1969.

[23] Arthur E. Bryson, Jr., "Applied Optimal Control", John Wiley & Sons, 1975.
[24] Donard E. Kirk, "Optimal Control Theory", Prentice-Hall, 1970.

[25] Campbell, S. L., "Singular Systems of Differential Equations", Vol.1, Pitman,
London, 1980.

[26] Campbelll, S. L, "Singular Systems of Differential Equations", Vol. 2, Pitman,
London, 1.982

[27] Cobb, D., "Descriptor Variable Systems and Optimal State Regulation", IEEE
Trans. Autom. Contr, Vol. AC—28, 601-611, May 1983.

[28] Verghese,G., Levy, B. C., and Kailath, T., "A generalized State- Space for Singu-
lar Systems", IEEE Trans.Autom. Contr, Vol.AC—26, 811-831, Aug.198l.

92

[29] Lovass-Nagy, V., Schilling, R., and Yan, H.-C., "A note an optimal control of
generalized state-space(descriptor) systems", Int. J. Control, Vol. 44, 613-624,
1986.

[30] Luenberger, D. G, "Dynamic equations in descriptor form," IEEE Trans.

Automat. Contr, Vol. AC— 22, 312- 321, Mar. 1977.

[31] Luenberger, D. G.,"Time-Invariant Descriptors Systems," Automatica, Vol. 14,
473-480, 1978.

[32] Wu, H.S., "Generalized maximum principle for optimal control of generalized
state—space systems," Int. J. control, Vol. 47, 373-380, 1988.

[33] Ibrahim, E.Y., Lovass-Nagy, V., Mukundan R. and Schilling R.J., "Free final-time
optimal control of descriptor systems," Int. J. Control, Vol. 48, 407-416, 1988.

[34] Pandolfi, L., "Controlability and Stabilization for Linear Systems of Algebraic and
Differential Equations," J. Optimization Theory Applic., Vol. 30, April 1980.

[35] Luenberger, D. G., "Optimization By Vector Space Methods," John Wiley &
Sons,Inc. New York 1969.

[36] Pontryagin, L. S., V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko,
"The Mathematical Theory of Optimal Processes(trans, by K.N. Trirogoff, ed. by
L. W. Neustadt)," John Wiley, New York, 1962.

[37] Luistemik, L. and V. Sobolev, "Elements of Functional Analysis," Frederick
Unger, New York, 1961.

[38] Gelfand, I. M., and S. V. Fomin, "Calculus of Variations," Englewood Cliffs,
NJ, Prentice-Hall, Inc., 1963.

[39] Rudin, W., "Functional Analysis," Mcgraw-Hill, Inc., 1973.

[40] Blum E. K. ,"The calculus of Variations, Functional Analysis, and Optimal Con-
trol Problems," Topics in Optimization, Ed. by G. Leitman, Academic Press,
New York, 417- 461, 1967.

[41] Swisher G. H.,"Introduction to Linear Systems Analysis," Matrix Publishers, Inc.

[42] Paiewonsky,"0ptimal Control: A Review of Theory and Practice," AIAA Journal,
Vol. 3, 1985-2006, 1965

[43] Athans M. and P. L. Falb, "Optimal Control: An Introduction to the Theory and
It’s Applications," McGraw—Hill, Inc., 1966.

93

[44] Nikravesh, P. E.,"Computer-Aided Analysis of Mechanical Systems," Prentice-
Hall, Englewood Cliffs, N.J., 1988.

[45] Paul, B. and Krajcinovic, D., "Computer Analysis of Machines with Planar
Motion-Part I: Kinematic ; Part H: Dynamics," Journal of Applied Mechanics,
Vol. 37, 697-712, 1970

[46] Chace, M.A. and Smith, D.A., "DAMN-A Digital Computer Program for the
Dynamic Analysis of Generalized Mechanical Systems," SAE paper 71024, Janu-
ary 1971.

[47] Sheth, P. N. and Uicker, J. J., Jr., "IMP(Integrated Mechanisms Program), A
Computer Aided Design Analysis System for Mechanisms and Linkages," Journal
of Engineering for Industry, Vol. 94, 454-464, 1972

[48] Athans, M., "The Status of Optimal Control Theory and Applications for Deter-
ministic Systems," IEEE Trans. On Automatic Control, Vol. AC-11, No. 3, 580-
596, 1966.

[49] Haug, E. J., "Elements and Methods of Computational Dynamics," : A Survey,
"Computer Aided Analysis and Optimization of Mechanical System Dynamic
Systems, E. J. Haug, ed., Springer-Verlag, Heidelberg, West Germany, 1984.

[50] Barman N. C., "Design Sensitivity Analysis and Optimization of Constrained
Dynamic Systems," Ph.D. Diss. The University of Iowa, 1979

[51] Kokotovic, Khalil, O’reilly, "Singular Perturbation Methods in Control,"
Academic Press, New York, 1986

[52] Joe H. Chow, John J. Allemong and Petar V. Kokotovic,"Singular Perturbation
Analysis of Systems with Sustained High Frequency Oscillations," Automatica,
Vol. 14, 271-279, 1978

 

 

 

HICH

   

3

11111111111111111111111111]1E5