How to use scipy.optimize.minimize to find travel time of a parameterized curve?

Question

I would like to determine the travel time of an object moving a certain distance along a parameterized curve. I already learned how to do this mathematically, but I think there should be a better way of implementing this in Python using scipy.optimize.minimize. However, for some reason, I cannot get it to work. The result always goes to +inf, even though my initial guess should be fairly close. What am I doing wrong?

The problem in more specific terms:

Given the curve (in 2d) parameterized by time, you choose an arbitrary point in time (t_end), which also specifies a specific point on the curve. From there you go back in time along the curve until the path traveled is equal to some arbitrary distance (d_min). What I want to know is the travel time along this segment of the curve, or in other words, given t_end and d_min, what is t_start so that the line integral along the curve from t_startto t_end is equal to d_min.

Below is an MWE that contains just the crucial part:

import numpy as np
from scipy.optimize import minimize
from scipy.integrate import quad


def speed(t, coeffs):
    vx = np.polyval(np.polyder(coeffs[:, 0]), t)
    vy = np.polyval(np.polyder(coeffs[:, 1]), t)
    return np.hypot(vx, vy)


def line_integral(t_start, t_end, coeffs):
    return quad(speed, t_start, t_end, args=(coeffs,))[0]


def travel_distance(t_start, t_end, coeffs, dist):
    return line_integral(t_start, t_end, coeffs) - dist


coeffs = np.array(
    [
        [-3.05831338e-10,  7.01828123e-10],
        [-2.47221941e-05,  2.38793003e-05],
        [-5.96508764e-02, -2.64191402e-01],
        [ 7.93182941e+02,  6.87323487e+02],
    ]
)
t_start = -2344.434017935311
t_end = 4211.317
d_min = 694.459957887133

res = minimize(
    travel_distance, t_start, args=(t_end, coeffs, d_min),
    method='Nelder-Mead', tol=1e-6, bounds=[(5*t_start, 0), ],
)
print(f'travel time guess: {t_start}, result: {res.x[0]}')
print(
    f'minimal distance: {d_min}\n'
    f'travel distance using guess: {line_integral(t_start, t_end, coeffs)}\n'
    f'travel distance using result: {line_integral(res.x[0], t_end, coeffs)}'
)

Background

I am trying to estimate how far back in time I should extrapolate a track to see if it intersects a specific area. So as an initial constraint I'm using the shortest distance from the earliest data point to the ellipse and estimate the travel time for that distance using the mean speed and acceleration of the object. Here is a larger MWE that provides some more context:

import numpy as np
from scipy.optimize import minimize
from scipy.integrate import quad
from shapely import affinity
from shapely.geometry import LineString, Point


# defining some track
x = [
    793.2978629,  792.79734879, 785.9774051,  785.36722833, 783.40637669,
    782.64515164, 780.3664593,  779.79485231, 777.48085524, 777.28616809,
    774.51204993, 774.11564649, 771.69632972, 771.02165865, 768.48989798,
    767.99362458, 765.28968539, 764.53914535, 762.12842011, 761.41434577,
    758.79663562, 758.18888603, 748.03518007, 747.31157785, 685.70431359,
    684.04913833, 604.28080405, 602.44069938, 503.82330074, 501.54973592,
    383.30603688, 380.51423267, 243.34684644, 239.99098054,  84.7050222,
    80.69489434,
]
y = [
    687.20838319, 683.22855991, 656.29479042, 654.72684098, 645.88731676,
    644.44054023, 635.79673583, 634.20437526, 625.60077485, 623.90628204,
    615.55599689, 614.04715941, 605.48088519, 604.07937581, 595.71970051,
    594.0450717,  586.02784429, 584.17400625, 576.13634702, 574.55298797,
    566.3923443,  564.97163436, 537.48129155, 534.29439651, 406.25295174,
    402.89392009, 293.21722135, 290.93423101, 200.88320562, 198.76239079,
    128.84454584, 127.20496349,  78.34941384,  77.91473383,  50.81598911,
    50.69878393,
]
t = [   
    0.,      12.171,  119.019,  125.316,  159.007,  165.307,  199.008,  
    205.317, 239.024, 245.333,  279.013,  285.314,  319.031,  325.333,
    359.024, 365.32,  399.069,  405.369,  439.051,  445.355,  479.064,
    485.364, 599.99,  612.165, 1199.136, 1212.173, 1800.038, 1811.285,
    2400.027, 2412.165, 3000.033, 3012.197, 3600.007, 3612.168, 4199.096,
    4211.317,
]
w = [
    0.9287239,  0.9287239,  1.,         1.,         0.887263,   0.887263,
    0.9605867,  0.9605867,  0.95934916, 0.95934916, 0.96882457, 0.96882457,
    0.9372613,  0.9372613,  0.963637,   0.963637,   0.91846114, 0.91846114,
    0.9222912,  0.9222912,  0.9175395,  0.9175395,  0.94061875, 0.94061875,
    0.8428271,  0.8428271,  0.894225,   0.894225,   0.8533954,  0.8533954,
    0.84179366, 0.84179366, 0.7897369,  0.7897369,  0.8689509,  0.8689509, 
]

# fitting the track
order = 3
coeffs = np.polyfit(t, np.array([x, y]).T, order, w=w)
mean_speed = np.hypot(*coeffs[-2])
mean_acceleration = np.hypot(*coeffs[-3])

# defining area of interest, following https://gis.stackexchange.com/a/243462
e_props = ((850, 1450), (70, 140))
ellipse = affinity.scale(Point(e_props[0]).buffer(1), *e_props[1])
ellipse = affinity.rotate(ellipse, 90)

# find the shortest distance from the track to the ellipse
track_origin = Point((x[0], y[0]))
d_min = track_origin.distance(ellipse)

# using pq formula to make a first estimate of when the curve will have
# traveled d_min
p = mean_speed/mean_acceleration
q = 2*d_min/mean_acceleration
dt_min = 100
t_start = - (dt_min - p + np.sqrt(p**2 + q))
t_end = t[-1]


""" start of crucial part """


def speed(t, coeffs):
    vx = np.polyval(np.polyder(coeffs[:, 0]), t)
    vy = np.polyval(np.polyder(coeffs[:, 1]), t)
    return np.hypot(vx, vy)


def line_integral(t_start, t_end, coeffs):
    return quad(speed, t_start, t_end, args=(coeffs,))[0]


def travel_distance(t_start, t_end, coeffs, dist):
    return line_integral(t_start, t_end, coeffs) - dist


res = minimize(
    travel_distance, t_start, args=(t_end, coeffs, d_min),
    method='Nelder-Mead', tol=1e-6, bounds=[(5*t_start, 0), ],
)
print(f'travel time guess: {t_start}, result: {res.x[0]}')
print(
    f'minimal distance: {d_min}\n'
    f'travel distance using guess: {line_integral(t_start, t_end, coeffs)}\n'
    f'travel distance using result: {line_integral(res.x[0], t_end, coeffs)}'
)


""" end of crucial part """


# use the the start time the estimate how far back the track should be
# extrapolated to see if it intersects with the ellipse


def position(t, coeffs):
    return np.array([np.polyval(coeffs[:, 0], t), np.polyval(coeffs[:, 1], t)])


n_fit = 1000
t_fit = np.linspace(t_start, t_end, n_fit)
x_fit, y_fit = position(t_fit, coeffs)
curve = LineString(np.column_stack((x_fit, y_fit)))
intersection = curve.intersection(ellipse)

# do some plotting down here
if False:
    from descartes import PolygonPatch
    import matplotlib.pyplot as plt
    from matplotlib.patches import Ellipse
    from matplotlib.collections import LineCollection

    fig, ax = plt.subplots()
    ellipse_artist = PolygonPatch(ellipse, fc='none', ec='gray', lw=2)
    ax.add_patch(ellipse_artist)
    ellipse_center_artist = Ellipse(
        e_props[0], width=8, height=8, color='black', zorder=5
    )
    ax.add_artist(ellipse_center_artist)
    points = np.array([x_fit, y_fit]).T.reshape(-1, 1, 2)
    segments = np.concatenate([points[:-1], points[1:]], axis=1)
    norm = plt.Normalize(t_fit[0], t_fit[-1])
    lc = LineCollection(
        segments, cmap='autumn', norm=norm, alpha=1,
        linewidths=2, capstyle='round', joinstyle='round'
    )
    lc.set_array(t_fit)
    ax.add_collection(lc)
    ax.plot(*intersection.coords[-1], 'x', c='black')
    ax.set_xlim(0, 2048)
    ax.set_ylim(2048, 0)
    fig.show()

Solution

I will only address your crucial part. Here are a few points that crossed my mind:

According to the comments, you basically want to solve an equation F(t) = d by minimizing the objective q(t) = F(t)-d. However, mathematically, this is not the same in general. To see why, let's consider the quadratic function F(t) = t^2 and d = 1. Then, t = 1 solves the equation F(t) = d. However, minimizing the objective q(t) = t^2 - 1 yields the local minimum t = 0 with objective function value q(0) = -1. And obviously, 0*0 ≠ 1. Note that F(t) = d if and only if q(t) = 0, i.e. the objective q has the objective value 0. Therefore, you need a minimum with an objective value of 0. And since any norm is non-negative by definition, we just minimize the euclidean norm of your function, i.e. we minimize p(t) = ||q(t)|| = ||F(t) - d||.
Because you have a (probably non-convex) scalar optimization problem of one variable, it's highly recommended to use the specialized algorithms behind scipy.optimize.minimize_scalar:

from scipy.optimize import minimize_scalar

res = minimize_scalar(lambda t_start: np.sqrt((line_integral(t_start, t_end, coeffs) - dist) ** 2))

This yields t_start = 1353.71. Note that I ignored your initial bounds since there's no solution to your equation within that interval.