Minimisation en Python pour trouver le plus court chemin entre deux points

Question

J'essaie de trouver le chemin le plus court entre deux points, (0,0) et (1000, -100). Le chemin doit être défini par une fonction polynomiale septième ordre:

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

Pour ce faire, J'essayé de minimiser la fonction qui permet de calculer la longueur totale du trajet à partir de la fonction polynomiale:

longueur = int de 0 à 1000 de {sqrt (1 + (dp (x)/dx)^2)}

Évidemment, la solution correcte sera une ligne linéaire, mais plus tard je veux ajouter des contraintes au problème. Celui-ci était censé être une première approche.

Le code I a été mis en œuvre:

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) 

Quand je l'ai couru mais la routine de minimisation échoue et renvoie les paramètres initiaux inchangés. La sortie est:

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]) 

Est-ce que je fais quelque chose de mal ou est-ce que la routine n'est pas appropriée pour résoudre ce genre de problème? Si oui, existe-t-il une alternative en Python?

Answer 1

Vous pouvez utiliser cette routine, mais il y a des problèmes avec votre approche:

• Le domaine du polynôme doit être normalisé à quelque chose de raisonnable, comme [0, 1]. Cela rend l'optimisation beaucoup plus facile. Vous pouvez revenir cette fois que vous avez terminé avec l'optimisation

• Vous pouvez simplifier le code en utilisant polyval et fonctions connexes

• La solution à cette question est bien évidemment -0.1 x, donc je ne sais pas pourquoi vous ressentir le besoin d'optimiser.

Une solution qui fonctionne est

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) 

Là où nous pouvons tracer le résultat

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