Matlab粒子群算法搜索费马点-求到多个点距离之和最小的点

0. 前言

本文是在我之前一篇博客的基础上进行了扩展和延申，原文使用了matlab自带的优化函数，而本文采用了一个比较经典的随机优化算法——粒子群算法，对给定的目标函数进行求解。

Matlab实现搜索费马点-求到多个点距离之和最小的点_ 一只博客-CSDN博客_到多个点距离之和最小的点https://blog.csdn.net/qq_42276781/article/details/118356037

1. 粒子群算法介绍

粒子群优化算法(Particle Swarm optimization,PSO)又翻译为粒子群算法、微粒群算法、或微粒群优化算法。是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。它可以被纳入多主体优化系统(Multiagent Optimization System, MAOS)。粒子群优化算法是由Eberhart博士和kennedy博士发明。                                                                                                     ——摘自百度百科

 基本步骤如下：

Step1  初始化所有粒子的位置和速度，计算其适应度，并以适应度最佳的粒子作为最佳个体，将当前这批粒子作为最佳群体。

Step2  更新各个粒子的速度，再以各个粒子的速度更新其位置，重新计算每个粒子的适应度。

Step3  在新生成的这批粒子中：如果有粒子的适应度优于最佳个体，则更新最佳个体；如果第个粒子的适应度优于最佳群体中第个粒子的适应度，更新最佳群体中第个粒子的适应度。

Step4  重复Step2和Step3，直到满足优化目标或迭代次数达到上限。

2. 粒子群算法搜索费马点

2.1 适应度函数

粒子群算法中的适应度函数就是目标函数，即对于费马点，其到其余点的距离之和最小

该适应度函数写作matlab的表达式如下，对该代码有疑问的可以参考前言中提到的博客。

%% I. 导入数据，设置适应度函数
% 数据
data = [1 8;3 7;6 4;7 5;9 2;5 8];
% 适应度函数，点x到所有点的距离之和
fun = @(x) sum(sqrt(sum((data-x).^2, 2)));
% 如果出现矩阵维度报错，请使用下方的fun
%fun = @(x) sum(sqrt(sum((data-repmat(x,length(data),1)).^2, 2)));

2.2 参数配置

设置学习因子和，最大进化次数，粒子群规模， 中横纵坐标的最大最小值，和方向的最大速度和最小速度。

%% II. 离散粒子群算法 参数配置 
c1 = 1.49445;        %个体学习因子 
c2 = 1.49445;        %群体学习因子 
evolution_max = 200; %进化次数
pop_size = 1000;     %种群规模
%最大最小种群位置，设置最大最小速度
x_max = max(data(:,1));
x_min = min(data(:,1));
y_max = max(data(:,2));
y_min = min(data(:,2));
v_x_max = x_max-x_min;
v_x_min = -v_x_max;
v_y_max = y_max-y_min;
v_y_min = -v_y_max;

2.3 初始化粒子

%% III. 产生初始粒子和速度,并计算其适应度
pop = zeros(pop_size, 2);
v = zeros(pop_size, 2);
fitness = zeros(1, pop_size);
pop(:,1) = x_min + (x_max-x_min) * rand(pop_size,1);
pop(:,2) = y_min + (y_max-y_min) * rand(pop_size,1);
v(:,1) = v_x_min + (v_x_max-v_x_min) * rand(pop_size,1);
v(:,2) = v_y_min + (v_y_max-v_y_min) * rand(pop_size,1);
for i = 1: pop_size
    fitness(i) = fun(pop(i,:));
end

是一个pop_size*2的矩阵，表示一群粒子的位置，第行表示第个粒子的位置，第一列表示粒子在方向的位置，第二列表示粒子在方向的位置，如表示第个粒子的横坐标。

是一个pop_zize*2的矩阵，表示一群粒子的速度，第行表示第个粒子的速度，第一列表示粒子在方向的速度，第二列表示粒子在方向的速度，如表示第个粒子在方向的速度。

是一维行向量，表示pop_size个粒子的适应度，如表示第个粒子的适应度。

这里使用了一个生成指定范围随机数的语句，用以生成区间[a,b]内的随机浮点数，返回形式是m*n的矩阵，里面的元素不完全相同。

a + (b-a) * rand(m, n)

2.4 计算最佳个体和最佳群体

%% IV. 初始化最佳个体和最佳群体
[fitness_best, index_best] = min(fitness); % 获取最佳适应度及其索引
I = pop(index_best,:);                     % 初始化最佳个体
P = pop;                                   % 初始化最佳群体
I_fitness_best = fitness_best;             % 初始化最佳个体适应度
P_fitness_best = fitness;                  % 初始化最佳群体适应度

在2.3节中我们计算了所有粒子的适应度，现在需要从中选出最优的适应度和其索引（它是的第几个粒子）。matlab提供了min函数，可以搜索一个向量中的最小值及其位置。

接着把有最优适应度的粒子的位置保存到，把当前的这群粒子暂时作为最佳群体，保存到，最后分别保存和的适应度和，用以之后的对比。

2.5 迭代寻优

2.5.1 速度更新

在Step2中，粒子需要更新其速度，第个粒子速度需要向最佳个体和最佳群体中的第个粒子学习。更新后其在方向的速度可表示为

表示第个粒子方向的速度，和分别表示第个粒子向最佳个体和最佳群体中的第个粒子的学习因子，是一个常数，和表示两个随机生成的0-1之间的随机数，表示当前所有粒子，表示当前第个粒子在方向的位置，表示最佳个体在方向的位置，表示最佳群体中的第个粒子在方向的位置。同理，更新后其在方向的速度可表示为

 

表示第个粒子方向的速度，和分别表示第个粒子向最佳个体和最佳群体中的第个粒子的学习因子，是一个常数，和表示两个随机生成的0-1之间的随机数(不同的下标用来示意这四个随机数并不一定完全相同），表示当前所有粒子，表示当前第个粒子在方向的位置，表示最佳个体在方向的位置，表示最佳群体中的第个粒子在方向的位置。

为了提高效率，可以一次性更新所有粒子在方向和方向的速度

这里的符号意义如下

            

           

         

 matlab相应的代码如下

%更新速度
v(:,1) = v(:,1) + c1*rand*(I(1)-pop(:,1)) + c2*rand*(P(:,1)-pop(:,1));
v(:,2) = v(:,2) + c1*rand*(I(2)-pop(:,2)) + c2*rand*(P(:,2)-pop(:,2));

2.5.2 位置更新

每一次迭代所需的时间为单位1，则

这里，

更新后的粒子可能会超出预设的边界，需要对超出边界的粒子进行修正。

如图，对于费马点的搜索，费马点的横坐标必定不会小于中的任何点的横坐标，必定不会大于任何点的横坐标，即

纵坐标必定不会小于中的任何点的纵坐标，必定不会大于任何点的纵坐标，即

 找出中横纵坐标的最大最小值，组合后得到4个点，将其连接得到上述矩形。不难发现，如果费马点在红色矩形框外，记距离最近的矩形框上的点为，显然到6个蓝色点的距离之和小于到6个蓝色点的距离之和。

故我们只需要使用构成这个矩形框的和对超出边界的粒子进行修正即可。单个粒子的修正逻辑如下

 

 

同样，仿照2.5.1节，我们也同时更新所有粒子的位置

         

更新所有粒子位置以及修正越界粒子位置的代码如下

pop(:,1) = pop(:,1) + v(:,1);
pop(:,2) = pop(:,2) + v(:,2);
%调整越界的粒子
pop(pop(:,1) > x_max, 1) = x_max;
pop(pop(:,1) < x_min, 1) = x_min;
pop(pop(:,2) > y_max, 2) = y_max;
pop(pop(:,2) < y_min, 2) = y_min;

在修正越界粒子的代码中，嵌套了matlab的索引代码，如果有不理解之处可以阅读我的这篇博客，或在评论区讨论。

matlab核心知识点-索引index运用_ 一只博客-CSDN博客_matlab中index索引https://blog.csdn.net/qq_42276781/article/details/121715378

2.5.3 迭代逻辑

%% V. 迭代寻优
fitness_evolution = zeros(evolution_max, 1);
position_evolution = zeros(evolution_max, 2);
for i = 1: evolution_max
    %记录最佳个体
    fitness_evolution(i) = I_fitness_best;
    position_evolution(i,:) = I;
    %更新速度
    v(:,1) = v(:,1) + c1*rand*(I(1)-pop(:,1)) + c2*rand*(P(:,1)-pop(:,1));
	v(:,2) = v(:,2) + c1*rand*(I(2)-pop(:,2)) + c2*rand*(P(:,2)-pop(:,2));
    %更新位置
    pop(:,1) = pop(:,1) + v(:,1);
    pop(:,2) = pop(:,2) + v(:,2);
    %调整越界的粒子
    pop(pop(:,1) > x_max, 1) = x_max;
    pop(pop(:,1) < x_min, 1) = x_min;
    pop(pop(:,2) > y_max, 2) = y_max;
    pop(pop(:,2) < y_min, 2) = y_min;
    %更新适应度
    for j = 1: pop_size
        fitness(j) = fun(pop(j,:));
    end
    %更新最佳群体
    index = fitness < P_fitness_best;
    P_fitness_best(index) = fitness(index);
    P(index,:) = pop(index,:);
    %更新最佳个体
    [fitness_best, index_best] = min(fitness);
    if fitness_best < I_fitness_best
        I_fitness_best = fitness_best; %更新最佳个体适应度
        I = pop(index_best,:);    %更新最佳个体
    end
end

2.6 输出结果

%% VI.输出结果
figure, plot(fitness_evolution), title('最优个体适应度')
xlabel('进化代数'), ylabel('适应度');
fprintf('最佳费马点为(%.2f, %.2f)\n',I(1), I(2))
fprintf('最小距离和为%.2f\n', I_fitness_best)

绘制每一次迭代的最优适应度，输出最后求解得到的费马点和最小距离和

2.7 模拟动画

通过动画来模拟每一次迭代得到的费马点，即费马点的搜索过程。因为本文使用的数据量较少，费马点很快就搜索到了，第一次迭代的费马点和最终的费马点相距较近（第一次迭代得到的最短距离之和为18.736左右，最终的距离之和为18.724左右），难以发现明显的移动。如果点较多，就可以发现明显移动。

figure
for i = 1: evolution_max
    plot(data(:,1), data(:,2), 'bo');
    hold on
    plot(position_evolution(i,1), position_evolution(i,2), 'rp');
    hold off
    m(:,i) =getframe;
end

3. 完整代码

clear, close all
%% I. 导入数据，设置适应度函数
%数据
data = [1 8;3 7;6 4;7 5;9 2;5 8];
%适应度函数，点x到所有点的距离之和
fun = @(x) sum(sqrt(sum((data-x).^2, 2)));
%如果出现矩阵维度报错，请使用下方的fun
%fun = @(x) sum(sqrt(sum((data-repmat(x,length(data),1)).^2, 2)));
%% II. 离散粒子群算法 参数配置 
c1 = 1.49445;        %个体学习因子 
c2 = 1.49445;        %群体学习因子 
evolution_max = 200; %进化次数
pop_size = 1000;     %种群规模
%最大最小种群位置，设置最大最小速度
x_max = max(data(:,1));
x_min = min(data(:,1));
y_max = max(data(:,2));
y_min = min(data(:,2));
v_x_max = x_max-x_min;
v_x_min = -v_x_max;
v_y_max = y_max-y_min;
v_y_min = -v_y_max;
%% III. 产生初始粒子和速度,并计算其适应度
pop = zeros(pop_size, 2);
v = zeros(pop_size, 2);
fitness = zeros(1, pop_size);
pop(:,1) = x_min + (x_max-x_min) * rand(pop_size,1);
pop(:,2) = y_min + (y_max-y_min) * rand(pop_size,1);
v(:,1) = v_x_min + (v_x_max-v_x_min) * rand(pop_size,1);
v(:,2) = v_y_min + (v_y_max-v_y_min) * rand(pop_size,1);
for i = 1: pop_size
    fitness(i) = fun(pop(i,:));
end
%% IV. 初始化最佳个体和最佳群体
[fitness_best, index_best] = min(fitness); %获取最佳适应度及其索引
I = pop(index_best,:);                     %初始化最佳个体
P = pop;                                   %初始化最佳群体
I_fitness_best = fitness_best;             %初始化最佳个体适应度
P_fitness_best = fitness;                  %初始化最佳群体适应度
%% V. 迭代寻优
fitness_evolution = zeros(evolution_max, 1);
position_evolution = zeros(evolution_max, 2);
for i = 1: evolution_max
    %记录最佳个体
    fitness_evolution(i) = I_fitness_best;
    position_evolution(i,:) = I;
    %更新速度
    v(:,1) = v(:,1) + c1*rand*(I(1)-pop(:,1)) + c2*rand*(P(:,1)-pop(:,1));
	v(:,2) = v(:,2) + c1*rand*(I(2)-pop(:,2)) + c2*rand*(P(:,2)-pop(:,2));
    %更新位置
    pop(:,1) = pop(:,1) + v(:,1);
    pop(:,2) = pop(:,2) + v(:,2);
    %调整越界的粒子
    pop(pop(:,1) > x_max, 1) = x_max;
    pop(pop(:,1) < x_min, 1) = x_min;
    pop(pop(:,2) > y_max, 2) = y_max;
    pop(pop(:,2) < y_min, 2) = y_min;
    %更新适应度
    for j = 1: pop_size
        fitness(j) = fun(pop(j,:));
    end
    %更新最佳群体
    index = fitness < P_fitness_best;
    P_fitness_best(index) = fitness(index);
    P(index,:) = pop(index,:);
    %更新最佳个体
    [fitness_best, index_best] = min(fitness);
    if fitness_best < I_fitness_best
        I_fitness_best = fitness_best; %更新最佳个体适应度
        I = pop(index_best,:);         %更新最佳个体
    end
end
%% VI.输出结果
figure, plot(fitness_evolution), title('最优个体适应度')
xlabel('进化代数'), ylabel('适应度');
fprintf('最佳费马点为(%.2f, %.2f)\n',I(1), I(2))
fprintf('最小距离和为%.2f\n', I_fitness_best)
%% VII.模拟动画
figure
for i = 1: evolution_max
    plot(data(:,1), data(:,2), 'bo');
    hold on
    plot(position_evolution(i,1), position_evolution(i,2), 'rp');
    hold off
    m(:,i) =getframe;
end

 

 有问题可以进群交流讨论