Matlab实现搜索费马点-求到多个点距离之和最小的点

目录

1. 初始化

2. 设置目标函数

step1 data-x

step2 (data-x).^2

step3 sum((data-x).^2, 2)

step4 sqrt(sum((data-x).^2, 2))

step5 sum(sqrt(sum((data-x).^2, 2)))

3. 初值x0的选取 

4. 求解最优点

5. 输出数据和画图

6. 完整代码和推广

7. 交流讨论

---

费马问题(Fermat problem)是著名的几何极值问题。费马(Fermat , P. de)曾提出一问题征解：“已知一个三角形，求作一点，使其与这个三角形的三个顶点的距离之和为极小。”

费马点即到三个点距离之和最小的点，推而广之，到n个点距离之和最小的点是不是也可以叫做费马点呢？

这个费马点非常具有现实意义，比如：

某物流公司准备为n个城市修建一个物流中心，为了效益最大化，要求物流中心到这n个城市的距离之和最小（物流中心不局限在这n个城市中）。

为简化问题，将n设置为6，随机生成6个点的坐标位置如下：

p1(1, 8)       p2(3, 7)       p3(6, 4)
p4(7, 5)       p5(9, 2)       p6(5, 8)

本文使用matlab自带的函数fminsearch来搜索费马点 

fminsearch函数使用无导数法来计算无约束的多变量函数的最小值，非常适合费马点的搜索。

当然，还可以使用遗传算法、模拟退火算法、粒子群算法等随机优化算法进行求解，现给出粒子群算法求解费马点的方法如下Matlab粒子群算法搜索费马点-求到多个点距离之和最小的点_ 一只博客-CSDN博客https://blog.csdn.net/qq_42276781/article/details/121699607

1. 初始化

%% 初始化
data = [1 8;3 7;6 4;7 5;9 2;5 8];
n = length(data);
distance = zeros(1, n);

使用length函数获取data行数和列数中最大的那个，即输入的点的数量，生成一个一行n列的一维向量，用于储存n个点到其他点的距离。

2. 设置目标函数

%% 目标函数
fun = @(x) sum(sqrt(sum((data-x).^2, 2)));
% 如果出现矩阵维度报错，请使用下方的fun
%fun = @(x) sum(sqrt(sum((data-repmat(x,n,1)).^2, 2)));

这里用到了匿名函数

下面简单举例介绍一下匿名函数的用法，这里的x不是横坐标的意思，x可以是一维的，也可以是二维的。f([4 7]) = ([4 7]-[2 4]).^2 = [2 3].^2 = [4 9]

下面来解读目标函数的运作逻辑，设这里的x表示的一维向量为[m n]，其中m表示点的横坐标，n表示点的纵坐标。

step1 data-x

版本较低的matlab不支持该语法，会提示“矩阵维度必须一致”。

可以将 data-x替换为 data-repmat(x,n,1)，repmat函数可以将变量向多维扩充，下面的示例是将横向量[1 2]向纵向重复4次，横向重复2次。

所以repmat(x,n,1)就是将x向纵向重复n次（n为data的长度），横向重复1次，即

注意这里repmat(x,n,1)里的n是6，因为data的长度为6

step2 (data-x).^2

需要注意，matlab中 A.*B 和 A*B是完全不同的两个概念

A.*B计算时会将A中的每个元素分别乘以B中的元素，如

而A*B计算时服从矩阵乘法运算，要求矩阵A的列数必须等于矩阵B的行数。

step3 sum((data-x).^2, 2)

这里简单介绍一下sum函数：

sum(A,2)会累加矩阵A的每一行，得到一个列向量；

sum(A,1)会累加矩阵A的每一列，得到一个行向量。

step4 sqrt(sum((data-x).^2, 2))

step5 sum(sqrt(sum((data-x).^2, 2)))

这里也可以用sum(sqrt(sum((data-x).^2, 2)),1)

因为一维向量是一种特殊的矩阵，所以：

对于一维列向量A来说，sum(A)等价于sum(A,1)

对于一维行向量A来说，sum(A)等价于sum(A,2)

3. 初值x0的选取 

fminsearch函数的输入参数有两个，分别是目标函数fun和变量初值x0。

变量初值x0并不是说是横坐标的值，它可以是二维的，三维的甚至n维的，本文中的x0是一个二维变量，第一项是横坐标初值，第二项是纵坐标初值。

初值的选取也比较讲究，选的好可以加快搜索速度节约时间。

一个很自然的想法就是在这6个点中寻找一个点P，要求这个点到其他5个点的距离之和最小（到它自己的距离是0，即到6个点的距离之和最小）

最后搜索得到的费马点Q必然落在P附近，这样就能提高搜索效率，避免盲目搜索。

初值x0的搜索代码如下

%% 在已知点中所搜费马点
for i = 1:n
    distance(i) = fun(data(i,:));
end     
[~,index] = min(distance);

这里我们只需要获取距离和最小的点的编号即可（通过编号来获取该点的横纵坐标），不需要知道距离和的具体数值。

[~, index] = min(distance)返回最小距离的索引index, 具体用法见min

4. 求解最优点

%% 将该点作为初始值，使用fminsearch搜索
x0 = data(index,:);
[best_point, min_dis] = fminsearch(fun, x0);

这里的best_point就是最优点，向量best_point的第一个元素是最优点的横坐标，第二个元素是最优点的纵坐标，mis_dis是最短距离之和。

5. 输出数据和画图

%% 输出数据
disp(['最短路径之和为',num2str(min_dis),...
      ',横坐标为',num2str(best_point(1)),',纵坐标为',num2str(best_point(2))]);
%% 作图
figure
plot(data(:,1),data(:,2),'bo');
axis equal
hold on
plot(best_point(1),best_point(2),'rp');

结果如下：

最短路径之和为18.7238,横坐标为5.5534,纵坐标为5.4711

6. 完整代码和推广

%% 初始化
data = [1 8;3 7;6 4;7 5;9 2;5 8];
n = length(data);
distance = zeros(1, n);
%% 目标函数
fun = @(x) sum(sqrt(sum((data-x).^2, 2)));
% 如果出现矩阵维度报错，请使用下方的fun
%fun = @(x) sum(sqrt(sum((data-repmat(x,n,1)).^2, 2)));
%% 在已知点中所搜费马点
for i = 1:n
    distance(i) = fun(data(i,:));
end     
[~,index] = min(distance);
%% 将该点作为初始值，使用fminsearch搜索
x0 = data(index,:);
[best_point, min_dis] = fminsearch(fun, x0);
%% 输出数据
disp(['最短路径之和为',num2str(min_dis),...
      ',横坐标为',num2str(best_point(1)),',纵坐标为',num2str(best_point(2))]);
%% 作图
figure
plot(data(:,1),data(:,2),'bo');
axis equal
hold on
plot(best_point(1),best_point(2),'rp');

简单将代码封装后，调用方法有如下两种，代码下载链接

7. 交流讨论