The U n i v e r s i t y of B r i t i s h Columbia FACULTY OF GRADUATE STUDIES PROGRAMME OF THE FINAL ORAL EXAMINATION FOR THE DEGREE OF DOCTOR OF PHILOSOPHY , of JAMES WILLIAM SUTHERLAND B.A.Sc., Un i v e r s i t y of B r i t i s h Columbia, 1964 MONDAY, MAY 1, 1967, AT 10:30 A.M. IN ROOM 402, MacLEOD BUILDING COMMITTEE IN CHARGE Chairman: B. N. Moyls E. V. Bohn' Y. N. Yu V. J. Modi A. C. Soudack R. W. Donaldson M. Davies External Examiner: Dr. I. H. Mofti Analysis Laboratory National Research Council, Ottawa, Ontario. Research Supervisor: E. V. Bohn THE SYNTHESIS OF OPTIMAL CONTROLLERS FOR A CLASS OF AERODYNAMICAL SYSTEMS, AND THE NUMERICAL SOLUTION OF NONLINEAR OPTIMAL CONTROL PROBLEMS ABSTRACT This d i s s e r t a t i o n i s divided into two parts. In Part I, a method i s developed for determining the o p t i -mal c o n t r o l laws fo r - a c l a s s of aerodynamical systems whose dynamics are l i n e a r i n the thrust and nonlinear in the l i f t and thrust angle. Conditions under which the adjoint variables can be eliminated from the control equations are derived, and expressions for the thrust and rate of change of l i f t and thrust angle are obtained which depend ,only'on state variables and a small number of time invariant parameters. The optimal values of the unknown parameters are determined by a d i r e c t search i n parameter space. I t i s shown that the proposed tech-nique i s considerably'simpler than standard gradient techniques which require a separate search i n function space for each component of the c o n t r o l vector. Further-more, since the controls are generated by the d i r e c t s o l u t i o n of d i f f e r e n t i a l equations, the method appears suitable f or use with i n - f l i g h t guidance computers. In Part I I , a three stage numerical algorithm i s developed for a general cl a s s of optimal c o n t r o l prob-lems. The f i r s t two stages of the algorithm are based on a gradient search i n the parameter space of i n i t i a l Lagrange m u l t i p l i e r s . The f i r s t stage.attempts to s a t i s f y the given end constraints without regard to system performance, and the second stage attempts to im-prove the system performance while simultaneously main-taining the end constraints set by the f i r s t stage. The f i n a l stage of the algorithm i s based on e i t h e r a modi-f i e d method of matching end points, or a method of determining the optimal step size for the gradient method of the second stage. E i t h e r combination r e s u l t s in a three stage algorithm which has good i n i t i a l con-vergence, good f i n a l convergence, and which requires storage at terminal points only. GRADUATE STUDIES F i e l d of Study: E l e c t r i c a l Engineering Analogue Computers Nonlinear System Theory Network Theory Servomechanisms Applied Electromagnetic Theory Communications Theory. D i g i t a l Computers E.V. Bohn A.C. Soudack L. Mason E.V. Bohn M.M. Kharadly A.D. Moore H.P. Zeiger Related Studies: The Calculus of Variations Space Dynamics J.F. Scott-Thomas V.J. Modi PUBLICATIONS Bohn, E.V., Chan, W.C., and Sutherland, J.W., "The synthesis of r e a l time optimal c o n t r o l l e r s for rockets", Astronautica Acta, 12, pp 26-32, 1966. Sutherland, J.W., and Bohn, E.V., "A numerical t r a j e c -tory optimization method suitable for a computer of l i m i t e d memory", IEEE Trans, on Automatic Control, Vol. AC-11, Number 3, July 1966, pp 440-447. THE SYNTHESIS OP OPTIMAL CONTROLLERS POR A CLASS OP AERODYNAMICAL SYSTEMS, AND THE NUMERICAL SOLUTION OP NONLINEAR OPTIMAL CONTROL PROBLEMS by JAMES WILLIAM SUTHERLAND B.A.Sc, University of B r i t i s h Columbia, 1964 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS POR THE DEGREE OP DOCTOR OP PHILOSOPHY i n the Department of E l e c t r i c a l Engineering We accept t h i s thesis as conforming to the required standard Research Supervisor Members of the Committee . . . . . a o o o a o . f i . . 0 0 e o o o t o o o o o o o o o o o o Head of the Department o o ~ o » o » o o o o o o o o o » o o o Members of the Department of E l e c t r i c a l Engineering THE UNIVERSITY OF BRITISH COLUMBIA . May, 1967 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r an a d v a n c e d d e g r e e a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e Head o f my D e p a r t m e n t o r by h i s r e p r e s e n t a t i v e s . . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . D e p a r t m e n t o f T h e U n i v e r s i t y o f B r i t i s h C o l u m b i a V a n c o u v e r 8, C a n a d a ABSTRACT In Part I, a method i s developed for determining the optimal control laws for a class of aerodynamical systems whose dynamics are l i n e a r i n the thrust and nonlinear i n the l i f t and thrust angle. Due to the presence of the l i n e a r thrust control, a singular subarc exists along which i t i s often possible to eliminate the Lagrange m u l t i p l i e r s from the control equations. Conditions under which t h i s elimination i s possible are derived, and expressions for thrust and the rate of change of l i f t and thrust angle are obtained that depend only on state variables and a small number of time-invariant parameters. The optimal values of the unknown parameters are determined by a direct search i n parameter space for that set which minimizes the system performance function. As a r e s u l t , the proposed method i s considerably simpler than standard numerical techniques that require a separate search i n function space for each com-ponent of the control vector. Furthermore, since the control vector i s generated by the direct solution of d i f f e r e n t i a l equations, the method appears suitable for use with i n - f l i g h t guidance computers. Several numerical examples are presented consisting of one, two, and three dimensional control. In each case, i t i s shown that the search i n multi-dimensional function space can be replaced by an equivalent search i n the parameter space of i n i t i a l conditions. In Part I I , a three stage numerical algorithm i s developed for a general class of optimal control problems. The technique i s es s e n t i a l l y a combination of the direct and i i i n d i r e c t approaches. Like the i n d i r e c t approach, the control law equations are used to eliminate the control vector from the system and adjoint equations. However, instead of trying to solve the two point boundary-value problem d i r e c t l y , the aug-mented performance function i s f i r s t considered to be a func-t i o n of the unknown i n i t i a l conditions and i s minimized by a gradient search i n the i n i t i a l condition space. I t i s shown that i t i s s u f f i c i e n t to search over the surface of any sphere * * for the intersection of the l i n e u\ . where A. i s the c l a s s i c a l ^ o 7 o solution of i n i t i a l values. As a r e s u l t , t h i s f i r s t approach i s not dependent on a good i n i t i a l estimate of the optimal t r a -jectory, and i s therefore used i n the f i r s t two stages of the proposed algorithm to provide the property of rapid i n i t i a l con-vergence,. The property of rapid f i n a l convergence i s obtained by employing either a modified method of matching end points, or a method of determining the optimal step size for the gradient method of the f i r s t two stages. Either combination results i n a three stage numerical algorithm that has good i n i t i a l con-vergence, good f i n a l convergence, and which requires storage at terminal points only. Several examples are presented consisting of both bounded and unbounded control. i i i TABLE OF CONTENTS Page LIST OP ILLUSTRATIONS v i i i ACKNOWLEDGEMENTS . ... .... ....... x 1. INTRODUCTION ., ., :. ,, . , ., . , , ,.' „ .. ; „ , . ., ,. „•„ , ., , . . ., 1.1 The Indirect Approach 1 1.2 The Direct Approach 2 1.3 The Dynamic Programming Approach 3 1.4 The Proposed Techniques 4 PART I OPTIMAL CONTROL LAWS POR A CLASS OP AERODYNAMICAL SYSTEMS 2. GENERAL ANALYSIS 2.1 Optimal Control of Aerodynamical Systems 7 2.2 Problem Statement 7 2.3 Discussion of the Necessary Conditions for an Ojp"bXlTlSll Tx*£t J S C "b 02?y • o e o o o 0 t > * o o » * o 0 0 0 * O 0 o o o o « * > o « o o 1 "1 2.4 The Control Equations for Thrust, Thrust Angle, 5. THE SOUNDING ROCKET PROBLEM 3.1 Derivation of Optimal Control Law 19 4. THE MAXIMUM RANGE PROBLEM 4.1 Derivation of Optimal Control Laws for the case 4.2 Derivation of Optimal Control Law for the case L = 0 34 4.3 The Maximum Range Problem with Control of F i r i n g 4.4 The Fixed End Point Problem 47 4.5 The Direct Search by the Modified Relaxation 5. PROBLEM FOR WHICH c x = 0 5.2 The Three Dimensional Control Problem 57 5.3 The Class of Problems for which 0-^ = 0 60 i v Page 5.3.1 The Case of Thrust Control Only 61 5.3.2 The Case of Thrust and Thrust Angle Control ... 62 5.3.3 The Case of Thrust, Thrust Angle, and L i f t Control 62 5.4 The Maximum Altitude Problem for the Case of Thrust and Thrust-Angle Control 63 5.4.1 Approximate Solution for u m o„ F i n i t e 65 5.4.2 The Exact Solution for u o v F i n i t e 65 max. 5.4.3 Comparison of Exact and Approximate Solutions 66 5.5 Testing the Necessary Conditions of the Calculus of Variations 66 5.6 The Problems of Maximizing the Function G^ at 5.6.1 Computing Algorithm for Thrust Control Only ... 75 5.6.2 Computing Algorithm for Thrust and Thrust-Angle Control 75 5.6.3 Computing Algorithm for Thrust, Thrust-Angle, and L i f t Control 77 5.6.4 A Comparison of the Three Cases 79 6. SUB0PTIMA1 CONTROLS 6a 1 Iii.ijX,oc3.'U-C"tion • © © < » a « © © o © © © © * o « o © © © © o o o « * © © « o * « » © * * 32 6.2 Eliminating ^ from the Control Equations 82 6.3 Using (3 = 0 i n the Two Dimensional Control Case.. 83 6.4 Comparison of Suboptimal and Optimal Controls ... 85 7. OPTIMAL CONTROLLERS DURING- SINGULAR SUBARC Y o 1 In~b3TOCLULC"fclOri © © O 8 © o o f t » « « © © o o o © » » o » © e « « o » o » » « © © e # » 86 7.2 Direct Feedback Control 86 7.3 Hybrid Optimal Controllers 86 7.4 Impl i c i t Function Generation 88 7.5 Conclusion 88 PART I I NUMERICAL ALGORITHMS 8. NUMERICAL TECHNIQUES r • 8 a 1 In~fc3? OCLULC "t i O H « « o © 6 0 O « « o o o o o o » © © o o o © o © * » « © « o « « a o o © 91 8.1.1 The Direct and Indirect Approaches 93 8.2 Linear Time-Varying D i f f e r e n t i a l Systems 93 v Page 8.2.1 The Zero Input Response 94 8.2.2 The Forced Response 94 8.2.3 The Adjoint System 95 8 . 3 Linearization About a Nominal Trajectory 96 8.4 The Optimal Control Problem 97 8.4.1 The Necessary Conditions for a Local Extremum.. 97 8.5 The Method of Steepest Descent 101 8 o & Til.© Mill -H S " b 3 T 3 , " f c 6 ^ y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 102 8.7 The Newton-Raphson Technique 105 8.8 The Method of Matching End Points 107 9. AN ALGORITHM BASED ON FIRST VARIATION 9.1 Introduction a . . o o » . o o . o o « o o « . « o » » o « a . . . . « « . . » . . . I l l 9.2 The Proposed Algorithm I l l 9 .3 Extension of the Proposed Technique 115 9.3.1 An Algorithm for Numerical Computation 120 9.4 The Arbitrary Scaling of the Lagrange M u l t i - .' ] 0 l X© Z * S o o o o o o o o o o o o o o o o o o o o o o o o o o o a o ' o o o o o o o * * * * * 122 5 0 5 EX£lHl]pl© 1 0 0 0 0 0 0 O 0 0 0 O O O 0 O 0 O O O 0 O O 0 O O O 0 0 0 0 O 0 0 0 0 0 0 0 0 12 5 9 » 6 EXcHlipl© 2 O O O O O O O O O O O O O O O O A O O O O O O O O O O O A O O O O O O O O O * 129 9 ° T E X c l l T l J O l © 3 O O O 0 0 O O O O O O 0 O O O O 0 0 O O O O O O O O O O O 0 O O O O 0 O O O 0 13 4" 10. TECHNIQUES FOR IMPROVING FINAL CONVERGENCE 10.1 Matching End Points Using an Optimal Scale FStCijODT 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 0 1^^ 10.1.1 Extension of the Method of Matching End P O X H ~ f c S O O O O 0 O O O O O O 0 O 0 O O O 0 O 0 O O 0 O O O O O O O O O O O 0 O O O O 1^ "1 10.1.2 Computing Algorithm F l Using Matching End POXH"fcS 0 O O 0 0 O O 0 O O O O O 0 O O O O O 0 O 0 O O O 0 0 O O O O O 0 O O O O 0 O 14" 5 10.2 Determining Optimal Step Size for the F i r s t Variation Approach 146 10.2.1 The Second Variation Method for Determining the Optimal Step Size 150 10.2.2 Computing Algorithm F2 Using Second Vstx*Xci~fcxo3x o o o o o o o o o o e o o o o o o o o o o o o o o o o o o o o o o o t o 154* 10.2.3 A Curve F i t t i n g Technique to Determine the Ojpijimstl S~fc©jo Sxz© 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15A 10.2.4 Computing Algorithm P3 Using a Curve F i t t i n g T© C l X U X C ^ l i © 0 0 0 O 9 9 0 0 O 0 O 0 O 0 O 0 O 0 0 0 0 O 0 0 0 0 0 0 O 0 0 0 0 0 0 0 15^ v i Page 10.3 The Combined Algorithm 157 10.3-1 Example 1 157 10•3°2 Exsmipls 2 • « o « * * * » * o o © © « o « « © « © © * » o » » o « < > o « « * » o o 160 10.4 The Second Variation Technique Used With the Method of Steepest Descent 166 10.5 Conclusions 169 Appendix A . . . . . . o o © e o 9 « o . o a © . o * « . « « o « . . o . . . « . . . . . . . . . 172 AppendiX B o o o . . . e . « o o e o 0 0 . « o o o o e o o » • . . o o . . . . 175 Appendix C 178 Appendix D . a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . • • • • • • 183 Appendix E 189 Appendix P 190 Appendix G* . . o o . . . . . . . . . ' . . . . o . . . . . . . . . . . . . . . . . . . . . . . . . 191 Appendix H 191 REFERENCES 195 v i i LIST OF ILLUSTRATIONS Figure Page 2 . 1 The Motion of a Rocket Within the Earth's 8 3 . 1 Sequence Diagram for the Sounding Rocket with Impulsive Boosting 22 3o2 Sequence Diagram for the Sounding Rocket with F i n i t e Maximum Thrust 25 4 . 1 The F i n a l Range as a Function of v'g and X g 32 4 .2 u(t) and $(t) f or the Optimal Trajectory 33 4 .3 The F i n a l Range as a Function of v g and (3Q 38 4o4 The Optimal Controls u and 3 for the maximum "RSOLI^Q J?X*OlolSni e o o o o o o o o 9 « 9 o o o o o * o o 0 o « * o « e o o o o o 9 o 39 4 . 5 The Optimal Trajectories for the u Control, the (u,p) Control, and the B a l l i s t i c Trajectory .... 4 1 4 . 6 44 4 . 7 The F i n a l Range as a Function of X v>. 45 4 . 8 The Optimal Control u for the! Maximum Range Problem with Free F i r i n g Angle 46 4 . 9 49 4 . 1 0 One Dimensional Search Employing a Parabolic 51 4 . 1 1 A. Standard Relaxation Search i n the ( (3 Q , v ) •Plctl!.© • e o o o 0 * « o o o * « e o e e « o o « e o e o « o o o e o o « o e e « « « o * « 53 4 . 1 2 The Modified Relaxation Search i n the (p , v ) Pistil© • Q 0 o o o o 6 o e « . o o o « e o * * o e « « o « o o « o o * o * « o o e o « * « * 55 4 . 1 3 Contour Map About the Optimum Found by the Modified Relaxation 56 5 . 1 The F i n a l Altitude as a Function of |3 Q 70 5.2 70 5.3 The Optimal Lagrange M u l t i p l i e r s for the Maximum Altitude Problem . . . . . . . . a . . . . . . . . . . . . . . . . . . . . . . 71 5.4 The Control Law Equations for the Maximum A l t i -"fcULCL© PX*0"bl©Hl » » « © © » o © 0 0 0 O O 0 0 « » 0 O 0 O f l 0 0 0 0 0 © 0 0 O 0 0 e 0 72 5.5 The Performance G, as a Function of v„ ......... D S 76 5.6 The Optimal Control u, (max. energy) 76 5.7 The Performance as a Function of (3Q 78 5.8 The Optimal Controls u and (3, (max. energy) .... 78 v i i i Figure Page 5.9 The Performance as a Function of P Q 80 5.10 The Optimal Controls u, |3, and L, (max. energy). 80 7.1 The Direct Feedback Optimal Controller 87 7.2a The Hybrid Optimal Controller - Maximum Range with Thrust Control 87 7.2b The Hybrid Optimal C o n t r o l l e r - Maximum Range with Thrust and Thrust Angle Control 87 7.3 Implicit Solution Using a High Gain Amplifier .. 89 7.4a Implicit Solution Controller - Two Dimensional Control with c 1 = 0 89 7.4b Implicit Solution Controller - Three Dimensional Control with c, =0 89 9.1 The Solution Line \x\Q 124 9.2a Transition Flow Graph for Example 1 125 9.2b Transition Flow Graph for Example 3 - 125 9.3 Iterative Path i n I n i t i a l Condition Space of Lagrange M u l t i p l i e r s 128 9.4 Temperature P r o f i l e s for Various Iteration ; Cycles 130 9.5 Xg P r o f i l e for Various Iteration Cycles 130 9.6 The Time Variations for the Optimal Path (E x . l ) . 131 9.7 Temperature P r o f i l e s for Cases a to b (Ex. 2)... 131 9.8 The Time Variations for the Optimal Path (Ex. 2) 133 9.9 The Time Variations for the Optimal Path ( EX o 3 ) a « o « 0 o o » > o o o o o e o o o e o e o a 0 o o o s « o « e o e e o o e o * a 1 3 3 10.1 The Solution Line \i\Q and the Cone of Conver-j££©l'lC 6 -R. 9 o a 0 o o 0 o e 0 O O 0 e 0 0 e o o o e 0 O 0 o e o 0 O 0 0 O 0 0 0 i > a « o 3-39 c 10.2 Convergence to \ by the Modified Method of Matching End Points 140 10.3 The Regions Rg, RJ , and Ra About the Solution . # a Line u-\^ . o . . . . . . . . . . . . © o . . . . . . . . . . . . . . o e o . . . . . . 149 10.4 The Graph of J (Soc) i n the Neighbourhood of the n 10.5 Flow Chart for the Combined Algorithm 158 10.6 The Optimal Trajectory for Example 2 162 10.7 The Optimal Lagrange M u l t i p l i e r s for Example 2.. 163 10.8 J i n the.Neighbourhood of the n-th Step 165 El i x ACKNOWLEDGEMENT I wish to express my gratitude to the Royal Canadian Ai r Force for making th i s graduate programme possible. Also, I wish to thank my supervisor, Dr. E. V. Bohn, and the head of thi s department, Dr. F. Noakes, for t h e i r continued interest and guidance throughout t h i s course of study. Appreciation i s extended to Dr. Y. N. Yu, Dr. V. J. Modi, and Mr. B. Wilbee for reviewing t h i s thesis and l e v e l l i n g constructive c r i t i -cisms. A special debt of gratitude i s held for my wife Susan and my son James Dean for t h e i r patient understanding and unselfish support over the many years of study leading towards the preparation of t h i s thesis. Acknowledgement i s g r a t e f u l l y given to the National Research Council for providing an N.R.C. Bursary, 1964-1965, and two N.R.C. Studentships, 1965-1967. x 1. I. INTRODUCTION 1.1 The Indirect Approach The c l a s s i c a l theory of the calculus of variations was developed nearly two-hundred years ago by Euler and Lagrange. Recently, a more complete and mathematically rigorous treatment of optimal control theory was presented by Pontryagin i n the form of the maximum pri n c i p l e [ l ^ j . Both techniques are essen-t i a l l y equivalent and form the basis of the indirect approach to the solution of optimal control problems. The indirect approach i s based on establishing a set of conditions which are necessary but not s u f f i c i e n t for an extremum. The solution to the optimal control problem i s then taken as any trajectory along which these necessary conditions are s a t i s f i e d . As a r e s u l t , the second v a r i a t i o n may have to be computed to insure a true extremum of the solution and not merely a stationary point. Although this approach often allows some a n a l y t i c a l information to be obtained about the optimal trajectory, a complete a n a l y t i c a l solutions i s usually not pos-s i b l e and; therefore, numerical techniques are employed i n associ-' ation with the indirect approach. However, these numerical techniques usually require a high-capacity d i g i t a l computer and usually need a good i n i t i a l estimate of the optimal trajectory before the solution w i l l converge. On the other hand, once a trajectory i s found which i s i n the neighbourhood of the extre-mum, very rapid f i n a l convergence i s realized. Typical examples of numerical techniques based on the i n d i r e c t approach are the two-point boundary-value problem |J2, 3] 9 "the successive sweep method and the min-H strategy 1.2 The Direct Approach A further computational scheme available to solve optimal control problems i s the method of steepest descent (or ascent) which i s based on a standard h i l l - c l i m b i n g approach. In contrast to the ind i r e c t approach, which s a t i s f i e s the necessary conditions for an extremum, th i s h i l l - c l i m b i n g tech-nique seeks d i r e c t l y that trajectory along which the system performance i s an extremum. As a r e s u l t , the method of steepest descent i s known as a direct approach to the solution of optimal control problems. However, unlike most h i l l - c l i m b i n g techniques, which are based on a gradient search i n parameter space, t h i s approach i s based on an i t e r a t i v e scheme for improving the con-t r o l function by means of a gradient search i n function space. Consequently, as these control functions must be stored at many points along the trajectory, the r e s u l t i n g memory requirements of the computer may become excessively large. Furthermore, since the technique i s predicated on a gradient method, the convergence slows down as the extremum i s approached and i t i s not known when the search should be terminated. On the other hand, a main advantage of the direct approach i s that the solu-t i o n does not depend upon a good i n i t i a l estimate of the optimal trajectory. In f a c t , once a nominal solution i s obtained, i n i t i a l convergence i s guaranteed. direct method are the method of gradients by Kelly [_6j and the method of steepest descent by Bryson and Dehham 7 . Typical examples of numerical techniques based on the 3. 1.3 The Dynamic Programming Approach The t h i r d approach to the solution of optimal control problems i s the dynamic programming approach recently developed by Bellman [ 8 j . This approach i s based on the pr i n c i p l e of optimality which states that "an optimal policy has the property that whatever the i n i t i a l state and i n i t i a l decision are, the remaining decisions must constitute an optimal policy from the state r e s u l t i n g from the f i r s t decision". Employing t h i s p r i n c i p l e , the dynamic programming technique works backward from the desired f i n a l conditions and evaluates the optimal controls at discrete points i n the entire space of permissible states. This flooding of the solution throughout the state space has, i n p r i n c i p l e , three main advantages. F i r s t , the technique i s capable of handling a very general control problem including problems with bounded state and/or bounded control variables. Second, the solution i s good for any set of allowable i n i t i a l conditions, and; t h i r d , as the optimal controls are known at a l l points i n state space, the solution i s useful for r e a l -time optimal control. However, for a l l but the simplest cases, the computer memory that i s required to store the complete solution i s p r o h i b i t i v e l y large j^9j. This severe r e s t r i c t i o n i s what Bellman c a l l s "the curse of dimensionality". Some techniques have been recently developed that s i g n i f i c a n t l y reduce this problem of dimensionality; however, the dynamic programming approach i s s t i l l l i m i t e d to r e l a t i v e l y simple problems Jioj . Due to this l i m i t a t i o n , the emphasis i n this thesis w i l l be placed on the use of the direct and the indirect 4V approaches with an aim of developing techniques that reduce some of the d i f f i c u l t i e s currently experienced with these methods. 1.4 The Proposed Techniques The optimization techniques developed i n t h i s thesis are primarily a.combination of the direct and i n d i r e c t approaches. Two classes of problems are studied, and since the results are e s s e n t i a l l y unique, the material i s presented i n two parts. In Part I, a method i s developed for determining the optimal control laws for a class of aerodynamical systems of which the dynamics are l i n e a r i n the thrust and nonlinear i n the l i f t and thrust angle. Due to the presence of the l i n e a r con-t r o l , a variable thrust or singular subarc exists along which the maximum prin c i p l e cannot be applied. To overcome the d i f -f i c u l t y associated with the singular control, a method i s devel-oped to eliminate the unknown Lagrange m u l t i p l i e r s along t h i s variable thrust subarc. Conditions under which t h i s elimination i s possible are derived and expressions for thrust and rate of change of l i f t and thrust angle are obtained which depend only on state variables and a small number of scalar time-invariant parameters. The unknown parameters are then determined by a' search i n parameter space, based on a direct approach, for the set which minimizes the performance function. The proposed method i s considerably simpler than the standard numerical techniques which require a separate search i n function space for each component of the control vector. Furthermore, since the control vector i s generated by the direct solution of d i f f e r e n t i a l equations, the method appears suitable for i n - f l i g h t 5. guidance computers; that i s , either the control law i s obtained as a feedback law involving state variables only, or i t can be generated from state variables and a small number of scalar parameters. Several numerical examples are given i l l u s t r a t i n g t h is technique. In Part I I , a numerical algorithm i s developed to solve optimal control problems which do not contain singular control. The technique i s based on replacing a gradient search i n func-t i o n space by an equivalent search i n parameter space through the use of the necessary conditions of the indirect approach. To accomplish t h i s , the control law equations of the calculus of variations are used to eliminate the control vector from the system and the adjoint*!equations. However, instead of trying to solve the r e s u l t i n g two point boundary-value problem, the aug-augmented performance function i s considered to be a function of a l l unknown i n i t i a l conditions and i s minimized by a gradient search i n the parameter space of these i n i t i a l conditions. It i s shown that i t i s s u f f i c i e n t to search over any sphere for the intersection of the l i n e u.\Q, where \ Q i s the c l a s s i c a l solution of i n i t i a l values. As a r e s u l t , the proposed method i s not dependent on a good i n i t i a l estimate of the optimal trajectory. However, since the technique i s based on a gradient search, the f i n a l convergence slows down as the optimum i s approached. To provide improved f i n a l convergence, three techniques are devel-oped. The f i r s t i s based on the method of matching end points which uses an optimal scale factor for the Lagrange m u l t i p l i e r s such that the error i n transversality i s a minimum at each step 6. i n the i t e r a t i o n . The r e s u l t i n g algorithm i s independent of the i n i t i a l scale factor f or the Lagrange m u l t i p l i e r s and, hence, can be conveniently used with the proposed gradient technique. The other two method are based on determining the optimal step size for the gradient technique once i n the v i c i n i t y of the optimum. One approach uses a second v a r i a t i o n of the augmented performance function, and the other uses a method of curve f i t t i n g . It i s shown that by combining these techniques, a three stage algorithm can be developed which has good i n i t i a l convergence, good f i n a l convergence, and which requires storage at terminal points only. It i s also shown that a s i m i l a r approach can be used to improve the f i n a l convergence properties of the gradient search i n function space. Several examples are presented, consisting of both bounded and unbounded control. PART I OPTIMAL CONTROL LAWS FOR A CLASS AERODYNAMICAL SYSTEMS 7. 2. GENERAL ANALYSIS 2.1 'Optimal Control of Aerodynamical Systems Miele has given a general v a r i a t i o n a l theory for o p t i -mal f l i g h t paths of aerodynamical systems and derived general expressions for the Euler-Lagrange equations [ l l | . In special cases, these equations are useful i n deriving a n a l y t i c a l results concerning the nature of optimal f l i g h t paths Jl2j . However, with the exception of these few p a r t i c u l a r cases, numerical techniques are required to solve optimization problems and, as mentioned e a r l i e r , these techniques usually require the Storage capability of a large-size general purpose, d i g i t a l computer j l3-14j • However, i t has been shown that i n the special case of the sounding rocket, the Euler-Lagrange.equations can be used to solve the synthesis problem and the optimal thrust i s expressed i n closed form as a function of state variables only [12]. A more complex and more interesting problem i s the general case of a missile moving within the earth's atmosphere under the control of thrust, thrust angle, and l i f t such that a specified performance function i s a minimum. This i s an example of multi-dimensional control, and the standard numerical techniques are time consuming to apply since a separate search i n function space must be performed for each component of the control vec-tor. I t i s the purpose of Part I of t h i s thesis to extend the techniques, used to solve the sounding..rocket problem, to the case of multi-dimensional control. 2.2 Problem Statement The class of f l i g h t systems to be discussed are those 8. which can be represented by the following state vector differen-t i a l equation: ( 2 . 1 ) x = G Q(X . L) + uG x(x,p) where x = X = n - vector of state variables dx x = TJTT denotes the time derivative of x G-Q(X.L) = n - vector of functions gQ^ . of x and L G-^(x.p) = n - vector of functions g^^ of x and p and where the dynamics are l i n e a r i n the thrust u and non-l i n e a r i n the l i f t L and thrust angle |3. (See Figure 2 . 1 . ) -> X Figure 2 . 1 The Motion of a Rocket Within the Earth's Atmosphere 9. The derivation of (2.1) for the case of a missile moving i n the earth's atmosphere i s given i n Appendix A. Consider now the problem of determining the set of controls (u, L, |3) which takes the system (2.1) from some i n i t i a l manifold defined by the k-vector constraint H ( x ( t Q ) , t Q ) = 0 (2.2) to some f i n a l manifold defined by the p-vector constraint G ( x ( t f ) , t f ) = 0 (2.3) where u i s subject to the constraint u ( umax. " u ) " a" = 0 ( 2 ' 4 ) and where the performance function P = P ( x ( t f ) , t f ) (2.5) i s to be minimized.. The constraints (2.2) and (2.3) are vec-tor equations whose dimensions are equal to or less than n. The stated problem i s of the Mayer type and can be solved by introducing the augmented function [ l l j P k X T[x - G Q - uG^| + X n + 1 £ ( u m a X o - u) - a*] (2.6) T Here \ i s an n-vector of Lagrange m u l t i p l i e r s , \ i s the trans-pose of X a n d ^ n + 1 i s the (n+l) - th Lagrange m u l t i p l i e r . The Euler-Lagrange equations obtained from (2.6) are \ = - (2.7) 0 = k L (2.8) 0 = uk p (2.9) 0 = k u + K+l ( 2 u " V u J ( 2' 1 0> where . 0 = a ? w L 3L A 3G p 3(3 A T k = G-A u 1 10. (2.11) (2.12) (2.13) (2.14) and where the following abbreviated notation i s used: 3 g 0 1 . . 3 x 0 A 3 x n 3 g 0 1 5 x — n 3 f 0 n , . . . § f O n 3 x - , 3 X 3 G 0 A 3L~ = n (2.15) (2.16) Since GQ and G ^ are formally independent of t , a f i r s t i n t e g r a l ( x ) T \ = c (2.17) ex i s t s , where c i s a constant of integration. Substituting (2.1) into (2.17) yields (2.18) The transversality condition for the stated problem i s it. T T G Q \ + uG x \ = c [dP + VTdG + (dx) T\ - c d t l " f = 0 J t Q (2.19) o 11. where dx and d t must be c o n s i s t e n t w i t h t he t e r m i n a l c o n s t r a i n t s (2.2) and (2.3) and where V i s a p - v e c t o r o f c o n s t a n t Lagrange m u l t i p l i e r s w h i c h a re i n t r o d u c e d t o accoun t f o r (2 .3) . 2.3 D i s c u s s i o n o f t he Necessary C o n d i t i o n s f o r an O p t i m a l T r a j e c t o r y Because the sys tem (2.1) i s l i n e a r i n u , and because o f t h e c o n s t r a i n t (2 .4) , t he o p t i m a l t r a j e c t o r y will i n g e n e r a l c o n s i s t o f maximum t h r u s t subarcs where u = , v a r i a b l e t h r u s t subarcs where O^u-^^ax » a n d minimum thrust or c o a s t i n g subaros where u = 0. The variable thrust subarcs are also known as singular subarcs slnoe the Hamiltonian is then independent o f u and, consequently, the maximum prinoiple oannot be applied to determine u. It follows from ( 2 . 4 ) that a a 0 along maximum and minimum thrust eubaros and henoe (2.11) is satisfied. Along the variable thrust eubaro i t follows from (2.11) that \ ^ •'0 and henoe (2.10) yields \ m 0 (2.20) The method to be dieoussed requires a knowledge of the sequenoe of subaros. This oan be determined from the Legendre-Olebsoh oondition and the Erdmann-Weierstrass corner oonditions. The Legendre-Clebsch oondition applied to (2.6) requires that § £ | ( SL ) 2 + ( Sp ) 2 + ^ | ( Su ) 2 + 1 2 ! ( S « ) 2 * 0 3L 3p 3 u ' 3 t x (2.21) where S u and S a a re r e l a t e d by (2 .4) . S u b s t i t u t i n g (2.6) and (2.10) i n t o (2.21) y i e l d s GOLL X ( § L ) 2 ( S P ) 2 + u 2u-u, max. 12. [(Su)2+(S:a)2]^ 0 (2.22) where the following abbreviated notation i s used: OLL, ( 2 . 2 3 ) ' As §L, S(3 and Su are independent, i t follows from (2.22) that GOLL A = 0 (2.24) (2.25) everywhere along the optimal trajectory and that k > 0 , (u=u ) u ' max. k = 0, (0<u<u ) u ' ^ max. (2.26) k u<0, (u=0) Due to the Erdmann-Weierstrass corner conditions, the Lagrange m u l t i p l i e r s and the f i r s t i n t e g r a l (2.18) are continuous along the optimal trajectory. Substituting (2.14) into (2.18) ;yields GQ\ + ;uk u ( 2 . 2 7 ) I t follows from (2.27)- and-the Erdmann-Weierstrass corner conditions that a discontinuity i n u i s possible provided that (k ) = (k ) = 0 u - u' + (2.28) where the minus and plus subscripts denote evaluation of the 13. brackets just before and just after each switching point res-pectively. Hence, k u i s a continuous function which vanishes at each switching point (see (2.26)). Tor this reason, k^ w i l l be defined as the switching function. The switching func-t i o n i s a fundamental importance i n determining the allowable sequence of subarcs. To investigate the properties of k u, (2.14) i s diff e r e n t i a t e d with respect to time. Using (2.1) and (2.7) to eliminate x and X yields k = u T T T T G,. G-, — G, G^ 0 l x . 1 Ox X + u „ T T o,To, T G-. Gn - G-, G, 1 l x 1 l x + Pk p (2.29) It i s seen that the coeff i c i e n t of u i n (2.29) i s i d e n t i c a l l y zero. Furthermore, the term pk^ i s zero, since either u^ O and kp=0 (see (2.9)) or i f u=0, the thrust angle can be defined to be constant so that (3=0. Hence (2.29) yields ku.= z\ (2.30) where K = G l x G0 " G0x G l ( 2-5l) It follows from (2.30) and the Erdmann-Weierstrass corner con-ditions that at the corners of the optimal trajectory ( k u ) _ = ( i u ) + (2.32) The condition (2.32) can be used to determine.the possible sequences of subarcs. Expanding k^ i n a Taylor series i n t at the switching instant t yields s k (t + At) = k (t ) + k (t )At + (2.33) u s u s u s 14. Substituting (2.28) and (2.33) into the condition (2.26), and choosing At so that second and higher order terms are ne g l i g i b l e , yields the conditions (2.34) At = 0 which apply for maximum thrust, variable thrust, and zero thrust respectively. Table 2.1 u(t - At) s u ( t s + At) (1) u = O k u < ° (2) u = u k u > o (la) (lb) K TX>0 Kl\ = 0 max. (3) u = u(t) k = ,0 u (2a) Z X\<0 (2b) K T\ = 0 (3a) K T\ = 0 (3b) K T\ = 0 ( l a> u = umax,' ku>° (lb) u = u ( t ) , k u = 0 (2a) u = 0, k u<0 (2b) u = u ( t ) , k u = 0 (3a) u = 0 (3b) u = u .. k u < ° 'max. l> i f u ( t s ) = u max. Table 2.1 i l l u s t r a t e s the possible sequences of sub-arcs which s a t i s f y the Legendre-Olebsch condition. The sign of k defines the state of u(see (2.26)). The instants of time u where k vanishes defines the switching points t . The sign u. s T ' of K \ at t = t determines the state of u after switching. s The symbolism (la),....,(3b), used to denote these states w i l l be' made use of l a t e r i n obtaining a sequence diagram. It w i l l be shown by means of examples that an a n a l y t i c a l study of the system and i t s constraints can often provide the necessary 1 5 . information to determine a unique sequence of subarcs by use of ( 2 . 2 8 ) and Table 2 . 1 . 2.4 The Control Equations for Thrust. Thrust Angle and L i f t To determine suitable control equations, i t i s , desirable to eliminate the Lagrange m u l t i p l i e r s and obtain equations r e l a t i n g the state variables and u, p and L only. I f these equations do not involve time derivatives of u, 0 and L, a feedback control law i s obtained. However, th i s may not always be possible. I t i s the purpose of t h i s section to determine con-ditions under which a feedback law can be d i r e c t l y obtained and to develop suitable alternatives when these conditions are not s a t i s f i e d . The elimination of the Lagrange m u l t i p l i e r s i s advantageous since the magnitudes of the control variables are usually known approximately while the magnitude of the Lagrange m u l t i p l i e r s are unknown. To eliminate the Lagrange m u l t i p l i e r s , ( 2 . 8 ) and ( 2 . 9 ) are f i r s t d i f f erentiated with respect to time. Eliminating x and \ by means of ( 2 . 1 ) and ( 2 . 7 ) yields L = L Q + u L ± ( 2 . 3 5 ) and p = p Q + u p i ( 2 . 3 6 ) where T A ( G Q L G0x " G Q T G Q L X ) X ' (0 , 7 N LQ = : rji a { 2 . 3 / ) & O L L X L 1 = ^ 1 ^ 1 ^ 4 1 1 ( 2 . 3 8 ) G O L L X a A ( G13 G0x " GQ G 1 8 x U G l p p X a A ( G13 G l x - G i T G 1 8 x ) X G18S X 16. (2.39). (2.40) and where the following abbreviated notation i s used: >01 . . . . G, OLx 3 Sr 3 2 g o i n d 2gpn -2 d S 0n 3Ld" x n (2.41) I f the inequality sign holds i n (2.21), i t follows from (2.24) T T and (2.25) that the d i v i s i o n by \ and, G-^ p X i s permissible. Along the variable thrust subarc k u i s i d e n t i c a l l y zero (see 2.20). Thus the time derivatives of are also zero. I t follows from (2.30) that during variable thrust k u = Kl\ = 0 u • rp rp • k 1 X + KX \ = 0 (2.42) (2.43) Using (2.31) to evaluate K and eliminating x, \, L and S by means of (2.1), (2.7), (2.35) and (2.36) yields u = if G0x 0 0 i ^ r n r t i m m m m r n r P Gi, K__ - KG-, _ + K/ + L, K^1 (2.44) 1 x l x T ^1 p 1 If \ can be eliminated from (2.35), (2.36) and (2.44), • « equations are obtained for u, L and p i n terms of state variables 17. only. These equations can then be used to investigate the p o s s i b i l i t y of optimal and sub-optimal feedback laws which hold during the variable thrust subarc. The elimination of \ i s pos-s i b l e i f a s u f f i c i e n t number of l i n e a r independent equations between the components of X can be found. There are f i v e such l i n e a r equations given by (2.8), (2.9), (2.18), (2.20) and (2.4-2). I t may be possible to augment these f i v e equations by further l i n e a r equations obtained by a direct integration of the Euler-Lagrange equations. Let m be the t o t a l number of such l i n e a r equations. These equations can be written i n the matrix form A\ = b (2.45) where A i s a mxn matrix, where c "1 c l b = c m (2.46) and where the c^ . are integration constants. I f b / 0, (2.45) can be used to determine X uniquely ; i n terms of the pro-vided that rank (A) = n = m (2.47) (It should be noted that for the case m>n, the integration constants cannot be independently specified.) If b = 0, (2.45) i s a set of l i n e a r homogeneous equations i n X and i f rank (A) = n - 1 =. m (2.48) 18. i t i s possible to obtain a unique n o n - t r i v i a l solution of (2.45) i n the form X± = Q i \ > ( i ?^ k ) ( 2 * 4 9 ) so that a l l \^ (i^k) are expressed i n terms of one Lagrange m u l t i p l i e r Introducing (2.49) into (2.37), (2.38), (2.39), (2.40) and (2.44) then yields the control equations of u, L and p i n terms of state variables alone. Conditions (2.47) and (2.48) therefore serve as a test to see i f i t i s possible to eliminate the Lagrange m u l t i p l i e r s during a variable thrust subarc. Several examples w i l l now be given to i l l u s t r a t e var-ious p o s s i b i l i t i e s which can occur. 3. THE SOUNDING ROCKET PROBLEM 19. 3.1 Derivation of Optimal Control Law The derivation of the feedback control law for the variable thrust subarc has been given i n |l2 . However, the use of the theory outlined i n Chapter 2 can be used to prove the existence of feedback control laws for more general cases and, also, to provide a more systematic means for obtaining analytic expressions for control laws which are a function of the state variables only. The terminal conditions for the sounding rocket problem are t Q = 0 y (o) = 0 v (o ) = v n 0 ( 3 - D m(0) = mQ v ( t f ) = 0 m(t f) = mf and the f i n a l a l t i t u d e y^ i s to be maximized for a given amount of f u e l . The performance function P = -y^ i s to be minimized. Appendix B gives the analysis associated with this problem. The transversality condition ( B - l l ) yields c = 0 , \ 2 f = 1 (3.2) It follows that b = 0 (see (2.46) and (B-13)), and since n = m = 3, condition (2.48) requires that rank (A) = 2 i f a feedback control law i s to exist. To determine rank (A), (B-12) can 20. be triangularized y i e l d i n g the matrix where B = 1 0 0 £ _ Pi. v mv 1 0 0 SL v. f. (3.3) mg - (D + 2L- D) e (3.4) and where rank (A) = rank (B). It i s seen from (3«3) that the condition (2.48) i s s a t i s f i e d along the variable thrust subarc i f (3.5) Substituting the f i r s t i n t e g r a l (B-10) into (B-9) yields the time-derivative of the switching function. k = u mv X3 ' uk 2 f u r n s (3.6) During a variable thrust subarc, i t i s seen from (3.6), (3«5) and (2.26) that k = 0:as required. The function f , which n s i s a function of state variables only, can be used to determine the optimal feedback control law. The time derivative of f i s f = -uM + N s (3.7) where 21. M A g + £ ( 3 + ^ ) (3.8) N k (g + + 22) + a D v ( 1 + Z - ) ( 3 . 9 ) e e (For the d e f i n i t i o n of a see (B-4).) It follows from (3.5) that f = 0 along a variable thrust subarc. Equation (3*7) can then be solved for the optimal control law. u = | (3.10) To determine the control law for the complete t r a -jectory requires additional information about the sequence of subarcs. Consider f i r s t the case where the maximum thrust subarc i s one of impulsive boosting (u„__ = co). Since M>0, i t follows from (3.7) that f < 0 along the maximum thrust s subarc. The sequence of possible subarcs can be obtained with the aid of the sequence diagram i l l u s t r a t e d i n Figure 3.1 which i s a graphical representation of the f i v e possible states associated with the signs of k u, k and condition f =0. Eight gates, ( l a ) , (lb) to (4a), (4b) are provided s to indicate the allowable change i n state. The position of gates ( l a ) , (lb) to ( 3 a ) , (3b) for a l l problems of the type discussed i n Section 2.1 are determined by the use of Table 2.1. From (2.26) and (2.42), i t i s seen that regions 1 .:. and I I i n the sequence diagram are regions of maximum thrust, regions'III.and. IV: are regions.of zero thrust, and'region V i s a region of variable- thrust where Figure 3.1 Sequence Diagram for the Sounding Rocket with Impulsive Boosting 23. u i s given by ( 3 - 1 0 ) . The positions of gates (4a) and (4b) are not given by Table 2.1 and must be determined by further analysis of each par t i c u l a r problem. Consider, for example, the case of the sounding rocket. To determine the position of gate (4b) the sign of k must be determined at the instant where :k ••= 0. ° u u * * D i f f e r e n t i a t i n g (3-6) with respect to time and evaluating k when k = 0 and u = 0 yields \,v k u = - - i - s f s (3.1D m v From (3.7) i t i s seen that during coasting (u = 0) ta = N> 0 (3.12) In Appendix B i t i s shown that 0 along the i n t e r i o r of the optimal trajectory. I t then follows from (3.12) and (3.11) that k u< 0 when k = 0 and gate (4b) must therefore open down as shown i n Figure 3 . 1 . The determination of the position of gate (4a) proceeds i n a s i m i l a r manner. D i f f e r e n t i a t i n g (3»£.) with respect to time and using u = u yields u 2 m v p, 2v f D-, D ( _ e v _ } _ N _ B _ -max. m v v ' mv e e (3.13) I t follows from (3.13) and 0 that during maximum thrust, the condition for k^> 0 i s f D s mv % a x . > D ' 2v * , (5,14) e In t h i s example, i t i s assumed that u —** 0 0 , and hence (3.14) 24. i s s a t i s f i e d . Gate (4a) must therefore be directed up as shown. The sequence diagram for the sounding rocket i s now complete and the optimal sequence for any set of terminal con-st r a i n t s can be determined. It i s proven i n Appendix B that for the end constraints ( 3 . 1 ) the f i n a l subarc must be a coasting subarc. Also, from ( 3 . 7 ) , i t i s seen that for the cases of a coasting, a variable thrust, and a maximum thrust subarc respectively. Combining th i s information with the sequence diagram permits the evaluation of an acceptable sequence of subarcs. Several p o s s i b i l i t i e s are i l l u s t r a t e d i n Figure 3 . 1 . The type of i n i t i a l subarc i s determined from the i n i t i a l value of f as shown. This i n i t i a l subarc i s maintained s u n t i l either f vanishes or u n t i l a l l f u e l i s consumed. The s vanishing of f indicates a switch to the variable-thrust sub-to s arc where u i s programmed according to ( 3 . 1 0 ) , and the instant of burn-out indicates a switch to the f i n a l coasting subarc. It i s seen, therefore, that ( 3 . 4 ) and (3.10) completely define the optimal control as a function of state. Also, f or any system s a t i s f y i n g condition ( 3 . 1 4 ) , there can be at most three subarcs i n the optimal trajectory. If (3.14) i s not s a t i s f i e d at the instant when k vanishes, u ' then i t follows from (3.13) that gate (4a) must be changed to the down position as i l l u s t r a t e d i n the sequence diagram shown i n Figure 3«2. Assuming that t h i s be the case, i t i s seen that a closed loop i n the sequence diagram can exist whose (3.15) The case of u max. < ° o i s of p a r t i c u l a r i n t e r e s t . 25. Figure 3.2 Sequence Diagram for the Sounding Rocket with F i n i t e Maximum Thrust 26. sequence i s ..... u = u m & x > , ° < u < u m a x . » u = u m a x . ' «(case l a ) . This closed loop presents no a n a l y t i c a l d i f f i c u l t y as switching to a maximum thrust subarc i s defined when u given by (3«10) equals u , and switching back to a variable thrust max. subarc occurs when f vanishes. However, d i f f i c u l t y does arise i f the i n i t i a l subarc i s a coasting subarc. For t h i s case, the f i r s t switching may be to a maximum thrust subarc (case I-IIb) or to a variable thrust subarc (case I l i a ) , depending on the 9 value of k^ when k u = 0 . To avoid the use of k , which con-tains unknown Lagrange m u l t i p l i e r s , the f i r s t switching instant t , where k vanishes, can be introduced as an unknown parameter. S "U. Let T be the instant where f vanishes which determines the s s switching instant to a variable thrust s u b a r c .If t <T the s s control switches to. ..a maximum thrust:, sub arc and i f t = TT the s s control switches to a variable thrust subarc. To determine t , s the performance function P can be considered to be a function of t and a search over the i n t e r v a l 0 = t ^ % can be performed s s s to determine the minimum of P. A numerical example of this type of search i s given i n [12] . * o = 0 x(0) = 0 y ( o ) = 0 v(0) = v , 4. THE MAXIMUM RANGE PROBLEM 4.1 Derivation of Optimal Control Laws for the Case L = 0. "P = Q Consider the problem of maximizing the range for a given amount of f u e l when the terminal constraints are ©(0) = © 0 m(0) = m0 (4.1) y ( t f ) = o 0 m(t f) = mf where the control constraints 6 = 0, L = 0 are applied, and where the performance function P = -x f (4.2) i s to be minimized. Appendix C gives the analysis associated with t h i s problem. It i s assumed that the i n i t i a l conditions for the rocket are obtained by means of a launching platform. Also, i t i s assumed that u = CO Is a good approximation to I H £ I X a the maximum thrust subarc. These assumptions do not impose any r e s t r i c t i o n s on the control laws derived for the variable thrust subarc. However, they do simplify the numerical compu-tation and attention can thus be focused on the derivation of the control laws and t h e i r application i n determining optimal t r a j e c t o r i e s . The f i r s t step, i n obtaining these control laws, i s to evaluate the transversality condition (C-14) according to (4.1) and (4.2) which yields \ 3 f =0 \ l f = C l = 1 \ 4 f = 0 c = 0 (4.3) 28. Taking c = 0 i n the f i r s t i n t e g r a l (C-13) and substi-tuting into (C-12) yields v \ , f 2\. k = __& ( 2_s 1 g cos 9 + uk ) (4.4) u mv m v ° u ' where f o = mg s i n 9 - D(l+v/vJ (4.5) Substituting (4.3) into (C-16) i t i s seen that b ^ 0. As (C-15) i s a 4 x 5 matrix where m = 4, n = 5, the condition (2.47) i s violated. Hence, a feedback control law, i n terms of state variables only, cannot be obtained d i r e c t l y and a modified approach must be adopted. During a variable thrust subarc, equation (2.26) must be s a t i s f i e d . Substituting ( C - l l ) into (2.26) yields k = - \_ = 0 (4.6) u m 5 Substituting (4.6) into (C-33) and equating k to zero gives 2\ g cot 9 + \,N * = — 0 5 ^ (4.7) 3 where M and N are given by (C-34) and (C-35) respectively. As \j £ o along a variable thrust subarc (see Appendix C), d i v i s i o n by \^ i s permissible. With the aid of (4.6), therefore, i t i s possible to write (4.7) i n the form 2ov o g cot 9/m + U u = - M where A •e-— 1 (4.8) V= (4.9) D i f f e r e n t i a t i n g (4.9) with respect to time and using (4.8) and 29. (C-10) yields X = — - u) (4.10) e which i s v a l i d for a l l values of Equations (4.8) and (4.10) give the optimal control law for th i s problem.. As a result of (4.10), an unknown parameter X i s introduced which i s the i n i t i a l value of $ at the beginning of the variable thrust subarc. One other unknown parameter e x i s t s , however, and that i s the instant of switching to the variable thrust subarc. s This instant i s characterized by. the vanishing of the switching function k . However, using k to define involves the use of u ~' ° u s the unknown Lagrange m u l t i p l i e r s . To avoid determining the Lagrange m u l t i p l i e r s , an alternate approach i s adopted. The unknown switching instant T l i s replaced by an unknown switching v e l o c i t y v . Consequently, the unknown parameters for this problem are X g and v g. The proposed technique, for obtaining the value of any unknown parameter a^ ., i s to consider the per-formance function as a function of thex parameters and then solve the minimization problem Min f Pi ^ (4.11) k by a direct search i n parameter space. Using (C-3) and a sequence diagram, i t can be shown that acceptable sequences for th i s problem are given by < 1 } v 0 < v s ' u = umax. ' °< u< umax. ' u = ° (2) v Q > v s , u = 0 : 0 < u < u m a x o , u = 0 (3) v Q = v s , 0 < u < u m a x > , u = 0 / 30. In case ( l ) , v n < v and hence the i n i t i a l subarc must be a maxi-mum thrust subarc along which v n increases to v . In case (2) v A > v and hence the i n i t i a l subarc must be a coasting subarc 0 s ° along which v n decreases to v„,. When v~. = v , the thrust ° 0 s 0 s switches to variable thrust u n t i l a l l the f u e l i s consumed. As an i l l u s t r a t i v e numerical example, consider the following data IIIQ = 35.0 slugs m^ = 10.0 slugs «o - 4 5 ° v Q = 1000 ft/sec (4.12) v = 5500 ft/sec a = (22000 f t ) " 1 K = 10" 4 slug/ft (see (B-4)) In t h i s example, VQ i s a r e l a t i v e l y small i n i t i a l v e l o c i t y and hence the sequence associated with case ( l ) r e s u l t s . As u ^ ^ = O o i s assumed, (C-l) to (C-5) can be integrated over the maximum thrust subarc using (4.1) from v = v n to v = v to y i e l d x s = 0 y s = o 9 S = 9 Q (4.13) ms = m 0exp((v 0-v s)/v e) The integration along the variable thrust subarc i s now per-formed using (4.13) as i n i t i a l values. The unknowns are v and 31. as m can be obtained from v through (4.13). The performance S S S function to be minimized i s P = -x„(v , ^ ). The computation i s s was performed on an IBM 7040 d i g i t a l computer using the following algorithm: 1» Select a set of values (v , X ) . s s 2. Solve (4.13) to obtain m . s 3. With u and X given by (4.8) and (4.10), integrate the system equations (C-l) to (C-5) and determine x f ( v s J s ) . 4. Repeat 1, 2, 3 "and carry out a direct search i n parameter space (v , ^ g) for the maximum x^. The values for the state variables at the switching points along the optimal trajectory are shown i n Table 4.1. Table 4.1 I n i t i a l Point End of Max. Thrust Burnout P i n a l Point t(sec.) 0 0+ 26.53 288.1 x(ft.) 0 0 90560 1470000 y(ft.) 0 0 74870 0 v(ft/sec.) 1000 3045 6670 4970 9(deg.) 45 45 36.6 -41.0 m(slug) 35 24.13 10.0 10.0 Figure 4.1 i l l u s t r a t e s x« as a function of v o for the optimal value of 6 , and x„ as a function of $ for the optimal value s i s of v . The graphs of u(t) and $(t) f or the optimal trajectory are shown i n Figure 4.2, 32. 1470020 14700Q0Y 1469960 1469920V 1469880V 1469840V \1469800 3025-0 3045-0 3 065-0 Vs (FEET/SEC) «/> 0-000005 0 00001 0-000015 0-0000 Figure 4.1 The F i n a l Range as a Function of V and "6, o CO \ CO CD =3 ~ J CO 0-55 0-54 0-53 0-52 0-50 33 . $ 10 12 14 16 TIME (SECONDS) 18 20 22 24 26\ 28 30 26-53 o-ooooi II 0-000005 0-0000000 c_> § rco CD CO o o § 03 5 Q 8 10 12 14 16 TIME (SECONDS) 18 20 22 24 26j 28 i X '53 Figure 4.2 u(t) and'u'(t) for the Optimal Trajectory JO 34. 4.2 Derivation of Optimal Control Law for the Case L = 0 The zero l i f t case where u and 3 are the control variables i s a two dimensional control problem. Appendix D gives the analysis associated with t h i s case. I f P i s inde-pendent of time, the transversality condition (D-15) yields c = 0. The class of problems for which c = 0, ^ 0, w i l l be discussed i n this section as the maximum range problem i s a pa r t i c u l a r example. For th i s class of problems, i t follows that b ^ 0 (see D-18). Furthermore, i f the f i r s t row of matrix A given by (D-16) i s interchanged successively with the second and then with the t h i r d row of A, the re s u l t i n g matrix can be triangularized y i e l d i n g the matrix B 1 0 0 0 0 0 1 0 0 0 0 g s i n 9 + D/m v s i n 9 1 0 0 0 g cot 9 2 v cot 8 V 1 0 0 0 0 mv s i n 0 v (4.14) which has the same rank as A and where f = mg sin 2(9+8) - D s i n 9(1 + + D c o s Q. P Vo •sin 8 cos 8 (4.15) For the case where f^ ^ 0, (2.47) i s s a t i s f i e d (rank (A) = rank (B) = 5 = m). It then follows that the Lagrange m u l t i p l i e r s can be eliminated from (2.39), (2.40), and (2.44), y i e l d i n g equations for u and 8 i n terms of state variables only. Substituting (D-16) and (D-18) into (2.45) yields a set of 35. l i n e a r non-homogeneous equations i n A... This set of equations can be used to express a l l Lagrange m u l t i p l i e r s i n terms of c-^ for the case that c^ ^ 0. However, to obtain a control law which i s more generally applicable, and which includes the special case c-^ = 0 for which b = 0, i t i s more convenient to express a l l Lagrange mu l t i p l i e r s i n terms of Conse-quently, a set of equations of the form (2.49) i s obtained. The triangularized form given by (4.14) i s useful i n obtaining these relationships. Substituting (D-19), (D-25) and the relations (2.49) into (2.36), (2.39), and (2.40) yields 8 = mv 2mg cos 9 - cot 3 <(D(l + ^— cos 3) - mg s i n 9^ > - v u s i n 3 e r (4.16) Comparing (2.30) and (D-12) yields K = — m s i n 9 s i n 3 - cos 9 cos 3 -cos 3 s i n 9 - s i n 3 cos 9 2D mv c o s B + ^ + f s i n 3 cos 9 e _ g_ The matrices K = x 2 cos 9 cos 3 + —2 s ^ n 3 v mv 0 3 k., (4.17) 0 0 9 0 9 0 8 k E dx x 9 k . rjx^ (4.18) 36. and K 6 = dk2 3T (4.19) 3k are now evaluated where k^, ..... k^ are the elements of K. Substituting (4-17), (4.18), (4.19), and (D-19) into (2.44), eliminating the Lagrange m u l t i p l i e r s by use Of (2.49), and per-forming a series of algebraic manipulations yields u = N/M (4.20) where M = D 2 v ( c o s 8 + c 7 i ~ p ) + v e ( 2 + s i n 2 p ) + v7 r D ( l + v ° 0 S 3) - mg s i n 9 (4.21) 2 N = - cos B si h 8 mg cos 9 cos 8 + mD(l + — cos 8) e g(2 cos 9 s i n 6 + 3 s i n 9) av s i n 9 + cos 8 + . 2 D m cos 8 D 2 + + (g - av ) cos 9 s i n 8 e 2g s i n 8 (v 9 s i n 8 - s i n 9) cos 8 v K r e (4.22) Equations (4.16) and (4.20) are the desired control laws which are v a l i d for a l l problems where c = 0. An i l l u s t r a t i v e numerical example of (u, 0) control w i l l now be presented. To allow a comparison with the case 37. where 3 was constrained to be zero, the maximum range problem of Section 4.1 w i l l be solved, with the addition of 8 control. For the maximum thrust subarc, equations (D-34), (D-37), and (D-40) apply for v n ^ v ^ v . (The subscript s denotes the \J S values of the variables at the instant which terminates the s maximum thrust subarc.) 'The unknown parameters i n t h i s problem are 3^ and v . Note that although there has been an increase 0^ s to i n control dimension, there has been no increase i n the number of unknown parameters. The computation was carried out on an IBM 7040 d i g i t a l computer using an algorithm similar to that i n Section 4.1. The values for the state variables at the switching points along the optimal trajectory are given i n Table 4.2. Table 4.2 I n i t i a l Point End of Max. Thrust Burnout Pi n a l Point t(sec.) 0 0 + 30.32 331.9 x(ft.) 0 0 80533 1545400 y(ft.) 0 0 83440 0 v(ft./sec.) 1000 2357.4 6727.0 5252 9(deg,) 45 45 43.1 -46.95 m(slug) 35 27.29 10.0 10.0 T Figure 4.3 i l l u s t r a t e s x„ as a function of v for the optimal value of 3Q> a n (^ x f a s a function of p A for the optimal value of v„. The optimal controls r e s u l t i n g from these parameters are shown i n Figure 4.4. For purposes of comparison, a ba l -l i s t i c trajectory from the same i n i t i a l conditions i s presented for which a l l the f u e l i s consumed during the boosting subarc. .1 38. 1525 00 0^ ' 1 — I 1—4 I 1 U 1 L_ 2335-0 234 5-0 2355-0 2365-0 0-18 0-19 0'20 0-21 o o Vs (FEET/SEC) (l0'3) @o (RADIANS) ('2'0) (DEGREES) Figure 4.3 The F i n a l Range as a Function of v • and 39. U MAX BOOS TING # U MAX (NOT TO SCALE) VARIABLE THRUST SUBARC-4 8 12 16 20 24 TIME (SECONDS)- VARIABLE THRUST SUBARC Figure 4.4 The Optimal Controls u and 6 for the Maximum Range Problem 40, The graphs of the optimal t r a j e c t o r i e s for the u control, (u, 8) control, and the b a l l i s t i c trajectory are shown i n Figure 4.5. The f i n a l ranges for these three cases are: 1. u = 0, 8 = 0 ( b a l l i s t i c ) , x f = 1,050,000 f t . 2 . 8 = 0 , x f = 1,470,000 f t . 3. M 0 , x f = 1,545,400 f t . - •', (4.23) The sequence, boundary,: shown in. Figure, 4.3 and subsequent figures, represents the locus of a l l points at which u = 0 at burnout. For points on the other.side, of th i s boundary, the variable thrust goes to zero before a l l the fuel i s con-sumed, and a sequence i s required of the form u = ~^ m a x > 0>$Cu<Cu , u = 0, u = u , O^u^Cu , u = 0, .... . ^ max.' ' max. ^ ^ max.' ' However, as a true extremum was obtained for the acceptable sequence u = u , 0<^ u><Cu , u = 0, other possible sequences were not investigated. 4.3 The Maximum Range Problem with Control of F i r i n g Angle Consider the maximum range problem of Section (4.1) when the terminal constraints are m(0) = 50.0 slug y ( t f ) = 0 m(t f) = 10.0 slug (4.24) In t h i s problem, the f i r i n g angle, 9(0) = QA, i s free and hence, the transversality condition (C-14) yields V 0 ) = ^40 = 0 ( 4 e 2 5 ) *0 = 0 x(0) = 0 y(o ) = 0 v(0) = £>0 5 / 2 © (U,fi) CONTROL 8 U CONTROL BALLISTIC TRAJECTORY 12B 256 384 512 64 0 768 096 1024 1152 1280 1408 1536 RANGE- (FEET x 10~3) Figure 4.5 The Optimal Trajectories for the u Control, the (u,8) Control, and the B a l l i s t i c Trajectory 4^ H 42. The i n i t i a l v e l o c i t y i s assumed to be zero. However, to keep o 0 f i n i t e during the boosting stage, the i n i t i a l v e l o c i t y i s taken equal to a small non-zero value e. For this case, as e—•O, the f i r s t subarc must be a maximum thrust subarc. It cannot be a coasting subarc for the system would remain at rest for v(0) =0. For a variable thrust subarc, k (see C - l l ) ' u and k u (see (4.4)) are both zero. Therefore, i f the conditions (4.24), (4.4). and (4.5) are evaluated at t = 0, and i f con-d i t i o n (4.25) i s used, i t i s seen that k u(0) = 0 only i f \^(0) = 0. However, i t i s seen i n Appendix C. that \^ ^ 0 along a variable thrust subarc. Consequently, the only remaining p o s s i b i l i t y i s that the f i r s t subarc be a maximum thrust subarc. It then follows from (4.25) and (C-9) that at the end of the maximum thrust subarc for u „ — 0 0 max. X 4 s * X 4 Q = 0 (4.26) Substituting (4.26) into (4.4) and evaluating k u at the start of the variable thrust subarc yields k ( r ) = - ^ s f s = 0 ( 4 o 2 7 ) u s in v s s Hence, as ^ / 0 along a variable thrust subarc (see Appendix C) f i s determined as the instant f vanishes, which from (4.5) s s yields ni g s i n 9 - D (l+v / v j = 0 (4.28) By use of a sequence diagram, an acceptable sequence i s found to be 4 3 . u = u , max.' 0<u<u , ^ ^ max.' u = 0 , o<t<-£ where f i s defined by ( 4 . 2 8 ) and f, i s the instant of burnout. The control laws for the problem are s t i l l given by ( 4 . 8 ) and ( 4 . 1 0 ) . Furthermore, i t i s interesting to note that the free f i r i n g angle does not result i n an increase i n the number of unknown parameters as ^ g and © g can be taken as the two unknown parameters with 9Q, v g , and mg being obtained from ( 4 . 2 8 ) and ( 4 . 1 3 ) . The performance function to be minimized i s then P = -x„(X , © ) and an algorithm s i m i l a r to Section ( 4 . 1 ) X s s i s used. The re s u l t i n g values for the state variables along the optimal trajectory are shown i n Table 4 . 3 . Table 4 . 3 I n i t i a l Point End of Max. Thrust Burnout F i n a l Point t(sec.) 0 0 + 3 0 . ' 0 3 5 3 7 8 . 8 x(ft.) 0 0 8 6 0 0 0 . 0 2 0 1 4 2 1 2 . 0 y ( f t . ) 0 0 9 2 0 0 0 . 0 - 0 v(ft,/sec.) e ^ 0 2 4 1 0 7 7 3 1 . 0 5 9 9 5 . 6 9(deg.) 5 3 . 2 5 3 . 2 4 3 . 9 - 4 7 . 1 m(slug) 5 0 . 0 3 2 . 2 6 1 0 . 0 1 0 . 0 The v a r i a t i o n of x^ i n the neighbourhood of the optimum i s shown i n Figure 4 . 6 for x„ as a function of 9 with o* optimum and i n Figure 4 . 7 for x~ as a function of with 9 optimum. The re s u l t i n g optimal control, u ( t ) , i s shown i n Figure 4 . 8 . 20U200 2014000 2013800 2013600 - 2013400 U. 2013200 2013000 0-9250 (53°) NT IT) o a? = 24/0 FT/ SEC Ys = 0-000008 0-S275 O-.9300 0S325 0-9350 (53-5°) 0-9375 0-9400 0-9425 0-9450 (54-f) 9S - RADIANS (DEGREES) Figure 4.6 The Fi n a l Range as a Function of 9 4^ Figure 4.7 The F i n a l Range as a Function of $ g TIME ( SECONDS) Figure 4.8 The Optimal Control u for the Maximum Range Problem with Free F i r i n g Angle 47. It i s interesting to note that the optimal f i r i n g angle for t h i s case i s about 53°; whereas, the optimal f i r i n g angle of a b a l l i s t i c missile i n a vacuum i s 45°• This same problem was also studied with the addition of thrust angle control s i m i l a r to the problem i n Section 4.2. However, i t was found that with the freedom to select an o p t i -mal f i r i n g angle, the 6 control was extremely small and that the r e s u l t i n g increase' i n range was i n s i g n i f i c a n t . (1) 6 = 0 , x f = 2,014,272 f t . (2) 8 ^ 0 , x f = 2,014,425 f t . The i n s i g n i f i c a n t increase i n range indicates that, i f the f i r i n g angle i s free, i t i s impractical to employ 8 control to maximize range. However, i f the f i r i n g angle i s f i x e d , a s i g n i -f icant improvement i n range can result through 8 control. Furthermore, the use of thrust angle control could be extremely important for optimal t r a j e c t o r i e s requiring some maneuverability of the rocket during f l i g h t . 4.4 The Fixed End Point Problem An interesting v a r i a t i o n of the type of problem handled i n Section 4=3 i s that of minimizing the f u e l to deliver a rocket between two fixed points i n space. For this example i t i s assumed that the i n i t i a l mass i s given''and that the f i n a l \ i mass i s to be maximized. As a r e s u l t , there' exists an additional unknown parameter which i s the instant of burnout. However, the addition of 7^ does not result i n a three dimensional search for the extremum of the performance function as one of the unknown parameters must be used to insure that the desired 48. f i n a l conditions are met. To i l l u s t r a t e t h i s example, consider the following data: t 0 = 0 m(0) =50 slugs x(0) = 0 x ( t f ) = 2,000,000 f t . y ( o ) = 0 y ( t f ) = 0 (4.29) v(0) = £>0 The performance function to be maximized i s P =m„(o , © ) and * Is s the p a r a m e t e r i s used to insure that (4.29) i s s a t i s f i e d for each set of parameters C&s> ®s)« ^ n e computational pro-cedure i s otherwise the same as Section 4«3« The r e s u l t i n g state variables at the switching points along the optimal trajectory are given i n Table 4.4. The var i a t i o n of m„ i n the neighbourhood of the optimum i s shown i n Table 4.4 I n i t i a l Point End of Max. Thrust Burnout P i n a l Point t(sec.) 0 0 + 29.99 376.3 x( f t . ) 0 0 85889 2000000 y( f t . ) 0 0 91630 0 v(ft./sec.) E £ 0 2408.0 7703.0 5977.0 ©(deg.) 53.3 53.3 43.3 -46.8 m(slug) 50.0 32.27 10.05 10.05 Figure 4.9 for m„ as a function of © with X optimum, and for i s s m_p as a function of X with © optimum. The optimal trajectory i s s i m i l a r i n form to that of Section 4.3° 4.5 The Direct Search by the Modified Relaxation Method The direct search employed i n these examples deals FINAL MASS Mt (SLUG) CQ FINAL MASS Mf (SLUG) — r ^ r-O CD '6t 50. with the problem of finding the extremum of a function P = P(a k) (4.30) over a set of parameters 'o^ for which the functional r e l a t i o n -ship (4.30) i s not known e x p l i c i t l y . In t h i s part of the thesis, the search i s accomplished through the use of a modi-fi e d relaxation method. In the standard relaxation technique, one parameter at a time i s varied while the remaining parameters are held fixed. The optimal value of the parameter being varied i s determined by a one dimensional search procedure as that value which maximizes P. The parameter then maintains t h i s optimal value while the next parameter i n the sequence i s varied and so on. The technique i s i t e r a t i v e and the cycle i s repeated u n t i l a l l the parameters converge to an optimal value. The one dimensional search procedure used at each step i n the i t e r a -t i o n i s based on finding a parabolic" approximation to the curve P(a^.) i n the v i c i n i t y of the optimum (see Figure 4 .10). From an i n i t i a l guess a^, the parameter <xk i s varied i n steps of Acck i n the d i r e c t i o n of increasing P u n t i l a value of <xk i s found which yields a larger value of P than for points on either side of oc^ . Assume that the following ine q u a l i t i e s hold P ( a k n - 1 ) < P ( a k n ) > P ( a k n + 1 ) ( 4 . 3 D As shown i n Appendix E, the three points with coordinates (P ,, a , 1 1 - 1 ) , (P , a, 1 1), and (P , , a, n + 1) can be used to n-1' k ' ' n' k ' n+1' k determine a parabolic approximation to the curve of the form P = a + ba k + c a k 2 (4.32) 51. Figure 4. 10 One Dimensional Search Employing a Parabolic Approximation The optimal value of oc^. i s taken to be the value of ex^. that maximizes (4.32). It i s shown i n Appendix E that t h i s value i s Vopt) - "k + ^ ( P ^ - V l ) For most problems, t h i s standard relaxation approach i s a convenient means of accomplishing the direct search i n parameter space. However, the speed of convergence i s depen-dent on the nature of the surface ?-(a^.) a n d often the conver-gence slows down long before the true optimum i s reached. A good example of th i s d i f f i c u l t y i s i l l u s t r a t e d by the search over the (3 , v ) plane for the maximum range problem of 52. Section 4.2. It i s seen i n Figure 4.11 that from an i n i t i a l guess (3 , v ) = (0.3, -2400), the standard relaxation technique o s appears to converge to a solution i n the v i c i n i t y of (0.24 , 2400). However, the true extremum i s at (0.1908, 2357.4). The cause for th i s "apparent" convergence i s that the contours of P(6 , v ) i n the (S , v ) plane are very nearly e l l i p s e s O S O S with a common major axis t i l t e d at approximately 45° to the coordinate axis. As the extremum i s approached, the eccentri-c i t y of the e l l i p s e s becomes increasingly smaller which causes the solution obtained by the relaxation method to o s c i l l a t e very rapidly, thus giving the impression of convergence. However, i t can be observed that i f the major axis were p a r a l l e l to one of the coordinate axes, the search technique would be f a i r l y independent of the ecc e n t r i c i t y . Indeed, i f the con-tours were true e l l i p s e s with a major axis p a r a l l e l to one of the coordinate axes, then one step convergence would r e s u l t . To benefit from this property, a modified relaxation method was developed. The technique i s based on rotating the coordi-nate axes such that the search i s along a new coordinate which i s p a r a l l e l to the di r e c t i o n of the major axis at that point. As t h i s "major axis" i s the projection of a ridge on the •P(ock) surface onto the hyperplane, i t i s generally not a straight l i n e and hence the direction^of the new coordinate axis has to be recomputed several times during the search procedure. In essence, the procedure i s as follows: (l ) Carry out one complete cycle of the standard relaxation technique to establish a point <^ a )^"'" on the ridge of the P(a, ) surface. 54. (2) Starting from <C°tjc/,"'"> carry out one more complete cycle of the standard relaxation method to establish another point on t h i s ridge. (3) Using <^a^.^> and ^x-^ determine the di r e c t i o n of t h i s ridge. Carry out a one dimensional search along t h i s d i r e c t i o n u n t i l an extremum of P i s found a t <ak>opt' (4) Repeat steps ( l ) , (2), and (3) from <(«k^ jp^ u n t i l the solution converges. (5) Investigate the surface ^ (o^) i n "the neighbour-hood of the solution found by the above procedure to prove that a true extremum has been obtained. The above technique i s called the modified relaxa-t i o n method and the r e s u l t i n g improved convergence i s i l l u s -trated i n Figure 4.12. It can be observed that since the search i s carried out i n a f i n i t e dimensional parameter space, i t i s possible to prove whether or not the solution i s a true extre-mum (see Figure 4.11 and 4.12). This fact i s further i l l u s -trated by Figure 4.13 which i s the contour map r e s u l t i n g from step (5) of the above method, Note that such a contour map cannot be obtained when the search i s carried out i n function space since a function can only be exactly represented i n a space of i n f i n i t e dimensions. 2330 2340 2350 2360 2370 2380 2390 Vs (FEET/SEC) Figure 4.13 Contour Map About the Optimum Found by the Modified Relaxation 57. 5. PROBLEMS EOR WHICH = 0 5.1 Introduction In Chapter 4, the problem of maximum range was studied using combinations of thrust and thrust-angle controls. In thi s chapter the three dimensional control problem i s i n t r o -duced which consists of thrust, thrust-angle, and l i f t control. The analysis associated with t h i s case i s given i n Appendix A. Problems consisting of one, two, and three dimensional control are studied for which the terminal conditions are independent of range and, hence, the transversality conditions y i e l d \^ = c^ = 0. The maximum alti t u d e problem, for the two dimen-sional control of thrust and thrust-angle, i s studied i n d e t a i l to show that the assumption u = cx>is v a l i d for most cases, r max. and to i l l u s t r a t e that optimal solutions obtained by the pro-posed techniques do indeed s a t i s f y a l l the necessary conditions of the calculus of variations. Subsequently, a problem :of maximizing a performance function at burnout i s investigated using one, two, and three dimensional control. It i s shown that t h i s problem i s equivalent to the maximum alti t u d e pro-blem i f burn-out occurs outside the earth:' s atmosphere and that, i n a l l cases, the search i n multi-dimensional function space i s reduced to a search over one time-invariant parameter. 5.2 The Three Dimensional Control Problem The development of the optimal control laws for the three dimensional control problem -proceeds .in-a; s i m i l a r manner to that of the two dimensional control case i n Section 4.2. Consider the case where the f i n a l time i s unspecified. The 58. t ransversality condition (A.21) yields c =0 (5.1) Using (A.16) and (A.17) i n (A .22), i t i s seen that for a non-t r i v i a l solution for \, and X. 3 4 tan 0 = DL (5.2) Due to (5.2), the f i r s t two rows :of matrix A are not indepen-dent and hence the f i r s t row of A may be eliminated. Rearranging the remaining rows of (A .25), the matrix A can be triangularized to y i e l d the matrix 0 B = 1 0 0 1 0 0 0 0 0 0 0 D (_ £ _ v v mv s i n 9 1 0 ) (k . _ % c o s 9 ) / 'mv v " v s i n 9 - cot 0/v 1 0 0 0 mv s i n 0 ve • (5.3) where B has the same rank as A and where .2, f T •= mg s i n ^ Q + 0) - D s i n 9 (sin 0 - cos^0 + v cos 0/v + D cos 9 s i n 0 cos 0 - v D^ s i n 9 cos 0 - L sin0 • •(2 s i n 9 cos 0 + s i n 0 cos 9 - v s i n 9/v ) (5.4) For the case f^ / 0 and ' c-^ ^ 0, (A.24) provides b ^ 0. It then follows from (5-3) that (2.47) i s s a t i s f i e d (rank (A) = rank (1) = 5 = m) and hence,the Lagrange mu l t i p l i e r s can be eliminated from (2.37), (2.38), and (2.44) y i e l d i n g equations for u and L i n terms of state variables only. Equation (5.2) i s used to define 0. Substituting (A.23) and (A.24) into 59.. (2.45) yields a set of l i n e a r , non-homogeneous equations i n \. This set of equations can "be used to express a l l the Lagrange mu l t i p l i e r s i n terms of c-^ for the case c^ ^ 0. However, to obtain a solution which i s v a l i d for c^ . = .0, the Lagrange mu l t i p l i e r s are expressed i n terms of \^ to obtain a set of equations of the form (2.49). Consequently, the optimal con-t r o l laws derived from t h i s set of equations w i l l be v a l i d for a l l c^. Following the procedure of Section (4.2), the desired optimal control laws obtained are L = X: + Y/sin 3 - v e u mD. LL s i n 3 Lv cos 3 D + vL >LL cos 3 0 = tan 1(D L) u = N/M (5.5) (5.6) (5.7) where X 1 D. f ( 2 ^ o s _ ^ : + . : s i n Q t a n p ) + D_Jan_^ cos 3 LL m 1_ LL h ( v mv W e cos 3 c Q s 2 p ) " D L y v s i n 9 + D L v(g s i n 9 - g) Y = IT =• cos 3 DLL • 2 Q sxn 3 Lv D D v s i n 9 + mv D L L cos^B YP . 2 mv cos 8 e r + v s i n 9 v m v D v y cos 8 - D (cos 8 - s i n g - v/v ) -r • Q y v p cos 8 ' eJ_] s i n 8 v Y D vL 1 cos + P(£ s i n 9» S + P ( 2 8 , O O B 9 + I ( M B _ f t . 2 ) ) e c o s s + D_c£sJ3 _ 5L4!L ! ) . ; + § S ( S I N Q +- 2 C O S 9 . mv m v 60. D . t a n p ) + Eg(D v _ 2 L tan_8.} _ ( + D}< H v mv m mv & m — + — _ L s i n P. + D ( S i ^ - f i + 3L.) + v D c o s ? v v v cos 3 v cos 8 v ' vv K e e r , e (L_ _ fi cos 9)Lg(2. s i n Q s i n 6 + cos 9(fjg_^ ^TTITT Y • o > . r COS p mv- cos , ) ) + D e i n p - - ^ ^ cos 8 •cos 8) + Y/cos 6 - D v s i n 8 cos 9 M =' — m 2 s i n 2 8 3 v e s i n 0 _ v cos 8 v v v D. v m 2 cos p + ^ (3 s i n ^ p)| + £ m v s i n 3 v cos 8 u e v , - 2 s i n 3 + (2 s i n 3 cos 3 - s i n ^ 3/cos 3) LI v cos p • ( ^ V : * - D v l oos 3) ' COS 3 . 2 E = mg s i n (9 + 8) + D(cos 8 - ^os 3^ " v / v e } - v V cos 3 + ' L ( ^ - f | t - | " 2 s i n 0) e O T TT P•= mg s i n 9 + D/cos 3 - D cos p/v - v D - — • • 3 6 , s i r r 3 2 cos 3 5«3 The Class of Problems for which c-^ = 0 For a l l rocket problems considered i n this thesis, i t 61. i s seen that the f i n a l transversality condition yields \^ = c^ = 0 whenever the f i n a l range i s unspecified. Two examples of t h i s type of problem are considered i n th i s chapter. The f i r s t i s that of maximizing the f i n a l a l t i t u d e for a given amount of f u e l . The second, which i s an approximation to the f i r s t and which eliminates the f i n a l coasting subarc, i s that of maxi-mizing the function at burnout. This function i s derived from the energy equation by di v i d i n g by the constant m^ g. (The subscript b denotes evaluation at the burnout instant The k i n e t i c energy term i n (5.9) i s associated with the y component of ve l o c i t y . Assuming that burnout occurs outside the earth's atmosphere, a l l aerodynamical forces are zero and the conservation of energy must apply. Consequently, (5.8) i s constant during the coasting subarc and maximizing (5.8) at burnout i s equivalent to maxi-mizing the f i n a l a l t i t u d e at (v s i n 9)^ _ t =0. Therefore, whenever burnout occurs at a r e l a t i v e l y high altitude (say above 75,000 f t . for t h i s case), the problem of maximizing-:^: at burnout i s a good approximation to the maximum alt i t u d e problem.. Both of these problems have c-^ = 0 and the resulting" analysis for the various cases i s given i n the following sections. 5.5.1 The Case of Thrust Control Only (5.8) ! y b = m bgy b + i m b ( v b s i n 9 f e) 2 (5.9) For the thrust control case of Section (4.1) (5.10) 62. However, as \^ = c^ = 0, (5.10) yields X = 0 and hence by equation (4.8) the optimal control law reduces to u = U/M (5.11) where M and N are given by (C-34) and (C-35) respectively. It i s seen, therefore, that for c^ = 0 the optimal control law (5.11) i s a function of state variables only. I t can also be shown that for the special case 9(0) = 90°, (5.11) reduces to the optimal control law for the sounding rocket |l2J . 5.3.2 The Case of Thrust and Thrust Angle Control The optimal control laws developed i n Section (4.2) are v a l i d for c = 0 and a l l values of c-^ . However, note that (Drl8) demands that b = 0 when c = c^ = 0. For th i s case, condition (2.48) requires that rank (A) = rank (B) = 4 and hence from (4.14) i t i s seen that f s = 0 (5.12) i s required everywhere along the variable thrust subarc where fg i s defined i n (4.15). This function fg i s analogous to the switching function f for the sounding rocket case, (see (3.4)). Also note, that for the special case 9(0) = 90°, f Q reduces to f s . 5 » 3 ° 3 The Case of Thrust, Thrust-Angle, and L i f t Controls The optimal control laws for t h i s problem are given by (5.5), (5.6) and (5.7). However, for the case c = c^ = 0, (A.24) gives b = 0 and condition (2.48) demands that rank (A) = rank (B) = 4 . ""From the d e f i n i t i o n of matrix B i n equation (5.3), i t i s seen that condition (2.48) i s s a t i s f i e d i f every-where along the variable thrust subarc 63. f L = 0 (5.13) where f^ i s defined i n (5«4). The function f ^ i s the switching function for the three dimensional control problem. Note that for the case L = 0, f^, (5.4)> reduces to f 0 , (4.15), and for the case L = 0 and ©(0) = 90°, f T reduces to f , (3.4). The function f^ i s , therefore, the general switching function from which the other switching functions can be derived. 5.4 The Maximum Altitude Problem for the Case of Thrust and Thrust-Angle Control Consider the two dimensional control problem, with terminal constraints t = 0 o x(0) = 0 y ( o ) = o v(0) = 1000 ft. / s e c . ©(o) = 70° e ( t f ) = o m(0) = 41.69 slugs m(t f) = 10 slugs (5.14) The performance function to be minimized i s P = - y ( t f ) (5.15) Substituting (5.14) and (5.15) into the transversality equation (DT15) yields X± = C l = 0 \ 2 ( t f ) = 1 \ 3 ( t f ) = 0 0 (5.16) 64. As c = = 0, the optimal control laws (4.16) and (4.20) apply, and condition (5.12) must be s a t i s f i e d everywhere along the singular subarc. Using the assumption that u = oo, equations (D-34), (1-37), and (DT40) completely define the system during the maximum thrust subarc. Hence the problem i s solved up to the unknown parameter 3(0) = 8 q. I t i s desired, however, to test the accuracy of the 3 control generated by (D-37) against the exact solution for the case of ^ m a x f i n i t e . To obtain the exact solution i t i s necessary to generate the Lagrange m u l t i p l i e r s during the maximum thrust subarc and to programme 8 according to (D-31) which i s 8 = t a n " 1 (^|-) (5.17) 3 Let IT be the instant of switching to the variable thrust sub-s arc and l e t x = x ( T ) where x i s any variable which depends on s s time. Then using the matrix B i n (4.14) and the fact that b = f f l '= 0, i t i s seen that at "X P s X2s (g + D/m s i n 9 + g cot 9 tan 3) * =^s (5.18) and X 3 e 5s - |m cos 8_|, _ (5.19) t = T S For t h i s problem i t i s convenient to scale the Lagrange multi-p l i e r s such that \^(0) = = 1 # After the optimal solution i s obtained, the Lagrange m u l t i p l i e r s can be rescaled to y i e l d the c l a s s i c a l solution \Q(t„) = 1. The computing algorithms for the 65. approximate and exact solutions are given i n the following sec-tions. 5.4.1 Approximate Solution for u. F i n i t e (1) For a given value u = u a n d the i n i t i a l conditions (5.14), select a value for 3(0) = 8 Q and using (D-57) to define 3, integrate (D-l) to (D-5) from t = 0 u n t i l f Q vanishes which defines T • P S (2) Use (4.16) and (4.20) to define 6 and u res-pectively, and continue integrating u n t i l mCt^) = m f (3) l e t u = 6 =0, and continue integrating u n t i l 9 ( t f ) = 0 which defines t f . Record y ( t ^ ) . (4) Return to ( l ) and perform a one dimensional search over B Q for the maximum y ( t ^ ) . 5.4.2 The Exact Solution for U. „ F i n i t e nictx • (1) For each value of u , use the value of 3 HL0.X o O found by the approximate technique as the i n i t i a l e s t i -mate. With = 1, select values for \ 2 Q a n d ^50° Using these values, solve (5.17) for \^Q. Integrate (D-l) to (D-10) from (5.14) u n t i l f^ vanishes which defines T . s (2) In general, \ 0 ( T ) and X^-iV) w i l l not s a t i s f y (5.1§) and (5.19) respectively. Return to ( l ) and select as an improved set of values X2G)(new) = \ 2 Q ( o l d ) + ^2s ~ ^ 2 ^ ^ \ 5 Q(new) = \ 5 Q ( o l d ) + \ 5 s - ^ ( l ^ ) 66. where \ 0 and \_ are found from (5.18) and (5.19), ahd: .\ 0 (Tl) and \ C-(T') are the actual values obtained. Repeat ( l ) and (2) u n t i l (5.18) and (5.19) are s a t i s -f i e d . (3) Use (4.6) and (4.20) to generate B and u respec-t i v e l y , and continue integrating u n t i l = m f (4) Set u = 6 = 0 and continue u n t i l 9(t^) =0 which defines t ^ . (5) Return to ( l ) and perform a one dimensional search over 6 Q for the maximum y(t^) s t a r t i n g with the values of and found by the previous i t e r a t i o n . (6) Rescale the Lagrange m u l t i p l i e r s such that \ 2 ( t f ) = 1. 5.4.3 Comparison of Exact and Approximate Solutions Table 5-1 i l l u s t r a t e s the exact and approximate values of the variables at the end of the maximum thrust subarc for u ranging from 2.0 slugs/sec. to i n f i n i t y . I t i s seen that the approximate solution for 8 given by (D-37) i s accurate to within one percent and that the difference i n f i n a l range i s i n s i g n i f i c a n t for a l l values of ^ - m a x greater than 10 slugs/sec. It can be concluded, therefore, that the use of (D-37) to pro-gramme B during the maximum thrust subarc i s j u s t i f i e d for a l l medium and high thrust rocket engines. 5.5 Testing the Necessary Conditions of the Calculus of Variations It was previously shown that the calculus of v a r i a -tions could be used to obtain: (a) a n a l y t i c a l forms for the Table 5-1 u max. (Slug/Sec) Exact Approx. $ Error (Sec) x s (feet) y s (feet) V s (ft/sec) , 9 s , (.rad.) m s (slug) Ps (rad.) y f (feet) OO E A foe 0 0 0 .2621,32 1.4453 30.82 0.1339 1026330 10 6 - E • A $e .11x10"! .11x10 4 0 3.6x10"^ 3.6x10 ^ 0 0.018E 0.0186 0 2621.28 2621.28 0 1.4453 1.4453 0 30V82 30.82 0 0.1339 0.1339 0 1026330 1026330 0 10 5 E A % e .•11x10"! .11x10 y 0 3.6xlO" 2 3.6x10 d 0 0.188 0.188 0 2621.28 2621.28 0 1.4453 1.4453 0. 30.82 30.82 0 0.1339 0.1339 0 1026330 1026330 0 10 4 •" E A foe • l l x l O " 2 .11x10 0 .36 .36 0 - 1.88 1.88 0 2621.38 2621.38 0 1.4453 1.4453 0 30.82 30.82 0 0.13389 0.13389 0 1026330 1026330 0 10? ' E A foe . l l x l O " 1 .11x10 ± 0 -3.6 3.6 0 18.81 18.81 0 2622.0 2622.0 . 0 1.4454 1,4454 0 30.815 30.815 0 0.13385 0.13385 0 1026330 1026330 0 1 0 2 " " - E • - -A foe - - .11 - • -.11 0 36.68 36.67 0.027 189.5 189.3 0.100 2628.44 2628.50 0.006 1.4458 1.4460 0.004 30.756 30.775 0.003 0.13338 0.13355 0.13 1026328 1026328 0 50 - • - E A fe • .22 .22 0 -73.74 73.71 0.04 - -381.7 381.7 0 ' 2635.78-2636.00 0.01 1.4462 1.4460 0.012 30.69 30.69 0 0.13287 0.13320 0.28 1026322 1026320 0.0002 10 • - E • - -A fe 1.157 ' 1.157 0 - 386.8 386.1 0.-12 2040. 20.40 0 2699-09 2699.20 - 0.004 1.4498 1.451 0.06 30.12 30.12 0 0.12857 0.12999 1.09 1026320 1026319 0.0001 2 - - E - -A foe 8.05 8.05 0 2647.0 2599.0 1.8 16230. 16250. 0.125 3301.74 3301.0 0.022 1.475 1.481 0.38 25.585 25.580 0.02 0.09886 0.1602 7.3 1026317 1026315 0.0002 68. optimal control laws, (b) the correct sequences of subarcs, (c) switching functions of state variables only, and (d) i n f o r -mation about the i n i t i a l values of the state variables. In most cases, once this information was obtained, the adjoinj; sys-tem could be completely disregarded and the optimal solution could be found by a search i n the parameter space of i n i t i a l conditions. I t i s now desired to show that the solution ob-tained i n t h i s manner i s a true extremal of the calculus of variations as manifested by the fact that along t h i s trajectory a l l the necessary conditions are s a t i s f i e d . Using the two dimensional control problem of the previous sections as an example, these necessary conditions are (1) from (2.9) ^ r\a.;(t)-X5) v X. cos 6 = jS. (x 3 s i n 8 - -) = 0, for u £ 0 (2) from the Legendre-Clebsch conditions (2.26) and (D - l l ) . v XA s i n 8 > K = (\, cos 0 + -2 : ) <0 u m • 3 v f o r u = umax.' 0 < u < u m a x . ' u = 0 respectively. (3) from the transversality condition (5.16) X± = c 1 = 0 A 2 ( t f ) = 1 \ 5 ( t f ) = 0 c = 0 (4) from (5.16) and the f i r s t i n t e g r a l (D-14) \ 4 ( t f ) =0 69. (5) from equations (5.12) and (5.4) V = .,0 ' f o r 0 < u < u m a x . i (6) and, from the Erdmann-Weierstrass corner condi-tions, a l l Lagrange m u l t i p l i e r s , state variables and the constant of integration are continuous functions of time. To test these conditions, an approximate solution was f i r s t obtained for u = 50 slugs/sec, (see Section 5.4.1). Using t h i s trajectory, equations (D-6) to (D-10) were integrated and the correct i n i t i a l conditions for the Lagrange mu l t i p l i e r s were determined as i n Section 5.4.2. The r e s u l t i n g search over P i s shown i n Figure 5-1 and the associated optimal controls for u and p are i l l u s t r a t e d i n Figure 5.2. I t i s seen from Figure 5.3 and Figure 5.4 that a l l the: ..necessary; conditions\(.i) : to (6) are s a t i s f i e d along t h i s trajectory and, hence, i t can be concluded 1 that the solutions obtained by the proposed tech-nique are true extremals of the calculus of variations. 5.6 The Problem of Maximizing the Function G^ at Burnout In this section, the problem of maximizing 2 r (,r j . (v s i n Q) x Gb = ^ + — \ „ i s investigated using one, two, and three dimensional control. For this problem the f i n a l time i s the instant of burnout t ^ . The terminal conditions are t = 0 9(0) = 70° x(0) - 0 m(0) = 41.69 slugs (5.20) 1026350 W26250\ 1026150\ 2; 1026050 20-0 20-66 /3Q (DEGREES) Figure 5 .1 The Final Altitude as a Function of 0 28-65 \1-1168\ co Uj Uj ct CD UJ Q V-566\ U=U MAX MAX THRUST SUBARC * NOT. TO SCALE VARIABLE THRUST SUBARC ~ 5 ~VMAX 0 + 6 9 TIME (SECONDS) 12 Figure 5.2 Optimal Controls u and 8, (Max. alt.) 8-0 s-oh 240 260 T5-~ 10-9 UMAX TIME-(SECONDS)- \jlME AFTER MAXIMUM THRUST~ Figure 5.3 The Optimal Lagrange M u l t i p l i e r s for the Maximum Altitude Problem ° r ~ 10>9 >s 5 UMAX % b 8 19 69 149 229 261 TIME- (SECONDS) - [TIME AFTER MAXIMUM THRUST] Figure 5 . 4 The Control Law Equations for the Maximum Altitude Problem 75. y(0) = 0 m(tf) = mf v(0) = 1000 ft/sec. and the performance function to be minimized i s ,2 1 P = - (v s i n 9)' y + 2g t = t f (5.21) A value of u =50 slugs/sec. i s used during the maximum max. D ° thrust subarc. The computing algorithms for the three cases studied are given i n the following sections and the acceptable sequence u = u , 0 < ^ u ^ u , u = 0 i s assumed, max. ^ max • 5.6.1 Computing Algorithm for Thrust Control Only (1) Select a value for v_ which represents the velo-s c i t y at the switching instant f ' between maximum thrust and variable thrust. Starting from (5.21) integrate (C-l) to (C-5) with u = u u n t i l the vel o c i t y increases to v which defines T. J s s (2) Programme u according to (5.11) and continue integrating u n t i l m = m^ . Record G^ . (5) Return to ( l ) and perform a one dimensional search over v for the maximum G, . s b Figure 5-5 i l l u s t r a t e s the search over v for the s maximum G^ and Figure 5.6 i l l u s t r a t e s the associated optimal thrust. 5.6.2 Computing Algorithm for Thrust and Thrust-Angle Control ( l ) Using u = u and defining 8 by (D-37), select a value for 8 Q and integrate (D-l) to (D-5) from (5.21) u n t i l fp vanishes which defines TT. 899000 Lu 889000 2 ct e Q , 879000 2500 3000 3500 4000 SWITCHING POINT VELOCITY (VS) - ( FEET/SEC) Figure 5.5 The Performance as a Function of v g l-40r o Uj to to CD ^) to 0-75 U=U MAX MAX —THRUST-SUB ARC * NOT TO SCALE •VARIABLE THRUST SUBARC-t o ll 03 |— To* 15-7 r 0 4 6 8 TIME-( SECONDS) 10 12 UMAX Figure 5.6 The Optimal Control u, (mgx. energy, 14 16 77. (2) Programme 8 and u according to (4.16) and (4.20) and continue integrating u n t i l m = m^ . Record G^ . (3) Return to ( l ) and perform a one dimensional search over 8 for the maximum G, . ro b Figure 5.7 i l l u s t r a t e s the search over 8 Q for the maximum G^ and Figure 5.8 i l l u s t r a t e s the re s u l t i n g optimal con-t r o l s for u and 8. 5.6.3 Computing Algorithm for Thrust. Thrust-Angle. and L i f t Control It i s assumed that for the case where l i f t i s non-zero, the drag force D(y,v,L) i s of the form |l7J D = K v V ^ + K j - e ^ l V 2 (5.22) where K = 10~ 4 a = 500 and a = (22,000 f t ) " 1 During the maximum thrust subarc u = ©o i s assumed. The equation for l i f t i s found by substituting (5.22) into (5.2) L = W[ ( v 2 e _ a y t a n P) (5.23) Using (A-l) to (A-5), (A-10) to (A-14), and (5-23), i t i s seen that 8 = s i n " 1 ( v Q s i n 8Q/v) (5.24) along the maximum thrust subarc. Under these assumptions, the res u l t i n g computing algorithm i s (l) For u = u , select a value for 8 . Define L ' max.' o^ and 8 by (5.23) and (5.24) respectively, and integrate 1033500T 78. 1031500\ 1029500 20 f30-DEGREES Figure 5.7 The Performance as a Function of 8( 23 2066- 1-45x 055 * NOT TO SCALE _L 0 % \ T ~ 1±6 s r's~UMAX~* 4 6 8 TIME SECONDS 10 12 14 16 Figure 5.8 The Optimal Controls u and 8, (max. energy) 18 79. (A-l) to (A-5) from (5.14) u n t i l f^ vanishes which defines f . s (2) Programme L, 0, and u according to (5-5), (5.6), and (5.7) respectively, and continue integrating u n t i l m = m^ . Record G^ . (3) Return to ( l ) and perform a one dimensional search over 8 Q for the maximum G^ . Figure-5*9 i l l u s t r a t e s the one dimensional search over S Q for the maximum of G^ , and the optimal controls are shown i n Figure 5.10. 5.6.4 A Comparison of the Three Cases Comparing these three cases with the case of a b a l -l i s t i c trajectory for which a l l the f u e l i s consumed during boosting, i t i s found that (1) B a l l i s t i c , Gfe = 681,000 f t . (2) u control , Gfe = 895,000 f t . (3) (u,S) control , Gfe = 1,033,400 f t . (4) (u,S,L) control , Gfe = 1,036,400 f t . As predicted by theory, the value of G^ increases with an i n -crease i n the control vector. In p a r t i c u l a r , s i g n i f i c a n t increases are realized between cases ( l ) and (2), and between cases (2) and (3). The increase between cases (3) and (4) i s r e l a t i v e l y small; however, i t can be observed from Figure 5.8 and Figure 5.10 that the demand on 8 control i s reduced when l i f t control i s added. This may be an important feature since, i n a l l p r a c t i c a l cases, there w i l l be an upper and lower bound 705 8 0 . 900-\ fb0 (DECREES) Figure 5-9 The Performance as a "Function of 3 Q *N0T TO SCALE 4 6 8 K7^ ' L C A X ^ T I M E (SECONDS) Figure 5.10 The Optimal Controls u, p, and L. (iwax. energy) on 3 which could be violated by case (3) and not by case (4). Furthermore, i t may be economically advantageous to reduce 0 control at the expense of l i f t control. The answers to such problems of course w i l l depend on the p a r t i c u l a r system under study and the type of trajectory desired. Certainly those tra j e c t o r i e s which require much maneuverability would favour the use of both (3 and l i f t controls. 82. 6. SUBOPTIMAL CONTROLS 6.1 Introduction It i s sometimes economically advantageous to trade off a loss i n system performance for a s i m p l i f i c a t i o n i n the design of the optimal controller. This type of consideration leads to the area of suboptimal control. By d e f i n i t i o n , a suboptimal con-t r o l i s any control, other than the optimal, which takes the system from a given i n i t i a l manifold to a given f i n a l manifold. This type of control always experiences a loss i n system per-formance . However, as shown i n t h i s chapter, under certain conditions t h i s loss may be i n s i g n i f i c a n t . For an example, the maximum range problem i s studied and two means of generating suboptimal controls which are functions of state variables only are presented. 6.2 Eliminating ^ from the Control Equations Consider the optimal control problem of Section 4«3. For this problem, a value Y = 0.8x10"^ i s found to be the s optimum value of at t = T^. During the remainder of the trajectory ^ i s generated from the d i f f e r e n t i a l equation (4.10) using X as the i n i t i a l value. However, th i s solution i s s v a l i d only i f no disturbances occur during f l i g h t . Should a disturbance occur at t = t ^ , the optimal trajectory for T T - ^ t ^ t ^ requires that a new optimal value of ^ ^ ^ ( t ^ ) be found. To i l l u s t r a t e t h i s property, disturbances to the v e l o c i t y and path i n c l i n a t i o n are provided at various times during the optimal trajectory found i n Section 4.3. Three cases are studied. F i r s t , the new optimal value of ^ i s found', second, no change i n the value of $ i s made.; and t h i r d , 0* i s 83. kept equal to zero throughout the trajectory. The results are shown i n Table 6.1 and i t can be observed that for perturbations i n state variables up to 10$, the r e s u l t i n g difference i n f i n a l range for the l a s t two cases i s less than 1% of the optimum. Furthermore, since X i s l o c a l i z e d to zero for a l l disturbances, the case for which $ i s kept equal to zero appears to be a good choice for the thrust control equation. Using X = 0 in^4-.8), a suboptimal control of the type u s o l = N/M . (6.1) i s developed where N and M are given by (C-34) and (C-35) respectively. Note that the suboptimal control (6.1) i s a function of state variables only. 6.3 Using 0 = 0 i n the Two Dimensional Control Case In Section 4.2, the optimal controls for the two dimensional control problem were found to be (see (4.16) and (4.20) ) i = ^ [j> mg cos 9 - §g_| <D(1 + 2- C O S P) . M G S I N Q> (6.2) and u = N/M (6.3) where N and M are functions of state variables only as given by (4.21) and (4.22) respectively. The object here i s to reduce this two dimensional control problem to a one dimensional con-t r o l problem by forcing 0 to be i d e n t i c a l l y zero. This requires that (6.2) vanish. Defining / s i n 0 (6.4) v u s i n 0 e r R t D(l + — cos 0) - mg s i n 9 v e and evaluating (6.2) for 0 = 0, i t i s seen that 0 i s zero only i f SCALE: V(ft/sec) , Q(rad), y ( f t ) , and Eime *1 Values Before Values After Percent Change S v , s© xlO 5 y ( t f ) With X0(tJ) y ( t f ) With y(t f) With io Error; With S o ^ % Error With X = o V\f-j_) © ( T - L ) v(rj) o(rj) 0 2410. .8960 — — — .8 • — 2014212 - — 2014027 . — .009 5 2752.4 — . 3000 — +9.1, - .7345 .IxlO" 4 2139186 2140092 2139638 .0425 .0212 5 2752.4 — 2500 — -9.1, - .734 .5xlO" 4 1886340 1888260 1885098 .1020 .168 5 — .8960 — .92 -, 3.4 .734 -.IxlO" 4 2009796 2010262 2010239 .0425 .0204 5 — .8960 — .86 -, -3-4 .734 .5xlO" 4 2006504 2009293 2005132 .139 .207 10 3 2 0 0.6 — 3500 — 10, 0. .662 . I x l O " 1 0 2164927 2165166 216-5166 .011 .00 10 3 2 0 0.6 — 3000 — -6.66, - .662 .5xlO - 4 1915000 1915257 1914524 .010- .0250 10 — .8613 — .93 -, 8.1 .662 -.IxlO" 4 1986387 1987732 1986931 .068 .0405 10 — .8613 — .85 -, -1.17 .662 .IxlO" 4 2013959 2014105 2013721 .0073 .0241 10 3200.6 .8113 3000 .9 -6.7,5.0 .662 .1x10 1908300 1908576 1908576 .010 .000 20 4646.8 — 5100 — 10.08, - .487 -.IxlO" 4 2240427 2240541 2240541 .005 .000 20 4646.8 — 4200 — -9.8, - .487 .5xlO" 4 1802958 1803487 1802773 .029 .038 20 — .8043 — .88 -, 10.0 .487 -.5xlO" 4 1987371 1988354 1987438 .0495 .0461 20 — , .8043 — .72 -, -10.0 .487 .IxlO" 5 1988129 1990283 1987753 .108 .127 Table 6.1 85. 2mg cos 9 = R (6.5) Substituting (6.5) and 0 = 0 i n (6.3) yields uso2 = V M o ( 6 - 6 ) where IVT 2 2 2 p. D p. N Q = -m g cos 9 + + mD e •sin 9 + —) m ( l + ^-)(3g s i n 9 + av-e (6.7) and M = D o 4v + 2v e + v2/v^ j (6.8) Equation (6.6) i s the desired suboptimal control. However, for 8 = 0 , (6.4) yields a f i n i t e R i f and only i f D(l + 2_) - mg s i n 9 = f n = 0 (6.9) V ^ B U e Di f f e r e n t i a t i n g (6.9) with respect to time and using (D-l) to (D-5), i t can be shown that (6,9) i s s a t i s f i e d i f u = V M o = uso2 ( 6 ' 1 0 ) which i s consistent with (6.6). 6.4 Oomparison of Suboptimal and Optimal Controls The maximum range problem of Section 4.3 was solved using the suboptimal controls (6.1) and (6.6). Comparing with the optimal control, the results are: (1) u optimal , x f = 2,014,425 f t . (2) U S Q 2 , x f = 2,014,212 f t . (3) u g o l , x f = 2,014,156 f t . It i s seen that the suboptimal controls are excellent approxi-mations to the optimal. Also, the advantage of the suboptimal controls, u -, and u 0, i s that they are functions of state ' s o l so2 17 86. variables only and hence can be generated by standard feedback techniques. u s 02 has the added advantage that the function f ' i n (6.9) must be i d e n t i c a l l y zero during variable thrust and hende the thrust can be generated by (6.9) and a high gain amplifier j^9j. I t can be concluded, therefore, that since the suboptimal controls are simpler to implement, and since the loss i n system performance i s n e g l i g i b l e , the use of (6.1) or (6.6) to generate the thrust i s j u s t i f i e d i n t h i s case. In a si m i l a r manner, suboptimal controls for other cases could be generated and tested. As the true optimum for each case is' known, a performance measure can be assigned to any candidate for sub-optimal control. 86 7. OPTIMAL CONTROLLERS DURING SINGULAR SUBARC 7.1 Introduction Por the purpose of synthesizing optimal cont r o l l e r s , the various optimal control laws developed i n the previous chapters can be separated into three cases: ( l ) the case for which the control i s a function of state variables only, (2) the case for which one control i s integrated from a d i f -f e r e n t i a l equation, and (3) the case for which a function of state and control variables exists which i s to be zero through-out the variable thrust subarc. The form of the optimal con-t r o l l e r for each of the three cases i s different as i l l u s t r a t e d i n the following sections. 7.2 Direct Feedback Control For the optimal control, law i n (5.1l), and for the suboptimal controls (6.1) and (6.6), the controller i s obtained by standard feedback techniques as i l l u s t r a t e d i n Figure 7.1. 7.3 Hybrid Optimal Controllers For the second type of cont r o l l e r , one parameter exists whose value must be up-dated along the trajectory to account for disturbances to the system during f l i g h t . The class of problems which require this type of controller are those for which c^ i s not equal to zero and one of the controls i s generated from a d i f f e r e n t i a l equation, (see Chapter 4). The r e s u l t i n g controller i s of the hybrid computer variety which uses an analog simulator to perform high-speed trajectory com-putations and a h i l l - c l i m b i n g d i g i t a l computer to carry out the one dimensions search 9 • Figure 7.2(a) shows the controller U=U(X) Figure 7.1 The Direct Feedback Optimal Controller HILL CLIMBING DIGITAL COMPUTER X(0) Y(0) X=f(X.U) u=U(X. X) TRAJECTORY COMPUTER Figure 7.2a The Hybrid Optimal Controller - Maximum Ran ere with Thrust Control HILL CLIMBING DIGITAL COMPUTER I X(0) X=f(X.U) U=U(X.fr) /3=/3(X./3) TRAJECTORY COMPUTER Figure 7.2b The Hybrid Optimal Controller - Maximum Range With Thrust and Thrust Angle Control 88. for the case of thrust only, and Figure 7.2(b),: shows the con-t r o l l e r for the two dimensional control of u. arid 6. 7.4- Implicit Function Generation The t h i r d type of controller deals with those cases for which a function of state variables and at most one control variable exists which must be zero throughout the singular subarc. In general, the control cannot be solved e x p l i c i t l y from this function and an i m p l i c i t solution i s required. F i g -ure 7.3 shows the case of the sounding rocket and the subopti-mal control of Section 6 . 3 . Figure 7 .4(a) shows the case of the two dimensional control of Section 5 . 3 . 2 , and Figure 7 .4(b) shows the case of the three dimensional control of Section 5 . 3 « 3 . 7.5 Conclusions It has been shown that for systems whose dynamics are l i n e a r i n control u, i t i s possible to derive control equations for u, (3, and L which are functions of state variables only for a variety of optimization problems... Furthermore, these control equations are convenient f or the study of optimal and suboptimal feedback control laws which can be implemented by direct feedback, up-dating the parameters through hybrid computation, or by i m p l i c i t solution of a switching function to obtain the desired control. A l l unknown parameters oc^. which enter into the problem are found by a dir e c t search in.a parameter space for the minimum of the .performance, function-.;.^.)••• It was shown that a modified relaxation method i s a suitable technique for accomplishing th i s search and that, since the search i s carried out i n a u PLANT Figure 7.3 Implicit Solution Using a High Gain Amplifier U PLANT ft U=U(X,ft) ft ft IMPLICIT SOLUTION f(X./3) =0 Figure 7.4a Implicit Solution Controller - Two Dimensional Control With ^ = 0 X U-U(X.L) L IMPLICIT SOLUTION f(X,L)~0 Figure 7.4b Implicit Solution Controller - Three Dimensional Control with c-, = 0 90. parameter space of f i n i t e dimensions, the solution can be eas i l y tested to insure that i t i s a true extremum and not merely a stationary point... The control components u, 0, and L can then be generated from the state variables and the optimal parameters. For the class of systems given by (2.1), the proposed technique i s considerably more convenient than standard numeri-c a l procedures which require not only a search i n multi-dimensional function space but are also' unsuitable for real-time control by i n - f l i g h t guidance computers. PART I I NUMERICAL ALGORITHMS 8. NUMERICAL TECHNIQUES 91. 8.1 Introduction In Part I of th i s thesis, the optimal control laws for a class of aerodynamical systems were obtained as a func-t i o n of state variables and, at most, one time-invariant para-meter. These control laws provided an e f f i c i e n t means of gener-ating the optimal t r a j e c t o r i e s and allowed the p o s s i b i l i t y of implementing real-time control. In general, however, such a n a l y t i c a l forms for the optimal controls cannot be obtained and i t becomes necessary to employ numerical techniques for the solution of the optimization problem. As mentioned previously, these numerical methods are b a s i c a l l y i t e r a t i v e schemes which require the use of large size d i g i t a l or hybrid computers. Although some success has been realized with these techniques, problems s t i l l e x i s t - i n the areas of i n i t i a l and f i n a l conver-gence, computer storage requirements, and computational algor-ithms. In th i s part of the thesis.,. numerical algorithms are discussed which are es s e n t i a l l y a combination of the direct and indir e c t approaches and which a l l e v i a t e some of these present d i f f i c u l t i e s . I t i s shown that the concepts used to develop these new algorithms can also be used to improve the properties.of exi s t i n g techniques. E s s e n t i a l l y , there are three basic con-cepts used: (l ) a f i r s t v a r i a t i o n approach applied to the aug-mented performance function which results i n a gradient search i n the parameter space of i n i t i a l Lagrange m u l t i p l i e r s , 92. (2) a second v a r i a t i o n approach applied to the aug-mented performance function which determines the o p t i -mal step size for the gradient approach i n ( l ) , and (3) an approach which determines the optimal scale factor for the Lagrange m u l t i p l i e r s such that the error i n f i n a l transversality i s a minimum at each step i n the i t e r a t i o n . It i s shown that for algorithms based on ( l ) , the scale factor f or the Lagrange m u l t i p l i e r s i s arb i t r a r y , and instead of searching over the entire \ Q-space, i t i s s u f f i c i e n t to determine the intersection of a l i n e with any sphere T \ \ = constant. Consequently, the i n i t i a l convergence does not depend on a good estimate of the optimal trajectory. Furthermore, as the gradient search i n (l) i s performed i n parameter space, computer storage i s required at the terminal points only. A disadvantage to ( l ) i s that, since i t i s a gradient technique, the convergence slows down as the optimum i s approached and i t i s not known when the search should be terminated. To overcome t h i s d i f f i c u l t y , (2) and (3) are used to determine the optimal step size f or the gradient technique i n the v i c i n i t y of the extremum. Concepts (2) and (3) are also applied to the method of steepest descent and the ind i r e c t methods based on matching end points. I t i s shown that some of the undesirable properties previously associated with these tech-niques can be s i g n i f i c a n t l y reduced. In t h i s chapter, the fundamental concepts of numerical methods w i l l be discussed and some of the ex i s t i n g numerical 93. techniques w i l l he presented. 8.1.1 The Direct and Indirect Approaches The direct and in d i r e c t approaches are i t e r a t i v e schemes which start from some i n i t i a l estimate of the optimal trajectory and generate a series of tra j e c t o r i e s that eventually converges to the optimum. Each trajectory i n the series i s ob-tained by a search i n the neighbourhood of the previous t r a j e c -tory, called the nominal trajectory, for that trajectory which best s a t i s f i e s the search c r i t e r i o n . As a r e s u l t , the new trajectory i n the series i s , i n some sense, "better" than i t s predecessor and i s therefore closer to the optimum. This procedure i s repeated u n t i l the search c r i t e r i o n i s s a t i s f i e d . Thus, i t can be observed that there are three basic features of these numerical techniques. These features are based on the manner of defining ( l ) the nominal trajectory about which the search i s conducted, ( 2 ) the space i n uhich the search i s carried out, and ( 3 ) the c r i t e r i o n upon which the search i s based. The method of generating the neighbouring t r a j e c t o r i e s i s common to a l l techniques and i s based on a technique of l i n e a r i z a t i o n about a nominal trajectory. This technique w i l l be discussed i n the following sections beginning with a review of l i n e a r system theory. 8.2 Linear Time-Varying D i f f e r e n t i a l Systems jl6J To begin, consider the zero input response and the forced response of systems described by x(t) = A(t)x(t) + B(t)u(t) (8.1) where A(t) i s an n x n matrix of scalar functions assumed to 94. be continuous for a l l t ; s i m i l a r l y , B(t) i s assumed to be an n x m continuous matrix, x(t) i s the state vector, and u(t) the input. 8.2.1 The Zero Input Response Theorem 1: Let 'I ( t , t Q ) be an n x n matrix which i s the solution of the matrix equation d fi(t,t ) o A(t) 5 ( t , t ) (8.2) dt - u y ~ v o where (^"^ 0» t Q ) = I Then the zero-input response of (8.1) x(t) = A(t)x(t) , x ( t Q ) = Xq (8.J) i s given by x(t) = 5 ( t , t o ) x ( t Q ) V t , V x 0 (8-4) Proof: By the d e f i n i t i o n of $ ( t , t Q ) , observe that (8.4) reduces the x at t = t Q . F i n a l l y , (8.3) i s s a t i s f i e d by d i f f e r e n t i a t i n g (8.4). The matrix $( t , t ) i s called the state t r a n s i t i o n matrix for the ' o system (8.3)• 8.2.2 The Forced Response Theorem 2: Let $ ( t , t Q ) be defined by (8.2). Then the forced response of (8.1) which goes through x Q at t Q i s given by t x(t) = s ( t , t o ) x o + i ( t , a ) B ( a ) u ( a ) d a (8.5) V t , V x Q 95. Proof: The proof i s immediate by direct v e r i f i c a t i o n of i n i t i a l conditions and direct substitution of (8.5) into (8.1). To effect t h i s substitution note that t t | ^ 5(t,CQB(a)u(CL)dCL= B(t)u(t) + A(t) ^ / ^ S ( t , a > t t o o •B (a)u (a)da (8.6) In practice, the part i c u l a r solution which i s given by the in t e g r a l i n (8.5) i s obtained by solving 1 th (8.1) with x Q i d e n t i c a l l y zero. Also, the i column of l ( t , t ) i s obtained by finding the zero input solution to (8.1) with xj_(t 0) equal to unity and the remaining i n i t i a l conditions equal to zero. 8.2.3 The Ad .joint, System Por the input response of (8.1) x(t) = A(t)x(t) (8.7) the system defined by Z(t) = -A T ( t ) Z ( t ) (8.8) i s called the adjoint system. According to Theorem 1, a state t r a n s i t i o n matrix SP(t,t Q) exists for the adjoint system such that | ^ T ( t , t Q ) =- A T ( t M t , t o ) (8.9) where «(t o,t Q) = I As a r e s u l t , the solution of (8.8) which passes through Z Q at t i s given by •Z(t) =«(t,t )Z Vt,\/Z Q (8.10) Lemma 1 Lemma 2 Lemma 3 x T ( t ) Z ( t ) = constant (8.14) 96. The following lemmas, which relate the system and adjoint system, w i l l be stated without proof: (see |l6j) vP(t 1,t Q) = g T 1 ( t o , t 1 ) (8.11) $ T ( t 1 , t o ) v i ' ( t 1 , t o ) =1 (8.12) ^ ( t Q , t 1 ) = 5 T ( t 1 , t 0 ) (8.13) Lemma 4: 8.3 Linearization About a Nominal Trajectory In optimization problems, the system to be controlled i s generally described by a set of nonlinear d i f f e r e n t i a l equa-tions of the form x(t) = f(x,u) (8.15) where x i s the state vector of n components, u i s a control vector of m components, and f(x,u) i s an n x 1 vector whose components are continuous functions of x and u. L4i-"a',;nominal trajectory x(t) of (8.15) be defined by some control u(t) and a set of i n i t i a l conditions x(t ). It i s desired to examine the effect of perturbing the i n i t i a l state by SXq and perturbing the control bySu(t). Equation (8.15) can be expanded i n a Taylor series about the nominal trajectory to y i e l d Sx = A(t)Sx + B(t)Su (8.16) where A(t) = ^ , B(t) = *g and Sx(t ) =Sx o o and where the partial derivatives are evaluated along the 97. nominal trajectory. Applying Theorem 2 to the l i n e a r time-varying system (8.16) yields t Sx(t) = s(t,t )8± + C ffi(t,a)B(a)8 u(Q.)d(H (8.17) t o where $( t , t ) i s the state t r a n s i t i o n matrix of the zero input 'o response of (8.16). The tr a j e c t o r i e s given .by x(t) = x(t) +Sx(t) (8.18) are the tr a j e c t o r i e s i n the neighbourhood of the nominal t r a j e c -tory and they are functions of SXq andSu(t) as given by (8.17). 8.4 The Optimal Control Problem To i l l u s t r a t e the basic principles of the various numerical techniques and because there i s no loss i n generality, a somewhat si m p l i f i e d optimal control problem w i l l be used as an example. Extension of the techniques to problems involving free f i n a l time, additional end constraints, bounded control, etc., may be obtained by consulting the references given i n each section. The problem w i l l be to fi n d that set of controls u(t) which w i l l minimize the system performance function J = 0(x(T)) (8.19) subject to the constraints x = f(x,u) (8.20) x(0) = x (8.21) o over ~the fixed time i n t e r v a l O ^ t ^ T . (8.22) 8.4.1 The Necessary Conditions for a Local Extremum,: j^ 21 The constraints (8.20) are adjoined to the perfor-mance function ( 8 . 1 9 ) by means of an n x 1 vector of Lagrange mu l t i p l i e r s X to y i e l d the augmented performance function T J a = 0(x(T)) + J ~ (\ Tx -H)dt ( 8 . 2 3 ) 0 where H = \ T f ( 8 . 2 4 ) i s the v a r i a t i o n a l Hamiltonian. The problem of minimizing J i s thus transformed to a problem of minimizing J . Taking the f i r s t v a r i a t i o n of ( 8 . 2 3 ) yields T S j a =S0(x(T)) (S\Tx + \^Sx - § H ) d t ( 8 . 2 5 ) 0 §0(x(T)) = 8 x f T 0 x f ( 8 . 2 6 ) §H =8U T H + 8X T H + SA T H . ( 8 . 2 7 ) "LI X A, where and ' ~lT -T xdt =Sx T\J 0 - ^ S x T \ d t ( 8 . 2 8 ) 0 0 Substituting ( 8 . 2 1 ) , ( 8 . 2 6 ) , ( 8 . 2 7 ) and ( 8 . 2 8 ) into ( 8 . 2 5 ) yields T S j a = S x / [0 x f + Xf] - J ( S u \ + Sx T(H x, + X) 0 + S\ T(H X - x)} dt ( 8 . 2 9 ) where the subscript f denotes evaluation at the f i n a l time T. Hence, for J to be a minimum, i t i s necessary that S j be a a zero. Equating ( 8 . 2 9 ) to zero, for independent variations i n Su, Sx, and Sx, provides the following set of necessary con-d i t i o n s : 99. (l) the system equations, x = H. = f (x,u) (8.30) (2) the Euler-Lagrange equations, (8.31) (3) the gradient condition (control equation), 0 = H u = V * (4) the i n i t i a l conditions, x(0) = x„ (8.32) (8.33) (5) and, the f i n a l conditions ( t r a n s v e r s a l i t y ) , .\f = " 0 x f (8.34) where the short-hand notation - A 3 f *x = 5x = 'lx, "nx. and 0 X A aei(T) f - 3 x 0 X (T) X l 0 X (T) n 'lx.. n "nx n (8.35) (8.36) 100. i s used to denote p a r t i a l derivatives of vectors and scalars respectively. A solution which s a t i s f i e s the above necessary condi-tions i s called an extremal solution. Such a solution i s a candidate for optimality but i s not necessarily the global optimum since the conditions (8.30) to (8.34) are ( l ) l o c a l i n nature and"(2) generally not s u f f i c i e n t . An additional test, such as a second v a r i a t i o n t e s t , i s required to separate the l o c a l optima from the extremals, and, subsequently, a search over a l l l o c a l optima i s needed to determine the global optimum. However, for the present argument, i t i s t a c i t l y assumed that the l o c a l optimum, and the extremal are unique so that only the conditions (8.30) to (8.34) need be considered. Based on these assumptions, the general approach to obtain a numerical solution i s as follows: (1) select a nominal trajectory which s a t i s f i e s as many of the necessary conditions as possible, (2) determine the space over which the search i s to be conducted by selecting those parameters and/or functions which w i l l be perturbed to generate the neighbouring t r a j e c t o r i e s , and, (3) select as a search c r i t e r i o n a direct approach (most improvement i n systems performance) or an i n d i -rect approach (most improvement i n meeting the neces-sary conditions not s a t i s f i e d i n ( l ) ) . To i l l u s t r a t e t h i s approach, several of the more common numeri-ca l techniques w i l l be presented i n the following sections. 101. 8.5 The Method of Steepest Descent |6, ^ For the method of steepest descent, the search i s carried out i n the function space of the control vector u ( t ) . An i n i t i a l control u(t) i s selected to define the nominal trajectory. The system equations (8.30) are integrated forward from (8.33) at t = t u n t i l t = T. The Euler-Lagrange equations (8.31) are integrated backward from the f i n a l conditions (8.34) using the values of x(t) obtained i n the forward integration. As a r e s u l t , only condition (8.32) i s not s a t i s f i e d and hence equation (8.29) reduces to T S j a = - J ^ (Su THu)dt (8.37) 0 The method of steepest descent i s based on the direct approach and a neighbouring trajectory i s sought which results i n a minimum of (8.37). Using Schwarz's i n e q u a l i t y , 8 J i s a minimum when Su(t) = kH u . (8.38) A where k i s a positive constant and H is, evaluated along the nominal trajectory. To insure that the l i n e a r i t y requirements are not vio l a t e d , a constraint on Su(t) i s imposed such that ^JI)u T(t)Su(t)dt =Sl 2 (8.39) 8:2 0 where 01*" i s chosen a r b i t r a r i l y small. The new nominal t r a j e c -tory i s -., ,:u(t) = u +Su(t) (8.40) 102. where u(t) i s defined by (8.38) and i s subject to the constraint (8.39). This process i s repeated u n t i l Su(t) goes to zero. It i s seen from (8.38) that t h i s condition w i l l occur when i s i d e n t i c a l l y zero and hence the remaining necessary condition (8.32) w i l l be s a t i s f i e d . The main advantage of this method i s that i n i t i a l con-vergence does not depend on a good i n i t i a l estimate of the I optimal control u * ( t ) . The disadvantages are that computer storage i s required at many points along the trajectory and that the convergence slows down as the optimum i s approached. In Section 1 (10.4), a second v a r i a t i o n technique i s developed :. which determines the optimal value or the parameter k i n the v i c i n i t y of the extremum. This modification to the steepest descent technique provides a means of improving f i n a l conver-gence without s i g n i f i c a n t l y increasing.the computational require-I ments. 8.6 The Min-H Strategy QT|* The min-H strategy i s s i m i l a r to the method of steep-est descent except for the c r i t e r i o n upon which Su(t) i s selected. In t h i s process, a Su(t) i s found which drives H closer to zero. Consequently, as the technique i s based on s a t i s f y i n g the remaining necessary condition (8.32), the min-H strategy i s an ind i r e c t approach.' The name "min-H" i s derived from the fact that H i s zero when H i s a mimimum (condition A (8.32)). For the nominal trajectory, H^ w i l l i n general not be zero1. Expanding H^ i n a Taylor series y i e l d s H = H + H S u + H ^ S x + H Sx (8.41) u u uu . u\ ux *Correctly,i the choise of sign i n (8.34) makes this "Max-H". 103. where the p a r t i a l derivatives are evaluated along the nominal trajectory. Equating (8.41) to zero such that condition (8.32) i s s a t i s f i e d , the desired v a r i a t i o n i s 8u(t) becomes S u ( t ) ^ H ^ " 1 $J>x + SnXBx + E j (8.42) Expanding (8.30) and (8.31) i n a Taylor series and replacing Su(t) by (8.42) yields S i ( t ) S\(t) S x ( t ) S\(t) (8.43) where C = o f - f H ^ H X u uu ux A A T l A _ T A -H + H H H XX ux uu ux d(t) e(t) A TA _T A - f H X f u uu u A T I A T l A — T A HP - f 1 + H 1 H • -1- f X ux uu u and , . A A _"| A d t) '= - f H H ' u uu u , , A A _ - i A e t = •+ H H H ' xu uu u (8.44) (8.45) (8.46) and where a l l p a r t i a l derivatives are evaluated along the nominal trajectory. As §x Q = 0, aSx. o must be determined that w i l l preserve the desired f i n a l condition (8.34). Expanding (8.34) about the nominal trajectory yields the desired change i n S \ . p = - 3 - Sx^. (8.47) This desired change i n f i n a l value can be transferred (or "swept") back to the i n i t i a l point by means of the R i c c a t i 104. transformation. 8\(t) •= P(t) Sx(t) + W(t) (8.48) Substituting (8.48) into (8.43) yields 8 i = (A+BP)Sx + BW + d - (8.49). &{•= (C-ATP)Sx - ATW + e (8.50) where A = f x - f u - ^ H^" 1 H ^ (8.51) B = - f H ± f 8.52) u uu ; u A A TJ A _ 1 A . . and C = -H + H H E v (8.53) XX ux uu ux " D i f f e r e n t i a t i n g (8.48) with respect to time and using (8.49)' gives S\ = (P + PA + PBA) Sx + PBW + Pd + W (8.54) Equating (8.54) to (8.50) for Sx arbitrary yields P = -PA - PBP - A TP + C , P(T) = - 0VY.(T) (8.55) ¥ = -PBW-Pd-ATW. + e , W(T) = 0 (8.56) Substituting (8.48) into (8.42) gives S . . A n p A A A C* A A A - , . . where Sx i s determined from Sx = x(t) rr x(t) ! (8.58) The procedure i s repeated u n t i l i s driven to zero and hence a l l the necessary conditions are s a t i s f i e d . The i n i t i a l conver-gence for t h i s technique i s not as good as the method of steepest descent.. However, t h i s approach offers good f i n a l convergence; i n f a c t , the speed of convergence increases as the extremum i s approached. The computational algorithm, however, i s much more 105. complex and computer storage i s required at many points along the trajectory. 8.7 The Newton-Raphson Technique j^ 2 J For the Fewton-Raphson technique, the search i s carried out i n the function space of the state variables x(t) and the adjoint variables \ ( t ) . An i n i t i a l guess at the time functions x(t) and \( t ) i s made such that the boundary conditions (8.33) and (8.34) are s a t i s f i e d . The control u(t) i s determined by solving the control equation (8.34). I t i s assumed here that (8.32) can be solved e x p l i c i t l y f o r u i n the form u = u(x,.\) (8.60) This assumption however, i s not a r e s t r i c t i o n on the numerical techniques and the case where (8.60) cannot be found e x p l i c i t l y i s covered i n Section (9.7). As a resu l t of (8.60), the ,only necessary conditions which are not s a t i s f i e d are the state equations (8.30) and the Euler-lagrange equations (8.31). The search c r i t e r i o n i s to fin d the neighbouring trajectory for which (8.30) and (8.31) are more closely s a t i s f i e d . Substituting (8.60) into-(8.30) and (8.31) yields two functions of the form h(t) = x - f ( x , u(x,;X) = x - r( x , \) (8.61) p(t) = X + f T ( x , u(x,\))\ = X * s(x,\) (8.62) and hence the necessary conditions (8.30) and (8.31) can be replaced by the conditions h(t) = 0 (8.63) and p(t) = 0 (8.64) everywhere along the extremal trajectory. Expanding (8.6l) and (8.62) i n a Taylor series about the nominal trajectory yields-h(t) = h.+Sx - r x S x - r^S\ (8.65) p(t) = p + S\ - i s Sx - :s,S\ (8,66) X A. where h(t) and p(t) are the values of h and p along the neigh-bouring t r a j e c t o r i e s . Substituting (8.63) and (8.6.4) into (8.65) and (8.66) yields S x = A ( t ) S x + B ( t ) Sx - h (8.67) S \ = C ( t ) S x - A T ( t ) S \ - p where A , B , and C are defined i n (8.51), (8.52) and (8.53) respectively. In an i d e n t i c a l manner to that used i n the Min-H strategy, the desired changes i n f i n a l values are swept back to the i n i t i a l . ' point by means of the R i c c a t i transformation (8.48). The f i n a l r e s u l t i s S i = ( A + B P ) S x + B W - h, 0 x Q = 0 (8.69) S\ = P5X + W (8.70) where P = - P A - P B P - A T P + C, P(T) = -0 (T) (8.71) W = , - P B W + P£ - A T W - p,: W(T) = 0 (8.72) Hence, the desired neighbouring trajectory i s x(t) = x(t) +Sx(t) (8.73) \(-t) = \(t) +S\(t) (8.74) where S x and S \ are defined by (8.69) and (8.70) respectively. The process i s continued u n t i l (8.63) and (8.64) are s a t i s f i e d . The characteristics of th i s method are very s i m i l a r to those of the Min-H Strategy i n that the i n i t i a l convergence i s f a i r and the f i n a l convergence i s very good (quadratic). 107. Also, the computing algorithm i s f a i r l y complex and storage i s required at many points along the trajectory. 8.8 The Method of Matching End Points 3„,, 13. In the previous techniques, the i t e r a t i v e search procedure was. performed i n function space and, as a r e s u l t , the desired changes i n these functions had to be computed and stored at many points along the trajectory. The present method i s an example of a technique for which the i t e r a t i v e search procedure i s performed i n a parameter space and storage i s required at the terminal points only. Eor th i s technique, the i n i t i a l trajectory i s determined by selecting a set of i n i t i a l Lagrange m u l t i p l i e r s \(0). The control equation (8.32) i s used to obtain u i n the form of (8.60). The system and Euler-Lagrange equations (8.30) and (8.31) are integrated forward from t = 0 with x(0) s a t i s f y i n g (8.33) and the assumed values for \ . Consequently, the only necessary condition which i s not s a t i s f i e d along t h i s nominal trajectory i s the f i n a l condition (8.34). Expanding (8.30) and (8.31) i n a Taylor series about the nominal trajectory yields (8.75) where C q i s defined i n (8.44). Equation (8.75) can be solved by means of a state t r a n s i t i o n matrix 2> such that (see Section (8.2.1)) Si" ~8*~ •- oo(t) Sx(t) S\(t) a(t,o) (8.76) Evaluating (8.76) at t = T and SXQ = 0 yields S x f = ffi12(T,0)S\Q and where (B 2 2(T,0 )3\ ( 108. (8.77) (8.78) & = 11 12 (8.79) £21 "22 and where the subscript f denotes evaluation of time t^= T. The necessary condition (8.34) can be expressed i n the form E f =,\f + 0 x f (8.80) such that E f = 0 (8.81) for the optimal trajectory. Expanding (8.80) i n a Taylor series about the nominal trajectory and using (8.81) yields E^ . = E, x - = 0 (8.82) J f - T ~ ' v f F x x f ^ A f Substituting (8.77) and (8.78) into (8.83) and solving for §\ yields SxQ = - [ t , 2 f + ^ x f a i 2 f J - l Q i f + ? x fJ • (8.83) As. a r e s u l t , the desired neighbouring trajectory becomes \(0) = iQ + S\Q (8.84) where i s obtained from (8.83). Using (8.84) to define the new nominal trajectory, the procedure i s repeated u n t i l S\ goes to zero. This condition occurs when (A--f+0xf) ± n (8.83) vanishes,, and hence condition (8.34) i s s a t i s f i e d . 109 In the equation (8.83), ^]_2f a n d ^22f m a y k e o b" t ; a i n e d by n forward integrations of the l i n e a r i z e d system of equations (8.75) (see Section (8.2.2)). Alternatively, n systems of (8.75) may be used i n p a r a l l e l , from appropriate sets of i n i t i a l con-d i t i o n s , to provide ^ j ^ f a n <^ ®22f a^"^er o n e ^ o r w a r d integration. In either case, computer storage i s required at terminal points only. The f i n a l convergence properties of t h i s approach are exceptionally good i n the v i c i n i t y of the'optimum; however, the i n i t i a l convergence depends on a good estimate of the optimal trajectory. The computing algorithm associated with t h i s technique i s r e l a t i v e l y complex and can involve much matrix inversion. A modification to t h i s method, proposed by Knapp and Frost [jf], suggests placing the desired end constraints i n a penalty function of the form (8.85) i = l and using a direct approach to f i n d the minimum of P. Comparing (8.34) and (8.85), i t i s seen that, when P attains i t s minimum value of zero, the desired f i n a l conditions of the Lagrange mu l t i p l i e r s are s a t i s f i e d . The technique used to minimize P i s based on a gradient search i n the parameter space of i n i t i a l Lagrange m u l t i p l i e r s . As a r e s u l t , the low memory requirement of the computer i s preserved and the computational algorithm i s r e l a t i v e l y simple. However, as the technique i s 110 based on a gradient method, the f i n a l convergence slows down as the optimum i s approached; and, as the technique i s based on matching the f i n a l values of Lagrange m u l t i p l i e r s , the i n i t i a l convergence depends on a good estimate of the optimal trajectory. I t i s shown i n Section (10.l) that the i n i t i a l convergence properties for these techniques can be improved by selecting an optimal scale factor f or the Lagrange m u l t i p l i e r s such that the error i n f i n a l transversality i s kept at a minimum. 1 1 1 . 9. AN ALGORITHM BASED ON A EIRST VARIATION 9 . 1 Introduction The method to be discussed i s ess e n t i a l l y a combina-tio n of the direct and ind i r e c t approaches. The state, co-tate, and control variables are generated for each trajectory from the necessary conditions developed for the ind i r e c t approach. However, instead of attempting to match end condi-tions, the augmented performance function J i s considered to be a function of the unknown i n i t i a l values of the Lagrange m u l t i p l i e r s , and a direct search for the minimum of J i s c l carried out i n the i n i t i a l conditions space. The result i s a technique which has good i n i t i a l convergence and which i s suitable for a d i g i t a l or hybrid computer of lim i t e d memory. The method also brings out an interesting point concerning the arbit r a r y scaling of the Lagrange m u l t i p l i e r s and helps to explain the d i f f i c u l t i e s , often encountered i n applying the ind i r e c t method, i f the i n i t i a l estimate for the optimal t r a j e c -tory i s not a good one. The a n a l y t i c a l relations necessary for formulating a computational algorithm are derived i n the next section. To avoid unnecessary d i f f i c u l t i e s the control problem i n Section ( 8 . 4 ) i s used. In subsequent sections, i t i s shown that the technique can be extended to problems with bounded control, fixed terminal constraints, and free f i n a l time. 9.2 The Proposed Algorithm 1 8 J Prom (8.29) the f i r s t v a r i a t i o n of the augmented performance function i s 112, S j o = S: X f T 0 x f + S x f T \ f a + SxT(H x +\) + S\ T(H x-x) j dt ( 9 . 1 ) The conditions necessary for an extremum are given by equations (8 .30) to (8 . 3 4 ) . In this technique, the state equations (8.3O), and the Euler-Lagrange equations (8.3 1 ) , are integrated from the i n i t i a l point given by (8 .33) and an assumed set of values \ ( 0 ) . During this integration, the controls are generated according to (8 .60 ) which s a t i s f i e s the necessary condition (8 .32)o The neighbouring tr a j e c t o r i e s are then generated by perturbing the i n i t i a l Lagrange mu l t i p l i e r s by SA.q. The res u l t i n g l i n e a r i z e d equation f or (8.30) and (8.31) are given by ( 9 . 2 ) where C o i s defined by ( 8 . 4 4 ) . The solution of ( 9 . 2 ) with 0 • Sx Ox = c (t) 0 Sx ox Sx Q = 0 i s given i n ( 8 . 7 6 ) to be Sx(t) = B 1 2(t,0)S\„ S \ ( t ) = ffi22(t,0)S\{ ( 9 . 3 ) ( 9 . 4 ) where a?12 and ^ 2 2 are defined i n ( 8 ^ 7 9 ) . The li n e a r i z e d form of ( 8 . 3 0 ) can also be written as Sx = f Sx + f Su Taking the transpose of ( 8 . 3 1 ) « m rn ~X = X f. 'x ( 9 . 5 ) ( 9 . 6 ) i t follows from ( 9 . 5 ) and ( 9 - 6 ) that 113. |^ (V^x) = X Tf u8u = B^Su (9 .7) However, as (8.32) i s s a t i s f i e d along the nominal trajectory, then X f = 0 and (9.7) reduces to d ( ^ S x ! = 0 ( 9 < 8 ) and hence \ T8x = constant = 0 (9-9) where the constant i n ( 9 - 9 ) i s zero since 8x^ = 0 by (8.33)» Substituting ( 8 . 3 0 ) , ( 8 . 3 1 ) , ( 8 - 3 2 ) , (9-3) and (9-9) into (9.1) yields - : S j a = S \ G T ffi12T(T,0) g x(T) (9.10) along the neighbouring tr a j e c t o r i e s generated by this approach. Equation (9.10) i s the desired expression r e l a t i n g the incre-mental change i n system performance to the incremental change i n i n i t i a l Lagrange m u l t i p l i e r s . As mentioned previously, i t i s desired to f i n d that S\ Q which yields the greatest decrease i n J : that i s , which makes 8 J a minimum. Using Schwarz1 ine-a a quality ( x T y ) 2 - (x Tx)(y Ty) and noting that equality holds only i f y i s proportional to x, i t i s seen that i s a minimum when a S\ o =- -k & l 2 T(T , 0)$ x(T) (9.11) where k> 0 i s a constant which determines the step s i z e . Equation (9.11) gives the incremental change i n X which produces the largest decrease i n J . The matrix (EL^ (T,0) 114. can be obtained by n forward integrations of (9.2) (see Section (8.2)) or i t can be obtained by one backward integration of the adjoint system to (9.2) z = - c T ( t ) z (9.12) Solving (9.12) backwards i n time for the state t r a n s i t i o n matrix defined by Z(0) = ^ (0,T)Z(T) (9.13) and p a r t i t i o n i n g SP: 11 12 21 22 (9.14) the relationship W2-^(0,T) = &^_2 (T>0) i s obtained by Lemma 3, Section (8.2). Substituting t h i s relationship into (9.11) yields 8 \ „ = -k^ 2 1(0,T) $ v(T) (9.15) x The vector *' 2 1(0,T)^(T) i n (9.15) can be obtained from one x backward integration of (9.12). with Z(T) = ? X(T) (9.16) 0 The values of Z at t = 0 are defined as Y, Z(0) = •0 '0 (9.17) where the l a s t n components of (9.17) are given by (9.13) to be Z 0 = ^ 2 1 ( T ' 0 ) 0 x ( T ) (9.18) 115. Substituting (9-18) into (9-15) yields the desired change i n the i n i t i a l Lagrange m u l t i p l i e r s S.\Q = -kZ^ (9.19) To insure that the l i n e a r i t y requirements are not viola t e d , a constraint on S\.Q i s imposed of the form S \ Q T S \ o = 8 l 2 (9.20) where S l 2 i s chosen a r b i t r a r i l y small. The desired neighbouring trajectory i s X o = £ o + S x o (9.21) where S\ i s given by (9.19) subject to (9.20). Equation (9.21) becomes the new nominal trajectory and the procedure i s repeated u n t i l S\ Q goes to zero. A comparison of th i s technique with the techniques discussed i n Chapter 8 i s given i n Table (9.1). 9.3 Extension of the Proposed Technique In the previous section, the basic p r i n c i p l e of the f i r s t v a r i a t i o n approach was presented for a somewhat simple-f i e d problem. This approach w i l l now be extended to include problems with fixed terminal constraints and free f i n a l time. The assumptions are s t i l l made that the control i s unbounded and that an e x p l i c i t solution for the control can be obtained i n the form of (8.60). However, i t w i l l be shown i n subsequent examples that these assumptions are not a r e s t r i c t i o n on the a p p l i c a b i l i t y of the proposed method. Consider the problem of determining a control vector. u(t) i n the free time i n t e r v a l 0 ^ t ^ t ^ , so that the per-formance function Table 9.1 Necessary Conditions of the Calculus of Variations Steepest Descent Min - H Strategy Newton -Raphson Matching End Points Proposed Techniques F i r s t Variation Combined Algorithms x z= f(x,u) X X X i X X • -3* T X = -g^ \ = h(x,\,u) X X X X X 3H 3 f T , ' du du - X X X X X X X Sx T(0)\(0) = 0 X X X X X X x(0) = x Q X X x X X X g ( x ( t f ) , t f ) = 0 X S ( x ( t f ) , t f ) =0 X X X X X X Search Over u(t) u(t) x(t) x(t) \(o) \(o) Search C r i t e r i o n S j a < o 3H -du g-0 (±-f)->0 (\-h)^0 * g-*o S J a < ° g-0 § J a < 0 or X f - A . J g-*o Computer Storage along trajectory along trajectory along trajectory at end points at end points at end points I n i t i a l Convergence good f a i r f a i r poor good good F i n a l Convergence poor good good good poor good 117. J = 0 ( x f , t f ) (9.22) i s a minimum. The state vector x i s subject to the constraints X = f(x,u) (9.23) x(Q) = X . 0 (9.24) g ( x f , t f ) - 0 (9.25) S ( x f , t f ) = 0 (9.26) where x i s an n x 1 vector of state variables, f i s an n x 1 vector of continuous functions, g i s a p x 1 vector of terminal constraints, and S i s a scalar function of terminal values. The terminal constraints (9.25) and the dynamical constraints (9.23) are adjoined to the system performance function (9.22) to y i e l d t J a = 0 ( x f , t f ) + g T ( x f , t f ) ^ + J (\Tx - H)dt (9.27) 0 where H = X f(x,u) and V i s a p x .1 vector of constant Lagrange m u l t i p l i e r s . The function (9.26) i s used as a stopping function to define t ^ . Taking the f i r s t v a r i a t i o n of (9.27) yields 8 j a =S0(x f,t f) + S ( g T ( x f , t 1 ) y ) (9.28) +J (-SH +8\TX + X^S^dt + (\Tx - H)Std 0 where S0 f = d x f T 0 x + 0 t8t f (9.29) 118. S(g Tl^) = dx f ! ISg x T V + sll&-f v (9.30) a X f = 8x f + x f S t f (9.31) SH = Su TH U + Sx TH X + S\TH^ (9.32) t f t f J ~ x^sidt = S x V f * - ^ T S x T \ d t (9.33) 0 0 Using (9.26) to define S t f yields Ss = d x f T S x + S t S t f = 0 (9.34) After substituting (9.29) to (9.34) i n (9.28), the expression for SJ i s a S j a =8x f T ( 0 s f + g s f V ) - J ( S u \ + S x T ( H x + i ) 0 + S\ T(H X - x))dt +8xVf' f + \ f T ( x - f ) f S t f (9.35) where As i n the previous case, the necessary conditions (8.30), (8 .31), (8.32) and (8.33) are s a t i s f i e d along the nominal trajectory and X^Sx = 0 along a l l neighbouring paths. As a re s u l t , using (9.3) to define 8Xg, equation (9.35) reduces to 119. (9.38) for the neighbouring t r a j e c t o r i e s . In a manner i d e n t i c a l to that of the previous section, i t can be shown that the minimum of S j occurs when a S\0 = -k ( z 0 + Z gp) (9.39) where A £ 0 - ffl12f )-sf A * A T Z - 3> T a-- ~ a12f g s f ' The .vector Z 0 i s found from the ,last n components of Z(0) when Z(T) = th sf (9.40) A and the i column, Z . , of the matrix Z i s found from the g i g l a s t n components of.'.Z(O) when hat Z.(t) =' where g^ i s the i ^ * 1 component of g. In general, the boundary constraints (9.25) w i l l not be s a t i s f i e d along the nominal trajectory. I f g were small for the nominal trajectory, a SXq could be chosen such that Sg = -g.. Thus, X. + § ^ 0 would result i n g + Sg = 0 as required. However, th i s choice could v i o l a t e the l i n e a r i t y requirement and thus, i n order to keep the error small, 120 g = -ag (9.41) i s chosen where a i s an arb i t r a r y small positive quantity, 0 - a - l . Equation (9.41) imposes a constraint onS\ Q which can be determined from the incremented equation Sg = ( g x - (S )§ / ) f 8 x f (9.42) S Substituting (9-3), (9.41), and using the d e f i n i t i o n of Z from (9.39), equation (9.42) becomes Sg = Z g TS\ Q = -ag (9.43) To insure that the li n e a r i z e d equations are v a l i d , a further constraint S\ 0 TS\ 0 = S l 2 (9.44) i s imposed where S l 2 can be chosen a r b i t r a r i l y small. The evaluation of V, k and SXq i s carried out i n Appendix G and yields S i 2 - . ^ ( 0 / 8 -\ —m a r - * — * m T—m— (9.45) 9.3.1 An Algorithm for Numerical Computation An algorithm for numerical computation based on (9.45) can now be formulated. ( l ) . Select an i n i t i a l \ and integrate (8.30) and A (8.31) from (8.33) at t = 0 u n t i l S = 0, which defines t ^ . During t h i s integration (8.32) i s used to define u. 121. (2) With the aid of the nominal trajectory determined i n step ( l ) , compute (9.40) and (9.41). Integrate (8.30) and (8.31) and (9.12) backwards p + 1 times and f i n d ZA and Z . ' 0 g (3) Select a and S l 2 and calculate O l 2 - a 2g T ( Z TZg)~xg. I f t h i s quantity i s less than zero, adjust a to make i t zero. I f th i s quantity i s greater than zero, no change need be made i n a. (Note: 0 ^ a - l ) (4) Compute S\ Q using (9.45). (5) Select a new trajectory using \ + 8A . q and repeat steps ( l ) to (5). A considerable s i m p l i f i c a t i o n i n the computation can be achieved i f . the. end constraints are considered i n a penalty function of the form P g = ) K ± g i 2 (9.46) 1=1 where the IL are assigned weighting factors. In th i s case, Z i s a vector and can be determined by one backward integra-t i o n . In theory, the weighting factors approach i n f i n i t e values as the g^ go to zero. However, i n numerical computa-tions, i t i s not possible to have i n f i n i t e values for the and, as a result,,spurious extremals may be introduced. Por most cases reported, however, th i s characteristic has not caused any r e s t r i c t i o n on the use of (9.46). It can be noticed that t h i s proposed algorithm re-quires a minimal amount of computer memory since only the i n i t i a l and f i n a l values need be stored. The elements of 122. T C Q can be determined by function generation using the values of x and X found from '(8.30) and (8.31). 9.4 The Arbitrary Scaling of the Lagrange M u l t i p l i e r s In the c l a s s i c a l theory, the transversality condition for the problem of Section (9.2) i s 0 (9.47) t=t f 0 (9.48) and ( - f \ * + 0 t + g t T^+ S ti7) t=t f for an unspecified t ^ . However, i f S,- i s used to define t ^ , then the f i n a l transversailty conditions can be modified. Solving (9.17) and (9-48) for V Q yields VL = - i (0 + g T ^ ) ~ (9.49) 0 S 1 Substituting (9.49) into (9.47) yields Xf*'= ( 0 s f + g s f T ^ } ' ( 9 ' 5 0 ) where 0 g f k ( 0 x - *x4))t S , T A / T _ ( £ T ) } = Sf - K S X ^ X V g U f •st-and where X i s the Lagrange m u l t i p l i e r s for the c l a s s i c a l case. In the present approach, however, i t i s seen from (9.38) that a trajectory i s sought for which any increment 8A.q results i n S j a =Sx f T ( 0 s f + g s f Tp) = 0 (9.51) Hence, the vector (0 - + g ~ V) i s normal to the hyperplane SI S X 123. formed by the Sx^. From (9.9) S x f T \ f = 0 (9.52) and hence X^ i s also perpendicular to the v a r i a t i o n Sx^. Since Sx^ i s an a r b i t r a r y vector i n the hyperplane, i t can be concluded from (9.51) and (9.52) that X^ i s colinear with (0O-P + g f );) for the optimal trajectory. Hence at optimality SI SI x f = V(0sf + S s f V ) (9.53) Note that by (9.50), \ f = \ f* only i f u. = -1.' This choice for u i s actually an unnecessary re s t r i c t i o n . . Since (8.31), (8.32), (9.47), and (9.48) are l i n e a r i n X, i t i s seen that i f X results i n a trajectory which extremizes 0 with * .... g = S = 0, then, w i l l : r e s u l t i n a trajectory which extremizes |i0 with ug = |i;S = 0. However, the control variable u i s the same i n both cases and i t i s u which i s desired. Thus, instead of searching for a point X , as i s done i n the c l a s s i c a l approach, i t i s s u f f i c i e n t to search' for the l i n e XQ = \xX (see Figure 9.1). It i s shown i n Appendix H, that a value of J i s associated with each r a d i a l l i n e through the a 0 o r i g i n i n the X -space. I f the i n i t i a l estimate X results 0 o r o 1 i n J = J > J . , a search can.be carried out over the sur-a a mm' face of a sphere for the minimum J (J . = 0 . ) which occurs when the l i n e a amm. ^mxn. * \iX intersects the sphere. As the search i s conducted over a sphere, i t i s convenient to define S l 2 i n (9.44) to be 124. 125. S l 2 =S<x 2\ o\ (9.55) where S o c i s the angular incremental rotation of the X vector r e s u l t i n g from the incremented displacement SA.q. For t h i s choice of S l 2 i t i s seen i n Appendix H that f o r - \ = c\ , where c i s any non zero scalar, then Sx. = cSx results from o o (9.45). Hence the speed of convergence for t h i s technique i s independent of any i n i t i a l scale factor. I t i s evident from t h i s result and Figure 9.1, that the i n i t i a l convergence i s not dependent on a good f i r s t estimate of X . 9.5 Example 1 Consider a system of the type i l l u s t r a t e d by the state t r a n s i t i o n flow graph of Figure 9.2a. The reaction k i n e t i c s are x l + k " l x l = ^ ^1 "^ 1 ^2 "^"15 o * o *- -o (a) Figure 9.2 (a) Transition flow graph for example 1. (b) Transition flow graph for example 3. 1 2 6 . , A -E,/RT k x = G-^e . 1 ' k 2 k & 2 e - E 2 / R T where the absolute temperature T(t) i s to be found such that x 2 ( t ^ ) i s a maximum, where t ^ i s fixed and where E. ^ = 18,090 cal/mole t f = 8 min E g = 3 0 , 0 0 0 cal/mole R = 2 cal/mole/° Q± = 0 . 5 3 5 x 1 0 1 I L/min G-2 = 0 . 4 6 1 x 1 0 1 8 / m i n Por convenience, the control variable i s taken as u = k^ and T i s considered unbounded. The matrix -C i s (see ( 8 . 4 4 ) and ( 9 . 1 2 ) ) -C T = u 1+p -px 1/x 2 -.px1/ (A . 2 - A . - L ) -(l+p\ 1/\ 2) p. ( ^ 2 "X 1 ) x. i-d+p^/X,) -p 1 (Xg-X-^ x. (X.2-\-L) x 1 X, U ^ ) X 1 A 2 " ( L + P ) px-^/X,(X,.^) 2 l+pX-^X, ^ ( where p ='E^/(Eg-E^) and n = Eg/E^. An i n i t i a l estimate \ 2 Q = 1 . 0 and X ^ Q = 0 . 1 was a r b i t r a r i l y selected. Table 9 . 2 i l l u s t r a t e s the numerical r e s u l t s . Note that the c l a s s i c a l theory requires that X-^ = 0 , X2£ = 1 . 0 (since 0 =-x^). The f i n a l values f or the i n i t i a l trajectory are - 8 2 1 0 and 1 3 3 , respectively, and are grossly i n error. However, after f i v e i t e r a t i o n s using ( 9 . 4 5 ) and the proposed'algorithm, i t i s 'X -X-, pX-Table 9.2 3TBP . F i r s t Variation (9.45) Combined Algorithm with F l Combined Algorithm with F2 X l f x 2 f ht X2f x 2 f * l f X2f 0 -8210. 133-3 0.02005 -8210. 133.3 0.02005 -8210. 133.3 0.02005 1 -7407. 122.3 0.02144 -1019.0 26.58 0.072622 -1019.0 26.58 0.072622 2 -1825. 40.58 0.05160 -64.92 4.882 0.274662 -64.92 , 4.882 0.274662 3 -317. 12.10 0.1348 -0.2057 1.103 0.679878 -0.2057 1.103 0.679878 4 -33. 3-52 0.3472 ::0. 01645 1.002 0.681691 -0.003409 0.9253 0.681707 5 -0.874 1.214 0.6626 -0.0001526 1.000 0.681707 -0.000005 1.0000 0.681707 6 -0.3896 1.123 0.6761 0.0000004 1.000 0.681707 7 -0.1665 l."075 0.6304 8 -0.6952 1.053 0.6815 9 -0.02866 : 1.043 0.6817 10 -0.01174 1.039 0.6817 11 -0.00930 1.037 0.6817 12 -0.00196 1.037 0.6817 13 -0.00079 1.036 0.6817 128. seen i n Figure 9.3 that the Lagrange mu l t i p l i e r s have converged # very close to l i n e \i\ . As the optimum i s approached, the 01 0-2 03 0-4 0-5 W 0-7 0 8 OS 10 Figure 9-3 Iter a t i o n Path i n the I n i t i a l Conditions Space of Lagrange M u l t i p l i e r s rate of convergence of the gradient method slows down and i t i s not known when the search should be terminated. In the next chapter, methods of improving t h i s f i n a l convergence are discussed. The f i r s t approach i s based on the method of matching end points i n which the Lagrange mul t i p l i e r s are continually re-scaled to maintain a minimum error i n f i n a l 129. transversality. As a r e s u l t , i t i s shown that t h i s method of matching end points can be used i n the f i n a l stages of the proposed technique to provide, the property of rapid f i n a l convergence. Another technique discussed i s a second v a r i a -t i o n method which determines the optimal step size for SA.Q as the extremum i s approached. The effect of using these two techniques i s shown i n Table 9.2.J Note that for these l a s t two methods, the solution converges to the c l a s s i c a l solution, yet for the f i r s t v a r i a t i o n approach,' the scale factor |i = 1.036 r e s u l t s . Figure 9.4 and Figure 9.5 i l l u s t r a t e the T and x 2 p r o f i l e s for various i t e r a t i v e cycles. The i t e r a t i o n path i s i l l u s t r a t e d i n Figure 9-3> and Figure 9.6 i l l u s t r a t e s the time variations of x.^ x 2, X-^, \ 2, and T for the optimal path. To plot these quantities on one graph, the following ordinate scales are used: y = (T-326)/2 y = 10 (-^ + l ) y = 10x 1 y = 10(-\ 2 + 1) y = 10x 2 9.6 Example 2 To i l l u s t r a t e how the; proposed method can be used i f u i s bounded, consider the equation of constraint (u-u . )(u -u) = a 2 (9.57) x mm.' v max ' Due to t h i s constraint, (8.32) has the form 0 = H + \-(u Q Y + u . - 2u) (9.38) u 3 max. mm. 7 2-0 40 SO 80 TIME (MINUTES) Temperature for Various Iteration Cycles TIME (MINUTESf Figure 9.6 The time variations for the optimal path (Ex.1) Figure 9.7 Temperature p r o f i l e s for oases a to d (Ex .2 ) 132. where \^ must be introduced because of (9.57), and i s determined by 0 = \5<x (9.59) Thus, when u = u . or u = u , a = 0 and (9.59) i s s a t i s -' mm. max. f i e d . I f u . < u<u , \, = 0 and (9.58) reduces once mm. ^ max. 3 more to (8.32). I f 0 i s to be a minimum, the Legandre Clebsch condition yields 2H (§u2 + S a 2 ) u u +u . -2u max. mm. ^ 0 ( 9 . 6 0 ) and hence H < 0 when u = u , H > 0 when u = u . and H = 0 u max.* u ^ mm. u i f u . <u<u mm^ ^ max. Consider the case of Example 1 where T has the following upper bounds: (a) T m = 345 (b) T m = 342 (c) T m = 340 (d) T m = 338 In t h i s problem i t i s known that u = u _ for 0 ^ t.<t . and \ msix c s u<u, for t < t ^ t o . The instant t i s determined when H max. s ^ f s u vanishes. After t h i s instant, u i s computed as i n the unbounded case. F o r ' t ^ t , the elements i n C (see ( 8 . 4 4 ) ) s o change, since ir£ = §^ = 0 when u = u . However, this change o x o A. max • i s r e a d i l y carried out i n the backward sweep by storing the value of t found during the forward sweep. With t h i s s l i g h t s modification, the computation i s the same as i n Example""'1. The temperature p r o f i l e s for cases (a) to (d) are i l l u s t r a t e d i n Figure 9.7. In Figure 9.8, the ordinate scales "used are: y = T - 331 y = -j-^ + 10 134. 10\, 100H + 10 y = Y y = 5 (" x2 + 2 A ) y = 10x 2 y = 10x 1 9.7 Example 3 It i s not always possible to use (8.32) to obtain an e x p l i c i t a n a l y t i c a l expression of the form (8.60) for u as a function of x and \. This would complicate the computa-ti o n since (8.32) would have to be used to determine u i m p l i c i t l y . Note, however, that (8.32) yields the l i n e a r i z e d equation Su = - ^ ( H ^ S x + Ejx) (9.61) To avoid the i m p l i c i t computation of u, note that (8.32) can be differentiated with respect to time to y i e l d u = -H - 1 ( H x + H , \) (9.62) uu v ux u\ where the right hand side i s a.function x, \, and u. The i n i t i a l value of u Q can be found by an i m p l i c i t solution of (8.32) at t = 0. With u Q known, (8.30), (8.31) and (9-62) can be used to compute.u for the nominal trajectory during the forward and backward sweeps. The procedure i s otherwise the same as before. To i l l u s t r a t e t his modification, consider the batch reactor problem of Example 1 where there i s an extra unwanted by-product x^ (see Figure 9.2b). The equations are x l = - ( k i + k 3 ) x i (9.63) 135, Let u = k^. In this case (8.32) has the form , n,-a-;: , n 2 - l H u = X-^x^d+YijGcj u D .;) + \ 2 ( " x i + n 2 G ' 2 X 2 U ^ where n 2 = \ / \ , k , &2' = 2 , ( 9 . 6 4 ) G ' = G 3/G ± n3. Hence u = N/M , i ( 9 . 6 5 ) where N = -(\^(l+n^G^ u "') - ^ 2 ) x 1 - \ 2 n2^2 u x 2 , n,-i . , ' n 2 - l . -x 1(l+n^G^ u J ~ (-x-L+n2&2 x 2u ) \ 2 x ' n 3 " 2 / \ ' n2"2 M = X.1x1n^(n^-1) G^ u ^ " + \ 2 n 2 ( n 2 - l ) G 2 u x 2 Figure 9 .9 i l l u s t r a t e s the computed results for the case where B 1 = 18,000 cal/mole, G1 = 0.535 lO^/min E 2 = 30,000 cal/mole, &2 = 0.401 10 1 8/min = 27,000 cal/mole, G^ = 0.500 10 1 6/min The ordinate scales used i n Figure 9 .9 are y = T-330, y = 10(1-^) y = 10x 1, y = 10(1.l-\ 2) y =liOx2., 1 3 6 . 1 0 . TECHNIQUES FOR IMPROVING- PINAL CONVERGENCE 1 0 . 1 Matching End Points Using, an Optimal Scale Factor J20 Consider the optimal control problem i n section ( 8 . 4 ) which i s to determine the control u(t) over the fixed time i n t e r v a l 0 - ^ t ^ T such that the system performance function J = 0(x(T)) ( 1 0 . 1 ) i s a minimum. The constraints on the state variables are x = f(x,u) ( 1 0 . 2 ) x ( 0 ) = x Q ( 1 0 . 3 ) As i n section ( 8 . 4 ) , the constraint ( 1 0 . 2 ) i s adjoined to the system performance function (lO.l) by an n x 1 vector of Lagrange m u l t i p l i e r s \. However, since i t was shown i n Chapter 9 that the scale factor for the Lagrange m u l t i p l i e r s i s a r b i -trary, the augmented performance function i s taken as T J a = 0(x(T)) + H / U T x - H)dt ( 1 0 . 4 ) 0 A T where H = \ f and where \i i s introduced as the arbitrary scale factor. For th i s case, the transversality condition ( 8 . 3 4 ) i s uX^ + 0 x f = 0 ( 1 0 . 5 ) Substituting ( 1 0 . 5 ) i n ( 8 . 8 0 ) yields E f = $ f + $xf ( 1 0 . 6 ) where E^ = 0 on the optimal trajectory. Expanding ( 1 0 . 6 ) i n a Taylor series about the nominal trajectory yields SE„ = JJS\ . + &~ 8 xxf ^ X f 137. Using SE^ = -E^ , and substituting for Sx^ and from (9 .3) and (9.4-) yields Sx o = - [ jA, 2 f + 0 x x f fi12fj + ? x f J (10.7) which d i f f e r s from (8 .33) by the scale factor u. As Sx IS 0 proportional to E^, p, i s selected such that the square of the error i n f i n a l transversality i s a minimum and hence SA. Q^SA. 0 i s minimized for each i t e r a t i o n . Erom ( 1 0 . 6 ) , the square of the error i s ,2 TC , TS ft T * , A o A r r i A A r n A A r n A A m A |Bf r = E f X E f = u. \ f x \ f + 2 i i X f 1 0 x f + 0 x f x 0 2 (10.8) D i f f e r e n t i a t i n g (10.8) with respect to \x and equating the derivative to zero yields 2p\ f x\ f + 2 \ f x 0 x f = 0 (10.9) Hence the value of \x which minimizes the f i n a l error i n trans-v e r s a l i t y i s ^ p t = - f r f <io-io> f f Substituting (10.10) into (10.7) yields the desired incremental change -1 ^J^Kf A A A 1 r~ A T A - 2 2 f T ) M x x f « 1 2 f A.^ A.^ A-f ) " x - p A — — + * T(v f ~ / x f f f (10.11) It i s shown i n Appendix H that, using (10.11) to define §A. Q , the rate of convergence i s independent of the i n i t i a l scale factor, and that the procedure converges to the solution l i n e 138. . |JA O This fact i s i l l u s t r a t e d i n Figure 10.1, for the problem of Example 1 i n Section ( 9 . 5 ) . I t i s seen that a "cone of convergence ", region RQ, exists about the solution l i n e |i\ . Within R , any i n i t i a l estimate for A. w i l l converge to the solution l i n e independent of the i n i t i a l scale factor. The number of steps required for convergence i s indicated on the r a d i a l l i n e s . For an i n i t i a l estimate which f a l l s outside the cone of convergence, region R , an unacceptable trajectory nc r e s u l t s . For t h i s p a r t i c u l a r example, matrix C Q i n (9.56) contains terms which are divided by (X^-Xj) and which become i n f i n i t e i f X^ = X^. Therefore, should t h i s s i t u a t i o n e x i s t , the trajectory i s unacceptable. To f i n d an i n i t i a l estimate which l i e s inside the cone of convergence, a random search can be ejnployed or a method of relaxing the f i n a l constraints, as done by Isaacs et a l |l5J , can be used. I f i t i s desired to have the; procedure converge to the c l a s s i c a l solution XQ , then i t i s shown i n Appendix H that the neighbouring trajectory should be taken as xo = %A +s\>} (10-l2) A where XQ i s defined by the nominal trajectory, \i-Q^ i s defined by (10.10) andS\ Q by (10.11). The effect of using (10.12) for the solution of Example 1 of Section (9.5) i s i l l u s t r a t e d i n Figure 10.2. For any i n i t i a l estimate of X which l i e s inside the region R , the f i r s t step i n the i t e r a t i o n estab-° ca * li s h e s the i n i t i a l conditions on the solution curve C . a Once on this solution curve, X moves along the curve u n t i l X i s reached. I f the i n i t i a l " e s t i m a t e l i e s i n region R , , O CD 139. Figure 10.1 The Solution Line u\. and the Cone of Convergence R Q 0-8 Figure 10.2 Convergence to X by the Modified Methoi of Matching End Points 141. the f i r s t step i n the i t e r a t i o n establishes the i n i t i a l condi-tions on the l i n e C, which i s i n region R . The procedure D ca i s then the same as before. For the technique developed by Znapp and Frost, (Section (8.8)), a si m i l a r improvement i n i n i t i a l convergence can be obtained when the penalty function i n (8.85) i s replaced by n p=y~K ± 2(\i\ ± t + 0 x i f ) 2 (10.13) i = l In t h i s case u i s selected to minimize P and hence the optimal value for u i s / (K. 2 \.« 0 . J u , = - ( l 0 . i 4 ) 0 p t V ^ K i ^ i f ) Using (10.14) i n (10.13), i t can be shown that the technique i s then independent of the i n i t i a l scale factor, and hence the rate of i n i t i a l convergence has been improved. 10.1.1 Extension of the Method of Matching End Points Consider the control problem i n Section (9.3) and l e t the augmented functional be T I T ' J & = 0 ( x f , t f ) + g (x f,t f)i/> + \x J (X x - H)dt 0 (10.15) where \i i s the scale factor for the Lagrange m u l t i p l i e r s . As before, the f i n a l time i s defined by S ( x f , t ) = 0 (10.16) 142. and the variation in f i n a l time for the neighbouring trajectory-is defined from (10.16) to be St- = - S x / (10.17) S f Using (10.17) to define St^, the f i r s t v a r i a t i o n of J i n i a (10.15) i s S j a = S x f T ( 0 g f + g s f \ 0 + pSx f\ f + (X Tx - H) fSt f •(\ + I \ ) j d t (10.18) where 0 s f £ 0 r f - 4)f and * T ^ * T - S (£"") na g g f _ g x f E > x f A. ; F s The transversality condition for t h i s case i s uA f + 0 s f + g ^ y =0. • (10-. 19) Therefore, for the nominal trajectory, the error i n transver-s a l i t y i s defined as Ef. = uA f + (# s f + g s f T^) (10.20) Expanding E f i n a Taylor series about the nominal trajectory yields §E f = pS\f + 0 s s f Sx f + gggf^Sxf + g s f TcSp (10.21) where 143. e S S f v = ( g x x i> - — — + — p — ) f 3XX = ( x 1 + 2X 1 g x t V + g ^ V ) / g >tt Substituting for S x f and S \ f from (9-3) and (9-4), and using S E f . = - E f yields (10.22) where A A = ^ + ( l s + S s s ^ ) J 1 2 Solving (10.21) for SX q yields ^ = - V 1 < § > + i f ) - - , 0 . (10.23) 0 A. °s w i Prom (9.25), i t i s required that g be zero on the neighbouring trajectory. Expanding (9.25) i n a Taylor series about the nominal trajectory and using (9.3) yields S g = g„ o S s - 1 2 ^ 0 (10.24) Substituting .(10.23) into (10.24) and using S g = -g yields -khz v 1 ( § > + v = -g Equation (10.24) can be solved for to provide (10.25) A A -1 A T S S ffi12 *\ g s g " g s ffi12 V -1 'f (10.26) Replacing i n (10.23) by (10.26), the desired incremental change becomes 144. 0 0 L o l (10.27) where 8 -1 A T X = g r\r\ A o 0 0 X s s f «sf ffi12f ffl\ g s f . " I A T / A A \ g s f ^ g s f "12f "X 5 s f -1 A T\ A A £ g o H, (ga.p ^ g ^ ) gg f ffi12f \ -1 -1' Note that for SA.q-= §\ , equation (10.24) becomes -1 £ • A A -1 A T ,-A £ - -1 A T-| A . 6 g = " g s f ffi12f \ g s f p s *12 \ g s ' ] f g f A (10.28) and hence §\ q o i s that component of SX q which attempts to' s a t i s f y the end constraints (9.25). Using&X Q =SXQ-^ i n (10.24) yields " Sg = '* * -1 A T g =: gsf ffi12f \ g s f g s ffl12 ffl\ -1 A T -1 ^sf *12f - 1 A A A - 1 A E f " g s f * 1 2 f \ E f \= 0 ( 1 0 . 2 9 ) and hence SA.q^ i s that component of &XQ which attempts to s a t i s f y ' the transversality condition without affec t i n g the end constraints set byS\ 0 0. Equation (10.27),. therefore, provides the desired incremental change i n X .. However, to evaluate (equation (10 . 2 2 ) ) , values of \x and V are: required. As i n the previous section, these values are selected to minimize the square of the error i n transversality. By (10.20) th i s error i s 145, JE f \< = E / E f = 0 r f- 0 s f + 2 P T Vj£ 0af + ^ q£ ^ r-,T ,TIT. (10.30) where mT A T "sf ' A f and D i f f e r e n t i a t i n g ( 1 0 . . 3 0 ) with respect, t o y , and equating the derivative to zero yields , / 2 % ^sf + V ^ 0 ( 1 0 . 3 1 ) The optimal values for the augmented vector V i s found from ( 1 0 . 3 1 ) to be-',' opt - 1 ^ s f ( 1 0 . 3 2 ) As a r e s u l t , the desired neighbouring trajectory becomes V = ^opt A +S\>) ( 1 0-33) where SA. i s defined by (10.27), u„ + and V^, are defined by o opu opu . . A (10.32), and \ i s defined by the nominal trajectory. 1 0 . 1 . 2 Computing Algorithm F l using Matching End Points ( 1 ) From \ Q, which defines the nominal trajectory, and ( 8 . 3 3 ) , integrate ( 8 . 3 0 ) , ( 8 . 3 1 ) and ( 9 . 3 2 ) from appropriate i n i t i a l conditions u n t i l S •= 0 , which defines t ^ . (Note: n systems of ( 9 . 3 2 ) are i n t e -grated i n p a r a l l e l from the i n i t i a l conditions given i n Section ( 8 . 2 . 2 ) ) . ( 2 ) Test |8jal<£-[_ and/or |g exi"k> where e-^ and £g are chosen to provide the desired degree of accuracy. 146. (3) Compute \i t, v %t andSxQ by (10.33) and (10.27) respectively. (4) Replace XQ by (10.33) and repeat from ( l ) . 10.2 Determining Optimal Step Size for the F i r s t Variation Approach In Chapter 9,. the incremental change i n \ based on a f i r s t v a r i a t i o n approach was shown to be (see (9.45)) SX q = a S \ Q 0 + kS\ o l (10.34) where Sx = -2 (Z TZ ) _ 1 g w oo g g g s S\ n = \Z (Z TZ )" 1Z TZ f l < - ZA ^ o l g g g g 0 0 W z 0 - z 0 T z g ( z A ) " l z g T z 0 and O ^ a ^ l It was further shown that since the r e s u l t i n g search i s carried out over the surface of a sphere, i t i s convenient to define S l 2 as (see(9»55)) S l 2 = S & \ H 0 (10.35) where Sa i s the angular incremental rotation of the X vector r e s u l t i n g from the incremental displacement SXq. The vector Sx i n (10.34) i s a l i n e a r combination of the vectors SX q 0 and S x Q l . It i s shown i n Appendix F that SX q 0 i s the component of SX q which attempts to s a t i s f y the end conditions (9.25), and that SXq-^ i s that component of SX q which attempts to mini-mize Ja without affecting the end conditions set by SX O Q. The r e l a t i v e emphasis placed on these two effects depends on 147. the values selected for the parameters a and k. Prom (10-34), i t i s seen that these values are subject to the constraints St XQ \ (1-a /b ) B S " f« — T fn (10.36) ZA~ ZA - ZA L Z (Z 1 Z ) x Z ZA 0 0 0 g g g g 0' where b = g A 2 A T A o o (10.37) As a r e s u l t of (10..36) and (10.37), two sets of values can be assigned to a and k depending on the value of b.l!. The two cases are: (1) for b ^ l , a = b and k = 0 (2) for b > l ,-. a = 1 and k i s defined i n (10.36) For the f i r s t case, the error i n the f i n a l end conditions i s large and f u l l emphasis i s placed on minimizing t h i s error. For the second case, the error i n end constraint i s small and, within the step Set, i t i s possible to reduce Ja while also s a t i s f y i n g the desired end conditions. To take advantage of the g o o d - i n i t i a l convergence properties of t h i s gradient tech-nique, a three stage algorithm can be developed that uses the gradient technique i n the f i r s t two stages to rapidly locate the region of the optimum. However, as the f i n a l convergence of the gradient method i s r e l a t i v e l y poor, a t h i r d stage i s used which has good f i n a l convergence properties and which can be , used with the gradient method. An example of such a technique i s the modified method of matching end points of Section 10.1. Also, i n the next sections, two techniques are developed which 148. determine'the optimal step size for the gradient method as the optimum i s approached. The f i r s t i s based on a second va r i a t i o n approach, and the second i s based on a method of curve f i t t i n g . Combining these techniques, therefore, the three stage algorithm i s as follows: (see Figure 10.3) (a) F i r s t stage: i f b<^l, a and k are defined by-Case ( l ) and hence K = ^ 0 0 = This region i s called the Rg region since f u l l emphasis i s placed on s a t i s f y i n g the end conditions, g = 0. The search i s carried out over the surface of a sphere with a constant rotation Sa u n t i l b >1 which defines the second stage. .(b) Second stage: i f b > l , a and k are defined i n Case (2) and hence S\rt = S\ + kS\ n 0 0 0 ol This region i s called the RJa region since the emphasis i s now on reducing Ja without affec t i n g the end conditions set by Sx . The search continues ° oo on the surface of the sphere with a constant rotation Sa u n t i l Ja increases which defines the t h i r d stage, (c) Third stage: At the point before Ja increases, i t can be concluded that the rotation Sec was too large and the region of the optimum was overstepped. This region, called the Ret region, l i e s i n the i n t e r i o r of a cone-shaped surface which has a maximum angular width Oct and which contains the solution l i n e \xk . CONE OF CONVERGENCE .A 02 10. "3 The Regions Rg, RJ , and Ra Abou.t the Solution * a L i n e U A 150. Within this region one^ of the techniques with good final convergence properties is used to complete the search procedure. 10.2.1 The Second Variation Method of Determining the Optimal Step Size Taking the variation of Ja defined in (10.15) and keeping all terms up to second order yields S J & = S x F T ( 0 s f + g s f a v ) + & F T ( 0 8 8 f + S S s f T ^ x f ,.. " fjT(Su TH u u8u+^Su T5 u xSx -r8xTHx x8x4SNTHX^S\)dt 0 (8\ T(i-H x) +8xT(\+Hx) +:8uTHu)dt + [\ Tx-HJ f8t f + n^/8\ T(Sx-H X xSx-H X l^u)dt f 0 + Sx T\ f (10.38) ^ 0 where SsfT» 0AF> £Bef^ ^ssf a r e a e f i n e d i n (10.18) and (10.21) respectively. Using the first variation approach developed in Chapter 9, the following relations hold for the nominal trajectory: (1) x - = 0 , from (8.30) (2) X + H = 0 , from (8.31) •A. (3) Hu = 0 , from (8.32) (4) ,H X X = 0 , from (8.24) 151 (5) S i - H^Sx - H^uSu = 0, from (9-5) (6) x(0) = x Q , from (8.33) (7) S x f \ f =-0 , from (9-9) Substituting ( l ) to (7) into (10.38) yields S j a = S x f T ( 0 s f + g s f T y ) - ^ f T ( 0 S 8 f + g s f f T y ) S ^ f " \ J ( S u \ u S u + 2 S u T H u x S x 4 S x T H x Sx)dt 0 (10.39) Prom (9.61) the v a r i a t i o n i n u for the neighbouring trajectory i s given by Su = -E - 1 (H Sx + H .S\) (10.40) uu ux u\ Por Sx = 0, the value of Sx(t) and S\(t) can be obtained from o (9.3) and (9.4) i n terms of the submatrices ffi^ and ®22' Using (9.3), (9.4), and (10.40) i n (10.39), S j becomes > S j a = S C ( Z 0 + Zg V ) + ^ o T ( ' a i 2 f ^ B B f ^ B B A ) S 1 2 f ^ o " %W[f ( D \ u D + 2 D T g U x % 2 A 2 \ x B 1 2 ) d t l S ^ 0 o (10.41) A A T where Z 0 = fl?12f ^ • „ A £ T A T Z g = *12f g s f a n d D = " i u " 1 ( Sux °12 + *uX ®22> By the f i r s t v a r i a t i o n approach, equation (10.34), the value of S\n that minimizes J subject to the constraint Sg = -g V cl 152. (a = l ) i s S\ = S\ n + kS\ , (10.42) O 0 0 o l By (Gf-2), the value for y-was found to be V = ^ + \ (10.43) where V = (Z TZ J " 1 ^ 0 g g ^ = - ( z T z ) _ 1 z T z M 1 g g' g 0 For. t h i s problem, the error i n f i n a l t r a n s v e r s a l i t y i s , B f = u\ f + (^ s f+g s f T^)' (10.44) Substituting (10.43) into (10.44), and determining a as that T value which minimizes E^ E^ yields H = ^ + ^ (10.45) where u Q = - \ f g g f X f ^ I = " x f ( 0 s f + g s f * i ) A f x f For this approach, i t i s desired to fi n d the value for k i n (10.42) which causes §J to be a minimum. Substituting (10.42), (10.43) , and (10.45) into (10.41) yields SJ = f S\ TW.S\ + US\ TW,S\ ^ T SA ) a k 0 0 4 0 0 0 0 4 o l 0 0 2 0 0 2 + t ( S \ 0 l \ S \ 0 l - f S \ 0 l T w ^ S \ 0 0 ) + | - S x 0 l T v 5 S x 0 l (10.46) where W]_ = 1 £ 1 2/(£ s s f + gssf V * s s f ^ s s f "1' 12f A . A <n , A rn A i-p , rpA rnA A A rpA A W3 £ - * T ' - ( D \ u I « 1 ) \ i f l 1 2 ^ 1 2 % x « l 2 > " 153. W4 k - I + V 2 + u o¥ 3 ¥5 k 2(W 1 +p 1W 3) As (10.46) i s a function of k only, the minimum of §J occurs when the derivative of (10.46) with respect to k i s zero. This yields the cubic equation k 5 + pk 2 + r = 0 (10.47) where p A ^ o l \ S ^ i ^ o l W ^ o o ) and Using Cardan's cubic formula for the r e a l root of (10.47) yields 1/3 kopt = " 3 ( 2 + ( l + ^ } ( 1 0 ' 4 8 ) I t i s shown i n Appendix H that the rate of convergence for the approach based on (10.48) i s independent of the i n i t i a l scale factor. As a r e s u l t , t h i s approach can be used within the Roc region as a means of providing improved f i n a l convergence. The associated neighbouring trajectory i s given by K = ^ o P t &0 + S V ( 1 0 ° 4 9 ) where \ i s defined on the nominal trajectory, S^-Q i s defined by ( 1 0 . 4 2 ) , k i s defined by (10.48) and u t by ( 1 0 . 4 5 ) . 154. 10.2.2 Computing Algorithm F2 using Second Variation (1) From the nominal trajectory defined by *k and (8.33), integrate (8.30), (8.31) and (9.32) from appropriate i n i t i a l conditions u n t i l S = 0, which defines t ^ . (Note: n systems of (9.32) are i n t e -grated i n p a r a l l e l from the i n i t i a l conditions given i n Section (8.2.2)). (2) Test J"a and/or jg |<e2 ^ o r e x i t , where and £ 2 a r e selected to provide the desired degree of accuracy. (3) Calculate k Q p t from (10.41) and use (10.42) to define S \ . o (4) With \ defined by (10.49) return to ( l ) and repeat. 10.2.3 A.' Curve F i t t i n g Technique to Determine the Optimal Step. Size In the f i r s t v a r i a t i o n procedure of Section 10.2, the A, vector i s swept through the i n i t i a l condition space at a constant angular rotation Set u n t i l J increases. At th i s ° a f i n a l step ((n + l ) - t h ) , two points on the curve J- = J (Set) a a have been established, where Set i s defined as the angular rotation i n the gradient d i r e c t i o n from the n-th step. The two points are (J ,0) and (J -, , S e t ) , where J i s the value of * ao' a l ' , ao J at the n-th step, and J , i s the value of J at the (n + l ) - t h Si c l -L £1 step. As the rotation from the n-th step i s i n the gradient d i r e c t i o n , an optimal rotation Soc ^ must exist which provides the maximum decrease i n J and whichvis-'within the, range a 155. O^Sa o p t<Sa ( (10.50) An approximation to 8* .|_ can be conveniently obtained by approximating the curve J (Soc) by a parabola of the form J = a + bSct + c Set2 (10.51) To solve for the coefficients a, b, and c i n (10.51), only one Figure 10.4 The Graph of Ja(Sct) i n the Neighbourhood of the n-th Step Substituting these three points into (10,51) yields the matrix equation 156, J ao J 1 at J . a l 1 0 0 1 i S a 2 1 Sa S a 2 a (10.52) The value for Sa , i s taken as.that value of Sa which mini-opt mizes (10.51). D i f f e r e n t i a t i n g (10.51) with respect to Sa and equating the derivative to zero yields (10.53) Solving for b and c from (10.52) yields 8 % t = - b / 2 0 g £ (10.54) ao' s ? (10.55) Substituting (10.54) and (10.55) into (10.53) yields the desired value Sa. Sfc <Jal - 4 J a 4 - + 3 J a o ) ^ 4 < TJ - 2J 1 + J T at ao' opt . ,, ,. a l The value for Sa i s replaced by Sa (10.56) opt and the procedure i s repeated u n t i l the desired degree of accuracy i s obtained. 10o2.4 Computing Algorithm F5 Using a Curve P i t t i n g Technique (1) Let J be a large positive number. £10 (2) Use the f i r s t v a r i a t i o n approach of Section 10.2 wi th Sa = Sci. / (3) If l ^ a l < ^ e l a n (^/ o r |-S I <Ce2 ex^-^» where e-^ and e 2 are chosen to provide the desired accuracy. (4) If J-<J n, set J ' a ao' ao J and store the f i n a l a values for this trajectory. ComputeiSx. by (10.34) 157. and return to (2). (5) I f J n, replace Sa by -gSa and, using the f i n a l a ^ ao 1 values of the previous trajectory, compute kj_ from 2 (10.36) and store (Note: a = 1 i n Ra) . (6) Replace \ by (^. + (ki-k)Sx„ i ) and integrate O O 2 O-L (8.30) and (8.31) forward to obtain J i . (7) Compute O a ^ from (10.53). Replace Sa by Sa .|_ and, using the stored f i n a l values, compute k .(. from (10.36). (8) Replace % by (% + (k ,-ki )S^„n) and return O O OpX 2 O-L to (2). 10.3 The Combined Computing Algorithm In t h i s section a combined algorithm i s presented which uses the f i r s t v a r i a t i o n approach to i n i t i a t e the search procedure and to locate the Ra region. Once i n the Ra region, one of the computing algorithms PI, F2, or F3 i s used to pro-vide the f i n a l convergence. A flow graph of the combined algorithm i s shown i n Figure 10.5. 10.3.1 Example 1 As a f i r s t example, the problem of Section 9.5 i s solved using the combined algorithm with F l and F2 as the f i n a l stage. The resul t s are shown i n Table 9.2. I t i s seen that the combined algorithms provide the desired improvement i n f i n a l convergence over the algorithm based on f i r s t v a r i a t i o n alone. A further comparison between F l and F2 i s i l l u s t r a t e d i n Table 10.1 and Table 10.2 for which the search procedure i s i n i t i a t e d from different points i n the i n i t i a l conditions space. It can be observed that, for t h i s example, the computing 158. START \ - R E i D / PERFORM INTEGRATIONS TO OBTAIN B0 AND Zg I CALCULATE b EQN. <70-3i) Ao =Ao + W A< J 1 a =; I COMPUTE K EON (lo-ib) NO Jao = Ja 6 Xo =ef Xoo+kdtol I I C/SE F / , F2,0/? F J 7^ 0 CALCULATE AO = no + </ A~o NO STOP Figure 10.5 Flow Chart for the Combined Algorithm 159. Table 10.1 STEP P2 - Second Variation on k PI - Matching End Points X 2 f .. X l f . X 2 f aS 0 1 2 0.166528 0.56796 0.659119 -197.1 -4.459 0.4128 8.869 1.627 0.8804 0.166528 0.56796 0.659119 -197.1 -4.459 0.4128 8.869 1.627 0.8804 B 3 4 ' 5 6 0.676695:: 0.681707 0.681707 -0.3273-.-~-0.00262 -0.0000043 -1.-006 0.8379 1.001 0.6 77784-0.681636 0.681707 0.681707 -0.2944 -0.03482 -0.0006664 0.0000003866 1.039 1.004 1.000 1.000C Table 10.2 P2 - Second Variation on k E l - Matching End points STEP X 2 f X U X2f X 2 f . ... X l f X2f . 0 0.006379 -85020. 1113. 0.006379 -.85020. 1113. 1 0.011924 -36610.' 529. 0.011924 -36610. 529. as 1 2 0,029196 -9714. 175.3 0.029196 -9714. 175.3 PA 3 0.094367 -134.8 50.61 0.094367 -1348. 40.61 4 0.339618 -76.88 • 7.783 0.339618 -76.88 7.783 5 0.670167 0.7321 2.036 0.670167 0.7321 2.036 6 0,680692 -0.1371 1.007 0.680782 -0.1318 1.017 a 7 0,681707 -0.002836 0.9509 0.681703 -0.008396 1.001 PA 8 0.681707 -0.00000391 0.9993 0.681707 -0.00003993 1.000 9 0.681707 0.00000057 1.0000 1 160. algorithm F2, which i s based on a second v a r i a t i o n approach, appears to provide the best f i n a l convergence. For both F l and F2, however, the convergence i s esse n t i a l l y quadratic. 10.3.2 Example 2 Consider the problem for which the control u(t) i s to be found over the time i n t e r v a l 0 ^ t ^ 4 . 2 , such that the system performance J = (2x 2 + x~) (10.57) x 0 t=4.2 i s a minimum. The equations of constraint are x-j^ = ( l - x 2 ) x 1 - x 2 + u , x 1(0) = 0 X n = X x o(0) = 1.0 (10.58) 2 - A i • » A2 x 3 = x l 2 + x 2 2 + u 2 > x^(0) = 0 and g = x 2 ( 4 . 2 ) = 0 (10.59) The associated Euler-Lagrange equations are \ ± = - ^ ( l - X g 2 ) - \ 2 - 2 \ ^ L ± (10.60) X 2 = A.1( 1+2x^2) - A,2 - 2\^x 1 \, = 0 3 H. = -X, - 2\,u = 0 u J. .5 For the c l a s s i c a l approach, the f i n a l transversality conditions are X 1 ( 4 . 2 ) = - 4 x 1 ( 4 . 2 ) (10.61) \ 2(4.2) = -V \ 3 ( 4 . 2 ) = -1.0 161. where V i s defined by the augmented performance function J = (2x-.2 + x 0i>+ x,) (10.62) a 1 2 3 t = 4 > 2 This problem i s an example of a problem with f i n a l end con-s t r a i n t s . The solution i s obtained using the combined algorithm with PI, F2, and F3 and the f i n a l stage. The optimal t r a j e c -tory and optimal control u(t) are shown i n Pigure 10.6. The associated Lagrange m u l t i p l i e r s are i l l u s t r a t e d i n Pigure 10.7. The values X-^O) = 0.5, A-2(0) = 5.0, and \^(0) = 1.0 were a r b i t r a r i l y selected as an i n i t i a l estimate. I t can be ob-served i n Table 10.3 that t h i s i n i t i a l estimate i s i n the region Rg since b = 0.3*0• An incremental rotation Sex = 0.005 radians i s used during the f i r s t v a r i a t i o n approach. Table 10.3 i l l u s t r a t e s very c l e a r l y the region Rg (steps 0 to 7) i n which a = b<|l and k = 0, the region RJ (steps 8 to 15) a i n which a = 1 and k / 0, and the region Ra ( f i n a l steps) i n which Sa !js Sa, • Notice that i n region Rg, f u l l emphasis i s placed on reducing g, and i n region RJ the emphasis i s trans- . c l ferred to reducing J while maintaining g close to zero. The intersection of the Ra region i s manifested by the increase i n •J at the 16-th step. Once i n Ra, PI, P2, and P3 are used to provide the f i n a l convergence. It can be observed that, for th i s example, F l and F3 provide more rapid convergence than P2. The convexity of J_ i n the neighbourhood of the 15-th step i s c l shown i n Figure 10.8. Note that the method of curve f i t t i n g , F3, estimates quite accurately the minimum of J , and that the c l corresponding estimate by the second v a r i a t i o n approach, F2, i s s l i g h t l y i n error. It i s believed that t h i s error i s caused Figure 10.7 The Optimal Lagrange M u l t i p l i e r s for Example 2 Table 10.3 STEP J J a g b " k V : o- -'• 41.353554 « -0.879200 0.302 0 — 1 27.258678 — -0.687343 0.207 0 — — 2 19.727938 — -0.588630 0.186 0 — — 3 20.036128 — -0.86092 0.223 0 — — 4 17.422018 — -0.382174 0.281 0 — — 5 17.310112 — -0.278221 0.382 0 — . — 6 14.493840 — -0.174874 0.601 0 — — 7 14.140386 :' 17-540385 -0.072850 1.421 0.000062 -46.67 2244.0 8 8.642633 - 10.201530 -0.071281 2.700 0.000195 -21.87 512.0 9 7.006617 7.058201 -0.004519 50.0 0.000290 -11.42 154.0 10 5.183621 : "5.305702 -0.021428 14.1 0.000514 -5.697 46.7 11 4.193642 4.219338 -0.007691 47.0 0.000770 -3.341 19.64 RJ 12 3.508284 - 3.'521721 -0.006890 58.8 0.001218 -1.950 8.33 a 13 3.096564'' 3.101889 -0.00480 90.7 0.002195 -1.109 3.403 14 2.902321 2.904174 -0.00339 134.0 0.006832 -0.5452 1.402 15 2.899740 . 2.900013 -0.002182 214.0 0.007452 -0.1250 1.189 16 2.903164 • 2.904198 -0.001860 — — — — E2 — Second -Variation on k P3 - Curve P i t t i n g PI - Matching End Points SIEI a STEP J g a g R a 17 18 19 20 21 22 23 24 25 2.884124 2.880123 2.879259 2.879046 2.878997 2.878985 2.878982 2.878981 2.878981 -0.000907 -0.000257 -0.000055 -0.000014 -0.000003 -0.000001 -0.000001 -0.000000 -0.000000 0.005334 0.004583 0.005135 0.005179 0.005160 0.005173 0.005168 0.005170 17 18 19 20 21 2.879017 2.878982 2.899972 2.884135 "2.878981 2.878982 2.878981 -0.000420 -0.000387' -0.000506 -0.000122_ 0.000000 0.000001 .0.000000 0.003551 0.003986 2.879007 2.878981 2.878981 ;o. 000371 -0.00002 0-000000 0.001054 0.000001 0.000000 166. by the uncertainty i n the value for V on the nominal trajectory of the F2 approach. By (10.43), i t i s seen that ^ i s a function of the future step s i z e , and hence J & on the nominal trajectory can only be computed after the future step has been specified. This type of problem i s not present with PI and P3 since for F l , V i s precisely specified on the nominal trajectory as that value which minimizes the error i n transversality (10 . 3 2 ) , and for P3, the actual curve J (6a) i s used. ' a 10.4 The Second Variation Technique Used with the Method of Steepest Descent A second v a r i a t i o n approach can also be used to determine the optimal parameter k for the method of steepest descent i n function space. Consider the optimal control problem i n Section (8.4)• Using (8.23), the v a r i a t i o n on J up to terms a of second order i s T S j a = S x f T 0 x f + & f \ x i S x f - i f ( Su TH u u Su)dt 0 ^ ( 2 S u T H u x S x + S x T H x x S x ) d t -i J ~ ( S\ TH x ? S\)dt . 0 . 0 T T -f (8u TH u+Sx T(H x+\) + S\ T(H x - x)dt - i J ~ (S\T-' 0 0 •(Sx-B\ 8x-H\ Su)dt (10.63) A.X A . U For the method of steepest descent developed i n Section (8.5), the following relations hold for the nominal trajectory: 167. (1) x - = 0 (2) \ + H x = 0 (3) H u , 0 (4) x(0) = x„ (5) \ f + 0 x f = 0 from ( 8 . 3 0 ) from ( 8 . 3 1 ) from ( 8 . 2 4 ) from ( 8 . 3 3 ) from ( 8 . 3 4 ) (6) Sx - H. Sx - H. Su = 0 , from (9.25) A.X A . U Using ( l ) to (6) i n (10.63) provides T S j a = - ^ f S u \ d t + i S x f T 0 x x ^ x f 0 T - i J ~ (Su TH u uSu + 2Su TH u xSx +Sx TH x xSx)dt 0 (10.64) Prom (8.38) the desired value for Su i s Su = kH u (10.65) The incremental v a r i a t i o n Sx(t) which results from this incre-mental change Su, i s given by (9.5) to be S . A p> A O -x = H X x6x + H X uOu Substituting (10.65).. into (10.66) yields S « A C" A A x = H. Ox + LL kH A.x \u u Let Sx =r kZ and hence, using (10.68) i n (10.67) results i n (10.66) (10.67) (10.68) • A A A Z = LL Z + H \x \u u (10.69) 168. which i s independent of the parameter k. The solution to (10.69) for Sx Q = Z(0) = 0 i s given i n Section (8.2.2) i n terms of the state t r a n s i t i o n matrix 3>: z(t) = C S(t , a ) H X u ( a ) H u ( a ) d a ( 1 0 . 7 0 ) Using (10.65) and (10.68) i n (10.64) provides T s J a = . k f \ \ ^ + 4- z f \ x t zi 0 T \r2 (~* , A rpA A A rnA TA N - f - / H XH H + 2H XH Z + Z \ Z dt 2 / u uu u u ux xx •0 (10.71) For §J to be a minimum with respect to k, §J must vanish and hence T § 2J = - / H TH + k ¥ = 0 (10.72) ^ a J u u o 0 where A rpA / , A rpA A A rnA rpA ¥ = z/0 -Z- - / (H H H +2H LE Z + ZXH Z)dt o f F x x f f J K U U U U U U X X X 0 Solving (10.72) for .yields T A rpA dt kopt = °" ¥ ( l 0 ' 7 5 ) r o and hence, from (10.65) and (10.73), 169. T H XE dt) u u Su(t) 0: A , H t) u ' (10.74) ¥ o i s the desired v a r i a t i o n for 6u i n the v i c i n i t y of the extremum.. Using a preselected value for Slz i n (8.39), the value of k for the f i r s t v a r i a t i o n approach i s The method of steepest descent i s carried out as before with k defined by (10.75) u n t i l the s i t u a t i o n exists where k ^ ^ ^ k . When t h i s occurs, the second v a r i a t i o n approach of (10.74) i s used for the remainder of the search to determine the va r i a t i o n Su(t). In a manner si m i l a r to that of the previous sections, this approach can be extended to cover problems with free f i n a l time and additional terminal constraints. 10.5 Conclusions An algorithm for the numerical solution of optimal control, problems has been developed which i s based on a combination of the direct and in d i r e c t approaches. The method i s s i m i l a r to the indi r e c t method i n that tra j e c t o r i e s are com-puted using d i f f e r e n t i a l equations and the known and computed i n i t i a l values. However, instead of matching end conditions as i s done i n the c l a s s i c a l i n d i r e c t approach, the augmented performance function J i s considered to be a function of the a unknown i n i t i a l values. The minimum of J i s found by gradient 9.-search i n the i n i t i a l condition space based on the f i r s t v a r i a t i o n . Unlike the in d i r e c t approach, convergence does not (10.75) 170. depend on a good i n i t i a l estimate of A.Q. It has been shown that the normalization of the Lagrange m u l t i p l i e r s \ , as carried out i n the c l a s s i c a l approach, i s not essential. Thus, instead of a search over the complete A, Q-space, i t i s s u f f i c i e n t to determine the intersection of the l i n e \x\ with the sphere A TA \ XQ = constant. Also, i t has been shown by means of examples that the method i s applicable to the case of bounded control, and that i t can be applied without computational d i f f i c u l t y to the case where u cannot be determined as an e x p l i c i t function of x and \U In the case of bounded control, some prior knowledge of the sequence of arcs i s required. Information of t h i s type can be determined with the aid of the Legendre-Clebsch condition. A disadvantage of the gradient technique, based on the f i r s t v a r i a t i o n only, i s that the convergence slows down as the optimum i s approached. It i s desirable, therefore, to use the gradient technique to i n i t i a t e the search, and then to use a technique with good f i n a l convergence properties to complete the search. In t h i s respect, a three stage computing algorithm was developed which i s based on a systematic search i n the i n i t i a l condition space of Lagrange m u l t i p l i e r s . The f i r s t two stages are steepest descent techniques which re s u l t i n a search over the surface of a sphere, at a constant sweep angle, u n t i l the region of the optimum i s overstepped. At t h i s point, the t h i r d stage comes into effect to provide a rapid f i n a l convergence. For t h i s t h i r d stage, three algorithms were developed. The f i r s t approach i s based on a method of matching end points i n which the Lagrange m u l t i p l i e r s are continuously re-scaled to 171. provide a minimum error i n f i n a l transversality. The other two approaches are based on finding the optimal step size for the method of steepest descent used i n the second stage. One approach uses a second v a r i a t i o n of the augmented performance function and the other uses a curve f i t t i n g technique to estimate the optimal step s i z e . I t i s e a s i l y observed that these techniques vary i n computational d i f f i c u l t y and i n the rate of convergence. For the method of matching end points, some matrix inversion i s required which may introduce computa-t i o n a l d i f f i c u l t i e s . Por t h i s technique, however, i t was demonstrated by examples that the rate of convergence i s essen-t i a l l y quadratic. The method of curve f i t t i n g , on the other hand, requires no extra equipment when used with the steepest descent approach. However, one extra forward.integration i s required f or each i t e r a t i o n . For this approach the rate of convergence was shown to be very satisfactory. In the second v a r i a t i o n approach, no matrix inversion and no extra forward integrations are required, and the rate of convergence appears to range between l i n e a r and quadratic. The pa r t i c u l a r approach to use, therefore, w i l l depend on the problem under study and the size and type of computing f a c i l i t i e s available. A common feature of a l l these techniques isvthat storage i s required at the, end points only and, as such, these, techniques are suitable for use with d i g i t a l or hybrid computers of lim i t e d memory. 172, Appendix A The equations of motion for a missile moving i n the earth's atmosphere are (see Figure 2.1, for s i m p l i c i t y , a f l a t earth and motion i n the xy-plane i s assumed). x = v cos 9 e y = v s i n 9 D v e u C 0 S ^ v = - g B i n 0 - ; + ~r-' _ g cos 9 + L + v e u s i n P v mv mv m -u (A-l) (A-2) (A-3) (A-4) (A-5) Here D = D(y, v, L) i s the drag and i t i s assumed that the engine thrust i s given by the i d e a l equation T = v u; Equations (A-l) to (A-5) can be written i n the form (2.1) by choosing x-^ = x, x 2 = y, x^ variables and by taking v cos 9 v s i n 9 Gf, A 0 -g s i n 9 - -g cos 9 L_ v mv 0 v, x. = 9, x,- = m as the state 4 5 0 0 v cos 8 e m v s i n 8 e m -1 (A-6) — 1 L —' It follows from (A-6) that 0 0 cos 9 -v s i n 9 0 0 0 s i n 9 V cos 9 0 Ox 0! D m D V m • -g COS 9 D/m2 0 0 (mg cos 9-L)/mv g s i n 9/v -L/m 0 0 0 0 0 (A-7) 173. Gr. l x G OL 0 0 0 0 0 0 0 0 0 0 0 0 -DL/m l/mv 0 0 0 0 -v s i n S/mV e 0 0 0 0 0 2 0 -v cos fi/m e 0 -v s i n p/m v 0 0 (A-8) 10 0 0 -v s i n 6/m e v cos 6/mv e r 0 (A-9) The Euler-Lagrange equations are X± = 0 ^2 - m (A-10) ( A - l l ) Xr \ , = COS 9 - \ 0 s i n © + ~ D, - A.„ ^ cos 9 3 1 2 m v 4 y2 X + — ^ (L + v u s i n p) mv (A-12) X^ = A^v s i n 9 - \ 2 v c o s e + ^ S c o s e ~ ^4 f s i n e (A-13) S = " 2 ( D " v e u ^ m Equation (A-10) yields X± = C-L L cos p) + —2^ — (L + v gu s i n p) m v (A-14) (A-15) where i s a constant of integration. Substituting (A-6) and (A-9) into (2.12), (2.13) and (2.14) yields 174. k T = i (\, D - ^ ) (A-16) L m 3 L v k8 = ST U 3 s i n S - ^ cos 8) (A-17) k u = 5s" U 3 cos p + ^ s i n p - \ 5) (A-18) Substituting (A-6) into (2.18) yields \,v cos 9 + \ 0v s i n 9 :- \- (g s i n 9 + —) i c J ni + \ A ( r - - & cos 9) + u k = c (A-19) 4\mv y u Evaluating (2.30) with the aid of (A-6), (A-7) and (A-8) yields (sin 9 s i n 8 - cos 9 cos 8 ) ^ - (cos 8 s i n 9 + s i n 8» cos Q)\ 2 + <^(Dv cos 8 + ~-' )jjj. + f s i n 8 cos 9>\, e k u = m +<^-^- cos 9 cos 8 + •^ ~p s i n 8 + ""^ p ( c o s 8 _ T - ) ^ x v mv mv e' (A-20) The transversality condition i s t t dP + A^dx + A-2dy + ^3 < i v + ^4^® +A.^ dm - cdt f = 0 (A-21) 0 During a variable thrust subarc (2.45) represents the following system of equations k L = 0 x l = c l G-0T\ = c (A-22) k = 0 u k = 0 u I t follows that 0 A- = 0 1 v cos 9 0 l_ a 5 1 0 0 0 v s i n 9 0 a 52 and b = 0 0 c-c 0 0 where a A m m s i n 3 0 -g s i n 9--v — cos 8 m • a 53 -v mv _ 1_ mv cos 8 0 L_ _ £ mv v v — s i n 8 mv K a 54 51 a 52 s i n 9 s i n 8 - cos 9 cos 8 -(cos 8 s i n 9 + s i n 8 cos 9) a 5 3 = m ^ Dv c o s 6 + + f s i n 8 cos 9 v 175, 0 0 0 cos 9 0 -1 0 (A-23) (A-24) (A-25) a 54 = - ^ cos 9 cos 8 + — 2 s i n P + — 2 ' c o s P ~ ~ v mv mv e Appendix B In the case of v e r t i c a l f l i g h t , the missile i s con-176. strained so that 9 = (3 = 0, L = 0 and the system dynamics (A-l) to (A-5) simplify to the form (B-l) y = v D . e V = - g - — + s m m m = -u Using a standard drag function of the form D = K V a y a (B-2) (B-3) (B-4) where K and a are constants, results i n the following Euler-• a & Lagrange equations a\,D m 2\_D X3 = " X2 + mv X5 = _ 2 ( D" v e u ) m (B-5) (B-6) (B-7) The switching function and i t s time derivative are given by (B-8) \.~v k__ = -4-JL _ v u m \ 0v^ \, 2v m (B-9) The f i r s t i n t e g r a l i s \ 2v - \ 5 ( g + -) + uk u = c The transversality condition i s (B-10) dP + X^dy + \^dv + X,^ dm - cdt = 0 ( B - l l ) J0 The matrix A and the vector b are given by (see (2.45) and (2.46)) 177. A = b = V - g -V 0 e m He D m 2 m c 0 0 D m 2v 0 -1 0 (B-12) (B-13) In Chapter 3,/use was made of the fact that X^ = 0 i n deriving the sequence diagram. To prove t h i s assertion,note that during a variable thrust subarc (k^ = 0) or during a coasting subarc (u = 0), (B-10) reduces to X2 V ~ X 3 ( g + m} 0 (B-14) after substituting the condition (3.2). It follows from (3.1) and (B-14) that at the f i n a l time 3f 0 (B-15) If \ j = 0 at any instant during the variable thrust or coasting subarc i t follows from (B-14) that ,\2• •= 0, and from (B-5) and (B-6) i t can then be concluded that X^ = 0, X^ = 0 everywhere, v i o l a t i n g the terminal condition for A,^ (see 3*2). Consider now a maximum thrust subarc. Substituting c = 0 and (B-10) into (B-6) yields X, \_ = uk + (D - mg) u mv (B-16) If \_ = 0 at any instant t during a maximum thrust subarc, i t follows from k u> 0 and (B-16) that \^>0 and consequently \,>0 for t < t = t„. From (B-5) i t i s then seen that 3 s ^ f 178. \ 2<0 (B-17) for t <t^t.p. Evaluating (B-10) at t = t yields S I s uk X2 = - < 0 (B-18) Conditions (B-17) and (B-18) require that A_2< 0 for t g < t ^ t f , which vi o l a t e s the terminal condition \ ^ - !• Thus X^ £ 0 for tQ<t<Ct^. To prove that \^>0 i t i s s u f f i c i e n t to note that i f (B-4) i s used to evaluate (B-6) at t = t ^ , i t yields \ 5 f = -1 (B-19) and thus the f i n a l value of X^ i s approached through positive values, since X^ = 0 (see B-15). Prom (B-2) i t follows that v>0 during an impulsive thrust subarc. The f i n a l subarc, therefore, cannot be an im-pulsive thrust subarc. Prom (3.1), (3.2) and (3.6) i t follows that v k (t J = - -£<0 (B-20) u f m^ which violates the condition (2.42) for a variable thrust sub-arc. Thus the f i n a l arc i s a coasting subarc. Appendix C The system dynamics for f l i g h t with the control con-st r a i n t s L = 0, 8 0 are x = v s i n 9 (C-l) y = v cos 9 (C-2) 9 = - £ cos 9 m = -u The Euler-Lagrange equations are \ 1 = 0 aA,D X2 = " m 179. (C-4) (C-5) (C-6) (C-7) A,2D \ 5 = -\1 cos 9 - \ 2 s i n 9 + X^2 cos 9 (C-8) A^ =• X^v s i n 9 - A_2 v cos 9 + X^g cos 9 - X^ s i n 9 (C-9) (C-10) X. X5 = " "2 ( D" v e u ) m The switching function and i t s time derivatives are ( 0 - 1 1 ) k = -±-2- - X u m v k = u m - X-^ cos 9 - \ 2 s i n © + -3=-(2-£- ) mv v L ^ cos 9 v The f i r s t i n t e g r a l i s (C-12) D 4. X-jV cos 9 + \ 2v s i n 9 - X- (g s i n 9 + -) - ^ g cos 9 + uk u dP The transversality condition i s + X-^ dx + \ 2dy + X^dv + A^d9 + A^ dm - cdt t (C -13) 0 (C-14) The matrix A and the vector b are given by (see (2.45) and (2.46)) 180. A = 1 v cos 9 0 - cos 9 0 v s i n 9 0 - s i n 9 '1 0 -g s i n 9— v m £-(2 + 2-) mv v ' e 0 cos 9 0 -^2 cos 9 v 0 0 -1 0 (C-15) (C-16) 0 0 where A-^ = c^ = constant follows from (C-6). A proof w i l l now be given of the fact that A^ £ 0, on a variable thrust subarc. During variable thrust (4.6) i s v a l i d . Substituting (4.6) into (C-10) yields Ar \ 5 = - ( D - V J A ) mv (C-17) If A^ becomes zero at any instant t ^ on a variable thrust sub-arc, i t follows from (4-6) and (C-17) that A^ and then remain zero for the remaining time i n t e r v a l t ^ =• t ^ t 2 °^ the variable thrust subarc. From (C-8) and (C-13) i t follows that since A.., and k are zero 3 u A-^ v cos 9 + \ 2 V S x n ® ~ ^4 Y c o s ® ~ 0 \-L cos 9 + A 2 S ± N ® + A4 ^2 c o s ® = ^ for t ^ = t =" t,,. Hence A^ = 0 and A^ cos 9 + A 2 s i n 9 = 0 (C-18) (C-19) From (C-6) and (C-7) i t i s seen that A-^ = const., A 2 = const, 181. Thi is implies that 9 = const, for t ^ ^ t ^ t 2 which contradicts (C-4). Hence \^ / 0 I 0 r " t n e variable thrust subarc. To determine an expression for k , (4.4) i s f i r s t written i n the form v k = — f (C-20) u mv ! where f = u K - " *4 c o s 9 ( 0 " 2 1 ) Substituting (C-13) with c = 0 into (0-8) and using (C-20) yields \_ = ^ -..k - -2- (C-22) v3 v u mv e e Substituting (C-13) with c = 0 into (0-9) yields D i f f e r e n t i a t i n g (4.5) with respect to time and using the system equations (C-l) to (C-5) yields f = uM, + Nn (0-24) s 1 1 where Mx k -g s i n 9 - ^ (2 + |3L) (C-26) e N A _ mg2 cos 2 9 D ( D ) ( 2 3v } ( c _ 2 6 ) 1 v v to m v e Di f f e r e n t i a t i n g (C-21) with respect to time yields °P ' i i' % f s , f s * 2 g cos 9 . 2g . Q A f = uk u + uk u - \ 5 — - \ 3 — - \ 4 — ^ + \ 4 ^ s i n 9 9 + X A 0 0 3 9 (C-27) ^ v 182. The term uk^ i n (C-27) i s always zero since u = 0 when u i s at i t s bounds of k = 0 when u i s variable. Eliminating the time derivatives i n (C-27) with the aid of (0-3), (C-4), (0-22), (0-23) and (C-24) yields f = u 2gv . 2k g 2 „ k k, + \,M„ + A. j ,—^ cos © - - r * - :?g n " - ~ f u -3-2 -4 m v2 v s i n 9' v s e + + XJNJ - \±2g cot 9 (C-28) where m e (0-29) f D N 2 =~T~ + m v A ^ s ^ . 2Dg cos 9 g C 0 S 9 _ ( G S I N Q + P.). mv s i n 9 v mv 63 nr (2 + 22.) _ aDv s i n 9 / , V _ N v v m v v A 2g cos 9 j ~ g cos^ 9 - s i n ^ 9 _ D si n 9 m Evaluating k u when k u equals zero yields (0-30) (0-31) v . k = — f u mv (C-32) Substituting (C-28) into (Cr-32), using (4.5) and taking k u = 0 yields v u mv 2, u .k v 1 u e mv uk ( S—E + —) + \,uM - \,N - \ n2g cot 9 u vv sm 9 mv 3 3 1 (C-33) where M A D ( ^ e v_ } _ Z£ s i n Q 2 w v v mv m e (C-34) N A Bg s i n 9 mv 2 + |v _ (3 +S_.) C o t 2 9 + V e e m v 183: 2 .(•* . 5 l . lL_) _ g 2 s i n 2 9 aDv s i n 9 ( 1 +v_\ ^ 5 + v e + y 1] v + m U + v e J , e (C-35) Appendix D The system dynamics for f l i g h t with the control con-s t r a i n t L = 0 are x = v s i n 9 (D-l « y = v cos 9 (D-2 D v u v = - g s i n 9 - - + cos 8 (D-3 ° . m m v u 9.= - £ cos 9 + - ~ — s i n 8 (D-4 v mv m = -u The Euler-Lagrange equations are (D-5 \ 1 - 0 (D-6 i2 = -X 3 f (L-7 •p. X. v u -\1 cos 9 - X2 s i n 9 + 2 \ ^ - | (g cos 9 |- s i n p) (D-8 X^ ~ XjY s i n 9 - X.2V c o s e + c o s e ~ ^4 y s i n Q A.., \ 5 = - -| (D - v gu cos p) + \ 4 -|- s i n 8 (D-10) V U —2 _ 8) /^ ~2~" m m u The equations for the switching function, i t s time derivative and k^ are ku = IT ( X 3 cos 3 + ^ sin 8 - ^ - ^ ) 184. (D-ll) k = -a u m (sin 9 s i n 8-cos 9 cos 3)A.^ - (cos 3 s i n 9 -h'sin 8 cos 9 ) \ 2 +<^~ (2 cos 3 + J-) + f sinp' v. V • cos 9 \ \ , + <- cos 9 cos 3 + — ^ s i n 3NX, / 5 \ v 2 mv2 7 1 (D-12) = -s (\ s i n 3 - - r 1 cos 3) 3 m 3 v The f i r s t i n t e g r a l i s (D-13) A^v s i n 9 + \ 2 V c o s ® ~ ^3 (fi s i n ® + JJj) ~ ^4 y c o s 9 + u k u = c (D-14) and the transversality condition i s - i t -dP + X^dx + X 2dy + X-^ dv + X^d9 + X^ dm - cdt = 0 (D-15) J0 The matrix A and the vector b are given by (see (2.45) and (2.46)) A -0 1 .0 0 ve ' — s i n 3 m 0 — cos p mv ^ 0 D g v cos 9 v s i n 9 -g s i n 9 - — - & cos 9 ° m v 0 0 v — cos 3 m ^ — s i n 3 mv M a 51 a 52 a 53 a 54 0 0 0 -1 0 (D-16) where a51 = s ^ n 9 s x n P ~ c o s 9 c o s P 185, a^ 2 = - (cos 0 s i n © + s i n 0 cos) a53 = mv ^ 2 c o s P + v~) + v~ s i n 13 c o s 9 2_, . £_ (D-17) a 54 A s : D = - ^ cos 9 cos 0 + — 2 s : L n 0 v mv and b = 0 °] c 0 0 (D-18) where \^ = c^ constant. For z e r o - l i f t f l i g h t , (A-6) to (A-9) take the forms °o -'Ox Gr. l x v cos 9 0 v s i n 9 0 -g sin 9 - D m » s i -V, i f c o s 0 (B-19) - cos 9 0 0 -1 0 cos 9 -v s i n 9 0 0 0 si n 9 -v cos 9 0 0 aD m -2D mv' - g cos 9 D 2 m (D-20) 0 0 ^2 cos:9 V ^ s i n 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ve 2 cos 0 (D-21) 0 0 _ ve . — 2 sm mv 0 0 \ 2 m si n 0 0 0 0 0 0 186. G0L = 0 (D-22) G. 10x G, op 0 0 — s i n p o (D-23) Di f f e r e n t i a t i n g (D-21) with respect to p yields 0 0 0 0 0 0 0 0 0 0 0 0 0 v. mv '— cos p 0 0 0 0 0 0 0 0 v mv e mv '— s i n p 2" cos P 0 (D-24) Di f f e r e n t i a t i n g (D-23) with respect to p yields 0 0 G opp v cos p m r v m s i n p (D-25) 0 For the case where u =co> i t i s possible to max -derive a control law for p for the maximum (impulsive) thrust subarc. Equations (D-3), (D-4), (D-8), (D-9) and (D-10) for 187. u—•oo, take the form v m v £? — cos 8 (D-26) m M v m 9 * - -2- s i n 8 (D-27) mv ^ \,v m \, = - * | s i n 8 (D-28) mv ! \_v u \ * c o s 0 (D-29) 5 m The remaining state variables and Lagrange mu l t i p l i e r s do not change during the i n f i n i t e s i m a l l y small i n t e r v a l of time i n which the impulsive thrust occurs. Thus \ 4 = const. (D-30) during the impulsive thrust subarc. I t follows from (2.9) and (2.13) that k^ = 0 and hence tan 8 = ^ (D-31) 188. D i f f e r e n t i a t i n g (D-31) with respect to time and noting (D-28) and (D-30) yields 3 v m s i n 8 r\j • e ~ mv Adding (D-27) and (D-32) yields 3 + 9 ^ 0 Integrating (D-33) yields 8 + 9 ^ 8 0 + e o = const. Substituting (D-32) into (D-26) yields v ^ - v 8 cos 3 si n 8 Integrating (D-35) yields c 2 = v Q s i n 8Q - v s i n 8 where c 2 i s an integration constant. Hence 8 ~ arc s i n v — s i n 8 A v M0 0 (D-32) (D-33) (D-34) (D-35) (D-36) (D-37) Equation (D-37) gives the control law for 8 during impulsive boosting. Substituting (D-36) into (D-32) yields 2 v = - u - -3-)* (D-38) m c 0 1 (1 - - § - ) T m v 2 Separating the variables i n (D-38) and integrating yields / 2 2N* , 2' 2 J ( c 2 - v ) - ( c 2 - v Q ) -1 m -v In :— e m0 (D-39) Substituting (D-36) into (D-39) and solving for m yields m = m^ exp v n COS 8Q ~ v cos 8 v_ (D-40) Equations (D-37) and (D-40) hold during the impulsive boosting 189. subarc. Appendix E After the n + 1-st step, in..the .=Qn.e:;dimensi.Qnal. search procedure, assume a condition exists of the form P ( a ^ " 1 ) < P ( a k n ) >. P ( a k n + 1 ) (E-l) and, hence, the optimal value of ct^ . i s located i n the region (see Figure 4°10) (E-2) n-1/ n+1 kopt ^ k To determine a^ . ^ approximately, the region i n the v i c i n i t y of the optimum i s represented by a parabola of the form 2 P = a + bct k + ca^ (E-3) Using the three coordinates (^n_]_, \ ^ > (^ n> a\ 1^ and (P. n+1 , a vn + 1 ) i n (E-3) yields 'n+1 • n 'n-1 1 1 1 a. a. h+1 n K n + 1 > a. n-1 ( a k } a (E-4) from which the values of a, b, and c can be established. From ordinary calculus (^JJ- = 0), the extremum of (E-3) i s located at a kopt k b_ 2c Substituting for b and c from (E-4), and noting that n-1 n . a, = a, Act, k k k n+1 a n k + Aa r (B-5) yields 190, Act, (P _ - P -, ) , k n-1 n+1 / - c i c \ «*opt = « , . + — ( P n + 1 - P n ) + ( P ^ - P n) (E-«) Appendix P The optimal value of S\. i s given i n (10.34) to be S^o = a S\>o + k S \ > l ( F " 1 } where ^ o l - Z g < Zg T Z / Z 0 - 2 0 ( p - 3 ) Prom (9-43), the va r i a t i o n Sg for the neighbouring trajectory i s Sg = Z g [ IS\ o (P-4) I t i s desired i n th i s appendix to determine the effect of SX q o and S\ o l on Sg. Using (P-2) i n (P-4) yields Sg = Z g T ( - Z g ( Z g T Z g ) _ 1 g ) = -g (P-5) Hence, by (P-5) i t can be concluded that S\ Q 0 i s the component of S\ concerned with s a t i s f y i n g the desired end conditions. Using (P-3) i n (P-4) yields Sg = Z T(Z (Z T Z ) 1 Z T Zfi - ZD) = 0 (F-6) g g g g g 0 0 Therefore, by (P-6), i t can be concluded that S\ ^ attempts to minimize the system performance without affecting the end conditions set by SA.Q0» 191. Appendix Gf To- f i n d i?, substitute (9.39) into (9.43). This yields a g = k zg L(z0 + ZJS) g (G-l) Hence 1P= (Z T Z )~1(§ g - Z T Zw) g g k g 0 To fi n d k, substitute (9.39) into (9.44). Noting that T (0-2) (Z 1 Z ) g g -1 = (Z T Z ) , and using (G—l) and (G--2) yields g g -1 k = £ n2 2AT /_ T „ N _ 1A O l - a g , ( Z g Z g)/ g Z,/ Zw - ZwT z (z T z ) 0 0 0 g g g -1 T g 0 (G-3) The res u l t (9.45) i s obtained by substituting (Gf-3) and (Gf-2) into (9.39). Appendix H -For the f i r s t v a r i a t i o n approach, the, d i f f e r e n t i a l equations used are x = f(x,u) , from (8.30) (H-l) \.= f„\ , from (8.31) (H-2) (H-3) The solution of (H-3) i s given by (9.3) and (9.4) to be Sx(t) = S 1 2(t,0) S\„ (H-4) Sx" Sx S\(.t) = ffi22(t,o) S\ ( (H-5) 192. It i s desired i n this appendix to determine the effect of multiplying a l l Lagrange m u l t i p l i e r s by a non-zero scale factor c. Let the primed quantities represent the values obtained when \ = c\Q, and l e t the unprimed quantities represent the values obtained when X = X . Prom (H-l) and (H-2) i t i s seen that o o x'(t) = x(t) (H-6) and x'(t) = c\(t) (H-7) rn Prom (H-6) and (H-7), and the d e f i n i t i o n H = \ f, i t i s seen that H' = cH •H1 = cH u u < = \ H x = cH x (H-8) \ H = cH uu uu Hux~' - - c H u x H' = cH XX XX (H" 1 ) ' = — H-"*" uu c uu Prom equations (9.2), (9.3), (9.4-), and the relations (H-8), i t can be determined that d^ 2(t,0) = J SB 1 2(t,0) S 2 2(t,0) = ffi22(t,0) (H-9) Sx'(t) =Sx(t) and S\' (t) = c S\(t) Using (H-9) i n (9.39) yields 193. 2 0 = Z 0 (H-10) and z' = - Z g c g From the d e f i n i t i o n of S l 2 i n (9.55), and (H -10), i t i s found from (G-3) and (G--2) that (Si 2)' = o 2 S i 2 k' = c 2k ; (H-ll) V = V Hence, from (H-ll) and the d e f i n i t i o n of J i n (9.27) and S\ El O i n (9.45), i t i s seen that j ' = J (H-12) a a and S\Q = oS\Q (H-13) The result of (H-12) i s that a value of J& i s associated with A each r a d i a l l i n e c\Q i n the i n i t i a l condition space of Lagrange m u l t i p l i e r s , and the res u l t of (H-13) i s that the rate at which the f i r s t v a r i a t i o n technique converges to the l i n e i s independent of the i n i t i a l scale factor. Using (H-6) to (H-10) i n (10.47) of the second v a r i -ation approach, i t i s found that 6 r = c r and p' = c2p' (H-14) from which i t can he determined by (10.48) and (10.42) that ' ? k . = cTk . opt opt and S\Q = cS\Q (H-15) Note that by (H-14) and (H-15), i f the neighbouring trajectory 194. X i s taken to be o XQ = \i\\Q +S\Q); (H-16) then XQ = XQ (H-17) Hence, instead of converging to the l i n e \iX , the sequence w i l l converge to the c l a s s i c a l solution \ since t h i s point has the minimum error i n transversality on the solution l i n e . The use of (H-16), however, i s merely a means of obtaining the c l a s s i -cal solution and does not provide an improvement i n convergence. Using (H-6), (H-7) and (H-ll) i n (10.27), i t can be shown that by using the modified method of matching end points then S\o. = c S \ Q (H-18) Prom (H-18) i t i s seen that the rate of convergence for the; modified method of Section 10.1 i s independent of the i n i t i a l scale factor for the Lagrange m u l t i p l i e r s . As a r e s u l t , t h i s modified approach has improved i n i t i a l convergence, and can be conveniently used with the combined algorithm of Section 10 .3 to provide the property of rapid f i n a l convergence. 195 = REFERENCES 1. Pontyagin, L.S., et a l , "The mathematical theory of optimal processes", Interscience Publishers, John Wiley, (1962). 2. McG r i l l , Ro, and Kenneth, P., "Solution of v a r i a t i o n a l problems by means of a generalized Newton-Raphson operation", AIAA Journal, Vol. 2, No. 10, October, 1964. 3. Knapp, C.H., and Frost, P.A., "Determination of optimal control and t r a j e c t o r i e s using the maximum prin c i p l e i n association with a gradient technique", JACC, A p r i l , 1965. 4. McReynolds, S.R., and Bryson, A.E., "A successive sweep method for solving optimal control problems", 1965, Proc. JACC, pp. 551-555. 5. Gottlieb, R.G., "Rapid convergence to optimal solutions using Min-H strategy", 1966 Proc. JACC, pp. 167-176. 6. K e l l y , H.J., "Method of gradients", i n "Optimization techniques", edited by Lietmann, G., Academic Press, New York, (1964). 7. Bryson, A.E., and Denham, W.F., "A steepest ascent method for solving optimal programming problems", J. Appl. Mech. 29, pp. 247-257, (1962). 8. Bellman, R.E., and Dreyfus, S.E., "Applied dynamic program-ming", Princeton University Press, 1962. 9. Chan, W.C., "The derivation of optimal control laws and the synthesis of real-time optimal controllers for a class of dynamic systems", Ph.D. dis s e r t a t i o n , University of B r i t i s h Columbia, January, 1965. 10. Larson, R.E., "Dynamic programming with reduced computa-t i o n a l requirements", IEEE Trans, on Automatic Control, Vol. AC-10, pp. 135-143, A p r i l , 1965. 11. Miele, A., "General v a r i a t i o n a l theory of f l i g h t paths of rocket-powered a i r c r a f t , missiles, and s a t e l l i t e carriers' 1, Astronautica Acta 4, pp. 264-288, (1958). 12. Bohn, E.V., Chan, W.C., and Sutherland, J.W. , "The synthe-s i s of real-time optimal controllers for rockets", Astro-nautica Acta 12, pp.' 26-32, (1966). 13. Breakwell, J.V., et a l , "Optimization and control of non-l i n e a r systems using the second v a r i a t i o n " , J.S.I.A.M. Control, Ser. A, Vol. 1, No. 2, 1963. 196. 14. Balakrishnan, A.V., and Neustadt, L.W., "Computing methods i n optimization problems", Academic Press, New York, (1964). 15. Isaacs, D., et a l , "On a sequential optimization approach i n nonlinear control", 1966 Proc. JACC, pp. 158-166. 16. Zadeh, L.A., and Desoer, C.A., "Linear system theory", McGraw-Hill Book Col., Inc., Naw^York, (1963). 17. Leitmann, G., "Optimization techniques", Chapter 6, Academic Press, New York,.(1962). 18. Sutherland, J.W., and Bohn, E.V., "A numerical trajectory optimization method suitable for a computer of li m i t e d memory", IEEE Trans, on Automatic Control, Vol. AC-11, Number 3,' July, 1966, pp. 440-447.-19. Sutherland, J.W., and Bohn, E.V., "Optimal control laws for a class of "aerodynamical systems", (Submitted for publication to Astronautica Acta). 20. Sutherland, J.W., and Bohn, E.V., "The numerical solution of optimization problems by a systematic search i n i n i t i a l condition space", (Submitted for publication to IEEE Trans... on Automatic Control). 21. Akhiezer, N.I., "The calculus of variations", B l a i s d e l l Publishing Company, New York, 1962.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Synthesis of optimal controllers for a class of aerodynmical...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Synthesis of optimal controllers for a class of aerodynmical systems, and the numerical solution of nonlinear… Sutherland, James William 1967
pdf
Page Metadata
Item Metadata
Title | Synthesis of optimal controllers for a class of aerodynmical systems, and the numerical solution of nonlinear optimal control problems |
Creator |
Sutherland, James William |
Publisher | University of British Columbia |
Date Issued | 1967 |
Description | In Part I, a method is developed for determining the optimal control laws for a class of aerodynamical systems whose dynamics are linear in the thrust and nonlinear in the lift and thrust angle. Due to the presence of the linear thrust control, a singular subarc exists along which it is often possible to eliminate the Lagrange multipliers from the control equations. Conditions under which this elimination is possible are derived, and expressions for thrust and the rate of change of lift and thrust angle are obtained that depend only on state variables and a small number of time-invariant parameters. The optimal values of the unknown parameters are determined by a direct search in parameter space for that set which minimizes the system performance function. As a result, the proposed method is considerably simpler than standard numerical techniques that require a separate search in function space for each component of the control vector. Furthermore, since the control vector is generated by the direct solution of differential equations, the method appears suitable for use with in-flight guidance computers. Several numerical examples are presented consisting of one, two, and three dimensional control. In each case, it is shown that the search in multi-dimensional function space can be replaced by an equivalent search in the parameter space of initial conditions. In Part II, a three stage numerical algorithm is developed for a general class of optimal control problems. The technique is essentially a combination of the direct and indirect approaches. Like the indirect approach, the control law equations are used to eliminate the control vector from the system and adjoint equations. However, instead of trying to solve the two point boundary-value problem directly, the augmented performance function is first considered to be a function of the unknown initial conditions and is minimized by a gradient search in the initial condition space. It is shown that it is sufficient to search over the surface of any sphere for the intersection of the line μλ*₀ where λ*₀ is the classical solution of initial values. As a result, this first approach is not dependent on a good initial estimate of the optimal trajectory, and is therefore used in the first two stages of the proposed algorithm to provide the property of rapid initial convergence. The property of rapid final convergence is obtained by employing either a modified method of matching end points, or a method of determining the optimal step size for the gradient method of the first two stages. Either combination results in a three stage numerical algorithm that has good initial convergence, good final convergence, and which requires storage at terminal points only. Several examples are presented consisting of both bounded and unbounded control. |
Subject |
Aerodynamics Automatic control Systems engineering |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-09-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0104714 |
URI | http://hdl.handle.net/2429/37266 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1967_A1 S88.pdf [ 8.97MB ]
- Metadata
- JSON: 831-1.0104714.json
- JSON-LD: 831-1.0104714-ld.json
- RDF/XML (Pretty): 831-1.0104714-rdf.xml
- RDF/JSON: 831-1.0104714-rdf.json
- Turtle: 831-1.0104714-turtle.txt
- N-Triples: 831-1.0104714-rdf-ntriples.txt
- Original Record: 831-1.0104714-source.json
- Full Text
- 831-1.0104714-fulltext.txt
- Citation
- 831-1.0104714.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0104714/manifest