Advanced Construction of the Dynamic Matrix in Numerically Efficient Fuzzy MPC Algorithms

Abstract

A method for the advanced construction of the dynamic matrix for Model Predictive Control (MPC) algorithms with linearization is proposed in the paper. It extends numerically efficient fuzzy algorithms utilizing skillful linearization. The algorithms combine the control performance offered by the MPC algorithms with nonlinear optimization (NMPC algorithms) with the numerical efficiency of the MPC algorithms based on linear models in which the optimization problem is a standard, easy-to-solve, quadratic programming problem with linear constraints. In the researched algorithms, the free response obtained using a nonlinear process model and the future trajectory of the control signals is used to construct an advanced dynamic matrix utilizing the easy-to-obtain fuzzy model. This leads to obtaining very good prediction and control quality very close to those offered by NMPC algorithms. The proposed approach is tested in the control system of a nonlinear chemical control plant—a CSTR reactor with the van de Vusse reaction.

1. Introduction

The Model Predictive Control (MPC) algorithms use the model of the control plant to predict the behavior of the process. Thanks to such an approach, the MPC algorithms can be successfully applied in control systems of processes with delays, with the inverse response, with constraints, and for MIMO (Multiple-Input Multiple-Output) processes; see, e.g., [1,2,3,4,5,6,7,8,9,10,11]. This is because the future control signal trajectories are generated by the MPC algorithms in such a way that the predicted, many sampling instants ahead, future behavior of the control plant and the shape of control signals fulfill the assumed criteria. Thus, the optimization problem, solved in each time step of the MPC algorithm to obtain these control signal trajectories, has the following form [1,4,7,9,12]:

$$arg\underset{\Delta \mathit{u}}{min}\left\{J_{\mathrm{MPC}}={(\overline{\mathit{y}}-\mathit{y})}^T·\kappa ·(\overline{\mathit{y}}-\mathit{y})+\Delta {\mathit{u}}^T·\mathsf{\Lambda }·\Delta \mathit{u}\right\}$$

(1)

subject to:

$$\Delta {\mathit{u}}_{min}\le \Delta \mathit{u}\le \Delta {\mathit{u}}_{max},$$

(2)

$${\mathit{u}}_{min}\le \mathit{u}\le {\mathit{u}}_{max},$$

(3)

$${\mathit{y}}_{min}\le \mathit{y}\le {\mathit{y}}_{max},$$

(4)

where:

$$\mathit{y}=\left[\begin{array}{c}{\mathit{y}}^1\\ ⋮\\ {\mathit{y}}^{n_y}\end{array}\right],{\mathit{y}}^j=\left[\begin{array}{c}{y}_{k+1|k}^j\\ ⋮\\ y_{k+p|k}^j\end{array}\right],$$

(5)

$$\overline{\mathit{y}}=\left[\begin{array}{c}{\overline{\mathit{y}}}^1\\ ⋮\\ {\overline{\mathit{y}}}^{n_y}\end{array}\right],{\overline{\mathit{y}}}^j=\left[\begin{array}{c}{\overline{y}}_{k+1|k}^j\\ ⋮\\ {\overline{y}}_{k+p|k}^j\end{array}\right],$$

(6)

$$\Delta \mathit{u}=\left[\begin{array}{c}\Delta {\mathit{u}}^1\\ ⋮\\ \Delta {\mathit{u}}^{n_u}\end{array}\right],\Delta {\mathit{u}}^m=\left[\begin{array}{c}\Delta u_{k+1|k}^m\\ ⋮\\ \Delta u_{k+s-1|k}^m\end{array}\right],$$

(7)

$$\mathit{u}=\left[\begin{array}{c}{\mathit{u}}^1\\ ⋮\\ {\mathit{u}}^{n_u}\end{array}\right],{\mathit{u}}^m=\left[\begin{array}{c}{u}_{k+1|k}^m\\ ⋮\\ u_{k+s-1|k}^m\end{array}\right],$$

(8)

$$\kappa =\left[{\kappa }_1,\dots ,{\kappa }_{n_y}\right]·\mathit{I},{\kappa }_i=\left[{\kappa }_i,\dots ,{\kappa }_i\right],$$

(9)

$$\mathsf{\Lambda }=\left[{\mathsf{\Lambda }}_1,\dots ,{\mathsf{\Lambda }}_{n_u}\right]·\mathit{I},{\mathsf{\Lambda }}_i=\left[{\lambda }_i,\dots ,{\lambda }_i\right],$$

(10)

where $y_{k+i|k}^j$ denotes a value, predicted using a process model the MPC algorithm is based on, of the jth output for the $(k+i)$th sampling instant from the prediction horizon, derived at the kth sampling instant, and ${\overline{y}}_{k+i|k}^j$ is an element of the reference trajectory for the jth output and for the $(k+i)$th sampling instant from the prediction horizon; if reference trajectories constant on the prediction horizon are used, then:

$${\overline{y}}_{k+i|k}^j={\overline{y}}_k^j$$

(11)

where ${\overline{y}}_k^j$ denotes a setpoint value for the jth output; $\Delta u_{k+i|k}^m$ are the decision variables of the optimization problem, being future changes in manipulated variables; ${\kappa }_i$ contains p elements, and ${\kappa }_j\ge 0$ denote the weighting coefficients for the predicted control errors of the jth output; ${\mathsf{\Lambda }}_i$ contains s elements, and ${\lambda }_m\ge 0$ denote the weighting coefficients for the changes of the mth manipulated variable; p is the prediction horizon; s is the control horizon; $n_y$ is the number of output variables; $n_u$ is the number of manipulated variables; $\Delta {\mathit{u}}_{min}$, $\Delta {\mathit{u}}_{max}$ denote the vectors defining the lower and upper bounds of the changes of the control signals; ${\mathit{u}}_{min}$, ${\mathit{u}}_{max}$ denote vectors defining the lower and upper bounds of the values of the control signals; and ${\mathit{y}}_{min}$, ${\mathit{y}}_{max}$ denote the vectors defining the lower and upper bounds of the values of the output variables.

The optimization problem (1)–(4) is formulated and solved in each sampling instant, yielding the optimal vector of future control action. From the vector $\Delta \mathit{u}$, the $\Delta u_{k|k}^m$ elements are extracted and applied in the control system. Then, the described procedure is repeated in the next time step.

1.1. MPC Algorithms Based on Linear Models

The simplest MPC algorithms are based on linear process models; see, e.g., [1,13]. In such a case, the superposition principle holds; therefore, the vector of predicted output values $\mathit{y}$ can be decomposed into two parts:

$$\mathit{y}=\tilde{\mathit{y}}+\mathit{A}·\Delta \mathit{u},$$

(12)

$$\tilde{\mathit{y}}=\left[\begin{array}{c}{\tilde{\mathit{y}}}^1\\ ⋮\\ {\tilde{\mathit{y}}}^{n_y}\end{array}\right],{\tilde{\mathit{y}}}^j=\left[\begin{array}{c}{\tilde{y}}_{k+1|k}^j\\ ⋮\\ {\tilde{y}}_{k+p|k}^j\end{array}\right],$$

(13)

where $\tilde{\mathit{y}}$ is called the free response of the control plant describing the influence of the past values of control signals on the process and $\mathit{A}·\Delta \mathit{u}$ is called the forced response, depending on the future changes of the control signals $\Delta \mathit{u}$; the matrix $\mathit{A}$ is called the dynamic matrix and has the following form:

$$\mathit{A}=\left[\begin{array}{cccc}{\mathit{A}}^{11}& {\mathit{A}}^{12}& \dots & {\mathit{A}}^{1n_u}\\ {\mathit{A}}^{21}& {\mathit{A}}^{22}& \dots & {\mathit{A}}^{2n_u}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathit{A}}^{n_{y}1}& {\mathit{A}}^{n_{y}2}& \dots & {\mathit{A}}^{n_y{n}_u}\end{array}\right],$$

(14)

$${\mathit{A}}^{jm}=\left[\begin{array}{ccccc}{a}_1^{j,m}& 0& \dots & 0& 0\\ a_2^{j,m}& a_1^{j,m}& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_p^{j,m}& a_{p-1}^{j,m}& \dots & a_{p-s+2}^{j,m}& a_{p-s+1}^{j,m}\end{array}\right],$$

(15)

where $a_i^{j,m}$ denote the step response coefficients of the control plant, describing the influence of the mth control on the jth output; for details, see, e.g., [9].

After applying the prediction based on a linear model (12), the performance function (1) is transformed into a function that depends quadratically on the decision variables $\Delta \mathit{u}$:

$$J_{LMPC}={(\overline{\mathit{y}}-\tilde{\mathit{y}}-\mathit{A}·\Delta \mathit{u})}^T·\kappa ·(\overline{\mathit{y}}-\tilde{\mathit{y}}-\mathit{A}·\Delta \mathit{u})+\Delta {\mathit{u}}^T·\mathsf{\Lambda }·\Delta \mathit{u}.$$

(16)

The prediction (12) depends linearly on the decision variables; therefore, after using it in the constraints on the output values (4), the optimization problem (1)–(4) becomes an easy-to-solve, standard quadratic optimization problem, with linear constraints. Unfortunately, the control performance offered by an LMPC algorithm, applied to a nonlinear process, may be unsatisfactory. In order to improve it, one can use an MPC algorithm based on a nonlinear model.

1.2. MPC Algorithms Based on Nonlinear Models

Assume that we have a nonlinear process model:

$${\widehat{\mathit{y}}}_k=\mathit{f}({\mathit{y}}_{k-1},{\mathit{y}}_{k-2},\dots ,{\mathit{y}}_{k-n_a},{\mathit{u}}_{k-1},{\mathit{u}}_{k-2},\dots ,{\mathit{u}}_{k-n_b}),$$

(17)

where ${\mathit{y}}_{k-i}={\left[y_{k-i}^1,\dots ,y_{k-i}^{n_y}\right]}^T$ is a vector that contains output values measured at the $(k-i)$th sampling instant and ${\mathit{u}}_{k-i}={\left[u_{k-i}^1,\dots ,u_{k-i}^{n_u}\right]}^T$ is a vector that contains control values applied at the $(k-i)$th sampling instant; denote the output values generated by the model in the $(k+i)$th sampling instant as ${\widehat{\mathit{y}}}_{k+i}={\left[{\widehat{y}}_{k+i}^1,\dots ,{\widehat{y}}_{k+i}^{n_y}\right]}^T$; $n_a,n_b$ determine how many past output and control values the model needs. If one wants to use the model (17) directly in the optimization problem (1)–(4), then the prediction takes the form of the following formulas, passed to the optimization problem as a set of equality constraints:

$$\begin{array}{ccc}\hfill {\mathit{y}}_{k+i|k}=& \mathit{f}& ({\widehat{\mathit{y}}}_{k+i-1},\dots ,{\widehat{\mathit{y}}}_{k+1},{\mathit{y}}_k,\dots ,{\mathit{y}}_{k-n_a+i-1},\hfill \\ & & {\mathit{u}}_{k+i-1|k},\dots ,{\mathit{u}}_{k|k},{\mathit{u}}_{k-1},\dots ,{\mathit{u}}_{k-n_b+i-1})+{\mathit{d}}_k,\hfill \end{array}$$

(18)

where ${\mathit{u}}_{k+i|k}={\left[u_{k+i|k}^1,\dots ,u_{k+i|k}^{n_u}\right]}^T$ is a vector containing future control values, depending on the decision variables from the vector $\Delta \mathit{u}$ and:

$${\mathit{d}}_k={\mathit{y}}_k-{\widehat{\mathit{y}}}_k,$$

(19)

where ${\mathit{y}}_k={\left[y_k^1,\dots ,y_k^{n_y}\right]}^T$ is a vector of recently measured output values; it is assumed that ${\mathit{d}}_k$ is the same for all instants in the prediction horizon—an approach proposed in the Dynamic Matrix Control (DMC) algorithm and therefore called the DMC-type disturbance model; see, e.g., [9].

Unfortunately, the optimization problem (1)–(4) with the prediction (18), based on a nonlinear model, is, in general, a non-convex nonlinear optimization problem; see, e.g., [14,15,16]. Thus, it is a hard-to-solve, time-consuming computational task, and there is no guarantee of finding a global solution, while the time needed to obtain the solution cannot be foreseen in advance. One of the methods to overcome this problem is to use so-called fast NMPC algorithms in which a suboptimal solution is generated faster than in the standard approach; see, e.g., [17,18,19]. To this group of methods belongs also the explicit approach in which most of the calculations are done off-line; see, e.g., [20,21,22] (it is interesting that this approach is optimal if the linear model is used [23,24,25]). Unfortunately, in this approach, the complexity of the controller grows significantly with the number of constraints taken into consideration.

If the process behavior is described by means of a fuzzy model, then the standard NMPC approach can be used, but also the structure of the model can be exploited to formulate algorithms that are easier to solve. One group of such algorithms is based on Linear Matrix Inequalities (LMIs); see, e.g., [26,27,28,29,30]. The other group is based on the classical fuzzy Takagi–Sugeno approach in which a few algorithms in the form of control laws, based on linear models, are used to obtain the fuzzy controller; see, e.g., [31,32].

In the other method, the linearization of the process model for the MPC algorithm is obtained in each time step, and the linear prediction relative to control changes is formulated; see, e.g., [9,33,34,35,36,37,38,39]. As a result, the optimization problem solved by the MPC algorithm in each time step is formulated as the quadratic one (like in LMPC algorithms). In the linearization-based algorithms, the prediction can be based on a classical linearization or a method of prediction generation can exploit the structure of the nonlinear model on which it is based. In the algorithms using the fuzzy Takagi–Sugeno model, described in [36,39], both the free response and the dynamic matrix are obtained using the model obtained after the linearization. In the algorithms detailed in [33,37], the (classical) free response is calculated using the nonlinear model. In [12], the advanced free response, calculated using the nonlinear model (which can have any form of the model generating outputs on the basis of input signals), takes into consideration the previously calculated trajectory of the future control signals (it can be improved iteratively if needed; the approach is similar, though slightly different in details, to the iterative prediction improvement in the iterative learning-based approaches to batch control described in [40,41,42]); the dynamic matrix is generated using the easy-to-obtain fuzzy model.

The approach proposed in the article is an extension of the algorithms presented in [12]. It is done in such a way that the optimization problem solved by the MPC algorithm in each time step is the quadratic one. The modification is introduced in the method of dynamic matrix construction. In [12], elements in the dynamic matrix were obtained taking into consideration the current operation point. Changes of the operating point on the prediction horizon are taken into consideration in the proposed approach. This is done using the free response, generated using the nonlinear control plant model. The obtained prediction is, however, still linear with respect to control changes. Thus, the algorithms utilizing the proposed approach combine the computational simplicity of the LMPC algorithms with the control performance offered by the NMPC algorithms.

The next section details the formulation of the Fuzzy MPC (FMPC) algorithms based on fuzzy and nonlinear models, with the advanced construction of the dynamic matrix. In Section 3, the operation of the FMPC algorithms exploiting the proposed advanced dynamic matrix is tested in a simulation example of the control system of a nonlinear chemical reactor with the inverse response. Conclusions are presented in the last section.

2. Efficient Fuzzy MPC Algorithms with the Advanced Construction of the Dynamic Matrix

The complications resulting from the need to solve a nonlinear optimization task by the MPC algorithm in each time step can be avoided by using an approximation of the process model carried out in each time step. Then, a nonlinear model is used to obtain the free response in such a way that the prediction is linear relative to decision variables, but mimics the nonlinearity of the process very well. An easy-to-obtain fuzzy model is used to get the dynamic matrix.

The dynamic matrix can be constructed in different ways. In Section 2.2, it is described how to do it in such a way that it fits the nonlinearity of the process better than in the standard approach. The method is based on a skillful utilization of the free response (generated using the nonlinear process model). Then, using both elements needed to obtain the prediction, the free response and the dynamic matrix, the optimization problem (1)–(4), solved in each time step, is formulated as a quadratic optimization problem. Such a problem is easy to solve using commonly available optimization routines. Moreover, the simplified versions of the algorithms, in which a control law is obtained, are also discussed.

The proposed algorithms are the generalization of their counterparts proposed in [12]. They use the same free response generation (reminded in Section 2.1) and the same formulation of the optimization problems solved at each time step (reminded in Section 2.3). The difference is in the way the dynamic matrix is constructed (the topic detailed in Section 2.2.3).

2.1. Generation of the Free Response

The method of advanced free response generation, based on a nonlinear process model, proposed in [12], is reminded in this subsection. Define the following vectors:

$${\mathit{u}}_{k+i|k-1}={\left[u_{k+i|k-1}^1,\dots ,u_{k+i|k-1}^{n_u}\right]}^T,$$

(20)

where $u_{k+i|k-1}^m$ are elements of control signal trajectories obtained in the last $(k-1)$st time step.

First, the components of the free response are derived iteratively using the nonlinear model (17) and the vectors (20):

$$\begin{array}{c}{\widehat{\mathit{y}}}_{k+1}=\mathit{f}({\mathit{y}}_k,{\mathit{y}}_{k-1},\dots ,{\mathit{y}}_{k-n_a+1},{\mathit{u}}_{k|k-1},{\mathit{u}}_{k-1},\dots ,{\mathit{u}}_{k-n_b+1});\\ {\widehat{\mathit{y}}}_{k+2}=\mathit{f}({\widehat{\mathit{y}}}_{k+1},{\mathit{y}}_k,\dots ,{\mathit{y}}_{k-n_a+2},{\mathit{u}}_{k+1|k-1},{\mathit{u}}_{k|k-1},{\mathit{u}}_{k-1},\dots ,{\mathit{u}}_{k-n_b+2});\\ ⋮\\ {\widehat{\mathit{y}}}_{k+i}=\mathit{f}({\widehat{\mathit{y}}}_{k+i-1},{\widehat{\mathit{y}}}_{k+i-2},\dots ,{\mathit{y}}_{k-n_a+i},{\mathit{u}}_{k+i-1|k-1},{\mathit{u}}_{k+i-2|k-1},\dots ,{\mathit{u}}_{k-n_b+i});\\ ⋮\end{array}$$

(21)

Note that the values ${\widehat{\mathit{y}}}_{k+1}$ are used when calculating the output values ${\widehat{\mathit{y}}}_{k+2}$, and in general, in the ith iteration, the values ${\widehat{\mathit{y}}}_{k+1},\dots ,{\widehat{\mathit{y}}}_{k+i-1}$ are used to obtain ${\widehat{\mathit{y}}}_{k+i}$.

Next, the final form of the free response is obtained, after taking into account the estimation of unmeasured disturbances ${\mathit{d}}_k={\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-1}$:

$${\tilde{\mathit{y}}}_{k+i|k}={\widehat{\mathit{y}}}_{k+i}+{\mathit{d}}_k.$$

(22)

Note that the free response can be iteratively improved using the nonlinear model. Moreover, it can be modified to include information about measured disturbances. These topics are detailed in [12,43].

2.2. Generation of the Dynamic Matrix

2.2.1. Fuzzy Model Used to Obtain the Dynamic Matrix

The fuzzy Takagi–Sugeno model used to generate the dynamic matrix has the following form:

$$\begin{array}{c}\hfill \mathrm{Rule}f:\\ \hfill & \mathrm{if}{y}_{k-1}^{j_y}\mathrm{is}{B}_1^{f,j_y}\mathrm{and}\dots \mathrm{and}{y}_{k-n}^{j_y}\mathrm{is}{B}_n^{f,j_y}\mathrm{and}\hfill \\ \hfill & u_{k-1}^{j_u}\mathrm{is}{C}_1^{f,j_u}\mathrm{and}\dots \mathrm{and}{u}_{k-m}^{j_u}\mathrm{is}{C}_m^{f,j_u}\hfill \\ \hfill & \mathrm{then}{\widehat{y}}_k^{j,f}=\sum _{m=1}^{n_u}\sum _{n=1}^{p_d-1}{a}_n^{j,m,f}·\Delta u_{k-n}^m+a_{p_d}^{j,m,f}·u_{k-p_d}^m,\hfill \end{array}$$

(23)

where $y_{k-i}^{j_y}$ is the value of the $j_y$th output variable at the $\left(k-i\right)$th sampling instant, $u_{k-i}^{j_u}$ is the value of $j_u$th manipulated variable at the $\left(k-i\right)$th sampling instant, $B_1^{f,j_y},\dots ,B_n^{f,j_y}$, $C_1^{f,j_u},\dots ,C_m^{f,j_u}$ are fuzzy sets, $a_n^{j,m,f}$ are the coefficients of step responses in the fth local model, $j_y=1,\dots ,n_y$, $j_u=1,\dots ,n_u$, $f=1,\dots ,l,$ and l is number of fuzzy rules. The model is composed of local models in the form of step responses [39]; thus, it can be obtained relatively easily, because it is sufficient to collect a few sets of step responses, near a few operating points, e.g., using a nonlinear model of the process. Next, the premises of the model should be formulated using expert knowledge and then tuned, e.g., by means of a fuzzy neural network [9].

2.2.2. Standard Dynamic Matrix

In order to calculate the output values of the model (23), the normalized firing strengths of the fuzzy rules should be calculated using fuzzy reasoning; see, e.g., [44,45]. These values are calculated using the previous values of the output and control signals. In the kth sampling instant, the normalized firing strengths ${\tilde{w}}_k^f$ are obtained, then the output values of the model are as follows:

$${\widehat{y}}_k^j=\sum _{m=1}^{n_u}\sum _{n=1}^{p_d-1}{\tilde{a}}_n^{j,m}·\Delta u_{k-n}^m+{\tilde{a}}_{p_d}^{j,m}·u_{k-p_d}^m,$$

(24)

where:

$${\tilde{a}}_n^{j,m}=\sum _{f=1}^l{\tilde{w}}_k^f·a_n^{j,m,f}.$$

(25)

In the standard approach to the dynamic matrix generation, exploited, e.g., in [12,46,47], the parameters ${\tilde{a}}_n^{j,m}$ are used in the construction of the dynamic matrix in the same way as in the LMPC algorithm, i.e., they are used at each sampling instant from the prediction horizon despite being calculated using the firing strengths ${\tilde{w}}_k^f$ obtained in the current (kth) sampling instant. Thus, they are obtained for the current operating point, which in general, changes on the prediction horizon. Therefore, now the improved version of the dynamic matrix, adopting the nonlinearity of the process on the prediction horizon, will be proposed. It can be relatively easily obtained using the elements of the free response and of the trajectory of future control signals. The method is detailed below.

2.2.3. Advanced Dynamic Matrix Generation

In the first time step from the prediction horizon, the fuzzy model (generating outputs for the $(k+1)$st time step) has the following form:

$$\begin{array}{c}\hfill \mathrm{Rule}f:\\ \hfill & \mathrm{if}{y}_k^{j_y}\mathrm{is}{B}_1^{f,j_y}\mathrm{and}{y}_{k-1}^{j_y}\mathrm{is}{B}_2^{f,j_y}\mathrm{and}\dots \mathrm{and}{y}_{k-n+1}^{j_y}\mathrm{is}{B}_n^{f,j_y}\mathrm{and}\hfill \\ \hfill & u_k^{j_u}\mathrm{is}{C}_1^{f,j_u}\mathrm{and}{u}_{k-1}^{j_u}\mathrm{is}{C}_2^{f,j_u}\mathrm{and}\dots \mathrm{and}{u}_{k-m+1}^{j_u}\mathrm{is}{C}_m^{f,j_u}\hfill \\ \hfill & \mathrm{then}{\widehat{y}}_{k+1}^{j,f}=\sum _{m=1}^{n_u}\sum _{n=1}^{p_d-1}{a}_n^{j,m,f}·\Delta u_{k-n+1}^m+a_{p_d}^{j,m,f}·u_{k-p_d+1}^m,\hfill \end{array}$$

(26)

However, $u_k^{j_u}$ are not known yet, but the approach used in the free response generation can be applied here, i.e., the appropriate element of the future trajectory of the control signals can be used, namely: $u_{k|k-1}^{j_u}$. Therefore, the premises used to obtain the firing strengths of the fuzzy rules will be as follows:

$$\begin{array}{c}\hfill \mathrm{Rule}f:\\ \hfill & \mathrm{if}{y}_k^{j_y}\mathrm{is}{B}_1^{f,j_y}\mathrm{and}{y}_{k-1}^{j_y}\mathrm{is}{B}_2^{f,j_y}\mathrm{and}\dots \mathrm{and}{y}_{k-n+1}^{j_y}\mathrm{is}{B}_n^{f,j_y}\mathrm{and}\hfill \\ \hfill & u_{k|k-1}^{j_u}\mathrm{is}{C}_1^{f,j_u}\mathrm{and}{u}_{k-1}^{j_u}\mathrm{is}{C}_2^{f,j_u}\mathrm{and}\dots \mathrm{and}{u}_{k-m+1}^{j_u}\mathrm{is}{C}_m^{f,j_u}\hfill \end{array}$$

(27)

As a result, the normalized firing strengths of the fuzzy rules ${\tilde{w}}_{k+1|k}^f$ will be obtained for the ${(k+1)}^{\mathrm{st}}$ time step from the prediction horizon (calculated at the kth sampling instant, using the fuzzy reasoning); the appropriate elements of the dynamic matrix will be then calculated using the formula:

$${\tilde{a}}_{n,k+1}^{j,m}=\sum _{f=1}^l{\tilde{w}}_{k+1|k}^f·a_n^{j,m,f}.$$

(28)

In the second time step from the prediction horizon, the fuzzy model has the following form:

$$\begin{array}{c}\hfill \mathrm{Rule}f:\\ \hfill & \mathrm{if}{y}_{k+1}^{j_y}\mathrm{is}{B}_1^{f,j_y}\mathrm{and}{y}_k^{j_y}\mathrm{is}{B}_2^{f,j_y}\mathrm{and}\dots \mathrm{and}{y}_{k-n+2}^{j_y}\mathrm{is}{B}_n^{f,j_y}\mathrm{and}\hfill \\ \hfill & u_{k+1}^{j_u}\mathrm{is}{C}_1^{f,j_u}\mathrm{and}{u}_k^{j_u}\mathrm{is}{C}_2^{f,j_u}\mathrm{and}\dots \mathrm{and}{u}_{k-m+2}^{j_u}\mathrm{is}{C}_m^{f,j_u}\hfill \\ \hfill & \mathrm{then}{\widehat{y}}_{k+2}^{j,f}=\sum _{m=1}^{n_u}\sum _{n=1}^{p_d-1}{a}_n^{j,m,f}·\Delta u_{k-n+2}^m+a_{p_d}^{j,m,f}·u_{k-p_d+2}^m\hfill \end{array}$$

(29)

This time, however, not only $u_k^{j_u}$ are not known yet, but also $u_{k+1}^{j_u}$ and $y_{k+1}^{j_y}$. One can use the approach from the previous step, i.e., instead of $u_{k+1}^{j_u}$, $u_{k+1|k-1}^{j_u}$ is used. The other problem is that also the values $y_{k+1}^{j_y}$ are not known yet. However, one can use the appropriate elements of the free response here, namely ${\tilde{y}}_{k+1|k}^{j_y}$. Therefore, the premises used to obtain the firing strengths in the next $(k+2)$nd time step from the prediction horizon will have the following form:

$$\begin{array}{c}\hfill \mathrm{Rule}f:\\ \hfill & \mathrm{if}{\tilde{y}}_{k+1|k}^{j_y}\mathrm{is}{B}_1^{f,j_y}\mathrm{and}{y}_k^{j_y}\mathrm{is}{B}_2^{f,j_y}\mathrm{and}\dots \mathrm{and}{y}_{k-n+2}^{j_y}\mathrm{is}{B}_n^{f,j_y}\mathrm{and}\hfill \\ \hfill & u_{k+1|k-1}^{j_u}\mathrm{is}{C}_1^{f,j_u}\mathrm{and}{u}_{k|k-1}^{j_u}\mathrm{is}{C}_2^{f,j_u}\mathrm{and}\dots \mathrm{and}{u}_{k-m+2}^{j_u}\mathrm{is}{C}_m^{f,j_u}\hfill \end{array}$$

(30)

the normalized firing strengths ${\tilde{w}}_{k+2|k}^f$, obtained for the $(k+2)$nd time step from the prediction horizon, using (30), are then utilized to calculate the elements of the dynamic matrix corresponding to the $(k+2)$nd time step from the prediction horizon:

$${\tilde{a}}_{n,k+2}^{j,m}=\sum _{f=1}^l{\tilde{w}}_{k+2|k}^f·a_n^{j,m,f}.$$

(31)

This procedure is repeated iteratively in the next time steps from the prediction horizon, resulting in obtaining the normalized firing strengths of the fuzzy rules ${\tilde{w}}_{k+i|k}^f$ in each, $(k+i)$th, time step from the prediction horizon. Next, the appropriate elements of the dynamic matrix, for each $(k+i)$th time step from the prediction horizon, will be calculated using the formula:

$${\tilde{a}}_{n,k+i}^{j,m}=\sum _{f=1}^l{\tilde{w}}_{k+i|k}^f·a_n^{j,m,f}.$$

(32)

Now, the dynamic matrix can be generated (and updated in each sampling instant):

$${\mathit{A}}_k=\left[\begin{array}{cccc}{\tilde{\mathit{A}}}_k^{11}& {\tilde{\mathit{A}}}_k^{12}& \dots & {\tilde{\mathit{A}}}_k^{1n_u}\\ {\tilde{\mathit{A}}}_k^{21}& {\tilde{\mathit{A}}}_k^{22}& \dots & {\tilde{\mathit{A}}}_k^{2n_u}\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{\mathit{A}}}_k^{n_{y}1}& {\tilde{\mathit{A}}}_k^{n_{y}2}& \dots & {\tilde{\mathit{A}}}_k^{n_y{n}_u}\end{array}\right],$$

(33)

$${\tilde{\mathit{A}}}_k^{jm}=\left[\begin{array}{ccccc}{\tilde{a}}_{1,k+1}^{j,m}& 0& \dots & 0& 0\\ {\tilde{a}}_{2,k+2}^{j,m}& {\tilde{a}}_{1,k+2}^{j,m}& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\tilde{a}}_{p,k+p}^{j,m}& {\tilde{a}}_{p-1,k+p}^{j,m}& \dots & {\tilde{a}}_{p-s+2,k+p}^{j,m}& {\tilde{a}}_{p-s+1,k+p}^{j,m}\end{array}\right].$$

(34)

Note that in the ith row of each matrix ${\tilde{\mathit{A}}}_k^{jm}$, the elements calculated for the $(k+i)$th time step from the prediction horizon are used.

Assume that the firing strengths for a given fuzzy rule are grouped in the following vector:

$${\tilde{\mathit{w}}}^f=\left[{\tilde{w}}_{k+1|k}^f,{\tilde{w}}_{k+2|k}^f,\dots ,{\tilde{w}}_{k+p|k}^f\right]$$

(35)

and define:

$${\mathit{A}}^{jmf}=\left[\begin{array}{ccccc}{a}_1^{j,m,f}& 0& \dots & 0& 0\\ a_2^{j,m,f}& a_1^{j,m,f}& \dots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_p^{j,m,f}& a_{p-1}^{j,m,f}& \dots & a_{p-s+2}^{j,m,f}& a_{p-s+1}^{j,m,f}\end{array}\right].$$

(36)

Then, the matrices ${\tilde{\mathit{A}}}_k^{jm}$ can be calculated using the relatively simple formula:

$${\tilde{\mathit{A}}}_k^{jm}=\sum _{f=1}^l{\tilde{\mathit{w}}}^f·\mathit{I}·{\mathit{A}}^{jmf},$$

(37)

where $\mathit{I}$ is the identity matrix of dimension $p\times p$, and the matrices ${\mathit{A}}^{jmf}$ remain the same in each time step.

2.3. Optimization Problem in the Numerical and Analytical Versions of the Algorithms

Now, assume that future control values are decomposed as follows:

$$u_{k+i|k}^m={\check{u}}_{k+i|k}^m+u_{k+i|k-1}^m,$$

(38)

where ${\check{u}}_{k+i|k}^m$ can be interpreted as the corrections of the control signal $u_{k+i|k-1}^m$ obtained in the last $(k-1)$st time step. Thus, future control increments are described by the following, similar formula:

$$\Delta u_{k+i|k}^m=\Delta {\check{u}}_{k+i|k}^m+\Delta u_{k+i|k-1}^m.$$

(39)

Note that the values $u_{k+i|k-1}^m$ and $\Delta u_{k+i|k-1}^m$ are known, as they were calculated in the last time step, and their influence on the output variables is contained in the free response (22), described in Section 2.1. Therefore, now the dynamic matrix (33) will be used to predict the influence of control corrections on the control plant outputs, and after using the free response (22) and the dynamic matrix (33), one obtains the following prediction:

$$\mathit{y}=\tilde{\mathit{y}}+{\mathit{A}}_k·\Delta \check{\mathit{u}}.$$

(40)

where:

$$\Delta \check{\mathit{u}}=\left[\begin{array}{c}\Delta {\check{\mathit{u}}}^1\\ ⋮\\ \Delta {\check{\mathit{u}}}^{n_u}\end{array}\right],\Delta {\check{\mathit{u}}}^m=\left[\begin{array}{c}\Delta {\check{u}}_{k+1|k}^m\\ ⋮\\ \Delta {\check{u}}_{k+s-1|k}^m\end{array}\right].$$

(41)

Note that the prediction depends linearly on corrections $\Delta \check{\mathit{u}}$.

2.3.1. Optimization with the Classical Performance Index

Optimization task (1)–(4), which is solved by the control algorithm in each time step, changes now to the following form with corrections $\Delta \check{\mathit{u}}$ being the decision variables:

$$arg\underset{\Delta \check{\mathit{u}}}{min}\left\{J_{MPC+}={(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{A}}_k·\Delta \check{\mathit{u}})}^T·\kappa ·(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{A}}_k·\Delta \check{\mathit{u}})+\Delta {\mathit{u}}^T·\mathsf{\Lambda }·\Delta \mathit{u}\right\}$$

(42)

subject to:

$$\Delta {\mathit{u}}_{min}\le \Delta \mathit{u}\le \Delta {\mathit{u}}_{max},$$

(43)

$${\mathit{u}}_{min}\le \mathit{u}\le {\mathit{u}}_{max},$$

(44)

$${\mathit{y}}_{min}\le \tilde{\mathit{y}}+{\mathit{A}}_k·\Delta \check{\mathit{u}}\le {\mathit{y}}_{max},$$

(45)

where $\Delta \mathit{u}=\overrightarrow{\Delta \mathit{u}}+\Delta \check{\mathit{u}}$, $\mathit{u}=\overrightarrow{\mathit{u}}+\check{\mathit{u}}$, and:

$$\overrightarrow{\Delta \mathit{u}}=\left[\begin{array}{c}{\overrightarrow{\Delta \mathit{u}}}^1\\ ⋮\\ {\overrightarrow{\Delta \mathit{u}}}^{n_u}\end{array}\right],{\overrightarrow{\Delta \mathit{u}}}^m=\left[\begin{array}{c}\Delta u_{k+1|k-1}^m\\ ⋮\\ \Delta u_{k+s-1|k-1}^m\end{array}\right],\overrightarrow{\mathit{u}}=\left[\begin{array}{c}{\overrightarrow{\mathit{u}}}^1\\ ⋮\\ {\overrightarrow{\mathit{u}}}^{n_u}\end{array}\right],{\overrightarrow{\mathit{u}}}^m=\left[\begin{array}{c}{u}_{k+1|k-1}^m\\ ⋮\\ u_{k+s-1|k-1}^m\end{array}\right]$$

contain elements of the future control increments’ trajectory and of the future control trajectory.

Note that the performance function in (42) depends quadratically on decision variables $\Delta \check{\mathit{u}}$, and all constraints depend linearly on decision variables. Thus, a standard, easy-to-solve linear-quadratic optimization problem is obtained.

If in each time step, the optimization problem with performance function from (42) is solved without constraints, then it has the following solution given by the analytical formula:

$$\Delta \check{\mathit{u}}={\left({\mathit{A}}_k^T·\kappa ·{\mathit{A}}_k+\mathsf{\Lambda }\right)}^{-1}·\left({\mathit{A}}_k^T·\kappa ·\left(\overline{\mathit{y}}-\tilde{\mathit{y}}\right)+\overrightarrow{\Delta \mathit{u}}\right).$$

(46)

2.3.2. Optimization with the Modified Performance Index

In the optimization task (42)–(45), a slightly modified performance index can be used:

$$J_{MPC+v2}={(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{A}}_k·\Delta \check{\mathit{u}})}^T·\kappa ·(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{A}}_k·\Delta \check{\mathit{u}})+\Delta {\check{\mathit{u}}}^T·\mathsf{\Lambda }·\Delta \check{\mathit{u}}.$$

(47)

The change consists of the modification of the second component of the performance index, which now depends only on the corrections of the control signals $\Delta \check{\mathit{u}}$. This modification leads to algorithms that generate faster responses, but the influence of the $\lambda$ parameter on system robustness becomes significantly limited.

If the performance index (47) is minimized without constraints, then the following analytical solution is obtained:

$$\Delta \check{\mathit{u}}={\left({\mathit{A}}_k^T·\kappa ·{\mathit{A}}_k+\mathsf{\Lambda }\right)}^{-1}·{\mathit{A}}_k^T·\kappa ·\left(\overline{\mathit{y}}-\tilde{\mathit{y}}\right).$$

(48)

2.3.3. Utilization of Analytical Versions of the Algorithms

Due to the fact that the dynamic matrix ${\mathit{A}}_k$ changes, in general, in every time step, then in each time step, the formula (46) or (48) allowing calculating the control action should be used by the algorithm. As a consequence, appropriate rows of the matrix ${\left({\mathit{A}}_k^T·\kappa ·{\mathit{A}}_k+\mathsf{\Lambda }\right)}^{-1}$ (corresponding to the control increments for the current sampling instant $\Delta {\mathit{u}}_{k|k}^m$) should be calculated by the algorithm in each time step.

3. Example

3.1. Control Plant

The control plant is a nonlinear, isothermal CSTR in which the van de Vusse reaction takes place. These CSTRs are a popular benchmark due to their nonlinearity and difficult dynamics, willingly used to test newly developed control algorithms; see, e.g., [48,49,50,51,52,53]. The process model of the reactor is composed of two composition balance equations; see, e.g., [54]:

$$\begin{array}{c}\frac{d{C}_{\mathrm{A}}}{dt}=-k_1·C_{\mathrm{A}}-k_3·C_{\mathrm{A}}^2+\frac{F}{V}\left(C_{\mathrm{Af}}-C_{\mathrm{A}}\right),\\ \frac{d{C}_{\mathrm{B}}}{dt}=k_1·C_{\mathrm{A}}-k_2·C_{\mathrm{B}}-\frac{F}{V}{C}_{\mathrm{B}},\end{array}$$

(49)

where $C_{\mathrm{A}}$, $C_{\mathrm{B}}$ are the concentrations of Components A and B, respectively, F is the inlet (and also outlet) flow rate, V is the volume in which the reaction takes place (assumed constant and $V=1\mathrm{L}$), and $C_{\mathrm{Af}}$ is the concentration of Component A in the inlet flow stream (if not declared otherwise, it is assumed that $C_{\mathrm{Af}}=10\mathrm{mol}/\mathrm{L}$). The values of the parameters are: $k_1=501/\mathrm{h}$, $k_2=1001/\mathrm{h}$, $k_3=10\mathrm{l}/(\mathrm{h}·\mathrm{mol})$. The output variable is the concentration $C_{\mathrm{B}}$ of Substance B. The manipulated variable is the inlet flow rate F. $C_{\mathrm{Af}}$ is the disturbance variable.

The described control plant was used during the tests also in [12]. Its nonlinear steady-state characteristic is reminded in Figure 1. The control plant has the inverse response; thus, it is natural to use an MPC algorithm in this case.

The fuzzy model used to generate the dynamic matrix, in the researched algorithms, is the same as in [12]. It is composed of three step responses obtained near the following operating points:

• R1 $C_{B0}=$ 0.91 mol/L, $C_{A0}=$ 2.18 mol/L, $F=$ 20 L/h;
• R2 $C_{B0}=$ 1.12 mol/L, $C_{A0}=$ 3 mol/L, $F=$ 34.3 L/h;
• R3 $C_{B0}=$ 1.22 mol/L, $C_{A0}=$ 3.66 mol/L, $F=$ 50 L/h.

The assumed membership functions are reminded in Figure 2.

3.2. Experiments

For the considered control plant, a few MPC algorithms were designed:

-

NMPC—based on the nonlinear model and nonlinear optimization and numerically efficient FMPC algorithms with advanced free response generation, proposed in [12], and:

-

FMPC1 with the conventional dynamic matrix and the classical performance index,

-

FMPC2 with the conventional dynamic matrix and the modified performance index,

-

FMPC1a with the advanced dynamic matrix and the classical performance index,

-

FMPC2a with the advanced dynamic matrix and the modified performance index.

The simulation experiments were done using MATLAB. During the experiments, the operation of the control systems with the NMPC and FMPC algorithms was compared. The FMPC algorithms use the nonlinear model in the form of state equations, to generate the free response and the fuzzy model (23) to obtain the dynamic matrix. A sampling time equal to $T_{\mathrm{s}}=3.6$ s was assumed; the values of tuning parameters were as follows (if not declared otherwise): prediction horizon $p=70$, control horizon $s=35$, and weighting coefficient $\lambda =0.001$.

The responses of the control systems to changes in the setpoint and to the change of the disturbance by 10% in the 6th minute of the experiment are shown in Figure 3 and Figure 4. If the setpoint was changed to ${\overline{C}}_{B1}=1$ mol/L, the responses obtained in the control system with the FMPC1a algorithm (magenta lines in Figure 3) and with the FMPC1 algorithm (blue lines in Figure 3) were very close to those obtained with the NMPC algorithm with nonlinear optimization (red lines in Figure 3). In all cases, there was almost no overshoot. The FMPC algorithms were slightly faster than the NMPC algorithm.

In the case when the setpoint changed to ${\overline{C}}_{B2}=1.25$ mol/L, the response obtained with the FMPC1 algorithm was the fastest one. The one generated with the NMPC algorithm was the slowest, and the response obtained with the FMPC1a algorithm was between these two responses; it was closer to the response generated with the NMPC algorithm than the response obtained with the FMPC1 algorithm. This is because the prediction used in the FMPC1a algorithm, thanks to using the advanced construction of the dynamic matrix, is more accurate. In the case of all the algorithms, the setpoint was reached without overshoot.

Disturbance responses obtained with the FMPC1 and the FMPC1a algorithms were practically the same. Near $C_B=1$ mol/L, the response generated with the NMPC algorithm was very close to the other responses. In all three cases, there was no overshoot. In the case of operation near $C_B=1.25$ mol/L, the NMPC was faster in the compensation of the disturbance change than the FMPC algorithms.

The FMPC2 and FMPC2a algorithms worked faster than FMPC1 and FMPC1a. However, the same relations between the algorithms could be observed when the responses generated with FMPC2 and FMPC2a algorithms were compared with the ones obtained with the NMPC algorithm (Figure 4). If the setpoint was changed to ${\overline{C}}_{B1}=1$ mol/L, the responses obtained with all the algorithms (FMPC2, FMPC2a, and NMPC) were very close to each other. There was almost no overshoot, and the FMPC2 algorithm was slightly slower than the other algorithms.

In the case when the setpoint changed to ${\overline{C}}_{B2}=1.25$ mol/L, the response obtained with the FMPC2 algorithm was the fastest one. The one generated with the NMPC algorithm was the slowest, and the response obtained with the FMPC2a algorithm, for the same reason as in the case of the FMPC1a algorithm, was between these two responses. In the case of all the algorithms, the setpoint was reached almost without overshoot.

The disturbance responses obtained with the FMPC2 and FMPC2a algorithms were practically the same (the same phenomenon was observed in the case of the FMPC1 and FMPC1a algorithms). Near $C_B=1$ mol/L and $C_B=1.25$ mol/L, the NMPC algorithm generated a bigger maximal error than the FMPC algorithms. In all three cases, there was no overshoot.

It can be noticed that the FMPC1a and FMPC2a algorithms generated responses closer to the ones generated by the NMPC algorithm than the FMPC1 and FMPC2 algorithms. This is because in the algorithms with the advanced construction of the dynamic matrix, the prediction was more accurate, thanks to using the proposed mechanism. It should be, however, emphasized that all the FMPC algorithms used a reliable quadratic programming routine to generate the control action instead of the non-convex, nonlinear optimization utilized in the NMPC algorithm. When comparing the SSE (Sum of Squared Errors), one can notice that the smallest value was obtained with the FMPC2 and FMPC2a algorithms. This is because they used the modified performance index in the optimization problem solved by the algorithm in each time step.

We also performed experiments with the changes of the parameters ($\lambda$ coefficient and the control horizon) of the algorithms, like the ones in [12]. First, the $\lambda$ coefficient was changed to $0.0001$ (Figure 5 and Figure 6). All the algorithms now worked faster than in the previous case (for $\lambda =0.001$); the control action was more aggressive, and the maximal errors at the beginning of the experiment were now bigger. However, the relations between the FMPC and NMPC algorithms remained unchanged. It can be also noticed that the differences between the FMPC1a and FMPC2a (and also FMPC1 and FMPC2) algorithms became much smaller compared to the previous experiments. The FMPC1a algorithm was practically as fast as the FMPC2a one. This is because for $\lambda =0$, both performance indexes (42) and (47) were the same. Thus, the closer the value of the $\lambda$ coefficient to zero is, the closer the responses generated with the FMPC1a and FMPC2a algorithms should be. The SSE values obtained for all the algorithms were smaller than in the previous case, i.e., for $\lambda =0.001$.

To simplify the optimization problem by reducing the number of decision variables, one can decrease the control horizon s. The responses obtained for different values of the control horizon are shown in Figure 7 and Figure 8 for the FMPC1a and FMPC2a algorithms, respectively. In the case of both algorithms, the control horizon can be shortened significantly, because the obtained responses do not change too much when the control horizon is between $s=10$ and $s=35$. However, reduction of the control horizon should be done carefully, because for a too short control horizon, the control performance can get worse.

The FMPC1a algorithm, for $s=5$, compensated the disturbance near $C_B=1.25$ mol/L faster, but achieved the setpoint ${\overline{C}}_{B1}=1$ mol/L slightly more slowly. A further decrease of the control horizon to $s=1$ brought faster disturbance compensation near $C_B=1.25$ mol/L, but worsened the compensation near $C_B=1$ mol/L. Reaching the setpoint ${\overline{C}}_{B2}=1.25$ mol/L was also faster, but at the cost of a slower response to the setpoint change to ${\overline{C}}_{B1}=1$ mol/L. The value of the SSE decreased with the decrease of the control horizon, but the value of the SSE obtained for $s=1$ was noticeably smaller than in the other cases.

The FMPC2a algorithm, for $s=5$ and $s=1$, compensates the disturbance near $C_B=1.25$ mol/L slightly faster, but unfortunately reaching the setpoint ${\overline{C}}_{B1}=1$ mol/L took longer, especially for $s=1$. Disturbance compensation near $C_B=1$ mol/L was also visibly slower for $s=1$. The smallest value of the SSE was obtained for $s=10$; however, the differences in the SSE values in all cases were very small.

4. Conclusions

The mechanism for the advanced construction of the dynamic matrix proposed in the article uses the easy-to-obtain fuzzy model and advanced free response to generate the dynamic matrix skillfully. Thanks to such an approach, the changes of the operating point on the prediction horizon are taken into consideration during the construction of the dynamic matrix. This leads to obtaining a very good prediction, which is, however, linear with respect to the control changes. As a result, the FMPC algorithms can offer control performance very close to the one characteristic for the NMPC algorithms (based on the nonlinear optimization) and, at the same time, numerical efficiency resulting from the fact that the proposed mechanism is designed in such a way that the FMPC algorithms are formulated as the easy-to-solve quadratic optimization problems.

The proposed mechanism can be applied by a control system designer also in the case of simpler FMPC algorithms, e.g., with the classical free response. It should improve the operation of such a controller, but a better prediction can be obtained in the case considered in the article when the advanced free response is used. The quality of the advanced dynamic matrix and of the prediction can be increased even more if iterative improvement of the free response is applied.