SIGNOMIAL PROGRAMS WITH EQUALITY CONSTRAINTS: NUMERICAL SOLUTION AND APPLICATIONS b y J A M E S Y A N B . A . S c . , U n i v e r s i t y o f B r i t i s h C o l u m b i a , V a n c o u v e r , 1 9 6 9 - M . A . S c . , U n i v e r s i t y o f B r i t i s h C o l u m b i a , V a n c o u v e r , 1 9 7 1 •A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n the Department of E l e c t r i c a l Engineering We accept t h i s t h e s i s as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1 9 7 6 © James Yan, 1 9 7 6 In presenting th i s thes is in pa r t i a l fu l f i lment of the requirements for an advanced degree at the Univers i ty of B r i t i s h Columbia, I agree that the L ibrary shal l make it f ree ly ava i l ab le for reference and study. I further agree that permission for extensive copying of th i s thes is for scho lar ly purposes may be granted by the Head of my Department or by his representat ives. It is understood that copying or pub l i ca t ion of th is thes is for f inanc ia l gain sha l l not be allowed without my written permission. Depa rtment The Univers i ty of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date ABSTRACT Many design problems from diverse engineering d i s c i p l i n e s can be formulated as signomial programs with both e q u a l i t y and inequa-l i t y c o n s t r a i n t s . However, e x i s t i n g computational methods of signomial programming are applicable to programs with i n e q u a l i t y constraints only. In t h i s thesis an algorithm i s proposed and implemented to solve signo-m i a l programs with mixed i n e q u a l i t y and e q u a l i t y c o n s t r a i n t s . The algo-rithm requires no manipulation of the constraints by the user and i s compatible with the r e s u l t s of signomial programming. The proposed algorithm i s a synthesis shaped by three concepts: the method of m u l t i p l i e r s viewed as a primal-dual method, p a r t i a l dua-l i z a t i o n , . a n d the r e t e n t i o n of the structure of a signomial program. The strategy of the algorithm i s to replace the o r i g i n a l problem with a se-quence of subproblems each subject to the o r i g i n a l signomial i n e q u a l i t y constraints only. The subproblem's objective function i s the o r i g i n a l cost augmented with the equality constraints and s p e c i f i e d by a parameter vector X_ and a penalty constant K. The algorithm then alternates between s o l v i n g the subproblem and updating X_ and K. The convergence of the algorithm i s , under s u i t a b l e assumptions, guaranteed by the convergence of the method of m u l t i p l i e r s and the method used to solve the subproblem. Because the subproblem has the form of a regular signomial program, i t can i n p r i n c i p l e be solved by the techniques of signomial programming. A new numerical method implementing a v a r i a n t of the A v r i e l -Williams algorithm for signomial programs i s suggested for s o l v i n g the subproblem. The method r e l i e s on monomial condensation, and the combined use of the reduced gradient method and a c u t t i n g plane scheme. D e t a i l s i i of the method's implementation are considered, and some computational experience has been acquired. The proposed method also has the advan-tageous f l e x i b i l i t y of being able to handle non-signomial d i f f e r e n t i a b l e objective functions. Four updating schemes for X_ and K are formulated and evaluated i n a s e r i e s of numerical experiments. In terms of the rate of convergence, the most promising scheme tested i s the use of the Hestenes-Powell r u l e f o r updating A_ and the moderate monotonic increase of K a f t e r the comple-t i o n of each subproblem. Convergence can also be considerably accele-rated by properly s c a l i n g the e q u a l i t y constraints and performing only inexact minimization i n the f i r s t few subproblems. The a p p l i c a b i l i t y of the algorithms developed i n t h i s thesis i s i l l u s t r a t e d with the s o l u t i o n of three design examples drawn from s t r u c t u r a l design and chemical process engineering. i i i TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS . . . . . . . . . . . . i v LIST OF ILLUSTRATIONS v l LIST OF TABLES . . v i i ACKNOWLEDGEMENT , v i i i . I. INTRODUCTION 1 1.1 Nonlinear Programming and Engineering Design 1 1.2 A Preview of the Thesis 2 I I . SIGNOMIAL PROGRAMMING 6 2.1 Introduction" . . . 6 2.2 Posynomlals and Their Condensation . . . ^ 2.3 Signomial Programming 9 2.4 Geometric Programming H 2.5 Solution of Geometric Programs 2.6 Solution of Signomial Programs 15 I I I . SIGNOMIAL PROGRAMS WITH EQUALITY CONSTRAINTS . . . . . . . . 1 9 3.1 Motivation i 9 3.2 Previous Work 2 0 3.3 Problem Statement 2 1 3.4 Proposed Algorithm 2 2 3.4.1 Background . . ; 2 2 3.4.2 Development of the Algorithm 2 3 3.4.3 Statement.of the Algorithm . 27 3.5 Convergence Considerations 2 9 3.6 Solution of the Subproblem 31 3.7 Conclusion 37 IV. A PROPOSED ALGORITHM FOR INEQUALITY-CONSTRAINED SIGNOMIAL PROGRAMS 39 4.1 Introduction . . ^9 4.2 Problem Statement 4 ° 4.3 Proposed Algorithm 4 1 4.3.1 Preview 4 1 4.3.2 The Reduced Gradient Method . 41 4.3.3 Algorithm Development 43 4.3.4 The Algorithm 46 4.3.5 The Algorithm i n Perspective 46 4.4 Considerations f o r Implementation 4 8 4.4.1 Finding a Feasible Point of a Signomial Program. 48 4.4.2 Getting an I n i t i a l Basic Feasible Solution of Problem (LST) - No Cut Added . . . . . . . . . . 49 4.4.3 Getting an I n i t i a l Basic Feasible Solution of Problem (LST) - A f t e r the Addition of a Cut . . 5 1 i v 9 Page 4.4.4 A Linear Search Method . . . . . . 53 4.4.5 Optimality C r i t e r i a , F e a s i b i l i t y Tolerances, and Updating the Basis' Inverse . . . . . . . . 56 4.5 Computational Experience . . . . . . . . . . . . . . . 58 4.5.1 Posynomial Representation . . . . 58 4.5.2 Software Implementation 59 4.5.3 Computational Results 6 0 4.6 Presence of Simple Upper Bounds 66 4.7 Minimizing Algebraic Functionals of Signomials . . . . 70 4.7.1 Introduction . . . . . . . 7 0 4.7.2 Problem Formulation and Solution 7 i 4.7.3 A p p l i c a t i o n to Constrained L o c a t i o n - A l l o c a t i o n Problems 72 V. NUMERICAL SOLUTION OF SIGNOMIAL PROGRAMS WITH EQUALITY CONSTRAINTS . 7 7 7 7 5.1 Introduction 5.2 Solution of the Subproblem . 7 8 5.3 Schemes f o r Updating X_ and K 7 9 5.3.1 Introduction 7 9 5.3.2 Choice of Step Size a 8 1 5.3.3 Updating the Penalty Constant K . . . . . . . . 8 4 5.4 Numerical Experiments ^ 5.4.1 Introduction 8 5 5.4.2 Software D e t a i l s . . . . . . . . . . . 8 5 5.4.3 Test Problems . . . . . . . . . . 8 7 5.4.4 Experiment Plan . 9 0 no 5.4.5 Results and Discussion ^ VI. APPLICATIONS IN ENGINEERING DESIGN . . . . 1 0 0 6.1 Optimum Design of S t a t i c a l l y Indeterminate Pin-Jointed Structures . . . 1 0 0 6.1.1 Introduction . . . . . 1 0 0 6.1.2 Problem Formulation 1 0 0 6.1.3 Example 1 0 3 6.2 Optimization Examples from Chemical Process Engineering 107 6.2.1 Introduction 1 0 7 6.2.2 A l k y l a t i o n Process Optimization 1 0 8 6.2.3 Design of a Heat Exchanger Network H 6 VII. CONCLUSION . 124 7.1 Synopsis' . . . 124 7.2 Contributions of the Thesis 128 7.3 Suggestions f o r Further Research 129 REFERENCES 131 v LIST OF ILLUSTRATIONS Figure Page 5.1 T y p i c a l p l o t of l o g 6(N) against the dual i t e r a t i o n count N , 97 6.1 The h y p e r s t a t i c p i n - j o i n t e d structure of the example, (a) Frame and loading (b) Basic frame with external load (c) Basic frame subject to redundant force . . . . 104 6.2 S i m p l i f i e d a l k y l a t i o n process diagram . 109 6.3 A small heat exchanger network 117 v i LIST OF TABLES Table Page 4.1 Coefficients for Problem C . . . . . . . . . 63 4.2 Statistics of the test problems for the CRGCP algorithm. 63 '4.3 Starting point and solution of (a) Problem A, (b) Pro-blem B, (c) Problem C . . 64 4.4 Coefficients of example problem with simple upper bounds. 69 4.5 Starting points and solutions of the example problem with simple upper bounds . 69 4.6 Locations and weights for major Ukraine cities . . . . . 76 4.7 Summary of the solution of the constrained location pro-blem 76 5.1 Summary of the characteristics of the test problems . . 87 5.2 Performance results of Schemes I-IV 93 6.1 The alkylation process problem's variables and their bounds 110 6.2 Definitions and values of the cost coefficients of the alkylation process problem . . . . . 110 6.3 The i n i t i a l and optimal solution of the alkylation process problem i X 5 6.4 Data for the heat exchanger network design problem . . . 121 6.5 Starting point and optimal solution of the heat ex-changer network problem 123 v i i ACKNOWLEDGEMENT I would l i k e to thank Dr. E r i k V. Bohn f o r introducing me to geometric programming and for supervising my doctoral research. My f i r s t encounter with the method of m u l t i p l i e r s was due to Dr. Mordecai A v r i e l of Technion, I s r a e l . I am also g r a t e f u l to Dr. A v r i e l f o r h i s encouragement and f o r the many valuable discussions we had during h i s stay at U.B.C. i n 1974-75. I wish to acknowledge the use of Dr. Ron Dembo's code GGP made av a i l a b l e to me by Dr. A v r i e l . The timely assistance of my f r i e n d , Rose Cam, during the pre-paration of the f i r s t d r a f t provided a much needed boost. She deserves a bouquet of roses. Many thanks go to Misses Mary-Ellen Flanagan and Sannifer Louie for t h e i r e x c e l l e n t work i n typing the f i n a l copy. F i n a l l y , the graduate fellowships from UBC and a National Research Council of Canada Postgraduate Scholarship are g r a t e f u l l y acknowledged. P a r t i a l support of the project also came from NRC Grant No. A3134. v i i i 1. I. INTRODUCTION 1 . 1 Nonlinear Programming and Engineering Design Design i s the e s s e n t i a l task of engineering. The design pro-cess usually begins with the recognition of a need. The process then proceeds with the analysis of the need, the synthesis of s o l u t i o n con-cepts l i k e l y to s a t i s f y the need, the evaluation of the proposed s o l u -t i o n s , and the s e l e c t i o n of a. f i n a l proposal. The selected concept i s then fur t h e r r e f i n e d and made s u f f i c i e n t l y s p e c i f i c f o r prototype eva-l u a t i o n and p o s s i b l e r e v i s i o n . F i n a l l y , the design plan i s implemented. In most stages of the design process opportunities f o r optimi-zation e x i s t . However, the disparate problems faced i n each stage d i c -t a te the use of r a d i c a l l y d i f f e r e n t optimization methods. One important step of the design process i n which the e f f i c a c y of optimization has been amply demonstrated i s the s p e c i f i c a t i o n of a selected s o l u t i o n con-cept. At t h i s point of the design process, the thorough a n a l y s i s of the perceived need and the evaluation of the set of suggested s o l u t i o n s have usually produced a mathematical model. The model characterizes the structure of the accepted concept, i d e n t i f i e s the d e c i s i o n v a r i a b l e s and t h e i r f u n c t i o n a l r e l a t i o n s h i p s , and l i s t s the constraints that have to be imposed because of p h y s i c a l , economic, or s o c i a l considerations. The task the designer now faces i s to choose the values of the design v a r i a b l e s i n order to impart to the selected s o l u t i o n concept a quanti-t a t i v e i n d i v i d u a l i t y . Most often, there i s a gamut of acceptable choices, with each choice implying a c e r t a i n l e v e l of performance by the model. I f the designer's expectations from the model can be stated mathematically, as they often can, then each chosen set of s a t i s f a c t o r y design v a r i a b l e s can be judged r e l a t i v e to the mathematical norm of performance. This i s z. p r e c i s e l y a s i t u a t i o n r i p e f o r systematic optimization. In many problems of engineering design, the mathematical model i s described by a set of continuous functions of a Euclidean space. For such models, the search f o r the optimal set of design parameters i s nothing e l s e but a nonlinear program. For the realm of nonlinear programming i s to solve the optimization problem minimize (min) f(x) (1.1) subject to (s.t.) g^(x) <_ 0, j = 1, 2, ..., p (1.2) h ^ x ) =0, k = 1, 2, q (1.3) where x i s an m-dimensional vector, and f, g^ ., and h^ are a l l continuous functions of the Euclidean space R m. The designer therefore has at h i s d i s p o s a l the powerful methods of nonlinear programming to seek the o p t i -mum design. He can abandon the t r a d i t i o n a l approach of obtaining only f e a s i b l e s o l u t i o n s . With the a i d of nonlinear programming he can achieve b e t t e r performance with economical savings, improve modeling by s e n s i t i -v i t y a n a l y s i s , and even generate new s o l u t i o n concepts. 1.2 A Preview of the Thesis A signomial i s a nonlinear function defined by the d i f f e r e n c e of two p o s i t i v e sums of power functions. I f a nonlinear program i n -volves only signomials, i t i s known as a signomial program. In engin-eering design, many models from diverse engineering d i s c i p l i n e s can be described by signomial programs. This observation i s confirmed by the i n c r e a s i n g number of published a r t i c l e s applying signomial programming to design. However, while engineering design models generally have both i n e q u a l i t y and e q u a l i t y c o n s t r a i n t s , the theory of signomial programming has been developed i n terms of i n e q u a l i t y constraints only. Applying the 3. theory to solve design problems, therefore, often requires the transforma-t i o n of e q u a l i t y constraints to i n e q u a l i t i e s . While such a transformation can be conveniently c a r r i e d out i n some cases, i t i s not c l e a r how the transformation should be undertaken i n general. The aim of t h i s t h e s i s i s to develop a numerical method which the designer can use to solve signomial programs with mixed constraints without having to change one type to the other. The method i s s t i l l rooted i n the theory of signomial programming. Consequently, the method can make f u l l use of the r e s u l t s of the theory. In order to provide the background necessary f o r l a t e r discuss-ions, the main r e s u l t s of signomial programming are summarized i n Chapter II. Geometric programming i s reviewed as a s i g n i f i c a n t s p e c i a l case of signomial programming. Methods f or s o l v i n g both geometric programs and signomial programs are discussed. In Chapter III, the ce n t r a l problem of signomial programs with i n e q u a l i t y and e q u a l i t y constraints i s formulated and solved by a pri m a l -dual method. The development of the proposed algorithm i s based on three ideas: the method of m u l t i p l i e r s viewed as a primal-dual algorithm, p a r t i a l d u a l i z a t i o n , and the retention of a signomial program's structure i n the primal problem. The algorithm replaces the o r i g i n a l problem with a sequence of ine q u a l i t y - c o n s t r a i n e d signomial programs dependent on a set of parameters defined as the vector \_ and the s c a l a r K. Between the s o l u t i o n of successive signomial programs, the parameters are updated i t e r a t i v e l y according to some r u l e s . The convergence of the algorithm i s guaranteed by that of the method of m u l t i p l i e r s and by the convergence of the algorithm used to solve the sequence of primal signomial programs. The implementation of the algorithm proposed i n Chapter III requires the s e l e c t i o n of an algorithm f o r s o l v i n g the primal signomial programs and the s p e c i f i c a t i o n of how the parameters A _ and K should be updated. In Chapter IV, a new numerical method based on the A v r i e l -Williams algorithm for s o l v i n g signomial programs i s proposed. In t h i s method the o r i g i n a l nonconvex f e a s i b l e set defined by signomial i n e q u a l i -t i e s i s approximated by a convex set obtained by monomial condensation. The nonconvex objective function i s then minimized over the convex appro-ximant by a combined reduced gradient and cutti n g plane algorithm. The minimization's s o l u t i o n i n turn determines the next convex approximant. Considerations f o r implementing the new method are discussed i n d e t a i l . Numerical experience with the method shows that i t i s comparable to, at times b e t t e r than, another recent implementation of the A v r i e l - W i l l i a m s algorithm. The proposed method, however, has the f l e x i b i l i t y of handling a wider class of objective functions such as algebraic f u n c t i o n a l s of signomials. F i n a l l y , the proposed algorithm i s r e f i n e d to t r e a t separately simple upper-bound constraints. The work i n Chapter V focuses on the performance of a set of updating rules f or the parameters A _ and K. A t o t a l of four combinations of updating rules f o r A . and K i s tested i n a s e r i e s of numerical e x p e r i -ments. In the experiments, a s e n s i t i v i t y study i s also performed to help understand the impact of the range of values of K. The r e s u l t s of the experiments serve as the bas i s f o r the choice of the s u i t a b l e updating rules f o r s o l v i n g the design problems of Chapter VI. To i l l u s t r a t e how the algorithms discussed i n the preceding chapters can be applied to engineering design, three selected design problems are presented i n Chapter VI. The f i r s t i s on obtaining the minimum-cost design of a h y p e r s t a t i c p i n - j o i n t e d s t r u c t u r e . The second i s concerned with achieving the optimum operating conditions of an a l k y l a t i o n process. The l a s t example involves the optimal trade-off between i n s t a l l a t i o n and operating costs of a small heat exchanger network. In each example, the problem i s formulated with s u f f i c i e n t d e t a i l to allow appreciation of the p h y s i c a l basis of the objective function and the con-s t r a i n t s . Numerical data are then s u b s t i t u t e d i n t o the problem to cast i t i n t o a format compatible with t h i s t h e s i s ' algorithms. The numerical s o l u t i o n i s then obtained and interpreted. Chapter VII concludes the thesis by g i v i n g a synopsis of the t h e s i s and some suggestions f o r further research. I I . SIGNOMIAL PROGRAMMING 2.1 Introduction The basic theory of signomial programming i s reviewed i n t h i s chapter. The fundamental notion of a posynomial i s f i r s t defined,, and i t s r e l a t i o n s to convexity and the geometric i n e q u a l i t y are summarized. The formulation of a signomial program i s then given, and the import-ant properties of the formulation are examined. The stated problem i s also compared to other equivalent formulations. I f a signomial program has no negative terms, then a geometric program r e s u l t s . The theory of geometric programming, i n p a r t i c u l a r , i t s d u a l i t y p r o p e r t i e s , are d i s -cussed along the l i n e s of Duffin, Peterson and Zener [ 1 ] . F i n a l l y , algorithms f o r so l v i n g geometric and signomial programs are considered. 2.2 Posynomials and T h e i r Condensation A posynomial P^ Cx) i s a real-valued function defined as r k m •p. (x) = r c . n x. i j k k - i£i Jk i = l i (2.1) where x = [x^, x^, x ^ ] 1 . The exponents a ^ j ^ are a r b i t r a r y r e a l numbers while the c o e f f i c i e n t s C j ^ and the v a r i a b l e s x, are r e s t r i c t e d to being p o s i t i v e . Each term of (2.1) i s c a l l e d a monomial. The form of a posynomial appears i n many engineering design equations. Hence an optimization theory such as signomial programming that i s based on posynomials i s a powerful a i d to the designer i n h i s search for optimal designs. In general, a posynomial i s nonconvex i n the space of x. But i f t. = In x. and t = [ t 1 , t O J t ]', then l I — 1 2' ' m 7. x = exp(t.) 1 1 (2.2) and r k m \ ( t ) - c . k e x p C ^ a . . k t . ) (2.3) Since the exponential function is convex and a l l the c o e f f i -cients c ^ a r e positive, i t follows [1, p 54] that P^ Ct) i s convex i n the space of t^ , which is just the Euclidean space Rm. Hence any posy-nomial can be transformed into a convex function. A consequence of this convexity property of posynomials i s the bounding from below of any posynomial by a monomial. Let ^ ( t ) be defined as P v(t) = In P, (t) k ~ k ~ (2.4) It can be shown that P^ Ct' 1 S a continuously differentiable convex function. Then for a given point t > 0, \ ( t ) >P k(t) + I (t. - t.) 3P k(t)/3t ± (2.5) The inequality (2.5) is simply the definition of a continuously d i f f e r -entiable convex function [2, p 28]. Exponentiating both sides of (2.5) gives ^ <\, „ m ^ „ exp(P k(t)) > exp(Pk(t_) . n exp((t i - t ± ) 3P k(t)/3t ± (2.6) Substituting t ± = In x± and t± = In x ± simplifies (2.6) to m b ~ \(*> > P k(x) * . n/x./x.) x k P k ( i , x) ( 2 > 7 ) where = 3 P k ( t ) / 3 t = [(3P (x)/3x.)(x./P v(x))] 1 K- 1 it — 1 1 K. — X = X (2.8) 8. O b s e r v e t h a t P ( x , x) i s a m o n o m i a l . T h i s t e c h n i q u e o f a p p r o x -K. i m a t i n g a p o s y n o m i a l P k ( x ) a t x . = J£ w i t h a m o n o m i a l P k (x> £ ) i s c a l l e d t h e c o n d e n s a t i o n o f P , (x . ) a t x = x . P. ( x , & ) i s k n o w n a s t h e d o n d e n s e d p o s y -n o m i a l o f P^Cx) a t x-T h e r e l a t i o n s h i p b e t w e e n P^OO a n d P^ Cx.* g i v e n h y L e m m a 2 . 1 . L e m m a 2 . 1 A. L e t ^ ( x * 2E.) D e t n e m o n o m i a l c o n d e n s e d f r o m t h e p o s y n o m i a l P k ( x ) a t x . T h e n ( a ) P k ( x ) > i > (x, x ) f o r a l l x > 0 ( b ) P k ( £ ) = P k ( £ , x ) ( 2 . 9 ) ( c ) [ 3 P . / 3 x . ] - = [3P./3X.] ^ ' k x x = x k i x = 2E P a r t ( a ) o f L e m m a 2 . 1 f o l l o w s d i r e c t l y f r o m t h e d e f i n i t i o n o f t h e c o n d e n s e d p o s y n o m i a l . P a r t ( b ) c a n b e v e r i f i e d b y d i r e c t s u b s t i -t u t i o n . P a r t ( c ) c a n b e s h o w n b y d i f f e r e n t i a t i n g t h e r i g h t s i d e o f ( 2.7) a n d d i r e c t s u b s t i t u t i o n . A n o t h e r a p p r o a c h t o c o n d e n s i n g p o s y n o m i a l s i s t o u s e t h e g e o m e t r i c i n e q u a l i t y b a s e d o n t h e f o l l o w i n g l e m m a . L e m m a 2 . 2 I f u . > 0 a n d e.. _> 0 f o r j = 1 , 2 , . . . , N , t h e n ' j i ^ ^ A W 3 x* (2.io) where N x = I e ( 2 . 1 1 ) j = l 3 a n d ( u . / e j r J = 1 i f E j = 0 ( 2 . 1 2 ) Moreover t h e i n e q u a l i t y becomes a n e q u a l i t y i f a n d o n l y i f N N = 1 , 2 , N ( 2 . 1 3 ) =j I u = u I e j 1 J i = l J i = l 9. The proof of Lemma 2.2 may be found i n [3] Given a posynomial P v(x) defined by r k m a \ ^ = j l i U j k ( ^ ' u j k C - } = c j k ± S i y be u point &_ by s e t t i n g x i j k (2.14) Lemma 2.2 may be used to derive the condensed posynomial P^X—' JO a t t n e e . = u M (x) /P. (x) j 3k~ k v - y (2.15) and X = 1.0 The condensed posynomial ' jO becomes (2.16) p k ( x , x) - ek(i) x . ^ ^ ( 2 > 1 7 ) where r V*> " j l i V/3 (2.18) and k a i k ( ^ = £ j a i j k . (2.19) I t can be e a s i l y v e r i f i e d that (2.17) s a t i s f i e s Lemma 2.1. 2.3 Signomial Programming A signomial g^Oi) l s a real-valued f u n c t i o n defined as g^x) = P k(x) - 0 k(x) m . ... .., , „m (2.20) where x e R +~ A p o s i t i v e orthant of R"'. P k and Q k are posynomials of the form given i n (2.1). A signomial program(SP) with i n e q u a l i t y con-s t r a i n t s only i s the optimization problem defined as follows (SP) min x^ s.t. g k(x) < 1, k = 1, 2, .... p ( 2 - 2 1 ) x E R + r a A signomial program with a signomial objective function g 00 can be e a s i l y put in t o the format of (2.21) by adding the i n e q u a l i t y (g Q(x) + C) <_ x^ to the constraint set. C i s a constant added to insure the p o s i t i v i t y of the l e f t side f o r a l l f e a s i b l e x. Note that (2.21) also allows simple upper and lower bound c o n s t r a i n t s . Hence i n t h i s t h e s i s , the formulation of (2.21) i s considered as the standard problem statement of a signomial program. D i f f e r e n t but equivalent formulations had been introduced by Passy and Wilde [4], A v r i e l and Williams [5], and D u f f i n and Peterson [6]. Because the d i f f e r e n c e of two convex functions i s generally nonconvex, a signomial program i s nonconvex. The s o l u t i o n of a s i g -nomial program cannot therefore guarantee a g l o b a l minimum. Only l o c a l minima are assured under c e r t a i n r e g u l a r i t y conditions. Another con-sequence of the nonconvexity of signomial programs i s the loss of a strong d u a l i t y theory that i s associated with convex programs. There i s , however, a weak d u a l i t y theory r e l a t i n g the l o c a l s o lutions of (SP) (the primal problem) to those of a l i n e a r l y constrained dual problem. The p r e c i s e d e f i n i t i o n of the dual problem depends on how the primal signomial program i s formulated. The essence of the d u a l i t y theory, however, remains the same, namely under some dual constraint q u a l i f -i c a t i o n s a dual Kuhn-Tucker point can be derived from a known primal Kuhn-Tucker point, and v i c e versa. In e i t h e r case, the primal and dual o b j e c t i v e functions are equal when evaluated at the respective Kuhn-Tucker points. More d e t a i l s on the d u a l i t y properties of s i g -nomial programs may be found i n [4] - [6]. X I . 2.4 Geometric Programming If no negative term i s present i n the signomial program of (2.21), then the problem becomes a geometric program involving posy-nomials only. The problem statement may, for later convenience, be rephrased as (GP) min P Q(x) s.t. P k(x) 1 1.0, k = l , 2, ...,p (2.22) x e R+m where M K m a . . V*) = I u (x), u (x) = c n x. ^ ( 2 > 2 3 ) and MD = 1, ^ - + 1, N k = Mk_! + n k ^ ^ n k i s the number of terms in P k(x). (2.22) i s called the primal geometric program and x^ are the primal variables. If the logarithmic transformation t^ = l n x^ i s used and (2.22) i s expressed i n terms of t_^, i = 1, 2, ..., m, then from the discussion of Section 2.2, the minimization problem reduces to a convex program since P ^ t ) > k = 0, 1, p, are convex functions of t^ . By exploiting the monotone increasing property of the natural logarithmic function, the primal geometric program can be cast into the following transformed primal problem associated with the primal geo-metric program. min In P Q(t) < T G p) ~ (2.25) s.t. In P k(t) _S °> k = 1, 2, . .., p where N k m \ & • V ^ J i a ^ ) ] (2-26) and and are given by (2.24). Another important problem that can be associated with a given primal geometric program i s the dual geometric program stated as follows: 12. N 6 P A (DGP) max v(S) = 1^ ( ^ / o ^ \ k (2.27) N k s' t > V = I 6 + » k = 1, 2, (2.28) and N £ P a 6 = 0 , i = 1, 2, m (2.29) j=l 3 ' N I 6 i = 1 (2.30) j=l 2 •Sj 1 0, j (2.31) In this program a ^ , C ^ , and are defined i n the same way as i n the primal programs. (see (2.23) and (2.24).) The function v(6_) is called the dual function and the variables 6^ are the dual variables. The relations (2.29), (2.30), and (2.31) are known as, respectively, the orthogonality condition, the normality condition, and the positivity condition. Note that each dual variable 6_j i s linked to a monomial term in the primal problem (GP), while each X^ i s paired with a posynomial inequality constraint of (GP). It i s shown in [1] that lcgv(6_) i s concave. Since max log (v (6)) = max v(6) over any feasible set of _5_, with log v(5) replacing v(6) as the objective function, (DGP) i s a concave program with linear constraints. The signi-ficance of this transformation w i l l be seen in the following discussion of the duality properties of geometric programs. It has been shown that the solution of a geometric program (GP) i s equivalent to the solution of a convex program (TGP) derived from (GP). Since convex programs have strong duality properties, the primal problem (GP) is expected to satisfy a strong duality theorem. This indeed i s the case as shown by Theorem 2.1 due to Duffin, Peterson and Zener [1]. 13. A program (primal or dual) i s consistent i f there e x i s t s a vector s a t i s f y i n g the program's constr a i n t s . (GP) i s superconsistent i f there e x i s t s a vector x such that P, (&)< 1 , k = 1 , 2 , p, and * > o. Theorem 2 . 1 Suppose that a primal geometric program (GP) i s supercon-s i s t e n t and p 0 ( x ) i s minimum at x*> a f e a s i b l e vector of (GP). Let Fp and F^ be, r e s p e c t i v e l y , the f e a s i b l e sets of (GP) and i t s a s s o c i -ated dual program (DGP). Then ( i ) F D i s non-empty, and f o r a l l x £ F p and S_ e F^, P n(x) > min P Q(x) 4 P 0(x*) = v(6*) A max v(6) >_ v(6) ( 2 . 3 2 ) k * e F P A e F D ^ ( i i ) there e x i s t non-negative Lagrange m u l t i p l i e r s y, associated with P, ( J C * ) , k = 1 , 2 , p, such that * u ( X * ) / P Q ( X * ) , j= 1 , 2 , N Q ( 2 . 3 3 ) [ pj*u ( x * ) / P k ( x * ) f j = M K , N F C ; k = 1 , 2 , p ( 2 . 3 4 ) Furthermore, X k* = y k * / \ ( x * ) ( 2 . 3 5 ) ( i i i ) i f 5* i s the maximizer of (DGP), then x* i s the s o l u t i o n to any system of m l i n e a r l y independant equations s e l e c t e d from the following equations: m I a ± j log x. = log (6 * v ( J * ) c ), j = 1 , 2 , N ( 2 > 3 6 ) i = l J J l^a±i log x. = log ( M / c . y ) , j = M K , N k , ( 2 * 3 7 ) It. - 1 ) 2 y 9 • • y P 14. Theorem 2.1 has two important practical implications. From a computational point of view, i t is indeed attractive that the solution to the primal program can be obtained via the solution of a concave program with linear constraints, and the solution of a system of linear algebraic equations. From an engineering design viewpoint, (2.32) i s appealing since i t allows the designer to use the inequality PQ(x.) 21 v ( § ) for x E Fp and 5_ e to decide whether the current feasible design i s acceptable, depending on the gap P Q(x) - v(6). 2.5 Solution of Geometric Programs Theorem 2.1 offers two options for solving a geometric program -a direct primal solution or one via the dual problem (DGP). The main advantage of the dual method is the structure of (DGP) - the maximization of a concave function over a linear manifold. Algorithms tailored for this class of nonlinear programs may be used. (e.g. see references [7] -[10].) A prime consideration in the dual approach i s the linear manifold's dimension usually called the degree of d i f f i c u l t y d. Since there are dual variables and m + 1 dual linear equality constraints, d = - (m + 1) for the case when the m x N (m < N ) exponent matrix A = [a. .1 has f u l l P P rank. If the rank of A i s less than m, the primal variables can be rede-fined to insure that A has f u l l rank [1, pp 82-83]. If d = 0, the dual optimal solution &f< i s conveniently obtained as the unique solution of a system of linear algebraic equations. But i f d > 0, iterative methods must be used. Employing the dual approach has Its d i f f i c u l t i e s . The degree of d i f f i c u l t y d can be quite large, although the primal space has low dimensionality. Hence the use of the dual method involves a trade off 15. between the case of linear constraints and the possible burden of high dimensionality. Another serious potential d i f f i c u l t y i s due to the existence of slack primal constraints. From the necessary complementary slackness condition for optimality, the optimal Lagrange multiplier u^* associated with the kth slack primal, inequality constraint must be zero. Hence the kth set of equations defined by (2.37) cannot be used. It- i s possible that not enough equations are available to yie l d the optimal primal vector x*. To overcome the d i f f i c u l t y posed by slack primal constraints, Kochenberger [11] proposed the addition of slack variables to loose constraints. The drawback of this method i s that the loose constraints have to be known a p r i o r i or else slack variables need to be added to each constraint. In the latter case, both the dimension of the primal problem and the degree of d i f f i c u l t y are i n -creased. A primal approach would solve the transformed primal geometric program (TGP) given i n (2.25). The motivation i s obvious. While (GP) is generally nonconvex, i t s equivalent problem (TGP) i s a nice convex program that can be solved by any of the algorithms developed for convex programs. For example, Avriel et a l [12] has applied a cutting plane algorithm to solve geometric programs. Beside exploiting convexity, the primal approach can handle loose constraints with no d i f f i c u l t y . Simple upper and lower bound constraints on x can also be easily taken into account, while the same type of constraints would mean increased degree of d i f f i c u l t y i n a dual method. 2.6 Solution of Signomial Programs Signomial programs are nonconvex programs and their solutions are usually at best local minima only. As in the case of geometric pro-grams, the solution to a signomial program can be obtained by a primal or a dual method. The dual approach solves a dual program that s t i l l has only linear constraints, but the logarithm of i t s dual objec-tive function is not concave [4] [6]. Furthermore, the dual solution no longer satisfies a strong duality theory. However, the potential d i f f i c u l t i e s of a dual approach, as discussed i n Section 2.5, remain. The primal approach to solving signomial programs i s based on convex approximation. In l i e u of solving the original nonconvex program, a sequence of convex programs i s solved such that the cor-responding sequence of solutions approaches a local minimum of the non-convex program. Charnes and Cooper [13] f i r s t suggested such a comp-utational strategy for signomial programs. Later several authors [5], [14], [15] adopted the same strategy but differed in the specifics of how to convexify the original nonconvex program. A practical algorithm that has been used for this thesis is that proposed by Avriel and Williams [5] and further developed by Avriel et a l [12]. A computer code implementing the algorithm was written by Dembo [16]. The Avriel - Williams algorithm solves the following sig-nomial program ( S « ' (2.38) s. t. x_ e X where x = [x]_, x 2 ', X = fx : x e Rm; P k(x) - Q k(x) <_ 1, k = 1,2,... ,p; 0 < x L <. x < x11} Note that (2.38) i s almost identical i n formulation to (2.21) except in the former the simple upper and lower bounds on x. are e x p l i c i t l y stated. Consider the kth signomial inequality constraint i n the feasible set X. The constraint can be rewritten as P k(x)/(1 + Qk(x)) < 1 (2.39) 1 7 . (2.39) i s elearly a ratio of two posynomials. Assume that a point x ( i _ 1 ) £ X is known. The denominator (1 + Qk (x)) can be condensed at ^ "^to yield a monomial Qk(x_5 ~ ^ ) - If the denominator of each constraint in X is replaced by their respective condensed posynomial, the following problem, called (GP^^), results. ( G P ^ ) min x^^ s.t. P k(x)/Q kU. 2£ ( i _ : L )) 11. k = 1, 2, p C2.40) 0 < x L 1 x < x U The above problem has the following properties: (a) It i s a geometric program since a posynomial divided by a monomial remains a posynomial. (b) Its feasible set is contained in the feasible set X of (2.38) because by (2.9) P k(x)/(1 + Q k(x)) < Pk(x)/Qk(2£» 1 (2.41) tor any _x feasible for (GP ( l )) (c) By (2.41), i t s optimal solution i s feasible, but not necessarily optimal, for (SP) (2.38). The Avriel - Williams algorithm solves a sequence of problems of the form specified by (2.40). It is shown in [5] that the algorithm produces a sequence of points converging to a local minimum of the original problem (SP), provided that the signomial program has a regular feasible set. A feasible set X is regular i f i t i s non-empty and compact, and the gradients of i t s tight constraints generate a pointed cone, .(i.e., the origin i convex h u l l of the gradients mentioned). In summary, the algorithm may be stated in the following manner Algorithm 2.1 Step 0: Given x ^ e X and let x ^ be the solution to G P ^ . Set i = 1 18. Step 1: Construct G P ^ according to (2.40). Step 2: Solve G P ^ and obtain the point . Step 3; Stop i f the convergence c r i t e r i o n i s s a t i s f i e d . Otherwise, i = i + 1 and go to Step 1. I t i s worth noting that t h i s algorithm may use any primal or dual algorithm f o r s o l v i n g the geometric programs. The A v r i e l -Williams algorithm i s a primal method only with respect to the nonconvex program (SP). I I I . SIGNOMIAL PROGRAMS WITH EQUALITY CONSTRAINTS 3.1 Motivation Since i t s f i r s t introduction in 1961, geometric programming, and later i t s generalization signomial programming, have been formulated in terms of inequality constraints only. Such a formulation allows the dev-elopment of a duality theory that i s , in the case of geometric programming, elegant and computationally attractive. However, an examination of many engineering design problems shows that there is a strong need for ex-tending the problem formulation to include e x p l i c i t l y equality constraints i n the f e a s i b i l i t y set. From the theoretical and computational viewpoints, there are also compelling reasons for an algorithm that can handle equality constraints and s t i l l retain the basic structure of the original signomial programming formulation. The presence of equality constraints in engineering design can be attributed to numerous reasons. Some of these are: 1) the conservation laws of mass, energy, momentum, and charge, 2) the input-output relation of a transformation, such as a stage in a process, 3) the need to main-tain equilibrium conditions, original topology, etc., 4) the relation of a variable in the cost function to other design parameters, 5) the need to satisfy boundary conditions. Some specific examples that are amenable to the signomial programming formulation and that i l l u s t r a t e well some of the mentioned reasons are chemical reactor design [17], steady-state process optimization [18], and optimal structure design [19]. A nonlinear program with nonlinear equality constraints is always nonconvex. In the case of geometric programming, the strong duality relations of Theorem 2.1 no longer apply. In fact, i t can easily be shown t h a t a g e o m e t r i c p r o g r a m w i t h e q u a l i t y c o n s t r a i n t s i s e q u i v a l e n t t o a n i n -e q u a l i t y - c o n s t r a i n e d s i g n o m i a l p r o g r a m , j u s t a s a r e s i g n o m i a l p r o g r a m s w i t h e q u a l i t y c o n s t r a i n t s . T h e o r e t i c a l l y , t h e s e e q u i v a l e n t s i g n o m i a l p r o -g r a m s c a n b e s o l v e d b y a l g o r i t h m s f o r s o l v i n g s i g n o m i a l p r o g r a m s . B u t t h i s a p p r o a c h r e q u i r e s t h e r e p l a c e m e n t o f e a c h e q u a l i t y c o n s t r a i n t g^(x)=l w i t h t h e p a i r o f i n e q u a l i t i e s g k ( x ) < £ a n d g k ( x ) r l x * ^ h e m e t h o d h a s t w o c o m p u t a t i o n a l d r a w b a c k s : 1 ) t h e s i z e o f t h e p r o b l e m i s i n c r e a s e d , 2 ) t h e p r i m a l f e a s i b l e s e t h a s n o i n t e r i o r a n d h e n c e i s n o t r e g u l a r . A n o t h e r p o s s i b l e s c h e m e i s t o r e p l a c e e a c h e q u a l i t y c o n s t r a i n t g k ( x ) = l w i t h e i t h e r g ( x ) < l o r g . ( x ) > l . T h e f i r s t o b j e c t i o n t o t h i s s c h e m e i s t h a t i f t h e n u m b e r q o f e q u a l i t y c o n s t r a i n t s i s l a r g e , t h e r e a r e 2 ^ c o m b i n a t i o n s t o c o n s i d e r , a n d t h e t r i a l - a n d - e r r o r c o m p u t a t i o n c a n b e p r o h i b i t i v e u n l e s s s o m e h e u r i s t i c s e a r c h p l a n i s f o l l o w e d . A m o r e i m p o r t a n t o b j e c t i o n t o t h i s s c h e m e i s t h a t t h e r e i s n o g u a r a n t e e t h a t t h e o r i g i n a l s o l u t i o n w i l l b e o b t a i n e d . I t i s c o n c e i v a b l e , f o r e x a m p l e , t h a t e a c h o f t h e 2 ^ c o m b i n -a t i o n s h a s a s l a c k c o n s t r a i n t , y e t t h e s o l u t i o n s o u g h t r e q u i r e s a l l t h e c o n s t r a i n t s t o b e t i g h t s i n c e t h e y a r e e q u a l i t y c o n s t r a i n t s . I n t h i s c h a p t e r a n a l g o r i t h m i s d e v e l o p e d t o h a n d l e s i g n o m i a l e q u a l i t y c o n s t r a i n t s w i t h o u t t h e r e p l a c e m e n t b y e q u i v a l e n t i n e q u a l i t y c o n -s t r a i n t s . T h e a l g o r i t h m i n c o r p o r a t e s t h e e q u a l i t y c o n s t r a i n t s i n t o t h e o b j e c t i v e f u n c t i o n a n d s o l v e s a s e q u e n c e o f i n e q u a l i t y - c o n s t r a i n e d s i g -n o m i a l p r o g r a m s i n l i e u o f t h e g i v e n p r o b l e m . T h e o r i g i n a l i n e q u a l i t y c o n s t r a i n t s a r e r e t a i n e d a s t h e y a r e . 3 . 2 P r e v i o u s W o r k T h e p u b l i s h e d l i t e r a t u r e o n s i g n o m i a l p r o g r a m m i n g c o n t a i n s v e r y l i m i t e d w o r k o n t h e s o l u t i o n o f s i g n o m i a l p r o g r a m s w i t h e q u a l i t y c o n s t r a i n t s . Blau and Wilde [ 1 7 ] reported an a p p l i c a t i o n of signomial programming i n v o l v i n g e q u a l i t y c o n s t r a i n t s . They, however, used p h y s i c a l reasoning and some algebra to transform the c o n s t r a i n t s i n t o i n e q u a l i t i e s . The same approach was adopted by R i j c k a e r t [ 1 8 ] . In applying signomial programming to nonlinear assignment problems, Passy [ 2 0 ] used the stand-ard technique of s u b s t i t u t i n g each equality c o n s t r a i n t with a p a i r of i n -e q u a l i t i e s opposite i n sense. Later, Blau and Wilde [ 2 1 ] considered s i g -nomial programs with e q u a l i t y constraints only. They proposed a Newton-Raphson method to solve f o r the stationary points of the primal's Lagran-gian. I t seems that the algorithm reported i n t h i s t h e s i s and i n [ 2 2 ] i s the f i r s t e f f o r t to solve signomial programs with mixed i n e q u a l i t y and e q u a l i t y c o n s t r a i n t s , i n a manner that preserves the key properties of signomials without r e q u i r i n g any constraint transformation by the user. 3 . 3 Problem Statement The problem that i s of i n t e r e s t i n t h i s chapter may be stated formally as follows (SPE) min X l . s . t . x e X (3 h k(x) A g p + k ( x ) - 1 = 0 , k = 1 , 2 , q where X = {x : x e R m; g k(x) 1 1 , k = 1 , 2 , p; 0 < x L <^ x £ x U} A l l the functions g : R+m-> R are signomials as defined i n ( 2 . 2 0 ) . The d i f ference between problem (SP) ( 2 . 3 8 ) and (SPE) ( 3 . 1 ) i s s o l e l y due to the q signomial equality c o n s t r a i n t s . The next s e c t i o n d e t a i l s the development of an algorithm that solves (SPE) v i a a sequence of (SP)-type problems which can be conveniently solved by the various algorithms r e f e r r e d to i n Section 2 . 6 . 3.4 Proposed Algorithm 3.4.1 Background Constraints i n nonlinear programming are generally treated i n four ways, each being the basis of a family of algorithms. The f i r s t i s to transform the constrained problem into an unconstrained one by i n c o r -porating a l l the constraints i n t o the objective function. While such an approach e x p l o i t s the e f f i c i e n c y and ease of unconstrained optimization algorithms, i t s major disadvantage i s that whatever s p e c i a l structure the c o n s t r a i n t s may have, i s l o s t . The second strategy i s the primal f e a s i b l e method i n which unconstrained optimization algorithms are modified to i n -sure that the sequence of points generated are a l l f e a s i b l e f o r the given problem. The primal f e a s i b l e method works w e l l with l i n e a r c o n s t r a i n t s , but with nonlinear constraints, i t i s plagued with the d i f f i c u l t i e s of g e t t i n g an i n i t i a l f e a s i b l e point and remaining w i t h i n the f e a s i b l e set. Approximation leading to the use of l i n e a r programming or simplex-type operations i s the t h i r d way to solve constrained nonlinear programs. This approach can be quite e f f i c i e n t , depending on how w e l l the f e a s i b l e set i s approximated. The fourth approach i s the Lagrangian or dual method which i s motivated by the viewpoint that the Lagrange m u l t i p l i e r s are the fundamental unknowns associated with a constrained optimization problem. Hence, the method does not solve d i r e c t l y the o r i g i n a l problem but instead attacks an alternate problem, the dual problem, whose unknowns are the Lagrange m u l t i p l i e r s of the o r i g i n a l problem. In almost a l l the algorithms of each family, both e q u a l i t y and i n e q u a l i t y constraints are treated i n the same manner. They are a l l included i n the objective function of a new problem, or a l l are used to maintain f e a s i b i l i t y , or a l l are approximated i n some form, or a l l are dualized to define the dual problem. The algorithm discussed in this chapter is an exception to this pattern of handling identically both equality and inequal-it y constraints. The algorithm is designed to be a synthesis of the four basic solution strategies in order to preserve the useful properties of a signomial program. 3.4.2 Development of the Algorithm The d i f f i c u l t y of Problem (SPE) defined in (3.1) is posed by the signomial equality constraints (3.1b). Without these constraints the problem simply reduces to a signomial program that can be conveniently s o l -ved by a variety of algorithms discussed in Section 2.6. It is clear then that the inequality constraints of (3.1) should remain untouched. The immediate question that has to be confronted i s : How should the equality constraints be manipulated? The techniques discussed in Section 3.1 have been ruled out for the reasons mentioned there. An examination of the basic solution strategies in nonlinear programming suggests that the trans-formation method, otherwise known as the penalty function method, i s perhaps the suitable approach. Using an external penalty formulation [23], the problem (SPE) can be rewritten as min X ; L + r ( 1 ) £ |h k(x)| B k = 1 (3.2) s. t. x e X In (3.2), 3 _> 1 and the sequence { r ^ } is an unbounded increasing positive sequence. If 3 is an even integer ( i t i s usually 2), then the objective func-tion of (3.2) is a signomial since any even integral power of a signomial is i t -self a signomial. Thus in this case (3.2) i s a signomial program as defined as problem (SP) . Let x ^ be the solution to the i t h subproblem defined by (3.2). Also let be the minimum value of x^ in (SPE). The convergence of the sequence {x^"^ } to a local minimum of (SPE) is obtained according to the following result due to Fiacco and McCormick [23]: If the local minima of problem (SPE) forma nonempty compact set A and r ^ — * • « • as i—•-<*>, then there exists a compact set S such that A <zz interior of S and (i) * for sufficiently large i , x e interior of S. Moreover, x^—»-x^ , and every limit point of any convergent subsequence of {x_^} is in A, While the convergence of the exterior penalty formulation is satisfactory, the transformation, indeed a l l penalty function methods, possess the unfavorable property of exhibiting ill-conditioned Hessians. It is a fundamental prop-erty [24] of penalty function methods that as r ^ — * • «, the Hessian of the transformed objective function is equal to the sum of the Hessian of the original problem's Lagrangian and a matrix of rank R, whose eigenvalues approach i n f i n i t y . The number R gives the number of active constraints. To overcome the numerical d i f f i c u l t i e s caused by the ill-conditioning of the Hessian as r^"""^-*» °°, the requirement that r ^ — * - «° in order to achieve convergence must be relaxed. This goal is realized by the method of multipliers f i r s t independently proposed by Hestenes [25] and Powell [26]. Because the method i s a modified penalty function method, i t preserves the attractive property that the transformed problem is an inequality-constrained signomial program. In i t s original formulation the method of multipliers solves the nonlinear program min f (x) o — s.t. f k(x) = 0, k = 1, 2, q ( 3 . 3 ) x e R by solving a sequence of unconstrained subproblems each defined as 25. •3 q min Hx, X ( i ) , K ( i ± ) ) = f (x) + £ A ^ f , (x) + K ( i ) ]> f, 2(x) ° k=l k k=l s.t. x e R m A e (3.4) It i s understood that X_ = [ X ^ X ^ ...,X ]'. In (3.4), \ ^ is the i t h * (i) estimate of the optimal Lagrange multipliers X^ , k = 1,2, ...,q and K is a sufficiently large but f i n i t e positive constant. The general outline of the algorithm i s as follows: Algorithm 3.1 Step 1 Set i = 1. Assume X _ ^ and are specified. Step 2 Solve (3.4). Let the solution be x ^ . Step 3 If ( x ^ ) =0, #k, stop. Otherwise, set i = i +1. Update X _ ^ and/or K . ^ \ and go to step 2. The updating rules for X^"^ and w i l l be discussed later. Conditions for the convergence of the algorithms w i l l be considered i n the next sec-tion. Now what is of interest i s the adaptation of the method of multi-pliers to accomodate inequality constraints. Besides being considered as a modified penalty method, the method of multipliers can also be viewed as a modified primal-dual method. This view was sketched by Luenberger [24] and more recently discussed in detail by Bertsekas [27]. Several other researchers (e.g. [28] - [32]) have adopted the same point of view to study the theoretical properties of the method of multipliers. Consider the nonlinear program given by (3.3). An equivalent problem, in the sense of having the same solution and the same minimum objective function, is the following problem min f Q(x) + K 1^ fk2O0 q I k=l (3.5) s.t. f (x) = 0, k = 1, 2, ..., q x E R m for any K > 0. The Lagrangian of (3.5) i s q q Z(x, X, K) = f D(x) + J x X kf k(x) + K J i f k 2 ( x ) (3.6) Let x be a local minimum of the original problem (3.3). Associated with this solution i s the optimal Lagrange multiplier _X . Then for suff i c i e n t l y large but f i n i t e K, say K >_ K for some K , the Hessian of J2.(x,X_,K) at * * ft (x ,X^ ) i s positive definite. As a result, for every K >_ K , problem (3.5) has a local l y convex structure, and local duality theory i s applicable. Let the function'd: Rq-^»-R be defined as d^CX) = min S,(x, X_, K) x e R m (3.6) Then themethod of multipliers can be cast as a primal-dual algorithm alter-nating between the primal problem (3.6) and the dual problem defined as max d k(\) (3-7) 2.eRq Note that both problems are unconstrained. Now suppose that a set of i n -equality constraints 6j(x) ^ . 0 » j = 1»2, p, i s appended to the feasible set of (3.5). These additional constraints can be easily handled by the concept of partial duality [24], [32]. The definition of the dual function cL^X) need not include the Lagrange multipliers of a l l the primal constraints. If local convexity applies, dualization can be with respect to any subset o f the primal constraints. Hence, the dual function ^(X) c a n be redefined as <L (X) = max £(x, X, K) X E X , " ~ < 3 ' 8 > where X•= {x : x e Rm, 6 ^ (x) <_ 0, j = 1, 2 , p} The dual problem i s the same as that of (3.7). Since the method of multi-pliers i s a primal-dual method, i t follows that one way to apply the method of multipliers to nonlinear programs with inequality constraints i s simply to use pa r t i a l duality and incorporate the inequality constraints into the definition of the dual function. As a result, the original problem is replaced with a sequence of less constrained subproblems. Clearly such an approach i s advantageous only i f the subproblem's constraints exhibit a special structure worthy of special consideration. This i s precisely the case with signomial programs with equality constraints. Applying the method of multipliers to the signomial program stated i n (3.1) yields the following i t h subproblem q q ( T T , ( I ) ) min x. + I A, ( i )h, (x) + K ( i ) £ h 2 ( x ) 1 1 k=l k=l s.t. x e X (3.9) where X = {x : x E Rm; g k(x) <_ 1, k = 1, 2, .. ., p; 0 < x L <_ x < x u} It can be seen that the minimization problem (3.9) is an inequality-con-strained signomial program. If a sufficiently large constant i s added to the objective function, (3.9) can be easily cast into the standard format of (2.21) for solution by the Avriel-Williams algorithm. 3.4.3 Statement of the Algorithm The algorithm that has just been developed for solving signomial programs with equality constraints may now be formally stated. Given the problem (SPE) posed in (3.1) a solution can be ob-tained iteratively by the following algorithm: Algorithm 3.2 Step 1 Set i = 1. Assume that arid are specified. Obtain a point x/°^ that i s feasible for problem (ir^-)) as stated in (3.9). Step 2 Solve (T T^ 1)) by the Avriel-Williams algorithm with x^"1^ as the starting point. Let the solution be x ^ • Step 3 If the equality constraints h ^ C x ^ ) = 0, k = 1, 2 , . . . , q, are satisfied, stop. Otherwise, i = i + 1. Update A^"^ and/or Go to Step 2. The parameters A_ and K may be updated in various ways, and the updating rules w i l l be discussed later. The operation of the proposed algorithm may be illustrated by considering the following simple signomial program with both equality and inequality constraints. 2 min x^ (3.10) s.t. 3 - x ^ - x 2 < _ l x 2 = l 0 < x 1 > x 2 Assuming that K i s sufficiently large and fixed, the kth subproblem i s min x 2 + ( A ( k ) - 2K) x 2 + Kx^ (3.11) s . t . 3 - x ^ - x 2 £ l 0 < x^, x2 From the Kuhn-Tucker conditions of the subproblem, the solution to (3.11) ± S (k) 2K + A ( k ) (k) 2K + 3 - A ( k ) X l 2(K + 1) » x2 2(K + 1) ( 3 < 1 2 ) One effective updating rule for A_ i s the Hestenes-Powell formula [25], [26] given as A_ ( k + 1 ) = A^ k ) + 2Kh(x ( k^). Using this formula yields X(k+D = ( 1 ) x(k) + 2K K + 1 + K + 1 ( 3 ' 1 3 ) = a A ( k ) + b, a, b > 0 From the theory of difference equations, the closed-form solution to (3.13) is k+1 i 1 - a A ( k ) = c a k + b ( - v a" ) (3.14) for some arbitrary constant c. Set A^°^ = 0. Then c = -b and 29. X ( k ) = 2[1 - ( K I 1 ) k ] (3.15) I t f o l l o w s t h a t x . ( k ) = 1 - 1 / ( K + l ) k + 1 , x 0 0 = 1 + 1 / ( K + l ) k + 1 (3.16) A s k °», c o n v e r g e n c e t o t h e u n i q u e - s o l u t i o n ( x ^ = 1, x^ = 1) o c c u r s f o r a l l K > 0 . H e n c e , K n e e d n o t a p p r o a c h i n f i n i t y . H o w e v e r , t h e r a t e o f c o n v e r g e n c e i n c r e a s e s r a p i d l y w i t h m o d e s t i n c r e a s e i n K . T h i s f a c t c a n b e v e r i f i e d b y s u b s t i t u t i n g v a l u e s o f K i n t o (3.15) - (3.16). I f t h e o b j e c t i v e f u n c t i o n i s j u s t x ^ , i t c a n b e s h o w n t h a t t h e r e i s a o n e - s t e p c o n v e r g e n c e i n d e p e n d e n t o f t h e i n i t i a l v a l u e o f A.* T h e g e n e r a l c o n v e r g -e n c e p r o p e r t i e s o f A l g o r i t h m 3.2 a r e d i s c u s s e d i n t h e n e x t s e c t i o n . 3.5 C o n v e r g e n c e C o n s i d e r a t i o n s T h e p r o p o s e d a l g o r i t h m i s a s p e c i a l i z a t i o n o f t h e m e t h o d o f m u l t i p l i e r s t o s i g n o m i a l p r o g r a m s w i t h e q u a l i t y c o n s t r a i n t s . T h e c o n -v e r g e n c e o f t h e a l g o r i t h m i s t h e r e f o r e d e t e r m i n e d b y t h e c o n v e r g e n c e o f t h e m e t h o d o f m u l t i p l i e r s . I n t h i s s e c t i o n , t h e t h e o r e t i c a l r e s u l t s p e r t a i n i n g t o t h e c o n v e r g e n c e p r o p e r t i e s o f t h e m e t h o d o f m u l t i p l i e r s a r e s t a t e d . T h e t h e o r e m s q u o t e d h e r e a r e d u e t o P o l y a k a n d T r e t ' y a k o v [34], a l t h o u g h r e c e n t l y B e r t s e k a s [35] i n d e p e n d e n t l y p r o v e d s i m i l a r r e -s u l t s . T h e p r o b l e m o f i n t e r e s t i s m i n f ( x ) o — s . t . h ^ x ) = 0 , k = 1, 2, q ( 3 . 1 7 ) x e X c r R m w h e r e f a n d h , a r e f u n c t i o n s m a p p i n g X t o R. N o t e t h a t i n (3.17) x i s r e s t r i c t e d t o X . T h e a u g m e n t e d L a g r a n g i a n £ ( x , _X> K ) c a n b e c o n s t r u c t e d a s £ ( x , X, K ) = f Q ( x ) + X' h ( x ) + K | | h ( x ) | | 2 ( 3 . 1 8 ) where h(x) = r ^ ) , .... h ^ x ) ] ' and | | . | | is the Euclidean norm. Then a primal-dual pair can be defined as follows (P) min £(x, X, K) s.t. x E X (3.19) and (D) max d^CX) (3.20) s.t. A e Rq where the dual functional i s defined as d v(A) = min £(x, X, K) (3.21) S . t. X E X Suppose that at (x*, A*), x* i s a local minimum of (3.17) and together with Xj* satisfies the second-order sufficient optimality conditions. 2 2 Also assume that the Hessian matrices V f(x) and V h^(x) are Lipschitz i n the neighborhood of x*. Then the following two theorems are v a l i d . Theorem 3.1 For every A_ in the bounded set Y = {A_ : | |X_ - X* \ \ <_ p}, there exists a number K*(p) such that for K >_ K*(p), the following are true: a) (P) has a unique local minimum x in the neighborhood of x*; b) &(x, >L» K ) i s locally convex about x*; c) For some scalar c > 0, o | | x - x * [ | < C q I I A. - A*| |/K, # K > K* (3.22) 111 + % ( x ) - X*j I £ c J |_A - X*\|/K, ¥• K > K* (3.23) Theorem 3.2 For every p, there exists a number K*(p) such that i f K > K*, the dual functional dR(A_) is twice dif ferentiable and s t r i c t l y concave with respect to X . Furthermore, the f i r s t and second derivatives of d K(A) are given by 31. V d K a ) = h(x) (3.24) V 2 d R(X) = -J(x) [V 2 £(x, X, K)] J(x) (3.25) where J is the Jacobian of the equality constraints. Theorem 3.1 proves not only the existence of local convexity at x*> but also the convergence of x to x* as A + 1*. Furthermore, the rate of convergence i s estimated by (3.22) - (3.23). If K is f i n i t e * a linear rate i s i n effect. If K -> 00, the rate i s superlinear. Theorem 3.2 gives the di f f e r e n t i a b i l i t y properties of the dual functional d^CA). The theorem also bears relevance to the choice of the updating rule for X_. This point is discussed further i n Chapter V. The convergence of the overall algorithm relies on reaching a local minimum of the primal problem (P). This i s equivalent to the ex-ecution of Step 2 of the proposed algorithm i n which the Avriel-Williams algorithm i s used. That the latter algorithm yields a local minimum was proved by Avriel and Williams [5] using Zangwill's Convergence Theorem with some regularity condition imposed on the inequality constraints. It follows that subject to the suitable assumptions, the proposed algo-rithm converges. 3 . 6 Solution of the Subproblem The proposed algorithm, as stated in Section 3.4.3, leaves un-specified two major steps: the solution of the subproblem (3.9) and the updating formulas for _A and K. In this section the discussion i s focused on some ways of solving the subproblem. Considerations on how A_ and K should be updated are deferred to Chapter V. A significant property of the signomial program defined by (3.9) i s that the program is l i k e l y to have a large number of terms. Specifically, q i f g p + k ( x ) has Nfc terms, the subproblem (3.9) has N k(N k+l)/2 terms more 32. than the original problem (3.1). A large number of primal terms means that the degree of d i f f i c u l t y is large, or equivalently that the dual problem's dimensionality i s high. For example, a three-variable geometric program with a monomial objective function and two three-term equality constraints would have a signomial subproblem whose dual problem has a dimensionality of 16. In view of the potential d i f f i c u l t y posed by the dual problem's high dimension-a l i t y a primal approach has been selected to solve the subproblem (3.9). expansion of the quadratic terms or the avoidance of such expansion by using some indirect method. The former course of action i s tedious and entails cumbersome preparation work for data entry. It would indeed be de-sirable to just have to enter the terms as stated in the equality constraints. In the rest of this section, three techniques which satisfy this requirement and which have been numerically tested are discussed and compared. The f i r s t technique, referred to as Method I, is the condensa-tion scheme of Avriel and Gurovich [36] for algebraic programs (programs involving algebraic functions of signomials). Consider problem (IT, ^ ) . Applying a primal algorithm to (3.9) would require either the L e t \ ± I v ^ 1 = ' V k ( ^ - V k ( ^ l -Then the following problem is equivalent to (3.9) min x. s . t . (3.26) (b) 0 < z k = 1, 2, ••••». q (c) x e X 33. Constraint (3.26a) (omitting the index subscript) i s equivalent to the following pair of inequalities: P(x)/G(x, z) A P(x)/(Q(x) + z) < 1 , Q(x)/H(x, z) A Q(x)/(P(x) + z) < 1 Condensing G and H at (x, z) gives m P(x)/G(x, z) <P(x)/(G(£, z) n (x./x.) a3 (z/z)" 1^ 1) < 1 j=l 2 2 Q(x)/H(x,: z) < Q(x)/(H(x, £) n (x./x.) 8J (z/2) B m f*) < 1 j=l J 2 where a. = (1/G(x, z)) [x. 3G/9x.] - , 3 ~ 3 3 x = x 3. - (1/H(x, z)) [x. 3H/3x.] J - 3 3 x = x j - 1» 2, ... , m am+l = * / G & Sm+1 = z/H(x, z) (3.27) (3.28) (3.29) Since and 8^-^ a r e both positive, the rlgh_t inequalities of (3.28) can be written as the following single inequality: max { E 2 / 0 W l , F 2/ em+l } < z 2 (3.30) From (3.28), satisfying (3.30) implies that (3.26a) is satisfied. Further-more, (3.30) has to be tight at the minimum. Hence, (3.26) i s solved by solving a sequence of problems each of the form min ion s.t. (a) q u co2 + I max {E. 2/°fc,nrf-l, F. 2 / \ ^ } < co0 1 k=l k k q (b) x± + C ( i ) + I - 2K ( i )) g + k U ) £ K ( i ) co2 (3.33) k=l " (c) x e x» u i > °> u o > 0 In (3.33) I O q and to^ are positive auxiliary variables satisfying (3.33a) and (3.33b). They are introduced to put the subproblem into the standard format and to eliminate the parameters A^"^ and from constraint (3.33a). i s a sufficiently positive constant to insure that the l e f t side of (3.33b) i s positive. The minimization problem posed in (3.33) i s a piecewise "pseudo" signomial program. One way to solve i t i s as follows: (i) Use the scheme proposed in [36] to condense (3.33a) to a monomial, ( i i ) Simultaneously condense the signomials in (3.33b) and (3.33c) to posynomials; (ii i ) S o l v e the resulting geometric program, (iv) Repeat steps (i) - ( i i i ) u n t i l some convergence criterion is satisfied. The advantage of this technique i s that the size of the original problem remains the same. A major weakness is the numerical d i f f i c u l t y that can result from either or $ m + 1 being very small. When this 35. condition occurs, exponent overflow or gross clipping errors destroy the convergence of the iterative solution. The second method, named Method II, exploits the quadratic character of the penalty terms in the subproblem's objective function. The subproblem (3.9) i s recast as follows: (^2) min u Q p+k* q s.t. (a) 4 + ^ 4*<2> ± "0 (b) X ; L + C + I (A k - 2 K ) g p + k ( x ) < K c o 2 (3.34) ( c) x e X , o>1 > 0, CDQ > 0 where C i s large enough to make the l e f t side of (3.34b) positive. Since a primal approach to solving problem (T^ ) has been adopted, the positive and negative terms of a l l the signomial inequality constraints must be d i s t i n -guished. Such a distinction among the terms in (3.34b) and (3.34c) i s straightforward. In (3.34a) because of the quadratic exponents, the separation of positive and negative terms can also be easily carried out by replacing gp + k(x) with Pp + k(x) - Q + k 0 0 a n c* expanding the squares. Then the constraint (3.12a) can be written as 4+ 1 f p p + k ( ^ + Q p + k ^ ) ] - 2 V k ^ v ^ 1 "0 ( 3 - 3 5 ) k=l k _ i 2 Being posynomials, P + k ( x ) and Q p + k(x) are positive for x e X. Hence, 10^ and the f i r s t summation constitute the positive terms while the rest are the negative terms. The specific primal algorithm selected to solve (3.34) i s the Avriel-Williams algorithm explained in Section 2.6. The algorithm requires that each signomial inequality constraint of (3.34) be rewritten as a ratio of positive terms to one plus the negative terms, with the ratio bounded from above by unity. The solution to problem (TT^) (3.34) is then obtained by solving a sequence of geometric programs each of which is obtained by condensing the denominator of each ratio constraint at a feasible point of ( T T 2 ) . The condensation of the ratio version of (3.35) merits special attention. Written as a ratio of posynomials, (3.35) becomes ^ < 1 (3.3. If condensation is viewed as an application of the arithmetic-geometric inequality, then the multiplication in each squared or product term has to be carried out. But this i s the precise procedure that i s unwanted. How-ever, i f condensation is interpreted as exponentiation following the f i r s t -order Taylor approximation of the logarithm of a posynomial, then either the numerator or the denominator of (3.36) can be condensed directly in the manner of (2.7) - (2.8). The geometric programs approximating the problem (T^) have been solved by Dembo's cutting-plane method [12], [16]. The numerical exper-ience acquired suggests that convergence is attained rather slowly. A closer analysis of the results shows that a good part of the computing time i s spent on approximating the constraint (3.34a). This discovery leads to the adoption of the third method detailed in Chapter V following the dis-cussion of an intimately related algorithm in the next chapter. Method III differs from Method II in the following ways: 1) The objective function in (3.9) is not incorporated into the constraint set, as i t i s done in (3.34) 2) Only the set X is convexified by monomial condensation. 3) The nonconvex objective function of (3.9), not i t s approximation, is directly minimized over the convex approximant of X. The numerical results obtained with Method III are detailed in Chapter VI. As a whole, the results indicate that Method III is preferred. 3.7 Conclusion In this chapter an algorithm has been synthesized to solve signomial programs with equality constraints. The development of the algo-rithm has been shaped by three main concepts: the method of multipliers viewed as a primal-dual algorithm, par t i a l duality, and the retention of the structure of signomial programs for convexification. The algorithm replaces the original problem with a sequence of inequality-constrained signomial programs. However, i f each subproblem is written out e x p l i c i t l y in the standard format of a signomial program, inconvenience in data preparation and computational d i f f i c u l t y due to large problem size are both very l i k e l y . To circumvent these problems, three indirect methods based on monomial condensation have been numerically explored, and one has been singled out as computationally promising. These methods require no or limited increase in variables and constraints, and the effort for data entry i s the necessary minimum. Finally, the convergence of the proposed algorithm is related to that of the method of multipliers and the A v r i e l -Williams algorithm. The proposed algorithm, as stated in Section 3.4.3, serves as a broad framework for specific methods. The conditions for the algorithm's convergence permit wide latitude of arbitrariness in casting the algorithm into a more concrete form for implementation in a computer. One major arbitrary area, the algorithm for solving the subproblem (3.9), has been discussed and w i l l be continued in Chapter V. The other important un-specified step in the proposed algorithm is the updating rules of X. and K. The discussion of these rules i s deferred to Chapter V. In both parts, the selected procedures not only must satisfy the basic conditions for convergence, but also should exhibit experimental efficiency in solving problems. In Section 3.6, the computational performance of the three indirect methods has been assessed on the basis of numerical experience acquired in the process of selecting a suitable method. In the subsequent chapters, further computational considerations are discussed and more numerical experimental results are offered. IV. A PROPOSED ALGORITHM FOR INEQUALITY-CONSTRAINED SIGNOMIAL PROGRAMS 4.1 Introduction The Method II discussed in Section 3.6 for solving the subproblem (3.9) i s effectively the application of the Avriel-Williams algorithm to transform a nonconvex signomial program to a sequence of convex geometric programs via posynomial condensation. Each of the convex geometric programs, in turn, is solved by Dembo's [12], [16] cutting plane algorithm. In this method, the signomial objective function is treated as a constraint through the use of an additional variable. The computational experience associated with Method II indicated that convergence is slow, that many cuts may be re-quired, and that most of required cuts are to satisfy the constraint derived from the objective function. A suspected reason why the objective function turned constraint requires so many linearizations i s the two-stage approxi-mation of i t s many terms with a single monomial term. To overcome this probable block on achieving acceptable convergence, an algorithm i s proposed in which only the nonconvex feasible set of an inequality-constrained signomial program is convexified and then outer linearized. The nonconvex objective function remains unchanged in the whole iterative process. As in any approximation-based strategy, the approximating problem is useful only i f i t can be handled with ease. In the problem of interest, an effective algorithm for minimizing a nonconvex function subject to linear con-straints must be selected. A promising candidate i n this regard i s the reduced gradient method f i r s t suggested by Wolfe [37]. Its generalized version [38] for nonlinear constraints was ranked as among the best of the techniques tested in Colville's study [39]. The reduced gradient method i s a hybrid of 40. simplex-type algorithms and those that are gradient-based. Thus, i f the solu-tion point i s an interior point, the method zeroes i n towards the solution point with the characteristic rate of the steepest descent method. This feature is a potential aid to accelerate the solution of the approximating convex program either by reaching an interior solution without too many cuts, or by ending at a point at which a deep cut can be made. The algorithm proposed in this chapter is, i n principle, applicable to a l l inequality-constrained signomial programs. For this reason the problem solved by the algorithm is formulated in Section 4.2 as a general signomial program. In Section 4.3, the algorithm i s developed and compared with related algorithms. Important computational considerations are elaborated in Section 4.4. Section 4.5 reports the computational experience with the algorithm. In the last two sections, the inclusion of simple upper-bound constraints and-the extension to non-signomial objective functions are considered. 4.2 Problem Statement The problem to be solved i s (T) min g Q(x) ( 4 > 1 ) S . t . X E F where F = {x: x_ E Rm; 0 < £ x; g i , ^ £ 1> k = 1, 2, ...,p} and g (x) are signomials. Note that F differs from X of (3.1) only i n that there i s no simple upper bound constraint on x. in F. Consideration of such type of constraint i s taken up in Section 4.6. A primal approach i s adopted to solve (4.1). Such an approach i s favorable i f the problem has a high degree of d i f f i c u l t y and/or slack primal 41. constraints. The algorithm to be discussed is particularly useful i f the objective function g 0(x) possesses so many terms that approximating i t with a monomial i s not preferred. 4.3 The Proposed Algorithm 4.3.1 Preview Consider the problem (T) in (4.1). It is a nonconvex program whose nonconvex feasible set F can be convexified by posynomial condensation. In the Avriel-Williams algorithm, the objective function g Q(x) is f i r s t trans-formed into a constraint and added to the set F before convexification. The whole problem is thus replaced with an approximating convex program. In the algorithm proposed in this section, the set F is not augmented with gQ(x.). Hence, only F i s convexified while g Q(x) remains unchanged. The original prob-lem i n this case i s approximated by a subproblem, say (ST), involving the min-imization of the original nonconvex cost over an approximating convex feasible set. The subproblem i s then solved by a combined reduced gradient and cutting plane (CRGCP) algorithm to be detailed later in the chapter. The proposed method has two major steps: the convexification of the set F and the solution of the resultant subproblem (ST) by a CRGCP algorithm. The discussion of how F can be convexified by posynomial condensation may be found in Section 2.6, and is not dealt with here. The reduced gradient method as proposed by Wolfe [37] is stated in the next section for completeness and for explaining subsequent algebraic manipulations and computational considera-tions. 4.3.2 The Reduced Gradient Method Assume that the problem to be solved is 42. min M(x) s.t. Ax = b, x _> 0 (4.2) where A is an m x n matrix and m < n. Let x. be partitioned into basic (depen-dent) variables x^ = (x^ , ..., \ )' and nonbasic (independent) variables —B B m x^ = ( x ^ , .. . , x^ )'. Without any loss of generality, a nonsingular square m x m matrix B can be defined as that consisting of the columns of A associated with Xg> Then = -B _ 1Cx T + B _ 1b (4.3) 0 r i % = ~ ^ % where C are the columns of A after those belonging to B have been deleted. The reduced gradient, i.e., dK(x)/d^, i s dM(x)/dx^ - V M(x) - C^B')" 1 V x M(x) (4.4) The strategy of the reduced gradient method is to decrease M(x) by maneuvering only in the subspace of x^ subject to the simple constraint x^ >_ 0. The equality constraint Ax = b is satisfied by (4.3) and the positivity of x —B is assured by suitable restriction on the step size i n x^. Furthermore, whenever a basic variable vanishes, a simplex pivot operation is done to inter-change the vanishing basic variable with a nonzero nonbasic variable. Algo-rithm 4.1 states formally the reduced gradient method for solving (4.2). Algorithm 4.1 (Wolfe's Reduced Gradient Mehtod) Step 1: Set k = 0; x° satisfies Ax° = b_. k Step 2: Compute VM(x ) and v k =V M(xk) + C Tu k k T -1 k where u = - (B ) V M(x ) ^B k Step 3: Define the descent direction s as i) s k. = 0 i f v =0 and v k > 0 Ni Wi l 43. ' ' P I K K Else, s... = -v. Nx x xx) s B - -B C % Step 4: If ||s^|| < t , e > 0 and small, stop. Else go to Step 5. k k k k Step 5: Find i) A„ - min {-x./s. : s. < 0, i = 1, 2, ..., n} c — M . x i l V k k k k i i ) X A minimizing M(x + Xs_ ) , .X e (0,X ] Step 6: Set X k = X k k+1 k . ,k k x_ = x + X s_ k+1 If X g > 0, set k = k+1 and go to Step 2. Else go to Step 7. Step 7: Perform a pivot operation on B * to interchange the vanishing basic variable with a nonvanishing nonbasic variable. Update the indices of the basic and nonbasic variables. Set k = k+1 and go to Step 2. From the foregoing summary of the reduced gradient method, four important requirements have to be met i f Algorithm 4.1 is to be applied. They are 1) the constraint structure of (4.1) has to be satisfied; 2) an i n i t i a l feasible point has to be available; 3) an i n i t i a l basis matrix and i t s i n -verse has to be determined; 4) the gradient of the objective function can be calculated. The next section gives the necessary algebraic manipulations for casting the outer-linearized version of the subproblem (ST) into the form of (4.1). Consideration of the other requirements i s found in Section 4.4. 4.3.3 Algorithm Development Consider the problem (4.1). Let g^ Cx) = ^ ( i ) ~ » k = 1» 2> p, where P^ and are posynomials. The feasible set F can be rewritten as F = {x : x e Rm; 0 < x L £ x; P (x)/(l + Q k(x)) < 1, k = 1, 2, ..., p} (4.5) 44. Assuming that a point x E F is known, the set F can be approximated with a smaller set F by applying posynomial condensation (see Section 2.2) to the denominator of Qfc(x) at x to yield the monomial Q k(x,x). The set F is then defined as 7 = {x : x e Rm; 0 < x L < x; P k(x)/Q (x, x) < 1, k = 1, 2, ...,p} (4.6) The approximate problem (ST) is (ST) min g (x) s.t. x E F o - (4.7) Both problems (T) and (ST) are nonconvex programs because F and F are gener-all y nonconvex for x j> x L > 0, x E R W . But i f the logarithmic transformation z^= In x^ is used to replace x with £, the problem equivalent to (T) and (ST), respectively, may be stated as follows: (T) min g Q(z) s.t. ze Z (4.8) where Z - {z : z £ Rm; z L < z; P k ( z ) / U + Q k(z)) < 1, k = 1, 2, ..... p}(4.9) (4.10) rm L c — \z_ : z_ £ f and (ST) Z min g 0(z) s.t. z_ E Z where Z = (z_ : z E Rm; z L £ z^ P k (O/Qk(z_, z) <_ 1, k = 1, 2, .. ., p} (4.11) Note that now problem (ST) Z is convex. Another formulation of (ST)^ equiva-lent to (4.10) - (4.11) i s (ST) min g (z) s.t. z £ \ (4.12) where Z L = ^ : £ e R 5 £ 1 G k(z, z) £ 0, k = 1, 2, .. . , p} (4.13) a n d % (z, £) = ^(Py^O/QkCz, z)) (4.14) In (4.12) - (4.14), the following statements are true: 1. G,(jz>J) i s a convex function V _z E Rm, <_ ^ 4 > . 2 . ZT = Z 3. ZT is a convex set. ti In Problem (ST) , a nonconvex function i s minimized over a convex set. To solve this problem, a combined reduced gradient and cutting plane algorithm i s used. The algorithm's strategy requires 1) the outer li n e a r i z -ation of the set Z , 2) the casting of the problem formulation into that of (4.2), 3) the use of Algorithm 4.1, 4) i f required, the improvement of the ap-proximation of Z with cutting planes. To outer linearize Z^ > replace the convex function Gk(z,2) with i t s m first-order Taylor approximation Gk(a,2) + (z - a) * VG^(a_,z) for some a_ e R satisfying a_ _> . Define the resultant polyhedron Z^ C R™ as Z L = (z : z e Rm; z L < z./^Ca.z) + (z-a) * VGj^a.z) £ 0, k=l,...,p} (4.15) Because of the convexity of G, (z ^|) , Z cr ZT . Next let = £ - and £ = [C £ ] be the vector of slack variables. Then define the hyperplane s i sp H c R ^ a s H = {(?', c ') : (?', c ') E R ^ ; 0 < (?', c ') ; [J : I](C\ ? ') - b}(4.16) — —S — — S — ~ — 8 S where J i s the p x m Jacobian matrix [ 3 G , (a,j|)/9z. ], I i s the pth-order ident-i t y matrix, and b_ is the p x 1 vector - [ ( z L - a.)' VGk(a,^) + Gk(a,2) J. The problem suitable for solution by Algorithm 4.1 i s (1ST) min gQ(0 (4.17) s.t. (?, ? ) e H — —s * * Let the solution to (LST) be (or correspondingly, some z_ ) . If * — z_ s Z- J proceed to the next convex approximation of Z. Otherwise, add a cutting plane to Z^. The cutting plane i s obtained by linearizing the most — * violated constraint of Z at jz . Li 46. 4.3.4 The Algorithm In this section the proposed algorithm for solving inequality-con-strained signomial programs of the type (4.1) is formally stated as Algorithm 4.2. Algorithm 4.2 (k) Iteration 0: Set k = 1. Assume that a point z e Z is known (k) Iteration k:; Step 1: Condense at 2 to produce problem (ST) Z Step 2: Set i = 1. Set a ( k ) = z^ k ) — (k) Step 3: Linearize at a . and set up problem (LST). Step 4: Solve (LST) by Algorithm 4.1. Let the solution (in the z - domain) be z ^ . Step 5: If z ^ e Z go to Step 7. Otherwise determine the most violated constraint, linearize i t at z , and use the linearization as the next cutting plane. Proceed to the next step. Step 6: Determine a new basis and a new feasible point for prob-lem (LST). Set i = i+1 and go to Step 4. Step 7: If z ^ satisfies the termination criterion, stop. Else, set k = k+1, z ( k ) = z ( l ) , and go to Step 1. Step 6 has to be included because before Algorithm 4.1 can be applied a feasible point and a basis must be available. How Step 6 is carried out is among the practical computational considerations taken up i n Section 4.4 4.3.5 The Algorithm in Perspective Algorithm 4.2 has been developed to solve signomial programs character-ized by a high degree of d i f f i c u l t y and by an objective function with many terms. The main basis of the algorithm is that a set defined by signomial inequalities 47. can be readily approximated by a convex set via posynomial condensation. In this regard, the algorithm is similar to the Avriel-Williams algorithm. There i s , however, a major difference between the two. Unlike the Avriel-Williams algorithm, the approach taken here does not incorporate the objective func-tion into the constraint set. It is hoped that by minimizing the original generally nonconvex objective function over the approximating convex set, convergence towards the original solution can be attained faster. The decision to approximate only the feasible set and leave the objective function untouched means that an efficient method for minimizing a nonconvex function over a convex set must be employed. A promising method i s the novel algorithm proposed in this chapter. The algorithm merges the a b i l i t y of Kelly's [40] cutting plane method to handle convex constraints with the efficiency of the reduced gradient method to minimize a general nonlinear func-tion subject to linear constraints. To tackle the issue of f e a s i b i l i t y , the algorithm exploits the con-vexity property by adopting the strategy of relaxing the constraints by outer linearization and progressively tightening the constraints with cutting planes. The reduction of the nonconvex cost over the set defined by the relaxed con-straint (a convex polyhedron) is achieved by the reduced gradient method. Thus the algorithm translates into a sequence of programs each solved by the reduced gradient method. In contrast, Kelley's cutting plane method for convex programs requires the solution of a sequence of linear programs. The l i k e l y benefit of Algorithm 4.2 is that because maneuvers into the interior of the feasible set are allowed, an interior solution point can be reached sooner, or i f the s o l -ution i s a boundary point, f e a s i b i l i t y can be achieved with deeper cuts. Another interesting comparison that can be made here is with Abadie and Carpentier's [38] generalized reduced gradient (GRG) algorithm for handling 48. nonlinear constraints. Because GRG is designed for general nonlinear constraints, f e a s i b i l i t y i s maintained in each iteration by a restoration phase i n which a system of nonlinear equations (the constraints) is solved by Newton's method. In the problem of interest in this chapter, because of the convexity of the feasible set, the need to maintain f e a s i b i l i t y in each iteration i s dispensed with and i s replaced with a relaxation strategy [41]. In GRG, the inverse of the Jacobian of the constraints with respect to the basic variables has to be either exactly or approximately updated at each point or when the basic variables are changed. This step can become a heavy computational load. In the suggested combined approach, because the relaxed constraints are just linear equations, the basis' inverse i s constant and i s easily updated by simple pivot operations when the basis i s changed. 4.4 Considerations for Implementation Together, Sections 4.3.2 and 4.3.4 spell out clearly the major steps of the proposed technique for solving signomial programs. However, i f the proposed method i s to be implemented i n a computer, some important practical questions need to be answered. How can the point c Z be found? How can the i n i t i a l basis and the i n i t i a l feasible point be derived before a c a l l to the reduced gradient method (Algorithm 4.1)? How can the basis' inverse be updated? How should the one-dimensional minimization required in Step 5 of Algorithm 4.1 be carried out? What are the termination c r i t e r i a and toler-ances? These are the questions to which this section is addressed. 4.4.1 Finding a feasible Point of a Signomial Program A requisite for approximating the set F defined in (4.1) with a smaller disguised convex set F is a point x e F. To find x, Dembo's [16] Phase I method i s used. In this method, a sequence of programs each of the 49. form p (PHSI) min ^ o>k s . t . P k ( x ) / ( 1 + Q k(x)) < (4.18) 0 < x J J _ < x , l _ < t o k , k = l , . . . , p i s solved. Each problem (PHSI) i s solved by Avriel-William's algorithm in. conjunction with Kelley's c u t t i n g plane method. (See [12], [16]). The se-quence i s terminated e i t h e r i f w, = 1, V" k, or some w, - ^ 1 when the optimal s o l u t i o n to (PHSI) has been obtained. 4.4.2 Getting an I n i t i a l Basic Feasible Solution of Prolbem (LST) - No Cut Added From (4.16), the problem (LST) i s defined as min g D(£) t s . t . (a) [ J : I ] [—]•=' b s (4.19) (b) 0 < c, 0 < c — --s where j; = z - z", J = [ 9G k(a,£) /dz ± ] , b = - [ ( z L - a)' VGfc (a ,J) (a, z) ], z e Z, and a^ e R m s a t i s f y i n g a _ 2 l ^ * A convenient i n i t i a l f e a s i b l e s o l u t i o n to (4.19) can be obtained by the following lemma. Lemma 4.1 Assume that jz e Z and l e t £ = _z - . I f a = z and = b - j£_ f then the point (|_, £^) i s a f e a s i b l e s o l u t i o n of problem (LST) . Proof: By the hypothesis and Lemma 2.1, Gk(!>D = 1" (Pk(A)/Qk(i»A)> = l n ( P k ( l ) / ( 1 + Qk(D) ) . n (4.20) But _z e Z. I t follows that ln(P k(£)/(! + Qk(z)))< 0 k^ (4.21) Or for some £ > 0, —s — G k(z, £) + = 0 • k^ (4.22) —s 50. Adding to both sides (4.22) the quantity j£ leads to the desired conclusion. Q.E.D. While the point ) is feasible for (4.18), i t s nonzero basic v a r i -s ables s t i l l have to be identified. A basis of (4.19a) has to be specified. Or equivalently, a set of m linearly independent columns of the composite mx(m+p) matrix [J:I] has to be known. A procedure to identify a set of such columns is explained in the next paragraph. Let the index sets V and W be defined as (4.23) W= {k^gk > °) = a32' UL2* (4.24) Clearly, L^ + L2 = P> Reorder the rows of (4.19a) such that the f i r s t L^ rows are indexed by V and the last L 2 rows, by W. Let the reordered equation with V-'{k:C 8 k-'0> = {v 2, v 2, v L i> C =? and C = ?, be written as s s Q. s ^ 1 i 0 0 i ^2 = b (4.25) It is assumed that the components of £ have been appropriately recorded. Then a basic solution to (4.25) consists of L^ variables selected from the nonzero elements of and the last L 2 nonzero slack variables. This basic solution i s feasible by Lemma 4.1. The L^ variables from must correspond to l i n -early independent columns selected from [-~] or simply from Q. An algo-us rithm such as Gaussian elimination used for the determination of a matrix's rank can be used to pick the desired L^ columns. Let the L^ linearly independent columns of [-—] be [-§-] where 0 i s a L x L n matrix and S i s s S ' 1 i a L£ x L j matrix. Then the i n i t i a l basis B of (4.25) i s I I Q B = Q s L 2 (4.26) 51. Since Q i s nonsingular by construction, the i n i t i a l inverse of the basis i s 1 _ 5 - i 0 ~ ~_1 1 (4.27) - S Q | Hence, in general, only a smaller matrix 6 need to be inverted to obtain the f u l l basis inverse. To summarize, the following steps are undertaken to obtain an i n i -t i a l basic feasible solution to the problem (LST) with no cut added: 1. Obtain a point z e Z and set a = z, ? = z^ - _zL . 2. Generate the matrix J and the vector b. Set £ = Jc, - b . — ~s — — 3. Reorder J, b_ and C such that the f i r s t L 1 rows are indexed by V while the last rows are indexed by W. 4. Identify L.^ linearly independent columns in Q. Let these columns form the L^ x L.^ matrix Q . 5. Invert Q and generate B ^ according to (4.27). 4.4.3 Getting an I n i t i a l Basic Feasible Solution of Problem (LST) - After the Addition of a Cut At the end of one application of the reduced gradient algorithm (i.e. the completion of Step 4 of Algorithm 4.2), l e t the solution be (i,£ ). —s Suppose that at this solution point, at least one of the constraints of the convex set Z^ is violated. Let the most violated constraint's index be r. Then with z = £ + z^, G (z,z) > 0. The cut to be added to the polyhedron ZL i s the supporting plane of Gr(z,J) at z_ = z. Expressed in terms of the equation of the new cut is 6 ' ? + Y (4.28) 52. where = !G^(.z,z) , c = -\{z - z) ' VGr(z,-2) + G r(z , 2 ) ] and y is the new slack variable associated with the new cut. Equation (4.25) becomes • - r ~ S i 0 - r -l* i 0 L2 (4.29) Let y = c - j? Clearly, (_£, ? , y) is a basic solution to (4.29). However, i t i s not feasible since by the assumption that Gr(z,_2) > 0, y < 0. A method to drive Y into the feasible set defined by (4.29) and the positivity constraint on (£, Z ,y) i s to solve the following auxiliary problem (AUX) (AUX) min -y (4.30) s.t. (4.29) , ? > 0, c > 0 —s — As i t is formally stated, the problem (AUX) is a linear program i n a form ame-nable to the use of a composite algorithm [42] to obtain an i n i t i a l basic feasible solution. But because the reduced gradient method allows nonzero nonbasic variables, the point X ? , , ^ , ? ) cannot be used as the i n i t i a l basic feasible point for the LP solution of (AUX). Using a completely different starting point would not only require a Phase I operation but also would des-troy the continuity of the original iterative process. A very convenient answer to this d i f f i c u l t y i s to solve problem (AUX) by a slightly modified ver-sion of the reduced gradient method. The only modifications are the objective function and the definition of the bound on the step size variable A . The usual bound i s defined as A . , = min {-x. /s. : s. < 0, Vs i} (4.31) For solving problem (AUX), (4.31) is changed to the following 53. X M = min [min {-x±/s± : x± > 0, s± < 0 , ^ i > , {-Y/S^: y < 0, > 0 } ] ( 4 . 3 2 ) where is the component of the direction vector s_ along the new dimension of y. The criterion specified by ( 4 . 3 2 ) ensures no new i n f e a s i b i l i t y , and elimin-ates the i n f e a s i b i l i t y due to y < 0 as soon as i t is possible. If both A and s^ are negative, then the constraints are inconsistent. Once y >_ 0, the bound given by ( 4 . 3 1 ) is used, and the original nonconvex objective function g Q(£) i s i n force. The i n i t i a l basic feasible solution to the problem (AUX) is simply (_?,£,y). The starting basis' inverse is ,-1 B -1 -1 B V B-l 0 1 ( 4 . 3 3 ) where B " i s the basis' inverse at the termination of Step 4 in Algorithm 4.2, and w' is the vector of the components of 6' that are associated with the basic variables in Note that once y > 0, B serves as the starting inverse for the next execution of Step 4 in Algorithm 4.2. Thus, from a computational viewpoint, using the modified reduced gradient method does not perturb too much the sequence of solution points obtained by Algorithm 4.2 . Furthermore, the additional programming effort required to implement the method i s very small since the method is s t i l l within the algorithmic framework of the main reduced gradient method. 4.4.4 A Linear Search Method Like other gradient-based algorithms, the reduced gradient method requires in each of i t s iterations the solution of the following one-dimensional problem: f(A k) = min f(A) = min g o U k + As k) s.t. A (0, A M] ( 4 . 3 4 ) 54. The conceptual algorithm calls for an exact linear search. In practice, how-ever, exact linear searches often consume an exorbitant part of the total com-puting time of a problem. Generally, approximate linear minimization i s found adequate to achieve convergence at a reasonable rate. The inexact approach i s adopted in this section. Inexact linear search methods generally f a l l under the following categories: direct searches involving function values only, low-order curve f i t t i n g requiring both function and/or gradient values, and nonoptimal step methods. The f i r s t two types are effectively iterative schemes which produce a sequence converging to the true minimizer of g(A). The nonoptimal step method subscribes to a different viewpoint, namely, that linear search i s a step imbedded in a larger algorithm. The convergence of the larger algorithm can be, and i s , maintained as long as the cost function is decreased s u f f i c -iently in each linear search. Nonoptimal step methods are preferred because they are operationally well defined and are e f f i c i e n t . Examples of nonoptimal step methods are the step size rules of Goldstein [43] and Armijo [44]. Other strategies were outlined by Bard [45]. In the problem (LST) that i s of interest, the signomial objective function gQ(.C) is continuous and d i f f e r e n t i a b l e over H. Furthermore, the evaluation of the derivative at a t r i a l point requires very l i t t l e additional . effort once the function has been evaluated. This observation suggests that a low-order polynomial f i t , such as the cubic f i t scheme derived by Davidon [46] is suitable. However, a curve f i t t i n g scheme by i t s e l f may s t i l l require numerous repetitions i f the termination criterion i s the usual relative change k k in f ( X ) or i n X . An appealing approach is to couple Armijo's test for s u f f i c -ient descent to the cubic f i t method. The strategy is to use Armijo's step size 55. rule except when i t i s clear that a relative minimum exists in the interval (0, X M ] . This condition i s confirmed, when df/dX < 0 at X = 0 and df/dX > 0 at X^. The proposed linear search method, stated as Algorithm 4.3, is thus a synthesis of the cubic f i t method and Armijo's step size rule. Algorithm 4.3 Step 0: Given a e (0,1), g e (0,1), y e (0,1), and X q > 0. Let X^ be determined by Step 5 o f Algorithm 4.1. Set X^ = min [X 0,X M] and X^ = 0. Calculate f(X ) and f';(X ). Step 1: Set X = X^, and evaluate f(X) and f'(X). A Step 2: If f*(X) > 0, go to Step 4. Else go to Step 3. Step 3: Set X = X i f f(X ) < f(X.) and go to Step 5. Else, X = YX and go m m x, r m m to Step 1. Step 4: Apply Algorithm 4.4 to f i t a cubic through (X„, f(X.)) and (X , f(X ) ) . x x m m Let the minimum estimate be X. Step 5: If f(X) - f U ^ ) < -a|x - xj |f\(X )|, Xfc = X stop. Else continue. Step 6: If X = X , go to Step 7. Else, set X = X i f f'(X) < 0, or X = X m Z m i f f'.(X) > 0. Go to Step 4. Step 7: X = g A . Go to Step 1. c — m m ^ Algorithm 4.4 (Cubic Fit) Step 0: Given X , A , f(X ), f(X ), f'(X.) and f'(X ). Assume f'(X„) < 0 and f'(X ) > 0. m — Step 1: Calculate u = f (A ) + f (X ) - 3(f(X.) - f(X )) / (X - X ) J- x m x m x m u 2 = [ u 2 - f ' ( X , ) f ' ( X m ) ] 1 / 2 Step_2: X = X m - (X^ - ( f ' ( X J + - u j / (f' ( X J - f ' (X £) + 2u£) 56. Algorithm 4.3 combines the useful features of a cubic f i t and Armijo's step size rule. The rule assures convergence by stipulating sufficient descent before a t r i a l point i s accepted as the next solution point. But the ith t r i a l point i s not always obtained by the relation = P** m^ »• a s x t i s the case in Armijo's algorithm. When the situation warrants i t , the a b i l i t y of a cubic f i t to rapidly pick a good interior estimate is exploited to produce a t r i a l point. 4'.4.5 Optimality Cri t e r i a , Feasibility Tolerance, and Updating the Basis' Inverse • . The main problem (T) posed i n (4.1) i s a nonconvex program. Solving (T) numerically with an algorithm generally leads to at best a local minimum at which the Kuhn-Tucker optimality conditions must be satisfied. In practice, unless a Lagrangian-type algorithm is used, testing the first-order necessary optimality conditions i s l i k e l y not possible since the optimal Lagrange multi-pliers may not be known. Even i f the multipliers are known, checking the Kuhn-Tucker conditions at the end of each iteration may create too large a compu-tational load. The whole question of knowing when optimality i s reached be-comes a trade-off between r e l i a b i l i t y and extravagance. For the sake of com-putational efficiency, simple but imperfect stopping rules for terminating the iterative process have to be adopted. The rule adopted in this chapter i s : Terminate i f the relative change in g Q (C.) £ £ and the relative change in £ <. e, where e is a small preassigned positive number. If either g Q or £ -*• 0, replace relative change with absolute change. Depending on which quantity is of i n -terest, the rule may be relaxed in two ways. If only g Q is important (e.g. only cost, not the design parameters, is sought), the test on the relative change in _£ may be dropped. Analogously, i f only the accuracy of i s emphasized, the test on g Q may be avoided. An excellent operational check on the answer obtained 57. by the previous stopping rules i s to perturb the s o l u t i o n and resolve the problem with the perturbed point as the s t a r t i n g point. Another important issue that must be faced i n implementing a con-ceptual algorithm i s the d e f i n i t i o n of f e a s i b i l i t y , or more generally, the d e f i n i t i o n of i n e q u a l i t y and e q u a l i t y . T h e o r e t i c a l l y , whether a and b are s c a l a r s or vectors, the expression a <^ b, a _> b, and a = b are unambiguous. However, i n the implemented ve r s i o n of an algorithm, because numbers are repre-sented i n f i n i t e lengths, i n t e r p r e t i n g i n e q u a l i t i e s and e q u a l i t i e s i n t h e i r s t r i c t mathematical sense could j u s t spawn a myriad of nuisances such as un-expected e a r l y program termination, undesired negative q u a n t i t i e s , d i v i s i o n by zero and other d i f f i c u l t i e s . For t h i s reason, i n the implementation of t h i s t h e s i s ' algorithms, i n e q u a l i t i e s and e q u a l i t i e s are replaced with s - i n e q u a l i -t i e s and e - e q u a l i t i e s . This means that given s c a l a r s a and b, and a small £ > 0 , a £ b + e replaces a <_ b; a>_b - E, a>_b; and |a - b| <_ e replaces a = b. If a. and b_ are vectors, then a__<b_ + e ^ a ± — ^± + e " S i m i l a r l y , a_ > b_ - e ^ a^ >_ b^- £ , and | | a_ - b_| | <_ t =^ a_ = b_. Note that the Euclidean norm i s used. The f i n a l point considered i n t h i s s e c t i o n i s the updating of the b a s i s ' inverse B ^ = [b. ]. The procedure followed i s very s i m i l a r to that of X j the revised simplex method. Let r be the index of the leaving basic v a r i a b l e . The key steps are: (1) Determine a nonzero nonbasic v a r i a b l e , say, Let i t s associated column be a . (2) Calculate v = B "*"a . I f y =0, return to (1) to determine another v a r i a b l e . If y^^ = 0 for a l l s, the system i s i l l -defined. (3) Assume that y r g 4 0. Let B = [b_^ _. ] be the new inverse. Then b. . = b. . - (y. b .)/y , i 4 r (4.35) b . = b . /y , i = r r j r j ' r s ' '58. 4.5 Computational Experience 4.5.1 Posynomial Representation A posynomial P(x) i s usu a l l y written as (4.36) orthant x > 0 (x-domain) to the r e a l space R (z-domain) by the l o g a r i t h -mic transformation z. = In x.. As a function of z_, P can be w r i t t e n as The r i g h t side of (4.37) i s defined as the exponential form of P(x). While the expression i n (4.36) i s generally nonconvex, the function P(z_) i s always a convex function of z_. The hidden convexity of (4.36) over the whole r e a l space i s thus unveiled by the logarithmic transformation. Further discussion on convex-transformable nonconvex functions may be found i n [47]. t h i s chapter's primal approach to s o l v i n g signomial programs, the expon-e n t i a l form i s favored for two important p r a c t i c a l reasons: 1) The prob-lem (LST) solved by the proposed combined reduced gradient-cutting plane algorithm i s defined i n the z-domain. 2) The CPU time required by the algebraic operations i n the exponential form i s much l e s s than that of the ordinary form. The second reason can be e a s i l y explained by consid-e r i n g the computational e f f o r t required to compute (4.36) and (4.37). The f i r s t equation requires mn m u l t i p l i c a t i o n s , nm exponentiations, and nm-1 additions , w h i l e the second requires nm-f-n m u l t i p l i c a t i o n s , n ex-ponentiations, and nm-1 additions. A comparison between the two shows J n m p ( z ) = 1 c. exp( I a z ), z e R.V^j i = l j = l (4.37) The exponential form has more than a formal s i g n i f i c a n c e . In 59. that the exponential form requires n multiplications and n(m-l) additions more, but n(m-l) exponentiations less than the ordinary form. Because an exponentiation usually requires much more time than an addition or a multiplication \ using the exponential form in any repetitive calcula-tion of posynomials should definitely be preferred. It i s true that i f the i n i t i a l starting point and the f i n a l solution are given i n the x-domain, there i s the additional overhead of taking m logarithms and m exponentiations. This overhead, however, becomes insignificant when compared to the savings in the iterations. In an attempt to gain some idea on how much saving can be obtained, Dembo's [16] program was com-piled in two versions and then used to solve a four-variable, two-con-straint problem. The improvement in total CPU time of using (4.37) i n l i e u of (4.36) was about 26%. 4.5.2 Software Implementation In this section the software implementation in FORTRAN IV of the proposed CRGCP algorithm is b r i e f l y described. The complete software package, hereafter referenced as SP, ' consists of a master program MAIN which oversees the operations of 16 subroutines. As a whole, the computer program SP has the following highlighting features: 1. It contains a l l the capabilities of Dembo's program GGP [16]. 2. It i s compatible with the I/O format of GGP except for minor changes i n the control cards and in the order of reading the data of the objective function. For example, in the IBM 370/168, the CPU time required by an expon-entiation operation i s approximately 17 and 26 times more than that of multiplication and addition, respectively. 60. 3. It uses the exponential form (4.37) for computing posynomials. 4. It offe rs a choice of whether GGP or CRGCP should be used. The choice i s indicated by the value of a single control parameter. 5. If CRGCP i s selected, a function OBJFCN declared in an EXTERNAL statement has to be supplied for the purpose of computing the function and gradient values of the objective function. The last two features provide the f l e x i b i l i t y of using CRGCP to solve signomial programs whose objective function is more complex than an ex p l i c i t l y stated signomial. This aspect w i l l be thoroughly discussed i n Section 4.7. It w i l l also prove useful for solving the subproblem (3.9) of Algorithm 3.2 designed to solve signomial programs with equality constraints. 4.5.3 Computational Results In order to acquire some experience with the performance of the CRGCP algorithm, the code SP was used to solve a set of signomial programs. The test problems together with the starting points are given below. Note that a l l the starting points are infeasible. Problem A (A quadratically constrained quadratic program) minimize g n(x) = 3x 2 + 2x? + x 2 + x 2 + 7x7 + 565 - 39x„ - 17x_ - x. 0 — 1 2 3 4 1 2 3 4 s . t . g;L(x) = (l/6)x 2 + (l/18)x 2 + (l/9)x 2 + (1/I8)x1 - (l/9)x 2 - (l/18)x 2 < 1 ? 2 2 2 g 2(x) = x~ + x^ + 2.5x^ + x^ + 9.5 - 4x 2 - <_ 1 61. 2 2 2 2 g 3(x) = x 2 + 3x^ + x^ + x 2 - x 2 - x^- x^ - 2 <_ 1 0.01 < x <. 10.0, j = 1, 2, 3, 4 Starting Point: x = (2.0, 2.0, 2.0, 2.0)' Problem B minimize g Q(x) = 0.000392xJ + 5.46x^ + O.OSx'^x" 1 + 6 . 0x 2x 3 + 2.5856X"1 - 0.125x2 - 1.8x l X 3 s. t. g 1(x) = 0.03x2 + O.lx 2 - 0.15x2 <_ 1 g 2(x) = 5.0x^x2 + 3.6x 2'\ 3 - 0.04x 2x 3 1 <_ 1 0.01 £ X j £ 10.0, j = 1, 2, 3 Starting Point: x = (2.0, 2 .0 , 2.0)' Problem C (Heat Exchanger Design [48]) minimize g Q(x) = + c 2 x 2 + c ^ s. t. g l ( x ) = c 4 x x x ^ " 1 + c ^ " 1 + c ^ X g 1 < 1 •r \ _ 1 _ 1 -1 ' -1 -1 g 2(x) = c 7 x 2 x 5x^ + c gx^x 7 + C g X 2 x^x 7 £ 1 83(2l) = ^ o^l1 + 'uXs**1 + C 1 2 X 3 l x 5 x 8 1 - 1 84(*> = c 1 3 x 4 + c 1 4 x 6 < 1 g 5(x) = c 1 5 x 5 + c 1 6 x ? + c 1 7 x A < 1 86(2L) = c 1 8 x g + c 1 9 x 5 < 1 100.0 <_• x 1 <_ 10000.0 1000.0 < x. < 10000.0, j = 2, 3 10.0 <_ x <_ 1000.0, j = 4, 8 62. Starting Point: x = (5000.0, 5000.0, 5000.0, 200.0, 350.0, 150.0, 225.0, 425.0) The coefficients c^, are given in Table 4.1. The s t a t i s t i c s of the test problems are summarized in Table 4.2. In this table the degree of d i f f i c u l t y D i s defined as the difference be-tween the total number of terms and the number of variables plus one. The terms due to the a r t i f i c i a l simple lower and upper bounds are neglected. If a dual approach were used to solve the problems, then the dimension of the respective dual problem would be D. A l l the problems were solved on an IBM 360/168 computer using -12 double-precision arithmetic. The sign tolerance was 10 while the feas-i b i l i t y tolerance was 10 The linear search was inexact and was term-inated either when the norm of the directional derivative along the search -3 direction was less than 10 or when a specified number of searches N s had been completed. The f i n a l termination criterion was relaxed to be -4 that the relative change in g Q(x) be less than 10 . The solutions obtain-ed by SP are given i n Tables 4.3 (a)-(b) for N = 5 and N = 10. Also s s quoted in the tables are the elapsed GPU time for obtaining each solution. The elapsed time i s the time between the moment when a feasible point has been obtained and the time when the termination criterion i s satis-fied. For the purpose of comparison, the solutions obtained by Dembo's G G P code under the same f e a s i b i l i t y and optimality tolerances are includ-ed in the last column of each table. 63. j C j 1 1.0 11 1.0 2 1.0 12 -2500.0 3 1.0 13 0.0025 4 833.33252 14 0.0025 5 100.0 15 0.0025 6 -83333.333 16 0.0025 7 1250.0 17 -0.0025 8 1.0 18 0.01 9 -1250.0 19 -0.01 10 1250000.0 Table 4.1 Coefficients for Problem C A B C No. of Variables 4 3 8 No. of Constraints 3 2 6 No. of Terms 30 13 19 Degree of D i f f i c u l t y 25 9 10 Table 4.2 Statistics of the test problems for the CRGCP algorithm. 64. Starting SP SP GGP Point N =5 s N =10 s x l 2.0 2.05478 2.05510 2.04978 X2 2.0 2.45474 2.45420 2.46178 x3 2.0 0.48977 0.59120 0.57044 X4 2.0 0.44928 0.44892 0.44517 Optimal Cost - 498.44053 498.44000 498.45041 CPU Time 1.190 sec (a) 0.933 sec 0.531 sec Starting SP SP GGP Point N = 5 s N =10 s X l 2.0 2.92012 2.92008 2.92042 X2 2.0 2.72930 2.72931 2.72920 X3 2.0 0.06929 0.06929 0.06930 Optimal — 0.40205 0.40205 0.40205 Cost CPU Time - 0.264 sec 0.109 sec 0.307 sec (b) Table 4.3 Starting point and solution of (a) Problem A, (b) Problem B 65. The numerical results of Table 4.3 indicate that while the code SP i s not consistently faster than GGP, the rate of convergence i s i n the same order of magnitude as that of GGP. The solutions obtained by the two codes are reasonably close. On a percentage basis, the variation in the optimal cost i s less than that in the norm of x. The smaller difference in the optimal cost i s perhaps due to the bias introduced by the choice of the termination criterion. Stopping the inexact linear search earl i e r by reducing N may hasten convergence. However, in certain problems the s gain i n time due to less linear searches may be offset by the need to solve more approximating geometric programs. Starting SP SP GGP Point N = 5 s N = 10 s x l 5.D3 553.58153 569.28585 568.55211 X2 5.D3 1367.87035 1369.15269 1372.33969 X3 5.D3 5128.21662 5110.96764 5108.52350 X4 200. 179.82943 181.17369 181.11155 x 5 350. 294.87134 295.56165 295.65912 X6 150. 220.17088 218.82631 218.88862 x 7 225. 284.95809 285.61097 285.45120 x8 425. 394.87134 395.56161 395.65911 Optimal Cost 7049.6685 7049.4062 7049.4107 CPU Time - 1.025 sec 1.268 sec 1.501 sec (c) Table 4.3 (cont'd) Starting point and solution of (c) Problem C 66. 4.6 Presence of Simple Upper Bounds Simple upper bound constraints of the form x < x U frequently occur in engineering design problems. In the problem statement of (4.1), such type of constraint is not e x p l i c i t l y expressed since in principle u x £ x can be included i n the feasible set F by writing m monomial con-straints of the form X j / X j £ !• Such a straightforward approach, however, i s computationally very i n e f f i c i e n t . In this section the problem formu-lation of (4.1) i s extended to allow explici t simple upper bounds on JC. The special simple structure of these constraints is exploited to specify modifications in Algorithm 4.1 so that the upper bounds are indirectly taken into account without any increase in the number of variables or constraints. The problem to be solved i s min g Q(z) (4.37) s.t. x e X where X = {x : x e Rm; 0 < x~" £ x £ x U; £ x> k = 1, 2, . . . p} and g (jx) are signomials. As in Section 4.3.3, i f x. = exp(z.), (4.37) 3 3 can be expressed in terms of z_ to yield the following problem min g 0(z) (4.38) s.t. z_ e Z where Z = {z_ : z_ e Rm; £ L £ z £ z_U; p k(^)/(1+Q k(z)) £ 1, k = 1» 2, ... p}. It follows that (4.38) can be solved by Algorithm 4.2 of Section 4.3.4 except now the problem corresponding to problem (LST) of (4.17) i s defined as min g0(jL) 8- t"- [J J I] [---] - b ( 4 ' 3 9 ) ~s 0 ± JL £ LU, 0 £ ^ 67. where J and b are as defined in (4.16), and = ln(Xj/x_. L) , j = 1, 2, m. Thus, the presence of simple upper bounds on x leads to simple upper bounds on but not on the slack variables j ^ . It follows that simple upper bounds in the original problem necessitate changes only in the re-duced gradient method used to solve (4.39). Consider the reduced gradient method summarized in Algorithm 4.1. For convenience of notation, define the upper bound on c; as a s vector ? U of a r b i t r a r i l y large numbers. Then the required changes in Algorithm 4.1 are in the definition of the search direction s^ , the de-k termination of the maximum allowable step size X^, and the criterion for the replacement of a basic variable. The revised algorithm capable of solving (4.39) is lik e Algorithm 4.1 except in the following steps i n the kth iteration: Step 2 i ) s N = 0 i f a) ^ = 0 and > 0 or b) ^ = ^ u and v N < 0 i i i I i i Else, s = -v ; i l Step 5 X = min [min{-S./s.Is. < 0 ,V- j } , mini— 1 1 I s. > 0}]-C — M 3 1 3 1 3 ' . Step 6 £ becomes nonbasic i f ? B = 0 or c,^1 . Note that a basic v a r i -i i i able has to be between the two bounds while a nonbasic variable can be nonzero. For reference purposes, the reduced gradient method with the above rules i s called Algorithm 4.1(a). A numerical example solved by Algorithm 4.1(a) i s Colville's [39] third test problem given below. 2 min g Q(x) = ^ x 3 + c 2 x 1 x 5 + c 3 X j ^ + c^ s.t. g-^x) = x 3 x 5 + c g ^ X j + c ? <_ 1 6 8 . S 2 <£> = c 8 x2 x5 + Cg X l X4 + C10 X3 X5 - 1 g 3(x) = c n x" 1 x' 1 + c 1 2 x x x" 1 + c i 3 x" 1 x 2 x" 1 < 1 84 (£) = c u x 2 x 5 + X i X 2 + c i 6 x 2 < 1 g 5(x) = c 1 ? x" 1 x^ 1 + c l g X ; L x" 1 + c i g x 4 x^ 1 < 1 g 6(x) - c 2 Q x 3 x 5 + c 2 1 x ± x3•+ c 2 2 x 3 x 4 < 1 7 8 1 x 1 i 1 0 2 3 3 1 x 2 - 4 5 27 < X j < 45 , j = 3, 4, 5 The coefficients c^, ... , c 2 2 are given in Table 4.4. A l l the lower and upper bounds on x^, j = 1, 5, are part of the problem form-ulation. Hence the degree of d i f f i c u l t y of the problem i s 26. Two feas-ible starting points and their respective solutions are presented in Table 4.5. The solutions were obtained with a f e a s i b i l i t y tolerance of -6 -4 10 and an optimality tolerance of 10 . Note that the solution i n -volves lower and upper bounds. The CPU elapsed times quoted in Table 4.5 are as defined in Section 4.5.3. 69. j C. 3 j C. 1 5.35785470 12 -0.4200261 2 0.83568910 13 -0.30585975 3 37.239239 14 0.00024186 4 -40792.141 15 0.00010159 5 0.00002584 16 0.00007379 6 -0.00006663 17 2275.132693 7 -0.00000734 18 -0.26680980 8 0.00853007 19 -0.40583930 9 0.00009395 20 0.00029955 10 -0.0087777 21 0.00007992 11 1330.32937 22 0.00012157 Table 4.4 Coefficients of example problem with simple upper, bounds. 1st 2nd Starting Pt. Solution Starting Pt. Solution x l 78.62 78.0 80.0 78.0 X2 33.44 33.0 35.0 33.0 X3 31.07 29.99307 35.0 29.99307 X4 44.18 45.0 35.0 45.0 X5 35.22 36.78135 35.0 36.78135 Optimal Cost - -30670.093 - -30670.093 CPU Time - 0.086 sec - 0.094 sec Table 4.5 Starting points and solutions of example problem with simple upper bounds. 70. 4.7 Minimizing Algebraic Functionals of Signomials 4.7.1 Introduction The algebraic operations are addition, subtraction, multipli-cation, division, and exponentiation. An algebraic functional of signomi-als i s a real-valued functional generated by performing specified algebraic operations on a f i n i t e set of signomials. The addition, subtraction, multiplication, and the positive integral powers of signomials can i n principle be formally expressed as signomials. Ratios and real powers of signomials generally cannot be reduced to the standard form of a s i g -nomial. The problem of interest i n this section is the minimization of algebraic functionals of signomials subject to signomial inequality con-straints. This type of problem is termed as a semi-algebraic program. The reasons for considering this class of problems are: 1. The subproblem (3.9) required in Step 1 of Algorithm 3.2 for sig -nomial programs with equality constraints belongs to this class of problems. 2. A wider range of applications can be modeled by signomial programs. 3. Considerable saving in the effort for data preparation is l i k e l y as the algebraic operations on the signomials need not be carried out. In Section 4.7.1, the problem of interest i s formulated. The same section also discusses how the_ combined reduced gradient - cutting plane algorithm can be used to solve semi-algebraic programs. A numerical example of semi-algebraic programs i s also presented in Section 4.7.1. In Section 4.7.2 the important practical constrained location-allocation problem is shown as an i l l u s t r a t i v e example of semi-algebraic programs. Although some constrained location-allocation problems can be formu-lated as standard signomial programs, such a formulation usually has a 71. large number of terms, hence, a high degree of d i f f i c u l t y . The advantage of solving constrained location-allocation problems by the CRGCP algo-rithm (Algorithm 4.2(a)) i s demonstrated via an example. 4.7.2 Problem Formulation and Solution Let s-^(x) , .... , s^ Cx) D e signomials and define f(s^(x), ... , s^(x)) as a real-valued functional generated by executing a set of speci-fied algebraic operations on s.. s . Then a semi-algebraic pro-l XJ gram may be formally defined as min f (s (x) , . . . , sCx-)) ( 4' 1 - L s . t . x. e X where X = {x : x e Rm; 0 < x L £ x £ x"; S^ Cx) £ 1, k = 1, 2, ... , p} and g^OO > k = 1, 2, ... ,p, are signomials. The only difference between (4.40) and (4.37) is the form of the objective function. But since Algorithm 4.2(a) does not really assume that the objective function has to be a signomial, i t follows that the algorithm remains applicable i f i t is used to solve a semi-algebraic program. This is precisely the approach adopted in this section. In fact, the only necessary restriction of f as a function of x is that f be continuously differentiable over X. To il l u s t r a t e the problem just formulated and the proposed method for solving the problem, consider the following two-variable numerical example min f(x) = ( s 1 ( x ) ) ( s 2 ( x ) ) " 1 - 5 - (s 3(x)) 2'°(s 4(x))°" 5 5 s. t. ( a ) x2/2.25 + x2/16 + 5 - 4x^1.5 - 0.5x2 £ 1 (4.41) (b) x±/5 + x2/15 £ 1 72. (c) 0.947x°\7 + 2.0 - x 2 1. 1 (d) 0.1 £ x 1 100.0, j = 1, 2 where s^, s^, s^, and are signomials given as t v 0.5 0.4 s l ( x ) = x l - x2 s 2(x) = x 2 + x 2 , s 1.25 s 3(x) = f ^ -0.4 . -1.0 s^(x) = x 2 + x i x2 (4.42) Note that the objective function i s a signomial-like func-tional of the signomials s^, S 2 » s^ and s^. The feasible set defined by (4.41) (a) - (d) is a nonconvex set. Example (4.41) was solved by Algo-rithm 4.2 on the IBM 370/168 computer. Beginning with the infeasible i n -i t i a l point (1.5, 1.5) the algorithm obtains via Dembo's Phase I proced-ure the feasible starting point (1.67477, 2.45071). From this feasible point, the algorithm obtains, after four approximations of the noncon-vex feasible set, the f i n a l solution of (2.43246, 7.70263). At the sol-ution, constraints (4.41) (a) and (b) are tight while constraint (4.41) (c) is loose. The overall CPU time required i s 1.01 sec. 4.7.3 Application to Constrained Location-Allocation Problems The optimal location-allocation problem may be generally stated as follows: Given the location and the requirements of a set of fixed known destinations, determine (a) the number of sources, (b) the location and the capacity of each source, (c) the source-to-source and source-to-destination shipment pattern such that the destinations' re-quirements are satisfied and that the total set-up, operation and ship-73. ping costs are minimized. This problem arises in the planning of ware-houses, distribution centers, communication switching centers, production and municipal services f a c i l i t i e s . the Generalized Weber problem, variants of the location-allocation prob-lem have been considered by a host of researchers. A recent bibliography [49] l i s t s over 200 papers in the decade before 1973. Most of the propos-ed methods, however, cannot handle constraints restricting the sources' coordinates. In real situations, some constraints need to be imposed on the location of the sources for a variety of reasons such as geographical barriers, zoning laws, or safety requirements. In recent years some theor-e t i c a l [50] and algorithmic [51] - [53] results related to the constrained location problem have been reported. In this section, a restricted ver-sion of the general problem i s formulated as a semi-algebraic program that can be solved by the approach discussed in this chapter. Algorithm. 4.2 i s demonstrated to be another tool available for solving certain types of constrained location-allocation problems. The version of the constrained location-allocation problem to be solved in this section i s posed as follows: Suppose that the co-ordinates a. and the demand r-of the jth sink, j = 1, 2, ... , m, are 3 J known. Let the i t h source, i = 1, 2, ... , n have a location at x^ and a capacity c^. If i t i s assumed that the shipping cost from one point to another i s proportional to the distance between the points, then the mathematical formulation i s Under the guise of such names as the Fermat, the Steiner, or min m n C(x) = .E. .1, w._,d(x., a.) + - 1=1 j=l j i — j — j n-l n 74. m n-l s.t. ( a) j l i w j i _ j£lvji - c i ' 1 = 1>2>---> n (4.43) n (b) i l ^ - j i l r j > j=l,2,...,m (c) g k ( x 1 5 . . . .x^) <_ 1.0, k=l,2,...,p where d(x, y) i s a distance function of x and y, and g^ Cx^ , ... , x^), k = 1, 2, ... , p are signomials expressed in terms of the sources' co-ordinates. w_.^ i s the amount shipped from source j to sink i while i s the shipment from source j to source k. Note that source n i s assum-ed not to supply anything to other sources in order to model a f a c i l i t y whose output does not serve as the raw material of other sources. In (4.43), i f the variables w.. and v., are a l l held fixed, the problem 3 K becomes a multi-source location problem. The common distance functions used are the Euclidean and the rectangular ("city block") distance func-tions. In the example of this section, the Euclidean distance is used. In (4.43), i f each d is replaced by the appropriate Euclidean distance expression, i t becomes clear that the constrained location -allocation problem as formulated in (4.43) i s a semi-algebraic program. It follows that Algorithm 4.2 can be used in a straightforward manner. This point is il l u s t r a t e d by the solution of the following single-source constrained location problem adapted from [51]. 24 „ • : v x t i \ 2 _ L / i N 2 , l / 2 mm E w^i^ - a][) + (x 2 - a 2) ] i=l s. t. (a) ( X ] L - 44.0) 2 + (x 2 - 36.2) 2 <_ 2.0 (4.44) 75. (b) x 2 + 1.25x < 92.75 (c) 0.0 < x l 5 x 2 The constraints in (4.44) require the source to be located with-in 2 units of (44, 36.2) and to be below the line x 2 = -1.25x '+.92.75. The data of the problem are given in Table 4.6. Originally provided by Kuhn and Kuenne [54], the data are those of the 24 cities in the Ukraine which rank among the 100 most populous c i t i e s in the U.S.S.R. The weight w^ of the i t h city i s taken as the proportion of the city's population of the total population of the 100 largest Russian ci t i e s . The coordin-ates (a*, a 2) of the i t h city are the city's north latitude and east longitude correct to the nearest f u l l degree. The results of solving (4.44) by Algorithm 4.2 are summarized in Table 4.7. In [51], Gurovich did not enforce the two constraints of (4.44) simultaneously. Hence i t is not possible to make a direct compar-ison between her results obtained by algebraic programming (AP) [36], [51] and those reported here. However, the results in Table 4.7 s t i l l are very much superior to Gurovich*s results obtained when only either of the constraints i s enforced. For example, in [51] enforcing (4.44a) only requires 6 convex approximations and about 6 CPU sees, while enforcing (4.44 b) only requires eight convex approximations and about 13 CPU sees. (IBM 370/168). With the use of Algorithm 4.2, enforcing both constraints requires only four convex approximations and 0.874 CPU sec. A plausible explanation for the wide margin in performance i s that in AP the objective function is incorporated into the constraint set. This not only increases the number of variables by one but also requires more approximations. It was precisely this d i f f i c u l t y that led to the development of the main algorithm proposed in this chapter. i i a l i a2 i w 1 54 28 0.012 2 50 24 0.010 3 47 28 0.005 4 . 46 31 0.014 5 47 32 0.005 6 47 33 0.004 7 45 34 0.004 8 45 34 0.004 9 46 . 34 0.001 10 48 35 0.004 11 48 35 0.015 12 48 36 0.010 Table 4.6 Locations and weights I n i t i a l Point Optimal Point Optimal Cost No. of Convex Approx Constraint Tolerance Optimality Tolerance Total CPU Time 76. i i a l i a2 i w 13 47 38 0.007 14 48 37 0.016 15 48 38 0.008 16 48 38 0.007 17 47 39 0.005 18 45 39 0,007 19 47 40 0.014 20 48 40 0.004 21 49 39 0.004 22 41 37 0.021 23 50 31 0.026 24 52 31 0.004 for major Ukraine c i t i e s . (40.0, 20.0) (45.3268, 36.0105) 1.0896 4 0.874 sec Table 4.7 Summary of the solution of the constrained location problem. 77. V. NUMERICAL SOLUTION OF SIGNOMIAL PROGRAMS WITH EQUALITY CONSTRAINTS 5.1 Introduction In Chapter III an algorithm based on the multiplier method and the concept of par t i a l duality has been developed to solve signomial programs with equality constraints. The algorithm i s so structured that a sequence of inequality-constrained signomial programs i s solved. This approach allows the exploitation of a key property of signomial programs, namely, their easy approximation by convex programs. As discussed i n Section 3.5, the convergence of the proposed method i s assured by the convergence of the method of multiplier provided that the assumptions for convergence hold and the following conditions are satisfied: 1) the penalty constant K i s sufficiently large after a f i n i t e number of iterations, 2) a local minimizing point of the primal problem (ir^) i s obtained, and 3) the unconstrained dual problem i s maximized. Except for the discussion on how ( T T J ) could be solved, no specific methods for achieving conditions (1) and (3) are mentioned in Chapter III. The f i r s t purpose of this chapter i s to delineate clearly how each of the three conditions can be met. Secondly, the proposed alter-natives are tested in a numerical experiment and their computational efficiency i s assessed. How the subproblem (TTJ) i s solved i s important not only be-cause the solution method determines the accuracy of the answer wanted but also because i t consumes a major part of the total computing time. As summarized in Section 3.6, earlier computational experience revealed 78. the numerical d i f f i c u l t i e s of the two indirect methods, Method I and Method II, that had been considered for solving ( 1 ^ ) . In Section 5.2, the Algorithm 4.2 developed in the previous chapter is adopted to solve (TTJ). The numerical results presented in Chapter IV suggest that Algori-thm 4.2 can handle well signomial programs whose objective function may have many terms or may not even be pure signomials. This a b i l i t y of the algorithm is confirmed in the numerical tests of this chapter. In Sec-tion 5.3, the discussion focuses on the important question of how to update the multiplier estimate vector X_ and the penalty constant K, given the i n i t i a l values of these parameters. Section 5.4 contains the des-cription and the results of the numerical experiments conducted to test the proposals formulated in Section 5.3. 5.2 Solution of the Subproblem (TTJ) The subproblem (TTJ) of Algorithm 3.2 i s (TTJ) min £(x, X, K) = X j + I X,h, (x) + K £ h 2,(x) (5.1) k=l k=l k s.t. x e X where \ ^ = g p + k ^ ~ 1 ' k = 1 ' 2 ' *** q > and X = .{x : x e Rm; ^ ( x ) < 1, k = 1, 2, ... p; 0 < x L < x < x u} It i s pointed out in Section 3.6 that the signomial £(x, X_, K) i s l i k e l y to have many terms. This observation suggests that any method which approximates £(x, X_, K) grossly may f a i l to exhibit good convergence, i f there is any at a l l . Algorithm 4.2 developed in the last chapter i s designed to solve inequality-constrained signomial programs without approximating the objective function. Later in Section 4.7, the algori-thm i s extended to solve programs whose objective functions are algebraic 79. functionals of signomials. Thus, Algorithm 4.2 i s just suitable for (TT^) since a quick examination of (5.1) shows that £(x, K) i s a qua-dratic form of the vector [hi, b.2, ... h^]'. I n t n e numerical tests reported in this chapter, the solution of the subproblem ( I T ^ ) i s accom-plished with Algorithm 4.2. This i s the Method III alluded to at the end of Chapter III. The application of Algorithm 4.2 requires the calculation of .34 •3x, £both the function &(x, X, K) and i t s gradient V &(x, _X, K) = [ — , 3£ - r — ] ' . The additional effort to calculate V £ i s l i t t l e indeed. Let 3 xm x m a i j k m b d JJ x i j k . [k] and {k} are appropriately defined index sets. Then 3fc i=i 1 for j > 2 where 3 1 8 X J je[k] i 3 k 3 k je{k} i j k j k ' Since in the process of calculating h^, uj k» je[k], and v ^ t je{k} can be temporarily stored, calculating ^ nk simply needs the additional arithmetic operations specified by (5.3). 5.3 Schemes for Updating _A and K 5.3.1 Introduction Algorithm 3.2 developed to solve signomial programs with equa-l i t y constraints i s an adaptation of the method of multipliers as modified 80. by p a r t i a l duality. For a sufficiently large penalty constant K, the algorithm can be viewed as a primal-dual method that solves the original problem (SPE) (3.1) by solving the unconstrained dual problem (D) max d^U) (5.4) X £ Rq where the dual functional (^A) 1 S defined as (P) dK(_A) = min £(x, \y K) (5.5) x e X The objective function £(x., _X» K) and the feasible set X are as defined in ( 3 . 9 ) . It i s stated in Section 3.5 that subject to certain assump-tions, the dual functional ^(X) i s twice continuously dif ferentiable wi th VdjrCX) = h(x) (5.6) v2.dK(l> = " ^ ( x ) [^2^(x, A» K ) ] " 1 Vh(x) (5.7) where h_(x) i s the vector of signomial equality constraints, Vh i s the Jacobian of h and V 2£ is the Hessian of £(x, X^, K) with respect to x« £ i s the solution of (5.5) given _X and K. It follows that the problem (D) can, in principle, be solved by any unconstrained gradient-based ascent algorithm. In such type of algorithm, the ( k + l ) s t iterate i s obtained according to the following general form X< k + 1 ) = V k> + a ( k V k ) h (x"(X ( k ) ,K)) (5.8) (k) (k) where a i s the step size and M i s a q x q positive definite matrix. A special case of (5.8) is the fixed-step steepest ascent formula of Hestenes and Powell in which = 2K and = I, the identity matrix. Note that in ( 5 . 8 ) , K remains constant in each iteration. Also the 81. dependence of x on 1 and K is emphasized in the equation. In Section (k) (k) 5.3.2, assuming that M =1, two choices of the step size a are con-sidered. In Section 5.3.3, alternative ways of specifying K to ensure, albeit probabilistically, the global convergence of the algorithm. 5.3.2 Choice of Step Size q Assume that = I,Vk and that = K > K*, for k > k, where K* is the minimum value of K to assure global convergence and k i s some f i n i t e integer. Then (5.8) simplifies to the iteration equation of the method of steepest ascent given by ACk+D = x(k) + a ( k ) k ^ ( A ( k ) > K ) ) ( 5 > 9 ) (k) The remaining question i s : How should the step size a be selected? Before answering this question, i t must f i r s t be noted that while <i^(X) and VdK(X) = h(x(A, K)) can be calculated, the effort required to compute <L^ (^ ) amounts to solving a minimization problem. The gradient Vd^OO i s essentially a free byproduct of the minimization. For this reason, the (k) step size a should be determined as simply as i t i s possible. Schemes requiring gradient information are acceptable, provided they do not re-quire too many t r i a l points since each point calls for a minimization problem. Two step size choices are considered. The f i r s t i s the simple but effective Hestenes-Powell rule stipulating that in general, a(k) = 2 K ( k ) j K(k) ^ K ( k - l ) j s q ^ X(k+L) = x(k) + 2 K(k) h ( 5 C k ) ) (5.10) (k) Note that the penalty constant K may change from iteration to itera-(k) tion. If K = KQ,Vk, then (5.10) becomes the fixed-step gradient 82. method. The local convergence rate of (5.10) was studied by Bertsekas [27] who proved that i f (x*, X*) satisfies the second-order sufficient (k) — optimality conditions and K = K*, k £ k , for some positive number K* and integer k , then A A* at a linear rate with a convergence ratio essentially inversely proportional to provided A^°^ i s sufficiently close to A*' If = K Q ) V k , the Hestenes-Powell updating formula (5.10) (k) uses only the gradient Vd R (A) at the present point X_ to estimate the o next point A^"*^ • A possible acceleration scheme that requires a modest increase in computational effort i s to use a low-order polynomial inter-polation scheme such as f i t t i n g a cubic. However, K q may have to be increased after some iterations in order to ensure convergence. Suppose that K ( k + 1^ = K , and K ( K ^ = K , K , > K . Then Vd„ (X) and Vd„ (X) are 1 o l ~ o K l — N> ~ gradients of different dual functionals. Clearly, they cannot be used to interpolate either d^ (A) or dj, (A) • A common functional i s needed i f interpolation i s to be used and the penalty constant K can s t i l l be changed. A solution i s to use AQ(A)> i.e., d Q(X) = min £(x, A' °) (5.11) x e X (5.11) i s just the minimization of the ordinary Lagrangian. Observe that d Q a + 2Kh(x(A,K))) = d ^ a ) (5.12) Hence to approximate dg(A) ^y passing a low-order polynomial f i t through two points, dg(A_ + 2Kh(x(A>*0)) has to be evaluated for any tT.?o pairs of A and K. This i s the rationale of the following algorithm whose local convergence was proved by Bertsekas [27]. However, no comparative 83. numerical evidence was provided. Algorithm 5.1 Step 0; Set k = 1. Assume A^ 2 k ^ ' the sequence {K^} are specified. Step 1; Obtain x&k-l) by minimizing i l(x, _A ( 2 k _ 1 ) , K ( 2 k - 1 ) ) s.t. x e X. Step 2: Pk)= P k - 1 )+ 2 K ( 2 k - 1 ) h ( x ( 2 k _ 1 ) ) . It follows from (5.12) and Step 2 that d Q(A* 2 k )) and Vd Q(A ( 2 k )) are known. Step 3: Obtain x/ 2 k^ by minimizing JL(x, A_^2k\ K^ 2 k )) s.t. x e X. Hence d 0 ( A ( 2 k ) + 2K ( 2 k ) h ( x ( 2 k ) ) ) and i t s gradient are known. (2k) (2k) Step 4; Approximate dQ(A) by f i t t i n g a cubic through (Av , dQ(A. )) and ( A ( 2 k> + 2K ( 2 k> h ( x ( 2 k ) ) , d0(A<2k> + 2K ( 2 k ) h ( x ( 2 k > ) ) ) . Let the estimated step size be a Step 5; Set the f i n a l step choice a^ 2 k^ according to the following rule: \ K ( 2 k ) . f 4 K ( 2 k ) < o(2k) -(2k) I a ((2k) . f 2 K(2k) < a(2k) < 4 K ( 2 k ) 2 K(2k) . f o(2k) K 2 K(2k) c a c '„ .(2k+l) ,(2k) . ^ (2k) (2k). . . . , . Step 6: Set _A = _A + a h(x ) , k = k + 1, and go to Step 1. In Algorithm 5.1 , cubic f i t i s used every second iteration, i.e., after every even iteration. The Hestenes-Powell step size rule is (2k) used in every odd iteration. In Step 5, the estimated step size a (2k) (2k) is bounded within [2K , 4K ] because the rate of convergence ana-ly s i s by Bertsekas shows such an interval provides linear local conver-gence. 84. 5.3.3 Updating the Penalty Constant K As discussed in Section 3.6, the penalty constant K must be sufficiently large to assure convergence. Another view of this condi-tion on K i s that K must be greater than some lower bound K* in order that a local convex structure exists. The latter condition, in turn, implies that local duality theory can be used. Furthermore, as K increases, the dual convergence rate improves, provided the i n i t i a l value of X_ is close enough to the optimal Lagrange multiplier vector Xf1. In practice, however, i t i s d i f f i c u l t , i f not impossible, to determine K*. If K i s too large, the primal problem may be d i f f i c u l t to solve because of i l l -conditioning. This situation is particularly l i k e l y i f the i n i t i a l X_ i s not in the neighborhood of If seems that the only practical recourse is to start with a modest value of K, and let i t increase in some way so that both the convergence and the fast rate of convergence promised by the method of multipliers can be attained with high probability. This approach seems to be not different from the ordinary penalty method. There i s , however, a crucial difference - the consolation that K does not have to approach a for convergence to occur. (k) Two rules for choosing the sequence {K } are considered here. Their performance in numerical experiments is detailed in Section 5 . 4 . The two rules are stated below. Rule 1; K ( k + 1 ) = K ( k ) i f || h ( x ( k + 1 ) ) | | < 9 ||h(x ( k ))|| , 9 < 1.0. Else, K ( k + 1 ) = cj) K ( k ), eb > 1.0. Rule 2: K ( k + 1 ) = <J>k K, <j> •> 1.0. , Rule 1 increases the penalty constant only when a prescribed rate of con-vergence is not met. Rule 2 divorces the increase in K from the current 8 5 . condition of the problem. In both rules, the penalty constant i s i n -creased in a geometric fashion. The geometric increase, however, is a convenient arbitrary realization of the requirement that there should be an increase. 5 . 4 Numerical Experiments 5 . 4 . 1 Introduction The schemes discussed in Section 5 . 3 for updating A. and K are combined and tested in a series of numerical experiments. The experi-ments are designed to compare the convergence and the rate of convergence of each of the four possible combinations of updating rules. The numeri-cal study also looks at the sensitivity of each scheme's performance to the key varying parameter of the updating rule for K. Section 5 . 4 . 2 explains the software implementation and other related details. The test functions used in the study are given in Section 5 . 4 . 3 while the experiment plan is described in Section 5 . 4 . 4 . The results of the experiments are presented and discussed in Section 5 . 4 . 5 . 5 . 4 . 2 Software Details The software needed to perform the numerical study requires the computer code SP developed in Chapter IV and an additional subroutine incorporated into the code. The functions of this subroutine are to monitor the violations of the equality constraints and, i f required, to update _A and/or K. At the entry to this subroutine, the subproblem (irj) has been approximately solved. Since the objective function of (TTJ) i s the augmented Lagrangian £(x, _A> K) which involves a l l the equality con-straints, these constraints need not be calculated inside the subroutine i t s e l f . They are a l l implicitly calculated, in the last evaluation of £(x, _X, K) prior to the c a l l of the subroutine. The violations of the equality constraints h^.(x) = P -Qp+k(x) - 1, k = 1, 2, ..., q, can be measured by various distance func-tions. For the study of this chapter, the measure of the equality con-straint violations i s given by the Euclidean distance function <5(N) de-fined as q 5(N) = [ I h \ ( x ^ ) ] 1 / 2 (5.11) k=l K (N) where N i s the dual iteration count and x i s the solution of the Nth subproblem (T^). The equality constraints are said to be sat i s f i e d i f 6(N) < e , where e is a prescribed small positive number. Another possible measure is the Chebyshev distance max I^ QE)! u s e d by k=l,...,q fc Powell [26]. However, the Euclidean distance yields an error region smaller than that of the Chebyshev distance. Two other important tolerances are the tolerance for termina-ting the subproblem (TTj) and the f e a s i b i l i t y tolerance of the inequality constraints. Both are defined as in Section 4.4.5. Computational ex-perience acquired during the preparatory phase of the experiments re-vealed that i f the threshold for terminating the subproblem is small -4 (e.g. < 10 ) and uniform for a l l (dual) iterations, an Inordinate amount of computing time would be consumed, especially in the f i r s t couple of iterations, without improving much the overall performance. Faster con-vergence could be achieved with earlier updating of X_ and/or K. Hence a heuristic rule to accelerate convergence through controlling the thres-hold for terminating a subproblem is included in the subroutine. The rule can be described as follows: Let the threshold for terminating a subproblem be e r r p and the iteration count be k. Step 0: Set = 10 2 . 87. Step 1; If S(k) < 10 or k > 2, Anl = e^/lO.O. Step 2: If ""* Gux Lux 6(k) < 1 0 " ( 1 + J ) , 1 < j < 5, e j j * = ejJ^/lO^. The rule links the reduc-tion of the threshold to the closeness of the approximate solution to the true solution. 5.4.3 Test Problems Four signomial programs of varying sizes serve as the test pro-blems in the numerical experiments. The s t a t i s t i c s of the problems are summarized in Table 5.1. The details of each problem are described in this section. Note that the starting value of _X in a l l cases i s the zero vector with the appropriate dimension. Test Problem No. of Variables No. of Ineq. Const. No. of Eq. Const. Total No. of Terms* Degree of D i f f i c u l t y A 2 1 3 5 12 B 4 1 1 9 4 C 3 2 2 20 16 D 4 1 2 30 25 * Excludes lower and upper bounds on variables. Table 5.1 Summary of Characteristics of Test Problems. 88. The following are the four test problems: Test Problem A -2 , v -2 mm g Q(x) = x 1x 2 - 0.5 x 2 s.t. g-^x) = 0.5 x 1 1 , 8 x 2 < 1.0 g 2(x) = 2x 1 2 + 2x £ 2 + 5.0 - 6x x - 2x £ = 1.0 g 3(x) = (2/3)x± + (1/3)x2 =1.0 g 4(x) = 4x 1 2 + x 2 2 + 10.0 - 12x1 - 2x 2 = 1.0 0.01 < x < 4.0 ~ o -0.01 < x_. < 10.0, j= 1, 2 Solution: X q = 0.5, x^ ^ = 1.0, - 1.0 X± = 0.48793, X2 = 0.075643, X3 = 0.0189 Starting point: X q = 0.5, x^ = 0.5, x 2 = 1.0 Test Problem B mxn g Q(x) = 2.0 - x 1x 2x 3 s.t. g (x) = 0.25 x. + 3.75 x 0x„ + 0.375 x„x, < 1.0 J - — j 15 3 4 " g„(x) = x.. + 2x. + 2x„ + 1.0 - x, = 1.0 ^ — 1 / 3 4 0.01 < x < 2.0 o ~ 0.01 < X j < 1.0, j = 1, 2, 3 0.01 < x. < 2.0 4 -Solution: X q = 52/27, x± = 2/3, x £ = 1/3, x 3 = 1/3, x 4 = 2.0 \ = 1/9 Starting point: X q = 2.0, x± = 0.5, x 2 = 0.5, x g = 0.25, x 4 = 1.5 89. Test Problem C min g Q(x) = x J ' 6 x 2 + X 2 X 3 ° ' 5 + 1 5 - 9 8 x i + 9.0824X2, - 60.72625x, s.t. -2 -2 g 1(x) = x 1x 2 + 1.48 - x 2 x3<_1.0 g 2(x) = ^ ' 2 5 x 3 + x 2 + 6.75 - xJ' 5x 3-° <_1.0 g 3(x) = x 2 + 4.0x2 + 2.0x3 - 57.0 = 1.0 -12 5 2 g^(x) = x 1x 2 x 3' + x 2x 3- x 2 - 15.55 =1.0 1.0 < x < 1000.0 _ o _ 0.1 < x.. < 10.0, j = 1, 2, 3 Solution: X q = 319.4, x^ = 1.0, x 2 = 2.5, x 3 = 4.0 X± = -3.0, A2 = 1.5 Starting Point: X q = 558.0, x± =2.0, x 2 = 2.0, x 3 = 8.0 Test Problem D m m s.t. g Q(x) = 3x1 + 2x 2 + x 3 + x 4 + 7x 1 + 565 - 39x2 - 17x 3 § 1 ( x ) = (l/6)x 2 + (l/18)x 2 + (l/9)x 2 + (1/18)^ -(l/9)x 3 - (l/18)x 2 < 1.0 2 2 2 2 g„(x) = x.. + x„ + 2.5x, + x„ + 9.5 - 4x. - x, = 1.0 2 — 1 3 4 3 Z 4 2 2 A 2 g 3(x) = x 2 + 3x 3 + x^ + x 2 - x^ - x^ - x^ - 2.0 = 1.0 1.0 < x < 1,000.0 - o -0.1 < x j < 10.0, j = 1, ..., 4 Solution: X q = 505.0, x± = 2.0, x 2 = 2.0, x 3 = 1.0, x 4 =1.0 A 1 = -1.0, A 2 = 3.0 Starting Point: X q = 700.0, x± = 1.0, x 2 - 1.0, x 3 = 3.0, x 4 = 2.0 90. In Test Problem A, the inequality constraint g^(x) < 1.0 i s loose at the solution. The same i s true with the constraint g^(x) < 1.0 of Test Problem D. A l l the other inequality constraints of a l l four test pro-blems are tight at the quoted solutions. 5.4.4 Experiment Plan Four schemes are considered in the experiments. The f i r s t two, Schemes I and II, result from combining the Hestenes-Powell step size rule and Rules 1 and 2, respectively, for updating K. Schemes III and IV are the combination of the same two rules for updating K and the use of Algorithm 5.1 to obtain the step size. When Rule 1 i s used to update K, the i n i t i a l values K/°^ tested are 0.1, 1.0, and 10.0 . The reduction factor 8 i s 1.0 while the expansion factor <f> is set to 2.0. When Rule 2 i s used, the values of <j> are 2.0, 4.0, and 8.0. The parameter K q is set to 1.0. Of the four schemes, Scheme I is the most similar to the algori-thm proposed by Powell [ 2 6 ] . There are, however, significant differences. The f i r s t i s that while only either A or K, but not both, is updated in Powell's algorithm, both can be updated in the same iteration in Scheme I. The second difference i s the definition of the norm of h(x) as explained in Section 5.4.2. Finally, Powell assumes that M i s a diagonal matrix not necessarily equal to I, the identity matrix. Because of this more general assumption, not one but a maximum of q > 1 numbers have to be updated. To carry out the numerical experiments, two performance mea-sures need to be defined. The f i r s t i s the convergence measure while the second i s for the rate of convergence. In this study, overall con-vergence is attained i f 6(N) < e p n j . provided the Nth subproblem satisfies 91. the convergence criterion of Section 4.4.5 as modified by the heuristic rule explained i n Section 5.4.2. An alternative convergence measure i s to stop i f fi(N) < s E Q and |d } ( X ( N > ) - d ( X ^ ) | /1 d ( X ^ ) | K k K < e . Both e and E are small positive numbers. Choosing a performance measure to reflect the rate of conver-gence requires some care. A f i r s t candidate is the number of dual i t e r a -tions (i.e. updatings of _X) before convergence i s attained. This measure can be viewed as the number of evaluations of the dual functional d^(X) and i t s gradient. While this measure may indicate the dual rate of convergence, i t i s not a f a i r indicator of the computational effort re-quired to achieve a given rate. In a primal unconstrained optimization problem, a function evaluation requires the same effort at any point of the iterative process. This is not true with the evaluation of the un-constrained dual functional. In the dual case, each evaluation is a mini-mization problem whose completion depends on many factors such as the proximity of the current test point to the solution, the accuracy required, and the tolerance levels. A second candidate is the number of linearly constrained nonlinear programs solved. This measure, while i t is an improvement over the f i r s t , nevertheless i s s t i l l unacceptable because each such problem requires different computational effort. The f i n a l choice adopted i s the CPU time measured from the confirmation of a feasi-ble solution to the satisfaction of the convergence criterion''-. The time used for accepting problem data and checking/achieving f e a s i b i l i t y i s common to a l l schemes and hence is ignored. In the literature on numerical studies of optimization algorithms, the chief objection to the ''"For the rest of this chapter, this time period is simply referred to as CPU time. 92. use of CPU time is that i t s significance varies with the computer code used and with the computer system. But since these two factors are con-stant in this study, the use of CPU time seems to be the most meaningful measure for rate of convergence. Throughout the experiments the following tolerance levels are observed: e_>. = 10~4; ej££ - 10~ 2; c-. = 10~6. (Recall that eT.7 i s the EQ CGP IN IN f e a s i b i l i t y tolerance for inequality constraints.) Nonconvergence is defined by the following condition: CPU time > 3 sec or number of dual iterations > 20. 5.4.5 Results and Discussion The performance results of solving a l l four test problems by each of the four schemes are presented in Tables 5.2(a)-(d). In each box of each table, the top number is the CPU time in seconds; the middle number, the number of linearly constrained nonlinear programs solved; and the bottom one gives the number of dual iterations needed to satisfy the convergence criterion. The latter two c r i t e r i a are included for completeness. Their lack of consistent correlation with computational effort i s confirmed by the numerical data. Overall, the experiments verify the convergence of the algori-thms proposed in Chapters III and IV, irrespective of which scheme i s used to update _A and K. In a l l cases except those noted i n Table 5.2, con-vergence is' attained i n the sense that the satisfaction of the conver-gence measure implies a f i n a l primal solution of reasonable accuracy. -3 In this study, the accuracy is within at least 10 of the true solution. Even in those few cases that were terminated because of l i k e l y excessive computation, the approximate solution progresses steadily but too slowly towards the theoretical solution. Similarly, when Problem D is solved Problem^. 0.1 1.0 10.0 A 0.651 (i) 0.559 30 10 0.436 18 5 B 2.019 89 21 1.217 46 11 0.820 48 8 C 1.009 49 9 0.392 21 5 1.637 79 8 D ( i i ) 1.434 33 15 2.418 68 12 (a) Scheme I Problem^. 2.0 4.0 8.0 A 0.423 18 5 0.448 18 4 0.727 24 5 B 0.756 28 7 0.579 23 4 0.455 15 3 C 0.363 17 3 0.366 17 3 0.371 16 3 D 1.156 29 12 0.675 17 6 0.681 17 6 (b) Scheme II (i) Terminated after 20 dual iterations. ( i i ) Terminated after 3 seconds of CPU time. Table 5.2 Performance results of Schemes I.- TV 94. Problem"--^ 0.1 1.0 10.0 A 1.081 (i) 0.621 30 10 0.446 18 5 B ( i i ) 1.877 61 18 0.936 54 5 C 1.046 51 10 0.393 21 5 ( i i ) D ( i i ) 1.43 37 15 2.251 62 9 (c) Scheme III Problem^\^ 2.0 4.0 8.0 A 0.541 20 5 1.06 30 6 0.722 24 5 B 0.961 26 5 0.781 24 4 0.543 15 3 C 0.350 17 3 0.354 17 3 0.361 16 3 D 0.605 16 * 6 0.494 13 * 5 0.415 13 * 5 (d) Scheme IV (i) Terminated after 20 dual iterations. ( i i ) Terminated after 3 seconds of CPU time. (*) Approaching but did not reach primal solution when dual convergence was satisfied. Table 5.2 Performance results of Schemes I - IV 95. by Scheme IV, termination occurs because the equality constraints are satisfied, although the primal solution i s s t i l l approaching but not sufficiently near the true solution. In this case, the dual rate i s much faster than the primal rate of convergence. Some l i k e l y reasons for the discrepancy between the two rates are the inexactness in solving the subproblem and the poor performance of the primal algorithm in this particular problem. The global convergence of A to a region near the optimal Lagrange multiplier _A* is satisfactory in a l l cases. However, the local convergence to X.* i t s e l f is sensitive to the tolerances, hence accuracy, of the primal computation as well as to the step size a. While i t is essential to confirm experimentally the convergence of the proposed algorithms, a major purpose of the experiments i s to assess the effect of the updating rules of _X and K on the rate of con-vergence. To this end, the data of Table 5.2 reveal significant informa-tion. The foremost is that by far, the preferred scheme is Scheme II using the Hestenes-Powell choice for a and a monotonic geometric increase of K. From the convergence viewpoint, both step size choices are almost equally effective. However, the 'Hestenes-Powell updating formula i s preferred for i t s simplicity, less demand for computation, and indepen-dence of the numerical accuracy of other quantities. Choosing the step size by cubic f i t does not yield the benefit commensurate with the extra computational work. As for the updating rule of K, Rule 2 relying on increasing K according to a geometric progression is definitely superior to Rule 1 in which K i s increased only i f ||h_(x)|| f a i l s to decrease at a specified rate. The choice of Rule 2 i s based on two related reasons. The f i r s t i s that independent of how the step size i s selected, Rule 2 generally needs far less CPU time. This observation can be checked by 96. comparing Tables 5 . 2(a) and (b) , and Tables 5 . 2 ( C ) and (d). The second reason is that Rule 2 almost invariably yields a fast monotonic decrease in 6 ( N ) (or || h(x)||) . On the other hand, using Rule 1 seems more l i k e l y to have nondecrease or increase in 6 ( N ) . When such a condition exists, i t necessitates solving more subproblems. Figure 5 . 1 gives a specific but typical example of the behavior of 6(N) as a function of the dual iteration count N. The figure, in which log 6(N) is plotted against N, i s derived from the results of solv-ing the test problems by Scheme I ( K ^ = 1 . 0 , cj> - 2 . 0 , 6 = 1 . 0 ) and Scheme II (K Q = 1 . 0 , <|> = 2 . 0 ) . The corresponding results from Schemes III and IV are not shown because updating _X by cubic f i t i s not as promising. Note that for the indicated parameters of Schemes I and II, each problem has the same starting subproblem. The step size choice i s common to both. Also, whenever K has to be updated, i t i s increased by the same factor <i>- 2 . 0 . Hence there is a f a i r basis for comparing the two updating rules for K. Let 5j.(N) and <5I:J.(N) be the function 6(N) obtained when Scheme I and Scheme II are used respectively. Several useful observa-tions about <5j(N) and 00 can be deduced from Fig. 5.1. The f i r s t i s that in a l l cases, ^ ( N ) < 6-j-(N) . Furthermore, SJJOO tends to decrease faster than S T(N). The two points together suggest that not only does 1 -4 Scheme II lead to an earlier convergence when e c r. = 1 0 (as i t i s the case in the experiments), the scheme also consistently finishes ea r l i e r than Scheme I for a l l values of This means that whether an accurate EQ and not too accurate solution i s sought, Rule 2 should s t i l l be preferred. In both Schemes I and III, the experimental data attest to the expected sensitivity of the rate of convergence to the value K^°^. As K^°^ i s larger, the rate generally improves. The improvement is reflected in less dual iterations and a faster reduction of i|h_(x)|| . The same '.Or 0.0 \ NUMBER OF ITERATIONS N (b) Fig. 5.1 Typical plot of log 6(N) against the dual iteration count N. 9 8 . effect is true for a larger .cb in Schemes I I and I V . In both instances, the CPU time usually decreases too. There are cases, however, in which the CPU time increases. This phenomenon i s not caused by ill-conditioning in the primal problem. Instead i t can be explained as follows. When K^°^ (or cb) i s selected to be too large, the dual convergence rate is higher andj|h_(x)|| tends to decrease much faster. This means that || h_(x)|| becomes small early in the iterative process. Because of the heuristic rule controlling the tolerance e£gp> a small value of ||h(x)|| effectively implies accurate minimization of the primal subproblem. Hence more com-putational work is needed. To avoid poor convergence or excessive effort to solve subpro-blems, prudent choices of the i n i t i a l value (K in Schemes II and IV) and the factor cb must be made. It seems that for suitably scaled con-straints, K Q = 1 . 0 and cb = 4 . 0 are quite acceptable. The emphasis here i s on the suitable scaling of the equality constraints. That these cons-(k) traints should be properly scaled i s implied by assuming M = I i n Section 5.3.2. Hence the scaling of the equality constraints and that of A a r e equivalent. Since K Q = 1 . 0 , i t is advisable to scale h^(x), Vk, such that i n i t i a l l y , |h^(x)| ^ 1 . 0 , where 'V means " i n the order of". Similarly, the original objective function should not be too much larger (o) 2 or smaller than K Y l v ( x ) . k In summary, the Hestenes-Powellstep size rule and the updating of K by a geometric progression form an effective scheme to implement Step 3 of Algorithm 3.2 for solving signomial programs with equality constraints. To achieve the high convergence rate promised by this scheme, care must be taken to scale the equality constraints and the o r i -ginal objective function, Furthermore, approximate minimization of the 99. f i r s t few subproblems can save considerable computing time without diminish-ing the power of the algorithm. 100. VI. APPLICATIONS IN ENGINEERING DESIGN 6.1 Optimum Design of Statically Indeterminate Pin-Jointed Structure 6.1.1 Introduction In the last decade the use of optimization techniques i n struc-tural design has been recognized as a more powerful approach than tradi-tional direct design methods. The optimization approach can handle more design constraints and s t i l l produce the most economic solution within the limits set by the constraints. The applicability of the optimization approach spans the whole structural design process ranging from the choice of the topology of the structure to the selection of the member sizes of structures with a defined shape. This section i s devoted to the optimum design of a class of structures with a fixed shape. As i t w i l l be made clear i n the problem formulation, signomial equality constraints arise naturally in this class of structure design problems. Furthermore, both the objective function and the inequality constraints are signomials. It follows that the algorithms developed in this thesis can be conveni-ently employed to yie l d the optimum design. The usefulness of the proposed algorithms extends beyond the group of structures considered in detail in this section. There are structures of other types (e.g. r i g i d l y jointed structures) whose design requires the satisfaction of signomial equality constraints. Constraints of similar type also appear in different formulations of optimum struc-tural design. It i s , however, beyond the scope of this section to con-sider them. 6.1.2 Problem Formulation Consider a sta t i c a l l y indeterminate pin-jointed N-member struc-ture'with a specified geometry. I t is desired to choose the most economical cross-sectional area of each member such that some structural constraints are satisfied. In order to obtain the optimum design wanted, both the objective function and the constraints need to be established. The objective function i s the cost of the structure. Let SL^ and be, respectively, the length and the area of the ith member. Then, i f is the weight density of the i t h member, the member's weight i s Y.^.A.. Suppose that the cost of the it h member is c./unit weight. 1 1 1 1 The desired objective function becomes N C = I c.y.i.A. (6.1) . L . 1'1 i l i=l The structural constraints are of three types: stress con-straints, deflection constraints, and compatibility constraint. The f i r s t two types are imposed for safety reasons and are in the form of inequalities. The compatibility constraints ensure that the topology of the structure remains the same during and after the application of ex-ternal forces. These constraints are equalities. To derive a l l the constraints, i t i s f i r s t necessary to determine the member forces P_ = [P-^ » •••». A s s u m e c ^ a t there are m loading points with the external load vector L^ = [L^, L , L^]' acting at these points. In a s t a t i c a l l y determinate structure, P_ i s linearly related to and uniquely determined by L^. That i s , £ = B b 4 ( 6- 2 ) where i s known as the load transformation matrix. However, i n a sta t i c a l l y indeterminate or hyperstatic structure, L, is not sufficient —b to specify P. The forces in the redundant members have to be included. In this case, i = \ h + B A <6-3> 102. where = [R^, R^ , R Q]'» being the force in the redundant member k, and is the force transformation matrix relating P_ to L^. The stress i n a member i s defined as the ratio of the member's force to i t s cross-sectional area. For hyperstatic structures, the stress constraints are therefore given by AP = AB, L, + AB L < o — b^b r r - — (6.4) where A is an NxN diagonal matrix whose i t h diagonal entry i s 1/A^ From (6.4), i t can be seen that the stress constraint on member i i s a signomial i n A., L , L R . 1 1 2., Q To derive the deflection constraints, the deflection vector X must f i r s t be related to the load vector L_. It can be shown [19] that in general for elas t i c structures, X = FL where F, known as the overall f l e x i b i l i t y matrix, is given by B'fB. B is the force transformation matrix and f is an NxN diagonal matrix whose ith entry i s ^ -J/E-JA^, being the Young's modulus of e l a s t i c i t y of member i . In the case of hyperstatic structures, X can be partitioned into X^, the deflections due to external loads and the deflections due to the redundant forces L_, Hence (6.5) The deflection constraints impose an upper bound on some or a l l of the components of X^. Mathematically, these constraints are B b , f B b ^ ^ v f B A ^ ( 6 - 6 ) 4 V f B b B, 'fB b r h X B ' f B, B ' fB L r b r r — r 103. Again, the i t h inequality is a signomial inequality constraint in A^, R^ , • V In (6.4) and (6.6), the redundant forces L are unknown. For this reason, X^must be required to be n u l l . Otherwise, the structure would assume a different topology after the application of . Thus, B 'fB tL, + B 'fB L = 0 (6.7) r t>~t> r r—r This equality constraint i s called the compatibility constraint. (6.7) is clearly a set of signomial equality constraints. In summary, the design problem is to minimize the cost function given by (6.1) subject to the constraints (6.4), (6.6), and (6.7). 6.1.3 Example,' Consider the pin-jointed structure shown in Fig. 6.1(a) and adapted from [19]. Suppose that the horizontal and v e r t i c a l deflections at D are both limited to 4 mm, and that the numerical value of the 6 2 stress in any member is not to exceed 10 kN/m . A l l three members are one meter in length. Let A^ be the cross-sectional area of member i . It i s required to choose the cross-sectional areas A^, A^s and A^ such that the structure has minimum volume. Assume that Young's modulus of 2 el a s t i c i t y E i s 207 kN/mm . The f i r s t step is to derive the constraints (6.4), (6.6), and (6.7). To construct the matrix B, , the redundant force is removed. In b this example, member 3 may be designated as the redundant member. With member 3 removed, the structure becomes that shown in Fig. 6.1(b). Resolution of the forces at D gives Efe = [Till Jill 0]' (6.8) Next remove the external load H and subject joint D to force P^ as shown Fig. 6 . 1 The hyperstatic pin-jointed structure of the example. (a) Frame and loading (b) Basic frame with external load (c) Basic frame subject to redundant force. 105. in Fig. 6.1(c). Again, resolving the forces at D yields B =11 0 1]' (6.9) The f l e x i b i l i t y matrix f is 1/EA1 0 0 0 1/EA2 0 1/EA, (6.10) while the load vector and the redundant force vector L are simply L = P, (6.11) It follows from (6.4) that the stress constraints are -1 ,,,,.-1 (a) (Jm/2a)k~L (l/a)A 1 J- P 3 < 1 (b) (y^H/2a)A2 1 < 1 (c) ( l / a ) A 3 - 1 P 3 ^ 1 (6.12) 6 2 where 0 = 10 kN/m . From (6.6), (6.8)-(6.11), the deflection constraint may be written as (H/2E5)A1~1 + (H/2E6)A 2 - 1 + (/2/2E6)A~1 P3< 1 (6.13) 2 Young's modulus of e l a s t i c i t y E i s assumed to be 207 kN/mm and the deflection tolerance 5 i s set at 0.004 m. The compatibility constraint defined i n (6.7) i s given by (>/2H/2)A~1 + k~X P 3 + A 3 _ 1 P 3 = 0 Finally, the areas A^, A 2, and A 3 must be nonnegative. (6.14) The optimization problem i s to minimize the volume V = A^ + A^ , + A 3 subject to the inequality constraints (6.12)-(6.13) , the equality constraint (6.14), and the nonnegativity constraint on the areas. To put the problem into a form compatible with the formulation of signomial programs defined in (3.1), several steps need to be undertaken. F i r s t , from (6.14), i t can be easily shown that P 3 i s nonpositive. But signo-mial programs have only positive variables. Hence P^ has to be replaced with -Py After this substitution, (6.12c) can immediately be dropped since the constraint i s loose for a l l feasible values of A^ and P^. The f i n a l manipulation i s to impose bounds on a l l the variables. Since signomials are defined over the positive orthant only, a l l the v a r i -ables except need to be bounded from below by a small positive quan-t i t y . With the a r t i f i c i a l upper bounds imposed and the problem data entered, the f i n a l problem is min A^ + + A^ s.t. (7.0711 x 10 _ 4)A ~ 1 - 10~6A ~1A. < 1 1 1 4 ~ (6.0385 x 10" 5)A 1~ 1 + (6.0385 x l O " 5 ^ - 1 -(8.54 x 10~ 7)A 1" 1 P 3 < 1 70.7107A1~1 - A ^ 1 P 3 - A 3 _ 1 P 3 + 1 = 1 (6.15) 10~ 8 < A 1 < 1.0 7.0711 x 10" 4 < A 2 < 1.0 10~ 8 < A 3 < 1.0 10~ 8 < P 3 < 1.0 2 Note that the areas are in m while i s i n kN. The numerical example was solved with the computer code SP using Scheme II of Section 5.4.4 to update A and K. The feasible primal 2 starting point i s at ^ = A 2 = A 2 = 0.001 m and P 3 = 0.001 kN, while A ^ ° \ = 0.0. K q and <}>,. the two parameters needed for updating K, are 1.0 and 4.0, respectively. Within 7 dual iterations (CPU time = 0.66 sec), the f i n a l answer obtained is V = 0.00141412 m3, A^ = 0.000707 m2, 107. A 2 = 0.00070711 in 2, A 3 = 10 8 m2, ? 3 = 0.11159 kN, and X = 0.6. The tolerance used are = 10 2 c^. = 10 e_. = 10 4 . Since A 0 ^.0.0, CGP IN EQ J i t can be deduced from the computed solution that the true solution i s A±* = 0.00070711, A ^ = 0.00070711, A 3* = 0.0, and P 3* = 0.0. The last relation is based on the fact that P 3 = A^/(A^ + A^). The discrepancy between P 3 and P 3* is due to the finite error e in the equality constraint and the smallness of A^. In fact, i t can be shown that for negligible A 3, P 3 - e/A^. In the example, e = 7.8 x 10 ^. From the design point of view, how should the optimal solution be interpreted? The most important result i s that the proposed design concept should be changed. Instead of a three-member structure, a two-member structure is adequate (and optimal) to meet the design require-ments. At the computed solution, the stress constraints are tight while the deflection constraint is loose. This result may prompt a re-examina-tion of the allowable stress and deflection because usually the deflec-tion constraint is more significant. 6.2 Optimization Examples from Chemical Process Engineering 6.2.1 Introduction The history of applying signomial programming to steady-state chemical process design is almost as old as the optimization method's. Rijckaert [18] gives a good summary of this area's key application papers published in the period 1963-1972. A distinctive feature of the mathe-matical models used in process design is the frequent presence of signo-mials in both the cost (or profit) function and the design constraints. This a f f i n i t y of the process models to signomials is the main reason why process engineering is such a f e r t i l e f i e l d for possible application of signomial programming. Another notable feature of the models is that 108. equality constraints are very common. These constraints arise because of mass and energy balance, as well as functional relationships derived empirically or from f i r s t principles. In a l l the related papers referenced by Rijckaert [18], the models containing both equality and inequality constraints are f i r s t transformed on the basis of "engineering intuition" into one with either only equality constraints or only inequality cons-traints. The design problems are then solved by methods accommodating one type of constraints only. In this section, two examples drawn from chemical process design demonstrate how the algorithms of Chapters III and IV can be used to solve the design problems based on the original models with mixed equality and inequality constraints. Because no con-straint transformation i s required, the physical arguments for equality or inequality constraints are not blurred by algebraic manipulations. 6.2.2 Alkylation Process Optimization The example used in this section i s typical of the models used in seeking the optimum operating conditions for a chemical process i n order to maximize a pro f i t function. This section's model of an alkyla-tion process, a common process in the petroleum industry, was f i r s t described by Sauer, Co l v i l l e , and Burwick [55] on the basis of the pro-cess relationships given by Payne [56]. The optimization problem was solved by Sauer et a l . [55] via the solution of a sequence of linear pro-grams. Later Bracken and McCormick [57] solved the problem using the penalty function method (SUMT). Recently, the same model, after i t had been transformed into an inequality-constrained signomial program, was solved by Avriel et a l . [12] using the computer code GGP developed by Dembo [16]. The mathematical formulation given below follows that of Bracken and McCormick. Figure 6.2 shows a simplified process diagram of an alkylation process. The system's main units are the reactor and the fractionator. Olefin feed and isobutane make-up are pumped into the reactor. To cata-lyze the reaction, fresh acid i s added and the spent acid i s withdrawn. The hydrocarbon product from the reactor is fed into a fractionator, from the top of which isobutane is recycled back to the reactor. Alky-late product i s withdrawn from the bottom of the fractionator. It is assumed that the olefin i s pure butylene, that both the isobutane make-up and the isobutane recycle are pure butane, and that the fresh acid strength i s 98% by weight. The process variables and their upper and lower bounds are defined in Table 6.1. The bounds on the variables x^ are due to the limitations imposed by the capability of the plant and/or the economic situation under analysis. For example, only 2,000 barrels per ISOBUTANE RECYCLE OLEFIN FEED ISOBUTANE MAKE-UP FRESH ACID REACTOR HYDROCARBON PRODUCT .SPENT ACID FRACTIONATOR ALKYLATE PRODUCT Fig. 6.2 Simplified alkylation process diagram. 110. VARIABLE LOWER BOUND UPPER BOUND x l o lefin feed (barrels per day) 1.0 2000.0 isobutane recycle (barrels per day) 1.0 16000.0 X3 acid dilution rate (1000 lb/day) 1.0 120.0 X4 alkylate yield (barrels per day) 1.0 5000.0 x5 isobutane make-up (barrels per day) 1.0 2000.0 X6 acid strength (weight per cent) 85.0 93.0 X7 motor octane number 90.0 95.0 X8 external isobutane-to-olefin ratio •3.0 12.0 X9 acid dilution factor 1.2 4.0 xio F-4 performance number 145.0 162.0 Table 6.1 The alkylation process problem's variables and their bounds. PROFIT AND COST PARAMETER VALUE alkylate product value olefin feed cost isobutane recycle cost acid addition cost isobutane make-up cost $0,063 per octane-barrel $5.04 per barrel $0,035 per barrel $10.00 per thousand pounds $3.36 per barrel Table 6.2 Definitions and values of the cost coefficients of the alkylation process problem. 111. d a y o f o l e f i n f e e d m a y b e a v a i l a b l e . T h e o t h e r v a r i a b l e s h a v e b o u n d s d i r e c t l y r e l a t e d t o t h e p r o c e s s . T h e p r o c e s s c o n s t r a i n t s c o n s i s t o f t w o g r o u p s : s t a t i s t i c a l r e l a t i o n s a m o n g s o m e v a r i a b l e s w i t h i n c e r t a i n o p e r a t i n g r a n g e s , a n d r e l a t i o n s h i p s d u e t o m a t e r i a l b a l a n c e s . T h e f i r s t g r o u p , d e s c r i b e d b y l i n e a r o r n o n l i n e a r r e g r e s s i o n s , l e a d s t o i n e q u a l i t y c o n s t r a i n t s w h i l e t h e l a t t e r g r o u p i s j u s t a s e t o f e q u a l i t y c o n s t r a i n t s . E a c h r e g r e s s i o n r e l a t i o n s h i p i s r e p l a c e d w i t h t w o i n e q u a l i t y c o n s t r a i n t s , w h i c h s p e c i f y t h e r a n g e f o r w h i c h t h e r e g r e s s i o n r e l a t i o n s h i p i s v a l i d . C o n s i d e r t h e r e g r e s s i o n e q u a t i o n Y = f ( z ) . T h i s r e l a t i o n s h i p w o u l d b e e x p r e s s e d i n t h e m o d e l a s d „ Y < f ( z ) < d Y ( 6 . 1 6 ) J6 XX O r f ( z ) < d Y — u - , f ( z ) < - d ^ Y ( 6 . 1 7 ) T h e d e v i a t i o n p a r a m e t e r s d ^ a n d d ^ e s t a b l i s h t h e p e r c e n t a g e d i f f e r e n c e o f t h e e s t i m a t e d v a l u e f r o m t h e t r u e v a l u e . I n t h i s m o d e l , t h e v a l u e s o f t h e d a n d d f o r e a c h r e g r e s s i o n r e l a t i o n a r e s e t t o t h e s a m e v a l u e s I u a s t h o s e g i v e n i n [ 5 7 ] ' . A s s u m i n g t h a t t h e r e a c t o r t e m p e r a t u r e i s b e t w e e n 8 0 ° F a n d 9 0 ° F , a n d t h a t t h e r e a c t o r a c i d b y w e i g h t p e r c e n t s t r e n g t h i s 8 5 - 9 3 % , n o n l i n e a r r e g r e s s i o n a n a l y s i s r e l a t e s t h e a l k y l a t e y i e l d x ^ t o t h e o l e f i n f e e d x ^ a n d t h e e x t e r n a l i s o b u t a n e - t o - o l e f i n r a t i o n x D b y t h e f o l l o w i n g p a i r o f o i n e q u a l i t i e s : (99/100)x, < 1.12X, + 0 . 1 3 1 6 7 x 1 x o - 0 . 0 0 6 6 7 x . , x 2 < ( 1 0 0 / 9 9 ) x , ( 6 . 1 8 ) H — 1 1 O 1 O — 4 S i m i l a r l y , u n d e r t h e s a m e r e a c t o r c o n d i t i o n s , t h e m o t o r o c t a n e n u m b e r x ^ i s r e l a t e d t o X g a n d t h e a c i d s t r e n g t h b y w e i g h t p e r c e n t x ^ b y 112. (99/100)x-, < 86.35 + 1.098xo - 0.038x2 + 0.325(x.-89) < (100/99)x_ (6.19) / — O O D — / A linear regression relation expresses the acid dilution factor X g as a linear function of the F-4 performance number X ^ Q , which in turn i s linearly related to x^ by regression analysis. (9/10)x9 £ 35.82 - 0.222x1Q < (10/9)xg (6.20) (9/100)x 1 Q < -133 + 3xy < (100/99)x 1 Q (6.21) Mass balance requires equality constraints for the isobutane make-up x^, the acid dilution factor x^, and the external isobutane-to-olefin ration x Q. Assuming that the volumetric shrinkage i s 0.22 volume o per volume alkylate yield, the isobutane make-up x,_ may be determined by a volumetric reactor balance given by x. = x 1 + x r - 0.22x. (6.22) 4 1 5 4 or ' xr-= 1.22x, - x. (6.23) 5 4 1 The acid dilution factor X g may be derived from an equation expressing acid addition rate x^ as a function of alkylate yield x^, acid dilution factor X g , and acid strength by weight percent x^. Thus, l,000x 3 = x 4x 9x 6/(98 - x 6) (6.24) or x g = gS.OOOx^^x"1 - l.OOOx^"1 (6.25) Finally, by definition, x g = (x 2 + x 5)/x L (6.26) Combining (6.23) and (6.26), and rearranging yields X2 = X1 X8 " 1 , 2 2 x 4 + x i (6-27) The profit function is taken as the value of the output, the 113. alkylate yi e l d , minus the feed and recycle costs. Other operating costs are assumed to be constant. The total daily profit to be maximized is P = C . X . X . . - C _ X T - c„x„ - c,x„ - c cx_ (6.28) 1 4 7 2 1 3 2 4 3 5 5 where the cost coefficients c^, c,. are as defined in Table 6.2. Also given in the same table are values of the coefficients used in the numerical solution. To cast the design problem into the format compatible with the formulation of Chapter III, the following are done: 1) The maximization of P i s replaced with the minimization of x o subject to x o^(-P + c) _< 1, where c i s a sufficiently large positive constant (e.g. 3,000) to insure that X q _< 0. 2) Each of the constraints (6.18)-(6.21) i s rewritten as a pair of signomial inequality constraints. 3) Each of (6.23), (6.25), and (6.27) i s stated as a signomial equality constraint. The resulting signomial program i s min x o s.t. (a) x~ 1(5.04x 1 + 0.035x2 + 10x3 + 3.36x5 + 3,000 - 0.063x4x7 (b) 0.00667x x2+ 0.99x( + 1 - 1.12x - 0.13167x x < 1 1 8 4 1 1 8 -(c) 1.12x1 + 0.13167x1xg+ 1 - 0.00667xjXg - l.OlOlx^ < 1 (d) 0.038X2, + 0.99X., - 0.325x£ - 1.098xo - 56.425 < 1 o / o o — (e) 0.325x£ + 1.098xo + 58.425 - 0.038x2 - l.OlOlx, < 1 o o o / — (f) 0.222x + 0.9xn - 34.82 < 1 1 0 9 - ( 6 > 2 9 ) (g) 36.82 - 0.222x 0 - l . l l l l x < 1 (h) 134 + 0.99xlQ - 3:<? < 1 (i) 3x ? - l.OlOlx 0 - 132 < 1 114. (j) 1.22x^x5 - x xx 5 =1 (k) 98,000x 3x~ 1x~ 1 Xg 1 - 1.22x21x4 =1 (1) x ^ ^ X g - x 1x~ 1 - 1.22x 2 1x 4 = 1 (m) the bounds specified in Table 6.1. The code SP was used to solve (6.29) with X_ and K updated accor-ding to Scheme II. As in Section 6.1, = (0.0, 0.0, 0.0)', K q = 1.0, and cb = 4.0. The optimality and f e a s i b i l i t y tolerances are also the same as those in Section 6.1. In the f i r s t few computation runs, i t was discovered that the dual convergence was not satisfactory in the sense that 6(N) was either constant for several dual iterations or decreasing very slowly in spite of K being very large (e.g. 10^). A modification that corrected the situation is to rewrite the equality constraints as (j) 1.22x4 - x x - x 5 + 1.0 = 1.0 (k) 98000.0x_x, - 1000.0x„ - x.x. + 1.0 = 1.0 3 6 3 4 9 (1) X j X g - X ; L - 1.22x4 - x 2 + 1.0 = 1.0 The above equations, however, are poorly scaled. The f i n a l scaled equa-tions used are (j) 1.22x4 - x± - x 5 + 1.0 = 1.0 (k) gso.ox-x"1 - 10.Ox. - O.Olx.xn + 1.0 = 1.0 3 6 3 4 9 (1) O . l X j X g - 0.1X;L - 0.122x4 - 0.1x2 + 1.0 = 1.0 Table 6.3 gives the starting point and the optimal solution computed. Note that the i n i t i a l point Is feasible with respect to the inequality constraints. The solution obtained here is very close but not identical to that of Bracken and McCormick [57]. The slight dis-crepancy is perhaps due to the primal algorithm used. As can be expected from the nature of the constraints (6.18)-(6.21), four of the inequality constraints of (6.29) are loose at the optimum. STARTING POINT OPTIMAL SOLUTION xo 2200.0 1232.88 x l 1745.0 1696.78 X2 12000.0 15805.58 X3 110.0 54.07 x4 3048.0 3028.87 X5 1974.0 1998.44 x6 89.5 90.11 X7 92.8 95.0 X8 8.0 10.49 X9 3.6 1.56 x10 145.0 153.53 Table 6.3 The i n i t i a l and optimal solution of the alkylation process problem. Since the profit P is C - x^, the maximum prof i t i s deduced to be 1767.11 dollars/day. From Table 6.3, i t can be seen that one of the inputs, the isobutane make-up x^, is at i t s upper bound. This suggests that the profit may be increased further by allowing more iso-butane make-up, the maximum profit is raised by about 5% to 1857.23 dollars/day. To accomplish this, the capacity of the fractionator has to be expanded to yield the isobutane recycle flow of 16593.27 barrels/ day. 116. The computational experience with the alkylation problem offers some valuable insight into how an equality constraint of the form y = f(x), where f(x) i s a signomial, should be converted to a standard signomial equality constraint g(x) = 1.0. Such type of equation often arises in input-output relations and in functional relations defining one design variable in terms of other variables. The f i r s t way to transform the equation is to write y ^f(x) = 1.0. This format i s usually well scaled provided y i s i n the proper range. However the factor y ^'^ appears in every term of f(x) and the augmented objective function i s very nonlinear in y. The second way i s to rewrite the equation as f(x) - y + 1.0 = 1.0. The advantage of this format is that the augmented cost i s quadratic in y. This format, however, is susceptible to poor scaling. If the latt e r form i s to be used, a scale factor 3 must be introduced such that, for example, i n i t i a l l y B|f(x) - y[ ^ 1.0. The f i n a l form becomes 8f(x)- 8y + 1.0 = 1.0. As long as 3 is suitably selected, the second format seems to be preferable. 6.2.3 Design of a Heat Exchanger Network The design of the heat transfer capability of a process varies very widely in complexity. The design task may range from the size specification of a simple heat exchange unit to the optimal structural synthesis of a heat exchanger network. The example in this section i s concerned with the optimal choice of the heat transfer areas and the temperature variables in a small heat exchanger network with a known configuration. Suppose that in an industrial chemical process, a flow denoted as Stream 2 has to be cooled from T 2 1°F to T 2 2°F. It i s found that the required heat loss by Stream 2 cannot be completely absorbed by an 117. existing cooler C through which cooling water i s pumped. There are, however, two streams in the process, Streams 1 and 3, which need to be heated and which can be conveniently routed to the v i c i n i t y of Stream 2. The two streams are therefore potential sinks for absorbing via heat exchangers part of the heat of Stream 2. Such a scheme, as shown dia-grammatical l y in Fig. 6.3, not only makes i t feasible to reduce the tem-perature of Stream 2 to "Y.^^* t> u t also may be cheaper than cooling by water alone, assuming the water flow i s fast enough. The design goal i s to satisfy system constraints and s t i l l achieve the optimal trade-off between the exchangers' installation costs and the operating cost of using cooling water. The exchangers are assumed to be the double-pipe type with countercurrent flow. The design variables are the output temperatures t^, t^ as indicated in Fig. 6.3, and the heat transfer areas and A ^ of exchangers XI and X2, respectively. The annual cost i s the sum of the exchanger's annual depreciation and the annual cost of cooling by water. The cost of in s t a l l i n g a heat exchanger with a heat transfer area A i s STREAM 2 HEAT EXCHANGER XI STREAM 1 © T, HEAT EXCHANGER X2 u STREAM 3 Q \ i COOLER C 1 COOLING WATER w Fig. 6.3 A small heat exchanger network'. 118. i s generally taken as hA dollars [58] for some positive constants h and y. If the cost i s spread over i t s expected service lifetime of n years, the annual depreciation cost is cA dollars/year where c = h/n. Because Streams 1 and 3 have different corrosion properties, Xl and X2 need to be constructed with different materials, and hence have different cost parameters c and y. If the annual cost rate of cooling by water i s c^ dollars/°F, then i t follows that the total annual cost f(A^, A^, t^) fCA^ A 2, t 4) = c ^ J l + c 2A 2 y2 + c 3 ( t 4 - T 2 2) (6.30) The constraints on the system of Fig. 6.3 consist of four groups: 1) Bounds on the minimum temperature difference (called the approach) between the hot and cold streams of an exchanger; 2) Bounds on the output temperature of the cold streams and on the areas of the heat exchangers; 3) Energy balance equalities relating the input and output tem-peratures of both streams; 4) Energy balance equalities relating the heat transferred and the heat transfer area of an exchanger. Since the temperature profile of the hot and cold streams of an exchanger is not known, the approach T cannot be e x p l i c i t l y bounded. An approximate bound on T may be obtained by bounding the differences T. - T . and T, . - T . The subscripts o and i refer to the output ho c i h i co v v temperature and input temperature respectively while the subscripts h and c refer to the hot and cold streams. Hence, from Fig. 6.3, the follow-ing inequalities must hold. 119. (a) T 2 1 ' - fcl > V (b) C3 " T l l > V (c) C3 " '2 > V (d) fc4 - T31 > V ( e ) fc4 " fc5 > T c (6.31) T i s the minimum allowable approach for the heat exchanger while T i s that for the cooler. Because Streams 1 and 3 are required for other purposes, the temperatures t^ and t^ may have to be bounded. It i s assumed that t^ <_ and a 2 £ t 2 £ b 2« Furthermore, because of heat pollution control, the cooling water cannot be heated beyond T W°F. Space limitation also demands that the heat exchangers' areas must be within 100 sq. f t . Let and C p^ be, respectively, the flow rate (lb/hr) and the heat capacity (BTU/lb-°F) of Stream i . Also, l e t the water's flow rate be Ww and i t s heat capacity be C • Then energy balance requires ( a ) W 2 C P 2 ( T 2 1 " ^ " V p l ^ l - V (b) W 2 C p 2 ( t 3 - t 4 ) - W 3C p 3(t 2 - T 3 1) (6.32) (c) W 2C p 2(t 4 - T 2 2) - W wC w(t 5 - T J Furthermore, the heat transfer rate Q of a heat exchanger can also be expressed as Q = UA AT. (6.33) log 2 where U i s the overall coefficient of heat transfer (BTU/hr-ft -°F), A i s the area available for heat transfer, and AT, is the logarithmic log 6 mean temperature difference. In terms of the input and output tempera-tures of both hot and cold streams, the logarithmic mean temperature 120. difference is defined as (T, . - T ) - (T - • T .) T _ h i co ho c i ((. . log ~ lnI(T, . - T )/(T, - T .)] K X > ' ^ } & hi co ho c i J where and T ^ are the input temperatures of the hot and cold streams, respectively, and T^ and are the output temperatures. Frequently AT, i s replaced with the arithmetic mean [(T, . - T ) - (T, - T .)]/2.0. log ^ h i co ho c i This approximation has less than 1% error i f ( T h i - T c o ) / ( T h o - T c i ) i s > 0.7 and less than 10% i f (T, . - T )/(T, - T .) is > 0.33 [59]. For — h i co ho c i — the purpose of this example, the arithmetic mean i s used to y i e l d (a) WjC ^ - T u ) = 0.5U 1A 1[(T 2 1 - t ^ + ( t 3 - T n > ] (b) W 3C p 3(t 2 - T 3 1) = 0.5U 2A 2[(t 3 - t 2) + ( t 4 - T 3 1 ) ] (6.35) In summary, the design problem calls for the minimization of (6.30) subject to the inequality constraints (6.31), the bounds on t^, t 2 > t^, A^ and A^, and the equality constraints (6.32) and (6.35). Table 6.4 gives the required data used to solve the numerical example presented i n this section. To simplify and put the optimization problem into the proper format, the following are done: 1) The constant term _°3^22 ^ S ^ r oPP e <* r r o m t n e objective function. 2) Constraint (6.31c) i s ignored since T 2 2 - x c = 260°F > b w = 180°F _> ty 3) Constraints (6.32a) and (6.32c) are used to eliminate t 3 and t^. The resulting signomial program is min 350A°_*6 + 275A°' 8 + 142.5t5 s.t. (a) (9.634669 x 10" 4)t 1 + (1.67124 x 1 0 _ 3 ) t 2 < 1.0 (b) 0.72251t2 + 0.95t5'•+ 0.5765t]L - 633.21639 = 1.0 (c) t± + 0.010255A1t1- 5.58343A1 - 239.0 = 1.0 121. Stream 1 Input temperature Flow rate Heat capacity Upper bound of output temp. Stream 2 Input temperature Output temperature Flow rate Heat capacity Stream 3 Input temperature Flow rate Heat capacity Range of output temp. Codling Water Input temperature Flow rate Heat capacity Upper bound of output temp. Annual cooling cost Heat Exchanger Xl Installation cost Overall heat transfer coefficient U Heat Exchanger X2 Installation cost c^k Overall heat transfer coefficient U L l l • W, "pi L21 r 22 w. 31 P3 a2^2^ b2 w W w pw w C 1 A Minimum approach for exchangers Minimum approach for cooler 240 °F 23,060 Ib/hr. 0.5 BTU/lb °F 550 °F 480 °F 280 °F 25,000 lb/hr. 0.8 BTU/lb °F 278 °F 20,643 lb/hr. 0.7 BTU/lb °F 300 °F < t 2 <_ 350 °F 100 °F 19,000 lb/hr. 1 BTU/lb °F 180 °F 150$/°F yr. 350A°*6$/yr. 150 BTU/hr. f t . ' 275A°'8$/yr. 2 oT x 150 BTU/hr. f t . " °F 20 °F 20 °F Table 6.4 Data for the heat exchanger network design problem. 122. (d) t £ + 0.0051903A2t2 + 0.0029922A2t1 - 2.72676A2 -0.0049308A2t5 - 277.0 = 1.0 (e) 0 <_ t± < 500.0 (f) 300.0 < t 2 <. 350.0 (g) 119.0 < t 5 <_ 180.0 (h) 0 < A. < 100.0, i = 1, 2 (6.36) As in the previous two design problems, the optimal solution i s obtained by the code SP using Scheme II to update X. and K. The parameters of Scheme II are as before. Table 6.5 gives both the starting point and the optimal solution (6.36). The optimal values of t^ and t^ are obtained from t^ and t^ according to the following linear relations derived from (6.32a), (6.32c) and the data of Table 6.4. t 3 = -0.5765t1 + 618.0 (6.37) t 4 = 0.95t5 + 185.0 Thus at the optimum t 3 = 391.13°F and t $ = 355.9°F. It can be easily verified that the error due to replacing the logarithmic mean with the arithmetic mean is less than 5% in exchanger XI and less than 1% i n exchanger X2. Thus the approximation by the arithmetic mean is reason-able. The energy balance equations (6.32) and (6.35) are a l l satisfied within a relative discrepancy of only 0.02%. From the optimal solution, i t can be seen that t,-, the cooling water's output temperature, is almost at the maximum allowable value. This confirms that even at the optimum, the cooler by i t s e l f cannot absorb the heat load. The solution also reveals the interesting result that i t i s more economical to have A^ at i t upper bound than to have a higher output temperature t 1 in exchanger Xl. 123. STARTING POINT OPTIMAL SOLUTION Cost 313226.83 24787.63 h 300.0 394.15 T 2 325.0 326.76 150.0 179.9 A l 50.0 100.0 A 50.0 66.04 2 Table 6.5 S t a r t i n g point and optimal s o l u t i o n of the heat exchanger network problem. VII.CONCLUSION 7.1 Synopsis A wide range of engineering design problems can be modeled as signomial programs with both equality and inequality constraints. The theory of signomial programming, however, has been developed in terms of inequality constraints only. In this thesis, an algorithm is proposed and implemented to solve signomial programs with mixed inequality and equality constraints. The algorithm does not require i t s user to trans-form the equality constraints to inequalities. The equality constraints are automatically incorporated into the original objective function to reduce the original problem into a sequence of less constrained signomial programs with inequality constraints only. Because i t is within the framework of signomial programming, the proposed algorithm i s computatio-nally practical. Because no problem manipulation i s needed prior to data entry, the physical significance of the original problem formulation i s preserved. In Chapter II, a concise review of signomial programming i s given. Signomial programs are defined and shown to be nonconvex. Because of nonconvexity, the solution of a signomial program can at best be gua-ranteed as a local minimum. Another consequence of nonconvexity is the existence of only a weak duality theorem relating the Kuhn-Tucker points of a primal signomial program to the Kuhn-Tucker points of a linearly constrained dual problem. At these Kuhn-Tucker points, the primal and dual objective functions are equal. An important special class of signo-mial programs is that of geometric programs. It is shown that geometric programs can be changed via the logarithmic transformation into convex 125. programs. As a result, geometric programming possesses an elegant strong duality theorem which identifies the global minimum of a geometric program as the global maximum of a concave dual function restricted to a linear manifold. The duality theorem thus offers two options, the primal and the dual approach, for computing the solution of a geometric program. Nei-ther option can be categorically stated as being superior in a l l cases. However, given an adopted method to solve geometric programs, the numeri-cal solution of a signomial program can be obtained by successively appro-ximating the nonconvex signomial program with convex geometric programs* each of which can be conveniently solved. This i s the strategy of the promising Avriel-Williams algorithm used in a slightly different form i n the later chapters of the thesis. In the Avriel-Williams algorithm, the approximation i s accomplished by condensing the signomials' negative terms into monomials. The main problem of signomial programs with equality and inequa-l i t y constraints i s formulated and solved in Chapter III. The proposed a l -gorithm, Algorithm 3.2, is the synthesis of three basic guiding concepts: the multiplier method viewed as a primal-dual method, partial dualization, and preservation of compatibility with signomial programming. In l i e u of confronting the original problem, the algorithm solves a sequence of i n -equality-constrained signomial programs each specified by the Lagrange multiplier estimate vector X_ and the penalty constant K. Alternating with the solution of each signomial program of the sequence i s the updating of X'and K according to some rules. The convergence of Algorithm 3.2 i s as-sured by that of the method of multipliers and by the convergence of the Avriel-Williams algorithm used to solve the sequence of signomial programs. It is observed that expressing each problem of the sequence in the stan-126. dard format of a signomial program would l i k e l y entail inconvenience in data preparation and computational d i f f i c u l t y due to large problem size. To circumvent these problems, three indirect methods of applying the Avriel-Williams algorithm to the subproblem (3.9) have been numerically explored. A l l three are based on monomial condensation, require no or small increase in variables and constraints, and need only the minimum effort for data entry. Of the three methods tested, one is selected as promising and dis-cussed i n detail i n Chapter IV. The consideration of the other major step of Algorithm 3.2, the updating of X_ and K, is deferred to Chapter V. Chapter IV presents a new numerical method to implement a proposed variant of the Avriel-Williams algorithm for solving signomial programs. In this method the original nonconvex feasible set defined by signomial inequalities is approximated by a convex set obtained by monomial conden-sation. The nonconvex objective function is then minimized over the convex approximant by a combined reduced gradient and cutting plane algorithm. The minimization's solution in turn determines the next convex approximant. In contrast to this method is the original version of the Avriel-Williams algorithm in which the signomial objective function is f i r s t incorporated into the feasible set prior to the convex approximation. The proposed method seems to match better the need of signomial programs characterized by a high degree of d i f f i c u l t y with most of the terms being in the objec-tive function. Unlike other methods reported in the literature, the method here readily admits extension to problems with a nonsignomial objective function but signomial inequality constraints. A practical example of such a problem is the constrained location problem solved in Section 4.7.3. The software implementation of the combined reduced gradient and cutting plane method requires detailed considerations of several aspects 127. of the algorithm. Section 4.4 considers the questions relating to how to obtain an i n i t i a l basic feasible solution, how to generate an i n i t i a l basis of the linear constraints, how to augment the basis' inverse and get a basic feasible solution after adding a cut, how to perform the linear search, and how to update the basis. In the same section, the optimality c r i t e r i a and the f e a s i b i l i t y tolerances are defined. Computational expe-rience with the software implemented shows that the proposed numerical method compares favorably with another well tested recent implementation. However, the method suggested in this thesis has the added advantageous f l e x i b i l i t y of handling general nonconvex objective functions such as the algebraic functionals of signomials demonstrated in Section 4.7. The formulation and testing of the updating rules for X_ and K are reported in Chapter V. The two rules for updating X_ are the Hestenes-Powell rule and one based on f i t t i n g a cubic to the ordinary Lagrangian. The parameter K i s increased by a factor <j> > 1 either only when the norm of the equality constraints does not decrease sufficiently, or in every dual iteration independent of the behavior of the equality constraints. Altogether four schemes are tested i n a series of numerical experiments. In terms of rate of convergence, the most promising scheme i s that using the Hestenes-Powell rule for A_ and increasing K according to the second criterion. The impact of the range of values of K on the rate of conver-gence i s also studied. It appears that best results are obtained when K is i n i t i a l l y small and then is moderately increased, provided the equa-l i t y constraints are i n i t i a l l y scaled proportionately. Furthermore appro-ximate minimization in the f i r s t few iterations can reduce considerably computing time without any serious loss of convergence and the required accuracy. 1 2 8 . In Chapter VI, the algorithms discussed in this thesis are applied to three engineering design problems. The f i r s t design problem involves choosing the cross-sectional areas of a hyperstatic pin-jointed structure such that the total volume is minimized subject to stress, deflection and topology constraints. The second example required the specification of the most profitable operating condition of an alkylation process commonly found in the petroleum industry. The last problem aims at achieving the optimal trade-off between installation and operating costs of a small heat exchanger network that is typical of those used in process engineering. Each design problem is f i r s t formulated from i t s engineering point of view before the numerical solution i s obtained and interpreted. 6.2 Contributions of the Thesis The research reported in this thesis is the f i r s t known effort to develop an algorithm for solving signomial programs with mixed inequa-l i t y and equality constraints without having to transform one type of cons-traints to the other. Furthermore, the algorithm i s compatible with the existing methods of signomial programming. The specific original contri-butions of this thesis may be lis t e d as follows: 1. the proposal of an algorithm that can automatically handle both inequality and equality signomial constraints within the frame-work of signomial programming and without any need for user mani-pulation of the constraints; 2. the development of a new numerical method to realize a variant of the Avriel-Williams algorithm for signomial programs whose objective function may have many terms or is an implicit signomial; 3. the extension of the numerical method of Item 2 to include non-signomial objective functions such as algebraic functionals of signomials, and the application of the method to constrained location-allocation problems; 4. the software implementation of the aforementioned algorithms, and the acquisition of computational experience with the imple-mentation; 5. the numerical study of the updating schemes for the parameters X_ and K yielding useful guidelines for later applications of the algorithms; 6. the formulation and numerical solution of selected r e a l i s t i c design problems cast as signomial programs with mixed types of constraints. 6.3 Suggestions for Further Research Further research based on the results of this thesis may be pursued along several directions. The f i r s t i s to investigate al t e r -native more efficient ways of solving the primal subproblem (3.9). The solution of the subproblem consumes a good part of the total computing time. Hence the a b i l i t y to solve the subproblem faster i s definitely desirable. Faster convergence may be achieved In various ways. Within the present framework of the combined reduced gradient and cutting plane (CRGCP) algorithm, effort may be focused on improving the inexact linear search, the procedure of obtaining simultaneously a basis and i t s inverse in the presence of tight constraints, and the method of obtaining a fea-sible and better point after the addition of a cut. A more involved study i s to consider incorporating a quasi-Newton feature into the CRGCP algorithm so that superlinear convergence may replace the present algorithm's linear convergence. 130. Another item worthy of further research is to consider using the CRGCP algorithm for solving the class of problems characterized by a nonconvex objective function and convex constraints. This study would include devising procedures for obtaining an i n i t i a l feasible (possibly interior) point and comparing different cutting plane schemes so as to achieve f e a s i b i l i t y earlier. More investigation is also needed on the specification and the updating of X_ and K. For example, the i n i t i a l values of A. and K may be more judiciously selected on the basis of primal f e a s i b i l i t y or optimality. K may be allowed to be different with each of the q equality constraints, and then the q penalty constants may be automatically adjusted in an adaptive fashion to maintain proper scaling in the dual space. The algorithms proposed in this thesis may also be used to solve wider classes of problems to extend the applicability of signomial programming. One such class is that of "pseudo" signomial programs a few terms of which have transcendental functions. Clearly such problems are not signomial programs. However, through signomial approximation of the transcendental functions and the use of new variables and equality constraints, the programs can be closely approximated by true signomial programs with equality constraints. Hence they can be solved by the a l -gorithms of this thesis. The approximation schemes need to be explored, and i t would be very convenient i f the approximation and the required new variables and equality constraints could be internally generated by the software package. 131. REFERENCES 1. R.J. Duffin, E.L. Peterson, C. Zener, Geometric Programming, New York: Wiley, 1967. 2. W.I. Zangwill, Nonlinear Programming: A Unified Approach, Englewood C l i f f s , N.J.: Prentice-Hall, 1969. 3. R.J. Duffin, "Linearizing geometric programs", SIAM Review, 12, 211-227, Apr i l , 1970. 4. U. Passy and D.J. Wilde, "Generalized polynomial optimization", SlAM J. Appl. Math., 15, 1344-1356, September, 1967. 5. M. Avriel and A.C. Williams, "Complementary geometric programming", SIAM J. Appl. Math., 19, 125-141, July, 1970. 6. R.J. Duffin and E.L. Peterson, "Geometric programming with signomials", J. Optim. Th. Applic., 11, 3-35, 1973. 7. G.W. Westley, "A geometric programming algorithm", Report No. ORNL-4650, Oak Ridge National Laboratory, Oak Ridge, Tennessee, July, 1971 8. T. Jefferson, "Geometric programming with an application to transporta-tion planning", Ph.D. Dissertation, Dept. of Operations Research, Northwestern Univ., Evanston, I l l i n o i s , August, 1972. 9. P.A. Beck and J.G. Ecker, "Some computational experience with a modi-fied convex simplex algorithm for geometric programming", Armament Development and Test Center, Elgin A.F.B., Florida, ADTC-72-20. 10. G.A. Kochenberger, R.E.D. Woolsey, and B.A. McCarl, "On the solution of geometric programs via separable programming", Operational Research Quarterly, 24, 285-294, June, 1973. 11. G.A. Kochenberger, "Geometric programming - extension to deal with degrees of d i f f i c u l t y and loose constraints", Ph.D. Thesis, School of Bus. Admin., Univ. of Colorado, 1969. 12. M. Avriel, R. Dembo and U. Passy, "Solution of generalized geometric programs", Int'l J. for Numerical Methods in Eng., 9_, 149-168, January, 1975. 13. A. Charnes and W.W. Cooper, "A convex approximant method for nonconvex extensions of geometric programming", Proc. Nat'l Academy of Sciences, 56, 1361-1364, 1966. 14. R.J. Duffin and E.L. Peterson, "Reversed geometric programs treated by harmonic means", Indiana Univ. Math.J., 22, 531-550, December, 1972. 15. U. Passy, "Condensing generalized polynomials", J. Optim. Th. A p p l i c , 9_, 221-237, 1972. LSI. 16. R.S. Dembo, "Solution of complementary geometric programming problems", M.Sc. Thesis, Dept. of Chem. Eng., Technion, Haifa, Israel, 1972. 17. G.E. Blau. and D.J. Wilde, "Generalized polynomial programming", Canad. J. Chem. Eng., 47, 317-326, June 1969. 18. M.J. Rijckaert, "Engineering applications of geometric programming", Optimization and Design, ed. by M. Avriel et a l . , Englewood C l i f f s , N.J.: Prentice-Hall, 196-220, 1973. 19. K.I. Majid, Optimum Design of Structures, London-Newner-Butterworths, 1974. 20. U. Passy, "Nonlinear assignment problems treated by geometric programming", Operations Research, 19, 1675-1689, 1971. 21. G. Blau and D.J. Wilde, "A Lagrangian algorithm for equality constrained generalized polynomial optimization", A.I.Ch.E.Journal, 17, 235-240, January, 1971. 22. J. Yan and M. Avriel, "Signomial programs with equality constraints", Proc. 1975 Int ' l Conf. Cybernetics and Society, 450-452, September, 1975. 23. A.V. Fiacco and G.P. McCormick, Nonlinear Programming: Sequential Un- constrained Minimization Techniques, New York, Wiley, 1968. 24. D.G. Luenberger, Introduction to Linear and Nonlinear Programming, Reading, Mass.: Addison-Wesley, 1973. 25. M.R. Hestenes, "Multiplier and gradient methods", J. Optim. Th. Applic., 4_, 303-320, 1969. 26. M.J.D. Powell, "A method for nonlinear constraints i n minimization prob-lems", in Optimization, ed. by R. Fletcher, New York: Academic Press, 283-297, 1969. 27. D.P. Bertsekas, "Combined primal-dual and penalty methods for constrained minimization", SIAM J. Control, 13, 521-544, 1975. : 28. R.T. Rockafellar, "A dual approach to solving nonlinear programming prob-lems by unconstrained optimization". Math. Programming, 5_, 354-373, 1973. 29. R.T. Rockafellar, "Augmented Lagrange multiplier functions and duality in nonconvex programming", SIAM J. Control, 12, 268-285, 1974. 30. O.L. Mangasarian, "Unconstrained Lagrangians in nonlinear programming", SIAM J. Control, 13, 772-791, July 1975. 31. H. Nakayama, H. Sayama, and Y. Sawaragi, "A generalized Lagrangian fun-ction and multiplier method", J. Optim. Th. A p p l i c , 17, 211-227, 1975. 32. H. Nakayama, H. Sayama, and Y. Sawaragi, "Multiplier method and optimal control problems with terminal state constraints", Int. J. Sys. Sci. , 6_, 465-477, 1975. 133. 33. A.M. Geoffrion, "Duality in nonlinear programming: a simplified appli-cations-oriented development", SIAM Review, 13, No. 1, 1-37, 1971. 34. V.T. Polyak and N.V. Tret'yakov, "The method of penalty estimates for conditional extremum problems", USSR Computational Math and Math Phys. 13, 42-58, 1973. 35. D.P. Bertsekas, "On penalty and multiplier methods for constrained mini-mization", SIAM J. Control and Optim., 14, 216-235, 1976. 36. M. Avriel and V. Gurovich, "A condensation algorithm for a class of algebraic programs", submitted to Operations Research. 37. P. Wolfe, "Methods of nonlinear programming", in Recent Advances in Mathematical Programming, ed. by R.L. Granes and P. Wolfe, New York: McGraw-Hill, 1963. 38. J. Abadie and J. Carpentier, "Generalization of the Wolfe reduced grad-ient method to the case of nonlinear constraints", in Optimization, ed. by R. Fletcher, New York: Academic Press, 1969. 39. A.R. Colville,"A comparative study on nonlinear programming codes", Report No. 320-2949, IBM Scientific Center, 1968. 40. J.E. Kelley, "The cutting plane method for solving convex programs", SIAM J. Appl. Math., 8 703-712, 1960. 41. A.M. Geoffrion, "Elements of large-scale mathematical programming. Part I: concepts", Man. Sci., 14, 652-675, 1970. 42. M. Simonnard, Linear Programming Transl. by W.S. Jewell, Englewood C l i f f s , N.J. : Prentice-Hall, 1966. 43. A.A. Goldstein, Constructive Real Analysis, New York: Harper, 1967. 44. L. Armijo, "Minimization of functions having continuous p a r t i a l deri-vatives" , Pa£i^ic_2i_MathJL, 16, 1-3, 1966. 45. Y. Bard, "Comparison of gradient methods for the solution of nonlinear parameter estimation problems", SIAM J. Numer. Anal., J7, 157-186, 1970. 46. W.C. Davidon, "Variable metric method for minimization", AEC Research and Development Report ANL-5990, Nov. 1959. 47. M. Avriel, Nonlinear Programming: Analysis and Methods, Prentice-Hall, in press. 48. M. Avriel and A.C. Williams, "An extension of geometric programming with applications in engineering optimization", J. of Eng. Math, 5_, 187-194, July 1971. 49. R.L. Francis and J.M. Goldstein, "Location theory: a selective b i b l i -ography", Op. Res., 22, 400-410, 1974. 134. 50. A.P. Hurter, Jr., M.K. Schaefer, and R.E. Wendell, "Solutions of con-strained location problems", Man. Sci., 22, 51-56, 1975. 51. V. Gurovich, "Extensions of signomial programming: numerical solution", M.Sc. Thesis, Technion, Haifa, Israel, 1974. 52. M.K. Schaefer and A.P. Hurter, Jr., "An algorithm for the solution of a location problem with metric constraints", Naval Research Logistic Quart., 21, 625-636, 1974. 53. R.F. Love and J.G. Morris, "Solving constrained m u l t i - f a c i l i t y location problems involving &p distances using convex programming", Op. Res., 23, 581-587, 1975. 54. H.W. Kuhn and R.E. Kuenne, "An efficient algorithm for the numerical solution of the generalized Weber problem in spatial economics, J. of Regional Science, 4, 21-33, 1962. 55. R.N. Sauer, A.R. C o l v i l l e , Jr., and C.W. Burwick, "Computer points the way to more profits", Hydrocarbon Process. Petrol. Refiner., 43, 84-92, 1964. 56. R.E. Payne, "Alkylation - what you should know about this process", Petrol. Refiner. , 3_7, 316-329, 1958. 57. J. Bracken and G.P. McCormick, Selected Applications in Nonlinear Pro- gramming, New York: Wiley, 1968. 58. J.R. Backhurst and J.H. Harker, Process Plant Design, London, U.K.: Heinemann Education Books, 1973. 59. R.A. Greenkorn and D.P. Kessler, Transfer Operations, New York: McGraw-H i l l , 1974.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Signomial programs with equality constraints : numerical...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Signomial programs with equality constraints : numerical solution and applications Yan, James 1976
pdf
Page Metadata
Item Metadata
Title | Signomial programs with equality constraints : numerical solution and applications |
Creator |
Yan, James |
Date Issued | 1976 |
Description | Many design problems from diverse engineering disciplines can be formulated as signomial programs with both equality and inequality constraints. However, existing computational methods of signomial programming are applicable to programs with inequality constraints only. In this thesis an algorithm is proposed and implemented to solve signomial programs with mixed inequality and equality constraints. The algorithm requires no manipulation of the constraints by the user and is compatible with the results of signomial programming. The proposed algorithm is a synthesis shaped by three concepts: the method of multipliers viewed as a primal-dual method, partial dualization, and the retention of the structure of a signomial program. The strategy of the algorithm is to replace the original problem with a sequence of subproblems each subject to the original signomial inequality constraints only. The subproblem's objective function is the original cost augmented with the equality constraints and specified by a parameter vector λ and a penalty constant K. The algorithm then alternates between solving the subproblem and updating λ and K. The convergence of the algorithm is, under suitable assumptions, guaranteed by the convergence of the method of multipliers and the method used to solve the subproblem. Because the subproblem has the form of a regular signomial program, it can in principle be solved by the techniques of signomial programming. A new numerical method implementing a variant of the Avriel-Williams algorithm for signomial programs is suggested for solving the subproblem. The method relies on monomial condensation, and the combined use of the reduced gradient method and a cutting plane scheme. Details of the method's implementation are considered, and some computational experience has been acquired. The proposed method also has the advantageous flexibility of being able to handle non-signomial differentiable objective functions. Four updating schemes for λ and K are formulated and evaluated in a series of numerical experiments. In terms of the rate of convergence, the most promising scheme tested is the use of the Hestenes-Powell rule for updating λ and the moderate monotonic increase of K after the completion of each subproblem. Convergence can also be considerably accelerated by properly scaling the equality constraints and performing only inexact minimization in the first few subproblems. The applicability of the algorithms developed in this thesis is illustrated with the solution of three design examples drawn from structural design and chemical process engineering. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0065625 |
URI | http://hdl.handle.net/2429/20802 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1977_A1 Y35.pdf [ 6.98MB ]
- Metadata
- JSON: 831-1.0065625.json
- JSON-LD: 831-1.0065625-ld.json
- RDF/XML (Pretty): 831-1.0065625-rdf.xml
- RDF/JSON: 831-1.0065625-rdf.json
- Turtle: 831-1.0065625-turtle.txt
- N-Triples: 831-1.0065625-rdf-ntriples.txt
- Original Record: 831-1.0065625-source.json
- Full Text
- 831-1.0065625-fulltext.txt
- Citation
- 831-1.0065625.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065625/manifest