Ｍatlab 非线性有约束规划的粒子群算法「建议收藏」

clear;close;clc
%% 约束条件和目标函数构建
fun = @(x) x(1)^2 + x(2)^2 + x(3)^2 + 8;
bind1 = @(x) x(1)^2 - x(2) + x(3)^2 >= 0;
bind2 = @(x) x(1) + x(2)^2 + x(3)^2 <= 20;
% 不太适合等式约束
ekc = 1e-10;
bind3 = @(x) abs(-x(1) - x(2)^2 + 2) <= ekc;
bind4 = @(x) abs(x(2) + 2*x(3)^2 - 3) <= ekc;

%% 初始化
popsize = 500; % 粒子个数
dim = 3; % 维度
max_iter = 100; % 最大迭代次数
xlimit_max = [2 3 20]'; % 由等式约束推出位置边界
xlimit_min = zeros(dim,1); 
vlimit_max = 1*ones(dim,1);
vlimit_min = -1*ones(dim,1);
w = 0.6; % 惯性权重
c1 = 0.5;c2 = 1.5;
pr = 0.4; % 变异率
pop_x = zeros(dim,popsize);  % 当前粒子位置
pop_v = zeros(dim,popsize); % 当前粒子速度
fitness_pop = zeros(1,popsize); % 粒子群当前位置适应度函数
fitness_lbest = zeros(1,popsize); % 个体粒子的历史最优极值
rand('state',sum(clock));

for j = 1:popsize 
    % 位置初始化
    pop_x(1,j) = xlimit_min(1) + rand*(xlimit_max(1) - xlimit_min(1));
    pop_x(2,j) = sqrt(2-pop_x(1,j));
    pop_x(3,j) = sqrt((3 - pop_x(2,j))/2);
    % 速度初始化
    for  i = 1:dim
        pop_v(i,j) = vlimit_min(i) + rand*(vlimit_max(i) - vlimit_min(i));
    end
end
%% 初始化个体极值
lbest = pop_x; % 个体历史最佳极值记录
for j =1: popsize 
    if bind1(pop_x(:,j))
        if bind2(pop_x(:,j))
            fitness_lbest(j) = fun(pop_x(:,j));
        else fitness_lbest(j) = 500;
        end
    else fitness_lbest(j) = 500;
    end
end

%% 初始化全局极值
popbest = pop_x(:,1);
fitness_popbest = fitness_lbest(1);
for j = 2:popsize 
    if fitness_lbest(j) < fitness_popbest
        fitness_popbest = fitness_lbest(j);
        popbest = pop_x(:,j);
    end
end
tic
%% 粒子群迭代
iter = 1; % 当前迭代次数
record = zeros(max_iter,1); % 记录每次迭代的全局极小值
format long;
while iter <= max_iter
    for j = 1:popsize 
        % 更新速度 边界处理
        pop_v(:,j) = w*pop_v(:,j) + c1*rand*(lbest(:,j) - pop_x(:,j)) +...
            c2*rand*(popbest - pop_x(:,j));
        for i = 1:dim 
            if pop_v(i,j) > vlimit_max(i)
                pop_v(i,j) = vlimit_max(i);
            elseif pop_v(i,j) < vlimit_min(i) 
                pop_v(i,j) = vlimit_min(i);
            end
        end
        % 更新位置 边界处理 修正位置 （等式约束）
        pop_x(:,j) = pop_x(:,j) + pop_v(:,j);
        for i = 1:dim 
            if pop_x(i,j) > xlimit_max(i)
                pop_x(i,j)  = xlimit_max(i);
            elseif pop_x(i,j) < xlimit_min(i)
                pop_x(i,j) = xlimit_min(i);
            end
        end
        
        % 进行自适应变异
        if rand < pr 
            i = ceil(dim*rand);
            pop_x(i,j) = xlimit_min(i) + rand*(xlimit_max(i) - xlimit_min(i));
        end
        % 约束条件限制 类似罚函数法
        if bind1(pop_x(:,j))
            if bind2(pop_x(:,j))
                if bind3(pop_x(:,j))
                    if bind4(pop_x(:,j))
                        fitness_pop(j) = fun(pop_x(:,j));
                    else fitness_pop(j) = 500;
                    end
                else fitness_pop(j) = 500;
                end
            else fitness_pop(j) = 500;
            end
        else fitness_pop(j) = 500;
        end
        % 当前适应度与个体历史最佳适应度作比较
        if fitness_pop(j) < fitness_lbest(j)
            lbest(:,j) = pop_x(:,j);
            fitness_lbest(j) = fitness_pop(j);
        end
        % 个体历史最佳适应度与种群历史最佳适应度作比较
        if fitness_popbest > fitness_lbest(j)
            fitness_popbest = fitness_lbest(j);
            popbest = lbest(:,j);
        end
    end
    record(iter) = fitness_popbest;
    iter = iter + 1;
    
end
toc
%% 输出解
minx = popbest
miny = fitness_popbest
plot(record,'r-');
title('粒子群算法迭代过程');
xlabel('迭代次数');
ylabel('当前迭代最佳函数值');