Below is the full problem formulation
I like your recommendation, but I think the issue with only using MTZ is that I can not enforce sequentialness.
if u[i] <= u[j] then that doesn’t prevent u[i] <= u[x] <= u[y] < u[j] where nodes x and y are not meant to interrupt the sequence of i → j
To put another way, I don’t know if you can prevent the tour from doing a massive loop going from node i, to node j, passing through any number of additional nodes, when you want to enforce node i → node j with only a few optional candidate nodes.
# Tie point defects to Linear jobs
# point_i can be done between line_i_s and line_i_e
# hack it to be more effecient
using Pkg
println("----------------------------------------")
#if isfile("Project.toml") && isfile("Manifest.toml")
# Pkg.activate(".")
# Pkg.instantiate()
#end
using JuMP
using Random
using PlotlyJS
# using SCIP
# using Gurobi
using GLPK
using Pipe
# MAKE INPUT DATA
function basic_shortest_path_run()
n = 3 # number of nodes
lines = 3 # number of linear inspections
starts = 2 # number of start nodes (1 or 2)
extra_layers = [1, 2]
random_seed = 1
start_node = 1
figure_name = "third_one.png"
# make a random set of X-Y coordinates sized n
# return x vector, y vector, euclidian distance matrix x_i, y_i -> x_j, y_j for all pairwise combos
println("Generating Data")
rng = Random.MersenneTwister(random_seed)
X = 100 * rand(rng, n)
Y = 100 * rand(rng, n)
# d = [sqrt((X[i] - X[j]) ^ 2 + (Y[i] - Y[j])^2) for i in 1:n, j in 1:n]
# make random set of X-Y coordinate pairs sizes lines
line_Xs = 100 * rand(rng, Float64, (lines, 2))
line_Ys = 100 * rand(rng, Float64, (lines, 2))
row = line_Xs[1, :]
# make random start point
start_Xs = 100 * rand(rng, starts)
start_Ys = 100 * rand(rng, starts)
# plot input data
point_scatter = scatter(x = X, y = Y, mode = "markers", showlegend = false)
start_scatter = scatter(x = start_Xs, y = start_Ys, mode = "markers", showlegend = false)
lines_scatter = Vector{GenericTrace}(undef, lines)
for i in 1:lines
lines_scatter[i] = scatter(x = line_Xs[i, :], y = line_Ys[i, :], mode = "lines", showlegend = false)
end
savefig(plot([[point_scatter];lines_scatter; [start_scatter]]), "input_data.png")
# ORGANISE DATA
# need a list of nodes
# nodes are the points, the start/end points, the start/end of the lines
nodes_data = Dict()
# add in points
for i in 1:n
nodes_data["point_" * string(i)] = (X[i], Y[i])
end
# add in lines
for i in 1:lines
nodes_data["line_" * string(i) * "_start"] = (line_Xs[i, 1], line_Ys[i, 1])
nodes_data["line_" * string(i) * "_end"] = (line_Xs[i, 2], line_Ys[i, 2])
end
# add in start/end
if starts == 1
# duplicate single point as start/end
nodes_data["start"] = (start_Xs[1], start_Ys[1])
nodes_data["end"] = (start_Xs[1], start_Ys[1])
elseif starts == 2
# use first one as start and second on as end
nodes_data["start"] = (start_Xs[1], start_Ys[1])
nodes_data["end"] = (start_Xs[2], start_Ys[2])
else
# flag issue and stop
println("Bad start number")
exit()
end
nodes = string.(collect(keys(nodes_data)))
# track names of line segments
segment = ["line_" * string(i) for i in 1:lines]
# track names of segment start and ends
segment_starts = ["line_" * string(i) * "_start" for i in 1:lines]
segment_ends = ["line_" * string(i)* "_end" for i in 1:lines]
segment_points = [segment_starts;segment_ends]
# track names of point defects
defects = ["point_" * string(i) for i in 1:n]
# track names depots, start/end
depots = ["start", "end"]
# track jobs, jobs is segment + defects + depots
jobs = [segment;defects;depots]
# BUILD Distance matrix for all node pairs
# d = [sqrt((X[i] - X[j]) ^ 2 + (Y[i] - Y[j])^2) for i in 1:n, j in 1:n]
d = [sqrt((nodes_data[node_i][1] - nodes_data[node_j][1])^2 + (nodes_data[node_i][2] - nodes_data[node_j][2])^2) for node_i in nodes, node_j in nodes]
distance_matrix = @pipe d |> JuMP.Containers.DenseAxisArray(_, nodes, nodes)
println(nodes)
# build optimiser for TSP
# model = Model(Gurobi.Optimizer)
# model = Model(SCIP.Optimizer)
model = Model(GLPK.Optimizer)
println("Building Optimizer")
# simple adjacency matrix marking node joins, if x[i,j] = 1 means node_i <-> node_j
# TODO update this to handle points in lines requirement
# need to add dimension to matrix
@variable(model, x[nodes, nodes], Bin)
# time vector (distance) to track distance traveled to job
@variable(model, t[nodes] >= 0)
# constraint, inspection start / ends need to be joined
# TODO update this to handle points in lines requirement
# need to update constraints for each layer in the matrix
# each layer in the matrix will have their own requirement
for s in segment
local s_start = join([s, "_start"])
local s_end = join([s, "_end"])
local point = join(["point_", s[lastindex(s)]])
println(s)
println(point)
# Assume bi drectional
# slightly hard coded
# (x[s_start, s_end] + x[s_end, s_start])/1 + (x[s_start, point] + x[point, s_end]) + (x[s_end, point] + x[point, s_start])/2 = 0
@constraint(
model,
(x[s_start, s_end] + x[s_end, s_start])/1 + # start and end sequentially
(x[s_start, point] + x[point, s_end]) + (x[s_end, point] + x[point, s_start])/2 # start, middle point, end sequentially
== 1
)
#@constraint(
# model,
# (x[s_start, s_end])/1 +
# (x[s_start, point] + x[point, s_end])/2
# == 1
#)
#@constraint(
# model,
# (x[s_end, s_start])/1 +
# (x[s_end, point] + x[point, s_start])/2
# == 0
#)
# Assume bi directional
# @constraint(model, x[s_start, s_end] + x[s_end, s_start] == 1)
# Assume single directional
# @constraint(model, x[s_start, s_end] == 1)
# @constraint(model, x[s_end, s_start] == 0)
end
# constraint, each node can be left once
# each node can only have 1 exit, otherwise you would do the node more than once
# TODO update to handle points in lines
# adjust to blanket apply on all layers
for i in nodes
@constraint(model, sum(x[i, j] for j in nodes) <= 1)
end
println("test")
# constraint, each node can be entered onc,e
# each node can only have 1 entrance, else more than once
# TODO update to handle points in lines
# adjust to blanket apply on all layers
for j in nodes
@constraint(model, sum(x[i, j] for i in nodes) <= 1)
end
# constraint, you want to start at start and end at start and end
# start to any next node needs to be <= 1
# any prev node to end needs to be <= 1
# TODO handle points in lines
# adjust to blanket apply on all layers
@constraint(model, sum(x["start", j] for j in nodes) <= 1)
@constraint(model, sum(x[i, "end"] for i in nodes) <= 1)
# constraint, no self loops
# TODO apply on all point in line layers
for i in nodes
@constraint(model, x[i,i] == 0)
end
# constraint, subtour contraints
# if a node is connected to the end of a chain, then all the previous nodes in the chain need to have a smaller total distance traveled
# TODO need to handle point in lines layers
# hack to make it very beneficial to be on layer 2
for i in nodes, j in nodes
@constraint(model, t[j] >= t[i] + distance_matrix[i,j] - (1000 * (1 - x[i,j]))) # big M is related to max distance between any two points, corners
end
# constraint, each pair of nodes can only be joined maximum one time
# TODO blanket apply on all layers
for i in nodes, j in nodes
@constraint(model, x[i,j] + x[j,i] <= 1)
end
# constraint, each node can only have 1 exit
# except for end which has 0 exits
# TODO blanket apply on all layers
for i in nodes
if i != "end"
@constraint(model, 1 == sum(x[i,j] for j in nodes))
end
end
# constraint, each node can only have 1 entrance
# except for start which has 0 entrants
# TODO blanket apply on all layers
for i in nodes
if i != "start"
@constraint(model, 1 == sum(x[j,i] for j in nodes))
end
end
# objective
# minimise total of product of adjacency matrix and distance matrix
# TODO apply on all layers
@expression(model, tot_dist, sum(coalesce(distance_matrix[i, j], 10000) * x[i,j] for i in nodes, j in nodes))
@objective(model, Min, tot_dist)
# Model summary
show(model)
# reset attributes and run model
#MOIU.reset_optimizer(model) # Resets attributes?
#set_optimizer_attribute(model, "SolutionLimit", 10000)
#set_optimizer_attribute(model, "TimeLimit", 1020)
#set_optimizer_attribute(model, "NoRelHeurTime", 1020)
# set_optimizer_attribute(model, "limits_time", 10000)
println("Running Optimizer")
optimize!(model)
function selected_edges(x::Matrix{Float64}, n)
return Tuple{Int,Int}[(i, j) for i in 1:n, j in 1:n if x[i, j] > 0.5]
end
println("Ploting output")
# TODO update to handle multiple layers
println(value.(model[:x]))
function plot_tour_v2(nodes_data, sol_mat, nodes)
# plot input data
plot_list = Vector{GenericTrace}(undef, 0)
for node in nodes
data = nodes_data[node]
point_scatter = scatter(x = [data[1]], y = [data[2]], text = [node], mode = "markers+text", showlegend = false)
plot_list = [plot_list;[point_scatter]]
end
# go through tour from start node
node_counter = 0
node_counter = node_counter + 1
row = sol_mat["start", :]
current_index = "start"
println(current_index)
next_index = 0
for node in nodes
if row[node] == 1
next_index = node
break
end
end
# next_index = findfirst(item->item==1, row)
println(next_index)
line_plot = scatter(
x = [nodes_data[current_index][1], nodes_data[next_index][1]],
y = [nodes_data[current_index][2], nodes_data[next_index][2]],
mode = "lines",
marker=attr(
symbol = "arrow"
),
showlegend = false
)
text_plot = scatter(
x = [nodes_data[next_index][1]],
y = [nodes_data[next_index][2]],
text = [node_counter],
mode = "text", showlegend = false
)
plot_list = [plot_list;[line_plot];[text_plot]]
while next_index != "end"
node_counter = node_counter + 1
row = sol_mat[next_index, :]
current_index = next_index
for node in nodes
if row[node] == 1
next_index = node
break
end
end
#next_index = findfirst(item->item==1, row)
println(next_index)
line_plot = scatter(
x = [nodes_data[current_index][1], nodes_data[next_index][1]],
y = [nodes_data[current_index][2], nodes_data[next_index][2]],
mode = "lines",
marker=attr(
symbol = "arrow"
),
showlegend = false
)
text_plot = scatter(
x = [nodes_data[next_index][1]],
y = [nodes_data[next_index][2]],
text = [node_counter],
mode = "text", textposition = "bottom", showlegend = false
)
plot_list = [plot_list;[line_plot];[text_plot]]
end
return plot(plot_list)
end
function plot_tour(X, Y, x, start_index, n)
# plot points
point_scatter = scatter(x = X, y = Y, mode = "markers", showlegend = false)
# get first row
node_counter = 1
row = x[start_index, :]
current_index = start_index
next_index = findall(x->x==1, row)[1]
lines_and_texts = Vector{GenericTrace}(undef, n * 2)
# plot line and text
lines_and_texts[node_counter] = scatter(
x = [X[current_index], X[next_index]],
y = [Y[current_index], Y[next_index]],
mode = "markers+lines",
marker=attr(
symbol = "arrow"
),
showlegend = false
)
lines_and_texts[node_counter + n] = scatter(
x = [X[current_index], X[next_index]],
y = [Y[current_index], Y[next_index]],
text = [node_counter],
mode = "text", showlegend = false
)
node_counter = node_counter + 1
while next_index != start_index
row = x[next_index, :]
current_index = next_index
next_index = findall(x->x==1, row)[1]
# plot line and text
lines_and_texts[node_counter] = scatter(
x = [X[current_index], X[next_index]],
y = [Y[current_index], Y[next_index]],
mode = "markers+lines",
marker=attr(
symbol = "arrow"
),
showlegend = false
)
lines_and_texts[node_counter + n] = scatter(
x = [X[current_index], X[next_index]],
y = [Y[current_index], Y[next_index]],
text = [node_counter],
mode = "text", showlegend = false
)
node_counter = node_counter + 1
end
return plot([[point_scatter];lines_and_texts])
end
path_plot = plot_tour_v2(nodes_data, value.(model[:x]), nodes)
savefig(path_plot, figure_name)
println("Done!")
end
# RUN CODE
basic_shortest_path_run()
# plot_tour(X, Y, value.(iterative_model[:x]))