0 回顾
回顾 M i n i m u m Minimum Minimum S n a p Snap Snap当中的闭式求解的问题,闭式求解可以把轨迹变的更加顺滑,但是由于没有过多考虑障碍物的信息,可能会与障碍物发生碰撞。所以现在软约束的轨迹优化原理就是在原本顺滑性的基础上考虑代价碰撞函数,以及速度和加速度的限制。
1 软约束的轨迹优化
1.1 方法描述
优化方程:
min
λ
1
f
s
+
λ
2
f
o
+
λ
3
(
f
v
+
f
a
)
\min \quad \lambda_{1} f_{s}+\lambda_{2} f_{o}+\lambda_{3}\left(f_{v}+f_{a}\right)
minλ1fs+λ2fo+λ3(fv+fa),前面的这个顺滑性
f
s
f_{s}
fs就是之前在
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap当中介绍过的表达式,又加入了
f
o
f_{o}
fo碰撞目标函数,以及
f
v
f_{v}
fv和
f
a
f_{a}
fa速度和加速度的目标函数,
λ
1
、
λ
2
、
λ
3
\lambda_{1}、\lambda_{2}、\lambda_{3}
λ1、λ2、λ3是他们的权重系数。求解目标函数的最小值,可以采用梯度下降法。
轨迹依然是分成与
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap当中的多段轨迹。
p
μ
(
t
)
=
{
∑
j
=
0
N
η
1
j
(
t
−
T
0
)
j
T
0
≤
t
≤
T
1
∑
j
=
0
N
η
2
j
(
t
−
T
1
)
j
T
1
≤
t
≤
T
2
⋮
⋮
∑
j
=
0
N
η
M
j
(
t
−
T
M
−
1
)
j
T
M
−
1
≤
t
≤
T
M
p_{\mu}(t)=\left\{\begin{array}{lc} \sum_{j=0}^{N} \eta_{1 j}\left(t-T_{0}\right)^{j} & T_{0} \leq t \leq T_{1} \\ \sum_{j=0}^{N} \eta_{2 j}\left(t-T_{1}\right)^{j} & T_{1} \leq t \leq T_{2} \\ \vdots & \vdots \\ \sum_{j=0}^{N} \eta_{M j}\left(t-T_{M-1}\right)^{j} & T_{M-1} \leq t \leq T_{M} \end{array}\right.
pμ(t)=⎩
⎨
⎧∑j=0Nη1j(t−T0)j∑j=0Nη2j(t−T1)j⋮∑j=0NηMj(t−TM−1)jT0≤t≤T1T1≤t≤T2⋮TM−1≤t≤TM
1.2 顺滑性目标函数的构建及其梯度
软约束轨迹优化方法是用来“完善”
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap方法的,也用到了
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap当中的一些理论,有关
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap的推导可以去看之前的博客。首先就是把参数的映射成每段轨迹两个端点的
P
V
A
PVA
PVA(位置,速度,加速度),
η
=
M
−
1
C
[
d
F
d
P
]
\eta=\mathbf{M}^{-1}\mathbf{C}\left[\begin{array}{l}\mathbf{d}_{F} \\\mathbf{d}_{P}\end{array}\right]
η=M−1C[dFdP]
这样做的目的就是,让系数突出其物理含义。其中,
d
F
d_{F}
dF表示端点需要固定的一些值,
d
P
d_{P}
dP表示需要优化的一些值。
M
−
1
M^{-1}
M−1表示在
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap当中的
A
−
1
M
A^{-1}M
A−1M。
C
C
C是类似的。将上面的
η
\eta
η带入
M
i
n
i
m
u
m
Minimum
Minimum
S
n
a
p
Snap
Snap的顺滑性目标函数,即可得到以下表达式。
顺滑性的目标函数:
f
s
=
[
d
F
d
P
]
T
C
T
M
−
T
Q
M
−
1
C
[
d
F
d
P
]
f_{s}=\left[\begin{array}{l} \mathbf{d}_{F} \\ \mathbf{d}_{P} \end{array}\right]^{T} \mathbf{C}^{T} \mathbf{M}^{-T} \mathbf{Q} \mathbf{M}^{-1} \mathbf{C}\left[\begin{array}{l} \mathbf{d}_{F} \\ \mathbf{d}_{P} \end{array}\right]
fs=[dFdP]TCTM−TQM−1C[dFdP]
把
C
T
M
−
T
Q
M
−
1
C
\mathbf{C}^{T} \mathbf{M}^{-T} \mathbf{Q} \mathbf{M}^{-1} \mathbf{C}
CTM−TQM−1C记作
R
R
R矩阵,如下。
f
s
=
[
d
F
d
P
]
T
[
R
F
F
R
F
P
R
P
F
R
P
P
]
[
d
F
d
P
]
f_{s}=\left[\begin{array}{l} \mathbf{d}_{F} \\ \mathbf{d}_{P} \end{array}\right]^{T}\left[\begin{array}{ll} \mathbf{R}_{F F} & \mathbf{R}_{F P} \\ \mathbf{R}_{P F} & \mathbf{R}_{P P} \end{array}\right]\left[\begin{array}{l} \mathbf{d}_{F} \\ \mathbf{d}_{P} \end{array}\right]
fs=[dFdP]T[RFFRPFRFPRPP][dFdP]
R
R
R阵通过
d
F
、
d
P
d_F、d_P
dF、dP的下标被分为四块,分界的地方分为四块。之后对上述表达式进行求导处理,对
d
p
η
d_{{p}_{\eta}}
dpη进行求导,
η
\eta
η表示的是维度,因为在
x
、
y
、
z
x、y、z
x、y、z三个方向都是用分段的多项式进行表示的,所以三个方向就得考虑
d
p
x
、
d
p
y
、
d
p
z
d_{{p}_{x}}、d_{{p}_{y}}、d_{{p}_{z}}
dpx、dpy、dpz,就可以得出一个通式,如下方最后一行的式子,主要是为了简化符号,表示对一个轴的计算:
J
s
=
[
∂
f
s
∂
d
P
x
,
∂
f
s
∂
d
P
y
,
∂
f
s
∂
d
P
z
]
,
H
s
=
[
∂
2
f
s
∂
d
P
x
2
,
∂
2
f
s
∂
d
P
y
2
,
∂
2
f
s
∂
d
P
z
2
]
,
∂
f
s
∂
d
P
μ
=
2
d
F
T
R
F
P
+
2
d
P
T
R
P
P
,
∂
2
f
s
∂
d
P
μ
2
=
2
R
P
P
T
,
\begin{aligned} &\mathbf{J}_{s}=\left[\frac{\partial f_{s}}{\partial \mathbf{d}_{P_{x}}}, \frac{\partial f_{s}}{\partial \mathbf{d}_{P_{y}}}, \frac{\partial f_{s}}{\partial \mathbf{d}_{P_{z}}}\right], \mathbf{H}_{s}=\left[\frac{\partial^{2} f_{s}}{\partial \mathbf{d}_{P_{x}}^{2}}, \frac{\partial^{2} f_{s}}{\partial \mathbf{d}_{P_{y}}^{2}}, \frac{\partial^{2} f_{s}}{\partial \mathbf{d}_{P_{z}}^{2}}\right], \\ &\frac{\partial f_{s}}{\partial \mathbf{d}_{P_{\mu}}}=2 \mathbf{d}_{F}^{T} \mathbf{R}_{F P}+2 \mathbf{d}_{P}^{T} \mathbf{R}_{P P}, \quad \frac{\partial^{2} f_{s}}{\partial \mathbf{d}_{P \mu}^{2}}=2 \mathbf{R}_{P P}^{T}, \end{aligned}
Js=[∂dPx∂fs,∂dPy∂fs,∂dPz∂fs],Hs=[∂dPx2∂2fs,∂dPy2∂2fs,∂dPz2∂2fs],∂dPμ∂fs=2dFTRFP+2dPTRPP,∂dPμ2∂2fs=2RPPT,
这就是一个对矩阵的求导,想象成数的求导就行。
得到导数后,因为是凸函数,可以得到最优解,因为此时又考虑了碰撞的代价函数以及速度和加速度的代价函数,所以现在考虑顺滑性的代价函数所得到的最优解并不一定是全局的最优解。所以接下来考虑碰撞的目标函数。
1.3 碰撞目标函数的构建及其梯度
因为接下来的速度与加速度目标函数的构建与碰撞的目标函数构建方法大同小异,所以只列举碰撞目标函数的构建方法。
- 引入 P e n a l t y Penalty Penalty F u n c t i o n Function Function,当 U A V UAV UAV距离障碍物很近的时候,将这个函数值设置的非常大,使得 U A V UAV UAV远离这个障碍物;但是当 U A V UAV UAV出于一个非常安全的情况下,就把这个函数值设置的小一点,以便快速通过空旷区域。
- 惩罚函数是可微的,可以选择如下式子所示的指数函数: c ( d ) = α ⋅ exp ( − ( d − d 0 ) / r ) c(d)=\alpha \cdot \exp \left(-\left(d-d_{0}\right) / r\right) c(d)=α⋅exp(−(d−d0)/r) d 0 d_0 d0表示的是一个安全距离, d d d表示的是无人机当前在一个位置(与障碍物之间的距离),假如这个 d d d非常接近障碍物了,所以此时的 d d d肯定是要小于安全距离 d 0 d_0 d0的,所以 d − d 0 d-d_0 d−d0就是一个负数,则 − ( d − d 0 ) -(d-d_0) −(d−d0)就是一个正数,而 r r r是个系数可以自定的,假设 r = 0.1 r=0.1 r=0.1,此时的这个指数 e x p exp exp的数值就会很大,所以整个惩罚函数的值就会非常大,这样就使得无人机会远离障碍物,其实基于软约束的这种方法和人工势场法的思想会不会有点类似?
- 轨迹是一条曲线,轨迹上的点都必须去计算这个惩罚函数。每个点都考虑的话就差不多是个线积分,就沿着轨迹
d
s
d_s
ds进行积分,积分其实就是个累加的想法,把每个点加起来,因为是
d
s
d_s
ds积分,所以可以用每一小段的速度乘以时间
d
t
d_t
dt进行代替
d
s
d_s
ds,这个速度
v
(
t
)
v(t)
v(t)要注意一下,这是
x
、
y
、
z
x、y、z
x、y、z三个方向的速度合并后的一个速度,也就是无人机真实的速度。
f
o
=
∫
T
0
T
M
c
(
p
(
t
)
)
d
s
=
∫
T
0
T
M
c
(
p
(
t
)
)
∥
v
(
t
)
∥
d
t
=
∑
k
=
0
τ
/
δ
t
c
(
p
(
T
k
)
)
∥
v
(
t
)
∥
δ
t
\begin{aligned} f_{o} &=\int_{T_{0}}^{T_{M}} c(p(t)) d s \\ &=\int_{T_{0}}^{T_{M}} c(p(t))\|v(t)\| d t=\sum_{k=0}^{\tau / \delta t} c\left(p\left(\mathcal{T}_{k}\right)\right)\|v(t)\| \delta t \end{aligned}
fo=∫T0TMc(p(t))ds=∫T0TMc(p(t))∥v(t)∥dt=k=0∑τ/δtc(p(Tk))∥v(t)∥δt
其中 ∥ v ( t ) ∥ = ( v x ) 2 + ( v y ) 2 + ( v z ) 2 \|v(t)\| = \sqrt{(v_x)^2 + (v_y)^2 + (v_z)^2} ∥v(t)∥=(vx)2+(vy)2+(vz)2。
上述式子因为没有一个理论的表达式,所以进行离散化处理。 - 求导,过程类似顺滑性约束的目标函数求导,对三个方向分别求导,通式为对
d
p
μ
d_{{p}_{\mu}}
dpμ求导的式子,因为三个轴的求法都是一样的。
J o = [ ∂ f o ∂ d P x , ∂ f o ∂ d P y , ∂ f o ∂ d P z ] ∂ f o ∂ d P μ = ∑ k = 0 τ / δ t { ∇ μ c ( p ( T k ) ) ∥ v ∥ F + c ( p ( T k ) ) v μ ∥ v ∥ G } δ t , ( 9 ) \begin{aligned} &\mathbf{J}_{o}=\left[\frac{\partial f_{o}}{\partial \mathbf{d}_{P_{x}}}, \frac{\partial f_{o}}{\partial \mathbf{d}_{P_{y}}}, \frac{\partial f_{o}}{\partial \mathbf{d}_{P_{z}}}\right] \\ &\frac{\partial f_{o}}{\partial \mathbf{d}_{P_{\mu}}}=\sum_{k=0}^{\tau / \delta t}\left\{\nabla_{\mu} c\left(p\left(\mathcal{T}_{k}\right)\right)\|v\| \mathbf{F}+c\left(p\left(\mathcal{T}_{k}\right)\right) \frac{v_{\mu}}{\|v\|} \mathbf{G}\right\} \delta t,(9) \end{aligned} Jo=[∂dPx∂fo,∂dPy∂fo,∂dPz∂fo]∂dPμ∂fo=k=0∑τ/δt{∇μc(p(Tk))∥v∥F+c(p(Tk))∥v∥vμG}δt,(9)
对 c ( p ( T k ) ) ∥ v ( t ) ∥ δ t c\left(p\left(\mathcal{T}_{k}\right)\right)\|v(t)\| \delta t c(p(Tk))∥v(t)∥δt这部分进行求导,可以使用乘法求导的法则,即 c ( p ( T k ) ) c\left(p\left(\mathcal{T}_{k}\right)\right) c(p(Tk))的导数乘以 ∥ v ( t ) ∥ \|v(t)\| ∥v(t)∥再加上 c ( p ( T k ) ) c\left(p\left(\mathcal{T}_{k}\right)\right) c(p(Tk))乘以 ∥ v ( t ) ∥ \|v(t)\| ∥v(t)∥的导数。上方求导后的 F \mathbf{F} F就是对第一项求导后所产生的项, G \mathbf{G} G就是对第二项求导后产生的项。 ( p ( T k ) ) \left(p\left(\mathcal{T}_{k}\right)\right) (p(Tk))就是 η ∗ T \eta*\mathbf{T} η∗T得来的。 η = M − 1 C [ d F d P ] \eta=\mathbf{M}^{-1}\mathbf{C}\left[\begin{array}{l}\mathbf{d}_{F} \\\mathbf{d}_{P}\end{array}\right] η=M−1C[dFdP],现在引进一个新的符号 L d p L_{dp} Ldp来指代 η \eta η当中对 d p d_p dp求导后的项,所以可以得到 F = T L d p \mathbf{F} = \mathbf{T}\mathbf{L}_{dp} F=TLdp, G = T V m L d p \mathbf{G} = \mathbf{T}\mathbf{V}_{m}\mathbf{L}_{dp} G=TVmLdp。然后对速度求导的时候,因为有根号,所以套用有根号的函数求导公式。 ∇ μ \nabla_{\mu} ∇μ就是碰撞代价在某个轴( x 、 y 、 z x、y、z x、y、z)上的梯度。 - 软约束优化中还得到了 H e s s i a n Hessian Hessian矩阵(一阶导数再求导),了解即可如下: H o = [ ∂ 2 f o ∂ d P x 2 , ∂ 2 f o ∂ d P y 2 , ∂ 2 f o ∂ d P z 2 ] ∂ 2 f o ∂ d P μ 2 = ∑ k = 0 τ / δ t { F T ∇ μ c ( p ( T k ) v μ ∥ v ∥ G + F T ∇ μ 2 c ( p ( T k ) ) ∥ v ∥ F + G T ∇ μ c ( p ( T k ) ) v μ ∥ v ∥ F + G T c ( p ( T k ) ) v μ 2 ∥ v ∥ 3 G } δ t \begin{aligned} &\mathbf{H}_o=\left[\frac{\partial^2 f_o}{\partial \mathbf{d}_{P_x}^2}, \frac{\partial^2 f_o}{\partial \mathbf{d}_{P_y}^2}, \frac{\partial^2 f_o}{\partial \mathbf{d}_{P_z}^2}\right] \\ &\frac{\partial^2 f_o}{\partial \mathbf{d}_{P \mu}^2}=\sum_{k=0}^{\tau / \delta t}\left\{\mathbf { F } ^ { T } \nabla _ { \mu } c \left(p\left(\mathcal{T}_k\right) \frac{v_\mu}{\|v\|} \mathbf{G}+\mathbf{F}^T \nabla_\mu^2 c\left(p\left(\mathcal{T}_k\right)\right)\|v\| \mathbf{F}\right.\right. \\ &\left.+\mathbf{G}^T \nabla_\mu c\left(p\left(\mathcal{T}_k\right)\right) \frac{v_\mu}{\|v\|} \mathbf{F}+\mathbf{G}^T c\left(p\left(\mathcal{T}_k\right)\right) \frac{v_\mu^2}{\|v\|^3} \mathbf{G}\right\} \delta t \end{aligned} Ho=[∂dPx2∂2fo,∂dPy2∂2fo,∂dPz2∂2fo]∂dPμ2∂2fo=k=0∑τ/δt{FT∇μc(p(Tk)∥v∥vμG+FT∇μ2c(p(Tk))∥v∥F+GT∇μc(p(Tk))∥v∥vμF+GTc(p(Tk))∥v∥3vμ2G}δt
- 其他的速度约束与加速度约束差不多构建,就不细说了。
- 这样就有了三部分的雅克比矩阵了, J = λ 1 J s + λ 2 J o + λ 3 ( J v + J a ) \mathbf{J}=\lambda_1 \mathbf{J}_s+\lambda_2 \mathbf{J}_o+\lambda_3\left(\mathbf{J}_v+\mathbf{J}_a\right) J=λ1Js+λ2Jo+λ3(Jv+Ja) , H = λ 1 H s + λ 2 H o + λ 3 ( H v + H a ) \mathbf{H}=\lambda_1 \mathbf{H}_s+\lambda_2 \mathbf{H}_o+\lambda_3\left(\mathbf{H}_v+\mathbf{H}_a\right) H=λ1Hs+λ2Ho+λ3(Hv+Ha),以上就是碰撞速度加速度目标函数及其梯度的构建步骤。
2 软约束轨迹优化方法的优缺点
地图上保留了当前点和距离最近的障碍物及其导数,安全距离
d
d
d和
d
0
d_0
d0之类的信息,看图说明也是和人工势场法比较相似。
如果比较近(与障碍物),就会有力把无人机往外推,防止碰撞;如果比较远,力就不会那么大。
这个方法也有一些确定,从名字上就可以体现出来,因为是软约束的方法,有好几个优化目标,可能会导致函数取到局部最优解,当地形比较复杂的时候,可能会陷入局部最优解反而出不来,但是大多数情况下还是能找到最优解的。
3 Code
Code下次来仿真。