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

Figure 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. , , ,nst

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 g1 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.)

VR = 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

VT = 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)

Greek 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) .

0 = objective function; may be superscripted as 0 to represent its value at point i. May also be used as 0 to denote a vector of 01!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

Mathematical 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 0 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 0 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 Specification

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 - 0 4

For Stream III Tja - 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 C1 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)

F3 = (1 - a)F1 (2)

h2 = hx (3)

h3 = h1 (4)

F4 « F2 (5)

(6) h2F2 + hllFll = h4F4 + h12F12 (?)

(8)

(9) h3F3 + h14F14 = h5F5 + h15F15 (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 Fg - Fy (13)

F13 - F12 (14)

h7F7 +h12Fl2 = h8F8+hl3F13 (15)

(16)

(17) h9F9+h15F15 = h10F10 + h16F16 (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 . T2 > T3 + 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 = fMfP(aAm)

Where f is a materials factor, fD 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

T5 Tnrc T« Ti

Scg 3 F3.T3.H3

Cold Str

rmn " Dew polnt temrcrat:ure. hoC »trem

r - Bubble point teinperature , hot strca

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

T5

=

l^T'1"

urc of

the

cold

stre.-E

i at a

locaClo

3 where

the

hot a

tream

teraper-

aCure 1

TDFM

T6

-

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

T8

-

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

VT2

(NODE

• 10) + 3 * 33

T^T3

(NODE

* 10) + U = 34

F1IR = lower bound LUi for How F

F12FLB-

(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 r3 3 LB

(NODE

* 10) + 7 = 37

F3DB = 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

WT6 + 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\

K3iH3LB

- (SEG « + 1)

1000 + NODE » 10 ^-3031

H^rm = uPPer 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 . Deriving ^-h£ Solution Procedure and_ Solving Model Equations

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

F7 = F7LB + °55

T3 - T5 + a33

(19)

(20)

where F is the lower bound on F7. /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

T3 - f(h3) (21)

T5 = f(h5) (22)

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

h3 = f(h5,o33) . (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 Fi'Fii'Fi4'F7LB

Enthalpies ^V'? 'h8'Vh10'hll 'h13>h14'h16

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 Point

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\

xf1

-*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. h5

4. Guess Fp(=F. )

5. a

6. F (=F '

7. h,

8. F,

'• *-2<=y

10. H

15

ii. f9(=f10:

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 feasible 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: (Ec - EQ).

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 0 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 0 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.)

EQ = [(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 0 c

2) Solution procedure derived

4) Procedure 'A' executed

Constraint Violation = 33.

26

E = [En, 33] c U

Solution procedure derived

Procedure 'A' executed

Constraint Violation = 25

Degree of freedom = 1

Ec = [EQ5 33, 25]

Solution procedure derived

Procedure 'A' executed

Constraint Violation = 22

Degree of freedom = 0

H = [33, 25, 22], V = [33, 25], V2 = [22, 33], V3 = [22, 25;

Ec = [EQ, V2] - [EQ, 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 0 H = '0'

Degree of freedom = 0

H = [22, 33, 55], V = [22, 33], V2 = [55, 22], V3 = [55, 33 j

Ec = [EQ, V2] = [EQ, 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 = [Vm,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 g1 is detected as being violated. The algorithm will respond by changing the solution procedure so that the slack variable 0 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 Algorithm

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 VT = V , VR = '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 Vc, 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 VR 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 VT. Set V = [VD)Vn>], 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 VR

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 VR = '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, (NT7 + 1 constraints) are all satisfied in a c 1 Vc

subspace of dimension NTI . 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 VR 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 VT = VTold, 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) VR = [22]

3 and 4) Optimization with the complex method.

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

8) Replace constraint 22 with constraint 35.

9) VR = [35], VT = [55]. Ec = [EQ,VR,VT].

3) Optimization with the complex method.

4) Optimization complete.

5 and 6) Vm 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„ = Fc = 2650, T, „ = 500°, T, c = 300°, 2 4 3 5 12 15

F7 = 0, Fn = 2500, T. = 392°, T = 423° and 0 = 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 .

Reliability 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 0 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 0

-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 0 is given by

0 = ch1fh1 + ch;fh; + ch-'fh^ + cuc"Fuc" +

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 0 . 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[(A2,)U'6 + (A2 - A2,)°-6 + (A2„ - A2)°-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 0

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° 150C

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

Cost Coefficient

0.10 0.015 0.015 0.008

Unspecified.

Let U. represent the heat transfer coefficient for exchanger I

ui ; V = 7CJ0

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

U3 = 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 0 is defined as

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

13

+ 35 I (A. + 10)°-6 - ID0"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

0 = $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 - 100,6] (1)

is used in preference to

cost = 35 A0"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 0 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 0 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

0 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 . Statement of tlie_ Problem

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 Lower Bound

To permit decomposition of a system structure, the primal or

overall system optimization problem may be written as

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

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

o1 ± 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 Area0'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 .

0

100

ibo

200 T3

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 = 0 - I j, A.[x. - t.(q.)]

j=l ie0(j)

which may be rewritten as

n n

L= J ff.(q.) + /V *4t4(q,) ~ I AixJ = LfMi'h)

I [f.CqJ + I \h(^ ~ I Aixi] = .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 0 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 0 at the intersection of this polytope with the axis £ = 0 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 0 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 = 0 in its interior or at its boundaries, a hyperplane

intersecting the axis g = 0 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 = 0 on or below the point

the polytope intersects this axis. Consequently the dual cannot exceed

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

Ill .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, 0 can be represented by the linear relationship

T

0 = £ a + b

(1)

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

= b

(2)

then

UB = b

(3)

assuming that g = 0 lies inside or at the boundaries.

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

0 =

,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

aX0 = b (5)

and from equation (3)

UB = aT0 (6)

Thus the "price" associated with each point (01,g1) is 01.

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

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

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

s.t. [gJ

= 0

1 2 r

a + a + . . . + a1

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 , . . .g1) do not surround the origin.

III. 5. Improvement of the Upper 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 0 on H at g_ = _g_ is 0 , then

0S > 0n8W

rtH rtnew , _ ,

0=0 implies that no further improvement can be ob- tained on the upper bound. Thus if 0 is greater than 0 , 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 ynew

n+1 , r new

i = l y • i=i x

(10)

(ii)

Points on H may be expressed by equation (1)

+ b

[gT,l]

0

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 0 = y G

Using equations (10) and (11)

new _, new y 0 = 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 mH ,, , N

y 0 = 0 (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 = 0 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 Investigate Global Qptimality

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 0 .

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(AJ), 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 0 < LB + e (c is a small positive quantity, say .02 0"), stop. 0 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 0 > 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 rt2 A

(g - g ) A = 0-0

(15)

,3 l.T ,p+l ,A3 A (g - g ) X1 =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 ) AF =0 - 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~,...gm. 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 rt2 .1 (g - g ) Ap =0-0

(g3 - gV AP+1 - 03 - 01 (16)

(gm - gV AP+L = 01" - 01

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 0P+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 .

81 = (x2 - yl} On optimization of the original system the two possible solutions are

Case: 1 T. = 400° T2 = 900° T3 = 300° 0 = 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 01 = 476.2

g^ = +200 h(A^) = -882.0 02 = 0

LB = -224.6

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

Step 4 UB = 238.1 0 > 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, 0 189.3 Steps 1, 2 and 3 are the same as in Case 1.

Step 4 UB = 238. 1 0 < UB + E

Step 5 A3 = -1.19

h(X?) =189.3 , gj = 0 , 03 = 189.3 P = 3 LB = 189.3 Step 3 0 = 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„ + A9y - A x ;) yll yi2,X22

4- Min (Area., - Ax,,,,) X23

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

th

T . 2

81 = (x22 ~ yll) and 82 = (x23 " y12^

On optimization

Tx = 181.9° T2 = 295° t = 218.1° t? = 286.4° t3 = 395.5° and 0 = 7049.

68

T = 100" /^-\

-u — - i

tcJ

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: 0 = 7049

Step 1

A1 =

-10 0

2

0 10

A3 =

-15 _-30_

Step 2

1

S =

200 100

2 S =

0 300

3 1 =

-94.6 -188

h(A ) = 500 01 = 2500

Step 3 0 > LB + e Step 4 UB = 9001.46

■k

0 < UB + e

h(A2) = -500 h(A3) = 5787

Step 5 A4 =

-21.67 -21.67

0 = 2500 LB = 5787 Let z = 140

■112.3 132

= 12,846

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

0 < UB + e

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

-11.04 -24.65

122.87 7J

h(A5) = 6640.3 05= 3533.67 p = 5

LB = 6640.3

P = 3

Step 3 0 > LB + E

71

Step 4 UB = 6928.94 0 < 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(A6) = 6916.26 06 = 3668.11

P = 6

LB = 6916.26

Step 3 0 < 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 , ) 712 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 - ^ix?o + ^o^l 2 + ^-^22^ +

X22'y12'y22

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

'll>A24"l3

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

~25

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

^26

and

72

? • lO.OL'O

r - 200°

Subsystem 2

yil

T

X24

F - 10,000

1

10,000

200°

*

T - 800

T

■ 40U^

2500 <- *22'yi2 s "00

300 £ y S 500°

2500 i y1,,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

gl = (x22 ~ yll}

82 = (x24 " y12}

g3 = (x25 - y22)

84 = (x26 " y13}

On optimization there are two possible solutions:

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

12

= 300°

T15 = 500° Fn = 2500 Fm

4154

Case ;

2 Y = 7350 F3

2650 T = 500'

T15 = 300c

F1I = °

FIII = 25°° T5 = /(23°

T, = 392°

4

3538

Application of the Algorithm to Case 1 0 4154

Step 1 A =

0

0

-10

•10

A2 =

0 0.000555 0 0

0

0

10

0

A4 =

-0.1

0.001

-10

-10

A5 =

0 0.001 -10 - 5

P = 5

Step 2 gJ

5000

-2,000,000 ■200 ■200

2 g =

-5000

2,000,000

200

200

2345 -1,876,000 ■200

200

■5000 2,000,000

0

0

-2345

1,876,000 ■200

200

74

hCA1) = 2350 01 = 6350

h(A2) = 230 02 = 1340

h(A3) = 1917 03 = 3917

LB = 2350 h(A4) = 2118 0 = 4618

h(A5) = 2188 05 = 5064

Step 3

0" > LB + e Let E = 80

Step 4

UB = 3845

0 > UB + e

Hence Case 1 corresponds to a local optimum.

*

Application of the Algorithm to Case 2 0 3538

Steps 1 to 3 same as in Case 1 .

Step 4

UB = 3845

0 < 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 .

A6 =

-0.065

6 S =

0

0.000225

0

-9.81 -6.58

-200 200

h(A6) = 3338 0 = 4037

p = 6 LB = 3338

Step 3

0 > LB + e

Step 4

UB = :

5845

0 < 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 .

A7 =

0.128

0.000709

-8.68

-7.72

7 8 =

-2345

4690

0

0

h(A7) = 3534

P = 7

0 = 4564 LB = 3534 Step 3 0 < 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 .

a0/. 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/m2 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 = 0 if a starting point is not provided

= number of constraints specified for a starting point

= 0 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

Variable Card Format Description

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) .

= 0 if the option is not used

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

= 0 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) .

= 0 if the option is not used

= 0 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 Any 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 Format 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

= 0 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 = 0 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 = 0 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 Format Description

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

= 0 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 < 0 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 0 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

0

LL

lu f\i rn — <

-1 sT

m -1 -1

m <\i

tyi

■>r

Q.

ITi CM

t—

000000000

0

►—

>

LL

-—

X

0

h-

0

•-<

*— •

rH -1 3 O O O O O O

0

0

3C

<

ci.

"0

2 u

H

h-

OOO^OOO-nO

0

<t

0 -J

1/1

5:

OO-iOOOr-^oo

0

2

ILI — • — ^ !\j

01 sf

<t lA «o

n ^

l~

2

0

0

.4

O O

Hi ..J

O O O ~H O ^1 D O O

0

-H

m

0

fNJ

0 0 -1 0 0 -< 0 0 0

_,

•4

u. •

UJ

C\J

i-j

X.

— 0

>

Q

000 000 0 ~* 0

0

0

► 2

2

(—

U

0OO0OO-HO.H

0

a.

0 -* <^J en

vT LTA

0 1*- -X)

0^ 0

...j 0

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

.— 0 "1 0 — 1 — ' 0 0 0

0 0 0

■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 0 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 JOOvO--<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 <-0i-.-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