mtravelling-salesmen-variant.jl

#=
main:
- Julia version:
- Author: quangio
- Date: 2019-05-06
=#



using JuMP
using SCIP
# include("input.jl")

model = Model(with_optimizer(SCIP.Optimizer))


dist_mat_gr17 = [
  0 633 257  91 412 150  80 134 259 505 353 324  70 211 268 246 121
633   0 390 661 227 488 572 530 555 289 282 638 567 466 420 745 518
257 390   0 228 169 112 196 154 372 262 110 437 191  74  53 472 142
 91 661 228   0 383 120  77 105 175 476 324 240  27 182 239 237  84
412 227 169 383   0 267 351 309 338 196  61 421 346 243 199 528 297
150 488 112 120 267   0  63  34 264 360 208 329  83 105 123 364  35
 80 572 196  77 351  63   0  29 232 444 292 297  47 150 207 332  29
134 530 154 105 309  34  29   0 249 402 250 314  68 108 165 349  36
259 555 372 175 338 264 232 249   0 495 352  95 189 326 383 202 236
505 289 262 476 196 360 444 402 495   0 154 578 439 336 240 685 390
353 282 110 324  61 208 292 250 352 154   0 435 287 184 140 542 238
324 638 437 240 421 329 297 314  95 578 435   0 254 391 448 157 301
 70 567 191  27 346  83  47  68 189 439 287 254   0 145 202 289  55
211 466  74 182 243 105 150 108 326 336 184 391 145   0  57 426  96
268 420  53 239 199 123 207 165 383 240 140 448 202  57   0 483 153
246 745 472 237 528 364 332 349 202 685 542 157 289 426 483   0 336
121 518 142  84 297  35  29  36 236 390 238 301  55  96 153 336   0] # shoud return 2085
dist_mat_5 = [
0.0  3.0  4.0  2.0  7.0
3.0  0.0  4.0  6.0  3.0
4.0  4.0  0.0  5.0  8.0
2.0  6.0  5.0  0.0  6.0
7.0  3.0  8.0  6.0  0.0
] # shoud return 19.0

n = 17 # number of vertices
c = zeros(n)   # cost assigned to each vertex (weight needed to be picked up)
d = dist_mat_gr17 # TestSet.dist_mat_48 # distance matrix of the graph

w0 = 1
no_salesman = 1
cap = 100

range = collect(1:n)

for i = 1:n
    d[i, i] = 0
end


@variable(model, y[1:n, 1:n])
@variable(model, x[1:n, 1:n], Bin)

@objective(model, Min, sum(y[i, j] * d[i, j] for i = 1:n, j = range[1:n .!= i]))

@constraint(model, sum(x[1, i] for i = 2:n) == no_salesman)
@constraint(model, sum(x[i, 1] for i = 2:n) == no_salesman)

for i = 2:n
    @constraint(model, y[1, i] == w0 * x[1, i])
    @constraint(model, sum(x[i, j] for j = range[1:n .!= i]) == 1)
    @constraint(model, sum(x[j, i] for j = range[1:n .!= i]) == 1)
    @constraint(model, c[i] == sum(y[i, j] - y[j, i] for j = range[1:n .!= i]))
end

for i = 1:n
    for j = 1:n
        if i == j continue end
        @constraint(model, y[i, j] >= (w0 + c[i]) * x[i, j])
        @constraint(model, y[i, j] <= (w0 + cap - c[j]) * x[i, j])
    end
end


function solved(m)
    x_val = JuMP.value.(x)
    cycle_idx = [1] # we can adapt this to support mTSP, just keep this for simplicity for now
    while true
        v, idx = findmax(x_val[cycle_idx[length(cycle_idx)], :])
        if idx == 1
            break
        else
            push!(cycle_idx, idx)
        end
    end
    if length(cycle_idx) < n
        @constraint(m, sum(x[i, j] for i = cycle_idx, j = cycle_idx[cycle_idx .!= i]) <= length(cycle_idx) - 1)
        return false
    end
    return true
end

if no_salesman == 0
    JuMP.optimize!(model)
    while !solved(model)
        JuMP.optimize!(model)
    end
end

if no_salesman >= 1
    @variable(model, 0 <= u[1:n] <= n - 1)
    for i = 2:n
        for j = collect(2:n)[2:n .!= i]
            @constraint(model, u[i] - u[j] + n * x[i, j] <= n - 1)
        end
    end
    @time JuMP.optimize!(model)
end

println(string("Objective value:", JuMP.objective_value(model)))

println("Arcs: ")

for i = range
    for j = range
        if JuMP.value(x[i, j]) > 0.0
            println(string(i, "->", j, "[label=", d[i, j], "]")) # you can use graphviz to visualize this
        end
    end
end

open("model2.lp", "w") do f
    print(f, model)
end

travelling-salesmen-dp.cpp

#include <iostream>
#include <vector>
#include <cstdio>
#include <algorithm>

using namespace std;

unsigned n, src = 0;
vector<vector<unsigned> > graph, dp;
vector<int> c;

const unsigned INIT_SIZE = 1;

inline void init() {
    for (size_t i = 0; i < n; ++i)
        dp[1u << i][i] = graph[src][i] * INIT_SIZE;
}

int TSP(unsigned status, unsigned x) {
    if (dp[status][x] != -1)
        return dp[status][x];

    int mask = 1u << x;
    dp[status][x] = int(1e9);
    for (size_t i = 0; i < n; ++i) {
        if (i != x && (status & (1u << i))) {
            unsigned prev = status - mask;
            int pickup = 0;
            for (size_t j = 0; j < n; ++j)
                if (prev & (1u << j))
                    pickup += c[j];
            dp[status][x] = min(dp[status][x], TSP(prev, i) + (INIT_SIZE + pickup) * graph[i][x]);
        }
    }
    return dp[status][x];
}

int main() {
    cin >> n;
    graph = vector<vector<unsigned>>(n, vector<unsigned>(n));
    dp = vector<vector<unsigned>>(1u << n, vector<unsigned>(n, -1));
    c = vector<int>(n);
    for (int i = 0; i < n; ++i) cin >> c[i];
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            cin >> graph[i][j];
    init();
    int r = TSP((1u << n) - 1, src);
    cout << r << endl;
    return 0;
}

For the first file, please install Julia and SCIP in your system. After installing Julia, open your terminal, type julia, then ], add JuMP, and add SCIP. Type julia mtravelling-salesmen-variant.jl in your command line to execute the code. You can change number of cities n, the weights c, distance matrix d, no_salesman... (from line 43 to 49) For instance,

n = 4
c = [1 1 1 1]
d = [0 1 9 9; 1 0 9 9; 9 9 0 1; 9 9 1 0]
# 
w0 = 1
no_salesman = 1
cap = 100

The C++ should be compiled using g++ -Ofast travelling-salesmen-dp.cpp -o tsp.out. After executing ./tsp.out, you can input the number of cities (integer n), weights of vertices (n integers), and distance matrix (n x n integers). Sample input

3
1 1 1
0 1 1
1 0 1
1 1 0

Simulated Annealing is implemented here (based on https://github.com/evanfields/TravelingSalesmanHeuristics.jl)
https://github.com/BoyanXu/simulated_annealing