# Modeling Examples¶

This chapter includes commented examples on modeling and solving optimization problems with Python-MIP.

## The 0/1 Knapsack Problem¶

As a first example, consider the solution of the 0/1 knapsack problem: given a set \(I\) of items, each one with a weight \(w_i\) and estimated profit \(p_i\), one wants to select a subset with maximum profit such that the summation of the weights of the selected items is less or equal to the knapsack capacity \(c\). Considering a set of decision binary variables \(x_i\) that receive value 1 if the \(i\)-th item is selected, or 0 if not, the resulting mathematical programming formulation is:

The following python code creates, optimizes and prints the optimal solution for the 0/1 knapsack problem

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | ```
from mip import Model, xsum, maximize, BINARY
p = [10, 13, 18, 31, 7, 15]
w = [11, 15, 20, 35, 10, 33]
c, I = 47, range(len(w))
m = Model('knapsack')
x = [m.add_var(var_type=BINARY) for i in I]
m.objective = maximize(xsum(p[i] * x[i] for i in I))
m += xsum(w[i] * x[i] for i in I) <= c
m.optimize()
selected = [i for i in I if x[i].x >= 0.99]
print('selected items: {}'.format(selected))
``` |

Line 3 imports the required classes and definitions from Python-MIP. Lines 5-8 define the problem data. Line 10 creates an empty maximization
problem `m`

with the (optional) name of “knapsack”. Line 12 adds the binary decision variables to model `m`

and stores their references in a list `x`

. Line 14 defines the objective function of this model and line 16 adds the capacity constraint. The model is optimized in line 18 and the solution, a list of the selected
items, is computed at line 20.

## The Traveling Salesman Problem¶

The traveling salesman problem (TSP) is one of the most studied combinatorial optimization problems, with the first computational studies dating back to the 50s [Dantz54], [Appleg06]. To to illustrate this problem, consider that you will spend some time in Belgium and wish to visit some of its main tourist attractions, depicted in the map bellow:

You want to find the shortest possible tour to visit all these places. More formally, considering \(n\) points \(V=\{0,\ldots,n-1\}\) and a distance matrix \(D_{n \times n}\) with elements \(c_{i,j} \in \mathbb{R}^+\), a solution consists in a set of exactly \(n\) (origin, destination) pairs indicating the itinerary of your trip, resulting in the following formulation:

The first two sets of constraints enforce that we leave and arrive only once at each point. The optimal solution for the problem including only these constraints could result in a solution with sub-tours, such as the one bellow.

To enforce the production of connected routes, additional variables \(y_{i} \geq 0\) are included in the model indicating the sequential order of each point in the produced route. Point zero is arbitrarily selected as the initial point and conditional constraints linking variables \(x_{i,j}\), \(y_{i}\) and \(y_{j}\) are created for all nodes except the the initial one to ensure that the selection of the arc \(x_{i,j}\) implies that \(y_{j}\geq y_{i}+1\).

The Python code to create, optimize and print the optimal route for the TSP is included bellow:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 | ```
from itertools import product
from sys import stdout as out
from mip import Model, xsum, minimize, BINARY
# names of places to visit
places = ['Antwerp', 'Bruges', 'C-Mine', 'Dinant', 'Ghent',
'Grand-Place de Bruxelles', 'Hasselt', 'Leuven',
'Mechelen', 'Mons', 'Montagne de Bueren', 'Namur',
'Remouchamps', 'Waterloo']
# distances in an upper triangular matrix
dists = [[83, 81, 113, 52, 42, 73, 44, 23, 91, 105, 90, 124, 57],
[161, 160, 39, 89, 151, 110, 90, 99, 177, 143, 193, 100],
[90, 125, 82, 13, 57, 71, 123, 38, 72, 59, 82],
[123, 77, 81, 71, 91, 72, 64, 24, 62, 63],
[51, 114, 72, 54, 69, 139, 105, 155, 62],
[70, 25, 22, 52, 90, 56, 105, 16],
[45, 61, 111, 36, 61, 57, 70],
[23, 71, 67, 48, 85, 29],
[74, 89, 69, 107, 36],
[117, 65, 125, 43],
[54, 22, 84],
[60, 44],
[97],
[]]
# number of nodes and list of vertices
n, V = len(dists), set(range(len(dists)))
# distances matrix
c = [[0 if i == j
else dists[i][j-i-1] if j > i
else dists[j][i-j-1]
for j in V] for i in V]
model = Model()
# binary variables indicating if arc (i,j) is used on the route or not
x = [[model.add_var(var_type=BINARY) for j in V] for i in V]
# continuous variable to prevent subtours: each city will have a
# different sequential id in the planned route except the first one
y = [model.add_var() for i in V]
# objective function: minimize the distance
model.objective = minimize(xsum(c[i][j]*x[i][j] for i in V for j in V))
# constraint : leave each city only once
for i in V:
model += xsum(x[i][j] for j in V - {i}) == 1
# constraint : enter each city only once
for i in V:
model += xsum(x[j][i] for j in V - {i}) == 1
# subtour elimination
for (i, j) in product(V - {0}, V - {0}):
if i != j:
model += y[i] - (n+1)*x[i][j] >= y[j]-n
# optimizing
model.optimize(max_seconds=30)
# checking if a solution was found
if model.num_solutions:
out.write('route with total distance %g found: %s'
% (model.objective_value, places[0]))
nc = 0
while True:
nc = [i for i in V if x[nc][i].x >= 0.99][0]
out.write(' -> %s' % places[nc])
if nc == 0:
break
out.write('\n')
``` |

In line 10 names of the places to visit are informed. In line 17 distances are informed in an upper triangular matrix. Line 33 stores the number of nodes and a list with nodes sequential ids starting from 0. In line 36 a full \(n \times n\) distance matrix is filled. Line 41 creates an empty MIP model. In line 44 all binary decision variables for the selection of arcs are created and their references are stored a \(n \times n\) matrix named `x`

. Differently from the \(x\) variables, \(y\) variables (line 48) are not required to be
binary or integral, they can be declared just as continuous variables, the default variable type. In this case, the parameter `var_type`

can be omitted from the `add_var`

call.

Line 51 sets the total traveled distance as objective function and lines 54-62 include the constraints. In line 66 we call the optimizer specifying a time limit of 30 seconds. This
will surely not be necessary for our Belgium example, which will be solved instantly, but may be important for larger problems: even though high quality solutions may be found very quickly by the MIP solver, the time required to *prove* that the current solution is optimal may be very
large. With a time limit, the search is truncated and the best solution found during the search is reported. In line 69 we check for the availability of a feasible solution. To repeatedly check for the next node in the route we check for the solution value (`.x`

attribute) of all variables of outgoing arcs of the current node in the route (line 73). The optimal solution for our trip has length 547 and is depicted bellow:

## n-Queens¶

In the \(n\)-queens puzzle \(n\) chess queens should to be placed in a board with \(n\times n\) cells in a way that no queen can attack another, i.e., there must be at most one queen per row, column and diagonal. This is a constraint satisfaction problem: any feasible solution is acceptable and no objective function is defined. The following binary programming formulation can be used to solve this problem:

The following code builds the previous model, solves it and prints the queen placements:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | ```
from sys import stdout
from mip import Model, xsum, BINARY
# number of queens
n = 40
queens = Model()
x = [[queens.add_var('x({},{})'.format(i, j), var_type=BINARY)
for j in range(n)] for i in range(n)]
# one per row
for i in range(n):
queens += xsum(x[i][j] for j in range(n)) == 1, 'row({})'.format(i)
# one per column
for j in range(n):
queens += xsum(x[i][j] for i in range(n)) == 1, 'col({})'.format(j)
# diagonal \
for p, k in enumerate(range(2 - n, n - 2 + 1)):
queens += xsum(x[i][j] for i in range(n) for j in range(n)
if i - j == k) <= 1, 'diag1({})'.format(p)
# diagonal /
for p, k in enumerate(range(3, n + n)):
queens += xsum(x[i][j] for i in range(n) for j in range(n)
if i + j == k) <= 1, 'diag2({})'.format(p)
queens.optimize()
if queens.num_solutions:
stdout.write('\n')
for i, v in enumerate(queens.vars):
stdout.write('O ' if v.x >= 0.99 else '. ')
if i % n == n-1:
stdout.write('\n')
``` |

## Frequency Assignment¶

The design of wireless networks, such as cell phone networks, involves assigning communication frequencies to devices. These communication frequencies can be separated into channels. The geographical area covered by a network can be divided into hexagonal cells, where each cell has a base station that covers a given area. Each cell requires a different number of channels, based on usage statistics and each cell has a set of neighbor cells, based on the geographical distances. The design of an efficient mobile network involves selecting subsets of channels for each cell, avoiding interference between calls in the same cell and in neighboring cells. Also, for economical reasons, the total bandwidth in use must be minimized, i.e., the total number of different channels used. One of the first real cases discussed in literature are the Philadelphia [Ande73] instances, with the structure depicted bellow:

Each cell has a demand with the required number of channels drawn at the center of the hexagon, and a sequential id at the top left corner. Also, in this example, each cell has a set of at most 6 adjacent neighboring cells (distance 1). The largest demand (8) occurs on cell 2. This cell has the following adjacent cells, with distance 1: (1, 6). The minimum distances between channels in the same cell in this example is 3 and channels in neighbor cells should differ by at least 2 units.

A generalization of this problem (not restricted to the hexagonal topology), is the Bandwidth Multicoloring Problem (BMCP), which has the following input data:

- \(N\):
set of cells, numbered from 1 to \(n\);

- \(r_i \in \mathbb{Z}^+\):
demand of cell \(i \in N\), i.e., the required number of channels;

- \(d_{i,j} \in \mathbb{Z}^+\):
minimum distance between channels assigned to nodes \(i\) and \(j\), \(d_{i,i}\) indicates the minimum distance between different channels allocated to the same cell.

Given an upper limit \(\overline{u}\) on the maximum number of channels \(U=\{1,\ldots,\overline{u}\}\) used, which can be obtained using a simple greedy heuristic, the BMPC can be formally stated as the combinatorial optimization problem of defining subsets of channels \(C_1, \ldots, C_n\) while minimizing the used bandwidth and avoiding interference:

This problem can be formulated as a mixed integer program with binary variables indicating the composition of the subsets: binary variables \(x_{(i,c)}\) indicate if for a given cell \(i\) channel \(c\) is selected (\(x_{(i,c)}=1\)) or not (\(x_{(i,c)}=0\)). The BMCP can be modeled with the following MIP formulation:

Follows the example of a solver for the BMCP using the previous MIP formulation:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | ```
from itertools import product
from mip import Model, xsum, minimize, BINARY
# number of channels per node
r = [3, 5, 8, 3, 6, 5, 7, 3]
# distance between channels in the same node (i, i) and in adjacent nodes
# 0 1 2 3 4 5 6 7
d = [[3, 2, 0, 0, 2, 2, 0, 0], # 0
[2, 3, 2, 0, 0, 2, 2, 0], # 1
[0, 2, 3, 0, 0, 0, 3, 0], # 2
[0, 0, 0, 3, 2, 0, 0, 2], # 3
[2, 0, 0, 2, 3, 2, 0, 0], # 4
[2, 2, 0, 0, 2, 3, 2, 0], # 5
[0, 2, 2, 0, 0, 2, 3, 0], # 6
[0, 0, 0, 2, 0, 0, 0, 3]] # 7
N = range(len(r))
# in complete applications this upper bound should be obtained from a feasible
# solution produced with some heuristic
U = range(sum(d[i][j] for (i, j) in product(N, N)) + sum(el for el in r))
m = Model()
x = [[m.add_var('x({},{})'.format(i, c), var_type=BINARY)
for c in U] for i in N]
z = m.add_var('z')
m.objective = minimize(z)
for i in N:
m += xsum(x[i][c] for c in U) == r[i]
for i, j, c1, c2 in product(N, N, U, U):
if i != j and c1 <= c2 < c1+d[i][j]:
m += x[i][c1] + x[j][c2] <= 1
for i, c1, c2 in product(N, U, U):
if c1 < c2 < c1+d[i][i]:
m += x[i][c1] + x[i][c2] <= 1
for i, c in product(N, U):
m += z >= (c+1)*x[i][c]
m.optimize(max_nodes=30)
if m.num_solutions:
for i in N:
print('Channels of node %d: %s' % (i, [c for c in U if x[i][c].x >=
0.99]))
``` |

## Resource Constrained Project Scheduling¶

The Resource-Constrained Project Scheduling Problem (RCPSP) is a combinatorial optimization problem that consists of finding a feasible scheduling for a set of \(n\) jobs subject to resource and precedence constraints. Each job has a processing time, a set of successors jobs and a required amount of different resources. Resources are scarce but are renewable at each time period. Precedence constraints between jobs mean that no jobs may start before all its predecessors are completed. The jobs must be scheduled non-preemptively, i.e., once started, their processing cannot be interrupted.

The RCPSP has the following input data:

- \(\mathcal{J}\)
jobs set

- \(\mathcal{R}\)
renewable resources set

- \(\mathcal{S}\)
set of precedences between jobs \((i,j) \in \mathcal{J} \times \mathcal{J}\)

- \(\mathcal{T}\)
planning horizon: set of possible processing times for jobs

- \(p_{j}\)
processing time of job \(j\)

- \(u_{(j,r)}\)
amount of resource \(r\) required for processing job \(j\)

- \(c_r\)
capacity of renewable resource \(r\)

In addition to the jobs that belong to the project, the set \(\mathcal{J}\) contains the jobs \(x_{0}\) and \(x_{n+1}\). These jobs are dummy jobs and represent the beginning of the planning and the end of the planning. The processing time for the dummy jobs is zero and does not consume resources.

A binary programming formulation was proposed by Pritsker et al. [PWW69]. In this formulation, decision variables \(x_{jt} = 1\) if job \(j\) is assigned a completion time at the end of time \(t\); otherwise, \(x_{jt} = 0\). All jobs must finish in a single instant of time without violating the relationships of precedence and amount of available resources. The model proposed by Pristker can be stated as follows:

An instance is shown below. The figure shows a graph where jobs \(\mathcal{J}\) are represented by nodes and precedence relations \(\mathcal{S}\) are represented by directed edges. Arc weights represent the time-consumption \(p_{j}\), while the information about resource consumption \(u_{(j,r)}\) is included next to the graph. This instance contains 10 jobs and 2 renewable resources (\(\mathcal{R}=\{r_{1}, r_{2}\}\)), where \(c_{1}\) = 6 and \(c_{2}\) = 8. The time horizon \(\mathcal{T}\) can be estimated by summing the duration of all jobs.

The Python code for creating the binary programming model, optimize it and print the optimal scheduling for RCPSP is included below:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 | ```
from itertools import product
from mip import Model, xsum, BINARY
p = [0, 3, 2, 5, 4, 2, 3, 4, 2, 4, 6, 0]
u = [[0, 0], [5, 1], [0, 4], [1, 4], [1, 3], [3, 2], [3, 1], [2, 4],
[4, 0], [5, 2], [2, 5], [0, 0]]
c = [6, 8]
S = [[0, 1], [0, 2], [0, 3], [1, 4], [1, 5], [2, 9], [2, 10], [3, 8], [4, 6],
[4, 7], [5, 9], [5, 10], [6, 8], [6, 9], [7, 8], [8, 11], [9, 11], [10, 11]]
(R, J, T) = (range(len(c)), range(len(p)), range(sum(p)))
model = Model()
x = [[model.add_var(name="x({},{})".format(j, t), var_type=BINARY) for t in T] for j in J]
model.objective = xsum(x[len(J) - 1][t] * t for t in T)
for j in J:
model += xsum(x[j][t] for t in T) == 1
for (r, t) in product(R, T):
model += (
xsum(u[j][r] * x[j][t2] for j in J for t2 in range(max(0, t - p[j] + 1), t + 1))
<= c[r])
for (j, s) in S:
model += xsum(t * x[s][t] - t * x[j][t] for t in T) >= p[j]
model.optimize()
print("Schedule: ")
for (j, t) in product(J, T):
if x[j][t].x >= 0.99:
print("({},{})".format(j, t))
print("Makespan = {}".format(model.objective_value))
``` |

The optimal solution is shown bellow, from the viewpoint of resource consumption:

## Job Shop Scheduling Problem¶

The Job Shop Scheduling Problem (JSSP) is an NP-hard problem defined by a set of jobs that must be executed by a set of machines in a specific order for each job. Each job has a defined execution time for each machine and a defined processing order of machines. Also, each job must use each machine only once. The machines can only execute a job at a time and once started, the machine cannot be interrupted until the completion of the assigned job. The objective is to minimize the makespan, i.e. the maximum completion time among all jobs.

For instance, suppose we have 3 machines and 3 jobs. The processing order for each job is as follows (the processing time of each job in each machine is between parenthesis):

Job \(j_1\): \(m_3\) (2) \(\rightarrow\) \(m_1\) (1) \(\rightarrow\) \(m_2\) (2)

Job \(j_2\): \(m_2\) (1) \(\rightarrow\) \(m_3\) (2) \(\rightarrow\) \(m_1\) (2)

Job \(j_3\): \(m_3\) (1) \(\rightarrow\) \(m_2\) (2) \(\rightarrow\) \(m_1\) (1)

Bellow there are two feasible schedules:

The first schedule shows a naive solution: jobs are processed in a sequence and machines stay idle quite often. The second solution is the optimal one, where jobs execute in parallel.

The JSSP has the following input data:

- \(\mathcal{J}\)
set of jobs, \(\mathcal{J} = \{1,...,n\}\),

- \(\mathcal{M}\)
set of machines, \(\mathcal{M} = \{1,...,m\}\),

- \(o^j_r\)
the machine that processes the \(r\)-th operation of job \(j\), the sequence without repetition \(O^j = (o^j_1,o^j_2,...,o^j_m)\) is the processing order of \(j\),

- \(p_{ij}\)
non-negative integer processing time of job \(j\) in machine \(i\).

A JSSP solution must respect the following constraints:

All jobs \(j\) must be executed following the sequence of machines given by \(O^j\),

Each machine can process only one job at a time,

Once a machine starts a job, it must be completed without interruptions.

The objective is to minimize the makespan, the end of the last job to be executed. The JSSP is NP-hard for any fixed \(n \ge 3\) and also for any fixed \(m \ge 3\).

The decision variables are defined by:

- \(x_{ij}\)
starting time of job \(j \in J\) on machine \(i \in M\)

- \(y_{ijk}=\)
\(\begin{cases} 1, & \text{if job } j \text{ precedes job } k \text{ on machine } i \text{,}\\ & i \in \mathcal{M} \text{, } j, k \in \mathcal{J} \text{, } j \neq k \\ 0, & \text{otherwise} \end{cases}\)

- \(C\)
variable for the makespan

Follows a MIP formulation [Mann60] for the JSSP. The objective function is computed in the auxiliary variable \(C\). The first set of constraints are the precedence constraints, that ensure that a job on a machine only starts after the processing of the previous machine concluded. The second and third set of disjunctive constraints ensure that only one job is processing at a given time in a given machine. The \(M\) constant must be large enough to ensure the correctness of these constraints. A valid (but weak) estimate for this value can be the summation of all processing times. The fourth set of constrains ensure that the makespan value is computed correctly and the last constraints indicate variable domains.

The following Python-MIP code creates the previous formulation, optimizes it and prints the optimal solution found:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 | ```
from itertools import product
from mip import Model, BINARY
n = m = 3
times = [[2, 1, 2],
[1, 2, 2],
[1, 2, 1]]
M = sum(times[i][j] for i in range(n) for j in range(m))
machines = [[2, 0, 1],
[1, 2, 0],
[2, 1, 0]]
model = Model('JSSP')
c = model.add_var(name="C")
x = [[model.add_var(name='x({},{})'.format(j+1, i+1))
for i in range(m)] for j in range(n)]
y = [[[model.add_var(var_type=BINARY, name='y({},{},{})'.format(j+1, k+1, i+1))
for i in range(m)] for k in range(n)] for j in range(n)]
model.objective = c
for (j, i) in product(range(n), range(1, m)):
model += x[j][machines[j][i]] - x[j][machines[j][i-1]] >= \
times[j][machines[j][i-1]]
for (j, k) in product(range(n), range(n)):
if k != j:
for i in range(m):
model += x[j][i] - x[k][i] + M*y[j][k][i] >= times[k][i]
model += -x[j][i] + x[k][i] - M*y[j][k][i] >= times[j][i] - M
for j in range(n):
model += c - x[j][machines[j][m - 1]] >= times[j][machines[j][m - 1]]
model.optimize()
print("Completion time: ", c.x)
for (j, i) in product(range(n), range(m)):
print("task %d starts on machine %d at time %g " % (j+1, i+1, x[j][i].x))
``` |

## Cutting Stock / One-dimensional Bin Packing Problem¶

The One-dimensional Cutting Stock Problem (also often referred to as One-dimensional Bin Packing Problem) is an NP-hard problem first studied by Kantorovich in 1939 [Kan60]. The problem consists of deciding how to cut a set of pieces out of a set of stock materials (paper rolls, metals, etc.) in a way that minimizes the number of stock materials used.

[Kan60] proposed an integer programming formulation for the problem, given below:

This formulation can be improved by including symmetry reducing constraints, such as:

The following Python-MIP code creates the formulation proposed by [Kan60], optimizes it and prints the optimal solution found.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | ```
from mip import Model, xsum, BINARY, INTEGER
n = 10 # maximum number of bars
L = 250 # bar length
m = 4 # number of requests
w = [187, 119, 74, 90] # size of each item
b = [1, 2, 2, 1] # demand for each item
# creating the model
model = Model()
x = {(i, j): model.add_var(obj=0, var_type=INTEGER, name="x[%d,%d]" % (i, j))
for i in range(m) for j in range(n)}
y = {j: model.add_var(obj=1, var_type=BINARY, name="y[%d]" % j)
for j in range(n)}
# constraints
for i in range(m):
model.add_constr(xsum(x[i, j] for j in range(n)) >= b[i])
for j in range(n):
model.add_constr(xsum(w[i] * x[i, j] for i in range(m)) <= L * y[j])
# additional constraints to reduce symmetry
for j in range(1, n):
model.add_constr(y[j - 1] >= y[j])
# optimizing the model
model.optimize()
# printing the solution
print('')
print('Objective value: {model.objective_value:.3}'.format(**locals()))
print('Solution: ', end='')
for v in model.vars:
if v.x > 1e-5:
print('{v.name} = {v.x}'.format(**locals()))
print(' ', end='')
``` |