【Matlab代码】蚁群算法解得最小值

clc;
clear;
f = inline('20 + x.^2 + y.^2 -10*(cos(2*pi*x) + cos(2*pi*y))');
x = -5:0.01:5;
y = -5:0.01:5;
[X,Y] = meshgrid(x,y);
F = f(X,Y);
figure(1);
mesh(X,Y,F);
xlabel('横坐标x'); ylabel('纵坐标y'); zlabel('空间坐标z');
hold on;
lower_x = -5;
upper_x = 5;
lower_y = -5;
upper_y = 5; 
ant = 1000;    
times = 300;   
rou = 0.9;     
p0 = 0.2;      
ant_x = zeros(1,ant);  
ant_y = zeros(1,ant);  
tau = zeros(1,ant);    
Macro = zeros(1,ant); 
for i=1:ant
    ant_x(i) = (upper_x-lower_x)*rand() + lower_x;
    ant_y(i) = (upper_y-lower_y)*rand() + lower_y;
    tau(i) = f(ant_x(i),ant_y(i));          
    % plot3(ant_x(i),ant_y(i),tau(i),'k*'); 
    hold on;
    Macro = zeros(1,ant);
end
fprintf('蚁群搜索开始(找最小值,绿色):\n');
T = 1;
trap = 0; 
tau_best = zeros(1,times);  
p = zeros(1,ant);   
while T < times
    lamda = 1/T;
    [tau_best(T),bestindex] = min(tau);
    plot3(ant_x(bestindex), ant_y(bestindex), f(ant_x(bestindex),ant_y(bestindex)), 'r*');
    hold on;
    for i = 1:ant
        p(i) = (tau(bestindex) - tau(i))/tau(bestindex); 
    end
    for i = 1:ant
        if p(i) < p0
            tempx = ant_x(i) + 0.5*(2*rand-1)*lamda;  
            tempy = ant_y(i) + 0.5*(2*rand-1)*lamda;
        else
            tempx = ant_x(i) + (upper_x-lower_x)*(rand-2.5);  % 5/2=2.5
            tempy = ant_y(i) + (upper_y-lower_y)*(rand-2.5);
        end
        if tempx < lower_x
            tempx = lower_x;
        end
        if tempx > upper_x
            tempx = upper_x;
        end
        if tempy < lower_y
            tempy = lower_y;
        end
        if tempy > upper_y
            tempy = upper_y;
        end
        if f(tempx,tempy) < tau(i)
            ant_x(i) = tempx;
            ant_y(i) = tempy;
            Macro(i) = f(tempx,tempy);
        end
    end
    for i = 1:ant
        tau(i) = (1-rou)*tau(i) + Macro(i);
    end
    T = T + 1;
    if T >= times
        fprintf('蚁群搜索到的最小值点:(%.5f,%.5f,%.5f)\n',...
        ant_x(bestindex), ant_y(bestindex), f(ant_x(bestindex),ant_y(bestindex)));
        fprintf('搜索次数:%d\n',T)
        if f(ant_x(bestindex),ant_y(bestindex)) > 0.8
            trap = trap + 1;  % 进入陷阱次数+1
            fprintf('发生大概率事件: 陷入局部极小值,自动再次执行程序!\n\n')
            T = 1;  
            tau_best = zeros(1,times);
            p = zeros(1,ant);
            for i=1:ant
                ant_x(i) = (upper_x-lower_x)*rand() + lower_x;
                ant_y(i) = (upper_y-lower_y)*rand() + lower_y;
                tau(i) = f(ant_x(i),ant_y(i));        
                Macro = zeros(1,ant);
            end
        end
    end
end
hold on;
trap = trap + 1;  
syms x y;
f = 20 + x^2 + y^2 -10*(cos(2*pi*x) + cos(2*pi*y));
fx = diff(f,x);
fy = diff(f,y);
acc = 0.0001;          
study_step = 0.001;    
x = ant_x(bestindex); 
y = ant_y(bestindex);
k = 0; 
fprintf('走出局部极值,梯度下降精确搜索开始(黄线):\n');
while eval(fx)~=0 | eval(fy)~=0 
 ans_tmp = [x,y] - study_step*[eval(fx),eval(fy)];
 acc_tmp = sqrt((ans_tmp(1)-x)^2 + (ans_tmp(2)-y)^2);
 if acc_tmp <= acc
 fprintf('精确极值坐标为:(%.5f,%.5f,%.5f)\n',ans_tmp(1),ans_tmp(2),f_tmp);
        fprintf('迭代次数:%d\n\n',k);
        plot3(ans_tmp(1),ans_tmp(2),f_tmp,'y.');
        hold off
 break;
 end
 x = ans_tmp(1);
 y = ans_tmp(2);
 f_tmp = eval(f);
    plot3(x,y,f_tmp,'y.')
    hold on;
    k = k + 1;  
end
if trap > 0
    fprintf('本模型走出陷阱概率:%.1f%%\n',1/trap*100);
end