如何使用 scipy.optimize.minimize 求出参数化曲线的走时？
How to use scipy.optimize.minimize to find travel time of a parameterized curve?
问题
我想确定一个物体沿参数化曲线移动一定距离的旅行时间。我已经学会了how to do this mathematically, but I think there should be a better way of implementing this in Python using scipy.optimize.minimize。但是，由于某种原因，我无法让它工作。结果总是转到 +inf，即使我最初的猜测应该相当接近。我做错了什么？
更具体的问题：
给定按时间参数化的曲线（在 2d 中），您选择任意时间点 (t_end)，这也指定了曲线上的特定点。从那里你沿着曲线回到过去，直到行进的路径等于某个任意距离（d_min）。我想知道的是沿着这条曲线的旅行时间，或者换句话说，给定 t_end 和 d_min，t_start 是多少，这样沿着曲线的线积分从t_start到t_end等于d_min.
下面是一个仅包含关键部分的 MWE：
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)}'
)
背景
我正在尝试估计我应该推断出一条轨迹的时间回溯多远，以查看它是否与特定区域相交。因此，作为初始约束，我使用从最早数据点到椭圆的最短距离，并使用物体的平均速度和加速度来估计该距离的行进时间。这是一个更大的 MWE，它提供了更多上下文：
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()
我只会解决你的关键部分。以下是我想到的几点：
根据评论，您基本上想通过最小化 objective q(t) = F(t)-d 来求解方程 F(t) = d。但是，从数学上讲，这与一般情况下 不同 。要了解原因，让我们考虑二次函数 F(t) = t^2 和 d = 1。然后，t = 1 求解方程 F(t) = d。但是，最小化 objective q(t) = t^2 - 1 会产生局部最小值 t = 0，函数值为 objective q(0) = -1。显然，0*0 ≠ 1。请注意，当且仅当 q(t) = 0 时，F(t) = d，即 objective q 的 objective 值为 0。因此，您需要objective 值为 0 的最小值。并且由于根据定义任何规范都是 non-negative，我们只是最小化函数的 euclidean norm，即我们最小化 p(t) = ||q （吨）|| = ||F(t) - d||.
因为你有一个（可能non-convex）一个变量的标量优化问题，强烈建议使用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))
这会产生 t_start = 1353.71。请注意，我忽略了您的初始界限，因为您的方程在该区间内没有解。
问题
我想确定一个物体沿参数化曲线移动一定距离的旅行时间。我已经学会了how to do this mathematically, but I think there should be a better way of implementing this in Python using scipy.optimize.minimize。但是，由于某种原因，我无法让它工作。结果总是转到 +inf，即使我最初的猜测应该相当接近。我做错了什么？
更具体的问题：
给定按时间参数化的曲线（在 2d 中），您选择任意时间点 (t_end)，这也指定了曲线上的特定点。从那里你沿着曲线回到过去，直到行进的路径等于某个任意距离（d_min）。我想知道的是沿着这条曲线的旅行时间，或者换句话说，给定 t_end 和 d_min，t_start 是多少，这样沿着曲线的线积分从t_start到t_end等于d_min.
下面是一个仅包含关键部分的 MWE：
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)}'
)
背景
我正在尝试估计我应该推断出一条轨迹的时间回溯多远，以查看它是否与特定区域相交。因此，作为初始约束，我使用从最早数据点到椭圆的最短距离，并使用物体的平均速度和加速度来估计该距离的行进时间。这是一个更大的 MWE，它提供了更多上下文：
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()
我只会解决你的关键部分。以下是我想到的几点：
根据评论，您基本上想通过最小化 objective q(t) = F(t)-d 来求解方程 F(t) = d。但是，从数学上讲，这与一般情况下 不同 。要了解原因，让我们考虑二次函数 F(t) = t^2 和 d = 1。然后，t = 1 求解方程 F(t) = d。但是，最小化 objective q(t) = t^2 - 1 会产生局部最小值 t = 0，函数值为 objective q(0) = -1。显然，0*0 ≠ 1。请注意，当且仅当 q(t) = 0 时，F(t) = d，即 objective q 的 objective 值为 0。因此，您需要objective 值为 0 的最小值。并且由于根据定义任何规范都是 non-negative，我们只是最小化函数的 euclidean norm，即我们最小化 p(t) = ||q （吨）|| = ||F(t) - d||.
因为你有一个（可能non-convex）一个变量的标量优化问题，强烈建议使用
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))
这会产生 t_start = 1353.71。请注意，我忽略了您的初始界限，因为您的方程在该区间内没有解。