SHARE
TWEET

tsp

a guest Dec 11th, 2019 86 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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top