234 lines
9.3 KiB
Org Mode
234 lines
9.3 KiB
Org Mode
---
|
||
title: "Solving a problem with mathematical programming"
|
||
date: 2021-04-02
|
||
tags: optimization, operations research, programming, julia
|
||
toc: false
|
||
---
|
||
|
||
Every month, IBM Research publish an interesting puzzle on their
|
||
[[https://www.research.ibm.com/haifa/ponderthis/index.shtml][Ponder This]] page. Last month [[https://www.research.ibm.com/haifa/ponderthis/challenges/March2021.html][puzzle]] was a nice optimization problem
|
||
about a rover exploring the surface of Mars.
|
||
|
||
In this post, I will explore how to formulate the problem as a
|
||
mixed-integer linear program (MILP)[fn:operations-research], and how
|
||
to solve it with Julia's [[https://jump.dev/][JuMP]] package.
|
||
|
||
[fn:operations-research] {-} See [[./operations-research-references.html ][my previous post]] for additional
|
||
background and references on operations research and optimization.
|
||
|
||
|
||
* The problem
|
||
|
||
The surface of Mars is represented as a $N \times N$ grid, where each
|
||
cell has a "score" (i.e. a reward for exploring the cell), and a
|
||
constant exploration cost of 128. The goal is to find the set of cell
|
||
which maximizes the total score. There is an additional constraint:
|
||
each cell can only be explored if all its upper neighbors were also
|
||
explored.[fn:problem]
|
||
|
||
This problem has a typical structure: we have to choose some variables
|
||
to maximize a specific quantity, subject to some constraints. Here,
|
||
JuMP will make it easy for us to formulate and solve this problem,
|
||
with minimal code.
|
||
|
||
[fn:problem] {-} The full problem statement is [[https://www.research.ibm.com/haifa/ponderthis/challenges/March2021.html][here]], along with an
|
||
example on a small grid.
|
||
|
||
|
||
* Solution
|
||
|
||
The grid scores are represented as a 20 × 20 array of hexadecimal
|
||
numbers:
|
||
#+begin_src txt
|
||
BC E6 56 29 99 95 AE 27 9F 89 88 8F BC B4 2A 71 44 7F AF 96
|
||
72 57 13 DD 08 44 9E A0 13 09 3F D5 AA 06 5E DB E1 EF 14 0B
|
||
42 B8 F3 8E 58 F0 FA 7F 7C BD FF AF DB D9 13 3E 5D D4 30 FB
|
||
60 CA B4 A1 73 E4 31 B5 B3 0C 85 DD 27 42 4F D0 11 09 28 39
|
||
1B 40 7C B1 01 79 52 53 65 65 BE 0F 4A 43 CD D7 A6 FE 7F 51
|
||
25 AB CC 20 F9 CC 7F 3B 4F 22 9C 72 F5 FE F9 BF A5 58 1F C7
|
||
EA B2 E4 F8 72 7B 80 A2 D7 C1 4F 46 D1 5E FA AB 12 40 82 7E
|
||
52 BF 4D 37 C6 5F 3D EF 56 11 D2 69 A4 02 0D 58 11 A7 9E 06
|
||
F6 B2 60 AF 83 08 4E 11 71 27 60 6F 9E 0A D3 19 20 F6 A3 40
|
||
B7 26 1B 3A 18 FE E3 3C FB DA 7E 78 CA 49 F3 FE 14 86 53 E9
|
||
1A 19 54 BD 1A 55 20 3B 59 42 8C 07 BA C5 27 A6 31 87 2A E2
|
||
36 82 E0 14 B6 09 C9 F5 57 5B 16 1A FA 1C 8A B2 DB F2 41 52
|
||
87 AC 9F CC 65 0A 4C 6F 87 FD 30 7D B4 FA CB 6D 03 64 CD 19
|
||
DC 22 FB B1 32 98 75 62 EF 1A 14 DC 5E 0A A2 ED 12 B5 CA C0
|
||
05 BE F3 1F CB B7 8A 8F 62 BA 11 12 A0 F6 79 FC 4D 97 74 4A
|
||
3C B9 0A 92 5E 8A DD A6 09 FF 68 82 F2 EE 9F 17 D2 D5 5C 72
|
||
76 CD 8D 05 61 BB 41 94 F9 FD 5C 72 71 21 54 3F 3B 32 E6 8F
|
||
45 3F 00 43 BB 07 1D 85 FC E2 24 CE 76 2C 96 40 10 FB 64 88
|
||
FB 89 D1 E3 81 0C E1 4C 37 B2 1D 60 40 D1 A5 2D 3B E4 85 87
|
||
E5 D7 05 D7 7D 9C C9 F5 70 0B 17 7B EF 18 83 46 79 0D 49 59
|
||
#+end_src
|
||
|
||
We can parse it easily with the [[https://docs.julialang.org/en/v1/stdlib/DelimitedFiles/][DelimitedFiles]] module from Julia's
|
||
standard library.
|
||
[fn::{-} [[file:../images/ponderthis_202103_grid.svg]]]
|
||
|
||
#+begin_src julia
|
||
using DelimitedFiles
|
||
|
||
function readgrid(filename)
|
||
open(filename) do f
|
||
parse.(Int, readdlm(f, String); base=16) .- 128
|
||
end
|
||
end
|
||
|
||
grid = readgrid("grid.txt")
|
||
#+end_src
|
||
|
||
We now need to define the actual optimization problem. First, we load
|
||
JuMP and a [[https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers][solver]] which supports MILP (for instance [[https://www.gnu.org/software/glpk/][GLPK]]).
|
||
|
||
#+begin_src julia
|
||
using JuMP, GLPK
|
||
#+end_src
|
||
|
||
Defining a model consists of three stages[fn:jump]:
|
||
- declare some variables, their types, and their bounds,
|
||
- add some constraints,
|
||
- specify an objective.
|
||
|
||
[fn:jump]{-} Check out the [[https://jump.dev/JuMP.jl/stable/quickstart/][Quick Start Guide]] for more info.
|
||
|
||
|
||
In our case, we have a single binary variable for each cell, which
|
||
will be 1 if the cell is explored by the rover and 0 otherwise. After
|
||
creating the model, we use the ~@variable~ macro to declare our
|
||
variable ~x~ of size ~(n, n)~.
|
||
|
||
#+begin_src julia
|
||
n = size(grid, 1)
|
||
model = Model(GLPK.Optimizer)
|
||
@variable(model, x[1:n, 1:n], Bin)
|
||
#+end_src
|
||
|
||
The "upper neighbors" of a cell ~(i, j)~ are ~[(i-1, j-1), (i-1, j),
|
||
(i-1, j+1)]~. Ensuring that a cell is explored only if all of its
|
||
upper neighbors are also explored means ensuring that ~x[i, j]~ is 1
|
||
only if it is also 1 for all the upper neighbors. We also have to
|
||
check that these neighbors are not outside the grid.
|
||
|
||
#+begin_src julia
|
||
for i = 2:n, j = 1:n
|
||
if j > 1
|
||
@constraint(model, x[i, j] <= x[i-1, j-1])
|
||
end
|
||
@constraint(model, x[i, j] <= x[i-1, j])
|
||
if j < n
|
||
@constraint(model, x[i, j] <= x[i-1, j+1])
|
||
end
|
||
end
|
||
#+end_src
|
||
|
||
Finally, the objective is to maximize the total of all rewards on explored cells:
|
||
|
||
#+begin_src julia
|
||
@objective(model, Max, sum(grid[i, j] * x[i, j] for i = 1:n, j = 1:n))
|
||
#+end_src
|
||
|
||
We now can send our model to the solver to be optimized. We retrieve
|
||
the objective value and the values of our variable ~x~, and do some
|
||
additional processing to get it in the expected format (0-based
|
||
indices while Julia uses 1-based indexing).[fn:check]
|
||
|
||
[fn:check]{-} In practice, you should also check that the solver
|
||
actually found an optimal solution, didn't find that the model is
|
||
infeasible, and did not run into numerical issues, using
|
||
~termination_status(model)~.
|
||
|
||
|
||
#+begin_src julia
|
||
optimize!(model)
|
||
obj = Int(objective_value(model))
|
||
indices = Tuple.(findall(value.(x) .> 0))
|
||
indices = sort([(a-1, b-1) for (a, b) = indices])
|
||
#+end_src
|
||
|
||
[fn::{-} [[file:../images/ponderthis_202103_explore.svg]]]
|
||
|
||
The resulting objective value is 1424, and the explored indices are
|
||
|
||
#+begin_src txt
|
||
[(0, 0), (0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9),
|
||
(0, 10), (0, 11), (0, 12), (0, 13), (0, 14), (0, 15), (0, 16), (0, 17), (0, 18),
|
||
(0, 19), (1, 0), (1, 1), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8),
|
||
(1, 9), (1, 10), (1, 11), (1, 12), (1, 13), (1, 14), (1, 15), (1, 16), (1, 17),
|
||
(1, 18), (1, 19), (2, 0), (2, 1), (2, 2), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7),
|
||
(2, 8), (2, 9), (2, 10), (2, 11), (2, 12), (2, 13), (2, 14), (2, 15), (2, 16),
|
||
(2, 17), (2, 18), (2, 19), (3, 0), (3, 1), (3, 2), (3, 3), (3, 4), (3, 5), (3, 6),
|
||
(3, 7), (3, 8), (3, 9), (3, 10), (3, 11), (3, 12), (3, 13), (3, 14), (3, 15),
|
||
(3, 16), (3, 17), (3, 18), (4, 0), (4, 1), (4, 2), (4, 3), (4, 4), (4, 5), (4, 6),
|
||
(4, 7), (4, 8), (4, 9), (4, 10), (4, 11), (4, 12), (4, 13), (4, 14), (4, 15),
|
||
(4, 16), (4, 17), (5, 0), (5, 1), (5, 2), (5, 3), (5, 4), (5, 5), (5, 6), (5, 7),
|
||
(5, 8), (5, 9), (5, 10), (5, 11), (5, 12), (5, 13), (5, 14), (5, 15), (5, 16),
|
||
(6, 0), (6, 1), (6, 2), (6, 3), (6, 4), (6, 5), (6, 6), (6, 7), (6, 8), (6, 9),
|
||
(6, 12), (6, 14), (6, 15), (7, 0), (7, 1), (7, 2), (7, 4), (7, 7), (8, 0), (8, 1),
|
||
(9, 0)]
|
||
#+end_src
|
||
|
||
* Exporting the model for external solvers
|
||
|
||
JuMP supports a wide variety of solvers, and this model is quite small
|
||
so open-source solvers are more than sufficient. However, let's see
|
||
how to use the [[https://neos-server.org/neos/][NEOS Server]] to give this problem to state-of-the-art
|
||
solvers!
|
||
|
||
Depending on the solver you plan to use, you will have to submit the
|
||
problem in a specific format. Looking at the [[https://neos-server.org/neos/solvers/index.html][solvers page]], we can use
|
||
[[https://www.gurobi.com/documentation/9.1/refman/mps_format.html][MPS]] or [[https://www.gurobi.com/documentation/9.1/refman/lp_format.html][LP]] format to use CPLEX or Gurobi for instance. Luckily, JuMP
|
||
(or more accurately [[https://github.com/jump-dev/MathOptInterface.jl][MathOptInterface]]) supports these formats (among
|
||
[[https://jump.dev/MathOptInterface.jl/stable/apireference/#File-Formats][others]]).
|
||
|
||
#+begin_src julia
|
||
write_to_file(model, "rover.lp") # or "rover.mps"
|
||
#+end_src
|
||
|
||
We can now upload this file to the NEOS Server, and sure enough, a few
|
||
seconds later, we get Gurobi's output:
|
||
|
||
#+begin_src txt
|
||
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (linux64)
|
||
Thread count: 32 physical cores, 64 logical processors, using up to 4 threads
|
||
Optimize a model with 1102 rows, 400 columns and 2204 nonzeros
|
||
Model fingerprint: 0x69169161
|
||
Variable types: 0 continuous, 400 integer (400 binary)
|
||
Coefficient statistics:
|
||
Matrix range [1e+00, 1e+00]
|
||
Objective range [1e+00, 1e+02]
|
||
Bounds range [1e+00, 1e+00]
|
||
RHS range [0e+00, 0e+00]
|
||
Found heuristic solution: objective 625.0000000
|
||
Presolve removed 116 rows and 45 columns
|
||
Presolve time: 0.01s
|
||
Presolved: 986 rows, 355 columns, 1972 nonzeros
|
||
Variable types: 0 continuous, 355 integer (355 binary)
|
||
|
||
Root relaxation: objective 1.424000e+03, 123 iterations, 0.00 seconds
|
||
|
||
Nodes | Current Node | Objective Bounds | Work
|
||
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
|
||
|
||
* 0 0 0 1424.0000000 1424.00000 0.00% - 0s
|
||
|
||
Explored 0 nodes (123 simplex iterations) in 0.01 seconds
|
||
Thread count was 4 (of 64 available processors)
|
||
|
||
Solution count 2: 1424 625
|
||
|
||
Optimal solution found (tolerance 1.00e-04)
|
||
Best objective 1.424000000000e+03, best bound 1.424000000000e+03, gap 0.0000%
|
||
|
||
********** Begin .sol file *************
|
||
|
||
# Solution for model obj
|
||
# Objective value = 1424
|
||
[...]
|
||
#+end_src
|
||
|
||
We get the same solution!
|
||
|
||
* Code
|
||
|
||
My complete solution is available [[https://github.com/dlozeve/ponder-this/blob/master/202103/rover.jl][on GitHub]].
|