I was first faced with a problem of grouping 1000 pairs into unique pairs to minimise the total distance - since this was a small number I solved it using linear optimisation:

Part 1 - Basic/Naive approach:

A bit more detail on it here but the gist is:

Create an objective function which wants to minimise all the entries in the lower triangle of a distance matrix (multiplied by the respective dummy). For example v1_10 means that point 1 and point 10 are connected. Such that each store is connected to only one (other store).

enter image description here

The result of the optimisation takes us from the first picture to the second. The structure of the code is thus:

shops = [[id, latitude, longitude],...]

The helper functions to create the optimisation:

def store_connect(shops):
    Cross stores and return connection e.g. "v1-2" by ID
    for i in range(len(shops)-1):
        for j in range(i+1,len(shops)):
            yield 'v'+str(shops[i][0])+'-'+str(shops[j][0])

def distances(shops):
    Return distance in miles for crosses
    for i in range(len(shops)-1):
        for j in range(i+1,len(shops)):
            yield haversine([shops[i][1],shops[i][2]],

def all_constraints(shops):
    Return constraint row
    for a in shops:
        cons = []
        for b in shops:
            if a[0]<b[0]:
            elif a[0]>b[0]:
        yield cons

And then the optimisation itself:

prob = LpProblem("Minimise distance",LpMinimize)

# Variables
x = []
xnames = []
for one in store_connect(shops):
print("Created %d variables" % len(x))

# Objective
prob += pulp.lpSum([i[1]*x[i[0]] for i in enumerate(distances(shops))])
print("Created objective")

# Constraints
counter = 0
for constraint_rw in all_constraints(shops):
    prob += pulp.lpSum(x[xnames.index(one_cn)] for one_cn in constraint_rw) == 1
    counter += 1
    if counter % 10 == 0:
        print("Added %d contraints out of %d" % (counter,len(shops)))

Which is straight-forward to solve directly or upload the .lp file to CPLEX:

#Using PULP - if not too complicated ... otherwise CPLEX:
# The status of the solution is printed to the screen
print("Status:", LpStatus[prob.status])

# Optimum variable values
sol = []
for v in prob.variables():
    if v.varValue:
print("Total minimised distance: ", value(prob.objective))  

shp_names = [row[0] for row in shops]
edges = []
for x in sol:
    v_list = (x)[1:].split("_") #remove the first letter 'v' and split
    conn = [it for sub in [shops[shp_names.index(float(v_list[0]))][1:3],
                           shops[shp_names.index(float(v_list[1]))][1:3]] for it in sub]

# Print diagram

Part 2 - Better approach (for lots of points):

If I expand the problem to 10,000 points I get roughly 50,000,000 variables. I was thinking of the following ways of simplifying the problem:

1) Use clustering to group the 10,000 points into something like 100 groups and then solve 100 times within those groups (however I'm not sure how to check whether the solution would be optimal).

2) Project the points; put them into a quad-tree and then instead of allowing connections between one point and any other point -> only allow connections to the closest 100 points (e.g.): however, I am not sure how to properly set this (the cut-off). I thought perhaps if a point is chosen that is the cut-off point then the algorithm will restart with a larger cut-off.

3) Run something like a minimum spanning tree or concorde tsp algorithm and then back-track the solution from that?

  • Are you looking for a Code Review or is there a specific something that is not working? The latter would suit the focussed Q&A that we do here but I think the former is too broad.
    – PolyGeo
    Commented Feb 5, 2016 at 10:24
  • @PolyGeo I don't think code review is appropriate because I have something that 'works' only for a small number of points. However I would like an approach that works for 10k points (or maybe even 100k points); which I think turns this into a different beast entirely.
    – mptevsion
    Commented Feb 5, 2016 at 10:26
  • It still seems to me that you are describing something that works for a small number of points but has performance that degrades severely as larger numbers of points are thrown at it. Code Review says that "If you are looking for feedback on a specific working piece of code from your project in the following areas… •Performance" then your question would be on-topic there.
    – PolyGeo
    Commented Feb 5, 2016 at 10:32
  • @PolyGeo I think the problem is that this works in n^2 time. So after 1000 points it doesn't just 'degrade' but stops working at all. So I was looking for a completely different approach which would scale to 10,000 points (it take a long time but it would work ... not like 100 years). However, if you still think that this should be Code Review I can repost - or can this be moved?
    – mptevsion
    Commented Feb 5, 2016 at 10:49

1 Answer 1


My suggestion is based on a method I apply to group smaller subcatchments into larger.

There are no streams, thus 1st step to compute Euclidean minimum spanning tree:

enter image description here

I guess it is some sort of optimisation already.

Pick your sink and create directed graph. Picture below shows "Flow Accumulation" in links, i.e. count of nodes discharging into it:

enter image description here

Pick links with flow greater or equal to desired "subcatchment" size, 10 in this case. Prune graph (remove nodes) at "source" branches. Picture show remaining nodes:

enter image description here

Nodes of each pruned subgraph form new "subcatchments". As one can see at 1st iteration 4 larger subcatchments were formed.

Repeat calculation of "Flow Accumulation" and pruning until no nodes left. Nodes at the picture below labelled by their parent subcatchment size. As one can see it ranges from 5 to 15, although target was set to 10.

Start making your pairs, going downstream. If count of nodes is even you are lucky, if not just add downstream subcatchment to the group.

enter image description here

Anyway it seem reasonable to use triangulation inside groups to compute distance matrix:

enter image description here

Hopefully it will reduce number of parameters to optimise


It took me 15 minutes to regroup original split of 500 points:

enter image description here

Into series each containing even count of points:

enter image description here

  • This seems to rely on luck, guesswork, and intuition. It also seems it could leave an awful lot of orphan (unpaired) points behind. Could you supply any reasons why it should work, if only approximately?
    – whuber
    Commented Feb 6, 2016 at 16:30
  • @whuber yes it will most likely produce what they call sub-optimal solution. Some manual intervention required to eliminate orphans using simple rule (reason you want?) that is total (or average) of 2 odd numbers result in even number. See updated answer
    – FelixIP
    Commented Feb 7, 2016 at 4:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.