Robust Model Predictive Controller Using Recurrent Neural Networks for Input–Output Linear Parameter Varying Systems

Abstract

This paper develops a model predictive controller (MPC) for constrained nonlinear MIMO systems subjected to bounded disturbances. A linear parameter varying (LPV) model assists MPC in dealing with nonlinear dynamics. In this study, the nonlinear process is represented by an LPV using past input–output information (LPV-IO). Two primary objectives of this study are to reduce online computational load compared with the existing literature of MPC with an LPV-IO model and to confirm the robustness of the controller in the presence of disturbance. For the first goal, a recurrent neural network (RNN) is employed to solve real-time optimization problems with lower online computation. Regarding robustness, a new control law is developed, which comprises a fixed control gain ($K$) and a free perturbation ($C)$. The proposed method enjoys a shrunken conservatism owing to the finding of a larger possible terminal region and using free control moves. The strategy is examined in an alkylation of benzene process and displays outstanding performance in both setpoint tracking and disturbance rejection problems. Moreover, the superiority of RNN over three conventional optimization algorithms is underlined in terms of MSE, the average time for solving the optimization problem, and the value of the cost function.

1. Introduction

Over the past few decades, model predictive controllers (MPCs) have attracted a great deal of attention in industry and academia [1,2,3,4]. This considerable prominence of MPC is rooted in several factors, including simple generalization to a MIMO system, applicability to non-minimum phase and unstable processes, capability to compensate delay, and imposing output/input constraints [1]. While influential theories exist in MPC for linear systems, some systems are nonlinear and operate under large operating conditions, or some external parameters can fundamentally change the process response [5]. Nonlinear MPC (NMPC) has been developed to deal with systems with nonlinear dynamics. A higher level of algorithm complexity and excessive computational load are two main drawbacks of NMPC.

A linear parameter varying (LPV) framework is established in this paper to mitigate these problems by assisting the MPC to model dynamic and static nonlinearities [6,7,8,9]. MPC founded on LPV has aroused enthusiasm among researchers to ensure practical implementation to bridge the gap between linear MPC and NMPC [10,11,12] (less complicated than nonlinear models and more accurate than linear models). With the help of LPV, MPC can be employed for nonlinear and/or time-varying systems [6,13,14,15]. LPV offers MPC a linear model instead of nonlinear models, where the scheduling variable denotes nonlinear dynamic parameters, changing environmental conditions, and operating points. It is also noteworthy that, in an LPV model, current values of scheduling variables can be measured or estimated, but the future values are unknown, so a robust MPC is required, or the future values can be estimated.

Most MPC-LPV literature is centered around the state-space model, while states might not be measurable in many industrial applications. Moreover, using an observer can complicate the design. In Reference [16], MPC-LPV with SS and IO presentations are compared for a control moment gyroscope (CMG). The results showed that the execution time would be reduced when using the IO framework. To this end, MPC is established based on an LPV-IO framework in this paper.

Only a few MPCs have been documented for LPV systems under the I/O formulation, mainly focused on stability, computation, and conservatism [5]. One of the initial research works on MPC with the LPV-IO model can be found at [10], where min-max (worse-case) cost functions are regarded. The cost function includes a stage cost for control performance and a terminal cost to grant stability. The optimization problem of MPC is written in bilinear matrix inequality (BMI), which can be non-convex and challenging. This approach was improved in [17] in order to reduce the complexity level markedly. The authors designed a controller based on a linear matrix inequality (LMI) convex problem, successfully assessed in an ideal continuous stirred tank reactor (CSTR). They also enlarge the terminal region to moderately reduce conservatism. This work employed offline controllers and terminal regions to reduce the online computational burden moderately. However, the controller still suffers from substantial online computational capacity, stemming from solving LMI equations and calculating the Markov coefficient of the LPV system [18]. Above all, it is assumed that the future values of scheduling variables can vary inside a convex polytope. This assumption is less conservative than a frozen or fixed scheduling variable over the prediction horizon while making the optimization problem nonlinear or unreachable [5,17]. The online computational complexity of [17] has been slightly ameliorated in a sequence of quadratic programs. The main drawback of this work is that it lacks a provision for recursive feasibility.

The first novelty of our paper is to make the controller, developed in [10], more computationally efficient. To the best of our knowledge, for the first time, a computationally efficient MPC with an LPV-IO model was presented in this study while keeping the conservatism degree as low as reasonably achievable. A recurrent neural network (RNN) is proposed to deal with the online optimization problem in this paper. Thanks to the global convergence and low computational burden, neural network-based optimizations have been applied to linear, quadratic [19,20], nonlinear programming [21,22], and variational inequality problems [23,24]. The RNN is pointed out to display convincing performance in either MPC or NMPC [25,26,27], in which nonlinear optimization can be converted to quadratic programming. The motive behind selecting RNN is its good adaptability as well as lower complexity, even in real-time nonlinear and large-scale optimization. That is why the studied method was assessed in the alkylation of benzene process, a large-scale nonlinear process with 25 state variables, 5 inputs, and 5 outputs.

The second novelty of this paper is to guarantee the robustness of the system when facing bounded disturbance. The problem of disturbance rejection has not been studied so far for MPC with LPV-IO models. A stabilizing control law based on the input–output measurements is proved to ensure the superior closed-loop performance of the controller in the presence of additive disturbance. The control law includes a fixed control gain ($K$) and a free perturbation ($C$), the former of which is designed offline with the largest possible terminal region to reduce conservatism. The free control moves are determined online to ensure input-state stability in the presence of bounded disturbances. Adding control moves provides a less conservative framework than those studied in [17,18], where only offline control gain is employed.

To sum up, the contributions of this paper are as follows:

• An RNN-based optimization algorithm is developed to offer global convergence and lower the online computational load.
• Free control moves are added to the constant control gain to maintain the closed-loop stability when facing bounded disturbances.
• Concerning previous studies for MPC with LPV, the proposed method inherently enjoys a shrunken conservatism degree as a result of finding the larger possible terminal region, using free control moves, and the global solution of the optimization problem.

The rest of the paper is organized as follows. The standard nonlinear model is reformulated as an LPV-IO model in Section 2. The proposed RMPC-LPV structure is developed in Section 3 and proved to be stable and feasible. The MPC-RNN is presented in Section 4. Next, the state-space equations, inputs, outputs, and critical parameters of the alkylation of benzene process are wholly outlined in Section 5, followed by a collection of simulations to evaluate the proposed control methodology. Key findings are finally summarized in Section 6.

2. Problem Statement

Assume the following discrete-time MIMO linear parameter-varying (LPV) transfer function with constraints:

$$y\left(k\right)=-\left({\displaystyle \sum }_{i=1}^{n_a}{A}_i\left(p\left(k\right)\right)q^{-i}\right)y\left(k\right)-\left({\displaystyle \sum }_{j=1}^{n_b}{B}_j\left(p\left(k\right)\right)q^{-j}\right)u\left(k\right)+d\left(p\left(k\right)\right)$$

(1)

Subject to the following constraints

$$\begin{gathered} u\left(k\right)\in U\equiv \left\{u\in {ℝ}^{n_u}| \left|u\left(k\right)\right|\le u_{max}\right\} \\ \Delta u\left(k\right)\in V\equiv \left\{\Delta u\in {ℝ}^{n_u}| \left|\Delta u\left(k\right)\right|\le \Delta u_{max}\right\} \\ y\left(k\right)\in Y\equiv \left\{y\in {ℝ}^{n_y}| \left|y\left(k\right)\right|\le y_{max}\right\} \end{gathered}$$

(2)

where $y\left(k\right)$ are outputs; $u\left(k\right)$ are inputs; $p\left(k\right)$ are scheduling variables; $q^{-i}$ is backshift operator; $n_a$ is the degree of the output polynomial $;n_b$ is the degree of the input polynomial;$n_u$ is the number of input; $n_y$ is the number of outputs; and $u_{max}, \Delta u_{max}$, and $y_{max}$ are boundaries. The state-space representation of the dynamic model given by Equation (1) is as follows [10]:

$$x\left(k+1\right)=\underset{A\left(p\left(k\right)\right)}{\underbrace{\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{ccc}-a_1& \dots & -a_{n_a-1}\end{array}\\ \begin{array}{ccc}{I}_{n_y}& \dots & 0\end{array}\\ \begin{array}{ccc}\begin{array}{c}⋮\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}& \begin{array}{c}\ddots \\ \dots \\ \begin{array}{c}\dots \\ \dots \\ \begin{array}{c}\ddots \\ \dots \end{array}\end{array}\end{array}& \begin{array}{c}⋮\\ I_{n_y}\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}-\begin{array}{cc}{a}_{n_a}& b_0+b_1\end{array}\\ \begin{array}{cc}0& 0\end{array}\\ \begin{array}{cc}\begin{array}{c}⋮\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}& \begin{array}{c}⋮\\ 0\\ \begin{array}{c}{I}_{n_u}\\ I_{n_u}\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}& \begin{array}{c}\begin{array}{ccc}\dots & b_{n_b-1}& b_{n_b}\end{array}\\ \begin{array}{ccc}\dots & 0& 0\end{array}\\ \begin{array}{ccc}\begin{array}{c}\ddots \\ \dots \\ \begin{array}{c}\dots \\ \dots \\ \begin{array}{c}\ddots \\ 0\end{array}\end{array}\end{array}& \begin{array}{c}⋮\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}⋮\\ I_{n_u}\end{array}\end{array}\end{array}& \begin{array}{c}⋮\\ 0\\ \begin{array}{c}0\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}}x\left(k\right)+\underset{B\left(p\left(k\right)\right)}{\underbrace{\left[\begin{array}{c}\begin{array}{c}\begin{array}{c}{b}_0\\ 0\end{array}\\ ⋮\\ 0\end{array}\\ \begin{array}{c}{I}_{n_y}\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\end{array}\right]}}u\left(k\right)+d\left(k\right)$$

(3)

where $x\left(k\right)={\left[y\left(k-1\right)\dots y\left(k-n_a\right) u\left(k-1\right)\dots u\left(k-n_b\right)\right]}^T$. For the sake of simplicity, $p\left(k\right)$ is removed from all coefficients $a_i$ and $b_i$. The following assumptions are made:

A1: The parameter-varying matrix is $\left[A\left(k\right) B\left(k\right)\right]\in \mathsf{\Omega }=Co\left[\left[A_1 B_1\right],\cdots ,\left[A_m B_m\right]\right]$, where $\mathsf{\Omega }$ is the polytope,$Co$ is the convex hull, and $\left[A_i B_i\right]$ are vertices corresponding to scheduling variable $p_i$:

$$p=Co\left\{p_1,\dots p_m\right\}$$

(4)

A2: The profile of changes in scheduling variables is predetermined from a safety and economic point of view with the following upper and lower bound:

$$p\in P\equiv \left\{p\in {ℝ}^{n_p}| \left|p\left(k\right)\right|\le p_{max}\right\}$$

(5)

A3: The disturbance is bounded:

$$d\left(k\right)\in D\equiv \left\{d\in {ℝ}^{n_d}| \left|d\left(k\right)\right|\le d_{max}\right\}$$

(6)

A4: System 3 is controllable, i.e., the controllability matrix (M) has a full rank:

$$M=\left[B AB A^{2}B\dots A^{n-1}B\right]$$

(7)

3. Robust Model Predictive Controller

This paper aims to define a control law that ensures the input-to-state practical stability (ISpS) of the system (1) with constraints (2). The control law is calculated by minimizing the cost function at any time instance k:

$$\begin{gathered} V_N\left(x_k,u_k,p_k\right)={\displaystyle \sum }_{k=0}^{N-1}[e{\left(k+i\right)}^{T}Qe\left(k+i\right)+\Delta u{\left(k+i-1\right)}^{T}R\Delta u\left(k+i-1\right)] \\ +[x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right)] \\ \mathrm{Subject} \mathrm{to} \\ y\left(k+j\right)=-\left({{\displaystyle \sum }}_{i=1}^{n_a}{A}_i{q}^{-i}\right)y\left(k+j\right)-\left({{\displaystyle \sum }}_{j=1}^{n_b}{B}_i{q}^{-i}\right)u\left(k+j\right)+d\left(k\right), \mathrm{j}=1,2,\dots ,\mathrm{N}-1 \\ u\left(k+j\right)\in U \\ \Delta u\left(k+j\right)\in V \\ y\left(k+j\right)\in Y \\ x\left(k+N\right)\in X_f \end{gathered}$$

(8)

where $e=r-y$, the deviation of output from the reference trajectory, and $\Delta u$ is the control increment. Q and R are positive definite state and input weighting matrices in charge of closed-loop achievements, $Q_p=Q_p{}^T>0$ is the terminal penalty matrix designed to ensure stability, and $X_f$ is the terminal region. The first and second part of $V_N$ are named stage cost ($V_s$) and terminal cost ($V_T$), respectively. Disturbances are supposed to be fixed over the prediction horizon.

According to [28], a state feedback controller can represent the system’s stability with the LPV-IO form. The control law ($u\left(k\right)$) is comprised of a fixed state feedback $K$ and a free control move ($c$), in which $K$ is computed offline and $c$ is obtained from the minimization of (8):

$$u\left(k\right)=-Kx\left(k\right)+c\left(k\right)$$

(9)

The state-feedback gain accounts for maintaining the final state variable $x\left(k+N\right)$ in the terminal region while meeting constraints and keeping the controller as less conservative as possible. In fact, a constant state feedback controller confirms the stability when there is no disturbance. Otherwise, free control moves (c) will be determined online in a min-max problem to guarantee the ISpS. Thus, a part of the controller will be calculated offline for conditions in which there is not any disturbance, and a part of the controller will be set online to deal with existing disturbances or violation of $x\left(k+N\right)\in X_f$.

3.1. Offline Controller

In this section, the stability of the system will be represented when $x\left(k\right)\in X_f$ through finding control gain $K$. After calculating the offline controller, a procedure to specify the terminal region is defined.

Theorem 1.

Taking a disturbance-free state-space model of Equation (3), there exists a control law $u=-Kx$that asymptotically stabilizes the system if$V_T$ is a positive definite Lyapunov function such that

• $V_T\left(x\left(k+1\right)\right)-V_T\left(x\left(k\right)\right)<0$ for $\forall x\in X_f$ And $\forall p\in P$.
• If $x\left(k\right)\in X_f$ then $x\left(k+1\right)\in X_f$
• $\left|u\right|=\left|Kx\right|\le u_{max} for \forall x\in X_f$.

Proof. 

Considering $V_T=x{\left(k\right)}^T{Q}_{p}x\left(k\right)>0$ as the candidate Lyapunov function, the controller $u=-Kx$ exists, if the equation below satisfies the following:

$$\begin{gathered} V_T\left(x\left(k+1\right)\right)-V_T\left(x\left(k\right)\right)<0 \\ x{\left(k+1\right)}^T{Q}_{p}x\left(k+1\right)-x{\left(k\right)}^T{Q}_{p}x\left(k\right)\le -x{\left(k\right)}^{T}Qx\left(k\right)-\Delta u{\left(k\right)}^{T}R\Delta u\left(k\right) \end{gathered}$$

(10)

It is assumed that there exists a K such that $A-BK$ is stable for all possible pairs of $\left[A\left(k\right) B\left(k\right)\right]\in \mathsf{\Omega }$. By substituting $\Delta u=-Kx$ and $x\left(k+1\right)=\left(A-BK\right)x\left(k\right)$, it yields

$$(\left(A-BK\right)x\left(k\right)){ }^T{Q}_p\left(\left(A-BK\right)x\left(k\right)\right)-x{\left(k\right)}^T{Q}_{p}x\left(k\right)\le -x{\left(k\right)}^{T}Qx\left(k\right)-{\left(-Kx\left(k\right)\right)}^{T}R\left(-Kx\left(k\right)\right)$$

(11)

It can be rewritten in a more compact form:

$$\begin{gathered} M^T{Q}_{p}M-Q_p+Q+K^{T}RK\le 0 \\ M=A-BK \end{gathered}$$

(12)

Therefore, $V_T$ is a Lyapunov function, and the stabilizing state-feedback gain K can be found from (12). After showing the controller’s stability, the terminal region is required to be specified, reducing the conservatism level. The terminal region $X_f$ is considered to be an ellipsoidal invariant set:

$$X_f=\left\{x\in R^n|x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right)<\theta \right\}$$

(13)

While $X_f$ is intended to be as broad as possible to reduce conservatism, it might lead to extensive control law and violation of input constraints. The maximum value of $\gamma$ can be derived from an optimization problem such that input limitations are met:

$$\begin{gathered} \max \theta \\ \mathrm{Subject} \mathrm{to} \\ x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right)<\theta \\ \left|Kx\right|\le u_{max} \end{gathered}$$

(14)

In Reference [10], it is proved that $\gamma$ can equivalently be derived from the following optimization problem:

$$\begin{gathered} \underset{\gamma }{\min }\theta \\ {\theta }^2{A}_T{P}^{-1}{A}_T{}^T\le B_T{B}_T{}^T \end{gathered}$$

(15)

where $A_T={\left[-I_n K I_n-K\right]}^T, B_T={\left[x_{max}-x_s \Delta u_{max} x_{max}-x_s \Delta u_{max}\right]}^T$, and $x_s$ is the desired state. In this section, an offline state-feedback controller ($K$) and the corresponding terminal region ($X_f$) that stabilize the system are introduced. Unlike [17], where the future profile of scheduling variables $\left\{p\left(k+1\right),\dots ,p\left(k+N\right)\right\}$ is assumed to be unknown, in this paper, a least-square algorithm was used to find scheduling variables over the prediction horizon. This means that the controller must not be robust against the possible uncertainties in the scheduling variables. For further information about the procedure for predicting the scheduling variables, please refer to [29]. □

3.2. Online Controller

The free control moves will be determined online in the case wherein there is a disturbance or the condition $x\left(k\right)\in X_f$ is not meet. The online controller steers the states towards $X_f$, where offline control asymptotically stabilizes the system. Two theorems are defined here to demonstrate the controller validity. Theorem 2 shows that the system remains ISpS when there is a bounded disturbance. Subsequently, the recursive feasibility of the system will be proved in Theorem 3. Because of additive disturbances, a min-max (worst-case) optimization problem is defined to cope with uncertainties. The optimization problem (8) can be rewritten considering the bounded disturbance:

$$\begin{gathered} \underset{c\left(k\right)}{\min }\underset{d\left(k\right)\in D}{\max }{V}_N\left(x\left(k\right),u\left(k\right),p\left(k\right),K,c\left(k\right)\right) \\ y\left(k+j\right)=-\left({{\displaystyle \sum }}_{i=1}^{n_a}{A}_i{q}^{-i}\right)y\left(k+j\right)-\left({{\displaystyle \sum }}_{j=1}^{n_b}{B}_i{q}^{-i}\right)u\left(k+j\right)+d\left(k\right), \mathrm{j}=1,2,\dots ,\mathrm{N}-1 \\ u\left(k+j\right)=-Kx\left(k+j\right)+c\left(k+j\right) \\ u\left(k+j\right)\in U \\ \Delta u\left(k+j\right)\in V \\ y\left(k+j\right)\in Y \\ d\left(k+j\right)\in D \\ x\left(k+N\right)\in X_f \end{gathered}$$

(16)

The optimal solution of problem (8), subjected to the system (1) with constraints (3), is the sequence $\left[c^{*}\left(k\right) c^{*}\left(k+1\right)\dots c^{*}\left(k+N-1\right)\right]$ corresponding to $\left[u^{*}\left(k\right) u^{*}\left(k+1\right)\dots u^{*}\left(k+N-1\right)\right]$. Two theorems are defined in this section; Theorem 2 for showing stability and Theorem 3 for showing recursive feasibility. To begin with, a robust positively invariant (RPI) set is needed to be defined as bellow:

A set $X$ is RPI for $x^{+}=Ax+Bu+d$ with $u=-Kx$ if $\in X$ $for \forall x^{+}\in X$,$d\in D$, $\left[A\left(k\right) B\left(k\right)\right]\in \mathsf{\Omega }$

Theorem 2.

There exists asequence of optimal control input ($c^{*}$) and$\alpha ,\beta ,\gamma$$\lambda$that ensures the ISpS of the system (1) with bounded disturbances (6) with the following assumptions:

• $\left|u\right|=\left|Kx\right|\le u_{max} for \forall x\in X_f$.
• $X_f$ is an RPI set of the system (1) with $u=-Kx$.
• $\alpha x^{\lambda }\le Px\le \beta x^{\lambda } for \forall x\in X_f.$
• $P{x}^{+}-Px\le -\left(Qx+Rkx\right)+\mu d$ for $x\in X$,$d\in D$, $\left[A\left(k\right) B\left(k\right)\right]\in \mathsf{\Omega }$
• $Qx+Ru\ge \gamma x^{\lambda } for \forall x\in X_f$

The underlying causes of these assumptions can be found in [30].

Proof. 

The definition of ISpS for the system (1) is

$$\begin{gathered} {\alpha }_1\left(\right|\left|x\right|\left|\right)\le V\left(x\right)\le {\alpha }_2\left(\right|\left|x\right|\left|\right)+s_1 \\ V\left(x^{+}\right)-V\left(x\right)\le -{\alpha }_3\left(\right|\left|x\right|\left|\right)+{\alpha }_4\left(\right|\left|d\right|\left|\right)+s_2 \end{gathered}$$

(17)

where $x^{+}$ RPI set in which $x^{+}=Ax+Bu+d$ with $u=-Kx$; ${\alpha }_1,{\alpha }_2,$ and ${\alpha }_3$ are ${\mathcal{K}}_{\infty }$-function; ${\alpha }_4$ is $\mathcal{K}$-function; and $s_1$ and $s_2$ are positive numbers such that $x\in X,$ $x^{+}\in X$, and $d\left(k\right)\in D$. It is also noteworthy that ISpS is equivalent to ISS for $s_1=s_2=0$. The $V_N\left(x\right)>0$ should first be proved to be bounded so that Equation (17) can be expressed as follows:

$$V_i\left(x\left(k+i\right)\right)=\underset{c_k}{\min }\underset{d\left(k\right)\in D}{\max }\left[x{\left(k+i\right)}^{T}Qx\left(k+i\right)+u{\left(k+i\right)}^{T}Ru\left(k+i\right)+V_{i+1}\left(x\left(k+i+1\right)\right)\right]$$

(18)

Equation (18) is derived from a mathematical induction technique such that the following holds:

$$\begin{gathered} V_0\left(x\left(k\right)\right)=\underset{c_k}{\min }\underset{d\left(k\right)\in D}{\max }\left[x{\left(k\right)}^{T}Qx\left(k\right)+u{\left(k\right)}^{T}Ru\left(k\right)+V_1\left(x\left(k+1\right)\right)\right] \\ {V}_1\left(x\left(k+1\right)\right)=\underset{c_k}{\min }\underset{d\left(k\right)\in D}{\max }\left[x{\left(k+1\right)}^{T}Qx\left(k+1\right)+u{\left(k+1\right)}^{T}Ru\left(k+1\right)+V_2\left(x\left(k+2\right)\right)\right] \\ \dots \\ {V}_N\left(x\left(k+N\right)\right)=[x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right) \end{gathered}$$

(19)

Equation (10), when the disturbance is regarded to be non-zero, results in the following:

$$x{\left(k+1\right)}^T{Q}_{p}x\left(k+1\right)-x{\left(k\right)}^T{Q}_{p}x\left(k\right)\le -x{\left(k\right)}^{T}Qx\left(k\right)-u{\left(k\right)}^{T}Ru\left(k\right)+d_{max}$$

(20)

By developing this equation for the next value of state variables, the following equations will be reached:

$$\begin{gathered} x{\left(k+2\right)}^T{Q}_{p}x\left(k+2\right)-x{\left(k\right)}^T{Q}_{p}x\left(k\right)\le -x{\left(k+1\right)}^{T}Qx\left(k+1\right)-u{\left(k+1\right)}^{T}Ru\left(k+1\right)+d_{max} \\ \dots \\ x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right)-x{\left(k\right)}^T{Q}_{p}x\left(k\right)\le -{\displaystyle \sum }_{i=0}^{N-1}x{\left(k+i\right)}^{T}Qx\left(k+i\right)+u{\left(k+i\right)}^{T}Ru\left(k+i\right)+d_{max} \end{gathered}$$

(21)

From Equations (18) and (20), and substituting $u\left(k\right)=-Kx\left(k\right)$, the following is obtained.

$$V_N\left(x\left(k\right)\right)\le \underset{d\left(k\right)\in D}{\max }\left[x{\left(k\right)}^{T}Qx\left(k\right)-x{\left(k\right)}^{T}KRx\left(k\right)\right]+x{\left(k+1\right)}^T{Q}_{p}x\left(k+1\right)$$

(22)

where $V\left(x_{k+1}\right)\le x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right)$ and, therefore,

$$V_N\left(x\right)\le x{\left(k\right)}^T\left(Q-KR+Q_p\right)x\left(k\right)+d_{max}=\left(Q-KR+Q_p\right)\left|\right|x\left|\right|+d_{max}$$

(23)

It can be concluded that $V_N\left(x\right)$ has an upper bound and ${\alpha }_1=Q-KR+Q_p$, $s_1=d_{max}$ in the right-hand side of the first equation of (23). The difference between two sequential Lyapunov functions is

$$V\left(x\left(k+1\right)\right)-V\left(x\left(k\right)\right)\le \underset{d\left(k\right)\in D}{\max }-\left[x{\left(k\right)}^{T}Qx\left(k\right)+u{\left(k\right)}^{T}Ru\left(k\right)\right]+x{\left(k+N\right)}^T{Q}_{p}x\left(k+N\right)$$

(24)

This can be expressed by

$$V\left(x\left(k+1\right)\right)-V\left(x\left(k\right)\right)\le x{\left(k\right)}^T\left(-Q+KR+Q_p\right)x\left(k\right)+d_{max}=-\left(Q-KR-Q_p\right)x+d_{max}$$

(25)

Accordingly, the difference of Lyapunov function is bounded such that ${\alpha }_3=Q-KR-Q_p$, $s_2=d_{max}$. It has been proven that the proposed controller is ISpS, and the optimal solution, satisfying input, output, and terminal constraints, can be derived from a min-max problem (16). In the next stage, it has been shown that the solution is feasible. □

Theorem 3.

Given system (1) with constraints (2) and bounded disturbance$\left(6\right)$, the closed-loop system is ISpS for$x_0\in x_N$, and the system is feasible.

Proof. 

Suppose that there exists a terminal region, derived from (15), and a control gain K, calculated by Equation (12); then, the system is stable and $x_0\in x_N$. Consequently, if the optimization problem (8) is feasible at time instance k, it remains feasible at the subsequent instances. After finding the optimal K, $c^{*}$, and the corresponding $u^{*}$ at time instance k, the system of Equation (3) at the next instance can be rewritten.

$$x\left(k+1\right)=A\left(p\left(k\right)\right)x_k+B\left(p\left(k\right)\right)u_k{}^{*}+d\left(k\right)$$

(26)

Moreover,

$$x\left(k+2\right)=A\left(p\left(k+1\right)\right)x_1+B\left(p\left(k+1\right)\right)u_0{}^{*}+d\left(k+1\right)$$

(27)

In general, the sequence of system states can be described as follows:

$$x\left(k+i+1\right)=A\left(p\left(k+i\right)\right)x_{k+i}+B\left(p\left(k+i\right)\right)u_{k+i}{}^{*}+d\left(k+i\right)$$

(28)

By substituting $u\left(k\right)=-Kx\left(k\right)+c\left(k\right)$, Equation (28) can be derived:

$$x\left(k+i+1\right)=\left(A\left(p\left(k+i\right)\right)-B\left(p\left(k+i\right)\right)K\right)*x_{k+i}+B\left(p\left(k+i\right)\right)c_{k+i}{}^{*}+d\left(k+i\right)$$

(29)

The above equation for $i=N$ results in

$$x\left(k+N+1\right)=\left(A\left(p\left(k+N\right)\right)-B\left(p\left(k+N\right)\right)\right)*x_{k+N}+B\left(p\left(k+N\right)\right)c_{k+N}{}^{*}+d\left(k+N\right)$$

(30)

This implies that the optimal solution is feasible because

$$\begin{gathered} for i=0, X_{N-1}\in X_N \\ for i=1, X_{N-2}\in X_{N-1}\in X_N \\ for i=N, X_0\in X_N and x\left(k+N+1\right)\in X_N \end{gathered}$$

(31)

According to Theorems 1 to 3, the developed controller has two components. The first is determined offline, which included a control gain $K$ and terminal region $X_f$. When $x\in X_f$, the offline controller asymptotically stabilizes the system. Otherwise, when there is uncertainty, an online controller verifies the ISpS. The online optimization problem is shown to be feasible. In the next section, a dynamic neural network is constructed to solve the online optimization problem for the first time. □

4. Real-Time Optimization Problem Using RNN

The global convergence and low complexity of RNN for optimization of linear models, constrained linear models, linear models with uncertainty, and various nonlinear models are represented in the literature. In this study, RNN optimizes a real-time QP problem enjoying parallel computation. To begin with, the original optimization problem is required to transform into a standard form. The input-state relationship of Equation (1) can be expressed in the following form:

$$x\left(k+1\right)=A\left(p\left(k\right)\right)x\left(k\right)+B\left(p\left(k\right)\right)u\left(k\right)+d\left(p\left(k\right)\right)$$

(32)

The vectors of predicted outputs, inputs, and disturbance are as follows:

$$\begin{gathered} X\left(k\right)={\left[X\left(k+1\right)\dots X\left(k+N\right)\right]}^T\in {ℛ}^{m\times N} \\ U\left(k\right)={\left[u\left(k\right)\dots u\left(k+N-1\right)\right]}^T\in {ℛ}^{n\times N} \\ D\left(k\right)={\left[d\left(k+1\right)\dots d\left(k+N\right)\right]}^T\in {ℛ}^{o\times N} \end{gathered}$$

(33)

where $m, n,$ and $o$ are the number of states, inputs, and disturbances, respectively. The predicted states can be easily shown to be in the following form:

$$\begin{gathered} X\left(k+j\right)=G\left(p\left(k+j-1\right)\right)X\left(k+j-1\right)+F\left(p\left(k+j-1\right)\right)U\left(k+j-1\right)+D\left(k+j-1\right),j=1,\dots ,N \\ G\left(p\left(k\right)\right)={\left[A\left(p\left(k\right)\right) A{\left(p\left(k\right)\right)}^2\dots A{\left(p\left(k\right)\right)}^N\right]}^T \\ F\left(p\left(k\right)\right)=\left[\begin{array}{c}\begin{array}{c}B\left(p\left(k\right)\right)\\ A\left(p\left(k\right)\right)B\left(p\left(k\right)\right)+B\left(p\left(k\right)\right)\end{array}\\ \begin{array}{c}⋮\\ A{\left(p\left(k\right)\right)}^{N-1}B\left(p\left(k\right)\right)+\dots +A\left(p\left(k\right)\right)B\left(p\left(k\right)\right)+B\left(p\left(k\right)\right)\end{array}\end{array}\right] \end{gathered}$$

(34)

The constraints (2) can be represented as follows:

$$\begin{gathered} -u_{max}\le u\left(k\right) \le u_{max} \\ -\Delta u_{max}\le u\left(k\right)-u\left(k-1\right) \le \Delta u_{max} \\ -x_{max}\le Gu+F+d \le x_{max} \end{gathered}$$

(35)

The optimization problem (16) can be written as follows:

$$\begin{gathered} \underset{U}{\min }{X}^T\overline{Q}X+U^{T}RU=\underset{U}{\min }{\left(GX+FU+D\right)}^T\overline{Q}\left(GX+FU+D\right)+U^{T}RU \\ \underset{U}{\min }{U}^T\times ({\left(FU\right)}^T\overline{Q}\left(FU\right)+R)\times U+{\left(GX+FU\right)}^T\overline{Q}\left(GX+FU\right) \end{gathered}$$

(36)

where

$$\overline{Q}=\left[\begin{array}{cc}Q& 0\\ 0& Q_p\end{array}\right]$$

According to Equations (35) and (36), the standard form can be expressed as follows:

$$\begin{gathered} \underset{v}{\min }\frac{1}{2}{v}^{T}Hv+b^{T}v \\ \mathrm{Subject} \mathrm{to} Tv\le q \end{gathered}$$

(37)

where

$$\begin{gathered} H=2\times ({\left(FU\right)}^T\overline{Q}\left(FU\right)+R) \\ b=2\times {\left(GX+FU\right)}^T\overline{Q}\left(GX+FU\right) \\ T=\left[\begin{array}{c}\begin{array}{c}{I}_{n\times n}\\ -I_{n\times n}\end{array}\\ I_{n\times n}\\ \begin{array}{c}-I_{n\times n}\\ \begin{array}{c}G\\ -G\end{array}\end{array}\end{array}\right], q=\left[\begin{array}{c}\begin{array}{c}{u}_{max}\\ -u_{max}\end{array}\\ \Delta u_{max}+u\left(k-1\right)\\ \begin{array}{c}\Delta u_{max}-u\left(k-1\right)\\ \begin{array}{c}{X}_{max}-F-d\\ -X_{max}+F+d\end{array}\end{array}\end{array}\right] \end{gathered}$$

Remark 1.

The vector$q$is unknown as$d$is not measurable, while it can be estimated by the difference between the measured output and model output. In Reference [31], a simple model to find unknown parameters is proposed.

Remark 2.

In References [25,31], a simplified form of a dual neural network is described to ensure the lower computational burden by defining the dual form of the optimization problem (37) as follows:

State equation:

$$\frac{d\omega }{dt}=\lambda \left(-T{H}^{-1}{T}^T\omega +M\left(T{H}^{-1}{T}^T\omega -T{H}^{-1}c-\omega \right)+T{H}^{-1}c\right)$$

(38)

where $\omega$ is the state variable of the network, $\lambda >0$ adjusts the convergence rate of RNN, and M is a piecewise linear function as follows:

$$M\left(z\right)=\{\begin{array}{c}{q}_{min} z<q_{min}\\ z q_{min}<z<q_{max}\\ q_{max} z>q_{max}\end{array}$$

(39)

Output equation:

$$v=H^{-1}{T}^T\omega -H^{-1}c$$

(40)

In the dual form of the optimization problem, constraints are added to the primary cost function as a penalty term. This means that, if $Tv<q,$ the cost function decreases by a factor of $\alpha$; otherwise, for $Tv>q$, the cost function is penalized. The single-layer RNN is given in Figure 1. The global convergence of the proposed network is verified in [31].

Finally, the proposed MPC algorithm can be summarized in the following steps:

• The value of $Q, R,\lambda ,N,u_{max},y_{max},d_{max},p$ and model $\left(A,B\right)$;
• Specify K using (12) and $\mathsf{\theta }$ using (15);
• Repeat the procedure of finding $u\left(k\right)$ by solving Equations (38)–(40).

The whole procedure is also depicted in Figure 2.

5. Case Study

Alkylation of benzene with ethylene, a principal process in the petrochemical industry, produces ethylbenzene, which is widely used as a large-scale benchmark. As shown in Figure 3, the process studied in this work consists of four continuously stirred tank reactors (CSTRs) and a flash tank separator. A detailed description of the process is described in [32,33]. The process’s manipulated variables are heat inputs to five vessels shown by Q1, Q2, Q3, Q4, and Q5. The temperatures of five vessels T1, T2, T3, T4, and T5 are considered the process’s outputs. The process’s dynamic behavior is set out in detail in Equations (41)–(45), and all parameters are named in

Table 1. Process variables of the alkylation of benzene process.

Variables	Definition
$C_{A1},C_{B1},C_{C1},C_{D1}$	Concentrations of A, B, C, D in CSTR-1
$C_{A2},C_{B2},C_{C2},C_{D2}$	Concentrations of A, B, C, D in CSTR-2
$C_{A3},C_{B3},C_{C3},C_{D3}$	Concentrations of A, B, C, D in CSTR-3
$C_{A4},C_{B4},C_{C4},C_{D4}$	Concentrations of A, B, C, D in Separator
$C_{A5},C_{B5},C_{C5},C_{D5}$	Concentrations of A, B, C, D in CSTR-1
$C_{Ar},C_{Br},C_{Cr},C_{Dr}$	Concentrations of A, B, C, D in $F_r,F_{r1},F_{r2}$
$T_1,T_3,T_3,T_4,T_5$	Temperatures in each vessel
$T_{ref}$	Reference temperature
$F_3,F_5,F_7,F_8,F_9$	Effluent flow rates from each vessel
$F_1,F_2,F_4,F_6,F_{10}$	Feed flow rates to each vessel
$F_r,F_{r1},F_{r2}$	Recycle flow rates
$H_{vapA},H_{vapB},H_{vapC},H_{vapD}$	Enthalpies of vaporization of A, B, C, D
$H_{Aref},H_{Bref},H_{Cref},H_{Dref}$	Enthalpies of A, B, C, D at $T_{ref}$
$\Delta H_{r1},\Delta H_{r2},\Delta H_{r3}$	Heat of reactions 1, 2, and 3
$V_1,V_2,V_3,V_4,V_5$	Volume of each vessel
$Q_1,Q_2,Q_3,Q_4,Q_5$	External heat/coolant inputs to each vessel
$C_{pA},C_{pB},C_{pC},C_{pD}$	Heat capacity of A, B, C, D at liquid phase
${\alpha }_A,{\alpha }_B,{\alpha }_C,{\alpha }_D$	Relative volatilities of A, B, C, D
$C_{A0},C_{B0},C_{C0},C_{D0}$	Molar densities of pure A, B, C, D
$T_{A0},T_{B0},T_{D0}$	Feed temperatures of pure A, B, D
$k$	Fraction of overhead flow recycled to the reactors
$C_{A1},C_{B1},C_{C1},C_{D1}$	Concentrations of A, B, C, D in CSTR-1
$C_{A2},C_{B2},C_{C2},C_{D2}$	Concentrations of A, B, C, D in CSTR-2
$C_{A3},C_{B3},C_{C3},C_{D3}$	Concentrations of A, B, C, D in CSTR-3
$C_{A4},C_{B4},C_{C4},C_{D4}$	Concentrations of A, B, C, D in Separator
$C_{A5},C_{B5},C_{C5},C_{D5}$	Concentrations of A, B, C, D in CSTR-1
$C_{Ar},C_{Br},C_{Cr},C_{Dr}$	Concentrations of A, B, C, D in $F_r,F_{r1},F_{r2}$
$T_1,T_3,T_3,T_4,T_5$	Temperatures in each vessel
$T_{ref}$	Reference temperature
$F_3,F_5,F_7,F_8,F_9$	Effluent flow rates from each vessel
$F_1,F_2,F_4,F_6,F_{10}$	Feed flow rates to each vessel
$F_r,F_{r1},F_{r2}$	Recycle flow rates
$H_{vapA},H_{vapB},H_{vapC},H_{vapD}$	Enthalpies of vaporization of A, B, C, D
$H_{Aref},H_{Bref},H_{Cref},H_{Dref}$	Enthalpies of A, B, C, D at $T_{ref}$
$\Delta H_{r1},\Delta H_{r2},\Delta H_{r3}$	Heat of reactions 1, 2, and 3
$V_1,V_2,V_3,V_4,V_5$	Volume of each vessel
$Q_1,Q_2,Q_3,Q_4,Q_5$	External heat/coolant inputs to each vessel
$C_{pA},C_{pB},C_{pC},C_{pD}$	Heat capacity of A, B, C, D at liquid phase
${\alpha }_A,{\alpha }_B,{\alpha }_C,{\alpha }_D$	Relative volatilities of A, B, C, D
$C_{A0},C_{B0},C_{C0},C_{D0}$	Molar densities of pure A, B, C, D

. For more information and the value of fixed parameters, refer to [34].

$$\frac{d{T}_1}{dt}=\frac{Q_1+F_1{C}_{A0}{H}_A\left(T_{A0}\right)+F_2{C}_{B0}{H}_B\left(T_{B0}\right)+{{\displaystyle \sum }}_i^{A,B,C,D}{F}_{r2}{C}_{ir}{H}_i\left(T_4\right)-F_3{C}_{i1}{H}_i\left(T_1\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i1}{C}_{pi}{V}_1}+\frac{-\Delta H_{r1}{r}_1\left(T_1,C_{A1},C_{B1}\right)-\Delta H_{r2}{r}_2\left(T_1,C_{B1},C_{C1}\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i1}{C}_{pi}}$$

(41)

$$\frac{d{T}_2}{dt}=\frac{Q_2+F_4{C}_{B0}{H}_B\left(T_{B0}\right)+{{\displaystyle \sum }}_i^{A,B,C,D}{F}_3{C}_{i1}{H}_i\left(T_1\right)-F_5{C}_{i2}{H}_i\left(T_2\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i2}{C}_{pi}{V}_2}+\frac{-\Delta H_{r1}{r}_1\left(T_2,C_{A2},C_{B2}\right)-\Delta H_{r2}{r}_2\left(T_2,C_{A2},C_{B2}\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i2}{C}_{pi}}$$

(42)

$$\frac{d{T}_3}{dt}=\frac{Q_3+F_6{C}_{B0}{H}_B\left(T_{B0}\right)+{{\displaystyle \sum }}_i^{A,B,C,D}{F}_5{C}_{i2}{H}_i\left(T_2\right)-F_7{C}_{i3}{H}_i\left(T_3\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i3}{C}_{pi}{V}_3}+\frac{-\Delta H_{r3}{r}_1\left(T_3,C_{A3},C_{B3}\right)-\Delta H_{r3}{r}_2\left(T_3,C_{B3},C_{C3}\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i3}{C}_{pi}}$$

(43)

$$\frac{d{T}_4}{dt}=\frac{Q_4+{{\displaystyle \sum }}_i^{A,B,C,D}{F}_7{C}_{i3}{H}_i\left(T_3\right)+F_9{C}_{i5}{H}_i\left(T_5\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i4}{C}_{pi}{V}_4}+\frac{{{\displaystyle \sum }}_i^{A,B,C,D}-M_i{H}_i\left(T_4\right)-F_8{C}_{i4}{H}_i\left(T_4\right)-M_i{H}_{vapi}}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i4}{C}_{pi}{V}_4}$$

(44)

$$\frac{d{T}_5}{dt}=\frac{Q_5+F_{10}{C}_{D0}{H}_D\left(T_{D0}\right)+{{\displaystyle \sum }}_i^{A,B,C,D}{F}_{r1}{C}_{ir}{H}_i\left(T_4\right)-F_9{C}_{i5}{H}_i\left(T_5\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i5}{C}_{pi}{V}_5}+\frac{-\Delta H_{r2}{r}_2\left(T_5,C_{B5},C_{C5}\right)-\Delta H_{r3}{r}_3\left(T_5,C_{A5},C_{D5}\right)}{{{\displaystyle \sum }}_i^{A,B,C,D}{C}_{i5}{C}_{pi}}$$

(45)

The steady-state values of inputs and outputs and initial conditions are indicated in

Table 2. The steady-state and initial values.

Steady-State Temperatures of Vessels (K)	Steady-State Inputs (J/K)	Initial Temperatures of Vessels (K)
$T_{1s}=477.24$	$Q_{1s}=-4.4\times {10}^6$	$T_{1,0}=443.02$
$T_{2s}=476.97$	$Q_{2s}=-4.6\times {10}^6$	$T_{2,0}=437.12$
$T_{3s}=473.47$	$Q_{3s}=-4.7\times {10}^6$	$T_{3,0}=428.37$
$T_{4s}=470.60$	$Q_{4s}=9.2\times {10}^6$	$T_{4,0}=433.15$
$T_{5s}=478.28$	$Q_{5s}=5.9\times {10}^6$	$T_{5,0}=457.55$

. The controller desires to steer the system from this initial condition to the steady- states condition while satisfying the constraints:

$$\begin{gathered} \left|Q_1-Q_{1s}\right|<7.5\times {10}^5, \left|Q_2-Q_{2s}\right|<5\times {10}^5, \left|Q_3-Q_{3s}\right|<5\times {10}^5 \\ \left|Q_4-Q_{4s}\right|<6\times {10}^5, \left|Q_5-Q_{5s}\right|<5\times {10}^5 \end{gathered}$$

In contrast to previous works that developed nonlinear MPC and/or noncentralized MPC, an MPC setup built on the LPV model, improved by RNN, is proposed in this paper, having lower complexity and better performance. More significantly, $F_r$, $C_{C0}$, and $C_{D0}$ are considered to be time-varying, for the first time, to the best of the authors’ knowledge. From a practical point of view, these variables are selected to be scheduling variables with $\pm 10\%$ variation around their nominal values. Before applying the controller, the LPV-IO framework is required to be identified so that ten levels are defined within the range of 90% to 110% of the nominal value for each scheduling variable. The global LPV model originates from a polynomial interpolation of ten local LTI models, identified with 500 data points in each level (operating condition). Each local model has five inputs and five outputs. All manipulated inputs are chosen to be pseudorandom binary signals to extract data with a sampling time of $T=10 s$. With the measured data, the LTI models, afterward, are identified using MATLAB system identification Toolbox. With several simulations, it has been realized that a third-degree polynomial interpolation produced a good fit for constructing the global LPV-IO model, while the order of local transfer functions varies from second to fifth.

The scheduling variables and the responses that come from the global LPV-IO model and nonlinear model are presented in Figure 4 and Figure 5, respectively. The figures supplied show that the identified LPV model can predict the nonlinear system’s behavior. Moreover, mean square error (MSE) and Akaike’s final prediction error (FPE) for the global model were 0.032 and 0.058 in succession, and for the local LPV models, the range of MSE was between 0.004 and 0.062, and FPE changed from 0.040 to 0.081. Figure 5 and the numerical results show that the LPV-IO model is sufficiently reliable for estimating responses.

6. Results and Discussion

In the next step, the performance of the controller in terms of setpoint tracking and disturbance rejection is investigated, in which a step change at $t=1500$ s with amplitude $2\times {10}^{-3}$ is made in $F_1$ and $F_2$. The controller parameters were $N=6, Q=1000\times I_{5\times 5}, R={10}^{-7}\times I_{5\times 5}$, and $\lambda =3.5$. All controllers employed in the performed comparison used the same parameters. The proposed approach is compared with linear RMPC [35], and the LPV controller studied in [17] is shown in Figure 6, where the dashed red line denotes the references, the blue line is the proposed method, the green line is linear RMPC, and the black line is LPV-IO RMPC. In all methods, MPC parameters are chosen according to [36,37].

Looking firstly at reference tracking, the proposed controller’s rise time is slightly larger than other methods, while being fast enough. On the other hand, the proposed controller has a lower settling time. In contrast to LPV-MPCs that reach the reference, the linear RMPC has a non-zero steady-state error, especially in $T_4$ and $T_5$. The proposed controller reached the setpoint without any oscillations, while linear RMPC and LPV-RMPC suffer from large-amplitude oscillations. The reason behind the fluctuation in $T_5$ for linear RMPC can be found in Equation (35), such that $C_{D0}$ directly affects the output, and LRMPC failed to tackle the changes in molar densities of pure D. In general, regardless of rise time, settling time, and overshoot, all three methods followed the setpoints for $T_1$ to $T_4$.

Turning to disturbance rejection, the studied approach outperforms others with removing the disturbance in a short time ranging from 120 s to 355 s for $T_1$ to $T_5$. After time t = 1500 s, when the disturbance was applied to the process, the LRMPC violated both inputs and outputs constraints, especially in $T_1$, although it can overcome the disturbances after a while; on the other hand, LPV-RMPC became highly unstable. Therefore, LPV-RMPC failed to cope with disturbances imposing the process in contrast to other methods, and high amplitude of oscillation in the LRMPC response can be devastating.

Regarding the cost function, the wasted resources by LPV-RMPC are 5.04 times more than LRMPC and 1.59 × 10 times greater than those of the proposed method, which proved that the studied approach is significantly cost-effective. More importantly, the cost function of LRMPC and LPV-RMPC did not converge to zero, so that they cannot deal with changes in three scheduling variables. The MSE of all three methods for different outputs is reported in

Table 3. The steady-state and initial values.

Method	${\mathit{T}}_1$	${\mathit{T}}_2$	${\mathit{T}}_3$	${\mathit{T}}_4$	${\mathit{T}}_5$	Total
LRMPC	391.40	387.25	346.72	331.76	408.73	347.56
LPV-RMPC	$5.33\times {10}^3$	$5.31\times {10}^3$	$5.07\times {10}^3$	$4.90\times {10}^3$	$5.41\times {10}^3$	$5.09\times {10}^3$
Proposed method	67.32	65.30	52.34	60.00	76.45	43.84

. Using the proposed method leads to a sharp decline in error compared with other methods, as the MSE for LRMPC, LPV-RMPC, and proposed method were 447.56, 5.09 × 103, and 43.84, respectively. In short, the proposed method had an acceptable speed (rise time, settling time, and time required to remove disturbance), as well as lower MSE and cost.

After evaluating the performance of the proposed controller, the importance of RNN needs to be emphasized. Finding the global minimum in a shorter time was the primary reason behind using RNN in this study. To this end, RNN is compared with three widely used optimization algorithms for MPCs, namely, sequential quadratic programming (SQP), genetic algorithm (GA), and singular value decomposition (SVD). It has been proven that using RNN reduces the average time required for computing the control action in each sampling in contrast with NMPC based on SQP, GA, and SVD. The average run-time and MSE of these methods is reported in

Table 4. A comparison of different optimization algorithms for finding the control signals.

Optimization Algorithm	MSE	Average Time	Cost
RNN	43.84	0.033	$2.03\times {10}^4$
SQP	153.12	0.76	$8.04\times {10}^6$
GA	76.55	0.81	$5.01\times {10}^4$
SVD	241.22	0.013	$2.05\times {10}^8$

and Figure 7, Figure 8 and Figure 9. The RNN-based method experienced far less MSE (43.84) and cost ($2.03\times {10}^4$) than other optimization algorithms and found the optimal control actions remarkably quick (0.033). GA ranked second in cost and MSE, while being slow, resulting in instability and significant fluctuations in some simulations. In stark contrast, however, SVD has faster responses, although it failed to converge the minimum cost function at times and saw monumental errors. SQP, on the other side, had a smaller MSE and cost than SVD and solved the optimization problem in a shorter time than GA did. According to Figure 7, RNN and SVD have the lowest and highest MSE, respectively. Figure 8 depicts that SQP and GA are sluggish as opposed to RNN and SVD. As described in Figure 9, SVD is the only method with a significant error and fails to find the proper control actions from the optimization problem.

7. Conclusions

An RMPC with an LPV-IO model is investigated in this paper, in which an RNN algorithm solves the real-time optimization problem. The study had three contributions. Firstly, an RNN-based optimization algorithm was developed to offer global convergence and lower the online computational load. Secondly, free control moves were added to the constant control gain to maintain the closed-loop stability when facing bounded disturbances. Lastly, concerning previous studies for MPC with LPV, the proposed method inherently enjoys a shrunken conservatism degree owing to finding the larger possible terminal region, using free control moves, and the global solution of the optimization problem.

This approach’s effectiveness was monitored and compared with LRMPC [35] and LPV-RMPC [17] in an alkylation of benzene process, nonlinear, and large-scale with three scheduling variables. The proposed method was astonishingly successful in both setpoint tracking and disturbance rejection having a reasonable rise time, settling time, MSE, and response amplitude. The LRPMC had a similar speed response while suffering from adverse oscillation. The LPV-RMPC can partially track the predefined reference outputs, whereas it failed to remain stable when facing disturbances. Meanwhile, four optimization algorithms are utilized to solve the online optimization problem of the proposed controller. The results showed that RNN and GA outdid in error reduction and finding the optimal solution; nonetheless, GA was sluggish, and thus detrimental to stability. SVD had the fastest convergence rate with the highest MSE and cost value, and SQP, with a stable response, had the worst execution with an excellent average time and MSE. Future works will focus on showing the stability of the LPV-IO system with multiplicative disturbances. Moreover, the conservatism can be further lowered by ellipsoidal invariant sets with polyhedral invariant sets.