Objectives
-
Understand how to solve integer linear optimization problems using Julia.
Section 5.5 Solving Integer LOPS with Julia
In this chapter, we will look at how to solve Integer Linear Optimization Problems (ILOPs) and Binary Linear Optimization Problems (BLOPs) using Julia.
We solved Problem 4.1.3 in Section 4.1 two different ways. The first involved cutting planes to trim the feasible set without eliminating integer solutions and the second was branch and bound, which successfully split the feasible set into pieces by eliminating non-integer solutions.
These were quite complex but important to understand how they work. However, when solving problems using software, we don’t need the details. As in the previous section, we will solve this using the
JuMP and HiGHS packages of Julia.
Subsection 5.5.1 Solving Integer Problems with Julia
Setting up an integer problem to solve with
JuMP and HiGHS is very similar to non-integer problems. Problem 4.1.3 can be written as
m = Model(HiGHS.Optimizer)
@variable(m, x₁ ≥ 0, Int)
@variable(m, x₂ ≥ 0, Int)
@constraint(m, 17x₁ + 32x₂ ≤ 136)
@constraint(m, 4x₁ + 4x₂ ≤ 25)
@objective(m, Max, 4x₁ + 5x₂)
print(m)
and note that when defining variables on lines 2 and 3, that we add
Int to the end indicating that the variables are integers. This returns
Max 4 x₁ + 5 x₂ Subject to 17 x₁ + 32 x₂ ≤ 136 4 x₁ + 4 x₂ ≤ 25 x₁ ≥ 0 x₂ ≥ 0 x₁ integer x₂ integer
We can solve this exactly the same way as before with
set_silent(m)
optimize!(m)
is_solved_and_feasible(m)
which returns
true so we have an optimal feasible solution. Then finally, we can find the values and the objective with
value(x₁), value(x₂)
which returns
(4.0, 2.0)
and
objective_value(m)
returns
24.0.
Checkpoint 5.5.1.
Solve the ILOP defined in Checkpoint 4.3.10 using the
JuMP and HiGHS packages in Julia. You can use the code above as a template.
Subsection 5.5.2 Solving Binary LOPs using Julia
Another common constraint on Linear Optimization problems is that variables may take on either 0 or 1. For example, the \(n\)-queens problem in Problem 2.1.10 is such an example.
In Subsubsection 2.1.2.1, the problem is written in mathematical form. First, if we look at this for \(n=4\text{,}\) the problem is a bit more reasonable.
First, the variable \(x_{i,j}\) will be a binary variable that will take on the value 1 if there is a queen on the \(i\)th row and \(j\)th column. To create the variables in Julia, we first create a model, then the variables, which we can do in bulk as the following:
m2 = Model(HiGHS.Optimizer)
n = 4
@variable(m2, x[1:n,1:n], Bin)
where the last term in the
@variable macro shows the variable is binary. Also, note that this make a array variable that will match the size of the board.
There are a number of ways to do set up the constraints including typing all of the constraints out. However, the following is a nice way to do this in a very general way. First, let’s do the row and column sums with the
mapslices command. Consider mapslices(sum, x, dims = [1]), which returns:
1×4 Matrix{AffExpr}:
x[1,1] + x[2,1] + x[3,1] + x[4,1] … x[1,4] + x[2,4] + x[3,4] + x[4,4]
Note that this is a row matrix. We then bound each of these by 1 with the following.
@constraint(m2, c1, mapslices(sum, x, dims = [1]) .<= ones(1,n))
@constraint(m2, c2, mapslices(sum, x, dims = [2]) .<= ones(n,1))
where the second line above is the column sums (note that the difference is the
dims = [2]).
To get the diagonal constraints, we will use the
diagm function, which produce a diagonal matrix. Consider diagm(1 => ones(3)) which returns
4×4 Matrix{Float64}:
0.0 1.0 0.0 0.0
0.0 0.0 1.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.0 0.0
and this will give us the diagonal constraint. The form of the
diagm function uses the argument of the form i = > v, where i is an integer that shows how far from the main diagonal the vector v should exist. If the integer is positive, then the diagonal will be on the upper right part of the matrix and negative values on the lower left. If we multiply this by x with diagm(1 => ones(3)).*x, the result is
4×4 Matrix{AffExpr}:
0 x[1,2] 0 0
0 0 x[2,3] 0
0 0 0 x[3,4]
0 0 0 0
And if this is summed, then we get the sum of the off diagonals. That is
sum(diagm(1 => ones(3)).*x) returns
\begin{equation*}
x_{1,2} + x_{2,3} + x_{3,4}
\end{equation*}
The diagonal constraints can be built with the following lines.
P = Matrix(I,n,n)[n:-1:1, :]
@constraint(m2, c3, [sum(diagm(i= > ones(n-abs(i))).*x) for i=-(n-2):(n-2)] .<= ones(2n-3))
@constraint(m2, c4, [sum(P*diagm(i= > ones(n-abs(i))).*x) for i=-(n-2):(n-2)] .<= ones(2n-3))
where the first line is a permutation matrix that will flip all of the rows of a matrix upon right multiplication. The structure
[. for i= -(n-2):(n-2)] is called a comprehension and builds a vector. See Section 5.2 for more examples.
Lastly, the objective needs to be defined and can be with
@objective(m2, Max, sum(x[:,:]))
and the 3rd argument sums over all elements of
x. If we put all of these commands in the same cell, then
m2 = Model(HiGHS.Optimizer)
n = 4
P = [ i+j == n+1 ? 1 : 0 for i=1:n, j=1:n]
@variable(m2, x[1:n,1:n], Bin)
@constraint(m2, c1, mapslices(sum, x, dims = [1]) .<= ones(1,n))
@constraint(m2, c2, mapslices(sum, x, dims = [2]) .<= ones(n,1))
@constraint(m2, c3, [sum(diagm(i=>ones(n-abs(i))).*x) for i=-(n-2):(n-2)] .<= ones(2n-3))
@constraint(m2, c4, [sum(P*diagm(i=>ones(n-abs(i))).*x) for i=-(n-2):(n-2)] .<= ones(2n-3))
@objective(m2, Max, sum([x[i,j] for i=1:n, j=1:n]))
Again, not necessary, but a good double check to do a
print(m2) which will print out the LOP as
Max x[1,1] + x[2,1] + x[3,1] + x[4,1] + x[1,2] + x[2,2] + x[3,2] + x[4,2] + x[1,3] + x[2,3] + x[3,3] + x[4,3] + x[1,4] + x[2,4] + x[3,4] + x[4,4]
Subject to
x[1,1] + x[2,1] + x[3,1] + x[4,1] ≤ 1
x[1,2] + x[2,2] + x[3,2] + x[4,2] ≤ 1
x[1,3] + x[2,3] + x[3,3] + x[4,3] ≤ 1
x[1,4] + x[2,4] + x[3,4] + x[4,4] ≤ 1
x[1,1] + x[1,2] + x[1,3] + x[1,4] ≤ 1
x[2,1] + x[2,2] + x[2,3] + x[2,4] ≤ 1
x[3,1] + x[3,2] + x[3,3] + x[3,4] ≤ 1
x[4,1] + x[4,2] + x[4,3] + x[4,4] ≤ 1
x[1,1] + x[2,2] + x[3,3] + x[4,4] ≤ 1
x[1,2] + x[2,3] + x[3,4] ≤ 1
x[1,3] + x[2,4] ≤ 1
x[2,1] + x[3,2] + x[4,3] ≤ 1
x[3,1] + x[4,2] ≤ 1
x[4,1] + x[3,2] + x[2,3] + x[1,4] ≤ 1
x[3,1] + x[2,2] + x[1,3] ≤ 1
x[2,1] + x[1,2] ≤ 1
x[4,2] + x[3,3] + x[2,4] ≤ 1
x[4,3] + x[3,4] ≤ 1
x[1,1] binary
x[2,1] binary
x[3,1] binary
x[4,1] binary
x[1,2] binary
x[2,2] binary
x[3,2] binary
x[4,2] binary
x[1,3] binary
x[2,3] binary
x[3,3] binary
x[4,3] binary
x[1,4] binary
x[2,4] binary
x[3,4] binary
x[4,4] binary
and if carefully look at all of the constraints that they are all that there is at most one queen (sum of the
x values along each row, column and diagonal), that all of the constraints are there.
Solve the problem with
set_silent(m2)
optimize!(m2)
is_solved_and_feasible(m2)
and again, since
true is returned that we have a feasible optimal solution. We can use the following shortcut (called broadcasting) to output the values with
value.(x)
where recall that the
. does broadcasting, that is, the function value is applied to all elements of x term by term. This returns
4×4 Matrix{Float64}:
-0.0 -0.0 1.0 0.0
1.0 0.0 0.0 -0.0
0.0 -0.0 -0.0 1.0
0.0 1.0 0.0 -0.0
Although we can read the solution from this, it seems odd that this is not an integer matrix. We can round this with
round.(Int, value.(x)) which returns
4×4 Matrix{Int64}:
0 0 1 0
1 0 0 0
0 0 0 1
0 1 0 0
which will show the locations of the queens. On a 4 by 4 chessboard, this would look like:
And in Figure 5.5.3, one can see that no queen can attack another one in a single move. Note that in general you get a single solution to a LOP when solving either via the simplex method or using software tools. There is another solution to this problem in which the board above can be flipped either vertically or horizontally. The software will not generate this other solution.
This problem can be scaled up to larger chessboards by careful adjusting the constraints above.
Subsection 5.5.3 Does the 3-queens problem have a solution?
It’s probably not too hard to see that if you have a 3 by 3 chess board, that you cannot get 3 queens on the board. Play with it for a few minutes and you’ll probably see this.
So what happens with the ILOP? The model can be created with:
q3 = Model(HiGHS.Optimizer)
@variable(q3, x[1:3,1:3], Bin)
for j=1:3
@constraint(q3, sum(x[:,j]) ≤ 1)
end
for i=1:3
@constraint(q3, sum(x[i,:]) ≤ 1)
end
@constraint(q3, sum( x[i,i] for i=1:3) ≤ 1)
@constraint(q3, sum( x[i,4-i] for i=1:3) ≤ 1)
@constraint(q3, sum( x[1,2]+x[2,1]) ≤ 1)
@constraint(q3, sum( x[3,2]+x[2,3]) ≤ 1)
@constraint(q3, sum( x[1,2]+x[2,3]) ≤ 1)
@constraint(q3, sum( x[3,2]+x[2,1]) ≤ 1)
@objective(q3, Max, sum(x[:,:]))
print(q3)
and this prints out the model:
Max x[1,1] + x[2,1] + x[3,1] + x[1,2] + x[2,2] + x[3,2] + x[1,3] + x[2,3] + x[3,3] Subject to x[1,1] + x[2,1] + x[3,1] ≤ 1 x[1,2] + x[2,2] + x[3,2] ≤ 1 x[1,3] + x[2,3] + x[3,3] ≤ 1 x[1,1] + x[1,2] + x[1,3] ≤ 1 x[2,1] + x[2,2] + x[2,3] ≤ 1 x[3,1] + x[3,2] + x[3,3] ≤ 1 x[1,1] + x[2,2] + x[3,3] ≤ 1 x[3,1] + x[2,2] + x[1,3] ≤ 1 x[2,1] + x[1,2] ≤ 1 x[3,2] + x[2,3] ≤ 1 x[1,2] + x[2,3] ≤ 1 x[2,1] + x[3,2] ≤ 1 x[1,1] binary x[2,1] binary x[3,1] binary x[1,2] binary x[2,2] binary x[3,2] binary x[1,3] binary x[2,3] binary x[3,3] binary
and the find the optimal solution with
set_silent(q3)
optimize!(q3)
is_solved_and_feasible(q3)
so
true is returned. You might think that this won’t have a solution, however, it seems to and if we print out the solution with
round.(Int, value.(x))
the result is
3×3 Matrix{Int64}:
0 0 0
1 0 0
0 0 1
This shows that the optimal solution has 2 queens on the board. Recall that the goal is to maximize the objective function with the given constraints. This means that putting a queen on an edge and an opposite corner is the largest number of queens allowed. You probably noticed this if you played with it a bit.
