蚁群聚类算法

1.问题描述:

       为了提高蚁群化学聚类算法的聚类效果 ,通过引入径向基核函数改进蚁群化学聚类算法的相似度度量 式 ,提出了核蚁群化学聚类算法。将核蚁群化学聚类算法用于三个标准数据集合 ,计算结果表明与蚁群化学聚 类算法相比 ,核蚁群化学聚类算法聚类效果提升显著。将核蚁群化学聚类算法、核自组织神经网络映射算法和 基于多项式核的结构化有向树数据聚类算法同时用于 Iris数据集合 ,结果显示三种核聚类算法聚类效果相当。

2.部分程序:

 

clc;
clf;
clear;
 % X = 测试样本矩阵;
X1 = load('data.txt');
X = X1(:,1:2);

[N,n]=size(X);      % N =测试样本数;n =测试样本的属性数;
K = 4;              % K = 组数; 
R = 100;            % R = 蚂蚁数; 
t_max = 1000;       % t_max =最大迭代次数;           
% 初始化
c = 10^-2;
tau = ones(N,K) * c;    %信息素矩阵,初始值为0.01的N*K矩阵(样本数*聚类数)
q = 0.9;                % 阈值q
rho = 0.1;              % 蒸发率
best_solution_function_value = inf; % 最佳路径度量值(初值为无穷大,该值越小聚类效果越好)
tic                     %计算程序运行时间
t = 1; 
%=======程序终止条件(下列两个终止条件任选其一)======
 while ((t<=t_max))                             %达到最大迭代次数而终止
% while ((best_solution_function_value>=19727))  %达到一定的聚类效果而终止
%while ((best_solution_function_value>=1000))    
%========================= 
    %路径标识字符:标识每只蚂蚁的路径
    solution_string = zeros(R,N+1);     
    for i = 1 : R                       %以信息素为依据确定蚂蚁的路径
        r = rand(1,N);    %随机产生值为0-1随机数的1*51的数组
        for g = 1 : N
            if r(g) < q     %如果r(g)小于阈值
                tau_max = max(tau(g,:));
                Cluster_number = find(tau(g,:)==tau_max);   %聚类标识数,选择信息素最多的路径
                solution_string(i,g) = Cluster_number(1);   %确定第i只蚂蚁对第g个样本的路径标识
            else            %如果r(g)大于阈值,求出各路径信息素占在总信息素的比例,按概率选择路径
                sum_p = sum(tau(g,:)); 
                p = tau(g,:) / sum_p; 
                for u = 2 : K 
                    p(u) = p(u) + p(u-1); 
                end
               rr = rand;          
                for s = 1 : K 
                    if (rr <= p(s)) 
                       Cluster_number = s;
                       solution_string(i,g) = Cluster_number;  
                    break; 
                end 
            end
        end
    end

    % 计算聚类中心
    weight = zeros(N,K);
       for h = 1:N              %给路径做计算标识
           Cluster_index = solution_string(i,h); %类的索引编号          
           weight(h,Cluster_index) = 1;          %对样本选择的类在weight数组的相应位置标1
       end

       cluster_center = zeros(K,n);  %聚类中心(聚类数K个中心)
       for j = 1:K
           for v = 1:n
               sum_wx = sum(weight(:,j).*X(:,v));   %各类样本各属性值之和
               sum_w = sum(weight(:,j));            %各类样本个数
               if sum_w==0                          %该类样本数为0,则该类的聚类中心为0
                 cluster_center(j,v) =0
                  continue;
               else                                 %该类样本数不为0,则聚类中心的值取样本属性值的平均值
               cluster_center(j,v) = sum_wx/sum_w;
               end
            end
       end

    % 计算各样本点各属性到其对应的聚类中心的均方差之和,该值存入solution_string的最后一位
      F = 0;
      for j= 1:K
          for ii = 1:N
              Temp=0;
              if solution_string(i,ii)==j;                
                  for v = 1:n
                      Temp = ((abs(X(ii,v)-cluster_center(j,v))).^2)+Temp;
                  end
                  Temp = sqrt(Temp);
              end
            F = (Temp)+F;
          end        
      end

       solution_string(i,end) = F;                      

    end 
    %根据F值,把solution_string矩阵升序排序
    [fitness_ascend,solution_index] = sort(solution_string(:,end),1);
    solution_ascend = [solution_string(solution_index,1:end-1) fitness_ascend];

    pls = 0.1;        %局部寻优阈值pls(相当于变异率)
    L = 2;              % 在L条路径内局部寻优
    % 局部寻优程序
    solution_temp = zeros(L,N+1);
    k = 1;
    while(k <= L)
           solution_temp(k,:) = solution_ascend(k,:);
           rp = rand(1,N);     %产生一个1*N(51)维的随机数组,某值小于pls则随机改变其对应的路径标识
           for i = 1:N
               if rp(i) <= pls
                   current_cluster_number = setdiff([1:K],solution_temp(k,i));
                   rrr=randint(1,1,[1,K-1]);
                   change_cluster = current_cluster_number(rrr);
                   solution_temp(k,i) = change_cluster;
               end
           end

        % 计算临时聚类中心   
           solution_temp_weight = zeros(N,K);
           for h = 1:N
               solution_temp_cluster_index = solution_temp(k,h);           
               solution_temp_weight(h,solution_temp_cluster_index) = 1;
           end

           solution_temp_cluster_center = zeros(K,n);
           for j = 1:K
               for v = 1:n
                   solution_temp_sum_wx = sum(solution_temp_weight(:,j).*X(:,v));
                   solution_temp_sum_w = sum(solution_temp_weight(:,j));
                   if solution_temp_sum_w==0
                   solution_temp_cluster_center(j,v) =0;
                   continue;
                   else
                       solution_temp_cluster_center(j,v) = solution_temp_sum_wx/solution_temp_sum_w;
                   end
               end
          end
          % 计算各样本点各属性到其对应的临时聚类中心的均方差之和Ft;
          solution_temp_F = 0;
          for j= 1:K
              for ii = 1:N
                  st_Temp=0;
                  if solution_temp(k,ii)==j;                               
                      for v = 1:n
                          st_Temp = ((abs(X(ii,v)-solution_temp_cluster_center(j,v))).^2)+st_Temp;
                      end
                      st_Temp = sqrt(st_Temp);
                  end
                  solution_temp_F = (st_Temp)+solution_temp_F;
              end
          end
        solution_temp(k,end) = solution_temp_F;   
        %根据临时聚类度量调整路径
        % 如果 Ft<Fl 则 Fl=Ft , Sl=St
          if solution_temp(k,end) <= solution_ascend(k,end)              
              solution_ascend(k,:) = solution_temp(k,:);               
          end

          if solution_ascend(k,end)<=best_solution_function_value
              best_solution = solution_ascend(k,:);
          end
      k = k+1;
      end   

    % 用最好的L条路径更新信息数矩阵
    tau_F = 0;
    for j = 1:L    
       tau_F = tau_F + solution_ascend(j,end);
    end
    for i = 1 : N        
       tau(i,best_solution(1,i)) = (1 - rho) * tau(i,best_solution(1,i)) + 1/ tau_F; 
    %1/tau_F和rho/tau_F效果都很好
    end 
    t=t+1
    best_solution_function_value =  solution_ascend(1,end);
    best_solution_function_value
end
time=toc;       %输出程序运行时间
clc

time
cluster_center
best_solution = solution_ascend(1,1:end-1);
IDY=ctranspose(best_solution)
best_solution_function_value =  solution_ascend(1,end)
%分类结果显示
plot(cluster_center(:,1),cluster_center(:,2),'o');grid;box
title('蚁群聚类结果(R=100,t=1000)')
xlabel('X')
ylabel('Y')
YY=[1 2 3 4];
index1 = find(YY(1) == best_solution)
index2 = find(YY(2) == best_solution)
index3 = find(YY(3) == best_solution)
index4 = find(YY(4) == best_solution)
line(X(index1,1),X(index1,2),X(index1,3),'linestyle','none','marker','*','color','g');
line(X(index2,1),X(index2,2),X(index2,3),'linestyle','none','marker','*','color','r');
line(X(index3,1),X(index3,2),X(index3,3),'linestyle','none','marker','+','color','b');
line(X(index4,1),X(index4,2),X(index4,3),'linestyle','none','marker','s','color','b');
rotate3d

3.仿真结论:

D18

%-- Unknown date --% else p(:,j)=0; end; if maxp(1)<p(1,j) maxp(1)=p(1,j); end; linear_index=find(maxp(1)==p(1,:)); size1=[1,n]; [r_index,c_index]=ind2sub(size1,linear_index(1)); solution_medium(k,1)=distance(g(NC,k),c_index(1)); route(k,1)=c_index(1); tabu(k,c_index(1))=1; tao(g(NC,k),c_index(1))=(1-rou)*tao(g(NC,k),c_index(1))+rou*initao(g(NC,k),c_index(1));%local updating for the first itertion tao(c_index(1),g(NC,k))=(1-rou)*tao(c_index(1),g(NC,k)))+rou*initao(c_index(1),g(NC,k)); solution(k)=solution_medium(k,1); for s=2:(n-1) c_index(s)=0; r_index(s)=0; linear_index(s)=0; maxp(s)=0; for j=1:n if tabu(k,j)==0 psum_medium0(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta); else psum_medium0(s,j)=0; end; psum_medium=psum_medium0.'; psum(k,s)=sum(psum_medium(:,s)); for j=1:n if tabu(k,j)==0; p(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta)/psum(k,s); else p(s:(n-1),j)=0; end; if maxp(s)<p(s,j) maxp(s)=p(s,j); end; linear_index=find(maxp(s)==p); size2=[n-1,n]; [r_index(s),c_index(s)]=ind2sub(size2,linear_index(1)); solution_medium(k,s)=distance(c_index(s-1),c_index(s)); route(k,s)=c_index(s); tabu(k,c_index(s))=1; tao(c_index(s-1),c_index(s))=(1-rou)*tao(c_index(s-1),c_index(s))+rou*initao(c_index(s-1),c_index(s)); tao(c_index(s),c_index(s-1))=(1-rou)*tao(c_index(s),c_index(s-1)))+rou*initao(c_index(s),c_index(s-1)); solution(k)=solution_medium(k,s); end; tao(c_index(n-1),g(NC,k))=(1-rou)*tao(c_index(n-1),g(NC,k)))+rou*initao(c_index(n-1),g(NC,k)); tao(g(NC,k),c_index(n-1))=(1-rou)*tao(g(NC,k),c_index(n-1))+rou*initao(g(NC,k),c_index(n-1)); solution_medium(k,n)=distance(c_index(n-1),g(NC,k)); solution(k)=solution(k)+solution_medium(k,n); solution_NC(NC,k)=solution(k); besolution_NC(NC,k)=solution_NC(NC,1); end; for k=1:m if besolution_NC(NC)=solution_NC(NC,k) besolution_NC(NC)=solution_NC(NC,k); end; %% ******** In this phase globle updating occurs ********%% linear_index1=find(besolution_NC(NC)==solution_NC(NC,:)); size3=[NC,m]; [r_index1(NC),c_index1(NC)]=ind2sub(size3,linear_index(1)); bestroute_NC(NC,:)=route(c_index(NC),:); [aa,bb]=size(linaer_index1); for i=1:bb [r_index1_tao(NC,i),c_index1_t(NC,i)]=ind2sub(size3,linear_index(i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=Q/solution(c_index1_t(NC,i)); detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1)); tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i))); detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i))); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1)); for s=2:n-1 detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=Q/solution(c_index1_t(NC,i)); detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=(1-rou).*tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))+rou.*detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1)); tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=(1-rou),*tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))+rou.*detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s)); end; if bestsolution>bestsolution_NC(NC) bestsolution=bestsolution_NC(NC); end; linear_index2=find(bestsolution==bestsolution_NC); size4[1,NC]; [r_index2,c_index2]=ind2sub(size4,linear_index2(1)); bestroute(1,:)=bestroute_NC(c_index2,:); initao=tao; end; %% ******** Otuput the best rusult of TSP ******** %% for NC-1:NCmax bestroute_NC(NC,n)=g(NC,c_index1(NC)); t(NC)=NC; end; bestroute(1,n)=bestroute_NC(c_index2,n); polt(x(bestroute(n)),y(bestroute(n)),'*'); hold on; end; line([x(bestroute(n)),x(bestroute(1))],[y(bestroute(n)),y(bestroute(1))]); hold on; for j=1:(n-1) line([x(bestroute(j)),x(bestroute(j+1))],[y(bestroute(j)),y(bestroute(j+1))]); hold on; end; hold off; xlabel('X coordinate'); ylabel('Y coordinate'); title('best tour in NCmax iterations'); %% ********Funtion: Open and improt file of city coordinate in TSP %% ******** %% [fname,pname,fiterindex]=uigetfile('*.*','All Files(*.*)'); set(handles.text_filename,'string',filename); fid=fopen(filename,'r')' if fid==-1; warndlg('can't open the file','WARN'); fclose(fid); end; [A0,COUNT]=fscanf(fid,'%g'); n=COUNT/3; set(handles.edit_citysum,'string',n); fclose(fid); m=str2double(get(handles.edit_citysum,'string')); NCmax=str2double(get(handles.edit_ncmax,'string')); for NC=1:NCmax for k=1:m g(NC,k)=fix((n-1)*rand(1))+1; end; end %-- 10-9-1 下午4:54 --% %-- 10-9-1 下午8:09 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\matlab.mat') %-- 10-9-2 上午9:01 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\matlab.mat') %% ******Initialization phase ****** %% global A0 global n; %city number global g; m=str2double(get(handles.edit_antsum, 'String')); %set ant number by using Matlab GUI initao=str2double(get(handles.edit_tao, 'String')); alpha=str2double(get(handles.edit_alpha, 'String')); beta=str2double(get(handles.edit_beta, 'String')); Q=str2double(get(handles.edit_q, 'String')); rou=str2double(get(handles.edit_rou, 'String')); NCmax=str2double(get(handles.edit_ncmax, 'String')); A=reshape(A0,3,n); %get the city coordinates in TSP x=A(2,:); %get x=coordinate y=A(3,:); %get y=coordinate for i=1:n for j=1:n distance(i,j)=sqrt((x(i)-x(j)+(y(i)-y(j))*y(i)-y(j))); end; for i=1:n for j=1:n if j!=i tao(i,j)=initao; yita(i,j)=1/distance(i,j); end; initao=initao.*ones(n,n); detatao=zeors(n,n); bestsolution=10000000000000; %% ******This is in phase in which ants build tours ****** %% for NC=1:NCmax tabu=zeros(m.n)l for k=1:m tabu(k,g(NC,k))=1; maxp(1)=0; for j=1:n if tabu(k,j)==0 psum_medium0(1,j)=(tao(g(NC,k)^alpha).*(yita(g(NC,k),j)^beta); else psum_medium0(1,j)=0; end; psum_medium=psum_medium0.'; psum(k,1)=sum(psum_medium(:,1)); for j=1:n if tabu(k,j)==0 p(1,j)=(tao(g(NC,k)^alpha).*(yita(g(NC,k),j)^beta/psm(k,1); else p(:,j)=0; end; if maxp(1)<p(1,j) maxp(1)=p(1,j); end; linear_index=find(maxp(1)==p(1,:)); size1=[1,n]; [r_index,c_index]=ind2sub(size1,linear_index(1)); solution_medium(k,1)=distance(g(NC,k),c_index(1)); route(k,1)=c_index(1); tabu(k,c_index(1))=1; tao(g(NC,k),c_index(1))=(1-rou)*tao(g(NC,k),c_index(1))+rou*initao(g(NC,k),c_index(1));%local updating for the first itertion tao(c_index(1),g(NC,k))=(1-rou)*tao(c_index(1),g(NC,k)))+rou*initao(c_index(1),g(NC,k)); solution(k)=solution_medium(k,1); for s=2:(n-1) c_index(s)=0; r_index(s)=0; linear_index(s)=0; maxp(s)=0; for j=1:n if tabu(k,j)==0 psum_medium0(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta); else psum_medium0(s,j)=0; end; psum_medium=psum_medium0.'; psum(k,s)=sum(psum_medium(:,s)); for j=1:n if tabu(k,j)==0; p(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta)/psum(k,s); else p(s:(n-1),j)=0; end; if maxp(s)<p(s,j) maxp(s)=p(s,j); end; linear_index=find(maxp(s)==p); size2=[n-1,n]; [r_index(s),c_index(s)]=ind2sub(size2,linear_index(1)); solution_medium(k,s)=distance(c_index(s-1),c_index(s)); route(k,s)=c_index(s); tabu(k,c_index(s))=1; tao(c_index(s-1),c_index(s))=(1-rou)*tao(c_index(s-1),c_index(s))+rou*initao(c_index(s-1),c_index(s)); tao(c_index(s),c_index(s-1))=(1-rou)*tao(c_index(s),c_index(s-1)))+rou*initao(c_index(s),c_index(s-1)); solution(k)=solution_medium(k,s); end; tao(c_index(n-1),g(NC,k))=(1-rou)*tao(c_index(n-1),g(NC,k)))+rou*initao(c_index(n-1),g(NC,k)); tao(g(NC,k),c_index(n-1))=(1-rou)*tao(g(NC,k),c_index(n-1))+rou*initao(g(NC,k),c_index(n-1)); solution_medium(k,n)=distance(c_index(n-1),g(NC,k)); solution(k)=solution(k)+solution_medium(k,n); solution_NC(NC,k)=solution(k); besolution_NC(NC,k)=solution_NC(NC,1); end; for k=1:m if besolution_NC(NC)=solution_NC(NC,k) besolution_NC(NC)=solution_NC(NC,k); end; %% ******** In this phase globle updating occurs ********%% linear_index1=find(besolution_NC(NC)==solution_NC(NC,:)); size3=[NC,m]; [r_index1(NC),c_index1(NC)]=ind2sub(size3,linear_index(1)); bestroute_NC(NC,:)=route(c_index(NC),:); [aa,bb]=size(linaer_index1); for i=1:bb [r_index1_tao(NC,i),c_index1_t(NC,i)]=ind2sub(size3,linear_index(i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=Q/solution(c_index1_t(NC,i)); detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),1)); tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),1),g(NC,c_index1_t(NC,i))); detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=Q/solution(c_index1_t(NC,i)); detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))=(1-rou).* tao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i)))+rou*detatao(route(c_index1_t(NC,i),n-1),g(NC,c_index1_t(NC,i))); tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))=(1-rou).* tao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1))+rou*detatao(g(NC,c_index1_t(NC,i)),route(c_index1_t(NC,i),n-1)); for s=2:n-1 detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=Q/solution(c_index1_t(NC,i)); detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=Q/solution(c_index1_t(NC,i)); tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))=(1-rou).*tao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1))+rou.*detatao(route(c_index1_t(NC,i),s),route(c_index1_t(NC,i),s-1)); tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))=(1-rou),*tao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s))+rou.*detatao(route(route(c_index1_t(NC,i),s-1),c_index1_t(NC,i),s)); end; if bestsolution>bestsolution_NC(NC) bestsolution=bestsolution_NC(NC); end; linear_index2=find(bestsolution==bestsolution_NC); size4[1,NC]; [r_index2,c_index2]=ind2sub(size4,linear_index2(1)); bestroute(1,:)=bestroute_NC(c_index2,:); initao=tao; end; %% ******** Otuput the best rusult of TSP ******** %% for NC-1:NCmax bestroute_NC(NC,n)=g(NC,c_index1(NC)); t(NC)=NC; end; bestroute(1,n)=bestroute_NC(c_index2,n); polt(x(bestroute(n)),y(bestroute(n)),'*'); hold on; end; line([x(bestroute(n)),x(bestroute(1))],[y(bestroute(n)),y(bestroute(1))]); hold on; for j=1:(n-1) line([x(bestroute(j)),x(bestroute(j+1))],[y(bestroute(j)),y(bestroute(j+1))]); hold on; end; hold off; xlabel('X coordinate'); ylabel('Y coordinate'); title('best tour in NCmax iterations'); %% ********Funtion: Open and improt file of city coordinate in TSP %% ******** %% [fname,pname,fiterindex]=uigetfile('*.*','All Files(*.*)'); set(handles.text_filename,'string',filename); fid=fopen(filename,'r')' if fid==-1; warndlg('can't open the file','WARN'); fclose(fid); end; [A0,COUNT]=fscanf(fid,'%g'); n=COUNT/3; set(handles.edit_citysum,'string',n); fclose(fid); m=str2double(get(handles.edit_citysum,'string')); NCmax=str2double(get(handles.edit_ncmax,'string')); for NC=1:NCmax for k=1:m g(NC,k)=fix((n-1)*rand(1))+1; end; end %-- 10-9-2 下午4:47 --% %-- 10-9-2 下午8:02 --% %{ ? 当前蚂蚁移到邻近区域内的没有被其他蚂蚁占据的节点 a) 输入参数 ? 蚂蚁编号ant_no ? S局部查找范围 ? ant_matrix(蚂蚁的窗格矩阵) ? Z平面窗格总区域 b) 输出参数 ? 蚂蚁前往的结点坐标 %} function ant_matrix=ant_move(ant_no,S,ant_matrix,Z) %1、获得ant_no点坐标 %[x,y,z]=ant_matrix(ant_no,:); one_ant=ant_matrix(ant_no,:); x=one_ant(1); y=one_ant(2); z=one_ant(3); %2、获得上下界限 x_low=x-S/2; x_high=x+S/2; if(x_low<0) x_low=0; end if(x_high>Z) x_high=Z; end y_low=y-S/2; y_high=y+S/2; if(y_low<0) y_low=0; end if(y_high>Z) y_high=Z; end %获得所有邻接结点下标数组 %{ allneighbours=[]; [row,col]=size(ant_matrix); for i=1:row if i~=ant_no [x,y,z]=ant_matrix(i,:); if x>=x_low && x<=x_high && y>=y_low && y<=y_high allneighbours=[allneighbours i]; end end end allneighbours=[allneighbours ant_no]; %} [row,col]=size(ant_matrix); positions=[]; %遍历所有上下界之间的点 for i=x_low:x_high for j=y_low:y_high status=0; for k=1:row % [x,y,z]=ant_matrix(k,:); one_ant=ant_matrix(k,:); x=one_ant(1); y=one_ant(2); z=one_ant(3); if x==i && y==j status=1; end end if status==0 positions=[positions;i j]; end end end %随机产生这个位置 [row,col]=size(positions); i=ceil(rand*row); position=positions(i,:); ant_matrix(ant_no,1)=position(1); ant_matrix(ant_no,2)=position(2); %{ ? 求Oi点在平面窗格内S*S区域内的邻接点的下标矩阵 a) 输入参数: ? Oi点 ? S局部查找范围 ? Item_Window(结点的窗格矩阵) ? Z平面窗格总区域 b) 输出参数 ? 所有在S*S区域内的结点的下标 具体算法: 1、S区域的上界 2、S区域的下界 3、比较一下点在哪些区域内 %} function allneighbours=get_allneighbours(Oi,S,item_window,Z) %1、获得Oi点坐标 %[x,y]=item_window(Oi,:); one_item=item_window(Oi,:); x=one_item(1); y=one_item(2); %2、获得上下界限 x_low=x-S/2; x_high=x+S/2; if(x_low<0) x_low=0; end if(x_high>Z) x_high=Z; end y_low=y-S/2; y_high=y+S/2; if(y_low<0) y_low=0; end if(y_high>Z) y_high=Z; end %获得所有邻接结点下标数组 allneighbours=[]; [row,col]=size(item_window); for i=1:row if i~=Oi % [x,y]=item_window(i,:); one_item=item_window(i,:); x=one_item(1); y=one_item(2); if x>=x_low && x<=x_high && y>=y_low && y<=y_high allneighbours=[allneighbours i]; end end end %{ ? 求两点空间的距离函数 a) 输入参数: ? Oi点 ? Oj点 ? 点的空间矩阵Item_Space b) 输出参数: ? 两点之间的距离 %} function distance=get_distance(Oi,Oj,item_space) distance=0.0; X1=item_space(Oi,:); X2=item_space(Oj,:); size=length(X1); for i=1:size distance=distance+(X1(i)-X2(i))^2; end distance=sqrt(distance); %{ ? 设计函数计算群体相似度计算 a) 输入参数: ? Oi点(需要计算的点)、 ? S(S*S区域内)、 ? Alpha、 ? Item_Window(结点的窗格矩阵) ? Z(Z*Z区域) ? item_space(结点的空间矩阵) b) 输出参数: ? Oi点的相似度 程序流程: 1、计算Oi的所有邻居 2、所有邻居的距离相应的和 3、判断和值和0的关系 4、最后决定值的大小 %} function fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space) %1、计算Oi的所有邻居 allneighbours=get_allneighbours(Oi,S,item_window,Z); %2、所有邻居的距离相应的和 len=length(allneighbours); fi=0; for i=1:len distance=get_distance(Oi,allneighbours(i),item_space); fi=fi+(1-distance)/Alpha; end if(fi>0) fi=fi/(S^2); else fi=0; end %{ ? 放下概率的计算 a) 输入参数: ? Oi点(需要计算的点)、 ? S(S*S区域内)、 ? Alpha、 ? Item_Window(结点的窗格矩阵) ? K2 ? Z(Z*Z区域) ? item_space(结点的空间矩阵) b) 输出参数 ? 放下概率 程序流程 1、先计算群体相似概率 2、计算放下概率 %} function Pd=get_Pd(Oi,S,Alpha,item_window,K2,Z,item_space) %1、计算群体相似概率 fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space); %2、计算放下概率 if fi %{ ? 拾起概率的计算 a) 输入参数: ? Oi点(需要计算的点)、 ? S(S*S区域内)、 ? Alpha、 ? Item_Window(结点的窗格矩阵) ? K2 b) 输出参数 ? 拾起概率 程序流程: 1、首先计算群体相似概率 2、计算拾起概率 %} function Pp=get_Pp(Oi,S,Alpha,item_window,K1,Z,item_space) %1 fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space); %2 Pp=(K1/(K1+fi))^2; %{ ? 判断蚂蚁所在处是否有点 a) 输入参数 ? x蚂蚁横坐标 ? y蚂蚁纵坐标 ? item_window所有点的坐标矩阵 b) 输出参数 ? 蚂蚁所在处是否有点,有点返回这个点在item_window中的行号,无点0 %} function position=has_item(x,y,item_window) %{ [row,col]=size(item_window); for i=1:row end %} %position=find(item_window(find(item_window(:,1)==x),2)==y) position=0; [row,col]=size(item_window); for i=1:row if ~isempty(find(x==item_window(i,1))) && ~isempty(find(y==item_window(i,2))) position=i; end end %{ ? 初始化函数 a) 输入参数 ? 蚂蚁的数目ant_number ? 点的数目item_number ? 空间尺寸Z,空间大小Z*Z b) 输出参数 ? 蚂蚁的平面窗格矩阵 ? 点的平面窗格矩阵 %} function [ant_matrix,item_window]=initialize(ant_number,item_number,Z) ant_matrix_x=randperm(Z); ant_matrix_y=randperm(Z); item_matrix_x=randperm(Z); item_matrix_y=randperm(Z); %生成蚂蚁矩阵 for i=1:ant_number ant_matrix(i,1)=ant_matrix_x(i); ant_matrix(i,2)=ant_matrix_y(i); ant_matrix(i,3)=0; end for i=1:item_number item_window(i,1)=item_matrix_x(i); item_window(i,2)=item_matrix_y(i); end %{ test 测试程序 整个程序流程 1、初始化: 2、迭代tmax次 3、所有的蚂蚁运动一次 4、产生一个0-1之间的随机数R 5、如果当前蚂蚁处于未负载状态,而且当前蚂蚁所在处的有点Oi 5.1、计算群体相似度f(Oi)和拾起概率Pp(Oi) 5.2、如果拾起概率Pp(Oi)》R 5.2.1、当前蚂蚁拾起点Oi(注意Oi在窗格中的位置是不断变动的) 5.3 5.2结束 6、如果条件5不成立,如果当前蚂蚁处于负载状态,持有点Oi,而且当前位置没有其他点 6.1计算群体相似度f(Oi)和放下概率Pd(Oi) 6.2如果放下概率Pd(Oi)》R 6.2.1放下节点Oi(注意点Oi在窗格中的位置是不断变动的) 6.2.2修改蚂蚁状态为卸载 6.3 6.2结束 7、5结束 8、当前蚂蚁移到邻近区域内的没有被其他蚂蚁占据的节点 9、所有的蚂蚁运动一次结束 10、迭代tmax次结束 输入参数: 1、ant_number 蚂蚁数 2、item_number 点数 3、Z 区域范围 4、tmax 最大迭代次数 5、S 邻接区域 6、Alpha 相异度参数 7、item_space 点的空间矩阵 8、K1 拾起概率参数 9、K2 放下概率参数 %} clear; clc; K1=0.1; K2=0.15; Alpha=0.15; %S=6; S=30; tmax=300; ant_number=20; Z=50; axis([0 Z 0 Z]); item_space=[ 5.1 3.5 1.4 0.2; 4.9 3.0 1.4 0.2; 4.7 3.2 1.3 0.2; 4.6 3.1 1.5 0.2; 5.0 3.6 1.4 0.2; 5.4 3.9 1.7 0.4; 4.6 3.4 1.4 0.3; 5.0 3.4 1.5 0.2; ]; [row,col]=size(item_space); item_number=row; %1、初始化: [ant_matrix,item_window]=initialize(ant_number,item_number,Z); item_window %2、迭代tmax次 for i=1:tmax %3、所有的蚂蚁运动一次 for j=1:ant_number %4、产生一个0-1之间的随机数R r=rand; %5、如果当前蚂蚁处于未负载状态,而且当前蚂蚁所在处的有点Oi %[x,y,z]=ant_matrix(j,:); one_ant=ant_matrix(j,:); x=one_ant(1); y=one_ant(2); z=one_ant(3); %rows1=find(item_window(:,1)==x); %rows2=find(item_window(:,2)==y); %Oi=find(item_window(find(item_window(:,1)==x),2)==y); Oi=has_item(x,y,item_window); if z==0 && Oi>0 %5.1、计算群体相似度f(Oi)和拾起概率Pp(Oi) %fi=get_Fi(Oi,S,Alpha,item_window,Z,item_space); Pp=get_Pp(Oi,S,Alpha,item_window,K1,Z,item_space); Pp; %5.2、如果拾起概率Pp(Oi)》R if Pp>=r %5.2.1、当前蚂蚁拾起点Oi(注意Oi在窗格中的位置是不断变动的) ant_matrix(j,3)=Oi; end elseif ant_matrix(j,3)>0 && Oi==0 %6、如果条件5不成立,如果当前蚂蚁处于负载状态,持有点Oi,而且当前位置没有其他点 %6.1计算群体相似度f(Oi)和放下概率Pd(Oi) Oi=ant_matrix(j,3);%代表当前蚂蚁所携带的点的行号 Pd=get_Pd(Oi,S,Alpha,item_window,K2,Z,item_space); Pd; %6.2如果放下概率Pd(Oi)》R if Pd>=r %6.2.1放下节点Oi(注意点Oi在窗格中的位置是不断变动的) item_window(Oi,1)=x; item_window(Oi,2)=y; ant_matrix(j,3)=0; end end %8、当前蚂蚁移到邻近区域内的没有被其他蚂蚁占据的节点 ant_matrix=ant_move(j,S,ant_matrix,Z); end end item_window axis([0 Z 0 Z]); plot(item_window(:,1),item_window(:,2),'--rs'); %-- 10-9-2 下午9:11 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\matlab.mat') %-- 10-9-3 下午2:54 --% load('E:\苏雪敬个人资料\苏雪敬-开题报告\程序代码\lf蚁群算法.mat') matlab7
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

fpga和matlab

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值