蚁群聚类算法matlab源代码
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
1.function Q = FINDabc(A) %给定 B1,B2 两点坐标计算平面参数 a,b,c v = 1e-3*[1 1 1 1];r=1.3;l=0; K=6; a = RAND; a = RESHAP_MAX(A,a,K); while(sum(STAND(A,a(K,:))>v)~=0) l=0; [Q,w]=REF(a,r,K); while((~D_ZONE1(Q))|(COMP(A,Q)>=COMP(A,a(K,:)))) [Q,w]=REF(a,w,K); if (w<=eps^12) a(K,:)=[]; K=K-1; l = 1; end if(l==1) break; end end if(~l) a(K,:)=Q; end if(K<=1) a=FINDabc(A); break; end a = RESHAP_MAX(A,a,K); end if(K<=1) Q = a(1,:); else Q=sum(a)./K; end 11 2.function P=TEXT(i) %对 i 组[a,b,c]检验,返回[a,b,c,a',b',c'] j = 1; while(j<=i) A(1,1:2)= CREAT_RAND12; A(1,3)=CREAT_RAND3; if(D_ZONE(A(1,:))==1) b(j,:)=A(1,:);j=j+1; end
if x(i)==j,p=i; end end a=A(p,:); A(p,:) = A(K,:); A(K,:)=a; Q_M=A; End
12.function w=STAND(A,R) D=-7*1e-1;h=3*1e-1; B = -.253; F = -.259; a = R(1); b = R(2); c = R(3); B11 = A(1); B12 = A(2); B21 = A(3); B22 = A(4); s3 = 2*(c+D)/(1+a^2+b^2); s2=-b*s3; s1=-a*s3; c13 = (c+D)/(1-b*B); c12 = B*c13; c23 = (c+D)/(1-b*F); c22 = F*c23; F11 = (h-s3)*(-s1)/(c13-s3)+s1; F12 = (h-s3)*(c12-s2)/(c13-s3)+s2; F21 = (h-s3)*(-s1)/(c23-s3)+s1; F22 = (h-s3)*(c22-s2)/(c23-s3)+s2; w(1)=abs(B11-F11); w(2)=abs(B12-F12); w(3)=abs(B21-F21); w(4)=abs(B22-F22);
3.function p = recognize(i) %分辨率的影响 j = 1; while(j<=i) a(1,1:2) = CREAT_RAND12 (); a(1,3) = CREAT_RAND3(); if(D_ZONE(a(1,:))) b(j,:) = a(1,:); j = j+1; end end for j = 1:i c(j,:) = Fx(b(j,:)); e(j,:)=1e-4*round(1e4*c(j,:)); end for j = 1:i f(j,1:3) = FINDabc(e(j,:)); end hl = (b-f)./b; for j = 1:3 h(:,j)=sum(abs(hl(:,j)))/i*ones(i,1); end plot(h(:,1),'--r');hold on;plot(h(:,2),'-y');plot(h(:,3),':k'); p = h; 12 4.function Q_C12 =CREAT_RAND12 () if(y_C12<=0.5) Q_C12(1) =0.0017 *rand(1); else Q_C12(1) = -0.0017*rand(1); end y_C12 = rand(1);
end for j=1:i c(j,:)=Fx(b(j,:)); d(j,:)=FINDabc(c(j,:)); e(j,:)=(d(j,:)-b(j,:))./b(j,:); end P=[b,d,e]; plot(e(:,1),'b--o');hold on; plot(e(:,2),'g-.x'); plot(e(:,3),'r:+');
==================== %% ACATSP.m %% Ant Colony Algorithm for Traveling Salesman Problem %% ChengAihua,PLA Information Engineering University,ZhengZhou,China %% Email:aihuacheng@gmail.com %% All rights reserved %%------------------------------------------------------------------------%% 主要符号说明 %% C n 个城市的坐标,n×2 的矩阵 %% NC_max 最大迭代次数 %% m 蚂蚁个数 %% Alpha 表征信息素重要程度的参数 %% Beta 表征启发式因子重要程度的参数 %% Rho 信息素蒸发系数 %% Q 信息素增加强度系数 %% R_best 各代最佳路线 %% L_best 各代最佳路线的长度 %%===================================================== ==================== %%第一步:变量初始化 n=size(C,1);%n 表示问题的规模(城市个数) D=zeros(n,n);%D 表示完全图的赋权邻接矩阵 for i=1:n for j=1:n if i~=j D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5; else D(i,j)=eps; end D(j,i)=D(i,j); end end Eta=1./D;%Eta 为启发因子,这里设为距离的倒数 Tau=ones(n,n);%Tau 为信息素矩阵 Tabu=zeros(m,n);%存储并记录路径的生成 NC=1;%迭代计数器 R_best=zeros(NC_max,n);%各代最佳路线 L_best=inf.*ones(NC_max,1);%各代最佳路线的长度 L_ave=zeros(NC_max,1);%各代路线的平均长度 while NC<=NC_max%停止条件之一:达到最大迭代次数 %%第二步:将 m 只蚂蚁放到 n 个城市上 Randpos=[]; for i=1:(ceil(m/n))
function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta, Rho,Q) %%=====================================================
Randpos=[Randpos,randperm(n)]; end Tabu(:,1)=(Randpos(1,1:m))'; %%第三步:m 只蚂蚁按概率函数选择下一座城市,完成各自的周游 for j=2:n for i=1:m visited=Tabu(i,1:(j-1));%已访问的城市 J=zeros(1,(n-j+1));%待访问的城市 P=J;%待访问城市的选择概率分布 Jc=1; for k=1:n if length(find(visited==k))==0 J(Jc)=k; Jc=Jc+1; end end %下面计算待选城市的概率分布 for k=1:length(J) P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta); end P=P/(sum(P)); %按概率原则选取下一个城市 Pcum=cumsum(P); Select=find(Pcum>=rand); to_visit=J(Select(1)); Tabu(i,j)=to_visit; end end if NC>=2 Tabu(1,:)=R_best(NC-1,:); end %%第四步:记录本次迭代最佳路线 L=zeros(m,1); for i=1:m R=Tabu(i,:); for j=1:(n-1) L(i)=L(i)+D(R(j),R(j+1)); end L(i)=L(i)+D(R(1),R(n)); end L_best(NC)=min(L); pos=find(L==L_best(NC)); R_best(NC,:)=Tabu(pos(1),:);
10.function [Q,w] = REF(C,r,k) B=C(k,:);C(k,:)=[];A=zeros(1,3); for i=1:(k-1) A=A+C(i,:); end if k>=2 A=A/(k-1); end Q = A+r*(A-B); w=r/2; end 14 11.function Q_M = RESHAP_MAX(B,A,K) p=0;a=0; for i =wenku.baidu.com:K x(i)=COMP(B,A(i,:)); end j=max(x); for i=1:K
L_ave(NC)=mean(L); NC=NC+1 %%第五步:更新信息素 Delta_Tau=zeros(n,n); for i=1:m for j=1:(n-1) Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i); end Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i); end Tau=(1-Rho).*Tau+Delta_Tau; %%第六步:禁忌表清零 Tabu=zeros(m,n); end %%第七步:输出结果 Pos=find(L_best==min(L_best)); Shortest_Route=R_best(Pos(1),:) Shortest_Length=L_best(Pos(1)) subplot(1,2,1) DrawRoute(C,Shortest_Route) subplot(1,2,2) plot(L_best) hold on plot(L_ave) function DrawRoute(C,R) %%===================================================== ==================== %% DrawRoute.m %% 画路线图的子函数 %%------------------------------------------------------------------------%% C Coordinate 节点坐标,由一个 N×2 的矩阵存储 %% R Route 路线 %%===================================================== ==================== N=length(R); scatter(C(:,1),C(:,2)); hold on plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)]) hold on for ii=2:N
13 8.function Q_F = Fx(x) D=-7*1e-1;h=3*1e-1; B = -0.253; F = -0.259; a = x(1); b = x(2); c = x(3); s3 = 2*(c+D)/(1+a^2+b^2);
s2=-b*s3; s1=-a*s3; c13 = (c+D)/(1-b*B); c12 = B*c13; c23 = (c+D)/(1-b*F); c22 = F*c23; B11 = (h-s3)*(-s1)/(c13-s3)+s1; B12 = (h-s3)*(c12-s2)/(c13-s3)+s2; B21 = (h-s3)*(-s1)/(c23-s3)+s1; B22 = (h-s3)*(c22-s2)/(c23-s3)+s2; Q_F=[B11,B12,B21,B22]; 9.function A = RAND() i = 1; while(i<=6) A(i,1:2)= CREAT_RAND12; A(i,3)=CREAT_RAND3; if(D_ZONE(A(i,:))==1) i=i+1; end end
if(y_C12<=0.5) Q_C12(2) = 0.0017*rand(1); else Q_C12(2) = -0.0017*rand(1); End
5.function Q_C = CREAT_RAND3() y_C = rand(1); if(y_C<=0.5) Q_C = 3*1e-6*rand(1); else Q_C = -3*1e-6*rand(1); end 6.function y = D_ZONE(Q_D) p = pi/180/10; a=Q_D(1);b=Q_D(2);c=Q_D(3); if((abs(c)<3*1e-6)&((a^2+b^2)<(1/cos(p))^2-1)) y =1; else y=0; end 7.function y = D_ZONE1(Q_D) p = pi/180/10; a=Q_D(1);b=Q_D(2);c=Q_D(3); if((abs(c)<=3*1e-6)&((a^2+b^2)<=4*(1/cos(p))^2-1)) y =1; else y=0; end
if x(i)==j,p=i; end end a=A(p,:); A(p,:) = A(K,:); A(K,:)=a; Q_M=A; End
12.function w=STAND(A,R) D=-7*1e-1;h=3*1e-1; B = -.253; F = -.259; a = R(1); b = R(2); c = R(3); B11 = A(1); B12 = A(2); B21 = A(3); B22 = A(4); s3 = 2*(c+D)/(1+a^2+b^2); s2=-b*s3; s1=-a*s3; c13 = (c+D)/(1-b*B); c12 = B*c13; c23 = (c+D)/(1-b*F); c22 = F*c23; F11 = (h-s3)*(-s1)/(c13-s3)+s1; F12 = (h-s3)*(c12-s2)/(c13-s3)+s2; F21 = (h-s3)*(-s1)/(c23-s3)+s1; F22 = (h-s3)*(c22-s2)/(c23-s3)+s2; w(1)=abs(B11-F11); w(2)=abs(B12-F12); w(3)=abs(B21-F21); w(4)=abs(B22-F22);
3.function p = recognize(i) %分辨率的影响 j = 1; while(j<=i) a(1,1:2) = CREAT_RAND12 (); a(1,3) = CREAT_RAND3(); if(D_ZONE(a(1,:))) b(j,:) = a(1,:); j = j+1; end end for j = 1:i c(j,:) = Fx(b(j,:)); e(j,:)=1e-4*round(1e4*c(j,:)); end for j = 1:i f(j,1:3) = FINDabc(e(j,:)); end hl = (b-f)./b; for j = 1:3 h(:,j)=sum(abs(hl(:,j)))/i*ones(i,1); end plot(h(:,1),'--r');hold on;plot(h(:,2),'-y');plot(h(:,3),':k'); p = h; 12 4.function Q_C12 =CREAT_RAND12 () if(y_C12<=0.5) Q_C12(1) =0.0017 *rand(1); else Q_C12(1) = -0.0017*rand(1); end y_C12 = rand(1);
end for j=1:i c(j,:)=Fx(b(j,:)); d(j,:)=FINDabc(c(j,:)); e(j,:)=(d(j,:)-b(j,:))./b(j,:); end P=[b,d,e]; plot(e(:,1),'b--o');hold on; plot(e(:,2),'g-.x'); plot(e(:,3),'r:+');
==================== %% ACATSP.m %% Ant Colony Algorithm for Traveling Salesman Problem %% ChengAihua,PLA Information Engineering University,ZhengZhou,China %% Email:aihuacheng@gmail.com %% All rights reserved %%------------------------------------------------------------------------%% 主要符号说明 %% C n 个城市的坐标,n×2 的矩阵 %% NC_max 最大迭代次数 %% m 蚂蚁个数 %% Alpha 表征信息素重要程度的参数 %% Beta 表征启发式因子重要程度的参数 %% Rho 信息素蒸发系数 %% Q 信息素增加强度系数 %% R_best 各代最佳路线 %% L_best 各代最佳路线的长度 %%===================================================== ==================== %%第一步:变量初始化 n=size(C,1);%n 表示问题的规模(城市个数) D=zeros(n,n);%D 表示完全图的赋权邻接矩阵 for i=1:n for j=1:n if i~=j D(i,j)=((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5; else D(i,j)=eps; end D(j,i)=D(i,j); end end Eta=1./D;%Eta 为启发因子,这里设为距离的倒数 Tau=ones(n,n);%Tau 为信息素矩阵 Tabu=zeros(m,n);%存储并记录路径的生成 NC=1;%迭代计数器 R_best=zeros(NC_max,n);%各代最佳路线 L_best=inf.*ones(NC_max,1);%各代最佳路线的长度 L_ave=zeros(NC_max,1);%各代路线的平均长度 while NC<=NC_max%停止条件之一:达到最大迭代次数 %%第二步:将 m 只蚂蚁放到 n 个城市上 Randpos=[]; for i=1:(ceil(m/n))
function [R_best,L_best,L_ave,Shortest_Route,Shortest_Length]=ACATSP(C,NC_max,m,Alpha,Beta, Rho,Q) %%=====================================================
Randpos=[Randpos,randperm(n)]; end Tabu(:,1)=(Randpos(1,1:m))'; %%第三步:m 只蚂蚁按概率函数选择下一座城市,完成各自的周游 for j=2:n for i=1:m visited=Tabu(i,1:(j-1));%已访问的城市 J=zeros(1,(n-j+1));%待访问的城市 P=J;%待访问城市的选择概率分布 Jc=1; for k=1:n if length(find(visited==k))==0 J(Jc)=k; Jc=Jc+1; end end %下面计算待选城市的概率分布 for k=1:length(J) P(k)=(Tau(visited(end),J(k))^Alpha)*(Eta(visited(end),J(k))^Beta); end P=P/(sum(P)); %按概率原则选取下一个城市 Pcum=cumsum(P); Select=find(Pcum>=rand); to_visit=J(Select(1)); Tabu(i,j)=to_visit; end end if NC>=2 Tabu(1,:)=R_best(NC-1,:); end %%第四步:记录本次迭代最佳路线 L=zeros(m,1); for i=1:m R=Tabu(i,:); for j=1:(n-1) L(i)=L(i)+D(R(j),R(j+1)); end L(i)=L(i)+D(R(1),R(n)); end L_best(NC)=min(L); pos=find(L==L_best(NC)); R_best(NC,:)=Tabu(pos(1),:);
10.function [Q,w] = REF(C,r,k) B=C(k,:);C(k,:)=[];A=zeros(1,3); for i=1:(k-1) A=A+C(i,:); end if k>=2 A=A/(k-1); end Q = A+r*(A-B); w=r/2; end 14 11.function Q_M = RESHAP_MAX(B,A,K) p=0;a=0; for i =wenku.baidu.com:K x(i)=COMP(B,A(i,:)); end j=max(x); for i=1:K
L_ave(NC)=mean(L); NC=NC+1 %%第五步:更新信息素 Delta_Tau=zeros(n,n); for i=1:m for j=1:(n-1) Delta_Tau(Tabu(i,j),Tabu(i,j+1))=Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i); end Delta_Tau(Tabu(i,n),Tabu(i,1))=Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i); end Tau=(1-Rho).*Tau+Delta_Tau; %%第六步:禁忌表清零 Tabu=zeros(m,n); end %%第七步:输出结果 Pos=find(L_best==min(L_best)); Shortest_Route=R_best(Pos(1),:) Shortest_Length=L_best(Pos(1)) subplot(1,2,1) DrawRoute(C,Shortest_Route) subplot(1,2,2) plot(L_best) hold on plot(L_ave) function DrawRoute(C,R) %%===================================================== ==================== %% DrawRoute.m %% 画路线图的子函数 %%------------------------------------------------------------------------%% C Coordinate 节点坐标,由一个 N×2 的矩阵存储 %% R Route 路线 %%===================================================== ==================== N=length(R); scatter(C(:,1),C(:,2)); hold on plot([C(R(1),1),C(R(N),1)],[C(R(1),2),C(R(N),2)]) hold on for ii=2:N
13 8.function Q_F = Fx(x) D=-7*1e-1;h=3*1e-1; B = -0.253; F = -0.259; a = x(1); b = x(2); c = x(3); s3 = 2*(c+D)/(1+a^2+b^2);
s2=-b*s3; s1=-a*s3; c13 = (c+D)/(1-b*B); c12 = B*c13; c23 = (c+D)/(1-b*F); c22 = F*c23; B11 = (h-s3)*(-s1)/(c13-s3)+s1; B12 = (h-s3)*(c12-s2)/(c13-s3)+s2; B21 = (h-s3)*(-s1)/(c23-s3)+s1; B22 = (h-s3)*(c22-s2)/(c23-s3)+s2; Q_F=[B11,B12,B21,B22]; 9.function A = RAND() i = 1; while(i<=6) A(i,1:2)= CREAT_RAND12; A(i,3)=CREAT_RAND3; if(D_ZONE(A(i,:))==1) i=i+1; end end
if(y_C12<=0.5) Q_C12(2) = 0.0017*rand(1); else Q_C12(2) = -0.0017*rand(1); End
5.function Q_C = CREAT_RAND3() y_C = rand(1); if(y_C<=0.5) Q_C = 3*1e-6*rand(1); else Q_C = -3*1e-6*rand(1); end 6.function y = D_ZONE(Q_D) p = pi/180/10; a=Q_D(1);b=Q_D(2);c=Q_D(3); if((abs(c)<3*1e-6)&((a^2+b^2)<(1/cos(p))^2-1)) y =1; else y=0; end 7.function y = D_ZONE1(Q_D) p = pi/180/10; a=Q_D(1);b=Q_D(2);c=Q_D(3); if((abs(c)<=3*1e-6)&((a^2+b^2)<=4*(1/cos(p))^2-1)) y =1; else y=0; end