禁忌搜索算法的TSP问题求解

Halcom

禁忌搜索算法的TSP问题求解：

1. clc,clear,close all
   
2. warning off
   
3. x=[82 91 12 92 63 9 28 55 96 97 15 98 96 49 80 12 80 55];
   
4. y=[14 42 92 80 96 66 3 85 94 68 76 75 39 66 17 78 80 45];
   
5. [Position, BestCost]= Tabu_Search_TSP(x,y,50);
   
6. 
   
7. figure(1);
   
8. hold on
   
9. plot(x,y,'o')
   
10. for i=1:length(x)
    
11.     if i<length(x)
    
12.         plot( [x(Position(i)),x(Position(i+1))], [y(Position(i)),y(Position(i+1))],'-','linewidth',2 )
    
13.     else
    
14.         plot( [x(Position(i)),x(Position(1))], [y(Position(i)),y(Position(1))],'-','linewidth',2 )
    
15.     end
    
16. end
    
17. hold off

复制代码

函数如下：

1. function [Position, BestCost]= Tabu_Search_TSP(x,y,MaxIt)
   
2. % x=[82 91 12 92 63 9 28 55 96 97 15 98 96 49 80 12 80 55];
   
3. % y=[14 42 92 80 96 66 3 85 94 68 76 75 39 66 17 78 80 45];
   
4. 
   
5. model = CreateModel_ysw(x, y);                  % Create TSP Model
   
6. CostFunction=@(tour) TourLength(tour, model);   % Cost Function
   
7. ActionList=CreatePermActionList(model.n);       % Action List
   
8. nAction=numel(ActionList);                      % Number of Actions
   
9. %% Tabu Search Parameters
   
10. % MaxIt=50;                   % Maximum Number of Iterations
    
11. TL=round(0.5*nAction);      % Tabu Length
    
12. %% Initialization
    
13. % Create Empty Individual Structure
    
14. empty_individual.Position=[];
    
15. empty_individual.Cost=[];
    
16. % Create Initial Solution
    
17. sol=empty_individual;
    
18. sol.Position=randperm(model.n);
    
19. sol.Cost=CostFunction(sol.Position);
    
20. % Initialize Best Solution Ever Found
    
21. BestSol=sol;
    
22. % Array to Hold Best Costs
    
23. BestCost=zeros(MaxIt,1);
    
24. % Initialize Action Tabu Counters
    
25. TC=zeros(nAction,1);
    
26. %% Tabu Search Main Loop
    
27. for it=1:MaxIt
    
28.     bestnewsol.Cost=inf;
    
29.     % Apply Actions
    
30.     for i=1:nAction
    
31.         if TC(i)==0
    
32.             newsol.Position=DoAction(sol.Position,ActionList{i});
    
33.             newsol.Cost=CostFunction(newsol.Position);
    
34.             newsol.ActionIndex=i;
    
35. 
    
36.             if newsol.Cost<=bestnewsol.Cost
    
37.                 bestnewsol=newsol;
    
38.             end
    
39.         end
    
40.     end
    
41.     % Update Current Solution
    
42.     sol=bestnewsol;
    
43.     % Update Tabu List
    
44.     for i=1:nAction
    
45.         if i==bestnewsol.ActionIndex
    
46.             TC(i)=TL;               % Add To Tabu List
    
47.         else
    
48.             TC(i)=max(TC(i)-1,0);   % Reduce Tabu Counter
    
49.         end
    
50.     end
    
51.     % Update Best Solution Ever Found
    
52.     if sol.Cost<=BestSol.Cost
    
53.         BestSol=sol;
    
54.     end
    
55.     % Save Best Cost Ever Found
    
56.     BestCost(it)=BestSol.Cost;
    
57.     % Show Iteration Information
    
58. %     disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
59.    
    
60. %     % Plot Best Solution
    
61. %     figure(1);
    
62. %     PlotSolution(BestSol, model);
    
63. %     pause(0.01);
    
64.     % If Global Minimum is Reached
    
65.     if BestCost(it)==0
    
66.         break;
    
67.     end
    
68. end
    
69. % BestCost=BestCost(1:it);
    
70. Position = BestSol.Position;
    
71. BestCost = BestSol.Cost;
    
72. 
    
73. end
    
74. 
    
75. function model=CreateModel_ysw(x,y)
    
76. %     x=[82 91 12 92 63 9 28 55 96 97 15 98 96 49 80 12 80 55];
    
77. %     y=[14 42 92 80 96 66 3 85 94 68 76 75 39 66 17 78 80 45];
    
78.     n=numel(x);
    
79.     d=zeros(n,n);
    
80.     for i=1:n-1
    
81.         for j=i+1:n
    
82.             d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
    
83.             d(j,i)=d(i,j);
    
84.         end
    
85.     end
    
86.     xmin=0;
    
87.     xmax=100;
    
88.     ymin=0;
    
89.     ymax=100;
    
90.     model.n=n;
    
91.     model.x=x;
    
92.     model.y=y;
    
93.     model.d=d;
    
94.     model.xmin=xmin;
    
95.     model.xmax=xmax;
    
96.     model.ymin=ymin;
    
97.     model.ymax=ymax;
    
98. end
    
99. 
    
100. function L=TourLength(tour,model)
     
101.     n=numel(tour);
     
102.     tour=[tour tour(1)];
     
103.     L=0;
     
104.     for k=1:n
     
105.         i=tour(k);
     
106.         j=tour(k+1);
     
107.         L=L+model.d(i,j);
     
108.     end
     
109. end
     
110. 
     
111. function ActionList=CreatePermActionList(n)
     
112.     nSwap=n*(n-1)/2;
     
113.     nReversion=n*(n-1)/2;
     
114.     nInsertion=n^2;
     
115.     nAction=nSwap+nReversion+nInsertion;
     
116.     ActionList=cell(nAction,1);
     
117.     c=0;
     
118.     % Add SWAP
     
119.     for i=1:n-1
     
120.         for j=i+1:n
     
121.             c=c+1;
     
122.             ActionList{c}=[1 i j];
     
123.         end
     
124.     end
     
125.     % Add REVERSION
     
126.     for i=1:n-1
     
127.         for j=i+1:n
     
128.             if abs(i-j)>2
     
129.                 c=c+1;
     
130.                 ActionList{c}=[2 i j];
     
131.             end
     
132.         end
     
133.     end
     
134.     % Add Insertion
     
135.     for i=1:n
     
136.         for j=1:n
     
137.             if abs(i-j)>1
     
138.                 c=c+1;
     
139.                 ActionList{c}=[3 i j];
     
140.             end
     
141.         end
     
142.     end
     
143.     ActionList=ActionList(1:c);
     
144. end
     
145. 
     
146. function q=DoAction(p,a)
     
147.     switch a(1)
     
148.         case 1
     
149.             % Swap
     
150.             q=DoSwap(p,a(2),a(3));
     
151.         case 2
     
152.             % Reversion
     
153.             q=DoReversion(p,a(2),a(3));
     
154.         case 3
     
155.             % Insertion
     
156.             q=DoInsertion(p,a(2),a(3));     
     
157.     end
     
158. end
     
159. 
     
160. function q=DoSwap(p,i1,i2)
     
161.     q=p;
     
162.     q([i1 i2])=p([i2 i1]);
     
163. end
     
164. 
     
165. function q=DoReversion(p,i1,i2)
     
166.     q=p;
     
167.     if i1<i2
     
168.         q(i1:i2)=p(i2:-1:i1);
     
169.     else
     
170.         q(i1:-1:i2)=p(i2:i1);
     
171.     end
     
172. end
     
173. 
     
174. function q=DoInsertion(p,i1,i2)
     
175.     if i1<i2
     
176.         q=p([1:i1-1 i1+1:i2 i1 i2+1:end]);
     
177.     else
     
178.         q=p([1:i2 i1 i2+1:i1-1 i1+1:end]);
     
179.     end
     
180. end

复制代码

why2018470

谢谢楼主分享