Ind. Eng. Chem. Res. 1998, 37, 1341-1350
1341
PROCESS DESIGN AND CONTROL Optimization of a Semibatch Distillation Process with Model Validation on the Industrial Site Pu Li, H. Arellano Garcia, and Gu 1 nter Wozny* Institut fu¨ r Prozess-und Anlagentechnik, Technische Universita¨ t Berlin, Skr. KWT-9, 10623 Berlin, Germany
Erich Reuter Abt. CPV, Henkel KGaA, 40191 Du¨ sseldorf, Germany
An industrial semibatch distillation process with a chemical reaction in the reboiler is optimized to minimize the batch operation time. A detailed dynamic model is employed to describe the process and incorporated in the optimization problem. This model is validated with the measured data from the experiment directly on the industrial site. The optimal reflux ratio and the feed rate policies of the process are developed with an efficient optimization approach. It can be shown that up to 30% of the operation time can be saved with the optimal policies developed in comparison to the conventional operation. The modeling and optimization approach proposed can be readily applied to the optimization of general batch distillation processes. 1. Introduction Batch distillation processes become more and more attractive to the chemical, food, and pharmaceutical industries because of their operational flexibility. A charge of a multicomponent mixture can be separated to achieve multiple products (fractions) by operating a single column. In particular, the amounts as well as the kinds of both the feed charge and the product of an industrial process may be frequently changed corresponding to the changes of the market conditions. Large-scale continuous processes may not endure such drastic changes due to the fact that the redesign work required and the operation costs will be much larger than those of batch processes. Therefore, more and more attention has been directed to the studies of batch distillation in the past few years. Previous and recent studies of batch distillation have focused on simulation (Ahmad and Barton, 1996; Skogestad et al., 1997; Wajge et al., 1997; Li et al., 1997a) and optimization (Logsdon and Biegler, 1989; Fahart et al., 1990; Diwekar and Madhavan, 1991; Mujtaba and Macchietto, 1992a) of such processes. To facilitate both the distillation and the reaction, reactive batch distillation has been proposed and investigated (Albet et al., 1991; Mujtaba and Macchietto, 1992b; Sørensen et al., 1996), in which usually a reversible reaction, which takes place in the reboiler, is considered. Semibatch distillation has also been studied because in some reactive and azeotropic batch distillation operations it is necessary to introduce a feed into the reboiler during the batch (Reuter et al., 1989; Dussel and Stichlmair, 1994). The rigorous modeling of batch distillation processes * To whom correspondence is to be addressed. Fax: 004930-31426915. Phone: 0049-30-31423893. E-mail: wozny@ dynamik.fbl0.TU-Berlin.DE.
formulates a large-scale nonlinear differential and algebraic equation (DAE) system. Therefore, in the previous studies of batch distillation optimization, simplified models have often been used as a compromise between the model accuracy and the solution possibility of the optimization problem (Li et al., 1997b). In addition, in these studies the reflux ratio profile is optimized either to maximize the operating profit, to maximize the amount of product or to minimize the batch time. For semibatch distillation, the feed policy plays an important role in the process operation and is also to be optimized which, however, leads to a more complicated problem. The optimization of both the reflux ratio and the feed flow rate profiles has not been considered in the previous studies. Besides the above theoretical investigations, a limited number of experimental works have been conducted. Nad and Spiegel (1987) reported their experimental study of a pilot batch column to validate the model through a comparison of the measured and simulated profiles. More recently, Barolo et al. (1996) studied, through experiment, a pilot batch column with a middle vessel, in which some practical issues of the batch operation are provided. However, to date no reports on experimental studies of batch distillation processes optimization of an industrial scale have been found in the open literature. In this study, we consider an industrial semibatch distillation process with a transesterification taking place in the reboiler. An experiment of a conventional batch run of the process was conducted directly on the industrial site. To describe the process more accurately, we use a detailed dynamic model which is validated with the measured data from the experiment. Then, by using a modified optimization approach including successive quadratic programming (SQP) and orthogonal collocation on finite elements, both the reflux ratio and the
S0888-5885(97)00695-7 CCC: $15.00 © 1998 American Chemical Society Published on Web 03/10/1998
1342 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
Figure 1. The semibatch distillation process.
feed rate policies are optimized to minimize the batch time. The approach is especially suitable for large-scale systems because only the independent variables are treated by SQP. The dependent variables are solved through integration of the DAEs and the necessary sensitivities are computed through transmission of the gradients from interval (element) to interval. With the developed optimal policy, the reaction rate will be accelerated and thus about 30% of the batch time can be saved in comparison to that with the conservative conventional operation. 2. Process Description An industrial semibatch distillation process (Figure 1) composed of a total condenser, a column with 30 bubble-cap trays, a reboiler, and two accumulators (for a main cut and an offcut) is considered. A slightly endothermic transesterification of two esters and two alcohols takes place in the reboiler, which can be described as follows:
educt ester (A) + educt alcohol (B) T product ester (C) + product alcohol (D) The pure educt ester, educt alcohol, a mixture of the educt and product alcohols from the off-cut of the last batch, and the homogeneous catalyst will be charged to the reboiler at the beginning of a batch operation. Then, the reboiler will be heated and during the startup period the liquid condensed in the condenser will be totally fed back to the column. An inert gas is introduced to keep the column at a constant operating pressure. During the batch, a limited amount of educt alcohol will be fed to the reboiler to accelerate the reaction rate. The product alcohol (the lightest component) is distillated from the reboiler, through which the reaction will
be shifted toward the product direction. The use of excess educt alcohol and the simultaneous removal of the product alcohol through distillation are highly effective in displacing the reaction equilibrium toward the right and hence to a high reaction rate. In the maincut period, the product alcohol is accumulated in the first distillate accumulator with a given purity specification. In the off-cut period, the reaction temperature will be raised to react the educt ester off to an acceptable extend so as to prevent an additional special separation step. At the end of the batch, a mixture of the product ester and the educt alcohol will be achieved in the reboiler and then separated by a recovery column. The off-cut will be recharged to the reboiler of the next batch. The temperatures, flow rates, and pressures of the different positions shown in Figure 1 can be measured on-line, while the composition measurement is provided by off-line chromatographic analysis. Several control loops are implemented to ensure the stable operation of the process. Since the major costs are the working time of the operators and the consumption of the heating steam, it is desired to reduce the batch time through optimization of the operation policies. For this semibatch distillation process, we consider the profiles of the reflux ratio and the feed flow rate within the batch time as the operation policies to be optimized. The constraints will be the purity specifications of both the product alcohol in the distillate and of the educt ester in the reboiler at the end of the batch as well as the feed amount and flow rate limitations. To formulate a mathematical optimization problem and solve it with a model-based approach, it is imperative to set up a mathematical model that can accurately describe the process. 3. Modeling and Model Validation 3.1. Detailed Modeling of the Process. In this study, a detailed tray-by-tray model is used to describe the batch column. Included are component balance, energy balance and vapor-liquid equilibrium of each tray in the unit model. The holdup, the pressure drop of each tray, and non-equimolar flow in the column are taken into consideration. The vapor-liquid equilibrium is described with nonideal liquid phase computed by the NRTL model. The reaction kinetics are added to the model to depict the chemical reaction in the reboiler. The component compositions of both the liquid and the vapor phases, the liquid and vapor flow rates, and the temperature are the variables of each tray. That is to say, the total number of the variables will be the following: (N-component × 2 + 3)N-tray. Thus, a large complex dynamic system composed of nonlinear DAEs will be formulated. A detailed description of the model is given in Appendix A. The parameters required in the model are the constants of the Antoine component vapor pressure equation, of the component heat capacity correlation, and of the NRTL model which are taken from Reid et al. (1987) and Gmehling et al. (1977), respectively. The kinetic parameters of the transesterification are taken from the experimental study by Reuter (1994). The control (independent) variables include the reflux ratio, the distillate, and the feed flow rate. The profiles of these variables within the batch operation should be given prior to the simulation. It should be noted that the start-up period with total reflux is not included in the model since this period is relatively small compared to
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1343
Figure 2. Measured feed and distillate flow rate profiles.
Figure 3. Measured reflux ratio profile.
the main- and off-cut periods. We therefore consider the process from the starting point of distillation. 3.2. Experiment. To validate the model, several experiments were carried out directly on the industrial site. We show the results of one operation. The conventional operation of the process was performed during the validation experiment. At the beginning, 6.9 m3 of pure educt ester, 2.6 m3 of a mixture of the educt alcohol and the product alcohol of the off-cut from the last batch, and 0.4 m3 of pure educt alcohol were charged to the reboiler, which led to a starting composition (molar fraction) with 〈0.2195, 0.2944, 0.1329, 0.3532〉 for the components A, B, C, and D, respectively. Also, 80 kg of catalyst were charged into the reboiler. The reboiler was then stirred and heated with the pressurized steam. After a half-hour of total reflux, the process was switched to the main-cut period to produce the product alcohol from the distillate. When the composition of the product alcohol in the distillate sank below its specification, the process was switched to the off-cut period. The total batch time of the experiment was about 26 h. During the batch, the empirical polices of the feed, distillate, and reflux flow rates were realized through manipulating the set points of the existing control loops. The level of the reflux drum was kept within its allowable range during the experimental run. It can be noted here the column was not always operated at the maximum vapor load because both the distillate and the reflux flow rate were control variables. Figures 2 and 3 show the measured profiles of the feed and the distillate flow rate as well as the reflux ratio policy. The
Figure 4. Simulated and measured profiles of distillate composition.
samples from the condenser and the reboiler were taken per hour and analyzed to obtain the composition profiles. 3.3. Simulation Results. To efficiently simulate this batch operation, we use collocation on finite elements to discretize the dynamic model equations (i.e., in each time interval (length of the element) the dependent variables are approximated with Lagrange polynomials with shifted roots of Legendre polynomials (Finlayson, 1980)). To keep the continuity of the variables between intervals, we use the last collocation point as the starting point of the next interval. The Newton-Raphson algorithm is then applied to solve the algebraic equations derived from the discretization, so that the operation of the process can be simulated from interval to interval. To take advantage of the dominant block tridiagonal structure of the Jacobian matrix, a Gauss elimination approach for sparse matrices is employed. A FORTRAN program was used to perform the simulation. This simulation program can be directly incorporated into the optimization approach. Prior to the simulation, the pressure, holdup, and tray efficiency parameters of the model must be computed. The pressure of each tray is calculated through interpolation between the top and bottom operating pressures. The holdup of each tray is determined from the tray volume. Since the value of the Murphree tray efficiency is usually in the range between 0.7 and 1.0, it can be determined through trial-and-error by comparing the simulated results with the measured data from the experiment. Figures 4-7 show the comparison of the simulation results to the measured experimental profiles. As desired, the product alcohol (Figure 4) is received during the main-cut period from the distillate. Near the end, the educt ester is almost completely eliminated through the reaction and the product left in the reboiler is a mixture of the product ester and the remaining educt alcohol (Figures 5 and 6). The temperature at the column top shown in Figure 7 remains constant in the main-cut period and rises in the off-cut period with the change in its composition, while the reboiler temperature increases steadily during the batch along with the distillation of the light components. It can be seen from these results that satisfactory agreement between the simulated and the measured profiles have been achieved. This implies that the model considered and the parameters employed have been validated and are suitable to be used in the optimization of the process.
1344 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
distillated from the reboiler due to the higher temperature. Thus, it is not desirable at this time to feed the educt alcohol to the reboiler. 4. Optimization
Figure 5. Simulated alcohol profiles in the reboiler.
Figure 6. Simulated and measured ester profiles in the reboiler.
Figure 7. Simulated and measured temperature profiles.
In addition, we can notice from these results the optimization potential to reduce the batch time. For instance, the reflux ratio in the main-cut period shown in Figure 4 is fairly large, which leads to a conservative operation: the purity of the product alcohol is much greater (average value: 0.993 mol/mol) than its specification (0.98 mol/mol). Since the product alcohol cannot be removed quickly from the reboiler with this reflux ratio profile, the reaction will be lowered. Another fact that points to the inefficiency of the conventional operation is the improper distribution of the limited amount of the feed of the educt alcohol. In the off-cut period, for example, a part of the educt alcohol will be
4.1. The Optimization Approach. Recent developments in dynamic optimization tend to discretize the DAEs within the time period considered, usually with orthogonal collocation on finite elements, and then to simultaneously solve the static optimization problem and integrate the DAEs with a nonlinear programming (NLP) algorithm (e.g., the SQP method (Cuthrell and Biegler, 1987: Renfro et al., 1987; Vassiliadis et al., 1994a,b)). The simultaneous approach has been modified (Logsdon and Biegler, 1992) for large-scale systems by dividing the variable space into a state and a control space. Thus, only the control or independent variables will be treated by SQP. The dependent variables will be solved in the integration of the DAEs with the Newton method and the necessary sensitivities computed from the interval gradient information. Another advantage of the modified approach over the simultaneous approach is that the solution is always feasible, even if the search is interrupted at a suboptimum or near the optimum due to a strong nonlinearity and a nonconvexity of the system. This modified approach has been successfully applied to the optimization of batch distillation for the development of optimal reflux ratio policies (Logsdon and Biegler, 1989; Li et al., 1997b). The modified approach can be structured as a twolayer solution approach. Beginning from an estimated profile for the independent variables, in the simulation layer the process in terms of DAEs will be discretized with orthogonal collocation and simulated using the Newton method. The sensitivities of the dependent variables to the independent variables of each time interval, which are required from the optimization layer, will be transmitted through the continuity relation from interval to interval. In the optimization layer, in which only the inequality constraints are included, the new solution of the optimal profiles of the independent variables will be found using SQP, on the basis of the sensitivities given from the simulation layer. There is a detailed derivation of the approach in Appendix B. 4.2. Problem Formulation. In this study, we use this modified approach to optimize the industrial semibatch distillation process by optimizing both the reflux ratio and the feed flow rate policies. The aim of the optimization is to minimize the batch operation time under the equality constraints of the model equations and the inequality constraints including the product purity specifications and the physical restriction of the feed. It should be noted that the model should be slightly modified for the optimization: the vapor flow from the column to the condenser will be fixed at the maximum vapor load, which means operating at a constant, namely the capacity of the existing equipment. According to the structure of the column considered, the maximum vapor load was computed by its F-factor (Kister, 1992) and is 21.6 kmol/h. Thus, the independent variables of the problem are the feed flow rate and the reflux ratio. The nonlinear dynamic optimization problem can now be described as follows:
min tf(F(t),RV(t),tu,tf) subject to the following equality constraints: component
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1345
Figure 8. Optimal reflux ratio profile (scenario 1). Figure 9. Optimal feed flow rate profile (scenario 1).
balance for each tray and each component; phase equilibrium for each tray and each component; summation equations for each tray; energy balance for each tray; total mass balance for the condenser; and total mass balance for the reboiler. It is also subject to the following inequality constraints: product alcohol specification, xD,1(t) g 0.98 mol/mol; t0 e t e tu; purity specification in the reboiler, xA,NST(tf) e 0.002 mol/mol; feed flow rate restriction. 0 e F(t) e 150 1/h; limitation of the total feed amount:
∫tt F(t) dt e 20 kmol f
0
A three-point collocation was used to discretize the system and divide the total batch time into 100 time intervals. Piecewise constant forms of the independent variables (the reflux ratio and the feed flow rate) are considered. In addition, the lengths of different intervals are also regarded as independent variables in order to conveniently handle the fraction switching time tu and the total batch time tf. As mentioned before, the simulation approach proposed in section 3 can be easily incorporated into the optimization. For the optimization layer, the SQP subroutine from IMSL (1987) was used. The computation of the sensitivities is coded in FORTRAN. 4.3. Results and Analysis. The conventional operation profiles of the reflux ratio and the feed flow rate measured in the validation experiment were chosen as the initial estimates for the optimization. This means that the first iteration of the computation represents the conventional operation of the process. The computation will then proceed to reduce the batch time by improving the profiles of the two independent variables from iteration to iteration until an optimum is reached or at least a local optimum is achieved. For comparison, two scenarios have been computed: in scenario 1, the constraints of the conventional operation will be kept; in scenario 2, all the constraints are the same except that the purity specification will be lowered. Scenario 1. Figure 8 shows the resultant optimal reflux ratio policy. We see that, compared with the conventional profile shown in Figure 3, the reflux ratio is substantially lowered in the main-cut period to remove the product alcohol from the reboiler as far as possible and thus to accelerate the reaction rate. In contrast, the reflux ratio is raised during the off-cut period in order to prevent the distillation of the educt alcohol from the reboiler.
Figure 10. Reaction rate profile (scenario 1).
The policy of the feed flow rate shown in Figure 9 illustrates an optimal distribution of the educt alcohol within the batch, presenting an obvious improvement over the conventional profile shown in Figure 4. That is to say, the educt alcohol should be fed into the reboiler at the maximum flow rate during the main-cut period to “push up” the reaction, while in the off-cut period the feed should not be introduced because a part of the educt alcohol will be distillated due to the rising temperature in the reboiler. However, there are two peaks of the feed rate at the beginning and the end of this period. The first peak is to prevent the sharp reduction of the reaction rate and the second is to eliminate the educt ester left in the reboiler in order to fulfill the strict purity requirement, which can be more clearly illustrated with the reaction rate profile shown in Figure 10. From the distillate composition profiles (Figure 11), we see that the composition of the product alcohol is maintained at the purity specification value during the main-cut period, which means the product alcohol is “pulled out” at the maximum possible rate in order to facilitate the reaction. This is known as a singular optimal control problem or path constraint problem (Logsdon and Biegler, 1989) that may present difficulties to be solved with other optimization methods. With the approach used in this case, satisfactory results have been achieved through the combination of the optimal policies for the reflux ratio and for the feed flow rate. Figures 12 and 13 show the reboiler composition profiles, indicating that the educt components in the reboiler are not dominant within the off-cut period.
1346 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
Figure 11. Distillate composition profiles (scenario 1). Figure 14. Temperature profiles (scenario 1).
Figure 12. Alcohol profiles in the reboiler (scenario 1). Figure 15. Optimal reflux ratio profile (scenario 2).
Figure 13. Ester profiles in the reboiler (scenario 1).
Figure 16. Optimal feed flow rate profile (scenario 2).
Therefore, another possibility to increase the reaction in this period should be considered (i.e., to increase the reaction temperature). Figure 14 illustrates that the temperature in the reboiler is raised as desired in this period by the optimal policies. Scenario 2. The aim of the second scenario is to test the influence of the changes of the product specification on the optimum for the process, namely, to test the flexibility and adaptivity of the optimization approach. We assume in this scenario that the purity of the product alcohol is lowered from 0.98 to 0.96 mol/mol, but all other conditions remain the same as in the first scenario. Figure 15 gives the optimal policy of the reflux ratio for this scenario. Because of the lower purity of the distillate during the main-cut period, the reflux ratio can be decreased even further, compared with that of
scenario 1, shown in Figure 8. Therefore, the product alcohol will be removed more quickly and the reaction rate can be accelerated, so that the batch time can be further reduced. The flow rate profile shown in Figure 16 again illustrates that a great part of the feed should be fed during the main-cut period. But to meet the purity requirement of the educt ester, an amount of feed is needed to react at the end of the batch. The composition profile of the product alcohol (Figure 17) is successfully maintained at the new specification during the main-cut period, which again demonstrates the effectiveness of this optimization approach. Table 1 summarizes the different operating conditions under the different purity constraints. It is shown that over 30% operation time can be saved with the optimal policies compared to the conventional operation. Additionally, if the product specification lowered, which
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1347
Figure 17. Distillate composition profiles (scenario 2). Table 1. Conventional and Optimal Operation of the Process
conventional scenario 1 scenario 2
specification (mol/mol)
switching time (h)
batch time (h)
x1,D(t) ) 0.98 x*NST,A(tf) ) 0.002 x1,D(t) ) 0.98 x*NST,A(tf) ) 0 002 x1,D(t) ) 0.96 x*NST,A(tf) ) 0.002
tu ) 16.15
tf ) 25.52
tu ) 9.22
tf ) 17.39
tu ) 7.28
tf ) 16.16
may occur due to the market changes, the batch time can be further reduced through the optimization of the reflux ratio and the feed rate policies. 5. Conclusions Although many studies of batch distillation optimization have been conducted, the optimization on industrial processes remained a challenging problem. In this study, a detailed model of a semibatch distillation process was validated directly using the industrial equipment on site, and then this model was used as the equality constraints to optimize the process. To carry out the experiment, a conventional run of the process was chosen and measured using the existing apparatus. Modeling the tray-by-tray column with a reversible reaction in the reboiler leads to a set of nonlinear DAEs of a high dimension. For simulating the process, collocation on finite elements was used to discretize the system and the DAEs can be integrated with the Newton method from interval to interval. Good agreement between the measured and the simulated profiles of the process have been achieved, which implies that the model can accurately describe the process and be used as part of the optimization. Additionally, through analysis of the experimental and simulation results, it is possible to obtain useful information as to the potential of the optimization being performed. Minimization of the batch time was considered as the aim of the optimization. Together with the equality and inequality constraints a complex dynamic optimization problem has been formulated. Both the reflux ratio and the feed flow rate policies were optimized with a modified approach containing the SQP algorithm for finding the policies and collocation method for simulating the process as well as computing the necessary sensitivities. From the optimization results of two operating scenarios, we demonstrate the effectiveness and usability of the approach. Generally speaking, during the main-cut period the reflux ratio should be reduced as far as possible to quickly remove the product
alcohol so as to increase the reaction rate, while in the off-cut period it should be increased to keep more of the educt alcohol in the reboiler. The feed of educt alcohol should be primarily introduced into the reboiler in the main-cut period because it will be distillated out of the reboiler due to the rising temperature during the offcut period. But at the end of the batch, a peak of the feed rate is necessary to eliminate the educt ester left in the reboiler, thus satisfying the strict purity requirement on the educt ester. The optimal policies developed can be used to guide and improve the operation of the process. However, the modeling and optimization approaches are not limited to the process considered; they can be applied to solve general dynamic optimization problems of large-scale systems described in terms of DAEs. Despite the fact that the solution found by this approach is feasible, it is quite easy to locate a suboptimum or near optimum due to the nonconvexity and complexity of the problem. In some cases, we had to restart the computation by initializing the approach with a computed local optimum as new estimated profiles for the independent variables. In addition. because of the demand of process simulation and the inversion of the Jacobian at each integration step for calculating the sensitivities. the computational expense is large: the CPU time of a SUN SPARC 20 for an optimization run of the considered problem amounts to about 3 h. Therefore, development of a global and efficient searching method remains an important task for batch distillation optimization. In addition, the implementation of the optimal policies developed for the industrial process and development of a proper approach for the on-line reoptimization, in situations where the process is disturbed and diverges from the desired trajectories, are areas for future work. Acknowledgment Financial support from the scholarship program to Chinese researchers by Technische Universita¨t Berlin is gratefully acknowledged. Nomenclature a ) heat capacity parameters A, B, C ) Antoine constants C ) molar concentration, mol/m3 cp ) heat capacity, kJ/kmol K D ) distillate flow rate, kmol/s E ) activity energy, J/mol f ) normal function F ) feed flow rate, kmol/s g ) array of equality constraints h ) enthalpy of a component, kJ/kmol H ) enthalpy of a mixture, kJ/kmol HU ) holdup, kmol J ) objective function k ) reaction rate constant k0 ) frequency factor K ) phase equilibrium constant L ) liquid flow rate, kmol/s M ) molecular weight, g/mol NC ) number of collocation points NK ) number of components NL ) number of time intervals NST ) number of trays p ) pressure, bar p0 ) component partial pressure, bar
1348 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 p ) array of inequality constraints r ) reaction rate, kmol/s m3 R ) normal gas constant, J/mol K RV ) reflux ratio t ) time, s T ) temperature, K u ) array of control (independent) variables V ) vapor flow rate, kmol/s VR ) reaction volume, m3 x ) liquid composition, mol/mol x ) array of dependent variables y ) vapor composition, mol/mol
NK
Ki,1xi,1 ) 1 ∑ i)1 L1 )
RV )
L1 D
(A4)
Internal Trays: j ) 2, NST - 1.
component balance:
Greek Letters
HUj
η ) tray efficiency γ ) activity coefficient F ) density of a component, g/cm3
dxi,j ) Lj-1xi,j-1 + Vj+1yi,j+1 - Ljxi,j - Vjyi,j (A5) dt
vapor-liquid equilibrium:
Subscripts 0 ) start f ) end H ) forward reaction i ) index of components or collocation points j ) index of trays or collocation points l ) index of time intervals R ) back reaction
yi,j ) ηjKi,jxi,j + (1 - ηj)yi,j+1
(A6)
summation equations: NK
∑ i)1
Superscripts k ) index of iteration L ) liquid V ) vapor
NK
xi,j ) 1;
yi,j ) 1 ∑ i)1
(A7)
energy balance:
Appendix A: Model of the Semibatch Distillation Process with a Chemical Reaction in the Reboiler The model considered in the simulation and the optimization is based on the following assumptions: (1) constant molar holdup in the condenser and on the trays; (2) constant tray pressure and tray efficiency; (3) total condenser without subcooling; (4) ideal vapor phase on each tray; (5) negligible vapor holdup; (6) ideal heat exchange in the condenser and reboiler. The last assumption means that the heat needed for the vapor condensation and boilup can be obtained from the heat exchange. The trays are numbered from the condenser (j ) 1) to the reboiler (j ) NST). The component compositions of both liquid and vapor phases, the liquid and vapor flow rates, and the temperature are the variables of each tray. With the above assumptions, we have the following model in terms of DAEs. Total Condenser: j ) 1.
component balance: dxi,1 HU1 ) V2yi,2 - (L1 + D)xi,1 dt
(A1)
summation equation: xi,1 ) 1 ∑ i)1
dHLj L V + Vj+1Hj+1 - LjHLj - VjHVj HUj ) Lj-1Hj-1 dt (A8) Reboiler: j ) NST.
component balance: dHUNSTxi,NST ) dt LNST-1xi,NST-1 - VNSTyi,NST + FNSTxFi + ri,NSTVR,NST (A9) vapor-liquid equilibrium: yi,NST ) Ki,NSTxi,NST
(A2)
Instead of using the energy balance equation, for a total condenser at the dew point it can be replaced with the following relation:
(A10)
summation equations: NK
xi,NST ) 1; ∑ i)1
NK
yi,NST ) 1 ∑ i)1
(A11)
Similar to the total condenser, the energy balance of the reboiler can be replaced as follows: with the assumption of an ideal heat exchange the energy needed for the evaporation will be obtained from the exchanger. This means
VNST - LNST-1 ) D
NK
total mass balance:
RVV2 ; 1 + RV
(A3)
(A12)
total mass balance: dHUNST ) LNST-1 - VNST + FNST dt volume of the reaction:
(A13)
(∑ ) i,NSTMi
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1349
NK x
VR,NST ) HUNST
i)1Fi(TNST)
(A14)
min J[u(t),tf] ) G1[x(tf),tf] +
Ki,j )
pj
γi,j(Tj,xi,j)
(A15)
where the partial pressure of a component will be computed with the Antoine equation:
ln p0i,j ) Ai +
Bi Ci + Tj
p(x3 ,x,u,t) g 0
and the activity coefficient γi,j can be calculated by either the NRTL, UNIQUAC, or UNIFAC model. Enthalpy. NK
HVj )
∑ i)1
NK
yi,jhVi,j
HLj )
xi,jhLi,j ∑ i)1
(A17)
umin e u e umax x(0) ) x0 For discretization, with the Lagrange polynomials with shifted roots of Legendre polynomials to represent each dependent variable in each time interval
∫
i
- tj
xi,j
(A23)
To keep the continuity of the variables between intervals, we use the last collocation point as the starting point for the next interval. Thus, the problem (A22) can now be described with an NLP problem:
min J(xl,iul,i) s.t.
gl,i(xl,i,ul,i) ) 0 (A24)
xmin e xl,i e xmax umin < ul,i < umax
T
c dT 298 pi,j
cpi,j ) ai,0 + ai,1Tj + ai,2Tj2 + ai,3Tj3
(A18)
(A19)
Reaction Rate. The rate of the reversible reaction in the reboiler can be described as
ri,NST ) kHCA,NSTCB,NST - kRCC,NSTCD,NST
(A20)
)
E RTNST
where the variables are the values of the state and control variables of (A22) at the collocation points. Since there are many variables in the optimization problem of a batch distillation, the SQP algorithm will become inefficient. To reduce the number of variables, we include only the control (independent) variables in the SQP algorithm, while the dependent variables are solved through integration of the model equations, through which the sensitivities of the objective function and the constraints are computed. Now the optimization problem can be simplified as
min J[f(ul,i),ul,i]
where the reaction constants can be calculated using the Arrhenius equation:
(
t - tj
∑ ∏ i)1 j)0 t
pl,i(xl,i,ul,i) g 0
and
k ) k0 exp -
[ ]
NC NC
xl(t) )
where the vapor and liquid molar enthalpy of a component can be computed respectively with
hi,j ) ∆H298i +
(A22)
xmin e x e xmax
j*i
(A16)
f
g(x3 ,x,u,t) ) 0
s.t.
In the above model equations, there are several important auxiliary variables which depend on the dependent variables and should be calculated from the following known physical and chemical relations. Phase Equilibrium.
p0i,j
∫0t G2(x,u,t) dt
s.t.
pl,i[f(ul,i),ul,i] g 0
(A25)
umin e ul,i e umax (A21)
There are (2NK + 3) variables and the same number of equations per tray. Thus the above DAE system can be discretized (e.g., an implicit Euler or a collocation method) and solved by using a global Newton method at each integration step. Appendix B: The Modified Dynamic Optimization Approach Based on SQP and Collocation on Finite Elements A general dynamic optimization problem can be described with
Now only the control variables and the inequality constraints appear in the problem to be solved. For the sensitivity computation, the gradients of the dependent variables to controls need to be computed. Since the relation x ) f(u) is not known, the gradients are implicitly computed with the following method. In the kth iteration of SQP and in the lth interval the model equations of the collocation points can be described as k gl(xl,0 ,xkl ,ukl ) ) 0
(A26)
where x is dependent on u and x0 (initial value of the variables at the interval). Through the first-order Taylor expansion of (A26), we obtain
1350 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 k D0kl (xl,0 - xl,0 ) + DXkl (xl - xkl ) + DUkl (ul - ukl ) ) 0 (A27)
where D0kl , DXkl , and DUkl are the gradient matrices of the related variables. From (A27) we can determine the sensitivities of the dependent variables to the controls and the start values:
∂xl ) -(DXkl )-1DUkl ∂ul
(A28)
∂xl ) -(DXkl )-1D0kl ∂xl,0
(A29)
Since the last collocation point is used as the starting point of the next interval, the sensitivities of the initial values of the next interval are part of those of the dependent variables. Thus, the gradients of the variables over the batch time can be computed and moved from interval to interval:
∂xl+1,0 ∂xl+1,0 ∂xl,0 ) ∂ul-1 ∂xl,0 ∂ul-1
(A30)
Literature Cited Ahmad, S. B.; Barton, P. I. Homogeneous multicomponent azeotropic batch distillation. AIChE J. 1996, 42, 3419. Albet, J.; Le Lann, J. M.; Joulia, X.; Koehret, B. Rigorous simulation of multicomponent multisequence batch reactive distillation. In Proceedings of COPE’91, Barcelona, Spain, 1991; pp 75-80. Barolo, M.; Guarise, G. B.; Rienzi, S. A.; Trotta, A.; Macchietto, S. Running batch distillation in a column with a middle vessel. Ind. Eng. Chem. Res. 1996, 35, 4612. Cuthrell, J. E.; Biegler, L. T. On the optimization of differentialalgebraic process systems. AIChE J. 1987, 33, 1257. Diwekar, U. M.; Madhavan, K. P. Batch-Dist: a comprehensive package for simulation, design, optimization and optimal control of multicomponent, multifraction batch distillation columns. Comput. Chem. Eng, 1991, 15, 833. Dussel, R.; Stichlmair, J. Zerlegung azeotroper Gemische durch Batch-Rektifikation unter Verwendung eines Zusatzsstoffes. Chem.-Ing-Tech. 1994, 66, 1061. Fahart, S.; Czernicki, M.; Pibouleau, L.; Domenech, S. Optimization of multiple-fraction batch distillation by nonlinear programming. AIChE J. 1990, 36, 1349. Finlayson, B. A. Nonlinear Analysis in Chemical Engineering; McGraw-Hill: New York, 1980. Gmehling, J.; Onken, U.; Arlt, W. Vapor-Liquid Data Collection; DECHEMA: Frankfurt, 1977. IMSL, MATH/LIBRARY User’s Manual; IMSL Inc.: Houston, TX, 1987. Kister, H. Z. Distillation Design; McGraw-Hill: New York, 1992.
Li, P.; Hoo, H. P.; Wozny, G. Efficient simulation of batch distillation processes by using orthogonal collocation. AIChE Annual Meeting, Los Angeles, CA, Nov 1997a; Paper 204g. Li, P.; Wozny, G.; Reuter, E. Optimization of multiple-fraction batch distillation with detailed dynamic process model. Distillation and Absorption ‘97, Maastricht, The Netherlands, September 1997; Darton, R., Ed.; IChemE: Rugby, 1997b; IChernE Symposium Series No. 142; pp 289-300. Logsdon, J. S.; Biegler, L. T. Accurate determination of optimal reflux policies for the maximum distillate problem in batch distillation. Ind. Eng. Chem. Res. 1989, 28, 1628. Logsdon, J. S.; Biegler, L. T. Decomposition strategies for largescale dynamic optimisation problems. Chem. Eng. Sci. 1992, 47, 851. Mujtaba, I. M.; Macchietto, S. An optimal recycle policy for multicomponent batch distillation. Comput. Chem. Eng. 1992a, 16, 273. Mujtaba, I. M.; Macchietto, S. Optimal operation of reactive batch distillation. AIChE Annual Meeting, Miami Beach, FL, November 1992b; Paper 135g. Nad, M.; Spiegel, L. Simulation of batch distillation by computer and comparison with experiment. Comput. Chem. Eng/EFCE Giardini Naxos, Italy, 1987. Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids; McGraw-Hill: New York, 1987. Renfro, J. G.; Morshedi, A. M.; Asbjornsen, O. A. Simultaneous optimization and solution of systems described by differential/ algebraic equations. Comput. Chem. Eng. 1987, 11, 503. Reuter, E. Simulation und Optimierung einer chemischen Reaktion mit uberlagerter Rektifikation. VDI Fortschritt-Berichte; Reihe 3: Verfahrenstechnik, Nr. 394, 1994. Reuter, E.; Wozny, G.; Jeromin, L. Modeling of multicomponent batch distillation processes with chemical reaction and their control systems. Comput. Chem. Eng. 1989, 13, S499. Skogestad, S.; Wittgens, B.; Litto, R. Multivessel batch distillation. AIChE J. 1997, 43, 971. Sørensen, E.; Macchietto, S.; Stuart, G.; Skogestad, S. Optimal control and on-line operation of reactive batch distillation. Comput. Chem. Eng, 1996, 20, 1491. Vassiliadis, V. S.; Pantelides, C. C.; Sargent, R. W. H. Solution of a class of multistage dynamic optimization problems. 1. Problems without path constraints. Ind. Eng. Chem. Res. 1994a, 33, 2111. Vassiliadis, V. S.; Pantelides, C. C.; Sargent, R. W. H. Solution of a class of multistage dynamic optimization problems. 2. Problems with path constraints. Ind. Eng. Chem. Res. 1994b, 33, 2123. Wajge, R. M.; Wilson, J. M.; Pekny, J. F.; Reklaitis, G. V. Investigation of numerical solution approaches to multicomponent batch distillation in packed beds. Ind. Eng. Chem. Res. 1997, 36, 2287.
Received for review September 30, 1997 Revised manuscript received January 2, 1998 Accepted January 5, 1998 IE970695L
