
% 主调用函数(求最大值)
clc;
clear;
close all;
% 初始化种群
N = 100; % 初始种群个数
D = 4; % 空间维数
iter = 50; % 迭代次数
x_limit = [10, 13; 14, 25; 3, 8; 0.3, 0.5]; % 位置限制
v_limit = [-1, 1; -1, 1; -0.1, 0.1; -0.1, 0.1]; % 速度限制
x = zeros(N, D);
for i = 1:100
x(i,1) = randi([10,13],1,1);%初始种群的位置
x(i,2) = randi([14,25],1,1);%初始种群的位置
x(i,3) = x_limit(3, 1) + (x_limit(3, 2) - x_limit(3, 1)) * rand();%初始种群的位置
x(i,4) = x_limit(4, 1) + (x_limit(4, 2) - x_limit(4, 1)) * rand();%初始种群的位置
end
v(:,1) = randi([-1,1],N,1); % 初始种群z1方向的速度
v(:,2) = randi([-1,1],N,1); % 初始种群z2方向的速度
v(:,3) = rands(N, 1) * 1; % 初始种群3方向的速度
v(:,4) = rands(N, 1) * 1; % 初始种群phi_R方向的速度
p_best = x; % (初始化)每个个体的历史最佳位置
f_best = zeros(1, D); % (初始化)种群的历史最佳位置
fp_best = zeros(N, 1) - inf; % (初始化)每个个体的历史最佳适应度为负无穷 -
fg_best = -inf; % (初始化)种群历史最佳适应度为负无穷 -
w = 1; % 惯性权重
c1 = 1; % 自我学习因子
c2 = 1; % 群体学习因子
i = 1;
record = zeros(iter, 1); % 记录器
while i <= iter
fx = f_xy(x(:,1), x(:,2), x(:,3), x(:,4)); % 个体当前适应度
for j = 1:N
if fp_best(j) < fx(j) % 记录最大值
fp_best(j) = fx(j); % 更新个体历史最佳适应度
p_best(j,:) = x(j,:); % 更新个体历史最佳位置
end
end
if fg_best < max(fp_best)
[fg_best, ind_max] = max(fp_best); % 更新群体历史最佳适应度
f_best = p_best(ind_max, :); % 更新群体历史最佳位置
end
v = v * w + c1 * randi([0,1],1,1) * (p_best - x) + c2 * randi([0,1],1,1) * (repmat(f_best, N, 1) - x); % 速度更新
% 速度处理
for t=1:N
for k=1:D
if v(t,k) > v_limit(k,2) % 超速处理
v(t,k) = v_limit(k,2);
elseif v(t,k) < v_limit(k,1) % 慢速处理
v(t,k) = v_limit(k,1);
end
end
end
x = x + v; % 位置更新
% 边界处理
for t=1:N
for k=1:D
if x(t,k) > x_limit(k,2) % 超过边界上限
x(t,k) = x_limit(k,2);
elseif x(t,k) < x_limit(k,1) % 超过边界下限
x(t,k) = x_limit(k,1);
end
end
end
record(i) = fg_best; % 最大值记录
i = i + 1;
if mod(i,10) == 0
i % 收敛进度输出
end
pause(0.1)
end
figure(1)
plot(-record);
xlabel('迭代次数');
ylabel('适应度值');
title('收敛过程');
disp(['最小值:',num2str(-fg_best)]);
disp(['变量取值:z1=',num2str(f_best(1)),' z2=',num2str(f_best(2)),...
' m=',num2str(f_best(3)),' phi_R=',num2str(f_best(4))]);
function f = f_xy(z1, z2, m, phi_R)
% 适应度函数
n = 2;
b = phi_R.*m.*sqrt(z1.^2+z2.^2)./2;
theta1 = atan(z1./z2);
V1 = b.*pi.*cos(theta1).*(3.*m.^2.*z1.^2-6.*m.*z1.*b.*sin(theta1)+4.*b.^2.*(sin(theta1)).^2)./12;
V2 = b.*pi.*sin(theta1).*(3.*m.^2.*z1.^2-6.*m.*z1.*b.*cos(theta1)+4.*b.^2.*(cos(theta1)).^2)./12;
f = -(n.*V1+2.*V2);
end
参考:https://blog.csdn.net/Wang_Dou_Dou_/article/details/119417078
浙公网安备 33010602011771号