2017-09-01 116 views
2

我試圖找到兩點之間(0,0)和(1000,-100)之間的最短路徑。路徑是由一個7階多項式函數來定義:在Python中最小化以找到兩點之間的最短路徑

P(X)= A0 + A1 * X + A2 *的x^2 + ... + a7中* X^7

爲此我試圖最小化,其計算從所述多項式函數的總路徑長度的函數:

長度= INT從0到1000 {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中是否有其他選擇?

回答

0

你可以使用這個程序,但也有一些問題,你的方法:

  • 多項式的域進行歸一化的東西合理,喜歡[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() 

enter image description here

+0

謝謝,這工作。爲什麼它找不到最佳解決方案? 我認爲這是一個很容易解決的問題。關於這個問題有一個明顯的解決方案,我意識到這一點。我希望從最簡單的情況開始,讓它變得更加複雜(我將有4條路徑,曲率半徑必須始終低於某個閾值,初始點和終點的導數必須爲零)。 –

相關問題