diff --git a/images/ponderthis_202103_explore.svg b/images/ponderthis_202103_explore.svg
new file mode 100644
index 0000000..254cb43
--- /dev/null
+++ b/images/ponderthis_202103_explore.svg
@@ -0,0 +1,536 @@
+
+
diff --git a/images/ponderthis_202103_grid.svg b/images/ponderthis_202103_grid.svg
new file mode 100644
index 0000000..6565310
--- /dev/null
+++ b/images/ponderthis_202103_grid.svg
@@ -0,0 +1,585 @@
+
+
diff --git a/posts/ponder-this-2021-03.org b/posts/ponder-this-2021-03.org
new file mode 100644
index 0000000..acd2ba0
--- /dev/null
+++ b/posts/ponder-this-2021-03.org
@@ -0,0 +1,234 @@
+---
+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]].