Matlab：4维、单目标、约束、粒子群优化算法

% 主调用函数（求最大值）
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