Influence of Estimators and Numerical Approaches on the Implementation of NMPCs

Abstract

Advanced control strategies, together with state-estimation methods, are frequently applied to nonlinear and complex systems. It is crucial to understand which of these are the most efficient methods for the best use of these approaches in a chemical process. In the current work, nonlinear model predictive control (NMPC) approaches were developed that considered three numerical methods: single shooting (SS), multiple shooting (MS), and orthogonal collocation (OC). Their performance was compared against the Van de Vusse reactor benchmark while considering set-point changes, unreachable set-point, disturbances, and mismatches. The results showed that the NMPC based on OC presented less computational cost than the other approaches. The extended Kalman filter (EKF), constrained extended Kalman filter (CEKF), and the moving horizon estimator (MHE) were also developed. The estimators’ performance was compared for the same benchmark by considering the computational cost and the mean squared error (MSE) for the estimated variables, thereby verifying the CEKF as the best option. Finally, the performance of the nine combinations of estimators and control approaches was compared to consider the Van de Vusse reactor and the same scenarios, thereby verifying the best performance of the CEKF with the OC. The present work can help with choosing the numerical method and the estimator for controlling chemical processes.

1. Introduction

Modern technological advances have enabled the evolution of predictive strategies in several sectors of the global economy. Within this context, model predictive control (MPC) methodologies have emerged as a real possibility for various industry fields due to progress from dynamic optimization and computational areas. MPC is a widespread technique that consists of creating a mathematical model responsible for predicting process output variables at future times through system inputs/outputs and future control signals [1,2,3].

Within chemical engineering, advanced control strategies are widely studied in different fields of knowledge and have been increasing in popularity since the beginning of this century [3,4]. In this context, it is possible to highlight a diversity of works in traditional areas, such as continuous and batch chemical reactors [5,6,7,8], catalytic cracking units [9,10], distillation columns [11,12], and polymerization reactors [13,14,15]. In addition, there are also applications of the MPC for controlling crystallization processes [16,17,18,19,20], and artificial lift gas injection for oil production [21,22,23].

Therefore, it is important to highlight that MPC controllers can be divided into two groups: linear (LMPC) and nonlinear (NMPC). LMPC refers to a set of schemes in which linear models are used to predict the process dynamics and consider linear constraints on the states and inputs. The nonlinear counterpart describes a set based on nonlinear models and/or considers a non-quadratic cost function and general nonlinear constraints [24,25].

Usually, chemical engineering systems rely on mechanistic or empirical models, which present nonlinear equations and constraints in their formulations [3]. Thus, LMPC controllers might have certain limitations in dealing with nonlinear systems, thereby making NMPC controllers more suitable to deal with these characteristics.

Large-scale nonlinear systems pose challenges when it comes to obtaining reliable dynamic models. Furthermore, the online estimation of unmeasured state variables is also a complicating factor due to the need for a nonlinear estimator. In this sense, advanced control strategies combined with state-estimation techniques are very important tools in the control of nonlinear and complex systems [26,27,28,29].

To deal with the nonlinearities in each process, it is important to provide a robust control approach to handle uncertainties, wherein the control trajectory can be optimized to satisfy process constraints [4]. Process monitoring and control require that reliable information be obtained about state variables in real-time to ensure the proper operation of the industrial plant. The controller determines its actions based on these measurements, and some important variables for the process may not be provided in real time [30].

EKF, CEKF, and MHE are three of the most important applied state estimators. The first two are based on the traditional Kalman filter (KF) and have been adapted to nonlinear systems using linearization in the estimation process [31,32,33]. Unlike these recursive estimators, (EKF and CEKF), the MHE is a batch estimator based on a past moving horizon. The MHE is well-suited for systems with many states and/or measurements, but it can be computationally expensive compared to the first two counterparts [33,34]. Other state estimators are also employed, such as CEKF&S (CEKF with a smoother) [35,36], the unscented Kalman filter (UKF) [37], the cubature Kalman filter (CKF) [38], the particle filter (PF) [39], the unscented particle filter (UPF) [40], the receding nonlinear Kalman filter (RNK) [41], and variations of the EKF [42].

There are many applications of state estimators in the literature, such as the work of Tonel et al. [30] that developed a comparative study based on Kalman filter estimators (EKF and CEKF) and the MHE for a classical four-cylindrical-tanks system and nonlinear reaction systems. The authors verified a greater computational effort for processing the moving horizon estimator when compared to the EKF and CEKF due to the nonlinear conditions and disturbances attributed to the system. Furthermore, Gesthuisen et al. [31] investigated the application of a simultaneous and a sequential MHE for the polycondensation of PET and compared their performance to a CEKF.

Bohler et al. [43] successfully applied an EKF for a polymer electrolyte membrane fuel cell system. Spivey et al. [44] used the MHE and UKF for industrial process data provided by the ExxonMobil Chemical Company. They compared the performance of these estimators to the traditional EKF, thereby demonstrating better performance from the MHE and UKF. Soares et al. [45] proposed an NMPC for controlling a gas-lift oil well combined with a machine learning approach to estimate the internal state from sensor data. The proposed strategy was applied against slugging, set-point changes, and unmeasurable disturbances, and its performance was compared to an NMPC combined with an EKF. Tuveri et al. [46] did simultaneous parameter and state estimations using a MHE in a microbial experimental fed-batch process.

The choice of numerical method is also an important aspect of solving the dynamic optimization problem implied by NMPC controllers. SS, MS, and OC on finite element methods can be used for solving NMPC problems in chemical engineering with high accuracy and computational efficiency. The SS is a traditional approach in which only the control inputs are parameterized or discretized over the prediction horizon [47]. On the other hand, the MS divides the prediction horizon into multiple intervals. The initial value problem of each interval is solved independently, and matching conditions are added in the optimization problem, thereby leading to a piecewise continuous representation of the control inputs and system states [48]. The OC also discretizes the state variables based on polynomial approximation, thereby converting the dynamic optimization problem into nonlinear programming. It provides better stability compared to other approaches, particularly for open-loop unstable systems [49].

For instance, the SS method was implemented in an NMPC controller combined with an EKF estimator to operate with flash separators and solve a dynamic optimization problem [50]. Jordanou et al. [51] studied the implementation of echo state networks (ESN) for artificial lifting in oil production using electrical submersible pumps (ESP). The authors designed one strategy based on an SS-NMPC, which included the ESN model inside its structure. Diehl et al. [26] proposed an NMPC implementation based on the phenomenological modeling of a copolymerization process using the MS method. They combined a predictive controller and an MHE to address the estimation of states and parameters for the proposed dynamic system. Similarly, Roman et al. [52] used the MS method to study the combination of an NMPC and MHE controller for catalytic cracking. They only used available variables measured by the industrial plant, while other states were estimated with uncertain model parameters. Furthermore, Kirches et al. [53] also studied the application of MS methods and the NMPC for three different benchmarks by considering the existence of long prediction horizons. Rodriguez [54] implemented an NMPC controller for a pulp reactor using the OC for the discretization problem. The predictive controller was directly compared to the proportional integral (PI) controller, thereby obtaining a faster response and better performance from the quantified metrics. Furthermore, Kim et al. [28] studied the application of an NMPC and MHE structure, thereby promoting the discretization of the differential equations by the OC method to optimize the operation of batch cultures of E. coli from collected experimental data.

Even with the relevance of numerical methods’ choice for the NMPC, the studies in the literature focus on tuning, such as those works by Giraldo et al. [55], Fontes et al. [56], and Tran et al. [57]. Moreover, some studies compare the estimators’ performance, such as those works by Tonel et al. [30], Valipour et al. [34], Haseltine et al. [58], and Salau et al. [59]. However, the analysis of the combination between the NMPC numerical methods and the estimators was not considered. In this sense, this work evaluated numerical aspects of modeling a complex multivariate system by applying NMPC algorithms, such as the SS, MS, and OC, combined with state estimation techniques, such as the EKF, CEKF, and MHE. The goal of the current work is to analyze the performance and computational efficiency of these combinations of estimators and NMPC algorithms while considering set-point changes, process disturbances, and model mismatch. For this analysis, the MSE and the computational cost were considered for the six scenarios studied and all approaches used. The case studied is the Van de Vusse reactor, which is a complex system due to nonlinearities, challenging process dynamics, and multiple steady states. This work contributes to the process of selecting the most efficient numerical methods for associating NMPCs with state estimators for complex chemical engineering systems.

2. Methodology

2.1. Problem Statement

Each configuration component shown in Figure 1 was evaluated in the present paper. The dynamic optimization of the process (online) with an optimizer based on the NMPC is proposed. Here, this receding horizon strategy, implemented in a closed-loop framework, was evaluated for numerical methods based on the direct approach (direct single shooting, direct multiple shooting, and direct collocation) and state estimators (the EKF, CEKF, and MHE). In total, nine configurations were obtained by combining each technique. In the next sections, a brief methodological review of the methods is carried out. The NMPC is transformed into nonlinear programming (NLP) by parameterizing the states and controls for the direct approach. Each numerical method defines the cost function, constraints, and future errors differently. The process noise and model uncertainty are approached differently for each state estimator.

The classic time-invariant tracking control problem with a discrete-time system ($x(i+1)=f\left(x\right(i),u(i\left)\right)$) was formulated in Equation (1) and taken as a reference to address predictive control according to the objective pointed out in Section 1. The problem implemented in a closed-loop framework employs the receding control strategy to determine the vector of actions ($u\left(0\right),\cdots ,u(p-1)$). The first action ($u\left(0\right)\in U$) is implemented in the process, and the output (controlled variable) is feedback as the initial state after filtering by the estimator ($x\left(0\right)\in X$). Furthermore, the cost function comprises the Lagrange term weighted by the matrix $Q_{con}$ and control actions weighted by the matrix $R_{con}$ [60].

(1)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{u\left(0\right),\cdots ,u(p-1)}{min}& \sum _{i=1}^p\left[{(y\left(i\right)-y^{sp}\left(i\right))}^T{Q}_{con}(y\left(i\right)-y^{sp}\left(i\right))+\Delta u^T\left(i\right)R_{con}\Delta u\left(i\right)\right]\hfill \\ \hfill & s.t.x(i+1)=f(x\left(i\right),u\left(i\right)),x\left(0\right)=x_0\hfill \\ \hfill & y\left(i\right)=h\left(x\right(i),u(i\left)\right)\hfill \\ \hfill & u\left(i\right)=\Delta u\left(i\right)+u(i-1)\in U,i=1,\cdots ,p\hfill \\ \hfill & x\left(i\right)\in X,y\left(i\right)\in Y,i=1,\cdots ,p\hfill \end{array}\end{array}$

2.2. Direct Approach Strategy

2.2.1. Single Shooting

After presenting the classic MPC, the numerical method based on the direct approach (single shooting [61]) was formulated for the problem defined by Equation (2) (an adapted version of Equation (1)), which considers a control horizon (H) and a finite prediction horizon (p), as well as process constraints on states and control.

(2)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{{\left\{x\left(i\right)\right\}}_{i=k+1}^{k+p},{\left\{u\left(i\right)\right\}}_{i=k}^{k+H-1}}{min}& \sum _{i=k+1}^{k+p}\left[{(y\left(i\right)-y^{sp}\left(i\right))}^T{Q}_{con}(y\left(i\right)-y^{sp}\left(i\right))+\Delta u^T\left(i\right)R_{con}\Delta u\left(i\right)\right]\hfill \\ \hfill & s.t.x(i+1)=f\left(x\right(i),u(i\left)\right),\forall i=k,\cdots k+p-1\hfill \\ \hfill & y\left(i\right)=h\left(x\right(i),u(i\left)\right),\forall i=k+1,\cdots k+p\hfill \\ \hfill & g\left(x\right(i),u(i),y(i\left)\right)\le 0,\forall i=k+1,\cdots k+p\hfill \\ \hfill & x^l\le x\left(i\right)\le x^u,\forall i=k+1,\cdots k+p\hfill \\ \hfill & y^l\le y\left(i\right)\le y^u,\forall i=k+1,\cdots k+p\hfill \\ \hfill & u^l\le u\left(i\right)\le u^u,\forall i=k,\cdots k+p\hfill \\ \hfill & u\left(i\right)=u(i-1),\forall i=k+H,\cdots k+p\hfill \\ \hfill & \Delta u\left(i\right)=u\left(i\right)-u(i-1),\forall i=k+1,\cdots k+p\hfill \\ \hfill & x\left(k\right)=x_0,\hfill \end{array}\end{array}$

The SS method first discretizes the problem defined by Equation (2) by using a piecewise constant control representation $u(t\in [t\left(k\right),t(k+1)\left]\right)=u\left(k\right)$ over a time grid $i=k,\cdots k+p$. Then, it optimizes the problem defined by Equation (3) by treating the values of the control variables as a set of NLP variables ($w=u\left(k\right),\cdots u(k+p)$) with inequality constraints enforced on the shooting nodes (c) [60,62,63]. In addition, the model ($m(w,x_0,t)$) is solved by an integration routine for each control in w over the entire time grid, where the trajectories $x\left(k\right),\cdots x(k+p)$ are dependent variables.

(3)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{w}{min}& \varphi \left(J(m(w,x_0,t),w)\right)\hfill \\ \hfill & s.t.c(m(w,x_0,t\left(k\right)),w\left(k\right))\le 0\hfill \end{array}\end{array}$

SS is a method that results in low-dimension NLP. However, it strongly depends on the initial values and initial guesses for the control parameters, which can be challenging to obtain for open-loop unstable systems [60]. In addition, in the present paper, the open-source CasADi [64] software for nonlinear optimization and algorithmic differentiation was employed with the initial value problems (ODE) solved with the classic Runge–Kutta method (RK4) instead of the solver interface CVODES/IDAS, which presented a prohibitive computational time. In addition, the open-source primal-dual interior point method (IPOPT) was implemented in CasADi and used as optimization algorithm to solve the NLP problem.

2.2.2. Multiple Shooting

MS is an extension of SS that directly adds process state variables at the matching conditions to the NLP ($w=x\left(k\right),\cdots x(k+p),u\left(k\right),\cdots u(k+p)$), thereby breaking down the system integration into short time intervals. It is a lifted SS through the recursion in Equation (4). The NLP decision variables are reformulated in Equation (5), where equality (g) and inequality (c) constraints are enforced on the shooting nodes [65].

(4)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & x\left(1\right)-f\left(x\right(0),u(0\left)\right)=0\hfill \\ \hfill & x\left(2\right)-f\left(x\right(1),u(1\left)\right)=0\hfill \\ \hfill & \cdots \hfill \\ \hfill & x(k+p)-f\left(x\right(k+p-1),u(k+p-1\left)\right)=0\hfill \end{array}\end{array}$

(5)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{w}{min}& \varphi \left(J\right(w\left)\right)\hfill \\ \hfill s.t.& g\left(w\right)=\left[\begin{array}{c}x\left(0\right)-x_0\\ x\left(1\right)-f\left(x\right(0),u(0\left)\right)\\ \cdots \\ x(k+p)-f\left(x\right(k+p-1),u(k+p-1\left)\right)\end{array}\right]=0\hfill \\ \hfill & c\left(w\right)=\left[\begin{array}{c}c\left(x\right(0),u(0\left)\right)\\ \cdots \\ c\left(x\right(k+p),u(k+p\left)\right)\\ c\left(x\right(p\left)\right)\end{array}\right]\le 0\hfill \end{array}\end{array}$

This method significantly improves the distribution of nonlinearity and numerical stability over SS methods, especially for ill-conditioned initial value problems, thereby accelerating the convergence for line search and trust region methods [63]. Ultimately, the open-source CasADi software was employed with the same specifications applied to the SS to solve this problem.

2.2.3. Direct Collocation

The direct collocation method simultaneously performs the discretization and optimization of the control problem seen in Equation (2). Unlike SS and MS, it approximates state trajectories via polynomials [66]. Therefore, they do not directly employ an integration routine. For example, the Gauss–Radau collocation defined by Equation (6) adjusts the polynomial coefficients ($\psi {\left(k\right)}_i$) so that they match the state derivative at the collocation points and the initial state of each interval. Hence, the NLP decision variables for each interval of integration are reformulated in Equation (7), with equality (g) constraints enforced on the collocation and shooting nodes (i.e., $[t{\left(k\right)}_0,\cdots ,t{\left(k\right)}_d]\in [t\left(k\right),t(k+1)]$) in addition to the inequalities constraints.

(6)$\begin{array}{c}\hfill \begin{array}{c}\hfill x(\psi \left(k\right),t)=\sum _{i=0}^d\psi {\left(k\right)}_i·P{\left(k\right)}_i\left(t\right)\\ \hfill P{\left(k\right)}_i\left(t\right)=\prod _{j=0,j\ne i}^d\frac{t-t{\left(k\right)}_j}{t{\left(k\right)}_i-t{\left(k\right)}_j}\\ \hfill P{\left(k\right)}_i\left(t{\left(k\right)}_j\right)=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}i=j\hfill \\ 0,\hfill & \mathrm{Otherwise}\hfill \end{array}\right.\end{array}\end{array}$

(7)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{w}{min}& \varphi \left(J\right(w\left)\right)\hfill \\ \hfill s.t.& g\left(w\right)=\left[\begin{array}{c}\psi {\left(0\right)}_0-x_0\\ x\left(\psi (0,t_1)\right)-\psi {\left(1\right)}_0\\ f(\psi {\left(0\right)}_i,u\left(0\right))-\sum _{j=0}^d\psi {\left(0\right)}_j·\stackrel{·}{P}{\left(0\right)}_j\left(t{\left(0\right)}_i\right)\\ \cdots \\ x(\psi \left(k\right),t(k+1))-\psi {(k+1)}_0)\\ f(\psi {\left(k\right)}_i,u\left(k\right))-\sum _{j=0}^d\psi {\left(k\right)}_j·\stackrel{·}{P}{\left(k\right)}_j\left(t{\left(k\right)}_i\right)\\ \cdots \end{array}\right]=0\hfill \\ \hfill & c\left(w\right)=\left[\begin{array}{c}c\left(x\right(0),u(0),\psi (0\left)\right)\\ \cdots \\ c\left(x\right(k+p),u(k+p),\psi (k+p\left)\right)\\ c\left(x\right(p),\psi (k+p\left)\right)\end{array}\right]\le 0\hfill \end{array}\end{array}$

Despite its inherent higher dimension, this simultaneous optimization and discretization approach is considered easier to solve than boundary value problems. NLP’s sparsity leads to potentially faster computations compared to shooting techniques [60]. In addition, the problem was solved via the open-source CasADi software with a solver interface based on the simultaneous approach as in Equation (7).

2.3. State Estimators

State estimators act in Equation (2) as state observers to predict $x_0$ from noisy measurements and models with uncertainty. Its correct estimation influences the solution of Equations (3), (5) and (7). In the present article, the difference equation model and observation model described by Equation (2) were reformulated to Equations (8) and (9), respectively, wherein the model parameters and perturbations could also be estimated. However, they do not change inside the time interval during the evolution of the model (i.e., ${\theta }_{k+1}={\theta }_k$) [67].

(8)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \tilde{x}(k+1)=\left[\begin{array}{c}x(k+1)\\ \theta (k+1)\end{array}\right]=\tilde{f}(x\left(k\right),u\left(k\right),\theta \left(k\right))+\tilde{w}\left(k\right)\hfill \\ \hfill & \tilde{f}(x\left(k\right),u\left(k\right),\theta \left(k\right))=\left[\begin{array}{c}f\left(x\right(k),u(k),\theta (k\left)\right)+w\left(k\right)\\ \theta \left(k\right)+w_{\theta }\left(k\right)\end{array}\right]\hfill \end{array}\end{array}$

(9)$y\left(k\right)=\left[\begin{array}{cc}h\left(x\right(k),u(k\left)\right)& 0\end{array}\right]\left[\begin{array}{c}x\left(k\right)\\ \theta \left(k\right)\end{array}\right]+v\left(k\right)$

where ${\theta }_k$ is the vector of the model parameters, w and $w_{\theta }$ are process noise and parameter uncertainty, respectively, and v is the vector of the measurement noise.

2.3.1. Extended Kalman Filter

The Kalman filter assumes that the random vectors $w\left(k\right)$ and $v\left(k\right)$ are defined as white, zero-mean, uncorrelated noises with known covariance matrices $Q\left(k\right)$ and $R\left(k\right)$, respectively. Thus, regarding these vectors and the initial state $\tilde{x}\left(0\right)$, it has the following properties: $E\left[\tilde{w}\left(k\right)\right]=0$ and $E\left[\tilde{w}\left(k\right)\tilde{w}{\left(k\right)}^T\right]=Q$; $E\left[v\right(k\left)\right]=0$ and $E\left[v\left(k\right)v{\left(k\right)}^T\right]=R$, $E\left[\tilde{x}\left(0\right)\right]=\tilde{\mu }\left(0\right)$, and $P\left(0\right)=E[(\tilde{x}\left(0\right)-\tilde{\mu }\left(0\right))(\tilde{x}\left(0\right)-\tilde{\mu }{\left(0\right)}^T]$ [68], where $E\left[a\right]$ is the expected value of the random variable a [68]. Since the first formulation of the Kalman filter was based on linear discrete-time systems, the extended version was adapted to non-linear dynamical systems. According to Terejanau et al. [68] and Bohler et al. [43], the estimate for the predictor (the model forecast step) that is linearized around the optimal solution (based on the Taylor series) is given in Equation (10). In Equation (11) the estimate covariance depends on the Jacobian matrix A, which is successively updated based on model predictions.

(10)$\stackrel{-}{\tilde{x}}(k+1)=\tilde{f}(\stackrel{-}{x}\left(k\right),u\left(k\right),\theta \left(k\right))$

(11)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \stackrel{-}{P}(k+1)=& A\left(k\right)\stackrel{-}{P}\left(k\right)A\left(k\right)+Q\left(k\right)\hfill \\ \hfill & A\left(k\right)=\frac{\delta \tilde{f}(\stackrel{-}{x}\left(k\right),u\left(k\right),\theta \left(k\right))}{\delta (\stackrel{-}{\tilde{x}}\left(k\right))}{|}_{\stackrel{-}{\tilde{x}}\left(k\right),u\left(k\right)}\hfill \end{array}\end{array}$

After estimating the covariance, the data assimilation step with the corrector calculates the posterior estimate of the states and parameters based on the plant’s available measurements ($y\left(k\right)$). Unlike the prediction step, the posterior estimate is determined in Equation (12) as the linear combination of $\stackrel{-}{\tilde{x}}\left(k\right)$ and $y\left(k\right)$ with the Kalman gain ($K\left(k\right)$) as the constant of proportionality.

(12)$\stackrel{\wedge }{\tilde{x}}\left(k\right)=\stackrel{-}{\tilde{x}}\left(k\right)+K\left(k\right)(y\left(k\right)-h(\stackrel{-}{x}\left(k\right),u\left(k\right))$

The covariance is determined based on the error between the steps of correction and prediction ($\stackrel{\wedge }{\tilde{x}}\left(k\right)-\stackrel{-}{\tilde{x}}\left(k\right)$), thereby resulting in Equation (13) [68].

(13)$\begin{array}{c}\hfill \begin{array}{cc}\hfill P\left(k\right)=(I-K(k)& H\left(k\right))\stackrel{-}{P}\left(k\right){(I-K\left(k\right)H\left(k\right))}^T+K\left(k\right)RK\left(k\right)\hfill \\ \hfill & H\left(k\right)=\frac{\delta h(\stackrel{-}{x}\left(k\right),u\left(k\right))}{\delta \left(\stackrel{-}{x}\left(k\right)\right)}{|}_{\stackrel{-}{x}\left(k\right),u\left(k\right)}\hfill \end{array}\end{array}$

The optimal gain defined by Equation (14) is then determined by minimizing the trace of $P\left(k\right)$ w.r.t. $K\left(k\right)$ and substituting in Equation (13) the results in Equation (15).

(14)$K\left(k\right)=\stackrel{-}{P}\left(k\right)H\left(k\right){(H\left(k\right)\stackrel{-}{P}\left(k\right)H{\left(k\right)}^T+R)}^{-1}$

(15)$P\left(k\right)=(I-K\left(k\right)H\left(k\right))\stackrel{-}{P}\left(k\right)$

In terms of implementing this method, the open-source CasADi software was also used. However, the solver interface IDAS was employed instead of the classic Runge–Kutta (RK4), because it is necessary to approximate the Jacobian matrix for the difference equation and observation model.

2.3.2. Moving Horizon Estimation

The definition of moving horizon is derived from the problem seen in Equation (2). However, according to Binder et al. [60], instead of a receding horizon strategy implemented for a regulation problem, it is an estimation problem implemented in an open loop with a finite-time horizon of size N (past measurements). This problem is defined by Muske et al. [33] as a receding horizon recursive state estimation, currently known as the MHE. In the present work, the optimization problem was reformulated in Equation (16), and the same properties for $\tilde{w}\left(k\right)$, $v\left(k\right)$, and $x\left(k\right)$ mentioned in the previous section were considered.

(16)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{{\{\tilde{w}\left(i\right),v\left(i\right)\}}_{i=k-N}^{k-1}}{min}& x\left(k\right)-\tilde{x}{\left(k\right))}^{T}P(x\left(k\right)-\tilde{x}\left(k\right))+\sum _{i=k-N}^{k-1}(\tilde{w}\left(i\right)Q\tilde{w}{\left(i\right)}^T+v\left(i\right)Rv{\left(i\right)}^T)\hfill \\ \hfill & s.t.\tilde{x}(i+1)=\tilde{f}(x\left(i\right),u\left(i\right),\theta \left(i\right)),\forall i=k-N,\cdots k-1\hfill \\ \hfill & y\left(i\right)=h\left(x\right(i),u(i\left)\right)+v\left(i\right),\forall i=k-N,\cdots k-1\hfill \\ \hfill & {\tilde{x}}^l\le \tilde{x}\left(i\right)\le {\tilde{x}}^u,\forall i=k-N,\cdots k\hfill \end{array}\end{array}$

This optimization problem is solved via the MS method, which transforms the problem defined by Equation (16) into an NLP in which the values of the manipulated variable are known for the finite-time horizon N, and the perturbations and parameters are also determined (implicitly in the form of $\tilde{x}$). However, they do not change during model evolution. Furthermore, the parameters Q and R are known a priori, and the covariance (P) is successively updated based on a priori estimates obtained from an EKF filter with the calculation of the covariance matrix updated from the Ricatti formula as defined by Equation (17) and proposed in Salau et al. [59] and Valipour and Ricardez-Sandoval [34]. In this equation, the exponential matrix is represented by $expm$ and computed with the open-source SciPy software [69].

(17)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & P\left(k\right)=Q+\varphi \left(k\right)\stackrel{-}{P}\left(k\right)\varphi \left(k\right)-\varphi \left(k\right)\stackrel{-}{P}\left(k\right)H{\left(k\right)}^T\hfill \\ \hfill & \times {[H\left(k\right)\stackrel{-}{P}\left(k\right)H{\left(k\right)}^T+R\left(k\right)]}^{-1}H\left(k\right)\stackrel{-}{P}\left(k\right)\varphi {\left(k\right)}^T\hfill \\ \hfill & \varphi \left(k\right)=expm\left(A\right(k)\Delta k)\hfill \\ \hfill & H\left(k\right)=\frac{\delta h(\stackrel{-}{x}\left(k\right),u\left(k\right))}{\delta \left(\stackrel{-}{x}\left(k\right)\right)}{|}_{\stackrel{-}{x}\left(k\right),u\left(k\right)}\hfill \end{array}\end{array}$

Ultimately, the same implementation details for the MS were used to solve this problem.

2.3.3. Constrained Extended Kalman Filter

The CEKF is an MHE with $N=0$. Then the problem defined by Equation (16) was reformulated in Equation (18), which was updated based solely on the last measurement for the EKF. As proposed in Gesthuisen et al. [31], this problem can be solved via quadratic programming (QP) when the observation model is linear.

(18)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\{\tilde{w}\left(k\right),v\left(k\right)\}}{min}& (\tilde{w}\left(k\right)Q\tilde{w}{\left(k\right)}^T+v\left(k\right)Rv{\left(k\right)}^T)\hfill \\ \hfill & s.t.\tilde{x}(k+1)=\tilde{f}(x\left(k\right),u\left(k\right),\theta \left(k\right))\hfill \\ \hfill & y\left(k\right)=h\left(x\right(k),u(k\left)\right)+v\left(k\right)\hfill \\ \hfill & {\tilde{x}}^l\le \tilde{x}\left(k\right)\le {\tilde{x}}^u\hfill \end{array}\end{array}$

This problem is formulated in Equation (19), which does not depend on a numerical integration method to compose the QP, thereby resulting in a well-conditioned problem with a reduced computational burden compared to the MHE. The state equation is solved after the optimization problem to make the state prediction for the next time step.

(19)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\{\tilde{w}\left(k\right),v\left(k\right)\}}{min}& (\psi \left(k\right)S{\left(k\right)}^{-1}\psi \left(k\right)+d^T\psi \left(k\right))\hfill \\ \hfill & s.t.\left[HI\right]\psi \left(k\right)=y\left(k\right)-h\left(x\right(k),u(k\left)\right)\hfill \\ \hfill & {\tilde{x}}^l\le \tilde{x}\left(k\right)\le {\tilde{x}}^u\hfill \\ \hfill & \psi \left(k\right)=\left[\begin{array}{c}\tilde{w}\left(k\right)\\ v\left(k\right)\end{array}\right],S\left(k\right)=\left[\begin{array}{cc}\stackrel{-}{P}\left(k\right)& 0\\ 0& R\end{array}\right],d=0\hfill \end{array}\end{array}$

The covariance was estimated based on the formulation of Equation (17). In addition, the open-source qpsolver software with the interface solver OSQP was used to solve the QP [70], wherein the covariance and Jacobian matrix were determined with CasADi as was done for the EKF.

2.4. Case Study: The Van de Vusse CSTR

The frameworks detailed in the previous subsections were evaluated using the Van de Vusse reactor as the illustrative model example, as shown in Figure 2.

It consists of a continuous stirred-tank reactor (CSTR) where the catalytic synthesis of cyclopentanol (B) from the electrophilic addition of water to cyclopentadiene (A) takes place. The reaction is said to be a cycloaddition (or Diels–Alder reaction), since it is catalyzed by an acid solution with the generation of cyclopentenediol (C) and dicyclopentadiene (D) as side products due to the high reactivity of the reagents and products [71]. The simplified scheme with the respective kinetic constants is shown in (Equation (20)):

(20)$\begin{array}{c}\hfill \begin{array}{c}\hfill A\stackrel{k_1}{⟶}B\stackrel{k_2}{⟶}C\\ \hfill 2A\stackrel{k_3}{⟶}D\end{array}\end{array}$

Some hypotheses are often considered concerning this process modeling:

• Perfect mixing;

• Constant volume;

• Constant heat capacity;

• Incompressible fluid;

• Elementary kinetics for all reactions;

• Reaction rate constants described by the Arrhenius Law;

• No inlet flow of B, C, and D.

Thus, one can obtain the following ODEs from the material balances of A and B and from the energy balance in the reactor and in the jacket (Equation (21)):

(21)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{d{C}_A}{dt}=\frac{F}{V}(C_{A,in}-C_A)-k_1\left(T_r\right)C_A-k_3\left(T_r\right)C_A^2\hfill \\ \hfill & \frac{d{C}_B}{dt}=-\frac{F}{V}{C}_B+k_1\left(T_r\right)C_A-k_2\left(T_r\right)C_B\hfill \\ \hfill & \frac{d{T}_k}{dt}=\frac{\stackrel{.}{Q}+K_w{A}_r(T_r-T_k)}{m_k{C}_{P,k}}\hfill \\ \hfill & \frac{d{T}_r}{dt}=\frac{F}{V}(T_{in}-T_r)+\frac{K_w{A}_r}{\rho C_{P}V}(T_k-T_r)+\hfill \\ \hfill & \frac{1}{\rho C_P}[k_1\left(T_r\right)C_A(-\Delta H_1)+k_2\left(T_r\right)C_B(-\Delta H_2)+k_3\left(T_r\right)C_A^2(-\Delta H_3)]\hfill \end{array}\end{array}$

where $C_A$ and $C_B$ are the concentrations of components A and B in the reactor, respectively; $C_{A,in}$ and $T_{in}$ are the concentration of A and the temperature in the inlet stream, respectively; $F/V$ is the reactor spatial velocity; and $T_r$ is the reactor temperature. $\rho$ and $C_P$ are the reaction mixture’s specific mass and heat capacity. $T_k$, $A_r$, $K_w$, $m_k$, and $C_{P,k}$ are the temperature, the heat exchange area, the heat transfer coefficient, the mass, and the heat capacity of the cooling jacket, respectively. $\Delta H_i$ and $k_i$ are the molar enthalpy and the rate constant of reaction i, of which the latter is dependent on the reactor temperature and expressed by the Arrhenius Law (Equation (22)):

(22)$k_i\left(T_r\right)=k_{0i}\exp \left[\frac{-E_{A,i}}{T_r+273.15}\right]$

where $k_{0i}$ and $E_{A,i}$ are the pre-exponential frequency factor and the activation energy for reaction i, respectively.

This system combines nonlinear behavior in several operational ranges with a relatively simple modeling, which is the reason that makes it an established benchmark for multivariable control algorithms and nonlinear optimization in the literature. The complete list of the physical–chemical parameters and their respective values is given by Table 1 [72]:

Summarizing, the present work considers the measurement vector $y={[C_A,C_B,T_r,T_k]}^T$, the input vector $u=[F/V,Q_k/(K_w$$A_r{\left)\right]}^T$, $C_{A,in}$ as an unmeasured disturbance, $T_{in}$ as a measured disturbance, and $C_B$ and $T_r$ as the controlled variables. The remainder of the physical–chemical parameters, except for $C_P$ and $k_{01}$, are considered known and certain. $C_P$ and $k_{01}$, considering uncertainty parameters, are approximated by $C_P=3.01\Delta C_P$ and $k_{01}=1.287\times {10}^{12}\Delta k_{01}$, with $\Delta C_P\in [0.5,1.5]$ and $\Delta k_{01}\in [0.5,1.5]$, which are estimated based on the methodology described in Section 2.3. Hence, the augmented state vector is $\tilde{x}\left(k\right)={[C_A,C_B,T_r,T_k,C_{A,in},\Delta k_{01},\Delta C_P]}^T$. In addition, the upper and lower bounds of $\tilde{x}\left(k\right)$ considered to compose the NMPC, CEKF, and MHE are given by Equation (23).

(23)$\left[\begin{array}{c}0\\ 0\\ 50\\ 50\\ 0.1\\ 0.5\\ 0.5\end{array}\right]\le \tilde{x}\left(k\right)\le \left[\begin{array}{c}4\\ 3\\ 200\\ 200\\ 6.1\\ 1.5\\ 1.5\end{array}\right]$

All simulations were done in a computer in LASAP/UFRJ with the following specifications: Intel Core i7-3770K, CPU 3.5 GHz, and 16 GB of RAM. A timestep of $2.5\times {10}^{-3}$h was adopted, and the simulations were done for two hours of the process. The initial applied steady state is presented in Table 2.

All scenarios considered were analyzed for set-point change, unreachable set-point, disturbances in $C_{A,in}$ and $T_{in}$, and for mismatches in the pre-exponential constant of $k_1$ and the specific heat of the mixture. Table 3 presents details of all six cases considered in the study using the different methods.

3. Results

This section shows the results with all the discussion related to it. Initially, the numerical methods used to develop the NMPC are presented, thereby comparing their computational cost and efficiency. Then, a comparison between the three estimators is made. Finally, the results for the nine combinations of estimators and numerical methods are shown, thereby demonstrating which combination is the best.

3.1. Numerical Methods Analysis

The absolute and relative tolerances were ${10}^{-8}$ for the SS, MS, and OC. For the OC, three collocation points were used in each time interval. The prediction and control horizons were tuned as 40 and 10 for all cases analyzed in this study based on Seborg et al. [73]. Furthermore, the weights of the controlled variable deviations were tuned as 5 for $C_B$ and 0.05 for $T_r$, while the suppression factors for manipulated variables were defined as ${10}^{-5}$ for $F⁄V$ and $9\times {10}^{-3}$ for $Q⁄(K_w$ $A_r)$. The difference the manipulated variables in each time step was constrained as represented in Equations (24) and (25), where $u_1$ is $F⁄V$ and $u_2$ is $Q⁄(K_w$ $A_r)$.

(24)$-50h^{-1}\le \Delta u_1\left(t\right)\le 50h^{-1}$

(25)$-0.15K\le \Delta u_2\left(t\right)\le 0.15K$

In the current analysis, the simulations were done using the CEKF estimator, and all results are presented in the Supplementary Material. As examples, Figure 3 exposes the results for all three NMPC approaches considering changes in the set-points of 1.1 to 0.8 mol/L and 140 to 135 °C. Moreover, Figure 4 illustrates the performance of these methods considering two pulses in $T_{in}$ equal to 150 °C in 0.5 h and 120 °C in 1.5 h. In all cases, the set-points could be reached according to all the constraints imposed on all variables. Furthermore, the different approaches analyzed presented similar results, which is expected, because, even though different numerical methods were considered, the same problem was solved.

Table 4 presents the MSE calculated for all cases analyzed and considers the difference between the controlled variables and the set-points. These results confirm that all three approaches presented similar performances in all cases studied.

The numerical methods used to develop the NMPC presented similar control results, but the major difference between these approaches is related to the computational cost, as is shown in Table 5. The results presented in Table 5 show that the controller based on SS took around 70 min to run each simulation. In comparison, this was done in about 4 min for the controller based on the OC. For the controller based on the MS, each simulation ran in around 20 min, which was an intermediate case compared to the previous ones. Therefore, the OC showed a considerable advantage compared to the other methods considered in this study.

3.2. Estimator Analysis

In the current work, the performance of the EKF, CEKF, and MHE was also studied to estimate $C_A,C_B,T_r,T_k,C_{A,in}$, the pre-exponential constant of $k_1$, and the specific heat of the mixture. The capacity of these estimators was analyzed for the same cases presented in Table 3 and for the controller based on the OC. For the MHE, the horizon was tuned to three. The initial covariance ($P_0$), model uncertainties (Q), and the measurement uncertainties (R) were defined according to Equations (26)–(28) for all cases.

(26)$P_0=[{10}^3,{10}^3,{10}^3,{10}^3,{10}^3,{10}^3,{10}^3]$

(27)$Q=[3\times {10}^{-6},5\times {10}^{-6},{10}^{-3},{10}^{-3},{10}^{-3},{10}^{-2},{10}^{-2}]$

(28)$R=[3\times {10}^{-2},5\times {10}^{-3},3\times {10}^{-1},8\times {10}^{-1}]$

The CEKF and MHE consider constraints in their formulation. The same limits previously defined for these estimators were imposed for the state and controlled variables. Moreover, the limits of −1 and 1 were imposed for the parameters $\tilde{w}(k+1)$ and $v\left(k\right)$.

Figure 5 and Figure 6 show results for the estimation of the pre-exponential constant of $k_1$ ($k_{01}$) and the specific heat of the mixture ($C_p$). The results from Figure 5 are from Case 5, which considered a difference of 5% in the value of the pre-exponential constant of $k_1$. The results from Figure 6 are from the scenario that considered a mismatch in $C_p$. All three methods could efficiently estimate the desired parameters.

Figure 7 presents the MSE values for estimations of all seven variables and the three methods used. These results show that the three estimators presented similar performances for all variables. However, the CEKF and MHE presented the advantage of considering constraints in their formulation, while the EKF did not consider them. For the cases analyzed in this work, the variables did not violate the constraints, thereby making estimations close to the simulated values of all the variables.

The time demanded to run each simulation was very close, as seen in Table 6. The CEKF took less time to run most simulations, while the MHE took more time than all three. The MHE has the most complex structure that considers a horizon to make the estimations, which demands more computational cost. Therefore, the estimators showed very similar performances for the predicted values and the computational cost, thereby observing the necessity to evaluate the estimators’ formulation in addition to the computational cost and the MSE. Considering theses points, the CEKF is the best option because it presents a more straightforward formulation than the MHE while demanding less computational effort. Moreover, the CEFK considers constraints in its formulation, while the EKF does not. Alternatively, the CEKF&S [35] could be used to obtain the advantages of the MHE regarding the filtering and the advantages of the CEKF regarding the low computational cost.

3.3. Combinations Analysis

The last analysis consisted of comparing each combination of the estimator and NMPC approach. The nine cases considered are shown in Table 7.

When comparing the computational cost, it is evident that combinations considering the NMPC based on the SS took much more time to run than the other numerical methods, as seen in Table 8. Therefore, Combinations 1, 2, and 3 took around 70 min to run each simulation, while the combinations considering the OC ran it in about 4 min. Combinations 4, 5, and 6, which applied the NMPC based on the MS, took around 20 min to run each simulation, which was an intermediate time between the other two. Because of the disadvantage of Combinations 1, 2, and 3, the combinations considering the SS presented the worst performances compared to the others.

Considering the MSE values for the estimators’ predictions and the performance of the NMPC approaches to achieve the set-points, the results for all nine combinations are very similar, as seen in Figure 8 and Figure 9. In most cases analyzed, the MSE values were almost the same, wherein we observed a few cases presenting small differences in the MSE values. Therefore, all combinations presented efficient performances when considering mismatches, set-point changes, disturbances, and unreachable set-point situations.

All combinations could efficiently achieve these goals by considering the estimated values and reaching the set points. Accounting for the computational cost, the NMPC based on the OC presented a significant advantage, as it ran the simulations in less time than the other methods. Moreover, the CEKF had some advantages compared to the other two estimators, thereby making it the best choice. Therefore, the best combination used the CEKF estimator with the NMPC based on the OC. An example of the efficiency of this approach can be observed in Figure 10, where it was applied to maintain the set-points of 140 °C and 1.1 mol/L for a change in $C_{a,in}$.

4. Conclusions

In the current work, NMPCs were developed using SS, MS, and OC numerical methods. Their performance was compared using the Van de Vusse benchmark while considering set-point changes, unreachable set-point, disturbances in $C_{A,in}$ and $T_{in}$, and mismatches in $k_1$ and $C_P$. The performance of the three control approaches for reaching the set-points was evaluated for the MSE, thereby verifying a similar behaviour for all scenarios. The computational cost was also evaluated for the controllers’ performance, wherein we observed that the NMPC based on the OC took around 4 min to run the simulations, while the controllers based on the SS and MS ran it in about 70 and 20 min, respectively.

The MHE, EKF, and CEKF estimators were also developed, and their performance was compared for the same cases and analysis of computational cost and MSE. All three approaches could efficiently estimate $C_A,C_B,T_r,T_k,C_{A,in}$, the preexponential constant of $k_1$, and the specific heat of the mixture, thereby achieving acceptable and similar values of MSE. Moreover, all three methods presented similar computational costsand took around 4 min to run each simulation. For the Van de Vuse reactor and the analyzed scenarios, the estimators presented similar values of computational cost and MSE. Therefore, it was necessary to evaluate the estimators while also considering their formulation. The EKF presented the disadvantage of not considering constraints in its formulation, which differed from the CEKF and MHE. Furthermore, the CEKF was a better option than the MHE due to its simpler formulation.

Finally, nine different combinations of estimators and NMPC approaches were compared for the developed techniques. Their performance was also tested for the Van de Vusse reactor for the same previous scenarios and analysis. The combinations could efficiently achieve the set-points and estimate the desired variables in all scenarios, thereby presenting acceptable MSE values for the estimated and controlled variables. However, the combination of the CEKF and OC showed the best performance due to the lower computational cost of the OC and the better CEKF formulation. These results indicate a strong potential for the use of the OC combined with the CEKF in chemical processes. The current work presents a detailed performance analysis of the numerical methods for NMPC and their combinations with estimators. This study can help when it comes to choosing the approach for implementing NMPCs.

Author Contributions

As the article’s main author, F.A.R.D.L. participated in developing the CEKF and MHE codes, running all simulations, and writing in all sections. R.d.R.F. participated in developing single shooting, multiple shooting, CEKF, and MHE libraries, writing the Methodology section, and paper formatting. R.C. developed the base code in CasADi, which includes the ODE model, NMPC, and EKF codes, and contributed to the Methodology section. M.C.F.C. participated in the bibliographic review and development of the Introduction. C.A.G.E. participated in the bibliographic review and development of the Introduction. For review and editing, all authors participated; conceptualization and supervision, M.B.d.S.J. and A.R.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Predictive control configuration for the direct approach and state estimators.

Figure 2. Van de Vusse reactor scheme.

Figure 3. Results for the three NMPC approaches for Case 1, considering set-point changes.

Figure 4. Results for the three NMPC approach for Case 4, considering changes in [Forumla omitted. See PDF.].

Figure 5. [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] estimations for Case 5, considering mismatch in [Forumla omitted. See PDF.].

Figure 6. [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] estimations for Case 6, considering mismatch in [Forumla omitted. See PDF.].

Figure 7. MSE values for estimations of the seven variables for the three approaches.

Figure 8. Mean squared error values calculated for estimations of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], and the set-points of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for the studied cases and the nine combinations.

Figure 9. Mean squared error values calculated for estimations of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] for the studied cases and the nine combinations.

Figure 10. Results for Combination 8 for maintaining the set-points of 140 °C and 1.1 mol/L, as well as a change in [Forumla omitted. See PDF.].

Table 1

Van de Vusse CSTR parameters.

Parameters	Value
$k_{01}$ (h−1)	$1.287\times {10}^{12}$
$k_{02}$ (h−1)	$1.287\times {10}^{12}$
$k_{03}$ (L/mol L)	$9.043\times {10}^9$
$E_{A,1}$ (K)	$9758.3$
$E_{A,2}$ (K)	$9758.3$
$E_{A,3}$ (K)	8560
$-\Delta H_1$ (kJ/mol A)	$-4.20$
$-\Delta H_2$ (kJ/mol B)	11
$-\Delta H_3$ (kJ/mol C)	$41.85$
$\rho$ (kg/L)	$0.9342$
$C_P$ (kJ/Kg K)	$3.01$
$C_{P,k}$ (kJ/Kg K)	2
$K_w$ (kJ/h K m2)	$4032.0$
$A_r$ (m2)	$0.215$
$m_k$ (kg)	5
V (L)	$10.01$

Table 2

Initial conditions used in the simulations.

${\mathit{C}}_{\mathit{A}}$	${\mathit{C}}_{\mathit{B}}$	${\mathit{T}}_{\mathit{r}}$	${\mathit{T}}_{\mathit{k}}$	$\mathit{F}/\mathit{V}$	$\mathit{Q}/\left({\mathit{K}}_{\mathit{w}}{\mathit{A}}_{\mathit{r}}\right)$	${\mathit{C}}_{\mathit{A},\mathbf{in}}$	${\mathit{T}}_{\mathbf{in}}$
(mol/L)	(mol/L)	(°C)	(°C)	(h−1)	$\left(\mathit{K}\right)$	(mol/L)	(°C)
2.4	1.1	140	140	85	−0.04	5.1	130

Table 3

Cases considered to analyze the different approaches proposed.

Case 1	Set-point changes from 1.1 to 0.8 mol/L and 140 to 135 °C in one hour of process.
Case 2	Unreachable set-point, maintaining the set-point for $T_r$ in 140 °C, and changing the set-point for $C_B$ from 1.1 to 1.3 mol/L in one hour of process.
Case 3	Application of a disturbance in $C_{A,in}$, maintaining the set-points of 1.1 mol/L and 140 °C. This value was changed from 5.1 to 6.0 mol/L in one hour of process.
Case 4	Application of disturbances in $T_{in}$, maintaining the set-points of 1.1 mol/L and 140 °C. For this case, two pulses were considered, one equals to 150 °C in 0.5 h, and the other equals 120 °C in 1.5 h. Each pulse lasted one sampling time.
Case 5	Application of a mismatch, imposing a difference of 5% in the value of the pre-exponential constant of $k_1$.
Case 6	Application of a mismatch, imposing a difference of 20% in the value of the specific heat of the mixture.

Table 4

Mean squared error values for the three NMPC methods and all cases analyzed.

Case	SS		MS		CO
	${\mathit{C}}_{\mathit{B}}$	${\mathit{T}}_{\mathit{r}}$	${\mathit{C}}_{\mathit{B}}$	${\mathit{T}}_{\mathit{r}}$	${\mathit{C}}_{\mathit{B}}$	${\mathit{T}}_{\mathit{r}}$
1	$3.44\times {10}^{-4}$	$6.57\times {10}^{-2}$	$3.55\times {10}^{-4}$	$6.57\times {10}^{-2}$	$3.16\times {10}^{-4}$	$6.87\times {10}^{-2}$
2	$1.65\times {10}^{-2}$	$1.03\times {10}^{-2}$	$1.66\times {10}^{-2}$	$9.42\times {10}^{-3}$	$1.68\times {10}^{-2}$	$3.01\times {10}^{-2}$
3	$1.16\times {10}^{-4}$	$8.58\times {10}^{-3}$	$9.80\times {10}^{-5}$	$7.59\times {10}^{-3}$	$1.16\times {10}^{-4}$	$9.29\times {10}^{-3}$
4	$1.98\times {10}^{-3}$	$1.79$	$1.85\times {10}^{-3}$	$1.82$	$1.65\times {10}^{-3}$	$1.71$
5	$2.37\times {10}^{-3}$	$2.14\times {10}^{-2}$	$2.43\times {10}^{-3}$	$2.09\times {10}^{-2}$	$2.43\times {10}^{-3}$	$2.45\times {10}^{-2}$
6	$9.99\times {10}^{-4}$	$1.30\times {10}^{-1}$	$1.01\times {10}^{-3}$	$1.28\times {10}^{-1}$	$1.23\times {10}^{-3}$	$1.12\times {10}^{-1}$

Table 5

Time needed to realize each simulation for NMPCs based on the different numerical methods.

Case	${\mathit{t}}_{\mathit{SS}}$ (min)	${\mathit{t}}_{\mathit{MS}}$ (min)	${\mathit{t}}_{\mathit{OC}}^{*}$ (min)
1	69.18	20.13	4.02
2	68.97	19.18	3.54
3	68.89	19.37	3.90
4	68.98	19.82	3.97
5	75.13	23.63	3.61
6	73.05	19.55	4.02

^* Best scenario.

Table 6

Time spent to run each simulation for the three estimators.

Case	${\mathit{t}}_{\mathit{EKF}}$ (min)	${\mathit{t}}_{\mathit{CEKF}}$ (min)	${\mathit{t}}_{\mathit{MHE}}$ (min)
1	4.26	4.02	4.52
2	3.66	3.54	3.94
3	3.96	3.90	4.27
4	3.93	3.97	4.12
5	3.83	3.61	3.94
6	4.20	4.02	4.31

Table 7

Combination of estimators and NMPC methods considered in the analysis.

Combination 1	SS and EKF
Combination 2	SS and CEKF
Combination 3	SS and MHE
Combination 4	MS and EKF
Combination 5	MS and CEKF
Combination 6	MS and MHE
Combination 7	OC and EKF
Combination 8	OC and CEKF
Combination 9	OC and MHE

Table 8

Computational time demanded in minutes to run each simulation for the nine combinations and the six cases.

Case	1	2	3	4	5	6
$t_{Comb,1}$	71.01	67.88	67.72	70.61	74.95	71.93
$t_{Comb,2}$	69.18	68.97	68.89	68.98	75.13	73.05
$t_{Comb,3}$	69.57	67.34	70.16	70.61	74.25	73.39
$t_{Comb,4}$	19.94	18.89	21.05	20.41	22.57	19.38
$t_{Comb,5}$	20.13	19.18	19.37	19.82	23.63	19.55
$t_{Comb,6}$	20.30	20.44	20.91	20.90	23.79	19.79
$t_{Comb,7}^{*}$	4.26	3.66	3.96	3.93	3.83	4.20
$t_{Comb,8}^{*}$	4.02	3.54	3.90	3.97	3.61	4.02
$t_{Comb,9}^{*}$	4.52	3.94	4.27	4.12	3.94	4.31

^* Best scenarios.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/pr11041102/s1, Figure S1: Simulation for combination 1 and case 1, considering set-point changes; Figure S2: Simulation for combination 1 and case 2, considering unreachable set-point; Figure S3: Simulation for combination 1 and case 3, considering a change in $C_{A,in}$; Figure S4: Simulation for combination 1 and case 4, considering changes in $T_{in}$; Figure S5: Simulation for combination 1 and case 5, considering mismatch in $k_{01}$; Figure S6: Simulation for combination 1 and case 6, considering mismatch in $C_P$; Figure S7: Simulation for combination 2 and case 1, considering set-point changes; Figure S8: Simulation for combination 2 and case 2, considering unreachable set-point; Figure S9: Simulation for combination 2 and case 3, considering a change in $C_{A,in}$; Figure S10: Simulation for combination 2 and case 4, considering changes in $T_{in}$; Figure S11: Simulation for combination 2 and case 5, considering mismatch in $k_{01}$; Figure S12: Simulation for combination 2 and case 6, considering mismatch in $C_P$; Figure S13: Simulation for combination 3 and case 1, considering set-point changes; Figure S14: Simulation for combination 3 and case 2, considering unreachable set-point; Figure S15: Simulation for combination 3 and case 3, considering a change in $C_{A,in}$; Figure S16: Simulation for combination 3 and case 4, considering changes in $T_{in}$; Figure S17: Simulation for combination 3 and case 5, considering mismatch in $k_{01}$; Figure S18: Simulation for combination 3 and case 6, considering mismatch in $C_P$; Figure S19: Simulation for combination 4 and case 1, considering set-point changes; Figure S20: Simulation for combination 4 and case 2, considering unreachable set-point; Figure S21: Simulation for combination 4 and case 3, considering a change in $C_{A,in}$; Figure S22: Simulation for combination 4 and case 4, considering changes in $T_{in}$; Figure S23: Simulation for combination 4 and case 5, considering mismatch in $k_{01}$; Figure S24: Simulation for combination 4 and case 6, considering mismatch in $C_P$; Figure S25: Simulation for combination 5 and case 1, considering set-point changes; Figure S26: Simulation for combination 5 and case 2, considering unreachable set-point; Figure S27: Simulation for combination 5 and case 3, considering a change in $C_{A,in}$; Figure S28: Simulation for combination 5 and case 4, considering changes in $T_{in}$; Figure S29: Simulation for combination 5 and case 5, considering mismatch in $k_{01}$; Figure S30: Simulation for combination 5 and case 6, considering mismatch in $C_P$. Figure S31: Simulation for combination 6 and case 1, considering set-point changes; Figure S32: Simulation for combination 6 and case 2, considering unreachable set-point; Figure S33: Simulation for combination 6 and case 3, considering a change in $C_{A,in}$; Figure S34: Simulation for combination 6 and case 4, considering changes in $T_{in}$; Figure S35: Simulation for combination 6 and case 5, considering mismatch in $k_{01}$; Figure S36: Simulation for combination 6 and case 6, considering mismatch in $C_P$; Figure S37: Simulation for combination 7 and case 1, considering set-point changes; Figure S38: Simulation for combination 7 and case 2, considering unreachable set-point; Figure S39: Simulation for combination 7 and case 3, considering a change in $C_{A,in}$; Figure S40: Simulation for combination 7 and case 4, considering changes in $T_{in}$; Figure S41: Simulation for combination 7 and case 5, considering mismatch in $k_{01}$; Figure S42: Simulation for combination 7 and case 6, considering mismatch in $C_P$; Figure S43: Simulation for combination 8 and case 1, considering set-point changes; Figure S44: Simulation for combination 8 and case 2, considering unreachable set-point; Figure S45: Simulation for combination 8 and case 3, considering a change in $C_{A,in}$; Figure S46: Simulation for combination 8 and case 4, considering changes in $T_{in}$; Figure S47: Simulation for combination 8 and case 5, considering mismatch in $k_{01}$; Figure S48: Simulation for combination 8 and case 6, considering mismatch in $C_P$; Figure S49: Simulation for combination 9 and case 1, considering set-point changes; Figure S50: Simulation for combination 9 and case 2, considering unreachable set-point; Figure S51: Simulation for combination 9 and case 3, considering a change in $C_{A,in}$; Figure S52: Simulation for combination 9 and case 4, considering changes in $T_{in}$; Figure S53: Simulation for combination 9 and case 5, considering mismatch in $k_{01}$; Figure S54: Simulation for combination 9 and case 6, considering mismatch in $C_P$.

References

1. Rawlings, J.B.; Mayne, D.Q.; Diehl, M. Model Predictive Control: Theory, Computation, and Design; Nob Hill Publishing: Madison, WI, USA, 2017; Volume 2.

2. Morato, M.M.; Normey-Rico, J.E.; Sename, O. Model predictive control design for linear parameter varying systems: A survey. Annu. Rev. Control; 2020; 49, pp. 64-80. [DOI: https://dx.doi.org/10.1016/j.arcontrol.2020.04.016]

3. Kumar, A.S.; Ahmad, Z. Model Predictive Control (MPC) and Its Current Issues in Chemical Engineering. Chem. Eng. Commun.; 2012; 199, pp. 472-511. [DOI: https://dx.doi.org/10.1080/00986445.2011.592446]

4. Qin, S.J.; Badgwell, T.A. A survey of industrial model predictive control technology. Control Eng. Pract.; 2003; 11, pp. 733-764. [DOI: https://dx.doi.org/10.1016/S0967-0661(02)00186-7]

5. Joy, P.; Schultz, E.S.; Ebrahimi, F.; Turan, U.; Casteel, S.; Schaffrath, T.; Hammen, R.; Mhamdi, A. Dynamic optimization and nonlinear model predictive control of a semi-batch epoxidation process. J. Process Control; 2021; 108, pp. 55-67. [DOI: https://dx.doi.org/10.1016/j.jprocont.2021.10.013]

6. Che, Y.; Zhao, Z.; Wang, Z.; Liu, F. Iterative learning model predictive control for multivariable nonlinear batch processes based on dynamic fuzzy PLS model. J. Process Control; 2022; 119, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.jprocont.2022.09.005]

7. Chen, X.; Heidarinejad, M.; Liu, J.; Christofides, P.D. Distributed economic MPC: Application to a nonlinear chemical process network. J. Process Control; 2012; 22, pp. 689-699. [DOI: https://dx.doi.org/10.1016/j.jprocont.2012.01.016]

8. Ramaswamy, S.; Cutright, T.; Qammar, H. Control of a continuous bioreactor using model predictive control. Process Biochem.; 2005; 40, pp. 2763-2770. [DOI: https://dx.doi.org/10.1016/j.procbio.2004.12.019]

9. Roman, R.; Nagy, Z.K.; Allgöwer, F.; Agachi, S.P. Dynamic modeling and nonlinear model predictive control of a fluid catalytic cracking unit. Proceedings of the European Symposium on Computer-Aided Process Engineering-15, 38th European Symposium of the Working Party on Computer Aided Process Engineering; Barcelona, Spain, 29 May–1 June 2005; Computer Aided Chemical Engineering Volume 20, pp. 1363-1368. [DOI: https://dx.doi.org/10.1016/S1570-7946(05)80069-0]

10. Cristea, M.V.; Agachi, Ş.P.; Marinoiu, V. Simulation and model predictive control of a UOP fluid catalytic cracking unit. Chem. Eng. Process. Process Intensif.; 2003; 42, pp. 67-91. [DOI: https://dx.doi.org/10.1016/S0255-2701(02)00017-X]

11. Abou-Jeyab, R.; Gupta, Y.; Gervais, J.; Branchi, P.; Woo, S. Constrained multivariable control of a distillation column using a simplified model predictive control algorithm. J. Process Control; 2001; 11, pp. 509-517. [DOI: https://dx.doi.org/10.1016/S0959-1524(00)00029-9]

12. Alpbaz, M.; Karacan, S.; Cabbar, Y.; Hapoğlu, H. Application of model predictive control and dynamic analysis to a pilot distillation column and experimental verification. Chem. Eng. J.; 2002; 88, pp. 163-174. [DOI: https://dx.doi.org/10.1016/S1385-8947(01)00304-7]

13. Alhamad, B.; Romagnoli, J.A.; Gomes, V.G. On-line multi-variable predictive control of molar mass and particle size distributions in free-radical emulsion copolymerization. Chem. Eng. Sci.; 2005; 60, pp. 6596-6606. [DOI: https://dx.doi.org/10.1016/j.ces.2005.05.053]

14. Bindlish, R. Nonlinear model predictive control of an industrial process with steady-state gain inversion. Comput. Chem. Eng.; 2020; 135, 106739. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2020.106739]

15. Mokhtarname, R.; Safavi, A.A.; Urbas, L.; Salimi, F.; Zerafat, M.M.; Harasi, N. Model development and control of an auto-refrigerated polystyrene polymerization reactor. Trans. Inst. Meas. Control; 2021; 43, pp. 3456-3472. [DOI: https://dx.doi.org/10.1177/01423312211025125]

16. Damour, C.; Benne, M.; Boillereaux, L.; Grondin-Perez, B.; Chabriat, J.P. NMPC of an industrial crystallization process using model-based observers. J. Ind. Eng. Chem.; 2010; 16, pp. 708-716. [DOI: https://dx.doi.org/10.1016/j.jiec.2010.07.014]

17. Moraes, M.G.; Souza, M.B., Jr.; Secchi, A.R. Dynamics and MPC of an Evaporative Continuous Crystallization Process. Proceedings of the 28th European Symposium on Computer Aided Process Engineering; Graz, Austria, 10–13 June 2018; Computer Aided Chemical Engineering Volume 43, pp. 997-1002. [DOI: https://dx.doi.org/10.1016/B978-0-444-64235-6.50175-3]

18. Lima, F.A.R.D.; de Moraes, M.G.F.; Secchi, A.R.; de Souza, M.B., Jr. Development of a recurrent neural networks-based NMPC for controlling the concentration of a crystallization process. Digit. Chem. Eng.; 2022; 5, 100052. [DOI: https://dx.doi.org/10.1016/j.dche.2022.100052]

19. Hermanto, M.W.; Braatz, R.D.; Chiu, M.S. Integrated batch-to-batch and nonlinear model predictive control for polymorphic transformation in pharmaceutical crystallization. AIChE J.; 2011; 57, pp. 1008-1019. [DOI: https://dx.doi.org/10.1002/aic.12331]

20. Nagy, Z.K.; Braatz, R.D. Robust nonlinear model predictive control of batch processes. AIChE J.; 2003; 49, pp. 1776-1786. [DOI: https://dx.doi.org/10.1002/aic.690490715]

21. Delou, P.A.; de Azevedo, J.P.; Krishnamoorthy, D.; de Souza, M.B.; Secchi, A.R. Model Predictive Control with Adaptive Strategy Applied to an Electric Submersible Pump in a Subsea Environment. IFAC PapersOnLine; 2019; 52, pp. 784-789. [DOI: https://dx.doi.org/10.1016/j.ifacol.2019.06.157]

22. Dias, A.C.S.R.; Soares, F.R.; Jäschke, J.; de Souza, M.B.; Pinto, J.C. Extracting Valuable Information from Big Data for Machine Learning Control: An Application for a Gas Lift Process. Processes; 2019; 7, 252. [DOI: https://dx.doi.org/10.3390/pr7050252]

23. Miyoshi, S.C.; Nunes, M.; Salles, A.; Secchi, A.R.; de Souza, M.B.; Brandão, A.L. Nonlinear model predictive control application for gas-lift based oil production. Proceedings of the 28th European Symposium on Computer Aided Process Engineering; Graz, Austria, 10–13 June 2018; Computer Aided Chemical Engineering Volume 43, pp. 1177-1182. [DOI: https://dx.doi.org/10.1016/B978-0-444-64235-6.50205-9]

24. Allgower, F.; Findeisen, R.; Nagy, Z. Nonlinear model predictive control: From theory to application. J. Chin. Inst. Chem. Eng.; 2004; 35, pp. 299-315.

25. Shettigar, J.P.; Lochan, K.; Jeppu, G.; Palanki, S.; Indiran, T. Development and Validation of Advanced Nonlinear Predictive Control Algorithms for Trajectory Tracking in Batch Polymerization. ACS Omega; 2021; 6, pp. 22857-22865. [DOI: https://dx.doi.org/10.1021/acsomega.1c03386]

26. Diehl, M.; Kühl, P.; Bock, H.; Schlöder, J.; Mahn, B.; Kallrath, J. Combined nonlinear model predictive control and moving horizon estimation for a copolymerization process. Proceedings of the 16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process Systems Engineering; Garmisch-Partenkirchen, Germany, 9–13 July 2006; Computer Aided Chemical Engineering Volume 21, pp. 1527-1532. [DOI: https://dx.doi.org/10.1016/S1570-7946(06)80264-6]

27. Hebing, L.; Tran, F.; Brandt, H.; Engell, S. Robust Optimizing Control of Fermentation Processes Based on a Set of Structurally Different Process Models. Ind. Eng. Chem. Res.; 2020; 59, pp. 2566-2580. [DOI: https://dx.doi.org/10.1021/acs.iecr.9b05504]

28. Kim, J.W.; Krausch, N.; Aizpuru, J.; Barz, T.; Lucia, S.; Neubauer, P.; Cruz Bournazou, M.N. Model predictive control and moving horizon estimation for adaptive optimal bolus feeding in high-throughput cultivation of E. coli. Comput. Chem. Eng.; 2023; 172, 108158. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2023.108158]

29. de Azevedo Delou, P.; Curvelo, R.; de Souza, M.B.; Secchi, A.R. Steady-state real-time optimization using transient measurements in the absence of a dynamic mechanistic model: A framework of HRTO integrated with Adaptive Self-Optimizing IHMPC. J. Process Control; 2021; 106, pp. 1-19. [DOI: https://dx.doi.org/10.1016/j.jprocont.2021.08.013]

30. Tonel, G.; Salau, N.P.; Trierweiler, J.O.; Secchi, A.R. Comprehensive evaluation of EKF, CEKF, and Moving Horizon estimators for on-line processes applications. Proceedings of the 18th European Symposium on Computer Aided Process Engineering; Lyon, France, 1–4 June 2008.

31. Gesthuisen, R.; Klatt, K.U.; Engell, S. Optimization-based state estimation — A comparative study for the batch polycondensation of polyethyleneterephthalate. Proceedings of the 2001 European Control Conference (ECC); Porto, Portugal, 4–7 September 2001; pp. 1062-1067. [DOI: https://dx.doi.org/10.23919/ECC.2001.7076055]

32. Julier, S.J.; Uhlmann, J.K. New extension of the Kalman filter to nonlinear systems. Proceedings of the Signal Processing, Sensor Fusion, and Target Recognition VI; Orlando, FL, USA, 21 April 1997; SPIE: New York, NY, USA, 1997; Volume 3068, pp. 182-193.

33. Muske, K.R.; Rawlings, J.B.; Lee, J.H. Receding Horizon Recursive State Estimation. Proceedings of the 1993 American Control Conference; San Francisco, CA, USA, 2–4 June 1993; pp. 900-904. [DOI: https://dx.doi.org/10.23919/ACC.1993.4792993]

34. Valipour, M.; Ricardez-Sandoval, L.A. Assessing the Impact of EKF as the Arrival Cost in the Moving Horizon Estimation under Nonlinear Model Predictive Control. Ind. Eng. Chem. Res.; 2021; 60, pp. 2994-3012. [DOI: https://dx.doi.org/10.1021/acs.iecr.0c06095]

35. Salau, N.P.; Trierweiler, J.O.; Secchi, A.R. State estimators for better bioprocesses operation. Proceedings of the 22nd European Symposium on Computer Aided Process Engineering; Eindhoven, The Netherlands, 16–19 June 2012; Computer Aided Chemical Engineering Bogle, I.D.L.; Fairweather, M. Elsevier: Amsterdam, The Netherlands, 2012; Volume 30, pp. 1267-1271. [DOI: https://dx.doi.org/10.1016/B978-0-444-59520-1.50112-3]

36. Silva, D.M.; Secchi, A.R. Recursive state and parameter estimation of COVID-19 circulating variants dynamics. Sci. Rep.; 2022; 12, 15879. [DOI: https://dx.doi.org/10.1038/s41598-022-18208-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36151226]

37. Julier, S.; Uhlmann, J. Unscented filtering and nonlinear estimation. Proc. IEEE; 2004; 92, pp. 401-422. [DOI: https://dx.doi.org/10.1109/JPROC.2003.823141]

38. Arasaratnam, H.; Haykin, S. Cubature Kalman Filters. IEEE Trans. Autom. Control; 2009; 54, pp. 1254-1269. [DOI: https://dx.doi.org/10.1109/TAC.2009.2019800]

39. Rawlings, J.B.; Bakshi, B.R. Particle filtering and moving horizon estimation. Comput. Chem. Eng.; 2006; 30, pp. 1529-1541. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2006.05.031]

40. Van Der Merwe, R.; Doucet, A.; De Freitas, N.; Wan, E. The unscented particle filter. Adv. Neural Inf. Process. Syst.; 2000; 13, pp. 1529-1541.

41. Rengaswamy, R.; Narasimhan, S.; Kuppuraj, V. Receding-Horizon Nonlinear Kalman (RNK) Filter for State Estimation. IEEE Trans. Autom. Control; 2013; 58, pp. 2054-2059. [DOI: https://dx.doi.org/10.1109/TAC.2013.2253271]

42. Salau, N.P.G.; Secchi, A.R.; Otávio, J. Five Formulations of Extended Kalman Filter: Which is the best for D-RTO?. Proceedings of the 17th European Symposium on Computer Aided Process Engineering—ESCAPE17; Bucharest, Romania, 27–30 May 2007.

43. Böhler, L.; Ritzberger, D.; Hametner, C.; Jakubek, S. Constrained extended Kalman filter design and application for on-line state estimation of high-order polymer electrolyte membrane fuel cell systems. Int. J. Hydrog. Energy; 2021; 46, pp. 18604-18614. [DOI: https://dx.doi.org/10.1016/j.ijhydene.2021.03.014]

44. Spivey, B.J.; Hedengren, J.D.; Edgar, T.F. Constrained Nonlinear Estimation for Industrial Process Fouling. Ind. Eng. Chem. Res.; 2010; 49, pp. 7824-7831. [DOI: https://dx.doi.org/10.1021/ie9018116]

45. Soares, F.D.R.; Secchi, A.R.; Souza, M.B., Jr. Development of a Nonlinear Model Predictive Control for Stabilization of a Gas-Lift Oil Well. Ind. Eng. Chem. Res.; 2022; 61, pp. 8411-8421. [DOI: https://dx.doi.org/10.1021/acs.iecr.1c04728]

46. Tuveri, A.; Nakama, C.S.; Matias, J.; Holck, H.E.; Jäschke, J.; Imsland, L.; Bar, N. A regularized Moving Horizon Estimator for combined state and parameter estimation in a bioprocess experimental application. Comput. Chem. Eng.; 2023; 172, 108183. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2023.108183]

47. Morari, M.; Garcia, C.E.; Prett, D.M. Model predictive control: Theory and practice. IFAC Proc. Vol.; 1988; 21, pp. 1-12. [DOI: https://dx.doi.org/10.1016/B978-0-08-035735-5.50006-1]

48. Rao, C.V.; Wright, S.J.; Rawlings, J.B. Application of interior-point methods to model predictive control. J. Optim. Theory Appl.; 1998; 99, pp. 723-757.

49. Carey, G.; Finlayson, B.A. Orthogonal collocation on finite elements. Chem. Eng. Sci.; 1975; 30, pp. 587-596. [DOI: https://dx.doi.org/10.1016/0009-2509(75)80031-5]

50. Ritschel, T.K.; Jørgensen, J.B. Nonlinear Model Predictive Control for Disturbance Rejection in Isoenergetic-isochoric Flash Processes. IFAC-PapersOnLine; 2019; 52, pp. 796-801. [DOI: https://dx.doi.org/10.1016/j.ifacol.2019.06.159]

51. Jordanou, J.P.; Osnes, I.; Hernes, S.B.; Camponogara, E.; Antonelo, E.A.; Imsland, L. Nonlinear Model Predictive Control of Electrical Submersible Pumps based on Echo State Networks. Adv. Eng. Inform.; 2022; 52, 101553. [DOI: https://dx.doi.org/10.1016/j.aei.2022.101553]

52. Roman, R.; Nagy, Z.K.; Cristea, M.V.; Agachi, S.P. Dynamic modelling and nonlinear model predictive control of a Fluid Catalytic Cracking Unit. Comput. Chem. Eng.; 2009; 33, pp. 605-617. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2008.08.007]

53. Kirches, C.; Wirsching, L.; Bock, H.; Schlöder, J. Efficient direct multiple shooting for nonlinear model predictive control on long horizons. J. Process Control; 2012; 22, pp. 540-550. [DOI: https://dx.doi.org/10.1016/j.jprocont.2012.01.008]

54. Rodríguez, D.E.S. Controle Preditivo não Linear para Sistemas de Parâmetros Distribuídos. Master’s Thesis; Universidade Federal de São Carlos: São Paulo, Brazil, 2014.

55. Giraldo, S.A.C.; Melo, P.A.; Secchi, A.R. Tuning of Model Predictive Controllers Based on Hybrid Optimization. Processes; 2022; 10, 351. [DOI: https://dx.doi.org/10.3390/pr10020351]

56. Fontes, R.M.; Martins, M.A.F.; Odloak, D. An Automatic Tuning Method for Model Predictive Control Strategies. Ind. Eng. Chem. Res.; 2019; 58, pp. 21602-21613. [DOI: https://dx.doi.org/10.1021/acs.iecr.9b03502]

57. Tran, Q.N.; Özkan, L.; Backx, A. Generalized predictive control tuning by controller matching. J. Process Control; 2015; 25, pp. 1-18. [DOI: https://dx.doi.org/10.1016/j.jprocont.2014.10.002]

58. Haseltine, E.L.; Rawlings, J.B. Critical Evaluation of Extended Kalman Filtering and Moving-Horizon Estimation. Ind. Eng. Chem. Res.; 2005; 44, pp. 2451-2460. [DOI: https://dx.doi.org/10.1021/ie034308l]

59. Salau, N.P.G.; Trierweiler, J.O.; Secchi, A.R. State estimation of chemical engineering systems tending to multiple solutions. Braz. J. Chem. Eng.; 2014; 31, pp. 771-785. [DOI: https://dx.doi.org/10.1590/0104-6632.20140313s00002625]

60. Binder, T.; Blank, L.; Bock, H.G.; Bulirsch, R.; Dahmen, W.; Diehl, M.; Kronseder, T.; Marquardt, W.; Schlöder, J.P.; von Stryk, O. Introduction to Model Based Optimization of Chemical Processes on Moving Horizons. Online Optimization of Large Scale Systems; Grötschel, M.; Krumke, S.O.; Rambau, J. Springer: Berlin/Heidelberg, Germany, 2001; pp. 295-339. [DOI: https://dx.doi.org/10.1007/978-3-662-04331-8_18]

61. Pollard, G.; Sargent, R. Off line computation of optimum controls for a plate distillation column. Automatica; 1970; 6, pp. 59-76. [DOI: https://dx.doi.org/10.1016/0005-1098(70)90075-0]

62. Betts, J.T. Practical Methods for Optimal Control and Estimation Using Nonlinear Programming; SIAM: Philadelphia, PA, USA, 2010.

63. Biegler, L.T. Efficient solution of dynamic optimization and NMPC problems. Nonlinear Model Predictive Control; Springer: Cham, Switzerland, 2000; pp. 219-243.

64. Andersson, J.A.E.; Gillis, J.; Horn, G.; Rawlings, J.B.; Diehl, M. CasADi—A software framework for nonlinear optimization and optimal control. Math. Program. Comput.; 2019; 11, pp. 1-36. [DOI: https://dx.doi.org/10.1007/s12532-018-0139-4]

65. Bock, H.; Plitt, K. A Multiple Shooting Algorithm for Direct Solution of Optimal Control Problems*. IFAC Proc. Vol.; 1984; 17, pp. 1603-1608. [DOI: https://dx.doi.org/10.1016/S1474-6670(17)61205-9]

66. Cuthrell, J.E.; Biegler, L.T. On the optimization of differential-algebraic process systems. AIChE J.; 1987; 33, pp. 1257-1270. [DOI: https://dx.doi.org/10.1002/aic.690330804]

67. Blanchard, E.D.; Sandu, A.; Sandu, C. Parameter Estimation Method Using an Extended Kalman Filter. 2007; Available online: https://www.researchgate.net/publication/254410401_Parameter_Estimation_Method_using_an_Extended_Kalman_Filter (accessed on 10 January 2023).

68. Terejanu, G.A. Extended Kalman Filter Tutorial; University at Buffalo: Buffalo, NY, USA, 2008.

69. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J. et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nat. Methods; 2020; 17, pp. 261-272. [DOI: https://dx.doi.org/10.1038/s41592-019-0686-2] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32015543]

70. Stellato, B.; Banjac, G.; Goulart, P.; Bemporad, A.; Boyd, S. OSQP: An operator splitting solver for quadratic programs. Math. Program. Comput.; 2020; 12, pp. 637-672. [DOI: https://dx.doi.org/10.1007/s12532-020-00179-2]

71. Vusse, J.V. A new model for the stirred tank reactor. Chem. Eng. Sci.; 1962; 17, pp. 507-521. [DOI: https://dx.doi.org/10.1016/0009-2509(62)87002-X]

72. Engell, S.; Klatt, K.U. Nonlinear control of a non-minimum-phase CSTR. Proceedings of the 1993 American Control Conference; San Francisco, CA, USA, 2–4 June 1993; IEEE: New York, NY, USA, 1993; pp. 2941-2945.

73. Seborg, D.E.; Edgar, T.F.; Mellichamp, D.A.; Doyle, F.J., III. Process Dynamics and Control; John Wiley & Sons: New York, NY, USA, 2016.