tsp问题遗传模拟退火算法.txt

TSP问题遗传退火算法

【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:[e-mail]
%%      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   路线
%%=========================================================================
function Draw_route(C,R)
N=length(R);
scatter(C(:,1),C(:,2));
hold on
plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)])
hold on
for ii=2:N
    plot([C(R(ii-1),1),C(R(ii),1)],[C(R(ii-1),2),C(R(ii),2)])
    hold on
end