Contact: iaindunning at gmail, @iaindunning, github/IainNZ
The PuzzlOR is a website connected to INFORMS that publishes operations research-related problems bimonthly. In a series of posts I’m going to solve some of the problems to demonstrate JuMP, an algebraic modeling language for optimization embedded in Julia.
(Link to full problem statement)
In this problem we have a 5x5 grid with two possibilities for each lot: residential or commercial. Let’s define a binary variable \( x_{ij} \) that is 1 if and only if lot \( (i,j) \) is residential. The challenge now is to express the objective. For each row and column there are the following possibilities:
We will need additional variables that don’t directly correspond to a decision in the original problem, commonly referred to as auxiliary variables, to model this problem. Before we talk about that let’s first review some integer programming modelling “tricks”.
In linear integer programming, all constraints must be linear in the decision variables. Solving non-linear integer programming problems is possible, but is significantly more difficult. The challenge of integer programming modeling then is to “linearize” all the logic of the problem. Consider the following constraints:
Because \( y \) is a binary variable, the effect of these two constraints is to force \( y = 1 \). Consider the following extension, where \( z_1, z_2, z_3 \) are all binary variables:
If any of the three \( z \) variables is 1, \( y \) is forced to equal one. This is equivalent to the logical statement y = z1 or z2 or z3. Another way to write this for binary \( z \) is
Let’s apply this to the problem at hand.
Define new variables
and equivalent \( y_{C,i} \) for columns. The objective function of our problem, for the positive parts at least, can be written as
Note that \( x \) does not appear in the objective, and that the optimizer will want to set these \( y \) variables to one if at all possible to drive the objective function value higher. Knowing that, we simply need to restrict these variables from being 1 if the condition is not met. The following constraints achieve this goal:
To understand why, consider the case where \( \sum_{j=1}^5 x_{ij} = 4 \). If we plug that in to the right-hand-side of each constraint we obtain:
When combined with the objective sense, the \( \leq \) is effectively an equality. Note also that if the sum of the \( x \) variables is zero on the right-hand-side, the problem is still feasible - that is, we don’t force the \( y \) variables to satisfy mutually exclusive constraints like
We can use a similar approach for the negative point cases, but the constraints are perhaps less intuitive. As before, define variables
The objective function “contribution” from these variables can be written as
The optimizer tries to drive these variables towards zero. Before we were restricting variables away from one, but now we will restrict them away from zero with the following constraints:
As before, let’s consider an example. Suppose \( \sum_{j=1}^5 x_{ij} = 1 \):
You can also check that all variables will set to 1 in the case where the sum is zero, and that if the sum is 3 or higher then all can be set to zero.
Apart from a constraint to set the total number of residential lots, we have everything we need. Let’s build the model in JuMP and Julia.
# Assuming a solver has been previously installed, e.g. Cbc
using JuMP
function SolveUrban()
m = Model()
# x is indexed by row and column
@defVar(m, 0 <= x[1:5,1:5] <= 1, Int)
# y is indexed by R or C, and the points
# JuMP allows indexing on arbitrary sets
rowcol = ["R","C"]
points = [+5,+4,+3,-3,-4,-5]
@defVar(m, 0 <= y[rowcol,points,i=1:5] <= 1, Int)
# Objective - combine the positive and negative parts
@setObjective(m, Max, sum{
3*(y["R", 3,i] + y["C", 3,i])
+ 1*(y["R", 4,i] + y["C", 4,i])
+ 1*(y["R", 5,i] + y["C", 5,i])
- 3*(y["R",-3,i] + y["C",-3,i])
- 1*(y["R",-4,i] + y["C",-4,i])
- 1*(y["R",-5,i] + y["C",-5,i]), i=1:5})
# Constrain the number of residential lots
@addConstraint(m, sum{x[i,j], i=1:5, j=1:5} == 12)
# Add the constraints that link the auxiliary y variables
# to the x variables
# Rows
for i = 1:5
@addConstraint(m, y["R", 5,i] <= 1/5*sum{x[i,j], j=1:5}) # sum = 5
@addConstraint(m, y["R", 4,i] <= 1/4*sum{x[i,j], j=1:5}) # sum = 4
@addConstraint(m, y["R", 3,i] <= 1/3*sum{x[i,j], j=1:5}) # sum = 3
@addConstraint(m, y["R",-3,i] >= 1-1/3*sum{x[i,j], j=1:5}) # sum = 2
@addConstraint(m, y["R",-4,i] >= 1-1/2*sum{x[i,j], j=1:5}) # sum = 1
@addConstraint(m, y["R",-5,i] >= 1-1/1*sum{x[i,j], j=1:5}) # sum = 0
end
# Columns
for j = 1:5
@addConstraint(m, y["C", 5,j] <= 1/5*sum{x[i,j], i=1:5}) # sum = 5
@addConstraint(m, y["C", 4,j] <= 1/4*sum{x[i,j], i=1:5}) # sum = 4
@addConstraint(m, y["C", 3,j] <= 1/3*sum{x[i,j], i=1:5}) # sum = 3
@addConstraint(m, y["C",-3,j] >= 1-1/3*sum{x[i,j], i=1:5}) # sum = 2
@addConstraint(m, y["C",-4,j] >= 1-1/2*sum{x[i,j], i=1:5}) # sum = 1
@addConstraint(m, y["C",-5,j] >= 1-1/1*sum{x[i,j], i=1:5}) # sum = 0
end
# Solve it with the default solver (CBC)
status = solve(m)
if status == :Infeasible
error("Solver couldn't find solution!")
end
# Print results
println("Best objective: $(round(getObjectiveValue(m)))")
println(getValue(x))
end
SolveUrban()
Try it yourself to see what the solution is. Do you believe the solution is correct, does it match your intuition?
While you can probably solve this problem just by inspection and trying a few ideas, its interesting to consider how the an integer programming solver would handle this problem. The deep details of this are way beyond the scope of this blog post, but there are two main issues:
Symmetry: the rows and columns are all effectively interchangeable. For example, if row 1 has three residential lots and row 2 has four, we could swap these rows without changing the objective. This can be an issue for MIP solvers as there is a chance the solution will “bounce” between many equivalent possibilities before reaching optimality, although many techniques and heuristics are built into these solvers to counteract this and even harness it.
Fractionality: MIP solvers actually solve a succession of linear problems where the constraints \( x \in {0,1} \) are relaxed to \( 0 \leq x \leq 1 \). The result of this is that when we initially solve the problem, many of \( y \) variables will be fractional. MIP solvers add cuts (constraints that push the values towards the integer values) and branch (solve two new problems with a fractional variables fixed to either 0 or 1) to eventually find an integer solution. They may also employ heuristics that attempt to construct an integer solution from the fractional solution.
The problem size here is so small that a solution is found almost instantly. However these are real issues for “real-world” problems and a lot of research is devoted to both improving solvers and studying the implications of different formulations.