1. No category

Computing the Maximum Volume Inscribed
Ellipsoid of a Polytopic Projection
Jianzhe Zhen∗ and Dick den Hertog†
Department of Econometrics and Operations Research
Tilburg University
Tilburg, The Netherlands
Abstract
This paper introduces a method for computing the maximum volume inscribed ellipsoid and
k-ball of a projected polytope. It is known that deriving an explicit description of a projected
polytope is NP-hard. By using adjustable robust optimization techniques, we construct a computationally tractable method that does not require an explicit description of the projection.
The obtained centers of the projected polytope are considered as the robust solutions, e.g., for
design centering problems. We perform numerical experiments on a simple polytope and a color
tube design problem. The color tube design problem demonstrates that the obtained solutions
are much more robust than the nominal approach with a slight compromise on the objective
value. Some other potential applications include ellipsoidal approximations to polytopic sets,
nominal scenario recovery, and cutting-plane method.
Keywords: Maximum volume inscribed ellipsoid; chebyshev center; polytopic projection; adjustable
robust optimization.
JEL Classification: C44, C61, C63
1
Introduction
Let P be a polytope in Rn (not necessarily full-dimensional) defined by m linear inequalities
P=
x
x
∈ Rn | A
≤b ,
y
y
(1)
where x ∈ Rn1 , y ∈ Rn2 , n = n1 + n2 , A ∈ Rm×n , and b ∈ Rm . In this paper, we develop a method
to solve the following problem:
Problem. Given polytope P , find the centers of maximum volume inscribed ellipsoid (abbr., MVE
center) and k-ball (a.k.a., Chebyshev center, where k ∈ R denotes the dimension of the ball) of the
projected P onto the x-space without deriving an explicit description of the projection.
∗ Corresponding author. Email: [email protected]. The research of this author is supported by NWO Grant
613.001.208.
† Email: [email protected].
1
One possible application of our method is to find robust solutions for design centering problems.
Consider a manufacturing process in which the characteristics of a product are represented by a
vector x. The product is acceptable if x lies in the “region of acceptability" P defined by (1). Due
to unavoidable disturbances, e.g., implementation errors, in the manufacturing process, the actual
value of x will deviate from the nominal value in the design. Moreover, because of the nature of
the problem at hand or the result of reformulations (such as slack or surplus variables, or variables
introduced in order to eliminate piecewise-linear functions like | xi | or max{ xi , 0}, see Example
2.2), there exists a auxiliary variable y in P . The auxiliary variable y is certain. Since only the
variable x suffers from uncertainties, a product design x that is “most interior" or “most robust" in
the x-space of P is desired. We consider the MVE center and Chebyshev center of the projected P
onto the x-space as our robust solutions. In this paper, we develop a new method for computing
the MVE center and Chebyshev center of P only with respect to a subset of variables, e.g., the
variable x. This problem can be seen as a generalized version of the classical design centering or
tolerance design problems (see Graeb [2007], Hendrix et al. [1996]). The main difference between
our version and the classical one is that the classical version considers the case where all the
variables are uncertain, whereas, we deal with the case where only some of the variables are
uncertain. In Example 2.1 and 2.2, it is shown that the classical approach produces poor quality
robust solutions due to the existence of the auxiliary variables. We refer to §6 for more potential
applications of our method.
The reason we are focusing on the MVE center and Chebyshev center is that both centers possess
many appealing properties. Firstly, both of the centers are (relative) interior points and they
are invariant of the representation of the given convex body. A Chebyshev center is a point
that is farthest from the exterior of the polytope; the MVE center is unique (see Danzer et al.
[1957], Zaguskin [1958]), affine invariant and more centralized than the Chebyshev center. Other
popular centers are the geometric center and analytic center. Rademacher [2007] shows that
computing the geometric center is #P-hard. Although the analytic center is easy to compute, it
is not invariant of the representation of the given set. For more details on the geometric center and analytic center, we refer the interested reader to the book by Boyd and Vandenberghe [2004].
The classical method for finding the MVE center and Chebyshev center (with respect to all the
variables) of a full-dimensional polytope is well-established. These two centers can be computed in polynomial time by solving a linear programming (LP) problem and a semidefinite
programming (SDP) problem, respectively (see Boyd and Vandenberghe [2004]). The aim of this
paper is to find these two centers under a more general framework. Instead of finding the MVE
center and Chebyshev center of a given polytope (i.e., with respect to all the variables), we compute the centers of the projection of the polytope (i.e., only with respect to a subset of the variables).
One obvious way of determining the centers in the x-space of a given polytope P is as follows.
Firstly, we project the auxiliary variables onto the x-space, in other words, eliminate the auxiliary
variable y from the set. Then, we can use the classical method to find the centers in the projected
polytope. The elimination of the auxiliary variables can equivalently be seen as orthogonal projections of the original polytope. The pioneer method Fourier-Motzkin elimination, a.k.a. the FME
method, was first introduced in Fourier [1824] in 1827, and was rediscovered in Motzkin [1936]. In
each step of the iterative algorithm the dimension of the polytope is reduced by projecting it onto
a hyperplane. The FME method has been improved many times since. The numerical performance
of a variant of the FME method, and two other projection methods, i.e., the Extreme Point Method
2
(see Lassez [1990]) and Convex Hull Method (see Lassez and Lassez [1992]), are evaluated in
Huynh et al. [1992]. In Jones et al. [2004], a new algorithm for obtaining the projection of polytopes
is developed. Their algorithm is suited for problems in which the number of vertices far exceeds
the number of facets. A more recent development on the method for variable elimination in
systems of inequalities can be found in Chaharsooghi et al. [2011]. For all the existing methods, the
size of the representation of the projection and of the intermediate computations can be intractable
(see Example 2.2). Tiwary [2008] shows that deriving an explicit description of a projected polytope
is NP-hard.
Alternatively, we develop a new method for computing the MVE center and Chebyshev center of
the projected polytope without deriving an explicit description of the projection. To this end, we
employ adjustable robust optimization techniques introduced in Ben-Tal et al. [2004]. The main
advantage of our method compared to all the existing methods is that it does not require the
explicit description of the projection. The maximum volume inscribed ellipsoid and k-ball found
by our method may not be optimal. In order to evaluate the quality of our (under-)approximations,
we adapt the idea introduced in Hadjiyiannis et al. [2011] to obtain upper bounds on the optimal
solutions. Note that the Chebyshev center is a special case of the MVE center. Without loss
of generality, we will focus on elaborating our method via the MVE center and only treat the
Chebyshev center separately when it is necessary.
The rest of this paper is organized as follows. In §2, we illustrate the main difficulties of computing
the MVE center in the projections. In §3, we introduce our method for full-dimensional polytopes
and extend our method to general polytopes. We adapt the method of Hadjiyiannis et al. [2011] to
obtain upper bounds of the optimal ellipsoid or k-ball in §4. In §5, we evaluate the quality of our
solutions via numerical experiments. In §6, we conclude with some possible applications of our
new method.
2
Why is it difficult to compute the MVE center in a projection?
There are two main sources of difficulties to compute the MVE center and Chebyshev center in a
polytopic projection. The first source is that the size of the representation of the projection and of
the intermediate computations can be intractable (see Tiwary [2008]). The second source of the
difficulties is that the projected space may not be full-dimensional. Hence, the full-dimensional
ellipsoid in the projected space simply does not exist and the classical method cannot be applied.
In this section, we first briefly review the classical method for computing the two centers in
a full-dimensional polytope. Then, it is shown that deriving the explicit representation of the
projections can be computational intractable.
Suppose polytope P defined by (1) is full-dimensional, which means that there is no implicit
equality constraint. We can find the MVE center of P by solving the following semi-infinite
programming problem (see Boyd and Vandenberghe [2004]):
( P)
max
x,y,E
s.t.
log detE
aiT·
x
+ Eζ ≤ bi
y
3
∀ζ : ||ζ ||2 ≤ 1, ∀i,
Figure 1: The maximum volume inscribed ellipsoid of the set P in Example 2.1.
where ai· denote the i-th row of matrix A, E ∈ Sn with implicit constraint E 0, where Sn is the
set of n × n symmetric matrices. The matrix E models the shape and volume of a n-dimensional
ellipsoid with center ( x, y). This problem can also be seen as a robust optimization problem under
ellipsoidal uncertainty set (see Ben-Tal et al. [2009]). It can be reformulated into an equivalent SDP
problem (see Boyd and Vandenberghe [2004]):
( RP)
max
x,y,E
s.t.
log detE
aiT·
x
+ || Eai· ||2 ≤ bi
y
∀i.
The optimal (x, y) from ( RP) corresponds to the unique MVE center of the polytope P . In the
following example, we apply ( RP) to find the MVE center of a specific polytope.
Example 2.1. The maximum volume inscribed ellipsoid of a full-dimensional polyhedron. Suppose P is the full-dimensional polytope that is formed by the intersection of three half-planes:
P = {( x, y)| − 0.5x − y ≤ −9, 0.6x + y ≤ 10, − x − y ≤ −10} .
(2)
This polytope is plotted in Figure 1. The MVE center of P is obtained by solving the optimization
problem
2.70
−
1.43
( P). The MVE center of P is at (4, 7 31 ) and the optimal matrix E is
.
−1.43 1.04
Now, let us assume that we are only interested to find the MVE center of P with respect to a
subset of variables, e.g., the variable x. The projection of P onto the x-space can be represented as
follows:
x
n1
n2
Projx (P ) = x ∈ R | ∃y ∈ R : A
≤b ,
(3)
y
where Projx (·) denotes the projection of (·) onto the space that variable x resides in. One obvious
attempt to determine the MVE center of Projx (P ) only with respect to the variable x, may be as
4
follows
(mP)
max
x,y,E
s.t.
log detE
aiT·
x + Eζ
y
≤ bi
∀ζ : ||ζ ||2 ≤ 1, ∀i,
where E ∈ Sn1 . Problem (mP) maximizes the volume of the ellipsoid around x with respect to the
best possible y. The next example shows that both ( RP) and (mP) fail to find the optimal MVE
center of Projx (P ).
EXAMPLE 1 (continued). The maximum volume inscribed ellipsoid of a full-dimensional polyhedron. We apply (mP) to the polytope P defined by (2) and find a new center at (2.22, 8). The best
possible y is at 8, which corresponds to the longest horizontal line within the feasible region (see Figure 2).
Our interest only resides in the x-space of P . From Figure 1 and 2, it can readily be seen that the feasible x
lies in the interval [0, 10]. Clearly, its MVE center coincides with its Chebyshev center. The MVE center
(or Chebyshev center) is x = 5, because it is the value that is farthest away (with distance 5) from the
boundaries 0 and 10. Both ( RP) and (mP) fail to find the optimal x at 5.
The centers found by solving ( RP) and (mP) corresponds to a fixed feasible y. However, it is not
necessary to restrict to a fixed y. One way of relaxing the restrictions imposed by the variable y
is to eliminate y from the system by using elimination methods. This is the same as deriving an
equivalent representation of Projx (P ) without the variable y. Tiwary [2008] shows that deriving an
explicit description of a projected polytope is NP-hard. For all the existing projection methods, the
size of the representation of the projection and of the intermediate computations can be intractable.
This is illustrated via the following example.
Example 2.2. The complexity of variable elimination/polytopic projection. Let us consider the
following polytope:
(
)
P=
x ∈ Rn |
n
∑ | xi | ≤ 1
.
i =1
We convert P into the following polytopic representation:
(
)
n
n
n
e
P = P = x ∈ R | ∃y ∈ R : ∑ yi ≤ 1, xi ≤ yi , − xi ≤ yi , ∀i .
i =1
e . If one employs elimination methods (e.g.,
In total there are 2n extra constraints and n extra variables in P
Fourier-Motzkin elimination, see Fourier [1824], Motzkin [1936]) to eliminate the auxiliary variable y, the
resulting polytope would contain 2n constraints.
3
MVE Center of Projected Polytope
In this section, we first consider the case where the polytopic projection is full-dimensional. The
adjustable robust optimization techniques is employed to compute the MVE center and Chebyshev
center without deriving the explicit description of the projections. Then, we extend our new
method to general polytopes.
5
Figure 2: The center of P obtained from (mP).
3.1
Full-dimensional Case
Suppose the polytope P defined by (1) is full-dimensional. The variable y is an auxiliary variable.
There exists a n1 -dimensional ellipsoid in the projected polytope Projx (P ) defined by (3). Finding
the MVE center of Projx (P ) can be formulated as follows:
x + Eζ
max log detE | ∀(ζ : ||ζ ||2 ≤ 1), ∃y(ζ ) : A
≤b .
(4)
y(ζ )
x,E
This can be interpreted as an adjustable robust optimization problem (see Ben-Tal et al. [2004]),
in which x is a main variable (non-adjustable variable) and y is an auxiliary variable (adjustable
variable). The vector function y(ζ ) is called a decision rule. Problem (4) can be reformulated as an
adjustable robust optimization problem:
( AP)
max
x,y(ζ ),E
s.t.
log detE
aiT·
x + Eζ
y(ζ )
≤ bi
∀ζ : ||ζ ||2 ≤ 1, ∀i,
where E ∈ Sn1 and ζ ∈ Rn1 . Ben-Tal et al. [2004] show that in general, solving such a problem is
NP-hard. In order to be computationally tractable, we need to restrict the vector function y(ζ ) to a
given class. Despite that linear functions may not be optimal, it appears that such decision rules
perform well in practice (see examples in Ben-Tal et al. [2004]). Let the decision rule y(ζ ) be linear,
i.e.,
y = u + Vζ,
where the coefficients u ∈ Rn2 and V ∈ Rn2 ×n1 are optimization variables. Let us substitute the
linear decision rule of y in ( AP) to get the affinely adjustable robust formulation:
( AAP)
max
x,u,E,V
s.t.
log detE
x + Eζ
aiT·
≤ bi
u + Vζ
6
∀ζ : ||ζ ||2 ≤ 1, ∀i.
Problem ( AAP) is a semi-infinite programming problem that approximates the MVE center of
Projx (P ). The tractable robust reformulation of ( AAP) is as follows (see Boyd and Vandenberghe
[2004]):
( AARP)
max
x,u,E,V
s.t.
log detE
aiT·
T x
E
+ a i · ≤ bi
u
V
∀i,
2
E
where
∈ Rn×n1 . Problem ( AARP) is a SDP problem. In case of Chebyshev center, E can
V
be replaced with ρI in all the constraints and maximize ρ in ( AARP), where ρ ∈ R and I is the
identity matrix. Then, the problem becomes a conic quadratic programming (CQP) problem. The
optimal ρ is the approximated radius of the maximum volume inscribed n1 -ball. In the following
example, we apply the new formulation ( AARP) to the polytope in Example 2.1 to determine the
MVE center in the projected polytope.
Example 3.1. The maximum volume inscribed ellipsoid of a full-dimensional polytope with
auxiliary variables. Let us again consider the full-dimensional polytope (2). By solving ( AARP), we
find the optimal MVE center at x = 5 with a radius 5 in Projx (P ). The solution for the adjustable variable
y = 7 − 3ζ for ζ ∈ [−1, 1], which coincides with the line that covers all the feasible values of x in Projx (P ),
i.e., 0.6x + y = 10.
Note that since we restrict y(ζ ) to a linear decision rule instead of a general vector function, the
optimal value of ( AARP) is merely a lower bound on the optimal value of ( AP). This can be
observed from the following two examples.
Example 3.2. Using a linear decision rule results in an approximated MVE. Let us consider the
tetrahedron
P = {( x, y) | x1 + x2 + y ≤ 1, −2x1 − y ≤ 0, x1 − x2 + y ≤ 1, y ≥ 0 } .
(5)
We are interested in the MVE center in the x-space. The approximated MVE of the projection Projx (P )
is a 2-ball with logdet( E) √
= −1.39 centered at x∗ = (0, 0). However, the projection Projx (P ) forms
a square with side length 2. The MVE of the square (i.e. the optimal MVE of Projx (P )) is a disc
with logdet( E) = −0.69 also centered at x∗ = (0, 0). Even though the MVE of Projx (P ) is an underapproximation, its center is the same as the optimal solution in this example.
Example 3.3. Using a linear decision rule results in an approximated MVE center of the projection. Let us consider the tetrahedron
P = {( x, y) | x1 + x2 + y ≤ 1, −1.5x1 − y ≤ 0, x1 − x2 + y ≤ 1, y ≥ 0 } .
(6)
We are interested in the MVE center in the x-space. The approximated MVE of the projection Projx (P )
∗
from ( AARP) is an ellipsoid with logdet(√
E) = −0.92
√ centered at x = (−0.5, 0). However, the projection
Projx (P ) forms a kite with side length 2 and 5. The MVE of the kite (i.e. the optimal MVE of
Projx (P )) is an ellipsoid with logdet( E) = −0.32 centered at x∗ = (−0.26, 0). The approximated MVE
center of Projx (P ) is not equal to the optimal MVE center.
7
Ben-Tal et al. [2009] show that the tractable robust counterpart can be derived for fixed resource
adjustable robust optimization problem with quadratic decision rule under ellipsoidal uncertainty
set. Let the decision rule y(ζ ) be quadratic, i.e.,
y j (ζ ) = ζ T Wj ζ + v Tj· ζ + u j ,
for j = 1, ..., n2 ,
(7)
where u j ∈ R, v j· ∈ Rn1 and Wj ∈ Rn1 ×n1 . The following theorem provides a SDP formulation of
( AP) with a quadratic decision rule.
Theorem 3.1. The solution ( x∗ , u∗ , E∗ , V ∗ , Wj∗ ) is optimal for ( AP) with quadratic decision rule
(7) if and only if it is an optimal solution of the following SDP problem:
( QARP)
max
x,u,τ,λ,E,V,Wj
s.t.
log detE
τ≤b
λ≥0





x
τi − λi − aiT·
u
T
E
− 12
ai ·
V
− 12 aiT·
E
V
λi I − ∑nj=n1 +1 aij Wj−n1



0

∀i.
Proof. See Appendix A. The proof is adapted from Ben-Tal et al. [2009].
In the following example, we apply ( QARP) to approximate the MVE of Projx (P ) in Example 3.2
and 3.3.
Example 3.4. Using quadratic decision rule may result in a better approximated MVE center of
the projection. We use ( QARP) to compute the MVE of the projected tetrahedron Projx (P ) of Example
3.2 and 3.3. In case of Example 3.2, the approximated MVE coincides with the optimal MVE. The approximated MVE for Example 3.3 is an ellipsoid with logdet( E) = −0.35 centered at x∗ = (−0.250, 0), which
is very close to the optimal MVE, i.e., an ellipsoid with logdet( E) = −0.32 centered at x∗ = (−0.257, 0).
The quality of the approximations from ( QARP) is at least as good as the approximations from
( AARP). In §5.1, we provide a more detailed comparison of the performance of ( AARP) and
( QARP). Analogously, in case of Chebyshev center, E can be replaced with ρI in all the constraints
and maximize ρ in ( QARP).
3.2
General Case
Suppose the x-space of the polytope P defined by (1) is not full-dimensional, i.e., there does not
exist a n1 -dimensional ellipsoid in the x-space. Hence, there are (hidden) equality constraints in
polytope P . Let us again assume that x is the main variable (non-adjustable variable) and y is the
auxiliary variable (adjustable variable) in P .
8
We consider two classes of (hidden) equality constraints in P . The first class consists of equality
constraints that contain the auxiliary variable y. The second class consists of equality constraints
that only contain the main variable x. The equality constraints formed by inequalities can be
detected and separated efficiently (see Fukuda [2013]). After separating all the hidden equality
constraints from the inequalities, we can represent the set as follows:
n
o
Projx (P ) = x | ∃y : A1x x + A1y y ≤ b1 , A2x x = b2 , A3x x + A3y y = b3 ,
(8)
where

 1
A1x A1y
b
A =  A2x 0  and b = b2  ,
A3x A3y
b3

where the set without (hidden) equalities { ( x, y) ∈ Rn | A1x x + A1y y ≤ b1 } is full-dimensional.
Firstly, we propose to use constraints that contain auxiliary variable, i.e., A3x x + A3y y = b3 , to
eliminate (part of) the auxiliary variable y (e.g., by Gaussian elimination), until there are no
equality constraints containing y. In Gorissen et al. [2013], the authors briefly discuss two ways
of dealing with uncertain equality constraints with auxiliary variables. Either one can eliminate
the equality constraints with auxiliary variables, or one can apply linear decision rules. Theorem
2 in Appendix B shows that both procedures are equivalent. However, using linear decision
rules requires extra variables. It is therefore preferred to eliminate the auxiliary variables in the
uncertain equality constraints rather than using linear decision rules. In the following example,
the auxiliary variables elimination is applied to a given polyhedron.
Example 3.5. Auxiliary variables elimination by using equality constraints.
Suppose we have the set
Projx (P ) = { x | ∃y ∈ R3 : 2x1 + 3x2 + y1 + y2 + 2y3 ≤ 5, x1 + x2 + y1 = 3,
x2 + y1 = 2.5, y2 + y3 = −2}.
The constraint x1 + x2 + y1 = 3 (or x2 + y1 = 2.5) and y2 + y3 = −2 can be used to eliminate y1 and y2 .
Then, we have a reformulation of the polyhedron set Projx (P ):
Projx (P ) = { x | ∃y3 ∈ R, x1 + 2x2 + y3 ≤ 4, x1 = 0.5}.
There are no equality constraints containing y in the reformulated Projx (P ).
Note that a different choice in the auxiliary variable elimination does not influence the feasible
region of x in Projx (P ). If the elimination procedure results in some equality constraints that
only contain variables x, we can of course add those to the constraints A2x x = b2 . Hence, we can
represent the new (equivalent) set Projx (P ) with no equality constraints containing y as follows:
n
o
Projx (P ) = x | ∃y
e : Dx1 x + Dye1 y
e ≤ c1 , Dx2 x = c2 ,
(9)
where the variable y
e ∈ Rne2 denotes a vector of remaining auxiliary variables after the elimination.
1
1
The matrices Dx , Dye and Dx2 , and the vectors c1 and c2 are the resulting matrices and vectors after
9
the elimination.
Obviously, the equality constraints only in the variable x, i.e., Dx2 x = c2 , reduce the dimension of
the set Projx (P ). For Dx2 6= 0, there does not exist a positive definite matrix E such that,
Dx2 ( x + Eζ ) = c2
∀ζ : ||ζ ||2 ≤ 1.
Hence, the method introduced in §2 cannot be applied directly. We propose to use the following
constraints elimination technique (see Boyd and Vandenberghe [2004]). Let us denote x0 as
a solution of the linear equations Dx2 x = c2 and define the columns of the full rank matrix
F ∈ Rn1 ×(n1 −l ) to be the basis of the null space of Dx2 , where l = rank ( Dx2 ). Then, we have
n
x ∈ Rn1 | Dx2 x = c2
o
n
o
= Fz + x0 | z ∈ Rn1 −l ,
(10)
where the new set (on the right hand side of the equality) is full-dimensional in the variable z.
Hence, the polytope defined by (9) can be reformulated as a full-dimensional set in the variable z:
Pz =
z∈R
n1 − l
| ∃y
e:D
1
Fz + x0
y
e
≤c
1
,
where D1 = ( Dx1 Dye1 ). The MVE center of the full-dimensional set Pz (with auxiliary variables)
can be obtained by solving the following problem:
( MVE − P)
max
z,e
y(ζ ),E
s.t.
log detE
d1i·
T F (z + Eζ ) + x 0
y
e(ζ )
≤ c1i
∀ζ : ||ζ ||2 ≤ 1, ∀i,
where c1i is the i-th element of the vector c1 , and the vector d1i· is the i-th row of the matrix D1 .
Since ( MVE − P) is computationally intractable for a general vector function y
e(ζ ), we restrict the
vector function y
e(ζ ) to a linear decision rule, i.e.,
y
e = u + Vζ.
Then, the MVE center of Pz can be approximated by solving the SDP problem:
( MVE − LP)
max
z,u,V,E
s.t.
log detE
d1i·
T Fz + x 0
u
T
FE
+ d1i· ≤ c1i
V
∀i.
2
Note that the optimal value of ( MVE − LP) is a lower bound on the optimal value of ( MVE − P).
Since the MVE is affine invariant (so is the MVE center), the MVE center of Projx (P ) can be
obtained via an affine transformation of z∗ , i.e.,
x∗ = Fz∗ + x0 .
The Chebyshev center is not affine invariant. A modification of our method is necessary for finding
the Chebyshev center of a given set with auxiliary variables. Let k = dim( Projx (P )) and k ≤ n1 .
10
Table 1: The complexities for computing the Chebyshev center and MVE center.
Classical method
N.A.
LP
SDP
Decision rules
Chebyshev center
MVE center
Our method
Linear Quadratic
CQP
SDP
SDP
SDP
A largest k-ball in the set Projx (P ) defined by (9) can be determined by solving the following
semi-infinite programming problem
(CC − P)
max
x,e
y(ζ ),ρ
s.t.
ρ
d1i·
T x + ρζ y
e(ζ )
≤ c1i
∀ζ : ||ζ ||2 ≤ 1, Dx2 ζ = 0, ∀i
Dx2 x = c2 ,
where y
e(ζ ) is a vector function of ζ, and the extra constraints Dx2 ζ = 0 reduce the dimension of
the ball. Hence, the constraints Dx2 ζ = 0 ensure the ball to be in Projx (P ). We again restrict the
vector function y
e(ζ ) to a linear decision rule, i.e.,
y
e = u + Vζ.
A Chebyshev center of Projx (P ) can be approximated by solving the following CQP problem:
(CC − LP)
max
ρ
s.t.
x,u,ρ,V
d1i·
T
ρI
d1i· ≤ c1i
+ P
u
V
T x ∀i
2
Dx2 x = c2 ,
where I ∈ Rn1 ×n1 is the identity matrix, P = I − ( Dx2 ) T [ Dx2 ( Dx2 ) T ]−1 Dx2 . The matrix P is the
projection matrix that projects vector ζ ∈ Rn1 onto the null space of Dx2 . This technique cannot be
used for computing MVE centers, since it leads to unbounded objective value if the x-space is not
full-dimensional. Note that this method does not reduce the (redundant) dimension of ζ before the
optimization; the dimensions of ζ, u and V in (CC − LP) are higher than the ones in ( MVE − LP).
Both ( MVE − QP) and (CC − QP) (i.e., with quadratic decision rule) can be derived similarly by
using the technique in Appendix A. From Table 1, one can observe that the theoretical complexity
for computing MVE center does not change with respect to different methods.
4
HGK Upper Bounding Method
In this section, we employ the idea in Hadjiyiannis et al. Hadjiyiannis et al. [2011] to compute
upper bounds on the optimal value of ( MVE − P). Hadjiyiannis et al. Hadjiyiannis et al. [2011]
propose to solve adjustable robust optimization problems only with respect to a finite subset of
uncertain parameters, i.e., the critically binding scenarios. The main advantage of their method
is that it does not require to impose any restriction (e.g., linear or quadratic decision rule) on
11
the adjustable variable y. The proposed optimization problem is less restrictive than ( MVE − P).
Hence, it provides an upper bound on the optimal value of ( MVE − P).
In order to obtain a tight upper bound, the subset of scenarios is carefully selected. Hadjiyiannis
et al. Hadjiyiannis et al. [2011] propose to consider the critically binding scenarios (CBS) in Ub =
{ζ 1 , ..., ζ K }, where the CBSs, i.e., ζ 1 , ..., ζ K , are determined by solving the auxiliary optimization
problem:
T FE∗ ζ kL = argmax d1k·
ζ,
(11)
V∗
ζ:||ζ || ≤1
2
where E∗ and V ∗ denote the optimal solution from ( MVE − LP), and k = 1, ..., m. Similarly, in
case of quadratic decision rule, the CBS of ( MVE − P) are defined as follows:
"
#
e2
n1 + n
T FE∗ k
1
T
1
∗
ζ Q = argmax dk·
ζ+ζ
(12)
∑ dkj Wj−n1 ζ,
V∗
ζ:||ζ || ≤1
j = n +1
2
1
where E∗ , V ∗ and W ∗ denote the optimal solution of ( MVE − P) with quadratic decision rule, and
k = 1, ..., m. The CBSs defined by (12) can be obtained by solving a SDP problem (see Appendix
A). In case more than one CBS is determined from the k-th constraint, an interesting research topic
is to determine the CBS that yields the best upper bound. We do not discuss this topic further in
the present paper. Here, we arbitrarily choose one of the CBSs and include it in the CBS set.
The scenario counterpart of ( MVE − P) with respect to the CBS set UbL , where UbL = {ζ 1L , ..., ζ m
L },
is given by the following optimization problem:
( MVE − ubL)
max
z,e
yk ,E
s.t.
log detE
d1i·
T
F (z + Eζ kL ) + x0
y
ek
!
≤ bi1
∀i, ∀k = 1, ..., m.
For the kth scenario ζ kL ∈ UbL , we only need a feasible y
ek to exist. Hence, the auxiliary variables y
ek
are not restricted to any decision rule. Problem ( MVE − ubL) provides an upper bound on the
optimal value of ( MVE − P).
Note that although the HGK upper bounds can be computed in polynomial time, for large
instances, the quadratic growth in the number of constraints of the scenario counterpart can be
problematic. Therefore, we propose the following iterative procedure for solving ( MVE − ubL).
Firstly, we solve the initialization problem
( It − ubL)
max
z,e
yk ,E
s.t.
log detE
d1k·
T
F (z + Eζ kL ) + x0
y
ek
!
≤ bk1
∀k = 1, ..., m.
For each constraint, only its corresponding CBS is considered in ( It − ubL) in the initial iteration.
Then, for k = 1, ..., m, we check if the kth constraint satisfies all the CBSs in UbL with respect to the
optimal z∗ , E∗ and y
ek∗ from ( It − ubL). If there are violating constraints detected, we add those
12
constraints to ( It − ubL) and solve it again. We repeat this procedure until no violating constraints
can be found. The resulting z∗ , E∗ and y
ek∗ from this iterative procedure is an optimal solution for
( MVE − ubL).
Analogously, the upper bounds ( MVE − ubQ) and ( MVE − ubLQ) can be obtained by solving the
b
b
scenario counterpart of ( MVE − P) with respect to UbQ = {ζ 1Q , ..., ζ m
Q } and ( U L ∪ U Q ), respectively.
The iterative procedure only adds the violating constraints in each iteration. It efficiently avoids to
include the redundant constraints. In §5, we applied this procedure to compute the upper bounds
from ( MVE − ubL), ( MVE − ubQ) and ( MVE − ubLQ). The upper bounds are obtained within
5 iterates. In all the cases, this iterative procedure only uses less than 10% of the total constraints
to obtain the upper bounds. More than 90% of the constraints appeared to be redundant. This
procedure significantly improves the computation capability.
The techniques that are discussed in the present section can be easily adapted to compute upper
bounds for (CC − P).
5
Numerical Results
In this section, we evaluate the performance and applicability of our method. In §5.1, we conduct
a simple experiment to examine the behavior of our lower bound method and HGK upper bound
method. In §5.2, our method is applied to find a robust temperature profile for a color picture
tube design problem.
5.1
A Simple Experiment
Let us consider the full-dimensional tetrahedron
Pi = {( x1 , x2 , y) | x1 + x2 + y ≤ 1, x1 − x2 + y ≤ 1, −(1 + i ) x1 − iy ≤ 0, y ≥ 0 } ,
(13)
where x1 , x2 are main variables, y is an auxiliary variable and i = 0, 1, ..., 5. We are interested in
the MVE center and Chebyshev center in the projected Pi onto the x-space. As shown
√ in Example
3.2 and 3.3, when i = 0, the projection Projx (P0 ) forms a square with side length 2; when i > 1,
the Projx (Pi ) forms a kite. In fact, the explicit description of the projected Pi onto the x-space is
as follows:
Projx (Pi ) = {( x1 , x2 ) | x1 + x2 ≤ 1, x1 − x2 ≤ 1, − x1 + ix2 ≤ i, − x1 − ix2 ≤ i } .
We can easily compute the optimal MVE center and Chebyshev center of Projx (Pi ). In this section,
we apply our method to approximate the MVE center and Chebyshev center in Projx (Pi ) for
i = 0, 1, 2, 3, 5. The mean squared deviation (abbr. MSD) of 105 uniformly sampled points in
Projx (Pi ) from the optimal and approximated centers is reported. The MSD is considered as a
measure of the robustness. The same samples are used for generating both Table 2 and Table
3. Note that the geometric center is the most robust solution with respect to the MSD measure,
but it is difficult to compute (see Rademacher [2007]). Hence, we consider the MVE center and
Chebyshev center as our robust solutions.
13
Table 2: The optimal and approximated MVE centers of the projected Pi given in (13) for i = 0, 1, 2, 3, 4, 5.
i
0
1
2
3
4
5
Value
x1
x2
Volume
MSD
x1
x2
Volume
MSD
x1
x2
Volume
MSD
x1
x2
Volume
MSD
x1
x2
Volume
MSD
x1
x2
Volume
MSD
Optimal
0.33
0
-1.65
0.22
0
0
-0.69
0.33
-0.26
0
-0.32
0.56
-0.54
0
-0.06
0.91
-0.83
0
0.13
1.34
-1.14
0
0.29
1.92
Lower bounds
MVE-LP MVE-QP
0.33
0.33
0
0
-1.65
-1.65
0.22
0.22
0
0
0
0
-1.39
-0.69
0.33
0.33
-0.5
-0.25
0
0
-0.92
-0.35
0.58
0.56
-1
-0.58
0
0
-0.55
-0.12
1
0.9
-1.33
-0.92
0
0
-0.26
0.07
1.44
1.32
-1.67
-1.25
0
0
-0.04
0.24
1.99
1.88
14
MVE-ubL
0.33
0
-1.65
0.22
0
0
0
0.33
-0.12
0
-0.11
0.6
-0.33
0
0.03
1
-0.67
0
0.18
1.42
-1
0
0.32
1.99
Upper bounds
MVE-ubQ MVE-ubLQ
0.33
0.33
0
0
-1.65
-1.65
0.22
0.22
0.5
0.09
-0.5
-0.09
0
-0.53
0.84
0.35
0.06
-0.05
-0.01
-0.01
-0.12
-0.13
0.71
0.63
-0.27
-0.33
-0.01
0
0.04
0.03
1.05
1
-0.6
-0.67
0
0
0.19
0.18
1.47
1.42
-0.93
-1
0
0
0.33
0.32
2.04
1.99
In Table 2, the optimal and approximated MVE centers of Projx (Pi ) with their corresponding
Volume are reported. The Volume in Table 2 denotes the objective value logdet( E), which is
proportional to the volume of the ellipsoid. The ( MVE − LP) and ( MVE − QP) denote our (lower
bound) approximation method for MVE centers with respect to linear decision rule and quadratic
decision rule. Suppose UbL and UbQ are the CBS sets obtained from (11) and (12) (see §4). The
( MVE − ubL), ( MVE − ubQ) and ( MVE − ubLQ) are computed from the HGK upper bound
method with respect to different CBS sets, i.e., UbL , UbQ and (UbL ∪ UbQ ), respectively. Table 2 can
be read as follows. For i = 1, the optimal MVE center is at ( x1 , x2 ) = (0, 0) with (proportional)
volume −0.69 and MSD 0.33. The approximated MVE centers from ( MVE − LP) and ( MVE − QP)
are also at (0, 0) and (0, 0) with volume −1.39 and −0.69 (i.e., lower bounds), and MSD 0.33. The
approximated MVE centers from ( MVE − ubL), ( MVE − ubQ) and ( MVE − ubLQ) are at (0, 0),
(0.5, −0.5) and (0.09, −0.09) with volume 0, 0 and −0.53 (i.e., upper bounds) and MSDs 0.33, 0.84
and 0.35, respectively.
From Table 2, one can easily observe that ( MVE − QP) produces the best approximation of
the optimal MVE centers among all the methods. For i = 0, 1, ( MVE − QP) gives the optimal
volume of the MVE. In case of ( MVE − LP), it obtains the optimal volume only when i = 0.
The approximated volumes from ( MVE − QP) are closer to the optimal solutions than those
from ( MVE − LP). The HGK upper bound method gives poorer approximations than those
from ( MVE − QP). The quality of the approximated volumes from ( MVE − ubLQ) is at least as
good as the ones from ( MVE − ubL) and ( MVE − ubQ). For i = 1, 2, the approximated volume
from ( MVE − ubLQ) is strictly better than those from ( MVE − ubL) and ( MVE − ubQ). The
approximated MVE centers from ( MVE − ubLQ) fall in between the approximated centers from
( MVE − ubL) and ( MVE − ubQ).
With respect to the robustness of the centers (i.e., see the MSDs), for most of the cases, our lower
bound method produces more robust centers than the ones from HGK method. The ( MVE − QP)
gives the most robust center with respect to the MSD measure. This is because, for this example,
the centers computed from ( MVE − QP) are the closest to the geometric centers of Projx (Pi ) (i.e.,
the true optimal center with respect to the MSD measure).
In Table 3, the optimal and approximated Chebyshev centers of the projected Pi are presented.
Again, the (CC − QP), i.e., our lower bound method with quadratic decision rule, gives the best
approximations of the optimal Chebyshev centers. The (CC − QP) obtains the optimal radii
for i = 0, ..., 3, whereas, (CC − LP) only produces the optimal radius when i = 0. Our lower
bound method produces more robust centers than the ones from HGK upper bound method. The
(CC − LP) gives the most robust centers with respect to the MSD measure. This is because the
approximated centers from (CC − LP) are the closest to the true optimal centers, i.e., the geometric
centers.
Our new method produces high quality lower bounds, especially, in case of quadratic decision
rule, i.e., ( MVE − QP) and (CC − QP). Using quadratic decision rule (instead of linear decision
rule) can significantly improve the quality of the approximations. The approximated centers
produced by our method are very robust with respect to the MSD measure. From our lower
bounds and the HGK upper bounds, the intervals that contain the optimal volume (or radius)
of the ellipsoids (or balls) can be constructed. This provide us valuable information about the
optimal solutions. Furthermore, the quality of the upper bounds from the HGK method can be
15
Table 3: The optimal and approximated Chebyshev centers of the projected Pi given in (13) for i = 0, 1, 2, 3, 4, 5.
i
0
1
2
3
4
5
Value
x1
x2
Radius
MSD
x1
x2
Radius
MSD
x1
x2
Radius
MSD
x1
x2
Radius
MSD
x1
x2
Radius
MSD
x1
x2
Radius
MSD
Optimal
0.41
0
0.41
0.23
0
0
0.71
0.33
-0.16
0
0.82
0.58
-0.24
0
0.87
1.07
-0.28
0
0.9
1.83
-0.3
0
0.92
2.94
Lower bounds
CC-LP CC-QP
0.41
0.41
0
0
0.41
0.41
0.23
0.23
0
0
0
0
0.5
0.71
0.33
0.33
-0.5
-0.17
0
0
0.62
0.82
0.58
0.58
-0.72
-0.25
0
0
0.72
0.87
0.89
1.06
-0.78
-0.3
0
0
0.78
0.89
1.36
1.79
-0.82
-0.33
0
0
0.82
0.91
2.14
2.88
16
CC-ubL
0.41
0
0.41
0.23
0
0
1
0.33
0.07
0
0.93
0.71
0.04
0
0.96
1.39
0.02
0
0.98
2.35
-0.19
0
0.98
3.7
Upper bounds
CC-ubQ CC-ubLQ
0.41
0.41
0
0
0.59
0.41
0.23
0.23
0.36
0.07
-0.36
-0.07
0.71
0.71
0.6
0.35
0
0.07
0
0
1.22
0.93
0.66
0.71
0
0.04
0
0
1.15
0.96
1.33
1.39
0
0.02
0
0
1.12
0.98
2.3
2.35
0
-0.19
0
0
1.1
0.98
3.66
3.7
improved by considering the CBSs from other decision rules.
Note that the quality of our lower bound approximations obtained by using quadratic decision
rules are at least as good as the ones from linear decision rules. However, this is not the case for
the upper bounds produced by the HGK method. In Table 3, one can observe that it is not clear
whether the upper bounds from (CC − ubL) or the ones from (CC − ubQ) are better. Hence, we
can conclude that the CBSs from UbQ (i.e., by using quadratic decision rules) do not necessarily
yield better upper bounds than the ones from UbL (i.e., by using linear decision rules).
For small i, the MSD of the Chebyshev centers (see Table 3) and MVE centers (see Table 2) are
almost identical. For i = 3, 4, 5, the MVE centers are noticeably more robust than the Chebyshev
centers (as the MSDs are smaller in Table 2). This is because the MVE centers are more centralized
than Chebyshev centers.
5.2
Robust Color Picture Tube Design
In the manufacturing process of a CRT TV, the color picture tube is assembled to the mask-screen
and the inner shield by a industrial oven at a high temperature. The oven temperature causes
thermal stresses on the tube and if the temperature is too high, the pressure will scrap the tube
due to implosion. In den Hertog and Stehouwer [2002], the authors present a temperature profile
of a color picture tube, see Figure 3. The profile is characterized by temperature at 20 locations. To
minimize the production cost and hence the number of scraps, an optimal temperature profile
that satisfies the following criteria:
• the maximal stress for the TV tube is minimal
• the temperature differences between near locations are not too high
• the temperatures are in the specified range.
The authors formulate the associated problem as follows:
( TV − P)
min
x,smax
s.t.
smax
ak + bkT x ≤ smax
∀k ∈ {1, ..., K }
| Ax| ≤ ∆Tmax
l ≤ x ≤ u,
where smax ∈ R is the maximal stress, ak + bkT x ∈ R is the stress at location k of the temperature
profile x ∈ Rn . There are n = 20 temperature points on the tube (see Figure 3). Furthermore,
these temperatures result in K = 216 thermal stresses on different parts of the tube. The 216
response functions of the thermal stresses, ak + bkT x, are derived by using the finite element method
simulator and regression. The vectors l ∈ Rn and u ∈ Rn are the lower and upper bounds of the
temperature profile. The parameter ∆Tmax ∈ Rd represents the maximal allowed temperature on d
location combinations; A ∈ Rd×n is the coefficient matrix that enforce the specified temperatures
do not differ more than ∆Tmax . For example, the temperatures at locations 2 and 5 in Figure 3
cannot differ more than ∆Tmax . By solving the nominal problem ( TV − P), the unique minimum
17
Figure 3: Temperature Profile.
smax = 14.15 is found (see Table 4).
Suppose the maximum tolerance of smax is at most 15. We are interested to find a temperature
profile that is robust against implementation error and not far from the optimal profile x¯ of
( TV − P), say, || x − x¯ ||1 ≤ 15. The feasible set for the robust temperature profile is as follows:
n
o
P = x ∈ Rn | ak + bkT x ≤ 15, | Ax| ≤ ∆Tmax , l ≤ x ≤ u, || x − x¯ ||1 ≤ 15, ∀k .
(14)
Due to the extra constraints, the feasible set P is much smaller than the feasible region of ( TV − P).
Our version of robust solution is the MVE center of the feasible set. As shown in Example 2.2,
the set P can be easily rewritten into a polytopic representation with some auxiliary variables.
However, the explicit description of the projected P onto the x-space requires exponentially many
constraints. We apply our method to find the robust profile, i.e., the approximated MVE centers of
the projected P onto the x-space.
In Table 4, the profiles that are obtained from ( TV − P) (the nominal profile), our method and
HGK method are presented, respectively. The smax denotes the maximal stress of the profiles. The
Volume denotes the objective value logdet( E), which is proportional to the volume of the ellipsoid.
The mean squared deviation (MSD) of the profiles is reported in Table 4. The MSDs are computed
from 105 uniformly sampled profiles in the set (14). One can observe that the nominal profile has
a slightly lower maximal stress smax , i.e., less than 1%, but its MSD is more than 30% higher than
the other profiles. With less than 1% compromise on the smax , one can increase the robustness of a
profile by more than 30%. The profiles produced by using our method are more robust than other
profiles. The profile obtained from ( MVE − LP) is more robust than the one from ( MVE − QP).
In Table 5, the computation time of all the methods are given. All the procedures are performed by
using SDPT3 (see Toh et al. [1999]) within Matlab R2012a on an Intel Core i5 CPU running at 2.9
GHz with 4 GB RAM under Windows 7 operating system. One can observe that the HGK upper
bound method takes much more time than our method. Although the theoretical complexity of
( MVE − LP) and ( MVE − QP) are equivalent, the required computation time for ( MVE − LP) is
18
Table 4: The robust color picture tube design: the approximated MVE centers.
Value
x1
x2
x3
x4
x5
x6
x7
x8
x9
x10
x11
x12
x13
x14
x15
x16
x17
x18
x19
x20
smax
Volume
MSD
Nominal
27.02
16.47
26.68
-20
27.17
25.76
16.86
20.66
-4.02
50
50
50
16.12
42.39
50
50
50
50
-20
-20
14.15
N.A.
19.16
Lower bounds
MVE-LP MVE-QP
27.02
27.02
16.47
16.45
26.68
26.63
-19.29
-19.26
27.17
27.17
25.76
25.80
16.86
16.89
20.66
20.68
-4.02
-4.02
49.29
49.26
49.29
49.24
49.29
49.25
16.12
16.13
42.39
42.39
49.29
49.26
49.29
49.26
49.29
49.26
49.29
49.26
-19.29
-19.26
-19.29
-19.26
14.24
14.24
-6.09
5.35
14.24
14.28
MVE-ubL
27.09
16.30
26.49
-19.19
27.22
25.93
16.90
20.82
-4.08
49.20
49.19
49.22
16.03
42.26
49.24
49.24
49.23
49.22
-19.22
-19.22
14.24
10.92
14.53
Upper bounds
MVE-ubQ MVE-ubLQ
27
26.99
16.51
16.5
26.55
26.54
-19.19
-19.19
27.14
27.17
25.03
26.03
16.92
16.94
20.57
20.55
-3.98
-4
49.15
49.15
49.16
49.17
49.16
49.16
16.09
16.11
42.3
42.31
49.18
49.19
49.17
49.18
49.17
49.18
49.17
49.18
-19.16
-19.17
-19.16
-19.17
14.25
14.26
6.54
6.51
14.52
14.51
Table 5: The robust color picture tube design: the computation time (in seconds) of different methods.
Computation time
MVE-LP
14.6
MVE-QP
244.9
19
MVE-ubL
792.9
MVE-ubQ
1331.8
MVE-ubLQ
3881.9
much less. In case of HGK method, the more CBSs are considered, the more time is needed for
computing the corresponding upper bound.
6
Conclusions and Further Research
In this paper, we have developed a new method to approximate two centers, i.e., the MVE center
and Chebyshev center, of a projected polytope. The main advantage of our method is that it does
not require an explicit description of the projection. Our method is computationally tractable. The
numerical results suggest that our method provides high quality approximations.
Besides finding robust designs, our method can be applied for many other purposes. In the
following, we list some possible applications.
Application 1. Ellipsoidal approximations to polytopic sets. For polytopes with many redundant
constraints and/or facets, Zhang and Gao [2003] argue that ellipsoids are much easier to handle, both
theoretically and computationally. For instance, the global minimum of any quadratic in an ellipsoid can be
located in polynomial time while finding such global minimum in a polytope is generally a NP-hard problem.
Since the description of the polytopic feasible set for power system problems often contains prohibitively
many constraints, Sari´c and Stankovi´c [2008] propose to use the MVE to approximate the polytopic sets.
Our new method allows auxiliary variables in the polytope which enhances the modeling flexibilities of the
proposed approach in Sari´c and Stankovi´c [2008].
Application 2. Nominal scenario recovery in polytopic uncertainty set. In robust optimization,
the uncertainty set that contains all possible scenarios is given. Suppose the uncertainty set is a polytope.
One may be interested to find a nominal scenario that is not far from the (later) actual realization. For
example, the approximated nominal scenario can be used to evaluate the price of robustness. However, there
may exist auxiliary variables in the description of the uncertainty set. The MVE center or Chebyshev center
of the projected set can be considered as approximations of the nominal scenario.
Application 3. Cutting-plane method. The goal of cutting-plane methods is to find a point in a convex
set X . The general procedure of cutting-plane methods is as follows. Suppose a polytope P that contains the
set X is provided, i.e., X ⊆ P . We first select a point x(0) in P . If x(0) ∈ X , then we are done; otherwise,
a separating hyperplane between x(0) and X can be constructed. This hyperplane cuts off the halfspace that
contains no points of X in P . Let us denote the resulting set P after the first cut as P1 . Then, we select
another point x(1) in P1 and repeat the procedures until we locate a point in X . A centralized point x(k) of
Pk is preferred as it leads to a deeper cut of Pk . The location of the selected point x(k) ∈ Pk in each iteration
determines the convergence rate of the method. Many choices of x(k) ∈ Pk are proposed in the literature
(see Boyd and Vandenberghe [2007]), e.g., the geometric center, the Chebyshev center, the MVE center, the
analytic center. Our new method can be used as an extension of the MVE center cutting-plane method and
Chebyshev center cutting-plane method. Suppose the set X is described by two kinds of variables, i.e., the
variable y appears linearly in all the constraints and the variable x appears nonlinearly in some constraints.
The complexity of locating a point in X arises by the variable x, since if all the variables appear linearly in
the description of the set X , finding a point in it can be done by solving a LP problem. Therefore, we focus
on cutting the projection of P where the variable x resides in. The cutting-plane can be determined by the
Chebyshev center or MVE center of the projected polytope (i.e., the given polytope that is projected onto the
20
space where the nonlinear variables reside in). Our method may significantly speed up the convergence rate
of the cutting-plane method in those cases.
One of the possible extensions of our method is that instead of only considering (centers of)
ellipsoids and balls, one can also consider (centers of) boxes. This can be done by replacing
the 2-norm function of ζ with 1-norm or ∞-norm in, for example, ( AP). Another interesting
research topic is to find robust solutions for system of uncertain linear equations. Given a system
of uncertain linear equations Ax = b, where the coefficient matrix A and right-hand side vector b
reside in a uncertainty set U. Our method can be applied to find the MVE center and Chebyshev
center of the solution set (i.e., { x | ∃( A, b) ∈ U : Ax = b}) as the robust solutions. The extension
of our method to find robust solutions for system of uncertain linear equations will be further
analyzed in future research.
Appendices
A
Tractable Robust Counterpart of ( AP) with Quadratic Decision Rule
Proof. This proof is adapted from Ben-Tal et al. [2009]. The quadratic decision rule is defined by
(7):
y j (ζ ) = ζ T Wj ζ + v Tj· ζ + u j
for j = 1, ..., n2 ,
where u j ∈ R, v j· ∈ Rn1 and Wj ∈ Rn1 ×n1 . Let us consider the constraint i in (AP), i.e.,
aiT·
x + Eζ
y(ζ )
≤ bi
∀ζ : ||ζ ||2 ≤ 1.
(15)
We have
"
#
n
x
T
T E
max
+ ai ·
ζ+ζ
∑ aij Wj−n1 ζ
u
V
ζ:||ζ ||2 ≤1
j = n1 +1
| {z } | {z }
|
{z
}
aiT·
pi
= max
ζ:||ζ ||2 ≤1
n
= min τi
τi
n
= min τi
τi
n
= min τi
τi
= min τi
τi
pi + 2qiT ζ
2qiT
Ri
T
+ ζ Ri ζ
| ∀(ζ : ||ζ ||2 ≤ 1) : τi ≥ pi + 2qiT ζ + ζ T Ri ζ
o
| ∀((ζ, t) : ||ζ ||2 ≤ t2 ) : (τi − pi )t2 − 2tqiT ζ − ζ T Ri ζ ≥ 0
o
o
| ∀(ζ, t), ∃λi ≥ 0 : (τi − pi )t2 − 2tqiT ζ − ζ T Ri ζ − λi (t2 − ζ T ζ ) ≥ 0 [S − Lemma]
−qiT
τi − pi − λi
0, λi ≥ 0 .
|
− qi
λi I − Ri
21
Hence, the tractable robust counterpart of (15) can be formulated as follows

τi ≤ bi




λi ≥ 0




 
x
1 T E
T
− 2 ai ·

 τi − λi − ai· u
V




T
 
 0.




E


− 21
ai ·
λi I − ∑nj=n1 +1 aij Wj−n1

V
B
Uncertain Equality Constraints with Auxiliary Variables
Let us consider the following uncertain linear equality constraint:
a(ζ )T x + y = b
∀ζ ∈ U ,
(16)
where a(ζ ) = a + Pζ is linear in ζ ∈ Rm , a ∈ Rn and P ∈ Rn×m . The variable y ∈ R is an
auxiliary variable (without loss of generality we assume the coefficient for y is 1) and U denotes
the uncertainty set.
Theorem B.1. Assume that 0 ∈ U . Then using a linear decision rule for the auxiliary variable y in
(16) is equivalent to eliminating y.
Proof. Let the linear decision rule be as follows:
y = u + v T ζ,
where u ∈ R and v ∈ Rm . Then we have
( a + Pζ )T x + u + v T ζ = b
T
T
T
⇔ a x + u + ( P x + v) ζ = b
∀ζ ∈ U
∀ζ ∈ U .
The equality constraint is satisfied for all ζ in U if and only if
T
a x+u = b
( P T x + v)T ζ = 0 ∀ζ ∈ U .
With u = b − a T x and v T ζ = − x T Pζ, we have,
y = u + v T ζ = b − a T x − x T Pζ = b − a(ζ ) T x.
Hence, substituting y = b − a(ζ ) T x everywhere in the problem is equivalent to using a linear
decision rule for y.
22
References
A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski. Adjustable robust solutions of
uncertain linear programs. Mathematical Programming, 99-2:351–376, 2004.
A. Ben-Tal, L. El Ghaoui, and A. Nemirovski. Robust Optimization. Princeton Series in Applied
Mathematics. Princeton University Press, Princeton, NJ, 2009.
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, Cambridge, UK,
2004.
S. Boyd and L. Vandenberghe. Localization and Cutting-Plane Methods. From Stanford EE 364b
lecture notes, 2007.
F.S. Chaharsooghi, M.J. Emadi, M. Zamanighomi, and M.R. Aref. A new method for variable
elimination in systems of inequations. Proceedings IEEE International Symposium on Information
Theory (ISIT), :1215–1219, 2011.
L. Danzer, D. Laugwitz, and H. Lenz. Über das Löwnersche Ellipsoid und sein Analogon unter
den einem Eikörper einbeschriebenen Ellipsoiden. Archiv der Mathematik, 8:214–219, 1957.
D den Hertog and HP Stehouwer. Optimizing color picture tubes by high-cost nonlinear programming. European Journal of Operational Research, 140(2):121–197, 2002.
J.B.J. Fourier. reported in: Analyse des travaux de l’Académie Royale des Sciences, pendant l’année
1824, Partie mathématique (1827). Histoire de l’Academie Royale des Sciences de l’Institut de France,
7:47–55, 1824.
K. Fukuda. Polyhedral Computation. Lecture notes, ETH Zurich, 2013.
B. Gorissen, I. Yanikoglu, and D. den Hertog. Hints for practical robust optimization. Technical
report, CentER Discussion Paper No. 2013-065, 2013.
H.E. Graeb. Analog Design Centering and Sizing. Springer-Verlag, 2007.
M.J. Hadjiyiannis, P.J. Goulart, and D. Kuhn. A scenario approach for estimating the suboptimality
of linear decision rules in two-stage robust optimization. Proceedings IEEE Conference on Decision
and Control and European Control Conference (CDC-ECC), :7386–7391, 2011.
E.M.T. Hendrix, C.J. Mecking, and T.H.B. Hendriks. Finding robust solutions for product design
problems. European Journal of Operational Research, 92:28–36, 1996.
T. Huynh, C. Lassez, and J.-L. Lassez. Practical issues on the projection of polyhedral sets. Annals
of Mathematics and Artificial Intelligence, 6:295–316, 1992.
C.N. Jones, E.C. Kerrigan, and J.M. Maciejowski. Equality set projection: A new algorithm for
the projection of polytopes in halfspace representation. CUED Technical Report CUED/FINFENG/TR.463, 2004.
C. Lassez and J.-L. Lassez. Symbolic and Numerical Computation for Artificial Intelligence, chapter 4,
pages 103–119. Academic Press, 1992.
J.-L. Lassez. Query constraints. Proceedings ACM Conference Principles of Database Systems, 4:288–298,
1990.
23
T.S. Motzkin. Beiträge zur Theorie der linearen Ungleichungen. University Basel Dissertation.
Jerusalem, Israel, 1936.
L.A. Rademacher. Approximating the centroid is hard. In Proceedings of the Twenty-third Annual
Symposium on Computational Geometry, SCG ’07, pages 302–305, New York, NY, USA, 2007. ACM.
ISBN 978-1-59593-705-6. doi: 10.1145/1247069.1247123. URL http://doi.acm.org/10.1145/
1247069.1247123.
A.T. Sari´c and A.M. Stankovi´c. Applications of ellipsoidal approximations to polyhedral sets in
power system optimization. IEEE Transactions on Power Systems, 23-3:956–965, 2008.
H.R. Tiwary. On computing the shadows and slices of polytopes. CoRR, abs/0804.4150, 2008. URL
http://arxiv.org/abs/0804.4150.
K.C. Toh, M.J. Todd, and R.H. Tütüncü. Sdpt3 – a matlab software package for semidefinite
programming. Optimization Methods and Software, 11:545–581, 1999.
V.L. Zaguskin. Circumscribed and inscribed ellipsoids of extremal volume. Uspehi Matematicheskih
Nauk, 13:89–93, 1958.
Y. Zhang and L. Gao. An interior-point algorithm for the maximum-volume ellipsoid problem.
SIAM Journal on Optimization, 14:53–76, 2003.
24