---
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]].
