ABSTRACT COMPUTATION OF OPTIMLL CONTROLS FOR NONLINEAR SYSTEMS VIA GEOIETRIC SEARCH By Richard Blain Stratton This thesis develops computational procedures for determining optimal controls for a rather general class of nonlinear systems. The procedures combine the general applicability of search techniques with the more rapid con, vergence of reachable set-oriented methods. Example prob- lems and numerical results are given for the algorithms developed. The minimum distance problem is principally considered although the time optimal control problem is also discussed. Extensions to other control problems are also possible. For the minimum distance problem, all optima lie on the boundary of the reachable set—-the collection of attainable system states for a specified final time. Reachable sets resulting from linear systems are con, vex, and the optimum final state is well-defined and in most cases unique. For nonlinear systems, however, the resulting reachable sets are, in general, nonconvex. As a result, there may be many boundary points on the reachable set which. are optima in a local sense. The global optima or optimum,a if unique, are found within this collection of local optima. Since the reachable set may not be convex, many of the pre- Richard Blain Stratton viously developed reachable set techniques are not easily applied. Thus a relatively new approach is taken. To evaluate each control which is considered, several error functions are developed which depend on the collin~ earity of the final adjoint and the final state at an opti- mum. In as much as an explicit expression for the boundary of the reachable set is not available, principles from dif- ferential geometry are used to define a path on the boundary of the reachable set. A sequence of final system states for which the error functions decrease monotonically to the op- timum value characterizes the path. Because of the fact that the reachable set is defined implicitly through the system differential equation, it is not possible to write an explicit equation for this path. Hewever, an algorithm to determine an optimum final state is developed utilizing an approximation to the path. Because of the approximate nature of this path, several alternative decisions relating to the algorithm are considered and their relationship to the error functions are investigated. Some special problems pertaining to reachable set char— acteristics are discussed and shown to be related to the global problem--that of find a global optimum. To treat the global jproblem, a random sequence of starting points are generated as a basis for each determination of a local opti- Richard Blain Stratton Several example nonlinear systems are considered and algorithm alternatives are compared. Example computational results for a variety of applications are given as are exam- ple reachable sets and trajectories. A summary of the theo- retical and computational results for the algorithms devel- Oped in this thesis is presented in the concluding section. COMPUTATION OF OPTIMAL CONTROLS FOR NONLINEAR SYSTEMS VIA GEOMETRIC SEARCH By Richard Blain Stratton A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHI Department of Electrical Engineering and Systems Science 1969 To My Wife And Our Parents 11 ACKNOWLEDGMENTS The author wishes to express his appreciation to the members of his committee for their interest, encouragement and advice. In particular, the significant and thoughtful suggestions, guidance and criticism of the committee chairman, Dr. Robert 0. Barr, are gratefully acknowledged. Special gratitude is due Dr. Leroy ll. Kelly for the inform- ative and helpful discussions of differential geometry. The author also deeply appreciates the assistance of his wife, Teddy, in typing the manuscript and her con- tinuing support and encouragement. 111 TABLE OF CONTENTS Page List of Tables vii List of Figures viii Chapter 1 Introduction 1 Chapter 2 Problem Formulation and Reachable Set Concepts 6 2.1 Notation.and Terminology 6 2.2 System Definition 10 2.3 Optimal Control Problem Definition lb 2.h Minimum Regulator Problem 15 2.5 Reachable Set Definition. 15 2.6 Pontryagin's Maximum Principle 16 2.7 Normality 22 2.8 Properties of the Reachable Set 25 2.9 Problem.Statements 27 Chapter 3 The Local Optimum Procedure 29 3.1 A Discussion of Iterative Methods 30 3.2 Initial Costate Iterations and Error Functions 32 3.3 Determination of a Local Optimmn via Direct 39 Search ‘ 3J4 Essential Concepts from Differential Geometry ’45 3.5 Convergence to a Local Optimum via Lines of 50 Curvature iv 3.6 3.7 Chapter h.l h.2 b.3 h.h h.5 h.6 h.7 h.8 Chapter 5.1 5.2' 5.3 Page Effective Boundary Path.Selection 3.6.1 Curvature Algorithm.A1ternatives-- Final Time 3.6.2 Curvature Algorithm.Alternatives-- Initial Time The Local Optimum Procedure--Composite Method h The Global Optimum and Related.Problems Special Problem 1: Nonconvex regions on.aR(T) Special Problem 2: Flats and Corners Special Problem 3: Extremum but Not Local Optimum Special Problem h: Internal Boundaries Special Problem 5: False Boundary Points Special Problem 6: Origin Interior to R(T) The Global Optimum Application of GOP to Time Optimal Control Problems 5 Compumational Results and Conclusions Example Systems .An.Introduction.to Computational Examples COmputational Comparisons of Algorithm.llter- natives 5.3.1 Study of Perturbation.Relationships-— Final Time 5.3.2 Comparison of Curvature Values 5.3.3 Effect of the Basic Curvature Formula Choice 5.3.” Comparison.of LOP-CM Subalgorithms ha and hb 5.3.5 Comparison.of Perturbation Direction Alternatives 56 59 72 7b 80 87 9o 91 93 95 97 99 103 lou 108 110 111 113 123 12h 125 Page 5.3.6 Comparison of Perturbation Step Size 126 Alternatives 5.3.7 Analog Error and the Integration 12? Correction.Routine 5.3.8 Comparison of Analog vs Digital 131 Integration.in.LOP-CM 5.h Comparison.of LOP-CM with LOP-DS 13h 5.5 Global Optimization Examples ' 13h 5. 5.1 Examples for ES—l: 2nd-Order System, 136 Scalar Control 5.5.2 Examples for Higher Order Systems 138 5.6 Time Optimization.Examples 138 5.7 Summary and Conclusions lb? Bibliography 156 Appendix A Analog Diagram for Example Problem 1 159 Appendix B Example Computer Programs 160 Appendix 0 Input Data Sets 171: Appendix D IPossible Extensions and Future Investiga- 175 tions vi Table 5.1 5.2 I.IEiT (1F TALBliEE§ Comparison of various Analog Estimates of Curvature Digital Estimates of Curvature for a 2nd-Order .System 5.3 5.1: 5.5 5.6 5.7 5.8 5.9 Digital Estimates of Curvature for a 3rdp0rder System Evaluation of the Significance of on Forward—Reverse Time Integration.Error Evaluation Application of GOP'to ES—l with x0 = (-1o,_5)T Application of GOP to Higher Order Systems Application of TOP to ES-1-, xo Application of TOP to ES-l, xO vii (-109'5)T (5"5)T Page 116 118 120 122 13C 137 1’40 1143 lh5 L128]? ()F FTIGIJRIBS Figure 2.1 2.2 2.3 3.1 3.2 3.5 3.1: 3:5 3.6 14.1 h.2 u.3 11.14 MS u.6 5.1 5.2 5.3 5.1: 5.5 5:6 5.7 Outward Normals for Various Subsets and Sets Endpoint Terminology and the Reachable Set R(T) Singular Trajectories, Flow Chart for LOP-D8, Final Time Perturbation Decisions Final State Perturbations: éxz-cx. Final State Perturbations near x+(T) Flow Chart for LOP-CR Flow Chart for Subalgorithm h.b Optima on Concave Boundary Regions Flats and Corners on Surfaces Extremum‘but Not a Local Optimum for Cos y Development of an.Internal Boundary A Local Optimum on.an Internal Boundary Reachable Sets with Interior Origins Example System #2—-3rd Order, Nonlinear Joint VS Single (IT only) R(T) Perturbations Example Curvature Ealues in.an Iteration Sequence Example Curvature values for R(20l, xb=(-10,-5)T Comparison of Perturbation.Step Size Alternatives Analog Error Evaluation and Correction.Flow Chart Example Trajectories for a Second-Order System viii Page 11 20 2h 1&3 60 65 66 76 78 82 87 91 92 93 96 106 112 11L» 117 ‘ 127 129 132 Figure Page 5.8 LOP-CM Convergence for Analog or Digital Integra- 133 tion 5.9 Comparison of LOP-CM and LOP—DS Convergence 135 5.10 Reachable Set and Example Iterations for LOP-CM 139 5.11 Example Iterations of GOP for Example Problem 2 lhl 5:12 R(t) for ES-l, x0 = (—10,-5)T, Various 1; Inn 5.13 R(t) for Es-l, x0 = (5,—5)T, various t 1H6 ix CHAPTER 1 INTRODUCTION The traditional approach to the investigation of con- trol systems has evolved significantly in the past twenty years. The criteria previously used for evaluating systems have changed as have the approaches used to insure that a system has desirable performance characteristics. Concepts such as rise time and steady-state error are being sup- planted with performance functional, target set, control constraints, etc. Bode, Nyquist and root-locus procedures are being supplemented by various other theoretical and computational methods. The newer concepts of modern opti- mal control theory are used in addition to the classical control principles. Together with this evolution in concepts, terminology and methods, there has been a change in the computational methods utilized to solve problems. These are becoming much.more computer-oriented because of the advent of com- puters which are faster and have greater capacity. As a result of this increased use of computers in solving optimal control problems, many computational tech- niques have been presented [Tl]. These methods, though related, have some significant differences. Dynamic pro- gramming, which was developed by Bellman [B3] and others is one approach. The methods of linear and nonlinear pro- gramming [H1,Zl] are appropriate for certain static opti- mization.problems. Another major class of procedures in, cludes gradient methods [B7,Nl] and search methods [H6] which.attempt to "directly" minimize the performance func- tional. In 1958 Soviet Union mathematicians led by Pontryagin [P2] presented the important maximum principle which is a necessary condition for optimality. Its impact on control theory has been great. Many computational methods involve the determination of a solution to the two-point boundary value problem generated by the maximum principle [D1,K1,Pl]. Closely related to these approaches are those which utilize properties of the reachable system states and the related adjoint system vectors [B2,E1,G1,H5,N2]. These concepts are basic to the work of this thesis. Initially, optimal control problems and the associ- ated computational methods were applicable only to low or- der, linear systems for which time optimal or minimum error regulator solutions were desired. Recently, however, empha- sis has been.placed on more complex systems--systems of higher order and systems which are nonlinear or stochastic in nature. Hereover, an additional objective has been to increase the computation speed for determining an.Optimal control. The subject area for this thesis is the computation of optimal controls for nonlinear systems. Within this general area, several approaches have been suggested. By far the most Obvious is a linearization [H7,Ll] of the systems such that existing linear methods can be utilized. Here subtle linearizations include the methods of successive approxima- tions applied to the system, to the control or to the reach. able set [H2,K2]. Other techniques for nonlinear systems include direct, random and pattern search. The method of quasilinearization in.which time-independent nonlinear systems become time- dependent linear systems is sometimes employed [Bh,B5]. The approach to be used in this thesis combines direct search methods with reachable set techniques. A feature of direct search procedures is general applicability to a wide class of problems. By utilizing the geometrically-oriented reachable set concepts, more efficient computational pro- cedures can be obtained. The application of reachable set techniques to the optimization of nonlinear systems is relatively new. The resulting reachable sets are usually not convex, hence existing techniques for convex sets must be significantly modified. Another difficulty is the lack of useful examples and computational data for comparison. Thus, one contribu- tion of this thesis is providing data for example nonlinear problems and reachable sets. The minimum error regulator problem is primarilyexamined but, as has been shown [B1,Fl], solution to more complex problems can.be based upon the b successive solution of this basic problem. Often thesis tOpics are generated by an attempt to find the solution to a specific physical problem or by restrict- ing the system considered to a very specialized class. While such a restriction often.yields a definite mathematical structure, it results in a method which is, as expected, restricted in the area of application. On the contrary, the approach of this thesis is general in the physical applica- tions which can be treated and in the rather general form of the nonlinear equations. It should be noted that such a general approach does not preclude specific applications-- as is evident from the examples which are included. This dissertation may be outlined in the following man- ner. Chapter 2 includes comments relative to notation and basic definitions. The systems and problems to be consid- ered are also defined. The important maximum principle is introduced as are the concepts of normality, extremality and optimality. Also in Chapter 2 the reachable set is defined and related properties are given. Finally, the modified (for nonconvex reachable sets) minimum distance problem (HP) and the associated local optimum problem (LP) are defined. Reachable set computational methods are discussed in Chapter 3,_particularly as they apply to problems involving nonlinear systems. Preparatory to the introduction of an algorithm to solve the local minimum distance problem, sev- eral error functions are developed. They emphasize the collinearity of the Optimum final state and adjoint vectors. A direct search algorithm is given.which utilizes these error functions to evaluate convergence. Also in Chapter 3 various concepts from differential geometry are presented and are used in the development of a geometric search procedure for computing optima. The previ- ously defined error functions are shown to decrease mono- tonically along paths defined on the boundary of the reachp able set. Various algorithm alternatives are discussed as they relate to effective boundary path selection. Con» cluding Chapter 3 is an algorithm based on geometric search, to determine the optimum in a locally convex subset of the reachable set. In Chapter h, the global optimization problem is OOH! sidered. In this case there may be many local optima. Re- lated to the global problem are several special cases which are also discussed. A global minimum distance procedure is then presented and its application to time optimal control problems is demonstrated. In.Chapter 5 example nonlinear problems are given. Reachable sets and typical trajectories are also presented. Within the algorithms, alternative choices are compared and general computational results are given for the local and the global problem. Finally, conclusions relating to the computational results and the developments of previous chapters are discussed. ' I CHAPTER 2 PROBLEM FORMULATION AND REACHABLE SET CONCEPTS A discussion of systems and optimal control problems can.be approached in any one of several ways. One approach is to start with a basic, simple system and later extend the discussion to more general systems. In this chapter, however, the more general formulation is first introduced and then specialized as needed. System and control assump- tions necessary for future development are also introduced. Knowledge of many common optimal centrol concepts, such as target set, performance functional, etc., is as- sumed and only discussed as deemed necessary or instructive. The reachable set is defined and important related results are summarized. In the last section of this chapter, an introduction is given to the basic problems to be solved. 2.1 Notation and Terminology Let En denote npdimensional Euclidean.space. No spe- cial notation is used to distinguish between scalars and vectors. lest symbols in this thesis are vectors (for ex- ample, x, u and p); scalars are so designated as they are introduced. The components of any vector are denoted by subscripts, namely, x1, 1 a l,:-:, n. Let t be a scalar, time. Let [to,.T] denote a general time interval where to is initial time and T is final time. lhere any vector, for example x, is a function, x(t), of 6 time, the following abbreviated forms are frequently used: x(t) = xt’ (2.1) x(to) . 10 (2.2) and x(T) = IT. (2.3) Where both component subscripts and the time subscripts of Equations 2.1 through 2.3 are simultaneously used, the com- ponent subscripts are placed first and the time subscripts last (for example, xZT). To represent the vector function x(t) over an interval of time, x(o) is used. The time de- rivative of x(t), dx/dt, is denoted x(t) and 3x/ at is used to denote the partial derivative. Let [a] denote the absolute value of the scalar a. Let denote the inner product of two vectors, x and y: n <1 , y) 3 Elxiyi. (20“) 1:: Let llxl] denote the Euclidean norm of x: llxu = < )i. (2.5) and thus ”x“ represents the length of the vector x. Since x may be viewed as either the vector x or the point x in E”, “x“ also represents the distance of the point x to the origin. Similarly, [Ix-y” denotes the distance between the points x and y. Superscripts are used to denote iteration indices, 1 namely x denotes the x vector for the 1th iteration. To 8 denote Optima, superscripts are also used: x+ denotes a local optimum and x. denotes a global optimum. Optima. are also indexed, if necessary, using pro-superscripts. For 1x+ represents the it11 local Optimum. example, Standard set notation is used (I) ,C, 3, etc.) with one exception: brackets are used to define a set. For example, Y=Echn=ll7|l<1L (2.6) denotes the set of all y in En such that norm y is less than 1. The boundary of a set V is denoted aY, the com- plement is denoted Yo and the closure is denoted 1'. A neighborhood, or open sphere, with center x and radius e, is denoted N(x;e): N(x:e) = [y z "7-!“ < e]. (2.7) A set K in En is convex if for any x1 and any x2 in K, the point x3 = nxl + (1-n)x2, 051751, x3 belongs to K. A set K in En is strictly convex if for any x1 and x2 in K, the point x3 = rrx1 + (l-n)x2, O - ], p ,4 O. (2.11) The closed half-space bounded by 91x52) _w_i_t_h outward nOrmal p is defined as: ‘ Q-(xgp) a: [y e En : S ], pg! 0. (2.12) Let K be a closed, convex set in E”. A hyperplane such that K ”an) is nonempty and KCQ-(np) is called a support gypegplane 3.3 K 1111; outward normal p. DEFINITION 2.1 L_e1 p e En, let R = R(T)C En be a nonempty, compact, reachable set (defined in Section 2.5) l and let x e E1'1 be such that x e DR. Case 1 (convex surface): Li: R(x;e) for some e > O is a convex set, then p lg _e_._n outward norml to R at x, i_f_ p is the outward normal to a support hyperplane to the con- vex set R(xge) at x. Case 2 (concave surface): If; N(x:e)n R3 for some x i; -p is the outward normal to a support hyperplane to the convex set R(xge) n R5 at x. 10 Case mixed surface : .I_f_ N(x:e)n 9R is a mixed is a final adjoint corresponding to the extremal endpoint x(T) = x (extremal endpoints are subsequently defined in Definition 2.b). Note that there may be many such.p for one x or many x for one p. See Figure 2.1 for examples of these cases. The scalar signum function, sgn, is defined by: sgn (a) = 1 a > O. (2.13) [sgn (a)| S 1 a = O, (2.1h) sgn (a) = -1 a < O. (2.15) The vector signum function is denoted SGN and is defined as: sgn'(x1) SGN (x) = Z . (2.16) Ls gno' (xn L 2.2 System Definition Consider a system whose state at any time t is de- scribed by the solution.x(t), to S t S T, to the following nonhomogeneous,.nonlinear, vector differential equation: . A X(t) = f(t,x(t),v(t),W(t)), 1(to) = x0, (2.17) where: t represents the independent variable, time, x(t) e En is the state vector, x(t) is its time derivative, 11 11 p1 a. Case 1 - Strictly Convex b. Case 1 - Convex but has ”flats" and a ”corner” (See Section.h.2) i” N(x;e) ‘/ x1 p1 \R 1’2 a \ x2 / 4’ \\_/, c. Case 2 - Concave d. Case 2‘- Concave with "flats“ and a ”corner" e. Case 3 - Mixed with corners f. Case 3 - Mixed FIGURE 2.1 Outward Normals for Various subsets and Sets 12 x0 is the initial state v(t) e Em is the control vector defined on a compact interval of E1, namely, I = [to,T], w(t) e Eq is the parameter vector defined on I, f('.°,',') is an npdimensional vector function defined on I x En'x Em x Eq. In the development to follow, v(t) and w(t) are treat- ed as one composite control vector, u(t), i.e., v(t) u(t) = (2.18) w(t) where u(t) is an m + q = r-dimensional vector defined on 1. Thus Equation 2.17 becomes: x(t) = f(t,x(t),u(t)), x(to) = x0, (2.19) where f(-,-,-) is an npdimensional vector function defined on I x En‘x Er. Unless otherwise specified, the term "con. trol" will hereafter refer to the composite vector, u(t). Let U’be a nonempty compact set in Er. A measurable func- tion u(-), defined on I with range space U is said'to be an admissible control and F is used to denote the family of admissible controls. In order that the solution exists, is unique and con-‘ tinuous [A1] for all u(-) in F, additional assumptions are introduced: H1) f(t,x(t),u(t)) is continuous on I x o x U, where e is a nonempty open set in E“, 13 H2) for any x in e, u in U and t in I, f(t,x,u) e Cl, (i.e., the first partial derivative with re- spect to x is also continuous). For linear systems it can be shown that a unique, global solution exists, but for nonlinear systems only a local (unique) solution can be proven. Note that the assumption of measurability for u(-) is often replaced with the strong- er assumption that u(-) is piecewise continuous. Note also that H2) is often replaced with the Lipschitz Condition: H2‘) There exists an integrableAfunction K(-) on I such that: [I f(t,x,u)-f(t.y.u) [I S R(t) [Ix-y”, (2.20) for any x and y in 9, u in U and t in I. The stronger assumption H2) is used since it is necessary for proving convergence to a local optimum. Additional assumptions are required for some of the results. In another section of this chapter the reachable set R(t) is defined. To guarantee that R(t) is compact and varies continuously with time (hence to guarantee the gens eral existence of the optimal control) the following two conditions (boundedness and convexity) are necessary [L2]: H3) III”) II S B for any t in 1, any u(-) in F (uniform bound). Rh) V(t,x) = [f(t,x,u) : u e U] is convex for each fixed x and t. q The assumptions listed above are not excessively re— strictive. For instance, consider one of the most often used families of admissible controls, F1, corresponding to 1“ the range space U1 in which each component of u(-) has absolute value less than or equal to l on 1: U1 = [y e Er : |y1| S 1, i = l,'-',r]. (2.21) Certainly this control satisfies the requirement that U'be compact and simplies the fulfillment of Rh). Further, if the control is a vector signum function of a continuous argument, (then u(-) ‘is certainly measurable. 2‘3 Optimal Contgol‘Problem Definition In addition to the system description and the class of available control functions, the optimal control problem also includes a prescribed set of conditions (final and sometimes initial) and a performance functional to be opti- mized. The initial conditions for time and state were pre- viously given as to and x(to) = x0, respectively. The final conditions are often determined by a target set G for the problem. For example, xT e G(T) may be required where G(o) is a nonempty set in En for each t e 1. Thus the general Optimal control problem is as follows: PROBLEM’2,1: Qiygg: the system (Equation 2.19), the class F of admissible controllers, and the performance func— tional T t 0 where K(o,-) is a continuous function from E1 x 9 to E1, and L(-,-,-) is a continuous function from I x e x U to E1. 15 Find: a control function, u*(o) in F which optimizes (maxi- mizes or minimizes) the performance functional while satis- fying Equation 2.19 and the prescribed set of conditions. It should be emphasized that an optimal control u‘(o) which generates an optimal trajectory x*(-) need not be unique. 2.h Minimum Regulator Problem In preparation for more complex problems, initial attention is given to the minimum regulator problem. For this problem, given a specific final time T, a control which drives the state xT closest to the origin is an optimal control. Specifically, the problem is defined as follows: PROBLEM 2,2: Q1122: the system (Equation 2.19), the class F of admissible control functions, final time T and the performance functional J(to,T,x,u) = K(xT) = 1le“, (2.23) ,Eigg: a control function u'(-) in.F which.mdnimizes K(xT) while satisfying Equation 2.19. 2.5 Reachable Set Dgfinition In much of the discussion to follow, the concept of the reachable set is important. For example, many system characteristics are directly related to properties of the reachable set. In fact, the search for a solution to Prob- lem 2.2 may be viewed as a search along the boundary of a reachable set. 16 For each u(o) belonging to F there corresponds a trajectory, xn(o) (a solution to Equation 2.19), which originates at x0 and terminates at xu(T). DEFINITION 2.2 The reachable (attainable or ob- tainable) gel at time t e I, denoted R(t), is the set of all states which can be reached at time t utilizing admis- sible controls, i.e., R(t) a [x 3 En : x a xu(t), u(-) e F]. (2.21.) Let R(-) designate the reachable set as a function of time on the interval 1. As previously indicated, c)R(t) represents the boundary of the reachable set at time t and let DR(o) designate' the boundary of the reachable set on the time interval 1. For most problems it is nearly impossible to give an explicit formula for c)R(t). Certain general properties of R(t) are known and are described in SoCtion 2. 8 e 2.6 Pontryagin's Maximum Principle Consider now, Pontryagin's Maximum Principle [P2] and its relationship to the optimal control problems previously introduced. The statement of the maximum principle varies with the nature of the problem, specifically with the na— ture of the prescribed conditions and the performance func- tional. There are, however, several essential concepts in the description of the maximum'principle regardless of the nature of the problem. These include the Hamiltonian function: 17 H(t,1(t) .u(t) .p(t)) 1' L(t,x(t) w(t)) + (2.25) with the associated Hamiltonian differential system: i“) g ggt,x(t),u(t).p(t)) (2.26) {w(t) = - gflmflfluhpm), (2.27) where p(t) is a nontrivial solution of the differential systems called the adjoint or the costate response. In the event that L(t,x(t),u(t)) is independent of x(t) (for instance, constant, as in the case of time optimal control problems), the adjoint equation.(Equation 2.27) becomes: in.) = - gghflflmfin p(t) (2.2a). because the partial derivative of L is zero. In fact, this is the same adjoint equation which one would Obtain if L(t,x(t),u(t)) = O, i.e. if the Hamiltonian function were ”unaugmented”: H(t,x(t),u(t)) = . (2.29) Utilizing this unaugmented Hamiltonian function and the adjoint equation (Equation 2.28) above, the following theorem results [L2]: THEOREM 2,1 Consider the process given in Equation 2.19 with assumptions H1) through.H3). 1L2: u'(-) belong to F and have the response x'(o) with x'(T) on the boundary of the reachable set, R(T). ‘ghgg there exists a nontrivial adjoint response p'(o) of Equation.2.28 such that the maximum condition holds almost l8 everywhere: R(t,x'(t).n'(t).p'(t)) = u(t,x'(t).p'(t)). (2.30) where u(t.x(t).p(t)) a max R(t.x(t),y,p(t)). (2.31) yeU This theorem is proved for autonomous systems in.Loe and Markus [L2], page 25“ and is extended to nonautonomous systems on page 318 and following pages. Before discussing this theorem in relation to the op- timal control, consider the following theorem which is a general existence theorem for optimal controllers [L2]. THEOREM 2,2 Consider Problem 2.1. ‘22; the target set G(t) in En'be a nonempty, compact set which varies conp tinuously for all t in 1. ‘L21 the family of admissible controllers, F, be nonempty. Further, 121.Hypotheses H1) and H2) aPPIY. ‘Thgg there exists an optimal control, u'(-), in F, on I minimizing J(tO,T,xu,u). Before relating the results of the preceding two theo- rems, consider the following terminology. DEFINITION 2,} Controls which satisfy the maximum principle (Equation 2.30) are called maximal controls. The resulting trajectories are maximal trajectories and termi- nate at maximal endpoints. DEFINTION 2.b Controls which result in trajectories terminating on the boundary of the reachable set are called extremal controls and the corresponding trajectories are l9 extremal trajectories. The boundary of the reachable set thus consists of extremal epgpoints. DEFINITION 2,5. Controls which minimize J(to,T,xu,u), as stated in Theorem 2.2, are optimal controls. The corre- sponding trajectories are called optimal trajectories and terminate at optimal epgpoints. Theorem 2.1 asserts that state trajectories terminate on the boundary of the reachable set (extremal trajecto— ries)‘pplz.1£ the Hamiltonian is maximized. For nonlinear systems trajectories corresponding to controls which sat- isfy the maximum principle (maximal controls) do not neces- sarily terminate on the boundary of the reachable set. For linear systems, however, maximal controls are also extremal controls. A third important and well-known theorem is the follows ing which asserts that all optimal controls are extremal controls. . THEOREM 2.2 ‘Lpt the hypotheses of Theorem 2.2 apply. Let u‘(.) be an admissible control with corresponding tra- jectory x‘(.) fromxO to G(T). gppp the control u*(.) is optimal pplz‘ig it is extremal. This theorem is proven.in Lee and Markus [L2,page 310] and in.Athans and Falb [A1, page 305], among others. Theorem 2.2 states the existence of an optimal control but does not guarantee that such a control is uniquw. Be— cause of Theorem 2.3, Theorem 2.2 also indicates the 20 existence of at least one extremal control. In general, of course, there are many extremal controls. It is also pos- sible that several distinct extremal controls could generate the same extremal endpoint. The terminology relating opti- mal, extremal and maximal endpoints is illustrated in Fig— ure 2.2. Maximal Endpoints (/ \ ‘ Origin \_/ wuml 3R(T)=Extrema1 Endpoints FIGURE 2.2 Epgpgint Terminologz and the Reachable Set R(T) 21 Since the reachable set represents all possible tra- jectory endpoints, one possible means of locating an opti- mal trajectory would be to examine the entire reachable set--a prohibitive procedure in the case of higher order systems. It should be noted, however, that it is possible to examine the boundary points of the reachable set by con_ sidering all maximal endpoints. Certainly this consider- ably reduces the computations necessary to determine an optimal control, but for higher order systems, such an.ex- amination would still include a prohibitive number of pos- sible trajectory endpoints. It should also be repeated that interior points of R(T) might also exist as maximal endpoints (See Figure 2.2). Finally, several important additional facts relating to extremal controls should be stated. The adjoint or co- state variable p(T) is an outward normal to the reachable set R(T) at x(T) [L2]. This fact is important in.the de— velopment of several error functions later to be considered. It is also equivalent to Equation.2.30. In addition it is related to the transversality condition which is an addi- tional necessary condition in the event that the target set G(t) is a convex set. The transversality condition states that the final adjoint p(T) is normal (inward) to the tar- get set at x(T). 22 2. Normalit The concept of normality is briefly discussed in this section due to its close relationship to the character of the reachable set and because of its effect upon the avail- ability and difficulty in computing the optimal control. It has previously been stated that extremal controls belong to the more general class of nximal controls. Hence extreml controls must mimize the Hamiltonian. Such mimization would be straightforward if for each adjoint variable there exists exactly one corresponding control. Unfortunately this is not always the case. In fact, there my be several (say two: 310) and ui(-) ), corresponding to the same ad- joint p1(o), which maximize the Hamiltonian for an arbi- trary time interval IB-C I. Were the two controls equal almost everywhere: fil(t) '1 ui(t) 3.0., (2032) no singularity would result. If, however, 'filu) ,4 u1'(t), t e la, (2.33) where Is is finite or countably infinite, then a singu— larity occurs. DEFINITION 2,6 lgfor any solution to Equation 2.28, there are two or more controls (u1(t) :4 u2(t), t e la, I,I finite or countably infinite) which differ yet which max- imize the Hamiltonian on 1s, 113313 the problem is sipgular. DEFINITION 2,2 iffor each solution of the adjoint equation (Equation 2.28) there is one unique mimal 23 control 31222 the system (problem) is normal. If a problem is singular, optimal controls may exist, but may not be uniquely defined on.some interval Is. Such nonuniqueness may have a variety of effects on the reachp able set. For linear systems, normality is equivalent to strict convexity of the reachable set, while singularity causes “flats” on.the boundary of R(t). With a singular- ity, i.e. u1(t) # u2(t), t in Is, it is still possible that x1(T) = x2(T) (that the same maximal point is attained). For more insight into this problem, consider the fol- lowing nonlinear state equation in.which the control be- longs to the control family F1 and is separable: i(t) = i(t,x(t)) + R(t,x(t)) u(t). (2.3u) The corresponding unaugmented Hamiltonian function.is: R(t,x(t).u(t),p(t)) = + . (2.35) At any instant in time, the maximum with respect to u(t) is attainted if u(t) = 3an BT(t,x(t)) p(t)] . (2.36) If, however, any component of BT(',x(-))p(t) is zero for a finite interval of time, the corresponding control component is indeterminate, taking on any value or variation allowed for F1. Each of these various controls may, in.turn, lead to different maximal endpoints. This possibility is illus- trated in Figure 2.3. Note that even.though the adjoint trajectories all start at 36, the final adjoints may vary 2U ( ) indicates the value of the costate corresponding to the state vector. FIGURE 2 Si ular Tra ectories greatly since the partial derivative in Equation 2.28 is control dependent . According to Theorem.2.3, an optimal control is one of the extremal controls. Associated with each of the ex- tremal controls is an extremal endpoint which.lies on the boundary of the reachable set, and a final value for the adjoint variable which is normal to the boundary of the reachable set at the extremal endpoint. For such a control to be extremal, it must also be maximal. Note that once the initial state and initial adjoint have been.specified, the final state and final adjoint are also specified (through Equation 2.30, the maximum.condition, and inte- gration of Equations 2.19 and 2.28) unless the problem is singular. In this case, a whole section of the boundary of the reachable set might correspond to the same initial 25 adjoint (See Figure 2.3), but the extremal controls differ over the singularity interval, Is“ while such a “singularity gap” can occur and would complicate the determination of the optimal control, it would be readily recognizable. That is, any procedure yielding a series of extremal endpoints would encounter a “jump” between successive final states when the system singularity is encountered. Finally, it is possible that the singularity may not affect the determination of the optimal control-all extremal endpoints in a large region around the optimal endpoint are the result of normal conp trols. As a result of the above consideration, normality is not a requirement imposed upon the problems to be con- sidered. Even if it were desired, proving normality for a general class of nonlinear systems would be extremely dif- ficult, if possible at all. 2.8 Properties of the Rgachable Set . . The concopt of the reachable‘set, R(t), was needed for earlier discussion, hence was defined in Section 2.5. It is the purpose of this section to list and discuss some of the important properties of R(t). Necessary assumptions ‘ for proving these properties have already been.given. As expected, fewer results are available in.nonlinear system study than in the study of linear systems. Reachable sets resulting from linear systems can.be proved to be continue ous, compact, convex and to have a computable contact 26 function [Bl,L2]. For nonlinear systems, however, the reachable set is generally not convex (although there may be locally convex regions) and the computation of a contact function is complicated. For most nonlinear system,.how- ever, compactness and continuity have been proven [L2,Rl] : THEOREM 2.14 Consider the process given in linuatio'n 2.19 with assumptions Hl) through Hit). Let F be defined as in Section 2.2. 1339;; R(t) is compact and varies con- tinuously with time on I. Both compactness and continuity are important in prov- ing the existence of a unique optimal control. For in- stance, if R(t) is not continuous then a unique time op- timal control is not guaranteed for time optiml control problems. .. Another important observation is given in the follow- ing theorem [RI]: THEOREM 2,: L23 R(t). be the reachable set for the process given in Equation 2.19 with assumptions H1) and H2). if, IT a EMT) 392p 1,6311”) for any t in I. Stated differently, all points on an extremal trajectory will belong to the boundary of the reachable set of cor- responding time, t. Several other important properties. of the boundary of the reachable set have already been given, including: 1) x, e blur) (extremal trajectory)” xT is a maximal endpoint. 27 2) xT c R(T) —> pT (corresponding to x1.) is nor- mal to the boundary of R(t) at IT. 2,2 Problem Statements Gilbert and Barr have shown that a useful approach to optimal control problems centers around the solution or repetitive solutions to the following basic problem (BP) [G1]: PROBLEM 2,} (BP) m: K, a compact, convex set in E"; Fig}: A point x* e X such that IIx‘ [I a :1; ”x“. In the case of linear systems, an obvious candidate for K is the reachable set R(T) and the problem essentially be- comes a minimum error regulator problem. Assuming that the origin is external to R(T), the solution lieslon the bound- ary of the reachable set. Since R(T) is not, in general, convex for nonlinear system, the following modified prob- lem (MP) must be solved: PROBLEM 2,1) (MP) W: R, a compact set in E“: mg: a point r" e n, such that no" u - 113:: ”r”. As part of this modified problem, it is possible that several subproblem similar to Problem 2.3 (HP) but M in nature mustbe solved. Consider the local problem (LP): .PROBLE! 2,5 (LP) 9.1.7.4.“: 1R, a convex subsetof R, a compact set in. En: Lipid: a point 1r+ o 1R such that ‘1 + n r n = mm1 u 1.”. re R 28 ‘ Since R and 1R are defined to be compact, and since ”r“ is a continuous function of r, solutions always exist to Problems 2.“ and 2.5. It may be euphasized that R is not necessarily convex. .In addition, the following pro- perties can be shown [B1,Gl] for LP: 1) 1r+ is unique, 2) llir || .. o if and only if 0 e in. 3) For ||1r+|| ? 0, ir+ e 31R. For Problem 2.“ (MP), properties 2) and 3) are true but r* is not necessarily unique. In, the next chapter, one additional important property is proven-the collinearity of r,“ (MP) (or 1r+ for LP) with the normal to clR (or 31R) at r' (1r+). CHAPTER 3 THE LOCAL OPTIMUM PROCEDURE In this chapter various iterative methods and their limitations are discussed. Consideration is given to the important function which the initial adjoint has in the determination of the extremal endpoint (for a given.xo) and the associated normals to the boundary of the reachable set (final adjoints). It is shown that a straightforward solution to LP can be implemented by a direct search on the initial adjoints. Utilizing properties of reachable sets and principles from differential geometry, a more sephistioated iterative“ procedure is developed for solution of LP. In as much.as reachable sets for nonlinear problems are generally not explicitly defined, part of the development of this algo- rithm is based on geometrical considerations of the reach— able set and perturbation analysis. Special consideration is given to the choice of error functions which correctly indicate a solution to LP and which are based on significant properties of reachable sets. Finally, convergence is considered and the solution algo- rithm is given. Special problems arising as a result of nonlinearity will be deferred until Chapter 0. 29 3° 2.1 A Discussion of Iterative Methods The set R of the modified.problem (MP) given in Chapter 2 is not convex, hence the global optimum (r‘, the point closest to the origin) is generally difficult to compute. Methods developed for obtaining such an optimum depend on the nature of the system, of the admissible controls and of the reachable set. For linear systems, hence convex reachable sets (assuming certain conditions on the control, etc.), the iterative procedures given by Neudstadt, Gilbert and Barr [Nl,Gl and B1] are effective in solving the problem. For nonlinear systems, however, convergence to the global optimum is not guaranteed. If this were the only handicap, such.methods would still have direct application in determining local optima. Possibly a linearization of the system or a 'convexization' of the reachable set could be used to implement these methods: In the event, however, that one desires to retain the nonlinear equations describ— ing the system, difficulties are encountered in.applying the above-mentioned methods. These methods require the determination of a contact point corresponding to each fi- nal costate selected [B2]. A contact point is any point in DR(T) which maximizes the projection onto the final co- state. In typical methods for linear systems, this deter- mination is relatively easy (considering present-day compup tational equipment) to implement by the following steps: 1) Consider the desired outward normal which is in the direction of the final value of the costate. 31 2) Since Equation 2.19 reduces to: i(t) . I(t)x(t) + R(t)u(t), (3.1) I where A(t) is an‘n x n matrix and R(t) is an n.x r matrix, the costate differential equap tion is also linear and homogeneous: fi(t) . -iT(t) P(t). (3.2) Thus, given one bOundary point, pT, p(:) 18 defined on I. 3) Once p(o) is defined, u(.) (a maximal control) and x(.) are also defined by Equations 2.30 and 3.1. Hence xT(the contact point) is computable.o Not only is the resulting set usually nonconvex for nonlinear systems, but also the above listed method to solve for contact points is not directly applicable. An iterative method to solve MP would thus have 319 levels of iteration instead of one, with the additional level result- ing from the difficulty in obtaining contact points. To demonstrate this, consider the adjoint equation for nonlinear systems (Equation 2.28): T é(t) = -’§§—1"“"'““” p(t). (3.3) Since contact points are on the boundary of R(T), maximal controls are employed. Certainly the determination of p(o), u(-) and hence x(o) is possible once po is known: but it is not possible to solve for p(-) from the final adjoint (as suggested above) since the adjoint differential equation also depends on the yet unknown.x(o). In.aummary, two lev- els of iteration would be necessary: 32 1) An iterative solution of the two-point boundary value problem to determine each (p ,SI) pair corresponding to each given.(xo,pT? ir. 2) Some type of iteration.(such as the Basic Iter- ative or Improved Iterative Procedure of Gibert and Barr, respectively [G1 and B1], on.xT and pT to determine x as defined in LP. 2.2 Initial Costate Iterations and Error Functions. Consideration of the above discussion suggested tothe author that an iterative method based on.the initial co— state would be simple, most direct, yet effective. Since x0 is given and maximal controls are utilized, po is suffi— cent to yield x(o) and.p(o). Starting with an arbitrary initial adjoint, a sequence of initial adjoints can be de- termined such that the resulting sequence of final states converges to a local optimum. An evaluation.of whether the local optimum is being approached or has been.reached is based on error functions to be later discussed. Both digital and hybrid computation.methods are feasi- ble: Each approach.has advantages and disadvantages. Nu- merical integration.methods (particularly for ill-behaved nonlinear systems) are sometimes slow, but consistent and accurate if sufficient computation time is available. Hy- brid computation.techniques have improved significantly in recent years and thus represent another effective approach. Since analog components are used to determine x(-) and p(o), hybrid computations are usually faster. 0n the other hand; less accuracy and consistency can be expected from hybrid 33 computers. Overall control of the iterations, of course, is a digital function for either approach, as is the eval- uation of the error functions and selection of the next initial adjoint. In this thesis, both hybrid and complete digital computation are employed. Based on these computa- tions, conclusions comparing their relative value are given in Chapter 5. The choice of the error function is critical and deter- mines the convergence and efficiency of the method. Con. sider the local problem.(LP); one obvious error function would simply be: ’ E1(Po) = H IT H- (3.u)' If a small change inpo results in a correspondingly small change in xT (proof to be given later), a method using E1 would generally converge if the step size of changesin.po were reasonably chosen. The error function E is certainly not the only avail- l able test for optimality. The following development, based upon a theorem proven for convex sets provides an effective error function, E2, which is later used in conjunction with E1. Since the reachable set is not, in general, convex, for nonlinear problems, it is necessary to select a subset of the reachable set in the following manner. Consider the compact, but not necessarily convex, set R(T) in En. Let 1,. be on the boundary of R(T) and let R(T) be convex in a local region of xT. That is, R(xT;e) defined 3.1! R(xT3e) a: N(xT;e)nR(T), (3.5) is convex. Then M(xT,p) :2 [y in En :

s c], (3.6) where c :2 <11. ,p> is a scalar constant, can be called a local support hyperplane of R(T) at xT. The following the- orem can be applied in this case with K a R(nge). THEOREM 3,1 L_e_t K be a compact, convex set in E”, 0 e I; The); IT (3 DH is the closest point of K to the origin: ' . ‘ I] . ”IT“ <||x|l, for anyx e'K, xixT (3.7) if, m; gn_l_z i_f, there exists a support hyperplane M(xT,-‘pT) of K through the point IT such that “IT is norml to M(IT,pT), i.e., , lull-.131.) = Mar-IT). (3.8) Or, stated differently, xT and pr are collinear and oppo- sitely directed. 23391: This theorem has been proven by Gilbert [G1]; how;- ever, a different proof which provides additional insight is presented here. First sufficiency is shown. If M(xT,pT). is asupport hyperplane to K, then either Z c for any x in K, (3.9) or ' f c for any x in K. . (3.10) Since 2 0 for any x in K, (3.11) 35 or . S 0 for any x in K. (3.12) Thus if -xT is to be norml to the support hyperplane, then either . _ (“T , (x-xT) > Z 0 for any x in R(szc), (3.13) or _<-xT , (x—xT) > S 0 for any x in R(szc). (3.11)) In fact, itwill be shown that Equation 3.114 is true: Z 0 for any 'x in R(ngc). (3.15) Since x.r is the closest point to the origin, 1|qu - IIITIIZ : o. (3.1.) Since R(nge) is convex, it is known that for any x1 and x2 in R(nge), x3 defined by: x3 = nxl + (l - n)x2, 0 S 11 S l (3.17) is also in R(xT;e). In particular, let x1 be an arbitrary x 5‘ IT and let x2 be LP. The resulting x3 is in R(xTz’e), but is not (of 0) the closest point to the origin. Define, for0$u51,_ _ . . ., . f(n) éllnx + (l—n)xT|I2 - nxTu2, (3.18) or _ . . _ _ _ ._ 1(a) = <(rrx+(l-rrxT)) ,(nx+(1-nx.r))> - nxTnZ. (3.19) Note that ‘ . f(fl') Z 0 (3.20) and that , f(n) = o if and only if 11 = o. (3.21) 36 Therefore,f(u)'has a minimum on the interval [0,1] at the endpoint, n a 0, consequently: iii-ELL), a 0 z o. (3.22) Differentiating Eqmtion 3.19 yields (if '11 ' ' - * 73111—1 =< (113: + (1-")xT) , (x-xT) 1’ + < (x-xT) . (nx+ (l-n'kT) > , (3.23) or 9,11%“)- : 2< (nx + (l-n:)xT), (x-xT)> , (3.214) and _ df‘ ' .. .. . , .2 #Ifl280_2 (3 5) But by Equation 3-. '22. . . 2 (IT , (1.1T) > Z O, (3026) or .3 . , .. Z 0. (3.27) which is precisely Equation 3.15, thus the proof of sufficiency is complete. Necessity is easily proved. Since M(xT,pT) is a sup- port hyperplane and since the origin is a point external to the hyperplane, there is exactly one line through the origin which is normal to the hyperplane [B6]. But this is the shortest path to the hyperplane, hence to the set I, since M(xT,pT) is a support hyperplane. Thus the proof of Theorem 3.1 is conplete. Using the results given in Theorem 3.1, several error 37 functions can be developed. They are all based on the fact that at a local optimum.(closest point‘to the origin.in a convex subset of R(T)), IT and pr are collinear and oppo- sitely directed. One obvious candidate for an error func- tion isthe following: ‘< .P >‘ I '. E2(po) = 008 Y 3 IT , T (3028) Hz." up." 1 It should be noted that at a local optimum, cos y is -1 since x1. and pT are collinear. Much attention is later given to this error function; in particular, it is shown that cos y monotonically decreases. to -1 along a (yet to be defined) path.on.the boundary of a locally convex sub- set of the reachable set. It is also possible to consider other error functions. They are, however, just modifications and combinations of E1 and E2. Since both E1 and E2 are a minimum at the opti- mum, one possible combination would be the product of the two, thus compounding their convergence: 3 < p > E3(po) " ET, I ’ (3.29) ” PT (I with the optimum occurring at the minimum negative number. While E3 should compound the convergence rate, it is less desirable than 1!:2 from one standpointu-it does not approach a specific value at the minimum. Converting E2 such that it approaches zero gives: 38 < 13 > ' ' typo) =~ 1” L H, (3.30) andthe corresponding compound error fmction would be: < ,p > . . Esme) .. 33—J— + ”IT". (3.31) . . . “PT“ Multiplying by ”pr”, another version is: E693) " “T :PT> + “IT” llpfllo . (3.32) Of course, many other combinations are possible, but those listed above represent the most convenient forms. Examination of the error functions indicates tint there is a possibility of erroneous minima if either xT or. p1, equals 0. If xT‘is zero then E2 and El: have the inde- terminant 0/0 form. If pT is zero then E2 through ES all have the indeterminant 0/0 form and 156 is zero even though 1:1. is not necessarily an optimum. Thus E1 must be. used in the. event that PT is 0. 0f the compound error functions (E3, E5 and E6), E5 is preferable to E3 because it, equals mro at the optimum and to E6 because E5 is not dependent upon the magnitude of the final adjoint.“ In most of the computations described in Chapter 5, the error function E5 is used.- In summary, several error functions have been intro- duced which can be used in an algorithm for solving LP. The suggested iterative procedure is: 1) Pick an initial p0. 39.. 2) Integrate Equations 2.19 and 2.28 using a maximl control (Equation 2.30) until xT and pr are available. , . 3) Compute the value of the error ftmction, E. ’4) Change p0 (by some method yet to be determined) . such that E is improved. Continue until the error function indicates (forexample, E5 =- 0) that an optimum has been determined. _ Now that error flmctions are available to differen- tiate between local optinm. and other points, it is imor- tant to consider the precise method for changing Po (item ’4 of the suggested method) such that E is improved. Two methods are consideredr One is based upon a direct search- on the initial adjoints, the other upon geometric consider- ations of the reachable set. These. methods and their con- vergence are discussed in‘the following sections of this chapter. 2,} Determination of a Local Optimum via Direct Search Any direct (or pattern) search a'technique presupposes that the change in the controlled parameter (Do) can be made sufficiently small such that the resulting change in the evaluation parameter (E), is correspondingly small. For linear systems it. is evident that a sufficiently snmll . change in E can be obtained. For nonlinear systems, how- ever, it is not so evident. ‘ Given a specified bound on the change in PT and LT, it must be shown that the change in pa can be made suffi- ciently small to keep the change in PT and IT within this no bmmd. To prove this result for nonlinear systems, an embedding theorem by Hestenes' and Guinn is utilized [H3]. The theorem is somewhat more general than the embedding theorem given in Hestenes' book [Hit]. THEOREM 2,2 ,_I._g;t_ x( ‘), p(o .) and u(-) be defined by Equations 2.19, 2.28 and 2.30 withcorresponding initial adjoint p0. 933-213(1)“ represent anyof the error fimctions previously defined (Equations 3.’-), 3.28-3.32). The}; for any a > 0, there exists a 6 > 0 such that I E(po) - R(pg) I < (a. (3.33) whenever _ ' II p, -.p,', ll < 6. (3.31)) 259231“ e be given. Since each cf the error fimctions being considered is continuous in IT and pr, there exists a 51 and a 52 such that . . ll pT - pT ll < 61 + (3.35) lle .. 1., n < :2“ (3.36) imply Ineqml-ity 3.33 is satisfied. Now consider Inequal- ity 3.35. The relationship between po and pr can be writ- ten as p,-'l._-=:‘I'(T,t0)p0 (3. 37) where 9(T, t o) is the fundamental untrix for the adjoint system (Equation 2.28) which satisfies: $11" to )-- 13?“. t ). (3.38) Since ‘1’(t,to) is nonsingula‘r, there exists a 63 such that #1 ll po-pf, H < 63 (3.39) implies . . . n pT-pé n < 61. (3.1m) Since hypotheses H1) through H3) are sufficient assumptions for the embedding theorem [113], there exists a 514 such that H goo-p; II < a), (3.1.1) implies . . I] xT—xé. I] < 62. (3.152) Now let 6 be the. snaller of 63 and 6,4 and the theorem is proved. The fact that snail changes in p0 result in small cmnges in IT is necessary for the following direct search algorithm:- LOCAL OPTIMUM PROCEDm-QIRECT SEARCH (LQP-IB) A. f Mpg). Evaluation: Whenever Mpg) is to be evaluated, then these steps are followed: ‘Given p}, and x0, ' initial conditions, integrate Equations 2.19 and . 2.28 utilizing a naxinal control. From the values obtained for z} and p1, evaluate Mpg) by means of Equation 3.31 (E5(Po) unless otherwise indicated. B. Initialization: Choose a step size h for the'com- ponents of p , a final stopping tolerance E , an improvement factor Rt) 1 and a mximum number of allowed iterations, I . Set 1 a 0 and select an arbitrary pg; then.evaluate E(pg). If E(pg) s E ~ then x51!- 1? is "a suitable approximtion to an opti- mum 14;. C. Iterations: Define the vector 6 whose 3‘13 com— ponent is h, all other componen s being zero. Let 3 ‘3 10 ‘ ‘ 1. Evaluate E(pg + i3) and E(pg — 63). 2. Let p:(i) = 3. 5. 7. 1.2 p: + 63 if E(p1 + a j) < E(p1) pg-a, 11' 3(1):“ )>E(p1) and E(p1-6 ) - 0. (3.149) The preceding definitions are given.for curves in.genp eral; now consider curves on.a hypersurface. The hyper- surface itself has a tangent hyperplane (assuming smooth. ness) and an.associated unit normal N to the hypersurface. The two normals (n, to the curve and N, to the hypersur- face) do not necessarily coincide. In.fact, two of the most important curves on.a hypersurface are defined by the behavior of the normal to the hypersurface as related to the tangent and to the normal to a curve on.the'hypersur- face. DEFINITION 3,2 For any curve on a hypersurface the curvature vector is d2x(s)/ds2, which can further be #8 represented as a combination.of the hypersurface normal and the vector w: dzxés) g kn(s) g kn N +1. (3.50) ds where kn is the normal curvature. One important curve on.a hypersurface is the geodesic. In differential geometry it is defined in the following manner: DEFINITION 3,8 A geodesic is a curve on.a hypersur- face for which w s 0, i.e. for which the principal normal to the curve coincides with the normal to the hypersurface. Thus the curvature k and the normal curvature kh are equal. Generally speaking, geodesics (on.hypersurfaces) are analogous to straight lines in Euclidean.space and are - curves of shortest distance. Since geodesics are curves of shortest distance, it was desired to seek a modificaa' tion of LOP-DB so that the curve L(xg,xg) would be a geo- desic. Thus 144,4) would represent a shortest path from the arbitrary starting point xg on.the boundary of the reachable set to :3. Attempts of the author to achieve such a modification have thus far been unsuccessful. A second important class of curves on.a hypersurface are the lines of curvature. Such.a line is determined by considering the rate of change of the hypersurface normal. In.general, one can.write '%g a "k %% +1; v (3.51) “9 where v is a unit vector orthogonal to dx/ds 9.2.4. contained in the tangent hyperspace at the point under consideration. The factor “rs is called the torsion of the hypersurface in the direction of the tangent (dx/ds) to the curve. DEFINITION 3,2 Consider curves which lave a direction such that “(as 0: such directions are called principal directions (of curvature) and the associated curvatures, k, are called princiml curvatures. DEFINITION 3,10 A curve on a hypersurface whose tan- gent at every point is along a principal direction is a line 2; curvature. Lines of curvature thus have the property that the rate of change of hypersurface normal coincides (in direc- tion) with the tangent to the curve on the hypersurface. This is expressed in Rodriques' formula: dN +. k dx = o. A (3.52) Further, concerning lines of curvature and principal directions, it has been shown [G2,N1] that: l) A point where the principal directions are wholly indeterminate (all principal curvatures have the same value) is called an umbilical point. A hyperplane and a hypersphere (or portion thereof) are the only hypersurfaces whose points are all umbilical points. ‘ 2) Except for the two instances mentioned in item 1), the principal directions always exist and define an orthogonal system if, directions. 3) The principal directions represent directions of extreme values 2; curvature. 50 Conver ence to aLQcal Optimum via Lines of Curvature In Section 3.3 a direct search technique on the ini- tial adjoints was shown to result in a sequence of extremal endpoints (defining a curve, L, on the boundary of R(T) ) which converges to an.approximation.to a minimum of E(p°). In this section, the definitions of the previous section are combined with the concepts of extremal endpoints and orthogonal final adjoints to formulate a different converg- ing sequence. It is shown that for a convex set, there exists a well-defined curve, L(xg,x;), consisting of lines of curvature, along which the error function monotonically decreases to the optimum. This is shown using E2(p°) or cos y. THEOREM’3,3 ‘Lgl R inEn be a strictly convex compact set with 0 s B. Let x: be an arbitrary boundary point of R and let_x; be such that llell > lleII for any xT .4 x; in n. (3.53) ‘Thgp.there exists a.path L(xg,x;) on the boundary of R, consisting entirely of lines of curvature, such that the error function <1) > . E2(po) a cosy: ' T 1:11 ' (3.511) llpTll ”Ii-ll monotonically decreases to -l and such that ”sq” meno- tonically decreases to its minimum.va1ue x}. Proof: Note first that B may represent the entire reach- able set R(T) if the system being considered is linear, or it may represent a convex subset of R(T). The boundary 51 points in these cases would represent extremal endpoints and the corresponding outward normals, pT, would represent final adjoints. As previously introduced, ds represents the in. cremental arc length along L(xg,x;). To show the monotonic— ity of cos y, it is sufficient to show that LII—3&1; < O (3o55) along the specified path. Since pT is an outward normal to the surface at xT, utilizing the terminology of the previous section yields ' , (3.56) ”PT“ As a notational simplification, the subscript T will be Ni deleted. Thus Equation 3.5b becomes 1:,N3» E2(po) 3 OOBY '13 frr'x—rr—‘o (3o57) ) Taking the derivative of Equation 3.57 one obtains: dcosxagllx‘ dxds ds .. . .60 as as 'U‘TII ’ ((xu (3 1 To this point, the path has not been identified more pre- cisely than.belonging to the boundary of B3thus it is new 52 assumed that L is composed entirely of segments of lines of curvature. Hence imitation 3.52 applies and Equation 3.60 becomes _ g—c-gg—l z <..k-d—1— x >-o (3o61) a. 18’qu ’ len Now consider Thus (1 2 Or Md d: a rr—nv ; (X’s-J?) . (3o6’4) Substituting Equation 3.6’-) into Equation 3.61 gives: Leg—1-1: 41.111- da< ,m> (3.65) Hence , d cos 1=—(k + cos ) glx || . (3.66) (18 Since the curve, L, is on the surface of a convex set, the curvature k is always negative. At or near the optimum on any reasonably behaved surface cos y is negative; thus for the derivative of cos y to always satisfy Equation 3.5!). it Elf—Ll < o. (3.67) Not only is this requirement necessary for the proof of the is necessary that theorem, it is also desirable from an understanding of the problem since the optimum point is the point of minimum 53 [I x I]. It is now shown that Equation 3.6? can be satisfied while simultaneously renining on lines of curvature. The hypersurface an is m-dimensional where m < n- 1. Through each point x on an there are m orthogonal direc- tions defined by the tangents to the lines of curvature through x. Although these vectors are orthogonal, in the proof to follow, it is' only necessary to require that they are linearly independent. Consider the arbitrary starting point x0 ,4 x+ and let L(x°,x1) denote a path on all composed of lines of cur-_ vature, along which llxlI decreases monotonically from [I x0“ to ||x1 II. It is first shown that there must exist one such path, i.e. that it is possible to move from x0 along a line of curvature such that “x” is decreased. Define level hypersurfaces, P(r), on 9R as the inter- section of all with m—dimensional hyperspheres S(r) of radi- us r and centered at the origin. The resulting level hyper- surfaces are of dimension m—l, thus at least one of the tangent vectors for any xix... belonging to one of these level hypersurfaces must intersect that level hypersurface. Thus it is possible to move along the lines of curvature corresponding to this tangent such that [I x [I is decreased. In the event that r = ||x+ll, the intersection of S(r) with :33 yields a point, x+. Consider the collection of all paths, L(xo,x), formed of lines of curvature for which ”x” is decreased. If 514 L(xo,x+) is in this collection, then the theorem is proved. If not, then there must be a lower bound greater than ||x+|| for the “x“ obtainable, i.e., a level hypersurface P(h) must bound all L(xo,x) from below. There may be many such paths which approach P(b). For each x1 there is a corresponding orthogonal system through x1 consisting of lines of curvature. Along some of these lines, ”x” is decreased. Let D(x1) represent the point x1 together with those arcs of its lines of curvature, containing x1, along which llxll decreases and for which llxll < um). Since b is the lower bound of the norms of points which are joinable to x0 by admissible paths, there exists a se- quence of points {x1} such that {II x1 H} converges to b. Since (DH is compact, {x1} contains a subsequence converging to a point 3', whose norm is b. For notational simplicity, {x1} is also used to denote the convergent subsequence. As previously mentioned for an arbitrary x 1‘ 1+, at least one member of the orthogonal system must penetrate the level hypersurface to which the point x belongs. Thus there exists an 3' e DGE) with “x" [I < b. From considerations of continuity, it follows that {D(x1)} approaches DCx'). Thus there exists a sequence {x'i}. where x'1 e D(x1), which converges to 32'. Also {llx'1ll} 1Ephroaches Hi" [I < b. Hence for some k, le'kll 0, (3o73) or . 1 . . k < . (3.7“) For the convex set of Theorem 3.3 this inequality certainly holds since k is negative. Questions relating to non- convex sets are considered in Chapter ’4. Note that Eqm— tion 3.71) allows this possibility since dE5/ds my be nega- tive even though k is positive. 3.6 Effective Beundapz Path Selection For a convex region of R(T) which includes a local optimum, the results of the previous section show tint at least one and probably many acceptable paths exist on the boundary of the reachable set. Each of these paths run 57 from the arbitrary starting point to the local optimum. These paths represent collections of extremal trajectory endpoints. Any method based upon these results would have as its goal the successive determination of these end— points, hence of these trajectories. Certainly it is not desirable to identify all of the endpoints since the path consists of an infinite number, but it is desirable to ‘ identify a sufficient number to define the path and hence locate the optimum endpoint. Were the boundary of the reachable set explicitly de- fined by a vector function, the path would likewise be pre- cisely defined. If this were the case, however, determina- tion of the path would be unnecessary since such an explic- it definition of aR(T) would easily lead to location of the local optimum either by direct calculation or through func- tional minimization. In any practical nonlinear problem, however, it is not possible to explicitly define the reach- able set. Based on hypotheses and theorems given earlier in this thesis, only the following general observations are available: 1) R(t) is of dimension r S n, where n is the dimension of the state equations. 2) 'The boundary of R(t) consists of extremal trajec— tory endpoints resulting from extremal controls. 3) R(t) is compact and varies continuously with time. 1)) R(t) is, in general, not globally convex, but will have regions of local convexity. 58 Howthen, does one, experimentally, determine the path from the arbitrary starting point to the local optimum? Because of the fact that so little is known of the shape of the reachable set, it was the author's original decision to implement some type of search method on aB(T) (i.e. among the initial adjoints corresponding to these extremal end- points). One such method is LOP-IB given in Section 3.3. Theorem 3.3, however, yields some additional insight into the nature of an effective path to the optimum. Some of these results were previously apparent from the nature of the problem: l) d“ lel/ds< O, i.e. “IT” must successively decrease. 2) The decrease in || [I should be made as large as possible (i.e., dl LIT” / ds should be minimized). Others are provided as a result of Theorem 3.3: 3) Since the path consists of lines of curvature, 91.15. 4, k 1?. = o. (3.75) '4) If there is a choice of lines of curvature, as exists in most cases, theone for which I: is minimum should be selected at each decision point along the path. While arbitrary changes in' the. pd's might. be acceptable in lower order systems, higher order systems require more sophisticated methods; hence the insight provided through Theorem 3.3 is important and should be utilized. The basic structure of such a method would be as follows: 59 O . 1) Try an arbitrary po , note E(po), p: and 1:. 2) Emlcying insight into the relationship of ’1‘ and/or pT to IT and/or pT and the interrelation- ship between 51.1. and 4691. (perturbations), change 0 IT and pT yielding IT and pr. 3) Run the system in reverse time from .0 and '5: using extremal switching. This yields an i: and a '53. Note that 320 probably differs from xo (the specified initial state). ’4) Consider a new initial adjoint p1 which is related to 53 and perhaps to p00 and to the change inxo o(x0 -xo.) 5) Assuming that p10 results1 in a better value of the error function, E(po), then the method is repeated. If the result is not better, some alternative approach must be taken. Within these five steps, there are strategies which must be chosen on the basis of good judgement as well as math- ematical development. Basically these alternatives can be classed into two subdivisions, according to the time at which they occur: in step 2), after forward integration of the system and in step ’4), after reverse time integra- tions 3 , 6 , l Curvature Algorithm Al ternat ives-F1 nal Time Consider first the alternatives available at the final time, i.e. on the boundary of the reachable set. The fol- lowing choices must be nude: 1) Should just IT or pT be perturbed or both IT and pr? 6O 2) What should be the size of the perturbation? 3) In.what direction should the perturbation.be made? b) If both xT andpT are perturbed, are the perturba- tions performed dependently or independently? 5) If xT andpT are perturbed dependently, which.one is perturbed independently and what is the rela- ticnship between.the perturbations? These decisions can be represented as shown in the decision flow chart given in Figure 3.2. It should be emphasized that these decisions must be made on the basis of little knowledge of the nature of the reachable set in the neigh- borhood of xT-any additional information.must be experi- mentally determined. The goal, of course, is a perturba- tion of xT andpT such that the error function is decreased and such that the perturbed final state 'iTebnu') and the perturbed final adjoint is normal to¢8B(T) at if. Perturb. . . ‘T “0R:.AND p'1' Dependp ' I ‘dH‘Lt nt '~v 'n‘ifly erturbation. - - ‘ - Direction Which :elationp , . . ' Independ- 8111p _i._._g Perturbation 6pT * Magnitude “IT FIGURE 2 F nal Time Perturbation Decisions 61 Consider some aspects relative to these perturbations and the above listed alternatives. In the discussion to follow, the perturbed final state and costate are desig- nated if and if: ‘ i1. = IT + 61.1. (3.76) ii, = Pr + apT. . (3.77) For any perturbation in xT, there is a good possibility that 3% does not lie on the boundary of the reachable set. Of course, the goal of any perturbation.is to minimize this possibility. In.addition.tc this fact, however, other as— pects must be considered. If only IT is perturbed, stpT probably does not represent a correct outward normal for the new boundary point. Similarly if only pT is perturbed, 55 does not represent a correct outward normal for the unperturbed final state. If both.xT andpT are perturbed independp ently or if an incorrect perturbation.interrelationship is utilized, again the resulting perturbed normal may not be accurate for the resulting boundary point. Thus a judi- ciously chosen, related perturbation is most desirable. Even.in.this case, some departure from the reachable set andthe proper outward normal can be expected, particularly if the step size is too large. With careful control of the step size, however, a joint, related perturbation may re- sult in a new ET close to aR(T) yet allow 6:1. to be large enough to represent a considerable improvement in.the final 62 state. The interrelationship between IT and pr is chosen to coincide with the theoretical developments for lines of curvature, namely, Equation 3.52 applies, or, in incre- mental steps: 6N + k 6xT .. o. (3.78) Obviously, I: must experimentally be estimted and updated as the path moves towards the optimum. The final state, 1.1,, is chosen for the independent perturbation since the proof, of Theorem 3.3 requires (and common sense dictates) that ”IT” continuously decreases. Thus, choosing ml. to be independent allows easy verification of this requirement. The crucial question remining unanswered is number , 3)-—in what direction should the perturbation be made? While this decision must be made without an explicit ex— pression for the reachable set, some facts are available: 1) ”le] must decrease. 2) E(po) must decrease. 3) The final perturbed state, 21., should rennin on the boundary of the reachable set. To give justification for the direction of perturbation to be selected, consider the following perturbation analyses of E(po) which are presented as theorems. The first per- turbation direction is based mostly on fact number l)-- decreasing || IT [I . 63 THEOREM 3,14 Consider the error function E2(po) \(the final time subscript, T, has been deleted for notational simplicity): E2(po) = 0087 - We Sfi-ffiz (3.79) where x belongs to a convex region of an”) (thus the cur:- vature, k, is negative) and N is norml to R(T) at x. p23 ox be a perturbation in x with corresponding 6N given by Equation 3.78. .1; 6x is defined by 6x a -cx (3.80) where c is a positive constant, 0 < c < 1, than E(x+ 61.p+6p) S E(x,p) a E2(po) (3.81) with equality only at the optimum. Proof: Let as . sum) - E(x+6x,p+6p). (3.82) It must be shown that 6E > 0. (3.83) Consider <‘x N) 6E3 x - 1+6! , (Beau) or, =__<:x 11> 6E ”x” - ”xi-bx” x-‘+6x _ _ <5: 6N> nit-3?” "fin-[1' (3'85) Substituting Equation 3.78 into 3.85 gives: '1' ' l ' <6x N> ‘E = (x ’N’ (U31? ururu)‘ [FT-07.5?“ + ("x k6x> + <5: Rex) . (3.86) x+6x ' ||x+ox|| 614 Note that . _ A A ”1+5!" = III-01H = ll-Ol Hill ‘3 ”-0) “III (3.87) because of Equation 3.80 and since O< c< 1. Substituting Equations 3.87 and 3.80 into 3.86 results in: . _ ' ' 2 ' ‘ . 6E as cos y (€130) + I—f—ocos y - 19%|!” + fill”; (3.83) or, or = afigflwn -.- -ckl|x|l . (3.89) Since 0 is positive and k is negative, Inequality 3.83 is satisfied and the proof is complete. Note that maximizing [kl will mximize 6E. , It should be noted that Equation 3.80 disregards the fact that x+ 6x should also lie on the boundary of the reachable set. If cos y is near zero (i.e., far from the optimum), ex is approximately tangential and Equation 3.80 is a good approximtion (See Figure 3.3a). On the other hand, consider an 1 near the optimum (Figure 3.3b). In this event, any 6x defined by Equation 3.80 is nearly nor- ml to aB(T) and thus is a poor choice. Hence, some other choice must be made for ox based on the available informa— tion which includes past knowledge of 6B(T'), x(T) and p(T) and current knowledge of p(T) and x(T). When IT is near x; 6x is nearly orthogonal to IT. For second order systems, this would be sufficient infor— mation for calculating 6x, but for higher dimensional sys- tem, this still does not adequately define 6x. Considera- tion of pT and IT near an optimum shows that they are 65 origin or gin a. Far from Optimum b. Near the Optimum FIGURE Final State Perturbations: bxs-cx nearly collinear, thus the veetor u a -0. ‘fi'fi Am) , (3.90) is nearly tangential near the optimum. This fact is illus- trated in Figure 3J4, and leads to the next perturbation analysis as given in Theorem 3.5. THEOREM 3,5 Let the error function be defined by Equation‘ 3.79 and 1;; Equation 3.78 apply. Le; x and N be as previously defined. l1 6x “is defined by Equation 3.90 with o < c' < ”x" /2 , j_l_1_(_5_1_1 Equation 3.81 is satisfied. 3391: Consider 6E as defined in Equation 3.82, and as ex- pressed in Equation 3.86. Substituting Equation 3.90 into parts of Equation 3.86, one obtains: ' 66 ‘l' ' l ' 6E= (le||-||x+3x||) + ||x||||x+3x||+ ”xi-Ex” _ '' _ <‘x kc'N> IIx+6x|||Ix|I x+ x _ 2 ’ _ ' . . kc' I. 2 I >+‘(|lel'llx+6xll)+|Ix~'+61ii+l_lx+oxll' x+ x _ ke' x 2k 0' 2 .. . , 2 nx—‘S‘L—x-HOOB Y 4' 1+5! (1 008 Y) (3 9 ) Rearranging Equation 3.92 and collecting term yields: or . a... u-" ,3, )+°—'é}117f,%firfl(l-Idlxll+2h')o (3.93) To relate the magnitudes of the bx's in this theorem and the previous theorem, define a new constant 0' 3 fi" e (3o9u) FIGURE h Final State Perturbations near :3 T 67 Note that the hypothesis :for the theorem requires that 0 < c' < i. With substitution of c', Equation 3.93 becomes: 1‘) 6E .n'ces y (l - 1+; + 0" gmxlms (1-1.1121) +2kc'llxlI ). (3.95) Now it must be proven that GE is positive. Note tint the second term is positive if . . l .. kllx||+ 2kc'llx[|>0, (3.96) since cos y 2 -1 with equality only at the optimum (hence the second term is zero at the optimum). Inequality 3.96 can be rearranged as follows: . 2kc" [I x” > k” x" - I, (3.97) or since k is negative (convex surface), at l‘é‘fi'lrfi'r (3.98) Certainly if m is less than t, Inequality 3.98 is satis- fied since -l/2k|| x” is positive; hence the second term of Equation 3.95 is positive except at the optimum where it is zero. Now consider the first term of Equation 3.95, for which it must be shown that . x . _ cos Y (1- 1.“: )5 O- - (3.99) It has previously been stated that cos y is negative in any well-behaved region near a local optimum: thus one must show that . (1 -' x ' )503 (3.100) 68 or . . . llxflxll - “III S 0. (3.101) which was already the objective of this method. With the introduction of c", Equation 3.90 has become 6x :2 -c"x - c'lltxll N. (3.102) Consider ‘ _ . . . . llxflxl) =ll(1- 9")“ O'IIIIIN)". (3.103) or . e . . . [I x-I-bxll =||'(x(1-c") - c" [I xll N)“. (3.101!) Using the triangular inequality one obtains: III“!!! 5 "1(1-0') II + ll(°"llxllN)l|: (3.105) or . . . . . . [Ix-tax“ S (l-c')[l:x|| + c'llxll _ _ (3.106) since a“ < é (thus l-c' is positive) and since UN“ is 1. Thus Inequality 3.101 is satisfied and the proof of the theorem is conplete. Sumrizing, two candidates have been selected for perturbing the final state, as given in Equations 3.80 and 3.90. Equation 3.80 is most effective at points where cos y is near 0 (far from optimum) and Equation 3.90 is most ef- fective for points on bR(T) which are near the optimum (cos 7 near -l). Letting 0 5 c". (3.107) an obvious candidate for a cogosite choice for ax is 6x = -c(x-[Ig.[|N cos y )‘, (3.108) because it reduces to Equation 3.80 when cos y is O and to 69 Equation 3.90 when cos y is -1. This perturbation choice is considered in the next theorem. THEOREM 3,6 1:21 the error function be defined by Eqmtion 3.79 and 131 Equation 3.78 apply. _1_._e_t_ x and N. be as previously defined. .I_!_ 6x is defined by Equation 3.108, 0 < c < 1, M Equation 3.81 is satisfied. 6 PM: Again consider 6E as defined in Equation 3.82. By partial substitution of Equation 3.108 into Equation 3.86, one obtains: 1. . l . o. GEI (m-m)+m .. 0 t ‘_ ck' amn<-l"””°°"'“’ m“"’ ck + W< x , [I xHN cosy > + "-31%” GIII2. (3.109) Since ' 3 [lel cos y a +flx|lzcosz Y) (3.112) ’1' A . _ . . . . , u '6: u? . «:20! x ()2- u x llzooezv) - «221)qu (1- coszy). (3.113) Substituting this result into Equation 3.111 gives: 6E=cosy(1- )- x (1- cos 27) x+Ox ext-x kc + 1+5: 22(1-cos v): (3.111)) 70 or or: a cosy (l-- )+ 2(1-0082y)(o-1). (3.115) x+6x c,x-I-6x While the first term of this equation is the same as de- rivod in Equation 3.95, it must again be considered since 6x differs. For that term to befpositive (or zero) it must be shown that . . . . IIxII 2 11:11:” . (3.116) From Equation 3.108 we can write, [I x+6x II - |I(x-cx+c II x“ Recs 7 )II =IKx(1-c)+oIIxIINcos 7)". (3. 117) Using the triangular inequality one obtains: llx+61llllx(1-°)ll + II (cllrllNoos (Noll .(3-118) Since Icos yI< .. l, c and l-c are positive; and IINII-l; IIx+6xII S (1-c) ”x“ + c IIxII a IIxII. (3.119) Thus the first term in Equation 3.115 is positive or zero (at the optimum). Now consider the second term. The factor (c-'l) is negative as' is It: all other factors are positive hence the term is positive except at the optimum where it is zero; hence the' proof of the theorem. Experimental results are given in Chapter 5 to compare the three possible methods of perturbing, 6x as given by Equations 3.80, 3.90 and 3.108. Once 6: has been defined, then on is determined through Equation 3.78. Once k is determined then 6N is specifically defined. The deter- mination of curvature, k, however, is not an easy task for higher order system. Since R(T) is not explicitly 71 defined, 1: must be experimentally approximted. Given two neighboring points on a line of curvature, (x,p) and (-x',p'), several means are available for calcu- lating an estimte of k. The first estimates are apparent from Equation 3.78 where n (dimension of the state) approx- imations are possible E1 1: — ———L:i- :' , i-l,n~,n. (3.120) 1" i any of these estimates could be used or the average: 11 KT ‘3 (1211(1) / no (3o121) A second estimate of k is available from Equation 3.66 which, solving for k gives ks-W-fi. (3.122) Or, in terms of perturbations Finally, a combination of k and Eav say be used: iavz = i (PM, + E) . (3.1214) For two final statee x and x' close to each other (and on a line of curvature), all of the estimates of k (Equations 3.121, 3.123 and 3.121)) are near the actual value. For a second order system, with 6x sufficiently shall, k can easily be estimted since the boundary of the reachable set is the line of curvature. In higher order systems, how- ever, a perturbation of (x,p) yielding ,72 E, s E s E, i=1,”‘,n—1 (3.125) itl would indeed be fortuitous. In general one can expect that the ki's and k have differing values and differing signs I since the. perturbation of (x,p) my not be along a line of curvature through x. For small perturbations in x, the variance of the values for k and the 151's is certainly an indication of whether or not the perturbation is along a line of curvature.” If the variance is large, several op- tions are available in selectingk: an average my be chosen, one particular estimte formula may be relied upon or more perturbations may be taken until a better estinte of k re- sults. The methods for estimating. k are experimentally compared in the next chapter. 3,6,2 Curvature Algorithm Alternativeso-Initial Time As indicated in Section 3.6, after reverse time inte- gration from (2:55;), there are several alternative means “1', the initial adjoint for the next itera- o tion. Since 6x}. and 6p&. are only estimates for an accurate of choosing p perturbation on 3R(T), the resulting 1; found by reverse time integration of the estate and costate equations using 1 " maximal switching often differs from the specified initial state, x0. In addition, computational errors may develop which also contribute to the difference between :3, Indeed, experimentation indicates that significant errors and :0. do develop and that a ccmputated (iT’fiT) pair when used, without perturbations, as initial points in a reverse time 73 integration, may not produce x0 and p0. The origin of this computational error is discussed in éhapter 5. For the de- velopments in this section, it is sufficient to observe that a computed x; may not be the initial state because of: 1) Computational errors. 2) Failure of i}. to lie on bB(T) and/or failure of 5.}. to take the correct direction. Regardless of the source of this difference, when a new pi“ is selected, based upon '13:, the validity of 3: should be examined as well as the possible benefit of adjusting for 1 -1 6x0 = 10 - Io. , (3.126) Disregarding these differences, a new approximation for Po is i+l -i‘ po 3 Poe (30127) Consider the computational errors in integration which might develop. Ons’ possible means of reducing their effect would be to consider the following error correction vectors: (:1 a i: — 10 (3.128) 31 a p: - pg, (3.129) where 201 and 15; are initial points obtained by reverse time 1 T and 31 thus re- integration with extremal switching from. p i and I; (1.0., without perturbation). The differences a sult entirely from computational errors. These correction factors may then be used to give new estimates of ‘0 and pa: a}, = if; - «1 (3.130) and . 7'4 1 ‘-i i 50 = p0 - a . (3.131) These new, hopefully better estimates can then be used to determine pg“. The effect of these correction factors is evaluated through experimentation in the next chapter. It should be noted that their computation at each step re- quires an extra integration of the state/costate equation pair. Once i; and $3 (or E; and 5;) have been computed, then a new pgfl'must be calculated. is mentioned, a straight- forward method is to pick pf,“ = 15; (or 33,). (3.132) Other alternatives are available such as: pf,“ = p}, + ow; - pg), o > o (3.133) where c < 1 if ox: is large and c z 1 if the estimates are EOOdo 2,2 The Lgcal mtimum Procedure-moogposite Method Incorporating the results obtained in the previous sections of this chapter, an algorithm can be given to solve Problem 2’.5 (LP). Included in this algorithm are several alternatives. Some are given in the form of sub- algorithms, others are evident by the choice between sev- eral alternative equations. In Figure 3.5 a flow chart is given for LOP-CM and in Figure 3.6 a flow chart is given for an example subalgorithm (hb). 75 L00 OPTIMUM PROCEDURE-~COMPOSITE METHOD LOP-CH l. 2. 3. h. 5. 9. Select an error function Select an arbitrary pg, (1 = O ). Determine x;, p; and E(p:) corresponding to xb, p: and extremal switching. Determine an estimate for k through Subalgo— rithm b.a or h.b. Perturb x% by using Equation 3.80, 3.90 or 3.108 yielding 'x',}_.), where the magnitude of 0 (step size) may be dependent on the error of the i-l'st iter- ation. Perturb p: in.a corresponding manner (through Equation 3.78), giving '51. Determine E: and.§: corresponding to extremal switching,'§%, E} and reverse time integration. Evaluate 32'1 and ii: and determine (through one of c the Subalgorithms 8.a through 8.0) a new'p:+l. Determine x$+l, p%+l and.e(p:+l). Test to see if the error value is decreased. If so, test to determine if the error value is sufficiently close to the optimum value. ‘If it is sufficient- ly close, then a local optimum has been.found: if not, go to step b and repeat. If E(p:+1) is net an.improvement, and if the step size has not been reduced past.its reduction limits, decrease c and go to step number 5. If the step size can no longer be reduced, go to LOP-D8 for final de- termination of the optimum. ' SUBALGORITHMS h——CURV1TURE DETERMINATION NOTE: For either of the following, the value finally see lected for curvature may be given.by Equations 3.120, 3.121, 3o123 or 3o12’4o 76 Select E1; , Ae£ggtagiyi .t ‘ fi‘ 1 = 0 (:START TRY AGAIN? w indicates: initial condi- tions p3 and x0; integra- ubalgorit. tion of Equations.2.l9 and H:Curvatur 2.28, utilizing a maximal cona terminat. trol, to yield 1 and.p}; then . ‘ evaluation of E( 3). - . @ indicates. ° initial conditions p ; integrat¥on.of Equations 2.19 and 2. 28 in reverse time to give and pure xd :- pubalgcrit fiEvalugfet and ,0 Digit GO TO ' LOP-BS {iii} 3 < FIGURE ‘ Flow Chart for LOPuCMI 77 1 resulting in 5:. o 11,. Determine is}. and it}. from x0, 13.}. and extremal switching. Lha. I . Perturb p III. If i}. a 1.}, increase the perturbation size and go to step I. If not, go'to step IV. IV. Calculate an estimate of I: from i}, 13.}, i 1.}. and pr. h.b.1. 'Perturb p: resulting in 39:! J :- 1. II. Determine in}. and 3p; from x0, 1p; and extremal switching. III. If 31.}. a in}, increase the perturlntion size and go to ste'p II. If not, go to step IV. IV. Determine the standard deviation jan' ”(B/Fa. where jam) is the standard deviation of the 51's and Ea is the average of their absolute values. If on is less-‘than ”a (allowable limitfor the standard deviation), then calculate an estimate of k. If Jon is greater than c then go to step V. V. Perturb pg in a random manner from the pre- vious perturbations, yielding 3+lpg. Unless 3+1 is greater“- than a preset limit on the-‘number of. allowed initial costate perturbations, go to step II and repeat steps II through V. In the event that 1+1 exceeds the preset limit, use-the value of p: for which o'n is a minimum to estimate k. SUBLLGORITHKS B—NEI INITIAL COSTLTE IETEBMINATION 8.a. Utilize Equation 3.127, i.e. p3“ - 3:. do not correct for 61:. 78 Select pz, A and”? Given Po START For meaning of the symbol refer to Fig-.130 3o5o i . I ; N11" jNiTul) i IiT" 3x11. 1 Generate a random ve ‘ tor p —' i=1+1 0 T0 MAI PROGm m Calcul # [ 10112?) J 0 TO m1 * PROGRAM 1 FIGURE * 6 F ow Chart for Subal orithm it b 79 8.11.1.1 Determine 21 and p1 corresponding to p; and IT (without perturbations) and extreml switching in reverse time integration. ,II. Calculate the correction factors of Equations 3.128 and 3.129. 111. Correct 3: 3.131. IV. Let pz+1 ' figo Ignore 61:0 and 2: using Equations 3.130 and -i -l 8.0. Using either x0 and p0 or i1 and p1(as determined in Subalgorithmo 8 .b),0 define a new opifl by using Equation 3.133. Examination of the above algorithm and subalgorithms points out a number of possible alternatives within the basic algorithm. By way of sumry these alternatives occur at: 1) Step 1. -- Error Function Selection 2) Step 11'. - Curvature Determination 3) Step 5. - Final State Perturbation 14) Step 8. -- New pg“ Estimtion. In addition to these, there are alternative step sizes to be chosen and also step 6 could be altered so that m}. would not be perturbed, or some of the other unrelated perturbations as discussed earlier (See Figure 3.2) could be selected. In Chapter 5 comparisons are made using these various a1 ternat ives . C H AP T E R ’4 THE GLOBAL OPTIMUM AND RELATED PROBLEIB In the previous chapter, a method (including several options) was developed for determining the local optimum of a convex subset of R(T) (i.e. solving LP). Direct applica- tion of this method to an arbitrary R(T) might encounter one of several special cases for which the results of the previous chapter need to be reassessed. In this chapter, means of identifying and treating these are discussed. lost of these special cases are actually part of the more gener- al problem of locating the global optimum; hence it is not necessary to implement specialized techniques for their solution. The more general problem of determining the global optimum for an arbitrary R(T) which my. have. several or even numerous local optima is most difficult. In this the- sis a random approach is taken for finding the local optima. and thus identifying a global optimum. This global optimum procedure solves Problem 2.14 (NP). Using some of the con-3 oepts of Fadden and Barr [F1 and BI], other types or opti- mal control problems can also be solved. Of these,.on1y the time optimal control problem is. considered here. All experimental results and example problems are presented in the next chapter. 80 81 14,1 Special Pgoblem l: Nonconzex regions on Q3”) The algorithm presented in Chapter 3 assumes that the region on 3R(T) boing considered is convex. For an arbi- trary reachable set there may be mnylocal optim (i.e. locally closest points to the origin). Much of the surface my not be convex and, in fact, some of these local optimal and even the global optimum my lie on a non-convex region. Such regions my be concave or my be “mixed" (saddle points, neither concave nor convex). Upon first consideration it my seem that it is not possible for an cptimmn to noon a concave region of DIRT). This, however, is not .the case. Consider, for example, the reachable set shown in Figure 11.1a, with the origin located as indicated. Since all of the surface (in this case a ' curve--R(T) is 2-dimensiona'l) near the origin is concave, the Optimum is obviously on a concave re'gion, namely at 1;. A very special case is shown. in Figure 15179, where the op- timum is not unique: infact, where there are an infinite number of minim, all equally close .to the origin. As a final example, Figure 14.10 is given in which the optimum _‘ is located at a "corner". .Note that the determining fac- tor in all three examples is' the‘relationship between the radius of curvature of the concave surface and the, distance of that surface from the origin. This observation leads to the following theorem which is given afterlsome preliminary definitions. 0 .1 82 = (T) ' o , origin origin origin . # R(T) R(T) 1% a. f be I. _ co FIGURE 14,1 gptim on Concave Boimg_a_r_z Regiw , I . , DEFINITION “,1 Let DIRT), x(T) and I: be aspreviously defined. Let' theprincipal curvatures 1:1, 1 - 1,0vf,n-l, be nonzero, then the principal. radii g_f_ curvature, fi’ at x(T) are defined as 1 f1 B-IE? o (uol) IEFINITION h 2 Let EMT), x(T) and It. be as previous- 1y defined, then the mximum convex c ture, Kn, at x(T) is defined as , ‘ x” . min ki , if any k1< o, (14.2) or . . Kn = o, if all 1:1 2 o. (11.3) mrml'rlon 14,; Let em), x(T) and 1111 be as previous- 2 ly defined, then the mximum concave curvature, K0,, at x(T) is defined as 83 I" a mix k1, if any 1:1 > O, (14.1!) or Kev e o, if all 1:1 5 0. (M5) THEORE]! 11,1 .112! R(T) be a compact set in ER, on R(T). Consider a point x(T) e 3R(T) for which. cos y = -l, where p(T) is defined as the outward normal to R(T) at x(T). .1451 No(x,c) be a neighborhood of x(T) on 3R(T). _1._e_t_ k1, ;' i a 1,0",n-1 be the principal curvatures at x(T). a. _I_£ . m > {—1—— (M6) . . ’ 111‘“ then x(T) is not a local optimum: i.e., there exists a y c Nr)(x,e) such that “7” < III”) II- (in?) b. l; . .1 _ ‘ - K07 B m, (14.3) then x(T) is not a unique optimum; i.e. there exist many 2 c No(x,e) such that IIZII '3 ”x(T) ll- (M9) c. Finally, if _ 1 . _ K < 04.10) W llx'rll ’ then x(T) is a local optimum: i.e. for any w c Na(x,c), .u an >. n and) H. (14.11) Proof: Consider part a. first. There exists at least one principal curvature, namely RC? 3 Kbv (“.12) for which 1 . k" ) W o (“313) Consider the line of curvature, 1.", through x(T) cor- responding to kcv' Since ch is a line of curvature, Equa- tion 3.66 applies. Since cos 7 at x(T) is -1, Equation , 3.66 reduces to: d 008 x = -(kcv _ A x¥r ') d “IdssTl II . (1.1.114) By Equation (4.13 this reduces. to H31? =--(“fin—(fa—Lflj‘” . (“-15) where l is a positive number. If an infinitesiml change is made in x(T) along I‘cv’ i-CT’i‘iE-l > O (14.16) since cos y is a minimum at x(T). Therefore by Equation b.15, d 451.1 < 0. (£1.17) Let y be a point on ch infinitesimlly close to x(T); thus by Equation l$.17, Equation ’4.7 is true and part a. is proven. Now consider part b. If for one curvature, k a K 8 A 1 ' “.18 then at x(T), d—gm a: O, (14.19) ds 85 over a segment of Lev’ hence the optimum is not unique. Finally consider part c. In this case, the number A of Inequality 14.15 is negative, regardless of the choice of principal curvature. Since the principal curvatures repre- sent the extreme values which curvature can take, K rep- cv resents the maximum value for curvature on the surface in the neighborhood of x(T). Since Equation 14.16 holds and l is negative, ‘1 'IsT > 0 (11.20) and x(T) is the local minimum. Hence the proof of the theo- rem. It is interesting. to note that the important relation- ship between k and 1 /|| x(T) [I for concave surfaces is con- firmed by Equation 3.72 for the composite error flmctien E5, In this case the requirement on k for monotonicity of the error flmction is the same as that for the existence of a local optimum on a concave surface (part c of Theorem 14.1) . Since the object of this chapter is to eventially adapt LOP-Cl! to the global problem, consider the results of Theorem “.1 as they apply to LOP-CM. First of all, it is evident that any minimum of cos 7 must be verified if even one of the principal curvatures is positive. In this event, one of the alternatives suggested by Theorem 14.1 applies. Thus, consider the following decision scheme for computa- tion based on these alternatives: " 86 1) If Equation 11.6 applies, the apparent 'cptimum must be a‘ ”false" optimum. Switching to E (p ) = ”x(T) [I effectively overcomes this difficulty. Note that this circumstance is improbable since the nature of LOP-CM, par- ticularly the method of choosingox, contra- dicts Equation 11.7. The chance of an arbi- trary trajector, x0(o) yielding this special case is likewise remote. ‘ 2) If Equation 14.8 applies, the optimum is not unique, bu computationally x(T) is one of the (infinite number of) local optima. 3) If Equation ’4.10 applies, then no computational change need be .mde since x(T) is the local Optimum. It is appropriate to mks an observation at this point. Considering the fact that a convex reachable set results from a linear system and considering the very nature of the optimization problem with O 6 R(T), it is reasonable to ex- pect that the global optimum usually lies on a convex, or at worst, a mixed (saddle-point) surface. Thus if an init- tial guess for po yields a boundary point where Equation 14.6 applies, a mJor change or jump in the initial costate might be most effective rather than a routine application of LOP-CM. In the discussion given above and in Chapter 3 only convex or concave surfaces have been cons idered. The re- sults, however, apply to mixed surfaces which are neither concave nor convex. For example, if Equation 15.10 holds (which allows the possibility of negative principal cur- vatures) then it is quite possible that the optimum lies on a mixed region. For such a region, examination of 87 Equations 3.66 and 3.72 indicates that if a choice of path is available, based on the value of curvature, then the principal curve corresponding to the minimum value for k should be chosen. Since the actual computational method of perturbation from iteration to iteration is based on -the ox's, the value of k only enters in the computation of 6p. 14,2 Special Pgoblem 2: Flats and Corners g. In the application of mny reachable set-oriented com— putational methods, two surface characteristics are partic-g ularly troublesome. These are the £13.55, for which one out- ward normal‘ corresponds to more than one adjacent boundary point, and the corner, for which a unique norml plane is not defined to the surface at that point. These are illus- trated in Figure 11.2. The norml plane to a flat is said to be nonregular and the corner is termed a nonregular ; point. p c(pT at corner) pf (pT at flat) -' p c .0 a a. Flat only b. Corner only 0. Flat and Corners FIGURE 11,2 Flats and Corners on Surfaces 88 Consider first the flat. In mny computational meth- ods which determine contact points [132] corresponding to a specific p(T), sets such as that shown in Figure h.2a and 14.20 are especially difficult to handle since mny boundary points correspond to the outward normal pf. In Section 2.7 it was indicated that the flats result from singular con- trols in linear systems. For nonlinear systems, however, flats my occur even though the problem is nonsingular. If the flat is not the result of singular controls, then it does not present special problems for the computational method introduced in Section 3.7. Although 'the direction of p(T) is constant for a region on 3R(T) (i.e. k is zero), there is a change in the final state corresponding to a change in the initial adJo int, hence any of the previously introduced error functions differentiate between boundary points even though the outward normls are the same. If the flat results from a singular control, then gr_1_e_ initial adjoint corresponds to all of the points on the flat. Thus if only the mximl values of the control are assigned at singular points, only the boundary points of the flat are determined, resulting in the singularity gap mentioned in Section 2.7. This gap is only significant if the global optimum is on the flat. In the event that this optimum is the only optimum, then the failure of LOP-Cl! to converge would identify this special case. If, on the other hand, other optim exist, then a global technique 89 would incorrectly identify one of these optim as the global optimum. The possibility of such an occurrence is slight but must be considered. Short of restricting LOP-CM and its global modificaa. tion only to norml problem, there is an additional alter- native which can be taken. Any minim (but not optima) located by a global procedure would have to be carefully examined to verify that they are not boundary points of a flat, and that. lle" for these points is not, in fact, less than “IT" for the supposedly global optimum. While this approach does not guarantee the identification of a global Optimum on a flat, it does reduce the already small possibility that such an optimum is overlooked. One final observation should be made. For linear systems the adjoint equation is fixed regardless of the state trajectory; thus to each p0 there corresponds one pT even though the state trajectories my vary (due to singularities). This, of course, results in a flat on 3R(T). For nonlinear systems, however, the adjoint equa- tion is dependent on the state trajectory, thus an initial adjoint which results in singular controls usually would produce differing final adjoints since the state trajec- tories vary. While no conclusion can be drawn, this would seem to indicate that flats caused by singular controls are less likely for nonlinear systems than for linear systems . 90 The second special case to consider is the corner. 'hereas the flat spot results from a boundary section for which p(T) is constant in direction, a corner results from a boundary point which is constant for a number of final adjoints (hence initial adjoints). In this event no par- ticular computational difficulty is introduced for LOP-CM if any of the error functions including the cos 7 factor are used. Since the final adjoints vary, the error func- tion oranges even though the final state does not clmnge. Note that the effective curvature in this case is infinite, thus indicating a perturbation in p(T) but no corresponding perturbation in x(T). 14.3 Special Problem 3: Ememum but Not hocal Optimum If error function E1 = ”x(T) [I is used, any minimum of the error function must also indicate a local optimum. For the other error functions, however, a minimum my be ob— tained but my not be the optimum value. For instance, cos y may reach a minimum but not the optimum value of -1. This possibility? is illustrated in FigurezbJ. ’r‘lhile LOP-IE or LOP-CM may converge to, such .a point, no particular com- putational problem would result since any global procedure excludes the choice of such an extremum as the global opti- mum. In the event that convergence to such an extremum is to be prevented, E1(p°) is utilized. 91 origin R(T) x(T) = extremum of E1 but not a ' local optimum; in fact cosyh is greater than -1. FIGURE h,3 Extremum.but Net a Local Optimum for Cos 1 b,h §pecia1 Problem.h: Internal Boundaries Consider the complement of R(T): Q(T) e [ x e En : x e R(T) J. (n.21) Usually Q(T) is connected; however, a situation as shown in Figure h.hc is certainly possible. Designate the por- tion of G(T) which contains the origin as QO(T) and index the remaining disjoint subsets of G(T) as Q1(T), i=l,oo-,k. IEFINITION Let) The 332 internal bounda , diner), of R(T) is the subset of 9301') defined by oinm = menses). (14.22) Several observations can.be made concerning Q1 of Figure h.#c. Because of Theorem 2.5, it is known that any x c J&R(T) was a boundary point for all reachable sets R(tL t < T: thus one can.assume that tho reachable set shown in Figure h.hc evolved in a manner similar to that illustrated in.Figures h.ha and h.hb. 92 R(tl) R(tz) o O origin origin origin FIGURE '4 14 Develo ment of an Internal Boundar The global optimum (or optim, if not unique) must lie on 30R(T) : W L_e_t_ R(T) be a compact set for which Q(T) is not connected. 1.21 Q0(T) be the connected subset of Q(T) which contains the origin. 1119}; any global optimum belongs to 0011”). Proof: Let x* c 3301') be a global optimum. Assume that x“ e 3312(1)), 3 .L 0. Thus x“ c 6301-), :j ,4 0. But the ori- gin belongs to QO(T), thus there is no path lying entirely in 153”) from x* to the origin. Consider the line from x* to the origin. Since for any i i j, Q1(T) and Q3(T) are disjoint, part of the line must lie in R(T). But this con- tradicts the assumption that x* is a global optimum, hence x* e ()JMT), j :4 O. The only portion of R(T) remining is JOINT) and since x‘ 63R”), x* c 30R“); hence the proof. 93 Note that the theorem is proven for a global optimum, x*. A local optimum may lie on.an internal boundary (See Figure “.5) and thus may be determined by LOP-CM. (Although a means of identifying such an.eptimum is available (i.e., examination.of the line segment fromx+ to the origin), such an.approach is not necessary since a global technique cannot select such local optima as global optima. Thus, no specialized algorithm is given to differentiate these local optima from other local Optima. “.5 Special Ppoblem 5: [also Boundagz;Points As has previously been indicated, all extremal coup trols are maximal controls but the converse does not hold. As a result, ”false boundary points" (fbp) are generated. 0 origin FIGURE “.5 A Local Optimum on.an Internal Boundar: 9“ DEFINITION “.5 1 false bound_a_1_‘_z. point is a mximl endpoint which is not extremal (i.e., x e OR(T)). It can be seen that the situation depicted in.Figure “.“ can result in false boundary points. The evolution.of the reachable set is such that the points lying on the line F(T) in Figure “.“c are false boundary points. Since LOP-CM utilizes maximal controls with.the an- ticipated result that x(T) e DENT), false boundary points could cause computational difficulties. If LOP-CM were re- quired to converge only to true (i.e., not fbp's) local optima, some specialized technique would be needed to idenp tify and treat the case of false boundary points. Such a method could be based on the following facts: 1) Each fbp (e.g., x e F(T) in.Figure “.“c) previous- ly existed as a boundary point for some t1 < T. 2) At some time t1, for each x1 e F(T), the state— costate pairs (x1(t1),p1(t1)) and (x1(t1),-Pi(t1)) represent true boundary points. 3) For each x e F(T) there’exists a y,' y =- cx, O = -1. (“.2“) a +1, (“.25) mg x“(T) is not the optimum for MP app 0 e R(T). Proof: COnsider the first part of the theorem. Since x*(T) is the global boundary optimum, it must only be shown that o d R(T) to establish that x*(T) is the global optimum. Since Equation “.2“ holds, x*(T) and p*(T) must be collinear but oppositely directed. Let L(x*(T),O) denote the line segment from x"'(T) to the origin. Since the final adjoint, p*(T), is directed outward from R(T) (i.e., from x*(T) ) along L(x*(T),O), points on L near x*(T) must lie exterior 97 to R(T). If this line were to intersect R(T) then it would include at least one point on BR(T), closer to the origin than x*(T). Since this contradicts Equation “.23, the ori- gin.must be external to R(T) and thus the global boundary optimum.x*(T) is the optimum for HP, Now consider the second part of the theorem. If Equa- tion.“.25 holds true, then xf(T) and.p*(T) must be collinear and similarly directed. Since p‘(T) is outward to the . reachable set, points on the line L(x*(T),O) near x*(T) must lie in R(T). Thus the origin must lie internal to R(T) or the line L(xf(T),O) would intersect DR(T). If the latter were the case, then.Equation “.23 would be contra- dicted, thus the origin must belong to the reachable set. Hence there is a trajectory endpoint (namely x(T) a O) which is closer to the origin than 1*(T); hence x‘(T) is not the optimum for MP. “.7 The Global Optimum In the previous sections, special problems have been suggested as possible difficulties in determining an.epti- mum. Careful consideration, however, has demonstrated that these special problems are really part of the global prob— lem. In as much as an explicit equation is not available for.aR(T), the global problem is very difficult, especially for high order, nonlinear systems. Several possible approaches can be taken to extend LOP-CM to determination of global optima. Once a local 98 optimum, x+, has been located, then the global optimum must be located within a hypersphere of radius [I x+|I. One meth- od of determining a global optimum would thus be to inves- tigate all possible final states within this radius. Since extremal endpoints, and thus global optima, are associated with initial adjoints (assuming normality), any other approach is to consider the set of all possible ini- tial adjoints. Since only the magnitude of the adjoint (hence initial adjoint-~see Equation 3.37) is important: max=max, O. 0. Consider the alternatives given.above. As was the case for GOP, there is a tradeoff between the reliability of the method and the amount of allowable computational time. If GOP is utilized at each time increment and if the time increments are made sufficiently small, then.there is confidence in the choice of t‘. On the other hand, to be realistic, some compromise in computations mmst be made. For example, GOP may just be utilized at the initial itera- tion and as a check of the final result. In.as much as the step size alternatives (step 3) are not difficult to implement, the choice should'be‘based entirely on.the effectiveness of the two alternatives. For the alternative given in.3), it is logical to use prior information, rather than starting with an arbitrary initial 102 costate. In fact, if GT1 is kept reasonably small, the new optimum initial costate “11:; my be very close to 1p;, in— dicating that LOP-IE is a good candidate for use in alter- native h). Caution, however, must be exercised in attemting to reduce the computation time. Even though a T:l has been determined for which 0 o DMTJ) and even though 0 o R(Ti-l), it is possible that there existed a Tk < T3 for which 0 e 3R(Tk)‘.. If the step size were too large or if an in- correct, local optimum was considered the global optimum at each time increment, then this situation could occur. One final observation should be made. The successive values of final states at each iteration, 1x(Ti) do not necessarily approach the origin monotonically. It is quite possible that an optimum final state 3x(TJ) has the prop- erty that . . . . u Jx(ri) u > u J-lx(r1-1) u, (n.29) even though the sequence of .1x(T1)'s converges to 0. Examples of this as well as other aspects of TOP are given in the following chapter. CHLPTER,5 COMPUTATIONAL RESULTS AND CONCLUSIONS In.the preceding chapters, general problems, methods and alternative choices have been given which can.effec—‘ tively be evaluated utilizing computational results. It is the purpose of this chapter to provide these results and to form.conclusions based on the data obtained. Specifically, example systems are detailed and dis- cussed in the first section. These systems are then used as examples in.the remaining sections of the chapter. These computational examples are introduced with several purposes in mind: 1) To give insight into the nature of the problems. 2) To present the resulting reachable sets and give example extremal trajectories. 3) To use these examples in comparing the various computational alternatives previously introduced. b) To demonstrate that the previously introduced algorithms are effective in.computing Optima. In particular, the author feels, based on personal experi- ence, that example nonlinear problems with resulting reach, able sets should be given for the benefit of others wishing to consider this type of problem. As the computational results are presented and analy- zed, comparisons are made and conclusions are drawn. These conclusions are listed in the various sections of this chapter and are summarized in the final section. 103 10b .1 Exa le S stems, In this section several example systems are intro- duced. For each of these systems, a differential equa— tion is given, state and costate equations are introduced and maximal control is considered. For simplicity, these results are presented in outline‘form. EXAMPLE SYSTEM £1 (ES-1) a. Differential Equation (13$th .01 [x(t) [x(t)] 1+ .1[%Q§u(t)] 55.1) b. State Equations ( x = x1 i1 (5 2) a u .» o 5:2 -.1| x1|.1x23 c. Hamiltonian H = .lplxz- .lp2[ x1 I x1 + .lpzxzu _ (5.3) d. Ldjoint System N ' a (50") pz -01 ' “.1111 e. Control 1. Restraint Set: [u] 5 ii. Maximal Switching: u a sgn (.lxépz) a sgn (xzpz) (5.5) 1’. An analog diagram for this system is given in Appendixwl. 105 EXAMPLE SYSTEM £2 (ES-2) a. Differential Equation: 3 2 du dx=_ dx de __?. E3 153+(at) +“14"”‘2J’ dt (5’6) This system equation was derived from a well- known equation of fluid dynamics known as Schlichting's Equation (Equation 5.6 becomes Schlichting's Equation if u a l and =0). In general, equations of this type r resent fluid flow across surfaces. With the addition of u , flow across a wedge is indicated. [13, 1]. b. State Equations (x1 2 x) 121’ *0 1. o— z: 70 er i2 8 0 0 1 12 + O l “1 (507) u 0 2 :34 f3 12 (3. IL 1 o c. Hamiltonian 2 . ' H 3 13112 + p213 + p2u2 - DBIBII + p312 + p3ul (5o8) d. Adjoint System '2“[ 1" "' ”"1 pl 0 0 x3 1’1 . fl, = -1 0 -2x2 p2 (5.9) p 0 -d x P L34 L. ' 1.... L3. e. Control 1. Restraint Set: [nil 5 1, i=l,2 ii. Maximal Control: [11] 3 [p2] (Solo) 1’3 106 EXAMPLE SYSTEI, ES- a e a. Differential Equation (for the system shown.in Figure 5.1) d3! 3 -«QE§ - 12,2; + u(t) (5.11) dt3 dt2 ' at ' ‘ b. State Equations (xwx1 ) F. T — —- gl— .1 I-O-q x1 0 l 0 x1 :12 s o 02 1 x2 + o u (5.12) i O -x -J l c. Hamiltonian 2 H = Izpi + are - 1112p: - 131°: + up: (W) d. Adjoint System ‘ h. 1' P _ _ _ pl 0 0 2"1“2 I’1 p O -l l P .L. 3. L ._ i. :1. 9. Control _ i. Restraint Set: In} f 1 ii. Maximal control: 11 = sgn (p3) (5.15) u(t) .1 + s2 (s + 1) ’1”) dx(t earity FIGURE 1 Exa 1e stem. 2- rd 0 der Nonlinear 1 107 EXAMPLE SYSTEM f4 (ES-J4) IE}, page 226] a. Differential Equation U 3 d x t g _2 d x t _ t t ( .16) —§-)-dt -——§—)-dt s[x( )] + u( ) 5 2 2-~ l h S[X(t)] 3 "5’1“” “1'3 x(t) *5-5-6 x(t) (5.17) b. State Equations: (x exl ) E17 — o 1 o 6- E1— 0‘ ’ a o o o 1 + (5.18) T3 t 21:1 11 x3 “2 x ~+F‘m O O -2 xu “l _L.L ' J J- J. _. .4 c. Hamil—tonian 2 2 l '4 H = plx2 + pr3 «I-pr,4 + fpuxl +1-5pux1 - W311 - qupg, + 133,111 + p3u2 (5.19) d. __A_d:loint System .1 131 {—0 o 0 .. fi- «in—5- + g; Fill— 152 g -l o o 0 p2) (5.20) 133 o -l o 0 pg 31“} _o o -1 2 .1 ”pu-J e. Control , . i. Restraint Set: lull s 1, [112‘ S 1, ii. Maximal Control: u2 sgn p3 108 5.2 An Introduction to Computational Examples - The computations summarized in the following sections were performed in the Hybrid Simulation.and Control Labora- tory at Michigan State University. In this lab an IBM+1800 digital computer and an.ADbh analog computer are linked to- gether, thus allowing digital control of analog operation and also hybrid computation. In all computational examples in this thesis, the dig- ital computer was used to provide overall control (direct the optimization routines, implement alternative compari- sons, etc.) and to input/output data. Depending on the example, either the analog or the digital computer was used to integrate the state and costate equations at each itera— . tion. While the capability exists, the digital was not used in oneline integration linkage (i.e., true hybrid cp- eration) with the analog computer. The digital integrations were performed using a hth order Bunge-Kutta method. While this method is relatively slow on the IBM-1800, it is reasonably accurate. It is presumed that the most significant source of inaccuracies is introduced through.the control, which in.all example systems is only piecewise continuous (signum.funotion COD! trol--signum switching). To partially overcome this dif- ficulty, whenever a discontinuity in any component of the control occurs, the integration step size is reduced by a factor of ten for the interval containing the control dis- continuity. 109 The analog computer, on the other hand, is generally faster for equivalent accuracy. Its accuracy, however, cannot readily be improved by reducing step size, as is the case for the digital.computer. It also has the tendency to be less consistent. This is especially apparent in the forwardpreverse time integrations as is indicated in Sec- tion 3.6.2. Example computer programs and subroutines are given in.tppendix B. It should be noted that these programs are more complex than.the basic algorithms because capabilities (data and sense switch options) for quick alternative com- parisons are included. host of the alternative method comparisons were made using ES-l and.ES-2 for various reachable sets. This is possible since the nature of the reachable set significant- ly changes, for a given.system, with changes in initial state, initial time and final time. For ES—l the analog- digital computer combination was utilized. Otherwise, to- tal digital methods were employed. The comparisons are given in.two sections. The first contains comparisons for the various algorithm alternatives. The second presents example optimization.problems and in, cludes example trajectories, reachable sets and comparisons with a totally direct search method. As previously menp tioned, conclusions are included within these sections and are clearly identified. 110 5,} Computational Copparisons of Algorithm Alternatives In Chapter 3, algorithms have been.presented which include alternative choices and subalgorithms. Although theoretical considerations indicate, in.most instances, which alternatives are best, it is the purpose of this section to substantiate (or refute) these choices on the basis of actual computational results. For each comparison, only one change or alternative is evaluated, Taking a specified group of reachable sets, LOP-CM is started at an arbitrary (but fixed) po and the number of iterations required for convergence is measured for the first of the two alternatives. Then the second alternative is incorporated into LOP-CM and again the spec- ified reachable sets and.po's are used. A comparison of the average number of iterations required for the two al- ternatives or of the average percent decrease in the error function.per iteration for each method gives an estimate of their relative value. All results are identified by the alternatives to be compared and the basis on which the comparisons are made. These results are presented in.semi-outline form and are arranged in.an order corresponding to the discussion of the alternatives in Chapter 3. It is only possible to present a summary of the computational work and example results. In.most cases, data, other than that listed, also confirm the results presented. 111 5.3.1 Stud: of Perturbation Rglationships--Fina1 Time In this section, comparisons are made of some of the possibilities introduced in Section 3.6.1 (See Figure 3.2). Specifically, the relative effectiveness of a joint per- turbation (both IT and.pT) versus a single perturbation.(xT and PT) is considered. a. Joint Perturbation.(;d and_pT via Eguation.2.28 1g 3,1. perturbation only. i. Copparison.Basis: ES-l, Analog Integration, t0 = O, T = 20 secs.,lnput Data Set 1 (See Appendix C). ii. Illustrative Exa 1e: In Figure 5.2 the reachable set for x = 13,-35T is given. On this set the sequence of final agates generated for each approach of this comparison are given (joint-~13 , xT only- 0). iii. Summa ‘2; Co arison Results: The average number of iterat one requires io aitain a Iocal optimum were compared for the joint perturbation and for the per- turbation of xT only. The joint perturbation.approach re- quired an average of h.7 iterations whereas perturbing xT only resulted in an average of 6.0 iterations. iv. CONCLUSION: The joint perturbation.approach is more effective than an xT only perturbation approach. b. Joint Perturbation _Y_S_ p,r perturbation only. i. Copparison.Basis: ES-l, Analog Integration, t0 = O, T = 20 secs., Input Data Set 1: ES—2, Digital ID! tegration, to = O, T = .5 sec., Input Data Set 3. ii. Summa ‘2; Copparison.Resu1ts: The average number of iterat ons require to attain a Iocal optimum for the joint perturbation.was 8.1, for perturbation.of pT only, 8.5. Another means of comparison is the average percent improvement for each curvature move (perturbation of xT and PT). On this basis, a joint perturbation yielded an aver- age 56% improvement while the pT only approach gave an average 51% improvement. iii. CONCLUSION: The joint perturbation approach is slightly better than a PT only approach. 112 x2 1 l 9 R(ZO) for i . T r _ 3 :0: (59-5) 7 origin 6 . g "' .9... ‘ =2. 2...... b .---_-2 6 §r--- 5 ‘- R(20) for T 10: (59‘1) , Dindicates 3 int perturbat tons Oindicatos onl 5. perturba lion 1 numbered .-L \N ‘— -7‘ill 2 x1 FIGURE 5.2 Joint vs Sipgle (IT only) R(T) Perturbations 113 0. GENERAL CONCLUSIONS: On the basis of the two come parisons made above, it can be concluded that the perturba- tion of the final adjoint has the most significant effect in improving the error function. The improvement due to perturbations in the final state are probably diminished due to the fact that slightly inaccurate perturbations in xT result in perturbed.points lying off the boundary of the reachable set. In either case the joint perturbation ap- proach is the most effective. 5.2,2 Copparison of Curvature Values In.Chapter 3 several alternative means of estimating curvature are given. It is the purpose of this section to compare these means and the resulting curvature values. In as much as the estimates of the curvature are obtained by perturbing the final state-final costate pair, one can ex- pect that small errors in the determination of the final state and final costate can introduce significant errors in curvature determination. Because of this, the analog and digital means of integration are analyzed separately since significant inconsistencies do result from analog computation. a. Copparison_ of Curvature Values as Determined on the Analo Copputer (using Equations 3.120, 3.121 and 3. 123%. 1. Co arison Basis: ES-1,tO:-O, T==20 secs., Input Data Set 1 ESee Appendix C). ii. IllustrativeE mppl : In Figure( 5. 3 the reachable set correspondingE to T= 20 andx 5,0)T is given. Curvature values for several extregal endpoints in an iteration sequence are given. 11h x2 x1 0 1 2 j 3 L1 _2_ k s: :37 (convex) 2 k s —,.50 (convex) it. “ . Ecomx) ! f _i icates E. ittration number. I R(20) for -5 -6 / -7 l/k s .0036 (concave) A A / FIGURE Exa 1e Curvature Values in an Iteration Sg'uence 115 iii. Summarz of Copparison Results: In.as much as ES—l is a second order_ system, the boundary of the reachp able set is a line of curvature. If the perturbations of and PT are sufficiently small and if no computational inaccuracies are encountered, all estimates of curvature as defined in Section 3. 6. 1 should coincide. This is not, however,the case for higher order systems as the perturbed final state may not be on a line of curvature through the original final state. Comparison results are given in Table 5.1. There is reasonable correspondence between'k and except for in_ put data numbers 2, 6,10 and 11. Since gel is a second order system, these inaccuracies must result from too large of perturbations or from integration.inaccuracies in the analog compuper. It is interesting to note that in.these four cases, k1 and'kz have Opposite signs, indicating that the estimate of k is not accurate. Several methods of cor— relating the various curvature estimates are available. One is the standard deviation of the'k 's (denoted u(k)), another the standard deviation of k an (denoted u(k, r? ) and a third, the ratio of u(k) and the average 1) (denoted an). The importance of these various comparisons are discussed further in comparison 5. 3. 2. d of this section. iv. CONCLUSIONS: Experimental curvature esti- mates agree with those expected from reachable set geometry. Perhaps the best measure of the accuracy of the curvature determination (also the measure of whether or not the per- turbation is on.a line of curvature for greater than sec- ond order systems) is o = a(k)/ . If this value is too large, then.the perihrbation [(ioi‘v determining curvature) is too large, computational inaccuracies have resulted , and/or the perturbation is not along a line of curvature. b. Copparison of Curvature Values as Determined on the Di ita Gogputer‘l’using Equations 5.120. 3.121 .. and 3.1 ) for a Second Order sttem. 1. Co arison Basis: ES-l, to -O, T= 20 sec., Input Data Set i (See Appendix C). ii. Illustrative Exapple: In Figure 5. h, the reachable set corresponding to x0 = (- lO ,-5)T is given. Curvature values for several extremal endpoints are in- dicated on the boundary. iii. Summarz of Copparison Results: In Table 5.2 estimates of curvature— are given. It is readily apparent that there is good agreement between these estimates. 116 TABLE 5.1 Copparison of Various Analog Estimates of Curvature INPUT DATA f* 1 \OmVO‘xUl-F'N 10 11 12 13 10 15 3.1.. ‘2 .5:— E .13.)- £55.!) .31.. .3630 .1217 ..2023 .2175 .1206 £0120 .0977 .0809 —.1111 -.0157 1.g311 .0966 .6230 1.0000 -.0185 —.0159 -.0172 .0123 .0013 10108 .0756 .1037 .1302 .1389 .1390 .0005 .0001 .0338 —.2101 .0870 -.0616 2.2019 .1086 1.1518 1.0000 .0799 .0975 .0887 .0821 .0088 .0033 .0992 -.0256 —.0552 —.o000 -.1670 .0108 10633 .3663 .0563 .0520 .0502 .0055 .0022 .0000 .0006 _45915 -.6657 -.0371 1.0009 .6286 .5010 1.0000 -.0008 polg1 ».0006 .1208 .0095 .0607 1.0000 .2186 .1980 .2085 .2096 .0101 .0005 .0089 .0650 .0538 .0590 .0921 .0056 .0163 L0903 -.0168 -.0159 -.0160 .0358 .0005 .0261 .0305 .0160 .0120 .0100 .0090 .0020 .0158 .1389 * determined for Data Set 1 (See Appendix C). No data was obtained for #3 since the initial guess of p '. was approximately the optimum; hence no curvaturg estimates were calculated. 117 J‘2 60 -.'01F "’ .16 50 00 30 on —"VVJ 20 lO 0 10 20 -1. -M -0O FIGURE 5.0.Exagple Curvature Values for R(20), xn= (-—10,...5)T 118 TABLE 5.2 Digital Estimates of Curvature for a 2nd-Order System INPUT DATA f* _El_ .1500 -.7205 .0002 ~1559 .0195 .0685 .0868 -.0292 -.2000 .0697 .0003 -.0019 ‘iev .__EL__ .JiSEl_‘9§EiEaI).__SR__ ._:1539 .1501 .0005 .0001 .0033 -.7166 -.7117 .0039 .0025 .0055 .0002 .0001 .0000 .0020 .0020 .1561 .1557 .0002 .0003 .0010 .0195 .0210 .0000 .0010 .0009 .0686 .0670 .0001 .0008 .0016 .0866 .0806 .0001 .0010 .0018 -.0292 -.0301 .0001 .0005 .0020 -.2027 -.2062 .0013 .0018 .0055 .0696 .0690 .0001 .0003 .0013 .0003 .0086 .0000 .0022 .0009 -.0019' -.0023 .0000 .0002 .0051 I"determined for Data Set 1 (See Appendix C). 6 and 13 are left off since the original guess was so near the optimum that LOP-BS was employed rather than LOP-CM; hence no estimates of curvature were calculated. Data numbers 3, 119 iv. CONCLUSIONS: For the second order system be- ing considered, perturbations of the digitally integrated differential equations gave very good estimates of the cur- vature. As expected, all variances were low: thus verify- ing that the perturbation size was sufficiently small to give a good estimate of curvature. c. Copparison of Curvature Values as Determined on the Di ital Co uter—(using Equations 3.120, 3.121, SoEZS; f 51 and or a T rd Order System. i. Copparison Basis: ES-2, t°=O, T: .5 sec., Input Data Sets 3 and . ii. Summar ‘21 Copparison Results: In Table 5.3 various estimates 0 curvature are given. Unlike the pre- vious comparison, made for a second order system, most data points have rather large variance of curvature values. This is, of course, expected where the boundary of the reachable set is not implicitly a line of curvature. Er- ror function values (E2) are given.at the final state be- fore perturbation and as the result of the next itera-w tion (which is based on the estimate of curvature value). iii. CONCLUSIONS: Since the same methods were utilized in this example and in comparison 5.3.2.b, it is reasonable to expect that the large variances as shown in Table 5.3 result from perturbations which are not generally on lines of curvature. There is no reason to believe that these large variances result from perturbations which are too large.or from computational inaccuracies. d. GENERAL CONCLUSIONS: Examination of the three comparisons made in this section indicate several impor- tant facts pertaining to curvature determination as it re- lates to computational efficiency. First of all, compari- sons 5.3.2.a and 5.3.2.b indicate the relative consistency and accuracy of the digital integration approach.as com- pared to that of analog integration. Since a second order system represents a very special case, the important OOH! clusions relative to curvature determination rests on third and higher order comparisons. 120 *determined for Data Sets 3 and 0 (See Appendix C). TABLE . Di ital Estimates of Curvature fora rd—Orders stem — 1th 1+1“ _iniiim“ _1‘1 _‘2 _‘3 Lav _5 __°n 0.22.1 m 3o1.1 -u82o6 -9o232 -0059 -163o9 23.679 1o37u -o111 -O962 3.1.2 .375 -.253 -.026 .032 -.070 1.191 -.962 -.979 3.3.1 -.l00 —1.288 -1.220 -.880 —l.513 .596 -.190 -.688 30302 "30292 ”0370 ‘10531 '10731 ‘025 069“ ‘o688 ”0971 3.0.1 3.803 —2.000 -2.590 -.262 -2.590 1.031 -.070 -.896 3.0.2 -2.801 —5.075 -0.328 -0.215 -12.80 .256 -.896 -.999 3.5.1 1.055 -.061 -.332 .097 -.899 1.106 -.561 -.797 3.5.2 1.195 -.793 -l. 6 —.361 -1.512 .981 -.797 -.918 3.6.3 -.259 -.062 -3.013 71.378 -.090 1.006 -.918 -.975 3o7o1 -uo633 ‘0597 -31o23 -12o15 -u0867 lo118 -0u10 -0315 3.7.3 -7.708 .352 -.996 -2.797 8.375 1.169 -.691 -. 82 3.9.1 -0.716 -.709 .057 -1.656 -. 05 1.130 -.29 -.300 3.9.2 -7.500 -1.057 .092 -2.823 -1. 90 1.159 —.30 ..330 3.9.3 —7.001 -1.385 -.367 -2.920 -1.720 1.000 -.ugo —. 5 309o ‘Eo656 '1o139 ‘03u0 ”20378 ‘1o389 ‘.980 -o 5 'o629 3o905 - 0303 ”0912 “0555 “10923 '100u1 0878 '0629 -0887 3.10.1 22.013 -1.718 .837 7.177 .001 1.300 4.115 .186 0.1.1 -.389 .153 -.l07 -.127 -.169 .966 —.769 -.901 0.1.2 -.555 -1.371 -.506 -.812 -.939 .091 -.901 -.976 0.2.1 -0.766 —.063 -.095 -l.77l -.860 1.195 -.222 -.176 nouol 10.097 “3080“ -o173 2.0u0 o810 1.255 -o666 -.003 0.5.1 -0.091 10.213 -.060 1.887 -.102 1.251 -.899 -.771 uo6o1 -3o153 -10u03 ‘ollu -1o557 -20502 oBOO -o06l -o973 The input data # indicates the data set, data point and the iteration number for that point. For example, 3.9.3 indicates Data Set 3, 110 number 9 and LOP-CM iteration number 3. 121 Examination of comparison 5.3.2.0 (third order) yields several important conclusions. Obviously, if it were nec- essary to-aocm'a-te‘lyidenttify a line of curvature, it would be necessary to attempt several, perhaps even many, perturba- tions. In fact, of all the perturbations represented in Table 5.3, only several give definite indication of being near a line of curvature (data number 3.0.2, for example). Perhaps the single most important number given for each data.point is on, the ratio of the standard deviation of the curvature estimates to the average of the absolute values of these estimates. To demonstrate the significance of this value and to indicate its use in LOP-CM, consider Table 5.0. In this table, the data is arranged according to increasing magnitude of on. 'Except for one data point (3.1.1), there is a close correspondence between.the value for on and the percent decrease of the error function‘de- noted by. The data point 3.1.1 which does not follow this general observation is assumed to be an anomaly-resu1ting from the extreme values taken.by the curvature estimates. The value for on which seems to represent the cut-off point (i.e., the point at which no improvement in cos y is noted) is approximately at on = 1.2. Stated differently, for any curvature estimate with a value of on below 1.159, the error function.is improved in the iteration based upon that estimate. It would certainly be desirable to obtain.a very good 122 TABLE 5.0 Evaluation of the Significance of Cg 1th 1+1“ "a? 3 DATA f" 0:1 cos 1' cos I -§cos y 51* ”We 59 3.”.2 .256 -0896 -0999 .103 .99 uo1o2 ou91 'o901 '0976 007 0758 817 3.3.1 .596 -.190 -.688 .09 1.613 ' 30302 069“ -0688 ‘0971 0283 0908 0.6.1 .800 -.061 -.973 .912 . 71 3.9.5 .878 -.629 -.887 .158 . 27 601 0.1.1 .966 -.769 -.901 .132 .572 ’ 3.5.2 .981 —.797 -.918 .121 .596 3.9.0 .980 -.005 -.629 .180 .332 , 3.B.3 1.000 -.350 -.005 .095 .106 .516 3. .1 1.031 -.O70 -.896 .826 .892 3o6o3 1.006 ' 'o918 -0975 0057 0695 3o5o1 1o106 -o 61 -0797 0236 0537 3.7.1 1.118 -. 10 —.51 .105 .178) .197 3.9.1 1.130 -.29 -.3 .005 .007 3.9.2 1.159 -.3 -.350 .006 .066 3.7.3 1.169 -.691 -.082 -.209 -.678 3.1.2 1.191 -.962 -.979 .017 .008 _ 272 0.2.1 1.195 -.222 -.176 -.o06 -.059 ‘ uoSol 1o251 -o899 ‘o771 -0128 -1.268 0.0.1 1.255 -.666 --°”3 -o623 ‘1-368) -1-100 0.10.1 1.300 -.115 -.186 -. Ol -. 00 3.1.1 1.37“ -.111 —.962 o 57 o 51 * determined for Data Sets 3 and 0 (See Appendix C). The data number corresponds to those in Table 5.3. * CS?) 0 = n ki av * _ -gcos 1' 6y ‘ + cos y 123 estimate of curvature, but there is a tradeoff between the computational time required for obtaining that estimate and the time required for iterating to a final state near- er the optimum. For this reason, the factor on is partic- ularly useful in indicating when the curvature estimate is close enough to use in LOP-CM. Another observation is apparent from Tables 5.3 and 5.0. As the optimum is approached, the curvature estimate generally improves. This is especially apparent in cone sidering on for data points 3.1, 3.0, 3.5 and 3.9. In.summary, the experimentation.of this section indi- cates: 1) For second order systems (where the boundary is the line of curvature), the curvature estimates correspond closely to that expected from reach, able set geometry. 2) Estimates of curvature obtained by digital inte- gration.are more consistent and more accurate than those determined by the analog computer. 3) The factor a is a good measure of both the loca- tion of a fim 1 state perturbation relative to a line of curvature and the usefulness of the cur- vature estimate. 0) As the optimum is approached, the curvature esti- mate generally improves. 5.3.3 Effect of the Basic Curvature Formula Choice In this section the effect of the basic curvature for- mula choice (Equations 3.120, 3.121, 3.123 and 3.120) on the rate of optimization convergence is considered. Com- parisons were made for LOP—CM using both Subalgorithms 0 (i.e. with and without standard deviation evaluation of 120 the curvature). a. Go arison.Basis: ES—l, to = O, T = 20 secs., Input Data Sets 1 and 2 (See Appendix C): ES—2, to = O T = .5 sec., Input Data Sets 3 and 0 (See Appendix C). b. Summa ‘2; Copparison.Results: The average num— ber of itera ions requ red to attain a local optimum were compared for each of the basic curvature formulas. The following results were obtained: 12 (Equation L123)-----7.7 iterations Eav2(Equation 3.120)---7.7 iterations Tc, (Equation 3.120)----—8.3 iterations 'Eav (Equation 3.121)----9.3 iterations. In as much as and k.had the same average number of iterations, a fur her comparison of these two estimates was made. This comparison was based on ES-2, a third order system. In this case the percent improvement in the error function for each curvature move was consider- ed. For E the percent improvement was 27.9 while it was c. CONCLUSIONS: If the estimates of the curvature were accurate, the choice of basic curvature formula would have little effect on the algorithm efficiency. In the case of third and higher order systems, where only approximations to lines of curvature are achieved, the choice of basic curvature formula does alter the algo- .rithm convergence. The data of this section indicates that 2, the overall average of curvature estimates, is the est choice with'E'being the second choice. 5,3,0, Copparison of ng-CM'Subalgprithms 0a and 0b, Within the structure of LOP-CM, two potential sub— algorithms are introduced to determine curvature. sub— algorithm 0.a provides an estimate of k without consider- ing the accuracy of the estimate. On the other hand, Subalgorithm 0.b considers the normalized standard devi- ation on and rejects any estimate for which on is greater than a specified value. In the computations of this 125 section, the limit is specified as 1.2 corresponding to the results indicated in Section 5.3.2. a. Copparison Basis: ES—Z, t0 = O, T = .5 sec., Ins put Data Set 3. ' b. Summapy p£.Comparison Results: In as much as there is tradeoff between computation time spent estimating cur- vature and time spent iterating, only a maximum of n (3 for this third order.system) attempts at obtaining a suit- able estimate of curvature are allowed. For this case, the average number of iterations to achieve a local opti- mum without considering cn_was11.0 whereas utilizing a to evaluate the curvature estimate improved this averagg number to 9.3. Consideration of the average percent improvement in the error function.per curvature move (LOP- CM iteration), substantiates these figures. If a is not considered, this percent is 30.7 whereas utilizatibn of the on factor gives an improvement of 03.3%. c. CONCLUSION: The data of this section indicates that rejection 0f poor curvature estimates (high a ) im- proves the convergence of LOP-CM. This effect woufd be even.more marked were the allowable a smaller and the allowed number of attempts at obtaining a good curvature estimate 1arger.. 5.3.5 Cepparison of Perturbation Direction Alternatives In step 5 of LOP-CM, three alternative equations (3.80, 3.90 and 3.108) are listed for determining the di- rection of the perturbation of the state at the final time. It is the purpose of this section to present ex— perimental results comparing these three choices. a. Go arison.Bpsis: ES-l, to = O, T s 20 sec., Input Data Sets and 2; ES-2, to = 0, T = .5 sec., Input Data Set 3. ' h. Summapy‘pg Copparison.Resu1ts: The average num— ber of iterations required to attain a local optimum for these three alternatives were as follows: i. Equation.3.80 ------ 7.20 iterations i i . Equat ion 3. 1.08------ 7 . 31 iterations our I\‘ j.“ .. .r 126 iii. Equation.3.90 7.81 iterations. Similar results were obtained for the average percent de- crease in the error function for each LOP-CM iteration: i. 02.0%, ii. 00.7% and iii. 20.9% 0. CONCLUSIONS: Very little difference was noted between the effect on.LOP-CM.of using Equations 3.80 and 3.108. On the other hand, Equation 3.90 gives less sat- isfactory results. It should be recalled that Equation 3.80 was developed for regions far from an.eptimum, Equa- tion 3.90 for regions near an.0ptimum and.Equation 3.108 for general applicability. In the computations of this section, LOP-CM was used until E reached .1 and then LOP-DS was employed. Had LOP-CM een.used for greater convergence, the utility of Equation 3.108 would have been more readily apparent. 5.3,6 Copparison of Perturbation.Step Size Alternatives In.the previous subsection, the direction of the final state perturbation was considered. In this section, its relative magnitude, c, is discussed. The important com- parison.to be made is that of a constant step size as com— pared to a step size dependent upon the error of the unper- turbed boundary point. Specifically, the two candidates selected are c = .2 and c = (l + cos y). a. Copparison.Basis: ES—l, t = O, T = 20 secs, Ins put Data Sets 1 and.2 (See AppendixOC); ES-2, t0 = 0, T = 20 sec., Input Data Set 3. b. Illustrative 1e: In Figure 5.5 an example iteration sequence is g ven. or both the constant and variable step factors. c. Summary‘p; Copparison Results: The average number of iterations required to attain.a iocal Optimum for cone stant step factor (c = .2) was 6.68, while the error de- dendent factor yielded a lower average number of itera— tions--5.92. The average percent decrease in.the error for each LOP-CM iteration. ave corresponding results: constant step factor-03.0 , error dependent step ' factor-09.l%. 127 10.,5 “ U independent step factor j” 0 error dependent step factor 1+ . = 01"? I -<>1.'. .- .L eOOI-fi 1 ”1.9.1. _____ ._ _. - _.._ ._ _. _.,_. _. _. _. _ , ._ “ “ bound :- i 1 ‘ ' :1 .1 3‘ : 1 2 3 0 5 6 7 a iteration'number FIGURE 5,5 Cogparison of Perturbation Stgp Size Alternatives d. CONCLUSION: The error dependent step factor (c = l + cos y) is more satisfactory in.achieving coup vergence than the constant step factor. 5.3.2 Analog Error and the Integration Correction Routine Section 3.6.2 pointed out that significant errors are encountered in analog integration of the differential equa- tions of Example Problem 1. As a consequence, correction factors a and a are introduced in Step 8 of LOP-CM. It is the purpose of this section to give‘examples of these er— rors, consider their origin and evaluate the effectiveness of the correction factors. a. Go arison Basis: ES—l, t =0, T = 20 secs., Ins put Data Set!“ !'1""'(See ip"'p'e'nd'"ir C), Roth Analog and Digital Integration. 'L 9129' 128 b. Illustrative Examples: Figure 5.6 is given to re- late the various state and adjoint vectors and correction factors. In Table 5.5, examples of these various vectors are given for both.analog and digital integration. values of the various vectors depicted in Figure 5.6 are presented as are error function.values at the start and at the conp clusion of the iteration. c. Summarz‘gg Results: The computations of this sec- tion yield several important r sults. First of all, COD! sideration.of the values for a and 31 in Table 5.5 demon— ‘ strates that analog computation errors are indeed important. Furthermore, at points'on R(T) far from.an Optimum, these. errors appear to be larger than.for points near an optimum. Two digital iterations are also presented and it is readily apparent that the digital integration routine does not pro- duce these large errors. Comparison of LOP~CM with and without the use of the error correction factors a1 and 51 demonstrates their use- fulness. 'Without the error corrections, the average number of iterations required to determine a local Optimum is 9.75 whereas the introduction of the correction factors lowered this to 7.00. Comparison of the relative efficiency of LOP-CM with analog integration and LOP-CM’with digital integration is given in the next section. d. CONCLUSIONS: It is, of course, obvious that ana- log computation errors do occur and that the correction factors a1 and 31 do help compensate for these difficul- ties. The significant item to determine is the cause of these errors. At first the author considered the possi- bility of incorrect patching on the analog computer or incorrect conversion to reversed time integration. After eliminating these as possibilities, analog-digital con- version and digital-analog conversion.were considered. These too, do not appear to be significant since the error is only on the order of one-half of one percent. Further consideration indicates that the errors enp countered are a result of the interaction.between several factors: 1) The nature of the differential equations and the large time interval involved in these comparisons (20 seconds). 2) The sensitivity of the trajectories to switching time. 3) The errors created.by analog equipment-- comparitors, integrators, multipliers, electronic switches, etc. 129 initial final states states —i Pg forward time .1); perturbation PT I ,, . i i ‘1 x0 integ’ration IT IT re rse rev rse t me ti e inte ation integ ation ERROR _1 EVALUATION 13 x0 1 -i p. P. i =-~,'i_._ L.“ 3 to l i .. i ._ i B - fio po L—i>CORRECTION «f .. ‘i 1 1(1) - IO - a. 'i - "i _ 1 5° - po . B ———-—>COMPARISON 1: 1- 6x0 20 x0 6p: = 33 - p3 ‘ DECISION Pg+l For Instance: 1+1 _ i i + 26pi orp3+1 = Di 0 0 FIGURE 5,6 Analog Error EvaluatiOn and Correction Flow Chart 130 TABLE 5,: ForwardpReverse Time Integration.Error Evaluation IE: 1141 A1-2 D1-1 D1—2 iu-l Au—z Au—a Ah—h 115-1 115-2 AIS-3 AlS-h IO -10. _ 5, -10. -5, -10. -5, -10. -5. 6. 2. 6. 2. pf. 0.000.. 10.000 l.h7h 6.779 0.000 10.000 2.166 9.769 10.000 “'3 o 000 -6.932 -6.6h2 -9.5 6 “207 3 ”90883 -l.uOO '50000 -5.000 “90389 “50995 -9.862 -1.892 ‘9o9h6 0.271 error 1.087 .600 1.277 1.219 3.1962 1 .052 .393 .286 1.862 .UBO .293 .056 Hi .131 n.292 .063 3.218 .001 .002 .001 .002 -.360 “0755 '0360 -1.h87 -0360 -1.63h -ouo9 -l.2h3 —1.167 .630 -.630 .386 ”0532 -0396 -0581 ‘0200 i 8.212 ““0293 .021 ”0358 ‘0029 -.000 -.026 ~.006 -0139 “0&29 -0271 .890 “0115 1.218 -.060 .927 n.036 -2.h57 .156 -.51h .002 .780 .373 1.567 new i error .600 n.619 1.219 .689 1.052 .393 o 286 .167 .uao 2293 .056 .030 *qui indicates the ith iteration for the 3th data.point using analog integration. the same data point is indicated by Di—i. Digital_integration for 131 The first two factors multiply the effect of the er— rors introduced by the analog equipment. Take, for exam- ple, the state and costate trajectories plotted in.Figure 5.7. While the initial and final points are not far re- moved from.one another, the trajectories traced out do by no means represent a shortest path from the initial to the final pOints. Also, four switching times occur for this extremal trajectory even though only a second order system is involved. Even slight differences in.awitching time due to comparitor hysteresis, slight multiplier inaccura- cies, etc. can have a significant effect on trajectory be- havior. 5.2,8 Comparison of Analog vs Digital Integration in LOP-CM The previous subsection strongly indicates that the er- rors due to inaccurate analog computation have a signifi- cant effect on the algorithm convergence. The introduction of the correction factors a1 and 31 improve this conver- gence. It is the purpose of this section to determine if the resulting improved LOP-CM with analog integration is as effective as a completely digital LOP-CM. a. Comparison.Bgsis: ES—l, t0 = 0, T = 20 secs., Input Data Set 1. _ b. Summa ‘2; Results: The average number of itera— tions required to attain.a local Optimum for LOP-CM with corrected analog integration was 6.25..‘With.digital inte- gration, the average number was only 5.33 for the same data. In addition, LOP-CM with analog integration failed to con, verge in two cases because of analog inconsistencies where- as convergence was always achieved using digital integraa t on. c. Illustrative E les: Two typical examples of the convergence of LOP-CMZwith.analog and LOP-CM with dig- ital integration are shown in.Figure 5.8. d. CONCLUSION: Use Of digital integration in LOPaCM is far more consistent and accurate. As a result, it con- verges in less iterations than.LOP-CM with analog integra- tion. It has the disadvantage of requiring more digital computer time. A judgmentl as to which.approach is best would depend on the equipment available, the accuracy of 132 P2 1" -'42-6-0-‘t’4-I8- \ \ ‘;\~‘~._ h x = (-lo,—5)T p = (0,5)T t -t indicate switching ° ° 1 14 times 1:1: .118 sec. t2-3.35 sec. t3=ll.26 sec. tu=13.08 sec. FIGURE 5,2 Exagle Trajectories for a Second-Order sttem 133 0 digital D analog .OIJ' p =[5] .007” O 3 0002" . 001 He» N UM C“ up. 0\ \l'” on T 0 digital n anlaog . \l I M O H .1; fi‘ ’°7T . t ' h - *- .02o x g [g] ‘, 1 ”01% O p. {3] o 0021' v u [I £175 6 7‘ 8% iterations . 001 q. \o.. .... Ci” :A+ {5'1 1 2 FIGURE 5,8 LOP-CM Convergence for Analog or Digital Integration 13h the desired results and the computation.time available. It is the author's opinion that," in general, LOP-CM with digital integration.wou1d.be the better choice. 5,h Comparison of LOP~CM with LOP-BS In the preceding section many LOP-CM alternatives were considered. Choice of the most suitable of these alterna- tives yields an effective optimization algorithm. It is the purpose of this section to summarize the comparisons made between.the resulting LOP-CMIand LOP-BS, a direct search algorithm. . ~ The comparisons were made for Example Systems 1 and 2, with selected data points from Appendix C. The results of several sequences Of comparisons showed an.average of h.17 iterations required to achieve each local optimum us- ing LOP-CM. The direct search routine, LOPeDS, required an average of 7.12. Thus LOP-CM represents a definite imp provement over a direct search routine. In Figure 5.9 an example of their relative convergence is given. 5,5 Global thimization.Examples All of the preceding developments and comparisons were aimed at producing an efficient algorithm for solving the modified problem HP. As a result Of the comparisons summarized thus far in this chapter, an efficient form of LOP-CM is identified. Using this ”Optimized" LOP-CM, the Global Optimum.Procedure GOP given in Section.h.7 is now applied to example problems. This procedure is based on 135 5 7. 2 - D — LOP—CM ' o - LOP-DS l.,-- 07 it _ ~10 _ .2 o n ‘0 “ _ 5 p0 ‘ g .1 . l .07t .02? “ ‘. - 001-Ji- . ' .007. 0002'“ '0010 i 2 i 5 n ‘ 9 10 ll 12 13 5 l 7 iterations FIGURE 5,2 Comparison of LOP-CM.and LOP-D8 Convergence the generation of N local minima Of the error function.by starting at N arbitrary initial adjoints and utilizing LOP-CM until a local minimum is obtained. Of course, the probability that a global Optimum is included in this re- sulting group of local minima increases as the integer N is increased. In this section, the Global Optimum Procedure is ap- plied tO the example problems Of Section 5.1. The integer N is fixed at 10, thus ten arbitrary initial adjoints are specified and ten.minima of the error function.iE5 are lo- cated. Many or all of these minima may coincide. In fact, if only one minimum Of E5 exists, i.e. the global optimum is unique, then all ten.must coincide. 136 Some Of the special problems discussed in.Chapter h are encountered in the computations of this section. They are discussed as they occur. In.addition, typical itera- tions for LOP-CM are presented for some of the gIObal prob— lems considered. 5.5.1 Examples for ES—l: 2nd—0rder sttem, Scalar Control Example Problem 1 is particularly useful as an example Of global optimization since it was specifically developed to produce significantly nonconvex reachable sets. The shape of these-sets, of course, depends on the initial state selected and on the time interval of integration. The effect of increasing the time interval is demonstrated in the next section (time optimal control). In this sec- tion the global problem is considered for several reachp able sets for ES-l corresponding to different initial states. The time interval is fixed at T = 20 seconds. The first initial state to be considered is ’ x0 = (-10,-5)T. The reachable set has alreadybeen shown. in Figure 5.0. The ten resulting minima of GOP are given in Table 5.6. Consideration of these data points shows that GOP selects three possibilities for the global Opti- mum--data points 1, Band-8: data points 2, ’4, 5, 7 and 9; and data.points 6 and 10. Of these, E1 is minimum for data number 1, 3 and 8. Hence they represent the global Optimum as consideration of the reachable set demonstrates. 137 TABLE 5.6 Application of GOP to ES—1.with 1gp: (:lo,—51T 22M .311. gm. him! E2 .21. 1 .526 -8.268 8.281: -.9993 .0057 2 44.10“ -9.106 9.989 r -.9999 .0013 3 .526 -8.268 8.28u -.9990 .0078 u -u.2ua -9.039 9.987 -.9995 .0052 5 -u.lou -9.108 9.989 -.9998 .0021 6 -8.151 -2.lul 8.u27 -.9999 .0009 7 -u.lou +9.108 9.989 -.9998 .0022 8 .586 -8.26l: 8.285 -.9997 .0027 9 A: . 2'48 -9 . 039 9 . 987 - . 999’4 . 0056 10 —8.039 -2.539 . 8.180 -.9997 .0023 Consideration Of the other two minima of E5 shows that one is a local Optimum (Data pOints 6 and 10) and that the other is a false optimum as discussed in SectiOn 15.1. Both of these minima occur. on concave surfaces. Application Of Theorem.u.1 to the first point (6 and 10) demonstrates that it is a true local Optimum. The curvature at this point is .02; thus K0,, is .02 since the boundary .is $2 line Of cur- vature. But [Ilel is 8.11 or l/[Ilel is .119. Thus Equa- tion “.10 applies: _ _ no, = .02 < .119 ,s 1/||x.,.||. (5.22) and thus 1:1. is a local Optimum. For data points 2, '4, 5, 7 and 9, the curvature is. .22; hence Kcv is also .22.. But “IT” is 9.99 13111.18 l/llell is approximtely .1, thus Equa- tion 11.6 applies: 138 Kev = .22 > .1 = 1/[IxT|I, (5.23) or xT in this case is not a local Optimum. The second initial state considered for ES—l is x0 = (5.5)T. The resulting reachable set is shown in Figure 5.10. While this set is nonconvex, it does not produce as interesting results as the previous example. In this case, there is only one local Optimum (located at xT = (-2.’+,.5)T ) and thus each iteration of GOP se- lected it. It, of course, represents the global Optimum. Also illustrated in Figure 5.10 is an example se- quence of iterations. In this case, it took 5 iterations to converge to the Optimum. The final state perturbations for the first four steps are shown. 5,5,2 Emmmples for-Higher Order sttems While application of GOP to ES-l resulted in.aeveral minima of E5 being achieved, application to the other ex- ample problems yielded unique optima. The resulting Opti- mum final states are given in.Tab1e 5.7 as are the corre- sponding values of the error functions. Given in Figure 5.11 are several examples of the change of error with each 0 . iteration for ES-2 for various p0. 5,6 Time thimization Examples In Section.h.8 the application of GOP to time optimal control problems was introduced. It is the purpose of this section to apply the resulting procedure to several example 139 R(20 ‘ '“W‘i 6x orir in R(20) for E:3—l , \ I 0 x1 -6 -5 .. .1 o FIGURE,5.10 Reachable Set and Example Iterations for LOP-CM lho TABLE 5.2 -Application of GOP to Higher Order Systems Example Problem ES-3 ES-h .5 .5 .5 .5 .5 # xT Illill -.9998 -.998 -.998 -.9996 —.996 -.9990 '0962 .001 .009 .001 .0011 .00“ .001 .007 .007 .062 1'41 T = .5 seconds A g 1. Io = (-3,-2,-1)T V i ' various p0 O 'A 001+ a? n I A “i .I l A l. A L J: °°1 1239.?6778910 iterations E 10. T = .5 seconds “ is. (2"2’0) ” “ 0 various pO ’1 U H II . . - - p o°01L——L 1 1 1 L l 2 3 h 5 -6 .7 8 9 10 ‘Jiterations FIGURE 5,11 Example Iterations of GOP for Example Problem 2 1H2 reachable sets for ES-l. Example Problem 1 is chosen since the resulting reachable sets are particularly nonconvex and interesting in their behavior. Unfortunately, however, the sequence of optimum final states for these reachable sets only asymptotically approaches the origin as' :‘time is in- creased. ‘ The first example considered is for the reachable set corresponding to the initial’state xb = (-10,-5)T. In Ta- ble 5.8 the results of some time increments are given. It should be observed that [I IT” does approach the origin, but the convergence is asymptotic. Likewise, it should be not- ed that the sequence Of llell's is not monotonic. These results are illustrated in the reachable sets corresponding to various times as given in Figure 5.12. From this figure it is easy to see the spiraling which R(t) does around the origin and its increasingly nonconvex nature. The locus of optimum final states (for the minimum distance problem) is also plotted in Figure 5.12. The second example corresponds to the initial state x0 = (5,-5)T. In Table 5.9 the results for various final ,times are given. The reachable set at various times is plotted in.Figure 5.13. Again, there is the spiraling ef- fect noted in.the previous example, with the reachable set becoming increasingly nonconvex as time increases. Consideration of the initial adjoints listed in Tables 5.8 and 5.9 shows that there is a consistency in the manner 1“3 TABLE 5,8~ Application 0: ‘TOP to ES-gl. xx = (-10.-5)E Time 0.00 5.00 7.10 9.65 12.06 15.75 19.17 25.00 28.03 30.09 33.53 36.79 38.“? 00.00 01.09 “2.86 00.26 05.60 06.92 08.22 “9.51 50.8“ JET.”— 11.180 20.980 18.009 12.060 6.639 9.899 9.792 0.321 3.955 “.551 “.309 2.812 2.000 1.809 1.659 1.620 1.519 1.002 1.391 1.33“ 1.303 1.358 xlT‘ -3.5“0 0.757 “.810 6.519 5.“?3 1.001 -“.O77 -3.881 -3.662 -l.2“5 -O.l70 0.268 0.659 0.903 1.098 1.098 1.196 1.098 1.098 1.196 1.293 321 -5.000 20.679 15.99“ 11.099 1.205 -8.250 -9.7“2 -1.“89 0.756 2.710 0.126 2.807 ’2.026 1.68“ .1.391 1.196 1.0“9 0.805 0.85“ 0.756 0.610 0.015 Plo 5.000 5.929 6.78“ 5.955 5.972 5.998 5.999 5.900 5.882 6.000 3.333 0.117 3.817 3.902 “.003 “.003 0 0.980 5.““9 6.377 6.038 5.962 5.962 P20 -10.000 0.071 0.836 0.730 0.999 0.110 0.618 0.100 -0.007 0.000 3.533 “.360 3.070 0.322 0.003 “.“33 3.500 3.2“9 1.072 0.670 0.666 0.666 1““ 80 6O J7 6; /Locus of Opti Fi'nall State :Tm “O V FIGURE 5.12 R(t) for E84, 5,, = g-10.-5)T. Various t 1“5 TABLE 5.9 Applicatiom_of,TOPLtOJES-l, IQ = (5. -5)T Time. 0.00 5.00 5.78 6.12 7.50 8.02 9.38 10.33 12.80 1“.“9 15.29 16.13 17.05 19.16 21.50 23.7“ 25.65 27.23 28.66 30.0“ 32.02 3“.12 35.37 36.03 llgmll 7.071 7.758 7.202 6.716 6.102 5.679 5.081 0.022 3.089 2.760 2.701 2.803 2.960 3.000 2.702 2.107 1.587 1.3“8 1.202 1.103 1.038 0.911 0.929 0.950 5.000 1.3“2 0.756 0.170 -0.“15 -1.001 --I.“89 -l.928 -2.758 -2.758 -2.710 -2.661 -2.563 —2.075 -1.538 —0.366 -0.020 0.122 0.122 0.317 0.015 0.317 0.063 0.512 -5.000 -7.601 -7.202 -6.71“ -6.128 -5.590 -“.858 -3.979 -1.391 —O.170 0.015 1.001 1.089 2.172 2.221 2.075 1.586 1.302 1.196 1.098 0.952 0.85“ 0.805 0.805 5.000 -5.992 -5.983 -5.99“ -5.990 -5.99“ -5.988 -5.990 -5.035 -6.302 -5.3“9 -5.577 -5.577 -6.“53 -6.21“ -5.955 -5.893 -5.881 -5.950 -5.950 -5.968 -6.679 -5.999 25.999 .3111__£2L”_I:10__1120_ 5.000 0.395 0.251 0.251 0.251 0.251 0.188 0.002 -0.13“ 1.556 2.211 2.211 2.211 2.011 1.081 01.137 -1.125 -0.995 -0.766 -0.766 -0.612 —0.036 -0.032 _0.032 l“6 R(20) R(15)7 -21 ~10 20 Path 01 mum FinTl eState ? Opti- ~30 .00 FIGURE 5.13 R(t)fifor ES.-1. 11:,, = (5.-5)T. various t l“? in which.the optimum initial costate changes as time is incremented. In.the examples of this section, use is made of this fact. For the first local Optimum (T = 5 seconds) LOP-CM is employed. fFor succeeding Optima, the change in the initial adjoints is sufficiently small to optimize using LOP-D8. In fact, in several instances the 1th p; is also the 1+18t p3 (See Table 5.9, time = 6.12, 7.50 and 8.02). 5.2 Summar: and Conclusions It is the purpose of this section to present a general summary of the theoretical developments, computational re- sults and comparisons of the previous sections of this the- sis. Conclusions which.are given in previous sections of this chapter are also briefly summarized in this section. Additional comments or conclusions which seem appropriate are also made. As indicated by the title, the purpose Of this thesis is to consider the computation of Optimal controls for DOD! linear systems. Rather general nonlinear systems are al- lowed. The restrictions placed on these systems are given '3 ‘ww.w:fiQ-FI . .s... “-4 in Section 2.2. The concept of the reachable set is intro— duced and investigated. A number of related definitions and results are givenxin.Chapter 2. The minimum-error regulator prO‘Bie-sism principally treated and computational methods are developed to solve this problem. The problem is somewhat difficult when ‘ .slmw’" ' ”1.-A xv. 1“8 nonlinearities are allowed in the state equations. The nonlinearities have a two-fold effect on the computation of an optimal control for the minimum distance, problem. First Of all, they introduce the possibility of the occur- rence of several or even many local optim since the re- sulting reachable set is not convex. Secondly, the deter- mination of a local optimum is made more difficult since the adjoint system becomes dependent on the state. Although many reachable set-oriented methods have pre- viously been utiiized, their application to nonlinear sys- tems has been limited. One of the most restrictive limita- tions is a result of the dependence of the adjoint system of equations on the state variables. As a result, the di- rect determination of an initial adjoint given a final ad- joint is not possible without knowing the corresponding state trajectory. To overcome this difficulty, the author decided to place the emphasis of his computational approach on the initial adjoint. Once the initial adjoint is specified, it is possible to compute the state response, a maximal trajectory, and to identify a boundary point on the reach- able Bate An intrinsic part of any effective computational meth- od is a means Of evaluating each iteration and of deter-. mining when an optimum has been achieved. This is the pur- pose for the error function. One obvious error function 109 for the minimum distance problem is E1, the norm Of the final state. Other error functions are available and are introduced in Chapter 3. It is shown that for the minimum distance problem, the final state and the final adjoint are collinear at the Optimlmi. This leads to a secondverror g =' ‘ E2 003 Y I p 9 (502,4) where x is the final state and p is the final adjoint. fmetion: This error function is more sensitive to changes in extre-g mal trajectories than is IIxII. Also, [lel may just be ap- proaching an unspecified minimum whereas E2 approaches -1 at the Optimum. For these reasons, E2 is usually employed in combination with E1: E5 = (1 + E2) E1. (5025) One approach to the solution of the minimum error reg- ulator problem is a direct search (related to gradient tech- niques) method based on varying the initial adjoints such that the error flmction is improved. A direct search pro- cedure is given.on pages “1 and.“2 and a flow chart of this method is given in Figure 3.1. Since the direct 'search technique is inefficient, other possible methods are considered. To develOp another method, the determination of the optimal control can be viewed as the determination of .a sequence of initial ad- joints which, in turn, determines a sequence of extremal endpoints starting at an arbitrary final state and 150 terminating at an optimum. Thus the problem is one of de- fininga path on the boundary of the reachable set, specif- ically a path for which the error flmctions decrease mono- tonically. To show that such a path exists and to give insight into the nature of this path, various principles and facts from differential geometry are utilized. Specifically, it is possible to prove that on a locally convex surface, a satisfactory path can be constructed from lines of curva- ture. Requiring that the path consist of lines of curva- ture gives a relationship between the perturbations of the final state and the final adjoint on the boundary of the reachable set: I: 6x + 0p ,= O, (5.26) where It represents the curvature of the surface at the pOint'x on 'the line of curvature. This aids in the imple- mentation of a procedure using this geometric approach. Chapter 3 develops this procedure and considers the al- ternative choices which are encountered. Several interesting facts develop as a result of these geometric considerations. To prove the monotonicity Of the error function E2, it is necessary to also prove the mono- tonicity of the error function E1,.along the desired path. As a result of this proof the importance of the curvature of the surface is demonstrated (see Equations 3.66 and 3.72). These facts lead to Theorem “.1 which analyzes the 151 effect of concave surfaces on.the procedure. The alternative choices available in the procedure can. generally be classified as those which are available on.the boundary of the reachable set (perturbation of the final state and the final adjoint to yield a better estimate of the initial state, hence Of an.eptimum) and those which are important at the initial time (perturbation of the initial adjoint). The computations summarized in Section.5.3.l indicate that the final time perturbation.is most effective if both the final state and the final adjoint are perturbed as re- lated through Equation.5.26. The perturbation.af xT is considered to be the independent perturbation, hence it is important to consider the direction.and magnitude of this perturbation. Sections 5.3.5 and 5.3.6 consider these choices and indicate that the best means of perturbing the final state IT is: . OxT = -c(xT-||x.r|| "—15%" E2) (5.27) where c represents the step size. Consideration of the effect Of the step size demonstrates that the method is most effective if c is made dependent on the error at each iteration. The perturbation of the final adjoint is related to the perturbation.of the final state through Equation 5.26, hence determination.of the curvature k is important. Sec- tions 5.3.2 through 5.3.“ consider the various methods 152 of estimating curvature and of evaluating the validity of the estimate. Examples Of estimates Of curvature are given in these sections and shown to correspond to the values ex- pected from reachable set geometry. The analog correction factors a and 8 introduced in Section 3.6.2 are shown to improve the convergence Of LOP-Cliusing analog integration. While the experimentation of this thesis was not intended to evaluate the relative merits Of analog or hybrid computation versus digital com- putation, the necessity of these correction factors intro- duces this subject. The advantages and disadvantages of both are apparent in.the computations of this thesis. No decision is clearly Obvious concerning their relative mer- its, but digital computation is more reliable. While the hybrid computer has the advantage of Speed of integrations, it has disadvantages also. Accuracy is limited due to the analog elements employed and to analog- digital and digital-analog conversion.limits. Because of the nonlinear systems considered and the discontinuities introduced by the controls, the trajectories are very same sitive to hystersis of comparitors and switches, inaooue racies in.multipliers, integrators, etc. Thus there is a large amount Of sensitivity to the control switching times. On the other hand, the digital computer is consistent and repeatable in its results. The accuracy of the inte- gration can be controlled by controlling the step size. 153 However, the major drawback of the digital. computer is its speed of integration of the differential equations. Thus, as previously indicated,.the decision as to whether to use analog or digital computers for the differential equation integration must be based on the nature of the system be- ing considered, the accuracy desired, the computational equipment available and the speed of computation desired. For the equipment available to the author, total digital computation was preferred because of its accuracy. Both local Optimum procedures, LOP-CM and LOP-DS, were utilized and performed as expectedm-converging to lo- cal optim. The )convergence of LOP-CM, as expected, was more rapid than that Of LOP-DS, In none of the experimen- tation performed for this thesis were any limit points of LOP-CM or LOP-DS encountered. The only failures of these methods to converge can be traced to analog integration inconsistencies. Some of the special problems introduced in Chapter “ were encmmtered. Special Problem 1, relating to concave curvatures, was encountered in Example 5.5.1 and the re- sults verify those predicted by Theorem “.1. NO flats or singularities were noted, nor were false bomidary points and interior boundaries. Corners were produced in several examples, but caused no special problem. In most cases, in fact, the corner is the optimum. Generating a sequence of local minima. did locate the 15“ global optimum in the examples considered (See Section 5.5). Only one optimum was computed for the reachable sets cor— responding to Example Problems 2 through “. While this is a disappointment from the standpoint of effectively testing the Global Optimum Procedure, it is encouraging in that it indicates that the reachable sets corresponding to many nonlinear problems are not extremely badly behaved. In Section 5.6, the time Optimal control problem is , considered and example computation is done. Examination of this computation demonstrates the mobility and changing shape of the reachable sets with increasing time. It is important to note that for the examples considered, the initial adjoint (corresponding to the optimal control) did not significantly change as time was" increased. In sumry it can be concluded that the theory of Chapters. 2 through “ provides the basis for an effective computational procedure. The consideration of the algo- rithm alternatives in Chapter 5 verifies the choices made in Chapter 3. The resulting LOP-CM is effective in the de- termination Of local minima and, when employed as part of a global procedure, determines the global optimum. Certainly there are many areas Open for future in- vestigation. Possibly other principles Of differential geometry can be brought to bear on the Optimization prob— lem (utilization of geodesics, for instance). Other appli- cations Of LOP-CM, GOP and TOP could be considered as well 155 as comparisons with other existing methods. Possible future investigations and extensions are suggested in Appendix D. BIBLIOGRAPHY “ ‘0‘- _L-A rr A1 B1 BZ B3 BS B6 B7 D1 E1 BIBLIOGRAPHY Athans, M. and P.L. Falb, "Optimal Control," McGraw-Hill Book Co., New York, (1966). Barr, R. 0., "Computation of Optiml Controls by Quadratic Programming on Convex Reachable Sets ," Ph.D. thesis, University of Michigan, 1966. Barr, R. O. and E. G. Gilbert, Some Efficient Al O- rithms for 9, Class 9; Abstract thimization Prob ems Arisimg-fp thiml Control, To Appear in IEEE Trans- actions on Automatic Control, 1969. Bellmn, R. "Dynamic Programming," Princeton University Press, Princeton, N.J., (1957). Bellman, R. and R. Kalaba, gffiamgProgrammimg, Invariant impeddimg 9mm uas nearizat on: omar- Isons _a_ng lpjterconnections, in ”Computing Methods In Optimization Problems," Academic Press, New York, (196“), pp l35—1“5.. Bellmn, R. and R. E. Kalaba, "Quasilinearization and Boundary Value Problems," American Elsevier Publishing Co., New York, (1965). Benson, R., "Eucl idean Geometry and Convexity,” McGraw-Hill Book Company, New York, (1966). Bryson, A. E. and W. F. Denhamri glzegpest-Ascent Method for Soivimg timum Pro rammimg Problems . J. Appl, Mech. Ser E. 01 29, 519620.131). 257-257. de Jong, J. L., “Application of Picard's and Newton's Methods to the Solution of Two-Point Boundary-Value Problems in Optiml Control Theory," Ph.D. Thesis, University of Michigan, 1967‘. i , ' Eaton J . H. Anlterative Solution to Time-gp‘timi ContrOl‘, J. 11555. AnaI and Ipp1., .1701— ' p"_3"‘—p 29.300; Errpta and Addenda, J. Math. Anal. 11nd App1., Vol. , 617): pp 157-152. 156 F1 G1 G2 H1 H3 H“ H5 H6 H7 157 Fadden, E.J., ”Computational Aspects of a Class Of Optiml Control Problems,“ Ph.D. Thesis, University of Michigan, 1965. Gilbert, E. E., An Iterative procedure for CO uti , the Minimum if, _a_ ufiratIc Form pm 9, Eonvex ,et, AM. J. on Control, S'er. L, VO ’3, NO. 1, (1966), pp 61-80. Gerretsen, J. C. H., "Lectures on Tensor Calculus and Differential Geometry,” P. Noordhoff, Groningen, The Netherlands, (1962). Hadley G., "Nonlinear and Dynamic Programming," Addison-Wesley, Reading Mass., (196“). Halkin, H., 1_A_e_j;hod _9_f_ Convex Ascent, in ”Computing Methods .in Optimization Problems ," Academic Press, New York, (196“), pp.‘211-239- Hestenes, M.R. and T. Guinn, g Embeddimg Theorem _fpm Differential E uations, J. of Opt mization Theory and AppIications, Vol. 2,‘ NO. 2 (March 1969), pp 87-101. Hestenes, M. R., ”Calculus of Variations and Optimal Control Theory," John Wiley and Sons, New York, (1967). HO, Y. E. A Successive Approximation Technigue for thimal ControI S stems Subject 1159 Imp ut Saturation, Trans. fiSME,‘ J. Res c Eng., Series , . 82, (1960), pp. 33" 00 Hooke, R. and T. A. Je‘eves, "Direct Search" Solution of Numeric l and" StatisticaI PrObIems J. Assoc. 5311p. Mac ., 701 8, NO. 2, (Aer 1961), pp 212-229. Hsu, J.C, and A. U. Meyer, "Modern Control Principles and Its Applications,” McGraw Hill Book Co., New York, (1968). Knudsen, H. R., All Iterative .P ocedure for Co uti thimal Controls IEEE Trans. on Au 0. Contro . V01. ‘10-9’ NO. 1: (19.6“), pp 23-300 , Kepp, R. E.‘ and R. McGill, _S__;vera1 Trajectomz $timiza- tion Technigues, in ”Computing Methods in Opt mizat on Problems," Academic Press, New York, (196“), pp. 65-89”. Ku, Y. H., ”Analysis and Control of Nonlinear Systems,” The Ronald Press Co., New York, (1958). i L1 L2 N1 112 P1 P2 R1 81 T1 W1 21 158 Lapidus, L. and R. Luus, ”Optiml Control of Engineer- ing Processes," Blaisdell Publishing Co., Waltham, Mass., (1967). Lee, E. B. and L. Markus, "Foundations of Optiml Theory,“ John Wiley and Sons, New York, (1967). Neustadt, L. W., Szm‘thepizimg Time timal Control §zstems,uJ. Math. Anal. and App1., VO 1, (1960), Pp. 92. Neustadt, L. W. and B. Paiewonsky, 9m mthesimg thimal Controls, in "Proceedings of the 'Second Con- gress of the Internatiomzl Federation Of Automtic Control (IFAC), Basel, 1963', Butterworth, London. Plant, J'. R., “Some Iterative Solutions in Optimal Control," The M,LI.T. Press, Cambridge, Mass., (1968). Pontryagin, L. S,, V. G. Boltyanskii, R. V, Gam- krelidze and E. F. Mischenko, ”The Mathemti-cal Theory of Optiml Processes,” John Wiley and Sons, Inc., New York, (1962). Roxin, E., A Gpometric Intempretation pg Pontmzmgin's Maximum Principle, in ”Nonlinear Differentia Equat one and Nonlinear Mechanics,“ Academic‘ Press, 'New York, Schliching, H., "Boundary Layer Theor ", 6th Ed., McGraw Hill Book Co., New York, (1968). Tapley, B.D‘. and J. M. Lewallen, Comparison of Several Ng‘erical thimizmtion etho s, J. Of Optimization , Theory and App cat ons,‘VO ' , No. 1, (1967), pp‘1—32. Willmore, T. J., o'An Introduction to Differential Geometry,“ The Clarendon Press, Oxford, England, (1959). Zukhovitskiy, S. I. and L. I. Avdeyera, "Linear and (gongg? Programming," W, B. Saunders Co., Philadelphia, 19 . ‘ APPENDICES AI’PIENIJIXC A ANALOG DIAGRAM FOR EXAMPLE PROBLEM 1 T3 T10 T20 T35 T15 T20 .01 11 ' 13 R3 03 2 1’4 R0 A0 P1 D2 02 A M A7: A81 3 .2 both .01 3 M 5 S6 Al ‘1 “ A 6 13 «0 output T16 [b- output indicates an integrator, A indicates an amplifier, indicates a relay; P indicates a potentiometer; ' indicates a multiplier; C indicates a comparitor; indicates a trunk line (to the digital computer); indicates an electronic switch. » A11 amplifiers and integrators consisted of two operational amplifiers, hence both a positive and negative output were available. Trunk 20 was used to switch the system for reverse time integration. . 001 (0'93de 159 APPENDIX B EXAMPLE COMPUTER PROGRAMS 1. GLOBAL OPTIMUM PROCEDURE AND OPTIONS“ // FOR GLOP2 *IOCS(CARD,1““3 PRINTER) *LIST SOURCE PROGRAM *NONPROCESS PROGRAM *ONE WORD INTEGERS DATSW 0 PROVIDES FOR XM=-XM*XM DATSw 1 PROVIDES KICXOUT OF ANY INEFFECTIVE ITERATIONS DATSW 2 PROVIDES DIRECT ACCESS TO DIRECT SEARCH ROUTINE DATSW PROVIDES STEP SIZE INIEPENIENCE 0F ERROR DATSW PROVIDES FOR LESS SUBCURVATURE MOVES, DEL(KK)=1 DATSw 5 DETERMINES THE STARTING VALUE OF KK,NORMAL-10 DATSW 6 BYPASSES THE ORTHOGONALIZATION IMPROVEMENT ROUTINE DATSw 7 CORRECTS FOR ANALOG INTEGRATION INACCURACIES DATSW 8 DETERMINES KM, NORMAL = COS, UP =- -1 DATsw 9 DETERMINES THE CURVATURE CHOICE (SINGLE 0R COMP. ) DATSW 10 STOPS ALL GOAN. PRINTOUTS EXCEPT ERROR PRINTOUT DATSW 11 SKIPS THE RANDOM INPUT AND SEARCH 0ch mm 12 PROVIDES FOR XICIIOUT OF THE RANDOM INITIAL ROUTINE DATSWl PROVIDES FOR PRINTING. OUT TRAJECTORY POINTS DATSWl SETS THE BASIC NMER OF INTEGRATION STEPS AT 10 DATSW 15 IETERMINES-THE INTEGRATION STEP MULTIPLIER SSWTCH 0 PROVIDES FOR TYPEWRITER INPUT AND OUTPUT SSWTCH 1 PROVIDES FOR THE CURVATURE AVERAGE 0F xx AND m SSWTCH 2 PROVIDES FOR ST. IEV. DEPENIENT ESTIMATE OF OURV. SSWTCH PROVIIES A BASIS OF 100 R(T) BOUNDARY PTS, N-25. SSWTCH PROVIDES ANR(T) DATA POINT MULT. 0F 5, N-2. SSWTCH 5 PROVIIES FOR RANDOMLY PICEING POINTS ON BOUND R(T) SSWTCH 6 PROVIDES FOR MULTIPLYING VALUE OF XD BY 2 ITI (FIRST READ) SPECIFIES THE NUMBER OF STANDARD DEV. . POSSIBILITIES OR CORNER CURVATURE ATTEMPTS TOLl (2ND READ) SPECIFIES THE KICKOUT ERROR TO THE DIRECT SEARCH ROUTINE TOL2 ( RD READ) SPECIFIES THE FINAL ERROR TOLERANCE TOL3 ( TH READ) SPECIFIES THE STANDARD IEYIATION MAXIMUM DIMENSION X(0)P P,(0),XI(0),PI(0) XJ(0),PJ(0), 0(0) ,XL(0), x0(0),Pc(0), XS(“, 10),PS(“, 10), XERS(.11) IY=5309 wRITE (1,2180) -' 218“ FORMAT((' ()', ,'() ) )','( )','( )', 0 READ (6,2183) ITI TOL1,T0L2, TOL3 2183 FORMAT (12, 3F10.“) CALL DIGO (20,—1) READ (2, 2)T, N 2 FORMAT (F10. 0 ,11) 00 O OOOOOOOOOOOQOOOQOQOQOQOO 160 “:00 “01 “02 “03 “06 “O7 “08 1000 1117 1001 1002 1003 100“ 1005 1006 1007 1020 1“75 1010 1011 090 091 8 161 FORMAT (8F10. 3) READ (2, 3) (X1(I),I=1,N), (PI(I),I=1,N) WRITE (3,001) FORMAT (' THE FOLLowI'm DATA SWITCHES wERE UP ') D0 ’40le 1:1,16 I"=I-1 1 CALL DATSw(Iww 1w) GO TO (002, “0“),1 WRITE (3,003) IWW' CONTINUE FORMAT (I00) , WRITE (3,006) FORMAT (v THE FOLLOWING SENSE SWITCHES WERE UP ') D0 008 I=1,8 1WW=I-1 CALL SSWTCH (IWW, 1w) GO TO (007, 008), 1W ‘WITE (3, “03)IWW CONTINUE KY=1771 ITRMsO CALL SSWTCH (5, IR) GO TO (1000, 1011), IR READ (6,1117) IRSY .. FORMAT (I0) CALL SSWTCH (3,11) GO TO (1001,1002), II INRS= 100 GO TO 1003 INRS=25 CALL SSWTCH(“, JJ) GO TO (1000, 1005), JJ INRS=1NRS*5 GO TO 1006 INRS=INRS*2 DO 1010 I=1,INRS D0 1007 J=1,N IRSX=IRSY CALL RANDU (IRSX, IRSY, RSY) PI(J)=RSY-.5 CALL SSWTCH (0,100) GO TO (1020,1010), IOO WRITE (1,1“75) FORMAT ('( )( )( )') READ (6, 3) (PI(I),I-1,N) ' CALL GOAN(N,T, XI, PI, X, P ,,X2,P2 XP ,ER) CALL DATSw (111m)' GO TO (090, 091), IRX IZA=1O GO TO 092 IZA = 0 1Y=53“9 DO 882 I=1,N IX=1Y 111‘ I .*-\. 1335;"? )' 3' inn 1:,0118) .’ .1‘ , , Ll.) E ~ 1' ,n: 1 ) .,.'.<"0{5.EOA|) :- «iitrod,c)- .0 l. " o J .2 :Jb—q _ J- rigs: k ‘ no 0" ”Hi I .~“-,-.TA-_1_ x 1. A A” 882 “92 101 11 379 13 1“ “923 1113 1322 12 555 556 101 1“11 1321 1212 2181 2182 218 219 162 CALL RANDU (IX, IY, Y) CALL DATSW (12,1RK) GO TO (9,882L IRX PI(I)=Y‘O WRITE (3,101)(XI(IL I-I,N) FORMAT (l/l0F20 .8) CALL GOAN(N, T, XI, PI, X, P, X2, P2, XP, ER) ERsaER IF(ER-TOL2)5, 5, 379 ERx=ER/X2)-1. IF (ERX) 0923,13,13 DO 10 I=1N PI(I)=-PI(I) CALL GOAN (N, T, XI, PI, X, P ,X2,P2, XP ,ER) ERssER ERX=(ER/X2)-1. IF(ERX) 0923, 8, 8 CALL DATSW(2, IE) GO TO (1113 12L IE HH:S0RT(ERS +. 5’ DO 1322 1= 1, N 0(I)sPI(I) ERDsERS GO TO 1662 PZ=.05 XPI=0. m555 1.1 N XPI=XPI+PI(I)*PI(I) XPI=SQRT(XPI) 556 I=1N PI(I)=PI(1)*5. /XPI PCO=O DO 101 I=1,N KX=KY CALL RANDU (XX, KY, YK) PC(I)=YK905 PCO=PCO+PC(I)*PC(I) PCO=SQRT(PCO) Do 1011 I=1, N PC(I)=PC(1)/PCO LL=0 IF(ERS—TOL1)20,20,1212 PZ=_5.*PZ LL=LL+1 lF(LL—“) 217, 217, 2181 IF(ITKM-IT1)2182, 2182, 218 ITKMsITXMAI GO TO 12 D0 219 I=1,N XL(I)=-1000 xxe-IOO. XXO=-IOO. 216 2160 2162 2161 2163 '180 180 .181 160 162 161 220 102 170 171 1791 179 20 163 GO TO 220 Do 216 I=L N Q(I)sPI(I)+Pz*20(I) CALL SSWTCH (0,100) GO TO (2160, 2162), 100 wRITE (1,1075) READ (6,3) (0(1), I=1.N) CALL GOAN (N, T, XI, .0, XJ, PJ, X2D, P2D, XPD, ERD) IF(ERD-TOL2) L 1,2161 IF(ERD—TOLI) 1662, 1662, 2163 CALL DATSw (6 I0) GO TO (160,180), IO IF(.95-ERD7ERs)16O,l60, 180 D0 181 I=1,N PI(I)=Q(I) x(I)=XJ(I) P(I)=PJ(I) ER=ERD ERS=ER X2=X2D P2=P2D GO TO 217 xxeo. 0 DO 161 I=L N IF (ABS(XJ(I)-X(I))-.OOOOI) 1212,1212, 162 XL(I)=((PJ(I)/P2D)- (P(I)/P2))/(X(I)-XJ(I)) XK=H+XL(I) xxsXX/N WRITE (3,102) (nu), 14,15), Xx C1=XP/;P2*x2) C2=XPD (22D*X2D) XKQ=(Cl-CZ)/(X2D-X2)-C1/XZ XXAV=(XX0XXN)/2. WRITE (3,102) XEQ, XKAV FORMAT (7F15. 5) CALL SSWTCH (2 IV) GO TO (170,179 ,IV AAer VAR=0 DO 171 1:1, N VARuXL(I)*XL(I)+VAR AAXsAAX+ABs(XL(I)) VAR=VAB/N-XK#XK STDV=SQRT(VAR) AAstAK/N TEST=STDV7AAX WRITE (3.102) VAR ,STDV,AAX,TEST ITXM:ITXM+1 IF(TEST-TOL3)179,179 1791 IF (ITKMAITI) 12,12,179 IF(ERD—ER) 30,30, 20 DO 21 I=1N Q(I)=PI(I)N 21 30 3111 3112 3113 3211 1201 1200 1202 1203 120k 32 3291 3292 329 329 1663 166b 1662 1665 1666 1661 166 1671 1672 C 1673 1670 1675 16k XJ(I)=X(I) FJ(I):P(I) ‘ERDéER X2D=X2 Eassnn :F2D:F2 CALL DATsw (5,1MX) GO TO (3111,3112), IMX KK:-3 . GO To 3113 KK=O HH:SQRT(ERS)+.5 IF (ERS-TOLZ) 1,1,3211 CALL DATsw (3,18)‘ G0 T0 (120L 1200), IS XD=. 2 GO To 1202 XD=ENSIX2D CALL SSWTCH (6,1DM) GO T0 (1203,1201), 1111 XD=XD"2. WRITE (3,102) XD IDITE (3,102)(XL(I),I=1,N),XK IF(ERS-TOLI) 1662 1662,3291 - CALL DATSW'(u, IXX GO TO (3293,3292), IXX KK#KK+1 G0 T0 32911 KK:KK+2 IF(KK-2)1661, 1661,1662 HH=HH*. 5 IF(HH+.005)1 166b,166b WRITE (3,102 HH CALL DATsw (1 MM) GO To (1,1662, CALL EXPLR(N T, XI Q, XJ ,PJ ,,X2D P2D, XPD, ERD, HH) IF (END-T0L231, L 166 IF (ERD—ERS)1666, 1663,1663 ERS=EBD GO To 1661; CALL DATsw (1,11) G0 T0 (L 166), 11 CALL DATsw(8, IXM) GO TC (1671, i672), IXX ”“1. GO TO 1673 XM:(DDs/X2D)~1. XM=0 CALL DATsw (0,1IM) G0 T0 (167u,1675), IYM XM=—XM¥XM WRITE (3,102) XM 165 DX=- XD*(XJ(I)- (X2D*PJ(I)*XM)/P2D) CALL SSWTCH (L IXAv) GO TO (1210,1213), IXAv 1210 F(I):PJ(I)+F2D*DX*XXAV GO To 35 1213 CALL DATSW (9,1X) GO To (1210,1211L IX C F(I)=FJ(I) C P(I)=PJ(I)-P2D*DX*XK 1210 F(I)=FJ(I)4?2D*DX*XL(I) G0 T0 35 1211 P(I)=PJ(I)-P2D*11X*IKQ C X(I)=XJ(I) 35 X(I)=XJ(I)+DX CALL DIGO (2L 1) WRITE (3,102) (X(IL I=1,N) WRITE (3,102) (P(I),1=1,N) CALLGCAN (N,,,TXP X,F, X2,P2,,XP ER) CALL DATSW’(7, JL) G0 T0 (891, 886), JL 891 CALL GOAN (N, T, XJ, PJ, XC, PC, X20, 220, XPC, ERC) D0 888 I=1N X0(I)=XC(I)-XI(1) 888 PC(I)=PC(I)-Q(I) WRITE (3,102) (XC(1), I=1,N) WRITE (3:102) (FC(1L I=1,N) D0 889 I=L N X(I)=X(I)-XC(I) 889 PI(I)=P(I)-PC(I) WRITE (3,102) (X(l) 1:1,N) WRITE (3,102) (FI(1,I=1,N) ~Do 890 l=lN 89o FI(I)=:F1(I)N GO TO 887 886 D0 885 I=L N 885 RI(I)=T(I) 887 CALL DIGO (2L -1) . _, CALL CCAN (N T, XI, PI, X, P ,X2,P2,XF,ER) IF (ER.ERs) 0L 01, 01 XD=XD*.2 012 WRITE (3,102)XD GO T0 32 02 ERs=ER WRITE (3,102) ERs ,ERs,ERs, ERs ITXM=0 IF (ERs-ToLz) 5, 5.12 D0 6 I= L N Q(I)=PI(I) XJ(I)=X(I) X2D=X2 D0 7 1= 1 N PS(I,IZA1=Q(I) XS(I, IZA)=XJ(I) XERs(IZA)=x2D IZA=IZA+1 IF(IZA~11)8,997,997 \IHQUI 1 2114 1213 1210 1211 35 891 888 889 890 886 885 887 01 012 \li-‘O\Ut 165 DX=~XD*(XJ(I)- (xszJ(I)*XM)/P20) CALL SSWTCH (L IXAv) G0 T0 (1210,1213), IXAv P(I)=PJ(I)—P2D*DX*XZKLV GO TO 35 CALL DATSW (9,1K) GO T0 (1210,1211), IK P(I)=PJ(I) P(I)=PJ(I)-P2D*DX*XK P(I)=PJ(I)-P2D*DX*XL(I) 00 T0 35 P(I)=PJ(I)-P2D*DX*IKQ X(I)=XJ(I) X(I)=XJ(I)+DX CALL 0100 (20,1) WRITE (3,102) (X(I),I=1,N) WRITE (3I102) (P(I):I=1:N) CALLGOAN (N,,,TXP X,,P X2,,P2XP, ER) CALL DATsw (7, JLI GO To (891, 886), JL CALL GOAN (N, T, XJ, PJ, XL PC, X20, P2C, XRL ERC) no 888 I=1N XC(I)=XC(I$-XI(1) PC(I)=PC(I)-Q(I) WRITE (3,102) (XC(IL I=1,N) WRITE (3,102) (PC(I), I=1IN) DO 889 I=1,N x(I)=X(I)-XC(I) PI(I):P(I)-PC(I) WRITE (3,102) (X( PI ) I=1,N) WRITE (3I102) ( I I ( ,I=1,N) ’DO 890 I=1N PI(I)=:PI(I$N 00 T0 887 no 885 I=L N PI(I)=P(I) CALL 0100 (20, -1) CALL GOAN (N T, XI, PI, X, r, X2, r2, 3?, ER) IF (ER.ERs) 02I 0L XD=XD*.2 WRITE (3,102)XD GO TO 32 ERs=ER WRITE (3,102) ERs ,ERs,ERs, ERs ITXM= 0 IF (ERS-TOLZ) 5, 5.12 no 6 I=1, N Q(I)=PI(I) XJ(I)=X(I) X2D=X2 DO 7 I=1N PS(I, IZA3=Q(I) XS(I, IZA)=XJ(I) XERs(IzA)=X20 IZA=IZA+1 IF(IZA911)899979997 166 99? CALL DATS"(11, IRX) GO To (u, 9), 9 no 999 J=1,1oIRX 999 ‘WRITE £3,102)XEBS(J), ,(xs(1, J), I=1,N), (PS(I, J), I-1,N) GO TO 99 CALL EXIT END *As previously mentioned, several options or alternatives are provided as a result of the CALL DATSW and.CALL SSWTCH subrouiines. 'In.addition, unnecessary printout routines are included. 2. DIGITAL INTEGRATION AND ERROR FUNCTION EVALUATION ROUTINE // FOR *LIST SOURCE PROGRAM *ONE WORD INTEGERS *NONPROCESS PROGRAM SURRCUTINE GoAN (N, T, n, PI, 11,? ,,XZ,P2 11? ,ER) EXTERNAL F, F0 DINEN310N(§1(9),,X(u), R(u), x1(u), 1(9), nr(9), Aux (1o ,9), PRMT COMMON 1TRNK(8),0LD(9),IYZ(5)' NAX=2*N no 1 1:1, N Y(2*I-l)=XI(I) 1 Y(2*I)éPI(I) ITRNK(3)=O E=.5/N CALL RATsw(1h,11) GO To (591,592),11 591 HT=10. GO To 593 592 HT=50. 593 CALL EATsw (15 JJ) GO To (59h, 595$.JJ 59h HT=HT*2. 595 HT=HT 596 IF(ITRNK(1)) 3, 3, u 3 PBMT(1)=O TEaT XYZ(1)=TF PRMT(2)=T PRMT(3)=T/HT GO TO 5 u PRMT (l)=T TF=O XYZ(1)=TF PBMT(2)=0 PRMT(3)=—T/HT NU‘ 525 528 526 527 77 30 308 110 112 111 18 10 20 31 32’ 38 167 no 2 I-I,NAx nW(I):E ITRNK(2)=O CALL HPCG (RRMT, Y, DY, 2*N, IHLF, F, no ,Anx) IF(PRMT(5)) )525 7 525 1R(ITRNK(3)) 7, 52é,7 no 526 KJ=1,9 Y(KJ)=OLD(KJ) PRMT(1)=PRMT(1)-PRMT(3) PRME(3)=PRMT(3)/10. RRMT (2):?RMT(1)+IO.%PRMT(3) no 527 KJ=1,NAX . DY(KJ)=E ITRNK(2)=1 CALL HPOG (PRMT, Y, ET, 2*N, IHLF, F, no, AUX) RRMT(1):PRMT(2) RRMTB):PRMT(3)*10. IF(ITRNK(3)) 7, 77, 7 PBMT(2)=TF GO TO 5 DO 6 I=l, N X(I):Y(2*I-1) P(I)=Y(2*I) n2=o. no 30 I=1 N 22:P2+P(I)*P(I) P2=SQRT(P2) no 308 1:1, N P(I)=T(I)*5 /P2 R2=5. CALL nATsw'(lo, IPX) I GO TO (20,110), IPX WRITE (3,112)T FORMAT (F15. 5) FORMAT (' THE CURRENT ERROR 1: ',F15. 5,' NORM x = ', F15 5.‘ XP= ',F10 5) no 18 L=1, N WRITE (3,10) x1(L), PI(L), P(LL X(L) FORMAT (6F20. 5) x2=o. o . no 31 I=l N x2=xz+X(I$*X(I) X2=SQRT(XZ) xn=o.‘ no 32 I=l N xr=xn+x(1)*r(1) ER=(XP/P2)+x2 an=ER/x2—1. ‘ WRITE (3,111) ER, 12, an CALL SSWTCH (0,100) GO TO (38, no), 100 no 39 J1=i, N 168 39 WRITE (1,10) XI(JI) ,PI(JI) ,P(JI), x(JI) WRITE (1, ’111) ER, x2, XPN no RETURN END 3. DIGITAL INTEGRATION FORWARD-REVERSE TIME SURROUTINE" // POR *LIST SOURCE PROGRAM *NONPROCESS PROGRAM *ONE WORn INTEGERS SURROUTINE mm (M, N) COMMON ITRNK(8), OLn(9) ITRNK (M-19)=N ' RETURN ENn *Thil aubrOutine supplements the digital integration sub- routine to replace a hybrid subroutine called in.the main program. ’4. HYBRID INTEGRATION AND ERROR FUNCTION EVALUATION ROUTINE // FOR *NONPROCESS PROGRAM *LIST SOURCE PROGRAM *ONE WORn INTEGERS SURROUTINE GOAN (N T, XI PI, 1, P x2, P2 KP ,ER) DIMENSION PI(u) ,X(fi),P(fiL x1(u$,JI(u3, K1(u) ,KJ(h),JJ(h) CALL LOGEX(1) no 5 J=1, N u JI(J)=XI(J)*3276. 7 5 CALL ANOUT (31+J, JI(J)) XPI=O. W555 I=lg N 555 XPI=XPI+PI(I)*PI(I) XPI=SQRT(XPI) no 556 121 N 556 PI(I):PI(I$*IO./XPI no 8 x:1, N JJ(K):PI(K)*3276. 7 8 CALL ANOUT (33+E:,JJ(R)) NT=T*h0. +30. CALL LOAn CALL DELAY (8000) CALL RUN no 7&2 I=1, 10 7h2 CALL DELAY (NT) CALL STOP! CALL AINP (12 ,N,KI(1), KI(2)) CALL AINP (1'4 ,N,KJ(1), KJ(2)) no 19 L=1,N 19 30 308 110 111 112 18 20 31 32 5. 169 x(L)=-K1(L)/327.67 P(L)=-KJ(L)/327. 67 CALL AINP (16,1 TK#-KL/327. 67 P2=0. D030 I=1 N ' P2#P2+P(I$*P(I) F2=s0RT(F2) D0 308 I=1, N P(I)=P(I)*10. [F2 F2=10. CALL DATSW (10,1Fx) GO To (20,110), EFx WRITE (3,112)TK FORMAT (' THE CURRENT ERROR IS ',F15.5) FORMAT (F15. 5) - DO 18 L=1, N ‘WBITE (3,10) x1(L), P1(L), P(L), x(L) FORMAT (6F20. 5) CALL LOAD X2=0. DO 31 I=1 N x2=x2+X(1$*X(I) x2=sORT(x2) xF=0.0 D0 32 I=1N XPBXP+X(I‘*P(I) ER=(XP/P2)+X2 WRITE (3,111) ER RETURN END DIRECT SEARCH SUBROUTINE // FOR *LIST SOURCE FROGRAM *ONE WORD INTEGERS . *NONPROCESS PROGRAM 205 210 215 SUBROUTINE EXPLR(N, Tx ,FI, x,F DIMENSION x1(u) ,F1(u$ ”x(h F(u3 Do 260 I=1N FI(I)¢FI(I$-HR . 011.1. GOAN(N, T, H F1, 1:, F, x2, F2, xF ,Es) IF (ES-EB) 205,210, 210 ET=ER EBFEB IF( .80—ER/ET)260, 260, 265 FI(i)=FI(I)+2.*HH CALL GOAN (N, T, XI ,FI, x, P ,,12,F2 1F ,Es) IF(Es-ER) 215,220, 220 . ET=ER ERES x2,F2, 1F, ER, HR) 170 IF(.95—ER/ET)260, 260, 265 220 PI(I)#PI(I)-HH 260 CONTINUE 265 NRITE (3,270) (PI(II) ,II-1,N) 270 FORMAT (9F10. 5) 99 RETURN END 6. DIGITAL INTEGRATION OUTPUT AND CONTROL DISCONTINUITY SENSOR SURROUTINE //'FOR *NONPROCESS PROGRAM *ONE WORD INTEGERS *LIST SOURCE PROGRAM 100 101 111 113 10h 103 10 1 2 22 109 108 107 110 7. SUEROUTINE FO(T, T ,DY,IHLF, NN,PRMT) DIMENSION 1(9) PRMTCS) COMMON ITRNE(85, OLn(9‘), XYZ(5) IF(T-PRMT(1)) 101,100,101 SIGN=Y(9) SIGNV=Y(8) GO TO 103 IF(ITRNK(2)) 111,111,103 IF(Y(9)*SIGN)IOu 1131 IF(Y(8)*SIGNV(1M 03,103) FRMT(5)=1. ‘ PRMT(1)=T CALL DATSw (13, Ex) GO TO (1, 2), ER FORMAT(12F10. 3) ‘wRITE E3Y10)8 T ,PRMT(u), (Y(I), I=1, NN, 2), (Y(I), 1:2, NN, 2), Y 9 IF(PRMT(53)108,22,108 DO 109 I=19 OLD(I)=Y(I$9 IF(ARS(XY2(1)-T)—. 0001)107,107, 110 ITRNK(3)=1 PBMT(5)=1. RETURN END RANDOM NUMBER GENERATOR SUEROUTINE // FOR *NONPROCESS PROGRAM *ONE WORD INTEGERS O\UK SUEROUTINE RANDU(Ix, IY, Y) IY=IX*899 IF(IY)5, 6, 6 IY=IY+32767+1 Y=IY Y=Y/32767 RETURN END 8. 171 BASIC DIGITAL INTEGRATION SUEROUTINE“ // FOR *LIST SOURCE PROGRAM *ONE WORD INTEGERS 22 23 2h 25 2? SUBROUTINE HPCG (PRMT, Y IERL NDIM IHLF, FCT, OUTP ,AUX) DIMENSION PRMT(5), 1(2003, DERY(200LAUX(10, 200) IHLF=O x:PRMT(1) H=PRMT(3) PRMT(u)=o, PRMT(5)=O'. DO 1 I=1 NDIM AUX(10, 15=0. AUX(9, I)=IERY(I) IF (H*(PRMT(2)-X)) 3, 2, 5 IHLF=12 GO TO 1; IHLF = 13 RETURN DO 10 N=1,3 CALL FCT (x, Y, DERI) IF (N—l) 1L 11,12 CALL 0UTP(x, L DERY ,IRLF,NDIM, PRMT) IF(PRMT(5)) ’11; 6,14 Do 9 I=1N,DIM AUX(N, I)=Y(I) AUX(N+h, I)=DERY(I) D0 101 I=L NDIM F(I)=AUX(N, I)+H*AUM(N+u, I) x=x+n CALL FCT (x, L IERY) DO 102 1:1, NDIM Y(I)=AUX(N, I) +. 5*H*(AUX(N+14, I)+IERY(I)) CONTINUE. N31 CALL FCT (x, Y ,IEBY) x=PRMT(1) DO 22 I=L NDIM AUX(8,I)=1ERY(I) Y(I)=AUX(1, I)+H*(. 375*A,Ux(5 I)+.791667*AUX(6, I) -.2083533*AOX(7, I )tOll166667" IERY(1) ) x=x+H N=N+1 CALL FCT(x, L IERY) CALL OUTP(x, L mRL IHLF, NDIM, PRMT) IF(PRMT(5)) ”u ,200 IF(N;h)25, 20k, 206 D0 26 1:1 NDIM AUX(N+LL I3: -IERY(I) IF(NA3)é7, 29, 20h DO 28 I=1, NDIM DELT=AUx(6, I)+AUX(6, I) 28 29 30 200 383 205 207 208 209 215 210 212 213 21b 172 IELT=IELT+IELT _ Y(I)=LUX(1,I)+.3333333*H*(AUX(5,11+DELT+AUX(79I1) GO TO 23 DO 30 I=L NDIM DELTaAUX(6, I)+AUX(7,I) DELT =IELT+IELT+1ELT Y(I)=AUX(1,I)+.375*H*(AUX(5,I)+DELT+AUX(8,I)) GO TO 23 Do 203 N=2, 8 DO 203 I=1,NDIM AUX(N-L I)=AUX(N,I) DO 205 I=L NDIM AUX(u,I)=F(I) AUX(8, I)= DERY(I) =x+R DO 207 I=L NDIM DELT=AUX(1,I)+1.333333*R*(AUX(8,1)+AUx(8,1)-AUx(7,1) +AUX(6, I)+AUX(6, 1)) 2(1)=DELT-.9256198*AUX(10, I) AUX(10, I)=DELT CALL FCT(x, Y, DERY) DO 208 I=1, NDIM DELT=.125*(9.*AUX(u, I)—AUX(2, I)+3.*H*(DER2(I)+AUX(8, I) +AUX(8, I)-AUX(7,I))) AUX(10, I)=AUX(10, I)-UELT 2(1)=DELT+.07O38017*AUX(10, I) DELT=). DO 209 I=1, NDIM RELT=DELT+AUX(9, I)*ABS(LDX(10, 1)) IF (PRMT(u)-DELT) 215,210, 210 PRMT(u) = DELT CALL FCT(x, L DERY) CALL OUTP(L 2,11ERLIRLLNDIM,FRMT) IF(PRMT(5)) 212,213,212 RETURN IF (H*XAPRMT(2))) 21h, 212, 212 IF (ABS(XAPRMT(2))-.1*ABS(H)) 212,200,200 END I"Thus: integration.routine is part of the Scientific Sub- 9. routine Set for the IBNZSyStem/36 O. FUNCTION EVALUATION SURROUTINE FOR.ES-2 ' I l // FOR *LIST SOURCE PROGRAM *NONPROCESS PROGRAM *ONE WORD INTEGERS SURROUTINE F(T, L DR) DIMENSION 2(9), NIB) %F{Y(6)), 2, - O\U\ F'uJN 173 G0 T03 U=ABS(Y Y(6)) Y(6) gin” )) 5., GO T06 Y=ABS(Y(#))/Y(U) Y(9)=U Y(8)=v DY(1)=Y(3) DY(3)=Y(5 )+v DY(5)=—Y(5)*Y(1)+Y(3)*Y(3)+U DY(2)=Y(5)*Y(6 ) mm )=-Y(2)-2.*Y(3)*Y(6) m(6)=-Y(h)+Y(1)*Y(6) RETURN END APPENDIX C INPUT DATA SETS 1. lgput Data Set 1 11210.21; _fl (1%)" Me; 41523. .2223. 1 (-109 ‘5) ( O, 5) .9 ( “59‘5) ( 59-5) 2 ( 2,-10) ( -5, -E) 10 (-10,-5) ( 5,-1) 3 ( 2, 6) ( 5, ) 11 ( 5, 5) ( 5,-1; h ( 6, 2) ( 10, -3) 12 ( -5, 5) ( 5. 5» 5 ( -7. -2) ( -2, -2) 13 (210,-5) ( 5, 3) 6 ( '19 9) ( 39 ‘u) I” ( 59 5) (“39‘1) 7 ( 5. 0) ( 5, -E) 15 ( 5,-5) ( 5. 5) 8 ( -5, 0) ( -3, ) 2. I ut Data Set 2: x0 = (-10,-5)T, po 19 randomly genp eratefi and normalized to 10. Number (PO)T NuMber (P9)T 1 ( 9.879, 1.551) 6 (44.919, 8.706 2 ( 9.997.-0.25'4) 7 ( 9.82 ,-1.871 a ( 3.671,-9.302) e 8 ( 9.2 ,-3.766) {-6.959.-7.181) 9 ( 9.970,-0.769) 5 (-8.I71, 5.76h) 10 (-0.235,-9.997) 3- .1222£_2232J§2142 Number (xb)T (PO)T 1 (“39‘29'1) ( 19‘29-3) 2 ( 2,-1, I) (—1, 0, 0) a ( 2"‘2, 2) ( 1, 1,-1) ( 2,12, 0) ( 1,-I,-g) 5 ( 2,-2,-2) (-2’-2, ) 6 ( , o, 2) (—E,—h,—2) 7 ( 2, 0, 0) ( 9'39 5) 8 ( 29 0"2) ( 79'79'7) 9 ( 29 29 2) ('39 1, 2) 10 ( 2, 2, O) ( 5, 2, 5) h. Lgput Data Set h: 10 = (-3,-2,-1)T, p0 is randmmfly generated: Number (PO)T Number ‘ (P0)T 1 ‘- ( .251“, .039“, ,ou733) '4 ( .1732,-;2'.497, .01-120) 2 (_.0120, .1630,-.h131) 5 ( .uOO3,-.O762, .0392) 3 ("ou‘659-0u6079-02’4‘56) 6 ('01785, .u760,-.0366) 17h 1. 2. 3; 9. 10. 11. APPENDIX D POSSIBLE EXTENSIONS AND FUTURE INVESTIGATIONS Further examination of time optiml central problem. Extension to other Optiml control problems such as min-,- immn fuel, convex function minimum error regulator,etc. Determination of computational results for the special cases given in Chapter it. Another approach to the solution of MP based on LOP, other than simply the generation of random initial adjoints. DevelOpment of higher order systems with nonconvex reachable sets and mny local optim. Application to parameter Optimization problems. Comparison of GOP and LOP with other nonlinear methods. Further investigations of the nature of the reachable sets as they relateto the nonlinear problem involved: a. Classification of reachable sets in some manner. b. Transformation of reachable sets as a result of the transformation of the origin. Effect of singularities and oontrollabilityon the general problems considered in this thesis. Relationship of the path composed of lines of curvature to the path determined by gradient methods. Use of geodesics instead of lines of curvature as the means of defining the path on the boundary of the reachable set. 175