CHEMICAL ENGINEERING TRANSACTIONS  
 

VOL. 56, 2017 

A publication of 

 
The Italian Association 

of Chemical Engineering 
Online at www.aidic.it/cet 

Guest Editors: Jiří Jaromír Klemeš, Peng Yen Liew, Wai Shin Ho, Jeng Shiun Lim 
Copyright © 2017, AIDIC Servizi S.r.l., 

ISBN 978-88-95608-47-1; ISSN 2283-9216 

Auto-Regressive with Exogenous Input Model Predictive 
Controller for Water Activity in Esterification 

Siti Asyura Zulkeflee, Suhairi Abdul Sata, Norashid Aziz* 
School of Chemical Engineering, Engineering Campus, Universiti Sains Malaysia, Seri Ampangan, 14300 Nibong Tebal, 
Seberang Perai Selatan, Malaysia  
chnaziz@usm.my 

In this work, the Auto-Regressive with Exogenous Input Model Predictive Controller (ARX-MPC) was designed 
and implemented to control the water activity of the lipase-catalysed esterification process. The empirical 
model, which was embedded in the MPC was developed using the Autoregressive with Exogenous input 
(ARX) model. The parameter estimation and model validation for the ARX model were carried out using the 
recursive least squares estimation (RLSE) system identification toolbox in MATLAB®. The capability of the 
models to capture the dynamics of the input and output variables was also verified. ARX models were solved 
using the quadratic programming (QP) method in the MPC toolbox in MATLAB®/Simulink. The ARX-MPC 
parameters were tuned to determine the best controller performance. Then, the performances of the best-
tuned ARX-MPC were evaluated in terms of set point tracking and disturbance rejection. According to results, 
the ARX-MPC was considered suitable and reliable for tracking the set-point changes of the controlled 
process variable and able to eliminate the presence of disturbance in the process. 

1. Introduction 
The rising demand for the natural products which are environmentally friendly leads the application of lipases 
as catalyst in many chemical processes. The rate and equilibrium yield as well as the stability of lipases can 
be influenced by many factors such as temperature, enzyme concentration, water activity, biocatalyst loading, 
initial molar substrate ratio and different solvents. Among these factors, the most critical is water activity 
(Wehtje et al., 1997). Most lipase catalysed need a certain amount of water activity in order to work. 
Continuous control of water activity for esterification process is a must. Control of water activity in the lipase-
catalysed esterification process is very important to achieve high yields and rates, as well as to reduce side 
products. A small number of studies has been done in controlling the esterification process, and most of them 
implemented conventional controllers, such as Proportional-Integral (PI) controller (Won and Lee, 2002). This 
type of controller is expected to perform efficiently only when the process operates in a genuinely linear 
behaviour. The control system has to cope with the process uncertainties, as well as the changes in 
operations conditions, in the presence of disturbances. Due to these difficulties, studies of advanced control 
strategies during batch processes have received considerable interest in the past decades (Hosen et al., 
2011). Previous studies have proven that the Model Predictive Control (MPC) is a good control for batch 
reactor processes because its algorithms can cope with multivariable systems and can be formulated to 
handle process constraints explicitly (Santos et al., 2013).  
The MPC can be defined as a control system in which the control algorithm optimises the manipulated 
variable profile over a finite future time horizon in order to maximise an objective function subjected to plant 
models and constraints (Eaton and Rawlings, 1992). The MPC has been successfully implemented in a 
number of chemical processes. The implementation of the MPC controller in an esterification system is still 
scarce. Among the major issues in the MPC development are the developments of a suitable internal model 
that can represent the real process (Mohd and Aziz, 2015). A proper internal model with an efficient 
optimisation algorithm must be embedded in the MPC in order to achieve better performance in controlling the 
esterification process. ARX modelling was regularly used as time series models. Just like any other empirical 
modelling approaches pertaining to the given time series data, the key problem is to determine the most 

                               
 
 

 

 
   

                                                  
DOI: 10.3303/CET1756037

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Please cite this article as: Zulkeflee A.A., Sata S.A., Aziz N., 2017, Auto-regressive with exogenous input model predictive controller for water 
activity in esterification, Chemical Engineering Transactions, 56, 217-222  DOI:10.3303/CET1756037   

217



suitable model complexity that can also depict the dynamic structure of the objective system. The complexity 
of the ARX modelling is determined by the number of model orders which is the number of the past 
consecutive measurable input and output response under a constant time interval. 
In this work ARX model has been developed and used as internal model that can capture the dynamics 
relating the inputs to the output of the esterification process. For control purposes, the developed ARX model 
is embedded in the MPC and known as ARX-MPC. Finally, ARX-MPC performance is evaluated for set-point 
tracking, set-point change and disturbance rejection. 

2. Water activity control system in batch esterification reactor 
In this work, Citronellyl laurate ester was synthesised from Citronellol and lauric acid using immobilised 
Candida Rugosa lipase as a biocatalyst. The esterification process was carried out in a jacketed batch reactor 
with the presence of iso-octane as a solvent. In the water activity control loop system, water activity was 
measured in the gaseous phase with a humidity sensor and connected to a computer for monitoring. The 
MATLAB® software was used to design the program for the controller to control the inflow of dry air, which 
was needed to ensure the water activity was retained at the set point value (MATLAB, 2016). The measured 
value in the reactor was compared with the set point and the air flow was controlled via a proportional valve, 

and directed into a desiccating column to dry the air, respectively. The air was bubbled into the reactor 
through small holes in a narrow steel tube. Details on the process can be found in Zulkeflee et al. (2013). 

3. Auto-regressive with exogenous input model predictive control (ARX-MPC) 
The control objective for the lipase-catalysed esterification process is to control the water activity (aw) by 
manipulating the air flowrate (Fair). The idea of the MPC is to attain the current control action by providing an 
optimal solution in every sampling time through a finite horizon open-loop optimal control problem based on 
the plant’s current state as the primary state. Based on the difference between the required and actual output 
response, there is a reduction of the required objective function within the optimisation method that is linked to 
an error function. The earliest optimal input is attached to the plant at time t, while other optimal inputs that 
come later will be rejected. In the meantime, another size of the optimal control problem is solved at time t + 1, 
while the regressing horizon mechanism offers the preferred feedback mechanism for the controller (Nagy et 
al., 2007). A formulation of the MPC online optimisation used in this work is as follows: 

min
u[t|t],…u[m+p|t]i=0

J(y(t), u(t)) (1) 

min
u[t|t],…u[m+p|t]k=0

∑ wk(y[t + k|t] − ysp)
2 + ∑ rk∆u[t + k|t]

2

M

k=1

P

k=1

 (2) 

where u(t) and y(t) are the input and output variables. P and M are the length of the process output prediction 
and the manipulated process input horizons. u[t + k|t]k=0,…P is the set of future process input values. The 
vector wk and rk represent the weight vector for output and input control response. The notation shows the 
value of the variable at instant t + k calculated at instant t.  
At each time step t, the current state y(t) of the process and the reference input over the finite horizon were 
provided to the controller which then computed the optimal input for the process response. 

ysp(t) = [ysp(t + 1)ysp(t + 2) … . ysp(t + P)]
T (3) 

In this work, the ARX model was developed used as the internal model embedded in the MPC. This structure 
is known as ARX-MPC. Figure 1 depicts the structure of the proposed ARX-MPC. 
ARX-MPC controller was formulated by considering the objective function, the constraints imposed by the 
empirical model, and input and output variables. The ARX model error, e(t) is: 

e[t|t] = y(t) − ∑ a(i). u(t −

nu

i=0

i) − ∑ b(j). y(t −

ny

j=1

j) (4) 

The prediction of future outputs is: 

y(t + k) = ∑ a(i). u(t −

nu

i=0

i + k) + ∑ b(j). y(t −

ny

j=1

j + k) + e(t) (5) 

Substitution of Eq(4) and Eq(5) into Eq(2) yields Eq(6): 

218



min
u[t|t],…u[m+p|t]k=0

∑ wk (∑ a(i). u(t −

nu

i=0

i + k + ∑ b(j). y(t −

ny

j=1

j + k + y(t) − ∑ a(i). u(t −

nu

i=0

i)

P

k=1

− ∑ b(j). y(t −

ny

j=1

j) − ysp)

2

+ ∑ rk∆u[t + k|t]
2

M

k=1

 

(6) 

 

Figure 1: Structure of the ARX-MPC 

The issue of online optimisation as shown above can also incorporate a few limitations. In this work, the 
normal operating condition for the air valve was at 50 %. The valve opening for the jacket and air was set to 
be in the range between 15 % and 85 %. The minimum and maximum constraints of input changes [∆umin, 
∆umax] can be set as [−35 %, 35 %]. The process constraints for this work can be seen in the following 
equations: 

16.67 L/s ≤  Fair[t + k|t]  ≤  83.30  L/s 

0 ≤  𝑎𝑤[𝑡 + 𝑘|𝑡]  ≤  0.8 

(7) 

(8) 

The tuning of the MPC parameter is based on the rule of thumb suggested by Bemporad et al. (2010) for 
application in the MPC toolbox in the MATLAB software. The sampling period, 𝑇𝑠, is approximated to one fifth 
of the process time constant, the small control horizon is selected with the default; M ≥ 1 3⁄ (Tsett Ts⁄ ) where 
Tsett is the settling time and the default prediction horizon is initially based on, P ≥ M. The weighting factor (wk 
and rk) is set to default value, 0.1. The value of the controller parameters is increased until further increases 
have a minor impact on the controller performance. The controller performance of the process response in 
following and tracking the set-point is evaluated based on the error value, rise-time, overshoot, and settling 
time. In this work, ISE is used to evaluate the performance of the developed controller. The ISE criterion 
formula is written as follows: 

ISE = ∫ |ysp(t) − y(t)|
2
dt

∞

0

 (9) 

4. Identification of ARX model 
As a linear recurrence equation, the ARX model is used to join the current value of an objective variable 
together with its past finite time series, as well as other exogenous input variable’s finite time series. This can 

be written as follows:  

y(t)ARX  = ∑ a(i). u(t −

nu

i=0

i) + ∑ b(j). y(t −

ny

j=1

j) + e(t) (10) 

where a(i) is the coefficient of the originating exogenous terms; b(i) is the coefficient of the autoregressive 
terms and e(t) is the process noise. Here, nu and ny are the input and output orders of the dynamical model in 
which nu  ≥  0, ny  ≥  1. The matrix form of the ARX model in Eq(10) can be written as follows:  

[
 
 
 

y(t)

y(t + 1)
⋮

y(t + ny)]
 
 
 

ARX 

= a. uT + b. yT + e(t) (11) 

where: 
a = [a(0) a(1) … a(nu)]

T (12) 
b = [b(1) b(2) … b(ny)]

T (13) 
u = [u(t) u(t − 1) … u(nu)] (14) 

 
OPTIMISER 

PLANT ysp (t) 
y(t) 

z-1 
ARX model 

z-1 

y(t-1) 
constraints 

ai , bi, 
f(t-1) 

y(t-1) 

u(t-1) 

u(t-1) 

219



y = [y(t − 1) u(t − 2) … u(ny)] (15) 
Eq(12) can be alternatively expressed as: 
y(t)ARX = [u

TyT] [
a
b
] (16) 

The constants 𝑎 and 𝑏 are the coefficients of the term cluster which contains the form unu(t − i)yny(t − j) for 
nu + ny ≤ nl, where 𝑖 and 𝑗 are any time lag. Such coefficients are called cluster coefficients and are 
represented as θnuny. Eq(18) can be represented in a simpler form:  
Y̅ = U̅. θnuny  (17) 
where: 
Y̅ = y(t) 
U̅ =  [uTyT] 
θnuny = [a b]

T 

(18) 
(19) 
(20) 

The estimation of the model constant in Eq(19) can be solved by:  
θnuny =  U̅/ Y̅ (21) 
With full state measurements, choosing the number of input and output lags is crucial and the process 
conditions can be initially described by the current states and inputs. The state and input transition matrices 
are estimated by regressing a matrix of concatenated states and inputs on a matrix of forward shifted states. 
In this work, the input of pseudo random binary signal (PRBS) has been implemented in the process 
simulation system and the output signal has been verified to be used in the ARX model identification. The 
minimum and maximum of air flowrate is equal to 16.67 L/s and 83.30 L/s.  
In ARX model identification, the model matrices function is estimated using the recursive least squares 
estimation (RLSE) function block in the MATLAB®/Simulink. The Kalman filter is chosen as the ARX model 
estimation method since it allows a prediction of the output response and then updates the output based on 
the available measurement. The parameters obtained from the model is exported in the MATLAB® workspace 
and used in the model validation procedure. For validation, the multiple step input signal is introduced into the 
process system to generate the validation data set. In developing the ARX models, the important phase that 
should be taken into account is the choice of the model order. (Eq)22 represents the calculation of the Mean 
Square Error (MSE) between the model output response and the real process output response.  

MSE ⟹ ∀u ∈ 𝕌: ξ {|y(t, u, y) − ym (t, u, y, nu, ny, θnuny)|
2
} (22) 

The model ym(t, u, y, nu, ny) converges in mean square sense to a system y(t, u, y), if all the column vectors of 
the model residuals, ε > 0, ∃Mε, is independent of θnuny . 

5. Results and discussion 
Different orders of the ARX models, which involved mapping the past input (nu) and output (ny) terms to future 
outputs, were tested and the best were selected according to the MSE criterion. The results for the variations 
in the model order and the corresponding MSE values for the ARX model are tabulated in Table 1. 

Table 1: MSE values for different numbers of 𝑛𝑢 and 𝑛𝑦 

(nu, ny)  MSE (x 10
-2) 

 Training Validation 
0,1 3.0849 7.9382 
1,1 2.6833 5.5213 
2,1 2.2910 5.0042 
2,2 2.0628 4.9017 
2,3 2.1735 4.9834 
 
For the first identification, all the model orders were set to nu = 0 and ny = 1 to keep the models as simple as 
possible. As can be seen from the results, the ARX model orders, nu = 2 and ny = 2 was selected as the best 
models with MSE = 4.9017 for validation. It can be observed that the MSE value decreased by increasing the 
model order. The MSE values started to increase after the model reached a certain model order, which 
showed a larger value of error. Based on this result, it is important to consider the input terms in the model 
since the variations of the input determine the direction of the output.  
Then, the best developed ARX model was used as the predictive models in the MPC. Table 2 shows the 
compilation results for the tuning parameters and the ISE criteria for the developed controllers in the set point 
tracking. Based on the results, the value of M = 4 gave the smallest error output response with ISE values = 

220



195.43. The results indicated that smaller numbers of the M would lead to a better control performance if 
compared to a long control horizon. This was because the small value of the M had the tendency to increase 
the robustness of the system and reduce computational burden, whereas a long control horizon would result in 
unnecessary control action and a long settling time. The results also show that increasing the P would 
decrease the value of ISE. Further increase of the P led to higher ISE values. The best performances for the 
ARX-MPC was when P = 10. The longer P was required to determine the number of the output predictions 
that were used in the optimisation calculation and caused the control system to be less sensitive to model 
error. As can be seen, the controller performed the best at rk = 0.1. For the  wk, controllers gave the lowest 
values of error at wk= 1. The weighting factor penalised the tracking errors and guided the servo performance 
of the control system. Next, ARX-MPC with the best tuning was evaluated based on their performances on set 
point tracking, set point change and disturbance rejection.  

Table 2: Tuning parameters and ISE criteria for set point tracking 

Tuning M ISE Tuning P ISE Tuning wk ISE Tuning rk ISE 
M = 2 
M = 3 
M = 4 
M = 5 
 
with  
P = 5; wk = 1; 
rk = 1 

385.04 
200.87 
195.43 
231.57 

P = 7 
P = 8 
P = 9 
P = 10  
P = 11 
with 
M = 4; wk = 
1; rk = 1 

185.02 
181.58 
176.92 
174.75 
190.11 
 

wk = 10 
wk = 1 
wk = 0.1 
 
 
with 
P = 10; M = 
4; rk = 1 

318.02 
174.75 
190.04 

rk = 10 
rk = 1 
rk = 0.1 
rk = 0.01 
 
with 
P = 10; M = 
4; wk = 1 

581.01 
174.75 
138.17 
152.22 

 

  

Figure 2: (a) Water activity control responses of the ARX-MPC controllers for set-point tracking and (b) their 

respective manipulated variable actions 

  

Figure 3: (a) Water activity control responses of the ARX-MPC controllers for set-point changes and (b) their 

respective manipulated variable actions 

  

Figure 4: (a) Water activity control responses of the ARX-MPC controllers for disturbance rejection and (b) 

their respective manipulated variable actions. 

, 

, , 

221



The water activity was controlled at its set point value of 0.3. Based on Figure 2(a), the ARX-MPC presented a 
large output overshoot with the ISE value of 46.27. Based on Figure 2(a), the ARX-MPC presented a large 
output overshoot (ISE = 46.27) as the water activity increased rapidly at the beginning of the reaction caused 
from the water production from the esterification process. With some overshoot and fluctuation, the results 
showed that the ARX-MPC controller can drive the water activity to the desired set points. As shown in Figure 
3, the ARX-MPC controllers were able to bring the water activity to the new set point. The ARX-MPC showed 
small fluctuating responses around the set point trajectory. In this test, the ISE value of the the ARX-MPC was 
41.49. A step input with 2 % increase of the jacket temperature variable is introduced after the measured 
process response becomes steady which at allocated time the cooling line was closed and the heating system 
was allowed to run for a few minutes until the jacket temperature has reached a 2 % increase from the initial 
value. In Figure 4(a), it is clear that when the process was disturbed, water activity responses by the ARX-
MPC deviated from the original set-point. A similar behaviour response can be observed from the manipulated 
profiles in Figure 4(b). For all tests, it can be observed that the manipulated variable profiles for ARX-MPC 
controllers fluctuated within the input constraints. In this test, the ISE value for the ARX-MPC was 31.79. As 
observed, the ARX-MPC showed a good performance in rejecting the process disturbance 

6. Conclusions 
The ARX-MPC controller had been proposed in this work. The identification of the ARX models for lipase-
catalysed esterification process was presented. The ARX model with model order nu = 2 and ny = 2 was 
selected as the best predictive models for the water activity control loop system. The identified model served 
as the basis for the design and evaluation of the MPC controllers. The important MPC were tuned to obtain 
the best controller performance. The ARX-MPC with controller parameters: P = 10, M = 4, wk = 1, and rk = 0.1 
was selected since these controllers showed the best performances in controlling water activity. Finally, all the 
controllers were evaluated based on their performances on set point tracking, set point change, and 
disturbance rejection. The performance results of the ARX-MPC controllers revealed that the closed-loop 
responses were satisfactory and effective in tracking the optimal set point. Aim of the future work is to embed 
the nonlinear empirical model in the MPC which can provide better prediction on the reaction process and thus 
will be able to bring the process variable to the set point quickly with smaller amount of overshoot. 

Acknowledgments  

The financial support from University Saints Malaysia through Research University Grant No. 
1001/PJKIMIA/814237 and NSF scholarship from Ministry of Science, Technology & Innovation to the first 
author are greatly acknowledged. 

References  

Bemporad A., Morari M., Ricker N.L., 2010, Model predictive control toolbox, user’s guide, The Mathworks 
Inc., Massachusetts, USA. 

Eaton J.W., Rawlings J.B., 1992, Model-predictive control of chemical processes, Chemical Engineering 
Science 47 (4), 705-720. 

Hosen M.A., Hussain M.A., Mjalli F.S., 2011, Control of polystyrene batch reactors using neural network 
based model predictive control (NNMPC): An experimental investigation, Control Engineering Practice 19, 
454-467. 

Matlab, 2016, The Mathworks Inc., Massachusetts, USA. 
Mohd N., Aziz N., 2015, Control of Bioethanol Fermentation Process: NARX-based MPC (NARX-MPC) versus 

Linear-based MPC (LMPC), Chemical Engineering Transactions 45, 1297-1302. 
Nagy Z.K., Mahn B., Franke R., Allgower F., 2007, Evaluation study of an efficient output feedback nonlinear 

model predictive control for temperature tracking in an industrial batch reactor, Control Engineering 
Practice 15, 839-850. 

Santos R.R.C., Santos B.F., Fileti A.M.F, Da Silva F.V., Zemp R.J., 2013, Application of Artificial Neural 
Networks in an Experimental Batch Reactor of Styrene Polymerization for Predictive Model Development. 
Chemical Engineering Transactions 32, 1399-1404. 

Wehtje E., Kaur J., Adlercreutz P., Chand S., Mattiasson B., 1997, Water activity control in enzymatic 
esterification processes, Enzyme and Microbial Technology 21 (7), 502-510. 

Won K., Lee S.B., 2002, On-line conversion estimation for solvent-free enzymatic esterification systems with 
water activity control, Biotechnology and Bioprocess Engineering 7 (2), 76-84. 

Zulkeflee S.A., Abd Sata S., Aziz N., 2013, Kinetic Model with Effect of Water Content for Enzyme-Catalyzed 
Citronellyl Laurate Esterification Process, Applied Mechanics and Materials 284 -287, 423-428. 

222