tsp

NS = 80;
xmax = 20; ymax = 20; delta = 1.0;
xs = zeros(NS, 1);
ys = zeros(NS, 1);
xs(1) = xmax*rand;
ys(1) = ymax*rand;
n = 1;
while n < NS
    x = xmax * rand;
    y = ymax * rand;
    if all(hypot(x-xs,y-ys)> delta)
        n = n+1;
        xs(n) =x; ys(n) = y;
    end
end



dist = zeros(NS, NS);
for i= 1:NS
    for j = 1:NS
        dist(i, j) = hypot(xs(i)-xs(j), ys(i)-ys(j));
    end
end

f = dist(:);
intcon = 1:length(f);
lb = zeros(length(f), 1);
ub = ones(length(f), 1);
i = 1:NS;
ub(i + (i-1)*NS) = 0.

A = zeros(NS*NS, length(f));
b = zeros(size(A, 1), 1);
n = 1;
for i = 1:NS
    for j = 1:NS
        A(n, i+(j-1)*NS) = 1;
        A(n, j+(i-1)*NS) = 1;
        b(n) = 1;
        n = n+1;
    end
end
Aeq = zeros(2*NS, length(f));
beq = zeros(2*NS, 1);
n = 1;
for i = 1:NS
    for j = 1:NS
        Aeq(n, i+(j-1)*NS) = 1;
        Aeq(n+1, j+(i-1)*NS) = 1;
    end
    beq(n) = 1;
    beq(n+1) = 1;
    n = n + 2;
end

while 1
    [xy, fval, exitflag] = intlinprog(f, intcon, A, b, Aeq, beq, lb, ub);
    xy = round(xy);

    % draw
    plot(xs, ys, 'r*')
    hold on
    for j = 1:NS
        for i = 1:NS
            if xy(i+(j-1)*NS) == 1
                plot([xs(i) xs(j)], [ys(i), ys(j)], 'b-')
            end
        end
    end
    hold off
    pause
    % draw end

    xy = reshape(xy, NS, NS);

    stops = zeros(NS, 1);
    for n = 1:NS
        subtour = find(stops==0, 1);
        while 1
           i = subtour(end);
           stops(i) = 1;
           j = find(xy(i, :), 1);
           if j == subtour(1), break; end
           subtour = [subtour j];
        end
        if n == 1 && all(stops==1)
            break;
        end

        A = [A; zeros(1, length(f))];
        for i = subtour
            for j = subtour
                A(end, i+(j-1)*NS) = 1;
            end
        end
        b = [b; length(subtour)-1];
        if all(stops==1) break; end
    end

    if n == 1 break; end;
end

plot(xs, ys, 'r*')
hold on
for j = 1:NS
    for i = 1:NS
        if xy(i+(j-1)*NS) == 1
            plot([xs(i) xs(j)], [ys(i), ys(j)], 'b-')
        end
    end
end
hold off
pause