m基于ACO蚁群优化的FCM模糊聚类算法matlab仿真

目录

1.算法概述

2.仿真效果预览

3.核心MATLAB程序

4.完整MATLAB程序


1.算法概述

       蚁群算法是通过对自然界中真实蚂蚁的集体行为的观察,模拟而得到一种仿生优化算法,它具有很好的并行性,分布性.根据蚂蚁群体不同的集体行为特征,蚁群算法可分为受蚂蚁觅食行为启发的模型和受孵化分类启发的模型,受劳动分工和协作运输启发的模型.本文重点研究了前两种蚁群算法模型. 受蚂蚁觅食行为启发的模型又称为蚁群优化算法(ACO),是继模拟退火算法,遗传算法,禁忌搜索等之后又一启发式智能优化算法.目前它已成功应用于求解TSP问题,地图着色,路径车辆调度等优化问题.本文针对蚁群算法收敛时间长,易陷入局部最优的缺点,通过对路径上信息素的更新方式作出动态调整,建立信息素平滑机制,进而使得不同路径上的信息素的更新速度有所不同,从而使改进后算法能够有效地缩短搜索的时间,并能对最终解进行优化,避免过早的陷入局部最优. 聚类是数据挖掘的重要技术之一,它可按照某种规则将数据对象划分为多个类或簇,使同一类的数据对象有较高的相似度,而不同类的数据对象差异较大.       

       用FCM算法实现基于目标函数的模糊聚类又称交替的迭代优化法。迭代优化本质上属于局部搜索的爬山法,很容易陷入局部极值点,因此对初始化很敏感。通常是根据一定的经验准则选取初始参数,这样计算结果与初始参数设置是否恰当密切相关。特别是在数据量较大和高维情况下,设置合理的参数非常困难,只能通过多次实验比较选定。由于初始聚类中心和样本的输入次序对最终的结果有重大影响,往往是用若干不同的初始中心和聚类数目分别聚类,然后选择最满意的聚类作为最终的结果。 通过蚁群算法,我们可以得到最佳的初始聚类中心,然后进行快速的聚类。

2.仿真效果预览

matlab2022a仿真结果如下:

 

 

 

3.核心MATLAB程序

function cluster_center = func_ant_opt(data,Class_Num);

X     = data;
%N =测试样本数;n =测试样本的属性数;
[N,n] = size(X); 
%K = 组数; 
K     = Class_Num; 
%R = 蚂蚁数; 
R     = 200;      
%t_max =最大迭代次数;  
t_max = 100;               
% 初始化
c     = 10^-1;
tau   = ones(N,K) * c;    %信息素矩阵,初始值为0.001的N*K矩阵(样本数*聚类数)
q     = 0.90;                % 阈值q
rho   = 0.1;              % 蒸发率
best_solution_function_value = inf; % 最佳路径度量值(初值为无穷大,该值越小聚类效果越好)
tic                     %计算程序运行时间
t = 1; 

%=======程序终止条件(下列两个终止条件任选其一)======
 while ((t<=t_max))                             %达到最大迭代次数而终止
%========================= 
    %路径标识字符:标识每只蚂蚁的路径
    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=1+floor((K-2)*rand(1,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);
end
02-011M

4.完整MATLAB程序

matlab源码说明_我爱C编程的博客-CSDN博客

V

  • 4
    点赞
  • 20
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
%-- 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

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

打赏作者

我爱C编程

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

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

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

打赏作者

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

抵扣说明:

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

余额充值