Tuning and Feasibility Analysis of Classical First-Order MIMO Non-Linear Sliding Mode Control Design for Industrial Applications

Abstract

Model-based control techniques have been gaining more and more interest these days. These complex control systems are mostly based on theories, such as feedback linearization, model predictive control, adaptive and robust control. In this paper the latter approach is investigated, in particular, sliding mode (SM) control is analyzed. While several works on the description and application of SM control on single-input single-output systems can easily be found, its application on multi-input multi-output systems is not examined in depth at the same level. Hence, this work aims at formalizing some theoretical complements about the necessary conditions for the feasibility of the SM control for multi-input-multi-output systems. Furthermore, in order to obtain the desired performance from the control system, a method for parameter tuning is proposed in the particular case in which the relative degree of the controlled channels is equal to one. Finally, a simple control problem example is shown with the aim of stressing the benefits derived from the application of the theoretical complements described here.

1. Introduction

Thanks to the increasing spread and development of performing microprocessors, the computation time of controller processors has been drastically reduced over the past decade, allowing the implementation of more and more complex control architectures. This fact, combined with the current more and more pressing issue of increased efficiency in each industrial field (e.g., in power generation [1,2,3,4], the transport sector [5,6,7], energy utilization [8,9], renewable energy [10,11,12], and so on), strongly favors the development and implementation of more sophisticated control theories, such as model-based ones. Among model-based control techniques experiencing wide popularity, model predictive control (MPC) [13,14,15], feedback linearization (FBL) [16,17], and sliding mode (SM) control [18] stand out. In particular, the SM control theory has been receiving growing interest since the early 1970s [19] and has found application in numerous industrial applications; to name a few, it was employed in References [20,21] for photovoltaic systems, in References [22,23] for gas turbines, in Reference [24], and Reference [25] for hybrid electric vehicles, and the list goes on. Many studies have been developed on SM, each one addressing a particular limitation, such as chattering [26], asymptotical convergence of the state variables to the desired value [27], insensitivity with respect to matched uncertainties [28], higher order derivative requirements for sliding surface design [29], and so on. Many other approaches can be found; for instance, a method to design integral sliding manifolds in the presence of additive unmatched uncertainties was discussed in Reference [30], while an adaptive second order sliding mode was proposed in Reference [31] in order to avoid the necessity of any a priori knowledge of the uncertainty upper bounds. In Reference [32] a peculiar approach was presented, where the physical model of the controlled system is not needed for the controller design, and the chattering behavior (determined by SM control) affecting the inverter switching frequency was analyzed and successfully addressed in Reference [33]. This paper has a slightly different goal than those previously mentioned. As a matter of fact, the attention here is focused on SM theory application to the non-linear dynamic systems presented in Reference [34] and to its further development and clarification for multi-input-multi-output systems provided in References [35,36,37]. In particular, a novel tuning method and the theoretical feasibility condition of a multi-input-multi-output (MIMO) classical first-order SM controller are discussed. With specific reference to theoretical assessments of SM application, Reference [34] provides all the necessary information for non-linear single-input-single-output (SISO) systems, but the extension of the theory to non-linear MIMO systems has not been deeply analyzed. Moreover, no hint on how to obtain the majorant matrix of the errors introduced by the multiplicative uncertainties was given in References [34,37].

In addition, for second order (or higher) SISO systems (and for MIMO system channels of which the order is at least two), tuning the SM controller parameters is a simple and straightforward action. Nevertheless, if the sliding variable corresponds to the tracking error (i.e., the dynamic order of the controlled channel is equal to one, as in Reference [23]) the classical SM tuning parameters disappear. Hence, in order to guarantee the desired performance of the control system, an alternative tuning method of the controller has to be investigated.

Therefore, motivated by what is previously mentioned, it appears from References [34,37] that a deeper analysis on how to obtain the majorant matrix of the errors, introduced by the multiplicative uncertainties, is necessary to provide a useful tool able to make the SM controller design more straightforward. The research developed in References [34,35,36,37,38] discuss the tuning method of the sliding function coefficients; however, they do not take into account the possibility of exploiting the coefficient compensating for the uncertain behavior of the system for further controller tuning, hence, allowing a higher performance to be obtained. As a consequence, all the first-order SM controllers could benefit from this new insight. For example, the performance determined by the controller in Reference [22] could be further enhanced.

In view of this state-of-the-art, this paper has three goals: (i) examining in depth the SM theory for MIMO systems presented in Reference [34] extending the definition of the necessary conditions to be verified in the control system matrix in order to successfully design an SM controller; (ii) providing a theoretical approach for the definition of the majorant matrix of multiplicative uncertainties; and (iii) giving an outline of how to handle controller tuning when the relative degree of the controlled channel is equal to one.

The main contributions of this paper can be summarized as follows:

• The classical first-order SM control for MIMO system feasibility conditions are deeply analyzed and correlated to the entity of the uncertainties affecting the real system.

• A straightforward procedure to obtain the majorant matrix of the errors introduced by the multiplicative uncertainties is given.

• A novel method to properly tune SM controllers exploiting the coefficients which guarantee the sliding condition to be verified is proposed.

Finally, the aforementioned points are illustrated by their application to a simple industrial-like test case.

2. Theoretical Remarks

2.1. Necessary Conditions for Non-Linear MIMO Sliding Mode Control

Consider a non-linear MIMO system of the form:

(1)${\overline{x}}^{\left(N\right)}=\overline{f}\left(\overline{x}\right)+[B\left(\overline{x}\right)]\overline{u}$

where:

(2)${\overline{x}}^{\left(N\right)}={\left[x_1^{\left(n_1\right)},\dots ,x_m^{\left(n_m\right)}\right]}^T$

(3)$\overline{u}={\left[u_1,\dots ,u_m\right]}^T$

(4)$\overline{f}\left(\overline{x}\right)={\left[f_1\left(\overline{x}\right),\dots ,f_m\left(\overline{x}\right)\right]}^T$

(5)$[B\left(\overline{x}\right)]=\left[\begin{array}{cccc}{b}_{11}\left(\overline{x}\right)& b_{12}\left(\overline{x}\right)& \cdots & b_{1m}\left(\overline{x}\right)\\ ⋮& ⋮& \ddots & ⋮\\ b_{m1}\left(\overline{x}\right)& b_{m2}\left(\overline{x}\right)& \cdots & b_{mm}\left(\overline{x}\right)\end{array}\right]$

and state vector $\overline{x}$ is defined as:

(6)$\overline{x}={\left[x_1,{\dot{x}}_1,\dots ,x_1^{\left(n_1-1\right)},x_2,{\dot{x}}_2,\dots ,x_2^{\left(n_2-1\right)},\dots ,x_m,{\dot{x}}_m,\dots ,x_m^{\left(n_m-1\right)}\right]}^T$

Assume vector $\overline{f}$ and matrix $\left[B\right]$ are not exactly known, however, their estimates, $\widehat{\overline{f}}$ and $\left[\widehat{B}\right]$, are known, such that the error introduced by the so-called additive uncertainties is bounded by one known function, $F_i$, namely:

(7)$F_i\ge \left|f_i-{\widehat{f}}_i\right|$

while it is possible to link $[\widehat{B}]$ and $\left[B\right]$ as follows:

(8)$[B]=\left([I]+[\Delta ]\right)[\widehat{B}]$

where $\left[\mathsf{\Delta }\right]$ is an unknown matrix whose elements will be indicated with ${\delta }_{ij}$. The sliding variable vector is then defined:

(9)$\overline{s}={\left[\begin{array}{cccc}{s}_1& s_2& \cdots & s_m\end{array}\right]}^T$

where:

(10)$s_i={\left(\frac{d}{dt}+{\lambda }_i\right)}^{n_i-1}{\tilde{x}}_{i}i=1,2,\dots ,m$

In which ${\tilde{x}}_i$ is the ith tracking error. The sliding variables can be rewritten for convenience:

(11)$\overline{s}={\overline{x}}^{\left(N-1\right)}-{\overline{x}}_r^{\left(N-1\right)}$

where:

(12)${\overline{x}}_r^{\left(N-1\right)}={\left[x_{r1}^{\left(n_1-1\right)},x_{r2}^{\left(n_2-1\right)},\dots ,x_{rm}^{\left(n_m-1\right)}\right]}^T$

and:

(13)$x_{ri}^{\left(n_i-1\right)}=x_{di}^{\left(n_i-1\right)}-\left[{\displaystyle \sum _{k=1}^{n_i-1}\left(\begin{array}{c}{n}_i-1\\ k\end{array}\right)}\frac{d^{n_i-k-1}}{d{t}^{n_i-k-1}}{\lambda }_i^k\right]{\tilde{x}}_{i}i=1,2,\dots ,m$

where $x_{di}$ is the ith state desired trajectory. The control problem can be solved by choosing the command laws as [34]:

(14)${\overline{u}}_{SM}={[\widehat{B}]}^{-1}\left({\widehat{\overline{u}}}_{eq}-{\overline{k}}_{sat}\right)$

in which:

(15)${\widehat{\overline{u}}}_{eq}={\overline{x}}_r^{\left(N\right)}-\widehat{\overline{f}}$

(16)${\overline{k}}_{sat}=[\Gamma ]\overline{k}$

and:

(17)$[\Gamma ]=diag\left[sat\left(s_i/{\varphi }_i\right)\right]i=1,\dots ,m$

where ${\varphi }_i$ is the boundary layer introduced to avoid input chattering. Please note that, in order for Equation (14) to exist, $[\widehat{B}]$ must be full-rank. The sliding condition on each $s_i$ is defined as:

(18)$\frac{1}{2}\frac{d}{dt}{s}_i^2\le -{\eta }_i\left|s_i\right|$

Which can be split into:

(19)$\{\begin{array}{l}{\dot{s}}_i\le -{\eta }_{i}for{s}_i\ge 0\\ {\dot{s}}_i\ge {\eta }_{i}for{s}_i<0\end{array}$

The expression of the vector $\dot{\overline{s}}$ is given by:

(20)$\dot{\overline{s}}=\overline{f}+[B][\widehat{B}{]}^{-1}\left({\widehat{\overline{u}}}_{eq}-{\overline{k}}_{sat}\right)-{\overline{x}}_r^{\left(N\right)}$

Thus, Equation (20) becomes the following after a few passages:

(21)$\{\begin{array}{l}\left(1+{\delta }_{ii}\right)k_i\ge {\eta }_i+\left(f_i-{\widehat{f}}_i\right)+{\displaystyle \sum _{j=1}^m{\delta }_{ij}{\widehat{u}}_{jeq}-{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^m{\delta }_{ij}{k}_{j}sat\left(\frac{s_j}{{\varphi }_j}\right)}}{s}_i>{\varphi }_i\\ \left(1+{\delta }_{ii}\right)k_i\ge {\eta }_i-\left(f_i-{\widehat{f}}_i\right)-{\displaystyle \sum _{j=1}^m{\delta }_{ij}{\widehat{u}}_{jeq}+{\displaystyle \sum _{\begin{array}{l}j=1\\ j\ne i\end{array}}^m{\delta }_{ij}{k}_{j}sat\left(\frac{s_j}{{\varphi }_j}\right)}}{s}_i<-{\varphi }_i\end{array}$

These conditions can be verified by solving for $\overline{k}$:

(22)$\{\begin{array}{l}\left([I]-[D]\right)\overline{k}\ge \left(\overline{\eta }+\overline{F}+[D]\left|{\overline{x}}_r^{\left(N\right)}-\widehat{\overline{f}}\right|\right)\\ \overline{k}>\overline{0}\end{array}$

where the inequality is to be verified, component by component, and:

(23)$\overline{k}={\left[k_1,k_2,\dots ,k_m\right]}^T$

(24)${[D]}_{ij}\ge \left|{[\Delta ]}_{ij}\right|$

where ${\left[D\right]}_{ij}$ and ${\left[\mathsf{\Delta }\right]}_{ij}$ are the elements of $\left[D\right]$ and $\left[\mathsf{\Delta }\right]$, respectevely. It follows from Equation (24) that the D-matrix is a majorant matrix of the uncertainties affecting the B-matrix. Please note that the $\overline{k}$, $\overline{F}$, $\widehat{\overline{f}}$ vectors and the D-matrix are, in general, state-dependent. In order to achieve the solution of Equation (22) one has to verify that:

(25)$\{\begin{array}{l}{[D]}_{ii}<1\\ {\lambda }_{max}[D]<1\end{array}$

where ${\lambda }_{max}\left[D\right]$ is the maximum eigenvalue of $\left[D\right]$. Indeed, it is clear from Equations (24) and (25) that matrix $\left[D\right]$ cannot assume any form. In particular, it is possible to obtain Equation (22) from Equation (21) only if the elements on the diagonal of matrix $\left[D\right]$ are less than 1, in addition system (22) admits a solution if matrix $\left(\left[I\right]-\left[D\right]\right)$ is invertible and its inverse matrix has all positive elements.

As a matter of fact, not satisfying Equation (25) implies the impossibility of sliding mode control. A detailed analysis about the origin of the second condition of Equation (25) is reported in the Appendix.

2.2. Procedure for Multiplicative Uncertainty Upper Bound Matrix Definition

In this section, a method to shape the [D] matrix is described. Referring to Equation (8), one has:

(26)$[B]=[\widehat{B}]+[\Delta ][\widehat{B}]$

Thus:

(27)$[\Delta ]=([B]-[\widehat{B}]){[\widehat{B}]}^{-1}$

where one can define the quantity:

(28)${[E]}_{ij}\ge \left|{[B]}_{ij}-{[\widehat{B}]}_{ij}\right|$

As regards additive uncertainties, it is necessary to know the function ${\left[E\right]}_{ij}$ which bind the errors introduced by the multiplicative uncertainties. Remembering Equation (24), after some simple passages, the final form of [D] can be written as:

(29)${[D]}_{ij}=\frac{1}{\left|det[\widehat{B}]\right|}{\displaystyle \sum _{k=1}^m{[E]}_{ik}\cdot \left|adj{[\widehat{B}]}_{kj}\right|}i,j=1,2,\dots ,m$

Therefore, the conditions of Equation (25) can be verified only a posteriori. It is possible to notice from Equation (29) that the entity of the error increase, ${\left[E\right]}_{ik}$, caused by the multiplicative uncertainties, negatively affects the possibility of successfully designing the controller.

2.3. Controller Tuning Method for Systems with First-Order Channel

Consider a non-linear system of the form:

(30)$\left[\begin{array}{l}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{l}{f}_1\left(x_1,x_2\right)\\ f_2\left(x_1,x_2\right)\end{array}\right]+\left[\begin{array}{cc}{b}_{11}\left(x_1,x_2\right)& b_{12}\left(x_1,x_2\right)\\ b_{21}\left(x_1,x_2\right)& b_{22}\left(x_1,x_2\right)\end{array}\right]\left[\begin{array}{l}{u}_1\\ u_2\end{array}\right]$

According to Equation (14), the command laws can be written as:

(31)$u_{1SM}={[\widehat{B}]}_{11}^{-1}\left({\dot{x}}_{1d}-{\widehat{f}}_1-k_{1}sat\left(\frac{s_1}{{\varphi }_1}\right)\right)+{[\widehat{B}]}_{12}^{-1}\left({\dot{x}}_{2d}-{\widehat{f}}_2-k_{2}sat\left(\frac{s_2}{{\varphi }_2}\right)\right)$

(32)$u_{2SM}={[\widehat{B}]}_{21}^{-1}\left({\dot{x}}_{1d}-{\widehat{f}}_1-k_{1}sat\left(\frac{s_1}{{\varphi }_1}\right)\right)+{[\widehat{B}]}_{22}^{-1}\left({\dot{x}}_{2d}-{\widehat{f}}_2-k_{2}sat\left(\frac{s_2}{{\varphi }_2}\right)\right)$

where ${[\widehat{B}]}_{ij}$ are the elements of $[\widehat{B}]$.

It is possible to notice from Equations (31) and (32) that the only parameters available to modify the control laws are $k_1$ and $k_2$. However, these parameters cannot be chosen deliberately, indeed they have to satisfy Equation (22) in order to verify Equation (18). For Equation (30) conditions, Equation (22) can be written as:

(33)$\{\begin{array}{c}\left(1-{\left[\mathrm{D}\right]}_{11}\right)k_1-{\left[\mathrm{D}\right]}_{12}{k}_2\ge h_1\\ -{\left[\mathrm{D}\right]}_{21}{k}_1+\left(1-{\left[\mathrm{D}\right]}_{22}\right)k_2\ge h_2\\ k_1>0\\ k_2>0\end{array}$

where:

(34)$h_1={\eta }_1+F_1+{\left[\mathrm{D}\right]}_{11}\left|{\dot{x}}_{1d}-{\widehat{f}}_1\right|+{\left[\mathrm{D}\right]}_{12}\left|{\dot{x}}_{2d}-{\widehat{f}}_2\right|$

(35)$h_2={\eta }_2+F_2+{\left[\mathrm{D}\right]}_{21}\left|{\dot{x}}_{1d}-{\widehat{f}}_1\right|+{\left[\mathrm{D}\right]}_{22}\left|{\dot{x}}_{2d}-{\widehat{f}}_2\right|$

The solution of Equation (33) can be graphically obtained, as depicted in Figure 1, where one can observe that, in order to satisfy the sliding conditions, $k_1$ and $k_2$ can be deliberately chosen within the dashed area. One usually lets the first inequality of Equation (22) degenerate into an equation, hence the solution is given by the “cone vertex”, i.e., the point ($k_{20}$, $k_{10}$). This solution area allows one to choose the appropriate values for $k_1$ and $k_2$ (within certain boundaries) in order to obtain the desired performance from the system. In addition, $k_1$ and $k_2$ can be “moved” inside the cone depending on the state measurement, further improving the controller performance.

It is clear that both $k_1$ and $k_2$ have to be increased in order to leave the cone vertex. However, minimal “moving” of one of the two coefficients also allows a wider regulation of the other. From an operative point of view, one can approach the problem as follows. The two straight lines which individuate the cone are defined as:

(36)$\{\begin{array}{l}{k}_{1,max}\left(k_2\right)=A{k}_2+B\\ k_{1,min}\left(k_2\right)=a{k}_2+b\end{array}$

where:

(37)$\{\begin{array}{l}A=\frac{1-{\left[D\right]}_{22}}{{\left[D\right]}_{21}}\\ B=-\frac{h_2}{{\left[D\right]}_{21}}\\ a=\frac{{\left[D\right]}_{12}}{1-{\left[D\right]}_{11}}\\ b=\frac{h_1}{1-{\left[D\right]}_{11}}\end{array}$

Hence, for a fixed value of $k_2$, Equation (36) individuates the boundary values $k_1$ can assume. Inside this range, one can choose the desired value for $k_1$ by properly setting the parameter $\epsilon$∈[0; 1], according to:

(38)$k_1\left(k_2,\epsilon \right)=\left[\left(A-a\right)\epsilon +a\right]k_2+\left(b-B\right)\epsilon +b$

That is, the equation of the segment depicted in red in Figure 2. Finally, $k_1$ and $k_2$ can be obtained according to the following relations:

(39)$\{\begin{array}{c}{k}_2=\sigma k_{20}\\ k_1\left(k_2,\epsilon \right)=\left[\left(A-a\right)\epsilon +a\right]k_2+\left(b-B\right)\epsilon +b\end{array}$

With $\sigma >1$. The situation is clarified by Figure 2.

Thus, regulation is achieved by properly choosing parameters σ and ε, which can also be dynamical variable functions. The performance obtained by the $k_1$ and $k_2$ motion inside the cone are compared with the ones obtained by simply setting $k_{10}$ and $k_{20}$ in the test case shown in the next section.

3. Test Case

3.1. System Description and Modeling

Let us consider a system composed of two kilns, schematically represented in Figure 3.

Where $\Delta T_1$ and $\Delta T_1$ are the relative temperatures of the kilns with respect to ambient temperature $T_{amb}$ in K, $R_w$ is the interior kiln thermal resistance in m²K/W while $T_{f1}$ and $T_{f2}$ are the thermal flows generated by the two heat sources in W. The behavior of the system can be simply described as follows:

(40)$\{\begin{array}{l}\frac{T_{f_1}}{2}-T_{f_{ex}}=c{d}_1{V}_1\frac{d\Delta T_1}{dt}+\frac{S_1}{R_W}\Delta T_1\\ \frac{T_{f_1}}{2}+T_{f_2}+T_{f_{ex}}=c{d}_2{V}_2\frac{d\Delta T_2}{dt}+\frac{S_2}{R_W}\Delta T_2\end{array}$

In which c is the air specific heat inside the kilns in J/(kgK), $V_1$, $V_2$, $S_1$ and $S_2$ are the kiln volumes and surfaces, respectively, in m³ and in m², while $d_1$ and $d_2$ are the air densities inside the two kilns in kg/m³. In addition, the term $T_{fex}$ takes into account the exchanged thermal flow between the two kilns according to:

(41)$T_{f_{ex}}=\frac{S_{12}}{2R_W}\left(\Delta T_1-\Delta T_2\right)$

considering $s_{12}$ as the heat exchange surface in m². Furthermore, a non-linear relation between the air density and the temperature inside each kiln is defined:

(42)$d_i=\frac{T_{amb}}{\Delta T_i+T_{amb}}{d}_0$

Finally, the equations of the system written in the normal form are:

(43)$\left[\begin{array}{l}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{c}\frac{\left(T_{amb}+x_1\right)}{a_{1}c}\left(-A{x}_1+C{x}_2\right)\\ \frac{\left(T_{amb}+x_2\right)}{a_{2}c}\left(C{x}_1-B{x}_2\right)\end{array}\right]+\left[\begin{array}{cc}\frac{\left(T_{amb}+x_1\right)}{2a_{1}c}& 0\\ \frac{\left(T_{amb}+x_2\right)}{2a_{2}c}& \frac{\left(T_{amb}+x_2\right)}{a_{2}c}\end{array}\right]\left[\begin{array}{l}{u}_1\\ u_2\end{array}\right]$

where the state vector is defined as:

(44)$\overline{x}={\left[x_1{x}_2\right]}^T={\left[\Delta T_1\Delta T_2\right]}^T$

the inputs are:

(45)$\overline{u}={\left[u_1{u}_2\right]}^T={\left[T_{f_1}{T}_{f_2}\right]}^T$

and the output vector is:

(46)$\overline{y}={\left[y_1{y}_2\right]}^T={\left[x_1{x}_2\right]}^T$

In addition, some constant terms have been collected in:

(47)$a_1=d_0{T}_{amb}{V}_1;a_2=d_0{T}_{amb}{V}_2$

(48)$A=\frac{S_{12}+2S_1}{2R_w};B=\frac{S_{12}+2S_2}{2R_w};C=\frac{S_{12}}{2R_W}$

Hence, Equation (43) describes a MIMO non-linear system in which channels are of the 1st order. Finally, a parametric uncertainty on term c is taken into consideration, supposing its exact (unknown) value to belong to the range:

(49)$c\in \left[c_{min}{c}_{max}\right]$

3.2. Sliding Mode Controller Design

According to the theory shown in Section 2, one can proceed to design the SM controller. The control objectives are the two relative temperatures, $\Delta T_1$ and $\Delta T_2$, corresponding to the state variables of the system. Therefore, the sliding variables can be defined as:

(50)$s_1=x_1-x_{1d}$

(51)$s_2=x_2-x_{2d}$

In order to handle the uncertainties introduced by the not-exactly-known parameter $c$, an estimation of the system functions needs to be defined. In particular defining:

(52)$\widehat{c}=\frac{c_{max}+c_{min}}{2}$

one can write:

(53)${\widehat{f}}_1={f_1|}_{c=\widehat{c}}=\frac{\left(T_{amb}+x_1\right)}{a_1\widehat{c}}\left(-A{x}_1+C{x}_2\right)$

(54)${\widehat{f}}_2={f_2|}_{c=\widehat{c}}=\frac{\left(T_{amb}+x_2\right)}{a_2\widehat{c}}\left(C{x}_1-B{x}_2\right)$

(55)${\widehat{b}}_{11}={b_{11}|}_{c=\widehat{c}}=\frac{\left(T_{amb}+x_1\right)}{2a_1\widehat{c}}$

(56)${\widehat{b}}_{12}={b_{12}|}_{c=\widehat{c}}=0$

(57)${\widehat{b}}_{21}={b_{21}|}_{c=\widehat{c}}=\frac{\left(T_{amb}+x_2\right)}{2a_2\widehat{c}}$

(58)${\widehat{b}}_{22}={b_{22}|}_{c=\widehat{c}}=\frac{\left(T_{amb}+x_2\right)}{a_2\widehat{c}}$

and consequently:

(59)$F_1=\left|\frac{1}{c_{min}}-\frac{1}{\widehat{c}}\right|\left|\frac{T_{amb}+x_1}{a_1}\left(-A{x}_1+C{x}_2\right)\right|$

(60)$F_2=\left|\frac{1}{c_{min}}-\frac{1}{\widehat{c}}\right|\left|\frac{T_{amb}+x_2}{a_2}\left(-B{x}_2+C{x}_1\right)\right|$

(61)$E_{11}=\left|\frac{1}{c_{min}}-\frac{1}{\widehat{c}}\right|\left|\frac{T_{amb}+x_1}{2a_1}\right|$

(62)$E_{12}=0$

(63)$E_{21}=\left|\frac{1}{c_{min}}-\frac{1}{\widehat{c}}\right|\left|\frac{T_{amb}+x_2}{2a_2}\right|$

(64)$E_{22}=\left|\frac{1}{c_{min}}-\frac{1}{\widehat{c}}\right|\left|\frac{T_{amb}+x_2}{a_2}\right|$

According to Equation (29), the D-matrix assumes the form:

(65)$[D]=\left[\begin{array}{cc}\frac{E_{11}}{\left|{\widehat{b}}_{11}\right|}& 0\\ \frac{E_{21}\left|{\widehat{b}}_{22}\right|+E_{22}\left|{\widehat{b}}_{21}\right|}{\left|{\widehat{b}}_{11}{\widehat{b}}_{22}\right|}& \frac{E_{22}}{\left|{\widehat{b}}_{22}\right|}\end{array}\right]$

Finally, one can write the control laws as shown in Equations (31) and (32). In particular, in this case, one has:

(66)$u_{1SM}={[\widehat{B}]}_{11}^{-1}\left({\dot{x}}_{1d}-{\widehat{f}}_1-k_{1}sat\left(\frac{s_1}{{\varphi }_1}\right)\right)$

(67)$u_{2SM}={[\widehat{B}]}_{21}^{-1}\left({\dot{x}}_{1d}-{\widehat{f}}_1-k_{1}sat\left(\frac{s_1}{{\varphi }_1}\right)\right)+{[\widehat{B}]}_{22}^{-1}\left({\dot{x}}_{2d}-{\widehat{f}}_2-k_{2}sat\left(\frac{s_2}{{\varphi }_2}\right)\right)$

where in the first moment, $\overline{k}$ is chosen as equal to the second side of the inequality of Equation (22):

(68)$k_{10}={[I-D]}_{11}^{-1}\left({\eta }_1+F_1+d_{11}\left|{\dot{x}}_{1d}-{\widehat{f}}_1\right|\right)$

(69)$k_{20}={[I-D]}_{21}^{-1}\left({\eta }_1+F_1+d_{21}\left|{\dot{x}}_{1d}-{\widehat{f}}_1\right|\right)+{[I-D]}_{22}^{-1}\left({\eta }_2+F_2+d_{22}\left|{\dot{x}}_{2d}-{\widehat{f}}_2\right|\right)$

3.3. Considerations on the Entity of Uncertainties

In this section, the effect of the error increase on the system, generated by parameter uncertainties, is investigated. The system test case is quite simple, and it is possible to make some considerations starting from the first condition of Equation (25). Looking at Equation (65) one has:

(70)$d_{11}=\frac{E_{11}}{\left|{\widehat{b}}_{11}\right|}<1$

which, remembering Equations (55) and (61), becomes:

(71)$\frac{\left|\frac{1}{c_{min}}-\frac{1}{\widehat{c}}\right|\left|\frac{\left(T_{amb}+x_1\right)}{2a_1}\right|}{\left|\frac{\left(T_{amb}+x_1\right)}{2a_1\widehat{c}}\right|}<1$

Finally, considering $I=c_{max}-c_{min}$ and after a few simple passages, one can obtain:

(72)$\frac{I}{\widehat{c}}<1$

This relation means that a limit on the uncertainty range extension exists, beyond which SM control is not feasible. The same analysis can be carried out for $d_{22}$, obtaining the same result. The overall control scheme is depicted in Figure 4.

The system model and the controllers were both implemented in the MATLAB and Simulink environments, as it can be appreciated from Figure 5:

3.4. Simulation Set I

This first set of simulations aims at assessing the correct SM controller synthesis and its robustness with respect to parametric and modelling uncertainties. To this extent, together with the uncertainty on the knowledge of parameter $c$, two actuators (neglected in the controller synthesis) have been inserted in the system, characterized by a time constant ${\tau }_a$ representing the delay time between the instant in which the controller orders the suitable values for $T_{f_1}$ and $T_{f_2}$ (labeled in Figure 4 as $T_{f_1}^{*}$ and $T_{f_2}^{*}$) and the one in which the two heat sources actually generate them.

The numerical values of the constant parameters used in the simulations are reported in Table 1.

As previously mentioned, the value of parameter c is not exactly known, thus it is necessary to provide its estimation and range of belonging to the SM controller. In view of what has been illustrated in the previous section, to guarantee the control feasibility, these values are chosen as:

(73)$c_{min}=772\left[\mathrm{J}/\mathrm{kg}\cdot \mathrm{K}\right];c_{max}=1127\left[\mathrm{J}/\mathrm{kg}\cdot \mathrm{K}\right]$

while the value of $\widehat{c}$ derives from Equation (52). In this situation one has:

(74)$d_{11}=d_{22}={\lambda }_{max}=0.23$

Thus the conditions of Equation (25) are verified.

The simulation, as shown in Figure 6, considers a step variation in the $\Delta T_1$ reference 10 s after the simulation start, from an initial condition of 500 °C to a steady state working point at 650 °C, while $\Delta T_2$ reference is kept constant at 700 °C. Please remember that in this first set of simulations, $k_1=k_{10}$ and $k_2=k_{20}$.

The simulation results show the correct operation of the SM controller, which is able to take the controlled variables as their references. However, from Figure 7, one can notice a non-negligible dip in the $\Delta T_2$ profile during the transient. Figure 8 shows the time profiles of the two thermal flows coming from the two actuators.

It is interesting to investigate the system behavior with the increasing of the entity of uncertainties. By choosing a wider uncertainty range of parameter c, such that Equation (72) is violated (e.g., $c_{min}=472$ J/kgK; $c_{max}=1427$ J/kgK), the system becomes unstable, as shown in Figure 9.

Indeed, the conditions necessary for SM control feasibility are not complied with, as now:

(75)$d_{11}=d_{22}={\lambda }_{max}=1.01$

3.5. Simulation Set II

In this section, the controller tuning proposed in Section 3 is applied in order to improve the performance of the system. In particular, in the following attention is focused on keeping $\Delta T_2$ as constant as possible throughout the entire transient, by dynamically regulating the terms $k_1$ and $k_2$. For this specific test case, as $D_{12}=0$, the cone assumes the shape in Figure 10.

As one can intuitively understand, an increase or decrease of $k_i$ involves a greater or smaller “strength” of the ith channel. As a consequence, in order to pursue the aforementioned goal, it is necessary to make the 2nd channel of the system “stronger” than the 1st. Therefore, during the transient, $k_1$ must be kept at its minimum value, while $k_2$ has to be dynamically increased within the dashed area in Figure 10. This means that regulation parameter ε must be zero value throughout the transient (which implies that $k_1=k_{1,min}=k_{10}$), while $\sigma$ is defined on the basis of the tracking error of $\Delta T_2$ such that, if $e_2$ increases, $k_1$ must also increase accordingly, e.g., the functional relation depicted in Figure 11.

The comparison between the performance obtained, with and without the proposed controller tuning, are presented in Figure 12, Figure 13, Figure 14, Figure 15 and Figure 16. The variable profiles obtained through the proposed tuning method are reported by red dash-dot lines, while the ones related to the conventionally tuned control system are reported by continuous blue lines.

The results of this second simulation set clearly point out the performance improvement obtained thanks to the tuning method proposed. Indeed, with an almost identical behavior of the variable $\Delta T_1$ (Figure 12), one can notice in Figure 13 how the regulation leads to a much better profile for $\Delta T_2$, which does not deviate more than 5 °C from the constant reference.

Figure 14 highlights that the thermal flow injected by the first heat source is substantially the same with and without the proposed controller tuning approach, while the second is forced by the controller to produce a greater amount of heat in order to prevent the second kiln temperature from decreasing too much.

Finally, Figure 15 and Figure 16 show the regulation profiles of $k_1$ and $k_2$ compared to the limit values they can dynamically assume, highlighting that, while $k_1$ is always set at its minimum value, $k_2$ is increased during the transient to keep the corresponding temperature as constant as possible.

4. Conclusions

Some analytical remarks on the SM control theory for MIMO systems have been reported in this paper. In particular, the conditions necessary for controller feasibility on the entity of uncertainties have been thoroughly investigated. In addition, a novel controller tuning method has been proposed in order to successfully obtain the desired performance, even when the typical SM regulation parameters are not available. More precisely, it is obtained that first order SM controllers for MIMO systems can be designed as long as the so-called multiplicative uncertainties are sufficiently contained; additionally, a procedure able to exploit the uncertainty compensating terms ($k_i$) both for satisfying the sliding condition and guaranteeing high controller performance is outlined. Finally, the validity of the theoretical considerations is proved through a simple test case. In order to show the effect of increasing entities of uncertainty and to illustrate the effective performance improvement obtained thanks to the proposed tuning method, a MIMO SM controller is designed and applied to the temperature regulation problem of two adjacent industrial kilns. Future developments will consider the implementation and validation of the proposed approach in a real test case or in a real time simulation environment in hardware in a loop configuration.

Author Contributions

A.P. and R.P. studied and developed the proposed tuning criteria, A.B. and M.B. provided its implementation and validation in the considered test case. A.M. supported the software implementation. I.M. and B.V. supervised the work and collaborated in the writing of the paper, together with A.M.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

In mathematics, the class of Z-matrices are those matrices whose off-diagonal entries are less than or equal to zero; that is, a Z-matrix A satisfies:

(A1)

$A=\left[a_{ij}\right]with{a}_{ij}\le 0foralli\ne j$

Definition A1.

Let A be a n × n Z-matrix. That is $A=\left(a_{ij}\right)$ where $a_{ij}\le 0$ for all $i\ne j$, $1\le i,j\le n.$ Then matrix A is also a M-matrix if it can be expressed in the form:

(A2)

$A=sI-D$

where $D=\left(d_{ij}\right)$ with $d_{ij}\ge 0$, for all $1\le i,j\le n$, s is greater than the maximum of the module of the eigenvalues of D and I is the identity matrix.

Definition A2.

Below, $\ge$ denotes the element-wise-order. That is, for any real matrices A, B of size m × n, we write $A\ge B$ if $a_{ij}\ge b_{ij}$ for all $i,j$.

Proposition A1.

Let $D=\left(d_{ij}\right)$ be such that $d_{ij}\ge 0$, for all $1\le i,j\le n$ with maximum eigenvalue ${\lambda }_{max}\le 1$, then $I-D$ is a M-matrix.

Proof. 

Obviously each off-diagonal term is less than or equal to 0, moreover, by choosing $s=1$ we get the definition. □

Proposition A2.

If A is a M-matrix then A is inverse-positive. That is, $A^{-1}$ exists and $A^{-1}\ge 0$.

Proof. 

It can be found in Reference [39]. □

Figures and Table

Figure 1. Graphical solution for k1 and k2.

Figure 2. Regulation of k1 and k2 coefficients inside the cone.

Figure 3. Test case system scheme.

Figure 4. Control system scheme.

Figure 5. MATLAB and Simulink implementation.

Figure 6. ΔT1 time profile.

Figure 7. ΔT2 time profile.

Figure 8. Profiles of the thermal flows injected in the two kilns.

Figure 9. Unstable behavior of the system due to large uncertainties.

Figure 10. Graphical solution for k1 and k2.

Figure 11. Regulation function of k2.

Figure 12. Performance comparison—ΔT1.

Figure 13. Performance comparison—ΔT2.

Figure 14. Performance comparison—injected thermal flows.

Figure 15. Performance comparison—k1 profile.

Figure 16. Performance comparison—k2 profile.

Table 1

Parameters definition.

Symbol	Value	Symbol	Value
$V_1$	8 m3	$V_2$	15.6 m3
$S_1$	16 m2	$S_2$	25 m2
$S_{12}$	4 m2	$T_{amb}$	15 °C
$d_0$	1.225 kg/m3	$c_{real}$	1005 J/kg·K
$R_W$	0.2 m2·K/W	${\tau }_a$	5 s