Tuning of Multivariable Model Predictive Control for Industrial Tasks

Abstract

This work is concerned with the tuning of the parameters of Model Predictive Control (MPC) algorithms when used for industrial tasks, i.e., compensation of disturbances that affect the process (process uncontrolled inputs and measurement noises). The discussed simulation optimisation tuning procedure is quite computationally simple since the consecutive parameters are optimised separately, and it requires only a very limited number of simulations. It makes it possible to perform a multicriteria control assessment as a few control quality measures may be taken into account. The effectiveness of the tuning method is demonstrated for a multivariable distillation column. Two cases are considered: a perfect model case and a more practical case in which the model is characterised by some error. It is shown that the discussed tuning approach makes it possible to obtain very good control quality, much better than in the most common case in which all tuning parameters are constant.

1. Introduction

There are a few advanced control methods, including model reference adaptive control [1], fault-tolerant control [2], stochastic control [3], fuzzy control [4] and Model Predictive Control (MPC) [5,6]. In particular, MPC algorithms make it possible to obtain excellent control quality in the case of Multiple-Input Multiple-Output (MIMO) processes with constraints. As a result, MPC methods have been used to numerous industrial processes [7], e.g., chemical reactors [8], distillation columns [9], waste water treatment plants [10], solar power stations [11], cement kilns [12], pasteurization plants [13] and pulp digesters [14]. In addition to that, MPC algorithms are more and more popular in other areas; example applications are: fuel cells [15], active vibration attenuation [16], combustion engines [17], robots [18], synchronous motors [19], mechanical systems [20], freeway traffic congestion control [21] and autonomous driving [22].

Tuning of MPC is an important issue since a good choice of parameters is likely to significantly increase control quality. The prediction and control horizons may be determined taking into account the process dynamics, the possible sampling period of the controller and the computational performance of the hardware platform used [23,24]. Additionally, it is possible to use the automatic tuning methodology [25] in which a step response model is experimentally obtained from the process, and the horizons are chosen using some hands-on tuning rules. A thorough review on how to find proper MPC horizons was given in [26]. On the other hand, the choice of numerous coefficients of the minimised MPC cost function is always an issue. This task is particularly important in the case of MIMO processes with strong interactions between the consecutive inputs and outputs. A few tuning methods have been discussed in the literature. In some simple cases, it is possible to explicitly calculate the tuning coefficients [27,28]. A tuning method based on the Robust Performance Number (RPN) was described in [29]. It is also possible to dynamically calculate set-points for MPC to accommodate user-defined output importance [30]. Multi-objective performance optimisation using the goal attainment approach was considered in [31]. A thorough comparison of a few heuristic optimisation algorithms was reported in [32] (the Particle Swarm Optimisation (PSO) method, the firefly algorithm, the grey wolf optimiser and the Jaya algorithm were used). An application of the genetic algorithm was described in [12], whereas the PSO algorithm used to find the parameters of MPC with model uncertainty was considered in [33]. In general, the optimisation-based tuning methods are very computationally demanding as there are numerous decision variables, which means that a large number of simulations must be performed. As a practical solution, a relatively computationally simple simulation optimisation tuning method presented in [34] may be used. It needs only a very limited number of simulations. The considered tuning method was discussed in [35] for the MPC employed for vehicle obstacle avoidance.

The tuning method discussed in [34] typically uses step set-point changes. In this work, a more industrially practical scenario is considered, i.e., compensation of disturbances that affect the process. The disturbances considered include process uncontrolled inputs and measurement noises. For a MIMO distillation column, the tuning procedure is detailed. Two cases are considered: a perfect model case and a more practical case in which the model is characterised by some error. It is shown that the discussed tuning approach makes it possible to use very good control quality, much better than in the most common case in which all tuning parameters are constant. The multicriteria control assessment is used since the control quality is assessed taking into account three factors: the sum of squared errors, the Huber standard deviation and the entropy [36,37].

The article is structured in the following way. First, in Section 2, the MPC task and its tuning parameters are recalled. Next, in Section 3, the tuning procedure is detailed, and the indicators used for control quality assessment are defined. Section 4 reports an application of the considered method for a MIMO distillation process with industrial disturbances. Finally, Section 5 summarises the article.

2. MPC Problem Formulation

In this work, we deal with MIMO processes that have as many as $n_{\mathrm{u}}$ manipulated variables (inputs) and $n_{\mathrm{y}}$ controlled variables (outputs). Hence, the following vectors are used: $u={\left[u_1\dots u_{n_{\mathrm{u}}}\right]}^{\mathrm{T}}$ and $y={\left[y_1\dots y_{n_{\mathrm{y}}}\right]}^{\mathrm{T}}$. At each discrete sampling instant, $k=0,1,\dots$, the MPC algorithm calculates on-line the vector of decision variables, which is defined by the increments of the future manipulated variables [6]:

$$▵\mathit{u}\left(k\right)=\left[\begin{array}{c}▵u\left(k\right|k)\\ ⋮\\ ▵u(k+N_{\mathrm{u}}-1|k)\end{array}\right]$$

(1)

where $N_{\mathrm{u}}$ is the control horizon. The MPC decision variables (1) are found as a result of a computational procedure in which the predicted control quality is optimised. An example of the minimised MPC cost function is:

$$J\left(k\right)=\sum _{p=1}^N\sum _{n=1}^{n_{\mathrm{y}}}{\psi }_{p,n}{\left(y_n^{\mathrm{sp}}(k+p|k)-{\widehat{y}}_n(k+p|k)\right)}^2+\sum _{p=0}^{N_{\mathrm{u}}-1}\sum _{n=1}^{n_{\mathrm{u}}}{\lambda }_{p,n}{\left(▵u_n(k+p|k)\right)}^2$$

(2)

The first part of the cost function measures the control errors predicted over the prediction horizon N. The set-points and predicted values for the future sampling instant $k+p$ known at the sampling instant k are denoted by $y_n^{\mathrm{sp}}(k+p|k)$ and ${\widehat{y}}_n(k+p|k)$, respectively, for all process outputs, i.e., $n=1,\dots ,n_{\mathrm{y}}$. The role of the second part of the cost function is to minimise unwanted big changes of the manipulated variables. In some applications for the cost function (2), the predicted squared values of the future manipulated variables may be also taken into account. The tuning coefficients are: ${\psi }_{p,n}\ge 0$ for $p=1,\dots ,N$, $n=1,\dots ,n_{\mathrm{y}}$ and ${\lambda }_{p,n}>0$ for $p=0,\dots ,N_{\mathrm{u}}-1$, $n=1,\dots ,n_{\mathrm{u}}$. Typically, the penalties ${\lambda }_{p,n}$ are chosen in such a way that the manipulated variables do not change rapidly. As far as the coefficients ${\psi }_{p,n}$ are concerned, the choice is more difficult. This is because these coefficients prioritise the predicted control errors of the consecutive controlled variables. Additionally, predicted control errors for different sampling instants over the prediction horizon may be taken into account in a different way. In practice, however, for simplicity, all coefficients ${\psi }_{p,n}$ are set to a constant value. Unfortunately, as is shown in Section 4, this may deteriorate the resulting control quality. The cost function (2) may be conveniently rewritten in a compact form:

$$\begin{array}{cc}\hfill J\left(k\right)& =\sum _{p=1}^N{∥y^{\mathrm{sp}}(k+p|k)-\widehat{y}(k+p|k)∥}_{{\mathbf{\Psi }}_p}^2+\sum _{p=0}^{N_{\mathrm{u}}-1}{∥▵u(k+p|k)∥}_{{\mathbf{\Lambda }}_p}^2\hfill \end{array}$$

(3)

where the set-point vector for the sampling instant $k+p$ is denoted by $y^{\mathrm{sp}}(k+p|k)$, the predicted vector of the output variables for the sampling instant $k+p$ is denoted by $\widehat{y}(k+p|k)$ and both vectors are of length $n_{\mathrm{y}}$. The semidefinite-positive matrix ${\mathbf{\Psi }}_p=\mathrm{diag}({\psi }_{p,1},\dots ,{\psi }_{p,n_{\mathrm{y}}})$ is of dimensionality $n_{\mathrm{y}}\times n_{\mathrm{y}}$, and the definite-positive matrix ${\mathbf{\Lambda }}_p=\mathrm{diag}({\lambda }_{p,1},\dots ,{\lambda }_{p,n_{\mathrm{u}}})$ is of dimensionality $n_{\mathrm{u}}\times n_{\mathrm{u}}$. In practice, we usually consider some constraints put on process variables. The typical rudimentary MPC optimisation problem is:

$$\begin{array}{cc}\hfill & \underset{▵\mathit{u}\left(k\right)}{min}\left\{J\left(k\right)={\displaystyle \sum _{p=1}^N}{∥y^{\mathrm{sp}}(k+p|k)-\widehat{y}(k+p|k)∥}_{{\mathbf{\Psi }}_p}^2+{\displaystyle \sum _{p=0}^{N_{\mathrm{u}}-1}}{∥▵u(k+p|k)∥}_{{\mathbf{\Lambda }}_p}^2\right\}\hfill \\ \hfill & \mathrm{subject}\mathrm{to}:\hfill \\ \hfill & u^{min}\le u(k+p|k)\le u^{max},p=0,\dots ,N_{\mathrm{u}}-1\hfill \\ \hfill & ▵u^{min}\le ▵u(k+p|k)\le ▵u^{max},p=0,\dots ,N_{\mathrm{u}}-1\hfill \\ \hfill & y^{min}\le \widehat{y}(k+p|k)\le y^{max},p=1,\dots ,N\hfill \end{array}$$

(4)

where the vectors $u^{min}$ and $u^{max}$ define the minimal and maximal values of the manipulated variables, respectively, the vectors $▵u^{min}$ and $▵u^{max}$ define the minimal and maximal changes of the manipulated variables, respectively, and the vectors $y^{min}$ and $y^{max}$ define the minimal and maximal predicted values of the controlled variables, respectively. In spite of the fact that as many as $n_{\mathrm{u}}{N}_{\mathrm{u}}$ decision variables (1) are computed at every sampling instant k, only the first $n_{\mathrm{u}}$ of them are actually applied to the process; in the consecutive sampling instants, the whole calculation procedure is repeated.

In this work, the Dynamic Matrix Control (DMC) algorithm [38] is considered. A characteristic feature of the DMC algorithm is the fact that for prediction, i.e., to calculate the scalars ${\widehat{y}}_n(k+p|k)$ or the vectors $\widehat{y}(k+p|k)$, a discrete-time step-response model of the controlled process is used. The model may be obtained in a very simple way from the real process; no complicated identification algorithms need be used, which is a huge advantage. Since the step-response model is linear in terms of the manipulated variables, minimisation of the MPC optimisation task (4) is, in fact, a quadratic optimisation problem. Details related to the implementation of the DMC algorithm for MIMO processes may be found in [6].

3. MPC Tuning Procedure

When all parameters ${\psi }_{p,n}$ for $p=1,\dots ,N$, $n=1,\dots ,n_{\mathrm{y}}$ are optimised at the same time, the optimisation problem has as many as $n_{\mathrm{y}}N$ decision variables [12,32,33]. For typical lengths of the prediction horizon and a few controlled variables, the resulting optimisation problem is quite difficult. In order to simplify the computational task, the discussed tuning method is based on the following two rules:

• The tuning coefficients for only one manipulated variable are calculated at the same time, i.e., ${\psi }_{1,n},\dots ,{\psi }_{N,n}$, for the consecutive process outputs $n=1,\dots ,n_{\mathrm{y}}$.
• If all the coefficients ${\psi }_{1,n},\dots ,{\psi }_{N,n}$ were optimised directly at the same time by a numerical optimisation procedure, we still would have as many as N decision variables. As an alternative, that sequence of coefficients is parameterised using Gauss-like functions [34]. At first, a very general approximation of the Gauss function is used in which for one point, the function has a value of $K_n$; for all other points, it has a default value equal to one:

  $${\psi }_{p,n}=\left\{\begin{array}{cc}{K}_n\hfill & \mathrm{if}p=m_n\hfill \\ 1\hfill & \mathrm{if}p\ne m_n\hfill \end{array}\right.$$

  (5)
  The function (5) is characterised by two parameters: $m_n$ and $K_n$. The first one defines the chosen sampling instant within the prediction horizon for which the function has the value of $K_n$. When the parameters $m_n$ and $K_n$ are selected, the trajectory of the weights is calculated precisely from the Gauss function:

  $${\psi }_{p,n}=K_{n}exp\left(-{\left(\frac{p-m_n}{a_n}\right)}^2\right)$$

  (6)

  where $a_n$ defines the spread.

All things considered, irrespective of the number of process outputs and the prediction horizon, in the discussed approach, there are always only three decision variables: $m_n$, $K_n$ and $a_n$. What is important is that they are not calculated at the same time, but always, only one parameter is found. Calculations are repeated for all process outputs, $n=1,\dots ,n_{\mathrm{y}}$. Unlike fully-fledged optimisation-based tuning methods, in our simulation-based approach, we simply assess how the consecutive parameters influence the value of the performance indices, and we choose the parameters for which the best results are possible. No numerical optimisation is used, but we simply evaluate the values of control the performance indices for a chosen set of tuning parameters. Of course, our procedure is suboptimal, i.e., one may expect that full numerical optimisation of all tuning parameters at the same time may give better results, but our procedure is deliberately rather uncomplicated and is recommended to be used in industrial applications in which all tuning parameters are usually constant.

The tuning procedure is comprised of the following steps:

• The trajectory of the N coefficients ${\psi }_{1,n},\dots ,{\psi }_{N,n}$ is found for the controlled variable n. All other coefficients are tuned or have their default value.

  (a)

  The best trajectory (5) is found, i.e., its parameters $m_n$ and $K_n$ are determined.

  • The constant value of the parameter $K_n$ is assumed. Its value has to be larger than one, and it should result in a noticeable change in control quality in comparison with the control quality achieved for the scenario in which all coefficients ${\psi }_{p,n}=1$.
  • The performance indices are calculated from simulations of the MPC algorithm for a few values of the parameter $m_n$. It is recommended to start the tests from $m_n=N/2$ and analyse the results obtained in some neighbourhood of this value first.
  • The value of the parameter $m_n$ is chosen for which the performance indices are the best.
  • For the chosen value of $m_n$, the performance indices are calculated from the simulations of the MPC algorithm for a few values of the parameter $K_n$. This step should start with the values of the parameter that are lower than those assumed in the first step. Next, it is increased. It is recommended not to choose too large values of $K_n$ as this may result in dangerously large values of the manipulated variables.
  • The value of the parameter $K_n$ is chosen for which the performance indices are the best.

  (b)

  The best trajectory (6) is found, i.e., its parameter $a_n$ is determined.

  • The selected parameters $m_n$ and $K_n$ are used.
  • The performance indices are calculated from the simulations of the MPC algorithm for a few values of the parameter $a_n$. It is recommended to perform simulations starting from relatively low values of the parameter, e.g., $a_n=1$, and slowly increase it until it is clear that any further increment results in no improvement of control quality.
  • The value of the parameter $a_n$ is chosen for which the performance indices are the best.

• Tuning is repeated for the consecutive manipulated variables, for $n=1,\dots ,n_{\mathrm{y}}$, i.e., the algorithm goes to Step 1.

Having completed the above procedure, the values of the coefficients ${\psi }_{1,n},\dots ,{\psi }_{N,n}$ are calculated from Equation (6) for the found parameters $a_n$, $K_n$ and $m_n$, for all $n=1,\dots ,n_{\mathrm{y}}$. The flowchart of the tuning procedure is depicted in Figure 1.

Let us define the control error for the sampling instant k and process output n:

$$e_n\left(k\right)=y_n^{\mathrm{sp}}\left(k\right)-y_n\left(k\right)$$

(7)

Although different performance criteria may be used to assess control accuracy, in this work, we use the multicriteria approach based on the following three indicators [36,39]:

• the Sum of Squared Errors (${\mathrm{SSE}}_n$) of the control error $e_n$,
• the Huber standard deviation ${\sigma }_n$ of the control error $e_n$,
• the rational entropy $H_n$ of the control error $e_n$.

The above indicators are calculated for all $n=1,\dots ,n_{\mathrm{y}}$ controlled variables.

In addition to the Gauss function (6), the application of other functions such as bell-shaped, triangular and trapezoidal ones was studied in [34]. The tuning procedure is similar, i.e., we do not optimise all parameters at the same time, but the influence of the consecutive parameters is analysed, while the best ones are chosen. It turns out that the proposed Gauss function results in the best control quality. This conclusion is based on the authors’ experiences with various MPC control algorithms, applied to a few Single-Input Single-Output (SISO) and MIMO processes.

4. Simulation Results

The discussed tuning method is verified in the control system of a simulated binary distillation column, which separates a two component mixture of water and methanol. The linear approximation of the process (scaled) is described by the following transfer function model introduced by Wood and Berry [40]:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{Y}_1\left(s\right)\\ Y_2\left(s\right)\end{array}\right]& =\left[\begin{array}{cc}{\displaystyle \frac{12.8e^{-s}}{16.7s+1}}& {\displaystyle \frac{-18.9e^{-3s}}{21s+1}}\\ {\displaystyle \frac{6.6e^{-7s}}{10.9s+1}}& {\displaystyle \frac{-19.4e^{-2s}}{14.4s+1}}\end{array}\right]\left[\begin{array}{c}{U}_1\left(s\right)\\ U_2\left(s\right)\end{array}\right]+\left[\begin{array}{c}{\displaystyle \frac{3.8e^{-8s}}{14.9s+1}}\\ {\displaystyle \frac{4.9e^{-3s}}{13.2s+1}}\end{array}\right]F\left(s\right)\hfill \end{array}$$

(8)

The controlled variables, $Y_1$ and $Y_2$, are: the top product (distillate) composition and the bottom product composition, respectively. The manipulated variables, $U_1$ and $U_2$, are: the reflux flow rate and the vapour flow, respectively. The uncontrolled process input (the disturbance), F, is the flow rate of the input stream. All variables are deviations from a typical operating point, and all variables are dimensionless.

The sampling time of MPC is 1 min. To obtain a discrete process representation that is used to determine the step-response model necessary for prediction in the DMC algorithm [6], we use the forward Euler method and the zero-order holder. As a result, the following transfer function representation of the model (8) is obtained:

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{Y}_1\left(z\right)\\ Y_2\left(z\right)\end{array}\right]& =\left[\begin{array}{cc}{\displaystyle \frac{0.7440z^{-2}}{1-0.9419z^{-1}}}& {\displaystyle \frac{-0.8789z^{-4}}{1-0.9535z^{-1}}}\\ {\displaystyle \frac{0.5786z^{-8}}{1-0.9123z^{-1}}}& {\displaystyle \frac{-1.3015z^{-3}}{1-0.9329z^{-1}}}\end{array}\right]\left[\begin{array}{c}{U}_1\left(z\right)\\ U_2\left(z\right)\end{array}\right]+\left[\begin{array}{c}{\displaystyle \frac{0.2467z^{-9}}{1-0.9351z^{-1}}}\\ {\displaystyle \frac{0.3575z^{-4}}{1-0.9270z^{-1}}}\end{array}\right]F\left(z\right)\hfill \end{array}$$

(9)

A resulting step-response model is obtained and next used for prediction in the DMC algorithm [6].

The horizons used in the DMC algorithm are: the prediction horizon $N=30$, the control horizon $N_u=5$ and the horizon of process dynamics $D=100$ (the horizon of process dynamics is the number of step-response coefficients taken into account by the model used in the DMC algorithm) [39]. All values of the weights that penalise the excessive increments of the manipulated variables, i.e., ${\lambda }_{p,n}$, are equal to one. In this work, the objective of the controller is to effectively compensate the influence of three types of disturbances: changes of the flow rate of the input stream F, as well as the measurement noise of the first and the second process outputs, denoted as ${\delta }_1$ and ${\delta }_2$, respectively. As many as 600,000 min (samples) are considered, i.e., approximately 417 days. It must be emphasised that in this work, the disturbances are not generated artificially, but real industrial disturbances recorded in a distillation process are used [39]. All 600,000 samples of the disturbances are shown in Figure 2.

All coefficients ${\lambda }_{p,n}=1$ for $p=0,\dots ,N_{\mathrm{u}}-1$, $n=1,\dots ,n_{\mathrm{u}}$ [39]. The problem to solve is to carefully select the tuning parameters ${\psi }_{p,n}$ for $p=1,\dots ,N$, $n=1,\dots ,n_{\mathrm{y}}$. For this purpose, the tuning procedure presented in Section 3 is used.

Let us first concentrate on tuning the coefficients ${\psi }_{p,1}$, i.e., for the first process output, $p=1,\dots ,N$. Initially, the trajectory of the weights ${\psi }_{p,1}$ is parameterised by means of Equation (5), which is characterised by the parameters $m_1$ and $K_1$. In the first step of the tuning procedure, the influence of the coefficient $m_1$ on the performance criteria ${\mathrm{SSE}}_1$, ${\mathrm{SSE}}_2$, ${\sigma }_1$, ${\sigma }_2$, ${\mathrm{H}}_1$ and ${\mathrm{H}}_2$ is analysed, and the results are shown in Figure 3. We take into account all three criteria for the first process output, i.e., the indices ${\mathrm{SSE}}_1$, ${\sigma }_1$ and ${\mathrm{H}}_1$, and we choose the value $m_1=2$ because it gives the best control quality. Next, for the chosen parameter $m_1=2$, in the second step, we assess the influence of the parameter $K_1$, and it is shown in Figure 4. It is clear that the best results are obtained for $K_1=100$. Finally, the trajectory of the weights ${\psi }_{p,1}$ is parameterised by means of Equation (6), which is characterised by the parameter $a_1$ as well as the fixed parameters $m_1$ and $K_1$. Figure 5 depicts the influence of the parameter $a_1$ on the considered performance indices. The best set of control quality criteria is obtained for $a_1=2$.

Next, having tuned the coefficients ${\psi }_{p,1}$, let us concentrate on tuning the coefficients ${\psi }_{p,2}$, i.e., for the second process output, $p=1,\dots ,N$. In the first step of the tuning procedure, the influence of the coefficient $m_2$ on the performance criteria is analysed, and the results are shown in Figure 6. We take into account all three criteria for the first process output, i.e., the indices ${\mathrm{SSE}}_2$, ${\sigma }_2$ and ${\mathrm{H}}_2$. We easily find the best value $m_2=3$. Next, for the chosen parameter $m_2=3$, in the second step, we assess the influence of the parameter $K_2$, and it is shown in Figure 7. It is clear that the best results are obtained for $K_1=100$. Finally, Figure 8 depicts the influence of the parameter $a_1$ on the considered performance indices. The best results are obtained for $a_2=1$.

All things considered, the coefficients ${\psi }_{p,1}$ are characterised by the set of parameters $m_1=2$, $K_1=100$ and $a_1=2$, while the coefficients ${\psi }_{p,2}$ by the parameters $m_2=3$, $K_2=100$ and $a_1=1$. The actual values of the coefficients used in MPC optimisation are computed from Equation (6).

To evaluate the efficiency of the chosen parameters ${\psi }_{p,1}$ and ${\psi }_{p,2}$, let us consider the values of all six performance criteria. We not consider only the found coefficients, but additionally, we take into account the most practically cases in which all coefficients are constant, i.e., ${\psi }_{p,1}={\psi }_{p,2}$. In the latter case, we consider six values of the coefficients: 0.1, 1, 10, 100, 1000 and 10,000. The values of the obtained performance criteria are given in

Table 1. The values of the considered control quality indices for different constant parameters ${\psi }_{p,1}$, ${\psi }_{p,2}$, as well as for the tuned one (for all 600,000 samples).

Parameters	${\mathbf{SSE}}_1$	${\mathbf{SSE}}_2$	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	${\mathit{H}}_1$	${\mathit{H}}_2$
${\psi }_{p,1}={\psi }_{p,2}=0.1$	$3.91\times {10}^6$	$2.69\times {10}^6$	$1.36$	$1.0489$	$6.58\times {10}^1$	$5.80\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=1$	$2.21\times {10}^6$	$1.73\times {10}^6$	$1.02$	$8.23\times {10}^{-1}$	$4.94\times {10}^1$	$4.63\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=10$	$1.52\times {10}^6$	$1.39\times {10}^6$	$8.48\times {10}^{-1}$	$7.32\times {10}^{-1}$	$4.08\times {10}^1$	$4.14\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=100$	$1.31\times {10}^6$	$1.27\times {10}^6$	$7.87\times {10}^{-1}$	$6.92\times {10}^{-1}$	$3.85\times {10}^1$	$3.95\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=1000$	$1.28\times {10}^6$	$1.25\times {10}^6$	$7.79\times {10}^{-1}$	$6.85\times {10}^{-1}$	$3.83\times {10}^1$	$3.94\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=10,000$	$1.28\times {10}^6$	$1.25\times {10}^6$	$7.78\times {10}^{-1}$	$6.85\times {10}^{-1}$	$3.82\times {10}^1$	$3.92\times {10}^1$
${\psi }_{p,1}$, ${\psi }_{p,2}$ tuned	$9.62\times {10}^5$	$1.23\times {10}^6$	$6.74\times {10}^{-1}$	$6.99\times {10}^{-1}$	$3.32\times {10}^1$	$3.91\times {10}^1$

. For almost all considered indices (the only exception is ${\sigma }_2$, but the difference is not significant), the discussed tuning method gives the best results. For constant weights ${\psi }_{p,1}={\psi }_{p,2}$ greater than 10,000, the results do not change.

Let us now consider histograms of the control errors for constant values of the coefficients ${\psi }_{p,1}$, ${\psi }_{p,2}=0.1,1,10$ and for the tuned ones. They are shown in Figure 9. Two observations are possible. Firstly, our tuning method gives the most narrow histograms. Secondly, they do not have long “tails”, which means that our method eliminates large negative and positive control errors. Figure 10 depicts the first 1000 samples of the control errors for constant values of the coefficients ${\psi }_{p,1}$, ${\psi }_{p,2}=0.1,1,10$ and for the tuned ones. The obtained trends confirm the observations made on the basis of the histograms, i.e., our tuning method gives the best control quality.

In the second part of the simulations, we use the obtained set of tuning coefficients, but all gains of all input-output process channels are increased by 20%. Such a case is particularly interesting since the model used in MPC is (practically) never perfect.

Table 2. The values of the considered control quality indices for different constant parameters ${\psi }_{p,1}$, ${\psi }_{p,2}$, as well as for the tuned one (the imprecise model case, for all 600,000 samples).

Parameters	${\mathbf{SSE}}_1$	${\mathbf{SSE}}_2$	${\mathit{\sigma }}_1$	${\mathit{\sigma }}_2$	${\mathit{H}}_1$	${\mathit{H}}_2$
${\psi }_{p,1}={\psi }_{p,2}=0.1$	$2.73\times {10}^6$	$1.88\times {10}^6$	$1.14$	$8.77\times {10}^{-1}$	$5.55\times {10}^1$	$4.84\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=1$	$2.21\times {10}^6$	$1.73\times {10}^6$	$1.02$	$8.23\times {10}^{-1}$	$4.94\times {10}^1$	$4.63\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=10$	$1.06\times {10}^6$	$9.67\times {10}^5$	$7.07\times {10}^{-1}$	$6.10\times {10}^{-1}$	$3.43\times {10}^1$	$3.48\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=100$	$9.12\times {10}^5$	$8.85\times {10}^5$	$6.56\times {10}^{-1}$	$5.77\times {10}^{-1}$	$3.25\times {10}^1$	$3.35\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=1000$	$8.94\times {10}^5$	$8.73\times {10}^5$	$6.50\times {10}^{-1}$	$5.71\times {10}^{-1}$	$3.23\times {10}^1$	$3.34\times {10}^1$
${\psi }_{p,1}={\psi }_{p,2}=10,000$	$8.92\times {10}^5$	$8.72\times {10}^5$	$6.49\times {10}^{-1}$	$5.71\times {10}^{-1}$	$3.19\times {10}^1$	$3.27\times {10}^1$
${\psi }_{p,1}$, ${\psi }_{p,2}$ tuned	$6.68\times {10}^5$	$8.59\times {10}^5$	$5.62\times {10}^{-1}$	$5.83\times {10}^{-1}$	$2.75\times {10}^1$	$3.27\times {10}^1$

shows the values of all six performance criteria. Additionally, we consider six cases in which all coefficients are constant and have the values 0.1, 1, 10, 100, 1000 and 10,000, respectively. For almost all considered indices (the only exception is ${\sigma }_2$, but the difference is not significant), the discussed method gives the best results. For constant weights ${\psi }_{p,1}={\psi }_{p,2}$ greater than 10,000, the results do not change.

The histograms of the control errors are given in Figure 11, and Figure 12 depicts the first 1000 samples of the control errors. Similar to the perfect model case (Figure 9 and Figure 10), our tuning method gives the best results.

5. Conclusions

For MIMO processes, the classical optimisation approaches to finding the tuning parameters of MPC algorithms require very long vectors of decision variables. The resulting parameter optimisation task may be solved in simulations using some heuristic optimisation algorithms [12,32,33], but the huge computational complexity of such an approach in practice is an issue (the number of decision variables is large). This work presents a computationally simple approach to tuning of MPC parameters. The coefficients are not optimised directly, but they are parameterised. Moreover, the parameters are not calculated at the same time, but always, only one parameter is optimised. Multicriteria control quality assessment is possible. In this work, the following three indices are used: the sum of squared errors, the Huber standard deviation and the rational entropy. Unlike fully-fledged optimisation-based tuning methods, in our simulation-based approach, we simply assess how the consecutive parameters influence the value of the performance indices, and we choose the parameters for which the best results are possible. No numerical optimisation is used. The resulting tuning procedure is suboptimal, but our procedure is deliberately rather uncomplicated and is recommended to be used in industrial applications in which all tuning parameters are usually constant.

The efficiency of the discussed tuning method is demonstrated for a simulated MIMO distillation column. The classical process control task is considered, which is the compensation of the influence of the disturbances. It is shown that the discussed tuning method makes it possible to find the parameters that result in very good control quality, better than in the case of classical equal and large parameters. It is also shown that the determined set of parameters yields good control also in the imperfect model case.

Finally, let us stress the fact that our tuning method may be used to find parameters of numerous types of MPC algorithms. Classical MPC schemes based on linear models [6], as well as nonlinear ones [41] may be taken into account. Additionally, the constraints of different kinds may be included in the MPC optimisation problem.