Python中的最小化以找到两点之间的最短路径

Question

我试图找到两点 (0,0) 和 (1000,-100) 之间的最短路径。路径将由 7 阶多项式函数定义：

p(x) = a0 + a1*x + a2*x^2 + ... + a7*x^7

为此，我尝试最小化从多项式函数计算总路径长度的函数：

length = int 从 0 到 1000 of { sqrt(1 + (dp(x)/dx)^2 ) }

显然，正确的解决方案将是一条直线，但稍后我想为问题添加约束。这应该是第一种方法。

我实现的代码是：

import numpy as np
import matplotlib.pyplot as plt
import math
import sys
import scipy

def path_tracer(a,x):
     return a[0] + a[1]*x + a[2]*x**2 + a[3]*x**3 + a[4]*x**4 + a[5]*x**5 + a[6]*x**6 + a[7]*x**7


def lof(a):
     upper_lim = a[8]

     L = lambda x: np.sqrt(1 + (a[1] + 2*a[2]*x + 3*a[3]*x**2 + 4*a[4]*x**3 + 5*a[5]*x**4 + 6*a[6]*x**5 + 7*a[7]*x**6)**2)
     length_of_path = scipy.integrate.quad(L,0,upper_lim)

     return length_of_path[0]

a = np.array([-4E-11, -.4146,.0003,-7e-8,0,0,0,0,1000]) # [polynomial parameters, x end point]

xx = np.linspace(0,1200,1200)
y = [path_tracer(a,x) for x in xx]

cons = ({'type': 'eq', 'fun': lambda x:path_tracer(a,a[8])+50})
c = scipy.optimize.minimize(lof, a, constraints = cons)
print(c)

当我运行它时，最小化例程失败并返回初始参数不变。输出是：

fun: 1022.9651540965604
     jac: array([  0.00000000e+00,  -1.78130722e+02,  -1.17327499e+05,
        -7.62458172e+07,   9.42803815e+11,   9.99924786e+14,
         9.99999921e+17,   1.00000000e+21,   1.00029755e+00])
 message: 'Singular matrix C in LSQ subproblem'
    nfev: 11
     nit: 1
    njev: 1
  status: 6
 success: False
       x: array([ -4.00000000e-11,  -4.14600000e-01,   3.00000000e-04,
        -7.00000000e-08,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   1.00000000e+03])

是我做错了什么还是例行公事不适合解决这类问题？如果是这样，Python中是否有替代方案？

Answer 1

您可以使用此例程，但您的方法存在一些问题：

• 多项式的域应该归一化为合理的值，例如 [0, 1]。这使得优化更容易。完成优化后，您可以恢复它

• 你可以使用polyval和相关函数来简化代码

• 这个问题的最佳解决方案很明显是-0.1 x，所以我不确定你为什么觉得需要优化。

一个可行的解决方案是

import numpy as np
import scipy.optimize

x = np.linspace(0, 1, 1000)

def obj_fun(p):
    deriv = np.polyval(np.polyder(p), x)
    return np.sum(np.sqrt(1 + deriv ** 2))

cons = ({'type': 'eq', 'fun': lambda p: np.polyval(p, [0, 1]) - [0, -100]})

p0 = np.zeros(8)
c = scipy.optimize.minimize(obj_fun, p0, constraints = cons)

我们可以在哪里绘制结果

import matplotlib.pyplot as plt
plt.plot(np.polyval(c.x, x), label='result')
plt.plot(-100 * x, label='optimal')
plt.legend()