A 优化实现细节
对于轨迹生成和优化,我们选择了一个五次多项式。
平滑性项:我们通过最小化轨迹的三阶导数的积分来确保轨迹的平滑性。
Jm=∫0Ti∥pi(3)(t)∥2dt,∂Jm∂ci=2(∫0Tiβ(3)(t)β(3)(t)Tdt)ci,∂Jm∂Ti=ciTβ(3)(Ti)β(3)(Ti)Tci. \begin{aligned} J_{m} &= \int_{0}^{T_{i}} \left\| p_{i}^{(3)}(t) \right\|^{2} dt, \\ \frac{\partial J_{m}}{\partial c_{i}} &= 2 \left( \int_{0}^{T_{i}} \beta^{(3)}(t) \beta^{(3)}(t)^{\mathrm{T}} dt \right) c_{i}, \\ \frac{\partial J_{m}}{\partial T_{i}} &= c_{i}^{\mathrm{T}} \beta^{(3)}(T_{i}) \beta^{(3)}(T_{i})^{\mathrm{T}} c_{i}. \end{aligned} Jm∂ci∂Jm∂Ti∂Jm=∫0Tipi(3)(t)2dt,=2(∫0Tiβ(3)(t)β(3)(t)Tdt)ci,=ciTβ(3)(Ti)β(3)(Ti)Tci.
时间惩罚项:我们通过最小化总时间来提高轨迹的激进性。
Jt=∑i=1NTi,∂Jt∂ci=0,∂Jt∂Ti=1. J_{t} = \sum_{i=1}^{N} T_{i}, \quad \frac{\partial J_{t}}{\partial c_{i}} = 0, \quad \frac{\partial J_{t}}{\partial T_{i}} = 1. Jt=i=1∑NTi,∂ci∂Jt=0,∂Ti∂Jt=1.
动态可行性项:我们以多旋翼无人机的动力学为例,说明轨迹生成方法的细节。为了确保动态可行性,我们对轨迹的导数进行约束。这包括在整个轨迹上限制速度、加速度和急动度的大小,这一原则也适用于其他类型机器人的动力学。
Gd(t)=λvGv+λaGa+λjGj,Gv=p˙(t)2−vm2,Ga=p¨(t)2−am2,Gj=v¨(t)2−jm2. \begin{aligned} \mathcal{G}_{d}(t) &= \lambda_{v} \mathcal{G}_{v} + \lambda_{a} \mathcal{G}_{a} + \lambda_{j} \mathcal{G}_{j}, \\ \mathcal{G}_{v} &= \dot{p}(t)^{2} - v_{m}^{2}, \\ \mathcal{G}_{a} &= \ddot{p}(t)^{2} - a_{m}^{2}, \\ \mathcal{G}_{j} &= \ddot{v}(t)^{2} - j_{m}^{2}. \end{aligned} Gd(t)GvGaGj=λvGv+λaGa+λjGj,=p˙(t)2−vm2,=p¨(t)2−am2,=v¨(t)2−jm2.
其中vmv_{m}vm、ama_{m}am和jmj_{m}jm分别是允许的最大速度、加速度和急动度,λv\lambda_{v}λv、λa\lambda_{a}λa和λj\lambda_{j}λj是相应的权重。对应的梯度为:
∂Qx∂ci=2β(n)(t)p(n)(t)T,∂Qx∂Ti=2β(n+1)(t)Tcip(n)(t), \frac{\partial Q_{x}}{\partial \mathbf{c}_{i}} = 2 \beta^{(n)}(t) p^{(n)}(t)^{\mathrm{T}}, \quad \frac{\partial Q_{x}}{\partial T_{i}} = 2 \beta^{(n+1)}(t)^{\mathrm{T}} \mathbf{c}_{i} p^{(n)}(t), ∂ci∂Qx=2β(n)(t)p(n)(t)T,∂Ti∂Qx=2β(n+1)(t)Tcip(n)(t),
其中x={v,a,j}x = \{v, a, j\}x={v,a,j},n={1,2,3}n = \{1, 2, 3\}n={1,2,3}。
位置残差项:
Gp(t)=Lμ[∥p(t)−pi(t)∥2], \mathcal{G}_{p}(t) = \mathcal{L}_{\mu} \left[ \| p(t) - p^{i(t)} \|^{2} \right], Gp(t)=Lμ[∥p(t)−pi(t)∥2],
∂Gp∂ci=2L˙μβ(t)(p(t)−pi(t))T,∂Gp∂Ti=2L˙μβ(1)(t)Tcip(t), \frac{\partial \mathcal{G}_{p}}{\partial \mathbf{c}_{i}} = 2 \dot{\mathcal{L}}_{\mu} \beta(t) (p(t) - p^{i(t)})^{\mathrm{T}}, \quad \frac{\partial \mathcal{G}_{p}}{\partial T_{i}} = 2 \dot{\mathcal{L}}_{\mu} \beta^{(1)}(t)^{\mathrm{T}} \mathbf{c}_{i} p(t), ∂ci∂Gp=2L˙μβ(t)(p(t)−pi(t))T,∂Ti∂Gp=2L˙μβ(1)(t)Tcip(t),
其中L˙μ\dot{\mathcal{L}}_{\mu}L˙μ可通过对第4.1节中公式(23)求导得到。
姿态残差项:
GR(t)=Lμ[∥R(t)−1Ri(t)−I∥F2] \mathcal{G}_{R}(t) = \mathcal{L}_{\mu} \left[ \left\| R(t)^{-1} R^{i(t)} - I \right\|_{F}^{2} \right] GR(t)=Lμ[R(t)−1Ri(t)−IF2]
其中∥R(t)−1Ri(t)−I∥F2\left\| R(t)^{-1} R^{i(t)} - I \right\|_{F}^{2}R(t)−1Ri(t)−IF2可简化为tr(6−2R(t)−1Ri(t))\mathrm{tr}(6 - 2 R(t)^{-1} R^{i(t)})tr(6−2R(t)−1Ri(t))。
为了优化目的,我们使用归一化的四元数q=[w,x,y,z]T\mathbf{q} = [w, x, y, z]^{\mathrm{T}}q=[w,x,y,z]T来表示旋转。对应的旋转矩阵RRR定义为:
R=[1−2(y2+z2)2(xy−wz)2(xz+wy)2(xy+wz)1−2(x2+z2)2(yz−wx)2(xz−wy)2(yz+wx)1−2(x2+y2)]. R = \begin{bmatrix} 1 - 2(y^{2} + z^{2}) & 2(xy - wz) & 2(xz + wy) \\ 2(xy + wz) & 1 - 2(x^{2} + z^{2}) & 2(yz - wx) \\ 2(xz - wy) & 2(yz + wx) & 1 - 2(x^{2} + y^{2}) \end{bmatrix}. R=1−2(y2+z2)2(xy+wz)2(xz−wy)2(xy−wz)1−2(x2+z2)2(yz+wx)2(xz+wy)2(yz−wx)1−2(x2+y2).
特别地,R−1(t)=R(t)TR^{-1}(t) = R(t)^{\mathrm{T}}R−1(t)=R(t)T。利用这一性质,我们可以轻松计算旋转矩阵R(t)−1R(t)^{-1}R(t)−1关于四元数元素w,x,y,zw, x, y, zw,x,y,z的偏导数。最终的梯度∂GR∂ci\frac{\partial \mathcal{G}_{R}}{\partial \mathbf{c}_{i}}∂ci∂GR和∂GR∂Ti\frac{\partial \mathcal{G}_{R}}{\partial T_{i}}∂Ti∂GR可以通过结合 [Faessler et al. 2017] 中的微分平坦性属性以及链式法则得到。
避障项:
Jo(t)=Lμ[sthr−SVSDF(p)], J_{o}(t) = \mathcal{L}_{\mu} \left[ s_{\mathrm{thr}} - \mathrm{SVSDF}(\pmb{p}) \right], Jo(t)=Lμ[sthr−SVSDF(p)],
其中sthrs_{\mathrm{thr}}sthr是安全阈值。对应的梯度与 [Zhang et al. 2023] 中研究人员使用的技术一致。
为了得到一个较为准确的近似值,可以方便地计算时间积分惩罚J⋆J^{\star}J⋆及其梯度:
Ii⋆=Tiκ∑j=0kiωˉjG⋆(jκTi),J⋆=∑i=1NTi⋆,∂J⋆∂ci=∂Iˉi⋆∂G⋆∂G⋆∂ci,∂J⋆∂Ti=∂Ii⋆∂Ti+∂Ii⋆∂G⋆∂G⋆∂Ti=Iˉi⋆Ti+jκ∂Ii⋆∂G⋆∂G⋆∂t⋆. \begin{aligned} \boldsymbol{I}_{i}^{\star} &= \frac{T_{i}}{\kappa} \sum_{j=0}^{k_{i}} \bar{\omega}_{j} \mathcal{G}_{\star} \left( \frac{j}{\kappa} T_{i} \right), \quad \mathcal{J}^{\star} = \sum_{i=1}^{N} T_{i}^{\star}, \\ \frac{\partial \mathcal{J}^{\star}}{\partial \mathbf{c}_{i}} &= \frac{\partial \bar{\boldsymbol{I}}_{i}^{\star}}{\partial \mathcal{G}_{\star}} \frac{\partial \mathcal{G}_{\star}}{\partial \mathbf{c}_{i}}, \\ \frac{\partial \mathcal{J}^{\star}}{\partial T_{i}} &= \frac{\partial \mathcal{I}_{i}^{\star}}{\partial T_{i}} + \frac{\partial \mathcal{I}_{i}^{\star}}{\partial \mathcal{G}_{\star}} \frac{\partial \mathcal{G}_{\star}}{\partial T_{i}} = \frac{\bar{\boldsymbol{I}}_{i}^{\star}}{T_{i}} + \frac{j}{\kappa} \frac{\partial \mathcal{I}_{i}^{\star}}{\partial \mathcal{G}_{\star}} \frac{\partial \mathcal{G}_{\star}}{\partial t_{\star}}. \end{aligned} Ii⋆∂ci∂J⋆∂Ti∂J⋆=κTij=0∑kiωˉjG⋆(κjTi),J⋆=i=1∑NTi⋆,=∂G⋆∂Iˉi⋆∂ci∂G⋆,=∂Ti∂Ii⋆+∂G⋆∂Ii⋆∂Ti∂G⋆=TiIˉi⋆+κj∂G⋆∂Ii⋆∂t⋆∂G⋆.
其中(ωˉ0,ωˉ1,…,ωˉκi−1,ωˉκi)=(12,1,⋯ ,1,12)\left( \bar{\omega}_{0}, \bar{\omega}_{1}, \ldots, \bar{\omega}_{\kappa_{i}-1}, \bar{\omega}_{\kappa_{i}} \right) = \left( \frac{1}{2}, 1, \cdots, 1, \frac{1}{2} \right)(ωˉ0,ωˉ1,…,ωˉκi−1,ωˉκi)=(21,1,⋯,1,21)是梯形法则的系数 [Press 2007]。i=(1,2,⋯ ,N)i = \left(1, 2, \cdots, N\right)i=(1,2,⋯,N)表示第iii段,j=(1,2,⋯ ,κ)j = \left(1, 2, \cdots, \kappa \right)j=(1,2,⋯,κ)。我们将原始问题重新表述为一个无约束的非线性优化问题,并在获得上述优化变量的必要梯度后,使用 L-BFGS 算法 [Liu and Nocedal 1989] 高效求解。
B 消融研究
为了验证我们分层规划器各个组成部分的有效贡献,我们在三维场景中进行了消融实验。实验在两种不同环境中进行:一种是随机障碍物环境(Random Obstacles),另一种是墙壁上有缝隙的环境(Wall with gaps)。对于每种场景,我们比较了分层规划器的不同组成部分被替换时的规划结果。总共进行了四种类型的实验:1)使用标准A算法(在R3\mathbb{R}^3R3空间中)代替非对称A作为前端;2)没有中间端;3)在后端使用离散采样评估方法(如 [Geng et al. 2023; Wang et al. 2022] 中的方法)代替我们的连续SVSDF;4)使用包含所有组成部分的完整分层规划器。消融实验的结果如图2所示。非对称A*算法在导航复杂环境时表现出色,尤其是在存在类似墙壁的结构时。当机器人需要在由三面墙壁包围的狭窄缝隙中连续机动时(Wall with gaps),这一点尤为明显。由于CCA(连续碰撞避免)问题的高度非凸性,优化器在没有良好初始路径的情况下,通常难以找到无碰撞的轨迹。当不同形状的机器人在越来越复杂的环境中运行时,这种复杂性被进一步放大,这也凸显了前端的重要性。中间端主要帮助减轻后端优化的负担,尽管它对轨迹质量有一定影响,但不如前端显著。当障碍物位于初始轨迹形成的SV(扫掠体积)内部时,基于SVSDF的后端是不可或缺的。其对CCA的有效梯度贡献帮助SV高效地避开障碍物。从前端到中间端,再到后端,每个部分在任意形状机器人的连续无碰撞规划中都发挥着不可或缺的作用,每个组成部分对整个框架的有效性都至关重要。
我们运行了蒙特卡洛实验,计算机器人在其最终轨迹上沿最近障碍物的平均最小有符号距离。这一指标用于评估避障的有效性,具体是SV与障碍物之间的最小距离,正值表示有间隙,负值表示碰撞。我们还评估了CCA成功率,定义为SVSDF返回正值的试验次数与总试验次数的比率。
消融研究的统计结果如图2所示,清晰地表明,使用标准A算法作为前端无法在复杂场景中实现无碰撞规划。这并不令人意外,因为标准A算法没有考虑机器人的形状和姿态信息。因此,它生成的前端轨迹质量较低,使得优化器难以在这一高度非凸问题中找到满足约束的轨迹。移除中间端对轨迹质量有一定影响,但由于我们仍然有基于SVSDF的后端和考虑任意形状机器人的前端,轨迹的避障性能优于使用标准A*替换前端的情况。此时,我们观察到后端优化的迭代次数增加,轨迹质量略有下降。因此,我们可以得出结论,中间端的主要作用是减轻后端的优化负担。当使用离散后端时,在复杂场景中实现CCA也具有挑战性。离散后端评估方法可能导致优化器忽略可能导致碰撞的障碍物,从而生成穿过障碍物的轨迹。此外,由于缺乏连续碰撞评估,机器人的复杂形状可能导致避障梯度出现振荡,从而导致CCA的成功率相对较低。
C GSIP收敛性证明
在研究广义半无限规划(GSIP)时,离散化方法的收敛性起着关键作用。我们的证明基于以下原理 [López and Still 2007]:
引理:在连续函数和Y(x)Y(x)Y(x)的紧致性假设下,如果YkY_kYk和YYY之间的Hausdorff距离随着k→∞k \to \inftyk→∞趋于零,则近似解序列将收敛到GSIP的解。
证明:在二维情况下,YYY是一个圆,而在三维情况下,YYY是一个球体。离散集合YkY_kYk由在给定分辨率ϵres\epsilon_{\text{res}}ϵres下均匀分布在圆(二维)或球体(三维)内的离散点组成。YkY_kYk和YYY之间的Hausdorff距离小于半径rkr_krk和rSVSDFr_{\text{SVSDF}}rSVSDF之间的差值。我们的算法可以实现任意精度要求ϵ\epsilonϵ,即rSVSDF−r∗≤ϵr_{\text{SVSDF}} - r^* \leq \epsilonrSVSDF−r∗≤ϵ,其中r∗r^*r∗是最终迭代输出。只要在最终迭代中,圆周/球面上任意两个均匀分布的采样点之间的距离小于ϵ/2\epsilon/2ϵ/2即可。鉴于评估函数ggg的空间连续性,从离散点采样评估的g∗g^*g∗与实际非离散LP(⋅\cdot⋅)解之间的差异在ϵ\epsilonϵ以内。因此,从离散方法得到的解rkr_krk与真实值rSVSDFr_{\text{SVSDF}}rSVSDF之间的差异也在ϵ\epsilonϵ以内。以二维为例进行说明,我们只需设计一个均匀分布点的采样策略,具体是一个均匀分布分辨率ϵres\epsilon_{\text{res}}ϵres的更新策略。简单地应用ϵresk=max(ϵresk−1/2,ϵ/2)\epsilon_{\text{res}}^k = \max(\epsilon_{\text{res}}^{k-1}/2, \epsilon/2)ϵresk=max(ϵresk−1/2,ϵ/2)即可确保采样密度增加并满足收敛时的精度要求。在三维场景中,我们可以设计类似的采样策略。
因此,我们基于离散化方法的算法(第3.2节)具有理论上的收敛保证,该方法受到 [Blankenship and Falk 1976; Remez 1962] 的启发。