Novel Strategy of Adaptive Predictive Control Based on a MIMO-ARX Model

Abstract

Many industrial processes include MIMO (multiple-input, multiple-output) systems that are difficult to control by standard commercial controllers. This paper describes a MIMO case of a class of SISO-APC (single-input, single-output adaptive predictive controller) based upon an ARX (autoregressive with exogenous variable) model. This class of SISO-APC based on ARX models has been successfully and extensively used in many industrial applications. This approach aims to minimize the barriers between the theory of predictive adaptive control and its application in the industrial environment. The proposed MIMO-APC (MIMO adaptive predictive controller) performance is validated with two simulated processes: a quadrotor drone and the quadruple tank process. In the first experiment the proposed MIMO APC shows ISE-IAE-ITAE performance indices improvements of up to 25%, 25.4% and 38.9%, respectively. For the quadruple tank process the water levels in the lower tanks follow closely the set points, with the exception of a 13% overshoot in tank 1 for the minimum phase behavior response. The controller responses show significant performance improvements when compared with previously published MIMO control strategies.

1. Introduction

The industrial process control sector has undergone a significant change in recent years with the incorporation of more complex, faster and multivariable processes. It is necessary to incorporate advanced control techniques that allow these processes to be controlled. Advanced process controls, now, are not only based on mathematical modeling and optimization objectives; it is important that new control systems should consider the interactions between hardware, software and mathematical modeling. One of several advanced control techniques is predictive control.

Model predictive control is becoming one of the most popular advanced control technique and has been use for more than four decades. This control strategy is based on numerical optimization [1]. This strategy considers the future values of a variable based on existing information on the process and the use of the explicit form of a mathematical model of the process’ internal dynamics, which is used to predict the evolution of the controlled variables over a prediction’s time horizon defined by the user. Thus it is possible to calculate the future manipulated variables to ensure that, in this horizon, the controlled variables converge to their reference values. If the predictive model is able to predict the behavior of the system, the variable under control will match the desired variable [2]. Predictive control can be used to control processes with relatively simple dynamics, as well as complex processes having long delay times, having non-minimum phase response or even having open-loop instabilities. For those applications, a well implemented predictive control technique can provide a high degree of robustness [3].

The unsuspected perturbations present in the industrial application field may be strong enough to alter the process nature; in this case, a mathematical model may not be able to accurately represent the process dynamic. By including an adaptation mechanism in the predictive controller, the APC can adjust the parameters of the predictive model and minimize the impact of the process dynamics related to the unknown disturbances in and time-varying nature of industrial processes. This approach has been implemented and documented in different application fields, such as automobile suspesion systems, oil field, aerospace and others [4,5,6,7,8,9,10,11].

The adaptive predictive control strategy offers an effective way to tackle the adaptation problems in multivariable control systems by including an estimated process model in the computation of control actions. The process model estimation is achieved through the incorporation of an identification mechanism on the adaptive control strategy. In this work, multivariable recursive least square (MRLS) is used as an adaptive mechanism tool. This paper considers linear difference equation models, an important class of stochastic models for describing dynamic behavior, which has received a great deal of attention and includes stationary models. One can find different structures of this kind of model in the literature, such as ARX/CAR models. On the other hand, a class of models for nonstationary dynamic behaviors are known as ARIMAX/CARIMA models. This work presents an extension to the (ARX) MIMO case of the previous (ARX) SISO predictive controller proposed in [5]. This SISO controller, based on the ARX model, is successfully and extensively used in many industrial applications [5,7,8,9,10,11]. Additionally, this SISO algorithm has been implemented and tested in embedded systems, showing improvements in terms of execution time compared with other MPC strategies found in the literature [12].

Although the APC algorithm has been successfully applied in the industry, to date its practical implementation has been mostly limited to the SISO case due to the fact that it can replace actual PID control loops. As such, industrial processes often present interactions between the input and output variables. The input–outputs interaction makes the design of the control strategy even more complex, due to the fact one several outputs can be affected by the change in one input; this interaction is also knows as a “coupled” system [13]. The presence of this coupling typically limits the performance of the control tools developed for SISO processes, because these controllers do not take this coupling into account. Therefore, the need for a performance improvement in the field of MIMO processes control has motivated the research of control techniques specifically developed for multivariable processes including linear and nonlinear systems, different control approaches like fuzzy controllers, neural network strategies, MPC and adaptive predictive controllers [14,15,16,17,18].

The structure of the paper is the following: Section 2 shows the general control scheme and the formulation, synthesis and recursive algorithm of the proposed MIMO-APC strategy. Section 3 briefly describes the adaptation mechanism. In Section 4, the proposed MIMO-APC is validated by implementing the proposed controller in two simulated processes. Finally, Section 5 presents the conclusions gathered and suggestions for further research.

2. Predictive Control for MIMO-ARX Process

2.1. General Control Strategy

Adaptive predictive control was first introduced as a new solution in the context of adaptive control [6]. The adaptive predictive control consist of a driver block, a predictive model and an adaptation mechanism. Figure 1 shows a global scheme of the APC strategy.

The Driver Block (DB) calculates the trajectory that guides the controlled variables to their references in an optimal way, based on the driver block equation [5,19]. The Predictive Model (PM) computes the control signals at every instant k; its objective is to lead the process output to meet the desired trajectory estimated by the DB. The predictive model will estimate the process output to calculate the prediction error every instant k using the estimated ARX model parameters at the instant k and the applied control signal $u(k-1)$. The ARX equation described in the predictive model is, at the beginning, an approximation model based on the process knowledge; this equation will be adjusted by the Adaptation Mechanism (AM), which will adjust the predictive model’s parameters based on prediction error until this tends asymptotically towards zero [20]; simultaneously, the current value of the values under control will feedback to the DB through the AM.

2.2. Representation of a MIMO-ARX Process

In MIMO systems, having more than one variable to be controlled has causal relationships with more than one variable to be manipulated, unlike SISO systems, wherein the input–output causal relationship is straightforward. The performance of the MIMO control system depends not only on the applied control strategy, but also upon a mathematical model that includes the interconnection between the variables of a coupled system. Based on the SISO-ARX model with wide application in real time control systems, this paper considers a multivariable, with autoregressive and exogenous variables, MIMO-ARX predictive model.

A MIMO process with a number of controlled variables r and a number of control signals s is considered. The set of input and output variables can be represented in a vector form as $\mathbf{y}\in {\mathbb{R}}^r$ and $\mathbf{u}\in {\mathbb{R}}^s$ respectively.

The output vector $\mathbf{y}$ is defined as follows:

$$\mathbf{y}\left(k\right)=\sum _{i=1}^n{\mathbf{A}}_i\mathbf{y}(k-i)+\sum _{i=1}^m{\mathbf{B}}_i\mathbf{u}(k-i)$$

(1)

$${\mathbf{A}}_i=\left[\begin{array}{cccc}{a}_i^{11}& a_i^{12}& \cdots & a_i^{1r}\\ a_i^{21}& a_i^{22}& \cdots & a_i^{2r}\\ ⋮& ⋮& \ddots & ⋮\\ a_i^{r1}& a_i^{r2}& \cdots & a_i^{rr}\end{array}\right]$$

(2)

$${\mathbf{B}}_i=\left[\begin{array}{cccc}{b}_i^{11}& b_i^{12}& \cdots & b_i^{1s}\\ b_i^{21}& b_i^{22}& \cdots & b_i^{2s}\\ ⋮& ⋮& \ddots & ⋮\\ b_i^{r1}& b_i^{r2}& \cdots & b_i^{rs}\end{array}\right]$$

(3)

where ${\mathbf{A}}_i\in {\mathbb{R}}^{r\times r}$ and ${\mathbf{B}}_i\in {\mathbb{R}}^{r\times s}$, n and m respectively represent the real order of the autorregresive and exogenous variables of the MIMO process.

2.3. Multivariable Predictive Control Strategy

In order to obtain a prediction model, a prediction horizon $N_p\ge 1$ is used. Thus, if the model (1) is extended to a prediction horizon $j=1,2,\cdots ,N_p$, a prediction model is obtained as follows:

$$\widehat{\mathbf{y}}(k+j|k)=\sum _{i=1}^{\widehat{n}}{\widehat{\mathbf{A}}}_i\widehat{\mathbf{y}}(k-i+j|k)+\sum _{i=1}^{\widehat{m}}{\widehat{\mathbf{B}}}_i\widehat{\mathbf{u}}(k-i+j|k)$$

(4)

where $\widehat{n}$ and $\widehat{m}$ represent the estimated order of the autorregresive and exogenous variables of the MIMO ARX Model, $\widehat{\mathbf{A}}$ and $\widehat{\mathbf{B}}$ are the estimated parameters of the MIMO model and the input–output vectors $\widehat{\mathbf{u}}-\widehat{\mathbf{y}}$ contain real and estimated values.

It is possible to rewrite (4) in order to split the terms of ${\widehat{\mathbf{A}}}^{\prime }s$ and ${\widehat{\mathbf{B}}}^{\prime }s$ associated with the terms $\mathbf{y}$ and $\mathbf{u}$ in the future and past, respectively.

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+j)=& \sum _{i=1}^{j-1}{\widehat{\mathbf{A}}}_i\widehat{\mathbf{y}}(k-i+j)+\sum _{i=j}^{\widehat{n}}{\widehat{\mathbf{A}}}_i\mathbf{y}(k-i+j)+\hfill \\ \hfill & \sum _{i=1}^{j-1}{\widehat{\mathbf{B}}}_i\widehat{\mathbf{u}}(k-i+j)+\sum _{i=j}^{\widehat{m}}{\widehat{\mathbf{B}}}_i\mathbf{u}(k-i+j)\hfill \end{array}$$

(5)

To conduct a prediction, Equation (5) at instant k will be extended j time steps as follows:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+1)=& {\widehat{\mathbf{A}}}_1\mathbf{y}\left(k\right)+{\widehat{\mathbf{A}}}_2\mathbf{y}(k-1)+\cdots +{\widehat{\mathbf{A}}}_{\widehat{n}}\mathbf{y}(k-\widehat{n}+1)+\hfill \\ \hfill & {\widehat{\mathbf{B}}}_1\mathbf{u}\left(k\right)+{\widehat{\mathbf{B}}}_2\mathbf{u}(k-1)+\cdots +{\widehat{\mathbf{B}}}_{\widehat{m}}\mathbf{u}(k-\widehat{m}+1)\hfill \end{array}$$

Now define $\widehat{\mathbf{y}}(k+1)$ as

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+1)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(1\right)}\mathbf{y}(k-i+1)+\sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(1\right)}\mathbf{u}(k-i+1)\hfill \end{array}$$

(6)

where ${\widehat{\mathbf{A}}}_i={\mathbf{E}}_i^{\left(1\right)}$ and ${\widehat{\mathbf{B}}}_i={\mathbf{G}}_i^{\left(1\right)}$.

Then, $\widehat{\mathbf{y}}(k+2)$ is obtained as follows:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+2)=& {\widehat{\mathbf{A}}}_1\widehat{\mathbf{y}}(k+1)+\sum _{i=2}^{\widehat{n}}{\widehat{\mathbf{A}}}_i\mathbf{y}(k-i+2)+\hfill \\ \hfill & {\widehat{\mathbf{B}}}_1\widehat{\mathbf{u}}(k+1)+\sum _{i=2}^{\widehat{m}}{\widehat{\mathbf{B}}}_i\mathbf{u}(k-i+2)\hfill \end{array}$$

(7)

Substituting Equations (6) in (7) obtains:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+2)=& {\widehat{\mathbf{A}}}_1\left[\sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(1\right)}\mathbf{y}(k-i+1)\right]+\hfill \\ & {\widehat{\mathbf{A}}}_1\left[\sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(1\right)}\mathbf{u}(k-i+1)\right]+\hfill \\ \hfill & \sum _{i=2}^{\widehat{n}}{\widehat{\mathbf{A}}}_i\mathbf{y}(k-i+2)+{\widehat{\mathbf{B}}}_1\widehat{\mathbf{u}}(k+1)+\hfill \\ & \sum _{i=2}^{\widehat{m}}{\widehat{\mathbf{B}}}_i\mathbf{u}(k-i+2)\hfill \end{array}$$

(8)

Associating terms in the function of $\mathbf{y}(k-i+1)$ and $\mathbf{u}(k-i+1)$, we get:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+2)=& \sum _{i=1}^{\widehat{n}}\left({\widehat{\mathbf{A}}}_1{\mathbf{E}}_i^{\left(1\right)}+{\widehat{\mathbf{A}}}_{i+1}\right)\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}\left({\widehat{\mathbf{A}}}_1{\mathbf{G}}_i^{\left(1\right)}+{\widehat{\mathbf{B}}}_{i+1}\right)\mathbf{u}(k-i+1)+\hfill \\ \hfill & {\widehat{\mathbf{B}}}_1\widehat{\mathbf{u}}(k+1)\hfill \end{array}$$

(9)

where ${\widehat{\mathbf{A}}}_{\mathbf{i}}=0$ for $i>\widehat{n}$ and ${\widehat{\mathbf{B}}}_i=0$ for $i>\widehat{m}$.

As in (6), Equation (9) can be expressed in terms of $\mathbf{E}$ and $\mathbf{G}$:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+2)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(2\right)}\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(2\right)}\mathbf{u}(k-i+1)+\hfill \\ \hfill & {\mathbf{G}}_1^{\left(1\right)}\widehat{\mathbf{u}}(k+1)\hfill \end{array}$$

(10)

Equation (10) shows that the future terms have disappeared; in this case, the output at instant $j=2$ depends only on past terms.

Now, from Equation (5), a horizon $j=3$ is used:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+3)=& \sum _{i=1}^2{\widehat{\mathbf{A}}}_i\widehat{\mathbf{y}}(k-i+3)+\sum _{i=3}^{\widehat{n}}{\widehat{\mathbf{A}}}_i\mathbf{y}(k-i+3)+\hfill \\ \hfill & \sum _{i=1}^2{\widehat{\mathbf{B}}}_i\widehat{\mathbf{u}}(k-i+3)+\hfill \\ & \sum _{i=3}^{\widehat{m}}{\widehat{\mathbf{B}}}_i\mathbf{u}(k-i+3)\hfill \end{array}$$

(11)

Substituting (6) and (10) in (11)

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+3)=& {\widehat{\mathbf{A}}}_1\left[\sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(2\right)}\mathbf{y}(k-i+1)\right]+\hfill \\ & {\widehat{\mathbf{A}}}_1\left[\sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(2\right)}\mathbf{u}(k-i+1)\right]+\hfill \\ \hfill & {\widehat{\mathbf{A}}}_1{\mathbf{G}}_1^{\left(1\right)}\widehat{\mathbf{u}}(k+1)+\hfill \\ \hfill & {\mathbf{A}}_2\left[\sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(1\right)}\mathbf{y}(k-i+1)\right]+\hfill \\ & {\mathbf{A}}_2\left[\sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(1\right)}\mathbf{u}(k-i+1)\right]+\hfill \\ \hfill & \sum _{i=3}^{\widehat{n}}{\widehat{\mathbf{A}}}_i\mathbf{y}(k-i+3)+\sum _{i=1}^2{\widehat{\mathbf{B}}}_i\widehat{\mathbf{u}}(k-i+3)+\hfill \\ \hfill & \sum _{i=3}^{\widehat{m}}{\widehat{\mathbf{B}}}_i\mathbf{u}(k-i+3)\hfill \end{array}$$

(12)

Associating terms of $\widehat{\mathbf{y}}(k-i+1)$ and $\widehat{\mathbf{u}}(k-i+1)$:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+3)=& \sum _{i=1}^{\widehat{n}}\left({\widehat{\mathbf{A}}}_1{\mathbf{E}}_i^{\left(2\right)}+{\widehat{\mathbf{A}}}_2{\mathbf{E}}_i^{\left(1\right)}+{\widehat{\mathbf{A}}}_{i+2}\right)\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}\left({\widehat{\mathbf{A}}}_1{\mathbf{G}}_i^{\left(2\right)}+{\widehat{\mathbf{A}}}_2{\mathbf{G}}_i^{\left(1\right)}+{\widehat{\mathbf{B}}}_{i+2}\right)\mathbf{u}(k-i+1)+\hfill \\ \hfill & \left({\widehat{\mathbf{A}}}_1{\mathbf{G}}_1^{\left(1\right)}+{\widehat{\mathbf{B}}}_2\right)\widehat{\mathbf{u}}(k+1)+{\widehat{\mathbf{B}}}_1\widehat{\mathbf{u}}(k+2)\hfill \end{array}$$

(13)

where, again ${\widehat{\mathbf{A}}}_i=0$ for $i>\widehat{n}$ and ${\widehat{\mathbf{B}}}_i=0$ for $i>\widehat{m}$.

Therefore, Equation (13) is rewritten as follows:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+3)=& \sum _{i=1}^{\widehat{n}}\left[\sum _{h=1}^2\left[{\widehat{\mathbf{A}}}_h{\mathbf{E}}_i^{(3-h)}\right]+{\widehat{\mathbf{A}}}_{i+2}\right]\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}\left[\sum _{h=1}^2\left[{\widehat{\mathbf{A}}}_h{\mathbf{G}}_i^{3-h}\right]+{\widehat{\mathbf{B}}}_{i+2}\right]\mathbf{u}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^2{\mathbf{G}}_1^{\left(i\right)}\widehat{\mathbf{u}}(k+3-i)\hfill \end{array}$$

(14)

As in (6) and (10) the prediction equation for $\widehat{\mathbf{y}}(k+3)$ is:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+3)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(3\right)}\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(3\right)}\mathbf{u}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^2{\mathbf{G}}_1^{\left(i\right)}\widehat{\mathbf{u}}(k+3-i)\hfill \end{array}$$

(15)

In general, the prediction equation for any horizon j can be expressed as follows:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+j)=& \sum _{i=1}^{\widehat{n}}\left[\sum _{h=1}^{j-1}\left[{\widehat{\mathbf{A}}}_h{\mathbf{E}}_i^{(j-h)}\right]+{\widehat{\mathbf{A}}}_{j+i-1}\right]\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}\left[\sum _{h=1}^{j-1}\left[{\widehat{\mathbf{A}}}_h{\mathbf{G}}_i^{j-h}\right]+{\widehat{\mathbf{B}}}_{j+i-1}\right]\mathbf{u}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{j-1}{\mathbf{G}}_1^{\left(i\right)}\widehat{\mathbf{u}}(k+j-i)\hfill \end{array}$$

(16)

as in previous developments, Equation (16) can be rewritten as

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+j)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(j\right)}\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{\widehat{m}}{\mathbf{G}}_i^{\left(j\right)}\mathbf{u}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^{j-1}{\mathbf{G}}_1^{\left(i\right)}\widehat{\mathbf{u}}(k+j-i)\hfill \end{array}$$

(17)

In (17) the second summation contains $\widehat{\mathbf{u}}\left(k\right)$, for practical purposes, this term is separated:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+j)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(j\right)}\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=2}^{\widehat{m}}{\mathbf{G}}_i^{\left(j\right)}\mathbf{u}(k-i+1)+\hfill \\ \hfill & {\mathbf{G}}_1^{\left(j\right)}\widehat{\mathbf{u}}\left(k\right)+\hfill \\ \hfill & \sum _{i=1}^{j-1}{\mathbf{G}}_1^{\left(i\right)}\widehat{\mathbf{u}}(k+j-i)\hfill \end{array}$$

(18)

Reordering the terms of ${\mathbf{G}}_1^{\left(i\right)}$ associated with the control actions in the present and future, Equation (18) is expressed as:

$$\begin{array}{cc}\hfill \widehat{\mathbf{y}}(k+j)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(j\right)}\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=2}^{\widehat{m}}{\mathbf{G}}_i^{\left(j\right)}\mathbf{u}(k-i+1)+\hfill \\ \hfill & \sum _{i=1}^j{\mathbf{G}}_1^{\left(i\right)}\widehat{\mathbf{u}}(k+j-i)\hfill \end{array}$$

(19)

Equation (19) shows the predicted output at instant $k+j$ based on the process parameters, inputs and outputs at every instant k. The recursiveness of the parameters matrices are defined as follow:

$$\begin{array}{cc}\hfill {\mathbf{E}}_i^{\left(j\right)}& =\sum _{h=1}^{j-1}\left[{\widehat{\mathbf{A}}}_p{\mathbf{E}}_i^{(j-h)}\right]+{\widehat{\mathbf{A}}}_{i+j-1}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\mathbf{G}}_i^{\left(j\right)}& =\sum _{h=1}^{j-1}\left[{\widehat{\mathbf{A}}}_p{\mathbf{G}}_i^{(j-h)}\right]+{\widehat{\mathbf{B}}}_{i+j-1}\hfill \end{array}$$

(21)

under the following initial conditions:

$$\begin{array}{cc}\hfill {\mathbf{E}}_i^1=& {\widehat{\mathbf{A}}}_{i}i=\{1,2,\dots ,\widehat{n}\}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\mathbf{E}}_i^{\left(j\right)}=& {\mathbb{O}}_{r\times r}\widehat{m}>\widehat{n}\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\mathbf{G}}_i^1=& {\widehat{\mathbf{B}}}_{i}i=\{1,2,\dots ,\widehat{m}\}\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathbf{G}}_i^{\left(j\right)}=& {\mathbb{O}}_{r\times s}\widehat{n}>\widehat{m}\hfill \end{array}$$

(25)

2.4. Control Law

Based on the MIMO predictive control law (26):

$$\widehat{\mathbf{y}}(k+N_p|k)={\widehat{\mathbf{y}}}_d(k+N_p|k)$$

(26)

Assuming that the value of the control signals of the system remains constant over the entire prediction horizon ($\widehat{\mathbf{u}}\left(k\right|k)=\widehat{\mathbf{u}}(k+1|k)=\cdots =\widehat{\mathbf{u}}(k+N_p-1|k)$) and considering the predictive control law (26), Equation (19) leads to the control law expression as follows:

$$\begin{array}{cc}\hfill {\widehat{\mathbf{y}}}_d(k+N_p)=& \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(N_p\right)}\mathbf{y}(k-i+1)+\hfill \\ \hfill & \sum _{i=2}^{\widehat{m}}{\mathbf{G}}_i^{\left(N_p\right)}\mathbf{u}(k-i+1)+\hfill \\ \hfill & {\mathbf{H}}^{\left(N_p\right)}\widehat{\mathbf{u}}\left(k\right|k)\hfill \end{array}$$

(27)

where ${\mathbf{H}}^{\left(N_p\right)}={\sum }_{i=1}^{N_p}{\mathbf{G}}_1^{\left(i\right)}$.

From (27) the vector containing the control signals $\widehat{\mathbf{u}}\left(k\right|k)$ for the entire prediction horizon $N_p$ at time k is:

$$\begin{array}{cc}\hfill \widehat{\mathbf{u}}\left(k\right|k)={\left({\mathbf{H}}^{\left(N_p\right)}\right)}^{-1}& [{\mathbf{y}}_d(k+N_p|k)-\hfill \\ & \sum _{i=1}^{\widehat{n}}{\mathbf{E}}_i^{\left(N_p\right)}\mathbf{y}(k-i+1)-\hfill \\ & \sum _{i=2}^{\widehat{m}}{\mathbf{G}}_i^{\left(N_p\right)}\mathbf{u}(k-i+1)]\hfill \end{array}$$

(28)

The simplicity of the calculation of this control law and its easy implementation on different hardware/software platforms makes this approach attractive for real industrial applications.

2.5. Desired Trajectory for the Multivariable Predictive Control

The control law (28) is used to implement the new proposed multivariable predictive control. Since the desired output trajectories are independent, it is not necessary to apply a multivariable recursion algorithm to calculate them. The desired trajectory for each controlled variable is generated using the recursive SISO algorithm, described in [20] as follows:

$$\begin{array}{cc}\hfill y_d^j(k+N_p|k)=& \sum _{i=1}^p{\phi }_i^{\left(N_p\right)}y(k+1-i)+\hfill \\ & \sum _{i=2}^q{\delta }_i^{\left(N_p\right)}{y}_{sp}(k+1-i)+{\mu }^{\left(N_p\right)}{y}_{sp}\left(k\right)\hfill \end{array}$$

(29)

where

$$\begin{array}{cc}\hfill {\phi }_i^{\left(j\right)}& ={\phi }_1^{(j-1)}{\alpha }_i+{\phi }_{i+1}^{j-1}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\delta }_i^{\left(j\right)}& ={\phi }_1^{(j-1)}{\beta }_i+{\delta }_{i+1}^{j-1}\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\mu }^{\left(j\right)}& ={\delta }_1^{\left(N_p\right)}+{\delta }_1^{(N_p-1)}+\cdots +{\delta }_1^{\left(1\right)}\hfill \end{array}$$

(32)

$\forall i\in \{1,2,\dots ,p\},\forall j\in \{2,3,\dots ,N_p\}$.

The coefficients $\alpha$ and $\beta$ are the coefficients of the driver block, the selection of these parameters is thoroughly discussed in [5]. Note that if the same dynamic behavior for each output is considered, the same coefficients are used. If different dynamics are necessary, then the desired coefficients must be defined for each output independently.

Having generated the desired trajectory for each of the controlled variables, the vector containing those variables for each instant k is defined as:

$$\begin{array}{c}\hfill {\mathbf{y}}_d(k+N_p|k)=\left[\begin{array}{c}{y}_d^1(k+N_p|k)\\ y_d^2(k+N_p|k)\\ ⋮\\ y_d^r(k+N_p|k)\end{array}\right]\end{array}$$

(33)

Up to this point we have presented the extension of the recursive algorithm for a multivariable APC strategy. As mentioned earlier in this section, the main advantage of the proposed algorithm is that it is based on a multivariable system with a coupling between all input and output variables of the system. Systems based upon output decoupling algorithms are shown in [3], from which we can assert that the proposed case is special.

3. Adaptation Mechanism Based on the Least Squares Algorithm

The method of least squares identification is a strategy for identifying dynamical systems by which an ideal model with an ARX structure is adjusted from an initial estimation; a recursive algorithm updates the model parameters when new data from the inputs and outputs of the system are available, thus making the identification of a process over time possible. This recursive property allows the least squares strategy to be used as the adaptation mechanism in the proposed adaptive control strategy.

The multivariable recursive least Squares (MRLS) equations for identifying the parameters of a multivariable system with s number of inputs and r number of outputs can be summarized as follows:

$$\begin{array}{cc}\hfill {\mathbf{e}}_k& ={\mathbf{y}}_k-{\Theta }_{k-1}^T{\Psi }_k\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill C& =\left[\frac{1}{a}+{\Psi }_k^T\frac{P_{k-1}}{\gamma }{\Psi }_k\right]\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill L_k& =\frac{1}{\gamma C}{P}_{k-1}{\Psi }_k\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\widehat{\Theta }}_k& =\left[\begin{array}{c}{\widehat{\mathbf{A}}}_1^T\\ {\widehat{\mathbf{A}}}_2^T\\ ⋮\\ {\widehat{\mathbf{A}}}_{\widehat{n}}^T\\ {\widehat{\mathbf{B}}}_1^T\\ {\widehat{\mathbf{B}}}_2^T\\ ⋮\\ {\widehat{\mathbf{B}}}_{\widehat{m}}^T\end{array}\right]={\widehat{\Theta }}_{k-1}+L_k{\mathbf{e}}_k^T\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill P_k& =\frac{1}{\gamma }\left[I-{\Psi }_k^T{L}_k\right]P_{k-1}\hfill \end{array}$$

(38)

where at instant k;

${\mathbf{e}}_k\in {\mathbb{R}}^r$ is the prediction error vector, ${\mathbf{y}}_k\in {\mathbb{R}}^r$ is the output vector, ${\Psi }_k\in {\mathbb{R}}^{r·\widehat{n}+s·\widehat{m}}$ is themultivariable regressive input-output vector, ${\widehat{\Theta }}_k\in {\mathbb{R}}^{(r·\widehat{n}+s·\widehat{m})\times \left(max(r,s)\right)}$ defines the parameters matrix, $L_k\in {\mathbb{R}}^{r·\widehat{n}+s·\widehat{m}}$ represents the gain vector, $P_k\in {\mathbb{R}}^{(r·\widehat{n}+s·\widehat{m})\times (r·\widehat{n}+s·\widehat{m})}$ is the inverse covariance matrix and $\gamma \in \mathbb{R}$ is the forgetting factor.

As can be seen from the proposed algorithm for the adaptation mechanism, it is necessary to initialize the covariance matrix, $P_k$, and the estimated parameters vector ${\widehat{\Theta }}_k$. The covariance matrix, $P_k$, could be initialized with high values on the diagonal and the vector, ${\widehat{\Theta }}_k$, with zeros.

In order to optimize the adaptation mechanism, a forgotten factor $0.98\ge \gamma <0.995$ is considered. Additionally, a reset strategy for the covariance matrix is applied when the estimation error between the estimated outputs and measured outputs exceeds a defined limit value. This restart is translated as a diagonal matrix with high values in the covariance matrix and has the effect of increasing the speed of adaptation.

The computation of the driver block, the prediction model and the adaptive mechanism can be summarized in the algorithm shown in Algorithm 1.

Algorithm 1: MIMO-APC algorithm

1.    Initialization 1.1. Define the sampling time T 1.2. Set the initial adaptive model (${\widehat{\mathbf{A}}}^{\prime }s,{\widehat{\mathbf{B}}}^{\prime }s$) 1.3. Define the PDT parameters ($\alpha ,\beta$) 1.4. Define the prediction horizon $N_p$ 2.    Driver Block 2.1. Update ${\mathbf{y}}_{sp}\left(k\right)$ and measure $\mathbf{y}\left(k\right)$ 2.2. Compute ϕ, δ and μ 2.3. Compute the DDT at $N_p$ prediction horizon 3.    Control Law 3.1. Compute ${\mathbf{E}}^{\prime }s$, ${\mathbf{G}}^{\prime }s$ and $\mathbf{H}$ coefficients 3.2. Compute of $\mathbf{u}\left(k\right)$ 3.3. Apply the control vector $\mathbf{u}\left(k\right)$ 4.    Adaptation Mechanism 4.1. Compute ${\Psi }_k$ 4.2. Compute the estimation error ${\mathbf{e}}_k$ 4.3. Compute ${\widehat{\Theta }}_k$ 4.4. Compute $P_k$ 4.5. Update the predictive model parameters ($\widehat{\mathbf{A}},\widehat{\mathbf{B}}$) 4.6. Repeat from 1.3

4. Results and Discussion

In the following section the results of two cases of study wherein the APC was implemented are shown. The simulation and some performance indexes under different conditions are documented, as well.

4.1. Quadrotor Flight Control

Interest in unmanned aerial vehicles (UAVs), also call drones, has grown over the last decade. The reason can be associated to many factors, but most important among them is reaching places out of the visual line of sight with a certain level of autonomy. It is increasingly common to see UAVs applications in the fields of agriculture, military, 2D and 3D mapping, logistics, medical, surveillance and many others.

A quadrotor drone is an small aircraft that has four rotors with rotating blades that enables the UAV to take off and tilt forward for propulsion while in flight. The flight movement and speed can be changed by varying the speeds of each independent blade, giving the drone six degrees of freedom (DoF).

The Drone parameters and details used in this experiment are well-described in [21]. Nevertheless, a brief description is presented in this section in order to clarify the predictive model used in the proposed MIMO-APC strategy.

4.1.1. Quadrotor Dynamics

For this simulation we considered the same simulation parameters as reported in [22]. The four degrees of freedom of the simulated drone provide attitude and position. Movements are thus achieved:

• Pitch: By rotational movement along the transverse axis y, translational movement on x axis is made.
• Roll: By rotational movement along longitudinal axis x, translational movement on y axis is made.
• Yaw: Rotational movement along the z axis.
• Throttle: Translational movement on the z axis.

The control parameters’ range of $[-1, 1]$ represents the min and max movement percentage value and this is given to the internal control. Additionally, $\varphi$ represents the roll angle in radians, $\theta$ the pitch angle in radians, $\dot{z}$ the vertical speeds in $m/s$ and $\dot{\psi }$ the angular speeds in $m/s$.

For this simulation the quadrotor model is considered a linear time-invariant system, as described in [22] a sampling time of 66 $ms$ and transfer functions equations are defined by Equations (39)–(42).

$$H_x\left(q^{-1}\right)=\frac{x\left(q^{-1}\right)}{u_x\left(q^{-1}\right)}=\frac{0.01477q^{-1}+0.01446q^{-2}}{1-1.939q^{-1}+0.9391q^{-2}}$$

(39)

$$H_y\left(q^{-1}\right)=\frac{y\left(q^{-1}\right)}{u_y\left(q^{-1}\right)}=\frac{0.01477q^{-1}+0.01446q^{-2}}{1-1.939q^{-1}+0.9391q^{-2}}$$

(40)

$$H_z\left(q^{-1}\right)=\frac{z\left(q^{-1}\right)}{\dot{z}\left(q^{-1}\right)}=\frac{0.00621q^{-1}+0.005644q^{-2}}{1-1.751q^{-1}+0.7505q^{-2}}$$

(41)

$$H_{yaw}\left(q^{-1}\right)=\frac{\psi \left(q^{-1}\right)}{\dot{\psi }\left(q^{-1}\right)}=\frac{0.1137q^{-1}+0.05722q^{-2}}{1-1.119q^{-1}+0.119q^{-2}}$$

(42)

Unlike the strategy reported in [21], here, the MIMO-APC is used to control the whole quadrotor behavior, instead of using a SISO controller for each DoF of the drone. Thus, MIMO-ARX model of the quadrotor is given by Equation (43).

$$\widehat{Y}\left(k\right)={\widehat{A}}_1\widehat{Y}(k-1)+{\widehat{A}}_2\widehat{Y}(k-2)+{\widehat{B}}_{1}U(k-1)+{\widehat{B}}_{2}U(k-2)$$

(43)

where

$$\begin{array}{cc}\hfill {\widehat{A}}_1=& \left[\begin{array}{cccc}1.9390& 0& 0& 0\\ 0& 1.9390& 0& 0\\ 0& 0& 1.7510& 0\\ 0& 0& 0& 1.1190\end{array}\right],\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill {\widehat{A}}_2=& \left[\begin{array}{cccc}-0.9391& 0& 0& 0\\ 0& -0.9391& 0& 0\\ 0& 0& -0.7505& 0\\ 0& 0& 0& -0.1190\end{array}\right],\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\widehat{B}}_1=& \left[\begin{array}{cccc}0.0148& 0& 0& 0\\ 0& 0.0148& 0& 0\\ 0& 0& 0.0062& 0\\ 0& 0& 0& 0.1137\end{array}\right],\hfill \end{array}$$

(46)

$$\begin{array}{cc}\hfill {\widehat{B}}_2=& \left[\begin{array}{cccc}0.0145& 0& 0& 0\\ 0& 0.0145& 0& 0\\ 0& 0& 0.0056& 0\\ 0& 0& 0& 0.0572\end{array}\right]\hfill \end{array}$$

(47)

The model described by Equation (43) has been selected as the drone time-invariant model, for this reason, the analysis of the adaptation mechanism performance is not considered in this experiment case.

4.1.2. Simulation Results

The first test consists of a tracking experiment of four points in a 3D space. In this way, the drone will simulate free flight-path control. This test is used to evaluate the performance of the proposed MIMO-APC, described in Section 3. The result obtained from this test for a prediction horizon $N_p=4$ is shown in Figure 2. This result indicates that a better reference tracking can be achieved using the MIMO-APC controller proposed in this paper in comparison with the multiple SISO controller implemented in [22].

Now, in order to have a better comparison of both controllers, the performance of all degrees of freedom are analyzed separately. Figure 2 shows the performance of the proposed MIMO-APC controller according to translational movement changes in the x, y and z axes. Notable improvements were achieved using the proposed MIMO-APC. Figure 3 and Figure 4 show a reduction of the settling time for the translation movement in axes x and y respectively, and the saturation time has decreased in comparison with the results presented in [21,22]. Figure 5 shows the performance of the elevation controller (z) where a reduction of the settling from 4 to 2 s and the saturation time for the control action can be observed. In both cases, a desired behavior without overshoot is presented.

The performance of the APC is compared with the controllers reported in [21,22] using the same performance indexes as ISE, IAE and ITAE.

Table 1. Performance index for path-following strategies.

CONTROLLER	ISE	Δ	IAE	Δ	ITAE	Δ
PD X	118.12	0	159.99	0	1862.70	0
EPSAC X	100.81	*	134.20	16.1	1279.60	31.3
MIMO-APC X	25	88.53	25.4	119.33	1138.00	38.9
PD Y	87.28	0	164.65	0	2360.50	0
EPSAC Y	68.89	21	131.68	20	1990.90	15.6
MIMO-APC Y	59.60	31.7	101.32	38.4	1438.70	39

Δ are improvements shown in percentages with respect to the PD..X or PD..Y case. (*) here the EPSAC controller did not improve the performance index criteria.

illustrates the X and Y performance comparisons for the path-following task. The drone´s performance shows improvements of 25%, 25.4% and 38.9% in X for the ISE, IAE and ITAE indices respectively, while improvements of 31.7%, 38.4% and 39% are shown in Y for the ISE, IAE and ITAE indices, respectively, in comparison with the reported PD controller. The proposed MIMO-APC, applied to the quadrotor drone, shows a remarkable reduction in the error of the path-following flight behavior, in comparison with the result reported with the decentralized MPC controllers presented in [21].

4.2. Quadruple-Tank Process Control

In this example the adaptive predictive controller is used to control the quadruple tank process. This process is considered a multivariable and coupled system that consists of four interconnected tanks that can configured to show the effect of multivariable zero (minimum and non minimum phase) and nonlinear behavior [17,23,24]. Figure 6 shows a schematic diagram of the plant.

The regulation problem in this system focuses on the lower tanks, A2 and B2, that are directly filled by two pumps, BP and AP, and indirectly by Tank A1 and tank B1, respectively. Pump BP supplies with fluid tanks A2 and B1 using a three-way valve, while pump AP supplies B2 and A1. On each tank (A1, A2, B1 and B2) a flow meter, regulated by a pneumatic valve and level sensors, have been included. The three-way valves are manually emulated; the positions of these two valves determine the non-minimum and minimum phase behaviors for the linear system, due to the location of the zero in the transfer function matrix.

The control inputs $u_1$, $u_2$ are the (0–10 V) simulated signal voltages applied to the two pumps. The simulated control outputs $y_1$, $y_2$ are (0–10 V) signal voltages representing the levels in the lower tanks.

4.2.1. Linear Model

Consider the following state space model of the quadruple tank process described in [17]:

$$\begin{array}{cc}\hfill \dot{x}& =Ax+Bu\hfill \\ \hfill y& =Cx\hfill \end{array}$$

(48)

where

$$\begin{array}{cc}\hfill A& =\left[\begin{array}{cccc}-\frac{1}{T_1}& 0& \frac{A_3}{A_1{T}_3}& 0\\ 0& -\frac{1}{T_2}& 0& \frac{A_4}{A_2{T}_4}\\ 0& 0& -\frac{1}{T_3}& 0\\ 0& 0& 0& -\frac{1}{T_4}\end{array}\right],\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill B& =\left[\begin{array}{cc}\frac{{\gamma }_1{k}_1}{A_1}& 0\\ 0& \frac{{\gamma }_2{k}_2}{A_2}\\ 0& \frac{(1-{\gamma }_2)k_2}{A_3}\\ \frac{(1-{\gamma }_1)k_1}{A_4}& 0\end{array}\right],\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill C& =\left[\begin{array}{cccc}{k}_c& 0& 0& 0\\ 0& k_c& 0& 0\end{array}\right],\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill D& =0.\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill u& =\left[\begin{array}{c}▵u_1\\ ▵u_2\end{array}\right],\hfill \end{array}$$

(53)

$$\begin{array}{cc}\hfill x& =\left[\begin{array}{c}▵h_1\\ ▵h_2\\ ▵h_3\\ ▵h_4\end{array}\right],\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill y& =\left[\begin{array}{c}▵y_1\\ ▵y_2\end{array}\right]\hfill \end{array}$$

(55)

and

$$\begin{array}{c}\hfill T_i=\frac{A_i}{a_i}\sqrt{\frac{2h_i^0}{g}}\end{array}$$

(56)

The transfer function matrix from u to y is given by

$$\begin{array}{c}\hfill G\left(s\right)=\left[\begin{array}{cc}\frac{{\gamma }_1{c}_1}{1+s{T}_1}& \frac{(1-{\gamma }_2)c_1}{(1+s{T}_1)(1+s{T}_3)}\\ \frac{(1-{\gamma }_1)c_2}{(1+s{T}_2)(1+s{T}_4)}& \frac{{\gamma }_2{c}_2}{1+s{T}_2}\end{array}\right]\end{array}$$

(57)

where $c_1=\frac{T_1{k}_1{k}_c}{A_1}$ and $c_2=\frac{T_2{k}_2{k}_c}{A_2}$

In the system, the manual position adjustment of the valve determines the minimum phase or non-minimum phase behavior of the system. These valve adjustments can be written with respect to flow ratios ${\gamma }_1$ and ${\gamma }_2$, as shown in

Table 2. Location of zeros on the linearized system as a function of the flow ratios ${\gamma }_1$ and ${\gamma }_2$.

	${\mathit{z}}_1$	${\mathit{z}}_2$	System Behavior
$1<{\gamma }_1+{\gamma }_2\le 2$	negative	negative	minimum phase
${\gamma }_1+{\gamma }_2=1$	zero	negative	boundary
$0<{\gamma }_1+{\gamma }_2\le 1$	positive	negative	non-minimum phase

.

The quadruple tank has approximately the physical constants shown in

Table 3. Parameter values for the quadruple tank simulation.

Parameter	Units	Value
$A_i$	${\mathrm{cm}}^2$	4.9
$a_i$	${\mathrm{cm}}^2$	0.03
$k_i$	V/cm	1.6
$k_c$	V/cm	0.5
g	${\mathrm{cm}/\mathrm{s}}^2$	981

.

The following parameters have been used for all experiments: $\widehat{n}=\widehat{m}=2$, $s=r=2$. The driver block for both outputs is represented by q first-order equation with a $\tau =3$.

4.2.2. Experiment Result

The first part of the experiment consists of adjusting the flow valves so that the system behavior is its minimum phase. In this case, the value of ${\gamma }_1$ has been changed from 0.7 to 0.6, in online mode (30 s). Based on this change, it is necessary to update the model parameters with the adaptive mechanism and to adjust the MIMO-APC settings. Figure 7 and Figure 8 show the performance of the MIMO-APC for the system operating in both minimum and non-minimum phases. Figure 7 shows that the level controllers closely follow the set points with the exception of a 13% overshoot in the tank A2 response, while Figure 8 illustrates very close set-point tracking throughout the experiment.

The adaptation mechanism starts adjusting parameters after a change in the setpoint. As in [23], in the non-minimum phase behavior, the value of ${\gamma }_1+{\gamma }_2<1$. In this experiment, the cases of ${\gamma }_1=0.5$ and ${\gamma }_1=0.4$ are considered. In Figure 8 the behavior of the controller is shown when it is operating in the non-minimum phase. In both experiment cases, a forgetting factor = 1 in the adaptation mechanism and a prediction horizon $N_p=4$ in the MIMO-APC strategy are considered.

Once the predictive model is adjusted, the effect of the adaptation mechanism is based upon RLS and the performance of the MIMO-APC, in both the minimum or non-minimum phase experiment cases, can be observed.

5. Conclusions and Future Work

This paper showed a formulation of a novel MIMO-APC approach for MIMO-ARX models. The development focused on an APC strategy applied to multivariable processes based on a MIMO-ARX model as a predictive model without restrictions in the number of inputs/outputs, of the system order or prediction horizon. The application of the obtained MIMO-APC algorithm in two benchmark application examples was presented, showing important improvements of as much as 39% in the quadrotor drone experiment case and excellent set-point tracking behavior in the quadruple-tank experiment case, both in comparison with the control strategies previously reported in the literature.

There are two aspects from which our future work will be considered. One is that, during the calculation of the desired control path, in the driver block $Y_d\left(k\right)$ will be formulated considering the actuator’s constraints; this helps to reduce or eliminate the saturation of the control action of the APC. Second, our objective is to achieve the application of an embedded APC controller and/or implement this APC in programmable logic controllers.