Guest User

tsp

a guest
Dec 11th, 2019
96
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. NS = 80;
  2. xmax = 20; ymax = 20; delta = 1.0;
  3. xs = zeros(NS, 1);
  4. ys = zeros(NS, 1);
  5. xs(1) = xmax*rand;
  6. ys(1) = ymax*rand;
  7. n = 1;
  8. while n < NS
  9.     x = xmax * rand;
  10.     y = ymax * rand;
  11.     if all(hypot(x-xs,y-ys)> delta)
  12.         n = n+1;
  13.         xs(n) =x; ys(n) = y;
  14.     end
  15. end
  16.  
  17.  
  18.  
  19. dist = zeros(NS, NS);
  20. for i= 1:NS
  21.     for j = 1:NS
  22.         dist(i, j) = hypot(xs(i)-xs(j), ys(i)-ys(j));
  23.     end
  24. end
  25.  
  26. f = dist(:);
  27. intcon = 1:length(f);
  28. lb = zeros(length(f), 1);
  29. ub = ones(length(f), 1);
  30. i = 1:NS;
  31. ub(i + (i-1)*NS) = 0.
  32.  
  33. A = zeros(NS*NS, length(f));
  34. b = zeros(size(A, 1), 1);
  35. n = 1;
  36. for i = 1:NS
  37.     for j = 1:NS
  38.         A(n, i+(j-1)*NS) = 1;
  39.         A(n, j+(i-1)*NS) = 1;
  40.         b(n) = 1;
  41.         n = n+1;
  42.     end
  43. end
  44. Aeq = zeros(2*NS, length(f));
  45. beq = zeros(2*NS, 1);
  46. n = 1;
  47. for i = 1:NS
  48.     for j = 1:NS
  49.         Aeq(n, i+(j-1)*NS) = 1;
  50.         Aeq(n+1, j+(i-1)*NS) = 1;
  51.     end
  52.     beq(n) = 1;
  53.     beq(n+1) = 1;
  54.     n = n + 2;
  55. end
  56.  
  57. while 1
  58.     [xy, fval, exitflag] = intlinprog(f, intcon, A, b, Aeq, beq, lb, ub);
  59.     xy = round(xy);
  60.    
  61.     % draw
  62.     plot(xs, ys, 'r*')
  63.     hold on
  64.     for j = 1:NS
  65.         for i = 1:NS
  66.             if xy(i+(j-1)*NS) == 1
  67.                 plot([xs(i) xs(j)], [ys(i), ys(j)], 'b-')
  68.             end
  69.         end
  70.     end
  71.     hold off
  72.     pause
  73.     % draw end
  74.    
  75.     xy = reshape(xy, NS, NS);
  76.  
  77.     stops = zeros(NS, 1);
  78.     for n = 1:NS
  79.         subtour = find(stops==0, 1);
  80.         while 1
  81.            i = subtour(end);
  82.            stops(i) = 1;
  83.            j = find(xy(i, :), 1);
  84.            if j == subtour(1), break; end
  85.            subtour = [subtour j];
  86.         end
  87.         if n == 1 && all(stops==1)
  88.             break;
  89.         end
  90.        
  91.         A = [A; zeros(1, length(f))];
  92.         for i = subtour
  93.             for j = subtour
  94.                 A(end, i+(j-1)*NS) = 1;
  95.             end
  96.         end
  97.         b = [b; length(subtour)-1];
  98.         if all(stops==1) break; end
  99.     end
  100.    
  101.     if n == 1 break; end;
  102. end
  103.  
  104. plot(xs, ys, 'r*')
  105. hold on
  106. for j = 1:NS
  107.     for i = 1:NS
  108.         if xy(i+(j-1)*NS) == 1
  109.             plot([xs(i) xs(j)], [ys(i), ys(j)], 'b-')
  110.         end
  111.     end
  112. end
  113. hold off
  114. pause
RAW Paste Data