New Approach to Robust-Stability Analysis of Linear Time-Invariant ...
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
New Approach to Robust -Stability Analysis of Linear Time-Invariant Systems With Polytope-Bounded Uncertainty Eduardo N. Gonçalves, Reinaldo M. Palhares, Ricardo H. C. Takahashi, and Renato C. Mesquita
Abstract—This note presents a new approach to robust -stability analysis of linear time-invariant systems with polytope-bounded uncertainty. The proposed approach combines sufficient conditions for robust -stability in terms of feasibility problems with linear matrix inequalities (LMI) constraints and a new polytope partition technique. If the initial polytope does not attain the robust -stability sufficient condition, the polytope is successively subdivided until all subpolytopes attain the sufficient condition, in the case of robustly -stable uncertain system, or it is found a subpolytope vertex that does not attain the regional pole-placement constraints, in the case of an uncertain system that is not robustly -stable. It is presented a new general format polytope partition technique that allows the implementation of the proposed approach. The efficiency of the proposed approach is verified by means of illustrative examples and three different LMI-based analysis formulations. Index Terms—Polytope-bounded uncertainty, robust stability, simplex subdivision.
1709
D-stable. This branch-and-bound like strategy has been already considered for robust stability analysis, e.g., [7], [8]. The partition technique considered in these previous works is based on the bisection of hyperrectangles. Because of this, these works are restrict to the case of systems with uncertain parameters that range between extremal values. The major contribution of this work is not only to extend this strategy to robust D -stability analysis, but it is how this strategy is implemented based on up-to-date LMI-based analysis formulations and a new polytope partition technique based on simplicial meshes. The advantage of considering simplicial meshes is the capability for dealing with polytopes with any shape, not restricted to the hyperrectangle case as in [7], [8]. This feature allows the proposed analysis procedure to be applied to both affine parameter-dependent and polytopic models, what is not possible with the aforementioned works. A complete discussion about advantages of simplicial meshes over box meshes and other mesh generation approaches can be found in [9]. Without the conservatism of LMI-based sufficient conditions, the proposed strategy allows the stability analysis of closed-loop robust systems designed through approaches not based only on LMI framework, as illustrated in the examples.
-stability, robust
NOTATION AND DEFINITIONS
I. INTRODUCTION The Lyapunov theory has been frequently employed to robust stability analysis due to the fact that the resulting linear matrix inequality (LMI)-based feasibility problem is easy to solve by one of the available LMI-solvers. The quadratic stability condition based on a single Lyapunov function is the simplest formulation but it can be too conservative in the case of time-invariant systems. To reduce the conservativeness, one can use parameter-dependent Lyapunov functions and slack variables [1]–[3]. Unfortunately, these LMI-based analysis formulations are only sufficient conditions. As verified in [4], the reliability of the LMI-based analysis formulations decreases with the increasing of the number of polytope vertices and the order of the system. It is possible to formulate a more reliable formulation increasing the number of decision variables but at the cost of a higher computational cost. The Lyapunov theory can be extended to deal with the problem of robust D-stability to guarantee that all closed-loop poles are robustly placed at convex regions of the complex plane, denoted LMI regions [5], [6]. In this note, it is presented an analysis approach that finds out if a linear time-invariant system is robustly D -stable or not. To accomplish this task, the proposed analysis approach is based on a procedure that combines a D -stability sufficient condition and a polytope partition technique. The proposed procedure subdivides the polytope iteratively until it is verified that all generated subpolytopes are robustly D-stable or it is found a system instance of the polytope that is not Manuscript received June 3, 2005; revised December 7, 2005 and March 29, 2006. Recommended by Associate Editor Y. Ohta. This work was supported in part by the Brazilian agencies CNPq and FAPEMIG. E. N. Gonçalves is with the Department of Electrical Engineering, Federal Center of Technological Education of Minas Gerais, 31510-470 Belo Horizonte, MG, Brazil (e-mail: [email protected]). R. M. Palhares and R. C. Mesquita are with the Department of Electronics Engineering, Federal University of Minas Gerais, 31270-010 Belo Horizonte, MG, Brazil (e-mail: [email protected]; [email protected]). R. H. C. Takahashi is with the Department of Mathematics, Federal University of Minas Gerais, 31270-010 Belo Horizonte, MG, Brazil (e-mail: [email protected]). Digital Object Identifier 10.1109/TAC.2006.883061
The transpose of a matrix A is denoted by AT and the complex conjugate of the complex variable z = + j! is denoted by z 3 . For symmetric matrices, A 0 (A 0) means that matrix A is definite positive (negative). The operator denotes the Kronecker product of two matrices, C = A B is a block matrix with block entry Cij = Aij B . The matrix Im is the identity matrix of size m 2 m. In this work, it is considered the definition of LMI region presented in [6]. Let R 2 2m22m be a symmetric matrix partitioned as
R=
R11 R12 T R22 R12
T 2 m2m R11 = R11 T 2 m2m ; R22 = R22
R22 0:
(1)
T z 3 + R22 zz 3 0 : R11 + R12 z + R12
(2)
:
The LMI -region is defined as
D
z2
:
The matrix A 2 n2n is said to be D -stable if and only if all its eigenvalues lie in the region D defined by (2). Consider the uncertain linear time-invariant system described by
[x(t)] = Ax(t); A 2 P
(3)
where [x(t)] dx(t)=dt and t 2 for continuous-time systems, or [x(t)] x(t + 1) and t 2 for discrete-time systems, x 2 n is the vector of state variables, and the matrix A is not precisely known, but
it belongs to a polytope type uncertain domain P defined as the set of all matrices obtained with the convex combination of its vertices
P
A 2 n2n : A =
N i=1
i Ai ;
N i=1
i = 1; i 0; 8i
(4)
where N is the number of polytope vertices and = [ 1 . . . N ]T is the polytope coordinate vector.
0018-9286/$20.00 © 2006 IEEE
1710
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
The system described by (3) is robustly D -stable if and only if A is
D-stable for all A 2 P .
Theorem 1 in [6] states that a system described by (3) is quadratically D-stable if and only if there exists a single symmetric positive definite matrix P , such that for all i = 1; . . . ; N
T (A T P ) R11 P + R12 (P Ai ) + R12 i T +R22 (Ai P Ai ) 0:
(5)
II. THE PROPOSED ROBUST D -STABILITY ANALYSIS APPROACH
Let LMI (P ) be a function that implements the LMI-based feasibility problem that characterizes the robust D -stability, where the polytope P is represented by its set of vertices fA1 ; . . . ; AN g. Given a feasibility problem of the form L(x) 0, where x is the vector of decision variables, LMI (P ) solves the auxiliary convex program: Minimize t subject to L(x) tI . This function returns the scalar tmin achieved in the minimization problem. The problem is feasible if and only if tmin 0, which means that the system is robustly D -stable for A 2 P. The main contribution of this work is the implementation of the analysis procedure presented in the sequel, where L is the set of simplices that does not attain the sufficient robust D -stability condition and (p) is either empty or the coordinate of the first system found in the polytope that is not D -stable. Robust D -Stability Analysis Procedure
;, p ;. Step 2. Compute t = LMI (P ), if t 0, then go to step 9. Step 3. Verify if all vertices of P are D -stable, if it is found a vertex whose eigenvalues do not lie in the D region, then ( )
min
(p)
min
(vertex coordinate)
and go to step 9.
Step 4. If P is a simplex, then set L = P , else use the Delaunay triangulation to decompose P into a set of simplices S = fS1 ; . . . ; Sr g and go to step 8. Step 5. Find the simplex Sm
2 S , if t
min
=
LMI (Si ) > 0, then L
L [ Si .
;, then go to step 5, else stop. At the end of the proposed analysis procedure, if L and p are both empty, then the uncertain system is robustly D -stable, else p 6 Step 9. If L =
; and p
( )
=
( )
( )
The main result in [6] is the Theorem 4 reproduced here: if there exists F 2 mn2mn , G 2 mn2n , and N symmetric positive–definite matrices Pi such that, for all i 2 f1; . . . ; N g; see (6) as shown at the bottom of the page; then, (3) is robustly D -stable. The proof can be found in [6]. Theorem 1 in [4] introduces more complex LMI conditions based on [6, Th. 4], but with the variables F and G replaced by parameterdependent variables Fi and Gi 2 mn2mn , i = 1; . . . ; N .
Step 1. Initialize L
Step 8. For all Si
2 L with the maximum value of t Step 6. Create the new vertices in Sm required by the edgewise
min .
subdivision technique.
Verify if all new vertices are D -stable, if it is found a vertex whose (vertex coordinate) and eigenvalues do not lie in D , then (p) go to step 9. Step 7. Split Sm in the set of simplices S = fS1 ; . . . ; Sr g using the edgewise subdivision technique and exclude Sm from L.
contains the polytope coordinate of the first system instance found in P that is not D-stable. If the LMI-based formulation is feasible for the whole polytope, the proposed analysis approach ends after step 2 without additional computational time. The efficiency and applicability of this analysis procedure depend on the appropriate choice of the sufficient conditions and polytope partition technique. The sufficient condition can be characterized by any LMI-based formulation but the best choice is not obvious since it is necessary to consider the tradeoff between complexity and conservatism. Remark: The strategy of polytope partition to reduce the conservatism of an LMI sufficient condition can be also applied to the computation of H2 or H1 "-guaranteed costs with a prescribed relative accuracy ". By means of the polytope partition technique presented in the next section, it is possible to implement an efficient branch-and-bound algorithm where the upper bound function is the maximum value of the guaranteed costs computed for each subpolytope P 2 L and the lower bound function is the maximum value of the norms computed for each vertex of each subpolytope. In this case, L is the set of polytopes that have guaranteed costs greater than the lower bound function. At each iteration, the subpolytope with the maximum guaranteed cost is subdivided until the difference between the upper and lower bound functions reach the relative prescribed accuracy ". III. THE PROPOSED POLYTOPE PARTITION TECHNIQUE The major problem in the implementation of the analysis procedure is the polytope subdivision approach to deal with polytope with any shape in a space of any dimension. In the case of polytopic models, the uncertain domain is a polytope with N vertices in the space of dimension (N 0 1), i.e., a simplex, parameterized by = [ 1 . . . N01 ]T . In the case of affine parameter-dependent models, the uncertain domain is a polytope in the uncertain parameter space of dimension d, parameterized by p = [p1 . . . pd ]T , generally a hyperrectangle with 2d vertices. The partition of P can be accomplished by several different ways. In this work, the partition is implemented considering simplicial meshes. This choice is due to: 1) in the case of polytopic models, the uncertain domain is already a simplex; 2) any polytope can be exactly decomposed in a set of simplices by a triangulation technique; and 3) the simplex is the polytope with the minimum number of vertices in any dimensional space which reduces the computational effort required by the LMI-based sufficient conditions. In the proposed implementation, in the case of parameter-dependent models, it is considered the well known Delaunay triangulation (MATLAB function delaunayn(1)) to decompose a two-dimensional polytope in triangles (tetrahedra in three-dimensional or simplices in d-dimensional spaces). It is necessary a simplex subdivision technique to implement the simplicial mesh refinement. It is essential to the convergence of branchand-bound like algorithms to avoid degenerated simplices, i.e., simplices with interior angles near to zero. Because of this, it is interesting that the simplex subdivision techinque generates a finite number of congruence classes. Two simplices ; 0 2 d are called congruent to each
R11 Pi + F (Im Ai ) + (Im ATi )F T R12 Pi + (Im ATi )G 0 F T P i + G T (I m A i ) 0 F T R22 Pi 0 G 0 GT R12
0
(6)
GONÇALVES et al.: LEE: LIANG AND XU: HAN: KWON AND PARK: AUTHORS’ REPLY
Fig. 1. Two steps of three possible subdivisions of a triangle (second step in dashed lines). (a) Division by two, or bisection. (b) Division by three. (c) Division by four, or edgewise subdivision.
order if there exists a translation vector v 2 d , a scaling factor c > 0, and an orthogonal matrix Q 2 d2d such that 0 = v + cQ [10]. Fig. 1 shows two subdivision steps by three possible triangle subdivision methods. The numbers inside the triangles in Fig. 1 correspond to the congruence classes, where the initial triangle belongs to class 1. The bisection subdivision method is simple to implement in any dimension but there is no guarantee that the number of congruence classes is limited for dimensions greater than 2. The division by three is obviously inappropriate even for two-dimensional space. The edgewise subdivision is the only one that generates triangles of the same class 1, regardless the number of refinements. It is proved in [10] that the number of congruence classes generated by the 2d simplex edgewise subdivision is limited to d!=2, which corresponds to the optimal value. This means that the maximum distance between vertices tends to zero faster than other techniques. Besides, all subsimplices have the same volume, which guarantees that the simplex volumes tend to zero with rate 1=2d . Due these two features, in this work, the simplicial mesh refinement is done by the simplex edgewise subdivision presented in [11] and implemented by the algorithm proposed in the next section. IV. EDGEWISE SUBDIVISION OF A SIMPLEX The main novelty of the proposed D -stability analysis procedure implementation is the formalization of a simple algorithm to implement the edgewise subdivision of a simplex. The proposed algorithm is based on the abacus model of the simplex presented in [11] which can be applied to subdivide a simplex in the d-dimensional space into kd simplices. It is proved in [11] that the number of congruence classes generated by the kd simplex edgewise subdivision is also limited to d!=2 independently of k . Although examples of simplex subdivision are presented in [11], it is not formalized a “ready-to-implement” algorithm as the one proposed here. The same notation used in [11] is applied here. Consider a d-simplex defined as a set of d + 1 points, P0 ; P1 ; . . . ; Pd , that are affinely independent in d . Consider the notation
P ...
=
(
P
+
P
+
k
111
+
P
)
1711
Fig. 2. New vertices in the partition of a tetrahedron by means of a 2 edgewise subdivision and the sub-tetrahedron achieved with .
M
TABLE I NUMBER OF SIMPLICES GENERATED BY THE SUBDIVISION TECHNIQUE
rows of M start with “0” instead of “1” as presented in [11]. One contribution of this work is just to formulate an algorithm to compute the kd color schemes to achieve the subdivision. Let n i;j be the entry of the ith row and j th column of the nth color scheme M n . Each color scheme M n , n = 0; 1; . . . ; kd 0 1, will be created row by row starting with 0n;0 = 0. To find out if the next entry will be kept or increased by one, it is necessary to represent the index n of the new simplex in the numerical system with base k and d digits
n = xd01 2 kd01 + xd02 2 kd02 + 1 1 1 + x0 2 k0 :
The values of the digits xd0i , i = 1; 2; . . . ; d, will determine which row of column i 0 1 will be increased by one to generate the column i. The following row starts with the last color in the preceding row, i.e., i;n0 = in01;d . This idea is implemented by the following algorithm. Algorithm Color scheme for n
; ;
= 0 1 ...
xd01 . . . x0 color
0
; kd 0 1 do conversion of n to base k ;
;
for i = 0; 1; . . . ; k 0 1 do i;n0 color; for j
(7)
where i 2 [0; d], is a vertex index. The edgewise subdivision of in kd will be derived from the vertices P0 ; P1 ; . . . ; Pd and new vertices P ... as defined in (7). The vertices which define each new simplex will be obtained from a matrix M 2 k2(d+1) , called the color scheme, whose entries, i;j , are integer numbers in the range [0, d], called the colors, that represent the indices of the vertices P0 ; P1 ; . . . ; Pd . The ith column of M will define the ith vertex P ... of the new simplex. In this work, the indices of the
(8)
;
= 1 ...
; d do
if xd0j = i then color n i;j color; end end end end
color + 1;
1712
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
With k = 2, the edgewise subdivision will be accomplished by introducing one new vertex over each edge of the d-simplex P , computed as Pij = (Pi + Pj )=2, i = 0; . . . ; d 0 1, j = i + 1; . . . ; d. Consider the subdivision of a tetrahedron in kd = 23 subtetrahedrons. The color scheme is formatted as
M
=
TABLE II PERFORMANCE OF THE PROPOSED ANALYSIS PROCEDURE TO IDENTIFY A ROBUSTLY -STABLE SYSTEM (EXAMPLE 1)
D
0;0 0;1 0;2 0;3 : 1;0 1;1 1;2 1;3
For example, to compute the sub-tetrahedron with n = 3, the change of colors will be specified by 3 in the base 2, i.e., x = 0112 , which means that the first color change occurs in row 0 and the next two in row 1
M3 =
0
1
1
1
1
1
2
3
:
The columns of this color scheme show that the subtetrahedron is defined by the vertex set fP01 ; P11 = P1 ; P12 ; P13 g as shown in Fig. 2, where the point P is computed by (7): Pij = (Pi + Pj )=2. In the proposed analysis approach it is considered only k = 2. The computational time required by the proposed analysis approach depends mainly on the number of evaluations of the sufficient robust stability condition for each sub-polytope generated in each step of the polytope partition. It also depends on the complexity of the LMI-based analysis formulation, i.e., the number of decision variables and constraints. Table I shows how the number of required evaluations, in each algorithm step, increases with the space dimension.
V. NUMERICAL ILLUSTRATIVE EXAMPLES The following examples will be presented to illustrate the effectiveness of the proposed approach to robust D -stability analysis. The proposed analysis procedure was implemented in both LMI Control Toolbox (LMIC) [12] and SeDuMi (SDMi) [13] for use with MATLAB. The computational complexity of LMIC is O(K 3 L) while SDMi has a complexity O(K 2 L2:5 + L3:5 ), where K is the total number of scalar decision variables and L is the total row size of the LMI system. In both cases the feasibility radius is set to 109 and the maximum number of iteration is set to 500. The implementation allows the partition of the uncertain domain based on the polytope coordinates (polytopic models) or uncertain parameter vector (affine parameter-dependent models). The following results have been computed in a Pentium IV 2.8 GHz with 1 Gb of system RAM and 4096 MB of virtual memory. The examples will consider the formulations based on: Quadratic stability characterization, with single Lyapunov function, as presented in [6, Th. 1]; [6, Th. 4]; and [4, Th. 1], both with parameter-dependent Lyapunov functions. The total row size L and the number of decision scalar variables K are as follows: For [6, Th. 1], L = n(N + 1) (half-plane and disk), L = n(2N + 1) (sector), K = n(n + 1)=2; for [6, Th. 4], L = 3nN , (half-plane and disk), L = 5nN (sector), K = 2n2 + Nn(n + 1)=2 (half-plane and disk), K = 2(2n)2 + Nn(n + 1)=2 (sector); for [4, Th. 1], L = n(N 3 + 3N 2 + 5N )=3, (half-plane and disk), L = n(4N 3 + 12N 2 + 14N )=6 (sector), K = N [2n2 + n(n + 1)=2] (half-plane and disk), K = N [2(2n)2 + n(n + 1)=2] (sector). Instead of analyzing all region constraints simultaneously, it is considered one constraint at each time. This strategy results in smaller computational times.
TABLE III PERFORMANCE OF THE PROPOSED ANALYSIS PROCEDURE TO IDENTIFY A NOT ROBUSTLY -STABLE SYSTEM (EXAMPLE 1)
D
1) Example 1: Consider the stabilization and disturbance attenuation problem of the longitudinal short period mode of the F4E fighter aircraft analyzed in [14]
dx dt
=
y=
a11 a12 a13 b1 a21 a22 a23 x + 0 u 0
x1 x2
0
030
30
where the uncertain parameters, determined in four operation points, are presented in [14, Table I]. Based on the approach presented in [15], a static output feedback controller is designed to robustly place the closed-loop eigenvalues into the intersection of the half-plane Real(z ) < 01:25, the disk with radius r = 25 centered at the origin, and the conic sector with inner angle 2 = =2. It is considered the polytope defined by the four vertices related to the four operation points. The design procedure computes the controller K = [0:0638 0:8034]. In this example, the LMI-based analysis formulations based on [6, Th. 1] (quadratic stability) and on [6, Th. 4] are infeasible for all regions. Only the formulation based on [4, Th. 1] is capable to identify the robust D -stability, but the required total computational time is greater than the proposed analysis procedure implemented with [6, Th. 4], for both solvers, as presented in Table II. Considering the robust D -stability analysis of the three regions simultaneously, it is required 213, 25, and 171 s, for the proposed approach (LMIC solver) implemented with [6, Th. 1], [6, Th. 4], and [4, Th. 1], respectively. These computational times are higher than the ones presented in Table II. Consider now a controller that does not result in a robustly D -stable system: K = [0:0639 0:8044]. The proposed analysis procedure finds the polytope coordinates (p) = [0 0:6875 0 0:3125] where the closed-loop system is not D -stable (there are poles outside the disk region). Table III compares the performance of the proposed analysis
GONÇALVES et al.: LEE: LIANG AND XU: HAN: KWON AND PARK: AUTHORS’ REPLY
1713
Fig. 3. Partition of the uncertain domain (Example 1).
procedure, implemented with the three LMI sufficient conditions, to identify a not robustly D -stable system for the case of the disk region. Fig. 3 shows the uncertain domain partition when considering [6, Th. 4]. The arrow indicates the first polytope coordinate found where the system is not D -stable. The uncertain domain is a tetrahedron, in the three-dimensional space, parameterized by the first three polytope coordinates since 4 = 1 0 3i=1 i . 2) Example 2: Consider the state–space model of the linearized vertical-plane dynamics of an aircraft (example (AC1) in [16]), with five uncertain parameters included here
0
dx dt
0
=
0
0
0:0485
0
0
p1
0 0
0
00:12
0
0
1
0
0
0
0
+
y=
1:132
0
0
00:0538 00:1712
0
p4 p5
0 0
01:665 00:0732
1
0
0
0
0
0
1
0
0
0
0
0
1
0
0
[0:2618; 0:3200]
p
0:0705
1
0
09:8556 01:013 p2 0 p3
x
disk centered in the origin with radius r = 3:6, and the conic sector with internal angle 2 = 2=4. The controller matrices are given by
u
= [ p1 ; . . . ; p 5 ]
D
01
0
x:
The uncertain parameter vector perrectangle
TABLE IV PERFORMANCE OF THE PROPOSED ROBUSTLY -STABILITY ANALYSIS APPROACH (EXAMPLE 2)
varies in the hy-
2 [0:9479; 1:1585] 2 [0:6173; 0:7545] 2[3:9771; 4:8609] 2 [1:4175; 1:7325]
with 32 vertices in the five-dimensional space. Based on the design procedure presented in [15], it was designed a first-order dynamic output controller to robustly place the closed-loop poles in the intersection of the half-plane region Real(z ) < 01:5, the
Ac = 0 5:1437 Bc = [ 5:2419 0:0817 00:2764 ] 05:7436 7:4180 00:1314 0:4392 Cc = 4:7128 Dc = 04:2613 1:3786 6:0419 : 07:2602 17:5314 9:7853 25:8029 In this example, with the formulation based on [4, Th. 1], the LMIC solver generates a memory error for the the initial polytope analysis which does not occur with the SDMi solver, but the latter requires prohibitive computational time in relation to the other LMI formulations. Table IV compares the performance of the proposed analysis approach implemented with the other two LMI analysis formulations. Without the polytope partition, both LMI formulations are not capable to inform if the uncertain system is robustly D -stable. Based on the proposed analysis approach implemented with [6, Th. 4], it is verified that the system is robustly D -stable. Considering [6, Th. 1] (quadratic stability), the proposed procedure reaches the maximum number of iter-
1714
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
ations without identifying the robust stability for the two first regions due to the high conservatism of this sufficient condition. As presented in Table I, in this case, 103 evaluations of the sufficient stability conditions are necessary in the first step of the polytope partition (Delaunay triangulation) and 32 evaluations in each further refinement (simplex subdivision). To verify the three regions simultaneously, the proposed procedure spends 23 subdivisions and 52989 s considering [6, Th. 4] and SDMi solver. One can observe in Table IV the reduction of the complexity (K, L) when working with simplex instead hyperrectangle. Considering the results in Table IV and the memory problem with the analysis formulation based on [4, Th. 1], even for higher dimensional spaces, the proposed analysis approach, with the appropriate choice of the sufficient condition, can be more computationally efficient than LMI-based formulations alone.
VI. CONCLUSION It has been proposed an approach to establish if a linear time-invariant system, with polytope-bounded uncertainty, is robustly D-stable or not. Due to the conservativeness of the LMI-based analysis formulations, when the LMI-solvers do not find a solution to the feasibility problem it is not possible to state that the system is not robustly D -stable. If the system is D -stable, the proposed approach iteratively subdivides the uncertain domain until the LMI sufficient condition is verified for all subpolytopes. In the contrary case, the proposed approach terminates when any system related to a subpolytope vertex is verified to be not robustly D -stable. In spite of the implementation complexity, the proposed approach can spend less computational time than LMI-based formulations with a large number of decision variables and constraints that are more influenced by the number of polytope vertices and system order. The main difficulty in the implementation of the proposed approach is the polytope partition. To overcome this problem, it has been proposed a new approach that combines the Delaunay triangulation, to decompose a polytope with any shape in a set of simplices, and a new simplex edgewise subdivision algorithm.
REFERENCES [1] M. C. de Oliveira, J. Bernussou, and J. C. Geromel, “A new discretetime robust stability condition,” Syst. Control Lett., vol. 37, no. 4, pp. 261–265, 1999. [2] D. C. W. Ramos and L. D. Peres, “A less conservative LMI condition for the robust stability of discrete-time uncertain systems,” Syst. Control Lett., vol. 43, pp. 371–378, Apr. 2001. [3] D. C. W. Ramos and P. L. D. Peres, “An LMI condition for the robust stability of uncertain continuous-time linear systems,” IEEE Trans. Autom. Control, vol. 47, no. 4, pp. 675–678, Apr. 2002. [4] V. J. S. Leite and P. L. D. Peres, “An improved LMI condition for robust -stability of uncertain polytopic systems,” IEEE Trans. Autom. Control, vol. 48, no. 3, pp. 500–504, Mar. 2003. [5] M. Chilali and P. Gahinet, “ design with pole placement constraints: An LMI approach,” IEEE Trans. Autom. Control, vol. 41, no. 3, pp. 358–367, Mar. 1996. [6] D. Peaucelle, D. Arzelier, O. Bachelier, and J. Bernussou, “A new robust -stability condition for real convex polytopic uncertainty,” Syst. Control Lett., vol. 40, pp. 21–30, 2000. [7] C. DeMarco, V. Balakrishnan, and S. Boyd, “A branch and bound methodology for matrix polytope stability problems arising in power systems,” in Proc. IEEE Conf. Decision and Control, Honolulu, HI, Dec. 1990, pp. 3022–3027. [8] V. Balakrishnan, S. Boyd, and S. Balemi, “Branch and bound algorithm for computing the minimum stability degree of parameter-dependent linear systems,” Int. J. Robust Nonlinear Control, vol. 1, no. 4, pp. 295–317, Oct.–Dec. 1991.
D
H
D
[9] D. W. Moore, Simplicial mesh generation with applications. Ithaca, NY, Dept. Comput. Sci., Cornell Univ., 1992, Ph.D. dissertation. [10] J. Bey, “Simplicial grid refinement: On Freudenthal’s algorithm and the optimal number of congruence classes,” Numerische Mathematik, vol. 85, pp. 1–29, 2000. [11] H. Edelsbrunner and D. R. Grayson, “Edgewise subdivision of a simplex,” Discrete Comput. Geom., vol. 24, pp. 707–719, 2000. [12] P. Gahinet, A. Nemirovski, A. J. Laub, and M. Chilali, LMI Control Toolbox: For Use with MATLAB. Natick, MA: The MathWorks, Inc., 1995. [13] J. Sturm, “Using SeDuMi 1.02, A MATLAB toolbox for optimization over symmetric cones,” Optim. Meth. Software vol. 11–12, pp. 625–653, 1999 [Online]. Available: http://sedumi.mcmaster.ca, ed. 1.1 and static [14] U. Shaked, “An LPD approach to robust output-feedback design,” IEEE Trans. Autom. Control, vol. 48, no. 5, pp. 866–872, May 2003. [15] E. N. Gonçalves, R. M. Palhares, and R. H. C. Takahashi, “Improved = control problem for linear optimisation approach to robust systems,” Proc. Inst. Elect. Eng. Control Theory Appl., vol. 152, no. 2, pp. 171–176, 2005. [16] F. Leibfritz and W. Lipinski, Description of the benchmark examples in COMP l ib 1.0 Univ. Trier, Trier, Germany, 2003, Tech. Rep.
H
H
H H
On Uniform Stabilization of Discrete-Time Linear Parameter-Varying Control Systems Ji-Woong Lee
Abstract—For discrete-time polytopic linear parameter-varying systems, the uniform exponential stability, which is equivalent to the asymptotic stability and includes the quadratic stability as a special case, is characterized by the union of an increasing family of linear matrix inequality conditions. On the other hand, in certain cases, nonconservative synthesis of uniformly stabilizing controllers is achieved via a system of finite number of linear matrix inequalities if the controller is allowed to have finite memory of past parameters. Two such cases are robust state feedback (in a relaxed sense) against polytopic parameter variations, and multiple output injections under polytopic fusion rules. Index Terms—Discrete linear inclusion, linear matrix inequality (LMI), linear parameter-varying (LPV) system, output injection, state feedback.
I. INTRODUCTION A common approach to analysis and control of linear systems whose parameters vary with uncertainty or in a complicated way (possibly depending on the state) is first to approximate the set of coefficients by a convex polytope and then to abstract the system as the family of all linear time-varying systems whose coefficients vary within this polytope. The resulting abstractions are called (polytopic) linear parameter-varying (LPV) systems. An LPV system is said to be stable if
Manuscript received September 9, 2005; revised February 3, 2006. Recommended by Associate Editor Y. Ohta. The author is with the Department of Electrical and Computer Engineering, the University of Florida, Gainesville, FL 32611 USA (e-mail: jiwoong@ufl. edu). Digital Object Identifier 10.1109/TAC.2006.880804
0018-9286/$20.00 © 2006 IEEE
New Approach to Robust -Stability Analysis of Linear Time-Invariant Systems With Polytope-Bounded Uncertainty Eduardo N. Gonçalves, Reinaldo M. Palhares, Ricardo H. C. Takahashi, and Renato C. Mesquita
Abstract—This note presents a new approach to robust -stability analysis of linear time-invariant systems with polytope-bounded uncertainty. The proposed approach combines sufficient conditions for robust -stability in terms of feasibility problems with linear matrix inequalities (LMI) constraints and a new polytope partition technique. If the initial polytope does not attain the robust -stability sufficient condition, the polytope is successively subdivided until all subpolytopes attain the sufficient condition, in the case of robustly -stable uncertain system, or it is found a subpolytope vertex that does not attain the regional pole-placement constraints, in the case of an uncertain system that is not robustly -stable. It is presented a new general format polytope partition technique that allows the implementation of the proposed approach. The efficiency of the proposed approach is verified by means of illustrative examples and three different LMI-based analysis formulations. Index Terms—Polytope-bounded uncertainty, robust stability, simplex subdivision.
1709
D-stable. This branch-and-bound like strategy has been already considered for robust stability analysis, e.g., [7], [8]. The partition technique considered in these previous works is based on the bisection of hyperrectangles. Because of this, these works are restrict to the case of systems with uncertain parameters that range between extremal values. The major contribution of this work is not only to extend this strategy to robust D -stability analysis, but it is how this strategy is implemented based on up-to-date LMI-based analysis formulations and a new polytope partition technique based on simplicial meshes. The advantage of considering simplicial meshes is the capability for dealing with polytopes with any shape, not restricted to the hyperrectangle case as in [7], [8]. This feature allows the proposed analysis procedure to be applied to both affine parameter-dependent and polytopic models, what is not possible with the aforementioned works. A complete discussion about advantages of simplicial meshes over box meshes and other mesh generation approaches can be found in [9]. Without the conservatism of LMI-based sufficient conditions, the proposed strategy allows the stability analysis of closed-loop robust systems designed through approaches not based only on LMI framework, as illustrated in the examples.
-stability, robust
NOTATION AND DEFINITIONS
I. INTRODUCTION The Lyapunov theory has been frequently employed to robust stability analysis due to the fact that the resulting linear matrix inequality (LMI)-based feasibility problem is easy to solve by one of the available LMI-solvers. The quadratic stability condition based on a single Lyapunov function is the simplest formulation but it can be too conservative in the case of time-invariant systems. To reduce the conservativeness, one can use parameter-dependent Lyapunov functions and slack variables [1]–[3]. Unfortunately, these LMI-based analysis formulations are only sufficient conditions. As verified in [4], the reliability of the LMI-based analysis formulations decreases with the increasing of the number of polytope vertices and the order of the system. It is possible to formulate a more reliable formulation increasing the number of decision variables but at the cost of a higher computational cost. The Lyapunov theory can be extended to deal with the problem of robust D-stability to guarantee that all closed-loop poles are robustly placed at convex regions of the complex plane, denoted LMI regions [5], [6]. In this note, it is presented an analysis approach that finds out if a linear time-invariant system is robustly D -stable or not. To accomplish this task, the proposed analysis approach is based on a procedure that combines a D -stability sufficient condition and a polytope partition technique. The proposed procedure subdivides the polytope iteratively until it is verified that all generated subpolytopes are robustly D-stable or it is found a system instance of the polytope that is not Manuscript received June 3, 2005; revised December 7, 2005 and March 29, 2006. Recommended by Associate Editor Y. Ohta. This work was supported in part by the Brazilian agencies CNPq and FAPEMIG. E. N. Gonçalves is with the Department of Electrical Engineering, Federal Center of Technological Education of Minas Gerais, 31510-470 Belo Horizonte, MG, Brazil (e-mail: [email protected]). R. M. Palhares and R. C. Mesquita are with the Department of Electronics Engineering, Federal University of Minas Gerais, 31270-010 Belo Horizonte, MG, Brazil (e-mail: [email protected]; [email protected]). R. H. C. Takahashi is with the Department of Mathematics, Federal University of Minas Gerais, 31270-010 Belo Horizonte, MG, Brazil (e-mail: [email protected]). Digital Object Identifier 10.1109/TAC.2006.883061
The transpose of a matrix A is denoted by AT and the complex conjugate of the complex variable z = + j! is denoted by z 3 . For symmetric matrices, A 0 (A 0) means that matrix A is definite positive (negative). The operator denotes the Kronecker product of two matrices, C = A B is a block matrix with block entry Cij = Aij B . The matrix Im is the identity matrix of size m 2 m. In this work, it is considered the definition of LMI region presented in [6]. Let R 2 2m22m be a symmetric matrix partitioned as
R=
R11 R12 T R22 R12
T 2 m2m R11 = R11 T 2 m2m ; R22 = R22
R22 0:
(1)
T z 3 + R22 zz 3 0 : R11 + R12 z + R12
(2)
:
The LMI -region is defined as
D
z2
:
The matrix A 2 n2n is said to be D -stable if and only if all its eigenvalues lie in the region D defined by (2). Consider the uncertain linear time-invariant system described by
[x(t)] = Ax(t); A 2 P
(3)
where [x(t)] dx(t)=dt and t 2 for continuous-time systems, or [x(t)] x(t + 1) and t 2 for discrete-time systems, x 2 n is the vector of state variables, and the matrix A is not precisely known, but
it belongs to a polytope type uncertain domain P defined as the set of all matrices obtained with the convex combination of its vertices
P
A 2 n2n : A =
N i=1
i Ai ;
N i=1
i = 1; i 0; 8i
(4)
where N is the number of polytope vertices and = [ 1 . . . N ]T is the polytope coordinate vector.
0018-9286/$20.00 © 2006 IEEE
1710
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
The system described by (3) is robustly D -stable if and only if A is
D-stable for all A 2 P .
Theorem 1 in [6] states that a system described by (3) is quadratically D-stable if and only if there exists a single symmetric positive definite matrix P , such that for all i = 1; . . . ; N
T (A T P ) R11 P + R12 (P Ai ) + R12 i T +R22 (Ai P Ai ) 0:
(5)
II. THE PROPOSED ROBUST D -STABILITY ANALYSIS APPROACH
Let LMI (P ) be a function that implements the LMI-based feasibility problem that characterizes the robust D -stability, where the polytope P is represented by its set of vertices fA1 ; . . . ; AN g. Given a feasibility problem of the form L(x) 0, where x is the vector of decision variables, LMI (P ) solves the auxiliary convex program: Minimize t subject to L(x) tI . This function returns the scalar tmin achieved in the minimization problem. The problem is feasible if and only if tmin 0, which means that the system is robustly D -stable for A 2 P. The main contribution of this work is the implementation of the analysis procedure presented in the sequel, where L is the set of simplices that does not attain the sufficient robust D -stability condition and (p) is either empty or the coordinate of the first system found in the polytope that is not D -stable. Robust D -Stability Analysis Procedure
;, p ;. Step 2. Compute t = LMI (P ), if t 0, then go to step 9. Step 3. Verify if all vertices of P are D -stable, if it is found a vertex whose eigenvalues do not lie in the D region, then ( )
min
(p)
min
(vertex coordinate)
and go to step 9.
Step 4. If P is a simplex, then set L = P , else use the Delaunay triangulation to decompose P into a set of simplices S = fS1 ; . . . ; Sr g and go to step 8. Step 5. Find the simplex Sm
2 S , if t
min
=
LMI (Si ) > 0, then L
L [ Si .
;, then go to step 5, else stop. At the end of the proposed analysis procedure, if L and p are both empty, then the uncertain system is robustly D -stable, else p 6 Step 9. If L =
; and p
( )
=
( )
( )
The main result in [6] is the Theorem 4 reproduced here: if there exists F 2 mn2mn , G 2 mn2n , and N symmetric positive–definite matrices Pi such that, for all i 2 f1; . . . ; N g; see (6) as shown at the bottom of the page; then, (3) is robustly D -stable. The proof can be found in [6]. Theorem 1 in [4] introduces more complex LMI conditions based on [6, Th. 4], but with the variables F and G replaced by parameterdependent variables Fi and Gi 2 mn2mn , i = 1; . . . ; N .
Step 1. Initialize L
Step 8. For all Si
2 L with the maximum value of t Step 6. Create the new vertices in Sm required by the edgewise
min .
subdivision technique.
Verify if all new vertices are D -stable, if it is found a vertex whose (vertex coordinate) and eigenvalues do not lie in D , then (p) go to step 9. Step 7. Split Sm in the set of simplices S = fS1 ; . . . ; Sr g using the edgewise subdivision technique and exclude Sm from L.
contains the polytope coordinate of the first system instance found in P that is not D-stable. If the LMI-based formulation is feasible for the whole polytope, the proposed analysis approach ends after step 2 without additional computational time. The efficiency and applicability of this analysis procedure depend on the appropriate choice of the sufficient conditions and polytope partition technique. The sufficient condition can be characterized by any LMI-based formulation but the best choice is not obvious since it is necessary to consider the tradeoff between complexity and conservatism. Remark: The strategy of polytope partition to reduce the conservatism of an LMI sufficient condition can be also applied to the computation of H2 or H1 "-guaranteed costs with a prescribed relative accuracy ". By means of the polytope partition technique presented in the next section, it is possible to implement an efficient branch-and-bound algorithm where the upper bound function is the maximum value of the guaranteed costs computed for each subpolytope P 2 L and the lower bound function is the maximum value of the norms computed for each vertex of each subpolytope. In this case, L is the set of polytopes that have guaranteed costs greater than the lower bound function. At each iteration, the subpolytope with the maximum guaranteed cost is subdivided until the difference between the upper and lower bound functions reach the relative prescribed accuracy ". III. THE PROPOSED POLYTOPE PARTITION TECHNIQUE The major problem in the implementation of the analysis procedure is the polytope subdivision approach to deal with polytope with any shape in a space of any dimension. In the case of polytopic models, the uncertain domain is a polytope with N vertices in the space of dimension (N 0 1), i.e., a simplex, parameterized by = [ 1 . . . N01 ]T . In the case of affine parameter-dependent models, the uncertain domain is a polytope in the uncertain parameter space of dimension d, parameterized by p = [p1 . . . pd ]T , generally a hyperrectangle with 2d vertices. The partition of P can be accomplished by several different ways. In this work, the partition is implemented considering simplicial meshes. This choice is due to: 1) in the case of polytopic models, the uncertain domain is already a simplex; 2) any polytope can be exactly decomposed in a set of simplices by a triangulation technique; and 3) the simplex is the polytope with the minimum number of vertices in any dimensional space which reduces the computational effort required by the LMI-based sufficient conditions. In the proposed implementation, in the case of parameter-dependent models, it is considered the well known Delaunay triangulation (MATLAB function delaunayn(1)) to decompose a two-dimensional polytope in triangles (tetrahedra in three-dimensional or simplices in d-dimensional spaces). It is necessary a simplex subdivision technique to implement the simplicial mesh refinement. It is essential to the convergence of branchand-bound like algorithms to avoid degenerated simplices, i.e., simplices with interior angles near to zero. Because of this, it is interesting that the simplex subdivision techinque generates a finite number of congruence classes. Two simplices ; 0 2 d are called congruent to each
R11 Pi + F (Im Ai ) + (Im ATi )F T R12 Pi + (Im ATi )G 0 F T P i + G T (I m A i ) 0 F T R22 Pi 0 G 0 GT R12
0
(6)
GONÇALVES et al.: LEE: LIANG AND XU: HAN: KWON AND PARK: AUTHORS’ REPLY
Fig. 1. Two steps of three possible subdivisions of a triangle (second step in dashed lines). (a) Division by two, or bisection. (b) Division by three. (c) Division by four, or edgewise subdivision.
order if there exists a translation vector v 2 d , a scaling factor c > 0, and an orthogonal matrix Q 2 d2d such that 0 = v + cQ [10]. Fig. 1 shows two subdivision steps by three possible triangle subdivision methods. The numbers inside the triangles in Fig. 1 correspond to the congruence classes, where the initial triangle belongs to class 1. The bisection subdivision method is simple to implement in any dimension but there is no guarantee that the number of congruence classes is limited for dimensions greater than 2. The division by three is obviously inappropriate even for two-dimensional space. The edgewise subdivision is the only one that generates triangles of the same class 1, regardless the number of refinements. It is proved in [10] that the number of congruence classes generated by the 2d simplex edgewise subdivision is limited to d!=2, which corresponds to the optimal value. This means that the maximum distance between vertices tends to zero faster than other techniques. Besides, all subsimplices have the same volume, which guarantees that the simplex volumes tend to zero with rate 1=2d . Due these two features, in this work, the simplicial mesh refinement is done by the simplex edgewise subdivision presented in [11] and implemented by the algorithm proposed in the next section. IV. EDGEWISE SUBDIVISION OF A SIMPLEX The main novelty of the proposed D -stability analysis procedure implementation is the formalization of a simple algorithm to implement the edgewise subdivision of a simplex. The proposed algorithm is based on the abacus model of the simplex presented in [11] which can be applied to subdivide a simplex in the d-dimensional space into kd simplices. It is proved in [11] that the number of congruence classes generated by the kd simplex edgewise subdivision is also limited to d!=2 independently of k . Although examples of simplex subdivision are presented in [11], it is not formalized a “ready-to-implement” algorithm as the one proposed here. The same notation used in [11] is applied here. Consider a d-simplex defined as a set of d + 1 points, P0 ; P1 ; . . . ; Pd , that are affinely independent in d . Consider the notation
P ...
=
(
P
+
P
+
k
111
+
P
)
1711
Fig. 2. New vertices in the partition of a tetrahedron by means of a 2 edgewise subdivision and the sub-tetrahedron achieved with .
M
TABLE I NUMBER OF SIMPLICES GENERATED BY THE SUBDIVISION TECHNIQUE
rows of M start with “0” instead of “1” as presented in [11]. One contribution of this work is just to formulate an algorithm to compute the kd color schemes to achieve the subdivision. Let n i;j be the entry of the ith row and j th column of the nth color scheme M n . Each color scheme M n , n = 0; 1; . . . ; kd 0 1, will be created row by row starting with 0n;0 = 0. To find out if the next entry will be kept or increased by one, it is necessary to represent the index n of the new simplex in the numerical system with base k and d digits
n = xd01 2 kd01 + xd02 2 kd02 + 1 1 1 + x0 2 k0 :
The values of the digits xd0i , i = 1; 2; . . . ; d, will determine which row of column i 0 1 will be increased by one to generate the column i. The following row starts with the last color in the preceding row, i.e., i;n0 = in01;d . This idea is implemented by the following algorithm. Algorithm Color scheme for n
; ;
= 0 1 ...
xd01 . . . x0 color
0
; kd 0 1 do conversion of n to base k ;
;
for i = 0; 1; . . . ; k 0 1 do i;n0 color; for j
(7)
where i 2 [0; d], is a vertex index. The edgewise subdivision of in kd will be derived from the vertices P0 ; P1 ; . . . ; Pd and new vertices P ... as defined in (7). The vertices which define each new simplex will be obtained from a matrix M 2 k2(d+1) , called the color scheme, whose entries, i;j , are integer numbers in the range [0, d], called the colors, that represent the indices of the vertices P0 ; P1 ; . . . ; Pd . The ith column of M will define the ith vertex P ... of the new simplex. In this work, the indices of the
(8)
;
= 1 ...
; d do
if xd0j = i then color n i;j color; end end end end
color + 1;
1712
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
With k = 2, the edgewise subdivision will be accomplished by introducing one new vertex over each edge of the d-simplex P , computed as Pij = (Pi + Pj )=2, i = 0; . . . ; d 0 1, j = i + 1; . . . ; d. Consider the subdivision of a tetrahedron in kd = 23 subtetrahedrons. The color scheme is formatted as
M
=
TABLE II PERFORMANCE OF THE PROPOSED ANALYSIS PROCEDURE TO IDENTIFY A ROBUSTLY -STABLE SYSTEM (EXAMPLE 1)
D
0;0 0;1 0;2 0;3 : 1;0 1;1 1;2 1;3
For example, to compute the sub-tetrahedron with n = 3, the change of colors will be specified by 3 in the base 2, i.e., x = 0112 , which means that the first color change occurs in row 0 and the next two in row 1
M3 =
0
1
1
1
1
1
2
3
:
The columns of this color scheme show that the subtetrahedron is defined by the vertex set fP01 ; P11 = P1 ; P12 ; P13 g as shown in Fig. 2, where the point P is computed by (7): Pij = (Pi + Pj )=2. In the proposed analysis approach it is considered only k = 2. The computational time required by the proposed analysis approach depends mainly on the number of evaluations of the sufficient robust stability condition for each sub-polytope generated in each step of the polytope partition. It also depends on the complexity of the LMI-based analysis formulation, i.e., the number of decision variables and constraints. Table I shows how the number of required evaluations, in each algorithm step, increases with the space dimension.
V. NUMERICAL ILLUSTRATIVE EXAMPLES The following examples will be presented to illustrate the effectiveness of the proposed approach to robust D -stability analysis. The proposed analysis procedure was implemented in both LMI Control Toolbox (LMIC) [12] and SeDuMi (SDMi) [13] for use with MATLAB. The computational complexity of LMIC is O(K 3 L) while SDMi has a complexity O(K 2 L2:5 + L3:5 ), where K is the total number of scalar decision variables and L is the total row size of the LMI system. In both cases the feasibility radius is set to 109 and the maximum number of iteration is set to 500. The implementation allows the partition of the uncertain domain based on the polytope coordinates (polytopic models) or uncertain parameter vector (affine parameter-dependent models). The following results have been computed in a Pentium IV 2.8 GHz with 1 Gb of system RAM and 4096 MB of virtual memory. The examples will consider the formulations based on: Quadratic stability characterization, with single Lyapunov function, as presented in [6, Th. 1]; [6, Th. 4]; and [4, Th. 1], both with parameter-dependent Lyapunov functions. The total row size L and the number of decision scalar variables K are as follows: For [6, Th. 1], L = n(N + 1) (half-plane and disk), L = n(2N + 1) (sector), K = n(n + 1)=2; for [6, Th. 4], L = 3nN , (half-plane and disk), L = 5nN (sector), K = 2n2 + Nn(n + 1)=2 (half-plane and disk), K = 2(2n)2 + Nn(n + 1)=2 (sector); for [4, Th. 1], L = n(N 3 + 3N 2 + 5N )=3, (half-plane and disk), L = n(4N 3 + 12N 2 + 14N )=6 (sector), K = N [2n2 + n(n + 1)=2] (half-plane and disk), K = N [2(2n)2 + n(n + 1)=2] (sector). Instead of analyzing all region constraints simultaneously, it is considered one constraint at each time. This strategy results in smaller computational times.
TABLE III PERFORMANCE OF THE PROPOSED ANALYSIS PROCEDURE TO IDENTIFY A NOT ROBUSTLY -STABLE SYSTEM (EXAMPLE 1)
D
1) Example 1: Consider the stabilization and disturbance attenuation problem of the longitudinal short period mode of the F4E fighter aircraft analyzed in [14]
dx dt
=
y=
a11 a12 a13 b1 a21 a22 a23 x + 0 u 0
x1 x2
0
030
30
where the uncertain parameters, determined in four operation points, are presented in [14, Table I]. Based on the approach presented in [15], a static output feedback controller is designed to robustly place the closed-loop eigenvalues into the intersection of the half-plane Real(z ) < 01:25, the disk with radius r = 25 centered at the origin, and the conic sector with inner angle 2 = =2. It is considered the polytope defined by the four vertices related to the four operation points. The design procedure computes the controller K = [0:0638 0:8034]. In this example, the LMI-based analysis formulations based on [6, Th. 1] (quadratic stability) and on [6, Th. 4] are infeasible for all regions. Only the formulation based on [4, Th. 1] is capable to identify the robust D -stability, but the required total computational time is greater than the proposed analysis procedure implemented with [6, Th. 4], for both solvers, as presented in Table II. Considering the robust D -stability analysis of the three regions simultaneously, it is required 213, 25, and 171 s, for the proposed approach (LMIC solver) implemented with [6, Th. 1], [6, Th. 4], and [4, Th. 1], respectively. These computational times are higher than the ones presented in Table II. Consider now a controller that does not result in a robustly D -stable system: K = [0:0639 0:8044]. The proposed analysis procedure finds the polytope coordinates (p) = [0 0:6875 0 0:3125] where the closed-loop system is not D -stable (there are poles outside the disk region). Table III compares the performance of the proposed analysis
GONÇALVES et al.: LEE: LIANG AND XU: HAN: KWON AND PARK: AUTHORS’ REPLY
1713
Fig. 3. Partition of the uncertain domain (Example 1).
procedure, implemented with the three LMI sufficient conditions, to identify a not robustly D -stable system for the case of the disk region. Fig. 3 shows the uncertain domain partition when considering [6, Th. 4]. The arrow indicates the first polytope coordinate found where the system is not D -stable. The uncertain domain is a tetrahedron, in the three-dimensional space, parameterized by the first three polytope coordinates since 4 = 1 0 3i=1 i . 2) Example 2: Consider the state–space model of the linearized vertical-plane dynamics of an aircraft (example (AC1) in [16]), with five uncertain parameters included here
0
dx dt
0
=
0
0
0:0485
0
0
p1
0 0
0
00:12
0
0
1
0
0
0
0
+
y=
1:132
0
0
00:0538 00:1712
0
p4 p5
0 0
01:665 00:0732
1
0
0
0
0
0
1
0
0
0
0
0
1
0
0
[0:2618; 0:3200]
p
0:0705
1
0
09:8556 01:013 p2 0 p3
x
disk centered in the origin with radius r = 3:6, and the conic sector with internal angle 2 = 2=4. The controller matrices are given by
u
= [ p1 ; . . . ; p 5 ]
D
01
0
x:
The uncertain parameter vector perrectangle
TABLE IV PERFORMANCE OF THE PROPOSED ROBUSTLY -STABILITY ANALYSIS APPROACH (EXAMPLE 2)
varies in the hy-
2 [0:9479; 1:1585] 2 [0:6173; 0:7545] 2[3:9771; 4:8609] 2 [1:4175; 1:7325]
with 32 vertices in the five-dimensional space. Based on the design procedure presented in [15], it was designed a first-order dynamic output controller to robustly place the closed-loop poles in the intersection of the half-plane region Real(z ) < 01:5, the
Ac = 0 5:1437 Bc = [ 5:2419 0:0817 00:2764 ] 05:7436 7:4180 00:1314 0:4392 Cc = 4:7128 Dc = 04:2613 1:3786 6:0419 : 07:2602 17:5314 9:7853 25:8029 In this example, with the formulation based on [4, Th. 1], the LMIC solver generates a memory error for the the initial polytope analysis which does not occur with the SDMi solver, but the latter requires prohibitive computational time in relation to the other LMI formulations. Table IV compares the performance of the proposed analysis approach implemented with the other two LMI analysis formulations. Without the polytope partition, both LMI formulations are not capable to inform if the uncertain system is robustly D -stable. Based on the proposed analysis approach implemented with [6, Th. 4], it is verified that the system is robustly D -stable. Considering [6, Th. 1] (quadratic stability), the proposed procedure reaches the maximum number of iter-
1714
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 51, NO. 10, OCTOBER 2006
ations without identifying the robust stability for the two first regions due to the high conservatism of this sufficient condition. As presented in Table I, in this case, 103 evaluations of the sufficient stability conditions are necessary in the first step of the polytope partition (Delaunay triangulation) and 32 evaluations in each further refinement (simplex subdivision). To verify the three regions simultaneously, the proposed procedure spends 23 subdivisions and 52989 s considering [6, Th. 4] and SDMi solver. One can observe in Table IV the reduction of the complexity (K, L) when working with simplex instead hyperrectangle. Considering the results in Table IV and the memory problem with the analysis formulation based on [4, Th. 1], even for higher dimensional spaces, the proposed analysis approach, with the appropriate choice of the sufficient condition, can be more computationally efficient than LMI-based formulations alone.
VI. CONCLUSION It has been proposed an approach to establish if a linear time-invariant system, with polytope-bounded uncertainty, is robustly D-stable or not. Due to the conservativeness of the LMI-based analysis formulations, when the LMI-solvers do not find a solution to the feasibility problem it is not possible to state that the system is not robustly D -stable. If the system is D -stable, the proposed approach iteratively subdivides the uncertain domain until the LMI sufficient condition is verified for all subpolytopes. In the contrary case, the proposed approach terminates when any system related to a subpolytope vertex is verified to be not robustly D -stable. In spite of the implementation complexity, the proposed approach can spend less computational time than LMI-based formulations with a large number of decision variables and constraints that are more influenced by the number of polytope vertices and system order. The main difficulty in the implementation of the proposed approach is the polytope partition. To overcome this problem, it has been proposed a new approach that combines the Delaunay triangulation, to decompose a polytope with any shape in a set of simplices, and a new simplex edgewise subdivision algorithm.
REFERENCES [1] M. C. de Oliveira, J. Bernussou, and J. C. Geromel, “A new discretetime robust stability condition,” Syst. Control Lett., vol. 37, no. 4, pp. 261–265, 1999. [2] D. C. W. Ramos and L. D. Peres, “A less conservative LMI condition for the robust stability of discrete-time uncertain systems,” Syst. Control Lett., vol. 43, pp. 371–378, Apr. 2001. [3] D. C. W. Ramos and P. L. D. Peres, “An LMI condition for the robust stability of uncertain continuous-time linear systems,” IEEE Trans. Autom. Control, vol. 47, no. 4, pp. 675–678, Apr. 2002. [4] V. J. S. Leite and P. L. D. Peres, “An improved LMI condition for robust -stability of uncertain polytopic systems,” IEEE Trans. Autom. Control, vol. 48, no. 3, pp. 500–504, Mar. 2003. [5] M. Chilali and P. Gahinet, “ design with pole placement constraints: An LMI approach,” IEEE Trans. Autom. Control, vol. 41, no. 3, pp. 358–367, Mar. 1996. [6] D. Peaucelle, D. Arzelier, O. Bachelier, and J. Bernussou, “A new robust -stability condition for real convex polytopic uncertainty,” Syst. Control Lett., vol. 40, pp. 21–30, 2000. [7] C. DeMarco, V. Balakrishnan, and S. Boyd, “A branch and bound methodology for matrix polytope stability problems arising in power systems,” in Proc. IEEE Conf. Decision and Control, Honolulu, HI, Dec. 1990, pp. 3022–3027. [8] V. Balakrishnan, S. Boyd, and S. Balemi, “Branch and bound algorithm for computing the minimum stability degree of parameter-dependent linear systems,” Int. J. Robust Nonlinear Control, vol. 1, no. 4, pp. 295–317, Oct.–Dec. 1991.
D
H
D
[9] D. W. Moore, Simplicial mesh generation with applications. Ithaca, NY, Dept. Comput. Sci., Cornell Univ., 1992, Ph.D. dissertation. [10] J. Bey, “Simplicial grid refinement: On Freudenthal’s algorithm and the optimal number of congruence classes,” Numerische Mathematik, vol. 85, pp. 1–29, 2000. [11] H. Edelsbrunner and D. R. Grayson, “Edgewise subdivision of a simplex,” Discrete Comput. Geom., vol. 24, pp. 707–719, 2000. [12] P. Gahinet, A. Nemirovski, A. J. Laub, and M. Chilali, LMI Control Toolbox: For Use with MATLAB. Natick, MA: The MathWorks, Inc., 1995. [13] J. Sturm, “Using SeDuMi 1.02, A MATLAB toolbox for optimization over symmetric cones,” Optim. Meth. Software vol. 11–12, pp. 625–653, 1999 [Online]. Available: http://sedumi.mcmaster.ca, ed. 1.1 and static [14] U. Shaked, “An LPD approach to robust output-feedback design,” IEEE Trans. Autom. Control, vol. 48, no. 5, pp. 866–872, May 2003. [15] E. N. Gonçalves, R. M. Palhares, and R. H. C. Takahashi, “Improved = control problem for linear optimisation approach to robust systems,” Proc. Inst. Elect. Eng. Control Theory Appl., vol. 152, no. 2, pp. 171–176, 2005. [16] F. Leibfritz and W. Lipinski, Description of the benchmark examples in COMP l ib 1.0 Univ. Trier, Trier, Germany, 2003, Tech. Rep.
H
H
H H
On Uniform Stabilization of Discrete-Time Linear Parameter-Varying Control Systems Ji-Woong Lee
Abstract—For discrete-time polytopic linear parameter-varying systems, the uniform exponential stability, which is equivalent to the asymptotic stability and includes the quadratic stability as a special case, is characterized by the union of an increasing family of linear matrix inequality conditions. On the other hand, in certain cases, nonconservative synthesis of uniformly stabilizing controllers is achieved via a system of finite number of linear matrix inequalities if the controller is allowed to have finite memory of past parameters. Two such cases are robust state feedback (in a relaxed sense) against polytopic parameter variations, and multiple output injections under polytopic fusion rules. Index Terms—Discrete linear inclusion, linear matrix inequality (LMI), linear parameter-varying (LPV) system, output injection, state feedback.
I. INTRODUCTION A common approach to analysis and control of linear systems whose parameters vary with uncertainty or in a complicated way (possibly depending on the state) is first to approximate the set of coefficients by a convex polytope and then to abstract the system as the family of all linear time-varying systems whose coefficients vary within this polytope. The resulting abstractions are called (polytopic) linear parameter-varying (LPV) systems. An LPV system is said to be stable if
Manuscript received September 9, 2005; revised February 3, 2006. Recommended by Associate Editor Y. Ohta. The author is with the Department of Electrical and Computer Engineering, the University of Florida, Gainesville, FL 32611 USA (e-mail: jiwoong@ufl. edu). Digital Object Identifier 10.1109/TAC.2006.880804
0018-9286/$20.00 © 2006 IEEE
