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( 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]. Mii — qTx = f(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 (q.t) = rem). 2(q.r),....<1>m 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 = —, (2.12) Since equation (2.11) can be expressed as , (2.13) Since the coefficient matrix wa of w in Eq. (2.13) is nonsingular, the dependent velocity w can be expressed as iv = -,,-1,& — (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 — 2mw — 2¢,,9 - (1),, where we,» a awe + car 19 Since the (1),, is nonsingular, the dependent acceleration at can be expressed as 66 = —,,-1[ cm + ywv'v + 7,6 + 2,,,€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 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, dr5q + n 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> (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.r) 1 } dr + htqoxioro ) (3.6) ‘0 so that the performance measure becomes : ’/ J(um) = j{ g.<1(z).1,t) + $144611) 1 } at + h ’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.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 = — 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) + 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_§(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) = _; 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 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 + £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) = —(qq = [<¢,pq>q4]qd + <¢qpq>qii + 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 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