# Solving a problem with mathematical programming

Every month, IBM Research publish an interesting puzzle on their Ponder This page. Last month 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) See my previous post for additional
background and references on operations research and optimization.

, and how
to solve it with Julia’s JuMP package.

## 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. The full problem statement is here, along with an
example on a small grid.

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.

## Solution

The grid scores are represented as a 20 × 20 array of hexadecimal numbers:

```
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
```

We can parse it easily with the DelimitedFiles module from Julia’s
standard library.

```
using DelimitedFiles
function readgrid(filename)
open(filename) do f
parse.(Int, readdlm(f, String); base=16) .- 128
end
end
= readgrid("grid.txt") grid
```

We now need to define the actual optimization problem. First, we load JuMP and a solver which supports MILP (for instance GLPK).

`using JuMP, GLPK`

Defining a model consists of three stages Check out the Quick Start Guide for more info.

:

- declare some variables, their types, and their bounds,
- add some constraints,
- specify an objective.

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)`

.

```
= size(grid, 1)
n = Model(GLPK.Optimizer)
model @variable(model, x[1:n, 1:n], Bin)
```

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.

```
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
```

Finally, the objective is to maximize the total of all rewards on explored cells:

`@objective(model, Max, sum(grid[i, j] * x[i, j] for i = 1:n, j = 1:n))`

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). 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)`

.

```
optimize!(model)
= Int(objective_value(model))
obj = Tuple.(findall(value.(x) .> 0))
indices = sort([(a-1, b-1) for (a, b) = indices]) indices
```

The resulting objective value is 1424, and the explored indices are

```
[(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)]
```

## 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 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 solvers page, we can use MPS or LP format to use CPLEX or Gurobi for instance. Luckily, JuMP (or more accurately MathOptInterface) supports these formats (among others).

`write_to_file(model, "rover.lp") # or "rover.mps"`

We can now upload this file to the NEOS Server, and sure enough, a few seconds later, we get Gurobi’s output:

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

We get the same solution!

## Code

My complete solution is available on GitHub.