# Robust Local Stabilization of Discrete-Time Systems with Time-Varying State Delay and Saturating Actuators

Authors: J. V. V. Silva; L. F. P. Silva; I. Rubio Scola; V. J. S. Leite
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The robust local stabilization of uncertain discrete-time systems with time-varying state delayed and subject to saturating actuators is investigated in this work. A convex optimization method is proposed to compute robust state feedback control law such that the uncertain closed-loop is locally asymptotically stable if the initial condition belongs to an estimate of the region of attraction for the origin. The proposed procedure allows computing estimates of the region of attraction through the intersection of ellipsoidal sets in an augmented space, reducing the conservatism of the estimates found in the literature. Also, the conditions can handle the amount of delay variation between two consecutive samples, which is new in the literature for the discrete-time case. Although the given synthesis conditions are delay dependent, the proposed control law is delay independent, yielding to easier real time implementations. A convex procedure is proposed to maximize the size of the set of safe initial conditions. Numerical examples are provided to illustrate the effectiveness of our approach and also to compare it with other conditions in the literature.


---

## Body

## 1. Introduction

We are concerned with the robust local stabilization of discrete-time systems with time-varying delay in the state subject to saturating actuators. Note that discrete-time systems with delayed state have been widely investigated in the literature and, in most of the cases, the approach based on Lyapunov-Krasovskii function has been used [1, 2]. Recent advances on discrete-time delayed state systems include several works where an improved (double) summation inequality is employed or some tighter sum over bounds (for instance, Wirtinger’s inequalities) yielding less conservative analysis conditions [3–6]. Although handling time-varying delay, none of these conditions are able to take into account bounds on the amount of variation of the delay between two successive instants. The estimate of region of safe initial conditions is also considered in [7] through augmented delay-free state representation by means of D-invariant and D-contractive sets. Besides this work, we can find similar studies in [8–11], but the saturating actuator is not addressed. However, it is interesting to reinforce that the delay is constant and known. In fact, differently from the continuous-time case, all conditions we have found in the literature assume that the time-varying delay, τk, can jump from its lower bound, usually τk=τ_=1 at one instant, to its maximum value, τ-, in the next instant, i.e., τk+1=τ-. Therefore, it is expected that such an assumption leads to conservative analysis and synthesis conditions.Another relevant issue in control systems concerns saturating actuators, which can reduce performance or even leads the closed-loop system to unstable behavior [12, 13]. Therefore, a fundamental (and difficult) task is to estimate the region of safe initial conditions, i.e., an estimate of the set of initial conditions such that the closed-loop trajectories starting from there remain bounded and go asymptotically to the origin (the local equilibrium point). Such a region is not easily characterized specially for delayed state systems, because it is required to characterize a sequence of vectors and not only a vector for a complete initial condition (thus ensuring uniqueness of solutions). A few works are found in the literature handling saturating actuators in the context of discrete-time delayed systems such as [14–19]. The description of the saturating function can be done by a polytopic based representation [15–17] or through a sector condition approach as in [14, 17–19] in the context of fuzzy systems. The estimates of the region of initial conditions may be characterized in different ways as we can verify in literature. One of the first approaches was to represent such a region by ball in Rn. In this case, both the current state and the delayed ones are required to belong to the ball [15, 16]. A more general approach has been considered in [14, 17], where two balls in Rn have been used: one where the current and delayed states must belong and the other to bound the norm of the time variation of the sequence of initial states. By a special choice of the radius of later ball, one can recover the representation in [15, 16]. This is discussed later in this work. A more recent approach used in [18, 19] also employs two sets to characterize the region of safe initial conditions: an ellipsoidal set in Rn that contains the current state vector and a ball in Rn, where all delayed states must belong. However, none of these works can handle bounds on the amount of variation of the delay between two successive samples.Inspired by [20], we introduce a new characterization of the region of initial conditions by using an augmented space where each point corresponds to a sequence of vectors corresponding to an initial condition. The delayed system is represented by an augmented delay-free switched system with the delay as the switch function. The saturating actuator is handled by the generalized sector condition [21]. As argued by the authors of [20], this approach has a great appealing because an augmented Lyapunov candidate function is employed, which is equivalent to a quite general set of Lyapunov-Krasovskii candidate functions. Thanks to this approach, we can handle the case with null delay, which is usually avoided in the literature. See [3–6, 14–20, 22]. Another highlight of our work is the possibility of feedback of the current state or the current state with all delayed ones or with only the delayest one. These control laws are all delay independent in the sense that they dismiss the knowledge of the current state value, which facilitate their use in control applications. On the other hand, the proposed conditions are delay-dependent ones. Two numerical examples are presented to illustrate our approach. Also, we compare our approach with others found in the literature showing that our conditions lead to less conservative results.Notation. The symbol ⋆ represents the transpose of the blocks symmetric with respect to the main diagonal of square symmetric matrices with real entries. Matrices I and 0 stand for identity and null matrices, respectively, with appropriate dimensions. The sets of real and integer numbers are noted by R and Z, respectively. The set of matrices with real entries and with dimensions m×n is Rm×n. A block diagonal matrix with the matrices M1 and M2 in the diagonal blocks is noted as block-diag{M1,M2}. The ith line of the matrix (or vector) M is denoted by M(i). If φ is a sequence of vectors, then its jth element is noted by [φ]j. The sequence φd,k∈Eϕ, with Eϕ=E1×E2×⋯×Ed+1, Ej⊆Rn, has its jth element given by [φd,k]j=xk+j-(d+1)∈Ej, j=1,…,d+1. Thus, φd,k={xk-d,xk-(d-1),…,xk}. φd,k=supj=1,…,d+1φd,kj and · refers to the Euclidean norm.


## 2. Problem Statement

Consider the class of uncertain discrete-time systems with state delay and saturating actuators described by(1)xk+1=Aαxk+Adαxk-τk+Bαsatvk,where xk=x(k)∈Rn is the state vector, vk is the linear control signal, sat(vk) is a decentralized vectorial saturation function: sat(vk)(i)=vk(i), if vki≤v-i, sat(vk)(i)=sign(vk)v-(i), otherwise, with v-(i), i=1,…,p, being the symmetrical bound on the values achieved by the ith actuator. The delay in the state is time-varying and represented by τk verifying(2)0≤τk≤τ-.The matrices of the system are denoted by Υ(α)≡A(α)Ad(α)B(α)∈Rn×(2n+p)=AAdB(α) with the time-invariant uncertain parameter α∈Ξ where (3)Ξ=α∈RN:∑i=1Nαi=1,αi≥0,i=1,…,N,the matrices Υ(α) belonging to the polytope S with N known vertices (4)S=Υα∈Rn×2n+p∣Υα=∑i=1NαiΥi,α∈Ξ,and Υi≡AiAd,iBi=AAdBi. The uniqueness of the solution of (1) is to ensure assuming an initial condition given by the sequence φτ-,0. Note that, due to the saturating actuator, an estimate of the initial conditions is also required to ensure the closed-loop stability of (1) with any control law yielding vk, as discussed later in this work. A contribution of this approach is to consider the time variation of the delay τk bounded as follows: (5)τk+1-τk≤Δτmax≤τ-.Therefore, the transitions of τk between two successive instants k and k+1 are such that if τk=τ, then 0≤τ-Δτmax≤τk+1≤τ+Δτmax≤τ-. We define C(τ) as the set that contains all possible values that can be assumed by the delay in the next sample, given its current value τk=τ: (6)Cτ=τ+∈I⊂Z:I≜max⁡τ_,τ-Δτmax,min⁡τ-,τ+Δτmax,where τ+=τk+1. Then, for τ-=5, τk=1, and Δτmax=2, we have C(τ)={0,1,2,3}, which is the set of all possible values that can be assumed by τk+1.The following delay independent control law is investigated to robust stabilize system (1): (7)vk=Kx-k,where K=K0K1⋯Kτ-∈Rp×(τ-+1)n and x-k=xkTxk-1T⋯xk-τ-TT∈R(τ-+1)n. Two particular structures of interest for (7) are(i)Ki=0, i=1,…,τ-, where only the current state is used in the feedback, and(ii)Ki=0, i=1,…,τ--1, where only the current state and the delayed one are employed to stabilize the closed-loop system. Note that control law (7) dismisses the real time knowledge of τk and, thus, it is called delay independent. Moreover, its application in real process is easier than control laws depending on the delay, as, for instance, vk=K0xk+Kτ-xk-τk. Observe that the notion of stability in state delayed systems requires that all elements of the sequence φτ-,k goes to the origin as k→∞ [1, 26].Because of the saturation of the actuator, there are initial conditionsφτ-,0 that may diverge from the origin even if the nonsaturated closed-loop system is stable; i.e., the elements of the sequence φτ-,k do not go to zero as k→∞. Therefore, we adapt the definition of region of attraction, RA, from [13] as follows.Definition 1 (region of attraction,RA).
The region of attractionRA of the origin for system (1) is the set of all sequences φτ-,0 whose respective trajectories converge asymptotically to the origin, i.e., limk→∞⁡φτ-,k={0,…,0}.Therefore, under saturated actuators it is necessary to determine the setRA or a subset of it, LV⊂RA, as large as possible, such that the local stability can be ensured. As it is widely known in the literature, the determination of RA is very difficult task [13, p. 15] and the determination of an estimate LV⊂RA usually leads to a numerical tractable conditions. In this sense, the following problem is proposed.Problem 2 (robust local stabilization).
Determine the robust control gain matrixK and characterize an estimate LV of the region of attraction RA such that the trajectories of system (1)–(5) under control law (7) and initial conditions given by all sequences φτ-,0 with elements taken from LV∈RA converge to the origin, for all α∈Ξ.Note that it is embedded in Problem2 the issue of consider a bound on the time variation of the delay, which is also a contribution of this work.


## 3. Auxiliary Results

In what follows, we present some auxiliary results used to obtain our main contributions. Inspired by [20], the system (1) is represented as an augmented switched system, where the delay τk is the switching function: (8)x-k+1=Aα,τkx-k+Bαsatvk,where (9)
A
α
,
τ
k
=
A
α
+
Γ
0
,
α
,
τ
k
Γ
1
,
α
,
τ
k
⋯
Γ
τ
-
-
1
,
α
,
τ
k
Γ
τ
-
,
α
,
τ
k
I
n
τ
-
0
n
τ
-
×
n
,Γi,α,τk=Ad(α) if i=τk and Γi,α,τk=0n×n, otherwise, B(α)=BT(α)0p×nτ-T, and x-k given after (7). Using the control law (7) into system (8), we can get x-k+1=A(α,τk)x-k+B(α)sat(Kx-k). With the aid of the dead zone function, Ψ(vk)=vk-sat(vk), we sum and subtract B(α)Kx-k in the right side of (8) to get x-k+1=(A(α,τk)+B(α)K)x-k-B(α)Ψ(Kx-k). To simplify the notation, we assume A~(α,τk)=A(α,τk)+B(α)K and the closed-loop system can be rewritten as (10)x-k+1=A~α,τkx-k-BαΨKx-k.Lemma 3 (generalized sector condition [21]).
Assumevk given by (7), v-=-v_>0 and a matrix G∈Rp×(τ-+1)n, such that (11)S=x-k∈Rτ-+1n:Kl-Glx-k≤v-l,forl=1,…,p..If x-k∈S, then Ψ(Kx-k)TT[Ψ(Kx-k)-Gx-k]≤0 is satisfied for any positive definite diagonal matrix T∈Rp×p.The stability here is investigated by means of Lyapunov’s Theorem [1]. A candidate function V(φτ-,k,α,τk)>0 is called a Lyapunov function if there exist positive constants ϵi, for i=0,1,2, such that ϵ0x-k2≤V(φτ-,k,αk,τk)≤ϵ1x-k2, and ΔV(φτ-,k,αk,τk)=V(φτ-,k+1,α,τk+1)-V(φτ-,k,α,τk)≤-ϵ3x-k2<0. In this work, we use the following augmented Lyapunov candidate function: (12)Vφτ-,k,α,τk=x-kTWα,τk-1x-k,where (13)Wα,τk=∑i=1NαiWiτk,α∈Ξ.It is worth noting that (12)-(13) encompass several Lyapunov-Krasovskii candidate functions (see [20]) and, thus, can yield less conservative stability and synthesis conditions. Also note that the mapping provided by (12)-(13) slightly differs from that proposed in [20] because of the inversion of W(α,τk).The estimate ofRA is performed here through the level set LV associated with a contractive set determined by the augmented Lyapunov function.Definition 4 (level set,LV).
Suppose that the function given by (12)-(13) is a Lyapunov function for system (8) with (3)-(5). The level set LV associated with (12)-(13) is defined as the intersection of the ellipsoidal sets E(W(α,τk)={x-k∈Rn(τ-+1):x-kTW(α,τk)-1x-k≤1}, for all α∈Ξ and τk∈{0,…,τ-}.It is worth saying that the computation of the setLV can be done as indicated in the next lemma, which is adapted from [27].Lemma 5.
Suppose a L-K function given by (12)-(13). The level set LV can be calculated as (14)LV=⋂τ∈0,…,τ-α∈ΞEWα,τk-1=⋂τ∈0,…,τ-i∈1,…,NEWiτk-1,with EWi(τ)-1, for i={1,…,N} and τ∈{0,…,τ-}.The level setLV from Lemma 5 is positively invariant and, thus, a trajectory started inside it goes to the origin. Therefore, LV is an estimate of RA, providing a set of safe initial conditions for the elements of the sequence φτ,0: each point in the space of LV corresponds to a complete sequence of initial conditions. Such a set characterization differs from those used in [14–19], resulting in a more general region. As shown in the numerical examples, our approach leads to less conservative estimates of the region of attraction, thus yielding a larger set of admissible initial conditions with respect to other conditions in the literature.


## 4. Robust Local Stabilization

In what follows, we present our main contribution, providing a solution to Problem2.Theorem 6.
For a given scalarv-, if there exist symmetric and positive definite matrices 0<Wi,τ=Wi,τT∈R(τ-+1)n×(τ-+1)n, τ=0,…,τ-, i=1,…,N, a diagonal positive definite matrix 0<S∈Rp×p, and matrices U∈R(τ-+1)n×(τ-+1)n, Y∈Rp×(τ-+1)n, and Z∈Rp×(τ-+1)n such that the following LMIs are satisfied: (15)-Wi,jAi,τU+BiY-BiS⋆Wi,τ-U-UTZT⋆⋆-2S<0,i=1…,N;τ=0,…,τ-;j∈Cτ,(16)Wi,τ-UT-UYlT-ZlT⋆-v-l2≤0,i=1,…,N;τ=0,…,τ-;l=1,…,p,with (17)
A
i
,
τ
=
A
i
+
Γ
0
,
i
,
τ
Γ
1
,
i
,
τ
⋯
Γ
τ
-
-
1
,
i
,
τ
Γ
τ
-
,
i
,
τ
I
n
τ
-
0
n
τ
-
×
n
,
B
i
=
B
i
0
n
τ
-
×
p
,and, for m=0,…,τ-, Γm,i,τ=Ad,i if m=τ otherwise Γm,i,τ=0n×n, then the control law (7), with (18)K=YU-1,ensures that the origin of the closed-loop system consisting of (1)-(7) is robustly asymptotically stable for every initial condition φτ-,0 whose elements are taken in the set LV⊂RA computed as in Lemma 5. Moreover, (12)-(13) are a Lyapunov function for this system.Proof.
If (15) is verified, then, by its block (1,1) we have the positivity of Wi,j ensured and thus, by (12)-(13) we have that V(x-k)>0 for x-k≠0. Moreover, if we multiply the left-hand side of (15) by αi and sum it up for i=1,…,N, replace Y and Z by KU and GU, respectively, and rename j as τ+ (meaning the value of τk in the next instant), we get (19)-Wα,τ+Aα,τU+BαKU-BαS⋆Wα,τ-U-UTGUT⋆⋆-2S<0.By block (2,2) of (19) we have ensured the regularity of U and from [28], we use the fact -UTW(α,τ)-1U≤W(α,τ)-UT-U; pre- and postmultiply the resulting inequality by T1=block-diag{I,U-1,S-1} and its transpose to get (20)-Wα,τ+Aα,τ+BαK-Bα⋆-Wα,τ-1GTS-1⋆⋆-2S-1<0.We call A~(α,τ)=A(α,τ)+B(α)K and apply Schur’s complement to obtain (21)-Wα,τ-1GTS-1⋆-2S-1+A~α,τT-BαTWα,τ+-1A~α,τ-Bα<0.By pre- and postmultiplying (21) by the augmented state XkT=x-kTΨ(Kx-k)T and Xk, where ΨKx-k is the decentralized dead zone function, we can ensure that –x-kTW(α,τ)-1x-k + 2Ψ(Kx-k)TS-1Gx-k − 2Ψ(Kx-k)TS-1Ψ(Kx-k) + x-kTA~(α,τ)T-Ψ(Kx-kT)B(α)TW(α,τ+)-1 [A~(α,τ)x-k-B(α)Ψ(Kx-k)]<0. By (10), we can replace A~(α,τ)x-k-B(α)Ψ(Kx-k) by x-k+1 and reorganize the expression 2Ψ(Kx-k)TS-1Gx-k-2Ψ(Kx-k)TS-1Ψ(Kx-k) as -2Ψ(Kx-k)TS-1(Ψ(Kx-k)-Gx-k), which allows to verify Ωk≡x-k+1TW(α,τ+)-1x-k+1-x-kTW(α,τ)-1x-k-2Ψ(Kx-k)TS-1(Ψ(Kx-k)-Gx-k)<0. By (12)-(13), we can write ΔV(φτ-,k,α,τk)=x-k+1TW(α,τk+1)-1x-k+1-x-kTW(α,τk)-1x-k and calling T=S-1 we can conclude that ΔV(φτ-,k,α,τk)-2Ψ(Kx-k)TT(Ψ(Kx-k)-Gx-k)≤Ωk<0. Therefore, we conclude that the feasibility of (15) ensures the positivity of (12) with (13), the negativity of ΔV(φτ-,k,α,τk), and, by Lemma 3, we have that x-k∈S. By Lyapunov’s Theorem [1], the local stability of the uncertain discrete-time system with state delay and bounded delay variation on time (1)-(7) is ensured with the control gain given by (18) and always the state trajectories of the closed-loop system evolve inside the set S with ϵ0=σx-k2, where σ is small enough constant, ϵ1=maxτ=0,…,τ-i=1,…,N⁡λmaxWi,τ-1, and for some ϵ3>0.
If (16) is also verified, then we multiply it by αi and sum it up for i=1,…,N and replace the block (1,1) of the resulting matrix by -UTW(α,τ)-1U which is less than or equal to W(α,τ)-UT-U. Then, we pre- and postmultiply the resulting inequality by T2=block-diag{U-T,1} and its transpose, respectively, and we use (18) to get (22)-Wα,τ-1Kl-GlT⋆-v-l2≤0.By using Schur’s complement and pre- and postmultiplying the resulting expression by x-kT and x-k we get (23)-x-kTWα,τk-1x-k+v-l-2x-kTKl-GlTKl-Glx-k≤0,whose consequence is, by Lemma 5, that LV⊆S and any trajectory of the closed-loop system starting in LV remains in S.Remark 7.
Theorem6 provides a solution for Problem 2 and up to the authors’ knowledge this is the first time that the variation of the delay is taken into account in the context of discrete-time systems with delay. Conditions (15)-(16) encompass the case when Δτmax⁡=τ--1, i.e., the case where it is admissible that the delay value can jump from its minimal value (τ=0) to its maximal value (τ=τ-) from one instant k to another k+1. Also note that Theorem 6 can handle null delays and, thus, it is more general than the conditions proposed in [3–7, 7–19], which work only for delays equal to or greater than 1.Remark 8.
The precisely known case can be addressed by conditions in Theorem6 just by suppressing α and, as a consequence, also the index i on LMIs (15)-(16). So, it is possible to show that the Lyapunov candidate function (12)-(13) can be rewritten in a form closer to that proposed by [20], i.e., without the inversion of W(α,τk). The numerical complexity in this case is, obviously, smaller than that of the conditions of Theorem 6.Remark 9.
Some structures of interest for the gain matrixK in (7) can be recovered from (15)-(18) by imposing some special structures on matrices Y and U. To get K=K00⋯0, it is enough to impose (24)Y=Y1,p×n0,U=U1,n×n0U2,nτ-×nU3,nτ-×nτ-.On the other hand, to get the gain K=K00⋯0Kτ-, it is necessary to constrain the structure of Y and U as (25)Y=Y1,p×n0Y2,p×n,U=U1,n×n00U2,nτ--1×nU3,nτ--1×nτ--1U4,nτ--1×nU5,n×n0U6,n×n.An interesting issue is that, in all cases, no structure is imposed to the Lyapunov candidate matrices Wi,τ, i=1,…,N and τ=0,…,τ-, allowing a larger set of search for gain K.Remark 10.
Another way to estimate the region of attraction is through the setLV1={φτ-,0:φτ-,0<ρ0} as it is done in [15, 16]. In this case, the estimate of RA corresponds to τ+1 balls in Rn with radius ρ, one for each delayed state x-k, k=0,…,τ, and LV1⊆LV⊆RA. In [14, 17–19], the estimate of RA is given in terms of LV2={φτ-,0:φτ-,0<ρ1,Δφτ-,0<ρ2}, with Δφτ-,0=φτ-,1-φτ-,0. Note that the set LV1 implies Δφτ-,0<φτ-,1+φτ-,0<2ρ0. In all of them, the parameters ρi, for i=0,1,2, are computed from maximum eigenvalues of the Lyapunov-Krasovskii functions matrices and upper value of the time-varying delay.


### 4.1. Maximizing the Estimate ofRA

An underlying issue on Problem2 is to maximize the size of the estimate of the region of attraction RA, i.e., to maximize the size of LV. Some methods to this end are reported, for instance, in [13, 29]. One possibility is to inflate the set (26)EH=x-k∈Rnτ-+1:0<x-kTHx-k≤1,that is included in LV. Such an inclusion is ensured by (27)HIIWi,τ≥0,τ=0,…,τ-,i=1,…,N.Thus, a convex optimization procedure can be formulated as (28)mintraceHsujectto15,16,and27.If the optimization procedure (28) has a feasible solution, EH is an approximation for LV⊂RA. Naturally a more precise determination is performed through (14) given in Lemma 5.Remark 11.
The numerical burden associated with the convex optimization procedure (28) depends on the number of scalar variables V and the number of LMI rows R. In our approach we have V=0.5(τ-+1)n[(τ-+1)n+1](N+τ-+1)+2(τ-+1)n[((τ-+1)n)+p]+p and R=(τ-+1)n+p+1. Although these numbers may be larger than some other approaches, the achieved results are much less conservative as it is clear in the next section.


## 4.1. Maximizing the Estimate ofRA

An underlying issue on Problem2 is to maximize the size of the estimate of the region of attraction RA, i.e., to maximize the size of LV. Some methods to this end are reported, for instance, in [13, 29]. One possibility is to inflate the set (26)EH=x-k∈Rnτ-+1:0<x-kTHx-k≤1,that is included in LV. Such an inclusion is ensured by (27)HIIWi,τ≥0,τ=0,…,τ-,i=1,…,N.Thus, a convex optimization procedure can be formulated as (28)mintraceHsujectto15,16,and27.If the optimization procedure (28) has a feasible solution, EH is an approximation for LV⊂RA. Naturally a more precise determination is performed through (14) given in Lemma 5.Remark 11.
The numerical burden associated with the convex optimization procedure (28) depends on the number of scalar variables V and the number of LMI rows R. In our approach we have V=0.5(τ-+1)n[(τ-+1)n+1](N+τ-+1)+2(τ-+1)n[((τ-+1)n)+p]+p and R=(τ-+1)n+p+1. Although these numbers may be larger than some other approaches, the achieved results are much less conservative as it is clear in the next section.


## 5. Numerical Examples

All calculations performed in this section were done with Python language and the library PICOS [30] with the solver SDPA [31] that has, in general, solved the optimization procedures with one-twentieths of the time spent by Matlab under the same computer, an Intel i53470@3.2GHz with 4Gb of RAM and running Ubuntu 16.04. In what follows, the average solving time among simple and more complex cases was 9s, including presolver processing. We present two examples. In the first one, a scalar case is addressed allowing the reader to get some insights on this new characterization of the estimate of the region of attraction and to observe the small conservatism of our approach with respect to two works in the literature. In the second example, we compare our approach with other three works in the literature to illustrate its efficacy.


### 5.1. Example 1

We study the scalar case withτ-=2 to provide a better understanding of the proposed methods. Thus, consider system (1) with A=An(1+δ), Ad=Adn(1+δ) and B=Bn(1+δ), and |δ|≤0.3 an uncertain parameter. The nominal parameters are An=2, Adn=-0.1, and Bn=1, and the symmetric saturation limit v-=0.7. Therefore, we have a polytope with two vertices computed at the extreme values of δ.Firstly, we have designed a controller for both the nominal system and the uncertain system, usingΔτmax=1 through the optimization procedure (28). The following gains were obtained: Kn=-1.01640.03590.0284 for the nominal case and Kr=-1.27200.04190.0336 for the uncertain case. Although very close from each other, these gains lead to different estimates of LV⊂RA. In Figure 1, we present the intersection of such regions with the planes xk×xk-1, xk×xk-2, and xk-1×xk-2. The gray dot marks refer to the cut on the estimate LV of the nominal system and the blue dot marks refer to the cut of the estimate LV related to the uncertain case. It is clear that the uncertainty reduces the size to the safe initial conditions. Through a spacial grid, the volumes of the estimates of safe initial conditions are for both the nominal case and the uncertain one. We verified that the estimate of region of attraction for the nominal is about 95% greater than that we got for the uncertain system. To illustrate the low conservatism of our approach, we have picked two possible systems and presented their phase portrait: one corresponding to the nominal system and the other taken inside the polytope, with α1=2/9, thus Υ(2/97/9T)=2/9A1Ad1B1+7/9A2Ad2B2. The projections of parametrized response for four different initial conditions are presented with the ellipsoidal cuts shown in Figure 1. The lines in cyan and in yellow correspond to the nominal system: one starting in the estimate of the region of attraction of the nominal system, {x0,x-1,x-2}={0.25,-2.1,-3.6}, gray dots, trajectory in cyan line; and another starting close, but outside it, {x0,x-1,x-2}={0.435,-3.2,-4}, trajectory in yellow line. Note that the initial sequence taken outside the estimated region of attraction yields a sequence that does not converge to the origin, although its trajectory passes very close to the estimated region. With regard to the uncertain system, we have taken two initial sequences: one starting in the estimate region of attraction of the uncertain system close to its boundary, {x0,x-1,x-2}={-0.554,0,0}, blue dots region, trajectory in red line; and another starting outside the estimated region of the uncertain system, but inside the one for the nominal system, {x0,x-1,x-2}={-0.621,0.0355,0.012}, gray dot, trajectory in black line. It is clear that our approach achieved a very nice estimate of RA. As expected trajectories in red and cyan lines go to the origin. On the other hand, trajectories in black and yellow lines, even starting very close to their respective safe initial conditions, diverge from the origin. Moreover, if one considers the approach in [15, 16] the bound on the estimation is, in the best case, limited by the black divergent trajectory yielding a small radius for LV1 (less than 0.62) and thus is much more conservative than ours which is about 8.7 times bigger than the optimal corresponding sphere.Figure 1
Intersections of the estimates of the region of attraction for the nominal system (gray mark) and for the uncertain system (blue dots) with the planesxk×xk-1, xk×xk-2, and xk-1×xk-2.In what follows, we evaluate the size ofLV for the uncertain system under different values of τ- and Δτmax, through the optimization procedure (28). For each value of τ-∈{3,…,7}, we have varied Δτmax from 0 up to τ--1. We have compared the percentage of augmentation in the volume of the estimate LV with respect to the case of maximum delay variation: (29)R%=VolLVΔτmaxVolLVΔτmax=τ--1-1×100%,where Vol(LV•) is volume of the estimate of the region of attraction obtained under the condition •. This volume has been verified through a regular spacial grid. The results are shown in Figure 2, where the ordinate axis is the lnR% (thus, 8 in such a figure is about 2120%). It is clear that the smaller the Δτmax, the bigger the size of the estimate of the region of attraction. Also note the case where Δτmax=0, i.e., the delay is uncertain but time-invariant. In this case the estimates of the region of attraction are much bigger than the case typically handled in the literature, where Δτmax=τ--1. This clearly shows the relevance of our approach in considering the bounds on the delay variation to get a less conservative estimate of the region of attraction.Figure 2
Percentage augmentation of the estimate of the region of attraction (R%, see the text) as a function of the maximum time variation of the delay, for τ-=3 (red square), τ-=4 (blue cross), τ-=5 (green circle), τ-=6 (purple star), and τ-=7 (▼).


### 5.2. Example 2

Consider a continuous-time system treaded in Example 2 of [23] that has been discretized here with sample time Ts=0.24s and we got(30)A=1+Ts0.6Ts-Ts1+0.8Ts,Ad=Ts0.5Ts-0.5Ts0,B=TsTs.Initially assume the delay of the continuous-time system as 1.2s which leads to τ-=5 samples. In the sequel, we describe two numerical experiments performed to compare our approach with others found in the literature.Firstly, we look for robust gains to provide local stability and that maximize the estimate of the region of attraction,RA. Because of discretization, the variation of one sample is already larger than the bounded considered in continuous-time case by [24] and thus, its conditions cannot be applied even for Δτ=1. The conditions found in [23] are applied only for time-invariant delay and those in [25] can handle the continuous-time version of the time-varying delay. Note that the conditions in [23–25] apply only to precisely known systems and yield a region of attraction described by a ball in Rn with radius r. The respective computed values are presented in Table 1.Table 1
Comparison of the radius of the ball inscribed in the estimate region of attraction.
Delay [23, Corollary 1] [24, Theorem 1] [25, Corollary 1] [25, Theorem 1] Opt. Proc. (28) Δτ=0 0.0806 0.0742 0.0772 — 0.1694  Δτ=3 — — — Infeasible 0.1289  Δτ=5 — — — Infeasible 0.1279On the other hand, our approach tracts time-varying delay with bounded time variation and the achieved results are also shown in the last column of the same table forΔτ∈{0,3,5}. Remembering that Δτ=0 implies uncertain and time-invariant delay. Following Remark 10, the achieved radius value in our approach corresponds to the biggest ball inscribed in estimated RA computed by the optimization procedure (28) and taking the smallest axis length among matrices W1,j-1,j=0,…,τ-. It is worth mentioning that even for the worst case of time-delay variation, Δτ=5, the obtained radius by our approach is more than 58.6% greater than the other approaches, at least.As we have mentioned in Remark10, the real estimated region given by our approach is much larger than this ball. For instance, with Δτ=3 we can found the sequence: (31)φ5,0=0.02480.2274,0.07240.0026,0.0490.0098,0.01360.0039,-0.0483-0.0129,-0.3795-0.1603which has φ5,0=0.2287 that is about 77% larger than the estimated norm achieved by [23], in the time-invariant delay case.The second numerical experiment performed consists in considering the uncertain time sampling interval which leads to an uncertain model representation. To this end we suppose that the sampling time belongs to the interval0.2s≤Ts≤0.3s. This leads system (1) with (30) to be represented in a two-vertex polytope. Furthermore, the delay in continuous-time becomes represented in discrete-time version by 4≤τ≤6 samples. We have performed the optimization procedure (28) for Δτ=0 and got an estimate of the RA where the bigger ball inscribed in it has a radius of 0.1870. This illustrates that, even with an uncertain sampling time, our approach yields better results than those achieved with [23–25] for precisely known system with time-invariant delay (see Table 1). As we mentioned before, our approach can find sequences such as (32)φ6,0=0.13120.3585,0.18230.0613,0.11390.0383,0.02460.0082,-0.0987-0.039,-0.1731-0.0615,-0.3618-0.1317which has φ6,0=0.3817. Note that despite the uncertainty in the matrices of the system, the delay interval is shorter and, as a consequence, the estimate for RA is even greater than the previous (precisely known) case.


## 5.1. Example 1

We study the scalar case withτ-=2 to provide a better understanding of the proposed methods. Thus, consider system (1) with A=An(1+δ), Ad=Adn(1+δ) and B=Bn(1+δ), and |δ|≤0.3 an uncertain parameter. The nominal parameters are An=2, Adn=-0.1, and Bn=1, and the symmetric saturation limit v-=0.7. Therefore, we have a polytope with two vertices computed at the extreme values of δ.Firstly, we have designed a controller for both the nominal system and the uncertain system, usingΔτmax=1 through the optimization procedure (28). The following gains were obtained: Kn=-1.01640.03590.0284 for the nominal case and Kr=-1.27200.04190.0336 for the uncertain case. Although very close from each other, these gains lead to different estimates of LV⊂RA. In Figure 1, we present the intersection of such regions with the planes xk×xk-1, xk×xk-2, and xk-1×xk-2. The gray dot marks refer to the cut on the estimate LV of the nominal system and the blue dot marks refer to the cut of the estimate LV related to the uncertain case. It is clear that the uncertainty reduces the size to the safe initial conditions. Through a spacial grid, the volumes of the estimates of safe initial conditions are for both the nominal case and the uncertain one. We verified that the estimate of region of attraction for the nominal is about 95% greater than that we got for the uncertain system. To illustrate the low conservatism of our approach, we have picked two possible systems and presented their phase portrait: one corresponding to the nominal system and the other taken inside the polytope, with α1=2/9, thus Υ(2/97/9T)=2/9A1Ad1B1+7/9A2Ad2B2. The projections of parametrized response for four different initial conditions are presented with the ellipsoidal cuts shown in Figure 1. The lines in cyan and in yellow correspond to the nominal system: one starting in the estimate of the region of attraction of the nominal system, {x0,x-1,x-2}={0.25,-2.1,-3.6}, gray dots, trajectory in cyan line; and another starting close, but outside it, {x0,x-1,x-2}={0.435,-3.2,-4}, trajectory in yellow line. Note that the initial sequence taken outside the estimated region of attraction yields a sequence that does not converge to the origin, although its trajectory passes very close to the estimated region. With regard to the uncertain system, we have taken two initial sequences: one starting in the estimate region of attraction of the uncertain system close to its boundary, {x0,x-1,x-2}={-0.554,0,0}, blue dots region, trajectory in red line; and another starting outside the estimated region of the uncertain system, but inside the one for the nominal system, {x0,x-1,x-2}={-0.621,0.0355,0.012}, gray dot, trajectory in black line. It is clear that our approach achieved a very nice estimate of RA. As expected trajectories in red and cyan lines go to the origin. On the other hand, trajectories in black and yellow lines, even starting very close to their respective safe initial conditions, diverge from the origin. Moreover, if one considers the approach in [15, 16] the bound on the estimation is, in the best case, limited by the black divergent trajectory yielding a small radius for LV1 (less than 0.62) and thus is much more conservative than ours which is about 8.7 times bigger than the optimal corresponding sphere.Figure 1
Intersections of the estimates of the region of attraction for the nominal system (gray mark) and for the uncertain system (blue dots) with the planesxk×xk-1, xk×xk-2, and xk-1×xk-2.In what follows, we evaluate the size ofLV for the uncertain system under different values of τ- and Δτmax, through the optimization procedure (28). For each value of τ-∈{3,…,7}, we have varied Δτmax from 0 up to τ--1. We have compared the percentage of augmentation in the volume of the estimate LV with respect to the case of maximum delay variation: (29)R%=VolLVΔτmaxVolLVΔτmax=τ--1-1×100%,where Vol(LV•) is volume of the estimate of the region of attraction obtained under the condition •. This volume has been verified through a regular spacial grid. The results are shown in Figure 2, where the ordinate axis is the lnR% (thus, 8 in such a figure is about 2120%). It is clear that the smaller the Δτmax, the bigger the size of the estimate of the region of attraction. Also note the case where Δτmax=0, i.e., the delay is uncertain but time-invariant. In this case the estimates of the region of attraction are much bigger than the case typically handled in the literature, where Δτmax=τ--1. This clearly shows the relevance of our approach in considering the bounds on the delay variation to get a less conservative estimate of the region of attraction.Figure 2
Percentage augmentation of the estimate of the region of attraction (R%, see the text) as a function of the maximum time variation of the delay, for τ-=3 (red square), τ-=4 (blue cross), τ-=5 (green circle), τ-=6 (purple star), and τ-=7 (▼).


## 5.2. Example 2

Consider a continuous-time system treaded in Example 2 of [23] that has been discretized here with sample time Ts=0.24s and we got(30)A=1+Ts0.6Ts-Ts1+0.8Ts,Ad=Ts0.5Ts-0.5Ts0,B=TsTs.Initially assume the delay of the continuous-time system as 1.2s which leads to τ-=5 samples. In the sequel, we describe two numerical experiments performed to compare our approach with others found in the literature.Firstly, we look for robust gains to provide local stability and that maximize the estimate of the region of attraction,RA. Because of discretization, the variation of one sample is already larger than the bounded considered in continuous-time case by [24] and thus, its conditions cannot be applied even for Δτ=1. The conditions found in [23] are applied only for time-invariant delay and those in [25] can handle the continuous-time version of the time-varying delay. Note that the conditions in [23–25] apply only to precisely known systems and yield a region of attraction described by a ball in Rn with radius r. The respective computed values are presented in Table 1.Table 1
Comparison of the radius of the ball inscribed in the estimate region of attraction.
Delay [23, Corollary 1] [24, Theorem 1] [25, Corollary 1] [25, Theorem 1] Opt. Proc. (28) Δτ=0 0.0806 0.0742 0.0772 — 0.1694  Δτ=3 — — — Infeasible 0.1289  Δτ=5 — — — Infeasible 0.1279On the other hand, our approach tracts time-varying delay with bounded time variation and the achieved results are also shown in the last column of the same table forΔτ∈{0,3,5}. Remembering that Δτ=0 implies uncertain and time-invariant delay. Following Remark 10, the achieved radius value in our approach corresponds to the biggest ball inscribed in estimated RA computed by the optimization procedure (28) and taking the smallest axis length among matrices W1,j-1,j=0,…,τ-. It is worth mentioning that even for the worst case of time-delay variation, Δτ=5, the obtained radius by our approach is more than 58.6% greater than the other approaches, at least.As we have mentioned in Remark10, the real estimated region given by our approach is much larger than this ball. For instance, with Δτ=3 we can found the sequence: (31)φ5,0=0.02480.2274,0.07240.0026,0.0490.0098,0.01360.0039,-0.0483-0.0129,-0.3795-0.1603which has φ5,0=0.2287 that is about 77% larger than the estimated norm achieved by [23], in the time-invariant delay case.The second numerical experiment performed consists in considering the uncertain time sampling interval which leads to an uncertain model representation. To this end we suppose that the sampling time belongs to the interval0.2s≤Ts≤0.3s. This leads system (1) with (30) to be represented in a two-vertex polytope. Furthermore, the delay in continuous-time becomes represented in discrete-time version by 4≤τ≤6 samples. We have performed the optimization procedure (28) for Δτ=0 and got an estimate of the RA where the bigger ball inscribed in it has a radius of 0.1870. This illustrates that, even with an uncertain sampling time, our approach yields better results than those achieved with [23–25] for precisely known system with time-invariant delay (see Table 1). As we mentioned before, our approach can find sequences such as (32)φ6,0=0.13120.3585,0.18230.0613,0.11390.0383,0.02460.0082,-0.0987-0.039,-0.1731-0.0615,-0.3618-0.1317which has φ6,0=0.3817. Note that despite the uncertainty in the matrices of the system, the delay interval is shorter and, as a consequence, the estimate for RA is even greater than the previous (precisely known) case.


## 6. Conclusions

We have presented new convex conditions formulated in terms of LMIs to design state feedback control law that robustly locally stabilizes uncertain discrete-time systems with time-varying delay actioned by saturating actuators and under bounded time variation delay. All matrices of the systems are supposed to be affected by polytopic uncertainty. It has been shown that taking into account a constraint on the time variation of the delay can lead to less conservative analysis and synthesis conditions. We have given a new characterization of the region of attraction and the proposed conditions allow computing much less conservative estimates of such a region. Although the proposed conditions are delay dependent—which is recognized as a nice property in the literature—the designed control law is delay independent, which facilitates the real time implementation. We have compared our approach with others found in the literature through two numerical examples, showing the efficiency of our proposal.



---

*Source: 5013056-2018-07-29.xml*