粒子动力学

Particle Dynamics

我正在为 3D 粒子建模 space。

{0} 粒子在时间 t0 从已知位置 P0 开始,速度为 V0。使用已知的 P-1 在 t-1 的先前位置计算速度。

{1} 粒子的目标是在 t1 以已知速度 V1 到达 P1。

{..} 粒子尽可能快地移动,没有受到一组约束的冲击(C1 连续),这些约束独立地限制沿 x、y 和 z 的加速度。沿 x、y 和 z 的最大值 acceleration/deceleration 是已知的,分别是 Xa、Ya 和 Za。加速度沿 x、y 和 z 的最大变化率由 Xr、Yr 和 Zr 定义。

{n} 经过未知数量的时间步后,它在某个时间(比如 tn)以速度 Vn 到达 Pn。

{n+1} 在tn+1移动到Pn+1。

我的问题是计算从P0 到Pn 转换的最短时间并生成其中间位置和速度方向。次要目标是平稳加速,而不是应用导致颠簸的加速。

当前方法:

  1. 找到从起点 P0 到终点 Pn 对齐所需时间最长的维度 {x、y 或 z}。这将是关键维度并将决定总时间。这相当简单,我可以为此写一些东西。

  2. 在所有维度上从 P0 到 Pn 平滑无抖动地插值,使得 Pn 处的速度符合预期。我不确定,如何解决这个问题。

任何已经这样做的 inputs/physics 引擎都会有用。这是一个商业项目,我不能依赖于具有限制性许可的大型 3rd 方库。

注意:P0 和 Pn 处的粒子几乎没有加速度。

如果我没理解错的话,你有一个点(P0, V0),有V0 = P0 - P-1,还有一个点(Pn, Vn),有Vn = Pn - Pn-1,你想找到最少的通过调整每个时间步的加速度来获得中间点。

让我们定义 ti 处的加速度:Ai = Vi - Vi-1,使用 abs(Ai) <= mA。在这里,由于问题是 axis-independant,abs 是 member-wise 绝对值而不是范数(或矢量幅度),mA 是最大加速度矢量,在每个尺寸。我们还考虑 Pn > P0 (member-wise).

据此,我们得到 Vi = Vi-1 + Ai,因此 Pi = Pi-1 + Vi-1 + Ai

如果您需要以最快的速度从某个点到达另一个点,那么显而易见的事情是,无论初始速度如何,都尽可能地加速,直到达到目标。但是,由于您的问题是离散的,并且您有一个终端速度 Vn,因此使用该方法可能会走得太远并且具有不同的终端速度。

然而,你可以反过来做同样的事情,从终点开始。如果你同时从两个点开始,你将在每个维度上创建两条路径相互交叉(不一定在 3D 中交叉,但是,在每个维度中,两条路径的相对方向在某个 "crossing" 点发生变化)。

让我们举一个one-dimensional的例子。 (P0, V0) = (0, -2)(Pn, Vn) = (35, -1),以及 mA = 1
第一条路径 Ai = mA 是这样的:

(0, -2) -> (-1, -1) -> (-1, 0) -> (0, 1) -> (2, 2) ->
(5, 3) -> (9, 4) -> (14, 5) -> (20, 6) -> (27, 7) -> ...

第二条路径 Ai = -mA 但相反,是这样的:

(35, -1) <- (36, 0) <- (36, 1) <- (35, 2) <- (33, 3) <-
(30, 4) <- (26, 5) <- (21, 6) <- (15, 7) <- ...

您可以看到路径以相同的速度在 20 和 21 之间的某处交叉。这为您提供了所需路径中最快的加速和减速部分,但是这两个部分没有连接。但是,通过找到速度相同的最近点很容易将它们连接起来;我们称这些点为 PqPr。这里,Pq = (20, 6)Pr = (21, 6)。由于该速度是在当前点和先前点之间计算的,因此取 Pq 之前的点(示例中的 Pq-1(14, 5))和点 Pr,然后尝试将它们连接起来.

如果Pq >= Pr >= Pq - 2mA,则Pq-1不变,PrVr = Pr - Pq-1直接连接即可。

否则,采用 Pq-2Pr-1(其中 Vr-1 = Vr - mA,因为它是反向的)并尝试通过添加中间点来连接它们。由于这些点的速度差 mA,您只能搜索具有相同速度 Vs 的中间点,这样 Vq-2 <= Vs <= Vr-1.

如果仍然找不到解决方案,则采用 Pq-3Pr-2 并使用更多中间点重复该过程。

在我举的例子中,Pq < Pr,所以我们必须尝试使用​​ Pq-2 = (9, 4)Pr-1 = (26, 5)。我们可以用 3 个点的序列连接它们,例如 (9, 4) -> (13, 4) -> (17, 4) -> (21, 4) -> (26, 5).

无论如何,此方法将为您提供最少数量的中间点,这意味着 P0Pn 之间的最快路径。

如果您随后想减少抖动,那么您可以忘记之前计算的点,并使用您现在知道的最小的 数量 个点进行插值。

在尝试了一些想法之后,我想到了另一个解决方案,如果正确完成,它比我之前的答案更准确,也可能更快。然而,它非常复杂并且需要相当多的数学知识,尽管不是很复杂的数学知识。此外,这是一项正在进行的工作:我仍在调查某些领域。尽管如此,根据我的尝试,它确实已经产生了非常好的结果。

问题

定义和目标

在整个答案中,p[n] 指的是第 n 个点的位置,v[n] 指的是它的速度,a[n] 指的是它的加速度,j[n] 指的是它的加加速度(加速度的导数)。第 n 个点的速度只取决于它的位置和前一点的位置。与加速度和加加速度类似,但分别具有点速度和加速度。

我们有一个起点和一个终点,分别为 p[0]p[n],两者都具有相关的速度 v[0]v[n]。目标是在中间放置 n-1 个点,任意 n,这样,沿着 X、Y 和 Z 轴,加速度和加加速度的绝对值在任何这些点(和在 p[n]) 低于某些限制,分别是 aMaxXaMaxYaMaxZ 加速度,以及 jMaxXjMaxYjMaxZ 混蛋。

我们要查找的是所有 i ∈ [1; n-1]p[i] 的值。因为 p[i] = p[i-1] + v[i],这与查找 v[i] 相同。同理,用v[i] = v[i-1] + a[i]a[i] = a[i-1] + j[i],找a[i]j[i]也是一样的。

a[0]a[n+1] 假定为零。

观察和简化

因为问题的约束与维度无关,所以我们可以分别求解三个维度中的每一个,只要每种情况下得到的点数相同即可。因此,我只打算解决 one-dimensional 版本的问题,使用 aMaxjMax,与轴无关。

*[WIP]* 确定最坏的情况先解决,再解决其他的,知道点数。

两个给定点的实际位置无关紧要,重要的是它们之间的相对距离,我们可以定义为P = p[n] - p[0]。我们还定义范围 R = [1; n]R* = [1; n+1].

由于问题的离散性,我们可以得到以下方程。请注意,∑{i∈R}(x[i])i∈R 的所有 x[i] 的总和。

Ⓐ ∑{i∈R}(v[i]) = P
Ⓑ ∑{i∈R}(a[i]) = v[n] - v[0]
Ⓧ ∑{i∈R*}(j[i]) = 0

Ⓧ来自于a[0] = a[n+1] = 0.
的假设 由Ⓐ和v[i] = v[i-1] + a[i], i∈R,我们可以推导出:

Ⓒ ∑{i∈R}((n+1-i)*a[i]) = P - n*v[0]

同理,由Ⓑ、Ⓒ、a[i] = a[i-1] + j[i], i∈R,可以推导出:

Ⓨ ∑{i∈R}((n+1-i)*j[i]) = v[n] - v[0]
Ⓩ ∑{i∈R}(T[n+1-i]*j[i]) = P - n*v[0]

这里,T[n]是第n个三角数,定义为T[n] = n*(n+1)/2

方程式Ⓧ、Ⓨ、Ⓩ与接下来的部分相关

方法

为了最小化n,我们可以从n(1, 2?)的小值开始求解。然后,如果 max{i∈R}(abs(a[i])) > aMaxmax{i∈R}(abs(j[i])) > jMax,我们可以递增 n 并重复该过程。

*[WIP]* 找到 n 的下限以避免从 n 的小值进行不必要的计算。或估计 n 的正确值并通过测试解决方案查明它。

找到解决方案需要为所有 i∈R* 找到 j[i] 的值。我还没有找到 j[i] 的最佳形式,但是定义 j*[i]r[i]s[i] 这样
j[i] = j*[i] + r[i]v[0] + s[i]v[n]
效果很好。

*[WIP]* 为 j[i]

寻找更好的形式

通过这样做,我们将 n-1 个未知数(j[i], i∈R,请注意 j[n+1] = -∑{i∈R}(j[i]))转换为 3(n-1) 个更容易找到的未知数。以下是我们现在可以从 Ⓧ、Ⓨ 和 Ⓩ 中推断出的一些事情。

∑{i∈R*}(r[i]) = 0
∑{i∈R*}(s[i]) = 0
∑{i∈R}((n+1-i)*r[i]) = -1
∑{i∈R}((n+1-i)*s[i]) = 1
∑{i∈R}(T(n+1-i)*r[i]) = -n
∑{i∈R}(T(n+1-i)*s[i]) = 0

提醒一下,这里是Ⓧ、Ⓨ和Ⓩ。

Ⓧ ∑{i∈R*}(j[i]) + j[n+1] = 0
Ⓨ ∑{i∈R}((n+1-i)*j[i]) = v[n] - v[0]
Ⓩ ∑{i∈R}(T[n+1-i]*j[i]) = P - n*v[0]

现在的目标是找到足够的特例来帮助我们确定这些未知数。

特殊情况

v[0] = v[n] = 0

通过研究加加速度的值,我观察到将所有 j[i], i∈R* 作为抛物线的一部分可以在最小化加加速度和加速度方面产生出色的结果。虽然它不是最合适的,但我还没有找到更好的。

来自抛物线的 jerk 值背后的直觉是,如果 position 的值要遵循多项式,那么它的次数必须至少为 5,并且可以是 5。如果你这样做就更容易理解考虑 4 次多项式之后的速度值。设置 v[0]v[n] 的约束条件 a[0] = a[n+1] = 0,并且其在 [0; n] 上的积分必须等于 P,此多项式的次数必须至少为4. 这适用于连续和离散的情况。最后,似乎取最小的度数会导致更平滑的混蛋,也更容易计算。

这是一个连续情况的例子,紫色是位置,蓝色是速度,黄色是加速度,加加速度n 红色.

如果你想玩这个,这里是如何根据 np[0]p[n]v[0]v[0] 定义位置曲线v[n](其他的只是导数)

a = (-3(v[n]+v[0]) + 6(p[n]-p[0])) / n^5
b = (n(7v[n]+8v[0]) - 15(p[n]-p[0])) / n^4
c = (-n(4v[n]+6v[0]) + 10(p[n]-p[0])) / n^3
p[x] = ax^5 + bx^4 + cx^3 + v[0]x + p[0]

如果v[0] = v[n] = 0,则j[i] = j*[i], i∈R*。这意味着值 j*[i] 遵循二次多项式。所以我们想要找到 αβγ 使得 Ⓟ 成立。

Ⓟ j*[i] = αi^2 + βi + γ, i∈R*

从 Ⓧ、Ⓨ 和 Ⓩ 遵循这些等式。

α*∑{i∈R*}(i^2) + β*∑{i∈R*}(i) + c*∑{i∈R*}(1) = 0
α*∑{i∈R}((n+1-i)*i^2) + β*∑{i∈R}((n+1-i)*i) + c*∑{i∈R}(n+1-i) = 0
α*∑{i∈R}(T(n+1-i)*i^2) + β*∑{i∈R}(T(n+1-i)*i) + c*∑{i∈R}(T(n+1-i)) = P

求解该系统得到αβγ,可与Ⓟ一起计算j*[i], i∈R*。注意j*[i] = j*[n+2-i],所以只需要计算上半部分

v[0] = v[n] = 1/n

如果v[0] = v[n] = 1/n,则j[i] = 0, i∈R*。这意味着 Ⓠ 成立。

Ⓠ r[i] + s[i] = -n*j[i], i∈R*

v[0] = 0, j[i∈L] = J, j[h] = 0, j[i∈U] = -J

LU分别是R*的下半部分和上半部分,h是中间的值,如果n+1是奇数.也就是说:

if n is odd:
L = [1; (n+1)/2]
U = [(n+3)/2; n+1]
if n is even:
L = [1; n/2]
h = n/2+1
U = [n/2+2; n]

这种特殊情况对应于 p[0]p[n] 之间的最大整体加速度,同时最小化 abs(j[i]), i∈R*。在这里,Ⓩ 给出了以下等式。

∑{i∈R}(T[n+1-i]*j[i]) = P
∑{i∈L}(T[n+1-i])*j[1] + ∑{i∈U}(T[n+1-i])*j[n+1] = P
j[1] = P / [ ∑{i∈L}(T[n+1-i]) - ∑{i∈U}(T[n+1-i]) ]

这给出 j[1],因此每个 j[i], i∈R*。然后我们可以使用 Ⓨ.

计算 v[n]

拼凑

对于 v[0]v[n]P 的某些值,每个特例都为我们提供了形式为
的关系 αj*[i] + βr[i] + γs[i] = δ.
通过处理三个特殊情况(假设它们不相似,意味着不给出相同的关系),我们有一个由三个方程组成的系统,一旦求解,就会给出 j*[i]r[i]s[i] 所有 i∈R*.

因此,对于n的每个值,我们可以根据v[0]v[n]P计算j[i]的值。它们可以预先计算,这意味着测试它们的任何 n 值都可以非常快。因此,只要我们预先计算的值达到足够大的值 n.

回答

我建议你使用以下功能:

X(n) = Xstart + Vxstart n+ (-6xstart +3Vx开始+6x结束-3Vx结束+c/2 ) n^2 + (8xstart+3Vxstart-8xend+5Vxend-c) n^3 + (-3Xstart-Vxstart+3xend-2Vxend+c/2) n^4

(对于每个坐标 X、Y、Z)

这里有一些图表,我为每个样本取了 c=3:

对于 xstart=1、vstart=1、xend=3、vstart=-2,这给出:

X(n)= 1 + n + 16 n^2 -25 n^3 + 10 n^4

对于 xstart = -4、vstart =-4、xend = 4、vend = 0,这给出:

(-4 -4n +61n^2 -78n^3 + 29yn^4)

其中c是0.1到5之间的数字,由你决定,c越高,函数运行到那个点的速度越快(但如果c > 4,它可能必须返回). (见下图)。

多项式来自以下计算:其中 a=x0,b=v0,c=xe,d=v2,e=魔法常数

说明

根据 Nelfeal 的回答,我的想法是尝试用多项式解决给定的问题。

我们可以改变问题,定义一个新的 Axis 进入 P[last]-P[0],将问题减少到维度 1。

我们可以用连续数学而不是离散数学来思考问题(比如用函数代替序列),回到离散世界,离散世界只是连续的一个特例。

我们可以把时间和space的单位换一下,让时间为1,距离为1,这样问题就简化为

Find a function which satisfies the following :

  1. (0) = 0 and (1) = 1
  2. '(0) = 0 and '(1) = 0
  3. For x∈ℝ |''(x)| < c, where c is the max speed

我们有

P(X) = ∑{i∈ℕ} Ai Xi

P'(X) = ∑{i∈ℕ} (i+1) Ai+1 Xi

P''(X) = ∑{i∈ℕ} (i+2)(i+1) Ai+2 Xi

我们需要:

  1. P(0) = 0
  2. P(1) = 1
  3. P'(0) = 0
  4. P'(1) = 0
  5. -c <= P''(x) <= c

这意味着:

a0 = 0(从 1.) a1 = 0 (来自 3.)

P(1) = ∑{i∈ℕ} Ai = 1

P'(1) = ∑{i∈ℕ} (i+1) Ai = 0

P''(x) = ∑{i∈ℕ} (i+2)(i+1) Ai Xi 在 [-c,c]

第三个方程是最复杂的一个,可以简化为 P(1) = c。

我们将让 c vary 看看有什么变化。

反转 3x3 矩阵后,我们得到以下结果:

P(x) = (c/2+6) x^2 - (c+8) x^3 + (c/2+3) x^4

对于 c=0.15,这给出:

对于 c=1,这给出:

对于 c=4,我们看到反弹:

如果我们将 c 从 0.1 取到 6,我们将得到以下 3d 图:

请注意,我们已经解决了 4 次多项式的问题,但您可以在更高的次数(如果您愿意,最多可以达到 10 次)上执行相同的操作,以便在您的函数中获得更多可能性。