Optimal Steady-State Design of Reactive Distillation Processes Using Simulated Annealing

37 Pages·2008·0.48 MB·English

Optimal Steady-State Design of Reactive Distillation Processes
Using Simulated Annealing
Jian-Kai Cheng, Hao-Yeh Lee, Hsiao-Ping Huang, Cheng-Ching Yu*
Department of Chemical Engineering, National Taiwan University, Taipei 106,
Taiwan
Abstract - In this paper, we investigate on a derivative-free optimization approach,
Simulated Annealing (SA), for the optimization of the reactive distillation (RD)
column design. Because RD systems exhibit non-monotonic behavior for key design
variables, flowsheet optimization using simulator is difficult when derivative-based
approach is employed. The SA-based optimization procedure gives an equally good
or better design than the optimal flowsheet obtained from the sequential design
approach. More importantly, this is achieved with much more efficient computing.

Keywords: reactive distillation; simulation, simulated annealing; process design;
optimization

* corresponding author, e-mail: [email protected]; tel: [phone], fax: [phone]
1. Introduction
Reactive distillation (RD) offers an attractive alternative for the production of
many important industrial chemicals. As pointed out by Hung et al. (2006), the
multifunctional nature leads to input multiplicity which makes the simulation and
convergence of the RD flowsheet difficult. Unlike conventional separation units, there
are much more design variables, e.g., number of separation trays, number of reactive
trays, feed tray location etc., associated with the RD column. The design of the RD
column becomes a combinatorial optimization problem where the objective function
is the total annual cost (TAC). Conventionally, a sequential design approach is taken
where the variable is changed one at a time to minimize the design objective (Kaymak
and Luyben, 2004a, 2004b and 2004c; Tang et al., 2005, Hung et al., 2006; Lai et al.,
2007; Tung and Yu, 2007; Lin et al., 2008). However, as the dimensionality of design
variables increases, the sequential method approach imposes heavy computing load,
especially for the coupled reactive distillation/ separation columns processes (Lai et
al., 2007; Lin et al., 2008). Mathematical optimization algorithms provide an
attractive alternative to improve the computational efficiency.
Numerical optimization algorithms have been studied extensively in the
literature and optimization problems are categorized into two categories: convex and
non-convex. The former can be solved by the mathematical programming to obtain
the global optimum (Grossmann, 2002; Biegler and Grossmann, 2004). However, the
derivative-based mathematical programming cannot guarantee the global solution for
the non-convex problem. Unfortunately, the optimization problems of the RD design
are non-convex, and derivative-free algorithms, e.g., Simulated Annealing (SA),
Genetic Algorithm (GA), Tabu search (TS) (Linke et al., 2003) etc, are preferred.
Among all of them, SA has the advantage that it is simple to program. The concept of
SA was originally proposed by Metropolis et al. (1953). It uses Monte Carlo
simulation to calculate the energy distribution of molecules and simulate the cooling
of material in a heat bath. Kirkpatrick et al. (1983) is among the first to apply this idea
to optimization and a traveling salesman problem was solved successfully using SA.
After that, many applications of this method were proposed for the optimization of
chemical processes. Examples include: heat exchanger networks (Dolan et al., 1989),
process synthesis (Lo, 2006), the hydrodealkylation process (Chaudhuri and Diwekar,
1996), and the reactive distillation (Cardoso et al., 2000; Lima et al., 2006; Gomez et
al., 2006). Adaptive SA-based methods have also been proposed (Cardoso et al., 2000;
Lima et al., 2006). Suman and Kumar (2006) gives a comprehensive review of
SA-based optimization algorithms which solve single and multiple objective
optimization problems.
In recent years, process simulators have been used extensively for the simulation
of chemical processes. Because of the “black box” nature of these simulators,
derivative-free algorithms can easily be implemented for optimization problems
(Diwekar et al., 1992; Dantus and High, 1999; Leboreiro and Acevedo, 2004; Lo,
2006). In this paper, we investigate the application of a SA algorithm to obtain the
optimal design of two reactive distillation systems. We adopt a conceptually simple
cooling schedule that is frequently used in the literature. The SA algorithm is
implemented in the Visual Basic Application (VBA) which interfaces with the
process simulator, Aspen Plus.

2. Simulated Annealing
The concept of the simulated annealing is based on the manner in which melts
freeze or metals grow crystals in the process of annealing. Initially, molecules of
melts are disordered at high temperature and they are slowly cooled so that the system
at each time instance can be viewed as in thermodynamic equilibrium. As the
temperature decreases, the system becomes more ordered and approaches a "frozen"
ground state at T=0. If the system is quenched too rapidly, a metastable structure
forms. At each temperature T, if the system is in thermodynamic equilibrium, it can be
described by the probability (P) of being in an energy state (E) at a temperature (T)
which is given by the Boltzmann energy distribution:
1
P{E}= exp(−E/K T) (1)
Z(t) b
where K is the Boltzmann’s constant and 1/Z(t) is a normalization factor.
b
In statistical mechanical simulations of annealing, molecules are allowed to
move randomly to a new position. Based on the Boltzmann energy distribution, we
calculate the probabilities of the old and new configurations and the ratio of the
probabilities is defined as the acceptance probability.
P{E } ⎛ E −E ⎞ ⎛ ΔE⎞
P = new =exp⎜− new old ⎟=exp⎜− ⎟ (2)
acc P{E } ⎝ K T ⎠ ⎝ T ⎠
old b
where ΔE is the difference of the normalized energy between the new and the
old states. The new energy state is accepted if one of the following conditions is met.
If the energy of the new configuration is less than that of the original, the system
accepts this new configuration. If the energy of the new configuration is greater than
the old one, the system accepts the new configuration only if the probability,
exp(-ΔE/T), is greater than a random value generated from a uniformly distributed
random number. This is based on the Metropolis criterion and it implies that the
system may accept high energy states at a higher temperature. As the temperature is
lowered, the system tends to prefer lower energy states. The mathematical algorithm
for optimization is analogous to the annealing process. An important feature of the SA
method is “probabilistic uphill moves”. In such a manner, the SA algorithm reduces
the chance of getting trapped in local minimum with a wide random move. Detailed
steps of the SA algorithm are described in Figure 1. In terms of process optimization
with a set of design variables (d), the energy state (E) is equivalent to the objective
function to be minimized. That is E is a function of d, i.e., E(d), and the energy state
will be lowered as the annealing starts. Initially, we pick a random vector of the
design variable d, the initial system temperature (T ), the final temperature (T),
i 0 f
criteria for thermodynamic equilibrium (e.g., N ), and the cooling schedule (α).
T
1. Perturb d to obtain a neighboring design vector d .
i i+1
2. Evaluate E(d ) using a simulation model and calculate ΔE= E(d ) - E(d).
i+1 i+1 i
3. (a) If ΔE < 0, then accept d as the new current solution.
i+1
(b) If ΔE > 0 and exp(-ΔE/T) > U(0,1), then accept d as the new current
i+1
solution. Here U(0,1) is a random number between 0 and 1.
4. Return to 2 until a thermodynamic equilibrium at temperature T is achieved.
5. Reduce the system temperature according to the cooling schedule (i.e., T=
αT) and reset i=1 and then go back to 1.
6. Terminate the algorithm when the frozen temperature is reached, i.e., T≦T.
f
Before implementing the SA algorithm for the optimization problem, we need to
determine the annealing schedule, which affects the convergence of the algorithm and
CPU time. Aarts and Krosr (1989), Patel et al. (1991) and Painton and Diwekar (1994)
discussed about the parameter settings when using the simulated annealing algorithm.
Due to its simplicity, some of basic parameter settings are described next.
2.1 Initial Temperature (T )
0
In the application of the SA algorithm, energy (E) and temperature (T) has the
same unit because we have combined Boltzmann’s constant and temperature shown in
the equation of Boltzmann energy distribution. The purpose of testing the initial
temperature is to make sure that there is a high enough acceptance ratio at the start of
cooling. There are two criteria for the determination of the initial temperature:
(a) Start at an arbitrary initial temperature T and attempting few hundred moves. If
0
the fraction of the moves accepted is greater than 80%, T is suitable for the system. If
0
not, double T and repeats the procedure.
0
(b) Take a few hundred random moves and find the largest difference of the energy
state function (ΔE) within these moves, and set 10 times this ΔE as the value of T .
0
2.2 Final Temperature
The final (frozen) temperature is the termination of an annealing process. If T is
f
too small, it will give rise to significant CPU time. There are two common criteria for
the determination of the final temperature:
(a) Set the final temperature to a given value. The algorithm stops when the system
reaches T.
f
(b) The final temperature is chosen such that the algorithm stops after 10 successive
temperature decrements with no change in the optimal configuration.
2.3 Temperature Decrement
If the temperature decreases too fast, the energy state will be too low for the
algorithm to escape from a local minimum. There are two common criteria for the
determination of the temperature decrement:
(a) T = αT , where 0.8 ≤ α≤ 0.9
n+1 n
(b) T = T [1+ln(1+γ)T /3σ]-1, where σ is the standard deviation of the cost at T 
n+1 n n n
and γ is a decreasing parameter (cid:4666)the smaller γ the slower cooling (cid:4667).
2.4 Equilibrium Detection
We assume that cooling process reaches quasi-equilibrium for every temperature
throughout the annealing process. There are two criteria:
(a) Reach N moves at each temperature.
T
(b) Achieve an accept/reject ratio of 1:10.
2.5 Move Generator
The limits of the search space are defined by the user. In this work, integer and
continuous design variables have maximum step sizes.
(a) Integer variables: dnew =dold +INT{(2⋅U(0,1)−1)N }
step
(b) Continuous variables: dnew =dold +{(2⋅U(0,1)−1)N }
step
where N is the vector of maximum step sizes for each move and INT means to
step
round off the number.
In this work, the first criterion of each parameter setting is applied.
3. Applications
In this section we present two examples. One is the reactive distillation to
produce methyl acetate (MeAc) and the other is the production of butyl acetate
(BuAc). Both of them have been studied by Tang et al. (2005) using the sequential
optimization approach. All the simulations are carried out in Aspen Plus using the
RADFRAC module with FORTRAN subroutines for the activity-based reaction
kinetics.

3.1 Example 1 – Methyl Acetate System
In the methyl acetate (MeAc) system, we assume a “neat” design (Tung and Yu,
2007) where exact stoichiometric feeds are introduced to the column. Thus the
dominant design variables are: number of rectifying trays (N ), number of reactive
R
trays (N ), number of stripping trays (N ), feed location of acetic acid (NF ) and
RX S HAc
feed location of methanol (NF ), as shown Figure 2. In order to reduce the
MeOH
occurrence of infeasible solutions, we give constraints on the design variables: 0≦N
R
≦10, 25≦N ≦45, 0≦N ≦10, (NF －NF )≧10. Reflux ratio and rebolier
RX S MeOH HAc
duty are adjusted automatically in Aspen Plus to achieve the product purity
specifications (x = 98 mol% and x = 98 mol%). Note that the tray number is
MeAc H2O
counted from the up- bottom.
3.1.1 Implementation
In our simulation, VBA is used to implement the simulated annealing algorithm,
as shown Figure 3. The optimization procedure consists of the following steps. First,
the VBA generates a set of design variables randomly, and then they are passed to
Aspen Plus as design variables for the process simulation. At the completion of the
simulation, the results are transferred back to the VBA program which calculates the
energy state (Total Annual Cost, the TAC calculation is listed in Appendix A) and the
design variables are updated according to the SA algorithm. The steps are repeated
until minimal energy state is reached. The VBA code is listed in the Appendix B.
Before getting into the SA optimization, we explore the effects of initial
temperature settings. The results of Table 1 clearly indicate that as the initial
temperature increases, the acceptance ratio increases. At T = 60, the acceptance ratio
0
reaches 80% as shown in Table 1. Note that T has the unit of $/yr. Thus, we choose
T =60 as the initial temperature.
0
3.1.2 Results and Discussion
Three different cases (case A, B, C and D) are used to illustrate the effects of the
SA parameter settings, such as α, T, N and N , to the optimization results. In order
f T step
to achieve a more general view on the SA algorithm, for each case, the SA schedule is
executed 3 times (3 runs in case D are tested by different initial conditions of design
variable set) and the results are given in Table 2. The results show that, for all cases, 3
runs in case D obtain the global optimum. It clearly indicates that reliable design can
be achieved using the SA algorithm. Figure 4 reveals that the objective function (TAC)
tends to fluctuate at high temperature state (to avoid being trapped in a local minima),
and improves steadily at low temperature states. The effects of specific parameter
values are discussed next.
(a) Effect of N (case A and B): When fewer moves are allowed at each
T
temperature, the result gives a larger deviation from the global optimal.
Acceptable quasi-equilibrium is achieved with 15 moves for the MeAc system.
(b) Effect of α and T (case B and C): If the rate of cooling is too fast, the result is
0
a suboptimal metastable configuration. As mentioned earlier, the outcome with a
low T is similar to that of a small α.
0
(c) Effect of T: We use a very small value of T (0.00001) in the simulation. Figure
f f
4 indicates that the TAC stabilizes when T reaches 0.001. This suggests that the
f
termination temperature of T=0.001 is a suitable choice.
f
(d) Effect of N (case B and D): Comparing N =1 with N =3, the variable
step step step
generation with N =1 has better results. In the low temperature, the algorithm
step
has more chances to find smaller TAC in the finite moves (N ), when the new
T
design variables are perturbed with the smaller neighborhood (N =1).
step
(e) Effect of initial condition (case D): 3 runs starting with different initial design
variables have the same optimal results. It turns out that the initial condition
doesn’t affect the outcome of the SA algorithm because of the characteristics of
randomly generated variables.
We compare these results with two other optimization approaches. One method
is evaluating all possible flowsheets in the design space and the other is a sequential
optimization procedure of Tang et al. (2005). The drawback of the former is that it
requires significant computing time. The drawback of the later method is that the
outcome can be affected by the initial condition and sequence in executing
optimization with one variable at a time. The results indicate that the SA method uses
less computing time while giving a satisfactory design. Table 3 shows that the
proposed SA algorithm is more efficient than the exhausted search by almost 10 folds.
3.1.3 Cases for distributed feeds
In this section, we explore the effects of feed splits (splitting the acetic acid feed
and methanol feed to different trays) to the designs. Figure 5 show the feed-split
configuration for the heavy and the light components. Thus, we have two more
optimization variables (SR and SR ), and they are continuous variables. Similar
HAc MeOH
to the earlier case, we test different initial temperatures and choose T = 120 as the
0
initial temperature. The other settings for the SA algorithm are: α=0.8, T =0.0001 and
f
N =15. The SA schedule is executed 3 times and the results are given in Table 4. It is
T
interesting to note that no economic incentive is gained by splitting the acetic acid
feed, and a slight TAC decrease is observed by feeding methanol to different trays.
3.2 Example 2 – Butyl Acetate System
In the BuAc system, we also assume the “neat” design where exact
stoichiometric feeds are introduced to system. Thus the dominant variables for
optimization are the same as that of the MeAc system (Figure 6). In order to reduce
the occurrence of infeasible solutions, the following constraints are imposed on the
design variables: (1≦N ≦10; 15≦N ≦35; 1≦N ≦10; NF －NF ≧2). The
R RX S MeOH HAc
organic phase is totally refluxed and the the rebolier duty is adjusted to maintain the
product purity specification (x = 99 mol% and x ≦ 50 ppm). As the previous
BuAc HAc
example, we test different initial temperatures and choose T = 240 as the initial
0
temperature. The other settings for the SA algorithm are: α=0.8, T =0.0001 and N =15.
f T
The SA schedule is executed 3 times and the results are given in Table 5. Again,
Figure 7 reveals that the objective function (TAC) tends to fluctuate at high
temperature state, and improves steadily at low temperature states. The best case
results in 25% improvement in the TAC as compared to that of Tang et al. (2005)
using the sequential optimization.
Because of the nonlinearity, the BuAc system has different results for 3 runs in
SA algorithm. Figure 8 shows that the profiles of TAC as the design variables (N and
S
N ) vary for the MeAc and the BuAc systems. Note that the TAC shown in Fig. 8
R
corresponding to the optimized feed locations. Figure 8 also reveals that, the MeAc
system has a much smooth surface as compared to that of the BuAc system. This
explains why it is more likely to be trapped at local minimum for the BuAc system. In
terms of the converging a process flowsheet, Table 6 shows that the BuAc system is
more difficult to converge as compared to the MeAc system. Thus, efficient
optimization method is crucial for the design BuAc types of processes.
4. Conclusion
Although the simulated annealing method cannot guarantee a globally optimal
design, it is computationally efficient and provides a reasonable engineering design.
There are two advantages: (1) Compared with the approach evaluating the entire
solution space, it takes less CPU time and (2) Compared with the sequential design
approach, it has a smaller probability of getting trapped in a local minimum. We have
used an SA algorithm for the design of two reactive distillation systems. The results
indicate that improved design can be achieved with relatively efficient computing. It
is recommended for highly nonlinear processes with non-monotonic characteristics
such as reactive distillation systems.

Acknowledgement
This work was supported by the Ministry of Economic Affairs and National
Taiwan University Excellence Research Project (NTU-ERP-95R0066-BE04-04). We
thank Professor S. H. Cheng of Tunghai University for helpful comments.