1

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

enter image description here

mapf
  • 1,906
  • 1
  • 14
  • 40
  • 1
    Can you come up with a minimalistic example? If this question is answered like this it is unlikely to benefit other people with a similar problem – Marat Mar 11 '22 at 15:01
  • Hey, that's fair of course, I just don't know how I can make it more abstract. I can tell you that this hasn't worked for not just that one track, but all of the hundreds of tracks I tested so far, so I'm pretty sure it must be a systematic error. On the other hand though, determining the travel time for a specific distance along an arbitrary parameterized curve is a pretty broad application that should be interesting to many people. – mapf Mar 11 '22 at 15:38
  • @Marat so I just separted the crucial part of the MWE for a better overview. But I don't how what else to reduce. – mapf Mar 11 '22 at 15:51
  • Sorry, I still don't get what is fixed and what is an independent variable. From the code it looks like we choose the best `t_start`, which doesn't make a lot of sense to me – Marat Mar 11 '22 at 16:03
  • Sorry for the confusion! Yes, I want to find "the best" `t_start`. The idea is this: Imagine you have a curve (in 2d) parameterized by time. Now 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_start` to `t_end` is equal to `d_min`. – mapf Mar 11 '22 at 18:31
  • It's unfortunately a bit awkward to explain without being able to use properly formatted formulas. It's maybe a little bit easier to understand if you look at the [other post on the mathematics](https://math.stackexchange.com/a/4398673/760701). – mapf Mar 11 '22 at 18:33

1 Answers1

3

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.

joni
  • 6,840
  • 2
  • 13
  • 20
  • Hi, thank you! "According to the comments, you basically want to solve an equation F(t) = d by minimizing the objective F(t)-d" yes, that's exactly it. Your approach sounds very promising, however, I'm surprised by the result, as `t_start` should be smaller than `t_end`. Of course there are always two solutions, one in each direction of the curve. But I don't understand why it would yield the positive solution. – mapf Mar 11 '22 at 22:13
  • Ok, nvm my brain fart. `t_start = 1353.71` is of course smaller than `t_end = 4211.317`. I made a mistake and gave the wrong `t_end`, it should of course be `t[0]` (which is 0) instead of `t[-1]`. But in my mind I was still thinking of `t_end` as something very close to 0, and thus expected a negative value for `t_start`. Sorry, my bad! With that, it all makes sense now, and the result seems correct: `t_start = -2204.49`. Thanks a lot. – mapf Mar 11 '22 at 23:22
  • Hi, sorry, I was thinking a bit about your first point "However, mathematically, this is not the same in general. Instead, you've got to minimize the euclidian norm of the objective, i.e. ||F(t)-d|| in case you really want to solve an equation by a minimization problem." Could you maybe elaborate a bit why that is? Also, it sounds like minimalization is not the best way to go about it. What would be a better way? – mapf Mar 12 '22 at 15:17
  • And solving equations (of multiple variables) as minimization problems is the standard approach. In the scalar case, it's a bit overkill, though. – joni Mar 13 '22 at 16:32