Global optimization approach to unequal sphere packing problems in 3D.
The problem of the unequal sphere packing in a 3-dimensional polytope is analyzed. Given a set of unequal spheres and a polytope, the double goal is to assemble the spheres in such a way that (i) they do not overlap with each other and (ii) the sum of the volumes of the spheres packed in the polytope is maximized. This optimization has an application in automated radiosurgical treatment planning and can be formulated as a nonconvex optimization problem with quadratic constraints and a linear objective function. On the basis of the special structures associated with this problem, we propose a variety of algorithms which improve markedly the existing simplicial branch-and-bound algorithm for the general nonconvex quadratic program. Further, heuristic algorithms are incorporated to strengthen the efficiency of the algorithm. The computational study demonstrates that the proposed algorithm can obtain successfully the optimization up to a limiting size.
Introduction
The optimization of the packing of unequal spheres in a 3-dimensional polytope is ana-
lyzed. Given a set of unequal spheres and a polytope, the objective is to assemble them
in such a way that (1) the spheres do not overlap with each other and (2) the sum of the
volumes of the packed spheres is maximized. We note that the conventional 2-dimensional
and 3-dimensional packing problems (also called the bin-packing problem), which have
been extensively studied [2, 7, 11], are fundamentally di#erent from the problem considered
in the present work.
The unequal sphere packing problem has important applications in automated radio-
surgical treatment planning [12, 14]. Stereotactic radiosurgery is an advanced medical
technology for treating brain and sinus tumors. It uses the Gamma knife to deliver a set
of extremely high dose ionizing radiations, called "shots", to the target tumor area [13]. In
good approximation, these shots can be considered as solid spheres. For large or irregular
target regions, multiple shots are used to cover di#erent parts of the tumor. However,
this procedure usually results in (1) large dose inhomogeneities, due to the overlap of the
di#erent shots, and (2) the delivery of large amount of dose to normal tissue arising from
enlargement of the treated region when two or more shots overlap.
Optimizing the number, the position, and the individual sizes of the shots can significantly
reduce both the inhomogeneities and the dose to normal tissue while simultaneously
achieving the required coverage. Unfortunately, since the treatment planning process is
tedious, the quality of the protocol depends heavily on the experience of the users. There-
fore, an automated planning process is desired. To achieve this goal, Wang and Wu et
al. [12, 14] mathematically formulated this planning problem as the packing of spheres
into a 3D region with a packing density greater than a certain given level. This packing
problem was proved to be NP-complete and an approximate algorithm was proposed [12].
In this work we formulate the question as a nonconvex quadratic optimization and present
solution methods based on the branch-and-bound technique.
Let K be the number of di#erent radii of the spheres in the given set and r k
1, ., K) the corresponding radii. There are L available spheres for each radius. Therefore,
the total number of the spheres in the set is KL. Here, we use a single value of L for
simplicity of the presentation. However, this model can be easily modified for the case
where di#erent numbers L k of spheres are available for di#erent radii r k and the total
number of spheres in the given set is
Let the polytope be given by { (x, y, z) # R 3
L be the maximum number of the spheres to be packed.
We designate variables as the location of sphere i in a packing. For
each sphere in the packing, a radius has to be assigned. The variables t ik
1, ., K) are used to handle this task:
sphere i has radius r k ,
With these preliminaries, the optimization can be formulated as follows.
.3
s.t.
# a 2
Constraints (1) and (3) respectively ensure that no two spheres overlap with each other
and that each sphere is centered within the polytope. Constraints (4) and (5) guarantee
that at most one radius is chosen for each sphere, i.e., if t then the sphere i is
packed with radius r k , and if t then the sphere i with radius r k is not packed.
Together with (3), (4) and (5), constraints (2) state that the distance between the center
of a sphere and the boundary of the polytope is at least as large as the radius of that
sphere. Hence, (2) and (3) force all the spheres to be packed inside the polytope.
. By (3), constraint (2) can be rewritten as
for each i, m. (6)
Since the right-hand side of (6) is nonnegative, (6) implies (3). Moreover, the binary 0-1
variables t ik in (5) can be replaced by the inequalities t ik (t ik -
Note that t ik # 1 is implied by the constraint (4). These steps allow restatement of
Problem (P1) as follows.
Note that the constant 4# in the original objective function is omitted here. The numbers
of the variables and the quadratic constraints are (3+K)L and L(L-1)+LK, respectively.
The quadratic function in each constraint (7) is neither convex nor concave.
There exist several algorithms [1, 9] that have been developed for solving the general
nonconvex quadratic program (NQP in short). A NQP can be transformed into a
semidefinite programming problem (SDP) with an additional rank-one constraint. Dropping
the rank-one constraint, one obtains a SDP relaxation problem, which is the tightest
relaxation among all others [3]. Hence, (P2) can be solved by using a branch-and-bound
type method based on the SDP relaxation. However, due to the relative slowness and the
instability of the SDP software, this conceptual algorithm is impractical even for problem
having a small size, since it requires many iterations of the SDP relaxation.
Commonly, methods of solution for the NQP class are designed through linear programming
(LP) relaxation, an approach known as the reformulation-linearization technique
[10]. Under the assumption that additional box constraints for the variables are
new nonlinear constraints are generated from each pair of the box constraints in
order to construct a relaxed convex area of the original nonconvex feasible region. The
nonlinear terms in each nonlinear constraint are accordingly replaced by new variables.
The result of this recasting is a linear program. Based on this technique of linearization,
Al-khayyal et al. [1] proposed a rectangular branch-and-bound algorithm to solve a class
of quadratically constrained programs. Raber [9] proposed another branch-and-bound
algorithm for the same problem based on the use of simplices as the partition elements
and the use of an underestimate a#ne function, whose value at each vertex of the simplex
agrees with that of the corresponding nonconvex quadratic function. This work [9] demonstrated
that the simplicial algorithm often has a better performance over the rectangular
algorithm with respect to the computational time. Interestingly, Raber [9] mentioned
that both the simplicial and the rectangular algorithms exhibit poor performance for a
packing problem :
giving the experimental details. It is obvious that their packing problem is similar to
ours, but has much simpler structures for the constraints.
In this paper, we examine a customization of Raber's simplicial branch-and-bound
algorithm [9] tailored to our problem. It is well known that the underlying simplicial sub-division
is a key factor influencing the quality of the relaxation, therefore, the e#ciency of
the branch-and-bound algorithm. The investigation of the structure of the optimization
suggests (i) an e#cient simplicial subdivision and (ii) di#erent underestimations of the
nonconvex functions. Based on these observations, three variants of the algorithm have
been constructed. The discrete nature of the packing enables a heuristic design for obtaining
good feasible solutions, an outcome which leads to savings on both computational
time and memory size.
The remainder of this paper is organized as follows. In Section 2 we derive the LP
relaxation of the problem with respect to the simplicial subdivision. The simplicial branch-
and-bound algorithm is presented in Section 3. Section 4 gives the heuristic algorithms
and Section 5 presents three variations of the previous algorithm based on the use of
special structures. Section 6 reports the computational results of the proposed branch-
and-bound algorithm. The conclusions of the work are presented in Section 7.
Programming Relaxation
The construction of the LP relaxation of Problem (P2) is the same as the one developed
in [9]. For the sake of a complete description, we outline this procedure below.
the transpose of a vector a. First, we write Problem (P2) in the following form.
d
d ik # R n , and c # R n is the coe#cient vector of the objective
d
correspond to the constraints (7) and (8),
respectively. Av # b represents all linear constraints of (9), (10) and (11). Furthermore,
the matrix Q il can be specified as follows.
il O
O
where O is a matrix having zero for all entries with appropriate size,
.
.
. 1 . -1
. 1 . -1 .
corresponds to the coe#cients of in (7) for i and l, and
. r 2
. r 1 r K r 2 r K . r 2
. r 2
. r 1 r K r 2 r K . r 2
K .
corresponds to those of (t 11 , ., t 1K , ., t L1 , ., t LK ) T in (7) for i and k. Similarly,
R 3L-3L and
d ik # R KL-KL can be written as follows:
O
. O. 0
. O
d
m be the number of the linear constraints in (P2). Denote the polytope defined by
these linear constraints as
m-n and b # R -
.
To construct the LP relaxation problem, we need to represent the matrix Q il (resp.
by the sum of a positive semidefinite matrix C il (resp.
negative semidefinite
matrix D il (resp.
# D ik ). Usually, a spectrum decomposition achieves this goal.
However, we do not need to perform such a task, since the matrices Q il and
possess
special structures that give the decomposition immediately. It is readily seen that the
decompositions
O
il O
O O
and
satisfy the desired property.
Now we consider how to construct our linear programming relaxation problem. Let
be an n-simplex (U #S #), where v i are its vertices. Then
Let W S # R n-n be a matrix which consists of columns (v
each point v in S can be represented as
Through substitution of v in (P2') by (14), we obtain the following equivalent problem.
s.t.
d
By replacing Q il and
with (12) and (13), the quadratic term of the left-hand side
of each constraint of (15) and (16) is divided into a convex and a concave functions
by replacing Q il and
with (12) and (13), respectively. The relaxation of the above
problem is constructed by ignoring the convex part and replacing the concave part with a
linear underestimate function. For such an underestimation, we use the convex envelope
of a concave function f with respect to the simplex S, which is an a#ne function whose
value at each vertex of S coincides with that of f . More precisely, for the quadratic
constraints (15) we have
where
(v
is the convex envelope of (W S #) T D il W S # with respect to S. In a similar fashion, for the
quadratic constraints (17), we have
d
d
where
(v
is the convex envelope of (W S #) T
# D ik W S # with respect to S.
Obviously, an upper bound of the objective function of (P3) can be obtained by solving
the following LP relaxation problem.
d
3 The Simplicial Branch-and-Bound Algorithm
The simplicial branch-and-bound algorithm is presented in this section. As mentioned
above, we use the algorithm in [9] as our prototype. However, two heuristics designed for
obtaining feasible solutions are embedded.
The branching operation is carried out by dividing the current simplex S into two
simplices. Let v i # and v j # be two vertices of S satisfying
where #a# 2 denotes the 2-norm of a vector a. Define
The simplex S is split into two simplices
and
The splitting of the simplices has the property that for each nested sequence {S q } of
simplices,
For further details, see Horst [4, 5]. The
resulting algorithm is presented below.
Branch and Bound Algorithm :
Step 1. Start Heuristic-1 to calculate a possible feasible solution v f and the objective
function value f(v f ). If successful, set
Step 1.1. Let Construct a simplex S 0 which contains the polytope U . Set
Solve the problem (LPR) S0
. If the problem is infeasible, then the
original problem has no solution, stop. Otherwise, let the optimal solution be
and the optimal value be -(S 0 ). Set is a feasible solution of
(P2), stop. Otherwise, start Heuristic-2 with v 0 to calculate a feasible solution v # 0
and the value f(v # 0 ). If LB < f(v # 0 ), then set
Step 1.2. If (UB - LB)/UB #, then stop.
k and S 2
k according to (25) and (28). Set
k }. For
Step 2.1. Solve the problem (LPR) S j
. If it is infeasible, set
let the optimal solution be
k and the optimal value be -(S j
Step 2.2. If v j
k is a feasible solution of (P2), set
Step 2.3. If v j
k is not feasible, then run Heuristic-2 with v j
k to calculate a feasible solution
(v j
and the value f((v j
Otherwise, select a simplex -
-(S). Set S
2.
The details of Heuristic-1 and Heuristic-2 will be given in Section 4. The parameter # gives
the tolerance of the solution obtained by the algorithm. We call the solution obtained
from the above algorithm an #-optimal solution. The convergence of the algorithm is
guaranteed as follows.
Theorem 1 ([6, 8, 9]) If the algorithm generates an infinite sequence {v k }, then every
accumulation point v # of this sequence is an #-optimal solution of Problem (1).
The algorithms Heuristic-1 and Heuristic-2 are described in this section. Recall that
Heuristic-1 finds a feasible solution at the beginning of the branch-and-bound algorithm
and Heuristic-2 generates a feasible solution from an infeasible solution obtained from a
relaxation subproblem.
First, we give the algorithm of Heuristic-1. The basic idea is to place as many spheres
as possible having relatively large radii in the polytope. Let C be a list of the given
set of KL candidate spheres, which are ordered such that r 1 # r KL . Consider a
3-D triangle defined by four constraints arbitrarily chosen from those of the polytope P .
Designate
triangles. Fixing a P i , the algorithm starts
by picking a sphere from the top of C. Then it checks whether the sphere can be located
at one of the four corners of the triangle P i . This procedure is continued for the rest of
the spheres on the list until four spheres are placed or the search of the spheres in the list
is exhausted. All the spheres packed have to satisfy (1)that they touch exactly three sides
of the triangle, (2) that no mutual intersection occurs between each pair of the spheres
(see

Figure

1), and (3) that they satisfy all other constrains on P . After obtaining the
initial packing, the algorithm attempts to insert more spheres between each pair of the
spheres in P i without violating the packing condition (see Figure 2).
Phase I

Figure

1: Figure 2:
Step 1. Determine the center of the sphere l so that it is tangent to three planes corresponding
to the three constraints of P s . If the sphere l satisfies the constraints of
polytope P and does not overlap with the spheres 1, ., k - 1, then pack the sphere
l and set
Step 2. If (i)k = L or (ii) l = KL and k > 0, go to Step 5.
Step 3. If stop. Otherwise set l=1 and
go to Step 1.
Step 4. Set l to Step 1.
Phase II
Step 5. If Otherwise, for each pair of the packed spheres i, j, set l = 1,
repeat Steps 6 and 7.
Step 6. If l > KL stop. Otherwise set
Step 7. Locate the center of the sphere l at M . If it satisfies the constraints P and does
not overlap with the spheres 1, ., k, then pack the sphere l and set
go to Step 6.
By the termination of the algorithm if k > 0 a feasible packing is obtained.
Next, we describe Heuristic-2. Let be the solution
of a relaxation subproblem. Suppose that it is not feasible. Then either t # ik is not integral
or some spheres overlap with each other if the radius is decided by t # each i,
i.e., use r k as the radius of the sphere i (see Figure 3). In the following, we fix the
centers of these spheres and determine their radii so that they do not mutually overlap
(see

Figure

4). Note that for a triplet
means
that no sphere is placed there. Let r 1 > - > r K .

Figure

3: Figure 4:
Heuristic-2
Step 1. Set
Step 2. If the sphere l centered at satisfies the polytope constraints
and does not overlap with spheres 1, ., l - 1, then set t #
(the sphere l centered at
l , z #
l ) has radius r k ), go to Step 4.
Step 3. If k < K, set 2. Otherwise, set t # lk
(no sphere is centered at
l , z #
l )).
Step 4. If l < L, set l 2. Otherwise, stop.
5 The Improvement of the Algorithm
In this section we discuss the special structures of the optimization and present the improvements
on the previous algorithm based on these structures. First, we describe how
to construct the initial simplex.
5.1 Generating the Initial Simplex
To start the branch-and-bound algorithm, we need an initial simplex S 0 which contains
the polytope U . It is generated by the following algorithm.
Generate Initial Simplex:
Step 1. Choose a nondegenerate vertex v 0 of the polytope U .
Step 2. Construct a matrix A # R n-n from the n binding constraints
Step 3. Compute
Step 4. Set S
Remark 1: The n constraints which decide the nondegenerate vertex v 0 can be selected
as follows. For each i we choose three constraints from (9) and all constraints of (11). It
is not di#cult to see that n coe#cient vectors from these constraints are
linearly independent. Set these linear inequality constraints as linear equations, i.e.,
The above system of linear equations determines an unique solution which can be considered
as v 0 .
Remark 2: The vertex v 0 has KL zero entries, which are determined by (31). Further-
more, if an equation in (30) is replaced by (-
#, then the solution obtained
maintains t k. Repeating this process for all equations in (30) yields
3L solutions, which we index as v 1 , ., v 3L . All the other solutions generated by replacing
are denoted by v 3L+1 , ., v n . More precisely, we have
and
5.2 Splitting a Simplex in a Di#erent Way
It is well known that in the simplicial branch-and-bound paradigm both the computational
time and the usage of the memory grow extremely fast as the dimension of the polytope
U increases.
When a simplex S is split into two simplices according to (25) - (28), the length of
is the longest among all other pairs of the vertices of the simplex S. If the
vertex v 0 is not replaced by the vertex v M , then the matrix W S 1 is the same as W S except
the column v i # - v 0 , which is replaced by the column v M - v 0 . Similarly, all columns
in the matrix W S 2 are the same as W S except the column v j # - v 0 , which is replaced
by v M - v 0 . Recall that the entries (i, are zeros in the matrices
(12). Hence, the coe#cients of # j in # Sil (#) of (22),
only depends on the first 3L coordinates of v 0 ,
values for the first 3L coordinates, then # Sil (#) remains unchanged in the constraints
(22) of the subproblem with respect to the simplices S j (j = 1, 2), and the quality of the
relaxation would not be significantly improved. Such a splitting leads to the computation
of subproblems which provide neither good lower nor useful upper bounds. To avoid this
situation, we choose v i # , v j # such that
(v
is the maximum among all other Consequently, the coordinates corresponding
to t ik will not be considered.
One negative e#ect of this choice is that the quality of the relaxation of the constraint
necessarily improved. However, since the linear constraint
is always considered, this defect is not expected to be serious. Accordingly, the number
of subproblems can be decreased in the whole branch-and-bound process.
There is another potential di#culty with the convergence of the algorithm which must
be addressed. Let -
(v
are vertices of S q }. Since we
divide the simplex to minimize the largest value
(v vertices, it is
possible that a nested sequence {S q } of simplices with -
does not satisfy (29). In this case, we are not guaranteed that the accumulation point of
an infinite sequence obtained by the algorithm will be an optimal solution.
5.3 Another Decomposition of the Matrix Q il
As shown in Section 2, to make the underestimation of the quadratic function of the
left-hand side of the quadratic constraint (7), we use the decomposition of (12). From
Replacing t 2
ik by t ik in (7) could result in a di#erent matrix
il O
O
.
. r 2
. r 1 r k r 2 r k . r 2
and
A spectral decomposition of the matrix Q # il provides positive semidefinite and negative
semidefinite matrices which are di#erent from those given in (12). Following a procedure
similar to that in Section 2, we have
where C # il and D # il are positive semidefinite and negative semidefinite, respectively. For
example, if only two magnitudes of radii r 1 and r 2 are considered, then C # t
il and D # t
R KL-KL can be constructed as follows. Let
2 . Then
O
il O
O
where
. r 2(r 2+3r 2+s)
. r 2(r 2+3r 2+s)
. r 1 r 2 r 2
. r 2
. r 2
. r 1 r 2 r 2
and
il .
It can be shown that the eigenvalues of matrix C # t
il are nonnegative (resp.
nonpositive). Hence, the two matrices have the desired properties. Consequently, another
LP relaxation can be constructed based on the new spectral decomposition of
that
remains the same here, so there is no change in its decomposition.
5.4 Another Form of the Relaxation
In this section we focus on a di#erent form of relaxation of problem (P2). Let us omit
the quadratic constraints (8), which are corresponding to the 0 - 1 condition of t ik .
Furthermore, we relax the condition that two spheres can not overlap with each other.
Here the distance between two spheres only needs to be greater than a small positive
value -#. More precisely, consider the problem
s.t.
Remark. One can take the smallest value among all radii given in the set as the magnitude
of - #.
Since only the variables are quadratic in (P4), the convex
envelope of the concave function -(x i - x l
can be determined
by the first 3L coordinates of the vertices of the corresponding simplex as given above.
This implies that it is su#cient to construct simplices in the 3L-dimensional space. Let
be the initial simplex in the 3L-dimensional space that contains the
polytope P . The ith entry of v # j is identical to that of v j for
# be defined
similarly as Q il and A, respectively. Let A # x and A # t be the submatrices of A # , which are
corresponding to
we have the following quadratic program.
s.t.
negative semidefinite matrix, the convex envelope of the quadratic term
(v
Therefore, we obtain the relaxation of (P4) as follows.
Since the constraints (8) are ignored and the constraints (7) are relaxed as (42) in the
above problem, the quality of this relaxation may be inferior to the previous methods.
However, the dimension of the simplices kept in the memory is only 3L. Therefore, the
total memory used may be far smaller than that required in the other methods.
6 The Computational Study
In this section we discuss details of the implementation of the algorithm and report the
experimental results.
Reoptimization : Suppose that S is the simplex at a branching node. The simplex S
is divided into S 1 and S 2 by (25) to (28). If the vertex v 0 is not replaced by v M , then
only columns of matrices W S 1 and W S 2 , are replaced by v M - v 0 ,
respectively. All the other columns remain unchanged. Therefore, in this case, when
two new subproblems are generated by splitting S into S 1 and S 2 , we could use the
information of the optimal solution of the LP relaxation problem associated with S. By
the construction of the LP relaxation problem with respect to S 1 and S 2 , only those
coe#cients corresponding to # i # and # j # in the constraints (22) and (23) are changed due
to the substitutions of the column v M - v 0 in W S 1 and W S 2, respectively. Consequently,
by changing these coe#cients, we could use the reoptimization technique to solve the LP
relaxation problem starting from the optimal solution corresponding to the simplex S.
Test Problems : The test problems are generated as follows. First, construct a simplex
with vertices (0, 0, 10), (10, 0, 0), (0, 10, 0) and (10, 10, 10). Calculate the maximum
inscribed sphere of the simplex (the radius is 2.88675). Then randomly generate
m points on the surface of the sphere. Construct tangent planes {(x, y, z) | a
which pass each of those m points respectively. Let
y, z) | a be the halfspace containing the inscribed sphere.
The intersection of H i with the simplex is the 3-dimensional polytope in which the sphere
packing problem is considered. Hence, the total number of the linear constraints of the
polytope is In our test, m was set at 4. Di#erent pairs of radii of spheres
that we used are shown in Tables 1-2.

Table

1: Combinations of radii for
C r radius
6 1.00 0.50

Table

2: Combinations of radii for
C r radius
The computational experiments were conducted on a DEC Alpha 21164 Workstation
(600MHz). We used CPLEX 6.5.1 as an LP solver for the relaxation problems. The
limits of the memory and the computational time were set at 512MB and 3600 seconds,
respectively. Besides the prototype algorithm given in Section 3, three variations were
also implemented. They were (1) Algorithm NSS, which uses the simplex splitting given
in Section 5.2, (2) Algorithm NMD, which uses the matrix decomposition given in Section
5.3, and (3) Algorithm XYZ, which uses the simplex in the XYZ space given in Section
5.4. For each algorithm, we tested 5 instances for each and each pair C r of radii
from the Table 1. Since the other three algorithms reached either the limit of memory
or the limiting computational time for most of the instances for present
the results of Algorithm XYZ. The data shown in Tables 4-7 and Table 10 represent the
average of the results. The legends used in the tables are given in Table 3.
From

Tables

4-7 we observe that the computational time grows drastically as the
number of spheres packed is increased, since the dimension of the problem is (3
This outcome is consistent with the observation in [1, 9] that the computational time

Table

3: Legends used in the tables
legend meaning
L the maximum number of the spheres packed
C r the combination of the radii
#.LP the number of linear programs solved
#.T the number of cases terminated caused by the limit of CPU
time (3600 seconds)
#.E the number of cases terminated when
#.M the number of cases terminated caused by the limit of memory (512 MB)
Time CPU time (seconds)
UB-LB
the ratio of the values (UB - LB) and UB, where UB and LB is the
upper and lower bounds of the objective function value, respectively
Algorithm ORG the prototype algorithm given in Section 3
Algorithm NSS the algorithm using the simplex splitting given in Section 5.2
Algorithm NMD the algorithm using the spectral decomposition given in Section 5.3
Algorithm XYZ the algorithm using the simplex in the XYZ space in Section 5.4
increases exponentially as the dimension increases. There is no big di#erence between
Algorithm ORG and Algorithm NMD, in the sense of their gross behaviors. Algorithm
NSS solves fewer LPs than the two previously mentioned algorithms do, therefore, it needs
less time. This leads to the solution of a greater number of instances without violating
the limits on either memory or time. Since the simplex splitting is based on the length of
the first 3L coordinates, it avoids solving unnecessary subproblems.
Among all the algorithms, Algorithm XYZ shows the best performance. All instances
for are solved except the case C It even successfully obtained
the solutions for about half of the instances for 5. Two reasons for this behavior
are considered. First, since the dimension of simplices is only 3L, much less memory is
needed to keep the information on the simplices. Second, since the simplices involve only
the coordinates of variables (x, y, z), the simplex splitting has the favorable characteristic
of Algorithm NSS that fewer LPs are needed to be solved.
It should be noted that (i) the dimension of the instances solved in the previous
papers [1, 9] range up to 16 and (ii) the computational time and memory demand of the
simplicial branch-and-bound algorithm usually increases exponentially with the dimension
of the problem. Therefore, a modest increase in the dimension could put the solution out
of range. In contrast, for the sphere packing problem discussed above, the dimension of
the instances solved by Algorithm XYZ extended up to
Hence, we conclude that Algorithm XYZ gains e#ciency, since it takes advantage of the
intrinsic structure of the problem.
Next, let us focus on Algorithm ORG and Algorithm XYZ and consider the influence
of the polytope P on their performance. Five instances of the polytope are generated
randomly and they have same number of constraints. Taking the results are
shown in Tables 8-9. It is observed that even if L and C r are identical, the comparative
behaviors of the algorithms are very sensitive to the shape of the polytope. This behavior
arises since (i) the quality of solutions generated by the heuristics depends on the shape
of the polytope and (ii) the initial simplex S 0 depends on the polytope. For a skinny
polytope, the volume of S 0 \ P 0 would be large, a fact which means that the algorithms
copiously waste e#ort on solving LPs defined on this zone before reaching the optimal
solution.

Table

shows the results when three values of radii are considered, i.e.,
that (i) the value of K e#ects the number of the quadratic constraints in (8), and (ii)
increasing the number of the quadratic constraints enhances the di#culty of the problem.
Comparing this case with the results for (Table 7), we observed that there is no
big change in the numbers of LPs solved. The algorithm is less sensitive to the value of K
than to the size of of L. This occurs, since the numbers of both the quadratic constraints
(7) and (8) grow with increasing of L.
Finally, all of the algorithms exhibited reduced performance when C namely,
when the di#erences of the radii are large.
Conclusions
In this paper, we considered the optimization of unequal sphere packing and demonstrated
the improvements over the existing simplicial branch-and-bound algorithm through advantageous
use of the intrinsic structure of the problem. Specially, the computational
study showed that the improved Algorithm XYZ could solve instances with much larger
size. We observed that optimal solutions were found for many instances when the algorithm
reached the limitations of either computational time or memory. This signals that

Table

4: Results of Algorithm ORG
Time #.E #.T #.M
the algorithm spends a large amount of e#ort verifying the optimality of the solution.
In turn, this behavior indicates that (1) developing an improved method of relaxation is
necessary for solving problems with a larger size, (2) when the algorithm is terminated,
the solution obtained then can be of high quality, and (3) a better heuristic method, which
could start with a feasible solution obtained from the branch-and-bound algorithm, would
be very important for obtaining the approximate solution to larger problems.

Table

5: Results of Algorithm NMD
Time #. E #. T #.M

Table

Results of Algorithm NSS

Table

7: Results of Algorithm XYZ
Time #.E #.T #.M

Table

8: The results of Algorithm ORG
Time

Table

9: The results of Algorithm XYZ
Time

Table

10: The results of Algorithm XYZ
Time #.E #.T #.M



--R

A relaxation method for nonconvex quadratically constrained quadratic programs

Semidefinite programming relaxations for nonconvex quadratic programs
On generalized bisection of n-simplices
Handbook of Global Optimization
Global Optimization - Deterministic Approaches
SIAM Journal on Comput- ing
Global Optimization in Action
A simplicial branch-and-bound method for solving nonconvex all- quadratic programs
A new reformulation-linearization technique for bilinear programming problems
A strip-packing algorithm with absolute performance bound 2
Packing of unequal spheres and automated radiosurgical treatment plan- ning
Physics and dosimetry of the gamma knife

--TR
On three-dimensional packing
On generalized bisection of <italic>n</italic>-simplices
A Strip-Packing Algorithm with Absolute Performance Bound 2
A Simplicial Branch-and-Bound Method for Solving Nonconvex All-Quadratic Programs
LP relaxation;nonconvex quadratic programming;simplicial branch-and-bound algorithm;heuristic algorithms;unequal sphere packing problem