Nonlinear Model Predictive Control of Reactive Distillation Based on

Ind. Eng. Chem. Res. 2008, 47, 6949–6960

6949

Nonlinear Model Predictive Control of Reactive Distillation Based on Stochastic Optimization Ch. Venkateswarlu* and A. Damodar Reddy Process Dynamics and Control Group, Chemical Engineering Sciences DiVision, Indian Institute of Chemical Technology, Hyderabad-500 007, India

Stochastic optimization algorithms such as genetic algorithm (GA) and simulated annealing (SA) are combined with a polynomial-type empirical process model to develop nonlinear model predictive control (NMPC) strategies, namely, GANMPC and SANMPC, in the perspective of control of a nonlinear reactive distillation column. In these strategies, the nonlinear input-output process model is cascaded itself to generate future predictions for the process output based on which the control sequence is computed by stochastic optimizers while satisfying the specified performance criteria. The performance of the proposed controllers is evaluated by applying to single input-single output (SISO) control of an ethyl acetate reactive distillation column with double-feed configuration involving an esterification reaction with azeotropism. The results demonstrate better performance of the stochastic optimization based NMPCs over a conventional proportional-integral (PI) controller, a linear model predictive controller (LMPC), and a NMPC based on sequential quadratic programming (SQP) in tracking the setpoint changes as well as stabilizing the operation in the presence of input disturbances. Although both the GANMPC and SANMPC are found to exhibit almost equal performance, the easier tuning and the lower computational effort suggests the better suitability of SANMPC for the control of a nonlinear reactive distillation column. 1. Introduction Model predictive control (MPC) is an important branch of automatic control theory that has been widely applied in industry.1 MPC relies on a process model that predicts the behavior of the process over some future time interval, and an optimal sequence of control actions is determined based on these model predictions by minimizing an objective function. In order to realize optimal performance, the process model used in MPC needs to be fairly accurate and the optimization algorithm should be efficient. Linear MPC employs linear or linearized models to obtain the predictive response of the controlled process. Linear MPC has proved useful for controlling processes that exhibit even some degree of nonlinear behavior.2,3 However, as the mismatch between the actual process and the representative model increases, the degree of deterioration in the control performance increases. Thus, the control of a highly nonlinear process by MPC requires a suitable model that represents the salient nonlinearities of the process. Basically, two different approaches are used to develop nonlinear dynamic models. These approaches are developing a first-principle model using available process knowledge and developing an empirical model from input-output data. The first-principle modeling approach results in models in the form of coupled nonlinear ordinary differential equations, and various model predictive controllers based on this approach have been reported for nonlinear chemical processes.4,5 The first-principle models will be larger in size for high-dimensional systems, thus limiting their usage for MPC design. On the other hand, the input-output modeling approach can be conveniently used to identify nonlinear empirical models from plant data, and there has been a growing interest in the development of different types of MPCs based on this approach.6,7 The other nonlinear model predictive control approaches include an approach8 involving a nonlinear model for free response prediction and a linear model for force * Corresponding author. Tel: +91-40-27193121. Fax: +91-4027193626. E-mail: [email protected].

prediction, an extended Kalman filter based strategy,9 and a quasi-infinite horizon nonlinear model predictive control scheme.10 One more important aspect in model predictive control of highly nonlinear systems is the optimization algorithm. Efficient optimization algorithms exist for convex optimization problems. However, the optimization problem often becomes nonconvex in the presence of nonlinear characteristics/ constraints and is usually more complex than convex optimization. According to Camacho and Bordons,8 the practical usefulness of nonlinear predictive control is hampered by the unavailability of suitable optimization tools. Sequential quadratic programming (SQP) is a widely used classical optimization algorithm to solve nonlinear optimization problems. However, for the solution of large problems, it has been reported that gradient-based methods like SQP require more computational efforts.11 Moreover, classical optimization methods are more sensitive to the initialization of the algorithm and usually lead to unacceptable solutions due to convergence to local optima. These methods are not expected to provide the global optimum for a multimodal function. Consequently, new optimization techniques are being proposed to achieve efficient control performance. Stochastic search and optimization algorithms such as genetic algorithm (GA) and simulated annealing (SA) derived from the principles of natural phenomena are useful to find the global optimum of complex engineering problems. These algorithms are attractive because of their flexibility, ease of operation, and global perspective. In a recent work on dynamic optimization, Li et al.12 have shown that an SA-based stochastic optimizer requires much less computational effort than SQP. Reactive distillation, the process of performing chemical reaction and multistage distillation simultaneously, has emerged as a favorable alternative to conventional reactor-separator configurations. This combined operation is especially useful for the chemical reactions in which chemical equilibrium limits the conversion and can be advantageously used for those reactions that occur at temperatures and pressures suitable to the distillation of components. When reactions are equilibrium limited,

10.1021/ie070972g CCC: $40.75  2008 American Chemical Society Published on Web 08/13/2008

6950 Ind. Eng. Chem. Res., Vol. 47, No. 18, 2008

reactive distillation continuously removes products from the reaction zone, which significantly increases the overall conversion. Reactive distillation may also increase selectivity in certain competing reaction systems by continuously separating products from reactants. Another advantage is that separation by reactive distillation can overcome azeotropic limitations on achievable separation. These advantages of reactive distillation over conventional configurations of reactors followed by separators have motivated a renewed interest in the use of reactive distillation technology for the production of important chemicals.13 However, in reactive distillation, the interaction between the simultaneous reaction and distillation introduces a much more complex behavior compared to conventional processes and leads to challenging problems in design, optimization, and control. Despite increasing attention toward dynamic modeling and simulation, the literature available on control of reactive distillation is relatively small. Recently, various studies examining the column design and control structures for reactive distillation have been reported. Sneesby et al.14 have studied several control configurations with proportional-integral (PI) controllers for ethyl tert-butyl ether (ETBE) reactive distillation. Among these configurations, LV and LB configurations with single-composition control have shown to provide better performance. Al-Arfaj and Luyben15 have considered different designs for a single-feed pentene metathesis reactive distillation with various single input-single output (SISO) control structures involving PI controllers. Simulation results have shown that a control structure that uses two temperatures to maintain the desired purities of top and bottom products was found to be more effective. Further studies on control of reactive distillation have shown the inadequacy of conventional controllers and highlight the need for more advanced controllers. Toward the application of advanced controllers, a nonlinear controller based on a rigorous mathematical model has been derived and used for an ethylene glycol reactive distillation column with kinetically limited reversible reactions.16 Simulation results have shown the better performance of a nonlinear controller over linear PI controllers. However, Al-Arfaj and Luyben15 have shown that an ethylene glycol reactive distillation column can be effectively controlled by a simple PI control scheme. Vora and Daoutidis17 have studied the application of model-based linear and nonlinear state feedback controllers and SISO controllers with PI configuration for ethyl acetate reactive distillation. The simulation results of these controllers have shown the better performance of the nonlinear model-based controller. This work is focused toward the development of nonlinear model predictive control strategies for reactive distillation based on stochastic optimization techniques. Stochastic optimization algorithms such as genetic algorithms (GAs) and simulated annealing (SA) are combined with the polynomial-type nonlinear empirical process model identified from input-output data of the process to develop nonlinear model predictive control (NMPC) strategies, namely, GANMPC and SANMPC. In these strategies, the nonlinear input-output process model is cascaded itself to generate future predictions for the process output, based on which the control sequence is computed by stochastic optimizers while satisfying the specified performance criteria. The performance of the proposed controllers is evaluated by applying to single input-single output (SISO) control of an ethyl acetate reactive distillation column with double-feed configuration involving an esterification reaction and compared with other conventional, linear, and nonlinear model predictive controllers.

2. Stochastic Optimization Algorithms Stochastic approach based genetic algorithms (GAs) and simulated annealing (SA) are potential tools for optimization because of their ability to handle constrained, nonlinear, and nonconvex optimization problems. These methods have the capacity to escape local optima and find solutions in the vicinity of the global optimum. They have the ability to use the values from the model in a black box optimization approach without requiring the derivatives. Various studies have been reported to demonstrate the ability of these methods in providing efficient optimization solutions.18,19 2.1. Genetic Algorithms. Genetic algorithms (GAs) are stochastic search algorithms that are based on the mechanics of natural selection and genetics.20 A GA aims at evolving an optimal solution by a randomized information exchange through a sequence of probabilistic transformations governed by a selection scheme biased toward high-quality solutions. Since GAs use a population of potential solutions and the probabilistic transition rules to create new solutions, they usually have a higher probability of finding the global optimum.20 Many applications of GAs have been reported to solve different optimization problems.21,22 GAs encode the candidate solutions of an optimization algorithm as a string of characters, which are usually binary digits. The string is called a chromosome, and the variables that are coded into a chromosome form substrings. The length of the string is usually determined according to the solution accuracy. The genetic algorithm considers a number of random strings for the variables, which together form a population. GA modifies and updates the population iteratively, searching for good solutions of the optimization problem. Each iteration step is called a generation. GA begins with a population of random strings representing the variables and then evaluates the strings fitness by using a fitness function, F(x), defined by F(x) )

1 1 + f(x)

(1)

where f(x) is an objective function. The specification of the fitness function is an important aspect of GA design, because the solution of the optimization problem and the performance of the algorithm depend on this function. Inspired by the “survival of the fittest” idea, the GA maximizes/minimizes the fitness value. The fitness function evaluates the fitness of every individual in the population until the end condition, such as the generation number or degree of convergence is satisfied. If the end condition is not satisfied, the next step will be population evolution using genetic operators. 2.1.1. Population Evolution. The algorithm starts with the generation of an initial population. This population contains individuals that represent initial estimates for the optimization problem. The fitness of every individual in the population is obtained through a cost function evaluation and new individuals are generated for the next generation. The population evolution for the next generation is performed using genetic operators such as selection, crossover, and mutation. Selection chooses individuals from the previous population for reproduction according to their fitness values. Individuals with better fitness values survive to form a mating pool. Different selection operators are used in selection. The essential idea of these operators is to pick up the above-average strings from the current population and insert their multiple copies in the mating pool in a probabilistic manner. Crossover is applied after selection, and this produces new individuals. Crossover

Ind. Eng. Chem. Res., Vol. 47, No. 18, 2008 6951

consists of merging the chromosomes of two individuals (parents) to obtain two other individuals (children). Parent pairs are randomly chosen from the selected population, and a crossover operator is used for merging. Mutation is the last operation, in which a mutation operator is applied to perform bitwise mutation with a specified mutation probability. In mutation, better strings are created by carrying out a local search around the current solution and diversity in the population is maintained. Selection, crossover, and mutation are repeated for a fixed number of generations. 2.2. Simulated Annealing. Simulated annealing (SA) is analogous to the process of atomic rearrangement of a substance into a highly ordered crystalline structure by way of slowly cooling-annealing the substance through successive stages. This method is found to be a potential tool to solve a variety of optimization problems.23–25 Crystalline structure with a high degree of atomic order is the purest form of the substance, indicating the minimum energy state. The principle of SA mimics the annealing process of slow cooling of molten metal to achieve the minimum function value. The cooling phenomena is simulated by controlling a temperature-like parameter introduced with the concept of the Bolzmann probability distribution, which determines the energy distribution probability, P, of the system at the temperature, T, according to the equation P(E) ) exp(-E ⁄ kBTA)

(2)

where kB is the Bolzmann constant. The Bolzmann distribution concentrates around the state with the lowest energy. For T f 0, P(E) f 0 and only the state with the lowest energy can have a probability greater than zero. However, cooling the system too fast could result in a higher state of energy and may lead to freezing the system to a metastable state. SA is a point-by-point method based on Monte Carlo approach. The algorithm begins at an initial random point called u, and a high temperature T, and the function value at this point is evaluated as E(k). A second point is created in the vicinity of the initial point u, and the function value corresponding to this point is obtained as E(k + 1). The difference in function values at these points ∆E is obtained as ∆E ) E(k + 1) - E(k) (3) If ∆E e 0, the second point is accepted; otherwise, the point is accepted probabilistically, governed by the temperature-dependent Bolzmann probability factor: Pr ) exp(-∆E ⁄ kBTA)

(4)

The annealing temperature, TA, is a parameter in the optimization algorithm and is set by a predefined annealing schedule, which starts at a relatively high temperature and steps slowly downward at a prescribed rate in accordance with the equation TAm+1 ) RTAm, 0 < R < 1

(5)

As the temperature decreases, the probability of the acceptance of the point u will be decreased according to eq 4. The parameter R is set such that, at the point of convergence, the temperature TA reaches a small value. The procedure is iteratively repeated at each temperature with the generation of new points, and the search is terminated when the convergence criterion set for the objective is met. 3. Nonlinear Modeling and Model Identification Various model structures, such as Volterra series models, Hammerstein and Wiener models, bilinear models, state affine

Figure 1. Structure of stochastic optimization based NMPC.

models, and neural network models have been reported in the literature for identification of nonlinear systems. Haber and Unbehauen26 presented a comprehensive review on these model structures. The model considered in this study for identification of a nonlinear reactive distillation has a polynomial ARMA structure of the form ny

ˆy(k) ) θ0 +

∑θ i)1

ny+1

∑θ

nu

1,iy(k - i) +

∑θ

2,iu(k - i) + i)1 nu i+1

∑∑θ

4,iju(k - i)u(k - j) + ...

(6)

ˆy(k) ) f[y(k - 1), ..., y(k - ny), u(k - 1), ..., u(k - nu)]

(7)

3,iy(k - i)u(k - i) +

i)1

i)1 j)1

or simply

Here, k refers to the sampling time, y and u are the output and input variables, and ny and nu refer to the number of output and input lags, respectively. These types of polynomial model structures have been used by various researchers for process control.27,28 The main advantage of this model is that it represents the process nonlinearities in a structure with linear model parameters, which can be estimated by using efficient parameter estimation methods such as recursive least-squares. Thus, the model in eq 7 can be rearranged in a linear regression form as ˆy(k) ) θT(k - 1)φ(k - 1) + ε(k)

(8)

where θ is a parameter vector, φ represents input-output process information, and ε is the estimation error. The parameters in the model can be estimated by using recursive leastsquares based on a priori process knowledge representing the process characteristics over a wide range of operating conditions. 4. Design of Nonlinear Model Predictive Controllers In NMPC design, the identified input-output nonlinear process model is explicitly used to predict the process output at future time instants over a specified prediction horizon. A sequence of future control actions over a specified control horizon is carried out using stochastic optimizers that minimize


the objective function under given operating constraints. In this receding horizon approach, only the first control action in the sequence is implemented. The horizons are moved toward the future. The structure of the stochastic optimizer based NMPC is shown in Figure 1. 4.1. Predictive Model Formulation. The primary purpose of MPC is to deal with complex dynamics over an extended horizon. Thus, an MPC model must predict the process dynamics over a prediction horizon, enabling the controller to incorporate future setpoint changes or disturbances. The polynomial input-output model provides one-step-ahead prediction for process output. By feeding back the model outputs and control inputs, the one-step-ahead predictive model can be recurrently cascaded to itself to generate future predictions for process output. The N step predictions can be obtained as follows: ˆy(k + 1 ⁄ k) ) f[y(k), ..., y(k + 1 - ny), u(k), ..., u(k + 1 - nu)] ˆy(k + 2 ⁄ k) ) f[yˆ(k + 1 ⁄ k), ..., y(k + 2 - ny), u(k + 1), ..., u(k + M - nu)] l l ˆy(k + N ⁄ k) ) f[yˆ(k + N - 1 ⁄ k), ..., ˆy(k + N - ny ⁄ k), u(k + M - 1), ..., u(k + M - nu)] (9) where N is the prediction horizon and M is the control horizon. 4.2. Objective Function. The optimal control input sequence in NMPC is computed by minimizing an objective function based on a desired output trajectory over a prediction horizon, N

Min J )

∑ λ[w(k + i) - ˆy (k + i)]

2

p

+

i)1

M

∑ γ[∆u(k + i - 1) ∆u(k + i - 1)] T

i)1

u(k), u(k + 1), ..., u(k + M - 1)

(10)

subject to constraints ymin e ˆyp(k + i) e ymax (i ) 1, ..., N) umin e u(k + i) e umax (i ) 0, ..., M - 1) ∆umin e ∆u(k + i) e ∆umax (i ) 0, ..., M - 1) where yˆp(k + i), i ) 1,..., N, are the future process outputs predicted over the prediction horizon; wk+i, i ) 1,..., N, are the setpoints; and u(k + i), i ) 0,..., M - 1, are the future control signals. The λ and γ represent the output and input weightings, respectively. The umin and umax are the minimum and maximum values of the manipulated inputs, respectively, and ∆umin and ∆umax represent their corresponding changes. Computation of future control signals involves the minimization of the objective function so as to bring and keep the process output as close as possible to the given reference trajectory, even in the presence of load disturbances. The control actions are computed online by solving an optimization problem while taking into consideration of constraints on the outputs and inputs. The control signal, u, is manipulated only within the control horizon and remains constant afterward, i.e., u(k + i) ) u(k + M - 1) for i ) M,..., N - 1. Only the first control move of the optimized control sequence is implemented on the process, and the output measurements are obtained. At the next sampling instant, the prediction and control horizons are moved ahead by one step, and the optimization problem is solved again using the updated measurements from the process. The mismatch dk between the process y(k) and the model yˆ(k) is computed as

dk ) b(y(k) - ˆy(k))

(11)

where b is a tunable parameter lying between 0 and 1. This mismatch is used to compensate the model predictions in eq 8: ˆyp(k + i) ) ˆy(k + i) + dk (for i ) 1 to N)

(12)

These predictions are incorporated in the objective function defined by eq 10 along with the corresponding setpoint values. 4.3. GA-Based Nonlinear Model Predictive Controller (GANMPC). In GANMPC design, the principles for fitness evaluation and the representation of the variables in genes are specified. The control input, u, is normalized and constrained within the specified limits. The GA coding for ∆u has to ensure that its decoded value should lie within the specified constraints. The number of input variables coded in a string equals the control horizon, and the length of the string increases as the control horizon increases. If each of the ∆u over the control horizon, M, is represented by a substring of s bits, then the string length over the control horizon becomes sM. A penalty function approach is considered to satisfy the constraints on the input variables. In this approach, a penalty term corresponding to the penalty violation is added to the objective function defined in eq 10. Thus, the violation of the constraints on the variables is accounted for by defining a penalty function of the form N

P)

∑ µ(∆u(k + i))

2

(13)

i)1

where the penalty parameter, µ, is selected as a high value. The penalized objective function is then given by f(x) ) J + P (14) The fitness function is formed by the usual transformation as in eq 1. The operation of GA begins with a population of random strings representing the manipulated input variables, ∆uk+i for i ) 0,..., M - 1. The population is then operated through selection, crossover, and mutation to create a new population of points representing the future control actions over a specified control horizon. The updated values of uk+i are used to compute the model predictions explicitly at future time instants over a specified prediction horizon in order to evaluate the objective function. The procedure is continued iteratively until the objective function is minimized while satisfying the specified termination criteria. The first value of the uk+i is implemented on the process, and the procedure repeated for the next process measurement condition. 4.3.1. Implementation Procedure. The implementation of GANMPC proceeds with the following steps. (1) Initialize the algorithm with the parameters used for processing the information. Choose the initial control vector for u and ∆u. Specify the generation number. (2) Find the process outputs and compute the model predictions using eq 9. (3) Perform GA search to find the ∆u that optimizes the cost function and satisfies the constraints. This is accomplished by performing selection, crossover, and mutation on the random population while evaluating the fitness of the solution. (4) Update the control vector as u(k + i) ) u(k) + ∆u(k) (15) (5) Terminate GA if the specified number of generations is reached. (6) Go to step 2 and repeat the procedure for every sampling point based on the updated control vector and its corresponding process output using the random population for ∆u.


4.4. SA-Based Nonlinear Model Predictive Controller (SANMPC). SANMPC design needs to specify the energy function and the random number selection for control input calculation. The normalized control input and its constraints used in SANMPC are the same as in GANMPC. The random numbers used for the control input, ∆u, equal the length of the control horizon, and these numbers are generated so as to satisfy the constraints. Constraint violation is accounted for by using the penalty function procedure as in the case of GA controller. At any instant, the current control signal, uk, and the prediction output based on this control input, yˆp(k + i), are used to compute the objective function J as the energy function, E(k + i). The E(k + i) and the previously evaluated E(k) provide the ∆E as ∆E(k) ) E(k + i) - E(k)

(16)

The comparison of the ∆E with the random numbers generated between 0 and 1 determines the probability of acceptance of u(k). If ∆E e 0, all u(k) are accepted. If ∆E > 0, u(k) are accepted with a probability of exp(-∆E/TA). If nm is the number of variables, nk is the number of function evaluations, and nT is the number of temperature reductions, then the total number of function evaluations required for every sampling condition are (nT × nk × nm). 4.4.1. Implementation Procedure. The implementation of SANMPC proceeds with the following steps. (1) Set TA as a sufficiently high value and let nk be the number of function evaluations to be performed at a particular TA. Specify the termination criterion, ε. Choose the initial control vector, u, and obtain the process output predictions using eq 9. Evaluate the objective function, eq 14, as the energy function E(k). (2) Compute the incremental input vector ∆uk stochastically and update the control vector as in eq 15. Calculate the objective function, E(k + i), as the energy function based on this vector. (3) Accept u(k + i) unconditionally if the energy function satisfies the condition E(k + i) e E(k)

(17)

Otherwise, accept u(k + i) with the probability according to the Metropolis criterion

(

)

(E(k + i) - E(k)) gr exp TA′

(18)

where TA′ is the current annealing temperature and r represents a random number. This step proceeds until the specified function evaluations, nk, are completed. (4) Carry out the temperature reduction in the outer loop according to the decrement function TA′ ) RTA

(19)

where R is the temperature reduction factor. Terminate the algorithm if all the differences are less than the prespecified ε. (5) Go to step 2 and repeat the procedure for every measurement condition based on the updated control vector and its corresponding process output. 5. Process Representation Ethyl acetate is produced through an esterification reaction between acetic acid and ethyl alcohol: H+

CH3COOH + C2H5OH 798 H2O + CH3COOC2H5 (20)

Figure 2. Control structure of two-feed ethyl acetate reactive distillation column.

The achievable conversion in this reversible reaction is limited by the equilibrium conversion. This quaternary system is highly nonideal and forms binary and ternary azeotropes, which introduce complexity to the separation by conventional distillation. Reactive distillation can provide a means of breaking the azeotropes by altering or eliminating the conditions for azeotrope formation. Thus, reactive distillation becomes an attractive alternative for the production of ethyl acetate. The rate equation of this reversible reaction in the presence of a homogeneous acid catalyst is given by Alejski and Duprat29 as k1 CC Kc 3 4 k1 ) (4.195Ck + 0.08815) exp(-6500.1⁄T) Kc ) 7.558 - 0.012T r1 ) k1C1C2 -

(21)

All the plates in the reactive distillation column are considered to be reactive in the sense that reaction takes place on all plates including the condenser and the reboiler. Vora and Daoutidis17 have presented a two-feed column configuration for ethyl acetate reactive distillation and found that feeding the two reactants on different trays countercurrently allows one to enhance the forward reaction on trays and results in higher conversion and purity over the conventional column configuration, which involves feeding the reactants on a single tray. In this work, the design and performance evaluation of the proposed stochastic optimization based NMPCs is carried out by using the doublefeed column configuration of ethyl acetate reactive distillation with the control structure shown in Figure 2. The dynamic model representing the process involves mass and component balance equations with reaction terms, as well as algebraic energy equations supported by vapor-liquid equilibrium and physical properties.29 The assumptions made in the formulation of the model include adiabatic column operation, negligible heat of reaction, negligible vapor holdup, liquid-phase reaction, and physical equilibrium in streams


leaving each stage. The assumption of negligible energy accumulation on the plates changes differential energy balance equations to the algebraic equations. The equations representing the process are given in Appendix A. 6. Results and Discussion The ethyl acetate reactive distillation column consists of 13 stages including reboiler and condenser, in which the less volatile acetic acid enters the 3rd tray and the more volatile ethanol enters the 10th tray. Vora and Daoutidis17 have shown that the column with the double-feed configuration has achieved a conversion of 76.8%, which is much higher than the equilibrium conversion of 66%, and a product purity of 65%, which is higher than the azeotropic composition of 54%. The conversion and the product purity obtained with the single-feed column configuration are found to be lower than the equilibrium conversion and azeotropic composition. The steady-state data considered for the column with the double-feed configuration are D ) 6.68 mol/s, B ) 7.0825 mol/s, FAc ) 6.9 mol/s, FEth ) 6.865 mol/s, and Lo) 13.51 mol/s. The stochastic optimization based NMPCs of this study are applied to control the ethyl acetate reactive distillation column with its double-feed configuration. Since ethyl acetate is produced significantly and is withdrawn as a product in the distillate stream, controlling the purity of this main product is important in spite of disturbances in the column operation. This becomes the main control loop for GANMPC/SANMPC in which reflux flow rate is used as a manipulated variable to control the purity of ethyl acetate. Since reboiler and condenser holdups act as pure integrators, they also need to be controlled. These become the auxiliary control loops and are controlled by conventional PI controllers. The distillate flow rate is considered as a manipulated variable to control the condenser molar holdup, and the bottom flow rate is used to control the reboiler molar holdup. The tuning parameters17 used for the PI controllers of both the reflux drum and the reboiler holdups are kc ) -0.001 and τI ) 1.99 × 104 s. The SISO control scheme for the column with its double-feed configuration is shown in Figure 2. The input-output data to construct the nonlinear empirical model is obtained by solving the model equations using Euler’s integration with a step size of 2.0 s. A PI controller with a series of step changes in the setpoint of ethyl acetate composition is used for data generation. The input data (reflux flow) is normalized and used along with the outputs (ethyl acetate composition) in model building. The reflux flow rate is constrained within the limits of 20 and 5 mol/s. A total number of 25 000 data sets is considered to develop the model. The model parameters are determined by using the well-known recursive least-squares algorithm,30 the application of which has been discussed elsewhere.31 After evaluating the model structure in eq 6 for different orders of ny and nu, the model with ny ) 2 and nu ) 2 is found to be more appropriate to design and implement the stochastic optimization based NMPCs. The structure of such a model has the form ˆy(k) ) θ0 + θ1y(k - 1) + θ2u(k - 1) + θ3y(k - 1)u(k - 1) + θ4y(k - 2)u(k - 2) + θ5u(k - 1)u(k - 2) (22) The parameters of this model are determined as θ0 ) -0.000774, θ1 ) 1.000553, θ2 ) 0.002943, θ3 ) -0.003828, θ4 ) 0.000766, and θ5 ) -0.000117. This identified model is then used to derive future predictions for the process output by cascading the model to itself, as in eq 9. These model predictions are added with the modeling error, d(k), defined by eq 11, which is considered to be constant for the entire prediction horizon.

Table 1. Effect of Prediction and Control Horizons on the Performance of GANMPC control horizon prediction horizon

1

2

3

4

5

1 2 3 5 10 20 30 50

94.8940 20.9834 17.3545 16.6682 16.5520 16.6265 16.7322 17.0724

21.5340 18.9550 16.9950 16.6249 16.7052 16.8056 17.1443

21.4920 16.6620 16.6710 16.7648 16.8750 17.2120

16.6750 16.6870 16.8168 16.9110 17.2120

16.6980 16.6987 16.8510 17.0015 17.3125

Table 2. Effect of Prediction and Control Horizons on the Performance of SANMPC control horizon 3

prediction horizon

1

2

1 2 3 5 10 20 30 50

94.6020 20.7010 17.2495 16.6218 16.5355 16.6006 16.7108 17.0440

20.9520 17.3350 16.6699 16.5980 16.6440 16.7750 17.0790

17.2790 16.6536 16.5938 16.6625 16.8273 17.1294

4

5

16.7391 16.6569 16.7066 16.8420 17.2038

16.7380 16.6464 16.7053 16.8827 17.1706

The weightings λ and γ in the objective function, eq 10, are set as 1.0 × 107 and 7.5 × 104, respectively. The penalty parameter, µ, in eq 13 is assigned as 1.0 × 105. The cost function used in GANMPC is the penalized objective function, eq 14, based on which the fitness function, eq 1 in GA search is computed. The incremental input, ∆u in GA search is constrained within the limits -0.0025 and 0.0025. The actual input, u, involved with the optimization scheme is a normalized value and is constrained between 0 and 1. A crossover probability of 0.8 and a mutation probability of 0.05 are employed for genetic operators. Roulette wheel selection, single-point crossover, and bitwise mutation are used for GA implementation. Each input variable used in GA is coded as a substring of length 10. The length of the string depends on the control horizon and increases with the increase of input variables. The size of the population is considered equal to the string length. The allowable generation number is fixed as 100. The objective function and the constraints used in SANMPC are the same as in GANMPC. The energy function at each instant is evaluated as the objective function, eq 10. The initial temperature T in SANMPC is chosen as 500, and the number of iterations at each temperature is set as 250. The temperature reduction factor, R in eq 19, is set as 0.5. The control input determined by the stochastic optimizer is denormalized and implemented on the process. A sample time of 2 s is considered for implementation of the controllers. The performances of GANMPC and SANMPC are evaluated by applying for the servo and regulatory control of the ethyl acetate reactive distillation column. The effect of the length of prediction and control horizons of these controllers is expressed by means of integral squared error (ISE) as shown in Tables 1 and 2. These ISE results are evaluated for a period of 30 h simulated operation by considering a step change in ethyl acetate composition setpoint from 0.6827 to 0.75 at time zero with no changes in load variables. From the ISE results of Tables 1 and 2, it is observed that increasing the prediction horizon with a fixed control horizon gradually improves the control performance up to a certain length of prediction horizon and afterward shows a marginal degradation. The GANMPC with a prediction horizon of around 10 and a control horizon of around 1 or 2 presents better performance. The SANMPC with a prediction horizon of


around 10 and a control horizon of around 1-3 provides enhanced performance. A similar trend is observed in the ISE results of GANMPC and SANMPC for a step change in reboiler heat input. However, these ISE results indicate almost equal performance of GANMPC and SANMPC. The nonlinear model predictive controllers with stochastic optimization are also compared with a PI controller, a linear model predictive controller (LMPC), and a nonlinear model predictive controller (NMPC) based on sequential quadratic programming (SQP). The PI controller17 with the parameters kC ) 10.0 and τI ) 1.99 × 104 s is used to control the ethyl acetate composition in the reactive distillation column. The LMPC formulation of this work is based on an auto-regressive exogeneous (ARX) model structure that relates the predicted plant output with the past plant input-output data using the relation

Y(k) ) FX(k) + GU(k) + Ee(k)

(30)

Y(k) ) [yp(k + 1)...yp(k + N)]T

(31)

where

X(k) ) [y(k)y(k - 1)...y(k - ny + 1)u(k - 1)... u(k - nu + 1)]T

(32)

U(k) ) [u(k)...u(k + m - 1)] (33) The elements in F, G, and E have the following matrix notation: T

y(k + 1) ) a1y(k) + ... + anyy(k - ny + 1) + b1u(k) + ... + bnuu(k - nu + 1) (23) where y(k) and u(k) are the process and controller outputs at time k, respectively, and y(k + 1) is the one-step-ahead model prediction at time k. The nu and ny are input and output orders of the system, respectively. The ARX model parameters are updated at each sampling time by an adaptive mechanism. The one-step-ahead model prediction in eq 23 has the form ym(k + 1) ) φxm(k)

(24)

φ ) [R1...Rny β1...βnu]

(25)

xm(k) ) [y(k)...y(k - ny + 1) u(k)...u(k - nu + 1)]T

(26)

where

and

with R and β as identified model parameters. The model parameters are adaptively updated using the popular recursive least-squares (RLS) algorithm.30,31 The N-time-steps-ahead predicted output over a prediction horizon is given by yp(k + N) ) R1y(k + N - 1) + ... + Rnyy(k - ny + N) + β1u(k + N - 1) + ... + βnuu(k - nu + N) + e(k)

(27)

where yp(k + N) represents the model predictions for N steps and e(k) is an estimate of the modeling error, which is assumed to be constant for the entire prediction horizon. If the control horizon is m, then the controller output after m time steps can be assumed to be constant. An internal model is used to eliminate the discrepancy between model and process outputs, e(k), at each sampling instant: e(k) ) y(k) - ym(k)

(28)

where ym(k) is the one-step-ahead model prediction at time (k - 1). Back-substitutions transform the prediction model equations into the following form: yp(k + N) ) fN,1y(k) + ... + fN,nyy(k - ny + 1) + fN,ny+1u(k - 1) + ... + fN,ny+nu-1u(k - nu + 1) + gN,1u(k) + ... + gN,mu(k + m - 1) + eNe(k) (29) The elements f, g, and e are recursively calculated using the parameters R and β of eq 25. The above equations can be written in a condensed form as

The notation Y(k) represents the model predictions over the prediction horizon, X(k) is a vector of past plant and controller outputs, and U(k) is a vector of future controller outputs. If the coefficients of F, G, and E can be determined, then the transformation is completed. The number of columns in F is determined by the ARX model structure used to represent the system, whereas the number of columns in G is determined by the length of the control horizon. The number of rows is fixed by the length of the prediction horizon. The optimal control input sequence can be generated by considering a cost function of the form J ) [yp(k+i)-w(k+i)]2 + [u(k+i-1)]2 ) [Y(k)-W(k)]T[Y(k) - W(k)] + U(k)TU(k) (34) where W(k) is a setpoint vector over the prediction horizon W(k) ) [w(k + 1)...w(k + N)]T (35) The minimization of the cost function, J, gives the optimal controller output sequence: U(k) ) [GTG + γI]-1GT[W(k) - FX(k) - Ee(k)] (36) The vector U(k) generates a control sequence over the entire control horizon. However, the first component of U(k) is actually implemented and the whole procedure is repeated again at the next sampling instant using the latest measured information. The predictive model in LMPC is assigned with the input and output orders as nu ) 2 and ny ) 2. The diagonal elements of the initial covariance matrix, P(0) in the RLS algorithm, are selected as 10.0, 1.0, 0.01, 0.01, respectively. The forgetting


Figure 3. Output and input profiles for step increase in ethyl acetate composition.

factor used in recursive least-squares is chosen as 5.0. The tuning parameter γ in the control law is set as 0.115 × 10-6. The LMPC is evaluated toward the effect of different prediction and control horizons. It is observed that there is no significant change in the results of LMPC for the prediction horizons in the range of 1-10 with the control horizons of 1-3. Increasing the prediction and control horizons beyond these ranges gradually degraded the controller performance. It is observed that the LMPC with a prediction horizon of around 5 and a control horizon of 2 has shown reasonably better control performance. The stochastic optimization based NMPCs and LMPC with the better selected prediction and control horizons are evaluated for their servo and regulatory performance, and the results of these controllers along with those of the PI controller are presented in Figures 3–6. Figures 3 and 4 compare the input and output profiles of GANMPC and SANMPC with those of LMPC and PI controllers for a step change in the ethyl acetate composition setpoint from 0.6827 to 0.75. The results in Figure 5 represent a 20% step decrease in acetic acid feed flow rate, and the results of Figure 6 correspond to a 20% step decrease in ethanol feed flow rate. Figure 7 compares the responses of SANMPC with those of a PI controller and LMPC for a 20% step increase in reboiler heat load. These results exhibit the superior performance of GANMPC and SANMPC over a conventional PI controller. These results also indicate the enhanced performance of NMPCs over LMPC. Figure 8 compares the responses of GANMPC and SANMPC against LMPC for tracking multiple step changes in the setpoint of the controlled variable. These results indicate the better setpointtracking efficiency of NMPCs. The ISE values of stochastic optimization based NMPCs and LMPC are also evaluated for different load disturbances by considering a 20% step change in each disturbance condition, and the results obtained for a period of 30 h duration are shown in Table 3. These results indicate the better regulatory performance of stochastic optimization based NMPCs over LMPC. Though NMPCs provide better results than LMPC, the execution of LMPC requires less computational effort. The GANMPC and SANMPC are also found to be effective in

Figure 4. Output and input profiles for step increase in ethyl acetate composition.

Figure 5. Output and input profiles for step decrease in acetic acid feed flow rate.

tracking setpoint changes in the presence of load disturbances. The results thus show the stability and robustness of stochastic optimization based NMPCs toward the control of reactive distillation column, but only a marginal difference is observed between the results of GANMPC and SANMPC. In stochastic optimization, the problem formulation along with the constraints is simple and straightforward and can be easily extended to many prediction and control horizons.


Minimize f(x)

Subject to h(x) ) 0

g(x) e 0 (37) where f(x) represents the objective function and h(x) and g(x) are the equality and inequality constraints, respectively. SQP forms a quadratic programming approximation to a nonlinear programming problem defined in eq 37 and solves the following quadratic programming subproblem: 1 Minimize ∇ f(x)T∆x + ∆xT∇2f(x)∆x 2 Subject to h(x) + ∇ h(x)T∆x ) 0 g(x) + ∇ g(x)T∆x e 0 (38) The process output predictions of eq 12 are incorporated in the objective function in eq 10 to compute the control sequence at every sampling instant subject to the following equality and inequality constraints: ˆyp(k + i) ) 0, i ) 1 to N 0 e u(t) e 1

0 e u(t + i) e 1 (39) Slack variables are introduced for the inequality constraints, and the subproblem in eq 39 is transformed to Lagrangian form to derive the following Kuhn-Tucker optimality conditions: Qx + ATλ + CTV + c ) 0 Figure 6. Output and input profiles for step decrease in ethanol feed flow rate.

Ax ) b Cx + s ) d SVe ) 0

(s, V) g 0 (40) where c ) 3f(x), A ) 3h(x), C ) 3g(x)T, b ) -h(x), and d ) -g(x). The λ and V are the multipliers of equalities and inequalities, and S and V are diagonal matrices of s and V, respectively. The QP subproblem is solved by introducing a central path concept and defining a complementary gap and by integrating an interior

Figure 7. Output and input profiles for step increase in reboiler heat load.

The NMPCs with stochastic optimization are further compared with a sequential quadratic programming based nonlinear model predictive controller (SQPNMPC). The controller is designed by combining the identified polynomial input-output model in eq 22 with the SQP optimizer. The design is based on considering a general nonlinear programming problem as

Figure 8. Output responses for multiple changes in ethyl acetate composition setpoint.

6958 Ind. Eng. Chem. Res., Vol. 47, No. 18, 2008 Table 3. Performance of Model Predictive Controllers ISE disturbance

LMPC

SQPNMPC

GANMPC

SANMPC

acetic acid feed flow decrease acetic acid feed flow increase ethanol feed flow decrease ethanol feed flow increase reboiler heat input decrease reboiler heat input increase

0.559167

0.004180

0.003253

0.001671

0.113510

0.006457

0.005327

0.002763

0.570151

0.006511

0.006117

0.003427

0.106454

0.005651

0.001979

0.000988

0.511803

0.334305

0.319454

0.146713

0.691946

0.396545

0.389169

0.271435

point method known as predictor-corrector algorithm within the SQP. A line search is imposed on the algorithm to determine the step length at every instant. The Hessian of the Lagrangian ∇2f(x) is approximated by a positive definite matrix Q, which is considered as an identity matrix initially. At every sampling instant, the QP subproblem is solved iteratively for the prediction, centering, and correction steps of the predictor-corrector algorithm. A line search is imposed for the new set of values, and the Hessian of the Lagrangian is updated using the Broyden-Fletcher-GoldfarbShanno (BFGS) method. After attaining convergence, the first of the input variables of the specified control horizon is implemented on the process. More details concerning the implementation of interior -point SQP are available in the literature.32 The weightings λ and γ in the objective function in eq 10 are considered to be unity for SQPNMPC. The controller is evaluated toward the effects of different identical prediction and control horizons. The SQPNMPC with a prediction horizon of 2 and a control horizon of 2 is found to be effective over other horizons. The ISE values of the SQPNMPC are evaluated toward the effect of different load disturbances, each with a 20% step change, while keeping the composition setpoint at its base value. The ISE results shown in column 3 of Table 3 indicate the better performance of SQPNMPC over LMPC. From these results, it is observed that the GANMPC and SANMPC provides slightly better performance over SQPNMPC. It is found that there is no significant change in the results of SQPNMPC for prediction and control horizons beyond 2, but the computational burden increases considerably with the increase of these horizons. Figure 9 shows improved results of SANMPC over SQPNMPC for the cases of (i) a step change in ethyl acetate composition from 0.6827 to 0.75 while keeping the load changes at their base values and (ii) a 20% step decrease in acetic acid feed flow rate while keeping the ethyl acetate composition setpoint at its base value. The ISE value of SQPNMPC that is evaluated for a period of 30 h duration for a step increase in ethyl acetate composition setpoint is found to be 769.59. This value appeared to be significantly higher to those of GANMPC and SANMPC evaluated for the same case with their better selected prediction and control horizons as given among the values in Tables 1 and 2. These results thus indicate the better performance of stochastic optimization based NMPCs over SQPNMPC. The computational efficiencies of stochastic optimization based NMPCs are evaluated in terms of CPU time. The computational times of GANMPC and SANMPC for different prediction and control horizons on a X86 computer with 127 MB RAM using C Programming are shown in Table 4. These times correspond to a complete simulation period of 10 h with a sample time of 2 s when a step change in ethyl acetate composition from 0.6827 to 0.75 is

Figure 9. Output responses for step increase in ethyl acetate composition setpoint and step decrease in acetic acid feed flow rate. Table 4. Computational Efficiency

prediction horizon

control horizon

time of execution of GANMPC (min)

10 10 10 10 10 1 2 3 4 5 6 7 8 9

1 2 3 4 5 1 1 1 1 1 1 1 1 1

66.15 92.11 118.27 144.43 170.06 51.26 52.16 54.37 55.30 56.73 59.04 61.33 63.10 64.43

time of execution of SANMPC (min) 21.41 22.73 24.16 25.83 27.31 5.62 7.91 9.53 11.37 13.73 15.17 16.67 18.19 19.76

considered. The CPU times of SQPNMPC evaluated for the same operation are found to be 31.74, 67.57, and 151.62 min, respectively, for identical prediction and control horizons of 2, 3, and 4. From these results, it is observed that the execution time required for SQPNMPC is higher than that for SANMPC and lower than that for GANMPC. However, the ISE results of regulatory and servo operations indicate improved performance of stochastic optimization based NMPCs over SQPNMPC. GANMPC involves a higher number of tuning parameters and random numbers, thus requiring more tuning effort. Although GANMPC and SANMPC are found to exhibit robust performances, SANMPC is preferred for the control of a reactive distillation column because of its better performance, easier tuning, and lower computational effort. 7. Conclusions Stochastic methods are being increasingly used for solving complex, nonlinear, and nonconvex optimization problems. In this study, stochastic optimization based nonlinear model predictive controllers (NMPCs), namely, GANMPC and SANMPC, are proposed and applied to control a reactive


distillation column. The results demonstrate the better performance of the stochastic optimization based NMPCs over a conventional PI controller, a linear model predictive controller (LMPC), and a nonlinaer model predictive controller (NMPC) based on sequential quadratic programming (SQP). The stochastic optimizers, because of their stochastic nature, may not provide a rigorous guaranty for global optimality, but the problem formulation along with the constraints is simple and straightforward. With the availability of low-cost computing, these methods can be applied to solve even larger optimization problems. Although both the GANMPC and SANMPC of this study are found to exhibit almost equal performance, SANMPC involves a smaller number of tuning parameters and requires less execution time. Hence, SANMPC is better suited for the control of reactive distillation. The method proposed in this study can be easily extended to multiple input-multiple output (MIMO) systems. However, its implementation to such systems requires additional computational effort because of the involvement of multistep predictive models and control horizons. However, with the enormous computation power available in present times, the application of this methodology to such systems will not pose any difficulty. Appendix A: Mathematical Model of Reactive Distillation The equations representing the process are given as follows. Total mass balance Total condenser: dM1 ) V2 - (L1 + D) + ∆R1 dt

(A1)

Reboiler: dEn ) Ln-1hn-1 - VnHn - Lnhn + Qn dt LeVel of liquid on the tray: Lliq )

dMn ) Ln-1 - Vn - Ln + ∆Rn dt Component mass balance Total condenser:

If Lliq < hweir then Ln ) 0

d(M1xi,1) ) V2yi,2 - (L1 + D)xi,1 + Ri,1 dt

Ln ) 1.84

Plate j: d(Mjxi,j) ) FVjyfi,j + Vj+1yi,j+1 + FLjxfi,j + Lj-1xi,j-1 - Vjyi,j dt Ljxi,j + Ri,j (A5) Reboiler: d(Mnxi,n) ) Ln-1xi,n-1 - Vnyi,n - Lnxi,n + Ri,n dt Total energy balance Total condenser: dE1 ) V2H2 - (L1 + D)h1 + Q1 dt

LweirFav liq (L - hweir)3⁄2 MWav

(A6)

(A7)

Plate j: dEj ) FVjHfj + Vj+1Hj+1 + FLjhfj + Lj-1hj-1 - VjHj - Ljhj + Qj dt (A8)

(A12)

Mole fraction normalization: NC

NC

∑ ∑y )1 xi )

i)1

i

(A13)

i)1

VLE Calculations. For the column operation under moderate pressures, the VLE equation assumes the ideal gas model for the vapor phase, thus making the vapor-phase activity coefficient equal to unity. The VLE relation is given by yiP ) xiγiPisat (i ) 1, 2, ... , NC)

(A14)

The liquid-phase activity coefficients are calculated using UNIFAC method.33 Enthalpies Calculations. The relations for the liquid enthalpy h, the vapor enthalpy H, and the liquid density F are as follows: h ) h(P, T, x) H ) H(P, T, y) Fliq ) Fliq(P, T, x) Notation

(A4)

(A11)

else

(A2)

(A3)

(A10)

Flow of liquid oVer the weir:

Plate j: dMj ) FVj + Vj+1 + FLj + Lj-1 - Vj - Lj + ∆Rj dt Reboiler:

MnMWav AtrayFav

(A9)

Atray ) tray area, m2 B ) bottom flow rate, mol s-1 C ) concentration, mol m-3 Ck ) catalyst concentration, vol % D ) distillate flow rate, mol s-1 E ) total enthalpy of liquid on plate, kJ FL ) liquid feed flow rate on plate, mol s-1 FV ) vapor feed on plate, mol s-1 FAc ) acetic acid feed flow rate, mol s-1 FEth ) ethanol feed flow rate, mol s-1 H ) molar enthalpy of vapor stream, kJ mol-1 h ) molar enthalpy of liquid stream, kJ mol-1 k1 ) reaction rate constant, m3 mol-1 s-1 hweir ) weir height, m KC ) constant of reaction equilibrium L ) molar liquid flow rate, mol s-1 Lweir ) weir length, m Lliquid ) liquid level on tray, m M ) molar holdup on plate, m MWav ) average molecular weight, g mol-1 P ) pressure on plate, pascal Q ) heat exchange, kJ R ) number of moles reacted, mol s-1 r ) rate of reaction, mol/s · m3 Fav ) average density, g m-3 T ) temperature, K V ) molar vapor flow rate, mol s-1

(A15)

6960 Ind. Eng. Chem. Res., Vol. 47, No. 18, 2008 x ) mole fraction in liquid phase y ) mole fraction in vapor phase

Literature Cited (1) Qin, J.; Badgwell, T. An overview of industrial model predictive control technology. In Vth International Conference on Chemical Process Control; Kantor, J. C., Garcia, C. E., Carnhan, B., Eds.; AIChE Symposium Series 93; AIChE: New York, 1997; pp 232-256. (2) Eaton, J. W.; Rawlings, J. B. Model predictive control of chemical processes. Chem. Eng. Sci. 1992, 47, 705–720. (3) Venkateswarlu, Ch.; Gangiah, K. Constrained generalized predictive control of unstable nonlinear processes. Trans. Inst. Chem. Eng. 1997, 75, 371–376. (4) Wright, G. T.; Edgar, T. F. Nonlinear model predictive control of a fixed-bed water-gas shift reactor: An experimental study. Comput. Chem. Eng. 1994, 18, 83–102. (5) Ricker, N. L.; Lee, J. H. Nonlinear model predictive control of the Tennessee Eastman challenging process. Comput. Chem. Eng. 1995, 19, 961–981. (6) Hernandez, E.; Arkun, Y. On the global solution of nonlinear model predictive control algorithms that use polynomial models. Comput. Chem. Eng. 1994, 18, 533–536. (7) Venkateswarlu, Ch.; Venkat Rao, K. Dynamic recurrent radial basis function network model predictive control of unstable nonlinear processes. Chem. Eng. Sci. 2005, 60, 6718–6732. (8) Camacho, E. F.; Bordons, C. Model PredictiVe Control in the Process Industry; Springer Verlag: Berlin, Germany, 1995. (9) Ahn, S. M.; Park, M. J.; Rhee, H. K. Extended Kalman based nonlinear model predictive control of a continuous polymerization reactor. Ind. Eng. Chem. Res. 1999, 38, 3942–3949. (10) Chen, H.; Allgower, F. A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. Automatica 1998, 34, 1205–1218. (11) Faber, R.; Jockenhovel, T.; Tsatsaronis, G. Dynamic optimization with simulated annealing. Comput. Chem. Eng. 2005, 29, 273–290. (12) Li, P.; Lowe, K.; Arellano-Garcia, K. H.; Wozny, G. Integration of simulated annealing to a simulation tool for dynamic optimization of chemical processes. Chem. Eng. Process. 2000, 39, 357–363. (13) Degarmo, J. L.; Parulekar, V. N.; Pinjala, V. Consider reactive distillation. Chem. Eng. Prog. 1992, 88 (3), 43–50. (14) Sneesby, M. G.; Tade, M. O.; Datta, R.; Smith, T. N. ETBE synthesis via reactive distillation. 2. Dynamic simulation and control aspects. Ind. Eng. Chem. Res. 1997, 36, 1870–1881. (15) Al-Arfaj, M. A.; Luyben, W. L. Design and control of olefin metathesis reactive distillation column. Chem. Eng. Sci. 2002, 57, 715– 733. (16) Kumar, A.; Dauotidis, P. Modeling, analysis and control of ethylene glycol reactive distillation column. AIChE J. 1999, 45, 51–68.

(17) Vora, N.; Daoutidis, P. Dynamics and control of ethyl acetate reactive distillation column. Ind. Eng. Chem. Res. 2001, 40, 833–849. (18) Shopova, E. G.; Vaklieva-Bancheva, N. G. BASIC-A genetic algorithm for engineering problems solution. Comput. Chem. Eng. 2006, 30, 1293–1309. (19) Hanke, M.; Li, P. Simulated annealing for the optimization of batch distillation process. Comput. Chem. Eng. 2000, 24, 1–8. (20) Goldberg, D. E. Genetic Algorithms in Search, Optimization and Machine Learning; Addison-Wesley: Reading, MA, 1989. (21) Chen, L.; Nguang, S. K.; Chen, X. D.; Li, X. M. Modelling and optimization of fed-batch fermentation processes using dynamical neural networks and genetic algorithms. Biochem. Eng. J. 2004, 22, 51–61. (22) Onnen, C.; Babuska, R.; Kaymark, J. M.; Sousa, J. M.; Verbruggen, H. B.; Isermann, R. Genetic algorithms for optimization in predictive control. Control Eng. Pract. 1997, 5 (10), 1363–1372. (23) Kirkpatrick, S.; Gelatt, C. D., Jr.; Veccchi, M. P. Optimization by simulated annealing. Science 1983, 220, 671–680. (24) Vanderbilt, D.; Louie, S. G. A Monte Carlo simulated annealing approach to optimization over continuous variables. J. Comput. Phys. 1984, 56, 259–271. (25) Dolan, W. B.; Cummings, P. T.; Le Van, M. D. Process optimization via simulated annealing: Application to network design. AIChE J. 1989, 35, 725–736. (26) Haber, R.; Unbehauen, H. Structure identification of nonlinear dynamical systemssA survey on input/output approaches. Automatica 1990, 26, 651–677. (27) Morningred, J. D.; Paden, B. E.; Seborg, D. E.; Mellichamp, D. A. An adaptive nonlinear predictive controller. Chem. Eng. Sci. 1992, 47, 755– 762. (28) Hernandez, E.; Arkun, Y. Control of nonlinear systems using polynimial ARMA models. AIChE J. 1993, 39, 446–460. (29) Alejski, K.; Duprat, F. Dynamic simulation of the multicomponent reactive distillation. Chem. Eng. Sci. 1996, 51, 4237–4252. (30) Goodwin, G. C.; Sin, K. S. AdaptiVe Filtering, Prediction and Control; Prentice Hall: Englewood Cliffs, NJ, 1984. (31) Venkateswarlu, Ch.; Naidu, K. V. S. Adaptive fuzzy model predictive control of an exothermic batch chemical reactor. Chem. Eng. Commun. 2001, 186, 1–23. (32) Albuquerque, J.; Gopal, V.; Staus, G.; Biegler, L. T.; Ydstie, B. E. Interior point SQP strategies for large-scale, structured process optimization problems. Comput. Chem. Eng. 1999, 23, 543–554. (33) Smith, J. M.; Van Ness, H. C.; Abbot, M. M. A Text Book on Introduction to Chemical Engineering Thermodynamics, 5th ed.; McGraw Hill International: New York, 1996.

ReceiVed for reView July 18, 2007 ReVised manuscript receiVed March 24, 2008 Accepted May 29, 2008 IE070972G

Nonlinear Model Predictive Control of Reactive Distillation Based on

Recommend Documents