Thank you again for suggestions Oscar and Victor for the example. After a lot of trying different things I’ve come at it from the opposite direction of solving the spatial problem of multiple simultaneous TSP’s first then worrying about other objectives.
I used your suggested school bus routing code as a base to calculate a series of networks with unique starting points.
Essentially instead of several buses from one point, a single bus or route runs from each of the designated ‘source’ points.
Since each demand is associated with only one node, the y and z variables and associated constraints were consolidated.
Two objectives maximise demand and minimise distance.
This runs successfully on a toy example (10 nodes, 3 source + 7 demand) with either HiGHS or GLPK.
However, on a slightly more meaningful example (70 nodes, 20 source + 50 demand) the run time quickly became intractable.
I tried adding lazy constraints to improve the efficiency but haven’t got the syntax quite right as we want to read in the 3-dimensional x matrix but use a 2-dimensional matrix for the subtour elimination.
From what I can see, lazy constraints aren’t supported by HiGHS so have tried switching to GLPK. I assume I will also need to turn off the MTZ constraints.
Any suggestions or corrections welcome.
Below is revised code
using HiGHS
using GLPK
import MultiObjectiveAlgorithms as MOA
distance = [
0 3.16 1.41 1 2.24 2.44 2.64 2.84 3.04 3.24
3.16 0 4 3.61 4.12 4.32 4.52 4.72 4.92 5.12
1.41 4 0 2.24 1 1.2 1.4 1.6 1.8 2
1 3.61 2.24 0 3.16 3.36 3.56 3.76 3.96 4.16
2.24 4.12 1 3.16 0 0.2 0.4 0.6 0.8 1
2.44 4.32 1.2 3.36 0.2 0 0.2 0.4 0.6 0.8
2.64 4.52 1.4 3.56 0.4 0.2 0 0.2 0.4 0.6
2.84 4.72 1.6 3.76 0.6 0.4 0.2 0 0.2 0.4
3.04 4.92 1.8 3.96 0.8 0.6 0.4 0.2 0 0.2
3.24 5.12 2 4.16 1 0.8 0.6 0.4 0.2 0
]
demands = [
0 0 0 2 2 1 4 5 3 4
] #First 3 nodes are source nodes, last 7 nodes are demand nodes
a = [
0 1 1
1 0 1
1 1 0
0 0 0
0 0 0
1 0 0
0 0 0
0 0 0
0 0 0
0 0 0
] #Allowed node combinations matrix (sources x demands): Use 1 to disallow nodes
V = 1:size(distance, 1)
supply_capacity = [4,5,7] #Allows unique capacity for each bus/source
n = length(supply_capacity) #Modify to set number of routes based on supply capacity vector
#Lazy constraints functions
#https://jump.dev/JuMP.jl/stable/tutorials/algorithms/tsp_lazy_constraints/
function subtour(edges::Vector{Tuple{Int,Int}}, n)
shortest_subtour, unvisited = collect(1:n), Set(collect(1:n))
while !isempty(unvisited)
this_cycle, neighbors = Int[], unvisited
while !isempty(neighbors)
current = pop!(neighbors)
push!(this_cycle, current)
if length(this_cycle) > 1
pop!(unvisited, current)
end
neighbors =
[j for (i, j) in edges if i == current && j in unvisited]
end
if length(this_cycle) < length(shortest_subtour)
shortest_subtour = this_cycle
end
end
return shortest_subtour
end
#Let us declare a helper function `selected_edges()` that will be repeatedly
#used in what follows.
function selected_edges(x_k::Matrix{Float64}, n)
return Tuple{Int,Int}[(i, j) for i in 1:n, j in 1:n if x_k[i, j] > 0.5]
end
#Other helper functions for computing subtours:
subtour(x::Matrix{Float64}) = subtour(selected_edges(x_k, size(x_k, 1)), size(x_k, 1))
subtour(x::AbstractMatrix{VariableRef}) = subtour(value.(x_k))
#model = Model(() -> MOA.Optimizer(HiGHS.Optimizer)) #For multiple objectives
model = Model(() -> MOA.Optimizer(GLPK.Optimizer)) #For multiple objectives
set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ALL)
@variables(model, begin
x[i in V, j in V, k in 1:n], Bin #A set of symmetric matrices, one for each potential route/bus. k is the number of school buses/sources
y[i in V, k in 1:n], Bin #1 if bus k visits stop i, 0 otherwise
2 <= u[i in V, k in 1:n] <= length(V) #A route vector of 5 nodes for each school bus
end)
fix.(u[1, :], 1; force = true) #Set first node of route vector u to always be 1
@expression(model, p4_expr, sum(demands[i] * y[i, k] for i in V, k in 1:n))
@expression(model, dist_expr,-sum(distance[i, j] .* x[i, j, k] for i in V, j in V, k in 1:n))
@objective(
model,
Max, [p4_expr,-dist_expr]
) #Minimise distance and maximise demands supplied
@constraints(model, begin
#Spatial network constraints
# For all nodes and routes, flow out == flow in
[i in V, k in 1:n], sum(x[i, :, k]) == sum(x[:, i, k])
# Flow in and out of all nodes other than current source node is 1
[i in V, k in 1:n; i != k], sum(x[i, :, :]) <= 1
[i in V, k in 1:n; i != k], sum(x[:, i, :]) <= 1
# MTZ constraints revised for multiple sources
[i in V, j in V, k in 1:n; i != k && j != k],
u[i, k] - u[j, k] + 1 <= (n - 1) * (1 - x[i, j, k])
# y[i, k] is 1 if node i is in route k
[i in V, k in 1:n], sum(x[i, :, k]) == y[i, k]
# Each node can be in at most one route
[i in V, k in 1:n], sum(y[i, :]) <= 1
#Demand and supply constraints
[k in 1:n], sum(demands[:] .* y[:, k]) <= supply_capacity[k] #Check demands on each route do not exceed supply capacity of source
[i in V], sum(y[i, :]) <= 1 #Each node on only one route
#Constrain which source-demand combinations are allowed or not allowed
[i in V, k in 1:n], y[i, k] * a[i, k] == 0
end)
#Lazy constraints
function subtour_elimination_callback(cb_data) #Not working
status = callback_node_status(cb_data, model)
if status != MOI.CALLBACK_NODE_STATUS_INTEGER
return # Only run at integer solutions
end
# x_all = callback_value.(cb_data, model[:x])
x_all = callback_value.(Ref(cb_data), model[:x])
# m = callback_value.(cb_data, model[:n])
m=3 #hard-coded for testing, bring in through callback
for k = 1:m
x_k = x_all[:,:,k] #Create a matrix based on x for each source node or k
cycle = subtour(x_k)
if !(1 < length(cycle) < n)
return # Only add a constraint if there is a cycle
end
println("Found cycle of length $(length(cycle))")
S = [(i, j) for (i, j) in Iterators.product(cycle, cycle) if i < j]
con = @build_constraint(
sum(x_k[i, j, k] for (i, j) in S) <= length(cycle) - 1, #Modified to use x_k
)
MOI.submit(model, MOI.LazyConstraint(cb_data), con)
end #k loop
return
end
set_attribute(model, MOI.LazyConstraintCallback(), subtour_elimination_callback)
optimize!(model)
solution_summary(model)
#Print route results
results = result_count(model)
if results >= 1
for r in 1:results
println("Result $r routes")
for k in 1:n
print("Route $k: $k")
i = k
while true
for j in V
if value(x[i, j, k];result=r) > 0.5
print(" => $j")
i = j
break
end
end
if i == k
println()
break
end
end #while
end #for k
end #for r
end #if
println()