使用Matlab群智能优化算法求解Shubert函数最优解

%% 定义Shubert函数 function y = shubert(x) n = size(x,2); sum1 = 0; sum2 = 0; for i=1:n for j=1:5 sum1 = sum1 + j*cos((j+1)x(i)+j); sum2 = sum2 + jcos((j+1)x(i)+j); end end y = sum1sum2; end

%% 绘制Shubert函数图形 x = -10:0.1:10; y = -10:0.1:10; [X,Y] = meshgrid(x,y); Z = zeros(size(X)); for i=1:size(X,1) for j=1:size(X,2) Z(i,j) = shubert([X(i,j),Y(i,j)]); end end figure;surf(X,Y,Z);

%% 绘制Shubert函数的二维图形 x = -10:0.1:10; y = shubert([x;zeros(size(x))]); figure;plot(x,y);

%% 使用PSO算法求解Shubert函数最优解 nVar = 2; % 两个变量：x和y VarSize = [1 nVar]; % 每个变量的大小 VarMin = -10; % 变量下限 VarMax = 10; % 变量上限

% PSO参数 MaxIt = 100; % 最大迭代次数 nPop = 50; % 群体大小 w = 1; % 惯性权重 wdamp = 0.99; % 惯性权重衰减系数 c1 = 1.5; % 个体学习因子 c2 = 2; % 群体学习因子

% 初始化粒子群 empty_particle.Position = []; empty_particle.Velocity = []; empty_particle.Cost = []; empty_particle.Best.Position = []; empty_particle.Best.Cost = [];

particle = repmat(empty_particle, nPop, 1);

% 初始化群体最优解 GlobalBest.Cost = inf;

for i=1:nPop % 初始化位置 particle(i).Position = unifrnd(VarMin, VarMax, VarSize); % 初始化速度 particle(i).Velocity = zeros(VarSize); % 计算成本 particle(i).Cost = shubert(particle(i).Position); % 更新个体最优解 particle(i).Best.Position = particle(i).Position; particle(i).Best.Cost = particle(i).Cost; % 更新群体最优解 if particle(i).Best.Cost < GlobalBest.Cost GlobalBest = particle(i).Best; end end

% 初始化历史最优解 BestCosts = zeros(MaxIt, 1);

% 主循环 for it=1:MaxIt for i=1:nPop % 更新速度 particle(i).Velocity = wparticle(i).Velocity ... + c1rand(VarSize).(particle(i).Best.Position - particle(i).Position) ... + c2rand(VarSize).*(GlobalBest.Position - particle(i).Position); % 限制速度范围 particle(i).Velocity = max(particle(i).Velocity, VarMin); particle(i).Velocity = min(particle(i).Velocity, VarMax); % 更新位置 particle(i).Position = particle(i).Position + particle(i).Velocity; % 限制位置范围 particle(i).Position = max(particle(i).Position, VarMin); particle(i).Position = min(particle(i).Position, VarMax); % 计算成本 particle(i).Cost = shubert(particle(i).Position); % 更新个体最优解 if particle(i).Cost < particle(i).Best.Cost particle(i).Best.Position = particle(i).Position; particle(i).Best.Cost = particle(i).Cost; % 更新群体最优解 if particle(i).Best.Cost < GlobalBest.Cost GlobalBest = particle(i).Best; end end end % 保存历史最优解 BestCosts(it) = GlobalBest.Cost; % 输出迭代信息 disp(['迭代次数 ' num2str(it) '，最优解为 ' num2str(GlobalBest.Cost)]); % 更新惯性权重 w = w * wdamp; end

% 绘制成本函数历史记录 figure;plot(BestCosts);

% 输出最优解 disp(['最优解：(' num2str(GlobalBest.Position(1)) ',' num2str(GlobalBest.Position(2)) '), f=' num2str(GlobalBest.Cost)]);

%% 计算Shubert函数的18个最优解 % 利用全局最优解和局部最优解 % 局部最优解通过在全局最优解周围随机生成新的解并进行优化得到

% 计算全局最优解 x0 = GlobalBest.Position; f0 = GlobalBest.Cost;

% 计算局部最优解 nLocal = 5; % 生成5个局部最优解 LocalBest = repmat(empty_particle, nLocal, 1); for i=1:nLocal % 在全局最优解的周围随机生成新的解 x = GlobalBest.Position + unifrnd(-1, 1, size(GlobalBest.Position)); % 优化新的解 [x, f] = particleswarm(@shubert, nVar, VarMin, VarMax, optimoptions('particleswarm', 'Display', 'off'), x); % 保存局部最优解 LocalBest(i).Position = x; LocalBest(i).Cost = f; end

% 将全局最优解和局部最优解合并 AllBest = [GlobalBest; LocalBest];

% 对所有最优解按照f值进行排序 [~, SortOrder] = sort([AllBest.Cost]); AllBest = AllBest(SortOrder);

% 输出18个最优解 disp('18个最优解：'); for i=1:18 disp(['(' num2str(AllBest(i).Position(1)) ',' num2str(AllBest(i).Position(2)) '), f=' num2str(AllBest(i).Cost)]); end

%% 绘制18个最优解的位置 figure;surf(X,Y,Z); hold on; for i=1:18 plot3(AllBest(i).Position(1), AllBest(i).Position(2), AllBest(i).Cost, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); end hold off;

%% 绘制18个最优解的二维图形 x = -10:0.1:10; y = shubert([x;zeros(size(x))]); figure;plot(x,y); hold on; for i=1:18 x = AllBest(i).Position(1); y = AllBest(i).Cost; plot(x, y, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); end hold off;

%% 绘制18个最优解的三维图形 figure;plot3(AllBest.Cost, AllBest.Position(:,1), AllBest.Position(:,2), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); xlabel('f'); ylabel('x'); zlabel('y');