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

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.


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][4][5][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][9][10][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,   , can jump from its lower bound, usually   =  = 1 at one instant, to its maximum value, , in the next instant, i.e.,  +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][15][16][17][18][19].

Mathematical Problems in Engineering
The description of the saturating function can be done by a polytopic based representation [15][16][17] or through a sector condition approach as in [14,[17][18][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 R  .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 R  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 R  that contains the current state vector and a ball in R  , 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  ×  is R × .A block diagonal matrix with the matrices  1 and  2 in the diagonal blocks is noted as block-diag { 1 ,  2 }.The  ℎ line of the matrix (or vector)  is denoted by  () .If  is a sequence of vectors, then its  ℎ element is noted by []  .The sequence  , ∈   , with

Problem Statement
Consider the class of uncertain discrete-time systems with state delay and saturating actuators described by  +1 =  ()   +   ()  −  +  () sat (V  ) , (1) where   = () ∈ R  is the state vector, V  is the linear control signal, sat(V  ) is a decentralized vectorial saturation function: , otherwise, with V () ,  = 1, . . ., , being the symmetrical bound on the values achieved by the  ℎ actuator.The delay in the state is time-varying and represented by   verifying 0 ≤   ≤ . ( The matrices of the system are denoted by Υ() ≡ [()   () ()] ∈ R ×(2+) = [   ]() with the time-invariant uncertain parameter  ∈ Ξ where the matrices Υ() belonging to the polytope S with  known vertices and Υ  ≡ [   ,   ] = [   ]  .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 V  , as discussed later in this work.A contribution of this approach is to consider the time variation of the delay   bounded as follows: Therefore, the transitions of   between two successive instants  and +1 are such that if   = , then 0 ≤ −Δ max ≤  +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   = : where  + =  +1 .Then, for  = 5,   = 1, and Δ max = 2, we have C() = {0, 1, 2, 3}, which is the set of all possible values that can be assumed by  +1 .
The following delay independent control law is investigated to robust stabilize system (1): where . Two particular structures of interest for (7) are (i)  = 0,  = 1, . . ., , where only the current state is used in the feedback, and (ii)  = 0,  = 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   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, V  =  0   +    −  .Observe that the notion of stability in state delayed systems requires that all elements of the sequence  , goes to the origin as  → ∞ [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  , do not go to zero as  → ∞.Therefore, we adapt the definition of region of attraction, R A , from [13] as follows.
Definition 1 (region of attraction, R A ).The region of attraction R A of the origin for system (1) is the set of all sequences  ,0 whose respective trajectories converge asymptotically to the origin, i.e., lim →∞  , = {0, . . ., 0}.
Therefore, under saturated actuators it is necessary to determine the set R A or a subset of it, L V ⊂ R A , as large as possible, such that the local stability can be ensured.As it is widely known in the literature, the determination of R A is very difficult task [13, p. 15] and the determination of an estimate L V ⊂ R A 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 matrix K and characterize an estimate L V of the region of attraction R A such that the trajectories of system (1)-( 5) under control law (7) and initial conditions given by all sequences  ,0 with elements taken from L V ∈ R A converge to the origin, for all  ∈ Ξ.
Note that it is embedded in Problem 2 the issue of consider a bound on the time variation of the delay, which is also a contribution of this work.

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   is the switching function: where Γ ,,  =   () if  =   and Γ ,,  = 0 × , otherwise, , and   given after (7).Using the control law (7) into system (8), we can get  +1 = A(,   )  + B()sat(K  ).With the aid of the dead zone function, Ψ(V  ) = V  − sat(V  ), we sum and subtract B()K  in the right side of (8) to get  +1 = (A(,   ) + B()K)  − B()Ψ(K  ).To simplify the notation, we assume Ã(,   ) = A(,   )+B()K and the closed-loop system can be rewritten as Lemma 3 (generalized sector condition [21]).Assume V  given by ( 7), V = −V > 0 and a matrix  ∈ R ×(+1) , such that The stability here is investigated by means of Lyapunov's Theorem [1].A candidate function ( , , ,   ) > 0 is called a Lyapunov function if there exist positive constants   , for In this work, we use the following augmented Lyapunov candidate function: where 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 (,   ).The estimate of R A is performed here through the level set L V associated with a contractive set determined by the augmented Lyapunov function.Definition 4 (level set, L V ).Suppose that the function given by ( 12)-( 13) is a Lyapunov function for system (8) with ( 3)- (5).The level set L V associated with ( 12)-( 13) is defined as the intersection of the ellipsoidal sets E((, , for all  ∈ Ξ and   ∈ {0, . . ., }.
It is worth saying that the computation of the set L V can be done as indicated in the next lemma, which is adapted from [27].
The level set L V from Lemma 5 is positively invariant and, thus, a trajectory started inside it goes to the origin.Therefore, L V is an estimate of R A , providing a set of safe initial conditions for the elements of the sequence  ,0 : each point in the space of L V corresponds to a complete sequence of initial conditions.Such a set characterization differs from those used in [14][15][16][17][18][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.
Remark 7. Theorem 6 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  to another  + 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 Theorem 6 just by suppressing  and, as a consequence, also the index  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 (,   ).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 matrix K in (7) can be recovered from ( 15)-( 18) by imposing some special structures on matrices  and .
On the other hand, to get the gain K = [ 0 0 ⋅ ⋅ ⋅ 0   ], it is necessary to constrain the structure of  and  as An interesting issue is that, in all cases, no structure is imposed to the Lyapunov candidate matrices  , ,  = 1, . . .,  and  = 0, . . ., , allowing a larger set of search for gain K.

Maximizing the Estimate of R
A .An underlying issue on Problem 2 is to maximize the size of the estimate of the region of attraction R A , i.e., to maximize the size of L V .Some methods to this end are reported, for instance, in [13,29].One possibility is to inflate the set that is included in L  .Such an inclusion is ensured by Thus, a convex optimization procedure can be formulated as min trace () suject to (15) , (16) , and ( 27) .
If the optimization procedure (28) has a feasible solution, E  is an approximation for L V ⊂ R A .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)[( + 1) + 1]( +  + 1) + 2( + 1)[(( + 1)) + ] +  and R = ( + 1) +  + 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.

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 5 3470@3.2GHzwith 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 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 planes   ×  −1 ,   ×  −2 , and 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.
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: K  = [−1.01640.0359 0.0284] for the nominal case and K  = [−1.27200.0419 0.0336] for the uncertain case.Although very close from each other, these gains lead to different estimates of L V ⊂ R A .In Figure 1, we present the intersection of such regions with the planes   ×  −1 ,   ×  −2 , and  −1 ×  −2 .The gray dot marks refer to the cut on the estimate L V of the nominal system and the blue dot marks refer to the cut of the estimate L V 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/9 7/9]  ) = (2/9) [ 1  1  1 ] + (7/9) [ 2  2  2 ].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, { 0 ,  −1 ,  −2 } = {0.25,−2.1, −3.6}, gray dots, trajectory in cyan line; and another starting close, but outside it, { 0 ,  −1 ,  −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, { 0 ,  −1 ,  −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, { 0 ,  −1 ,  −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 R A .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 L V1 (less than 0.62) and thus is much more conservative than ours which is about 8.7 times bigger than the optimal corresponding sphere.
In what follows, we evaluate the size of L V 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 L V with respect to the case of maximum delay variation: where Vol(L V• ) 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 ln  % (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.
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, R A .Because of discretization, the variation of one sample is already larger than the bounded considered in continuoustime 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][24][25] apply only to precisely known systems and yield a region of attraction described by a ball in R  with radius .The respective computed values are presented in Table 1.
On 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 R A computed by the optimization procedure (28) and taking the smallest axis length among matrices  −1 1, ,  = 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 Remark 10, the real estimated region given by our approach is much larger than this ball.For instance, with Δ = 3 we can found the sequence: which 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 interval 0.2 ≤   ≤ 0.3.This leads system (1) with (30) to be represented in a twovertex 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 R A 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][24][25] for precisely known system with time-invariant delay (see Table 1).As we mentioned before, our approach can find sequences such as which 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 R A is even greater than the previous (precisely known) case.

Conclusions
We have presented new convex conditions formulated in terms of LMIs to design state feedback control law that Mathematical Problems in Engineering 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.

Table 1 :
Comparison of the radius of the ball inscribed in the estimate region of attraction.Percentage augmentation of the estimate of the region of attraction ( % , 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 ().
[23]Example 2.Consider a continuous-time system treaded in Example 2 of[23]that has been discretized here with sample time   = 0.24 and we got