粒子群算法(PSO)的改进以及Matlab自带的PSO算法

时间:2021-08-08
本文章向大家介绍粒子群算法(PSO)的改进以及Matlab自带的PSO算法,主要包括粒子群算法(PSO)的改进以及Matlab自带的PSO算法使用实例、应用技巧、基本知识点总结和需要注意事项,具有一定的参考价值,需要的朋友可以参考一下。

粒子群算法的改进

线性递减惯性权重

我们的目的就是要在开始的时间内尽量的多在全局搜索,在后期的时间内不断缩短搜索范围,从而找到全局最大值。

线性递减惯性权重就是为了平衡全局搜索和局部搜索的地位。

%% 线性递减惯性权重的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示)
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比,使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形,继续在上面画图

%% 粒子群算法中的预设参数(参数的设置不是固定的,可以适当修改)
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子,也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子,也称为社会加速常数
w_start = 0.9;  % 初始惯性权重,通常取0.9
w_end = 0.4; % 粒子群最大迭代次数时的惯性权重,通常取0.4
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度(这里我们设置为[-vmax,vmax])

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群,计算每一个粒子的适应度
    fit(i) = Obj_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置(是一个n*narvs的向量)
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置(是一个1*narvs的向量)

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r');  % scatter3是绘制三维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代,一共迭代K次
    for i = 1:n   % 依次更新第i个粒子的速度与位置       
        w = w_start - (w_start - w_end) * d / K;  
        v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:));  % 更新第i个粒子的速度
        % 如果粒子的速度超过了最大速度限制,就对其进行调整
        for j = 1: narvs
            if v(i,j) < -vmax(j)
                v(i,j) = -vmax(j);
            elseif v(i,j) > vmax(j)
                v(i,j) = vmax(j);
            end
        end
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域,就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x(:,1);  % 更新散点图句柄的x轴的数据(此时粒子的位置在图上发生了变化)
    h.YData = x(:,2);   % 更新散点图句柄的y轴的数据(此时粒子的位置在图上发生了变化)
    h.ZData = fit;  % 更新散点图句柄的z轴的数据(此时粒子的位置在图上发生了变化)
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是:'); disp(gbest)
disp('此时最优值是:'); disp(Obj_fun2(gbest))


%% 三种递减惯性权重的对比
w_start = 0.9;  w_end = 0.4; 
K = 30;  d = 1:K; 
w1 = w_start-(w_start-w_end)*d/K;
w2 = w_start-(w_start-w_end)*(d/K).^2;
w3 = w_start-(w_start-w_end)*(2*d/K-(d/K).^2);
figure(3)
plot(d,w1,'r',d,w2,'b',d,w3,'g'); 
legend('w1','w2','w3')

非线性递减惯性权重

三种递减惯性权重对比

自适应惯性权重

与原来的相比,现在惯性权重和迭代次数与每个粒子的适应度有关

我们引入了适应度概念,适应度越小,说明距离最优解越近,我们需要局部搜索;适应度越大,说明我们距离最优解越远,我们需要全局搜索。

最大值问题与求解最小值问题相反。适应度越小,说明距离最优解越远,我们需要全局搜索;适应度越大,说明我们距离最优解越近,我们需要局部搜索。

%% 自适应权重的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示)
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比,使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形,继续在上面画图

%% 粒子群算法中的预设参数(参数的设置不是固定的,可以适当修改)
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子,也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子,也称为社会加速常数
w_max = 0.9;  % 最大惯性权重,通常取0.9
w_min = 0.4; % 最小惯性权重,通常取0.4
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度(这里我们设置为[-vmax,vmax])

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群,计算每一个粒子的适应度
    fit(i) = Obj_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置(是一个n*narvs的向量)
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置(是一个1*narvs的向量)

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r');  % scatter3是绘制三维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代,一共迭代K次
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        f_i = fit(i);  % 取出第i个粒子的适应度
        f_avg = sum(fit)/n;  % 计算此时适应度的平均值
        f_min = min(fit); % 计算此时适应度的最小值
        if f_i <= f_avg  
            if f_avg ~= f_min  % 如果分母为0,我们就干脆令w=w_min
                w = w_min + (w_max - w_min)*(f_i - f_min)/(f_avg - f_min);
            else
                w = w_min;
            end
        else
            w = w_max;
        end
        v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:));  % 更新第i个粒子的速度
        % 如果粒子的速度超过了最大速度限制,就对其进行调整
        for j = 1: narvs
            if v(i,j) < -vmax(j)
                v(i,j) = -vmax(j);
            elseif v(i,j) > vmax(j)
                v(i,j) = vmax(j);
            end
        end
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域,就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x(:,1);  % 更新散点图句柄的x轴的数据(此时粒子的位置在图上发生了变化)
    h.YData = x(:,2);   % 更新散点图句柄的y轴的数据(此时粒子的位置在图上发生了变化)
    h.ZData = fit;  % 更新散点图句柄的z轴的数据(此时粒子的位置在图上发生了变化)
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是:'); disp(gbest)
disp('此时最优值是:'); disp(Obj_fun2(gbest))

随机惯性权重

使用随机的惯性权重,可以避免在迭代前期局部搜索能力的不足;也可以避免在迭代后期全局搜索能力的不足。

红色方框部分应是标准差,一般取0.2~0.5之间的数

%% 随机权重的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示)
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比,使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形,继续在上面画图

%% 粒子群算法中的预设参数(参数的设置不是固定的,可以适当修改)
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2;  % 每个粒子的个体学习因子,也称为个体加速常数
c2 = 2;  % 每个粒子的社会学习因子,也称为社会加速常数
w_max = 0.9;  % 最大惯性权重,通常取0.9
w_min = 0.4; % 最小惯性权重,通常取0.4
sigma = 0.3; % 正态分布的随机扰动项的标准差(一般取0.2-0.5之间)
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度(这里我们设置为[-vmax,vmax])

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群,计算每一个粒子的适应度
    fit(i) = Obj_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置(是一个n*narvs的向量)
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置(是一个1*narvs的向量)

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r');  % scatter3是绘制三维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代,一共迭代K次
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        w = w_min + (w_max - w_min)*rand(1) + sigma*randn(1);
        v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:));  % 更新第i个粒子的速度
        % 如果粒子的速度超过了最大速度限制,就对其进行调整
        for j = 1: narvs
            if v(i,j) < -vmax(j)
                v(i,j) = -vmax(j);
            elseif v(i,j) > vmax(j)
                v(i,j) = vmax(j);
            end
        end
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域,就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x(:,1);  % 更新散点图句柄的x轴的数据(此时粒子的位置在图上发生了变化)
    h.YData = x(:,2);   % 更新散点图句柄的y轴的数据(此时粒子的位置在图上发生了变化)
    h.ZData = fit;  % 更新散点图句柄的z轴的数据(此时粒子的位置在图上发生了变化)
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是:'); disp(gbest)
disp('此时最优值是:'); disp(Obj_fun2(gbest))

压缩(收缩)因子

讨论完前面的速度权重之后,我们来讨论个体学习因子和社会学系因子前的系数\(c_1\)\(c_2\)

压缩因子法中应用较多的个体学习因子\(c_1\)和社会学习因子\(c_2\)均取2.05,用我们自己的符号可以表示为:

%% 带有压缩因子(收缩因子)的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示)
clear; clc

%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比,使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形,继续在上面画图

%% 粒子群算法中的预设参数(参数的设置不是固定的,可以适当修改)
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2.05;  % 每个粒子的个体学习因子,也称为个体加速常数
c2 = 2.05;  % 每个粒子的社会学习因子,也称为社会加速常数
C = c1+c2; 
fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子
w = 0.9;  % 惯性权重 
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度(这里我们设置为[-vmax,vmax])

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群,计算每一个粒子的适应度
    fit(i) = Obj_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置(是一个n*narvs的向量)
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置(是一个1*narvs的向量)

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r');  % scatter3是绘制三维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代,一共迭代K次
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)));  % 更新第i个粒子的速度
        %  有时候会看到直接写成:0.657*v(i,:) +  1.496*rand(1)*(pbest(i,:) - x(i,:)) + 1.496*rand(1)*(gbest - x(i,:)));
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域,就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x(:,1);  % 更新散点图句柄的x轴的数据(此时粒子的位置在图上发生了变化)
    h.YData = x(:,2);   % 更新散点图句柄的y轴的数据(此时粒子的位置在图上发生了变化)
    h.ZData = fit;  % 更新散点图句柄的z轴的数据(此时粒子的位置在图上发生了变化)
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是:'); disp(gbest)
disp('此时最优值是:'); disp(Obj_fun2(gbest))

非对称学习因子

简而言之,相对于之前的对称学习因子的\(c_1\)\(c_2\)的固定取值,非对称学习因子引入了公式,使\(c_1\)逐渐减小,\(c_2\)逐渐增大

%% 非对称学习因子的粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(动画演示)
clear; clc
% 毛开富, 包广清, 徐驰. 基于非对称学习因子调节的粒子群优化算法[J]. 计算机工程, 2010(19):188-190.(公式10中的减号应该改为加号)


%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1');  ylabel('x2');  zlabel('y');  % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比,使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on  % 不关闭图形,继续在上面画图

%% 粒子群算法中的预设参数(参数的设置不是固定的,可以适当修改)
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1_ini = 2.5;  % 个体学习因子的初始值
c1_fin = 0.5;  % 个体学习因子的最终值
c2_ini = 1;  % 社会学习因子的初始值
c2_fin = 2.25;  % 社会学习因子的最终值
w = 0.9;  % 惯性权重
K = 100;  % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度(这里我们设置为[-vmax,vmax])

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群,计算每一个粒子的适应度
    fit(i) = Obj_fun2(x(i,:));   % 调用Obj_fun2函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置(是一个n*narvs的向量)
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置(是一个1*narvs的向量)

%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r');  % scatter3是绘制三维散点图的函数(这里返回h是为了得到图形的句柄,未来我们对其位置进行更新)

%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
for d = 1:K  % 开始迭代,一共迭代K次
    %    w = 0.9 - 0.5*d/K;  % 效果不是很好的话可以同时考虑线性递减权重
        c1 = c1_ini + (c1_fin - c1_ini)*d/K;
        c2 = c2_ini + (c2_fin - c2_ini)*d/K;
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:));  % 更新第i个粒子的速度
        % 如果粒子的速度超过了最大速度限制,就对其进行调整
        for j = 1: narvs
            if v(i,j) < -vmax(j)
                v(i,j) = -vmax(j);
            elseif v(i,j) > vmax(j)
                v(i,j) = vmax(j);
            end
        end
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域,就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun2(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun2(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun2(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun2(gbest);  % 更新第d次迭代得到的最佳的适应度
    pause(0.1)  % 暂停0.1s
    h.XData = x(:,1);  % 更新散点图句柄的x轴的数据(此时粒子的位置在图上发生了变化)
    h.YData = x(:,2);   % 更新散点图句柄的y轴的数据(此时粒子的位置在图上发生了变化)
    h.ZData = fit;  % 更新散点图句柄的z轴的数据(此时粒子的位置在图上发生了变化)
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是:'); disp(gbest)
disp('此时最优值是:'); disp(Obj_fun2(gbest))

自动退出迭代循环

当粒子已经找到最佳位置后,再增加迭代次数只会浪费计算时间,那么我们能否设计一个策略,能够自动退出迭代呢?

我们给出如下思路供大家研究:

%% 自适应权重且带有收缩因子的粒子群算法PSO: 求解四种不同测试函数的最小值(动画演示)
clear; clc

%% 粒子群算法中的预设参数(参数的设置不是固定的,可以适当修改)
tic % 开始计时
n = 1000; % 粒子数量
narvs = 30; % 变量个数
c1 = 2.05;  % 每个粒子的个体学习因子,也称为个体加速常数
c2 = 2.05;  % 每个粒子的社会学习因子,也称为社会加速常数
C = c1+c2; 
fai = 2/abs((2-C-sqrt(C^2-4*C))); % 收缩因子
w_max = 0.9;  % 最大惯性权重,通常取0.9
w_min = 0.4; % 最小惯性权重,通常取0.4
K = 1000;  % 迭代的次数
% Sphere函数
% vmax = 30*ones(1,30); % 粒子的最大速度
% x_lb = -100*ones(1,30); % x的下界
% x_ub = 100*ones(1,30); % x的上界
% Rosenbrock函数
vmax = 10*ones(1,30); % 粒子的最大速度
x_lb = -30*ones(1,30); % x的下界
x_ub = 30*ones(1,30); % x的上界
% Rastrigin函数
% vmax = 1.5*ones(1,30); % 粒子的最大速度
% x_lb = -5.12*ones(1,30); % x的下界
% x_ub = 5.12*ones(1,30); % x的上界
% Griewank函数
% vmax = 150*ones(1,30); % 粒子的最大速度
% x_lb = -600*ones(1,30); % x的下界
% x_ub = 600*ones(1,30); % x的上界

%% 改进:自动退出迭代循环
Count = 0; % 计数器初始化为0
max_Count = 30;  % 最大计数值初始化为30
tolerance = 1e-6;  % 函数变化量容忍度,取10^(-6)

%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
    x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1);    % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs);  % 随机初始化粒子的速度(这里我们设置为[-vmax,vmax])

%% 计算适应度
fit = zeros(n,1);  % 初始化这n个粒子的适应度全为0
for i = 1:n  % 循环整个粒子群,计算每一个粒子的适应度
    fit(i) = Obj_fun3(x(i,:));   % 调用Obj_fun3函数来计算适应度
end 
pbest = x;   % 初始化这n个粒子迄今为止找到的最佳位置(是一个n*narvs的向量)
ind = find(fit == min(fit), 1);  % 找到适应度最小的那个粒子的下标
gbest = x(ind,:);  % 定义所有粒子迄今为止找到的最佳位置(是一个1*narvs的向量)


%% 迭代K次来更新速度与位置
% fitnessbest = ones(K,1);  % 初始化每次迭代得到的最佳的适应度
% 注意我把上面这行注释掉了,因为我们不知道最后一共要迭代多少次,可能后面会自动跳出循环
% 初始化的目的一般是节省运行的时间,在这里这个初始化的步骤我们可以跳过
for d = 1:K  % 开始迭代,一共迭代K次
    tem = gbest; % 将上一步找到的最佳位置保存为临时变量
    for i = 1:n   % 依次更新第i个粒子的速度与位置
        f_i = fit(i);  % 取出第i个粒子的适应度
        f_avg = sum(fit)/n;  % 计算此时适应度的平均值
        f_min = min(fit); % 计算此时适应度的最小值
        if f_i <= f_avg  
            w = w_min + (w_max - w_min)*(f_i - f_min)/(f_avg - f_min);
        else
            w = w_max;
        end
        v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)));  % 更新第i个粒子的速度
        x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
        % 如果粒子的位置超出了定义域,就对其进行调整
        for j = 1: narvs
            if x(i,j) < x_lb(j)
                x(i,j) = x_lb(j);
            elseif x(i,j) > x_ub(j)
                x(i,j) = x_ub(j);
            end
        end
        fit(i) = Obj_fun3(x(i,:));  % 重新计算第i个粒子的适应度
        if fit(i) < Obj_fun3(pbest(i,:))   % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
           pbest(i,:) = x(i,:);   % 那就更新第i个粒子迄今为止找到的最佳位置
        end
        if  fit(i) < Obj_fun3(gbest)  % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
            gbest = pbest(i,:);   % 那就更新所有粒子迄今为止找到的最佳位置
        end
    end
    fitnessbest(d) = Obj_fun3(gbest);  % 更新第d次迭代得到的最佳的适应度
    delta_fit = abs(Obj_fun3(gbest) - Obj_fun3(tem));   % 计算相邻两次迭代适应度的变化量
    if delta_fit < tolerance  % 判断这个变化量和“函数变化量容忍度”的相对大小,如果前者小,则计数器加1
        Count = Count + 1;
    else
        Count = 0;  % 否则计数器清0
    end   
    if Count > max_Count  % 如果计数器的值达到了最大计数值
        break;  % 跳出循环
    end
end

figure(2) 
plot(fitnessbest)  % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是:'); disp(gbest)
disp('此时最优值是:'); disp(Obj_fun3(gbest))

toc % 结束计时

Matlab自带的粒子群函数

我们之前已经实现过了粒子群函数,也都对粒子群函数做了大量改进。实际上,matlab2014后的版本已经提供给了我们现成的粒子群函数。我们之前所有的改进都是基于最原始的粒子群算法公式提出的,但是Matlab自身的粒子群函数是截然不同的实现方式。

Matlab中particleswarm函数采用的是自适应的邻域模式

Q:什么是自适应呢?

A:自适应体现在: 如果适应度开始停滞时,粒子群搜索会从邻域模式向全局模式转换,一旦适应度开始下降,则又恢复到邻域模式,以免陷入局部最优。 当适应度的停滞次数足够大时,惯性系数开始逐渐变小,从而利于局部搜索

Q:什么是邻域模式呢?

A:首先我们来了解一下全局模式,我们之前的粒子群算法是基于全局模式的,全局模式就是将周围所有的粒子均看成邻域粒子(就是邻居),是邻域模式的极端情况。会有较快的收敛速度,不过容易陷入局部最优。而邻域模式只将周围的几个粒子视为邻域粒子,有利于在多个区域进行搜索,避免进入局部最优。

预设参数选取

函数内部的实现

这里不要求掌握,仅了解即可:

自动退出迭代循环

第一个方法需要设置Maxstallterations和functiontolerance两个参数,功能是当迭代大于Maxstallterations时就开始检查funchange(最优值减去当前最优值),当funchange小于functiontolerance时就自动退出循环。

第二个方法直接设置最大循环次数,当大于循环次数时就退出循环。

其余的方法不常用,就不多介绍了。

如果是第一种方式退出迭代的话,我们可以将粒子群算法得到的解作为初始值,继续调用其他的函数来进行混合求解,例如我们熟悉的fmincon函数

代码实现

matlab代码如下:

%% Matlab自带的粒子群函数 particleswarm
% particleswarm函数是求最小值的
% 如果目标函数是求最大值则需要添加负号从而转换为求最小值。
clear;clc
%  Matlab中粒子群算法函数的参考文献
%       Mezura-Montes, E., and C. A. Coello Coello. "Constraint-handling in nature-inspired numerical optimization: Past, present and future." Swarm and Evolutionary Computation. 2011, pp. 173–194.
%       Pedersen, M. E. "Good Parameters for Particle Swarm Optimization." Luxembourg: Hvass Laboratories, 2010.
%       Iadevaia, S., Lu, Y., Morales, F. C., Mills, G. B. & Ram, P. T. Identification of optimal drug combinations targeting cellular networks: integrating phospho-proteomics and computational network analysis. Cancer Res. 70, 6704-6714 (2010).
%       Liu, Mingshou , D. Shin , and H. I. Kang . "Parameter estimation in dynamic biochemical systems based on adaptive Particle Swarm Optimization." Information, Communications and Signal Processing, 2009. ICICS 2009. 7th International Conference on IEEE Press, 2010.

%% 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值(最小值为8)
narvs = 2; % 变量个数
x_lb = [-15 -15]; % x的下界(长度等于变量的个数,每个变量对应一个下界约束)
x_ub = [15 15]; % x的上界
[x,fval,exitflag,output] = particleswarm(@Obj_fun2, narvs, x_lb, x_ub)   

%% 直接调用particleswarm函数进行求解测试函数
narvs = 30; % 变量个数
% Sphere函数
% x_lb = -100*ones(1,30); % x的下界
% x_ub = 100*ones(1,30); % x的上界
% Rosenbrock函数
x_lb = -30*ones(1,30); % x的下界
x_ub = 30*ones(1,30); % x的上界
% Rastrigin函数
% x_lb = -5.12*ones(1,30); % x的下界
% x_ub = 5.12*ones(1,30); % x的上界
% Griewank函数
% x_lb = -600*ones(1,30); % x的下界
% x_ub = 600*ones(1,30); % x的上界
[x,fval,exitflag,output] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub)  


%% 绘制最佳的函数值随迭代次数的变化图
options = optimoptions('particleswarm','PlotFcn','pswplotbestf')   
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 展示函数的迭代过程
options = optimoptions('particleswarm','Display','iter');
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 修改粒子数量,默认的是:min(100,10*nvars)
options = optimoptions('particleswarm','SwarmSize',50);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 在粒子群算法结束后继续调用其他函数进行混合求解(hybrid  n.混合物合成物; adj.混合的; 杂种的;) 
options = optimoptions('particleswarm','HybridFcn',@fmincon);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 惯性权重的变化范围,默认的是0.1-1.1
options = optimoptions('particleswarm','InertiaRange',[0.2 1.2]);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 个体学习因子,默认的是1.49(压缩因子)
options = optimoptions('particleswarm','SelfAdjustmentWeight',2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 社会学习因子,默认的是1.49(压缩因子)
options = optimoptions('particleswarm','SocialAdjustmentWeight',2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 最大的迭代次数,默认的是200*nvars
options = optimoptions('particleswarm','MaxIterations',10000);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 领域内粒子的比例 MinNeighborsFraction,默认是0.25 
options = optimoptions('particleswarm','MinNeighborsFraction',0.2);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 函数容忍度FunctionTolerance, 默认1e-6, 用于控制自动退出迭代的参数
options = optimoptions('particleswarm','FunctionTolerance',1e-8);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 最大停滞迭代数MaxStallIterations, 默认20, 用于控制自动退出迭代的参数
options = optimoptions('particleswarm','MaxStallIterations',50);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)

%% 不考虑计算时间,同时修改三个控制迭代退出的参数
tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',100,'MaxIterations',100000);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)
toc

%% 在粒子群结束后调用其他函数进行混合求解
tic
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',50,'MaxIterations',20000,'HybridFcn',@fmincon);
[x,fval] = particleswarm(@Obj_fun3,narvs,x_lb,x_ub,options)
toc

求解的函数要单独放在新的函数文件中,在主函数文件(上)中用函数句柄的形式引用!

Matlab输出停止搜索的原因


原文地址:https://www.cnblogs.com/Cassiopeia/p/15114510.html