4158
Ind. Eng. Chem. Res. 2008, 47, 4158–4169
Time-Optimal Control of Startup Distillation Columns by Iterative Dynamic Programming Alexandru Woinaroschy* Department of Chemical Engineering, “Politehnica” UniVersity of Bucharest, Polizu Str. 1-5, Bucharest, Romania
Iterative dynamic programming employing randomly chosen candidates for admissible control is applied to the minimization of distillation startup time. The solution of this high-dimensional problem was facilitated mainly by avoiding iterative algebraic equation solving inside a rigorous dynamic distillation model, thereby reducing the computation time. Several illustrative applications of various complexity degrees are presented. The control variables are the reflux ratio, the reboiler heat duty, or both. Both piecewise constant control and piecewise linear control strategies were employed. The decrease of the startup times and heat consumption is high, having an appreciable economic significance. Introduction Time-optimal control (TOC) is a challenge in engineering: the general problem of reaching a desired final state from a set of given initial conditions in minimum time by an adequate control policy occurs very often in chemical engineering, especially in process startup. This transition should be accomplished as fast as possible because the intermediate products have no economic value. In many cases, the values of utilities and the time consumed during the startup period are important because, in fact, they are wasted. The first solutions for TOC problems were obtained for linear or linearized systems.1–3 For high-dimensional nonlinear systems, TOC problems are extremely difficult to solve. In recent decades, a promising method, respectively iterative dynamic programming (IDP), was used in TOC of nonlinear systems.4–8 Usually, methods based on dynamic programming were limited to systems with low dimensionality. By introduction of IDP by Luus9 that employs systematic contraction of the search region and thus dispenses the need for a fine grid, the dimensionality of the problem can be significantly increased. However, the applications presented in previous citations were mainly from the field of chemical reactors, involving a reduced number of state variables. TOC of startup distillation columns involves a large number of state variables, being a very high-dimensional problem. Fikar et al.10 have studied the application of IDP on optimal control of a distillation column for separation of the binary mixture methanol-isopropanol. They have used the top product molar flow rate and the reflux molar flow rate as control variables. IDP was applied to minimize the duration of transition between two stationary states defined byproduct compositions. They have used a distillation column model with several strong simplifications: negligible pressure variation through the column, constant relative volatility, constant molar liquid volumes, constant vaporization enthalpies, etc.
umini e ui e umaxi i ) 1, 2, ..., nu
The continuous dynamic system is described by differential equations (1)
* To whom correspondence should be addressed. Phone: +40-214023902. Fax: +40-21-3185900 E-mail: A_WOINAROSCHY@ chim.upb.ro.
(2)
a. Piecewise Constant Control. To transform the continuous control policy into a piecewise constant control problem, the time interval 0 e t e tf is divided into S time stages of variable length V(k) ) tk - tk-1 k ) 1, 2, ..., S
(3)
where the values of control variables are kept constant during each of these time intervals. The reason of variable length stages is to accurately determine switching times. By introduction5,6 of the normalized time variable τ, such that τ ) t/tf, with discrete values τk ) k/S, k ) 0, 1, ..., S in the transformed time domain, each stage is of equal length 1/ S. Equation 1 in the time interval τk e τ e τk+1 becomes dsi ) V(k)Sfi(s1, s2, ..., sns, u1, u2, ..., unu) i ) 1, 2, ..., ns; dτ k ) 1, 2, ..., S (4) The TOC problem is to minimize the final time tf and to determine the corresponding control variables ui (τk), i ) 1, 2, ..., nu ; k ) 1, 2, ..., S and length of time stages V(k), k ) 1, 2, ..., S restricted to si(tf) ) ssi i ) 1, 2, ..., ns
(5)
where ssi represents the desired final stationary state values. To minimize the final time tf subject to constraint eq 5, the penalty performance index was formulated as
∑ |1 ns
I ) tf + ω
Mathematical Formulation of TOC
dsi ) fi(s1, s2, ..., sns, u1, u2, ...unu) i ) 1, 2, ..., ns dt
where si respresents the state variables at a given initial state si(0), i ) 1, 2, ..., ns, and ui represents control variables subject to the bounds
i)1
ns
si(tf) ssi
|
(6)
where ω > 0 is a penalty coefficient and the final values of state variables si (tf), i ) 1, 2, ... ns are calculated by integration of eq 4. The TOC algorithm proposed by Bojkov and Luus5 was used in the following applications. This algorithm is based on the IDP procedure given by Bojkov and Luus,11,12 which employs randomly chosen candidates for the admissible control. The stop criterion is to execute a specified number of IDP iterations. The
10.1021/ie0713745 CCC: $40.75 2008 American Chemical Society Published on Web 05/22/2008
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4159
second term in the objective function (eq 6) ensures that the desired final stationary state (restrictions defined by eq 5). For this aim, the selected value of the penalty coefficient ω is very important. For a small value of this coefficient, the final stationary state will not be reached, whereas a high value will seriously increase the computing time. In the examples presented in this paper, the values of the penalty coefficient were between 5 × 105 and 106. b. Piecewise Linear Control. In some cases, a continuous control policy is preferred rather than a policy that requires sudden switching from one level to another. According to the procedure proposed by Luus,13 a piecewise linear control policy ui in the time interval (tk, tk+1) is given by ui(t) ) ui(k) +
(
)
ui(k + 1) - ui(k) (t - tk) tk+1 - tk
(7)
where ui(k) is the value of ui at tk and ui(k + 1) is the value of ui at tk+1. Then, the optimal control problem is to find ui(k), i ) 1, 2, ..., nu; k ) 0, 1, ..., S - 1, minimizing the performance index in eq 6. For the last stage ui is kept constant at the value ui(S - 1). By iterative dynamic programming, eq 7 does not introduce any difficulties because ui(k + 1) is known from the previous stage, and the remaining task is to find the best value for ui(k) at the beginning of the interval. To provide a continuous control policy, a single point is used for the s-grid at each time stage. Luus13 considered that, despite the convergence difficulties that may arise, it is a good compromise from the point of view of simplicity. Startup Distillation Operational Procedure Startup of distillation columns is a very challenging control and simulation problem because of both theoretical and practical aspects. Several researchers have studied this subject.14–19 A general sequence of actions which forms the basis for different startup procedures was formulated by Ruiz et al.14 At the end of several actions, all plates have enough liquid holdup so that the liquid can start to fall down the downcomers. The downcomers are sealed, and no vapor can go up through them. The liquid composition is the same on all plates, being equal with the feed composition. In the frame of the present work, these conditions define the initial state from where the effective startup transient operating regime procedure begins. Traditionally, the column is operated at constant values of control parameters (reflux ratio, reboiler heat duty, etc). Usually, these are the prescribed values for the desired stationary regime. In an optimal transient regime procedure, the column will be operated at prescribed time-distributed values of control parameters to minimize the duration of the transient regime. To establish the optimal time-distributed values of control parameters besides a TOC algorithm, an adequate dynamic distillation model (DDM) is needed. This model corresponds to the specific formulation of the general eqs 1 and 2. Dynamic Distillation Model The dynamic simulation of the distillation columns was very extensively studied. Hundreds of papers have been written on this subject, and it now can be considered almost closed. Therefore, this is not the place nor is it the aim of this paper to make an incursion into this field. Nevertheless, some considerations regarding DDM for TOC startup are required. During the application of TOC algorithm based on IDP, the DDM will be solved ten thousand times. Therefore, the model must be a low computer-time consumer. But, if the model is
too simplified, the results will be insignificant. In a rigorous DDM at each time stage representing an integration step of differential mass and energy balances, the calculation of temperature and vapor composition on the column plates is made by an iterative procedure solving algebraic nonlinear equations. For optimization purposes, because of computer time reasons, these models are not suitable. The DDM proposed by the author20,21 represents a good compromise between the degree of complexity and correctness. Because these references are published in Romanian language, the corresponding DDM is detailed in the following paragraphs. The advantage and originality of this model result from the iterative algebraic equations being avoided. The following simplification assumptions are present in this model: (i) The molar vapor holdup is negligible compared to the molar liquid holdup. (ii) Interphase heat transfer is considered much more intense than interphase mass transport; consequently, the liquid and vapor leaving each plate are in thermal equilibrium at the boiling temperature, corresponding to the liquid composition. (iii) On each plate the liquid and vapor are perfectly mixed and Murphree plate efficiency is applied. (iv) Entrainment and weeping rates, flooding of the plates, downcomer holdup, and delay time between plates are neglected. With exception of the fourth assumption, these hypotheses are present in the most of rigorous DDM procedures. For a typical plate in a distillation column, neglect of the molar vapor holdup, mass and energy balance equations is thus defined: Component mass balance around plate j for component i d(Njxi,j) ) Lj-1xi,j-1 + Vj+1yi,j+1 - Ljxi,j dt Vjyi,j ( FL,jxF,i,j ( FV,jyF,i,j
(8)
Total mass balance around plate j dNj ) Lj-1 + Vj+1 - Lj - Vj ( FL,j ( FV,j dt
(9)
Total energy balance around plate j d(Njhj) ) Lj-1hj-1 + Vj+1Hj+1 - Ljhj dt VjHj ( FL,jhF,j ( FV,jHF,j - qj
(10)
The left side of eq 8 is developed to dxi,j dNj d(Ni,jxi,j) ) Nj + xi,j (11) dt dt dt and after substitution of eqs 9 and 11 in eq 8 and rearrangement, the applied form of component mass balance around plate j for component i is dxi,j 1 ) [Lj-1(xi,j-1 - xi,j) + Vj+1(yi,j+1 - xi,j) ( FL,j(xF,i,j dt Nj xi,j) ( FV,j(yF,i,j - xi,j)] (12) The original feature of the model (resulting from the iterative solving of nonlinear algebraic equations is avoided) consists of the calculation of the temperature on each plate. Starting from the constraints m
∑y
i,j ) 1
i)1
that is
(13)
4160 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 m
∑κ x
ij ij ) 1
ln γ1,j )
(14)
i)1
The differentiation with respect to time of eq 14 gives
∑( m
i)1
)
dκi,j dxi,j + κi,j xi,j )0 dt dt
ln γ2,j )
(15)
For ideal and azeotropic mixtures
[
(16)
dγ1,j )dx1,j
that is κi,j ) κi,j(x1,j, x2,j, ..., xm,j, Tj, pj)
pj ) pj+1 - ∆pj
A21
(17)
dγ2,j )dx2,j
(18)
where the pressure drop is calculated on the base of hydraulic correlations, specific for the plate type. With this assumption, the differentiation with respect to time of eq 16 gives dκi,j Pi,j dγi,j γi,j dPi,j ) + dt pj dt pj dt
(19)
It can be written that dγi,j dγi,j dxi,j ) dt dxi,j dt
(20)
Vj )
dPi,j dPi,j dTj ) dt dTj dt
(21)
Substitution of eq 19 into 15 and then substituting eqs 16, 20, and 21, after rearranging, an original equation for dynamic calculation of temperature on plate j results
dTj i)1 )dt
(
1+
m
∑ i)1
)
xi,j dγi,j dxi,j γi,j dxi,j dt
xi,jγi,j dPi,j pj dTj
[
4.606γ1,jA212 3 A12x1,j 1+ (1 - x1,j)2 A21(1 - x1,j)
]
(
)
(22)
[
(26b)
dγ2,j )0 dx1,j
(26c)
4.606γ21,jA221 3 A21x2,j 1+ (1 - x2,j)2 A12(1 - x2,j)
]
(26d)
(
1 L h + Vj+1Hj+1 - Ljhj ( FL,jhF,j ( FV,jHF,j - qj Hj j-1 j-1 dNj dhj dTj - hj (27) dTj dt dt
)
The liquid flow rate is obtained on the base of Francis’ correlation for a plate weir22 (28)
where the height of the liquid over weir is calculated as the difference between the average height of liquid on the plate and the height of the weir Zj )
WL,j - zW Ap
(29)
with m
Nj WL,j )
(23)
∑Mx
i i,j
i)1
(30)
εL,jFL,j
After replacement of the volumetric flow rate with molar flow rate and substitution of eqs 29 and 30 in eq 28, the resulting equation for molar flow rate is
then Pi,jBi dPi,j ) dTj (Tj + Ci)2
A12
(26a)
dγ1,j )0 dx2,j
QL,j ) 1.84Z1.5 j
In eq 22, dxi,j/dt is replaced from eq 12, while dPi,j/dT and dγi,j/ dxi,j are computed according with the selected equations for Pi,j and γi,j. If the vapor pressure of component i on pate j is expressed by Antoine equation Bi Pi,j ) exp Ai Tj + Ci
]
(25b)
2
The core of the model consists of the system of ordinary differential eqs 9, 12, and 22 that represent the particularization of the general continuous dynamic system (eq 1), the state variables being Nj, j ) 1,2, ..., n; xi,j, i ) 1, 2, ..., m - 1, j ) 1, 2, ..., n; and Tj, j ) 1, 2, ..., n. In the mass and energy balance equations, some variables, such as reflux ratio or reboiler duty, can be assigned as control variables u. The rest of the variables involved in these equations are explicitly calculated as follows: Vapor flow rate is obtained from total energy balance (eq 10)
and
i,j
A21 A21x2,j 1+ A12(1 - x2,j)
[
The variation of total pressure during one time integration step is much smaller than the variations of composition and temperature. To simplify the procedure, the pressure is considered constant along the time integration step, but it will be recomputed at the beginning of each new time step
m
]
(25a)
2
results
κij ) γijPij ⁄ pj
∑κ
A12 A12x1,j 1+ A21(1 - x1,j)
(24)
The derivatives dγi,j/dxi,j obviously depend on the activity coefficients expressions. For an ideal mixture, all γi,j are unity, and the derivatives vanish. For a binary azeotropic mixture, using van Laar equations
Lj )
1.84lwFL,j m
∑Mx
i i,j
i)1
(
m
Nj
∑
Mixi,j
i)1
εL,jApFL,j
- zW
)
1.5
(31)
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4161
The vapor composition according to Murphree efficiency is yi,j ) yi,j+1 + Ei,j(y/i,j - yi,j+1)
(32)
Ei,j ) 1 - e-[(KG,i,jAa)/Vj]
(33)
Table 1. General Data for Example 1
where
and KG,i,j )
1 1
+
kG,i,j
(34)
κi,j kL,i,j
The mass transfer coefficients in vapor film kG,i,j and in liquid film kL,i,j on the plate j are calculated using the following criterial equations23 0.9 0.25 ShG,i,j ) 2ReG,j ScG,i,j
(35)
1.1 0.5 0.24 -1 ScL,i,jGaL,j Γ ShL,i,j ) 0.23ReL,j
(36)
Replacement of the criteria with their classical expressions results in kG,i,j ) 95013
0.75 0.65 0.05 0.9 DG,i,j FG,j FL,j wG,j 0.05 0.65 (Tj + 273)σj µG,j
(37)
and kL,i,j ) 1423.3
0.5 0.82 2.08 1.1 DL,i,j zst,j FL,j wL,j
(∑ ) m
(38)
1.08 xi,jMi εL,jµL,j
i)1
The physical properties corresponding to the vapor phase DG,i,j, µG,j, and FG,j are calculated by usual relations in function of vapor composition y1,j, y2,j, ..., ym,j, temperature Tj, and total pressure pj on plate j. Similarly, the physical properties corresponding to the liquid phase DL,i,j, µL,j, FL,j, and σL,j are calculated in function of liquid composition x1,j, x2,j, ... xm,j, and temperature Tj . From hydrodynamic considerations, the static height of the liquid on the sieve tray j is calculated with the relation
( ) m
Lj
zst,j ) 0.65zW + 0.8
∑x
0.667
i,jMi
i)1
lWFL,j
(39)
In eq 37, the vapor velocity is related to the entire cross section of the column m
Vj wG,j )
∑y
i,jMi
i)1
FG,j(πd2c ⁄ 4)
(40)
and in eq 38, the liquid velocity is related to active area per plate m
Lj wL,j )
∑x
i,jMi
i)1
FL,jAa
(41)
The equilibrium, thermodynamic data, and other physical properties correlations are selected according to the mixture nature. The system of differential equations was numerically integrated by fourth-order Runge-Kutta-Gill method. Because of the strong nonlinearity of eq 22, a small initial integration step
a
data
value
number of plates feed plate position plate diameter (m) weir length (m) weir height (m) hole diameter (m) holes area (m2) type of condenser type of reboiler condenser pressure (mm Hg) bottom liquid volume (m3) feed specification feed composition (mole fractions) feed flow rate (kmol s-1) reflux ratioa reboiler heat dutya (kW)
10 4 0.8 0.48 0.025 0.002 0.04 total total 760 0.146 liquid at bubble point benzene 0.6, toluene 0.3, ethylbenzene 0.1 0.0053 3 417
Prescribed values for the desired stationary regime.
Table 2. Initial Condition and Desired State for Example 1 initial condition tray j x1,j x2,j x3,j Tj pj 1 2 3 4 5 6 7 8 9 10 B
0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6
0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
90.81 90.94 91.07 91.19 91.32 91.45 91.58 91.71 91.83 91.96 92.08
770.7 773.7 776.6 779.6 782.5 785.5 788.4 791.3 794.2 797.1 800.0
desired state x1,j
x2,j
x3,j
Tj
pj
0.8887 0.8197 0.7314 0.6192 0.5313 0.4250 0.3127 0.2109 0.1310 0.0755 0.0324
0.1075 0.1719 0.2513 0.3308 0.4137 0.5123 0.6128 0.6969 0.7500 0.7653 0.7026
0.0038 0.0084 0.0173 0.0500 0.0550 0.0627 0.0745 0.0922 0.1190 0.1592 0.2650
82.87 84.56 85.86 90.16 92.79 96.23 100.27 104.43 108.26 111.55 116.04
762.2 765.8 769.5 773.1 776.9 780.7 784.6 788.4 792.2 796.1 800.0
must be used, 3 s. In fact the reliability of the entire numerical procedure depends only on the value of this integration step. We appreciate that the corresponding increasing of the computer time resulting from the small value of the integration step is justified by avoiding the iterative solution of algebraic equations, which are more computer time consumptive. The parameters of differential equations which depend on the integrated variables are calculated from algebraic equations. These parameters are expressed explicitly in algebraic equations and are computed at each integration step. To avoid iterative calculations, the vapor flow rates (eq 27) and vapor composition (eq 32) are computed in the order j ) n, n - 1, ..., 1, while the liquid flow rates (eq 31) are calculated in the order j ) 1, 2, ..., n. This DDM was favorably verified by several theoretical and experimental24 tests. For example 1, the final stationary state given by this model (considering theoretical plates) was confronted with steady-state results of several commercial simulators. The average percentage relative errors of all state variables between the DDM model and the results given by the following simulators are AspenMAX 1.42%, HYSYS 2.31%, PROII 1.12%, and CHEMCAD 2.06%. A short presentation of the experimental test is given in Appendix 1. The computer programs were coded in FORTRAN. The applications were executed on a notebook computer ACER TravelMate 8102WLMi with Intel Pentium M740 (1.73 GHz) processor, 512 MB DDR2 533 memory, and Microsoft Windows XP Professional SP2 operating system. Example 1: A Ternary Ideal Mixture. A mixture of benzene, toluene, and ethylbenzene is distilled in a sieve plates column with lateral downcomers. The general data of this example are given in Table 1, and the initial condition and desired state are indicated in Table 2.
4162 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008
Figure 1. Optimal reflux control for example 1 case a.
Figure 4. Temperature evolutions in time on plate 1 for example 1 case a.
Figure 2. Optimal reboiler heat duty control for example 1 case a.
Figure 5. Temperature evolutions in time on plate 10 for example 1 case a. Table 3. Startup Times and Heat Consumption for Example 1 control policy
Figure 3. Optimal reflux (---) and reboiler heat duty (s) control for example 1 case a.
The vapor pressures of components were calculated by Antoine equation. a. Piecewise Constant Control. The TOC algorithm was applied for S ) 5 time stages. Three types of optimal control policies were established: optimal reflux control (Figure 1), optimal reboiler heat duty control (Figure 2), and the combination of optimal reflux and reboiler heat duty control (Figure 3). The temperature evolutions in time, on plate 1 (top plate), and on plate 10 (bottom plate), for each of the three optimal control policies and for the “traditional” procedure (mentioned above) are shown in Figure 4 and, respectively, Figure 5. In Figure 4 the final values of temperatures on plate 1 are 82.87 °C for traditional procedure, 82.86 °C for optimal reflux control, 82.69
traditional procedure piecewise constant control (5 time stages) optimal reflux optimal reboiler heat duty optimal reflux and reboiler heat duty piecewise linear control (5 time stages) optimal reflux optimal reboiler heat duty piecewise linear control (10 time stages) optimal reflux optimal reboiler heat duty
startup heat time (min) consumption (kWh) 100
695
45 38 24
313 293 191
61 36
424 286
37 28
257 256
°C for optimal reboiler heat control, and 82.87 °C for combined optimal reflux and reboiler control. The deviation from final stationary state in the case of optimal reboiler heat control corresponds to a relative error of 0.217%. This result can be explained by the nonuniform distribution of the state variables relative deviations in the average penalty term of the objective function (eq 6). An argument of this explication is the fact that all the values of final temperature on plate 10 (Figure 5) are in very good agreement (111.56, 111.63, 111.65, and 111.69 °C).
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4163
Figure 6. Optimal reflux control (S ) 5) for example 1 case b.
Figure 7. Optimal reflux control (S ) 10) for example 1 case b.
The comparison of startup times and heat consumption (Table 3) indicates the TOC benefits. The number of grid points for optimal reflux control and for optimal reboiler heat duty control was 5 for s-grid, 9 for u-grid, and 9 for V-grid. For combined optimal reflux and reboiler heat duty control, due to the increased dimensionality, the number of grid points was 3 for s-grid, 5 for u-grid, and 5 for V-grid. According to the recommendations of Bojkov and Luus,4,5,12 the region contraction factor was set at 0.8. The computer execution time on the mentioned notebook computer was 6 h 54 min for optimal reflux control, 5 h 58 min for optimal reboiler duty control, and 2 h 43 min for combined optimal reflux and reboiler duty control. In the first two cases, the total number of IDP iterations was 20, whereas in the last case, there were only 10, and as mentioned before, in this case, the number of grid points was smaller. b. Piecewise Linear Control. The TOC algorithm was applied for S ) 5 and for S ) 10 time stages. Two types of optimal control policies were established: optimal reflux control (Figures 6 and 7) and optimal reboiler heat duty control (Figures 8 and 9). The temperature evolutions in time, on plate 1 and on plate 10, for each of the two optimal control policies and for the traditional procedure for 10 time stages are shown in Figures 10 and 11. The number of grid points was 9 for u-grid, and 9 for V-grid. As in the case of piecewise constant control the region contraction factor was set at 0.8, and the total number of IDP iterations was 20. The comparison of startup times and heat consumption are given in Table 3. Because a single point is used for s-grid at each time stage, increasing the number of time stages improves the value of control policy. The computer
Figure 8. Optimal reboiler heat duty control (S ) 5) for example 1 case b.
Figure 9. Optimal reboiler heat duty control (S ) 10) for example 1 case b.
Figure 10. Temperature evolutions in time on plate 1 (S ) 10) for example 1 case b.
execution time was 2 h 12 min for optimal reflux control with S ) 5 and 2 h 38 min for S ) 10; 1 h 21 min for optimal reboiler duty control with S ) 5 and 1 h 54 min for S ) 10. For piecewise linear control the computer effort is appreciably smaller because a single point for s-grid is used. The value of penalty coefficient was 106 for all applications of this example. This example could become a standard test problem for TOC of distillation columns. To help the reader interested to test
4164 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008
Figure 12. Optimal reflux control for example 2, case a. Figure 11. Temperature evolutions in time on plate 10 (S ) 10) for example 1 case b. Table 4. Initial Condition and Desired State for Example 2 initial condition
desired state
tray j
x1,j
x2,j
Tj
pj
x1,j
x2,j
Tj
pj
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 B
0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
88.82 88.89 88.97 89.05 89.13 89.21 89.28 89.36 89.44 89.51 89.59 89.67 89.74 89.82 89.89 89.97 90.04 90.12 90.19 90.27 90.34
772.9 775.3 777.6 780.0 782.4 784.8 787.2 789.5 791.9 794.3 796.6 799.0 801.3 803.7 806.1 808.4 810.8 813.2 815.5 817.9 820.2
0.4578 0.4651 0.4733 0.4825 0.4931 0.5050 0.5187 0.5342 0.5520 0.5723 0.6006 0.6239 0.6518 0.6846 0.7216 0.7612 0.8008 0.8376 0.8694 0.8955 0.9551
0.5422 0.5349 0.5267 0.5175 0.5069 0.4950 0.4813 0.4658 0.4480 0.4277 0.3994 0.3761 0.3482 0.3154 0.2784 0.2388 0.1992 0.1624 0.1306 0.1045 0.0449
88.38 88.45 88.53 88.64 88.74 88.83 88.94 89.06 89.20 89.37 89.59 89.82 90.12 90.51 91.03 91.69 92.48 93.35 94.21 95.03 97.17
760.0 762.9 765.8 768.7 771.6 774.5 777.4 780.2 783.1 786.0 788.9 791.9 795.0 798.1 801.2 804.3 807.4 810.6 813.7 817.0 820.2
different optimization methods, more information about this case study is given in Appendix 2. Example 2: A Binary Azeotropic Mixture. A binary azeotropic mixture of n-propanol and water is distilled in a column with 20 sieve plates and lateral downcomers. The feed is introduced on the plate 11 with the mole fraction composition 0.6 n-propanol and 0.4 water. The prescribed control values of the desired stationary regime are 2.8 for reflux ratio and 600 kW for reboiler heat duty. The rest of general data are the same as in the previous example shown in Table 1. The initial condition and the desired state are given in Table 4. The vapor pressures of components were calculated on the basis of Antoine equation, and the activity coefficients were calculated with van Laar equation. a. Piecewise constant control. As in the previous example, the TOC algorithm was applied for S ) 5 time stages. Optimal reflux control (Figure 12) and combined optimal reflux and reboiler heat duty control (Figure 13) were established. The temperature evolutions in time, on plate 20 (bottom plate), for these optimal control policies and for traditional procedure are shown in Figure 14.
Figure 13. Optimal reflux (---) and reboiler heat duty (s) control for example 2, case a.
Figure 14. Temperature evolutions in time on plate 20 for example 2, case a.
The number of grid points and the region contraction factor were the same as in the previous example. The total number of IDP iterations was 10. The comparison of startup times and heat consumption, shown in Table 5, indicates the very favorable effect of optimal control policies. For combined optimal reflux and reboiler heat duty control, the startup time must be equal to or lower than startup time in the case of optimal reflux control. Because of the lower number of grid points used in the case of combined optimal
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4165 Table 5. Startup Times and Heat Consumption for Example 2 control policy traditional procedure piecewise constant control (5 time stages) optimal reflux optimal reflux and reboiler heat duty piecewise linear control (5 time stages) optimal reflux piecewise linear control (10 time stages) optimal reflux
startup heat consumption time (min) (kWh) 210
2100
33 33
329 316
54
537
49
491
control policies, this expected result was not obtained. The computer execution time was 4 h 55 min for optimal reflux control and 9 h 27 min for combined optimal reflux and reboiler duty control. b. Piecewise linear control. The TOC algorithm was applied for 10 time stages to establish the optimal reflux control (Figure 15).The temperature evolutions in time, on plate 20 for this optimal control policy and for traditional procedure are shown in Figure 16. The corresponding startup time and heat consumption are given in Table 5. The number of grid points was 9 for u-grid, and 9 for V-grid. As in the case of piecewise constant control, the region contraction factor was set at 0.8, and the total number of IDP iterations was 10. Because the time stages for piecewise constant control and for piecewise linear control are different, a comparison of corresponding policies is not relevant. But, it is unexpected that piecewise linear control established for a higher number of stages is an inferior result compared to piecewise constant control. A possible explanation of this fact can be the use of a single point for s-grid in the case of piecewise linear control. The computer execution time was 1 h 44 min. The value of penalty coefficient was 5 × 105 for all applications of this example. Example 3: An Industrial Scale Process. A mixture of four components (benzene, toluene, ethylbenzene, and o-xylene) is distillated in an industrial scale sieve plates column with lateral downcomers. The general data for this example are given in Table 6, and the initial condition and the desired state are given in Table 7. The vapor pressures of components were calculated by Antoine equation. The TOC algorithm was applied for S ) 5 time stages to establish the optimal reflux control for piecewise constant case (Figure 17) and for piecewise linear case (Figure 18). The temperature evolutions in time on plate 1 (top plate) and plate 30 (bottom plate) are shown in Figures 19 and 20, respectively. The number of grid points for piecewise constant control was 5 for s-grid, 9 for u-grid, and 9 for V-grid, and for piecewise
Figure 16. Temperature evolutions in time on plate 20 for example 2, case b. Table 6. General Data for Example 3 data
value
number of plates feed plate position plate diameter (m) weir length (m) weir height (m) hole diameter (m) holes area (m2) type of condenser type of reboiler condenser pressure (mm Hg) bottom liquid volume (m3) feed specification feed composition (mole fractions) feed flow rate (kmol s-1) reflux ratio a reboiler heat duty a (kW)
30 12 3.5 2.1 0.025 0.005 0.862 total total 760 2.79 liquid at bubble point benzene 0.6; toluene 0.2; ethylbenzene 0.1; o-xylene 0.1 0.03472 2 1251
a
Prescribed values for the desired stationary regime.
linear control, it was 9 for u-grid and 9 for V-grid. For both policies, the region contraction factor was set at 0.8, the total number of IDP iterations was 10, and the value of penalty coefficient was 106. The comparison of startup times and heat consumption (Table 8) indicates the advantages of both control policies. Because of the high dimension of the system of differential equations, the computer effort was appreciable: the execution time was 19 h 1 min for piecewise constant control and 4 h 3 min for piecewise linear control. This example demonstrates the ability of the proposed TOC procedure to be used for industrial scale processes. Discussion
Figure 15. Optimal reflux control for example 2, case b.
The selected value for the number of time stages S ) 5 seems to be too small. A higher value of this number strongly increases the CPU time without substantially improving the performance index. This fact is sustained by several reference examples.9,25 Unlike the case of reflux control when the switching of the reflux values can be made instantaneously, the modification of the reboiler heat duty is made with a time delay, and the startup time and state variables will be affected. These effects are not included in the above applications.
4166 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 Table 7. Initial Condition and Desired State for Example 3 initial condition
desired state
tray j
x1,j
x2,j
x3,j
x4,j
Tj
pj
x1,j
x2,j
x3,j
x4,j
Tj
pj
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 B
0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6
0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
92.81 92.85 92.89 92.94 92.99 93.04 93.08 93.13 93.17 93.22 93.27 93.31 93.36 93.41 93.45 93.50 93.54 93.59 93.64 93.68 93.72 93.77 93.82 93.86 93.91 93.96 94.00 94.04 94.09 94.14 94.18
788.1 789.2 790.2 791.3 792.4 793.4 794.5 795.6 796.6 797.7 798.8 799.8 800.9 801.9 803.0 804.1 805.1 806.2 807.3 808.3 809.4 810.4 811.5 812.6 813.6 814.7 815.8 816.8 817.9 818.8 820.0
0.9975 0.9957 0.9931 0.9893 0.9838 0.9756 0.9636 0.9461 0.9204 0.8826 0.8268 0.6955 0.6956 0.6956 0.6956 0.6956 0.6956 0.6956 0.6956 0.6955 0.6953 0.6949 0.6940 0.6924 0.6892 0.6832 0.6718 0.6505 0.6110 0.5406 0.3540
0.0025 0.0043 0.0068 0.0105 0.0159 0.0237 0.0348 0.0504 0.0715 0.0991 0.1323 0.1714 0.1714 0.1714 0.1714 0.1714 0.1714 0.1714 0.1714 0.1715 0.1717 0.1721 0.1729 0.1743 0.1770 0.1820 0.1910 0.2070 0.2339 0.2750 0.3229
0.0000 0.0000 0.0001 0.0001 0.0002 0.0005 0.0010 0.0023 0.0051 0.0111 0.0236 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0678 0.0679 0.0680 0.0683 0.0689 0.0703 0.0734 0.0805 0.0963 0.1616
0.0000 0.0000 0.0000 0.0001 0.0001 0.0002 0.0006 0.0012 0.0030 0.0072 0.0173 0.0653 0.0652 0.0652 0.0652 0.0652 0.0652 0.0652 0.0652 0.0652 0.0652 0.0652 0.0652 0.0653 0.0655 0.0659 0.0669 0.0691 0.0746 0.0881 0.1615
81.22 81.23 81.29 81.40 81.55 81.76 82.07 82.50 83.14 84.11 85.64 89.87 89.92 89.97 90.01 90.06 90.11 90.16 90.21 90.26 90.31 90.38 90.45 90.54 90.68 90.90 91.28 91.96 93.25 95.71 104.07
763.9 765.6 767.2 768.8 770.4 772.0 773.6 775.2 776.8 778.4 780.0 781.6 783.6 785.7 787.7 789.7 791.7 793.7 795.8 797.8 799.8 801.8 803.9 805.9 807.9 809.9 811.9 814.0 816.0 818.0 820.0
The improvement obtained by TOC of startup distillation shown in Tables 3, 5, and 8 in comparison with the traditional startup procedure seems to be artificially amplified, but this traditional startup procedure is frequently used. Of course, it is possible that the comparison with other particular (empirically improved) startup procedures to give lower differences in performance. Because of bang-bang control (especially for piecewise constant control), the practical implementation of the optimal policies obtained here can lead to some hydrodynamic problems (flooding, weeping, liquid entrainment, etc). Therefore, it is useful to test these control policies by simulations, using a DDM with more detailed plate hydraulics. In this way, some suboptimal control policies, avoiding wrong hydrodynamic regimes, can be found by suitable corrections. Of course, a very useful, but a little bit difficult way to validate these control policies, based exclusively on simulation, is to test the approach on a pilot plant.
The DDM used here was applied for ideal and azeotropic mixtures. For nonideal mixtures (distillation at low temperatures, high pressures, or both), the calculation of equilibrium constants are made with implicit methods (e.g., Soave, Peng-Robinson, etc.). In this case, the application of numerical differentiation involved in the plate temperature equation (eq 22) is expected to fail. A more promising way to extend the use of present DDM for this case is to obtain previously correlation relations for equilibrium constants and to derive them analytically with respect to temperature and composition.
Figure 17. Optimal reflux control (piecewise constant) for example 3.
Figure 18. Optimal reflux control (piecewise linear) for example 3.
Conclusions In this work, the ability of IDP to solve high dimensional TOC problems, involving complex models, was demonstrated. Optimal control policies for startup distillation columns were obtained on the base of a DDM representing a good compromise between correctness and complexity. This result was possible mainly because of the transformation of the differential algebraic equations model into an ordinary differential equations model, which is
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4167
Figure 19. Temperature evolutions in time on plate 1 for example 3.
Figure 20. Temperature evolutions in time on plate 30 for example 3. Table 8. Startup Times and Heat Consumption for Example 3
control policy traditional procedure piecewise constant control (5 time stages) optimal reflux piecewise linear control (5 time stages) optimal reflux
In a first stage, the distillation was conducted at total reflux and constant heating of the still to stationary regime, respectively, time-constant values for top and bottom temperatures and pressures. Some experimental values obtained in stationary conditions were used in the dynamic simulation. These values are still heat flux 22 W, average plate efficiency 0.424, and pressure drop per a tray 0.5 mmHg. In a second stage, a step signal for the bottom composition was generated by a very quick injection of one mole of toluene into the still. The temperature of injected toluene was the boiling temperature of the new bottom composition. The variations of top and bottom temperatures were continuously registered up to their constant values, at total reflux and constant bottom heating. The top temperature variation during the transient regime is insignificant. The explanation consists of the total separation of the light component because of the high number of trays in the laboratory column. The simulated dynamic response of bottom temperature according to bottom composition modification is very close to the experimental dynamic response (Figure A1).
Figure A1. The simulated (s) and experimental (---) dynamic response of bottom temperature.
startup time (min)
heat consumption (kWh)
138
2877
57
1188
Component index designation: 1, benzene; 2, toluene; 3, ethylbenzene. a. Equilibrium Data. Antoine equation with
48
1001
A1 ) 5.881966 A2 ) 16.010667 A3 ) 16.011394
Appendix 2
B1 ) 2777.724 B2 ) 3094.543 B3 ) 3274.078
computationally much more efficient. For three illustrative examples, the startup times and heat consumption were reduced substantially, providing an appreciable economic benefit. Acknowledgment The research presented here has no kind of financial support. This paper is dedicated to Professor Rein Luus for his brilliant works in the field of iterative dynamic programming.
C1 ) 220.237 C2 ) 219.337 C3 ) 212.931 b. Thermodynamic Data. Liquid enthalpy of components on plate j (A1)
cp,i,j ) ai + biTj
(A2)
where
Appendix 1 A simple experimental test was conducted in a laboratory bubble cap trays column. The column characteristics are as follows: well insulated, 17 trays, and 0.023 m internal diameter. The top and the bottom temperatures and pressures were measured online. A binary benzene-toluene mixture prepared from pure components was distilled. The initial mole ratio was 1:1.
hi,j ) cp,i,jTj
with a1 ) 130897.1 a2 ) 153005 a3 ) 173104.5 b1 ) 326.5 b2 ) 346.6 b3 ) 406.9 Saturated vapor enthalpy of components on plate j
4168 Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008
Hi,j ) hi,j + rV,i,j
DL,k,i,j ) DL,i,k,j
(A3)
(A14)
Diffusion coefficient in liquid phase of component i on plate j
where rV,i,j ) aV,i + bV,iTj
DL,i.j )
with
∑ i)1
aV,1 ) 35752626 aV,2 ) 38954079 aV,3 ) 52791147 bV,1 ) -62063.6 bV,2 ) -51219.9 bV,3 ) -127510.5
µG,i,j ) µ0,i
3
∑x
i,jhi,j
(A5)
(A15)
3
xk,j D k)1;k*i L,i,k,j
∑
Vapor viscosity of component i on plate j
Liquid mixture enthalpy on plate j hj )
1 - xi,j
3
(A4)
(
273 + Cµ,i Tj + 273 Tj + 273 + Cµ,i 273
)
1.5
(A16)
with
i)1
µ0,1 ) 0.68 × 10-5 µ0,2 ) 0.66 × 10-5 µ0,3 ) 0.59 × 10-5
Saturated vapor mixture enthalpy on plate j
Cµ1 ) 281 Cµ2 ) 269 Cµ3 ) 323
3
Hj )
∑
yi,jHi,j
(A6)
i)1
( )
Vapor mixture viscosity on plate j 3
c. Hydrodynamics. Pressure drop on sieve tray j 3
∆pj ) -
0.00683 FG,j
QG,j
∑y
i,jMi
µG,j ) σj - 0.02175 - 0.07352FL,jzst,j do (A7)
i)1
Ao
∑ i)1
(A17)
yi,jMi µG,i,j
log µL,i,j ) aµi +
(A8)
bµi Tj + 273
(A18)
with
and
aµ1 ) -5.096225 aµ2 ) -4.90586 aµ3 ) -4.799094
3
σj )
∑
xi,jσi,j
(A9)
bµ1 ) 563.8389 bµ2 ) 498.3618 bµ3 ) 478.5327
i)1
Liquid mixture viscosity on plate j
with
(
σ0,1 ) 0.03108 σ0,2 ) 0.0299 σ0,3 ) 0.0315
µL,j ) exp
R0,1 ) 0.00394 R0,2 ) 0.00351 R0,3 ) 0.0021 d. Other Physical Data. Diffusion coefficient in vapor phase of component i into component k on plate j -7
DG,i,k,j )
4.3 · 10 (Tj + 273)1.5 pj 0.33 (υ + υ0.33 k ) 760 i
1 1 + Mi Mk
DG,k,i,j ) DG,i,k,j
i)1
3
yk,j D k)1;k*i G,i,k,j
∑
M0.5 k (Tj + 273) µL,k,jυ0.6 i
i,jMi
Tj + 273
FL,i,j ) F0,i(1 - βiTj) (A11)
(A20)
(A21)
with F0,1 ) 878.6 F0,2 ) 866 F0,3 ) 866.9 β1 ) 0.00106 β2 ) 0.001109 β3 ) 0.00096 Liquid density on plate j 3
FL,j )
∑ i)1
(A12)
xi,jMi 3
∑x
FL,i,j
i,jMi
i)1
Diffusion coefficient in liquid phase of component i into component k on plate j DL,i,k,j ) 7.4 × 10-14
(A19)
∑y i)1
FG,j )
Diffusion coefficient in vapor phase of component i on plate j
∑
)
Liquid density of component i on tray j
M1 ) 78.11 M2 ) 92.13 M3 ) 106.16
DG,i.j )
ln µL,i,j
3
(A10)
υ1 ) 96.118 υ2 ) 118.2 υ3 ) 140.4
1 - yi,j
i,j
i)1
0.01603pj
with
3
3
∑x
Vapor density on plate j
and
and
i)1 3
Liquid viscosity of component i on plate j
where σi,j ) σ0,i(1 - RiTj)
∑y
i,jMi
2
(A13)
Nomenclature Aa ) active area per plate (m2) A0 ) holes area per plate (m2) Ap ) total area per plate (m2) A, B, C ) Antoine equation constants A12, A21 ) van Laar equations constants D ) diffusion coefficient (m2 s-1)
(A22)
Ind. Eng. Chem. Res., Vol. 47, No. 12, 2008 4169 dc ) column diameter (m) d0 ) hole diameter (m) E ) Murphree efficiency F ) feed or output side streamflow rate (kmol s-1) Ga ) Galileu number H ) vapor enthalpy (J kmol-1) h ) liquid enthalpy (J kmol-1) I ) performance index K ) overall mass transfer coefficient (kmol m-2 s-1) k ) film mass transfer coefficient (kmol m-2 s-1) L ) liquid flow rate (kmol s-1) lW ) weir length (m) M ) molecular weight (kg kmol-1) m ) number of components N ) liquid holdup per plate (kmol) n ) number of plates ns ) dimension of state vector nu ) dimension of control vector P ) vapor pressure (mmHg) p ) total pressure on plate (mmHg) Q ) volumetric flow rate (m3 s-1) q ) heat rate losses per plate (W) R ) reflux ratio Re ) Reynolds number S ) number of time stages s ) state vector Sc ) Schmidt number Sh ) Sherwood number T ) temperature (°C) t )time (s) u ) control vector V ) vapor flow rate (kmol s-1) V ) time stage (s) W )liquid volume on plate (m3) w ) velocity (m s-1) x ) liquid mole fraction on plate y ) vapor mole fraction on plate Z ) height of the liquid over weir (m) zst ) static height of the liquid (m) zW ) weir height (m) Greek Letters Γ) geometric simplex γ ) activity coefficient ∆p ) pressure drop (mmHg) ε ) volumetric fraction κ ) equilibrium constant µ ) viscosity (kg m-1 s-1) F ) density (kg m-3) σ ) superficial tension (kg s-2) τ ) normalized time variable ω ) penalty coefficient Subscripts F ) feed f ) final G )gas i ) component j ) plate k ) time stage L ) liquid s ) stationary V ) vapor
Superscript * ) equilibrium
Literature Cited (1) Lapidus, L.; Luus, R. Optimal Control of Engineering Processes; Blaisdel: Waltham, MA, 1967. (2) Bashein, G. A Simplex Algorithm for On-line Computational of Time-Optimal Controls. IEEE Trans. Autom. Control 1971, AC-16, 479. (3) Rosen, O.; Imanudin, I.; Luus, R. Sensitivity of the Final State Specification for Time-Optimal Control Problems. Int. J. Control 1987, 45, 1371. (4) Bojkov, B.; Luus, R. Application of Iterative Dynamic Programming to Time Optimal Control. Chem. Eng. Res. Des. 1994, 72, 72. (5) Bojkov, B.; Luus, R. Time-Optimal Control by Iterative Dynamic Programming. Ind. Eng. Chem. Res. 1994, 33, 1486. (6) Bojkov, B.; Luus, R. Time Optimal Control of High Dimensional Systems by Iterative Dynamic Programming. 1995. Can. J. Chem. Eng. 1995, 73, 380. (7) Dadebo, S. A.; McAuley, K. B. A Simultaneous Iterative Solution Technique for Time-Optimal Control Using Dynamic Programming. Ind. Eng. Chem. Res. 1995, 34, 2077. (8) Luus, R.; Okongwu, O. N. Towards Practical Optimal Control of Batch Reactors. Chem. Eng. J. 1999, 75, 1. (9) Luus, R. Optimal Control by Dynamic Programming Using Accessible Grid Points and Region Reduction. Hung. J. Ind. Chem. 1989, 17, 523. (10) Fikar, M.; Latifi, A. M.; Fournier, F. ; Creff, Y. Application of Iterative Dynamic Programming to Optimal Control of a Distillation Column. Can. J. Chem. Eng. 1998, 76, 1110. (11) Bojkov, B.; Luus, R. Use of Random Admissible Values for Control in Iterative Dynamic Programming. Ind. Eng. Chem. Res. 1992, 31, 1308. (12) Bojkov, B.; Luus, R. Evaluation of the Parameters Used in Iterative Dynamic Programming. Can. J. Chem. Eng. 1993, 71, 451. (13) Luus, R. Piecewise Liniar Continuous Control by Iterative Dynamic Programming. Ind. Eng. Chem. Res. 1993, 32, 859. (14) Ruiz, C. A.; Cameron, I. T.; Gani, R. A Generalized Model for Distillation Columns III. Study of Startup Operations. Comput. Chem. Eng. 1988, 12, 1. (15) Yasuoka, H.; Nakanishi, E.; Kunikita, E. Design of an On-line Startup System for Distillation Column Based on a Simple Algorithm. Int. Chem. Eng. 1987, 27, 466. (16) Barolo, M.; Guaries, G. B.; Rienzi, S.; Trotta, A. On-line Startup of a Distillation Column Using Generic Model Control. Comput. Chem. Eng. 1993, 17 (S), 349. (17) Barolo, M.; Guaries, G. B.; Rienzi, S.; Trotta, A. Nonlinear ModelBased Startup and Operational Control of a Distillation Column: An Experimental Study. Ind. Eng. Chem. Res. 1994, 33, 3160. (18) Fieg, G.; Wozny, G.; Kruse, C. Experimental and Theoretical Studies of the Dynamics of Startup and Product Switchover Operations of Distillation Columns. Chem. Eng. Process. 1993, 32, 283. (19) Han, M.; Park, S. Startup of Distillation Columns Using Profile Position Control Based on a Nonlinear Wave Model. Ind. Eng. Chem. Res. 1999, 38, 1565. (20) Woinaroschy, A. A New Model for the Dynamic Simulation of the Rectification Processes. I. Development of the Calculation Model and Algorithm. ReV. Chim. 1986, 37, 697. (21) Woinaroschy, A. A New Model for the Dynamic Simulation of the Rectification Processes. V. Simulation of Non-ideal Mixtures Separation. ReV. Chim. 1987, 38, 701. (22) Perry, R. H., Ed. Perry’s Chemical Engineering Handbook., 7th ed; Mc Graw-Hill: New-York, 1999. (23) Kafarov, V. Fundamentals of Mass Transfer; Mir Publishers: Moscow, 1975. (24) Woinaroschy, A. A New Model for the Dynamic Simulation of the Rectification Processes. II. Theoretical and Experimental Testing of the Model. ReV. Chim. 1986, 37, 879. (25) Luus, R.; Galli, M. Multiplicity of Solutions in Using Dynamic Programming for Optimal Control. Hung. J. Ind. Chem. 1991, 19, 55.
ReceiVed for reView October 13, 2007 ReVised manuscript receiVed March 4, 2008 Accepted March 5, 2008 IE0713745
