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
Live Music Archive
Librivox Free Audio
Metropolitan Museum
Cleveland Museum of Art
Internet Arcade
Console Living Room
Books to Borrow
Open Library
TV News
Understanding 9/11