THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING Information and Control Engineering Program COMPUTATION OF OPTIMAL CONTROLS BY QUADRATIC PROGRAMMING ON CONVEX REACHABLE SETS Robert Ortha, Barr, Jr. ORA Projects 06999 and 08180 Supported by: AIR FORCE OFFICE OF SCIENTIFIC RESEARCH OFFICE OF AEROSPACE RESEARCH ARLINGTON, VIRGINIA Grant No. AF-AFOSR-814-65 and AF-AFOSR-814-66 Administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBIOR September 1966

via \N, 1-INA YV\ L 1; -ooo

This report was also a dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan 1966

Abstract COMPUTATION OF OPTIMAL CONTROLS BY QUADRATIC PROGRAMMING ON CONVEX REACHABLE SETS by Robert Ortha Barr, Jr. This thesis develops iterative procedures applicable to a wide variety of optimal control problems. The computational algorithms are quite different from these which have been described previously and apply under more general conditions, Convergence of the procedures is proved and numerical results for example problems are given. Essential to all the procedures is a method for solving the basic n problem BP: given K a compact, convex set in E; find a point z* e K such that |z2 = mi I z 12z ( 1. denotes Euclidean norm). The constraint set K in this quadratic programming problem, in contrast to the quadratic programming problems usually treated in the literature, need not be specified by some set of functional inequalities. It is required only that there be a known contact function of K, i.e., a function s(.) from En to K such that y. s(y) = mKx y' z for y / 0. The method given by E. G. Gilbert for solving BP using only a contact function of K is extended so that on each iteration a quadratic minimization problem is solved on a convex polyhedron instead of on a line segment. As with Gilbert's method, computable error bounds are available. Techniques for readily solving the minimization problem on a convex polyhedron are discussed and extensive computational results for a rather general example using both Gilbert's method and its extension are presented. These results indicate that the extended procedure has a much improved rate of convergence. Furthermore, it is proved that if K is a convex polyhedron and the range of s(y), y E En, is a finite set of points, then the extended procedure exhibits finite convergence.

Certain optimal control problems are related directly to BP. However, it is possible to state an abstract problem which has application to a much broader class of optimal control problems. This problem, called the general problem GP, is: given a compact interval Q = [0o, a], > 0, a family of sets K(w) in En which are compact, convex, and continuous on Q2, and the fact that there exists o" E 02 such that 0 a K(w), 0 w < w-', and 0 E K(wo); find o:-. An iterative procedure for solving GP, which on each iteration involves the minimization problem BP, is described and shown to converge. Numerical results for this procedure applied to a minimum fuel example are given. Optimal control problems are considered in which there are. i) a system differential equation of the form x(t) - A(t) x(t) + f(u(t), t), x(0), where x is the m-dimensional state vector, x(0) is the initial state, u(.) is an r-dimensional control vector, admissible if measurable with range in a compact set U, the matrix function A(t) and vector function f(u, t) are continuous in their arguments; ii) compact, convex target sets W(t) in EM which are continuous in t; iii) a cost functional t Jt(u) = [a(cr) x () + f0(u(o-),s)]dr- + h~(x (t)), where x (t) is the solution of the system equation corresponding to the u control u( ), h~( ) is a convex function from Em to E1, the vector function a(t) and scalar function fo(u, t) are continuous in their arguments; iv) an optimization objective corresponding to fixed or free terminal time. For a fixed terminal time T >0 the objective is: find an admissible control u*(-) such that x,(T) E W(T) and the cost for u'*() at T < the cost for u(e) at T where u( ) is any admissible control for which x (T) e W(T); if the terminal time is free: find an admissible control u*(. ) and an optimal time t* such that x,(t*:) E W(t*) and the cost for

u:( ) at t<' the cost for u(-) at t where u'(.) is any admissible control for which x,( t) E W( t), any t. The iterative procedures apply to many u problems of this class and to other optimal control problems as well. The set of all solutions of the system equation at a particular time t which are generated by admissible controls is the reachable set. This set is (1) compact, (2) convex, (3) continuous in t, and (4) has a contact function which can be evaluated. These four properties are the essential features which permit application of the iterative procedures to optimal control problems.

Acknowledgements The author wishes to express his gratitude to each committee member, especially to Dr. E. G. Gilbert who provided thoughtful criticisms, encouragement, and suggestions throughout the course of this research. The helpful comments of Dr. P. Falb are also gratefully acknowledged. The author is indebted to Mr. John A. Parsons for his excellent work in typing the manuscript. Sponsorship for the research was provided by the Air Force Office of Scientific Research, Office of Aerospace Research, United States Air Force, under AFOSR grant numbers AF-AFOSR-814-65 and AF-AFOSR-814- 66.

Table of Contents Page Acknowledgements ii List of Tables v List of Figures ix Chapter 1 Introduction 1 Chapter 2 The Basic Iterative Procedure BIP 4 2. 1 Preliminary Comments on Notation 4 2. 2 Contact Function; the Basic Problem BP 5 2.3 The Basic Iterative Procedure BIP 6 2.4 Convergence Theorem for BIP 9 2.5 Nature of Convergence; Examples 13 2.6 Numerical Results for BIP 18 Chapter 3 The Improved Iterative Procedure II P 30 3. 1 Convex Sets and Polyhedra 30 3.2 The Subproblem 1 34 3.3 The Improved Iterative Procedure IIP 35 3.4 Convergence Theorem for II P 37 3.5 Selection Rules for IIP 39 3.6 Solution of Subproblem 1 48 3.7 A Finite Convergence Theorem 52 3.8 Numerical Results for IIP 55 Chapter 4 The General Iterative Procedure GIP 80 4.1 The General Problem GP 80 4. 2 The General Iterative Procedure GIP 82 4.3 Convergence Theorem for GIP 86 Chapter 5 Application of BIP, II P, and GIP to Optimal Control 90 Problems 5. 1 Some Optimal Control Problems Solvable by GIP 90 iii

Page 5. 2 The Reachable Set; Determination of a Contact Function 93 5. 3 Solution of Problems 1 and 2 by BIP or IIP Alone 99 5.4 The Sets K(w) = X(X) - Y(w) 102 5. 5 Solution of Problems 1, 2, and 3 by GIP 103 5.6 Solution of Problems 4, 5, and 6 by GIP 105 5.7 Numerical Results for GIP Applied to a Minimum-Fuel 108 Example List of Symbols 128 Bibliography 133 iv

List of Tables Table Page 2. 6. 1 Number of iterations to satisfy IZk | Iz:r I; 23 n = 2, z0 = (6, 2), BIP. 2.6.2 Number of iterations to satisfy Izk - z': $ s E; 23 n = 2, z0 = (6, 2), BIP. 2.6. 3 Number of iterations to satisfy IZ': I zk | Y(zk) 24 n = 2, z0 = (6, 2), BIP. 2. 6. 4 Number of iterations to satisfy Js(-zk) - z' I<; 24 n = 2, zo = (6, 2), BIP 2. 6. 5 Number of iterations to satisfy zkl- z': $ e; 25 n = 2, 2 = 100, BIP. 2. 6. 6 Number of iterations to satisfy I k |- |z': e; 28 n = 3, z0 = (6, 2, 2), BIP. 2. 6. 7 Number of iterations to satisfy Izk - z:-' I e; 28 n = 3, zo = (6, 2, 2), BIP. 2. 6. 8 Number of iterations to satisfy j z:' - lzk y( zk) e 29 n = 3, z0 = (6, 2, 2), BIP. 2. 6. 9 Number of iterations to satisfy Is(-zk) - z;: I; 29 n = 3, z0 = (6, 2, 2), BIP. 3. 8. 1 Number of iterations to satisfy zk- Iz': J E; 61 n = 2, p = 2, z0 = (6, 2), IIP Selection Rule A. 3. 8. 2 Number of iterations to satisfy I zk - z: ]< $ E; 61 n = 2, p = 2, z0 = (6, 2), IIP Selection Rule A. 3. 8. 3 Number of iterations to satisfy I - I Zzk I Y(z k) E; 62 n = 2, p = 2, z0 = (6, 2), IIP Selection Rule A.

Table Page 3. 8. 4 Number of iterations to satisfy I s(-zk) - z* E; 62 n = 2, p = 2, z0 = (6, 2), II P Selection Rule A. 3. 8. 5 Number of iterations to satisfy I - e; 63 n = 2, p = 2, X2 = 100, IIP Selection Rule A. 3. 8. 6 Number of iterations to satisfy I k I - lz' I s E 66 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 7 Number of iterations to satisfy |zk - zI< <e; 66 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 8 Number of iterations to satisfy z' -J zk JY( zk) - E; 67 n = 3, p = 3, zO = (6, 2, 2), IIP Selection Rule A. 3. 8. 9 Number of iterations to satisfy Is(-zk) - zIj: <; 67 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 10 Number of iterations to satisfy IZkI lIzk Y(zk) - E; 68 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 11 Number of iterations to satisfy l/1 -y(Zk) ZI kl 68 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 12 Number of iterations to satisfy I s(-zk) I - I I E; 69 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 13 Number of iterations to satisfy I Zk l- z' I 69 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. 3. 8. 14 Number of iterations to satisfy Izk I -lz I I < 70 n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule B. 3. 8. 15 Number of iterations to satisfy I Zk I- zI I; 71 n = 2, p = 1, z0 = (6, 2), IIP Selection Rule A. 3. 8. 16 Number of iterations to satisfy I Zk - I zI E; 71 n = 2, z0 = (6, 2), IIP Selection Rule A and BIP. vi

Table Page 3. 8. 17 Number of iterations to satisfy I zkl- z* < E 72 n = 3, ZO = (5, 4, 2), X2 = 1000, X3 = 100, IIP Selection Rule A and BIP. 3. 8. 18 Number of iterations to satisfy Izk -z' | E 73 n = 4, p = 4, z0 = (6, 2, 2, 1), II P Selection Rule A. 3. 8. 19 Number of iterations to satisfy I zk - z E; 74 n = 4, p = 4, zo = (6, 2, 2, 1), IIP Selection Rule A. 3. 8. 20 Number of iterations to satisfy I Zk I- z' |l E; 7 5 n = 4, z0 = (6, 1, 2, 2), X2 = 1000, X3 = 500, X4 = 100, II P Selection Rule A. 3. 8. 21 Number of iterations to satisfy Izk l- z*j 76 n = 5, p = 5, z0 = (5, 3, 1, 1. 8, 2. 6), IIP Selection Rule A. 3. 8. 22 Number of iterations to satisfy zk - z'l: E; 76 n = 5, p = 5, z0 (5, 3, 1, 1. 8, 2. 6), IIP Selection Rule A. 3. 8. 23 Number of iterations to satisfy error criteria; 77 n = 6, p = 6, z0 = (4, 3, 2. 6, 2.6, 1. 8,1. 8), IIP Selection Rule A. 3. 8. 24 Number of iterations to satisfy zkl- zl I E; 7 8 II P Selection Rule C. 3. 8. 25 Number of iterations to satisfy I Zk - z'l -; 79 II P Selection Rule C. 5.7. 1 The first j for which 6 s E and p s E; x(O) = (2, -1). 114 J J 5. 7. 2 The first j for which 6.j E and pj. E; x(O) = (2, 0). 114 J J 5. 7. 3 The first j for which 6. < E and p.' E; x(O) = (7, -3); J J I I P Selection Rule A with p = 3 used on each iteration of GIP. vii

Table Page 5. 7.4 The first j for which 6. < E and pj - E; x(0) = (7. 5, -3); 115 J j II P Selection Rule A with p = 3 used on each iteration of GIP. 5. 7. 5 The functions x. and (i) for x(0) = (2, -1): IIP Selection 116 Rule A with p = 3 used on each iteration of GIP. 5. 7. 6 The functions x. and 2(i) for x(0) = (2, -1); BIP used on 117 1 each iteration of GIP. 5. 7. 7 The functions w. and 2(i) for x(0) = (2, 0); IIP Selection 118 Rule A with p = 3 used on each iteration of GIP. 5. 7. 8 The functions w. and 2(i) for x(0) = (2, 0); BIP used on 119 each iteration of GIP. 5. 7. 9 The functions wi and 2(i) for x(0) = (7, -3); IIP Selection 120 Rule A with p = 3 used on each iteration of GIP. 5.7. 10 The functions x. and B(i) for x(0) = (7. 5, -3); IIP 121 Selection Rule A with p = 3 used on each iteration of GIP.

List of Figures Figure Page 2. 3. 1 Geometric significance of p(') and y(). 7 2. 3.2 Geometry of BIP for 6 = 1. 9 2.6.1 IZkl- IIZ.I for n = 2, zo = (6, 2): A) X2 = 10, B) X = 100, 21 C) X2 = 200, D) X2 = 500, E) X2 = 1000; BIP. 2.6. 2 P(zk) for n = 2, z0 = (6, 2), X3 = 100; BIP. 22 2.6.3 Results for n = 3, z = (6, 2, 2), X2 = 100, X3 = 10: 26 max A) z-z*| i<k IzilY(zi), B)!Zkl-lZ* ], C) zk - Z*l; BIP. For k < 14, Izk - z:l Izkl-: I. 2.6.4 Results for n = 3, z0 = (6, 2, 2), X2 = 100, 3 10: 27 A) Iz*l-lzk Y(Zk), B) Is(-zk) - z'; BIP. 3.5. 1 Example 3 of Section 2. 5, n = 2, v = 1, k2 = 1. 40 3.5. 2 Example 1 of Section 2. 5, v = 1. 41 3.7. 1 A Kc EZ which does not satisfy assumption i) of 54 Comment 3.7.2. 3.8. 1 Zk I-Iz':I for n = 2, p = 2, z0 = (6, 2): A) X2 = 100, 59 B) X2 = 1000; IIP Selection Rule A. 3.8.2 Results for n = 2, p = 2, zo = (6, 2), X2 = 100: 60 A)ik zi(zi), B) IZk - Z*; IIP Selection Rule A. 3.8.3 Results for n = 3, p = 3, z0 = (6, 2, 2), X2 = 100, X3 = 10: 64 A) z- i<k Izil(zi), B) Izk l-lz'I, C) IZk - Zl:; IIP Selection Rule A. For k - 9, zk- z |l-zkl-lz'l.

Figure Page 3.8.4 Results for n = 3, p = 3, z0 = (6, 2, 2) X2 = 100, 65 X3 = 10: A) zI'I-IzklY(zk), B) Is(-zk) - z*l; liP Selection Rule A. 4. 2. 1 Notation for the case p(w)' E in the proof of Theorem 84 4. 2. 1. 4. 3. 1 A configuration for which Jzi+l |> Izil in GIP. 88 5. 7. 1 Initial conditions and isofuel contours for the minimum- 109 fuel example. 5. 7. 2 Geometry for minimum-fuel example. 112 5. 7. 3 6. and p for x(0) = (2, -1), 0 =. 4; BIP used on each 122 iteration of GIP. 5. 7. 4 6j and Pj for x(0) = (2, -1), 0 -. 4; IIP Selection Rule 123 J3 J A with p = 3 used on each iteration of GIP. 5. 7. 5 6j and pj for x(0) = (2, 0), 0 =. 4; BIP used on each 124 iteration of GIP. 5. 7. 6 6. and.j for x(0) = (2, 0), 0 =. 4; IIP Selection Rule A 125 with p = 3 used on each iteration of GIP. 5. 7. 7 4. for x(0) = (2, -1), 0 =.4: A) IIP Selection Rule A 126 with p = 3 used on each iteration of GIP; B) BIP used on each iteration of GIP. 5.7. 8 for x(0) = (2, 0), 0 =.4: A) IIP Selection Rule A with 127 p = 3 used on each iteration of GIP; B) BIP used on each iteration of GIP.

CHAPTER 1 INTRODUCTION Recent advances in engineering and science, especially in space technology, have stimulated much interest in optimal controls. The availability of modern computers makes it feasible to consider the actual computation of these controls. The central theme of this dissertation is the development of new and improved computational methods for a broad class of optimal control problems. Some of the early work in this area concerned rather specialized low-order problems in which an explicit solution could be obtained, For example, Bushaw [B7] solved the second-order minimum time regulator by geometrical constructions in the phase plane. More recently, because of the need for solving complex problems, a great variety of general computational methods have appeared in the literature. The method of dynamic programming introduced by Bellman [B8] has application to many optimization problems. However, the large storage capacity and long computational times that are usually required limit the practical usefulness of the technique. Gradient methods in function space also have been proposed [e. g., B 1, K1] and are quite commonly used, but are difficult to apply in problems with control constraints and usually exhibit slow convergence. The two-point boundary value problems which arise from the necessary conditions of Pontryagin [P2] or variational calculus provide the starting point for many other computational methods. For example, there are techniques [B9, B10, K2, P1] which consider certain mappings of the boundary values at the initial or final time. There are also methods which are based on the convexity of sets related to the reachable set of system states. Such convexity methods are the principal concern in this thesis.

Most of the convexity methods involve gradient-type minimization of scalar functions of n variables [e.g., B1, El, Fl, F2, N1, N3, N4]. These techniques possess inherent difficulties of step-size determination and exhibit slow convergence in the presence of poor problem conditioning. Furthermore, they require an assumption that the reachable sets be strictly convex. This assumption is not satisfied in certain important problems; even if it is satisfied, it may be difficult to verify. Another convexity method, which provided the motivation for this dissertation, is that due to Gilbert [Gl]. Gilbert's method and the procedures developed here are not gradient-type techniques and do not require strict convexity. Other desirable features are: rapid monotone convergence and the availability of computable error bounds. It should be noted that Gilbert's investigations in this area began by proving convergence for an extension by Fancher [F3] of a procedure due to Ho [H3]. The outline of the dissertation is as follows. In Chapter 2 a basic quadratic programming problem BP is stated and Gilbert's algorithm for solving it is presented in detail, The results of many numerical computations with this algorithm are also given. Gilbert's method requires a quadratic minimization problem on a line segment to be solved on each iteration. Chapter 3 contains an extension of Gilbert's method in which quadratic minimization problems are solved on convex polyhedra. The theory is carefully developed and includes a convergence proof and a sufficient condition for convergence in a finite number of iterations, Extensive computational results for a rather general example are presented. These results indicate that the extended procedure for solving BP is a significant improvement over Gilbert's method. In Chapter 4 a general problem GP which has application to a wide variety of optimal control problems is formulated. The basic idea of the formulation is due to Fadden [F1], but the

iterative procedure for solving GP and the proof of convergence are new, Finally, in Chapter 5 the iterative procedures are applied to optimal control problems of a broad class. Numerical results for a minimum fuel example are also given.

CHAPTER 2 THE BASIC ITERATIVE PROCEDURE BIP Gilbert [G1] has presented an iterative procedure for computing the minimum of a quadratic form on a compact convex set in finite dimensional space. This algorithm will be called the basic iterative procedure since it is fundamental to all the computing procedures developed in this thesis. The first five sections of this chapter are devoted to the theoretical details of Gilbert's method and except for minor changes are taken directly from his paper [Gl]. In Section 2. 6 the results of the author's extensive computer experimentation with Gilbert's procedure are presented. 2.1 Preliminary Comments on Notation Let x = (xl, xZ,..., x), y = (y,y,...,y ), z = (z, z,..., zn) be vectors in n-dimensional Euclidean space E and let wc be a scalar. n i i The following notation is employed: x y =.Z x y, the scalar product i 1-1 of x and y; Ix| = (x, x)2, the Euclidean norm of x, and Ix - y, the distance between points x and y; L(x;y) = {z: z = x + w(y - x), -oo < w < oo}, x ~ y, the line passing through x and y; N(x;w) = {z: z - x < W}, w > 0, the open sphere with center at x and radius w; N(x;w) = {z o Iz - x I W}, the corresponding closed sphere; Q(x;y) = {z: z. y = x. y}, y 1 0, the hyperplane (dimension n - 1) through x with normal y; Q (x;y) = {z z -y < x- y}, y 1 0, the open half-space bounded by Q(x;y) with outward normal y; Q (x;y) {z: z.y > x. y}, y i 0, the closed half-space bounded by Q(x;y) with inward normal y. In addition, if X and Y are arbitrary sets in E: aX denotes the boundary of X; X +~Y is the set z: z = x + y, xEX, y E Y}; X - Y is the set {z z= x -y, xE X, yY}.

2. 2 Contact Function; the Basic Problem BP Consider a set Kc E which is compact and convex. The realvalued function r(y) = max z. y is called the support function of K. zEK Since K is compact, (b(y) is defined for all y. Furthermore, it can be shown that rI(y) is a convex function on En, a result which implies that r{(y) is continuous on En [E2]. Let P(y), y i 0, be the hyperplane {x: x- y = r(y)}. Since z. y < j(y) for all z E K and P(y) nK is not empty, P(y) is the (unique) support hyperplane of K with outward normal y. For each y i 0 the set S(y) = P(y) nK is called the contact set of K and its elements are called contact points of K. It follows that S(y) is not empty, S(y) c aK, S(wy) = S(y) for o > 0. If for every y i 0, S(y) contains only a single point, K is strictly convex. A function s(y) mapping E into E is defined to be a contact function of K if s(y) E S(y), y / 0, and s(O) E K. From the preceding it may be concluded that: s(.) is bounded; s(y) = s(wy), w > 0; vb(Y) = s(y). y. Furthermore, on the set {y: y I > 0} each of the following is true if and only if K is strictly convex: s(.) is uniquely determined, s(.) is continuous. The continuity result is proved in [N1]. If for every y there is a method for determining a point x(y) E K such that x(y) * y = zEK z. y = flry), then it is said that a contact function of K is available. Such a determination, which corresponds to the solution of a linear programming problem on K, is possible for the convex sets which arise in a variety of optimal control problems (see Chapter 5). This availability of a contact function is essential to all the computing procedures which follow. Consider now the basic problem: BP Given: K a compact, convex set in En Find: a point z: E K such that I zI = min I zi zEK Since K is compact and I z is a continuous function of z, a solution z* exists. The following additional results hold:

THEOREM 2. 2.1 (Solution Properties for BP). i) z*' is unique; ii) I z' I- 0 if and only if 0 E K; iii) for z* j > 0, z: E aK; iv) for I z > 0 z = z*' if and only if z E P(-z) n K = S(-z). Proof: Properties ii) and iii) are obvious. Property i) is proved by contradiction. Suppose z"' and z" are distinct solutions. Then by convexity -z + 1 2z' K, which means lz - I| zII I =Z' I. But this implies I z' + 12 Z[2 2 I - Z +Z I Z2 which can be written z - Z 0 O an inequality which is only true for z'" = z'. Consider iv). The condition z E P(- az) nOK implies z C P(-z) - Q(z;z). But Q(z;z) is the support hyperplane for the closed sphere N(O; I z) whose outward normal is z and whose contact point is z. Therefore, Q(z;z) is a (separating) support hyperplane for K and N(0; Iz ). Thus K nN(0; z I) is empty. Since zE KnN(0; Izl), this implies z = z'o The steps of this argument may be reversed to obtain the converse result. Since I z is the Euclidean norm of z, the above minimization problem is equivalent to finding z* e K such that z;'. | = mIin ZZ zEK Thus BP is a quadratic programming problem. It differs from the quadratic programming problems which are frequently treated in the literature [e. g., Al, B2, H1, V1] in that the constraint set is not described by some set of known algebraic equations or inequalities. Information about the set K is obtainable only through a contact function s(-)o In actual computations it is generally not possible to obtain z*-. However, it is important to know if z E K is such that I z - zt_ I 5 E1 and z - z;:' < E2 where cl and Ez are specified positive constants. Expressions for determining if these inequalities are satisfied are given as part of the iterative procedures which follow. 2. 3 The Basic Iterative Procedure EBIP In this section Gilbert's algorithm for computing the solution to

BP is described, As a first step, let s( o) be a specific contact function of K and consider: P(z) = lz-s(-Z) z - (z-s(-z)), z-s(-z) 0 (2 3.1) =0 z - s(-z) =0 and y(z) = zi z. S(-Z), IzI > 0, z s(-z) > o (2. 3. 2) =0, -0 or Izi > 0, z s(-z)' 0. Thus (. ) and y( ) are functions which are defined on K. Figure 2.3.1 indicates their geometric significance: x = z + f(z)(s(-z) - z) is the point on the line L(z;s(-z)) with minimum Euclidean length; y( z)z is either the point L(O; z) n P(-z) or the origin, depending on whether or not L(O;z) n P(-z) is on the line segment connecting z and the origin. L(O; z);~igin) i L(z; s( - z)) P(-z) x = z +p(z)(s(-z)-z) Figure 2.3.1 Geometric significance of 3(.) and ~y( ). The functions P( ) and y( ) have the following properties. THEOREM 2. 3.1 Let K be the set described in BP and restrict z to K. Then: i) ((z) > 0; ii) P(z) = 0 if and only if z = z-; iii) O < y(z) < 1; iv) if O E K, y(z) - O; v) if O d K, y(z) = 1 if and only if

z = z'.; vi) y(z) is continuous. Proof: In this paragraph z always denotes a point in K. Later in this chapter (inequality (2. 4. 5)) it is shown that O < z - (z - s(-z)). Hence, i) and iii) follow from (2. 3. 1) and (2. 3. 2). For the time being assume Iz* I| > 0. The conditions P(z) = 0 and y(z) = 1 both imply z. (-z) = s(-z). (-z) = ri(-z) which requires z E P(-z). Since z E K, part iv) of Theorem 2. 2.1 yields z = z-o*. Reversing these arguments completes the proof of ii) for I z:' I| > 0 and of v). Now take I z*: I = 0. Inequality (2. 4. 4) then implies s(-z) * z s O which by (2.3. 2) yields iv). If P(z) = Othen it must follow from (2.3.1) that s(-z) ~ z = I zl Because of s(-z). z - 0 this implies z = 0 = z*. Since z = z* - 0 also yields P(z) = 0, the proof of ii) is complete. For jz |I' |z: I| > 0, the continuity of y(z) follows from (2. 3. 2) and the continuity of the support function rl(y) = s(y) - y. For I zt': = 0, it is trivially true from iv). It is of interest to note that p( ) may be discontinuous on K even though s(-) is continuous on K. See Example 3, Section 2. 5. For 6 a fixed number, 0 < 6 < 1, and z E K, define the closed interval I(z) = [min{6p(z), l}, min{(2 - 6)(z), 1}]. (2.3.3) Then Gilbert's algorithm may be stated as follows: The Basic Iterative Procedure BIP Let s( ) be an arbitrary contact function of the set K specified in BP. Take z0 E K and choose 6 in the interval 0O< 6 < 1. Then a sequence of vectors {Zk} in En is generated according to the rule k+ = k + ak((-zk) k) Z) k = 0,., Z (2. 3.4) where scalars a k are selected arbitrarily from I( zk). For the case 6 = 1 the selection of ak E I(zk) reduces to ak = sat (Zk)v where sat to is o for 0 _ X o 1 and is 1 for X > 1. Figure

2. 3. 2 gives the geometric interpretation of BIP for this case. N(_; I Zk I Y(Zk)) Origin. 2 Geometry of BIP for = s(-zk+ < z; if s(-zk) z Zk+1procedure is finite _ ~(zk)k1 zk is also clear that yzk (zk) >, k on each step upper Figure 2.3.2 Geometry of BIP for 6 - 1o if P(zk) > 0 an improvement is obtained on the kth step, i.e., Izk+ I < Iz k' if P(zk) = O0 zk = z*: and the iterative procedure is finite, i.e., the solution has been obtained in k steps. From Figure 2. 3. 2 it is also clear that Iz*[ < I Zk[. Thus on each step upper and lower bounds on zj' I may be computed. Notice that in applying the iterative procedure it is not necessary to know beforehand whether or not 0 e K. A more precise and complete statement of results is contained in the theorem of the next section. 2.4 Convergence Theorem for BIP THEOREM 2.4.1 Consider the sequence {Zk} generated by BIP. For k > O and k -- oo: i) zk E K; ii) the sequence {|zk |}is decreasing (IZk -> [zk+ l), I Zk -Z 1IzI, and IZkl = Zk+ I implies zk = z* iii) zk -- z; iv) IzkIY(zk) < Iz [ and IzklY(zk) -> Iz*I; v) IZk - zk I < /1 - Y(Zk) IZk and /1 - Y(Zk) IZkI 0~; vi) Is(-zk) - z':- - Is(-Zk) - (zk)zkl.

10 Since the bounds given in parts iv), v), and vi) are computable as the iterative process proceeds, they may be used to generate stopping criteria for the termination of the iterative process. Example problems show { IZk Y(Zk)} is not necessarily increasing. Thus Izkl - max I zi Y(Zi) is more satisfactory as an upper bound for I Zk Iz - z than lzkl - IZklI(zk). Since examples also show that { Izk z- *} and { I (-zk) - Z: |} are not necessarily decreasing, it is not possible to improve similarly the bounds given in v) and vi). Suppose IZ''l > 0 and s(.) is continuous in a neighborhood of -z:(the latter is certainly implied if K is strictly convex). Then it follows from the continuity of y(.) and part iii) that the upper bound in part vi) converges to zero. Thus {s(-zk)} may be used as an approximating sequence, an approach which may be advantageous in some situations. In addition it is clear from part iv) that Is(-zk) | - II Z': I s(-zk) I - max I Zi Y(zi), where the right side converges to zero. Therefore meaningful stopping criteria are available. Proof of Theorem 2.4.1 (due to Gilbert [G1]): First some basic inequalities are stated. From z*" E P(-z*), 0 # K, and s(-y) E P(-y), y # 0, it follows by the definition of P(.) that Z*. e Z* - Z. <z*, zE K, 0 K; (2.4.1) s(-y) ~ y <z zK y z, y EE (2. 4. 2) These inequalities lead to: a n I z' 12 < s(-y) ~ z*, 0 4 K, y E; (2.4. 3) n s(-y)- y s z.y, y E E; (2. 4. 4) s(-z) - - zE K; (2. 4. 5) ly- z:-' z+z':-. ~(y-z)y <y- (y-s(-y)), y-E; (2.4.6) Iz -z' < z- (z - s(-z)), z E K, 0 K. (2.4.7)

Inequalities (2.4. 3), (2. 4. 4), and (2.4. 5) are deduced from (2.4.1) and (2.4. 2 by obvious substitutions. From the identity ly - Z I + Z''. (y - z'') + y'(z*- - s(-y)) = y (y - s(-y)) inequality (2.4. 6) follows by (2.4. 4). Inequality (2. 4.7) follows from (2.4. 6) by use of (2. 4.1). Part i) of the theorem depends on ak E I( zk) which insures o sak < 1. Thus from (2.3.4), s(-zk) E K, and the convexity of K, Zk E K implies Zk+1 E K. Consider now the inequalities in parts iv), v), and vi). From (2.4.4) and the Schwartz inequality s(-y) - y l Y Iz I. Thus iv) follows from (2. 3. 2). The proof of the inequalities in parts v) and vi) makes use of z = zk E K. For s(-z). z > 0, z. (z - s(-z))= Iz z(1- (z)) and from (2.4.7) the inequality in v) is true. Now consider s(-z). z < 0, which corresponds toy(z) = 0. For z*: = 0, y(z) = 0 (Theorem 2.3.1) and v) holds as an equality; for z* / 0, the inequality in v) follows from (2. 4.1) which insures l - z-l =12 zl - 2z - z.: + 21z2>' Iz - I z <I l <' lz. If z' = 0 the inequality in vi) is trivially true. Consider now z* $ 0. If s(-z) ~ z < 0, part vi) reduces to -2s(-z) ~ z" + I z* 12 0 which is true by (2.4. 3). The following identity is easily verified I s(-z) - z" 1 = s(-z) - Y(z)z2 + z -2(s(-z). )2 + IZ': 2 _ 2s(-z). z*. Assuming s(-z)' z > 0 and using s(-z). z < Izl Iz*:- yields Iz-Z2(s(-z) 2 s z2 T s( z*_2 < Is 2 I z I. Thus Js(-z) - Z*J < Is(-Z) - y(z)ZI + 2( z2 - s(-z). zi) and by (2. 4. 3) the inequality in vi) follows. In order to complete the proof of the theorem, the function 1r() = z2 - Iz':Z = Iz - z':2 + 2z'. (z - z*:) (2.4.8) is introduced. For 0 B K inequality (2.4.1) gives 0 < Iz - z:12 r(z), z E K, (2.4.9) a result which is obviously true for 0 E K. In the following paragraphs it will be shown that {F(zk)} is decreasing and F(zk) -* 0. By (2.4.8) and (2.4.9) this proves the first two results in ii) and iii). From

12 (2.3.4) Zk+ I < IZk + ak s(-Zk) - zki so IZkI = Zk+ limplies ak s(-zk) - Zk O0. Thus either ak = 0 or zk = s(-zk), both of which yield P(zk) = 0 by (2.3.3) and (2. 3.1). Then part ii) of Theorem 2.3.1 implies zk = z4, which completes the proof of ii). The remaining results in parts iv) and v) follow from the known value of y(z'), the continuity of y(-), and part iii). For simplicity let A(z;C) = r(z) - r(z + a(s(-z) - z)) (2.4.10) and assume tacitly in what follows that z E K. Then from (2.4.8) z( z; ) = 2a( I z s( -z) ) _ a I z - s( z) Z (2. 411) Because the coefficient of a 2 is not positive, min) (z;a) is attained at the end points of I(z). It is readily shown that A(z; 6 P(z)) = ( z;( 2 - 6) P(z)). Thus from the definition of I(z), min A(z;a) = z(z;6p(z)), P(z)' 6 ~~~~~~~aEI(z) -1 (2.4.12) = (z; 1), 3(z) - 6 Equation (2.4.12) is now used to obtain a lower bound on A(z;), a E I(z). From (2.4.11) and (2. 3.1) it follows that t(z;6p(z)) - z - s(-z)l [z.(z- s(-Z))](26 - 6z). (2. 4.13) Let max KIZI - z21 (2. 4. 14) - z, zzEK2 denote the diameter of K and recall that 0 < 6 < 1. Then A(z;63(z)) > -2 6[z. (z - s(-z))]z (2.4.15) From (2.4.8) and (2.4.6) r( z) 21z - +2z (z-z) 2z (z - s(-z)) (2.4.16) (for z': = O this may be sharpened to 1(z) < z' (z - s(-z))). Thus

13 A(z;6 (z)) > 4[ 26r (Z). (2.4.17) For P(z) > 1, z (z - s(-z)) z - s(-z)I and consequently A(z;l) = 2zz (z - s(-z))- Iz-s(-z)2 > z. (z- s(-z)). Therefore (2. 4.16) yields (z;l) - i2r(z), P(z)' 1. (2. 4.18) Finally, utilizing (2.4.17) and (2. 4.18) in (2.4.12) yields -'1 -2 2 1 A(z; c) > min{ 4 26r z), r(z) (2. 4.19) tI( z) Letting z = zk in (2. 4.19), using (2. 3.4), and returning to (2.4.10), it is seen that Fk k+ min{ 26F(z k) 2 (Zk) (2.4.20) Therefore the sequence {f(zk} is decreasing and, since it is bounded from below by zero, has a limit point. Thus passing to the limit on the left side of (2. 4. 20) gives zero and therefore from the right side r(zk) - 0. 2. 5 Nature of Convergence; Examples This section gives some further results of Gilbert on the convergence of BIP. Theorem 2. 5.1 establishes upper bounds on the elements of the sequences { Izkl} and { Izk - z*I} Three example problems are analyzed to demonstrate still more fully the nature of convergence. Emphasis is on the case 0 B K, since it is most important in applications. THEOREM 2.5.1 Let ok = o0(1 + 4~ 60ok) 1, Oo = z02 _- Iz 12 (2. 5.1) and assume that Izo02 < zs:'2 + 26L26l, Then if {zk} is generated by

14 BIP, the following inequalities hold for k > 0: | Z k or B~lk | | (2. 5. 2) Izk- z, |< oNF. (2. 5. 3) The assumption on Izo0 is often met in practice. For example, it is easily shown that it must be satisfied if iz*j l -<(26 - 1)1. In any case, z0 may be interpreted as a suitable intermediate point in the iterative process, and inequalities (2. 5. 2) and (2. 5.3) may be used to estimate the subsequent rate of convergence. For I z*' > 0 and k > 1 inequalities (2. 5. 2) and (2.5. 3) imply zk I - Iz*I < Z.Iz"':l- l' k1, (2. 5.4) 1 1 lZk - Zj < 26- 2k- 2 (2. 5. 5) results which conform closely to (2. 5. 2) and (2. 5. 3) for k sufficiently large. In Examples 1 and 2, which appear later in this section, it is demonstrated that within a constant multiplicative factor it is impos - sible to obtain bounds on I Zk - IZ' and I Zk - z'* which approach zero more rapidly than those given in (2. 5. 4) and (2. 5. 5) Proof of Theorem 2. 5.1: Since Izol 0. From (2.4. 20) this implies F(zk+ ) (zk) - 6 (Zk) k > 0. Since 1 - 1 - 2 1 - ~4 b27r (1 + 42 6Yr)- for all F > 0, it is possible to write 7(Zk+l) $ r(zk)(l + 4 26r(zk)) k > 0. (2.5.6) But substitution shows that 0k is the solution of = - (1+: -26 )-1 (2.5.7)

15 with O0 = Izso - Iz 12 = r(z0). Thus comparison of (2. 5. 6) and (2. 5. 7) yields r(zk) < Ok, k > 0. Finally, (2. 5. 2) and (2. 5. 3) follow from (2.4. 8) and (2.4. 9). The complexity of the difference equation (2.3.4) makes it difficult to obtain more specific analytic results than those obtained in Theorem 2.5.1. Thus the remainder of this section is limited to the presentation and discussion of three example problems. Example 1. Take 6 = 1 and let K be the convex hull of three points in 2-space, (1, v), (-1, v), (0,1 + v), where v > 0. Clearly z- =(O,v) and I z;' = v. Simple inspection shows that the iterative process is finite (zl = z*) if and only if z0 is on the line segment connecting (1, v) and (-1, v). Moreover when the process is not finite Zk, k > 1, is determined by the scalar qk = IZk jl(zQ. Thus the second order nonlinear difference equation (2. 3. 4) may be replaced by a first order difference equation in Wk. It is not difficult to show that Tlk+ = <k( 1- V k)(1 + Vk + 2Zk) k > 1 (2. 5. 8) For qJk << 1 this equation is approximated by =qk+l = k(l + 2vk)-, an equation of the same form as (2. 5.7). These observations and some tedious, but straightforward, computations lead to (the notation o(o) means lime0 o(W) = 0) Zk -|z'I (Zvk) ~- + o(k ), (2.5.9) Iz -zk I -(2vk)-i + V2 + o(k1) (2.5.10) Equation (2. 5.9) demonstrates that is is impossible to obtain an upper bound on I Zk I- I zI which approaches zero more rapidly than (const.) k. For large k the upper bound in inequality (2. 5. 4) is conservative by a factor of sixteen. This factor can be traced to two sources each of which contributes a factor of four: in equation (2.4.15)

16. is an unsatisfactory estimate of Izk - s(-zk)I, in the derivation of inequality (2. 4. 6) the term y - (z*- - s(y)) has been omitted. For this example the upper bound in inequality (2. 5. 5) is a poor estimate be1 -1 cause it is order k 2 rather than order k It is also possible to show that z, I (zk)/Izkl (Zvk) + o(k) (2. 5.11) 1 1 ly(zk) Ikl k 2 +o(k 2) (2.5.12) By comparing (2. 5. 11) with (2. 5.9) and (2.5.12) with (2. 5.10) it is seen that in Theorem 2. 4.1, part iv) provides a reasonably good stopping criterion while part v) does not. Example 2. Take 6 = 1 and let K be the convex hull of three points in 3-space, (1, 0, v), (-1, 0, v), (0, 1, v), where v > 0. Thus z* = (0, 0, v) and Iz l = v. The iterative process is much the same as in Example 1, the points zk E K, k - 1 being determined by the scalar =k = |kZ(zk. The first order difference equation for Lk is (2. 5.8) with v = 0. By using the fact that k % ~ (1 + 4q0ok) 2 is the solution of _,. 1 _ 1 k+l = pk(l + 4pk4) 2 and that (1 + 4 ) 2 _1+ 2 for << 1, the following results can be derived: Izkl I- IZ1 = (8vk) + o(k ), (2. 5.13) k = (2k2) + o(k 2), (2. 5.14) >z- - Y((zk)lzk = 3(8vk) + o(k ), (2.5.15) 1 1 N/1 - (zk) zk (2k) 2+ o(k 2). (2. 5.16) Equation (2. 5.14) shows that the asymptotic behavior of I zk - z matches the bound given in inequality (2. 5. 5), except for a multiplicative factor of eight. The bound given in inequality (2. 5.4) is conservative by a multiplicative factor of 64. Comparison of (2. 5.15) with

17 (20 5.13) and (2. 5.16) with (2. 5.14) shows that parts iv) and v) of Theorem 2.4.1 both provide reasonable stopping criteria. Example 3 (The Hyperparaboloid Problem) Take 6 = 1 and in nspace let 2 n K =zlz' v + 2 F, (ZYiX2 Z' Z 2v}, VX2,XI, X3 n> ~i=2 (2. 5.17) In the neighborhood of the optimum z,* = (v, 0,..., 0), K is the elliptic n - 1 hyperparaboloid z1 = v + a F (z )Zi_ where X k are the principal 2 1 i- 2 n2 the principal radii of curvature at the vertex z*. Many other convex sets K have a boundary surface which in the neighborhood of z^' can be closely approximated by such an elliptic hyperparaboloid (see Chapter I of Busemann [B6]). Thus this example is of general interest. For y1 < 0 and - (y) (y)2 <K v it is easy to show that 1=2 n s'(y) = v+ (yl)i(yl)z s (y) = (yl)l, i=...n. (2. 5.18) i=2 Let X = imax {X.} and assume the conditions = v Ni +2vX > IyI - >y v y' (2. 5.19) are satisfied, which in turn imply z i - Y i(Y )2 < v. Thus (2. 5.19) defines a set on which (2. 5.18) is valid. Using this fact, z1 > v for z E K, and (Z.3.1) gives n (z' - v)2 + V(zW- V)+ Fa (I+ I 1Z)- X )(z i) )i-Z( 2 ). )() n 2 -v)2 (z' - v)2 Z C (Z1+(z X. +(z +(z ) i)( zl) +)22zC (Z)_-2 X.zz)2)2 1i 1 i=2 1 (2. 5.20) 7,; 7,v'v 7. K, IZI < r

18 n Because z E K, |zj <, imply z1 > v and - (z)2 (zl)2 <v it follows i12 i that n (z1 -v)2 +; (z1)2 ((z) - i= 2(5)'_-Z -2 n 1 - (z V)2 +( + 2Xv +Xv )(z) 1+ 5 -+X v 12 (2. 5. 21) z / z, E K, Izl < f. Because P(z':) = 0 this inequality implies that P(z) is discontinuous on K at z*. By starting with (2. 4.13) and repeating the derivation of Section 2.4 with [z. (z - s(-z))] z - s(-z)j2 = P(z) > ( it can be shown that (Zk+z ) F(z)(l - - 6), k # zK (2. 5. 22) For zk =, F(Zk+ ) = 0 and (2. 5,g22) is trivially true. Thus r(zk) < (ZO)(l -j 2) I k' O, zo E K, jzo| <. Using (2.4. 8) and (2.4. 9) this leads to IZkI -| |zsvlo- 0 (I - I6)k, (2.5. 23) IZk - Z,-I< %0 (1 6) where 00 is given as before in (2. 5.1). Since ( > 0 inequalities (2. 5. 23) and (2. 5. 24) guarantee that the convergence of { Zk } and { Izk - z' } is geometric. However, the guaranteed rate of convergence is not rapid if P << 1, i.e., vZ << X2 2. 6 Numerical Results for BIP Actual computations with BIP were carried out for Example 3 of Section 2. 5 at the University of Michigan Computing Center. The programs were written in the MAD language and an IBM 7090 digital computer was used. In all the experiments v = 1, i. e., z' =(1, 0,...,0)

Also, the extent of K was increased beyond z1 = 2v so that (2. 5.18) is valid even though (2. 5.19) is violated. Data was obtained for n = 2 and 1 3figures and 3 and various combinations of X2, 3, Z z0, zand z0. The figures and tables included in this section contain an important part of these results. In optimal control applications (see Chapter 5) the evaluation of a contact function is the most time-consuming part of the iterative procedure. Since such an evaluation is required on each iteration of BIP, the number of iterations to satisfy certain error criteria is used as a measure of the speed of convergence. The results for n = 2 are presented first. Figure 2. 6. 1 shows the strong dependence of BIP on the parameter X2v, convergence becoming very slow as X2 v -- oo. Additional data for various performance measures and X2 values, given in Tables 2. 6. 1 through 2. 6. 4, illustrate further this dependence. It is interesting to note that 3(zk) shown in Figure 2. 6. 2 is large immediately preceding major improvements in I ZkI |- z I. Table 2. 6. 5, which displays the case X2 = 100 for different points z0, shows a somewhat random dependence of BIP on the initial point z0. As a rough average, about 16 iterations are required for a decade of improvement in zk, -Iz:j (n = 2, X2 = 100). The results for n = 3 are presented in Figures 2. 6. 3, 2. 6. 4 and Tables 2. 6. 6 through 2. 6. 9. The behavior of I Zkl- z, I zk - z |, Iz-* - ZkJl(Zk), and ls(-zk) - Z" shown in Figures 2. 6. 3 and 2. 6. 4 is typical of that observed in all computations with BIP. It can be stated that I z' zk Y(zk) decreases most rapidly, followed in order by IZkl-IIzI, IZk - z*, and Is(-zk) - z*j. The results indicate that convergence is quite slow when v << X. Furthermore, compared with the parameter Xv, the ratio X2 /13 has little effect. By way of contrast, for iterative procedures of the gradient type [e. g., B1, El, F1, F2, N1, N3, N4], Xz /X3 has an important influence.

20 Attempts to improve convergence for this example by taking 6 i 1 and choosing ak E I(Zk) in a variety of ways were unsuccessful.

. 1 1-4 10 " — 10-5 lo-S I -L —Z O 40 50 10 l- 20 30_~= 10, B) 3z= 100 Figure Z.6.1 zkl- zAl forn =2, zo=(6, ): A) Xz, B) Z= C) Xz = 200, D) Xz = 500, E) Xz = 1000; BIP.

22 10 1-. 1.01 i i I / iiI I I I k 0 10 20 30 40 50 Figure 2.6.2 (zk) for n - 2, z0 = (6, 2), k2 = 100; BIP.

23 Table 2. 6. 1 Number of iterations to satisfy Zkl- lzI:n = 2, z0 = (6, 2), BIP. 0 1. 1.01 10-3 104 10-5 10 10 4 13 17 18 27 29 40 100 2 9 25 42 60 60 75 200 18 49 98 147 148 167 216 500 36 113 126 134 150 158 241 1000 69 96 200 215 1 281 376 443 Table 2. 6. 2 Number of iterations to satisfy I Zk - z I C n = 2, zo = (6, 2), BIP. E x 21 1.01 10-3 10-4 10-5 10-6 10 4 14 19 26 39 >40 100 4 11 29 46 64 83 >98 200 18 49 99 151 166 213 >269 500 38 115 129 139 157 165!>300 i 1000 71 98 203 224 301 443 >500

24 Table 2. 6, 3 Number of iterations to satisfy j z" | - k I t zk) (z n = 2, zo = (6, 2), BIP,,,~1 1 0 1 10 3 1004 -105 10 6 10 0 9 16 16 23.26 33 100 0 1 8 8 24 41 59 200 0 1 3 17 49 99 146 146 500 0 2 34 112 112 125 i41 1000 0 4 68 85 95 199 21I Table 2. 6. 4 Number of iterations to satisfy Is(-zk) - z' I n = 2, zo = (6, 2), BIP, IS;Tt 1 1.01 103 10 4105 10' 10 9 16 26 39 >40 100 8 24 59 84 >98 200 17 146 1 166 253 >269 500 112 125 157 299 >300 1000 85 199 392 487 >500

25 Table 2. 6. 5 Number of iterations to satisfy Izkl I I < n = 2, X2 = 100, BIP. 1 -3 -4 -5 0-6 z C 1 o 1 010 10 10 10 6 1 18 19 44 44 62 102 137 6 1. 5 5 36 57 79 90 100 129 6 2 2 9 25 42 60 60 75 6 2. 5 13 24 34 34 44 49 49 6 3 4 11 30 30 51 69 92 6 3. 5 11 34 34 34 34 49 49 6 4 18 22 34 34 59 81 81 6 4. 5 10 10 22 37 38 71 89 6 5 3 15 47 47 47 47 47 2 6 5 14 33 61 62 115 136 3 5. 57 25 67 72 89 122 160 160 4 4. 9 6 16 28 59 I 93 106 134 4.5 4. 44 10 20 62 73 88 113 137 5 3. 87 3 23 50 50 71 92 121 5. 5 3. 12 8 j40 68 102 120 136 175 6. 2 1.25 10 21 23 36 36 36 36 Note: For z0 = (6, 2) and the last seven cases in the table, IzoI2 40

26 1.01 10-3 L lX A -4 L-5 LO k 0 10 20 30 40 50 Figure 2.6.3 Results for n = 3, zO = (6, 2, 2), X2 = 100, X3 = 10: A) /i I" z i l (zi ), B) IZkl lZ' lC) Ilzk- z'l; BIP, For k' 14, 1zk - z:kI I-kl: |.

27 10 1.1I.01 10 -3 A 1o-4 L 0 10 20 30 40 50 k Figure 2.6.4 Results for n = 3, zo = (6, 2, 2), X2 = 100, X3 = 10: A) I z" IzkY(Zk), B) Is(-zk) - z': I; BIP.

28 Table 2. 6. 6 Number of iterations to satisfy I Zkl-z:- l -< E n = 3, z0 = (6, 2, 2), BIP. E X2 X 1. 1. 0 1 10 10 1 1 1 2 3 3 5 5 5 10 10 4 9 21 28 31 38 41 100 100 10 19 31 59 59 74 111 1000 1000 82 82 126 216 250 290 340 100 1 2 6 24 35 46 46 57 100 10 11 15 21 27 52 73 81 100 50 20 33 78 124 124 151 151 100 90 1 5 38 38 58 112 137 137 1000 1 14 198 218 218 274 321 i 321 1000 10 83 174 207 229 267 29 8 3 59 1000 100 81 88 108 197 221 273 358 Table 2. 6. 7 Number of iterations to satisfy I Zk - zI < E n = 3, zo = (6, 2, 2), BIP. -3 - -6 X2 XL 1.1.01 103 10-4 10 10 1 1 1 3 4 4 5 5 >5 10 10 4 11 25 30 40 > 41 100 100 10 20 35 63 73 >111 1000 1000 82 84 130 247 309 >340 100 1 2 15 36 45 > 57 100 10 11 17 26 47 64 92 >101 100 50 20 35 81 127 146 >151 100 90 17 38 42 67 132 >137 1000 1 16 199 232 238 3 18 >321 1000 10 84 177 216 249 291 >3 59 = 1000 100 81 90 112 203 259 337 >358

29 Table 2. 6. 8 Number of iterations to satisfy J z zk(zk) - E n = 3, zo = (6, 2, 2), BIP. k3 " " 1 1 1... Z..01......10 10 _ 1 1 0 2 2 2 4 4 4 10 10 0 5 8 20 27 30 37 100 100 0 9 9 18 35 58 58 1000 1000 0 26 47 81 81 125 165 100 1 0 1 19 3 1 38 45 45 100 10 0 1 1 14 14 26 51 51 100 50 0 13 1 9 67 77 123 123 100 90 0 14 37 37 37 57 127 1000 1 0 13 41 146 197 241 255 1000 10 0 34 82 88 173 223 264...:r "_ _w _> * * v -. _, r -._ 10000 100 00 0 9 80 80 87 107 196 Table 2. 6.9 Number of iterations to satisfy Is(-zk) - zf E; n = 3, zo = (6, 2, 2), BIP.' 1 i 3 1.01 3 10-4 10-5 10 i 1 1_ i _O i 2 2 4 4 4 4 > 5 10 10 l 5 [ 20 30 40 >41 _ 1 00 i_ 100 9 1305 73 > 1 1 1; 1000 1000 81 125 1 165 29 21>340 100 1 1_ 1 | 23 | 45 56 >57 - t____r -t* —- -- t - - -- 10 10 1 4 26 51 921>101t100 50 19 77 123 150 1>151 t__ 4 _ _ _- + i <. 100 90 37 I 37 1 27 i 136 >137 1000 1 96 i 243 1 286 i>321 1 1000 10 159 241 1295 >359 1000 100 80 107 1232 _

C HA PTER 3 THE IMPROVED ITERATIVE PROCEDURE II P The examples of Section 2. 5 and results of Section 2. 6 indicate that for many problems BP the convergence of BIP is not rapid. In Example 3, Section 2. 5 slow convergence is obtained with BIP for cases in which the surface aK at z* has at least one principal radius of curvature large compared with I z |. For problems in which S(-z'*) contains more than one point such as Examples 1 and 2, Section 2. 5, the convergence of BIP is especially poor. Since a large number of convex sets K which occur in practical problems either have a boundary surface near z* approximately like that of Example 3 or have contact sets containing more than one point (K not strictly convex), it is important to seek ways to improve BIP. By using more than one contact point at each iteration to gain information about AK, it is possible to develop an iterative procedure for solving BP which exhibits much more rapid convergence than BIP. Such an improved iterative procedure II P is discussed in this chapter. First some background material on convex polyhedra is given after which II P is stated and shown to converge. In the last two sections a sufficient condition for finite convergence of II P and the results of many numerical computations are presented. 3. 1 Convex Sets and Polyhedra It is convenient to introduce here some basic definitions and results for convex sets and polyhedra which will be needed in the presentation of II P. The proofs of all the results except Theorem 3. 1. 2 are omitted since they can be found in standard references such as Eggleston [E2] and Valentine [V2]. 30

31 If X is an arbitrary set of points in E, then the convex hull of X, written AX, is the set of points which is the intersection of all the convex sets that contain X. AX is convex, and a necessary and sufficient condition that X be convex is that X = AX. Furthermore, if X is compact, then AX is compact. The convex hull of a finite number of points from En is called a convex polyhedron. This set is compact and may be viewed as the intersection of a finite number of closed half -spaces. Thus it has a representation in terms of a finite set of linear equations and/or inequalities. Given yl, Y2z, y m E E n the convex polyhedron A{Y1, Y2,. *, Ym} is identical with the set of points of the form y = m m LCirYi, where. - 1 and cr-. O (i = 1, 2,...,m) Each such point y i=1 i 1=1 1 1 is said to be a convex combination of yl, Yz,, Ymo A point y E En is said to be an extreme point of a convex set X if y E X and there do not exist two distinct points yl, y2 E X such that Y e A{Yl, y}, y / Yl, Y / yz. The set of extreme points of a convex polyhedron is finite and is called the skeleton of the convex polyhedron. Furthermore, a convex polyhedron is the convex hull of its skeleton. n The dimension of a convex set Xc E, written dim X, is the largest integer m such that X contains m + 1 points y1 y, ~. e. Ym+l for which the collection of vectors {y1 - Ym+i ~Ym - Ym+i } is linearly independent. It follows that dim X c no A non-empty set in E is called a linear manifold if it consists of a single point or if for every Yl, yz in the set, yl i Y2, the line L(y-ly2) is in the set. If dim X = m, then X is a subset of an m-dimensional linear manifold (i. e., a point if m = 0, a line if m = 1, a plane if m = 2, a hyperplane if m = n - 1, the whole space En if m = n) The relative interior and relative boundary of a convex set X having dimension m are defined to be the interior and boundary of X relative to the m-dimensional linear manifold which contains X.

32 A convex polyhedron in En having dimension m is called an mpolyhedron. The relative boundary of an m-polyhedron is the union of a finite number of (m - l)-polyhedra. Furthermore, the skeleton of each (m - 1)-polyhedron is a subset of the skeleton of the m-polyhedron. A simplex is a special case of a convex polyhedron and is defined as follows: an m-dimensional simplex, or more briefly an msimplex, in En (m < n) is the convex hull of m + 1 points from En which do not lie on a linear manifold of dimension m - 1. The set of m + 1 points is the skeleton of the m-simplex. THEOREM 3.1.1 (Carath6codory [C1]) Let Xc En If y e AX, there is a set of m points yl, y2,...m all belonging to X with m s n + 1 such that y is contained in the (m - 1)-simplex A{yl, Yz,..o, Ym) The following theorem is a collection of some rather obvious results on convex polyhedra which will be useful in the sequel. THEOREM 3.1.2 Given m points yl, Y,...,ym e E, consider the convex polyhedron H = A{y1, Yz,. Ym) Let the dimension of H be q. Then: i) q < min{m - 1, n}) ii) the skeleton of H, which is a subset of {Y1, Y2..., Ym), contains at least q + 1 points; iii) if q = n, AH is the union of a finite number of (n - 1) -polyhedra; iv) if q < n, aH= H (a q-polyhedron); v) if m - n and Al a2,... _, m rrlon n( - n)!, denote the convex polyhedra formed b the convex hull of n points chosen from y,,y..., Ym}, then aHc U vi) if m n and m m - y e aH, then y has the representation y = ~ iYi- where ( = 1, ai - 0 1=1 1 1i =1 1 1 (i = l, 2..., m), and at least m - n values of -i, 1 i < m, are 0; vii) if m n +1 and 1, 2,. = (n l. n-l)' denote the convex polyhedra formed by the convex hull of n + 1 points chosen from {Y1Yz, ~ ~ Ym}, then H = Ao j=1 J Proof: Clearly q < n so i) is true for n < m - 1o Consider n > m - lo Ifq > m - 1, then there exist xx2,. oO., xm+ E H such that

33 the set of vectors {xl - xm+,l... x - Xm+l} is linearly independent. m m Every y E H has the representation y = o.y., cr-. 1. That is, m m m 1=1 1 1i= 1 Y = yi(l - -.i)+Z C.iy Yi +~ i(Yi Y1i) Thus each vector x, - i=Z2 i=z i2 1 3 x (j = 1, 2,..,m) can be written as a linear combination of yi - Yi (i = 2, 3,..., m), contradicting the fact that the x. - x+ are linearly independent. Hence, q < m - 1 which completes the proof of i). Consider ii). Suppose there is a point x # y., 1 $ i < m, which is mn 1 m an extreme point of H. Since x E H, x -= iy where C r. = 1, 1 i i=1 1 - > 0 (i = 12, 2..,m). But x f Yi3 1 0. For definiteness, suppose c. > 0. From -m 1 -1 1 (1 ~- c1) =2]i = 1 and (1 - >l). >O(i = 2, 3,..., m), it follows that i=2.lm 1 y = (1 - E1)-1?OiYi C H. Then x = -lyI + (1 - c1 )y where x t Yl and x I y, violating the assumption that x is an extreme point. Let p be the number of points in the skeleton of H. Suppose p < q. Since q < n, p < n. Then H is the convex hull of p points, and part i) implies dim H < p - 1. But p < q yields q = dim H < q - 1, a contradiction. Thus ii) is established. If q = n, aH equals the relative boundary of H and thus is the union of a finite number of (n - l)-polyhedra. This proves iii). In iv) the fact that H is a subset of a q-dimensional linear manifold, q < n, implies that if y E H and E > 0 then N(y;E) contains a point not in H. Thus aH = H. Consider v). By iii) and iv) aH is the union of a finite number of convex polyhedra each having dimension < n - 1. From the remarks following the definition of m-polyhedron, the set of extreme points of each of these convex polyhedra is a subset of {yl, yz,.., Ym Thus given y E aH there is a set X c {y1,yz,...,ym} such that y E AX and X is contained in a linear manifold of dimension n - 1. Since dim Xs n - 1 Theorem 3.1. 1 applied to X yields: there exist y, yz,....,yp e X with p < n such that y E a{y,z.y..., }. But for p < n and any set

34 {Yi-,-.-, Yp} of points from {y, y2,, ym} A{Yi, Yz,oo, p} -.j A and therefore v) is true. Part vi) follows at once from v) since y E aH implies y is contained in at least one A., 1' j < m. Thus y may be written as a convex comJ m bination of n points from {Yi, Y2,.. o Ym} and the remaining oi in -.iYi may be set equal to O. Consider vii). Let X = {Yl, Yz,.., Ym}. Then Theorem 3. 1.1 yields: given y E AX = H there exist y.Y-2,.. Ep ~ X with p < n + 1 such that y E A{ Y-1, Y2., }. If A,, 1 s j s -, are defined as in vii), - 3 then {Yl,yz,...,yp } c Aj for any set lY P n+l, of points from X. Hence Hcj U A.. Now suppose x E U A.. Then x E A j=1 j j-1 j j for some j, 1 < j < 2, which means x is a convex combination of n + 1 points from X. By giving the remaining m - n - 1 points in X a zero coefficient, x can be written as a convex combination of all the points in X. Hence x E H and A. c H. This completes the proof. j=1 J 3. 2 The Subproblem 1 On each iteration of II P the following minimization problem, which is a special case of BP, will occur. SUBPROBLEM 1 Given: H, the convex hull of m points Y Yz.,y frEnomE Find: a point y' Hsuchthat y': = yEHYJ Since H is compact and y Il is a continuous function of y, a solution y*: exists. Let PH(x), x / 0, be the support hyperplane of H with outward normal x. Then: THEOREM 3.2.1 (Solution Properties for Subproblem 1) i) y*4 is unique; ii) | Y':= 0 if and only if 0 E H; iii) for Iy':- j > 0, y:- E aH; iv) for Iy': - > 0, y = y if and only if y E P(-y) nH; v) if m n and m m y*E aH, y': has a representation y'' =iC.iYi where ~ c. = 1, -= ji=1 1

35.-*> 0 (i 1, 2,..., m), and at least m - n values of -. 1 i < m, are 0; 1n- n vi) if m > n, Jy' I > 0, and y* = Tlici where H - m, Z F. = 1, g. > 0 =1= i = 1 1 (i = 1, 2,...,n), and Yi (i = 1, 2,...,n are vectors from {Y1,yz,...ym}I then y. (i = 1, 2,..., n) are contained in PH(-y*). - 1 - H Proof: Parts i) through iv) follow by the same arguments used in proving Theorem 2. 2.1 and v) comes directly from part vi) of Theorem 3.1. 2. Consider vi). By iv) y' E PH(-y*') so vi) is clearly true for n = 1. Thus take 1 < in m. Suppose any Yi, say yl, PH(-y:) Then (yl - y*") Y* y i 0. This, Yj E H, and the fact that PH(-y*) is a support hyperplane of H with outward normal -y*' imply (Yl - y*)' Y*- > 0, n- n Since y'* = i-i where Z.i = 1, - 2, and >0 (i 1, 2,..,n-, i-1Yi i=1 1 1 - 0(i = 2,3, 3 E...,n), it follows that y E H which implies (y - y':) ~ y > 0. But y;:- = Fly l + (1 - F, )Y, 0 < rl < 1, only if there exists E > 0 such that y*' - y = E(yl - Y) Then (Yl - y*") * y:" > 0 yields 1(Y' - y)' y* > 0. This leads to (Y - y'). y' < 0, which contradicts the earlier result and completes the proof. It is important to note the distinction between BP and Subproblem 1, both of which are quadratic programming problems on a compact, convex constraint set. The set K in BP is described only by a contact function s( ) of K whereas the convex polyhedron H in Subproblem 1 is the convex hull of m known points. Thus Subproblem 1 is much simpler than BP. It is shown in Section 3. 6 that Subproblem 1 is amenable to solution by standard quadratic programming techniques. 3. 3 The Improved Iterative Procedure IIP In this section the improved iterative procedure for solving BP is described.

36 First let s( ) be a specific contact function of the set K in BP and consider: ( z) = Iz| z. s(-z), z 0 (3.3.1) =0, z =0 Thus the function ('), which is closely related to y( ) of (2. 3. 2), is defined on K. Geometrically, 1I4(z) |, z J 0, is the Euclidean distance from the origin to the support hyperplane P( -z) of K. The following properties also hold. THEOREM 3.3.1 Let K be the set described in BP and restrict z to K. Then: i),(z) z zY(z) -< z*' I; ii) if [(z) > 0, 1X(z) =- Iz (z); iii) I[(z;') = Iz;: J; iv) if 0 B K and Ep(z) = Iz;: J, S(-z) = S(-z*:); v) if 0 4 K, [Jl(Z) is continuous on K. Proof: In this paragraph z always denotes a point in K. The inequality Iz Iy(z)' Jz' I was shown in the proof of Theorem 2. 4.1. Part ii) and the remaining result in i) follow from equations (2. 3. 2) and (3.3.1). Consider iii). If 0 E K, z;: = 0 and clearly p(z;,) =- z* = 0. If 0 K, part v) of Theorem 2. 3.1 yields y(z;'') = 1. But y(z) > 0 implies k(z) = IZ Y(z) > 0 and thus 4L(z*) = Iz~ IY(Z) = Iz* I. In iv) 0 q K implies z z-* > 0 so P(-z), the support hyperplane of K with outward normal -z, is defined. But P(-z) = {x: x- Iz-l (-z) = s(-z) ~ Izl- (-z) = -[(z)} is also the support hyperplane of N(O; L(z)) with outward normal z and contact point p.(z) Izl-lz. Thus if p.(z) = z' i, P(-z) is a (separating) support hyperplane for K and N(O; I z* ). This implies P(-z) = P(-z*) and S(-z) = S(-z*'). Consider v). From (3.3.1) and the continuity of the support function i](y) = y. s(y), y E E, it follows that L(z) is continuous except possibly at z = 0. For 0 4 K, [z| > 0 for all z E K so the proof is complete. Now consider: The Improved Iterative Procedure II P Let s( ) be an arbitrary

37 contact function of the set K specified in BP. Take zo E K and choose a positive integer p. Then a sequence of vectors {Zk}, k = 0,1, 2,..., in En is generated as follows: Step 1 Select any p vectors Y1(k), yz(k),.. yp(k) in K and let Yk {yl(k),yz(k),..., yp(k)}, (3.3.2) Hk {yl(k),yz(k),...,yp(k), s(-zk), k} (3. 3.3) Step 2 Solve Subproblem 1 on the convex polyhedron Hk: find zk+ E Hk such that Izkl + m I zEHk Steps 1 and 2 constitute one iteration, called iteration k, of II P. Note that II P differs from BIP (with 6 = 1) only in the fact that Zk+ is obtained by minimizing over Hk instead of A{s(-zk), Zk}. There are a great variety of selection rules for choosing the elements of Yk in Step 1 and each choice of a selection rule gives a different version of IIP. As is discussed in Section 3.5 the function 1j(.) forms the basis of selection for several versions of II P which exhibit good convergence. 3. 4 Convergence Theorem for II P THEOREM 3.4.1 Consider the sequence {Zk} generated by IIP. For k > 0 and k - co: i) zk E K; ii) the sequence { zk } is decreasing (IZk - I Zk+lI)' IZkil l[ZI and Zk = Izk+l implies zk = iii) zk - z:'; iv) zzklY(zk) Iz'' I and ZkY(Zk) z'; v) zk - z < -y(zk) IZkl and 1 - Y(Zk)Izk > 0; vi) Is(-zk) - z':' Is(-z ) - Y(zk)Zkl. Furthermore: vii) 4L(zk) < IzkIY(zk); vii) if 0 f K, (zk) zIt should be observed that parts i) through vi) are identical to the corresponding parts of Theorem 2. 4.1 (convergence theorem for BIP).

38 The remarks in the two paragraphs following the statement of Theorem 2.4.1 also apply here. Proof of Theorem 3.4. 1: Since Yk cK and s(-zk) E K, it follows from the convexity of K and the definition of convex hull that zk E K implies Hkc K and Zk+1 E K. Thus by induction z0 E K proves part i). As in the proof of Theorem 2.4.1, the inequalities (2.4.1) through (2. 4. 7) and zk E K yield the inequalities in iv), v), and vi). Part i) of Theorem 3.3.1 and i) prove vii). Consider ii). Since II P differs from BIP (with 6 = 1) only in the fact that zk+l is obtained by minimizing over Hk instead of A{s(-zk), Zk} it follows from A{s(-zk), Zk} cHk that i Zk+i I using II P | Zk+ I using BIP (6 = 1). By comparison of the sequences { zkJ}, the first two results in ii) are a consequence of the corresponding results in part ii) of Theorem 2.4.1. Now suppose I kl = Zk+ I and zk f z*. Since izkl >- k z*> 0O Zk Z z*; implies zk I > 0 so that support hyperplanes P(-zk) and PHk(-z k) of K and Hk respectively are defined. Clearly s(-zk) E P(-z k) Q(s(-zk);zk) and part iv) of Theorem 3.2.1 yieldszk+lE PHk(-zk+ ) Q( k+;z. From Zk, Zk+1 E Hk I k I = I minizI, and part k+i k+i kk1 kIzHk k i) of Theorem 3. 2.1, it follows that zk = k+ and Q(Zk; Zk) =Q(z k+; Zk+l) Then Q(Zk;zk) is the support hyperplane of Hk with outward normal -zk and s(-zk) E Hkc Q-(zk;zk). Furthermore, Hkc K implies Q(Zk;zk) c Q (s(-zk); Zk). The last two statements can be true only if Q(Zk;zk) = Q(s(-zk);zk). Hence Zk E P(-Zk) which by i) and part iv) of Theorem 2. 2.1 implies zk = z*. This contradiction completes the proof of ii). Let the function F(z) satisfying (2. 4. 8) and (2. 4. 9) be introduced again. Part ii) and (2. 4. 8) imply that {r(zk)} is decreasing and r(Zk) - 0. By (2.4.9) this proves iii). The remaining results in parts iv) and v) follow from the known value of y(z*I), the continuity of

39 )y(), and iii) Parts iii) and v) of Theorem 3. 3.1 and iii) give viii). Thus the proof is complete. 3. 5 Selection Rules for IIP In this section several rules for selecting the vectors y (k), y (k),..,y p(k) E K in Step 1 of iteration k of II P are presented. Each selection rule yields a particular version of II P and it is desired that the iterative procedure converge much more rapidly than BIP. Consider the case z' E AK which is most important in applications. If after a few iterations of II P the surface aHk in the vicinity of Zk+l closely approximates AK in the vicinity of z-:', then it is likely that I I P will exhibit improved convergence. If aHk is to approximate WK, the dimension of Hk must be sufficiently large, namely n, and Yk must include boundary points of K. To illustrate these remarks, consider Figure 3. 5.1 which is Example 3 of Section 2. 5 with n = 2, v = 1, kz = Io In Figure 3.5.1(a) BIP is shown, where dimHk = 1 and convergence is slow. In Figure 3.5.1(b) IIP with Selection Rule A (to be described subsequently) is shown, where dim Hk s- 2 and convergence is notably improved. An even more startling improvement is exhibited in Figure 3. 5.2 which is Example 1 of Section 2. 5 with v = 1. Theorem 3.7.1 shows that when K is a convex polyhedron, IIP (with a suitable selection rule and contact function) converges in a finite number of iterations Furthermore, the extensive numerical results of Section 3. 8 provide strong evidence that II P is far superior to BIP. Let the p points in Yk be contact points of K. Observe that &Hk is a better local approximation to aK for larger values of p. However, the larger p is, the more difficult it is to solve the Subproblem 1 in Step 2 of II P. The computational results of Section 3.8 indicate that convergence is good for p = n and little improvement is obtained for p > n. The desirability of choosing p = n is also evident from the

40 0 5(-Zi) z3 s\s(- ZO) H1 Z3 S(-z2) (a) BIP. Hz z' O Z 3 (-bP c R A p (b) IIP Selection Rule A, p = 2. Figure 3. 5.1 Example 3 of Section 2 5, n = 2, v g - 1. finite convergence material in Section 3 7. In optimal control applications (see Chapter 5) it is advantageous to limit the number of times the contact function is evaluated. Thus in the selection rules which follow Yk+' k: 0, contains every point in Yk except perhaps one. There remains the question of how to reject one contact point in favor of another. The approach in the selection rules given here is to use i(z) as an indication of the quality of the contact point s(-z), z E K. Roughly speaking, contact points corresponding to larger values of (.) are preferred. Other quantities, e.g.u Is(-z), i z, Y(z), may be suggested for judging the merit of s( z). Howe vser,

41.0.0 s(-z1) Z' S(=- s( -ZZ) ZO Ho (a) BI P -0 s(-zl) z':-' = Zz s(-Zo) (b) II P Selection Rule A, p = 2. Figure 3.5.2 Example 1 of Section 2.5, v = 1.

42 careful examination of example problems such as Example 3, Section 2. 5 with n = 3, Xz > > X3 shows that these quantities are less desirable than 4(z). It is convenient to view each selection rule as two phases, 1 and 2. Phase 1 is to be used for the first few iterations of IIP until a certain condition is satisfied, and Phase 2 is to be used thereafter. Selection Rule A Phase Al: For k = 0 set Yi(O) = s(-zo), i = 1, 2,.., p, and define scalars ~1,,..., Epp equal to 4(Zo). Whenever 0 < k < p, set Yi(k) = Yi(k - 1), i = 1, 2,.... P Then set Yk(k) = s(-zk_ ) and =k'(z1). When k > p is satisfied, begin using Phase A2. Phase A2: For each k set yi(k) = yi(k - 1), i =, 2,..., p. Then let _ = min{ 1,. P 2,. and let j be the smallest integer in [1, p] for which j = _A. Whenever 4j < 4(Zk ), replace yj(k) by s(-zk 1) and [j by 4(Zkl). Selection Rule B Phase B 1: For k = 0 set Yi(O) = s(-z0), i = 1, 2,..., p, and define scalars [1, z2,..., B p equal to 4(z0). When k > 0 is satisfied, begin using Phase B2. Phase B2: Same as Phase A2. COMMENT 3. 5.1 In Phase A2(B2) there is nothing crucial about the way of handling the possibility of two or more Fig 1 c i < p, being equal to _. Moreover, Bj < (Zk l) may be used as the condition for replacement instead of j < (Zk 1). For selection rules such as A and B in which Yk is a subset of the set {yl(k), yz(k),...,yp(k), s(-zk)}, all k' 0, it is possible to state additional results. Let

43 Sk = {yl(k), yz(k2 k) s(-zk) k= 0, 1, 2,.... (3. 5.1) and consider the determination of Zk+~ in Step 2 of iteration k of II P. p+z i Since Zk+ E Hk, it has a representation: z+ = x Yi where kk+1 i-1 i p+2 i ~~i~~~~~~~i x 0 (i = l, 2..., p + 2). Let Ik be the set of superscripts i, 1 i 0. If in k+1 k - Step 2 of iteration k, k > 0, the integer p + 1 I, then zk k+ = Z' kk k+ik ii) Suppose p is chosen so that q = p + 2- n > 0. Then in Step 2 of iteration k, k > 0, either Zk+l = (which implies zk+ = z' = 0) or there exists a vector x = (x', xZ xp+) such that the corresponding Ik contains at least q elements. Now assume z- / z':- for some k > 0. Then for all k, 0 < k < k: k++ iii) Yk Skc k c QQZk+;Z+); iv) if Y k+lc Sk ) E Q (z. k k k k+i k+i k~i k k+ k+ Zk+ ); v) if Ylc So, S1 contains two distinct points; vi) if Y k+ Sk, Y contains every point in Sk except one, and S contains p distinct k+i k k points (1 < p < p + 1), it follows that there are at least P distinct points in S k+; vii) if Phase Al is used as the first phase in any selection rule and z I z:', S contains p + 1 distinct points. P P Proof: Consider i). Note that k > 0 and let H = a{y,(k-1), Yz(k - 1)...,y (k - 1), s(-zk ), Z k} Since H c H and zk E H, IZk = min lz But p + 1E Ik andY imply zk+ E H and zEH k k k-i k-ii > min I z 1i Thus by part ii) of Theorem 3. 4. 1 = Zk = IZk+ and zEH k zk~I Zk = Zk+i = Z' In ii) suppose Zk+l 0. Part iii) of Theorem 3.2.1 yields Zk+ E aHk. Then part v) of the same theorem and p + 2 - n > 0 imply that a vector x with the desired property exists, Consider iii)o

44 Since z-+l - Z*, part ii) of Theorem 3.4.1 implies IZkI > Iz' > O for k+1 I z I all k, 0 - k $ k + 1. Hence, for 0 < k < k, the hyperplane Q(Zk+;zk+) is defined. For w i 0 let PHk(w) be the support hyperplane of Hk with outward normal w. By part iv) of Theorem 3.2.1 it follows that zk+~ e PHk(-Zk+l) = Q(Zk+;zk+). Thus Hk is contained in the closed halfspace Q (Zk+;zk+l) and iii) is true. In iv) suppose for some k, 0 < k k, s(-z k+l) Q+(z k+l;zk+ ) Then s(-z k+) E Q-(zk+l; z k+) k+i k+i k k+i ki k' and by iii) Q (Zk+l;zk+1 ) also contains Sk and Zk+ Since Yk+ C Sk H Ck+Q k+i; z ) and iz k+ zk+. By part ii) of Theorem k+i k+i k+i k+i k+z 3.4.1 this implies Zk+1 = z* for some k, 0 < k s k, which contradicts Izj > |Iz, 0 < j -<k + 1, and thus established iv). From iii), iv) and + Y1 c So it follows that Y1 E Q-(zl;zl) and s(-zl) E Q (zl;zl). Hence S5 must contain two distinct points. Similarly iii), iv), and Yk+1 Sk imply Yk+ E Q-(Zk+;Zk+ ) and s(-zk+ ) E Q (Z;+) which means s(-z k+) is distinct from the points in Yk+ If Y contains every k+i k+ii point in Sk except one and Sk contains p distinct points, then there are at least p - 1 distinct points in Y +. Since S is the union of Y p k+i' k+1 k+i and s(-z k+), the conclusion in vi) is true. Consider vii). From z i z*, it follows that iii) and iv) hold for 0 < k < p - 1. Thus Yk+lcSk and 0 < -k < p - 1 imply s(-z k+) is distinct from the points in Y k+ For Phase Al Y1 and So contain only s(-zo), Y2 and S1 contain only s(-z0) and s(-zl),...,Y and S contain only s(-zo), s(-zl),..., s(-z ) p p-i p-i Consequently the fact that s(-zk+ ) Y k+ applied successively for k = 0, 1,..., p - 1 yields the result: s(-z0), s(-zl),..., s(-z ) are disp tinct. Since S is the union of these p + 1 points, vii) holds and the proof is complete. COMMENT 3. 5. 2 Observe that if either Selection Rule A or Selection Rule B is used in II P, Yk+ l Sk and Yk+i contains every point in Sk except one for all k > 0. However, Rule A is generally better than Rule B because only p iterations are required with A to ensure that S contains p + 1 distinct points (assuming z i z':). k p

45 Now consider two additional selection rules. Selection Rule C (Assume that p > n and that in Step 2 of every iteration k of I I P a vector x is determined such that the corresponding Ik contains at least p + 2 - n integers. By part ii) of Theorem 3. 5.1 such an x must exist or Zk+ = 0. Terminate IIP in Step 2 of iteration k if Zk+1 = 0 (which implies Zk+1 = Z''' = 0) or if p +1 E Ik (which implies k Zk+l =)) Phase C1: For k = 0 set yi(0) = s(-z), i = 1, 2,...,p, and define scalars j1,,..., p equal to ~L(ZO). When k > 0 and in Step 2 of iteration k - 1 the condition p + 2 E Ik-1 is satisfied, begin using Phase C2. For O < k -< p, set y i(k) yi(k - 1), i = 1, 2,...p. Then set Yk(k) = s(-Zk-) and -k ='(Zk-) For k > p, set Yi(k) = Yi(k - 1), i = 1, 2,..., p. Then let _ = min{irl, I z,..., p} and let j be the smallest integer in [1, p] for which kj = -. Whenever.j- II(zk ), replace y.(k) by s(-z ) and 4j by (Z k-1 ) k-i Phase C2: For 0 < k < p, set y i(k) = yi(k - 1), i = l 2,..9 p Then setyk(k) = s(-zk_l) and k= (zk). For k > p, set Yi(k) = Yi(k - 1), i = 1 2 p Then 1 1 let 4' min mi and let j be the smallest - iEI k-1 ijp+l, p+z integer in {i: i E iIk i p 1, p + 2} for which j = i'. Replace yj(k) by s(-zk ) and.j by4(zk )o — j k-i kCOMMENT 3. 5. 3 For a particular problem it is possible that the condition p + 2 E Ik, k > 0, for entering Phase C2 may never be satisfied. In that case Selection Rules A and C are identical. The computational experience with II P indicates, however, that this

46 condition is satisfied after only a few iterations for a broad class of problems. Selection Rule D (Assume that p - n, that z0 = s(-z ) for some z E K, and that in Step 2 of every iteration k of II P a vector x is -1 determined such that the corresponding Ik contains at least p + 2 - n integers. By part ii) of Theorem 3. 5.1 such an x must exist or Zk+1 = 0. Terminate IIP in Step 2 of iteration k if Zk+ = 0 (which implies Zk+i = Z' = 0) or if p + 1E Ik (which implies zk = Zk+l = Phase D 1: For k = 0 set yj(0) = z0 = s(-z ), set Yi(0) =s(-z0), i = 2 3, 3,...,p, and define ll = -d(Z 1) i = ) (Z~) i = 2, 3,..., p. When k > 0 is satisfied, begin using Phase D2. Phase D2: For 0 < k < p, set Yi(k) = Yi(k - 1), i, 2,...,p. Then set yk+ (k) = s(-zk ) and pk+l (zk ) 1 1 For k > p, set Yi(k) = yi(k - 1), i = 1, 2,. *. p. Then let min i and let j be the smallest k-1 isp+i, p+2 integer in {i: i E Ik_ i # p + 1, p + 2} for which j = i'. Replace y.(k) by s(-zk ) and 4j by j(zk ). COMMENT 3. 5. 4 The assumption p > n is required for Selection Rule C and Selection Rule D so that the set {i: i E Ik i p + 1 p + 2} which occurs in Phases C2 and D2 is not empty. Since p + 2 - n > 2, Ik contains at least 2 integers in [1, p + 2]. Moreover, p+ 1'Ik-i or II P would have terminated in Step 2 of iteration k - 1. COMMENT 3. 5. 5 Observe that if either Selection Rule C or Selection Rule D is used in II P, Yk+ l Sk and Yk+l contains every point in Sk except one for all k > 0. Furthermore, arguments like those used for Selection Rule A show that if z # z*, S (with Rule C) P P and S (with Rule D) contain p + 1 distinct points. p-i

47 THEOREM 3. 5. 2 Consider II P and assume that Selection Rule A C or Selection Rule D is used. Let k be the first k > 0 for which p + 2 E I if Selection Rule C is used and let k = 0 if Selection Rule D is used. Then in Step 2 of iteration k, all k > I, Subproblem 1 can be solved on ASk instead of Hk That is, let k+ satisfy z k k. k+i k+ EAS z k+ = min Iz I and find a vector x = (x',x2,...,xp+) such that k+p ZeAS p+2 i k =k E xy. where y. =yi(k) (i = 1,.*, P), Yp+l s(kz p+z k p+2 i i p+2 P x (i 10 2,...,p +l), xp - 0. i=1 A Proof: Since p + 2 E IA (with Rule C) and yl(k) = z (with Rule A k D), on iteration k Subproblem 1 can certainly be solved on S SA instead k of Ho. Thus zA+ E 4AS. By comment 3. 5. 5 YA c Si' Furthermore, Selection Rules C and D are such that YA contains every point k+i in S. except one which has a coefficient of 0 in the convex combination expression for z. Hence zA E AY and ASA H + so that A k+1 +1 k+ + on iteration k + 1 Subproblem 1 can be solved on AS +1 to yield z E AS + By induction Subproblem 1 can be solved on ASk instead of Hk for all iterations k, k > k. This completes the proof. Note that it is simpler to solve Subproblem 1 on ASk rather than on Hk: the constraint set for the quadratic programming problem is the convex hull of only p + 1 points instead of p + 2. It will henceforth be assumed that whenever Selection Rule C or Selection Rule D is used A in II P, Subproblem 1 is solved on ASk for all k > k. Section 3. 8 contains computational results for II P with Selection Rules A, B and C. These results indicate that it is good to choose p = n and that for a broad class of problems, IIP (with p = n and any of the Selection Rules A, B, C) exhibits much more rapid convergence than BIP. Selection Rule D is identical with Selection Rule C for k > max{p, k}, where k is the first k > O for which p + 2 E Ik when

48 Rule C is used. For all the computations in which Selection Rule C was used, k was observed to be very small. Thus it can be stated that II P with Selection Rule D also converges much more rapidly than BIP. The computational results show, furthermore, that IIP converges at about the same rate with Selection Rule A or Selection Rule C. Thus the decision of which selection rule to use may be based on other considerations. By comment 3. 5. 2 Selection Rule B may be rejected in favor of Selection Rule A. From Theorem 3. 5 2 Selection Rules C and D have A an advantage over Rule A in that for iterations k, k - k, Subproblem 1 can be solved on the convex hull of p + 1 points instead of p + 2. However, the requirement with Rules C and D that Ik contain at least p + 2 - n integers adds complexity to the solution of Subproblem 1 in Step 2 of every iteration k, k > 0 (see Section 3 6). Selection Rule D is most desirable for guaranteeing finite convergence in certain problems (see Section 3.7) and Rules C and D are advantageous for certain optimal control applications (see Chapter 5). 3. 6 Solution of Subproblem 1 Since each iteration of II P requires the solution of Subproblem 1, it is important that methods exist for readily computing its solution. As mentioned in Section 3, 2 Subproblem 1 is a quadratic programming problem on a convex polyhedron constraint set. This is the type of problem that is usually described in the literature [e. g., Al, B2, Hl, V1] under the heading "quadratic programming". However, the computational algorithms which are suggested always begin by assuming the constraint set is described by a set of linear equations and/or inequalities rather than by the points whose convex hull is the constraint polyhedron. Thus to apply the standard quadratic programming techniques directly to Subproblem 1 it is necessary to first determine

49 from these points a description of the constraint set in terms of linear equations and/or inequalities. Such a determination presents nearly insurmountable computational difficulties. There is an alternative method of attacking Subproblem 1 which makes possible the use of the standard algorithms. It is shown in this section that the solution to Subproblem 1 is given by the solution to another quadratic programming problem, Subproblem 2, which has a constraint set described by linear equations and inequalities. Subproblem 2 is solvable by any of the well-known quadratic programming techniques such as those due to Frank and Wolfe [F4], Wolfe [WI], Beale [B4, B5], Houthakker [H4], Hildreth [HZ], Markowitz [M1], and Lemke [L2]. Let Yl, Yz, *- -, m E E be the known points whose convex hull is the convex polyhedron H = A{y1, Yz,..., y } specified in Subproblem 1. Note that since | is the Euclidean norm, an equivalent statement of Subproblem 1 is: find y-* E H such that jY:'2 = =in Jy Y m i m i Each y E H has the representation y -= x y. where Z x = 1, i i= 1 i= x > O (i = l, 2. m) Thus m m vm jy~ I = i EE x yi = E E xYxjY y (3. 6.1) i=1 i=1 j=l If x is the m-vector (x, xZ,.., x ) and D is the m X m symmetric matrix with elements dij = yi yj. then y = x Dx. (3.6. 2) Since Iy 2 _ 0, the quadratic form xo Dx is non-negative definite, a fact which implies it is a convex function of x on E. Consider now the following quadratic programming problem. SUBPROBLEM 2 Given: D, an m X m symmetric non-negative definite matrix, and the constraint set

50 m m~ x = 1 x -> 0 (i = 1, 2,..., m)}. i=1 rain Find: a point x*>' EX such that x*'. Dx* = xE Dx. If D = [dij] = [yi. yj], Subproblem 2 is said to be associated with Subproblem 1. For this case (3.6. 2) implies that minimization of x ~Dx on X is equivalent to minimization of ly I on H. Thus, if x* solves Subproblem 2, y*" is given by m y =, xYi'y (3.6.3) i-= A solution x' E E to Subproblem 2 may have as many as m non, zero elements. However, if Selection Rule C or Selection Rule D is used in Step 1 of II P, it is required that a solution to the Subproblem 2 associated with Subproblem 1 in Step 2 of II P be obtained which has at least m - n zero elements. If no such solution exists, then II P terminates and y:', the solution to Subproblem 1, and z', the solution to BP, both equal 0. Assume m > n and consider Subproblem 1 and the associated Subproblem 2. Suppose a solution x' to Subproblem 2 is obtained for which more than n elements, say n, are nonzero. It will now be shown how a second solution x' to Subproblem 2 which has at least m - n zero elements can be determined, provided it exists. m.m i The point y' is given by y* = xl Y Thus y = x yi reprei=~1 1 i 1=i=1 sents a system of n equations in m unknowns: x=', x2Z x, m The 2 2 T n X m matrix of coefficients [yl, y ym] has rank s n so that at least m - n of the unknowns x'~i may be given arbitrary values. There are m - n zero elements of x' so set the corresponding elements of x' ~~~~~~~~~~1 2 equal to zero. It is required to find n - n of the remaining n~ unknowns which can be set equal to zero.

51 This can be determined by trial-and-error where the number of trials is at most ( — -n)' n'. Suppose n - n of the remaining unknowns are selected arbitrarily and set equal to zero. There results a set of n equations in nunknowns. These equations may present computational difficulties because a solution does not necessarily exist, and even if it exists it may not be unique. If there is a solution such that the resultingx' satisfies i xz 1, xzi > 0 (i = 1, 2,...,m), then this is the de2 g x i=1 2 sired x4. Otherwise another set of n - n unknowns are selected and the procedure is repeated. If the solution x;' indicates that y*r 1 0, there is another way of m m obtaining y*' in the form y* = lcry. where i al. =1 i > 0 i=1 1 1 j=1 1 1 (i = 1, 2,...,m), and at least m - n of the -., 1' i m, are zero. 1 Parts iii) and v) of Theorem 3. 2.1 and y* 1 0 imply that this form exists and that y* E aH. The method is to consider the m = M; (m - n) n. convex polyhedra Al, za,. _ formed by the convex hull of n points chosen from {y1,yz,.*.y }m. For each A. (j = 1, Z o..m) a Subproblem 1 and its associated Subproblem 2 (D an n X n matrix here) can be solved to yield a point y" E A. such that j I min Y If y min (1: j min < m) satisfies y i| =main m Ym: I then part v) of Theorem j min -)sis Y j 3.1. 2 and y*: E AH imply that y* = y. Thus y*- can be written as a convex combination of n points from {Yl,Y2,..., Ym} which is the desired form. Of the two methods just described the former is simpler because it requires a Subproblem 2 to be solved just once instead of mH + 1 times. It may be possible to avoid entirely the possibility of trial-anderror or solution of subproblems on Aj (j = 1, 2,...T ). This is the case if a quadratic programming technique can be found for Subproblem 2 which yields a solution with a maximum number of zero elements. Considering the simple nature of the constraint set X it is likely that this would not be difficult to do.

52 3.7 A Finite Convergence Theorem The following theorem gives a sufficient condition for I I P to exhibit finite convergence. THEOREM 3. 7.1 Let s(.) be an arbitrary contact function of the set K specified in BP, choose p > n, and consider IIP with Selection Rule D. Assume that K is a convex polyhedron and that the range of s(y) for y E En is a finite set of points. Then the sequence {Zk} generated by II P converges in a finite number of iterations. If K is a convex polygon in E (n = 2), then the range of s(y) for y E En is a finite set of points. This conclusion is not necessarily true, however, for a convex polyhedron K in En, n > 3. For example, in E3 an entire edge of K could lie in the range of s(y), yE E3. Nevertheless, for many convex polyhedra K such as those which arise in optimization problems for linear sampled-data systems, it may be possible to choose a contact function s(y) of K whose range for y E E is a finite set of points. Proof of Theorem 3.7.1: Consider first the following result. LEMMA 3.7.1 If K is a convex polyhedron in En and s(-) is an arbitrary contact function of K, then every extreme point of K lies in the range of s(y), y E En Proof: Let x be any extreme point of K. It can be shown that there exists a support hyperplane of K that contacts K in the single point x. Therefore, if y is the outward normal to this support hyperplane, x = s(y) and the lemma is proved. Now let dS, sz,.9.s denote the points in the range of s(y) y e E and define S = {sl,sZ,...,s}. By Lemma 3.7.1 every extreme point of K E S. This and Sc K imply K = AS. Thus z'EE AS. Suppose that the convergence of II P is not finite. From

53 Comment 3.5. 5 and part vi) of Theorem 3. 5.1, Sk, k > p - 1, contains p + 1 distinct contact points when Selection Rule D is used. Each of these points is in S. If 2 p + 1. Since p > n, this implies 2 > n+l. Let A1, A2, 2,, ( n+l)(2 -n-l) denote the convex polyhedra formed by the convex hull of n + 1 points chosen from S. By part vii) of Theorem 3.1. 2, AS = A. Thus z*' E AS implies z:- E A. - j=j - j for at least one j, i' j < 2. Let 6.(j = 1, 2,... ) be the distance of minI each set Aj from z*; that is, 6j zE - z* At least one 6., 1 - j -, is 0 and 6S = 0 if and only if z:- E Aj Since Sk, k > p - 1, contains p + 1 distinct points from S and p > n, there is at least one a., 1 < j < 2, such that A c AS Furthermore, part vii) of Theorem 3.1. 2 implies U Aj = ASk. Let 6(k) A-cASk mn in min be defined for k > p - 1 by: 6(k) -= A.SS That is, 6(k) = A.cK J zEASk lj k z -z*. By Theorem 3.5. 2 it follows that zk E ASk = Hk, k > 0, when Selection Rule D is used. Hence Izk - z*|> 6(k) > 0, all k > p - 1. Part iii) of Theorem 3.4.1 implies izk - z'l -* 0 as k - oo, so 6(k) -+ 0 as k -* 00o But 6(k) has 2 values at most. Consequently there exists k > p - 1 for which 6(k) = 0. This implies z*' E AS_ and thus z- = z k k+ contradicting the supposition that convergence is not finite. This completes the proof. COMMENT 3.7.1 If there exists k > 0 such that the condition p + 2 E Ik is satisfied in Step 2 of iteration k, then Theorem 3.7. 1 holds for II P with Selection Rule C. COMMENT 3. 7. 2 If 0 K, the assumption in Theorem 3. 7. 1 may be replaced by the following two less restrictive assumptions:

54 i) there exists E > 0 such that for z E K and I z - z*- < E, s(-z) E S(-z*), ii) S(-z-') is a convex polyhedron and there are only a finite number of contact points contained in S(-z':-') which lie in the range of s(-z), z E Kr Then the sequence {Zk} generated by IIP with Selection Rule D converges in a finite number of iterations. The proof of this result is based on the fact that for 0 q K, 1(Zk) -+ IZ:": i By supposing convergence is not finite it can be shown that there exists k' such that for k > k':- Skc S(-z*'), zkE S(-z*:), and Sk contains p + 1 distinct points. Then in much the same manner as in the proof of Theorem 3. 7. 1 a contradiction can be established. Figure 3. 7. 1 illustrates a set K in EZ which satisfies assumption ii) but not i). By inspection it is clear that unless zo E K is such that zo = wz':-, w > 0, convergence is not finite. O = Origin S(-z ) I K Figure 3.7.1 A Kc EZ which does not satisfy assumption i) of Comment 3o7.2.

55 3. 8 Numerical Results for II P This section contains results of extensive computations with II P for Example 3 of Section 2. 5. As in Section 2. 6 v = 1, z0 = (1, 0,., 0), and the extent of K is increased beyond z1 = 2v so that (2. 5.18) is valid. Data is presented for n = 2, 3,.. 6; various p, X2, X,3.. X6, z0; and Selection Rules A, B, and C. The quadratic programming technique of Frank and Wolfe [F4] was used to solve the Subproblem 2 which occurs on each iteration of II P. This technique, which exhibits finite convergence, makes possible the use of the simplex method of linear programming [D1] as a subroutine. Briefly, the technique is as follows. If D is the m X m matrix in Subproblem 2, define the (m + 1) X (2m + 2) matrix D by 1 1...1 0 0 0 0...0 -1 1 D= -1 1 (3.8.1) -2D I -l 1 Furthermore, if r is the column vector (x,,,,, v), where x, v E Em and s,, E El, let r+ denote the row vector (v, 0, 0, x). Note that r+r = 2x E v and consider SUBPROBLEM 3 Given: D as in (3. 8.1). Find: r E E such that Dr = (1, 0,..., 0), r > 0, and r+r is a minimum. By employing the Kuhn-Tucker conditions [K3] Frank and Wolfe show that if a solution x:- to Subproblem 2 exists then there is a solution r to Subproblem 3 for which r+r = 0. Moreover, if such an r is found, its first m components are a solution to Subproblem 2. Subproblem 3 is solved in the following two phases (a feasible vector r E E satisfies the constraints Dr = (1, 0,...,0), r > 0; a

56 basic vector r E E contains no more than m + 1 nonzero elements): Phase I A basic feasible vector ro is found with which to begin Phase II (the simplex method with one artificial vector [HI] is easy to use here). Phase II There are a finite number of iterations in this phase. On the th h iteration a feasible vector wh and a basic feasible vector rh-1 are available (on the first iteration use wl = r0). Employ the simplex method to minimize the linear form w+tr, obtaining the sequence of basic feasible vectors { rj, j =1, Z... r~ - rh such that w+ 0 wr > wr2 >. h-1l h.h h-. Stop at the first rJ such that either (rJ)+rJ = 0 or w+rJ <ww h. h 2 h h' If the first condition is satisfied, rj solves Subproblem 3. w+ (wh - rh) Otherwise let rh = r, h = min h, 1}, and w = ah + -; and (rh wh )(r Wh h+l h + h(rh - Wh); then repeat Phase II using Wh+ and rh' It is interesting to note that in all the computations with II P, every solution x* to Subproblem 2 obtained by Frank and Wolfe's technique had at least m - n zero elements. Thus the difficulties (trial-anderror or solution of additional subproblems) mentioned in the latter part of Section 3. 6 were not encountered. Tables 3.8.1 through 3.8.9 (for IIP Selection Rule A and p = n) correspond directly to the nine tables in Section 2. 6 (for BIP). Furthermore, Figures 3.8.1, 3.8. 3, and 3.8.4 are similar to Figures 2. 6.1, 2. 6,3, and 2. 6.4. These tables and figures show the marked improvement of IIP over BIP. This improvement has three facets: (1) much more rapid convergence, (2) very little dependence on the parameter -i k v, and (3) no noticeable influence of the initial point z0.

57 The behavior of lZkj-iz'., I k - z, Iz:i l-lzkY (zk), and I s(-zk) - z'': I shown in Figures 3. 8. 3 and 3.8.4 is typical of that observed in all computations with II P0 As with BIP it can be stated that z I- Izk lY(zk) decreases most rapidly, followed in order by IZkl- IZk I I zk - z* I, and is(-zk) - z': I. Some other quantities related to the sequence {Zk} are displayed in Tables 3.8.10, 3.8.11, and 3 8.12 (for n = p = 3 and II P Selection Rule A), and Table 3. 8.13 gives more data on IZk - I z-*' for different Nz/X3 ratios. As with BIP the parameter 2 /3 has very little effect. It is clear from Table 3. 8.14 that Selection Rule A is superior to Selection Rule B. Note that with Selection Rule B it may happen that a large number of iterations are required to satisfy I Zk - l zl: 1, after which convergence is quite rapid. This shows the significance of Sk, k > p, containing p + 1 distinct points (see Comment 3. 5. 2). Table 3.8.15 differs from Table 3.8.1 in that p equals 1 instead of 2. Tables 3.8.16 (for n = 2) and 3.8.17 (for n = 3) compare results for a variety of p values. This data and Table 3. 8. 20 (for n = 4) show that convergence is good for p = n and little improvement is obtained for p > n. The desirability of choosing p = n is also indicated by the actual computing time required to satisfy I Zk - z-'< 10- However, as mentioned in Section 2 6 the evaluation of a contact function in optimal control applications is the most time-consuming part of the iterative procedure. Thus the number of iterations to satisfy given error criteria is a better measure of the performance of the procedure. Results for n = p = 4, 5, 6 and IIP Selection Rule A are given in Tables 3.8.18, 3.8.19, 3.8.21, 3.8.22, and 3.8.23. Roughly speaking, the rate of decrease of Zk - I z is dependent on n alone. The number of iterations per decade after a few initial iterations, is approximately 2 for n = 2, 4 for n = 3, 6 for n = 4, 9 for n = 5, and 13 for n = 6.

58 Finally, Tables 3. 8. 24 and 3. 8. 25 give results for IIP Selection Rule C. Observe that the rate of convergence in all cases is about the same as that for II P Selection Rule A. The sequences { I Zk - z: }, {Iz:"I-lIzkIY(zk)}, { Is(-zk) - z':'I}, as well as {IZkl-IZ":} I behave similarly for II P Selection Rules A and C.

59 1 A ~ 1 —- B.01 10I I k 0 10 20 30 40 50 Figure 3.8.1!zk/I IZj for n = 2, p = 2, zO =(6, 2): A) X2 = 100, B) Xz = 1000; IIP Selection Rule A.

60 1 A B -4 10 10 k 0 10 20 30 40 50 Figure 3.8.2 Results for n = 2, p = 2, zo = (6, 2), X2 - 100: A) i<k z ly(zi), B) i zk -z:'; IIP Selection Rule A.

61 Table 3. 8. 1 Number of iterations to satisfy I Zk Z- J <E; n = 2, p = 2, z0 = (6, 2), IIP Selection Rule A. k2 1. 1.01 10 10-4 10- 10-5 6 10 3 4 6 7 10 11 13 100 2 7 9 11 12 14 15 200 4 5 7 9 10 10 12 500 3 5 7 7 10 12 14 1000 5 6 9 9 12 13 15 104 8 11 14 16 18 18 20 105 8 9 10 12 14 16 18 106 10 11 14 15 17 18 20 Table 3. 8. 2 Number of iterations to satisfy I zk - z ) $ E; n = 2, p = 2, z0 = (6, 2), IIP Selection Rule A. [I 1.1.01 10o3 10-4 10-5 i0-6 10 3 5 6 9 12 >13 100 4 7 9 11 14 > 15 200 4 5 8 9 > 12 500 3 5 7 9 11 >14 1000 5 6 9 10 12 14 > 15 104 8 11 14 16 18 >20 105 8 9 10 12 14 16 18 106 10 11 14 15 17 18 >20

62 Table 3. 8. 3 Number of iterations to satisfy Iz:i:- zk y(zk) ke n - 2, p -, zo - (6, 2), IIP Selection Rule A. E~z X Z i 1. 1.01 10-3 10-4 1o 10 0 3 3 6 8 9 12 100 0 1 8 9 11 11 14 200 0 4 4 7 9 9 9 500 0 2 6 6 6 11 11 1000 0 5 5 8 8 12 14 104 0 4 4 15 17 17 17 10o 0 7 9 9 13 15 15 106 0 10 10 13 16 16 19 Table 3. 8. 4 Number of iterations to satisfy ts(-,Zk) - z:' I s E; n - 2, p = 2, zo = (6, 2), IIP Selection Rule A. 2l 1 1.01 1 0-' 108 10 10 10 3 6 9 12 > 13 100 8 11 14 >15 200 7 9 9 >12 500 6 11 14 >14 1000 8 12 14 >15 104 17 17 20 >20 10 15 18 >18 106 19 >2O

63 Table 3. 8. 5 Number of iterations to satisfy I Zk -:: -< n = 2, = 2 2, X = 100, IIP Selection!Rule A. 1 2 0 1 1 0 -3 1 0-4 - -6 zoB. ZO2 1 1.01 10 1 10 10 6 1.5 2 4 8 11 14 16 16 6 2 2 7 9 11 12 14 15 6 2. 5 3 6 7 10 11 12 13 6 3 2 7 10 12 14 15 17 6 3. 5 2 6 8 9 11 13 14 6 4 3 5 6 8 10 11 13 6 4.5 4 8 10 11 12 14 15 6 5 3 7 8 10 12 13 14 2 6 2 6 9 12 12 13 15 3 5.57 2 3 7 8 10 10 12 4 4.9 3 7 9 10 12 14 16 4.5 4. 44 2 6 8 1 1 13 16 18 5 3. 87 2 7 9 10 10 12 13 5.5 3. 12 3 6 9 1 1 1 3 14 15 6. 2 1. 25 3 5 7 8 10 12 14 Note: For ZO = (6, 2) and the last seven cases in the table, zo I. - 40.

64.01 10-3 10- A B 100 10 20 30 40 50 Figure 3.8.3 Results for n = 3, p = 3, zo = (6, 2, 2), X2 = 100, X3 = 10: A) IJz,-' i] k Iy(zl), B) IJzk1 - z, C) Ik Z*; I IP Selection RuleA. For k 9, IZk - z*I ~ IZk-IZI*.

65 10.01 10 - 10 4 0 10 20 30 40 50 Figure 3. 8.4 Results for n - 3, p = 3, zo = (6, 2, 2), X2 = 100, 3 - 10: A) Iz I -lzkY(zk), B)!s(-zk) - zI; II P Selection Rule A.o

66 Table 3. 8. 6 Number of iterations to satisfy | zk - z E < E n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. C E ( -3 1 4 -5 -6 N2 X3 I. 1. 1 001 10 1 1 1 2 3 3 3 4 4 10 10 3 3 6 7 9 10 12 100 100 4 5 7 9 11 13 13 1000 1000 4 6 8 10 11 12 12 100 1 2 4 8 10 12 14 17 100 10 6 11 14 17 21 26 32 100 50 7 11 15 18 22 26 30 100 90 6 9 12 17 23 28 3 1 1000 1 6 7 8 11 13 15 17 1000 10 7 9 13 1 20 23 26 29 1000 100 8 14 17 23 26 29 32 Table 3. 8. 7 Number of iterations to satisfy zk - z* <: n = 3, p = 3, zo = (6, 2, 2), IIP Selection Rule A. E e 3 -4 - -|X2,1. 1. 01 10-3 10 10 10 1 1 1 2 3 4 4 >4 10 10 3 5 6 8 >12 100 100 4 5 8 10 12 13 1000 1000 4 6 8 10 11 >12 100 1 2 6 10 16 >17 100 10 6 11 16 20 31 >32 100 50 7 11 15 21 27 >30 100 90 6 9 13 23 3i 31 1000 1 6 7 12 12 16 >17 1000 10 7 10 17 24 28 >29 1000 100 8 14 18 25 30 32

67 Table 3o 8. 8 Number of iterations to satisfy j z>- -Zk I(z) E n = 3, p = 3, zo = (6, 2, 2), IIP Selection Rule A. XX 1 3 1 1.01 1- 10 10 -6 _ _ _1 o 1 10 1 1 0 2 2 2 2 3 3 10 10 0 2 2 6 8 8 11 100 100 0 4 4 8 9 10 12 1000 1000 0 5 7 8 9 11 11 100 1 0 1 9 100 10 0 8 12 16 20 20 31 100 50 0 9 13 15 19 24 27 100 90 0 7 11 13 18 25 27 1000 1 0O 5 7 10 12 14 16 1000 10 0 7 12 14 22 25 28 1000 100 0 12 16 16 25 27 30 Table 3, 8. 9 Number of iterations to satisfy I s(-Zk) - z*: < E n = 3, p = 3, zo = (6, 2, 2), II P Selection Rule A. L3E3 0 01 |i.4 10i 1 0 [ 1 |- 1 0 2 3 3 4 4 >4 10 10 2 6 8 >2 100 100 4 10 12 >13 1000 1000 9 11 11 1 1 > 12 _ 100 1 1 11 13 >17 __ 100 10 10 18 31 >32 100 50 13 19 27 >30 100 90 11 1 8 >31 1000 1000 7 14 >17 1000 10 12 22 >29 1000 100 16 27 >32.., _

68 Table 3. 8. 10 Number of iterations to satisfy I Zk lIzk Y(zk) <) n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule A. _ _ Xi 1.1 1.01 109 l06 10 I 10 1 1 1 2 i 3 3 3 41 4 10 10 3 3 6 8 9 11 >12 100 100 4 5 7 9 11 >13 1000 1000 4 6 8 10 1 1 12 > 12 100 1 2 6 8 10 12 14 17 100 10 7 11 14 17 2! 27 32 100 50 8 11 15 19 22 27 t 30 100 90 6 9 13 18 i 24 28 >31 1000 1 6 7'8 11 13 1,6 1 7 1000 10 7 9 14 20 1 24 27 29. 1000 100 9 14 i 17 24 27 30 32 Table 3. 8. 11 Number of iterations to satisfy - yZk) zk E n = 3, p = 3, zo = (6, 2, 2), II P Selection Rule A. _____~ X, 1.1 0i 10 1- 0 10- 10-6 1 1 2 3 4 >4 10 10 3 6 9 > 12 100 100 5 8 11 >13 1000 1000 5 8 11 12 100 1 3 9 12.1 100 10 8 14 24 32 100 50 9 15 23 >30 100 90 7 13 25 > 3 1 1000 1 6 9 14 > 17 i t 1000 10 8 14 24 29 1 1000 100 10 17 27 32,

69 Table 3. 8. 12 Number of iterations to satisfy s(-z k)- |I z' E; n = 3, p = 3, z0 = (6, 2, 2), IIP Selection Rule Ao., E2 0 -3 0-4 10-' -., 1 01 10 1 1 1 1 1 2 3 4 4 >4 10 10 2 2 6 8 8 11 12 100 100 4 8 9 10 12 12 >13 1000 1000 8 9 11 11 11 11 11 100 I 1 9 11 11 13 13 >17 100 10 10 14 16 20 20 31 31 100 50 12 15 19 24 27 > 30 100 90 9 13 18 25 27 >31 1000 1 7 12 14 16 >17 1000 10 12 12 22 22 28 29 1000 100 16 23 27 30 30 >32 Table 3. 8. 13 Number of iterations to satisfy I Zkl- |z:'I <; n = 3, p = 3, zo = (6, 2, 2) II P Selection Rule A. r E [ I i I F... I. I.7 Wn6 xkz | 1. 1 0. 0 i 11 10 104 1 5 8 12 14 16 17 20 104 10 8 11 14.19 21 26 31!04 100 13 17 19 23 3 1 34 40.............................. 104 1000 14 19 23 26 30 34 37 104 5000 13 18 21 30 33 36 40 104 9000 13 15 21 25 28 37 41 10 1 8 8 12 14 17 19 23 105 10 9 13 16 20 25 29 33 10S 100 10 19 23 25 32 35 40 10 loo1000 12 15 18 22 27 33 38 L10 104 15 17 23 28 32 33 38

70 Table 3. 8. 14 Number of iterations to satisfy I k z E; n = 3, p = 3, z, = (6 2, 23 IIP Selection Rule B. E -A -5 -3 A 1. 1.01 10. 10 10 10 k 1 1 1 2 3 3 3 4 4 4 10 10 3 3 6 7 9 10 12 12 100 100 4 5 7 9 11 13 13 13 1000 1000 6 9 10 12 14 14 15 12 100 1 2 4 8 10 12 14 17 17 100 10 6 11 14 17 21 26 32. 32 100 50 8 11 16 20 23 27 30 30 100 90 9 12 19 26 30 33 36 31 1000 1 6 7 8 11 13 1i5 7 17 1000 10 10 12 17 20 23 27 30 29 1000 100 23 39 42 46 47 54 57 32 104 1 8 11 14 17 19 2 23 20 104 10 10 16 18 25 28 30 34 31 104 100 54 57 61 63 67 71 73 40 104 1000 26 29 33 37 41 44 51 37 10_ 4 5000 52 63 67 70 75 83 88 40 104 9000 155 160 164 173 177 180 185 41 105 1000 165 168 17 1 174 179 s185 9 1 38 105 1 104 55 59 62 65 68 76 79 38 A -- Note: k = the first k for which Izkl - J Z 10 with II P Selection Rule Ao

71 Table 3. 8. 15 Number of iterations to satisfy | zk -z <E n = 2, p = 1, z0 = (6, 2). liP Selection Rule A. 2 [E 1 o 1.0 01 10 10 10 3 5 8 2 16 19 22 100 2 7 12 i5 18 23 27 200 5 9 13 16 21 24 27 500 6 9 12 16 19 23 27 1000 6 1 0 15 19 23 28 32 104 8 11 14 20 22 26 30 105 10 15 19 21 25 29 32 10 10 14 18 21 25 28 31 Table 3 8. 16 Number of iterations to satisfy I zk -z z* I sn = 2, z0 = (6, 2), II P Selection Rule A and BIP. -3 -4 -, 6 \201. 1 01 10 10 10 t 0 2 9 25 42 60 60 75 25.5 1 2 7 12 15 18 23 27 29. 1 100 2 2 7 9 11 12 14 15 17.2 3 2 7 9 11 12 14 1 5 18 9 0 36 113 126 134 150 158 241 82, 6 1 6 9 12 16 19 23 27 28. 5 500 2 3 5 7 7 10 12 14 15.2 3 3 5 8 10 1! 13 15 18 4 4 3 5 8 10 12 13 15 21.5 Note: 1) p = O corresponds to BIP. 2) t = actual computing time (seconds) for IBM 7090.

72 Table 3. 8. 17 Number of iterations to satisfy I Zkl- Z-' - < n = 3, z0 = (5, 4, 2), X2 = 1000, X3 = 100, II P Selection Rule A and BIP. p I I 0 1. 0 1. 0 10- 1 t 0 53 77 138 202 241 298 366 125. 3 1 30 52 97 134 171 249 293 319. 1 2 9 15 27 36 46 52 58 7 5. 9 3 9 16 23 26 29 33 36 50. 1 4 8 14 17 24 27 32 37 77.8 5 8 13 16 23 26 31 35 90. 3 Notes: 1) p = 0 corresponds to BIP. 2) t = actual computing time (seconds) for IBM 7090.

73 Table 3. 8. 18 Number of iterations to satisfy I Zk |- I I -< E; n = 4, p = 4, zo = (6, 2, 2, 1), I P Selection Rule A. kZ kB3 \. 1. 1. 01 10- 0-43 10-0 110- 100 10 1 5 9 15 19 23 27 29 100 10 10 6 10 15 20 25 28 33 100 50 10 7 14 23 29 36 43 49 100 50 50 7 10 13 17 20 23 27 100 90 10 7 13 19 28 34 41 48 100 90 30 9 14 21 28 37 42 49 100 90 50 8 14 19 23 3 1 37 43 100 90 70 9 17 23 30 35 43 52 100 100 50 6 9 13 16 20 22 26 1000 50 10 9 15 21 26 32 38 43 1000 100 1 7 9 14 18 22 25 29 1000 100 10 8 14 21 25 29 40 46 1000 500 100 10 15 26 3 1 38 42 48 1000 700 100 13 19 24 30 38 45 51 1000 900 500 13 19 30 41 47 57 63 1000 950 900 8 19 23 42 47 53 58 1000 995 990 4 15 21 3 1 44 49 54 104 5000 1000 18 28 35 44 48 59 64 105 5 X 104 104 16 23 33 39 46 52 57

74 Table 3. 8& 19 Number of iterations to satisfy I zk - z'' < E; n = 4, p 4, zo = (62, 2, 21) i I P Selection Rule A. k3 X4 1 1.01110 10 10 10 100 10 1 5 10 18 22 >29 100 10 10 6 10 19 27 32 >33 100 50 10 7 14 26 [40 >49 100 50 50 7 1 10 16 19 > 27 100 90 10 7 13 24 34 46 >48 100 90 30 9 15 25 t 36 >49 100 90 50 8 14 21 35 >43 100 90 70 9 17 26 i37 51 >52 100 100 50 6 9 15 19 > 26 1000 50 10 9 15 2i 32 >43 1000 100 1 7 9 16 21 27 >29 1000 100 10 8 15 23 28 > 46 1000 500 100 10 15 27 33 41 >48 1000 700 100 13 19 25 34 48 >51 1000 900 500 13 20 30 44 56 > 63 1000 950 900 8 19 24 42 52 > 58 1000 99 5 990 4 16 21 38 49 > 54 104 5000 1000 18 28 3 5 44 51 >64 i0O 5 X 104 iO 16 23 33 39 46 55 >57

75 Table 3. 8. 20 Number of iterations to satisfy Zk I - z': E n = 4, z0 = (6, 1,, 2 2), Xk = 1000, \3 = 500, X4 =100, II P Selection Rule A. 1 1.01 1 10- 10-4 10- 10 t 1 23 67 107 152 185 226 264 331. 1 2 19 38 79 101 127 141 158 211. 4 3 15 27 36 46 62 74 78 123. 1 4 12 18 25 29 37 48 53 86. 6 5 11 17 22 30 39 49 57 113.3 Note: t = actual computing time (seconds) for IBM 7090.

76 Table 3. 8. 21 Number of iterations to satisfy I Zk- - I n - 5, p = 5, z = (5, 3,1,1. 8, 2. 6), II IP Selection Rule A. -3 -4 -5 -6 x2 X 3 4 1 1.01 10 10 10 10 100 50 30 10 8 19 30 39 42 54 64 100 70 50 10 9 20 29 39 48 58 66 100 90 50 10 10 20 3 1 39 47 58 67 1000 700 500 100 19 27 33 42 53 68 I 82 1000 900 500 100 17 25 34 43 51 72 79 Table 3. 8. 22 Number of iterations to satisfy | zk - z':- c, n = 5, p = 5, zo =(5,3,1,1. 8, 2. 6) _I P Selection Rule A. E 2 X3 X4 X5 1.1 01. 1 3 10-4 10-5 1 0-' 100 50 30 10 8 21 34 52 > 64 100 70 50 10 9 20 33 48 >66 100 90 50 10 10 20 34 54 > 67 _,,....... -..._..... ~_ 1000 700 500 100 20 29 39 62 80 > 82 1000 900 500 100 17 26 37 58 >79

77 Table 3. 8. 23 Number of iterations to satisfy error criteria; n = 6, p = 6, z0 = (4, 3, 2. 6, 2. 6,1L 8,1. 8), II P Selection Rule Ac Error E Criterion A1.1 01 10 01 Al 10 24 34 53 64 79 93 zk|-| z*:- E -.......,2 1 Z 26 38 45 63 76 9 1 A1 10 24 41 73 > 93 IZk-Z E f A2 13 25 39 60 >91 AlZ-kY(Z) A 0 23 23 41 61 73 90 kzr k zki Az 0 20 31 39 60 74 85 Al 23 61 87 >93 Is(-zk) - z.A2 33 60 90 >91 A1 10 24 37 56 67 82 >93 lZk- IIzk ) I E zk -< Az 14 26 39 47 65 79 >91 Al 16 39 69 > 93 A2 16 39 68 >91 A1 23 4 1 59 66.82 92 >93 Is(-zk)I-IZ I -E A 26 39 39 68 5 >9 1 Note: 11: X2 = 100, X3 =70, X4 = 50, \5 =30, X6 =10. A2: x2 = 100, x3 = 90, x4 = 70, \s = 50s X6 = 10

78 Table 3. 8. 24 Number of iterations to satisfy Zk l-Z' t; II P Selection Rule C. n=p k X3 X4 15 1. 1. 10 10-410-5 10 k k'. _, ____......._ j -. 2 100 2 7 9 11 12 14 15 15 2 3 100 10 5 13 16 20 24 28 31 32 3 3 100 50 7 14 i 17 21 24 29 31 30 3 3 100 90 6 12 15 19 22 26 29 31 3 3 1000100 9 15 18 23 27 30 32 32 7 4 100 50 10 8 15 125 30 38 45 51 49 3 4 100 90 10 7 14 21 28 34 42 47 48 3 4 100 90 50 9 18 24 31 35 43 48 43 3 4 100 90 70 7 13 25 33 38 47 51 52 3 4 10001500 100 10 17 28 35 40 44 50 48 6 5 0O 70 50 10 8 17 26 36 42 54 63 66 5 5 i100 90 80 70 9 16 21 37 48 58 69 70 4 6 o100 90 70 50 1 10 11 24 37 49 62 78 92 91 5 6 11000 90 70 50 10 15 29 45 60 73 84 94 95 7 Notes: 1) k = the first k for which IZk| -z': 10 with II P Selection Rule A. 2) k' = the first k for which Phase C2 is used. 3) For n = 2, zo = (6, 2); for n = 3, zo = (6, 2, 2); for n = 4, zo = (6, 2, 2,1); for n = 5, zo (5, 3,l,l. 8, 2. 6); for n = 6 Zo (4, 3, 2. 6, 2. 6, 1. 8, 1. 8).

79 Table 3. 8. 25 Number of iterations to satisfy Zk J - Iz' < E II P Selection Rule C. -3 -4 -5 -6 n p 1.1.01 10 10 10 10 t 3 2 7 15 26 34 43 52 61 72. 4 3 3 5 13 16 20 24 28 31 44. 3 3 4 5 14 18 21 23 27 i 30 61. 7 3 5 5 12 15 20 24 26 29 80.6o 4 3 i 9 22 37 51 61 72 88 137. 1 4 4 11 7 13 25 33 38 47 51 88. 3 4 5 7 16 25 32 39 44 50 134. 3 4 6 1 7 13 23 31 38 44 50 f1866,O 4 7 7 14 22 31 39 45 51 1284, 4 Notes: 1) t = actual computing time (seconds) for IBM 7090. 2) For n = 3: z = (6, 2, 2), X2 = 100, X3 = 10. 3) For n = 4: zo = (6, 2, 2, 1), 2 =- 100, X3 = 90, X4 = 70.

CHAPTER 4 THE GENERAL ITERATIVE PROCEDURE GIP By following some of the ideas introduced by Fadden [F1], it is possible to state a general problem which has application to a wide variety of optimal control problems. This general problem GP is formulated in this chapter and a method for solving it, called the general iterative procedure GIP, is described and shown to converge, Each iteration of GIP involves the minimization problem BP, which can be solved by II P or BIP. 4. 1 The General Problem GP Let 2 = [0,, A],,> 0, be a compact interval in El. For w E Q consider sets K(w) c E which are compact and convex. Also let K(o) be continuous on 2, i. e., for every E > 0 there is a 6 = 6(E, w) such that K() c K(w + -) + N(0;E ) and K(w + a) c K(w) + N(0;E) whenever JUI < 6 wherew, w + WE e Before stating GP the support and contact functions for K(w), w E 2, are introduced. Let rj(w, y) =zEK(w) z y denote the support function of K(w). Since K(w) is compact, r(w, y) is defined on ~ X Eo As in Section 2. 2, for fixed w E 2 I(w, ~) is a convex continuous funcn tion on E. The following result also holds. THEOREM 4. 1. 1 Given that the sets K(w) are compact, convex, and continuous on 2. The support function rj(w, y) is continuous as a function of w and y on ~ X En Proof (due to Fadden [F1] o It is necessary to show that for every E > 0 there exists a 6 = 6(E,, y) such that r(, y) - r(W, y) I < E whenever I(, y) - (w, y) I < 6 where w, E and y, y E En, Now 80

81 n(@ Yy) -( y)J = n(y -( Y) -n( y)l < rl(~, y -rl(o, y' -+1 ri(, Y) -r(w, Y)Jo The continuity of (q(w, y) for fixed w implies that given e > 0 there is a 61 = 61(E,W,y) such that Ir(Wo, y) -n(w,y) < 2 whenever Y - y < 6, That 6t(1, y) is a continuous function of w for fixed y follows from the continuity of K(w). For every T > 0 there is a 8 = 6(w, w) such that K(w) c K(-) + N(O; ) and K(s) cK(w) + N(O;,) whenever | - I | <. Consider fixed y and note that the support function of N(0O; ) is y|. Hence, i(w, y-) < r(i, y) + Tj|Y and r(U, y-) < 3(w), y-) + Ely. Thus r(3, y) - (wo, y-) < e JyJ whenever J~ - | < 6. Clearly for every E > 0 there exist a 62 = 62(E,, y) and a 63 > 0 such that ir(u, y-)-( -(,y- Y7) < - for all y E N(y;63) whenever J - w J < 62. Therefore, since I(W, y) - (w, y)J I|3 - w I + Jy - y J, the desired property is established if 6 = 62 + min{6l, 63}. For c E 2 let P(w, y), y i 0, be the hyperplane {x E E: x y = I(wy y)} Since z o y -< r(w, y) for all z e K(w) and P(w, y) nK(w) is not empty, P(w, y) is the support hyperplane of K(ow) with outward normal y. For each y i O the set S(w, y) = P(wy)nK(w) is called the contact set of K(w) and its elements are called contact points of K(w). It follows that S(w, y) is not empty, S(w, y) c aK(), S(w, Xy) = S(w, y) for X > 0. A function s(w, y) defined on 2 X En is a contact function of K(w) if s(o, y) E S(w, y), y # 0, and s(w, 0) E K(W). Thus for w E 2- s(w, o) is bounded; s(o, y) = s((w, Xy), X > 0; r(wo, y) - s(w, y) - y. If for every w E Q2 and every y E En there is a method for determining a point x(o, y) E K(w) such that x(w, y) ~ y = K() z - y rl(, y), then it is said that a contact function of K(w) is available. This availability is essential to the computing procedure GIP presented in the next section. Consider now the general problem:

82 GP Let 2 = [O, A], W > 0, be a compact interval in E1 and for w E 2 let sets K(C) c En be compact, convex, and continuous in w. Assume that there exists o*i E 2 such that 0 4 K(@), 0 _ w < 0<, and 0 E K(0w0*). F'ind: c'. Note that because K(X) is continuous it follows that 0 E aK('o*). The relationship of GP to problems in optimal control is discussed in Chapter 5. The computational problem posed by Fadden [F1] differs from GP in that the sets K(w) are required to be strictly convex. Furthermore, in addition to wo" Fadden seeks to find the outward normal to a support hyperplane of K(w^c) which contacts K(^-*) at the point 0. 4. 2 The General Iterative Procedure GIP In this section the iterative procedure for computing the solution to GP is described. First some preliminary definitions and results are presented. Let Q2 and K(w) be as specified in GP. For w E i2 define functions z (w) and p(w) = IJ(cw)I such that Z(w) E K(X) and p(w) = z(w) = mmin Izi zEK(w) The compactness of K(c) and continuity of I z| ensure the existence of these functions. Geometrically,'Z(w) is the point in K(49) closest to the origin and p(w) is the distance of'Z(w) from 0. It follows from Theorem 2. 2.1 that: z'(w) is unique; p(w) = 0 if and only if 0 E K(W); for p(o) > 0, z(w) E aK(w); for p(w) > 0, z ='(Co) if and only if z E P(o, -z) n K(w) - S(w, -z). Moreover, the following result holds. THEOREM 4. 2.1 The functions p(w) and'(w) are continuous on Proof: Consider p(o) first and let E > 0 be given. It is required to find 6, = 61(E, o) such that J p(w) - p(cc + )l <E whenever )~, < 6, where a, X + X E 2. Since K(X) is continuous, there exists 6 - 6(E, )

83 such that K(W) K(o + W) + N(OE) and K(w + r) c K(w) + N(O E) whenever IW1 < 6f It will be shown that 61 = 6. The first inclusion and z(w) E K(X) imply that there exist z1 E K(w + &3) and xl E N(OE ) such that z(w) = zl + xl. This yields Izl I = I(w) - xl z - I|(W)I+ixi |<( i(4+) + E. Similarly the second inclusion and z(f + W) E K(w + a-) imply that there exist z2 E K(XC) and x2 E N(O;E) such that Z( +X) = Z2 + X2, 1Z2 < Iz(w +&3) +E. But Iz(o + ) I I Jzl and Iz(w2) - i zz Therefore, jz(o +WC) < IZ(w)l +E and Iz(w)I 0 be given and let 6 = 6(E, ) be chosen as in the preceding paragraph, The fact that K(X) is bounded for all w E 2 implies that there is an a > 0 such that p(w) < a, w E Q. Since z () - z ( + w) |- p(c) + p(eo + -) < 2r, E may be restricted to the interval 0 < E < Zc. Two cases will be considered, corresponding to p(w) < E and p(w) > E Suppose first of all that p(w) < E For | <K 6(E,) I P(') - P(W +-) < E which yields p(w + a) < p() + E < 2 E Then I (w) - (i v + ) | - p(w) + p(w + CJ) < 3E Now consider p( c) > E and I) | < 6(<: w4) Let Ng = N(O. p(w) + E) Q1 = Q(r c'(W).' - z~(W)), Q Q(- G(W) - E p ()() - Z()) and Q+ Q +(z(W) - Ep-(X) z(W); - M(X)). As illustrated in Figure 4. Z. I Q, is the support hyperplane of K(w) at z (w) with outward normal -'z(w). It can be shown that the parallel hyperplane Q2 is the support hyperplane of K(w) +N(O;E) at Z (w) - Ep1l(w)'z(w). Since Ip(G) - p(0 + )J < E, z(w + X) EN1o In addition, z(w +X) E K(w) + N(O;E) implies'(w +~) E Qz2 the open half-space bounded by Qz with outward normal -z(i)- Thus if x is any point in Qz n AN1 it follows that |z(c) - z(zc+C-)l < |i(z) - x | Furthermore,

84 /, / / / K~w) p(W) + - E -1, Z(W)-Ep (Z)Z(O )' /K(w) +N(0f) / Q2 Q1 Figure 4. 2, 1 Notation for the case p(w) > E in the proof of Theorem 4.2.1 z(w) - XI = E + [p(W) +] _ [p(W) - E = E2 + 4p(w)E = E [E + 4p(w)]. Since p(w) < ac and E < 2a, the quantity in the brackets is less than 6ac. 1 Thus Iz(W) -'z(w + )I < (6aE) 2. Clearly the above results imply Ij(w) - z(w+&)I <E if [j] < min{6( -, w), 6( -, W)}. Therefore ~(w) is continuous on Q2 and the proof is complete.

85 Now let w be an element of 2 and s(w, ) a specific contact function of K(w). Consider -2 y(w, z)= Iz z z S (w, -Z), Iz| > 0, 9 s((, -Z) > 0 (4. 2.1) = 0, z = o or zi >, z s(w, -z) < 0. Thus y(o, ) is a function defined on K(o). It follows from Theorem 2. 3.1 that for w E ~ and z E K(wo) 0 -< y(w, z) < 1; if 0 E K(w), y(O, z) = 0 for all z E K(w); if 0 q K(O), y(w, z) = 1 if and only if z = z'(co); for fixed w E ~2, y(w, z) is continuous in z. The following results also hold. THEOREM 4. 2. 2 Let Q2, K(w), and': be as specified in GP. Then for 0 - w< wd - and z E K(@): i) Iz IY(O, Z) < P(W); ii) rj(o:', -y) > 0, y EE n; iii) if y(wo, z) >, r( z) I ) < 0. Proof: Part i) follows from the definition of p(w) and part i) of Theorem 3. 3. 1. Inequality (2. 4. 4) applied to the set K(w*) yields s(o*, -y) - y'z ('*) y, ye E. Since z (cW*) = 0, i(b:- -y) = -y. s(:.-, -y) > 0 which proves ii). Consider iii). The fact that 0 ~ K(w) for 0 -<w<* implies Iz l > 0, all z E K(c). This, y(w, z) > 0, and (4. 2. 1) yield z ~ s(w, -z) > 0. So iii) follows from (w, -z) = -z ~ s(, -z) and (4. 2. 1). Now GIP can be stated. The General Iterative Procedure GIP Let ~2, K(o), and o" be as described in GP. For w E ~2 let (w,. ) be the support function of K(w), let s(w, ) be an arbitrary contact function of K(X), and define y(w, ~ ) by (4. 2. 1). Set wo = 0 and i = 0. The following two steps constitute one iteration, called iteration i, of GIP. Step 1 Consider.i fixed, 0 - i. < a*, and the corresponding set K( i). Apply II P or BIP to the minimization of z, z E K(wi). If W. < *, IIP or BIP may be continued (Theorems 3o 4. 1 and 2. 4. 1) until a point z. E K(oi) is obtained such that y(i~, zi)> 0i, where 0 < ~ < Oi < 1 and 0, Oi are preselected numbers. When 1 1~~~

86 this happens proceed to Step 2. If w. = Io>*, I P or BIP will not produce a z E K(ci) satisfying y(i., z) O0i but will generate (Theorems 3.4.1 and 2.4.1) a sequence of points in K(wi) = K(w*) which converges to 0. Step 2 Consider z. fixed and look at ri(, -zi) for X E [i., X*]. Since 1 1 1 y(Wi' zi) - 0, r(i, -Zi) < 0 (Theorem 4.2. 2). This, rii(':, -zi) 0~ and the continuity of q(w, z) (Theorem 4. 1.1) imply there exists X E (wi.,'*] such that r(, -z.) = 0. Let X increase from x. and ~~~1 1~i 1 let Wi+. be the first w for which (w, -z.) = 0o Then increase i by 1 one and return to Step 1. It should be noted that in Step 1 with i > 0 the point s(wi. -z, ) lying on the hyperplane Q(0O-z. ) can be used to initiate II P or BIP. For i - 0 the 1-1 point s(w0., x) for any x can be used. If the condition x.i = act in the latter part of Step 1 occurs, then the solution to GP has been obtained in a finite number of iterations of GIP. If this condition does not occur, then GIP generates infinite sequences {wi} and {zi} where cw. E [0, w) and zi E K(i) To show convergence of wi to a'- and zi to 0 it is necessary to make an additional assumption on 11(Wa, Y) q r(b-9 Y rj(w, y), namely, that the difference quotient be bounded a b from above for all a,a b E f, Wa f b bounded y E En This is treated in the next section. 4.3 Convergence Theorem for GIP THEOREM 4. 3.1 Consider GIP and assume the condition to. - * in the latter part of Step 1 does not occur for i < ooo Furthermore, assume ](wa, Y) - 7(wbP Y) that is bounded from above for all co, c E, wX, _ W a bb a b a b bounded y E Eno Then GIP generates sequences {i } and {z.} which for i > O and i -oo satisfy: i) o. e [0, o*) and z c K(w.); ii) {X } is strictly 1 1i i 1

87 increasing (wi < ai+ ) and w. -* w*; iii) z. -- 0. 1 j+1 1 1 Proof: The first half of i) follows from o0 E [0, co*) and the definition of ci+i in Step 2. Since all points z generated by II P or BIP applied to the set K(w ) are in K(.i), z. E K(wi). i i 1 i Consider ii). The property oi < Wi+ follows by the remarks in Step 2. Thus {Wi} is a strictly increasing sequence on a compact interval [0, X;:'] and must have a limit, say'. To complete the proof of ii) it suffices to show that cn* - E < X < w:', where E > 0 is chosen arbitrarily. Assume the contrary, i. e., s < * - E. Let supremum 71 ay) - (b yY) a b ab a a b -yE U K(w) which exists by hypothesis. Clearly v > 0. At the start of Step 2 Xi < W', zilJ > 0, and Y(oi, Zi) > 0. By part iii) of Theorem 4. 2. 2,,(zi1 - yzi) =0. But -ri(z, -Zi) = 0. Hence, i+1 V Jzi Y(i.9 z.) I V i 0 > 0. (4. 3.1) Now consider w E [0, o - E] and the function Z(w). Since z(w) is continuous on [0, we& - El (Theorem 4. 2. 1), it takes on its minimum, say Z. on this interval. Then. = m En [ mn in z]l> 0 Con- min' min 05i''e ZEKM sequently Izi I > |zein I for all i: 0 which by (4. 3. 1) implies.+l -. > v Iz0in I = constant >0 for all i ~ 0. This contradicts the fact that Wi has a limit on [0, a":< - E] and thus ii) is true. Now consider iii). From xi -) the continuity of'z(w), and Z(wo) = 0, it follows that Z(Oi) 0. Part i) of Theorem 4. 2. 2 and y(e., zi) > 0 yield 0 fz. r In z iY(oi, zi) i (hi). Since the right side of this inequality converges to 0, zi J - 0 and thus zi i - 0. That is, z. -0, which completes the proof. 1

88 Note that it is not possible to say that zi+1 I Izi | (see Figure 4. 3. 1). Q(O; ( i z) 0 Z. \iP~a~iX Zi 1 1 P(Wi0- zi) P(W+I'Zi+li+l z.). Figure 4. 3. 1 A configuration for which I Z+ > Z. in GIP. 1 ~a 1 Also observe that if on an iteration of II P or BIP within Step 1 of GIP a point z E K(w.) is obtained for which z | = 0, then x. = o:. However, 1 1 if I z < E (arbitrary E > 0) is used as a stopping condition for GIP, then xi may or may not be close to w:. In certain applications where w represents fuel, cost, effort, time, etc. (see Chapter 5) it may be satisfactory to obtain a, O 0 < w <d:, provided that a point z E K( ) sufficiently

89 near to the origin is found. In such cases J z| < E provides a stopping condition for GIP. Consider the following: COROLLARY 4. 3. 1 Let the hypotheses of Theorem 4. 3. 1 be satisfied, consider {wi} and {zi} generated by GIP, and let E > 0 be given. Then: i) ifD'X <evl Zi < (Ve E)2; ii) if I Zi < E and for P(Wa) - P(wb) wa, b E [0, w*], Wa Wb - < 0, it follows that'* - w. a - b -1 -1 a b -1 < -cr E. Proof: Inequality (4. 3. 1), w. < X. < w', and w -. < E imply -v zi zi ) <E which proves i).. Consider ) a w = w < Wi it 1<1 1<1 1 1 p(w'') - p(w ) follows: -. Since p(w'*:) = O and w. < w-' -p (wi) -(: Wi) andw - W - P(wi) But P(wi) I zil < E so - W. <~ Ecompleting the proof. Part ii) of the corollary has importance if the scalar C- > 0 exists and can be calculated. For in that case whenever the stopping condition |z < E is satisfied for a point z E K(wi) generated by IIP or BIP within Step 1 of GIP, ii) provides a measure of the error in w. Moreover, if z < w-E is used as the stopping condition for GIP, then w- - w.i < E. 1

CHAPTER 5 APPLICATION OF BIP, II P, AND GIP TO OPTIMAL CONTROL PROBLEMS 5.1 Some Optimal Control Problems Solvable by GIP In this section a class of optimal control problems is formulated using some ideas of Fadden [Fl]. Many problems of this class, and other optimal control problems as well, can be solved by GIP. Some can be solved by BIP or IIP alone. Six specific problems which illustrate a variety of optimization objectives are stated. The application of the iterative procedures to these six problems is discussed later in the chapter. Let U be a compact set in Er and let O = [0,t ], t > 0, be a compact interval in E1. A measurable function u( ) defined on O whose range is in U is said to be an admissible control. Consider dynamical systems of the form x(t) = A(t) x(t) + f(u(t), t), x(0) specified, ( 5.11) where x is the m-dimensional state vector, x is its time derivative, x(0) is the initial state, u(.) is an r-dimensional vector control function defined on 0, A(.) is an m X m matrix function defined and continuous on O, f(.,.) is an m-dimensional vector function defined and continuous on U X 0. For every admissible control u( ) there is an absolutely continuous solution function xu( )(t) = x (t), which satisfies (5.1.1) almost everywhere on O. An admissible u(.) is said to generate the states x( t). Let there be given target sets W(t) c E which are compact, convex, and continuous on 0. An admissible control u(o) is said to transfer x from x(O) to W(t) in time t E 0 if x (t) E W(t).o 90

91 Let there also be given a cost functional of the form t Jt(U) = [a() x (i) +f(u(o), cr)]dcr+h~(x (t)) (5.1.2) 0 where a( ) is a continuous function from e to Em, f (.,) is a continuous function from U X O to E1, and h~(.) is a convex (continuous) function from Em to El. For every admissible control u(.) Jt(u) yields a particular cost function, which if evaluated at t E 0 results in a number called the cost at time t for this admissible control. Consider now the following class of optimal control problems: Fixed Terminal Time Problems Let T > 0 be a fixed point in O. Find an admissible control u(-~), if one exists, such that: a) u-'(.) transfers x from x(O) to W(T) in time T; b) the cost at time T for u*"(.) is less than or equal to the cost at time T for any admissible control u( ) which transfers x from x(O) to W(T) in time T. (For these problems Jt(u) = JT(u) = J(u).) Free Terminal Time Problems Find an admissible control u*( ~) and an optimal time t:' E O, if these exist, such that: a) u:-'(-) transfers x from x(O) to W(t':') in time t"; b) the cost at time t' for u-(-) is less than or equal to the cost at time t for any admissible control u( ) and time t C 0 for which x(t) c W(t). u Any admissible control u'*(-) which satisfies conditions a) and b) in either a fixed or free terminal time problem is said to be an optimal control. The above formulation is also appropriate if there are constraints on the control function of the form s(i),rd i 513o jc(A -~rd Mi=,2,2 513

92 where the ( (,.) are continuous functions from U X O to E' and M are closed intervals in E'. By introducing I additional differential equations.(m~i) i.(m+i) x = (u(t), t), x (0) = 0, i =1, 2,...., (5.1.4) and letting (x(t), m(m+l) m+l and letting (x( t), x(l) (t),..., x( ) (t)) E E be the new state vector, these constraints can be incorporated into a system description of the form of equation (5.1.1). The target sets in this case are W(t) X M1 X -X M. Since there is no loss of generality, it is assumed in the sequel that no constraints of this type are present. Many optimal control problems of the class described above can be solved by GIP. Furthermore, it is often possible to use GIP on problems of an extended class which includes other cost functionals (e.g., equation (5.1.9)) and/or target sets which are closed but not bounded (e. g., half-spaces or linear manifolds of dimension > 1 which may or may not vary with t). Some of these problems can be solved by BIP or I I P alone. The following six example problems will be used later in the chapter to illustrate the application of the iterative procedures. The first five have a fixed terminal time T; the sixth has a free terminal time. Problem 1 (Minimum-Error Regulator) J(u) = IXU(T)I (5.1.5) For this problem, first posed by Ho [H3], the target set W(T) is equal to the whole space Em Problem 2 (Generalized Minimum-Error Regulator) J(u) = xu(T) Gx (T) + gxu(T) (5. 1.6) where G is a symmetric non-negative definite m>x m matrix and g is an

93 m-vector in the range of G. This problem, treated by Gilbert [G1], also has the target set W(T) = Em Problem 3 (Convex Function Minimum-Error Regulator) J(u) = h (x (T)), (5.1.7) where h~(.) is a convex (continuous) function from Em to E1 such that {x' h~(x) < constant} is compact, and W(T) = Em, Problem 4 (Minimum-Fuel Terminal Control) T J(u) = ka(t). x (t) + f0(u(t), t)]dt, (5.1. 8) 0 where a(-) is a continuous function from [0, T] to Em and f0(,) is a continuous function from U X [0, T] to E1. Problem 5 (Minimum Effort) j~u). 0 sT i ( l iu) tT lu(t)l. (5.1.9) For this problem, which was introduced by Neustadt [N3], the target set W(T) is a single point w(T) and f(u, t) = B(t)u where B(~) is an m X r matrix function continuous on [0, T]o. Problem 6 (Minimum Time) Jt(u) = t, (5.1.10) which corresponds to (5.1. 2) with a(.) = 0, f0(,) = 1 and h~( ) = 0. There are no further restrictions on W(t) for this problem. 5. 2 The Reachable Set; Determination of a Contact Function The problem of computing an optimal control can be approached by considering the set of all possible solutions of (5.1.1) using admis - sible controls, This set and some of its important properties are

94 discussed in this section. In addition, two closely-related sets are introduced. Consider the system (5.1.1) and let t e. The set R(t) = {x x = x (t), u(-) admissible} (5.2.1) is called the reachable set. It is the set of all possible states in Em which can be reached from x(O) at time t using an admissible control. Since the solution to (5.1.1) can be written as t x (t) = (t) [x(O) + i(o) f(u(cr), -) do] (5. 2. 2) 0 where P(t) is the m X m matrix solution of i - A(t)d), 4(0) = I, it is clear that R(t) = {~(t) [x(O)+ 5 l'()f(u(), -)d cr: u(*) admissible}. (5. 2.3) 0 Neustadt has shown [N2] that for each t E O the set R(t) is: (1) compact; (2) convex. From (5.2. 3) and the continuity of f(.,') and ~ (-) it is easy to prove that: (3) R(t) is continuous on O. Furthermore, for fixed t E ) it will be shown below that: (4) a contact function of the set R(t) is available. The properties (1), (2), (3), (4) are the essential features which permit application of the computing procedures of Chapters 2, 3, and 4 to optimal control problems. It should be observed that for all but the most elementary systems (5.1.1) the complexity of (5. 2.3) prohibits explicit calculation of the boundaries of R(t). Hence, the most satisfactory computing algorithms are those based solely on properties (1), (2), (3), (4). Let T E be given and consider the determination of a contact function s (T, Y), Y E M, of R(T). For y a specified m-vector it follows from (5.2.2) that

95 T Y X (T) =Y. (T) X(O) + 0 ([~(T) () ] *y)',r)d-, (5. 2. 4) 0 where the prime denotes matrix transpose. But the m-vector 4l(t, y) = [(T) (t)]'y, defined on [0, T] X E, is the solution of the adjoint differential equation @(t) = -A'(t) @l(t), 41(T) = y. (5.2. 5) Hence, xER() x = (0) x(O) + maxU [p(o-,y). f(v, x)]do(. (5. 2. 6) 0 It can be shown that there exists an admissible control u(t, y) defined on [0, T] X Em such that for almost all t E [0, T]: max 4(t, y) - f(u(t, y), t) = U 4(t, y) - f(v, t). (5.2. 7) Then from (5.2. 6) it is clear that y XU(ty) (T) > Y XU(t) ( T) for every admissible control u(t). It follows that a contact function of R(T) is SR( T, Y) x Xu ) (T). (5. 2 8) This result agrees with the well-known fact that boundary points of the reachable set must "satisfy" the Pontryagin maximum principle. In most practical problems it is not difficult to obtain a function u(t, y) which satisfies (5. 2.7). Consider, for example, the case where f(u, t) = B(t)u, B(t) is an m X r matrix function continuous on 70, T1-] and U is the unit hypercube {u' l u I < 1, i = 1, 2,.... r}. Notice that (5. 2. 7) may not uniquely define u(t, y) almost everywhere in [0, T]. For in the present example, (5. 2. 7) yields u(t, y) = sgn B'(t) L4(t, y) where the ith component of sgn v (i = 1, 2,...r) is 1 for vi > 0, is -1 for v < 0, and is arbitrary if v = 0o This is of no concern, however, since even if a component of B'(t) 4(t, y) is identically zero on [0, T] (io e., LaSalle's

96 normality condition [L1] is not satisfied) different choices for u(t, y) will at most lead only to different contact functions of R(T). Previous computational procedures [B1, El, Fl, F2, N1, N3, N4] have required assumptions which correspond to a unique determination of u(t, y) by (5.2o7). Such "unique maximum" assumptions imply strict convexity of R( T) o Computer evaluation of SR(T y) for a specified y E Em entails three steps: i) Evaluation of LP(t, y) by solving (5. 2. 5) backwards from t = T to t = 0; ii) Determination of u(t, y) from ](t, y) by (5. 2. 7); iii) Solution of (5.1.1) with u(t) = u(t, y) from t = 0 to t = T. As observed by Fancher [F3], storage of the function u(t, y) can be avoided by first obtaining i(0, y) from the backward integration of (5.2. 5) and then integrating (5. 2. 5) and (5.1.1) together from t = O to t = T while obtaining u(t, y) from (5. 2.7). Thus each evaluation of a contact function involves the sequential solution of two differential equations. This situation can be handled effectively by a digital or hybrid computer [G2]. It is convenient to introduce here two sets C(t) and R(t) which are closely related to R(t) and have application to a number of optimal control problems. Let t E e and define t C(t) - (Cl ) f(u(fo),) do: u(') admissible}. (5. 2.9) Note that C(t) C Em is related to R(t) by the continuous, non-singular affine transformation C(t) = m (t)R(t) - x(0). (5.2.10)

97 Consequently, C(t) is compact, convex, and continuous in t for all t c 0. Furthermore, a contact function s (t, y), y E E, of C(t) is available, as is shown in the following paragraph. Let T E O and y e Em be specified. If c (T) denotes the point in u C(T) generated by the admissible control u(t), then: T y' Ca(T) =S (['(-()]y) f(u(cr), c) dc (5. 2,11) 0 and T max - I cEC(T) y C = maVaU ([l ()] y) f(v, c) dcr (5. 2.12) 0 It can be shown that there exists an admissible control u(t, y) defined on [0, T] X Em such that for almost all t E [0, T]: ([m ~(t)] y) f(u(t, y), t) = maU ([ (t)] y). f(v,t). (5. 2.13) Then if T sC(T y) = ('-() f(u(ar-y), o)d, (5. 2.14) 0 it is clear that y - sC(T, y) A Ya C u(T) for every admissible control u(t). Thus C(T, y) is a contact function of C(T). An alternate formula for evaluating sc( T, y) is C (T, y) - (T)x (T) - X(O), (5.,. 15) where u(t, y) is obtained from (5.2.13). This follows directly from (5.2.10) and (5.2.14). In certain optimal control problems such as Problem 4 of Section 5. 1 it is useful to introduce a set R(t). For these problems there is a cost function x~(t) E E1 satisfying almost everywhere on O the differential equation

98 xo(t) = a(t) x (t) + f(U(t), t), x~(0) = 0, (5.16) where a(-) and fo(, *) are as in (5.1. 2). Then if x = (x~ x) denotes a m+i point in E, (5.1.1) and (5.2.16) are equivalent to the equation x(t) = A(t) x(t) + f(u(t), t), x(0) - (0, x(0)), (5. 2.17) where f = (f~, f) and A(t) is given by ~ al'(t) a (t)' a(t) A(t) = 0 A _ -., A(t) O I 0i Let x (t) be an absolutely continuous solution of (5. 217) corresponding -u to the admissible control u( ). Then for t E O define R(t) = {x: x = x (t), u(.) admissible. (5. 2.18) If 1(t) is the (m + 1) X (m + 1) matrix solution of b = A(t) T, T(0) = I, y is a specified (m + 1)-vector, and d(t, y), defined on [0, ] X E m+, is the solution of,(t) = -A (t) 4i(t), 41(T) = y (5z. 219) then equations (5. 2. 2) through (5. 2. 8) may be rewritten with sub-bars on each x, c, f, R, y, A, and dL. Thus R(t) is compact, convex, and continuous on 0. Moreover, for -E (O, a contact function of R(T) is sR(T, y) = X (t )(T) ( (5. 2. 20) where u(t, y) satisfies for almost all t E [0, T]: (t, y) fu(t y), t) = max 1(t, y)' f( u( t, y_), t) vu( tt, y) - f(v, t). (5.2.21)

99 5.3 Solution of Problems 1 and 2 by BIP or II P Alone Problems 1 and 2 of Section 5o1 can be solved by GIP or by BIP or II P alone. The method of solution by BIP or II P alone is treated in this section. In Problem 1 let K = R(T), where T > 0 is the fixed terminal time. The optimization objective is: find a point z*' E K and an admissible control u*4(.) such that Iz* m = mzJ z and z* = x,.(T). Since R(T) zEK U-;. is compact, convex, and has an available contact function (equation (5.2. 8)), this is clearly a variant of BP and BIP or I IP may be used to obtain a point'JE K such that z - z* < E (arbitrary E > 0)o There is no difficulty in initiating the iterative procedures: choose zo = x (T) E K, where uo(-) is any admissible control. (In IIP with Selection Rule D, an arbitrary admissible control say u (.), can be used to generate z E K.) Furthermore, for all k > O the contact point s(-zk) = SR(T, -Zk) is generated by an admissible control u(t, -zk) satisfying (5o 2.7). There remains the issue of finding admissible controls Uk(.) and u(-) which will generate zk and zo For BIP or II P (with Selection Rules A, B, C, or D) zk, k > 0, has a representation:z = s0-.(-zi) + -1 zO k-1 i= - where -> 2 0 (i = -1,,..., k - 1) and ] i. = 1. Suppose that for almost 1 i=_1 1 all t E [0, T] the sets f(U, t) = {f(v, t): v E U} are convex. Then for k> 0 there exists an admissible control uk( ) such that almost everywhere in [0, T] f(uk(t),t) = Z c f(u(t, -zi) t) t - f(u0(t), t), which implies by (. hk = iO - 1 (5.1.1) that zk = Xuk(t)(T). For example, if f(u, t) = B(t)u + b(t) and U is ~* ~ k-1 convex, then uk(t) = cr.iu(t, -Zi) + C uo(t). If the sets f (U, t) are not convex for almost all t E [0, T] and additional approximation process, the construction of a chattering control, is necessary [B3]o This chattering process yields an admissible control Uk(') which generates a point arbitrarily close to zk. As a brief illustration, for f(u, t) = B(t)u + b(t) and U not convex, the process is as follows.

100 k-i Consider u(t) =.I -. v.(t) where the c. are as above and the adj=-1 1 1 1 missible functions v.(t) are: v (t) = uo(t), v.(t) = u(t, -z), i = 0,1,... k -1. A 1 -1 1 k-i i Suppose u(t) does not have its range in U. Since. =u 1 and cr. ~ 0 1=-i 1 1 (i = -1, 0,.... k-l), the intervals r = (r - 1) T+ T -- T, (r - 1)T + E CT. rs i- s -i T = q T, r= 19 2...,q, s = -1,0...,k -l, are mutually disjoint and cover (0, T]. Thus uk(t) = vs(t), t EF, r = 1, 29.. q, s =-1,0 O,..., k -1, is an admissible control defined on (0, T] which "chatters among v (t), vo(t),... vk (t) the fractional time spent on vi(t) being.. It is not difficult to show that there exists a q (dependent on E, A(t), B(t), b(t), T, and U) such that IxA(T) - x_ (T)I <E for q > qE. Thus a sufU uk ficiently fine division of (0,T] implies IxA(T) - x_ (T)I < I' (arbitrary A u uk E > 0). It is clear that chattering can be omitted on any subinterval of A (0. T] where u(t) has range in U. By Theorems Z. 4. 1 and 3. 4. 1 zk -- z-* for BIP and IIP. Thus for k sufficiently large zk will satisfy specified stopping criteria. The approximate optimal state is z = z k and it is generated by the admissible control u( ) = Uk('). Observe that uk( *) is required only on the final iteration of the iterative procedures. In some situations, however, it may be desirable to compute Uk( *) for each k > 0 in a recursive manner, using uk ( * ) and those u(, -z.i) 0 i - k - 1 which occur on iteration k - 1. If K is strictly convex, then s(-zk) will satisfy specified stopping criteria for k sufficiently large. Hence' and i"'( -) can be taken to be s(-zk) and u(, -z k). This avoids the necessity for finding uk( ), but examples show that an increased number of iterations are required. There are several instances when it is possible to simplify the finding of an admissible uk( ) which generates zk. Consider II P with Selection Rule D. Then Zk, k > 0, has a representationmzk-= Ai s(-zi) iEI 1

101 A where I is a set of p + 1 or fewer integers from [0, k - 1], i 0(i E) and LA ci =1. It follows that Uk(*) can be constructed from u(t, -z.), i E I. Since it may occur that p + 1<< k - 1 where k denotes the iteration on which the stopping criteria are satisfied, the use of Selection Rule D simplifies the construction of u( *) = Uk( *). For example, f(u, t) = B(t)u + b(t) and U convex imply u(t) = uk(t) - iu(t' -zi). This same simplification is achieved in II P with Selection Rule C if on some iteration prior to the final one, Phase CZ is entered. Furthermore, even with II P (Selection Rules A or B) or BIP it may be possible to represent zk as E uis(-zi) where I is a proper subset of {0, 1,..., k - 1} iE I (note: I may contain more than p + 1 integers), c-. - 0 (i E I), and Cari.=lo iE~I This is the case if for some j < k z, has a zero coefficient in the convex combination expression for z j+ Now consider Problem 2. Let q(x), x E E B, be the quadratic form q(x) = x- Gx + g. x, (5.3.1) where G and g are as in (5.1. 6). It is easy to show that q(x) can be written as q(x) = Gx + gl + q (5 3 2) where G is an 2 X m matrix, 2 = rank G, G = G'G, g =!(G G') Gg (that is, g = 2GI g), and qO - - gl 2 min is g Gg and = xE q(x) Thus if K is defined by K = {z: z = Gx + g, x E R(T)} (5.3.3) it follows that K is compact, convex, and xER(T) q(x) = ztK 2 +z + q0. (5.3.4) Consequently Problem 2 can be stated as a variant of BP: find an x* E R(T) and an admissible control u*( I) such that

102 x*E{x: Gx + g= z::, x R(T)}, z,2 = minxI2 a =x 2T). Note that x' is not necessarily unique. A contact function s(y), y E En, of K is st(y) = GsR(T, G'y) + g. (5.3 5) This follows from max max zeK Y' z = xeR(T) Y. (Gx + g) = sR(T, G'y) ~ G'y+y- g = y (GsR(T, G'y) + g). (5.3.6) Since for every x E R(T) there is a corresponding point Gx + g = z E K, BIP or IIP can be extended by simple substitution to Problem 2. For BIP Gilbert [G1] has written out the complete details. The procedures yield a point x~E R(T) such that x - x' | < E (arbitrary e > 0), x>:* an optimal state. An admissible control u(' ) which generates x is determined in a way similar to that described in Problem 1. 5. 4 The Sets K(wc) = X(&.) - Y(w) Before applying GIP to Problems 1 through 6, some results on the difference of two sets are stated. Let Q2 = [0, i], A > 0, be a compact interval in E1 and consider sets X(w), Y(w)cEn which are compact, convex, and continuous on i2. Define K(o) = X(o) - Y(w) {z: z = x -y, x E X(w), y E Y(w)} (5.4.1) It is easy to show that the sets K(c) are compact, convex, and continuous on 12. Furthermore, the support function rl(w, *) of K(w) can be written in terms of the support functions X(w,') and rly(w. ) of X(w) and Y(w) as follows: (arbitrary yE E n)

103 max max (W, y-)- zEK(o) z y =xEX(W (x-Y) y yEY(w) max max xEX(c). y + yEY(o) y (-y IX-(W ) y7 + ry("' m~1. (5.4. 2) From (5.4. Z2) it follows that if sx(o, ) and sy(w,.) are contact functions of X(w) and Y(w), then a contact function of K(w) is s(W, y-) =s (X(, y) - Sy(W, -y), y En (5.4. Clearly 0 4 K(w) if and only if sets X(w) and Y(o) do not intersect. 5. 5 Solution of Problems 1, 2, and 3 by GIP Let X(w), X E 2 = [O, c], be the fixed set R(T) in each of Problems max 1, 2, and 3. A reasonable choice for A' is xER(T)f(x) where,(x) equals lx I, x Gx + g x, h~(x) in Problems 1, 2, 3 respectively. Define Y(w), w E 2, as follows: Y(o) = {x E Em: Ix i' } for Problem 1, (5.5.1) Y(w) = {x E E x Gx + g x < w} for Problem 2, (5.5. 2) Y(M) = {x E E: h~(x) s aw} for Problem 3. (5.5.3) Since Ixl, x. Gx + g. x, and h~(x) are convex functions on Em it follows that these sets Y(w) are convex and continuous on ~2. Clearly Y(w) in (5. 5.1) is compact, and if it is assumed that G is non-singular then Y(W) in (5.5. 2) is compact. The compactness of {x: h~(x) - constant} was assumed in (5.1.7). (Methods exist for treating convex functions h (.) which do not satisfy this assumption.) Thus if K(w) = X(W) - Y(w), K(w) is compact, convex, and continuous on 2 for Problems 1, 2, 3. Consequently each of these problems can be stated as a variant of GP: find'&- E 2, a point xv E R(T), and an admissible control u*( ) such that

104 0 4 K(w), 0 < X < <w*, 0 E K(co*"), x'-' E Y(w-o), and x*" = x (T). The quantity o.' is the minimum cost in each problem. To apply GIP it is necessary that a contact function s(w, y) of K(i) be available. Since X(X) = R(T), sx(w, y) = sR(T, y) and thus (5. 4. 3) yields s(w, y) if sy(o, y) can be determined. In Problem 1 sy(w, y) = oly l y for lYl > 0; sy(w, y) = arbitrary x E Y(w) for lY = 0. Finding sy(w, y) in Problems 2 and 3 is more difficult. Consider Problem 2. Since the gradient of x. Gx + g. x is ZGx + g, sy(w, y) must satisfy (for lyl > 0): 2Gsy(W, y) + g = y, > (5> 4) and sy(w, y) Gsy(w, y) + g sy(CO, y) = o. (5.5.5) These are a set of m + 1 equations in m + 1 unknowns: sy(CO, y) and At. Any solution for sy(w, y) will suffice but it may not be easy to determine. Similarly in Problem 3 if the gradiant of h~(x) exists, it follows that sy(w, y) must satisfy (for ly I > 0): grad h0(s(w, y)) = [Ly, 0 > O (5.5.6) and h~(s(w, y)) =w. (5.5.7) For convergence of GIP it is required that the support function rl(yas Y) - l(ob, Y) t](w, y) of K(o) satisfy: a bounded from above for all wo - b a b Woa, b E 2, w a # Wb' bounded by y e Em. Since X(w) does not vary with w in Problems 1, 2, 3, (5.4. 2) implies that this requirement is satisfy(Ca( Y) - ~y (Wb, Y) fied if Y is bounded from above (wa, Wb, y as before). a b For Problem 1 ly(w, y) = y. sy(wo, y) = wXY |, jY > O, so the desired condition is true.

105 Assuming that a contact function s(w, y) is available and that q(w, y) satisifes the requirement of the preceding paragraph, GIP may be used to generate coE [0, o:'-] and a point' E K(~') such that z = x - y, t_-.; Irae.m in x E X(w), y E Y(w~), where IzI < E and Y(: - y < E (arbitrary yEY( w E > 0). For Problems 1 and 2 {x~- x':> < E, where x:' is an optimal state (x': unique for Problem 1). An admissible control u( ) which will generate x can be found by a method similar to that described in Section 5. 3 since x can be written as a convex combination of an initial point xO E R(T) and a certain number of contact points of R(T). 5. 6 Solution of Problems 4, 5, and 6 by GIP A maximum O Consider Problem 4 and choose c -(x0 x)ER(t) x~, where R(t) is defined by (5. 2. 18). The function f0(., ) in (5. 1. 8) may be chosen so that the lower boundary of the closed interval 2 is 0. For w E 2 = [0, A], let X(w) = R(T) and Y(w) = {((w, w) w E W(T)}. Since W(T) is compact and convex, Y(w) has these properties for each C E 02. The fact that Y(c) is continuous on Q2 is obvious. Hence K(w) = X(w) - Y(w) is compact, convex, and continuous on ~2 and Problem 4 is a variant of GP: find E"' E 2, a point x': E R(T), and admissible control u':'(o) such that O. K(X), 0 c <:; d::, 0 E K(O'::), x:-' e Y(wO), and x': = x,.(T). The minimum cost (fuel) is wo,'. A contact function sx(w, y), y E E, is sR(T, y) given by (5.2.20). Suppose a contact function s (y), y E Em, of W(T) is known. Write y = (y0, y) where y0 E E1, yE Em. Then a contact function sy(o, y) of Y(w) is (w, SW(y)) for jYI > 0; sy( w y) = arbitrary x E Y(w) for lY = 0o If W(T) is a single point w(T), then sy(wg y) = (w, w(T)) for all y. A contact function of K(w) is s(w, y) = sx(w y) - s (w, - y). Since qX(w, y) does not vary with w and Tl y(o, y) = y s y(w y), it is clear that the condition on the difference quotient of il(o, y) in Theorem 4. 3. 1 is satisfied.

106 For Problem 5 define sets = {u: luIJ S X, i =1, 2,..,r}i, (5. 6.1) and T C( T) = ({ -(v) B(-) u(o) do: u(.) measurable with range in U } 0 (5. 6. 2) sup ui A control u(') on [0, T] is said to have effort A if 0<t< TIu (t)I = X. Sup-1 pose P (T) w(T) - x(0) E C (T). Then there exists a measurable u(.) X T -T with range in U such that w(T) = 4(T) x(0) + f)(T) - (a) B( ) u(-) do, Xp( T) (O T0 which implies that u(.) transfers x from x(0) to w(T) in time T. Observe CX(T) = KXC1(T) = {. c = Xc, c E C1(T)}. It follows that if X' is the smallest X such that - (T)w(T) - x(0) E C>(T) and w'* is the largest w such that w[4l (T) w(T) - x(0)] E Ci(T), then X* = (=,)-1 A -1 -1 max One choice for A is |I- (T) w(T) - x(0)| cEC(T) Ic. For w E Q2 = [0, A] let X(w) = C1(T) and let Y(w) be the point r4[ -(T)w(T) -x(0)]. Since C (T) is compact and convex (Section 5. 2), K(o) = X(w) - Y(o) has these properties. Clearly K(w) is continuous on i2 and has a contact function s(w, y) = sc (T,y) - o[C l(T) w(T) - x(0)j. where s (T, y) is given by (5. 2.15) and (5. 2.13) with U = U1 and f(u, t) = B(t)u. Thus Problem 5 is a variant of GP: find c.: E e2, o' > 0, and a measurable control u*(-) with range in U1 such that 0 d K(w3), 0:< < X:w-, O E K(':-), and x( )- -1,(T) - w(T). The function (wo*) u'( ) is a minimum effort control and (o*') is the minimum effort. Note that the support function ](w, y) = y - s (T, y) - wy- [, (T) w(T) - x(O)] has a difference quotient which satisfies the requirement in Theorem 4.3.1. Now consider Problem 6. Let X = t, ~ = O, X(X) - R(t), and Y(0) = W(t). Then K(w) - X(M) - Y(W) is compact, convex, and continuous on Q. A contact function of K(X) is s(h, y) - SR(~ y) - SW(, y), y E E, where sR(, ) is given by (5. 2. 8) and Sw(~ ) is a contact function of W(Q.~.

107 The optimization objective in Problem 6 is: find w' E S2, a point x' E R(w"), and an admissible control u*,(.) such that 0 4 K(W), o0 _ w < oae,- 0 E K(w*,)), x=< E W(&w*), and x*c = x w(:) The minimum time is c*. The support function of K(w) is q(ow y) = lR(O, Y) + W(, -Y). W (Wa, Y) - W(Wb, Y) Assume that is bounded from above for all w () - wb a a b Wb e, wa X W b9 bounded y E E. This is certainly true if the target set W(w) is a single point w(wo) which satisfies a Lipschitz condition on i. From (5. 2. 6) nR ( y ) = y' 4(W) X(0) + ymOaxU [(-, y). f(v, c-)]dc. (5.6.3) 0 Thus anlR(aW Y) max RI = y. A(w) (w)x(O) + vEU [Ww(, Y) f(vI W)] (5.6.4) aqR(w Y) Since A(.), ( ), i(-, y), and f(v,.) are continuous on., a exists'rlR(W a, Y) rl R( b Y) and is bounded for all bounded y. It follows that - a(wb~ a?(wa, Y) rl( b Y) a b and are bounded from above for all a WX b E a a - b Wa bi Wb bounded y EE Consequently for Problems 4, 5, and 6 GIP can be applied and convergence is guaranteed. An arbitrary admissible control is chosen to initiate the procedure. Then GIP is continued until a point z E K('tX). w E Q, satisfying specified terminal criteria is found. Finally an admissible control u(.) which generates an approximate optimal state is determined in a manner similar to that described in Section 5. 3.

108 5.7 Numerical Results for GIP Applied to a Minimum-Fuel Example Barr and Gilbert [B3] presented an iterative technique for solving minimum-fuel terminal control problems. The technique is a special case of GIP and uses BIP for the minimization problems which occur in Step 1 of each iteration. Numerical results for a minimum-fuel example using this technique were obtained by Hutcheson [H5]. Computations have been carried out for this same example using II P (with Selection Rule A and p = n) in Step 1 of each iteration. Results from these computations are presented in this section and are compared with Hutches on's data. The example treated is a particular case of Problem 4: A = 0 0 f =, U = {u: u < 1}, a = 0, f0 =!ui, T = 4, W(T) = 0. Because of the simplicity of this example, analytical results can be obtained [e. g., F5]. In particular, the set of all initial states x(0) which can be transferred to the origin in time T with admissible controls can be drawn in the 2-dimensional state space. This set is shown in Figure 5. 7. 1, which also indicates some isofuel contours and the initial states used for the numerical computations. Clearly the minimum fuel -* E [0, 4] and if x(0) 1 0, c': > 0. For an initial state x(0) lying in the shaded region there is no unique optimal control. Since the isofuel contours contain straight-line segments, it follows that the reachable set R(T) c E3 is not strictly convex. Thus the gradient-type convexity methods are not applicable to the example. A detailed description of the iterative procedure for the case when BIP is used in Step 1 of GIP is presented in [B3]. The procedure is much the same when IIP is used. For this simple example it is rot necessary to actually solve (5. 2.17), (5. 2.19), and (502.21) on the computer. As described by Hutcheson, a contact function s R(T, y) can be evaluatled in terms of certain switching times which are related algebraically to the vector y.

109 Fuel = 4' // Fuel 1 / // /,. l /.,I. // -/ / / I! (7, -31 Figure 5.7 1 Initial conditions and iscfuel contours for the minimum-fuel example.

110 Let e be the unit vector (1, 0, 0). Then according to Section 5. 6 K(wo), w E Q = [0, 4], is equal to R(T) - we. For x E R(T) define the scalar functions (x) = Ix - (x e)ej, (5.7. 1) <(X) = x. e. (5.7.2) If x = (x~, x), p(x) is the Euclidean distance of the state x from the origin and 7(x) = x~ is the fuel associated with the state x. Suppose an w E [0, wo"'] and a z E K(w) are generated by GIP. Since z + we E R(T) and ~(z + we) - wA* < ~(z + we) - w, it follows that 5(z + we) - w is an upper bound on the difference between the fuel corresponding to z = (x~, x) -we and the minimum fuel w'; furthermore, p(z + we) measures the error in satisfying the terminal constraint x = 0, Hence a reasonable computational objective is to find z E K(w) such that ( z + we) - w < E1 and T(z + we) < EZ (E1, Ez specified arbitrary positive constants). Consider GIP. Let 2(i) denote the number of points generated by BIP or II P (including the initial point) which occur on iteration i of GIP. The kth point generated by BIP or IIP on iteration i of GIP will be denoted z(j) where i=1 j = k +, (q), i # 0, k = 01,.. {,e(i) - 1 (5. 7. 3) q=o = k, i = 0, k = 0,1... (0) - 1. This notation avoids the use of a double index. In the notation of i Chapter 4 then, z. = z(i F (q) - 1). Note that j + 1 represents the number 1 q=0 of points generated up to and including z(j). Also j is the number of times the contact function has been evaluated. Define p- = -(z(j) + we), (5.7.4) -;z+ce 1 (5. 7. 5) J 1

111 and 6. = c(z(j) + w.e) -.. (5. 7.6) J3 ~ 1 1 The sequences {pj}, {cj}, and {6} are used to display the numerical results. A fixed 0. = 0 was used for each computation and various initial states x(O) (shown in Figure 5.7.1) were investigated. For iteration 0 of GIP a point z(0) in K(0) = R(T) with which to initiate BIP or IIP was found by solving (5. 2. 17) analytically using the admissible control u(t) = -, t c [0, 4]. For iteration i > 0 of GIP the point z. -(W. -a.W )e 1-1 1 1-1 was used to initiate BIP or IIP (see Figure 5.7. 2). In Table 5. 7. 1 the first j for which 6. s E and p < c is shown for J j the case x(0) = (2, -1) and a variety of 0 values. The same quantities are shown in Tables 5. 7. 2, 5. 7. 3, and 5. 7. 4 for the cases x(0) = (2, 0), x(O) = (7, -3), and x(0) = (7. 5, -3). A comparison of the data in A and B of Tables 5. 7. 1 and 5. 7. 2 makes clear the significant improvement achieved when II P is used rather than BIP. Observe that there is little dependence on 8; however, an intermediate value, say 0 =. 4, seems to be most desirable. Figures 5. 7. 3 through 5. 7. 6 show the details of the sequences {6 } and {PJ} for 0 =. 4, BIP and I I P, and initial states (2, 0) and (2, -1). Similar behavior was observed for other 0 values. Note the especially significant improvement which was achieved for the initial state (2, -1), which lies in the shaded region of Figure 5. 7. 1. That convergence is indeed rapid for initial states in this region is verified by the data for x(0) = (7, -3) and x(0) = (7. 5, -3) in Tables 5. 7. 3 and 5. 7. 4. Tables 5. 7. 5 through 5. 7. 10 show the functions xi and #(i) for the various 0, initial states, and BIP and II P. Furthermore, the sequence {4j} is illustrated in Figures 5. 7. 7 and 5.7. 8. This again shows the

112 Starting point for II P or BIP (iteration i - 1): i-2 z. -(Wo. -W )e z(3; I(q)) 1-2 1-1 1-2 q=0 e Terminal point for II P or X i - l )A BIP (iteration i - 1) i-1 z. = z( 2I(q) -1) 1-1 IZ q=o ~~K(wQS, ), CA). _ i). Starting point for II P or BIP (iteration i)o P/. -Z ). 7 2(,Prm~-rT f-r Mnil-fi 1 1-l Figure 5.7.2 Geometry forz minimum-fuel example.

113 improvement possible when IIP is used rather than BIP on each iteration of GIP. It can be seen that a sizeable decrease in 6. occurs on transitions J where i is increased by one (i. e., transitions from k f 0 to k = 0). On these transitions both pj and.j are constant.

114 Table 5. 7. 1 The first j for which 6.j E and, s E; x(O) -(2, -1). A B 1 1.5.05. 0110-3 1. 5.1.05.01 10-3 ~ 1 3 4 1 1 11 11 11 4 8 199 >499.3 3 4 11 11 11 11 4 8 1 53 >499 ~ 4 3 4 11 11 11 i 1 4 8 85>1499.~ 5 3 4 1 1 1 1 1 1 11 4 8 81 >499. 7 3 5 9 9 9 9 4 81> 19. 9 3 6 7 11 11 11 4 28> 19 A: II P Selection Rule A with p = 3 used on each iteration of GIP. B: BIP used on each iteration of GIP. Table 5. 7. 2 The first j for which 6. 5 E and pj E; x(O) =(2, 0). J j A B 1.5.1.05.0110 1.5.1.05.01 10 o 1 4 5 16 19 29 39 4 7 94 96 379 >499 ~3 2 5 12 15 21 27 4 9 76 91 244 >499. 4 2 5 14 14 14 24 3 13 25 40 129 436 o 5 2 5 14 14 14 24 3 13 25 40 349 >499 ~ 7 2 5 13 21 31 45 3 12 70 139 295 >499 ~9 2 7 18 24 34 41 3 17 85 99 >99 A: II P Selection Rule A with p = 3 used on each iteration of GIP. B: BIP used on each iteration of GHP.

Table 5. 7. 3 The first j tao which 5. and.p - ~ J J x(0) (=7, -3); IiF Selection Rule A with p = 3 used on each iteration of GIP. E -3 0 1.5.o.05.01 10 1 2 4 11 11 i ii.3 2 4 1 1 11 11 11. 4 2 4 1 1 iI 11 II.5 2 4 il 11 11 11. 7 2 5 7 7 27 27.9 4 6 14 15 15 15 Table 5. 7. 4 The first j for which 6. < E and j < x(0) - (7. 5, -3); i P Selection Rule A with p = 3 used on each iteration of GiF. I o 5. *.05.0 1 lo. I 2 5 9 9 1! 16.3 2 3 8 8 8 8.4 2 3 8 8 8 8.5 2 3 8 8 8 8.7 2 3 8 8 8 8 ~ 9 4 6 13 16 2! 26

116 Table 5. 7. 5 The functions.i and 2(i) for x(O) = (2, -1); II P Selection Rule A with p = 3 used on each iteration of GIP. 0e 0 1 2.~ 0. 527. 000 ~3 0 ~ 527 1. 000.4 0. 527 1. 000. 5 0. 527 1. 000 ~ 7 0 659 1. 000 ~ 9 0. 8 52 1. 000.1 4 2 5.3 4 2 5.4 4 2 5 2(i).5 4 2 5.7 5 2 2.9 6 3 2

117 Table 5.7. 6 The functions 0i and 2(i) for x(0O) = (2, Ki); BIP used on each iteration of GIP. o 0 1 2 3 4 5 6 7 8 ~ i 0. 599. 725.793 857.908.935.954.968.3 0. 599. 806.928 o 4 0 o 599. 823.945.5 0.776. 946.7 0.776.9 0.845.1 3 9 4 4 2 48 129 260 > 40.3 3 14 2 >480.4 4 16 62 >1417 2(i).5 7 68 >424.7 7 >193.9 27 >172

1 8 Table 5.7.7 The functions o. and f(i) for x(O) = (2, 0), 1 II P Selection Rule A with p = 3 used on each iteration of GIP. 0 1 2 3 4 5 ~ 1 0 [ 545 1. 032 1. 137 1. 161 1. 171. 3 0 1.000 1. 101 1. 168 1, 171 4 0 1. 000 1. 169 1. 171. 5 0 1.000 1. 169 1. 171.7 0 1. 000 1. 152 1. 171 ~9 0 1.071 1. 169 1. 171.1 2 5 5 10 4 13 J4(i)...5 4 9 10 1 ~ 7 4 10 18 13.9 5 21 14 1

119 Table 5. 7. 8 The functions w. and I(i) for x(0) = (2, 0); BIP used on each iteration of GIP. o 0 1 2 3 4 5. 1 0. 545 1. 113 1. 150 1. 166 1. 171 ~ 3 0. 560 1. 138 1. 171 ~ 4 0. 897 1. 126 1. 167 1. 171. 5 0. 897 1. 126 1. 171 ~ 7 0. 907 1. 117 1. 167 ~ 9 0.961 1. 166 1 1 4 7 4 26 45 >349 e3 3 2 96 >398 ~ 4 7 1 5 22 127 265 2(i).5 7 15 80 >397.7 9 57 99 >334.9 15 67 >17

120 Table 5. 7. 9 The functions wi and 1(i) for x(0) = (7, -3); II P Selection Rule A with p = 3 used on each iteration of GIP. i o 0 1 2 3. 1 0 1. 7 50 2. 647 3. 000.3 0 1. 750 2. 647 3. 000. 4 0 1.7 50 2. 647 3. 000 1 5 0 1. 7 50 2. 647 3. 000.7 0 2.002 3.000 ~ 9 0 2. 543 3. 000.1 1 2 7 1.3 1 2 7 1.4 1 2 7 1 P(i).5 1 2 7 1.7 2 4 21.9 3 7 5

121 Table 5. 7. 10 The functions o. and (i) for x(O) = (7. 5. -3); II I P Selection Rule A with p = 3 used on each iteration of GIP. 0 0 1 2 3 X1 0 2.125 2. 545 3.000 * 3 0 2. 125 3. 000. 4 0 2. 125 3. 000. 5 0 2. 125 3. 000.7 0 2. 125 3. 000 ~ 9 0 3. 000.1 1 2 3 10.3 1 3 4.4 1 3 4 (i).5 1 3 4.9 1 3 4.9 3 23

122 Pi 1 01_ II 6. _ 13 0 1 -50 2 50 3 50 I III A 1 1 ~ Note: k 0 0 for j = 0, 4, 20, 82; the arrows denote transitions from k 0 0 to k =O. 10-3 I I I I, 0 10 20 30 50 150 250 350 j Figure 5.7.3 6j and p for x(0) = (2, -1), 0 =.4; BIP used on each iteration of GIP.

123. 6. (6j negative for i j =6, 7) V I.1-.01- - Note: k = 0 for j = 0, 4, 6; the arrows denote transitions from k $ 0 to k 0. 10 I I I i I I I 0 10 20 30 40 50 60 70 80 90 100 Figure 5.7.4 6j and p. for x(0) = (2, -1), 0 =.4; IIP Selection Rule A with p = 3 used on each iteration of GIP.

124 6. (6j negative for j = 22, 23, 24, 44) Note: k = 0 for j = 0, 7, 22, 44, 171; the arrows denote transitions from k 0 O to k = 0. \An ~n! 1 -31 i 10-3 LI 1l 1 1 {l~ll 1 - 0 10 20 30 50 150 250 350 j Figure 5.7.5 6j and pj for x(O) = (2, O), E =.4; BIP used on each iteration of GIP.

125 |I- i. (6. negative for j = 4, 13) Note: k = 0 for j = 0, 4, 13, 23; the arrows denote transitions from k # 0 to k = 0. 11 3 v-I 10-3 0 10 20 30 40 50 60 70 80 90 100 j Figure 5.7.6 6j and pj for x(0) = (2, 0), 0 =.4; IIP Selection Rule A with p = 3 used on each iteration of GIP.

126 2.0 ------ A (A = B for j $ 2) B 1. 8 1. 6 | Note: With A, k = 0 for j = 0, 4, 6; upward arrows denote transitions from k $ 0 to k = 0. With B, k = 0 for j = 0,4, 20, 82; downward arrows denote transitions from k ~ 0 to k = 0. 1.4 1,2 I i':- - 1. 000 0.8tt I I I, I I I 0 10 20 30 50 1 50 250 3 50 j Figure 5.7,7. for x(O) =(2, -1), 0 =.4: A) IIP Selection Rule A with p = 3 used on each iteration of GIP; B) BIP used on each iteration of GIP.

127 2.,0 --- A (A = B for j < 2) - B 1. 8 Note: With A, k = 0 for j = 0, 4, 13, 23; upward arrows 1. 6 denote transitions from k i 0 to k = 0. With B, k = 0 for j = 0, 7, 22, 44, 171; downward arrows denote transitions from k f 0 to k = 0. 1.4 1. 1 0 ___ I I I I I _ "I 0 10 20 30 50 150 250 350 j Figure 5.7.8 j for x(0) -(2, 0), 0 =.4: A) IIP Selection Rule A with p = 3 used on each iteration of GIP; B) BIP used on each iteration of GIP.

List of Symbols A(t) continuous m X m matrix function. B(t) continuous m X r matrix function. t C(t) f{f -l() f(u(r), r) dcr: u(') admissible}. 0 D m X m symmetric non-negative definite matrix. D (m + 1) X (2m + 2) matrix defined by (3.8.1). En n-dimensional Euclidean space. G m X m symmetric non-negative definite matrix. H convex hull of m points; convex polyhedron. Hk convex hull of yl(k), y (k), s(-Zk)' Zk. I identity matrix. ik the set of i, 1 i is p + 2, for which xi = O. I(z) [min {6 5 z), 1}, min {(2 - 6) (z), 1 ]. t JY(u) f [a(fl) - Xu(a) + fo (u(a), a-)] du-+ h~(x (t)). K compact, convex set in En K(w) compact, convex sets in En continuous on 2. L(x;y) {z z = x +co(y-x), -o0< co < co}, x / y; line. Mi closed intervals in E' N(x;w) {z Jz-xj < w}, w > O; open sphere. N(x;w) {z Iz-xf < W}; closed sphere. O origin. P(y) the support hyperplane {x ~ x- y = r(y)}, y / 0, of K with outward normal y. 128

129 P(w, y) the support hyperplane {x: x y = (w y), y ) 0, of K(o) with outward normal y. P H(Y) the support hyperplane of H with outward normal y # 0. Q(x;y) {z ~ z. y = x y}, y # 0; the hyperplane through x with normal y. Q+(x;y) {z: z* y < x. y}, y # 0; the open half-space bounded by Q(x;y) with outward normal y. Q-(x;y) {z: z y > x y}, y # 0; the closed half-space bounded by Q(x;y) with inward normal y. R(t) {x: x = x (t), u(') admissible}; reachable set. Sk y), (k)k), s(-zk)}. S(y) P(y)nK; the contact set of K for outward normal y. S(o, y) P(w, y)n K; the contact set of K(c) for outward normal y. T fixed terminal time. U compact set in E. W(t) compact, convex, and continuous sets in Em X(t) compact, convex, and continuous sets in En Y(Wo) compact, convex, and continuous sets in En Y~k {y1 (k),..., yp (k)}. a(t) continuous function from, T] to E a(t) continuous function from [0, T] to Em b(t) continuous function from O to Em. cU(t) c (t) = P (t) x u(t) - x(O). e the 3-vector (1, 0, 0). f(u, t) continuous function from U X O to EM f~(u, t) continuous function from U X Oto E1'.

130 g m-vector in the range of G. h~(x) convex function from EM to E1 i iteration number for GIP. i-1 ~j k + Z f(q), i ~ 0, k, i = 0. q =o k iteration number for BIP or II P. e(i) number of points generated by BIP or IIP (including the initial point) which occur on iteration i of GIP. p number of points in Yk; parameter of II P. n max s(y) function from E to K such that y - s(y) = zEK Y z, y 0; contact function of K. s(w, y) function from E to K(w) such that y ~ s(w0, y) = zEK(aw) Y z y / 0; contact function of K(o). sX (W, y) contact function of X(X). t time. u(') r-dimensional control vector function defined on O. w(T) target set consisting of a single point. x(t) m-dimensional state vector. xI', solution of Subproblem 2. x (t) an absolutely continuous solution of (5.1.1) for the admissible control u('). y * solution of Subproblem 1. rain z(00) point in K(X) satisfying I'(tw) = zCK(w) Jz. z. terminal point generated by BIP or II P on iteration i of GIP. z(j) the kth point generated by BIP or II P on iteration i of GIP. z,- =solution of BP.

131 a k an element of I(zk). (z) Iz-s(-z) I z (z-s(-Z)), Z-S(-Z) /0; 0, Z-S(-Z) = 0. Y(z) zL z- s(-z), Izl > 0 and z. s(-z) > 0; 0, z = 0 or Izl > 0, z s(-Z) < 0, -2 y(w, z) Izl z s(, -z), Iz > 0 andz. s(w, -z) > 0; 0, z=0or z > 0, z. s(W, -z) -< 0. 6 fixed number in (0, 1] used for I(z); parameter for BIP.,5j o(z(j) + Wie) - wi. max ](Y) zEK Z. y; support function of K. max l(w, Y) zEK(co) z. y; support function of K(w). 0 parameter for GIP. Xk principal radii of curvature at z*' for Example 3 of Section 2.5. max X X.'R(z) Izl z. s(-z), z / 0; 0, z - 0. Pji PFz(j) + Woie). ~(x) x' e. -(z(j) + Wie). ~i(u, t) continuous functions from U X O to E1. P( t) m-dimensional adjoint vector. p(t, y) solution of (5, 2. 5) with boundary condition p( T) = y. W scalar E 2. o* ~solution of GP.

132 P(z) Izlz - Jz-, z E K. a denotes "the boundary of" denotes "the convex hull of' A( z;o) r(z) - r(z + a(s(-z) - z)). e compact interval in E'. ~ (t) solution of D = A(t)b,'(0) = I. compact interval in E1.

Bibliography Al Arrow, K.J., Hurwicz, L., and Uzawa, H., "Studies in Linear and Nonlinear Programming, "' Stanford University Press, Stanford, Calif., (1958). B1 Babunashvili, T. G., Synthesis of Linear Optimal Systems, (Russian) Doklady Akademii Nauk SSSR, Vol. 155, (1964), pp 295-298, (English translation in SIAM J. on Control, Ser. A, Vol. 2, (1964), pp 261-265). B2 Barankin, E. W. and Dorfman, R., "Towards Quadratic Programming, " Office of Naval Research Logistics Project at Columbia University and University of California, Berkeley, (19 55). B3 Barr, R. 0. and Gilbert, E. G., Some Iterative Procedures for Computing Optimal Controls, to appear in "Proceedings of the Third Congress of the International Federation of Automatic Control (IFAC), London, 1966, " Butterworth, London. B4 Beale, E. M. L., On Minimizing a Convex Function Subject to Linear Inequalities, J. of Royal Statistical Society, Ser. B, Vol. 17, (1955), pp 173-184. B5 Beale, E. M. L., On Quadratic Programming, Naval Research Logistics Quarterly, Vol. 6, (1959), pp 227-243. B6 Busemann, H., "Convex Surfaces, " Interscience Tracts in Pure and Appl. Math. No. 6, Interscience Publishers, Inc., New York, (19 58). B7 Bushaw, D., Optimal Discontinuous Forcing Terms, in "Contributions to the Theory of Nonlinear Oscillations, " Princeton University Press, Vol. 4, (1958). B8 Bellman, R., "Dynamic Programming, " Princeton University Press, Princeton, N. J., (19 57) B9 Breakwell, J. V., Optimization of Trajectories, J. SIAM, Vol. 7, No. 2, (1959), pp 215-247. B10 Breakwell, J. V., Speyer, J. L., and Bryson, A. E., Optimization and Control of Nonlinear Systems Using the Second Variation, SIAM J on Control, Ser. A, Vol. 1, No 2, (1963), pp 193 -223. 133

134 B ll Bryson, A. E., and Denham, W. FK, A Steepest-Ascent Method for Solving Optimum Programming Problems, J, Appl. Mech,, ASME Trans., Ser. E, Vol. 29, (1962), pp 247-257, C 1 Carathe'odory, C., Uber den Variabilitaitsbereich der Fourierschen Konstanten von positiven harmonischen Funktionen, Rend. Circ. Mat. Palermo, Vol. 32, (1911), pp 1P3-217, D1 Dantzig, G. B., Maximization of a Linear Function of Variables Subject to Linear Inequalities. in "Activity Analysis of Production and Allocation, " (T. C. Koopmans, editor), Wiley, New York, (19 51). El Eaton, J. H., An Iterative Solution to Time-Optimal Control, J. Math. Anal. and Appl., Vol, 5, (1962), pp 329-344, Errata and Addenda, J. Math. Anal. and Applo, Vol. 9, (1964), pp 147152. E2 Eggleston, H. G., "Convexity, "Cambridge Tracts in Math, and Math. Physics No. 47, Cambridge University Press, (1958). Fl Fadden, E. J,, "Computational Aspects of a Class of Optimal Control Problems,' Ph. D. thesis, University of Michigan, 1965. F2 Fadden, E. J. and Gilbert, E. G., Computational Aspects of the Time-Optimal Control Problem, in "Computing Methods in Optimization Problems, Academic Press, New York, (1964), pp 167192. F3 Fancher, P. S., Iterative Computation Procedures for an Optimum Control Problem, IEEE Trans. on Auto. Control, Vol, AC-10, No. 3, (1965), pp 346-348. F4 Frank, M. and Wolfe P., An Algorithm for Quadratic Programming, Naval Research Logistics Quarterly, Vol. 3, (1956), pp 9 5-110. F5 Flugge-Lotz, I., and Marbach, HF, The Optimal Control of Some Attitude Control Systems for Different Performance Criteria, J. Basic Engr., ASME Trans., Ser, D, Vol. 85, (1963), pp 16517 6 G1 Gilbert, E, G., An Iterative Procedure for Computing the Minimum of a Quadratic Form on a Convex Set, SIAM J. on Control, Ser. A, Vol. 4, No. 1, (1966), pp 61-800

135 G2 Gilbert, E. G., The Application of Hybrid Computers to the Iterative Solution of Optimal Control. Problems, in "Computing Methods in Optimization Problems " Academic Press, New York, (1964), pp 261-284. H1 Hadley, G., "Nonlinear and Dynamic Programming, " AddisonWesley, Reading, Mass., (1964). H2 Hildreth, C., A Quadratic Programming Procedure, Naval Research Logistics Quarterly, Volo 4, (1957), pp 79-85. H3 Ho, Y. C., A Successive Approximation Technique for Optimal Control Systems Subject to Input Saturation, J. Basic Engo, ASME Trans., Ser, D, Volo 84, No. 1, (1962). pp 33-40. H4 Houthakker, H. S., The Capacity Method of Quadratic Programming, Econometrica, Vol. 28, (1960), pp 62-87, H5 Hutcheson, R. D., The Computation of Optimal Controls by an Iterative Procedure, Professional Degree Thesis, University of Michigan, 1965. K1 Kelley, H. J., Gradient Theory of Optimal Flight Paths, J, American Rocket Society, Vol. 30, (1960), pp 947-953. K2 Knudsen, H. K,, An Iterative Procedure for Computing Optimal Controls, IEEE Trans. on Auto, Control, Vol. AC-9, No. 1, (1964), pp 23-30. K3 Kuhn, H. W., and Tucker, A. W,, Nonlinear Programming, Proceedings Second Berkeley Symposium on Mathematical Statistics and Probability, (1951), pp 481-492. L1 La Salle, J. P., The Time-Optimal Control Problem, in "Contributions to the Theory of Nonlinear Oscillations, " Volo 5, Princeton University Press, Princeton, N. J., (1960), pp 1-24. L2 Lemke, C. E., A Method for Solution of Quadratic Programs, Management Science, Vol. 8, (1962), pp 442-453. M1 Markowitz, H., The Optimization of a Quadratic Function Subject to Linear Constraints, Naval Research Logistics Quarterly, Vol. 3, (1956), pp 111-133. N1 Neustadt, Lo W. and Paiewonsky, Bo, On Synthesizing Optimal Controls, in "Proceedings of the Second Congress of the Inter

136 national Federation of Automatic Control (IFAC), Basel, 1963", Butterworth, London. N2 Neustadt, L. W., The Existence of Optimal Controls in the Absence of Convexity Conditions, J. Math. Anal. and Appl., Vol. 7, (1963), pp 110-117. N3 Neustadt, L. W., Minimum Effort Control Systems. SIAM J. on Control, Ser. A, Vol. 1, (1962), pp 16-31. N4 Neustadt, L. W., Synthesizing Time Optimal Control Systems, J. Math. Anal. and Appl., Vol. 1, (1960), pp 484-492. P1 Plant, J. B., and Athans, M., An Iterative Technique for the Computation of Time Optimal Controls, to appear in "Proceedings of the Third Congress of the International Federation of Automatic Control (IFAC), London, 1966, " Butterworth, London. P2 Pontryagin, L. S., Boltyanskii, V. G., Gamkrelidze, R. V., and Mishchenko, E. F., "The Mathematical Theory of Optimal Processes, " John Wiley and Sons, Inc., New York, (1962). Vl Vajda, S., "Mathematical Programming, " Addison-Wesley, Reading, Mass., (1961). V2 Valentine, F. A., "Convex Sets, " McGraw-Hill, New York, (1964). W1 Wolfe, P., The Simplex Method for Quadratic Programming, Econometrica, Vol. 27, (1959), pp 382-398.

UNIVERSITY OF MICHIGAN 1139015 02229 0129111111 3 9015 02229 0129