遺傳算法解決TSP問題
遺傳算法
遺傳算法的基本原理是通過作用於染色體上的基因尋找好的染色體來求解問題，它需要對算法所產生的每個染色體進行評價，並基於適應度值來選擇染色體，使適應性好的染色體有更多的繁殖機會，在遺傳算法中，通過隨機方式產生若干個所求解問題的數字編碼，即染色體，形成初始種群；通過適應度函數給每個個體一個數值評價，淘汰低適應度的個體，選擇高適應度的個體參加遺傳操作，經過遺產操作后的個體集合形成下一代新的種群，對這個新的種群進行下一輪的進化。
TSP問題
TSP問題即旅行商問題，經典的TSP可以描述為：一個商品推銷員要去若干個城市推銷商品，該推銷員從一個城市出發，需要經過所有城市后，回到出發地。應如何選擇行進路線，以使總的行程最短。從圖論的角度來看，該問題實質是在一個帶權完全無向圖中，找一個權值最小的哈密爾頓回路。
遺傳算法解決TSP問題
概念介紹：
種群 ==> 可行解集
個體 ==> 可行解
染色體 ==> 可行解的編碼
基因 ==> 可行解編碼的分量
基因形式 ==> 遺傳編碼
適應度 ==> 評價的函數值（適應度函數）
選擇 ==> 選擇操作
交叉 ==> 編碼的交叉操作
變異 ==> 可行解編碼的變異
遺傳操作：就包括優選適應性強的個體的“選擇”；個體間交換基因產生新個體的“交叉”；個體間的基因突變而產生新個體的“變異”。其中遺傳算法是運用遺傳算子來進行遺傳操作的。即：選擇算子、變異算子、交叉算子。
遺傳算法的基本運算過程
（1）種群初始化：個體編碼方法有二進制編碼和實數編碼，在解決TSP問題過程中個體編碼方法為實數編碼。對於TSP問題，實數編碼為1-n的實數的隨機排列，初始化的參數有種群個數M、染色體基因個數N（即城市的個數）、迭代次數C、交叉概率Pc、變異概率Pmutation。
（2）適應度函數：在TSP問題中，對於任意兩個城市之間的距離D(i,j)已知，每個染色體（即n個城市的隨機排列）可計算出總距離，因此可將一個隨機全排列的總距離的倒數作為適應度函數，即距離越短，適應度函數越好，滿足TSP要求。
（3）選擇操作：采用基於適應度比例的選擇策略，即適應度越好的個體被選擇的概率越大，同時在選擇中保存適應度最高的個體。
（4）交叉操作：對於個體，隨機選擇兩個個體，在對應位置交換若干個基因片段，同時保證每個個體依然是1-n的隨機排列，防止進入局部收斂。
例如：
對A的4 3 2進行和B的交叉時，將B的2 6 7換至A，同時將A中的6 7（B換上來但A已包含）和B中對應的數字進行交叉，6->3 7->2
（5）變異操作：對於變異操作，隨機選取個體，同時隨機選取個體的兩個基因進行交換以實現變異操作。
產生新的群體后再次進行評估，然后再選擇、交叉、變異，一直循環這幾步，最終找到一個近似最優解。
遺傳算法的終止條件
- 當最優個體的適應度達到給定的閾值
- 最優個體的適應度和群體適應度不再上升
- 迭代次數達到預設的代數時，算法終止
運行結果：
第一次測試
N=25; %城市的個數
M=100; %種群的個數
TER=2000; %迭代次數
Pc=0.8; %%交叉概率
Pmutation=0.05; %%變異概率
迭代2000次之后運行結果並不能算是很好，整體迭代了600次左右趨於平穩，但仍有部分細微波動，運行時間相較於蟻群算法更快。
第二次測試
N=25; %城市的個數
M=100; %種群的個數
TER=2000; %迭代次數
Pc=0.95; %交叉概率
Pmutation=0.1; %變異概率
提高了變異概率和交叉概率，雖然最終結果的准確性有略微提升，但是犧牲了運行時間和迭代次數，迭代次數遠遠大於第一次測試，可見變異概率和交叉概率過大對實驗結果存在一定影響。
第三次測試
N=25; %城市的個數
M=100; %種群的個數
TER=2000; %迭代次數
Pc=0.4; %交叉概率
Pmutation=0.01; %變異概率
結合第二次測試可以看出，交叉概率和變異概率不能過大，否則波動較大，迭代時間較長，同樣也不能過小，否則很難找到最優解
第四次測試
N=25; %城市的個數
M=200; %種群的個數
TER=2000; %迭代次數
Pc=0.8; %交叉概率
Pmutation=0.05; %變異概率
種群的個數增加成原來的一倍后，運行的結果更加精確，但是相應的迭代次數會增加，時間成本增加
第五次測試
N=25; %城市的個數
M=500; %種群的個數
TER=2000; %迭代次數
Pc=0.8; %交叉概率
Pmutation=0.05; %變異概率
可以發現當種群個數增加更多之后，迭代次數大幅減少，但是運行的時間與第四次測試相比增加了近20秒，時間成本增加了很多，所以種群個數同樣不應該過大。
總結
對同一個TSP問題，分析種群規模、交叉概率和變異概率對算法結果的影響：
（1）改變種群個數后的影響-->種群個數增大，算法結果的精確度會有一定的提升，但是運行時間相較之前更久了。
（2）改變交叉概率后的影響-->交叉概率過低將很難得不到最優解，交叉概率越高則平均適應度越好，但是也不能過高，否則會影響迭代時間。
（3）改變變異概率概率后的影響-->變異概率過高或者過低都會影響運行得到最優解
（4）遺傳算法得到的結果的精確度會受到交叉概率、變異概率、迭代次數的影響：迭代次數越大，變異概率越小，則遺傳算法的精確度越高。執行時間隨着迭代次數的增加而增加。當交叉概率為0.8，變異概率為0.5，結果相對來說會稍微好一些，整體來說運用遺傳算法來求解TSP問題的速度要比用蟻群算法來解決TSP問題快的多，所以適當將迭代次數增加一些也沒有什么影響。
代碼
%main clear; clc; %%%%%%%%%%%%%%%輸入參數%%%%%%%% N=25; %%城市的個數 M=100; %%種群的個數 ITER=2000; %%迭代次數 %C_old=C; m=2; %%適應值歸一化淘汰加速指數 Pc=0.8; %%交叉概率 Pmutation=0.05; %%變異概率 %%生成城市的坐標 pos=[-0.0596300734351115 -0.851297433985992; 0.438398458335735 0.841724248850041; -0.251804207069105 -1.30744728793171; 1.17790692164388 0.0952195708913433; 1.66427009988301 -0.133382476901910; -0.337308956003599 -0.407225198007313; 0.311862707690034 -0.431724333072457; 0.182396732713558 -0.157788333124223; 0.0378390295101906 -0.367920724129644; -0.853709951356673 0.0363449796675028; -1.64191242274810 -0.565028275383299; -0.821481309108439 -1.36166562718530; -0.0666502703726364 0.756667745411741; 2.25931856568426 -1.30140683361980; 0.759098750479421 2.81038357991309; -0.00925105144207478 0.672464995693229; 0.862239245550680 -0.576524359942781; -0.742512628202077 0.740723643506810; 0.371005938734921 -0.835740815137822; 0.677683424862960 -0.284758907693345; -1.48234180628974 0.955162776644004; 0.452806831664590 0.608019647737753; 1.22196689592397 -0.878318328597688; 0.352173262723895 1.60512181080303; 1.80729484625293 0.593061432086341 ]; %%生成城市之間距離矩陣 D=zeros(N,N); for i=1:N for j=i+1:N dis=(pos(i,1)-pos(j,1)).^2+(pos(i,2)-pos(j,2)).^2; D(i,j)=dis^(0.5); D(j,i)=D(i,j); end end %%生成初始群體 popm=zeros(M,N); for i=1:M popm(i,:)=randperm(N);%隨機排列，比如[2 4 5 6 1 3] end %%隨機選擇一個種群 R=popm(1,:); figure(1); scatter(pos(:,1),pos(:,2),'rx');%畫出所有城市坐標 axis([-3 3 -3 3]); figure(2); plot_route(pos,R); %%畫出初始種群對應各城市之間的連線 axis([-3 3 -3 3]); %%初始化種群及其適應函數 fitness=zeros(M,1); len=zeros(M,1); for i=1:M%計算每個染色體對應的總長度 len(i,1)=myLength(D,popm(i,:)); end maxlen=max(len);%最大回路 minlen=min(len);%最小回路 fitness=fit(len,m,maxlen,minlen); rr=find(len==minlen);%找到最小值的下標，賦值為rr R=popm(rr(1,1),:);%提取該染色體，賦值為R for i=1:N fprintf('%d ',R(i));%把R順序打印出來 end fprintf('\n'); fitness=fitness/sum(fitness); distance_min=zeros(ITER+1,1); %%各次迭代的最小的種群的路徑總長 nn=M; iter=0; while iter<=ITER fprintf('迭代第%d次\n',iter); %%選擇操作 p=fitness./sum(fitness); q=cumsum(p);%累加 for i=1:(M-1) len_1(i,1)=myLength(D,popm(i,:)); r=rand; tmp=find(r<=q); popm_sel(i,:)=popm(tmp(1),:); end [fmax,indmax]=max(fitness);%求當代最佳個體 popm_sel(M,:)=popm(indmax,:); %%交叉操作 nnper=randperm(M); % A=popm_sel(nnper(1),:); % B=popm_sel(nnper(2),:); %% for i=1:M*Pc*0.5 A=popm_sel(nnper(i),:); B=popm_sel(nnper(i+1),:); [A,B]=cross(A,B); % popm_sel(nnper(1),:)=A; % popm_sel(nnper(2),:)=B; popm_sel(nnper(i),:)=A; popm_sel(nnper(i+1),:)=B; end %%變異操作 for i=1:M pick=rand; while pick==0 pick=rand; end if pick<=Pmutation popm_sel(i,:)=Mutation(popm_sel(i,:)); end end %%求適應度函數 NN=size(popm_sel,1); len=zeros(NN,1); for i=1:NN len(i,1)=myLength(D,popm_sel(i,:)); end maxlen=max(len); minlen=min(len); distance_min(iter+1,1)=minlen; fitness=fit(len,m,maxlen,minlen); rr=find(len==minlen); fprintf('minlen=%d\n',minlen); R=popm_sel(rr(1,1),:); for i=1:N fprintf('%d ',R(i)); end fprintf('\n'); popm=[]; popm=popm_sel; iter=iter+1; %pause(1); end %end of while figure(3) plot_route(pos,R); axis([-3 3 -3 3]); figure(4) plot(distance_min);
%交叉操作函數 cross.m function [A,B]=cross(A,B) L=length(A); if L<10 W=L; elseif ((L/10)-floor(L/10))>=rand&&L>10 W=ceil(L/10)+8; else W=floor(L/10)+8; end %%W為需要交叉的位數 p=(L-W+1);%隨機產生一個交叉位置 %fprintf('p=%d ',p);%交叉位置 for i=1:W x=find(A==B(1,p+i-1)); y=find(B==A(1,p+i-1)); [A(1,p+i-1),B(1,p+i-1)]=exchange(A(1,p+i-1),B(1,p+i-1)); [A(1,x),B(1,y)]=exchange(A(1,x),B(1,y)); end end
%對調函數 exchange.m function [x,y]=exchange(x,y) temp=x; x=y; y=temp; end
%適應度函數fit.m，每次迭代都要計算每個染色體在本種群內部的優先級別，類似歸一化參數。越大約好！ function fitness=fit(len,m,maxlen,minlen) fitness=len; for i=1:length(len) fitness(i,1)=(1-(len(i,1)-minlen)/(maxlen-minlen+0.0001)).^m; end
%變異函數 Mutation.m function a=Mutation(A) index1=0;index2=0; nnper=randperm(size(A,2)); index1=nnper(1); index2=nnper(2); %fprintf('index1=%d ',index1); %fprintf('index2=%d ',index2); temp=0; temp=A(index1); A(index1)=A(index2); A(index2)=temp; a=A; end
%染色體的路程代價函數 mylength.m function len=myLength(D,p)%p是一個排列 [N,NN]=size(D); len=D(p(1,N),p(1,1)); for i=1:(N-1) len=len+D(p(1,i),p(1,i+1)); end end
%連點畫圖函數 plot_route.m function plot_route(a,R) scatter(a(:,1),a(:,2),'rx'); hold on; plot([a(R(1),1),a(R(length(R)),1)],[a(R(1),2),a(R(length(R)),2)]); hold on; for i=2:length(R) x0=a(R(i-1),1); y0=a(R(i-1),2); x1=a(R(i),1); y1=a(R(i),2); xx=[x0,x1]; yy=[y0,y1]; plot(xx,yy); hold on; end end