Numerically Efficient Fuzzy MPC Algorithm with Advanced Generation of Prediction—Application to a Chemical Reactor

Abstract

In Model Predictive Control (MPC) algorithms, control signals are generated after solving optimization problems. If the model used for prediction is linear then the optimization problem is a standard, easy to solve, quadratic programming problem with linear constraints. However, such an algorithm may offer insufficient performance if applied to a nonlinear control plant. On the other hand, if a model used for prediction is nonlinear, then non–convex optimization problem must be solved at each algorithm iteration. Then the numerical problems may occur during solving it and the time needed to calculate the control signals cannot be determined. Therefore approaches based on linearized models are preferred in practical applications. A fuzzy algorithm with an advanced generation of the prediction is proposed in the article. The prediction is obtained in such a way that the algorithm is formulated as a quadratic optimization problem but offers performance very close to that of the MPC algorithm with nonlinear optimization. The efficiency of the proposed approach is demonstrated in the control system of a nonlinear chemical control plant—a CSTR (Continuous Stirred–Tank Reactor) with van de Vusse reaction.

1. Introduction

The Model Predictive Control (MPC) algorithms during the calculation of the control signals use the model of the control plant directly, predicting the behavior of the process. Thanks to such an approach they can be successfully applied in control systems of control plants with difficult dynamics (delay, inverse response), with constraints and for MIMO (Multiple–Input Multiple–Output) processes, see for example, References [1,2,3,4,5,6,7].

The simplest MPC algorithms are based on linear control plant models; see for example, References [1,8]. Unfortunately, if a control plant is nonlinear and the MPC algorithm should operate for different operating points, as in control system structures with steady–state setpoint optimization (see e.g., References [7,9,10]), then the algorithm based on a linear model may be insufficient in terms of control performance. In such a case operation of the control system can be usually improved using the MPC algorithm based on a nonlinear model.

On the other hand, if a nonlinear process model is used in the MPC algorithm directly (such algorithm is called NMPC—Nonlinear MPC), then the optimization problem solved at each algorithm iteration to generate the control signals is a nonlinear optimization problem; see for example, References [11,12,13]. Thus, it is usually non–convex, hard to solve, and time–consuming computational task. Moreover, the time needed to find the solution cannot be foreseen in advance. One of the ways to cope with this problem is to apply so called fast NMPC algorithms in which the solution is generated faster than in the standard approach, but unfortunately, at the expense of suboptimality; see for example, References [14,15,16]. One of such approaches is the explicit one in which most of calculations are done off–line; see for example, References [17,18,19]. It is an interesting approach and it is optimal if the linear model is used [20,21,22]. However it has a serious drawback—the complexity of the controller grows significantly with the number of inequality constraints in the problem. When the fuzzy models are used, then of course the nonlinear optimization can be applied in the algorithm. However the form of the model is such that easier to solve algorithms can be formulated. Among the algorithms exploiting the structure of the fuzzy models there is a group which is based on Linear Matrix Inequalities (LMIs); see for example, References [23,24,25,26,27]. The other group is based on the classical fuzzy Takagi–Sugeno approach in which a few algorithms based on linear models are used to form the fuzzy controller; see for example, References [28,29].

The other approach consists in online linearization of the process model at each MPC algorithm iteration; see for example, References [7,30,31,32,33,34,35,36]. Then prediction is formulated in such a way that it is linear relative to decision variables in the optimization problem solved at each iteration of the algorithm. Such form of the prediction leads to the formulation of the optimization problem as the quadratic one. In the linearization–based algorithms the prediction can be based on a classical linearization or a clever method of prediction generation can be used, exploiting structure of the nonlinear model the algorithm is based on. In the algorithms using linearization of the fuzzy Takagi–Sugeno model, described in References [31,34], both—the free response and the dynamic matrix were obtained using the model obtained after the linearization. In References [32,36] the model obtained after linearization was used to obtain the dynamic matrix whereas the (classical) free response was calculated using the nonlinear model.

The MPC algorithms proposed in the article are an extension of the algorithm presented in Reference [36] and, like their predecessor, are based on two models—a nonlinear model and its, easy to obtain, fuzzy approximation. The nonlinear model can have any form of the model generating outputs on the base of input signals. This model is used to obtain an advanced (based on the previously calculated trajectory of the future control signals) free response of the process prediction. The easy–to–obtain fuzzy model is used to derive the dynamic matrix describing influence of the future, calculated control signals on prediction. The prediction has such a form that the proposed algorithms are formulated as the standard quadratic programming problem. Thus, the algorithms combine the computational simplicity of the MPC algorithms based on linear models with control performance almost as good as the one offered by the algorithm with nonlinear optimization. The prediction can be also improved iteratively, if needed. Iterative prediction improvement, though slightly different in details, is also used in the interesting iterative–learning–based approaches to batch control, see for example, References [37,38,39].

Standard MPC algorithms based on linear and nonlinear models are described in the next section. In Section 3 the MPC algorithms based on fuzzy and nonlinear models are proposed, and mechanisms of their improvements and simplifications are detailed. A simulation example of the application of the proposed algorithms to a nonlinear chemical reactor with the inverse response is presented in Section 4, illustrating the efficiency of the proposed approach. The article is concluded in the last section.

2. Model Predictive Control Algorithms

Future control signal trajectories are generated in the MPC algorithms in such a way that predicted, many sampling instants ahead, future behavior of the control plant fulfills assumed criteria. It is usually formulated as an optimization problem in which, it is demanded that a performance function is minimized subject to the constraints put on control and output signals. The performance function takes into consideration future, predicted control errors and the shape of control signals. Thus, the optimization problem solved by the MPC algorithm in each iteration has, most often, the following form [1,4,5,7]:

$$arg\underset{\Delta \mathit{u}}{min}\left\{J_{\mathrm{MPC}}={(\overline{\mathit{y}}-\mathit{y})}^T·\mathit{\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[{\mathit{y}}^1,\dots ,{\mathit{y}}^{n_y}\right]}^T,{\mathit{y}}^j=\left[y_{k+1|k}^j,\dots ,y_{k+p|k}^j\right],$$

(5)

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

(6)

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

(7)

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

(8)

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

(9)

$y_{k+j|k}^j$ is a value of the jth output for the $(k+i)$th sampling instant predicted at the kth sampling instant using a control plant model, $\Delta u_{k+i|k}^m$ are future changes in manipulated variables, ${\kappa }_j\ge 0$ and ${\lambda }_m\ge 0$ are weighting coefficients for the predicted control errors of the jth output and for the changes of the mth manipulated variable, respectively; p and s denote prediction and control horizons, respectively; $n_y$, $n_u$ denote number of output and manipulated variables, respectively; $\Delta {\mathit{u}}_{min}$, $\Delta {\mathit{u}}_{max}$, ${\mathit{u}}_{min}$, ${\mathit{u}}_{max}$, ${\mathit{y}}_{min}$, ${\mathit{y}}_{max}$ are vectors of lower and upper bounds of changes and values of the control signals and of the values of output variables, respectively; ${\mathit{\kappa }}_i$ has p elements, ${\mathsf{\Lambda }}_i$ has s elements. The reference trajectories are grouped in the following vector:

$$\overline{\mathit{y}}={\left[{\overline{\mathit{y}}}^1,\dots ,{\overline{\mathit{y}}}^{n_y}\right]}^T,{\overline{\mathit{y}}}^j=\left[{\overline{y}}_{k+1|k}^j,\dots ,{\overline{y}}_{k+p|k}^j\right].$$

(10)

It is often assumed that reference trajectories are constant on the prediction horizon, then the vectors ${\overline{\mathit{y}}}^j$ (each composed of p elements) have the following form:

$${\overline{\mathit{y}}}^j=\left[{\overline{y}}_k^j,\dots ,{\overline{y}}_k^j\right],$$

(11)

where ${\overline{y}}_k^j$ is a setpoint value for the jth output.

As a solution to the optimization problem ((1)–(4)) the optimal vector of control action is obtained in such a way that future, predicted control errors are minimized (the first component of the performance index), control signals do not vary too much (the second component) and the constraints on control and output signals are fulfilled. The $\Delta u_{k|k}^m$ elements from the vector $\Delta \mathit{u}$ are applied in the control system and the algorithm passes to the next iteration. The predicted values of output signals $y_{k+j|k}^j$ are derived using a control plant model the MPC algorithm is based on.

2.1. MPC Algorithms Based on Nonlinear Models (NMPC)

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}),$$

(12)

where ${\mathit{y}}_{k-i}={\left[y_{k-i}^1,\dots ,y_{k-i}^{n_y}\right]}^T$ is a vector which contains output values measured at the $(k-i)$th sampling instant, ${\mathit{u}}_{k-i}={\left[u_{k-i}^1,\dots ,u_{k-i}^{n_u}\right]}^T$ is a vector which contains control values applied at the $(k-i)$th sampling instant; denote 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]$, $n_a,n_b$ determine how many past output and control values the model needs. If one wants to use the model (12) 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}{\mathit{y}}_{k+i|k}=\hfill & \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}$$

(13)

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 decision variables from the vector $\Delta \mathit{u}$ and

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

(14)

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 for example, Reference [7].

Unfortunately, the optimization task with the prediction (13), based on a nonlinear model, is not a quadratic optimization task. The optimization problem ((1)–(4)) becomes, in general, non–convex nonlinear optimization problem. Solving such a problem is difficult and time–consuming task, there is no guarantee to find a global solution and time needed to obtain solution cannot be estimated in advance. Different kinds of nonlinear process models can be used in the algorithm, but they are exploited in a similar way, see for example, References [13,40,41,42]; such algorithms will be referred to as Nonlinear MPC (NMPC). Due to problems caused by formulation of NMPC algorithms, in practice MPC algorithms based on linear approximations of process models, obtained at each iteration of the algorithm are considered. Such linear approximation is used to formulate the, easy to solve, quadratic optimization problem [35]. Such an optimization problem is also used in the MPC algorithms based on linear models, therefore now their formulation will be reminded. In the next section, the proposed algorithms, based on an advanced approximation of the nonlinear model, will be proposed.

2.2. MPC Algorithms Based on Linear Models (LMPC)

In the case of linear models 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},$$

(15)

$$\tilde{\mathit{y}}={\left[{\tilde{\mathit{y}}}^1,\dots ,{\tilde{\mathit{y}}}^{n_y}\right]}^T,{\tilde{\mathit{y}}}^j=\left[{\tilde{y}}_{k+1|k}^j,\dots ,{\tilde{y}}_{k+p|k}^j\right],$$

(16)

where $\tilde{\mathit{y}}$ is called the free response of the plant, because it contains future values of output signals calculated assuming that the control signal does not change in the future — it describes influence of the past values of control signals; form of the free response depends on the model used for prediction, see for example, References [1,4,5,7]. The second component in the formula (15) depends on future (to be calculated) changes of the control signals and is called the forced response; the matrix $\mathit{A}$, called the dynamic matrix, is given by the following formula:

$$\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],$$

(17)

$${\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],$$

(18)

where $a_i^{j,m}$ are the control plant step response coefficients describing influence of the mth input on the jth output; it can be shown that the dynamic matrix has the same form regardless a type of linear models MPC algorithms are based on [7].

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

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

(19)

After using prediction (15) also in the constraints (3), all constraints depend linearly on decision variables. Thus, the optimization problem ((1)–(4)) becomes an, easy to solve, standard linear–quadratic optimization problem. Unfortunately, application of an LMPC algorithm to a nonlinear process may result in unsatisfactory control performance, especially if operation in different operating points is demanded. Therefore, the fuzzy MPC (FMPC) algorithm being a combination of the LMPC and NMPC algorithms was proposed in Reference [36]. Now it will be improved by appropriate usage of the nonlinear model. Formulation of the proposed algorithm is presented in the next section. It still uses two models the original, nonlinear one is used to calculate the free response in an advanced way. The fuzzy one, composed of a set of a few step responses, is applied to calculate the dynamic matrix, which is updated at each iteration of the algorithm.

3. Fuzzy MPC Algorithm with Advanced Generation of Prediction

In order to avoid algorithm complications resulting from the need to solve a nonlinear optimization task in each iteration, one can use approximation of the process model carried out in each iteration. Using an easy–to–obtain fuzzy model to get the dynamic matrix and a nonlinear model to obtain the advanced free response leads to prediction which can be generated in a relatively simple way, is linear relative to decision variables, and mimics the nonlinearity of the process very well. Next, in each iteration of the algorithm, an optimization problem ((1)–(4)) which uses this prediction is formulated as a quadratic optimization problem. Such a problem is easy to solve using commonly available optimization routines. Now it will be described how to obtain such an advanced prediction.

3.1. Generation of the Dynamic Matrix

The dynamic matrix, needed to obtain forced response in the prediction can be calculated using an easy–to–obtain fuzzy Takagi–Sugeno model which has local models in the form of step responses [34]. Design process of such a model is very simple, because it is sufficient to collect a few sets of step responses (near a few operating points). These step responses can be obtained directly on the control plant or using a nonlinear model of the process. Next, the premises of the model should be formulated using expert knowledge, and then tuned, for example, by means of a fuzzy neural network [7]. Thus, the model 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}\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}\dots \mathrm{and} u_{k-m+1}^{j_u} \mathrm{is} C_m^{f,j_u}\hfill \\ \hfill & \mathrm{then} {\widehat{y}}_k^{j,f}={\displaystyle \sum _{m=1}^{n_u}}{\displaystyle \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}$$

(20)

where $y_k^{j_y}$ is the value of $j_y$th output variable at the kth sampling instant, $u_k^{j_u}$ is the value of $j_u$th manipulated variable at the kth 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$, l is number of fuzzy rules. The output values of the model (20) can be calculated using the following formula:

$${\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,$$

(21)

where ${\tilde{a}}_n^{j,m}={\sum }_{f=1}^l{\tilde{w}}_f·a_n^{j,m,f}$ and ${\tilde{w}}_f$ are the normalized firing strengths of the fuzzy rules calculated using fuzzy reasoning, see for example, References [43,44]. Note that one can consider the model (21) as the step response describing behavior of the control plant near the current operating point. This model is then used at each sampling instant of the proposed algorithms to obtain a new, updated dynamic matrix:

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

(22)

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

(23)

3.2. Generation of the Free Response

The basic idea of the free response generation for an MPC algorithm is to use a non–linear process model assuming that the control signals do not change over the entire prediction horizon—a classic free response; see for example, Reference [36]. In the case of the proposed algorithm, future control increments calculated in the algorithm in the previous iterations are also used when the free response is calculated. Thus, more information is exploited to obtain the improved prediction. 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,$$

(24)

where ${\check{u}}_{k+i|k}^m$ can be interpreted as a corrections of the control signal $u_{k+i|k-1}^m$ obtained in the last $(k-1)$st iteration. Similarly, future control increments will be given by the formula:

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

(25)

Assume

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

(26)

In this case, using the nonlinear model, the free response is iteratively derived as follows:

• First, the nonlinear model (12) is used to obtain

  $${\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});$$

  (27)
• Next, the values ${\widehat{\mathit{y}}}_{k+1}$ are used when calculating the output values ${\widehat{\mathit{y}}}_{k+2}$ for the next sampling instant:

  $${\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});$$

  (28)
• For the ith iteration, using the values ${\widehat{\mathit{y}}}_{k+1},\dots ,{\widehat{\mathit{y}}}_{k+i-1}$ one obtains:

  $${\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}).$$

  (29)
• After taking into account the estimation of unmeasured disturbances ${\mathit{d}}_k={\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-1}$, containing also influence of modeling errors, the final form of the formula describing the elements of the free response is obtained:

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

  (30)

3.3. Formulation of the Optimization Problem

Now the dynamic matrix (22) will be used to predict the influence of control corrections on the control plant outputs. Therefore, after using the free response (30) and the dynamic matrix (22) one obtains the following prediction:

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

(31)

Due to the change of decision variables, optimization task ((1)–(4)), which is solved by the control algorithm in each iteration can be written as:

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

(32)

subject to:

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

(33)

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

(34)

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

(35)

where

$$\Delta \check{\mathit{u}}={\left[\Delta {\check{\mathit{u}}}^1,\dots ,\Delta {\check{\mathit{u}}}^{n_u}\right]}^T,\Delta {\check{\mathit{u}}}^m=\left[\Delta {\check{u}}_{k+1|k}^m,\dots ,\Delta {\check{u}}_{k+s-1|k}^m\right].$$

(36)

Note that the performance function in (32) 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 was obtained.

Remark 1.

In the proposed algorithm, one can also use a slightly modified form of the performance index, in which the second component depends only on the corrections $\Delta \check{\mathit{u}}$ of the control signals determined in the previous iteration of the algorithm, that is,

$$J_{MPC+v2}={(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{A}}_k·\Delta \check{\mathit{u}})}^T·\mathit{\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}}.$$

(37)

It should be emphasized, however, that such modification causes a change in the meaning of the λ parameter in relation to the classical formulation. Consequently, algorithms with a modified form of the performance function generate faster responses, but influence of the λ parameter on system robustness becomes significantly limited.

3.4. Iterative Improvement of the Prediction

In practice, the proposed approach gives results very similar to those offered by the NMPC algorithm. However it is possible to further improve the control quality, if necessary, thanks to iterative improvement of the prediction determined in the algorithm. Then smaller corrections of control trajectories in each subsequent iteration are obtained, what is beneficial because impact of them on the process outputs is predicted using the dynamic matrix instead of a more accurate non–linear process model.

The first, internal iteration of the algorithm is the same as in the basic version of the algorithm. After receiving the future control trajectories, one can either apply them in the process or try to improve the resulting solution. In such a case, the second internal iteration of the algorithm should be executed. Denote the values of future control signal trajectories obtained in the first iteration as $u_{k+i\left|k\right(1)}^m$ and future control increments as $\Delta u_{k+i\left|k\right(1)}^m$.

In subsequent iterations, because subsequent control signal trajectories are available, future control values are decomposed as follows:

$$u_{k+i\left|k\right(i_w)}^m={\check{u}}_{k+i\left|k\right(i_w)}^m+u_{k+i\left|k\right(i_w-1)}^m,$$

(38)

where $i_w$—index of the internal iteration of the algorithm, $i_w=2,3,\dots$. Similarly, future control increments are decomposed:

$$\Delta u_{k+i\left|k\right(i_w)}^m=\Delta {\check{u}}_{k+i\left|k\right(i_w)}^m+\Delta u_{k+i\left|k\right(i_w-1)}^m.$$

(39)

Assume

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

(40)

Thus, using the non–linear process model, the free response is iteratively derived as follows:

• First, the model is used to obtain values

  $${\widehat{\mathit{y}}}_{k+1}=\mathit{f}({\mathit{y}}_k,{\mathit{y}}_{k-1},\dots ,{\mathit{y}}_{k-n_a+1},{\mathit{u}}_{k\left|k\right(i_w-1)},{\mathit{u}}_{k-1},\dots ,{\mathit{u}}_{k-n_b+1});$$

  (41)
• Next, the values ${\widehat{\mathit{y}}}_{k+1}$ are used to calculate the output values ${\widehat{\mathit{y}}}_{k+2}$ for the next sampling instant:

  $${\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\left|k\right(i_w-1)},{\mathit{u}}_{k\left|k\right(i_w-1)},{\mathit{u}}_{k-1},\dots ,{\mathit{u}}_{k-n_b+2});$$

  (42)
• For the ith iteration, using the values ${\widehat{\mathit{y}}}_{k+1},\dots ,{\widehat{\mathit{y}}}_{k+i-1}$ one obtains:

  $$\begin{array}{cc}{\widehat{\mathit{y}}}_{k+i}=\mathit{f}(\hfill & {\widehat{\mathit{y}}}_{k+i-1},{\widehat{\mathit{y}}}_{k+i-2},\dots ,{\mathit{y}}_{k-n_a+i},\hfill \\ & {\mathit{u}}_{k+i-1\left|k\right(i_w-1)},{\mathit{u}}_{k+i-2\left|k\right(i_w-1)},\dots ,{\mathit{u}}_{k-n_b+i}).\hfill \end{array}$$

  (43)
• After taking into account the estimation of unmeasured disturbances ${\mathit{d}}_k={\mathit{y}}_k-{\widehat{\mathit{y}}}_{k|k-1}$, containing also influence of modeling errors, like in (30), the final formula describing the elements of the free response is obtained:

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

  (44)

As in the basic algorithm, the dynamic matrix is used to predict the influence of control corrections on the control plant outputs. Thus, after using the free response (44) and the dynamic matrix (22) one obtains the following prediction:

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

(45)

where

$$\Delta \check{\mathit{u}}={\left[\Delta {\check{\mathit{u}}}^1,\dots ,\Delta {\check{\mathit{u}}}^{n_u}\right]}^T,\Delta {\check{\mathit{u}}}^m=\left[\Delta {\check{u}}_{k+1\left|k\right(i_w)}^m,\dots ,\Delta {\check{u}}_{k+s-1\left|k\right(i_w)}^m\right].$$

(46)

This prediction should then be used to formulate the optimization task in the same way as in the basic algorithm. After solving this task, one can either use the result to control the process or to perform another internal iteration of the algorithm.

Note that the lower the values of the corrections of the future control trajectories are obtained, the better the prediction is obtained, because it will be closer to the one obtained fully using a nonlinear model. The number of internal iterations of the algorithm can therefore be selected depending on how large the corrections of the future control trajectories are (provided that the time for generating control signals has not expired). If these adjustments become small enough, one should stop doing further internal iterations.

Due to the lack of proof of the convergence of the described procedure, it is advisable to check after each internal iteration whether the performance function in the problem solved by the algorithm and control corrections decrease. If there is no improvement in the value of the performance function (especially when there is a large increase of it) or control adjustments begin to increase the internal iterations should be interrupted.

3.5. Fast Generation of the Control Action—Analytical Approach

Note that if in each iteration of the proposed approach the following optimization problem without constraints is solved:

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

then it has the following solution given by the analytical formula:

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

(47)

where $\Delta \overrightarrow{\mathit{u}}$ contains elements of future control increments trajectory, corrections of which are collected in the vector $\Delta \check{\mathit{u}}$.

Remark 2.

If one uses the performance index (37), that is, the modified version:

$$J_{MPC+v2}={(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{A}}_k·\Delta \check{\mathit{u}})}^T·\mathit{\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}},$$

(48)

then the following analytical formula should be used:

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

(49)

It should be once again emphasized, however, that this modification changes the meaning of the parameter matrix Λ comparing to the classical formulation and algorithms with a modified form of the performance index generate faster responses than the algorithm using the control law (47), but influence of the λ parameter on system robustness becomes significantly limited.

Due to the fact that the dynamic matrix ${\mathit{A}}_k$ changes, in general, in every iteration, then in each iteration of the algorithm the formulas (47) (or (49)) allowing to calculate the control action should be used. As a consequence, appropriate rows of the matrix ${\left({\mathit{A}}_k^T·\mathit{\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 in each iteration of the algorithm.

3.6. Disturbance Measurement Utilization

The proposed algorithms are formulated in such a way that the mechanisms of disturbance measurement utilization presented in Reference [45] can be applied in them easily. The proposed algorithms are based on two models, therefore two such mechanisms can be applied. They will be now shortly reminded.

3.6.1. Employing Fuzzy Model

The first method consists in extending local models of the fuzzy model (20) using the models of disturbance influence on the control plant in the form of step responses. It results then in extending the prediction (31) with a term dependent on measured disturbance; for details see References [34,45]:

$$\mathit{y}=\tilde{\mathit{y}}+{\mathit{E}}_k·\Delta {\mathit{u}}^d+{\mathit{A}}_k·\Delta \check{\mathit{u}},$$

(50)

where

$${\mathit{E}}_k=\left[\begin{array}{c}{\mathit{E}}_k^1\\ {\mathit{E}}_k^2\\ ⋮\\ {\mathit{E}}_k^{n_y}\end{array}\right],{\mathit{E}}_k^j=\left[\begin{array}{cccccc}{\tilde{e}}_1^j& {\tilde{e}}_2^j-{\tilde{e}}_1^j& {\tilde{e}}_3^j-{\tilde{e}}_2^j& \dots & {\tilde{e}}_{p_z-1}^j-{\tilde{e}}_{p_z-2}^j& {\tilde{e}}_{p_z}^j-{\tilde{e}}_{p_z-1}^j\\ {\tilde{e}}_2^j& {\tilde{e}}_3^j-{\tilde{e}}_1^j& {\tilde{e}}_4^j-{\tilde{e}}_2^j& \dots & {\tilde{e}}_{p_z}^j-{\tilde{e}}_{p_z-2}^j& {\tilde{e}}_{p_z}^j-{\tilde{e}}_{p_z-1}^j\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ {\tilde{e}}_p^j& {\tilde{e}}_{p+1}^j-{\tilde{e}}_1^j& {\tilde{e}}_{p+2}^j-{\tilde{e}}_2^j& \dots & {\tilde{e}}_{p_z}^j-{\tilde{e}}_{p_z-2}^j& {\tilde{e}}_{p_z}^j-{\tilde{e}}_{p_z-1}^j\end{array}\right],$$

and

$$\Delta {\mathit{u}}^d={\left[\Delta u_k^d,\Delta u_{k-1}^d,\Delta u_{k-2}^d,\dots ,\Delta u_{k-p_z+2}^d,\Delta u_{k-p_z+1}^d\right]}^T,$$

where ${\tilde{e}}_n^j={\sum }_{f=1}^l{\tilde{w}}_f·e_i^{j,f}$, $e_i^{j,f}$ ($i=1,\dots ,p_z$, $j=1,\dots ,n_y$, $f=1,\dots ,l$) are the coefficients of step responses to the change of disturbance, in the fth local model, $u_{k-i}^d$ are values of the disturbance measured in the $(k-i)$th sampling instant.

Then the following performance function in the optimization problem (1)–(4) should be now used:

$$\begin{array}{cc}\hfill J_{MPC+}=& {(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{E}}_k·\Delta {\mathit{u}}^d-{\mathit{A}}_k·\Delta \check{\mathit{u}})}^T·\mathit{\kappa }·(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{E}}_k·\Delta {\mathit{u}}^d-{\mathit{A}}_k·\Delta \check{\mathit{u}})\hfill \\ & +\Delta {\mathit{u}}^T·\mathsf{\Lambda }·\Delta \mathit{u}.\hfill \end{array}$$

(51)

or the modified one:

$$\begin{array}{cc}{J}_{MPC+v2}=\hfill & {(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{E}}_k·\Delta {\mathit{u}}^d-{\mathit{A}}_k·\Delta \check{\mathit{u}})}^T·\mathit{\kappa }·(\overline{\mathit{y}}-\tilde{\mathit{y}}-{\mathit{E}}_k·\Delta {\mathit{u}}^d-{\mathit{A}}_k·\Delta \check{\mathit{u}})\hfill \\ & +\Delta {\check{\mathit{u}}}^T·\mathsf{\Lambda }·\Delta \check{\mathit{u}}.\hfill \end{array}$$

(52)

In the example presented in the next section, this mechanism was simplified by assuming that only one model of disturbance influence on the process is obtained, that is, $e_i^{j,f}=e_i^j$ for $f=1,\dots ,l$, then the prediction transforms to the following formula:

$$\mathit{y}=\tilde{\mathit{y}}+\mathit{E}·\Delta {\mathit{u}}^d+{\mathit{A}}_k·\Delta \check{\mathit{u}},$$

(53)

where

$$\mathit{E}=\left[\begin{array}{c}{\mathit{E}}^1\\ {\mathit{E}}^2\\ ⋮\\ {\mathit{E}}^{n_y}\end{array}\right],{\mathit{E}}^j=\left[\begin{array}{cccccc}{e}_1^j& e_2^j-e_1^j& e_3^j-e_2^j& \dots & e_{p_z-1}^j-e_{p_z-2}^j& e_{p_z}^j-e_{p_z-1}^j\\ e_2^j& e_3^j-e_1^j& e_4^j-e_2^j& \dots & e_{p_z}^j-e_{p_z-2}^j& e_{p_z}^j-e_{p_z-1}^j\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ e_p^j& e_{p+1}^j-e_1^j& e_{p+2}^j-e_2^j& \dots & e_{p_z}^j-e_{p_z-2}^j& e_{p_z}^j-e_{p_z-1}^j\end{array}\right].$$

and the matrix $\mathit{E}$ is constant in each iteration, what simplifies the calculations. Thus, disturbance measurement utilization mechanism becomes actually the same as in the DMC algorithm.

3.6.2. Employing Nonlinear Model

The second approach assumes that the disturbance is one of the inputs of the model used to obtain the free response:

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

(54)

where $u_{k-i}^d$ is value of disturbance in the $(k-i)$th sampling instant. Then the nonlinear process model (12) becomes:

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

(55)

Then it is sufficient to use measured values of the disturbance during calculation of the free response to include this information in the prediction. During calculation of the free response, as future values of the disturbance its estimates can be used or, if one does not have any, they should be assumed unchanged, like in the case of classical free response.

4. Example

The control plant is a nonlinear, isothermal CSTR in which the van de Vusse reaction takes place. The steady–state characteristic of the control plant is depicted in Figure 1.

The CSTRs in which the van de Vusse reaction takes place are a popular benchmark due to their nonlinearity and difficult dynamics, willingly used to test newly developed control algorithms; see for example, References [46,47,48,49,50,51]. The process model of the reactor is composed of two composition balance equations; see for example, Reference [52]:

$$\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}$$

(56)

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}$), $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 parameters are: $k_1=50 1/\mathrm{h}$, $k_2=100 1/\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 control plant has the inverse response, thus it is a control plant difficult to control using a standard algorithm. Therefore, it is natural to use an MPC algorithm in this case.

Experiments

For the considered control plant a few MPC algorithms were designed—an NMPC one based on the nonlinear model and nonlinear optimization, an LMPC one based on a linear model in the form of a step response (of DMC type) and the proposed FMPC algorithms using a fuzzy model, different methods of prediction improvements, and different forms of the performance index (FMPC1—with classical performance index, and FMPC2—with the modified performance index). The fuzzy model used in the FMPC algorithms is composed of three step responses obtained near the following operating points:

$$\begin{gathered} \mathrm{R}1 C_{B0}=0.91 \mathrm{mol}/\mathrm{L}, C_{A0}=2.18 \mathrm{mol}/\mathrm{L}, F=20 \mathrm{L}/\mathrm{h}; \\ \mathrm{R}2 C_{B0}=1.12 \mathrm{mol}/\mathrm{L}, C_{A0}=3 \mathrm{mol}/\mathrm{L}, F=34.3 \mathrm{L}/\mathrm{h}; \\ \mathrm{R}3 C_{B0}=1.22 \mathrm{mol}/\mathrm{L}, C_{A0}=3.66 \mathrm{mol}/\mathrm{L}, F=50 \mathrm{L}/\mathrm{h}. \end{gathered}$$

The model used to obtain the dynamic matrix is thus composed of three rules; the membership functions, same as in Reference [34], chosen after analysis of the steady–state characteristic of the control plant, are shown in Figure 2.

The simulation experiments were done using Matlab. The sampling time equal to $T_{\mathrm{s}}=3.6$ s was assumed; the values of tuning parameters were as follows: prediction horizon $p=70$, control horizon $s=35$, and at the beginning weighting coefficient $\lambda =0.001$. During the experiments operation of control systems with NMPC, LMPC and FMPC algorithms was compared. The LMPC algorithm was designed using the step response obtained near the R2 point. The FMPC algorithms use two models: nonlinear one in the form of state equations, to generate a free response and the fuzzy one with local models in the form of step responses, to obtain the dynamic matrix.

Responses of the control system 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. If the setpoint was changed to ${\overline{C}}_{B1}=1$ mol/L, responses obtained in the control system with the FMPC1 algorithm (blue lines in Figure 3) are very close to those obtained with NMPC algorithm with nonlinear optimization (red lines in Figure 3). In both cases there is almost no overshoot and both algorithms work much better than their LMPC counterpart (magenta lines in Figure 3). In the latter case there is a significant overshoot. Moreover, settling time is much longer than in the case of NMPC and FMPC algorithms. The FMPC1 algorithm is slightly slower than the NMPC one, whereas the FMPC2 is practically as fast as NMPC.

In the case when setpoint changed to ${\overline{C}}_{B2}=1.25$ mol/L, the responses obtained with the FMPC algorithms are even faster than the one generated with NMPC algorithm. In the case of all the algorithms the setpoint is reached without overshoot, but LMPC works much slower than the algorithms based on nonlinear models.

Disturbance responses obtained near $C_B=1$ mol/L with FMPC and NMPC algorithms are almost the same. There is no overshoot which is present in the response of the LMPC algorithm. The algorithms based on nonlinear models slightly faster compensate the change in disturbance. In the case of operation near $C_B=1.25$ mol/L, the NMPC is the fastest in compensation of the disturbance change. The FMPC1 is slightly slower and the LMPC is much slower than its counterparts.

It can be noticed that the FMPC algorithms offer control performance close to the one which can be obtained with the NMPC algorithm and with nonlinear optimization. In one case (setpoint change to ${\overline{C}}_{B2}=1.25$ mol/L) they are even faster. It should be, however, emphasized that the FMPC algorithms use reliable quadratic programming routine to generate control action.

Duration of action of the researched algorithms in the experiments depicted in Figure 3 is compared in

Table 1. Average duration of action of the algorithms in the experiments depicted in Figure 3.

${\overline{\mathit{C}}}_{\mathit{B}}$	NMPC	FMPC1	FMPC2
$1.00$ mol/L	29.1671	1.8175	1.8118
$1.25$ mol/L	47.0518	1.8116	1.8137
sum	76.2189	3.6291	3.6254

. For each algorithm and for each setpoint five measurements of duration were done (time needed to simulate the control plant was subtracted) and averaged. Then these average durations were summed in the last row of the table. Both FDMC algorithms work much (around 20 times) faster than the NMPC algorithm.

There was also done an experiment with $\lambda =0.0001$ (10 times smaller than in the first experiment; Figure 4). All the algorithms work now faster; the control action in all cases is more aggressive, but the relations between FMPC, NMPC, and LMPC algorithms remain unchanged; they are the same as in the previous experiment, for $\lambda =0.001$. However, as the $\lambda$ parameter is closer to zero, the differences between FMPC1 and FMPC2 algorithms became much smaller comparing to the previous experiment. The FMPC1 algorithm is practically as fast as the FMPC2 one. Such a phenomenon is reasonable, because for $\lambda =0$ both performance indexes (32) and (37) are the same. Thus, the closer value of the $\lambda$ coefficient to 0 is, the closer responses generated with FMPC1 and FMPC2 are.

There were also done experiments with shortening the control horizon s (Figure 5 and Figure 6). Note, that the smaller the control horizon is the less number of decision variables in the optimization problem. Therefore, it is simpler to solve and can be solved faster. In the case of both algorithms FMPC1 and FMPC2, with different values of $\lambda$ for many values of the control horizon the responses generated with the algorithms do not change too much. Thus, the control horizon can be shorted significantly. In the considered examples for FMPC1 algorithm the responses generated for $s=35$ and for $s=10$ are almost the same; in the case of control system with FMPC2 algorithm—the responses generated for $s=35$ and for $s=5$ are almost the same. Further decreasing of the control horizon brings changes in the responses.

In the case of the control system with FMDC1 algorithm, for $s=5$ the disturbance compensation near $C_B=1.25$ mol/L becomes faster, but unfortunately achieving the setpoint ${\overline{C}}_{B1}=1$ mol/L lasts longer. Further decrease of the control horizon to $s=1$ causes that the disturbance compensation near $C_B=1.25$ mol/L becomes even faster as well as reaching the setpoint ${\overline{C}}_{B2}=1.25$ mol/L, however at the cost of slower response to the setpoint change to ${\overline{C}}_{B1}=1$ mol/L and worsened disturbance compensation near $C_B=1$ mol/L.

In the control system with FMDC2 algorithm, for $s=2$ and $s=1$ the disturbance compensation near $C_B=1.25$ mol/L becomes slightly faster, but unfortunately achieving both the setpoints ${\overline{C}}_{B1}=1$ mol/L and ${\overline{C}}_{B2}=1.25$ mol/L lasts longer. Moreover disturbance compensation near $C_B=1$ mol/L is now slower. Thus, in this case, as in the case of FMPC1 algorithm, decreasing the control horizon should be done also carefully and there is also a threshold which should not be crossed in order to avoid control performance deterioration.

To sum up, during designing of an algorithm it is advisable to check if the control horizon can be decreased (and optimization problem solved by the algorithm simplified). However, it should be done carefully, because for too short control horizon control performance can got worse.

The mechanisms of disturbance measurement utilization were used in the next experiments. Both FMPC algorithms generated similar responses. The first method of disturbance measurement utilization (red lines in Figure 7 and Figure 8) is based on the same approach as in the standard DMC algorithm, the second one employs the nonlinear process model (magenta lines in Figure 7 and Figure 8). Disturbance compensation near $C_B=1.25$ mol/L is almost the same in the case of both methods; differences in responses are really small. It is because the model of disturbance influence on the control plant for the first method was obtained near this output value.

Application of both methods of disturbance measurement utilization near $C_B=1$ mol/L causes reduction of the maximal control error by around 50% comparing to the case when disturbance measurement was not exploited (blue lines in Figure 7 and Figure 8). Both methods offer very good disturbance compensation. However, in the case of the one based on the nonlinear model (magenta lines in Figure 7 and Figure 8) the output value achieves the setpoint slightly faster than the one based on the linear model (red lines in Figure 7 and Figure 8) while changes of the control signal are smaller.

5. Conclusions

The numerically efficient FMPC algorithms proposed in the article use the easy–to–obtain fuzzy model to calculate the dynamic matrix, which describes the influence of future control action on prediction, and the nonlinear model to derive the free response of the control plant in a clever way which allows obtaining very good prediction. As a result, the algorithms offer control performance very close to the one offered by the NMPC algorithm, based on the nonlinear optimization. The FMPC algorithms, however, despite they are also based on nonlinear models, thanks to a clever way of prediction generation are formulated as the easy–to–solve quadratic optimization problem.

Modifications of the algorithms proposed to improve their operation can be applied by a system designer if needed. Thus the algorithm can be tailored to a particular application. Thus, it can offer a good balance between control performance which is needed in the problem at hand, and control algorithm complication. One can adjust the algorithm to a given problem so the computational demand is as small possible.