禁忌搜索算法的TSP问题求解
禁忌搜索算法的TSP问题求解：
clc,clear,close all
warning off
x=;
y=;
= Tabu_Search_TSP(x,y,50);
figure(1);
hold on
plot(x,y,'o')
for i=1:length(x)
if i<length(x)
plot( , ,'-','linewidth',2 )
else
plot( , ,'-','linewidth',2 )
end
end
hold off函数如下：
function = Tabu_Search_TSP(x,y,MaxIt)
% x=;
% y=;
model = CreateModel_ysw(x, y); % Create TSP Model
CostFunction=@(tour) TourLength(tour, model); % Cost Function
ActionList=CreatePermActionList(model.n); % Action List
nAction=numel(ActionList); % Number of Actions
%% Tabu Search Parameters
% MaxIt=50; % Maximum Number of Iterations
TL=round(0.5*nAction); % Tabu Length
%% Initialization
% Create Empty Individual Structure
empty_individual.Position=[];
empty_individual.Cost=[];
% Create Initial Solution
sol=empty_individual;
sol.Position=randperm(model.n);
sol.Cost=CostFunction(sol.Position);
% Initialize Best Solution Ever Found
BestSol=sol;
% Array to Hold Best Costs
BestCost=zeros(MaxIt,1);
% Initialize Action Tabu Counters
TC=zeros(nAction,1);
%% Tabu Search Main Loop
for it=1:MaxIt
bestnewsol.Cost=inf;
% Apply Actions
for i=1:nAction
if TC(i)==0
newsol.Position=DoAction(sol.Position,ActionList{i});
newsol.Cost=CostFunction(newsol.Position);
newsol.ActionIndex=i;
if newsol.Cost<=bestnewsol.Cost
bestnewsol=newsol;
end
end
end
% Update Current Solution
sol=bestnewsol;
% Update Tabu List
for i=1:nAction
if i==bestnewsol.ActionIndex
TC(i)=TL; % Add To Tabu List
else
TC(i)=max(TC(i)-1,0); % Reduce Tabu Counter
end
end
% Update Best Solution Ever Found
if sol.Cost<=BestSol.Cost
BestSol=sol;
end
% Save Best Cost Ever Found
BestCost(it)=BestSol.Cost;
% Show Iteration Information
% disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
% % Plot Best Solution
% figure(1);
% PlotSolution(BestSol, model);
% pause(0.01);
% If Global Minimum is Reached
if BestCost(it)==0
break;
end
end
% BestCost=BestCost(1:it);
Position = BestSol.Position;
BestCost = BestSol.Cost;
end
function model=CreateModel_ysw(x,y)
% x=;
% y=;
n=numel(x);
d=zeros(n,n);
for i=1:n-1
for j=i+1:n
d(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
d(j,i)=d(i,j);
end
end
xmin=0;
xmax=100;
ymin=0;
ymax=100;
model.n=n;
model.x=x;
model.y=y;
model.d=d;
model.xmin=xmin;
model.xmax=xmax;
model.ymin=ymin;
model.ymax=ymax;
end
function L=TourLength(tour,model)
n=numel(tour);
tour=;
L=0;
for k=1:n
i=tour(k);
j=tour(k+1);
L=L+model.d(i,j);
end
end
function ActionList=CreatePermActionList(n)
nSwap=n*(n-1)/2;
nReversion=n*(n-1)/2;
nInsertion=n^2;
nAction=nSwap+nReversion+nInsertion;
ActionList=cell(nAction,1);
c=0;
% Add SWAP
for i=1:n-1
for j=i+1:n
c=c+1;
ActionList{c}=;
end
end
% Add REVERSION
for i=1:n-1
for j=i+1:n
if abs(i-j)>2
c=c+1;
ActionList{c}=;
end
end
end
% Add Insertion
for i=1:n
for j=1:n
if abs(i-j)>1
c=c+1;
ActionList{c}=;
end
end
end
ActionList=ActionList(1:c);
end
function q=DoAction(p,a)
switch a(1)
case 1
% Swap
q=DoSwap(p,a(2),a(3));
case 2
% Reversion
q=DoReversion(p,a(2),a(3));
case 3
% Insertion
q=DoInsertion(p,a(2),a(3));
end
end
function q=DoSwap(p,i1,i2)
q=p;
q()=p();
end
function q=DoReversion(p,i1,i2)
q=p;
if i1<i2
q(i1:i2)=p(i2:-1:i1);
else
q(i1:-1:i2)=p(i2:i1);
end
end
function q=DoInsertion(p,i1,i2)
if i1<i2
q=p();
else
q=p();
end
end
谢谢楼主分享
页:
[1]