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

Specific Domains Optimization (Mathematical)
1

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.

2

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
3

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

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

5

Hey all

Apologies for leavening this post bare for so long

I managed to solve the issue for myself, making use of the AffExpr functionality in JUMP (Expressions · JuMP)

The basic idea is a progressively add to an AffExpr as I am looping through a for loop

the end result is of the AffExpr is constrained to to be 1, allowing for the xor nature of the addition operation to operate properly

A full example solution is pasted below, its a bit dense but hopefully easy to understand

# 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
using Combinatorics


# 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 = "fourth_one_basic_2.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
    @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
    segment_expressions = Dict()
    new_segment_expressions = Dict()
    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)

        # soft coded bi directional
        # small issues with middle points needing to be joined breaking the constraints
        # middle_point_list = ["point_1", "point_2"]
        middle_point_list = [point]
        println(middle_point_list)

        thingy = parse(Int, s[lastindex(s)])
        #@expression(model, holder_var, 
        #    (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
        #)
        #set_name(holder_var, join(["line_variable_", s[lastindex(s)]]))
        
        # build base ex for start and end
        combined_forward_ex = AffExpr(0.0)
        combined_backwards_ex = AffExpr(0.0)

        # by directional
        # add_to_expression!(base_ex, (x[s_start, s_end] + x[s_end, s_start])/1)
        for sub_list in powerset(middle_point_list)
            forward_ex = AffExpr(0.0)
            backwards_ex = AffExpr(0.0)

            # add_to_expression!(next_ex)
            if length(sub_list) == 0
                add_to_expression!(forward_ex, x[s_start, s_end])
                add_to_expression!(backwards_ex, x[s_end, s_start])
            else

                # build forward pass
                println(sub_list)
                # link start and first point of list
                cur_index = 1
                add_to_expression!(forward_ex, x[s_start, sub_list[cur_index]])
                while cur_index < length(sub_list)
                    add_to_expression!(forward_ex, x[sub_list[cur_index], sub_list[cur_index + 1]])
                    cur_index += 1
                end
                add_to_expression!(forward_ex, x[sub_list[cur_index], s_end])

                # build backwards pass
                cur_index = length(sub_list)
                add_to_expression!(backwards_ex, x[s_end, sub_list[cur_index]])
                while cur_index > 1
                    add_to_expression!(backwards_ex, x[sub_list[cur_index], sub_list[cur_index - 1]])
                    cur_index -= 1
                end
                add_to_expression!(backwards_ex, x[sub_list[cur_index], s_start])
            end
            
            # combined_ex = AffExpr(0.0)
            add_to_expression!(combined_forward_ex, forward_ex/(length(sub_list) + 1))
            add_to_expression!(combined_backwards_ex, backwards_ex/((length(sub_list) + 1)))
            # add_to_expression!(base_ex, combined_ex/(length(sub_list) + 1))
        end

        if false
            # bi directional
            base_ex = AffExpr(0.0)
            add_to_expression!(base_ex, combined_forward_ex)
            add_to_expression!(base_ex, combined_backwards_ex)
            new_segment_expressions[thingy] = base_ex
            @constraint(
                model,
                base_ex == 1
            )
        else
            # uni directional
            @constraint(
                model,
                combined_forward_ex == 1
            )
            @constraint(
                model,
                combined_backwards_ex == 0
            )
            new_segment_expressions[thingy] = [combined_forward_ex, combined_backwards_ex]
        end
    end

    # constraint, each node can be left once
    # each node can only have 1 exit, otherwise you would do the node more than once
    for i in nodes
        @constraint(model, sum(x[i, j] for j in nodes) <= 1)
    end

    println("test")
    # constraint, each node can be entered once
    # each node can only have 1 entrance, else more than once
    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
    @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
    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
    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
    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
    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
    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
    @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")
    #println(value.(model[:x]))
    #println(segment_expressions[1])
    #println(typeof(segment_expressions[1]))
    println(new_segment_expressions[1])
    println(typeof(new_segment_expressions[1]))


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