SD-ARX modeling and robust MPC with variable feedback gain for nonlinear systems

Abstract

As a generalized input-output model, the state-dependent exogenous variable autoregressive (SD-ARX) model has been intensively utilized to model complex nonlinear systems. Considering that more freedom can be provided by the state feedback control with variable feedback gain for constructing robust controllers, we propose a robust model predictive control (RMPC) algorithm with variable feedback gain on the basis of the SD-ARX model. First, the polytopic state space models (SSMs) of the system are constructed and the prediction accuracy of the SSMs is further improved by using the parameter variation rate information of the SD-ARX model. Then, an RMPC algorithm with variable feedback gain is synthesized for increasing the design freedom and enhancing the control performance. Two simulation examples, that is, the modeling and control of a continuous stirred tank reactor (CSTR) and a water tank system, are provided to demonstrate the feasibility and effectiveness of the proposed RMPC algorithm.

Introduction

Nonlinear systems with physical constraints can be controlled by a model-based control technique called model predictive control (MPC), which has been attracting widespread attention from both the academic sphere and the industrial domain, respectively. ¹ For the MPC, accurate system model is very important, because usually this kind of algorithm needs to solve an optimal control problem based on model prediction at each sampling interval. However, in actual industrial systems, it is difficult to obtain accurate explicit models due to various uncertainties. ² So far, dealing with systemic uncertainties such as modeling errors and external disturbances is one of the main topics in MPC research. ³

Since the dynamics of nonlinear and/or uncertain systems can be conveniently wrapped by the models in the polytopic form, the polytopic model has attracted much attention. ⁴ Moreover, based on this type of model, robust model predictive control (RMPC) algorithms for the system with disturbances and/or uncertainties can also be conveniently constructed to achieve a reasonable performance. For a system with uncertain parameters, Kothare et al. ⁵ proposed a linear matrix inequality (LMI)-based RMPC algorithm for state feedback control optimization of the controller. For input-constrained linear parameter varying (LPV) systems, Lee et al. ⁶ proposed an RMPC algorithm based on novel sufficient conditions for the cost monotonicity. Lu and Arkun ⁷ proposed two quasi-min-max MPCs for polytopic linear parameter varying (PLPV) systems. For the input-constrained LPV system with unknown but bounded disturbances, He et al. ⁸ proposed an input-to-state stable RMPC algorithm. For the constrained PLPV systems with bounded disturbances, Casavola et al. ⁹ used a novel RMPC for addressing feedback regulation problems. Li et al. ¹⁰ proposed a new RMPC algorithm with constrained feedback for uncertain systems. Polytopic form model-based RMPCs have produced effective results,^5–9 but the majority of the research has aimed to deal with the problems of system’s state or output tracking by assuming the given or accurately measurable steady states of a system. ¹¹ Nevertheless, the steady states of most systems remain unknown or beyond quantification in practicality, owing to extraneous disturbances or inescapable approximative disparities. ¹²

To solve these problems, attention has been focused on system state estimator-based RMPC synthesis schemes. ¹³ For polytopic systems, Ding et al. ¹³ proposed a novel feedback RMPC with a state estimator-based synthesis scheme. For LPV systems, Park et al. ¹² proposed a quasi-min–max RMPC algorithm. Song et al. ¹⁴ proposed a distributed RMPC with output feedback. For uncertain systems, Hu and Ding ¹⁵ proposed a dynamic output feedback RMPC. However, employing observers for the purpose of estimating the states may introduce additional complexity, and when the constraints of the variables such as the state, input, and output are activated, the closed-loop performance may be greatly degraded in terms of input disturbance rejection.^11,16

It is important to directly construct an output feedback RMPC with an input-output type model, since modeling an actual system with an input-output equation is generally easier. ¹⁶ As a generalized input-output model, the state-dependent exogenous variable autoregressive (SD-ARX) model is intensively applied in complex nonlinear system modeling,^17,18 and it belongs to a quasi-LPV model and its state-dependent coefficients and flexible model structure can be convenient for the design of subsequent robust controller.^19,20 By approximating the functional coefficients in an SD-ARX model with radial basis function (RBF) neural networks, the RBF-ARX model can be obtained.^21–23 The SD-ARX model-based RMPC for systems without disturbance was first studied in the work of Peng et al. ²⁴ Based on this type of model, an RMPC with unknown steady-state information was proposed in Zhou et al. ²⁵ To decrease the conservativeness of the LPV models in studies such as Peng et al, ²⁴ Zhou et al. ²⁵ Zhou et al. ²⁶ proposed an RMPC synthesis method using the model’s parameter boundary information, and based on this model, they ²⁷ further proposed a one-stage scheduling RMPC. So far, many meaningful researches of SD-ARX model-based RMPC have been reported,^24–27 whereas the research has mainly focused on studying single-constant feedback RMPC algorithms. According to the research results in Lu et al., ² Li et al., ¹⁰ the improvement of the performance and the enlargement of the feasible region with more freedom of the RMPC can be achieved if variable feedback gains are used, compared with the case of using constant single feedback gains.

Considering that more freedom can be provided by the state feedback control with variable feedback gain for constructing robust controllers, in this paper, we propose a SD-ARX model-based RMPC with variable feedback gain for nonlinear systems. The information of parameter variation rate in the process of constructing the system’s state space models is further explored to both decrease the conservatism and increase the predication accuracy. In this paper, the RMPC algorithm has been designed without using the steady-state information of the system, which differs from those with given by^5–9 or calculated with steady-state information.^12–14 The rest of the work is constructed as follows: the construction process of an identified SD-ARX model-based polytopic state space model (SSM) is presented in Section II, the design of the RMPC with variable feedback gain is introduced in Section III, the numerical example is illustrated in Section IV. Finally, the conclusions are described in Section V.

Design of polytopic state space models

Considering that the nonlinearity of many practical problems depends on one or more system variables,^17,18 the SD-ARX model ¹⁹ can be used to present the dynamics of the system. By approximating the functional coefficients with RBF networks, the following RBF-ARX model can be obtained ²⁸ :

$$\begin{array}{c}y(t+1)=\sum _{k=0}^{n_y-1}{\pi }_{y,k}(\mathbf{w}(t))y(t-k)+\sum _{k=0}^{n_u-1}{\pi }_{u,k}(\mathbf{w}(t))u(t-k)\hfill \\ +{\pi }_0(\mathbf{w}(t))+\xi (t+1)\hfill \\ \{\begin{array}{c}{\pi }_{j,k}(\mathbf{w}(t))=c_{k,0}^j+\sum _{m=1}^h{c}_{k,m}^j\exp \{-{\lambda }_m^j\Vert \mathbf{w}(t)-{\mathbf{z}}_m^j{\Vert }_2^2\},\hfill \\ \multicolumn{1}{c}{{\pi }_0(\mathbf{w}(t))=c_0^0+\sum _{m=1}^h{c}_m^0\exp \{-{\lambda }_m^y\Vert \mathbf{w}(t)-{\mathbf{z}}_m^y{\Vert }_2^2\},j=y,u.}\end{array}\hfill \end{array}$$

(1)

where $u(t)$ denotes the input of the system; $y(t)$ denotes the output; $\xi (t+1)$ is the white noise item; $n_u$ is the maximal input lag; $n_y$ is the maximal output lag; the process variable $\mathbf{w}(t)={[w(t),w(t-1),\cdots ,w(t-d+1)]}^{\mathrm{T}}$ is a measurable state vector that causes the operating point of the system to vary with time, in many cases, it contains maybe the output series and/or the input series^17,28, $\{{\pi }_0(\mathbf{w}(t)),{\pi }_{j,k}(\mathbf{w}(t))|j=y,u\}$ are the RBF network-type functional coefficients composed of multiple exp-functions; $\Vert \mathbf{x}{\Vert }_2^2={\mathbf{x}}^{\mathrm{T}}\mathbf{x}$ is the vector 2-norm; $\{c_0^0,c_m^0,c_{k,0}^j,c_{k,m}^j\}$ are linear weights of the RBF networks; and ${\mathbf{z}}_m^j$ and ${\mathrm{\lambda }}_m^j$ are the RBF network centers and scaling parameters, respectively. The optimization of the linear parameters ${\theta }_L=\{c_0^0,c_m^0,c_{k,0}^j,c_{k,m}^j\}$ and nonlinear parameters ${\theta }_N=\{{\mathbf{z}}_m^j,{\mathrm{\lambda }}_m^j\}$ of model (1) can be conducted by using the variable projection algorithm, ²⁰ the normal ¹⁷ or the regularized structured nonlinear parameter optimization method (SNPOM). ²⁸

If the bounded external disturbance of the system is considered, then the following more general system model can be obtained:

$$\begin{array}{c}y(t+1)=\sum _{k=0}^{k_n-1}(a_{k+1,t}y(t-k)+b_{k+1,t}u(t-k))\hfill \\ +{\pi }_{0,t}+\zeta (t+1)\hfill \\ \{\begin{array}{c}{a}_{k+1,t}=\{\begin{array}{c}{\pi }_{y,k}(\mathbf{w}(t)),0\le k\le n_y-1\hfill \\ \multicolumn{1}{c}{0,n_y\le k\le k_n}\end{array},\hfill \\ \multicolumn{1}{c}{b_{k+1,t}=\{\begin{array}{c}{\pi }_{u,k}(\mathbf{w}(t)),0\le k\le n_u-1\hfill \\ \multicolumn{1}{c}{0,n_u\le k\le k_n}\end{array},}\\ \multicolumn{1}{c}{{\pi }_{0,t}={\pi }_0(\mathbf{w}(t)),k_n=max(n_y,n_u),|\zeta (t+1)|\le \delta .}\end{array}\hfill \end{array}$$

(2)

where $\delta$ is a known upper bound and $\zeta (t+1)$ is the bounded external unknown disturbance, containing modeling errors.

To obtain the polytopic model of the system, we first define the input and output increments as follows:

$$\{\begin{array}{c}\overline{y}(t+i)=y(t+i)-y_r,\hfill \\ \multicolumn{1}{c}{\overline{u}(t+j)=u(t+j)-u(t+j-1).}\end{array}$$

(3)

where $i$ and $j$ are integers less than or equal to zero, and $y_r$ denotes the expected output reference of the system. Then, the one-step forward prediction $\overline{y}(t+1|t)$ can be derived as follows:

$$\overline{y}(t+1|t)=\sum _{k=0}^{k_n-1}(a_{k+1,t}\overline{y}(t-k)+b_{k+1,t}\overline{u}(t-k))+{\phi }_t$$

(4)

$$\begin{array}{c}\multicolumn{1}{c}{{\phi }_t=\sum _{k=0}^{k_n-1}(a_{k+1,t}{y}_r+b_{k+1,t}u(t-k-1))-y_r+{\pi }_{0,t}+\zeta (t+1|t)}\\ ={\psi }_t+\zeta (t+1|t)\hfill \end{array}$$

(5)

where $\zeta (t+1|t)$ denotes the bounded unknown disturbance with $|\zeta (t+1|t)|\le \delta$ . The term $\left|{\phi }_t\right|$ in (5) can be used to indicate the controlled system’s steady state information, since it equals zero under the steady state for a perfect input $u(t)$ and a stabilized output of $y_r$ . ²⁶ Therefore, the future controlled system dynamics can be specified by designing model (6), where the “perfect” input increment set $\{\overline{u}(t),\overline{u}(t+1|t),\dots \}$ will be constructed to achieve future items $\{\left|{\phi }_{t+j|t}\right|,j\ge 1\}=0$ .

$$\begin{array}{c}\overline{y}(t+j+1|t)\hfill \\ =\sum _{k=0}^{k_n-1}(a_{k+1,t+j}\overline{y}(t+j-k|t)+b_{k+1,t+j}\overline{u}(t+j-k|t))\hfill \end{array}$$

(6)

where $j\ge 1$ ; if $j\le k$ , $\overline{u}(t+j-k|t)=\overline{u}(t+j-k)$ , and $\overline{y}(t+j-k|t)=\overline{y}(t+j-k)$ .

To establish the state space model of the system, first, the state vector of the system is defined as follows:

$$\{\begin{array}{c}X(t+j|t)={[x_{1,t+j|t}{x}_{2,t+j|t}\cdots x_{k_n,t+j|t}]}^{\mathrm{T}},\\ x_{1,t+j|t}=\overline{y}(t+j|t),\\ x_{i,t+j|t}=\sum _{k=1}^{k_n+1-i}(a_{k+i-1,t+j-1}\overline{y}(t+j-k|t)\\ +b_{k+i-1,t+j-1}\overline{u}(t+j-k|t))\end{array}.$$

(7)

where $j\ge 0$ and $i\in \{2,3,\dots ,k_n\}$ .

Thus, the state space models corresponding to the polynomial models (4) and (6) can be obtained as follows:

$$\{\begin{array}{c}X(t+1|t)=A_{t}X(t|t)+B_t\overline{u}(t|t)+\mathrm{\Xi }(t),\hfill \\ \multicolumn{1}{c}{A_t=\left[\begin{array}{ccccc}{a}_{1,t}\hfill & 1\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{a_{2,t}}& 0\hfill & 1\hfill & \cdots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & 0\hfill \\ \multicolumn{1}{c}{a_{k_n-1,t}}& 0\hfill & 0\hfill & \cdots \hfill & 1\hfill \\ \multicolumn{1}{c}{a_{k_n,t}}& 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right],B_t=\left[\begin{array}{c}{b}_{1,t}\hfill \\ \multicolumn{1}{c}{b_{2,t}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{b_{k_n,t}}\end{array}\right],\mathrm{\Xi }(t)=\left[\begin{array}{c}{\phi }_t\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0}\end{array}\right].}\end{array}$$

(8)

and

$$\{\begin{array}{c}X(t+j+1|t)=A_{t+j|t}X(t+j|t)+B_{t+j|t}\overline{u}(t+j|t),\hfill \\ \multicolumn{1}{c}{A_{t+j|t}=\left[\begin{array}{ccccc}{a}_{1,t+j|t}\hfill & 1\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{a_{2,t+j|t}}& 0\hfill & 1\hfill & \cdots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & 0\hfill \\ \multicolumn{1}{c}{a_{k_n-1,t+j|t}}& 0\hfill & 0\hfill & \cdots \hfill & 1\hfill \\ \multicolumn{1}{c}{a_{k_n,t+j|t}}& 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right],B_{t+j|t}=\left[\begin{array}{c}{b}_{1,t+j|t}\hfill \\ \multicolumn{1}{c}{b_{2,t+j|t}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{b_{k_n,t+j|t}}\end{array}\right].}\end{array}$$

(9)

Based on model (1) and state vector $\mathbf{w}(t)$ at the current time, it is clearly that $X(t|t)$ , $A_t$ , and $B_t$ in the state space model (8) can be obtained. From (8) and (5), one can see that $\mathrm{\Xi }(t)$ is an unknown item. But it varies in a determined range. In this paper, a convex polytopic set ${\mathrm{\Omega }}_{\mathrm{\Xi }}$ is designed for wrapping the range of $\mathrm{\Xi }(t)$ , in which two vertices of ${\mathrm{\Omega }}_{\mathrm{\Xi }}$ are defined as follows:

$$\{\begin{array}{c}{\mathrm{\Xi }}_{t,1}={[{\psi }_t-\delta ,0,\cdots ,0]}^{\mathrm{T}},\hfill \\ \multicolumn{1}{c}{{\mathrm{\Xi }}_{t,2}={[{\psi }_t+\delta ,0,\cdots ,0]}^{\mathrm{T}}.}\end{array}$$

(10)

• Remark 1: $\{A_{t+j|t},B_{t+j|t}|j\ge 1\}$ in model (9) consists of the functional coefficients $\{a_{k+1,t+j|t},b_{k+1,t+j|t}|k=0,...,k_n-1\}$ , which depend on the future state vectors. However, more often than not, the future state vectors $\mathbf{w}(t+j|t)$ cannot be obtained. Therefore, using the special SD-ARX structure of model (1), we can achieve the ranges of all exp-functions in the functional coefficients, and then we can construct convex polytopic sets for wrapping the variation range of the set $\{A_{t+j|t},B_{t+j|t}|j\ge 1\}$ in model (9).

First, based on the current state matrix set $\{A_t,B_t\}$ in (8) and the information of the parameter variation rates in model (1) calculated with the global identification data, we can construct two convex polytopic sets to wrap the future state matrices $A_{t+j|t}$ and $B_{t+j|t}$ in (9) as follows:

$${\mathrm{\Omega }}_A^j:=\{A_{t+j|t}=\sum _{\alpha =1}^{L_h}{\lambda }_{\alpha ,j}^a{A}_{\alpha }^j,{\lambda }_{\alpha ,j}^a\ge 0,\sum _{\alpha =1}^{L_h}{\lambda }_{\alpha ,j}^a=1\}$$

(11)

$${\mathrm{\Omega }}_B^{j}j:=\{B_{t+j|t}=\sum _{\beta =1}^{L_h}{\lambda }_{\beta ,j}^b{B}_{\beta }^j,{\lambda }_{\beta ,j}^b\ge 0,\sum _{\beta =1}^{L_h}{\lambda }_{\beta ,j}^b=1\}$$

(12)

where $L_h=2^h$ , $h$ is the model order; and their vertices can be defined as follows:

$$A_{\alpha }^j=\left[\begin{array}{ccccc}{c}_{0,0}^y+\sum _{m=1}^h{c}_{0,m}^y({\overline{e}}_{y,m}^j|{\underset{¯}{e}}_{y,m}^j)\hfill & 1\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{c_{1,0}^y+\sum _{m=1}^h{c}_{1,m}^y({\overline{e}}_{y,m}^j|{\underset{¯}{e}}_{y,m}^j)}& 0\hfill & 1\hfill & \cdots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & ⋮\hfill & \ddots \hfill & 0\hfill \\ \multicolumn{1}{c}{c_{k_n-2,0}^y+\sum _{m=1}^h{c}_{k_n-2,m}^y({\overline{e}}_{y,m}^j|{\underset{¯}{e}}_{y,m}^j)}& 0\hfill & 0\hfill & \cdots \hfill & 1\hfill \\ \multicolumn{1}{c}{c_{k_n-1,0}^y+\sum _{m=1}^h{c}_{k_n-1,m}^y({\overline{e}}_{y,m}^j|{\underset{¯}{e}}_{y,m}^j)}& 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right]$$

(13)

$$B_{\beta }^j=\left[\begin{array}{c}{c}_{0,0}^u+\sum _{m=1}^h{c}_{0,m}^u({\overline{e}}_{u,m}^j|{\underset{¯}{e}}_{u,m}^j)\hfill \\ \multicolumn{1}{c}{c_{1,0}^u+\sum _{m=1}^h{c}_{1,m}^u({\overline{e}}_{u,m}^j|{\underset{¯}{e}}_{u,m}^j)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{c_{k_n-1,0}^u+\sum _{m=1}^h{c}_{k_n-1,m}^u({\overline{e}}_{u,m}^j|{\underset{¯}{e}}_{u,m}^j)}\end{array}\right]$$

(14)

where $\{c_{k,m}^j|j=y,u;k=0,...,k_n-1;m=0,...,h\}$ is the linear parameter set of model (1); the mark $({\overline{e}}_{y,m}^j|{\underset{¯}{e}}_{y,m}^j)$ represents ${\overline{e}}_{y,m}^j$ or ${\underset{¯}{e}}_{y,m}^j$ ; and the mark $({\overline{e}}_{u,m}^j|{\underset{¯}{e}}_{u,m}^j)$ represents ${\overline{e}}_{u,m}^j$ or ${\underset{¯}{e}}_{u,m}^j$ . Their specific structures are as follows:

$$\{\begin{array}{c}{\overline{e}}_{y,m}^j=\min \left\{e_m^y(t)+j\mathrm{\Delta }{\overline{e}}_{y,m},{\overline{e}}_{y,m}\right\},\\ {\underset{\_}{e}}_{y,m}^j=\max \{e_m^y(t)+j\mathrm{\Delta }{\underset{\_}{e}}_{y,m},{\underset{\_}{e}}_{y,m}\},\\ {\overline{e}}_{u,m}^j=\min \left\{e_m^u(t)+j\mathrm{\Delta }{\overline{e}}_{u,m},{\overline{e}}_{u,m}\right\},\\ {\underset{\_}{e}}_{u,m}^j=\max \{e_m^u(t)+j\mathrm{\Delta }{\underset{\_}{e}}_{u,m},{\underset{\_}{e}}_{u,m}\}.\end{array}$$

(15)

where $\{e_m^{\theta }(t),{\overline{e}}_{\theta ,m},{\underset{\_}{e}}_{\theta ,m},\mathrm{\Delta }{\overline{e}}_{\theta ,m},\mathrm{\Delta }{\underset{\_}{e}}_{\theta ,m}|\theta =y,u;m=1,\dots ,h\}$ are defined as follows:

$$\{\begin{array}{c}{e}_m^{\theta }(t)=\exp \{-{\lambda }_m^{\theta }\Vert \mathbf{w}(t)-{\mathbf{z}}_m^{\theta }{\Vert }_2^2\},\hfill \\ {\overline{e}}_{\theta ,m}=max\{e_m^{\theta }(t),\forall \mathbf{w}(t)\},\hfill \\ {\underset{¯}{e}}_{\theta ,m}=min\{e_m^{\theta }(t),\forall \mathbf{w}(t)\},\hfill \\ \mathrm{\Delta }{e}_m^{\theta }(t)=-2{\lambda }_m^{\theta }{(\mathbf{w}(t)-{\mathbf{z}}_m^{\theta })}^{\mathrm{T}}{e}_m^{\theta }(t)\mathrm{\Delta }\mathbf{w}(t),\hfill \\ \mathrm{\Delta }{\overline{e}}_{\theta ,m}=max\{\mathrm{\Delta }{e}_m^{\theta }(t),\forall \mathbf{w}(t)\},\hfill \\ \mathrm{\Delta }{\underset{¯}{e}}_{\theta ,m}=min\{\mathrm{\Delta }{e}_m^{\theta }(t),\forall \mathbf{w}(t)\}.\hfill \end{array}$$

(16)

where $\{{\mathbf{z}}_m^{\theta },{\lambda }_m^{\theta }\}$ are the nonlinear parameters of model (1); $\forall \mathbf{w}(t)$ denotes the entire historical $\mathbf{w}(t)$ in the system identification process; $\mathrm{\Delta }\mathbf{w}(t)$ denotes the varying rate of $\mathbf{w}(t)$ ; $e_m^{\theta }(t)$ denotes the exp-functions in model (8) at the current time; ${\overline{e}}_{\theta ,m}$ and ${\underset{¯}{e}}_{\theta ,m}$ denote the upper and lower boundaries of $e_m^{\theta }(t)$ , respectively; $\mathrm{\Delta }{e}_m^{\theta }(t)$ denotes the variation rate of $e_m^{\theta }(t)$ at the current time; $\mathrm{\Delta }{\overline{e}}_{\theta ,m}$ and $\mathrm{\Delta }{\underset{¯}{e}}_{\theta ,m}$ denote the upper and lower boundaries of the variation rate of $e_m^{\theta }(t)$ , respectively.

Remark 2: Using the information of the parameter variation rates and the boundaries of model (1), the variation ranges $[{\underset{¯}{e}}_{y,m}^j,{\overline{e}}_{y,m}^j]$ and $[{\underset{¯}{e}}_{u,m}^j,{\overline{e}}_{u,m}^j]$ can be calculated by (15) and (16), respectively. Then, the vertices ${\mathrm{\Omega }}_A^j$ or ${\mathrm{\Omega }}_B^j$ of the convex polytopic sets ${\mathrm{\Omega }}_A^j$ or ${\mathrm{\Omega }}_B^j$ that “increase gradually” with the number of prediction steps $j$ can be constructed.

Next, to facilitate the robust controller design, the convex polytopic set, which wraps the range of set $\{A_{t+j|t},B_{t+j|t}\}$ in (9) can be redesigned as follows:

$$\begin{array}{c}{\mathrm{\Omega }}_C^j:=\hfill \\ \{(A_{t+j|t},B_{t+j|t})=\sum _{\mu =1}^{n_p}{\lambda }_{\mu ,j}^c(A_{\mu }^j,B_{\mu }^j),{\lambda }_{\mu ,j}^c\ge 0,\sum _{\mu =1}^{n_p}{\lambda }_{\mu ,j}^c=1\}\hfill \end{array}$$

(17)

where $n_p=4^h$ and the state matrices $(A_{\mu }^j,B_{\mu }^j)$ denote the vertices of ${\mathrm{\Omega }}_C^j$ . From $2^h$ vertices $A_{\alpha }^j$ of ${\mathrm{\Omega }}_A^j$ in (11) and $B_{\beta }^j$ of ${\mathrm{\Omega }}_B^j$ in (12), we can build $2^h\times 2^h$ vertices $(A_{\mu }^j,B_{\mu }^j)$ in (17). Thus, the polytopic set ${\mathrm{\Omega }}_C^j$ , which $\{A_{t+j|t},B_{t+j|t}\}$ belongs to, can be obtained.

Remark 3: The limitation of our method is that it consumes relatively large online computation resources, especially with high model order $h$ . To balance the accuracy and the complexity of the polytopic SSM construction method, the marks in (15) can be simplified as ${\overline{e}}_{y,m}^j={\overline{e}}_{y,m}$ , ${\underset{¯}{e}}_{y,m}^j={\underset{¯}{e}}_{y,m}$ , ${\overline{e}}_{u,m}^j={\overline{e}}_{u,m}$ , and ${\underset{¯}{e}}_{u,m}^j={\underset{¯}{e}}_{u,m}$ for all $j\ge 2$ . So, the corresponding vertices $(A_{\mu },B_{\mu })$ of the convex polytopic set defined as ${\mathrm{\Omega }}_D^1$ , to which the set $\{A_{t+j|t},B_{t+j|t}\}$ in (9) belongs, can be obtained.

Overall, one can get the convex polytopic set ${\mathrm{\Omega }}_C^1$ that wraps the state matrix set $\{A_{t+1|t},B_{t+1|t}\}$ and the convex polytopic set ${\mathrm{\Omega }}_D^1$ that wraps the state matrix set $\{A_{t+j|t},B_{t+j|t}|j\ge 2\}$ . Then their corresponding polytopic LPV models can be obtained to wrap the dynamics of the SSMs (8–9).

Synthesis method of the robust MPC with variable feedback gain

In this section, an output tracking RMPC with variable feedback gain are designed based on the SSMs of the system which is established in the previous section. First, the following cost function is utilized in the RMPC synthesis:

$$\begin{array}{c}\underset{\overline{U}(t)}{min}\underset{\mathrm{\Xi }(t)\in {\mathrm{\Omega }}_{\mathrm{\Xi }},\{A_{t+1|t},B_{t+1|t}\}\in {\mathrm{\Omega }}_C^1,\{A_{t+j|t},B_{t+j|t}|j\ge 2\}\in {\mathrm{\Omega }}_D^1}{max}{J}_0^{\infty }(t)\hfill \\ s.t.(8-9),|\overline{u}(t+j|t)|\le {\overline{u}}_{max},j\ge 0.\hfill \end{array}$$

(18)

where $J_0^{\infty }(t)=\sum _{j=0}^{\infty }\{\Vert \overline{u}(t+j|t){\Vert }_R^2+\Vert X(t+j|t){\Vert }_W^2\}$ ; $R$ denotes the positive constant; $W$ denotes the positive definite diagonal matrix; $\Vert x{\Vert }_R^2=x^{\mathrm{T}}\mathrm{Rx}$ ; $X(t+j|t)$ denotes the state predication at time $t$ ; $\overline{u}(t+j|t)$ denotes the input increments constrained by $|\overline{u}(t+j|t)|\le {\overline{u}}_{max}$ ; and $\overline{U}(t)=\{\overline{u}(t|t),\overline{u}(t+1|t),...\}$ .

Then, by dividing function (18) into two parts, one can obtain:

$$J_0^{\infty }(t)=J_0^1(t)+J_1^{\infty }(t)$$

(19)

where $J_0^1(t)=\Vert \overline{u}(t|t){\Vert }_R^2+\Vert X(t|t){\Vert }_W^2$ and $\overline{u}(t|t)$ is the free decision variable.

In this paper, the control strategy for the RMPC algorithm with variable feedback gain is designed as follows:

$$\Lambda :\overline{u}(t+j|t)=\{\begin{array}{c}{F}_{j}X(t+j|t)1\le j<N-1\hfill \\ \multicolumn{1}{c}{F_{N-1}X(t+j|t)j\ge N-1}\end{array}$$

(20)

where $F_j$ is the variable feedback gain with the infinite control horizon of $j$ , and $N\ge 3$ . For $j\ge N-1$ , $F_j$ in the control strategy $\Lambda$ is always set to be $F_{N-1}$ .

Next, the Lyapunov function of the robust controller is designed as follows:

$$V(t+j)=\Vert X(t+j|t){\Vert }_{P_j}^2$$

(21)

where $j\ge 1$ , and if $j>N-1$ , $P_j$ is always set to be $P_{N-1}$ .

From the selected Lyapunov function $V(t+j)$ , the robust stability conditions of the RMPC can be designed as follows:

$$V(t+j+1)-V(t+j)<-\{\Vert \overline{u}(t+j|t){\Vert }_R^2+\Vert X(t+j|t){\Vert }_W^2\}$$

(22)

Add inequality (22) from $j=1$ to $\infty$ , one can obtain the upper limit of $J_1^{\infty }(t)$ :

$$\underset{\mathrm{\Xi }(t)\in {\mathrm{\Omega }}_{\mathrm{\Xi }},\{A_{t+1|t},B_{t+1|t}\}\in {\mathrm{\Omega }}_C^1,\{A_{t+j|t},B_{t+j|t}|j\ge 2\}\in {\mathrm{\Omega }}_D^1}{max}{J}_1^{\infty }(t)<V(t+1)$$

(23)

Therefore, the optimization problem (18) can be rewritten as follows:

$$\begin{array}{c}\underset{\overline{U}(k),P(t+j),j\ge 1}{min}{J}_0^1(t)+\Vert X(t+1|t){\Vert }_{P_1}^2\hfill \\ s.t.(8-9),(22),|\overline{u}(t+j|t)|\le {\overline{u}}_{max},j\ge 0.\hfill \end{array}$$

(24)

Thus, as shown in Theorem 1, the optimization problem (18) can be converted into a linear programing problem (LPP) for solution.

• Theorem 1: The optimizing problem (18) with the control strategy as in (20) is identical with the following LMI constrained convex optimization problem:

$$\underset{\gamma ,\overline{u}(t|t),\{Y_j,G_j,Q_j,Z_j|j=1,\dots ,N-1\}}{min}\gamma ,$$

(25)

subject to

$$\begin{array}{c}\left[\begin{array}{cccc}{G}_1+G_1^{\mathrm{T}}-Q_1\hfill & *\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{A_{\mu }^1{G}_1+B_{\mu }^1{Y}_1}& Q_2\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{W^{1/2}{G}_1}& 0\hfill & \gamma \mathrm{I}\hfill & *\hfill \\ \multicolumn{1}{c}{R^{1/2}{Y}_1}& 0\hfill & 0\hfill & \gamma \mathrm{I}\hfill \end{array}\right]>0,\hfill \\ (\forall \mu =1,...,n_p)\hfill \end{array}$$

(26)

$$\begin{array}{c}\left[\begin{array}{cccc}{G}_j+G_j^{\mathrm{T}}-Q_j\hfill & *\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{A_{\mu }{G}_j+B_{\mu }{Y}_j}& Q_{j+1}\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{W^{1/2}{G}_j}& 0\hfill & \gamma \mathrm{I}\hfill & *\hfill \\ \multicolumn{1}{c}{R^{1/2}{Y}_j}& 0\hfill & 0\hfill & \gamma \mathrm{I}\hfill \end{array}\right]>0,\hfill \\ (\forall j=2,...,N-1;N\ge 3)\hfill \end{array}$$

(27)

$$\begin{array}{c}\left[\begin{array}{cccc}1\hfill & *\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{A_{t}X(t|t)+B_t\overline{u}(t|t)+{\mathrm{\Xi }}_{t,g}}& Q_1\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{W^{1/2}X(t|t)}& 0\hfill & \gamma \mathrm{I}\hfill & *\hfill \\ \multicolumn{1}{c}{R^{1/2}\overline{u}(t|t)}& 0\hfill & 0\hfill & \gamma \mathrm{I}\hfill \end{array}\right]\ge 0,\hfill \\ (\forall g=1,2)\hfill \end{array}$$

(28)

$$-{\overline{u}}_{max}\le \overline{u}(t|t)\le {\overline{u}}_{max}$$

(29)

$$\begin{array}{c}\left[\begin{array}{cc}{Z}_j\hfill & *\hfill \\ \multicolumn{1}{c}{Y_j}& G_j+G_j^{\mathrm{T}}-Q_j\hfill \end{array}\right]\ge 0,Z_j\le {{\overline{u}}_{max}}^2.\hfill \\ (\forall j=1,...,N-1;N\ge 3)\hfill \end{array}$$

(30)

where $R$ and $W$ denote weights; $A_t$ , $B_t$ , $X(t|t)$ , ${\mathrm{\Xi }}_{t,1}$ , ${\mathrm{\Xi }}_{t,2}$ , and $\{A_{\mu }^1,B_{\mu }^1,A_{\mu },B_{\mu }|\mu =1,...,n_p\}$ can be computed based on the obtained historical system data and model (1). Therefore, the control input subjected to the control system at time $t$ is $u(t)=\overline{u}(t|t)+u(t-1)$ , where $\overline{u}(t|t)$ is the solution of the LPP in Theorem 1.

• Proof of Theorem 1: Substitute (20–21) and model (9) into (22), and one can obtain:

$$\Vert A_{t+j|t}+B_{t+j|t}{F}_j{\Vert }_{P_{j+1}}^2-P_j+{F_j}^{\mathrm{T}}R{F}_j+W<0$$

(31)

Considering that $\{A_{t+1|t},B_{t+1|t}\}\in {\mathrm{\Omega }}_C^1$ and $\{A_{t+j|t},B_{t+j|t}|j\ge 2\}\in {\mathrm{\Omega }}_D^1$ , one can see that (31) is equivalent to (32–33).

$$\Vert A_{\mu }^1+B_{\mu }^1{F}_1{\Vert }_{P_2}^2-P_1+{F_1}^{\mathrm{T}}R{F}_1+W<0$$

(32)

$$\Vert A_{\mu }+B_{\mu }{F}_j{\Vert }_{P_{j+1}}^2-P_j+{F_j}^{\mathrm{T}}R{F}_j+W<0$$

(33)

In (32–33), $\mu =1,2,\dots ,n_p$ ; $j=2,\dots ,N-1$ and $N\ge 3$ . Defining $\{F_j=Y_j{G}_j^{-1},Q_j=\gamma {P_j}^{-1}|j=1,\dots ,N-1\}$ and $\gamma \ge 0$ , and then using the Schur complement, one can obtain the following LMIs.

$$\left[\begin{array}{cccc}{G}_1^{\mathrm{T}}{Q}_1^{-1}{G}_1\hfill & *\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{A_{\mu }^1{G}_1+B_{\mu }^1{Y}_1}& Q_2\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{W^{1/2}{G}_1}& 0\hfill & \gamma \mathrm{I}\hfill & *\hfill \\ \multicolumn{1}{c}{R^{1/2}{Y}_1}& 0\hfill & 0\hfill & \gamma \mathrm{I}\hfill \end{array}\right]>0$$

(34)

$$\left[\begin{array}{cccc}{G}_j^{\mathrm{T}}{Q}_j^{-1}{G}_j\hfill & *\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{A_{\mu }{G}_j+B_{\mu }{Y}_j}& Q_{j+1}\hfill & *\hfill & *\hfill \\ \multicolumn{1}{c}{W^{1/2}{G}_j}& 0\hfill & \gamma \mathrm{I}\hfill & *\hfill \\ \multicolumn{1}{c}{R^{1/2}{Y}_j}& 0\hfill & 0\hfill & \gamma \mathrm{I}\hfill \end{array}\right]>0$$

(35)

Following Zhou et al., ²⁵ one can obtain $0<G_j+G_j^{\mathrm{T}}-Q_j\le G_j^{\mathrm{T}}{Q}_j^{-1}{G}_j$ . Given that $Q_j>0$ , conditions (26) and (27) imply (34) and (35), respectively.

Thus, we can convert the form of the minimization problem (24) into the following form:

$$\begin{array}{c}\underset{\gamma ,\overline{u}(t|t),\{Y_j,G_j,Q_j|j=1,\dots ,N-1\}}{min}\gamma \hfill \\ s.t.(8),(34-35),J_0^1(t)+\Vert X(t+1|t){\Vert }_{P_1}^2\le \gamma ,\hfill \\ |\overline{u}(t+j|t)|\le {\overline{u}}_{max},j\ge 0.\hfill \end{array}$$

(36)

Considering that the vector $\mathrm{\Xi }(t)$ in (8) belongs to the convex polytopic set ${\mathrm{\Omega }}_{\mathrm{\Xi }}$ and introducing (8) into (36), the inequality $J_0^1(t)+\Vert X(t+1|t){\Vert }_{P_1}^2\le \gamma$ can be finally expressed as LMIs (28).

Moreover, referring to Zhou et al., ²⁵ one can convert the input and input increment constraints imposed in (36) into the LMIs (29–30).

In the end, the min-max optimization problem shown in (18) is transformed into the minimization problem $\underset{\gamma ,\overline{u}(t|t),\{Y_j,G_j,Q_j,Z_j|j=1,\dots ,N-1\}}{min}\gamma$ that subjects to the LMIs (26–30).

The computation process of the proposed RMPC algorithm is summarized as follows.

• Step 1: SD-ARX model (1) is used to model the nonlinear system, and the regularized-SNPOM is used to optimize the parameters of the model (1).

• Step 2: Initialize controller parameters $W$ and $R$ , predictive feedback horizon $N$ , upper limit of disturbance $\delta$ and constraint value ${\overline{u}}_{max}$ in (29–30).

• Step 3: At sampling time $t$ , based on the identified model (1) and the current system state vector $\mathbf{w}(t)$ , calculate $A_t$ , $B_t$ , $X(t|t)$ according to (7–8) and calculate ${\mathrm{\Xi }}_{t,1}$ , ${\mathrm{\Xi }}_{t,2}$ using (10).

• Step 4: Based on the identified model (1) and the historical state vectors obtained from the global identification data, calculate vertices $\{A_{\mu }^1,B_{\mu }^1,A_{\mu },B_{\mu }|\mu =1,...,n_p\}$ of the polytopic sets ${\mathrm{\Omega }}_C^1$ and ${\mathrm{\Omega }}_D^1$ in (17) according to (11–16).

• Step 5: Solve the optimization problem (25) to get the input increment $\overline{u}(t|t)$ .

• Step 6: Calculate control law $u(t)=\overline{u}(t|t)+u(t-1)$ and implement it on the controlled system at the current time; Increase $t$ by 1 and go to Step 2.

• Theorem 2: When the robust MPC synthesis method is implemented with a rolling-horizon strategy, the stability can be guaranteed by the feasibility of the LPP (25).

• Proof of Theorem 2: For Models (8–9), because LPP (25) is feasible at time $t$ , the optimized solution can be assumed to be $\{{\gamma }^{*}(t),\overline{u}{(t|t)}^{*},$ ${\Phi }_1^{*},\dots ,{\Phi }_{N-1}^{*}\}$ where ${\Phi }_j^{*}=\{Y_j^{*},G_j^{*},Q_j^{*},Z_j^{*}\}$ . From (22), one can obtain that $V(t+2)+J_1^2(t)<V(t+1)$ . According to $J_0^1(t)+V(t+1)\le \gamma$ , one can obtain that ${\gamma }^{*}(t)\ge J_0^1(t)+$ $V(t+1)\ge V(t+1)>V(t+2)+J_1^2(t)$ . To prove the theorem, we construct a solution at time $t+1$ as $\{J_1^2(t)+V(t+2),Y_1^{*}{G}_1^{*-1}X(t+1|t),\sigma {\Phi }_2^{*},...,\sigma {\Phi }_{N-1}^{*},\sigma {\Phi }_{N-1}^{*}\}$ , where $\sigma =(V(k+2)+J_1^2(t))/{\gamma }^{*}(t)<1$ and $\sigma {\Phi }_j^{*}=\{\sigma Y_j^{*},$ $\sigma G_j^{*},\sigma Q_j^{*},\sigma Z_j^{*}\}$ . Because $0<\sigma <1$ , one can see that the new solution meets conditions (26–27) and (30). And, according to (30), $\overline{u}(t+1|t+1)=\overline{u}(t+1|t)=Y_1^{*}{G}_1^{*-1}X(t+1|t)$ satisfies condition (29). Moreover, based on ${\gamma }^{*}(t+1)=J_1^2(t)+V(t+2)\ge J_0^{\infty }(t+1)$ and $\overline{u}(t+1|t+1)=$ $\overline{u}(t+1|t)$ , the new solution meets condition (28). Therefore, it is definitely a practical solution to LPP (25) at $t+1$ . Besides, it is clearly that ${\gamma }^{*}(t)\ge V(t+1)+J_0^1(t)\ge V(t+1)>V(t+2)$ $+J_1^2(t)={\gamma }^{*}(t+1)$ and one can obtain ${\gamma }^{*}(t)>{\gamma }^{*}(t+1)$ , which completes the proof.

• Remark 4: Based on the uncertain SSMs (8–9) and Theorems 1 and 2, when there is a feasible solution for the LPP (25) at time t, one can achieve the output-tracking RMPC without steady-state offset to nonlinear system (2) because there is a $u_s$ satisfying $u(t)\to u_s$ as $t\to \infty$ since $y(t)\to y_r$ when $X(t)\to 0$ and $\overline{u}(t)\to 0$ as $t\to \infty$ .

A numerical example

Case I – CSTR process

Considering that continuous stirred tank reactor (CSTR) ²⁹ is very important in the core equipment for chemical production, a CSTR process is selected and utilized to verify the performance of the RMPC strategy. The selected CSTR process is an irreversible first-order exothermic reaction ( $A\to B$ ), and its characteristics described by continuous-time nonlinear differential equations are as follows ³⁰ :

$$\{\begin{array}{c}\stackrel{·}{C}a(t)=\frac{q}{v}C{a}_0-\mathrm{Ca}(t)-k_0\mathrm{Ca}(t)e^{(-\frac{E}{\mathrm{RT}(t)})},\hfill \\ \stackrel{·}{T}(t)=\frac{q}{v}{T}_0-\frac{q}{v}T(t)+k_1\mathrm{Ca}(t)e^{(-\frac{E}{\mathrm{RT}(t)})}\hfill \\ +k_2{q}_c(t)(T_{c0}-T(t))(1-e^{(-\frac{k_3}{q_c(t)})}).\hfill \end{array}$$

(37)

where $T(t)$ is the reaction temperature; $q_c(t)$ is the coolant flow rate; $\mathrm{Ca}(t)$ is the product concentration of $A$ ; $T_{c0}$ is the inlet coolant temperature; $T_0$ is the feed temperature; $q$ is the process flow rate; $C{a}_0$ is the feed concentration; $k_0$ is the reaction rate constant; $v$ is the reactor volume; $E/R$ is the activation energy term; $k_1=(\mathrm{\Delta }H{k}_0)/(\rho C_p)$ , $k_2=1/v$ and $k_3=\mathrm{hA}/(\rho C_p)$ are the constant coefficients of the chemical reaction; $\mathrm{\Delta }H$ is the heat of reaction; $\rho$ is the liquid density; $C_p$ is the specific heats; $\mathrm{hA}$ is the heat transfer term. The CSTR parameters in (37) are shown in Table 1. In this process, the reaction temperature $T(t)$ is measured in real time and controlled by the manipulated variable $q_c(t)$ optimized by the proposed robust MPC algorithm.

Table 1.

CSTR parameters in (37).

$C{a}_0$	$1\mathrm{mol}/\mathrm{L}$	$k_1$	$9\times {10}^{12}\mathrm{KL}/\mathrm{cal}/\min$
$q$	$100\mathrm{L}/\min$	$k_2$	$1\times {10}^{-2}/\mathrm{L}$
$T_0$	$350\mathrm{K}$	$k_3$	$7\times {10}^2\mathrm{L}/\min$
$T_{c0}$	$350\mathrm{K}$	$v$	$100\mathrm{L}$
$k_0$	$4.5\times {10}^{10}\mathrm{mi}{\mathrm{n}}^{-1}$	$E/R$	$1\times {10}^4\mathrm{K}$

Polytopic SSM construction process

First, the basic model structure of the process can be selected as:

$$\begin{array}{c}T(t+1)=\sum _{k=0}^{n_y-1}{\pi }_{y,k}(\mathbf{w}(t))T(t-k)+\sum _{k=0}^{n_u-1}{\pi }_{u,k}(\mathbf{w}(t))q_c(t-k)\hfill \\ +{\pi }_0(\mathbf{w}(t))+\xi (t+1)\hfill \\ \{\begin{array}{c}{\pi }_{j,k}(\mathbf{w}(t))=c_{k,0}^j+\sum _{m=1}^h{c}_{k,m}^j\exp \{-{\lambda }_m^j\Vert \mathbf{w}(t)-{\mathbf{z}}_m^j{\Vert }_2^2\},\hfill \\ \multicolumn{1}{c}{{\pi }_0(\mathbf{w}(t))=c_0^0+\sum _{m=1}^h{c}_m^0\exp \{-{\lambda }_m^y\Vert \mathbf{w}(t)-{\mathbf{z}}_m^y{\Vert }_2^2\},j=y,u.}\end{array}\hfill \end{array}$$

(38)

where the system input and the corresponding output are the coolant flow rate $q_c(t)$ and the reaction temperature $T(t)$ , respectively; the measurable state vector is chosen as $\mathbf{w}(t)={[T(t),T(t-1),\cdots ,T(t-d+1)]}^{\mathrm{T}}$ ; and $\{{\pi }_0(\mathbf{w}(t)),$ ${\pi }_{j,k}(\mathbf{w}(t))|j=y,u\}$  are the functional coefficients composed of multiple exp-functions. The linear parameters ${\theta }_L=\{c_0^0,c_m^0,c_{k,0}^j,c_{k,m}^j\}$ and nonlinear parameters ${\theta }_N=\{{\mathbf{z}}_m^j,{\mathrm{\lambda }}_m^j\}$ in model (38) are optimized by the regularized-SNPOM. ²⁸

In this paper, to get a set of globally varying input/output data, a set of sine signals and step signals are set to be the required outputs under the PID controller for the CSTR, which make the input and output of the system adequately vary within a large range, and a set of Gaussian white noise signal with small power is also added to the inputs in order to excite the systems’ various dynamic modes. The global identification data of the CSTR process with a sampling period of 0.05 min are shown in Figure 1. We use the first and the last 2500 data points for training and testing of the model, respectively. In the modeling process, the step response performance of the identified model is evaluated as the criterion for determining the model orders in (38) and the final selected model orders in (38) are $n_y=2$ , $n_u=2$ , $h=1$ and $d=1$ .

Figure 1.

Global identification data of the process.

Based on the state vector defined in (7) and the global identification data of the process, the SSMs of the system can be achieved at each sampling instant. Their poles, varying with the variable $\mathbf{w}(t)$ , are indicated by different colors, as shown in Figure 2. It is obvious that the characteristics of the process vary with variable $\mathbf{w}(t)$ . Therefore, the parameter changing rate and the parameter boundary of the system can be obtained through the identification process. Then, we can obtain the previously designed polytopic SSMs (8–9) in Section II to wrap the dynamics of the system at each sampling time. Finally, based on these models, the proposed RMPC algorithm can be constructed and used for the process control.

Figure 2.

Poles of SSMs varying with $\mathbf{w}(t)$ .

Control results and analysis of the RMPC algorithms

The numerical experiments are conducted under the MATLAB R2017b running on a computer equipped with an Intel Core i7-3770 CPU and 16 GB RAM. In this CSTR process, the constraints of the control variable $q_c(t)$ and the input increment are selected as $10\le q_c(t)\le 150$ and ${\overline{u}}_{max}=60\mathrm{L}/\min$ , respectively. Just as in (18), the constraints of the control variable $u(t)$ is not considered in the cost function of the proposed RMPC algorithm. So, the constraint on the input variable in this paper is the physical constraint of the system. After a series of trials and errors, the determined controller weights of our RMPC are $W=\mathrm{diag}(1,1)$ and $R=0.005$ .

Figure 3 shows the control results of our RMPC for several predictive feedback horizons $N$ (as in (20)). The pink dotted line in Figure 3 represents the expected output trajectory of the system. A square wave disturbance with a $2\mathrm{K}$ amplitude and a 0.5 min duration is added into the system output at the fourth min and ninth min for testing the anti-jamming performance of the proposed RMPC. Therefore, the upper limit $\delta$ , as in (2), is set to 2 for the uncertain system disturbance $\zeta (t+1)$ . Figure 4 shows the real-time computational burden corresponding to the control results in Figure 3.

Figure 3.

Control results of our RMPC for selected predictive horizons.

Figure 4.

Online computational burden corresponding to the results in Figure 3.

From Figure 3, it can be seen that the control results of our RMPC are improved with the increase of the predictive feedback horizon $N$ , especially the anti-interference performance of the system when an external interference is introduced. Moreover, when the feedback horizon $N$ is >5, the improvement of the control effect becomes weak. From Figure 4, one can see that the online computational cost corresponding to the results in Figure 3 also increases obviously with increasing the predictive feedback horizon $N$ . This paper comprehensively considers the control performance of the RMPC algorithm and the online computing burden of the algorithm, the predictive feedback horizon $N$ is finally selected as 7 for the CSTR.

In this subsection, two RBF-ARX model-based RMPC algorithms suggested in Peng et al., ²⁴ Zhou et al. ²⁶ (named RMPC1 and RMPC2, respectively) are also performed for comparative research. To be fair, the same RBF-ARX model obtained previously is used in both RMPC algorithms. After a series of trials and errors, the well-adjusted parameters of the three RMPCs are presented in Table 2.

Table 2.

Controller parameters.

Algorithm	Parameters
RMPC1	$R=0.01$ , $W=\mathrm{diag}(1,1)$
RMPC2	$R=0.008$ , $W=\mathrm{diag}(1,1)$ , $\delta =2$
RMPC(N = 7)	$R=0.005$ , $W=\mathrm{diag}(1,1)$ , $N=7$ , $\delta =2$

The control results of the 3 RMPCs are shown in Figure 5, and the corresponding online computational burden is shown in Figure 6.

Figure 5.

Control results of the three RMPCs for CSTR.

Figure 6.

Online computational burden corresponding to the results in Figure 5.

From Figure 5, one can see that the control results of the RMPC1 are worse than those of the others. The control results of the RMPC2 and RMPC(N = 7) are not significantly different when the system has no external interference. However, when the system is interfered with by the square wave disturbance, the RMPC(N = 7) has a better anti-interference performance than the RMPC2. From Figure 6, it is also clearly that the RMPC(N = 7) shows a heavier online computing burden compared with the other two algorithms. However, its maximal value is still far smaller than that of the control system with a 0.05 min sampling period.

Case II – water tank system

Level control of the water tank system is often used in literatures^30,31 for comparison study. The mass and energy balance equations for the water tank system are given in Zhou et al. ³⁰ In this case study, the control problem is to control the level $H_2(t)$ of the lower tank by adjusting the input voltage $U(t)$ of the pump. The detailed RBF-ARX modeling process and results of the water tank system have been given in Zhou et al., ³⁰ that is, the estimated RBF-ARX model orders are $n_y=5$ , $n_u=4$ , $h=1$ and $d=2$ . Based on this identified model, the proposed RMPC algorithm can be obtained and used for the level control of the water tank system.

The numerical experiments are conducted under the same software and hardware environment as in the previous subsection. The sampling period of the control system is 1.5 s. In this water tank system, the constraints of the control variable $H_2(t)$ and the input increment $U(t)$ are selected as $2\le H_2(t)\le 100$ and ${\overline{u}}_{max}=150$ , respectively. After a series of trials and errors, the determined controller weights of our RMPC are $W=\mathrm{diag}(1,1,1,1,1)$ and $R=0.001$ . Similar to the previous section, the RMPC1 and RMPC2 based on the same RBF-ARX model ( $n_y=5$ , $n_u=4$ , $h=1$ and $d=2$ ) are also performed for comparative research. After a series of trials and errors, the well-adjusted parameters of the three RMPCs are presented in Table 3.

Table 3.

Controller parameters.

Algorithm	Parameters
RMPC1	$R=0.0012$ , $W=\mathrm{diag}(1,1,1,1,1)$
RMPC2	$R=0.001$ , $W=\mathrm{diag}(1,1,1,1,1)$ , $\delta =3$
RMPC(N = 3)	$R=0.001$ , $W=\mathrm{diag}(1,1,1,1,1)$ , $N=3$ , $\delta =3$

The control results of the three RMPCs for the water tank system are shown in Figure 7, and the corresponding online computational burden is shown in Figure 8.

Figure 7.

Control results of the three RMPCs for water tank system.

Figure 8.

Online computational burden corresponding to the results in Figure 7.

From Figure 7, one can see that the control results of the RMPC(N = 3) are better than those of the others, especially in terms of the step response performance. From Figure 8, it is also clearly that the RMPC(N = 3) shows a heavier online computing burden compared with the others. However, its maximal value is still a little smaller than that of the system sampling period 1.5 s. It is worth noting that when $N$ is >3, the calculation time of the RMPC has exceeded the sampling time of the water tank system, so the predictive feedback horizon $N$ is finally selected as 3 for the water tank system. In general, for the actual system, it is recommended to choose the maximum value of $N$ that can complete the RMPC algorithm execution within a sample time of the system as the infinite control horizon.

Conclusions

In this paper, a SD-ARX model-based RMPC algorithm with variable feedback gain for nonlinear systems were proposed. First, the construction process of the polytopic SSMs of the system was introduced. Then, the synthesis method of the proposed RMPC was designed. In general, the proposed approach has two advantages: by using the information of the variation rate of the SD-ARX model parameters, our scheme has the ability to shrink the variation regions of the state matrices in the system’s polytopic SSMs, and by using the control policy with variable feedback gain, it also presents the ability to enlarge the robust controller design freedom. The feasibility and effectiveness of the algorithm were verified by two numerical examples. However, the limitation of the RMPC is its relatively large computation burden, especially as the predictive feedback horizon $N$ increases. In order to reduce the online computation burden, we shall try to study a partially offline computing and partially online synthesizing RMPC in the future research.