=Paper=
{{Paper
|id=Vol-1987/paper72
|storemode=property
|title=Semidefinite Programming Methods of Ascertaining Stability for a Nonlinear Delay System Subject to Constraints
|pdfUrl=https://ceur-ws.org/Vol-1987/paper72.pdf
|volume=Vol-1987
|authors=Natalya Sedova,Vadim I. Shmyrev,Ruslan Yu. Simanchev,Inna V. Urazova,Alexander K. Skiba,Stepan P. Sorokin,Maxim V. Staritsyn,Alexander S. Strekalovsky,Anna Tatarczak,Nikolay P. Tikhomirov,Dmitry V. Aron,Alexey A. Tret'yakov,Sergey Trofimov,Aleksey Ivanov,Yury Fettser,Tatiana S. Zarodnyuk,Alexander Yu. Gornov,Anton S. Anikin,Evgeniya A. Finkelstein,Elena S. Zasukhina,Sergey V. Zasukhin,Vitaly Zhadan,Anna V. Zykina,Olga N. Kaneva
}}
==Semidefinite Programming Methods of Ascertaining Stability for a Nonlinear Delay System Subject to Constraints==
https://ceur-ws.org/Vol-1987/paper72.pdf
Semidefinite Programming Methods of Ascertaining
Stability for a Nonlinear Delay System Subject to
Constraints
Natalya Sedova
Ulyanovsk State University
Leo Tolstoi str. 42,
432000 Ul’yanovsk, Russia.
sedovano@sv.ulsu.ru
Abstract
This paper deals with stability problem for nonlinear delay systems.
The proposed approach is based on the assumption that on a subset of
instant states the system is represented by a continuous-time Takagi–
Sugeno system with delay. The aim of the study is to present sufficient
conditions of stability for the system subject to the constraints defin-
ing the subset under consideration. Because of the constraints, another
important issue is how to estimate the set of initial points such that the
trajectories starting at those points remain to satisfy the constraints.
Using Lyapunov functions with the Razumikhin conditions, we obtain
stability conditions and an estimation of the domain of attraction in
terms of linear matrix inequalities with scalar parameters. The ob-
tained problems can be formulated by semidefinite programming and
solved numerically.
1 Introduction
This paper aims at establishing a framework to ascertain the stability and to estimate the domain of attrac-
tion (DA) of the equilibrium point (the origin) for nonlinear systems with time delay. It is assumed that
there is a representation of the system by a continuous-time Takagi–Sugeno (TS) system with delay (see e.g.
[Tanaka & Wang, 2001] and references therein). The original TS system was proposed as a fuzzy system and it
was described by fuzzy IF-THEN rules which represent local linear input-output relations of a nonlinear system.
The main feature of a TS fuzzy model is to express the local dynamics of each fuzzy implication (rule) by a linear
system model. The overall model of the system is achieved by fuzzy “blending” of the linear system models, so
it can be presented in a “non-fuzzy” form. A wide class nonlinear dynamic systems can be represented by TS
systems. In fact, it is proved that TS systems are universal approximators; some approaches, which derive a TS
model from given nonlinear dynamical models, are also known [Tanaka & Wang, 2001].
For such systems, the method of Lyapunov-Krasovskii functionals leads to stability conditions which can be
determined by solving some generalized eigenvalue minimization problems (GEVPs) or a system of linear matrix
Copyright ⃝
c by the paper’s authors. Copying permitted for private and academic purposes.
In: Yu. G. Evtushenko, M. Yu. Khachay, O. V. Khamisov, Yu. A. Kochetov, V.U. Malkova, M.A. Posypkin (eds.): Proceedings of
the OPTIMA-2017 Conference, Petrovac, Montenegro, 02-Oct-2017, published at http://ceur-ws.org
504
inequalities (LMIs). Such problems can be efficiently handled by means of convex optimization techniques
[Nesterov & Nemirovski, 1993]. From the other hand, when dealing with TS systems, what makes the problem
more challenging is that the system can be usually represented (exactly or approximately) in the TS form only
on some subset of the state space containing the origin (“modeling region”). Also, the system thus modeled may
have physical constraints precisely reflected in the states belonging to some modeling region. So the problem
arises to estimate the DA, that is the set of initial points such that the trajectories starting at those points
remain in the modeling region and tend to zero.
For a Lyapunov function which guarantees the local stability of the equilibrium of ODE, any sublevel set of the
Lyapunov function is an inner estimate of the DA if the set belongs to the region where the function is positive
definite and its derivative with respect to the system is negative definite. For delay equations, due to the infinite
dimensional property of the space of initial values, a proper definition of DA remains a fundamental problem.
Some authors use generalized concepts of DA. Such a DA is assumed to be a subset of a certain functional space.
However, it is a hard task to construct a set that estimates such DA, as well as to check in practice whether an
initial function is contained in the set. For instance, in [Leng et al., 2016] the volume of the DA is defined using
the projection of the initial state space to a finite-dimensional space by expanding the initial function along with
some basis. But the authors of [Leng et al., 2016] propose only a numerical procedure to estimate the volume via
numerical integrating the equation with random initial functions from the projection and counting the number
of the functions that finally arrive at the origin.
It is more convenient to define the so-called direct DA that is a subset of the finite dimensional space of instant
states x(t) (see, e.g., [Prasolov, 1998]). Then it would reasonable to use classical auxiliary functions rather than
Lyapunov–Krasovskii’s functionals which depend on the true states xt . Taking into account time delay, we
will use the method of Lyapunov functions with additional restrictions, namely, the Razumikhin conditions
[Razumikhin, 1956]. And, as for the ODEs, sublevel sets of Lyapunov functions will be employed as inner
estimates of the direct DA. The usage of a quadratic function for a TS system leads to stability conditions in
terms of LMIs. If the LMIs are feasible, the resulted quadratic function seems to be a global Lyapunov function
for the system. But the set of LMIs obtained depends not only on a matrix variable but also on some scalar
parameters. Therefore, the LMIs are not linear in all variables. In addition, asymptotic stability conditions
in this case stand good only within the modeling region. Thus, the problem is to obtain the largest possible
Lyapunov-based estimate of the DA of a nonlinear system with the asymptotically stable origin and subject to
the given constraints. In fact, we need to find invariant subsets of the DA that fit into the modeling region.
In this note, some linear matrix inequalities with parameters are presented which can be used to check a
delay-dependent sufficient condition for stability of the zero solution of a nonlinear delay differential system, and
to find an inner estimate of the DA of the solution. We consider some kinds of constraints that can be also
turned into LMIs. Then, semidefinite programs associated with those problems are formulated which can be
efficiently handled via available software.
2 Preliminaries
Notations used throughout the paper is fairly standard. Time delay is denoted by√r∑ (r > 0), R+ = [0, +∞), Rn
⊤ n
denotes the n-dimensional space of vectors x = (x1 , . . . , xn ) with the norm |x| = 2 n
i=1 xi , C = C([−r, 0], R )
is the Banach space with the supremum-norm ∥·∥. For a continuous function x(t) ∈ C([α−r, α+β), R ) (α ∈ R ,
n +
β > 0) an element xt ∈ C is defined for any t ∈ [α, α + β) by xt (s) = x(t + s), −r ≤ s ≤ 0, ẋ(t) stands for the
right-hand derivative.
Then, M < 0 (M ≤ 0) means that M is a real symmetric and negative definite (semidefinite) matrix. The
symbol ∗ within a matrix represents the symmetric of the matrix.
We will study systems dynamics of which are given by the following equation:
ẋ(t) = X(t, xt ), X(t, 0) ≡ 0, (1)
where X : R+ ×C → Rn . The function x(t) is said to be a solution of (1) if x(t) ∈ C([t0 −r, t0 +β], Rn ) for certain
t0 ∈ R+ , β > 0 and x(t) is an absolutely continuous function on [t0 , t0 + β] satisfying (1) almost everywhere on
[t0 , t0 + β]. The solution of (1) satisfying the given initial condition xt0 = φ is denoted by x(t; t0 , φ).
In the sequel, we assume that the system is considered on a set D ⊆ Rn and the origin is an internal point of
this set (the specification of such a set can be stipulated, e.g., by the physical meaning of state variables). For
convenience, suppose that the set D is defined by
D = {x ∈ Rn : gk (x) ≤ 0, k = 1, . . . , l} (2)
505
for some l ≥ 1. Also, the right-hand side of (1) is believed to satisfy the assumptions from [Sedova, 2011] which
ensure the existence, uniqueness, and continuability of solution for any initial function φ ∈ C([−r, 0], D). In
particular, if φ = 0 then X(t, 0) ≡ 0 implies x(t; t0 , 0) ≡ 0 for any t0 ∈ R+ , t ≥ t0 .
Let ξ(t) be a piecewise continuous vector function whose values at the current time t depend on t and on xt
and belong to the set Dξ ⊂ Rs whenever t ∈ R+ and xt (s) ∈ D for all s ∈ [−r, 0]. For example, ξ(t) can be a
vector (x⊤ (t), x⊤ (t − r))⊤ , then Dξ = D × D, s = 2n.
On the set Dξ we define continuous functions µk , k = 1, . . . , p, such that
∑
p
µk (ξ) ∈ [0, 1], µk (ξ) = 1 for all ξ ∈ Dξ . (3)
k=1
Suppose that on the set D equation (1) is represented in the form of TS system (see e. g.
[Tanaka & Wang, 2001] about ways of representation by TS systems):
∑
p
ẋ(t) = µk (ξ(t))(Ak x(t) + Akτ x(t − τ (t))). (4)
k=1
Here τ (t) : R+ → [0, r] is a piecewise continuous function, µk meet (3), Ak , Akτ are constant n × n-matrices
(k = 1, . . . , p).
3 Results
We use the following statement on asymptotic stability theorem with Razumikhin conditions (for a time-invariant
delay, see [Druzhinina & Sedova, 2017]):
Theorem 1 Assume that, for some positive numbers a and b, and symmetric positive definite matrix Q ∈ Rn×n
the(LMIs ) ( )
−aQ Ak Q −bQ Akτ Q
a) ≤ 0, ≤ 0,
∗ −Q ∗ −Q
( 1 )
r Φk + (a + b)Q Akτ Q
b) < 0,
∗ − 12 Q
hold for k = 1, . . . , p with Φk = Q(Ak + Akτ )⊤ + (Ak + Akτ )Q. Then the zero solution of (4) is uniformly
asymptotically stable.
It should be noted that under the conditions of Theorem 1 the function V (x) = x⊤ P x with P = Q−1 is the
Lyapunov function for (4).
Since the LMIs in Theorem 1 depend on some scalar parameters, they are not linear in all variables. Of
course, for each combination of those parameters we obtain the LMIs that are linear in the matrix variable Q.
But in any case the direct use of a LMI solver can be troublesome. A possible solution is to state a sequence of
minimization problems, each of which is the minimization of a linear function under constraints represented by
linear matrix inequalities. Such a problem can be therefore formulated as a semidefinite programming (SDP),
which is a convex optimization problem:
Step 1. Solve the following GEVP:
min α subject to
Q
Q > 0, Φk < αQ, k = 1, . . . , p;
∗
If the corresponding α is negative, we obtain matrix Q and go to Step 2.
Step 2. With the value Q, obtained in Step 1, compute the minimal positive a, b, meeting inequalities a) in
Theorem 1;
Step 3. With the value Q, obtained in Step 1, and a, b, obtained in Step 2, compute the maximal r = r∗ > 0
such that LMIs b) in Theorem 1 are feasible.
So Theorem 1 gives the sufficient stability conditions which can be efficiently checked by solving SDPs. Indeed,
if a solution Q of GEVP in Step 1 exists with α∗ < 0, then there exist a > 0, b > 0 meeting condition a) of
Theorem 1, as well as condition b) is also valid for sufficiently small r > 0. Moreover, the smaller a and b the less
506
restrictive condition b) becomes. If the maximal value of r is to find, we obtain the delay-dependent sufficient
stability conditions. If the value of delay is specified then step 3 should be replaced by checking whether LMIs
b) of Theorem 1 are feasible. The problem described above by steps 1–3 will be referred to as Problem 1.
Now let us discuss the problem of constructing the DA for system (4) with the asymptotically stable origin.
The goal is to obtain the largest possible estimate of the DA whose shape is fixed by some quadratic function.
For system (4) the problem has some special feature. Unless nonlinearities are naturally bounded for the
whole state space Rn of the nonlinear system, any TS representation will be only locally valid for a subset of the
state space. This is due to the fact that functions µk , which depend on the system nonlinearities, do not stand
for the convex-sum properties (3) outside the set D, thus invalidating any Lyapunov-based analysis beyond this
set. In particular, outside of the set D, the conditions of Theorem 1 do not ensure asymptotic stability of the
zero solution.
Moreover, the set D can have the sense of the domain of safe system operation. So, the problem is to construct
a set A ⊂ D such that the solutions beginning in the set should not leave the set D and must tend to equilibrium.
Finally, we use the following definition of DA (see the similar definitions in [Prasolov, 1998] and
[Kamenetskii, 1996]):
The direct domain of attraction (DA) for system (4) on a set D ⊂ Rn is the set A ⊂ D such that x(t0 + s) ∈ A
for all s ∈ [−r, 0] implies x(t; t0 , xt0 ) ∈ D for all t ≥ t0 and x(t; t0 , xt0 ) → 0 as t → +∞.
So, the solutions beginning at A: 1) must tend to equilibrium, and 2) should not leave the given set D.
Within the framework of the direct Lyapunov method, the DA is estimated by the sets bounded by the level
surfaces of the Lyapunov functions. The results on this topic obtained for ODE through the Lyapunov direct
method and other techniques are numerous (see e.g. [Kamenetskii, 1996] and references therein).
Here we use quadratic Lyapunov functions. Thus, we can find an invariant estimate A0 ⊂ A that defined in
the form A0 = B(P, c) = {x ∈ D : x⊤ P x ≤ c} with some c > 0 and a positive definite matrix P . Assuming that
the conditions of Theorem 1 are valid, we obtain that the set B(P, c) is invariant with respect to system (4) on
the set D, and B(P, c) ⊂ D implies therefore B(P, c) ⊂ A. That is the outermost set B(P, c) ⊂ D is the largest
valid estimation of the DA through the proposed method.
Furthermore, since B(P, c1 ) ⊂ B(P, c2 ) for c1 < c2 , that outermost set is defined by c∗ = max{c ∈ R+ :
B(P, c) ⊆ D}. The number c∗ can be numerically estimated via a gridding procedure by using the fact that
c∗ = min{V (x) : x ∈ ∂D}, where ∂D denotes the boundary of D.
Assume now that the set D is of the form
D = {x ∈ Rn : |xi | ≤ di , i ∈ I ⊂ {1, . . . , n}}. (5)
Notice that this is the typical form of a set on which the TS-representation corresponds to the original nonlinear
system.
Let us use the standard basis: Eij is the matrix with the i, j-th entry equal to 1 and all the rest are zero.
Then, we can translate the condition B(P, c) ⊂ D into LMIs. Indeed, it is easy to see that if for β > 0 we have
1
Eii ≤ βP for all i ∈ I, (6)
d2i
then B(P, 1/β) ⊂ D.
So, the method of computing an inner estimate of the DA is summarized in the following procedure:
Step 1. Solve Problem 1. If the problem is feasible, we obtain positive definite matrix P = Q−1 .
Step 2. With the value P , obtained in Step 1, compute
β ∗ = min β : (6)
(notice that β ∗ is positive). Then set c∗ = 1/β ∗ .
The problem described above by steps 1–2 will be refereed to as Problem 2.
As the result, we obtain the set A0 = B(P, c∗ ) ⊂ D which is invariant with respect to (4) and it is the largest
inner estimate of the DA whose shape is fixed by the function V (x).
A similar result can be obtained for the constraints of the form gk (x) = x⊤ Ck x − bk with a symmetric matrix
Ck and bk > 0: it is easy to check that any x ∈ B(P, 1/β) meets the inequality gk (x) 6 0 if b1k Ck 6 βP .
In accordance with the latter condition we can modify Step 2 in the procedure above.
507
4 Example
Consider the system
ẋ(t) = (3x2 (t − r) + 4x(t − r)y(t − r) + 2y 2 (t − r) + 0.5)x(t) + y(t − r)
(7)
ẏ(t) = −x(t − r) − y(t − r)
in the ellipse D = {(x, y) : 3x2 + 4xy + 2y ≤ 0.25}. System (7) can be represented in the form (4) with p = 2,
( ) ( ) ( )
1 0.5 0 2 0.75 0 1 2 0 1
A = , A = , Aτ = Aτ = .
0 0 0 0 −1 −1
Solving Problem 1, we obtain (with an accuracy of 4 decimal places):
( )
9.7387 −8.2377
Q= , α∗ = −0.1406;
−8.2377 9.7395
a = 1.9768, b = 3.3860; r∗ = 0.013.
Then, solving Problem 2 with ( )
−1 0.3608 0.3052
P =Q = ,
0.3052 0.3608
we obtain c∗ = 1/β ∗ = 0.0233. So for all r ∈ [0, r∗ ] the zero solution of (7) is asymptotically stable, and
A0 = B(P, c∗ ) ⊂ D is an inner estimate of the DA. The sets B(P, c∗ ) and D are presented in Fig. 1(left).
Remark 1. It is easy to see that for sufficiently large initial deviations from the equilibrium the corresponding
solutions tend to infinity, that is the zero solution of (7) has the bounded DA. Numerical experiments show
that the DA is not much larger than the estimate obtained. Some phase trajectories for r = 0 are presented in
Fig. 1(right).
Figure 1: (left) The sets B(P, c∗ ) (gray domain) and D for (7); (right) The set D and phase trajectories for (7)
(r = 0)
.
Remark 2. Numerical experiments show that as the number r > r∗ increases, the DA diminishes. Two phase
trajectories for (7) with r = 0.25 and initial conditions belonging to the set B(P, c∗ ) are presented on Fig. 2(left);
it can be seen from the figure that the trajectory beginning near the bound of B(P, c∗ ) does not tend to zero.
Moreover, for sufficiently large values of r (the numerical estimate of this value is approximately 0.3) the zero
solution of (7) is unstable. Two phase trajectories for (7) with the same small initial conditions and different
delays are presented on Fig. 2(right).
508
Notice that the estimation r∗ of the maximum allowable value of delay that preserves stability is far from
precise. This is due, among other things, to the rather conservative nature of the method of Razumikhin
functions. It is well known that the usage of the Lyapunov-Krasovskii functionals may yield less conservative
stability conditions. However, the functionals are inconvenient from the point of view of estimating the DA, as
noted above. In addition, when using the method of Razumikhin functions, the function τ (t) may vary quite
arbitrarily within a specified range.
Figure 2: (left) Two phase trajectories for (7) with r = 0.25 and initial conditions belonging the set
B(P, c∗ ); (right) Two phase trajectories for (7) with the same initial conditions and different delays ((x0 , y0 ) =
(0.01, −0.01), r = 0.1 (solid line) and r = 0.4 (dashed line)
5 Conclusion
In this paper, the direct Lyapunov method is used to obtain sufficient conditions for stability the zero solution of
nonlinear delay systems as well as to construct an estimation for the domain of attraction subject to constraints.
The assumption that on a subset of instant states the system is represented by a continuous-time Takagi–Sugeno
system and a special structure of a Lyapunov function allow the stability conditions to be expressed in terms
of SDP problems. Here, the simplest case of common quadratic Lyapunov function is considered. Similar ideas
can be used for other techniques of stability analysis, such as fuzzy Lyapunov function or piecewise Lyapunov
function approach (some recent results for ordinary differential and difference equations see e.g. [Xie et al., 2015,
González & Bernal, 2016]). These ways can lead to less conservative LMIs which express sufficient conditions
for stability.
References
[Druzhinina & Sedova, 2017] Druzhinina, O. V., & Sedova, N. O. (2017). Analysis of Stability and Stabilization
of Cascade Systems with Time Delay in Terms of Linear Matrix Inequalities, Journal of Computer and
Systems Sciences International, 56(1), 19–32. doi:10.1134/S1064230717010063
[González & Bernal, 2016] González, T., & Bernal, M. (2016). Progressively better estimates of the domain of
attraction for nonlinear systems via piecewise Takagi–Sugeno models: Stability and stabilization issues,
Fuzzy Sets and Systems, 297, 73–95. doi:10.1016/j.fss.2015.11.010
[Kamenetskii, 1996] Kamenetskii, V. A. (1996). Parametric Stabilization of Nonlinear Control Systems under
State Constraints. Autom. Remote Control, 57(10), 1427–1435.
509
[Khlebnikov et al., 2011] Khlebnikov, M. V., Polyak, B. T., & Kuntsevich, V. M. (2011). Optimization of Linear
Systems Subject to Bounded Exogenous Disturbances: The Invariant Ellipsoid Technique. Automation
and Remote Control, 72(11), 2227–2275. doi:10.1134/S0005117911110026
[Leng et al., 2016] Leng, S., Lin,W., & Kurths, J. Basin stability in delayed dynamics. Nature Scientific Re-
ports, 6. doi: 10.1038/srep21449
[Nesterov & Nemirovski, 1993] Nesterov, Y., & Nemirovski, A. (1993). Interior Point Polynomial Methods in
Convex Programming: Theory and Applications, Philadelphia: SIAM.
[Prasolov, 1998] Prasolov, A. V. (1998). On an Estimate of the Domain of Attraction for Systems with Aftereffect.
Ukrainian Mathematical Journal, 50(5), 761–769. doi:10.1007/BF02514329
[Razumikhin, 1956] Razumikhin, B. S. (1956). On stability of systems with delay. Prikl. Mat. Mekh., 20(4),
500–512 [in Russian].
[Sedova, 2008] Sedova, N. O. (2008). The global asymptotic stability and stabilization in nonlinear cascade
systems with delay. Russian Mathematics (Iz VUZ), 52(11), 60–69. doi:10.3103/S1066369X08110078
[Sedova, 2011] Sedova, N. O. (2011). On the principle of reduction for the nonlinear delay systems, Automation
and Remote Control, 72(9), 1864–1875. doi:10.1134/S0005117911090086
[Tanaka & Wang, 2001] Tanaka, K., & Wang, H. O. (2001). Fuzzy control systems design and analysis: a linear
matrix inequality approach. New York: Wiley.
[Xie et al., 2015] Xie, X.-P., Liu, Z.-W., & Zhu, X.-L. (2015). An efficient approach for reducing the conservatism
of LMI-based stability conditions for continuous-time TS fuzzy systems. Fuzzy Sets and Systems, 263,
71–81. doi:10.1016/j.fss.2014.05.020
510