Index Hybrid DifferentialâAlgebraic Equations Model Based on

Article pubs.acs.org/IECR

Index Hybrid Differential−Algebraic Equations Model Based on Fundamental Principles for Nonlinear Model Predictive Control of a Flash Separation Drum Federico Lozano Santamaría and Jorge M. Gómez* Grupo de Diseño de Productos y Procesos, Departamento de Ingeniería Química, Universidad de los Andes, Carrera 1 No. 18A-10, Bogotá, Colombia ABSTRACT: Model predictive control (MPC) has been widely used in the chemical process industry; however, the industrial applications use a simplified model to describe the time evolution of the system in order to make the dynamic optimization problem easier to solve. Models based on fundamental principles describe in detail the dynamics of the system, but they define a problem of differential−algebraic equations (DAE) that can present high-index difficulties related to the definition of consistent initial conditions. Index reduction techniques have been applied to solve this kind of problem. This study presents an alternative to using high-index DAE models in MPC. The problem is discretized using orthogonal collocation on finite elements, and the original high-index model is used in the numerical integration, but its equivalent reduced index one model is used to define efficiently consistent initial conditions. The set point tracking problem and the disturbance rejection problem for a flash separation drum using nonlinear model predictive control are considered. It is shown that the proposed model does not present a practical difference with respect to the index one reduced model and that it is more efficient in computational terms.

1. INTRODUCTION Model based predictive control (MPC) is an online feedback control strategy that solves an optimal control problem (OCP) at each instant of time. It has been widely used in the process industry because the main advantages of this strategy are prediction of the behavior of the process in a future period of time, inherent compensation to dead-time and time-delay phenomena, and the ability to handle constraint and the multivariable problem directly.1 Nevertheless, these implementations usually use simplified linear models of the process (nonlinear) to predict its future behavior taking into account that these kinds of models are valid only in the vicinity of the linearization point; beyond this region the prediction is not suitable to match the real process evolution. The advantage of working with linear models is that the global optimum of the OCP is guaranteed. Nonlinear models can be included in the MPC strategy (nonlinear model predictive control, NMPC) to predict the future behavior of the process, but it increases the difficulty of solving the OCP because the optimization problem could be not convex and the global optimum cannot be guaranteed; it also requires specialized nonlinear programming (NLP) algorithms. Thus, there is a compromise between the quality of the prediction and the difficulty of solving the OCP. Some empirical nonlinear models, which do not require an extensive knowledge of the system, have been used in the NMPC of different processes. For example, Hammerstein and Wiener models were used for the level control of a three tank interconnected system,2 artificial neural networks for the pH control of a chemical reactor,3 and wave models for the product composition control of an internal thermally coupled distillation column.4 However, because of their empirical characteristics, the quality of the prediction obtained with these models is also restricted to a certain region, but it is better © XXXX American Chemical Society

than the linear model’s prediction. The quality of the prediction is especially important in the process start-up and shut-down and when it is subject to drastic transitions or disturbances that move the system far from its nominal operation point. The models based on fundamental principles may predict better the behavior of the process over all the operation conditions because they describe how the system changes over time using conservation laws, but they usually define large-scale models with differential−algebraic equations (DAE), increasing the complexity of the dynamic optimization. The dynamic operation of a flash separation drum and a distillation column are represented by DAE models, in which the differential equations correspond to material and energy balances and the algebraic equations represent hydraulic and equilibrium relations. This kind of model cannot be solved using traditional solvers for ordinary differential equations (ODEs) as they require a special treatment for the algebraic equations. Throughout this study the semiexplicit (Hessenberg form) representation of a DAE system is considered; even though this representation does not cover all the DAE systems, it is suitable for most chemical process applications in which the algebraic equations arise naturally from physical constraints. The semiexplicit DAE system consists of a group of differential equations subject to a set of algebraic equations that represent equality constraints (eqs 1a and 1b).5

dy − f (y , z ) = 0 dt

(1a)

Received: September 15, 2014 Revised: January 30, 2015 Accepted: February 3, 2015

A

DOI: 10.1021/ie5036594 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research Table 1. Comparison of Different Studies about the Dynamics of Distillation Columns and Flash Separation year

problem

model used to describe the process

Flatby et al.12

author

1994

Han Kim13

1999

Kumar and Daoutidis8

1999

Raghunathan et al.9

2004

Index one and index two DAE based on fundamental principles. The index two DAE system is reduced to an index one DAE system using differentiation and substitution. Index one DAE system based on fundamental principles. It is obtained after applying an index reduction technique over an index two DAE system. Index one and index two DAE system based on fundamental principles. The index two problem is addressed using state space realization, which is equivalent to a direct differentiation index reduction technique. Index one DAE system based on fundamental principles obtain after applying an index reduction technique of differentiation and substitution over an index two DAE system.

Kawathekar and Riggs11 Gonçalves et al.14 Venkateswarlu and Reddy10 Jogwar and Daoutidis15 Sharma and Singh16

2007

general study about the dynamics of distillation columns optimal design and operation of a batch distillation column nonlinear control of a reactive distillation column using a state space representation dynamic optimization of distillation column with mathematical program with equilibrium constraints (MPEC) nonlinear model predictive control of a reactive distillation column dynamic simulation of a flash separation drum using equations of state nonlinear model predictive control of a reactive distillation column control of a vapor recompression distillation unit model predictive control of a reactive distillation column

Ramos et al.17

2013

g (y , z ) = 0

2007 2008 2009 2012

optimal control of an extractive distillation column

Index one DAE system based on fundamental principles obtain after applying an index reduction technique of differentiation and substitution over an index two DAE system. Index one DAE system based on fundamental principles that consider the vapor hold-up and the pressure variation with respect to time. Polynomial model autoregressive mean average (ARMA) adjusted from a detailed index one DAE model based on fundamental principles. Stiff differential equations representing the mass and energy balances. They use singular perturbations to formulate reduced order models valid over different time scales. Second-order transfer function model and artificial neural networks, both adjusted from an index one DAE model based on fundamental principles obtain after applying an index reduction technique of differentiation and substitution over an index two DAE system. Detailed index one DAE model based on fundamental principles.

algebraic variables (z0) determines easily consistent initial conditions for all the variables.5 On the other hand, the initial conditions of higher-index DAEs have to satisfy the algebraic constraint and other hidden constraints that are associated with the derivatives of the algebraic constraints. For instance, in an index two DAE system at least one algebraic variable is not defined by the subset of algebraic equations g(y0,z0) = 0 when the initial values of the differential variables are fixed; in this case, the subsystem of algebraic equations is singular for the algebraic variables. Consistent initial conditions for index two DAE systems must satisfy g(y0, z0) = 0 and (dg/dt)(y0, z0) = 0; this new set of algebraic equations is nonsingular for the algebraic initial conditions (z0). Higher-index systems include higher-order derivatives of the algebraic constraints in the hidden constraints.7 The derivatives of the constraints are not directly available with the DAE representation of eq 1, which makes consistent initial conditions difficult to define; thus, it compromises the correct numerical integration of the original high-index DAE system. Several authors have studied the dynamic behavior and control of distillation columns and flash separation drums using DAE models, as is shown in Table 1. In the case of offline studies such as optimal control or operational analysis under disturbances, in which the computational time required to solve the problem is not important, large-scale DAE systems based on fundamental principles have been used.8,9 On the other hand, in online applications like model predictive control, simplified models are used because of the high influence of the computational time in the controller performance; however, these models are adjusted from a more detailed model that is generally a DAE system based on fundamental principles.10,11 Notice in Table 1 that when a high-index DAE system arises in the modeling of the process, an index reduction technique is applied to obtain an equivalent index one DAE and, as far as we know, an index two or higher DAE has never been used directly in the dynamic study of distillation columns or flash separation drums. This study uses the direct numerical integration of an index two DAE system used as a prediction model of the future behavior of a flash separation drum in a NMPC controller. The

(1b)

The index of the model is generally an indicator of the difficulty of solving the DAE system. There is not a unique definition of index, and it is usual to define different types of it, for example, differential index, perturbation index, and tractability index, each one with a different mathematical interpretation.6 Generally, the differential index is the most studied and reported for large problems, and it is defined as the minimum number of times that all or part of the original DAE system must be differentiated with respect to t in order to determine (dy/dt) and (dz/dt) as continuous functions of y, z, and t, in other words, a system of ordinary differential equations.5 Among the different alternatives to solve high-index problems, the index reduction techniques are the most common, and they can be classified in three categories: (i) direct differentiation, which consists of differentiating the entire system of equations until it becomes an explicit ODE; (ii) differentiation and substitution, in which some of the algebraic constraints are differentiated to replace a differential equation by its algebraic equivalent; and (iii) forming an augmented DAE system, which is similar to direct differentiation but without differentiating all the equations of the system; for each additional equation a new variable, equivalent to a Lagrange’s multiplier, is added. In the DAE problems related with chemical processes, technique (ii) is the most applied.5,7 However, index reduction is a costly process in terms of model dimension because it introduces more equations to the original model and it requires the analytical calculation of many derivatives; thus, it is advantageous to work directly with the high-index DAE. To solve any index DAE system correctly, the initial conditions of all the variables must be consistent. For index one DAE systems this means that they have to satisfy the algebraic constraints at the initial point g(y0,z0) = 0. It is important to point out that for index one DAE systems the subsystem of algebraic equations is nonsingular for the algebraic variables, so setting arbitrary initial conditions on the differential variables (y0) and solving the subsystem of algebraic equations g(y0,z0) = 0 for the initial value of the B


Article

Industrial & Engineering Chemistry Research

vapor phase are faster than those of the liquid phase; there is a time-scale multiplicity. The vapor phase responds rapidly to perturbations, almost instantaneously; thus, taking this limit condition over the fast variables (vapor hold-ups) eliminates its dynamic response, and they have to be defined by the set of algebraic equations (e.g., equilibrium relations). However, because the pressure is fixed there are fewer algebraic states than algebraic equations, producing a high-index problem19,20 It is valid to apply these assumptions when the operation pressure is less than 10 bar, when the vapor hold-up is less than the 20% of the liquid hold-up, and when there is a very tight control of the pressure, as in distillation columns.20 For this study, ideal vapor phase is assumed and the nonrandom two liquids (NRTL) model is used to describe the behavior of the liquid phase. When these simplifications are applied, the model can be described by the following algebraic and differential equations. Total mass balance:

article is structured as follows: section 2 describes the index two DAE model for a flash separation drum and the index reduction technique applied to obtain the equivalent index one model. Section 3 presents the hybrid index model proposed to determinate efficiently consistent initial conditions and the use of the index two DAE system directly in the numerical integration. Section 4 presents the discretization used to transform the differential equations into a set of algebraic equations in order to solve the OCP as a NLP problem. Section 5 describes the NMPC strategy, assumptions, and configuration. Section 6 illustrates the application of the methodology proposed with a case study corresponding to the control of a flash drum used in the azeotropic separation of a water−ethanol mixture. Section 7 presents the results and discussion of the NMPC controller for the flash separation drum using the methodology proposed. The final section provides conclusions and future challenges.

2. DYNAMIC MODELS FOR THE FLASH SEPARATION DRUM The dynamic model for a flash separation drum consists of material and energy balances between the inlet and outlet process streams (Figure 1) that define differential equations,

dML = F − L − V, dt

M(t0) = M 0

(2)

Partial mass balance: F(zi − xi) − V (yi − xi) dxi = , ML dt xi(t0) = xi ,0 ∀ i ∈ {1, 2, ..., NC}

(3)

Energy balance: d(MLHL) = FHF − LHL − VHV + Q dt

(4)

Thermodynamic equilibrium: yi = K ixi =

Pisatγi P

xi

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

(5)

Phase equilibrium error: Figure 1. Flash separation representation.

NC

∑ (yi − xi) = 0

the thermodynamic equilibrium relations, and physical constraint (e.g., vessel capacity) that defines the algebraic equations. An index one DAE model that considers the dynamics of the vapor phase and all the physical constraints of the system is one of the closest representations of the real process, but this introduces additional algebraic equations and variables. Also, this DAE model is a stiff system because the vapor phase dynamics are faster than the liquid phase dynamics;18 as a result, an adequate numerical method with a step size control algorithm is required. These drawbacks affect the performance of the NMPC controller because they increase the complexity of the associated optimization problem, and computational delays could be observed. Finally, it is usual to apply some model simplifications to make the optimization problem easier to solve, but they increase the index of the DAE system. The next two subsections present the index two DAE system used to model this separation process under certain assumptions and the index reduction technique applied in order to obtain an equivalent index one DAE system. 2.1. Index Two DAE System (DAE2). To simplify the model, it is usual to assume that the pressure is constant and the vapor hold-up is negligible.12 Normally, the dynamics of the

i=1

(6)

Liquid flow rate: L = ρ(T , xi)Ao 2gh

(7)

Equations 1−7 lead to a DAE system, with differential variables ML, xi, and MLHL where i ∈ {1, 2, ..., NC}. This model is an index two DAE system, and it can be seen easily because the vapor flow rate, an algebraic variable, does not appear in any of the algebraic equations. The change of this variable over time cannot be known directly when a method to solve index one problems is applied because the subsystem of algebraic equations is singular. Also, the model has only one degree of freedom or manipulated variable in the case of NMPC (e.g., heat duty). Consistent initial conditions cannot be known directly for this problem because of the singularity of the subsystem of algebraic equations, which limits the application of this model. If arbitrary initial conditions are set for these algebraic variables, which cannot be determined, the solution of the index two DAE system will present impulse behavior and will not be a smooth function; thus, it will not be appropriate for predictive purposes in a NMPC scheme.21 C


Article

Industrial & Engineering Chemistry Research 2.2. Index One Reduced DAE System (DAE2r). An index reduction technique is applied over the index two DAE system (eqs 2−7) to obtain an equivalent representation in an index one DAE system. In the case of distillation columns and flash separation drums, the differential energy balance is replaced by an algebraic energy balance using differentiation and substitution while all the other equations are kept intact.22 First, the left-hand side of the differential energy balance (eq 4) is differentiated to obtain the expression ML

dHL dML + HL = FHF − LHL − VHV + Q dt dt

hand, the subsystem of algebraic equations of the index one reduced DAE system is nonsingular; thus, after the initial conditions for the differential variables (ML and xi) are defined, the square system of algebraic equations can be solved for the initial values of all the algebraic variables, guaranteeing consistent initial conditions and a smooth solution.

3. INDEX HYBRID DAE MODEL The model presented here is a combination of the index two DAE system (DAE2) and its equivalent reduced index one model (DAE2r). The index hybrid DAE system (DAE2h) uses the algebraic equations of the DAE2r system to define consistent initial conditions when the differential states (ML and xi) are fixed. Then, the numerical integration is done using all the equations of the original DAE2 system. The subsystem of algebraic equations of DAE2r defines two factors that allow a direct integration of the DAE2 system: (1) the definition of consistent initial conditions solving the square nonlinear algebraic equations when the differential states are fixed and (2) the definition of an invariant set for the algebraic variables of the DAE2 system.7,23 The first point guarantees a smooth solution without the presence of impulse behavior. In addtion, the initial point for the numerical integration obtained with this methodology belongs to the invariant set defined by the algebraic equations and hidden constraints of the DAE2 system.24 The second point ensures that the hidden constraints of the DAE2 system are satisfied along the numerical integration, even though they are not used directly, because the initial point belongs to the invariant set.7 Using the methodology proposed (DAE2h) there is no need to calculate the derivatives of the index reduction over the entire time interval. They are calculated just at the first point of integration, which reduces significantly the number of equations and variables after discretizing and implementing the differential and algebraic equations in an algebraic programming language, in comparison with the DAE2r system.25 The DAE2h system has a computational advantage over the DAE2r because it has fewer equations after discretization, which implies less computational time for the solving of the system. It is especially important in NMPC controllers where usually the computational time required for the solution of the dynamic optimization problem is a limitation of the controller performance. The most important part of this methodology is the direct numerical integration of the model DAE2 that can be done under consistent initial conditions, but using the equivalent DAE2r system is not the only alternative to define these conditions. For instance, a steady-state solution of the DAE2 system also guarantees consistent initial conditions imposing additional algebraic constraints in the differential equations; in this case, the analytical calculation of the derivatives involved in the index reduction is not necessary. However, not all the numerical methods used to integrate index one DAEs and ODEs are appropriate for high-index DAE systems because it is required that the numerical method have strong stability properties and control of the local integral error in order to avoid a high global error in the algebraic variables.26 It has been reported that collocation methods have the desired properties and that for index two DAEs the minimum requirement to obtain a solution with an acceptable global error is to use threepoint collocation.26 The collocation technique used to solve the high-index DAE system is explained in the next section.

(8)

Then eq 2 is replaced in eq 8: dHL F(HF − HL) − V (HV − HL) + Q = dt ML

(9)

∑i NC = 1xiHL,i,

Consider the liquid enthalpy as HL = where HL,i is the specific liquid molar enthalpy of each component and it is an algebraic function of temperature. This expression is differentiated with respect to time applying the chain rule N dHL dT ⎡⎢ C dHL, i ⎤⎥ = + ∑ xi dt dt ⎢⎣ i = 1 dT ⎥⎦

NC

∑ i=1

dxi HL, i dt

(10)

The expression (dHLi/dT) can be derived analytically knowing the algebraic function of HL,i, and (dxi/dt) is defined by eq 3; thus, the only definition left to construct is an algebraic energy balance (dT/dt), which is obtained by substituting eq 5 in eq 6 and differentiating the resulting equation (eq 11). N

C d {[∑ (K ixi)] − 1} = dt i = 1

NC

⎡⎛ dxi ⎞ ⎛ dK i ⎞⎤ ⎟ + ⎜xi ⎟⎥ ⎣ d t ⎠ ⎝ dt ⎠⎦

∑ ⎢⎜⎝K i i=1

⎛ 1 ⎡ dγ(x , T ) ⎛ dx ⎞ dP sat dT ⎤⎞ i j ⎥⎟ = ∑ ⎜K i i ⎟ + ∑ ⎜⎜xi ⎢Pisat + γi i ⎝ dt ⎠ i = 1 ⎝ P ⎣⎢ dt dT dt ⎥⎦⎟⎠ i=1 NC

NC

(11)

=0

After the chain rule is applied several times and the terms are rearranged, the final expression for the derivative of temperature with respect to time is obtained. dT = dt

dx 1 N ⎡ N ∂γ (xj , T ) dxj ⎤ N − P ∑i =C1 ⎢xiPisat ∑ j =C1 i ∂x dt ⎥ − ∑i =C1 K i dti ⎦ ⎣ j

(

1 P

∂γ (xj , T ) NC x P sat i ∂T i=1 i i

(∑

N

+ ∑i =C1 xiγi

dPisat dT

)

) (12)

where ∂γi(xj,T)/∂xj, ∂γi(xj,T)/∂T, and dPisat/dT can be calculated analytically knowing their respective functions of temperature and composition and (dxi/dt) is defined by eq 3. Finally, replacing eq 12 in eq 10 and the resulting expression in eq 9, the new energy balance is an algebraic equation. The index one reduced DAE system has the same number of equations as its equivalent index two DAE, but this index reduction process requires the analytical calculation of a significant number of derivatives. When the DAE model is transcribed to an algebraic programming language, these derivatives introduce dummy variables and equations for intermediate calculations, increasing the model dimensions. It is advantageous to use directly the index two model because it does not need additional calculations, but with this model, it is not easy to define consistent initial conditions that ensure the validity of the high-index DAE system integration. On the other D


Article


4. DISCRETIZATION OF THE DIFFERENTIAL EQUATIONS All the state and control variables present in the differential algebraic equations are discretized with the purpose of solving the problem using simultaneous optimization strategies that are compatible with NLP algorithms. The discretization used is orthogonal collocation on finite elements; it approximates the state and control variables into a family of orthogonal polynomials inside each finite element. A fixed step length of 0.5 min is used for each finite element, and inside each finite element the variables are discretized around three collocation Radau points. This type of discretization is chosen because (1) it is widely applied to solve index two problems;26,27 (2) it is a numerical method with L-stability, which means that when the integration step tends to zero the stability region tends to infinity;26 (3) it is capable of stabilizing high-index DAE systems;26 and (4) this approach allows a direct transcription of the optimization problem in order to be solved simultaneously in an equation oriented environment that can exploit the dispersity and structural properties of the algebraic system generated.28 Properties 2 and 3 hold only if the initial conditions are consistent, otherwise the numerical method can fail at the first steps of integration showing impulse behavior in the response of the state variables.7,29 Moreover, this type of collocation presents the lowest truncation error for stiff ODE problems after orthogonal collocation with Legendre points, but because of the characteristic of Radau collocation, the truncation error is still relative low for DAE problems (Table 2). However, the method is less accurate in the algebraic variables for index two problems.28,30

at each Radau root is evaluated, the term dli(τk)/dτ in eq 13 becomes the kth element of the collocation matrix. There is an additional constraint that represents the continuity of the differential variables between finite elements (eq 16): the function that describes the evolution of the differential variables has to be differentiable over the entire time interval. The manipulated variables are restricted to piece-wise step functions, and some algebraic variables that are directly related to them may present discontinuities between finite elements. The different behavior of the differential variable, algebraic variable, and manipulated variables is illustrated in Figure 2. K

yj + 1 =

∑ yjk lk(1),

j = 1, 2, ...J − 1 (16)

k=0

Table 2. Truncation Error Inside Each Finite Element for the Method of Orthogonal Colocation with Radau Point DAE index

error in algebraic variablesa

Error in differential variablesa

0 1 2

− O(h2s−1) O(hs)

O(h2s−1) O(h2s−1) O(h2s−1)

a

Figure 2. Orthogonal collocation on finite elements representation.29

s is the number of stages (collocation points) of the method.

Figure 2 also shows where each model presented in section 2 is used to construct the hybrid DAE model proposed. The DAE2r system is used at the first point of the integration to define consistent initial conditions when the differential states and manipulated variables are fixed, while the DAE2 system is used in the numerical integration along the internal points of each finite element. However, the discontinuities of the manipulated variables force the use of the DAE2r system periodically to define consistent initial conditions for each time interval.

When this method is applied, the differential equations become a weighting of the variables using the derivative of an orthogonal polynomial (Lagrange’s polynomial) (eq 13). These polynomials are represented in eq 14. Like any other method for integration of differential equations, it is necessary to define initial conditions for the differential variables at the first collocation point of the first finite element (eq 15). K

∑ yji i=0

dli(τk) = hjf (t jk , yjk ), dτ

k = 1, 2, ..., K , j = 1, 2, ...J (13)

K

lk(τ ) =

∏ n = 0, ≠ k

y1,0 = y(t0)

τ − τn , τk − τn

5. NMPC PROBLEM FORMULATION k = 0, 1, 2, ..., K

In a NMPC scheme at each instant of time (tk), an OCP has to be solved to obtain the control signal to be applied at the current time. Assuming that the current states of the process can be measured or estimated, the OCP is stated as a general dynamic optimization problem described by eq 17.28

(14) (15)

where K is the total number of collocation points, J the number of finite elements, hj the length of the finite element j, and f(tjk, yjk) the right-hand side of the differential equations presented in section 2; τn(k) are the Radau roots of the orthogonal polynomial. After the derivative of the orthogonal polynomial

min J = ϕ(t f ) + E

∫t

tf 0

φ(y(t ), z(t ), u(t )) dt

(17a)


Article

Industrial & Engineering Chemistry Research s.t.

dy = f (y(t ), z(t ), u(t )) dt

(17b)

y(0) = y0

(17c)

h(y(t ), z(t ), u(t )) = 0

(17d)

g (y(t ), z(t ), u(t )) ≤ 0

(17e)

6. CASE STUDY The case study selected corresponds to the dehydration of ethanol. Nowadays bioethanol is the most utilized biofuel as an alternative to oil fuels. It can be produced by two different routes: fermenting sugars and starches (first-generation bioethanol) or from biomass such as lignocellulose (secondgeneration bioethanol).33 Nevertheless, it is necessary that the ethanol has a purity higher than 99.5% molar to avoid phase separation when it is mixed with gasoline.34 Fuels that are blends of ethanol and gasoline contain between 51 and 83 volume percentage of ethanol, and even when the water content is low (less than 0.5 volume percentage) the mixture will separate in two phases, one hydrocarbon phase and the other water phase rich in ethanol; this phase separation may produce several operational problems in vehicle engines.35 Also, fuel grade ethanol is hygroscopic and can absorb water from the atmosphere during storage and transportation.34 For these reasons the water content in the fuel grade ethanol has to be minimized. The bioethanol produced has high amounts of water forming an azeotropic mixture of minimum boiling point (∼90% molar of ethanol at 1 bar) that makes the conventional distillation infeasible.17 There are many alternatives to overcome the azeotropic point, such as vacuum distillation, azeotropic distillation, and extractive distillation. Different possibilities of extractive distillation can be considered according to the entrainer used, and it has been reported that glycerol has a good performance for this system.17,36 Furthermore, glycerol has an economical advantage over other solvents (e.g., benzene) because it is a subproduct of the transesterification reactions in biodiesel production. This process ensures the desired product purity, but it also requires an advanced control strategy because the distillation process is highly nonlinear and there are several multivariable interactions that have to be considered. In addition, an advanced control strategy is required because the process is subject to constant disturbances, the most significant of which is the variability of the raw materials (e.g., biomass quantity processed according to seasons or harvest efficiency, water content of biomass, variations of the market prices)37 that may compromise the product quality and production rates. The study of a flash separation drum is a first approximation to the dynamic behavior of the distillation process because the model that describes the dynamics of an extractive distillation column is an extrapolation of the models presented in section 2. However, the NMPC controller for the extractive distillation column using the index hybrid DAE model to predict the process future behavior will be considered in future studies. The problem studied here consists of the flash separation of the mixture ethanol−water−glycerol. Two examples are analyzed to show the performance of the index hybrid DAE model in a NMPC scheme: 1. Process start-up. This is a set point tracking problem in which the system starts with zero production rates and it is desired that it reaches a specific production rate and purity. 2. Disturbance rejection. The system starts from the steady state reached in example 1, and a sinusoidal perturbation with 10% amplitude and a 1 h period is applied in the feed flow rate. The operation conditions for the flash separation drum are listed in Table 3.

where t0 is the initial time, tf the final time, y a vector of differential state variables, z a vector of algebraic state variables, y0 a vector of initial conditions of the differential variables, and u a vector of manipulated variables. The objective function is composed of two terms: the first one (ϕ(tf)) is a terminal cost that is a function of only 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 hybrid index DAE model described in section 3 contains 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 represents bounds for the states and manipulated variables (polytopic constraints) or specific operational constraints like product purity limits;31 however, they are not considered in this study. The objective function of this problem is formulated as the minimization of the quadratic error of the controlled variables with respect to a set point over a prediction horizon, while the drastic changes of the manipulated variable are penalized over a control horizon. It is important to point out that the control horizon length is always lower than or equal to the prediction horizon length. The first two terms of the objective function (eq 18) represent the set point tracking of the controlled variables that in this case are the vapor flow rate and the composition of the vapor stream, and the third term represents the penalization of the variation of the manipulated variable (heat duty). Np

min Jk =

∑ αV(Vk ,i − Vsp)2 + αy(ykEtOH − yspEtOH )2 ,i i=0

Nu

+

∑ β(Q k ,i − Q k ,i− 1)2 i=0

(18)

In eq 18, the subindex i represents the future prediction of the variable while the subindex k represents the current time instant. The ideal NMPC scheme is considered here. In this methodology, the OCPk is solved at each instant of time (tk) and its solution provides a control sequence over a future time horizon (uk,i), but only the first element of this sequence (uk,0) is feedback to the controller in order to define the control law, and the system moves to the next time instant (tk+1) where a new OPCk+1 is solved. A time step of 0.5 min equal to the finite element length is used in the NMPC scheme to facilitate the state feedback and the time evolution of the system. Two assumptions are made to simplify the control problem: (1) All the states are measurable and they are feedback to the controller. (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),32 whereas the second assumption implies that the states predicted are the same as the states measured. F


Article

Industrial & Engineering Chemistry Research Table 3. Operation Conditions for the Flash Separation Drum parameter operation pressure (bar) feed flow rate (mol/min) feed temperature (K) feed molar percentage water ethanol glycerol tank diameter (m) liquid initial level (m) tank height (m) set point vapor flow rate (mol/min) ethanol composition in vapor control parameters sampling time (min) prediction horizon length control horizon length αV αy β

value 1 1000 357.56 37% 37% 26% 1.6 1.6 3.2 209.5 75% 0.5 30 15 0.01 0.012 0.004

7. 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 of example 1 is solved for one simulation period with the three models: index two direct solution (DAE2), index two reduced to index one (DAE2r), and index hybrid (DAE2h). The response of the vapor flow rate presented in Figure 3a shows impulse behavior; therefore, the direct solution of the index two model is not valid. There is no equation to define the vapor flow rate at the initial point of each finite element, and the algorithm fixed this value at its set point in order to minimize the objective function; thus it is not possible to guarantee the consistency of the initial conditions for each time interval of the numerical integration. This impulse behavior of the response of the controlled variable affects directly the future prediction of the manipulated variable, as is shown in Figure 3b, which makes this model inappropriate for prediction purposes in a NMPC scheme. The impulse behavior and the high model−plant mismatch may introduce closed loop instabilities and significant steady-state errors. The models DAE2r and DAE2h are better options to predict the process behavior in a NMPC scheme. There is no practical difference between the solution of the DAE2r and DAE2h; the mean absolute error in the vapor flow rate for this case is 0.0031 mol/min, and for the heat duty it is 0.0137 kJ/min (Figure 3). 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, while 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,30 Table 4 shows the size of the models after discretization using orthogonal collocation on finite elements with the parameters listed in Table 3 to solve the optimal control

Figure 3. Comparison OCP solutions: (a) vapor flow rate and (b) heat duty.

Table 4. Size of the Models after Discretization number of equations number of variables

DAE2r

DAE2h

DAE2

16 215 16 242

8755 8783

6374 6402

problem at each instant of time in the NMPC scheme. The hybrid model has fewer equations than the index one reduced DAE system because the dummy variables and equations used in the index reduction process (e.g., eq 12) are included in the mathematical programming transcription. This factor affects the computational load and therefore the computational time required to solve the OCP. Note that the DAE2 system has fewer equations and variables as it does not apply any index reduction technique, but its solution is not suitable for prediction in a NMPC scheme. It is important to point out that the number of equations and variables presented in Table 4 includes intermediate programming calculations, that is, the decomposition of the original model equations in different parts to define the original equations in a more structured and organized way. These partitions are the same for both models, DAE2h and DAE2r; thus, the way they are programmed is the same, and the advantage of the DAE2h model is that it uses a simpler definition of the energy balance that does not require as many partitions as the algebraic energy balance used in the DAE2r model. The rest of this section shows the results of the NMPC strategy using the models DAE2r and DAE2h as predictors for the examples described in section 6. G


Article

Industrial & Engineering Chemistry Research 7.1. Example 1. Figure 4 shows the simulation results for the set point tracking problem using a NMPC controller with

does not present high-index integration difficulties (e.g., impulse behavior, failure of the numerical integration at the first steps) because it defines consistent initial condition, and it is computationally easier to solve. Regarding this last point, Figure 5 shows the computational time required to solve the

Figure 5. Computational results example 1: set point tracking, background task NMPC.

OCP at each time step. For the DAE2r, the mean time is 2.7 s with a maximum time of 5.3 s, whereas for the DAE2h the mean time is 1.2 s with a maximum time of 2.5 s. Note that the maximum time is reached close to the first time step because the initialization of the variables is far from its solution, then each new OPC is initialized with the results of the previous one that do not differ significantly. 7.2. Example 2. The presence of disturbances is common in the ethanol dehydration process. They drive the system away from its nominal operational point, and the control strategy should be capable of minimizing their influence. The NMPC controller allows a direct handling of the process disturbances; however, to exploit all the advantages of the predictive controller, it is necessary to know the future behavior of the disturbances so the controller can respond faster to their influence. It is not possible to know a priori the future behavior of the disturbances, but it is possible to use a forecasting model to estimate it.38 Usually the disturbances are represented as time-series models that use past data to forecast its future values,38 but here we consider only two types of forecasting: (1) a constant forecast (CF) of the disturbance equal to the measured value at the current time, that is, the most common forecasting strategy applied,38 and (2) a perfect forecast (PF) of the disturbance that assumes a perfect knowledge of the future behavior of the disturbance. Other types of forecasting models (e.g., autocorrelated models) will be studied in future work. Figure 6 shows the response of the controlled variables and the manipulated variables using the models DAE2h and DAE2r to predict the behavior of the process when a sinusoidal disturbance of the feed flow rate is applied on the flash separation drum. Also, the performance of the NMPC controller is compared under the two forecasting strategies for the disturbance. It is important to note that the manipulated variable (heat duty) responds in a manner consistent with the perturbation applied in order to minimize the variation of the controlled variables with respect to its set points. The process reaches a cyclic steady state because of the characteristic of the perturbation, as shown in Figure 6a. In the case of the perfect forecast, the process reaches this cyclic steady state

Figure 4. Example 1 set point tracking, NMPC control of a flash separation drum: (a) vapor flow rate, (b) ethanol composition in vapor, and (c) heat duty.

the DAE2r and DAE2h as predictor models. The system reaches the new set point in ∼30 min without steady-state error in the two controlled variables (vapor flow rate and ethanol mole fraction). Also, the behavior of the manipulated variable (heat duty) does not change drastically because of the penalization in the objective function, and it can be implemented in the real process without operational problems. As can be observed, there is no practical difference between the responses of the two models, and the mean absolute error in each variable presented in Figure 4 is (a) 1.62 × 10−4 mol/min, (b) 5.77 × 10−5%, and (c) 0.008 kJ/min. The principal advantages of the methodology proposed (DAE2h) are it conserves the same physical description of the process that is provided by the reduced index one model, it H


Article


Figure 7. Computational results example 2: disturbance rejection, background task NMPC.

solution of the OCP using the model DAE2r requires a computational time greater than that using the DAE2h model, regardless of the forecasting strategy used to estimate the future behavior of the disturbance.

8. CONCLUSIONS AND PERSPECTIVES We present a model based on fundamental principles defined by differential algebraic equations to describe the dynamic behavior of a flash separation drum. This model uses an index two DAE system in the numerical integration and its equivalent reduced index one DAE system to define efficiently consistent initial conditions at the first point of each time interval; thus, it is a hybrid index DAE model. It is shown that the solution of the index two DAE under consistent initial conditions its practically equal to the solution obtained using the reduced index one equivalent DAE system. There is a small difference between the two solutions produced by the precision of the numerical method used in the integration (orthogonal collocation). The model proposed is compared to the reduced index one model to show that there is no practical difference between them. This comparison is done analyzing the control of a flash separation drum using a NMPC controller. The set point tracking and disturbance rejection problems are considered separately, and both cases show that the two models (DAE2h and DAE2r) describe the same phenomena, but the solution of the dynamic optimization problem using the hybrid model requires less computational time. The difference in computational time is related with the computational load of each model; the DAE2r model has a higher number of variables and equations because of the introduction of dummy variables associated with the index reduction process. Finally, the methodology proposed may be improved if consistent initial conditions for the index two DAE are defined in a way that does not require the calculation of the derivatives of the algebraic constraints. In terms of modeling, future work deals with the extrapolation of the hybrid index model to a distillation column case, and in terms of NMPC, the effects of disturbance, robust control, state estimation, and NLP sensitivity has to be analyzed. Also, other types of forecasting models, like ARMA or artificial neural networks, may be considered in the NMPC scheme to have a more realistic description of the disturbance behavior over the prediction horizon.

Figure 6. Example 2 disturbance rejection, NMPC control of a flash separation drum: (a) vapor flow rate, (b) ethanol composition in vapor, and (c) heat duty.

instantaneously, whereas the constant forecast case introduces a delay response of the manipulated variable that causes a different behavior of the vapor flow rate during the first period of the disturbance (

Index Hybrid DifferentialâAlgebraic Equations Model Based on

Recommend Documents

Index Hybrid DifferentialâAlgebraic Equations Model Based on