zoukankan      html  css  js  c++  java
  • 蚁群算法

    蚁群算法解决TSP问题、二次分配问题、背包问题,是蚁群算法的经典应用。从mathworks上下载了三个代码,看了注释,对蚁群算法的有了更全面的了解。

    % Project Title: Ant Colony Optimization for Traveling Salesman Problem
    clc;
    clear;
    close all;
    
    % Problem Definition
    model=CreateModel();
    CostFunction=@(tour) TourLength(tour,model);
    nVar=model.n;
    
    % ACO Parameters
    MaxIt=300;      % Maximum Number of Iterations 最大迭代次数
    nAnt=40;        % Number of Ants (Population Size)  种群数量
    Q=1;
    tau0=10*Q/(nVar*mean(model.D(:)));	% Initial Phromone
    
    alpha=1;        % Phromone Exponential Weight
    beta=1;         % Heuristic Exponential Weight
    
    rho=0.05;       % Evaporation Rate蒸发率(信息素损失率)
    
    % Initialization
    eta=1./model.D;             % Heuristic Information Matrix  启发信息矩阵
    tau=tau0*ones(nVar,nVar);   % Phromone Matrix
    BestCost=zeros(MaxIt,1);    % Array to Hold Best Cost Values   保存最好的值的数组
    
    % Empty Ant
    empty_ant.Tour=[];
    empty_ant.Cost=[];
    
    % Ant Colony Matrix
    ant=repmat(empty_ant,nAnt,1);
    
    % Best Ant
    BestSol.Cost=inf;
    
    % ACO Main Loop
    for it=1:MaxIt    
        % Move Ants
        for k=1:nAnt        
            ant(k).Tour=randi([1 nVar]);        %在[1,nVar]中均匀随机产生一个数
            for l=2:nVar            
                i=ant(k).Tour(end);            
                P=tau(i,:).^alpha.*eta(i,:).^beta;            
                P(ant(k).Tour)=0;            
                P=P/sum(P);            
                j=RouletteWheelSelection(P);    %轮盘赌        
                ant(k).Tour=[ant(k).Tour j];            
            end
            
            ant(k).Cost=CostFunction(ant(k).Tour);        
            if ant(k).Cost<BestSol.Cost
                BestSol=ant(k);
            end        
        end
        
        % Update Phromones 更新信息素
        for k=1:nAnt        
            tour=ant(k).Tour;        
            tour=[tour tour(1)]; %#ok        
            for l=1:nVar            
                i=tour(l);
                j=tour(l+1);            
                tau(i,j)=tau(i,j)+Q/ant(k).Cost;     %可以用其他的信息素更新规则
                
            end        
        end
        
        % Evaporation
        tau=(1-rho)*tau;
        
        % Store Best Cost
        BestCost(it)=BestSol.Cost;
        
        % Show Iteration Information
        disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
        
        % Plot Solution
        figure(1);
        PlotSolution(BestSol.Tour,model);
        pause(0.01);    
    end
    
    % Results
    figure;
    plot(BestCost,'LineWidth',2);
    xlabel('Iteration');
    ylabel('Best Cost');
    grid on;
    
    % Project Title: Ant Colony Optimization for Traveling Salesman Problem
    function model=CreateModel()
        x=[82 91 12 92 63 9  28 55 96 97 15 98 96 49 80 14 42 92 80 96];    
        y=[66 3  85 94 68 76 75 39 66 17 71 3  27 4  9  83 70 32 95 3 ];    
        n=numel(x);     %元素个数
        D=zeros(n,n);   %计算邻接矩阵    
        for i=1:n-1
            for j=i+1:n            
                D(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);            
                D(j,i)=D(i,j);            
            end
        end
        
        model.n=n;
        model.x=x;
        model.y=y;
        model.D=D;
    end
    
    function PlotSolution(tour,model) 
        tour=[tour tour(1)];    
        plot(model.x(tour),model.y(tour),'k-o',...
            'MarkerSize',10,...
            'MarkerFaceColor','y',...
            'LineWidth',1.5);
        
        xlabel('x');
        ylabel('y');
        
        axis equal;
        grid on;
        
    	alpha = 0.1;
    	
        xmin = min(model.x);
        xmax = max(model.x);
        dx = xmax - xmin;
        xmin = floor((xmin - alpha*dx)/10)*10;
        xmax = ceil((xmax + alpha*dx)/10)*10;
        xlim([xmin xmax]);
        
        ymin = min(model.y);
        ymax = max(model.y);
        dy = ymax - ymin;
        ymin = floor((ymin - alpha*dy)/10)*10;
        ymax = ceil((ymax + alpha*dy)/10)*10;
        ylim([ymin ymax]);    
    end
    
    %轮盘赌选择
    function j=RouletteWheelSelection(P)
        r=rand;    
        C=cumsum(P);    
        j=find(r<=C,1,'first');
    end
    
    %当前解的长度
    function L=TourLength(tour,model)
        n=numel(tour);   %tour里面是TSP的遍历路径
        tour=[tour tour(1)];    %围成圈
        L=0;
        for i=1:n
            L=L+model.D(tour(i),tour(i+1));   %累加计算总长度
        end
    end
    

      

    % Project Title: Ant Colony Optimization for Quadratic Assignment Problem
    clc;
    clear;
    close all;
    
    % Problem Definition
    model=CreateModel();
    CostFunction=@(p) MyCost(p,model);
    nVar=model.n;
    
    % ACO Parameters
    MaxIt=500;      % Maximum Number of Iterations   迭代次数
    nAnt=50;        % Number of Ants (Population Size) 种群数量
    Q=1;
    tau0=10;        % Initial Phromone
    alpha=0.3;      % Phromone Exponential Weight--蚂蚁学习历史经验的权重
    rho=0.1;        % Evaporation Rate蒸发率
    
    % Initialization
    tau=tau0*ones(model.m,nVar);
    BestCost=zeros(MaxIt,1);    % Array to Hold Best Cost Values
    
    % Empty Ant
    empty_ant.Tour=[];
    empty_ant.Cost=[];
    
    % Ant Colony Matrix
    ant=repmat(empty_ant,nAnt,1);
    
    % Best Ant
    BestSol.Cost=inf;
    
    % ACO Main Loop
    for it=1:MaxIt    
        % Move Ants
        for k=1:nAnt        
            ant(k).Tour=[];        
            for l=1:nVar            
                P=tau(:,l).^alpha;            
                P(ant(k).Tour)=0;            
                P=P/sum(P);            
                j=RouletteWheelSelection(P);            
                ant(k).Tour=[ant(k).Tour j];            
            end
            ant(k).Cost=CostFunction(ant(k).Tour);  
            %蚂蚁方案的代价只取决于方案本身,不像TSP还取决于邻接边的长短
            
            if ant(k).Cost<BestSol.Cost
                BestSol=ant(k);
            end        
        end
        
        % Update Phromones
        for k=1:nAnt        
            tour=ant(k).Tour;        
            for l=1:nVar            
                tau(tour(l),l)=tau(tour(l),l)+Q/ant(k).Cost;            
            end        
        end
        
        % Evaporation
        tau=(1-rho)*tau;
        
        % Store Best Cost
        BestCost(it)=BestSol.Cost;
        
        % Show Iteration Information
        disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
    
        % Plot Solution
        figure(1);
        PlotSolution(BestSol.Tour,model);
        pause(0.01);
        
    end
    
    % Results
    figure;
    plot(BestCost,'LineWidth',2);
    xlabel('Iteration');
    ylabel('Best Cost');
    grid on;
    
    
    function model=CreateModel()
        x=[67 80 62 34 54 36 53 46 39 35 83 58 87 90 83 38 26 78 49 67];    
        y=[9  81 9  43 89 30 95 87 1  74 85 86 56 86 22 73 36 34 17 37];    
        m=numel(x);    
        d=zeros(m,m);
        for p=1:m-1
            for q=p+1:m
                d(p,q)=sqrt((x(p)-x(q))^2+(y(p)-y(q))^2);
                d(q,p)=d(p,q);
            end
        end
        
        w=[0    6    6    3    5    5    5
           6    0    6    4  -10    3    6
           6    6    0    4    5    8    6
           3    4    4    0    4    4  100
           5  -10    5    4    0    3    4
           5    3    8    4    3    0    4
           5    6    6  100    4    4    0];
    
        n=size(w,1);
        
        model.n=n;
        model.m=m;
        model.w=w;  %权重
        model.x=x;
        model.y=y;
        model.d=d;  %邻接矩阵
    end
    
    function z=MyCost(p,model)
        n=model.n;
        w=model.w;
        d=model.d;
        z=0;
        
        for i=1:n-1
            for j=i+1:n            
                z=z+w(i,j)*d(p(i),p(j));            
            end
        end
    end
    
    function PlotSolution(p,model)
        x=model.x;
        y=model.y;    
        plot(x,y,'o','MarkerSize',20,'Color',[0.4 0.5 0.9]);    
        hold on;
        
        plot(x(p),y(p),'sk','MarkerSize',16,'MarkerFaceColor','y');
        
        n=model.n;
        for i=1:n
            text(x(p(i)),y(p(i)),num2str(i),...
                'HorizontalAlignment','center',...
                'VerticalAlignment','middle');
        end
        
        title('Quadratic Assignment Problem');
        
        hold off;
        axis equal;
        grid on;
        
        alpha = 0.1;
        
        xmin = min(x);
        xmax = max(x);
        dx = xmax - xmin;
        xmin = floor((xmin - alpha*dx)/10)*10;
        xmax = ceil((xmax + alpha*dx)/10)*10;
        xlim([xmin xmax]);
        
        ymin = min(y);
        ymax = max(y);
        dy = ymax - ymin;
        ymin = floor((ymin - alpha*dy)/10)*10;
        ymax = ceil((ymax + alpha*dy)/10)*10;
        ylim([ymin ymax]);
    end
    
    function j=RouletteWheelSelection(P)
        r=rand;    
        C=cumsum(P);    
        j=find(r<=C,1,'first');
    end
    
    % Project Title: Ant Colony Optimization for Binary Knapsack Problem
    clc;
    clear;
    close all;
    % Problem Definition
    model=CreateModel();
    CostFunction=@(x) MyCost(x,model);
    nVar=model.n;
    
    % ACO Parameters
    MaxIt=300;      % Maximum Number of Iterations
    nAnt=40;        % Number of Ants (Population Size)
    Q=1;
    tau0=0.1;        % Initial Phromone
    alpha=1;        % Phromone Exponential Weight
    beta=0.02;      % Heuristic Exponential Weight启发式权重指数
    
    rho=0.1;        % Evaporation Rate蒸发率
    % Initialization
    N=[0 1];
    
    eta=[model.w./model.v
         model.v./model.w];           % Heuristic Information
    
    tau=tau0*ones(2,nVar);      % Phromone Matrix
    BestCost=zeros(MaxIt,1);    % Array to Hold Best Cost Values
    
    % Empty Ant
    empty_ant.Tour=[];
    empty_ant.x=[];
    empty_ant.Cost=[];
    empty_ant.Sol=[];
    
    % Ant Colony Matrix
    ant=repmat(empty_ant,nAnt,1);
    
    % Best Ant
    BestSol.Cost=inf;
    
    % ACO Main Loop
    for it=1:MaxIt    
        % Move Ants
        for k=1:nAnt        
            ant(k).Tour=[];        
            for l=1:nVar            
                P=tau(:,l).^alpha.*eta(:,l).^beta;            
                P=P/sum(P);            
                j=RouletteWheelSelection(P);            
                ant(k).Tour=[ant(k).Tour j];            
            end        
            ant(k).x=N(ant(k).Tour);        
            [ant(k).Cost, ant(k).Sol]=CostFunction(ant(k).x);        
            if ant(k).Cost<BestSol.Cost
                BestSol=ant(k);
            end        
        end
        
        % Update Phromones
        for k=1:nAnt        
            tour=ant(k).Tour;        
            for l=1:nVar
                tau(tour(l),l)=tau(tour(l),l)+Q/ant(k).Cost;            
            end        
        end
        
        % Evaporation
        tau=(1-rho)*tau;
        
        % Store Best Cost
        BestCost(it)=BestSol.Cost;
        
        % Show Iteration Information
        if BestSol.Sol.IsFeasible
            FeasiblityFlag = '*';
        else
            FeasiblityFlag = '';
        end
        disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it)) ' ' FeasiblityFlag]);
        
    end
    
    % Results
    figure;
    plot(BestCost,'LineWidth',2);
    xlabel('Iteration');
    ylabel('Best Cost');
    grid on;
    
    function model=CreateModel()
        v=[391 444 250 330 246 400 150 266 268 293 471 388 364 493 202 161 410 270 384 486];    
        w=[55  52  59  24  52  46  45  34  34  59  59  28  57  21  47  66  64  42  22  23];    
        n=numel(v);    
        W=500;
        
        model.n=n;
        model.v=v;
        model.w=w;
        model.W=W;
    end
    
    function [z, sol]=MyCost(x,model)
        v=model.v;
        w=model.w;
        W=model.W;
    
        V1=sum(v.*x);
        W1=sum(w.*x);
        V0=sum(v.*(1-x));
        W0=sum(w.*(1-x));
        
        Violation=max(W1/W-1,0);    
        z=V0*(1+100*Violation);
        
        sol.V1=V1;
        sol.W1=W1;
        sol.V0=V0;
        sol.W0=W0;
        sol.Violation=Violation;
        sol.z=z;
        sol.IsFeasible=(Violation==0);
    end
    
    function j=RouletteWheelSelection(P)
        r=rand;    
        C=cumsum(P);    
        j=find(r<=C,1,'first');
    end
    

      

  • 相关阅读:
    cuda npp库旋转图片
    Xml序列化 详解
    jsonp简介
    在centos7下安装.net core
    安装vs2017后造成无法打开xproj项目无法打开
    SqlServer 语法
    js自定义事件
    HttpWebResponse 解压gzip、deflate压缩
    centos7 安装.net core的方法
    帮助类-从tfs获取数据
  • 原文地址:https://www.cnblogs.com/liugl7/p/7856546.html
Copyright © 2011-2022 走看看