simulated annealing technique to design minimum cost exchanger

19 Pages·2011·0.52 MB·English

Available on line at
Association of the Chemical Engineers of Serbia AChE

www.ache.org.rs/CICEQ
Chemical Industry & Chemical Engineering Quarterly 17 (4) 409−427 (2011) CI&CEQ

NADEEM M. KHALFE SIMULATED ANNEALING TECHNIQUE TO
SANDIP KUMAR LAHIRI DESIGN MINIMUM COST EXCHANGER
SHIV KUMAR WADHWA
Jubail Industrial City, Saudi Arabia Owing to the wide utilization of heat exchangers in industrial processes, their
cost minimization is an important target for both designers and users. Tra-
SCIENTIFIC PAPER ditional design approaches are based on iterative procedures which gradually
change the design and geometric parameters to satisfy a given heat duty and
UDC 621.78:519.876.5
constraints. Although well proven, this kind of approach is time consuming and
may not lead to cost effective design as no cost criteria are explicitly accounted
DOI 10.2298/CICEQ110204027K
for. The present study explores the use of non-traditional optimization tech-
nique called simulated annealing (SA), for design optimization of shell and
tube heat exchangers from an economic point of view. The optimization pro-
cedure involves the selection of the major geometric parameters such as tube
diameters, tube length, baffle spacing, number of tube passes, tube layout,
type of head, baffle cut, etc., and minimization of total annual cost is consi-
dered as the design target. The presented simulated annealing technique is
simple in concept, few in parameters and easy for implementations. Further-
more, the SA algorithm explores the good quality solutions quickly, giving the
designer more degrees of freedom in the final choice with respect to traditional
methods. The methodology takes into account the geometric and operational
constraints typically recommended by design codes. Three different case stu-
dies are presented to demonstrate the effectiveness and accuracy of proposed
algorithm. The SA approach is able to reduce the total cost of the heat ex-
changer compared to cost obtained by the previously reported GA approach.
Keywords: simulated annealing; heat exchanger design; optimization;
heat transfer; mathematical modeling.

Chemical process industries are currently facing nities. Designers of STHE normally keep a design
the economic squeeze. Global demand for oils and margin to accommodate any uncertainties in design
chemical products is low, revenue has fallen quickly calculations and to ensure that the heat exchanger
and economic downturn has reduced the availability deliver its services in actual shop floor. The need of
of financing for working capital and investment. This the hour is to reduce the investment cost of STHE by
cut throat competition and shrinking profit margin trimming down the fat in design through more efficient
forced the process industries to introspect critically design strategies. The classical approach to STHE
the new investment decision. Shell and tube heat ex- design involves a significant amount of trial-and-error
changers (STHE) are the most common type of ther- because an acceptable design needs to satisfy a
mal equipment employed in chemical process Indus- number of constraints (e.g., fouling allowance and al-
tries and contribute a major portion of capital invest- lowable pressure drops). Computer software mar-
ment in new projects. Because of their sheer large keted by companies such as Heat Transfer Research,
numbers in any chemical plants, small improvements Inc. (HTRI), and Heat Transfer and Fluid Flow Service
in STHE design strategies offer big saving opportu- (HTFS) are used extensively in the thermal design
and rating of heat exchangers. These packages incor-
porate various design options for the heat exchangers
Correspondening author: N.M. Khalfe, , P.O. Box 10085, 31961
Jubail Industrial City, Saudi Arabia. including the variations in the tube diameter, tube
E-mail: [e-mail] pitch, shell type, number of tube passes, baffle spa-
Paper received: 4 February, 2011
cing, baffle cut, etc. Typically, a designer chooses va-
Paper revised: 1 July, 2011
Paper accepted: 3 July, 2011
409
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

rious geometrical parameters such as tube length, straint for achieving optimal design parameters. All of
shell diameter and the baffle spacing based on expe- the above researchers concluded that these algo-
rience to arrive at a possible design. If the design rithms result in considerable savings in computational
does not satisfy the constraints, a new set of geomet- time compared to an exhaustive search, and have an
rical parameters must be chosen to check if there is advantage over other methods in obtaining multiple
any possibility of reducing the heat transfer area while solutions of the same quality, thus providing more fle-
satisfying the constraints. Although well proven, this xibility to the designer. Fesanghary et al. [8] used
kind of approach is time consuming and may not lead global sensitivity analysis to identify the most influen-
to cost effective design as no cost criteria are expli- tial geometrical parameters (like tube diameter, shell
citly accounted for. Since several discrete combina- diameter, baffle spacing, etc.) that affect total cost of
tions of the design configurations are possible, the STHEs in order to reduce the size of optimization pro-
designer needs an efficient strategy to quickly locate blem and carried out optimization by applying harmo-
the design configuration having the minimum heat ex- nic search. Recently, Patel and Rao [9] have applied
changer cost. Thus the optimal design of heat ex- particle swarm optimization technique to design lowest
changer can be posed as a large scale, discrete, cost heat exchangers.
combinatorial optimization problem [1]. In view of the encouraging results found out by
In literature, attempts to automate and optimize the above researchers, an attempt has been made in
the heat exchanger design process have been pro- the present study to apply a strategy called simulated
posed for a long time and the subject is still evolving. annealing (SA), to the optimal heat exchanger design
The suggested approaches mainly vary in the choice problem. Simulated annealing (SA), a recent optimi-
of the objective function, in the number and kind of zation technique, is an exceptionally simple evolution
sizing parameters utilized and in the numerical opti- strategy that is significantly faster and robust at nu-
mization method employed. In commercial software, merical optimization and is more likely to find a func-
heat exchanger cost function has been recently incor- tion’s true global optimum. Simulated annealing (SA)
porated and cost minimization is performed by apply- algorithms that are members of the stochastic optimi-
ing mainly gradient based methods. Depending upon zation formalisms have been used with a great suc-
the degree of non-linearity and initial guess, most of cess in solving problems involving very large search
the traditional optimization techniques based on gra- spaces. This optimization technique resembles the
dient methods have the possibility of getting trapped cooling process of molten metals through annealing.
at local optimum. Hence, these traditional optimiza- The cooling phenomenon is simulated by controlling a
tion techniques do not ensure global optimum and also temperature like parameter introduced with the con-
have limited applications. In the recent past, some ex- cept of the Boltzmann probability distribution. Accord-
pert systems based on natural phenomena (evolu- ing to the Boltzmann probability distribution, a system
tionary computation) such as simulated annealing [2] in thermal equilibrium at a temperature T has its ener-
and genetic algorithms [3,4] have been developed to gy distributed probabilistically according to P(E) =
overcome this problem. = exp(–E/kT), where k is the Boltzmann constant. This
Chaudhuri et al. [1] used simulated annealing for expression suggests that a system at high tempera-
the optimal design of heat exchangers and developed ture has almost uniform probability of being in a high
a command procedure, to run the HTRI design pro- energy state. Therefore, by controlling the tempera-
gram coupled to the annealing algorithm, iteratively. ture T and assuming that the search process follows
They have compared the results of the SA program the Boltzmann probability distribution, the conver-
with a base case design and concluded that signi- gence of an algorithm can be controlled. There are
ficant savings in the heat transfer area and hence the numerous papers in literature discussing application
STHE cost can be obtained using SA. Manish et al. of SA in various problems. In a comprehensive study
[5] used a genetic algorithm framework to solve this of SA, Johnson et al. [10-12] discuss the performance
optimal problem of heat exchanger design along with of SA on four problems: the travelling salesman prob-
SA and compared the performance of SA and GAs in lem (TSP), graph partitioning problem (GPP), graph
solving this problem. They also presented GA strate- coloring problem and number partitioning problem. In
gies to improve the performance of the optimization general, the performance of SA was mixed – in some
framework. Recently, Caputo et al. [6] also used GA problems, it outperformed the best known heuristics
based optimization of heat exchanger design. Selbas for these problems, and, in other cases, specialized
et al. [7] used genetic algorithm for optimal design of heuristics performed better.
STHEs, in which pressure drop was applied as a con-
410
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

The main objective of this study is to explore the b. Evaluation of the capital investment, ope-
effectiveness of Simulated Annealing technique in the rating cost and the objective function.
design optimization of STHEs from economic point of c. Utilization of the optimization algorithm to se-
view. Ability of the SA based technique is demon- lect a new set of values for the design variables.
strated using different case studies and parametric d. Iterations of the previous steps until a mini-
analysis. mum of the objective function is found.
The second objective of the present study is to The entire process is schematized in Figure 1.
design the optimum heat exchanger, which would
User input. Following parameters are required as a
comply with tubular exchangers manufactures asso-
user defined input to calculate the heat exchanger
ciation (TEMA) standards and obey the industrial re-
area:
quirement of geometric, velocity and pressure drop
constraints. Most of the optimum solutions (say tube 1. Mass flow rate and inlet/outlet temperatures
diameter, tube length, baffle spacing, shell diameter, of shell side and tube side fluids.
etc.) found in the literature are not available as per 2. Thermophysical properties of both fluids, e.g.,
TEMA standard sizes and thus makes the fabrication density, viscosity, heat capacity, thermal conductivity.
of heat exchanger costly due to nonstandard sizes. 3. Fouling resistances, Rfoul, shell and Rfoul, tube.
Also, some practical design rules such as geometric Input optimization variables. The optimization
constraints, velocity and pressure drop constraints variables, with values assigned iteratively by the opti-
are usually ignored in the algorithms present in the mization techniques (e.g. SA) are given in Table 1.
literature, which restrains the effective application of Calculation sequence
design solutions found. In spite of the new algorithmic
1. Assume overall heat transfer coefficient based
developments applied to heat exchanger design in
on type of shell and tube side fluid.
literature, the complexity of the task allows some cri-
2. Based on the actual values of the design
ticism of the effectiveness of optimization procedures
specifications and on the current values of the opti-
for real industrial problems [13]. In this context of the
mization variables, the exchanger design routine de-
development of new design algorithms, this paper
termines the values of the shell side and tube side
presents an optimization procedure integrated with
heat transfer coefficients, the overall exchanger area,
practical design guidelines, aiming to provide a fea-
the number of tubes, the shell diameter and tube side
sible alternative in an engineering point of view.
and shell side flow velocities, thus defining all cons-
With the application of SA algorithm, it is found
tructive details of the exchanger satisfying the as-
in the present study that multiple heat exchanger con-
signed thermal duty specifications.
figurations are possible with practically same cost or
3. Overall actual heat transfer coefficient is then
with little cost difference. The lowest cost exchangers
calculated and compared with the assumed values of
are not always performing best in actual shop floor.
U.
Maintainability, ease of cleaning of tubes and shells,
4. If the difference between assumed and ac-
less fouling tendency, less flow induced vibrations,
tual values is within 0.5%, then algorithm goes to the
less floor space requirement, compactness of design
next step. Otherwise it assigns new U = U
assumed calc
etc are some of the criteria which must be considered
and iterates the step 2 until the U and U
assumed calc
in industrial scenario. All these solutions are feasible
values are within 0.5% agreement.
and user has flexibility to choose any one of them
5. The computed values of flow velocities and
based on his requirement and engineering judgment.
the constructive details of the exchanger structure are
This paper collects some practical guidelines from li-
then used to evaluate the cost objective function. The
terature regarding how to choose the best exchan-
optimization algorithm, based on the value of the ob-
gers among various alternatives.
jective function, updates the trial values of the optimi-
zation variables which are then passed on to the de-
THE OPTIMAL HEAT EXCHANGER DESIGN
sign routine to define a new architecture of the heat
PROBLEM
exchanger. The process is iterated until a minimum of
the objective function is found or a prescribed conver-
The procedure for optimal heat exchanger de-
gence criterion is met.
sign includes the following step:
a. Estimation of the exchanger heat transfer Heat exchanger design procedure
area based on the required duty and other design This section describes step by step calculation
specifications assuming a set of design variable. procedure to evaluate heat exchanger area, tube side
411
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

Figure 1. Flowchart for heat exchanger design using SA algorithm.
and shell side pressure drop and total cost of heat (T −T )−(T −T )
LMTD= hi co ho ci (1)
exchanger. Equations (1) to (28) can be found from T −T 
ln hi co
Kern [14], Sinnot [15], Caputo et al. [6] as well as Tho−Tci
Patel and Rao [9].
The logarithmic mean temperature difference For sensible heat transfer, the heat transfer rate
(LMTD) is determined by: is given by:
412
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

Table 1. Search optimization variables and their options
Optimization Variable Total options
Variable name Corresponding options
variable notation available
x d Tube diameter, m 12 [0.00635, 0.009525, 0.0127, 0.01905, 0.022225, 0.0254, 0.03175,
1 0
0.0381, 0.04445, 0.0508, 0.05715, 0.0635].
x L Tube length, m 8 [1.2192 1.8288 2.4384 3.048 3.6576 4.8768 6.096 6.7056 7.3152]
2
x R Ratio of baffle spacing to 17 [0.2 ,0.25, 0.3, 0.35, 0.4, 0.45, 0.5 ,0.55, 0.6 ,0.65, 0.7, 0.75, 0.8,
3 b
shell diameter, (-) 0.85, 0.9, 0.95, 1]
x N Number of tube pass, (-) 5 1-1, 1-2, 1-4, 1-6, 1-8
4 pass
x P Type of pitch, (-) 2 Triangular and square
5 type
x Headtype Head type, (-) 4 Fixed tube sheet or U tube, outside packed head,
6
split ring floating head, and pull through floating head
x Bafflecut Baffle cut, % 4 [15, 25, 35, 45]
7

Q=mC (T −T ) (2) m n
h ph hi ho v = t ( ) (9)
The correction factor, F, for the flow configuration t π4dt2ρt Nt
involved is found as a function of dimensionless tem-
perature ratio for most flow configurations of: The Darcy friction factor f is calculated as fol-
t
lows:
1−P
F= R2+1 ln1−PR (3) ft=(1.82log 10Ret − 1.64)−2 (10)
R−1 ln(2−P(R+1− R2+1))
where Re is the tube side Reynolds number and
t
where the correction coefficient R is given by: given by:
T −T ρvd
R= hi ho (4) Re =  t t i (11)
T −T t μ
co ci t
Efficiency P is given by: where d is the inside diameter of tube and for sim-
i
plicity kept constant as d = 0.8d in this work. How-
T −T i 0
P= co ci (5) ever, for more rigorous calculations, instead of keep-
T −T
hi ci ing it constant, variable pipe thickness can be used
Assuming the overall heat transfer coefficient from tables of tube manufacturers: schedules 40 or 80
U , the heat exchanger surface S is computed or standard tubing gages: B.W.G. and Stub’s gage.
assumed
by: Prt is the tube side Prandtl number given by:
S= Q (6) Pr = μtCpt (12)
U FLMTD t k
assumed t
Tube side calculations. Number of tubes, N, and According to flow regime, the tube side heat
t
tube bundle diameter, Db, are calculated as follows: transfer coefficient (ht) is computed from following
correlation for Re ≤ 10000:
t
S
Nt=πd0L (7) h =kt[ (ft/8)(Ret−1000)Prt (1+ di)0.67] (13)
t d 1+12.7(f /8)1/2(Pr2/3− 1) L
1 i t t
D =d Ntn1 and for R et > 100 00: (8)
b 0K1 k μ
h =0.027 t Re0.8Pr1/3( t )0.14 (14)
where K and n are coefficients that take values ac- t d t t μ
1 1 o wt
cording to flow arrangement and number of passes as
Shell side calculation
per table given by Sinnot [15].
Velocity through tubes is found by: Clearance between tube bundle diameter and
shell diameter is calculated from the figure given by
Sinnot [15] for different head types as follows:
413
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

Clearance=D −D =mD +c where coefficients m and c are evaluated from figure
s b b
given by Sinnot [15] for different baffle cuts, as fol-
where m and c are empirical constants and assume lows: 15% baffle cut (-6.7×10-8; 0.100067); 25% baffle
following values for different head type: cut (-4.5×10-8; 0.070045); 35% baffle cut (-4.4×10-8;
– fixed and U tube (0.01, 0.008), outside packed 0.063044); 45% baffle cut (-3.4×10-8; 0.050034).
head (0.0, 0.038), split ring floating head (0.027,
The overall heat transfer coefficient is calculated
0.0446), pull through floating head (0.009, 0.0862).
as follows:
Baffle spacing is calculated as:
1
B = R D (15) U = (23)
bs s calc  1  d  1 
Cross-sectional area normal to flow direction is  +R + ( o)(R +  )
h  fs d ft h 
determined by: s i t
d Compare and check the absolute error,%,
S = DB(1− 0) (16)
a s P between Uassumed and Ucalc as per following:
t
U −U
where P is tube pitch and given by P = 1.25 d. Error [%]= 100 assumed calc  (24)
t t 0 U
Flow velocity for shell side can be obtained from: assumed
m If error is more than 0.5%, assume new U as
vs = (ρSs ) (17) Uassumed = Ucalc and go back to Eq. (5) and recalculate
s a everything until error reaches less than 0.5%.
Reynolds number and Prandtl number for shell When error is less than 0.5%, i.e., Uassumed and
side can be calculated as: Ucalc is very near to each other proceeding for further
calculations as follows:
mDe
Re =  s s (18) – Exchanger area:
s μS
s a Q
S= (25)
μC U FLMTD
Pr =  s ps (19) calc
s k The tube side pressure drop includes distributed
s
pressure drops along the tube length and concen-
where De is the shell hydraulic diameter and com-
s trated pressure losses in elbows and in the inlet and
puted as:
outlet nozzles:
– For triangular pitch:
4(0.43S2− (0.5πd2/4)) ΔPt= ρtv2t2(dLft+2.5)n (26)
De =  t 0 (20) i
s 0.5πd
0 The shell side pressure drop is given by:
– For square pitch: ρv2 L D
ΔP =f s s s (27)
4(S2− (πd2 /4)) s s 2 BDes
De =  t 0 (21)
s πd where fs is given by:
0
f =2b Re−0.15 (28)
Kern’s formulation for segmental baffle shell- s 0 s
and-tube exchanger is used for computing shell side
where b = 0.72 valid for Re < 40000.
0 s
heat transfer coefficient h:
s
Objective function
j K RePr0.333
hs =  h s Des s (22) Total cost Ctot is taken as the objective function,
s which includes capital investment (Ci), energy cost
(C), annual operating cost (C) and total discounted
where coefficient j is calculated from the figure given e o
h operating cost (C ):
by Sinnot [15] for different baffle cuts. od
The following approximate equation was used: C =C +C (29)
tot i od
j =mRec
h s
414
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

Adopting Hall’s correlation, the capital invest- Constraints
ment (C) is computed as a function of the exchanger
i Though the lowest cost exchanger is the main
surface area:
selection criterion for STHEs, this is not the only cri-
Ci =a1+a2Sa3 (30) terion for commercial plants. The concept of a good
design involves aspects that cannot be easily des-
where a1 = 8000, a2 = 259.2 and a3 = 0.93 for ex- cribed in a single economic objective function e.g.
changer made with stainless steel for both shell and fouling suppression, maintenance ease, mechanical
tubes. resistance, simplicity, flow distribution, potential tube
The total discounted operating cost related to vibration etc. These criteria, though subjective, have a
pumping power to overcome friction losses is com- profound effect on exchanger performance in com-
puted from the following equation: mercial plants. These criteria are sometimes expres-
C =PC H (31) sed as geometric and hydraulic and service cons-
0 e traints [17].
C =ny C0 (32) Geometric constraints
od x=1(1+i)x The STHE candidate must respect a series of
geometric constraints, involving the following rules:
where pumping power P is computed from:
the ratio between tube length and shell diameter must
1 m m be between 3 and 15, the ratio between baffle spa-
P = ( t ΔP + s ΔP ) (33)
η ρ t ρ s cing and shell diameter must be between 0.2 and 1;
t s
the baffle spacing cannot be lower than 50 mm; the
Based on all above calculations, total cost is baffle spacing must obey the maximum unsupported
computed from Eq. (29) for case studies 1 and 2. The span. These constraints can be represented by the
cost objective function is kept exactly same as case following mathematical expressions:
study by GA approach by Caputo et al. [6] to enable
3≤L/D ≤5 (37)
the performance comparison of SA approach and GA s 
approach. 0.2≤R ≤1 (38)
However, for case study 3 the authors did not bs
use same cost function in their case studies. Some of 0.050D /≤R ≤0.5Lmax /D (39)
the authors [16] use total annual cost slightly diffe- s  bs b s
rently as follows: the total cost consists of five com- Velocity constraints
ponents: the capital cost of the exchanger, the capital
These constraints represent fluid velocity limits
costs for two pumps, and the operating (power) costs
in order to reduce fouling and erosion problems. The
of the pumps. The expression for the total annual cost
tube and shell side velocity must obey lower and up-
is of the form:
per bounds:
C i=Af(Cexc +Cpump,T)+Cpump,s)= νmin ≤ν ≤νmax (40)
t   t  t
m
=A [(C +C Ac)+(C +C ( t ΔP)e)+ (34)
f a b e f ρ t νmin ≤ν ≤νmax (41)
t s   s  s
m
+(C +C ( s ΔP )e)] High velocities will give high heat transfer coef-
e f ρ s
s ficients but also a high pressure drop. The velocity
must be high enough to prevent any suspended so-
where C , C and C are the capital costs for
exc pump,t pump,s
lids settling, but not so high as to cause erosion. High
the exchanger, tube side and shell side pumps,
velocities will reduce fouling. Plastic inserts are some-
respectively.
timesused to reduce erosion at the tube inlet. Typical
C = 1 (35) design velocities are given below [15].
od η(mt ΔP + ms ΔP )C H Liquids. Tube-side, process fluids: 1 to 2 m/s,
ρ t ρ s pow
t s maximum 4 m/s if required to reduce fouling; water:
1.5 to 2.5 m/s; shell-side: 0.3 to 1 m/s.
C =C +C (36)
tot i od Vapors. For vapors, the velocity used will de-
pend on the operating pressure and fluid density; the
This new cost calculation is used in case study 3
lower values in the ranges given below will apply to
to have a same basis for comparison of performance
high molecular weight materials; vacuum: 50 to 70
of SA approach and GA approach.
415
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

m/s; atmospheric pressure: 10 to 30 m/s; high pres- where x is the vector of optimization variables (Table
sure: 5 to 10 m/s. 1). The set of constraints g(x) corresponds to the
In this work, the following constraints were im- inequalities given by Eqs. (37)–(43).
posed on objective function: For implementation of the SA algorithm, we
used a penalty function in the objective function, to
1≤v ≤2 m/s and 0.3≤v ≤1 m/s
t s provide the following objective function to be mini-
Service constraints mized [16].
The hydraulic requirements of the service are Obj(x)=C (x)+penalty(x) (44)
tot
represented by upper bounds on the pressure drop of
The penalty function accounts for the violation of
both streams:
the constraints such that:
ΔP ≤ΔPmax (42)
t  t  0 if x is feasible

ΔP ≤ΔPmax (43) Penalty(x)=m (45)
s  s rg2(x) otherwise
 i i
i=1
The values suggested below can be used as a
general guide, and will normally give designs that are where r is a variable penalty coefficient for the ith
i
near the optimum. [15]. constraint, r varies according to the level of violation.
i
Liquids. For fluids having Viscosity <1 mN s/m2
Simulated annealing: at a glance [22]
maximum pressure drop: 35 kN/m2 [15]. For fluids
having viscosities 1 to 10 mN s/m2, maximum pres- What Is simulated annealing?
sure drop: 50-70 kN/m2 [15]. The simulated annealing method is based on
Gas and vapors. High vacuum: 0.4-0.8 kN/m2; the simulation of thermal annealing of critically heated
medium vacuum: 0.1×absolute pressure; 1 to 2 bar: solids. When a solid (metal) is brought into a molten
0.5×system gauge pressure; above 10 bar: 0.1×system state by heating it to a high temperature, the atoms in
gauge pressure. the molten metal move freely with respect to each
For the present case study following constraints other. However, the movements of atoms get restric-
were imposed: ted as the temperature is reduced. As the tempera-
ΔP ≤ 35000 Pa; ΔP ≤ 35000 Pa ture reduces, the atoms tend to get ordered and fi-
t s
nally form crystals having the minimum possible inter-
The prerequisite of a good design is to choose nal energy. The process of formation of crystals es-
the lowest cost exchanger with standard dimensions sentially depends on the cooling rate. When the tem-
(as per TEMA standard) while obeying the above perature of the molten metal is reduced at a very fast
constraints. Attempt has been made in this work to rate, it may not be able to achieve the crystalline state;
apply SA optimization technique to design a lowest instead, it may attain a polycrystalline state having a
cost heat exchanger with TEMA dimensions and sa- higher energy state compared to that of the crystalline
tisfying all of the above constraints. state. In engineering applications, rapid cooling may
However, the value of the constraints of pres- introduce defects inside the material. Thus, the tem-
sure drop and velocity is dependent on the detailed perature of the heated solid (molten metal) needs to
design and very much problem specific. In this work, be reduced at a slow and controlled rate to ensure
the values of constraints are selected as per general proper solidification with a highly ordered crystalline
guidelines given by Sinnot [15] and the user is not state that corresponds to the lowest energy state (in-
restricted to adhere this value. The value of these ternal energy). This process of cooling at a slow rate
constraints must be judiciously selected as they have is known as annealing.
a big impact on final solution and cost. In case the
Procedure [22]
user does not have specific restriction on these va-
lues, the constraints should be kept as broad as pos- The simulated annealing method simulates the
sible. This will facilitate the lowest cost heat exchanger. process of slow cooling of molten metal to achieve
the minimum function value in a minimization prob-
Handling the constraints
lem. The cooling phenomenon of the molten metal is
The original problem can be set as: simulated by introducing a temperature-like parame-
Minimize Ctot(x) ter and controlling it using the concept of Boltzmann’s
Subject to gi(x) ≥ 0 where i = 1, 2, m probability distribution. The Boltzmann’s probability
distribution implies that the energy (E) of a system in
416
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

thermal equilibrium at temperature T is distributed Note that the probability of accepting the point X is
i+1
probabilistically according to the relation: not same in all situations. As can be seen from Eq.
(50), this probability depends on the values of ∆E and
P(E)=e−kET (46) T. If the temperature T is large, the probability will be
high for design points X with larger function values
i+1
where P(E) denotes the probability of achieving the
(with larger values of ∆E = ∆f). Thus at high tempe-
energy level E, and k is called the Boltzmann’s cons-
ratures, even worse design points X are likely to be
i+1
tant. Equation (46) shows that at high temperatures
accepted because of larger probabilities. However, if
the system has nearly a uniform probability of being
the temperature T is small, the probability of accept-
at any energy state. However, at low temperatures,
ing worse design points X (with larger values of ∆E =
i+1
the system has a small probability of being at a high-
= ∆f) will be small. Thus as the temperature values
energy state. This indicates that when the search
get smaller (that is, as the process gets closer to the
process is assumed to follow Boltzmann’s probability
optimum solution), the design points X with larger
i+1
distribution, the convergence of the simulated anneal-
function values compared to the one at X are less
i
ing algorithm can be controlled by controlling the tem-
likely to be accepted.
perature T. The method of implementing the Boltz-
mann’s probability distribution in simulated thermody- P(Ei+1)=e−kΔTE (50)
namic systems, suggested by Metropolis et al. [18]
can also be used in the context of minimization of Algorithm [22]
functions. In the case of function minimization, let the
The SA algorithm [22] can be summarized as
current design point (state) be Xi, with the corres- follows. Start with an initial design vector X (iteration
1
ponding value of the objective function given by fi = number i = 1) and a high value of temperature T. Ge-
= f(Xi). Similar to the energy state of a thermodynamic nerate a new design point randomly in the vicinity of
system, the energy E at state X is given by:
i i the current design point and find the difference in
E =f =f(X) (47) function values:
i i i
ΔE=Δf −Δf =f(X )−f(X) (51)
Then, according to the Metropolis criterion, the pro- i+1 i i+1 i
bability of the next design point (state) X depends
i+1
If f is smaller than f (with a negative value of
on the difference in the energy state or function va- i+1 i
∆f), accept the point X as the next design point.
lues at the two design points (states) given by: i+1
Otherwise, when ∆f is positive, accept the point X
i+1
ΔE=ΔEi+1−ΔEi =f(Xi+1)−f(Xi) (48) as the next design point only with a probability
exp(−ΔE/kT). This means that if the value of a ran-
The new state or design point X can be found
i+1 domly generated number is larger than exp(−ΔE/kT),
using the Boltzmann’s probability distribution:
accept the point X ; otherwise, reject the point X .
i+1 i+1
P(Ei+1)=min1,e−kΔTE (49) Tthhei sp cooinmt pXlei+t1e sis orneeje citteerda, titohne no ft hthee p rSoAc easlsg oorfi thgmen. eI-f
rating a new design point X randomly in the vicinity
i+1
The Boltzmann’s constant serves as a scaling
of the current design point, evaluating the corres-
factor in simulated annealing and, as such, can be
ponding objective function value f , and deciding to
i+1
chosen as 1 for simplicity. Note that if ∆E ≤ 0, Eq. (49)
accept X as the new design point, based on the use
i+1
gives P(E ) = 1 and hence the point X is always
i+1 i+1 of the Metropolis criterion, Eq. (50), is continued. To
accepted. This is a logical choice in the context of mi-
simulate the attainment of thermal equilibrium at
nimization of a function because the function value at
every temperature, a predetermined number (n) of
X , f , is better (smaller) than at X, f, and hence the
i+1 i+1 i i new points X are tested at any specific value of the
i+1
design vector X must be accepted. On the other
i+1 temperature T. Once the number of new design points
hand, when ∆E > 0, the function value f at X is
i+1 i+1 X tested at any temperature T exceeds the value of
i+1
worse (larger) than the one at X. According to most
i n, the temperature T is reduced by a pre-specified
conventional optimization procedures, the point X
i+1 fractional value c (0 < c < 1) and the whole process is
cannot be accepted as the next point in the iterative
repeated. The procedure is assumed to have con-
process. However, the probability of accepting the
verged when the current value of temperature T is
point X , in spite of its being worse than X in terms
i+1 i sufficiently small or when changes in the function va-
of the objective function value, is finite (although it
lues (∆f) are observed to be sufficiently small. The
may be small) according to the Metropolis criterion.
choices of the initial temperature T, the number of
417
N.M. KHALFE, S.K. LAHIRI, S.K. WADHWA: SIMULATED ANNEALING TECHNIQUE… CI&CEQ 17 (4) 409−427 (2011)

iterations n before reducing the temperature, and the 1. The quality of the final solution is not affected
temperature reduction factor c play important roles in by the initial guesses, except that the computational
the successful convergence of the SA algorithm. For effort may increase with worse starting designs.
example, if the initial temperature T is too large, it re- 2. Because of the discrete nature of the function
quires a larger number of temperature reductions for and constraint evaluations, the convergence or
convergence. On the other hand, if the initial tempe- transition characteristics are not affected by the
rature is chosen to be too small, the search process continuity or differentiability of the functions.
may be incomplete in the sense that it might fail to 3. For problems involving behavior constraints
thoroughly investigate the design space in locating (in addition to lower and upper bounds on the design
the global minimum before convergence. The tempe- variables), an equivalent unconstrained function can
rature reduction factor c has a similar effect. Too be formulated as in the case of genetic algorithms.
large a value of c (such as 0.8 or 0.9) requires too
much computational effort for convergence. On the
other hand, too small a value of c (such as 0.1 or 0.2)
may result in a faster reduction in temperature that
might not permit a thorough exploration of the design
space for locating the global minimum solution. Si-
milarly, a large value of the number of iterations n will
help in achieving quasi equilibrium state at each tem-
perature but will result in a larger computational effort.
A smaller value of n, on the other hand, might result
either in a premature convergence or convergence to
a local minimum (due to inadequate exploration of the
design space for the global minimum). Unfortunately,
no unique set of values is available for T, n, and c that
will work well for every problem. As per recommen-
dation of Rao [22], the initial temperature T can be
chosen as the average value of the objective function
computed at a number of randomly selected points in
the design space. The number of iterations n can be
chosen between 50 and 100 based on the computing
resources and the desired accuracy of solution. The
temperature reduction factor c can be chosen be-
tween 0.4 and 0.6 for a reasonable temperature re-
duction strategy (also termed the cooling schedule).
More complex cooling schedules, based on the ex-
pected mathematical convergence rates, have been
used in the literature for the solution of complex prac-
tical optimization problems. In spite of all the research
being done on SA algorithms, the choice of the initial
temperature T, the number of iterations n at any spe-
cific temperature, and the temperature reduction fac-
tor (or cooling rate) c still remain an art and generally
require a trial-and-error process to find suitable va-
lues for solving any particular type of optimization
problems. The SA procedure is shown as a flowchart

in Figure 2 [22].
Figure 2. Simulated annealing procedure [22].
Features of the method
Some of the features of simulated annealing as SA implementation
found from literature [22] are as follows:
The objective function is the minimization of HE
cost Obj(x) given in Eq. (44) and x is a solution string
representing a design configuration.
418