Nonlinear Autoregressive with Exogenous Inputs Based Model Predictive Control for Batch Citronellyl Laurate Esterification Reactor

1. Introduction

Esterification is a widely employed reaction in organic process industry. Organic esters are most frequently used as plasticizers, solvents, perfumery, as flavor chemicals and also as precursors in pharmaceutical products. One of the important ester is Citronellyl laurate, a versatile component in flavors and fragrances, which are widely used in the food, beverage, cosmetic and pharmaceutical industries. In industry, the most common ester productions are carried out in batch reactors because this type of reactor is quite flexible and can be adapted to accommodate small production volumes (Barbosa-Póvoa, 2007). The mode of operation for a batch esterification reactor is similar to other batch reactor processes where there is no inflow or outflow of reactants or products while the reaction is being carried out. In the batch esterification system, there are various parameters affecting the ester rate of reaction such as different catalysts, solvents, speed of agitation, catalyst loading, temperature, mole ratio, molecular sieve and water activity (Yadav and Lathi, 2005). Control of this reactor is very important in achieving high yields, rates and to reduce side products. Due to its simple structure and easy implementation, 95% of control loops in chemical industries are still using linear controllers such as the conventional Proportional, Integral & Derivative (PID) controllers. However, linear controllers yield satisfactory performance only if the process is operated close to a nominal steady-state or if the process is fairly linear (Liu & Macchietto, 1995). Conversely, batch processes are characterized by limited reaction duration and by non-stationary operating conditions, then nonlinearities may have an important impact on the control problem (Hua et al., 2004). Moreover, the control system must cope with the process variables, as well as facing changing operation conditions, in the presence of unmeasured disturbances. Due to these difficulties, studies of advanced control strategy have received great interests during the past decade. Among the advanced control strategies available, the Model Predictive Control (MPC) has proved to be a good control for batch reactor processes (Foss et al., 1995; Dowd et al., 2001; Costa et al., 2002; Bouhenchir et al., 2006). MPC has influenced process control practices since late 1970s. Eaton and Rawlings (1992) defined MPC as a control scheme in which the control algorithm optimizes the manipulated variable profile over a finite future time horizon in order to maximize an objective function subjected to plant models and constraints. Due to these features, these model based control algorithms can be extended to include multivariable systems and can be formulated to handle process constraints explicitly. Most of the improvements on MPC algorithms are based on the developmental reconstruction of the MPC basic elements which include prediction model, objective function and optimization algorithm. There are several comprehensive technical surveys of theories and future exploration direction of MPC by Henson, 1998, Morari & Lee, 1999, Mayne et al., 2000 and Bequette, 2007. Early development of this kind of control strategy, the Linear Model Predictive Control (LMPC) techniques such as Dynamic Matrix Control (DMC) (Gattu and Zafiriou, 1992) have been successfully implemented on a large number of processes. One limitation to the LMPC methods is that they are based on linear system theory and may not perform well on highly nonlinear system. Because of this, a Nonlinear Model Predictive Control (NMPC) which is an extension of the LMPC is very much needed.

NMPC is conceptually similar to its linear counterpart, except that nonlinear dynamic models are used for process prediction and optimization. Even though NMPC has been successfully implemented in a number of applications (Braun et al., 2002; M’sahli et al., 2002; Ozkan et al., 2006; Nagy et al., 2007; Shafiee et al., 2008; Deshpande et al., 2009), there is no common or standard controller for all processes. In other words, NMPC is a unique controller which is meant only for the particular process under consideration. Among the major issues in NMPC development are firstly, the development of a suitable model that can represent the real process and secondly, the choice of the best optimization technique. Recently a number of modeling techniques have gained prominence. In most systems, linear models such as partial least squares (PLS), Auto Regressive with Exogenous inputs (ARX) and Auto Regressive Moving Average with Exogenous inputs (ARMAX) only perform well over a small region of operations. For these reasons, a lot of attention has been directed at identifying nonlinear models such as neural networks, Volterra, Hammerstein, Wiener and NARX model. Among of these models, the NARX model can be considered as an outstanding choice to represent the batch esterification process since it is easier to check the model parameters using the rank of information matrix, covariance matrices or evaluating the model prediction error using a given final prediction error criterion. The NARX model provides a powerful representation for time series analysis, modeling and prediction due to its strength in accommodating the dynamic, complex and nonlinear nature of real time series applications (Harris & Yu, 2007; Mu et al., 2005). Therefore, in this work, a NARX model has been developed and embedded in the NMPC with suitable and efficient optimization algorithm and thus currently, this model is known as NARX-MPC.

Citronellyl laurate is synthesized from DL-citronellol and Lauric acid using immobilized Candida Rugosa lipase (Serri et. al., 2006). This process has been chosen mainly because it is a very common and important process in the industry but it has yet to embrace the advanced control system such as the MPC in their plant operation. According to Petersson et al. (2005), temperature has a strong influence on the enzymatic esterification process. The temperature should preferably be above the melting points of the substrates and the product, but not too high, as the enzyme’s activity and stability decreases at elevated temperatures. Therefore, temperature control is important in the esterification process in order to achieve maximum ester production. In this work, the reactor’s temperature is controlled by manipulating the flowrate of cooling water into the reactor jacket. The performances of the NARX-MPC were evaluated based on its set-point tracking, set-point change and load change. Furthermore, the robustness of the NARX-MPC is studied by using four tests i.e. increasing heat transfer coefficient, increasing heat of reaction, decreasing inhibition activation energy and a simultaneous change of all the mentioned parameters. Finally, the performance of NARX-MPC is compared with a PID controller that is tuned using internal model control technique (IMC-PID).

2. Batch esterification reactor

The synthesis of Citronellyl laurate involved an exothermic process where Citronellol reacted with Lauric acid to produce Citronellyl Laurate and water.

Figure 1.

Schematic represent esterification of Citronellyl laurate

The esterification process took place in a batch reactor where the immobilized lipase catalyst was mixed freely in the reactor. A layout of the batch esterification reactor with associated heating and cooling configurations is shown in Fig.2.

Figure 2.

Schematic diagram of the batch esterification reactor.

Typical operating conditions were 310K and 1 bar. The reactor temperature was controlled by manipulating the water flowrate within the jacket. The reactor’s temperature should not exceed the maximal temperature of 320K, due to the temperature sensitivity of the catalysts (Yadav & Lathi, 2004; Serri et. al., 2006; Zulkeflee & Aziz, 2007). The reactor’s temperature control can be achieved by treating the limitation of the jacket’s flowrate, Fj, which can be viewed as a state of the process and as the constraint control problem. The control strategy proposed in this paper was designed to meet the specifications of the laboratory scale batch reactor at the Control Laboratory of School of Chemical Engineering, University Sains Malaysia, which has a maximum of 0.2 L/min limitation on the jacket’s flowrate. Therefore, the constraint of the jacket’s flowrate will be denoted as F_jmax= 0.2 L/min.

The fundamental equations of the mass and energy balances of the process are needed to generate data for empirical model identification. The equations are valid for all $t\in \left[0,\infty \right]$ . The reaction rate and kinetics are given by (Yadav & Lathi, 2004; Serri et. al., 2006; Zulkeflee & Aziz, 2007):

$$\frac{d{C}_{Al}}{dt}=\frac{\left[C_{Al}\right]r_{\max }}{\alpha K_{Al}\left(1+\left(\frac{K_{Ac}}{\left[C_{Ac}\right]}\right)+\left(\frac{K_{Ac}}{K_i}\right)+\left(\frac{\left[C_{Ac}\right]}{\beta K_i}\right)\right)}+\left[C_{Al}\right]\left(1+\left(\frac{\alpha K_{Ac}}{\left[C_{Ac}\right]}\right)\right)$$ E1

$$\frac{d{C}_{Ac}}{dt}=\frac{\left[C_{Ac}\right]r_{\max }}{\alpha K_{Ac}\left(1+\left(\frac{K_{Ac}}{\left[C_{Ac}\right]}\right)\right)}+\left[C_{Ac}\right]\left(1+\left(\frac{\alpha K_{Al}}{\left[C_{Al}\right]}\right)\right)+\frac{\alpha K_{Al}{K}_{Ac}}{K_i\left[C_{Al}\right]}+\frac{\alpha K_{Al}\left[C_{Ac}\right]}{\beta K_i\left[C_{Al}\right]}$$ E2

$$\frac{d{C}_{Es}}{dt}=-\frac{\left[C_{Al}\right]r_{\max }}{\alpha K_{Al}\left(1+\left(\frac{K_{Ac}}{\left[C_{Ac}\right]}\right)+\left(\frac{K_{Ac}}{K_i}\right)+\left(\frac{\left[C_{Ac}\right]}{\beta K_i}\right)\right)}+\left[Al\right]\left(1+\left(\frac{\alpha K_{Ac}}{\left[C_{Ac}\right]}\right)\right)$$ E3

$$\frac{d{C}_W}{dt}=-\frac{\left[C_{Ac}\right]r_{\max }}{\alpha K_{Ac}\left(1+\left(\frac{K_{Ac}}{\left[C_{Ac}\right]}\right)\right)}+\left[C_{Ac}\right]\left(1+\left(\frac{\alpha K_{Al}}{\left[C_{Al}\right]}\right)\right)+\frac{\alpha K_{Al}{K}_{Ac}}{K_i\left[C_{Al}\right]}+\frac{\alpha K_{Al}\left[C_{Ac}\right]}{\beta K_i\left[C_{Al}\right]}$$ E4

$$K_i=A_i{\exp }^{-E_i/R{T}_r}$$ E5

$$K_{Ac}=A_{Ac}{\exp }^{-E_{Ac}/R{T}_r}$$ E6

$$K_{Al}=A_{Al}{\exp }^{-E_{Al}/R{T}_r}$$ E7

where $C_{Ac},C_{Al},C_{Es}$ and $C_W$ are concentrations (mol/L) of Lauric acid, Citronellol, Citronellyl laurate and water respectively; r_max (mol l^-1 min^-1 g^-1 of enzyme) is the maximum rate of reaction,K _Ac (mol l^-1 g^-1 of enzyme), K _Al (mol l^-1 g^-1 of enzyme) and K _i (mol l^-1 g^-1 of enzyme) are the Michealis constant for Lauric acid, Citronellol and inhibition respectively; $A_i,A_{Ac}$ and $A_{Al}$ are the pre-exponential factors (L mol/s) for inhibition, Lauric acid and Citronellol respectively; $E_i,E_{Ac}$ and $E_{Al}$ are the activation energy (J mol/K) for inhibition, acid lauric and Citronellol respectively;R is the gas constant (J/mol K).

The reactor can be described by the following thermal balances (Aziz et al., 2000):

$$\frac{d{T}_r}{dt}=\Delta H_{rxn}{r}_{Ac}V+\frac{\dot{Q}}{\left[V\left(C_{Ac}C{p}_{Ac}+C_{Al}C{p}_{Al}+C_{Es}C{p}_{Es}+C_{W}C{p}_W\right)\right]}$$ E8

$$\frac{d{T}_j}{dt}=\frac{(F_{j}C{p}_w{\rho }_w\left(T_{jin}-T_j\right)-\dot{Q})}{V_{j}C{p}_w{\rho }_w}$$ E9

$$\dot{Q}=UA\left(T_j-T_r\right)$$ E10

where T _r (K), T _j (K) and T _jin is reactor, jacket and inlet jacket temperature respectively; $\Delta H_{rxn}$ (kJ/mol)is heat of reaction; V(l) and V _j (l) is the volume of the reactor and jacket respectively; $C{p}_{Ac},C{p}_{Al},C{p}_{Es}$ and $C{p}_W$ are specific heats (J/mol K) of Lauric acid, Citronellol, Citronellyl laurate and water respectively; ${\rho }_w$ is the water density (g/L) in the jacket; $F_j$ is the flowrate of the jacket (L/min); $\dot{Q}$ (kW) is the heat transfer through the jacket wall; A and U are the heat exchange area (m²) and the heat exchange coefficient (W/m²/K) respectively.

Eq. 1 - Eq. 10 were simulated using a 4^th/5^th order of the Runge Kutta method in MATLAB® environment. The model of the batch esterification process was derived under the assumption that the process is perfectly mixed where the concentrations of $[Ac]$ , $[Al]$ , $[Es]$ , $[w]$ and temperature of the fluid in the tank is uniform. Table 1 shows all the value of the parameters for the batch esterification process under consideration.The validations of corresponding dynamic models have been reported in Zulkeflee & Aziz (2007).

Parameters	Units	Values	Parameters	Units	Values
AAc	L mol/s	18.20871	Cpw	J/mol K	75.40
AAl	L mol/s	24.04675	V	L	1.5
Ai	L mol/s	0.319947	Vj	L	0.8
EAc	J mol/K	-105.405	Q	J/m^3	11.648
EAl	J mol/K	-66.093	ΔHrxn	kJ	16.73
Ei	J mol/K	-249.944	α	-	1
Tji	K	294	β	-	1
CpAc	J/mol K	420.53	U	J/s m^2 K	2.857
CpAl	J/mol K	235.27	A	m^2	0.077
CpEs	J/mol K	617.79	R	J/mol K	8.314

Table 1.

Operating Conditions and Calculated Parameters

3. NARX model

The Nonlinear Autoregressive with Exogenous inputs (NARX) model is characterized by the non-linear relations between the past inputs, past outputs and the predicted process output and can be delineated by the high order difference equation, as follows:

$$y\left(t\right)=f\left\{y\left(t-1\right),\dots y\left(t-n_y\right),u\left(t-1\right)\dots u\left(t-n_u\right)\right\}+e(t)$$ E11

where $u\left(t\right)$ and $y(t)$ represents the input and output of the model at time $t$ in which the current output $y(t)\in \Re$ depends entirely on the current input $u\left(t\right)\in \Re$ . Here $n_u$ and $n_y$ are the input and output orders of the dynamical model which are $n_u\ge 0$ , $n_y\ge 1$ . The function $f$ is a nonlinear function. $\overline{X}=[y\left(t-1\right)\dots y\left(t-n_y\right)u\left(t-1\right)\dots u{(t-n_{u)}]}^T$ denotes the system input vector with a known dimension $n=n_y+n_u$ . Since the function $f$ is unknown, it is approximated by the regression model of the form:

$$\begin{array}{l}y\left(t\right)={\displaystyle \sum }_{i=0}^{n_u}a\left(i\right).u(t-i)+{\displaystyle \sum }_{j=1}^{n_y}b\left(j\right).y(t-j)+{\displaystyle \sum }_{i=0}^{n_u}{\displaystyle \sum }_{j=i}^{n_u}a\left(i,j\right).u\left(t-i\right).u\left(t-j\right)+\\ {\displaystyle \sum }_{i=1}^{n_y}{\displaystyle \sum }_{j=i}^{n_y}b\left(i,j\right).y\left(t-i\right).y\left(t-j\right)+{\displaystyle \sum }_{i=0}^{n_u}{\displaystyle \sum }_{j=1}^{n_y}c\left(i,j\right).u\left(t-i\right).y\left(t-j\right)+e(t)\end{array}$$ E12

where $a\left(i\right)$ and $a\left(i,j\right)$ are the coefficients of linear and nonlinear for originating exogenous terms; $b\left(i\right)andb\left(i,j\right)$ are the coefficients of the linear and nonlinear autoregressive terms; $c(i,j)$ are the coefficients of the nonlinear cross terms. Eq. 12 can be written in matrix form:

$$\left[\begin{array}{c}y(t)\\ \begin{array}{c}y\left(t+1\right)\\ ⋮\\ y(t+n_y)\end{array}\end{array}\right]=a\text{.}{u}^T+b\text{.}{y}^T+A\text{.}{[U]}^T+B\text{.}{[Y]}^T+C\text{.}{[X]}^T$$ E13

where

$$a={[a\left(0\right)a\left(1\right)\dots a\left(n_u\right)]}^T$$ E14

$$b={[b\left(1\right)b\left(2\right)\dots b\left(n_y\right)]}^T$$ E15

$$A={[a\left(0,0\right)a\left(0,1\right)\dots a\left(0,n_u\right)a\left(1,1\right)\dots a\left(n_u,n_u\right)]}^T$$ E16

$$B={[b\left(1,1\right)b\left(1,2\right)\dots b\left(1,n_y\right)b\left(2,2\right)\dots b\left(n_y,n_y\right)]}^T$$ E17

$$C={[c\left(0,1\right)c\left(0,2\right)\dots c\left(0,n_y\right)c\left(1,1\right)\dots c\left(n_u,n_y\right)]}^T$$ E18

$$u=[u\left(t\right)u\left(t-1\right)\dots u\left(n_u\right)]$$ E19

$$y=[y\left(t-1\right)u\left(t-2\right)\dots u\left(n_y\right)]$$ E20

$$U=[u\left(t\right)\text{.}u\left(t\right)u\left(t\right)\text{.}u\left(t-1\right)\dots u\left(t\right)\text{.}u\left(t-n_u\right)u\left(t-1\right)\text{.}u\left(t-1\right)\dots u\left(t-n_u\right)\text{.}u\left(t-n_u\right)]$$ E21

$$Y=\left[\begin{array}{l}y\left(t-1\right)\text{.}y\left(t-1\right)y\left(t-1\right)\text{.}y\left(t-2\right)\dots \\ y\left(t-1\right)\text{.}y\left(t-n_y\right)y\left(t-2\right)\text{.}y\left(t-2\right)\dots y\left(t-n_y\right)\text{.}y\left(t-n_y\right)\end{array}\right]$$ E22

$$X=\left[\begin{array}{l}u\left(t\right)\text{.}y\left(t-1\right)u\text{.}y\left(t-2\right)\mathrm{..}u\left(t\right)\text{.}y\left(t-n_y\right)\\ u\left(t-1\right)\text{.}y\left(t-1\right)\mathrm{..}u\left(t-n_u\right)\text{.}u\left(t-n_y\right)\end{array}\right]$$ E23

The Eq. 13 can alternatively be expressed as:

$$y\left(t\right)=[u^T{y}^T{U}^T{Y}^T{X}^T]\left[\begin{array}{c}a\\ b\\ \begin{array}{c}A\\ \begin{array}{c}B\\ C\end{array}\end{array}\end{array}\right]$$ E24

and can be simplified as:

$$\stackrel{}{\overline{Y}}=\stackrel{}{U}\text{.}\stackrel{}{C}$$ E25

where

$$\stackrel{}{\overline{Y}}=y(t)$$ E26

$$\stackrel{}{U}=[u^T{y}^T{U}^T{Y}^T{X}^T]$$ E27

$$\overline{C}={[abABC]}^T$$ E28

Finally, the solution of the above identification problem is represented by

$$\overline{C}=\overline{U}/\overline{Y}$$ E29

The procedures fora NARX model identification is shown in Fig.3. This model identification process includes:

Figure 3.

NARX model identification procedure

• Identification pre-testing: This study is very important in order to choose the important controlled, manipulated and disturbance variables. A preliminary study of the response plots can also gives an idea of the response time and the process gain.

• Selection of input signal: The study of input range has to be done, to calculate the maximal possible values of all the input signals so that both inputs and outputs will be within the desired operating conditions range. The selection of input signal would allow theincorporationof additional objectives and constraints, i.e. minimum or maximum input event separations which are desirable for the input signals and the resulting process behavior.

• Selection of model order: The important step in estimating NARX models is to choose the model order. The model performance was evaluated by the Means Squared Error (MSE) and Sum Squared Error (SSE).

• Model validation: Finally, the model was validated with two sets of validation data which were unseen independent data sets that are not used in NARX model parameter estimation.

The details of the identification of NARX model for the batch esterification can be found at Zulkeflee & Aziz (2008).

4. MPC algorithm

The conceptual structure of MPC is depicted in Fig. 4. The conception of MPC is to obtain the current control action by solving, at each sampling instant, a finite horizon open-loop optimal control problem, using the current state of the plant as the initial state. The desired objective function is minimized within the optimization method and related to an error function based on the differences between the desired and actual output responses. The first optimal input was actually applied to the plant at time t and the remaining optimal inputs were discarded. Meanwhile, at time t+1, a new measurement of optimal control problem was resolved and the receding horizon mechanism provided the controller with the desired feedback mechanism (Morari & Lee, 1999; Qin & Badgwell, 2003; Allgower, Findeisen & Nagy, 2004).

Figure 4.

Basic structure of Model Predictive Control

A formulation of the MPC on-line optimization can be as follows:

$$\underset{u\left[t\text{|}t\right],\dots u{[m+p|t]}_{i=0}}{\min }\underset{u\left[t\text{|}t\right],\dots u{[m+p|t]}_{i=0}}{\min }J(y\left(t\right),u\left(t\right))$$ E30

$$\underset{u\left[t\text{|}t\right],\dots u{[m+p|t]}_{k=0}}{\min }\underset{u\left[t\text{|}t\right],\dots u{[m+p|t]}_{k=0}}{\min }{\displaystyle \sum }_{k=1}^P{w}_k{(y\left[t+k\text{|}t\right]-y_{sp})}^2+{\displaystyle \sum }_{k=1}^M{r}_k\Delta u{[t+k|t]}^2$$ E31

Where P and M is the length of the process output prediction and manipulated process input horizons respectively with P ≤ M. $u{[t+k|t]}_{k=0,\dots P}$ is the set of future process input values. The vector $w_k$ is the weight vector.

The above on-line optimization problem could also include certain constraints. There can be bounds on the input and output variables:

$$u_{\max }\ge u[t+k|t]\ge u_{\min }$$ E32

$$\Delta u_{\max }\ge \Delta u[t+k|t]\ge -\Delta u_{\min }$$ E33

$$y_{\max }\ge y[t+k|t]\ge y_{\min }$$ E34

It is clear that the above problem formulation necessitates the prediction of future outputs

$$y[t+k|t]$$ E35

In this NARX model, for k step ahead:

The error e(t):

$$\begin{array}{l}e\left[t\text{|}t\right]=y\left(t\right)-{\displaystyle \sum }_{i=0}^{n_u}a\left(i\right).u(t-i)-{\displaystyle \sum }_{j=1}^{n_y}b\left(j\right).y(t-j)-{\displaystyle \sum }_{i=0}^{n_u}{\displaystyle \sum }_{j=i}^{n_u}a\left(i,j\right).u\left(t-i\right).u\left(t-j\right)-\\ {\displaystyle \sum }_{i=1}^{n_y}{\displaystyle \sum }_{j=i}^{n_y}b\left(i,j\right).y\left(t-i\right).y\left(t-j\right)-{\displaystyle \sum }_{i=0}^{n_u}{\displaystyle \sum }_{j=1}^{n_y}c\left(i,j\right).u\left(t-i\right).y\left(t-j\right)\end{array}$$ E36

The prediction of future outputs:

$$\begin{array}{l}y\left(t+k\right)={\displaystyle \sum }_{i=0}^{n_u}a\left(i\right).u(t-i+k)+{\displaystyle \sum }_{j=1}^{n_y}b\left(j\right).y(t-j+k)\\ +{\displaystyle \sum }_{i=1}^{n_u}{\displaystyle \sum }_{j=i}^{n_u}a\left(i,j\right).u\left(t-i+k\right).u\left(t-j+k\right)\\ +{\displaystyle \sum }_{i=1}^{n_y}{\displaystyle \sum }_{j=i}^{n_y}b\left(i,j\right).y\left(t-i+k\right).y\left(t-j+k\right)\\ +{\displaystyle \sum }_{i=0}^{n_u}{\displaystyle \sum }_{j=1}^{n_y}c\left(i,j\right).u\left(t-i+k\right).y\left(t-j+k\right)+e\left(t+k\right)\end{array}$$ E37

Substitution of Eq. 35 and Eq. 36 into Eq 31 yields:

$$\begin{array}{l}\underset{u\left[t|t\right],\mathrm{...}u{\left[m+p|t\right]}_{k=0}}{\min }{\displaystyle \sum _{k=1}^P{W}_k}\left(\left({\displaystyle \sum _{i=0}^{n_u}a\left(i\right)\cdot u\left(t-i+k\right)}\right)+{\displaystyle \sum _{j=1}^{n_y}b\left(j\right)\cdot y\left(t-j+k\right)}\right)\\ +{\displaystyle \sum _{i=0}^{n_u}{\displaystyle \sum _{j=i}^{n_u}a\left(i,j\right)\cdot u\left(t-i+k\right)}}\cdot u\left(t-j+k\right)\\ +{\displaystyle \sum _{i=1}^{n_y}{\displaystyle \sum _{j=i}^{n_y}b\left(i,j\right)\cdot y\left(t-i+k\right)}}\cdot y\left(t-j+k\right)\\ +{\displaystyle \sum _{i=0}^{n_u}{\displaystyle \sum _{j=1}^{n_y}c\left(i,j\right)\cdot u\left(t-i+k\right)}}\cdot y\left(t-j+k\right)+y\left(t\right)-{\displaystyle \sum _{i=0}^{n_u}a\left(i\right)}\cdot u\left(t-i\right)\\ -{\displaystyle \sum _{j=1}^{n_y}b\left(j\right)}\cdot y\left(t-j\right)-{\displaystyle \sum _{i=0}^{n_u}{\displaystyle \sum _{j=i}^{n_u}a\left(i,j\right)\cdot u\left(t-i\right)}}\cdot u\left(t-j\right)\\ -{\displaystyle \sum _{i=1}^{n_y}{\displaystyle \sum _{j=i}^{n_y}b\left(i,j\right)\cdot y\left(t-i\right)}}\cdot y\left(t-j\right)-{\displaystyle \sum _{i=0}^{n_u}{\displaystyle \sum _{j=1}^{n_y}c\left(i,j\right)\cdot u\left(t-i\right)}}\cdot y\left(t-j\right)\\ -y_{sp}^2+{\displaystyle \sum _{k=1}^M{r}_k}\Delta u{\left[t+i|t\right]}^2\end{array}$$ E38

Where

$$y_{sp}\left(t\right)={[y_{sp}\left(t+1\right)y_{sp}\left(t+2\right)\dots \text{.}{y}_{sp}\left(t+P\right)]}^T$$ E39

$$\Delta u\left(t\right)={[\Delta u\left[t\text{|}t\right]\Delta u\left[t+1\text{|}t\right]\dots \Delta u\left[t+M-1\text{|}t\right]]}^T$$ E40

The above optimization problem is a nonlinear programming (NLP) which can be solved at each time t. Even though the input trajectory was calculated until M-1 sampling times into the future, only the first computed move was implemented for one sampling interval and the above optimization was repeated at the next sampling time. The structure of the proposed NARX-MPC is shown in Fig. 5.

In this work, the optimization problem was solved using constrained nonlinear optimization programming (fmincon) function in the MATLAB. A lower flowrate limit of 0 L/min and an upper limit of 0.2 L/min and a lower temperature limit of 300K and upper limit of 320K were chosen for the input and output variables respectively. In order to evaluate the performance of NARX-MPC controller, the NARX-MPC has been used to track the temperature set-point at 310K. For the set-point change, a step change from 310K to 315K was introduced to the process at t=25 min. For load change, a disturbance was implemented with a step change (+10%) for the jacket temperature from 294K to 309K. Finally, the performance of NARX-MPC is compared with the performance of PID controller. The parameters of PID controller have been estimated using the internal model based controller. The details of the implementation of IMC-PID controller can be found in Zulkeflee & Aziz (2009).

Figure 5.

The structure of the NARX-MPC

5. Results

5.1. NARX model identification

The input and output data for the identification of a NARX model have been generated from the validated first principle model. The input and output data used for nonlinear identification are shown in Fig. 6.The minimum-maximum range input (0 to 0.2 L/min) under the amplitude constraint was selected in order to achieve the most accurate parameter to determine the ratio of the output parameter. For training data, the inputs signal for jacket flowrate was chosen as multilevel signal. Different orders of NARX models which was a mapping of past inputs (n _u ) and output (n _y ) terms to future outputs were tested and the best one was selected according to the MSE and SSE criterion. Results have been summarized in Table 2. From the results, the MSE and SSE value decreased by increasing the model order until the NARX model with n _u =1 and n _y = 2. Therefore, the NARX model with n _u =1 and n _y = 2 was selected as the optimum model with MSE and SSE equal to 0.0025 and 0.7152 respectively. The respective graphical error of identification for training and validation of estimated NARX model is depicted in Fig. 7.

5.2. NARX-MPC

The identified NARX model of the process has been implemented in the MPC algorithm. Agachi et al., (2007) proposed some criteria to select the significant tuning parameters (prediction horizon, P; control horizon, M; penalty weight matrices w _k and r _k ) for the MPC controller. In many cases, the prediction (P) and control horizons (M) are introduced as P>M>1 due to the fact that it allows consequent control over the variables for the next future cycles. The value of weighting (w _k and r _k ) of the controlled variables must be large enough to minimize the constraint violations in objective function. Tuning parameters and SSE values of the NARX-MPC controller are shown in Table 3. Based on these results, the effect of changing the control horizon, M for M: 2, 3, 4 and 5 indicated that M=2 gave the smallest error of output response with SSE value=424.04. From the influence of prediction horizon, P results, the SSE value was found to decrease by increasing the number of prediction horizon until P=11 with the smallest SSE value = 404.94. SSE values shown in Table 3 demonstrate that adjusting the elements of the w _k and r _k weighting matrix can improve the control performance. The value of w _k = 0.1 and r _k = 1 had resulted in the smallest error with SSE=386.45. Therefore, the best tuning parameters for the NARX-MPC controller were P=11; M=2; w _k = 0.1 and r _k = 1.

Figure 6.

Input output data for NARX model identification

Model	Training		Validation1		Validation2
(nu,ny )	mse	sse	mse	sse	mse	sse
0,1	0.0205	6.1654	0.0285	8.5909	0.0254	7.6357
1,1	0.0202	6.0663	0.0307	9.2556	0.0251	7.5405
2,1	0.0194	5.8419	0.0392	11.8036	0.0266	8.0157
1,2	0.0025	0.7512	0.0034	1.0114	0.0059	1.7780
2,2	0.0026	0.7759	0.0029	0.8639	0.0038	1.1566
3,2	0.0024	0.7289	0.0035	1.0625	0.0097	2.9141
2,3	0.0024	0.7143	0.0033	0.9930	0.0064	1.9212

Table 2.

MSE and SSE values of NARX model for different number of n _u and n _y

Figure 7.

Graphical error of identification for the training and validation of estimated NARX model

Tuning Parameter	SSE	Tuning Parameter	SSE
M=2	424.04	wk = 10	410.13
M=3	511.35	wk = 1	404.94
M =4	505.26	wk = 0.1	386.45
M=5	509.95	wk = 0.01	439.37
with P= 7; wk= 1; rk= 1		with P= 11; M= 2; rk = 1
P =7	424.04	rk = 10	439.23
P =10	405.31	rk = 1	386.45
P =11	404.94	rk = 0.1	407.18
P =12	406.06	rk = 0.01	410.02
with M= 2; wk= 1; rk= 1		with P =11; M=2; wk =0.1

Table 3.

Tuning parameters and SSE criteria for applied controllers in set-point tracking

The responses obtained from the NARX-MPC and the IMC-PID controllers with parameter tuning, K _c =8.3; T _I =10.2; T _D =2.55 (Zulkeflee & Aziz, 2009) during the set-point tracking are shown in Fig. 8. The results show that the NARX-MPC controller had driven the process output to the desired set-point with a fast response time (10 minutes) and no overshoot or oscillatory response with SSE value = 386.45. In comparison, the output response for the unconstrained IMC-PID controller only reached the set-point after 25 minutes and had shown smooth and no overshoot response with SSE value = 402.24. However, in terms of input variable, the output response for the IMC-PID controller has shown large deviations as compared to the NARX-MPC. Normally,actuator saturation is among the most conventional and notable problem in control system designs and the IMC-PID controller did not take this into consideration. Concerning to this matter, an alternative to set a constraint value for the IMC-PID manipulated variable has been developed. As a result, the new IMC-PID control variable with constraint had resulted in higher overshoot with a settling time of about 18 minutes with SSE=457.12.

Figure 8.

Control response of NARX-MPC and IMC-PID controllers for set-point tracking with their respective manipulated variable action.

With respect to the conversion of ester, the implementation of the NARX-MPC controller led to a higher conversion of Citronellyl laurate (95% conversion) as compared to the IMC-PID, with 90% at time=150min (see Fig. 9).Hence, it has been proven that the NARX-MPC is far better thanthe IMC-PID control scheme.

Figure 9.

Profile of ester conversion for NARX-MPC, IMC-PID-Unconstraint and IMC-PIC controllers.

With a view to set-point changing (see Fig. 10), the responses of the NARX-MPC and IMC-PID for set-point change have been varied from 310K to 315K at t=25min. The NARX-MPC was found to drive the output response faster than the IMC-PID controller with settling time, t= 45min and had shown no overshoot response with SSE value = 352.17.On the other hand, the limitation of input constraints for IMC-PID was evidenced in the poor output response with some overshoot and longer settling time, t= 60min (SSE=391.78). These results showed that NARX-MPC response controller had managed to cope with the set-point change better than the IMC-PID controllers.

Fig. 11 shows the NARX-MPC and the IMC-PID responses for 10% load change (jacket temperature) from the nominal value at t=25min. The NARX-MPC was found to drive the output response faster than the IMC-PID controller. As can be seen in the lower axes of Fig 9, the input variable response for the IMC-PID had varied extremely as compared to the input variable from NARX-MPC. From the results, it was concluded that the NARX-MPC controller with SSE=10.80 was able to reject the effect of disturbance better than the IMC-PID with SSE=32.94.

Figure 10.

Control response of NARX-MPC and IMC-PID controllers for set-point changing with their respective manipulated variable action.

The performance of the NARX-MPC and the IMC-PID controllers was also evaluated under a robustness test associated with a model parameter mismatch condition. The tests were;

1. Test 1: A 30% increase for the heat of reaction, from 16.73 KJ to 21.75 KJ. It represented a change in the operating conditions that could be caused by a behavioral phase of the system.

2. Test 2: Reduction of heat transfer coefficient from 2.857 J/s m² K to 2.143 J/s m² K, which was a 25 % decrease. This test simulated a change in heat transfer that could be expected due to the fouling of the heat transfer surfaces.

3. Test 3: A 50% decrease of the inhibition activation energy, from 249.94 J mol/K to 124.97 J mol/K. This test representeda change in the rate of reaction that could be expected due to the deactivation of catalyst.

4. Test 4: Simultaneous changes in heat of reaction, heat transfer coefficient and inhibition activation energy based on previous tests. This test represented the realistic operation of an actual reactive batch reactor process which would involve more than one input variable changes at one time.

Figure 11.

Control response of NARX-MPC and IMC-PID controllers for load change with their respective manipulated variable action.

Fig.12- Fig.15 have shown the comparison of both IMC-PID and NARX-MPC control scheme’s response for the reactor temperature and their respective manipulated variable action for robustness test 1 to test 4 severally. As can be seen in Fig. 12- Fig. 15, in all tests, the time required for the IMC-PID controllers to track the set-point is greater compared to the NARX-MPC controller. Nevertheless, NARX-MPC still shows good profile of manipulated variable, maintaining its good performance. The SSE values for the entire robustness test are summarized in Table 4. These SSE values shows that both controllers manage to compensate with the robustness. However, the error values indicated that the NARX-MPC still gives better performance compared to the both IMC-PID controllers.

Controller	Test 1	Test 2	Test 3	Test 4
NARX-MPC	415.89	405.37	457.21	481.72
IMC-PID	546.64	521.47	547.13	593.46

Table 4.

SSE value of NARX-MPC and IMC-PID for robustness test

Figure 12.

Control response of NARX-MPC and IMC-PID controllers for robustness Test 1 with their respective manipulated variable action.

Figure 13.

Control response of NARX-MPC and IMC-PID controllers for robustness Test 2 with their respective manipulated variable action.

Figure 14.

Control response of NARX-MPC and IMC-PID controllers for robustness Test 3 with their respective manipulated variable action.

Figure 15.

Control response of NARX-MPC and IMC-PID controllers for robustness Test 4 with their respective manipulated variable action.

6. Conclusion

In this work, the NARX-MPC controller for the Batch Citronellyl Laurate Esterification Reactor has been developed. The validated first principle model was used as a process model to generate data required for NARX model identification. The NARX model with n _u =1 and n _y = 2 was chosen since it gave the best performance with MSE and SSE equal to 0.0025 and 0.7152 respectively. Finally, the performances of the NARX-MPC controller were evaluated for set-point tracking, set-point change, load change and robustness test. For all cases, the developed controller strategy (NARX-MPC) has been proven to perform well in controlling the temperature of the batch esterification reactor, as compared to the IMC-PID controllers.

Acknowledgments

The authors wish to acknowledge the financial support by Universiti Sains Malaysia (USM), for the research funding through the Research University Grant scholarship and Ministry of Science, Technology and Innovation, Malaysia (MOSTI) for the scholarship for the first author.