NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED
PROBLEMS : APPLICATIONS TO CONTROL DESIGN
Jean-Baptiste Thevenet
ONERA-CERT, 2 av. Edouard Belin, 31055 Toulouse, France
and UPS-MIP (Mathmatiques pour l’Industrie et la Physique), CNRS UMR 5640
118, route de Narbonne 31062 Toulouse, France
Dominikus Noll
UPS-MIP,
Pierre Apkarian
ONERA-CERT and UPS-MIP,
Keywords:
Bilinear matrix inequality, spectral penalty function, trust-region, control synthesis.
Abstract:
The purpose of this paper is to examine a nonlinear spectral semidefinite programming method to solve prob-
lems with bilinear matrix inequality (BMI) constraints. Such optimization programs arise frequently in au-
tomatic control and are difficult to solve due to the inherent non-convexity. The method we discuss here is
of augmented Lagrangian type and uses a succession of unconstrained subproblems to approximate the BMI
optimization program. These tangent programs are solved by a trust region strategy. The method is tested
against several difficult examples in feedback control synthesis.
1 INTRODUCTION
Minimizing a linear objective function subject to bi-
linear matrix inequality (BMI) constraints is a useful
way to describe many robust control synthesis prob-
lems. Many other problems in automatic control lead
to BMI feasibility and optimization programs, such as
filtering problems, synthesis of structured or reduced-
order feedback controllers, simultaneous stabilization
problems and many others. Formally, these problems
may be described as
minimize c
T
x, x R
n
subject to B(x) 0,
(1)
where 0 means negative semi-definite and
B(x)=A
0
+
n
i=1
x
i
A
i
+
1i<jn
x
i
x
j
B
ij
(2)
is a bilinear matrix-valued operator with data
A
i
,B
ij
S
m
.
The importance of BMIs was recognized over the
past decade and different solution strategies have been
proposed. The earliest ideas use interior point meth-
ods (imported from semidefinite programming), con-
cave programming or coordinate descent methods
employing successions of SDP subproblems. The
success of these methods is up to now very limited,
and even moderately sized programs with more than
100 variables usually make such methods fail.
Significant progress was achieved by two recent
approaches, both using successions of SDP tangent
subproblems. In (Fares et al., 2002), a sequen-
tial semidefinite is used, programming (S-SDP) al-
gorithm, expanding on the sequential quadratic pro-
gramming (SQP) technique in traditional nonlinear
programming. In (Noll et al., 2002; Apkarian et al.,
2003b; Fares et al., 2001) an augmented Lagrangian
strategy is proposed, which was since then success-
fully tested against a considerable number of difficult
synthesis problems (Fares et al., 2000; Apkarian et al.,
2002; Apkarian et al., 2003a; Noll et al., 2002). This
approach is particularly suited when BMI problems
are transformed to LMI-constrained programs with an
additional nonlinear equality constraints:
minimize c
T
x, x R
n
subject to A(x) 0
g(x)=0.
(3)
Here A is an affine operator (i.e., (2) with B
ij
=0)
with values in S
m
, and g : R
n
R
p
incorporates a
finite number of equality constraints. Notice that ev-
ery BMI problem may in principle be transformed to
the form (3), but the change will only be fruitful if it
is motivated from the control point of view. The lat-
ter is for instance the case in robust or reduced order
synthesis, where the projection lemma is successfully
brought into play. For semidefinite programming,
(Zibulevsky, 1996; Ben-Tal and Zibulevsky, 1997)
propose a modified augmented Lagrangian method,
237
Thevenet J., Noll D. and Apkarian P. (2004).
NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS : APPLICATIONS TO CONTROL DESIGN.
In Proceedings of the First International Conference on Informatics in Control, Automation and Robotics, pages 237-248
DOI: 10.5220/0001128402370248
Copyright
c
SciTePress
which in (Kocvara and Stingl, 2003) has been used to
build a software tool for SDP called PENNON. The
idea is further expanded in (Henrion et al., 2003) to
include BMI-problems.
The present paper addresses cases where the pas-
sage from (1) to (3) seems less suited and a direct
approach is required. Our algorithm solves (1) by a
sequence of unconstrained subproblems via a trust re-
gion technique. Using trust regions is of the essence,
since we have to bear in mind that (1) is usually highly
non-convex. In consequence, local optimal methods
risk failure by getting stalled at a local minimum of
constraint violation. This refers to points x satisfy-
ing B(x)  0, where the algorithm makes no further
progress. Such a case means failure to solve the un-
derlying control problem, which requires at least x
with B(x) 0. This phenomenon may be avoided
(or at least significantly reduced) with genuine second
order strategies. In particular, trust region techniques
do not require convexification of the tangent problem
Hessian, a major advantage over line search methods.
(Notice that in this situation, any approach to (1) us-
ing a succession of SDPs is bound to face similar dif-
ficulties due to the need to convexify the Hessian. On
the other hand, if LMI-constrained subproblems with
non-convex quadratic objectives f (x) can be handled,
such an approach is promising. This idea has been de-
veloped successfully in (Apkarian et al., 2003a).
While our technique could be applied basically to
any BMI-constrained program, we presently focus on
automatic control issues. In particular, we give ap-
plications in static output feedback control under H
2
,
H
or mixed H
/H
2
performance indices. Two
other difficult cases, multi-model and multi-objective
synthesis, and synthesis of structured controllers, are
presented.
The structure of the paper is as follows. In sec-
tion 2, we present our spectral SDP method and dis-
cuss its principal features. Section 3 presents the-
ory needed to compute first and second derivatives of
spectral functions. Finally, in section 4, numerical ex-
periments on automatic control issues are examined.
NOTATION
Our notation is standard. We let S
m
denote the set
of m × m symmetric matrices, M
T
the transpose
of the matrix M and Tr M its trace. We equip S
m
with the euclidean scalar product X, Y = X Y =
Tr (XY ). For symmetric matrices , M N means
that M N is positive definite and M N means
that M N is positive semi-definite. The symbol
stands for the usual Kronecker product of matrices
and vec stands for the columnwise vectorization on
matrices. We shall make use of the properties:
vec (AXB)=(B
T
A)vecX, (4)
Tr (AB)=vec
T
A
T
vec B, (5)
which hold for matrices A, X and B of compatible di-
mensions. The Hadamard or Schur product is defined
as
A B =((A
ij
B
ij
)). (6)
The following holds for matrices of the same dimen-
sion:
vec A vec B = diag(vec A)vec B, (7)
where the operator diag forms a diagonal matrix with
vec A on the main diagonal.
2 NONLINEAR SPECTRAL SDP
METHOD
In this chapter we present and discuss our method to
solve BMI-constrained optimization problems. Sec-
tion 2.1 gives the general out-set, while the subse-
quent sections discuss theoretical and implementa-
tional aspects of the method.
2.1 General outline
The method we are about to present applies to more
general classes of objective functions than (1). We
shall consider matrix inequality constrained programs
of the form
minimize f(x),x R
n
subject to F(x) 0,
(8)
where f is a class C
2
function, F : R
n
S
m
a class
C
2
operator. We will use the symbol B whenever we
wish to specialize to bilinear operators (2) or to the
form (1), and we use A for affine operators.
Our method is a penalty/barrier multiplier algo-
rithm as proposed in (Zibulevsky, 1996; Mosheyev
and Zibulevsky, 2000). According to that terminol-
ogy, a spectral penalty (SP) function for (1) is defined
as
F (x, p)=f(x)+TrΦ
p
(F(x)), (9)
where ϕ
p
: R R is a parametrized family of scalar
functions, which generates a family of matrix-valued
operators Φ
p
: S
m
S
m
upon defining:
Φ
p
(X):=Sdiag [ϕ
p
(λ
i
(X))] S
T
. (10)
Here λ
i
(X) stands for the ith eigenvalue of X
S
m
and S is an orthonormal matrix of associated
ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION
238
eigenvectors. An alternative expression for (9) using
(10) is:
F (x, p)=f(x)+
m
i=1
ϕ
p
(λ
i
(F(x)) . (11)
As ϕ
p
is chosen strictly increasing and satisfies
ϕ
p
(0) = 0, each of the following programs
(parametrized through p>0)
minimize f(x),x R
n
subject to Φ
p
(F(x)) 0
(12)
is equivalent to (8). Thus, F (x, p) may be under-
stood as a penalty function for (8). Forcing p 0,
we expect the solutions to the unconstrained program
min
x
F (x, p) to converge to a solution of (8).
It is well-known that pure penalty methods run into
numerical difficulties as soon as penalty constants get
large. Similarly, using pure SP functions as in (9)
would lead to ill-conditioning for small p>0. The
epoch-making idea of Hestenes (Hestenes, 1969) and
Powell (Powell, 1969), known as the augmented La-
grangian approach, was to avoid this phenomenon by
including a linear term carrying a Lagrange multiplier
estimate into the objective. In the present context, we
follow the same line, but incorporate Lagrange multi-
plier information by a nonlinear term. We define the
augmented Lagrangian function associated with the
matrix inequality constraints in (1) as
L(x, V, p)=f(x)+TrΦ
p
(V
T
F(x)V ) , (13)
= f (x)+
m
i=1
ϕ
p
(λ
i
(V
T
F(x)V )).
In this expression, the matrix variable V has the
same dimension as F(x) S
m
and serves as a fac-
tor of the Lagrange multiplier variable U S
m
, U =
VV
T
(for instance, a Cholesky factor). This has the
immediate advantage of maintaining non-negativity
of U 0. In contrast with classical augmented La-
grangians, however, the Lagrange multiplier U is not
involved linearly in (13). We nevertheless reserve
the name of an augmented Lagrangian for L(x, V, p),
as its properties resemble those of the classical aug-
mented Lagrangian. Surprisingly, the non-linear de-
pendence of L(x, V, p) on V is not at all troublesome
and a suitable first-order update formula V V
+
,
generalizing the classical one, will be readily derived
in section 3.2. We will even see some genuine advan-
tages of (13) over the case of a linear multiplier U .
Let us mention that the convergence theory for an
augmented Lagrangian method like the present one
splits into a local and a global branch. Global the-
ory gives weak convergence of the method towards
critical points from arbitrary starting points x, if nec-
essary by driving p 0. An important complement
is provided by local theory, which shows that as soon
as the iterates reach a neighbourhood of attraction of
one of those critical points predicted by global theory,
and if this critical point happens to be a local mini-
mum satisfying the sufficient second order optimality
conditions, then the sequence will stay in this neigh-
bourhood, converge to the minimum in question, and
the user will not have to push the parameter p below
a certain threshold ¯p>0. Proofs for global and lo-
cal convergence of the AL exist in traditional nonlin-
ear programming (see for instance (Conn et al., 1991;
Conn et al., 1996; Conn et al., 1993b; Conn et al.,
1993a; Bertsekas, 1982)). Convergence theory for
matrix inequality constraints is still a somewhat unex-
plored field. A global convergence result which could
be adapted to the present context is given in (Noll
et al., 2002). Local theory for the present approach
covering a large class of penalty functions ϕ
p
will be
presented in a forthcoming article. Notice that in the
convex case a convergence result has been published
in (Zibulevsky, 1996). It is based on Rockafellar’s
idea relating the AL method to a proximal point algo-
rithm. Our present testing of nonconvex programs (1)
shows a similar picture. Even without convexity the
method converges most of the time and the penalty
parameter is in the end stably away from 0, p ¯p.
Schematically, the augmented Lagrangian tech-
nique is as follows:
Spectral augmented Lagrangian algorithm
1. Initial phase. Set constants γ>0 < 1.
Initialize the algorithm with x
0
, V
0
and a penalty
parameter p
0
> 0.
2. Optimization phase. For fixed V
j
and p
j
solve
the unconstrained subproblem
minimize
xR
n
L(x, V
j
,p
j
) (14)
Let x
j+1
be the solution. Use the previous iterate
x
j
as a starting value for the inner optimization.
3. Update penalty and multiplier. Apply first-
order rule to estimate Lagrange multiplier :
V
j+1
V
T
j+1
= V
j
S
diagϕ
p
λ
i
V
T
j
·F(x
j+1
)V
j
))] S
T
V
T
j
, (15)
where S diagonalizes V
T
j
F(x
j+1
)V
j
.
Update the penalty parameter using :
p
j+1
=
ρp
j
, if λ
max
(0, F(x
j+1
))
λ
max
(0, F(x
j
))
p
j
, else
(16)
Increase j and go back to step 2.
NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS: APPLICATIONS TO CONTROL
DESIGN
239
In our implementation, following the recommenda-
tion in (Mosheyev and Zibulevsky, 2000), we have
used the log-quadratic penalty function ϕ
p
(t)=
1
(t/p) where
ϕ
1
(t)=
t +
1
2
t
2
if t ≥−
1
2
1
4
log(2t)
3
8
if t<
1
2
,
(17)
but other choices could be used (see for instance
(Zibulevsky, 1996) for an extended list).
The multiplier update formula (15) requires the full
machinery of differentiability of the spectral function
Tr Φ
p
and will be derived in section 3.2. Algorithmic
aspects of the subproblem in step 2 will be discussed
in section 3.3. We start by analyzing the idea of the
spectral AL-method.
2.2 The mechanism of the algorithm
In order to understand the rationale behind the AL-
algorithm, it may instructive to consider classical non-
linear programming, which is a special case of the
scheme if the values of the operator F are diagonal
matrices: F(x) = diag [c
1
(x),...,c
m
(x)]. Then
the Lagrange multiplier U and its factor V may also
be restricted to the space of diagonal matrices, U =
diag u
i
, V = diag v
i
, and we recover the situation of
classical polyhedral constraints. Switching to a more
convenient notation, the problem becomes
minimize f(x)
subject to c
j
(x) 0,j=1,...,m
With u
i
= v
2
i
, we obtain the analogue of the aug-
mented Lagrangian (13):
L(x, v, p)=f(x)+
m
i=1
1
(v
2
i
c
i
(x)/p).
Here we use (17) or another choice from the list in
(Zibulevsky, 1996). Computing derivatives is easy
here and we obtain
L(x, v, p)=f(x)
+
m
i=1
ϕ
1
(v
2
i
c
i
(x)/p)v
2
i
c
i
(x).
If we compare with the Lagrangian :
L(x, u)=f(x)+u
T
c(x),
and its gradient :
L(x, u)=f(x)+
m
i=1
u
i
c
i
(x),
the following update formula readily appears:
u
+
i
= v
+2
i
= v
2
i
ϕ
1
(v
2
i
c
i
(x
+
)/p),
the scalar analogue of (15). Suppose now that ϕ
1
has
the form (17). Then for sufficiently small v
2
i
c
i
(x) we
expect v
2
i
c
i
(x)/p >
1
2
,so
1
(v
2
i
c
i
(x)/p)=p
v
2
i
c
i
(x)/p
+
1
2
v
4
i
c
i
(x)
2
/p
2
= u
i
c
i
(x)+
1
2
(v
4
i
/p)c
i
(x)
2
.
If we remember the form of the classical augmented
Lagrangian for inequality constraints (cf. (Bertsekas,
1982))
L(x, u, µ)=f(x)
+
µ
2
m
i=1
max [0,u
i
+ c
i
(x)]
2
u
2
i
,
we can see that the term v
4
i
/p takes the place of
the penalty parameter 1as long as the v
i
s remain
bounded. This suggests that the method should be-
have similarly to the classical augmented Lagrangian
method in a neighbourhood of attraction of a local
minimum. In consequence, close to a minimum the
updating strategy for p should resemble that for the
classical parameter µ. In contrast, the algorithm may
perform very differently when iterates are far away
from local minima. Here the smoothness of the func-
tions ϕ
1
may play favorably. We propose the update
(16) for the penalty parameter, but other possibilities
exist, see for instance Conn et al., and need to be
tested and compared in the case of matrix inequality
constraints.
3 DERIVATIVES OF SP
FUNCTIONS
In this section we obtain formulas for the first-
and second-order derivatives of the augmented La-
grangian function (13). This part is based on the dif-
ferentiability theory for spectral functions developed
by Lewis et al. (Lewis, 1996; Lewis, 2001; Lewis and
S.Sendov, 2002). To begin with, recall the following
definition from (Lewis and S.Sendov, 2002).
Definition 3.1 Let λ : S
m
R
m
denote the eigen-
value map λ(X):=(λ
1
(X),...,λ
m
(X)). For a
symmetric function ψ : R
m
R, the spectral func-
tion Ψ associated with ψ is defined as Ψ:S
m
R,
Ψ=ψ λ.
First order differentiability theory for spectral func-
tions is covered by (Lewis, 1996). We will need the
following result from (Lewis, 1996):
Lemma 3.2 Let ψ a symmetric function defined on
an open subset D of R
n
. Let Ψ=ψλ be the spectral
function associated with ψ. Let X S
m
and suppose
ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION
240
λ(X) D. Then Ψ is differentiable at X if and only
if ψ is differentiable at λ(X). In that case we have
the formula:
Ψ(X)=(ψ λ)(X)=S diag ψ (λ(X)) S
T
for every orthogonal matrix S satisfying
X = S diagλ(X) S
T
.
The gradient Ψ(X) is a (dual) element of S
m
,
which we distinguish from the differential, dΨ(X).
The latter acts as a linear form on tangent vectors
dX S
m
as
dΨ(X)[dX]=Ψ(X) dX
=Tr
S diag ψ (λ(X))S
T
dX
.
It is now straightforward to compute the first deriva-
tive of a composite spectral functions Ψ ◦F, where
F : R
n
S
m
is sufficiently smooth. We obtain
◦F)(x)=dF(x)
[Ψ(F(x))] R
n
, (18)
where dF(x) is a linear operator mapping R
n
S
m
,
and dF(x)
its adjoint, mapping S
m
R
n
. With the
standing notation
F := F(x),F
i
:=
F(x)
∂x
i
S
m
,
we obtain the representations
dF(x)[δx]=
n
i=1
F
i
δx
i
,
dF(x)
[dX]=(F
1
dX,...,F
n
dX) .
Then the ith element of the gradient is
( ◦F)(x))
i
=Tr
F
i
Sdiagψ (λ(F )) S
T
. (19)
Observe that F
i
= A
i
+
j<i
B
ji
x
j
+
i<j
B
ij
x
j
if B is of the form (2), F
i
= A
i
for affine operators
A.
We are now ready for the first derivative of the aug-
mented Lagrangian function L(x, V, p). We consider
a special class of symmetric functions ψ : R
n
R,
ψ(x)=
n
i=1
ϕ(x
i
), (20)
generated by a single scalar function ϕ : R R.
Here the spectral function Ψ=ψ λ is of the form
Ψ=TrΦwith Φ:S
m
S
m
the matrix-valued
operator
Φ(X)=S diag ϕ (λ(X)) S
T
,
S being an orthogonal matrix diagonalizing X. Using
ϕ
p
, ψ
p
and Ψ
p
instead of ϕ, ψ and Ψ, the penalty term
in (13) takes the form
Ψ
p
(V
T
F(x)V )) = ψ
p
λ(V
T
F(x)V ))
=
m
i=1
ϕ
p
(λ
i
(V
T
F(x)V ))
=TrΦ
p
(V
T
F(x)V ),
where Φ
p
is given in (10). Now we apply Lemma
3.2, respectively formula (18) or even (19) to Ψ
p
=
Tr Φ
p
and the operator
F(x)=V
T
F(x)V , where
we observe en passant that
d
F(x)
[X]=d

V
T
FV
(x)
[X]
= dF(x)
[VXV
T
],
or put differently,
F
i
= V
T
F
i
V . This gives the fol-
lowing
Proposition 3.3 With the above notations the gradi-
ent of L(x, V, p) is given as
∂x
i
L(x, V, p)=
∂x
i
f(x)
+Tr
F
i
VSdiag [ϕ
p
(λ
i
(V
T
FV))] S
T
V
T
.
In vector notation
L(x, V, p)=f(x)+[vec (F
1
),...,vec (F
n
)]
T
· vec
VSdiag [ϕ
p
(λ
i
(V
T
FV))] S
T
V
T
,
where S is orthogonal and diagonalizes V
T
FV.
3.1 Second derivatives
The second order theory of spectral function is more
involved and presented in (Lewis, 2001; Lewis and
S.Sendov, 2002). We require the following central ob-
servation.
Lemma 3.4 With the same notations as above, Ψ=
ψ λ is twice differentiable at X S
m
if and only if
ψ is twice differentiable at λ(X) D.
We specialize again to the class of spectral func-
tions Ψ=TrΦgenerated by a single scalar function
ϕ, see (20). Then we have the following result im-
ported from (Shapiro, 2002):
Lemma 3.5 Let ϕ : R R be a function on the real
line, and let ψ be the symmetric function on R
n
gener-
ated by (20) above. Let Ψ=TrΦdenote the spectral
function associated with ψ. Suppose ϕ is twice differ-
entiable on an open interval I. Then Ψ is twice differ-
entiable at X whenever λ
i
= λ
i
(X) I for every i.
Moreover, if we define the matrix ϕ
[1]
= ϕ
[1]
(X)
S
m
as:
ϕ
[1]
ij
:=
ϕ
(λ
i
)ϕ
(λ
j
)
λ
i
λ
j
if λ
i
= λ
j
ϕ

(λ
i
) if λ
i
= λ
j
,
then the following formula defines the second order
derivative of Ψ at X
d
2
Ψ(X))[dX, dX]=Tr(S
T
dXSϕ
[1]
◦{S
T
dXS}),
where S is orthogonal and diagonalizes X.
NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS: APPLICATIONS TO CONTROL
DESIGN
241
In order to obtain second derivatives of composite
functions, Ψ◦F, we will require the chain rule, which
for the second differential forms d
2
Ψ(X)[dX, dX] of
Ψ and d
2
F(x)[δx,δx] of F amounts to:
d
2
◦F)(x)[δx, δx]
= dΨ(F(x))
d
2
F(x)[δx,δx]
+ d
2
Ψ(F(x)) [dF(x)[δx],dF(x)[δx]] . (21)
For bilinear B, the second order form d
2
B(x) is of
course independent of x and given as
d
2
B(x)[δx, δx]=
i<j
δx
i
δx
j
B
ij
.
In particular, d
2
A =0for affine operators A. During
the following we shall use these results to obtain sec-
ond derivatives of the augmented Lagrangian function
L(x, V, p).
As in the previous section, we apply the results to
Ψ
p
=TrΦ
p
generated by ϕ
p
and to the operator
F(x)=V
T
F(x)V . As we have already seen, for
F we have the formula
d
F(x)[δx]=
n
i=1
V
T
F
i
x
i
= V
T
n
i=1
F
i
δx
i
V
= V
T
dFV.
Therefore, using Lemma 3.5, the second term in (21)
becomes
d
2
Ψ
p
(V
T
FV)[V
T
dF V,V
T
dF V ]
=Tr
S
T
V
T
dF VSϕ
[1]
p
◦{S
T
V
T
dF VS}
,
where S is orthogonal and diagonalizes V
T
FV and
ϕ
[1]
p
= ϕ
[1]
p
(V
T
FV). This may now be transformed
to vector-matrix form using the operator vec and the
relations (4) - (7).
d
2
Ψ
p
(V
T
FV)[V
T
dF V,V
T
dF V ]
=vec
T
((VS)
T
dF (VS))
· vec (ϕ
[1]
p
◦{(VS)
T
dF (VS)})
=vec
T
(dF) K
p
vec (dF) , (22)
where the matrix K
p
is
K
p
=[(VS) (VS)]diag
vec ϕ
[1]
p
· [(VS)
T
(VS)
T
].
Using dF =
n
i=1
F
i
δx
i
, the last line in (22) gives
d
2
Ψ
p
(V
T
FV)[V
T
dF V,V
T
dF V ]
= δx
T
[vec (F
1
),...,vec (F
n
)]
T
K
p
· [vec (F
1
),...,vec (F
n
)] δx,
so we are led to retain the symmetric n × n matrix
H
1
=[vec(F
1
),...,vec (F
n
)]
T
K
p
· [vec (F
1
),...,vec (F
n
)] . (23)
Let us now consider the second part of the
Hessian
2
L(x, V, p), coming from the term
dΨ(F(x))
d
2
F(x)[δx,δx]
in (21). First observe
that trivially
d
2
F(x)[δx,δx]=d
2
{V
T
FV }(x)[δx, δx]
= V
T
d
2
F(x)[δx,δx] V.
Hence
dΨ
p
(V
T
FV)
V
T
d
2
F(x)[δx,δx]V
=Tr
S diag ϕ
p
λ
i
(V
T
FV)
S
T
·V
T
d
2
F(x)[δx,δx]V
=Tr
VSdiag ϕ
p
λ
i
(V
T
FV)
S
T
V
T
·d
2
F(x)[δx,δx]
. (24)
Introducing matrices F
ij
= F
ij
(x) S
m
to repre-
sent
d
2
F(x)[δx,δx]=
n
i,j=1
F
ij
δx
i
δx
j
,
relation (24) may be written as
dΨ
p
(V
T
FV)
V
T
d
2
F(x)[δx,δx]V
=
n
i,j=1
δx
i
δx
j
Tr
F
ij
VSdiag ϕ
p
λ
i
(V
T
·FV)) S
T
V
T
.
Hence the interest to define a second symmetric n ×n
matrix H
2
by
(H
2
)
ij
=Tr
F
ij
VSdiagϕ
p
λ
i
(V
T
FV)
S
T
V
T
,
(25)
where as before S is orthogonal and diagonalizes
V
T
FV. We summarize by the following
Proposition 3.6 The Hessian of the augmented La-
grangian function L(x, V, p) is of the form
2
L(x, V, p)=
2
f(x)+H
1
+ H
2
,
where H
1
, given by (23), is positive semidefinite. H
2
is given by (25).
Proof. The only element which remains to be proved
is non-negativity of H
1
, which hinges on the non-
negativity of K
p
, and hence on that of ϕ
[1]
p
. The latter
is a consequence of the convexity of ϕ
p
, hence the re-
sult.
ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION
242
3.2 Multiplier update rule
The first-order multiplier update rule is derived sim-
ilarly to the classical case of polyhedral constraints.
For our analysis we will require the traditional La-
grangian of problem (8), which is
L(x, U)=f (x)+U •F(x). (26)
Suppose x
+
is the solution of problem (14) in step
2 of our algorithm. Let U = VV
T
the current La-
grange multiplier estimate. A new multiplier U
+
and therefore a new factor estimate V
+
with U
+
=
V
+
V
+T
is then obtained by equating the gradients
of the augmented Lagrangian in (13) and the tradi-
tional Lagrangian of problem (1). Writing L(x, U )=
f(x)+Tr(U F(x)) = f(x)+Tr(VV
T
F(x)) =
f(x)+Tr(V
T
F(x)V ), and computing its gradient
at x
+
and U
+
= V
+
V
+T
implies
L(x
+
,U
+
)=f(x
+
)
+[vec(F
1
),...,vec (F
n
)]
T
vec (V
+
V
+
T
).
Comparing with
L(x
+
,V
+
,p)
= f(x
+
) + [vec (F
1
),...,vec (F
n
)]
T
· vec
VSdiag ϕ
p
(λ
i
(V
T
F(x
+
)V )S
T
V
T
suggests the update rule
V
+
V
+
T
= VS[diag ϕ
p
(λ
i
(V
T
F(x
+
)V ))] (VS)
T
,
(27)
where S is orthogonal and diagonalizes V
T
F(x
+
)V .
This was already presented in our algorithm. The ar-
gument suggests the equality to hold only modulo the
null space of [vec (F
1
),...,vec (F
n
)]
T
, but the limit-
ing case gives more information. Notice that equality
Lx,
¯
U)=Lx,
¯
V,p)
holds at a Karush-Kuhn-Tucker pair x,
¯
U) and
¯
U =
¯
V
¯
V
T
for every p>0, a formula whose analogue
in the classical case is well-known. That means we
would (27) expect to hold at x
+
x, V = V
+
=
¯
V . This is indeed the case since by complementarity,
¯
V
T
Fx)
¯
V =0. With ϕ
p
(0) = 1 for every p>
0, the right hand side in (27) is
¯
V
T
S
T
S
¯
V =
¯
V
T
¯
V ,
hence equals the left hand side. We may therefore
read (27) as some sort of fixed-point relation, which
is satisfied in the limit V = V
+
if iterates converge.
Notice, however, that convergence has to be proved by
an extra argument, because (27) is not of contraction
type.
Notice that by construction the right hand term in
(27) is indeed positive definite, so V
+
could for in-
stance be chosen as a Cholesky factor.
3.3 Solving the subproblem -
implementational issues
Efficient minimization of L(x, V
j
,p
j
) for fixed V
j
, p
j
is at the core of our approach and our implementation
is based on a Newton trust region method, following
the lines of Lin and Mor
´
e (Lin and More, 1998). As
compared to Newton line search algorithms or other
descent direction methods, trust regions can take bet-
ter advantage of second-order information. This is
witnessed by the fact that negative curvature in the
tangent problem Hessian, frequently arising in BMI-
minimization when iterates are still far away from any
neighbourhood of local convergence, may be taken
into account. Furthermore, trust region methods of-
ten (miraculously) find good local minima, leaving
the bad ones behind. This additional benefit is in
contrast with what line search methods achieve, and
is explained to some extent by the fact that, at least
over the horizon specified by the current trust region
radius, the minimization in the tangent problem is a
global one.
As part of our testing, and at an early stage of the
development, we studied the behavior of our code on
a class of LMI problems, where we compared line
search and trust region approaches. Even in that situa-
tion, where line search is supposedly at its best due to
convexity, the trust region variant appeared to be more
robust, though often slower. This is significant when
problems get ill-conditioned, which is often the case
in automatic control applications. Then the Newton
direction, giving descent in theory, often fails to do so
due to ill-conditioning. In that event recourse to first-
order methods has to be taken. This is where a trust
region strategy continues to perform reasonably well.
As currently implemented, both the trust region and
line search variant of our method apply a Cholesky
factorization to the tangent problem Hessian. This
serves either to build an efficient preconditioner, or
to compute the Newton direction in the second case.
When the Hessian is close to becoming indefinite, a
modified Cholesky factorization is used (
2
L + E =
JJ
T
, with J a lower triangular matrix and E a shift
matrix of minimal norm). This is a critical point of the
algorithm, particularly for our hard problems, because
the condition number of the shifted Hessian has to
be reasonable to avoid large computational roundoff
errors. However, trading a larger condition number
may be often desirable, for example when the Hes-
sian at the solution is ill-conditioned. We found that
the revised modified Cholesky algorithm proposed by
(Schnabel and Eskow, 1999) achieves this goal fairly
well.
There are a few more features of our approach,
which in our opinion merit special attention as they
seem to be typical for automatic control applications.
NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS: APPLICATIONS TO CONTROL
DESIGN
243
A first issue is the magnitude of the decision vari-
ables, which may vary greatly within a given prob-
lem. In particular the range of the Lyapunov variables
may show differences of several orders of magnitude,
a phenomenon which is difficult to predict in advance,
as the physical meaning of the Lyapunov variables
is often lacking. Properly scaling these variables is
therefore difficult in general.
A similar issue is to give prior bounds on the con-
troller gains. This is mandatory because these values
are to be useful in the hard-ware implementations of
the theoretically computed controllers. Exceedingly
large gains would bear the risk to amplify measure-
ment noise and thwart the process. This may also be
understood as a quest on the well-posedness of the
initial control problem, which should include explicit
gain bound constraints.
Incomplete knowledge about the numerical range
of the different variables is in our opinion the source
of most numerical difficulties and a method to balance
data is extremely important.
4 NUMERICAL EXAMPLES
In this section we test our code on a variety of difficult
applications in automatic control. The first part uses
examples from static output feedback control. We
then discuss applications like sparse linear constant
output-feedback design, simultaneous state-feedback
stabilization with limits on feedback gains and finally
multi-objective H
2
/H
-controller design. Here we
import examples from (Hassibi et al., 1999). In all
these cases our method achieves closed-loop stability,
while performance indices like H
, H
2
or mixed per-
formances previously published in the literature are
significantly improved. We always start our algorithm
at x
0
=0in order to capture a realistic scenario,
where no such prior information is available.
4.1 Static output feedback controller
synthesis
We first recall a BMI characterization of the classical
output feedback synthesis problem. To this end let
P (s) be an LTI (Linear Time Invariant) system with
state-space equations:
P (s):
˙x
z
y
=
AB
1
B
2
C
1
D
11
D
12
C
2
D
21
D
22

x
w
u
, (28)
where
x R
n
1
is the state vector,
u R
m
2
is the vector of control inputs,
w R
m
1
is a vector of exogenous inputs,
y R
p
2
is the vector of measurements,
z R
p
1
is the controlled or performance vector.
D
22
=0is assumed without loss of generality.
Let T
w,z
(s) denote the closed-loop transfer functions
from w to z for some static output-feedback control
law:
u = Ky . (29)
w
u
K(s)
P(s)
z
y
Figure 1: Output-feedback Linear Fractional Transform.
Our aim is to compute K subject to the following con-
straints:
internal stability: for w =0the state vector of the
closed-loop system (28) and (29) tends to zero as
time goes to infinity.
performance: the H
norm T
w,z
(s)
, respec-
tively the H
2
norm T
w,z
(s)
2
, is minimized
where the closed-loop transfer T
w,z
(s) is described
as :
˙x =(A + B
2
KC
2
)x +(B
1
+ B
2
KD
21
)w
z =(C
1
+ D
12
KC
2
)x +(D
11
+ D
12
KD
21
)w.
Recall that in the H
2
case, some feedthrough terms
must be nonexistent in order for the H
2
performance
to be well defined. We have to assume D
11
=
0,D
21
=0for the plant in (28). This guaran-
teed, both static H
and H
2
indices are then trans-
formed into a matrix inequality condition using the
bounded real lemma (Anderson and Vongpanitlerd,
1973). This leads to well-know characterizations:
Proposition 4.1 A stabilizing static output feedback
controller K with H
gain T
w,z
(s)
γ exists
provided there exists X S
n
1
such that:
(A+B
2
KC
2
)
T
X+ ∗∗
(B
1
+B
2
KD
21
)
T
X γI
(C
1
+D
12
KC
2
)(D
11
+D
12
KD
21
) γI
0
(30)
X
0. (31)
Proposition 4.2 A stabilizing static output feedback
controller K with H
2
performance T
w,z
(s)
2
γ
exists provided there exist X, M S
n
1
such that:
ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION
244
(A+B
2
KC
2
)
T
X+ ∗∗
(C
1
+D
12
KC
2
) I
0 (32)
M
X (B
1
+B
2
KD
21
) X
0 (33)
Tr (M )
γ. (34)
In order to convert these programs to the form (1), we
have to replace strict inequalities 0 by − for
a suitable threshold. The choice of >0 poses no
problem in practice.
4.1.1 Transport airplane (TA)
We consider the design of a static H
controller for
a transport airplane from (Gangsaas et al., 1986). Our
algorithm computes a static controller K after solving
29 unconstrained minimization subproblems (14):
K
=
0.69788
0.64050
0.83794
0.09769
1.57062
T
.
The associated H
performance is 2.22 and repre-
sents 30% improvement over the result in (Leibfritz,
1998).
4.1.2 VTOL helicopter (VH)
The state-space data of the VTOL helicopter are bor-
rowed from (Keel et al., 1988). They are obtained
by linearizing the helicopter rigid body dynamics at
given flight conditions. This leads to a fourth-order
model. Both H
and H
2
static controllers are com-
puted with the spectral SDP method. The algorithm
converges after solving 29 tangent problems (14):
K
=
0.50750
10.00000
.
The final H
performance is 1.59e-1, which is about
40% improvement over the performance obtained in
(Leibfritz, 1998).
K
2
=
0.13105
5.95163
,
which gives a H
2
performance of 9.54e-2 after solv-
ing 28 subproblems.
4.1.3 Chemical reactor (CR)
A model description for the chemical reactor can be
found in (Hung and MacFarlane, 1982). Both H
and H
2
static controllers are computed with the spec-
tral SDP method, which requires solving respectively
28 and 41 subproblems. The resulting gains are:
K
=
10.17039 31.54413
26.01021 100.00000
,
Table 1: Static output feedback synthesis. n and m
give size of (1), iter counts instances of (14), cpu in
seconds, index gives the corresponding H
or H
2
performance.
Problem n m iter cpu index
TA
(H
)
51 30 29 39 2.22
VH
(H
)
13 12 29 0.3 1.59e-1
VH
(H
2
)
16 13 28 0.5 9.54e-2
CR
(H
)
15 18 28 1.3 1.17
CR
(H
2
)
25 19 41 24.6 1.94
PA
(H
)
19 13 30 13.5 2.19e-3
with H
-performance T
w,z
(s)
=1.17.
K
2
=
0.35714 2.62418
2.58159 0.77642
,
with H
2
performance T
w,z
(s)
2
=1.94.
4.1.4 Piezoelectric actuator (PA)
In this section, we consider the design of a static con-
troller for a piezoelectric actuator system. A compre-
hensive presentation of this system can be found in
(Chen, 1998). Our algorithm computes a static H
controller after solving 30 instances of (14):
K
=[
0.09659 1.45023 100.00
] .
The computed H
-performance 2.19e-3 improves
significantly over the value 0.785 given in (Leibfritz,
1998).
4.2 Miscellaneous examples
4.2.1 Sparse linear constant output-feedback
design
In many applications it is of importance to design
sparse controllers with a small number of nonzero en-
tries. This is for instance the case when it is manda-
tory to limit the number of arithmetic operations in
the control law, or when reducing the number of sen-
sor/actuator devices is critical. In such a case one
may attempt to synthesis a controller with an imposed
sparsity structure. A more flexible idea is to let the
problem find its own sparsity structure by using a
1
-
norm cost function. We solve the following (BMI)
NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS: APPLICATIONS TO CONTROL
DESIGN
245
problem:
min
1in
1jn
, |K
ij
|,
s. t. P I
(A + BKC)
T
P + P (A + BKC)
−2 αP I
where is the usual threshold to adapt strict inequali-
ties to the form (1) or (8). This corresponds to design-
ing a sparse linear constant output feedback control
u = Ky for the system ˙x = Ax + Bu, y = Cx with
an imposed decay rate of at least α in the closed-loop
system.
Minimizing the
1
-norm of the feedback gain is a
good heuristic for obtaining sparse feedback gain ma-
trices. This may be explained by the fact that the func-
tion φ(t)=|t| is the natural continuous interpolant
of the counting function c(t)=t. Minimizing the
function
ij
K
ij
or even counting the nonzero en-
tries K
ij
in K is what we would really like to do, but
this is a discrete optimization problem subject to ma-
trix inequality constraints.
There are many more applications to this heuristic.
Finding sparse feedback gain matrices is also a good
way of solving the actuator/sensor placement or con-
troller topology design problems. In our example the
method works nicely and a sparsity pattern is appar-
ent. If small but nonzero values occur, one may set
these values to zero if below a certain threshold, de-
rive a sparsity pattern, and attempt a structured design
with that pattern imposed. This leads to another class
of BMI problems.
In our example, the system is initially unstable with
rate of α
0
= 0.2832. This value is obtained by com-
puting the smallest negative real part of the eigenval-
ues of A (data available from (Hassibi et al., 1999)).
In a first test similar to the one in (Hassibi et al.,
1999), we seek a sparse K so that the decay rate of the
closed-loop system is not less than 0.35. We obtain
the following gain matrix:
K
=
0.2278 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.3269 0.0000 0.0000
,
with only two non-zeros elements. The path-
following method presented in (Hassibi et al., 1999)
gives 3 non-zeros elements. The spectral method
gives once again a significant improvement, while
maintaining the closed-loop stability.
If now in a second test a higher decay rate (α
0.5) (and so a better damped closed-loop system), is
required, our method computes
K
=
0.1072 0.0000 0.0000 0.0000
0.0224 0.0000 0.0000 0.0000
0.3547 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.0186 0.1502 0.0000 0.0000
,
with 5 non-zeros entries. As expected, the sparsity
decreases when α increases.
4.2.2 Simultaneous state-feedback stabilization
with limits on feedback gains
One attempts here to stabilize three different lin-
ear systems using a common linear constant state-
feedback law with limits on the feedback gains.
Specifically, suppose that:
˙x = A
k
x + B
k
u, u = Kx,k =1, 2, 3.
The goal is to compute K such that |K
ij
|≤K
ij,max
and the three closed-loop systems:
˙x =(A
k
+ B
k
K ) x,k =1, 2, 3 ,
are stable. Under these constraints the worst decay-
rate of the closed-loop systems is maximized.
Proposition 4.3 The feedback gain K stabilizes the
above systems simultaneously if and only if there exist
P
k
0, α
k
> 0, k =1, 2, 3 such that
(A
k
+ B
k
K)
T
P
k
+ P
k
(A
k
+ B
k
K) ≺−2α
k
P
k
,
k =1, 2, 3.
Introducing our usual threshold parameter >0,
we have to check whether the value of the following
BMI program is positive:
maximize min
k=1,2,3
α
k
subject to |K
ij
|≤K
ij,max
,
(A
k
+ B
k
K)
T
P
k
+ P
k
(A
k
+ B
k
K)
−2 α
k
P
k
I
P
k
I, k =1, 2, 3.
By looking in (Hassibi et al., 1999) and computing
the corresponding eigenvalues, the reader will notice
that all three systems are unstable. We chose, as pro-
posed by Hassibi, How and Boyd, K
ij,max
=50. Ac-
tually this optimization problem turned out one of the
hardest for our algorithm. We obtain a solution after
68 unconstrained minimization subproblems, in 7.3
seconds. Actually a restart from a stabilizing point
was done here, which explains the unexpectedly high
number (68 = 32 + 36) of subproblems iterations.
On the other hand, this experiment gave the best
improvement among our present tests. The worst
decay-rate was increased to α =3.42, an improve-
ment of more than 200% over the value 1.05 obtained
in (Hassibi et al., 1999). Notice that none of the gain
bounds was active at the optimum:
K
=
42.4186 38.6457 17.4347
37.4186 15.9380 33.6942
.
ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION
246
Table 2: Miscellaneaous examples.
Problem n m iter cpu
Sparse design 1 56 51 29 2.9
Sparse design 2 56 51 19 7.3
Simult. stab. 28 34 68 7.3
H
2
/H
17 27 25 1.1
4.2.3 H
2
/H
controller design
Mixed H
2
/H
-performance indices represent one of
those prominent cases where the projection Lemma
fails and direct BMI-techniques are required. Here we
consider the case where a static state-feedback con-
troller u = Kx for the open loop system
˙x = Ax + B
1
u + B
2
w
z
1
= C
1
x + D
11
u
z
2
= C
2
x + D
21
u
is sought for. The purpose is to find K such that the
H
2
norm w z
2
is minimized, while the H
-norm
w z
1
is fixed below a performance level γ. After
introducing our usual threshold >0, the following
matrix inequality optimization problem arises:
minimize η
2
subject to
(A + B
1
K)
T
P
1
+P
1
(A + B
1
K)
+(C
1
+ D
11
K)
T
(C
1
+ D
11
K)
B
2
T
P
1
γ
2
I
−I
(A + B
1
K)
T
P
2
+ P
2
(A + B
1
K)
B
T
2
P
2
I
−I
P
2
C
2
Z
I, Tr (Z) η
2
,P
1
I , P
2
I.
The first block contains multilinear terms, which may
be converted to BMI-from using a Schur complement.
With γ =2, the algorithm needs 25 instances of (14)
to compute the compensator:
K
2,
=[
1.1500 0.4685 1.1500
] ,
with corresponding H
2
-norm 0.8384. In this case the
improvement over the result in (Hassibi et al., 1999)
is negligible, which we interpret in the sense that this
result was already close to the global optimum. (Nota
bene, one cannot have any guarantee as to the truth of
this statement. Indeed, in this example we found that
many other local minima could be computed).
5 CONCLUSION
A spectral penalty augmented Lagrangian method for
matrix inequality constrained nonlinear optimization
programs was presented and applied to a number of
small to medium-size test problems in automatic con-
trol.
The algorithm performs robustly if parameters are
carefully tuned. This refers to the choice of the initial
penalty value, p, the control of the condition number
of the tangent problem Hessian
2
L(x, V, p) in (14),
and to the stopping tests (complementarity, stationar-
ity), where we followed the lines of (Henrion et al.,
2003). For the updates p p
+
and V V
+
we
found that it is vital to avoid drastic changes, since the
new subproblem (14) may be much more difficult to
solve. This is particularly so for the nonconvex BMI
case.
We remind the reader that similar to our previous
approaches, the proposed algorithm is a local opti-
mization method, which gives no certificate as to find-
ing the global optimal solution of the control prob-
lem. If fact, such a method may in principle even
fail completely by converging to a local minimum of
constraint violation. For all that, our general expe-
rience, confirmed in the present testing, is that this
type of failure rarely occurs. We consider the local
approach largely superior for instance to a very typi-
cal class of methods in automatic control, where hard
problems are solved by a tower of LMI problems with
growing size, N, where a solution is guaranteed as
N →∞. Such a situation is often useless, since the
size of the LMI problem needed to solve the under-
lying hard problem is beyond reach, and since a prior
estimate is not available or too pessimistic. The major
advantage of local methods is that the subproblems,
on whose succession the problem rests, all have the
same dimension.
We finally point out that the local convergence of
the algorithm will be proved in a forthcoming article.
REFERENCES
Anderson, B. and Vongpanitlerd, S. (1973). Network analy-
sis and synthesis: a modern systems theory approach.
Prentice-Hall.
Apkarian, P., Noll, D., Thevenet, J. B., and Tuan, H. D.
(2003a). A Spectral Quadratic-SDP Method with Ap-
plications to Fixed-Order H
2
and H
Synthesis. sub-
mitted. Rapport Interne , MIP, UMR 5640, Maths.
Dept. - Paul Sabatier University.
Apkarian, P., Noll, D., and Tuan, H. D. (2002). Fixed-
order H
control design via an augmented La-
grangian method . Rapport Interne 02-13, MIP,
UMR 5640, Maths. Dept. - Paul Sabatier University .
http://mip.ups-tlse.fr/publi/publi.html.
Apkarian, P., Noll, D., and Tuan, H. D. (2003b). Fixed-
order H
control design via an augmented La-
grangian method . Int. J. Robust and Nonlinear Con-
trol, 13(12):1137–1148.
NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS: APPLICATIONS TO CONTROL
DESIGN
247
Ben-Tal, A. and Zibulevsky, M. (1997). Penalty/barrier
multiplier methods for convex programming prob-
lems. SIAM J. on Optimization, 7:347–366.
Bertsekas, D. P. (1982.). Constrained optimization and La-
grange multiplier methods. Academic Press, London.
Chen, B. M. (1998). H
Control and Its Applications, vol-
ume 235 of Lectures Notes in Control and Informa-
tion Sciences. Springer Verlag, New York, Heidelberg,
Berlin.
Conn, A. R., Gould, N., and Toint, P. L. (1991). A globally
convergent augmented Lagrangian algorithm for opti-
mization with general constraints and simple bounds.
SIAM J. Numer. Anal., 28(2):545 572.
Conn, A. R., Gould, N. I. M., Sartenaer, A., and Toint, P. L.
(1993a). Global Convergence of two Augmented La-
grangian Algorithms for Optimization with a Com-
bination of General Equality and Linear Constraints.
Technical Report TR/PA/93/26, CERFACS, Toulouse,
France.
Conn, A. R., Gould, N. I. M., Sartenaer, A., and Toint, P. L.
(1993b). Local convergence properties of two aug-
mented Lagrangian algorithms for optimization with a
combination of general equality and linear constraints.
Technical Report TR/PA/93/27, CERFACS, Toulouse,
France.
Conn, A. R., Gould, N. I. M., Sartenaer, A., and Toint, P. L.
(1996). Convergence properties of an augmented La-
grangian algorithm for optimization with a combina-
tion of general equality and linear constraints. SIAM
J. on Optimization, 6(3):674 703.
Fares, B., Apkarian, P., and Noll, D. (2000). An Augmented
Lagrangian Method for a Class of LMI-Constrained
Problems in Robust Control Theory. In Proc. Amer-
ican Control Conf., pages 3702–3705, Chicago, Illi-
nois.
Fares, B., Apkarian, P., and Noll, D. (2001). An Augmented
Lagrangian Method for a Class of LMI-Constrained
Problems in Robust Control Theory. Int. J. Control,
74(4):348–360.
Fares, B., Noll, D., and Apkarian, P. (2002). Robust Control
via Sequential Semidefinite Programming. SIAM J. on
Control and Optimization, 40(6):1791–1820.
Gangsaas, D., Bruce, K., Blight, J., and Ly, U.-L. (1986).
Application of modern synthesis to aircraft control:
Three case studies. IEEE Trans. Aut. Control,AC-
31(11):995–1014.
Hassibi, A., How, J., and Boyd, S. (1999). A pathfollowing
method for solving bmi problems in control. In Proc.
American Control Conf., pages 1385–1389.
Henrion, D., M.Kocvara, and Stingl, M. (2003). Solving
simultaneous stabilization BMI problems with PEN-
NON. In IFIP Conference on System Modeling and
Optimization, volume 7, Sophia Antipolis, France.
Hestenes, M. R. (1969). Multiplier and gradient method. J.
Optim. Theory Appl., 4:303 320.
Hung, Y. S. and MacFarlane, A. G. J. (1982). Multivari-
able feedback: A classical approach. Lectures Notes
in Control and Information Sciences. Springer Verlag,
New York, Heidelberg, Berlin.
Keel, L. H., Bhattacharyya, S. P., and Howze, J. W. (1988).
Robust control with structured perturbations. IEEE
Trans. Aut. Control, 36:68–77.
Kocvara, M. and Stingl, M. (2003). A Code for Convex
Nonlinear and Semidefinite Programming. Optimiza-
tion Methods and Software, 18(3):317–333.
Leibfritz, F. (1998). Computational design of stabiliz-
ing static output feedback controller. Rapports 1 et
2 Mathematik/Informatik, Forschungsbericht 99-02,
Universitt Trier.
Lewis, A. (1996). Derivatives of spectral functions. Math-
ematics of Operations Research, 21:576–588.
Lewis, A. (2001). Twice differentiable spectral functions.
SIAM J. on Matrix Analysis and Applications, 23:368–
386.
Lewis, A. and S.Sendov, H. (2002). Quadratic expan-
sions of spectral functions. Linear Algebra and Appl.,
340:97–121.
Lin, C. and More, J. (1998). Newton’s method for large
bound–constrained optimization problems. Techni-
cal Report ANL/MCS-P724–0898, Mathematics and
Computer Sciences Division, Argonne National Lab-
oratory.
Mosheyev, L. and Zibulevsky, M. (2000). Penalty/barrier
multiplier algorithm for semidefinite programming.
Optimization Methods and Software, 13(4):235–261.
Noll, D., Torki, M., and Apkarian, P. (2002). Partially
Augmented Lagrangian Method for Matrix Inequality
Constraints. submitted. Rapport Interne , MIP, UMR
5640, Maths. Dept. - Paul Sabatier University.
Powell, M. J. D. (1969). A method for nonlinear constraints
in minimization problem. In Fletcher, R., editor, Op-
timization. Academic Press, London, New York.
Schnabel, R. B. and Eskow, E. (1999). A revised modified
cholesky factorization algorithm. SIAM J. on Opti-
mization, 9(4):1135–1148.
Shapiro, A. (2002). On differentiability of symmetric ma-
trix valued functions. School of Industrial and Sys-
tems Engineering, Georgia Institute of Technology.
Preprint.
Zibulevsky, M. (1996). Penalty/Barrier Multiplier Methods
for Large-Scale Nonlinear and Semidefinite Program-
ming. Ph. D. Thesis, Technion Isral Institute of Tech-
nology.
ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION
248