粒子群算法求一元函数最值的MATLAB实现-一一网

在上一篇详细介绍粒子群算法实现分组背包的随笔中，已经详细介绍了粒子群算法的主要思想，如果掌握了用粒子群算法如何实现分组背包的话，那么要将其修改成一元函数求最值的应用简直易如反掌。这里如下先copy一份之前总结的用粒子群算法实现分组背包大致思想：

随机产生了一堆粒子，每个粒子代表背包的一种情况（选了哪3个物品），初始粒子全都是局部最优粒子。
计算每个粒子的适应度值（也就是每个粒子代表的背包的价值、体积、重量），然后每个粒子适应度按优劣选出非劣粒子。
之后进入迭代过程，每次迭代中，从非劣粒子中随机选一个作为全局最优粒子，按照公式计算粒子速度（跟较优粒子和全局最优粒子息息相关），使每个粒子产生移动（也就是有概率替换代表的物品）。
如果移动后的粒子比之前的粒子更优则替代原来的局部最优粒子
然后把非劣粒子和局部最优粒子放一起，再按优劣选出非劣粒子
去除重复的非劣粒子，进入下一次迭代。

现在针对求一元函数最小值对上述思想做出如下改动：

随机产生了一堆粒子，每个粒子代表一个横坐标，初始粒子全都是局部最优粒子。
计算每个粒子的适应度值（也就是每个粒子代表的纵坐标值），然后每个粒子适应度按优劣选出初始非劣粒子。
之后进入迭代过程，每次迭代中，从非劣粒子中随机选一个作为全局最优粒子，按照公式计算粒子速度（跟较优粒子和全局最优粒子息息相关），使每个粒子产生移动（也就是横坐标左移或右移）。
如果移动后的粒子比之前的粒子更优则替代原来的局部最优粒子
然后把非劣粒子和局部最优粒子放一起，再按优劣选出非劣粒子
去除重复的非劣粒子，进入下一次迭代。

一：一元函数最值问题

假设有如下函数：

$y(x) = sin(x^2) – 2sin(x) + x * sin(x)$

画出该函数的图像如图所示：

我们要如何找到函数在[-3, 9]之间的最小值呢？首先放一个展示粒子群如何找到最小值的动态图：

从图中可以看到，我们先生成了50个粒子，用蓝色小点表示，用红色小点表示非劣解粒子（最后的最小值点一定是包括在非劣解粒子中的），但这里可以看到非劣解粒子一直只有一个，所以它就是我们所求的最小值点。可以看到每次迭代粒子的位置都会变化，并不断向最小值位置处移动。

下面具体介绍下粒子群算法如何求一元函数最值。还是按照上一篇随笔一样将完整的MATLAB代码分解成七部分讲解。

二：粒子群算法求一元函数最值问题

2.1 输入参数、固定参数初始化

注意参数中的惯性权值可以自己修改成自己觉得合适的值，它会影响粒子每次运动的步长。wmax越大迭代初期粒子运动步长越大；wmin越小迭代末期粒子运动步长越小。

clear, clc, close all;

%% 输入参数、固定参数初始化
xsize = 50;         % 粒子数
ITER = 50;         % 迭代次数
c1 = 0.8; c2 = 0.8;      % 常数
wmax = 1.0; wmin = 0.01;  % 惯性权重相关常数
v = zeros(1, xsize); % 粒子速度初始化
复制代码

2.2 粒子群位置、适应度、最佳位置、最佳适应度初始化

随机产生粒子群 $x$ ，表示每个粒子的横坐标，注意想要产生 $[a, b]$ 之间的均匀分布随机数的话，就用 $a + (b – a) * rand(1,1)$ 这个表达式。

适应度其实就是用粒子群的纵坐标。

%% 粒子群位置、适应度、最佳位置、最佳适应度初始化
x = -3 + 12 * rand(1, xsize); % 随机粒子群位置生成（表示横坐标[-3, 9]）

% 粒子群适应度
y = zeros(1, xsize); % 粒子群纵坐标
for i = 1 : xsize
    y(i) = sin(x(i).^2) - 2*sin(x(i)) + x(i) * sin(x(i));
end

bestx = x; % 粒子群位置最佳值
besty = y; % 粒子群最佳适应度
复制代码

2.3 初始筛选非劣解

第一次筛选非劣解，之后每次迭代都会重新筛选一次。其中的判断条件很重要，可以根据问题的限制而改变。这里就是判断每个粒子是否比别的所有粒子都更符合要求（纵坐标更小）。

%% 初始筛选非劣解
cnt = 1;
for i = 1 : xsize
    fljflag = 1;
    for j = 1 : xsize
        if j ~= i
            if y(i) > y(j) % i为劣解
                fljflag = 0;
                break;
            end
        end
    end 
    if fljflag == 1
        flj(cnt) = y(i); % 非劣解适应度
        fljx(cnt) = x(i); % 非劣解位置
        cnt = cnt + 1;
    end
end
复制代码