(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Children's Library | Biodiversity Heritage Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

Full text of "Optimization of energy recovery systems"

THE OPTIMIZATION OF ENERGY RECOVERY SYSTEMS 



By 

JIGAR V. SHAH 



A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF 

THE UNIVERSITY OF FLORIDA 

IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE 

DEGREE OF DOCTOR OF PHILOSOPHY 



UNIVERSITY OF FLORIDA 
1978 



J 



Dedicated to 
my parents, Hemlata and Vinubhai Shah 



ACKNOWLEDGEMENTS 

The author wishes to express his sincere gratitude to Dr. 
A. W. Westerberg. Working as a graduate student under his super- 
vision was an invaluable experience. 

The author also wishes to express his appreciation to the 
department of chemical engineering at the University of Florida 
and to Dr. J. P. O'Connell and Dr. H. D. Ratliff of the supervisory 
committee. 

Thanks are due to the assistance provided by the department 
of chemical engineering at Carnegie-Mellon University and to the 
Computer-Aided Design Centre, Cambridge, U.K. 



TABLE OF CONTENTS 

Page 

ACKNOWLEDGEMENTS iii 

LIST OF TABLES vi 

LIST OF FIGURES vil 

LIST OF SYMBOLS ix 

ABSTRACT xii 

CHAPTERS: 

I INTRODUCTION 1 

II THE SYSTEM EROS 5 

II. 1. Background 5 

11.2 Data Specification 9 

11.3. Modeling Considerations 9 

11. 4. Deriving a Solution Procedure and 

Solving Model Equations 17 

11. 5. Starting the Problem: Finding a 

Feasible Point 20 

11.6. The Optimization Strategy 27 

II. 7. Special Uses of EROS 35 

II. 8. Results and Discussion 41 

III LOCAL AND GLOBAL OPTIMA 48 

III.l. Statement of the Problem 48 

II I. 2. The Lower Bound 49 

III. 3. The Upper Bound on the Lower (Dual) Bound... 54 



TABLE OF CONTENTS (Continued) 

Page 

III. 4. Finding the Best Upper Bound 57 

III. 5. Improvement of the Upper Bound 59 

III. 6. Procedure to Investigate Global 

Optimality 62 

III. 7. Algorithm 63 

III .8 . Examples 65 

III. 9. Discussion 75 

IV SYNTHESIS USING STRUCTURAL PARAMETERS: A PROBLEM 

WITH INEQUALITY CONSTRAINTS 76 

V CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER 

RESEARCH 83 

APPENDICES: 

A DESCRIPTION OF THE INPUT AND OUTPUT FEATURES 8 7 

B DESCRIPTION OF THE PROGRAM 102 

LITERATURE CITED Ill 

BIOGRAPHICAL SKETCH 114 



LIST OF TABLES 

Table Page 

1 Constraint: Representation for an Example Node 3 

in Figure 4 16 

2 The Incidence Matrix 21 

3 Solution Procedure for the Problem in Figure 2 22 

4 Additional User Specifications for the Problem 

in Figure 4 40 

5 Results for the Problem in Figure 9 42 

b Stream Specifications for the Example in Figure 10 44 

7 Results 46 

8 Stream Specifications and Problem Data 79 



LIST OF FIGURES 

Figur e Page 

1 Structure of an Optimizing System 7 

2 An Example Problem 10 

3 The Network Nodes 12 

4 Partitioning a Heat-Exchanger into Zones where 

Phase Changes Occur 15 

5 A Typical Cooling Curve 18 

fa Keeping Track of the Best Point 'e' 28 

7 Structure of an Optimization Algorithm 30 

8 Using an Existing Exchanger 36 

9 Reliability Analysis 38 

10 An Example Exchanger Network 43 

11 Example Problem 1 50 

12 The Behavior of the Objective Function for 

Example Problem 1 51 

13 Staged Processes 53 

14 Geometric Significance of the Dual 55 

15 Geometric Significance of the Upper Bound 56 

16 Subsystems for Example Problem 1 66 

17 Example Problem 2 68 

18 Subsystems for Example Problem 2 69 

19 Subsystems for Example Problem 3 72 



LIST OF FIGURES (Continued) 

Figure Page 

20 The Optimal Values of the Parameters 78 

21 The Discontinuity in Optimization 81 

22 The Optimal Network When Ignoring the Approach 
Temperature Requirement at Null Exchanger 2 of 

Figure 20 82 

23 Description of the Program EROS 103 



LIST OF SYMBOLS 

a = vector; constant for points on a hyperplane 

b = scalar; constant for points on a hyperplane 

C = the set of all newly violated constraints 

C = the violated constraint encountered first 

C = heat capacity 
P 

D = minimum allowable approach temperature 

E = the equation set without any inequality constraints present 
as equalities (heat and material balances only) 

f = scalar return function; may be subscripted 

f . = subproblem, subsystem, return function 
1 

F - flow rate 

o n ■ / 1 n+1. , , , n st 

G = matrix; first n rows are (g , . . .g ) and (n+l; row is 

(l,...l) 

g = vector constrained function; may be subscripted to repre- 
sent scalar elements. May also be superscripted as g 1 to 
represent the vector constraint function at a point i 

h = enthalpy 

H a set of (m+1 ) constraints where m is an integer 

H = a hyperplane 

I/.x = the set of i such that stream i is an input to subsystem j 

L = Lagrange function 

LB = lower bound 

0,.» = the set of i such that stream i is an output of subsystem j 



Q = duty of a heat-exchanger 

g = composite vector variable (x/u) : may be subscripted 

S = constrained variable set; may be subscripted 

T = temperature 

AT, = log mean temperature difference 

u = system decision variable: may be subscripted 

U = heat transfer coefficient 

UB = upper bound 

V = the current set of inequality constraints in the equation 

set (E - E„) 

c (J 



V .V . . V 
1' 2 



= subsets formed by removing one constraint at a 
time from H. (For example, V-^ is obtained by 
removing the first constraint in II.) 



V R = the set of constraints present in the equation set as 
equality constraints with the difference that their 
slack variables are used as search coordinates 

Vg = all the constraints in the system less the ones in V 
and V T 

V T = the set of constraints being held as equality constraints 

x = vector variable associated with interconnecting streams; 
may be subscripted 

y = multipliers to express a vector as a linear combination of 
other vectors; may be subscripted 

Z = objective function for the linear programming problem 
(maximization) 



Gree k Letters 

a = variable vector used in linear programming formulation; 
may be superscripted to represent scalar variables in 
linear programming formulation 

The scalar variable a and the subscripted variable 
u-jj may be used to represent a structural parameter 
(split fraction) . 



= objective function; may be superscripted as to represent 
its value at point i. May also be used as to denote a 
vector of 1! s 

A = vector of Lagrange variables; may be subscripted to repre- 
sent scalar elements. May also be superscripted as A to 
represent its value at point i 

o = slack variable; may be subscripted 

Mathemat ical Symbols 
H = intersection with 
'0' = null set 



Abstract of Dissertation Presented to the Graduate Council 

of the University of Florida in Partial Fulfillment of the Requirements 

for the Degree of Doctor of Philosophy 

THE OPTIMIZATION OF ENERGY RECOVERY SYSTEMS 

by 

Jigar V. Shah 

March, 1978 

Chairman: Arthur W. Westerberg 

Major Department: Chemical Engineering 

An emphasis on technical development oriented towards conservation 
of energy in the design of new chemical engineering processes and in the 
modification of the existing processes is as important as an emphasis 
on the search for alternate sources of energy. With the increasing 
sophistication in computer hardware, the design engineers are now able 
to assess innumerable options through the use of flowsheeting packages 
to maximize energy recovery. Most of the current flowsheeting packages 
essentially convert the well-defined input information about a process 
into a description about the output from the process. A next evolution 
in these packages is judged to be a system which will be very flexible 
in the input that is specified to it. In addition, it will possess 
an optimization capacity to yield the most desirable values for 
parameters left to it as degrees of freedom. This study describes 
a prototype for what a really convenient flowsheeting system ought to 
be. 



xii 



The program EROS (Energy Recovery Optimization System) is a flow- 
sheeting package for evaluating and optimizing the performance of simple 
networks of heat-exchangers. Using EROS one can set up an arbitrary 
structure of heat-exchangers, stream splitters and mixers. Stream 
flow-rates and entry and exit temperatures may be specified, free, or 
bounded above and/or below. Phase changes may be allowed to occur. 
Requiring no more user input (such as initial guesses), the program 
gathers together the modeling equations and appropriate inequality 
constraints. It then develops solution procedures repeatedly in the 
course of optimizing, initially to aid in locating a feasible point 
(which need not be specified by the user) and in the final stages to 
take advantage of tight inequality constraints to reduce the degree 
of freedom. The program is written in standard Fortran and is cur- 
rently operable on IBM 360/67. 

The development of this program is relevant from the point of view 
of evaluating and optimizing energy recovery systems because the subject 
of choosing an optimal structure has recently come under increasing 
attention. A tool such as EROS will prove very useful to make a more 
thorough analysis of the candidate structures chosen by the process of 
synthesis at the penultimate and final stages. Existing energy recovery 
schemes may be reevaluated with the intention of making changes to make 
them more efficient. EROS can also perform quick reliability studies 
accounting for fluctuations in stream flow-rates and temperatures, 
a common problem during start-up and periodic failures. 

As is the case in the optimization of most chemical engineering 
systems, a global optimum cannot always be guaranteed. An attempt 



to recognize whether the discovered optimum is a global one, in a 
large system, is not a trivial task. Usually the problem may be 
side-stepped by making several optimization runs from different 
starting points and considering the best result as the required 
optimum. A systematic procedure is presented in this study to 
investigate global optimality, and the results on application to 
simple examples prove it to be quick and effective. 

The use of EROS can be extended to perform process synthesis 
via structural parameters. However, the presence of inequality 
constraints inherent in the system can easily force the optimization 
to stop short of the desired result. This observation has not been 
reported before and attention is brought to it in this dissertation. 



CHAPTER I 
INTRODUCTION 

Most of the existing flowsheet ing packages for chemical engi- 
neering processes are based on Sequential or Simultaneous Modular 
Approaches. In these approaches each unit is modeled by writing a 
computer subroutine which converts the input stream and equipment 
parameter values into output stream values. Systems based on 
Sequential or Simultaneous Modular Approaches are relatively easy 
to build but the penalties paid are the lack of flexibility in the 
definition of the problem and the requirement for well-defined user- 
specifications . 

Equation solving approaches, as followed in EROS (and in Leigh, 
Jackson and Sargent (1974), Hutchison and Shewchuk (1974) and Kubicek, 
Hlavacek and Prochaska (1976)) present an alternative way to treat 
flowsheets. The flowsheet is represented as a collection of non-linear 
equations which must be solved simultaneously. In following an Equation 
Solving Approach, the user can specify many of the values for both unit 
inputs and outputs. The unit equipment parameters are then calculated 
to give these desired transformations of inputs to outputs by the unit; 
in other words, the unit is designed to meet these requirements. Since 
EROS also contains an optimization capability the rather striking 
advantage is that variables for whose values the user has no preference 
for are treated as the degree of freedom by the system. EROS incor- 
porates the general optimization strategy as outlined in Westerberg 



and deBrosse (1973) and demonstrates the applicability and effec- 
tiveness of their algorithm. 

The two important problems in the design of energy recovery 
systems are to choose the configuration, and, given a configuration, 
to choose the design parameters and operating variables. A recent 
review on the effort directed to choosing a configuration can be found 
in Nishida, Liu and Lapidus (1977). In choosing a suitable configura- 
tion the general trend has been to evaluate networks using the heuristic 
of setting the minimum allowable approach temperature to 20°F, whereas 
the economics often advocate a smaller value. In addition, when a 
stream is split, the need for finding the optimal value for the split 
fraction has been ignored. Grossman and Sargent (1977) optimized 
several heat-exchanger networks and found considerable savings 
(sometimes as much as 25%) . 

The problem of optimizing a heat-exchanger network to obtain 
the most suitable values for operating variables has been considered 
in the works of Westbrook (1961), Boas (1963), Fan and Wang (1964), 
Bragin (1966) and Avriel and Williams (1971). Typically each design 
problem is formulated and solved for as an optimization problem. 
Many investigators (Hwa (1965), Takamatsu, Hashimoto and Ohno (1970), 
Henley and Williams (1973), and Takamatsu et al. (1976)) have combined 
both the problems, choosing a configuration as well as the operating 
variables, and formulated it as an optimization problem. 

All the methods mentioned for optimizing over the operating 
variables require the problem to be cast into a mathematical format. 
EROS precludes this need because of its capability as a flowsheeting 
package . 



In Chapter II of the dissertation a description of the program 
EROS is provided. The data specification and modeling considerations 
for the system are discussed here while a user's guide to EROS is 
presented in the appendices. The derivation of a solution procedure 
is illustrated with the help of a simple example. If the user has not 
provided a feasible starting point, EROS has the capability of dis- 
covering one. The algorithm used is a modified version of that pre- 
sented in deBrosse and Westerberg (1973). The optimization strategy 
based on Westerberg and deBrosse also merits attention in this study 
because of its usefulness to this application. The different applica- 
tions of EROS and the results on optimization of 10 problems of varying 
sizes are also shown. 

The program EROS does not guarantee a global optimum because 
of the non-linear and multimodal nature of the problem. In Chapter 
III an algorithm is presented to establish, often quite quickly, 
whether the discovered optimum is global or not. The strategy is 
based on finding both a lower (dual) bound for the optimum of the 
structured system and an upper bound on this lower bound. In this 
chapter the effectiveness of the algorithm is then illustrated by 
applying it to three small heat-exchanger network problems. 

With a convenient evaluation and optimization package such as 
EROS, it is very tempting to use it to perform process synthesis via 
the use of structural parameters. Using this approach several 
alternate configurations are imbedded into one main structure 
which, on optimization, yields the required test configuration. 



Several authors, as mentioned earlier in this section, have already 

experimented with this idea. However, the problem to be solved has 

to be formulated very carefully, and not enough attention has been 

paid to this iact. In Chapter IV a discussion regarding this subject is 

presented. 

The conclusions from creating a system such as EROS and recom- 
mendations for further investigations form Chapter V of this disserta- 
tion. A guide to the use of the program is given in the appendices. 
The aspects of EROS presented are the input data format, the inter- 
pretation of the output, a typical computer printout, and a description 
of all the subroutines. 



CHAPTER II 
THE SYSTEM EROS 

II. 1. Background 

In order to evaluate and optimize heat-exchanger networks it 
is desirable to have a program which on being given information about 
the configuration and stream properties, yields all the required 
information about the optimal network. Since the program will perform 
several different tasks, it would be very attractive and in many 
instances necessary for i t to possess the following features. 

a) An ability to set up solution procedures. 

The program should eliminate or reduce com- 
putational recycles, choose the decision variables 
and discover the order for calculating the various 
unknown variables. 

b) An ability to obtain a feasible starting point. 

Often, locating a feasible starting point for 
optimization is not a simple task. In order to save 
users the time and trouble necessary to find a feasible 
starting point, the program should be capable of per- 
forming such a task on its own. 

c) Efficient optimization routines. 

These would be required for selecting the 
optimal values for the decision variables by 
searching over their full range. 

However, the optimization of a system such as this one raises 
certain problems. The objective function is highly non-linear and 
multimodal. Also, if phase changes are allowed, continuous derivatives 
cannot be obtained. These criteria force the use of a search algorithm 
such as the complex method. In having resigned to the use of the 



complex method for optimization, further demands can be made to improve 
the efficiency of the approach for optimization. If in the process of 
optimization several constraints are violated, one remedial action is 
the use of penalty functions, but this modification is inefficient and 
it increases the number of iterations required for convergence. 

Hence, with regard to optimization, few additional features would 
be deemed attractive. 

d) The ability to rederive a solution procedure. 

When a constraint is violated, the program 
should be able to modify the equation set and re- 
derive a solution procedure so that the optimization 
may be continued. 

This strategy will lessen the number of iterations required 
for convergence as compared with the penalty function method. However, 
it is essential that the saving in computer time thus incurred com- 
pensates for the extra time required in rederiving the solution 
procedure . 

e) The use of restriction as a solution strategy. 

In the presence of a large number of decision 
variables, some of them are set to zero and optimiza- 
tion is performed with only the remaining ones as 
search coordinates. Optimization is considered 
complete when the Kuhn-Tucker (1951) optimality 
conditions are satisfied with respect to the decision 
variables held at zero value. This strategy aids 
considerably in reducing the number of iterations 
during optimization. 

EROS is the author's attempt at developing an optimization 
program which incorporates all the features discussed above. Figure 1 
illustrates the general structure of the EROS system. 

The approach taken is to model each unit in a heat-exchanger 
network functionally by writing overall material and energy balances. 
The unit models themselves may be very complex internally, but the 




Objec tlve 

Function i <u,x(u)) 



Fig. 1. Structure of an Optimizing Syster 



net effect at the flowsheet level is that each unit satisfies overall 
heat and material balances. In the current version of EROS only simple 
models are used. Using the functional equations, the modeling of such 
a system is only at the flowsheet level, and a solution procedure or 
the order in which these equations are to be solved is found along with 
the degree of freedom to be chosen. The solution procedure sought is 
one that will eliminate, if possible, computational recycles at this 
level . 

The above approach is useful because the equations are solved 
repeatedly as an inner loop to an optimization program. As illustrated 
in Figure 1, the optimizer directs all the activity. Its primary func- 
tion is to adjust the decision variables to improve the objective 
function 0. For this system is the annualized cost of the equipment 
plus the annual cost of buying the utilities needed such as steam for 
the purpose of heating. 

To evaluate 0, the optimizer supplies the block labeled "Solve 
Model Equations" with the values it wishes to try for the decision 
variables u. The remaining problem variables x(u) are then obtained 
by solving the model equations. With u and x(u) values available, 
constraint violations are checked, and if some are violated, they are 
identified to the optimizer. Assuming none are violated, the units 
in the system are sized and an annualized cost is evaluated. The 
optimizer notes this cost and changes the decision variable values, 
with the aim of reducing 0. This calculation sequence is repeated 
many times during a typical optimization. If constraints are violated, 
special action is taken, which tor this system will result in a modi- 
fied set of model equations and a need to rederive a solution procedure 



for them (this approach is based on the optimization strategies in 
deBrosse and Westerberg (1973) and Westerberg and deBrosse (1973)). 
The modified complex optimization algorithm (Umeda and Ichikawa (1971)) 
is used in searching for the decision variable values, u. 

11. 2. Data Specifi cation 

Adequate data must be supplied to the computer to define the 
problem. A problem definition will require the following input. 

1) the flowsheet structure 

a) connectivity of segments and units 

2) the unit data 

a) unit type 

b) any equipment parameter specifications 

3) the stream data 

a) flow and temperature specifications 

b) physical property data 

c) film heat transfer coefficients 

d) materials specifications and other cost 
parameters 

4) the segment data 

a) associated stream identifier (i.e., what 
stream this segment is a part of) 

b) any specifications imposed on flow and 
temperature 

5) general user specifications 

6) guessed set of inequality constraints to be held 

11 .3. Modeling Considerations 

The modeling of a network will be illustrated with regard to 
the example in Figure 2. 

The network comprises a single hot process stream H which is 
split and used to heat two cold process streams C and C„. It then 
merges to its exit conditions. Streams C and C„ are heated further 
by steam utility streams S and S, ; . The network has four heat 



JO 



I (11,1 




11(5.) 



I - 500 



I - 500 



C - 1 U - 25 U, - 20 0. - 0, - 50 

p 2 3 5 6 

For Stream II T - T - 600° Ural of Vaporization ^800 

Cost multiplier per unit flov rate - 4 

For Stream III Tj a - T^ - 600° neat of Vapomation =• 800 

Cost multiplier per onlt flea; rate - 0.2 

Constraints: 2500 =: F,,F.. s 7500 , 300° £ T.-.T.. * 500° 



2' 3 



12' 15 



Aim: To minimize i. i "Y Cost, + 0.4 1 + 0.2 F 

1-2. I?M 

vhere Cost, - 35 (/area.) ' . Cost, and Area. 
1 - i 1 1 

axe associated u-lth exchanger 1. 



Fig. 2. An Example Problem 



11 



exchanger (or heater and cooler) units, 2, 3, 5 and 6, one stream 
splitter unit, 1, and one mixing unit, 4. These units may also be 
referred to as nodes. All the heat exchangers are assumed counter- 
current. The streams have been broken up into segments, of which 
there are 16 overall. For example, stream C 1 enters node 2 as 
segment 11. It exits and proceeds to unit 5 as segment 12, and 
finally leaves the system as segment 13. The naming scheme should 
now be evident. The 3 basic units used are shown in Figure 3. The 
unit models are written functionally by writing overall material 
and energy balances. 



Unit 



F = aF (I) 

F 3 = (1 - a)F 1 (2) 

h 2 = h x (3) 

h 3 = h 1 (4) 



F 4 « F 2 (5) 

(6) 
h 2 F 2 + h ll F ll = h 4 F 4 + h 12 F 12 (?) 



(8) 

(9) 
h 3 F 3 + h 14 F 14 = h 5 F 5 + h 15 F 15 (10) 



6 4b 

h,F, = h.F + h F, (12) 

6 6 4 4 5 5 



12 



cream VyHcat-Exclirfiige 



d 



Fig. 3. The Network Nodes 



13 



5 F g - F y (13) 

F 13 - F 12 (14) 

h 7 F 7 +h 12 F l2 = h 8 F 8 +h l3 F 13 (15) 

(16) 

(17) 
h 9 F 9 +h 15 F 15 = h 10 F 10 + h 16 F 16 (L8) 

In addition to these 18 equations, the associated inequality 
constraints and the equipment sizing and costing relations can be 
written. 

The basic inequality constraint is that at no point in the 
exchanger should the hot stream temperature equal or fall below that 
of the cold stream. Referring to the heat-exchanger in Figure 3a, 
this constraint is usually expressed as 

T > T, + D, where D is minimum allowable approach 

1 4 

temperature . 
T 2 > T 3 + D 

However the temperatures could cross over internally and these 
constraints may not be adequate to detect it, particularly when a 
stream passes through a phase change. A check should therefore be 
made at several points along the exchanger to prevent a "crossover" 
of temperatures. 

The final set of constraints indicate that a positive heat 
transfer must occur, that is, the hot stream must be cooled and the 
cold stream heated. These are 



14 



T„ < T. and T- < T. . 
2 1 3 4 

All Che constraints associated with a typical exchanger such 
as the one shown in Figure 4 are listed in Table 1 with their ap- 
propriate code so that they can be precisely identified by numbers. 
No constraints are written for the splitter and mixer units. 

The sizing calculation for an exchanger is to evaluate its 
area. This calculation can be very involved, but for design purposes 
may in fact be simplified by assuming film coefficients based on the 
fluid, whether it is heating or cooling and whether it is boiling or 
condensing (see Perry and Chilton (1973)). The exchanger may, again, 
for simplified design purposes, be considered to operate in zones as 
indicated in Figure 4. Each zone is then sized and costed using the 
appropriate film coefficients and temperature driving forces. 

To place a crude cost on the exchanger one might use an equation 
of the form (see Guthrie (1969)). 

COSt = f M f P (aAm) 

Where f is a materials factor, f D a pressure factor and A the 
M ' 

area of a zone within an exchanger. The zones are sized and costed 
differently because of the different types of heat-exchange duty. 
The terms a and m are constants, with m being about 0.6 to 0.8. 
Constant costs are assumed for the splitter and the mixer units. 

The last source of equations is the evaluation of physical prop- 
erties. The system must be able to convert from stream temperature (and 
vapor fraction for a pure component in two phase region) to enthalpy 
and vice versa. A cooling curve may be provided for the stream 



15 



Hoc Stream I 



II . Ill 



T 5 T nrc T « T i 



Scg 3 F 3 .T 3 .H 3 

Cold Str 



r mn " Dew polnt tem r crat:ure . hoC »trem 

r - Bubble point teinperature , hot strca 

T - new point temperature, cold stream 
T - But.ble point ceinperature , cold stre 



T 5 


= 


l^T' 1 " 


urc of 


the 


cold 


stre.-E 


i at a 






locaClo 


3 where 


the 


hot a 


tream 


teraper- 






aCure 1 


T DFM 










T 6 


- 


Terapera 


ure of 


Che 


cold 


streau 


ac a 




locatlo 


l where 


the 


hot s 


cream 


temper- 






ature 1 


T 
BPH 










I, 


- 


lenpcra 


ure of 


the 


hot s 


trcjn 


at a 






locatlo 


i where 


the 


cold 


strem. 


tempor- 






ature i 


T 
I) PC 










T 8 


- 


rempera 


ure of 


the 


hot s 


tre.TO 


al a 




locitloi 


where 


the 


cold 


strean 


temper- 






ature 1 


T 











Minimum allowable approach temperacu 



Fig. 4. Partitioning a Heat-exchanger into Zones where Phase 
Changes Occur 



16 



TABLE 1 
CONSTRAINT REPRESENTATION FOR AN EXAMPLE NODE 3 IN FIGURE 4 





Exterior 






Cggggent 


Constraints 


Code 




approach 


T s T + D 


(NODE 


* 10) + i = 31 


approach 


Tjit + D 


(NODE 


* 10) + 2 = 32 




V T 2 


(NODE 


• 10) + 3 * 33 




T^T 3 


(NODE 


* 10) + U = 34 


F 1IR = lower bound 
LUi for How F 


F 1 2F LB- 


(NODE 


* 10) + 5 = 35 


F. _. = upper bound 
1U ^ for flow F, 


F 2 F 

OB 1 


(MODE 


» 10) + 6 = 36 


?-., u - lower bound 
^ for flow F 


F 2 F 
r 3 3 LB 


(NODE 


* 10) + 7 = 37 


F 3DB = upper bound 
J for flow r 


F a F 
3 OB 3 

Interior 


(NODE 


» 10) -t 8 = 38 


Type 


Constraints 


Representation 


approach 


VW d 


(NODE 


* 1000) + 1 = 3001 


approach 


W T 6 + D 


(NODE 


* 1000) + 2 = 3002 


approach 


. VW ° 


(NODE 


* 1000) + 3 = 3003 


approach 


W° 


(NODE 


* 1000) + A = 3004 


In addition there could 


be constraints asr 


jciated with any stream seement. 


Example: Segaent 3, No 


le 3 






Coaiuent 


Constraints 


Code 




H,. R - lower bound 
J for enthalpy h\ 


K 3 iH 3LB 


- (SEG « 
+ 1) 


1000 + NODE » 10 
^-3031 


H ^rm = u PP er bound 


^OB*^ 


- (SEG " 

4. 2\ 


1OC0 + NODE * 10 



for enthalpy H 



17 



if the stream is assumed to be at a constant pressure. Figure 5 
illustrates a cooling, curve, where T and T are the dew and 
bubble point temperatures respectively. 

Properties such as thermal conductivities and densities 
should also be provided if the film coefficients are to be deter- 
mined from correlations. If for design purposes typical values for 
film coefficients are to be used, these properties will not be needed. 

II .4 . Deriv ing ^-h£ Sol uti on Pr oced ure and_ 
Solving Mod el Equat ions 

Consideration will now be given to developing a solution 
procedure and then solving the example problem. First, the system 
must gather together the necessary equations, or at least establish 
their structure, so that a solution procedure may be prepared. The 
desired solution procedure should eliminate all recycle loops in the 
computations if possible or minimize their number. 

The initial solution procedure will ignore all but the user 
specified inequality constraints. Thus the system sets up the 18 
heat and material relations shown in the last section. Assuming that 
the user has requested that constraints 55 and 33 be included, the 
additional relations required are 



F 7 = F 7LB + °55 



T 3 - T 5 + a 33 



(19) 

(20) 



where F is the lower bound on F 7 . 
/LB / 

The inequality constraints have been converted to equality by 
the slack variables o , and Q„„ which are required to be nonnegative. 



18 




T " Bubble point tecrpcratur 
T - Dew point temper a Cure 



Fig. 5. A Typical Cooling Curve 



19 



When the solution procedure is derived o and a will be required 

to be decision variables with an initial value of zero. In this way 

F-, and T„ will be forced to equal to F and T respectively. 
/ 3 /Lis j 

The relation (20) is in terms of temperatures rather than 
enthalpies. Hence the following relationships should also be added 

T 3 - f(h 3 ) (21) 

T 5 = f(h 5 ) (22) 

The system can implicitly account for equations (20) , (21) 
and (22) by the single expression 

h 3 = f(h 5 ,o 33 ) . (23) 

The incidence matrix can now be created. However its size can 
be significantly reduced. Note that a large number of equations simply 
equate one variable to another. Equations (3), (4), (5), (6), (8), (9). 
(13), (14), (16) and (17) are precisely of this form. These equations 
will automatically be satisfied if the values of variables so equated 
are stored in a common storage location. Hence these equations may 
be deleted and the two variables occurring in each one of them may be 
merged to a single one. 

Many of the variables in the incidence matrix are in fact 
specified. Let the following be specified in data input for the 
example in Figure 2 

Flows F i' F ii' F i4' F 7LB 

Enthalpies ^V'? ' h 8'V h 10' h ll ' h 13> h 14' h 16 



2U 



These specified variables along with the slack variables o and <J 
(required to be decision variables) may be eliminated in the incidence 
matrix. The resulting and much reduced matrix is illustrated in 
Table 2. 

A modification of Lee, Christensen and Rudd (1966) algorithm 
is applied to determine the solution procedure. The solution proce- 
dure that results on the application of this algorithm to the incidence 
matrix of Table 2 is shown in Table 3. The variables listed are cal- 
culated from the corresponding equations in the order indicated. Note 
that there is an iteration loop involving the single 'tear' variable 
F (from steps A to 9). F_(=F.) appears in equation (12) and the 
iterations between steps 4 and 9 are continued until the value of F. ; 
guessed at step 4 is essentially the same as the value of F^ calculated 
from equation (12) in step 9. 

The execution of the solution procedure, that is, calculation 
of the variables from the equations assigned to them, is termed 
"Solve Model Equations" in Figure 1. Corresponding to every unit, a 
subroutine is required to calculate any variable involved in the heat 
and mass balances of the particular unit. These subroutines may be 
supplied by the user for more sophisticated models. 

II. 5. Starting the Problem: Finding a Feasible Poin t 

If the user has not provided any information to aid in obtain- 
ing a feasible starting point, a modified version of an algorithm by 
deBrosse and Westerberg (1973) is used. As mentioned earlier, a signif- 
icant effort would be required on the part of a user to provide a 
feasible starting point if computational loops are involved in solving 



21 



TABLE 2 
THE INCIDENCE MATRIX 



Variables 

Equations 


tx 

;- 


(1) 




(2) 




(7) 




(10) 




(11) 




(12) 




(15) 


X 


(18) 




(19) 


X 


(23) 





-<f 


IT\ 


















c\j 






LPv 






t 1 




r-i 




IT\ 


xf 1 


-*0 


a 


X 














X 




X 












X 



X XX 

XX x 

XX XX X 



22 



TABLE 3 
SOLUTION PROCEDURE FOR TOE PROBLEM IN FIGURE 2 



Decision Variables 
Variable 



i. y=v 



n , _ . it 



2 - ^ 

3. h 5 



4. Guess Fp(=F. ) 



5. a 



6. F (=F ' 



7. h, 



8. F, 



'• *- 2 <=y 



10. H 



15 



ii. f 9 (=f 10 : 



55' 33 
Equation 

(19) 

(15) 
(23) 

(1) 
(2) 



HI) 

(12) 
(10) 
(18) 



23 



model equations. Computational loops are almost inevitable in complex 
networks. However, the user does have the option of providing a feasible 
starting point. 

The deBrosse and Westerberg (1973) algorithm uses an indirect 
approach. It hypothesizes that a subset of constraints has no feasible 
region and then attempts to verify the conjecture. If successful, the 
subset is identified as infeasible and obviously no feasible point 
exists. If unsuccessful, either a new hypothesis can be generated or 
the algorithm has indirectly found a feasible point. 

The feasi ble pq i n t __al go r i t hm 

In order to demonstrate the feasible point algorithm, the 

following definitions are necessary. 

E = The equation set without any inequality constraints present 
as equalities (heat and material balances only) . 

E = The current equation set (E., + inequality constraint (s) 
present as equalities with slack variables) . 

C = The set of all newly violated constraints. 

C. = The violated constraint encountered first. 

V = The current set of inequality constraints in the equation 
C set: (E c - E Q ). 

H = A set of (m + 1) constraints. 

V ,V , . . .V . = subsets formed by removing one constraint at 

a time from H (for example, V-i is obtained by 
removing the first constraint in H) . 

'0' = null set. 

Degree of freedom (as defined for this work) = The number of original 
variables (excluding slack variables) in problem less number of con- 
straints in E . 



24 



Procedure 'A' represents the following sequence of operations. 

1) Set all the decision variables to a value of 
zero . 

2) Solve the model equations. 

3) Check for constraint violations. 



Step 1 Set E = E„, V = '0' 
c c 

Step 2 Determine the solution procedure based on struc- 
tural considerations of the equation set. If a 
solution procedure can be determined, go to Step 4. 

Step 3 Attempt unsuccessful. Exit. 

Step 4 Execute Procedure 'A' . If any constraints are 
violated, go to Step 6. 

Step 5 Feasible starting point. Exit. 

Step 6 Set E = [E ,C, ] 
c c 1 

Step 7 Determine the solution procedure. If a solution 
procedure cannot be found, go to Step 11. 

Step 8 Execute Procedure 'A'. If the problem does not 
prove solvable, go to Step 11. If no constraints 
are violated, go to Step 5. 

Step 9 If the degree of freedom is not zero, go to 
Step 6. 

Step 10 Set H = [V ,C ], i = 1. Co to Step 12. 

Step 11 Set i = 1, H = V 

Step 12 Set E = [E„,V.], V. + V from Step 10. 
c l l c 

Step 13 Determine the solution procedure. If it cannot 
be determined, go to Step 17. 

Step 14 Execute Procedure 'A' . If the problem cannot be 

solved go to Step 17. If no constraints are violated, 
go to Step 5. 

Step 15 If C Pi H f '0', go to Step 17. 

Step 16 If the degree of freedom = 0, go to Step 10. If 
not, go to Step 6. 



25 



Step 17 Set i - i + I. If i > m + 1, go to Step 3 (the 
set of constraints in H do not allow for a feasible 
region). If i <• m + 1, go to Step 12. 



The main departures from the deBrosse and Westerberg (1973) 
algorithm are as follows: 

1) In this work, the set H is restricted to containing 
less than or equal to (n + 1) constraints where n 
represents the number of decision variables. 

2) Both the sets V and H are created in a different 
manner. In this work, one constraint is added at a 
time to build up the set V , and H is created either 
when there are already n constraints in V and an 
additional constraint is violated, or when a struc- 
tural inf easibility occurs. 

3) "No Solution" options (in Steps 8 and 14 where the 
problem is not solvable) are not treated in as 
rigorous a fashion as in deBrosse and Westerberg 
(1973) . The method here is the same unless several 
problems corresponding to the same "hypothesis" 

set H lead to no solution in Step 14. The algorithm 
here may terminate unsuccessfully at Step 17 whereas 
that of deBrosse and Westerberg might continue with 
alternate and reduced sets H. This option is quite 
complex and was never found necessary in EROS. 
Hence it was never included. 

The algorithm can now be applied to the example shown in Figure 

2. (Note that this problem involves 8 equations in 10 unknowns; thus 

two degrees of freedom exist at the start.) 



E Q = [(1), (2), (7), (10), (11), (12), (15), (18)] 

(all the equations from Table 2 except the last two, 
(19) and (23)) 

1) E = E.., V = '0* 

c c 

2) Solution procedure derived 

4) Procedure 'A' executed 

Constraint Violation = 33. 



26 



E = [E n , 33] 
c U 

Solution procedure derived 

Procedure 'A' executed 

Constraint Violation = 25 

Degree of freedom = 1 

E c = [E Q5 33, 25] 

Solution procedure derived 

Procedure 'A' executed 

Constraint Violation = 22 

Degree of freedom = 

H = [33, 25, 22], V = [33, 25], V 2 = [22, 33], V 3 = [22, 25; 

E c = [E Q , V 2 ] - [E Q , 22, 33] 

Note: V = [33, 25] is V in the previous iteration and 
it need not be considered again 

Solution procedure derived 

Procedure 'A' executed. Constraint Violation = 55 

C H = '0' 

Degree of freedom = 

H = [22, 33, 55], V = [22, 33], V 2 = [55, 22], V 3 = [55, 33 j 

E c = [E Q , V 2 ] = [E Q , 55, 22] 

Note: V = [22, 33] is V in the previous iteration and 
is eliminated from consideration 

13) Solution procedure derived 

14) Procedure 'A' executed 
No constraint violation 

5) Feasible starting point. Exit. 



2 7 



II . 6. The Optimization Strategy 

The optimization strategy is modeled after the algorithm 
presented in Westerberg and deBrosse (1973) . The algorithm is in- 
voked once a feasible point is available. 

The sets of inequality constraints are divided into three sets. 

V = The set of inequality constraints being held as 
equality constraints . 



V = The set of constraints present in the equation 

set as equality constraints with the difference 
that their slack variables are used as search 
coordinates. 

V = The set of all remaining constraints. 

V = [V m ,V„] (Note: V is as defined in the previous 
c T R c 

section. ) 

Solution procedures are modified as inequality constraints are 
moved from one set to another. Adding constraints in the set being 
held tends to aid the optimization process by reducing the dimension 
of search space for what is usually a marginal or no added burden in 
solving an enlarged set of equations. 

As the optimization proceeds, the values of all variables, in- 
cluding all slack variables, are stored for the point that yields the 
best value for the objective function. Hence, even when V^ is changed, 
optimization can be and is started at the best point discovered up to 
that moment. This modification makes a significant improvement to 
the Westerberg and deBrosse (1973) method. Figure 6 illustrates the 
typical dilemma faced by the Westerberg and deBrosse optimization 
algorithm when stepping from a current best point, point 'e', through 
one or more inequality constraints to point 'f'. At point 'f' the 



:-h 




Fig. 6. Keeping Track of the Best Point 'e' 



29 



constraint g 1 is detected as being violated. The algorithm will 
respond by changing the solution procedure so that the slack 
variable for g becomes a decision variable. The other decision 
variable will be either x or x . There are several options now as 
to where the optimization may be started. The algorithm could hold 
x. or x„ (whichever is selected as the decision variable) at its 
current value and find the point where a is zero, leading to point 
P, or P„ respectively. Alternatively it could attempt to locate P„ 
by searching along the direction leading from 'e' to 'f'. All of 
these options can, and often do, lead to a next point which has a 
higher and thus worse value for the objective function. By saving 
all the variable values for the best point, the search can always 
start, even after developing a new solution procedure, from that 
point, that is from point 'e'. This change reduces cycling because 
a change in the solution procedure cannot lead to a point that is 



The actual search strategy used is the modified complex 
method (Umeda and Ichikawa (1971)). The complex method is consid- 
ered suitable because gradients are not required. The treatment 
of phase changes creates discontinuities in first derivatives of 
functions . 

The optimization Al gori thm 

The notation followed is the same as that in the last section. 

Figure 7 gives the structure of this algorithm. 

Step 1 Obtain a feasible starting point. 

Step 2 Set V T = V , V R = '0'. Apply the Kuhn-Tucker (1951) 
conditions as follows. Perturb each slack variable 



30 




Unsuccessful Exit 



A - Horro/ll mode loop 

B - Find scorch coordin.ntrs, nbnonal node loop 

C - fVCul n to ii.-imal n iitc 



Fig. 7. Structure of the Optimization Algorithm 



31 



corresponding to constraints in V^, one at a 
time, away from zero. If an improvement in 
the objective function results, move the cor- 
responding constraint from V to V . 

A. If any solution procedure fails numerically during 
this step, go to Step 11. 

B. If on completion of this step, the number of con- 
straints in Vj is equal to the number of constraints 
in V c , and the degree of freedom is zero, go to Step 
18. 

C. Otherwise continue. 



The optimization is back in its normal mode if return is to Step 3 
or Step 7. 

Step 3 Perform an optimization (using the complex algo- 
rithm) . Use as search variables the slack variables 
corresponding to inequality constraints in V R along 
with any additional problem variables still retained 
as decision variables. Search only over nonnegative 
values of the slack variables. 

A. If at any point the solution procedure fails numeri- 
cally, go to Step 11. 

B. Otherwise, when the optimization algorithm exits 
normally, continue. 

Step 4 An optimization has just been completed. 

A. If one or more constraints are violated, go to 
Step 7. 

B. Otherwise continue. 



Step 5 Set V™ , , = V,,,. Apply the Kuhn-Tucker conditions 
as done in Step 2. 

A. If any solution procedure fails numerically, go to 
Step 11. 

B. Otherwise on completion continue. 



m 



Step 6 A. If V_ = V, T , .,, go to Step 18. 
1 Iold 

B. Otherwise repeat from Step 3. 

One or more constraints not in the current set V has been violated. 

c 

Step 7 A. If the degree of freedom is not zero, set 
E = [E ,C ,] and go to Step 9. 

B. Otherwise find the set of constraints R from V 
such that each constraint in R could, from struc- 
tural considerations of the equations, be traded 
for C . If R = '0', go to Step 10. 

C. Find the constraint C, in set R having the 
largest value for its corresponding slack 
variable. If that value is zero, go to 
Step 10. 

D. Otherwise continue. 

Step 8 Replace C with C in V . 

Step 9 For any slack variable corresponding to a constraint 
in V^ and having a value of zero at the current best 
point, transfer the corresponding constraint from Vp 
to V T . Set V = [V D) V n> ], E = [E„,V ] and develop 

1 . . c K , * C UC 

a new solution procedure. 

A. If a solution procedure cannot be developed, go 
to Step 11. 

B. Otherwise return to Step 3. 



The degree of freedom is zero and either (a) the values of all the 

slack variables in V.. are zero or (b) the set R is empty. 
K 

Step 10 Set H = [V ,C ], set i = 1, and go to Step 12. 

The current set of equations are (or appear to be) numerically 
inconsistent. 

Step 11 Set H = [V ] , set i = 1. 



33 



Step 12 Set E = [E„.V.], V. + V (this test is only relevant 
cOiiC 
if entry was originally from Step 10 above) . Set V R 

equal to the set of all constraints in V. whose slack 

variable values are greater than zero. Set V.,, equal 

to all remaining constraints in V,. 

A. If V R = '0' , go to Step 13. 

B. Otherwise, determine a new solution procedure. 

1. If one can be found, go to Step 14. 

2. Otherwise, go to Step 17. 

If entry to Step 12 was via Step 10 originally, a degenerate problem 

is at hand as V and C, (N T7 + 1 constraints) are all satisfied in a 
c 1 Vc 

subspace of dimension N TI . 
Vc 

Step 13 Determine a solution procedure. 

A. If one cannot be determined, go to Step 17. 

B. Otherwise, perturb, one at a time, the slack 
variables corresponding to the constraints in 

V 

1. If, on any perturbation, constraints are 
violated such that HOC i '0', go to Step 17. 

2. If, on any perturbation of a slack variable, 

an improvement is obtained in the objective 

function, put the corresponding constraint 

into V„ . 
K 

C. If, on completion of the perturbations in Step B, 
the set V_ = '0', go to Step 18. 

D. Otherwise continue. 

Step 14 Perform an optimization. Use as search variables 
the slack variables corresponding to inequality con- 
straints in V R along with any additional problem 
variables still retained as decision variables. Use 
the solution procedure derived in Step 13. Search 
only over nonnegative values of the slack variables. 



34 



A. If at any point the solution procedure fails 
numerically, go to Step 17. 

B. When the optimization algorithm exits normally 
and if a constraint violation occurs, then 

1. If HflC / '0', go to Step 17. 

2. Otherwise, HOC = '0'. Go to Step 7. 

C. Otherwise when the optimization algorithm exits 
normally and no constraints are violated, continue. 

Step 15 Set V,„ , , = V™. Apply the Kuhn-Tucker conditions 
, fo]d T ' 

(as done xn Step 2) . 

A. If any solution procedure fails numerically, go 
to Step 17. 

B. Otherwise continue. 

Step 16A. If V T = V Told , go to Step 18. 
B. Otherwise, return to Step 14. 

■k 

The current solution procedure failed or led to an immediate 
constraint violation of the constraint dropped in set H. 

Step 17 Set i = i + 1. 

A. If i > m + 1, go to Step 19. 

B. Otherwise repeat from Step 12. 

Normal exit. Optimization attempt apparently successful. 
Step 18 Optimization complete. Exit. 

Abnormal exit. Optimization attempt aborted. 

Step 19 Optimization failed. Exit. 

Application of the strategy to the Example Problem in 
Figure 2. 



3 5 



1) Feasible starting point. V = [55,22] (from the 
last section) . 

2) V R = [22] 

3 and 4) Optimization with the complex method. 

7) Constraint violated = 35, R = [22], and o J0 > 0. 

8) Replace constraint 22 with constraint 35. 

9) V R = [35], V T = [55]. E c = [E Q ,V R ,V T ]. 

3) Optimization with the complex method. 

4) Optimization complete. 

5 and 6) V m 1 , = [55]. V^ = [55] after Kuhn-Tucker test. 
iold T 

18) Optimization successful. Exit. 
At the optimum, the following values resulted: 

F = F. = 7350, F„ = F c = 2650, T, „ = 500°, T, c = 300°, 
2 4 3 5 12 15 

F 7 = 0, F n = 2500, T. = 392°, T = 423° and = 3538 

/ y 4 5 

II. 7. Special Uses of EROS 

Using existing Heat-Exchangers . The program may be used to 
analyze a network where some of the exchangers are already specified. 
A special effort must be made, however, to make use of these exchangers. 
It is assumed that these exchangers are available at no cost. In 
Figure 8, an exchanger with an area of A has been specified. On 
analysis, however, it is discovered that an exchanger with an area A 
is required at that particular site in the network. In the program, 
the costs assumed for different conditions of A and A are shown in 
Figure 8. The physical significance of 1) is that a by-pass will be 



3 fa 






Specified 
Exchanger 


Require 
Exchang 


1) 

2) 


If Aj« A, 

If A, > A, 


T03t 
Cost 



Cost associated 
with an exchange 
having the axes 



Fig. 8. Using an Existing Exchanger 



37 



used in the specified exchanger. 2) implies that an exchanger with 
area = A„ - A must be purchased in addition to the exchanger already 
present . 

Reliab ility analysis . The program has been extended to permit 
its use in quick reliability studies. The reliability studies will 
be demonstrated with the help of an example. Figure 9a represents 
a network that is operational under normal conditions. Cold process 
streams C, , C , and C are heated to their final temperatures with 
the help of a hot process stream H,, and a flue gas H . The flow 
rate and the outlet temperature of stream H. are undefined but are 
required to be within specified bounds. Assume that two abnormal 
occurrences take place separately, for certain periods of the year, 
namely , 

1) Stream H„ is unavailable. 

2) Heating of stream C is no longer required. 

The aim now is to find an optimal network such as the one 
shown in Figure 9a, fully provided for to meet the contingencies with 
the aid of by-passes in the exchangers or with the aid of auxiliary 
exchangers. The designer is permitted to allow a change in the 
flow of stream H and its outlet temperature provided they stay within 
their specified bounds. In the case of failure mode 2), a cooling 
utility stream may be used to cool stream C„ to 230° so that it may 
be recycled to be heated to 250° if it is desired to maintain a 
semblance to the normal operation. (Flow rate of C„ can be allowed 
to vary between and 70,000 in this instance.) 



38 



F - '•0,000 
T - 150° 



F » 20,000 
T - 300° 



1 > 



F * 70,00 

T - 2 30° 



F - 70,000 
I - 30° 



H3 O— Q 



T « 400 T - 250 T - 150 



o- — © 



(b) 



o 



o 






-6 



(c) 

* Uaspeclfled 

Fig. 9. Reliability Analysis 



39 



Figure 9b represents the network when stream H is unavailable, 
and, Figure 9c when stream C„ is not required to be heated. Additional 
user specifications to those shown in Figure 9 are presented in Table 4. 
In order to find the optimal operating system, networks in Figure 9a, 9b 
and 9c are optimized together. The objective function is given by 

= c h 1 f h 1 + c h; f h; + c h-' f h^ + c uc" F uc" + 

5 

I (cost of exchanger area at site i) 
i=l 

where C. and F. represent the cost coefficient and the flow rates of 
l l 

stream i, respectively. The cost coefficient C. should reflect the 

l 

expected fraction of the year that the network is in the particular 
state being represented. For example, for the problem in Figure 9 
it is assumed that the networks in Figure 9a, 9b and 9c are opera- 
tional 77%, 11.5%, and 11.5% of the time in a year, respectively. 

Hence if C is . 1 , then C , and C „ are 0.015 each. 
1 1 1 

The cost of exchanger area at a site will be illustrated with 

the help of an example. At the site 2, the exchangers of different 

sizes required in Figure 9 are A , A , and A ,,. Assume that A , is 

smallest and A ,, the largest area. 

The cost of exchanger area at site 2 is defined as 

35[(A 2 ,) U ' 6 + (A 2 - A 2 ,)°- 6 + (A 2 „ - A 2 )°- 6 ] 

This manner of costing areas in doing reliability analysis 
appears to be a good formulation of the real system. If, because of 



4(J 



TABLE 4 
ADDITIONAL USER SPECIFICATION FOR THE PKOBEEM IM FIGURE 9 



Stream 



UC" 



Flow Rate 



Lower Bound on Flow 



50,000 50,000 50,000 



Upper Bound on Flow 200,000 200,000 200,000 100,000 



Inlet Temp 



600 c 



600° 600° 100° 



Outlet Temp 



150° 



Lower Bound on Outlet Temp 190° 100° 190° 150 C 



Upper Bound on Outlet Temp 600° 600° 600° 150 c 



Cost Coefficient 



0.10 0.015 0.015 0.008 



Unspecified. 

Let U. represent the heat transfer coefficient for exchanger I 

u i ; V = 7CJ0 

U 2 = V " l, 2" " Ul1 • 2l 

U 3 = Uy U,» - 562.5 

U, -• U ' = U " --, 700 

U •' = 225 



4.1 



some departure from the normal mode, more exchanger area is required 
at a particular site, then one must pay for the auxiliary exchanger. 
If the exchanger area required is more for the normal mode, then an 
auxiliary exchanger is already in effect. The results obtained on 
optimization of the system in Figure 9 are shown in Table 5. 

II. 8. Results and Discuss ion 

A typical application of EROS to a heat exchanger network is 
illustrated by the example in Figure 10, the stream specifications 
for which are shown in Table 6. Stream I is a flue gas and streams 
IV and VIII are utility streams. The flow rates for these streams 
are not defined but are required to be within specified bounds. 
Some of the streams in the network are also multicomponent and are 
characterized by a dew and a bubble point. Thus the network is 
analyzed to ensure that minimum approach temperature violation does 
not occur inside the exchanger owing to discontinuities at the dew 
and bubble points. In fact, at the optimal solution for the example 
in Figure 10, the minimum approach temperature constraint at the bubble 
point of stream III is operative inside exchanger 5. The program can 
also evaluate variables that already exist, fully or in part. For 
example, exchanger 3 was assumed present with an area of 1500 units. 
The objective function is defined as 

= ^ cost coefficient * Flow Rate 
streams I, IV, VIII 

13 



+ 35 I (A. + 10)°- 6 - ID " 6 
1=1(^7, 8) 1 



42 

TABLE 5 
RESULTS FuR THE PROBLEM IN FIGURE 9 



Exchanger Area 

1 131.50 

2 30.01 

3 413.24 

4 103.24 
2 ' 41 . 80 
3* 462.31 
4* 105.06 
l" 131.50 

2 36.10 

3 0.00 
4" 64.67 
5" 0.00 



*'„ = 170,341 F , = 180,878 F - 50,148 

1 "l "l 

1 11 n 

T = T = T = 190 



= $23451. 25/yr. 



43 



T - 350 
F - 80,000 



Flue Gas I 
F-LB-1000 
I - 700 




ST - STKAM, CU - COOL It JO WATER 
LB - LCi.'R BOUND 
I) - till). ALLOWABLE Al'l'i.OACH TCMP. 

iat - I-'.t: .".:.■.!. ly iii-i.a :a:;. AiTRuACi ti:mp. 
* - Active constraints at oi-rrn': 



Fig. 10. An Example Exchanger Network 



44 





M 




£ 




"i 




3 




w 


sO 


a 




H 


w 




h-l 


Uh 


Sj 


y 



tt 



o o o 



•> O m O 



O O O 



o o o 



rf 



o o o 



o o o 



H O u^ 10 



o o o o 



o o o o 



O O O C) o 



r-i O O O O 



O O o 



o o o oo 





; -, 


o 


>i 




^~ 




O 


1 H 


1 H 


>H 


< ) 


N 


Ff 


C J 


1- 1 


t 1 
1 1 


n 
f ' 


u> 


( ) 


• ■: 


•4 


o. 


3 


IrH 


«; 


F 


£j 


■=■:• 


-<, 


E 


In 


n< 


n 


1:1 




U 


■i; 


en: 








t ) 


< ) 


i'i 


a; 


( ) 


H 




(-, 


fi 


o 

p-, 




d 




^ 








1 "* 


ft. 


n 




H 


M 




( ) 


W 


t i 


M 


t 1 


pi 





o o o o 



tr\ U^, U~\ 



u O O 



O O O O 



o o o 



rH C3 CO 



O O O O O o 



ti £ 



£ 3 



t-;' 




:-5 


s 


C 1 


f 1 


m 


m 


ic; 


f .' 


-i 


|:J 


1 •; 




(J 


»'; 



45 



where A is the area associated with exchanger i. In the calculation 
of cost associated with an exchanger, the relation 

cost = 35[(A + 10) 0,6 - 10 0,6 ] (1) 

is used in preference to 

cost = 35 A " 6 (2) 

because whenever the area of an exchanger currently at zero value, is 
increased, the cost for the exchanger as calculated from (2) increases 
abnormally compared to the change in cost for the rest of the system. 
The modification as shown in relation (1) dampens this ill-behaved 
effect. 

There are 8 decision variables for this problem and the optimum 
results after 328 iterations from an infeasible starting point. The 
stopping criterion is a 1 x 10 difference between the worst and best 
objective function values in the current set of points retained by 
the complex algorithm. The value of at the optimum is $152,109/yr. 

Results for 10 examples are shown in Table 7. In all the 
examples the feasible point results in very few iterations. It may be 
observed that the time required for the rewriting of solution proce- 
dures after finding a feasible point is relatively small as compared 
to the time taken for function evaluations during optimization. The 
maximum ratio of these two times occur in Example 6, but it is still 
less than 1/3. This observation indicates that the penalty paid for 
rewriting solution procedures whenever constraint violations occur is 
indeed very small. 



4h 



a-. 


in 


n 


-:r 


■"'< 


Xj 


O 

-J ' ". 


r 1 | _ 


-1- 
o 


;:; 


* 


CI 


en 
o 


H 


03 


CO 


I'"- f . 


OJ 

O) '*< 
I-- 


Oi 


cO 


00 


00 

* co 




o 

co 


o 


CJ 


m 


u"> 


m 


en 


CI 


r- CO 
O 

■-I O 
r-l 


r i 
CO 

r- 

c-j 


en 


cn 


CM 


m 


•4" 


-4 


o 

o ° 

o 


CJA 

a. co 


C) 

o 
o 

O) 


<t 


<r 


en 


i~- 


J 


2 


o 

o ° 

o 


2 J 


C l 
CO 

-rf 

^3 


en 


" 


- 


LO 


en 


Cl 


I i 
o 


c i m 
-cr 


OJ 


-J 


- 


<r 


c. 


CM 


o 

o ° 

o 


£ ^ 


CN 

CO 


- 


en 


c. 


-* 


Ol 


-.t 


o § 

o 


CM 

en 


r - 


f: 


6 

en 


n 

a 

C/i 


u 

c 

o 
w 


X> 

a 

O 

a 

n 


C 

o 

n. 

X) 
■H 

Ln 

a 

[-« 

O 
O 


IJl 

u -j, 
d -o 
•a c 
a o 
a u 
01 u 

%4 <--\ If: 
Pj Xl ^ 

C u) <u 
a a g 

u ^ H 
p 

O OJ 

in i-i 

CJ 


13 TJ 

a c 

f^ o 

u 

a ivj 
c h 


t-J 


U) 

-a 
a 
o 
u 

a 



^ 3 a 



47 



The figures shown for time in Table 7 are those required on 
IBM 360/67. The cost per second of CPU time is about 1.4 cents and 
the longest run (Example 10) costs $8.56 for complete execution while 
Example 1 costs $0.40. The size limitations are 50 process stream, 
25 nodes and 150 stream segments in the current version of the program. 
The program is fairly well tested and gave satisfactory results when 
used for sixteen different problems set up by students in a recent 
design course. 



CHAPTER III 
LOCAL AND GLOBAL OPTIMA 

III .1 . State ment of tlie_ Prob lem 

In the course of optimizing a chemical engineering system 
using a conventional minimum seeking algorithm, one should ask and 
then attempt to discover if the solution found is indeed a global 
one. If one is willing to try, then a common strategy is to re- 
start the optimization at several different initial points, and, 
if all or most lead to a single best point, that point is conjectured 
to be the global optimum. 

A second common strategy is to use one's intuition to claim 
that only one optimum is likely. This strategy is particularly 
dangerous when the optimum is at the boundary of the search region 
where portions of the system effectively disappear because the flow 
through them is zero. Often one models the capital cost of a process 
unit by an equation of the form 

cost = C (Throughput) 

This form is not convex and is particularly troublesome at zero 

throughput, where the slope is infinite for cost versus throughput. 

If a unit becomes zero in size, then a small positive perturbation 

in its throughput has an apparent positive infinite effect on cost, 

an effect usually large enough to trap the optimization algorithm. 

One can of course (and should) modify the cost equation to reduce 
this problem. 

48 



49 



These strategies may help but still do not guarantee that one 
has found the global optimum. The purpose of the work leading to this 
paper was to produce a reliable method to determine if an optimal solu- 
tion is either a global or a local optimal solution. It was hoped 
that the technique would be computationally effective, in particular 
for heat exchanger networks, a class of problems which commonly dis- 
plays local optima. 

The optimization problem relating to a heat-exchanger network 
is to minimize the annualized cost, 0, of heat-exchangers and utilities 
subject to the approach temperature inequality constraints and the 
heat and material balances representing equality constraints. 

The simple network in Figure 11 (due to Grossman and Sargent 
(1977)) has only one degree of freedom and the network cost can be 
explored by choosing different values for the temperature T... Figure 
12 represents the plot of costs, 0, vs. T„. At both the points A and 
B, the Kuhn-Tucker optimality conditions are satisfied as any slight 
perturbation away from each of these points results in an increase 
in 0. However, the point B clearly represents a local optimum. It 
is quite likely therefore that an optimization algorithm will not find 
the global optimum. 

III. 2. The L ower B ound 

To permit decomposition of a system structure, the primal or 

overall system optimization problem may be written as 

n 
Minimize = f(q) = \ f . (q . ) 
j=l J ] 

s.t. g. = x. - t.(q.) = ieO,.. j=l,2,...n 

o 1 ± XX (j) 

q.ES. j=l,2,...n (PI) 



50 



F - 10,000 
T - 600° 



F - 10,000 
T - 900° 



F - 10,000 



<^ — 6- 



C - 1 I) - V. - 200 

P 12 

Area - Q/u jlT , Cost ■ 35 x Area ' 6 

Cost, refers to exchanger 1. 
Aim: To minimize the objective function i, where i » Cost + Cos 



Fig. 11. Example Problem 1 



51 



i 














300 ■> 










8 


286.3 


200 . 














189. 3 


A 












100 . 





























100 




ibo 


200 
T 3 


250 




300 



Fig, 12. The Behavior of the Objective Function for 
Example Problem 1 



52 



(It is important to note that the constraints defining S must be 
complete enough to guarantee the boundedness of each subproblem j.) 

This representation of a system, such as the one in Figure 13, 
is based on Lasdon (1970) , but a procedure for equality rather than 
inequality constraints is stressed as shown in HcGalliard and 
Westerberg (1972). This stress follows from an interest in decomposed 
system structures, as against, say, allocation of resources. 

The problem is assumed to have an optimal feasible solution. 
The Lagrange function for problem (PI) is 



L = - I j, A.[x. - t.(q.)] 



j=l ie0 (j) 

which may be rewritten as 

n n 

L= J ff.(q.) + / V * 4 t 4 (q,) ~ I A i x J = LfMi'h) 



I [f.CqJ + I \h ( ^ ~ I A i x i ] = .V 

j-l 1 J ie0 (j) lel (j) 



J J 



This Lagrange function is separable in q, and a Lagrange problem may 
be defined which is equivalent to the following set of subproblems, 
one for each subsystem in the original system: 

Minimize f.(q.,A) j=l,2,...n 

q.ES. ' ] ( P2 > 

.] J 

The solution to the Lagrange problem equals the sum of the solutions 

of (P2), i.e., 

j=l J 



53 




_i 



Fig. 13. Staged processes 



54 



A dual function can now be defined as 

n . 
h(A) = I f (A) 
j = l J 

The geometric significance of the dual h(A) is illustrated in 
Figure 14 for a simple system like the one shown in Figure 13. The 
dual h(A) is the intercept of the line with slope A. Note that this 
line is a supporting hyperplane at the point (0°,g°). The geometric 
significance of the dual has been treated in Lasdon (1970).. The 
value of the dual is always less than or equal to the value of the 
primal optimum. Hence h(A) can be considered a lower bound on the 
global optimum. 

III. 3. The Upper Bound on the Lower (Dual) Bound 

Theorem 1 : If in the space of vs . g_, where £ is n- 
dimensional, a polytope (Geoffrion (1969)) is formed in the hyper- 
plane passing through n+1 support points, the value of at the 
intersection of this polytope with the axis £ = provides an upper 
bound on the lower bound if the point of intersection is contained 
inside or at the boundaries of the polytope. 

Discussion : Figure 15 illustrates the ideas underlying this 

1 2 
theorem. g and g are support points to the graph of vs. g. The 

1 2 
line connecting g to g^ is a polytope formed in the hyperplane pass- 
ing through these n+1 = 2 support points. This polytope or line 

12 
between g and g intersects the vertical axis g = 0, and thus the 



55 



PrLnAl Optical Valui 



Supporting hyporplane 
u-ltli slope X 




Support Point 



g - vVV 



Fig. 14. Geometric Significance of the Dual 



56 



Upper Bound (UB) 




Fig. 15. Geometric Significance of the Upper Bound 



57 



theorem says this point of intersection represents an upper bound on 

the lower bound. 

1 2 
Proof : All the support points of a graph (such as g and g ) 

are contained on the surface of a convex hull formed by the support 

hyperplanes for this graph. Assuming that the polytope intersects 

the axis g = in its interior or at its boundaries, a hyperplane 

intersecting the axis g = above this point must intersect a part 

of the graph, that is, it cannot be a support hyperplane. Thus every 

support hyperplane must intersect the axis g = on or below the point 

the polytope intersects this axis. Consequently the dual cannot exceed 

this value of (=VB) and UB is an upper bound to the dual. 

Il l . 4 . Finding th e Best Upper Bound 

Given support points (0 ,g ) i=l, 2 , . . . p, pin+1, the problem is 
to determine the following. 

(a) Find, if possible, a set of n+1 points that 
yields an upper bound. Posed differently, the 
problem is to find n+1 points sue!) that the 
polytope formed by them intersects the axis 

g = 0. 

(b) If several sets of n+1 points qualify to pro- 
vide an upper bound, find the set that yields 
the lowest value for the upper bound. 

Any point inside or at the boundaries of the convex polytope, 

specifically the point g = 0, can be obtained by a convex combination 

of (n+1) points forming the polytope. This result has been stated and 

applied in Director et al. (1978). Thus problem (a) can be formulated 

as a feasible point problem in linear programming. If a feasible 

solution exists (no artificial variables present in the solution at 

a nonzero level), then indeed the set of n+1 points does provide an 

upper bound . 



58 



In order Co solve problem (b) , if a "price" could be associated 
with each point such that the objective function reflects the value 
of the upper bound corresponding to a set of n+1 points, then it too 
is a linear programming problem. The solution to it, if one exists, 
would yield solutions to both problems (a) and (b) simultaneously. 
This "pricing" is indeed possible. 

For every point on the polytope formed by n+1 points, can 
be represented by the linear relationship 



T 



= £ a + b 



(1) 



where _a and b are constants. The value of at the intersection with 
the axis ^ = is 



= b 



(2) 



then 



UB = b 



(3) 



assuming that g = lies inside or at the boundaries. 

For n+1 points, equation (1) may be rewritten compactly as 



= 



,n+l 



n+1 



a + b 



(4) 



t _ 1' r 1 n+1, , r i .,12 n+1, _ 
Let a = [a , . . .a J such that [a = 1, [& ,£ ,...£ J a = 0, 

1=1 



and a > for all i. Premultiplying equation (4) with a gives 



59 



a X = b (5) 

and from equation (3) 

UB = a T (6) 

Thus the "price" associated with each point (0 1 ,g 1 ) is 1 . 

The linear programming formulation to solve both problems (a) 
and (b) simultaneously is as follows: 

Find u = [a , . . .a ] to 



Max Z = - [0 ,...0 F ] 



s.t. [g J 



= 



1 2 r 

a + a + . . . + a 1 



and a ,...,a , > = 0. The upper bound UB = - Z. 

The solution will contain (n+1) vectors in the basis. If a 
feasible solution does not exist, artificial variables will be present 
in the solution at a positive level and no upper bound will exist 
because the set of points in hand so far (g ,g , . . .g 1 ) do not surround 
the origin. 



III. 5. Impro vement of the Uppe r Bound 

Let N be the set of the current n+1 points forming the poly- 
tope P, H the hyperplane passing through these points with the slope 



60 



_A, and UB the value of the upper bound. The aim now is to find an 
improvement on UB . 

If a new point is introduced, it follows from the simplex 
method that, since the problem is bounded, a feasible solution and 
hence a new upper bound can be obtained by replacing one of the points 
in N by the new point. The question now is how can the new point be 
found so that the new upper bound (UB ' ) is an improvement on UB, 
i.e. , UB' < UB. 

Once again, the geometry of the problem and the theory of 
linear programming can be used advantageously. Assume for the moment 
that the solution to the linear programming problem obtained from the 
points in N is not degenerate. If a support point is determined for 
slope A_> this new point must lie on or below H. This result is ob- 
vious as a hyperplane with slope A_ cannot be a support hyperplane 

above H. If the coordinates of the new point are (0 ,g ), and 

*-u t r n It new . JA 

the value of on H at g_ = _g_ is , then 

S > n8W 

rt H rt new , _ , 

0=0 implies that no further improvement can be ob- 
tained on the upper bound. Thus if is greater than , it 
follows that UB' will be less than UB if the new point replaces one 
of the points in N by the simplex method. This result can be 
demonstrated as follows: 

Let the points in N be (0 ,g_ ) , . . . (0 ,£ ). The constraints 



/I 2 n+1. _ ,.,. 

(£ ,£ £ ) a = £ (7) 



61 



1,2 n+1 

a + a. . . . + a =1 



(8) 



Combining relations (7) and (8) 



G a 



(9) 



The new point can be represented as a linear combination of the points 
in the basis, that is, 



„, n st n . ~new 

ihe n+1 element of g gives 



a y new 



n+1 
, r new 

i = l y • 
i=i x 



(10) 



(ii) 



Points on H may be expressed by equation (1) 

+ b 



[g T ,l] 







For the (n+1) points in the basis, it follows from equation (4) 



that 



= G' 



+ b 



Premultiplying with y 



T T _, 
new ,, new ^1 
y = y G 



Using equations (10) and (11) 



new _, new 
y = g ,1 



+ b I y\ 
i = l " 



+ b 



(12) 



(13) 



b2 



However, the right hand side is a point on H at g and thus 

T 
new M m H ,, , N 

y = (14) 



In the simplex method, the criterion for an improvement in the 

objective function by replacing a point in the basis by a new one is 

T 
y 0-0 =0-0 >0 (see Hadley (1963)) . 

bmce si - v is greater than zero 

UB' < UB 

If the linear programming solution corresponding to the points 
in N is degenerate, that: is, the point g = lies on one of the bounda- 
ries of P, a difficulty arises. In this case, there is no unique 
hyperplane H and consequently there may be some degrees of freedom 
available in calculating A. There is then no guarantee that the new 
point will bring about an improvement in the upper bound. In such a 
case, several new points may have to be evaluated before an upper 
bound with a lower value is obtained. Degenerate linear programming 
problem solutions occur in Example Problem 3. 

III. 6. Procedure to Inve st igate Globa l Qpt imality 

The basic idea of this approach is that the value for the 
global optimum must lie between or at the boundary of the upper and 
lower bounds. The interval between the bounds is generally decreased 
at every iteration by the improvement on the upper bound and possibly 
an improvement on the lower bound. If at some point during this 



63 



procedure, an optimum is found to possess a value significantly 
greater than the upper bound, an inference can be made that the 
optimum is local and the procedure terminated. On the other hand, 
if a value for the optimum is very close to the lower bound, it may 
be inferred that the optimum is global and the process terminated. 
A problem arises when there is a dual gap in the optimiza- 
tion problem. It is possible that the global optimum may lie above 
the upper bound at some iteration. However, the optimum is not 
rejected as a local optimum if it lies within a certain interval 
above the upper bound, say 2% of the value of the optimum. Previous 
experience indicates that this modification is adequate for the type 
of problems considered in this study. This behavior can be observed 
in Example Problem 2. 

III. 7. Algorithm 



To investigate a primal optimal solution with an objective 
function . 

Step 1 Guess n+1 different A's (A where j=l to n+1). Attempt 
to select these AJ to obtain points g which surround the 
origin. 

Step 2 .For each X , evaluate the dual_h(A J ), the support point 
g , and the objective function 0J . Let p = the total 
number of support points (n+1) in this case) . Find LB 
where LB = Max (h(AJ), j = l to p) . 

Step 3 If < LB + e (c is a small positive quantity, say .02 
0"), stop. is the global optimum. Else go to Step 4. 

Step 4 Formulate an L.P. problem with all the available support 
points and solve for the optimal objective function (-UB) 
and the corresponding basis vector. If a feasible solu- 
tion does not exist, go to Step 7. If > UB + £, stop. 
The given solution is a local optimum. If the L.P. solu- 
tion is degenerate, go to Step 6. 



64 



Step 5 In the L.P. solution, let the vectors in the basis be 

g ,g ,g , . . . g . Set up n simultaneous equations in 
unknowns Ap+I as follows 

, 2 l.T .p+1 rt 2 A 

(g - g ) A = 0-0 



(15) 



,3 l.T ,p+l , A 3 A 
(g - g ) X 1 =0-0 

, 4 l.T .p+1 ,.4 -1 
(g - g ) Y =0-0 

, n+1 l.T ,p+l „,n+l A. 
(g - g ) A F =0 - 

Solve for A' . Find the dual h(A ), the support 
point gP , and the objective function 0P . Set 
p = p + 1. Find LB where 

LB = Max (h(A_J) , j=l,2, . . .p) 

Go to Step 3. 

Step 6 Select all vectors from the basis which are present 

in the L.P. solution in Step 4 at a nonzero level. Let 
these be points g',g~,...g m . Since only m (< n+1) such 
points exist, the AP +J - to be chosen can come anywhere 
from within an n+1 - m dimensional subspace. Either 
initiate (if first time through this step with above 
nonzero basis vectors) or continue a systematic search 
over the subspace to find the next A's to use. A A 
from the appropriate subspace is found by having it 
satisfy 

,2 l.T , p+1 rt 2 .1 
(g - g ) A p =0-0 

(g 3 - gV AP +1 - 3 - 1 (16) 

(g m - gV AP +L = 1 " - 1 

plus any additional n + 1 - m independent specifications. 
Go to Step 3. 

Step 7 Select a set of A to find a new support point g , 
the dual h(AP +1 ), and the objective function P+1 . Set 
p = p + 1. Let LB = Max (h(A.I), where j=l to p) . Go to 
Step 3. (Note this step is actually the same as Step 6 
but with no constraints of form (16) being written, i.e., 
one should search systematically over the entire space 
for the next A's.) 



65 



I I I. 8. Examples 

Example Problem 1 

Figure 11 represents the Example Problem 1, and its subsystems 
are shown in Figure 16. 

h(A) = Min (cost, + Ay) + Min (cost, - Ax) 



where A is the price associated with the variable T . 

8 1 = (x 2 - y l } 
On optimization of the original system the two possible 
solutions are 



Case: 1 T. = 400° T 2 = 900° T 3 = 300° = 286.93 
Case: 2 T = 600° 1' 2 = 700° T = 100° 0* = 189.3 

Apply the algorithm to Case 1,0 = 286.93 
Step 1 A = -4.14 ; A = 4.14 and p = 2 points. 
Step 2 gj = -200 h(xh = -224.6 1 = 476.2 

g^ = +200 h(A^) = -882.0 2 = 

LB = -224.6 

Step 3 0* > LB + e. (Let c = 0.02 0* = 5.74.) 

Step 4 UB = 238.1 
> UB + e 



Hence the solution in Case 1 is a local optimum. 



6b 



F » 10,000 

r - 6ou° 



F - 10,000 

r - 300° 



<> 



F - 10,000 



~ 0" 



T , 



100 s y S 300 



100 s x s J10 



FiR. 1*j. Subsystems for Example Problem 1 



67 



Apply the algorithm to Case 2, 189.3 Steps 1, 2 and 3 are 
the same as in Case 1. 

Step 4 UB = 238. 1 
< UB + E 

Step 5 A 3 = -1.19 

h(X?) =189.3 , gj = , 3 = 189.3 
P = 3 
LB = 189.3 
Step 3 = 189.3 < LB + e 

Hence the solution in Case 2 is a global optimum. 

Example Problem 2 

Figure 17 represents this problem and the subsystems are shown 
in Figure 18. 

h(A) = Min (Area + Ay) + Min (Area„ + A 9 y - A x ; ) 
y ll y i2 ,X 22 

4- Min (Area., - Ax,,,,) 
X 23 



A is the price associated with T and A is the price associated wi 



th 



T . 
2 



8 1 = (x 22 ~ y ll ) and 8 2 = (x 23 " y 12^ 



On optimization 



T x = 181.9° T 2 = 295° t = 218.1° t ? = 286.4° 
t 3 = 395.5° and = 7049. 



68 



T = 100" /^-\ 



-u — - i 



t c J 



FC - 10 tor all the stream 
P 

U, - 120 U„ - 80 IL ■ 40 



Arua refers to exchanger 1 
Area - Q/D il^ , d - Are^ 
Aira: To minimise 4 



>.rc», + Area, 



[•'ig. 17. Exsrunle Problem 



69 



Subsystem 1 

T > 300° 






Subsystem 2 

T - «iO° 



Subsystem 3 

I - 600° 



^A- *_A 



100 £ y s 300 
100° St, s 300° 



100 £ x„„ * 300 
22 

100 £ y £ 400° 

100° s t, £ 400° 



100 £ x £ 400 
200° st, £ 600° 



Fig. 18. Subsystems for Example Problem 2 



70 



Application of the Algorithm: = 7049 



Step 1 


A 1 = 


-10 



2 



10 


A 3 = 


-15 
_-30_ 


Step 2 


1 

S = 


200 
100 


2 
S = 



300 


3 
1 = 


-94.6 
-188 



h(A ) = 500 
1 = 2500 



Step 3 > LB + e 
Step 4 UB = 9001.46 

■k 

< UB + e 



h(A 2 ) = -500 h(A 3 ) = 5787 



Step 5 A 4 = 



-21.67 
-21.67 



= 2500 
LB = 5787 
Let z = 140 



■112.3 
132 



= 12,846 



h(A H ) = 5585 = 5158 p = 4 
LB = 5787 
Step 3 > LB + e 
Step 4 UB = 7173.34 

< UB + e 

13 4 
Step 5 The points in the basis are g , g and g . 



-11.04 
-24.65 



122.87 
7J 



h(A 5 ) = 6640.3 5 = 3533.67 p = 5 



LB = 6640.3 



P = 3 



Step 3 > LB + E 



71 



Step 4 UB = 6928.94 
< UB + e 



3 4 5 
Step 5 The points in the basis are g , g and g . 



-13.34 
-24.76 



111.79 
71 



h(A 6 ) = 6916.26 6 = 3668.11 



P = 6 



LB = 6916.26 

Step 3 < LB + e 

Hence the solution represents a global optimum. 

Example Problem 3 

Figure 2 represents this problem and the subsystems are 
shown in Figure 19. 

y , x and price A are associated with F ? 

y, n , x„, and price A„ are associated with (F. x T , ) 
7 12 24 2 4 4 

y„ 9 , x„ s and price A~ are associated with T ? 

y,.-,, x,., and price A, are associated with T, c 
1326 4 1j 



h(A) = Min (Cost - ^i x ? o + ^o^l 2 + ^-^22^ + 

X 22' y 12' y 22 

Min (Cost„ + Ay -. - A x„, + A, y „) + 



'll> A 24"l3 

Min (Cost- + 0.4F - ^X^) + 



~25 

Min (Cost. + 0.2F TTT - A.x„,) + 
4 III 4 2b 



^26 



and 



72 



? • lO.OL'O 

r - 200° 



Subsystem 2 









y il 


T 


X 24 










F - 10,000 


1 




10,000 

200° 


* 










T - 800 


T 


■ 40U^ 





















2500 <- * 2 2' y i2 s "00 



300 £ y S 500° 



2500 i y 1 ,,x 4 s 7500 
300° s y s 500° 




Sub ay stem t* 

F - 10,000 



300 s x„ c s 500 




Fig. 19. Subsystems for Example Problem 3 



73 



g l = (x 22 ~ y ll } 

8 2 = (x 24 " y 12 } 

g 3 = (x 25 - y 22 ) 

8 4 = (x 26 " y 13 } 



On optimization there are two possible solutions: 



Case: 1 F„ = 2500 F = 7500 T, = T c = 400° 
2 J 4 5 



12 



= 300° 



T 15 = 500° F n = 2500 F m 



4154 



Case ; 



2 Y = 7350 F 3 



2650 T = 500' 



T 15 = 300 c 



F 1I = ° 



F III = 25 °° T 5 = /(23 ° 



T, = 392° 

4 



3538 



Application of the Algorithm to Case 1 4154 



Step 1 A = 







-10 

•10 



A 2 = 




0.000555 









10 





A 4 = 



-0.1 


0.001 


-10 


-10 



A 5 = 




0.001 
-10 
- 5 



P = 5 



Step 2 g J 



5000 

-2,000,000 
■200 
■200 



2 
g = 



-5000 

2,000,000 

200 

200 



2345 
-1,876,000 
■200 

200 



■5000 
2,000,000 







-2345 

1,876,000 
■200 

200 





74 




hCA 1 ) = 2350 1 = 6350 




h(A 2 ) = 230 2 = 1340 




h(A 3 ) = 1917 3 = 3917 




LB = 2350 
h(A 4 ) = 2118 = 4618 






h(A 5 ) = 2188 5 = 5064 


Step 3 


0" > LB + e Let E = 80 


Step 4 


UB = 3845 




> UB + e 




Hence Case 1 corresponds to a local optimum. 


* 

Application of the Algorithm to Case 2 3538 




Steps 1 to 3 same as in Case 1 . 


Step 4 


UB = 3845 




< UB + e 




LP solution is degenerate 


Step 6 


Nonzero points in the basis of LP solution are g and g . 




Let the n+1 points for finding A be g , g , g , g , and g . 




A 6 = 


-0.065 


6 
S = 









0.000225 











-9.81 
-6.58 




-200 
200 






h(A 6 ) = 3338 = 4037 




p = 6 LB = 3338 


Step 3 


> LB + e 


Step 4 


UB = : 


5845 









< ub + e 

LP solution is degenerate 



75 



1 2 
Step 6 Nonzero points in the basis of the LP solution are g and g 

Let the n+1 points for finding A be g , g , g , g , and g . 



A 7 = 


0.128 


0.000709 




-8.68 




-7.72 



7 
8 = 



-2345 

4690 







h(A 7 ) = 3534 



P = 7 



= 4564 
LB = 3534 
Step 3 < LB + e. 

Hence Case 2 represents the global optimum. 

III. 9. Discussion 

In the first and the third examples, the local optima are 
recognized very quickly. The global optimum in each of these 
examples is confirmed albeit with a greater number of steps. In the 
second example, a dual gap (global optimum = 7049, upper bound = 6928.9) 
is present. However, the algorithm proves effective. It may also be 
noticed in the second example that an improvement results in the upper 
bound at every iteration. This is guaranteed in the absence of de- 
generacy. Degenerate solutions occur in the third example and several 
points may have to be evaluated before obtaining an improvement in the 
upper bound. 



CHAPTER IV 

PROCESS SYNTHESIS USING STRUCTURAL PARAMETERS: 

A PROBLEM WITH INEQUALITY CONSTRAINTS 

In this chapter a problem with the structural parameter method 
for process synthesis shall be discussed. The problem is almost 
obvious once stated, but it has not been emphasized in the literature. 
It can seriously affect the expected results. The discussion here does 
not imply that the method is valueless, it simply stresses that care 
must be taken to ensure that the method is really useful for a given 
problem. 

A system of interconnected process subsystems can be modeled 
using structural parameters which are defined by the following 

equations 

N 
x. = y a. .x! i=l,2,3, . . .N 

J = l 



£ a. . £ 1 , y a. . = 1 

1=1 



(1) 



where x. and x! are, respectively, the input and output variables of 
the i-th subsystem, N is the total number of subsystems in the entire 
system and the parameters a., are the structural parameters; that is, 
each is the fraction of the output stream of the j-th subsystem which 
flows into the i-th subsystem. 

By means of structural parameters, a system synthesis problem can 
be transformed into a non-linear programming problem with continuous 



76 



77 



decision variables. This idea has been stated and demonstrated in 

the studies by Umeda, Hirai and Ichikawa (1972), lchikawa and Fan 

(1973), Osakada and Fan (1973), Mishra, Fan and Erickson (1973), and 

Himmelblau (1975). In order to obtain an optimal structure, redundant 

subsystems are usually inserted into a "super" structure in which it 

is hoped that the optimal structure is imbedded. For example, Figure 

20 illustrates one use of the method for synthesizing a heat recovery 

network. The problem data and the stream specifications are shown 

in Table 8, and the problem is described in the tradition of problems 

like 4SP1 in Masso and Rudd (1969) and Lee et al. (1970). A decision 

is to be made on whether or not to use cold stream II to aid in 

cooling hot stream V. The problem is formulated with potentially 

redundant exchangers 2, 4 and 5 and the structural parameter a . 

a 0/ . is the fraction of cold stream II that goes through exchanger 
26 

2. The strategy is, in this formulation, to convert the discrete 
decision of whether or not to use cold stream II to aid in cooling 
hot stream V into a problem of optimizing over the continuous deci- 
sion variable a„, . To avoid a cost discontinuity caused by intro- 
Zb 

ducing an exchanger we must, of course, require the cost of an 
exchanger to approach zero as its area does. 

The notion that we have converted a discrete decision problem 
to a continuous one is invalid here. "It may be observed that hot 
stream V has sufficient heat content to drive cold stream I to its 
final temperature. However, the transfer of heat in exchanger I is 
limited because an inequality constraint in exchanger 2 has to be 
satisfied. This constraint is the approach temperature requirement 



78 



II, Cold 
t - 504.00 
T - 311.11 




Cost - 16,945 S/yx 



VII, Heating Utility 
T - 588.89 



Fig. 20. The Optimal Values of the Parameters 



79 



TABLE 8 
STREAM SPECIFICATIONS AND PROBLEM DATA 



'^J'-rimoN 


. is; is .... 








STREAM 








flow 


g/» 


I 


II 


III 


IV 


V 


vr 


¥11 


1260. 00 


504.00 


1260.00 


Unknown 


2520.00 


Unknown 


Unknown 


Iclet Teapeiature 


i 


311. a 


311.11 


422.22 


311.11 


644. U. 


644.44 


588.3-7 


Outlet Traperatire 


i 


588.89 


533.33 


477.78 


355.56 


477.78 


644.44 


588.89 


Bolllnj Point 


K 


755.56 


755.56 


755.56 


373.33 


755.56 


644. 44 


588.89 


Liquid Heat Capacity 


kJAg s 


A. 19 


4.19 


4.19 


4.19 


4.19 


4.19 


4.19 


Film lk.it Tr.w-.3fer 
Coefficient 


W/m 2 K 


1703.49 


1703.49 


1703.49 


1703.49 


1703 .49 


8517.45 


8517.45 


Coot 


8 


0.00 


0.00 


o.oo 


l.lOxlo" 7 


0.00 


22.05*10~ 7 


5.20xlO~ 7 


Heat of Vaporization 


kJAg 


6?/. 80 


697.60 


697.80 


1786.368 


697.80 


1628.20 


1395.60 


Inlet Vapor fraction 




0.00 


0.00 


0.00 


o.oo 


0.00 


1.00 


1.00 


Dutlet Vapor Fraction 
i — _ 




0.00 


0.00 


0.00 


0.00 


0.00 


o.co 


O.CO 



Approach Temperature = 2.78 K 

Beat exchanger coot equation = (a A ) 

where = arum-vl rata of return = 0.1 

a = 350 

b = 0.6 

Equijwiont down tins = 280 hr./yr. 
The cost of the network is the cost of utility streams +■ tho cost of the exchangers, In l/yr. 
Ill heat archangel's are ausuaed to bo counter-current 



80 



between the inlet hot stream and the outlet cold stream (T > 

533.33° + D, where D is a minimum allowed approach temperature of, say, 

2.78K). By numerical computation one finds that the overall cost can 

be lowered if a is increased (see Figure 21). This increase in a„, 
£D 26 

lowers the flows of cooling utilities IV and VII. Ultimately, when 

the flow of cooling utility IV is zero, an optimal structure is obtained 

with a cost of $16,945/yr. At this point a_, = 0.688. 

26 

Note, however, that when the flow of the cold stream through 
exchanger 2 is zero, and the heat exchanger 2 is totally neglected, 
a network with a cost of $6864/yr. is obtained as shown in Figure 22. 
Thus, a significant discontinuity, which was hoped to have been 
eliminated, has reappeared. It should be emphasized that the data 
for the problem and its formulation are important only in that they 
are plausible and that they help make this point. 

The main observation may be summarized as follows. When certain 
subsystems are rendered redundant during optimization (by the asso- 
ciated flow or a split fraction taking on a value of zero) and if 
inequality constraints are associated with these subsystems which 
force constrained behavior elsewhere, discontinuities very likely 
still exist in the problem. One must still make a discrete decision 
about whether or not to introduce that subsystem. This observation 
means that, if the structural parameters are to be used for a syn- 
thesis problem, and, if inequality constraints are involved, the 
problem has to be formulated very carefully, if indeed it can be, 
to make it a continuous one. 



81 




oi o.2 oT3 o~T 



0.5 0.6 0.7 



Fig. 21. The Discontinuity in Optimization 



, Cold 

- 1260.00 

- 311.11 



III, Cold 
F - 1260. CO 
T - 422.22 



*T - 477.78 




VTI, Heatlr^ Utility 
336.04 
538.89 



Fig. 22. The Optimal Network When Ignoring the Approach 
Temperature Requirement at Null Exchanger 2 of 
Figure 20 



CHAPTER V 
CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER RESEARCH 

The optimization strategy chosen for EROS proves to be efficient. 
It is the author's belief that the number of steps required for con- 
vergence is significantly lowered by rederiving a solution procedure 
every time a constraint violation occurs, and by the use of 'restriction' 
(Geoff rion (1970)). A very useful feature in EROS is its ability to 
find a feasible starting point, if none such is provided by the user. 
The solution yielded by EROS can account for portions of the network 
that already exist and for irregularities likely to occur in the 
process streams. Considering the nature of the problem treated, cost 
per a typical run of EROS seems small. The use of program may also 
be extended in carrying out synthesis via structural parameters but 
the problem must be formulated very carefully as mentioned in Chapter 
IV. In the results shown in Table 7, a feasible starting point is 
obtained in very few iterations. A feasible starting point may also 
be provided by the user as in example 10 of Table 7. As indicated 
in Chapter II, the time required for the rewriting of solution 
procedures after finding a feasible point is relatively small as 
compared to the function evaluations during optimization. Hence, 
the equation solving approach followed is fully justified. 



8 'J 



84 



The practical applicability of a program such as EROS can only 
be considered complete if a capability of handling power generation 
and power intensive units is incorporated. Hence there is a need for 
adding unit models such as turbines and compressors. Hereafter the 
pressure changes will be important and a pressure variable must be 
associated with each segment. These changes might also call for a 
better modeling of a heat-exchanger unit in which the heat-transfer 
coefficients may no longer be evaluated by simplistic relations. 
Currently the cooling curves for streams are approximated by three 
linear regions. A more accurate description could be provided for 
streams by specifying enthalpy versus temperature information for the 
range under consideration. A linear interpolation may be assumed 
between two adjacent points provided in this specification. 

In Chapter III an algorithm is presented which can effectively 
identify if an optimum point is a global optimum or a local one. It 
is practical as given only when applied to systems of interconnected 
subsystems, fortunately, a problem form which is common in engineer- 
ing design. The algorithm develops a lower (dual) bound using the 
ideas of Brosilow and Lasdon (1965) wherein the total system is de- 
composed into subsystems, each of which must be globally optimized. 
If the subsystems are small enough, this requirement is readily satis- 
fied. Frequently the optimal lower (dual) bound is equal to the 
optimal upper (primal) bound, but it also often lies strictly below 
because of nonconvexities in the problem. Previous experience on 
actual engineering problems indicates that this difference, caused by 
a duality gap, is small, i.e., only a small percent of the system cost. 



85 



The algorithm also develops an upper bound to this lower bound. 
Subject to the tolerances required because of the duality gap, an 
optimum more than a few percent above this upper bound is declared 
to be a local optimum, and one at or just above the lower or dual 
bound is considered to be a global optimum. If no discrimination is 
possible at the current step, the algorithm provides a means to 
guarantee improvement of the upper bound for nondegenerate problems, 
the lower bound usually being improved at the same time. The three 
example problems demonstrate that the algorithm is surprisingly 
effective. 

Relating to the algorithm presented in Chapter III, a study is 
recommended to discover the global optimum having discovered that 
the given primal optimum is local. It would be desirable to predict 
a value of A yielding the global optimum by perceiving a trend in 
the value of X from the information gathered. 



APPENDICES 



APPENDIX A 
DESCRIPTION OF THE INPUT AND OUTPUT FEATURES 



DATA INPUT FOR THE PROGRAM EROS 



Card 



Variable Card Format Description 



First Card 



FL 



SC 

Next Card MNODES 
MSEGS 

NS 
LS 

NESTK 



NSTOP 



4F8.1 Scaling factor for Flow Recommended 
FL = Minimum Flow 

Scaling factor for Temperature 
Recommended T = Minimum Temperature 

Scaling factor for split fraction 
Recommended A = 0.2 

Scaling factor for slack variables 
Recommended SC = Minimum Temperature 

11 14 Number of nodes 

Number of stream segments 

Number of streams 

Number of information items in a 
stream = 18 

(length of stream vector) NESTK 
= if a starting point is not 
provided 

= number of constraints specified 
for a starting point 

= if the program is to be executed 
until completion 

= 1 if the program is to be stopped 
after one iteration - to permit 
a check that the input has been 
read correctly 



87 



Card 



Var iable Car d Forma t Descri ption 



NSEGS 



Number of segments (note: not 
streams) on which an upper and a 
lower bound on enthalpy is to be 
specified 



INF 



1NH 



Total number of entries in the 
vector KNF (see the description 
of KNF) . 

= if the option is not used 

Total number of entries in the 
vector KNH (see the description 
of KNH) . 

= if the option is not used 



Optional 

Cards 
(when 
NESTK i 
0) 



1NA 



JWRIT 



KESTK (I) 918 
1=1, NESTK 



Total number of entries in the 
vector KNA (see the description 

of KNA) . 

= if the option is not used 

= if a complete output is to 
be printed 

= 1 if the output during optimi- 
zation is to be suppressed. 

List of codes for the initial 
constraints to be held (see the 
section on how to "code" a con- 
straint) 



Next set 
of Cards 
= 3xNS 
where NS 
is the 
number of 
streams 



STPAR (I, J) 8F8.1 
1=1, 18 
J=l, NS 



Stream Information 
A set of 3 cards is necessary for 
every stream, with up to 8 entries 
per card 



For A ny Stream .J 



1st Card 



STPAR (1,J) 



-2. for undefined flow of a 
process stream 

-1 . for undefined flow of a 
utility (steam or cooling water) 



89 



Card 



Variable Card Format Description 



2nd Card 



STPAR (2, J) 

STPAR (3, J) 

STPAR (4, J) 

STPAR (5, J) 

STPAR (6, J) 

STPAR (7, J) 

STPAR (8, J) 

STPAR (9, J) 

STPAR (10, J) 

STPAR (11, J) 

STPAR (12, J) 

STPAR (13, J) 

STPAR (14, J) 



The undefined flow is identified 
in two different ways because the 
costing of the streams is different 
in each case. 

The cost of a process stream need 
be found once regardless of the 
number of units it passes through. 
A utility stream cost has to be 
determined at every heater or 
cooler that exists in which it is 
a utility. 

Inlet temperature 
= -1. if undefined 

Output temperature 
= -1. if undefined 

= 0. Inlet pressure 
(Not used at present) . 

= 0. Outlet pressure 
(Not used at present). 

Inlet vapor fraction 
= -1. if the inlet temperature 
is undefined 

Outlet vapor fraction 
= -1. if the outlet temperature 
is undefined 

Liquid heat capacity 

Vapor heat capacity 

Heat of vaporization 

Dew point temperature 

(Must be at least 1° greater than 

the bubble point temperature) 

Bubble point temperature 

Liquid-side film heat transfer 
coefficient 

Vapor-side film heat transfer 
coefficient 



90 



Card 



Variable Card Format Description 



3rd Card 



STPAR (15, J) 
STPAR (16, J) 
STPAR (17, J) 



Next 
Card 



STPAR (18, J) 



Al 



4F16.8 



A 2 



Two phase film heat transfer 
coefficient 

Cost coefficient ($/lb.* no. of 
operating hrs. per yr.) 

Lower bound on flow. The lower 
bound must be set to quantity 
greater than zero (say 1% of total 
flow) if the stream is split. 
This helps avoid numerical dif- 
ficulties . 

Upper bound on flows 
(Can be very large) 

Step-size for Kuhn-Tucker test. 
Recommended value = 0.01 (change 
in scaled values for variables) . 

Initial step-size in scaled 
variable values for the complex 
method. Recommended value =0.1 



A3 



A 4 



Minimum allowable approach 
temperature in degrees 

Stopping criterion for the 
optimization (change in objec- 
tive function) . Recommended 
value = 0.00001 



Next 
Card 



CA 



3F16.8 Cost multiplier for an exchanger 
(heater or cooler) x annual rate 
of return. 



EXPO 



Exponent for the calculation of 
the cost of an exchanger, heater 
or cooler. Recommended value = 
0.6. 



CB 



Modification in the cost function 
of an exchanger, heater or cooler. 
Recommended value = 10. (Aids in 
avoiding the creation of local 
minima. ) 

Cost of an exchanger, heater or 
cooler = CA*((Area + CB)**0.6 - 
CB**0.6) 



91 



Card 



Variable 



Card Format Descrii 



Next Set of ARK(I) 
Cards I = 1.MN0DES 



Optional 
Cards 

(when INA / 
0) 



6F10.2 This vector contains the infor- 
mation whether area is specified 
for a given node. 

= 0. if area is not specified 
= value of area, if specified 
6 area entries per card. For 
mixers and splitters, a value 
of zero is specified. 



KNA(I) 1114 Gives information about ex- 
I = 1,INA changers, heaters and coolers 

that are to be related by areas 

Structure of KNA vector 
Entry 1 

Item 1: Number of nodes in 
entry (say, n,) 

Items 2,3,4, . . .rij + 1 

Node numbers in ascending order. 
Repeat for entry 2 etc. New 
entries do not start a new card. 



Optional KNF(l) 
Cards I = l.INF 
(when INF $ 
0) 



Optional KNH(I) 
Cards I = l.INH 
(when INH t 
0) 



1114 



Gives information about un- 
related segments which are to 
have equal flows. No segment 
flows are to be listed here as 
equal if they are already equal 
by the nature of the flowsheet. 

Structure of KNF vector 

Entry 1 

Item 1: Number of stream seg- 
ments in the entry (say, n-i) 

Item 2,3, .. .n + 1 

Stream segment numbers in 
ascending order. 

Repeat for entry 2 etc. New 
entries do not start a new card. 

The description of this vector 
is the same as that of KNF with 
the exception that enthalpies 
are involved instead of flows. 



92 



Card 



Variable Card F ormat Description 



No. of Cards NARRAY(I,J) 1114 

= MNODES J = 1,11 

(No. of 1=1, MNODES 

nodes) 



One card for each node. 
Every card has 11 entries. 



For node I which represents an exchanger heater or cooler 
N ARRAY (1,1) 



NARRAY(I,2) 
NARRAY(I,3) 

NARRAY (I, 4) 
NARRAY(I,5) 

NARRAY(1,6) 
NARRAY (1,7) 

N ARRAY (I, 8) 
NARRAY (I, 9) 

NARRAY (I, 10) 
NARRAY (1,11) 



Stream I.D. for the hot stream 
segments 

hot input segment number 

source node for hot input 
segment 

hot output segment number 

destination node for hot output 
segment 

cold input segment number 

source node for cold input 
segment 

cold output segment number 

destination node for cold 
output segment 

Type of node 

= 1 for a heat exchanger 

= 2 for a heater 

= 3 for a cooler 

Stream I.D. for the cold stream 
segments 



For node I which represents a splitter 
NARRAY (1,1) 



NARRAY (I, 2) 



Stream I.D. for the hot stream 

segment 

= for a cold stream 

hot input segment number 
= -1 for a cold stream 



9 3 



Card Variable Card Format Description 

NARRAY(I,3) source node for hot input segment 

= -1 for a cold stream 

NARRAY(I,4) the first output segment number 

NARRAY(I,5) destination node for the first 

output segment 

NARRAY(I,6) cold input segment number 

= -1 for a hot stream 

NARRAY(I,7) source node for cold input seg- 

ment 
= -1 for a hot stream 

NARRAY(I,8) the second output segment number 

NARRAY(I,9) destination node for the second 

output segment 

NARRAY(I,10) Type of node 

= 4 for a splitter 

NARRAY(I,11) Stream I.D. for the cold stream 

segments 
= for a hot stream 

Note : the order for the output 
segments can be arbitrary. 

For node I representing a mixer 

NARRAY(I,1) Stream I.D. for the hot stream 

segments 
= for a cold stream 



NARRAY(I,2) 
NARRAY(I,3) 

NARRAY(I,4) 

NARRAY(I,5) 

N ARRAY (I, 6) 
NARRAY(I,7) 



the first input segment number 

source node for the first in- 
put segment 

hot output segment number 
= -1 for a cold stream 

destination node for the hot 

output segment 

= -1 for a cold stream 

second input segment number 

source node for the second 
input segment 



94 



Card 



Variable Card For mat Descriptio n 



NARRAY(I,8) 

NARRAY(1,9) 

NARRAY(i,10) 
NARRAY(I,11) 



Optional 

Cards (when ICON (I) 

NSEGS + I = 1,NSEG 

0) 

Optional SEC0N(I,J) 

Cards J = 1,2 

(when NSECS I = 1, NSEGS 
4 0) 



1114 



8F8.1 



cold output segment number 
= -1 tor a hot stream 

destination node for the cold 

output segment 

= -1 for a hot stream 

type of node 
= 5 for a mixer 

Stream I.D. for the cold stream 

segments 

= for a hot stream 

Note : the order for the input 
segments can be arbitrary. 

This vector contains the segment 
numbers for which both the lower 
and upper bounds are to be 
specified for enthalpies. 

SEC0N(1,1) contains the lower 
bounds on enthalpies correspond- 
ing to the segment specified at 
the same 1 location in ICON. 

SEC0N(I,2) contains upper bounds. 

8 entries per card with all the 
lower bounds specified first. 
Start with a new card for upper 
bounds . 



Explanation About the Output From EROS 
The output provides the following information: 

1) A printout of the input. 

2) The progress of optimization (optional) . 

3) Incidence matrix. The rows represent functions and columns 
represent variables. 

4) Information about functions. The first column points to the row 
of the incidence matrix representing the function. The second 
column lists the node associated with the function. The third 
column indicates the type of equation involved. 



95 



TYPE = 1 Heat balance for an exchanger, heater or 

cooler 

TYPE = 2,3 Material balances for a splitter 

TYPE = 4 Material balance for a mixer 

TYPE = 5 Heat balance for a mixer 

TYPE > 10, TYPE < Type is an inequality constraint code. 
(See Table 1.) 

5) Identification of unknown variables. The first column stands for 
the segment involved. If the flow, temperature, or the split 
fraction associated with this segment is an unknown variable, the 
entry indicates the column number in the incidence matrix. Several 
variables may have the same column number if the system discovers 
they must be equal by heat or material balance. If the variable 

is known, a value of is used. 

6) Segment Information - Final values at optimal point. 

7) Precedence order list - Together with 8, indicates how the equa- 
tions were solved during the last optimization step. 

8) Decisions and tears. If a variable is a decision, the correspond- 
ing value for its Function and MLIST columns are -1 and 1 respec- 
tively. If a tear variable is involved, the function it is 
assigned to is indicated in the corresponding column. 

MLIST = 2 represents explicit assignment. 
MLIST = 3 represents explicit assignment. 

9) Objective function - Final value. 

10) Exchanger Areas - FinaL values. 

11) Information about time and number of iterations required. 



96 



O m 

t • 

O <M 



o o 

3 O 



o m 






< 


~ 


u 


CM 


a 


u 


>• 


iy 


^ 


D 




'.•> 


< 


•— 




U- 


a 




a 


2 


u. 


■— 


h- 


s: 


O LU 


U 


-j 



l/l CO I O 



c m 
o in 



o o 

o o 

o 



a 3 

o 



o n 

o 

■o 



X. U. 

o 



00 O 



o o o o 

o o o o 
o o o o 



o o o o O o o 

OOTOOOO 

O O O O 3 D O 

00 U"> O •O 00 o -o 

o o 

.-I J> 

o o o o O o 



o 



<A O 

Q. O 

H- O 

1/0 H 



"3 —I 

O I 



O rH 

I 



o o 

-I o 



97 



O 3 

O 3 



o o 

O O 



o o 
o 



o 


O O 


o 


o 


o 


r-i 


o 





O CO 

o o o 
o o 



m 


o 


a 


o 


«: 


"j 


i ) 


a 




o 




o 




o 




o 




o 




o 




o 




o 




o 




o 




o 




o 



o 


CO 


o 


o 


a. 


o 


o 


X 


o 



O -J- ^> o -t ^ 



vr — < w-* ir> eg cm 



ro m o -< o c 
l 



m m in h n <j 



— • o o ro <n m 



o o 


o 




~> 


o 
















o o 


o 




o 


o 
















■0 IT\ 


u> 




m 


^ 












o 

o 

o 


in -r ^r o o a 


o o o 


:j 


-> 


o 


o 


-1 


.— < 


O 


-r 


a 




M Nfin O n o 


• • • 


• 


• 


* 


• 


• 


< 


o 


u 


3 






O O D 


o 


"> 


~> 


n 


-i 




n 




3 






o z> a 


o 


D 


o 


a 


3 




o 




a 






C3CM 


lf\ 


-) 


!\J 


m 


3 




o 




o 




O H h M O ,T 


o 




a 






O 




.J 




o 






0^ 




--( 






*"* 




—1 




O 

o 


o 
o 




o o a 


o 


j 


c 


o 


D 




a 




tn 


o 


HftJ((l >f N ffv 

>- 





■ft 


ft ft 




«■ 


* 




•ft 


s.u -ft 


» 


ft 


i~ ft 


■ft 


* 


L.J -ft 


■ft 


•«■ 


-1 -ft 




«■ 


G. ft- 


h- 


•>» 


>: ft 


_> 


* 


i_) ft- 


a 


V 


o ft 


►— 


* 


* 


::> 


i< 


i* * 


a 


* 


u * 




•ft 


►-c tt 


■«■ 


•ft 


f- -ft 


■ft 


ft- 


< # 


* 


•ft 


M -ft 




•K- 


i- -ft 




-ft 


s: -ft 




# 


»-« ft- 




ft 


y- # 




■ft 


D. ft 




ft 


O * 




■ft 


ft 




* 


* * 




•ft 


# * 




•ft 


•ft V 



--( O O -H DO 



CNJ .1 



y..< 



u • 

UJ X) 



rr- 


0. 


<1 


L/-J 


> 




z 




_< 




O 


c 


2 


3; 


V 


UJ 


^: 


1- 







•» 


00000000c 





LL 


lu f\i rn — < 


-1 sT 


m -1 -1 


m <\i 




tyi 




■>r 








Q. 






ITi CM 




t— 






000000000 





►— 


> 








LL 


-— 




X 









h- 











•-< 




*— • 


rH -1 3 O O O O O O 
















3C 


< 




ci. 






"0 










2 u 


H 




h- 


OOO^OOO-nO 





<t 










-J 


1/1 




5: 


OO-iOOOr-^oo 





2 


ILI — • — ^ !\j 


01 sf 


<t lA «o 


n ^ 


l~ 


2 






















.4 


O 
O 




Hi 
..J 


O O O ~H O ^1 D O O 





- H 


m 













fNJ 


-1 -< 


_, 


•4 










u. • 


UJ 


C\J 


i-j 






X. 










— 


> 




Q 


000 000 ~* 
















► 2 

2 


(— 




U 


0OO0OO-HO.H 





a. 


-* <^J en 


vT LTA 


1*- -X) 


0^ 


...j 


O 




*-' " 






^L 


^ 






~( 


j.j i-j 


-i 


in 
in 




O <J O "3 r-l — 1 O 3 O 

O -4 O -1 -H ~< O O O 

.— "1 — 1 — ' 







■K 


Li. 








* 



99 



srcooooooocorjooo 



"n! r\i eo im -J- <T -J- in 



-i o o o m m o o o c o 

<i O O O CO O "DOo^O 



h- O O O >-< IN SOOOODOOO 
ZOOO^tMTO O C O T> O O O 

HJ U 2 'O :•! sf sf ■J fl -r -O r\) in i^ »-j 



-3 O0 0000000">00 "5 



U- o o o o o o 



O O O O "3 



OOOOOOOOOOOOOOOO 



I/IOOOOOOOOOOOOO 

^OOOOOOOOOOTOO 



OlOOOOOOOOOOOOOO 



COOOr-O-OQOOOCOOQO^O 



Q-OOOrO'-nOOOOOOOO"3 
SOOQCOmJOOOoOOOO "J 
UJ • • 

i— o o o <~* r\i ooooooooo 

oooo>tMOoooooct)o 
xi 35 ~0 -1 <r sj- o -o o ^ M win m 



O ■-« <NI — < f\l >'*! -J- if in in O O O o o o 



:s O m m .n m toq-joooot 

_!••••• 

O t in J- m -j o O O o o d 

o -n o ii o n m in o o o n 

O N rg s r^J j (*l N 3 O O O 






'N"! -t ifl >ON JOO v O--<rjnsJ- 



JJ 
•fr 1/1 



100 



o o 



o o 
o c 

o o 
o o 



o o 
o o 



o o 
o o 



o o 
o o 

o o 
o o 



o o 


LU 


o o 


Q 


• a 


X 


o c 


c 


o o 




o o 


LU 


o o 


o 


— « -I 


a: 




iU o 




a ? 



a- n o m h m m o 't xjooo 



o 


o 


in 


a 


o 


_l 


o 


— 1 


S\ 


o 


o 


r\J 


3 


— 1 


X) 


-3 


2 


TO 


o 


.0 


LTl 


O 


o 


o 




cr- 


-r 






N 




et- 


^-i 






og 



CM ,-0 -4" -A -O 



■it a: <r co o 



O (M i<1 N 



101 



— >• IT\ ^~i O 

l^- C3 f\j j« 

;u vj- (\> - 1 o 

s: r-i n r-j 



D C 




Q '.'J 




LU O 




u u 




a k 




..r a 




c_ 




~^ 












a ►- uj 




— f- T 




h- 3 — c 




3JH 




_J "3 




1) -O -0 




oo :j 




ii- o 




ll CJ iL 


LLl 


3 


T 


U") < 


— i 


15 -< _j 


H- 


^T .. ) „/ 




— i >- U 


_J 


f- —I U 


■-1 


»-- ^r oo 


►— 


.li :z »-i 


u 


io <r i: 


t~ 



o ■- 



APPENDIX B 
DESCRIPTION OF THE PROGRAM 

An abbreviated block, diagram of the program is illustrated in 
Figure 23. Only the subroutines performing important functions are 
shown. A typical run of EROS takes place in the following manner. 
The main program makes a call to the subroutine EXEC. EXEC reads in 
the input data. Next, the routine START is called to perform an 
initialization wherein all the undefined parameters are set to a value 
of zero. Once this has been accomplished, EXEC performs the task of 
identifying all the unknown variables and numbering them sequentially 
by using the routines FLOCTS, HCNTS , KNOWN and CLEAR. EXEC issues 
directives to find the first solution procedure for the system. With 
the aid of FMATRX, an incidence matrix is created with rows representing 
model equations and columns representing the unknown variables. PRECED 
uses this mati ix to find a solution procedure. At this point the 
control is transferred from EXEC to USPONT, a subroutine coordinating 
the two functions, optimization and creation of all subsequent solution 
procedures . 

Every time a new solution procedure is created, USPONT trans- 
fers control to the optimization routine USOPT. USOPT makes use of 
the complex method represented by OPT, and the initialization for the 
optimization is performed by ICOPT. OPT hands back different sets of 



102 



103 




Fig. 23. Description of the Program EROS 



104 



decision variable values to USOPT. For each set of decision variable 
values, model equations are solved and the objective function found. 
USOPT invokes the use of ANALYZ in order to solve the model equations. 
The convergence of the iteration variables in ANALYZ is aided by the 
Harwell Library subroutine VA05AD. Constraint violations are checked 
for in SUBANZ. If no constraints are violated, USOPT requests for 
the total system cost via the routine COST. If constraint violations 
occur, control is transferred back to USPONT. A decision regarding 
which inequality constraints are to be present as equalities in the 
equation set is made in POINT, and then a new incidence matrix and 
consequently a new solution procedure is found in USPONT. Optimiza- 
tion is started once again, using USOPT. 

When the stopping criteria for optimization are satisfied, 
control is transferred to EXEC from USOPT through USPONT. An output 
is requested by a call to PRINT and the run terminated. 

A brief description of all the subroutines and function 
subprograms used in EROS is presented next. The modules are presented 
in an alphabetical order and unless they are specified explicitly as 
function subprograms, they represent subroutines. 



ACYCLIC 

This routine helps in finding a solution procedure. It 
consists of the acyclic algorithm. 

ADDER 

When a new solution procedure is found, constraints in the 
equation set corresponding to slack variables whose values are zero, 
are held tight. 

AHEI P 

Using the area of a heat-exchanger zone, and the user 
specified constants, cost is calculated in this subroutine. 



105 



ANALYZ 

Groups of equations that can be solved without iterating and 
calculations involving an iteration loop are identified so that the 
model equations may be solved. 

Subroutines used: LOOP 
INIT1 
VA05AD 
SUBANZ 
FTLLX 

AREA 

Areas are calculated for different phase zones and subsequently 
the cost for an exchanger is evaluated. 

Subroutine used: AHELP 

BUBPT 

This is a function subprogram that discovers the bubble-point 
temperature for a stream segment. 

CALFUN 

This routine is called by VA05AD to provide the value of the 
error function. 

Subroutine used: SUBANZ 

CENTRE 

A centroid is found for the points generated by the complex 
method . 

CLEAR 

The unknown variables are numbered sequentially. 

Subroutine used: CLEAR1 

CLEAR1 

An aid to the routine CLEAR. 

C0LOC 

This subroutine is used to find the value of a slack variable 
from the current best point. 

C0N0DE 

The unknown variable involved in a constraint present as an 
equality is calculated by this routine. 

Subroutines used: FILLUP 
ENTHAL 
ENSLAK 
TEMP 
COLOC 



106 



Function Subprograms Used: DEWPT 

CONSTR 

A check is made to see if any constraints are violated. If 
none are violated the slack values for all the constraints are stored. 

Subroutines used: ENTHAL 
TEMP 
ENSLAK 

Function Subprograms used: DEWPT 

BUBPT 

CONVRT 

In this routine, an array that provides information about 
slack variables associated with constraints in the equation set, is 
created. 

COST 

Once all the model equations are solved, the costs associated 
with exchangers and the utility streams are added together in this 
subroutine to provide the value of the objective function. 

Subroutine used: AREA 

DEWPT 

This is a function subprogram that discovers the dew-point 
temperature for a stream segment. 

ENBALH 

During the setting up of an incidence matrix, this routine 
helps in creating a row for the heat-balance of an exchanger. 

ENSLAK 

This routine helps in treating the temperature constraints 
in terms of enthalpies 

Subroutines used: ENTHAL 
TEMP 

ENTHAL 

Used for calculating the enthalpy and vapor fraction of a 
segment given its temperature. For a pure component enthalpy is 
calculated given both temperature and vapor fraction. 

Function Subprograms used: BUBPT 

DEWPT 

EQS 

Helps relate flows, enthalpies or areas for reliability 
analysis. 

EXEC 

This is an executive routine to which all the input data are 
provided. Calls are made to several subroutines for the execution of 



107 


the program. The directive for providing the final output is also 


issued by EXEC. 


Subroutines used: START 


FLOCTS 


HCNTS 


KNOWN 


CLEAR 


FMATRX 


SETUP 


PRECED 


LINK 


ANALYZ 


COST 


USPONT 


PRINT 


FILLUP 


Once an unknown variable is calculated, its value is stored 


for all the segments it occurs in. 


FILLX 


The value of the decision variable is stored for all the 


segments this variable occurs in. 


Subroutine used: TEMP 


FLOCTS 


All the segments having a common flow-rate are identified 



and characterized by a unique number in this routine. 



Subroutine used: 



EQS 



FMATRX 

An incidence matrix is created in which the rows represent 
the heat and material balances and inequality constraints converted 
to equalities. The columns represent the unknown variables. 

Subroutine used: SETUP 

FPASS 

This routine is used when a new solution procedure is to be 
created by replacing one of the constraints present in the equation 
set by the most recently violated constraint. 



HCNTS 

All the segments having a common enthalpy are characterized 
by a unique number in this routine. 



Subroutine used; 



EQS 



108 



HTEX 

Any unknown variable in the heat-balance equation of the 
heat-exchanger is calculated by this subroutine. The error function 
associated with the heat-balance equation is also calculated. 

Subroutine used: FILLUP 
TEMP 

HTTC 

This function subprogram is used to calculate the overall 
heat-transfer coefficient based on the film heat-transfer coefficients 
of the streams involved. 

ICOPT 

Values for decision variables and other parameters involved 
in optimization are initialized. 

Subroutines used: C0L0C 
SORT 

INIT1 

Initializes parameters for the use of VA05AD. 

KNOWN 

In this routine the known values for flows and enthalpies 
are taken into account to help discover what the unknown variables 
are. 

Subroutine used: ENTHAL 

LINK 

An array is created with entries that keep track of the 
segments that have the same flow or enthalpy. 

LOOP 

Identifies all the iteration loops for ANALYZ . 

MB HAD 

A Harwell Library subroutine called by VA05AD. 

MIXER 

The unknown variables involved in the heat and material 
balances of a mixer and the error functions are calculated by this 
subroutine. 

Subroutines used: TEMP 

FILLUP 

OPT 

This routine consists of the modified complex algorithm 
for optimization. 

Subroutine used: CENTRE 



109 



POINT 

In this routine a decision is made regarding what constraints 
are to be incorporated into the equation set to obtain a new solution 
procedure. 

PRECED 

A solution procedure is discovered for the given set of 
equations represented in the matrix form. The decisions and tear 
variables are established. 

Subroutine used: ACYCLIC 

PRINT 

Provides an output. 

RECOG 

An array is created in which information concerning the un- 
known variables (whether they represent flows, enthalpies, or split 
fractions) is stored. 

RELEAS 

Slack variables corresponding to constraints in the equation 
set are identified if they are to be treated as search coordinates. 

SAFE 

The value of the objective function corresponding to the 
best current point is stored. 

SETUP 

Creates an incidence matrix with rows represented only by 
the heat and material balance equations. 

Subroutine used: ENBALH 

SORT 

A routine that creates an array containing information about 
slack variables that are to be used as search coordinates. 

SPLITR 

Any unknown variables in the material balance equations of 
a splitter, and the error functions involved, are calculated by this 
subroutine. 

Subroutine used: FILLUP 

START 

Initializes all the unspecified parameters at the start of 
execution. 

Subroutine used: EQS 



110 



SUBANZ 

For the calculation of an unknown variable, a call is made 
to the unit model routines associated with the variable. After all 
the calculations have been completed, CONSTR is called to check for 
constraint violations. 

Subroutines used: HTEX 

SPLITR 
MIXER 
C0N0DE 
CONSTR 

TEMP 

The temperature and vapor fraction of a segment are calcu- 
lated given its enthalpy. 

USOPT 

Optimization is performed using the strategy followed in 
this study. 

Subroutines called: CONVRT 
ADDER 
RELEAS 
1C0PT 
OPT 

ANALYZ 
COST 

USPONT 

The algorithm for finding a new solution procedure is co- 
ordinated with the optimization strategy. 

Subroutines used: SETUP 
FMATRX 
PRECED 
LINK 
RECOG 
USOPT 
FPASS 
POINT 
PRINT 

VA05AD 

A Harwell Library subroutine that performs a least-square 
fit on a set of non-linear equations. 

Subroutines used: MB11AD 

CALFUN 

In addition, a system subroutine (CPUTIM) on IBM 360/67 ar Carnegie- 
Mellon University was used. This routine provided the elapsed CPU 
time in milliseconds. It was called by the following subroutines: 

EXEC 

ANALYZ 

USPONT 



LITERATURE CITED 

Avriel, M., and A.C. Williams, "An Extension of Geometric Programming 
with Applications in Engineering Optimization," J. of Engr . Math., _5, 
187 (1971). 

Boas, A.H., "Optimization via Linear and Dynamic Programming," Chem. 
Eng., 2P_, 85 (1963). 

Bragin, M.S., "Optimization of Multistage Heat Exchanger Systems," 
Ph.D. dissertation, New York University (1966). 

Brosilow, C. and L. Lasdon, "A Two Level Optimization Technique for 
Recycle Processes," AIChE-IChemE Symp. Ser. No. 4, 75 (1965). 

deBrosse, C.J., and A.W. Westerberg, "A Feasible-Point Algorithm 
for Structured Design Systems in Chemical Engineering," AIChE J., 
19, 251 (1973). 

Director, S.W., Hachtel, CD. and L.M. Vidigal, "Computer Efficient 
Yield Estimation Procedure Based on Simplicial Approximation," IEEE 
Trans. Circuits and Systems, to be published (1978). 

Fan, L.T., and C.S. Wang, "The Discrete Maximum Principle," Wiley, 
N.Y. (1964). 

Geoffrion, A.M., "Elements of Large-Scale Mathematical Programming, 1" 
the Rand Corporation, R-481-PR, Santa Monica, Calif., November (1969). 

Grossman, I.E., and R.W.H. Sargent, a personal communication (1977). 

Guthrie, K.M. , "Capital Cost Estimating," Chem. Eng., pp 114-142, 
March 24 (1969) . 

Hadley, G., "Linear Programming, " Addison-Wesley, Reading, Massachusetts, 
pp 158-162 (1963) . 

Henley, E.J., and R.A. Williams, "Graph Theory in Modern Engineering," 
Academic Press, N.Y. (1973). 

Himmelblau, D.M., "Optimal Design via Structural Parameters and Non- 
linear Programming," U.S. -Japan Joint Seminar on Application of 
Process Systems Engineering to Chemical Technology Assessment, Kyoto, 
Japan, June 23-27, 1975. 



Ill 



112 



Hutchison, H.P., and C.F. Shewchuk, "Computational Method for Multiple 
Distillation Towers," Trans. Instu. Chem. Eng. , 52 (4) , p 325 (1974). 

Hwa, C.S., "Mathematical Formulation and Optimization of Heat Exchanger 
Networks Using Separable Programming," Symp . Series No. 4, pp 101-106, 
AIChE-IChemE Joint Meeting, London (1965). 

Ichikawa, A., and L.T. Fan, "Optimal Synthesis of Process Systems. 
Necessary Condition for Optimal System and its Use in Synthesis of 
Systems," Chem. Eng. Sci., _28, 357 (1973). 

Kubicek, M. , V. Hlavacek and F. Prochaska, "Global Modular Newton- 
Raphson Technique for Simulation of an Interconnected Plant Applied 
to Complex Rectification Columns," Chem. Engr. Sci., 31, 277 (1976). 

Kuhn, H.W., and A.W. Tucker, "Non-Linear Programming," in Proc. 2d 
Berkeley Symp. on Math. Statistics Prob., 481, Univ. of Calif. Press, 
Berkeley (1951). 

Lasdon, L.A., "Optimization Theory for Large Systems," Macmillan Co., 
New York (1970) . 

Lee, W. , J.H. Christensen and D.F. Rudd, "Design Variable Selection 
to Simplify Process Calculations," AIChE J., J_2, 1104 (1966). 

Lee, K.F., A.H. Masso and D.F. Rudd, "Branch and Bound Synthesis of 
Integrated Process Designs," Ind. Eng. Chem. Fundamentals, % 48 
(1970). 

Leigh, M.J., CD. Jackson and R.W.H. Sargent, "SPEED-UP — A Computer- 
Based System for the Design of Chemical Processes," paper presented 
at CAD-74, Imperial College, London, England, Sept. 24-27, 1974. 

Masso, A.M., and D.F. Rudd, "The Synthesis of System Designs - II 
Heuristic Structuring," AIChE J., L5_, 10 (1969). 

McGalliard, R.L., and A.W. Westerberg, "Structural Sensitivity Analysis 
in Design Synthesis," Chem. Eng. J., 4 (4), 127 (1972). 

Mishra, P.N., L.T. Fan and L.E. Erickson, "Biological Wastewater 
Treatment System Design," Can. J. Chem. Eng., 51_, 694 (1973). 

Nishida, N., Y.A. Liu and L. Lapidus, "Studies in Chemical Process 
Design and Synthesis: III. A Simple and Practical Approach to the 
Optimal Synthesis of Heat Exchanger Networks," AIChE J., 23, 77 (1977). 

Osakada, K. , and L.T. Fan, "Synthesis of an Optimal Large-Scale Inter- 
connected System by Structural Parameter Method Coupled with Multi- 
level Technique," Can. J. Chem. Eng., 51, 94 (1973). 

Perry, R.H., and C.H. Chilton, Chemical Engineering Handbook, 5th 
Edition, McGraw-Hill, New York, pp 10-39 to 10-42 (1973). 



113 



Takamatsu, T., I. Hashimoto, H. Nishitani and S. Tomita, "The 
Optimal Design of Large Complex Chemical Processes - Develop- 
ment of Max-Sensitive Method," Chem. Engng. Sci., 31, 705 
(1976). 

Takamatsu, T. , I. Hashimoto and H. Ohno, "Optimal Design of a 
Langle Complex System from the View Point of Sensitivity 
Analysis," Ind. Engng. Chem. Proc. Des. and Dev., 9, 368 (1970) 

Umeda, T. , A. Hirai and A. Ichikawa, "Synthesis of Optimal 
Processing System by an Integrated Approach," Chem. Eng . Sci., 
27, 795 (1972). 

Umeda, T., and A. Ichikawa, "A Modified Complex Method for 
Optimization," Ind. Engng. Chem. Proc. Des. and Dev., 1_0, 
229 (1971). 

Westbrook, G.T., "Use This Method to Size Each Stage for Best 
Operation," Hydrocarbon Processing, 40, 201 (1961). 

Westerberg, A.W. , and C.J. deBrosse, "An Optimization Algorithm 
for Structured Design Systems," AIChE J., 1_9, 335 (1973). 



BIOGRAPHICAL SKETCH 

Jigar Shah, the eldest son of Hernlata and Vinubhai Shah, was born 
on December 27, 1951, in Ahmedabad, India. He graduated from St. Peter's 
Boys' High School in Panchgani in December 1967. From July 1968 to May 
1973 he studied chemical engineering at Indian Institute of Technology, 
Madras. After the receipt of the B.Tech. degree, he was accepted as a 
graduate student at the University of Florida. He completed his M.S. in 
chemical engineering in December 1974 and proceeded to work for a doctoral 
degree. While working on doctoral research he accompanied his advisor, 
Dr. A. Westerberg, to Computer-Aided Design Centre, Cambridge, U.K., 
from January 1975 to June 1975, and to Carnegie-Mellon University, 
Pittsburgh, from September 1976 to January 1978. Jigar Shah is cur- 
rently a member of the American Institute of Chemical Engineers. 



114 



I certify that I have read this study and that in my opinion it 
conforms to acceptable standards of scholarly presentation and is fully 
adequate, in scope and quality, as a dissertation for the degree of 
Doctor of Philosophy. 



;ki- ■(■ 



VI t 



(/ 



/__ty 



Arthur W. Westerberg, Chairman 
Professor of Chemical Engineering 



1 certify that I have read this study and that in my opinion it 
conforms to acceptable standards of scholarly presentation and is fully 
adequate, in scope and quality, as a dissertation for the degree of 
Doctor of Philosophy. 



An f <- i-.-v..ct : 1 



John P. O'Connell 

Professor of Chemical Engineering 



I certify that I have read this study and that in my opinion it 
conforms to acceptable standards of scholarly presentation and is fully 
adequate, in scope and quality, as a dissertation for the degree of 
Doctor of Philosophy. 




Hugb/T). Rati iff 
Professor of Industrial ^and Systems 
Engineering 



This dissertation was submitted to the Craduate Faculty of the College 
of Engineering and to the Graduate Council, and was accepted as partial 
fulfillment of the requirements for the degree of Doctor of Philosophy. 



March 1978 



lULta 









Dean, College of Engineering 



Dean, Graduate School 



UNIVERSITY OF FLORIDA 

3 1262 08666 281 3 



I