

                          

! "##$%   ! &    ' (    ) *    ! )    +   ,  ,

Asymptotic Exactness of Parameter-Dependent
Lyapunov Functions:
An Error Bound and Exactness Verification∗
Yasuaki Oishi†This paper provides an approximate approach to a robust semidefinite programming problem with a functional variable and shows its asymptotic exactness. This problem covers a variety of control problems including a robust stability/performance analysis witha parameter-dependent Lyapunov function. In the proposed approach, an approximate
semidefinite programming problem is constructed based on the division of the set of parameter values. This approach is asymptotically exact in the sense that, as the resolutionof the division becomes higher, the optimal value of the constructed approximate problem
converges to that of the original problem. Our convergence analysis is quantitative. In particular, this paper gives an a priori upper bound on the discrepancy between the optimalvalues of the two problems. Moreover, it discusses how to verify that an optimal solutionof the approximate problem is actually optimal also for the original problem. Finally, theresults are generalized for robust stability/performance analysis against a time-varying
parameter and for stability analysis of a nonlinear system.
Keywords: parameter-dependent Lyapunov functions, robust semidefinite programming,
approximation error, exactness verification, matrix dilation, linear matrix inequalities.

1. Introduction
Search of a parameter-dependent Lyapunov function is an important control problem for lessconservative analysis and design of a parameter-dependent system. This problem is, however,
not easy due to its infinite-dimensional nature. A widely accepted approach is to assume
polynomial parameter dependence on a Lyapunov function and to apply recent techniques
based on positive polynomials [14, 5, 19, 21, 22]. As a result, we have a standard semidefinite
programming (SDP) problem, which is solvable with the efficient interior-point method. Thisapproach is conservative by two reasons. First, an infinite-dimensional problem to find an
unknown function is reduced to a finite-dimensional one. Second, a semi-infinite constraintnecessary to hold for infinitely many parameter values is replaced by a finitely many linear
matrix inequality (LMI) constraints.
Bliman [2] showed asymptotic exactness of a polynomially parameter-dependent Lyapunov
function. That is, there is no conservatism if we assume polynomial parameter dependence∗†

July 28, 2007
Department of Information Systems and Mathematical Sciences, Nanzan University, Seireicho 27, Seto

489-0863, Japan; Phone: +81-561-89-2000; Fax: +81-561-89-2082; Email: oishi@nanzan-u.ac.jp

of sufficiently high degree. This result combined with asymptotic exactness of the positivepolynomial techniques [10, 15, 20, 3, 24, 25] implies asymptotic exactness of the above approach.
This analysis is, however, only qualitative. Namely, it states that conservatism asymptoticallyvanishes as the degrees of the Lyapunov function and the associated positive polynomial increase. It does not tell how large degrees are required to suppress conservatism to a certain
level. This theoretical limitation prevents us from efficient use of the above approach to suppress
conservatism with small computational complexity.
In this paper, we give an approximate approach to a parameter-dependent Lyapunov function and show its asymptotic exactness in a quantitative manner. We generally formulate theproblem as a robust SDP problem with a functional variable. This problem covers a wide
range of problems such as robust stability/performance analysis and design of a gain-scheduled
state-feedback controller. The proposed approach is based on division of the parameter set. We
make the infinite-dimensional problem to a finite-dimensional one by considering a piecewise
polynomial consistent with the division. We make the semi-infinite constraint to a finite number of LMI constraints using the matrix-dilation technique in [17, 18]. Thus, a standard SDP
problem is constructed to approximate the original robust SDP problem. As the resolution of
the division goes up, the optimal value of the approximate problem converges to that of theoriginal problem. In particular, we give an a priori upper bound on the approximation error,that is, the discrepancy between the two optimal values. This bound shows a quantitative
relationship between the approximation error and the resolution of the division. This boundis also useful to make an efficient division, which gives good approximation with low computational complexity. Moreover, we discuss how to verify exactness of the approximation. With
the proposed method, we can detect that the obtained approximate result is actually exact andhigher resolution of the division is not necessary any more. At the end of the paper, the results
are generalized for robust stability/performance analysis against a time-varying parameter andfor stability analysis of a nonlinear system. There, the considered robust SDP problem includes
not only an unknown function but also its derivatives.
These results are extension of the results of the present author [17, 18], which is on a robust
SDP problem without a functional variable. It is notable that the extension is natural in spiteof apparently large difference due to the existence of a functional variable. Many authors haveconsidered the use of piecewise Lyapunov functions [4, 16, 9, 13, 23, 11, 6]. Although some ofthem discussed asymptotic exactness of their methods [9, 13, 23], their discussion was limited
to qualitative one. Coutinho–Danès [6] gave an approach close to the present one without
discussion on asymptotic exactness. Scherer [24] discussed verification of exactness in the caseof no functional variable.
Construction of the paper is as follows. After a problem is formulated in Section 2, our2approximate approach is presented in Section 3. Section 4 gives an upper bound on the approximation error and shows asymptotic exactness of our approach. Section 5 discusses verificationof exactness. Section 6 provides a numerical example. After the results are generalized in
Section 7, the paper is concluded in Section 8.
The symbol Rp stands for the set of p-dimensional real vectors while Zp+ for the set of pdimensional vectors of nonnegative integers. The symbol T denotes the transpose of a matrix
or a vector. For θ = [θ1 θ2 · · · θp ]T ∈ Rp and α = [α1 α2 · · · αp ]T ∈ Zp+ , the symbol θααmeans the product θ1α1 θ2α2 · · · θp p . The symbols Oq×r and Iq designate the q × r zero matrix andthe q × q identity matrix, respectively. The sizes of these matrices are omitted when they are
obvious from the context. The maximum singular value of a matrix A is written as σ(A). Fora real symmetric matrix A, the inequality A º O means that A is positive semidefinite, thatis, xT Ax is nonnegative for any real vector x. Similarly, A Â O expresses that A is positive
definite. For two real symmetric matrices A and B, the inequality A º B means A − B º O.
The Kronecker product of two (not necessarily symmetric) matrices A = (aij ) and B is definedas
a11 B · · · a1r B
 .
.. 
.
A ⊗ B := 
. 
 .
.
aq1 B · · · aqr B

For a set S, the symbol |S| denotes its cardinality.

2. Problem
The problem to be considered in this paper is the following:P : minimize

cT x

subject to E(x) º O,

F(x, φ(θ); θ) º O(∀θ ∈ Θ),where the optimization variables are a vector x ∈ Rn and a function φ ∈ Φ. The set Φ is alinear space of piecewise continuous functions mapping Θ to Rnφ . The set Θ is a given closed
polytope in Rp and works as the domain of an uncertain p-dimensional parameter θ. The
objective function is a linear function of x only, which is described as cT x with given c ∈ Rn . In
the constraints, the value of E(x) is an ` × ` symmetric matrix while the value of F(x, φ(θ); θ)is an m × m symmetric matrix. The function E(x) is affine in x ∈ Rn . The function F(x, a; θ)is affine in x ∈ Rn and a ∈ Rnφ while polynomial in θ. The optimal value of the problem P isdenoted by inf P with the attention that the minimum may not be attained.
The problem P is called a robust SDP problem because the second constraint has to besatisfied for all possible values of the uncertain parameter θ. It is also notable that P has a3

functional variable φ ∈ Φ. Many control problems are formulated into the form of P . Here isa simple example.
Example 1. Consider a parameter-dependent linear system
˙ = A(θ)ξ(t) + B(θ)u(t),ξ(t)y(t) = C(θ)ξ(t),where the matrices A(θ), B(θ), and C(θ) are polynomials of a time-invariant but unknown
parameter θ. Let the domain of θ be a closed polytope Θ. Suppose that A(θ) is Hurwitz forall θ ∈ Θ. Then, the maximum L2 -induced norm over Θ is computed through the problem:
minimizex−A(θ)R(θ) − R(θ)A(θ)T −R(θ)C(θ)T −B(θ)
ºO
subject to −C(θ)R(θ)
xI
O


T
−B(θ)
O
xI

(∀θ ∈ Θ),where the optimization variables are a real scalar x and a piecewise continuous symmetricmatrix-valued function R(θ). No generality is lost by the piecewise continuity of R(θ) because
A(θ) is continuous and Θ is compact [2]. This problem is in the form of P .

¤Various generalization is possible with this example. The matrices A(θ), B(θ), and C(θ)
can be rational functions of θ. In this case, multiplication of an appropriate polynomial makes
the constraint polynomial in θ. Design of a gain-scheduled state-feedback controller has asimilar structure to Example 1 and can be formulated into P . When the system depends on atime-varying parameter, we need to modify the problem P . This will be discussed in Section 7.
The problem P is difficult to solve by two reasons: it has a functional variable φ ∈ Φ andthe semi-infinite constraint F(x, φ(θ); θ) º O (∀θ ∈ Θ). In the succeeding sections, we propose
an approximate approach and discuss its properties.

3. Matrix-dilation approach
An approximate approach to the problem P is presented in this section. In the approach, we
require φ to be a low-order polynomial and make the problem finite-dimensional. In order toimprove the quality of approximation, we divide the parameter set Θ into several subregions
and allow φ to be a piecewise polynomial. To deal with the semi-infinite constraint, we employ
the matrix-dilation approach in [17, 18], which is again based on the division of the parameter
set. Hence, we use the division for two purposes: to approximate the functional variable
4and to approximate the semi-infinite constraint. As the resolution of the division goes up,
approximation is improved in the two senses. We will quantitatively investigate this property
in the next section.
For the function variable φ(θ), we use a fixed-order polynomial
set S ⊂ Zp+ .
φu (θ) =PPα∈S uα θ

αfor some finite

We use the coefficients u = (uα ) ∈ Rnu to characterize the polynomial and writeα∈S uα θα. Here, the dimension nu is equal to |S|nφ . Substitution of φu (θ) into F

makes F dependent on the finite-dimensional variables x and u. This is made explicit by thenotationF (x, u; θ) := F(x, φu (θ); θ).
Note that F is affine in x and u while polynomial in θ.
We next introduce a division ∆ of the parameter set Θ. This is a finite collection of pdimensional closed convex polytopes {Θ[1] , Θ[2] , . . . , Θ[J] } such that the intersection Θ[j] ∩ Θ[k]has no interior point for any j 6= k and the union ∪Jj=1 Θ[j] is equal to the whole parameter set
Θ. We allow u to take a different value u[j] depending on the subregion Θ[j] (j = 1, 2, . . . , J).
Hence, our φ is a piecewise polynomial. We consider the following approximate problem:
P0 (∆) : minimize

cT x

subject to E(x) º O,

F (x, u[j] ; θ) º O(∀θ ∈ Θ[j] , ∀j = 1, 2, . . . , J),

where the optimization variables are x ∈ Rn and u[1] , u[2] , . . . , u[J] ∈ Rnu . Since only a limitedclass of functions are considered for φ, we have inf P ≤ inf P0 (∆).
The problem P0 (∆) is still difficult to solve because it has a semi-infinite constraint. We
use the matrix-dilation approach to circumvent this difficulty. In particular, we replace thesemi-infinite constraint with its sufficient condition that can be expressed by a finite numberof LMIs. The resulting problem P (∆) is a standard SDP problem and thus solvable with
the interior-point method. We have inf P ≤ inf P0 (∆) ≤ inf P (∆) and expect that inf P (∆)
converges to inf P as the resolution of the division ∆ becomes higher.
For construction of the problem P (∆), we use the expansion F (x, u; θ) =Pα Fα (x, u)θ

α

. In

many practical situations, only a part of the coefficient matrices Fα (x, u) is nonzero. We would
like to exploit this sparse structure to construct P (∆) of small size. To this end, we consider adirected graph (V, A) having the following properties.
(i) The vertex set V is a finite subset of Zp+ . It consists of the origin, all α’s with nonzero
coefficient matrices Fα (x, u), and possibly some additional points.
(ii) The arc set A is a finite set of arcs (α, β) with α, β ∈ V such that βi = αi + 1 for one andonly one i = 1, 2, . . . , p and βi = αi for all the remaining i’s.

5

α2

6
w

α

(4)

6
s -

w-

α

(1)

α

(2)

α

w (3)

α1

Figure 1. A rectilinear Steiner arborescence in the case of p = 2 and {α | Fα (x, u) 6= O} ={[0 0]T , [2 0]T , [1 1]T }.
(iii) For any vertex α ∈ V , the directed graph (V, A) has one and only one directed path fromthe origin to α.
Figure 1 shows such a graph (V, A) in the case of p = 2 and {α | Fα (x, u) 6= O} = {[0 0]T , [2 0]T ,[1 1]T }. Note that the vertex [1 0]T is added for satisfaction of the conditions (ii) and(iii). A graph (V, A) satisfying (i)–(iii) is called a rectilinear Steiner arborescence for the set{α | Fα (x, u) 6= O}. When the set {α | Fα (x, u) 6= O} has small cardinality, we can construct a
rectilinear Steiner arborescence with small |V |.
We number the elements of V as α(1) , α(2) , . . . , α(|V |) . The numbering is arbitrary as far asα(1) = 0. We here define the matrices:"#2Fα(1) (x, u) F∗ (x, u)G(x, u) :=,F∗ (x, u)TO(1)F∗ (x, u) := [Fα(2) (x, u) Fα(3) (x, u) · · · Fα(|V |) (x, u)],
(1)

(2)

M (θ) := [θα Im θα Im · · · θα(|V |)Im ]T .

It is easy to see 2F (x, u; θ) = M (θ)T G(x, u)M (θ). We next define the matrix H(θ) bye
H(θ) := H(θ)
⊗ Im

(2)

e
with the (q, r)-element of H(θ)
being


−θ , if (α(q) , α(r+1) ) ∈ A and is parallel to the ith axis;


 i
e
(H(θ))
1,
if q = r + 1;
qr :=



0,
otherwisefor q = 1, 2, . . . , |V | and r = 1, 2, . . . , |V | − 1. It is straightforward to see M (θ)T H(θ) = O
and the matrix [M (θ) H(θ)] is square and nonsingular for any θ. Hence, H(θ) is called anorthogonal complement of M (θ). It is notable that this particular orthogonal complement is6

affine in θ. In the case of Figure 1, the matrices M (θ) and H(θ) look as−θ1 Im
Im




 Im
 θ1 Im −θ1 Im −θ2 Im 
.



,
H(θ) = M (θ) =  2


Im


 θ1 Im Imθ1 θ2 ImWith these preparations, we introduce the desired problem:
P (∆) : minimize

cT x

subject to E(x) º O,

G(x, u[j] ) + H(θ)(W [j] )T + W [j] H(θ)T º O(∀θ ∈ ver Θ[j] , ∀j = 1, 2, . . . , J),
where the optimization variables are x ∈ Rn , u[1] , u[2] , . . . , u[J] ∈ Rnu , and |V |m × (|V | − 1)mmatrices W [1] , W [2] , . . . , W [J] . The symbol ver Θ[j] denotes the set of the vertices of the convex
polytope Θ[j] . Since the second constraint is affine in θ, it actually holds for any θ ∈ Θ[j] . Then,
premultiplication of M (θ)T and postmultiplication of M (θ) give F (x, u[j] ; θ) º O. This leads
to the next proposition [17, 18].
Proposition 2. The feasible region of P (∆) projected onto the space of (x, u[1] , u[2] , . . . , u[J] )is contained in the feasible region of P0 (∆). In particular, inf P ≤ inf P0 (∆) ≤ inf P (∆).
We here have an approximate approach to P , which we call a matrix-dilation approach. Thatis, we first choose a division ∆ for the parameter set Θ and solve the approximate problem
P (∆). If the obtained optimal value inf P (∆) is satisfactory, we stop. Otherwise, we subdivide∆ and solve again the new approximate problem P (∆).
This is a natural approach to the problem P . However, it is not clear how inf P (∆) converges
to inf P . We quantitatively answer to this question in the next section.

4. An upper bound on the approximation error
4.1. Upper boundIn this section, we obtain an a priori upper bound on the approximation error | inf P (∆)−inf P |in terms of the resolution of the division. This is a generalization of the upper bounds in [17, 18],
which are for a robust SDP problem without a functional variable.
We make the following assumptions to have the result. The first assumption is a reasonable
one, which means strict feasibility of the problem P . The second assumption means thatsearch of bounded x and φ is sufficient for the problem P as well as its modification P² , to be7

introduced in Section 4.2. This assumption is mild enough because it is satisfied when P and
P² have bounded optimal solutions. The last assumption is again mild one, which means thatthe polynomial φu (θ) has a constant term.
Assumption 3.
(a) There exist x ∈ Rn and φ ∈ Φ such that E(x) Â O and F(x, φ(θ); θ) Â O (∀θ ∈ Θ).
(b) There exist three positive numbers ², x, and φ such that, for any 0 ≤ ² ≤ ² and any v ∈ R,
the set {(x, φ) ∈ Rn × Φ | cT x ≤ v, E(x) º O, F(x, φ(θ); θ) º ²I (∀θ ∈ Θ)} is either empty
or having an element (x, φ) with kxk < x and kφ(θ)k < φ (∀θ ∈ Θ).
(c) The set S in the polynomial φu (θ) =Pα∈S uα θ

αincludes the origin.

¤

We next need a measure of the resolution of the division. For a division ∆ = {Θ[1] , Θ[2] , . . . , Θ[J] }of the parameter set Θ, the radius of a subregion Θ[j] is rad Θ[j] := minθ∈Θ[j] maxθ0 ∈Θ[j] maxi=1,2,...,p|θi0 − θi |. A θ ∈ Θ[j] that attains the minimum is called a center of Θ[j] . The maximum radiusof the division ∆ is defined as rad ∆ := maxj=1,2,...,J rad Θ[j] , which measures the resolution of∆.
With these preparations, we have the desired upper bound. The proof will be given in
Section 4.2.
Theorem 4. Under Assumption 3, there exist positive numbers C and C 0 , which are independent of the division ∆, such that| inf P (∆) − inf P | ≤ C rad ∆for any ∆ satisfying rad ∆ ≤ C 0 .
The concrete forms of C and C 0 will be given in (6) and (7).
This theorem immediately implies the asymptotic exactness of our approach. To see this,consider a sequence of divisions whose maximum radii converge to zero. Then, the corresponding approximate problems P (∆) have the approximation errors converging to zero because C
and C 0 are independent of ∆. Note that Theorem 4 also gives the rate of convergence. Namely,
the convergence is at least in the linear order of the maximum radius.
Theorem 4 shows a tradeoff between the approximation error and the computational complexity of the present approach. Indeed, the theorem tells that the maximum radius should besmall for small approximation error. A small maximum radius leads to a large number of subregions and, then, to a large number of variables and constraints in the approximate problem
P (∆). Since the number of subregions is in the order of (rad ∆)−p , a small maximum radius
invites computational difficulty especially when p is large.
8In order to address this issue, we can improve the theorem by replacing the maximum
radius with a more sophisticated measure of the resolution. Consider an approximate problem
P (∆) for some division ∆ and suppose that its optimal value is attained. The constraint
G(x, u[j] ) + H(θ)(W [j] )T + W [j] H(θ)T º O is called active for an optimal solution if the lefthand-side matrix has a zero eigenvalue for this optimal solution. A subregion Θ[j] is called
active if it has an active constraint. We define the maximum active radius of the division∆, denoted by a-rad ∆, as the maximum radius over all active subregions. When P (∆) has
multiple optimal solutions, we define a-rad ∆ by taking the minimum of all possible values. In
this case, actual computation of a-rad ∆ may be difficult while computation of its upper boundis easy. When the optimal value is not attained in P (∆), we define a-rad ∆ := rad ∆. We here
have the following result. The proof is given in Section 4.2.
Corollary 5. Under Assumption 3, the inequality| inf P (∆) − inf P | ≤ C a-rad ∆holds for any ∆ satisfying rad ∆ ≤ C 0 , where C and C 0 are the same positive constants as in
Theorem 4.
This corollary states that we do not need to decrease the maximum radius but the maximum
active radius for the reduction of the approximation error. Although direct decrease of the
maximum active radius is difficult, a heuristic algorithm can be constructed.
Algorithm 6.
0. Choose a coarse division ∆.
1. Solve P (∆) for the current division ∆.
2. Stop if the obtained solution is satisfactory.
3. If an optimal solution is found, find an active subregion that attains a-rad ∆ and divide it
into two subregions so that they have small radii. Otherwise, find a subregion that attains
rad ∆ and divide it similarly.
4. Go back to Step 1 with the updated division ∆.

¤This type of algorithm was originally proposed for a robust SDP problem without a functional variable [17].

94.2. Proof
We prove Theorem 4 as well as Corollary 5 in this subsection. We present the full proof for thesake of completeness though it includes some discussion found also in [17, 18]. The key idea isto relate the following auxiliary problem with the approximate problem P (∆):
P² : minimize

cT x

subject to E(x) º O,

F(x, φ(θ); θ) º ²I(∀θ ∈ Θ),where ² is a nonnegative number. We will find below an ² that attains inf P ≤ inf P (∆) ≤ inf P² .
Note that inf P² is convex in ² and that P0 is identical with the original problem P . With g being
an upper bound on a subgradient of inf P² , we have | inf P² − inf P | = inf P² − inf P ≤ g². Since
this implies | inf P (∆) − inf P | ≤ g², we obtain the desired upper bound on the approximation
error.
We first prepare a |V |m × |V |m matrix L(θ), which will be used for simplification of thedilated LMI constraint G(x, u[j] ) + H(θ)(W [j] )T + W [j] H(θ)T º O. The definition of L(θ) isbased on the rectilinear Steiner arborescence (V, A). For α(q) , α(r) ∈ V , we say that α(q) is
reachable from α(r) in (V, A) if these two vertices are identical or (V, A) has a directed pathconnecting α(r) to α(q) . For the convenience of the proof, we assume that the numbering of
the vertices is consistent with the partial order defined by (V, A). In other words, we haveq ≥ r whenever α(q) is reachable from α(r) . This assumption does not contradict with α(1) = 0
because the origin is the root of the arborescence (V, A). Also, this assumption does not harm
the generality of the result because the result does not depend on the numbering.
Now, the matrix L(θ) is defined ase ⊗ ImL(θ) := L(θ)
with

(3)


θα(q) −α(r) , if α(q) is reachable from α(r) in (V, A);e
L(θ)qr =0,otherwise.

ee qr 6= 0 only if q ≥ r. Moreover, its diagonal
The matrix L(θ)
is lower triangular, i.e., L(θ)
eelements are all equal to unity. A consequence is nonsingularity of L(θ)
and also of L(θ).
Straightforward calculation gives the concrete forms of L(θ)T G(x, u)L(θ) and L(θ)T H(θ0 )for θ, θ0 ∈ Rp .
Lemma 7. For the matrices G(x, u) and L(θ) in (1) and (3), respectively, we have"#
2F (x, u; θ) F∗∗ (x, u; θ)TL(θ) G(x, u)L(θ) =F∗∗ (x, u; θ)T
O
10

with
F∗∗ (x, u; θ) =h X

Fα (x, u)θα−αα∈V (2)(2)X

Fα (x, u)θα−α

(3)

···

α∈V (3)X

Fα (x, u)θα−α(|V |)i,α∈V (|V |)where V (r) is the set of vertices reachable from α(r) in (V, A) for r = 2, 3, . . . , |V |.
Lemma 8. For the matrices H(θ) and L(θ) in (2) and (3), respectively, we can write∗ ∗ ··· ∗ 1 ∗ · · · ∗




T
0
L(θ) H(θ ) = 0 1 · · · ∗ ⊗ Im ,. . .
 .. .. . . ... 
0 0 ··· 1
where an element expressed by ∗ is either equal to zero or expressed as θα (θi − θi0 ) for someα ∈ Zp+ and i = 1, 2, . . . , p. When θ = θ0 in particular, the elements expressed by ∗ are all equal
to zero.
With these preparations, we relate P² and P (∆). By Assumption 3 (a), there exists ²0 > 0such that, for any 0 ≤ ² ≤ ²0 , the auxiliary problem P² is strictly feasible. This ²0 can be chosen
smaller than or equal to the ² in Assumption 3 (b). We consider the set of all u such that thecorresponding φu (θ) is constantly equal to some vector a with kak < φ. This set is nonempty
by Assumption 3 (c). We let F be an upper bound of σ[F∗∗ (x, u; θ)] over all kxk < x, all θ ∈ Θ,
and all u in the considered set.
We begin by the special case that Θ ⊆ [−1, 1]p . In this case, |θi | ≤ 1 for any θ ∈ Θ and any
i = 1, 2, . . . , p.
Lemma 9. Suppose that Θ ⊆ [−1, 1]p andnrad ∆ ≤ min

2²
1 o
p0
,
.
(F + |V |m)2 |V |Then, we have inf P ≤ inf P (∆) ≤ inf P² forp(F + |V |m)2²=
rad ∆.
2
Proof. Since inf P ≤ inf P (∆) by Proposition 2, we will show inf P (∆) ≤ inf P² . Let ∆ be[1]

{Θ , Θ[2] , . . . , Θ[J] }. The ² given in the lemma satisfies 0 ≤ ² ≤ ²0 . Let (x, φ) be a feasible
solution of P² satisfying kxk ≤ x and kφ(θ)k ≤ φ (∀θ ∈ Θ). Here, F(x, φ(θ); θ) º ²I for anyθ ∈ Θ. We show below that there exists (u[j] , W [j] ) for each j = 1, 2, . . . , J such that
G(x, u[j] ) + H(θ)(W [j] )T + W [j] H(θ)T º O
11(∀θ ∈ ver Θ[j] ),(4)

in other words, this x is feasible for P (∆). Then the proof is complete because (x, φ) can bechosen so that cT x is arbitrarily close to inf P² . In fact, the inequality (4) holds if we choose(u[j] , W [j] ) for each j so that φu[j] (θ) ≡ φ(θc ) and W [j] = (1/ rad Θ[j] )H(θc ) , where θc is a centerof Θ[j] . Note in this case that σ[F∗∗ (x, u[j] ; θ)] ≤ F for any j and θ ∈ Θ.
In order to show the desired inequality (4), we premultiply L(θc )T and postmultiply L(θc ) toit. Lemmas 7 and 8 give the concrete forms of L(θc )T G(x, u[j] )L(θc ) and (rad Θ[j] )(W [j] )T L(θc ) =H(θc )T L(θc ). By Lemma 8 again, the product (1/ rad Θ[j] )L(θc )T H(θ) has the form


∗
∗
···
∗


1/ rad Θ[j]

∗
·
·
·
∗




0
1/ rad Θ[j] · · ·∗ ⊗ Im ,


.
.
.
.


..
..
..
..


[j]
0
0
· · · 1/ rad Θ
where an element expressed by ∗ is either equal to zero or of the form(θc )α (θic − θi ),rad Θ[j]whose magnitude is at most one since |θic | ≤ 1 and |θic − θi | ≤ rad Θ[j] for i = 1, 2, . . . , p. Let us
write the product (1/ rad Θ[j] )L(θc )T H(θ) as [H1T H2T ]T with the m × (|V | − 1)m matrix H1and the (|V | − 1)m × (|V | − 1)m matrix H2 . Then, we have
"
#
O
H
m×m
1
L(θc )T H(θ)(W [j] )T L(θc ) =.
O(|V |−1)m×m H2
Since H1 has at most (|V | − 1)m nonzero elements whose magnitude is at most one, we havepσ(H1 ) ≤ (|V | − 1)m. On the other hand, H2 is upper triangular and each of its rows has the
diagonal element 1/ rad Θ[j] and at most |V |−2 nonzero off-diagonal elements, whose magnitude
is at most one. Hence, we have
H2 + H2T º

³

´
2
−
|V
|
+
2
I.
rad Θ[j]Now, the left-hand side of (4) multiplied by L(θc )T and L(θc ) has the upper-left m × m block
equal to2F (x, u[j] ; θc ) = 2F(x, φu[j] (θc ); θc ) = 2F(x, φ(θc ); θc ) º 2²I.
Its Schur complement is
H2 + H2T − [F∗∗ (x, u[j] ; θc )+H1 ]T [2F (x, u[j] ; θc )]−1 [F∗∗ (x, u[j] ; θc )+H1 ]
³ 2
´
p
ª2 1
©
[j] c
º
σ[F
(x,
u
;
θ
)]
+
(|V
|
−
1)m
I.
−
|V
|
+
2
I
−
∗∗
rad Θ[j]2²12Noting that 2/ rad Θ[j] − |V | + 2 > 1/ rad Θ[j] and σ[F∗∗ (x, u[j] ; θc )] ≤ F , we see the positive
semidefiniteness of the right-hand side matrix. This completes the proof.

¤

The general case that not necessarily Θ ⊆ [−1, 1]p can be reduced to the special case. Let
us writeθ := max{1, max max |θi |}(5)θ∈Θ i=1,2,...,pand |α| := α1 + α2 + · · · + αp . Since
F (x, u; θ) =XFα (x, u)θα =α∈V

we can regard Fα (x, u)θ

|α|XFα (x, u)θ

|α|³ θ ´αθα∈V,

as a coefficient and θ/θ as a parameter. Since θ/θ moves in [−1, 1]p ,the discussion in the special case can be applied. To present the result in this case, we define0F∗∗(x, u; θ) by
0F∗∗(x, u; θ) :=h X

Fα (x, u)θ

α∈V (2)|α|³ θ ´α−α(2)θXFα (x, u)θ

α∈V (3)···|α|³ θ ´α−α(3)θXFα (x, u)θ

α∈V (|V |)|α|³ θ ´α−α(|V |) iθ.

0(x, u; θ)] over all kxk < x, all θ ∈ Θ, and all u
We also define F 0 as an upper bound of σ[F∗∗corresponding to the polynomial φu (θ) constantly equal to a vector a with kak < φ.
Lemma 10. Suppose thatnrad ∆ ≤ min

2θ²0
θ o
p
,
.
(F 0 + |V |m)2 |V |Then, we have inf P ≤ inf P (∆) ≤ inf P² forp(F 0 + |V |m)2²=
rad ∆.
2θ
Proof. The replacement of Fα (x, u) by Fα (x, u)θ

|α|and θ by θ/θ does not change theproblems P and P² essentially. As is shown below, this replacement does not change either theapproximate problem P (∆). The parameter after the replacement, i.e., θ/θ, moves in [−1, 1]p .
Hence, the result of Lemma 9 is valid with F replaced by F 0 and rad ∆ by rad ∆/θ. This
completes the proof.
We show that the approximate problem P (∆) does not change by the replacement above.
Let G0 (x, u) and H 0 (θ) be the matrices obtained from G(x, u) and H(θ), respectively, by thisreplacement. It is routine to confirm that
G0 (x, u) = diag{Im , T }G(x, u) diag{Im , T },13

H 0 (θ) = diag{Im , T }H(θ)T −1 ,where T := diag{θ|α(2) |, θ|α(3) |, ..., θ|α(|V |) |} ⊗ Im and diag denotes a block-diagonal matrix.

Therefore, the existence of W satisfying
G(x, u) + H(θ)W T + W H(θ)T º O
is equivalent to the existence of W 0 satisfying
G0 (x, u) + H 0 (θ)(W 0 )T + W 0 H 0 (θ)T º O
with the correspondence W 0 = diag{Im , T }W T . This means that the approximate problem
P (∆) does not change essentially by the replacement.

¤

We now take the final step toward the proof of Theorem 4. Recall that we assume Assumption 3. The number ²0 is such that the auxiliary problem P² is strictly feasible for any
0 ≤ ² ≤ ²0 . It also satisfies ²0 ≤ ² for the ² in Assumption 3 (b). The number θ is as in (5)while F 0 is as in the same paragraph. Finally, let g be an upper bound on the left derivative ofinf P² at ² = ²0 . Then, withC=g(F 0 +

C 0 = minnp

|V |m)2

,

(6)

θ o
2θ²0
p
,
,
(F 0 + |V |m)2 |V |

(7)

2θ

we can prove the theorem.
Proof of Theorem 4. Lemma 10 implies that, when rad ∆ ≤ C 0 , we have inf P ≤ inf P (∆) ≤p
inf P² for ² = [(F 0 + |V |m)2 /2θ] rad ∆. Owing to the convexity of inf P² , the upper bound g
is greater than or equal to the left derivative of inf P² at this ². Hence, the convexity of inf P²
again implies inf P ≥ inf P² − g², from which inf P² − inf P ≤ g². Substitution of the concrete
form of ² gives the theorem.

¤

The proof of Corollary 5 follows from subdivision of ∆.
Proof of Corollary 5. The corollary is reduced to Theorem 4 when the optimal value ofP (∆) is not attained. Hence, we assume that the optimal value is attained.
Consider an optimal solution and subdivide each inactive subregion of ∆, if necessary, sothat each of the created subregion has the radius smaller than or equal to a-rad ∆. Then, thee = {Θ
e[k] } has rad ∆e = a-rad ∆. Note also inf P (∆)
e ≤ inf P (∆). We
resulting new division ∆next consider the SDP dual of P (∆). As is shown below, the dual problem of P (∆) has an
14

optimal solution that attains inf P (∆). Based on this solution, we construct a dual feasible
e In particular, if a subregion Θ
e[k] is a newly created subregion, assign zero
solution of P (∆).
e[k] coincides
matrices to the dual variables corresponding to this subregion; If a subregion Θwith one subregion in ∆, assign the same values as the dual optimal solution of P (∆). With
these assigned values, the dual objective function takes the same value as before, i.e., inf P (∆).
e ≤ inf P (∆). Hence,
By weak duality, we have inf P (∆) ≤ inf P (∆)
e − inf P | ≤ C rad ∆e = C a-rad ∆,| inf P (∆) − inf P | = | inf P (∆)
which shows the claim.

p
We show that the dual problem of P (∆) has an optimal solution. Let ² = [(F 0 + |V |m)2 /2θ]×rad ∆. Since 0 < ² ≤ ²0 ≤ ², the problem P² is strictly feasible and has a solution (x, φ) satisfying kxk < x and kφ(θ)k < φ for any θ ∈ Θ. It is possible to assume strict feasibility of this(x, φ) because it can be perturbed, if necessary, to the direction of a strictly feasible solution.
We construct from this (x, φ) a solution of P (∆) as in the proof of Lemma 9. We then obtain astrictly feasible solution of P (∆). On the other hand, the problem P is bounded from below by
Assumption 3 (b) and so is P (∆). The duality theorem on SDP (Theorem 2.4.1 of [1]) implies
that the dual problem of P (∆) has an optimal solution.

¤

5. Verification of exactness
Suppose that we solved an approximate problem P (∆) for some division ∆ and obtained an
c [j] }). In general, the attained value cT x
optimal solution (bx, {bu[j] }, {Wb = inf P (∆) is larger than
the true optimal value inf P . In some cases, however, these two values happen to be identical.
One of such cases is that there exist some j and θb ∈ Θ[j] such that (bx, ub[j] ) is optimal for the
problem:b : minimize
P (θ)

cT x

subject to E(x) º O,

b º O.
F (x, u[j] ; θ)Indeed, since only a particular parameter value θb ∈ Θ is considered in this problem, the attained
minimum value cT x
b is smaller than or equal to inf P . But we also have inf P ≤ inf P (∆) = cT xb,which implies that these three values are equal to each other. In this case, θb is understood as
the worst-case parameter value, which prevents inf P from becoming smaller than the present
value.
The purpose of this section is to give a condition for the above situation to occur. If thiscondition is satisfied, the obtained solution is optimal not only for the approximate problem
P (∆) but also for the original problem P . This result is useful because, in this case, we
15

do not have to improve the approximate problem any more. The condition is obtained bygeneralization of the result on a robust SDP problem without a functional variable [17]. See
also [24] for a related result.
A feasible solution of an SDP problem is optimal if there exists a dual feasible solution
satisfying the complementary slackness condition with the given (primal) solution (Section 4.1
of [1]). To make an explicit statement, we introduce the inner product hA, Bi := tr AB forsymmetric matrices A and B and let eni be the ith elementary n-dimensional vector (i =b is optimal if there exists a matrix pair
1, 2, . . . , n). Then, a feasible solution (bx, ub[j] ) of P (θ)
b where Yb is an ` × ` symmetric matrix and Zb is an m × m symmetric matrix, that
(Yb , Z),
satisfies the dual feasibility condition
Yb º O,
Zb º O,DE D
Eb − F (0, 0; θ)b = cib F (en , 0; θ)Yb , E(eni ) − E(0) + Z,iDEb − F (0, 0; θ)b = 0 (i = 1, 2, . . . , nu )b F (0, enu ; θ)
Z,
i

(8)
(i = 1, 2, . . . , n),

(9)
(10)

as well as the complementary slackness condition
hYb , E(bx)i = 0,(11)

b = 0.
b F (bhZ,x, ub[j] ; θ)i(12)b is equivalent to the existence of a dual feasible
We show below that the existence of such (Yb , Z)
solution of P (∆) having some special structure. Hence, if we solve the approximate problem
P (∆) and the obtained dual solution has this structure, the obtained solution is exact inthe sense that it optimizes the original problem P . Note that the primal-dual interior-point
method gives not only a primal solution but also a dual solution. Here, the dual variables ofP (∆) consist of an ` × ` symmetric matrix, which corresponds to E(x) º O, and |V |m × |V |msymmetric matrices, each of which corresponds to G(x, u[j] ) + H(θ)(W [j] )T + W [j] H(θ)T º O(θ ∈ ver Θ[j] , j = 1, 2, . . . , J).
c [j] }) be a primal feasible solution of the approximate problem
Theorem 11. Let (bx, {bu[j] }, {Wb has
P (∆) for a division ∆. Let θb be a point in some subregion Θ[j] . Then, the problem P (θ)
b satisfying the complementary slackness condition (11) and (12)a dual feasible solution (Yb , Z)with the primal feasible solution (bx, ub[j] ) if and only if the approximate problem P (∆) has ac [j] })
dual feasible solution satisfying the complementary slackness condition with (bx, {bu[j] }, {Wand having the following structure.
(a) The dual variable corresponding to the constraint E(x) º O is Yb .

16

(b) Each remaining nonzero dual variable corresponds to the constraint G(x, u[j] )+H(θ)(W [j] )T +W [j] H(θ)T º O for some vertex θ ∈ ver Θ[j] .
(c) Let the vertices in (b) be θ(1) , θ(2) , . . . , θ(Q) . Then, the set {θ(1) , θ(2) , . . . , θ(Q) } includes θb in
the relative interior of its convex hull.
b ZM
b T forb (θ)(d) The nonzero dual variable corresponding to θ(q) is expressed as (a(q) /2)M (θ)
each q = 1, 2, . . . , Q, where a(q) is the coefficient of a positive convex combination θb =a(1) θ(1) + a(2) θ(2) + · · · + a(Q) θ(Q) .
b ZM
b T can be checked by solving an SDP problem. Indeed,b (θ)
The structure (a(q) /2)M (θ)
b ZM
b T if and only if AH(θ)
b = O and theb (θ)
a symmetric matrix A is equal to (a(q) /2)M (θ)
b = O is affine in θ,b we canb Since AH(θ)
upper-left m × m submatrix of A is equal to (a(q) /2)Z.
search θb by solving an SDP problem. See Section 6 for an example.
Proof. Assume the existence of a dual feasible solution of P (∆) with the described structure.
Then, we have
Yb º O,D

Zb º O,(13)

Q D (q)
E X
E
a
n
b ZM
b T , G(en , 0) − G(0, 0) = cibb (θ)
Y , E(ei ) − E(0) +
M (θ)
i
2
q=1

(i = 1, 2, . . . , n),(14)

Q D (q)Xa
q=1

2

Q D (q)
X
a
q=1

2

E
nu
T
b
b
b
M (θ)ZM (θ) , G(0, ei ) − G(0, 0) = 0 (i = 1, 2, . . . , nu ),

(15)

E
b ZM
b T , H(θ(q) )W T + W H(θ(q) )T = 0
b (θ)
M (θ)
(for any |V |m × (|V | − 1)m matrix W ), (16)hYb , E(bx)i = 0,(17)
D a(q)
E
b ZM
b T , G(bb (θ)
c [j] )T + W
c [j] H(θ(q) )T = 0 (q = 1, 2, . . . , Q). (18)
M (θ)
x, u
b[j] ) + H(θ(q) )(W
2Here, (13)–(16) constitute the dual feasibility condition while (17) and (18) the complementaryslackness condition.
b T G(x, u)M (θ)
b =We consider (14) first. Application of a(1) + a(2) + · · · + a(Q) = 1 and M (θ)
b gives (9). Similarly, (15) implies (10). Summing up (18) for q = 1, 2, . . . , Q and
2F (x, u; θ)b we have (12). Hence, (Yb , Z)b satisfies (8)–(12),using a(1) θ(1) + a(2) θ(2) + · · · + a(Q) θ(Q) = θ,b satisfying the complementary slackness
in other words, it is a dual feasible solution of P (θ)
condition.
17

b satisfies (8)–(12). By the reversed reasoning, (14) and (15)Conversely, assume that (Yb , Z)
follow from (9) and (10), respectively. The equality (16) is a consequence of a(1) θ(1) + a(2) θ(2) +b T H(θ)
b = O. We finally derive (18) from (12). The equality (12)· · · + a(Q) θ(Q) = θb and M (θ)
implies
Q D (q)
X
a
q=1

2

E
b ZM
b T , G(bb (θ)
c [j] )T + W
c [j] H(θ(q) )T = 0.
M (θ)
x, u
b[j] ) + H(θ(q) )(W

Since each term is nonnegative, the equality (18) follows.

¤

6. Example
A numerical example is presented for illustration.
Consider a two-dimensional system
#
"
#
"
0
1
0
˙ =
u(t),
ξ(t)
ξ(t) +−1 −1 − 0.1θ1
1 + 0.2θ1 − 0.1θ12 − 0.1θ2y(t) = [1 0]ξ(t)depending on a two-dimensional parameter θ = [θ1 θ2 ]T ∈ Θ = [0, 1]2 . We compute its
maximum L2 -induced norm over Θ. This problem is formulated into a robust SDP problem as
in Example 1. Substituting a polynomial R(θ) = R00 +R10 θ1 +R20 θ12 +R01 θ2 into the constraint
there, we have an LMI with the terms of θ1 , θ12 , θ13 , θ2 , and θ1 θ2 as well as the constant term.
With the coarsest division ∆ = {Θ}, we constructed the approximate problem P (∆) andsolved it. The obtained optimal value was 1.20687. We used SeDuMi [26] for the SDP solver
with the help of the parser YALMIP [12] on a computer equipped with Pentium M of 1.10 GHz
and memory of 760 MBytes. The computational time was 7.25 seconds.
The region Θ has four vertices: θ(1) = [0 0]T , θ(2) = [1 0]T , θ(3) = [0 1]T , and θ(4) = [1 1]T .
Among them, θ(1) and θ(2) are interesting because the obtained dual variables corresponding tothese vertices are found nonzero. With the dual variable corresponding to θ(1) , we computed
the worst-case parameter θb as in Section 5. The result was θb = [0.694205 0]T . With thisparameter, the system has the L2 -induced norm equal to 1.20687, which shows the accuracyof our approximate optimal value. Moreover, the obtained dual variables corresponding to θ(1)b ZM
b T and (a(2) /2)M (θ)
b ZM
b T , respectively, up to the
b (θ)
b (θ)
and θ(2) coincide with (a(1) /2)M (θ)
b This verifies
order of 10−8 for a(1) = 1 − 0.694205, a(2) = 0.694205 and appropriately chosen Z.
the exactness of our approximation. We do not need to increase the order of R(θ) or divide Θ
any more.

18

It is interesting that an exact result is obtained with a low-degree R(θ) and a coarse division. This is consistent with the reported performance of other approaches based on positive
polynomials.

7. A robust SDP problem with function derivatives
In this section, we consider a robust SDP problem not only with a function but also with its
derivatives:
Pd : minimize

cT x

subject to E(x) º O,³
∂φ(θ) ∂φ(θ)∂φ(θ) ´F x, φ(θ),,,...,;θ º O∂θ1∂θ2∂θp(∀θ ∈ Θ).

Here, the optimization variables are a vector x ∈ Rn and a function φ ∈ Φd . The function set
Φd is a linear space of continuous and piecewise continuously differentiable functions definedon Θ. The function F is affine in x, φ(θ), (∂/∂θ1 )φ(θ), (∂/∂θ2 )φ(θ), . . ., (∂/∂θp )φ(θ) while
polynomial in θ. Other setting such as the size of E(x) is the same as in P .
Various control problems are formulated into Pd .
Suppose in Example 1 that the parameter θ varies with time in the set Θ and its time
derivative θ̇ belongs to the set Θ̇. Assume 0 ∈ Θ̇. The L2 -induced norm of the system can becomputed with an optimization problem similar to the one in Example 1 [27, 8, 28, 29]. The
difference is that the (1, 1)-block of the constraint is replaced by
p
X
∂R(θ)
i=1

∂θi

θ̇i − A(θ)R(θ) − R(θ)A(θ)T

with R(θ) being a continuous and piecewise continuously differentiable function. Moreover, theparameter of the problem is (θ, θ̇) ∈ Θ × Θ̇ here. This problem is in the form of Pd .
Another example is stability analysis of a nonlinear system.
˙ = f (ξ(t)) with a polynomial f (ξ) satisfying
Example 12. Consider a nonlinear system ξ(t)f (0) = 0. Then, the origin is asymptotically stable with this system if the optimum value ofthe following problem is negative:
minimizex

subject to φ(0) = 0,

1 − φ(ξ0 ) ≥ 0,
2

φ(ξ) + xkξk ≥ 0,2

xkξk −

19³ ∂φ(ξ) ´T∂ξf (ξ) ≥ 0(∀ξ ∈ Ξ).

Here, Ξ is some closed polytope containing the origin in its interior; ξ0 6= 0 is a point in Ξ. The
optimization variables are a scaler x and a Lyapunov function candidate φ, which is continuousand piecewise continuously differentiable. This problem is again in the form of Pd .
We consider an approximate approach to the problem Pd .
First, a general function φ(θ) is replaced by a polynomial φu (θ) =Pα∈S uα θα¤characterized

by u = (uα ) ∈ Rnu . The notation³∂φu (θ) ∂φu (θ)∂φu (θ) ´F (x, u; θ) := F x, φu (θ),
,
,...,
;θ
∂θ1
∂θ2
∂θp
is used for explicit representation of the finite dimensionality. Here, F is affine in x and u while
polynomial in θ.
Next, a division ∆ = {Θ[j] }Jj=1 of the parameter set Θ is prepared and a polynomial φu[j] (θ)
is considered for each subregion Θ[j] . Recall that the function φ(θ) has to be continuous over thewhole set Θ. This requirement is expressed by a linear equality constraint on u[1] , u[2] , . . . , u[J] ,which is denoted byd(u[1] , u[2] , . . . , u[J] ) = 0.
Finally, the matrices G(x, u), M (θ), H(θ) are defined in the same way as in P . Then, we
have an approximate problem:
Pd (∆) : minimize

cT x

subject to E(x) º O,

G(x, u[j] ) + H(θ)(W [j] )T + W [j] H(θ)T º O(∀θ ∈ ver Θ[j] , ∀j = 1, 2, . . . , J),d(u[1] , u[2] , . . . , u[J] ) = 0,
which can be solved with the standard interior-point method.
It is natural to expect an upper bound on the approximation error | inf Pd (∆) − inf Pd |.
For the present, such a bound is obtained in a rather special setting: the parameter set Θ is aQ
multi-dimensional interval pi=1 [θi , θi ], its division ∆ is a Cartesian division, and the polynomialφu (θ) is cubic in each θi . This is due to the requirement of continuity and the technical difficulty
to bound u[j] in the limit of rad ∆ → 0.
Q
A Cartesian division of Θ = pi=1 [θi , θi ] is defined as follows. For each i = 1, 2, . . . , p, divide
the interval [θi , θi ] as θi = θi0 < θi1 < · · · < θiQi −1 < θiQi = θi . Then, the collection of all theQmulti-dimensional intervals pi=1 [θiqi , θiqi +1 ] (qi = 0, 1, . . . , Qi − 1, i = 1, 2, . . . , p) is a Cartesiandivision of Θ.
When the polynomial φu (θ) is cubic in each θi (i = 1, 2, . . . , p), the dimension nu is equal to
4p nφ . Let us use this φu (θ) to interpolate a sufficiently smooth function φ(θ) on Θ. Pick up one
20subregion Θ[j] from a Cartesian division ∆. Then, we can choose φu[j] (θ) so as to interpolate thefunction value φ(θ) as well as the first derivatives (∂/∂θ1 )φ(θ), (∂/∂θ2 )φ(θ), . . ., (∂/∂θp )φ(θ)
at all the 2p vertices of Θ[j] . Since the required dimension is 2p (p + 1)nφ , the choice of u[j] ismade from the space of the dimension 4p nφ − 2p (p + 1)nφ . Whatever the particular choice is,if each u[j] is chosen as such the continuity condition d(u[1] , u[2] , . . . , u[J] ) = 0 is automaticallysatisfied. Moreover, it is possible to bound u[j] irrespectively of the maximum radius rad ∆ byexploiting the freedom in the choice of u[j] . More discussion can be found in Chapters 14–16 of[7].
We make the following two assumptions to obtain the bound. The first assumption is anatural generalization of Assumption 3 (a). Although the second assumption is a generalizationof Assumption 3 (b), it may look restrictive due to the bounds on the high-order derivatives of
φ(t). Note, however, that it is satisfied in a practically reasonable situation that the problem
Pd and its modification Pd,² , to be introduced, are optimized by some x and smooth enoughφ(θ).
Assumption 13.
(a) There exist x ∈ Rn and φ ∈ Φd such that E(x) Â O and F(x, φ(θ), (∂/∂θ1 )φ(θ), (∂/∂θ2 )φ(θ),. . . , (∂/∂θp )φ(θ); θ) Â O (∀θ ∈ Θ).
(b) There exist positive numbers ², x, and φα , where α = [α1 α2 · · · αp ]T and αi = 0, 1, 2, 3
for each i, having the following property: For any 0 ≤ ² ≤ ² and any v ∈ R, the set {(x, φ) ∈Rn × Φd | cT x ≤ v, E(x) º O, F(x, φ(θ), (∂/∂θ1 )φ(θ), (∂/∂θ2 )φ(θ), . . . , (∂/∂θp )φ(θ); θ) º
²I (∀θ ∈ Θ)} is either empty or having an element (x, φ) with kxk ≤ x and k(∂/∂θ1 )α1 (∂/∂θ2 )α2· · · (∂/∂θp )αp φ(θ)k ≤ φα for any θ ∈ Θ and αi = 0, 1, 2, 3 (i = 1, 2, . . . , p).

¤With this setting, we have the following bound on the approximation error of Pd (∆).
Theorem 14. Let Θ be a multi-dimensional interval and ∆ be its Cartesian division. Let also
φu (θ) be cubic in each θi . Then, under Assumption 13, there exist positive numbers Cd and Cd0 ,which are independent of the division ∆, such that| inf Pd (∆) − inf Pd | ≤ Cd rad ∆for any ∆ satisfying rad ∆ ≤ Cd0 .
The implication of this theorem is parallel to that of Theorem 4. This theorem immediately
implies the asymptotic exactness of Pd (∆), that is, the approximation error | inf Pd (∆) − inf Pd |converges to zero as the maximum radius rad ∆ of the division goes to zero. The convergence is
at least proportional to rad ∆. Exactness verification is possible in a similar fashion to Section 5.
21

Adaptive division as in Section 4 is not possible in the present setting because the division
is restricted to be Cartesian. Improvement on this point is a future research subject. In theliterature [4, 16, 9, 13, 23, 11, 6], a low-order polynomial is used for φu (θ) in place of the present
cubic function. An upper bound in this setting is another interesting subject.
The proof of Theorem 14 is similar to that of Theorem 4. We show only a sketch of theproof for avoidance of repetition.
Sketch of the proof of Theorem 14. We first consider the following auxiliary problem:
Pd,² : minimize

cT x

subject to E(x) º O,³
∂φ(θ) ∂φ(θ)∂φ(θ) ´,,...,; θ º ²I
F x, φ(θ),
∂θ1
∂θ2
∂θp

(∀θ ∈ Θ)for ² ≥ 0.
The critical step of the proof is to show inf Pd (∆) ≤ inf Pd,² for ∆ with small rad ∆ and forconsistently small ² > 0. To this aim, we construct a solution (x, {u[j] }, {W [j] }) of Pd (∆) froma solution (x, φ) of Pd,² .
Suppose Θ ⊆ [−1, 1]p first. If ² > 0 is small enough, the auxiliary problem Pd,² is feasible
due to Assumption 13 (a). Moreover, a feasible solution (x, φ) can be chosen so that kxk ≤ x
and k(∂/∂θ1 )α1 (∂/∂θ2 )α2 · · · (∂/∂θp )αp φ(θ)k ≤ φα , for any θ ∈ Θ and αi = 0, 1, 2, 3 owing
to Assumption 13 (b). Now, we choose a cubic function φu[j] (θ) for each j = 1, 2, . . . , J so
as to interpolate φ(θ) and its first derivatives at the vertices of Θ[j] . We also choose W [j] as(1/2 rad Θ[j] )H(θ0 ) for each j, where θ0 is the vertex of Θ[j] with the smallest coordinate ineach axis. Then, if the maximum radius of ∆ satisfies some linear relationship with ², this(x, {u[j] }, {W [j] }) forms a feasible solution of the approximate problem Pd (∆). This implies
inf Pd (∆) ≤ inf Pd,² .
The extension to the case of Θ 6⊆ [−1, 1]p is as in Section 4.2. Finally, the convexity of
inf Pd,² gives the desired bound.

¤

8. Conclusion
In this paper, we considered a robust SDP problem with a functional variable, which covers
many control problems with a parameter-dependent Lyapunov function. In particular, we
present an approximate approach to solve this problem and show its asymptotic exactness.
A quantitative error bound is given in terms of the resolution of the division. Verification of
exactness is also discussed. A numerical example shows the possibility of exact optimization
through a finite-dimensional approximate problem.
22

Acknowledgment
Part of this work is supported by a Grant-in-Aid of the Ministry of Education, Culture, Sports,
Science and Technology of Japan and the 2007 Pache Research Fund I-A-2 of Nanzan University,
Japan.

References[1] A. Ben-Tal and A. Nemirovski, Lectures on Modern Convex Optimization: Analysis, Algorithms,
and Engineering Applications. Philadelphia, USA: SIAM, 2001.
[2] P.-A. Bliman, “An existence result for polynomial solutions of parameter-dependent LMIs,” Systems & Control Letters, vol. 51, nos. 3–4, pp. 165–169, 2004.
[3] P.-A. Bliman, “On robust semidefinite programming,” in Proceedings of the 16th International
Symposium on Mathematical Theory of Networks and Systems (MTNS2004), Leuven, Belgium,
July 2004.
[4] R. Brayton and C. Tong, “Stability of dynamical systems: a constructive approach,” IEEE
Transactions on Circuits and Systems, vol. 26, no. 4, pp. 224–234, 1979.
[5] G. Chesi, A. Garulli, A. Tesi, and A. Vicino, “Polynomially parameter-dependent Lyapunov
functions for robust stability of polytopic systems: an LMI approach,” IEEE Transactions on
Automatic Control, vol. 50, no. 3, pp. 365–370, 2005.
[6] D. F. Coutinho and P. Danès, “Piecewise Lyapunov functions to the stability analysis of rational
systems subject to multiple state constraints,” in Proceedings of the 45th IEEE Conference on
Decision and Control, San Diego, USA, December 2006, pp. 5801–5806.
[7] G. Farin, Curves and Surfaces for CAGD: A Practical Guide, 5th ed. San Francisco, USA: Morgan
Kaufmann, 2002.
[8] P. Gahinet, P. Apkarian, and M. Chilali, “Affine parameter-dependent Lyapunov functions andreal parametric uncertainty,” IEEE Transactions on Automatic Control, vol. 41, no. 3, pp. 436–
442, 1996.
[9] M. Johansson and A. Rantzer, “On the computation of piecewise quadratic Lyapunov functions,”in Proceedings of the 36th IEEE Conference on Decision and Control, San Diego, USA, December
1997, pp. 3515–3520.
[10] J. B. Lasserre, “Global optimization with polynomials and the problem of moments,” SIAM
Journal on Optimization, vol. 11, no. 3, pp. 796–817, 2001.
[11] S. Lim and J. P. How, “Analysis of linear parameter-varying systems using a non-smooth dissipative systems framework,” International Journal of Robust and Nonlinear Control, vol. 12, no. 12,
pp. 1067–1092, 2002.
[12] J. Löfberg, “YALMIP: a toolbox for modeling and optimization in MATLAB,” inProceedings of the 2004 IEEE International Symposium on Computer Aided Control Systems Design, Taipei, Taiwan, September 2004, pp. 284–289, available fromhttp://control.ee.ethz.ch/∼ joloef/yalmip.php

23

[13] I. Masubuchi, A. Kume, and E. Shimemura, “Spline-type solution to parameter-dependent LMIs,”in Proceedings of the 37th IEEE Conference on Decision and Control, Tampa, USA, December
1998, pp. 1753–1758.
[14] A. T. Neto, “Parameter dependent Lyapunov functions for a class of uncertain linear systems: an
LMI approach,” in Proceedings of the 38th Conference on Decision and Control, Phoenix, USA,
December 1999, pp. 2341–2346.
[15] A. Ohara and Y. Sasaki, “On solvability and numerical solutions of parameter-dependent differential matrix inequality,” in Proceedings of the 40th IEEE Conference on Decision and Control,
Orlando, USA, December 2001, pp. 3593–3594.
[16] Y. Ohta, H. Imanishi, L. Gong, and H. Haneda, “Computer generated Lyapunov functions for aclass of nonlinear systems,” IEEE Transactions on Circuits and Systems I: Fundamental Theory
and Applications, vol. 40, no. 5, pp. 343–354, 1993.
[17] Y. Oishi, “A matrix-dilation approach to robust semidefinite programming and its error bound,”in Proceedings of the 2006 American Control Conference, Minneapolis, USA, June 2006, pp. 123–
129.
[18] Y. Oishi and Y. Isaka, “Exploiting sparsity in the matrix-dilation approach to robust semidefinite
programming,” in Proceedings of the 2007 American Control Conference, New York, USA, July
2007, pp. 6169–6176.
[19] R. C. L. F. Oliveira and P. L. D. Peres, “LMI conditions for robust stability analysis basedon polynomially parameter-dependent Lyapunov functions,” Systems & Control Letters, vol. 55,
no. 1, pp. 52–61, 2006.
[20] P. A. Parrilo, “Semidefinite programming relaxations for semialgebraic problems,” Mathematical
Programming, vol. 96, no. 2, pp. 293–320, 2003.
[21] D. Peaucelle, Y. Ebihara, D. Arzelier, and T. Hagiwara, “General polynomial parameterdependent Lyapunov functions for polytopic uncertain systems,” in Proceedings of the 17th International Symposium on Mathematical Theory of Networks and Systems (MTNS2006), Kyoto,
Japan, July 2006, pp. 2238–2242.
[22] M. Sato and D. Peaucelle, “Robust stability/performance analysis for linear time-invariant polynomially parameter-dependent systems using polynomially parameter-dependent Lyapunov functions,” in Proceedings of the 45th IEEE Conference on Decision and Control, San Diego, USA,
December 2006, pp. 5807–5813.
[23] C. W. Scherer, “Robust performance analysis for parameter dependent systems using tensor
product splines,” in Proceedings of the 37th IEEE Conference on Decision and Control, Tampa,
USA, December 1998, pp. 2216–2221.
[24] C. W. Scherer, “Relaxations for robust linear matrix inequality problems with verifications forexactness,” SIAM Journal on Matrix Analysis and Applications, vol. 27, no. 2, pp. 365–395, 2005.
[25] C. W. Scherer and C. W. J. Hol, “Matrix sum-of-squares relaxations for robust semi-definite
programs,” Mathematical Programming, vol. 107, nos. 1–2, pp. 189–211, 2006.
[26] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,”
Optimization Methods and Software, vols. 11–12, pp. 625–653, 1999.

24

[27] R. Watanabe, K. Uchida, M. Fujita, and E. Shimemura, “L2 gain and H∞ control of linear
systems with scheduling parameter,” in Proceedings of the 33rd IEEE Conference on Decision
and Control, Lake Buena Vista, USA, December 1994, pp. 1412–1414.
[28] F. Wu, X. H. Yang, A. Packard, and G. Becker, “Induced L2 -norm control for LPV systems withbounded parameter variation rates,” International Journal of Robust and Nonlinear Control,
vol. 6, nos. 9–10, pp. 983–998, 1996.
[29] J. Yu and A. Sideris, “H∞ control with parametric Lyapunov functions,” Systems & Control
Letters, vol. 30, nos. 2–3, pp. 57–69, 1997.

25