Minimization in Python to find shortest path between two points

Question

I'm trying to find the shortest path between two points, (0,0) and (1000,-100). The path is to be defined by a 7th order polynomial function:

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

To do so I tried to minimize the function that calculates the total path length from the polynomial function:

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

Obviously the correct solution will be a linear line, however later on I want to add constraints to the problem. This one was supposed to be a first approach.

The code I implemented was:

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)

When I ran it however the minimization routine fails and returns the initial parameters unchanged. The output is:

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

Am I doing something wrong or is the routine just not appropriate to solve this kind of problems? If so, is there an alternative in Python?

score 0 · Answer 1 · answered Sep 01 '17 at 13:40

You can use this routine, but there are some problems with your approach:

The domain of the polynomial should be normalized to something reasonable, like [0, 1]. This makes the optimization much easier. You can revert this after you are done with the optimization
You could simplify the code by using polyval and related functions
The optimal solution to this is quite obviously -0.1 x, so I'm not sure why you feel the need to optimize.

A solution that works is

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)

Where we can plot the result

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

Thanks, that worked. Why is it it couldn't find the optimal solution? I'd expect this to be an easy problem to solve for the routine.Regarding the fact that there's an obvious solution to the problem, I'm aware of that. I wanted to start with the easiest case before making it more complex (I'll have 4 paths, the radius of curvature must always be below a certain threshold, derivatives at the initial point and end point must be zero). — Daniel Jaló, Sep 01 '17 at 14:15

Minimization in Python to find shortest path between two points

1 Answers1