Solving a problem with mathematical programming

April 2, 2021

Tags: , , ,

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

grid = readgrid("grid.txt")

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 stagesCheck out the 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).

n = size(grid, 1)
model = Model(GLPK.Optimizer)
@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)
obj = Int(objective_value(model))
indices = Tuple.(findall(value.(x) .> 0))
indices = sort([(a-1, b-1) for (a, b) = 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.