Julia Community 🟣

Evan Wright
Evan Wright

Posted on

MILP with JuMP and HiGHS

In my previous post, I mentioned that the problem (Advent of Code 2018 day 23) can be reformulated as a mixed-integer linear program (MILP). In this post, we'll walk through a solution using JuMP.jl and HiGHS.jl. The formulation is based on this Reddit comment.

Input parsing is the same as last time. We set up the JuMP problem by defining variables x, y, and z as integers and the vector variable botsinrange as binary. This time, the value of botsinrange[i] is 1 if the point (x, y, z) is in range of bot i, otherwise 0.

using JuMP
using HiGHS

struct Nanobot
    x::Int
    y::Int
    z::Int
    radius::Int
end
function processinput(filename)
    puzin = [[parse(Int, m.match) for m in eachmatch(r"(\-?\d+)", line)]
             for line in readlines(filename)]
    bots = [Nanobot(l...) for l in puzin]
    return bots
end
bots = processinput("input201823.txt")

model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x, Int)
@variable(model, y, Int)
@variable(model, z, Int)
@variable(model, botsinrange[1:length(bots)], Bin)
Enter fullscreen mode Exit fullscreen mode

Next, define constraints of the form

abs(b.x - x) + abs(b.y - y) + abs(b.z - z) <= b.radius + (1 - botsinrange[i])*slack
Enter fullscreen mode Exit fullscreen mode

where i is the index of bot b. For the moment, assume slack is infinity. Then if we think bot i is out of range (botsinrange[i] = 0), the constraint holds (right hand side is Inf). Otherwise, if we think bot i is in range (botsinrange[i] = 1), the constraint holds if and only if abs(b.x - x) + abs(b.y - y) + abs(b.z - z) <= b.radius, i.e. the bot has to actually be in range. We will maximize the sum of botsinrange, so the optimizer will try to set botsinrange[i] = 1 for as many i as possible, and in those cases, it must also ensure the corresponding constraint holds.

So far so good, but there are two minor complications. First, the abs function is non-linear, so we can't have it in the constraint definition. However, note that if -x < c and x < c, then abs(x) < c. If x + y < c and -x + y < c and x - y < c and -x - y < c, then abs(x) + abs(y) < c. A similar argument holds for 3 variables, so we just need to create several linear constraints that imply the constraint defined above.

Second, we can't actually use Inf for the slack for numerical stability reasons; we'll have to use a large finite number. The slack can't be too low because we need the constraint to be relaxed if botsinrange[i] = 0. Consider the worst case, where the left hand side of the constraint is as large as possible. This may occur when (x, y, z) is at the "opposite" corner of 3d space from a bot i. The largest coordinate or radius in our puzzle input is

maxcoord = maximum(b -> max(abs(b.x), abs(b.y), abs(b.z), b.radius), bots)
Enter fullscreen mode Exit fullscreen mode

In the worst case, we are comparing a point like (maxcoord, maxcoord, maxcoord) to (-maxcoord, -maxcoord, -maxcoord), so the distance between those points is 2*3*maxcoord. We have additional distance for the bots' sensor radius, so we may use slack = 2*4*maxcoord.

Putting it all together, our constraints are

slack = 2*4*maximum(b -> max(abs(b.x), abs(b.y), abs(b.z), b.radius), bots)
for (i,b) in enumerate(bots)
    for absmult in Iterators.product((-1,1),(-1,1),(-1,1))
        @constraint(model, absmult[1]*(b.x - x) + absmult[2]*(b.y - y) +
                            absmult[3]*(b.z - z) <= b.radius + (1-botsinrange[i])*slack)
    end
end
Enter fullscreen mode Exit fullscreen mode

Finally, we can optimize the model.

@objective(model, Max, sum(botsinrange))
optimize!(model)
println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z)))
ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")
Enter fullscreen mode Exit fullscreen mode

For my input, we find the correct answer, but recall the problem statement asked if there are multiple points with maximum nanobot coverage, then we should choose the point closest to the origin. So, we can add a constraint that the objective function must be at least as high as we just found1

function coordisinrange(x, y, z, bot::Nanobot)
    return abs(bot.x - x) + abs(bot.y - y) + abs(bot.z - z) <= bot.radius
end
maxbots = count(b -> coordisinrange(round(Int, value(x)),
         round(Int, value(y)), round(Int, value(z)), b), bots)
@constraint(model, maxbots <= sum(botsinrange))
Enter fullscreen mode Exit fullscreen mode

and add a penalty term to our objective function to minimize abs(x) + abs(y) + abs(z).2 We have the same non-linear problem with abs as before, so we have to use a similar trick. The objective function becomes sum(botsinrange) - 0.00000001 * (xt + yt + zt) where the penalty is set by trial and error until the optimizer finds a solution. The coordinates are 8 digit numbers, so the penalty should be similar in scale to 1 unit of the objective---i.e. having one more bot in range.

@variable(model, xt, Int)
@variable(model, yt, Int)
@variable(model, zt, Int)
@constraint(model, x <= xt); @constraint(model, -x <= xt);
@constraint(model, y <= yt); @constraint(model, -y <= yt);
@constraint(model, z <= zt); @constraint(model, -z <= zt);

@constraint(model, maxbots <= sum(botsinrange))
@objective(model, Max, sum(botsinrange) - 0.00000001 * (xt + yt + zt))
optimize!(model)
println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z)))
ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")
Enter fullscreen mode Exit fullscreen mode

The solution is the same answer as before---which we know is correct thanks to Z3---but I'm curious if any other input needs the tie-breaking rule.

Discussion

The MILP solution feels a bit hacky compared to Z3, but it is faster at around 0.3-1.5 seconds without compilation, compared to 50+ seconds for Z3. It took trial and error to get the right slack, and even if the slack is high enough, the solver may provide an incorrect solution.3 For my input, slack of 1e11 results in an incorrect solution, while 1e10, 1e12-1e14 are ok.

Handling the L1 norm (abs) is annoying as well. I tried to use MathOptInterface.NormOneCone,

for (i,b) in enumerate(bots)
    @constraint(model2,
     [b.radius + (1-botsinrange[i])*Int(1e10); x - [b.x, b.y, b.z]] in MOI.NormOneCone(4))
end
Enter fullscreen mode Exit fullscreen mode

but there is a performance issue. With constraints from only 280 bots (the problem has 1000), the solver takes about 4 seconds, but with 290 bots it hangs (or takes longer than I care to wait).

In summary, MILP solvers can solve this problem more quickly than Z3 with a minor reformulation. There may be a better formulation that is less finicky. If you know of one, let me know.

Full Julia script

using JuMP
using HiGHS

struct Nanobot
    x::Int
    y::Int
    z::Int
    radius::Int
end

function processinput(filename)
    puzin = [[parse(Int, m.match) for m in eachmatch(r"(\-?\d+)", line)]
             for line in readlines(filename)]
    bots = [Nanobot(l...) for l in puzin]
    return bots
end

function coordisinrange(x, y, z, bot::Nanobot)
    return abs(bot.x - x) + abs(bot.y - y) + abs(bot.z - z) <= bot.radius
end

function main(filename)
    bots = processinput(filename)
    model = Model(HiGHS.Optimizer)
    slack = 2*4*maximum(b -> max(abs(b.x), abs(b.y), abs(b.z), b.radius), bots)
    set_silent(model)
    @variable(model, x, Int)
    @variable(model, y, Int)
    @variable(model, z, Int)
    @variable(model, botsinrange[1:length(bots)], Bin)
    for (i,b) in enumerate(bots)
        for absmult in Iterators.product((-1,1),(-1,1),(-1,1))
            @constraint(model, absmult[1]*(b.x - x) + absmult[2]*(b.y - y) +
                                absmult[3]*(b.z - z) <= b.radius + (1-botsinrange[i])*slack)
        end
    end
    @objective(model, Max, sum(botsinrange))
    optimize!(model)
    println("Step 1 solution:")
    println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z))),
     ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")

    # Step 2, with lower bound constraint for number of bots in range.
    maxbots = count(b -> coordisinrange(round(Int, value(x)),
         round(Int, value(y)), round(Int, value(z)), b), bots)
    @constraint(model, maxbots <= sum(botsinrange))

    @variable(model, xt, Int)
    @variable(model, yt, Int)
    @variable(model, zt, Int)
    @constraint(model, x <= xt); @constraint(model, -x <= xt);
    @constraint(model, y <= yt); @constraint(model, -y <= yt);
    @constraint(model, z <= zt); @constraint(model, -z <= zt);

    @objective(model, Max, sum(botsinrange) - 0.00000001 * (xt + yt + zt))
    optimize!(model)
    println("Step 2 solution:")
    println("""x=$(round(Int, value(x))), y=$(round(Int, value(y))), z=$(round(Int, value(z)))
     ans=$(round(Int, value(x)) + round(Int, value(y)) + round(Int, value(z)))""")

    return nothing
end

@time main("input201823.txt")
Enter fullscreen mode Exit fullscreen mode

Output (second run to avoid compilation time):

Step 1 solution:
x=44166763, y=43480550, z=38585775,
ans=126233088
Step 2 solution:
x=44166763, y=43480550, z=38585775
ans=126233088
  1.688667 seconds (590.15 k allocations: 32.233 MiB, 0.86% gc time)
Enter fullscreen mode Exit fullscreen mode

  1. You may be wondering why I didn't use round(Int, objective_value(model)). I found that for some values of slack, the objective_value is not the value of the objective at the solution (x, y, z). For example, the solution for step one results in an objective_value of 969 when the true number of bots in range of (x, y, z) is 972. 969 is still larger than the objective function evaluated at any other point, though. There may be a setting for HiGHS---that I'm unaware of---to resolve this issue. ↩

  2. We could do this in the first step as well, but I'm more confident in finding the global maxbots without the penalty in the first step, then trying to minimize the L1 norm of the solution. ↩

  3. If the slack is too low, then the problem is incorrectly specified. In the extreme case where slack is 0, the problem is to find a point within every bot's range, which isn't feasible. However, there shouldn't be an issue---mathematically---with a slack that is too large. ↩

Top comments (2)

Collapse
 
daniglez profile image
Daniel González Arribas

It's fun doing some AoC challenges with MILP. This year I did it for days 16 and 19, even if it's probably not the most efficient way.

Collapse
 
oheil profile image
oheil

Awesome!