使用灰狼算法求函数最小值

close all; 
clc; 
clear;
rand('seed',1);%保证每次的随机数结果一样
%------------------参数初始化------------------
wolf_num = 15;  % 狼群数量
max_iter = 150;  % 最大迭代次数
dim = 30;       % 变量维数
lb = -30*ones(1,dim);  % 自变量下边界
ub = 30*ones(1,dim);   % 自变量上边界
alpha_pos = zeros(1,dim); %"alpha" 狼只位置初始化
alpha_score = inf;        % "alpha"狼只适应度函数初始化
beta_pos = zeros(1,dim);  %  "beta" 狼只位置初始化
beta_score = inf;         %  "beta"狼只适应度函数初始化
delta_pos = zeros(1,dim); %  "delta"狼只位置初始化
delta_score = inf;        %  "delta"狼只适应度函数初始化
convergence_curve = zeros(1,max_iter);  % 每次迭代的最小适应度值初始化 
%------------------参数初始化------------------
count = 0;
position = init_pos(wolf_num,dim,ub,lb);  % 狼群位置初始化
while count < max_iter
    for i = 1:wolf_num
        flag_ub = position(i,:) > ub;%判断代表狼只位置的30维向量的第几维超出上边界，flag_ub为30维的布尔型，超出记为1，否则为0。
        flag_lb = position(i,:) < lb;%判断代表狼只位置的30维向量的第几维超出下边界，flag_lb为30维的布尔型，超出记为1，否则为0。
        %%%%%假设position(1：)=[40,-40,20,20,……],则flag_ub=[1,0,0,0……],flag_lb=[0,1,0,0……]
        position(i,:) = position(i,:).*(~(flag_ub+flag_lb))+flag_ub.*ub+flag_lb.*lb;  % 调整超出边界的狼群位置，将超出的部分限制在上边界和下边界
        fitness = f(position(i,:));   % 计算适应度值
        if fitness < alpha_score %更新alphe的适应度值
            alpha_score = fitness;
            alpha_pos = position(i,:);
        elseif fitness < beta_score%更新beta的适应度值
            beta_score = fitness;
            beta_pos = position(i,:);
        elseif fitness < delta_score%更新delta的适应度值
            delta_score = fitness;
            delta_pos = position(i,:);
        end
    end
    a = 2 - count*(2/max_iter);   % 更新a的值，a为线性
    for i = 1:wolf_num
        for j = 1:dim
            alpha = update_pos(alpha_pos(j),position(i,j),a);%更新alpha、beta、delta狼只位置
            beta = update_pos(beta_pos(j),position(i,j),a);
            delta = update_pos(delta_pos(j),position(i,j),a);
            position(i,j) = (alpha+beta+delta)/3;
        end
    end
    count = count + 1;
    convergence_curve(count) = alpha_score;
end
%-----------------------------绘图-------------------------------------
plot(1:max_iter,convergence_curve,'LineWidth',2,'LineStyle','-','Color','r');
xlabel('iteration'); ylabel('fitness'); title('GWO fitness curve');
grid on;
disp('The solution of GWO:'); disp(alpha_pos);
disp('The best fitness of GWO:'); disp(alpha_score);
%-----------------------------绘图-------------------------------------

%-----------------------------适应度函数-------------------------------
function res = f(x)
% f_min : 0
dim = 30; sum = 0;
for i = 1:dim-1
    sum = sum + 100*(x(i+1)-x(i)^2)^2+(x(i)-1)^2;
end
res = sum;
end
%-----------------------------适应度函数-------------------------------
%-----------------------------更新A、C和位置-------------------------------
function res = update_pos(v1,v2,a)
A = 2*a*rand()-a;
C = 2*rand();
temp = abs(C*v1-v2);
res = v1 - A*temp;
end
%-----------------------------更新A、C和位置-------------------------------
%-----------------------------位置初始化函数-------------------------------
function position = init_pos(wolf_num,dim,ub,lb)
position = zeros(wolf_num,dim);
for i = 1:wolf_num
    for j = 1:dim
        position(i,j) = rand( )*(ub(j)-lb(j)) + lb(j);
    end
end
end
%-----------------------------位置初始化函数-------------------------------