% 基于MATLAB的遺傳算法程序用于解決TSP難題
 
% 傳統遺傳算法：

M=36;popsize=36;Pm=0.15;Pc=0.9;Po=0.75;T=5000;Samen=500;
tic;
Smax=zeros(1,T+1);
num=0;
V=initiate(popsize,M,D);
S=fitness(popsize,M,V,D);
Smax(1)=max(S);
for i=1:T
    C=0;i
    V1=cross(V,popsize,M,Pc);
    V1=mutate(V1,popsize,M,Pm);
    V1=opt(V1,popsize,M,Po,D);
    S1=fitness(popsize,M,V1,D);
    S1m=mean(S1);
    G=find(S>=S1m);
    Vp=V(G,:);
    W=[V1;Vp];
    a=size(W,1);
    for j=1:a
        if W(j,1)==0
            continue;
        end
        for k=j+1:a
            if W(k,:)==W(j,:)
                W(k,1)=0;
            end
        end
    end
    W(find(W(:,1)==0),:)=[];
    d=size(W,1);
    Sw=fitness(d,M,W,D);
    [Smax(i+1),k]=max(Sw);
    t=W(k,:);Sw(k)=[];
    if Smax(i+1)==Smax(i)
        num=num+1;
    else 
        num=0;
    end
    if num==Samen
        break;
    end
    P=Sw./sum(Sw);
    Q=cumsum(P);
    R=rand(1,popsize-1);
    for j=1:popsize-1
        for k=1:d-1
            if R(j)<=Q(k)
                A(j)=k;
                break;
            end
        end
    end
    V=[t;W(A,:)];
    S=fitness(popsize,M,V,D);
    SNmax=max(S);
    for j=1:popsize
        if S(j)>=0.8*SNmax
            C=C+1;
        end
    end
    if C/popsize>=1;
        Pm=5*Pm;
    end
end
time=toc;Smin=1./Smax(1:i+1);
plot(1:i+1,Smin);
disp('最佳路線為'),genebest=t,
disp('其總路程為'),1/Smax(i+1),
disp('共用時'),time
disp('進化代數'),i,
[SMin,k]=min(Smin)

% 并行遺傳算法：

M=36;popsize=36;Pm=0.15;Pc=0.9;Po=0.75;T=5000;Samen=100;
tic;
Smax=zeros(1,T+1);
num=0;
V=initiate(popsize,M,D);
S=fitness(popsize,M,V,D);
Smax(1)=max(S);
for i=1:T
    C=0;
    Vc=cross(V,popsize,M,Pc);
    Sc=fitness(popsize,M,Vc,D);
    Smc=mean(Sc);Sc=[];
    Vm=mutate(V,popsize,M,Pm);
    Sm=fitness(popsize,M,Vm,D);
    Smm=mean(Sm);Sm=[];
    Vo=opt(V,popsize,M,Po,D);
    Vp=parent(V,S,Smc,Smm);
    W=[Vc;Vm;Vp;Vo];Vc=[];Vm=[];Vp=[];Vo=[];
    a=size(W,1);
    for j=1:a
        if W(j,1)==0
            continue;
        end
        for k=j+1:a
            if W(k,:)==W(j,:)
                W(k,1)=0;
            end
        end
    end
    W(find(W(:,1)==0),:)=[];
    d=size(W,1);
    Sw=fitness(d,M,W,D);
    [Smax(i+1),k]=max(Sw);
    t=W(k,:);Sw(k)=[];
    if Smax(i+1)==Smax(i)
        num=num+1;
    else num=0;
    end
    P=Sw./sum(Sw);
    Q=cumsum(P);P=[];
    R=rand(1,popsize-1);
    for j=1:popsize-1
        for k=1:d-1
            if R(j)<=Q(k)
                A(j)=k;
                break;
            end
        end
    end
    V=[t;W(A,:)];W=[];Q=[];R=[];
    S=fitness(popsize,M,V,D);SNmax=max(S);
    for j=1:popsize
        if S(j)>=0.8*SNmax
            C=C+1;
        end
    end
    if C/popsize>=1;
        Pm=5*Pm;
    end
end
time=toc;Smin=1./Smax(1:i+1);
plot(1:i+1,Smin);
disp('最佳路線為'),genebest=t,
disp('其總路程為'),Smin(i+1),
disp('共用時'),time
disp('進化代數'),i,
[SMin,k]=min(Smin)


% 貪婪算子初始化群體：

function V=initiate(popsize,M,D)
d=D;
for i=1:M
    d(i,i)=inf;
end
for t=1:popsize
    do=d;i=t;V(t,1)=t;
    for j=1:M-1
        [dmin,k]=min(do(i,:));
        do(1:M,i)=inf*ones(M,1);
        i=k;
        V(t,j+1)=k;
    end
end


% 雜交算子：

function Vc=cross(V,popsize,M,Pc)
R1=rand(1,popsize); % selectout font
a=find(R1<Pc);
b=size(a,2);
R2=randn(1,b);
for i=1:b-1
    for j=1:b-i
        if R2(j)<R2(j+1)
            t1=R2(j);R2(j)=R2(j+1);R2(j+1)=t1;
            t2=a(j);a(j)=a(j+1);a(j+1)=t2;
        end
    end
end 
R3=rand(floor(b/2),M);
[maxr,k1]=max(R3,[],2);
[minr,k2]=min(R3,[],2);
A=[k1,k2];
for i=1:floor(b/2)
    if A(i,2)<A(i,1)
        t1=A(i,1);A(i,1)=A(i,2);A(i,2)=t1;
    end
end % selectout end
Vc=V;
for i=1:floor(b/2)
    t3=Vc(a(2*i-1),[A(i,1):A(i,2)]);
    t4=Vc(a(2*i),A(i,1):A(i,2));
    if (A(i,2)==M)&(A(i,1)==1)
        continue;
    end
    C1=(A(i,1)~=1);C2=(A(i,2)~=M);
    t5=[C2*(A(i,2)+1:M),C1*(1:A(i,1)-1)];
    Vc1=[Vc(a(2*i),t5),t4];
    Vc12=[Vc(a(2*i-1),t5),t3];
    for j=1:A(i,2)-A(i,1)+1
        ao=find(Vc1==t3(j));
        a1=find(Vc12==t4(j));
        Vc1(ao)=[];
        Vc12(a1)=[];
    end 
    Vc(a(2*i-1),:)=[Vc1(M-A(i,2)+1:end),t3,Vc1(1:M-A(i,2))];
    Vc(a(2*i),:)=[Vc12(M-A(i,2)+1:end),t4,Vc12(1:M-A(i,2))];
end


% 變異算子：

function Vm=mutate(V,popsize,M,Pm)
R1=rand(popsize,M);
Vm=V;
for i=1:popsize
    a=find(R1(i,:)<Pm);
    b=size(a,2);
    if b==0
        continue;
    end
    if b==1
        Vm1=Vm(i,:);
        t=Vm1(a);
        Vm1(a)=[];
        R2=rand(1,M-1);
        [Rmax,k]=max(R2);
        Vm(i,:)=[Vm1(1:k-1),t,Vm1(k:M-1)];
    else
        R2=rand(1,b);
        for j=1:b-1
            for k=1:b-j
                if R2(k)<R2(k+1)
                    t1=R2(k);R2(k)=R2(k+1);R2(k+1)=t1;
                    t2=Vm(i,a(k));Vm(i,a(k))=Vm(i,a(k+1));Vm(i,a(k+1))=t2;
                end
            end
        end
    end
end


% 父代較優個體：

function Vp=parent(V,S,Smc,Smm)
Sp=Smc;
if Smm>Smc
    Sp=Smm;
end
a=find(S>=Sp);
Vp=V(a,:);

% 兩點法局部優化算子：

function Vo=opt(V,popsize,M,Po,D)
R=rand(1,popsize);
a=find(R<Po);
b=size(a,2);
Vo=V;
for k=1:b
    for i=1:M-3
        for j=i+2:M-1
            if D(i,j)+D(i+1,j+1)<=D(i,i+1)+D(j,j+1)
                Vo(a(k),i+1:j)=Vo(a(k),j:-1:i+1);
            end
        end
    end
end


% 計算適應值：

function S=fitness(L,M,V,D)
S=zeros(1,L);
V1=[V(:,2:M),V(:,1)];
for i=1:L
    for j=1:M
        S(i)=S(i)+D(V(i,j),V1(i,j));
    end
end
S=1./S;