STABILITY OF DIFFERENTIAL INCLUSIONS : A COMPUTATIONAL APPROACH

We present a technique for analysis of asymptotic stability for a class of differential inclusions. This technique is based on the Lyapunov-type theorems. The construction of the Lyapunov functions for differential inclusions is reduced to an auxiliary problem of mathematical programming, namely, to the problem of searching saddle points of a suitable function. The computational approach to the auxiliary problem contains a gradienttype algorithm for saddle-point problems. We also extend our main results to systems described by difference inclusions. The obtained numerical schemes are applied to some illustrative examples.


Motivation and overview
The method of Lyapunov functions is one of the most efficient methods for analyzing the stability of nonlinear dynamic systems (see, e.g., [12,22,41]).As in the case of ordinary differential equations, the main difficulty encountered in using this method for differential inclusions lies in constructing the Lyapunov function with requisite properties.The aim of our paper is to obtain an algorithm for choosing the Lyapunov function for a class of differential inclusions.Recall that differential inclusions are usual mathematical representations of control systems with ordinary differential equations (see, e.g., [3,10,28]).Let us consider the following initial-value problem [3,7,13,21]: x t 0 = 0, where x ∈ R n and Ꮾ q , q ∈ N, is the convex hull of the given real (n × n)-matrices A 1 ,..., A q , that is, For example, the set of the initial-value problems for linear nonstationary differential equations ẋ = A(t)x, x(0) = 0, A(t) = a i j (t) n i, j=1 , α i j ≤ a i j (t) ≤ β i j , i, j = 1,...,n, (1.3) where α i j , β i j are constants, can be reduced to problem (1.1).
In parallel with (1.1), we examine a more general type of the initial-value problem for differential inclusions ẋ ∈ F(x), x t 0 = 0, where Ꮾ is a compactum (in general, nonconvex) in the n 2 -dimensional space of real (n × n)-matrices A.
Recall that an absolutely continuous vector function x(•) satisfying the condition ẋ(t) ∈ F q (x(t)) (or the condition ẋ(t) ∈ F(x(t))) almost everywhere on a considered interval of time [t 0 ,t] is called a solution of the differential inclusion in (1.1) (a solution of the differential inclusion in (1.4)).Note that any solution x(•) of above inclusions can be continued on the whole semi-infinite axis [t 0 ,∞) (see [3]).The solution x(•) ≡ 0 of the initial-value problem (1.1) or (1.4) is called the zero solution (or trivial solution, see, e.g., [7]).The problem of asymptotic stability of the zero solution for systems of the type (1.1) and (1.4), to which many practically important control systems can be reduced, consists of choosing a Lyapunov function V : R n → R + .This function satisfies the following conditions (see [7,13]): (i) The Lyapunov functions V (•) for systems (1.1) and (1.4) can be chosen from the class of convex functions [27].For the differential inclusions in (1.1) and in (1.4), the role of usual derivatives is played by the functions where ∇V (x)y = lim h→0+ h −1 (V (x + hy) − V (x)) is the Gateaux derivative (see [13]).The bracket (•,•) denotes the inner product in R n .We call the function W q (•) or W(•) the derivative of the function V (•) along solutions of system (1.1) or (1.4), respectively.
In this paper, we also study the question of asymptotic stability of the zero solution for systems of difference inclusions.It is common knowledge that a difference inclusion can be obtained as discretization of the corresponding differential inclusion, and provides a basis for the numerical treatment of the given differential inclusion.We refer to [8,9,28,37] for details.Note that our paper is concerned only with a class of difference inclusions determined by the multifunction F q (•).
The problem of asymptotic stability of the zero solution for system (1.1) is closely related to the familiar problem of absolute stabilization of a nonstationary control system.Alternatively, the above problem can also be treated as a variant of the so-called robust stabilization problem.The concept of absolute stability of feedback systems dates back to Lourie and Postnikov [23].Consider a simple model of a plant with nonlinear gain u given in terms of the output σ as follows:  [40] and by Kalman [18].We refer to [30] for a historical survey of the problem.Some applications of the frequency techniques to the stability of delay systems are discussed in [15].For the extensions to distributed systems, see [6,39].The differential inclusion in (1.1) is equivalent to the following set of linear control systems (see [35] for details): where the control functions u ν (•), ν = 1,..., q, are Lebesque measurable on each finite segment of the t-axis and Ꮽ is an (n × n) matrix.Thus the problem of asymptotic stability of the zero solution for system (1.1) can be interpreted as the corresponding absolute stabilization problem for the presented set of linear control systems.Robust stability of time-varying systems with uncertainties has received much attention since the familiar paper of Kharitonov [20].A great amount of works is devoted to the theoretical and practical aspects of robust stability and robust stabilization of dynamical systems; see [1,16,25,29] and the references therein.The so-called quadratic stabilization is a powerful approach to the problem of robust stabilization for uncertain linear systems (see, e.g., [5,31]).The problem is to construct a common quadratic Lyapunov function for all possible uncertainties.In this context, an uncertainty is meant in a deterministic sense: it arises as a result of approximation, imprecision, or imperfect knowledge introduced during the modeling procedure.In the present paper, the uncertain system is defined by a differential inclusion, the right-hand side of which is a known multifunction F q or F. Quadratic stability dates back to the pioneering work of Lyapunov [24] who established that the existence of a quadratic Lyapunov function is a necessary and sufficient condition for asymptotic stability of a linear system.Quadratic Lyapunov functions are often the first resort in the analysis of nonlinear systems and much work on absolute stability is based on quadratic Lyapunov functions.We refer to [17,36] for some new directions in the area of quadratic stability and quadratic stabilization.The formulated problem of asymptotic stability of the zero solution for system (1.1) can also be considered in the framework of quadratic stabilization.
The aim of this paper is to propose a new numerical approach to the described problem of asymptotic stability for systems (1.1) and (1.4) and to the corresponding robust stabilization problem.Moreover, we also investigate the question of asymptotic stability for the linear "nonconvex" system (1.4) (Theorems 2.4 and 2.5).Using these new theoretical results, one can construct the Lyapunov function for a suitable system (1.1) in place of the corresponding problem for the more general system (1.4).Note that the problem of asymptotic stability of the zero solution for system (1.4) is not from the class of quadratic stabilization problems.In the present paper, we do not consider any standard techniques from the theory of absolute stability or robust stability (e.g., Popov-type criteria, Brockett technique, stability radius, LMI-based methods, Riccati equation approach, and the like).We apply some results from the classical stability theory for differential inclusions (see, e.g., [7,12,13]).Moreover, we discuss a relaxation of the usual stability concept [7].In contrast to the quadratic stabilization, we construct polynomial Lyapunov functions for system (1.1) and for the corresponding discrete system.This approach makes it possible to reduce the problem of choosing the Lyapunov function to an auxiliary problem of mathematical programming, namely, to a saddle points problem (Theorems 3.1 and 4.1).We solve this saddle points problem by application of a gradient-type method.
The remainder of the paper is organized as follows.Section 2 contains some basic facts about the asymptotic stability.In Section 3, we propose a constructive approach to the stability problem and establish a relation between stability of systems given above and an auxiliary saddle points problem.Section 4 contains the related results for a class of difference inclusions.In Section 5, we apply a gradient-type method to the auxiliary saddle points problem and present an algorithm for constructing the Lyapunov functions.Finally, we consider two illustrative examples.

Mathematical preliminaries
First let us discuss some definitions and theoretical results in connection with the stability theory for differential inclusions (see [3,7,13]).

Vadim Azhmyakov 5
The zero solution x(•) ≡ 0 is said to be weakly asymptotically stable if the above definition holds with "for each" replaced by "for some."For the classical Lyapunov-type stability theorems, see [7,13].Let us present two special stability theorems for the given system (1.4) [27].
Theorem 2.2.For the zero solution x(•) ≡ 0 of problem (1.4) to be asymptotically stable, it is necessary and sufficient that there exists a strictly convex, homogeneous (of second order) Lyapunov function V (•) of a quasiquadratic form, namely, whose derivative along solutions of system (1.4) satisfies the inequality Note that Theorem 2.2 is a variant of the converse theorems of Lyapunov's direct method for differential inclusions [22].In parallel with the quasiquadratic Lyapunov functions, one can also consider the Lyapunov functions from some other classes of functions, for one, from the class of homogeneous forms [27].

is asymptotically stable if and only if there exists a Lyapunov function in the class of homogeneous forms of order
where l i ∈ R n , i = 1,...,m are constant vectors with such that for its derivative along solutions of system (1.4), the inequality is satisfied with an integer p ≥ 1.
Since (1.1) is a special case of (1.4), Theorems 2.2 and 2.3 provide the necessary and sufficient conditions of the asymptotic stability of the zero solution x ≡ 0 for systems (1.1) and (1.4).Evidently, the parameters determining the class of Lyapunov functions V m,p (•,•) are the components of the vectors l i , i = 1,...,m, and the numbers m and p.As is shown in [27], the Lyapunov functions in Theorems 2.2 and 2.3 can be chosen from the class of convex (with respect to x) functions.
The next result establishes a useful link between problem (1.4) and an approximate initial-value problem.
Theorem 2.4.The zero solution x(•) ≡ 0 of problem (1.4) is asymptotically stable if and only if the zero solution of the problem x t 0 = 0, is asymptotically stable.
Proof.Sufficiency of the assertion follows from the inclusion F(x) ⊆ F c (x). Necessity follows from the equivalence of the closure of the solution set of problem (1.4) and the solutions set of problem (2.7) (see [3,13]).
Since a convex compactum co(Ꮾ) can be approximated by convex polyhedrons Ꮾ q , we have the following result.
Theorem 2.5.For the zero solution x(•) ≡ 0 of problem (1.4) to be asymptotically stable, it is necessary and sufficient that there exists a number q ≥ 1 and an initial-value problem (1.1) whose zero solution x(•) ≡ 0 is asymptotically stable and (2.8) Proof.Let q ∈ N be a number such that F c (x) ⊂ F q (x).We assume that the zero solution x(•) ≡ 0 of (1.1) is asymptotically stable.The asymptotic stability of the zero solution of the initial-value problem (2.7) is an immediate consequence of (2.8).By Theorem 2.4, the zero solution of problem (1.4) is asymptotically stable too.Now we assume that the zero solution x(•) ≡ 0 of problem (1.4) is asymptotically stable.It follows from the results of [3,12] that there exists some > 0 such that the zero solution of the initial-value problem ẋ ∈ F (x), x t 0 = 0, where Ꮾ is a compactum and Ꮾ ⊂ Ꮾ , is asymptotically stable.This implies that the zero solution of the problem is asymptotically stable (see Theorem 2.4).Then there exists a number q ∈ N such that co(Ꮾ) ⊂ co A 1 ,...,A q , A s ∈ co Ꮾ , s = 1,..., q. (2.11) We have co A 1 ,...,A q ⊂ co Ꮾ , A s ∈ co Ꮾ , s = 1,..., q. (2.12) Vadim Azhmyakov 7 This means that there is a number q ≥ 1 such that the solution x(•) ≡ 0 of (1.1) is asymptotically stable and condition (2.8) holds.
In a similar way, one can formulate the corresponding results for week asymptotically stable systems.
Theorem 2.6 can be proved in the same way as the above-mentioned Theorems 2.2 and 2.3.
We will touch briefly on the problem of stabilization for the control system where f : R × R n × R m , x is the state variable, and U is a control region.The asymptotic stability in this case means that if x(t 0 ) < δ, then the trajectory x(•) satisfies the inequality x(t) < (t 0 ≤ t < ∞) for every admissible control function u(•),u(t) ∈ U.
The weak asymptotic stability means the same inequality holds only for some admissible controls.

Lyapunov functions and saddle points problem
The Lyapunov function from Theorem 2.3 has the polynomial form or equivalently where ψ r (x), r = 1,...,N(p), are standard monomials of degree 2p ≥ n, p) are coefficients of monomials and N(p) = C 2p n+2p−1 is the number of monomials.The derivatives W p,q (z,x) and W − p,q (z,x) of the function V p (•,•) along solutions 8 Stability of differential inclusions of the given system (1.1) have the form where λ := (λ 1 ,...,λ q ) T , Θ := {λ ∈ R q : q j=1 λ j = 1, λ j ≥ 0}, and is an (N(p) × n)-matrix.We will use the following notation: x := 0. Now we consider the bounded region and investigate the asymptotic stability of system (1.1) for x ∈ G x .Note that under the above-presented condition (2.8) Theorems 2.4 and 2.5 make it possible to reduce the question of asymptotic stability of system (1.4) to the same question for system (1.1).Therefore, we give the constructive characterization of condition (2.6) in Theorem 2.3 only for the initial-value problem (1.1).Theorem 3.1.Let V p (•,•) be the Lyapunov function from the class of homogeneous forms of order 2p, p ∈ N, for system (1.1).The inequality has a solution z ∈ G z if and only if the inequalities hold for all z ∈ G z , x ∈ G x , and x = 0.
Proof.Assume that the inequality W p,q ( z,x) ≤ −γ x 2p holds for a vector z ∈ G z .Then we have and the necessity of the assertion is proved.Now we assume that (3.7) is true.Clearly, It is matter of direct verification to prove that the set F q (x), x ∈ G x , is a bounded convex closed (a convex compact) subset of R n .Moreover, we deal with continuous functions.
In other words, the pair ( z, x) ∈ G z × G x is a saddle point of the function W p,q (•,•).After this investigation, we come back to the question leading to the construction of the Lyapunov function V p (•,•) for the given system (1.1) (or system (1.4)).This problem is now reduced to an auxiliary problem of mathematical programming, namely, to the saddle points problem described by inequalities (3.7).Note that in practice we define The constant Γ z can be determined, for instance, by a maximal machine number.We can also obtain the similar result for weakly asymptotic stable systems.
Theorem 3.2.Let V p (•,•) be the Lyapunov function from the class of homogeneous forms of order 2p, p ∈ N, for system (1.1).The inequality has a solution z ∈ G z if and only if the inequalities hold for all z ∈ G z , x ∈ G x , and x = 0.
The proof of this theorem is analogous to the proof of Theorem 3.1.

On the asymptotic stability of difference inclusions
In this section, we study the initial-value problem for difference inclusions where k = 0,1,..., and the multifunction F q (•) is determined in Section 1.The definition of asymptotic stability of the zero solution x(•) ≡ 0 of (4.1) is analogous to Definition 2.1.
In particular, as in the case of the initial-valued problem (1.1), the asymptotic stability of the zero solution of (4.1) can be established by means of the Lyapunov function V m,p (•,•) which in this case will satisfy the inequality (see, e.g., [27]) V m,p (l, y) − αV m,p (l,x) ≤ 0, ( for some 0 < α < 1.For system (4.1),we introduce the Lyapunov function V p (•,•) (see Section 3) and consider the inequality for z ∈ G z and x ∈ G x .Evidently, relation (4.3) implies the inequality "along solutions" of the discrete system (4.1),namely, the inequality The difference inclusion in (4.1) can be considered as a discretization scheme for (1.1) (see [9]).
Theorem 4.1.Let V p (•,•) be the Lyapunov function from the class of homogeneous forms of order 2p, p ∈ N, for system (4.1).The inequality (4.3) has a solution z ∈ G z if and only if the inequalities W p,q ( z,x) ≤ W p,q ( z, x) ≤ W p,q (z, x) ( are fulfilled for all z ∈ G z and x ∈ G x , x = 0. Proof.Let z ∈ G z be a solution of (4.3).Since and F q (x) = {0}, we have W p,q ( z, x) = W p,q (z, x) = 0.This implies that the saddle points inequalities (4.5) are fulfilled.Now we assume that ( z, x) ∈ G z × G x is a solution of (4.5).Since we obtain W p,q ( z,x) ≤ 0 for x ∈ G x , x = 0.
As is evident from the foregoing, the problem of constructing the Lyapunov function for system (4.1) with difference inclusions is also reduced to an auxiliary saddle points problem.

Numerical aspects
This section is devoted to the numerical treatment of the auxiliary saddle points problems given above.We solve the saddle points problem by application of the gradient-type algorithm.Note that the gradient method is one of the most popular first-order methods for solving problems of mathematical programming (see, e.g., [11,32]).Let us consider with 0 < K s < ∞.Evidently, a function u s (x) := (c s ,x) in (5.4) can be interpreted as a linear feedback control function for the control system x(0) = 0. (5.6) In this case, the system (5.4) is a closed-loop control system.One can rewrite (5.4) in the following equivalent form [4]: where It can be shown (see [4]) that (5.7) is a special case of (1.1) if where h νs = 0 or h νs = K s , ν = 1,..., q, and s = 1,...,S.We now give an example of system (5.4).
Example 5.3.We examine the same system (5.10) and construct the Lyapunov function V 2 (•,•) from the class of the fourth order forms.Note that p = 2 and N(p) = 5.For K = 6.40 we have V 2 (z,x) = x 4 1 + 0.301x 3 1 x 2 + 0.620x 2 1 x 2 2 + 0.211x 1 x 3 2 + 0.101x 4 2 . (5.14) The gradient-type method (5.1) and the limit procedure from Theorem 5.1 in Examples 5.2 and 5.3 are implemented in C. We used the "numerical recipes in C" package [34] and some author programs for this purpose.

Concluding remarks
In this paper, we establish the important role of saddle points problems in computing polynomial Lyapunov functions for a class of initial-value problems for differential and difference inclusions.We investigate a new numerical schema derived from a variant of the gradient method.This schema can also be treated as a practicable computational approach to the related problem of robust stability.The presented algorithm has usual advantages and disadvantages of a first-order numerical method.Note that the main idea of the algorithm can be combined not only with the gradient-type method but also with some other techniques for saddle points problem.It is conceivable that one can use the proposed computational approach and the corresponding results in studies of the socalled receding horizon control (see, e.g., [14,26]).Moreover, this approach can also be applied to some other classes of continuous and discrete control systems governed by differential and difference inclusions.