19 Discrete Optimization
Previous chapters have focused on optimizing problems involving design variables that are continuous. Many problems, however, have design variables thatare naturally discrete, such as manufacturing problems involving mechanicalcomponents that come in fixed sizes or navigation problems involving choices between discrete paths. A discrete optimization problem has constraints such that thedesign variables must be chosen from a discrete set. Some discrete optimization
problems have infinite design spaces, and others are finite.1 Even for finite problems, where we could in theory enumerate every possible solution, it is generally
not computationally feasible to do so in practice. This chapter discusses both
exact and approximate approaches to solving discrete optimization problems
that avoid enumeration. Many of the methods covered earlier, such as simulated
annealing and genetic programming, can easily be adapted for discrete optimization problems, but we will focus this chapter on categories of techniques we have
not yet discussed.
Discrete optimization constrains the design to be integral. Consider the
problem:minimizex1 + x2
xExample 19.1. Discrete versionsof problems constrain the solution, often resulting in worse solutions than their continuous counterparts.

kxk ≤ 2x is integral√√
The optimum in the continuous case is x∗ = [− 2, − 2] with a value of√y = −2 2 ≈ −2.828. If x1 and x2 are constrained to be integer-valued, thenthe best we can do is to have y = −2 with x∗ ∈ {[−2, 0], [−1, −1], [0, −2]}.

2x2

subject toDiscrete optimization with a finite
design space is sometimes referred
to as combinatorial optimization. Fora review, see B. Korte and J. Vygen,
Combinatorial Optimization: Theory
and Algorithms, 5th ed. Springer,
2012.
1

0

−2
−2

0
x1

2

340 c ha p te r 19. d i sc re te op ti mi z a ti on

19.1 Integer Programs
An integer program is a linear program2 with integral constraints. By integral
constraints, we mean that the design variables must come from the set of integers.3
Integer programs are sometimes referred to as integer linear programs to emphasize
the assumption that the objective function and constraints are linear.
An integer program in standard form is expressed as:
minimizexc⊤ x

subject to Ax ≤ b
x≥0(19.1)

2

See chapter 11.

Integer programming is a very
mature field, with applicationsin operations research, communications networks, task scheduling, and other disciplines. Modern solvers, such as Gurobi and
CPLEX, can routinely handle problems with millions of variables.
There are packages for Julia thatprovide access to Gurobi, CPLEX,
and a variety of other solvers.
3x ∈ Znwhere Z n is the set of n-dimensional integral vectors.
Like linear programs, integer programs are often solved in equality form. Transforming an integer program to equality form often requires adding additionalslack variables s that do not need to be integral. Thus, the equality form for
integral programs is:
minimizec⊤ xx

subject toAx + s = bx≥0(19.2)s≥0x ∈ ZnMore generally, a mixed integer program (algorithm 19.1) includes both continuous and discrete design components. Such a problem, in equality form, isexpressed as:
minimizec⊤ xx

subject to

Ax = bx≥0(19.3)xD ∈ Z kDk

where D is a set of indices into the design variables that are constrained to bediscrete. Here, x = [xD , xC ], where xD represents the vector of discrete design
variables and xC the vector of continuous design variables.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 2. rou n d i n g 341

mutable struct MixedIntegerProgram
A
b
c
D
end

Algorithm 19.1. A mixed integer
linear program type that reflects
equation (19.3). Here, D is the setof design indices constrained to bediscrete.

19.2 Rounding

relax(MIP) = LinearProgram(MIP.A, MIP.b, MIP.c)function round_ip(MIP)x = minimize_lp(relax(MIP))
for i in MIP.D
x[i] = round(Int, x[i])endreturn xend

We can show that rounding the continuous solution for a constraint Ax ≤ bwhen A is integral is never too far from the optimal integral solution.4 If x∗c isan optimal solution of the LP with m × n matrix A, then there exists an optimal

Figure 19.1. Rounding can produce
an infeasible design point.

y

A common approach to discrete optimization is to relax the constraint that thedesign points must come from a discrete set. The advantage of this relaxation isthat we can use techniques, like gradient descent or linear programming, thattake advantage of the continuous nature of the objective function to direct the
search. After a continuous solution is found, the design variables are rounded to
the nearest feasible discrete design.
There are potential issues with rounding. Rounding might result in an infeasible design point, as shown in figure 19.1. Even if rounding results in a feasible
point, it may be far from optimal, as shown in figure 19.2. The addition of thediscrete constraint will typically worsen the objective value as illustrated in example 19.1. However, for some problems, we can show the relaxed solution is close
to the optimal discrete solution.
We can solve integer programs using rounding by removing the integer constraint, solving the corresponding linear program, or LP, and then rounding thesolution to the nearest integer. This method is implemented in algorithm 19.2.

x∗
x
Figure 19.2. The nearest feasible discrete design may be significantly worse than the best feasible
discrete design.
Algorithm 19.2. Methods for relaxing a mixed integer linear program
into a linear program and solving
a mixed integer linear program byrounding. Both methods accept amixed integer linear program MIP.
The solution obtained by rounding
may be suboptimal or infeasible.
4
W. Cook, A. M. Gerards, A. Schrijver, and É. Tardos, “Sensitivity Theorems in Integer Linear Programming,” Mathematical Programming,
vol. 34, no. 3, pp. 251–264, 1986.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

342 c ha p te r 19. d i sc re te op ti mi z a ti on

discrete solution x∗d with x∗c − x∗d ∞ less than or equal to n times the maximum
absolute value of the determinants of the submatrices of A.
The vector c need not be integral for an LP to have an optimal integral solution
because the feasible region is purely determined by A and b. Some approaches
use the dual formulation for the LP, which has a feasible region dependent on c,in which case having an integral c is also required.
In the special case of totally unimodular integer programs, where A, b, and c haveall integer entries and A is totally unimodular, the simplex algorithm is guaranteed
to return an integer solution. A matrix is totally unimodular if the determinant of
every submatrix5 is 0, 1, or −1, and the inverse of a totally unimodular matrix isalso integral. In fact, every vertex solution of a totally unimodular integer programis integral. Thus, every Ax = b for unimodular A and integral b has an integral
solution.
Several matrices and their total unimodularity are discussed in example 19.2.
Methods for determining whether a matrix or an integer linear program are totally
unimodular are given in algorithm 19.3.
Consider the following matrices:



1
0
1
0
1



0
0
0
 0
 0
1
0
1
0 −1

Example 19.2. Examples of totally
unimodular matrices.


1

0
0

The left matrix is not totally unimodular, as
1
1


−1 −1
0
0
0


0 −1 −1
0
 1
0
1
1
0 −1


1
= −2
−1

The other two matrices are totally unimodular.

19.3

A submatrix is a matrix obtainedby deleting rows and/or columnsof another matrix.
5

Cutting Planes

The cutting plane method is an exact method for solving mixed integer programs
when A is not totally unimodular.6 Modern practical methods for solving integer
programs use branch and cut algorithms7 that combine the cutting plane method

R. E. Gomory, “An Algorithm
for Integer Solutions to Linear
Programs,” Recent Advances inMathematical Programming, vol. 64,
pp. 269–302, 1963.

6

M. Padberg and G. Rinaldi, “A
Branch-and-Cut Algorithm for the
Resolution of Large-Scale Symmetric Traveling Salesman Problems,”
SIAM Review, vol. 33, no. 1, pp. 60–
100, 1991.
7© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 3. c u tti n g p l a n e s 343

isint(x, ϵ=1e-10) = abs(round(x) - x) ≤ ϵfunction is_totally_unimodular(A::Matrix)# all entries must be in [0,1,-1]
if any(a ∉ (0,-1,1) for a in A)
return falseend# brute force check every subdeterminantr,c = size(A)for i in 1 : min(r,c)
for a in subsets(1:r, i)for b in subsets(1:c, i)
B = A[a,b]if det(B) ∉ (0,-1,1)
return false
end
end
end
end
return true
end
function is_totally_unimodular(MIP)
return is_totally_unimodular(MIP.A) &&all(isint, MIP.b) && all(isint, MIP.c)
end

Algorithm 19.3. Methods for determining whether matrices A or
mixed integer programs MIP are
totally unimodular. The method
isint returns true if the given
value is integral.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

344 c ha p te r 19. d i sc re te op ti mi z a ti on

with the branch and bound method, discussed in the next section. The cutting
plane method works by solving the relaxed LP and then adding linear constraints
that result in an optimal solution.
We begin the cutting method with a solution x∗c to the relaxed LP, which must
be a vertex of Ax = b. If the D components in x∗c are integral, then it is also an
optimal solution to the original mixed integer program, and we are done. Aslong as the D components in x∗c are not integral, we find a hyperplane with x∗c on
one side and all feasible discrete solutions on the other. This cutting plane is anadditional linear constraint to exclude x∗c . The augmented LP is then solved for a
new x∗c .
Each iteration of algorithm 19.4 introduces cutting planes that make nonintegral components of x∗c infeasible while preserving the feasibility of the nearest
integral solutions and the rest of the feasible set. The integer program modified with these cutting plane constraints is solved for a new relaxed solution.
Figure 19.3 illustrates this process.

x2

Figure 19.3. The cutting plane
method introduces constraints until the solution to the LP is integral. The cutting plane is shown asa red line on the left. The feasible
region of the augmented LP is on
the right.

x2
x∗c

x∗c
x1

x1

We wish to add constraints that cut out nonintegral components of x∗c . For an
LP in equality form with constraint Ax = b, recall from section 11.2.1 that we canpartition a vertex solution x∗c to arrive at
AB x∗B + AV x∗V = b(19.4)

where x∗V = 0. The nonintegral components of x∗c will thus occur only in x∗B .
We can introduce an additional inequality constraint for each b ∈ B such thatxb∗ is nonintegral:8xb∗ − ⌊ xb∗ ⌋ − ∑ Ābv − Ābv xv ≤ 0(19.5)

Note that ⌊ x ⌋, or floor of x, rounds
x down to the nearest integer.

8v∈V© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 3. c u tti n g p l a n e s 345

frac(x) = modf(x)[1]function cutting_plane(MIP)
LP = relax(MIP)x, b_inds, v_inds = minimize_lp(LP)n_orig = length(x)D = copy(MIP.D)
while !all(isint(x[i]) for i in D)
AB, AV = LP.A[:,b_inds], LP.A[:,v_inds]
Abar = AB\AVb = 0
for i in D
if !isint(x[i])b += 1A2 = [LP.A zeros(size(LP.A,1));zeros(1,size(LP.A,2)+1)]A2[end,end] = 1
A2[end,v_inds] = (x->floor(x) - x).(Abar[b,:])b2 = vcat(LP.b, -frac(x[i]))c2 = vcat(LP.c, 0)
LP = LinearProgram(A2,b2,c2)
end
end
x, b_inds, v_inds = minimize_lp(LP)endreturn x[1:n_orig]
end

Algorithm 19.4. The cutting plane
method solves a given mixed integer program MIP and returns an
optimal design vector. An error isthrown if no feasible solution exists. The helper function frac returns the fractional part of a number, and the implementation forminimize_lp, algorithm 11.5, hasbeen adjusted to return the basic
and nonbasic indices b_inds andv_inds along with an optimal design x.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

346 c ha p te r 19. d i sc re te op ti mi z a ti on−1
where Ā = AB
AV . These cutting planes use only the V -components to cut out
the nonintegral components of x∗c .
Introducing a cutting plane constraint cuts out the relaxed solution x∗c , because
all xv are zero:xb∗ − ⌊ xb∗ ⌋ − ∑ Ābv − Ābv xv > 0(19.6)| {z } v∈V
>0
|
{z
}
0

A cutting plane is written in equality form using an additional integral slack
variable xk :
(19.7)
xk + ∑ Ābv − Ābv xv = ⌊ xb∗ ⌋ − xb∗v∈V

Each iteration of algorithm 19.4 thus increases the number of constraints andthe number of variables until solving the LP produces an integral solution. Only
the components corresponding to the original design variables are returned.
The cutting plane method is used to solve a simple integer linear program inexample 19.3.

19.4 Branch and Bound
One method for finding the global optimum of a discrete problem, such as aninteger program, is to enumerate all possible solutions. The branch and bound9
method guarantees that an optimal solution is found without having to evaluate
all possible solutions. Many commercial integer program solvers use ideas fromboth the cutting plane method and branch and bound. The method gets its namefrom the branch operation that partitions the solution space10 and the bound operationthat computes a lower bound for a partition.
Branch and bound is a general method that can be applied to many kinds of
discrete optimization problems, but we will focus here on how it can be used for
integer programming. Algorithm 19.5 provides an implementation that uses apriority queue, which is a data structure that associates priorities with elements
in a collection. We can add an element and its priority value to a priority queue
using the enqueue! operation. We can remove the element with the minimum
priority value using the dequeue! operation.

9
A. H. Land and A. G. Doig, “AnAutomatic Method of Solving Discrete Programming Problems,”
Econometrica, vol. 28, no. 3, pp. 497–
520, 1960.

The subsets are typically disjoint,
but this is not required. For branch
and bound to work, at least one
subset must have an optimal solution. D. A. Bader, W. E. Hart, and
C. A. Phillips, “Parallel Algorithm
Design for Branch and Bound,” inTutorials on Emerging Methodologies
and Applications in Operations Research, H. J. Greenberg, ed., Kluwer
Academic Press, 2004.
10© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 4 . bra n c h a n d bou n d 347

Example 19.3. The cutting plane
method used to solve an integer
program.

Consider the integer program:
2x1 + x2 + 3x3
#
#
"
2.5
−0.5
1.0
x=
−1.5
0.5 −1.5

minimizex

subject to

"

0.5
2.0

x ∈ Z3x≥0

The relaxed solution is x∗ ≈ [0.818, 0, 2.091], yielding:
"
#
"
#
"
#
0.5
1
−0.5
−0.091
AB =
AV =
Ā =
2 −1.50.5−0.455
From equation (19.7), the constraint for x1 with slack variable x4 is:
x4 + (⌊−0.091⌋ − (−0.091)) x2 = ⌊0.818⌋ − 0.818
x4 − 0.909x2 = −0.818

The constraint for x3 with slack variable x5 is:
x5 + (⌊−0.455⌋ − (−0.455)) x2 = ⌊2.091⌋ − 2.091
x5 − 0.545x2 = −0.091

The modified integer program has:


0.5
2

A=
0
0

−0.5
1
0.5
−1.5
−0.9090−0.545
0



0 0
0 0


1 0
0 1





2.5
−1.5 b=−0.818
−0.091

 
2
1
 
 
c = 3 00Solving the modified LP, we get x∗ ≈ [0.9, 0.9, 2.5, 0.0, 0.4]. Since this point
is not integral, we repeat the procedure with constraints:
x6 − 0.9x4 = −0.9

x8 − 0.5x4 = −0.5

x7 − 0.9x4 = −0.9

x9 − 0.4x4 = −0.4

and solve a third LP to obtain: x∗ = [1, 2, 3, 1, 1, 0, 0, 0, 0] with a final solutionof xi∗ = [1, 2, 3].
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

348 c ha p te r 19. d i sc re te op ti mi z a ti on

The algorithm begins with a priority queue containing a single LP relaxationof the original mixed integer program. Associated with that LP is a solution x∗cand objective value yc = c⊤ x∗c . The objective value serves as a lower bound on
the solution and thus serves as the LP’s priority in the priority queue. At each
iteration of the algorithm, we check whether the priority queue is empty. If it
is not empty, we dequeue the LP with the lowest priority value. If the solution
associated with that element has the necessary integral components, then wekeep track of whether it is the best integral solution found so far.
If the dequeued solution has one or more components in D that are nonintegral, we choose from x∗c such a component that is farthest from an integer value.
Suppose this component corresponds to the ith design variable. We branch by considering two new LPs, each one created by adding one of the following constraints
to the dequeued LP:11
 ∗
 ∗
(19.8)
orxi ≥ xi,cxi ≤ xi,cNote that ⌈ x ⌉, or ceiling of x,
rounds x up to the nearest integer.

11

as shown in figure 19.4. Example 19.4 demonstrates this process.
We compute the solution associated with these two LPs, which provide lower
bounds on the value of the original mixed integer program. If either solution
lowers the objective value when compared to the best integral solution seen sofar, it is placed into the priority queue. Not placing solutions already known to
be inferior to the best integral solution seen thus far allows branch and bound toprune the search space. The process continues until the priority queue is empty,
and we return the best feasible integral solution. Example 19.5 shows how branch
and bound can be applied to a small integer program.
Consider a relaxed solution x∗c = [3, 2.4, 1.2, 5.8] for an integer program withc = [−1, −2, −3, −4]. The lower bound is

Example 19.4. An example of asingle application of the branching
step in branch and bound.

y ≥ c⊤ x∗c = −34.6
We branch on a nonintegral coordinate of x∗c , typically the one farthest from∗ ,an integer value. In this case, we choose the first nonintegral coordinate, x2,cwhich is 0.4 from the nearest integer value. We then consider two new LPs,
one with x2 ≤ 2 as an additional constraint and the other with x2 ≥ 3 as anadditional constraint.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 4 . bra n c h a n d bou n d 349

x∗cFigure 19.4. Branching splits thefeasible set into subsets with anadditional integral inequality constraint.

x2

x1

x1 ≥ 2x1 ≤ 1x∗c,left
x2

x2

x1

x∗c,right

x1© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

350 c ha p te r 19. d i sc re te op ti mi z a ti on

function minimize_lp_and_y(LP)tryx = minimize_lp(LP)return (x, x⋅LP.c)
catch
return (fill(NaN, length(LP.c)), Inf)
end
end
function branch_and_bound(MIP)
LP = relax(MIP)x, y = minimize_lp_and_y(LP)n = length(x)x_best, y_best, Q = deepcopy(x), Inf, PriorityQueue()enqueue!(Q, (LP,x,y), y)while !isempty(Q)LP, x, y = dequeue!(Q)if any(isnan.(x)) || all(isint(x[i]) for i in MIP.D)
if y < y_best
x_best, y_best = x[1:n], yendelse
i = argmax([abs(x[i] - round(x[i])) for i in MIP.D])# x_i ≤ floor(x_i)
A, b, c = LP.A, LP.b, LP.c
A2 = [A zeros(size(A,1));[j==i for j in 1:size(A,2)]' 1]
b2, c2 = vcat(b, floor(x[i])), vcat(c, 0)
LP2 = LinearProgram(A2,b2,c2)
x2, y2 = minimize_lp_and_y(LP2)
if y2 ≤ y_best
enqueue!(Q, (LP2,x2,y2), y2)
end# x_i ≥ ceil(x_i)
A2 = [A zeros(size(A,1));[j==i for j in 1:size(A,2)]' -1]
b2, c2 = vcat(b, ceil(x[i])), vcat(c, 0)
LP2 = LinearProgram(A2,b2,c2)
x2, y2 = minimize_lp_and_y(LP2)
if y2 ≤ y_best
enqueue!(Q, (LP2,x2,y2), y2)
end
end
end
return x_best
end

Algorithm 19.5.
The branch
and
bound
algorithm
forsolving a mixted integer program MIP. The helper methodminimize_lp_and_y solves an LP
and returns both the solution
and its value. An infeasible LP
produces a NaN solution and an
Inf value. More sophisticatedimplementations will drop variables whose solutions are known
in order to speed computation.
The PriorityQueue type is provided by the DataStructures.jl
package.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 5 . d y n a mi c p rog ra mmi n g 351

We can use branch and bound to solve the integer program in example 19.3.
As before, the relaxed solution is x∗c = [0.818, 0, 2.09], with a value of 7.909.
We branch on the first component, resulting in two integer programs, onewith x1 ≤ 0 and one with x1 ≥ 1:
 




2
2.5
0.5 −0.5
1
0
1

 



cleft =  bleft = −1.5
Aleft =  2
0.5 −1.5 0
3
0
1
0
0
1
0


0.5

Aright =  2
1





−0.5
1
0

0.5 −1.5 0 
0
0
−1



2.5


bright = −1.5
1

 
2
1
 
cright =  30The left LP with x1 ≤ 0 is infeasible. The right LP with x1 ≥ 1 has arelaxed solution, x∗c = [1, 2, 3, 0], and a value of 13. We have thus obtained
our integral solution, xi∗ = [1, 2, 3].

19.5

Dynamic Programming

Dynamic programming12 is a technique that can be applied to problems with optimal
substructure and overlapping subproblems. A problem has optimal substructure if
an optimal solution can be constructed from optimal solutions of its subproblems.
Figure 19.5 shows an example.

ab

Example 19.5. Appying branch andbound to solve an integer programming problem.

c

12
The term dynamic programming
was chosen by Richard Bellman toreflect the time-varying aspect ofthe problems he applied it to andto avoid the sometimes negative
connotations words like research
and mathematics had. He wrote,‘‘I thought dynamic programming
was a good name. It was something
not even a Congressman could object to. So I used it as an umbrella
for my activities.’’ R. Bellman, Eyeof the Hurricane: An Autobiography.
World Scientific, 1984. p. 159.

Figure 19.5. Shortest path problems have optimal substructure because if the shortest path from any
a to c passes through b, then the
subpaths a → b and b → c are both
shortest paths.

A problem with overlapping subproblems solved recursively will encounterthe same subproblem many times. Instead of enumerating exponentially many potential solutions, dynamic programming either stores subproblem solutions, andthereby avoids having to recompute them, or recursively builds the optimal solu© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

352 c ha p te r 19. d i sc re te op ti mi z a ti ontion in a single pass. Problems with recurrence relations often have overlapping
subproblems. Figure 19.6 shows an example.
P9
P9

P4

P5
P3
P1

P2 P2

P7

P6

P7
P4
P1 P2

P4

P5

P3
P1 P1

P6

P0

P3

P2

Figure 19.6. We can computethe nth term of the Padovan
sequence, Pn = Pn−2 + Pn−3 , with
P0 = P1 = P2 = 1 by recursing
through all subterms (left). A
more efficient approach is tocompute subterms once and reuse
their values in subsequent calculations by exploiting the problem’s
overlapping substructure (right).

P0
P0

P1

Dynamic programming can be implemented either top-down or bottom-up, as
demonstrated in algorithm 19.6. The top-down approach begins with the desired
problem and recurses down to smaller and smaller subproblems. Subproblem
solutions are stored so that when we are given a new subproblem, we can either
retrieve the computed solution or solve and store it for future use.13 The bottomup approach starts by solving the smaller subproblems and uses their solutions
to obtain solutions to larger problems.
function padovan_topdown(n, P=Dict())
if !haskey(P, n)P[n] = n < 3 ? 1 :padovan_topdown(n-2,P) + padovan_topdown(n-3,P)
end
return P[n]endfunction padovan_bottomup(n)P = Dict(0=>1,1=>1,2=>1)
for i in 3 : n
P[i] = P[i-2] + P[i-3]end
return P[n]end

Storing subproblem solutions inthis manner is called memoization.

13

Algorithm 19.6. Computing the
Padovan sequence using dynamic
programming, with both the topdown and bottom-up approaches.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 5 . d y n a mi c p rog ra mmi n g 353

The knapsack problem is a well-known combinatorial optimization problem that
often arises in resource allocation.14 Suppose we are packing our knapsack for atrip, but we have limited space and want to pack the most valuable items. There
are several variations of the knapsack problem. In the 0-1 knapsack problem, we
have the following optimization problem:

The knapsack problem is an integer program with a single constraint, but it can be efficiently
solved using dynamic programming.

14n− ∑ vi xi

minimize
x

i =1

n

subject to∑ wi xi ≤ wmax

(19.9)

i =1

xi ∈ {0, 1} for all i in {1, . . . , n}

We have n items, with the ith item having integral weight wi > 0 and value vi .
The design vector x consists of binary values that indicate whether an item ispacked. The total weight cannot exceed our integral capacity wmax , and we seek
to maximize the total value of packed items.
There are 2n possible design vectors, which makes direct enumeration for large
n intractable. However, we can use dynamic programming. The 0-1 knapsackproblem has optimal substructure and overlapping subproblems. Consider having
solved knapsack problems with n items and several capacities up to wmax . The
solution to a larger knapsack problem with one additional item with weight wn+1
and capacity wmax will either include or not include the new item:• If it is not worth including the new item, the solution will have the same value
as a knapsack with n − 1 items and capacity wmax .
• If it is worth including the new item, the solution will have the value of aknapsack with n − 1 items and capacity wmax − wn+1 plus the value of the
new item.
The recurrence relation is:


0




knapsack(i − 1, wmax )knapsack(i, wmax ) =knapsack(i − 1, w


max )


max


knapsack(i − 1, wmax − wi ) + vi

if i = 0
if wi > wmax
(discard new item)(include new item)(19.10)otherwise© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

354 c ha p te r 19. d i sc re te op ti mi z a ti on

The 0-1 knapsack problem can be solved using the implementation in algorithm 19.7.
function knapsack(v, w, w_max)
n = length(v)y = Dict((0,j) => 0.0 for j in 0:w_max)
for i in 1 : n
for j in 0 : w_max
y[i,j] = w[i] > j ? y[i-1,j] :max(y[i-1,j],y[i-1,j-w[i]] + v[i])
end
end

end

Algorithm 19.7. A method for solving the 0-1 knapsack problem withitem values v, integral item weights
w, and integral capacity w_max. Recovering the design vector fromthe cached solutions requires additional iteration.

# recover solution
x, j = falses(n), w_max
for i in n: -1 : 1
if w[i] ≤ j && y[i,j] - y[i-1, j-w[i]] == v[i]# the ith element is in the knapsackx[i] = true
j -= w[i]endendreturn x19.6 Ant Colony Optimization
Ant colony optimization15 is a stochastic method for optimizing paths throughgraphs. This method was inspired by some ant species that wander randomly
in search of food, leaving pheromone trails as they go. Other ants that stumble
upon a pheromone trail are likely to start following it, thereby reinforcing thetrail’s scent. Pheromones slowly evaporate over time, causing unused trails to
fade. Short paths, with stronger pheromones, are traveled more often and thus
attract more ants. Thus, short paths create positive feedback that lead other ants
to follow and further reinforce the shorter path.
Basic shortest path problems, such as the shortest paths found by ants between
the ant hill and sources of food, can be efficiently solved using dynamic pro-

M. Dorigo, V. Maniezzo, and A.
Colorni, “Ant System: Optimization by a Colony of Cooperating
Agents,” IEEE Transactions on Systems, Man, and Cybernetics, Part B(Cybernetics), vol. 26, no. 1, pp. 29–
41, 1996.
15© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 6 . a n t c ol on y op ti mi z a ti on 355

gramming. Ant colony optimization has been used to find near-optimal solutions
to the traveling salesman problem, a much more difficult problem in which wewant to find the shortest path that passes through each node exactly once. Ant
colony optimization has also been used to route multiple vehicles, find optimal
locations for factories, and fold proteins.16 The algorithm is stochastic in natureand is thus resistant to changes to the problem over time, such as traffic delays
changing effective edge lengths in the graph or networking issues that remove
edges entirely.
Ants move stochastically based on the attractiveness of the edges available tothem. The attractiveness of transition i → j depends on the pheromone level andan optional prior factor:
A (i → j ) = τ (i → j )α η (i → j ) β(19.11)

where α and β are exponents for the pheromone level τ and prior factor η, respectively.17 For problems involving shortest paths, we can set the prior factor
to the inverse edge length ℓ(i → j) to encourage the traversal of shorter paths:
η (i → j) = 1/ℓ(i → j). A method for computing the edge attractiveness is given
in algorithm 19.8.
Suppose an ant is at node i and can transition to any of the nodes j ∈ J . The setof successor nodes J contains all valid outgoing neighbors.18 Sometimes edges
are excluded, such as in the traveling salesman problem where ants are prevented
from visiting the same node twice. It follows that J is dependent on both i andthe ant’s history.
function edge_attractiveness(graph, τ, η; α=1, β=5)A = Dict()
for i in 1 : nv(graph)neighbors = outneighbors(graph, i)
for j in neighborsv = τ[(i,j)]^α * η[(i,j)]^β
A[(i,j)] = v
end
end
return A
end

These and other applications
are discussed in these references:M. Manfrin, “Ant Colony Optimization for the Vehicle Routing
Problem,” Ph.D. dissertation,
Université Libre de Bruxelles,
2004. T. Stützle, “MAX-MIN Ant
System for Quadratic Assignment
Problems,” Technical University
Darmstadt, Tech. Rep., 1997.
and A. Shmygelska, R. AguirreHernández, and H. H. Hoos, “AnAnt Colony Algorithm for the 2D
HP Protein Folding Problem,”in International Workshop on Ant
Algorithms (ANTS), 2002.
16

Dorigo, Maniezzo, and Colorni
recommend α = 1 and β = 5.

17

The outgoing neighbors of a node
i are all nodes j such that i → j
is in the graph. In an undirected
graph, the neighbors and the outgoing neighbors are identical.

18

Algorithm 19.8. A method for computing the edge attractiveness table given graph graph, pheromone
levels τ, prior edge weights η,pheromone exponent α, and priorexponent β.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

356 c ha p te r 19. d i sc re te op ti mi z a ti on

The probability of edge transition i → j is:
P (i → j ) =A (i → j )∑ j′ ∈J A(i → j′ )(19.12)

Ants affect subsequent generations of ants by depositing pheromones. There
are several methods for modeling pheromone deposition. A common approach is
to deposit pheromones establishing after a complete path.19 Ants that do not finda path do not deposit pheromones. For shortest path problems, a successful ant
that has established a path of length ℓ deposits 1/ℓ pheromones on each edge it
traversed.
import StatsBase: Weights, sample
function run_ant(G, lengths, τ, A, x_best, y_best)
x = [1]
while length(x) < nv(G)i = x[end]neighbors = setdiff(outneighbors(G, i), x)
if isempty(neighbors) # ant got stuck
return (x_best, y_best)
end

end

end

as = [A[(i,j)] for j in neighbors]
push!(x, neighbors[sample(Weights(as))])M. Dorigo, G. Di Caro, and L. M.
Gambardella, “Ant Algorithms for
Discrete Optimization,” Artificial
Life, vol. 5, no. 2, pp. 137–172, 1999.
19

Algorithm 19.9. A method for simulating a single ant on a travelingsalesman problem in which the ant
starts at the first node and attempts
to visit each node exactly once.
Pheromone levels are increased at
the end of a successful tour. The
parameters are the graph G, edge
lengths lengths, pheromone levels
τ, edge attractiveness A, the best solution found thus far x_best, andits value y_best.

l = sum(lengths[(x[i-1],x[i])] for i in 2:length(x))
for i in 2 : length(x)τ[(x[i-1],x[i])] += 1/l
end
if l < y_best
return (x, l)
else
return (x_best, y_best)end

Ant colony optimization also models pheromone evaporation, which naturally
occurs in the real world. Modeling evaporation helps prevent the algorithm fromprematurely converging to a single, potentially suboptimal, solution. Pheromone
evaporation is executed at the end of each iteration after all ant simulations have© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 6 . a n t c ol on y op ti mi z a ti on 357

been completed. Evaporation decreases the pheromone level of each transition
by a factor of 1 − ρ, with ρ ∈ [0, 1].20
For m ants at iteration k, the effective pheromone update ismτ ( i → j ) ( k +1) = (1 − ρ ) τ ( i → j ) ( k ) + ∑a =1 ℓ1 ( a)(i → j ) ∈ P ( a)20

It is common to use ρ = 1/2.

(19.13)where ℓ( a) is the path length and P ( a) is the set of edges traversed by ant a.
Ant colony optimization is implemented in algorithm 19.10, with individual ant
simulations using algorithm 19.9. Figure 19.7 visualizes ant colony optimization
used to solve a traveling salesman problem.
function ant_colony_optimization(G, lengths;
m = 1000, k_max=100, α=1.0, β=5.0, ρ=0.5,
η = Dict((e.src,e.dst)=>1/lengths[(e.src,e.dst)]
for e in edges(G)))τ = Dict((e.src,e.dst)=>1.0 for e in edges(G))
x_best, y_best = [], Inf
for k in 1 : k_max
A = edge_attractiveness(G, τ, η, α=α, β=β)for (e,v) in ττ[e] = (1-ρ)*v
end
for ant in 1 : m
x_best,y_best = run_ant(G,lengths,τ,A,x_best,y_best)
end
end
return x_best
end

Algorithm 19.10.
Ant colony
optimization, which takes a directed or undirected graph G from
LightGraphs.jl and a dictionaryof edge tuples to path lengths
lengths. Ants start at the first nodein the graph. Optional parameters
include the number of ants per iteration m, the number of iterations
k_max, the pheromone exponent α,
the prior exponent β, the evaporation scalar ρ, and a dictionary ofprior edge weights η.

Figure 19.7. Ant colony optimization used to solve a travelingsalesman problem on a directed
graph using 50 ants per iteration.
Path lengths are the Euclidean distances. Color opacity corresponds
to pheromone level.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

358 c ha p te r 19. d i sc re te op ti mi z a ti on

19.7 Summary• Discrete optimization problems require that the design variables be chosenfrom discrete sets.
• Relaxation, in which the continuous version of the discrete problem is solved,
is by itself an unreliable technique for finding an optimal discrete solution but
is central to more sophisticated algorithms.
• Many combinatorial optimization problems can be framed as an integer program, which is a linear program with integer constraints.
• Both the cutting plane and branch and bound methods can be used to solve
integer programs efficiently and exactly. The branch and bound method isquite general and can be applied to a wide variety of discrete optimization
problems.
• Dynamic programming is a powerful technique that exploits optimal overlapping substructure in some problems.
• Ant colony optimization is a nature-inspired algorithm that can be used for
optimizing paths in graphs.

19.8 Exercises
Exercise 19.1. A Boolean satisfiability problem, often abbreviated SAT, requires determining whether a Boolean design exists that causes a Boolean-valued objectivefunction to output true. SAT problems were the first to be proven to belong to
the difficult class of NP-complete problems.21 This means that SAT is at least as
difficult as all other problems whose solutions can be verified in polynomial time.
Consider the Boolean objective function:

S. Cook, “The Complexity ofTheorem-Proving Procedures,” in
ACM Symposium on Theory of Computing, 1971.
21f (x) = x1 ∧ ( x2 ∨ ¬ x3 ) ∧ (¬ x1 ∨ ¬ x2 )Find an optimal design using enumeration. How many designs must be considered for an n-dimensional design vector in the worst case?
Exercise 19.2. Formulate the problem in exercise 19.1 as an integer linear program. Can any Boolean satisfiability problem be formulated as an integer linear
program?
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

19. 8 . e x e rc i se s 359

Exercise 19.3. Why are we interested in totally unimodular matrices? Furthermore,
why does every totally unimodular matrix contain only entries that are 0, 1, or−1?
Exercise 19.4. This chapter solved the 0-1 knapsack problem using dynamic
programming. Show how to apply branch and bound to the 0-1 knapsack problem, and use your approach to solve the knapsack problem with values v =[9, 4, 2, 3, 5, 3], and weights w = [7, 8, 4, 5, 9, 4] with capacity wmax = 20.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com