Economic Oriented NMPC for an Extractive Distillation Column Using

Jun 2, 2015 - Finally, we use the index hybrid DAE model in the economic oriented NMPC (EO-NMPC) of an extractive distillation column considering diff...
0 downloads 9 Views 887KB Size
Article pubs.acs.org/IECR

Economic Oriented NMPC for an Extractive Distillation Column Using an Index Hybrid DAE Model Based on Fundamental Principles Federico Lozano Santamaría and Jorge M. Gómez* Grupo de Diseño de Productos y Procesos, Departamento de Ingeniería Química, Universidad de los Andes, Carrera 1 No. 18A-10, Bogotá 111711, Colombia ABSTRACT: Nonlinear model predictive control (NMPC) is an advanced control strategy that uses a rigorous dynamic model of the process, usually described by differential-algebraic equations (DAE), to predict its future behavior. The economic optimization of the process can be included directly in the objectives of the NMPC controller, which make it attractive for application in highly nonlinear and energy-intensive processes that are subject to fluctuating operational conditions. Currently, the implementation of NMPC controllers with large-scale dynamic models is limited by the computational difficulty of solving the associated dynamic optimization problem at each period of time. In this work, we propose an index two DAE model based on fundamental principles that uses its equivalent reduced index one model to define consistent initial conditions (index hybrid DAE) with the aim of representing the dynamics of a distillation column. We demonstrate that this model describes the same physical behavior of its equivalent index one reduced model and that it decreases significantly the computational complexity of the online optimization associated with the NMPC controller. Finally, we use the index hybrid DAE model in the economic oriented NMPC (EO-NMPC) of an extractive distillation column considering different disturbances over the inlet conditions. Also, we compare the EO-NMPC with a classic PI (proportional and integral) controller to show that the first one has a better dynamic performance and improves the economic profit of the process.

1. INTRODUCTION Model predictive control (MPC) is an advanced feedback control strategy that solves an optimal control problem (OCP) at certain periods of time. It uses a detailed model of the process to predict its behavior in a future horizon in order to determine the control sequence that minimizes an objective function that usually is a tracking trajectory function. The control law is defined by the first element of the control sequence obtained as a solution of the OCP at each sampling time, which means that only the first element of the control sequence predicted is feedback to the system. The main advantages of this strategy over traditional controllers like PID (Proportional, Integral and Derivative) are its implementation does not need wide knowledge of control theory, its tuning is relatively intuitive, inherent compensation of dead-time and time delay phenomena and it can handle constraints and the multivariable problem directly.1 MPC has been widely used in the process industry due to their advantages and easy implementation;2,3 however, they normally use a linear model of the process (nonlinear) that is obtained by direct measurements of the system at its operation point under small disturbances. Thus, the dynamic behavior predicted in the future horizon does not match the real process when it is displaced far from its normal operation point. On the other hand, the use of linear models in MPC schemes simplified the OCP problem and guarantees the global optimum; also, there are specialized algorithms to solve linear and quadratic problems (e.g., CPLEX, Gurobi). Nonlinear models can be used in MPC strategies (nonlinear model predictive control, NMPC) so that the prediction of the dynamic behavior is closer to the real process, but it increases the complexity of the OCP that has to be solved at each sampling time. Using nonlinear models in MPC controllers defines an © 2015 American Chemical Society

optimization problem that may not be convex, thus the global optimum cannot be guaranteed; also, it requires the implementation of specialized nonlinear programming (NLP) algorithms. There is a compromise between the quality of the prediction over the future horizon and the difficulty to solve the OCP at each sampling time. Simplified or empirical nonlinear models have been used in the NMPC of different processes; for example, a nonlinear autoregressive with exogenous input (NARX) model is adjusted and applied to control a distillation column4 and an autoregressive mean average (ARMA) model is applied to a reactive distillation column used to produce ethyl acetate.5 Another alternative is to use a linearized model of a rigorous nonlinear description of the process; for example, Purohit et al.6 applied a successive linearization technique of a rigorous nonlinear DAE system describing the dynamics of a reactive distillation column and it shows a reduction of the computational time and a similar response of the controller compared with the NMPC using the nonlinear DAE model. Nevertheless, the prediction obtains with these models is still limited to certain regions due to their empirical characteristics or to the small regions around the linearization point. On the other hand, models based on fundamental principles can predict a better behavior of the dynamic process over all the operation conditions, even when the system is far away from its nominal operation point. However, the phenomenological equations that describe the process (e.g., thermodynamic relations, material and energy balances) are defined as large-scale differential algebraic Received: Revised: Accepted: Published: 6344

March 5, 2015 May 3, 2015 June 2, 2015 June 2, 2015 DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research equations (DAE) increasing the complexity of the NLP problem that has to be solved at each sampling time. Also, computational delays may be introduced. Distillation is an energy-intensive operation and it is also sensitive to disturbances that modify its dynamic behavior, thus any operational improvement produces significant economical savings. Processes economic optimization has traditionally been addressed in a two-layered structure. The upper layer defines the operation point of the system maximizing the economic profit of the process via real time optimization (RTO) using a detailed steady sate model. The lower layer receives the set-point and reference trajectories obtained by RTO to implement a feedback control strategy (e.g., MPC, PID).7 However, this two-layered structure is not suitable for cyclic processes (e.g., pressure swing adsorption) or processes subject to constant disturbances; it may also lead to suboptimal solutions or even to unreachable setpoints because the upper layer does not consider the dynamic behavior nor the disturbances of the process.8,9 Using (N)MPC, the economic objective of the process can be handled directly in the control strategy, integrating the two-layered structure and improving the process performance, which may produce significant savings for the operation of distillation columns. Economic oriented NMPC (EO-NMPC) changes the standard tracking objective function of the (N)MPC by a general economic function and it has been applied in many works;10,11 however, these works add a regularization term to the objective function because the original change does not guarantee neither the stability nor the robustness of the controller.12 In this work, the EO-NMPC using a model based on fundamental principles is applied to an extractive distillation column. The nonlinear model is an index two DAE with consistent initial conditions determined from its equivalent index one reduced model. In the next section, two different dynamic models of the distillation column are described, both based on the same principles. Section 3 describes the index hybrid methodology proposed that uses the index one reduced model to define consistent initial conditions and then integrate directly the index two DAE. Section 4 presents the EO-NMPC formulation adding regularization terms, the description of the dynamic optimization problem and the discretization used to transform the differential equations into a large-scale system of algebraic equations. Section 5 describes the ethanol dehydration that is the case of study selected. The results show that the index hybrid DAE model is equivalent to the index one reduced DAE and the EO-NMPC performance are presented in Section 6. Section 7 provides conclusions and future challenges.

Figure 1. Distillation column representation.

there is no pressure drop in the reboiler, (5) constant mass holdup in the condenser and reboiler, (6) the vapor hold-up is negligible with respect to the liquid hold-up for all the stages (Mv,n ≈ 0 ∀ n ∈ {1...NT}); the density of the liquid phase is much greater than the density of the vapor phase and (7) the pressure does not vary with time. It helps the problem to be less rigid. Additionally, for this study, it is assumed ideal behavior of the vapor phase and the liquid phase is modeled with NRTL (nonrandom two liquids). Under these assumptions, the model can be described by the following algebraic and differential equations. total mass balance: dML, n dt dML, n dt

⎛ 1⎞ = Vn + 1 − Ln⎜1 + ⎟ = 0, ⎝ R⎠

n=1

= Fn + Ln − 1 + Vn + 1 − Ln − Vn ,

n ∈ {2...NT − 1} dML, n dt

(1)

= Ln − 1 − Ln − Vn = 0,

(2)

n = NT

(3)

partial mass balance:

2. DYNAMIC MODELS FOR A DISTILLATION COLUMN The dynamic model for a distillation column is an extrapolation of the model for a flash separation drum and it is based on mass and energy balances between the process inlets and outlets at each stage (Figure 1). The thermodynamic equilibrium and hydraulic relations are additional constraints needed to define all the variables along the column; for example, tray temperature, tray hold-ups and liquid flow rates. The next subsections present an index two DAE system and its equivalent index one reduced DAE, both based on MESH (material, energy, summation and enthalpy) equations. 2.1. Index Two DAE System (DAE2). The model presented here uses the MESH equations to describe the dynamic behavior of a distillation column. The most important assumptions applied are13 (1) thermodynamic equilibrium at each stage, (2) adiabatic operation, (3) total condenser and partial reboiler, (4)

ML, n

⎛ 1⎞ = yn + 1, i Vn + 1 − xn , iLn⎜1 + ⎟ , ⎝ dt R⎠

dxn , i

i ∈ {1, 2, ..., NC}, n = 1 ML, n

dxn , i dt

(4)

= Fn(zn , i − xn , i) + Ln − 1(xn − 1, i − xn , i)

+ Vn + 1(yn + 1, i − xn , i) − Vn(yn , i − xn , i), i ∈ {1, 2, ..., NC}, n ∈ {2...NT − 1}

ML, n

dxn , i dt

(5)

= xn − 1, iLn − 1 − xn , iLn − yn , i Vn ,

i ∈ {1, 2, ..., NC}, n = NT

(6)

energy balance: 6345

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research d(HL, nML, n) dt

⎛ 1⎞ = Vn + 1HV, n + 1 − Ln⎜1 + ⎟HL, n − Q C , ⎝ R⎠

different studies related to the dynamic or control of distillation columns.16,18,19 It consists in differentiating the left-hand side of the differential energy balance (eqs 7−9) and replace it with an algebraic equation that lead to an algebraic energy balance. The new energy balance (eqs 13−15) is an algebraic equation that defines the vapor flow rate directly for each stage. The term dML,n/dt is known from eqs 1−3. Considering the liquid enthalpy as HL,n = ∑iNC = 1xi,nHL,i,n where HL,i,n is the specific liquid molar enthalpy of each component and it is an algebraic function of temperature, the term dHL,n/dt can be calculated analytically using the chain rule and the derivative of temperature with respect to time for each stage (eq 16). This derivative is obtained from the thermodynamic equilibrium relation under the assumption of ideal gas behavior (eq 10).

(7)

n=1

d(HL, nML, n) dt

= HF, nFn + HL, n − 1Ln − 1 + HV, n + 1Vn + 1

− HL, nLn − HV, nVn , d(HL, nML, n) dt

n ∈ {2...NT − 1}

(8)

= Q R + HL, n − 1Ln − 1 − HL, nLn − HV, nVn , (9)

n = NT

phase equilibrium error:



(yn , i − xn , i) = 0,

ML, n

n ∈ {1...NT}

i ∈ {1,2,...,NC}

Pnsat, i γn , i Pn

ML, n

n ∈ {2...NT − 1}

dHL, n dt

⎛ 1⎞ HL, n = Vn + 1HV, n + 1 − Ln⎜1 + ⎟HL, n ⎝ dt R⎠

dML, n

n=1 +

dML, n dt

(13)

HL, n = HF, nFn + HL, n − 1Ln − 1

+ HV, n + 1Vn + 1 − HL, nLn − HV, nVn ,

(11)

n ∈ {2...NT − 1} (14)

liquid flow rate: Ln = f (ML, n , Vn),

+

− Q C,

xn , i ,

i ∈ {1, 2, ..., NC}, n ∈ {1...NT}

dt

(10)

thermodynamic equilibrium: yn , i = K n , ixn , i =

dHL, n

ML, n

(12)

The hydraulic equations used to model the liquid flow rate at each stage of the distillation column are obtained from correlations and empirical relations. The equations were taken from Monteil 196514 and they consider the liquid flow rate as a function of the liquid hold-up, the vapor flow rate and the geometry of the tray (eq 12). Equations 1−12 lead to a DAE system, with differential variables ML,n, xn,i and ML,nHL,n where i ∈ {1,2,...,NC} and n ∈ {1...NT}. The number of differential equations is NT(NC + 2) − 2 and the number of algebraic equations is. This model has two degrees of freedom (e.g., the reflux ratio and the reboiler duty) that, in the case of a control problem, they correspond to the manipulated variables. In this case, the model is an index two DAE because the vapor flow rate of each stage is an algebraic variable and it is not represented in the subset of algebraic equations; thus, the change of this variable with respect to time cannot be known directly when a method to solve index one problems is applied. Nevertheless, the number of equations and variables is reduced by a factor of NT(NC + 3) compared with a model that considers the vapor hold-up,15 also eliminating the vapor phase dynamics reduces the stiffness of the model;16 as a result, this simplified model that does not consider the vapor hold-up can be solved easier in computational terms. Defining consistent initial conditions is the principal difficulty to solve high-index DAE systems. When the initial values of the differential variables are set, the subsystem of algebraic equations is singular for high index DAE systems. If arbitrary initial conditions are set for those algebraic variables that cannot be defined by solving the subsystem of algebraic equations, the numerical solution may present impulse behavior at the beginning of the integration and it will be displaced from the true manifold.17 2.2. Index One Reduced DAE System (DAE2r). To overcome the problem of defining consistent initial conditions for the previous DAE system (Section 2.1), an index reduction technique of differentiation and substitution has been applied in

dHL, n dt

+

− HV, nVn , dTn = dt

dML, n dt

HL, n = Q R + HL, n − 1Ln − 1 − HL, nLn

n = NT

(15)

dx n , i N ⎡ N ∂γ (xn , j , Tn) dxn , j ⎤ NC 1 − P ∑i =C1 ⎢xn , iPnsat, i ∑ j =C1 n , i∂x ⎥ − ∑i = 1 K n , i dt dt ⎦ ⎣ n n,j , dPnsat NC 1 ⎛ NC ,i ⎞ sat ∂γn , i(xn , j , Tn) ⎜∑ ⎟ ∑ + γ x P x n , i n , i n , i n , i dTn ⎠ i=1 Pn ⎝ i = 1 ∂Tn

n ∈ {1...NT }

(

)

(16)

The equations used in this study to calculate thermodynamic 20 properties (e.g., Psat i , γi) can be found in Aspen-Technology and their respective derivatives can be calculated. When the algebraic energy balance is used instead of the differential balance the variation of the vapor flow rate of each stage can be known in each instant of time, this reduces the problem to an index one DAE with differential variables ML,n and xn,i where i ∈ {1, 2, ..., NC} and n ∈ {1...NT}. Even though, the number of equations of the DAE2 system and the DAE2r system is the same, the index reduction process introduces a significant number of dummy variables and equations during the calculation of all the derivatives required; thus, when the model is translated into an algebraic programming language, the dimension of the model DAE2r is greater than the dimension of the model DAE2. For this reason, in computational terms, it is more difficult to solve the DAE2r system than the DAE2 system. The DAE2r system is index one; therefore, it has the advantage of defining a nonsingular subsystem of algebraic equations when the initial values of the differential variables are set. This guarantees the definition of consistent initial conditions and an adequate solution of the DAE system.

3. INDEX HYBRID DAE MODEL Another alternative to solve index two DAE systems, besides the index reduction, is defining consistent initial conditions and then integrating the model with an adequate integration algorithm (e.g., orthogonal collocation, multistep backward differentiation formula, Runge−Kutta methods).21 Consistent initial conditions 6346

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research are those that satisfy all the algebraic equations of the index two DAE and also the hidden constraints of the system, which means that they have to satisfy the additional algebraic equations obtained after the index of the model has been reduced to one using any reduction technique.22 Defining consistent initial conditions allows the direct integration of index two DAE systems because this point belongs to an invariant set defined by all the algebraic constraints, including the hidden constraints, thus during the integration these constraints will not be violated.23,24 A steady-state solution of the DAE system always guarantees consistent initial conditions, which is useful for certain types of studies such as start-up and shut-down of processes,25 but for process control under continuous disturbances, this is not feasible. However, as the index two DAE is equivalent to the index one model obtained after applying an index reduction technique, this index one model can be used to define consistent initial conditions for the index two DAE at any point and continue the integration with the high index model. This strategy that uses an index hybrid DAE model (DAE2h), integrating the index two DAE and defining consistent initial conditions with its equivalent index one reduced model, was previously tested and proved in the NMPC controller of a flash separation drum.26 The main advantages of the methodology proposed (DAE2h) are (1) the reduction of the number of equations and variables after discretization because the derivatives involved in the index reduction are only calculated at the first point of integration, and (2) the reduction of the computational time used for solving the DAE system: this is a direct implication of the previous advantage. These advantages allow an improvement of the performance of NMPC controllers where usually the computational time is a limitation. However, the direct integration of an index two DAE system proposed with this methodology requires a numerical method with strong stability properties and control of the local integral error to avoid large global errors in the algebraic variables.27 It has been reported that collocation methods with at least three internal points are suitable for the direct integration of index two DAE systems with an acceptable accuracy.27

g (zd(t ), za(t ), u(t )) ≤ 0

where t0 is the initial time, tf is the final time, zd is a vector of differential state variables, za is a vector of algebraic state variables, z0 is a vector of initial conditions and u is a vector of manipulated variables. The objective function is composed by two terms; the first one ϕ(tf) is a terminal cost that is only a function of the states at the final time and the second term is a moving cost that considers all the trajectory of the states and manipulated variables. The DAE models described in Section 2 and Section 3 contain the rigorous description of the process and they correspond to eqs 17b−17d in the dynamic optimization problem formulation. The inequality constraints (eq 17e) usually represent bounds for the states and manipulated variables (polytopic constraints) or specific operational constraints like product purity limits.29 The dynamic optimization problem stated in eqs 17a−17e can be transcribed into a large-scale NLP using a discrete time representation and particular strategies have been developed for its solution. Among these alternatives, the most common are single shooting, multiple shooting and direct transcription. The first two strategies only discretize the manipulated variables and all the DAE system is solved at each iteration of the optimization routine, with the difference that in multiple shooting the time space is divided into intervals and the problem is solved sequentially, whereas in single shooting, the problem is solved directly for all the time space.28 On the other hand, the direct transcription strategy discretizes all the state and manipulated variables transforming the differential equations into equality constraints. Important advantages of this strategy are it solves the DAE system only once at the optimal point, it involves the ability to exploit the sparse structure of the dynamic model represented in the Hessian and constraint matrices using the NLP solver and it transforms the process model in a completely algebraic representation that allows an easier calculation of first and second derivatives from modeling platforms (e.g., GAMS, AMPL).16,30,31 In this work, we implement the direct transcription strategy using orthogonal collocation on finite elements described in the following subsection. 4.2. Discretization of the Differential Equations. The differential equations are transformed into a set of algebraic equations, discretizing all the state and control variables, with the purpose of solving the optimization problem using simultaneous strategies that are compatible with NLP algorithms. This transformation is done using orthogonal collocation on finite elements with a fixed step length of 5 min and using three Radau points. The length of each finite element is equal to the time step of integration (h), note that also the length of each element is equivalent to the sampling time in order to simplify the implementation of the NMPC controller. This type of discretization is chosen because: (1) it is widely used to solve index two problems,32 (2) it is a numerical method with Lstability, which means that when the integration step tends to zero the stability region tends to infinity, (3) it is capable of stabilizing high index DAE systems. The previous properties hold only if the initial conditions are consistent, otherwise the numerical method can fail at the first steps of integration.23,31 However, the method is less accurate in the algebraic variables for index two problems.28,33 Applying this method, the differential equations become a weighting of the variables using the derivative of an orthogonal polynomial (Lagrange’s polynomial) (eq18). These polynomials are represented in eq 19.

4. EO-NMPC PROBLEM FORMULATION In this section, we first describe the general formulation and solution strategies of a dynamic optimization problem of the kind that is solved at each sampling time in a NMPC scheme. Then, we present the EO-NMPC formulation, its characteristics and implementation. 4.1. Dynamic Optimization Strategy. In a NMPC scheme, at each instant of time (tk), an OCP has to be solved in order to obtain the control signal that will be applied to the process. Assuming that the current states of the process can be measured or estimated, the OCP is defined as a general dynamic optimization problem described by eqs 17a−17e.28 min J = ϕ(tf ) +

s. t.

∫t

t0

φ(zd(t ), za(t ), u(t ))dt

f

dzd = f (zd(t ), za(t ), u(t )) dt

(17a)

(17b)

zd(0) = z 0

(17c)

h(zd(t ), za(t ), u(t )) = 0

(17d)

(17e)

6347

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research K

∑ zdjk′ k ′= 0

Np

dlk ′(τk) = hjf (t jk , zdjk , zajk), dτ

min JEO‐NMPC =

i=0

k = 1, 2, ..., K , j = 1, 2, ..., N K

lk(τ ) =

∏ n = 0, ≠ k

τ − τn , τk − τn

(22)

The main advantage of the EO-NMPC is that it integrates the two-layered structure used in the classic process economic optimization, which means that the set-points and control trajectories are determined simultaneously solving the OCP. Also, the EO-NMPC is suitable for process with constant disturbances, cyclic operations or batch operations.36 However, replacing the objective function of the NMPC formulation by a general economic function does not guarantee the stability and robustness of the controller because it may not be a bounded function on the infinite horizon or a monotonic function.12 If the objective function is a Lyapunov function, robustness and nominal stability are guaranteed for EO-NMPC applied to cyclic processes as it was demonstrated by,37 however their analysis is immediately extended and validated to general EONMPC considering only one period. Specifically, the moving cost of the objective function has to satisfy Lipschitz continuity, weak controllability; it has to be greater than a monotonic unbounded function (K∞) of the state variables and the optimization problem has to have a local isolated optimum.34,37,38 A modification of the objective function that guarantees the previous conditions for stability is to add a regularization term to the original moving cost and it is equivalent to the original NMPC tracking objective function:37

k = 1, 2, ..., K (19)

K

∑ zdjklk(1),

j = 1, 2, ..., N (20)

k=0

4.3. EO-NMPC Description. The conventional (N)MPC aims to minimize, at each time period k, the quadratic error of the controlled variables with respect to a set-point over a prediction horizon (N p ), whereas the quadratic deviation of the manipulated variables with respect to a reference trajectory or its instant change is penalized over a control horizon (Nu ≤ Np). This dynamic optimization problem is subject to the dynamic model of the process, the polytopic constraints of the variables and the operational constraints, and its solution is a control sequence over the prediction horizon. Nevertheless, the (N)MPC implements the feedback control law by applying only the first element of the control sequence and solving repeatedly the OCP while moving the horizon forward. Considering a discrete time representation of the model, the objective function of the (N)MPC is Np

min JNMPC =

min J = wJEO‐NMPC + (1 − w)JNMPC Np

= −∑ w(C DDi − CQ Q i) + (1 − w) i=0 Np

×

∑ αx(xi − xsp)

+

i=0

+ βR (R i − R i − 1)2

∑ αx(xi − xsp)2 + (1 − w) i=0 Nu

×

∑ βQ (Q i − Q i− 1)2 + βR (R i − R i− 1)2 i=1

Nu 2

i=0

(18)

where K is the total number of collocation points, N is the number of finite elements, hj is the length of the finite element j, f(tjk, zdjk, zajk) is the right-hand side of the differential equations presented in Section 2 and τ are the Radau roots of the orthogonal polynomial. There is an additional constraint that represents the continuity of the differential state variables between finite elements eq 20. The algebraic states and the manipulated variables may present discontinuities in the frontier of the finite elements. zdj + 1 =

Np

∑ li(xi , ui) = −∑ li(CDDi − CQ Q i)

∑ βQ (Q i − Q i− 1)

2

(23)

This modification ensures that the moving cost is always strongly convex due to the quadratic terms presented, which is a sufficient condition to guarantee the nominal stability of the closed-loop system. In eq 23 w ∈ [0,1), and it represents the trade-off parameter between the two objectives of the function, the economic term and the tracking term.10,37 The objective function can be rewrite combining the weights of the move suppression coefficients and the weight of the tracking objective in just one. Also, all the weights can be expressed relatively to the economic objective, so only αx′ = (αx(1 − w)/w), β′Q = (βQ(1 − w)/w) and β′R = (βR(1 − w)/w) have to be defined reducing the number of tuning parameters of the controller. Note that all the terms of the objective function are multiply by the same scalar factor, which changes the value of the objective function at the optimal point but does not change the optimal point. This modification is shown in eq 24.

i=1

(21)

In the case of a distillation column control, x is the desired product purity in distillate (state variable), QR and R are the manipulated variables that represent the reboiler duty and the reflux ratio, respectively. The parameters α and β are weights that act as tuning parameters and as linear scalar factors for the multiobjective problem. The set-point is determined from an upper layer, in the hierarchy structure, that considers the economic factors of the process. In other words, it is the results of an economic steady-state optimization. It can be shown that the closed-loop system controlled with this ideal NMPC formulation is asymptotically stable if the control law satisfies certain conditions that make eq 21 a Lyapunov function.12,34 The EO-NMPC has the same characteristic of the formulation and implementation of a NMPC, the only difference is that the objective function of the OCP is a general economic function of the state and manipulated variables instead of a tracking quadratic function. For the distillation column, this function is a trade-off of the energy cost, reboiler duty, and the production revenue, distillate flow rate (eq 22).35

Np

Np

min J = −∑ (C DDi − CQ Q i) + i=0 Np

+

i=0

∑ βQ′ (Q i − Q i− 1)2 + βR′ (R i − R i− 1)2 i=0

6348

∑ αx′(xi − xsp)2

(24)

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research For a better understanding of the EO-NMPC scheme presented here, its stability properties and implementation, the references 10, 36 and 37 could be useful. Two assumptions are made to simplify the control problem: (1) All the states are measurable and they are feedback to the controller and (2) there is no model-plant mismatch. The first assumption avoids the use of an observer such as extended Kalman filter (EKF), Smith estimator, Luenberger observer or moving horizon estimation (MHE),39 it also allows a direct composition control over an indirect control using the most sensitive tray temperature as controlled variable. The second assumption implies that the states predicted are the same as the state measured thus the robustness of the controller is not considered.

Table 1. Operation Condition for the Extractive Distillation Column parameter

value

condenser pressure (bar) equilibrium stages (NT) mixture ethanol−water feed flow rate (mol/min) feed stage of the mixture ethanol−water ethanol molar percentage in the feed stream water molar percentage in the feed stream glycerol feed flow rate (mol/min) glycerol feed stage molar purity of ethanol in distillate

1 19 1666.67 12 80% 20% 583.33 3 ≥99.5%

Table 2. Configuration of the NMPC Controller

5. CASE STUDY The ethanol dehydration is studied because this biofuel has gained importance as an alternative to substitute oil fuels. Also, the bioethanol can be used in mixtures with gasoline in order to increase the octane rating, but the purity of ethanol has to be higher than 99.5% molar to avoid phase separations. During the fermentation process or the biomass reactions used to produce bioethanol, high amounts of water are produced too, forming an azeotropic mixture (∼90% molar of ethanol) of minimum boiling point that make the conventional distillation infeasible.35,40 Alternatives to overcome the azeotropic point and satisfy the purity constraint of the product are vacuum distillation, azeotropic distillation and extractive distillation. There are many variations of extractive distillation according to the entrainer used, and it has been reported that the glycerol can have a good performance for this system.35,41 Furthermore, the glycerol has an economical advantage over other solvents (e.g., benzene) because it is produced as a subproduct of the transesterification reactions to produce biodiesel. The problem studied is the control of a distillation column used in the separation of a mixture ethanol−water using glycerol as entrainer when the system is subject to disturbances in the inlet flow conditions. The initial point of operation of the distillation column is considered to be a point of steady-state operation which is computed from simulation. Additionally, in this case, the initial point of operation is the result of a steadystate operation that minimizes the negative of the process economic profit that correspond to the first term of the EONMPC controller objective function (eq 24). If the system does not start form a steady state that is an economic optimum and if there are no disturbances, using the EO-NMPC controller, the system moves to a new steady state that maximizes the economic profit of the process. The disturbances are applied after the first 30 min of operation and they are only considered in the feed stream of the mixture ethanol−water: (1) +6.25% step in the ethanol composition (ZEtOH ) [D1], (2) −6.25% step in the ethanol composition m (ZEtOH ) [D2] and (3) +10% step in the inlet flow rate (Fm) [D3]. m In this case in particular, ideal behavior of the vapor phase is assumed and the nonideality of the liquid phase is modeled with NRTL (nonrandom two liquids). All the physical properties of the system are obtained from Aspen Properties.20 The operational conditions for the distillation column are the same as those reported by Ramos et al.,35 and they are listed in Table 1. Additionally, the configuration of the EO-NMPC controller is presented in Table 2. It is important to point out that the weights of the control objectives in the NMPC objective function are different from those reported by Ramos et al.35 for the OCP

sampling time (min) prediction horizon length control horizon length periods simulated α′R Weight for ethanol purity β′R Weight for reboiler duty βR′ Weight for reflux ratio reboiler duty cost (CQ) ($/MJ) distillate revenue (CD) ($/mol)

5 10 10 60 10 0.5 0.01 1 × 10−3 30 × 10−3

because the variation of the control effort is penalized instead of the deviation of the manipulated variables respect to a reference point, also this modification allows more rapid changes of the manipulated variables (reflux ration and reboiler duty) that can improve the economic profit of the process represented in the first part of the EO-NMPC objective function. The weights of the EO-NMPC objective function (eq 24) are selected such that the response of the manipulated variables is not too drastic and the control objectives can be achieved. Also, a sensitivity analysis for these parameters is performed and it is found that as βR′ or βQ′ decrease, the economic profit of the process increases, but too drastic of changes in the manipulated variables are observed. Increasing α′R produces a tighter control of the ethanol purity in distillate but induces large and drastic changes in the manipulated variables. The parameters selected represent a compromise good enough between the economic profit of the process and a response of the manipulated variables that is not too drastic and that can be implemented in a real process. The tuning parameters of the EO-NMPC controller used in this study correspond to a first approach to evaluate the performance of the DAE2h model in a NMPC scheme and the advantages of considering the economic optimization within the control loop. In a future work, the tuning parameters of the controller will be improved using an adequate methodology.

6. RESULTS AND DISCUSSION The control problem is solved using the software GAMS 23.9.5 with the CONOPT algorithm that is based on Generalized Reduced Gradient (GRG) in an Intel Core i5 2.70 GHz, 8.0 GB memory. The OCP is solved for one simulation period using the two models presented in Section 2 and the new model proposed in Section 3: index two−direct solution (DAE2), index one obtained after applying an index reduction technique (DAE2r) and index hybrid (DAE2h). A step disturbance in the feed ethanol composition (80% to 75%) is applied after 5 min of 6349

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research

compared between the two models.) is small enough, for example for the variables presented in Figure 2 the AAEs are (a) 0.1489 mol/min, (b) 1.11 × 10−4 and (c) 5.16 × 10−3 MJ/min, these differences between the two models are explained by the precision of the discretization method. Orthogonal collocation with Radau point has an error of O(h2s−1), where s is the number of collocation points, for the algebraic variables of an index one DAE, whereas for the algebraic variables of an index two DAE, the error is O(hs); the precision in the differential variables is the same for index one and index two DAEs.28,33 The index hybrid DAE model proposed is adequate to describe in detail the dynamic behavior of a distillation process because it is a combination of the index two model and its equivalent reduced index one model, both presented in Section 2. These two models are based on fundamental principles; therefore they do not apply drastic simplifications or adjustable parameters. It is called index hybrid DAE model because it uses the index one reduced model to define consistent initial conditions at the first point of each finite element, and then the integration along the interior points of the finite elements uses the index two model without any modification. Table 3 shows the dimensions of the

operation in order to evaluate the changes in the state and manipulated variables. Figure 2 presents the response of one

Table 3. Size of Models after Discretization

number of variables number of equations

index one reduced DAE

index hybrid DAE

index two DAE

101781 101740

55611 55570

40212 40161

three models considered here after discretization using the parameters presented in Table 1 and Table 2, it is important to point out that without consistent initial conditions for the index two model it is not possible to know if the solution is correct. Note that the total number of variables and equations obtained after discretization is significantly lower for the index hybrid DAE model than for the index one reduced DAE model. It is important to clarify that the total number of equations and variables presented in Table 2 are obtained after discretizing and implementing the model in an algebraic programming language, thus it includes intermediate programming calculations. However, the decomposition of the original model equations in different parts to obtain a more structured and organized program is the same for the three models. The way they are programmed is the same but the main advantage of the DAE2h model is that it uses a simpler definition of the energy balances that does not require as many intermediate calculations along the time domain as the energy balances of the DAE2r model. The model proposed (DAE2h) is also used in the EO-NMPC for the disturbance rejection study using the three disturbances and the configuration presented in Section 5, it is also compared with the EO-NMPC using the model DAE2r and with a standard PI controller to illustrate the general benefits of the NMPC. The PI control strategy uses two SISO PI, the first one manipulates the reflux flow rate whereas the second controller manipulates the reboiler duty; both are used to control the ethanol fraction in distillate. The IMC (internal model control) rules are used to identify the tuning parameters of the PI controllers. This tuning strategy uses an open-loop step disturbance of the manipulated variables to identify a first-order plus dead time model of the plant, then the parameters of the PI controller are determined by comparison, and a correction is done using the closed-loop time constant or filter parameter.42 Multiple step tests are done varying the reflux flow rate and the reboiler duty around the

Figure 2. Comparison OCP solutions. (a) Vapor flow rate of the stage 12. (b) Reflux ratio. (c) Reboiler duty.

state variable (Vapor flow rate, (a)) and manipulated variables (reflux ratio (b) and reboiler duty (c)), it is clear that the direct solution of the index two model (DAE2) is not valid because the consistency of the initial conditions is lost at the first point of the finite elements and there is no equation to define this variable at these points. Also, it is not possible to force continuity of the vapor flow rate because it depends directly on the reboiler duty (eq 15), which is a discontinuous function. The drastic changes of the vapor flow rate due to the problem described produce a bad prediction of the dynamic behavior as can be observed in the manipulated variables (Figure 2b,c). In contrast, the models DAE2r and DAE2h do not present the high index problems described and they make a good prediction of the dynamic evolution of the distillation process. Finally, there is no practical difference between the models DAE2r and DAE2h because the analytical solutions of both problems are the same, also the average absolute error (AAE) (The AAE is calculated as AAE = (1/(tf − t0))∫ tft0|yDAE2r − yDAE2h|dt where y is any variable to be 6350

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research nominal operation point of the column to identify a first-order approximate model of the ethanol composition in distillate. Using the gain and time constant of this approximate model the constants of each PI controller are calculated according to the IMC tuning procedure. The parameters obtained are listed in Table 4.

Figure 3 also shows the profiles for the disturbance no. 2 (−6.25% step in the ethanol composition of the feed stream). It is clear that the PI does not satisfy the product purity constraint, the ethanol percentage decreases below a significant level in the time interval of 30 to 120 min (Figure 3a), whereas the control with NMPC implements more drastic changes in the two manipulated variables to satisfy the purity constraint (Figure 3c,d). Also, the NMPC always satisfy this constraint for the fact that it considers a detailed model of the process dynamics and solves a constraint dynamic optimization problem at each sampling time. Moreover, the EO-NMPC has a final distillate production rate similar to the one obtained with PI controller, but it reduces the reboiler duty in order to increase the process profit that is represented in the first term of the OCP objective function (eq 24). The results for the disturbance no. 3 (+10% step in the inlet flow rate) are shown in Figure 4. When a PI control strategy is used the ethanol composition in distillate decreases significantly over the time interval of 30 to 100 min. Similar to the other two disturbances, in this case the NMPC satisfies the purity constraint and decreases the reboiler duty in order to increase the economic profit of the process. Note that despite the increase in the inlet flow rate the distillation production rate does not increase significantly at the final steady state, which indicates that the entrainer flow rate is not sufficient to increase the production; therefore, to improve the objective function (eq 24), the NMPC controller reduces the reboiler duty because it is associated with the energy cost. In other scenarios where it is expected a high variability of the inlet flow rate, the solvent flow rate may be considered as an additional manipulated variable in order to increase the distillate production rate and also the process economic profit. In contrast with disturbances 1 and 2, for this case a small difference between the two models used for NMPC is observed at ∼30 min. It is explained by the fact that the

Table 4. Tuning Parameters for PI Controllers reflux flow rate control reboiler duty control

gain

integral time constant [h]

0.6039 (kmol/h/%) −82.1246 (MJ/h/%)

0.0220 0.0302

The closed-loop profiles for the disturbance no. 1 (+6.25% step in the ethanol composition of the feed stream) are shown in Figure 3. In this case, the purity constraint is always satisfied for all the control strategies proposed due to the nature of the disturbance (Figure 3a). Increasing the ethanol composition of the feed stream does not decrease the product purity in an openloop configuration with constant values of the manipulated variables, thus with any control strategy, the purity constraint is satisfied unless the controller applies a drastic change of the manipulated variables generated by bad tuning parameters. The closed-loop response of both controllers is similar, but the NMPC controller response faster than the PI controller. The behavior of the manipulated variables is different because the PI does not consider the economics of the process (Figure 3c,d), it just focuses on maintaining the product purity at the desired setpoint. On the other hand, the EO-NMPC reaches a different steady state with a higher distillate production rate (Figure 3b), thus it has a better economic profit. The final steady state reached by the NMPC controller improves the economic profit of the process in 0.06% compared to the final steady state of the PI controller.

Figure 3. Closed-loop profiles for disturbance 1 (+6.25% ZEtOH ) [D1] and disturbance 2 (−6.25% ZEtOH ) [D2] using NMPC with two models (DAE2h m m and DAE2r) and PI control. (a) Ethanol composition in distillate. (b) Distillate flow rate. (c) Reflux ratio. (d) Reboiler duty. 6351

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research

Figure 4. Closed-loop profiles for disturbance 3 (+10% Fm) using NMPC with two models (DAE2h and DAE2r) and PI control. (a) Ethanol composition in distillate. (b) Distillate flow rate. (c) Reflux ratio. (d) Reboiler duty.

integration algorithm, orthogonal collocation on finite elements, is more accurate for index one problems than index two problems when the time step is lower than 1. However, this difference does not imply an issue for the control performance because the behavior predicted is good enough. Moreover, when the system is the real distillation column, the model does not match the plant measurements exactly; also, there are unmeasured noises, and these have a bigger impact on the control performance than the deviations observed in the prediction of the models DAE2h and DAE2r. Introducing an observer or predictor can improve the quality of the prediction obtained by the model. Figure 3 and Figure 4 show that the dynamic behavior of the system controlled with NMPC is independent of the model used. The models DAE2h and DAE2r describe the same phenomena under certain mathematical consideration as consistent initial conditions and proper discretization method, but the NMPC using the model DAE2h has a better computational performance because it has a less number of equation than the model DAE2r, thus it requires less time to solve each OCP reducing the computational delay of the feedback action. Furthermore, in some cases the time required to reach an optimal solution of the OCP using the DAE2r model within the EO-NMPC scheme is greater than one sampling time (300s), which introduces significant delays or even instabilities of the closed-loop controller in real implementations. Table 5 shows a comparison of the background task results for the NMPC using both models. Note that NMPC using the model DAE2h requires a shorter average computational time to solve the OCPs than the NMPC using the model DAE2r, and for the NMPC with DAE2h, the maximum time required to solve the OCP never exceed the sampling time established for the system. Table 5 shows the total profit of the process during 5 h operation, this profit is calculated as the integral over all the time of the first term of the EO-NMPC objective function (eq 24)

Table 5. Computational Results, Background Task NMPC average computational time [s]

maximum computational time [s]

model

DAE2h

DAE2r

DAE2h

DAE2r

disturbance 1 disturbance 2 disturbance 3

31.59 42.49 23.55

62.90 60.53 31.64

201.86 220.41 219.30

984.99 238.39 215.76

rejecting the time interval when the distillate product is out of specifications. The EO-NMPC strategy always increases the process profit in comparison with the classic PI controller, even when the product is always in specification (disturbance no. 1) the EO-NMPC increases the profit 0.655% that for long operation times, it represents a significant revenue. Note in Figure 3 and Figure 4, the EO-NMPC reaches a different steady state with a better economic profit from the PI final steady state because it changes the two manipulated variables faster after the disturbance is applied, whereas the PI controller responds slower and only focus on maintaining the ethanol composition in distillate. The tuning parameters of the PI limits the performance of the controller; for example, a faster response that increases the economic profit can be obtained increasing the gain but the integral action that guarantees zero steady-state error is banished. (Table 6).

7. CONCLUSIONS AND PERSPECTIVES In this study, an index hybrid DAE model based on fundamental principles is proposed to describe the dynamic behavior of an extractive distillation column used in the ethanol dehydration process; also, the EO-NMPC controller of the same unit is implemented using this model as a predictor of the process behavior. The proposed model uses an index two DAE to do the numerical integration and its equivalent reduced index one 6352

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research CQ = reboiler duty cost ($/MJ) F = feed streamflow rate (mol/min) HL = liquid specific molar enthalpy (MJ/mol) K = thermodynamic equilibrium constant (−) L = liquid flow rate (mol/min) ML = liquid molar hold-up (mol) NC = total number of components NP = prediction horizon length Nu = control horizon length P = pressure (Pa) QC = condenser duty (MJ/min) QR = reboiler duty (MJ/min) R = molar reflux ratio (−) t = time (min) T = temperature (K) V = vapor flow rate (mol/min) x = liquid molar fraction (−) y = vapor molar fraction (−) z = feed stream molar fraction (−)

Table 6. Economic Performance Control Strategies economic profit after 300 min of operation [$] control strategy

disturbance 1

disturbance 2

disturbance 3

EO-NMPC (DAE2h) EO-NMPC (DAE2r) PI NMPC economic advantage

1.25225 × 10 1.25227 × 104 1.24407 × 104 0.655%

1.11387 × 10 1.11387 × 104 0.78253 × 104 42.34%

1.16868 × 104 1.16873 × 104 0.87700 × 104 33.26%

4

4

model to define consistent initial conditions at the beginning of each finite element, thus it is an index hybrid DAE model. It is shown that the model proposed has the same solution of the reduced index one DAE model and that it decreases the computational time required to solve the optimization problem associated with the EO-NMPC controller. Also, with this model, the online computational time never exceeds the sampling time of the controller whereas the model based control using the index one reduced model exhibits a maximum time greater than the sampling time, which in a real implementation introduces delay and possible instabilities. Three different disturbances are applied over the inlet conditions of the ethanol−water stream and a PI controller is implemented in order to compare the performance of the EONMPC controller. The EO-NMPC improves the net economic profit of the process compare to the PI because it satisfies the product quality constraint during all the operation time. The EONMPC controller allows higher production rates with less energy consumption in comparison with the PI controller. The different disturbances demonstrate that the EO-NMPC controller can handle nonlinear dynamics over a wide range of operation conditions. Future work deals with the improvement of the tuning parameter of the EO-NMPC. This type of controller has too many tuning parameters that are generally defined by trial and error, and many simulations are required to obtain a configuration that satisfies the designer preferences for all the possible cases. Multiobjective optimization or stochastic optimization algorithms could be used to define them in a better way.



Greek Symbols



REFERENCES

(1) Camacho, E. F.; Bordons, C. Model Predictive Control; Springer: London, 1999. (2) Alpbaz, M.; Karacan, S.; Cabbar, Y.; Hapoğlu, H. Application of model predictive control and dynamic analysis to a pilot distillation column and experimental verification. Chem. Eng. J. 2002, 88 (1−3), 163−174. (3) Volk, U.; Kniese, D. W.; Hahn, R.; Haber, R.; Schmitz, U. Optimized multivariable predictive control of an industrial distillation column considering hard and soft constraints. Control Eng. Pract. 2005, 13 (7), 913−927. (4) Ramesh, K.; Shukor, S. R. A.; Aziz, N. Nonlinear model predictive control of a distillation column using NARX model. Comput.-Aided Chem. Eng. 2009, 27, 1575−1580. (5) Venkateswarlu, C.; Reddy, A. D. Nonlinear model predictive control of reactive distillation based on stochastic optimization. Ind. Eng. Chem. Res. 2008, 47 (18), 6949−6960. (6) Purohit, J. L.; Patwardhan, S. C.; Mahajani, S. M. DAE-EKF-based nonlinear predictive control of reactive distillation systems exhibiting input and output multiplicities. Ind. Eng. Chem. Res. 2013, 52 (38), 13699−13716. (7) Marlin, T. E. Process Control: Designing Processes and Control Systems for Dynamic Performance, 2 ed.; McGraw-Hill: Boston, 2000. (8) Ellis, M.; Christofides, P. D. Integrating dynamic economic optimization and model predictive control for optimal operation of nonlinear process systems. Control Eng. Pract. 2014, 22 (0), 242−251. (9) Rawlings, J. B.; Bonne, D.; Jorgensen, J. B.; Venkat, A. N.; Jorgensen, S. B. Unreachable setpoints in model predictive control. IEEE Trans. Autom. Control 2008, 53 (9), 2209−2215. (10) Amrit, R.; Rawlings, J. B.; Biegler, L. T. Optimizing process economics online using model predictive control. Comput. Chem. Eng. 2013, 58 (0), 334−343. (11) Porfírio, C. R.; Odloak, D. Optimizing model predictive control of an industrial distillation column. Control Eng. Pract. 2011, 19 (10), 1137−1146. (12) Biegler, L. T. A survey on sensitivity-based nonlinear model predictive control. In Proceedings of the 10th IFAC International Symposium on Dynamics and Control of Process Systems, Mumbai, India, December 18−20, 2013; International Federation of Automatic Control: Laxenburg, Austria, 2013; pp 499−510. (13) Flatby, P.; Skogestad, S.; Lundstrom, P. Rigorous dynamic simulation of distillation columns based on UV-flash. In Proceedings of

AUTHOR INFORMATION

Corresponding Author

*J. M. Gómez. E-mail: [email protected]. Notes

The authors declare no competing financial interest.



α = weight parameter for the controlled variables β = weight parameter for the manipulated variables γ = activity coefficient (−)

NOMENCLATURE

Subscript

0 = initial condition a = algebraic variable d = differential variable f = final condition i = set of components m = mixture ethanol−water feed to the column n = set of stages for the distillation column (1, condenser; NT, reboiler) Superscript

EtOH = ethanol sat = saturation conditions Latin Symbols

CD = distillate revenue ($/mol) 6353

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354

Article

Industrial & Engineering Chemistry Research the IFAC Symposium on Advanced Control of Chemical Processes, Kyoto, Japan, 1994; International Federation of Automatic Control: Laxenburg, Austria, 1994; pp 261−266. (14) (a) Monteil, C. Techniques de L’ingénieur. Colonnes à plateaux: Dimensionnement, J263. Paris, 1965; Vol. II. (b) Zuiderweg, F. J. Sieve trays: A view on the state of the art. Chem. Eng. Sci. 1982, 37 (10), 1441− 1464. (15) Raghunathan, A. U.; Soledad Diaz, M.; Biegler, L. T. An MPEC formulation for dynamic optimization of distillation operations. Comput. Chem. Eng. 2004, 28 (10), 2037−2052. (16) Huang, R.; Zavala, V. M.; Biegler, L. T. Advanced step nonlinear model predictive control for air separation units. J. Process Control 2009, 19 (4), 678−685. (17) Kumar, A.; Daoutidis, P. Control of nonlinear differential algebraic equation systems: An overview. In Nonlinear Model Based Process Control, Berber, R., Kravaris, C., Eds.; Springer: Dordrecht, The Netherlands, 1998; Vol. 353, pp 311−344. (18) Rahul, M.; Kumar, M. V. P.; Dwivedi, D.; Kaistha, N. An efficient algorithm for rigorous dynamic simulation of reactive distillation columns. Comput. Chem. Eng. 2009, 33 (8), 1336−1343. (19) Sharma, N.; Singh, K. Model predictive control and neural network predictive control of TAME reactive distillation column. Chem. Eng. Process.: Process Intensif. 2012, 59 (0), 9−21. (20) Aspen-Technology. Aspen ONE; Aspen Technology: Cambridge, MA. (21) Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations; SIAM, Society for Industrial and Applied Mathematics: New York, 1996. (22) Cash, J. R. Review paper: Efficient numerical methods for the solution of stiff initial-value problems and differential algebraic equations. Proc. R. Soc. London, Ser. A 2003, 459 (2032), 797−815. (23) Ascher, U. M.; Petzold, L. R. Computer Methods For Ordinary Differential Equations And Differential-Algebraic Equations; SIAM, Society for Industrial and Applied Mathematics: New York, 1998. (24) Estévez Schwarz, D. Consistent initialization for DAEs in Hessenberg form. Numerical Algorithms 2009, 52 (4), 629−648. (25) Barton, P. I. Dynamic Modeling and Simulation Notes; Massachusetts Institute of Technology: Cambridge, MA, 1997. (26) Lozano Santamaría, F.; Gómez, J. M. Index hybrid differential− Algebraic equations model based on fundamental principles for nonlinear model predictive control of a flash separation drum. Ind. Eng. Chem. Res. 2015, 54 (7), 2145−2155. (27) Logsdon, J. S.; Biegler, L. T. Accurate solution of differentialalgebraic optimization problems. Ind. Eng. Chem. Res. 1989, 28 (11), 1628−1639. (28) Biegler, L. T. Nonlinear Programming: Concepts, Algorithms and Applications to Chemical Process; SIAM, Society for Industrial and Applied Mathematics: Philadelphia, 2010. (29) Maciejowski, J. M. Predictive Control with Constraints; Prentice Hall: Upper Saddle River, NJ, 2002. (30) Biegler, L. T.; Cervantes, A. M.; Wächter, A. Advances in simultaneous strategies for dynamic process optimization. Chem. Eng. Sci. 2002, 57 (4), 575−593. (31) Biegler, L. T. An overview of simultaneous strategies for dynamic optimization. Chem. Eng. Process.: Process Intensif. 2007, 46 (11), 1043− 1053. (32) Lima, E. R. A.; Castier, M.; Biscaia, E. C. Differential-algebraic approach to dynamic simulations of flash drums with rigorous evaluation of physical properties. Oil Gas Sci. Technol. 2008, 63 (5), 677−686. (33) Aubry, A.; Chartier, P. On improving the convergence of Radau IIA methods applied to index 2 DAEs. SIAM J. Numer. Anal. 1998, 35 (4), 1347−1367. (34) Grune, L.; Jurgen, P. Nonlinear Model Predictive Control: Theory and Algorithms; Springer: New York, 2011. (35) Ramos, M. A.; García-Herreros, P.; Gómez, J. M.; Reneaume, J.M. Optimal control of the extractive distillation for the production of fuel-grade ethanol. Ind. Eng. Chem. Res. 2013, 52 (25), 8471−8487.

(36) Gopalakrishnan, A.; Biegler, L. T. Economic nonlinear model predictive control for periodic optimal operation of gas pipeline networks. Comput. Chem. Eng. 2013, 52 (0), 90−99. (37) Huang, R.; Biegler, L. T.; Harinath, E. Robust stability of economically oriented infinite horizon NMPC that include cyclic processes. J. Process Control 2012, 22 (1), 51−59. (38) Huang, R.; Harinath, E.; Biegler, L. T. Lyapunov stability of economically oriented NMPC for cyclic processes. J. Process Control 2011, 21 (4), 501−509. (39) Zavala, V.; Biegler, L. Nonlinear programming strategies for state estimation and model predictive control. In Nonlinear Model Predictive Control; Magni, L., Raimondo, D., Allgöwer, F., Eds.; Springer: Berlin/ Heidelberg, 2009; Vol. 384, pp 419−432. (40) Luque, R.; Herrero-Davila, L.; Campelo, J. M.; Clark, J. H.; Hidalgo, J. M.; Luna, D.; Marinas, J. M.; Romero, A. A. Biofuels: A technological perspective. Energy Environ. Sci. 2008, 1 (5), 542−564. (41) García-Herreros, P.; Gòmez, J. M.; Gil, I. n. D.; Rodríguez, G. Optimization of the design and operation of an extractive distillation system for the production of fuel grade ethanol using glycerol as entrainer. Ind. Eng. Chem. Res. 2011, 50 (7), 3977−3985. (42) Romagnoli, J. A.; Palazoglu, A. Introduction to Process Control; CRC Press, Taylor & Francis Group: Boca Raton, FL, 2006.

6354

DOI: 10.1021/acs.iecr.5b00853 Ind. Eng. Chem. Res. 2015, 54, 6344−6354