【GreenSim.C原创】TSP问题遗传模拟退火算法 2007-03-10 10:09:27
大 中 小
function [R_best,L_best,LC,ALC]=TSP_GSA(C,M,S,Age,T0,P,K,alpha)
%%=========================================================================
%% TSP_GSA.m
%% 旅行商问题遗传模拟退火算法
%% ChengAihua,PLA Information Engineering University,ZhengZhou,China
%% Email:aihuacheng@gmail.com
%% All rights reserved
%%-------------------------------------------------------------------------
%% 输入参数：
%% C Coordinate 节点坐标，由一个N×2的矩阵存储
%% M Adjacent Matrix 赋权邻接矩阵，由一个N×N的矩阵存储
%% S Size 种群的规模
%% Age 进化代数
%% T0 初始温度
%% P Proportion 等比降温的降温系数
%% K 状态转移次数
%% alpha 竞争力；两级分化调节参数
%% 输出参数：
%% R_best 最优路径
%% L_best 最优路径对应的路线长度
%% LC Learning Ｃurve 学习曲线（记录各代最优路径对应的路线长度）
%% ALC Average Learning Ｃurve 记录各代种群的平均路线长度
%% 主要中间变量：
%% N 问题规模，城市个数，基因长度
%% R Route 路线
%% L Length 路线长度
%% farm 种群，S×N的矩阵
%% FARM 扩容后的种群
%% Tem 温度
%% counter 代数
%%=========================================================================
%%-------------------------------------------------------------------------
%% 随机产生初始种群
N=size(M,1);
farm=zeros(S,N);
for ii=1:S
farm(ii,:)=randperm(N);
end
%%-------------------------------------------------------------------------
%% 变量初始化
counter=1;%进化代数计数器置零
Tem=T0;%温度置为初始温度
R_best=zeros(1,N);
L_best=0;
LC=zeros(1,Age);
ALC=LC;
%%-------------------------------------------------------------------------
%% 最外层循环（进化）开始
while counter<=Age
%%---------------------------------------------------------------------
%% 交叉
W=ceil(N/6);
if mod(S,2)==1
farm1=[farm;farm(unidrnd(S-2)+1,:)];
farm2=farm1;
for ii=1:2:S
A=farm1(ii,:);
B=farm1(ii+1,:);
[a,b]=PMX_intercross(A,B,N,W);%调用部分匹配交叉算法
farm2(ii,:)=a;
farm2(ii+1,:)=b;
end
FARM=[farm1(1:S,:);farm(1:S,:)];
elseif mod(S,2)==0
farm1=farm;
farm2=farm1;
for ii=1:2:(S-1)
A=farm(ii,:);
B=farm(ii+1,:);
[a,b]=PMX_intercross(A,B,N,W);%调用部分匹配交叉算法
farm2(ii,:)=a;
farm2(ii+1,:)=b;
end
FARM=[farm1;farm2];
else
end
%打乱新种群中个体排列秩序，保证随机组合交配
Order=randperm(2*S);
newFARM=FARM;
for ii=1:(2*S)
newFARM(ii,:)=FARM(Order(ii),:);
end
FARM=newFARM;
%%---------------------------------------------------------------------
%% 状态转移（通过变异产生新解和对新解按照概率接受准则进行处理）
W_up=ceil(N/20);
for ii=1:(2*S)
R_old=FARM(ii,:);
R_new=State_jump(R_old,K,W_up,Tem,M);%调用状态转移子程序
FARM(ii,:)=R_new;
end
%%---------------------------------------------------------------------
%% 计算新群体的适应值，记录输出量
Fit=zeros(2*S,1);
for ii=1:(2*S)
Fit(ii,1)=Compute_length(FARM(ii,:),M);
end
L_best=min(Fit);
LC(counter)=L_best;
ALC(counter)=mean(Fit);
P_best=find(Fit==L_best);
R_best=FARM(P_best(1),:);
%%---------------------------------------------------------------------
%% 含竞争力两极分化程度调节的选择复制操作
minL=min(Fit);
maxL=max(Fit);
Fit=(maxL.*ones(2*S,1)-Fit)/(maxL-minL);
Fit=Fit./(sum(Fit));
Fit=Fit.^alpha;
Fit=Fit./(sum(Fit));
%构造转轮赌
Wh=zeros(1,2*S);
for ii=1:(2*S);
aa=Fit(1:ii);
Wh(ii)=sum(aa);
end
%随机选择下一代个体
Flag=zeros(1,2*S);
Flag(P_best(1))=1;%保存最优个体
sumFlag=sum(Flag);
while sumFlag<S
Pos=find(Wh>=rand);
Flag(Pos(1))=1;
sumFlag=sum(Flag);
end
kk=0;
for ii=1:(2*S)
if Flag(ii)==1
kk=kk+1;
farm(kk,:)=FARM(ii,:);
end
end
farm(1,:)=R_best;
%%---------------------------------------------------------------------
%% 降温进化
counter=counter+1
Tem=Tem*P;
end
%%-------------------------------------------------------------------------
%% 绘图输出结果
subplot(1,2,1)
Draw_route(C,R_best)
subplot(1,2,2)
plot(LC)
hold on
plot(ALC)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%=========================================================================
%% PMX_intercross.m
%% 部分匹配交叉算法
%% PMX（部分匹配交叉）由Goldberg和Lingle于1985年提出
%%-------------------------------------------------------------------------
%% A 父代A（1×N的行向量，下同）
%% B 父代B
%% a 子代a
%% b 子代b
%% N 基因长度
%% W 交叉宽度（假设已在本程序外部随机取定）
%%=========================================================================
function [a,b]=PMX_intercross(A,B,N,W)
%在环状基因链上独立选择交叉点
p1=unidrnd(N);
if p1~=1
A=[A(p1:N),A(1:p1-1)];
end
p2=unidrnd(N);
if p2~=1
B=[B(p2:N),B(1:p2-1)];
end
%部分匹配交叉
a=[B(1:W),A(W+1:N)];
b=[A(1:W),B(W+1:N)];
f1=zeros(1,N);
f2=f1;
for ii=W+1:N
if length(find(a==a(ii)))==2
f1(ii)=1;
end
if length(find(b==b(ii)))==2
f2(ii)=1;
end
end
ff1=find(f1==1);
ff2=find(f2==1);
for ii=1:length(ff1)
temp=a(ff1(ii));
a(ff1(ii))=b(ff2(ii));
b(ff2(ii))=temp;
end
%恢复
if p1~=1
a=[a(N-p1+2:N),a(1:N-p1+1)];
end
if p2~=1
b=[b(N-p2+2:N),b(1:N-p2+1)];
end
%%=========================================================================
%% State_jump.m
%% 状态转移（通过变异产生新解和对新解按照概率接受准则进行处理）
%%-------------------------------------------------------------------------
%% R_old 状态转移前的旧状态
%% R_new 经状态转移达到的新状态
%% W_up 随机变异位宽上限（假设已在本程序外部随机取定）
%% Tem 温度
%% K 状态转移次数
%% M Adjacent Matrix 赋权邻接矩阵，由一个N×N的矩阵存储
%%=========================================================================
function R_new=State_jump(R_old,K,W_up,Tem,M)
N=length(R_old);
k=1;
while k<=K
%随机选择变异开始点和变异宽度
p=unidrnd(N);
if p~=1
R_old=[R_old(p:N),R_old(1:p-1)];
end
R_new=R_old;
W_mut=unidrnd(W_up-1)+1;
%采用变异局部打乱排列秩序的方法变异（有待改进成启发式变异，以提高搜索效率）
Order=randperm(W_mut);
for ii=1:W_mut
R_new(ii)=R_old(Order(ii));
end
if p~=1
R_old=[R_old(N-p+2:N),R_old(1:N-p+1)];
R_new=[R_new(N-p+2:N),R_new(1:N-p+1)];
end
%概率接受准则
Fit_old=Compute_length(R_old,M);
Fit_new=Compute_length(R_new,M);
if Fit_new<=Fit_old
R_old=R_new;
elseif exp((Fit_old-Fit_new)/Tem)>rand
R_old=R_new;
else
end
k=k+1;
end
%%=========================================================================
%% Compute_length.m
%% 计算路径长度
%%-------------------------------------------------------------------------
%% M Adjacent Matrix 赋权邻接矩阵，由一个N×N的矩阵存储
%% R Route 路线
%% L Length 路线长度
%%=========================================================================
function L=Compute_length(R,M)
N=size(M,1);
L=M(R(N),R(1));
for ii=1:(N-1)
L=L+M(R(ii),R(ii+1));
end
%%=========================================================================
%% Draw_route.m
%% 画路线图
%%-------------------------------------------------------------------------
%% C Coordinate 节点坐标，由一个N×2的矩阵存储
%% R Route 路线
%%============================================================
