% 遺傳算法在TSP中的應用
function [s,q]=yichuan(D,n)       % D為距離矩陣，難為TSP中的城市個數
tic
fprintf('注意：種群中個體數、遺傳代數越多，需要運行的時間越長,但是結果越好！\n');
m=input('請輸入種群大小,最好為偶數m=');  
p=input('請輸入遺傳代數p=');
% 產生初始種群
for i=1:m
    e=1:n;
    for j=1:n
        t=unidrnd(length(e));
        q(i,j)=e(t);
        e(t)=[];
    end
end
s=len(q,D,n,m); % 計算路程長度
for i=1:p
    q=newq(s,q,m);% 產生新種群
    q=change(q,m,n); % 交叉
    q=vf(q,n,m);   % 變異
    s=len(q,D,n,m); % 計算路程長度
end
[s0,b]=min(s);
fprintf('中總路程為s=%4.2f\n',s0);
fprintf('站點順序為\n');
fprintf('%d  ',q(b,:));
toc


% 計算路程即目標函數
function s=len(q,D,n,m)
for i=1:m
    s(i)=0;
    for j=1:n-1
        s(i)=s(i)+D(q(i,j),q(i,j+1));
    end
end


%  選擇目標函數小的個體，產生新種群
function q=newq(s,q,m)
t=0.20;  % 淘汰比例
% 降序排列
for i=1:m-1
    for j=i:m
        if s(i)<s(j)
            k=s(j);
            s(j)=s(i);
            s(i)=k;
            
            Q=q(j,:);
            q(j,:)=q(i,:);
            q(i,:)=Q;
        end
    end
end
f=fix(m*t);  % 淘汰個數
% 將目標函數值大的前f個個體用目標函數值小的后f個個體
for i=1:f
    q(i,:)=q(m-f+i,:);
end
% 打亂順序
for i=1:2*m
    N1=unidrnd(m);
    N2=unidrnd(m);
    Q=q(N1,:);
    q(N1,:)=q(N2,:);
    q(N2,:)=Q;
end



%  兩兩交叉操作
function q=change(q,m,n)
while 1
    N1=unidrnd(n-1);
    N2=unidrnd(n-1);
    if N2>N1+1 
        break;
    end
end

%  采用部分匹配交叉，且一對一對交叉
for i=1:2:m-1
    % 中間部分交換
    Q=q(i,N1+1:N2);
    q(i,N1+1:N2)=q(i+1,N1+1:N2);
    q(i+1,N1+1:N2)=Q;
    % 選出可以保留的部分
    t1=1;t2=1;
    for j=1:n        
        if (q(i,j)~=q(i,N1+1:N2))
            e1(t1)=j;
            t1=t1+1;
        end
        if (j>N1&j<=N2)
            e1(t1)=j;
            t1=t1+1;
        end
       if (q(i+1,j)~=q(i+1,N1+1:N2))
           e2(t2)=j;
           t2=t2+1;
       end
       if (j>N1&j<=N2)
           e2(t2)=j;
           t2=t2+1;
       end
    end
    % 替換重復部分
    for j=1:n
        if j~=e1
            for k=N1+1:N2
                if q(i+1,k)~=q(i,e1)
                    q(i,j)=q(i+1,k);
                    e1(t1)=j;
                    t1=t1+1;
                    break;
                end
            end
        end
        if j~=e2
           for k=N1+1:N2
                if q(i,k)~=q(i+1,e2)
                    q(i+1,j)=q(i,k);
                    e2(t2)=j;
                    t2=t2+1;
                    break;
                end
            end 
        end
    end
    for j=1:n
        if j~=e1
            for k=1:n
                if k~=q(i,e1)
                    q(i,j)=k;
                    e1(t1)=j;
                    t1=t1+1;
                    break;
                end
            end
        end
        if j~=e2
            for k=1:n
                if k~=q(i+1,e2)
                    q(i+1,j)=k;
                    e2(t2)=j;
                    t2=t2+1;
                    break;
                end
            end
        end 
    end
    clear e1 e2
end


%  變異
function q=vf(q,n,m)
F=0.01 ; %以1%的概率變異
N1=unidrnd(m);   % 種群中第N1個個體可能變異
N2=1+unidrnd(n-1); % 編碼中第N2個號可能變異
k=1+unidrnd(n-1);  % q(n1,n2)變為k
p0=unifrnd(0,1);   % 隨機產生變異概率
if p0<F
    for i=1:n
        if q(N1,i)==k
            q(N1,i)=q(N1,N2);
            break;
        end
    end
    q(N1,N2)=k;
end