Combining Augmented Error Modeling Technique and Block-Pulse Functions Method for Tracking Control Design of Nonlinear Polynomial Systems with Multiple Time-Delayed States Subject to Nonsymmetric Input Saturation

Abstract

The present research work is intended to synthesize a novel tracking control strategy for a class of nonlinear polynomial systems characterized by multiple well-defined delays in state variables under the presence of nonsymmetric input saturation. The design strategy makes full use of an associate’s memory nonlinear state feedback control with integral-based actions. An original control scheme joining block-pulse functions method combined with the augmented error modeling technique is used to infer the controller’s tracking gains. The objective is to convert the investigated nonlinear algebraic problem governed by specifying constraints into a constrained linear one that can be solved in the constrained least square methodology. Detailed novel sufficient conditions proving the closed-loop augmented system’s practical stability are elaborated. The instance of a twofold inverted pendulums benchmark is considered so as to exhibit the benefits of the proposed control approach.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The design of an output tracking control scheme for continuous-time nonlinear systems subject to time delays is one of the most investigated problems in nonlinear control research theory [1]. This topic is not a challenging one only because this class of time-delayed systems is highly confronted in the industrial process, but mainly due to the restrictions in proving the desired closed-loop performance objectives, in addition to the stability requirement. As needs be, the delicacy identified with the tracking control problem is intensively considered by researchers nowadays [2–5].

In the framework of the Lyapunov theory, there are mainly two fundamental methodologies to synthesize a tracking control law for a particular class of nonlinear time-delayed systems. The first one is based on the Lyapunov–Krasovskii functionals, while the subsequent one exploits the Lyapunov–Razumikhin functionals. Both methodologies can provide an asymptotic tracking performance. Based on these Lyapunov functionals, most kinds of tracking control methods for various classes of nonlinear time-delayed systems have been widely developed, such as sliding mode control [6], model predictive control [7], fault-tolerant control [8], backstepping control [9, 10], dynamic surface control [11], $H_{\infty }$ control [12], and adaptive fuzzy or network control [13, 14].

The significant drawbacks of those strategies originate from the fact that no generic algorithms exist to define a Lyapunov–Krasovskii/Razumikhin function candidate, as well as from the accompanied computational complexity of the existing tracking control methods for providing the controller’s gains.

It is worth mentioning that all of the aforementioned control approaches are applicable for particular classes of nonlinear systems under the severe Lipschitz assumption and other additional constraints [6, 7], just to mention a few, cases of nonlinear systems with time delays designated by a particular differentiable form [10], cases of nonlinear systems with time delays represented in pure feedback form [9] or a lower-triangular form [11] or strict-feedback form [8], and cases of nonlinear systems with time delays modeled by advanced technique [12–14].

On the other hand, it is recognized that the actuator saturation is a potential problem of all practical dynamic systems [15]. So, every control system has to face an input saturation, which can result in performance abasement such as oscillations, high overshoot, undershoot, lag, large settling time, and sometimes even instability. Thus, both actuator saturation and time delays are frequently encountered when modeling an engineering system. To design appropriate control laws for dealing with the presence of input saturation, three main strategies have been presented in the literature for continuous-time nonlinear systems affected by time delays, which can be cited as follows:

(a) The anti-windup compensation-based methods, in which an auxiliary design system is introduced to compensate the effect of input saturation. There are two schemes in this case, called one-step design scheme or two-step design scheme. For the last one, first a controller is designed without regarding control input nonlinearity, and then, an anti-windup compensator will be augmented in order to minimize the undesirable degradation of closed-loop performance caused by input saturation. The one-step design scheme deals with saturation effect by simultaneously designing the controller and its associated anti-windup compensation. Recently, a few exceptional works have been developed in which the tracking control design problems have been investigated. In [16], a full-order dynamic internal model control-based anti-windup compensation architecture was proposed for stable nonlinear time-delayed systems, satisfying the Lipschitz condition, under actuator saturation. A local decoupling anti-windup compensation design was suggested using the delay-range-dependent approach, for compensation of the undesirable saturation effects. In [17], by employing the concept of decoupled architecture-based design, a dynamic anti-windup compensation applicable to both stable or unstable nonlinear time-delayed systems was proposed, by reformulating the Lipschitz nonlinearities of the considered class of system via linear parameter-varying models. In [18], a novel approach was proposed for designing a static anti-windup compensation for two classes of delayed nonlinear systems under input actuator saturation. Based on the Wirtinger-based inequality, a delay-range-dependent technique was exploited to derive a condition for calculating the anti-windup compensation gains. By using the quadratic Lyapunov–Krasovskii functional, Wirtinger-based inequality, sector conditions, bound on delays, time-varying range of delay, and $H_{\infty }$ gain reduction, several conditions in terms of LMIs are derived to guarantee local and global stabilization and performance of the overall closed-loop system. It should be mentioned that all of these methods can only handle the nonlinear time-delayed systems under the severe Lipschitz condition.

(b) The convexity-based methods, in which the input constraints are considered at the controller design stage in order to adjust the controller gains according to the saturation levels. In this case, there are two approaches. The first one is to use a polytopic representation of the saturation function. Based on this representation and using the Takagi–Sugeno modeling, some interesting works were devoted only to the stabilization control problem for saturated T-S fuzzy time-delayed systems by means of delay-dependent Lyapunov–Krasovskii functional [19] or delay-independent Lyapunov–Krasovskii functional [20]. Another way to take the saturation effects into account relying on the sector nonlinearity approach is adopted to derive Popov and Circle criteria [21]. Lately, this approach was effectively applied to tackle the stabilization control problem for Lipschitz nonlinear time-delayed systems subject to symmetric input saturation [22], as well as the tracking control problem for lower-triangular nonlinear time-delayed systems subject to nonsymmetric input saturation [23, 24]. It is important to highlight that these methods, despite that they are restricted to specific classes of nonlinear time-delayed systems, provide an attractive region of initial states that ensure the asymptotic stability of the considered closed-loop system.

(c) The approximation-based intelligent control methods, in which the nonsmooth nonlinear input saturation is approximated by a smooth function. Generally, two commonly used saturated functions are the logistic sigmoid and the hyperbolic tangent. By using neural network or fuzzy approximation technique to relax the growth assumptions and employing Lyapunov–Krasovskii functionals to deal with time delay, an increasing attention has been paid to the adaptive tracking control design for saturated nonlinear time-delayed systems in the triangular form [25, 26] or strict-feedback form [27, 28]. A main limitation of the mentioned works is that the required number of neural-network rules or fuzzy rules to achieve given accuracy leads to increasing complexities of the procedure of computing tracking control law.

Up to our knowledge, the time-varying set point tracking control problem for a particular class of multi-input multi-output nonlinear time-delayed systems, namely nonlinear polynomial system with multiple time-delayed states, subject to nonsymmetric input saturation is not yet investigated until now. It must be noted that such polynomial model is considered as a single unified form to represent the nonlinear dynamics of analytical systems. In fact, any analytical nonlinear model can be represented by a polynomial one utilizing the Kronecker tonsorial power [29] and the Taylor series expansion as main powerful analytical tools. Therefore, this kind of model can accurately describe the dynamical behavior of a large set of physical processes and real industrial plants, including but not limited to the following: electrical machines [30], robot manipulators [31], power systems [32], and interconnected inverted pendulums [33]. This advantage is justified by its capacity (i) to preserve the high-order model as well as its original nonlinearity, which leads to an accurate controller’s gains when computing the control law, and (ii) to enlarge the operating region approximation around the equilibrium point.

Due to the necessity of bringing remedy to the drawbacks of the exiting control methods based on the Lyapunov theory, increasing interests are showed up as of late toward the improvement of numerical approximation methods based on the modeling concept of nonlinear polynomial systems. Some works are devoted to the tracking control problem for such undelayed system in the case of a step input [31] as well as in the case of a time-varying set point input [34]. For the subclass of the nonlinear polynomial system under the presence of a single time-delayed state, whose value is supposed to be constant and known, the tracking control problem has been addressed in [35] using block-pulse functions as a tool of approximation. Nevertheless, it is notable that up to now, none published research in which any one of the control design problems for the nonlinear polynomial system with or without time delays under input saturation is studied.

Motivated by the two concerns: firstly, to tackle the output tracking control problem for the nonlinear polynomial system with multiple time-delayed states subject to nonsymmetric input saturation; secondly, to surmount the drawbacks of previous studies building through a Lyapunov–Krasovskii or Razumikhin function approach, whose cannot tackle the mentioned class of nonlinear time-delayed systems, then it might be beneficial to exploit simultaneously the augmented error modeling technique and numerical method tools. This is justified by the fact that the static errors are reduced effectively by increasing the state vector through a new output tracking error vector, as well as by the conversion of the posed saturated tracking control problem into solving a constrained linear system of algebraic equations, depending on the parameters of the feedback regulator.

Proceeding from the above idea, a new algebraic scheme for the development of a memory polynomial state feedback controller with an integral-action-based compensator for continuous-time nonlinear polynomial system affected by multiple time delays in the states and under nonsymmetric input saturation is presented in the current study. Right off the bat, the modeling technique of augmented error is performed to determine an equivalent augmented form for the closed-loop controlled system, which permits to achieve high trajectory tracking performance [34, 36]. Then, the block-pulse functions basis as an approximation tool as well as their operational matrices are exploited to transform the augmented state-space model of controlled saturated process into an equivalent constrained linear model of algebraic equations, which depends only on the controller gains. Thus, the proposed computation method in this work is simple as well as swift.

Walsh functions [37], block-pulse functions [38], and Haar wavelet functions [39] are considered as the most impressive types of piecewise constant basis functions. Contrary to them, block-pulse functions are featured by their simple structure and direct implementation [40]. This is generally because of their significant properties, that is, labeled disjointedness, fulfillment, and orthogonality. Thus, they have been exploited widely as a valuable tool in the synthesis, analysis, estimation, identification, and different problems related to the area of control and systems [41–43].

The main purpose of this work is to substitute both coefficients of delayed and undelayed augmented state vectors by those of reference models issued from projecting the state-space model over a basis of block-pulse functions. By applying the operational matrices of the considered basis accompanied by the use of explicit arithmetical control operations, it comes out a constrained linear system of algebraic equations, which can be solved by the means of constrained least squares optimization. It is worth noting that the insertion of the imposed limits by the actuator saturation on the control input is translated here through the expansion of the saturated input control over the considered block-pulse functions basis. As an original result, novel sufficient conditions are determined for practical stability guarantees of the augmented closed-loop system under the presence of input saturation constraint. Consequently, the significant contribution of this study can be summed up as follows:

(1) This work is addressing for the first time, the output trajectory tracking control problem for the class of nonlinear polynomial systems with multiple time-delayed states, subject to nonsymmetric input saturation. For this reason, a new combined block-pulse functions method joined to the augmented error modeling technique is applied. This scheme helps in designing memory polynomial state feedback controllers with integral actions.

(2) To our knowledge, the proposed method for the introduction of the imposed limits on the control input by the actuator saturation in the procedure of computing tracking control law for the addressed class of saturated nonlinear time-delayed systems is developed for the first time in this study. This method has the ability to treat both nonsymmetric or symmetric input saturation, which can be evaluated as an important advantage.

(3) Dissimilar to the based Lyapunov–Krasovskii/Razumikhin function tracking method and its related computational intricacy, the proposed controller’s gains in this work are derived by means of constrained least squares optimization. This is established throughout algebraic operations where sufficient stability conditions of the controlled augmented system are defined explicitly. These conditions are defined without precedent for this study.

Following this introduction, a short description of block-pulse functions and their fundamentals are presented. Next, the closed-loop system model is described leading to the problem statement. The essential outcomes are given in the subsequent section. The efficiency and performance of the developed strategy are evaluated on a twofold inverted pendulums benchmark.

Notation. During the whole of this study, the $n$ dimensional Euclidean space is denoted by ${ℝ}^n$, while ${ℝ}^{n\times m}$ is referred to the set of all real matrices with $n$ rows and $m$ columns, whereas ${ℝ}^{n\times m{n}^{i+j}}$ stands for the set of all real matrices with $n$ rows and $m\times n^{i+j}$ columns, where $i$ and $j$ are two natural numbers. $I_n$ denotes the identity matrix of size $n\times n$, whereas $0_{n\times m}$ denotes the zero matrix of size $n\times m$. $A^T$ corresponds to the transpose of the matrix $A$. $e_i^p$ denotes the $p$ dimensional unit vector which has 1 in the $i^{th}$ element and zero elsewhere. The adopted vector norm is the Euclidean norm, and the matrix norm is the corresponding induced norm.

2. Block-Pulse Functions and their Properties

2.1. Block-Pulse Functions

Over the interval $0,T$, an $N$ set of block-pulse functions (BPFs) can be defined as follows [38]:$$\begin{array}{}\text{(1)}& {\phi }_{i}t=\begin{array}{cc}1,& \frac{iT}{N}\le t\le \frac{i+1T}{N},\hspace{1em}i=0,\dots ,N-1,\\ 0,& \text{else\hspace{0.17em}where,}\end{array}\end{array}$$where for each $i\in 0,\dots ,N-1$, ${\phi }_{i}t$ represents the i-th block-pulse functions.

The principal attributes of block-pulse functions are as follows.

2.1.1. Disjointness

The block-pulse functions are disjoint with each other in the interval $0,T$:$$\begin{array}{}\text{(2)}& \forall i,j\in 0,\dots ,N-1,{\phi }_{i}t{\phi }_{j}t=\begin{array}{cc}{\phi }_{i}t& fori=j\\ 0& fori\ne j\end{array}.\end{array}$$

2.1.2. Orthogonality

The block-pulse functions are orthogonal with each other in the interval $0,T$:$$\begin{array}{}\text{(3)}& \forall i,j\in 0,\dots ,N-1,{\displaystyle {\int }_0^T{\phi }_{i}t{\phi }_{j}tdt}=\begin{array}{cc}\frac{T}{N}& fori=j\\ 0& fori\ne j\end{array}.\end{array}$$

The orthogonal property of block-pulse functions is the basis of expanding functions into their block-pulse series. Thus, any absolutely integrable function $ft$ of the Lebesgue measure in the interval $0,T$ may be expanded over a BPFs basis as follows:$$\begin{array}{}\text{(4)}& ft\cong {\displaystyle \sum _{i=0}^{N-1}{f}_i{\phi }_{i}t}=F_N{S}_{N}t,\end{array}$$where the block-pulse coefficient vector $F_N$ is given by$$\begin{array}{}\text{(5)}& F_N=f_0\cdots f_{N-1},\end{array}$$and the block-pulse vector basis of order $N$ is given by$$\begin{array}{}\text{(6)}& S_{N}t={{\phi }_{0}t\cdots {\phi }_{N-1}t}^T,\end{array}$$with $f_i$ is the block-pulse coefficient computing as follows:$$\begin{array}{}\text{(7)}& \forall i\in 0,\dots ,N-1,f_i=\frac{N}{T}{\displaystyle {\int }_0^{T}ft{\phi }_{i}tdt}=\frac{N}{T}{\displaystyle {\int }_{iT/N}^{i+1T/N}ftdt}.\end{array}$$

2.1.3. Completeness

The block-pulse function set is complete when the order $N$ approaches infinity. This means that for any square-integrable function $ft$ of the Lebesgue measure in the interval $0,T$, we have$$\begin{array}{}\text{(8)}& {\displaystyle {\int }_0^T{f}^{2}tdt}=\underset{N⟶+\infty }{\lim }{\displaystyle \sum _{i=0}^{N-1}{f}_i^2{\displaystyle {\int }_0^T{\phi }_i^{2}tdt}}.\end{array}$$

2.2. Error in BPFs Approximation

Assuming that $ft$ is a differentiable function, where his first derivative on the time interval $0,T$ is bounded, that is to say$$\begin{array}{}\text{(9)}& \exists M\succ 0,such{f}^{\prime }t\le M.\end{array}$$

Over every subinterval $iT/N,i+1T/N$, the residual error is given by$$\begin{array}{}\text{(10)}& e_{i}t=f_i{\phi }_{i}t-ft=f_i-ft=\frac{N}{T}{\displaystyle {\int }_{iT/N}^{i+1T/N}f\tau d\tau }-ft.\end{array}$$

Furthermore, the total error can be defined as follows:$$\begin{array}{}\text{(11)}& et={\displaystyle \sum _{i=0}^{N-1}{e}_{i}t},\hspace{1em}t\in 0,T,\end{array}$$which can be bounded as follows [43]:$$\begin{array}{}\text{(12)}& et\le \frac{M}{N}\sqrt{\frac{T^3}{12}},\end{array}$$allowing to conclude that$$\begin{array}{}\text{(13)}& \underset{N⟶+\infty }{\lim }et=0.\end{array}$$

Consequently, over a BPFs basis with a high order $N$, we can represent exactly the function $ft$.

2.3. BPFs Expansion of a Vector Function

A vector function $xt$ of $n$ dimensional components which are square-integrable with respect to Lebesgue measure in the interval $0,T$ can be expressed practically by a finite series of block-pulse functions as follows:$$\begin{array}{}\text{(14)}& xt\cong {\displaystyle \sum _{i=0}^{N-1}{x}_i{\phi }_{i}t}=X_{1,N}{S}_{N}t,\end{array}$$with $X_{1,N}=x_0,\dots ,x_{N-1},$ where $x_i$ are the block-pulse coefficients of vector function $xt$, as computed through the orthogonality property of the block-pulse functions as follows:$$\begin{array}{}\text{(15)}& x_i=\frac{N}{T}{\displaystyle {\int }_{\frac{iT}{N}}^{\frac{i+1T}{N}}xtdt.}\end{array}$$

2.4. BPFs Expansion of a Time-Delayed Vector Function

Assume that the vector $xt$ has its initial as$$\begin{array}{}\text{(16)}& xt=\zeta t\hspace{1em}-\tau \le t\le 0,\end{array}$$with $0<\tau <T$. Then, for $t\in 0,T$, the time-delayed vector $xt-\tau$ is presented by$$\begin{array}{}\text{(17)}& xt-\tau =\begin{array}{c}\zeta t-\tau \hspace{1em}\hspace{1em}for\hspace{1em}0\le t\le \tau \hfill \\ xt-\tau \hspace{1em}\hspace{1em}for\hspace{1em}\tau <t\le T\hfill \end{array}.\end{array}$$

As a result, the expansion of time-delayed vector $xt-\tau$ over the BPFs basis is given by [44]$$\begin{array}{}\text{(18)}& xt-\tau \cong {\displaystyle \sum _{i=0}^{N-1}{x}_i^{\ast }\tau {\phi }_{i}t}=X_{1,N}^{\ast }\tau S_{N}t,\end{array}$$with$$\begin{array}{}\text{(19)}& X_{1,N}^{\ast }\tau =x_0^{\ast }\tau ,\dots x_{N-1}^{\ast }\tau ,\text{(20)}& x_i^{\ast }\tau =\frac{N}{T}{\displaystyle {\int }_{iT/N}^{i+1T/N}xt-\tau \mathrm{d}t}=\begin{array}{cc}{\zeta }_i^{\ast }\tau & \mathrm{for}\hspace{1em}i\le \mu \\ x_{i-\mu }& \mathrm{for}\hspace{1em}i>\mu \end{array},\end{array}$$where $\mu$ is the number of BPFs considered over $0\le t\le \tau$, and$$\begin{array}{}\text{(21)}& {\zeta }_i^{\ast }\tau =\frac{N}{T}{\displaystyle {\int }_{iT/N}^{i+1T/N}\zeta t-\tau \mathrm{d}t}.\end{array}$$

2.5. Operational Matrix of Integration

For the BPFs basis, an operational matrix of integration is defined in [44] by$$\begin{array}{}\text{(22)}& {\displaystyle \underset{0}{\overset{t}{\int }}{S}_{N}tdt=}{P}_N{S}_{N}t,\end{array}$$where $P_N=T/2N\begin{array}{cccc}1& 2& \dots & 2\\ 0& 1& \dots & 2\\ ⋮& \ddots & \ddots & 2\\ 0& \dots & 0& 1\end{array}.$

2.6. Operational Matrix of Product

The elementary matrix can be defined as follows:$$\begin{array}{}\text{(23)}& E_{i,j}^{p\times q}=e_i^p\otimes e_j^{q^T},\end{array}$$where $\otimes$ is the symbol of the Kronecker product [29].

Based on the disjointness property of BPFs, the operational matrix of product is given by [35]$$\begin{array}{}\text{(24)}& S_{N}t\otimes S_{N}t=\begin{array}{c}{E}_{\mathrm{1,1}}^{N\times N}\\ ⋮\\ E_{N,N}^{N\times N}\end{array}{S}_{N}t=M_N{S}_{N}t.\end{array}$$

3. Problem Formulation and System Description

Consider a continuous-time nonlinear polynomial systems with known and fixed multiple time delays in states $S$, which are described by the following state equation:$$\begin{array}{}\text{(25)}& S\begin{array}{c}\dot{x}t={\displaystyle \sum _{i=1}^{q_1}{A}_i{x}^{i}t}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{D}_{jl}{x}^{j}t-{\tau }_l}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{G}_{ijl}{x}^{i}t\otimes x^{j}t-{\tau }_l}}}\\ +{\displaystyle \sum _{i=1}^{q_2}{B}_{i}ut\otimes x^{i}t}+{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{C}_{jl}ut\otimes x^{j}t-{\tau }_l}}\\ +{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{Q}_{ijl}ut\otimes x^{i}t\otimes x^{j}t-{\tau }_l}}}+But\\ yt=Cxt\\ xt=\zeta t\hspace{1em}for\hspace{1em}t\in -\overline{\tau },0\end{array},\end{array}$$where $ut={u_{1}t,\dots ,u_{m}t}^T\in {ℝ}^m$ is the input vector, $xt\in {ℝ}^n$ is the nondelayed state vector, $xt-{\tau }_l\in {ℝ}^n$ is the delayed state vector with ${\tau }_l$ denoted time delay, and $yt\in {ℝ}^p$ is the output vector. The continuous vector-valued function $\zeta t$ denotes the initial data, where $\overline{\tau }=\max {\tau }_l$ for each $l\in 1\dots v$, with $v=\max v_i$ for each $i\in \mathrm{1,2}$.

In (25), $A_i\in {ℝ}^{n\times n^i}$, $D_{jl}\in {ℝ}^{n\times n^j}$, $G_{ijl}\in {ℝ}^{n\times n^{i+j}}$, $B_i\in {ℝ}^{n\times m{n}^i}$, $C_{jl}\in {ℝ}^{n\times m{n}^j}$, $Q_{ijl}\in {ℝ}^{n\times m{n}^{i+j}}$, $B\in {ℝ}^{n\times m}$, $C\in {ℝ}^{p\times n}$, and $x^i$ is the $i^{th}$ Kronecker power of the state vector $x$, defined by $\begin{array}{c}{x}^0=1,\hfill \\ x^{i+1}=x^i\otimes x\hspace{1em}\forall i>0.\hfill \end{array}$

3.1. Assumption

The consider system (25) is locally controllable around $x_0=\zeta 0$ where its $n$ state components are all physically measurable.

3.2. Recommended Strategy

Control objective involves the tracking of outputs of saturated nonlinear time delays system for given reference input vector $y_{c}t\in {ℝ}^p$. For this purpose, we intend to design a memory polynomial state feedback controller accompanied by compensator gains:$$\begin{array}{}\text{(26)}& ut=\overline{N}Rt-{\displaystyle \sum _{i=1}^{q_1}{K}_i{x}^{i}t}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{L}_{jl}{x}^{j}t-{\tau }_l}},\end{array}$$with$$\begin{array}{}\text{(27)}& Rt={e{t}^T,e^1{t}^T,\dots e^{w-1}{t}^T}^T,\end{array}$$where$$\begin{array}{}\text{(28)}& et=\underset{wtimes}{\underbrace{{\displaystyle {\int }_0^t\dots }{{\displaystyle \int }}_0^t}}{y}_c\sigma -y\sigma d{\sigma }^w,\end{array}$$with $w$ is a positive integer which depends on the type of reference input vector. In (26), $\overline{N}\in {ℝ}^{m\times wp}$, $K_i\in {ℝ}^{m\times n^i}$, $L_{jl}\in {ℝ}^{m\times n^j}$, and $e^{h}t$ are the $h-th$ derivative of $et$.

Each element of the actuator vector $u_{i}t$ should satisfy the following constraint:$$\begin{array}{}\text{(29)}& -{\overline{v}}_i\le u_{i}t\le {\overline{u}}_{i}with0\prec {\overline{u}}_i\text{\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}}0\prec {\overline{v}}_i,\end{array}$$which is equivalent to$$\begin{array}{}\text{(30)}& -v_{\max }\le ut\le u_{\max },\end{array}$$where $u_{\max }={{\overline{u}}_1,\dots ,{\overline{u}}_m}^T$ and $v_{\max }={{\overline{v}}_1,\dots ,{\overline{v}}_m}^T$ are given positive reals vectors.

The constrained inputs, denoted $sat{u}_{i}t$, are saturating functions and defined for each $i\in 1\dots m$ as follows:$$\begin{array}{}\text{(31)}& sat{u}_{i}t=\begin{array}{c}{\overline{u}}_i\hspace{1em}\hspace{1em}\hspace{1em}if\hspace{1em}{u}_{i}t\succ {\overline{u}}_i\hfill \\ u_{i}t\hspace{1em}if\hspace{1em}-{\overline{v}}_i\le u_{i}t\le {\overline{u}}_i\hfill \\ -{\overline{v}}_i\hspace{1em}\hspace{1em}if\hspace{1em}{u}_{i}t\prec -{\overline{v}}_i\hfill \end{array},\end{array}$$with$$\begin{array}{}\text{(32)}& satut={\begin{array}{ccc}sat{u}_{1}t& \dots & sat{u}_{m}t\end{array}}^T.\end{array}$$

3.3. Control Objective

The undelayed augmented state vector $Xt$ as well as the delayed augmented state vector $Xt-{\tau }_l$ can be expressed using the augmented error modeling technique, as follows:$$\begin{array}{}\text{(33)}& Xt={x{t}^T,R{t}^T}^T,\text{(34)}& Xt-{\tau }_l={x{t-{\tau }_l}^T,R{t-{\tau }_l}^T}^T.\end{array}$$

When taking into consideration that$$\begin{array}{}\text{(35)}& xt=\begin{array}{cc}{I}_n& 0_{n,wp}\end{array}Xt=\Psi Xt,\text{(36)}& x^{i}t={\Psi }^i{X}^{i}t,\text{(37)}& xt-{\tau }_l=\Psi Xt-{\tau }_l,\text{(38)}& x^{i}t-{\tau }_l={\Psi }^i{X}^{i}t-{\tau }_l.\end{array}$$

Then, the augmented system subject to the constrained input control $S_a$ can be described by the following state equation:$$\begin{array}{}\text{(39)}& S_a\begin{array}{c}\dot{X}t={\displaystyle \sum _{i=1}^{q_1}{\overline{A}}_i{X}^{i}t}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{D}}_{jl}{X}^{j}t-{\tau }_l}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{X}^{i}t\otimes X^{j}t-{\tau }_l}}}\hfill \\ +{\displaystyle \sum _{i=1}^{q_2}{\overline{B}}_{i}satut\otimes X^{i}t}+{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{C}}_{jl}satut\otimes X^{j}t-{\tau }_l}}\hfill \\ +{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{Q}}_{ijl}satut\otimes X^{i}t\otimes X^{j}t-{\tau }_l}}}+\overline{B}satut+\overline{E}{y}_{c}t\hfill \\ Yt=\overline{C}Xt\hfill \\ Xt=\overline{\zeta }t={\begin{array}{cc}\zeta t^T& 0_{1,wp}\end{array}}^T\hspace{1em}\text{for}\hspace{1em}t\in -\overline{\tau },0,\hfill \end{array}\end{array}$$with ${\overline{A}}_1=\begin{array}{cccccc}{A}_1& 0_{n,p}& 0_{n,p}& \dots & \dots & 0_{n,p}\\ 0_{p,n}& 0_{p,p}& I_p& 0_{p,p}& \dots & 0_{p,p}\\ 0_{p,n}& 0_{p,p}& 0_{p,p}& I_p& \ddots & ⋮\\ ⋮& ⋮& ⋮& \ddots & \ddots & 0_{p,p}\\ 0_{p,n}& 0_{p,p}& 0_{p,p}& \dots & 0_{p,p}& I_p\\ -C& 0_{p,p}& 0_{p,p}& \dots & \dots & 0_{p,p}\end{array},\hspace{1em}\overline{E}=\begin{array}{c}{0}_{n+w-1p,p}\\ I_p\end{array},\hspace{1em}\overline{C}=\begin{array}{cc}C& 0_{p,wp}\end{array},\hspace{1em}\overline{B}=\begin{array}{c}B\\ 0_{wp,m}\end{array}$, and for each $i\in 2\dots q_1$, $j\in 1\dots s_1$ and $l\in 1\dots v_1$$$\begin{array}{}\text{(39a)}& \begin{array}{c}{\overline{A}}_i=\begin{array}{c}{A}_i\\ 0_{wp,n^i}\end{array}{\Psi }^i,\hspace{1em}{\overline{D}}_{jl}=\begin{array}{c}{D}_{jl}\\ 0_{wp,n^j}\end{array}{\Psi }^j,\hspace{1em}{\overline{G}}_{ijl}=\begin{array}{c}{G}_{ijl}\\ 0_{wp,n^{i+j}}\end{array}{\Psi }^{i+j}\end{array},\end{array}$$and for each $i\in 1\dots q_2$, $j\in 1\dots s_2$ and $l\in 1\dots v_2$${\overline{B}}_i=\begin{array}{c}{B}_i\\ 0_{wp,m{n}^i}\end{array}{I}_m\otimes {\mathrm{\Psi }}^i,\hspace{1em}{\overline{C}}_{jl}=\begin{array}{c}{C}_{jl}\\ 0_{wp,m{n}^j}\end{array}{I}_m\otimes {\mathrm{\Psi }}^j,{\overline{Q}}_{ijl}=\begin{array}{c}{Q}_{ijl}\\ 0_{wp,m{n}^{i+j}}\end{array}{I}_m\otimes {\mathrm{\Psi }}^{i+j}$.

Assuming that the obtained augmented system $S_a$ is locally controllable around $X_0=\overline{\zeta }0$, then the dynamical behavior of the closed-loop augmented system $S_a$ must exactly match the behavior of a linear reference model $S_r$, which is represented also by its state equations:$$\begin{array}{}\text{(40)}& S_r\begin{array}{c}{\dot{X}}_{r}t=E{X}_{r}t+F{y}_{c}t\hfill \\ Y_{r}t=Z{X}_{r}t\hfill \\ X_{r}t=0_{n+wp,1}\hfill \end{array},\end{array}$$where $X_{r}t\in {ℝ}^{n+wp}$ and $Y_{r}t\in {ℝ}^p$.

The output vector $Y_{r}t$ of the above reference model, which should be adequately designed, must track perfectly the reference input vector $y_{c}t$.

Remark 1.

In order to make the adopted control strategy (26) in the adequate form, corresponding to the augmented system (39), it should be written as follows:$$\begin{array}{}\text{(41)}& ut=-{\displaystyle \sum _{i=1}^{q_1}{\overline{K}}_i{X}^{i}t}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{L}}_{jl}{X}^{j}t-{\tau }_l}},\end{array}$$with$$\begin{array}{}\text{(42)}& {\overline{K}}_1=K_1-\overline{N},\hspace{1em}{\overline{K}}_i=K_i{\Psi }^i,\hspace{1em}\forall i\in 2\dots q_1,\text{(43)}& {\overline{L}}_{jl}=L_{jl}{\Psi }^j\hspace{1em}\forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1.\end{array}$$

Remark 2.

The structure of the adopted strategy of control is perfectly adequate to the structure of the considered state model. This is justified by the following facts:

(i) The compensator gain $\overline{N}$ is designed to eliminate the steady-state error for a time-varying set point reference input vector:$$\begin{array}{}\text{(44)}& e\infty =y_c\infty -y\infty =0.\end{array}$$

(ii) The control gains $K_i$ and $L_{jl}$ for each $i\in 1,\dots ,q_1,j\in 1\dots s_1,l\in 1\dots v_1$ are designed to ensure the stability performance of the following closed-loop nonlinear polynomial time delays system under input control saturation $\overline{u}t$:$$\begin{array}{c}\text{(45)}& \dot{x}t={\displaystyle \sum _{i=1}^{q_1}{A}_i{x}^{i}t}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{D}_{jl}{x}^{j}t-{\tau }_l}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{G}_{ijl}{x}^{i}t\otimes x^{j}t-{\tau }_l}}}+{\displaystyle \sum _{i=1}^{q_2}{B}_{i}sat\overline{u}t\otimes x^{i}t}\\ +{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{C}_{jl}sat\overline{u}t\otimes x^{j}t-{\tau }_l}}\\ +{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{Q}_{ijl}sat\overline{u}t\otimes x^{i}t\otimes x^{j}t-{\tau }_l}}}+Bsat\overline{u}t,\end{array}$$

  with$$\begin{array}{}\text{(46)}& \begin{array}{c}\overline{u}t=-{\displaystyle \sum _{i=1}^{q_1}{K}_i{x}^{i}t}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{L}_{jl}{x}^{j}t-{\tau }_l}}\end{array}.\end{array}$$

4. Major Outcomes

4.1. Proposed Tracking Control Approach for Nonlinear Polynomial Time Delays System under Input Saturation

Let us define the following matrices which depend only on control parameters:$$\begin{array}{c}\text{(47)}& \forall i\in 1\dots q_1,\hspace{1em}\forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,M_i={\overline{A}}_i-\overline{B}{\overline{K}}_i,\hspace{1em}{H}_{jl}={\overline{D}}_{jl}-\overline{B}{\overline{L}}_{jl}.\\ \text{(48)}& \forall i\in 1\dots q_1,\hspace{1em}\forall j\in 1\dots q_2,{\overline{\mathcal{L}}}_{ij}={\overline{B}}_j{\overline{K}}_i\otimes I_{{n+wp}^j},\\ \text{(49)}& \forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,\hspace{1em}\forall i\in 1\dots s_2,\hspace{1em}\forall h\in 1\dots v_2,{\overline{\mathcal{V}}}_{jlih}={\overline{C}}_{ih}{\overline{L}}_{jl}\otimes I_{{n+wp}^i},\\ \text{(50)}& \forall i\in 1\dots q_1,\hspace{1em}\forall h\in 1\dots q_2,\hspace{1em}\forall j\in 1\dots s_2,\hspace{1em}\forall l\in 1\dots v_2,{\overline{\mathcal{U}}}_{ijl}={\overline{C}}_{jl}{\overline{K}}_i\otimes I_{{n+wp}^j},\hspace{1em}{\overline{\mathcal{Y}}}_{ihjl}={\overline{Q}}_{hjl}{\overline{K}}_i\otimes I_{{n+wp}^{h+j}},\\ \text{(51)}& \forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,\hspace{1em}\forall i\in 1\dots q_2,\hspace{1em}\forall h\in 1\dots s_2,\hspace{1em}\forall c\in 1\dots v_2{\overline{\mathcal{N}}}_{jli}={\overline{B}}_i{\overline{L}}_{jl}\otimes I_{{n+wp}^i},\hspace{1em}{\overline{\mathcal{Z}}}_{jlihc}={\overline{Q}}_{ihc}{\overline{L}}_{jl}\otimes I_{{n+wp}^{i+h}}.\end{array}$$

By substituting the constrained inputs control $satut$ by the unconstrained inputs control $ut$ in state equation (39) and by taking into account the relations (41) and (47)–(51), it comes out:$$\begin{array}{c}\text{(52)}& \dot{X}t={\displaystyle \sum _{i=1}^{q_1}{M}_i{X}^{i}t}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{H}_{jl}{X}^{j}t-{\tau }_l}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{X}^{i}t\otimes X^{j}t-{\tau }_l}}}+\overline{E}{y}_{c}t-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{X}^{i}t\otimes X^{j}t}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{X}^{j}t-{\tau }_l\otimes X^{i}t}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{X}^{i}t\otimes X^{j}t-{\tau }_l}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{X}^{j}t-{\tau }_l\otimes X^{i}t-{\tau }_h}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{X}^{i}t\otimes X^{h}t\otimes X^{j}t-{\tau }_l}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{X}^{j}t-{\tau }_l\otimes X^{i}t\otimes X^{h}t-{\tau }_c.}}}}}\end{array}$$

With respect to $t$ over the interval $0,T$, the integration of equation (52) leads to$$\begin{array}{c}\text{(53)}& Xt-X0={\displaystyle \sum _{i=1}^{q_1}{M}_i{\displaystyle \underset{0}{\overset{t}{\int }}{X}^i\sigma d\sigma }}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{H}_{jl}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^j\sigma -{\tau }_{l}d\sigma }}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^i\sigma \otimes X^j\sigma -{\tau }_{l}d\sigma }}}}+\overline{E}{\displaystyle \underset{0}{\overset{t}{\int }}{y}_c\sigma d\sigma }\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^i\sigma \otimes X^j\sigma d\sigma }}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^j\sigma -{\tau }_l\otimes X^i\sigma d\sigma }}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^i\sigma \otimes X^j\sigma -{\tau }_{l}d\sigma }}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^j\sigma -{\tau }_l\otimes X^i\sigma -{\tau }_{h}d\sigma }}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^i\sigma \otimes X^h\sigma \otimes X^j\sigma -{\tau }_{l}d\sigma }}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{\displaystyle \underset{0}{\overset{t}{\int }}{X}^j\sigma -{\tau }_l\otimes X^i\sigma \otimes X^h\sigma -{\tau }_{c}d\sigma }}}}}}.\end{array}$$

Over the basis of the block-pulse functions truncated to an order $N$, the expansion of the fixed reference input vector $y_{c}t$ can be done as follows:$$\begin{array}{}\text{(54)}& \begin{array}{c}{y}_{c}t=Y_{1,cN}{S}_{N}t\end{array},\end{array}$$where $Y_{1,cN}$ denote reference input coefficients, computing from the scalar product (15).

We indicate that the essential key is based on the equalization of the undelayed state vector of the controlled augmented system and the state vector of the reference model. By way of explanation $\forall t\in 0,T$:$$\begin{array}{}\text{(55)}& Xt=X_{r}t=X_{1,rN}{S}_{N}t,\end{array}$$where $X_{1,rN}$ denote state coefficients, computing from the scalar product (15).

As a result, the delayed state vector of the controlled augmented system $Xt-{\tau }_l$, for each $l\in 1\dots v$, can be expressed as follows:$$\begin{array}{}\text{(56)}& \begin{array}{c}Xt-{\tau }_l=\begin{array}{c}\overline{\zeta }t\hspace{1em}for\hspace{1em}0\le t\le {\tau }_l\hfill \\ Xt-{\tau }_l=X_{r}t-{\tau }_l\hspace{1em}for\hspace{1em}{\tau }_l<t\le T\hfill \end{array}\end{array},\end{array}$$and then, its expansion over the considered basis leads to$$\begin{array}{}\text{(57)}& \begin{array}{c}Xt-{\tau }_l=X_{1,rN}^{\ast }{\tau }_l{S}_{N}t\end{array},\end{array}$$where $X_{1,rN}^{\ast }{\tau }_l$ denote the delayed state coefficients, computing from the scalar product (20).

Using the operational matrix of product, the $i^{th}$ Kronecker power of $Xt$ could be expressed as follows:$$\begin{array}{c}\text{(58)}& X^{2}t=Xt\otimes Xt=X_{1,rN}{S}_{N}t\otimes X_{1,rN}{S}_{N}t=X_{1,rN}\otimes X_{1,rN}{S}_{N}t\otimes S_{N}t\\ =X_{1,rN}\otimes X_{1,rN}{M}_N{S}_{N}t=X_{2,rN}{S}_{N}t\\ ⋮\\ X^{q_1}t=X^{q_1-1}t\otimes Xt=X_{q_1-1,rN}{S}_{N}t\otimes X_{1,rN}{S}_{N}t\\ =X_{q_1-1,rN}\otimes X_{1,rN}{S}_{N}t\otimes S_{N}t\\ =X_{q_1-1,rN}\otimes X_{1,rN}{M}_N{S}_{N}t=X_{q_1,rN}{S}_{N}t.\end{array}$$

In the same way, the $j^{th}$ Kronecker power of $Xt-{\tau }_l$, for each $l\in 1\dots v_1$, could be expressed as follows:$$\begin{array}{c}\text{(59)}& X^{2}t-{\tau }_l=Xt-{\tau }_l\otimes Xt-{\tau }_l=X_{1,rN}^{\ast }{\tau }_l{S}_{N}t\otimes X_{1,rN}^{\ast }{\tau }_l{S}_{N}t=X_{1,rN}^{\ast }{\tau }_l\otimes X_{1,rN}^{\ast }{\tau }_l{S}_{N}t\otimes S_{N}t=X_{1,rN}^{\ast }{\tau }_l\otimes X_{1,rN}^{\ast }{\tau }_l{M}_N{S}_{N}t\\ =X_{2,rN}^{\ast }{\tau }_l{S}_{N}t\\ ⋮\\ X^{s_1}t-{\tau }_l=X^{s_1-1}t-{\tau }_l\otimes Xt-{\tau }_l\\ =X_{s_1-1,rN}^{\ast }{\tau }_l{S}_{N}t\otimes X_{1,rN}^{\ast }{\tau }_l{S}_{N}t\\ =X_{s_1-1,rN}^{\ast }{\tau }_l\otimes X_{1,rN}^{\ast }{\tau }_l{S}_{N}t\otimes S_{N}t\\ =X_{s_1-1,rN}^{\ast }{\tau }_l\otimes X_{1,rN}^{\ast }{\tau }_l{M}_N{S}_{N}t=X_{s_1,rN}^{\ast }{\tau }_l{S}_{N}t.\end{array}$$

In addition to this, for each $i\in 1\dots q_1$, $j\in 1\dots s_1$, and $l\in 1\dots v_1$, the Kronecker product terms $X^{i}t\otimes X^{j}t-{\tau }_l$ can be expressed as follows:$$\begin{array}{c}\text{(60)}& X^{i}t\otimes X^{j}t-{\tau }_l=X_{i,rN}{S}_{N}t\otimes X_{j,rN}^{\ast }{\tau }_l{S}_{N}t=X_{i,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\otimes S_{N}t=X_{i,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{M}_N{S}_{N}t.\end{array}$$

As well, for each $i\in 1\dots q_1$ and $j\in 1\dots q_2$, the Kronecker product terms $X^{i}t\otimes X^{j}t$ can be also expressed as follows:$$\begin{array}{c}\text{(61)}& X^{i}t\otimes X^{j}t=X_{i,rN}{S}_{N}t\otimes X_{j,rN}{S}_{N}t=X_{i,rN}\otimes X_{j,rN}{S}_{N}t\otimes S_{N}t=X_{i,rN}\otimes X_{j,rN}{M}_N{S}_{N}t.\end{array}$$

As well, for each $i\in 1\dots q_2$, $j\in 1\dots s_1$, and $l\in 1\dots v_1$, the Kronecker product terms $X^{j}t-{\tau }_l\otimes X^{i}t$ can be also expressed as follows:$$\begin{array}{c}\text{(62)}& X^{j}t-{\tau }_l\otimes X^{i}t=X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\otimes X_{i,rN}{S}_{N}t=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}{S}_{N}t\otimes S_{N}t=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}{M}_N{S}_{N}t.\hfill \end{array}$$

As well, for each $i\in 1\dots q_1$, $j\in 1\dots s_2$, and $l\in 1\dots v_2$, the Kronecker product terms $X^{i}t\otimes X^{j}t-{\tau }_l$ can be also expressed as follows:$$\begin{array}{c}\text{(63)}& X^{i}t\otimes X^{j}t-{\tau }_l=X_{i,rN}{S}_{N}t\otimes X_{j,rN}^{\ast }{\tau }_l{S}_{N}t=X_{i,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\otimes S_{N}t=X_{i,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{M}_N{S}_{N}t.\end{array}$$

As well, for each $i\in 1\dots s_2$, $j\in 1\dots s_1$, $l\in 1\dots v_1$, and $h\in 1\dots v_2$, the Kronecker product terms $X^{j}t-{\tau }_l\otimes X^{i}t-{\tau }_h$ can be also expressed as follows:$$\begin{array}{c}\text{(64)}& X^{j}t-{\tau }_l\otimes X^{i}t-{\tau }_h=X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\otimes X_{i,rN}^{\ast }{\tau }_h{S}_{N}t=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}^{\ast }{\tau }_h{S}_{N}t\otimes S_{N}t=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}^{\ast }{\tau }_h{M}_N{S}_{N}t.\hfill \end{array}$$

As well, for each $i\in 1\dots q_1$, $h\in 1\dots q_2$, $j\in 1\dots s_2$, and $l\in 1\dots v_2$, the Kronecker product terms $X^{i}t\otimes X^{h}t\otimes X^{j}t-{\tau }_l$ can be also expressed as follows:$$\begin{array}{}\text{(65)}& \begin{array}{c}{X}^{i}t\otimes X^{h}t\otimes X^{j}t-{\tau }_l\\ =X_{i,rN}{S}_{N}t\otimes X_{h,rN}{S}_{N}t\otimes X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\\ =X_{i,rN}\otimes X_{h,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\otimes S_{N}t\otimes S_{N}t\\ =X_{i,rN}\otimes X_{h,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{I}_N\otimes M_N{M}_N{S}_{N}t.\end{array}\end{array}$$

As well, for each $i\in 1\dots q_2$, $h\in 1\dots s_2$, $j\in 1\dots s_1$, $l\in 1\dots v_1$, and $c\in 1\dots v_2$, the Kronecker product terms $X^{j}t-{\tau }_l\otimes X^{i}t\otimes X^{h}t-{\tau }_c$ can be also expressed as follows:$$\begin{array}{c}\text{(66)}& X^{j}t-{\tau }_l\otimes X^{i}t\otimes X^{h}t-{\tau }_c=X_{j,rN}^{\ast }{\tau }_l{S}_{N}t\otimes X_{i,rN}{S}_{N}t\otimes X_{h,rN}^{\ast }{\tau }_c{S}_{N}t\\ =X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}\otimes X_{h,rN}^{\ast }{\tau }_c{S}_{N}t\otimes S_{N}t\otimes S_{N}t\\ =X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}\otimes X_{h,rN}^{\ast }{\tau }_c{I}_N\otimes M_N{M}_N{S}_{N}t.\end{array}$$

For simplicity, let us denote from relations (60)–(66) the following constants matrices:$$\begin{array}{c}\text{(67)}& \begin{array}{c}\forall i\in 1\dots q_1,\hspace{1em}\forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,\hspace{1em}{\mathcal{G}}_{ijl}=X_{i,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{M}_N\end{array},\text{(68)}& \begin{array}{c}\forall i\in 1\dots q_1,\hspace{1em}\forall j\in 1\dots q_2,\hspace{1em}{\mathrm{ℒ}}_{ij}=X_{i,rN}\otimes X_{j,rN}{M}_N\end{array},\text{(69)}& \forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,\hspace{1em}\forall i\in 1\dots s_2,\hspace{1em}\forall h\in 1\dots v_2,{\mathcal{V}}_{jlih}=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}^{\ast }{\tau }_h{M}_N,\\ \text{(70)}& \forall i\in 1\dots q_1,\hspace{1em}\forall h\in 1\dots q_2,\hspace{1em}\forall j\in 1\dots s_2,\hspace{1em}\forall l\in 1\dots v_2,{\mathcal{U}}_{ijl}=X_{i,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{M}_N,\hspace{1em}{\mathcal{Y}}_{ihjl}=X_{i,rN}\otimes X_{h,rN}\otimes X_{j,rN}^{\ast }{\tau }_l{I}_N\otimes M_N{M}_N,\\ \text{(71)}& \forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,\hspace{1em}\forall i\in 1\dots q_2,\hspace{1em}\forall h\in 1\dots s_2,\hspace{1em}\forall c\in 1\dots v_2{\mathcal{N}}_{jli}=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}{M}_N,\\ {\mathcal{Z}}_{jlihc}=X_{j,rN}^{\ast }{\tau }_l\otimes X_{i,rN}\otimes X_{h,rN}^{\ast }{\tau }_c{I}_N\otimes M_N{M}_N.\end{array}$$

Then, it results from the expansion of the whole state model (53) over the considered orthogonal functions basis:$$\begin{array}{c}\text{(72)}& X_{1,rN}{S}_{N}t-X_{0N}{S}_{N}t={\displaystyle \sum _{i=1}^{q_1}{M}_i{X}_{i,rN}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{H}_{jl}{X}_{j,rN}^{\ast }{\tau }_l{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{\mathcal{G}}_{ijl}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}}+\overline{E}{Y}_{1,cN}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{\mathrm{ℒ}}_{ij}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{\mathcal{N}}_{jli}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}}-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{u}}}_{ijl}{\mathcal{U}}_{ijl}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{\mathcal{V}}_{jlih}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}}}-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{\mathcal{Y}}_{ihjl}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{\mathcal{Z}}_{jlihc}{\displaystyle {\int }_0^t{S}_N\sigma d\sigma }}}}}}.\end{array}$$where $X_{0N}=\begin{array}{cccc}{x}_0& 0_{n,1}& \dots & 0_{n,1}\\ 0_{wp,1}& 0_{wp,1}& \dots & 0_{wp,1}\end{array}\in R^{n+wp\times N}.$

The application of the integration operational matrix $P_N$, defined by relation (22), leads to$$\begin{array}{c}\text{(73)}& X_{1,rN}{S}_{N}t-X_{0N}{S}_{N}t={\displaystyle \sum _{i=1}^{q_1}{M}_i{X}_{i,rN}{P}_N{S}_{N}t}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{H}_{jl}{X}_{j,rN}^{\ast }{\tau }_l{P}_N{S}_{N}t}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{\mathcal{G}}_{ijl}{P}_N{S}_{N}t}}}+\overline{E}{Y}_{1,cN}{P}_N{S}_{N}t\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{\mathcal{L}}_{ij}{P}_N{S}_{N}t}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{\mathcal{N}}_{jli}{P}_N{S}_{N}t}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{\mathcal{U}}_{ijl}{P}_N{S}_{N}t}}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{\mathcal{V}}_{jlih}{P}_N{S}_{N}t}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{\mathcal{Y}}_{ihjl}{P}_N{S}_{N}t}}}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{\mathcal{Z}}_{jlihc}{P}_N{S}_{N}t.}}}}}\end{array}$$

By removing the block-pulse vector basis $S_{N}t$ in both sides of the last relation (73), we attain the following algebraic equation:$$\begin{array}{c}\text{(74)}& X_{1,rN}-X_{0N}={\displaystyle \sum _{i=1}^{q_1}{M}_i{X}_{i,rN}{P}_N}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{H}_{jl}{X}_{j,rN}^{\ast }{\tau }_l{P}_N}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{\mathcal{G}}_{ijl}{P}_N}}}+\overline{E}{Y}_{1,cN}{P}_N-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{\mathcal{ℒ}}_{ij}{P}_N}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{\mathcal{N}}_{jli}{P}_N}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{\mathcal{U}}_{ijl}{P}_N}}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{\mathcal{V}}_{jlih}{P}_N}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{\mathcal{Y}}_{ihjl}{P}_N}}}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{\mathcal{Z}}_{jlihc}{P}_N.}}}}}\end{array}$$

By taking into account the relations (47)–(51), it results the following algebraic relation, where the only unknowns are the control law parameters:$$\begin{array}{c}\text{(75)}& X_{1,rN}-X_{0N}={\displaystyle \sum _{i=1}^{q_1}{\overline{A}}_i-\overline{B}{\overline{K}}_i{X}_{i,rN}{P}_N}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{D}}_{jl}-\overline{B}{\overline{L}}_{jl}{X}_{j,rN}^{\ast }{\tau }_l{P}_N}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{\mathcal{G}}_{ijl}{P}_N}}}+\overline{E}{Y}_{1,cN}{P}_N-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{B}}_j{\overline{K}}_i\otimes I_{{n+wp}^j}{\mathrm{ℒ}}_{ij}{P}_N}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{B}}_i{\overline{L}}_{jl}\otimes I_{{n+wp}^i}{\mathcal{N}}_{jli}{P}_N}}}-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{C}}_{jl}{\overline{K}}_i\otimes I_{{n+wp}^j}{\mathcal{U}}_{ijl}{P}_N}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{C}}_{ih}{\overline{L}}_{jl}\otimes I_{{n+wp}^i}{\mathcal{V}}_{jlih}{P}_N}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{Q}}_{hjl}{\overline{K}}_i\otimes I_{{n+wp}^{h+j}}{\mathcal{Y}}_{ihjl}{P}_N}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{Q}}_{ihc}{\overline{L}}_{jl}\otimes I_{{n+wp}^{i+h}}{\mathcal{Z}}_{jlihc}{P}_N.}}}}}\end{array}$$

Using the vectorization operator $vec$ (see Appendix B in [34]) allows us to develop the last relation (75) in the following form:$$\begin{array}{c}\text{(76)}& \beta =-{\displaystyle \sum _{i=1}^{q_1}{\chi }_{i}vec{\overline{K}}_i}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\eta }_{jl}vec{\overline{L}}_{jl}}}-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\xi }_{ij}vec{\overline{K}}_i\otimes I_{{n+wp}^j}}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\sigma }_{jli}vec{\overline{L}}_{jl}\otimes I_{{n+wp}^i}}}}-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\zeta }_{ijl}vec{\overline{K}}_i\otimes I_{{n+wp}^j}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\mu }_{jlih}vec{\overline{L}}_{jl}\otimes I_{{n+wp}^i}}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\varpi }_{ihjl}vec{\overline{K}}_i\otimes I_{{n+wp}^{h+j}}}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\kappa }_{jlihc}vec{\overline{L}}_{jl}\otimes I_{{n+wp}^{i+h}}}}}}},\end{array}$$with $\beta =vec{X}_{1,rN}-vec{X}_{0N}-{\displaystyle {\sum }_{i=1}^{q_1}{P}_N^T\otimes {\overline{A}}_{i}vec{X}_{i,rN}-P_N^T\otimes \overline{E}vec{Y}_{1,cN}}-{\displaystyle {\sum }_{j=1}^{s_1}{\displaystyle {\sum }_{l=1}^{v_1}{P}_N^T\otimes {\overline{D}}_{jl}vec{X}_{j,rN}^{\ast }{\tau }_l}}-{\displaystyle {\sum }_{i=1}^{q_1}{\displaystyle {\sum }_{j=1}^{s_1}{\displaystyle {\sum }_{l=1}^{v_1}{P}_N^T\otimes {\overline{G}}_{ijl}vec{\mathcal{G}}_{ijl}}}},$ where for each $i\in 1\dots q_1$, $j\in 1\dots s_1$, $l\in 1\dots v_1$: $\begin{array}{c}{\chi }_i={X_{i,rN}{P}_N}^T\otimes \overline{B},\hspace{1em}{\eta }_{jl}={X_{j,rN}^{\ast }{\tau }_l{P}_N}^T\otimes \overline{B}\end{array}$ and for each $i\in 1\dots q_1$ and $j\in 1\dots q_2$$\begin{array}{c}{\xi }_{ij}={{\mathrm{ℒ}}_{ij}{P}_N}^T\otimes {\overline{B}}_j\end{array}$ and for each $j\in 1\dots s_1$, $l\in 1\dots v_1$, $i\in 1\dots s_2$, and $h\in 1\dots v_2$$\begin{array}{c}{\mu }_{jlih}={{\mathcal{V}}_{jlih}{P}_N}^T\otimes {\overline{C}}_{ih}\end{array}$ and for each $i\in 1\dots q_1$, $h\in 1\dots q_2$, $j\in 1\dots s_2$, $l\in 1\dots v_2$: $\begin{array}{c}{\zeta }_{ijl}={{\mathcal{U}}_{ijl}{P}_N}^T\otimes {\overline{C}}_{jl},\hspace{1em}{\varpi }_{ihjl}={{\mathcal{Y}}_{ihjl}{P}_N}^T\otimes {\overline{Q}}_{hjl}\end{array}$ and for each $j\in 1\dots s_1$, $l\in 1\dots v_1$, $i\in 1\dots q_2$, $h\in 1\dots s_2$, $c\in 1\dots v_2$: $\begin{array}{c}{\sigma }_{jli}={{\mathcal{N}}_{jli}{P}_N}^T\otimes {\overline{B}}_i,{\kappa }_{jlihc}={{\mathcal{Z}}_{jlihc}{P}_N}^T\otimes {\overline{Q}}_{ihc}.\end{array}$

With the famous property of the vectorization operator, given for each matrix $A\in {ℝ}^{m\times n}$ by$$\begin{array}{}\text{(77)}& \begin{array}{c}vecA\otimes I_p=\begin{array}{c}vec{E}_{11}^{m\times n}\otimes I_p⋮\dots ⋮vec{E}_{m1}^{m\times n}\otimes I_p⋮\\ vec{E}_{22}^{m\times n}\otimes I_p⋮\dots ⋮vec{E}_{m2}^{m\times n}\otimes I_p⋮\dots ⋮\\ vec{E}_{1n}^{m\times n}\otimes I_p⋮\dots ⋮vec{E}_{mn}^{m\times n}\otimes I_p\end{array}={\Pi }_{m,n}{I}_{p}vecA\end{array},\end{array}$$where $E_{ij}^{m\times n}$ is the elementary matrix defined in relation (23), it results the following outcome from the previous relation (76):$$\begin{array}{c}\text{(78)}& \beta =-{\displaystyle \sum _{i=1}^{q_1}{\chi }_{i}vec{\overline{K}}_i}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\eta }_{jl}vec{\overline{L}}_{jl}}}-{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\xi }_{ij}{\Pi }_{m,{n+wp}^i}{I}_{{n+wp}^j}vec{\overline{K}}_i}}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\sigma }_{jli}{\Pi }_{m,{n+wp}^j}{I}_{{n+wp}^i}vec{\overline{L}}_{jl}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\zeta }_{ijl}{\Pi }_{m,{n+wp}^i}{I}_{{n+wp}^j}vec{\overline{K}}_i}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\mu }_{jlih}{\Pi }_{m,{n+wp}^j}{I}_{{n+wp}^i}vec{\overline{L}}_{jl}}}}}\\ -{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\varpi }_{ihjl}{\Pi }_{m,{n+wp}^i}{I}_{{n+wp}^{h+j}}vec{\overline{K}}_i}}}}\\ -{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\kappa }_{jlihc}{\Pi }_{m,{n+wp}^j}{I}_{{n+wp}^{i+h}}vec{\overline{L}}_{jl}.}}}}}\end{array}$$

It should be borne in mind that our goal is the gain determination of the adopted control law. Thus, from relations (42) and (43), we can infer that$$\begin{array}{}\text{(79)}& \begin{array}{c}vec{\overline{K}}_1=\begin{array}{c}{I}_{mn}\\ 0_{mwp,mn}\end{array}vec{K}_1-\begin{array}{c}{0}_{mn,mwp}\\ I_{mwp}\end{array}vec\overline{N}\end{array},\end{array}$$and for each $i\in 2\dots q_1$, $j\in 1\dots s_1$, $l\in 1\dots v_1$:$$\begin{array}{}\text{(80)}& \begin{array}{c}vec{\overline{K}}_i={{\Psi }^i}^T\otimes I_{m}vec{K}_i,\hspace{1em}vec{\overline{L}}_{jl}={{\Psi }^j}^T\otimes I_{m}vec{L}_{jl}\end{array}.\end{array}$$

Consequently, relation (78) can be transformed into the following form:$$\begin{array}{}\text{(81)}& \begin{array}{c}\beta ={\alpha }_{1}vec\overline{N}-{\displaystyle \sum _{i=1}^{q_1}{\overline{\chi }}_{i}vec{K}_i}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{\eta }}_{jl}vec{L}_{jl}}}\end{array},\end{array}$$with $\begin{array}{c}{\alpha }_1={\alpha }_1^{\ast }\begin{array}{c}{0}_{mn,mwp}\\ I_{mwp}\end{array},\hspace{1em}{\overline{\chi }}_1={\overline{\chi }}_1^{\ast }\begin{array}{c}{I}_{mn}\\ 0_{mwp,mn}\end{array}\end{array}$ where$$\begin{array}{c}\text{(82)}& {\alpha }_1^{\ast }={\chi }_1+{\displaystyle \sum _{j=1}^{q_2}{\xi }_{1j}{\Pi }_{m,n+wp}{I}_{{n+wp}^j}}+{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\zeta }_{1jl}{\Pi }_{m,n+wp}{I}_{{n+wp}^j}}}+{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\varpi }_{1hjl}{\Pi }_{m,n+wp}{I}_{{n+wp}^{h+j}}}}},\\ {\overline{\chi }}_1^{\ast }={\chi }_1+{\displaystyle \sum _{j=1}^{q_2}{\xi }_{1j}{\Pi }_{m,n+wp}{I}_{{n+wp}^j}}+{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\zeta }_{1jl}{\Pi }_{m,n+wp}{I}_{{n+wp}^j}}}\\ +{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\varpi }_{1hjl}{\Pi }_{m,n+wp}{I}_{{n+wp}^{h+j}}}}},\end{array}$$and for each $i\in 2\dots q_1$, $j\in 1\dots s_1$, $l\in 1\dots v_1$: $\begin{array}{c}{\overline{\chi }}_i={\overline{\chi }}_i^{\ast }{{\mathrm{\Psi }}^i}^T\otimes I_m,\hspace{1em}{\overline{\eta }}_{jl}={\overline{\eta }}_{jl}^{\ast }{{\mathrm{\Psi }}^j}^T\otimes I_m\end{array}$ where ${\overline{\chi }}_i^{\ast }={\chi }_i+{\displaystyle {\sum }_{j=1}^{q_2}{\xi }_{ij}{\mathrm{\Pi }}_{m,{n+wp}^i}{I}_{{n+wp}^j}}+{\displaystyle {\sum }_{j=1}^{s_2}{\displaystyle {\sum }_{l=1}^{v_2}{\zeta }_{ijl}{\mathrm{\Pi }}_{m,{n+wp}^i}{I}_{{n+wp}^j}}}+{\displaystyle {\sum }_{h=1}^{q_2}{\displaystyle {\sum }_{j=1}^{s_2}{\displaystyle {\sum }_{l=1}^{v_2}{\varpi }_{ihjl}{\mathrm{\Pi }}_{m,{n+wp}^i}{I}_{{n+wp}^{h+j}}}}}$${\overline{\eta }}_{jl}^{\ast }={\eta }_{jl}+{\displaystyle {\sum }_{i=1}^{q_2}{\sigma }_{jli}{\mathrm{\Pi }}_{m,{n+wp}^j}{I}_{{n+wp}^i}}+{\displaystyle {\sum }_{i=1}^{s_2}{\displaystyle {\sum }_{h=1}^{v_2}{\mu }_{jlih}{\mathrm{\Pi }}_{m,{n+wp}^j}{I}_{{n+wp}^i}}}+{\displaystyle {\sum }_{i=1}^{q_2}{\displaystyle {\sum }_{h=1}^{s_2}{\displaystyle {\sum }_{c=1}^{v_2}{\kappa }_{jlihc}{\mathrm{\Pi }}_{m,{n+wp}^j}{I}_{{n+wp}^{i+h}}}}}.$

Letting $\begin{array}{c}{\alpha }_2=-\begin{array}{cccc}{\overline{\chi }}_1& {\overline{\chi }}_2& \dots & {\overline{\chi }}_{q_1}\end{array}\end{array}$ and $\begin{array}{c}\forall \hspace{1em}j\in 1\dots s_1,\hspace{1em}{\gamma }_j=-\begin{array}{cccc}{\overline{\eta }}_{j1}& {\overline{\eta }}_{j2}& \dots & {\overline{\eta }}_{j{v}_1}\end{array}\end{array}$then, from equation (81), it comes out$$\begin{array}{}\text{(83)}& \begin{array}{c}\beta ={\alpha }_{1}vec\overline{N}+{\alpha }_2{\theta }_K+{\displaystyle \sum _{j=1}^{s_1}{\gamma }_j{\theta }_{Lj}}\end{array},\end{array}$$with $\begin{array}{c}{\theta }_K=\begin{array}{c}vec{K}_1\\ ⋮\\ vec{K}_{q_1}\end{array}\end{array}$ and $\begin{array}{c}\forall \hspace{1em}j\in 1\dots s_1,\hspace{1em}{\theta }_{Lj}=\begin{array}{c}vec{L}_{j1}\\ ⋮\\ vec{L}_{j{v}_1}\end{array}\end{array}$ which is equivalent to$$\begin{array}{}\text{(84)}& \begin{array}{c}\beta ={\alpha }_{1}vec\overline{N}+{\alpha }_2{\theta }_K+{\alpha }_3{\theta }_L\end{array},\end{array}$$where $\begin{array}{c}{\alpha }_3=\begin{array}{cccc}{\gamma }_1& {\gamma }_2& \dots & {\gamma }_{s_1}\end{array}\mathrm{and}{\theta }_L=\begin{array}{c}{\theta }_{L1}\\ ⋮\\ {\theta }_{L{s}_1}\end{array}.\end{array}$

It is appreciable to systematize such problem as follows:$$\begin{array}{}\text{(85)}& \mathcal{A}\theta =\mathcal{B},\end{array}$$where constant matrix $\mathcal{A}$, constant vector $\mathcal{B}$, and vector of searched parameters $\theta$ are given by$$\begin{array}{c}\text{(86)}& \mathcal{A}=\begin{array}{ccc}{\alpha }_1& {\alpha }_2& {\alpha }_3\end{array},\mathcal{B}=\beta ,\\ \theta =\begin{array}{c}\mathrm{vec}\overline{N}\\ {\theta }_K\\ {\theta }_L\end{array}.\end{array}$$

Now, by taking into account that the control law $ut$ has to respect upper and lower bounds (30), and then from relation (41), we get the two following inequalities:$$\begin{array}{}\text{(87)}& \begin{array}{c}-{\displaystyle \sum _{i=1}^{q_1}{\overline{K}}_i{X}^{i}t}-{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{L}}_{jl}{X}^{j}t-{\tau }_l}}\le u_{\max }\\ {\displaystyle \sum _{i=1}^{q_1}{\overline{K}}_i{X}^{i}t}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{L}}_{jl}{X}^{j}t-{\tau }_l}}\le v_{\max }\end{array}.\end{array}$$

Undertaken the equality (55), the expansion of the latter inequalities over the considered block-pulse functions basis and the use of the $vec$ operator yield to the following inequality system:$$\begin{array}{}\text{(88)}& \begin{array}{c}\mathcal{C}\theta \le d\end{array},\end{array}$$where constant matrix $\mathcal{C}$ and constant vector $d$ are given as follows:$$\begin{array}{c}\text{(89)}& \mathcal{C}=\begin{array}{ccc}-{\mathcal{D}}_1\hfill & -{\mathcal{D}}_2\hfill & -{\mathcal{D}}_3\hfill \\ {\mathcal{D}}_1\hfill & {\mathcal{D}}_2\hfill & {\mathcal{D}}_3\hfill \end{array},d=\begin{array}{c}vec{U}_{\max ,N}\\ vec{V}_{\max ,N}\end{array}.\end{array}$$with$$\begin{array}{c}\text{(90)}& {\mathcal{D}}_1=-{X_{1,rN}}^T\otimes I_m\begin{array}{c}{0}_{mn,mwp}\\ I_{mwp}\end{array},{\mathcal{D}}_2=\begin{array}{cc}{X_{1,rN}}^T\otimes I_m\begin{array}{c}{I}_{mn}\\ 0_{mwp,mn}\end{array}& {\overline{\mathcal{D}}}_2\end{array},\\ {\mathcal{D}}_3=\begin{array}{ccc}{d}_1& \dots & d_{s_1}\end{array},\end{array}$$where$$\begin{array}{c}\text{(91)}& {\overline{\mathcal{D}}}_2=\begin{array}{ccc}{{\Psi }^2{X}_{2,rN}}^T\otimes I_m& \dots & {{\Psi }^{q_1}{X}_{q_1,rN}}^T\otimes I_m\end{array},d_j=\begin{array}{ccc}{{\Psi }^j{X}_{j,rN}^{\ast }{\tau }_1}^T\otimes I_m& \dots & {{\Psi }^j{X}_{j,rN}^{\ast }{\tau }_{v_1}}^T\otimes I_m\end{array},\end{array}$$and $U_{\max ,N}$ and $V_{\max ,N}$ denote the upper bound coefficients, which are computed from the scalar product (15).

In the end, the posed control problem has been effectively reduced to an overdetermined linear system subject to an inequality constraint, which can be formulated from equality (85) and inequality (88) as follows:$$\begin{array}{}\text{(92a)}& \mathcal{A}\theta =\mathcal{B},\end{array}$$subject to$$\begin{array}{}\text{(92b)}& \mathcal{C}\theta \le d.\end{array}$$

The above constrained linear system is of $n+wpN$ algebraic equations subject to $2mN$ algebraic constraints with respect to $mwp+{\displaystyle {\sum }_{i=1}^{q_1}{n}^i}+v_1{\displaystyle {\sum }_{j=1}^{s_1}{n}^j}$ unknowns.

This kind of system can be easily solved in the constrained least square sense, by the formulation as a quadratic programming problem as follows:$$\begin{array}{}\text{(93a)}& \underset{\theta }{\min }\mathcal{A}\theta -\mathcal{B},\end{array}$$subject to$$\begin{array}{}\text{(93b)}& \mathcal{C}\theta \le d.\end{array}$$

4.2. Practical Stability Problem Analysis

Presently, the control parameters $\overline{N}$, $K_i$, and $L_{jl}$ can be found by solving the resulted quadratic programming problem (93a)-(93b) for fixed reference input $y_{c}t$. Thereafter, we intend to analyze the practical stability of the closed-loop augmented system (39).

For this purpose, let us consider the following constants matrices:$$\begin{array}{c}\text{(94)}& \forall i\in 1\dots q_1,\hspace{1em}\forall j\in 1\dots s_1,\hspace{1em}\forall l\in 1\dots v_1,{\mathrm{ℳ}}_i={\overline{A}}_i-\frac{1}{2}\overline{B}{\overline{K}}_i,\hspace{1em}{\mathrm{\mathscr{H}}}_{jl}={\overline{D}}_{jl}-\frac{1}{2}\overline{B}{\overline{L}}_{jl}.\end{array}$$

Definition 1.

The controlled system (39) is called practically stable if every bounded input vector produces a bounded states vector. In other words, there exist $0<R_0<r$, such that (see [45]):$$\begin{array}{}\text{(95)}& X0<R_0\hspace{1em}\Rightarrow \hspace{1em}Xt<r,\hspace{1em}\forall t\ge 0.\end{array}$$

Theorem 1.

The controlled system (39) is practical stable if all eigenvalues of matrix ${\mathrm{ℳ}}_1$ have a strictly negative real part and if$$\begin{array}{}\text{(96)}& \begin{array}{c}{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}+{\omega }_1<0\end{array},\text{(97)}& \begin{array}{c}X0<R_0\end{array},\end{array}$$where $R_0>0$, ${\lambda }_1>0$ and ${\omega }_1<0$ are given scalars satisfying:$$\begin{array}{}\text{(98)}& e^{{\mathrm{ℳ}}_{1}t}\le {\lambda }_1{e}^{{\omega }_{1}t},\forall t\ge 0,\text{(99)}& R_0=\frac{R_1}{{\lambda }_1},\end{array}$$with $R_1$ is the unique positive solution of the following equation:$$\begin{array}{}\text{(100)}& f{R}_1+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}+{\omega }_1=0,\end{array}$$where$$\begin{array}{c}\text{(101)}& f{R}_1={\lambda }_1{\displaystyle \sum _{i=2}^{q_1}{\mathcal{M}}_i{R}_1^{i-1}}+{\displaystyle \sum _{j=2}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{jl}{R}_1^{j-1}}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{R}_1^{i+j-1}}}}+\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{R}_1^{i+j-1}}}+{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{R}_1^{j+i-1}}}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{R}_1^{i+j-1}}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{R}_1^{j+i-1}}}}}+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{R}_1^{i+h+j-1}}}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{R}_1^{j+i+h-1}}}}}}+{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{q_1}{\overline{B}}_i{\overline{K}}_j{R}_1^{i+j-1}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{B}}_i{\overline{L}}_{jl}{R}_1^{i+j-1}}}}+{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{i=1}^{q_1}{\overline{C}}_{jl}{\overline{K}}_i{R}_1^{i+j-1}}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{i=1}^{s_1}{\displaystyle \sum _{h=1}^{v_1}{\overline{C}}_{jl}{\overline{L}}_{ih}{R}_1^{i+j-1}}}}}+{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{h=1}^{q_1}{\overline{Q}}_{ijl}{\overline{K}}_h{R}_1^{i+j+h-1}}}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{h=1}^{s_1}{\displaystyle \sum _{c=1}^{v_1}{\overline{Q}}_{ijl}{\overline{L}}_{hc}{R}_1^{i+j+h-1}}}}}}.\end{array}$$

Proof 1.

By taking into account that $satut=\overline{u}t+1/2ut$ with $\overline{u}t=satut-1/2ut$, the augmented closed-loop system can be written from relations (39), (41), (94), and (48)–(51) as follows:$$\begin{array}{c}\text{(102)}& \dot{X}t-{\mathcal{M}}_{1}Xt={\displaystyle \sum _{i=2}^{q_1}{\mathcal{M}}_i{X}^{i}t}+{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{1l}Xt-{\tau }_l}+{\displaystyle \sum _{j=2}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{jl}{X}^{j}t-{\tau }_l}}+\overline{B}\overline{u}t+{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{X}^{i}t{X}^{j}t-{\tau }_l}}}+\overline{E}{y}_{c}t-\frac{1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{X}^{i}t\otimes X^{j}t}}\\ -\frac{1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{X}^{j}t-{\tau }_l\otimes X^{i}t}}}-\frac{1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{X}^{i}t\otimes X^{j}t-{\tau }_l}}}\\ -\frac{1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{X}^{j}t-{\tau }_l\otimes X^{i}t-{\tau }_h}}}}\\ -\frac{1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{X}^{i}t\otimes X^{h}t\otimes X^{j}t-{\tau }_l}}}}\\ -\frac{1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{X}^{j}t-{\tau }_l\otimes X^{i}t\otimes X^{h}t-{\tau }_c}}}}}+{\displaystyle \sum _{i=1}^{q_2}{\overline{B}}_i\overline{u}t\otimes X^{i}t}\\ +{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{C}}_{jl}\overline{u}t\otimes X^{j}t-{\tau }_l}}+{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{Q}}_{ijl}\overline{u}t\otimes X^{i}t\otimes X^{j}t-{\tau }_l}}}.\end{array}$$

By respecting that ${\mathrm{ℳ}}_1{e}^{-{\mathrm{ℳ}}_{1}t}=e^{-{\mathrm{ℳ}}_{1}t}{\mathrm{ℳ}}_1$ which can be deduced from the series definition, we get the following:$$\begin{array}{}\text{(103)}& \begin{array}{c}\frac{d}{dt}{e}^{-{\mathrm{ℳ}}_{1}t}Xt=-{\mathrm{ℳ}}_1{e}^{-{\mathrm{ℳ}}_{1}t}Xt+e^{-{\mathrm{ℳ}}_{1}t}\dot{X}t=e^{-{\mathrm{ℳ}}_{1}t}\dot{X}t-{\mathrm{ℳ}}_{1}Xt\end{array}.\end{array}$$Thus, by integrating the last differential equation from 0 to $t$, it comes out:$$\begin{array}{}\text{(104)}& \begin{array}{c}{e}^{-{\mathrm{ℳ}}_{1}t}Xt-X0={\displaystyle {\int }_0^t{e}^{-{\mathrm{ℳ}}_1\sigma }\dot{X}\sigma -{\mathrm{ℳ}}_{1}X\sigma d\sigma }\end{array},\end{array}$$or equivalently$$\begin{array}{}\text{(105)}& Xt=e^{{\mathrm{ℳ}}_{1}t}X0+{\displaystyle {\int }_0^t{e}^{{\mathrm{ℳ}}_{1}t-\sigma }\dot{X}\sigma -{\mathrm{ℳ}}_{1}X\sigma d\sigma }.\end{array}$$

Since all eigenvalues of matrix ${\mathrm{ℳ}}_1$ must have the negative real part, then there exist two reals ${\lambda }_1>0$ and ${\omega }_1<0$ such that relation (98) holds on. Consequently, the state vector $Xt$ can be bounded as:$$\begin{array}{}\text{(106)}& Xt\le {\lambda }_1{e}^{{\omega }_{1}t}X0+{\lambda }_1{e}^{{\omega }_{1}t}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }\dot{X}\sigma -{\mathrm{ℳ}}_{1}X\sigma d\sigma },\end{array}$$which is similar to the following:$$\begin{array}{}\text{(107)}& e^{-{\omega }_{1}t}Xt\le {\lambda }_{1}X0+{\lambda }_1{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }\dot{X}\sigma -{\mathrm{ℳ}}_{1}X\sigma d\sigma }.\end{array}$$

Let us suppose that$$\begin{array}{}\text{(108)}& Xt<R.\end{array}$$

Through the following matrix norm properties:$$\begin{array}{c}\text{(109)}& X^{i}t<R^{i-1}Xt,\hspace{1em}{X}^{j}t-{\tau }_l<R^{j-1}Xt-{\tau }_l{X}^{i}t\otimes X^{j}t-{\tau }_l<R^{i+j-1}Xt-{\tau }_l,\hspace{1em}{X}^{i}t\otimes X^{j}t<R^{i+j-1}Xt\\ X^{j}t-{\tau }_l\otimes X^{i}t<R^{j+i-1}Xt\\ X^i\sigma \otimes X^h\sigma \otimes X^j\sigma -{\tau }_l<R^{i+h+j-1}Xt-{\tau }_l\\ X^j\sigma -{\tau }_l\otimes X^i\sigma \otimes X^h\sigma -{\tau }_c<R^{j+i+h-1}Xt-{\tau }_c.\end{array}$$then we can deduce, from inequality (107) and relation (102), the following expression:$$\begin{array}{c}\text{(110)}& e^{-{\omega }_{1}t}Xt<{\lambda }_{1}X0+{\lambda }_1{\displaystyle \sum _{i=2}^{q_1}{\mathcal{M}}_i{R}^{i-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma d\sigma }}+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{1l}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{l}d\sigma }}+{\lambda }_1{\displaystyle \sum _{j=2}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{jl}{R}^{j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{l}d\sigma }}}\\ +{\lambda }_1{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{l}d\sigma }}}}+{\lambda }_1\overline{E}\delta {\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }d\sigma }\\ +{\lambda }_1\overline{B}\overline{\delta }{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }d\sigma }+\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma }}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{R}^{j+i-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma }}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_l}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{R}^{j+i-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_h}}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{R}^{i+h+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_l}}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{R}^{j+i+h-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_c}}}}}d\sigma }\\ +{\lambda }_1{\displaystyle \sum _{i=1}^{q_2}{\overline{B}}_i{R}^i{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }\overline{u}\sigma d\sigma }}+{\lambda }_1{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{C}}_{jl}{R}^j{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }\overline{u}\sigma d\sigma }}}\\ +{\lambda }_1{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{Q}}_{ijl}{R}^{i+j}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }\overline{u}\sigma d\sigma }}}}.\end{array}$$where $\delta >0$ satisfying $y_{c}t<\delta$, and $\overline{\delta }=\max \begin{array}{cc}{v}_{\max },& u_{\max }\end{array}$.

Now by taking into account the following special property of function saturation (see [46]):$$\begin{array}{}\text{(111)}& \overline{u}t=satut-\frac{1}{2}ut\le \frac{1}{2}ut,\end{array}$$which can be written from relation (41) and condition (108) as follows:$$\begin{array}{c}\text{(112)}& \overline{u}t<\frac{1}{2}{\displaystyle \sum _{i=1}^{q_1}{\overline{K}}_i{X}^{i}t}+\frac{1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{L}}_{jl}{X}^{j}t-{\tau }_l}}<\frac{1}{2}{\displaystyle \sum _{i=1}^{q_1}{\overline{K}}_i{R}^{i-1}Xt}+\frac{1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{L}}_{jl}{R}^{j-1}Xt-{\tau }_l}}.\end{array}$$

It results the following inequality:$$\begin{array}{c}\text{(113)}& e^{-{\omega }_{1}t}Xt<{\lambda }_{1}X0+{\lambda }_1{\displaystyle \sum _{i=2}^{q_1}{\mathcal{M}}_i{R}^{i-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma d\sigma }}+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{1l}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{h}d\sigma }}+{\lambda }_1{\displaystyle \sum _{j=2}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\mathcal{H}}_{jl}{R}^{j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{h}d\sigma }}}\\ +{\lambda }_1{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{G}}_{ijl}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{h}d\sigma }}}}+{\lambda }_1\overline{E}\delta {\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }d\sigma }\\ +{\lambda }_1\overline{B}\overline{\delta }{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }d\sigma }+\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{q_2}{\overline{\mathcal{L}}}_{ij}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma }}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\overline{\mathcal{N}}}_{jli}{R}^{j+i-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma }}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{U}}}_{ijl}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_l}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{s_2}{\displaystyle \sum _{h=1}^{v_2}{\overline{\mathcal{V}}}_{jlih}{R}^{j+i-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_h}}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_1}{\displaystyle \sum _{h=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\overline{\mathcal{Y}}}_{ihjl}{R}^{i+h+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_l}}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{h=1}^{s_2}{\displaystyle \sum _{c=1}^{v_2}{\overline{\mathcal{Z}}}_{jlihc}{R}^{j+i+h-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_c}}}}}d\sigma }\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{q_1}{\overline{B}}_i{\overline{K}}_j{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma d\sigma }}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_1}{\displaystyle \sum _{l=1}^{v_1}{\overline{B}}_i{\overline{L}}_{jl}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{l}d\sigma }}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{i=1}^{q_1}{\overline{C}}_{jl}{\overline{K}}_i{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma d\sigma }}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{i=1}^{s_1}{\displaystyle \sum _{h=1}^{v_1}{\overline{C}}_{jl}{\overline{L}}_{ih}{R}^{i+j-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{h}d\sigma }}}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{h=1}^{q_1}{\overline{Q}}_{ijl}{\overline{K}}_h{R}^{i+j+h-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma d\sigma }}}}}\\ +\frac{{\lambda }_1}{2}{\displaystyle \sum _{i=1}^{q_2}{\displaystyle \sum _{j=1}^{s_2}{\displaystyle \sum _{l=1}^{v_2}{\displaystyle \sum _{h=1}^{s_1}{\displaystyle \sum _{c=1}^{v_1}{\overline{Q}}_{ijl}{\overline{L}}_{hc}{R}^{i+j+h-1}{\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }X\sigma -{\tau }_{c}d\sigma }}}}}}.\end{array}$$

Let us introduce the following function:$$\begin{array}{}\text{(114)}& \begin{array}{c}zt={\sup }_{\rho \in t-\overline{\tau },t}X\rho e^{-{\omega }_{1}t}\hspace{1em}\forall t\ge 0\end{array}.\end{array}$$then under considering that $\forall t\ge 0$ and $l\in 1\dots v$: $\begin{array}{c}Xt\le {\sup }_{\rho \in t-\overline{\tau },t}X\rho \end{array}$ and $\begin{array}{c}Xt-{\tau }_l\le {\sup }_{\rho \in t-\overline{\tau },t}X\rho \end{array}$ it turns out, from inequality (113), that$$\begin{array}{}\text{(115)}& \begin{array}{c}{e}^{-{\omega }_{1}t}Xt<{\lambda }_{1}X0+fR+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}{\displaystyle {\int }_0^{t}z\sigma d\sigma }+\Omega {\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }d\sigma }\end{array},\end{array}$$where function $f.$ is given by relation (101) and $\mathrm{\Omega }={\lambda }_1\overline{E}\delta +{\lambda }_1\overline{B}\overline{\delta }$.

Now, by introducing a new function $gt$ as the right-hand member of inequality (115):$$\begin{array}{}\text{(116)}& \begin{array}{c}gt={\lambda }_{1}X0+fR+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}{\displaystyle {\int }_0^{t}z\sigma d\sigma }+\Omega {\displaystyle {\int }_0^t{e}^{-{\omega }_1\sigma }d\sigma }\end{array}.\end{array}$$we can show that$$\begin{array}{}\text{(117)}& \forall t\ge 0,zt<gt.\end{array}$$

By deriving the function $gt$:$$\begin{array}{}\text{(118)}& \begin{array}{c}\frac{dgt}{dt}=fR+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}zt+\Omega e^{-{\omega }_{1}t}\end{array},\end{array}$$it comes out from inequality (117) that$$\begin{array}{}\text{(119)}& \frac{dgt}{dt}<fR+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}gt+\Omega e^{-{\omega }_{1}t}.\end{array}$$

By integrating the last inequality from 0 to t, it follows, based on the use of the Bellman-Gronwall inequality (see appendix D in [34]) that$$\begin{array}{}\text{(120)}& \begin{array}{c}gt<e^{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}t}g0+\frac{\Omega }{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}+{\omega }_1}}{e}^{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}t}-e^{-{\omega }_{1}t}\end{array},\end{array}$$where $g0={\lambda }_{1}X0.$

Inequalities (117) and (120) imply$$\begin{array}{}\text{(121)}& zt<{\lambda }_{1}X0e^{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}t}+\frac{\Omega }{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}+{\omega }_1}}{e}^{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}t}-e^{-{\omega }_{1}t}.\end{array}$$

From the definition of function $zt$, we can derived that$$\begin{array}{}\text{(122)}& Xt<{\lambda }_{1}X0e^{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}+{\omega }_{1}t}+\frac{\Omega }{fR+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}+{\omega }_1}{e}^{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}+{\omega }_{1}t}-1.\end{array}$$

If condition (96) holds on, then for$$\begin{array}{}\text{(123)}& \begin{array}{c}fR+{\lambda }_1{\displaystyle \sum _{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}}+{\omega }_1<0\end{array}.\end{array}$$which is verified for each radius $R$ satisfying$$\begin{array}{}\text{(124)}& \begin{array}{c}R<R_1\end{array},\end{array}$$where $R_1$ is the unique positive solution of the equality (100), we can make out the following result:$$\begin{array}{}\text{(125)}& Xt<r={\lambda }_{1}X0-\frac{\Omega }{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}+{\omega }_1}}.\end{array}$$

At this point, in order to make sure the hypothesis (108) for all $t$, it is enough to have$$\begin{array}{}\text{(126)}& {\lambda }_{1}X0<R_1.\end{array}$$

Consequent, from inequality (125) and condition (126), it allows that$$\begin{array}{}\text{(127)}& X0<R_0=\frac{R_1}{{\lambda }_1}\hspace{1em}\Rightarrow \hspace{1em}Xt<r={\lambda }_1{R}_0-\frac{\Omega }{fR+{\lambda }_1{\displaystyle {\sum }_{l=1}^{v_1}{\mathrm{\mathscr{H}}}_{1l}+{\omega }_1}}.\end{array}$$

Hence, the controlled system (39) is practically stable.

5. Study of an Illustrative Benchmark System

In this section, simulations will be performed on the model of double-inverted pendulums, which can be described concisely in Figure 1.

[figure(s) omitted; refer to PDF]

The parameters $m_1=2\mathrm{kg}$, $J_1=5\mathrm{kg}\mathrm{.}{\mathrm{m}}^{\mathrm{2}}$ and $m_2=2.5\mathrm{kg}$, $J_2=6.5\mathrm{kg}\mathrm{.}{\mathrm{m}}^{\mathrm{2}}$ are, respectively, the masses and the moments of inertia of the two pendulums, $r=0.5\mathrm{m}$ is the pendulum height, $k=100\mathrm{N}/\mathrm{m}$ is the constant of the connecting torsional spring, $b=0.5\mathrm{m}$ is the distance between the pendulum hinges, $l=0.5\mathrm{m}$ is the spring natural length, and $g=9.81\mathrm{m}/{\mathrm{s}}^2$ is the gravitational acceleration.

5.1. System Modeling

According to [47], the considered benchmark system is described as follows:$$\begin{array}{}\text{(128)}& \begin{array}{c}\begin{array}{c}{\dot{x}}_{1}t=x_{2}t\hfill \\ {\dot{x}}_{2}t=\frac{m_{1}gr}{J_1}-\frac{k{r}^2}{4J_1}\sin x_{1}t+\frac{kr}{2J_1}l-b+\frac{k{r}^2}{4J_1}\sin x_{3}t\hfill \\ +f_1{x}_{1}t-{\tau }_1,x_{2}t+g_1{x}_{1}t,x_{3}tsat{u}_{1}t\hfill \\ {\dot{x}}_{3}t=x_{4}t\hfill \\ {\dot{x}}_{4}t=\frac{m_{2}gr}{J_2}-\frac{k{r}^2}{4J_2}\sin x_{3}t+\frac{kr}{2J_2}l-b+\frac{k{r}^2}{4J_2}\sin x_{1}t\hfill \\ +f_2{x}_{3}t-{\tau }_2,x_{4}t+g_2{x}_{1}t,x_{3}tsat{u}_{2}t\hfill \end{array}\end{array},\end{array}$$where $x_{1}t={\theta }_{1}t$ and $x_{3}t={\theta }_{2}t$ denote the angular displacements of the pendulums from vertical, $f_1.$ and $f_2.$ are nonlinear time delays functions associated, respectively, with each pendulum, and $g_1.$ and $g_2.$ are nonlinear functions associated, respectively, with each input channel $u_{1}t$ and $u_{2}t$. The output systems are ${\theta }_{1}t$ and ${\theta }_{2}t$.

In addition, we assumed that$$\begin{array}{c}\text{(129)}& f_1{x}_{1}t-{\tau }_1,x_{2}t=\frac{x_{1}t-{\tau }_1}{1+x_1^{2}t-{\tau }_1}+\frac{1}{J_1}{x}_{1}t-{\tau }_1{x}_{2}t\cos x_{1}t-{\tau }_1,f_2{x}_{3}t-{\tau }_2,x_{4}t=\frac{x_{3}t-{\tau }_2}{1+x_3^{2}t-{\tau }_2}+\frac{1}{J_2}{x}_{3}t-{\tau }_2{x}_{4}t\cos x_{3}t-{\tau }_2,\\ g_1{x}_{1}t,x_{3}t=\frac{1}{J_1}-0.25x_1^{2}t+0.25x_3^{2}t,\\ g_2{x}_{1}t,x_{3}t=\frac{1}{J_2}-0.25x_1^{2}t+0.25x_3^{2}t.\end{array}$$

For simulation, we also give the values of the time delays, the initial conditions, and the two saturating actuators as follows:$$\begin{array}{c}\text{(130)}& {\tau }_1=2\mathrm{s},{\tau }_2=4\mathrm{s},\\ \text{(131)}& \zeta t=\frac{1}{2}{\sin t,1-\cos t,\sin t,1-\cos t}^{T}fort\in -\mathrm{4,0},\text{(132)}& sat{u}_{1}t=\begin{array}{cc}50& \mathrm{if}{u}_{1}t\succ 50,\\ u_{1}t& \mathrm{if}-50\le u_{1}t\le 50,\\ -50& \mathrm{if}{u}_{1}t\prec -50,\end{array}sat{u}_{2}t=\begin{array}{cc}80& \mathrm{if}{u}_{2}t\succ 80,\\ u_{2}t& \mathrm{if}-60\le u_{2}t\le 80,\\ -60& \mathrm{if}{u}_{2}t\prec -60.\end{array}\end{array}$$

5.2. Polynomial Model for Studied Benchmark System

Using Taylor series expansions, the nonlinear benchmark model (128), without actuator saturation, can be developed on a third-order polynomial system (25) which is described by the following state equations:$$\begin{array}{}\text{(133)}& \begin{array}{c}\dot{x}t={\displaystyle \sum _{i=1}^3{A}_i{x}^{i}t}+{\displaystyle \sum _{j=1}^3{\displaystyle \sum _{l=1}^2{D}_{jl}{x}^{j}t-{\tau }_l}}+{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{j=1}^3{\displaystyle \sum _{l=1}^2{G}_{ijl}{x}^{i}t\otimes x^{j}t-{\tau }_l+{\displaystyle \sum _{i=1}^2{B}_{i}ut\otimes x^{i}t}+But}}}\hfill \\ yt=Cxt\hfill \\ xt=\zeta tfort\in -\mathrm{4,0}\hfill \end{array},\end{array}$$with$$\begin{array}{c}\text{(134)}& xt=\begin{array}{c}{x}_{1}t\\ x_{2}t\\ x_{3}t\\ x_{4}t\end{array},A_1=\begin{array}{cccc}0& 1& 0& 0\\ \frac{m_{1}gr}{J_1}-\frac{k{r}^2}{4J_1}& 0& \frac{k{r}^2}{4J_1}& 0\\ 0& 0& 0& 1\\ \frac{k{r}^2}{4J_2}& 0& \frac{m_{2}gr}{J_2}-\frac{k{r}^2}{4J_2}& 0\end{array},D_{11}=\begin{array}{cccc}0& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array},D_{12}=\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 1& 0\end{array},B=\begin{array}{cc}\begin{array}{c}0\\ \frac{1}{J_1}\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ \frac{1}{J_2}\end{array}\end{array},\\ C=\begin{array}{c}\begin{array}{cccc}1& 0& 0& 0\end{array}\\ \begin{array}{cccc}0& 0& 1& 0\end{array}\end{array},\\ G_{111}\mathrm{2,5}=\frac{1}{J_1},G_{111}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^2;\\ G_{121}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 4^3;\\ G_{131}\mathrm{2,65}=-\frac{1}{2J_1},G_{131}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^4;\\ G_{112}\mathrm{4,15}=\frac{1}{J_2},G_{112}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^2;\\ G_{122}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 4^3;\\ G_{132}\mathrm{4,235}=-\frac{1}{2J_2},G_{132}i,j=0\hspace{1em}elsewherei=1\dots 4;j=1\dots 4^4;\\ B_{1}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 8;B_2\mathrm{2,1}=-0.25,B_2\mathrm{2,11}=0.25,\\ B_2\mathrm{4,17}=-0.25,B_2\mathrm{2,27}=0.25,B_{2}i,j=0\\ elsewhere\hspace{1em}i=1\dots 4;j=1\dots 32;\\ A_{2}i,j=D_{21}i,j=D_{22}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 4^2;\\ A_3\mathrm{2,1}=-\frac{1}{6}\frac{m_{1}gr}{J_1}-\frac{k{r}^2}{4J_1},A_3\mathrm{2,43}=-\frac{k{r}^2}{24J_1},A_3\mathrm{4,1}=-\frac{k{r}^2}{24J_2};\\ A_3\mathrm{4,43}=-\frac{1}{6}\frac{m_{2}gr}{J_2}-\frac{k{r}^2}{4J_2},A_{3}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^3;\\ D_{31}\mathrm{2,1}=-1,D_{31}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^3;\\ D_{32}\mathrm{4,43}=-1,D_{32}i,j=0\hspace{1em}elsewherei=1\dots 4;j=1\dots 4^3;\\ G_{111}\mathrm{2,5}=\frac{1}{J_1},G_{111}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^2;\\ G_{121}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 4^3;\\ G_{131}\mathrm{2,65}=-\frac{1}{2J_1},G_{131}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^4;\\ G_{112}\mathrm{4,15}=\frac{1}{J_2},G_{112}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 4;j=1\dots 4^2;\\ G_{122}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 4^3;\\ G_{132}\mathrm{4,235}=-\frac{1}{2J_2},G_{132}i,j=0\hspace{1em}elsewherei=1\dots 4;j=1\dots 4^4;\\ B_{1}i,j=0\hspace{1em}for\hspace{1em}i=1\dots 4;j=1\dots 8;B_2\mathrm{2,1}=-0.25,B_2\mathrm{2,11}=0.25,\\ B_2\mathrm{4,17}=-0.25,B_2\mathrm{2,27}=0.25,B_{2}i,j=0\\ elsewhere\hspace{1em}i=1\dots 4;j=1\dots 32;\end{array}$$and for each $i\in \mathrm{2,3}$, $j\in 1\dots 3$, and $l\in \mathrm{1,2}$, $G_{ijl}$ are null matrices with appropriate dimensions.

5.3. Fidelity of the Adopted Polynomial Model

Given the following inputs control:$$\begin{array}{}\text{(135)}& \begin{array}{c}{u}_{1}t=-36x_{1}t-16x_{2}t-20x_{1}t-{\tau }_1-10x_{2}t-{\tau }_1-16x_{1}t-{\tau }_1-12x_{2}t-{\tau }_2,\\ u_{2}t=-36x_{3}t-12x_{4}t-20x_{3}t-{\tau }_1-10x_{4}t-{\tau }_1-16x_{3}t-{\tau }_1-12x_{4}t-{\tau }_2.\end{array}\end{array}$$

Then, the evolution of the state variables in closed loop for both nonlinear (128) and polynomial (133) systems is shown in Figure 2.

[figure(s) omitted; refer to PDF]

From the simulation results, it is proved that the dynamical behavior of the real system is adequately described through the adopted polynomial model affected by multiple time delays in states. This is confirmed by the fact that the behavior of the considered polynomial model can replicate with a enough accuracy the real system one.

5.4. Control Qualification Statements

The aim of the present application is the control of the angular positions ${\theta }_{1}t$ and ${\theta }_{2}t$, in order to force them to track the following reference input vector:$$\begin{array}{}\text{(136)}& y_{c}t=\begin{array}{c}{y}_{c1}t\\ y_{c2}t\end{array}=\begin{array}{c}\sin t+\sin 0.5t\\ \sin t+\sin 0.5t\end{array}.\end{array}$$

For this purpose, we adopt the following parameters matrices of the reference model (40):$$\begin{array}{c}\text{(137)}& E=\begin{array}{cccccccccc}0& 1& 0& 0& 0& 0& 0& 0& 0& 0\\ -480.9& -35.8& -122.4& -4.5& 9518.9& 6513.6& 8825.3& 4682.2& 3019.3& 1175.4\\ 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ -110.8& -4& -477.8& -35.7& 5864.3& 9597.1& 4218.6& 8821.2& 1060.9& 3004\\ 0& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 1\\ -1& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0& 0& 0& 0& 0\end{array},F={\begin{array}{c}\begin{array}{cccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 1& 0\end{array}\\ \begin{array}{cccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\end{array}}^T,\\ Z=\begin{array}{c}\begin{array}{cccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\\ \begin{array}{cccccccccc}0& 0& 1& 0& 0& 0& 0& 0& 0& 0\end{array}\end{array}.\end{array}$$

It should be noted that the chosen number of integrators $w=3$ allows for the output vector $Y_{r}t$ of the reference model to track adequately the considered reference input vector $y_{c}t$.

5.5. Parameters Determination of Block-Pulse Functions Basis

For the selected reference input vector $y_{c}t$, the exact solution $X_{r}t={x_{1,r}t,x_{2,r}t,\dots ,x_{n+wp,r}t}^T$ of the adopted reference model can be obtained from the following relation:$$\begin{array}{}\text{(138)}& X_{r}t=e^{Et}{X}_{r}0+{\displaystyle {\int }_0^t{e}^{Et-\tau }F{y}_c\tau d\tau }.\end{array}$$

In order to determine the parameters of block-pulse functions basis, namely the time interval $0,T$ and the order $N$, it suffices to follow the following steps:

(i) The time interval $0,T$ is sufficiently long to achieve reference model behavior.

(ii) By taking into consideration that each element $x_{i,r}t$ of the state vector $X_{r}t$ is of smooth functions, so from what has been indicated in Section 2.2, it can be approximated efficiently by a finite number of block-pulse functions under the condition that it is high as enough. To satisfy this condition, we compare the exact solution of the reference model $X_{r}t$ with the approximate one $X_{1,rN}{S}_{N}t$, which leads to the suitable choice of the order $N$.

From Figure 3, we can deduce that the agreement is excellent for $N=2^8$ (number of BPFs) and $T=32\mathrm{s}$.

[figure(s) omitted; refer to PDF]

In Table 1, three numerical methods, namely, root mean square error (RMSE), mean squared error (MSE), and mean absolute error (MAE), are used to measure the errors between the exact solution $X_{r}t$ and the approximate one $X_{1,rN}{S}_{N}t$ for $N=2^8$ and $T=32\mathrm{s}$. The errors are defined as follows:$$\begin{array}{}\text{(139)}& E_r=X_{r}t-X_{1,rN}{S}_{N}t={e_{r,1},e_{r,2},e_{r,3},e_{r,4},e_{r,5},e_{r,6},e_{r,7},e_{r,8},e_{r,9},e_{r,10}}^T.\end{array}$$

Table 1

The errors between the exact solution and the approximate one.

Methods	Errors
	$e_{r,1}$	$e_{r,2}$	$e_{r,3}$	$e_{r,4}$	$e_{r,5}$	$e_{r,6}$	$e_{r,7}$	$e_{r,8}$	$e_{r,9}$	$e_{r,10}$
RMSE	0.029 3	0.041 1	0.029 4	0.040 9	0.001 0	0.001 0	0.000 9	0.001 0	0.001 2	0.001 2
MSE	0.000 9	0.001 7	0.000 9	0.001 7	0.000 0	0.000 0	0.000 0	0.000 0	0.000 0	0.000 0
MAE	0.021 6	0.024 0	0.021 6	0.024 0	0.000 7	0.000 7	0.000 7	0.000 7	0.000 8	0.000 8

It appears clear from Table 1 that the errors can be judged to be very satisfactory with respect to the effectiveness expansivity of the exact solution over the considered block-pulse functions basis for $N=2^8$ and $T=32\mathrm{s}$.

5.6. Results and Discussion

For $N=2^8$ (number of BPFs) and $T=32\mathrm{s}$, which permit us to represent adequately the exact solution over the considered BPFs basis of the adopted reference model $X_{r}t$ under the fixed input reference vector $y_{c}t$, the execution of the proposed saturated tracking control approach leads to the following numerical results of the controller’s parameters:$$\begin{array}{c}\text{(140)}& \overline{N}=\begin{array}{cccccc}36103\hfill & 43147\hfill & 33031\hfill & 33945\hfill & 12799\hfill & 8100\hfill \\ 26870\hfill & 68777\hfill & 18904\hfill & 62039\hfill & 6455\hfill & 18875\hfill \end{array},K_1=\begin{array}{cccc}2305.6\hfill & 174.7\hfill & 708.8\hfill & 26.3\hfill \\ 968.8\hfill & 38.5\hfill & 2704\hfill & 208.5\hfill \end{array},\\ K_2=K_{2a},K_{2b},K_{2c},K_{2d},\end{array}$$with:$$\begin{array}{c}\text{(141)}& K_{2a}=\begin{array}{cccc}-12.6331\hfill & -159.6439\hfill & 9.8312\hfill & 161.4933\hfill \\ -15.2547\hfill & -922.4262\hfill & 12.6334\hfill & 921.7037\hfill \end{array},K_{2b}=\begin{array}{cccc}-159.6439\hfill & 13.0846\hfill & 162.3417\hfill & -13.1298\hfill \\ -922.4262\hfill & 8.9307\hfill & 923.6831\hfill & -10.0853\hfill \end{array},\\ K_{2c}=\begin{array}{cccc}9.8312\hfill & 162.3417\hfill & -7.0379\hfill & -164.1832\hfill \\ 12.6334\hfill & 923.6831\hfill & -10.0214\hfill & -922.9536\hfill \end{array},\\ K_{2d}=\begin{array}{cccc}161.4993\hfill & -13.1298\hfill & -164.1832\hfill & 13.1800\hfill \\ 921.7037\hfill & -10.0853\hfill & -922.9536\hfill & 11.2408\hfill \end{array},\\ K_3=K_{3a},K_{3b},K_{3c},K_{3d},K_{3e},K_{3f},K_{3g},K_{3h},K_{3i},K_{3j},K_{3k},\end{array}$$with:$$\begin{array}{c}\text{(142)}& K_{3a}=\begin{array}{cccccc}-72& -883.2& 110.3& 401.6& 118.3& 124.4\\ -167& -832.7& 294.5& 1225.1& -621.2& 40.3\end{array},K_{3b}=\begin{array}{cccccc}272.8& -102.8& 71.1& 118.3& -37.9& 434\\ 623.4& -217.5& 27.1& 1205.2& -126.9& -861.4\end{array},\\ K_{3c}=\begin{array}{cccccc}-398& -102.8& 91.3& 79.9& 118.3& 124.4\\ 589& -150.5& -241.1& 305.9& -621.2& 40.3\end{array},\\ K_{3d}=\begin{array}{cccccc}272.8& -102.8& 124.4& 59.2& -3.4& -23.4\\ 623.4& -217.5& 40.3& -95.3& 152.3& 43.9\end{array},\\ K_{3e}=\begin{array}{cccccc}395.2& -3.4& 192.6& -18.7& -102.8& -23.4\\ -640.7& 152.3& -994.1& -20.8& -150.5& 43.9\end{array},\\ K_{3f}=\begin{array}{cccccc}-18.7& -13.9& 71.1& 118.3& -37.9& 434\\ -20.8& 5.7& 27.1& 1205.2& -126.9& -861.4\end{array},\\ K_{3g}=\begin{array}{cccccc}395.2& -3.4& 192.6& -18.7& -181.2& -1312.6\\ -640.7& 152.3& -994.1& -20.8& 112.6& 1686.2\end{array},\\ K_{3h}=\begin{array}{cccccc}74.9& 1436.2& -271.8& -18.7& -773.7& 42.2\\ -42.5& 1783.0& -414.5& -20.8& -575.4& -111.3\end{array},\\ K_{3i}=\begin{array}{cccccc}-398& -102.8& 91.3& 79.9& -102.8& -23.4\\ 589& -150.5& -241.1& 305.9& -150.5& 43.9\end{array},\\ K_{3j}=\begin{array}{cccccc}-18.7& -13.9& -271.8& -18.7& -773.7& 42.2\\ -20.8& 5.7& -414.5& -20.8& -575.4& -111.3\end{array},\\ K_{3k}=\begin{array}{cccc}79.9& -13.9& 42.2& 53\\ 305.9& 5.7& -111.3& -53.6\end{array},\\ \text{(143)}& L_{11}\mathrm{1,1}=10,L_{11}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 2;j=1\dots 4;L_{12}\mathrm{2,3}=12.5,L_{12}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 2;j=1\dots 4;\\ L_{21}i,j=L_{22}i,j=0\hspace{1em}for\hspace{1em}each\hspace{1em}i=1\dots 2;j=1\dots 16;\\ L_{31}\mathrm{1,1}=-10,L_{31}i,j=0\hspace{1em}elsewhere\hspace{1em}i=1\dots 2;j=1\dots 64;\\ L_{32}i,j=0\hspace{1em}for\hspace{1em}each\hspace{1em}i=1\dots 2;j=1\dots 64.\end{array}$$

Consequently, the following simulation results are obtained:

In Figure 4, the responses of controlled inverted pendulums system using the proposed saturated tracking control method and the considered reference model are plotted.

[figure(s) omitted; refer to PDF]

It can be seen that the time delay polynomial state feedback controllers with triple-action integral, applied to the inverted pendulums system, permit them to achieve their purpose. In fact, the performance of the controlled inverted pendulums system is judged to be very satisfactory with respect to the effectiveness trackability of the reference input vector.

Figure 5 displays the errors between the output vector of the controlled benchmark system $Yt$ and the input reference vector $y_{c}t$, which are denoted as follows:$$\begin{array}{}\text{(144)}& \tilde{e}t=Yt-y_{c}t=\begin{array}{c}{e}_{1}t\\ e_{2}t\end{array}=\begin{array}{c}{\theta }_{1}t-y_{c1}t\\ {\theta }_{2}t-y_{c2}t\end{array}.\end{array}$$

[figure(s) omitted; refer to PDF]

It is clear that the error signals are not significant. In fact, we can see that they converge toward a small neighborhood of the origin, which confirm that the proposed saturated control law computed using the developed approach, permits to have a good tracking performance. This is confirmed by a close to zero steady-state error for the angular positions ${\theta }_{1}t$ and ${\theta }_{2}t$ compared with its time-varying reference inputs $y_{c1}t$ and $y_{c2}t$.

According to the simulation result given in Figure 6, it is seen that the control signals $u_{1}t$ and $u_{2}t$ are always inside, respectively, the interval $-\mathrm{50,50}$ and $-\mathrm{60,80}$. Due to the fact that the imposed limits on the control input by the actuator saturation have been incorporated into the procedure of computing tracking control law, the potential ability of the proposed saturated tracking control method to withstand these constraints can be confirmed from the obtained simulation results.

[figure(s) omitted; refer to PDF]

From what has been stated above, it can be deduced that the adopted strategy of control can guarantee for the closed-loop saturated nonlinear polynomial system subject to multiple time-delayed states an excellently trackability of the reference input vector.

5.7. Simulation Results Given in a Practical Environment

To demonstrate the effectiveness in a practical environment, some measurement noises and external disturbances are imposed on the system. The measurement noises $m_n$ are assumed as $m_n=0.01\mathrm{rand}\mathrm{4,1}$, and the external disturbances $e_{d}t$ are set as$$\begin{array}{}\text{(145)}& e_{d}t={0.2\sin t,0.2\cos t,0.2\sin t,0.2\cos t}^T.\end{array}$$

The corresponding simulation results are shown in Figures 7 and 8. We can deduce that the proposed strategy of control permits to conserve an acceptable tracking performance even in the presence of some measurement noises and external disturbances. Hence, this can be evaluated as a second advantage of the investigated approach in this study.

[figure(s) omitted; refer to PDF]

5.8. Practical Stability Test

The obtained gains verified that all eigenvalues of matrix ${\mathrm{ℳ}}_1$ have a strictly negative real part, and$$\begin{array}{}\text{(146)}& e^{{\mathrm{ℳ}}_{1}t}<360e^{0.4t},\end{array}$$which corresponds to inequality (98) with ${\lambda }_1=360$ and ${\omega }_1=-0.4$. Then, Theorem 1 provides an estimation of the radius of practical stability of the studied system $R_0=0.00000000144$.

6. Conclusion

In this article, we have built up an original control scheme for a class of nonlinear polynomial systems subject to control input constraints. These systems are characterized by their complex and high nonlinearities as well as fixed and well-defined time delays. The design methodology makes full use of a memory nonlinear state feedback control formalism, accompanied by a compensator gain built through the integral actions technique. In the initial step, the closed-loop augmented form is determined using the augmented error modeling technique. In the subsequent step, an algebraic system with constraints of controller parameters is derived so as to solve the tracking control problem in the sense of constrained least squares optimization methodology. In the last step, the system’s closed-loop stability is assessed through providing new sufficient conditions for the model’s practical stability. As an original performance of the synthesized technique, it is straightforward to estimate the enlarged region of stability.

The simulation analysis has demonstrated the viability of the developed control method in enhancing the performance of the studied benchmark model output tracking. Indeed, it is plainly observed throughout the simulation study that the suggested memory polynomial controllers with integrals action can quickly attenuate the steady-state error, as a result of a good tracking performance, even in the presence of actuator saturation.

Future works will stretch out the developed technique to deal with saturated tracking control design for nonlinear polynomial time-varying systems affected with both unstructured uncertainties and time delays in states.

References

[1] I. Karafyllis, M. Malisoff, F. Mazenc, P. Pepe, "Erratum to: recent results on nonlinear delay control systems," Recent Results on Nonlinear Delay Control Systems, vol. 4,DOI: 10.1007/978-3-319-18072-4_21, 2016.

[2] Y. Tu, H. Fang, Y. Yin, S. He, Reinforcement Learning-Based Nonlinear Tracking Control System Design via LDI Approach with Application to Trolley System, 2021.

[3] P. Cheng, S. He, X. Luan, F. Liu, "Finite-region asynchronous Hcontrol for 2D Markov jump systems," Automatica, vol. 129,DOI: 10.1016/j.automatica.2021.109590, 2021.

[4] M. Cai, Z. Xiang, "Adaptive finite-time consensus tracking for multiple uncertain mechanical systems with input saturation," International Journal of Robust and Nonlinear Control, vol. 27, pp. 1653-1676, DOI: 10.1002/rnc.3741, 2017.

[5] S. Li, Z. Xiang, "Adaptive prescribed performance control for switched nonlinear systems with input saturation," International Journal of Systems Science, vol. 49 no. 1, pp. 113-123, DOI: 10.1080/00207721.2017.1390706, 2018.

[6] B. Mirkin, J. Haddad, Y. Shtessel, "Tracking with asymptotic sliding mode and adaptive input delay effect compensation of nonlinearly perturbed delayed systems applied to traffic feedback control," International Journal of Control, vol. 89 no. 9, pp. 1890-1903, DOI: 10.1080/00207179.2016.1174340, 2016.

[7] I. Estrada-Sánchez, M. Velasco-Villa, H. Rodríguez-Cortés, "Prediction-based control for nonlinear systems with input delay," Mathematical Problems in Engineering, vol. 2017,DOI: 10.1155/2017/7415418, 2017.

[8] D. Zhai, C. Xi, J. Dong, Q. Zhang, "Adaptive fuzzy fault-tolerant tracking control of uncertain nonlinear time-varying delay systems," IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 50 no. 5, pp. 1840-1849, DOI: 10.1109/tsmc.2018.2789441, 2020.

[9] W. Liu, C.-C. Lim, P. Shi, S. Xu, "Observer-based tracking control for MIMO pure-feedback nonlinear systems with time-delay and input quantisation," International Journal of Control, vol. 90 no. 11, pp. 2433-2448, DOI: 10.1080/00207179.2016.1250162, 2017.

[10] L. Xue, Z. Liu, Z. Sun, W. Sun, "New results on robust tracking control for a class of high-order nonlinear time-delay systems," International Journal of Systems Science, vol. 50 no. 10, pp. 2002-2014, DOI: 10.1080/00207721.2019.1646346, 2019.

[11] Y. Sun, Q. Zhao, C. Zhang, J. Wu, "Adaptive dynamic surface control design with asymptotic tracking performance for a class of uncertain input-delayed systems," Circuits, Systems, and Signal Processing, vol. 39 no. 7, pp. 3512-3531, DOI: 10.1007/s00034-019-01335-3, 2020.

[12] C.-S. Chiu, T.-S. Chiang, "H ∞ output-feedback fuzzy proportional-integral control of fully delayed input/output systems," ISA Transactions, vol. 66, pp. 22-31, DOI: 10.1016/j.isatra.2016.10.003, 2017.

[13] Y. H. Choi, S. J. Yoo, "Neural-networks-based adaptive quantized feedback tracking of uncertain nonlinear strict-feedback systems with unknown time delays," Journal of the Franklin Institute, vol. 357 no. 15, pp. 10691-10715, DOI: 10.1016/j.jfranklin.2020.08.046, 2020.

[14] Z. Rahmani, B. Baigzadehnoe, B. Rezaie, "Tracking control of a class of nonlinear systems with output delay based on adaptive fuzzy dynamic surface control," International Journal of Systems Science, vol. 51 no. 7, pp. 1280-1306, DOI: 10.1080/00207721.2020.1758827, 2020.

[15] Y. Li, Z. Lin, Stability and Performance of Control Systems with Actuator Saturation, 2018.

[16] M. Hussain, M. Rehan, "Nonlinear time-delay anti-windup compensator synthesis for nonlinear time-delay systems: a delay-range-dependent approach," Neurocomputing, vol. 186, pp. 54-65, DOI: 10.1016/j.neucom.2015.12.078, 2016.

[17] A. Akram, M. Hussain, N. us Saqib, M. Rehan, "Dynamic anti-windup compensation of nonlinear time-delay systems using LPV approach," Nonlinear Dynamics, vol. 90 no. 1, pp. 513-533, DOI: 10.1007/s11071-017-3678-8, 2017.

[18] M. Hussain, M. Rehan, C. K. Ahn, Z. Zheng, "Static anti-windup compensator design for nonlinear time-delay systems subjected to input saturation," Nonlinear Dynamics, vol. 95 no. 3, pp. 1879-1901, DOI: 10.1007/s11071-018-4666-3, 2019.

[19] P. Selvaraj, B. Kaviarasan, R. Sakthivel, H. R. Karimi, "Fault‐tolerant SMC for Takagi-Sugeno fuzzy systems with time‐varying delay and actuator saturation," IET Control Theory & Applications, vol. 11 no. 8, pp. 1112-1123, DOI: 10.1049/iet-cta.2016.0845, 2017.

[20] H. Gassara, A. El Hajjaji, M. Chaabane, "Control of time delay polynomial fuzzy model subject to actuator saturation," International Journal of Fuzzy Systems, vol. 18 no. 5, pp. 763-772, DOI: 10.1007/s40815-015-0133-1, 2016.

[21] H. Khalil, Nonlinear Systems, 2002.

[22] M. Rehan, N. Iqbal, K.-S. Hong, "Delay-range-dependent control of nonlinear time-delay systems under input saturation," International Journal of Robust and Nonlinear Control, vol. 26 no. 8, pp. 1647-1666, DOI: 10.1002/rnc.3372, 2016.

[23] H. Min, S. Xu, Q. Ma, B. Zhang, Z. Zhang, "Composite-observer-based output-feedback control for nonlinear time-delay systems with input saturation and its application," IEEE Transactions on Industrial Electronics, vol. 65 no. 7, pp. 5856-5863, DOI: 10.1109/tie.2017.2784347, 2018.

[24] Q. Liu, H. Li, X. Zhang, "Output feedback tracking control for lower‐triangular nonlinear time‐delay systems with asymmetric input saturation," International Journal of Robust and Nonlinear Control, vol. 30 no. 8, pp. 3471-3483, DOI: 10.1002/rnc.4953, 2020.

[25] H. Wang, B. Chen, C. Lin, Y. Sun, "Observer-based adaptive fuzzy tracking control for a class of MIMO nonlinear systems with unknown dead zones and time-varying delays," International Journal of Systems Science, vol. 50 no. 3, pp. 546-562, DOI: 10.1080/00207721.2018.1563224, 2019.

[26] H. Dastres, B. Rezaie, B. Baigzadehnoe, "Neural-network-based adaptive backstepping control for a class of unknown nonlinear time-delay systems with unknown input saturation," Neurocomputing, vol. 398, pp. 131-152, DOI: 10.1016/j.neucom.2020.02.070, 2020.

[27] A. Moradvandi, M. Shahrokhi, S. A. Malek, "Adaptive fuzzy decentralized control for a class of MIMO large-scale nonlinear state delay systems with unmodeled dynamics subject to unknown input saturation and infinite number of actuator failures," Information Sciences, vol. 475, pp. 121-141, DOI: 10.1016/j.ins.2018.09.052, 2019.

[28] X. Zhao, H. Yang, W. Xia, X. Wang, "Adaptive fuzzy hierarchical sliding-mode control for a class of MIMO nonlinear time-delay systems with input saturation," IEEE Transactions on Fuzzy Systems, vol. 25 no. 5, pp. 1062-1077, DOI: 10.1109/tfuzz.2016.2594273, 2017.

[29] J. Brewer, "Kronecker products and matrix calculus in system theory," IEEE Transactions on Circuits and Systems, vol. 25 no. 9, pp. 772-781, DOI: 10.1109/tcs.1978.1084534, 1978.

[30] M. M. Belhaouane, N. B. Braiek, "Design of stabilizing control for synchronous machines via polynomial modelling and linear matrix inequalities approach," International Journal of Control, Automation and Systems, vol. 9 no. 3, pp. 425-436, DOI: 10.1007/s12555-011-0301-5, 2011.

[31] B. Iben Warrad, M. K. Bouafoura, N. Benhadj Braiek, "Tracking control synthesis of nonlinear polynomial systems," .

[32] H. Rtibi, S. Elloumi, N. Benhadj Braiek, "Development of a decentralized nonlinear controller for a class of uncertain polynomial interconnected systems: application for a large scale power system," Transactions of the Institute of Measurement and Control, vol. 41 no. 8, pp. 2236-2249, DOI: 10.1177/0142331218796150, 2018.

[33] A. Feydi, S. Elloumi, C. Jammazi, N. Benhadj Braiek, "Decentralized finite-horizon suboptimal control for nonlinear interconnected dynamic systems using SDRE approach," Transactions of the Institute of Measurement and Control, vol. 41 no. 11, pp. 3264-3275, DOI: 10.1177/0142331218820916, 2019.

[34] B. Iben Warrad, M. K. Bouafoura, N. Benhadj Braiek, "Tracking control design for nonlinear polynomial systems via augmented error system approach and block pulse functions technique," Kybernetika, vol. 5 no. 5, pp. 831-851, DOI: 10.14736/kyb-2019-5-0831, 2019.

[35] B. Iben Warrad, M. K. Bouafoura, N. Benhadj Braiek, "Static output tracking control for non-linear polynomial time-delay systems via block-pulse functions," Journal of the Chinese Institute of Engineers, vol. 41 no. 3, pp. 194-205, DOI: 10.1080/02533839.2018.1459849, 2018.

[36] Y. Xu, J. Zhang, F. Liao, "Augmented error system approach to control design for a class of neutral systems," Advances in Difference Equations, vol. 270,DOI: 10.1186/s13662-015-0596-2, 2015.

[37] W. L. Yen-Ping Shih, Y. P. Shih, "Shift Walsh matrix and delay-differential equations," IEEE Transactions on Automatic Control, vol. 23 no. 6, pp. 1023-1028, DOI: 10.1109/tac.1978.1101888, 1978.

[38] Z. H. Jiang, W. Schaufelberger, Block-Pulse Functions and Their Applications in Control Systems, 1992.

[39] S.-L. Chen, H.-C. Lai, K.-C. Ho, "Identification of linear time varying systems by Haar wavelet," International Journal of Systems Science, vol. 37 no. 9, pp. 619-628, DOI: 10.1080/00207720600774255, 2006.

[40] A. Deb, S. Roychoudhury, G. Sarkar, "System identification: parameter estimation of transfer function," Analysis and Identification of Time-Invariant Systems, Time-Varying Systems, and Multi-Delay Systems using Orthogonal Hybrid Functions, vol. 46, pp. 331-355, DOI: 10.1007/978-3-319-26684-8_15, 2016.

[41] B. Iben Warrad, M. K. Bouafoura, N. Benhadj Braiek, "Combined constrained robust least squares approach and block-pulse functions technique for tracking control synthesis of uncertain bilinear systems with multiple time-delayed states under bounded input control," Mathematical Problems in Engineering, vol. 2020,DOI: 10.1155/2020/7186928, 2020.

[42] H. R. Marzban, F. Malakoutikhah, "Solution of delay fractional optimal control problems using a hybrid of block-pulse functions and orthonormal Taylor polynomials," Journal of the Franklin Institute, vol. 356 no. 15, pp. 8182-8215, DOI: 10.1016/j.jfranklin.2019.07.010, 2019.

[43] O. Balatif, M. Rachik, E. Labriji, Z. Rachik, "Optimal control problem for a class of bilinear systems via block pulse functions," ”Journal of Mathematical Control and Information, vol. 34 no. 3, pp. 973-986, DOI: 10.1093/imamci/dnw005, 2016.

[44] C. Hwang, Y. P. Shih, "Optimal control of delay systems via block pulse functions," Journal of Optimization Theory and Applications, vol. 45 no. 1, pp. 101-112, DOI: 10.1007/bf00940816, 1985.

[45] V. Lakshmikantham, S. Leela, A. A. Martynyuk, Practical Stability of Nonlinear Systems, 1990.

[46] W.-J. Wang, B.-S. Chen, "Stability of large-scale systems with saturating actuators," International Journal of Control, vol. 47 no. 3, pp. 827-850, DOI: 10.1080/00207178808906056, 1988.

[47] H. Sun, L. Hou, G. Zong, X. Yu, "Adaptive decentralized neural network tracking control for uncertain interconnected nonlinear systems with input quantization and time delay," IEEE Transactions on Neural Networks and Learning Systems, vol. 31 no. 4, pp. 1401-1409, DOI: 10.1109/tnnls.2019.2919697, 2020.