C.I A COMPARISON OF SOLUTION METHODS FOR THE CHEMICAL EQUILIBRIUM PROBLEM by MARGARITA M. RUDA Li e . en Quimica, Universidad de Buenos Aires,Argentina, 1973 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in .THE FACULTY OF GRADUATE STUDIES (Department of Chemical Engineering) We accept t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 22 A p r i l 1982 Q Margarita M. Ruda, 1982 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. I t i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of CH4EM l&A-L. £A>G I MEE & /A/£ The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date ^OGUc.T v 5 / IQPZ DE-6 (3/81) i i ABSTRACT This thesis deals with computing the equilibrium composition of a multiple species reacting mixture. When pressure and temperature are constant and the system i s ideal, t h i s i s the chemical equilibrium problem. It i s possible to approach this problem as the minimization of a non-linear objective function subject to linear equality constraints. The objective function represents the Gibbs' free energy of the system; the constraints refer to the conservation of the elements. Such a formulation corresponds to a "dual geometric program" which i s related to another optimization problem known as the "primal geometric program". In those chemical equilibrium problems with many species, the "primal geometric programming" formulation includes less variables and constraints ( inequality ones) than the dual formulation. We f i r s t compared the primal and dual formulations of the chemical equilibrium problem. Both formulations were solved with a Generalized Reduced Gradient code on seven examples. The primal formulation proved to be 30% faster than the dual for middle-sized problems (up to six simultaneous reactions). The code f a i l e d when trying to solve a dual problem of 24 species and 4 elements.; but thi s same problem was ea s i l y solved when formulated as a primal geometric program. As the geometric programming theory includes s e n s i t i v i t y analysis, and we . were, also interested in the e f f e c t s of small changes of pressure and temperature on the optimal solution, we compared s e n s i t i v i t y analysis with re-optimization of the p r o b l e m . S e n s i t i v i t y a n a l y s i s p r o v e d t o be b e t w e e n 30% t o 50% f a s t e r t h a n r e - o p t i m i z a t i o n . I t a l s o y i e l d e d a c c u r a t e r e s u l t s f o r t h e more a b u n d a n t s p e c i e s when r e l a t i v e l y s m a l l c h a n g e s o f t e m p e r a t u r e a n d p r e s s u r e were o p e r a t e d . H o w e v e r , t h e e q u i l i b r i u m c o n c e n t r a t i o n s o f t r a c e s p e c i e s h a r d l y m a t c h e d t h o s e c a l c u l a t e d by r e - o p t i m i z a t i o n . F rom t h e s e r e s u l t s we recommend t h e u s e o f t h e p r i m a l f o r m u l a t i o n a n d o f r e - o p t i m i z a t i o n t o s o l v e t h e c h e m i c a l e q u i l i b r i u m p r o b l e m , a n d we p r e s e n t a c o m p u t e r c o d e t h a t h a s b e e n t e s t e d on a v a r i e t y o f e x a m p l e s . i v TABLE OF CONTENTS ABSTRACT i i LIST OF TABLES v i i i LIST OF FIGURES x ACKNOWLEDGEMENTS xi Chapter I: Introduction 1 Chapter I I : Mathematical formulations of the problem 5 Thermodynamic relations 5 Formulations of the chemical equilibrium problem 10 The free energy minimization method 10 Reducing the free energy minimization problem 13 Formulation of the equilibrium constant method 16 Geometric programming theory 19 Analogies between GP and the chemical equilibrium problem 22 The reduced dual 26 Summary of the mathematical formulations 26 S e n s i t i v i t y analysis 26 Chapter I I I : Literature review 33 General aspects of the problem 33 The RAND method 34 The NASA method 35 Algorithms to solve geometric programs 36 Algorithms that solve the primal 37 Algorithms that solve the dual 38 V General purpose non-linear optimization algorithms .. 40 Computational experience on GP codes and general non-linear codes for GP 42 Geometric programmming and the chemical equilibrium problem 47 Chapter IV: Writing the program--Preliminaries 49 Mathematical formulations 51 GRG. Description of the code 53 Stating the problem with GRG 53 Description of the algorithm 54 Use of GRG 56 The scaling problem 57 Attempts to solve the scaling problem 58 Comparison between scaling the primal and using the transformed primal 61 Determination of working parameters for GRG 63 A dual starting point 63 Calculation of primal variables from the dual ones 67 S e n s i t i v i t y analysis 69 Calculation of stoichiometric c o e f f i c i e n t s 70 Comparison between the two methods of s e n s i t i v i t y analysis . 71 Conclusions of this chapter 74 Chapter V: Description of the programs 76 Diagrams 76 Input needed. Examples 81 Example 1. CH4-Water gas reaction 82 Example 2. Claus Furnace reaction 85 Calculation of free energy c o e f f i c i e n t s 88 Subroutine LIPSU2 . A dual st a r t i n g point 89 Subroutine SINGV 89 Using GRG as a subroutine 90 Giving input to GRG 90 Subroutine GCOMP 91 Output from GRG 92 Performing s e n s i t i v i t y analysis 93 Output from the programs 93 Chapter VI: Results and Discussion 95 Evaluation of the starting point routine 95 Primal-dual comparisons 98 S e n s i t i v i t y analysis and re-optimization 101 Effe c t s of incorporating new species to a chemical equilibrium model 107 Chapter VII: Conclusions and Recommendations 113 Conclusions 113 Recommendations 115 References c i t e d 116 Nomenclature 122 Appendix A: Computer Programs 124 Program JOTA 125 Program NPLUSK 127 Program COMP1 129 Program COMP2 , . 137 Program COMP3 ' 1 47 Appendix B: examples 153 v i i i LIST OF TABLES Table 11 — 1. GP and chemical equilibrium nomenclature 27 Table II-2. Mathematical formulations of the chemical equilibrium problem 28 Table 111 — 1 . Algorithms compared by Rijckaert and Martens . 43 Table III-2. Algorithms compared by Dembo 44 Table 111-3. Algorithms compared by Gochet, Loute and Solow 45 Table IV-1. Some examples of the chemical equilibrium problem 59 Table IV-2. Scaled primal vs. In-transformed primal with GRG 1 62 Table IV-3. Working parameters for GRG 64 Table IV-4. A comparison of two methods of s e n s i t i v i t y analysis 73 Table VI-1. Starting points calculated within the programs 97 Table VI-2. Time comparison between primal and dual formulations 98 Table VI-3. Time comparison between s e n s i t i v i t y analysis and re-optimization 102 Table VI-4. S e n s i t i v i t y analysis. Effects of the number of increments on the composition 106 Table VI-5. Adding new species to a model effects on the objective function and in the t o t a l number of moles ....108 Table VI-6. Adding new species to a model: effects on ix execution time 109 Table V I-7. Adding new species to a model: effects on the composition 112 X LIST OF FIGURES Figure 1. % Error vs. Variation of temperature for problem 6 104 Figure 2. Composition vs. Temperature for problem 6 110 Figure 3. Composition vs. Temperature for problem 7 111 xi i ACKNOWLEDGEMENTS I would l i k e to thank my advisor, Dr D. W. ' Thompson for his patient help and understanding. ' ; ; •• I am very grateful to Mike Patterson and the rest of the people of the Numerical Analysis group in the Computer Centre at UBC. They were very helpful with programming hints. I also want to thank my husband, Luis Sancholuz, for his encouragement and his typing. Grants from UBC-NAHS are g r a t e f u l l y aknowledged. 1 CHAPTER I: INTRODUCTION The problem of finding the chemical composition of a reacting mixture when i t reaches equilibrium has been c a l l e d the "chemical equilibrium problem" (Dantzig et a l . , 1958). In thi s problem, two thermodynamic variables, pressure and temperature, are fixed. Throughout th i s thesis,we w i l l assume that the system i s id e a l . Knowing the composition at equilibrium is important in applied f i e l d s as diverse as chemical reactors design, development of rocket propellants, evaluation of explosives or biochemistry. When few reactions take place,the calculations needed are quite t r i v i a l . They become very cumbersome, however, when there is a large number of chemical species and simultaneous reactions . If these calculations have to be repeated for s l i g h t changes in the fixed conditions ( i . e . , pressure and temperature), the problem can be quite impossible to solve in a reasonable period of time. A d i g i t a l computer i s then required. This thesis aims at producing a general computer program to solve the chemical equilibrium problem for ideal systems of many chemical species. To avoid re-solving the problem when small changes in pressure and temperature need to be considered, the program should also be capable of performing s e n s i t i v i t y analysis. Eventually, t h i s program w i l l be used for simulation and design of chemical reactors, mainly for gaseous reactions. To achieve our goal, we needed to address the following 2 points: i) the mathematical formulation of the problem, i i ) t h e numerical method of solution, and i i i ) the s i m p l i c i t y of the program for the user. In that order we now proceed. The mathematical formulation of the problem has to be general enough to apply to a wide range of examples. The chemical equilibrium problem has been posed in many ways, a l l of them mathematically equivalent a l b e i t numerically d i f f e r e n t . Van Zeggeren and Storey (1970) c l a s s i f i e d these formulations as the "free energy minimization" approach and the "equilibrium constant" method. The former i s an optimization problem with a non-linear objective function and linear equality constraints. The l a t t e r can be derived from the condition of minimum Gibbs free energy of the system, but consists of a set of simultaneous non-linear equations to be solved. Passy and Wilde (1968) found a relationship between the chemical equilibrium problem and Geometric Programming, a mathematical programming technique (GP hereinafter). The "free energy minimization approach" can be regarded as a "dual GP" problem to which corresponds a related "primal GP" formulation. The "primal GP" formulation has a smaller dimension than the "dual GP", and has been succesfully used to solve examples of the chemical equilibrium problem (Passy and Wilde, 1970; Dembo, 1976; Lidor, 1975). The numerical method of solution i s of course related to the mathematical formulation of the problem; however, the same formulation can be solved using d i f f e r e n t numerical methods. The selection of the numerical method of solution i s important 3 for the computational effectiveness of our code. In other words, we want to produce accurate results fa s t . Both the primal GP and dual GP formulations present numerical d i f f i c u l t i e s . There i s not enough computer evidence in the l i t e r a t u r e as to which method of solving the chemical equilibrium problem performs better. In t h i s work we s h a l l compare the perfomances of the primal and the dual GP formulations when solving the chemical equilibrium problem. A general purpose non-linear optimization code, which can handle both formulations, seems appropiate to make that comparison possible. GP theory provides ways of performing s e n s i t i v i t y analysis. This is p a r t i c u l a r l y useful to us since we are interested in comparing i t s speed and accuracy with that of actually re-solving the problem. The results of the comparisons (primal vs. dual and s e n s i t i v i t y analysis vs. repeating the method) w i l l permit us to choose the more appropiate mathematical formulation and numerical method of solution for a f i n a l computer code to solve the chemical equilibrium problem. We f i n a l l y come to the question of s i m p l i c i t y . We want to avoid the user's tendency to err while inputing data. The amount of information needed to run the program should be reduced to a minimum. A subroutine to calculate a f i r s t s t a r t i n g point for the optimization i s therefore necessary. The remainder of the thesis is organized as follows. Chapter II presents the d i f f e r e n t mathematical formulations of the chemical equilibium problem . Chapter III reviews and discusses 4 the l i t e r a t u r e on the subject. Chapter IV deals with the d i f f i c u l t i e s found in trying to set up the programs, while Chapter V describes the programs actually written for the proposed comparisons. Chapter VI presents and discusses the results obtained with our programs in a series of examples taken from the l i t e r a t u r e . F i n a l l y , conclusions and recommendations are included in Chapter VII. 5 CHAPTER I I : MATHEMATICAL FORMULATIONS OF THE PROBLEM This chapter w i l l serve as a means of defining the problem, as well as introducing the nomenclature. It i s based on a l i t e r a t u r e review, and i t i s organized as follows. Section 1 deals with the thermodynamic relations that describe the chemical equilibrium problem. In section 2 we introduce the two t r a d i t i o n a l mathematical formulations known as the "Gibbs free energy minimization method" and the "equilibrium constant method" (Van Zeggeren and Storey,1970). In section 3 we present the geometric programming theory, and we discuss the relati o n s h i p between th i s mathematical programming technique and the chemical equilibrium problem. In section 4 we discuss the basis for s e n s i t i v i t y analysis in geometric programming. Thermodynamic relations The thermodynamic basis of the "chemical equilibrium problem" has been extensively discussed in many textbooks (Denbigh, 1966, Kirkwood and Oppenheim, 1961). We w i l l review the thermodynamic relations that w i l l allow us to formulate mathematically our problem. When temperature (T) and pressure (P) are chosen as independent variables for a thermodynamic system, the appropriate fundamental r e l a t i o n that completely describes the system i s the one expressing the Gibbs free energy (G) in terms of P,' T, and the composition variables. We w i l l assume a system 6 of K phases and N c h e m i c a l s p e c i e s ; each c h e m i c a l s p e c i e s i s p o t e n t i a l l y p r e s e n t i n each phase. Then, G = G(T,P, 6j k ) 2.1 where 6 j = number of moles of s p e c i e s j i n phase k S i n c e Gibbs f r e e energy i s an e x t e n s i v e p r o p e r t y , we can r e w r i t e eq. 2.1 as : N E j = 1 G - 6 j k 2 • 2 where the c h e m i c a l p o t e n t i a l i s an i n t e n s i v e p r o p e r t y d e f i n e d as f o l l o w s : M j k - ) 2 - 3 O 6 j K T , P , 6 L K I f we n e g l e c t a l l i n t e r a c t i o n s among the phases, each phase c o n t r i b u t e s a d d i t i v e l y t o the Gibb s f r e e energy of the system. Then, 7 K G = E Gk ( T , P , 6 j k ) 2.4 k=1 Where Gk = Gibbs free energy of phase k Let us now take a look at the functional form of the chemical p o t e n t i a l . The chemical potential of the species j in the phase k can be written in the following form: (Denbigh, 1966) „ J l t (T,P,6 J l c ) = „ J k 0 (T,P)+RT In a j < 2.5 where M j 1 < : 0 (T,P)= reference value for the chemical potential of the species j in the phase k. R = universal gas constant a j K = a c t i v i t y of species j . We w i l l now assume that each phase is an ideal mixture. Moreover, we w i l l consider each species as belonging to one phase; i f one species belongs to more than one phase, i t w i l l be given a d i f f e r e n t number. Hence, k i s no longer needed as a subscript for the number of moles. In real l i f e , our problem i s r e s t r i c t e d to : a) gaseous reactions, b) reactions in pure condensed phases, c) some b i o l o g i c a l models. 8 With these assumptions, the a c t i v i t y of one species in the gaseous phase is equal to the molar fraction of t h i s species in the gas phase, times the t o t a l pressure. For the ideal condensed phases, the a c t i v i t y equals the molar f r a c t i o n . In equations, i f we define = I jek 2.6 Then the molar fraction of species j , Xj i s : Xj = 6 , / X Y 2.7 and, for an ideal gas phase = Xj P/1 atm. 2.8 for pure condensed phases a, = Xj 2.9 we can rewrite equation 2.5 for the gas and for the pure condensed phases 9 u j (T,P) = K J 0 (T) + RT In P + RT In (6j / X k ) 2.10 V J (T,P) = „j 0 (T,P)+ RT In (6j / U K ) 2.11 where the superscript 0 refers to standard state. For the gas, i t is the chemical potential of the species j as an ideal gas, at zero pressure. For the condensed phases, i t i s usually the chemical potential of the pure species j in the condensed phase at the same T and P. We w i l l now define the free energy c o e f f i c i e n t s Cj as: Cj = ( MJ° (T) /RT) + In P/1 atm. 2.12 for the gas phase, and Cj = Mj° (T,P) / RT 2.13 for the condensed phases. Replacing equations 2.10, 2.11, into 2.2, dividing by R T , and using 2.12 and 2.13, we get a convenient formulation of the Gibbs free energy for an ideal system of N species and K ideal phases, at T and P constant. N G/RT = E 6j j = 1 (In 6j + Cj) K I Xt; In X.« k=1 2.14 We have now a working equation that describes Gibbs free 10 energy in terms of the number of moles of the species present at equilibrium, the t o t a l number of moles per phase, the temperature, pressure, and the "free energy c o e f f i c i e n t s " Cj Let us consider now the equilibrium conditions. The condition for equilibrium in a closed system described by Gibbs free energy (G) i s that G i s a minimum. For the ideal case of equation 2.14, this minimum exists and is unique. (Denbigh, 1966). The fact that G is a minimum at equilibrium implies that the variations of G produced by independent variations must be zero. But not a l l the variations in the number of moles are independent. They must s a t i s f y the requirement that the t o t a l mass of each element is d i s t r i b u t e d among the d i f f e r e n t chemical species (Zeleznik and Gordon, 1968). We need a mathematical description of these constraints. Formulations of the chemical equilibrium problem The free energy minimization method The "free energy minimization " method i s just the mathematical formulation of the equilibrium conditions stated above. We s t i l l have to formulate the constraints, since the Gibbs free energy is described by equation 2.14. To do that, we define the "exponent matrix" AA . It is an M x N matrix, where M i s the number of elements in the system, and N is the number of chemical species in the system. Each column of the matrix corresponds to the chemical formula of one species: AAi,j gram atoms of the element i are present in one mole of species j . 11 If the system contains Bi gram atoms of element i , then the conservation of elements can be written as: K E E AAi j 6j = Bi i = 1 ,M 2 .15 k=1 jek If ionization i s considered, the conservation of charge can be also expressed as in equation 2.15 ; charge i s considered to be the M+1 element, with zero amounts; the corresponding row in the exponent matrix i s the charge of each species. Besides the conservative constraints, we have N p o s i t i v i t y constraints, since the number of moles of a chemical species i s either p o s i t i v e or zero. They are: 6 j ^ 0 j = 1,N 2.16 Combining the equations 2.14, 2.15, and 2.16, we have the "Gibbs free energy minimization " formulation of the chemical e q u i l i l b r i u m problem. We w i l l c a l l this formulation Problem A. •v. 1 2 Problem A Min G/RT N E K 6 j ( I n 6 j + Cj) - E X. K I n X t< j = 1 k=1 s. to: K E E k=1 jek AAij 6j = Bi i = 1 , M > 0 j =1 ,N 2.17 In short , Problem A i s an optimization problem of: a) non linear objective function of N variables (the number of chemical species present at equilibrium). b) M linear equality constraints (number of elements). c) N p o s i t i v i t y conditions. In spite of Gibbs developing his theory in the past century, i t was not u n t i l recently that Problem A could be solved e f f i c i e n t l y by computer optimization techniques. Among the d i f f i c u l t i e s , the constraints are equality ones; the objective function i s convex for ideal systems, but i t is non-di f f e r e n t iable i f one of the species i s zero, and i t i s not defined in such a case; the dimension of the problem is the number of species present at equilibrium, which can be a large number. The main advantage of Problem A is i t s s i m p l i c i t y of formulation; the only data needed are the exponent matrix , the amount of elements (B vector) and the free energy c o e f f i c i e n t s C, determined by the working pressure and temperature. 13 Reducing the free energy minimization problem Problem A was too large to be solved when computers were not available; stoichiometry was used to simplify numerically the problem. The idea was to write some species (constituents) as a function of others (components) (Brinkley, 1946). A system of many constituents was regarded as being formed by the components through a set of simultaneous l i n e a r l y independent chemical reactions. Just how many of these reactions should be considered, or which should be the components, has been the object of much research for the case of many constituents (Denbigh, 1966; Schubert and Hoffmann, 1976; Waller et a l • 1980). A well posed chemical equilibrium problem should have a l l i t s proposed species greater than zero at equilibrium; also the balances (eq. 2.15) should form a set of linear equations with rank = M. If this i s the case, there exists a matrix U c a l l e d the stoichiometric matrix, with dimensions N x D , and rank D, where D= N-M , is the number of independent chemical reactions. The stoichiometric matrix i s such that, in matrix notation, AA.U = 0 2.18 Then the composition of each species can be written as a linear combination of D independent parameters r c a l l e d the extent of the reactions. That i s , 1 4 6j (r) D = 6j° + E Ujd r d j = 1 ,N 2.19 d=1 where 6j 0 = i n i t i a l amount of species j Also see that: d 6 j ) = Ujd T,P 2.20 which i s the t r a d i t i o n a l way of defining the extent of reaction. To calculate the stoichiometric matrix i s easy when few reactions take place but i t is not so for the case of many reactions. Some systematic ways of c a l c u l a t i n g i t from the mass balances have been developed l a t e l y (Schneider & R e k l a i t i s , 1975). In fact, any solution of the homogeneous system of equations formed by setting the B vector equal to zero in equation 2.15, is a stoichiometric matrix , even i f i t does not look so nice to a chemist. (Aris, 1970). The selection of the components is the more d i f f i c u l t task in a l l the systematic procedures. The best way to do i t is to choose as components those species present in greater quantities at equilibrium, which is most of the time d i f f i c u l t to know beforehand. Isomers also present d i f f i c u l t i e s ( C a v a l l o t t i et a l . , 1980). If the number of moles of the species are written as in equation 2.18, then the free energy of the system can be stated as a reduced problem of D variables, (the extents of each reaction ) with D p o s i t i v i t y constraints. To simplify the 15 expression resulting by replacing 2.18 in 2.14 and 2.15, we introduce the following constant terms: N E 6f Cj j = 1 •In Ko = 2.21 R T N E Ujd Cj j = 1 •In Kd = 2.22 R T We now define: X k(r) = X k° + E rd X K d k = 1,K d=1 2.23 ° " E 6 J ° j C k k = 1 ,K 2.24 X KJ = E Ujd j€k d = 1,D k = 1,K 2.25 After some algebra we get the reduced free energy minimization formulation, which we w i l l c a l l Problem B. 16 Problem B D Min G/RT = - In Ko + I r6 ln Kd + d=1 .:-N K E (r) ln 6 j ( r ) - E X. K(r)ln X. K( r) j = 1 k=1 s. to D (r) = 6 j 0 + E r j Ujd > 0 d=1 j : = 1 ,N Problem B is an optimization problem, with a) non-linear objective function of D variables, which are the extents of each independent reaction. b) N p o s i t i v i t y conditions, one per species present at equilibrium. In t h i s formulation, the number of variables i s smaller than in the previous one; but when N>2M , this advantage is not so important. The objective function i s s t i l l undefined and non d i f f e r e n t i a b l e when a variable is equal to zero. The set of stoichiometric c o e f f i c i e n t s has to be calculated before the optimization procedure. Formulation of the equilibrium constant method Problems A and B are two forms of the optimization approach to the chemical equilibrium problem. In order to derive the "equilibrium constant" approach, we w i l l consider Problem B and the conditions for optimality. If G is a minimum at 1 7 equilibrium, with pressure and temperature constant, then the gradient of G has to be zero at t h i s point. Taking the gradient of G, expressed as a function of the extents of reactions, and setting i t equal to zero, we get a set of D non-linear equations in D unknowns. These equations are the well known "mass action law". "Problem C" below states the equations. Problem C for gaseous phase • • Kd = exp E Ujd Cj jek RT n jek ( 6j° + D E d=1 Ujd r d ) P Ujd E (6,° jek D + E d=1 Ujd rd ) k = 1 d=1 ,D 18 Problem C for condensed phases K = exp E Ujd Cj jek RT n j€k k = 2 , K d = 1 ,D D ( 6 j ° + E Ujd r 6 ) d=1 D E (6j 0 + E Ujd r d ) jek d=1 Ujd 2 . 2 8 Stoichiometric c o e f f i c i e n t s are needed in Problem C. Once the set of reaction variables r that s a t i s f i e s the equations i s known, the number of moles present at equilibrium i s calculated through equations 2 . 1 9 . Problem C involves solving a system of non linear equations. This is a t o t a l l y d i f f e r e n t numerical approach than problems A and B, although the mathematical formulations have been shown to be equivalent. Almost a l l formulations of the chemical equilibrium problem in the l i t e r a t u r e f a l l into Problems A, B, or C, with some algebraic modifications. There i s s t i l l a d i f f e r e n t formulation, and to introduce i t we need some background on the mathematical programming technique known as Geometric Programming. 19 Geometric programming theory Geometric Programming is the mathematical formulation of a special kind of optimization problem. It involves the minimization of a posynomial (positive polynomial) objective function g , subject to inequality constraints that are also posynomials. The problem described i s c a l l e d a "primal problem" ; the primal function g belongs to the Euclidean space R+m there exists a related "dual" maximization problem, involving a function v that belongs to the dual space R+n. The function v has the form of a product of non-linear terms, and has M+1 linear constraints. It has been found that the constrained maximum value of v is equal to the constrained minimum value of g (Duff i n , Petersen and Zener, 1966). We w i l l now present the equations that exemplify a l l t h i s . A Primal Geometric Program is of the form: Pr imal GP Min g„ (t) s. to g K (t) < 1 t i >o k=1 ,K i = 1 ,M where 9K.(t) = I C j jek M n t L Ai j i = 1 k = 0,K j = 1 ,N The corresponding maximizing Dual Geometric Program (Dual GP) i s 20 as follows: Dual GP N Max v(6) = n j = 0 (cj / 6j K ) n k = 1 s. to E 6 J = 1 jeO normality condition N E A i j 6J j = 0 = 0 i = 1 ,M orthogonality condi t ions 6j > 0 j = 0,N p o s i t i v i t y conditions where X.k = E 6 j jek k = 1 ,K The r e l a t i o n between the primal and dual problems i s c a l l e d the "geometric inequality" and states: g (t) > min g(t) = max v(6) ^ v(6) 2.31 When the equality holds at the optimum (*), the primal and the dual variables are related through the following set of equations (Duff'in, Petersen & Zener, 1966) : 21 M In cj + E A i j In t * = In (6j* g c ( t * ) ) jek = 0 2.32 i = 1 M In C j + E A i j In t* = In ( 6 j * / k K* ) jek = 1,N 2.33 i = 1 If in the primal problem the constraints are active, and i f i t is solved by means of a Lagrangian technique, there is a rela t i o n s h i p between the Lagrange m u l t i p l i e r for the constraint and the corresponding t o t a l number of moles of the phase : X k = n K / g o 2.34 where nx = Lagrange m u l t i p l i e r of the k-primal constraint We w i l l now summarize some important aspects: a) Each primal variable is associated to one of the orthogonal conditions in the dual. b) Each primal constraint is associated to one of the X values. c) Each primal term is associated to one dual variable. d) The posynomial terms in the primal objective function correspond to the dual variables subject to the normality condit ion. We w i l l not go any further into Geometric Programming theory for the moment; we w i l l discuss instead the rel a t i o n s h i p between 22 Geometric Programming and the chemical equilibrium problem. Analogies between GP and the chemical equilibrium problem In 1968 Passy and Wilde found that the chemical equilibrium problem, stated as our Problem A , could be regarded as a dual geometric program. Problem A was a l g e b r a i c a l l y transformed through the following d e f i n i t i o n s to f i t into the standard formulation of a dual geometric program. 60 = 1 c j = exp(-Cj) c 0 = 1 -Bi = Aio Since exp(-G/RT) is a monotonic function, finding i t s maximum is equivalent to minimizing the negative of the function; taking logarithms, we have: Min (G/RT) = Max exp(-.G/RT) 2.39 The Dual GP chemical equilibrium problem is then: 2.35 2.36 2.37 2.38 23 Chemical equilibrium as a dual GP Max v = exp(-G/RT) = N n j = (cj / 6 3 0 6 j ) K kK n X f e k=1 s. to 60 = 1 normality condit ion N Z A i j 6j = j = 0 = 0 i = 1 ,M orthogonality conditions 6j >- 0 j = 0,N p o s i t i v i t y conditions where X K = Z 6j jek k = 1 ,K If we maximize the logarithm of the previous problem, we obtain a "transformed dual" formulation which is similar to our "Problem A". The matrix A is the augmented exponent matrix (AA); 6 0 is a dummy species. 24 Transformed dual GP Max In v N = Min G/RT = E j=o K I X K In X« k=l (In 6j + Cj) -s. to 6 0 = 1 normality condition N E A i j 6J = 0 j = 0 i = 1 ,M orthogonality conditions where 6j ^ 0 X < = E 6J jek j k = 0,N = 1 ,K p o s i t i v i t y conditions From the theory of Geometric Programming, to the previous dual problem corresponds the following primal problem: Chemical equilibrium as primal GP M Aio Min go(t) = n ti i=1 s. to g * ( t ) = E c-jek M A i j n t i < 1 k = 1 ,K i = 1 ti ^ 0 i = 1 ,M where Aio=-Bi 25 To the "transformed dual" corresponds a "transformed primal" formulation; we define Z = ln t and then: Transformed primal M Min ln g 0 ( t ) = Min h(Z) = E Aio Zi i = 1 s. to M E exp(Cj + E A i j Zi) < 1 k = = 1,K j ek i = 1 The primal GP chemical equilibrium problem i s a new formulation of the problem, with: a) minimization of a non-linear objective function of M variables (the number of elements) b) K non-linear inequality constraints, one per phase. The constraints are active at the optimum. Each term of the constraints is equal to the molar fr a c t i o n of the corresponding species in the phase; hence the constraints state that the summation of the molar fractions per phase has to be one at equilibrium. The t o t a l number of moles per phase can be calculated from eqn. 2.34 The transformed primal involves: a) minimization of a linear objective function of M variables b) K non-linear inequality constraints, with the same physical meaning as in the primal problem. If the transformed primal i s solved using the Lagrange m u l t i p l i e r s technique, each Lagrange 26 m u l t i p l i e r i s equal to the t o t a l number of moles in that phase. The reduced dual The same methods described before to derive Problem B from Problem A are used in the Geometric Programming theory. Problem B i s thus equivalent to a "Reduced Dual Geometric Program" , except for some changes in nomenclature. These changes are shown in Table 11 — 1. Summary of the mathematical formulations Table II-2 summarizes the most important features of the mathematical formulations of the chemical equilibrium problem. We w i l l now take a look at one convenient derivation from the Geometric Programming theory : the post-optimal analysis known as s e n s i t i v i t y analysis. S e n s i t i v i t y analysis A l l the formulations of the chemical equilibrium problem are depending on the free energy c o e f f i c i e n t s C, the exponent matrix AA and the vector of the quantity of the elements, B. The exponent mateix w i l l not change i f the model i s well formulated. The free energy c o e f f i c i e n t s w i l l probably vary with the source from where they are obtained, but more important for p r a c t i c a l purposes is their v a r i a t i o n with the temperature and pressure of the system. The B vector may vary with d i f f e r e n t i n i t i a l compositions, for example. 27 Table 11 — 1 . A p a r a l l e l between Geometric Programming and Chemical Equilibrium nomenclatures. Symbol Chemical Nomenclature GP Nomenclature r 6 k = 0 6o Ujd Kd Extent of reaction Number of moles of species j Not defined Not defined Stoichiometric coef f ic ients Equilibrium constant for the d reaction Basic variable Dual variable Corresponds to primal objective funct ion Normality condition N u l l i t y vectors Basic constant S e n s i t i v i t y analysis is just a way of evaluating the changes in the equilibrium Gibbs free energy and in the composition, when any of the parameters mentioned above i s changed, without act u a l l y re-solving the optimization problem. The method involves a numerical c a l c u l a t i o n of the p a r t i a l derivatives at the optimum, that i s , the truncation of a Taylor series expanded around an already found optimum. Of course, i f the variations in the parameters are large enough, the problem has to be re-solved . The problem of obtaining a set of s e n s i t i v i t y equations has been approached in two ways in the l i t e r a t u r e . Both approaches 28 Table II-2. A summary of the more important c h a r a c t e r i s t i c s of the formulations of the chemical equilibrium problem. Name of the Eqn. Type of Number of Number and Need formulation numerical variables type of stoich. in t h i s work solution constraints coeff. Problem A Problem B Problem C GP dual and transformed dual GP reduced dual GP primal 2.17 Optimization N 2.26 Optimization 2.27 System of D D 2.28 nonlinear unknowns equat ions 2.40 Optimization N+1 2.41 2.26 Optimization D 2.42 Optimization M GP transformed 2.43 Optimization M primal M linear eq. const. N p o s i t i v i t y conditions D nonlinear iheq. const. No M+1 linear eq. const. N p o s i t i v i t y conditions D nonlinear constraints Yes Yes No Yes K nonlinear No ineq. const. M p o s i t i v i t y conditions K non linear No ineq. const. Eqn.:equation, stoich. coeff.: stoichiometric c o e f f i c i e n t s , . eq.:equality, ineq.: inequality, const.: constraints, N: number of species, M: of elements, K : of phases, D: of reactions. d i f f e r from the numerical point of view. The f i r s t approach was stated in appendix B of Duffin, Petersen and Zener's book (1966). They work with the Reduced Dual GP (eqn. 2.26). The Jacdbian matrix of the transformation 29 from the r variables to the Kd "variables"(matrix J) is formed as follows : N K Jqs(6) = E (Ujq Ujs /6j *) " Z U K q X K S /\K*) 2.44 j=1 k=1 where q,s = 1,D X M = E Ujq 2.45 jek X K S = E Ujs 2.46 jek Xfc* = E 6j * 2.47 jek Duff in et a l . (1966) proved that the matrix J is also the Hessian matrix for the function In v r e l a t i v e to the basic variables r (extents of reactions). So, i f the reduced dual GP equivalent to the Problem B i s solved by any method^involving derivatives, the matrix J i s readily a v a i l a b l e . The matrix J evaluated at the optimum point is used both to introduce c r i t e r i a for the existence of derivatives and to derive expressions for these derivatives. The d i f f e r e n t i a l changes are approximated l i n e a r l y . Dinkel and Lakshmanan (1976, 1977) used these expressions in the case of the chemical equilibrium problem when P, T, and the amount of elements changed. They applied their results to two examples, the problems 3 and 5 of appendix B. They also 30 suggested the use of an incremental procedure in order to control the error produced by the linear approximation of the d i f f e r e n t i a l changes. The expressions they obtained are as follows: The expressions are v a l i d i f J i s non singular at the optimum. As seen from equations 2.48 and 2.49, the inversion of the matrix J evaluated at the optimum i s a necessary step to perform s e n s i t i v i t y analysis with this approach. When the number of reactions D is not very large, t h i s is not a problem; but as the dimension of the matrix increases, so do the rounding errors and the time that is consumed in the inversion. The second approach taken was related to Generalized Geometric Programming. This is an extension of the Geometric Programming theory to "generalized polynomials", that i s polynomials with some negative terms. These researchers also t r i e d to obtain expressions for numerical derivatives, but their approach was s l i g h t l y d i f f e r e n t than the previous one (Rijckaert, 1974) . The M+1 dual constraints were written in e x p l i c i t form, over the N+1 variables , as well as the D Av v* 2.48 2.49 31 "equilibrium conditions" and the K equations r e l a t i n g X to the summation of the number of moles in each phase. A l l t h i s i s equivalent to writing the Kuhn Tucker conditions for the dual GP chemical equilibrium problem. They wrote the N+1+K set of equations for the optimum and for a small perturbation of the optimum. Each equation in the perturbed set was subtracted from the corresponding one in the optimal set, and a Taylor's approximation series was used to l i n e a r i z e the equations. The f i n a l system of N+K+1 linear equations in N+K+1 unknowns( the variations in the number of moles) i s as follows: A 6 0 = 0 2.50 K N E E A i j A6j = 0 k=0 j=0 i = 1 ,M 2.51 E A6 0 " A X K = 0 jek k = 1 ,K 2.52 K K E E (Ujd/6j*) A6 - E (Ujd /X ,< * ) A X K k=0 jek k=0 K N = E E Ujd(ln C J - ln cj*) k=0 j=1 d = 1,D 2.53 For any changes in the free energy c o e f f i c i e n t s , only the 32 right hand side of eq. 2.53 changes. For any changes in the B vector, only the f i r s t column of eq. 2.51 changes. No inversion of a matrix i s required, but the dimensions of the system are much larger than in the f i r s t approach. 33 CHAPTER I I I : LITERATURE REVIEW This review of the l i t e r a t u r e i s organized in three sections. F i r s t we take a look at the l i t e r a t u r e on general methods of solution for the chemical equilibrium problem. In the second section, we focus on the computer codes available for solving geometric programs. The t h i r d section deals with the geometric programming approach to the chemical equilibrium problem. General aspects of the problem In the previous chapter we have presented d i f f e r e n t formulations of the chemical equilibrium problem. They are a l l mathematically equivalent, but di f f e r e n t numerical methods of solution are used for each case. The methods of solution f a l l into two main categories: a) optimization methods b) solution of non-linear systems of equations. The problem now seems to be reduced to a choice of a mathematical formulation and a numerical method of solution. The l i t e r a t u r e on the possible combinations i s very wide; an extensive review was conducted by Van Zeggeren and Storey (1970) and by Zeleznik and Gordon (1968). However, the problem is not yet s e t t l e d ( C a v a l l o t t i et a l . , 1980). Two factors should be pondered when judging the goodness of a solution method: the purpose of the calculations, and the means that are available to carry them on. We are not as much 34 interested in the best solution for one s p e c i f i c case, as we are in a method f l e x i b l e enough to accommodate a wide range of chemical equilibrium examples. The calcu l a t i o n s w i l l be performed by a d i g i t a l computer. At the moment, there are two general computer codes available commercially that solve the chemical equilibrium problem. They belong to the RAND Corporation (1965,1970) and to the NASA (1971) respectively. We w i l l b r i e f l y describe their methods. The RAND method RAND researchers devised the f i r s t version of the RAND program in 1958 (Dantzig, Johnson, White and De Howen, 1958). The program was revised and completed by R.J.Clasen, N. Shapiro, M. Shapley and others over a period of 15 years. The chemical equilibrium problem i s formulated as a minimization of the Gibbs free energy subject to mass balance constraints. The formulation is similar to Problem A of Chapter II in t h i s thesis; i t deals with an ideal gas phase and pure condensed phases. The solution i s obtained in two steps. The " f i r s t order method" provides a f i r s t set of composition values through a si m p l i f i c a t i o n of the problem to a line a r program. The "second order method" uses the previous set of values as a f i r s t guess. The Gibbs free energy i s approximated at t h i s point using a second order Taylor's series. Then, the problem i s transformed to an unconstrained optimization, using Lagrange m u l t i p l i e r s . A f i n a l set of equations i s then solved using the Newton-Raphson 35 technique. The RAND program has proved successful over a number of years in solving the chemical equilibrium problem. The main disadvantages, however, are: a) slow convergence when many trace species are included in the model b) i n a b i l i t y to handle cases where the number of moles of a species i s zero. The NASA method The f i r s t NASA programs to calculate chemical equilibrium composition were based on the "chemical equilibrium constant" formulation. The method was f i r s t devised by Huff (1951) and was later modified by Zeleznik and Gordon (1960,1962). However, in 1971, they changed the method to a free energy minimization procedure;the reasons given for the change included ( NASA , Manual Report,1971): a) more bookkeeping necessary for the equilibrium constant method. b) numerical d i f f i c u l t i e s with the use of components (compared with keeping a l l the variables as such). c) more d i f f i c u l t y in extending the generalized method for non-ideal equations (but s t i l l the program handles only ideal gases and pure condensed phases). It was shown by Gautam and Sei.der (1979) that the RAND and NASA methods,although derived d i f f e r e n t l y , give nearly i d e n t i c a l equations and are both implementations of Newton's method. The 3 6 NASA code calculates thermodynamic derivatives as well as the chemical composition. We have already stated that the chemical equilibrium problem can be regarded as a geometric program. Let us now look at some algorithms devised to solve geometric programming problems. Algorithms to solve geometric programs Many algorithms to solve geometric programs have been proposed in the l i t e r a t u r e . A summary of some of these can be found in Bleightler and P h i l l i p s ' book(l976). We may c l a s s i f y the algorithms in three groups. a) algorithms that solve the primal geometric program (GP), b) algorithms that solve the dual GP, c) general non-linear optimization methods that can be applied to solve either the primal or the dual GP. For the general case, i t i s not possible to predict i f the dual or the primal problem are easier to solve. Empirical evidence based on computational comparison i s needed. Some comparisons among codes were performed by Dembo (1978),Sarma, Martens et a l (1978) and by Gochet, Loute and Solow (1978). We w i l l give a brief description of some of the GP and general purpose optimization codes compared and we w i l l then try to summarize the conclusions of the papers mentioned above. 37 Algorithms that solve the primal The primal GP was stated in equation 2.29 (chapter I I ) . The chemical equilibrium problem as a primal GP was described by equation 2.42. .. ; If we minimize instead the logarithm of the objective function, and i f we make the transformation of variables z = ln t, then the problem becomes a convex program. The t's have to be positive,but the z's are unrestricted in sign. The resulting problem is c a l l e d the transformed primal problem,and was s p e c i f i e d in equation 2.43 for the chemical equilibrium problem. Another way of solving the primal is via separable programming techniques; t h i s formulation increases the dimensionality of the problem from M to N. Based on the previous formulations, the algorithms that solve the primal can be c l a s s i f i e d in three groups. 1) Condensation : the primal geometric programming i s solved d i r e c t l y by condensation or l i n e a r i z a t i o n of posynomial functions (cutting planes algorithms). A well known code based on t h i s approach i s GGP. Written by Dembo (1974), i t is based on a Kelley's cutting plane algorithm. 2) Kuhn-Tucker conditions for optimality. The Kuhn Tucker conditions for primal geometric programming are solved i t e r a t i v e l y , using a condensation technique. One code based on th i s approach i s GPKTC, written by Rijckaert and Martens (1976). This method is e s s e n t i a l l y equivalent to a Newton-Raphson algorithm for the Kuhn-Tucker conditions expressed in terms of the variables z=ln t. The.code FP from Gochet, Loute and Solow, 38 i s based on the same approach 3) Separable Programming. A good code based on thi s formulation is DAP, by G.V.Reklaitis ; he uses the d i f f e r e n t i a b l e algorithm of Wilde and Bleightler (Dembo, 1978). Algorithms that solve the dual The dual geometric programming (DGP) is the l i n e a r l y constrained, non-linear programming problem stated on equation 2.30 (chapter I I ) ; equation 2.40 for the case of chemical equi1ibr ium. Taking logarithms of the objective function, the problem i s transformed to a concave objective function to be maximized, subject to a linear equality constraints. For the chemical equilibrium, that transformation i s seen in equation 2.41. The reduced dual geometric program (RDGP) i s obtained by eliminating M+1 basic variables from the program DGP and expressing them in terms of D= N+1-(M+1) nonbasic variables (see eqn. 2.26, 2.27 and 2.28 on chapter II) most of the algorithms that solve the dual make use of the reduction of dimensionality that the RDGP problem provides. The main problems related to solving the dual are (Gochet, Louter, and Solow, 1974): a) non d i f f e r e n t i a b i l i t y of the objective function with respect to the dual variables where they take on the value zero. b) the dimensionality of the dual problem w i l l always be larger than that of the primal, except for cases with N<2M solved using the reduced dual GP. 39 c) dual variables have to be determined with higher accuracy i f primal variables are to be calculated from them. The algorithms that solve the dual problem are based on the following p r i n c i p l e s : 1) Linear approximation methods The dual program in i t s logarithmic form i s approximated at a s p e c i f i c point by a linear program (LP). The solution of the LP provides a d i r e c t i o n for improving the value of the objective function. The step i s provided by a l i n e minimization. One code that follows t h i s approach i s LAM, from Rijckaert and Martens (1978): the LP routine i s based on the Revised Simplex Method with product form of the inverse. The l i n e minimization is based on cubic interpolation. 2) Separable programming. After a logarithmic transformation, the dual objective function is separable and can be approximated by a piecewise linear function. One code that follows t h i s approach i s SP, by Rijckaert and Martens (1978) 3) Gradient projection methods. This approach combined the gradient projection method due to Rosen with a variable metric method in order to approximate the inverse of the Hessian of the objective function. One example i s the code VMP (Sargent and Murtagh, 1973). 4) Newton-Raphson. The Kuhn-Tucker conditions for the reduced dual geometric program result in a system of non linear equations that are solved using a Newton-Raphson technique (see the s i m i l a r i t i e s with the "equilibrium constant" approach in the chemical equilibrium problem). Many codes are based on t h i s 40 approach, l i k e NEWTGP by J. Bradley (1973) (that code has the cap a b i l i t y of performing s e n s i t i v i t y a n a l y s i s ) , GEOGRAD by Dinkel and Kochenberger (1974), and NRF by Rijckaert and Martens, (1978). 5) Other methods The code GOMTRY by Blau and Wilde (1971) solves the Kuhn Tucker conditions for the separable dual program ; the code CSGP by Beck and Ecker (1978) applies the concave simplex method to the dual; the code MCS introduces a modification to CSGP that allows for blocks of variables to go to zero simultaneously. General purpose non-linear optimization algorithms Many important algorithms for solving non-linear programs have been developed and refined in the last 10 years (Lasdon,1981). The more successful are : 1) Penalty function methods. The essential idea i s to transform the general non-linear problem into a sequence of unconstrained problems. The more robust code i s : SUMT ( C a r r o l l , 1959, Fiacco and Mc Cormicke, 1966). The objective function and inequality constraints can be non-linear functions but i f there are equality constraints they have to be li n e a r . Hence the code may be used to solve either the primal or the dual GP (Himmelblau, 1972). Another code , CONMIN, transforms a non-linear program with inequality constraints into an unconstrained problem using a penalty function method, and then the unconstrained problem i s solved by Fletcher and Powell's algorithm (Haarhof and Buys, 1970). 41 2) Generalized Reduced Gradient methods(GRG). The GRG algorithm uses the linear equality constraints to accomplish the equivalent of al g e b r a i c a l l y eliminating an equal number of dependent variables from the problem. This reduces the problem to one with exclusively decision variables. This reduction a f f e c t s the evaluation of the gradient. The reduced gradient is the gradient of the objective function with respect to the independent variables, subject to moving the independent variables in such a way that the equality constraints are s a t i s f i e d (Westerberg, 1981). The reduced gradient i s used to determine the dir e c t i o n of search. Hence, when GRG solves the dual GP problem, i t r e a l l y works with the "reduced dual problem", although the user should present the dual GP in i t s standard form. When dealing with inequality constraints, GRG converts them into e q u a l i t i e s by introducing slack variables. If the constraints are non lin e a r , they are replaced by their second-order Taylor series approximates expanded at the point of inter e s t . So GRG may also be used to solve the primal problem. UBC has a 1975 version of the code available. 3) Succesive Linear Programming (SLP). The Successive Linear Programming codes l i n e a r i z e any non-linear objective or constraint functions around a point, and then use the resulting linear program using e f f i c i e n t LP codes as subroutines. Some codes have been proposed by G r i f f i t h and Stewart (1961) and by Busby (1974) 42 Computational experience on GP codes and general non-linear codes for GP To decide which is the best algorithm to solve geometric programs would involve trying a l l of them on a wide variety of examples; the algorithms should be c o d i f i e d by the same programmer, and they should be run with the same compiler on the same computer. The evaluation can then be made on the basis of time. This approach, however , seems quite impractical. The number of proposed algorithms is quite large; so i s the range of problems to be tested. To be able to compare codes in d i f f e r e n t computers and with d i f f e r e n t compilers, standard times were defined. They refer to the execution time of the problem divided by the time required to execute C o l v i l l e ' s timing program (Himmelblau, 1972). They s t i l l vary with the programmer, and sometimes with the compiler, and they are not always used in the l i t e r a t u r e . We w i l l now show the code comparisons done by a series of authors on d i f f e r e n t examples of geometric programs. We w i l l focus on their conclusions when solving geometric programs that are also chemical equilibrium problems. Sarma, Martens, R e k l a i t i s and Rijckaert (1978) tested 5 algorithms to solve 16 GP problems . The codes were: SUMT on primal, GGP on primal, CSGP on dual, MCS on dual,and DAP on transformed primal. The basis for each code was given in the previous section in this chapter. GGP got the best results as for CPU times. SUMT had the largest times. They concluded that there was no d e f i n i t e evidence that solving the dual GP was 43 computationally more a t t r a c t i v e than solving the primal as has been stated by Duff in and Petersen (1967). Rijckaert and Martens (1978) also performed computational comparisons on GP algorithms. They t r i e d 16 GP codes ( a l l the ones we have mentioned before and three more) and one general purpose non-linear program, CONMIN , on 24 problems . Their problem 4 i s an scaled version of a primal chemical equilibrium as stated in Dembo (l976)--see prob.4, appendix B. On this p a r t i c u l a r problem the results of CPU times of the best algorithms are reported in table 111 — 1. Table 111-1. Best CPU times (in sec.) for the algorithms compared by Rijckaert and Martens (1978) on the chemical equilibrium problem. Problem GPKTC GGP NRF CONMIN (primal) (primal) (red. dual) (primal) Problem 4 1.0 5.82 7.73 12.73 For the general case, in t h i s comparison, GGP was the best code. But GPKTC worked better for the chemical equilibrium problem. Dembo(l978) compared six GP codes and five general non-linear optimization codes on 8 problems. Problems 1A and 1B of his series are, respectively, the unsealed and scaled versions of the primal chemical equilibrium problem (prob. 4, appendix 44 B),already compared by Rijckaert and Martens. The results for the scaled problem are shown in table I I I - 2 . Table I I I - 2 . Best standard times (in sec.) for the algorithms compared by Dembo (1978) on the chemical equilibrium problem. Problem GPKTC GGP GEOGRAD GRG (primal) (primal) (dual) (primal) Problem 1B .0554 .271 .0565 .0560 The time relationship between GPKTC and GGP i s similar, but not exactly the same as in the previous comparison. There i s not much difference between GPKTC , the dual-based GP code GEOGRAD , and the general- purpose non-linear code GRG . Only GPKTC performed in the badly scaled problem. Gochet,Loute and Solow (1974) t r i e d one GP algorithm for solving the primal (FP), one GP algorithm for solving the dual (CSGP modified) and GRG on the primal GP on 16 problems. Problems 3 and 8 of their set are examples of the chemical equilibrium problem. Problem 3 i s our problem 1 in Appendix B. It has 3 elements, 10 species and one phase; problem 8 has 6 elements, 16 species and 3 phases. The results are shown in table 111-3 . 45 Table III-3 . Best CPU times (in sec.) for the algorithms compared by Gochet, Loute and Solow (1974) on two chemical equilibrium problems. Problem CSGP FP GRG (dual) (primal) (primal) 3 1.68 .98 2.90 8 4.55 4.36 3.34 On the smaller problem, the primal based GP code was faster. GRG performed better on the bigger problem. Another conclusion of their work was that the application of GRG to the primal t-form of geometric programs was more e f f i c i e n t than using the In-transformed variables. This conclusion was v a l i d when the number of primal variables was less than 5. Ratner, Lasdon and Jain (1978) compared the perfomance of GRG on the Dembo set of problems (1976) with that of Dembo's code. They also compared GRG on the 24 problems given by Rijckaert and Martens with the best times given in their paper. Standard times were used as a comparison. They reproduced the results on Dembo's paper already mentioned. They also pointed out the s e n s i t i v i t y of GRG to some tolerances that affect i t s perfomance. They are the tolerances for the objective function (stopping c r i t e r i a ) and for binding constraints, and the use of quadratic or tangent approximations of the i n i t i a l values of the basis. F i n a l l y , Eckes,Gochet and Smeers (1978) discussed the d i f f i c u l t i e s of solving the dual GP with GRG. They modified GRG 46 using Beck and Ecker modification for the concave simplex in CSGP to account for the n o n - d i f f e r e n t i a b i l i t y of the objective function at a point where some of the variables are zero, and to allow for some variables to become zero at optimality. The modified GRG was compared with the CSGP code. Half the number of iter a t i o n s were needed with GRG than with the CSGP in most of the cases. From these results, we may conclude that: a) It i s not always clear when the dual or the primal GP should be used. b) For general cases, the best GP codes seem to be GGP (primal) and CSGP (dual). However, for the chemical equilibrium problem, GPKTC (primal) and GEOGRAD (dual) give better results in terms of time. c) The general purpose non linear code GRG compares well with the best GP codes when solving geometric programming problems, s p e c i f i c a l l y chemical equilibrium's scaled primal. The code can be used to solve either the primal or the dual GP . d) In a general case, GRG seems to work better when solving primal geometric programs on the t-variables. However, the GP primal codes that performed better on the chemical equilibrium problems worked with transformed primal variables. 47 Geometric programmminq and the chemical equilibrium problem As we pointed out in the previous chapter, the chemical equilibrium problem reduces to a special case of geometric programming (GP). If we translate into GP nomenclature the t r a d i t i o n a l ways of formulating the chemical equilibrium problem, we can see that they a l l attempt to solve a dual GP, either as a problem with a l l i t s variables e x p l i c i t or as a reduced problem . Passy and Wilde. (1968) devised a primal algorithm to solve the chemical equilibrium problem with only one ideal phase. It involved the formation of an unconstrained problem using Lagrange's method. They t r i e d i t on a hydrazine combustion problem with one ideal phase, 10 species and 6 elements (see problem 1 , appendix C). The algorithm worked with the primal variables t, which were scaled between 0 and 1. The objective function was elevated to the one tenth. The RAND corporation model of the chemistry of the human respiratory system was used as a test problem by Dembo in a set of problems intended to compare GP codes (Dembo, 1976). We saw the results in the previous section. B l e i g h t l e r and P h i l l i p s give.an extense discussion of t h i s p a r t i c u l a r problem solved by GGP (1976). Dinkel and Lakshmanan (1975,1977) were interested on s e n s i t i v i t y analysis applied to the chemical equilibrium problem. They solved the dual problem using Dinkel and Kochenberger's GP algorithm (GEOGRAD code). F i n a l l y Lidor (1975) devised a modification of the GGP 48 algorithm to solve the primal chemical equilibrium problem. He worked with the transformed primal variables z= In t, and used a Zangwill's cutting plane algorithm. His code includes the generation of a f i r s t starting point. He t r i e d the algorithm on seven chemical equilibrium test problems, but his CPU timings compared disfavourably to the RAND code. The RAND method, as we pointed out e a r l i e r in thi s chapter, i s b a s i c a l l y a dual-based algorithm, and the RAND code has been perfected over many years. The question remains as i f Lidor's code was slow because of the algorithm i t s e l f , or because of i t s implementation. 49 CHAPTER IV: WRITING THE PROGRAM—PRELIMINARIES In order to write a general computer program to solve the chemical equilibrium problem we had to decide which mathematical formulation and which numerical technique of solution should be used. As for a mathematical formulation, we decided to use the geometric programming approach, because i t provided a new form of presenting the problem (the primal formulation), and i t gave the basis for a s e n s i t i v i t y analysis. We wanted to compare the primal and dual formulations of the geometric program, since no d e f i n i t e conclusions about the superiority of either one could be drawn from the l i t e r a t u r e . The dual GP is equivalent to the t r a d i t i o n a l free energy minimization approach to the chemical equilibrium problem. The dual variables have a straightforward physical meaning: they are the number of moles of the species present at equilibrium. The disadvantages are the bigger dimensionality of the dual, the n o n - d i f f e r e n t i a b i l i t y of the objective function at zero ,and the numerical problems that result when one variable tends to zero (Gochet et a l . , 1 9 7 4 ) . The primal problem has less variables and few constraints. The constraints are highly non-linear, but they are inequality ones, and they are active at the optimum. Two related problems remain: the scaling of the problem and the finding of a f i r s t primal s t a r t i n g point. We decided to try one computer code that could be used for both the primal and the dual formulations, in order to compare 50 the two cases. GRG had^ been proven in the l i t e r a t u r e as comparable, i f not better than most general-purpose GP codes. It could also be applied on the dual and the primal problems. A 1975 version of the code was available at UBC (Wales, 1977) and could be used as a self-contained program and as a subroutine. We decided to use GRG as our comparison code. In t h i s preliminary part of the work we were concerned with the accuracy of the r e s u l t s . We had to determine the best conditions for GRG to get results comparable to the l i t e r a t u r e . The general scaling of the problems was a d i f f i c u l t task, and three d i f f e r e n t versions of solving i t were t r i e d . A l l this was performed with GRG as a self-contained program. The rest of the chapter i s organized in the following way. In the f i r s t two sections we respectively repeat the chemical equilibrium formulations that we used, and give a detailed description of GRG. The t h i r d section deals with the scaling problem. The forth section states the working parameters for GRG. In the next section we present the method used to generate a s t a r t i n g point. The steps taken for performing s e n s i t i v i t y analysis are shown in the l a s t section. 51 Mathematical formulations For convenience we s h a l l repeat here the mathematical formulations of the primal and the dual chemical equilibrium problems as we used them. They correspond to equations 2.40, 2.41, 2.42 and 2.43 of chapter II. For programming reasons, we incorporated the normality condition to the exponent matrix and to the B vector, and wrote the B vector e x p l i c i t l y . A l l e subscripts are shifted one unit; A11=B1=1, Ai1=0 for a l l i , A1j=0 for a l l j>1 Dual problem Min (-v) N+1 exp(G/RT) = n ( C J / 6j ) j = 1 6j K ] n X k k=1 s. to N+1 E A i j 6j = Bi j = 1 i = 1,M+1 6j £ 0 j = 1,N+1 4.1 Transformed dual N+1 Min (-In v) = G/RT = E 6; (In 6 - + Cj) -j = 1 K I X. Kln X K k=1 s. to N+1 E A i j 6J = = Bi i = 1,M+1 j=1 - 6j£ 0 j • = 1,N+1 Primal problem M Min g 0 ( t ) = n t i . i = 1 -Bi+1 s. to g K ( t ) = E c j J €k M Ai-M^j n t; . < 1 k = i = 1 = 1 ,K tj >> o i = = 1 ,M 53 Transformed primal M Min In g c ( t ) = h (Z) = = E -Bi+1 Zi i=1 s. to M E exp(Cj + E Ai + 1 , j Zi) <? 1 k = = 1 ,K jek i = 1 GRG. Description of the code The program for the GRG code available at UBC was written at Cleveland State University. The.manual UBC GRG, written by K. Wales (1977) is largely taken from Cleveland State University technical memorandum C1S-75-02 (1975). We s h a l l repeat here some of the more important features of this code. Stating the problem with GRG Warning: the nomenclature used by GRG clashes with ours. Later in t h i s chapter we include a table to compare them (Table IV-3) GRG solves constrained optimization problems stated as follows: 54 Min g M + 1 ( x ) s. to 9f = 0 i = = 1,NEQ o -< g UB w + r ) i = = NEQ+1, M LBi X UBi i = = 1,N Where: x = (x.,, ... x N) vector of N real variables. G = (g 1 , ... gM^ .^) vector of real continuous functions of x , linear or non-linear. 9-1 • • • 9uE<a equality constraints. g N e „ + 1 ••• 9M inequality constraints. g M 4^ objective function. LBi lower bound on x; i = 1,.. N UBi upper bound on x i = 1,.. N UB upper bound on the inequality constraint q ; i = NEQ + 1, M Description of the algorithm We w i l l now give a brief description of the algorithm. M slack variables are added to the constraints of problem ( 4 . 5 ) . The previous N variables are c a l l e d natural variables. Assume x i s a feasible point, and NB of the constraints are binding at x . Two sets of variables are distinguished in the GRG algorithm, provided there i s no degeneracy : the NB basic variables (dependent) and the non-basic variables which are the N-NB remaining natural variables and the M slack variables 55 associated with the binding constraints. The binding constraints can be written as G (Y,Z)= 0 4.6 where Y i s the vector of the NB basic variables and Z is the vector, of non-basic variables. Then the binding constraints may be solved for Y in terms of Z, y i e l d i n g a function Y(Z), v a l i d for a l l Y,Z close to the f i r s t feasible point. This reduces the objective function to a function of only the Z variables, F(Z) which i s c a l l e d the reduced objective. The gradient of F(Z) is c a l l e d the reduced gradient. The o r i g i n a l problem is now a reduced problem. Min F (Z) s. to LB <^ Z ^ UB 4.7 Where LB and UB are respectively the lower and upper bounds for Z, and the vanishing of the reduced gradient is sought. GRG solves (4.5) (with the slack variables) by solving a sequence of (4.7) reduced problems. Each problem i s solved by a gradient type unconstrained non-linear optimization method. The reduced gradient gives a search d i r e c t i o n d for a one dimensional linear subsearch as to minimize F(Z+ ad) with respect to a . A set of f i n a l equations is solved using Newton's method. 56 If Newton does not converge, GRG reduces a and t r i e s again. Otherwise the search i s finished. If Newton converges, but a basic variable v i o l a t e s the bounds, a new set of basic variables is determined and the solution of a new problem s t a r t s . The search may continue u n t i l the objective i s found larger than the one in the previous i t e r a t i o n . Then a quadratic interpolation is done to the three o values bracketting the minimum, and the objective i s evaluated at t h i s point. If the i n i t i a l x does not s a t i s f y the constraints, GRG optimizes a "Phase 1" objective which i s the sum of the constraints v i o l a t i o n s , in order to obtain a f i r s t feasible point. Use of GRG GRG can be used as a self-contained program or as a subroutine. The self-contained program was used in the preliminary stages of this work. Details on the use of GRG can be taken from the UBC-GRG writeup (Wales, 1977). Some parameters have to be s p e c i f i e d for GRG to work. The values of these parameters depend on the mathematical formulation of the problem, on the scaling of the problem, and on the accuracy required for the solution, and w i l l be discussed l a t e r . Running GRG as a self-contained program , these parameters are prompted by the user in a terminal on conversational mode. The objective function and the constraints are calculated by a subroutine named GCOMP, provided by the user 57 in a f i l e attached to input unit 5. The scaling problem The fact that the chemical equilibrium i s an i l l - c o n d i t i o n e d problem was pointed out by many authors (Beightler and P h i l l i p s , 1976; Dembo, 1976; Lidor, 1975). One d i f f i c u l t y i s the huge value of the GP objective function; the range of the primal variables at the optimum i s quite wide. The free energy c o e f f i c i e n t s also vary a l o t between species. Table IV-1 exemplifies the need for a scaling with some problems from the l i t e r a t u r e . A l l the figures in table IV-1 are from the l i t e r a t u r e or they were calculated from l i t e r a t u r e data. Problem 1 was solved in the l i t e r a t u r e (Passy and Wilde, 1968) as a primal problem. The t variables were scaled so that their values at the optimum had a range of 100 between the lower and the higher values. No comments were made in the l i t e r a t u r e on the logic of the scaling. Problem 4 is one of Dembo's test problems for evaluating geometric programming codes. GRG could not solve the unsealed version when using the t variables. The problem was scaled by Dembo, and i t s results were reproduced by us, star t i n g from a point :closer to the optimum than t h e i r s . The transformed primal problem could be solved unsealed. The sources of problems 2,6 and 7 did not specify the value of the objective function, which we calculated from their composition and free energy data. Problem 5 was solved in the 58 l i t e r a t u r e as a dual GP, using the logarithm of the objective funct ion. As a f i r s t approach, we t r i e d to reproduce the problems from table IV-1 using GRG. We succeeded in the scaled and transformed problems, but not on the unsealed versions of the primal or dual geometric programs. See appendix B for our re s u l t s . Attempts to solve the scaling problem From what we could see in the l i t e r a t u r e , the attempts done to scale the primal problem were: a) when working with the primal variables t, the objective function was elevated to an exponent so as to diminish i t s absolute value, and the primal variables were scaled in obscure ways so that they would f a l l within a close range from their optimum value. b) when working with the transformed primal variables Z, there was no scaling problem. The objective function corresponded to the standard Gibbs free energy of the system and the transformed primal problem was a convex program. As can be observed in table IV-1, the range of variation of the transformed primal variables i s not too wide. We t r i e d both schemes a) and b) on problems 5 and 6, running T a b l e IV-1 . C h a r a c t e r i s t i c s of the Chemical E q u i 1 i b r i u m Problem' PROBLEM NUMBER SPECIES, ELEMENTS, PHASES. OBJECTIVE FUNCTION @ OF TIMUM RANGE OF VARIABLES <a OPTIMUM V * V • 0 1 1 n v t * * Zt 6* 1 10. 3, 1 6.2904 E20 1 .61431 47.8907 5.61 E-5 ?2.47 E-7 -9.7882 -15.214 1.48 E-1 6.93 E-4 2 4 . 3 , 1 1.9699 E39 2.47141 90.4788 4.93 E-2 2.86 E-14 -5.3120 -33.488 2.48 E-1 2.52 E-1 3 5 , 3 . 1 2.9231 E34 2.21136 79.3605 4.25 EO ' 1.49 E-10 -.85543 -24.927 1.72 E-1 5.79 EO 4 30,12, 3 7.8257 E796 9.30969 E7 1834.91 6.51 E-1 1.19 E-9 -.42960 -20.552 6.60 E-21 28.9 E'1„ 5 10, 4, 1 5.2990 E377 2.38365 E3 777.638 1.10 E-5 1.80 E-9 -11.415 -20.127 4.40 E-4 19.9 E1 6 8 , 4 , 1 1.4576 E52 3.3281 120.111 5.84 E-3 3.32 E-20 -5.1423 -44.851 1.4 1 E-4 1.88 EO 7 24, 4, 1 1.4448 E52 3.32352 120.102 6.07 E-3 3.32 E-20 -5.1037 -44.852 5.77 E-22 1.88 EO 1 The sou r c e d a t a f o r t h i s t a b l e a r e i n Appendix B * v = exp (G/RT) ** t = p r i m a l v a r i a b l e s t Z = Int * 6 = dual v a r i a b l e s (number of moles) 60 GRG as a self-contained program. When using the primal variables t , the objective function was elevated to 0.01 and the primal variables t were scaled based on a primal s t a r t i n g point close to the optimum. After the scaling, the variables were a r b i t r a r i l y bounded between 0 and 500. When using the transformed primal variables Z, we were faced with the problem of f i x i n g the boundaries for these variables. We put zero as upper bound for a l l the variables. The tranformed primal objective function is linear on the Z variables; each Z can be regarded as the contribution of the associated element to the t o t a l free energy of the system (Zeleznik and Gordon, 1968). Therefore, i f Zi i s equal to zero that means that the element i is not contributing to the t o t a l free energy of the system, regardless of the amount of element i . This i s quite an u n r e a l i s t i c situation i f the model i s well posed, and so zero looks l i k e an appropiate upper bound for the Z variables. The lower bound i s not so ea s i l y determined. A fixed large negative number l i k e -100 i s one p o s s i b i l i t y , but this approach has some problems: i t does not account for cases when Z is below that boundary and the range may be too wide for some variables. Also, in the evaluation of the terms of the primal constraints we can have exponentials of too large negative numbers, which are undefined. A second p o s s i b i l i t y to determine general lower bounds for the transformed primal variables is to calculate a fixed percentage of a f i r s t good s t a r t i n g point and use i t as a 61 boundary. But then there i s the question of how big t h i s percentage should be; we t r i e d a few numbers, and found variations from problem to problem. The t h i r d p o s s i b i l i t y i s to add a fixed negative number to a good sta r t i n g point. If we choose -30, that i s equivalent of a range of 10 1 3 in the untransformed primal variables t. If we choose -20, the range i s 10 8. We chose -30 a r b i t r a r i l y and i t worked well in a l l our examples . This selection of a boundary is then strongly dependent on the f i r s t s tarting point. Using an aproximation of the problem to a linear program allowed us to obtain dual starting points that are quite close to the optimal values. The corresponding primal points can be calculated from the dual. A description of the procedure is done later in t h i s chapter. See Chapter VI for the results of the closeness of these points to the optimum values. Comparison between scaling the primal and using the transformed primal Table IV-2 shows the results of a comparison between the two procedures explained in the previous section. The problems No. 5 and 6 of table IV-1 were scaled, and also posed as transformed primal problems. Both problems have 4 primal variables, and according to the l i t e r a t u r e ( Gochet et a_l. , 1974) for general GP the problem posed as t variables should be solved more e f f i c i e n t l y by GRG than the transformed one. The number of i t e r a t i o n s and of function evaluations needed to go from the same star t i n g point to the .same optimum were computed with both 62 Table IV-2. Performances of the scaled primal problem and the ln-transformed primal with GRG (*) Problem Scaled Ln-transformed Iterations F. evaluations Iterations F. evaluations 5 (**) 7 51 5 46 g (*** ) 20 240 13 183 Notes: (*) GRG parameters: EPSTOP=EPNEWT=EPSBOUND=ESPIV= 10"6 (**) Scaling: objective function elevated to .01 ; y1= 11.105 y2 = t2. 106 , y3 = t3 . 106 , y4 = t4 . 105 Starting point : Z1=-11 . 42, Z2 = -20.13, Z3 = -15.86, Z4 = -11.63 Optimum : Z1=- 1 1 .4315, Z2 = -20.13541, Z3 = -15.8013, Z4=-11.6964 (***) Scaling: objective function elevated to .01 ; y1=t1.l0 y2=t2. 10 2 3 , y3=t3. 107 , y4=t4 . 105 Starting point: Z1=-4.35, Z2 — 51.6, Z3 = -16.3, Z4=-11.9 Optimum: Z1=-4.3370, Z2=-52.2190, Z3=-15.9968, Z4 = - 1 1 .9352 methods. The same tolerances were used in a l l cases. It took GRG less i t e r a t i o n s and function evaluations to solve the transformed primal, the opposite of what the l i t e r a t u r e said for general GP . We believe t h i s is due to the p a r t i c u l a r i l l -conditioning of the chemical equilibrium problem, and therefore should need a special scaling for each example. From here on, we w i l l use the words "primal" and "transformed primal" i n d i s t i n c t l y , and we w i l l always refer to the ln-transformed problem.. For the dual problem, we found no special computing advantages for either using the objective function as in a 63 standard geometric program (but elevated to a certain exponent), or minimizing the negative logarithm of that function. The second approach was a t t r a c t i v e because then the objective function represented the Gibbs free energy of the system divided by R and by T, same as in the transformed primal, and as in the t r a d i t i o n a l "free energy minimization method". For the rest of the thesis , whenever "dual problem" is mentioned, we w i l l be ref e r r i n g to equations 4.2. Determination of working parameters for GRG GRG needs as an input a series of control cards that have to be provided by the user. We determined some of them from l i t e r a t u r e values; some others through numerical experience. The res u l t i n g values may not be the best for a s p e c i f i c problem; however, they allowed us to obtain results similar to the l i t e r a t u r e in a l l cases. Table IV-3 explains the meaning of the key words and exemplifies the values used in our programs. A dual sta r t i n g point A sta r t i n g point for the dual problem is quite straightforward, since the dual variables are the number of moles of the chemical species. When using GRG as a s e l f -contained program, we usually started with a point close to the l i t e r a t u r e optimum, and then we changed i t at random. For the user's convenience, the f i n a l code should include the generation of a st a r t i n g point. The NASA code uses an equal TABLE IV-3. Parameters for GRG GRG Transformed DUAL PARAMETERS PRIMAL N < 100 (Number of var iables) M + 1 (Number of elements + 1) N + 1 (Number of spec ies + 1) M < 100 (number of constraints) K . (number of phases) M + 1 (number of elements + 1) NEQ (number of equality constraints) o N + 1 (number of elements + 1) LBV (lower bound on variables Starting point +(-30) E-2 4 (E-10 - E-30) - t r i e d -UBV (upper bound on variables) 0 40 UBC (upper bound on inequality constraint) 1 0 X ( i n i t i a l vector of variables) Generated with LIPSU2 and SINGV Generated with LIPSU2 EPNEWT (tolerance for equality and binding constraints) E-6 (very important for accuracy) E-6 (very important for accuracy) EPSBOUND (tolerance for var iables @ bounds) E-6 E-6 //cont inues / / / T a b l e . I V - 3 (continued) GRG Transformed DUAL PARAMETERS PRIMAL EPSTOP E-6 E-6 (tolerance - for objective function stopping c r i t e r i a ) EPSPIV E-6 E-6 (tolerance for pivot element in the basis) QUAD used used (quadrat ic • -important- -important-extrapolation for estimating i n i t i a l basic • variables) PRINTCTL. 2 2 • (control for s e l f - for s e l f -amount of contained contained output) program;1 for program;1 for subrout ine 'subroutine 66 d i s t r i b u t i o n of species for that purpose. The RAND code approximates the dual problem to a linear program to obtain a f i r s t point. Lidor (1975), did a di f f e r e n t approximation than the NASA , and obtained another linear program for his starting point routine. Lidor compared his results to those of the RAND code, and got points closer to the optimum. We decided to implement Lidor's approach. His formulation of the linear program (LP) i s as follows: LP approximation of the dual N N Min f = E Cj yj + € E Cj j = 1 j = 1 s. to N E A i j yj + £ = Bi i = 1 ,M j = 1 £ > « Yj > 0 j = 1,N Then: <5j = yj + t 4.9 To solve the linear program, we used subroutine LIPSUB from UBC (Patterson, 1979) on the program LIPSU2 . The following steps were followed: 67 1) Read C, A and B. 2 Set the tableau in the manner spec i f i e d by the UBC routine LIPSUB . 3) C a l l LIPSUB . It solves the linear program using a primal-dual algorithm. 4) The results (the y variables and z ) are used to calculate the dual st a r t i n g point with equation 4.9. Calculation of primal variables from the dual ones The primal starting point i s more d i f f i c u l t to v i s u a l i z e and is more c r i t i c a l than the dual. We already mentioned that we are defining the lower bounds of the primal variables as a certain function of their s t a r t i n g value. Since we already have a program to generate a dual starting point, i t sounds reasonable to transform this point to primal variables. For a geometric program at the optimum , the primal and dual variables are related as follows: M E A i j Zi = Cj + In (6jA k) i = 1 where There are M variables and N equations. In reacting chemical equilibrium problems, N>M and we have an overdetermined system of linear equations. The e q u a l i t i e s hold only at the optimum; j = 1,N 4.10 k = 1 ,K 68 this i s a nuisance i f we are interested in the starting point which is not normally the optimum. A linear least squares approach should be the best solution. We t r i e d subroutine DBEST from UBC with very good r e s u l t s . It decomposes the transpose of our A matrix (the exponent matrix) following a Gram-Schmidt orthogonalization. The problem was that the UBC subroutine was limited to cases with less than 30 species, so we t r i e d another method. Subroutine DSLVD from UBC (UBC MATRIX, 1980) uses a singular value decomposition to solve overdetermined systems of linear equations. The transpose of the matrix A i s decomposed as follows: A' = U I V' 4.11 where U, V = orthogonal matrices E =diagonal matrix of the singular values of A' The solution of the system of equations A'Z = b is the vector Z Z = V E + U' b 4.12 Where E + is the pseudoinverse of E. We implemented the program SINGV as follows: 1) The exponent matrix A, the free energy c o e f f i c i e n t s matrix C, the dual sta r t i n g point and the values of M and N are read. 69 o 2) C a l l subroutine DLSVD from UBC. It calculates the matrices V, E and U'b. 3) Calculate E + and V E +. 3) subroutine DGMULT from UBC i s c a l l e d to multiply V E + and U'b and obtain the vector Z, which i s actually the vector of the transformed primal variables z. S e n s i t i v i t y analysis As we already discussed elsewhere (chapter I I ) , there are two approaches in the l i t e r a t u r e to s e n s i t i v i t y analysis in geometric programming. Both approaches are mathematically equivalent but the f i r s t involves the inversion of a matrix and the second solves a system of linear equations. We chose the combustion of propane (problem 5) as a means of comparison between the two approaches, because : a) It was used by Dinkel and Lakshmanan in their paper on s e n s i t i v i t y analysis for the chemical equilibrium problem. They used the method of the inversion of the Jacobian (Dinkel and Lakshmanan, 1977). b) i t i s a "middle s i z e " problem within the l i t e r a t u r e . It involves six simultaneous reactions ; the dimension of the matrix to be inverted in the Jacobian approach is DxD, where D is the number of independent reactions . When D= 1 or 2, there is no doubt that t h i s method w i l l be faster than solving an equivalent system of N+K+1 by N+K+1 linear equations, where N= number of species, K is the number of phases. We showed before (Chapter II) that D=N-M; i f we assume N=3, M=2, K=1, we have a 5 70 x 5 system of equations to be solved vs. elevating a number to -1. In the combustion of propane example, D=6, N=10, M=4, K=1 and then the comparison is between a 6 x 6 matrix to be inverted, or a 12 x 12 system of l i n e a r equations to be solved. For the purpose of the comparison the stoichiometric c o e f f i c i e n t s used were the ones calculated by Dinkel and Lakshmanan (1977). We w i l l now explain a systematic procedure to obtain these c o e f f i c i e n t s from the exponents matrix AA. Calculation of stoichiometric c o e f f i c i e n t s A systematic method to construct a maximal set of l i n e a r l y independent chemical reactions from the exponent's matrix AA i s as follows: ( C a v a l l o t t i |et a l , 1979) a) Given the matrix AA, factorize i t into two matrices E and F, where E i s an MxM matrix of rank M, and F has dimensions MxD, so that, in m a t r i t i a l notation: AA = | E | F | In chemical language, that is equivalent to reacc'ommodating the matrix AA, so that the f i r s t M columns correspond to the key components; the matrix F i s the matrix of the constituents. b) Reduce the factorized matrix by a Gauss Jordan procedure u n t i l i t becomes AA = | I | H | where I i s the MxM identity matrix. c) Form the array : 71 U = -H I Where I i s now the DxD identity matrix, and U stands for the NxD matrix of stoichiometric c o e f f i c i e n t s for the D reactions. With thi s construction, the product of matrices AA.U = 0 , which was the condition . required for the stoichiometric matrix (see chapter I I ) . Comparison between the two methods of s e n s i t i v i t y analysis We implemented Dinkel and Lakshmanan's method in the program JOTA (appendix B). The method based on the solution of the linear system of equations i s contained in program NPLUSK (appendix B). We w i l l now describe both algorithms. a) Algorithm of program JOTA 1) read optimal composition at optimum temperature and pressure (data from the l i t e r a t u r e ) , stoichiometric c o e f f i c i e n t s , t o t a l number of moles in the phase, free energy c o e f f i c i e n t s at the optimal and at the perturbed T,P. 2) form the J matrix according to equations 2.44, 2.45, 2.46, 2.47. 3) the J matrix i s inverted by the UBC subroutine INV (UBC Matrix, 1979). The execution time to invert the matrix i s printed. 4) for the new set of free energy c o e f f i c i e n t s , the new composition and free energy are computed using equations 2.49 72 and 2.48, and .these values are printed. B) Algorithm of program NPLUSK 1) read same data as in JOTA . 2) the matrix to be solved i s posed at the optimal conditions of P, T, and composition. The equations 2.50 to 2.53 are used for this purpose. 3) the right-hand side vector i s also calculated for the perturbed condition. 4) the system of equations i s solved by UBC subroutine DSLIMP (see UBC MATRIX, 1979). The execution time for the performance of DSLIMP i s printed. 5) the new composition values are printed. In table IV-4, we can see that both programs y i e l d figures for composition comparable to those obtained in the l i t e r a t u r e for s e n s i t i v i t y analysis, but the execution time for the JOTA program was a thousand times longer. As a consequence of t h i s , we decided to incorporate the NPLUSK program as a subroutine of our main program in order to perform s e n s i t i v i t y analysis (see Chapter 5). When the free energy c o e f f i c i e n t s of the species change, only the right hand side of the system of equations needs to be recalculated. Since the matrix of the equations i s solved by DSLIMP, i t is very easy to make a new perturbation and resolve the system, for the matrix i s stored . It i s not the same when the amounts of elements are changed. In thi s case, the f i r s t column of the matrix has to be recalculated. No attempts were 73 Table IV-4. Two approaches to s e n s i t i v i t y analysis. Composition values for problem 5 * , at two thermodynamic conditions, using optimization. Composition at the perturbed condition, calculated with two methods of S.A. Execution times for the two methods. Spec ies Start T=2200 K P=40 at. Perturbat ion 2500 K , P = 50 at. Optimiz.* S.A. * JOTA ** t=2.24008 s NPLUSK *** t=0.00161 s H2 .02007 .05492 .04020 .04020 .04017 H .00065 .00433 .00189 .00189 .00189 OH .01500 .05579 .03534 .03534 .03532 H20 3.9719 3.9150 3.9412 3.9412 3.9412 CO .08160 .24835 .17378 .17379 .17267 C02 2.9184 2.7516 2.8262 2.8262 2.8263 N2 19.987 19.959 19.971 19.972 19.972 NO .02668 .08618 .05680 .05679 .05679 02 .03358 .09614 .06954 .06954 .06955 O .00044 .00356 .00137 .00139 .00138 T= temperature, K= Kelvin, P= pressure, at.= atmospheres, SA.= s e n s i t i v i t y analysis, Optimiz.= optimization, t= execution time, s= seconds. * Dinkel and Lakshmanan, 1977. ** JOTA : S.A. program based on inverting the Jacobian. *** NPLUSK : S.A. program that solves a system of linear equations. done on solving t h i s case. From the data in table IV-4 we can see that the re l a t i v e errors of the composition calculated by s e n s i t i v i t y analysis, are quite large. These errors are r e l a t i v e to the value of the composition as calculated by a new optimization. The errors also vary from species to species. Almost a l l papers on s e n s i t i v i t y analysis in GP (Rijckaert, 1974; Dinkel and Lakshmanan, 1975,1977) propose an incremental procedure to diminish these errors. An evaluation of the size of the increment • was done on the example No. 5, and is 7 4 e x p l a i n e d i n Chapter VI. S e n s i t i v i t y a n a l y s i s seems a f a s t e r way to solve the chemical e q u i l i b r i u m problem i n a region c l o s e to a known optimun, than r e o p t i m i z a t i o n i t s e l f . J u s t how much f a s t e r i t can be, with s i m i l a r accuracy, has not yet been s t a t e d . We t r i e d to compare the s e n s i t i v i t y a n a l y s i s procedure vs. r e o p t i m i z a t i o n , and the r e s u l t s and c o n c l u s i o n s are i n Chapter VI. Co n c l u s i o n s of t h i s chapter The c o n c l u s i o n s of these p r e l i m i n a r y s t u d i e s a r e : 1) GRG s o l v e s both p r i m a l and dual GP chemical e q u i l i b r i u m problems, provided there i s a s c a l i n g , or a l o g a r i t h m i c t r a n s f o r m a t i o n . I t i s more e f f e c t i v e i f ( i ) the p r i m a l i s posed as a problem of the transformed p r i m a l v a r i a b l e s Z= In t , ( i i ) the l o g a r i t h m of the dual o b j e c t i v e f u n c t i o n i s o p t i m i z e d . 2) The working parameters f o r GRG were s t a t e d . 3) A program to generate dual s t a r t i n g p o i n t s was w r i t t e n . 4 ) A program to c a l c u l a t e p r i m a l p o i n t s from the dual ones was implemented. 5) S e n s i t i v i t y a n a l y s i s f o r the v a r i a t i o n s of the f r e e energy c o e f f i c i e n t s on the dual v a r i a b l e s was performed f a s t e r by s o l v i n g a system of l i n e a r equations r a t h e r than by i n v e r t i n g the Hessian matrix. The problem t r e a t e d had a degree of d i f f i c u l t y ( number of independent r e a c t i o n s ) equal to s i x . 6) For s e n s i t i v i t y a n a l y s i s we need the s t o i c h i o m e t r i c c o e f f i c i e n t s f o r the r e a c t i o n s . A l i t e r a t u r e based systematic procedure to o b t a i n them from the exponents matrix was 75 e x p l a i n e d . 76 CHAPTER V: DESCRIPTION OF THE PROGRAMS Diagrams We w i l l present here the diagrams of the four computer programs that were w r i t t e n . The idea was to compare pr i m a l vs. dual perfomance, and s e n s i t i v i t y a n a l y s i s vs. r e -o p t i m i z a t i o n . The diagrams are q u i t e rough; a more d e t a i l e d e x p l a n a t i o n of each subroutine w i l l be given i n the f o l l o w i n g s e c t i o n s of t h i s chapter. Program COMP 1 so l v e s the p r i m a l transformed problem with GRG. Program COMP 2 so l v e s a l s o the p r i m a l , and performs s e n s i t i v i t y a n a l y s i s as w e l l . Program COMP 3 s o l v e s the dual and f i n a l l y program COMP 4 s o l v e s the dual and performs s e n s i t i v i t y a n a l y s i s . The programs are w r i t t e n i n FORTRAN IV . A l l subroutines are r e l a t e d through l a b e l l e d COMMON b l o c k s . Each subroutine c a l l s at l e a s t one subroutine from UBC (except f o r subroutine FREEN ). The programs were run on an AMDAHL 470 V/8 computer, using the Michigan Terminal System. They were compiled with IBM FORTRAN com p i l e r s G and H at l e v e l 21.8. Data were read from a f i l e a t t a c h e d to input u n i t 5. 77 PROGRAM COMP1 READ DATA SUBROUTINE FREEN C a l c u l a t e s f r e e energy c o e f f i c i e n t s at d i f f e r e n t T and P TIME = O.DO FIX T and P SUBROUTINE L I P S U 2 p r o v i d e s dual s t a r t i n g p o i n t SUBROUTINE SINGV c a l c u l a t e s p r i m a l v a r i a b l e s from dual GRG s o l v e s the o p t i m i z a t i o n problem (GCCOMP) CALCULATE dual v a r i a b l e s from the optimal p r i m a l WANT VAR .T AND/OR P? -VARY T,P. -GET new f r e e energy c o e f f , -USE optimum as s t . Point ' PRINT RESULTS/ PROGRAM COMP2 READ DATA SUBROUTINE FREEN C a l c u l a t e s f r e e energy c o e f f i c i e n t s at d i f f e r e n t T AND P TIME = O.DO 1 FIX T AND P WANT VARY T AND/OR P? SUBROUTI Provides dual WE LIPSU2 s t a r t i n g point SUBROUTINE SINGV C a l c u l a t e s p r i m a l v a r i a b l e s from dual ones GRG Solves the o p t i m i z a t i o n problem CALCULATE du from the op a l v a r i a b l e s t i m a l p r i m a l YES -VARY T and P. -GET new fr e e energy c o e f f . SUBROUTI] Performs Sensi WE NPLUSK t i v i y A n a l y s i s 'PRINT RESULTS/ 79 PROGRAM COMP3 READ DAT SUBROUTINE FREEN C a l c u l a t e s f r e e energy c o e f f i c i e n t s at d i f f e r e n t T AND P TIME = 0.D0 FIX T AND P SUBROUTIb Provides dual s IE LIPSU2 s t a r t i n g p o i n t GF Solves the c Problem ( IG >ptimization GCCOMP) PRINT 'RESULTS/ WANT VARY ^ P AND/OR T? T YES -VARY P,T. -GET new f r e e energy c o e f f -USE Optimum as s t . Point NO TIME; STOP 80 PROGRAM COMP4 READ DATA/ SUBROUTINE FREEN C a l c u l a t e s f r e e energy c o e f f i c i e n t s at d i f f e r e n t T AND P TIME = O.DO FIX T AND P SUBROUTIls p r o v i d e s dual £ IE LIPSU2 s t a r t i n g p o i n t GF Solves t opt i m i z a t ic IG .he dual >n problem YES -VARY T, P. -GET new fr e e energy c o e f f SUBROUTINE NPLUSK Performs s e n s i t i v i t y a n a l y s i s NO TIME/ STOP PRINT RESULTS 81 Input needed. Examples For any of the programs mentioned above, the user should p r o v i d e : a) c o n t r o l c a r d s : N+1= number of spe c i e s present at e q u i l i b r i u m +1 M+1= number of elements +1 K= number of phases NN= number of times P and/or T are v a r i e d . NE= 1 i f the f r e e energy c o e f f i c i e n t s are a v a i l a b l e as such; i f they have to be c a l c u l a t e d from the Gordon and Mc. Bride c o e f f i c i e n t s , use any other i n t e g e r . NF=1 i f the dual s t a r t i n g p o i n t r o u t i n e i s to be used. Otherwise, use another i n t e g e r . b) f r e e energy data: f r e e energy c o e f f i c i e n t s matrix C or temperature c o e f f i c i e n t s matrix S f o r the Gordon and Mc. Bride polynomials, NASA SP-3001. c) exponent matrix A. An example of i t s c o n s t r u c t i o n w i l l f o l l o w . d) amount of elements :vector B. e) s t o i c h i o m e t r i c c o e f f i c i e n t s f o r the r e a c t i o n s i f C0MP2 or C0MP4 are to be used (matrix U) 82 e) dual s t a r t i n g p o i n t i f a v a i l a b l e . We s h a l l now g i v e two examples of the i n p u t . Example 1• CH4-Water gas r e a c t i o n P r e l i m i n a r y data: -Reaction of 2 moles of CH4 and 3 moles of water -Reach e q u i l i b r i u m at T= 1000K, P= 1atm. -Species present at e q u i l i b r i u m = CO, C02, H20, H2,CH4 - I d e a l i t y assumed -One phase (gas) -Free energy c o e f f i c i e n t s a v a i l a b l e from the l i t e r a t u r e f o r one set of T,P c o n d i t i o n s . -Do not have a s t a r t i n g p o i n t . Input a) C o n t r o l c a r d s : N+1 = 6 ( 5 s p e c i e s +1 ) M+1 = 4 (3 elements +1) K = 1 (1 phase) NN = 1 (one thermodynamic s t a t e ) NE = 1 (the f r e e energy c o e f f i c i e n t s are data) NF = 1 (no s t a r t i n g p o i n t i s provided) b) Free energy c o e f f i c i e n t s : C i s a 6x1 matrix, as f o l l o w s : 83 0.0 dummy -24.025 CO -47.413 C02 C = -23.067 H20 0.0 H2 2.0847 CH4 The f i r s t c o e f f i c i e n t corresponds to the n o r m a l i t y c o n d i t i o n . The r e s t of the c o e f f i c i e n t s are the Gibbs f r e e energy of formation from i t s elements at temperature = 1000K f o r each s p e c i e s , d i v i d e d by R (the u n i v e r s a l gas constant) and by the temperature, and added to the lo g a r i t h m of the pressure (1 atm.), s i n c e there i s only gas phase. c) Exponents matrix : For computing purposes , we added the nor m a l i t y c o n d i t i o n to the exponents matrix AA d e s c r i b e d i n the previous c h a p t e r . Hence the new matrix A i s a (4 x 6) matrix c o n s t r u c t e d i n the f o l l o w i n g way: A ( 1 , 1) = 1.DO A ( 1, j ) = O.DO j=2,N+1 A ( i , 1) = 0.DO i = 2,M+1 the remaining of each column accounts f o r the formula of each chemical s p e c i e s . The f i n a l matrix looks l i k e t h i s : dummy CO C02 H20 H2 CH4 1 0 0 0 0 1 1 0 A= 0 0 0 2 0 1 2 1 0 0 dummy 0 1 C 2 4 H 0 0 0 d) Vector of amount of elements B 84 Since we s t a r t e d with 2 moles of CH4 and 3 moles of H20, 2 atom-grams of C, 14 of H and 3 of 0 should be conserved throughout the r e a c t i o n s . Hence the B vecto r w i l l be: 1 dummy B= 2 C 1 4 H 3 0 e) s i n c e we only have one set of thermodynamic data, we are not concerned with s e n s i t i v i t y a n a l y s i s and C0MP2 and C0MP4 should not be a p p l i e d . f) we do not have a s t a r t i n g p o i n t to enter here. The programs w i l l generate one. g) i f we want to c o n s i d e r the formation of s o l i d C, then the c o n t r o l cards and the exponents matrix should be modified as f o l l o w s : For k=1, j var i e s between 1 and 6; for k = 2, j = 7 Dummy CO C02 H20 H2 CH4 C 1 0 0 0 0 0 0 Dummy 0 1 1 0 0 1 1 C 0 0 0 2 2 4 0 H 0 1 2 1 0 0 0 0 The B v e c t o r remains as be f o r e . The C matrix has to i n c o r p o r a t e the f r e e energy c o e f f i c i e n t f o r carbon (0) as C(7,1) 85 Example 2. Claus Furnace r e a c t i o n P r e l i m i n a r y data : - P a r t i a l l y r e a c t e d mixture of 0.1 moles S02, 0.2 moles H2S, 0.8 moles H20 and 1.88 moles of N2 ( i n e r t ) -Reach e q u i l i b r i u m composition at P = 1 atm., and 9 d i f f e r e n t T from 550 K to 1000 K (50 K i n t e r v a l s ) -Species present at e q u i l i b r i u m : S02, H2S, H20, S2, S4, S6, S8, N2. - I d e a l i t y assumed. -Elements : S, 0, H, N -One phase (gas) -Temperature c o e f f i c i e n t s f o r thermodynamic f u n c t i o n s of Gordon and Mc. B r i d e -Want to compare s e n s i t i v i t y a n a l y s i s versus r e - o p t i m i z a t i o n . I nput a) C o n t r o l c a r d s : N+1 = 9 (8 s p e c i e s ) M+1 = 5 (4 elements) K = 1 (one phase) NN = 9 (nine T,P c o n d i t i o n s ) NE= 2 ( f r e e energy c o e f f i c i e n t s should be c a l c u l a t e d ) NF =1 the s t a r t i n g p o i n t should be c a l c u l a t e d w i t h i n the program. b) matrix of Gordon and Mc.Bride c o e f f i c i e n t s . M atrix S i s an ( N x 7) matrix each row of the matrix corresponds to the seven 8 6 c o e f f i c i e n t s of the Gordon and Mc.Bride polynomials f o r a chemical s p e c i e s . The dummy s p e c i e s 1 i s not i n c l u d e d , so that the s p e c i e s here are s h i f t e d by one u n i t when compared to t h e i r numeration i n matrix A. With these data the subroutine FREEN w i l l c a l c u l a t e the matrix C. ' •-" c) Exponents matrix . Matrix A i s as f o l l o w s : dummy S02 H2S H20 S2 S4 S6 S8 N2 1 0 0 0 0 0 0 0 0 dummy 0 1 1 0 2 4 6 8 0 S 0 2 0 1 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 H 0 0 0 0 0 0 0 0 2 N d) Vector of elements content . The B ve c t o r i s as f o l l o w s : e) Matrix of s t o i c h i o m e t r i c c o e f f i c i e n t s : Since we want to apply s e n s i t i v i t y a n a l y s i s , we need to pro v i d e the matrix of s t o i c h i o m e t r i c c o e f f i c i e n t s U . I t i s a (N+1 x D) matrix ; each column c o n t a i n s the s t o i c h i o m e t r i c c o e f f i c i e n t of the spe c i e s d e s c r i b e d by the row number f o r the D independent r e a c t i o n s . To the dummy spe c i e s corresponds a s t o i c h i o m e t r i c c o e f f i c i e n t equal to zero. If M+1 i s the rank of the matrix A, then D=N - M. In t h i s case D= 4 . The U matrix i s b u i l t i n the f o l l o w i n g way : 1.0 0.3 1 .0 2.0 3.76 dummy S 0 H N 87 - i ) take matrix A, without the f i r s t row and column (they correspond to the n o r m a l i t y c o n d i t i o n of the geometric program). Exchange column 4 and 8, and p a r t i t i o n the r e s u l t i n g matrix so that 1 1 0 0 4 6 8 2 2 0 1 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 2 0 0 0 0 - i i ) work with elementary matrix o p e r a t i o n s to get 1 0 0 0 4/3 6/3 8/3 2/3 0 1 0 0 8/3 12/3 16/3 4/3 0 0 1 0 -4/3 -6/3 -8/3 -2/3 0 0 0 1 0 0 0 0 - i i i ) form the s t o i c h i o m e t r i c matrix. Interchange the rows 4 and 8 to get the o r i g i n a l set of s p e c i e s . -4/3 -6/3 -8/3 -2/3 -8/3 -12/3 -16/3 -4/3 4/3 6/3 8/3 2/3 0 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 To o b t a i n the U matrix f o r our programs , a f i r s t row of zeroes should be added. If we m u l t i p l y t h i s U matrix by 3, the set of independent r e a c t i o n s i s : S02 + 2 H2S = 3/j Sj + H20 j = 2, 4, 6, 8. Any other combination of r e a c t i o n s can be used, provided they 88 are l i n e a r l y independent. C a l c u l a t i o n of f r e e energy c o e f f i c i e n t s Subroutine FREEN c a l c u l a t e s the f r e e energy parameters of each s p e c i e s at a s p e c i f i e d T and P. It uses the polynomial approximation to the Gibbs f r e e energy determined by Gordon and Mc. B r i d e (1971). The polynomials are as f o l l o w s : > j 0 ( T ) / RT = S j , (1 - In T) - S j 2 T/2 - S j 3 T 2 /6 -Sj , T 3 /1 2 - Sj 5 T 4 ./20 + S j 6 A - S j 7 5.1 Where S i s a Nx7 matrix of the Gordon and Mc. Bri d e c o e f f i c i e n t s f o r the N s p e c i e s . Then the f r e e energy c o e f f i c i e n t s are c a l c u l a t e d through the d e f i n i t i o n s : C j l = /"-J(T) / RT + l n P 5.2 f o r gases and C j l = Gj(T) / RT 5.3 fo r condensed phases. The s u b s c r i p t 1 r e f e r s to the number of times T and P are v a r i e d . A l l the in f o r m a t i o n needed by the r o u t i n e , and produced by 89 i t , i s handled by means of l a b e l l e d COMMON b l o c k s . When other sources of fr e e energy data are a v a i l a b l e , the programs s k i p subroutine FREEN and read the matrix C s u p p l i e d by the user. Subroutine LIPSU2 . A dual s t a r t i n g p o i n t . If a good f i r s t guess of the composition i s a v a i l a b l e , the programs may take t h i s guess as a s t a r t i n g p o i n t . If t h i s i s not the case, the programs generate a dual s t a r t i n g p o i n t c a l l i n g the subroutine LIPSU2 . Th i s subroutine i s almost i d e n t i c a l to the program LIPSU2 e x p l a i n e d i n the previous chapter. The subroutine o b t a i n s the data from the main program through l a b e l l e d COMMON blocks and i t se t s the ta b l e a u i n the manner s p e c i f i e d by the UBC r o u t i n e LIPSUB . LIPSUB i s then c a l l e d and i t s o l v e s the LP using a p r i m a l - d u a l a l g o r i t h m . The r e s u l t s are s t o r e d i n the COMMON block FX. The main program then uses t h i s i nformation to c a l c u l a t e the dual s t a r t i n g p o i n t . Subroutine SINGV Subroutine SINGV transforms dual v a r i a b l e s i n t o p r i m a l ones. It i s very s i m i l a r to the program of the same name d e s c r i b e d on chapter 4 . The reading and output are re p l a c e d by l a b e l l e d COMMON b l o c k s . 90 Using GRG as a subroutine GRG i s c a l l e d from the main program as three s u b r o u t i n e s from UBC. GRGIN reads the c o n t r o l cards d e s c r i b e d i n chapter IV. GRG2 c a l l s the o p t i m i z a t i o n s u b r o u t i n e . The value of the o b j e c t i v e f u n c t i o n and of the v a r i a b l e s at the optimum are arguments of GRG2. GREG p r i n t s the Lagrange m u l t i p l i e r s of the b i n d i n g and e q u a l i t y c o n s t r a i n t s , and the c h a r a c t e r i s t i c s of the o p t i m i z a t i o n such as number of i t e r a t i o n s , number of f u n c t i o n e v a l u a t i o n s , e t c . G i v i n g input to GRG To read the c o n t r o l cards needed for the o p t i m i z a t i o n , subroutine GRGIN does i t from a f i l e a t t a c h e d to the l o g i c u n i t 5 . That i n t r o d u c e s two problems : a) The main program reads i t s own data from a d a t a f i l e attached to the input u n i t 5 . b) Some of the parameters that should be s p e c i f i e d in the c o n t r o l cards are c a l c u l a t e d in the f i r s t p art of the program, l i k e the s t a r t i n g p o i n t and the boundaries. Hence the c o n t r o l cards cannot be s e t t l e d i n advance. To solve these two p o i n t s , a s c r a t c h f i l e -DATA i s c r e a t e d from the main program. Input u n i t 5 i s r e a s s i g n e d to -DATA . The c o n t r o l cards are then w r i t t e n i n t o -DATA using the i n f o r m a t i o n a v a i l a b l e i n the main program , with the r e q u i r e d format. 91 Subroutine GCOMP GRG2 i s now ready to perform the o p t i m i z a t i o n . I t needs to evalu a t e the o b j e c t i v e f u n c t i o n and l e f t - h a n d s i d e s of the c o n s t r a i n t s f o r the given values of the v a r i a b l e s x . To do t h a t , i t c a l l s subroutine GCOMP. The form of t h i s subroutine w i l l depend on the mathematical f o r m u l a t i o n of the problem we are d e a l i n g with, but the v a r i a b l e s have to be c a l l e d x. Subroutine GCOMP gets the needed data through l a b e l l e d COMMON bl o c k s . For the dual case, the problem of logarithms with zero or negative arguments i s avoided by d e f i n i n g the l o g a r i t h m i c f u n c t i o n as equal to a very l a r g e negative number when negative or zero v a r i a b l e s are encountered i n the search. A l s o , the lower bound of the v a r i a b l e s i s set to a very low va l u e . A f t e r the o p t i m i z a t i o n i s f i n i s h e d , any number of moles smaller than the i n v e r s e of the Avogadro's number i s set to zero by the main program. In the pri m a l problem, the e x p o n e n t i a t i o n of very l a r g e negative numbers may occur, and t h i s i s a l s o undefined. We do not allow the lower boundaries of the pri m a l v a r i a b l e s to go below the d e f i n e d r e g i o n . To d e l e t e s p e c i e s from the model, the corresponding column of the exponent matrix i s set to zero. For the pr i m a l problem, that means that the molar f r a c t i o n of the corresponding s p e c i e s i s one. To av o i d the problem, we d e f i n e d DEXP(0)=0 when Aij=0. 92 Output from GRG The amount of output produced by GRG i s c o n t r o l l e d by the keyword PRINTCTL=number which should be s p e c i f i e d i n a c o n t r o l c a r d . When PRINTCTL=0, no i n f o r m a t i o n i s output. When i t i s set equal to 1, the r e s u l t s are p r i n t e d ; when working at a t e r m i n a l , the user i s prompted to w r i t e the word GO a f t e r the i n i t i a l c o n d i t i o n s f o r each o p t i m i z a t i o n are s e t . Succesive i n f o r m a t i o n on the p r o g r e s s i o n of the o p t i m i z a t i o n i s p r i n t e d when PRINTCTL values range from 2 to 4. The input l i n e s are echoed by d e f a u l t . If t h i s i s not d e s i r e d , NOECHO has to be s p e c i f i e d i n a c o n t r o l c a r d . The v a l u e s of the o b j e c t i v e f u n c t i o n and of the v a r i a b l e s at the optimum i s a l l the i n f o r m a t i o n needed from GRG when s o l v i n g the dual problem. They are arguments of GRG2 , hence f o r t h i s case we can a v o i d any p r i n t e d output from GRG i f d e s i r e d . When s o l v i n g the p r i m a l , we need a l s o the Lagrange m u l t i p l i e r s f o r the b i n d i n g c o n s t r a i n t s (they are the t o t a l number of moles f o r each c o n s t r a i n t when s o l v i n g the transformed p r i m a l problem), and they are not a v a i l a b l e as arguments. We had to get them through a r a t h e r t w i s t e d way. The main program c r e a t e s a s c r a t c h f i l e -GRGOUT, and the output u n i t 6 i s reassigned to t h i s f i l e . Subroutine GRGEG, from UBC, w r i t e s the output from the o p t i m i z a t i o n i n t o -GRGOUT, provided PRINTCTL i s set to a value d i f f e r e n t than zer o . Subroutine SKIP then scans through -GRGOUT u n t i l i t f i n d s the expr e s s i o n "Lagrange m u l t i p l i e r s " . At t h i s p o i n t , the main program reads the values of the m u l t i p l i e r s from the next l i n e 93 on the output s c r a t c h f i l e , and c a l c u l a t e s the number of moles f o r each s p e c i e s . -GRGOUT i s rewound a f t e r each o p t i m i z a t i o n . L i s t i n g the f i l e a l lows us to compute the number of i t e r a t i o n s and of f u n c t i o n e v a l u a t i o n s of the l a s t run. Performing s e n s i t i v i t y a n a l y s i s Subroutine NPLUSK i s very s i m i l a r to program NPLUSK d e s c r i b e d in chapter IV. The input of data i s done through l a b e l l e d COMMON b l o c k s . The s t o i c h i o m e t r i c c o e f f i c i e n t s have to be p r o v i d e d by the user, s i n c e we d i d not implement a subroutine to c a l c u l a t e them from the exponents matrix, as e x e m p l i f i e d e a r l i e r t h i s chapter ( s e c t i o n 5.2.) Subroutine NPLUSK formulates the l i n e a r system of equations from the data a v a i l a b l e . The r i g h t hand s i d e of the equations v a r i e s with pressure and temperature, and i s r e c a l c u l a t e d each time these c o n d i t i o n s change, but the l e f t hand sid e remains the same. Subroutine DSLIMP from UBC solves the system of e q u a t i o n s . Output from the programs If the keyword PRINTCTL i s set equal to 1, the output of the programs c o n t a i n s : a) Matrix of f r e e energy c o e f f i c i e n t s ( i f a p p l i c a b l e ) . b) Dual s t a r t i n g p o i n t . c) Primal s t a r t i n g p o i n t ( i f a p p l i c a b l e ) . d) Input f o r GRG 94 e) Pressure and temperature f) Number of moles of the s p e c i e s present at e q u i l i b r i u m at the prev i o u s P and T. g) T o t a l number of moles per phase at these c o n d i t i o n s . h) Corresponding molar f r a c t i o n s P o i n t s d) through h) are repeated f o r the d i f f e r e n t T and P v a l u e s . Point d) i s only repeated when there i s no s e n s i t i v i t y a n a l y s i s , f o r programs COMP1 and COMP3. At the end, the execution time i s p r i n t e d . 95 CHAPTER VI: RESULTS AND DISCUSSION The four computer programs d e s c r i b e d i n Chapter V were t r i e d on three problems : 5, 6 ,6B and 7 of. Appendix C. We c o u l d thus compare the perfomances of dual vs. p r i m a l f o r m u l a t i o n s , s e n s i t i v i t y a n a l y s i s vs r e - o p t i m i z a t i o n . The r e s u l t s are presented and d i s c u s s e d i n the present Chapter. T h i s Chapter i s organized i n four s e c t i o n s . In s e c t i o n one we t r y to evaluate the perfomance of the s t a r t i n g p o i n t r o u t i n e by comparing the p o i n t s generated by the r o u t i n e to the optimum v a l u e s . In the second s e c t i o n we compare the pri m a l and dual problems. In s e c t i o n three we compare s e n s i t i v i t y a n a l y s i s and r e - o p t i m i z a t i o n . In s e c t i o n four we de a l with the e f f e c t s of i n c o r p o r a t i n g new sp e c i e s to a chemical e q u i l i b r i u m model. E v a l u a t i o n of the s t a r t i n g p o i n t r o u t i n e If the user does not have a good f i r s t guess, a s t a r t i n g p o i n t i s c a l c u l a t e d w i t h i n our programs. The c a l c u l a t i o n of a dual s t a r t i n g p o i n t i s based on approximating the dual problem to a l i n e a r program. The method d e s c r i b e d by L i d o r (1975) i s used, and i t i s implemented i n subroutine LIPSU 2. An account of the method, and of the subroutine LIPSU 2 was made i n Chapters IV and V. When a pri m a l s t a r t i n g p o i n t i s needed, the dual s t a r t i n g p o i n t i s determined as above, and then i t i s transformed to a pri m a l p o i n t . Subroutine SINGV, d e s c r i b e d i n Chapter IV, performs t h i s task. 9 6 Table VI-1 shows the values of the s t a r t i n g p o i n t s c a l c u l a t e d by the programs f o r problems 5 and 6, and t h e i r optimum values . We w i l l now d i s c u s s the r e s u l t s of Table VI-1. a) The values of the o b j e c t i v e f u n c t i o n s at the s t a r t i n g p o i n t s are w i t h i n 2-4% of the optimal v a l u e s , which looks l i k e a q u i t e good approximation. b) The dual v a r i a b l e s of the s p e c i e s that are present i n g r e a t e r amounts at e q u i l i b r i u m are determined w i t h i n 10-20% of t h e i r optimal v a l u e . The r e s t of the v a r i a b l e s are wide a p a r t . The s t a r t i n g p o i n t r o u t i n e only g i v e s an estimate of as many s p e c i e s as the rank of the exponent matrix A, which should be M in a w e l l posed model. These s p e c i e s correspond to the b a s i c v a r i a b l e s of the l i n e a r program, and most of the time the dummy sp e c i e s introduced i n the l i n e a r program to account fo r non-l i n e a r i t i e s i s one of the b a s i c v a r i a b l e s . The r e s t of the s p e c i e s are given a f i x e d p o s i t i v e number, r e g a r d l e s s of t h e i r r e l a t i v e importance. So, i f the problem we are d e a l i n g with has one or two s p e c i e s present as t r a c e q u a n t i t i e s , the s t a r t i n g p o i n t does not allow us to d i s t i n g u i s h them from other secondary s p e c i e s . See problem 6 in t a b l e VI -1. As a consequence, the goodness of the s t a r t i n g p o i n t w i l l depend s t r o n g l y on the s p e c i f i c problem. c) The p r i m a l v a r i a b l e s f a l l w i t h i n 1.5% to 50% of t h e i r o ptimal v a l u e s . Again, the range i s a c h a r a c t e r i s t i c of the p a r t i c u l a r problem. d) The p r i m a l s t a r t i n g p o i n t s c a l c u l a t e d with t h i s method are 97 Table VI-1. E v a l u a t i o n of the s t a r t i n g p o i n t r o u t i n e s . Composition, o b j e c t i v e f u n c t i o n and pri m a l c o n s t r a i n t values at the s t a r t i n g p o i n t and at optimum f o r problems 5* and 6**. Prob. N. of moles Z G/RT P. const. S.P. Opt. S.P. Opt. S.P./ Opt. S.P./ Opt. 5* 3. 3. 3. 5.5 3. 3. 21.5 3. 5.25 3. .0204 .0007 .01 53 3.973 .0819 2.918 19.98 .0269 .0338 .0005 -1 0.0606 -22.3678 -13.5499 -11.581 9 -11.4153 -20.1267 -15.8052 -11.6971 746.3633 777.6382 230.6 1.00001 6** .2 .3 1 . .2 .2 .2 .2 1 .98 .0003 .0006 .9989 E-9 E-9 .0005 .0371 1 .880 -1.4607 -66.8466 -18.0296 -11.9268 -2.7906 -66.4363 -18.6591 -11.7611 148.1885 150.8133 49210. 1.0000 * D i n k e l and Lakshmanan, 1977. P= 40 a t . T = 2200 K ** Bonsu, 1981. P= 1 a t . T=355 K N. = Number Z = ln-t r a n s f o r m e d p r i m a l v a r i a b l e s G/RT = O b j e c t i v e f u n c t i o n S.P. = S t a r t i n g p o i n t Opt. = optimum P. Const. = pri m a l c o n s t r a i n t i n f e a s i b l e . In f a c t , the valu e s of the pri m a l c o n s t r a i n t s at the s t a r t i n g p o i n t s are q u i t e f a r apart from t h e i r optimal v a l u e s . GRG takes care of the i n f e a s i b i l i t y of the pri m a l s t a r t i n g p o i n t by o p t i m i z i n g the sum of the c o n s t r a i n t v i o l a t i o n s i n a "Phase I" procedure. We found that the combination of our s t a r t i n g p o i n t procedure p l u s GRG 's "Phase 98 I" produced b e t t e r r e s u l t s than a random p o i n t . The s t a r t i n g p o i n t r o u t i n e s performed w e l l i n problems 5,6,6B,7 (Appendix B) however, i f a b e t t e r s t a r t i n g guess i s a v a i l a b l e , the user i s encouraged to av o i d the s t a r t i n g p o i n t r o u t i n e . B e t t e r guesses shorten computation time • Pr i m a l - d u a l comparisons To compare the primal and dual f o r m u l a t i o n s of the chemical e q u i l i b r i u m problem,we used the codes COMP1 and COMP2 d e s c r i b e d i n the pr e v i o u s chapter. The problems t e s t e d were 5,6,6B, and 7 from appendix B. The value s of the o b j e c t i v e f u n c t i o n s and composition obtained are i n appendix B. The execution times are in t a b l e VI-2. . Table VI-2. Execution times ( i n sec.) of primal and dual codes fo r a f i x e d P and T. Problem Species Elements Primal Dual (sec . ) (sec .) 5 10 4 .1265 .1612 6 8 4 .1342 .1695 6b 8 4 .1355 .1720 7 24 4 .5886 f a i l e d To e v a l u a t e the r e l a t i v e e f f e c t i v e n e s s of programming codes we need some c r i t e r i a . Himmelblau (1972) proposed and d i s c u s s e d s e v e r a l e v a l u a t i o n c r i t e r i a . We w i l l repeat them here, and 99 then we w i l l t r y to see how these c r i t e r i a apply to our problem. The c r i t e r i a are as f o l l o w s : i ) s i z e of the problem, i i ) time r e q u i r e d to introduce data i n t o the program , i i i ) s i m p l i c i t y of the computer program, i v ) success i n s o l v i n g r e a l world problems most of the time, v) accuracy of the r e s u l t s , v i ) time. The s i z e of the problem favours the p r i m a l over the dual f o r m u l a t i o n f o r the chemical e q u i l i b r i u m problem with many r e a c t i o n s . The time r e q u i r e d to introduce data i n t o the program i s the same f o r both the pri m a l and dual problems. The dual program i s more simple to implement than the p r i m a l , because i t does not need subroutines SINGV and SKIP . The p r i m a l c o u l d be s i m p l i f i e d i f the subroutine GRG 2 would have the Lagrange m u l t i p l i e r s f o r the p r i m a l c o n s t r a i n t s as arguments, thus subroutine SKIP can be avoided. As f o r " s o l v i n g r e a l - l i f e problems most of the time", the pr i m a l always d i d so. The dual f a i l e d when t r y i n g a problem with 24 v a r i a b l e s (problem 7). The problem had many t r a c e v a r i a b l e s , which were sent to the lower boundaries a f t e r a few i t e r a t i o n s and would not move t h e r e a f t e r . Making the values of the boundaries smaller d i d not improve the s i t u a t i o n . A new computer code, r e c e n t l y a v a i l a b l e at UBC , MINOS, i s based on an a l g o r i t h m s i m i l a r to GRG f o r the case of l i n e a r e q u a l i t y c o n s t r a i n t s . Murtagh and Saunders (1978) c l a i m that MINOS so l v e d a badly s c a l e d dual chemical e q u i l i b r i u m problem of 45 v a r i a b l e s . We d i d not t r y t h i s code ,though. The s t a r t i n g p o i n t and the r e q u i r e d accuracy of the s o l u t i o n i n f l u e n c e g r e a t l y on the computation time. Since the p r i m a l 100 s t a r t i n g p o i n t i s c a l c u l a t e d from the d u a l , we are almost sure that our comparison between the pr i m a l and the dual f o r m u l a t i o n s w i l l not be a f f e c t e d by the s t a r t i n g p o i n t . T h i s i s only a good approximation because the r e l a t i o n s between p r i m a l and dual v a r i a b l e s that we used i n subroutine SINGV only h o l d as e q u a l i t i e s at the optimum p o i n t . As f o r the accuracy of the s o l u t i o n , we used the same parameters i n GRG f o r the pr i m a l and dual problems. The r e s u l t s were comparable with the l i t e r a t u r e , and between the two codes, w i t h i n .2%. For the pri m a l problem the accuracy f o r the composition values was determined by the t o l e r a n c e of the c o n s t r a i n t ( EPNEWT ). The p r i m a l c o n s t r a i n t s s t a t e t hat the sum of the molar f r a c t i o n of the spe c i e s present at a p a r t i c u l a r phase should be one at e q u i l i b r i u m . A value of EPNEWT of 10 " 6 i m p l i e s that the absolute e r r o r s of the molar f r a c t i o n s of each s p e c i e s w i l l be of t h i s order of magnitude. That means that the r e l a t i v e e r r o r s f o r the s p e c i e s present i n t r a c e c o n c e n t r a t i o n s w i l l be q u i t e important. We can change the value of EPNEWT to improve the accuracy; i n f a c t we t r i e d up to a value of 1 0 " 1 0 f o r problem 7. S t i l l , f o r t h i s p a r t i c u l a r problem, the c o n c e n t r a t i o n of the sp e c i e s d i f f e r as much as 1 0 2 2 . We do not think i t i s p o s s i b l e to o b t a i n such accuracy at the moment; we are al r e a d y working with double p r e c i s i o n . Besides, we are i n t e r e s t e d on a code that may solve a wide v a r i e t y of examples, r a t h e r than produce the best s o l u t i o n f o r a p a r t i c u l a r problem. The in c r e a s e i n accuracy a l s o means i n c r e a s i n g the execution time. We have shown that our pri m a l and dual codes both s t a r t from 101 approximately the same po i n t and get the same r e s u l t s . The time i n v e s t e d i n that procedure i s an important c r i t e r i o n of the e f f e c t i v e n e s s of a code. Time r e f e r s e i t h e r to the number of f u n c t i o n e v a l u a t i o n s or to the execution time. In our case, there i s no p o i n t in comparing the number of f u n c t i o n e v a l u a t i o n s , s i n c e the mathematical f o r m u l a t i o n s f o r the problems are q u i t e d i f f e r e n t . We have to compare execution times. The e x e c u t i o n times i n t a b l e VI-2 are not s t a n d a r i z e d , s i n c e we were concerned only with the comparison between our codes. Care was taken to run a set of p r i m a l - d u a l problems one a f t e r the o t h e r , i n order to av o i d t i m e - s h a r i n g problems. The p r i m a l f o r m u l a t i o n proved to be between 23 and 39% f a s t e r than the dual f o r the case of m i d d l e - s i z e problems (8-10 s p e c i e s ) . S e n s i t i v i t y a n a l y s i s and r e - o p t i m i z a t i o n To be abl e to evaluate the performance of the s e n s i t i v y a n a l y s i s , we compared the codes COMP1 and COMP2 (primal with and without s e n s i t i v i t y a n a l y s i s ) and COMP3 and COMP4 (dual with and without s e n s i t i v i t y a n a l y s i s ) . We t e s t e d the codes on problems 5 and 6. Problem 5 had been used i n the l i t e r a t u r e as an example f o r s e n s i t i v i t y a n a l y s i s ( D i n k e l and Lakshmanan, 1977), but only the accuracy of the method had been d i s c u s s e d . We repeated t h e i r composition v a l u e s (see chapter I V ) . The f r e e energy data f o r problem 6 are c a l c u l a t e d from the Mc Bri d e and Gordon c o e f f i c i e n t s ( NASA , 1971) as a f u n c t i o n of temperature. That allowed us to vary the temperature at small i n t e r v a l s , and 1 02 to compare the e f f e c t of these v a r i a t i o n s on the accuracy of the r e s u l t s and on the execution times. Table VI-3. Execution times for s e n s i t i v i t y a n a l y s i s (S.A.) and r e - o p t i m i z a t i o n (R.O) Problem Species N. of T R. 0. S. A. Form. S t . Elements P, T (K) T i m e ( s e c ) T i m e ( s e c ) cond. 5 10/4 5 - .9673 .4357 D .44 5 10/4 5 - .7299 .2885 P .50 6 8/4 9 50 1.4361 .3144 P .51 6 8/4 9 50 1.8665 . 4661 D .44 6 8/4 9 1 0 1.2710 .3011 P .47 6 8/4 9 5 1.2243 .3235 P .42 6 8/4 9 1 1.0373 .3226 P .35 T i n t e r v a l s of temperature between two o p t i m i z a t i o n s , at P=1at. fo r problem 6. in problem 5, T and P v a r i e d . P= p r i m a l D= dual Form.= f o r m u l a t i o n cond.= c o n d i t i o n s . St= (R.O time / SA. t i m e ) / n. of P, T c o n d i t i o n s ) We w i l l use the same comparison c r i t e r i a mentioned b e f o r e . The s i z e of the problem i s s m a l l e r f o r s e n s i t i v i t y a n a l y s i s , s i n c e i t s o l v e s a system of l i n e a r equations .instead of using a n o n l i n e a r o p t i m i z a t i o n method. S e n s i t i v i t y a n a l y s i s needs e i t h e r more time to i n t r o d u c e data, or more programming, because of the s t o i c h i o m e t r i c c o e f f i c i e n t s f o r the r e a c t i o n s . Re-o p t i m i z a t i o n i s e a s i e r to program : i t i s only q u e s t i o n of adding a DO loop to the o p t i m i z a t i o n program. Let us take a look at the accuracy of the s e n s i t i v i t y 103 a n a l y s i s . We w i l l r e f e r the composition values obtained by s e n s i t i v i t y a n a l y s i s to the values c a l c u l a t e d by the o p t i m i z a t i o n procedure. The l a t t e r compare w e l l with the l i t e r a t u r e . We w i l l d e f i n e now % r e l a t i v e e r r o r of sp e c i e s j as: % € j = { ( 6 j * - 6j 1 ) / 6j *} 100 6.1 Where, %€: : percentage e r r o r i n the dete r m i n a t i o n with s e n s i t i v i t y a n a l y s i s of the number of moles of species j , r e f e r r e d to r e - o p t i m i z a t i o n 6 j 1 : number of moles of s p e c i e s j c a l c u l a t e d by s e n s i t i v i t y a n a l y s i s 6 j * : number of moles of s p e c i e s j c a l c u l a t e d by r e -o p t i m i z a t i o n ) The values c a l c u l a t e d by r e - o p t i m i z a t i o n are, of course, s u b j e c t to e r r o r s . We have d i s c u s s e d i n the p r e v i o u s s e c t i o n that the ab s o l u t e e r r o r i n the dete r m i n a t i o n of the molar f r a c t i o n s by o p t i m i z a t i o n was around 10~ 6. That means a c o n s i d e r a b l e r e l a t i v e e r r o r f o r s p e c i e s present as t r a c e s . In f i g u r e 1, we p l o t t e d % e r r o r vs v a r i a t i o n of temperature fo r problem 6. The s p e c i e s number 4 was present i n small q u a n t i t i e s . Changing the temperature 10K introduced an e r r o r of 18% when using s e n s i t i v i t y a n a l y s i s as compared to a new F i g u r e 1. % E r r o r v s . V a r i a t i o n o f T e m p e r a t u r e F o r P r o b l e m 6 105 o p t i m i z a t i o n . I f we c a l c u l a t e the e r r o r f o r a s p e c i e s present in g r e a t e r q u a n t i t i e s , ( l i k e the summation of s p e c i e s 4 to 8, which corresponds to t o t a l sulphur) the value of the e r r o r w i l l be around .2% , s i m i l a r to the u n c e r t a i n t y of our method. For a v a r i a t i o n of 100K, the e r r o r f o r s p e c i e s 4 i s around 95%, while the summation of spe c i e s 4 to 8 only has 4% e r r o r . The v a l u e s mentioned above f o r the e r r o r s are v a l i d only f o r t h i s s p e c i f i c problem. The f r e e energy c o e f f i c i e n t s vary i n a d i f f e r e n t way f o r d i f f e r e n t s p e c i e s i n d i f f e r e n t thermodynamical c o n d i t i o n s . However, we can assume that the magnitudes of the e r r o r s are more i n f l u e n c e d by the r e l a t i v e importance of the sp e c i e s than by the r e l a t i v e changes of the f r e e energy c o e f f i c i e n t s . T h i s i s shown i n Table VT-4 . We can see there that f o r a s i m i l a r v a r i a t i o n i n the f r e e energy c o e f f i c i e n t s (water and S02) the s p e c i e s with lower c o n c e n t r a t i o n (S02) has 33% e r r o r when determined with s e n s i t i v i t y a n a l y s i s , vs. .2% e r r o r f o r the water. D i v i d i n g the temperature v a r i a t i o n i n smaller i n t e r v a l s doesn't seem to have great e f f e c t when the t o t a l v a r i a t i o n of the fr e e energy c o e f f i c i e n t s i s below 10%. It does reduce the e r r o r of S2, which has a t o t a l v a r i a t i o n of the Cj of 41.6%. Thus, as f a r as accuracy i s concerned, we conclude that the s e n s i t i v i t y a n a l y s i s method should only be used when we are not concerned with the accuracy of the sp e c i e s present i n small q u a n t i t i t e s , and when the v a r i a t i o n of the f r e e energy c o e f f i c i e n t s are of the order of 10% or l e s s . Table VI-3 shows the execution times f o r s o l v i n g problems 5 and 6 with the four programs COMP1, COMP2, COMP3 and COMP4 106 Table VI"4. S e n s i t i v i t y a n a l y s i s . E f f e c t of the number of temperature increments on the composition, f o r a t o t a l v a r i a t i o n of temperature of 50 K. Problem 6* Spec i e s x j % C j * * % E r r o r (re j f e r r e d to opt: imization) 1 i n c . 5 i n c . 10 i n c . N2 H20 S02 S2 .64411 .34154 .00049 2.e-8 .54 8.1 7.9 41 .6 .0001 .23 33.0 85.2 .00001 .24 30.4 74. 1 .00001 .12 30.2 73.0 * Bonsu, 1981 % C j * * % v a r i a t i o n of the f r e e energy c o e f f i c i e n t s f o r each s p e c i e s , from 410 K ( b a s i s ) to 460 K. i n c . i n c r e m e n t s of temperature. Xj= molar f r a c t i o n at e q u i l i b r i u m d e s c r i b e d i n the previous s e c t i o n . To t r y to account f o r the d i f f e r e n t number of P and T c o n d i t i o n s c o n s i d e r e d i n the two problems, we d e f i n e d a r a t i o St as f o l l o w s : St = (time f o r R.O. / time f o r S.A)/ n.of P,T c o n d i t i o n s where R.O. stands for r e - o p t i m i z a t i o n , and S.A. f o r s e n s i t i v i t y a n a l y s i s . The St r a t i o g i v e s a rough idea of the r e l a t i v e speed of the s e n s i t i v i t y a n a l y s i s method . We can thus say, very approximately, that the s e n s i t i v i t y a n a l y s i s method i s between 35% to 50% f a s t e r than the o p t i m i z a t i o n per number of d i f f e r e n t thermodynamic c o n d i t i o n s . For problem 6 we see that the times f o r the p r i m a l r e -1 07 o p t i m i z a t i o n s are sh o r t e r f o r sm a l l e r v a r i a t i o n s of temperature. T h i s can be e x p l a i n e d by the f a c t t h a t , f o r any new o p t i m i z a t i o n , GRG uses the p r e v i o u s optimal v a l u e s as s t a r t i n g points.- When the temperature d i f f e r e n c e s are s m a l l , the optimal values vary very l i t t l e . So f a r , the d e c i s i o n of i n c o r p o r a t i n g a subroutine to perform s e n s i t i v i t y a n a l y s i s to the f i n a l program, w i l l depend on a trade between l o s s of accuracy and gain of speed. The l o s s of accuracy seems q u i t e b i g , and i t i s not e a s i l y p r e d i c t e d . The gain of speed i s not very i m p r e s s i v e . F i n a l l y , one comment. So f a r we used the s e n s i t i v i t y a n a l y s i s r o u t i n e with s t o i c h i o m e t r i c c o e f f i c i e n t s provided by the user. To in c r e a s e the s i m p l i c i t y of the program to the user, the code should a l s o have a subroutine to c a l c u l a t e these s t o i c h i o m e t r i c c o e f f i c i e n t s , u s i n g the f i r s t o p t i m i z a t i o n r e s u l t s to chose the key components. But t h i s c a l c u l a t i o n w i l l need some time; and the d i f f e r e n c e s i n CPU times between s e n s i t i v i t y a n a l y s i s and r e - o p t i m i z a t i o n w i l l be s m a l l e r . E f f e c t s of i n c o r p o r a t i n g new s p e c i e s to a chemical e q u i l i b r i u m model Problems 6B and 7 are c l o s e l y r e l a t e d . Problem 6B d e s c r i b e s the o x i d a t i o n of H2S i n a Claus r e a c t o r with 100% s t o i c h i o m e t r i c a i r . E i g h t s p e c i e s are assumed to be present at e q u i l i b r i u m . Problem 7 i s the. same case, with 24 sp e c i e s p r e s e n t . The r e s u l t s of s o l v i n g the p r i m a l of both problems are shown i n Appendix B f o r one value of temperature and Pressure. The 108 r e s u l t s agree with the l i t e r a t u r e in each case. A l s o , the t h e o r e t i c a l composition was checked a g a i n s t experimental values i n the l i t e r a t u r e . We observe the same d i s c r e p a n c i e s , s i n c e our problems agree with t h e i r t h e o r e t i c a l v a l u e s . Both problems, 7 and 6B, have the same q u a n t i t y of pr i m a l Table VT-5. E f f e c t s of adding new sp e c i e s to a model i n the o b j e c t i v e f u n c t i o n and in the t o t a l number of moles, f o r d i f f e r e n t temperatures. Temperature (K) Problem 6b Problem 7 G/RT ' X. G/RT X. 600 247.24402 6. 09281 247.17428 6. 08926 650 241.51526 6. 1 3869 241.46424 6. 1 3431 700 236.86576 .6. 20562 236.83423 6. 20244 750 233.11049 6. 30549 233.10187 6. 31 000 800 230.12108 6. 44289 230.13785 6. 45207 X : T o t a l number of moles (1 phase, gas) G/RT: O b j e c t i v e f u n c t i o n . v a r i a b l e s . So, f o r the pri m a l problems, we have the same o b j e c t i v e f u n c t i o n but the c o n s t r a i n t w i l l have t r i p l e d the number of terms. For the dual problem, there are three times more v a r i a b l e s i n problem 7, but we w i l l have the same number of c o n s t r a i n t s as bef o r e . We a l r e a d y mentioned that GRG f a i l e d when t r y i n g t h i s dual problem. Table VT-5 shows how the o b j e c t i v e f u n c t i o n and the t o t a l number of moles change with temperature f o r both problems. We 109 Table VI-6. Adding new sp e c i e s to a model: e f f e c t s on the execution time Problem Execution time * Execution time* (sec.) Tolerances= E-6 (sec.) Tolerances= E-10 6b 1.6384 7 4.2558 5.9839 6 1.4361 7 mod** 4.1532-* Over 9 d i f f e r e n t P and T c o n d i t i o n s ** Problem 7 was m o d i f i e d to have the same amount of elements as problem 6. can see that the values are s i m i l a r ; the o b j e c t i v e f u n c t i o n of problem 7 i s smal l e r than that of problem 6B. F i g u r e s 2 and 3 present the v a r i a t i o n with T of the composition of the two problems at a range of temperatures. By watching the two f i g u r e s , we can see that the s p e c i e s present in bigger q u a n t i t i e s remain in approximately the same composition i n both cases. The main d i f f e r e n c e s arose from the s u l f u r compounds. I t i s not only a q u e s t i o n of not i n c l u d i n g s p e c i e s ; the free energy c o e f f i c i e n t s f o r S4 was d i f f e r e n t i n the two examples. The r e s t of the s p e c i e s were not important i n the range of temperatures c o n s i d e r e d . The computation times f o r nine d i f f e r e n t values of P and T are presented i n t a b l e VI-6. The execution time for the problem with 24 s p e c i e s i s 30% higher than f o r the problem with 8 s p e c i e s , based on a "per o p t i m i z a t i o n " b a s i s . When s o l v i n g problem 7 on the same c o n d i t i o n s as above, but i n c r e a s i n g the accuracy of the pri m a l F i g u r e 3. C o m p o s i t i o n vs. Temperature For Problem 7 T C K ) 1 12 Table VI-7. Adding new s p e c i e s to a model: e f f e c t s on the composition. Problems 6 and 7, same i n i t i a l composition, T=600K and P=1 atm. Spec i e s Molar f r a c t i o n s Spec i e s Molar f r a c t ions Prob. 6 Prob.7 Prob.6 Prob.7 S02 .0112133 .0116286 H2 - .0000069 H2S .0224051 .0232480 H - E-1 9 H20 .3173050 .3163954 SO - E-1 1 S2 .0003625 .0003916 OH - E-20 S4 .0000480 E-22 S03 - E-1 3 S6 .0057975 .0034403 SN - E-1 9 S8 .0040892 .0040944 S20 - .0000425 N2 ..6387804 .6387978 NO - E-22 NH3 - E-1 0 S3 - .0000425 S - E-1 5 S5 - .0003127 SH - E-1 2 S7 - .0016156 0 - 0 02 - 0 prob. 6 : problem 6 of Appendix B. prob. 7 : problem 7 of Appendix B. c o n s t r a i n t from 10" 6 to 10 " 1 0 , the execution time i n c r e a s e d to 5.9839, 16% more on a per o p t i m i z a t i o n b a s i s . Since problem 6 and 6B only d i f f e r e d on the amount of a i r used (the B v e c t o r of amount of elements), we compared problem 6 with problem 7 mo d i f i e d so that the i n i t i a l c o n d i t i o n s of .problem 7 were the same as f o r problem 6. The r e s u l t s are in t a b l e VI-7 f o r T=600K and P= 1at. We observed the same type of v a r i a t i o n s as with problems 6B and 7 (unmodified): s l i g h t l y l e s s value of the o b j e c t i v e f u n c t i o n when more s p e c i e s are added to the model ; no a p p r e c i a b l e v a r i a t i o n on the composition of the s p e c i e s present i n g r e a t e r amounts. 1 1 3 CHAPTER V I I : CONCLUSIONS AND RECOMMENDATIONS Conclusio n s When s o l v i n g the chemical e q u i l i b r i u m problem with the computer code G e n e r a l i z e d Reduced Gr a d i e n t , we found: 1) GRG cannot so l v e the pr i m a l geometric programming fo r m u l a t i o n of the problem d i r e c t l y . E i t h e r a s c a l i n g of the v a r i a b l e s and of the o b j e c t i v e f u n c t i o n ,or a l o g a r i t h m i c t r a n s f o r m a t i o n has to be performed before using GRG. We found that GRG needs l e s s i t e r a t i o n s and f u n c t i o n e v a l u a t i o n s to s o l v e the l o g a r i t h m i c -transformed problem, than to sol v e the s c a l e d v e r s i o n . We worked with the transformed p r i m a l t h e r e a f t e r . 2) When s o l v i n g the dual geometric programming with GRG, the o b j e c t i v e f u n c t i o n needs to be s c a l e d , or the lo g a r i t h m of the o b j e c t i v e f u n c t i o n should be minimized. We d i d not f i n d any p a r t i c u l a r computational advantage f o r any of the two approaches, but we chose the l a t t e r to compare the r e s u l t s to those of the transformed p r i m a l . 3) In both the pri m a l and the dual f o r m u l a t i o n s , the o b j e c t i v e f u n c t i o n i s not very s e n s i t i v e to changes i n the v a r i a b l e s . Two val u e s of the o b j e c t i v e f u n c t i o n may d i f f e r in t h e i r n i n t h s i g n i f i c a n t f i g u r e , yet the v a r i a b l e s d i f f e r i n t h e i r s i x t h s i g n i f i c a n t f i g u r e . 1 1 4 4) The accuracy of the p r i m a l s o l u t i o n depends mainly on the t o l e r a n c e s f o r the c o n s t r a i n t s . At the optimum, each c o n s t r a i n t accounts f o r the summation of the molar f r a c t i o n s of each s p e c i e s present i n the phase d e s c r i b e d by the c o n s t r a i n t . The ab s o l u t e e r r o r of the c o n s t r a i n t equals to the summation of the ab s o l u t e e r r o r of each term. Hence any s p e c i e s whose molar f r a c t i o n i s smaller than the t o l e r a n c e of the c o n s t r a i n t w i l l be su b j e c t to e r r o r s . 5) The accuracy of the dual v a r i a b l e s depends on the t o l e r a n c e s f o r the c o n s t r a i n t s ( c o n s e r v a t i o n of elements) and on the accuracy of the t o t a l amount of each element. 6) GRG so l v e s the pri m a l problem 30% f a s t e r than the d u a l , f o r problems with s i x r e a c t i o n s . GRG f a i l e d to so l v e the dual of a l a r g e problem, with twenty r e a c t i o n s . Adding new sp e c i e s to a model i s more e a s i l y done with the pri m a l f o r m u l a t i o n . The only advantage of the dual seems to be that i t i s more easy to implement than the p r i m a l . We concluded that the pri m a l i s s u p e r i o r to the dual when both are sol v e d with the GRG code. 7) The method to perform s e n s i t i v i t y a n a l y s i s suggested by D i n k e l and Lakshmanan (1977) proved a thousand times slower than the method d e r i v e d by R i j c k a e r t (1974). Both methods were t r i e d on a chemical e q u i l i b r i u m problem with s i x r e a c t i o n s . We used the R i j c k a e r t ' s method from then on. 1 15 8) S e n s i t i v i t y a n a l y s i s was between 30% to 50% f a s t e r than r e o p t i m i z a t i o n , per number of pressure and temperature changes. However, the r e s u l t s obtained with s e n s i t i v i t y a n a l y s i s showed d i s c r e p a n c i e s when compared to the o p t i m i z a t i o n ones. " The e r r o r s were more s e r i o u s f o r the l e s s important s p e c i e s . ' I t i s our o p i n i o n that the computational speed gained with s e n s i t i v i t y a n a l y s i s does not compensate f o r the l o s s of accuracy of the r e s u l t s . Recommendations Based on the c o n c l u s i o n s s t a t e d above, i t i s recommended the use of the code COMP1, l i s t e d on Appendix C. I t i s a p r i m a l based code, and i t does not perform s e n s i t i v i t y a n a l y s i s . I t i s a l s o recommended that f u r t h e r r e s e a r c h should be c a r r i e d on with a b e t t e r implementation of GRG f o r h a n d l i n g l a r g e dual problems. Murtagh and Saunders (1978) c l a i m that t h e i r code MINOS can do so. 1 16 REFERENCES CITED A r i s , R. 1970. American S c i e n t i s t , 58:419-427. B a l l a r d , D.H. J e l i n e k , C O . and S c h i n z i n g e r , R. 1974 . Computer J . 22:261-266. n.3 Beck, P.A. and Ecker, J.G. 1975 . J.O.T.A. Jj>: 189-201 . n.2 B e i g h t l e r , C.S. and P h i l l i p s , D.T. 1976. " A p p l i e d Geometric Programming " , John Wiley and Sons, N.Y. Bennett, H.A. 1979. Ph.D. T h e s i s , U n i v e r s i t y of B r i t i s h Columbia, Dept. of Chem. Eng. Vancouver, B.C. Bigelow, J.H. and Shapiro, N.Z. 1971. RAND C o r p o r a t i o n Paper P-4628 B i r d , C. 1979. "UBC-MATRIX. A Guide To S o l v i n g M a t r i x Problems", UBC Computing Centre, Vancouver, B.C. Blau, G.E. and Wilde, D.J. 1971 . AIChE J o u r n a l , 17:235-240. n. 1 Bonsu, A.K. 1981. Ph.D. T h e s i s , U n i v e r s i t y of B r i t i s h Columbia, Dept. of Chem. Eng. Vancouver, B.C. Bradley, J . 1973. "An Al g o r i t h m For The Numerical S o l u t i o n of Prototype Geometric Programming", Report of the I n s t i t u t e Of I n t e r n a t i o n a l Research And Standards, D u b l i n , I r e l a n d . B r i n k l e y , S.R. 1947. J . Chem. Phys. J_5: 107-1 10. Buzby,.B.R. 1974. "Techniques And Experience S o l v i n g R e a l l y Big Non-linear Programs" i n " O p t i m i z a t i o n Methods", C o t t l e , R.W. and Krarup, J . Ed. E n g l i s h Univ. Press, London. 1 1 7 C a v a l l o t t i , P. C e l e r i , G. and L e o n a r d i s , B. 1980. Chem. Eng, S c i . 35:2297-2304. Clasen, R.J, 1965. RAND C o r p o r a t i o n Memo RM-4345-PR (January, 1965) Dembo, R.S. 1 976 . Math. Prog. J_0 : 1 92-21 3 . Dembo, R.S. 1978 . J.O.T.A. , 26:149-183. n.2 Dembo, R.S. 1978 . J.O.T.A. , 26:243-252. n.2 Dembo, R.S. 1979 . Math. Prog. J_3: 156-175. Denbigh, K. 1966. "The P r i n c i p l e s Of Chemical E q u i l i b r i u m " , Cambridge U n i v e r s i t y Press, 493 pp. D i n k e l , J . and Kochenberger, G. 1977. Op. Res. 25:155-163. D i n k e l , J . and Kochenberger, G. 1979 . Math. Prog. J_7:109-1 13. D i n k e l , J . and Lakshmanan, R. 1975. J . of Eng. Math. 9:343-352, D i n k e l , J . and Lakshmanan, R. 1977. Computers and Chemical Eng. 1 -.41-47. D i n k e l , J . E l l i o t , W. and Kochenberger, G. 1977 . Math. Prog. 13:200-220. D i n k e l , J . Kochenberger, G. and Mc.Carl, B. 1974 . Math. Prog 7:181. D u f f i n , R. Peterson, E. and Zener, C. 1967. " Geometric Programming. Theory And A p p l i c a t i o n . " , John Wiley and Sons, N.Y. , 278 pp. 118 D u f f i n , R.J. And Zener, C. 1969. Proc. of Nat. Ac. of S c i . 63: 629-636. D u f f i n , R.J. and Zener, C. 1970. J . Phys. Chem. 74:2419-2423. n. 1 2 Ecker, J.G. 1972 . J.O.T.A. 9:176-181. n.3 G a r c i a , A. and Hogg, G. 1973 . "Computer Code For The Blau A l g o r i t h m " i n "Optimizing With FORTRAN ", McGraw H i l l , N.Y. pp 136-154. Gautam, S. and S e i d e r , W.D. 1979. AIChE J o u r n a l 25:991-998. n. 6. Gochet, W. and Smeers, Y. 1974 . C a h i e r s Du Centre D'Etudes De Recherche O p e r a t i o n n e l l e , 16:23-36. Gochet, W. Loute, E. Solow, D. 1974 . C a h i e r s Du Centre D'Etudes De Recherche O p e r a t i o n n e l l e , 16:469-486. Gordon, S. and McBride, B.J. 1971. NASA SP- 273, Washington, D.C. G r i f f i t h , R.E. and Stewart, R. A. 1961. Management Science, 7_: 379-381 . Haarhoff, P.C. and Buys, J.D. 1970. Computer J o u r n a l , 13:178-184. Himmelblau, D.M. 1972. "Applied Non L i n e a r Programming", McGraw-H i l l ed. Kochenberger, G. Woolsey, R. and McCarl, B. 1973 . Op. Res. Q. 24:285-293. n.2 Lasdon, L.S. 1981. "A Survey Of Nonlinear Programming Algorithms And Software", i n "Foundations Of Computer Aided Chemical Process Design", Mah, R. S.H. and Se i d e r , W. D. E d i t o r s , Pub. by the Nat. S c i . Foundation, NY,NY. pp.185-217. 119 L i d o r , G. and Wilde, D.J. 1978. J.O.T.A. 26:77-96. n.1 L i d o r , G. 1975. Ph.D T h e s i s , S t a n f o r d U n i v e r s i t y , Dept. Of Operations Research, Tech. Report 75-8. McBride, B.J. Heimel, S. E h l e r s , J.G. and Gordon, S. 1963. NASA SP-3001, Washington, D.C. McGreggor, D.E. 1971. Ph.D. T h e s i s , U n i v e r s i t y Of A l b e r t a , Dept. Of Chemical And Petroleum Eng. Edmonton, A l b e r t a . McNamara, J.R. 1976 . Op. Res. 2^:15-25. n.1 Murtagh, B.A. and Saunders, M. A. 1978. Math. Programming, 14: 41-72. Passy, U. and Wilde, D. 1968. SIAM J . Appl. Math. j_6:363-373. n.2 Passy, U. and Wilde, D. 1969. J . of Eng. Math. 3:325. n. 4 Passy, U. 1972. J.O.T.A. 9:221-234. n.4 Pa t t e r s o n , M. 1981. "UBC-LIPSUB. A L i n e a r Programming Subroutine", UBC Computing Centre, Vancouver, B.C. Perry, R.H. and C h i l t o n , C H . 1973. "Chemical E n g i n e e r i n g Handbook", McGraw-Hill Chemical Eng. S e r i e s , 5th. Ed. P h i l l i p s , D.T. 1974. AIIE T r a n s a c t i o n s 6:114-119. n.2 Ratner, M. Lasdon, L.S. and J a i n , A. 1978. J.O.T.A. 26:253-263. n. 2 R i c h t e r , E. and Hofmann, H. 1976. I n t e r n a t . Chem. Eng. 16:137-143. 120 R i j c k a e r t . M. 1974. " S e n s i t i v i t y A n a l y s i s In Geometric Programming" i n "Mathematical Programming For A c t i v i t y A n a l y s i s " , Ed. P. Van Moeseke, North H o l l a n d Pub.Co. Amsterdam, H o l l a n d , pp. 61-72 R i j c k a e r t , M. and Martens, X.M. 1976 . Math. Prog. J J _ : 8 9 - 9 3 . R i j c k a e r t , M.J. and Martens, X.M. 1978 . J.O.T.A., 26:325-340. n. 2 Sarma, P.V. Martens, X.M. R e k l a i t i s , G.V. and R i j c k a e r t M.J. 1978 . J.O.T.A. 26:185-242. n.2 Schneider, D.R. and R e k l a i t i s , G.V. 1975. Chem. Eng. S c i . 30: 243-247. Schubert, E. and Hofmann, H. 1976. I n t e r n a t . Chem. Eng. 16:132-136. n.1 S c u l l y , D.B. 1962. Chem. Eng. S c i . r7 :977-985. Shapiro, N.Z. 1964. RAND C o r p o r a t i o n Report RM-4128-PR Shapley, M. and C u t l e r , L. 1970. RAND C o r p o r a t i o n Report R-495-PR (June, 1971) Van Zeggeren, F. and Storey, S.H. 1970. "The Computation Of 1 Chemical E q u i l i b r i a " , Cambridge U n i v e r s i t y Press, 176 pp. Wales, K. 1977. "UBC-GRG. G e n e r a l i z e d Reduced Gradient Program", UBC Computing Centre, Vancouver, B.C. Wal l e r , K. and M a k l l a , P. 1981. Ind. Eng. Chem. Process Des. Dev. 20:1-11 . n.1 Warga, J . 1963. J . Soc. Indust. A p p l . Math. n . : ^ 4 " 6 0 6 * 121 Westerberg, A.W. 1981. "O p t i m i z a t i o n i n Computer Aided Design", in "Foundations of Computer Aided Chemical Process Design", Mah, "R.S.H. and S e i d e r , W.D. E d i t o r s . Pub. by the Nat. S c i . Foundation, NY,NY. Pp.149- 183. White, W. Johnston, S. and D a n t z i g , G. 1958. J . Chem. Phys. 28: 751-755. Z e l e z n i k , F . J . and Gordon, S. 1968. Ind. Eng. Chem. 60:27-57. n. 6 1 22 NOMENCLATURE aj = A c t i v i t y of s p e c i e s j . A = Augmented exponents matrix. AA = Exponents matrix. A A i j atom grams of element i are present i n one mole of s p e c i e s j . Bi = Amount of atom grams of element i . Cj = Free energy c o e f f i c i e n t of sp e c i e s j . For gases, Cj = >j°(T)/RT + In P/1 atm. ; f o r condensed phases, Cj = /4.j»(T,P)/RT cj = Geometric program (GP) c o e f f i c i e n t s ; Cj = exp(-Cj) D = Number of independent chemical r e a c t i o n s . g 0 = Primal GP o b j e c t i v e f u n c t i o n . g K = Primal GP c o n s t r a i n t s . G = Gibbs f r e e energy. Gj = Molar Gibbs f r e e energy. h = O b j e c t i v e f u n c t i o n of the transformed p r i m a l GP. h = In g G = G/RT I = I d e n t i t y matrix. J = Hessian of In v. K = Number of phases at e q u i l i b r i u m . Kd = E q u i l i b r i u m constant f o r r e a c t i o n d. M = Number of elements. N = Number of chemical s p e c i e s at e q u i l i b r i u m . P = Pressure. rj = Extent of r e a c t i o n d. R = U n i v e r s a l gas con s t a n t . Sj = Vector of Gordon and McBride c o e f f i c i e n t s f o r sp e c i e s j t; = Primal v a r i a b l e a s o c i a t e d with element i T = Temperature. Ujd = S t o i c h i o m e t r i c c o e f f i c i e n t of sp e c i e s j f o r 123 r e a c t i o n d. v = Dual GP o b j e c t i v e f u n c t i o n . -In v = G/RT Xj = Molar f r a c t i o n of s p e c i e s j at e q u i l i b r i u m . y = V a r i a b l e of the l i n e a r program approximation to the dual GP. Zi = Ln-transformed p r i m a l v a r i a b l e . Z = l n t GREEK LETTERS: a = Parameter i n uni d i m e n s i o n a l search f o r GRG. 6j = Number of moles of s p e c i e s j . 6 j 0 = I n i t i a l number of moles of sp e c i e s j . %e = P e r c e n t u a l e r r o r of the composition determined by s e n s i t i v i t y a n a l y s i s , r e f e r r e d to r e - o p t i m i z a t i o n . \ K= T o t a l number of moles i n phase k. tij - Chemical p o t e n t i a l of s p e c i e s j in phase k. 0 = Reference chemical p o t e n t i a l . FlK = Lagrange m u l t i p l i e r f o r the k-primal c o n s t r a i n t . When s o l v i n g the transformed p r i m a l problem, n * = X-K e. = V a r i a b l e of the l i n e a r program approximation to the dual GP. I t accounts f o r n o n - l i n e a r i t i e s . APPENDIX A : COMPUTER PROGRAMS 125 PROGRAM JOTA C C C C THIS PROGRAM PERFORMS SENSITIVITY ANALYSIS BY INVERTING THE C JACOBIAN MATRIX C C IMPLICIT REAL*8 (A-H,0-Z) DIMENSION DA(20,20),DT(20,20),I PERM(40) DIMENSION U(20,20),DELT(20),C(20),C1(20),SUM(20,20) DIMENSION SUMA(20),SUMB(20),SUMC(20) C READ THE DATA READ (5,10) N,NDIMA,NDIMT 10 FORMAT (313) READ (5,20) N1,N2,OPTIM 20 FORMAT ( 213,F10.0) READ (5,30) (DELT(I),1=1,N2) 30 FORMAT (6F10.0) READ (5,40) (C(I),1=1,N2) 40 FORMAT (6F10.0) READ (5,50) (C1(I),1=1,N2) 50 FORMAT (6F10.0) READ (5,60) ((U(I,J),J=1,N2),1=1,N) 60 FORMAT (12F4.0) C FORM THE J MATRIX DO 9 I=1,N DO 9 J=1,N SUM(I,J)=0.D0 DO 8 K=1,N1 8 SUM(I,J)=SUM(I,J)+U(I,K)*U(J,K)/DELT(K) WRITE (6,200) I,J,SUM(I,J) 200 FORMAT (' I=',I4,' J=',I4,' SUM(I,J)=',G20.12) DA(I,J)=SUM(I,J)-U(I,N2)*U(J,N2)/DELT(N2) 9 CONTINUE C WRITE THE J MATRIX WRITE (6,70) 7 0 FORMAT (' MATRIX J') WRITE (6,80)((DA(I,j),J=1,N),1=1,N) 80 FORMAT (1X,6G14.6) C CALCULATE THE TIME TIME=SCLOCK(0.0) C CALCULATE THE INVERSE CALL INV(N,NDIMA,DA,I PERM,NDIMT,DT,DDET,JEXP,DCOND) C WRITE THE EXECUTION TIME WRITE (6,900) TIME 900 FORMAT (' EXECUTION TIME IS',F10.5) IF (DDET) 1,2,1 C WRITE THE INVERSE 1 WRITE (6,90) DCOND 90 FORMAT (' COND NO.=',G20.12,' INVERSE') WRITE (6,100)((DT(I,J),J=1,N),1=1,N) 100 FORMAT (1X,6G14.6) C CALCULATE THE NEW OPTIMUM SUMD=0.DO DO 3 1=1,N1 3 SUMD=SUMD+DELT(I)*(C1(I)-C(l)) OPTIM=OPTIM+SUMD C WRITE THE NEW OPTIMUM WRITE (6,120) OPTIM 120 FORMAT (' NEW OPTIM.=',6G14.6) C CALCULATE THE NEW SOLUTION DO 12 1=1,N2 SUMC(I)=0.D0 DO 11 L=1,N SUMA(L)=0.D0 DO 21 K=1,N SUMB(K)=0.D0 DO 22 J=1,N1 22 SUMB(K)=SUMB(K)+U(K,J)*(C1(J)-C(J)) 21 SUMA(L)=SUMA(L)+DT(L,K)* SUMB(K) 11 SUMC(I)=SUMC(I)+SUMA(L)*U(L,I) 12 DELT(I)=DELT(I)+SUMC(I) C WRITE THE NEW SOLUTION WRITE (6,130) 130 FORMAT (' NEW SOLUTION') WRITE (6,140) (DELT(I),1=1,N2) 140 FORMAT (1X,6G14.6) STOP 2 WRITE (6,110) 110 FORMAT (' INVERSION FAILED') STOP END 1 27 Program NPLUSK C C C C THIS PROGRAM PERFORMS SENSITIVITY ANALYSIS BY SOLVING A C SYSTEM OF LINEAR EQUATIONS C c IMPLICIT REAL*8 (A-H,0~Z) DIMENSION DA(12,12), DT(12,12), DB(12), DX(12), 1 DRZ(12), IPERM(24), A(12,12), DELT(12) DIMENSION C(12), C I ( 1 2 ) r U(12,12), SUM(12),V(12) READ (5,10) N1,K,M 10 FORMAT (313) READ (5,20) N,NDIMAT,N2 20 FORMAT (313) DEPS=1.D-14 C READ IN DATA TO FORM THE (N1+KXN1+K) MATRIX READ (5,30) ((A(I,J),J=1,N),1=1,M) 30 FORMAT (12F4.0) READ (5,40) (DELT(I),1=1 , N ) 40 FORMAT (6F10.0) READ (5,50) (C(I),1=1,N) 50 FORMAT (6F10.0) READ (5,60) (C1(I),I=1,N) 60 FORMAT (6F10.0) READ (5,70) ((U(I,J),J=1,N),1=1,N2) 70 FORMAT (12F4.0) NRHS=1 ITMAX=14 C FORM THE FINAL MATRIX C FIRST THE NORMALITY CONDITION DA(1,1)=1.DO DO 1 J=2,N 1 DA(1,J)=0.D0 C NOW THE ORTHOGONALITY COND. DO 2 1=1,M DO 2 J=1 ,N L= 1 + 1 2 DA(L,J)=A(I,J) C NOW,SUMMATION OF NUMBER OF MOLES DO 4 1=1,2 L=I+M+1 DA(L,1)=0.D0 DO 5 J=2,N1 5 DA(L,J)=1.D0 6 DA(L,N)=-1.DO 4 CONTINUE C NOW, THE EQUILIBRIUM CONDITIONS DO 7 1=1,N2 DO 7 J=1,N L=I+M+1+K 7 DA(L,J)=U(I,J)/DELT(J) C NOW CALCULATE THE RIGHT HAND SIDE VECTOR 1 2 8 M3=1+M+K DO 8 1 = 1 R M 3 8 D B ( I ) = 0 . D 0 DO 11 1 = 1 , N 2 SUM(I)=0.D0 DO 9 J = 1 , N 1 9 S U M ( I ) = S U M ( I ) + U ( l , J ) * ( C 1 ( J ) - C ( J ) ) L=I+M+1+K 11 D B ( L ) = S U M ( I ) C W R I T E T H E DA M A T R I X AND T H E DB V E C T O R W R I T E ( 6 , 8 5 0 ) 8 5 0 FORMAT (' M A T R I X OF C O E F F I C I E N T S * ) W R I T E ( 6 , 9 0 ) ( ( D A ( I , J ) , J = 1 , N ) , 1 = 1 , N ) 90 FORMAT ( 1 X , 6 G 1 2 . 6 ) W R I T E ( 6 , 9 5 0 ) 9 5 0 FORMAT (' B V E C T O R ' ) W R I T E ( 6 , 1 0 0 ) ( D B ( I ) , 1 = 1 , N ) 100 FORMAT ( 1 X , 6 G 1 2 . 6 ) C C A L C U L A T E T H E T I M E OF E X E C U T I O N T I M E = S C L O C K ( 0 . 0 ) C S O L V E T H E S Y S T E M C A L L D S L I M P ( D A , D T , D B , D X , D R Z , I P E R M , N , N D I M A T , D E P S , N R H S , I T M A X ) T I M E = S C L O C K ( T I M E ) C W R I T E T H E E X E C U T I O N T I M E W R I T E ( 6 , 8 0 0 ) T I M E 8 0 0 FORMAT (' E X E C U T I O N T I M E I S ' , F 1 0 . 5 ) C W R I T E OUT R E S U L T S W R I T E (6,goo-goo FORMAT (' V A R I A T I O N OF NUMBER OF M O L E S ' ) W R I T E ( 6 , 1 1 0 ) ( D X ( I ) , 1 = 1 , N ) 1 1 0 FORMAT ( 1 X , 6 G 1 4 . 6 ) C C A L C U L A T E T H E NEW NUMBER OF M O L E S DO 15 I = 1 , N 15 V ( I ) = D X ( I ) + D E L T ( l ) W R I T E ( 6 , 1 2 0 ) 1 2 0 FORMAT (' S O L U T I O N ' ) W R I T E ( 6 , 1 3 0 ) ( V ( I ) , 1 = 1 , N ) 13 0 FORMAT ( 1 X , 6 G 1 4 . 6 ) S T O P E N D 1 29 Program C0MP1 C C C THIS PROGRAM SOLVES THE PRIMAL GP FOR THE CHEMICAL EQUILIBRIUM C PROBLEM. IT CALLS SUBROUTINE LIPSU2 TO GET A FIRST DUAL POINT C FOR THE OPTIMIZATION. THEN IT CALLS SUBROUTINE SINGV WHICH C CALCULATES THE FIRST PRIMAL STARTING POINT. SUBROUTINES GRGIN, C GRG2 AND GRGEG, FROM UBC,SOLVE THE OPTIMIZATION PROBLEM. C SUBROUTINE GCOMP CALCULATES THE OBJECTIVE FUNCTION AND THE C CONSTRAINTS, AND IS CALLED FROM GRG. THE PROCEDURE IS C REPEATED FOR DIFFERENT TEMPERATURES AND PRESSURES. C SUBROUTINE FREEN CALCULATES THE FREE ENERGY COEFFICIENTS C (G/RT), USING THE DATA FROM GORDON AND MC. BRIDE C C C IMPLICIT REAL*8 (A-H fO-Z) LOGICAL*1 BLANK /' '/, SEMIC /';'/ -C DIMENSION STATMENTS . LABELLED COMMONS ARE USED. C N=NUMBER OF CHEMICAL SPECIES +1 C M=NUMBER OF CHEMICAL ELEMENTS (ATOMS) + 1 C K=NUMBER OF PHASES .. C NN= NUMBER OF DIFFERENT TEMP. AND PRESSURE CONDITIONS. C NE=1 IF THE FREE ENERGY COEFFICIENTS ARE ALREADY AVAILABLE. C IF THEY ARE TO BE CALCULATED BY THE GORDON AND MC.BRIDE C COEFFICIENTS, USE ANY OTHER INTEGER. C NF=1 IF THE DUAL STARTING POINT ROUTINE IS TO BE USED. C OTHERWISE, EQUAL NF TO ANY OTHER INTEGER, AND PROVIDE C A DUAL STARTING POINT. DON'T FORGETTHAT C THE FIRST DUAL VARIABLE IS DUMMY AND EQUAL TO 1.D0. C A(I,J),I=1,M,J=1,N = EXPONENTS MATRIX. A(1,J) CORRESPONDS C TO THE NORMALITY CONDITION C C(JJ,L),JJ=1,N, L=1,NN FREE ENERGY COEFFICIENTS OF SPECIES J J C T,P=EQUILIBRIUM TEMPERATURE AND PRESSURE C X(I)=Z(I,1),I=1,M-1 TRANSFORMED PRIMAL VARIABLES.FOR C SUBROUTINE LIPSU2,X(J),J=1,N+1, CORRESPOND TO THE VARIABLES C LINEAR PROGRAM. C XNEW(J),J=2,N =NUMBER OF MOLES OF SPECIES J ; XNEW(1)=1.D0 C IS A DUMMY VARIABLE TO ACCOUNT FOR THE NORMALITY CONDITION. C XMF(J)=M0LAR FRACTION OF SPECIES J,.J=2, 3, N. C XLAG(I),I=1,K = LAGRANGE MULTIPLIER OF CONSTRAINT I C S(J,L),J=1,N-1,L=1,7 COEFFICIENTS FOR GORDON & MC.BRIDE C POLYNOMIALS TO CALCULATE THE FREE ENERGY COEFFICIENTS. C B(I),I=2,M = VECTOR OF AMOUNTS OF ELEMENT I.B(1)=1, DUMMY. C F = OBJECTIVE FUNCTION. COMMON/AX/A(10,30) COMMON/BX/B(10) COMMON/CX/C(30, 10) COMMON/DX/N,M/HX/K,NN/GX/F/XL/L COMMON/EX/T(30),P(30) COMMON/FX/X(31)/ZX/Z(10,5) DIMENSION XMF(30),VL(30),VH(30) COMMON/XNE/XNEW( 3 1 ) /SX/'S (30,10) DIMENSION XLAG(10),SUMC(30) 1 30 C C READ IN DATA C READ(5,10) N,M,K,NN,NE,NF 10 FORMAT (615) C C CALCULATE FREE ENERGY COEFFICIENTS BY GORDON AND MC.BRIDE C IF (NE.EQ.1) GO TO 13 SX =550.DO C C CALCULATE THE TEMPERATURE AT FIXED INTERVALS. C DO 20 L=1,NN SX=SX +5.D1 T(L)=SX 20 P(L)=1.D0 N1=N-1 READ (5,40) ((S(I,L),L=1,7),I=1,N1) 40 FORMAT (5G13.7) CALL FREEN GO TO 400 13 READ (5,311) (T(L),L=1,NN) C C IF THE C(I,L) MATRIX IS AVAILABLE BY OTHER SOURCES THAN C GORDON & MC.BRIDE, READ NOW TEMPERATURE, PRESSURE, AND THE C C(I,L) MATRIX. C READ (5,311) (P(L),L=1,NN) READ (5,311) ((C(I,L),L=1,NN),I=1,N) 311 FORMAT (5F10.0) C C READ THE EXPONENT MATRIX AND THE B VECTOR. C 400 READ (5,50) ((A(I,J),J=1,N),I=1,M) 50 FORMAT (9F4.0) READ (5,60) (B(I),1=1,M) 60 FORMAT (5F10.0) C C PUT NF=1 IF YOU WANT THE PROGRAM TO CALCULATE A DUAL STARTING C POINT. OTHERWISE, IT WILL NOW READ YOUR FIRST DUAL GUESS. C IF (NF.EQ.1) GO TO 500 READ (5,511) (XNEW(J),J=1,N) 511 FORMAT (5F10.0) TIME=SCLOCK(0.) GO TO 600 500 TIME=SCLOCK(0.) C C CALLS SUBROUTINE LIPSUB2.IT PROVIDES A STARTING POINT FOR THE C DUAL PROBLEM. CALL LIPSU2 XNEW(1)=1.D0 DO 100 J=2,N 131 100 XNEW(J) = X(J - 1) + X(N) WRITE (6,12) (XNEW(J),J=1,N) 12 FORMAT (' FIRST DUAL POINT 1/1X,6G16.8) C C SUBROUTINE SINGV WILL NOW TRANSFORM THE DUAL VARIABLES INTO C PRIMAL ONES. C 600 CALL SINGV M1=M-1 WRITE (6,16) (Z(I,1) ,1 = 1 ,M1) 16 FORMAT (' FIRST PRIMAL STARTING POINT' ,4G16.8 ) DO 11 1=1,M1 11 X(I)=Z(I,1) DO 202 L=1,NN CALL FTNCMD('ASSIGN 5=-DATA;') CALL FTNCMD ('ASSIGN 7=*SINK*;') C OPTIMIZATION IS PERFORMED FOR NN DIFFERENT P AND T CONDITIONS C A SCRATCH FILE IS CREATED TO WRITE DOWN.THE DATA NEEDED FOR C THE UBC SUBROUTINES GRGIN AND GRG TO PERFORM THE OPTIMIZATION C C FIRST ARE THE CONTROL CARDS C IG=0.D0 WRITE (5,15) Ml,K,IG 15 FORMAT (316) C C NOW WRITE THE LOWER BOUNDS OF THE VARIABLES.WE ADD -30 TO THE C PRIMAL STARTING POINT. C WRITE (5,25) 25 FORMAT ('LBV=') DO 115 1=1,M1 115 VL(I)=X(I)-3.D1 WRITE (5,35) (BLANK,I,VL(I),I=1,M1),SEMIC 35 FORMAT (6(A1,I3,G10.2)) C NOW WRITE THE UPPER BOUNDS OF THE VARIABLES.WE SET THEM AS 0.0 C WRITE (5,55) 55 FORMAT ('UBV=') DO 12 5 I=1,M1 125 VH(I)=0.D0 WRITE (5,35) (BLANK,I,VH(I),I=1,M1),SEMIC C C PRIMAL CONSTRAINTS. IF THERE ARE MORE THAN ONE, THE FORMAT C SHOULD BE :FORMAT('UBC= 1 I.DO 2 1.D0 K 1.DO ;') C WRITE (5,66) 66 FORMAT('UBC= 1 1.DO ;') WRITE (5,65) 65 FORMAT ('QUAD') C C EPNEWT IS THE TOLERANCE FOR THE PRIMAL CONSTRAINTS.THE ACCURACY C OF THE RESULTS IS VERY SENSITIVE TO THIS VALUE. C 1 32 WRITE (5,75) 75 FORMAT ('EPNEWT=1.D~10') WRITE (5,85) C C EPSTOP IS THE TOLERANCE FOR THE OBJECTIVE FUNCTION C 85 FORMAT ('EPSTOP=1.D-6') WRITE (5,95) 95 FORMAT ('EPSBOUND=1.D-6') WRITE (5,105) 105 FORMAT ('EPSPIV=1.D-6') WRITE (5,61) C C PRINTCTL MAY BE SET AT HIGHER VALUES (UP TO 4) FOR MORE C INFORMATION ON THE OPTIMIZATION PROGRESS. C 61 FORMAT (' PRINTCTL=1') WRITE (5,107) 107 FORMAT('X=') WRITE (5,109) (X(I) ,1=1,M1) 109 FORMAT (6G18.6) WRITE (5,108) 108 FORMAT ('OPTIMIZE*/'GO'/'STOP * ) REWIND 5 C NOW GRGIN CAN READ INPUT FROM -DATA;GRG2 AND GRGRES WILL C THEIR OUTPUT INTO THE SCRATCH FILE -GRGOUT CALL FTNCMD('ASSIGN 6=-GRGOUT;') CALL GRGINU1 ,&2) CALL GRG2(X,F,&1) CALL GRGRES DO 22 1=1,M1 22 Z(I,L)=X(I) GO TO 2 1 WRITE (7,3) 3 FORMAT (' GRG HAS FAILED') GO TO 202 2 REWIND 6 C C SUBROUTINE SKIP READS THE LAGRANGE MULTIPLIERS FOR EACH PRIMAL C CONSTRAINT. THEY ARE THE NEGATIVE OF THE TOTAL NUMBER OF MOLES C IN EACH CONSTRAINT. C CALL SKIPU202) READ (6,26) (XLAG(I),I=1,K) 26 FORMAT (1X,1 OF 13.0) WRITE (7,205) T ( L ) , P ( L ) , F 205 FORMAT ('T(K)=',F10.0,'P(AT)=',F10.0,'OBJ.F=',G16.8) WRITE (7,215) 215 FORMAT (' MOLAR FRACTIONS') DO 33 J=2,N SUMC(J)=0.D0 DO 31 1=1,M1 11=1+1 31 SUMC(J).= SUMC(J)+A(II , J ) * X ( I ) IF (SUMC(J).GE.0.D0) GO TO 34 R=DEXP(SUMC(J)-C(J,L)) GO TO 33 34 R=0.D0 33 XMF(J)=R WRITE (7 , 235) (XMF(J),J = 2,N) 235 FORMAT (1X,6G16.8) WRITE (7,238) (XLAG(I),I=1,K) 238 FORMAT (' TOTAL NUMBER OF MOLES',G16.8) WRITE (7,239) 239 FORMAT ('NUMBER OF MOLES') DO 237 J=2,N 237 . XNEW(J)=-XLAG(1)*XMF(J) WRITE (7,240) (XNEW(J),J=2,N) 240 FORMAT (1X,6G16.8) 202 CONTINUE 23 TIME=SCLOCK(TIME) WRITE (7,19) TIME 19 FORMAT (' EXECUTION TIME =',F6.4) STOP END C C SUBROUTINE LIPSU2 CALCULATES A DUAL STARTING POINT C BY APPROXIMATING THE DUAL PROBLEM TO A LINEAR PROGRAM. C THE LINEAR PROGRAM ISTHEN SOLVED BY SUBROUTINE C LIPSUB FROM UBC. C SUBROUTINE LIPSU2 IMPLICIT REAL*8 (A-H,0~Z) DIMENSION TABLO(30,30),NVIN(30),NVOUT(30),BBOBJ(20), 1 UBOBJ(20),BRHS(20),UBRHS(10) COMMON/AX/A(10,30) COMMON/BX/B(10) COMMON/CX/C(30,10)/GX/OPTIM COMMON/DX/NVARS,NCONST COMMON/FX/X(31) N1=NVARS-1 M1=NCONST-1 NEQUAL=M1 MAXIM=0 IFOBJ=0 IFRHS=0 NCOLS=NVARS+1 NROWS=NCONST+1 C ZERO THE TABLEAU DO 200 J=1,NCOLS DO 2 00 1=1,NROWS 200 TABLO(l,J)=0.D0 C FIRST ROW IN THE TABLEAU IS THE OBJECTIVE FUNCTION SUM=0.D0 DO 101 J=1,N1 JJ=J+1 TABLO(1,J)=C(JJ,1) 101 SUM=SUM+TABLO(1,J) TABLO(1,NVARS)=SUM TABLO(2,NVARS)=-1.DO TABLO(2,NCOLS)=-.0001 DO DO 203 I=2,NCONST DO 202 J=1,N1 202 TABLO(I + 1, J) = A(I, J + 1) TABLO(I + 1, NVARS) = 1.DO 203 TABLO(I + 1, NCOLS) = B ( l ) 2 60 NLVARS = NVARS CALL LIPSUB(TABLO, 30, NCONST, NLVARS, NEQUAL, MAXIM, 1 IFOBJ, IFRHS, TOL, NVIN, NVOUT,BBOBJ, 2 UBOBJ, BBRHS, UBRHS, &540) NP1=NLVARS+1 OPTIM =TABLO(1,NP1) IF (MAXIM.NE.1) OPTIM=-OPTIM DO 12 1=1,NVARS 12 X(I)=0.D0 DO 13 I=2,NROWS J = NVIN(I) 13 X(J)=TABL0(I,NP1) WRITE (6,320) (X(J),J=1,NVARS) 320 FORMAT (1X,5G18.6) RETURN 540 WRITE (6,560) 560 FORMAT ('0DUAL STARTING POINT ROUTINE FAILED') DO 580 1 = 1 , NVARS 580 X(I) = 10.0D0 RETURN END C C SUBROUTINE GCOMP CALCULATES THE VALUES OF THE OBJECTIVE C FUNCTION AND OF THE CONSTRAINTS C SUBROUTINE GCOMP(G,X) IMPLICIT REAL*8 (A~H,0"Z) COMMON/AX/A(10,3 0)/BX/B(10)/CX/C(30,10)/DX/N,M/XL/L COMMON/HX/K,NN DIMENSION G(1),X(1),SUMF(5),XF(30),XMF(30) DO 10 I=1,K SUMF(I)=0.DO DO 20 J=2,N XF(J)=0.D0 DO 30 11 = 2,M 11=11-1 30 XF(J)=XF(J)+A(II,J)*X(I1) IF (XF(J).GE.0.D0) GO TO 40 XMF(J)=XF(J)-C(J,L) R=DEXP(XMF(J)) GO TO 20 40 R=0.D0 20 SUMF(I)=SUMF(I)+R 10 G(I)=SUMF(I) SUM=0.D0 DO 50 J=2,M JJ=J-1 50 SUM=SUM+B(J)*X(JJ) G(K+1) = -SUM RETURN END C C SUBROUTINE SINGV CALCULATES A FIRST STARTING POINT C FOR THE PRIMAL PROBLEM.IT SOLVES AN OVERDETERMINED C SYSTEM OF EQUATIONS 'THAT RELATES PRIMAL AND DUAL VARIABLES. C THE METHOD USED IS A SINGULAR VALUE DECOMPOSITION. C SUBROUTINE SINGV IMPLICIT REAL*8(A-H fO-Z) DIMENSION A(30,30),V(30,30),S(30) COMMON/AX/AD(10,30)/CX/C(30,10)/XNE/XNEW(31)/DX/N1,M1 COMMON/ZX/X(10,5)/HX/K,NN NP=1.D0 N=N1-1 M=M1-1 MNP=M+NP NDIMAU=30 NDIMV=30 DO 10 I=1,N 11=1+1 DO 10 J=1,M JJ=J+1 10 A(I,J)=AD(JJ,II) SUMX=0.D0 DO 12 J=2,N1 12 SUMX=SUMX+XNEW(J) DO 11 I=1,N 11=1+1 DO 11 J=M1,MNP . IF (XNEW(II).LE.O.DO) XNEW(II)=1.D-16 11 A(I,J)=C(11,1 )+DLOG(XNEW(II))-DLOG(SUMX) WRITE (6,14) ((A(I,J),J=1,MNP),1=1,N) 14 FORMAT (1X,5G18.6) CALL DSLSVD(A,S,V,NDIMAU,NDIMV,N,M,NP,&140) EPS=1.D-6 SS=S(1)*EPS DO 60 J=1,M IF (S(J).LT.SS) GO TO 70 DO 50 1=1,M 50 V (I , J) = V (I , J) /S (J) 60 CONTINUE J=M+1 GO TO 90 70 WRITE (6,80) ((V(l,K),1=1,M),K=J,M) 80 FORMAT (1X,4G13.5) 90 IF (J.GT.M) GO TO 120 DO 110 K=J,M DO 110 1=1,M 110 V(I,K)=0.D0 120 MP1=M+1 136 CALL DGMULT(V,A(1,MP1),X,M,M,NP,NDIMV,NDIMAU,10) RETURN 140 WRITE (6,150) 150 FORMAT (' ERROR RETURN FROM DSLSVD') STOP END C C SUBROUTINE SKIP WILL READ THE LAGRANGIAN MULTIPLIERS C FOR THE PRIMAL CONSTRAINTS AS PRINTED BY GRG IN THE SCRATCH C FILE -GRGOUT. C SUBROUTINE SKIP(*) LOGICAL* 1 RECORDO50) INTEGER*2 LEN LOGICAL EQCMP 100 CALL READ(RECORD,LEN,0,LNR,6,&200) IF (LEN .LT. 70 .OR. LEN .GT. 75) GO TO 100 IF (EQCMP (21, RECORD, '0LAGRANGE MULTIPLIERS')) RETURN GO TO 100 200 RETURN 1 END C C SUBROUTINE FREEN CALCULATES THE FREE ENERGY PARAMETERS C FROM THE GORDON AND MC.BRIDE POLYNOMIALS COEFFICIENTS. C SUBROUTINE FREEN IMPLICIT REAL*8 (A-H,0-Z) COMMON/CX/C(30,10)/SX/S(30,10) COMMON/DX/N,M/HX/K,NN COMMON/EX/T(30),P(30) DO 10 L=1,NN C(1,L)=0.DO N1=N-1 DO 10 J=1,N1 JJ=J+1 10 C(JJ,L)=S(J,1)*(1.D0-DLOG(T(L)))~S(J,2)*T(L)/2.D0-1 (S(J,3)*T(L)**2.D0)/6.D0 -(S(J,4)*T(L)**3.DO)/1.2D1 2 -(S(J,5)*T(L)**4.D0)/2.D1+S(J,6)/T(L)-S(J,7) 3 +DLOG(P ( D ) WRITE(6,11) ((C(J,L),J=1,N),L=1,NN) 11 FORMAT ('MATRIX OF FREE ENERGY COEFICIENTS',/,1X,4F12.7) RETURN END 137 Program C0MP2 C C C C THIS PROGRAM SOLVES THE PRIMAL GP FOR THE CHEMICAL EQUILIBRIUM C PROBLEM. IT CALLS THE SUBROUTINE LIPSU2 TO GET A FIRST DUAL C POINT FOR THE OPTIMIZATION. THEN IT CALLS SUBROUTINE SINGV E C WHICH CALCULATES THE CORRESPONDING PRIMAL C STARTING POINT. SUBROUTINES GRGIN AND GRG2, FROM UBC, C SOLVE THE OPTIMIZATION PROBLEM.SUBROUTINE GCOMP CALCULATES THE C OBJECTIVE FUNCTION AND THE CONSTRAINTS, AND IS CALLED FROM GRG. C SUBROUTINE FREEN CALCULATES THE FREE ENERGY COEFFICIENTS (G/RT) C USING THE DATA FROM GORDON AND MC.BRIDE. C A SENSITIVITY ANALYSIS OF THE DUAL VARIABLES IS PERFORMED BY C SUBROUTINE NPLUSK C C IMPLICIT REAL*8 (A-H,0-Z) LOGICAL*1 BLANK /' '/, SEMIC /';'/ C DIMENSION STATMENTS . LABELLED COMMONS ARE USED. C N=NUMBER OF CHEMICAL SPECIES +1 C M=NUMBER OF CHEMICAL ELEMENTS (ATOMS) + 1 C K=NUMBER OF PHASES C NN= NUMBER OF DIFFERENT TEMP. AND PRESSURE CONDITIONS. C NE=1 IF THE FREE ENERGY COEFFICIENTS ARE ALREADY AVAILABLE. IF C THEY ARE TO BE CALCULATED BY THE GORDON AND MC.BRIDE C COEFFICIENTS, USE ANY OTHER INTEGER . C NF=1 IF THE DUAL STARTING POINT ROUTINE IS TO BE USED. C OTHERWISE, EQUAL NF TO ANY OTHER INTEGER, C AND PROVIDE A DUAL STARTING POINT. DON'T FORGETTHAT C THE FIRST DUAL VARIABLE IS DUMMY AND EQUAL TO 1.D0. C A(l,J),I=1,M,J=1,N = EXPONENTS MATRIX. A(J,1) CORRESPONDS C TO THE NORMALITY CONDITION. C C(JJ,L),JJ=1,N, L=1,NN FREE ENERGY COEFFICIENTS OF SPECIES J J C T,P=EQUILIBRIUM TEMPERATURE AND PRESSURE C X(I)=Z(I,1),I=1,M-1 TRANSFORMED PRIMAL VARIABLES. FOR C SUBROUTINE LIPSU2,X(J),J=1,N+1, CORRESPOND TO THE VARIABLES OF C THE LINEAR PROGRAM. C XNEW(J),J=2,N =NUMBER OF MOLES OF SPECIES J ; XNEW(1)=1.D0 IS A C DUMMY VARIABLE CORRESPONDING TO THE NORMALITY CONDITION. C XMF(J)=M0LAR FRACTION OF SPECIES J C XLAG(I),1=1,K = LAGRANGE MULTIPLIER OF CONSTRAINT I C S(J,L),J=1,N-1,L=1,7 COEFFICIENTS FOR GORDON & MC.BRIDE C POLYNOMIALS C U(I,J) 1=1,N-M, J=1,N IS THE MATRIX OF STOICHIOMETRIC COEFFICIE C B(I),I=2,M = VECTOR OF AMOUNTS OF ELEMENT I.B(l)=1, DUMMY. C F = OBJECTIVE FUNCTION. COMMON/AX/A(10,30) COMMON/BX/B(10) COMMON/CX/C(30, 1 0) COMMON/DX/N,M/HX/K,NN/GX/F/XL/L COMMON/EX/T(30),P(30) COMMON/FX/X(31)/ZX/Z(10,5) DIMENSION XMF(30),VL(30),VH(30) COMMON/XNE/XNEW(31)/SX/SUMX/SUX/S(30,10)/UX/U(30,30) DIMENSION XLAG(10),SUMC(30) C C READ IN DATA C READ(5,10) N,M,K,NN,NE,NF 10 FORMAT (615) C C CALCULATE FREE ENERGY COEFFICIENTS BY GORDON AND MC.BRIDE C IF (NE.EQ.1) GO TO 13 SX =550.DO C C CALCULATE THE TEMPERATURE AT FIXED INTERVALS. C DO 20 L=1,NN SX=SX +5.D1 T(L)=SX 20 P(L)=1.D0 N1=N-1 READ (5,40) ((S(I,L),L=1,7),I=1,N1) 4 0 FORMAT (5G13.7) CALL FREEN GO TO 400 13 READ (5,311) (T(L),L=1,NN) C C IF THE C(I,L) MATRIX IS AVAILABLE BY OTHER SOURCES THAN C GORDON & MC.BRIDE, READ NOW TEMPERATURE, PRESSURE, AND THE C C(I,L) MATRIX. C READ (5,311) (P(L),L=1,NN) READ (5,311) ((C(I,L),L=1,NN),I=1,N) 311 FORMAT (5F10.0) C C READ THE EXPONENT MATRIX AND THE B VECTOR. C 400 READ (5,50) ((A(I,J),J=1,N),I=1,M) 50 FORMAT (13F4.0) READ (5,60) (B(I),1=1,M) 60 FORMAT (5F10.0) C C READ THE MATRIX OF STOICHIOMETRIC COEFFICIENTS. C N2=N-M READ (5,70) ((U(I,J),J=1,N),1=1,N2) 7 0 FORMAT (10F4.0) C C PUT NF=1 IF YOU WANT THE PROGRAM TO CALCULATE A DUAL STARTING C POINT. OTHERWISE, IT WILL NOW READ YOUR FIRST DUAL GUESS. C IF (NF.EQ.1) GO TO 500 READ (5,511) (XNEW(J),J=1,N) 511 FORMAT (5F10.0) TIME=SCLOCK(0.) 1 39 GO TO 600 500 TIME=SCLOCK(0.) C C CALLS SUBROUTINE LIPSUB2.IT PROVIDES A STARTING POINT FOR THE C DUAL PROBLEM. CALL LIPSU2 XNEW(1)=1.DO DO 100 J = 2,N 100 XNEW(J) = X(J - 1) + X(N) WRITE (6,12) (XNEW(J),J=1,N) 12 FORMAT (' FIRST DUAL POINT'/1X,6G16.8) C C SUBROUTINE SINGV WILL NOW TRANSFORM THE DUAL VARIABLES INTO C PRIMAL ONES. C 600 CALL SINGV M1=M-1 WRITE (6,16) (Z(I,1),1=1,M1) 16 FORMAT (' FIRST PRIMAL STARTING POINT',4G16.8) DO 11 1=1,M1 11 X(I)=Z(I,1) CALL FTNCMD('ASSIGN 5=-DATA;') CALL FTNCMD ('ASSIGN 7=*SINK*;*) C OPTIMIZATION IS PERFORMED FOR NN DIFFERENT P AND T CONDITIONS C A SCRATCH FILE IS CREATED TO WRITE DOWN THE DATA NEEDED FOR THE C UBC SUBROUTINES GRGIN AND GRG2 TO PERFORM THE OPTIMIZATION C C FIRST ARE THE CONTROL CARDS C IG=0.D0 WRITE (5,15) M1,K,IG 15 FORMAT (316) C C NOW WRITE THE LOWER BOUNDS OF THE VARIABLES.WE ADD -30 TO THE C PRIMAL STARTING POINT. C WRITE (5,25) 25 FORMAT ('LBV=') DO 115 1=1,M1 115 VL(I)=X(I)-3.D1 WRITE (5,35) (BLANK,I,VL(I),I=1,M1),SEMIC 35 FORMAT (6(A1,I3,G10.2)) C NOW WRITE THE UPPER BOUNDS OF THE VARIABLES.WE SET THEM AS 0.0 C WRITE (5,55) 55 FORMAT ('UBV=') DO 125 1=1,M1 125 VH(I)=0.D0 WRITE (5,35) (BLANK,I,VH(I),I=1,M1),SEMIC C C PRIMAL CONSTRAINTS. IF THERE ARE MORE THAN ONE, THE FORMAT C SHOULD BE : FORMAT('UBC= 1 1.D0 2 1.DO K 1.DO ;') C WRITE (5,66) 1 40 66 FORMAT('UBC= 1 1.DO ;') WRITE (5,65) 65 FORMAT ('QUAD') C C EPNEWT IS THE TOLERANCE FOR THE PRIMAL CONSTRAINTS.THE ACCURACY C OF THE RESULTS IS VERY SENSITIVE TO THIS VALUE. C WRITE (5,75) 75 FORMAT ('EPNEWT=1.D~10') WRITE (5,85) C C EPSTOP IS THE TOLERANCE FOR THE OBJECTIVE FUNCTION C 85 FORMAT ('EPSTOP=1.D-6') WRITE (5,95) 95 FORMAT ('EPSBOUND=1.D"6') WRITE (5,105) 105 FORMAT ('EPSPIV=1.D-6') WRITE (5,61) C C PRINTCTL MAY BE SET AT HIGHER VALUES (UP TO 4) FOR MORE C INFORMATION ON THE PROGRESS OF THE OPTIMIZATION. C 61 FORMAT (' PRINTCTL=1') WRITE (5,107) 107 FORMAT('X=') WRITE (5,109) (X(I),1=1,M1) 109 FORMAT (6G18.6) WRITE (5,108) 108 FORMAT ('OPTIMIZE'/'GO'/'STOP') REWIND 5 C NOW GRGIN CAN READ INPUT FROM -DATA;GRG2 AND GRGRES WILL WRITE C THEIR OUTPUTS INTO THE SCRATCH FILE -GRGOUT CALL FTNCMD('ASSIGN 6=~GRGOUT;') CALL GRGINU1 ,&2) CALL GRG2(X,F,&1) " CALL GRGRES GO TO 2 1 WRITE (7,3) 3 FORMAT (' GRG HAS FAILED') GO TO 202 2 REWIND 6 C C SUBROUTINE SKIP READS THE LAGRANGE MULTIPLIERS FOR EACH PRIMAL C CONSTRAINT. THEY ARE THE NEGATIVE OF THE TOTAL NUMBER OF C MOLES IN EACH CONSTRAINT. C CALL SKIP(&202) READ (6,26) (XLAG(I),1=1,K) 26 FORMAT (IX,1 OF 13.0) WRITE (7,205) T ( L ) , P ( L ) , F 205 FORMAT ('T(K)=',F10.0,'P(AT)=',F10.0,'OBJ.F=',G16.8) WRITE (7,215) 215 FORMAT (' MOLAR FRACTIONS') 141 DO 33 J=2,N SUMC(J)=0.DO-DO 31 1=1,M1 11=1+1 31 SUMC(J)=SUMC(J)+A(II,J)*X(I) IF (SUMC(J).GE.O.DO) GO TO 34 R=DEXP(SUMC(J)-C(J,L)) GO TO 33 34 R=0.D0 33 XMF(J)=R WRITE (7,235) (XMF(J),J = 2,N) 235 FORMAT (1X,6G16.8) WRITE (7,238) (-XLAG(I),I=1,K) 238 FORMAT (' TOTAL NUMBER OF MOLES 1 ,G16.8) WRITE (7,239) 239 FORMAT ('NUMBER OF MOLES') DO 237 J=2,N 237 XNEW(J)=-XLAG(1)*XMF(J) WRITE (7,240) (XNEW(J),J=2,N) 240 FORMAT (1X,6G16.8) XNEW(1)=1.D0 SUMX=-XLAG(1) CALL NPLUSK GO TO 23 23 TIME=SCLOCK(TIME) WRITE (7,19) TIME 19 FORMAT (' EXECUTION TIME =',F6.4) STOP END C C SUBROUTINE LIPSU2 CALCULATES A DUAL STARTING POINT APPROXIMATIN C THE DUAL PROBLEM TO A LINEAR PROGRAM.THE LINEAR PROGRAM IS C SOLVED BY SUBROUTINE LIPSUB FROM UBC. C SUBROUTINE LIPSU2 IMPLICIT REAL*8 (A-H,0-Z) DIMENSION TABLO(30,30),NVIN(30),NVOUT(30),BBOBJ(20), 1 UBOBJ(20),BRHS(20),UBRHS(10) COMMON/AX/A(10,30) COMMON/BX/B(10) COMMON/CX/C(30,10)/GX/OPTIM COMMON/DX/NVARS,NCONST COMMON/FX/X(31) N1=NVARS-1 M1=NCONST-1 NEQUAL=M1 MAXIM=0 IFOBJ=0 IFRHS=0 NCOLS=NVARS+1 NROWS=NCONST+1 C ZERO THE TABLEAU DO 200 J=1,NCOLS DO 200 1=1,NROWS 1 42 200 TABLO(I,J)= 0.DO C FIRST ROW IN THE TABLEAU IS THE OBJECTIVE FUNCTION SUM=0.D0 DO 101 J=1,N1 JJ=J+1 TABLO(1,J)=C(JJ, 1 ) 101 SUM=SUM+TABLO(1 , J) TABLO(1,NVARS)=SUM TABLO(2,NVARS)=-1.DO TABLO(2,NCOLS)=-.0001 DO DO 203 I=2,NCONST DO 202 J=1,N1 202 TABLO(I + 1, J) = A ( I , J + 1 ) TABLO(I + 1, NVARS) = 1.D0 203 TABLO(I + 1, NCOLS) = B ( l ) 260 NLVARS = NVARS CALL LIPSUB(TABLO, 30, NCONST, NLVARS, NEQUAL, MAXIM, 1 IFOBJ, IFRHS, TOL, NVIN, NVOUT,BBOBJ, 2 UBOBJ, BBRHS, UBRHS, &540) NP1=NLVARS+1 OPTIM =TABLO(1,NP1) IF (MAXIM.NE.1) OPTIM=-OPTIM DO 12 1=1,NVARS 12 X(I)=0.D0 DO 13 I=2,NROWS J = NVIN(I) 13 X(J)=TABLO(l,NP1 ) WRITE (6,320) (X(J),J=1,NVARS) 320 FORMAT (1X,5G18.6) RETURN 540 WRITE (6,560) 560 FORMAT ('0DUAL STARTING POINT ROUTINE FAILED') DO 580 I = 1, NVARS 580 X(I) = 10.0D0 RETURN END C C SUBROUTINE GCOMP CALCULATES THE VALUES OF THE OBJECTIVE FUNCTIO C AND OF THE CONSTRAINTS C SUBROUTINE GCOMP(G,X) IMPLICIT REAL*8 (A-H,0~Z) COMMON/AX/A(10,30)/BX/B(10)/CX/C(30,10)/DX/N,M/XL/L COMMON/HX/K,NN DIMENSION G(1),X(1),SUMF(5),XF(30),XMF(30) DO 10 I=1,K SUMF(I)=0.D0 DO 20 J=2,N XF(J)=0.D0 DO 30 11=2,M 11=11-1 30 XF(J)=XF(J)+A(II,J)*X(I1) IF (XF(J).GE.0.D0) GO TO 40 XMF(J)=XF(J)-C(J,L) 143 R=DEXP(XMF(J)) GO TO 2 0 40 R=0.D0 20 SUMF(I)=SUMF(I)+R 10 G(I)=SUMF(I) SUM=0.DO DO 50 J=2,M JJ=J-1 50 SUM=SUM+B(J)*X(JJ) G(K+1) = -SUM RETURN END C C SUBROUTINE SINGV CALCULATES A FIRST STARTING POINT FOR THE PRIM C PROBLEM.IT SOLVES AN OVERDETERMINED SYSTEM OF EQUATIONS THAT C RELATES PRIMAL AND DUAL VARIABLES. C THE METHOD USED IS A SINGULAR VALUE DECOMPOSITION. C SUBROUTINE SINGV IMPLICIT REAL*8(A-H,0-Z) DIMENSION A(30,30),V(30,30),S(30) COMMON/AX/AD(10,3 0)/CX/C(30,10)/XNE/XNEW(31)/DX/N1 ,M1 COMMON/ZX/X(10,5)/HX/K,NN NP=1.DO N=N1-1 M=M1- 1 MNP=M+NP NDIMAU=30 NDIMV=30 DO 10 I = 1 , N 11=1+1 DO 10 J=1,M JJ=J+1 10 A(I,J)=AD(JJ,II) SUMX=0.D0 DO 12 J=2,N1 12 SUMX=SUMX+XNEW(J) DO 11 I=1,N 11=1+1 DO 11 J=M1,MNP IF (XNEW(ll).LE.0.D0) XNEW(11) = 1 .D-16 11 A(I,J)=C(II,1)+DLOG(XNEW(II))-DLOG(SUMX) WRITE (6,14) ((A(I,J),J=1,MNP),1=1,N) 14 FORMAT (1X,5G18.6) CALL DSLSVD(A,S,V,NDIMAU,NDIMV,N,M,NP,£c1 40) EPS=1.D-6 SS=S(1)*EPS DO 60 J=1,M IF (S(J).LT.SS) GO TO 70 DO 50 1=1,M 50 V(I,J)=V(I,J)/S(J) 60 CONTINUE J=M+1 GO TO 90 1 44 70 WRITE (6,80) ((V(I,K),1=1,M),K=J,M) 80 FORMAT (1X,4G13.5) 90 IF (J.GT.M) GO TO 120 DO 110 K=J,M DO 110 I = 1 , M 110 V(I,K)=0.D0 120 MP1=M+1 CALL DGMULT(V,A(1,MP1),X,M,M,NP,NDIMV,NDIMAU,10) RETURN 140 WRITE (6,150) 150 FORMAT (' ERROR RETURN FROM DSLSVD') STOP END C C SUBROUTINE SKIP WILL READ THE LAGRANGIAN MULTIPLIERS FOR C THE PRIMAL CONSTRAINTS AS PRINTED BY GRGEG C IN THE SCRATCH FILE -GRGOUT. C SUBROUTINE SKIP(*) LOGICAL*1 RECORDO50) INTEGER*2 LEN LOGICAL EQCMP 100 CALL READ(RECORD,LEN,0,LNR,6,&200) IF (LEN .LT. 70 .OR. LEN .GT. 75) GO TO 100 IF (EQCMP (21, RECORD, '0LAGRANGE MULTIPLIERS')) RETURN GO TO 100 200 RETURN 1 END C C SUBROUTINE FREEN CALCULATES THE FREE ENERGY PARAMETERS FROM C THE GORDON AND MC BRIDE COEFFICIENTS. C SUBROUTINE FREEN IMPLICIT REAL*8 (A-H,0~Z) COMMON/CX/C(30,10)/SUX/S(30,10) COMMON/DX/N,M/HX/K,NN COMMON/EX/T(30),P(30) DO 10 L=1,NN C(1,L)=0.DO N1=N-1 DO 10 J=1,N1 JJ=J+1 10 C(JJ,L)=S(J,1)*(1.D0-DLOG(T(L)))-S(J,2)*T(L)/2.D0~ 1 (S(J,3)*T(L)**2.D0)/6.D0 -(S(J,4)*T(L)**3.DO)/1.2D1 2 -(S(J,5)*T(L)**4.D0)/2.D1+S(J,6)/T(L)-S(J,7) 3 +DLOG(P ( D ) WRITE(6,12) 12 FORMAT ('MATRIX OF FREE ENERGY COEFFICIENTS') WRITE(6,11) ((C(J,L),J=1,N),L=1,NN) 11 FORMAT (1X,4F12.7) RETURN END C C SUBROUTINE NPLUSK PERFORMS SENSITIVITY ANALYSIS SOLVING A SYSTE 145 C SYSTEM OF LINEAR EQUATIONS WITH SUBROUTINE DSLIMP(FROM UBC).THE C U(I,J) THAT IS , THE STOICHIOMETRIC COEFFICIENTS FOR THE C REACTIONS THAT TAKE PLACE, ARE PROVIDED BY THE USER. C SUBROUTINE NPLUSK IMPLICIT REAL*8 (A-H,0~Z) COMMON/AX/A(10,3 0)/BX/B(10)/XNE/DELT(31)/UX/U(30,30) COMMON/CX/C(30,10)/DX/N1,M/HX/K,NN/EX/T(3 0),P(3 0)/SX/SUMX DIMENSION DA(50,50), DT(50,50), DB(50), DX(50), 1 DRZ(50), IPERMO00) DIMENSION SUM(30),V(30),SUMB(30) N=N1+K NDIMAT=50 N2=N1-M DEPS=1.D-8 ITMAX=0 C NOW WE NEED TO FORM THE (N1+KXN1+K) MATRIX DO 83 1=1,N DO 83 J=1,N 83 DA(I,J)=0.D0 C NORMALITY AND ORTHOGONALITY CONDITIONS DA(1,1)=1.D0 DO 84 1=2,M DO 84 J=2,N1 DA(I,1)=-B(l) 84 DA(I,J)=A(I,J) C NOW, SUMMATION OF NUMBER OF MOLES DO 4 I=1,K LL=I+M DO 5 J=2,N1 5 DA(LL,J)=1.DO 4 DA(LL,N)=-1.DO C NOW, EQUILIBRIUM CONDITIONS DO 8 1=1,N2 SUMB(I)=0.D0 DO 7 J=1,N1 LL=M+K+I DA(LL,J)=U(I,J)/DELT(J) 7 SUMB(I)=SUMB(I)+U(I,J) 8 DA(LL,N)=SUMB(I)/SUMX C NOW CALCULATE THE RIGHT HAND SIDE VECTOR DO 95 L=2,NN M3= M+K DO 10 1=1,M3 10 DB(I)=0.D0 DO 11 1=1,N2 SUM(I)=0.D0 DO 9 J=1,N1 9 SUM(I)=SUM(I)-U(l,J)*(C(J,L)-C(J, 1 )) LL=I+M+K 11 DB(LL)=SUM(I) NRHS=L-1 C SOLVE THE SYSTEM CALL DSLIMP (DA, DT,DB,DX,DRZ , I PERM, N , NDIMAT, DEPS ,-NRHS , ITMAX) C WRITE OUT RESULTS WRITE (7,900) 900 FORMAT (' VARIATION OF NUMBER OF MOLES') WRITE (7,110) (DX(I),1=1,N) 110 FORMAT (1X,6G14.6) C CALCULATE THE NEW NUMBER OF MOLES DO 15 I=1,N 15 V(I)=DX(I)+DELT(l ) W=P(L) Z=T(L) WRITE (7,120) W,Z 120 FORMAT (' P=',F10.0,' T(KELVIN) = ',F10.0,1X,'SOLUTION * ) WRITE (7,130) (V(I),I=1,N) 130 FORMAT (1X,6G16.8) 95 CONTINUE RETURN END 147 Program COMP3 C C C C THIS PROGRAM SOLVES THE DUAL G.P. FOR THE CHEMICAL EQUILIBRIUM C PROBLEM. IT CALLS THE SUBROUTINE LIPSU2 TO GET A FIRST POINT C FOR THE OPTIMIZATION. THEN IT CALLS SUBROUTINES GRGIN AND GRG2 C WHICH SOLVE THE PROBLEM ITSELF.SUBROUTINE GCOMP CALCULATES THE C OBJECTIVE FUNCTION AND THE CONSTRAINTS, AND IS CALLED FROM C SUBROUTINE GRG2.SUBROUTINES GRGIN AND GRG2 ARE FROM UBC. C THE PROCEDURE IS REPEATED FOR DIFFERENT TEMPERATURES AND C PRESSURES. C SUBROUTINE FREEN CALCULATES THE FREE ENERGY COEFFICIENTS (G/RT) C USING THE DATA FROM GORDON & MC.BRIDE C c IMPLICIT REAL*8 (A-H ,0-Z) LOGICAL*1 BLANK /» */, SEMIC /';'/ C DIMENSION STATMENTS . LABELLED COMMONS ARE USED. C N=NUMBER OF CHEMICAL SPECIES +1 C M=NUMBER OF CHEMICAL ELEMENTS (ATOMS) + 1 C K=NUMBER OF PHASES C NN= NUMBER OF DIFFERENT TEMP. AND PRESSURE CONDITIONS. C NE=1 IF THE FREE ENERGY COEFFICIENTS AT TEMPERATURE T ARE C AVAILABLE. IF THEY HAVE TO BE CALCULATED BY C THE GORDON AND MC.BRIDE POLYNOMIALS, SET NE EQUAL TO C ANY INTEGER DIFFERENT THAN ONE. C NF=1 IF THE PROGRAM HAS TO PROVIDE A FIRST STARTING POINT. C IF YOU PROVIDE THE FIRST STARTING POINT, SET NF EQUAL TO C ANY OTHER INTEGER. C A(I,J),1=1,M,J=1,N = EXPONENTS MATRIX.A(1,1)=1 CORRESPONDS C TO THE NORMALITY CONDITION. C C ( J , L ) , J=1,N, L=1,NN = FREE ENERGY COEFFICIENTS OF SPECIES J C AT T(L) ,P(L) .SET C ( 1 , D = 0 . C T,P=EQUILIBRIUM TEMPERATURE AND PRESSURE. C XNEW(J),J=2,N = NUMBER OF MOLES OF SPECIES J AT EQUILIBRIUM. C FOR F=1, XNEW =1 (DUMMY SPECIES CORRESPONDING TO THE NORMALITY C CONDITION) C X(J),J=1,N = FOR SUBROUTINE LIPSU2, CORRESPONDS TO THE LINEAR C PROGRAM VARIABLES. FOR GRG, DUAL VARIABLES, EQUAL TO XNEW(j). C XMF(J),J=2,N MOLAR FRACTION OF SPECIES J C S(J,L) J=1,N-1, L=1,7 = COEFFICIENTS FOR GORDON AND MC.BRIDE C POLYNOMIALS FOR SPECIES J+1. C B(I),I=2,M VECTOR OF AMOUNTS OF ELEMENT I, THAT ARE CONSERVED. C B(1)=1, DUMMY. C F= OBJECTIVE FUNCTION. C c COMMON/AX/A(10,30) COMMON/BX/B(10) COMMON/CX/C(30, 10) COMMON/DX/N,M/HX/K,NN/GX/F/XL/L COMMON/EX/T(30),P(30)/SX/S(30,10) COMMON/FX/X(31) 1 48 DIMENSION XMF(30),H(30),Q(30) DIMENSION XNEW(31) C C READ IN DATA C READ(5,10) N,M,K,NN,NE,NF 10 FORMAT (615) C . :, K C CALCULATE FREE ENERGY COEFFICIENTS BY GORDON AND MC.BRIDE C IF (NE.EQ.1) GO TO 13 C C CALCULATE THE TEMPERATURE AT FIXED INTERVALS. C SX=55.D1 DO 20 L=1,NN SX=SX+50.D0 T(L)=SX 20 P(L)=1.DO N1=N-1 READ (5,40) ((S(I,L),L=1,7),I=1,N1) 40 FORMAT (5G13.7) CALL FREEN GO TO 400 C C IF THE C(I,L) MATRIX IS AVAILABLE BY OTHER SOURCES OTHER C THAN GORDON AND MC.BRIDE, READ TEMPERATURES, PRESSURES, & C THE C(I,L) MATRIX. C 13 READ (5,311) (T(L),L=1,NN) READ (5,311) (P(L),L=1,NN) READ (5,311) ((C(I,L),L=1 ,NN),I=1,N) 311 FORMAT (5F10.0) C C READ THE EXPONENT MATRIX AND THE B VECTOR C 400 READ (5,50) ((A(I,J),J=1,N),1=1,M) 50 FORMAT (13F4.0) READ (5,60) (B(I),1=1,M) 60 FORMAT (5F10.0) IF (NF.EQ.1) GO TO 500 C C IF YOU HAVE A GOOD STARTING POINT, THE PROGRAM SHOULD READ IT N C OTHERWISE, IT CALCULATES ITS OWN STARTING POINT. C READ (5,510) (XNEW(J),J=1,N) 510 FORMAT (5F10.0) TIME=SCLOCK(0.) GO TO 600 500 TIME=SCLOCK(0.) C CALLS SUBROUTINE LIPSUB2.IT PROVIDES A STARTING POINT FOR THE C OPTIMIZATION ROUTINE. CALL LIPSU2 XNEW(1)=1.D0 149 DO 100 J=2,N 100 XNEW(J) = X ( J - 1) + X(N) WRITE (6,12) (XNEW(J),J=1,N) 12 FORMAT (' FIRST STARTING POINT'/1X,6G16.8) DO 200 L=1,NN C OPTIMIZATION IS PERFORMED FOR NN DIFFERENT P AND T CONDITIONS C A SCRATCH FILE IS CREATED TO WRITE DOWN THE DATA NEEDED FOR THE C UBC SUBROUTINES GRGIN AND GRG2 TO PERFORM THE OPTIMIZATION CALL FTNCMD('ASSIGN 5=~DATA;') C C FIRST COME THE CONTROL CARDS. C WRITE (5,15) N,M,M 15 FORMAT (316) C C NOW WRITE THE LOWER BOUNDS FOR THE VARIABLES. C WRITE (5,25) 25 FORMAT ('LBV=') DO 115 I=1,N 115 H(l)=1.D-60 WRITE (5,35) (BLANK,I,H(l),I=1,N),SEMIC 35 FORMAT (6(A1,1 3,G10.2)) C C NOW, THE UPPER BOUND FOR THE VARIABLES. C WRITE (5,55) 55 FORMAT ('UBV=') DO 125 1=1,N 125 Q(I)=40.D0 WRITE (5,35) (BLANK,I,Q(l),I=1,N),SEMIC WRITE (5,65) 65 FORMAT ('QUAD') C C EPNEWT IS THE TOLERANCE FOR THE CONSTRAINTS. C WRITE (5,75) 7 5 FORMAT ('EPNEWT=1.D-6') C C EPSTOP IS THE TOLERANCE FOR THE OBJECTIVE FUNCTION. C WRITE (5,85) 85 FORMAT ('EPSTOP=1.D-6') WRITE (5,95) 95 FORMAT ('EPSBOUND=1.D-6') C C EPSPiV SHOULD BE OF THE ORDER OF THE LOWER BOUNDARY. C WRITE (5,105) 105 FORMAT ('EPSPIV=1.D-16') C C PRINTCTL MAY BE SET AT HIGHER VALUES (UP TO 4) FOR MORE C INFORMATION ON THE PROGRESS OF THE OPTIMIZATION. IF 'NOECHO' C IS WRITEN ON FILE -DATA JUST BEFORE THE CONTROL.CARDS, AND 1 50 C IF PRINTCTL =0, NO OUTPUT IS WRITTEN FROM GRG. C WRITE (5,61) 61 FORMAT (' PRINTCTL=1 1) DO 101 I=1,N 101 X(I)=XNEW(I) WRITE (5,107) 107 FORMAT(' X= ') WRITE (5,109) (X(I),I=1,N) 109 FORMAT (6G18.6) WRITE (5,108) 108 FORMAT ( 1OPTIMIZE 1/'GO'/'STOP') REWIND 5 C NOW GRGIN CAN READ INPUT FROM -DATA;GRG WILL WRITE ITS OUTPUT C INTO THE SCRATCH FILE -GRGOUT. CALL GRGINU1 ,&2) CALL GRG2(XNEW,F,&1) GO TO 2 1 WRITE (6,3) 3 FORMAT (' GRG HAS FAILED') 2 WRITE (6,205) T ( L ) , P ( L ) , F 205 FORMAT ('T(K)=',F10.0,'P(AT)=',F10.0,'OBJ.F=',G16.8) WRITE (6,215) 215 FORMAT ('NUMBER OF MOLES') WRITE (6,225)(XNEW(I),1=1,N) 225 FORMAT (1X,6G16.8) SUM=0.D0 DO 8 J=2,N 8 SUM= SUM+XNEW(J) WRITE (6,235) SUM 235 FORMAT (' TOTAL NUMBER OF MOLES=',G16.8) DO 9 J=2,N 9 XMF(J)=XNEW(J)/SUM WRITE (6,245) (XMF(J),J = 2,N) 245 FORMAT (1X,' MOLAR FRACTIONS',1X,6G16.8) 200 CONTINUE GO TO 23 23 TIME=SCLOCK(TIME) WRITE (6,19) TIME 19 FORMAT (' EXECUTION TIME =',F6.4) STOP END C C SUBROUTINE LIPSU2 CALCULATES A STARTING POINT APPROXIMATING THE C DUAL PROBLEM TO A LINEAR PROGRAM. THE LINEAR PROGRAM IS SOLVED C BY SUBROUTINE LIPSUB FROM UBC. C • SUBROUTINE LIPSU2 IMPLICIT REAL*8 (A-H,0"Z) DIMENSION TABLO(30,30),NVIN(30),NVOUT(30),BBOBJ(20), 1 UBOBJ(20),BRHS(20),UBRHS(10) COMMON/AX/A(10,30) COMMON/BX/B( 1 0) COMMON/CX/C(30,10)/GX/OPTIM 151 COMMON/DX/NVARS,NCONST COMMON/FX/X(31) N1=NVARS-1 M1=NCONST-1 NEQUAL=M1 MAXIM=0 IFOBJ=0 IFRHS=0 NCOLS=NVARS+1 NROWS=NCONST+1 C ZERO THE TABLEAU DO 200 J=1,NCOLS DO 200 1=1,NROWS 200 TABLO(I,J) = 0.DO C FIRST ROW IN THE TABLEAU IS THE OBJECTIVE FUNCTION SUM=0.D0 DO 101 J=1,N1 . JJ=J+1 TABLO(1,J)=C(JJ,1) 101 SUM= SUM+TABLO(1 ,J) TABLO(1,NVARS)=SUM TABLO(2,NVARS)=~1.DO TABLO(2,NCOLS)=-.0001DO DO 203 1=2,NCONST DO 202 J=1,N1 202 TABLO(I + 1, J) = A ( l , J + 1) TABLO(I + 1, NVARS) = 1.D0 203 TABLO(I + 1, NCOLS) = B ( l ) 260 NLVARS = NVARS CALL LIPSUB(TABLO, 30, NCONST, NLVARS, NEQUAL, MAXIM, 1 IFOBJ, IFRHS, TOL, NVIN, NVOUT,BBOBJ, 2 UBOBJ, BBRHS, UBRHS, &540) NP1=NLVARS+1 OPTIM =TABLO(1,NP1) IF (MAXIM.NE.1) OPTIM=-OPTIM DO 12 1=1,NVARS 12 X(I)=0.D0 DO 13 I=2,NROWS J = NVIN(I) 13 X(J)=TABL0(I,NP1) WRITE (6,320) (X(J),J=1,NVARS) 320 FORMAT (1X,5G18.6) RETURN 540 WRITE (6,560) 560 FORMAT (' STARTING POINT ROUTINE FAILED') DO 580 1 = 1 , NVARS 580 X(I) = 10.0D0 RETURN END C C SUBROUTINE GCOMP CALCULATES THE VALUES OF THE OBJECTIVE C FUNCTION AND OF THE CONSTRAINTS. C SUBROUTINE GCOMP(G,X) 1 52 IMPLICIT REAL*8 (A-H,0~Z) COMMON/AX/A(10,30)/BX/B(10)/CX/C(30,10)/DX/N,M/XL/L DIMENSION G(1),X(1),SUM(30) DO 20 1=1,M SUM(I)=0.D0 DO 10 J=1,N IF (A(I,J) .LE. O.DO) X(J)=0.D0 10 SUM(I)=SUM(I)+A(l,J)*X(J) 20 G(I)=SUM(I)-B(l) M1=M+1 DELT=0.DO DO 30 1=2,N 30 DELT=DELT+X(I) SUMAC DO DO 40 1=1,N DLX = -1.E50 IF (X(I) .GT. 0.0D0) DLX = DLOG(X(l)) 40 SUMA=SUMA+X(I )*(DLX + C ( I , D ) DLD = -1.E40 IF (DELT .GT. 0.0D0) DLD = DLOG(DELT) G(M1)=SUMA-DELT*DLD RETURN END C C SUBROUTINE FREEN CALCULATES FREE ENERGY PARAMETERS C SUBROUTINE FREEN IMPLICIT REAL*8 (A-H,0"Z) COMMON/CX/C(30,10)/SX/S(30,10)/DX/N,M/HX/K,NN COMMON/EX/T(30),P(30) DO 10 L=1,NN C(1,L)=0.D0 N1=N-1 DO 10 J=1,N1 JJ=J+1 10 C(JJ,L)=S(J,1)*(1.D0-DLOG(T(L)))-S(J,2)*T(L)/2.D0-1 (S(J,3)*T(L)**2.D0)/6.D0-(S(J,4)*T(L)**3.D0)/1.2D1 2 -(S(J,5)*T(L)**4.D0)/2.D1+S(J,6)/T(L)-S(J,7) 3 +DLOG(P(D) WRITE (6,7) 7 FORMAT (' MATRIX OF FREE ENERGY COEFFICIENTS') WRITE(6,11) ((C(J,L),J=1,N),L=1,NN) 11 FORMAT (1X,4F12.7) RETURN END APPENDIX B: EXAMPLES 154 PROBLEM 1: HYDRAZINE COMBUSTION Passy and Wilde,1968 Species = 10 Elements = 3 Phases = 1 G/RT = 47.8907 (1) Temperature= 3500 K G/RT = 47.8907 (2) Pressure= 51.2 atm. i Elements Bi Zi (1) Zi (2) 1 H 2.0 -9. 78896 -9.78903 2 O 1 .0 -15 .2128 -15.2127 3 N 1 .0 -13 .1000 -13.1000 j Spec i e s cj 6 j (1 ) 6 j (2) 2 H -6 .089 4.06 E-2 3.78 E-2 3 H2 -17 . 1 64 1 . 48 E-1 1 .47 E-1 4 H20 -34 .054 7.83 E-1 7.85 E-1 5 N -5 .914 1 .41 E-3 1 . 24 E-3 6 N2 -24 .721 4.85 E- 1 4.86 E-1 7 NH -14 .986 6.93 E-4 6.06 E-4 8 NO -24 .100 2.72 E-2 2.43 E-2 9 O -10 .708 1 .79 E-2 1.81 E-2 1 0 02 -26 .662 3.73 E-2 3.80 E-2 1 1 OH -22 . 1 79 9.69 E-2 9.74 E-2 T o t a l number of moles = 1 . 638 (1) T o t a l number of moles = 1 . 63799 (2) Summation of molar f r a c t i o n s at optimum (1) = 0.995586 Summation of molar f r a c t i o n s at optimum (2) = 1.000001 NOTES: (1) = Passy and Wilde, 1968 (2) = T h i s t h e s i s 155 PROBLEM 2 : METHANE AND WATER REACTION Perry, 1968 Species = 5 G/RT = 79.3605 (1) G/RT = 79.3605 (2) Elements = 3 Phases = 1 Temperature = 1095 K Pressure = 1 atm. Elements Bi Zi ( 1 ) Zi(2) 1 2 3 C H 0 2.0 14.0 3.0 -24.9739 -0.7922 -0.2007 -24.9738 -0.7922 -0.2007 Spec i e s Cj 6 j ( D 6 j ( 2 2 3 4 5 6 CO C02 H20 H2 CH4 -24.025 -47.413 -23.067 0.0 2.0847 51 74 3107 861 2 7942 , 1 722 ,5174 ,31 07 ,861 3 ,7943 , 1 722 X. T o t a l number of moles = 8.6558 ( 1 ) V T o t a l number of moles = 8.6559 (2) Summation of molar f r a c t i o n s at optimum Summation of molar f r a c t i o n s at optimum 1.0000 (1) 1.0000 (2) Notes: (1) = Perry,1968 (2) = T h i s t h e s i s PROBLEM 3 : WATER GAS REACTION D i n k e l and Lakshmanan,1975 Species = 4 Elements = 3 Phases = 1 G/RT = 90.4788 (1) G/RT = 90.4789 (2) Temperature = 1000 K Pressure = 1 atm. Elements B i Zi (1) Zi (2) C 0 H 1 .0 2.0 2.0 N.A N.A N.A -5.31206 -33.4888 -9.09497 Spec i e s Cj 6 j (1) 6 j (2) 2 3 4 5 CO H20 C02 H2 •37.4239 •50.3023 •70.8924 •1 6.7936 .505025 .505025 .494975 .494975 505044 505088 494978 494978 X T o t a l number of moles = 3.000000 (1) X T o t a l number of moles = 3.000088 (2) Summation of molar f r a c t i o n s at optimum Summation of molar f r a c t i o n s at optimum 1.000000 (1) 1.000031 (2) Notes: (1) = D i n k e l and Lakshmanan, 1975 (2) = T h i s t h e s i s . 157 PROBLEM 4 : RESPIRATORY SYSTEM Dembo, 1976 Species =31 G/RT = (1) G/RT = 90.4789 (2) Elements =12 Phases = 3 Temperature = N.A. Pressure = N.A. T h i s problem was solved with the untransformed pr i m a l v a r i a b l e s t , s c a l e d as in the r e f e r e n c e . Elements B i s c a l e d t i (1) s c a l e d t i 1 2 3 4 5 6 7 8 9 1 0 1 1 1 2 02 CO 2 N2 H+ OH-C l -Na + K+ HB1-HPp-HPr-Z (cha 4, 4, rge) 0013317 0022709 0024855 67 671973 008140 008092 005 000909 00088 001 19 0 2.517968 2.539194 7.657035 1.221926 467072 291600 283088 781697 1.787073 2.001670 496106 449689 6, 6, 508494. 525456 664300 1 93585 645168 1.314099 4.384049 2.651293 1.761015 2.017308 487580 408977 6, 6, k j Species Cj 6 j (1) 6 j ( 2 ) 2 02 -10.89 4.46 E-3 4.40 E-3 3 C02 - 7.69 1 .83 E-3 1.81 E-3 4 N2 -11.49 2.47 E-2 2.45 E-2 5 H20 -36.44 2.02 E-3 2.00 E-3 6 02 0 7.32 E-5 7.05 E-5 7 C02 0 7.37 E-4 7.10 E-4 8 N2 0 2.22 E-4 2.15 E-4 9 H+ 0 2.24 E-8 3.35 E-8 10 OH- 0 3.43 E-7 2.15 E-7 1 1 CL- 0 5.83 E-2 5.69 E-2 1 2 Na + 0 8.09 E-2 1.23 E-1 1 3 H20 -39.23 2.89 E+1 2.80 E+1 1 4 HC03- -21 .20 1 .40 E-2 8.74 E-3 15 H2C03- 0 6.72 E-21 6.48 E-21 16 C03- 6.25 2.18 E-5 8.80 E-6 1 7 HPp- 0 8.75 E-3 5.67 E-3 18 02 0 4.52 E-5 4.73 E-5 19 C02 0 4.55 E-4 4.76 E-4 20 N2 0 1 .37 E-4 1 .44 E-4 21 H+ 0 2.13 E-8 2.25 E-8 158 PROBLEM 4 : RESPIRATORY SYSTEM Continued k j Spec i e s Cj 6 j (1) 6 j (2) 3 22 OH- 0 1 .38 E-7 1 .44 E-7 23 C l - 0 2.34 E-2 2.48 E-2 24 K+ 0 5.00 E-2 5.00 E-2 25 H20 -39.23 1 .78 E+1 1 .87 E+1 26 HC03- -21.20 5.64 E-3 5.86 E-3 27 H2C03 0 4.15 E-21 4.35 E-21 28 C03- 6.25 5.69 E-6 5.89 E-6 29 HB1 0 3.10 E-3 3.31 E-3 30 HB102- -16.23 8.74 E-3 9.34 E-3 31 HPr- . 0 1 .20 E-2 1 .22 E-2 X, T o t a l number of moles, phase 1 = N.A. (1) X ^ T o t a l number of moles, phase 2 = N.A. (1) X ^ T o t a l number of moles, phase 3 = N.A. (1) X, T o t a l number of moles, phase 1 = .032700 (2) X 2 T o t a l number of moles, phase 2 = 28.1013 (2) X 3 T o t a l number of moles, phase 3 = 15.8586 (2) Summation of molar f r a c t i o n s , phase 1 = 1.000000 (2) Summation of molar f r a c t i o n s , phase 2 = 1.000014 (2) Summation of molar f r a c t i o n s , phase 3 = .999991 (2) Notes: (1) = Dembo, 1976 (2) = T h i s t h e s i s 159 PROBLEM 5 : COMBUSTION OF PROPANE (1) Species =10 G/RT = 777.6387 ( 1 ) G/RT = 777.6387 (2) Elements = 4 Phases 1 Temperature = 2200 K Pressure = 40 atm. Elements Bi Zi (1) Z i (2) 1 H 8.0 N.A. -11.4153 2 C 3.0 N.A. -20.1267 3 O 10.0 N.A. -15.8052 4 N 40.0 N.A. -11.6971 j Spec i e s cj 6 j (1) 6 j 2 H2 -15.6191 .0200735 .0204434 3 H -.7824 .0006540 .0006811 4 OH -19.7527 .0154000 .0152673 5 H20 -36.7180 3.9718994 3.971582 6 CO -30 .1221 .0815971 .0819268 7 C02 -49.5104 2.9184028 2.9180737 8 N2 -23.0912 19.9866573 19.986533 9 NO -20.5868 .0266857 .0269337 10 02 -24.9310 .0335845 .0338415 1 1 O -4.7912 .0004428 .000459 X. T o t a l number of moles = 27.0553973 (1) X T o t a l number of moles = 27.057600 (2) Summation of molar f r a c t i o n s at optimum = 1.000001 (2) Notes: (1) = D i n k e l and Lakshmanan, 1977 (2) = T h i s t h e s i s N.A. = Not a v a i l a b l e 160 PROBLEM 6 : CLAUS FURNACE 1 Bonsu, 1981 Species = 8 Elements = 4 Phases = 1 G/RT = N.A. (1) Temperature = 800K G/RT = 110.48998 (2) Pressure = 1 atm. i Elements Bi Z i (1) Zi (2) 1 S 0.3 N. A. -7.01247 2 0 1.0 N. A. -36.5171 3 H 2.0 N. A. -12.6496 4 N 3.760 N. A. -12.3856 j Spec i e s Cj X j (1) x j (; 2 S02 -76.38437 .02568 .02565 3 H2S -29.34294 .05136 .05139 4 H20 -60.55412 .28304 .28304 5 S2 -9.51304 0.01098 0.01098 6 S4 -19.09567 0.00013 0.00013 7 S6 -33. 14049 .00013 .00013 8 S8 -43.41504 .00000 .00000 9 N2 -24.30711 .62868 .62869 X. T o t a l number of moles at e q u i l i b r i u m = N.A. (1) X T o t a l number of moles at e q u i l i b r i u m = 2.99003 (2) Summation of molar f r a c t i o n s at optimum = 1.00001 (1) Summation of molar f r a c t i o n s at optimum = 1.000001 (2) Notes: (1) = Bonsu, 1981 (2) = T h i s t h e s i s . N.A. = Not a v a i l a b l e . The r e f e r e n c e g i v e s the composition in molar f r a c t i o n s . 161 PROBLEM 6B :CLAUS FURNACE 1 M c G r e g o r , 1978 Species = 8 Elements = 4 " Phases = 1 G/RT = N.A. (1) Temperature = 800K G/RT = 230.12109 (2) Pressure = 1 atm. Elements Bi Zi (1) Zi (2) 1 S 2.0 N.A. -6.28329 2 0 2.0 N.A. -36.6006 3 H 4.0 N.A. -12.7332 4 N 3.76 N.A. -12.4229 Spec i e s Cj Xj (1) Xj (2) 2 3 4 5 6 7 8 9 S02 H2S H20 S2 H2 S6 S8 N2 •76. •29. •60, -9, -16, -32, -43, -24, 38437 34294 55412 51 304 96621 48436 75302 3071 1 0455 0908 2192 0483 0002 0112 001 1 5833 04505 09009 21922 04783 00020 ,01130 ,00120 ,5831 1 X. T o t a l number of moles at e q u i l i b r i u m = N.A. (1) X. T o t a l number of moles at e q u i l i b r i u m ^ 6.44320 (2) Summation of molar f r a c t i o n s at e q u i l i b r i u m = 1.0001 (1) Summation of molar f r a c t i o n s at e q u i l i b r i u m ^ 1.000001 (2) Notes: (1) McGregor, 1978 (2) T h i s t h e s i s . The r e f e r e n c e g i v e s the composition i n molar f r a c t i o n s . 1 62 PROBLEM 7 : CLAUS FURNACE 2 Bennett,1979 Species =24 Elements = 4 G/RT = N.A. (1) G/RT = 230.13785 Phases = 1 Temperature = 800K Pressure = 1 atm. ) i Elements Bi Zi (1 ) Zi (2) 1 S 2.0 N.A. -6.28617 2 0 4.0 N.A. -36.6064 3 H 4.0 N.A. -12.7318 4 N 7.52 N.A. -12.4236 j Spec i e s Cj Xj (1) Xj (2) 2 S02 -76.38437 . 4440 E-1 3 H2S -29.34294 .9010 E-1 4 H20 -60.55412 .21 97 E 0 5 S2 -9.51304 N.A. .4692 E-1 6 S4 21.72258 .4424 E-20 7 S6 -32.48436 .5339 E-2 8 S8 -43.75302 .14 50 E-2 9 N2 -24.30711 .5827 E 0 1 0 NH3 -31.83038 .6922 E-8 1 1 S 20.62205 .2060 E-1 1 1 2 SH -2.09771 .4484 E-7 1 3 H2 -16.96626 .2040 E-3 1 4 H 18.08661 . 41 28 E-1 3 1 5 SO -27.24001 . 1 593 E-6 1 6 HO -17.51132 . 1 506 E-1 3 1 7 S03 -92.84698 .7925 E- 1 0 18 SN 11.90983 .5036 E-1 3 19 S20 -42.68584 .1514 E-2 20 NO -13.05837 .2386 E-1 5 21 S3 -13.37118 .41 39 E-2 22 S5 -24.89094 . 1 445 E-2 23 S7 -37.89094 .2186 E-2 24 0 -17. 17544 .4394 E-23 25 02 -25.98145 . 1983 E-22 X T o t a l number of moles at e q u i l i b r i u m = 6 . 4 5 0 0 (1) (3) X T o t a l number of moles at e q u i l i b r i u m = 6 . 4 5 2 0 7 ( 2 ) Summation of molar f r a c t i o n s at e q u i l i b r i u m = N.A. (1) summation of molar f r a c t i o n s at equilibrium=1 . 0 0 0 0 0 01 (2) 163 PROBLEM 7 : CLAUS FURNACE 2 Continued Notes: (1) = Bennett, 1979 (2) = T h i s t h e s i s . (3) = C a l c u l a t e d from Table 6-1 , p.141 ,Bennett,1979 Moles of product formed at 800 K from 100 moles of SH2 238 moles of a i r : Species Moles of product (1) Moles of product (2) H2 H2S S02 N2 H20 Sj .06 .010 29. 1 6 1 4.54 187.91 70.46 19.52 29.063 14.323 187.976 70.864 1 9.821 Others T o t a l .85 322.50 .545 322.602 No t e s * (1) Table 6-1, p. 141 , Bennett, 1979 (2) T h i s t h e s i s .
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A comparison of solution methods for the chemical equilibrium...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A comparison of solution methods for the chemical equilibrium problem Ruda, Margarita M. 1982
pdf
Page Metadata
Item Metadata
Title | A comparison of solution methods for the chemical equilibrium problem |
Creator |
Ruda, Margarita M. |
Date Issued | 1982 |
Description | This thesis deals with computing the equilibrium composition of a multiple species reacting mixture. When pressure and temperature are constant and the system is ideal, this is the chemical equilibrium problem. It is possible to approach this problem as the minimization of a non-linear objective function subject to linear equality constraints. The objective function represents the Gibbs' free energy of the system; the constraints refer to the conservation of the elements. Such a formulation corresponds to a "dual geometric program" which is related to another optimization problem known as the "primal geometric program". In those chemical equilibrium problems with many species, the "primal geometric programming" formulation includes less variables and constraints (inequality ones) than the dual formulation. We first compared the primal and dual formulations of the chemical equilibrium problem. Both formulations were solved with a Generalized Reduced Gradient code on seven examples. The primal formulation proved to be 30% faster than the dual for middle-sized problems (up to six simultaneous reactions). The code failed when trying to solve a dual problem of 24 species and 4 elements; but this same problem was easily solved when formulated as a primal geometric program. As the geometric programming theory includes sensitivity analysis, and we were, also interested in the effects of small changes of pressure and temperature on the optimal solution, we compared sensitivity analysis with re-optimization of the problem. Sensitivity analysis proved to be between 30% to 50% faster than re-optimization. It also yielded accurate results for the more abundant species when relatively small changes of temperature and pressure were operated. However, the equilibrium concentrations of trace species hardly matched those calculated by re-optimization. From these results we recommend the use of the primal formulation and of re-optimization to solve the chemical equilibrium problem, and we present a computer code that has been tested on a variety of examples. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-13 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0058911 |
URI | http://hdl.handle.net/2429/23429 |
Degree |
Master of Applied Science - MASc |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of Chemical and Biological Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1982_A7 R94.pdf [ 6.9MB ]
- Metadata
- JSON: 831-1.0058911.json
- JSON-LD: 831-1.0058911-ld.json
- RDF/XML (Pretty): 831-1.0058911-rdf.xml
- RDF/JSON: 831-1.0058911-rdf.json
- Turtle: 831-1.0058911-turtle.txt
- N-Triples: 831-1.0058911-rdf-ntriples.txt
- Original Record: 831-1.0058911-source.json
- Full Text
- 831-1.0058911-fulltext.txt
- Citation
- 831-1.0058911.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-0058911/manifest