[分享]从网上找到的一些解决TSP问题的算法及源代码
<P>模拟退火算法模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，</P>
<P>内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据Metropolis</P>
<P>准则，粒子在温度T时趋于平衡的概率为e-ΔE/(kT)，其中E为温度T时的内能，ΔE为其改变量，k为Boltzmann常数。用固体退</P>
<P>火模拟组合优化问题，将内能E模拟为目标函数值f，温度T演化成控制参数t，即得到解组合优化问题的模拟退火算法：由初始</P>
<P>解i和控制参数初值t开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t值，算法终止时的</P>
<P>当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(Cooling </P>
<P>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步。
算法对应动态演示图：
模拟退火算法新解的产生和接受可分为如下四个步骤：
第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当</P>
<P>前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法</P>
<P>决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。
第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。</P>
<P>事实表明，对大多数应用而言，这是计算目标函数差的最快方法。
第三步是判断新解是否被接受,判断的依据是一个接受准则，最常用的接受准则是Metropo1is准则: 若Δt′<0则接受S′作</P>
<P>为新的当前解S，否则以概率exp(-Δt′/T)接受S′作为新的当前解S。
第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正</P>
<P>目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的</P>
<P>基础上继续下一轮试验。
模拟退火算法与初始值无关，算法求得的解与初始解状态S(是算法迭代的起点)无关；模拟退火算法具有渐近收敛性，已在</P>
<P>理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。 </P>
<P>3.5.2 模拟退火算法的简单应用
作为模拟退火算法应用，讨论货郎担问题(Travelling Salesman Problem，简记为TSP)：设有n个城市，用数码1,…,n代表</P>
<P>。城市i和城市j之间的距离为d(i，j) i, j=1,…,n．TSP问题是要找遍访每个域市恰好一次的一条回路，且其路径总长度为最</P>
<P>短.。
求解TSP的模拟退火算法模型可描述如下：
解空间 解空间S是遍访每个城市恰好一次的所有回路，是{1，……，n}的所有循环排列的集合，S中的成员记为(w1,w2 ,…</P>
<P>…，wn)，并记wn+1= w1。初始解可选为(1，……，n)
目标函数 此时的目标函数即为访问所有城市的路径总长度或称为代价函数： </P>
<P> 我们要求此代价函数的最小值。
新解的产生 随机产生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), 则代价函数差为： </P>
<P>根据上述分析，可写出用模拟退火算法求解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-)
S=S′;
IF the-halt-condition-is-TRUE THEN
termination=true;
End;
T_lower;
End;
End
模拟退火算法的应用很广泛，可以较高的效率求解最大截问题(Max Cut Problem)、0-1背包问题(Zero One Knapsack </P>
<P>Problem)、图着色问题(Graph Colouring Problem)、调度问题(Scheduling Problem)等等。 </P>
<P>3.5.3 模拟退火算法的参数控制问题
模拟退火算法的应用很广泛，可以求解NP完全问题，但其参数难以控制，其主要问题有以下三点：
(1) 温度T的初始值设置问题。
温度T的初始值设置是影响模拟退火算法全局搜索性能的重要因素之一、初始温度高，则搜索到全局最优解的可能性大，但</P>
<P>因此要花费大量的计算时间；反之，则可节约计算时间，但全局搜索性能可能受到影响。实际应用过程中，初始温度一般需要</P>
<P>依据实验结果进行若干次调整。
(2) 退火速度问题。
模拟退火算法的全局搜索性能也与退火速度密切相关。一般来说，同一温度下的“充分”搜索(退火)是相当必要的，但这</P>
<P>需要计算时间。实际应用中，要针对具体问题的性质和特征设置合理的退火平衡条件。
(3) 温度管理问题。
温度管理问题也是模拟退火算法难以处理的问题之一。实际应用中，由于必须考虑计算复杂度的切实可行性等问题，常采</P>
<P>用如下所示的降温方式： </P>
<P>T(t+1)＝k×T(t)
式中k为正的略小于1.00的常数，t为降温的次数 </P>
<P>使用SA解决TSP问题的Matlab程序：</P>
<DIV class=HtmlCode>
<P>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 (<a href="mailtjmyers@nwu" target="_blank" >jmyers@nwu</A>.
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 + ));
% 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(,)'.^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 = ;
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 = ;
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(' = [%f %f]\n', minE, maxE);
fprintf(' = [%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 </P></DIV>
<P>又一个相关的Matlab程序</P>
<DIV class=HtmlCode>
<P>function=zkp(w,c,M,t0,tf)
=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 =exchange_2(pi0,d)
=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; </P></DIV>
遗传算法GA
<P>遗传算法：</P><P>旅行商问题(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重复以下过程：从中产生一个随机数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的取值范围为,产生一个中随机数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的计算；是，结束遗传操作，跳出。</P>
<P>Matlab程序：</P>
<DIV class=HtmlCode>
<P>function =ga(d,termops,num,pc,cxops,pm,alpha)
%
%————————————————————————
%=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</P>
<P>=initializega(num,citynum);
for i=1:termops
=f(d,t);
=find(l==max(l));
trace(i)=-l(y(1));
bestpop=t(y(1),:);
=select(t,l,alpha);
=grefenstette(t);
=crossover(g,pc,cxops);
=mutation(g1,pm); %均匀变异
=congrefenstette(g);
end</P>
<P>---------------------------------------------------------
function =initializega(num,citynum)
for i=1:num
t(i,:)=randperm(citynum);
end
-----------------------------------------------------------
function =f(d,t)
=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 =select(t,l,alpha)
=size(l);
t1=t;
=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 =grefenstette(t)
=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 =crossover(g,pc,cxops)
=size(g);
ran=rand(1,m);
r=cxops;
=find(ran<pc);
if ru>=2
for k=1:2:length(ru)-1
g1(ru(k),:)=),g(ru(k+1),[(r+1):n])];
g(ru(k+1),:)=),g(ru(k),[(r+1):n])];
g(ru(k),:)=g1(ru(k),:);
end
end
--------------------------------------------
function =mutation(g,pm) %均匀变异
=size(g);
ran=rand(1,m);
r=rand(1,3); %dai gai jin
rr=floor(n*rand(1,3)+1);
=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 =congrefenstette(g)
=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
------------------------------------------------- </P></DIV>
<P>又一个Matlab程序，其中交叉算法采用的是由Goldberg和Lingle于1985年提出的PMX(部分匹配交叉)，淘汰保护指数alpha是我自己设计的，起到了加速优胜劣汰的作用。</P>
<DIV class=HtmlCode>
<P>%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 =geneticTSP(D,n,C,m,alpha)</P>
<P>=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;</P>
<P>while counter<C</P>
<P>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),:);%更新最短路径</P>
<P>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,:);</P>
<P>=size(FARM);%交叉和变异
while aa<n
if nn<=2
nnper=randperm(2);
else
nnper=randperm(nn);
end
A=FARM(nnper(1),:);
B=FARM(nnper(2),:);
=intercross(A,B);
FARM=;
=size(FARM);
end
if aa>n
FARM=FARM(1:n,:);%保持种群规模为n
end</P>
<P>farm=FARM;
clear FARM
counter=counter+1</P>
<P>end</P>
<P>Rlength=myLength(D,R);</P>
<P>function =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));
=exchange(a(1,p+i-1),b(1,p+i-1));
=exchange(a(1,x),b(1,y));
end
function =exchange(x,y)
temp=x;
x=y;
y=temp;</P>
<P>% 计算路径的子程序
function len=myLength(D,p)
=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</P>
<P>%计算归一化适应值子程序
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 </P></DIV>
<P>一个C++的程序：</P>
<DIV class=HtmlCode>
<P>//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* ;
for(int i=0;i<n;i++){
a=new T;
a=0;
for(int j=0;j<n;j++)if(i!=j)a=NoEdge;
}
}
template<class T>
MGraph<T>::~MGraph()
{
for(int i=0;i<n;i++)delete[]a;
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==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!=NoEdge){
cerr<<"BadInput!"<<endl;
return false;
}
a=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==NoEdge){
cerr<<"BadInput!"<<endl;
return false;
}
a=NoEdge;
e--;
return true;
}
template<class T>
void MGraph<T>::Floyd(T**& d,int**& path)
{
d=new T* ;
path=new int* ;
for(int i=0;i<n;i++){
d=new T;
path=new int;
for(int j=0;j<n;j++){
d=a;
if(i!=j&&a<NoEdge)path=i;
else path=-1;
}
}
for(int k=0;k<n;k++){
for(i=0;i<n;i++)
for(int j=0;j<n;j++)
if(d+d<d){
d=d+d;
path=path;
}
}
}
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<<' ';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<<' ';
if(j==vertices-1)cout<<endl;
}
}
int *V;
V=new int;
cout<<"请输入任意一个初始H-圈："<<endl;
for(int n=0;n<=vertices;n++){
cin>>V;
}
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]]+D]]<D]]+D]]){
int l;
l=V;V=V;V=l;
}
}
}
}
}
float total=0;
cout<<"最小回路："<<endl;
for(i=0;i<=vertices;i++){
cout<<V+1<<' ';
}
cout<<endl;
for(i=0;i<vertices;i++)
total+=D]];
cout<<"最短路径长度："<<endl;
cout<<total;
} </P></DIV>
<P>C语言程序：</P>
<DIV class=HtmlCode>
<P>#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></P>
<P>#define maxpop 100
#define maxstring 100</P>
<P>
struct pp{unsigned char chrom;
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;
unsigned char x,y;
float *dd,ff,maxdd,refpd,fm;
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();</P>
<P>main()
{unsigned int gen,k,j,tt;
char fname;
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.chrom='\0';
for(k=0;k<maxpop;k++) newpop.chrom='\0';
printf("Enter Result Data Filename:");
gets(fname);
if((fp=fopen(fname,"w+"))==NULL)
{printf("cannot open file\n");exit(0);}</P>
<P>
gen=0;
randomize();
initialize();</P>
<P>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);
}
fprintf(fp,"\n Y site:\n");
for(k=0;k<lchrom;k++)
{if((k%16)==0) fprintf(fp,"\n");
fprintf(fp,"%5d",y);
}
fprintf(fp,"\n");</P>
<P>
crtinit();
statistics(oldpop);
report(gen,oldpop);
getch();
maxold=min;
fm=100.0*oldpop.x/ff;
do {
gen=gen+1;
generation();
statistics(oldpop);
if(max>maxold)
{maxold=max;
co_min=0;
}
fm=100.0*oldpop.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=100.0*oldpop.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));</P>
<P>getch();
for(k=0;k<lchrom;k++)
{if((k%16)==0)fprintf(fp,"\n");
fprintf(fp,"%5d",oldpop.chrom);
}
fprintf(fp,"\n");</P>
<P>fclose(fp);
farfree(dd);
farfree(p1);
farfree(oldpop);
farfree(newpop);
restorecrtmode();
exit(0);
}</P>
<P>/*%%%%%%%%%%%%%%%%*/</P>
<P>float objfunc(float x1)
{float y;
y=100.0*ff/x1;
return y;
}</P>
<P>/*&&&&&&&&&&&&&&&&&&&*/</P>
<P>void statistics(pop)
struct pp *pop;
{int j;
sumfitness=pop.fitness;
min=pop.fitness;
max=pop.fitness;
maxpp=0;
minpp=0;
for(j=1;j<popsize;j++)
{sumfitness=sumfitness+pop.fitness;
if(pop.fitness>max)
{max=pop.fitness;
maxpp=j;
}
if(pop.fitness<min)
{min=pop.fitness;
minpp=j;
}
}</P>
<P>avg=sumfitness/(float)popsize;
}</P>
<P>/*%%%%%%%%%%%%%%%%%%%%*/</P>
<P>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.chrom,oldpop.chrom,j);
newpop.x=(float)decode(newpop.chrom);
newpop.fitness=objfunc(newpop.x);
newpop.parent1=mate1;
newpop.parent2=mate2;
newpop.xsite=jcross;
newpop.x=(float)decode(newpop.chrom);
newpop.fitness=objfunc(newpop.x);
newpop.parent1=mate1;
newpop.parent2=mate2;
newpop.xsite=jcross;
if(newpop.fitness>min)
{for(k=0;k<lchrom;k++)
oldpop.chrom=newpop.chrom;
oldpop.x=newpop.x;
oldpop.fitness=newpop.fitness;
co_min++;
return;
}</P>
<P> if(newpop.fitness>min)
{for(k=0;k<lchrom;k++)
oldpop.chrom=newpop.chrom;
oldpop.x=newpop.x;
oldpop.fitness=newpop.fitness;
co_min++;
return;
}
j=j+2;
}while(j<popsize);
}</P>
<P>/*%%%%%%%%%%%%%%%%%*/</P>
<P>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;
}</P>
<P>/*%%%%%%%%%%%%%%%%%%%%*/</P>
<P>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);
}</P>
<P>
void initpop()
{unsigned char j1;
unsigned int k5,i1,i2,j,i,k,j2,j3,j4,p5;
float f1,f2;
j=0;
for(k=0;k<lchrom;k++)
oldpop.chrom=k;
for(k=0;k<lchrom;k++)
p5=oldpop.chrom;
randomize();
for(;j<popsize;j++)
{j2=random(lchrom);
for(k=0;k<j2+20;k++)
{j3=random(lchrom);
j4=random(lchrom);
j1=p5;
p5=p5;
p5=j1;
}
for(k=0;k<lchrom;k++)
oldpop.chrom=p5;
}
for(k=0;k<lchrom;k++)
for(j=0;j<lchrom;j++)
dd=hypot(x-x,y-y);
for(j=0;j<popsize;j++)
{oldpop.x=(float)decode(oldpop.chrom);
oldpop.fitness=objfunc(oldpop.x);
oldpop.parent1=0;
oldpop.parent2=0;
oldpop.xsite=0;
}
}</P>
<P>/*&&&&&&&&&&&&&&&&&*/
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=rand();
if(x>maxx)maxx=x;
if(x<minx)minx=x;
y=rand();
if(y>maxy)maxy=y;
if(y<miny)miny=y;
}
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=hypot(x-x,y-y);
if(maxdd<dd)maxdd=dd;
}
refpd=dd;
for(k=0;k<lchrom;k++)
refpd=refpd+dd;
for(j=0;j<lchrom;j++)
dd=4.0*maxdd;
ff=(0.765*maxxy*pow(lchrom,0.5));
minpp=0;
min=dd;
for(j=0;j<lchrom-1;j++)
{if(dd<min)
{min=dd;
minpp=j;
}
}
initpop();
statistics(oldpop);
initreport();
}</P>
<P>/*&&&&&&&&&&&&&&&&&&*/</P>
<P>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.x/maxxy,pop.x,ff);
printf("Co_minpath:%6.4f Maxfit:%10.8f"
,100.0*pop.x/ff,pop.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.chrom];
iy=y.chrom]+110;
jx=x.chrom];
jy=y.chrom]+110;
line(ix,iy,jx,jy);
putpixel(ix,iy,RED);
}
ix=x.chrom];
iy=y.chrom]+110;
jx=x.chrom];
jy=y.chrom]+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/3;
jx=ix+1;
jy=366-fm/3;
line(ix,iy,jx,jy);
putpixel(ix,iy,RED);
}
printf("GEN:%3d",l);
printf("Minpath:%f Maxfit:%f",pop.x,pop.fitness);
printf("Clock:%2d:%2d:%4.2f\n",tt/60,tt%60,ttt-tt*60.0);
}</P>
<P>/*###############*/</P>
<P>float decode(unsigned char *pp)
{int j,k,l;
float tt;
tt=dd*lchrom+pp];
for(j=0;j<lchrom-1;j++)
{tt=tt+dd*lchrom+pp];}
l=0;
for(k=0;k<lchrom-1;k++)
for(j=k+1;j<lchrom;j++)
{if(pp==pp)l++;}
return tt+4*l*maxdd;
}</P>
<P>/*%%%%%%%%%%%%%%%%%%*/
void crtinit()
{int driver,mode;
struct palettetype p;
driver=DETECT;
mode=0;
initgraph(&driver,&mode,"");
cleardevice();
}</P>
<P>/*$$$$$$$$$$$$$$$$$$$$*/
int select()
{double rand1,partsum;
float r1;
int j;
partsum=0.0;
j=0;
rand1=random1()*sumfitness;
do{
partsum=partsum+oldpop.fitness;
j=j+1;
}while((partsum<rand1)&&(j<popsize));
return j-1;
}</P>
<P>/*$$$$$$$$$$$$$$$*/
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,ts2;
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=parent1;
ts2=parent2;
k++;
}
j3=k;
for(j=0;j<lchrom;j++)
{j2=0;
while((parent2!=ts1)&&(j2<k)){j2++;}
if(j2==k)
{ts1=parent2;
j3++;
}
}
j3=k;
for(j=0;j<lchrom;j++)
{j2=0;
while((parent1!=ts2)&&(j2<k)){j2++;}
if(j2==k)
{ts2=parent1;
j3++;
}
}
for(j=0;j<lchrom;j++)
{newpop.chrom=ts1;
newpop.chrom=ts2;
}
}
else
{for(j=0;j<lchrom;j++)
{newpop.chrom=parent1;
newpop.chrom=parent2;
}
mutate=flip(pmutation);
if(mutate)
{s1=1;
nmutation=nmutation+1;
for(j3=0;j3<200;j3++)
{j1=random(lchrom);
j=random(lchrom);
jj=newpop.chrom;
newpop.chrom=newpop.chrom;
newpop.chrom=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.chrom;
newpop.chrom=newpop.chrom;
newpop.chrom=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.chrom+newpop.chrom];
f1=f1+dd.chrom[(i1+1)%lchrom]+
newpop.chrom[(i2+1)%lchrom]];
f2=dd.chrom+
newpop.chrom[(i1+1)%lchrom]];
f2=f2+dd.chrom+
newpop.chrom[(i2+1)%lchrom]];
if(f1<f2){inversion(i1,i2,newpop.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.chrom+newpop.chrom];
f1=f1+dd.chrom[(i1+1)%lchrom]+
newpop.chrom[(i2+1)%lchrom]];
f2=dd.chrom+
newpop.chrom[(i1+1)%lchrom]];
f2=f2+dd.chrom+
newpop.chrom[(i2+1)%lchrom]];
if(f1<f2){inversion(i1,i2,newpop.chrom);}
}
return 1;
}</P>
<P>/*$$$$$$$$$$$$$$$*/</P>
<P>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;
ss=ss;
ss=tt;
}
}</P>
<P>/*%%%%%%%%%%%%%%%*/</P>
<P>void randomize1()
{int i;
randomize();
for(i=0;i<lchrom;i++)
oldrand=random(30001)/30000.0;
jrand=0;
}</P>
<P>/*%%%%%%%%%%%*/</P>
<P>float random1()
{jrand=jrand+1;
if(jrand>=lchrom)
{jrand=0;
randomize1();
}
return oldrand;
}</P>
<P>/*%%%%%%%%%%*/</P>
<P>int flip(float probability)
{float ppp;
ppp=random(20001)/20000.0;
if(ppp<=probability)return 1;
return 0;
}</P></DIV>
<P>改进后用来求解VRP问题的Delphi程序：</P>
<DIV class=HtmlCode>
<P>unit uEA;</P>
<P>interface</P>
<P>uses
uUtilsEA, uIEA, uITSP, Classes, GaPara, windows, SysUtils, fEA_TSP;</P>
<P>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;</P>
<P>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 : 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;</P>
<P>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;</P>
<P>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;</P>
<P>TParentSelectorTournament = class(TInterfacedObject, IParentSelector)
public
function SelectParent(Population: IPopulation): IIndividual;
end;</P>
<P>TTSPBreederCrossover = class(TInterfacedObject, IBreeder)
public
function BreedOffspring(PSelector: IParentSelector; Pop: IPopulation): IIndividual;
end;</P>
<P>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 Transposition: TFloat read GetTrans write SetTrans;
// Probability of doing an inversion
property Inversion: TFloat read GetInv write SetInv;
end;</P>
<P>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;</P>
<P>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 : 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;</P>
<P>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;</P>
<P>function GetWeightConstraint( Individual: IIndividual): TInt;</P>
<P>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 : TPoint2D read GetCity;
property NoVehicles : TInt read GetNoVehicle; ///////////////
end;</P>
<P>implementation</P>
<P>uses
Math;</P>
<P>{ TIndividual }</P>
<P>function TIndividual.GetFitness: TFloat;
begin
result := fFitness;
end;</P>
<P>function TIndividual.GetWeConstrain: integer;
begin
result := fWeConstrain;
end;</P>
<P>function TIndividual.GetBackConstrain: integer;
begin
result := fBackConstrain;
end;</P>
<P>function TIndividual.GetTimeConstrain: integer;
begin
result := fTimeConstrain;
end;</P>
<P>procedure TIndividual.SetBackConstrain(const Value: integer);
begin
fBackConstrain := Value;
end;</P>
<P>procedure TIndividual.SetFitness(const Value: TFloat);
begin
fFitness := Value;
end;</P>
<P>procedure TIndividual.SetWeConstrain(const Value: integer);
begin
fWeConstrain := Value;
end;</P>
<P>procedure TIndividual.SetTimeConstrain(const Value: integer);
begin
fTimeConstrain := Value;
end;</P>
<P>{ TTSPIndividual }</P>
<P>constructor TTSPIndividual.Create(Size: TInt);
begin
Inherited Create;
SetLength(fRouteArray, Size);
// fSteps := Size;
end;</P>
<P>destructor TTSPIndividual.Destroy;
begin
SetLength(fRouteArray, 0);
inherited;
end;</P>
<P>function TTSPIndividual.GetRouteArray(I: Integer): Integer;
begin
result := fRouteArray;
end;</P>
<P>function TTSPIndividual.GetSteps: Integer;
begin
result := Length(fRouteArray);
end;</P>
<P>procedure TTSPIndividual.SetSteps(const Value: Integer);
begin
SetLength(fRouteArray, Value);
end;</P>
<P>procedure TTSPIndividual.SetRouteArray(I: Integer; const Value: Integer);
begin
fRouteArray := Value;
end;</P>
<P>function TTSPIndividual.GetWeConstrain: integer;
begin
result := fWeConstrain;
end;</P>
<P>function TTSPIndividual.GetBackConstrain: integer;
begin
result := fBackConstrain;
end;</P>
<P>function TTSPIndividual.GetTimeConstrain: integer;
begin
result := fTimeConstrain;
end;</P>
<P>procedure TTSPIndividual.SetWeConstrain(const Value: integer);
begin
fWeConstrain := Value;
end;</P>
<P>procedure TTSPIndividual.SetBackConstrain(const Value: integer);
begin
fBackConstrain := Value;
end;</P>
<P>procedure TTSPIndividual.SetTimeConstrain(const Value: integer);
begin
fTimeConstrain := Value;
end;</P>
<P>{ TTSPCreator }</P>
<P>function TTSPCreator.CreateIndividual: IIndividual;
var
New: ITSPIndividual;
i, j, Top, Temp : Integer;
//trav:integer;
begin
// Get the number of cities
Top := fController.CityCount;
// Create the new individual
New := TTSPIndividual.Create(Top);
// Initialise it with a sequential route
for i := 0 to Top - 1 do
New.RouteArray := i;
// Shuffle the route
for i := Top - 1 downto 1 do
begin
j := Random(i);
Temp := New.RouteArray;
New.RouteArray := New.RouteArray;
New.RouteArray := Temp;
end;
result := New;
end;</P>
<P>function TTSPCreator.CreateFeasibleIndividual: IIndividual;
var
New: ITSPIndividual;
i, j, Top, Temp : Tint;
Msg:TMsg;
begin
// Get the number of cities
Top := fController.CityCount;
// Create the new individual
New := TTSPIndividual.Create(Top);
// Initialise it with a sequential route
repeat
begin//////////////////////////////////
for i := 0 to Top - 1 do
New.RouteArray := i;
// Shuffle the route
for i := Top - 1 downto 1 do
begin
j := Random(i);
Temp := New.RouteArray;
New.RouteArray := New.RouteArray;
New.RouteArray := Temp;
end;
//process message sequence//////////
while PeekMessage(Msg,0,0,0,1) do///
begin ///
if Msg.Message<>18 then ///
begin ///
TranslateMessage(Msg); ///
DispatchMessage(Msg); ///
end; ///
end; ///
////////////////////////////////////
end
until (fController.GetWeightConstraint(New)=0)and(fController.GetBackConstraint(New)=0);</P>
<P>result := New;
end;</P>
<P>function TTSPCreator.GetController: ITSPController;
begin
result := fController;
end;</P>
<P>procedure TTSPCreator.SetController(const Value: ITSPController);
begin
fController := Value;
end;</P>
<P>{ TKillerPercentage }</P>
<P>function TKillerPercentage.GetPercentage: TFloat;
begin
result := fPer;
end;</P>
<P>function TKillerPercentage.Kill(Pop: IPopulation): Integer;
var
KillCount, i : Integer;
begin
// Work out the number we have to kill
KillCount := Floor(Pop.Count * (fPer / 100));
// Delete the worst individuals - assuming the population is sorted
for i := 1 to KillCount do
Pop.Delete(Pop.Count - 1);
// Return the number killed
Result := KillCount;
end;</P>
<P>procedure TKillerPercentage.SetPercentage(const Value: TFloat);
begin
fPer := Value;
end;</P>
<P>{ TParentSelectorTournament }</P>
<P>function TParentSelectorTournament.SelectParent(
Population: IPopulation): IIndividual;
var
i1, i2 : Integer;
begin
// Select a random individual
i1 := Random(Population.Count);
// Select a *different* random individual
repeat
i2 := Random(Population.Count);
until i1 <> i2;
// Hold the tournament and return the fittest of the two
if Population.FitnessOf(i1) < Population.FitnessOf(i2) then
Result := Population
else
Result := Population;
end;</P>
<P>{ TTSPBreederCrossover }</P>
<P>function TTSPBreederCrossover.BreedOffspring(PSelector: IParentSelector;
Pop: IPopulation): IIndividual;
var
Child, Mom, Dad, Parent1, Parent2 : ITSPIndividual;
i, j, p : Integer;</P>
<P>function AlreadyAssigned(City, x : Integer) : Boolean;
var
y : Integer;
Found : Boolean;
begin
Found := False;
for y := 0 to x - 1 do
begin
if Child.RouteArray = City then
begin
Found := True;
Break;
end;
end;
Result := Found;
end;</P>
<P>begin
// Select a some parents...
Mom := PSelector.SelectParent(Pop) as ITSPIndividual;
Dad := PSelector.SelectParent(Pop) as ITSPIndividual;
// Create a child
Child := TTSPIndividual.Create(Mom.Steps);
// Copy the route from parents to child
for i := 0 to Child.Steps - 1 do
begin
// Choose a parent at random
p := Random(2);
if p = 0 then
begin
Parent1 := Mom;
Parent2 := Dad;
end else
begin
Parent1 := Dad;
Parent2 := Mom;
end;
if not AlreadyAssigned(Parent1.RouteArray, i) then
begin
// Use city from Parent 1 unless used already
Child.RouteArray := Parent1.RouteArray;
end else
if not AlreadyAssigned(Parent2.RouteArray, i) then
begin
// Otherwise use city from Parent 2 unless used already
Child.RouteArray := Parent2.RouteArray;
end else
begin
// If both assigned already then use a random city
repeat
j := Random(Child.Steps);
until not AlreadyAssigned(j, i);
Child.RouteArray := j;
end;
end;
// Return the child
Result := Child;
end;</P>
<P>{ TTSPMutator }</P>
<P>function TTSPMutator.GetInv: TFloat;
begin
result := fInv;
end;</P>
<P>function TTSPMutator.GetTrans: TFloat;
begin
result := fTrans;
end;</P>
<P>procedure TTSPMutator.Mutate(Individual: IIndividual);
var
P: Double;
i, j, t : Integer; Start, Finish : Integer;
begin
with Individual as ITSPIndividual do
begin
// Should we do an inversion?
P := Random * 100;
if P < FTrans then
begin
// Do an inversion (i.e. swap two cities at random)
// Choose first city
i := Random(Steps);
// Choose a second city
repeat
j := Random(Steps);
until i <> j;
// Swap them over
t := RouteArray;
RouteArray := RouteArray;
RouteArray := t;
end;
// Should we do a transposition?
P := Random * 100;
if P < FInv then
begin
// Do a transposition (i.e. reverse a sub-route)
// Choose random start and finish points
Start := Random(Steps - 1);
Finish := Start + Random(Steps - Start);
// Reverse the sub-route
for i := 0 to Floor((Finish - Start) / 2) do
begin
t := RouteArray;
RouteArray := RouteArray;
RouteArray := t;
end;
end;
end;
end;</P>
<P>procedure TTSPMutator.SetInv(const Value: TFloat);
begin
fInv := Value;
end;</P>
<P>procedure TTSPMutator.SetTrans(const Value: TFloat);
begin
fTrans := Value;
end;</P>
<P>{ TTSPExaminer }</P>
<P>function TTSPExaminer.GetController: ITSPController;
begin
result := fController;
end;</P>
<P>function TTSPExaminer.GetFitness(Individual: IIndividual): TFloat;
var
i , weightConstraint, backConstraint : TInt;
Distance , penaltyW, penaltyB : TFloat;
Indi : ITSPIndividual;
begin
Indi := Individual as ITSPIndividual;
Distance := 0;
penaltyW:=FormGaPara.EditWeightConstrain.Value;
penaltyB:=FormGaPara.EditBackConstrain.Value;
for i := 0 to Indi.Steps - 2 do
begin
Distance := Distance + fController.DistanceBetween(Indi.RouteArray, Indi.RouteArray);
end;
Distance := Distance + fController.DistanceBetween(Indi.RouteArray, Indi.RouteArray);
WeightConstraint:=fController.GetWeightConstraint(Indi);
backConstraint:=fController.GetBackConstraint(Indi);
Indi.WeConstrain:=WeightConstraint;
Indi.BackConstrain:=backConstraint;
Result := Distance+penaltyW*weightconstraint+penaltyB*backConstraint;
end;</P>
<P>procedure TTSPExaminer.SetController(const Value: ITSPController);
begin
fController := Value;
end;</P>
<P>{ TPopulation }</P>
<P>constructor TPopulation.Create;
begin
inherited;
fPop := TInterfaceList.Create;
end;</P>
<P>destructor TPopulation.Destroy;
begin
fPop.Free;
inherited;
end;</P>
<P>procedure TPopulation.Add(New: IIndividual);
begin
fPop.Add(New);
end;</P>
<P>procedure TPopulation.Clear;
begin
fPop.Clear;
end;</P>
<P>function TPopulation.CompareIndividuals(I1, I2: IIndividual): Integer;
var
A, B, D : TFloat;
begin
// Get the difference between the two individuals (real number)
A := I1.Fitness;
B := I2.Fitness;</P>
<P>D := A - B;</P>
<P>// Quickest way to convert that to an integer is...
if D > 0 then
Result := 1
else if D < 0 then
Result := -1
else
Result := 0;
end;</P>
<P>procedure TPopulation.Delete(I: Integer);
begin
fPop.Delete(I);
end;</P>
<P>function TPopulation.FitnessOf(I: Integer): TFloat;
begin
result := Pop.Fitness;
end;</P>
<P>procedure TPopulation.Change;
begin
if Assigned(fOnChange) then
FOnChange(Self);
end;</P>
<P>procedure TPopulation.Generation;
var
Replace, i : Integer;
New : IIndividual;
begin
// Kill some of the population
Replace := fKiller.Kill(Self);</P>
<P>for i := 1 to Replace do
begin
// Breed a new individual
New := fBreeder.BreedOffspring(fParentSelector, Self);
// Perform some mutation on the individual
FMutator.Mutate(New);
// Get the fitness of the new individual
New.Fitness := fExaminer.GetFitness(New);
// Add it to the population
Add(New);
end;
// Sort the population into fitness order where first <==> best
SortPopulation;</P>
<P>Change;
end;</P>
<P>function TPopulation.GetBreeder: IBreeder;
begin
result := fBreeder;
end;</P>
<P>function TPopulation.GetCount: Integer;
begin
result := fPop.Count;
end;</P>
<P>function TPopulation.GetCreator: ICreator;
begin
result := fCreator;
end;</P>
<P>function TPopulation.GetExaminer: IExaminer;
begin
result := fExaminer;
end;</P>
<P>function TPopulation.GetIndividual(I: Integer): IIndividual;
begin
result := (fPop as IIndividual);
end;</P>
<P>function TPopulation.GetKiller: IKiller;
begin
result := fKiller;
end;</P>
<P>function TPopulation.GetMutator: IMutator;
begin
result := fMutator;
end;</P>
<P>function TPopulation.GetOnChange: TNotifyEvent;
begin
result := fOnChange;
end;</P>
<P>function TPopulation.GetParentSelector: IParentSelector;
begin
result := fParentSelector;
end;</P>
<P>procedure TPopulation.Initialise(Size: Integer);
var
i,feasibleCount: Integer;
New: IIndividual;
begin
feasibleCount:=round(size*(FormGaPara.EditFeasible.Value)/100);
//feasibleCount:=1;
// Clear out the old stuff
Clear;
// Set the capacity first to save about 12 nanoseconds ;o)
fPop.Capacity := Size;
// Create the appropriate number of individuals
for i := 1 to feasibleCount do
begin
// Create the individual
New := fCreator.CreateFeasibleIndividual;
// Get the fitness of the new individual
New.Fitness := fExaminer.GetFitness(New);
// Add to the population
Add(New);
end;
for i := feasibleCount+1 to Size do
begin
// Create the individual
New := fCreator.CreateIndividual; ////////
// Get the fitness of the new individual
New.Fitness := fExaminer.GetFitness(New);
// Add to the population
Add(New);
end;
SortPopulation;
Change;
end;</P>
<P>procedure TPopulation.SetBreeder(const Value: IBreeder);
begin
fBreeder := Value;
end;</P>
<P>procedure TPopulation.SetCreator(const Value: ICreator);
begin
fCreator := Value;
end;</P>
<P>procedure TPopulation.SetExaminer(const Value: IExaminer);
begin
fExaminer := Value;
end;</P>
<P>procedure TPopulation.SetKiller(const Value: IKiller);
begin
fKiller := Value;
end;</P>
<P>procedure TPopulation.SetMutator(const Value: IMutator);
begin
fMutator := Value;
end;</P>
<P>procedure TPopulation.SetOnChange(const Value: TNotifyEvent);
begin
fOnChange := Value;
end;</P>
<P>procedure TPopulation.SetParentSelector(const Value: IParentSelector);
begin
fParentSelector := Value;
end;</P>
<P>procedure TPopulation.DanQuickSort(SortList: TInterfaceList; L, R: Integer;
SCompare: TInterfaceCompare);
var
I, J: Integer;
P: IIndividual;
begin
repeat
I := L;
J := R;
P := SortList.Items[(L + R) div 2] as IIndividual;
repeat
while SCompare(SortList.Items as IIndividual, P) < 0 do
Inc(I);
while SCompare(SortList.Items as IIndividual, P) > 0 do
Dec(J);
if I <= J then
begin
SortList.Exchange(I, J);
Inc(I);
Dec(J);
end;
until I > J;
if L < J then
DanQuickSort(SortList, L, J, SCompare);
L := I;
until I >= R;
end;</P>
<P>procedure TPopulation.Sort(Compare: TInterfaceCompare);
begin
if Assigned(fPop) and (Count > 0) then
DanQuickSort(fPop, 0, Count - 1, Compare);
end;</P>
<P>procedure TPopulation.SortPopulation;
begin
Sort(self.CompareIndividuals);
end;</P>
<P>{ TTSPController }</P>
<P>constructor TTSPController.Create;
begin
inherited;
end;</P>
<P>destructor TTSPController.Destroy;
begin
SetLength(FCities, 0);
SetLength(FNoVehicles, 0);
SetLength(FVehicles, 0);
inherited;
end;</P>
<P>{ Standard euclidian distance between two 2D vectors... }
function TTSPController.DistanceBetween( C1, C2: Integer): TFloat;
var
temp:TFloat;
i,j,intTemp,intTemp2:TInt;
begin
intTemp:=0;
temp:=FormGaPara.EditWrongConstrain.Value;</P>
<P>{if (Cities.id>=1)and(Cities.id<=fOldCityCount)and(Cities.id<>0) then
begin
fCities.serviceDepot:= fCities.serviceDepot;
end; //}
//8
if (Cities.id>=1)and(Cities.id<=fOldCityCount)and(Cities.id>=1)and(Cities.id<=fOldCityCount) then
begin
temp:=CostArray;
end;
//1
if Cities.id=0 then
begin
for i:=0 to fDepotCount-1 do
begin
intTemp:=intTemp+fNoVehicles;
if Cities.id =fOldCityCount + intTemp +1 then
temp:=0;
end;
intTemp:=0;
end;
//2
if Cities.id=0 then
begin
for i:=1 to fDepotCount do
begin
intTemp:=intTemp+fNoVehicles;
if Cities.id =fOldCityCount + intTemp then
temp:=0;
end;
intTemp:=0;
end;
//5
for i:=0 to fDepotCount-1 do
begin
intTemp:=intTemp+fNoVehicles;
{ if (Cities.id=fOldCityCount + intTemp +1)and(Cities.id=Cities.id+1) then
temp:=10; /////////////////////////// }
if (Cities.id>=fOldCityCount + intTemp +1)and(Cities.id<=fOldCityCount + intTemp+fNoVehicles)
and(Cities.id>=fOldCityCount + intTemp +1)and(Cities.id<=fOldCityCount + intTemp+fNoVehicles)
then
temp:=0;//}
end;
intTemp:=0;
//7
if (Cities.id=Cities.id)and(Cities.id > fOldCityCount) then
begin
temp:=0;
end;
//3
if (Cities.id > fOldCityCount)and(Cities.id>=1)and(Cities.id<=fOldCityCount) then
begin
//temp := Sqrt(sqr(Cities.X - Cities.X)+sqr(Cities.Y - Cities.Y));
temp:=CostArray;
end;
//4
if (Cities.id<=fOldCityCount)and(Cities.id>=1)and(Cities.id > fOldCityCount) then
begin
//if Cities.serviceDepot=Cities.serviceDepot then //back to the start point
temp:=CostArray;
end;
//6
intTemp:=0;
for i:=1 to fDepotCount do
begin
intTemp:=intTemp+fNoVehicles;
if Cities.id= fOldCityCount + intTemp then
begin
intTemp2:=0;
for j:=0 to fDepotCount-1 do
begin
intTemp2:=intTemp2+fNoVehicles;
if Cities.id=fOldCityCount + intTemp2 +1 then
if abs(Cities.id-Cities.id) <> fNoVehicles-1 then
temp:=0;
end; //}</P>
<P>end;
end;
intTemp:=0;
result:=temp;
end;</P>
<P>function TTSPController.CostBetween(C1, C2: Integer): TFloat; //matrix cij
var
distance:TFloat;
begin
distance := Sqrt(sqr(Cities.X - Cities.X)+ sqr(Cities.Y - Cities.Y));
//result:=distance+TimeCostBetween(C1,C2);
result:=distance;
end;</P>
<P>function TTSPController.TimeCostBetween(C1, C2: Integer): TFloat;
var
cost:TFloat;
i,j,penaltyW,penaltyD:TFloat;
startTime:TDateTime;
begin
startTime:=strToDateTime(FormGa.EditStartTime.Text);
penaltyW:=FormGaPara.EditWaitConstrain.Value;
penaltyD:=FormGaPara.EditDelayConstrain.Value;
if Cities.id>fOldCityCount then
fCities.totalTime:=0
else
fCities.totalTime:=Cities.totalTime+Cities.serviceTime+timeArray;</P>
<P>fCities.waitTime:= max(0,DateSpanToMin(startTime,Cities.early)-Cities.totalTime);
fCities.delayTime:=max(0,Cities.totalTime-DateSpanToMin(startTime,Cities.late));</P>
<P>if Cities.late<>0 then //consider time or not
begin
if Cities.early<>0 then //window or deadline
cost:=penaltyW*fCities.waitTime +penaltyD*fCities.delayTime
else
cost:=penaltyD*fCities.delayTime;
end
else
cost:=0;
result:=cost;
end;</P>
<P>function TTSPController.DateSpanToMin(d1,d2:TDateTime):integer;
var
span:TDateTime;
Year, Month, Day, Hour, Min, Sec, MSec: Word;
begin
span:=abs(d2-d1);
DecodeDate(span, Year, Month, Day);
DecodeTime(span, Hour, Min, Sec, MSec);
result:=Min;
end;</P>
<P>//return the position in the vehicles array
function TTSPController.GetVehicleInfo( routeInt:Tint):integer;
begin
result:=routeInt-fOldCityCount-1;
end;</P>
<P>function TTSPController.GetWeightConstraint( Individual: IIndividual): TInt;
var
Indi: ITSPIndividual;
totalCapacity,maxCapacity: TFloat;
i,j:TInt;
tempArray:array of TInt;
tempResult:TInt;
begin
Indi := Individual as ITSPIndividual;
SetLength(tempArray, fCityCount+1);
tempResult:=0;
/////////////////////////////////////////////////////////
for i:=0 to fCityCount-1 do
begin
if Indi.RouteArray=fOldCityCount+1 then
break;
end;
for j:=0 to fCityCount-i-1 do
begin
tempArray:= Indi.RouteArray;
end;
for j:=fCityCount-i to fCityCount-1 do
begin
tempArray:= Indi.RouteArray;
end;
tempArray:= tempArray;
//////////////////////////////////////////////////////////
//totalCapacity:=fCities].supply; //supply
maxCapacity:=fVehicles)].volume;
totalCapacity:=maxCapacity;
for i:=0 to fCityCount do
begin
if (FCities].id<=fOldCityCount)and(FCities].id>0) then
begin
totalCapacity:=totalCapacity+FCities].supply-FCities].demand;
if (totalCapacity>maxCapacity)or(totalCapacity<0) then
begin
tempResult:=tempResult+1;
//break;
end;
end;
if FCities].id>fOldCityCount then
begin
//totalCapacity:=fCities].supply; //supply
maxCapacity:=fVehicles)].volume;
totalCapacity:=maxCapacity;
end;
end;
SetLength(tempArray,0);
result:=tempResult;
end;</P>
<P>function TTSPController.GetBackConstraint( Individual: IIndividual): TInt;
var
Indi: ITSPIndividual;
i,j:TInt;
tempArray:array of TInt;
tempResult:TInt;
begin
Indi := Individual as ITSPIndividual;
SetLength(tempArray, fCityCount+1);
tempResult:=0;
for i:=0 to fCityCount-1 do
begin
if Indi.RouteArray=fOldCityCount+1 then
break;
end;
for j:=0 to fCityCount-i-1 do
begin
tempArray:= Indi.RouteArray;
end;
for j:=fCityCount-i to fCityCount-1 do
begin
tempArray:= Indi.RouteArray;
end;
tempArray:=tempArray;
{tempArray:=11;tempArray:=5;tempArray:=8;tempArray:=7;
tempArray:=9;tempArray:=6;tempArray:=12;tempArray:=10;
tempArray:=2;tempArray:=4;tempArray:=3;tempArray:=1;
tempArray:=0;tempArray:=11;tempArray:=3;tempArray:=1;
tempArray:=4;tempArray:=11;//10,2,2}
for i:=0 to fCityCount-1 do
begin
if (Cities].id<=fOldCityCount) then
begin
fCities].serviceDepot:= fCities].serviceDepot;
end;
if (Cities].id<=fOldCityCount)and(Cities].id>=1)and(Cities].id > fOldCityCount) then
begin
if Cities].serviceDepot<>Cities].serviceDepot then //back to the start point
begin
tempResult:=tempResult+1;
// break;
end;
end;
end;
SetLength(tempArray,0);
result:=tempResult;
end; </P>
<P>function TTSPController.GetTimeConstraint( Individual: IIndividual): TInt;
var
Indi: ITSPIndividual;
i,j:TInt;
totalTimeCost:TFloat;
tempArray:array of TInt;
tempResult:TInt;
begin
Indi := Individual as ITSPIndividual;
SetLength(tempArray, fCityCount+1);
tempResult:=0;
for i:=0 to fCityCount-1 do
begin
if Indi.RouteArray=fOldCityCount+1 then
break;
end;
for j:=0 to fCityCount-i-1 do
begin
tempArray:= Indi.RouteArray;
end;
for j:=fCityCount-i to fCityCount-1 do
begin
tempArray:= Indi.RouteArray;
end;
tempArray:=tempArray;</P>
<P>totalTimeCost:=0;
for i:=0 to fCityCount-1 do
begin
totalTimeCost:=totalTimeCost+timeCostBetween(tempArray,tempArray);
end;
if totalTimeCost<>0 then tempResult:=1;
SetLength(tempArray,0);
end;</P>
<P>function TTSPController.GetCity(I: Integer): TPoint2D;
begin
result := fCities;
end;</P>
<P>function TTSPController.GetNoVehicle(I: Integer): TInt;
begin
result := fNoVehicles;
end;</P>
<P>function TTSPController.GetCityCount: Integer;
begin
result := fCityCount;
end;</P>
<P>function TTSPController.GetOldCityCount: Integer;
begin
result := fOldCityCount;
end;</P>
<P>function TTSPController.GetTravelCount: Integer;
begin
result := fTravelCount;
end;</P>
<P>function TTSPController.GetDepotCount: Integer;
begin
result := fDepotCount;
end;</P>
<P>function TTSPController.GetXmax: TFloat;
begin
result := fXmax;
end;</P>
<P>function TTSPController.GetXmin: TFloat;
begin
result := fXmin;
end;</P>
<P>function TTSPController.GetYmax: TFloat;
begin
result := fYmax;
end;</P>
<P>function TTSPController.GetYmin: TFloat;
begin
result := fYmin;
end;</P>
<P>procedure TTSPController.RandomCities; //from database
var
i,j,k,m,intTemp,totalVehicleCount: Integer;
tempVehicle:TVehicle;
begin
//////////////////////////////////////////////////////////
fNoVehicles:=0;
totalVehicleCount:=0;
for i:=1 to fDepotCount do //from depots database
begin
fNoVehicles:=fTravelCount +1;
totalVehicleCount:=totalVehicleCount+ fNoVehicles; //real and virtual vehicles
end;
SetLength(fVehicles,totalVehicleCount);
intTemp:=0;
for i:=1 to fDepotCount do
begin
for j:=intTemp to intTemp+fNoVehicles-2 do
begin
fVehicles.index:=j+1;
fVehicles.id:='real vehicle';
fVehicles.volume:=50;
end;
with fVehicles-1] do
begin
index:=intTemp+fNoVehicles;
id:='virtual vehicle';
volume:=0;
end;
intTemp:=intTemp+ fNoVehicles;
end;</P>
<P>///////////////////////////////////////////////////////////
intTemp:=0;
for i:=1 to fDepotCount do //depot 1--value
begin
intTemp:=intTemp + fNoVehicles;
end;</P>
<P>for i := 0 to FOldCityCount do //from database
begin
FCities.id:= i;
FCities.X := Xmin+0.5+(Xmax-Xmin-1.0)*Random;
FCities.Y := Ymin+0.5+(Ymax-Ymin-1.0)*Random;
FCities.early:=0;
FCities.late:=0; //TDateTime
FCities.serviceTime:=0;
FCities.totalTime:=0;
FCities.waitTime:=0;
FCities.delayTime:=0;
end;
for i:=FOldCityCount+1 to FCityCount-1 do
begin
FCities.id:= i;
if fDepotCount=1 then
begin
FCities.X := Xmin+0.5+(Xmax-Xmin-1.0)*RandomRange(2,4)/5;
FCities.Y := Ymin+0.5+(Ymax-Ymin-1.0)*RandomRange(2,4)/5;
end
else
begin
FCities.X := Xmin+0.5+(Xmax-Xmin-1.0)*Random;
FCities.Y := Ymin+0.5+(Ymax-Ymin-1.0)*Random;
end;
FCities.early:=0;
FCities.late:=0; //TDateTime
FCities.serviceTime:=0;
FCities.totalTime:=0;
FCities.waitTime:=0;
FCities.delayTime:=0;
end;</P>
<P>for i := 0 to FOldCityCount do
begin
FCities.serviceDepot:=i;
end;</P>
<P>m:=FOldCityCount+1;
for k:=1 to fDepotCount do
begin
for j:=0 to fNoVehicles-1 do
begin
FCities.serviceDepot:= fOldCityCount+k;
m:=m+1;
end;
end;</P>
<P>//supply and demand //////////////////////////from database
FCities.demand:=0;
FCities.supply:=0;
for i:=1 to FOldCityCount do
begin
FCities.demand:=10;
FCities.supply:=0;
end;
for i:=FOldCityCount+1 to FCityCount-1 do
begin
FCities.demand:=0;
FCities.supply:=50;
end;
////////////////////////////////////////////////////////////</P>
<P>intTemp:=0;
for i:=0 to fDepotCount-1 do
begin
intTemp:=intTemp+fNoVehicles;
for j:=2 to fNoVehicles do
begin
FCities.X :=FCities.X;
FCities.Y :=FCities.Y;
end;
end;
writeTimeArray;
writeCostArray;
end;</P>
<P>procedure TTSPController.writeTimeArray; //database
var
i,j:integer;
begin
SetLength(timeArray,fCityCount,fCityCount);
for i:=0 to fCityCount-1 do
begin
for j:=0 to fCityCount-1 do
begin
if i=j then timeArray:=0
else timeArray:=10;
end;
end;
end;</P>
<P>procedure TTSPController.writeCostArray; //database
var
i,j:integer;
begin
SetLength(costArray,fCityCount,fCityCount);
for i:=0 to fCityCount-1 do
begin
for j:=0 to fCityCount-1 do
begin
if i=j then costArray:=0
else costArray:=costBetween(i,j);
end;
end;
end;</P>
<P>procedure TTSPController.SetCityCount(const Value: Integer);
begin
SetLength(fCities, Value);
fCityCount := Value;</P>
<P>RandomCities;
end;</P>
<P>procedure TTSPController.SetOldCityCount(const Value: Integer);
begin
fOldCityCount := Value;
end;</P>
<P>procedure TTSPController.SetTravelCount(const Value: Integer); ///////////
begin
fTravelCount := Value;
end;</P>
<P>procedure TTSPController.SetDepotCount(const Value: Integer); ///////////
begin
SetLength(fNoVehicles, Value+1); ///////////////
fDepotCount := Value;
end;</P>
<P>procedure TTSPController.SetXmax(const Value: TFloat);
begin
fXmax := Value;
end;</P>
<P>procedure TTSPController.SetXmin(const Value: TFloat);
begin
fXmin := Value;
end;</P>
<P>procedure TTSPController.SetYmax(const Value: TFloat);
begin
fYmax := Value;
end;</P>
<P>procedure TTSPController.SetYmin(const Value: TFloat);
begin
fYmin := Value;
end;</P>
<P>end.
</P></DIV>
[此贴子已经被作者于2005-4-27 15:51:02编辑过]
进化算法
<P>Evocosm encapsulates the fundamental principles of evolutionary algorithms in an extensible set of standard C++ templates and tools. Evolutionary algorithms come in a variety of shapes and flavors -- genetic algorithms, genetic programming, evolutionary computing -- but at their core, they all share certain characteristics: populations that reproduce and mutate through a series of generations, producing future generations based on some measure of fitness.Evocosm将进化算法的基本原则封装在一套可扩展的标准C++模板和工具中。进化算法可以各种各样--遗传算法、遗传编程、进化计算等等，但是它们的核心都共享一些特点：种群通过几代来繁殖和成熟，基于某种适当性量度来产生未来的代。</P>
<P>这是进化算法解决TSP问题的源代码：</P>
蚂蚁算法
<P>该算法是由意大利学者M.Dorigo、V.Maniez-zo、A.Colorini等人首先提出的,称之为蚁群系统(antcolonysystem), 该模型已成功应用于求旅行商问题(ＴＳＰ),二次指派问题,排序问题等ＮＰ-困难的组合最优化问题,结果可与模拟退火,遗传算法等通用的启发式算法相媲美.蚁群算法和局部搜索算法相结合(称为混合蚁群算法)应用于解二次指派问题和排序问题,得到的结果可以与专用算法相媲美].受其影响,蚁群系统模型逐渐引起了其它研究者的注意,D.Costa和A.Hertz.</P><P>在M.Dorigo等人研究成果的基础上,提出了一种求解分配类型问题(assignmenttypeproblem)的一般模型,并用来研究图着色问题.G.Bilchev、I.C.Parmee研究了求解连续空间优化问题的蚁群系统模型.。</P>
<P>蚁群算法是模仿蚂蚁工作方式的一种新的启发式算法.生物学研究表明一群互相协作的蚂蚁能够找到食物源和巢之间的最短路径,而单只蚂蚁则不能.蚂蚁间相互协作的方法是它们在所经过的路上留下一定数量的信息素(迹),该迹能被其它蚂蚁检测出来,一条路上的迹越多,其它蚂蚁将以越高的概率跟随此路径,从而该路径上的迹会被加强.</P>
<P>蚂蚁算法伪码
Begin
初始化:
t← 0
iteration← 0 (iteration为迭代步数)
将 m个蚂蚁随机置于 n个顶点上 ;
Loop:
将所有蚂蚁的初始出发点置于当前解集中 ;
for i← 0 to n-1 do
for k← 1 to m do
按概率 选择顶点 j;
移动蚂蚁 k至顶点 j;
将顶点j置于蚂蚁k的当前解集</P>
<P>end for
t←t＋1
end for
计算各蚂蚁的 L个目标函数值
更新当前的理想解
计算各个解的满意度 </P>
<P>置 t← t+1
重置所有 ← 0
iteration← iteration+1
若 iteration<预定的迭代次数
则 goto Loop
输出目前的最满意解
End</P>
<P>下面是蚂蚁算法解决TSP问题的C++程序：</P>
<DIV class=HtmlCode>
<P>/* >->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->
AntColonySystemTSP.cc, Heiko Stamer </P>
<P>Ant Colony System (ACS) for the Traveling Salesman Problem (TSP) </P>
<P>
by M. Dorigo, V. Maniezzo, A. Colorni
IEEE Transactions on Systems, Man and Cybernetics - Part B, Vol.26-1 1996 </P>
<P>
by M. Dorigo and L. M. Gambardella
IEEE Transactions on Evolutionary Computation, Vol. 1, No. 1, 1997 </P>
<P><a href="http://stinfwww.informatik.uni-leipzig.de/~mai97ixb" target="_blank" >http://stinfwww.informatik.uni-leipzig.de/~mai97ixb</A>
>->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->->-> </P>
<P>Copyright (C) 2001 - until_the_end_of_the_ants </P>
<P>This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version. </P>
<P>This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details. </P>
<P>You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ </P>
<P>#include
#include
#include
#include
#include
#include </P>
<P>#define N 70 </P>
<P>double C = { {64 , 96} , {80 , 39} , {69 , 23} , {72 , 42} , {48 , 67} , {58 ,
43} , {81 , 34} , {79 , 17} , {30 , 23} , {42 , 67} , {7 , 76} , {29 , 51} , {78
, 92} , {64 , 8} , {95 , 57} , {57 , 91} , {40 , 35} , {68 , 40} , {92 , 34} ,
{62 , 1} , {28 , 43} , {76 , 73} , {67 , 88} , {93 , 54} , {6 , 8} , {87 , 18} ,
{30 , 9} , {77 , 13} , {78 , 94} , {55 , 3} , {82 , 88} , {73 , 28} , {20 , 55}
, {27 , 43} , {95 , 86} , {67 , 99} , {48 , 83} , {75 , 81} , {8 , 19} , {20 ,
18} , {54 , 38} , {63 , 36} , {44 , 33} , {52 , 18} , {12 , 13} , {25 , 5} , {58
, 85} , {5 , 67} , {90 , 9} , {41 , 76} , {25 , 76} , {37 , 64} , {56 , 63} ,
{10 , 55} , {98 , 7} , {16 , 74} , {89 , 60} , {48 , 82} , {81 , 76} , {29 , 60}
, {17 , 22} , {5 , 45} , {79 , 70} , {9 , 100} , {17 , 82} , {74 , 67} , {10 ,
68} , {48 , 19} , {83 , 86} , {84 , 94} }; </P>
<P>typedef int Tour;
typedef double doubleMatrix; </P>
<P>doubleMatrix D; </P>
<P>double dist(int i, int j)
{
return sqrt(pow((C-C), 2.0) + pow((C-C), 2.0));
} </P>
<P>void calc_dist()
{
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
D = dist(i, j);
} </P>
<P>double max_dist()
{
double max_dist = 0.0;
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
if (dist(i, j) > max_dist)
max_dist = dist(i, j);
return max_dist;
} </P>
<P>double calc_length(Tour tour)
{
double l = 0.0;
for (int n = 0; n < N; n++)
{
int i = tour;
int j = tour;
l += D;
}
return (l);
} </P>
<P>void print_tour(Tour tour)
{
for (int n = 0; n < N; n++)
printf("( %d , %d ) ", tour, tour);
printf("\n");
} </P>
<P>int sum_sequence(int array[], int count)
{
int sum = 0;
for (int i = 0; i < count; i++)
sum += array;
return (sum);
} </P>
<P>/******************************************************************************/ </P>
<P>class Ant
{
protected:
int START_CITY, CURRENT_CITY;
int ALLOWED;
Tour CURRENT_TOUR;
int CURRENT_TOUR_INDEX;
public:
inline Ant(int start_city)
{
START_CITY = start_city;
}
inline void moveTo(int to_city)
{
ALLOWED = 0;
CURRENT_TOUR = CURRENT_CITY;
CURRENT_TOUR = to_city;
CURRENT_TOUR_INDEX++;
CURRENT_CITY = to_city;
}
}; </P>
<P>class NNAnt : Ant
{
public:
inline NNAnt(int start_city): Ant(start_city) { };
inline int choose()
{
double best_length = (double)N * max_dist();
int best_choose = -1; </P>
<P>for (int j = 0; j < N; j++)
{
if ((ALLOWED == 1) && (D < best_length))
{
best_choose = j;
best_length = D;
}
}
return best_choose;
}
inline Tour *search()
{
CURRENT_CITY = START_CITY;
CURRENT_TOUR_INDEX = 0;
for (int i = 0; i < N; i++)
ALLOWED = 1;
ALLOWED = 0;
while (sum_sequence(ALLOWED, N) > 0)
moveTo(choose());
ALLOWED = 1;
moveTo(START_CITY);
return &CURRENT_TOUR;
}
}; </P>
<P>class AntColonySystem;
class ACSAnt : Ant
{
private:
AntColonySystem *ACS;
public:
ACSAnt(AntColonySystem *acs, int start_city): Ant(start_city)
{
ACS = acs;
}
inline int choose();
inline Tour *search();
}; </P>
<P>class AntColonySystem
{
private:
double ALPHA, BETA, RHO, TAU0;
doubleMatrix TAU, dTAU;
static const int M = 420;
ACSAnt *ANTS; </P>
<P>public:
double Q0;
AntColonySystem(double alpha, double beta, double rho, double q0);
inline double calc_tau0();
inline void init_tau_by_value(double value);
inline void init_tau_by_matrix(doubleMatrix matrix);
inline void init_uniform();
inline void init_random();
inline void init_randomMOAPC();
inline double ETA(int i, int j);
inline double transition(int i, int j);
inline double sum_transition(int i, int allowed[]);
inline void local_update_rule(int i, int j);
inline void clear_global_update();
inline void add_global_update(Tour tour, double length);
inline void global_update_rule();
inline doubleMatrix *get_tau();
inline Tour *search(int T);
}; </P>
<P>inline int ACSAnt::choose()
{
double q = rand() / (double)RAND_MAX; </P>
<P>if (q <= ACS->Q0)
{
double best_value = -1.0;
int best_choose = -1;
for (int j = 0; j < N; j++)
{
if ((ALLOWED == 1) &&
(ACS->transition(CURRENT_CITY, j) > best_value))
{
best_choose = j;
best_value = ACS->transition(CURRENT_CITY, j);
}
}
return best_choose;
} </P>
<P>double sum = ACS->sum_transition(CURRENT_CITY, ALLOWED);
double p = rand() / (double)RAND_MAX;
double p_j = 0.0; </P>
<P>for (int j = 0; j < N; j++)
{
if (ALLOWED == 1) p_j += ACS->transition(CURRENT_CITY, j) / sum;
if ((p < p_j) && (ALLOWED == 1))
return j;
}
return -1;
} </P>
<P>inline Tour *ACSAnt::search()
{
CURRENT_CITY = START_CITY;
CURRENT_TOUR_INDEX = 0;
for (int i = 0; i < N; i++)
ALLOWED = 1;
ALLOWED = 0;
while (sum_sequence(ALLOWED, N) > 0)
{
int LAST_CITY = CURRENT_CITY;
moveTo(choose());
ACS->local_update_rule(LAST_CITY, CURRENT_CITY);
}
ALLOWED = 1;
ACS->local_update_rule(CURRENT_CITY, START_CITY);
moveTo(START_CITY);
return &CURRENT_TOUR;
} </P>
<P>/******************************************************************************/ </P>
<P>AntColonySystem::AntColonySystem(double alpha, double beta, double rho, double q0)
{
ALPHA = alpha;
BETA = beta;
RHO = rho;
Q0 = q0;
} </P>
<P>inline double AntColonySystem::calc_tau0()
{
double best_length = (double)N * max_dist(); </P>
<P>for (int n = 0; n < N; n++)
{
NNAnt *nnANT = new NNAnt(n);
Tour tour;
tour = *(nnANT->search());
double tour_length = calc_length(tour);
if (tour_length < best_length)
best_length = tour_length;
delete nnANT;
}
return 1.0 / ((double)N * best_length);
} </P>
<P>inline void AntColonySystem::init_tau_by_value(double value)
{
TAU0 = value;
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
TAU = TAU0;
} </P>
<P>inline void AntColonySystem::init_tau_by_matrix(doubleMatrix matrix)
{
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
TAU = matrix;
} </P>
<P>inline void AntColonySystem::init_uniform()
{
// uniformly distributed
for (int k = 0; k < M; k++)
ANTS = new ACSAnt(this, (k % N));
} </P>
<P>inline void AntColonySystem::init_random()
{
// randomly distributed
for (int k = 0; k < M; k++)
ANTS = new ACSAnt(this, (int)((double)N * (rand() / (double)RAND_MAX)));
} </P>
<P>inline void AntColonySystem::init_randomMOAPC()
{
// randomly distributed with MOAPC (most one ant per city)
bool MOAPCarray;
assert(M <= N); </P>
<P>for (int n = 0; n < N; n++)
MOAPCarray = false; </P>
<P>for (int k = 0; k < M; k++)
{
int c;
do
{
c = (int)((double)N * (rand() / (double)RAND_MAX));
}
while (MOAPCarray); </P>
<P>MOAPCarray = true;
ANTS = new ACSAnt(this, c);
}
} </P>
<P>inline double AntColonySystem::ETA(int i, int j)
{
return ( 1.0 / D );
} </P>
<P>inline double AntColonySystem::transition(int i, int j)
{
if (i != j)
return ( TAU * pow( ETA(i, j), BETA ) );
else
return(0.0);
} </P>
<P>inline double AntColonySystem::sum_transition(int i, int allowed[])
{
double sum = 0.0;
for (int j = 0; j < N; j++)
sum += ((double)allowed * transition(i, j));
return (sum);
} </P>
<P>inline void AntColonySystem::local_update_rule(int i, int j)
{
TAU = (1.0 - RHO) * TAU + RHO * TAU0;
// symmetric TSP
TAU = TAU;
} </P>
<P>inline void AntColonySystem::clear_global_update()
{
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
dTAU = 0.0;
} </P>
<P>
inline void AntColonySystem::add_global_update(Tour tour, double length)
{
for (int n = 0; n < N; n++)
{
int i = tour;
int j = tour;
dTAU += (1.0 / length);
// symmetric TSP
dTAU += (1.0 / length);
}
} </P>
<P>inline void AntColonySystem::global_update_rule()
{
for (int i = 0; i < N; i++)
for (int j = 0; j < N; j++)
TAU = (1.0 - ALPHA) * TAU + ALPHA * dTAU;
} </P>
<P>inline doubleMatrix *AntColonySystem::get_tau()
{
return &TAU;
} </P>
<P>inline Tour *AntColonySystem::search(int T)
{
Tour best_tour, tour;
double best_length = (double)N * max_dist(), tour_length;
clear_global_update(); </P>
<P>// do T iterations of ACS algorithm
int t;
for (t = 0; t < T; t++)
{
for (int k = 0; k < M; k++)
{
tour = *(ANTS->search());
tour_length = calc_length(tour);
if (tour_length < best_length)
{
best_tour = tour;
best_length = tour_length;
clear_global_update();
add_global_update(tour, tour_length);
//printf("[%d / %d]: %lf \n", t, T, tour_length);
}
}
global_update_rule();
} </P>
<P>//printf("[%d/%d] best tour (length = %f):\n", t, T, best_length);
//print_tour(best_tour);
//printf("[%d/%d] iterations done\n", t, T);
printf("%f\n", best_length);
return (&best_tour);
} </P>
<P>int main(int argc, char* argv[])
{
// PRNG initalisieren
time_t timer;
time(&timer);
pid_t pid = getpid() + getppid();
unsigned long seed = (timer * pid);
if (seed == 0)
{
time(&timer);
seed = 7 * timer * pid;
if (seed == 0) seed = pid; else seed = seed % 56000;
} else seed = seed % 56000;
srand((unsigned int)seed); </P>
<P>// EUC2D
calc_dist(); </P>
<P>// Ant Colony System
AntColonySystem *acs = new AntColonySystem(0.1, 2.0, 0.1, 0.9);
double tau0 = acs->calc_tau0();
acs->init_tau_by_value(tau0);
acs->init_uniform();
acs->search(1000); </P>
<P>return(0);
} </P></DIV>
<P>蚂蚁算法的一些文献：</P>
<P>
</P> 感谢共享!!,不过这样传上来看着不是很方便.不过还是非常感谢共享你的资源!! 多谢楼主,对我帮助很大. 谢谢楼主！！！ <font size="5">谢谢楼主！！！谢谢楼主！！！</font> <div>好</div> <a href="http://www.madio.net/bbs/dispbbs.asp?BoardID=107&ID=4088&replyID=32632&skin=1"><font color="#000000" size="+0">多谢楼主,对我帮助很大.</font></a>