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 semideﬁnite programming method to solve prob-

lems with bilinear matrix inequality (BMI) constraints. Such optimization programs arise frequently in au-

tomatic control and are difﬁcult 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 difﬁcult 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

ﬁltering problems, synthesis of structured or reduced-

order feedback controllers, simultaneous stabilization

problems and many others. Formally, these problems

may be described as

minimize c

x, x ∈ R

subject to B(x)  0,

(1)

where  0 means negative semi-deﬁnite and

B(x)=A



i=1



1≤i<j≤n

(2)

is a bilinear matrix-valued operator with data

∈ S

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 semideﬁnite 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.

Signiﬁcant progress was achieved by two recent

approaches, both using successions of SDP tangent

subproblems. In (Fares et al., 2002), a sequen-

tial semideﬁnite 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 difﬁcult

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

x, x ∈ R

subject to A(x)  0

g(x)=0.

(3)

Here A is an afﬁne operator (i.e., (2) with B

=0)

with values in S

, and g : R

→ R

incorporates a

ﬁnite 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 semideﬁnite programming,

(Zibulevsky, 1996; Ben-Tal and Zibulevsky, 1997)

propose a modiﬁed 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

 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 signiﬁcantly reduced) with genuine second

order strategies. In particular, trust region techniques

do not require convexiﬁcation 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-

ﬁculties 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

∞

or mixed H

∞

performance indices. Two

other difﬁcult 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 ﬁrst 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

denote the set

of m × m symmetric matrices, M

the transpose

of the matrix M and Tr M its trace. We equip S

with the euclidean scalar product X, Y  = X •Y =

Tr (XY ). For symmetric matrices , M  N means

that M − N is positive deﬁnite and M  N means

that M − N is positive semi-deﬁnite. 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

⊗ A)vecX, (4)

Tr (AB)=vec

vec B, (5)

which hold for matrices A, X and B of compatible di-

mensions. The Hadamard or Schur product is deﬁned

A ◦ B =((A

)). (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

subject to F(x)  0,

(8)

where f is a class C

function, F : R

→ S

a class

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 afﬁne 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 deﬁned

F (x, p)=f(x)+TrΦ

(F(x)), (9)

where ϕ

: R → R is a parametrized family of scalar

functions, which generates a family of matrix-valued

operators Φ

: S

→ S

upon deﬁning:

(X):=Sdiag [ϕ

(λ

(X))] S

. (10)

Here λ

(X) stands for the ith eigenvalue of X ∈

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)+



i=1

(λ

(F(x)) . (11)

As ϕ

is chosen strictly increasing and satisﬁes

(0) = 0, each of the following programs

(parametrized through p>0)

minimize f(x),x∈ R

subject to Φ

(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

F (x, p) to converge to a solution of (8).

It is well-known that pure penalty methods run into

numerical difﬁculties 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 deﬁne the

augmented Lagrangian function associated with the

matrix inequality constraints in (1) as

L(x, V, p)=f(x)+TrΦ

F(x)V ) , (13)

= f (x)+



i=1

(λ

F(x)V )).

In this expression, the matrix variable V has the

same dimension as F(x) ∈ S

and serves as a fac-

tor of the Lagrange multiplier variable U ∈ S

, U =

(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 ﬁrst-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 sufﬁcient 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 ﬁeld. 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 ϕ

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

, V

and a penalty

parameter p

> 0.

2. Optimization phase. For ﬁxed V

and p

solve

the unconstrained subproblem

minimize

x∈R

L(x, V

) (14)

Let x

j+1

be the solution. Use the previous iterate

as a starting value for the inner optimization.

3. Update penalty and multiplier. Apply ﬁrst-

order rule to estimate Lagrange multiplier :

j+1

= V



diagϕ





·F(x

j+1

))] S

, (15)

where S diagonalizes V

F(x

j+1

Update the penalty parameter using :

j+1



ρp

, if λ

max

(0, F(x

j+1

))

>γλ

max

(0, F(x

))

, 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 ϕ

(t)=

pϕ

(t/p) where

(t)=



t +

if t ≥−

−

log(−2t) −

if t<−

(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 Φ

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

(x),...,c

(x)]. Then

the Lagrange multiplier U and its factor V may also

be restricted to the space of diagonal matrices, U =

diag u

, V = diag v

, and we recover the situation of

classical polyhedral constraints. Switching to a more

convenient notation, the problem becomes

minimize f(x)

subject to c

(x) ≤ 0,j=1,...,m

With u

= v

, we obtain the analogue of the aug-

mented Lagrangian (13):

L(x, v, p)=f(x)+



i=1

pϕ

(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)



i=1



(x)/p)v

∇c

(x).

If we compare with the Lagrangian :

L(x, u)=f(x)+u

c(x),

and its gradient :

∇L(x, u)=∇f(x)+



i=1

∇c

(x),

the following update formula readily appears:

= v



)/p),

the scalar analogue of (15). Suppose now that ϕ

has

the form (17). Then for sufﬁciently small v

(x) we

expect v

(x)/p > −

,so

pϕ

(x)/p)=p



(x)/p

(x)



= u

(x)+

/p)c

(x)

If we remember the form of the classical augmented

Lagrangian for inequality constraints (cf. (Bertsekas,

1982))

L(x, u, µ)=f(x)



i=1



max [0,u

+ c

(x)/µ]

− u



we can see that the term v

/p takes the place of

the penalty parameter 1/µ as long as the v



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 ϕ

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 ﬁrst-

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

deﬁnition from (Lewis and S.Sendov, 2002).

Deﬁnition 3.1 Let λ : S

→ R

denote the eigen-

value map λ(X):=(λ

(X),...,λ

(X)). For a

symmetric function ψ : R

→ R, the spectral func-

tion Ψ associated with ψ is deﬁned as Ψ:S

→ 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 deﬁned on

an open subset D of R

. Let Ψ=ψ◦λ be the spectral

function associated with ψ. Let X ∈ S

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

for every orthogonal matrix S satisfying

X = S diagλ(X) S

. 

The gradient ∇Ψ(X) is a (dual) element of S

which we distinguish from the differential, dΨ(X).

The latter acts as a linear form on tangent vectors

dX ∈ S

dΨ(X)[dX]=∇Ψ(X) • dX

=Tr



S diag ∇ψ (λ(X))S



It is now straightforward to compute the ﬁrst deriva-

tive of a composite spectral functions Ψ ◦F, where

F : R

→ S

is sufﬁciently smooth. We obtain

∇(Ψ ◦F)(x)=dF(x)

∗

[∇Ψ(F(x))] ∈ R

, (18)

where dF(x) is a linear operator mapping R

→ S

and dF(x)

∗

its adjoint, mapping S

→ R

. With the

standing notation

F := F(x),F

∂F(x)

∂x

∈ S

we obtain the representations

dF(x)[δx]=



i=1

δx

dF(x)

∗

[dX]=(F

• dX,...,F

• dX) .

Then the ith element of the gradient is

(∇(Ψ ◦F)(x))

=Tr



Sdiag∇ψ (λ(F )) S



. (19)

Observe that F

= A



j<i



i<j

if B is of the form (2), F

= A

for afﬁne operators

We are now ready for the ﬁrst derivative of the aug-

mented Lagrangian function L(x, V, p). We consider

a special class of symmetric functions ψ : R

→ R,

ψ(x)=



i=1

ϕ(x

), (20)

generated by a single scalar function ϕ : R → R.

Here the spectral function Ψ=ψ ◦ λ is of the form

Ψ=TrΦwith Φ:S

→ S

the matrix-valued

operator

Φ(X)=S diag ϕ (λ(X)) S

S being an orthogonal matrix diagonalizing X. Using

, ψ

and Ψ

instead of ϕ, ψ and Ψ, the penalty term

in (13) takes the form

F(x)V )) = ψ

◦ λ(V

F(x)V ))



i=1

(λ

F(x)V ))

=TrΦ

F(x)V ),

where Φ

is given in (10). Now we apply Lemma

3.2, respectively formula (18) or even (19) to Ψ

Tr Φ

and the operator



F(x)=V

F(x)V , where

we observe en passant that



F(x)

∗

[X]=d





(x)



∗

[X]

= dF(x)

∗

[VXV

or put differently,



= V

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

L(x, V, p)=

∂

∂x

f(x)

+Tr



VSdiag [ϕ



(λ

FV))] S



In vector notation

∇L(x, V, p)=∇f(x)+[vec (F

),...,vec (F

)]

· vec



VSdiag [ϕ



(λ

FV))] S



where S is orthogonal and diagonalizes V

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

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

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 λ

= λ

(X) ∈ I for every i.

Moreover, if we deﬁne the matrix ϕ



[1]

= ϕ



[1]

(X) ∈

as:



[1]





(λ

)−ϕ



(λ

)

−λ

if λ

= λ



(λ

) if λ

= λ

then the following formula deﬁnes the second order

derivative of Ψ at X

Ψ(X))[dX, dX]=Tr(S

dXSϕ



[1]

◦{S

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

Ψ(X)[dX, dX] of

Ψ and d

F(x)[δx,δx] of F amounts to:

(Ψ ◦F)(x)[δx, δx]

= dΨ(F(x))



F(x)[δx,δx]



+ d

Ψ(F(x)) [dF(x)[δx],dF(x)[δx]] . (21)

For bilinear B, the second order form d

B(x) is of

course independent of x and given as

B(x)[δx, δx]=



i<j

δx

In particular, d

A =0for afﬁne 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

=TrΦ

generated by ϕ

and to the operator



F(x)=V

F(x)V . As we have already seen, for



F we have the formula



F(x)[δx]=



i=1

Vδx

= V



i=1

δx

= V

dFV.

Therefore, using Lemma 3.5, the second term in (21)

becomes

FV)[V

dF V,V

dF V ]

=Tr



dF VSϕ



[1]

◦{S

dF VS}



where S is orthogonal and diagonalizes V

FV and



[1]

= ϕ



[1]

FV). This may now be transformed

to vector-matrix form using the operator vec and the

relations (4) - (7).

FV)[V

dF V,V

dF V ]

=vec

((VS)

dF (VS))

· vec (ϕ



[1]

◦{(VS)

dF (VS)})

=vec

(dF) K

vec (dF) , (22)

where the matrix K

=[(VS) ⊗ (VS)]diag



vec ϕ



[1]



· [(VS)

⊗ (VS)

Using dF =



i=1

δx

, the last line in (22) gives

FV)[V

dF V,V

dF V ]

= δx

[vec (F

),...,vec (F

)]

· [vec (F

),...,vec (F

)] δx,

so we are led to retain the symmetric n × n matrix

=[vec(F

),...,vec (F

)]

· [vec (F

),...,vec (F

)] . (23)

Let us now consider the second part of the

Hessian ∇

L(x, V, p), coming from the term

dΨ(F(x))



F(x)[δx,δx]



in (21). First observe

that trivially



F(x)[δx,δx]=d

FV }(x)[δx, δx]

= V

F(x)[δx,δx] V.

Hence

dΨ

FV)



F(x)[δx,δx]V



=Tr



S diag ϕ





FV)



·V

F(x)[δx,δx]V



=Tr



VSdiag ϕ





FV)



·d

F(x)[δx,δx]



. (24)

Introducing matrices F

= F

(x) ∈ S

to repre-

sent

F(x)[δx,δx]=



i,j=1

δx

relation (24) may be written as

dΨ

FV)



F(x)[δx,δx]V





i,j=1

δx



VSdiag ϕ





·FV)) S



Hence the interest to deﬁne a second symmetric n ×n

matrix H

)

=Tr



VSdiagϕ





FV)



(25)

where as before S is orthogonal and diagonalizes

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

∇

L(x, V, p)=∇

f(x)+H

+ H

where H

, given by (23), is positive semideﬁnite. H

is given by (25).

Proof. The only element which remains to be proved

is non-negativity of H

, which hinges on the non-

negativity of K

, and hence on that of ϕ



[1]

. The latter

is a consequence of the convexity of ϕ

, hence the re-

sult. 

ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION

242

3.2 Multiplier update rule

The ﬁrst-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

the current La-

grange multiplier estimate. A new multiplier U

and therefore a new factor estimate V

with U

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

F(x)) =

f(x)+Tr(V

F(x)V ), and computing its gradient

at x

and U

= V

implies

∇L(x

)=∇f(x

)

+[vec(F

),...,vec (F

)]

vec (V

Comparing with

∇L(x

,p)

= ∇f(x

) + [vec (F

),...,vec (F

)]

· vec



VSdiag ϕ



(λ

F(x

)V )S



suggests the update rule

= VS[diag ϕ



(λ

F(x

)V ))] (VS)

(27)

where S is orthogonal and diagonalizes V

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

),...,vec (F

)]

, but the limit-

ing case gives more information. Notice that equality

∇L(¯x,

U)=∇L(¯x,

V,p)

holds at a Karush-Kuhn-Tucker pair (¯x,

U) and

U =

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,

F(¯x)

V =0. With ϕ



(0) = 1 for every p>

0, the right hand side in (27) is

V =

V ,

hence equals the left hand side. We may therefore

read (27) as some sort of ﬁxed-point relation, which

is satisﬁed 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 deﬁnite, so V

could for in-

stance be chosen as a Cholesky factor.

3.3 Solving the subproblem -

implementational issues

Efﬁcient minimization of L(x, V

) for ﬁxed V

, p

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) ﬁnd good local minima, leaving

the bad ones behind. This additional beneﬁt is in

contrast with what line search methods achieve, and

is explained to some extent by the fact that, at least

over the horizon speciﬁed 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 signiﬁcant 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 ﬁrst-

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 efﬁcient preconditioner, or

to compute the Newton direction in the second case.

When the Hessian is close to becoming indeﬁnite, a

modiﬁed Cholesky factorization is used (∇

L + E =

, 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 modiﬁed 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 ﬁrst 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 difﬁcult to predict in advance,

as the physical meaning of the Lyapunov variables

is often lacking. Properly scaling these variables is

therefore difﬁcult 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 difﬁculties and a method to balance

data is extremely important.

4 NUMERICAL EXAMPLES

In this section we test our code on a variety of difﬁcult

applications in automatic control. The ﬁrst 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 ﬁnally

multi-objective 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

or mixed per-

formances previously published in the literature are

signiﬁcantly improved. We always start our algorithm

at x

=0in order to capture a realistic scenario,

where no such prior information is available.

4.1 Static output feedback controller

synthesis

We ﬁrst 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









, (28)

where

• x ∈ R

is the state vector,

• u ∈ R

is the vector of control inputs,

• w ∈ R

is a vector of exogenous inputs,

• y ∈ R

is the vector of measurements,

• z ∈ R

is the controlled or performance vector.

• D

=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)

K(s)

P(s)

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 inﬁnity.

• performance: the H

∞

norm T

w,z

(s)

∞

, respec-

tively the H

norm T

w,z

(s)

, is minimized

where the closed-loop transfer T

w,z

(s) is described

as :



˙x =(A + B

)x +(B

+ B

z =(C

+ D

)x +(D

+ D

)w.

Recall that in the H

case, some feedthrough terms

must be nonexistent in order for the H

performance

to be well deﬁned. We have to assume D

0,D

=0for the plant in (28). This guaran-

teed, both static H

∞

and H

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

such that:



(A+B

)

X+ ∗∗ ∗

)

X −γI ∗

)(D

) −γI



≺ 0

(30)

 0. (31)

Proposition 4.2 A stabilizing static output feedback

controller K with H

performance T

w,z

(s)

≤

√

exists provided there exist X, M ∈ S

such that:

ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION

244



(A+B

)

X+ ∗∗

) −I



≺ 0 (32)



M ∗

X (B

) 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):

∞

⎡

⎢

⎣

0.69788

−0.64050

−0.83794

0.09769

1.57062

⎤

⎥

⎦

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 ﬂight conditions. This leads to a fourth-order

model. Both H

∞

and H

static controllers are com-

puted with the spectral SDP method. The algorithm

converges after solving 29 tangent problems (14):

∞



0.50750

10.00000



The ﬁnal H

∞

performance is 1.59e-1, which is about

40% improvement over the performance obtained in

(Leibfritz, 1998).



0.13105

5.95163



which gives a H

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

static controllers are computed with the spec-

tral SDP method, which requires solving respectively

28 and 41 subproblems. The resulting gains are:

∞



−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

performance.

Problem n m iter cpu index

∞

)

51 30 29 39 2.22

∞

)

13 12 29 0.3 1.59e-1

)

16 13 28 0.5 9.54e-2

∞

)

15 18 28 1.3 1.17

)

25 19 41 24.6 1.94

∞

)

19 13 30 13.5 2.19e-3

with H

∞

-performance T

w,z

(s)

∞

=1.17.



0.35714 −2.62418

2.58159 0.77642



with H

performance T

w,z

(s)

=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):

∞

0.09659 −1.45023 −100.00

] .

The computed H

∞

-performance 2.19e-3 improves

signiﬁcantly 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 ﬂexible idea is to let the

problem ﬁnd its own sparsity structure by using a 

norm cost function. We solve the following (BMI)

NON LINEAR SPECTRAL SDP METHOD FOR BMI-CONSTRAINED PROBLEMS: APPLICATIONS TO CONTROL

DESIGN

245

problem:

min



1≤i≤n

1≤j≤n

, |K

s. t. P  I

(A + BKC)

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 

-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



K

 or even counting the nonzero en-

tries K

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.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 ﬁrst 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:

∞

⎡

⎢

⎣

0.2278 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 signiﬁcant 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

∞

⎡

⎢

⎣

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.

Speciﬁcally, suppose that:

˙x = A

x + B

u, u = Kx,k =1, 2, 3.

The goal is to compute K such that |K

|≤K

ij,max

and the three closed-loop systems:

˙x =(A

+ B

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

 0, α

> 0, k =1, 2, 3 such that

+ B

+ P

+ B

K) ≺−2α

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

subject to |K

|≤K

ij,max

+ B

+ P

+ B

−2 α

− I

 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:

∞



−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

∞

17 27 25 1.1

4.2.3 H

∞

controller design

Mixed 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

u + B

= C

x + D

= C

x + D

is sought for. The purpose is to ﬁnd K such that the

norm w → z

is minimized, while the H

∞

-norm

w → z

is ﬁxed below a performance level γ. After

introducing our usual threshold >0, the following

matrix inequality optimization problem arises:

minimize η

subject to

⎡

⎢

⎣

⎛

⎝

(A + B

+(C

+ D

⎞

⎠

∗

− γ

⎤

⎥

⎦

−I



(A + B

+ P

(A + B

K) ∗

− I



−I



∗



 I, Tr (Z) ≤ η

 I , P

 I.

The ﬁrst 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:

2,∞

1.1500 −0.4685 −1.1500

] ,

with corresponding H

-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 ∇

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

found that it is vital to avoid drastic changes, since the

new subproblem (14) may be much more difﬁcult 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 certiﬁcate as to ﬁnd-

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, conﬁrmed 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 ﬁnally 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

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 Semideﬁnite 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 Semideﬁnite 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 semideﬁnite 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 modiﬁed

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 Semideﬁnite Program-

ming. Ph. D. Thesis, Technion Isral Institute of Tech-

nology.

ICINCO 2004 - INTELLIGENT CONTROL SYSTEMS AND OPTIMIZATION

248