Title: Optima
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00090046/00065
 Material Information
Title: Optima
Series Title: Optima
Physical Description: Serial
Language: English
Creator: Mathematical Programming Society, University of Florida
Publisher: Mathematical Programming Society, University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: May 2001
 Record Information
Bibliographic ID: UF00090046
Volume ID: VID00065
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.


This item has the following downloads:

optima65 ( PDF )

Full Text

Mathematical Programming Society Newsletter

MAY 2001

,. I pi

.. ..... ......
)-* it#4 -
i F

-- 1;
--I .
u'' h
i J':


The Seventeenth Symposium on Mathematical
Programming A Great Success
Karen Aardal

The Symposium was held at the Georgia
Institute of Technology, August 7-11, 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 semi-plenary 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, Jean-Philippe Vial. During the opening
session, the Beale Orchard-Hays, 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.







S 10 PA


MAY 2001


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 Orchard-F

The Beale Orchard-Hays
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



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

ard-Hays prize, new ideas in the field, from cutting planes to
on prize, and the branch-and-cut 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 non-expert can get a feel for the ideas
645-656. 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,509-city 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 of optimization software.
The authors have advanced the state of the art
to the point where solving 1,000-city instances is
routine. Indeed, in a recent test by Applegate (to
gather information on the Beardwood-Halton-
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
the way cutting-plane separation algorithms can
be viewed. Previously our main route to produc-
blem (TSP) has
ing facet-defining 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 problem-specific 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
state-of-the-art branch-and-cut 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 publicly-available
high-quality implementations of the Lin-
Kernighan and Chained-Lin-Kernighan heuris-
tics for finding near-optimal 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
An honorable mention for the Beale -
Orchard-Hays 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), 68-75.
William Gropp, Jorge More, Optimization
environments and the NEOS Server, in:



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-
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
self-concordance. 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.

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-
A balanced matrix is a 0-1 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
well-studied because packing and covering

Linear Programs with balanced constraint matri-
ces have integral solutions. Also, balanced matri-
ces are a natural generalization of bipartite
The authors give a polynomial time algorithm
for recognizing balanced matrices, settling a
long-standing open question. The algorithm is
based on a deep, beautiful decomposition theo-
rem whose arduous proof combines brilliant
insights and meticulous case-checking.
A second award went to Michel Goemans
and David Williamson for their paper
"Improved approximation algorithms for the
maximum-cut and satisfiability problems using
semi-definite programming," Journal of the
Association for Computing Machinery, 42
(1995), 1115-1145.
This paper gives a much better approximation
algorithm for the fundamental maximum-cut
problem in graphs and several other optimiza-
tion problems by an elegant application of semi-
definite programming. It is the first paper to
apply semi-definite 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 semi-definite
relaxation. The authors' techniques have since
found a host of other applications.
The prize jury: W. Cook, R. Graham, R.
Kannan (chair).




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
factor-2 approximation algorithm for the gener-
alized Steiner network problem.
The third finalist, and winner of the Tucker
prize, is Bertrand Guenin from Carnegie-Mellon
University. Dr. Guenin's paper gives a complete
characterization of graphs for which a natural
LP-relaxation 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
The prize jury: K. Anstreicher (chair), R.
Burkhard, D. den Hertog, D. Karger, J. Lee.


MAY 2001

Solution of a weighing


Cor A.J. Hurkens*
Gerhard J. Woegingert


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, NL-5600 MB Eindhoven, The Netherlands.
tgwoegi@opt.math.tu-graz.ac.at; Institut fiir Mathematik B, TU Graz, Steyrergasse 30, A-8010 Graz,
Austria. Supported by the START program Y43-MAT of the Austrian Ministry of Science.


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 +76-102
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
3 76 + 102
-1 12 + 43
-12 + 43
1 -12+ 43
-43 + 76
1 43 + 76
-7 + 43
7-12+ 43
-1 -3 + 43
-3 + 43
-1 + 43
1 +43
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 -3-7 +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


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) */


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
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 non-increasing
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 */



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

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

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 5-tuples (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 meta-language and then all the opti-
mization and all the searching are performed
automatically by the underlying software) can
compete against human specialization.


We thank Clive Tooth for stating this weighing
problem on his web page at



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 mid-1980s, 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 shared-memory model typ-
ified by the SGI Origin series and by computers
manufactured by Sun and Hewlett-Packard. 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 branch-and-bound tree, or different
candidates for the next iterate. Novel parallel
approaches were developed for global optimiza-
tion, network optimization, and direct-search
methods for nonlinear optimization. Activity
was particularly widespread in parallel branch-
and-bound 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
high-end 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 task-farming
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


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 master-worker
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

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, message-passing 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 glide-in, 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 master-worker
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 Jean-Pierre 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 user-defined 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 "bottom-level" inter-
face that allows it to be implemented in other
grid computing toolkits.

3 Integer Programming

Consider the linear mixed integer programming
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 branch-and-
bound procedure, in which the work of explor-
ing the branch-and-bound tree is distributed
among a fixed number of processors. These
approaches are differentiated by their use of vir-
tual-shared-memory vs. message-passing 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 "fault-tolerant
Condor PVM") is an enhanced branch-and-
bound procedure that utilizes (globally valid)
cuts, pseudocosts for branching, preprocessing at
nodes within the branch-and-bound tree, and
heuristics to identify integer feasible solutions
In FATCOP's master-worker algorithm, each
task consists of exploration of a subtree of the
branch-and-bound 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


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. (Depth-first 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.



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


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
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 well-known seymour prob-
aster process checkpoints" by
lem from the MIPLIB library of integer pro-
e current state of computation to
grams. This well-known integer-programming
;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 Four-Color 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 Four-Color
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


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 branch-and-bound
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-

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 matrix-form representation is as

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 second-largest 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 Jean-Pierre
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 branch-and-bound 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 Condor-based 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 Frank-Wolfe
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 Frank-Wolfe 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 branch-and-bound strate-
gy. At each node of the branch-and-bound 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
level-k 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 level-k 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


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 branch-and-
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 depth-first 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 "finish-up" 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 level-k 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 branch-and-bound tree, for a specific prob-
lem and specific choices of the algorithmic

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

5 Stochastic Programming

The two-stage stochastic linear programming
problem with recourse can be formulated as fol-

def N
min Q(x) = x+ piQ (x)
subject to Ax = b, x > 0,

Qi(x) = mingT y s.t. Wy = hi Tx, Yi 0.

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 two-stage problem with recourse repre-
sents a situation in which one set of decisions
(represented by the first-stage 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 second-stage sce-
nario is known. The objective function represents


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 second-stage 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
aQ(x)= c+ YpaQ (x).
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
trust-region 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
piecewise-linear model function M(x) satisfying
M(x) < Q (x). At the kth iteration, a candidate
iterate z is obtained by solving the following

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 trust-region 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 oo-norm, the sub-
problem can be formulated as a linear program.)
Most of the computational work in the algo-
rithm involves solution of the N second-stage
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
second-stage 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 others-its 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
second-stage 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 trust-region 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 private-line telecommunications
services. In this model, each of 86 parameters
representing demand can independently take on
3 to 7 values, giving a total of approximately


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
ATR-1 47 19 .35 130
ATR-10 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 L-shaped 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 trust-region
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 trust-region approaches were consider-
ably faster than the L-shaped 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 trust-region approach.
The real interest lies, however, not in solving
single sampled instances of SSN, but in obtain-
ing high-quality 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.)
ATR-1 25 67 .57 211
ATR-5 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 second-stage work can
be divided into a much larger number of chunks
than for SSN 125 chunks, rather than 10 -
the synchronous trust-region 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.


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, grid-enabled
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, higher-quality solutions can be found, and
improvements to sampling methodology can be
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. Data-intensive applications (from tomogra-
phy and data mining) represent a potentially
huge field of work, but these require a somewhat
different approach from the compute-intensive
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!


The metaneos project was funded by National
Science Foundation grant CDA-9726385 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 W-31-
For more information on the metaneos proj-
ect, see .
Information on Condor can be found at
and on Globus at

[1] I. Foster and C. Kesselman, edi-
tors. The Grid: Blueprintfor a
New Computing Infrastructure.
Morgan-Kaufmann, 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 4-18.
IEEE Press, 1998.
[5] J.-P Goux, J. T. Linderoth, andM.
Yoder. Metacomputing and the
master-worker paradigm. Preprint
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
master-worker 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 branch-and-
bound methods for mixed-integer
programming on the CM-5.
SIAMJournal on Optimization,
4(4):794-814, 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:19-43,
[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 RRR-40-2000, 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 99-11, Computer Sciences
Department, University of

Wisconsin at Madison, 1999. To
appear in Annals of Operations
[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:150-173,
[14] A. Marzetta and A. Brtinger. A
dynamic-programming 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 339-348.
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-
[18] S. Sen, R. D. Doverspike, and S.
Cosares. Network planning with
random demand. Telecommuni-
cations Systems, 3:11-30, 1994.
[19] J. R. Birge and E Louveaux.
Introduction to Stochastic
Programming. Springer Series in
Operations Research. Springer,
[20] J. M. Muley and A. Ruszczynski.
A new scenario decomposition
method for large scale stochastic
optimization. Operations Research,
43:477-490, 1995.

1 P Ii 5


MAY 2001



1st Cologne Twente Workshop on Graphs
and Combinatorial Optimization
June 6-8, 2001, University of Cologne,
Cologne, Germany.
URL: http://www.zaik.uni-koeln.de/AFS/

June 13-15, 2001, Utrecht,
The Netherlands.
URL: http://www.cs.uu.nl/events/ipco2001
email: ipco2001@cs.uu.nl

INFORMS International 2001
June 17-20, 2001, Maui, Hawaii, USA.
URL: http://www.informs.org/Conf/

SIAM Annual Meeting
July 9-13, 2001, San Diego,
California, USA.
URL: http://www.siam.org/meetings/an01/
Symposium on OR 2001
September 3-5, 2001, Gerhard-Mercator-
Universitdit, Duisburg, Germany.
URL: http://www.uni-duisburg.de/or2001/

INFORMS Annual Meeting
November 4-7, 2001, Miami Beach,
Florida, USA.
URL: http://www.informs.org/Conf/

9th International Conference on
Stochastic Programming
Berlin, Germany, August 25-31, 2001
(held under the auspices of the Committee on Stochastic Programming
of the Mathematical Programming Society)
Contact Information. Web site: http://www.mathematik.
hu-berlin.de/SPO1; E-mail: sp01@mathematik.hu-berlin.de;
Mailing Address: SP01, c/o Prof. W Roemisch, Institute of Mathematics,
Humboldt-University Berlin, 10099 Berlin, Germany

Support. Humboldt-Universitaet Berlin, Gerhard-Mercator-Universitaet
Duisburg, Weierstrass-Institut fuer Angewandte Analysis und Stochastik
(WIAS), Konrad-Zuse-Zentrum 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,
two-stage 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.
J-B Wets (USA), W.T Ziemba (Canada)

Invited Plenary Speakers. Hans Foellmer (Berlin), Suvrajeet Sen
(Tucson), Alexander Shapiro (Atlanta), Stein W Wallace (Trondheim),
Roger J-B 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. Groewe-Kuska, R.
Henrion, J. Kerger, H. Lange, A. Moeller, I. Nowak, M.C. Steinbach, I.

Location. Humboldt-University Berlin (Main Building), Unter den
Linden 6, 10117 Berlin

Schedule. Tutorials: August 25-26, 2001; International Conference:
August 27-31, 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



MAY 2001

EURO 2001: The XVIII-th Euro
Conference on Operations

Why Visit Euro 2001: Meeting OR friends in
convenient and well-equipped 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
When: July 9-11, 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 e-mail address (info@euro2001.org).
Fees: Early: Euro 300 for Non-Students, 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 e-mail address
Further Information: E-mail: 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

MING Series B

Issue on Algebraic and Geometric
Methods in Discrete Optimization


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!

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, derivative-free 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 non-traditional 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
Deadline for submission of full papers: April 30,
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 e-mail (erice@dis.uniromal.it).

Franco Giannessi, School Director
Gianni Di Pillo and Almerico Murli,
Workshop Directors

International Congress of
Mathematicians 2002

Beijing, China
August 20-28, 2002

The next International Congress of
Mathematicians will take place in Beijing,
People's Republic of China, August 20-28, 2002.
For pre-registration forms and additional
information, please visit the web site
(http://www.icm2002.org.cn), or e-mail the
organizers (icmsec@beijing.icm2002.org.cn).
Pre-registrants 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

Karen Aardal, Department of Computer Science,
Utrecht University, P.O. Box 80089, 3508 TB Utrecht,
The Netherlands; e-mail: aardal@cs.uu.nl and
Rekha R.Thomas, Department of Mathematics, Box
354350, University of Washington, Seattle, WA 98195-
4350, U.S.A.; e-mail:


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 paint-by-numbers puzzle consists of an m x n
grid of pixels (the canvas) together with m + n
cluster-size 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 cluster-size
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 paint-by-numbers puzzles are easy to
solve. One useful strategy is to alternate between
analyzing the row cluster-size sequences and the
column cluster-size 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 cluster-size 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 cluster-size 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 paint-by-numbers puz-
zles have unique solutions.

2 11 11 11111123
3 211 211123 8 9

3 3 ... ... ...

3 ... ... ... ...

Figure 2

1 1
1 1


2 5
2 5

Figure 3

11 1
1 1



Figure 4

An IP Formulation

In this section, we describe how integer pro-
gramming can be used to solve paint-by-num-
bers puzzles. Our approach is to think of a

paint-by-numbers 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.


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 cluster-size sequence
for row i,
sI, s. 2 sk= the cluster-size 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


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


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

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


impose, for each 1 < i < m and 1 following pair of inequalities:



'maxin{e ,

/'=min{ ,j-s +1l

k max{Z ,},
z, l
=1 i'=min{e ,,i-s',+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
and 7,5 Y7,2,4 + Y7,2,5
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
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
zt> 211
In our example problem, for row 14
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.


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/

and scroll down to the Optima
Mindsharpener section.

2. Is there a better method for solving paint-
by-numbers puzzles?

3. Solve the paint-by-numbers 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 Paint-by-
Numbers Web Page (much more difficult
than many of the 30 x 30 and 40 x 40 puz-
zles). To reach the Paint-by-Numbers Web
Page, point your browser to www02.so-net.

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



dear the next issue
teadnvjievt r the next issue
feiji-t6- *" -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

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





Mail to:
Mathematical Programming Society
3600 University City Sciences Center
Philadelphia PA 19104-2688 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 one-half the
above rate. Have a faculty member verify your
student status and send application with dues
to above address.

Faculty verifying status




i .illiaiii, ,11,,

Springer ad (insert film)




Center for Applied Optimization
371 Weil
PO Box 116595
Gainesville FL 32611-6595 USA


Karen Aardal
Department of Computer Science
Utrecht University
PO Box 80089
3508 TB Utrecht
The Netherlands
e-mail: aardal@cs.ruu.nl
URL: http://www.cs.ruu.nl/staff/aardal.html

Robert Weismantel
Universitat Magdeburg
Fakultat fir Mathematik
Universitatsplatz 2
D-39106 Magdeburg
e-mail: weismant@math.uni-magdeburg.de

Christina Loosli, DESIGNER
University of Florida

Journal contents are subject to change by the publisher

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs