Which JuMP.jl solver for this problem?

[mthelm85]

I’ve been struggling for a couple of days now to coerce JuMP to solve a particular problem. The goal is fairly simple and I’ll describe a simplified version of my real problem here, along with the code:

Given the latitudes, longitudes, and populations of 10 cities, I need to divide the cities into 3 groups in such a way that the geographical distance between all cities in a group is minimized, and I also need the total populations of each of the three groups to be as close as possible.

I’ve tried this with GLPK, Ipopt, and Juniper and I can’t get any of them to work so I’m hoping someone with more knowledge/experience than me can tell me which solver is actually capable of handling this problem.

using DataFrames
using Distances
using GLPK
using Ipopt
using JuMP
using Juniper
using LinearAlgebra
using Random
using Statistics

# DataFrame with 10 cities including their lat/lon and population
cities = DataFrame(
    city=[
        "Marysville",
        "Perris",
        "Cleveland",
        "Worcester",
        "Columbia",
        "Waterbury",
        "Eagan",
        "Southfield",
        "Lafayette",
        "Boise"
    ],
    lat=[
        48.0517637,
        33.7825194,
        41.49932,
        42.2625932,
        34.0007104,
        41.5581525,
        44.8041322,
        42.4733688,
        30.2240897,
        43.6187102
    ],
    lon=[
        -122.1770818,
        -117.2286478,
        -81.6943605,
        -71.8022934,
        -81.0348144,
        -73.0514965,
        -93.1668858,
        -83.2218731,
        -92.0198427,
        -116.2146068
    ],
    pop=[
        63269,
        72326,
        390113,
        182544,
        133358,
        109676,
        65453,
        73006,
        124276,
        214237
    ]
)

num_groups = 3

# Compute the pairwise distances between each city
dm = Distances.pairwise(Haversine(6372.8), Matrix(cities[:, [2,3]])', dims=2)

# Assign random group numbers to start
groups = rand(1:num_groups, size(cities,1))

# Filter the distance matrix by group, compute the norm, then take the sum of each group's norm
ϵ  = sum([norm(dm[groups .== i, groups .== i]) for i in 1:num_groups])

# Take the norm of each group's population - total population / num_groups (target pop per group)
σ = norm([sum(cities.pop[groups .== i]) for i in 1:num_groups] .- sum(cities.pop) / num_groups)

My approach thus far has been to create decision variables for each of the 10 cities @variable(model, x[1:10], Int) and restrict their value to the integers 1, 2, and 3. The objective function is to minimize σ + ϵ. The only additional constraint I have is that I must indeed end up with 3 groups but I’m not attempting to restrict the minimum number of cities per group at this point in time.

Will one of the available open-source solvers be able to handle this problem? If not, is there a better way that I can reformulate it?

I’m thinking about creating 3 vectors of length 10 and each value in the vector would be binary and then maybe there’s some way to solve the problem like that?

[joaquimg]

Your problem seems a more complex version of the p-median problem:

http://courses.ieor.berkeley.edu/ieor151/lecture_notes/ieor151_lec13.pdf

You can probably start from there and try to add you additional constraints.

Possibly a MIP solver (you can see a list on jump manual) will be enough, if you can restrict yourself to linear constraints and continuous+integer+binary variables.

[mthelm85]

I was able to figure out a very simple way to have it divide the cities into groups with fairly equal populations:

target_pop = sum(cities.pop) / num_groups

model = Model(Cbc.Optimizer)

@variable(model, x[1:num_groups, 1:size(cities, 1)], Bin)

for i in 1:size(cities,1)
    @constraint(model, sum(x[:,i]) == 1)
end

for i in 1:num_groups
    @constraint(model, sum(x[i, :]) >= 3)
    @constraint(model, (x * cities.pop)[i] - target_pop <= 45_000)
    @constraint(model, (x * cities.pop)[i] - target_pop >= -45_000)
end

optimize!(model)

cities.group = zeros(size(cities,1))

for i in 1:3, j in 1:10
    if round(value.(x)[i,j]) == 1.0
        cities.group[j] = i
    end
end

for group in groupby(cities, :group)
    println(sum(group.pop))
end

I just had to play with the threshold values in the constraints by decreasing the absolute value until it was no longer feasible. I still cannot figure out how to make it group cities such that the cumulative distances between cities in a group are minimized…

It seems that I am always going to need the values from x to be able to retrieve the appropriate elements in my distance matrix but there’s no way to set the problem up like that as far as I can tell…

[mtanneau]

@mthelm85

the geographical distance between all cities in a group is minimized, and I also need the total populations of each of the three groups to be as close as possible

This sounds like a bi-objective problem, unless you put one criterion as a constraint, e.g., all cities in a given group should have a population that differ from at most XXX, or optimize a weighted combination of the two.

That being said, if you only want to solve the problem for 10 cities, brute force might be simple enough.
There are 1024 possible groups of cities, which would give at most ~10^9 combinations to check (much less with a bit of smart elimination). That should take only a few seconds on a decent machine.

Obviously, it does not scale, but it’s only a few lines of code and it will give you all Pareto-optimal combinations.

[mthelm85]

Unfortunately, this is just a toy example to play with. My real problem involves grouping over 3,000 U.S. counties

[joaquimg]

You can probably do a multiplication trick instead of using x to access a matrix, just multiply x by the matrix entries so that it will be the matrix entry when x is one and zero otherwise.
The most straightforward way one would think it to do, for each group
sum( d[i,j] * x[i, group] * x[j, group] ) == sum_dist[group]
However you can’t have products of binaries (if you want performance).
Since they are binaries there are a bunch of formulations that can be used, if you problem is tough you will have to dig into it.
If you want something simple you can:
Replace x[i, group] * x[j, group] by z[i,j,group], also binary.
and add constraints like:
z[I,j,group] <= x[i, group] , so that z is zero if x[i] is zero (same for x[j])
and
-2 +x[i]+x[j] <= z[i,j] - 1 , so that z is one if both are one.

I might have been a bit sloppy… hope you can get the idea.

There is a lot of Art involved in writing good models, and better models lead to results in less time although their optimal solution is the same.

JuMP is here to give us enough flexibility to write all this crazy stuff.

[joaquimg]

I you can get something up and running consider a PR to https://github.com/jump-dev/JuMPTutorials.jl

[mthelm85]

Unfortunately, when I tried something like this before it didn’t throw any errors but it doesn’t work:

dm = Distances.pairwise(Haversine(6372.8), Matrix(cities[:, [2,3]])', dims=2)
n = size(cities,1)

model = Model(Cbc.Optimizer)

@variable(model, x[1:num_groups, 1:n], Bin)

for i in 1:n
    @constraint(model, sum(x[:,i]) == 1) # a city cannot be in more than 1 group
end

for i in 1:num_groups
    @constraint(model, sum(x[i, :]) >= 3) # at least 3 cities per group
end

julia> @objective(model, Min, sum(sum(dm[x[i,:] .== 1, x[i,:] .== 1]) for i in 1:num_groups))
0.0

optimize!(model)

julia> objective_value(model)
0.0

I don’t have any idea what it’s doing but the objective value is always just 0.0 and it’s not actually optimizing the problem.

[mtanneau]

dm[x[i,:] .== 1, x[i,:] .== 1])

You can’t do that.
Here x is an array of VariableRef, not of numbers, so x[i, :] .== 1 yields an array of falses, and the resulting sum is 0.

[mthelm85]

@mtanneau

Here x is an array of VariableRef , not of numbers, so x[i, :] .== 1 yields an array of false s, and the resulting sum is 0 .

Ah, yes. I’m struggling a bit to understand Joaquim’s suggestion but I’ll keep working at it

[mtanneau]

Finding feasible solutions is not hard (any partition of the cities in 3 groups will do).
The hard part is getting the objective right.

First, you need to find a well-defined expression for the objective:

@mthelm85

the geographical distance between all cities in a group is minimized

Does it mean:

• The average/median/maximum/total pair-wise distances between cities, among all pairs of cities in a given group?
• The average/median/maximum/total distance between two cities in a given group?
• The average/median/maximum/total distance between a city and the center of mass of the group?

@mthelm85

the total populations of each of the three groups to be as close as possible

Does it mean:

• The difference between the most and least populous groups?
• The average difference between each group’s population and the average group population?

Then, you still have the two-objective issue.
If you don’t want to compute a full Parto front, a common approach is to take a weighted sum of the two. Obviously, it requires you to choose suitable coefficients.

Once this part is sorted out, things will follow much more easily.

[mthelm85]

@mtanneau

Does it mean:

• The average/median/maximum/total pair-wise distances between cities, among all pairs of cities in a given group?
• The average/median/maximum/total distance between two cities in a given group?
• The average/median/maximum/total distance between a city and the center of mass of the group?

I think the 3rd option would result in the best outcome for my purpose but it seems to me that it’s also the hardest so I’ve just been using the sum of the pairwise distance matrix for now.

@mtanneau

Does it mean:

• The difference between the most and least populous groups?
• The average difference between each group’s population and the average group population?

I would like the standard deviation between each group’s population to be minimized but for the sake of simplicity I’m simply adding a constraint that the difference between each group’s population and the result of the quotient of the total population (of all 10 cities) and the number of groups must be between -45,000 and 45,000 (this seems to work really well).

I’ve actually separated the problem into two for now and I’m satisfied with the population portion, I just can’t figure out how to obtain a result for the distance portion.Once I figure out the distance piece, I’ll look at how to include both objectives

[mtanneau]

For the distance part, if you can express the cities’ locations in Cartesian coordinates (not sure if that’s already the case), you could use any clustering method (e.g., K-means) to get decent solutions.

It’s not exact, but it should be simple enough to code, and it’s reasonable to expect it would yield decent clusters.

[mthelm85]

I’m getting a little bit closer but I can’t figure out how to add a constraint so that you end up with the desired number of groups…the code that is commented out would have served that purpose but it has the same problem that a VariableRef cannot be passed to the count function.

num_groups = 3

n = size(cities,1)

model = Model(GLPK.Optimizer)

@variable(model, x[1:n, 1:n, 1:num_groups], Bin)

for i in 1:n, j in 1:n, k in 1:num_groups
    @constraint(model, x[i,j,k] == x[j,i,k])
end

for i in 1:n, j in 1:n
    @constraint(model, x[i,j,1] + x[i,j,2] + x[i,j,3] == 1)
end

# for k in 1:num_groups
#     @constraint(model, count(x -> x == 0, sum(x[:,j,k] for j in 1:n)) == 3)
# end

@objective(model, Min, sum(dot(x[:,:,i], dm) for i in 1:num_groups))

optimize!(model)

Basically, I need to be able to add a constraint such that each of the 3 10 x 10 matrices contains 3 columns of zeros and 3 rows of zeros…so a valid x[:,:,1] might look like this:

5×5 Array{Int64,2}:
 0  0  0  0  0
 0  1  0  1  0
 0  0  0  0  0
 0  1  0  1  0
 0  0  0  0  0

In this group, we would have cities 2 and 4.

[odow]

@mthelm85

@constraint(model, count(x -> x == 0, sum(x[:,j,k] for j in 1:n)) == 3)

Something like:

@constraint(model, sum(1 - x[:,j,k] for j in 1:n) == 3)

Mixed-integer linear programming uses a very specific subset of the things you can do (namely, linear expressions). You can’t throw arbitrary Julia functions at it.

[mthelm85]

@odow

Mixed-integer linear programming uses a very specific subset of the things you can do (namely, linear expressions). You can’t throw arbitrary Julia functions at it.

Yes, I’m still learning - slowly, but surely.

Mathematical optimization is very interesting and very satisfying but as I’ve been learning I’ve been thinking it’s a bit like tricking a child into doing something that they don’t want to do. I have young children and if you tell them to clean their room they don’t want to do it but if I tell them we’re going to play a game to see who can pick up the most toys in the shortest amount of time, they’re all for it.

Writing these JuMP models is a bit like that - you can’t just express the problem bluntly and let the model go to work. You have to trick it into solving the problem in a roundabout way This is the art that @joaquimg was speaking of.

[mtanneau]

Start with a binary variable x_{i, k} that takes value 1 if city i is in group k and 0 otherwise.
Each city must be in a group, so we add the constraint \sum_{k} x_{i, k} = 1 for every i.

The total population of group k is then Q_{k} = \sum_{i} x_{i, k} q_{i}.
Adding bounds Q_{min} \leq Q_{k} \leq Q_{max} (as you suggested initially) will ensure that groups have similar population levels.

Then, there’s the question of the distance. We want to minimize the sum of distances between pairs of cities in a same group.
Let z_{i, j} be a binary variable that takes value 1 if cities i and j are in the same group, and 0 otherwise. Then the total distance is just \sum_{i, j} d_{i, j} z_{i, j}. (note that I’m counting each pair twice here).

To ensure that z_{i, j} = 1 if and only if cities i and j are in the same group, one way is to add the constraints z_{i, j} \geq x_{i, k} + x_{j, k} - 1 for every pair i, j and every k.
It is easy to verify that, if cities i, j are in group k, then z_{i, j} >= 1+ 1 - 1 = 1.
Otherwise, x_{i, k} + x_{j, k} cannot exceed 1, thus z_{i, j} \geq 0 and, since we are minimizing and z_{i, j} has positive objective coefficient, we get z_{i, j} = 0 in any optimal solution.

The size of this model is thus:

• N \times k + k + N \times N variables
• N + N \times N \times k constraints

I tried this with 10 cities, CPLEX solved the model in a fraction of second.
For 30 cities, it took CPLEX 10s.

I let you imagine what happens with 3,000 cities…

[mthelm85]

This is absolutely brilliant - thank you very much . I thought for a while about out how to link x and z and wasn’t able to see this simple relationship:

@mtanneau

one way is to add the constraints z_{i,j}≥x_{i,k}+x_{j,k}−1 for every pair i,j and every k .

I’m currently running it for about 100 cities with the Cbc solver and it’s taking quite a long time, but that’s okay - I’m just thrilled to finally have something that works after working at it for several days

[mthelm85]

@mtanneau, I’d like to give back by taking @joaquimg’s suggestion and writing a tutorial for JuMPTutorials.jl so what’s the best way to credit you? The tutorials there all have an “Originally contributed by: author name” line at the top. You basically wrote the tutorial in Post #17

[mtanneau]

Putting my name as one of the contributors (you can find it in the Facility location tutorial), or just giving a pointer to this thread would be totally fine by me

According to Discourse’s terms of service, posts are licensed under a CC BY-SA 4.0 license, so just to be safe maybe add a link to this discussion.