模擬退火算法
模擬退火算法來(lái)源于固體退火原理，將固體加溫至充分高，再讓其徐徐冷卻，加溫時(shí)，固體內(nèi)部粒子隨溫升變?yōu)闊o(wú)序狀內(nèi)能增大，而徐徐冷卻時(shí)粒子漸趨有序，在每個(gè)溫度都達(dá)到平衡態(tài)，最后在常溫時(shí)達(dá)到基態(tài)，內(nèi)能減為最小。根據(jù)Metropolis準(zhǔn)則，粒子在溫度T時(shí)趨于平衡的概率為e-ΔE/(kT)，其中E為溫度T時(shí)的內(nèi)能，ΔE為其改變量，k為Boltzmann常數(shù)。用固體退火模擬組合優(yōu)化問(wèn)題，將內(nèi)能E模擬為目標(biāo)函數(shù)值f，溫度T演化成控制參數(shù)t，即得到解組合優(yōu)化問(wèn)題的模擬退火算法：由初始解i和控制參數(shù)初值t開(kāi)始，對(duì)當(dāng)前解重復(fù)“產(chǎn)生新解→計(jì)算目標(biāo)函數(shù)差→接受或舍棄”的迭代，并逐步衰減t值，算法終止時(shí)的當(dāng)前解即為所得近似最優(yōu)解，這是基于蒙特卡羅迭代求解法的一種啟發(fā)式隨機(jī)搜索過(guò)程。退火過(guò)程由冷卻進(jìn)度表(Cooling Schedule)控制，包括控制參數(shù)的初值t及其衰減因子Δt、每個(gè)t值時(shí)的迭代次數(shù)L和停止條件S。
3.5.1 模擬退火算法的模型
模擬退火算法可以分解為解空間、目標(biāo)函數(shù)和初始解三部分。
模擬退火的基本思想:
(1) 初始化：初始溫度T(充分大)，初始解狀態(tài)S(是算法迭代的起點(diǎn))， 每個(gè)T值的迭代次數(shù)L
(2) 對(duì)k=1，……，L做第(3)至第6步：
(3) 產(chǎn)生新解S′
(4) 計(jì)算增量Δt′=C(S′)-C(S)，其中C(S)為評(píng)價(jià)函數(shù)
(5) 若Δt′<0則接受S′作為新的當(dāng)前解，否則以概率exp(-Δt′/T)接受S′作為新的當(dāng)前解.
(6) 如果滿足終止條件則輸出當(dāng)前解作為最優(yōu)解，結(jié)束程序。
終止條件通常取為連續(xù)若干個(gè)新解都沒(méi)有被接受時(shí)終止算法。
(7) T逐漸減少，且T->0，然后轉(zhuǎn)第2步。
算法對(duì)應(yīng)動(dòng)態(tài)演示圖：
模擬退火算法新解的產(chǎn)生和接受可分為如下四個(gè)步驟：
第一步是由一個(gè)產(chǎn)生函數(shù)從當(dāng)前解產(chǎn)生一個(gè)位于解空間的新解；為便于后續(xù)的計(jì)算和接受，減少算法耗時(shí)，通常選擇由當(dāng)
前新解經(jīng)過(guò)簡(jiǎn)單地變換即可產(chǎn)生新解的方法，如對(duì)構(gòu)成新解的全部或部分元素進(jìn)行置換、互換等，注意到產(chǎn)生新解的變換方法
決定了當(dāng)前新解的鄰域結(jié)構(gòu)，因而對(duì)冷卻進(jìn)度表的選取有一定的影響。
第二步是計(jì)算與新解所對(duì)應(yīng)的目標(biāo)函數(shù)差。因?yàn)槟繕?biāo)函數(shù)差僅由變換部分產(chǎn)生，所以目標(biāo)函數(shù)差的計(jì)算最好按增量計(jì)算。
事實(shí)表明，對(duì)大多數(shù)應(yīng)用而言，這是計(jì)算目標(biāo)函數(shù)差的最快方法。
第三步是判斷新解是否被接受,判斷的依據(jù)是一個(gè)接受準(zhǔn)則，最常用的接受準(zhǔn)則是Metropo1is準(zhǔn)則: 若Δt′<0則接受S′作
為新的當(dāng)前解S，否則以概率exp(-Δt′/T)接受S′作為新的當(dāng)前解S。
第四步是當(dāng)新解被確定接受時(shí)，用新解代替當(dāng)前解，這只需將當(dāng)前解中對(duì)應(yīng)于產(chǎn)生新解時(shí)的變換部分予以實(shí)現(xiàn)，同時(shí)修正
目標(biāo)函數(shù)值即可。此時(shí)，當(dāng)前解實(shí)現(xiàn)了一次迭代?？稍诖嘶A(chǔ)上開(kāi)始下一輪試驗(yàn)。而當(dāng)新解被判定為舍棄時(shí)，則在原當(dāng)前解的
基礎(chǔ)上繼續(xù)下一輪試驗(yàn)。
模擬退火算法與初始值無(wú)關(guān)，算法求得的解與初始解狀態(tài)S(是算法迭代的起點(diǎn))無(wú)關(guān)；模擬退火算法具有漸近收斂性，已在
理論上被證明是一種以概率l 收斂于全局最優(yōu)解的全局優(yōu)化算法；模擬退火算法具有并行性。
3.5.2 模擬退火算法的簡(jiǎn)單應(yīng)用
作為模擬退火算法應(yīng)用，討論貨郎擔(dān)問(wèn)題(Travelling Salesman Problem，簡(jiǎn)記為TSP)：設(shè)有n個(gè)城市，用數(shù)碼1,…,n代表
。城市i和城市j之間的距離為d(i，j) i, j=1,…,n．TSP問(wèn)題是要找遍訪每個(gè)域市恰好一次的一條回路，且其路徑總長(zhǎng)度為最
短.。
求解TSP的模擬退火算法模型可描述如下：
解空間 解空間S是遍訪每個(gè)城市恰好一次的所有回路，是{1，……，n}的所有循環(huán)排列的集合，S中的成員記為(w1,w2 ,…
…，wn)，并記wn+1= w1。初始解可選為(1，……，n)
目標(biāo)函數(shù) 此時(shí)的目標(biāo)函數(shù)即為訪問(wèn)所有城市的路徑總長(zhǎng)度或稱為代價(jià)函數(shù)：
我們要求此代價(jià)函數(shù)的最小值。
新解的產(chǎn)生 隨機(jī)產(chǎn)生1和n之間的兩相異數(shù)k和m，若k<m，則將
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
變?yōu)椋?
(w1, w2 ,…，wm , wm-1 ,…，wk+1 , wk ,…，wn).
如果是k>m，則將
(w1, w2 ,…，wk , wk+1 ,…，wm ,…，wn)
變?yōu)椋?
(wm, wm-1 ,…，w1 , wm+1 ,…，wk-1 ,wn , wn-1 ,…，wk).
上述變換方法可簡(jiǎn)單說(shuō)成是“逆轉(zhuǎn)中間或者逆轉(zhuǎn)兩端”。
也可以采用其他的變換方法，有些變換有獨(dú)特的優(yōu)越性，有時(shí)也將它們交替使用，得到一種更好方法。
代價(jià)函數(shù)差 設(shè)將(w1, w2 ,……，wn)變換為(u1, u2 ,……，un), 則代價(jià)函數(shù)差為：
根據(jù)上述分析，可寫出用模擬退火算法求解TSP問(wèn)題的偽程序：
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); { 從當(dāng)前回路S產(chǎn)生新回路S′}
Δt:=f(S′))-f(S);{f(S)為路徑總長(zhǎng)}
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
模擬退火算法的應(yīng)用很廣泛，可以較高的效率求解最大截問(wèn)題(Max Cut Problem)、0-1背包問(wèn)題(Zero One Knapsack
Problem)、圖著色問(wèn)題(Graph Colouring Problem)、調(diào)度問(wèn)題(Scheduling Problem)等等。
3.5.3 模擬退火算法的參數(shù)控制問(wèn)題
模擬退火算法的應(yīng)用很廣泛，可以求解NP完全問(wèn)題，但其參數(shù)難以控制，其主要問(wèn)題有以下三點(diǎn)：
(1) 溫度T的初始值設(shè)置問(wèn)題。
溫度T的初始值設(shè)置是影響模擬退火算法全局搜索性能的重要因素之一、初始溫度高，則搜索到全局最優(yōu)解的可能性大，但
因此要花費(fèi)大量的計(jì)算時(shí)間；反之，則可節(jié)約計(jì)算時(shí)間，但全局搜索性能可能受到影響。實(shí)際應(yīng)用過(guò)程中，初始溫度一般需要
依據(jù)實(shí)驗(yàn)結(jié)果進(jìn)行若干次調(diào)整。
(2) 退火速度問(wèn)題。
模擬退火算法的全局搜索性能也與退火速度密切相關(guān)。一般來(lái)說(shuō)，同一溫度下的“充分”搜索(退火)是相當(dāng)必要的，但這
需要計(jì)算時(shí)間。實(shí)際應(yīng)用中，要針對(duì)具體問(wèn)題的性質(zhì)和特征設(shè)置合理的退火平衡條件。
(3) 溫度管理問(wèn)題。
溫度管理問(wèn)題也是模擬退火算法難以處理的問(wèn)題之一。實(shí)際應(yīng)用中，由于必須考慮計(jì)算復(fù)雜度的切實(shí)可行性等問(wèn)題，常采
用如下所示的降溫方式：
T(t+1)＝k×T(t)
式中k為正的略小于1.00的常數(shù)，t為降溫的次數(shù)
使用SA解決TSP問(wèn)題的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
又一個(gè)相關(guān)的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;
%下面的函數(shù)產(chǎn)生新解
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;
遺傳算法：
旅行商問(wèn)題(traveling saleman problem,簡(jiǎn)稱tsp)：
已知n個(gè)城市之間的相互距離，現(xiàn)有一個(gè)推銷員必須遍訪這n個(gè)城市，并且每個(gè)城市只能訪問(wèn)一次，最后又必須返回出發(fā)城市。如何安排他對(duì)這些城市的訪問(wèn)次序，可使其旅行路線的總長(zhǎng)度最短？
用圖論的術(shù)語(yǔ)來(lái)說(shuō)，假設(shè)有一個(gè)圖 g=(v,e)，其中v是頂點(diǎn)集，e是邊集，設(shè)d=(dij)是由頂點(diǎn)i和頂點(diǎn)j之間的距離所組成的距離矩陣，旅行商問(wèn)題就是求出一條通過(guò)所有頂點(diǎn)且每個(gè)頂點(diǎn)只通過(guò)一次的具有最短距離的回路。
這個(gè)問(wèn)題可分為對(duì)稱旅行商問(wèn)題(dij=dji,,任意i,j=1,2,3，…,n)和非對(duì)稱旅行商問(wèn)題(dij≠dji,,任意i,j=1,2,3，…,n)。
若對(duì)于城市v={v1,v2,v3,…,vn}的一個(gè)訪問(wèn)順序?yàn)閠=(t1,t2,t3,…,ti,…,tn),其中ti∈v(i=1,2,3,…,n)，且記tn+1= t1，則旅行商問(wèn)題的數(shù)學(xué)模型為：
min l=σd(t(i),t(i+1)) （i=1,…,n）
旅行商問(wèn)題是一個(gè)典型的組合優(yōu)化問(wèn)題，并且是一個(gè)np難問(wèn)題，其可能的路徑數(shù)目與城市數(shù)目n是成指數(shù)型增長(zhǎng)的，所以一般很難精確地求出其最優(yōu)解，本文采用遺傳算法求其近似解。
遺傳算法：
初始化過(guò)程：用v1,v2,v3,…,vn代表所選n個(gè)城市。定義整數(shù)pop-size作為染色體的個(gè)數(shù)，并且隨機(jī)產(chǎn)生pop-size個(gè)初始染色體，每個(gè)染色體為1到18的整數(shù)組成的隨機(jī)序列。
適應(yīng)度f(wàn)的計(jì)算：對(duì)種群中的每個(gè)染色體vi，計(jì)算其適應(yīng)度，f=σd(t(i),t(i+1)).
評(píng)
價(jià)函數(shù)eval(vi)：用來(lái)對(duì)種群中的每個(gè)染色體vi設(shè)定一個(gè)概率，以使該染色體被選中的可能性與其種群中其它染色體的適應(yīng)性成比例，既通過(guò)輪盤賭，適
應(yīng)性強(qiáng)的染色體被選擇產(chǎn)生后臺(tái)的機(jī)會(huì)要大，設(shè)alpha∈(0,1)，本文定義基于序的評(píng)價(jià)函數(shù)為eval(vi)=alpha*(1-alpha).^
(i-1) 。[隨機(jī)規(guī)劃與模糊規(guī)劃]
選擇過(guò)程：選擇過(guò)程是以旋轉(zhuǎn)賭輪pop-size次為基礎(chǔ)，每次旋轉(zhuǎn)都為新的種群選擇一個(gè)染色體。賭輪是按每個(gè)染色體的適應(yīng)度進(jìn)行選擇染色體的。
step1 、對(duì)每個(gè)染色體vi,計(jì)算累計(jì)概率qi，q0=0;qi=σeval(vj) j=1,…,i;i=1,…pop-size.
step2、從區(qū)間(0,pop-size)中產(chǎn)生一個(gè)隨機(jī)數(shù)r；
step3、若qi-1<r<qi,則選擇第i個(gè)染色體 ；
step4、重復(fù)step2和step3共pop-size次，這樣可以得到pop-size個(gè)復(fù)制的染色體。
grefenstette
編碼：由于常規(guī)的交叉運(yùn)算和變異運(yùn)算會(huì)使種群中產(chǎn)生一些無(wú)實(shí)際意義的染色體，本文采用grefenstette編碼《遺傳算法原理及應(yīng)用》可以避免這種情
況的出現(xiàn)。所謂的grefenstette編碼就是用所選隊(duì)員在未選（不含淘汰）隊(duì)員中的位置，如：
8 15 2 16 10 7 4 3 11 14 6 12 9 5 18 13 17 1
對(duì)應(yīng)：
8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1。
交叉過(guò)程：本文采用常規(guī)單點(diǎn)交叉。為確定交叉操作的父代，從 到pop-size重復(fù)以下過(guò)程：從[0，1]中產(chǎn)生一個(gè)隨機(jī)數(shù)r，如果r<pc ，則選擇vi作為一個(gè)父代。
將所選的父代兩兩組隊(duì)，隨機(jī)產(chǎn)生一個(gè)位置進(jìn)行交叉，如：
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
變
異過(guò)程：本文采用均勻多點(diǎn)變異。類似交叉操作中選擇父代的過(guò)程，在r<pm
的標(biāo)準(zhǔn)下選擇多個(gè)染色體vi作為父代。對(duì)每一個(gè)選擇的父代，隨機(jī)選擇多個(gè)位置，使其在每位置按均勻變異（該變異點(diǎn)xk的取值范圍為[ukmin,
ukmax],產(chǎn)生一個(gè)[0，1]中隨機(jī)數(shù)r，該點(diǎn)變異為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編碼之后進(jìn)行的，為了循環(huán)操作和返回最終結(jié)果，必須逆grefenstette編碼過(guò)程，將編碼恢復(fù)到自然編碼。
循環(huán)操作：判斷是否滿足設(shè)定的帶數(shù)xzome，否，則跳入適應(yīng)度f(wàn)的計(jì)算；是，結(jié)束遺傳操作，跳出。
Matlab程序：
function [bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%
%————————————————————————
%[bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%d:距離矩陣
%termops:種群代數(shù)
%num:每代染色體的個(gè)數(shù)
%pc:交叉概率
%cxops:由于本程序采用單點(diǎn)交叉，交叉點(diǎn)的設(shè)置在本程序中沒(méi)有很好的解決，所以本文了采用定點(diǎn)，即第cxops，可以隨機(jī)產(chǎn)生。
%pm:變異概率
%alpha:評(píng)價(jià)函數(shù)eval(vi)=alpha*(1-alpha).^(i-1).
%bestpop:返回的最優(yōu)種群
%trace:進(jìn)化軌跡
%------------------------------------------------
%####@@@##版權(quán)所有！歡迎廣大網(wǎng)友改正，改進(jìn)！##@@@####
%e-mail:tobysidney33@sohu.com
%####################################################
%
citynum=size(d,2);
n=nargin;
if n<2
disp('缺少變量！！')
disp('^_^開(kāi)個(gè)玩笑^_^')
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
-------------------------------------------------
又一個(gè)Matlab程序，其中交叉算法采用的是由Goldberg和Lingle于1985年提出的PMX(部分匹配交叉)，淘汰保護(hù)指數(shù)alpha是我自己設(shè)計(jì)的，起到了加速優(yōu)勝劣汰的作用。
%TSP問(wèn)題（又名：旅行商問(wèn)題，貨郎擔(dān)問(wèn)題）遺傳算法通用matlab程序
%D是距離矩陣，n為種群個(gè)數(shù),建議取為城市個(gè)數(shù)的1~2倍，
%C為停止代數(shù)，遺傳到第 C代時(shí)程序停止,C的具體取值視問(wèn)題的規(guī)模和耗費(fèi)的時(shí)間而定
%m為適應(yīng)值歸一化淘汰加速指數(shù) ,最好取為1,2,3,4 ,不宜太大
%alpha為淘汰保護(hù)指數(shù),可取為0~1之間任意小數(shù),取1時(shí)關(guān)閉保護(hù)功能,最好取為0.8~1.0
%R為最短路徑,Rlength為路徑長(zhǎng)度
function [R,Rlength]=geneticTSP(D,n,C,m,alpha)
[N,NN]=size(D);
farm=zeros(n,N);%用于存儲(chǔ)種群
for i=1:n
farm(i,:)=randperm(N);%隨機(jī)生成初始種群
end
R=farm(1,:);%存儲(chǔ)最優(yōu)種群
len=zeros(n,1);%存儲(chǔ)路徑長(zhǎng)度
fitness=zeros(n,1);%存儲(chǔ)歸一化適應(yīng)值
counter=0;
while counter<C
for i=1:n
len(i,1)=myLength(D,farm(i,:));%計(jì)算路徑長(zhǎng)度
end
maxlen=max(len);
minlen=min(len);
fitness=fit(len,m,maxlen,minlen);%計(jì)算歸一化適應(yīng)值
rr=find(len==minlen);
R=farm(rr(1,1),:);%更新最短路徑
FARM=farm;%優(yōu)勝劣汰，nn記錄了復(fù)制的個(gè)數(shù)
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,:);%保持種群規(guī)模為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);%隨機(jī)選擇交叉范圍，從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;
% 計(jì)算路徑的子程序
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
%計(jì)算歸一化適應(yīng)值子程序
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
一個(gè)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<<"請(qǐng)輸入該圖的節(jié)點(diǎn)數(shù)："<<endl;
int vertices;
cin>>vertices;
MGraph<float> b(vertices,noEdge);
cout<<"請(qǐng)輸入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<<"請(qǐng)輸入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<<"請(qǐng)輸入任意一個(gè)初始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<<"最短路徑長(zhǎng)度："<<endl;
cout<<total;
}
C語(yǔ)言程序：
#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;
}
改進(jìn)后用來(lái)求解VRP問(wèn)題的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;