p TIMA
Mathematical Programming Society Newsletter
MAY 2001
,. I pi
.. ..... ......
)* it#4 
i F
Ni
 1;
~~J
I .
u'' h
i J':
MA72001
GEORGIA ON MY MIND....
The Seventeenth Symposium on Mathematical
Programming A Great Success
Karen Aardal
The Symposium was held at the Georgia
Institute of Technology, August 711, 2000. The
main organizers were George Nemhauser,
Donna Llwellyn, and Martin Savelsbergh. The
atmosphere of the meeting was as warm as the
outdoor climate, and the quality of the technical
as well as the social program was supreme. There
were 1,055 registered participants at the meeting.
Seven plenary talks were given, including a
special lecture by Thomas Hales about his proof
of the Kepler Conjecture, and the Banquet pres
entation given by Harold Kuhn. The plenary
talks were focusing on the history of some of the
key areas of Mathematical Programming, where
as the ten semiplenary talks were considering
more recent developments. About 870 technical
talks were given in parallel sessions.
On Sunday evening there was a welcome
reception and registration at the Ferst Center for
the Arts the "heart" of the Symposium site.
The Symposium was opened on Monday with
welcome addresses given by: George Nemhauser;
the provost, Michael Thomas; and the MPS
chair, JeanPhilippe Vial. During the opening
session, the Beale OrchardHays, Dantzig,
Fulkerson, and Tucker prizes were awarded as
well. The winners and the prize citations can be
found on the following pages. Rita Graham con
tributed with a musical program that was clearly
appreciated by the mathematical programmers. I
have never seen so many people swinging in
their chairs at a conference! On Wednesday
there was a fantastic conference banquet. Harold
Kuhn gave an entertaining, interesting, and very
moving Banquet Presentation that will long remain
in our memories.
For the first time in the MPS history,
Founders Awards were given to some of our
prominent colleagues who have laid the founda
tion for the many areas of mathematical pro
gramming. The ten recipients were: William
Davidon, George Dantzig, Lester Ford, David
Gale, Ralph Gomory, Alan Hoffman, Harold
Kuhn, Harry Markowitz, Philip Wolfe, and
Guus Zoutendijk. Lester Ford and David Gale
could unfortunately not attend the meeting.
EIGHT OF THE FOUNDERS AWARDS RECIPIENTS, FRC
LEFT TO RIGHT: PHILIP WOLFE, HAROLD KUHN, HAF
MARKOWITZ, RALPH GOMORY, GEORGE DANTZIG, A
HOFFMAN, GUUS ZOUTENDIJK, AND WILLIAM DAVII
SOUTHERN HOSPITALITY! THE
HAPPY ORGANIZING TEAM:
GEORGE NEMHAUSER, DONNA
LLWELLYN, AND MARTIN
SAVELSBERGH.
WHICH SESSION TO GO TO
NEXT? GEORGE DANTZIG,
RALPH GOMORY, AND ELLIS
JOHNSON ARE CHECKING THE
PROGRAM.
MARTIN GROETSCHEL GIVING A PLENARY TALK
AT THE OPENING SESSION.
HAROLD KUHN DURING HIS BANQUET
PRESENTATION "IN THE RIGHT PLACE AND THE
RIGHT TIME."
PAST, PRESENT, AND FUTURE CHAIRMEN OF THE SOCIETY ENJOYING THE BANQUET GEORGE
NEMHAUSER, CHAIRELECT BOB BIXBY, JOHN DENNIS, GEORGE DANTZIG, JAN KAREL LENSTRA, AND
CURRENT CHAIR JEANPHILIPPE VIAL.
S 10 PA
PAGE 2
MAY 2001
THE MPS PRIZES
Listed below are the prize citat
MPS prizes (The Beale Orch
the Dantzig prize, the Fulkers
Tucker prize). Some of the cit
slightly edited by the OPTIM.
The Beale OrchardF
The Beale OrchardHays
lence in computational mathel
ming was awarded to David
Bixby, Vasek Chvital, and Wil
their paper, "On the Solution
Salesman Problems," Docume
Journal der Deutschen Mathel
Vereinigung, ICM III (1998),
Additional information and
the associated CONCORDE
are available online (http://ww
.edu/concorde.html also give
itself).
DAVID APPLEGATE
VASEK CHVATAL
The Traveling Salesman Pro
played a special role in the hist
cal programming. In addition
cations, from manufacturing t
computational biology, the TS
as the paradigmatic hard comb
ions for the four
mization problem and as a prime first testbed for
ardHays prize, new ideas in the field, from cutting planes to
on prize, and the branchandcut to simulated annealing and
nations have been other metaheuristics.
A Editor. The award winning paper gives an informa
tive summary of the history of cutting plane
ays Prize based approaches to the TSP and then describes
the ideas, both old and new, that the authors
prize for excel have incorporated in their own code, which is
matical program now by far the best optimization code for the
Applegate, Robert TSP that has been developed. In its 12 pages
liam Cook for (constrained by a page limitation for this special
of Traveling issue of the journal) the paper manages to
nta Mathematica include a large amount of material at a level such
matiker that the nonexpert can get a feel for the ideas
645656. involved and the expert can understand the key
source code for innovations involved in the code.
software package As with all experimental papers, however, it is
w.keck.caam.rice the code itself and the results obtained with it
;n in the paper that provide the measure of significance. The
authors' solution of a 13,509city instance from
the TSPLIB is arguably the most difficult opti
mization problem that has been solved to date.
The 13,509 solution was the culmination of 10
Years of work by the group, pushing the enve
lope of optimization computing in a wide range
of areas. The procedure that the authors
employed is a model for how to combine mathe
matical programming techniques with advanced
data structures and algorithms, together with
parallel computing, to greatly extend the reach
ROBERT BIXBY
ROBERT BIXBY of optimization software.
The authors have advanced the state of the art
to the point where solving 1,000city instances is
routine. Indeed, in a recent test by Applegate (to
gather information on the BeardwoodHalton
Hammersley constant) the authors' code was
used to solve 10,000 geometric instances each
containing 1,000 points.
The "local cuts" procedure described in the
paper (and implemented via a tour de force in
their computer code) is an important advance in
WILLIAM COOK .
the way cuttingplane separation algorithms can
be viewed. Previously our main route to produc
blem (TSP) has
ing facetdefining inequalities (so useful in prac
tory of mathemati
tice over the past twenty years) was a template
to its many appli
to its many appli based approach that required a great deal of
o scheduling to
scheduling to problemspecific knowledge, analysis, and algo
P has long served .
P has long served rithms. The new and complementary approach
inatorial opti
(and its natural extension to general integer pro
gramming) provides a general purpose routine
that can be specialized to produce facets auto
matically by a process of linear projection, given
only that one has a method for optimizing small
instances of the relevant subproblem.
Byproducts of the authors' TSP work have
found their way into the leading commercial
integer and linear programming codes. A promi
nent example of this is the strong branching rule
(for choosing a branching variable in mixed inte
ger programming branch and bound algorithms)
that is now the default choice when solving hard
MIP instances.
The authors' source code includes a little over
100,000 lines of C code. This is (to my knowl
edge) the first time the full source code to a
stateoftheart branchandcut algorithm has
been made publicly available, and it serves as a
model for the development of future efficient
implementations. Included with the cutting
plane routines are the first publiclyavailable
highquality implementations of the Lin
Kernighan and ChainedLinKernighan heuris
tics for finding nearoptimal tours, which them
selves are of great practical benefit when
attempting to get good solutions to instances
beyond the range of the optimization routines.
All this allows future research to start at the top,
rather than having to rebuild a code from the
ground up.
The two versions of the accompanying source
code (one from 1997 and one from 1999) have
had over 1,200 downloads (and the new version
continues to obtain about 40 new downloads
per week). Links to it have already been added
to over 25 web sites (outside of AT&T, Rice,
and Rutgers). This amount of activity is unusual
(to say the least) for a piece of advanced mathe
matical software, and a definite measure of its
impact.
An honorable mention for the Beale 
OrchardHays prize was given to the following
two papers: Joseph Czyzyk, Michael P Mesnier,
Jorge Mord, The NEOS Server, IEEE Journal on
Computational Science and Engineering, 5
(1998), 6875.
William Gropp, Jorge More, Optimization
environments and the NEOS Server, in:
SMP IM 65
PAGE 3
MAY 2001
Approximation Theory and Optimization, M.
D. Buhmann and A. Iserles, eds., pages 167
182, Cambridge University Press, 1997.
The prize jury: M. Ferris, B. Fourer, D.
Goldfarb (chair), M. Kojima, M. Saunders.
The George B. Dantzig Prize
The Dantzig prize is jointly administered by
the Mathematical Programming Society and the
Society for Industrial and Applied Mathematics
(SIAM). The prize is awarded to one or more
individuals for original research which, by virtue
of its originality, breadth and depth, is having a
major impact on the field of mathematical pro
gramming.
The prize was awarded to Yurii Nesterov.
Yurii Nesterov is probably best known for his
fundamental contributions to the theory of inte
rior point methods for convex programming,
including his introduction of the concept of
selfconcordance. However, his previous works
in convex optimization and in automatic differ
entiation, though perhaps less widely publicized,
were just as creative and authoritative. These
accomplishments have justifiably attracted inter
national recognition, and have significantly
affected the course of research in optimization
theory and algorithms.
The prize jury: W Cunningham, C.
Lemar&chal, S. M. Robinson (chair), L.A.
Wolsey.
The D. Ray Fulkerson Prize
The Fulkerson Prize for outstanding papers in
the area of discrete mathematics is sponsored
jointly by the Mathematical Programming
Society and the American Mathematical Society.
The jury awarded two prizes at the Symposium.
One award went to Michele Conforti, Gerard
Cornuejols and M. R. Rao for their paper
"Decomposition of balanced matrices," Journal
Combinatorial Theory, Series B 77 (1999), 292
406.
A balanced matrix is a 01 matrix with no
submatrix of odd order with precisely two l's
per row and per column. This simply defined
class of matrices has been of great interest and is
wellstudied because packing and covering
Linear Programs with balanced constraint matri
ces have integral solutions. Also, balanced matri
ces are a natural generalization of bipartite
graphs.
The authors give a polynomial time algorithm
for recognizing balanced matrices, settling a
longstanding open question. The algorithm is
based on a deep, beautiful decomposition theo
rem whose arduous proof combines brilliant
insights and meticulous casechecking.
A second award went to Michel Goemans
and David Williamson for their paper
"Improved approximation algorithms for the
maximumcut and satisfiability problems using
semidefinite programming," Journal of the
Association for Computing Machinery, 42
(1995), 11151145.
This paper gives a much better approximation
algorithm for the fundamental maximumcut
problem in graphs and several other optimiza
tion problems by an elegant application of semi
definite programming. It is the first paper to
apply semidefinite programming in an approxi
mation algorithm with a guaranteed perform
ance bound.
A crucial new element is the extraction of an
integer optimal solution from the semidefinite
relaxation. The authors' techniques have since
found a host of other applications.
The prize jury: W. Cook, R. Graham, R.
Kannan (chair).
M.R. RAO, MICHELE CONFORTI, GERARD
CORNUEJOLS, AND RAVI KANNAN (JURY CHAIR)
PAGE 4
DAVID WILLIAMSON, MICHEL GOEMANS AND RAVI
KANNAN
The A.W. Tucker Prize
The Tucker prize is awarded at each ISMP
meeting for the best paper authored by a student.
The committee chose three finalists for the
prize, and from among these a winner. The first
two finalists are Fabian Chudak from Cornell
University, and Kamal Jain from Georgia Tech.
Dr. Chudaks paper considers improved approxi
mation algorithms for the uncapacitated facility
location problem. Dr. Jain's paper describes a
factor2 approximation algorithm for the gener
alized Steiner network problem.
The third finalist, and winner of the Tucker
prize, is Bertrand Guenin from CarnegieMellon
University. Dr. Guenin's paper gives a complete
characterization of graphs for which a natural
LPrelaxation of the MAXCUT problem has all
integral vertices. Such graphs were termed
"weakly bipartite" by Grotschel and Pulleyblank.
Guenin's result verifies a substantial special case
of a conjecture by Seymour on ideal clutters dat
ing back to 1977, and represents the first major
progress on the resolution of Seymour's conjec
ture. The paper employs a sophisticated graph
theoretic analysis involving the use of Lehman's
results on minimally nonideal matrices, for
which Lehman was awarded the MPS Fulkerson
prize in 1991. In addition to its purely theoreti
cal interest, the study of ideal matrices has sig
nificant applicability in the design of algorithms
for partitioning and covering problems that arise
in many areas, such as vehicle routing and crew
scheduling.
The prize jury: K. Anstreicher (chair), R.
Burkhard, D. den Hertog, D. Karger, J. Lee.
10PTIMA65
MAY 2001
Solution of a weighing
problem
Cor A.J. Hurkens*
Gerhard J. Woegingert
PAGE 5
1 The seven weights problem
The following problem was published in the
Sunday Times on Sunday, 30 June 1963, and a
prize of three pounds was offered for the first
correct solution:
"On our last expedition to the interior," said
the famous explorer, "we came across a tribe
who had an interesting kind of Harvest Festival,
in which every male member of the tribe had to
contribute a levy of grain into the communal
tribal store. Their unit of weight was roughly
the same as our pound avoirdupois, and each
tribesman had to contribute one pound of grain
for every year of his age. The contributions were
weighed on the tribe's ceremonial scales, using a
set of seven ceremonial weights. Each of these
weighed an integral number of pounds, and it
was an essential part of the ritual that not more
than three of them should be used for each
weighing, though they need not all be in the
same pan.
We were told that if ever a tribesman lived to
such an age that his contribution could no
longer be weighed by using three weights only,
the levy of grain would terminate for ever. And
in the previous year, one old man had died only
a few months short of attaining the critical age,
greatly to the relief of the headman of the tribe.
The scientist with our expedition confirmed
that the weights had been selected so that the
critical age was the maximum possible."
What was the age of the old man when he
died? And what were the weights of the seven
ceremonial weights?
In a more mathematical formulation, the
problem asks for seven positive integers
W, ... W, that can represent consecutive
integers in the widest possible range 1 . L.
Here an integer n is said to be represented by the
numbers W .... W, if there exists an index
1 < a < 7 with n = W or if there exist two
indices a, b with 1 < a < b < 7 such that
n = + W + WY, or if there exist three indices
a, b, c with 1 < a < b < c < 7 such that
n = +W+ W, W.
In the solution published in the Sunday
Times of 7 July 1963, it was claimed that the
best possible range was for L =120, that the
weights of the seven ceremonial weights were 1,
3, 7, 12, 42, 75, and 100, and that the old man
had died just before his 121st birthday. In this
note we show that this claim is not justified:
The truth is that the old man had died just
before his 123rd birthday. In Section 2, we
describe a weight sequence that represents the
range 1... 122. In Section 3, we argue that a
range 1... 123 cannot be reached; the argu
ment relies on complete enumeration with the
help of a computer program. This program is
available on request from the first author (by
sending email to wscor@win.tue.nl).
2 Solutions to the problem
When we started our search, we used a primi
tive program that only searched through a very
limited number of feasible solutions. Never
theless, within two hours of running time this
program managed to find a solution for the
range 1... 122. Table 1 demonstrates that all
integers in the range 1... 122 indeed can be
represented by the seven weights 1, 3, 7, 12, 43,
76, and 102.
We then completely settled the problem by
performing an implicit complete enumeration
with the help of a computer program. For this
we needed a more sophisticated approach,
which is described in detail in Section 3. In
Table 2, we list in lexicographically increasing
order all possible sets of seven weights that can
represent ranges 1 . L with L > 120.
3 Implicit enumeration
We have implemented an enumeration scheme
to find the best possible sequence of weights W,
where W, W2< W3 W< W, < W6 W.
Our enumeration scheme searches through the
weight sequences in lexicographically increasing
order. A naive approach that simply enumerates
all possible sequences will fail since the number
of combinations to be inspected is fairly large:
Even if we could assume a priori that W, < 123
(and there is absolutely no reason for assuming
this), there still would be roughly 101 sequences
to consider; and for each single one of these
1011 sequences, we had to check whether it rep
resents the full range 1 .. 122 of integers.
*wscor@win.tue.nl; Department of Mathamatics and Computing Science, Eindhoven University of
Technology, P.O. Box 513, NL5600 MB Eindhoven, The Netherlands.
tgwoegi@opt.math.tugraz.ac.at; Institut fiir Mathematik B, TU Graz, Steyrergasse 30, A8010 Graz,
Austria. Supported by the START program Y43MAT of the Austrian Ministry of Science.
0PTIMA65
MAY 2001
1 = 1
2 1 +3
3 3
4 1+3
5 7 + 1:
6 =1 +7
8 1+7
9 3 + 12
10 3+7
11 1 + 12
12 12
13 1 + 12
14 1 + 3 + 12
15 3 +12
16 1 +3 +12
17 43 +76102
18 1 + 7 + 12
19 7 +12
20 1 +7 + 12
21 12 43 + 76
22 3+7+12
23 3 76 + 102
24 7 12 + 43
25 1 76 + 102
76 + 102
1 76 + 102
312+43
3 76 + 102
1 12 + 43
12 + 43
1 12+ 43
43 + 76
1 43 + 76
17+43
7 + 43
17+43
712+ 43
1 3 + 43
3 + 43
13+43
1 + 43
43
1 +43
1+3+43
3 +43
1+3+43
7 + 12 + 43
1 + 7 + 43
7 +43
51 + 7 + 43
52 3 + 12 + 43
53= 3 + 7 + 43
54 1 + 12 + 43
55 12 + 43
56 + 12 + 43
57= 7 12 + 76
58 1 43 + 102
59= 43 + 102
60 1 43 + 102
61 3 12 + 76
62= 3 43 + 102
63= 1 12 + 76
64= 12 + 76
65 = 12 + 76
66 3 7 + 76
67= 3 12 + 76
68  7 + 76
69= 7 + 76
70 1 7 + 76
71= 7 12 + 76
72  3 +76
73 3 +76
74= 1 3 + 76
75= 1 + 76
76 76
77 1 + 76
78= 1 + 3 + 76
79 3 + 76
80= 1+ 3 +76
81 7 + 12 + 76
82= 1 + 7 + 76
83= 7 + 76
84 + 7 + 76
85 3 + 12 + 76
86= 3 + 7 + 76
87= 1 + 12 + 76
88 12 + 76
89 = + 12 + 76
90= 12 + 102
91 1 12 + 102
92 37 +102
93 3 12 + 102
94= 1 7 +102
95= 7 + 102
96= 1 7 + 102
97 7 12 + 102
98 1 3 + 102
99 3 + 102
100 = 1 3 + 102
101= 1 + 102
102 102
103 = 1 102
104 1 3 + 102
105= 3 +102
106 = 1 3 + 102
107 7 + 12 + 102
108 1 +7 +102
109= 7 +102
110 =1 +7 + 102
111 3 +12 + 102
112 3 +7 + 102
113 1 + 12+ 102
114 12+ 102
115 + 12 + 102
116 3 +43 +76
117 3 + 12 + 102
118 1 +43 +76
119 43 + 76
120 = + 43 + 76
121 7 + 12 + 102
122 3 +43 + 76
PAGE 6
W, W W, W4 W, W, W, Range
1 2 12 17 33 74 99 1...120
1 3 7 12 42 75 100 1...120
1 3 7 12 43 76 102 1...122
1 3 9 14 40 70 103 1...120
1 3 10 15 38 72 103 1...121
2 4 7 22 35 78 90 1...121
2 5 6 15 43 74 103 1...120
2 6 13 23 33 83 84 1...120
3 5 15 16 33 75 99 1...120
Table 2: All possible solutions for ranges
1 ... L with L> 120
+4(73 ) /* combinations abs(W W
W) */
01
if 3W > /*compensating for too high
otherwise triples W +Wk+W */
Table 1: How to represent all the integers in the range 1. .. 122
by the seven weights 1, 3, 7, 12, 43, 76, and 102.
Hence, we implemented several tricks that
helped to considerably cut down the search tree.
These tricks investigate potential extensions of
partial solutions. We say that a partial solution
W, < .. < W (with 1 < k < 6) is interesting, if
there is still some possibility to extend it by
numbers W+,, ..., W, to a sequence that repre
sents all the numbers in the range 1 ... 0. The
value 0 is a threshold that will be set to some
value > 117; the choice of this threshold is
somewhat arbitrary, but the analysis below is
general. In order to find good bounds on the
number of distinct representable values in the
range 1 ... 0, we compute and maintain the fol
lowing sets and numbers for every partial solu
tion W, < ... W,.
X < 6 is the number of currently fixed
weights.
M(k, K) with 1< k < 3 is the set of all
integers that can be represented by using at
most k numbers from W, ... W. This set
includes 0 and negative numbers.
W is a (tentative) lower bound on the
numbers W,,..., W.
From now on we will write #(X) to denote the
cardinality of a set X. Note that the cardinalities
of the sets IN n M(1, X), IN n M(2, X), and
IN n M(3, X) can be bounded from above by the
values K, K2, and k2+ 4(3), respectively. Similarly,
the cardinalities of the sets M(1, X), M(2, X), and
M(3, X) can be bounded from above by 2X + 1,
2X2 + 1 and 2X2 + i i + 1, respectively. We will
use these bounds many times in our arguments.
X 1 2 3 4 5 6 7
2 1 4 9 16 25 36 49
2+4(Q) 1 4 13 32 65 116 189
3.1 Locally bounding the depth of
the search tree
The following expression A(L, W 0) is an
upper bound on the number of distinct positive
integers in 1 ... 0 that can be represented from
any extension of the current partial solution. For
ease of exposition, any number that can be cho
sen from the set M(., .) is denoted by A, and
any weights yet to be fixed are denoted by
W< W,_< W with < k< 1.
A(, W, 0)=
#( [1, 0] n M(3, )) ) /* combinations
without any W */
+( 7) #( ( e W] n M(2, X))
/* combinations abs(W f ) */
+(72 ) #(M(1, X)) /* combinations
abs(W W, ) */
+( ) # ( ( 0 2W ] M(1, X)) /*
combinations abs(W + W,+ M) */
2(73') if W > /*compensating
0 otherwise for W+W +W1 and
W WW+W "
E.g., the fifth additive term 4( ;") is justified
by the observation that at most 4 of the 8 com
binations +W Wk W, can be positive. Note
that the function A(X, W 0) is nonincreasing
in parameter W and that W > W. Hence, if
A(L, W, 0) < 0 holds, then the considered par
tial solution cannot be interesting. In this case,
we skip all extensions of this partial solution.
3.2 Locally bounding the width of
the search tree
If a partial solution survives the above feasibility
test and satisfies A(X, W, 0) > 0, then our goal is
to compute lower and upper bounds on W,,
the value of the next weight to be fixed. Since
we generate the weight sequences in lexicograph
ically increasing order, the value W yields a triv
ial lower bound for W +. Now let us look for
upper bounds. Intuitively it is clear that if we set
the weight W, too high, then we will not be
able to represent all the 'small' integers below
W +. To make this idea precise, we compute a
crude upper bound B(L, 0) on the number of
integers in the range 1... min {0, W+ 1} that
can be represented by any extension of the
considered partial solution:
B(L, 0)=
#( [1, 0] n M(3, )) ) /* combinations
without any W */
[S P TIMA Li5
I I I I
MAY 2001
+( ) #( [1, o) n M(2, X))
/* combinations W M */
+(72 ) #(M (1, X)) /* combinations
abs (W W+ M) *
+(7 3) /* combinations abs(W + Wk W) */
(1) First we discuss the case where B(X, 0) < 0
holds. Here we use the upper bound W+, <
B(L, 0) + 1. Note that we can bound B(L, 0)
from above by 2 + 4() + ( ) 2 2+(7 )
(21+ 1) + ( 3). Here L2 + 4(3) is an upper
bound on #( [1, 0] n (3, )) ), and (7 ') L2 is
an upper bound on (7 1) #( [1, oo) n M(2,
X)) and so on. This yields the following a priori
upper bounds on W+,:
X 0 1 2 3 4 5 6
B(X, 0)< 56 72 84 95 108 126 152
W,,+ 57 73 85 96 109
In fact if B(L, 0) < 0, we may compute an even
better upper bound on W, by excluding the
very high values in M(3, L). Define
x* = max{x I x < 0 and B(,, x) > x}.
Then the inequality W,,+ B(L, x* + 1) + 1
yields the improved upper bound on W,,.
(2) Next we discuss the remaining case where
B(L, 0) > 0. Here we define
y* = min{y A(L, y, 0) < 0}
By the definition of function A in the preceding
section, we get that for W >y* there is no feasi
ble extension of the considered partial solution.
Hence, in this second case we arrive at the upper
bound W,+
It remains to be shown that the above defined
value y* indeed exists, i.e., that the minimum in
the right hand side exists (which is not at all
obvious). By examining the function A, one
observes that A(L, , 0) is bounded from above
by L2 + 4(3) + (2L + 1) (72,) + (731): The right
hand side of the definition is the sum of seven
terms. The first term #( [1, 0]n M(3, M)) is
bounded from above by L2 + 4( ). The second
and the fourth term disappear for sufficiently
large values of W The third term
( 2 ) #M(1, X)) is bounded from above by
(/ ) (2X + 1). The sum of the last three terms
equals (7 ) for sufficiently large W Altogether,
this yields the following table.
X 1 2 3 4 5 6 7
A(0, 0,)< 66 64 59 60 76 116 189
Since < 6 and since 0 > 116, the value y*
indeed exists. Finally, we remark that the value
y* can be computed quickly by bisection search.
To summarize, for every partial solution we
first determine the value B( X, 0). Depending on
whether this value is < 0 or> t, we compute an
upper bound that is based on x* ory*, respec
tively. Then we only search extensions of the
partial solution for which W, 1 takes values
between W and this upper bound.
3.3 Fathoming a node by looking
ahead
Starting from an empty sequence of weights, we
recursively compute and extend partial solutions.
The bounds described above are used to limit
the number of combinations that have to be
considered. It is computationally cheap to evalu
ate the functions A and B, since we compute
and maintain all necessary auxiliary values and
sets on the fly.
Once five weights have been fixed (i.e., once
we are down to X = 5), it is possible to efficiently
determine whether the partial solution is inter
esting. By a slight adaptation of routine A, we
can compute even sharper bounds on the value
W,. This allows us to skip several additional val
ues of W6. Once six weights have been fixed, it is
an easy job to find the interesting values for W,.
At this stage, almost every partial solution
indeed leads to a sequence which represents each
integer in 1 ... 0.
3.4 Remarks on the computer
program
The program has been coded in C and took 475
minutes of CPU time on a 450 MHz PC run
ning under Linux. The final search tree has
approximately 68 million leaf nodes. The last 5
tuple (W,, W2, W, W, W5) that was investigat
ed is (57, 73, 84, 91, 95). The highest values
considered for the seven weights are 57, 73, 85,
96, 109, 110, and 103, respectively. There exist
116,315,910 5tuples (W, < W2 W3 < W4 <
W), with W below the respective upper
bounds, so at least half of them have been fath
omed by our bounding procedures. This shows
that the search tree has been kept under control,
as far as the higher weights are concerned.
We have not been able to cut down the range
for the first weights W, andW2. The reader may
want to observe from the entries in Table 2 that
good solutions are possible without using the
values 1 and 2 at all. Nevertheless, one would
expect some argument showing that the partial
solutions in which, say, W > 10 need not be
considered. Such a result would speed up the
algorithm considerably. The highest upper
bounds on W6 and W, (observed when the first
five weights have already been fixed) are 137 and
262, respectively. This indicates that it will be
rather difficult to bound the ranges of these vari
ables a priori.
We conclude that only the application of the
various bounds enabled us to implicitly enumer
ate all solutions.
4 Concluding remarks
Table 2 is the result of the search setting the
threshold 0 = 120. It therefore contains all possi
ble solutions from which a consecutive range
1 .. 120 can be represented.
This type of weighing problem should yield
nice challenges for the constraint programming
community. After all, we have invested a lot of
work and a lot of thinking into the solution of
this problem, and it is not clear whether the
generic constraint programming approaches
(where one only needs to formulate a problem
in some metalanguage and then all the opti
mization and all the searching are performed
automatically by the underlying software) can
compete against human specialization.
Acknowledgement
We thank Clive Tooth for stating this weighing
problem on his web page at
(http://www.pisquaredoversix.force9.co.uk/
3From7.htm).
10PTIMA65
PAGE 7
MAY 2001
Solving Optimization Problems on Computational Grids
Stephen J. Wright*
November 21, 2000
1 Multiprocessors and
Computational Grids
Multiprocessor computing platforms, which
have become available more and more widely
since the mid1980s, are now heavily used by
organizations that need to solve very demanding
computational problems. There has also been a
great deal of research on computational tech
niques that are suited to these platforms, and on
the software infrastructure needed to compile
and run programs on them. Parallel computing
is now central to the culture of many research
communities, in such areas as meteorology,
computational physics and chemistry, and cryp
tography. Nevertheless, fundamental research in
numerical techniques for these platforms
remains a major topic of investigation in numer
ical PDE and numerical linear algebra.
The nature of parallel platforms has evolved
rapidly during the past 15 years. The Eighties
saw a profusion of manufacturers (Intel,
Denelcor, Alliant, Thinking Machines, Convex,
Encore, Sequent, among others) with a corre
sponding proliferation of architectural features:
hypercube, mesh, and ring topologies; shared
and distributed memories; memory hierarchies
of different types; butterfly switches; and global
buses. Compilers and runtime support tools
were machine specific and idiosyncratic.
Argonne's Advanced Computing Research
Facility kept a zoo of these machines during the
late 1980s, allowing free access to many
researchers in the United States and giving many
of us our first taste of this brave new world.
By the early 1990s, the situation had started
to change and stabilize. Most of the vendors
went out of business, and their machines gradu
ally were turned off one processor at a time, in
some cases. Architectures gravitated toward two
easily understood paradigms that prevail today.
One paradigm is the sharedmemory model typ
ified by the SGI Origin series and by computers
manufactured by Sun and HewlettPackard. The
other paradigm, typified by the IBM SP series,
can be viewed roughly as a uniform collection of
processors, each with its own memory and all
able to pass messages to one another at a rate
roughly independent of the locations of the two
processors involved. On the software side, the
advent of software tools such as p4, MPI, and
*Mathematics and Computer Science Division,
Argonne National Laboratory, and Computer
Science Department, University of Chicago.
PVM allowed users to write code that could be
compiled and executed without alteration on the
machines of different manufacturers, as well as
on networks of workstations.
The optimization community was quick to
take advantage of parallel computers. In design
ing optimization algorithms for these machines,
it was best in some cases to exploit parallelism at
a lower level (that is, at the level of the linear
algebra or the function/derivative evaluations)
and leave the control flow of the optimization
algorithm essentially unchanged. Other cases
required a complete rethinking of the algo
rithms, to allow simultaneous exploration of dif
ferent regions of the solution space, different
parts of the branchandbound tree, or different
candidates for the next iterate. Novel parallel
approaches were developed for global optimiza
tion, network optimization, and directsearch
methods for nonlinear optimization. Activity
was particularly widespread in parallel branch
andbound approaches for various problems in
combinatorial and network optimization, as a
brief Web search can attest.
As the cost of personal computers and low
end workstations has continued to fall, while the
speed and capacity of processors and networks
have increased dramatically, "cluster" platforms
have become popular in many settings. Though
the software infrastructure has yet to mature,
clusters have made supercomputing inexpensive
and accessible to an even wider audience.
A somewhat different type of parallel com
puting platform known as a computational grid
(alternatively, metacomputer) has arisen in com
paratively recent times [1]. Broadly speaking,
this term refers not to a multiprocessor with
identical processing nodes but rather to a het
erogeneous collection of devices that are widely
distributed, possibly around the globe. The
advantage of such platforms is obvious: they
have the potential to deliver enormous comput
ing power. (A particular type of grid, one made
up of unused computer cycles of workstations
on a number of campuses, has the additional
advantage of costing essentially nothing.) Just as
obviously, however, the complexity of grids
makes them very difficult to use. The software
infrastructure and the applications programs
that run on them must be able to handle the
following features:
heterogeneity of the various processors in
the grid;
the dynamic nature of the platform (the
pool of processors available to the user may
grow and shrink during the computation);
the possibility that a processor performing
part of the computation may disappear
without warning;
latency (that is, time for communications
between the processors) that is highly vari
able but often slow.
In many applications, however, the potential
power and/or low cost of computational grids
make the effort of meeting these challenges
worthwhile. The Condor team, headed by
Miron Livny at the University of Wisconsin,
were among the pioneers in providing infra
structure for grid computations. The Condor
system can be used to assemble a parallel plat
form from workstations, PC clusters, and multi
processors and can be configured to use only
"free" cycles on these machines, sharing them
with their respective owners and other users.
More recently, the Globus project has developed
technologies to support computations on geo
graphically distributed platforms consisting of
highend computers, storage and visualization
devices, and other scientific instruments.
In 1997, we started the metaneos project as a
collaborative effort between optimization spe
cialists and the Condor and Globus groups. Our
aim was to address complex, difficult optimiza
tion problems in several areas, designing and
implementing the algorithms and the software
infrastructure needed to solve these problems on
computational grids. A coordinated effort on
both the optimization and the computer science
sides was essential. The existing Condor and
Globus tools were inadequate for direct use as a
base for programming optimization algorithms,
whose control structures are inevitably more
complex than those required for taskfarming
applications. Many existing parallel algorithms
for optimization were "not parallel enough" to
exploit the full power of typical grid platforms.
Moreover, they were often "not asynchronous
enough," in that they required too much com
munication between tasks to execute efficiently
on platforms with the heterogeneity and com
munications latency properties of our target
platforms. A further challenge was that, in con
trast to other grid applications, the computa
tional resources required to solve an optimiza
10S T I M 6
PAGE 8
MAY 2001
tion problem often cannot be predicted with
much confidence, making it difficult to assemble
and utilize these resources effectively.
This article describes some of the results we
have obtained during the first three years of the
metaneos project. Our efforts have led to devel
opment of the runtime support library MW, for
implementing algorithms with masterworker
control structure on Condor platforms. This
work is discussed below, along with our work on
algorithms and codes for integer linear program
ming, the quadratic assignment problem, and
stochastic linear programming. Other metaneos
work, not discussed below, includes work in
global optimization, integer nonlinear program
ming, and verification of solution quality for
stochastic programming.
2 Condor, Globus, and the MW
Framework
The Condor system [2, 3] had its origins at the
University of Wisconsin in the 1980s. It focuses
on collections of computing resources, known as
Condorpools, that are distributively owned. To
understand the implications of "distributed
ownership," consider a typical machine in a
pool: a workstation on the desk of a researcher.
The Condor system provides a means by which
other users (not known to the machine's owner)
can exploit some of the unused cycles on the
machine, which otherwise would sit idle most of
the time. The owner maintains control over the
access rights of Condor to his machine, specify
ing the hours in which Condor is allowed to
schedule processes on the machine and the con
ditions under which Condor must terminate any
process it is running when the owner starts a
process of his own. Whenever Condor needs to
terminate a process under these conditions, it
migrates the process to another machine in the
pool, guaranteeing eventual completion.
When a user submits a process, the Condor
system finds a machine in the pool that matches
the software and hardware requirements of the
user. Condor executes the user's process on this
machine, trapping any system calls made by the
process (such as input/output operations) and
referring them back to the submitting machine.
In this way, Condor preserves much of the sub
mitting machine's environment on the execution
machine. Users can submit a large number of
processes to the pool at once. Since each such
process maintains contact with the submitting
machine, this feature of Condor opens up the
possibility of parallel processing. Condor pro
vides an opportunistic environment, one that can
make use of whatever resources currently are
available in its pool. This set of resources grows
and shrinks dynamically during execution of the
user's job, and his algorithm should be able to
exploit this situation.
The Globus Toolkit [4] is a set of components
that can be used to develop applications or pro
gramming tools for computational grids.
Currently, the Toolkit contains tools for resource
allocation management and resource discovery
across a grid, security and authentication, data
movement, messagepassing communication,
and monitoring of grid components. The main
use of Globus within the metaneos project has
been at a level below Condor. By a Globus
mechanism known as glidein, a user can add
machines at a remote location into the Condor
pool on a temporary basis, making them accessi
ble only to his own processes. In this way, a user
can marshal a large and powerful set of resources
over multiple sites, some or all of them dedicat
ed exclusively to his job.
MW is a software framework that facilitates
implementation of algorithms of masterworker
type on computational grids. It was developed as
part of the metaneos project by Condor team
members Mike Yoder and Sanjeev Kulkarni in
collaboration with optimization specialists Jeff
Linderoth and JeanPierre Goux [5, 6]. MW
takes the form of a set of C++ abstract classes,
which the user implements to perform the par
ticular operations associated with his algorithm
and problem class. There are just ten virtual
functions, grouped into the following three fun
damental base classes:
MWDriver contains four functions that
obtain initial user information, set up the
initial set of tasks, pack the data required to
initialize each worker processor as it
becomes available, and act on the results
that are returned to the master when a task
is completed.
MWWorker contains two functions, to
unpack the initialization data for the work
er and to execute a task sent by the master.
MWTask contains four functions to pack
and unpack the data defining a single task
and to pack and unpack the results associ
ated with that task.
MW also contains functions that monitor per
formance of the grid and gather various statistics
about the run.
Internally, MW works by managing a list of
workers and a list of tasks. The resource manage
ment mechanisms of the underlying grid are
used to obtain new workers for the list and pro
vide information about each worker. The infor
mation can be used to order the worker list so
that the most suitable workers (e.g., the fastest
machines) are at the head of the list and hence
are the first to receive tasks. Similarly, the task
list can be ordered by a userdefined key to
ensure that the most important tasks are per
formed first. Scheduling of tasks to workers then
becomes simple: The first task on the list is
assigned to the first available worker. New tasks
are added to the list by the master process in
response to results received from completion of
an earlier task.
MW is currently implemented on two slightly
different grid platforms. The first uses Condor's
version of the PVM (parallel virtual machine)
protocol, while the second uses the remote I/O
features of Condor to allow master and workers
to communicate via series of shared files. In
addition, MW provides a "bottomlevel" inter
face that allows it to be implemented in other
grid computing toolkits.
3 Integer Programming
Consider the linear mixed integer programming
problem
min c x subject to Ax < b, I < x< u,
xE Z, for all E I,
where x is a vector of length n, Z represents the
integers, and IC {1,2,...,n}. Parallel algorithms
and frameworks for this problem have been
investigated by a number of authors in recent
times. The approaches described in [7, 8, 9, 10]
implement enhancements of the branchand
bound procedure, in which the work of explor
ing the branchandbound tree is distributed
among a fixed number of processors. These
approaches are differentiated by their use of vir
tualsharedmemory vs. messagepassing models,
their load balancing procedures, their choice of
branching rules, and their use of cuts.
By contrast, the FATCOP code of Chen,
Ferris, and Linderoth [11] uses the MW frame
work to greedily use whatever computational
resources become available from the Condor
pool. The algorithm implemented in FATCOP
(the name is a loose acronym for "faulttolerant
Condor PVM") is an enhanced branchand
bound procedure that utilizes (globally valid)
cuts, pseudocosts for branching, preprocessing at
nodes within the branchandbound tree, and
heuristics to identify integer feasible solutions
rapidly.
In FATCOP's masterworker algorithm, each
task consists of exploration of a subtree of the
branchandbound tree, not just evaluation of a
single node of the tree. Given a root node for
the subtree, and other information such as the
global cut and pseudocost pools, the task exe
cutes for a given amount of time, making its
own branching decisions and accumulating its
10PTS A65
PAGE 9
MAY 2001
own collection of cuts and pseudocosts. It may
also perform a "diving" heuristic from its root
node to seek a new integer feasible solution.
When the task is complete, it returns to the
master a stack representing the unexplored por
tions of its subtree. (Depthfirst search is used to
limit the size of this stack.) The task also sends
back any new cut and pseudocost information it
generated, which is added to the master's global
cut and pseudocost pools.
By processing a subtree rather than a single
node, FATCOP increases the granularity of the
task and improves utilization of the computa
tional power of each worker. The time for which
a processor is sitting idle while waiting for the
task information to be sent to and from the
master, and for the master to process its results
and assign it a new task, is generally small rela
tive to the computation time.
The master process is responsible for main
taining the task pool, as well as the pools of cuts
and pseudocosts. It recognizes new workers as
they join the computation pool, and sends them
copies of the problem data together with the
current cut and pseudocost pools. Moreover, it
sends tasks to these workers and processes the
results of these tasks by updating its pools and
possibly its incumbent solution. Another impor
tant function of the master is to detect when a
machine has disappeared from the worker pool.
Ta
F
minumum
average
maximum
In this case, th
machine is losi
another work
ture that make
stations, the ma
writing out th
disk. By doing
at the latest ch
ter processor.
To illustrate
sider the solut
the MIPLIB te
electricity gen
Islands. FATC
Condor pool a
Because of the
ment the si
available work
the network va
each run th
COP was quit
different from what one would obtain from a
serial implementation of the same approach.
However, using the statistical features of MW,
we can find the average number of workers used
during each run, defined as
P = k / T,
k=1
where tk is the total time during which the run
had control of k processors, and T is the wall
clock time for the run. The minimum, maxi
mum, and mean values of P over the ten runs
are shown in Table 1. This table also shows sta
tistics for the number of nodes evaluated by
FATCOP, and the wall clock times.
Figure 1: Number of Workers
during FATCOP on gesa2_o
Time
ble 1: Performance of Figure 1 profiles the size of the worker pool
during a particular run. Note the sharp dip,
4TCOP on gesa2_o
which occurred when a set of machines partici
P Nodes Time pated in a daily backup procedure, and the grad
43.2 6207993 6951 ual buildup during the end of the run, which
62.5 8214031 10074 occurred in the late afternoon when more
86.3 9693518 13198 machines became available as their owners went
home.
e task that was occupying that
For a detailed performance analysis of FAT
t, and the master must assign it to
COP, see [11].
r. (This is the "fault tolerant" fea C se
A separate but related activity involved solv
s FATCOP fat!) On long compu
e going to optimality the wellknown seymour prob
aster process checkpoints" by
lem from the MIPLIB library of integer pro
e current state of computation to
grams. This wellknown integerprogramming
;so, it can restart the computation
oit ca restart the com tat problem, posed by Paul Seymour, arises in a new
eckpoint after a crash of the mas
proof [12] of the famous FourColor Theorem,
which states that any map can be colored using
FATCOP's performance, we con
SFA s performance, con four colors in such a way that regions sharing a
ion of the problem gesa2_o from
Sp boundary segment receive different colors. The
;st set. This problem arises in an
Balearic seymour problem is to find the smallest set of
ration application in the Balearic
P was rn tn tien the configurations such that the FourColor
OP was run ten times on the
Theorem is true if none of these configurations
it the University of Wisconsin.
can exist in a minimal counterexample.
dynamic computational environ
dynamc compuational e Although Seymour claimed to have found a
ze and composition of the pool of
solution with objective value 423, nobody
ers and communication times on
d (including Seymour himself) had been able
tried between runs and during
rI, ed11 bete, rn r either to reproduce this solution, or prove a
el sadrlall owd.LL1 III WU Uy o 1 
e different in each instance, and
strong lower bound on the optimal value. There
PAGE 10
was some skepticism in the integer program
ming community as to whether this was indeed
the optimal value.
In July 2000, a team of metaneos researchers
found solutions with the value 423, and proved
its optimality. Gabor Pataki, Stefan Schmieta,
and Sebastian Ceria at Columbia used prepro
cessing, disjunctive cuts and branchandbound
to break down the problem into a list of 256
integer programs. Michael Ferris at Wisconsin
and Jeff Linderoth at Argonne joined the
Columbia group in working through this list.
The problems were solved separately with the
help of the Condor system, using the integer
programming packages CPLEX and XPRESS
MP. About 9,000 hours of CPU time was need
ed, the vast majority of it spent in proving opti
mality.
4 Quadratic Assignment Problem
The quadratic assignment problem (QAP) is a
problem in location theory that has proved to be
among the most difficult combinatorial opti
mization problems to solve in practice. Given
n x n matrices A, B, and C, where A, represents
the flow between facilities i andj, B is the dis
tance between locations i andj, and C is the
fixed cost of assigning facility i to location, the
problem is to find the permutation {r(1),
(2), . ,T(n)} of the index set {1, 2, . n}
that minimizes the following objective:
;=1 i=1 ;=1
An alternative matrixform representation is as
follows:
QAP(A,B,C,): min tr(AXB + C) XT,
s.t. X II,
where tr(.) represents the trace and H is the set
of n x n permutation matrices.
The practical difficulty of solving instances of
the QAP to optimality grows rapidly with n. As
recently as 1998, only the secondlargest prob
lem (n = 25) from the standard "Nugent"
benchmark set [13] had been solved, and this
effort required a powerful parallel platform [14].
In June 2000, a team consisting of Kurt
Anstreicher and Nate Brixius (University of
Iowa) and metaneos investigators JeanPierre
Goux and Jeff Linderoth solved the largest of
the Nugent problems the n = 30 instance
known as nug30 verifying that a solution
obtained earlier from a heuristic was optimal
[15]. They devised a branchandbound algo
rithm based on a convex quadratic programming
10PS IM 65
MAY 2001
relaxation of QAP, implemented it using MW,
and ran it on a Condorbased computational
grid spanning eight institutions. The computa
tion used over 1,000 worker processors at its
peak, and ran for a week. It was solving linear
assignment problems (the core computational
operation in the algorithm) at the rate of nearly
a million per second during this period.
In the remainder of this section, we outline
the various theoretical, heuristic, and computa
tional ingredients that combined to make this
achievement possible.
The convex quadratic programming (QP)
relaxation of QAP proposed by Anstreicher and
Brixius [16] yields a lower bound on the optimal
objective that is tighter than alternative bounds
based on projected eigenvalues or linear assign
ment problems. Just as important, an approximate
solution to the QP relaxation can be found at
reasonable cost by applying the FrankWolfe
method for quadratic programming. Each itera
tion of this method requires only the solution of
a dense linear assignment problem an inexpen
sive operation. Hence, the FrankWolfe method
is preferable in this context to more sophisticated
quadratic programming algorithms; its slow
asymptotic convergence properties are not
important because only an approximate solution
is required.
The QAP is solved by embedding the QP
relaxation scheme in a branchandbound strate
gy. At each node of the branchandbound tree,
some subset of the facilities is assigned to certain
locations in the nodes at level k of the tree,
exactly k such assignments have been made. At a
levelk node, a reduced QAP can be formulated
in which the unassigned part of the permutation
(which has n k components) is the unknown.
The QP relaxation can then be used on this
reduced QAP to find an approximate lower
bound on its solution, and therefore on the cost
of all possible permutations that include the k
assignments already made at this node. If the
bound is greater than the cost of the best permu
tation found to date (the incumbent), the sub
tree rooted at this node can be discarded.
Otherwise, we need to decide whether and how
to branch from this node.
Branching is performed by choosing a facility
and assigning it to each location in turn (row
branching) or by choosing a location and assign
ing each of the remaining facilities to it in turn
(column branching). However, it is not always
necessary to examine all possible n k children
of a levelk node; some of them can be eliminat
ed immediately by using information from the
dual of the QP relaxation. In fact, one rule for
deciding the next node from which to branch at
level k of the tree is to choose the node that
PAGE 11
yields the fewest children. A more expensive
branching rule, using a strong branching tech
nique, is employed near the root of the tree (k
smaller). Here, the consequence of fixing each
one of a collection of promising facilities (or
locations) is evaluated by provisionally making
the assignment in question and solving the cor
responding QP relaxation. Estimates of lower
bounds are then obtained for the grandchildren
of the current node, and these are summed. The
node for which this summation is largest is cho
sen as the branching facility (location).
The branching rule and the parameters that
govern the execution of the branching rule are
chosen according to the level in the tree and also
the gap, which measures the closeness of the
lower bound at the current node to the incum
bent objective. When the gap is large at a partic
ular node, it is likely that exploration of the sub
tree rooted at this node will be a costly process.
Use of a more elaborate (and expensive) branch
ing rule tends to ensure that exploration of
unprofitable parts of the subtree is avoided,
thereby reducing overall run time.
Parallel implementation of the branchand
bound technique uses an approach not unlike
the FATCOP code for integer programming.
Each worker is assigned the root node of a sub
tree to explore, in a depthfirst fashion, for a
given amount of time. When its time expires, it
returns unexplored nodes from its subtree to the
master, together with any new incumbent infor
mation. The pool of tasks on the master is
ordered by the gap, so that nodes with smaller
gaps (corresponding to subtrees that should be
less difficult to explore) are assigned first. To
reduce the number of easy tasks returned to the
master, a "finishup" heuristic permits a worker
extra time to explore its subtree if its gap
becomes small.
Exploitation of the symmetries that are pres
ent in many large QAPs is another important
factor in making solution of nug30 and other
large problems a practical proposition. Such
symmetries arise when the distance matrix is
derived from a rectangular grid. Symmetries can
be used, for instance, to decrease the number of
child nodes that need to be formed (to consider
ably fewer than n k children at a levelk node).
Prediction of the performance profile of a run
is also important in tuning algorithmic parame
ters and in estimating the amount of computa
tional resources needed to tackle the problem.
An estimation procedure due to Knuth was
enhanced to allow prediction of the number of
nodes that need to be evaluated at each level of
the branchandbound tree, for a specific prob
lem and specific choices of the algorithmic
parameters.
Figure 2 shows the number of workers used
during the course of the nug30 run in early June
2000. As can be seen from this graph, the run
was halted five times twice because of failures
in the resource management software and three
times for system maintenance and restarted
each time from the latest master checkpoint.
In the weeks since the nug30 computation,
the team has solved three more benchmark prob
lems of size n = 30 and n = 32, using even larger
computational grids. Several outstanding prob
lems of size n = 36 derived from a backboard
wiring application continue to stand as a chal
lenge to this group and to the wider combinato
rial optimization community.
Figure 2: Number of Workers during
nug30 Computation
1000 .
soo 0 i
600 Y
69 6 0 6/ 11 612 6/13 6/14 615
Time
5 Stochastic Programming
The twostage stochastic linear programming
problem with recourse can be formulated as fol
lows:
def N
min Q(x) = x+ piQ (x)
subject to Ax = b, x > 0,
where
def
Qi(x) = mingT y s.t. Wy = hi Tx, Yi 0.
Yi
The uncertainty in this formulation is modeled
by the data triplets (h, T, q,), i = 1, 2, .. ,N,
each of which represents a possible scenario for
the uncertain data (h, T q). Each p, represents
the probability that scenario i is the one that
actually happens; these quantities are nonnega
tive and sum to 1. The quantities p,, i =
1, 2 .. ., Nare nonnegative and sum to 1; p
is the probability that scenario i is the true one.
The twostage problem with recourse repre
sents a situation in which one set of decisions
(represented by the firststage variables x) must
be made at the present time, while a second set
of decisions (represented byy, i = 1, 2, . ., )
can be deferred to a later time, when the uncertain
ty has been resolved and the true secondstage sce
nario is known. The objective function represents
0PTIMA65
MAY 2001
the expected cost of x, integrated over the probability
distribution for the uncertain part of the model.
In many practical problems, the number of
possible scenarios Neither is infinite (that is, the
probability distribution is continuous) or is
finite but much too large to allow practical solu
tion of the full problem. In these instances, sam
pling is often used to obtain an approximate
problem with fewer scenarios.
Decomposition algorithms are well suited to
grid platforms, because they allow the computa
tions associated with the N secondstage scenar
ios to be performed independently and require
only modest amounts of data movement
between processors. These algorithms view Q as
a piecewise linear, convex function of the vari
ables x, whose subgradient is given by the
formula
N
aQ(x)= c+ YpaQ (x).
i=1
Evaluation of each function Q, i = 1, 2,..., N
requires solution of the linear program in y,
given above. One element of the subgradient
aQ (x) of this function can be calculated from
the dual solution of this linear program.
In the metaneos project, Linderoth and
Wright [17] have implemented a decomposition
algorithm based on techniques from non
smooth optimization and including various
enhancements to lessen the need for the master
and workers to synchronize their efforts. In the
remainder of this section, we outline in turn the
trustregion algorithm ATR and its convergence
properties, implementation of this algorithm on
the Condor grid platform using MW, and our
"asynchronous" variant.
The ATR algorithm progressively builds up a
piecewiselinear model function M(x) satisfying
M(x) < Q (x). At the kth iteration, a candidate
iterate z is obtained by solving the following
subproblem:
min, M(z) subject to Az = b, z > 0,
z x' < A,
where the last constraint represents a trust region
with radius A > 0. The candidate z becomes the
new iterate xk if the decrease in objective
Q(x) Q(z) is a significant fraction of the
decrease Q(x) M(z) promised by the model
function. Otherwise, no step is taken. In either
case, the trustregion radius A may be adjusted,
function and subgradient information about Q
at z is used to enhance the model M, and the
subproblem is solved again. The algorithm uses a
"multicut" variant in which subgradients for
partial sums of YXI (U) are included in the
model separately, allowing a more accurate
model to be constructed in fewer iterations.
In the MW implementation of the ATR algo
rithm, the function and subgradient information
defining M is accumulated at the master proces
sor, and the subproblem is solved on this proces
sor. (Since M is always piecewise linear and the
trust region is defined by an oonorm, the sub
problem can be formulated as a linear program.)
Most of the computational work in the algo
rithm involves solution of the N secondstage
linear programs in y, from which we obtain
Q(z) and aQ(z). This work is distributed
among Ttasks, to be executed in parallel, where
each task requires solution of a "chunk" of N/T
secondstage linear programs.
The use of chunking allows problems with
very large Nto make efficient use of a fairly
large number of processors. However, the
approach still requires evaluation of all the
chunks for Q(z) to be completed before deciding
whether to accept or reject z as the next iterate.
It is possible that one chunk will be processed
much more slowly than the othersits computa
tion may have been interrupted by the worksta
tion's owner reclaiming the machine, for
instance. All the other workers in the pool will
be left idle while waiting for evaluation of this
chunk to complete.
The ATR method maintains not just a single
candidate for the next iterate but rather a basket
1 containing 5 to 20 possible candidates. At any
given time, the workers are evaluating chunks of
secondstage problems associated with one or
other of these basket points. ATR also maintains
an "incumbent" x, which is the current best
estimate of the solution and is a point for which
Q(x') is known. When all the chunks for one of
the basket points z have been evaluated,Q(z) is
compared with the incumbent objective Q(x')
and with the decrease predicted by the model
function M at the time z was generated. As a
result, either z becomes the new incumbent and
x is discarded, or x remains the incumbent and
z is discarded. In either case, a vacancy is created
in the basket 1. To fill the vacancy, a new candi
date iterate z' is generated by solving a sub
problem with the trustregion constraint cen
tered on the incumbent, that is,
1z x' <_A.
We show results obtained for sampled
instances of problems from the stochastic pro
gramming literature, using the MW implemen
tation of ATR running on a Condor pool. The
SSN problem described in [18] arises in design
of a network for privateline telecommunications
services. In this model, each of 86 parameters
representing demand can independently take on
3 to 7 values, giving a total of approximately
PAGE 12
N= 1' scenarios. Sampling is used to obtain
problems with more modest values of N, which
are then solved with ATR.
Table 2: SSN, with N= 10,000
Run Iter. Procs. Eff. Time(min.)
L 255 19 .46 398
ATR1 47 19 .35 130
ATR10 164 71 .57 43
Results for an instance of SSN with
N= 10000 are shown in Table 2. When written
out as a linear program in the unknowns
(x y Y2 . YN), this problem has approxi
mately 1,750,000 rows and 7,060,000 columns.
Table 2 compares three algorithms. The first is
an Lshaped method (see [19]), which obtains
its iterates from a model function M but does
not use a trust region or check sufficient
decrease conditions. (The implementation
described here is modified to improve paral
lelism, in that it does not wait for all the chunks
for the current point to be evaluated before cal
culating a new iterate.) The second entry in
Table 2 is for the synchronous trustregion
approach (which is equivalent to ATR with a
basket size of 1), and the third entry is for ATR
with a basket size of 10. In all cases, the second
stage evaluations were divided into 10 chunks,
and 50 partial subgradients were added to M at
each iteration. The table shows the average num
ber of processors used during the run, the pro
portion of time for which these processors were
kept busy, and the wall clock time required to
find the solution.
The trustregion approaches were consider
ably faster than the Lshaped approach, indicat
ing that the need for sound algorithms remains
as keen as ever in a parallel environment; we
cannot rely on raw computing power to do all
the work. The benefits of asynchronicity can also
be seen. When ATR has a basket size of 10, it is
able to use a larger number of processors and
takes less time to complete, even though the
number of iterates increases significantly over the
synchronous trustregion approach.
The real interest lies, however, not in solving
single sampled instances of SSN, but in obtain
ing highquality solutions to the underlying
problem (the one with 1070 scenarios). ATR
gives a valuable tool that can be used in con
junction with variance reduction and verification
techniques to yield such solutions.
Table 3: storm, with N= 250,000
Run Iter. Procs. Eff. Time(min.)
ATR1 25 67 .57 211
ATR5 57 86 .96 229
10 T I M A
MAY 2001
Finally, we show results from the "storm"
problem, which arises from a cargo flight sched
uling application [20]. The ATR implementa
tion was used to solve a sampled instance with
N= 250,000, for which the full linear program
has approximately 132,000,000 rows and
315,000,000 columns. The results in Table 3
show that this huge linear program with nontriv
ial structure can be solved in less than 4 hours
on a computational platform that costs essential
ly nothing. Because the secondstage work can
be divided into a much larger number of chunks
than for SSN 125 chunks, rather than 10 
the synchronous trustregion algorithm is able to
make fairly effective use of an average of 67
processors and requires less wall clock time than
ATR with a basket size of 5.
Conclusions
Our experiences in the metaneos project have
shown that cheap, powerful computational grids
can be used to tackle large optimization prob
lems of various types. These results have several
interesting implications. In an industrial or com
mercial setting, the results demonstrate that one
may not have to buy powerful computational
servers to solve many of the large problems aris
ing in areas such as scheduling, portfolio opti
mization, or logistics; the idle time on employee
workstations (or, at worst, an investment in a
modest cluster of PCs) may do the job. For the
optimization research community, our results
motivate further work on parallel, gridenabled
algorithms for solving very large problems of
other types. The fact that very large problems
can be solved cheaply allows researchers to better
understand issues of "practical" complexity and
of the role of heuristics. In stochastic optimiza
tion, higherquality solutions can be found, and
improvements to sampling methodology can be
investigated.
Work remains to be done in making the grid
infrastructure robust enough for general use. The
logistics of assembling a grid issues of security
and shared ownership remain challenging. The
vision of a computational grid that is as easy to
tap into as the electric power grid remains far
off, though metaneos gives a glimpse of the way
in which optimizers could exploit such a system.
We have investigated just a few of the prob
lem classes that could benefit from solution on
computational grids. Global optimization prob
lems of different types should be examined fur
ther. Dataintensive applications (from tomogra
phy and data mining) represent a potentially
huge field of work, but these require a somewhat
different approach from the computeintensive
applications we have considered to date.
We hope that optimizers of all flavors, along
with grid computing experts and applications
specialists, will join the quest. There's plenty of
work for all!
Acknowledgements
The metaneos project was funded by National
Science Foundation grant CDA9726385 and
was supported by the Mathematical,
Information, and Computational Sciences
Division subprogram of the Office of Advanced
Scientific Computing Research, U.S.
Department of Energy, under Contract W31
109Eng38.
For more information on the metaneos proj
ect, see .
Information on Condor can be found at
and on Globus at
.
References
[1] I. Foster and C. Kesselman, edi
tors. The Grid: Blueprintfor a
New Computing Infrastructure.
MorganKaufmann, 1998.
[2] J. Basney, M. Livny, and T.
Tannenbaum. High throughput
computing with Condor. HPCU
News, 1(2), 1997.
[3] M. Livny, J. Basney, R. Raman,
and T. Tannenbaum. Mechanisms
for high throughput computing.
SPEEDUP, 11(1), June 1997.
[4] I. Foster and C. Kesselman. The
Globus project: A status report. In
Proceedings of the Heterogeneous
Computing Workshop, pages 418.
IEEE Press, 1998.
[5] J.P Goux, J. T. Linderoth, andM.
Yoder. Metacomputing and the
masterworker paradigm. Preprint
MCS/ANLP7920200,
Mathematics and Computer
Science Division, Argonne
National Laboratory, Argonne,
Ill., February 2000.
[6] J.P. Goux, S. Kulkarni, J. T.
Linderoth, and M. Yoder. An
enabling framework for
masterworker applications for the
computational grid. In Proceedings
of the Ninth IEEE Symposium on
High Performance Distributed
Computing (HPDC9), pages 43
50, Pittsburgh, Pennsylvania,
August 2000.
[7] J. Eckstein. Parallel branchand
bound methods for mixedinteger
programming on the CM5.
SIAMJournal on Optimization,
4(4):794814, 1994.
[8] R. E. Bixby, W Cook, A. Cox,
and E. K. Lee. Computational
experience with parallel mixed
integer programming in a distrib
uted environment. Annals of
Operations Research, 90:1943,
1999.
[9] J. T. Linderoth. Topics in Parallel
Integer Optimization. Ph.D. thesis,
School of Industrial and Systems
Engineering, Georgia Institute of
Technology, 1998.
[10] J. Eckstein, C. A. Phillips, and W
E. Hart. PPICO: An object
oriented framework for parallel
branch and bound. Research
Report RRR402000, RUTCOR,
August 2000.
[11] Q. Chen, M. C. Ferris, and J. T.
Linderoth. FATCOP 2.0:
Advanced features in an oppor
tunistic mixed integer program
ming solver. Technical Report
Data Mining Institute Technical
Report 9911, Computer Sciences
Department, University of
Wisconsin at Madison, 1999. To
appear in Annals of Operations
Research.
[12] N. Robertson, D. P. Sander, P. D.
Seymour, and R. Thomas. A new
proof of the four color theorem.
Electronic Research Announcements
ofthe American Mathematical
Society, 1996.
[13] C. E. Nugent, T. E. Vollman, and
J. Ruml. An experimental compar
ison of techniques for the assign
ment of facilities to locations.
Operations Research, 16:150173,
1968.
[14] A. Marzetta and A. Brtinger. A
dynamicprogramming bound for
the quadratic assignment problem.
In Computing and Combinatorics:
5th Annual International
Conference, COCOON '99, vol
ume 1627 of Lecture Notes in
Computer Science, pages 339348.
Springer, 1999.
[15] K. Anstreicher, N. Brixius, J.P.
Goux, and J. T. Linderoth.
Solving quadratic assignment
problems on computational grids.
Technical report, Mathematics
and Computer Science Division,
Argonne National Laboratory,
Argonne, Ill., October 2000.
[16] K. M. Anstreicher and N. Brixius.
A new bound for the quadratic
assignment problem based on con
vex quadratic programming.
Technical report, Department of
Management Sciences, The
University of Iowa, 1999.
[17] J. T. Linderoth and S. J. Wright.
Implementing decomposition
algorithms for stochastic program
ming on a computational grid.
Technical report, Mathematics
and Computer Science Division,
Argonne National Laboratory,
Argonne, 11, 2000. In prepara
tion.
[18] S. Sen, R. D. Doverspike, and S.
Cosares. Network planning with
random demand. Telecommuni
cations Systems, 3:1130, 1994.
[19] J. R. Birge and E Louveaux.
Introduction to Stochastic
Programming. Springer Series in
Operations Research. Springer,
1997.
[20] J. M. Muley and A. Ruszczynski.
A new scenario decomposition
method for large scale stochastic
optimization. Operations Research,
43:477490, 1995.
1 P Ii 5
PAGE 13
MAY 2001
CONFERENCE
Calendar
1st Cologne Twente Workshop on Graphs
and Combinatorial Optimization
June 68, 2001, University of Cologne,
Cologne, Germany.
URL: http://www.zaik.unikoeln.de/AFS/
conferences/CTW2001/CTW2001.html
IPCO VIII
June 1315, 2001, Utrecht,
The Netherlands.
URL: http://www.cs.uu.nl/events/ipco2001
email: ipco2001@cs.uu.nl
INFORMS International 2001
June 1720, 2001, Maui, Hawaii, USA.
URL: http://www.informs.org/Conf/
Hawaii2001
SIAM Annual Meeting
July 913, 2001, San Diego,
California, USA.
URL: http://www.siam.org/meetings/an01/
Symposium on OR 2001
September 35, 2001, GerhardMercator
Universitdit, Duisburg, Germany.
URL: http://www.uniduisburg.de/or2001/
INFORMS Annual Meeting
November 47, 2001, Miami Beach,
Florida, USA.
URL: http://www.informs.org/Conf/
Miami2001
ANNOUNCEMENT
9th International Conference on
Stochastic Programming
Berlin, Germany, August 2531, 2001
(held under the auspices of the Committee on Stochastic Programming
of the Mathematical Programming Society)
Contact Information. Web site: http://www.mathematik.
huberlin.de/SPO1; Email: sp01@mathematik.huberlin.de;
Mailing Address: SP01, c/o Prof. W Roemisch, Institute of Mathematics,
HumboldtUniversity Berlin, 10099 Berlin, Germany
Support. HumboldtUniversitaet Berlin, GerhardMercatorUniversitaet
Duisburg, WeierstrassInstitut fuer Angewandte Analysis und Stochastik
(WIAS), KonradZuseZentrum fuer Informationstechnik Berlin (ZIB),
Deutsche Forschungsgemeinschaft
Scope. Stochastic programming addresses the optimization of decision
making under uncertainty. Several branches of mathematics and its prac
tical applications inspire it. Studies in stochastic programming are direct
ed to model building, theoretical analysis of models, design of algo
rithms, software development, and practical implementation. The 9th
International Conference on Stochastic Programming will be a forum to
discuss recent advances, the state of the art, and future developments of
the field.
Topics of the conference. modeling techniques, scenario generation,
twostage and chance constrained models, dynamic stochastic optimiza
tion, integer and combinatorial optimization under uncertainty, risk
management, stability analysis, estimation and simulation, approximation
and bounding, decomposition methods, stochastic programming models
in finance, and practical applications (e.g., in energy and logistics).
International Scientific Committee. J.R. Birge (USA), M.A.H.
Dempster (Great Britain), J. Dupacova (Czech Republic), P. Kall
(Switzerland), A.J. King (USA), A. Prekopa (Hungary/USA), G. Pflug
(Austria), W. Roemisch (Germany), A. Ruszyriski (USA), R. Schultz
(Germany), S. Sen (USA), A. Shapiro (USA), S.W. Wallace (Norway), R.
JB Wets (USA), W.T Ziemba (Canada)
Invited Plenary Speakers. Hans Foellmer (Berlin), Suvrajeet Sen
(Tucson), Alexander Shapiro (Atlanta), Stein W Wallace (Trondheim),
Roger JB Wets (Davis)
Organizing Committee. K. Frauendorfer (St. Gallen), R. Henrion
(WIAS Berlin), P Kall (Zuerich), K. Marti (Muenchen), G. Pflug
(Wien), W Roemisch (HU Berlin), R. Schultz (Duisburg), M.C.
Steinbach (ZIB Berlin), S. Vogel (Ilmenau)
Local Organizers. W Roemisch, R. Schultz, N. GroeweKuska, R.
Henrion, J. Kerger, H. Lange, A. Moeller, I. Nowak, M.C. Steinbach, I.
Wegner
Location. HumboldtUniversity Berlin (Main Building), Unter den
Linden 6, 10117 Berlin
Schedule. Tutorials: August 2526, 2001; International Conference:
August 2731, 2001; Opening Session: August 27, 2001, 9:30 am
Fees. Conference fee: 400. DM (before May 31, 2001), 450. DM;
150. DM (student); Tutorial fee: 100. DM
Deadlines. February 15, 2001 Starting Online Registration; April 15,
2001 Abstract Submission; April 30, 2001 Notification of Acceptance;
May 31, 2001 Early Registration
I P T I
PAGE 14
MAY 2001
EURO 2001: The XVIIIth Euro
Conference on Operations
Research
Why Visit Euro 2001: Meeting OR friends in
convenient and wellequipped conference center;
participating in and contributing to high quality
theory and practice sessions; learning more
about smart logistics and innovative operations
during special conference sessions, company vis
its, and excursions to the largest port in the
world; reduced fee to a post conference seminar
on Quantitative Financial Risk Management in
Amsterdam; enjoying social events in Europe's
cultural capital of 2001
Where: Erasmus University Rotterdam, The
Netherlands
When: July 911, 2001
How to Register. The most convenient way to
register is via our web site (www.euro2001.org).
Alternatively, you may order a registration card
via our email address (info@euro2001.org).
Fees: Early: Euro 300 for NonStudents, Euro
200 for Students; Late: Euro 350 for Non
Students, Euro 250 for Students.
Abstract Submission: Please refer to our web site
for submission instructions, or ask for instruc
tions on paper via our email address
(info@euro2001.org).
Further Information: Email: info@euro2001.org;
web site: www.euro2001.org
Deadlines: The early registration deadline is
May 1, 2001; the abstract submission deadline is
March 1, 2001. Please note that submissions
that have been received after the deadline cannot
be accepted.
Call for Papers
MATHEMATICAL PROGRAM
MING Series B
Issue on Algebraic and Geometric
Methods in Discrete Optimization
PAGE 15
Abstracts can be submitted via the web site
(www.euro2001.org), which also contains a list of
topics and submission instructions.
We hope to welcome you at EURO 2001!
www.euro2001.org
Ettore Majorana Centre for
Scientific Culture International
School of Mathematics G.
Stampacchia Workshop
High Performance Algorithms and Software
for Nonlinear Optimization
Erice, Italy
June 30 July 8, 2001
Objectives: In the first year of a new century, the
Workshop aims to assess the past and to discuss
the future of Nonlinear Optimization. Recent
achievements and promising research trends will
be highlighted. Algorithmic and high perform
ance software aspects will be emphasized, along
with theoretical advances and new computation
al experiences. Contributions from different and
complementary standpoints, covering several
fields of numerical optimization and its applica
tions are expected.
Topics: Topics include, but are not limited to:
constrained and unconstrained optimization,
global optimization, derivativefree methods,
interior point techniques for nonlinear program
ming, large scale optimization, linear and non
linear complementarity problems, nonsmooth
optimization, neural networks and optimization,
applications of nonlinear optimization.
Lectures: The Workshop will consist of invited
lectures (1 hour) and contributed lectures (1/2
hour). Invited lecturers who have confirmed the
participation are: P. Boggs, A. R. Conn, J.
We invite research articles on algebraic and geomet
ric methods in discrete optimization and related
areas for a forthcoming issue of Mathematical
Programming, Series B. The goal is to provide a
common forum for nontraditional methods in the
theory and practice of discrete optimization. Of
special interest are methods from the geometry of
numbers, topology, commutative algebra, algebraic
geometry, probability theory, computer science and
combinatorics.
Deadline for submission of full papers: April 30,
2001.
We aim at completing a first review of all papers
by September 30, 2001.
Dennis, Y. G. Evtushenko, F. Facchinei, R.
Fletcher, J. C. Gilbert, P. E. Gill, D. Goldfarb, S.
Lucidi, J. J. More, S. G. Nash, E. J. Polak, E.
Sachs, D. Shanno, Ph. L. Toint, V. J. Torczon, J.
P. Vial, R. J. Vanderbei, M. H. Wright, S.
Wright, J. Z. Zhang.
Proceedings: Edited proceedings including the
invited lectures and a selection of contributed
lectures will be published.
Further information can be found online
(http://www.dis.uniromal.it/or/erice) or
requested via email (erice@dis.uniromal.it).
Franco Giannessi, School Director
Gianni Di Pillo and Almerico Murli,
Workshop Directors
International Congress of
Mathematicians 2002
Beijing, China
August 2028, 2002
The next International Congress of
Mathematicians will take place in Beijing,
People's Republic of China, August 2028, 2002.
For preregistration forms and additional
information, please visit the web site
(http://www.icm2002.org.cn), or email the
organizers (icmsec@beijing.icm2002.org.cn).
Preregistrants will automatically receive progress
reports on the upcoming Congress as well as
final registration materials.
The International Congress of Mathematicians
takes place every four years and is supported by
the International Mathematical Union (IMU).
For information regarding the IMU, please visit
its web site (http://mathunion.org/).
All submissions will be refereed according to the
usual standards of Mathematical Programming.
Information about this issue can be obtained from
the guest editors for this volume (or at
www.math.washington.edu/thomas).
Karen Aardal, Department of Computer Science,
Utrecht University, P.O. Box 80089, 3508 TB Utrecht,
The Netherlands; email: aardal@cs.uu.nl and
Rekha R.Thomas, Department of Mathematics, Box
354350, University of Washington, Seattle, WA 98195
4350, U.S.A.; email:
thomas@math.washington.edu
SoPTIMA65
MAY 2001
imn/Ind apeler We invite OPTIMA readers to submit solutions to the problems to Robert Bosch
l (bobb@cs.oberlin.edu). The most attractive solutions will be presented in a forthcoming issue.
Painting by Numbers
Robert A. Bosch
September 28, 2000
A paintbynumbers puzzle consists of an m x n
grid of pixels (the canvas) together with m + n
clustersize sequences, one for each row and col
umn. The goal is to paint the canvas with a pic
ture that satisfies the following constraints:
1. Each pixel must be black or white.
2. If a row or column has clustersize
sequence sp, S2 . sk, then it must contain
k clusters of black pixels the first with s,
black pixels, the second with s2 black pixels,
and so on.
It should be noted that "first" means "leftmost"
for rows and "topmost" for columns, and that
rows and columns need not begin or end with
black pixels.
1 4 .....
S1 3 ...
2 .5 U . .. .. ... ..
2 5 . . ...
3 ........
2 5 .. ... ... ... ... ................
Figure 1
Small paintbynumbers puzzles are easy to
solve. One useful strategy is to alternate between
analyzing the row clustersize sequences and the
column clustersize sequences. For the puzzle
displayed in Figure 1, a first pass through the
rows enables us to paint 15 pixels black and one
pixel white, as displayed in Figure 2. (Since row
l's clustersize sequence is 3, 6, its first three pix
els must be black, its fourth pixel must be white,
and its final six pixels must be black. Since rows
7 and 8 have clustersize sequence 2, 5, their
second clusters must occupy pixels 4 through 8,
5 through 9, or 6 through 10; in each case they
occupy pixels 6, 7, and 8.) A subsequent pass
through the columns enables us to paint 19
more pixels black and 19 more pixels white, as
displayed in Figure 3.
Four more passes through the rows and
columns enable us to paint the remaining 46
pixels. The unique solution, a picture of a space
man, is displayed in Figure 4. It should be noted
that only well designed paintbynumbers puz
zles have unique solutions.
2 11 11 11111123
3 211 211123 8 9
113
2
3 3 ... ... ...
14
25
25
3 ... ... ... ...
Figure 2
1 1
1 1
14
113
2 5
2 5
3
Figure 3
11 1
1 1
2
14
25
Figure 4
An IP Formulation
In this section, we describe how integer pro
gramming can be used to solve paintbynum
bers puzzles. Our approach is to think of a
paintbynumbers puzzle as a problem com
prised of two interlocking tiling problems: one
involving the placement of the row clusters on
the canvas, and the other involving the place
ment of the column clusters. Note that if a pixel
is painted black, it must be covered by both a
row cluster and a column cluster; if it is painted
white, it must be left uncovered by the row clus
ters and column clusters.
Notation
Let
m = the number of rows,
n = the number of columns,
k = the number of clusters in row i,
k = the number of clusters in column,
s 1, s,' .. s'k = the clustersize sequence
for row i,
sI, s. 2 sk= the clustersize sequence
for column,
In addition, let
e, = the smallest value ofj such that row i's th
cluster can be placed in row i with its left
most pixel occupying pixelj,
1,= the largest value ofj such that row i's th
cluster can be placed in row i with its left
most pixel occupying pixelj
e = the smallest value of i such that column
j's th cluster can be placed in column
with its topmost pixel occupying pixel i,
1, = the largest value of i such that column's
th cluster can be placed in column with
its topmost pixel occupying pixel i.
(The letters "e" and "l" stand for "earliest" and
"latest.") In our example puzzle, row 7's first
cluster must be placed so that its leftmost pixel
occupies pixel 1, 2, or 3; row 7's second cluster
must be placed so that its leftmost pixel occupies
pixel 4, 5, or 6. In other words, e, = 1, 1, = 3,
e, = 4 and = 6
Variables
Our formulation uses three sets of variables.
One set specifies which pixels are painted black.
For each 1
z 1 if row i's jth pixel is painted black,
Z,'J 0 if row i's jth pixel is painted white
0SP TI M 6 5
PAGE 16
MAY 2001
The other two sets of variables are concerned
with the placements of the row and column
clusters. For each 1 < i < m, 1 < t < k and
e,
1 if row i's tth cluster is placed in row
with its leftmost pixel occupying
itj = pixelj,
0 if not.
For each 1
e ,< i < I' let
1 if column's tth cluster is placed in
column with its topmost pixel
jti occupying pixel i,
xjti =
0 if not.
Cluster Constraints
To ensure that row i's th cluster appears in row i
exactly once, we impose
S.I i.
j=e,,t
The sum is over all pixels j such that row i 's
th cluster can be placed in row i with its left
most pixel occupying pixelj. For the two clus
ters of row 7 of our example puzzle, we impose
y',7,1 + Y7,1,2 + Y7,1,3 = 1
and
Y7,2,4+ Y7,2,5 + 7,2,6= 1
To guarantee that row i's (t + 1)'h cluster is
placed to the right of its th cluster, we impose
1
i + Yt+
j'=j+s: ,+1l
for each e, +1 < j < I:
For row 7 of our example puzzle, we impose
Y7,1,2 Y7,2,5 + Y7,2,6
and
y7,1,3 < Y,2,6'
Similar constraints (involving 1.. 's) are
needed for the column clusters.
Double Coverage Constraints
To guarantee that each black pixel is covered by
both a row cluster and a column cluster, we
PAGE 17
impose, for each 1 < i < m and 1
following pair of inequalities:
kh
t=K
'maxin{e ,
/'=min{ ,js +1l
k max{Z ,},
z, l
=1 i'=min{e ,,is',+1
The first inequality states that if row i'sjth
pixel is painted black, then at least one of row i's
clusters must be placed in such a way that it
covers row i'sjth pixel. (The upper and lower
limits of the second summation make sure that
j' satisfies the inequalities e,
+ 1
i's h cluster can be placed in row i with its left
most pixel occupying pixel j'. The last two hold
if and only if row i'sjth pixel is covered when
row i's th cluster is placed in row i with its left
most pixel occupyingj') The other one makes
sure that if row i's jth pixel is painted black, then
at least one of column's clusters covers it. In
our example problem, for row 7's fifth pixel, we
impose
and 7,5 Y7,2,4 + Y7,2,5
and
Z7,5 x5,3,7 + '
Finally, we include constraints that prevent
white pixels from being covered by the row clus
ters or column clusters. For each 1 < i < m,
each 1
each' that satisfies the inequalities
j s + 1
we impose
z.. >i ,. I I71
112
And for each 1 < i
1 < t k and each i' that satisfies 4 2 3
the inequalities e < i' < I and 3 1 4
, 3313
i s + 1 < i i, we impose 214
29
zt> 211
27
In our example problem, for row 14
82
7's fifth pixel, we impose 6 2 2
2 8 13
Z, Z y7,2,4, 1 5 5 2
13 2 1
,s Y7,2,5 3 1 2 41
Z 5,3,7 11313
z, x 5,4,7 2 112
Objective Function
There is no need for an objective function here.
Problems
Interested readers may enjoy trying to solve the
following problems:
1. Try to improve the IP formulation. To
obtain our code for constructing Cplex
".lp" files for the IP formulation, point
your browser to www.oberlin.edu/math/
faculty/people/bosch.html
and scroll down to the Optima
Mindsharpener section.
2. Is there a better method for solving paint
bynumbers puzzles?
3. Solve the paintbynumbers puzzle dis
played in Figure 5. This puzzle was
designed by Won Yoon Jo. Its solution is
quite lovely, and it is perhaps the most dif
ficult 20 x 20 puzzle on the Paintby
Numbers Web Page (much more difficult
than many of the 30 x 30 and 40 x 40 puz
zles). To reach the PaintbyNumbers Web
Page, point your browser to www02.sonet.
ne.jp/kajitani/cgibin/pbn.cgi/index.html
3 1
11411 5 I 1 13
121213141111 1 161 1322332
1 11112616 7171119111 1 233 2
2 111 4 121 1 2212 1213352 226
Figure 5
SOPTIMA65
PAGE 18
dear the next issue
teadnvjievt r the next issue
feijit6 *" i
This is the last issue by the
cuilii en ed Icoil eami! \\
v\ ih co c[link Ill oIt [ce
a1tl[1or, [I.u l1 ,re con
[ibutecd \irih .alricles duling
[lie pat yea.ll .
W'e also \\ish to thank EIW.
* Wtabrjkefor her [anudc, ob
', 4f de.ignei of OPTINIMA.
She hasbeen the designed
inlce.4i 85, helping make i[
one of the most important
and recognized neCv.letters
in [lie mi chcmle ti.ical l sciences
coinnuniyr. The qualinr
emphlasi oft [e ,ocic v has
been capmteid in hel cie
ati\e giaplnc, and lIVout.
The_ NPS c\piceec it,
appieciation for her xce\l
lenxc IoIk aind v aiihe' oti
\\ell. A v1, 1m tlunk \oIu 1"0
all of their assistance is also
directed to Don Hearn, the
founding editor, and Patsy
Messinger.
Finally, we wish the incom
ing Editor Jens Clausen and
his team GOOD LUCK!
 S *.*a* '
Application for Membership
I wish to enroll as a member of the Society.
My subscription is for my personal use and not for the benefit of any library or institution.
O I will pay my membership dues on receipt of your invoice.
O I wish to pay by credit card (Master/Euro or Visa).
CREDIT CARD NO.
EXPIRATION DATE
FAMILYNAME
MAILING ADDRESS
TELEPHONE NO. TELEFAX NO.
EMAIL
SIGNATURE 0
Mail to:
Mathematical Programming Society
3600 University City Sciences Center
Philadelphia PA 191042688 USA
Cheques or money orders should be made
payable to The Mathematical Programming
Society, Inc. Dues for 2001, including sub
scription to the journal Mathematical
Programming, are US $75.
Student applications: Dues are onehalf the
above rate. Have a faculty member verify your
student status and send application with dues
to above address.
Faculty verifying status
Institution
v
F
i .illiaiii, ,11,,
Springer ad (insert film)
O P T I M A
MATHEMATICAL PROGRAMMING SOCIETY
UNIVERSITY OF
SFLORIDA
Center for Applied Optimization
371 Weil
PO Box 116595
Gainesville FL 326116595 USA
FIRST CLASS MAIL
EDITOR:
Karen Aardal
Department of Computer Science
Utrecht University
PO Box 80089
3508 TB Utrecht
The Netherlands
email: aardal@cs.ruu.nl
URL: http://www.cs.ruu.nl/staff/aardal.html
BOOK REVIEW EDITOR:
Robert Weismantel
Universitat Magdeburg
Fakultat fir Mathematik
Universitatsplatz 2
D39106 Magdeburg
Germany
email: weismant@math.unimagdeburg.de
Donald W. Hearn, FOUNDING EDITOR
Christina Loosli, DESIGNER
PUBLISHED BY THE
MATHEMATICAL PROGRAMMING SOCIETY &
GATOREngineering, PUBLICATION SERVICES
University of Florida
Journal contents are subject to change by the publisher
