Homogeneous Lyapunov functions for systems ... - Semantic Scholar
Homogeneous Lyapunov functions for systems with structured uncertainties G. Chesi1 , A. Garulli1 , A. Tesi2 , A. Vicino1 1
Dipartimento di Ingegneria dell’Informazione, Universit`a di Siena Via Roma, 56 - 53100 Siena, Italia 2
Dipartimento di Sistemi e Informatica, Universit`a di Firenze Via di S. Marta, 3 - 50139 Firenze, Italia
Abstract The problem addressed in this paper is the construction of homogeneous polynomial Lyapunov functions for linear systems with time-varying structured uncertainties. A sufficient condition for the existence of a homogeneous polynomial Lyapunov function of given degree is formulated in terms of a Linear Matrix Inequalities (LMI) feasibility problem. This condition turns out to be also necessary in some cases depending on the dimension of the system and the degree of the Lyapunov function. The maximum `∞ norm of the parametric uncertainty for which there exists a homogeneous polynomial Lyapunov function is computed by solving a generalized eigenvalue problem. The construction of such Lyapunov functions is efficiently performed by means of popular convex optimization tools for the solution of problems in LMI form. Comparisons with other classes of Lyapunov functions in numerical examples taken from the literature show that homogeneous polynomial Lyapunov functions are a powerful tool for robustness analysis.
Keywords: Robustness, Lyapunov function, Homogeneous forms, Linear Matrix Inequalities.
1
1
Introduction
Lyapunov functions are a standard tool for tackling robust analysis of linear systems affected by time-varying structured uncertainties. Quadratic Lyapunov functions have been considered by many authors, since a long time (see e.g. [1, 2, 3, 4]). More recently, relationships between quadratic stability and Integral Quadratic Constraints have been also investigated [5]. It is however widely recognized that quadratic functions lead to conservative estimates of the robust stability margin. Hence, nonquadratic Lyapunov functions have been addressed in the literature. Piecewise quadratic Lyapunov functions have been considered for both linear systems with timevarying perturbations [6] and switching linear systems [7]. Polyhedral Lyapunov functions have been introduced in [8] and successively considered by several authors. It has been shown that they are not conservative for robust analysis [9] and synthesis [10] of linear systems with time-varying structured uncertainties. The main drawback of polyhedral functions is that the computational burden required by their construction dramatically increases with the dimension of the system and the number of vertices of the polytope of matrices describing the uncertainty. Homogeneous polynomial Lyapunov functions are a viable alternative to the above classes of Lyapunov functions. The fact that this class of Lyapunov functions can improve robust stability results provided by quadratic Lyapunov functions has been recognized since long time [11]. In [12], the use of homogeneous polynomial functions to prove robust stability of linear systems with time-varying uncertainties has been considered and an approach based on the S-procedure has been proposed to enhance the robustness degree achieved by such functions. More recently, it has been shown that for a wide class of uncertain systems, including linear systems with time-varying structured uncertainties, robust stability is equivalent to the existence of a smooth Lyapunov function [13], and moreover that this Lyapunov function can always be chosen as a homogeneous polynomial of arbitrary degree [14]. This means that the class of homogeneous polynomial Lyapunov functions is not conservative for proving robust stability of linear systems with time-varying uncertainties. In this paper, the problem of constructing homogeneous polynomial Lyapunov functions is tackled by means of recently developed convex optimization techniques based on Linear Matrix Inequalities (LMI) [15, 16] and on certain properties of homogeneous forms which lead to a suitable matricial representation [17, 18]. First, a sufficient condition for the existence of a homogeneous polynomial Lyapunov function of given degree is formulated in terms of an LMI feasibility problem. The condition turns out to be also necessary in some cases, depending on the system dimension and the degree of the Lyapunov function. Then, the problem of computing the maximum ` ∞ norm of the parametric uncertainty for which a homogeneous Lyapunov function exists is addressed and a lower bound is obtained by solving a generalized eigenvalue problem [15]. Tightness of the lower bound is discussed, and it is shown to be related to the necessity of the existence condition. Another problem of interest is that of computing the Lyapunov function which guarantees the best
2
transient performance of an uncertain system: it is shown that also this problem can be easily tackled in the considered framework. Finally, several examples are presented, which demonstrate that homogeneous polynomial Lyapunov functions are a powerful tool for robustness analysis. The paper is organized as follows. Section 2 contains the problem formulation and preliminary material on homogeneous forms. The existence conditions for homogeneous polynomial Lyapunov functions are provided in Section 3, where relationships with previous work are also discussed. Maximization of the `∞ norm of the uncertainty is treated in Section 4, while the computation of the homogeneous polynomial Lyapunov function achieving the optimal transient performance is addressed in Section 5. Several numerical examples are presented in Section 6 and some concluding remarks are given in Section 7.
2
Problem formulation and preliminaries
Let us introduce the notation used throughout the paper. In : n × n identity matrix; A0 : transpose
of matrix A; A > 0 (A ≥ 0): symmetric positive definite (semidefinite) matrix A; A ⊗ B: Kronecker
product of matrices A and B; kxk: euclidean norm of vector x; kxk∞ , maxi |xi |: `∞ norm of x;
co{w1 , . . . , wr }: convex hull of vectors w i , i = 1, . . . , r. A function fm (x) is a homogeneous form of degree m in x ∈ Rn if
fm (x) =
X
ci1 ,i2 ,...,in xi11 xi22 . . . xinn ,
i1 +i2 +...+in =m
where i1 , i2 , . . . , in are nonnegative integers, and ci1 ,i2 ,...,in ∈ R are weighting coefficients. The form
fm (x) is said positive if fm (x) > 0 ∀x 6= 0, and nonnegative if fm (x) ≥ 0 ∀x.
2.1
Problem formulation
Let us consider the uncertain linear system x(t) ˙ = A(w(t))x(t), where A(w(t)) = and A0 , . . . , As ∈
Rn×n
Ã
A0 +
s X i=1
wi (t)Ai
(1) !
,
(2)
are given matrices. The uncertain time-varying parameter vector w(t) =
0
(w1 (t), . . . , ws (t)) is assumed to satisfy for all t ≥ 0 the constraint w(t) ∈ W , co{w 1 , . . . , wr }
(3)
where w i ∈ Rs , i = 1, . . . , r, are given vectors.
The basic problem addressed in this paper is the construction of a Lyapunov function proving global asymptotic stability of system (1)-(3). The attention is restricted to a special class of Lyapunov 3
functions: the homogeneous polynomial forms of degree 2m, hereafter referred to as Homogeneous Polynomial Lyapunov Functions (HPLFs) and denoted by v2m (x). More specifically, the aim is to find v2m (x) such that (i) v2m (x) > 0 for all x 6= 0; (ii) v˙ 2m (x) < 0 for all x 6= 0 and for all w(t) ∈ W. Before proceeding, let us introduce the key representation of homogeneous forms that will be exploited throughout the paper.
2.2
Square Matricial Representation of homogeneous forms
Let g2m (x) be a homogeneous form of degree 2m. Then, the Square Matricial Representation (SMR) of g2m (x) is defined as 0
g2m (x) = x{m} Cg x{m} ,
(4)
where x{m} ∈ Rd is the base vector of homogeneous forms of degree m in x and Cg = Cg0 ∈ Rd×d is
a coefficient matrix. Let Ei ∈ R(n+1−i)×n be the matrix composed by the rows i, i + 1, . . . n of In .
Then, the base vector x{m} containing all monomials of x x1 x{m−1} x2 (E2 x){m−1} x{m} = .. . {m−1} xn−1 (En−1 x) xm n
degree m is defined as if m = 1,
otherwise.
It is easy to check that the dimension d of x{m} is given by d=
(n + m − 1)! . (n − 1)!m!
(5)
An important property of the SMR is that matrix Cg is not unique. Indeed, all the matrices Cg satisfying (4) can be parameterized as Cg + L, with L ∈ L, where L is the set of matrices describing
the null form, i.e.
n o 0 L = L = L0 ∈ Rd×d : x{m} Lx{m} = 0 ∀x ∈ Rn .
(6)
1 (n + 2m − 1)! dL = d(d + 1) − . 2 (n − 1)!(2m)!
(7)
In [18] it has been shown that L is a linear space whose dimension is given by
4
Then, the family of matrices Cg describing g2m (x) can be parameterized affinely as Cg (α) = Cg + L(α),
(8)
where α ∈ RdL is a vector of free parameters and L : RdL → L is a linear parameterization of L. Hence, the Complete SMR (CSMR) of g2m (x) is defined as 0
g2m (x) = x{m} Cg (α)x{m} .
(9)
In the sequel, we will refer to Cg and Cg (α) as SMR and CSMR matrix of g2m (x), respectively.
3
Existence conditions for homogeneous Lyapunov functions
In this section, conditions for the existence of a HPLF for system (1)-(3) are provided. Such conditions are based on the SMR of homogeneous forms previously introduced and are expressed as Linear Matrix Inequalities.
3.1
Sufficient condition
For a generic system x(t) ˙ = Ax(t), let us introduce the extended matrix A {m} ∈ Rd×d , defined by d {m} ∂x{m} x = Ax = A{m} x{m} . dt ∂x
(10)
An expression of the extended matrix is provided below in terms of Kronecker products. n Lemma 1 Let x[m] = x | ⊗x⊗ {z. . . ⊗ x} , and Km ∈ R
m ×d
be the matrix satisfying
m
x[m] = Km x{m} .
(11)
Then, the extended matrix A{m} of A ∈ Rn×n is given by A{m} =
0 0 (Km Km )−1 Km
Ãm−1 X i=0
Inm−1−i ⊗ A ⊗ Ini
!
Km .
Proof. See the Appendix. The following useful properties of the extended matrix A{m} hold [11, 12]: (I) Let A, B ∈ Rn×n and α, β ∈ R. Then (αA + βB){m} = αA{m} + βB{m} .
5
(12)
(II) Let A0,{m} and Ai,{m} denote respectively the extended matrices of A0 and Ai , according to (10). Consider the extended system ! Ã s X wi (t)Ai,{m} x{m} (t) = A(w(t)){m} x{m} (t). x˙ {m} (t) = A0,{m} + i=1
0
Then, for any v2m (x) = x{m} V x{m} , one has ¯ ¯ ¯ ¯ d d ¯ v2m (x)¯ = v2m (x)¯¯ . dt dt x˙ {m} =A(w(t)){m} x{m} x=A(w(t)) ˙ x
The following result provides a sufficient condition for establishing the existence of a HPLF of degree 2m for system (1)-(3). Theorem 1 Let A¯j = A(wj ), j = 1, . . . , r and let A¯j,{m} denote the extended matrix of A¯j . If the system of LMIs
(
V > 0, −V A¯j,{m} − A¯0
j,{m} V
− L(αj ) > 0,
(13)
j = 1, . . . , r 0
admits a feasible solution V = V 0 ∈ Rd×d , αj ∈ RdL , j = 1, . . . , r, then v2m (x) = x{m} V x{m} is a
HPLF for (1)-(3).
0
Proof. If (13) admits solution, then v2m (x) = x{m} V x{m} is positive definite. By differentiating v2m (x) along the trajectories of x˙ {m} = A¯j,{m} x{m} , one gets ¯ ¯ d v2m (x)¯¯ dt ¯j,{m} x{m} x˙ {m} =A
= x{m}
0
= x{m}
0
³
³
´ V A¯j,{m} + A¯0j,{m} V x{m}
´ V A¯j,{m} + A¯0j,{m} V + L(αj ) x{m} ,
(14)
which is negative definite by (13). Observe that the CSMR (8)-(9) has been used in (14). Now, ´ ³ P notice that, due to property (I), matrices A¯j,{m} = A0,{m} + si=1 wij Ai,{m} are the vertices of P the set of extended matrices A0,{m} + si=1 wi (t)Ai,{m} , w(t) ∈ W. Convexity of this set of matrices
implies that v2m (x) is a common HPLF for all matrices of the set. Hence, property (II) allows one to conclude that v2m (x) is also a HPLF for system (1).
¥
Remark 1. It is worth observing that (13) is an LMI feasibility problem [15], which can be solved by efficient computational tools based on convex optimization [16, 19]. Notice that the free variables in (13) are matrix V = V 0 ∈ Rd×d and vectors αj ∈ RdL , j = 1, . . . , r, and therefore their number
is equal to d(d + 1)/2 + rdL − 1, with d and dL given by (5) and (7), respectively (the term −1 is due to the fact that V can be properly scaled, for example by setting to 1 the first entry of the
matrix).
6
3.2
Necessary and sufficient condition
In some cases, depending on the dimension n of x and the degree 2m of v2m (x), the sufficient condition provided by Theorem 1 is also necessary for the existence of a HPLF for (1)-(3). The next result is central for establishing such a necessary and sufficient condition. Lemma 2 Consider the equivalence set E, defined as E = {(n, 2) ∀n} ∪ {(2, 2m) ∀m} ∪ {(3, 4)} .
(15)
Let x ∈ Rn and g2m (x) be a nonnegative homogeneous form of degree 2m. If (n, 2m) ∈ E, then g2m (x) can be written as the sum of squares of homogeneous forms of degree m. Moreover, if g2m (x) is positive, then it admits a positive definite SMR, i.e. there exists a positive definite matrix 0
V ∈ Rd×d such that g2m (x) = x{m} V x{m} . Proof. See the Appendix. Theorem 2 Let (n, 2m) belong to E in (15). Then, there exists a HPLF of degree 2m for system (1) if and only if the set of LMIs (13) admits a feasible solution.
Proof. Obviously, it must be proven only that if (n, 2m) ∈ E and there exists a HPLF v 2m (x)
for (1), then (13) admits a feasible solution. By assumption we have that v2m (x) is positive definite, and −v˙ 2m (x) is positive definite for any uncertainty w(t) ∈ W. Hence, due to Lemma 2, 0
there exists V > 0 such that v2m (x) = x{m} V x{m} , and there exist Vj > 0, j = 1, . . . , r, such that the time derivative of v2m (x) evaluated along the trajectories of x˙ {m} = A¯j,{m} x{m} satisfies 0
−v˙ 2m (x) = x{m} Vj x{m} . Since the parameterization of the SMR matrices given by (6) and (8) is complete, there exists αj such that Vj = −V A¯j,{m} − A¯0j,{m} V − L(αj ). Therefore, (13) admits a
¥
feasible solution.
Remark 2. It is worth observing that, though it is well known that there exist homogenous forms that cannot be written as sums of squares, experience shows that occurrences of such forms are quite rare in practice (see e.g. [18, 20]).
3.3
Relationships with previous work on HPLF
The use of HPLFs for robust analysis of system (1)-(3) has been addresses in [12]. In this paper, a result based on the S-procedure is exploited to provide a sufficient condition for the existence of a HPLF of degree 2m. Although the condition is formulated as the minimization of a nondifferentiable convex function, it can be easily rewritten as an LMI feasibility problem of the form V > 0, d−n X ¯j,{m} − A¯0 −V A V − βij Fi > 0, j = 1, . . . , r j,{m} i=1
7
(16)
where Fi ∈ Rd×d , i = 1, . . . , d − n are suitable matrices. The difference between (16) and (13) is
that the number of free variables βij for each LMI in (16) is dZ = d − n, which turns out to be much
smaller than dL , the dimension of vector αj in (13), as clarified by Table 1. The larger number n=2 m=2
3
4
5
dL = 1
dZ = 1
6
3
20
6
50
10
3
3
2
27
7
126
16
420
30
4
6
3
75
12
465
31
1990
65
5
10
4
165
18
1310
52
7000
121
Table 1: Values of dL and dZ for some n and m. of free variables in the LMIs makes the sufficient condition provided by Theorem 1 much more powerful. Moreover, the fact that the SMR of homogeneous forms is complete, allows to formulate the necessary and sufficient condition of Theorem 2 for the cases in which (n, 2m) ∈ E. On the contrary, the condition provided in [12] is necessary only for the case n = m = 2, in which the two
parameterization of homogeneous forms coincide (see Table 1).
4
Computation of the `∞ 2m-HPLF stability margin
In the analysis of uncertain systems of type (1)-(2), a key problem is that of finding the largest value of the positive scalar γ for which there exists a HPLF of degree 2m for all w(t) belonging to the scaled perturbation set γW. In order to simplify the presentation, we restrict our attention to the case in which the perturbation set is the `∞ box, i.e. Bγ = {q ∈ Rs : |qi | ≤ γ, i = 1, . . . , s} .
(17)
s
Clearly, Bγ is equal to γ co{u1 , . . . , u2 }, where uj , j = 1, . . . , 2s , are the vertices of the unit `∞ ball B1 . Said another way, the aim is to compute
∗ γ2m = sup{γ : ∃ v2m (x) for (1)-(2), with w(t) ∈ Bγ }.
(18)
∗ will be referred to as the ` In the following, γ2m ∞ 2m-HPLF stability margin for system (1)-(2).
We also consider the related problem of computing the following quantity κ∗2m = sup{κ : ∃ v2m (x) for (1)-(2), with w(t) ∈ B¯κ },
(19)
B¯κ = {q ∈ Rs : 0 ≤ qi ≤ κ, i = 1, . . . , s} .
(20)
where
∗ because the perturbation set is restricted to the positive orthant Notice that κ∗2m differs from γ2m
of the parameter space. For this reason, it will be referred to as the ` ∞ 2m-HPLF positive stability 8
∗ and κ∗ are bounded from above margin. It is worth recalling that the above stability margins γ2m 2m
by the well-known `∞ state space parametric stability margin (see for example [21, 22]). Let us first focus the attention on problem (18)-(17). Define the matrices A˜j = A(uj ) − A0 , where uj , j = 1, . . . , 2s are the vertices of the unit `∞ ball B1 . Moreover, let A˜j,{m} denote the extended ∗ matrix of A˜j . The following result shows that a lower bound γˆ2m of the `∞ 2m-HPLF stability ∗ can be computed by solving a quasiconvex optimization problem. margin γ2m ∗ be defined as Theorem 3 Let γˆ2m ∗ )−1 = (ˆ γ2m
inf z ∈ R; V = V 0 ∈ Rd×d ;
z
αj ∈ RdL , j = 0, . . . , 2s V > 0, z > 0 −V A0,{m} − A0,{m} 0 V − L(α0 ) > 0 s.t. ´ ³ z −V A 0 0 j ˜ ˜0 0,{m} − A0,{m} V − L(α ) > V Aj,{m} + Aj,{m} V + L(α )
(21)
j = 1, . . . , 2s .
∗ ≤ γ ∗ . Moreover, if (n, 2m) belongs to E, then γ ∗ = γ∗ . Then, γˆ2m ˆ2m 2m 2m 0
Proof. Let v2m (x) = x{m} V x{m} . The time derivative of v2m (x) evaluated for w(t) = z −1 uj , is given by ¯ · ³ ´ ³ ´0 ¸ ¯ d 0 v2m (x)¯¯ = x{m} V A0,{m} + z −1 A˜j,{m} + A0,{m} + z −1 A˜j,{m} V x{m} dt w(t)=z −1 uj
h ³ ´ i 0 = z −1 x{m} z V A0,{m} + A00,{m} V + L(α0 ) + V A˜j,{m} + A˜0j,{m} V + L(αj ) x{m} ,
where the CSMR has been exploited by adding the term zL(α0 ) + L(αj ) inside the square brackets. Hence, the constraint in (21) guarantees that v˙ 2m (x) is negative definite. Since this holds for all ∗ ≤ γ∗ . w(t) such that kw(t)k∞ ≤ z −1 , one has that γˆ2m 2m
∗ , Let us assume now that (n, 2m) belongs to E. From definition (18), for any z such that z −1 < γ2m
there exists a positive definite v2m (x) such that −v˙ 2m (x) is positive definite for w(t) : kw(t)k∞ ≤ z −1 . From Lemma 2, it follows that −v˙ 2m (x) admits a positive definite SMR matrix for any w(t) ∈ Bz −1 and therefore the constraint in problem (21) admits a feasible solution for any
∗ . Hence, γ ∗ = γ∗ . z −1 < γ2m ˆ2m 2m
¥
Remark 3. Problem (21) is a Generalized Eigenvalue Problem (GEVP), which has been proven to be a quasiconvex optimization problem and can be tackled by efficient optimization tools [15, 16]. Notice that the LMI − V A0,{m} − A0,{m} 0 V − L(α0 ) > 0 in (21) guarantees that the matrix multi-
plying the generalized eigenvalue z is positive, as required by the standard form of the GEVP (see
9
[15]). Observe that the number of free variables in the GEVP (21) is equal to d(d+1)/2+(2 s +1)dL . A result similar to Theorem 3 can be obtained for the `∞ 2m-HPLF positive stability margin κ∗2m in (19)-(20). Here, we report the result for the simpler case of a segment of matrices, because it has been widely addressed in the literature and it will be useful in the examples presented in the next section. The aim is to compute the largest value of κ such that there exist a HPLF of degree 2m for the matrices A0 + w(t)A1 , with 0 ≤ w(t) ≤ κ. Corollary 1 Let s = 1 in (19)-(20). Define κ ˆ ∗2m as (ˆ κ∗2m )−1 =
z min z ∈ R; V = V 0 ∈ Rd×d ;
α 1 , α 2 ∈ R dL V >0 −V A0,{m} − A00,{m} V − L(α1 ) > 0 s.t. ³ ´ z −V A 0 1) > V A 0 2 − A V − L(α 0,{m} 1,{m} + A1,{m} V + L(α ) 0,{m}
(22)
ˆ ∗2m = κ∗2m . Then, κ ˆ ∗2m ≤ κ∗2m . Moreover, if (n, 2m) belongs to E, then κ
5
Construction of the optimal performance HPLF
Another problem of interest in the robustness analysis of uncertain systems is that of determining the Lyapunov function that achieves the best transient performance (see, for example, the classical textbook [23]). For a given Lyapunov function v(x), one can define the transient performance index of system (1)-(3) as λ(v) = sup
sup
x∈Rn \0 w(t)∈W
v(x) ˙ . v(x)
(23)
From (23) one has v(x(t)) ≤ v(x(t0 ))eλ(v)(t−t0 ) , thus establishing the rate of decrease of v.
Therefore, it is natural to select among all feasible Lyapunov functions, the one that minimizes λ(v). For systems of the type (1)-(3), this problem has been addressed in [24] within the class of quadratic Lyapunov functions. Here, the aim is to select the optimal Lyapunov function among HPLFs of degree 2m, i.e. to compute λ∗2m = inf λ(v2m ). v2m
The following result shows that an upper bound of λ∗2m can be obtained by solving a GEVP.
10
(24)
ˆ ∗ be defined as Theorem 4 Let λ 2m ˆ∗ = λ 2m
z inf 0 d×d z ∈ R; V = V ∈ R ;
αj ∈ RdL , j = 1, . . . , r ( V > 0, s.t. zV > V A¯j,{m} + A¯0j,{m} V + L(αj )
(25) j = 1, . . . , r.
ˆ ∗ ≥ λ∗ . Moreover, if (n, 2m) belongs to E, then λ ˆ ∗ = λ∗ . Then, λ 2m 2m 2m 2m Proof. Let us assume that there exists z ∈ R and V = V 0 ∈ Rd×d satisfying the constraint in (25), 0
for some vectors αj ∈ RdL . Then, by setting v2m (x) = x{m} V x{m} and differentiating along the trajectories of x˙ {m} = A¯j,{m} x{m} , one has ³ ´ ¯ {m} 0 V A 0 ¯ ¯ x + A V x{m} j,{m} j,{m} v˙ 2m (x) ¯¯ = = {m} 0 V x{m} v2m (x) ¯x˙ {m} =A¯j,{m} x{m} x ³ ´ 0 x{m} V A¯j,{m} + A¯0j,{m} V + L(αj ) x{m} ≤z = 0 x{m} V x{m} for all j = 1, . . . , r. Therefore, by exploiting convexity of W, one can conclude that v˙ 2m /v2m < z
for all w(t) ∈ W and hence λ∗2m ≤ z. Minimizing with respect to z, one obtains the upper bound ˆ∗ . λ 2m Now, let (n, 2m) ∈ E. Definition (24) says that for any z > λ∗2m there exists a positive definite v2m (x) such that v˙ 2m (x)/v2m (x) < z, for all x and for all w(t) ∈ W. Then, Lemma 2 implies
that zv2m (x) − v˙ 2m (x) admits a positive definite SMR matrix for any w(t) ∈ W, and therefore the
constraint in problem (25) admits a feasible solution for any z > λ∗2m .
6
¥
Examples
The following examples illustrate the application of the methodology for the construction of HPLFs described in the previous section.
6.1
Example 1
This example is taken from [12]. Let us consider the matrices à ! à ! 0 1 0 0 A0 = , A1 = −2 −1 −1 0 and assume we want to compute the `∞ 2m-HPLF positive stability margin κ∗2m , defined in (19)(20). Since in this example n = 2, the lower bound provided by the GEVP (22) is tight, i.e. 11
κ ˆ ∗2m = κ∗2m , as assured by Corollary 1. In [12], it has been observed that a quadratic Lyapunov function exists only for κ < 3.82, and that this bound can be improved to κ < 5.73 by solving problem (16) with m = 2. Table 2 shows the values of κ∗2m obtained from (22), for some different m. m
2
3
4
5
6
7
8
κ∗2m
5.7393
6.2135
6.3982
6.6469
6.6581
6.7801
6.7962
Table 2: Values of κ∗2m for different m in Example 1. For the case m = 3, vector x{3} and extended matrices A0,{3} and A1,{3} are given by
x
{3}
x31
0
3
0
0
x2 x −2 −1 2 0 1 2 = , A0,{3} = x1 x22 0 −4 −2 1 3 x2 0 0 −6 −3
Matrix L(α) is equal to
0
0
0 2α1 L(α) = −α1 α2 −α2 −α3
0 0
α2
2α3
−α3 . 0 0
The matrix V obtained by solving (22) for m = 3 turns out to be 1.0000 0.1524 −0.1965 −0.0174 0.1524 1.2187 0.2086 −0.0280 V = −0.1965 0.2086 0.2938 0.0231 −0.0174 −0.0280 0.0231 0.0099
corresponding to the HPLF
0
−1 0 0 0 , A1,{3} = . 0 −2 0 0 0 0 −3 0
−α1 −α2
0
0
v6 (x) = x61 + 0.3048x51 x2 + 0.8258x41 x22 + 0.3825x31 x32 + 0.2378x21 x42 + 0.0463x1 x52 + 0.0099x62 . The same example has been considered also in [25] and [6]. In [25], a polyhedral Lyapunov function has been constructed, guaranteeing asymptotic stability for κ = 6. In [6], a piecewise quadratic Lyapunov function achieving stability for κ = 6.2, has been obtained via a sequence of LMI optimizations and a grid search over two free parameters. As it can be seen from Table 2, these bounds are improved by all the HPLFs of degree greater than 2. By the way, notice that the HPLF is obtained by solving one single GEVP, with 12 (3m2 + m + 2) free variables.
12
6.2
Example 2
This example shows that the gap in the Lyapunov stability margin between our technique and the one presented in [12] can be very large. Let us consider the matrices 0 1 0 −2 0 −1 , A1 = 1 −10 A0 = 0 0 1 3 −1 −2 −4 3 −4 2
(26)
and assume we want to solve the same problem as in Example 1. Quadratic stability is guaranteed only for κ ≤ 1.9042. Let us compute a HPLF of degree 4. Being n = 3 and m = 2 it follows that d = 6 and dL = 6. Then, vector x{2} and matrix L(α) are 0 0 x21 0 x1 x2 2α1 x1 x3 α2 , L(α) = 0 x{2} = 2 0 −α1 x2 −α2 −α4 x2 x3 2 −α3 −α5 x3
respectively given by: 0 α2
−α1 −α2 −α3 0
2α3
α4
α4
0
α5
0
0
−α6
−α4 −α5 α5 0 . 0 −α6 2α6 0 0 0
The GEVP (22) returns the lower bound κ ˆ ∗4 = 75.1071. Moreover, since (n, 2m) = (3, 4) ∈ E, we
have from Corollary 1 that κ ˆ ∗4 = κ∗4 .
Using the approach proposed in [12], one finds that the maximum κ for which robust stability is guaranteed is equal to 17.8347. The remarkable difference with respect to our approach is due to the fact that the parameterization of homogeneous forms adopted in [12] is not complete, as discussed in Subsection 3.3. Specifically, it is easy to see that in [12] only the parameters α 1 , α3 , α6 in the matrix L(α) are considered (see also Table 1).
6.3
Example 3
Let us consider the differential equation ˙ + k(t)ξ(t) = 0 ¨ + ξ(t) ξ(t) and assume that we want to compute the maximum κ such that the solution remain bounded, for all 0 ≤ k(t) ≤ κ. By adopting a strategy similar to that proposed in [26], one can determine numerically
the maximum value of κ. In particular, one has that such value must satisfy the equation ½ ¾ √ √ £ ¤ 1 k exp − √ π − arctan( 4κ − 1) − 1 = 0 4κ − 1 and is equal to κ∗ ≈ 3.0448.
13
Clearly, the problem can be easily tackled in the framework of HPLF by solving (19)-(20) with s = 1 and A0 =
Ã
0
1
0 −1
!
,
A1 =
Ã
0 0 −1 0
!
.
The values of κ∗2m obtained from (22) are reported in Table 3. Notice that the values exhibit a growth m
1
2
3
4
5
6
κ∗2m
1.0000
1.5000
1.9997
2.2952
2.4009
2.5053
m
7
8
9
10
11
12
κ∗2m
2.6102
2.6601
2.6954
2.7466
2.7726
2.7917
Table 3: Values of κ∗2m for different m in Example 3. towards the theoretical upper bound κ∗ . On the other hand, the results in [13, 14] guarantee that for any κ < κ∗ there exists a HPLF of suitable degree, hence it is expected that κ∗2m converges to κ∗ as m approaches infinity. Another interesting byproduct of the treatment in [26] is that one can calculate the worst-case sequence k(t) (which consists of suitable switchings between k(t) = 0 and k(t) = κ ∗ ) and the corresponding trajectory of the system, that represents the limit curve towards which the level surfaces of the Lyapunov functions are expected to tend. This is confirmed by Figure 1a, where the level curves of the obtained HPLFs for different values of m are depicted. Figure 1b clearly shows that the HPLF corresponding to m = 12 is very close to the limit trajectory predicted by the theory.
6.4
Example 4
This example has been considered by several authors [27, 24, 10]. Let us consider system (1)-(2) with s = 1 and matrices A0 =
Ã
0
1
−1 −1
!
,
A1 =
Ã
0 0 1 0
!
.
The maximum γ for which there exists a Lyapunov function for system (1)-(2), with w(t) ∈ B γ , has
been calculated for different classes of Lyapunov functions. Notice that the ` ∞ state space stability
∗ ≤ 1, for all margin is given by ρ = 1, because A0 + A1 is not asymptotically stable, and hence γ2m √ m. In [27], it has been proven that quadratic Lyapunov functions can achieve γ2∗ = 3/2. In [10], a
polyhedral Lyapunov function has been provided, which guarantees stability for |w(t)| ≤ γ = 0.98.
The level curves of this Lyapunov function have 30 vertices.
By applying Theorem 3 for different values of m, one can construct HPLFs of increasing degree, trying to improve the value of the 2m-HPLF stability margin. Observe that also in this example ∗ because n = 2. For m = 2, one obtains the lower bound provided by (21) always coincides with γ2m
14
1.5
1.5
1
1
0.5
0.5
x2
x2
limit trajectory 0
0
−0.5
−0.5
m=12
m=4 m=3 m=2 −1
−1.5 −1.5
−1
m=1
−1
−0.5
0
x
0.5
1
−1.5 −1.5
1.5
−1
−0.5
0
x
1
(a) level curves of v2m (x) for different m
0.5
1
1.5
1
(b) level curve of v24 (x) and limit trajectory (dashed)
Figure 1: Level curves of HPLFs in Example 3. ∗ = 1, which means that it is γ4∗ = 0.9771, but for m ≥ 3 the result returned by the GEVP is γˆ2m
possible to construct a HPLF of degree 6 (or more) for |w(t)| ≤ γ and γ arbitrarily close to 1.
In particular, for m = 3 the following matrix V has 1.0000 2.9987 2.9987 14.1343 V = 1.9974 12.4475 0.2491 2.1782 corresponding to γ6∗ = 1 and the HPLF
been obtained
1.9974 0.2491
12.4475 2.1782 14.5114 3.5064 3.5064 1.3106
v2m (x) = x61 + 5.9974x51 x2 + 18.1291x41 x22 + 25.3932x31 x32 + 18.8678x21 x42 + 7.0128x1 x52 + 1.3106x62 .
6.5
Example 5
This example concerns the computation of an optimal performance HPLF, in the sense explained in Section 5. Let us consider the helicopter model originally proposed in [28], and the robust controller designed in [29]. The resulting closed-loop uncertain system is given by −0.1027
0.1762
0
0
0.1995
−0.3446
−0.7275 − 0.1640q (t) −1.4280 + 0.2699q (t) 2.1852 + 0.4511q (t) 0.7218 + 0.4308q (t) 3 3 3 3 x(t) ˙ = 1.1689 −0.4446 + q1 (t) −3.6757 −3.0841 + q2 (t) 1
0
x(t)
with |q1 (t)| ≤ 0.2192, |q2 (t)| ≤ 1.2031, |q3 (t)| ≤ 2.0673. Clearly, the system can be easily written in
form (1)-(3) with r = 8.
15
In [24], the problem of computing the quadratic Lyapunov function achieving the best transient performance, defined as in (23), has been addressed. The value obtained was λ = −0.3839. By solving the GEVP problem (25) with m = 2, one obtains the HPLF
v4 (x) = x41 + 0.2184x31 x2 − 2.8083x31 x3 − 4.0483x31 x4 + 1.1458x21 x22 + 0.1831x21 x2 x3
+0.3233x21 x2 x4 + 9.7670x21 x23 + 16.4758x21 x3 x4 + 9.7066x21 x24 + 0.0283x1 x32
−2.5961x1 x22 x3 − 3.9673x1 x22 x4 − 0.1303x1 x2 x23 − 3.6305x1 x2 x3 x4 − 3.3848x1 x2 x24
−14.8574x1 x33 − 37.3189x1 x23 x4 − 31.4426x1 x3 x24 − 11.2239x1 x34 + 0.1064x42
+0.1657x32 x3 + 0.2465x32 x4 + 2.1908x22 x23 + 5.1952x22 x3 x4 + 4.2363x22 x24 −0.0017x2 x33 + 3.0142x2 x23 x4 + 6.4605x2 x3 x24 + 3.4452x2 x34 + 8.8681x43
+29.0364x33 x4 + 36.5118x23 x24 + 20.6917x3 x34 + 5.8479x44
ˆ ∗ = −0.8889. Therefore, once again it can be seen that the HPLFs outperform the achieving λ 4
classic approaches based on quadratic Lyapunov functions. Moreover, this example shows that the proposed LMI-based procedures are able to handle also quite complex systems (in this case, a fourth order uncertain system with 8 vertices).
7
Conclusions
The construction of homogeneous polynomial Lyapunov functions for linear systems with timevarying structured uncertainties has been addressed via LMI optimization techniques. Maximization of the size of the parametric uncertainty can also be efficiently tackled in the proposed framework. With respect to previous work on this class of Lyapunov functions, better results have been obtained by exploiting a complete parameterization of homogeneous forms of given degree. Moreover, this allows one to formulate necessary conditions for the existence of a homogeneous polynomial Lyapunov function in some cases (for example, for all second order systems). Comparisons with other classes of Lyapunov functions are also very promising. Quadratic Lyapunov functions are often outperformed by HPLFs. On the other hand, polyhedral Lyapunov functions seem to require level surfaces with a very high number of vertices, even in cases that can be handled by HPLFs of quite low degree (and hence with a limited number of coefficients). Further investigations on the potentialities of homogeneous polynomial Lyapunov functions will be the subject of future research.
References ˘ [1] D. D. Siljak. Nonlinear Systems: Parametric Analysis and Design. John Wiley & Sons, New York, 1969.
16
[2] K. S. Narendra and J. H. Taylor. Frequency Domain Criteria for Absolute Stability. Academic, New York, 1973. [3] E. S. Piatnitskiy and V. I. Skorodinskiy. Numerical methods of Lyapunov functions construction and their application to the absolute stability problem. Systems and Control Letters, 2:130–135, 1982. [4] K. Zhou and P. P. Khargonekar. Stability robustness bounds for linear state-space models with structured uncertainties. IEEE Transactions on Automatic Control, 32:621–623, 1987. [5] A. Megretski and A. Rantzer. System analysis via integral quadratic constraints. IEEE Transactions on Automatic Control, 42:819–830, 1997. [6] L. Xie, S. Shishkin, and M. Fu. Piecewise Lyapunov functions for robust stability of linear time-varying systems. Systems and Control Letters, 31:165–171, 1997. [7] M. Johansson and A. Rantzer. Computation of piecewise quadratic Lyapunov functions for hybrid systems. IEEE Transactions on Automatic Control, 43:555–559, 1998. [8] R. K. Brayton and C. H. Tong. Stability of dynamical systems: a constructive approach. IEEE Transactions on Circuits and Systems, 26:224–234, 1979. [9] A. P. Molchanov and E. S. Piatnitskiy. Lyapunov functions specifying necessary and sufficient conditions of absolute stability of nonlinear nonstationary control system. Automat. Remote Contr., 47:344–354 (part I), 443–451 (II), 620–630 (III), 1986. [10] F. Blanchini. Nonquadratic Lyapunov functions for robust control. Automatica, 31:451–461, 1995. [11] R. W. Brockett. Lie algebra and Lie groups in control theory. In D.Q. Mayne and R.W. Brockett, editors, Geometric Methods in Systems Theory, pages 43–82. Dordrecht, Reidel, 1973. [12] A. L. Zelentsovky. Nonquadratic Lyapunov functions for robust stability analysis of linear uncertain systems. IEEE Transactions on Automatic Control, 39:135–138, 1994. [13] W. P. Dayawansa and C. F. Martin. A converse Lyapunov theorem for a class of dynamical systems which undergo switching. IEEE Transactions on Automatic Control, 44:751–760, 1999. [14] D. Angeli. LMI conditions for the determination of polynomial Lyapunov functions. Technical report, DSI - University of Florence, Italy, 2001. [15] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix Inequalities in System and Control Theory. SIAM, Philadelphia, 1994. 17
[16] Y. Nesterov and A. Nemirovsky. Interior Point Polynomial Methods in Convex Programming: Theory and Applications. SIAM, Philadelphia, 1993. [17] G. Chesi, A. Tesi, A. Vicino, and R. Genesio. An LMI approach to constrained optimization with homogeneous forms. Systems and Control Letters, 42:11–19, 2001. [18] G. Chesi, A. Garulli, A. Tesi, and A. Vicino. LMI-based techniques for solving quadratic distance problems. In Proc. of 40-th IEEE Conf. on Decision and Control, Orlando, December 2001. [19] P. Gahinet, A. Nemirovski, A. J. Laub, and M. Chilali. LMI Control Toolbox. The Mathworks Inc., 1995. [20] P. Parrilo. Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization. Ph.D. Thesis, California Institute of Technology, 2000. [21] S. P. Bhattacharyya, H. Chapellat, and L. H. Keel. Robust Control: The Parametric Approach. Prentice Hall, NJ, 1995. [22] A. Vicino, A. Tesi, and M. Milanese. Computation of non conservative stability perturbation bounds for systems with nonlinearly correlated uncertainties. IEEE Transactions on Automatic Control, 35:835–841, 1990. [23] J. C. Hsu and A. U. Meyer. Modern Control Principles and Applications. McGraw-Hill Book Company, New York, 1968. [24] A. Olas. Construction of optimal Lyapunov function for systems with structured uncertainties. IEEE Transactions on Automatic Control, 39:167–171, 1994. [25] F. Blanchini and S. Miani. On the transient estimate for linear systems with time-varying uncertain parameters. IEEE Transactions on Circuits and Systems, Part I, 43:592–596, 1996. [26] R. W. Brockett. Finite Dimensional Linear Systems. John Wiley and Sons, NY, 1970. [27] B. Radziszewski. O Najlepszej Funkcji Lapunowa. IFTR Reports, 26, 1977. [28] K. S. Narendra and S. S. Tripathi. Identification and optimization of aircraft dynamics. Journal of Aircraft, 10:193–199, 1973. [29] Y. H. Chen and J. S. Chen. Robust control of uncertain systems with time-varying uncertainty: a computer-aided setting. In Lecture Notes in Control and Information Sciences, volume 151, pages 97–114. Springer Verlag, New York, 1991.
18
[30] A. Graham. Kronecker Products and Matrix Calculus with Applications. Ellis Horwood, Chichester, 1981. [31] G. Hardy, J.E. Littlewood, and G. P´olya. Inequalities: Second edition. Cambridge University Press, Cambridge, 1988.
Appendix Proof of Lemma 1 Let A[m] ∈ Rn
m ×nm
be the matrix satisfying ∂x[m] Ax = A[m] x[m] . ∂x
(27)
From the properties of the Knonecker product (see [30]) it turns out that m−1 X ∂x[m] = x[m−1−i] ⊗ In ⊗ x[i] ∂x i=0
and
¡
¢ x[m−1−i] ⊗ In ⊗ x[i] Ax = x[m−1−i] ⊗ Ax ⊗ x[i]
= (Inm−1−i ⊗ A ⊗ Ini ) x[m] .
Therefore, matrix A[m] admits the following expression A[m] =
m−1 X
Inm−1−i ⊗ A ⊗ Ini .
(28)
∂x{m} Ax = A[m] Km x{m} . ∂x
(29)
i=0
Now, from (11) and (27) one has that Km
It is not difficult to check that Km has full rank. Since Km has dimension nm × d and contains d linearly independent rows, expression (12) follows from (10) and (28), once that both sides of (29)
0 K )−1 K 0 . are multiplied by (Km m m
¥
Proof of Lemma 2 The first part is a well-known property of homogeneous forms [31]. In order to prove the second part, let us define the normalized minimum of g2m (x) as εg =
min g2m (x)
x∈Rn
subject to kxk = 1. 19
Obviously, εg > 0 since g2m (x) is positive definite. Let us introduce the homogeneous form h2m (x) = g2m (x) − εkxk2m . It turns out that h2m (x) is a nonnegative homogeneous form. Indeed, its normal-
ized minimum is equal to εh = εg − εg = 0. Hence, the first part of this lemma guarantees that
h2m (x) can be expressed as sum of squares of polynomials if (n, 2m) ∈ E, i.e. h2m (x) = =
P ¡
P
2 i i fm,i (x) = 0 P 0 {m} {m} x ( i ti ti ) x
t0i x{m}
¢2
0
= x{m} Hx{m} ,
where fm,i are homogeneous forms of degree m, ti ∈ Rd and H =
P
0 i ti ti
≥ 0. Therefore, there
exists a positive semidefinite SMR matrix for h2m (x). Let us define the following SMR matrix G
of g2m (x): G = H + εN where N is the diagonal SMR matrix of kxk2m (such N exists since the
homogeneous form kxk2m contains only monomials with even powers of the variables xi ). Then, let us observe that N ≥ Id since the coefficients of the monomials in kxk2m are greater or equal to 1.
Therefore, G ≥ H + εId > 0, that is G is positive definite.
20
¥
Dipartimento di Ingegneria dell’Informazione, Universit`a di Siena Via Roma, 56 - 53100 Siena, Italia 2
Dipartimento di Sistemi e Informatica, Universit`a di Firenze Via di S. Marta, 3 - 50139 Firenze, Italia
Abstract The problem addressed in this paper is the construction of homogeneous polynomial Lyapunov functions for linear systems with time-varying structured uncertainties. A sufficient condition for the existence of a homogeneous polynomial Lyapunov function of given degree is formulated in terms of a Linear Matrix Inequalities (LMI) feasibility problem. This condition turns out to be also necessary in some cases depending on the dimension of the system and the degree of the Lyapunov function. The maximum `∞ norm of the parametric uncertainty for which there exists a homogeneous polynomial Lyapunov function is computed by solving a generalized eigenvalue problem. The construction of such Lyapunov functions is efficiently performed by means of popular convex optimization tools for the solution of problems in LMI form. Comparisons with other classes of Lyapunov functions in numerical examples taken from the literature show that homogeneous polynomial Lyapunov functions are a powerful tool for robustness analysis.
Keywords: Robustness, Lyapunov function, Homogeneous forms, Linear Matrix Inequalities.
1
1
Introduction
Lyapunov functions are a standard tool for tackling robust analysis of linear systems affected by time-varying structured uncertainties. Quadratic Lyapunov functions have been considered by many authors, since a long time (see e.g. [1, 2, 3, 4]). More recently, relationships between quadratic stability and Integral Quadratic Constraints have been also investigated [5]. It is however widely recognized that quadratic functions lead to conservative estimates of the robust stability margin. Hence, nonquadratic Lyapunov functions have been addressed in the literature. Piecewise quadratic Lyapunov functions have been considered for both linear systems with timevarying perturbations [6] and switching linear systems [7]. Polyhedral Lyapunov functions have been introduced in [8] and successively considered by several authors. It has been shown that they are not conservative for robust analysis [9] and synthesis [10] of linear systems with time-varying structured uncertainties. The main drawback of polyhedral functions is that the computational burden required by their construction dramatically increases with the dimension of the system and the number of vertices of the polytope of matrices describing the uncertainty. Homogeneous polynomial Lyapunov functions are a viable alternative to the above classes of Lyapunov functions. The fact that this class of Lyapunov functions can improve robust stability results provided by quadratic Lyapunov functions has been recognized since long time [11]. In [12], the use of homogeneous polynomial functions to prove robust stability of linear systems with time-varying uncertainties has been considered and an approach based on the S-procedure has been proposed to enhance the robustness degree achieved by such functions. More recently, it has been shown that for a wide class of uncertain systems, including linear systems with time-varying structured uncertainties, robust stability is equivalent to the existence of a smooth Lyapunov function [13], and moreover that this Lyapunov function can always be chosen as a homogeneous polynomial of arbitrary degree [14]. This means that the class of homogeneous polynomial Lyapunov functions is not conservative for proving robust stability of linear systems with time-varying uncertainties. In this paper, the problem of constructing homogeneous polynomial Lyapunov functions is tackled by means of recently developed convex optimization techniques based on Linear Matrix Inequalities (LMI) [15, 16] and on certain properties of homogeneous forms which lead to a suitable matricial representation [17, 18]. First, a sufficient condition for the existence of a homogeneous polynomial Lyapunov function of given degree is formulated in terms of an LMI feasibility problem. The condition turns out to be also necessary in some cases, depending on the system dimension and the degree of the Lyapunov function. Then, the problem of computing the maximum ` ∞ norm of the parametric uncertainty for which a homogeneous Lyapunov function exists is addressed and a lower bound is obtained by solving a generalized eigenvalue problem [15]. Tightness of the lower bound is discussed, and it is shown to be related to the necessity of the existence condition. Another problem of interest is that of computing the Lyapunov function which guarantees the best
2
transient performance of an uncertain system: it is shown that also this problem can be easily tackled in the considered framework. Finally, several examples are presented, which demonstrate that homogeneous polynomial Lyapunov functions are a powerful tool for robustness analysis. The paper is organized as follows. Section 2 contains the problem formulation and preliminary material on homogeneous forms. The existence conditions for homogeneous polynomial Lyapunov functions are provided in Section 3, where relationships with previous work are also discussed. Maximization of the `∞ norm of the uncertainty is treated in Section 4, while the computation of the homogeneous polynomial Lyapunov function achieving the optimal transient performance is addressed in Section 5. Several numerical examples are presented in Section 6 and some concluding remarks are given in Section 7.
2
Problem formulation and preliminaries
Let us introduce the notation used throughout the paper. In : n × n identity matrix; A0 : transpose
of matrix A; A > 0 (A ≥ 0): symmetric positive definite (semidefinite) matrix A; A ⊗ B: Kronecker
product of matrices A and B; kxk: euclidean norm of vector x; kxk∞ , maxi |xi |: `∞ norm of x;
co{w1 , . . . , wr }: convex hull of vectors w i , i = 1, . . . , r. A function fm (x) is a homogeneous form of degree m in x ∈ Rn if
fm (x) =
X
ci1 ,i2 ,...,in xi11 xi22 . . . xinn ,
i1 +i2 +...+in =m
where i1 , i2 , . . . , in are nonnegative integers, and ci1 ,i2 ,...,in ∈ R are weighting coefficients. The form
fm (x) is said positive if fm (x) > 0 ∀x 6= 0, and nonnegative if fm (x) ≥ 0 ∀x.
2.1
Problem formulation
Let us consider the uncertain linear system x(t) ˙ = A(w(t))x(t), where A(w(t)) = and A0 , . . . , As ∈
Rn×n
Ã
A0 +
s X i=1
wi (t)Ai
(1) !
,
(2)
are given matrices. The uncertain time-varying parameter vector w(t) =
0
(w1 (t), . . . , ws (t)) is assumed to satisfy for all t ≥ 0 the constraint w(t) ∈ W , co{w 1 , . . . , wr }
(3)
where w i ∈ Rs , i = 1, . . . , r, are given vectors.
The basic problem addressed in this paper is the construction of a Lyapunov function proving global asymptotic stability of system (1)-(3). The attention is restricted to a special class of Lyapunov 3
functions: the homogeneous polynomial forms of degree 2m, hereafter referred to as Homogeneous Polynomial Lyapunov Functions (HPLFs) and denoted by v2m (x). More specifically, the aim is to find v2m (x) such that (i) v2m (x) > 0 for all x 6= 0; (ii) v˙ 2m (x) < 0 for all x 6= 0 and for all w(t) ∈ W. Before proceeding, let us introduce the key representation of homogeneous forms that will be exploited throughout the paper.
2.2
Square Matricial Representation of homogeneous forms
Let g2m (x) be a homogeneous form of degree 2m. Then, the Square Matricial Representation (SMR) of g2m (x) is defined as 0
g2m (x) = x{m} Cg x{m} ,
(4)
where x{m} ∈ Rd is the base vector of homogeneous forms of degree m in x and Cg = Cg0 ∈ Rd×d is
a coefficient matrix. Let Ei ∈ R(n+1−i)×n be the matrix composed by the rows i, i + 1, . . . n of In .
Then, the base vector x{m} containing all monomials of x x1 x{m−1} x2 (E2 x){m−1} x{m} = .. . {m−1} xn−1 (En−1 x) xm n
degree m is defined as if m = 1,
otherwise.
It is easy to check that the dimension d of x{m} is given by d=
(n + m − 1)! . (n − 1)!m!
(5)
An important property of the SMR is that matrix Cg is not unique. Indeed, all the matrices Cg satisfying (4) can be parameterized as Cg + L, with L ∈ L, where L is the set of matrices describing
the null form, i.e.
n o 0 L = L = L0 ∈ Rd×d : x{m} Lx{m} = 0 ∀x ∈ Rn .
(6)
1 (n + 2m − 1)! dL = d(d + 1) − . 2 (n − 1)!(2m)!
(7)
In [18] it has been shown that L is a linear space whose dimension is given by
4
Then, the family of matrices Cg describing g2m (x) can be parameterized affinely as Cg (α) = Cg + L(α),
(8)
where α ∈ RdL is a vector of free parameters and L : RdL → L is a linear parameterization of L. Hence, the Complete SMR (CSMR) of g2m (x) is defined as 0
g2m (x) = x{m} Cg (α)x{m} .
(9)
In the sequel, we will refer to Cg and Cg (α) as SMR and CSMR matrix of g2m (x), respectively.
3
Existence conditions for homogeneous Lyapunov functions
In this section, conditions for the existence of a HPLF for system (1)-(3) are provided. Such conditions are based on the SMR of homogeneous forms previously introduced and are expressed as Linear Matrix Inequalities.
3.1
Sufficient condition
For a generic system x(t) ˙ = Ax(t), let us introduce the extended matrix A {m} ∈ Rd×d , defined by d {m} ∂x{m} x = Ax = A{m} x{m} . dt ∂x
(10)
An expression of the extended matrix is provided below in terms of Kronecker products. n Lemma 1 Let x[m] = x | ⊗x⊗ {z. . . ⊗ x} , and Km ∈ R
m ×d
be the matrix satisfying
m
x[m] = Km x{m} .
(11)
Then, the extended matrix A{m} of A ∈ Rn×n is given by A{m} =
0 0 (Km Km )−1 Km
Ãm−1 X i=0
Inm−1−i ⊗ A ⊗ Ini
!
Km .
Proof. See the Appendix. The following useful properties of the extended matrix A{m} hold [11, 12]: (I) Let A, B ∈ Rn×n and α, β ∈ R. Then (αA + βB){m} = αA{m} + βB{m} .
5
(12)
(II) Let A0,{m} and Ai,{m} denote respectively the extended matrices of A0 and Ai , according to (10). Consider the extended system ! Ã s X wi (t)Ai,{m} x{m} (t) = A(w(t)){m} x{m} (t). x˙ {m} (t) = A0,{m} + i=1
0
Then, for any v2m (x) = x{m} V x{m} , one has ¯ ¯ ¯ ¯ d d ¯ v2m (x)¯ = v2m (x)¯¯ . dt dt x˙ {m} =A(w(t)){m} x{m} x=A(w(t)) ˙ x
The following result provides a sufficient condition for establishing the existence of a HPLF of degree 2m for system (1)-(3). Theorem 1 Let A¯j = A(wj ), j = 1, . . . , r and let A¯j,{m} denote the extended matrix of A¯j . If the system of LMIs
(
V > 0, −V A¯j,{m} − A¯0
j,{m} V
− L(αj ) > 0,
(13)
j = 1, . . . , r 0
admits a feasible solution V = V 0 ∈ Rd×d , αj ∈ RdL , j = 1, . . . , r, then v2m (x) = x{m} V x{m} is a
HPLF for (1)-(3).
0
Proof. If (13) admits solution, then v2m (x) = x{m} V x{m} is positive definite. By differentiating v2m (x) along the trajectories of x˙ {m} = A¯j,{m} x{m} , one gets ¯ ¯ d v2m (x)¯¯ dt ¯j,{m} x{m} x˙ {m} =A
= x{m}
0
= x{m}
0
³
³
´ V A¯j,{m} + A¯0j,{m} V x{m}
´ V A¯j,{m} + A¯0j,{m} V + L(αj ) x{m} ,
(14)
which is negative definite by (13). Observe that the CSMR (8)-(9) has been used in (14). Now, ´ ³ P notice that, due to property (I), matrices A¯j,{m} = A0,{m} + si=1 wij Ai,{m} are the vertices of P the set of extended matrices A0,{m} + si=1 wi (t)Ai,{m} , w(t) ∈ W. Convexity of this set of matrices
implies that v2m (x) is a common HPLF for all matrices of the set. Hence, property (II) allows one to conclude that v2m (x) is also a HPLF for system (1).
¥
Remark 1. It is worth observing that (13) is an LMI feasibility problem [15], which can be solved by efficient computational tools based on convex optimization [16, 19]. Notice that the free variables in (13) are matrix V = V 0 ∈ Rd×d and vectors αj ∈ RdL , j = 1, . . . , r, and therefore their number
is equal to d(d + 1)/2 + rdL − 1, with d and dL given by (5) and (7), respectively (the term −1 is due to the fact that V can be properly scaled, for example by setting to 1 the first entry of the
matrix).
6
3.2
Necessary and sufficient condition
In some cases, depending on the dimension n of x and the degree 2m of v2m (x), the sufficient condition provided by Theorem 1 is also necessary for the existence of a HPLF for (1)-(3). The next result is central for establishing such a necessary and sufficient condition. Lemma 2 Consider the equivalence set E, defined as E = {(n, 2) ∀n} ∪ {(2, 2m) ∀m} ∪ {(3, 4)} .
(15)
Let x ∈ Rn and g2m (x) be a nonnegative homogeneous form of degree 2m. If (n, 2m) ∈ E, then g2m (x) can be written as the sum of squares of homogeneous forms of degree m. Moreover, if g2m (x) is positive, then it admits a positive definite SMR, i.e. there exists a positive definite matrix 0
V ∈ Rd×d such that g2m (x) = x{m} V x{m} . Proof. See the Appendix. Theorem 2 Let (n, 2m) belong to E in (15). Then, there exists a HPLF of degree 2m for system (1) if and only if the set of LMIs (13) admits a feasible solution.
Proof. Obviously, it must be proven only that if (n, 2m) ∈ E and there exists a HPLF v 2m (x)
for (1), then (13) admits a feasible solution. By assumption we have that v2m (x) is positive definite, and −v˙ 2m (x) is positive definite for any uncertainty w(t) ∈ W. Hence, due to Lemma 2, 0
there exists V > 0 such that v2m (x) = x{m} V x{m} , and there exist Vj > 0, j = 1, . . . , r, such that the time derivative of v2m (x) evaluated along the trajectories of x˙ {m} = A¯j,{m} x{m} satisfies 0
−v˙ 2m (x) = x{m} Vj x{m} . Since the parameterization of the SMR matrices given by (6) and (8) is complete, there exists αj such that Vj = −V A¯j,{m} − A¯0j,{m} V − L(αj ). Therefore, (13) admits a
¥
feasible solution.
Remark 2. It is worth observing that, though it is well known that there exist homogenous forms that cannot be written as sums of squares, experience shows that occurrences of such forms are quite rare in practice (see e.g. [18, 20]).
3.3
Relationships with previous work on HPLF
The use of HPLFs for robust analysis of system (1)-(3) has been addresses in [12]. In this paper, a result based on the S-procedure is exploited to provide a sufficient condition for the existence of a HPLF of degree 2m. Although the condition is formulated as the minimization of a nondifferentiable convex function, it can be easily rewritten as an LMI feasibility problem of the form V > 0, d−n X ¯j,{m} − A¯0 −V A V − βij Fi > 0, j = 1, . . . , r j,{m} i=1
7
(16)
where Fi ∈ Rd×d , i = 1, . . . , d − n are suitable matrices. The difference between (16) and (13) is
that the number of free variables βij for each LMI in (16) is dZ = d − n, which turns out to be much
smaller than dL , the dimension of vector αj in (13), as clarified by Table 1. The larger number n=2 m=2
3
4
5
dL = 1
dZ = 1
6
3
20
6
50
10
3
3
2
27
7
126
16
420
30
4
6
3
75
12
465
31
1990
65
5
10
4
165
18
1310
52
7000
121
Table 1: Values of dL and dZ for some n and m. of free variables in the LMIs makes the sufficient condition provided by Theorem 1 much more powerful. Moreover, the fact that the SMR of homogeneous forms is complete, allows to formulate the necessary and sufficient condition of Theorem 2 for the cases in which (n, 2m) ∈ E. On the contrary, the condition provided in [12] is necessary only for the case n = m = 2, in which the two
parameterization of homogeneous forms coincide (see Table 1).
4
Computation of the `∞ 2m-HPLF stability margin
In the analysis of uncertain systems of type (1)-(2), a key problem is that of finding the largest value of the positive scalar γ for which there exists a HPLF of degree 2m for all w(t) belonging to the scaled perturbation set γW. In order to simplify the presentation, we restrict our attention to the case in which the perturbation set is the `∞ box, i.e. Bγ = {q ∈ Rs : |qi | ≤ γ, i = 1, . . . , s} .
(17)
s
Clearly, Bγ is equal to γ co{u1 , . . . , u2 }, where uj , j = 1, . . . , 2s , are the vertices of the unit `∞ ball B1 . Said another way, the aim is to compute
∗ γ2m = sup{γ : ∃ v2m (x) for (1)-(2), with w(t) ∈ Bγ }.
(18)
∗ will be referred to as the ` In the following, γ2m ∞ 2m-HPLF stability margin for system (1)-(2).
We also consider the related problem of computing the following quantity κ∗2m = sup{κ : ∃ v2m (x) for (1)-(2), with w(t) ∈ B¯κ },
(19)
B¯κ = {q ∈ Rs : 0 ≤ qi ≤ κ, i = 1, . . . , s} .
(20)
where
∗ because the perturbation set is restricted to the positive orthant Notice that κ∗2m differs from γ2m
of the parameter space. For this reason, it will be referred to as the ` ∞ 2m-HPLF positive stability 8
∗ and κ∗ are bounded from above margin. It is worth recalling that the above stability margins γ2m 2m
by the well-known `∞ state space parametric stability margin (see for example [21, 22]). Let us first focus the attention on problem (18)-(17). Define the matrices A˜j = A(uj ) − A0 , where uj , j = 1, . . . , 2s are the vertices of the unit `∞ ball B1 . Moreover, let A˜j,{m} denote the extended ∗ matrix of A˜j . The following result shows that a lower bound γˆ2m of the `∞ 2m-HPLF stability ∗ can be computed by solving a quasiconvex optimization problem. margin γ2m ∗ be defined as Theorem 3 Let γˆ2m ∗ )−1 = (ˆ γ2m
inf z ∈ R; V = V 0 ∈ Rd×d ;
z
αj ∈ RdL , j = 0, . . . , 2s V > 0, z > 0 −V A0,{m} − A0,{m} 0 V − L(α0 ) > 0 s.t. ´ ³ z −V A 0 0 j ˜ ˜0 0,{m} − A0,{m} V − L(α ) > V Aj,{m} + Aj,{m} V + L(α )
(21)
j = 1, . . . , 2s .
∗ ≤ γ ∗ . Moreover, if (n, 2m) belongs to E, then γ ∗ = γ∗ . Then, γˆ2m ˆ2m 2m 2m 0
Proof. Let v2m (x) = x{m} V x{m} . The time derivative of v2m (x) evaluated for w(t) = z −1 uj , is given by ¯ · ³ ´ ³ ´0 ¸ ¯ d 0 v2m (x)¯¯ = x{m} V A0,{m} + z −1 A˜j,{m} + A0,{m} + z −1 A˜j,{m} V x{m} dt w(t)=z −1 uj
h ³ ´ i 0 = z −1 x{m} z V A0,{m} + A00,{m} V + L(α0 ) + V A˜j,{m} + A˜0j,{m} V + L(αj ) x{m} ,
where the CSMR has been exploited by adding the term zL(α0 ) + L(αj ) inside the square brackets. Hence, the constraint in (21) guarantees that v˙ 2m (x) is negative definite. Since this holds for all ∗ ≤ γ∗ . w(t) such that kw(t)k∞ ≤ z −1 , one has that γˆ2m 2m
∗ , Let us assume now that (n, 2m) belongs to E. From definition (18), for any z such that z −1 < γ2m
there exists a positive definite v2m (x) such that −v˙ 2m (x) is positive definite for w(t) : kw(t)k∞ ≤ z −1 . From Lemma 2, it follows that −v˙ 2m (x) admits a positive definite SMR matrix for any w(t) ∈ Bz −1 and therefore the constraint in problem (21) admits a feasible solution for any
∗ . Hence, γ ∗ = γ∗ . z −1 < γ2m ˆ2m 2m
¥
Remark 3. Problem (21) is a Generalized Eigenvalue Problem (GEVP), which has been proven to be a quasiconvex optimization problem and can be tackled by efficient optimization tools [15, 16]. Notice that the LMI − V A0,{m} − A0,{m} 0 V − L(α0 ) > 0 in (21) guarantees that the matrix multi-
plying the generalized eigenvalue z is positive, as required by the standard form of the GEVP (see
9
[15]). Observe that the number of free variables in the GEVP (21) is equal to d(d+1)/2+(2 s +1)dL . A result similar to Theorem 3 can be obtained for the `∞ 2m-HPLF positive stability margin κ∗2m in (19)-(20). Here, we report the result for the simpler case of a segment of matrices, because it has been widely addressed in the literature and it will be useful in the examples presented in the next section. The aim is to compute the largest value of κ such that there exist a HPLF of degree 2m for the matrices A0 + w(t)A1 , with 0 ≤ w(t) ≤ κ. Corollary 1 Let s = 1 in (19)-(20). Define κ ˆ ∗2m as (ˆ κ∗2m )−1 =
z min z ∈ R; V = V 0 ∈ Rd×d ;
α 1 , α 2 ∈ R dL V >0 −V A0,{m} − A00,{m} V − L(α1 ) > 0 s.t. ³ ´ z −V A 0 1) > V A 0 2 − A V − L(α 0,{m} 1,{m} + A1,{m} V + L(α ) 0,{m}
(22)
ˆ ∗2m = κ∗2m . Then, κ ˆ ∗2m ≤ κ∗2m . Moreover, if (n, 2m) belongs to E, then κ
5
Construction of the optimal performance HPLF
Another problem of interest in the robustness analysis of uncertain systems is that of determining the Lyapunov function that achieves the best transient performance (see, for example, the classical textbook [23]). For a given Lyapunov function v(x), one can define the transient performance index of system (1)-(3) as λ(v) = sup
sup
x∈Rn \0 w(t)∈W
v(x) ˙ . v(x)
(23)
From (23) one has v(x(t)) ≤ v(x(t0 ))eλ(v)(t−t0 ) , thus establishing the rate of decrease of v.
Therefore, it is natural to select among all feasible Lyapunov functions, the one that minimizes λ(v). For systems of the type (1)-(3), this problem has been addressed in [24] within the class of quadratic Lyapunov functions. Here, the aim is to select the optimal Lyapunov function among HPLFs of degree 2m, i.e. to compute λ∗2m = inf λ(v2m ). v2m
The following result shows that an upper bound of λ∗2m can be obtained by solving a GEVP.
10
(24)
ˆ ∗ be defined as Theorem 4 Let λ 2m ˆ∗ = λ 2m
z inf 0 d×d z ∈ R; V = V ∈ R ;
αj ∈ RdL , j = 1, . . . , r ( V > 0, s.t. zV > V A¯j,{m} + A¯0j,{m} V + L(αj )
(25) j = 1, . . . , r.
ˆ ∗ ≥ λ∗ . Moreover, if (n, 2m) belongs to E, then λ ˆ ∗ = λ∗ . Then, λ 2m 2m 2m 2m Proof. Let us assume that there exists z ∈ R and V = V 0 ∈ Rd×d satisfying the constraint in (25), 0
for some vectors αj ∈ RdL . Then, by setting v2m (x) = x{m} V x{m} and differentiating along the trajectories of x˙ {m} = A¯j,{m} x{m} , one has ³ ´ ¯ {m} 0 V A 0 ¯ ¯ x + A V x{m} j,{m} j,{m} v˙ 2m (x) ¯¯ = = {m} 0 V x{m} v2m (x) ¯x˙ {m} =A¯j,{m} x{m} x ³ ´ 0 x{m} V A¯j,{m} + A¯0j,{m} V + L(αj ) x{m} ≤z = 0 x{m} V x{m} for all j = 1, . . . , r. Therefore, by exploiting convexity of W, one can conclude that v˙ 2m /v2m < z
for all w(t) ∈ W and hence λ∗2m ≤ z. Minimizing with respect to z, one obtains the upper bound ˆ∗ . λ 2m Now, let (n, 2m) ∈ E. Definition (24) says that for any z > λ∗2m there exists a positive definite v2m (x) such that v˙ 2m (x)/v2m (x) < z, for all x and for all w(t) ∈ W. Then, Lemma 2 implies
that zv2m (x) − v˙ 2m (x) admits a positive definite SMR matrix for any w(t) ∈ W, and therefore the
constraint in problem (25) admits a feasible solution for any z > λ∗2m .
6
¥
Examples
The following examples illustrate the application of the methodology for the construction of HPLFs described in the previous section.
6.1
Example 1
This example is taken from [12]. Let us consider the matrices à ! à ! 0 1 0 0 A0 = , A1 = −2 −1 −1 0 and assume we want to compute the `∞ 2m-HPLF positive stability margin κ∗2m , defined in (19)(20). Since in this example n = 2, the lower bound provided by the GEVP (22) is tight, i.e. 11
κ ˆ ∗2m = κ∗2m , as assured by Corollary 1. In [12], it has been observed that a quadratic Lyapunov function exists only for κ < 3.82, and that this bound can be improved to κ < 5.73 by solving problem (16) with m = 2. Table 2 shows the values of κ∗2m obtained from (22), for some different m. m
2
3
4
5
6
7
8
κ∗2m
5.7393
6.2135
6.3982
6.6469
6.6581
6.7801
6.7962
Table 2: Values of κ∗2m for different m in Example 1. For the case m = 3, vector x{3} and extended matrices A0,{3} and A1,{3} are given by
x
{3}
x31
0
3
0
0
x2 x −2 −1 2 0 1 2 = , A0,{3} = x1 x22 0 −4 −2 1 3 x2 0 0 −6 −3
Matrix L(α) is equal to
0
0
0 2α1 L(α) = −α1 α2 −α2 −α3
0 0
α2
2α3
−α3 . 0 0
The matrix V obtained by solving (22) for m = 3 turns out to be 1.0000 0.1524 −0.1965 −0.0174 0.1524 1.2187 0.2086 −0.0280 V = −0.1965 0.2086 0.2938 0.0231 −0.0174 −0.0280 0.0231 0.0099
corresponding to the HPLF
0
−1 0 0 0 , A1,{3} = . 0 −2 0 0 0 0 −3 0
−α1 −α2
0
0
v6 (x) = x61 + 0.3048x51 x2 + 0.8258x41 x22 + 0.3825x31 x32 + 0.2378x21 x42 + 0.0463x1 x52 + 0.0099x62 . The same example has been considered also in [25] and [6]. In [25], a polyhedral Lyapunov function has been constructed, guaranteeing asymptotic stability for κ = 6. In [6], a piecewise quadratic Lyapunov function achieving stability for κ = 6.2, has been obtained via a sequence of LMI optimizations and a grid search over two free parameters. As it can be seen from Table 2, these bounds are improved by all the HPLFs of degree greater than 2. By the way, notice that the HPLF is obtained by solving one single GEVP, with 12 (3m2 + m + 2) free variables.
12
6.2
Example 2
This example shows that the gap in the Lyapunov stability margin between our technique and the one presented in [12] can be very large. Let us consider the matrices 0 1 0 −2 0 −1 , A1 = 1 −10 A0 = 0 0 1 3 −1 −2 −4 3 −4 2
(26)
and assume we want to solve the same problem as in Example 1. Quadratic stability is guaranteed only for κ ≤ 1.9042. Let us compute a HPLF of degree 4. Being n = 3 and m = 2 it follows that d = 6 and dL = 6. Then, vector x{2} and matrix L(α) are 0 0 x21 0 x1 x2 2α1 x1 x3 α2 , L(α) = 0 x{2} = 2 0 −α1 x2 −α2 −α4 x2 x3 2 −α3 −α5 x3
respectively given by: 0 α2
−α1 −α2 −α3 0
2α3
α4
α4
0
α5
0
0
−α6
−α4 −α5 α5 0 . 0 −α6 2α6 0 0 0
The GEVP (22) returns the lower bound κ ˆ ∗4 = 75.1071. Moreover, since (n, 2m) = (3, 4) ∈ E, we
have from Corollary 1 that κ ˆ ∗4 = κ∗4 .
Using the approach proposed in [12], one finds that the maximum κ for which robust stability is guaranteed is equal to 17.8347. The remarkable difference with respect to our approach is due to the fact that the parameterization of homogeneous forms adopted in [12] is not complete, as discussed in Subsection 3.3. Specifically, it is easy to see that in [12] only the parameters α 1 , α3 , α6 in the matrix L(α) are considered (see also Table 1).
6.3
Example 3
Let us consider the differential equation ˙ + k(t)ξ(t) = 0 ¨ + ξ(t) ξ(t) and assume that we want to compute the maximum κ such that the solution remain bounded, for all 0 ≤ k(t) ≤ κ. By adopting a strategy similar to that proposed in [26], one can determine numerically
the maximum value of κ. In particular, one has that such value must satisfy the equation ½ ¾ √ √ £ ¤ 1 k exp − √ π − arctan( 4κ − 1) − 1 = 0 4κ − 1 and is equal to κ∗ ≈ 3.0448.
13
Clearly, the problem can be easily tackled in the framework of HPLF by solving (19)-(20) with s = 1 and A0 =
Ã
0
1
0 −1
!
,
A1 =
Ã
0 0 −1 0
!
.
The values of κ∗2m obtained from (22) are reported in Table 3. Notice that the values exhibit a growth m
1
2
3
4
5
6
κ∗2m
1.0000
1.5000
1.9997
2.2952
2.4009
2.5053
m
7
8
9
10
11
12
κ∗2m
2.6102
2.6601
2.6954
2.7466
2.7726
2.7917
Table 3: Values of κ∗2m for different m in Example 3. towards the theoretical upper bound κ∗ . On the other hand, the results in [13, 14] guarantee that for any κ < κ∗ there exists a HPLF of suitable degree, hence it is expected that κ∗2m converges to κ∗ as m approaches infinity. Another interesting byproduct of the treatment in [26] is that one can calculate the worst-case sequence k(t) (which consists of suitable switchings between k(t) = 0 and k(t) = κ ∗ ) and the corresponding trajectory of the system, that represents the limit curve towards which the level surfaces of the Lyapunov functions are expected to tend. This is confirmed by Figure 1a, where the level curves of the obtained HPLFs for different values of m are depicted. Figure 1b clearly shows that the HPLF corresponding to m = 12 is very close to the limit trajectory predicted by the theory.
6.4
Example 4
This example has been considered by several authors [27, 24, 10]. Let us consider system (1)-(2) with s = 1 and matrices A0 =
Ã
0
1
−1 −1
!
,
A1 =
Ã
0 0 1 0
!
.
The maximum γ for which there exists a Lyapunov function for system (1)-(2), with w(t) ∈ B γ , has
been calculated for different classes of Lyapunov functions. Notice that the ` ∞ state space stability
∗ ≤ 1, for all margin is given by ρ = 1, because A0 + A1 is not asymptotically stable, and hence γ2m √ m. In [27], it has been proven that quadratic Lyapunov functions can achieve γ2∗ = 3/2. In [10], a
polyhedral Lyapunov function has been provided, which guarantees stability for |w(t)| ≤ γ = 0.98.
The level curves of this Lyapunov function have 30 vertices.
By applying Theorem 3 for different values of m, one can construct HPLFs of increasing degree, trying to improve the value of the 2m-HPLF stability margin. Observe that also in this example ∗ because n = 2. For m = 2, one obtains the lower bound provided by (21) always coincides with γ2m
14
1.5
1.5
1
1
0.5
0.5
x2
x2
limit trajectory 0
0
−0.5
−0.5
m=12
m=4 m=3 m=2 −1
−1.5 −1.5
−1
m=1
−1
−0.5
0
x
0.5
1
−1.5 −1.5
1.5
−1
−0.5
0
x
1
(a) level curves of v2m (x) for different m
0.5
1
1.5
1
(b) level curve of v24 (x) and limit trajectory (dashed)
Figure 1: Level curves of HPLFs in Example 3. ∗ = 1, which means that it is γ4∗ = 0.9771, but for m ≥ 3 the result returned by the GEVP is γˆ2m
possible to construct a HPLF of degree 6 (or more) for |w(t)| ≤ γ and γ arbitrarily close to 1.
In particular, for m = 3 the following matrix V has 1.0000 2.9987 2.9987 14.1343 V = 1.9974 12.4475 0.2491 2.1782 corresponding to γ6∗ = 1 and the HPLF
been obtained
1.9974 0.2491
12.4475 2.1782 14.5114 3.5064 3.5064 1.3106
v2m (x) = x61 + 5.9974x51 x2 + 18.1291x41 x22 + 25.3932x31 x32 + 18.8678x21 x42 + 7.0128x1 x52 + 1.3106x62 .
6.5
Example 5
This example concerns the computation of an optimal performance HPLF, in the sense explained in Section 5. Let us consider the helicopter model originally proposed in [28], and the robust controller designed in [29]. The resulting closed-loop uncertain system is given by −0.1027
0.1762
0
0
0.1995
−0.3446
−0.7275 − 0.1640q (t) −1.4280 + 0.2699q (t) 2.1852 + 0.4511q (t) 0.7218 + 0.4308q (t) 3 3 3 3 x(t) ˙ = 1.1689 −0.4446 + q1 (t) −3.6757 −3.0841 + q2 (t) 1
0
x(t)
with |q1 (t)| ≤ 0.2192, |q2 (t)| ≤ 1.2031, |q3 (t)| ≤ 2.0673. Clearly, the system can be easily written in
form (1)-(3) with r = 8.
15
In [24], the problem of computing the quadratic Lyapunov function achieving the best transient performance, defined as in (23), has been addressed. The value obtained was λ = −0.3839. By solving the GEVP problem (25) with m = 2, one obtains the HPLF
v4 (x) = x41 + 0.2184x31 x2 − 2.8083x31 x3 − 4.0483x31 x4 + 1.1458x21 x22 + 0.1831x21 x2 x3
+0.3233x21 x2 x4 + 9.7670x21 x23 + 16.4758x21 x3 x4 + 9.7066x21 x24 + 0.0283x1 x32
−2.5961x1 x22 x3 − 3.9673x1 x22 x4 − 0.1303x1 x2 x23 − 3.6305x1 x2 x3 x4 − 3.3848x1 x2 x24
−14.8574x1 x33 − 37.3189x1 x23 x4 − 31.4426x1 x3 x24 − 11.2239x1 x34 + 0.1064x42
+0.1657x32 x3 + 0.2465x32 x4 + 2.1908x22 x23 + 5.1952x22 x3 x4 + 4.2363x22 x24 −0.0017x2 x33 + 3.0142x2 x23 x4 + 6.4605x2 x3 x24 + 3.4452x2 x34 + 8.8681x43
+29.0364x33 x4 + 36.5118x23 x24 + 20.6917x3 x34 + 5.8479x44
ˆ ∗ = −0.8889. Therefore, once again it can be seen that the HPLFs outperform the achieving λ 4
classic approaches based on quadratic Lyapunov functions. Moreover, this example shows that the proposed LMI-based procedures are able to handle also quite complex systems (in this case, a fourth order uncertain system with 8 vertices).
7
Conclusions
The construction of homogeneous polynomial Lyapunov functions for linear systems with timevarying structured uncertainties has been addressed via LMI optimization techniques. Maximization of the size of the parametric uncertainty can also be efficiently tackled in the proposed framework. With respect to previous work on this class of Lyapunov functions, better results have been obtained by exploiting a complete parameterization of homogeneous forms of given degree. Moreover, this allows one to formulate necessary conditions for the existence of a homogeneous polynomial Lyapunov function in some cases (for example, for all second order systems). Comparisons with other classes of Lyapunov functions are also very promising. Quadratic Lyapunov functions are often outperformed by HPLFs. On the other hand, polyhedral Lyapunov functions seem to require level surfaces with a very high number of vertices, even in cases that can be handled by HPLFs of quite low degree (and hence with a limited number of coefficients). Further investigations on the potentialities of homogeneous polynomial Lyapunov functions will be the subject of future research.
References ˘ [1] D. D. Siljak. Nonlinear Systems: Parametric Analysis and Design. John Wiley & Sons, New York, 1969.
16
[2] K. S. Narendra and J. H. Taylor. Frequency Domain Criteria for Absolute Stability. Academic, New York, 1973. [3] E. S. Piatnitskiy and V. I. Skorodinskiy. Numerical methods of Lyapunov functions construction and their application to the absolute stability problem. Systems and Control Letters, 2:130–135, 1982. [4] K. Zhou and P. P. Khargonekar. Stability robustness bounds for linear state-space models with structured uncertainties. IEEE Transactions on Automatic Control, 32:621–623, 1987. [5] A. Megretski and A. Rantzer. System analysis via integral quadratic constraints. IEEE Transactions on Automatic Control, 42:819–830, 1997. [6] L. Xie, S. Shishkin, and M. Fu. Piecewise Lyapunov functions for robust stability of linear time-varying systems. Systems and Control Letters, 31:165–171, 1997. [7] M. Johansson and A. Rantzer. Computation of piecewise quadratic Lyapunov functions for hybrid systems. IEEE Transactions on Automatic Control, 43:555–559, 1998. [8] R. K. Brayton and C. H. Tong. Stability of dynamical systems: a constructive approach. IEEE Transactions on Circuits and Systems, 26:224–234, 1979. [9] A. P. Molchanov and E. S. Piatnitskiy. Lyapunov functions specifying necessary and sufficient conditions of absolute stability of nonlinear nonstationary control system. Automat. Remote Contr., 47:344–354 (part I), 443–451 (II), 620–630 (III), 1986. [10] F. Blanchini. Nonquadratic Lyapunov functions for robust control. Automatica, 31:451–461, 1995. [11] R. W. Brockett. Lie algebra and Lie groups in control theory. In D.Q. Mayne and R.W. Brockett, editors, Geometric Methods in Systems Theory, pages 43–82. Dordrecht, Reidel, 1973. [12] A. L. Zelentsovky. Nonquadratic Lyapunov functions for robust stability analysis of linear uncertain systems. IEEE Transactions on Automatic Control, 39:135–138, 1994. [13] W. P. Dayawansa and C. F. Martin. A converse Lyapunov theorem for a class of dynamical systems which undergo switching. IEEE Transactions on Automatic Control, 44:751–760, 1999. [14] D. Angeli. LMI conditions for the determination of polynomial Lyapunov functions. Technical report, DSI - University of Florence, Italy, 2001. [15] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan. Linear Matrix Inequalities in System and Control Theory. SIAM, Philadelphia, 1994. 17
[16] Y. Nesterov and A. Nemirovsky. Interior Point Polynomial Methods in Convex Programming: Theory and Applications. SIAM, Philadelphia, 1993. [17] G. Chesi, A. Tesi, A. Vicino, and R. Genesio. An LMI approach to constrained optimization with homogeneous forms. Systems and Control Letters, 42:11–19, 2001. [18] G. Chesi, A. Garulli, A. Tesi, and A. Vicino. LMI-based techniques for solving quadratic distance problems. In Proc. of 40-th IEEE Conf. on Decision and Control, Orlando, December 2001. [19] P. Gahinet, A. Nemirovski, A. J. Laub, and M. Chilali. LMI Control Toolbox. The Mathworks Inc., 1995. [20] P. Parrilo. Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization. Ph.D. Thesis, California Institute of Technology, 2000. [21] S. P. Bhattacharyya, H. Chapellat, and L. H. Keel. Robust Control: The Parametric Approach. Prentice Hall, NJ, 1995. [22] A. Vicino, A. Tesi, and M. Milanese. Computation of non conservative stability perturbation bounds for systems with nonlinearly correlated uncertainties. IEEE Transactions on Automatic Control, 35:835–841, 1990. [23] J. C. Hsu and A. U. Meyer. Modern Control Principles and Applications. McGraw-Hill Book Company, New York, 1968. [24] A. Olas. Construction of optimal Lyapunov function for systems with structured uncertainties. IEEE Transactions on Automatic Control, 39:167–171, 1994. [25] F. Blanchini and S. Miani. On the transient estimate for linear systems with time-varying uncertain parameters. IEEE Transactions on Circuits and Systems, Part I, 43:592–596, 1996. [26] R. W. Brockett. Finite Dimensional Linear Systems. John Wiley and Sons, NY, 1970. [27] B. Radziszewski. O Najlepszej Funkcji Lapunowa. IFTR Reports, 26, 1977. [28] K. S. Narendra and S. S. Tripathi. Identification and optimization of aircraft dynamics. Journal of Aircraft, 10:193–199, 1973. [29] Y. H. Chen and J. S. Chen. Robust control of uncertain systems with time-varying uncertainty: a computer-aided setting. In Lecture Notes in Control and Information Sciences, volume 151, pages 97–114. Springer Verlag, New York, 1991.
18
[30] A. Graham. Kronecker Products and Matrix Calculus with Applications. Ellis Horwood, Chichester, 1981. [31] G. Hardy, J.E. Littlewood, and G. P´olya. Inequalities: Second edition. Cambridge University Press, Cambridge, 1988.
Appendix Proof of Lemma 1 Let A[m] ∈ Rn
m ×nm
be the matrix satisfying ∂x[m] Ax = A[m] x[m] . ∂x
(27)
From the properties of the Knonecker product (see [30]) it turns out that m−1 X ∂x[m] = x[m−1−i] ⊗ In ⊗ x[i] ∂x i=0
and
¡
¢ x[m−1−i] ⊗ In ⊗ x[i] Ax = x[m−1−i] ⊗ Ax ⊗ x[i]
= (Inm−1−i ⊗ A ⊗ Ini ) x[m] .
Therefore, matrix A[m] admits the following expression A[m] =
m−1 X
Inm−1−i ⊗ A ⊗ Ini .
(28)
∂x{m} Ax = A[m] Km x{m} . ∂x
(29)
i=0
Now, from (11) and (27) one has that Km
It is not difficult to check that Km has full rank. Since Km has dimension nm × d and contains d linearly independent rows, expression (12) follows from (10) and (28), once that both sides of (29)
0 K )−1 K 0 . are multiplied by (Km m m
¥
Proof of Lemma 2 The first part is a well-known property of homogeneous forms [31]. In order to prove the second part, let us define the normalized minimum of g2m (x) as εg =
min g2m (x)
x∈Rn
subject to kxk = 1. 19
Obviously, εg > 0 since g2m (x) is positive definite. Let us introduce the homogeneous form h2m (x) = g2m (x) − εkxk2m . It turns out that h2m (x) is a nonnegative homogeneous form. Indeed, its normal-
ized minimum is equal to εh = εg − εg = 0. Hence, the first part of this lemma guarantees that
h2m (x) can be expressed as sum of squares of polynomials if (n, 2m) ∈ E, i.e. h2m (x) = =
P ¡
P
2 i i fm,i (x) = 0 P 0 {m} {m} x ( i ti ti ) x
t0i x{m}
¢2
0
= x{m} Hx{m} ,
where fm,i are homogeneous forms of degree m, ti ∈ Rd and H =
P
0 i ti ti
≥ 0. Therefore, there
exists a positive semidefinite SMR matrix for h2m (x). Let us define the following SMR matrix G
of g2m (x): G = H + εN where N is the diagonal SMR matrix of kxk2m (such N exists since the
homogeneous form kxk2m contains only monomials with even powers of the variables xi ). Then, let us observe that N ≥ Id since the coefficients of the monomials in kxk2m are greater or equal to 1.
Therefore, G ≥ H + εId > 0, that is G is positive definite.
20
¥
