Lingo/scipy路径最短问题

*参考链接：http://blog.sina.com.cn/s/blog_a16714bf0101deuz.html*

对于在平面区域上随机分布的$N$个点，求一点$(a,b)$ 使该点到其他$N$ 个点的距离和最短的问题，用$MATLAB$ 和 $LINGO$ （http://pan.baidu.com/s/1nvjcdF7）可以简化求解，具体代码和思想如下：

---

$LINGO$ 程序

model:
sets:
	position/1..18/:x,y;!随机生成的18个点的坐标;
	ab/1/:a,b;!待求的坐标;
endsets

data:
	@text('G:\matlab programing\test\dataxy.txt')=x,y;!保存数据到G:\matlab programing\test\目录下;
	@text('G:\matlab programing\test\dataab.txt')=a,b;
enddata
x(1)=1+4*@rand(0.12345);
y(1)=1+4*@rand(0.25);
@for(position(i)|i #ge# 2 :x(i)=1+4*@rand(x(i-1)));!随机产生1~5中的18个点;
@for(position(i)|i #ge# 2 :y(i)=1+4*@rand(y(i-1)));
[obj]min=@sum(position(i):@sqrt((x(i)-a(1))^2+(y(i)-b(1))^2));!目标函数;
@bnd(1,a(1),5);!意思是待求的点（a,b）中a的范围(bound)在[1,5];
@bnd(1,b(1),5);
end

$MATLAB$ 程序

clear;
clc;
close all;
load('dataxy.txt'); %读取点的坐标
load('dataab.txt'); %读取所求点的坐标
hold on;
plot(dataab(1),dataab(2),'o','markersize',15,'markerfacecolor','r'); %绘制出所求点(a,b)
plot(dataxy(:,1),dataxy(:,2),'or','markersize',15,'markerfacecolor','b'); %绘制出其他各点
set(gcf,'Color','w'); %修改背景色为白色 Get Current Figure获取当前图像属性
set(gca,'fontsize',16); %修改图像字体大小
grid off; %去除网格
dataab=repmat(dataab,18,1); %把(a,b)坐标的值生成和(x,y)一样多的矩阵
P=[dataab(:,1)';dataxy(:,1)'];
Q=[dataab(:,2)';dataxy(:,2)'];
plot(P,Q,'g','LineWidth',2);
xlabel('$x$','interpreter','latex');
ylabel('$y$','interpreter','latex');

最终的结果如下图：

2019-04-14更新

该路径最短问题可以归结为最优化问题，即求解一个点到其他点距离的最小值。

$$minimizef(x,y)=\sum ((x-data{x}_i{)}^2+(y-data{y}_i{)}^2)$$
这里为了方便计算，使用了距离的平方，没有开根号计算。

根据前面的目标函数，我们写成符合scipy.optimzie.minimize的函数输入格式，最后完整的程序如下：

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

datax = 10*np.random.randn(15)
datay = 10*np.random.rand(15)

def f(x):
    return sum((x[0]-datax[i])**2+(x[1]-datay[i])**2 for i in range(len(datax)))

# 初值
x0 = [1,1]
# 计算结果
opt = minimize(f,x0=x0)
print(opt)

plt.figure()
plt.scatter(datax,datay,label='Original points')
plt.scatter(opt.x[0], opt.x[1], s=240, marker='h',label='Center point')
plt.legend()
# plt.show()
for i in range(len(datax)):
    plt.plot([opt.x[0],datax[i]],[opt.x[1],datay[i]])
plt.show()

结果如下：

      fun: 762.1598697487602
 hess_inv: array([[ 3.33359580e-02, -8.50862577e-06],
       [-8.50862577e-06,  3.33328322e-02]])
      jac: array([ 0.00000000e+00, -7.62939453e-06])
  message: 'Optimization terminated successfully.'
     nfev: 28
      nit: 5
     njev: 7
   status: 0
  success: True
        x: array([0.89233268, 5.09306968])

如果想用matlab作为个人学习使用，可以微信扫描下面公众号二维码，关注“向芽塔”，后台回复
“MATLAB”，就可以获取丰富的资源啦（支持WIN,MAC,Linux三大系统）