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