刚体动力学算法Chapter 4 刚体系统建模 Modelling Rigid Body Systems

Chapter 4 刚体系统建模 Modelling Rigid Body Systems

刚体系统是由组件部分组成的。这些组件包括刚体、关节和各种力元素。关节负责系统中的运动约束,因此我们广义地将术语“关节”解释为两个刚体之间任何可能的运动关系。在本章中,我们将从组件部分的角度来描述刚体系统的方法。特别是,我们试图描述仅由刚体和关节组成的系统。这样的描述可以称为系统模型,因为它描述的是系统本身,而不是其行为的某个方面。

系统模型的作用可以如下理解。假设我们希望在计算机上计算给定刚体系统的动力学。基本上有两种方法:

  1. 推导给定系统的符号运动方程,并将这些方程转化为计算机代码;或者
  2. 构建一个描述给定系统的系统模型,并将其作为参数提供给

(a) 基于模型的动力学计算程序;或者

(b) 基于模型的自动方程(和代码)生成器。

第二种方法更可取,部分原因是构建和修改系统模型比手工推导运动方程要简单得多,部分原因是基于模型的程序可以在被大量使用之前进行彻底的调试和优化。本书中的动力学算法是基于模型的,并且期望得到一个类似于本章描述的系统模型。方法2(b)在第10.4节中进行了讨论。
在这里插入图片描述

图4.1:两个刚体系统的连通性图

4.1 连通性

系统模型必须包含以下数据:组件部分的描述以及它们之间如何连接的描述。后者被称为系统的拓扑结构或连通性,可以通过连通性图的形式来表示,具有以下属性:

  1. 节点代表刚体。
  2. 弧代表关节。
  3. 恰好有一个节点代表固定刚体,称为基座。所有其他节点代表移动刚体。
  4. 图是无向的。
  5. 图是连通的(即,任意两个节点之间存在路径)。

术语“移动刚体”包括任何能够移动的刚体,无论实际上是否移动。因此,通过运动约束被固定的刚体也算作是移动刚体。除非另有说明,“图”和“连通图”应被理解为“无向图”。

图4.1 中显示了一些连通性图的例子。第一个是一个简单的平面机构,称为曲柄滑块连杆机构。它由三个移动刚体(B1到B3)、一个固定基座、三个旋转关节(J1 到 J3)和一个平移关节(J4)组成。刚体 B1(曲柄)的旋转使得 B3 相对于基座来回滑动。在连通性图中,每个节点代表一个刚体,每个弧代表一个关节。通常需要对至少一些节点和弧进行标记,以便清楚地显示图与原始机构之间的对应关系。此外,有时为了方便视觉识别,我们还会将基座节点与其他节点区分开来,通过给基座节点标记一个白色圆点来实现。请注意,机构中的刚体通常称为连杆。

第二个例子是一颗轨道上的卫星。这个卫星本身由两个移动刚体和一个关节组成。但是,我们不能直接从这个双刚体系统构建连通性图,因为 规则3 要求必须有一个节点表示固定基座。因此,我们添加了第三个刚体-在这种情况下是下方的行星作为固定基座。添加了第三个刚体后,规则5 要求我们还必须添加第二个关节,以便将行星连接到其他两个刚体之一。由于卫星和行星之间没有物理连接,我们使用了一种特殊的关节,它允许六个自由度运动-六自由度关节。该关节不施加运动约束,因此对系统没有物理影响。它的作用是定义额外的关节变量,以便系统整体的位置和速度变量可以从它所包含的关节的位置和速度变量中提取出来。在这种情况下,卫星具有总共 6+ns6+n_{s}6+ns 个自由度,其中 nsn_{s}ns 是卫星上关节允许的自由度。因此,我们需要 6+ns6+n_{s}6+ns 个速度变量,而这正好是两个关节定义的数量。

通过一个六自由度关节连接到固定基座的刚体称为浮动基座。在这个例子中,卫星的刚体被选择作为浮动基座,但我们也可以选择天线作为浮动基座。

4.1.1 运动树和回路 Kinematic Trees and Loops

如果一个图中任意两个节点之间存在唯一的路径,那么该图就是一个拓扑树。如果刚体系统的连通图是一个拓扑树,那么我们称该系统为运动树。图4.1中的卫星就是一个运动树,但曲柄滑块连杆机构不是。运动树具有特殊的性质,可以方便地高效计算它们的动力学。

GGG 是任意的连通图。GGG 的一个生成树,记为 GtG_{t}Gt,是 GGG 的一个子图,包含 GGG 中所有的节点,并包含 GGG 中的任意子集的弧,使得 GtG_{t}Gt 是一个拓扑树。每个连通图至少有一个生成树。如果 GGG 本身就是一棵树,那么 Gt=GG_{t}=GGt=G;否则,GtG_{t}GtGGG 的一个真子图,并且不唯一。图4.2 展示了一些生成树的例子。

一个有 nnn 个节点的拓扑树有 n−1n-1n1 条弧。因此,如果 GGG 包含 nnn_{n}nn 个节点和 nan_{a}na 条弧,那么 GtG_{t}Gt 将包含 nnn_{n}nn 个节点和 nn−1n_{n}-1nn1 弧,意味着在 GGG 中将会有 na−nn+1n_{a}-n_{n}+1nann+1 条弧不在 GtG_{t}Gt 中。这些不在拓扑树中的非树弧被称为弦。尽管 GGG 中弦的数量是固定的,但在生成树 GtG_{t}Gt 的选择下,作为弦的弧的集合会有所不同。这种变化在 图4.2 中有所体现。(我们通过虚线来区分弦和树弧。)

在这里插入图片描述

一个回路是图中的闭合路径。树中没有回路,所以图中的每个回路都必须经过至少一条弦。我们特别关注那些只经过一条弦的回路。对于图中的每条弦,总有且只有一个这样的回路,但它的路径取决于生成树的选择。例如,图4.2中的图形有两条弦和三个回路。对于任何生成树的选择,其中两个回路是单弦回路,第三个回路(未显示)经过两条弦。单弦回路的路径由它经过的一条弦以及生成树中连接该弦两个节点之间的唯一路径组成。

连通图中的每个回路都定义了一个运动回路。然而,我们只对独立的运动回路感兴趣,而不是所有的运动回路。这是在制定运动方程时必须考虑的最小回路集合。单弦回路的重要性在于它们确定了独立运动回路的集合。如果我们将刚体系统描述为具有nnn个运动回路,那么实际上意味着它具有nnn个独立的运动回路。

考虑一个由固定基座、NBN_{B}NB 个运动刚体和 NJN_{J}NJ 个关节组成的刚体系统,其中后者包括为了确保连通图连通而添加的任何 6自由度关节。设 GGG 是该系统的连通图,GtG_{t}GtGGG 的一棵生成树。这两个图都包含 NB+1N_{B}+1NB+1 个节点。然而,GGG 将包含 NJN_{J}NJ 条弧,而GtG_{t}Gt 只包含 NBN_{B}NB 条,因此系统中弦的数量为 NJ−NBN_{J}-N_{B}NJNB。系统中独立运动回路的数量 NLN_{L}NL 等于 GGG 中的弦的数量,因此由下式给出:

NL=NJ−NB(4.1) N_{L}=N_{J}-N_{B} \tag{4.1} NL=NJNB(4.1)

将该公式应用于 图4.1 中的曲柄滑块连杆机构可以发现它有一个运动回路。

如果一个系统具有运动回路,则被称为闭环系统,否则为运动树。与弦相关联的关节被称为回路关节或闭环关节,而那些属于生成树的关节被称为树关节。因此,刚体系统具有NBN_{B}NB个树关节和NLN_{L}NL个回路关节。

4.1.2 编号

刚体和关节可以用各种方式进行标识,但最方便的方法是进行编号。我们将使用一种适合描述和实施动力学算法的编号方案。根据这个方案,刚体从0到NBN_{B}NB编号,关节从1到NJN_{J}NJ编号,并且根据以下规则分配编号。

1.选择一个生成树 GtG_{t}Gt
2.将代表固定基座的节点赋予编号 0,并将该节点定义为 GtG_{t}Gt 的根节点。
3.以任意顺序将剩余节点从 1 到 NBN_{B}NB 编号,使得每个节点的编号都比它在 GtG_{t}Gt 中的父节点高。
4.将 GtG_{t}Gt 中的弧从 1 到 NBN_{B}NB 编号,使得弧iii连接节点iii和其父节点。
5.以任意顺序将所有剩余的弧从 NB+1N_{B}+1NB+1NJN_{J}NJ 编号。
6.每个刚体得到与其节点相同的编号,每个关节得到与其弧相同的编号。

这个方案称为常规编号。如果这些规则应用于无分支的运动学树(树中没有节点有多于一个子节点),那么刚体和关节将从基座开始按顺序编号,如图4.3所示。在所有其他情况下,编号不唯一。例如,图4.4显示了图4.1中曲柄滑块连杆机构的所有可能编号。这八种可能性来自于四个可能的生成树。两棵树没有分支,因此只有一种编号方式。这两种方式出现在图的左上角和右下角。另外两棵树在根节点处有分支,它们各自有三种编号方式。

4.1.3 Joint Polarity

在这里插入图片描述

每个关节连接两个刚体。我们将其中一个刚体称为关节的前驱,另一个称为关节的后继。关节本身被认为是从前驱到后继的连接。这个术语的目的有两个:定义关节和刚体变量之间的关系,并指定关节在系统中连接的方式。例如,关节上的速度定义为后继相对于前驱的速度,关节上传递的力定义为从前驱传递到后继的力。

确定关节在两个刚体之间连接的方式被称为它的极性。我们可以通过在连接图上的弧上放置箭头来确定关节的极性,箭头指向从前驱到后继。图4.5和图4.6显示了一些示例。为方便起见,我们采用以下约定关于树关节的极性:树关节的正常极性是从父节点连接到子节点。除非另有说明,所有树关节都被认为具有这种极性。

一些关节在物理上是对称的,意味着后继相对于前驱的运动与前驱相对于后继的运动在性质上是相同的。没有这个性质的关节是不对称的。改变对称关节的极性对刚体系统没有实质性影响,但改变不对称关节的极性会产生影响。例如,如果两个刚体通过球接头连接,它是对称的,那么无论哪个刚体有球,哪个有座,实际上并没有什么区别——如果我们交换它们的位置,运动仍然是相同的。相反,如果两个刚体由齿轮齿条连接,它是不对称的,那么哪个刚体具有齿条哪个具有齿轮就很重要。因此,不对称关节的极性由物理系统决定,但对称关节的极性不受限制,可以在系统模型中自由选择。

在这里插入图片描述

4.1.4 表示方法 Representation

在计算机中,有许多方法可以表示连接图。其中最简单的一种方法是使用一对数组 pppsss 记录关节的前驱和后继刚体编号,其中 p(i)p(i)p(i)s(i)s(i)s(i) 分别表示第 iii 个关节的前驱和后继刚体编号。这种方法在 图4.5 中进行了说明,图中重现了 图4.4 中的两个连接图,并添加了箭头以显示每个关节的极性。

通过 pppsss ,我们可以完整描述连接图,并可以计算其他有用的量。其中之一是父节点数组 λ\lambdaλ,用于标识跨度树中每个刚体的父节点:λ(i)\lambda(i)λ(i) 表示刚体 iii 的父节点。根据编号方案的规则3和规则4,树关节 iii 连接在刚体 iiiλ(i)\lambda(i)λ(i) 之间,且λ(i)<i\lambda(i)<iλ(i)<i,因此有

λ(i)=min⁡(p(i),s(i)).(1≤i≤NB)(4.2) \lambda(i)=\min (p(i), s(i)) . \quad\left(1 \leq i \leq N_{B}\right) \tag{4.2} λ(i)=min(p(i),s(i)).(1iNB)(4.2)

如果树关节满足 s(i)=is(i)=is(i)=i ,则它连接从父节点到子节点,并称为树中的正向关节。否则,它连接从子节点到父节点,并称为树中的反向关节。图4.5 中的树关节都是正向的(箭头指向远离根节点)。图4.6中的树关节也都是正向的,除了关节4。此图中省略了树关节编号和正向箭头:前者可以从编号方案中推导出来,后者是因为我们将“正向”视为树关节的标准方向。

还有三个量可以提供关于运动学树连接性的有用信息:对于除了基底之外的任何刚体 iiiκ(i)\kappa(i)κ(i) 是连接刚体 iii 和基底之间的所有树关节的集合;对于任何刚体 iii(包括基底),μ(i)\mu(i)μ(i) 是刚体 iii 的子节点集合,ν(i)\nu(i)ν(i) 是以刚体 iii 为起点的子树中的刚体集合。它们分别被称为支撑集、子节点集和子树集。(如果一棵树关节位于某个刚体与基底之间的路径上,则称该关节支撑该刚体。)ν(i)\nu(i)ν(i) 也可以视为由关节iii支撑的所有刚体的集合。对于 图4.6 中的跨度树,支撑集、子节点集和子树集如下:

κ(1)={1}μ(0)={1,4}ν(0)={0,1,2,…,8}κ(2)={1,2}μ(1)={2,3}ν(1)={1,2,3,5}κ(3)={1,3}μ(2)={}ν(2)={2}⋮⋮⋮κ(8)={4,8}μ(8)={}ν(8)={8} \begin{array}{lll} \kappa(1)=\{1\} & \mu(0)=\{1,4\} & \nu(0)=\{0,1,2, \ldots, 8\} \\ \kappa(2)=\{1,2\} & \mu(1)=\{2,3\} & \nu(1)=\{1,2,3,5\} \\ \kappa(3)=\{1,3\} & \mu(2)=\{\} & \nu(2)=\{2\} \\ \vdots & \vdots & \vdots \\ \kappa(8)=\{4,8\} & \mu(8)=\{\} & \nu(8)=\{8\} \end{array} κ(1)={1}κ(2)={1,2}κ(3)={1,3}κ(8)={4,8}μ(0)={1,4}μ(1)={2,3}μ(2)={}μ(8)={}ν(0)={0,1,2,,8}ν(1)={1,2,3,5}ν(2)={2}ν(8)={8}

显然,κ,λ,μ\kappa, \lambda, \muκ,λ,μν\nuν 都具有许多简单的性质,这些性质可以从它们的定义中得出。例如,μ(i)⊆ν(i),j∈κ(i)⇒i∈ν(j)\mu(i) \subseteq \nu(i), j \in \kappa(i) \Rightarrow i \in \nu(j)μ(i)ν(i),jκ(i)iν(j),等等。一个不太显然但非常有用的性质是下面的恒等式:

∑i=1NB∑j∈κ(i)(⋯ )=∑j=1NB∑i∈ν(j)(⋯ ).(4.3) \sum_{i=1}^{N_{B}} \sum_{j \in \kappa(i)}(\cdots)=\sum_{j=1}^{N_{B}} \sum_{i \in \nu(j)}(\cdots) . \tag{4.3} i=1NBjκ(i)()=j=1NBiν(j)().(4.3)

左边是对所有满足关节jjj支撑刚体iiii,ji, ji,j对的求和,而右边也是对所有满足关节jjj支撑刚体iiii,ji, ji,j对的求和,因此两者相等。

例4.1 展示了λ\lambdaλ 在动力学算法中的作用,以下是一个简单的例子。假设我们希望计算刚体系统中每个刚体的速度,已知树关节间的速度。设 vi\boldsymbol{v}_{i}vivJi\boldsymbol{v}_{\mathrm{J} i}vJi 分别为刚体 iii 的速度和过关节 iii 的速度,后者是相对于父节点的速度(即,关节是正向的)。现在,刚体 iii 的速度可以递归地定义为其父节点的速度与连接它到父节点的关节上的速度之和;因此,计算刚体速度的公式为:

vi=vλ(i)+vJi.(v0=0) \boldsymbol{v}_{i}=\boldsymbol{v}_{\lambda(i)}+\boldsymbol{v}_{\mathrm{J} i} . \quad\left(\boldsymbol{v}_{0}=\mathbf{0}\right) vi=vλ(i)+vJi.(v0=0)

计算这些速度的算法如下:

v0=0 for i=1 to NB do vi=vλ(i)+vJi end  \begin{aligned} & \boldsymbol{v}_{0}=\mathbf{0} \\ & \text { for } i=1 \text { to } N_{B} \text { do } \\ & \qquad \boldsymbol{v}_{i}=\boldsymbol{v}_{\lambda(i)}+\boldsymbol{v}_{\mathrm{J} i} \\ & \text { end } \end{aligned} v0=0 for i=1 to NB do vi=vλ(i)+vJi end 

性质 λ(i)<i\lambda(i)<iλ(i)<i 保证了在使用 vλ(i)\boldsymbol{v}_{\lambda(i)}vλ(i) 之前它总是被计算出来的。另一种表达刚体速度的方式是:

vi=∑j∈κ(i)vJj(4.4) \boldsymbol{v}_{i}=\sum_{j \in \kappa(i)} \boldsymbol{v}_{\mathrm{J} j} \tag{4.4} vi=jκ(i)vJj(4.4)

例4.2 例子2.3 介绍了一个称为刚体雅可比矩阵的矩阵,它将每个刚体的空间速度表示为关节速度变量的函数。如果 Ji\boldsymbol{J}_{i}Ji 是刚体 iii 的雅可比矩阵,则有:
vi=Jiq˙ \boldsymbol{v}_{i}=\boldsymbol{J}_{i} \dot{\boldsymbol{q}} vi=Jiq˙

方程2.18给出了 Ji\boldsymbol{J}_{i}Ji 的一个公式,适用于非分支的运动链的特殊情况,现在我们来得到一般的公式。设关节 jjj 上的速度为 Sjq˙j\boldsymbol{S}_{j} \dot{\boldsymbol{q}}_{j}Sjq˙j,并定义 ϵij\epsilon_{i j}ϵij 如下:

ϵij={1 if j∈κ(i)0 otherwise. (4.5) \epsilon_{i j}= \begin{cases}1 & \text { if } j \in \kappa(i) \\ 0 & \text { otherwise. }\end{cases} \tag{4.5} ϵij={10 if jκ(i) otherwise. (4.5)

因此,如果关节 jjj 支撑刚体 iii,则 ϵij\epsilon_{i j}ϵij 非零。刚体 iii 的速度可以表示如下:

vi=∑j∈κ(i)Sjq˙j=∑j=1NBϵijSjq˙j=[ϵi1S1⋯ϵiNBSNB][q˙1⋮q˙NB]. \boldsymbol{v}_{i}=\sum_{j \in \kappa(i)} \boldsymbol{S}_{j} \dot{\boldsymbol{q}}_{j}=\sum_{j=1}^{N_{B}} \epsilon_{i j} \boldsymbol{S}_{j} \dot{\boldsymbol{q}}_{j}=\left[\begin{array}{llll} \epsilon_{i 1} \boldsymbol{S}_{1} & \cdots & \epsilon_{i N_{B}} \boldsymbol{S}_{N_{B}} \end{array}\right]\left[\begin{array}{c} \dot{\boldsymbol{q}}_{1} \\ \vdots \\ \dot{\boldsymbol{q}}_{N_{B}} \end{array}\right] . vi=jκ(i)Sjq˙j=j=1NBϵijSjq˙j=[ϵi1S1ϵiNBSNB]q˙1q˙NB.

因此,刚体雅可比矩阵的一般公式为:

Ji=[ϵi1S1ϵi2S2⋯ϵiNBSNB](4.6) \boldsymbol{J}_{i}=\left[\begin{array}{llll} \epsilon_{i 1} \boldsymbol{S}_{1} & \epsilon_{i 2} \boldsymbol{S}_{2} & \cdots & \epsilon_{i N_{B}} \boldsymbol{S}_{N_{B}} \end{array}\right] \tag{4.6} Ji=[ϵi1S1ϵi2S2ϵiNBSNB](4.6)

4.2 几何性质 Geometry

如果两个刚体通过关节连接在一起,那么它们的相对运动是受限制的。完整描述约束需要了解两个方面:关节允许的运动以及关节相对于每个刚体的位置。前者由关节模型描述,后者由系统的几何性质描述。

刚体系统的几何模型描述了连接的位置。特别是,它描述了每个刚体上每个关节的相对位置。图4.7 显示了一个刚体系统的几何模型,包括一个固定基座 B0B_{0}B0,三个移动刚体 B1B_{1}B1B3B_{3}B3 ,三个树状关节 J1J_{1}J1J3J_{3}J3,以及一个回路关节 J4J_{4}J4。图中的虚线仅表示每个关节的连接性,而不表示其确切物理位置。

每个移动刚体内嵌有一个坐标系,为该刚体定义了一个局部坐标系统。这些坐标系标记为 F1F_{1}F1F3F_{3}F3,它们定义的坐标系称为刚体坐标或链接坐标。还有一个嵌入在基座中的坐标系 F0F_{0}F0,它为整个刚体系统定义了一个绝对、参考或基座坐标系。刚体 BiB_{i}Bi 的位置定义为 FiF_{i}Fi 相对于 F0F_{0}F0 的位置,可以用坐标变换 iX0{ }^{i} \boldsymbol{X}_{0}iX0 来表示。

在这里插入图片描述

这幅图还展示了几个坐标系,命名形式为 Fi,jF_{i, j}Fi,j,其中 iii 指代一个刚体,jjj 指代一个关节,描述了关节 jjj 在刚体 iii 中的位置。每个树状关节 (tree joint) iii连接从坐标系 Fλ(i),iF_{\lambda(i), i}Fλ(i),i 到坐标系 FiF_{i}Fi(如果关节是反向的,则相反),而每个回路关节 iii 连接从 Fp(i),iF_{p(i), i}Fp(i),iFs(i),iF_{s(i), i}Fs(i),i 。因此,总共有 NB+2NLN_{B}+2 N_{L}NB+2NL 个定位坐标系:每个树状关节一个,每个回路关节两个。坐标系 Fi,jF_{i, j}Fi,j 相对于 FiF_{i}Fi 的位置由变换 i,jXi{ }^{i, j} \boldsymbol{X}_{i}i,jXi 给出,它是一个常数。这些变换,或者可以计算它们的数据集,定义了刚体系统的几何模型。

为了方便起见,我们定义了一个树状变换(tree transforms)数组 XT\boldsymbol{X}_{\mathrm{T}}XT,使得 XT(i)=λ(i),iXλ(i)\boldsymbol{X}_{\mathrm{T}}(i)={ }^{\lambda(i), i} \boldsymbol{X}_{\lambda(i)}XT(i)=λ(i),iXλ(i)。这些是树状关节的定位变换,该数组提供了跨度树的完整几何描述。我们还定义了数组 XP\boldsymbol{X}_{\mathrm{P}}XPXS\boldsymbol{X}_{\mathrm{S}}XS,每个数组从 1 到 NLN_{L}NL 进行索引,使得对于每个回路关节 iiiXP(i−NB)=p(i),iXp(i)\boldsymbol{X}_{\mathrm{P}}\left(i-N_{B}\right)={ }^{p(i), i} \boldsymbol{X}_{p(i)}XP(iNB)=p(i),iXp(i),且 XS(i−NB)=s(i),iXs(i)\boldsymbol{X}_{\mathrm{S}}\left(i-N_{B}\right)={ }^{s(i), i} \boldsymbol{X}_{s(i)}XS(iNB)=s(i),iXs(i)。这些是回路关节的定位变换。

最后,每个关节定义了一个关节变换 XJi\boldsymbol{X}_{\mathrm{J} i}XJi,它是关节位置变量的函数。这些变换将在第4.4节中讨论。如图所示,如果 iii 是树状关节,则 $\boldsymbol{X}_{\mathrm{J} i} $ 等于 iXλ(i),i{ }^{i} \boldsymbol{X}_{\lambda(i), i}iXλ(i),i,如果 iii 是回路关节,则 XJi\boldsymbol{X}_{\mathrm{J} i}XJi 等于 s(i),iXp(i),i{ }^{s(i), i} \boldsymbol{X}_{p(i), i}s(i),iXp(i),i。例如,假设我们希望计算运动树中每个刚体的位置;也就是说,我们希望计算每个刚体的 iX0{ }^{i} \boldsymbol{X}_{0}iX0。完成此计算的算法如下所示:

在这里插入图片描述

在上述代码片段中,XJ\boldsymbol{X}_{\mathrm{J}}XJ 是一个局部变量,jtype (i)(i)(i)是描述关节iii的数据结构,jcalc 是一个从关节描述和关节位置变量向量计算关节变换的函数。关于 jtype (i)(i)(i) 和 jcalc 的详细信息可参考第4.4节。

4.3 Denavit-Hartenberg Parameters

在这里插入图片描述

通常情况下,为了指定一个坐标系相对于另一个坐标系的位置,需要六个参数。然而,对于某些刚体系统,可以以一种特定的方式定位刚体坐标系,从而只需要更少的参数。Denavit-Hartenberg(DH)方法是其中一种方法,最初由Denavit和Hartenberg(Denavit和Hartenberg,1955)描述,每个关节只需要四个参数。该方法利用了某些常见关节类型可以用空间中的一条直线来描述,并且只需要四个参数来定位该直线。

图4.8展示了如何将Denavit-Hartenberg(DH)方法应用于表示机器人手臂或类似装置的非分支运动树。该系统的重要几何特征有:

  • 基坐标系
  • nnn个关节轴
  • 嵌入在末端链接中的末端执行器坐标系

末端执行器坐标系用于机器人正确执行定位命令,但在动力学计算中并不需要。根据以下规则放置DH坐标系。

  1. z0z_0z0zn+1z_{n+1}zn+1与基坐标系和末端执行器坐标系的zzz轴对齐。

  2. z1z_1z1znz_nznnnn个关节轴对齐,使得ziz_izi与关节iii的轴对齐。(如果关节iii是移动关节,则ziz_izi可以是任何具有正确方向的轴;但在可行的情况下,选择与zi−1z_{i-1}zi1相交的轴是合理的。)

  3. xix_ixiziz_izizi+1z_{i+1}zi+1之间的公共垂线,从ziz_izi指向zi+1z_{i+1}zi+1。如果ziz_izizi+1z_{i+1}zi+1相交,则xix_ixi的方向存在180∘180^{\circ}180的歧义。如果ziz_izizi+1z_{i+1}zi+1平行,则公共垂线的位置不确定,合理的选择是与xi−1x_{i-1}xi1相交。

  4. 刚体坐标系FiF_iFi由轴xix_ixiyiy_iyiziz_izi定义,其中yiy_iyixix_ixiziz_izi派生,形成一个右手坐标系。

放置了坐标系之后,它们的相对位置由以下DH参数描述:

  • did_idi 是从xi−1x_{i-1}xi1沿ziz_izi测量到xix_ixi的距离。
  • θi\theta_iθi 是从xi−1x_{i-1}xi1xix_ixiziz_izi测量的角度。
  • aia_iai 是从zi−1z_{i-1}zi1沿xi−1x_{i-1}xi1测量到ziz_izi的距离。定义xi−1x_{i-1}xi1暗示着ai≥0a_i \geq 0ai0
  • αi\alpha_iαi 是从zi−1z_{i-1}zi1ziz_izixi−1x_{i-1}xi1测量的角度。

所有参数aia_iaiαi\alpha_iαi都是常数,但一些参数did_idiθi\theta_iθi用作关节变量。具体来说,如果关节iii是旋转关节,则θi\theta_iθi是其关节变量;如果关节iii是移动关节,则did_idi是其关节变量;如果关节iii是圆柱关节,则θi\theta_iθidid_idi都是关节变量;如果关节iii是螺旋关节,且螺距为hih_ihi,则θi\theta_iθi是关节变量,并且从xi−1x_{i-1}xi1xix_ixi的距离定义为di+hiθid_i+h_i \theta_idi+hiθi。如果每个关节都是旋转关节,则DH参数为图4.8中非分支树的以下几何模型定义:

i=1:XT(1)=xlt⁡([a10d1]T)rotx⁡(α1)rotz⁡(θ0)xlt⁡([00d0]T)i>1:XT(i)=xlt⁡([ai0di]T)rotx⁡(αi) \begin{array}{ll} i=1: & \boldsymbol{X}_{\mathrm{T}}(1)=\operatorname{xlt}\left(\left[\begin{array}{lll} a_{1} & 0 & d_{1} \end{array}\right]^{\mathrm{T}}\right) \operatorname{rotx}\left(\alpha_{1}\right) \operatorname{rotz}\left(\theta_{0}\right) \operatorname{xlt}\left(\left[\begin{array}{lll} 0 & 0 & d_{0} \end{array}\right]^{\mathrm{T}}\right) \\ i>1: & \boldsymbol{X}_{\mathrm{T}}(i)=\operatorname{xlt}\left(\left[\begin{array}{lll} a_{i} & 0 & d_{i} \end{array}\right]^{\mathrm{T}}\right) \operatorname{rotx}\left(\alpha_{i}\right) \end{array} i=1:i>1:XT(1)=xlt([a10d1]T)rotx(α1)rotz(θ0)xlt([00d0]T)XT(i)=xlt([ai0di]T)rotx(αi)

请参考表2.2中有关平移和旋转函数的定义。此外,请参阅附录中的§\S§ A.4,了解螺旋变换的详细信息。

对于具有nnn个关节的系统,需要4n+64n+64n+6个DH参数来描述,其中至少有nnn个关节变量而非参数。然而,最后四个参数仅用于定位末端执行器的坐标系,在不需要这些信息的情况下可以省略。如果基坐标系可以选择与z1z_{1}z1对齐,则还可以省略多达五个参数。

Denavit-Hartenberg方法的基本形式适用于任何刚体系统,其中每个关节可以由关节轴来描述,并且拓扑结构可以是非分支的运动树或单个封闭的运动回路。然而,该方法可以根据Khalil和Dombre(2002)的描述适应更广泛的系统类别。有关DH参数的描述可在各种教材中找到,例如Angeles(2003)和Craig(2005)。符号细节在不同的作者之间可能有所不同,特别是关于xix_{i}xiaia_{i}aiαi\alpha_{i}αi的编号。这里使用的符号遵循Springer机器人手册(Siciliano和Khatib,2008)。

在这里插入图片描述

4.4 Joint Models 关节模型

关节可以被看作是两个笛卡尔坐标系之间的运动约束,一个嵌入在关节的后继体中,另一个嵌入在其前驱中。我们将分别称这些坐标系为 FsF_{s}FsFpF_{p}Fp 。关节模型的目的是提供计算 FsF_{s}Fs 相对于 FpF_{p}Fp 的运动所需的信息,作为关节变量的函数。

关节允许的运动类型取决于其关节类型。例如,旋转关节允许围绕单个轴进行纯旋转,而球面关节允许围绕特定点进行任意旋转。因此,我们需要为每种类型的关节设计一个关节模型。关节模型的内容包括用于计算特定量的数据和公式的集合,例如:关节变换(sXp)\left({ }^{s} \boldsymbol{X}_{p}\right)(sXp)、关节速度(vJ)\left(\boldsymbol{v}_{\mathrm{J}}\right)(vJ)、关节的运动子空间(S)(\boldsymbol{S})(S)等等。表4.1显示了几种常见关节类型的公式,图4.9显示了其中一些关节的几何形状。

让我们更详细地研究一个模型。表4.1 中的第一个条目描述了一个允许 FsF_{s}Fs 相对于 FpF_{p}Fp 围绕它们的公共 zzz 轴旋转的旋转关节,并定义关节变量为旋转角度。标记为 E\boldsymbol{E}Er\boldsymbol{r}r 的两列描述了关节变换的旋转和平移分量,而变换本身由方程 sXp=rot⁡(E)xlt⁡(r){ }^{s} \boldsymbol{X}_{p}=\operatorname{rot}(\boldsymbol{E}) \operatorname{xlt}(\boldsymbol{r})sXp=rot(E)xlt(r) 给出。(有关rot、xlt和rz的定义,请参见表2.2。)注意,当关节变量 (q1)\left(q_{1}\right)(q1) 为零时,FsF_{s}FsFpF_{p}Fp 重合。我们尽可能地将此约定应用于所有关节模型。我们对所有关节模型应用的另一个约定是 S、T\boldsymbol{S}、\boldsymbol{T}ST 及相关量都是以 FsF_{s}Fs 坐标表示的。在旋转关节的情况下,这并不重要,因为 sS=pS{ }^{s} \boldsymbol{S}={ }^{p} \boldsymbol{S}sS=pS ,但对于其他关节来说这是重要的。

在设计良好的动力学模拟器中,很可能会有一个关节模型库,其中包含软件支持的每种关节类型的模型。然后,系统模型将包含系统中每个关节的关节类型描述符。我们将使用表达式 jtype (i)(i)(i) 来表示关节iii的关节类型描述符。最简单的情况下,关节类型描述符只是一个标识库中特定模型的类型代码。例如,类型代码 ‘R\mathrm{R}R’ 可以表示旋转关节,‘P\mathrm{P}P’ 表示平移关节。然而,某些关节既由类型又由一组参数来表征。例如,螺旋关节除了是螺旋形状外,还通过一个参数指定了螺旋线的螺距。刚体系统中的每个个体螺旋关节的螺距可能不同;因此,这种参数是关节的属性,而不是关节模型的属性。为了适应关节参数,我们定义关节类型描述符为一个数据结构,其中包含类型代码和所需的一组参数值(如果有)。

在这里插入图片描述

Notes: c1=cos⁡(q1),s1=sin⁡(q1)c_{1}=\cos \left(q_{1}\right), s_{1}=\sin \left(q_{1}\right)c1=cos(q1),s1=sin(q1), and sXp=rot⁡(E)xlt⁡(r){ }^{s} \boldsymbol{X}_{p}=\operatorname{rot}(\boldsymbol{E}) \operatorname{xlt}(\boldsymbol{r})sXp=rot(E)xlt(r) (see Table 2.2).

Table 4.1: Joint Model Formulae

Joint Model Calculation Routine 关节模型计算例程

第四章后面的章节描述了几种动力学算法,这些算法都使用关节模型数据。为了保持算法描述的简洁,我们假设存在一个名为jcalc 的单一函数,它可以在一次调用中计算所有请求的关节相关数据。该函数的示例调用如下所示:

[XJ,S,vJ,cJ]=jcalc⁡( type, q,α) (explicit joint constraint) [δ,T]=jcalc⁡( type ,sXp) (implicit joint constraint) [q˙]= jcalc ( type ,q,α) (for qdfn (Eq. 3.7)) \begin{array}{ll} {\left[\boldsymbol{X}_{\mathrm{J}}, \boldsymbol{S}, \boldsymbol{v}_{\mathrm{J}}, \boldsymbol{c}_{\mathrm{J}}\right]=\operatorname{jcalc}(\text { type, } \boldsymbol{q}, \boldsymbol{\alpha})} & \text { (explicit joint constraint) } \\ {[\boldsymbol{\delta}, \boldsymbol{T}]=\operatorname{jcalc}\left(\text { type },{ }^{s} \boldsymbol{X}_{p}\right)} & \text { (implicit joint constraint) } \\ {[\dot{\boldsymbol{q}}]=\text { jcalc }(\text { type }, \boldsymbol{q}, \boldsymbol{\alpha})} & \text { (for qdfn (Eq. } 3.7)) \end{array} [XJ,S,vJ,cJ]=jcalc( type, q,α)[δ,T]=jcalc( type ,sXp)[q˙]= jcalc ( type ,q,α) (explicit joint constraint)  (implicit joint constraint)  (for qdfn (Eq. 3.7))

在每种情况下,左侧的表达式是请求返回值的列表,参数类型指的是关节类型描述符(例如 jtype (i)(i)(i) )。关于使用 jcalc 计算隐式关节约束的方法将在第8章中介绍,而使用它计算 q˙\dot{\boldsymbol{q}}q˙ 的方法将在本章后面介绍。在算法列表中,我们通常假设 α=q˙\boldsymbol{\alpha}=\dot{\boldsymbol{q}}α=q˙ ,并在参数列表中使用 q˙\dot{\boldsymbol{q}}q˙ 代替 α\boldsymbol{\alpha}α。如果一个关节明确依赖于时间,那么时间必须作为列表中的最后一个参数添加进去。注意,vJ,cJ,S\boldsymbol{v}_{\mathrm{J}}, \boldsymbol{c}_{\mathrm{J}}, \boldsymbol{S}vJ,cJ,ST\boldsymbol{T}T 都用 FsF_{s}Fs 坐标表示。vJ\boldsymbol{v}_{\mathrm{J}}vJcJ\boldsymbol{c}_{\mathrm{J}}cJ 的表达式在第3.5节中给出。

示例4.4 示例4.3展示了如何计算运动树中每个刚体的坐标变换 iX0{ }^{i} \boldsymbol{X}_{0}iX0。现在我们扩展这个示例,计算 iX0{ }^{i} \boldsymbol{X}_{0}iX00vi{ }^{0} \boldsymbol{v}_{i}0vi(刚体 iii 在基座标系中的速度)。实现这一目标的算法如下所示:
0v0=0 for i=1 to NB do [XJ,vJ]=jcalc⁡(jtype⁡(i),qi,q˙i)iXλ(i)=XJXT(i) if λ(i)≠0 then iX0=iXλ(i)λ(i)X0 end 0vi=0vλ(i)+0XivJ end  \begin{aligned} & { }^{0} \boldsymbol{v}_{0}=\mathbf{0} \\ & \text { for } i=1 \text { to } N_{B} \text { do } \\ & \quad\left[\boldsymbol{X}_{\mathrm{J}}, \boldsymbol{v}_{\mathrm{J}}\right]=\operatorname{jcalc}\left(\operatorname{jtype}(i), \boldsymbol{q}_{i}, \dot{\boldsymbol{q}}_{i}\right) \\ & \quad{ }^{i} \boldsymbol{X}_{\lambda(i)}=\boldsymbol{X}_{\mathrm{J}} \boldsymbol{X}_{\mathrm{T}}(i) \\ & \quad \text { if } \lambda(i) \neq 0 \text { then } \\ & \quad{ }^{i} \boldsymbol{X}_{0}={ }^{i} \boldsymbol{X}_{\lambda(i)}{ }^{\lambda(i)} \boldsymbol{X}_{0} \\ & \quad \text { end } \\ & \quad{ }^{0} \boldsymbol{v}_{i}={ }^{0} \boldsymbol{v}_{\lambda(i)}+{ }^{0} \boldsymbol{X}_{i} \boldsymbol{v}_{\mathrm{J}} \\ & \text { end } \end{aligned} 0v0=0 for i=1 to NB do [XJ,vJ]=jcalc(jtype(i),qi,q˙i)iXλ(i)=XJXT(i) if λ(i)=0 then iX0=iXλ(i)λ(i)X0 end 0vi=0vλ(i)+0XivJ end 

这段代码片段包含一个新特性:它计算 iX0{ }^{i} \boldsymbol{X}_{0}iX0,然后使用 0Xi{ }^{0} \boldsymbol{X}_{i}0Xi 来变换 vJ\boldsymbol{v}_{\mathrm{J}}vJ。关于这种做法的解释可以在 附录AAA 中找到,在那里说明了如何直接从 X\boldsymbol{X}Xv\boldsymbol{v}v 计算X−1v\boldsymbol{X}^{-1} \boldsymbol{v}X1v,而不必先计算 X−1\boldsymbol{X}^{-1}X1。也可以直接从 X\boldsymbol{X}Xf\boldsymbol{f}f 计算 X∗f\boldsymbol{X}^{*} \boldsymbol{f}Xf(X∗)−1f\left(\boldsymbol{X}^{*}\right)^{-1} \boldsymbol{f}(X)1f,而不必先计算 X∗\boldsymbol{X}^{*}X(X∗)−1\left(\boldsymbol{X}^{*}\right)^{-1}(X)1。我们在几个动力学算法列表中利用了这个事实。

示例4.5 表4.1中平面关节和 6-DoF关节 的一个特点是,平移关节变量以 FsF_{s}Fs 坐标而不是 FpF_{p}Fp坐标表示。这样做的原因是为了使运动子空间矩阵成为一个简单的常数,以便 S∘=0,cJ=0\stackrel{\circ}{\mathbf{S}}=\mathbf{0}, \boldsymbol{c}_{\mathrm{J}}=\mathbf{0}S=0,cJ=0 ,并且 S\boldsymbol{S}S 的简单形式允许将各种优化应用于动力学算法。然而,这种安排可能存在一个很大的缺点:如果后继刚体与前驱相距很远,并开始旋转,则旋转将导致位置变量的变化速率很大(就像一个翻滚的卫星或飞行器)。在这种情况下,可以通过使用一组位置变量来表示以 FpF_{p}Fp 坐标为基础的平移,并结合一组速度变量,其中线速度以 FsF_{s}Fs 坐标表示,来兼顾两者。

pq1,pq2{ }^{p} q_{1},{ }^{p} q_{2}pq1,pq2pq3{ }^{p} q_{3}pq3 是平面关节的一组位置变量,其中 pq2{ }^{p} q_{2}pq2pq3{ }^{p} q_{3}pq3FpF_{p}Fp 坐标系中 xxxyyy 方向的平移;设 sq1,sq2{ }^{s} q_{1},{ }^{s} q_{2}sq1,sq2sq3{ }^{s} q_{3}sq3 是同一个平面关节的第二组位置变量,其中 sq2{ }^{s} q_{2}sq2sq3{ }^{s} q_{3}sq3FsF_{s}Fs 坐标系中 xxxyyy 方向的平移。第二组变量是 表4.1 中使用的变量并在 图4.9 中显示。这两组变量之间的关系是

[pq1pq2pq3]=[1000c1−s10s1c1][sq1sq2sq3]. \left[\begin{array}{l} { }^{p} q_{1} \\ { }^{p} q_{2} \\ { }^{p} q_{3} \end{array}\right]=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & c_{1} & -s_{1} \\ 0 & s_{1} & c_{1} \end{array}\right]\left[\begin{array}{c} { }^{s} q_{1} \\ { }^{s} q_{2} \\ { }^{s} q_{3} \end{array}\right] . pq1pq2pq3=1000c1s10s1c1sq1sq2sq3.

我们希望使用pqi{ }^{p} q_{i}pqi作为位置变量,但使用sq˙i{ }^{s} \dot{q}_{i}sq˙i作为速度变量。为了满足这一要求,我们需要一个公式来根据位置和速度变量计算位置变量的导数。这样的公式可以通过对上述方程进行求导获得,得到

[q˙1q˙2q˙3]=[100−q3c1−s1q2s1c1][α1α2α3], \left[\begin{array}{c} \dot{q}_{1} \\ \dot{q}_{2} \\ \dot{q}_{3} \end{array}\right]=\left[\begin{array}{ccc} 1 & 0 & 0 \\ -q_{3} & c_{1} & -s_{1} \\ q_{2} & s_{1} & c_{1} \end{array}\right]\left[\begin{array}{c} \alpha_{1} \\ \alpha_{2} \\ \alpha_{3} \end{array}\right], q˙1q˙2q˙3=1q3q20c1s10s1c1α1α2α3,

其中我们将 pqi{ }^{p} q_{i}pqipq˙i{ }^{p} \dot{q}_{i}pq˙i 的上标 ppp 省略,并用 αi\alpha_{i}αi 代替 sq˙i{ }^{s} \dot{q}_{i}sq˙i。这就是如果要求 jcalc 根据这个平面关节模型从 q\boldsymbol{q}qα\boldsymbol{\alpha}α 计算 q˙\dot{\boldsymbol{q}}q˙ 时所执行的计算。

例4.6 我们来构建一个示例,展示在图4.10中所示的齿轮齿条关节的关节模型。该关节允许继承体中的齿轮沿着前导体中的齿条滚动。齿条与FpF_{p}Fp坐标系的xxx轴平行,而齿轮位于FsF_{s}Fs坐标系的zzz轴上。齿轮具有半径rrr,这是关节模型的一个参数。该关节只允许单个自由度运动,因此只有一个关节变量,记作q1q_{1}q1,并且选择它作为齿轮的旋转角度。

在这里插入图片描述

如果当 q1=0q_{1}=0q1=0FpF_{p}FpFsF_{s}Fs 重合,则从 FpF_{p}FpFsF_{s}Fs 的变换包括沿 xxx 方向平移 rq1r q_{1}rq1,然后绕 zzz 轴旋转 q1q_{1}q1。因此,关节变换为:

XJ=sXp=[E00E][10−p×1] \boldsymbol{X}_{\mathrm{J}}={ }^{s} \boldsymbol{X}_{p}=\left[\begin{array}{cc} \boldsymbol{E} & \mathbf{0} \\ \mathbf{0} & \boldsymbol{E} \end{array}\right]\left[\begin{array}{cc} \mathbf{1} & \mathbf{0} \\ -\boldsymbol{p} \times & \mathbf{1} \end{array}\right] XJ=sXp=[E00E][1p×01]

其中,

E=rz⁡(q1)=[c1s10−s1c10001] and p=[rq100] \boldsymbol{E}=\operatorname{rz}\left(q_{1}\right)=\left[\begin{array}{ccc} c_{1} & s_{1} & 0 \\ -s_{1} & c_{1} & 0 \\ 0 & 0 & 1 \end{array}\right] \quad \text { and } \quad \boldsymbol{p}=\left[\begin{array}{c} r q_{1} \\ 0 \\ 0 \end{array}\right] E=rz(q1)=c1s10s1c10001 and p=rq100

下一步是找到 vJ\boldsymbol{v}_{\mathrm{J}}vJS\boldsymbol{S}S 的公式。我们首先将关节速度表示为关节位置和速度变量的函数。相对于 FpF_{p}Fp,齿轮以角速度ω=[00q˙1]T\boldsymbol{\omega}=\left[\begin{array}{lll}0 & 0 & \dot{q}_{1}\end{array}\right]^{\mathrm{T}}ω=[00q˙1]T绕通过 p\boldsymbol{p}p 的轴旋转,并且以线速度 p˙=[rq˙100]T\dot{\boldsymbol{p}}=\left[\begin{array}{lll}r \dot{q}_{1} & 0 & 0\end{array}\right]^{\mathrm{T}}p˙=[rq˙100]T 平移。这相当于在 FpF_{p}Fp 坐标系中表示的空间速度,即

在这里插入图片描述
由于只有一个关节变量,运动子空间矩阵的公式为 S=vJ/q˙1\boldsymbol{S}=\boldsymbol{v}_{\mathrm{J}} / \dot{q}_{1}S=vJ/q˙1。因此,

vJ=sXppv=[00q˙1rc1q˙1−rs1q˙10] and S=vJq˙1=[001rc1−rs10].  \boldsymbol{v}_{\mathrm{J}}={ }^{s} \boldsymbol{X}_{p}{ }^{p} \boldsymbol{v}=\left[\begin{array}{c} 0 \\ 0 \\ \dot{q}_{1} \\ r c_{1} \dot{q}_{1} \\ -r s_{1} \dot{q}_{1} \\ 0 \end{array}\right] \quad \text { and } \quad \boldsymbol{S}=\frac{\boldsymbol{v}_{\mathrm{J}}}{\dot{q}_{1}}=\left[\begin{array}{c} 0 \\ 0 \\ 1 \\ r c_{1} \\ -r s_{1} \\ 0 \end{array}\right] \text {. } vJ=sXppv=00q˙1rc1q˙1rs1q˙10 and S=q˙1vJ=001rc1rs10

最后,由于 S\boldsymbol{S}S 不是恒定的,我们需要一个计算速度乘积项 cJ\boldsymbol{c}_{\mathrm{J}}cJ 的公式。它可以表示为

cJ=S∘q˙1=∂S∂q1q˙12=[000−rs1q˙12−rc1q˙120] \boldsymbol{c}_{\mathrm{J}}=\stackrel{\circ}{\boldsymbol{S}} \dot{q}_{1}=\frac{\partial \boldsymbol{S}}{\partial q_{1}} \dot{q}_{1}^{2}=\left[\begin{array}{c} 0 \\ 0 \\ 0 \\ -r s_{1} \dot{q}_{1}^{2} \\ -r c_{1} \dot{q}_{1}^{2} \\ 0 \end{array}\right] cJ=Sq˙1=q1Sq˙12=000rs1q˙12rc1q˙120

(参见 方程3.42 中的 S∘\stackrel{\circ}{\mathbf{S}}S )。注意,如果 r=0r=0r=0,齿轮齿条关节简化为旋转关节。

Polarity Reversal 极性反转

如第4.1.4节所述,运动树中的关节在前向方向上意味着前驱是后继的父节点,否则就在反向方向上。如果一个关节是对称的,那么可以在定义系统模型时选择其极性,使其处于前向方向;但如果它是非对称的,则其极性由正在建模的刚体系统决定。因此,无法保证每个树关节都处于前向方向。

处理反向关节的最简单方法是将其转换为前向关节。我们可以按以下方式进行操作。首先,在关节描述符数据结构中添加一个布尔变量,如果关节是树关节且处于反向方向,则取值为true,否则为false。在给定这个额外的信息后,当这个布尔变量为true时,jcalc可以被修改为模拟相反极性的关节。实际上,jcalc 会颠倒前驱和后继的角色。为了完成这个任务,jcalc 只需要在返回之前修改一些计算得到的值。例如,它不再返回根据给定关节类型的公式计算得出的矩阵 sXp{ }^{s} \boldsymbol{X}_{p}sXp,而是返回该矩阵的逆;而且,它不再返回计算得到的 vJ\boldsymbol{v}_{\mathrm{J}}vJ 的值,而是返回 pXsvJ{ }^{p} \boldsymbol{X}_{s} \boldsymbol{v}_{\mathrm{J}}pXsvJ。修改后的返回值的完整集合如表4.2所示。这种方法的优点是,它允许动力学算法假定所有的树关节都处于前向方向。反过来,这使得算法更简单、更易于实现。

 normal  return value \begin{array}{l}\text { normal } \\ \text { return value }\end{array} normal  return value  simulating  opposite polarity \begin{array}{c}\text { simulating } \\ \text { opposite polarity }\end{array} simulating  opposite polarity 
sXp{ }^{s} \boldsymbol{X}_{p}sXppXs{ }^{p} \boldsymbol{X}_{s}pXs
SSS−pXsS-{ }^{p} \boldsymbol{X}_{s} \boldsymbol{S}pXsS
TTT−pXsT-{ }^{p} \boldsymbol{X}_{s} \boldsymbol{T}pXsT
vJ\boldsymbol{v}_{\mathrm{J}}vJ−pXsvJ-{ }^{p} \boldsymbol{X}_{s} \boldsymbol{v}_{\mathrm{J}}pXsvJ
cJc_{J}cJ−pXscJ-{ }^{p} \boldsymbol{X}_{s} \boldsymbol{c}_{\mathrm{J}}pXscJ
q˙\dot{\boldsymbol{q}}q˙q˙\dot{\boldsymbol{q}}q˙

Table 4.2: How jcalc simulates an opposite-polarity joint

4.5 Spherical Motion 球面运动

表示一般的3D旋转有许多方法(例如,参见Rooney,1977年),但主要的选择是欧拉角和欧拉参数。后者也被称为单位四元数。欧拉角的优点是位置变量是三个独立角度的集合,速度变量是它们的导数。这使得它们相对容易实现。它们的缺点是它们会出现奇点:随着方向接近奇点,表示变得越来越病态,并最终完全崩溃。欧拉参数的主要优点是它们不会出现奇点;但它们的缺点是它们是一组四个非独立的位置变量,这意味着速度变量不能是位置变量的导数。使用欧拉参数的软件必须编写成允许位置变量的数量与速度变量的数量不同。

欧拉角

我们将“欧拉角”一词分为广义和狭义两种意义。在广义上,它指的是通过一系列围绕指定坐标轴的三次旋转来表示一般旋转的任何方法。在狭义上,它指的是特殊情况,即第一次和第三次旋转围绕zzz轴进行,而第二次旋转围绕xxx轴或yyy轴进行。

每个版本的欧拉角都会出现奇点。当第一次和第三次轴对齐时,奇点出现。如果第一次和第三次旋转使用相同的坐标轴,则奇点出现在第二个角度为0∘0^{\circ}0180∘180^{\circ}180时;否则,在第二个角度为±90∘\pm 90^{\circ}±90时出现奇点。在接近奇点的附近,所表示的旋转的微小变化可以引起第一和第三角度的大变化。基于欧拉角的关节模型在奇点处无效。

让我们使用偏航角、俯仰角和滚转角构建一个球面关节模型。这些是围绕zzz轴、yyy轴和xxx轴的一系列旋转,按顺序执行。因此,相对于FpF_{p}FpFsF_{s}Fs的正确方向是先将FsF_{s}FsFpF_{p}Fp重合,然后围绕其zzz轴旋转一个量q1q_{1}q1(偏航角),然后围绕其yyy轴旋转一个量q2q_{2}q2(俯仰角),最后围绕其xxx轴旋转一个量q3q_{3}q3(滚转角)。因此,E\boldsymbol{E}E的公式为

E=rx⁡(q3)ry⁡(q2)rz⁡(q1)=[1000c3s30−s3c3][c20−s2010s20c2][c1s10−s1c10001]=[c1c2s1c2−s2c1s2s3−s1c3s1s2s3+c1c3c2s3c1s2c3+s1s3s1s2c3−c1s3c2c3]. \begin{aligned} \boldsymbol{E} & =\operatorname{rx}\left(q_{3}\right) \operatorname{ry}\left(q_{2}\right) \operatorname{rz}\left(q_{1}\right) \\ & =\left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & c_{3} & s_{3} \\ 0 & -s_{3} & c_{3} \end{array}\right]\left[\begin{array}{ccc} c_{2} & 0 & -s_{2} \\ 0 & 1 & 0 \\ s_{2} & 0 & c_{2} \end{array}\right]\left[\begin{array}{ccc} c_{1} & s_{1} & 0 \\ -s_{1} & c_{1} & 0 \\ 0 & 0 & 1 \end{array}\right] \\ & =\left[\begin{array}{ccc} c_{1} c_{2} & s_{1} c_{2} & -s_{2} \\ c_{1} s_{2} s_{3}-s_{1} c_{3} & s_{1} s_{2} s_{3}+c_{1} c_{3} & c_{2} s_{3} \\ c_{1} s_{2} c_{3}+s_{1} s_{3} & s_{1} s_{2} c_{3}-c_{1} s_{3} & c_{2} c_{3} \end{array}\right] . \end{aligned} E=rx(q3)ry(q2)rz(q1)=1000c3s30s3c3c20s2010s20c2c1s10s1c10001=c1c2c1s2s3s1c3c1s2c3+s1s3s1c2s1s2s3+c1c3s1s2c3c1s3s2c2s3c2c3.

S\boldsymbol{S}S 的公式可以通过将关节速度表示为 q˙\dot{\boldsymbol{q}}q˙ 的函数,并使用S=∂vJ/∂q˙\boldsymbol{S}=\partial \boldsymbol{v}_{\mathrm{J}} / \partial \dot{\boldsymbol{q}}S=vJ/q˙ 来获得。我们将跳过细节,只呈现最终结果:

S=[−s201c2s3c30c2c3−s30000000000]. \boldsymbol{S}=\left[\begin{array}{ccc} -s_{2} & 0 & 1 \\ c_{2} s_{3} & c_{3} & 0 \\ c_{2} c_{3} & -s_{3} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right] . S=s2c2s3c2c30000c3s3000100000.

S\boldsymbol{S}S的列包含相对于FsF_{s}Fs的坐标系,围绕其进行旋转的轴:第一列是FpF_{p}Fpzzz轴;第二列是第一次旋转后yyy轴所在的位置;第三列是前两次旋转后xxx轴所在的位置。由于S\boldsymbol{S}S不是常数,我们还需要一个关于cJ\boldsymbol{c}_{\mathrm{J}}cJ的公式。它可以表示为

cJ=S∘q˙=[−c2q˙200−s2s3q˙2+c2c3q˙3−s3q˙30−s2c3q˙2−c2s3q˙3−c3q˙30000000000][q˙1q˙2q˙3]=[−c2q˙1q˙2−s2s3q˙1q˙2+c2c3q˙1q˙3−s3q˙2q˙3−s2c3q˙1q˙2−c2s3q˙1q˙3−c3q˙2q˙3000]. \begin{aligned} \boldsymbol{c}_{\mathrm{J}}=\stackrel{\circ}{\boldsymbol{S}} \dot{\boldsymbol{q}} & =\left[\begin{array}{ccc} -c_{2} \dot{q}_{2} & 0 & 0 \\ -s_{2} s_{3} \dot{q}_{2}+c_{2} c_{3} \dot{q}_{3} & -s_{3} \dot{q}_{3} & 0 \\ -s_{2} c_{3} \dot{q}_{2}-c_{2} s_{3} \dot{q}_{3} & -c_{3} \dot{q}_{3} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right]\left[\begin{array}{c} \dot{q}_{1} \\ \dot{q}_{2} \\ \dot{q}_{3} \end{array}\right] \\ = & {\left[\begin{array}{c} -c_{2} \dot{q}_{1} \dot{q}_{2} \\ -s_{2} s_{3} \dot{q}_{1} \dot{q}_{2}+c_{2} c_{3} \dot{q}_{1} \dot{q}_{3}-s_{3} \dot{q}_{2} \dot{q}_{3} \\ -s_{2} c_{3} \dot{q}_{1} \dot{q}_{2}-c_{2} s_{3} \dot{q}_{1} \dot{q}_{3}-c_{3} \dot{q}_{2} \dot{q}_{3} \\ 0 \\ 0 \\ 0 \end{array}\right] . } \end{aligned} cJ=Sq˙==c2q˙2s2s3q˙2+c2c3q˙3s2c3q˙2c2s3q˙30000s3q˙3c3q˙3000000000q˙1q˙2q˙3c2q˙1q˙2s2s3q˙1q˙2+c2c3q˙1q˙3s3q˙2q˙3s2c3q˙1q˙2c2s3q˙1q˙3c3q˙2q˙3000.

方程4.7和4.8替代了表4.1中E\boldsymbol{E}ES\boldsymbol{S}S的公式,但T\boldsymbol{T}T的公式保持不变。因此,这是一个关节模型的例子,其中S\boldsymbol{S}S变化而T\boldsymbol{T}T保持不变。这种情况是可能的,因为范围(S)(\boldsymbol{S})(S)是一个常数,并且只需要T\boldsymbol{T}T满足范围(T)=range⁡(S)⊥(\boldsymbol{T})=\operatorname{range}(\boldsymbol{S})^{\perp}(T)=range(S)

欧拉参数 (四元数)

一般的旋转也可以通过单位向量和角度来描述。单位向量指定旋转轴,而角度指定旋转大小。欧拉参数与这种表示密切相关。如果u\boldsymbol{u}uθ\thetaθ是单位向量和角度,那么描述相同旋转的欧拉参数为

p0=cos⁡(θ/2)p1=sin⁡(θ/2)uxp2=sin⁡(θ/2)uyp3=sin⁡(θ/2)uz. \begin{aligned} & p_{0}=\cos (\theta / 2) \\ & p_{1}=\sin (\theta / 2) u_{x} \\ & p_{2}=\sin (\theta / 2) u_{y} \\ & p_{3}=\sin (\theta / 2) u_{z} . \end{aligned} p0=cos(θ/2)p1=sin(θ/2)uxp2=sin(θ/2)uyp3=sin(θ/2)uz.

They are not independent, but must satisfy the equation

p02+p12+p22+p32=1. p_{0}^{2}+p_{1}^{2}+p_{2}^{2}+p_{3}^{2}=1 . p02+p12+p22+p32=1.

球形关节必须具有三个独立的速度变量。因此,如果我们希望将四个非独立的欧拉参数用作位置变量,那么速度变量不能是位置变量的导数。速度变量的最佳选择是使用关节角速度矢量ω\boldsymbol{\omega}ωFsF_{s}Fs中的笛卡尔坐标。这个选择导致S\boldsymbol{S}S采用了表4.1中所示的简单且恒定的值。因此,我们只需要一个关于E\boldsymbol{E}E的公式和一个关于位置变量的导数的公式。它们为

E=2[p02+p12−1/2p1p2+p0p3p1p3−p0p2p1p2−p0p3p02+p22−1/2p2p3+p0p1p1p3+p0p2p2p3−p0p1p02+p32−1/2] \boldsymbol{E}=2\left[\begin{array}{ccc} p_{0}^{2}+p_{1}^{2}-1 / 2 & p_{1} p_{2}+p_{0} p_{3} & p_{1} p_{3}-p_{0} p_{2} \\ p_{1} p_{2}-p_{0} p_{3} & p_{0}^{2}+p_{2}^{2}-1 / 2 & p_{2} p_{3}+p_{0} p_{1} \\ p_{1} p_{3}+p_{0} p_{2} & p_{2} p_{3}-p_{0} p_{1} & p_{0}^{2}+p_{3}^{2}-1 / 2 \end{array}\right] E=2p02+p121/2p1p2p0p3p1p3+p0p2p1p2+p0p3p02+p221/2p2p3p0p1p1p3p0p2p2p3+p0p1p02+p321/2

and

[p˙0p˙1p˙2p˙3]=12[−p1−p2−p3p0−p3p2p3p0−p1−p2p1p0][ωxωyωz]. \left[\begin{array}{c} \dot{p}_{0} \\ \dot{p}_{1} \\ \dot{p}_{2} \\ \dot{p}_{3} \end{array}\right]=\frac{1}{2}\left[\begin{array}{ccc} -p_{1} & -p_{2} & -p_{3} \\ p_{0} & -p_{3} & p_{2} \\ p_{3} & p_{0} & -p_{1} \\ -p_{2} & p_{1} & p_{0} \end{array}\right]\left[\begin{array}{c} \omega_{x} \\ \omega_{y} \\ \omega_{z} \end{array}\right] . p˙0p˙1p˙2p˙3=21p1p0p3p2p2p3p0p1p3p2p1p0ωxωyωz.

A sp基于这些方程的球形关节模型在数值上优于基于方程4.7、4.8和4.9的模型。然而,这确实给周围的软件带来了负担,因为需要处理一个位置变量数量与速度变量数量不同的关节模型。还有三个复杂之处:位置变量的值必须通过积分方程4.13获得,而不是直接积分速度变量;零旋转由非零的欧拉参数集定义;在积分过程中需要频繁地对欧拉参数进行重新归一化。(积分过程中的截断误差导致参数不能完全满足方程4.11,从而导致线性缩放误差影响到由关节支持的整个子树。理想情况下,应该在每个积分步骤后重新归一化它们。)

 data for a  kinematic tree \begin{array}{l}\text { data for a } \\ \text { kinematic tree }\end{array} data for a  kinematic tree  extra data for a  closed-loop system \begin{array}{l}\text { extra data for a } \\ \text { closed-loop system }\end{array} extra data for a  closed-loop system 
NBN_{B}NBNLN_{L}NL
λ(1)⋯λ(NB)\lambda(1) \cdots \lambda\left(N_{B}\right)λ(1)λ(NB)p(NB+1)⋯p(NB+NL)p\left(N_{B}+1\right) \cdots p\left(N_{B}+N_{L}\right)p(NB+1)p(NB+NL)
jtype(1) ⋯\cdots jtype (NB)\left(N_{B}\right)(NB)s(NB+1)⋯s(NB+NL)s\left(N_{B}+1\right) \cdots s\left(N_{B}+N_{L}\right)s(NB+1)s(NB+NL)
XT(1)⋯XT(NB)\boldsymbol{X}_{\mathrm{T}}(1) \cdots \boldsymbol{X}_{\mathrm{T}}\left(N_{B}\right)XT(1)XT(NB)jtype(NB+1)⋯jtype⁡(NB+NL)\mathrm{jtype}\left(N_{B}+1\right) \cdots \operatorname{jtype}\left(N_{B}+N_{L}\right)jtype(NB+1)jtype(NB+NL)
I1⋯INB\boldsymbol{I}_{1} \cdots \boldsymbol{I}_{N_{B}}I1INBXP(1)⋯XP(NL)\boldsymbol{X}_{\mathrm{P}}(1) \cdots \boldsymbol{X}_{\mathrm{P}}\left(N_{L}\right)XP(1)XP(NL)
XS(1)⋯XS(NL)\boldsymbol{X}_{\mathrm{S}}(1) \cdots \boldsymbol{X}_{\mathrm{S}}\left(N_{L}\right)XS(1)XS(NL)

Table 4.3: Complete system model

完整的系统模型

表4.3展示了一个刚体系统的完整系统模型。对于一个运动链,所需的数据包括:NBN_{B}NB,父节点数组,关节类型描述数组,树形变换数组和刚体惯性数组。刚体的惯性以刚体坐标系表示,并且在该坐标系下是常数。对于闭环系统,额外需要的数据是:NLN_{L}NL,以及对于每个闭环关节,前驱和后继刚体编号、关节类型描述符、前驱和后继变换。

在最简单的情况下,软件包可能只针对一组受限制的关节类型进行设计,例如仅旋转关节和平移关节。在这种情况下,关节类型 jjj 可以简单地是一个数字类型码,我们可以将XJ\boldsymbol{X}_{\mathrm{J}}XJvJ\boldsymbol{v}_{\mathrm{J}}vJS\boldsymbol{S}S 的适当值直接硬编码到 jcalc 的源代码中,甚至直接硬编码到动力学算法的源代码中。另一个简化涉及访问关节变量向量(如 q\boldsymbol{q}qq˙\dot{\boldsymbol{q}}q˙)中的元素。在一般情况下,qi\boldsymbol{q}_{i}qiq\boldsymbol{q}q 的一个子向量,并且其在 q\boldsymbol{q}q 中的位置取决于所有前面子向量的大小。因此,访问和更新 qi\boldsymbol{q}_{i}qi 的代码应该具有一个偏移数组,每个关节一个,指示 q\boldsymbol{q}q 中的哪个元素是 qi\boldsymbol{q}_{i}qi 的第一个元素。然而,如果每个关节恰好有一个自由度,就像旋转关节和平移关节一样,那么 qi\boldsymbol{q}_{i}qi 就是 q\boldsymbol{q}q 的第 iii 个元素。

如果目标是开发一个通用的软件包,那么创建一个关节模型库,并让 jcalc 访问这个库会更有意义。其中一种特别有用的关节类型是用户自定义关节。通过使库具有可扩展性,可以提供这样的功能。或者,可以要求用户自定义关节的“参数”之一是该关节对应的用户提供的 jcalc 的地址,而内置的 jcalc 可以编程将控制权传递给用户提供的函数。

rbdl 代码注释

结构体 Joint

私有变量
  • mJointAxes : 关节的轴向
    • 这是一个数组,数组里的每一个元素为空间矢量
    • 数组里的元素个数等于关节的自由度 mDoFCount
  • mJointType:关节的类型
    • 旋转关节 Revolute、平移关节 Prismatic、绕 x 轴旋转的关节RevoluteX 等等
  • mDoFCount:关节的自由度
  • q_index:
    • 默认值为 0
  • custom_joint_index:
    • 默认值为 0

函数 jcalc

输入:机器人模型 model 、关节编号 joint_id、关节的位置 q、关节的速度 qdot

作用:计算出机器人模型里的 X_lambda[joint_id] 的 E 和 r ,以及计算 关节速度 v_J[joint_id][2]

计算方式

一堆 if else 实现根据不同的关节类型实现计算目标

以 绕z轴旋转为例子

	Scalar s, c;
	// 计算关节旋转角度的正弦值(s)和余弦值(c)。  
	// q[model.mJoints[joint_id].q_index]是关节的旋转角度。
    sincosp (q[model.mJoints[joint_id].q_index], &s, &c);

	// 计算关节的旋转矩阵。
    model.X_lambda[joint_id].E = Matrix3d (
         c * model.X_T[joint_id].E(0, 0) + s * model.X_T[joint_id].E(1, 0),
         c * model.X_T[joint_id].E(0, 1) + s * model.X_T[joint_id].E(1, 1),
         c * model.X_T[joint_id].E(0, 2) + s * model.X_T[joint_id].E(1, 2),

        -s * model.X_T[joint_id].E(0, 0) + c * model.X_T[joint_id].E(1, 0),
        -s * model.X_T[joint_id].E(0, 1) + c * model.X_T[joint_id].E(1, 1),
        -s * model.X_T[joint_id].E(0, 2) + c * model.X_T[joint_id].E(1, 2),

        model.X_T[joint_id].E(2, 0),
        model.X_T[joint_id].E(2, 1),
        model.X_T[joint_id].E(2, 2));
	// 将关节的初始位置赋值给关节的新位置。
	// 因为这是一个旋转关节,所以关节的位置没有改变,只是方向改变了。
    model.X_lambda[joint_id].r = model.X_T[joint_id].r;
	// 计算关节的角速度
    model.v_J[joint_id][2] = qdot[model.mJoints[joint_id].q_index];

绕 z 轴旋转,
Xj=rot⁡(E)xlt⁡(r) \boldsymbol{X}_{j}=\operatorname{rot}(\boldsymbol{E}) \operatorname{xlt}(\boldsymbol{r}) Xj=rot(E)xlt(r)
其中

E=rz⁡(θ)=[cs0−sc0001] \boldsymbol{E} = \operatorname{rz}(\theta)=\left[\begin{array}{ccc} c & s & 0 \\ -s & c & 0 \\ 0 & 0 & 1 \end{array}\right] E=rz(θ)=cs0sc0001
(r 表示位置,旋转关节只有旋转,没有位置变化)

因此,这里只计算 EEE 的变换,计算 第 iii 个关节的 旋转变化公式
Eλ=EJET E_{\lambda} = E_JE_T Eλ=EJET
最终计算出
iXλ(i)=XJXT(i) \quad{ }^{i} \boldsymbol{X}_{\lambda(i)}=\boldsymbol{X}_{\mathrm{J}} \boldsymbol{X}_{\mathrm{T}}(i) iXλ(i)=XJXT(i)
最后一行,是补充计算
0vi=0vλ(i)+0XivJ \quad{ }^{0} \boldsymbol{v}_{i}={ }^{0} \boldsymbol{v}_{\lambda(i)}+{ }^{0} \boldsymbol{X}_{i} \boldsymbol{v}_{\mathrm{J}} 0vi=0vλ(i)+0XivJ

函数 jcalc_XJ

返回值:Xj\boldsymbol{X}_{j}Xj

计算方式:根据不同的关节类型,用公式
Xj=rot⁡(E)xlt⁡(r) \boldsymbol{X}_{j}=\operatorname{rot}(\boldsymbol{E}) \operatorname{xlt}(\boldsymbol{r}) Xj=rot(E)xlt(r)
可以理解为是 函数 jcalc 的补充

函数 jcalc_X_lambda_S

作用:计算 iXλ(i)=XJXT(i)\quad{ }^{i} \boldsymbol{X}_{\lambda(i)}=\boldsymbol{X}_{\mathrm{J}} \boldsymbol{X}_{\mathrm{T}}(i)iXλ(i)=XJXT(i) 以及 第 iii 个关节运动子空间 Si\boldsymbol{S}_{i}Si

结构体Body

成员变量:

  • 质量大小 mMass
  • 描述刚体质心相对于刚体坐标系的三维向量 mCenterOfMass
  • 相对于刚体质心的惯性矩阵 mInertia

模型 Model

成员变量
  • lambda:vector 类型,存储了每个刚体的父刚体 id。
  • lambda_q:vector 类型,存储了直接影响当前自由度的父自由度的索引。
  • mu:这是一个二维向量,vector的vector类型,存储了给定刚体的所有子刚体的 id 。
  • dof_count:表示模型的自由度数量。
  • q_size:表示 q\mathbf{q}q 向量的大小。
  • qdot_size:表示q˙\mathbf{\dot{q}}q˙q¨\mathbf{\ddot{q}}q¨τ\mathbf{\tau}τ 向量的大小。
  • gravity:3维向量,使用笛卡尔坐标系表示重力的方向和大小。
  • v:SpatialVector 类型的vector,存储了刚体的空间速度
  • a:SpatialVector 类型的vector,存储了刚体的空间加速度。

关节相关:

  • mJoints :Joint 类型的 vector,存储了这个机器人模型的所有关节
  • S:SpatialVector 类型的vector,表示关节 i 的轴(关节运动子空间)。
  • v_J & c_J :SpatialVector 类型的vector,存储了关节的速度和加速度。
  • X_T:SpatialTransform 类型的 vector,存储了从父刚体坐标系到关节坐标系的变换。 它的表达基于父坐标系。

动力学相关:

  • c:速度相关的空间加速度

  • IA:空间惯量

  • pA:空间偏置力。

  • U、d、u :用于计算的临时变量

  • f:刚体的力,(只在 函数 InverseDynamics 使用)

  • X_lambda和X_base:这些是向量,分别存储了从父刚体到当前刚体的变换和从基座到刚体参考框架的变换。

  • mBodies:这是一个向量,存储了所有的刚体。

刚体相关:

  • X_lambda:SpatialTransform 的 vector 类型,存储 从父刚体到当前刚体的变换 。Xλ(i)=iXλ(i)X_{\lambda(i)} = {}^{i} X_{\lambda(i)}Xλ(i)=iXλ(i)
  • X_base:存储从基座标系到刚体坐标系的变换
  • mBodies:Body 类型的vector,存储所有刚体,包含 base 刚体在内。 0 … N_B 。
    • mBodies[0] - base body
    • mBodies[1] - 1st moveable body
    • mBodies[N_B] - N_Bth moveable body
成员函数
  • Model():构造函数
  • AddBody():这个函数用于将一个给定的刚体连接到模型中。
    • 添加方法:首先是要添加的父(或后继)刚体,其次是连接两个刚体的关节到原点的变换。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值