Struggling to define Jump constraints for modified TSP in a programmatic way

Hey All

I’m working on a traveling salesman like problem with a few additional constraints and I’m using the JUMP package.

One of the constraints is that certain nodes on the graph need to be performed sequentially

However, a some nodes can be done in-between the sequential ordered nodes if optimal

I have it working already with a slightly hard coded example. Below s_start and s_end are the nodes that need to be done sequentially, point is the node that can potentially be done in-between s_start and s_end.

@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
        )

However I need to recreate this for a given list of middle points sized N

So something like this for the specific example of 2 potential middle points

@constraint(
            model, 
            (x[s_start, s_end] + x[s_end, s_start])/1 + # start and end sequentially 
            (x[s_start, point_1] + x[point_1, s_end]) + (x[s_end, point_1] + x[point_1, s_start])/2  + # start, middle point_1, end sequentially
            (x[s_start, point_2] + x[point_2, s_end]) + (x[s_end, point_2] + x[point_2, s_start])/2 # start, middle point_2, end sequentially
            (x[s_start, point_1] + x[point_1, point_2] + x[point_2, s_end]) + (x[s_end, point_2] + x[point_2, point_1] + x[point_2, s_start])/2 # start, middle point_1, point_2, end sequentially
            == 1
        )

Do you have recommendations

The one idea I have is to build the constraint macro as a string and then use the Meta.parse(prog) command to run it. But I’m uncertain if scope rules will work for me.

Hi @Lee_Phillips, welcome to the forum.

Do you have the full formulation that you’re working with?

If i and j are sequential, this feels like you could use the MTZ formulation and constraint -1 <= u[i] - u[j] <= 1. If there is just a priority that i must be visited before j, then u[i] <= u[j]

1 Like

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]))

One option to add seqeuntialness in MTZ might be something like this:

model = Model()
@variable(model, 2 <= u[1:n] <= n, Int)
@variable(model, z, Bin) # if z = 0 then i => j else i => k => j
@constraint(model, u[i] - u[j] <= n * z)
@constraint(model, u[j] - u[i] <= n * z)
@constraint(model, u[i] - u[k] <= n * (1 - z))
@constraint(model, u[k] - u[i] <= n * (1 - z))
@constraint(model, u[k] - u[j] <= n * (1 - z))
@constraint(model, u[j] - u[k] <= n * (1 - z))

To simplify your constraint construction, you could otherwise do something like this:

term(x, i, j) = x[i, j] + x[j, i]
term(x, i, j, k) = (x[i, j] + x[j, k]) + (x[k, j] + x[j, i]) / 2
term(x, i, j, k, l) = (x[i, j] + x[j, k] + x[k, l]) + (x[l, k] + x[k, j] + x[j, i]) / 2
S = [
    (s_start, s_end),
    (s_start, point_1, s_end),
    (s_start, point_2, s_end),
    (s_start, point_1, point_2, s_end),
]
@constraint(model, sum(term(x, indices...) for indices in S) == 1)

(I haven’t tested, so there might be typos, etc, but you get the idea.)