模擬退火算法
模擬退火算法來源于固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時，固體內部粒子隨溫升變為無序狀內能增大，而徐徐冷卻時粒子漸趨有序，在每個溫度都達到平衡態，最后在常溫時達到基態，內能減為最小。根據Metropolis準則，粒子在溫度T時趨于平衡的概率為e-ΔE/(kT)，其中E為溫度T時的內能，ΔE為其改變量，k為Boltzmann常數。用固體退火模擬組合優化問題，將內能E模擬為目標函數值f，溫度T演化成控制參數t，即得到解組合優化問題的模擬退火算法：由初始解i和控制參數初值t開始，對當前解重復“產生新解→計算目標函數差→接受或舍棄”的迭代，并逐步衰減t值，算法終止時的當前解即為所得近似最優解，這是基于蒙特卡羅迭代求解法的一種啟發式隨機搜索過程。退火過程由冷卻進度表(Cooling Schedule)控制，包括控制參數的初值t及其衰減因子Δt、每個t值時的迭代次數L和停止條件S。
3.5.1 模擬退火算法的模型
模擬退火算法可以分解為解空間、目標函數和初始解三部分。
模擬退火的基本思想:
(1) 初始化：初始溫度T(充分大)，初始解狀態S(是算法迭代的起點)， 每個T值的迭代次數L
(2) 對k=1，……，L做第(3)至第6步：
(3) 產生新解S′
(4) 計算增量Δt′=C(S′)-C(S)，其中C(S)為評價函數
(5) 若Δt′<0則接受S′作為新的當前解，否則以概率exp(-Δt′/T)接受S′作為新的當前解.
(6) 如果滿足終止條件則輸出當前解作為最優解，結束程序。
終止條件通常取為連續若干個新解都沒有被接受時終止算法。
(7) T逐漸減少，且T->0，然后轉第2步。
算法對應動態演示圖：
模擬退火算法新解的產生和接受可分為如下四個步驟：
第一步是由一個產生函數從當前解產生一個位于解空間的新解；為便于后續的計算和接受，減少算法耗時，通常選擇由當
前新解經過簡單地變換即可產生新解的方法，如對構成新解的全部或部分元素進行置換、互換等，注意到產生新解的變換方法
決定了當前新解的鄰域結構，因而對冷卻進度表的選取有一定的影響。
第二步是計算與新解所對應的目標函數差。因為目標函數差僅由變換部分產生，所以目標函數差的計算最好按增量計算。
事實表明，對大多數應用而言，這是計算目標函數差的最快方法。
第三步是判斷新解是否被接受,判斷的依據是一個接受準則，最常用的接受準則是Metropo1is準則: 若Δt′<0則接受S′作
為新的當前解S，否則以概率exp(-Δt′/T)接受S′作為新的當前解S。
第四步是當新解被確定接受時，用新解代替當前解，這只需將當前解中對應于產生新解時的變換部分予以實現，同時修正
目標函數值即可。此時，當前解實現了一次迭代。可在此基礎上開始下一輪試驗。而當新解被判定為舍棄時，則在原當前解的
基礎上繼續下一輪試驗。
模擬退火算法與初始值無關，算法求得的解與初始解狀態S(是算法迭代的起點)無關；模擬退火算法具有漸近收斂性，已在
理論上被證明是一種以概率l 收斂于全局最優解的全局優化算法；模擬退火算法具有并行性。
3.5.2 模擬退火算法的簡單應用
作為模擬退火算法應用，討論貨郎擔問題(Travelling Salesman Problem，簡記為TSP)：設有n個城市，用數碼1,…,n代表
。城市i和城市j之間的距離為d(i，j) i, j=1,…,n．TSP問題是要找遍訪每個域市恰好一次的一條回路，且其路徑總長度為最
短.。
求解TSP的模擬退火算法模型可描述如下：
解空間 解空間S是遍訪每個城市恰好一次的所有回路，是{1，……，n}的所有循環排列的集合，S中的成員記為(w1,w2 ,…
…，wn)，并記wn+1= w1。初始解可選為(1，……，n)
目標函數 此時的目標函數即為訪問所有城市的路徑總長度或稱為代價函數：
我們要求此代價函數的最小值。
新解的產生 隨機產生1和n之間的兩相異數k和m，若k<m，則將
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
變為：
(w1, w2 ,…，wm , wm-1 ,…，wk+1 , wk ,…，wn).
如果是k>m，則將
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
變為：
(wm, wm-1 ,…，w1 , wm+1 ,…，wk-1 ,wn , wn-1 ,…，wk).
上述變換方法可簡單說成是“逆轉中間或者逆轉兩端”。
也可以采用其他的變換方法，有些變換有獨特的優越性，有時也將它們交替使用，得到一種更好方法。
代價函數差 設將(w1, w2 ,……，wn)變換為(u1, u2 ,……，un), 則代價函數差為：
根據上述分析，可寫出用模擬退火算法求解TSP問題的偽程序：
Procedure TSPSA:
begin
init-of-T; { T為初始溫度}
S={1，……，n}; {S為初始值}
termination=false;
while termination=false
begin
for i=1 to L do
begin
generate(S′form S); { 從當前回路S產生新回路S′}
Δt:=f(S′))-f(S);{f(S)為路徑總長}
IF(Δt<0) OR (EXP(-Δt/T)>Random-of-[0,1])
S=S′;
IF the-halt-condition-is-TRUE THEN
termination=true;
End;
T_lower;
End;
End
模擬退火算法的應用很廣泛，可以較高的效率求解最大截問題(Max Cut Problem)、0-1背包問題(Zero One Knapsack
Problem)、圖著色問題(Graph Colouring Problem)、調度問題(Scheduling Problem)等等。
3.5.3 模擬退火算法的參數控制問題
模擬退火算法的應用很廣泛，可以求解NP完全問題，但其參數難以控制，其主要問題有以下三點：
(1) 溫度T的初始值設置問題。
溫度T的初始值設置是影響模擬退火算法全局搜索性能的重要因素之一、初始溫度高，則搜索到全局最優解的可能性大，但
因此要花費大量的計算時間；反之，則可節約計算時間，但全局搜索性能可能受到影響。實際應用過程中，初始溫度一般需要
依據實驗結果進行若干次調整。
(2) 退火速度問題。
模擬退火算法的全局搜索性能也與退火速度密切相關。一般來說，同一溫度下的“充分”搜索(退火)是相當必要的，但這
需要計算時間。實際應用中，要針對具體問題的性質和特征設置合理的退火平衡條件。
(3) 溫度管理問題。
溫度管理問題也是模擬退火算法難以處理的問題之一。實際應用中，由于必須考慮計算復雜度的切實可行性等問題，常采
用如下所示的降溫方式：
T(t+1)＝k×T(t)
式中k為正的略小于1.00的常數，t為降溫的次數
使用SA解決TSP問題的Matlab程序：
function out = tsp(loc)
% TSP Traveling salesman problem (TSP) using SA (simulated annealing).
% TSP by itself will generate 20 cities within a unit cube and
% then use SA to slove this problem.
%
% TSP(LOC) solve the traveling salesman problem with cities'
% coordinates given by LOC, which is an M by 2 matrix and M is
% the number of cities.
%
% For example:
%
% loc = rand(50, 2);
% tsp(loc);
if nargin == 0,
% The following data is from the post by Jennifer Myers (jmyers@nwu.
edu)
edu)
% to comp.ai.neural-nets. It's obtained from the figure in
% Hopfield & Tank's 1985 paper in Biological Cybernetics
% (Vol 52, pp. 141-152).
loc = [0.3663, 0.9076; 0.7459, 0.8713; 0.4521, 0.8465;
0.7624, 0.7459; 0.7096, 0.7228; 0.0710, 0.7426;
0.4224, 0.7129; 0.5908, 0.6931; 0.3201, 0.6403;
0.5974, 0.6436; 0.3630, 0.5908; 0.6700, 0.5908;
0.6172, 0.5495; 0.6667, 0.5446; 0.1980, 0.4686;
0.3498, 0.4488; 0.2673, 0.4274; 0.9439, 0.4208;
0.8218, 0.3795; 0.3729, 0.2690; 0.6073, 0.2640;
0.4158, 0.2475; 0.5990, 0.2261; 0.3927, 0.1947;
0.5347, 0.1898; 0.3960, 0.1320; 0.6287, 0.0842;
0.5000, 0.0396; 0.9802, 0.0182; 0.6832, 0.8515];
end
NumCity = length(loc); % Number of cities
distance = zeros(NumCity); % Initialize a distance matrix
% Fill the distance matrix
for i = 1:NumCity,
for j = 1:NumCity,
distance(i, j) = norm(loc(i, - loc(j, );
distance(i, j) = norm(loc(i, - loc(j, );
end
end
% To generate energy (objective function) from path
%path = randperm(NumCity);
%energy = sum(distance((path-1)*NumCity + [path(2:NumCity) path(1)]));
% Find typical values of dE
count = 20;
all_dE = zeros(count, 1);
for i = 1:count
path = randperm(NumCity);
energy = sum(distance((path-1)*NumCity + [path(2:NumCity)
path(1)]));
new_path = path;
index = round(rand(2,1)*NumCity+.5);
inversion_index = (min(index):max(index));
new_path(inversion_index) = fliplr(path(inversion_index));
all_dE(i) = abs(energy - ...
sum(sum(diff(loc([new_path new_path(1)],)'.^2)));
end
dE = max(all_dE);
dE = max(all_dE);
temp = 10*dE; % Choose the temperature to be large enough
fprintf('Initial energy = %f"n"n',energy);
% Initial plots
out = [path path(1)];
plot(loc(out(, 1), loc(out(, 2),'r.', 'Markersize', 20);
axis square; hold on
h = plot(loc(out(, 1), loc(out(, 2)); hold off
MaxTrialN = NumCity*100; % Max. # of trials at a
temperature
MaxAcceptN = NumCity*10; % Max. # of acceptances at a
temperature
StopTolerance = 0.005; % Stopping tolerance
TempRatio = 0.5; % Temperature decrease ratio
minE = inf; % Initial value for min. energy
maxE = -1; % Initial value for max. energy
% Major annealing loop
while (maxE - minE)/maxE > StopTolerance,
minE = inf;
minE = inf;
maxE = 0;
TrialN = 0; % Number of trial moves
AcceptN = 0; % Number of actual moves
while TrialN < MaxTrialN & AcceptN < MaxAcceptN,
new_path = path;
index = round(rand(2,1)*NumCity+.5);
inversion_index = (min(index):max(index));
new_path(inversion_index) =
fliplr(path(inversion_index));
new_energy = sum(distance( ...
(new_path-1)*NumCity+[new_path(2:NumCity)
new_path(1)]));
if rand < exp((energy - new_energy)/temp), %
accept
it!
energy = new_energy;
path = new_path;
minE = min(minE, energy);
maxE = max(maxE, energy);
AcceptN = AcceptN + 1;
end
TrialN = TrialN + 1;
end
end
% Update plot
out = [path path(1)];
set(h, 'xdata', loc(out(, 1), 'ydata', loc(out(, 2));
drawnow;
% Print information in command window
fprintf('temp. = %f"n', temp);
tmp = sprintf('%d ',path);
fprintf('path = %s"n', tmp);
fprintf('energy = %f"n', energy);
fprintf('[minE maxE] = [%f %f]"n', minE, maxE);
fprintf('[AcceptN TrialN] = [%d %d]"n"n', AcceptN, TrialN);
% Lower the temperature
temp = temp*TempRatio;
end
% Print sequential numbers in the graphic window
for i = 1:NumCity,
text(loc(path(i),1)+0.01, loc(path(i),2)+0.01, int2str(i), ...
'fontsize', 8);
end
又一個相關的Matlab程序
function[X,P]=zkp(w,c,M,t0,tf)
[m,n]=size(w);
L=100*n;
t=t0;
clear m;
x=zeros(1,n)
xmax=x;
fmax=0;
while t>tf
for k=1:L
xr=change(x)
gx=g_0_1(w,x);
gxr=g_0_1(w,xr);
if gxr<=M
fr=f_0_1(xr,c);
f=f_0_1(x,c);
df=fr-f;
if df>0
x=xr;
if fr>fmax
fmax=fr;
xmax=xr;
end
else
p=rand;
if p<exp(df/t)
x=xr;
end
end
end
end
t=0.87*t
end
P=fmax;
X=xmax;
%下面的函數產生新解
function [d_f,pi_r]=exchange_2(pi0,d)
[m,n]=size(d);
clear m;
u=rand;
u=u*(n-2);
u=round(u);
if u<2
u=2;
end
if u>n-2
u=n-2;
end
v=rand;
v=v*(n-u+1);
v=round(v);
if v<1
v=1;
end
v=v+u;
if v>n
v=n;
end
pi_1(u)=pi0(v);
pi_1(u)=pi0(u);
if u>1
for k=1:(u-1)
pi_1(k)=pi0(k);
end
end
if v>(u+1)
for k=1:(v-u-1)
pi_1(u+k)=pi0(v-k);
end
end
if v<n
for k=(v+1):n
pi_1(k)=pi0(k);
end
end
d_f=0;
if v<n
d_f=d(pi0(u-1),pi0(v))+d(pi0(u),pi0(v+1));
for k=(u+1):n
d_f=d_f+d(pi0(k),pi0(k-1))-d(pi0(v),pi0(v+1));
end
d_f=d_f-d(pi0(u-1),pi0(u));
for k=(u+1):n
d_f=d_f-d(pi0(k-1),pi0(k));
end
else
d_f=d(pi0(u-1),pi0(v))+d(pi0(u),pi0(1))-d(pi0(u-1),pi0(u))-d(pi0(v),pi0(1));
for k=(u+1):n
d_f=d_f-d(pi0(k),pi0(k-1));
end
for k=(u+1):n
d_f=d_f-d(pi0(k-1),pi0(k));
end
end
pi_r=pi_1;
遺傳算法：
旅行商問題(traveling saleman problem,簡稱tsp)：
已知n個城市之間的相互距離，現有一個推銷員必須遍訪這n個城市，并且每個城市只能訪問一次，最后又必須返回出發城市。如何安排他對這些城市的訪問次序，可使其旅行路線的總長度最短？
用圖論的術語來說，假設有一個圖 g=(v,e)，其中v是頂點集，e是邊集，設d=(dij)是由頂點i和頂點j之間的距離所組成的距離矩陣，旅行商問題就是求出一條通過所有頂點且每個頂點只通過一次的具有最短距離的回路。
這個問題可分為對稱旅行商問題(dij=dji,,任意i,j=1,2,3，…,n)和非對稱旅行商問題(dij≠dji,,任意i,j=1,2,3，…,n)。
若對于城市v={v1,v2,v3,…,vn}的一個訪問順序為t=(t1,t2,t3,…,ti,…,tn),其中ti∈v(i=1,2,3,…,n)，且記tn+1= t1，則旅行商問題的數學模型為：
min l=σd(t(i),t(i+1)) （i=1,…,n）
旅行商問題是一個典型的組合優化問題，并且是一個np難問題，其可能的路徑數目與城市數目n是成指數型增長的，所以一般很難精確地求出其最優解，本文采用遺傳算法求其近似解。
遺傳算法：
初始化過程：用v1,v2,v3,…,vn代表所選n個城市。定義整數pop-size作為染色體的個數，并且隨機產生pop-size個初始染色體，每個染色體為1到18的整數組成的隨機序列。
適應度f的計算：對種群中的每個染色體vi，計算其適應度，f=σd(t(i),t(i+1)).
評
價函數eval(vi)：用來對種群中的每個染色體vi設定一個概率，以使該染色體被選中的可能性與其種群中其它染色體的適應性成比例，既通過輪盤賭，適
應性強的染色體被選擇產生后臺的機會要大，設alpha∈(0,1)，本文定義基于序的評價函數為eval(vi)=alpha*(1-alpha).^
(i-1) 。[隨機規劃與模糊規劃]
選擇過程：選擇過程是以旋轉賭輪pop-size次為基礎，每次旋轉都為新的種群選擇一個染色體。賭輪是按每個染色體的適應度進行選擇染色體的。
step1 、對每個染色體vi,計算累計概率qi，q0=0;qi=σeval(vj) j=1,…,i;i=1,…pop-size.
step2、從區間(0,pop-size)中產生一個隨機數r；
step3、若qi-1<r<qi,則選擇第i個染色體 ；
step4、重復step2和step3共pop-size次，這樣可以得到pop-size個復制的染色體。
grefenstette
編碼：由于常規的交叉運算和變異運算會使種群中產生一些無實際意義的染色體，本文采用grefenstette編碼《遺傳算法原理及應用》可以避免這種情
況的出現。所謂的grefenstette編碼就是用所選隊員在未選（不含淘汰）隊員中的位置，如：
8 15 2 16 10 7 4 3 11 14 6 12 9 5 18 13 17 1
對應：
8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1。
交叉過程：本文采用常規單點交叉。為確定交叉操作的父代，從 到pop-size重復以下過程：從[0，1]中產生一個隨機數r，如果r<pc ，則選擇vi作為一個父代。
將所選的父代兩兩組隊，隨機產生一個位置進行交叉，如：
8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1
6 12 3 5 6 8 5 6 3 1 8 5 6 3 3 2 1 1
交叉后為：
8 14 2 13 8 6 3 2 5 1 8 5 6 3 3 2 1 1
6 12 3 5 6 8 5 6 3 7 3 4 3 2 4 2 2 1
變
異過程：本文采用均勻多點變異。類似交叉操作中選擇父代的過程，在r<pm
的標準下選擇多個染色體vi作為父代。對每一個選擇的父代，隨機選擇多個位置，使其在每位置按均勻變異（該變異點xk的取值范圍為[ukmin,
ukmax],產生一個[0，1]中隨機數r，該點變異為x'k=ukmin+r(ukmax-ukmin)）操作。如：
8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1
變異后：
8 14 2 13 10 6 3 2 2 7 3 4 5 2 4 1 2 1
反grefenstette編碼：交叉和變異都是在grefenstette編碼之后進行的，為了循環操作和返回最終結果，必須逆grefenstette編碼過程，將編碼恢復到自然編碼。
循環操作：判斷是否滿足設定的帶數xzome，否，則跳入適應度f的計算；是，結束遺傳操作，跳出。
Matlab程序：
function [bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%
%————————————————————————
%[bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%d:距離矩陣
%termops:種群代數
%num:每代染色體的個數
%pc:交叉概率
%cxops:由于本程序采用單點交叉，交叉點的設置在本程序中沒有很好的解決，所以本文了采用定點，即第cxops，可以隨機產生。
%pm:變異概率
%alpha:評價函數eval(vi)=alpha*(1-alpha).^(i-1).
%bestpop:返回的最優種群
%trace:進化軌跡
%------------------------------------------------
%####@@@##版權所有！歡迎廣大網友改正，改進！##@@@####
%e-mail:tobysidney33@sohu.com
%####################################################
%
citynum=size(d,2);
n=nargin;
if n<2
disp('缺少變量！！')
disp('^_^開個玩笑^_^')
end
if n<2
termops=500;
num=50;
pc=0.25;
cxops=3;
pm=0.30;
alpha=0.10;
end
if n<3
num=50;
pc=0.25;
cxops=3;
pm=0.30;
alpha=0.10;
end
if n<4
pc=0.25;
cxops=3;
pm=0.30;
alpha=0.10;
end
if n<5
cxops=3;
pm=0.30;
alpha=0.10;
end
if n<6
pm=0.30;
alpha=0.10;
end
if n<7
alpha=0.10;
end
if isempty(cxops)
cxops=3;
end
[t]=initializega(num,citynum);
for i=1:termops
[l]=f(d,t);
[x,y]=find(l==max(l));
trace(i)=-l(y(1));
bestpop=t(y(1),:);
[t]=select(t,l,alpha);
[g]=grefenstette(t);
[g1]=crossover(g,pc,cxops);
[g]=mutation(g1,pm); %均勻變異
[t]=congrefenstette(g);
end
---------------------------------------------------------
function [t]=initializega(num,citynum)
for i=1:num
t(i,:)=randperm(citynum);
end
-----------------------------------------------------------
function [l]=f(d,t)
[m,n]=size(t);
for k=1:m
for i=1:n-1
l(k,i)=d(t(k,i),t(k,i+1));
end
l(k,n)=d(t(k,n),t(k,1));
l(k)=-sum(l(k,:));
end
-----------------------------------------------------------
function [t]=select(t,l,alpha)
[m,n]=size(l);
t1=t;
[beforesort,aftersort1]=sort(l,2);%fsort from l to u
for i=1:n
aftersort(i)=aftersort1(n+1-i); %change
end
for k=1:n;
t(k,:)=t1(aftersort(k),:);
l1(k)=l(aftersort(k));
end
t1=t;
l=l1;
for i=1:size(aftersort,2)
evalv(i)=alpha*(1-alpha).^(i-1);
end
m=size(t,1);
q=cumsum(evalv);
qmax=max(q);
for k=1:m
r=qmax*rand(1);
for j=1:m
if j==1&r<=q(1)
t(k,:)=t1(1,:);
elseif j~=1&r>q(j-1)&r<=q(j)
t(k,:)=t1(j,:);
end
end
end
--------------------------------------------------
function [g]=grefenstette(t)
[m,n]=size(t);
for k=1:m
t0=1:n;
for i=1:n
for j=1:length(t0)
if t(k,i)==t0(j)
g(k,i)=j;
t0(j)=[];
break
end
end
end
end
-------------------------------------------
function [g]=crossover(g,pc,cxops)
[m,n]=size(g);
ran=rand(1,m);
r=cxops;
[x,ru]=find(ran<pc);
if ru>=2
for k=1:2:length(ru)-1
g1(ru(k),:)=[g(ru(k),[1:r]),g(ru(k+1),[(r+1):n])];
g(ru(k+1),:)=[g(ru(k+1),[1:r]),g(ru(k),[(r+1):n])];
g(ru(k),:)=g1(ru(k),:);
end
end
--------------------------------------------
function [g]=mutation(g,pm) %均勻變異
[m,n]=size(g);
ran=rand(1,m);
r=rand(1,3); %dai gai jin
rr=floor(n*rand(1,3)+1);
[x,mu]=find(ran<pm);
for k=1:length(mu)
for i=1:length(r)
umax(i)=n+1-rr(i);
umin(i)=1;
g(mu(k),rr(i))=umin(i)+floor((umax(i)-umin(i))*r(i));
end
end
---------------------------------------------------
function [t]=congrefenstette(g)
[m,n]=size(g);
for k=1:m
t0=1:n;
for i=1:n
t(k,i)=t0(g(k,i));
t0(g(k,i))=[];
end
end
-------------------------------------------------
又一個Matlab程序，其中交叉算法采用的是由Goldberg和Lingle于1985年提出的PMX(部分匹配交叉)，淘汰保護指數alpha是我自己設計的，起到了加速優勝劣汰的作用。
%TSP問題（又名：旅行商問題，貨郎擔問題）遺傳算法通用matlab程序
%D是距離矩陣，n為種群個數,建議取為城市個數的1~2倍，
%C為停止代數，遺傳到第 C代時程序停止,C的具體取值視問題的規模和耗費的時間而定
%m為適應值歸一化淘汰加速指數 ,最好取為1,2,3,4 ,不宜太大
%alpha為淘汰保護指數,可取為0~1之間任意小數,取1時關閉保護功能,最好取為0.8~1.0
%R為最短路徑,Rlength為路徑長度
function [R,Rlength]=geneticTSP(D,n,C,m,alpha)
[N,NN]=size(D);
farm=zeros(n,N);%用于存儲種群
for i=1:n
farm(i,:)=randperm(N);%隨機生成初始種群
end
R=farm(1,:);%存儲最優種群
len=zeros(n,1);%存儲路徑長度
fitness=zeros(n,1);%存儲歸一化適應值
counter=0;
while counter<C
for i=1:n
len(i,1)=myLength(D,farm(i,:));%計算路徑長度
end
maxlen=max(len);
minlen=min(len);
fitness=fit(len,m,maxlen,minlen);%計算歸一化適應值
rr=find(len==minlen);
R=farm(rr(1,1),:);%更新最短路徑
FARM=farm;%優勝劣汰，nn記錄了復制的個數
nn=0;
for i=1:n
if fitness(i,1)>=alpha*rand
nn=nn+1;
FARM(nn,:)=farm(i,:);
end
end
FARM=FARM(1:nn,:);
[aa,bb]=size(FARM);%交叉和變異
while aa<n
if nn<=2
nnper=randperm(2);
else
nnper=randperm(nn);
end
A=FARM(nnper(1),:);
B=FARM(nnper(2),:);
[A,B]=intercross(A,B);
FARM=[FARM;A;B];
[aa,bb]=size(FARM);
end
if aa>n
FARM=FARM(1:n,:);%保持種群規模為n
end
farm=FARM;
clear FARM
counter=counter+1
end
Rlength=myLength(D,R);
function [a,b]=intercross(a,b)
L=length(a);
if L<=10%確定交叉寬度
W=1;
elseif ((L/10)-floor(L/10))>=rand&&L>10
W=ceil(L/10);
else
W=floor(L/10);
end
p=unidrnd(L-W+1);%隨機選擇交叉范圍，從p到p+W
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
function [x,y]=exchange(x,y)
temp=x;
x=y;
y=temp;
% 計算路徑的子程序
function len=myLength(D,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
%計算歸一化適應值子程序
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.000001))).^m;
end
一個C++的程序：
//c++的程序
#include<iostream.h>
#include<stdlib.h>
template<class T>
class Graph
{
public:
Graph(int vertices=10)
{
n=vertices;
e=0;
}
~Graph(){}
virtual bool Add(int u,int v,const T& w)=0;
virtual bool Delete(int u,int v)=0;
virtual bool Exist(int u,int v)const=0;
int Vertices()const{return n;}
int Edges()const{return e;}
protected:
int n;
int e;
};
template<class T>
class MGraph:public Graph<T>
{
public:
MGraph(int Vertices=10,T noEdge=0);
~MGraph();
bool Add(int u,int v,const T& w);
bool Delete(int u,int v);
bool Exist(int u,int v)const;
void Floyd(T**& d,int**& path);
void print(int Vertices);
private:
T NoEdge;
T** a;
};
template<class T>
MGraph<T>::MGraph(int Vertices,T noEdge)
{
n=Vertices;
NoEdge=noEdge;
a=new T* [n];
for(int i=0;i<n;i++){
a[i]=new T[n];
a[i][i]=0;
for(int j=0;j<n;j++)if(i!=j)a[i][j]=NoEdge;
}
}
template<class T>
MGraph<T>::~MGraph()
{
for(int i=0;i<n;i++)delete[]a[i];
delete[]a;
}
template<class T>
bool MGraph<T>::Exist(int u,int v)const
{
if(u<0||v<0||u>n-1||v>n-1||u==v||a[u][v]==NoEdge)return false;
return true;
}
template<class T>
bool MGraph<T>::Add(int u,int v,const T& w)
{
if(u<0||v<0||u>n-1||v>n-1||u==v||a[u][v]!=NoEdge){
cerr<<"BadInput!"<<endl;
return false;
}
a[u][v]=w;
e++;
return true;
}
template<class T>
bool MGraph<T>:delete(int u,int v)
{
if(u<0||v<0||u>n-1||v>n-1||u==v||a[u][v]==NoEdge){
cerr<<"BadInput!"<<endl;
return false;
}
a[u][v]=NoEdge;
e--;
return true;
}
template<class T>
void MGraph<T>::Floyd(T**& d,int**& path)
{
d=new T* [n];
path=new int* [n];
for(int i=0;i<n;i++){
d[i]=new T[n];
path[i]=new int[n];
for(int j=0;j<n;j++){
d[i][j]=a[i][j];
if(i!=j&&a[i][j]<NoEdge)path[i][j]=i;
else path[i][j]=-1;
}
}
for(int k=0;k<n;k++){
for(i=0;i<n;i++)
for(int j=0;j<n;j++)
if(d[i][k]+d[k][j]<d[i][j]){
d[i][j]=d[i][k]+d[k][j];
path[i][j]=path[k][j];
}
}
}
template<class T>
void MGraph<T>::print(int Vertices)
{
for(int i=0;i<Vertices;i++)
for(int j=0;j<Vertices;j++)
{
cout<<a[i][j]<<' ';if(j==Vertices-1)cout<<endl;
}
}
#define noEdge 10000
#include<iostream.h>
void main()
{
cout<<"請輸入該圖的節點數："<<endl;
int vertices;
cin>>vertices;
MGraph<float> b(vertices,noEdge);
cout<<"請輸入u,v,w:"<<endl;
int u,v;
float w;
cin>>u>>v>>w;
while(w!=noEdge){
//u=u-1;
b.Add(u-1,v-1,w);
b.Add(v-1,u-1,w);
cout<<"請輸入u,v,w:"<<endl;
cin>>u>>v>>w;
}
b.print(vertices);
int** Path;
int**& path=Path;
float** D;
float**& d=D;
b.Floyd(d,path);
for(int i=0;i<vertices;i++){
for(int j=0;j<vertices;j++){
cout<<Path[i][j]<<' ';
if(j==vertices-1)cout<<endl;
}
}
int *V;
V=new int[vertices+1];
cout<<"請輸入任意一個初始H-圈："<<endl;
for(int n=0;n<=vertices;n++){
cin>>V[n];
}
for(n=0;n<55;n++){
for(i=0;i<n-1;i++){
for(int j=0;j<n-1;j++)
{
if(i+1>0&&j>i+1&&j<n-1){
if(D[V[i]][V[j]]+D[V[i+1]][V[j+1]]<D[V[i]][V[i+1]]+D[V[j]][V[j+1]]){
int l;
l=V[i+1];V[i+1]=V[j];V[j]=l;
}
}
}
}
}
float total=0;
cout<<"最小回路："<<endl;
for(i=0;i<=vertices;i++){
cout<<V[i]+1<<' ';
}
cout<<endl;
for(i=0;i<vertices;i++)
total+=D[V[i]][V[i+1]];
cout<<"最短路徑長度："<<endl;
cout<<total;
}
C語言程序：
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<alloc.h>
#include<conio.h>
#include<float.h>
#include<time.h>
#include<graphics.h>
#include<bios.h>
#define maxpop 100
#define maxstring 100
struct pp{unsigned char chrom[maxstring];
float x,fitness;
unsigned int parent1,parent2,xsite;
};
struct pp *oldpop,*newpop,*p1;
unsigned int popsize,lchrom,gem,maxgen,co_min,jrand;
unsigned int nmutation,ncross,jcross,maxpp,minpp,maxxy;
float pcross,pmutation,sumfitness,avg,max,min,seed,maxold,oldrand[maxstring];
unsigned char x[maxstring],y[maxstring];
float *dd,ff,maxdd,refpd,fm[201];
FILE *fp,*fp1;
float objfunc(float);
void statistics();
int select();
int flip(float);
int crossover();
void generation();
void initialize();
void report();
float decode();
void crtinit();
void inversion();
float random1();
void randomize1();
main()
{unsigned int gen,k,j,tt;
char fname[10];
float ttt;
clrscr();
co_min=0;
if((oldpop=(struct pp *)farmalloc(maxpop*sizeof(struct pp)))==NULL)
{printf("memory requst fail!"n");exit(0);}
if((dd=(float *)farmalloc(maxstring*maxstring*sizeof(float)))==NULL)
{printf("memory requst fail!"n");exit(0);}
if((newpop=(struct pp *)farmalloc(maxpop*sizeof(struct pp)))==NULL)
{printf("memory requst fail!"n");exit(0);}
if((p1=(struct pp *)farmalloc(sizeof(struct pp)))==NULL)
{printf("memory requst fail!"n");exit(0);}
for(k=0;k<maxpop;k++) oldpop[k].chrom[0]='"0';
for(k=0;k<maxpop;k++) newpop[k].chrom[0]='"0';
printf("Enter Result Data Filename:");
gets(fname);
if((fp=fopen(fname,"w+"))==NULL)
{printf("cannot open file"n");exit(0);}
gen=0;
randomize();
initialize();
fputs("this is result of the TSP problem:",fp);
fprintf(fp,"city: %2d psize: %3d Ref.TSP_path: %f"n",lchrom,popsize,refpd);
fprintf(fp,"Pc: %f Pm: %f Seed: %f"n",pcross,pmutation,seed);
fprintf(fp,"X site:"n");
for(k=0;k<lchrom;k++)
{if((k%16)==0) fprintf(fp,""n");
fprintf(fp,"%5d",x[k]);
}
fprintf(fp,""n Y site:"n");
for(k=0;k<lchrom;k++)
{if((k%16)==0) fprintf(fp,""n");
fprintf(fp,"%5d",y[k]);
}
fprintf(fp,""n");
crtinit();
statistics(oldpop);
report(gen,oldpop);
getch();
maxold=min;
fm[0]=100.0*oldpop[maxpp].x/ff;
do {
gen=gen+1;
generation();
statistics(oldpop);
if(max>maxold)
{maxold=max;
co_min=0;
}
fm[gen%200]=100.0*oldpop[maxpp].x/ff;
report(gen,oldpop);
gotoxy(30,25);
ttt=clock()/18.2;
tt=ttt/60;
printf("Run Clock: %2d: %2d: %4.2f",tt/60,tt%60,ttt-tt*60.0);
printf("Min=%6.4f Nm:%d"n",min,co_min);
}while((gen<100)&&!bioskey(1));
printf(""n gen= %d",gen);
do{
gen=gen+1;
generation();
statistics(oldpop);
if(max>maxold)
{maxold=max;
co_min=0;
}
fm[gen%200]=100.0*oldpop[maxpp].x/ff;
report(gen,oldpop);
if((gen%100)==0)report(gen,oldpop);
gotoxy(30,25);
ttt=clock()/18.2;
tt=ttt/60;
printf("Run Clock: %2d: %2d: %4.2f",tt/60,tt%60,ttt-tt*60.0);
printf("Min=%6.4f Nm:%d"n",min,co_min);
}while((gen<maxgen)&&!bioskey(1));
getch();
for(k=0;k<lchrom;k++)
{if((k%16)==0)fprintf(fp,""n");
fprintf(fp,"%5d",oldpop[maxpp].chrom[k]);
}
fprintf(fp,""n");
fclose(fp);
farfree(dd);
farfree(p1);
farfree(oldpop);
farfree(newpop);
restorecrtmode();
exit(0);
}
/*%%%%%%%%%%%%%%%%*/
float objfunc(float x1)
{float y;
y=100.0*ff/x1;
return y;
}
/*&&&&&&&&&&&&&&&&&&&*/
void statistics(pop)
struct pp *pop;
{int j;
sumfitness=pop[0].fitness;
min=pop[0].fitness;
max=pop[0].fitness;
maxpp=0;
minpp=0;
for(j=1;j<popsize;j++)
{sumfitness=sumfitness+pop[j].fitness;
if(pop[j].fitness>max)
{max=pop[j].fitness;
maxpp=j;
}
if(pop[j].fitness<min)
{min=pop[j].fitness;
minpp=j;
}
}
avg=sumfitness/(float)popsize;
}
/*%%%%%%%%%%%%%%%%%%%%*/
void generation()
{unsigned int k,j,j1,j2,i1,i2,mate1,mate2;
float f1,f2;
j=0;
do{
mate1=select();
pp:mate2=select();
if(mate1==mate2)goto pp;
crossover(oldpop[mate1].chrom,oldpop[mate2].chrom,j);
newpop[j].x=(float)decode(newpop[j].chrom);
newpop[j].fitness=objfunc(newpop[j].x);
newpop[j].parent1=mate1;
newpop[j].parent2=mate2;
newpop[j].xsite=jcross;
newpop[j+1].x=(float)decode(newpop[j+1].chrom);
newpop[j+1].fitness=objfunc(newpop[j+1].x);
newpop[j+1].parent1=mate1;
newpop[j+1].parent2=mate2;
newpop[j+1].xsite=jcross;
if(newpop[j].fitness>min)
{for(k=0;k<lchrom;k++)
oldpop[minpp].chrom[k]=newpop[j].chrom[k];
oldpop[minpp].x=newpop[j].x;
oldpop[minpp].fitness=newpop[j].fitness;
co_min++;
return;
}
if(newpop[j+1].fitness>min)
{for(k=0;k<lchrom;k++)
oldpop[minpp].chrom[k]=newpop[j+1].chrom[k];
oldpop[minpp].x=newpop[j+1].x;
oldpop[minpp].fitness=newpop[j+1].fitness;
co_min++;
return;
}
j=j+2;
}while(j<popsize);
}
/*%%%%%%%%%%%%%%%%%*/
void initdata()
{unsigned int ch,j;
clrscr();
printf("-----------------------"n");
printf("A SGA"n");
printf("------------------------"n");
/*pause();*/clrscr();
printf("*******SGA DATA ENTRY AND INITILIZATION *******"n");
printf(""n");
printf("input pop size");scanf("%d",&popsize);
printf("input chrom length");scanf("%d",&lchrom);
printf("input max generations");scanf("%d",&maxgen);
printf("input crossover probability");scanf("%f",&pcross);
printf("input mutation prob");scanf("%f",&pmutation);
randomize1();
clrscr();
nmutation=0;
ncross=0;
}
/*%%%%%%%%%%%%%%%%%%%%*/
void initreport()
{int j,k;
printf("pop size=%d"n",popsize);
printf("chromosome length=%d"n",lchrom);
printf("maxgen=%d"n",maxgen);
printf("pmutation=%f"n",pmutation);
printf("pcross=%f"n",pcross);
printf("initial generation statistics"n");
printf("ini pop max fitness=%f"n",max);
printf("ini pop avr fitness=%f"n",avg);
printf("ini pop min fitness=%f"n",min);
printf("ini pop sum fit=%f"n",sumfitness);
}
void initpop()
{unsigned char j1;
unsigned int k5,i1,i2,j,i,k,j2,j3,j4,p5[maxstring];
float f1,f2;
j=0;
for(k=0;k<lchrom;k++)
oldpop[j].chrom[k]=k;
for(k=0;k<lchrom;k++)
p5[k]=oldpop[j].chrom[k];
randomize();
for(;j<popsize;j++)
{j2=random(lchrom);
for(k=0;k<j2+20;k++)
{j3=random(lchrom);
j4=random(lchrom);
j1=p5[j3];
p5[j3]=p5[j4];
p5[j4]=j1;
}
for(k=0;k<lchrom;k++)
oldpop[j].chrom[k]=p5[k];
}
for(k=0;k<lchrom;k++)
for(j=0;j<lchrom;j++)
dd[k*lchrom+j]=hypot(x[k]-x[j],y[k]-y[j]);
for(j=0;j<popsize;j++)
{oldpop[j].x=(float)decode(oldpop[j].chrom);
oldpop[j].fitness=objfunc(oldpop[j].x);
oldpop[j].parent1=0;
oldpop[j].parent2=0;
oldpop[j].xsite=0;
}
}
/*&&&&&&&&&&&&&&&&&*/
void initialize()
{int k,j,minx,miny,maxx,maxy;
initdata();
minx=0;
miny=0;
maxx=0;maxy=0;
for(k=0;k<lchrom;k++)
{x[k]=rand();
if(x[k]>maxx)maxx=x[k];
if(x[k]<minx)minx=x[k];
y[k]=rand();
if(y[k]>maxy)maxy=y[k];
if(y[k]<miny)miny=y[k];
}
if((maxx-minx)>(maxy-miny))
{maxxy=maxx-minx;}
else {maxxy=maxy-miny;}
maxdd=0.0;
for(k=0;k<lchrom;k++)
for(j=0;j<lchrom;j++)
{dd[k*lchrom+j]=hypot(x[k]-x[j],y[k]-y[j]);
if(maxdd<dd[k*lchrom+j])maxdd=dd[k*lchrom+j];
}
refpd=dd[lchrom-1];
for(k=0;k<lchrom;k++)
refpd=refpd+dd[k*lchrom+k+2];
for(j=0;j<lchrom;j++)
dd[j*lchrom+j]=4.0*maxdd;
ff=(0.765*maxxy*pow(lchrom,0.5));
minpp=0;
min=dd[lchrom-1];
for(j=0;j<lchrom-1;j++)
{if(dd[lchrom*j+lchrom-1]<min)
{min=dd[lchrom*j+lchrom-1];
minpp=j;
}
}
initpop();
statistics(oldpop);
initreport();
}
/*&&&&&&&&&&&&&&&&&&*/
void report(int l,struct pp *pop)
{int k,ix,iy,jx,jy;
unsigned int tt;
float ttt;
cleardevice();
gotoxy(1,1);
printf("city:%4d para_size:%4d maxgen:%4d ref_tour:%f"n"
,lchrom,popsize,maxgen,refpd);
printf("ncross:%4d Nmutation:%4d Rungen:%4d AVG=%8.4f MIN=%8.4f"n"n"
,ncross,nmutation,l,avg,min);
printf("Ref.cominpath:%6.4f Minpath length:%10.4f Ref_co_tour:%f"n"
,pop[maxpp].x/maxxy,pop[maxpp].x,ff);
printf("Co_minpath:%6.4f Maxfit:%10.8f"
,100.0*pop[maxpp].x/ff,pop[maxpp].fitness);
ttt=clock()/18.2;
tt=ttt/60;
printf("Run clock:%2d:%2d:%4d.2f"n",tt/60,tt%60,ttt-tt*60.0);
setcolor(1%15+1);
for(k=0;k<lchrom-1;k++)
{ix=x[pop[maxpp].chrom[k]];
iy=y[pop[maxpp].chrom[k]]+110;
jx=x[pop[maxpp].chrom[k+1]];
jy=y[pop[maxpp].chrom[k+1]]+110;
line(ix,iy,jx,jy);
putpixel(ix,iy,RED);
}
ix=x[pop[maxpp].chrom[0]];
iy=y[pop[maxpp].chrom[0]]+110;
jx=x[pop[maxpp].chrom[lchrom-1]];
jy=y[pop[maxpp].chrom[lchrom-1]]+110;
line(ix,iy,jx,jy);
putpixel(jx,jy,RED);
setcolor(11);
outtextxy(ix,iy,"*");
setcolor(12);
for(k=0;k<1%200;k++)
{ix=k+280;
iy=366-fm[k]/3;
jx=ix+1;
jy=366-fm[k+1]/3;
line(ix,iy,jx,jy);
putpixel(ix,iy,RED);
}
printf("GEN:%3d",l);
printf("Minpath:%f Maxfit:%f",pop[maxpp].x,pop[maxpp].fitness);
printf("Clock:%2d:%2d:%4.2f"n",tt/60,tt%60,ttt-tt*60.0);
}
/*###############*/
float decode(unsigned char *pp)
{int j,k,l;
float tt;
tt=dd[pp[0]*lchrom+pp[lchrom-1]];
for(j=0;j<lchrom-1;j++)
{tt=tt+dd[pp[j]*lchrom+pp[j+1]];}
l=0;
for(k=0;k<lchrom-1;k++)
for(j=k+1;j<lchrom;j++)
{if(pp[j]==pp[k])l++;}
return tt+4*l*maxdd;
}
/*%%%%%%%%%%%%%%%%%%*/
void crtinit()
{int driver,mode;
struct palettetype p;
driver=DETECT;
mode=0;
initgraph(&driver,&mode,"");
cleardevice();
}
/*$$$$$$$$$$$$$$$$$$$$*/
int select()
{double rand1,partsum;
float r1;
int j;
partsum=0.0;
j=0;
rand1=random1()*sumfitness;
do{
partsum=partsum+oldpop[j].fitness;
j=j+1;
}while((partsum<rand1)&&(j<popsize));
return j-1;
}
/*$$$$$$$$$$$$$$$*/
int crossover(unsigned char *parent1,unsigned char *parent2,int k5)
{int k,j,mutate,i1,i2,j5;
int j1,j2,j3,s0,s1,s2;
unsigned char jj,ts1[maxstring],ts2[maxstring];
float f1,f2;
s0=0;s1=0;s2=0;
if(flip(pcross))
{jcross=random(lchrom-1);
j5=random(lchrom-1);
ncross=ncross+1;
if(jcross>j5){k=jcross;jcross=j5;j5=k;}
}
else jcross=lchrom;
if(jcross!=lchrom)
{s0=1;
k=0;
for(j=jcross;j<j5;j++)
{ts1[k]=parent1[j];
ts2[k]=parent2[j];
k++;
}
j3=k;
for(j=0;j<lchrom;j++)
{j2=0;
while((parent2[j]!=ts1[j2])&&(j2<k)){j2++;}
if(j2==k)
{ts1[j3]=parent2[j];
j3++;
}
}
j3=k;
for(j=0;j<lchrom;j++)
{j2=0;
while((parent1[j]!=ts2[j2])&&(j2<k)){j2++;}
if(j2==k)
{ts2[j3]=parent1[j];
j3++;
}
}
for(j=0;j<lchrom;j++)
{newpop[k5].chrom[j]=ts1[j];
newpop[k5+1].chrom[j]=ts2[j];
}
}
else
{for(j=0;j<lchrom;j++)
{newpop[k5].chrom[j]=parent1[j];
newpop[k5+1].chrom[j]=parent2[j];
}
mutate=flip(pmutation);
if(mutate)
{s1=1;
nmutation=nmutation+1;
for(j3=0;j3<200;j3++)
{j1=random(lchrom);
j=random(lchrom);
jj=newpop[k5].chrom[j];
newpop[k5].chrom[j]=newpop[k5].chrom[j1];
newpop[k5].chrom[j1]=jj;
}
}
mutate=flip(pmutation);
if(mutate)
{s2=1;
nmutation=nmutation+1;
for(j3=0;j3<100;j3++)
{j1=random(lchrom);
j=random(lchrom);
jj=newpop[k5+1].chrom[j];
newpop[k5+1].chrom[j]=newpop[k5+1].chrom[j1];
newpop[k5+1].chrom[j1]=jj;
}
}
}
j2=random(2*lchrom/3);
for(j=j2;j<j2+lchrom/3-1;j++)
for(k=0;k<lchrom;k++)
{if(k==j)continue;
if(k>j){i2=k;i1=j;}
else{i1=k;i2=j;}
f1=dd[lchrom*newpop[k5].chrom[i1]+newpop[k5].chrom[i2]];
f1=f1+dd[lchrom*newpop[k5].chrom[(i1+1)%lchrom]+
newpop[k5].chrom[(i2+1)%lchrom]];
f2=dd[lchrom*newpop[k5].chrom[i1]+
newpop[k5].chrom[(i1+1)%lchrom]];
f2=f2+dd[lchrom*newpop[k5].chrom[i2]+
newpop[k5].chrom[(i2+1)%lchrom]];
if(f1<f2){inversion(i1,i2,newpop[k5].chrom);}
}
j2=random(2*lchrom/3);
for(j=j2;j<j2+lchrom/3-1;j++)
for(k=0;k<lchrom;k++)
{if(k==j)continue;
if(k>j){i2=k;i1=j;}
else{i1=k;i2=j;}
f1=dd[lchrom*newpop[k5+1].chrom[i1]+newpop[k5+1].chrom[i2]];
f1=f1+dd[lchrom*newpop[k5+1].chrom[(i1+1)%lchrom]+
newpop[k5+1].chrom[(i2+1)%lchrom]];
f2=dd[lchrom*newpop[k5+1].chrom[i1]+
newpop[k5+1].chrom[(i1+1)%lchrom]];
f2=f2+dd[lchrom*newpop[k5+1].chrom[i2]+
newpop[k5+1].chrom[(i2+1)%lchrom]];
if(f1<f2){inversion(i1,i2,newpop[k5+1].chrom);}
}
return 1;
}
/*$$$$$$$$$$$$$$$*/
void inversion(unsigned int k,unsigned int j,unsigned char *ss)
{unsigned int l1,i;
unsigned char tt;
l1=(j-k)/2;
for(i=0;i<l1;i++)
{tt=ss[k+i+1];
ss[k+i+1]=ss[j-i];
ss[j-i]=tt;
}
}
/*%%%%%%%%%%%%%%%*/
void randomize1()
{int i;
randomize();
for(i=0;i<lchrom;i++)
oldrand[i]=random(30001)/30000.0;
jrand=0;
}
/*%%%%%%%%%%%*/
float random1()
{jrand=jrand+1;
if(jrand>=lchrom)
{jrand=0;
randomize1();
}
return oldrand[jrand];
}
/*%%%%%%%%%%*/
int flip(float probability)
{float ppp;
ppp=random(20001)/20000.0;
if(ppp<=probability)return 1;
return 0;
}
改進后用來求解VRP問題的Delphi程序：
unit uEA;
interface
uses
uUtilsEA, uIEA, uITSP, Classes, GaPara, windows, SysUtils, fEA_TSP;
type
TIndividual = class(TInterfacedObject, IIndividual)
private
// The internally stored fitness value
fFitness: TFloat;
fWeConstrain: integer;
fBackConstrain: integer;
fTimeConstrain: integer;
procedure SetFitness(const Value: TFloat);
function GetFitness: TFloat;
function GetWeConstrain: integer;
procedure SetWeConstrain(const Value: integer);
procedure SetBackConstrain(const Value: integer);
function GetBackConstrain: integer;
function GetTimeConstrain: integer;
procedure SetTimeConstrain(const Value: integer);
public
property Fitness : TFloat read GetFitness write SetFitness;
property WeConstrain :integer read GetWeConstrain write SetWeConstrain;
property BackConstrain :integer read GetBackConstrain write SetBackConstrain;
property TimeConstrain :integer read GetTimeConstrain write SetTimeConstrain;
end;
TTSPIndividual = class(TIndividual, ITSPIndividual)
private
// The route we travel
fRouteArray : ArrayInt;
fWeConstrain: integer;
fBackConstrain: integer;
fTimeConstrain: integer;
function GetRouteArray(I: Integer): Integer;
procedure SetRouteArray(I: Integer; const Value: Integer);
procedure SetSteps(const Value: Integer);
function GetSteps: Integer;
function GetWeConstrain: integer;
procedure SetWeConstrain(const Value: integer);
procedure SetBackConstrain(const Value: integer);
procedure SetTimeConstrain(const Value: integer);
function GetBackConstrain: integer;
function GetTimeConstrain: integer;
public
// Constructor, called with initial route size
constructor Create(Size : TInt); reintroduce;
destructor Destroy; override;
property RouteArray[I : Integer] : Integer read GetRouteArray write SetRouteArray;
// The number of steps on the route
property Steps : Integer read GetSteps write SetSteps;
property Fitness : TFloat read GetFitness write SetFitness;
property WeConstrain :integer read GetWeConstrain write SetWeConstrain;
property BackConstrain :integer read GetWeConstrain write SetBackConstrain;
property TimeConstrain :integer read GetTimeConstrain write SetTimeConstrain;
end;
TTSPCreator = class(TInterfacedObject, ITSPCreator)
private
// The Control component we are associated with
fController: ITSPController;
function GetController: ITSPController;
procedure SetController(const Value: ITSPController);
public
// Function to create a random individual
function CreateIndividual : IIndividual;
function CreateFeasibleIndividual: IIndividual;
property Controller : ITSPController read GetController write SetController;
end;
TKillerPercentage = class(TInterfacedObject, IKillerPercentage)
private
fPer: TFloat;
procedure SetPercentage(const Value: TFloat);
function GetPercentage: TFloat;
public
function Kill(Pop : IPopulation): Integer;
// Percentage of population to be killed
property Percentage: TFloat read GetPercentage write SetPercentage;
end;
TParentSelectorTournament = class(TInterfacedObject, IParentSelector)
public
function SelectParent(Population: IPopulation): IIndividual;
end;
TTSPBreederCrossover = class(TInterfacedObject, IBreeder)
public
function BreedOffspring(PSelector: IParentSelector; Pop: IPopulation): IIndividual;
end;
TTSPMutator = class(TInterfacedObject, ITSPMutator)
private
fTrans: TFloat;
fInv: TFloat;
procedure SetInv(const Value: TFloat);
procedure SetTrans(const Value: TFloat);
function GetInv: TFloat;
function GetTrans: TFloat;
public
procedure Mutate(Individual: IIndividual);
published
// Probability of doing a transposition
property Trans read GetTrans write SetTrans;
// Probability of doing an inversion
property Inversion: TFloat read GetInv write SetInv;
end;
TTSPExaminer = class(TInterfacedObject, ITSPExaminer)
private
// The Control component we are associated with
fController: ITSPController;
function GetController: ITSPController;
procedure SetController(const Value: ITSPController);
public
// Returns the fitness of an individual as a real number where 0 => best
function GetFitness(Individual : IIndividual) : TFloat;
property Controller : ITSPController read GetController write SetController;
end;
TPopulation = class(TInterfacedObject, IPopulation)
private
// The population
fPop : TInterfaceList;
// Worker for breeding
fBreeder: IBreeder;
// Worker for killing
fKiller: IKiller;
// Worker for parent selection
fParentSelector: IParentSelector;
// Worker for mutation
fMutator: IMutator;
// Worker for initial creation
fCreator: ICreator;
// Worker for fitness calculation
fExaminer: IExaminer;
// On Change event
FOnChange: TNotifyEvent;
procedure Change;
// Getters and Setters
function GetIndividual(I: Integer): IIndividual;
function GetCount: Integer;
function GetBreeder: IBreeder;
function GetCreator: ICreator;
function GetExaminer: IExaminer;
function GetKiller: IKiller;
function GetMutator: IMutator;
function GetOnChange: TNotifyEvent;
function GetParentSelector: IParentSelector;
procedure SetBreeder(const Value: IBreeder);
procedure SetCreator(const Value: ICreator);
procedure SetExaminer(const Value: IExaminer);
procedure SetKiller(const Value: IKiller);
procedure SetMutator(const Value: IMutator);
procedure SetOnChange(const Value: TNotifyEvent);
procedure SetParentSelector(const Value: IParentSelector);
// not interfaced
procedure DanQuickSort(SortList: TInterfaceList; L, R: Integer; SCompare: TInterfaceCompare);
procedure Sort(Compare: TInterfaceCompare);
protected
// Comparison function for Sort()
function CompareIndividuals(I1, I2: IIndividual): Integer;
// Sort the population
procedure SortPopulation;
public
// The constructor
constructor Create;
// The destructor
destructor Destroy; override;
// Adds an individual to the population
procedure Add(New : IIndividual);
// Deletes an individual from the population
procedure Delete(I : Integer);
// Runs a single generation
procedure Generation;
// Initialise the population
procedure Initialise(Size : Integer);
// Clear ourselves out
procedure Clear;
// Get the fitness of an individual
function FitnessOf(I : Integer) : TFloat;
// Access to the population members
property Pop[I : Integer] : IIndividual read GetIndividual; default;
// The size of the population
property Count : Integer read GetCount;
property ParentSelector : IParentSelector read GetParentSelector write SetParentSelector;
property Breeder : IBreeder read GetBreeder write SetBreeder;
property Killer : IKiller read GetKiller write SetKiller;
property Mutator : IMutator read GetMutator write SetMutator;
property Creator : ICreator read GetCreator write SetCreator;
property Examiner : IExaminer read GetExaminer write SetExaminer;
// An event
property OnChange : TNotifyEvent read GetOnChange write SetOnChange;
end;
TTSPController = class(TInterfacedObject, ITSPController)
private
fXmin, fXmax, fYmin, fYmax: TFloat;
{ The array of 'cities' }
fCities : array of TPoint2D;
{ The array of 'vehicles' }
fVehicles : array of TVehicle;
{ The array of 'vehicle number' }
fNoVehicles : ArrayInt;/////////////////////
{ The number of 'new cities' }
fCityCount: Integer;
{ The number of 'old cities' }
foldCityCount: Integer;
{ The number of 'travelers' }
fTravelCount:Integer; ///////////////////////
{ The number of 'depots' }
fDepotCount:Integer; ///////////////////////
{ Getters... }
function GetCity(I: Integer): TPoint2D;
function GetNoVehicle(I: Integer): TInt;
function GetCityCount: Integer;
function GetOldCityCount: Integer;
function GetTravelCount:Integer;
function GetDepotCount:Integer;
function GetXmax: TFloat;
function GetXmin: TFloat;
function GetYmax: TFloat;
function GetYmin: TFloat;
{ Setters... }
procedure SetCityCount(const Value: Integer);
procedure SetOldCityCount(const Value: Integer);
procedure SetTravelCount(const Value: Integer); /////////////
procedure SetDepotCount(const Value: Integer); /////////////
procedure SetXmax(const Value: TFloat);
procedure SetXmin(const Value: TFloat);
procedure SetYmax(const Value: TFloat);
procedure SetYmin(const Value: TFloat);
function TimeCostBetween(C1, C2: Integer): TFloat;
function GetTimeConstraint(Individual: IIndividual): TInt;
function DateSpanToMin(d1, d2: TDateTime): integer;
function GetVehicleInfo(routeInt: Tint): integer;
procedure writeTimeArray;
procedure writeCostArray;
public
{ The constructor }
constructor Create;
{ The destructor }
destructor Destroy; override;
{ Get the distance between two cities }
function DistanceBetween(C1, C2 : Integer) : TFloat;
{ Get the cost between two cities }
function CostBetween(C1, C2: Integer): TFloat;
function GetWeightConstraint( Individual: IIndividual): TInt;
function GetBackConstraint( Individual: IIndividual): TInt;
{ Places the cities at random points }
procedure RandomCities;
{ Area limits }
property Xmin: TFloat read GetXmin write SetXmin;
property Xmax: TFloat read GetXmax write SetXmax;
property Ymin: TFloat read GetYmin write SetYmin;
property Ymax: TFloat read GetYmax write SetYmax;
{ Properties... }
property CityCount : Integer read GetCityCount write SetCityCount;
property OldCityCount : Integer read GetOldCityCount write SetOldCityCount;
property TravelCount : Integer read GetTravelCount write SetTravelCount; ///////////
property DepotCount : Integer read GetDepotCount write SetDepotCount; ///////////
{ Access to the cities array }
property Cities[I : Integer] : TPoint2D read GetCity;
property NoVehicles[I : Integer] : TInt read GetNoVehicle; ///////////////
end;
implementation
uses
Math;
{ TIndividual }
function TIndividual.GetFitness: TFloat;
begin
result := fFitness;
end;
function TIndividual.GetWeConstrain: integer;
begin
result := fWeConstrain;
end;
function TIndividual.GetBackConstrain: integer;
begin
result := fBackConstrain;
end;
function TIndividual.GetTimeConstrain: integer;
begin
result := fTimeConstrain;
end;
procedure TIndividual.SetBackConstrain(const Value: integer);
begin
fBackConstrain := Value;
end;
procedure TIndividual.SetFitness(const Value: TFloat);
begin
fFitness := Value;
end;
procedure TIndividual.SetWeConstrain(const Value: integer);
begin
fWeConstrain := Value;
end;
procedure TIndividual.SetTimeConstrain(const Value: integer);
begin
fTimeConstrain := Value;
end;
{ TTSPIndividual }
constructor TTSPIndividual.Create(Size: TInt);
begin
Inherited Create;
SetLength(fRouteArray, Size);
// fSteps := Size;
end;
destructor TTSPIndividual.Destroy;
begin
SetLength(fRouteArray, 0);
inherited;
end;
function TTSPIndividual.GetRouteArray(I: Integer): Integer;
begin
result := fRouteArray[I];
end;
function TTSPIndividual.GetSteps: Integer;
begin
result := Length(fRouteArray);
end;
procedure TTSPIndividual.SetSteps(const Value: Integer);
begin
SetLength(fRouteArray, Value);
end;
procedure TTSPIndividual.SetRouteArray(I: Integer; const Value: Integer);
begin
fRouteArray[I] := Value;
end;