DYNAMIC MODELING AND CONTROL OF MULTIBODY SYSTEMS USING DUAL QUATERNIONS
-
摘要: 空间机械臂的在轨作业是当前空间在轨服务中应用最为广泛的技术之一, 然而, 机械臂在操作过程中漂浮基座与臂体之间的位姿耦合效应极为显著, 这给控制系统设计带来了新挑战. 针对多刚体系统的位姿一体化建模与控制问题, 文章改进了基于对偶四元数的位姿一体化建模和控制方法, 使之可以应用于多刚体系统. 该方法不仅能够精确描述复杂的力学关系, 还能够在统一的数学框架中有效地处理位姿耦合问题, 这为后续设计姿轨一体化的控制系统提供了极大的便利. 首先, 基于铰链模型建立了对偶四元数形式的铰链和臂体之间的速度和加速度递推关系, 然后, 利用铰链和臂体间力-力矩传递关系建立了递推形式的逆向动力学方程, 为了便于控制系统设计与分析, 随后推导建立了矩阵形式的位姿一体化的正向动力学方程. 然后针对推进器和控制力矩陀螺讨论了具体的执行机构的动力学建模问题, 并分别讨论了机械臂和漂浮基座的位姿一体化控制问题. 最后, 对一个6自由度机械臂和漂浮基座的组合体进行了动力学建模和控制仿真, 动力学仿真的结果证明了所提动力学建模方法的正确性, 控制仿真表明控制系统能较快地抵消机械臂运动对基座产生的干扰力和干扰力矩, 证明了所提控制方法的有效性和可行性.Abstract: The operation of space robotic arms in orbit is one of the most extensively applied technologies in current space in-orbit services. However, the significant coupling effect of position and attitude between the floating base and the arm during operations presents new challenges for the design of control systems. To address the integrated modeling and control problems of position and attitude in multi-rigid body systems, this paper improves the dual quaternion-based integrated modeling and control method, making it applicable to multi-rigid body systems. This method not only accurately describes complex mechanical relationships but also effectively manages the coupling problems of position and attitude within a unified mathematical framework, greatly facilitating the subsequent design of integrated control systems for position and trajectory. Initially, leveraging the hinge model, the paper establishes recursive relationships for velocity and acceleration between the hinges and the arm in dual quaternion form. Then, using the force-torque transmission relationship between the hinges and the arm, a recursive form of the inverse dynamics equation is established to ease the design and analysis of control systems. Following this, a matrix form of the integrated position and attitude forward dynamics equation is derived. The paper then discusses the dynamics modeling issues related to actuators, such as thrusters and control moment gyroscopes. It addresses the integrated control problems of position and attitude for both the robotic arm and the floating base. Finally, dynamics modeling and control simulations for a composite entity comprising a six-degree-of-freedom robotic arm and a floating base are conducted. The dynamics simulation results confirm the correctness of the proposed dynamics modeling method, and the control simulation demonstrates that the control system can quickly counteract the disturbance forces and torques generated by the movement of the robotic arm on the base, proving the effectiveness and feasibility of the proposed control method.
-
引 言
在过去的几十年中, 多刚体系统在空间机械臂的建模和控制方面一直是航天工程领域的重要研究课题. 空间机械臂的应用场景广泛, 例如捕获空间碎片[1]、修理故障航天器[2]以及执行复杂的在轨组装和维护任务[3]. 这些操作在微重力环境下进行, 但空间机械臂的基座和机械臂之间存在位姿耦合, 加之空间环境中各种不可预测的扰动因素, 使得控制问题极具挑战性. 在这种背景下, 空间机器人系统的建模和控制需求迫切.
空间机械臂的传统控制策略面临着一定困难. 这些困难源于旋转群(SO(3))的拓扑结构与欧几里得空间的本质不同, 使得位置和姿态的描述及其控制更加复杂. 在处理空间机械臂的旋转和平移运动时, 传统方法通常将姿态和位置分开考虑[4-5], 此种分离的控制策略不仅增加了计算复杂性, 而且在某些情况下可能导致系统性能和稳定性的降低. 例如, 在执行精密操作或快速响应任务时, 姿态与位置之间的相互影响可能导致控制精度下降或响应迟缓[6].
对偶四元数作为一种有效的数学工具, 近年来在动力学建模领域显示出了巨大的潜力. 对偶四元数结合了四元数和对偶数的特点, 它的实部和对偶部分分别可以表示空间变换中的旋转和平移, 使其可以在统一的数学框架内表示SE(3)群的元素. 这种特性为描述和控制刚体运动提供了一种更加简洁和直观的方法. 与传统的欧拉角或轴角表示法相比, 对偶四元数提供了一种唯一且连续的空间运动表示方法. 它继承了四元数的优点, 克服了欧拉角中的万向锁现象, 并解决了轴角表示法存在的多解和不连续性问题. 此外, 对偶四元数在计算上显示出高效性, 因为它允许直接且无需转换为其他表示形式即可进行旋转和平移的组合操作[7], 这在机器人学、计算图形学和航空航天领域的应用中尤为重要. Dooley 等[8] 首次探索了将对偶四元数作为广义坐标用于建模的可能性. 随后, Brodsky等[9] 提出了一种基于对偶四元数的方法实现单刚体的动力学建模, 并发展了基于牛顿-欧拉法和拉格朗日法的动力学框架. Chen等[10]提出了一种基于对偶四元数的机械臂逆运动学问题的解决方法. 与传统的方法相比, 该方法需要更少的计算量, 并避免了奇异性问题.
在国内, 杨一岱等[11]利用对偶四元数, 推导了挠性航天器的姿轨一体化动力学模型并提出了一种自适应位置姿态跟踪控制器. 该模型可以完整地描述航天器平动、转动与挠性附件振动三者之间的关联耦合作用, 并使用自适应控制器补偿了挠性附件的耦合作用, 提高了控制系统的稳定性. 董宏洋[12]探讨了多航天器编队任务中的6自由度相对运动跟踪控制问题. 他采用对偶四元数描述航天器间的位姿一体化相对运动, 并设计了多种控制方法, 包括处理质量与惯量不确定性的自适应控制方法, 以及针对执行机构故障的容错控制方法. 张洪珠[13]进一步地提出了基于对偶四元数的滑模变结构控制器, 考虑了模型不确定性和外部干扰, 对比了不同控制算法的优势与不足, 并通过仿真验证了所提出控制器的有效性.
上述研究的重点多数是对偶四元数在单刚体动力学与控制中的应用. 在多体动力学领域, Benić等[14] 使用对偶四元数推导了双自由度机械臂的欧拉-拉格朗日形式的正向运动学和动力学方程. 这个方法仅仅对特定数目的刚体有效, 较难拓展至任意数量的多刚体系统. Antonello等[15] 和 Valverde等[16] 提出了一种牛顿-欧拉法的建模方法, 考虑了不同类型的关节, 但这一方法在多体系统结构变化时缺乏灵活性, 需要重新推导方程, 这在分析和计算上是不利的.
本文着重于利用对偶四元数进行多刚体系统的建模, 推导了链式多刚体位姿一体化建模问题和多刚体系统一体化跟踪控制问题. 针对多自由度冗余空间机械臂的位姿一体化建模问题, 首先推导了链式系统中的速度递推关系, 然后基于单刚体的动力学方程, 推导了链式系统的逆动力学递推方程和正向动力学的矩阵形式, 最后讨论了具体的执行机构建模. 通过这一过程, 本文建立了一个基于对偶四元数的递推形式的机械臂位姿一体化模型. 该模型不仅考虑了机械臂的动力学特性, 还涵盖了实际执行机构的建模细节. 最终, 通过数值仿真, 本研究验证了所提出方法的有效性和正确性.
1. 数学基础
1.1 对偶四元数
对偶四元数集具有如下形式
$$ {\boldsymbol{HD}} = \{ {\boldsymbol{\hat \sigma }} = {\boldsymbol{p}} + {{\varepsilon }}{\boldsymbol{q}},{\boldsymbol{p}},{\boldsymbol{q}} \in {\boldsymbol{H}}\} $$ (1) 其中, $ {\boldsymbol{p}}和{\boldsymbol{q}} $为四元数, 算子${{\varepsilon }}$满足${{\varepsilon }} \ne 0,{{{\varepsilon }}^2} = 0$. 单位对偶四元数可以表示空间任意6自由度变换, 已有较多文献[16-19]对对偶四元数做了详细介绍, 表1总结了本文用到的对偶四元数运算法则.
表 1 部分对偶四元数运算法则Table 1. Dual quaternion operation rulesOperation Formula expression multiplication ${\hat \sigma _1} \circ {\hat \sigma _2} = {{\boldsymbol{p}}_1} \circ {{\boldsymbol{p}}_2} + {{\varepsilon }}\left( {{{\boldsymbol{p}}_1} \circ {{\boldsymbol{q}}_2} + {{\boldsymbol{q}}_1} \circ {{\boldsymbol{p}}_2}} \right)$ conjugate ${\hat \sigma ^*} = {{\boldsymbol{p}}^*} + {{\varepsilon }}{{\boldsymbol{q}}^*}$ cross product ${\hat \sigma _1} \times {\hat \sigma _2} = {{\boldsymbol{p}}_1} \times {{\boldsymbol{p}}_2} + {{\varepsilon }}\left( {{{\boldsymbol{p}}_1} \times {{\boldsymbol{q}}_2} + {{\boldsymbol{p}}_2} \times {{\boldsymbol{q}}_1}} \right)$ frame transformation ${\boldsymbol{\hat r}}_b^b = {\boldsymbol{\hat q}}_{ba}^{a*} \circ {\boldsymbol{\hat r}}_a^a \circ {\boldsymbol{\hat q}}_{ba}^a$ logarithm $ \log \left( {{\boldsymbol{p}} + {{\varepsilon }}{\boldsymbol{q}}} \right) = \log {\boldsymbol{p}} + {{\varepsilon }}\left( {{{\boldsymbol{p}}^{ - 1}} \circ {\boldsymbol{q}}} \right) $ 其中对偶四元数的对数用到了四元数的逆, 其定义为[13]
$$ {{\boldsymbol{p}}^{ - 1}} = \frac{{{{\boldsymbol{p}}^ * }}}{{{{\left\| {\boldsymbol{p}} \right\|}^2}}} $$ (2) 对偶四元数的乘法又可以用矩阵形式表示为
$$ \left. \begin{split} & {{{\hat {\boldsymbol{\varLambda}} }}_{{1}}^{{ + }} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{1L}}}&{{\boldsymbol{0}}} \\ {{{\boldsymbol{q}}_{1L}}}&{{{\boldsymbol{p}}_{1L}}} \end{array}} \right] \in {\mathbb{R}^{8 \times 8}}} \\ & {{{\hat {\boldsymbol{\varLambda}} }}_{{2}}^{{ - }} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{1R}}}&{{\boldsymbol{0}}} \\ {{{\boldsymbol{q}}_{1R}}}&{{{\boldsymbol{p}}_{1R}}} \end{array}} \right] \in {\mathbb{R}^{8 \times 8}}} \\ & {{{\boldsymbol{P}}_L} = \left[ {\begin{array}{*{20}{c}} {{p_0}}&{ - {{{{\bar {\boldsymbol{p}}}}}^{\text{T}}}} \\ {{{\bar {\boldsymbol{p}}}}}&{{p_0}{{\boldsymbol{I}}_3} + {{\tilde {\bar {\boldsymbol{p}}}}}} \end{array}} \right] \in {\mathbb{R}^{4 \times 4}}} \\ & {{{\boldsymbol{Q}}_R} = \left[ {\begin{array}{*{20}{c}} {{q_0}}&{ - {{{{\bar {\boldsymbol{q}}}}}^{\text{T}}}} \\ {{{\bar {\boldsymbol{q}}}}}&{{q_0}{{\boldsymbol{I}}_3} - {{\tilde {\bar {\boldsymbol{q}}}}}} \end{array}} \right] \in {\mathbb{R}^{4 \times 4}}} \end{split} \right\}$$ (3) 其中, $\tilde \cdot $表示矢量的反对称阵形式, $\bar \cdot $表示取一个四元数的矢量部分.
1.2 对偶四元数形式的运动学和动力学
单位对偶四元数可以表示空间中的一个6自由度变换, 通常为了表示此变换它被写作
$$ {{\boldsymbol{\hat q}}_{ba}} = {{\boldsymbol{q}}_{ba}} + \frac{{{\varepsilon }}}{{\text{2}}}{{\boldsymbol{q}}_{ba}} \circ {\boldsymbol{r}}_{ba}^b $$ (4) 式中, $ {{\boldsymbol{q}}_{ba}} $为表示两个坐标系相对姿态的四元数, 下标ba表示坐标系b相对于坐标系a, 上标表示矢量投影所在的坐标系. $ {\boldsymbol{r}}_{ba}^b $是一个实部为0的四元数.
两个坐标系之间的对偶速度的定义为
$$ {\boldsymbol{\hat \omega }}_{ba}^b = {\boldsymbol{\omega }}_{ba}^b + {{\varepsilon }}\left( {{\boldsymbol{\dot r}}_{ba}^b + {\boldsymbol{\omega }}_{ba}^b \times {\boldsymbol{r}}_{ba}^b} \right) $$ (5) 刚体相对于任意一点C的动力学方程可以写作[18]
$$ {{{\hat {\boldsymbol{F}}}}_C} = {{{\hat {\boldsymbol{M}}}}_C}\frac{\partial }{{\partial t}}{{{\hat {\boldsymbol{\omega}} }}_C} + {{{\hat {\boldsymbol{\omega}} }}_C} \times {{{\hat {\boldsymbol{M}}}}_C}{{{\hat {\boldsymbol{\omega}} }}_C} $$ (6) 其中, ${{{\hat {\boldsymbol{F}}}}_C} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_C}}&{{{\boldsymbol{T}}_C}} \end{array}} \right]^{\text{T}}}$为对偶力, $\dfrac{\partial }{{\partial t}}$表示相对导数, ${{{\hat {\boldsymbol{M}}}}_C} = \left[ {\begin{array}{*{20}{c}} { - {{\boldsymbol{S}}_C}}&{m{\boldsymbol{I}}} \\ {{{\boldsymbol{J}}_C}}&{{{\boldsymbol{S}}_C}} \end{array}} \right]$为对偶质量阵.
2. 递推形式的对偶四元数动力学方程
2.1 坐标系定义
本节使用对偶四元数, 完成任意串联机械臂的正逆动力学的递推, 首先给出坐标系定义如下.
(1)链式多刚体系统(图1)由N个刚体组成, 从端部向基座编号, 依次为1, 2, ···, N.
(2)第k个刚体的质心记为${O_k}$, 第k个铰链的内节点记为${d_k}$,外节点记为${t_k}$.
(3)惯性坐标系为$ {\varTheta _O} $, 第k个刚体的随体坐标系${\varTheta _{{d_k}}}$位于内接铰点${d_k}$处, 第k铰链的随体坐标系${\varTheta _{{t_k}}}$位于第k + 1体的外接铰点${t_{k + 1}}$处.
2.2 逆向动力学递推
铰链的相对运动可以通过铰链的映射矩阵来表示, 该矩阵建立了铰链的广义速度与铰链的空间速度的关系
$$ {{\hat {\boldsymbol{\omega}} }} = {\boldsymbol{\hat H}}\dot {\boldsymbol{\theta}} $$ (7) 其中铰链映射矩阵$ {\boldsymbol{\hat H}} \in {\mathbb{R}^{8 \times n}} $, n为铰链自由度, $\dot \theta $为铰链的广义速度. 表2给出了一些常见的铰链所对应的映射矩阵和广义速度关系.
表 2 常见铰链的映射矩阵和广义速度Table 2. Mapping matrix and generalized velocity for common hingesHinge type Map matrix Generalized velocity revolute $\left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 3}}} \\ 1 \\ {{{{\boldsymbol{0}}}_{1 \times 4}}} \end{array}} \right)$ $ \dot \theta $ prismatic $\left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 7}}} \\ 1 \end{array}} \right)$ $\dot d$ spherical $\left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 3}}} \\ {{{\boldsymbol{I}}_3}} \\ {{{{\boldsymbol{0}}}_{4 \times 3}}} \end{array}} \right)$ $\left( {\begin{array}{*{20}{c}} {{\omega _x}} \\ {{\omega _y}} \\ {{\omega _z}} \end{array}} \right)$ 6-DoF joint $ \left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 3}}}&{{{{\boldsymbol{0}}}_{1 \times 3}}} \\ {{{\boldsymbol{I}}_3}}&{{{{\boldsymbol{0}}}_{3 \times 3}}} \\ {{{{\boldsymbol{0}}}_{1 \times 3}}}&{{{{\boldsymbol{0}}}_{1 \times 3}}} \\ {{{{\boldsymbol{0}}}_{3 \times 3}}}&{{{\boldsymbol{I}}_3}} \end{array}} \right) $ $\left( {\begin{array}{*{20}{c}} {{\omega _x}} \\ {{\omega _y}} \\ {{\omega _z}} \\ {{v_x}} \\ {{v_y}} \\ {{v_z}} \end{array}} \right)$ 第k + 1刚体到第k刚体间的速度递推可以由两部分组合: 第k + 1刚体体坐标系到铰链k中铰坐标系的速度传递, 铰链k中铰坐标系到第k刚体体坐标系的传递, 见图2.
设第k + 1刚体内接铰点对偶速度为${{\boldsymbol{\hat \omega }}_{{{{d}}_{{{k + 1}}}}}}$, 外接铰点对偶速度为${{\boldsymbol{\hat \omega }}_{{t_k}}}$, 第k + 1刚体内的速度传递可以表示为
$$ {{\boldsymbol{\hat \omega }}_{{t_k}}} = {\boldsymbol{\hat q}}_{{t_k}{d_{k + 1}}}^ * \circ {{\boldsymbol{\hat \omega }}_{{d_{k + 1}}}} \circ {{\boldsymbol{\hat q}}_{{t_k}{d_{k + 1}}}} $$ (8) 对偶四元数$ {\hat {\boldsymbol{q}}_{{t_k}{d_{k + 1}}}} = {{\boldsymbol{q}}_{{t_k}{d_{k + 1}}}} + \dfrac{\varepsilon }{2}{\boldsymbol{p}}_{{t_k}{d_{k + 1}}}^{{t_k}} \circ {{\boldsymbol{q}}_{{t_k}{d_{k + 1}}}} $, $ {\boldsymbol{p}}_{{t_k}{d_{k + 1}}}^{{t_k}} $表示 ${\varTheta _{{d_{k + 1}}}}$ 指向 ${\varTheta _{{t_k}}}$ 的矢量.
第k铰链间的速度传递可以表示为
$$ {{\boldsymbol{\hat \omega }}_{{d_k}}} = {\boldsymbol{\hat q}}_{{d_k}{t_k}}^ * \circ {{\boldsymbol{\hat \omega }}_{{t_k}}} \circ {{\boldsymbol{\hat q}}_{{d_k}{t_k}}} + {{\boldsymbol{\hat H}}_k}{\dot \theta _k} $$ (9) 根据式(8)和式(9)可以得到
$$ \begin{split} & {{{\boldsymbol{\hat \omega }}}_{{d_k}}} = {\boldsymbol{\hat q}}_{{d_k}{t_k}}^ * \circ \left( {{\boldsymbol{\hat q}}_{{t_k}{d_{k + 1}}}^ * \circ {{{\boldsymbol{\hat \omega }}}_{{d_{k + 1}}}} \circ {{{\boldsymbol{\hat q}}}_{{t_k}{d_{k + 1}}}}} \right) \circ {{{\boldsymbol{\hat q}}}_{{d_k}{t_k}}} + {{{\boldsymbol{\hat H}}}_k}{{\dot \theta }_k}= \\ &\qquad {\boldsymbol{\hat q}}_{{d_k}{d_{k + 1}}}^ * \circ {{{\boldsymbol{\hat \omega }}}_{{d_{k + 1}}}} \circ {{{\boldsymbol{\hat q}}}_{{d_k}{d_{k + 1}}}} + {{{\boldsymbol{\hat H}}}_k}{{\dot \theta }_k}\\[-1pt] \end{split}$$ (10) 为了方便起见, 在后文中统一用$ {{\boldsymbol{\hat \omega }}_k} $表示第k刚体的内接铰点速度, $ {{\boldsymbol{\hat q}}_{kk + 1}} $表示第k + 1刚体体坐标系到第k刚体体坐标系的对偶四元数. 于是
$$ {{\boldsymbol{\hat \omega }}_k} = {\boldsymbol{\hat q}}_{kk + 1}^ * \circ {{\boldsymbol{\hat \omega }}_{k + 1}} \circ {{\boldsymbol{\hat q}}_{kk + 1}} + {{\boldsymbol{\hat H}}_k}{\dot \theta _k} $$ (11) 第k刚体相对惯性系的对偶速度可以通过以下递推得到
$$\left. \begin{split} & {{{\boldsymbol{\hat \omega }}}_{n + 1}} = {0} \\ &\quad {{\mathrm{for}}}\quad {k = n: - 1:1} \\ &\quad {{{{\boldsymbol{\hat \omega }}}_k} = {\boldsymbol{\hat q}}_{kk + 1}^ * \circ {{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ {{{\boldsymbol{\hat q}}}_{kk + 1}} + {{{\boldsymbol{\hat H}}}_k}{{\dot \theta }_k}} \\ &\quad {{\mathrm{end}}} \end{split}\right\} $$ (12) 对式(8)求相对导数得到
$$ \begin{split} & \frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_k} = \left( {\frac{{{\partial ^k}}}{{\partial t}}{\boldsymbol{\hat q}}_{kk + 1}^ * } \right) \circ {{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ {{{\boldsymbol{\hat q}}}_{kk + 1}} + {\boldsymbol{\hat q}}_{kk + 1}^ * \circ \left( {\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_{k + 1}}} \right) \circ \\ &\qquad {{{\boldsymbol{\hat q}}}_{kk + 1}}+ {\boldsymbol{\hat q}}_{kk + 1}^ * \circ {{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ \left( {\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat q}}}_{kk + 1}}} \right) + \left( {\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat H}}}_k}} \right){{\dot \theta }_k} + {{{\boldsymbol{\hat H}}}_k}{{\ddot \theta }_k} \end{split} $$ (13) 由于
$$ \frac{{{\partial ^k}}}{{\partial t}}{{\boldsymbol{\hat \omega }}_{k + 1}} = \frac{{{\partial ^{k + 1}}}}{{\partial t}}{{\boldsymbol{\hat \omega }}_{k + 1}} + {{\boldsymbol{\hat \omega }}_{k + 1}} \times {{\boldsymbol{\hat \omega }}_{kk + 1}} $$ (14) 令
$$ \begin{split} & {{{\boldsymbol{\hat \alpha }}}_k} = \left( {\frac{{{\partial ^k}}}{{\partial t}}{\boldsymbol{\hat q}}_{kk + 1}^ * } \right) \circ {{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ {{{\boldsymbol{\hat q}}}_{kk + 1}} + {\boldsymbol{\hat q}}_{kk + 1}^ * \circ {{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ \left( {\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat q}}}_{kk + 1}}} \right)+ \\ &\qquad {\boldsymbol{\hat q}}_{kk + 1}^ * \circ \left( {{{{\boldsymbol{\hat \omega }}}_{k + 1}} \times {{{\boldsymbol{\hat \omega }}}_{kk + 1}}} \right) \circ {{{\boldsymbol{\hat q}}}_{kk + 1}} + \left( {\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat H}}}_k}} \right){{\dot \theta }_k} \end{split} $$ (15) 于是
$$ \frac{{{\partial ^k}}}{{\partial t}}{{\boldsymbol{\hat \omega }}_k} = {\boldsymbol{\hat q}}_{kk + 1}^ * \circ \frac{{{\partial ^{k + 1}}}}{{\partial t}}{{\boldsymbol{\hat \omega }}_{k + 1}} \circ {{\boldsymbol{\hat q}}_{kk + 1}} + {{\boldsymbol{\hat \alpha }}_k} + {{\boldsymbol{\hat H}}_k}{\ddot \theta _k} $$ (16) 根据文献[18], 刚体的6自由度动力学方程可以写作
$$ {{\boldsymbol{\hat f}}_B} = {{\boldsymbol{\hat M}}_B}\frac{\partial }{{\partial t}}{{\boldsymbol{\hat \omega }}_B} + {{\boldsymbol{\hat \omega }}_B} \times {{\boldsymbol{\hat M}}_B}{{\boldsymbol{\hat \omega }}_B} $$ (17) 将第k刚体内接铰点受到的对偶力记为${{\boldsymbol{\hat f}}_{{d_k}}}$, 外接铰点受到的对偶力记为${{\boldsymbol{\hat f}}_{{t_{k - 1}}}}$, 根据对偶形式的刚体动力学方程
$$ {{\boldsymbol{\hat f}}_{{d_k}}} - {\boldsymbol{\hat q}}_{{d_k}{t_{k - 1}}}^ * \circ {{\boldsymbol{\hat f}}_{{t_{k - 1}}}} \circ {{\boldsymbol{\hat q}}_{{d_k}{t_{k - 1}}}} = {{\boldsymbol{\hat M}}_{{d_k}}}\frac{{{\partial ^{{d_k}}}}}{{\partial t}}{{\boldsymbol{\hat \omega }}_{{d_k}}} + {{\boldsymbol{\hat \omega }}_{{d_k}}} \times {{\boldsymbol{\hat M}}_{{d_k}}}{{\boldsymbol{\hat \omega }}_{{d_k}}} $$ (18) 假设铰链没有质量, 则铰链k−1之间力传递为
$$ {{\boldsymbol{\hat f}}_{{t_{k - 1}}}} = {\boldsymbol{\hat q}}_{{t_{k - 1}}{d_{k - 1}}}^ * \circ {{\boldsymbol{\hat f}}_{{d_{k - 1}}}} \circ {{\boldsymbol{\hat q}}_{{t_{k - 1}}{d_{k - 1}}}} $$ (19) $$ \begin{split} & {{{\boldsymbol{\hat f}}}_{{d_k}}} = {\boldsymbol{\hat q}}_{{d_k}{d_{k - 1}}}^ * \circ {{{\boldsymbol{\hat f}}}_{{d_{k - 1}}}} \circ {{{\boldsymbol{\hat q}}}_{{d_k}{d_{k - 1}}}} + {{{\boldsymbol{\hat M}}}_{{d_k}}}\frac{{{\partial ^{{d_k}}}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_{{d_k}}}+ \\ &\qquad {{{\boldsymbol{\hat \omega }}}_{{d_k}}} \times {{{\boldsymbol{\hat M}}}_{{d_k}}}{{{\boldsymbol{\hat \omega }}}_{{d_k}}} \end{split} $$ (20) 同样地, 记第k刚体在内接铰点${d_k}$处所受的力为${{\boldsymbol{\hat f}}_k}$, 内接铰点处的对偶质量阵记为${{\boldsymbol{\hat M}}_k}$. 于是有
$$ {{\boldsymbol{\hat f}}_k} = {\boldsymbol{\hat q}}_{kk - 1}^ * \circ {{\boldsymbol{\hat f}}_{k - 1}} \circ {{\boldsymbol{\hat q}}_{kk - 1}} + {{\boldsymbol{\hat M}}_k}\frac{{{\partial ^k}}}{{\partial t}}{{\boldsymbol{\hat \omega }}_k} + {{\boldsymbol{\hat \omega }}_k} \times {{\boldsymbol{\hat M}}_k}{{\boldsymbol{\hat \omega }}_k} $$ (21) 因为铰k的映射矩阵为${{\boldsymbol{\hat H}}_k}$ , 所以铰k处主动力
$$ {{\boldsymbol{\hat \tau }}_k} = {{\boldsymbol{\hat H}}_k}{{\boldsymbol{\hat f}}_k} $$ (22) 逆向动力学中, 铰链的广义加速度$ {\ddot \theta _k} $已知, 可以通过以下由端到基的递推和基到端的递推计算广义力和铰链力
$$\left. \begin{split} & {{{{\boldsymbol{\hat \omega }}}_{n + 1}} = {0}} \\ &\quad {{\mathrm{for}}}\quad {k = n: - 1:1} \\ &\quad {{{{\boldsymbol{\hat \omega }}}_k} = {\boldsymbol{\hat q}}_{kk + 1}^* \circ {{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ {{{\boldsymbol{\hat q}}}_{kk + 1}} + {{{\boldsymbol{\hat H}}}_k}{{\dot \theta }_k}} \\ &\quad {\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_k} = {{{\boldsymbol{\hat \alpha }}}_k} + {\boldsymbol{\hat q}}_{kk + 1}^* \circ \frac{{{\partial ^{k + 1}}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_{k + 1}} \circ {{{\boldsymbol{\hat q}}}_{kk + 1}} + {{{\boldsymbol{\hat H}}}_k}{{\ddot \theta }_k}} \\ &\quad {{\mathrm{end}}} \\ & {{{{\boldsymbol{\hat f}}}_0} = {0}} \\ &\quad {{\mathrm{for}}}\quad {k = 1:1:n} \\ &\quad {{{{\boldsymbol{\hat f}}}_k} = {\boldsymbol{\hat q}}_{kk - 1}^* \circ {{{\boldsymbol{\hat f}}}_{k - 1}} \circ {{{\boldsymbol{\hat q}}}_{kk - 1}} + {{{\boldsymbol{\hat M}}}_k}\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_k} + {{{\boldsymbol{\hat \omega }}}_k} \times {{{\boldsymbol{\hat M}}}_k}{{{\boldsymbol{\hat \omega }}}_k}} \\ &\quad {{{{\boldsymbol{\hat \tau }}}_k} = {{{\boldsymbol{\hat H}}}_k}{{{\boldsymbol{\hat f}}}_k}} \\ &\quad {{\mathrm{end}}} \end{split} \right\}$$ (23) 2.3 正向动力学递推
正向动力学中, 铰链上施加的对偶力$ {{\boldsymbol{\hat \varGamma }}_k} $已知, 需要求解铰链对应的广义加速度$ {\ddot \theta _k} $. 由于铰链上的约束反力未知, 即对偶力$ {{\boldsymbol{\hat f}}_k} $ 未知, 这时无法通过递推进行直接求解. 可以通过一定变换, 把递推方程转换为矩阵形式进行求解.
根据对偶四元数的运算规则, 任意的对偶四元数${{\boldsymbol{\hat \sigma }}_{k + 1}}$, 式${\boldsymbol{\hat q}}_{kk + 1}^ * \circ {{\boldsymbol{\hat \sigma }}_{k + 1}} \circ {{\boldsymbol{\hat q}}_{kk + 1}}$对于的矩阵形式可以写作
$$ \begin{split} &{\boldsymbol{\hat q}}_{kk + 1}^ * \circ {{{\boldsymbol{\hat \sigma }}}_{k + 1}} \circ {{{\boldsymbol{\hat q}}}_{kk + 1}}{\text{ = }}{\left[ {{\boldsymbol{\hat q}}_{kk + 1}^ * } \right]^ + }{\left[ {{{{\boldsymbol{\hat q}}}_{kk + 1}}} \right]^ - }{{{\boldsymbol{\hat \sigma }}}_{k + 1}} {\overset{\Delta}{\mathop{=}}}\\ &\qquad {{\boldsymbol{Q}}_{kk + 1}}{{{\boldsymbol{\hat \sigma }}}_{k + 1}}\end{split} $$ (24) 于是相应递推公式可以写作
$$ \left.\begin{aligned} & {{{\boldsymbol{\hat \omega }}}_k} = {{\boldsymbol{Q}}_{kk + 1}}{{{\boldsymbol{\hat \omega }}}_{k + 1}} + {{{\boldsymbol{\hat H}}}_k}{{\dot \theta }_k} \\ & \frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_k} = {{\boldsymbol{Q}}_{kk + 1}}\frac{{{\partial ^{k + 1}}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_{k + 1}} + {{{\boldsymbol{\hat \alpha }}}_k} + {{{\boldsymbol{\hat H}}}_k}{{\ddot \theta }_k} \\ & {{{\boldsymbol{\hat f}}}_k} = {{\boldsymbol{Q}}_{kk - 1}}{{{\boldsymbol{\hat f}}}_{k - 1}} + {{{\boldsymbol{\hat M}}}_k}\frac{{{\partial ^k}}}{{\partial t}}{{{\boldsymbol{\hat \omega }}}_k} + {{{\boldsymbol{\hat \omega }}}_k} \times {{{\boldsymbol{\hat M}}}_k}{{{\boldsymbol{\hat \omega }}}_k} \\ & {{{\boldsymbol{\hat \tau }}}_k} = {{{\boldsymbol{\hat H}}}_k}{{{\boldsymbol{\hat f}}}_k} \end{aligned} \right\} $$ (25) 定义矩阵
$$ \left. \begin{aligned} & {{\hat {\boldsymbol{\varOmega}} }} = {\text{col}}\left\{ {{{{{\hat {\boldsymbol{\omega}} }}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n}} \\ & {{\hat {\boldsymbol{\varGamma}} }} = {\text{col}}\left\{ {{{{{\hat {\boldsymbol{\tau}} }}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n}} \\ & {{\hat { {\boldsymbol{A}}}}} = {\text{col}}\left\{ {{{{{\hat {\boldsymbol{\alpha}} }}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n}} \\ & {{\hat {\boldsymbol{F}}}} = {\text{col}}\left\{ {{{{{\hat {\boldsymbol{f}}}}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n}} \\ & {{\hat {\boldsymbol{M}}}} = {\text{diag}}\left\{ {{{{{\hat {\boldsymbol{M}}}}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n \times 8n}} \\ & {{\hat {\boldsymbol{H}}}} = {\text{diag}}\left\{ {{{{{\hat {\boldsymbol{H}}}}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n \times N}} \\ & {\boldsymbol{\varTheta}} = {\text{col}}\left\{ {{\theta _k}} \right\}_{k = 1}^n \in {\mathbb{R}^N} \\ & {{\hat { {\boldsymbol{E}}}}} = {\text{col}}\left\{ {\frac{{{\partial ^k}}}{{\partial t}}{{{{\hat {\boldsymbol{\omega}} }}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n}} \\ & {{\hat { {\boldsymbol{B}}}}} = {\text{col}}\left\{ {{{{{\hat {\boldsymbol{\omega}} }}}_k} \times {{{{\hat {\boldsymbol{M}}}}}_k}{{{{\hat {\boldsymbol{\omega}} }}}_k}} \right\}_{k = 1}^n \in {\mathbb{R}^{8n}} \end{aligned} \right\} $$ (26) $$ {\boldsymbol{\hat Q}} = \left[ {\begin{array}{*{20}{l}} {{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}} \\ {{{\boldsymbol{Q}}_{{{12}}}}}&{{\boldsymbol{0}}}&{\boldsymbol{L}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}} \\ {{\boldsymbol{0}}}&{{{\boldsymbol{Q}}_{{{23}}}}}&{\boldsymbol{L}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}} \\ {\boldsymbol{M}}&{\boldsymbol{M}}&{{\boldsymbol{O}}}&{\boldsymbol{M}}&{\boldsymbol{M}} \\ {{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{\boldsymbol{L}}&{{{\boldsymbol{Q}}_{{{n - 1n}}}}}&{{\boldsymbol{0}}} \end{array}} \right] \in {\mathbb{R}^{8n \times 8n}} $$ (27) 其中n为刚体个数, N为铰链自由度之和.
于是递推公式可以写作
$$ \left. \begin{split} & {{\hat {\boldsymbol{\varOmega}} }} = {{{{\hat {\boldsymbol{Q}}}}}^{\text{T}}}{{\hat{\boldsymbol{ \varOmega }}}} + {{\hat {\boldsymbol{H}}}}\dot \varTheta \\ & {{\hat {{\boldsymbol{E}}}}} = {{{{\hat {\boldsymbol{Q}}}}}^{\text{T}}}{{\hat {{\boldsymbol{E}}}}} + {{\hat {{\boldsymbol{A}}}}} + {{\hat {\boldsymbol{H}}}}\ddot \varTheta \\ & {{\hat {\boldsymbol{F}}}} = {{\hat {\boldsymbol{Q}}\hat {\boldsymbol{F}}}} + {{\hat {\boldsymbol{M}}\hat { {\boldsymbol{E}}}}} + {{\hat { {\boldsymbol{B}}}}} \\ & {{\hat{\boldsymbol{ \varGamma}} }} = {{\hat {\boldsymbol{H}}\hat {\boldsymbol{F}}}} \end{split} \right\} $$ (28) 从矩阵Q结构可知为幂零矩阵, 其幂指数为n, 根据幂零矩阵性质可知$ \left( {{\boldsymbol{I}} + {\boldsymbol{\hat Q}}} \right) $是可逆的, 可以证明
$$ {\boldsymbol{\varPhi }} {\overset{\Delta}{\mathop{=}}} {\left( {{\boldsymbol{I}} + {\boldsymbol{\hat Q}}} \right)^{ - 1}} = \left[ {\begin{array}{*{20}{l}} {{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}}&{{\boldsymbol{0}}} \\ {{{\boldsymbol{Q}}_{12}}}&{{\boldsymbol{0}}}& \cdots &{{\boldsymbol{0}}}&{{\boldsymbol{0}}} \\ {{{\boldsymbol{Q}}_{13}}}&{{{\boldsymbol{Q}}_{23}}}& \cdots &{{\boldsymbol{0}}}&{{\boldsymbol{0}}} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ {{{\boldsymbol{Q}}_{1n}}}&{{{\boldsymbol{Q}}_{2n}}}& \cdots &{{{\boldsymbol{Q}}_{n - 1n}}}&{{\boldsymbol{0}}} \end{array}} \right] $$ (29) 于是递推公式(28)可以写作
$$ \left.\begin{split} & {{\hat {\boldsymbol{\varOmega}} }} = {{\boldsymbol{\varPhi }}^{\text{T}}}{{\hat {\boldsymbol{H}}}}\dot {\boldsymbol{\varTheta}} \\ & {{\hat { {\boldsymbol{E}}}}} = {{\boldsymbol{\varPhi }}^{\text{T}}}\left( {{{\hat { {\boldsymbol{A}}}}} + {{\hat {\boldsymbol{H}}}}\ddot {\boldsymbol{\varTheta}} } \right) \\ & {{\hat {\boldsymbol{F}}}} = {\boldsymbol{\varPhi }}\left( {{{\hat {\boldsymbol{M}}\hat { {\boldsymbol{E}}}}} + {{\hat { {\boldsymbol{B}}}}}} \right) \\ & {{\hat {\boldsymbol{\varGamma }}}} = {{\hat {\boldsymbol{H}}\hat {\boldsymbol{F}}}} \end{split} \right\} $$ (30) 上式可以推导为以下形式的动力学方程
$$ {{\boldsymbol{M}}}(\theta )\ddot \theta + {{\boldsymbol{C}}}(\theta ,\dot \theta ) = {\boldsymbol{\hat \varGamma }} $$ (31) 其中${{\boldsymbol{M}}}(\theta ) = \left( {{{\hat {\boldsymbol{H}}{\boldsymbol{\varPhi}} \hat {\boldsymbol{M}}}}{{\boldsymbol{\varPhi }}^{\text{T}}}{{\hat {\boldsymbol{H}}}}} \right)$, ${{\boldsymbol{C}}}(\theta ,\dot \theta ) = {{\hat {\boldsymbol{H}}{\boldsymbol{\varPhi}} \hat {\boldsymbol{M}}}}{{\boldsymbol{\varPhi }}^{\text{T}}}{{\hat {{\boldsymbol{A}}}}} + {\boldsymbol{\hat H\varPhi \hat { {\boldsymbol{B}}}}}$,于是可以得到正向动力学求解
$$ \ddot{\boldsymbol{ \varTheta}} = {\left( {{{\hat {\boldsymbol{H}}{\boldsymbol{\varPhi}} \hat{\boldsymbol{ M}}}}{{\boldsymbol{\varPhi }}^{\text{T}}}{{\hat{\boldsymbol{ H}}}}} \right)^{ - 1}}\left( {{{\hat {\boldsymbol{\varGamma}} }} - {{\hat {\boldsymbol{H}}{\boldsymbol{\varPhi}} \hat {\boldsymbol{M}}}}{{\boldsymbol{\varPhi }}^{\text{T}}}{{\hat {{\boldsymbol{A}}}}} - {{\hat {\boldsymbol{H}}{\boldsymbol{\varPhi}} \hat { {\boldsymbol{B}}}}}} \right) $$ (32) 3. 执行机构建模
3.1 推进器
推进器在卫星本体上产生的力和力矩为
$$\qquad\qquad {\boldsymbol{f}}_{}^B = \sum\limits_{i = 1}^S {{\boldsymbol{f}}_i^B} $$ (33) $$\qquad\qquad {\boldsymbol{t}}_{}^B = \sum\limits_{i = 1}^S {{{\boldsymbol{R}}_i}} \times {\boldsymbol{f}}_i^B $$ (34) 式中S为推进器数量, ${{\boldsymbol{R}}_i}$为推进器相对于本体系原点的安装位置. 根据文献[20-21], 推进器产生的动力学效应可以表示为
$$ {\boldsymbol{\hat f}} = \left. {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{f}}_{}^B} \\ {{\boldsymbol{\tau }}_{}^B} \end{array}} \right.} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{H}}_f}} \\ {{{\boldsymbol{H}}_t}} \end{array}} \right]{\boldsymbol{\sigma }} = {\boldsymbol{H\sigma }} $$ (35) 其中${{\boldsymbol{H}}_f},{{\boldsymbol{H}}_t} \in {\mathbb{R}^{3 \times S}}$分别为推进器与安装体上的力和力矩映射矩阵, 它与推进器相对于载体质心的位置有关. 矢量${\boldsymbol{\sigma}} = {({\sigma _1}, \cdots ,{\sigma _i}, \cdots ,{\sigma _S})^{\mathrm{T}}} \in {\mathbb{R}^S}$为各个推进器沿推力器安装方向施加的推力.
3.2 动量交换装置
控制力矩陀螺(VSCMGs)是一种用于航天器姿态控制的惯性执行机构, 它主要通过交换动量来产生控制力矩, 从而实现对航天器姿态的控制或稳定[22-23]. 这种装置具有结构简单和高输出扭矩等特点. 控制力矩陀螺应用时需要考虑奇异回避问题, 并对此设计回避奇异点的控制率[24-26]. 本文重点介绍动力学建模过程, 控制律设计采用奇异鲁棒逆算法[27-28]. 考虑执行机构包含有K个变速控制力矩陀螺, 可以将每个控制力矩陀螺等效为独立的运动附件, 每个VSCMG产生的力和力矩用对偶四元数的形式可以表示为[18]
$$ \begin{split} & {{{{\hat {\boldsymbol{F}}}}}_{PB}} = {{{{\hat {\boldsymbol{M}}}}}_{PB}}\frac{{{\partial ^B}}}{{\partial t}}{{{{\hat {\boldsymbol{\omega}} }}}_B} + \\ &\qquad {{{{\hat {\boldsymbol{D}}}}}_{BP}}\left[ {{{{{\hat {\boldsymbol{\omega }}}}}_p} \times {{{{\hat {\boldsymbol{M}}}}}_P}{{{{\hat {\boldsymbol{\omega}} }}}_p} + {{{{\hat {\boldsymbol{M}}}}}_P}\left( {{{{{\hat {\boldsymbol{\omega}} }}}_p} \times {{{{\hat {\boldsymbol{\omega}} }}}_{PB}} + {{{{\dot {\hat {\boldsymbol{\omega}}} }}}_{PB}}} \right)} \right] \end{split} $$ (36) 其中, $ {{\boldsymbol{\hat M}}_{PB}} $为控制力矩陀螺相对于中心刚体本体坐标系的对偶质量阵, $ {{{\hat {\boldsymbol{\omega}} }}_p} $为控制力矩陀螺相对惯性系的对偶速度, $ {{{\hat {\boldsymbol{\omega}} }}_{PB}} $为控制力矩陀螺相对中心刚体的相对对偶速度, $ {{{\hat {\boldsymbol{D}}}}_{BP}} $表示了控制力矩陀螺相对中心刚体的安装位置与姿态, $ {{{\hat {\boldsymbol{F}}}}_{PB}} $表示VSCMG在中心刚体体坐标系上作用的力和力矩.
中心刚体与控制力矩陀螺组成的系统动力学方程为
$$ \begin{split} & {{{{\hat {\boldsymbol{F}}}}}_B} = \left( {{{{{\hat {\boldsymbol{M}}}}}_B} + \sum\limits_{i = 1}^K {{{{{\hat {\boldsymbol{M}}}}}_{P{B_i}}}} } \right)\frac{{{\partial ^B}}}{{\partial t}}{{{{\hat {\boldsymbol{\omega}} }}}_B} + {{{{\hat {\boldsymbol{\omega}} }}}_B} \times {{{{\hat {\boldsymbol{M}}}}}_B}{{{{\hat {\boldsymbol{\omega}} }}}_B} + \\ &\quad \sum\limits_{i = 1}^K {\left\{ {{{{{\hat {\boldsymbol{D}}}}}_{B{P_i}}}\left[ {{{{{\hat {\boldsymbol{\omega}} }}}_{{p_i}}} \times \left( {{{{{\hat {\boldsymbol{M}}}}}_{{p_i}}}{{{{\hat {\boldsymbol{\omega}} }}}_{{p_i}}}} \right) + {{{{\hat{\boldsymbol{ M}}}}}_{{P_i}}}\left( {{{{{\hat {\boldsymbol{\omega}} }}}_{{p_i}}} \times {{{{\hat {\boldsymbol{\omega}} }}}_{P{B_i}}} + {{{{\dot {\hat{\boldsymbol{ \omega}}} }}}_{P{B_i}}}} \right)} \right]} \right\}} \end{split} $$ (37) 由式(37)可知, 每个控制力矩陀螺产生的控制力矩为
$$ \begin{split} &{{\boldsymbol{\tau }}_i} = {\text{dual}}\left\{ {{{{\hat {\boldsymbol{D}}}}}_{B{P_i}}}\left[ {{{{\hat {\boldsymbol{\omega}} }}}_{{p_i}}} \times \left( {{{{{\hat {\boldsymbol{M}}}}}_{{p_i}}}{{{{\hat {\boldsymbol{\omega}} }}}_{{p_i}}}} \right) +\right.\right.\\ &\qquad \left.\left.{{{{\hat {\boldsymbol{M}}}}}_{{P_i}}}\left( {{{{{\hat {\boldsymbol{\omega}} }}}_{{p_i}}} \times {{{{\hat {\boldsymbol{\omega}} }}}_{P{B_i}}} + {{{{\dot \hat {\boldsymbol{\omega}} }}}_{P{B_i}}}} \right) \right] \right\}\end{split} $$ (38) 其中dual(·)表示取对偶四元数的对偶部, 设VSCMG的框架角速度和转子角速度分别为$\dot \gamma 和\varOmega $, 则
$$ {{\boldsymbol{\omega }}_{P{B_i}}} = {{\boldsymbol{A}}_i}\left[ {\begin{array}{*{20}{c}} {{{\dot \gamma }_i}} \\ {{\varOmega _i}} \end{array}} \right] $$ (39) 式中$ {{\boldsymbol{A}}_i} $是角速率与角速度变换矩阵.
将式(39)代入式(38)可以化简得
$$\left. \begin{split} & {{\boldsymbol{\tau }}_i} = {{\boldsymbol{B}}_i}{{\ddot \gamma }_i} + {{\boldsymbol{C}}_i}{{\dot \gamma }_i} + {{\boldsymbol{E}}_i}{{\dot \varOmega }_i} \\ & {{\boldsymbol{B}}_i} = {{\boldsymbol{A}}_g}{{\boldsymbol{I}}_{cg}} \\ & {{\boldsymbol{C}}_i} = {\boldsymbol{w}} + {{\boldsymbol{A}}_t}{{\boldsymbol{I}}_{ws}}{[\varOmega ]^d} + {\boldsymbol{\tilde \omega }}{{\boldsymbol{A}}_g}{{\boldsymbol{I}}_{cg}} \\ & {{\boldsymbol{E}}_i} = {{\boldsymbol{A}}_s}{{\boldsymbol{I}}_{ws}} \end{split} \right\}$$ (40) 式中, 下标s, t和g分别代表任一VSCMG自旋轴、输出力矩相反方向轴和框架轴. ${{\boldsymbol{A}}_ * } = \left[ {{{\boldsymbol{e}}_{ * 1}}}\;\;{{{\boldsymbol{e}}_{ * 2}}}\;\; \cdots \;\;{{{\boldsymbol{e}}_{ * k}}} \right]$, ${{\boldsymbol{e}}_{ * k}}$为第k个控制力矩陀螺的s或t或g轴在中心刚体本体系中的方向余弦阵. $ {{\boldsymbol{I}}_{{{ws}}}} = {\mathrm{diag}}\left[ {{{\boldsymbol{I}}_{{{ws1}}}}}\;\;{{{\boldsymbol{I}}_{{{ws2}}}}}\;\; \cdots \;\;{{{\boldsymbol{I}}_{{{wsk}}}}} \right] $为控制力矩陀螺转子相对于s轴的转动惯量, $ {{\boldsymbol{I}}_{{{cg}}}} = {\mathrm{diag}}\left[ {{{\boldsymbol{I}}_{{{cg1}}}}}\;\;{{{\boldsymbol{I}}_{{{cg2}}}}}\;\; \cdots \;\;{{{\boldsymbol{I}}_{{{cgk}}}}} \right] $为控制力矩陀螺相对于g轴的转动惯量. 根据文献[25], $ {\boldsymbol{B}}{\ddot \gamma _i} $的影响较小可以忽略, 于是K个VSCMG产生的合力矩可以写作
$$ {\boldsymbol{\tau }} = [\begin{array}{*{20}{c}} {\boldsymbol{C}}&{\boldsymbol{E}} \end{array}]\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\dot \gamma }}} \\ {{\boldsymbol{\dot \varOmega }}} \end{array}} \right] {\overset{\Delta}{\mathop{=}}} \;{\boldsymbol{G}}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\dot \gamma }}} \\ {{\boldsymbol{\dot \varOmega }}} \end{array}} \right] $$ (41) 其中, $ {\boldsymbol{C}} = \left[ {{{\boldsymbol{C}}_1}}\;\;{{{\boldsymbol{C}}_2}}\;\; \cdots \;\;{{{\boldsymbol{C}}_K}} \right],{\boldsymbol{E}} = \left[ {{{\boldsymbol{E}}_1}}\;\;{{{\boldsymbol{E}}_2}}\;\; \cdots \;\;{{{\boldsymbol{E}}_K}} \right] $, $ {\boldsymbol{\dot \gamma }} = \left[ {{{\dot \gamma }_1}}\;\;{{{\dot \gamma }_2}}\;\; \cdots {{{\dot \gamma }_K}} \right]^{\text{T}}, {\boldsymbol{\dot \varOmega }} = {\left[ {{{\dot \varOmega }_1}}\;\; {{{\dot \varOmega }_2}}\;\; \cdots \;\;{{{\dot \varOmega }_K}} \right]^{\text{T}}} $.
式(41)的加权范数解为
$$ \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\dot \gamma }}} \\ {{\boldsymbol{\dot \varOmega }}} \end{array}} \right] = {{\boldsymbol{G}}^{\text{T}}}{[{\boldsymbol{G}}{{\boldsymbol{G}}^{\text{T}}} + \lambda {\boldsymbol{R}}]^{ - 1}}{\boldsymbol{\tau }} $$ (42) 其中, $\lambda = {\lambda _0}\exp [ - \mu \det ({\boldsymbol{G}}{{\boldsymbol{G}}^{\text{T}}})]$, $ {\boldsymbol{R}} = \left[\begin{array}{ccc}1& {\varepsilon }_{3}& {\varepsilon}_{2}\\ {\varepsilon}_{3}& 1& {\varepsilon}_{1}\\ {\varepsilon}_{2}& {\varepsilon}_{1}& 1\end{array}\right] $, $ {\varepsilon}_{i}和 {\lambda }_{0} $为正实数.
4. 控制实现
在本方法中, 解耦了机械手的控制和卫星基座的控制. 机械臂将使用对偶PD控制方法进行控制, 而卫星基座将使用推进器和VSCMG的组合进行控制.
4.1 机械臂控制
安装在卫星基座上的有N个关节的机械臂动力学方程如式(27)所示. 机械臂的控制目标是控制每个关节跟踪任务空间的参考轨迹$ {{\boldsymbol{\hat q}}_{di}}(t) \in \boldsymbol{H} $和参考速度${{\boldsymbol{\hat \omega }}_{di}}(t)$, 其中$ {{\boldsymbol{\hat q}}_{di}}(t) $是每个臂的参考位姿, ${{\boldsymbol{\hat \omega }}_{di}}(t)$为每个臂的参考广义速度.
若臂i的位姿用对偶四元数描述${{\boldsymbol{\hat q}}_i}$, 速度用描述${{\boldsymbol{\hat \omega }}_i}$, 定义位姿误差为
$$ {\boldsymbol{\hat e}} = {\boldsymbol{\hat q}}_{di}^* \circ {{\boldsymbol{\hat q}}_i} $$ (43) 根据对偶四元数的定义展开可以得到
$$ {\boldsymbol{\hat e}} = {{\boldsymbol{e}}_q} + \frac{{{\varepsilon }}}{{\text{2}}}{{\boldsymbol{e}}_q} \circ {{\boldsymbol{e}}_r} $$ (44) 可以看出位姿误差对偶四元数的实部为姿态误差$ {{\boldsymbol{e}}_q} = {\boldsymbol{q}}_{di}^* \circ {{\boldsymbol{q}}_i} $, 对偶部包含了位置误差$ {{\boldsymbol{e}}_r} = {{\boldsymbol{r}}_{di}} - {{\boldsymbol{r}}_i} $.
定义速度误差为
$$ {{\boldsymbol{\hat \omega }}_e} = {{\boldsymbol{\hat \omega }}_{di}} - {{\boldsymbol{\hat \omega }}_i} $$ (45) 根据文献[13],可以通过广义PD控制对每个关节进行控制, 定义控制律为
$$ {{\boldsymbol{\hat u}}_i} = - {{\boldsymbol{k}}_{pi}}\lg {{\boldsymbol{\hat e}}_i} - {{\boldsymbol{k}}_{di}}{{\boldsymbol{\hat \omega }}_{ei}} $$ (46) 其中增益系数$ {{\boldsymbol{k}}_{pi}}和{{\boldsymbol{k}}_{di}} $都是正定的.
4.2 基座控制
在机械臂运动时, 使用VSCMG和推进器混合控制基座的位姿. 推进器用于控制基座的位置, VSCMG控制基座姿态. 控制律设计为
$$ {\boldsymbol{\hat f}} = - {{\boldsymbol{k}}_p}\lg {{\boldsymbol{\hat e}}_i} - {{\boldsymbol{k}}_d}{{\boldsymbol{\hat \omega }}_{ei}} - {{\boldsymbol{\hat f}}_0} $$ (47) 其中前馈项$ {{\boldsymbol{\hat f}}_0} $为机械臂运动时产生的力和力矩, 增益系数$ {{\boldsymbol{k}}_p}和{{\boldsymbol{k}}_d} $都是正定的. 执行机构产生的对偶主动力为
$$ {\boldsymbol{\hat f}} = {\boldsymbol{H\sigma }} + \left. {\left[ {\begin{array}{*{20}{c}} 0 \\ {{{\boldsymbol{\tau }}^B}} \end{array}} \right.} \right] = \left. {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{H}}_f}{\boldsymbol{\sigma }}} \\ {{{\boldsymbol{H}}_t}{\boldsymbol{\sigma }} + {{\boldsymbol{\tau }}^B}} \end{array}} \right.} \right] $$ (48) 由式(47)与式(48)可以先求推力器的解, 该问题可以转化为如下的最优求解问题[29]
$$\left.\begin{split} & \min \sum\limits_{i = 1}^n {{\sigma _i}} \\ & {\text{s}}{\text{.t}}{\text{. }}\left\{ \begin{aligned} & {{{\boldsymbol{H}}_f}{\boldsymbol{\sigma }} = {\boldsymbol{f}}} \\ & {{\sigma _{\min }} \leqslant {\sigma _i} \leqslant {\sigma _{\max }},\quad i = 1,2, \cdots ,S} \end{aligned} \right. \end{split} \right\}$$ (49) 求出推力器的控制率后, 代入式(48)的对偶部分结合式(42)即可求得每个VSCMG的控制规律.
5. 算例
本节对动力学建模和控制律设计进行仿真. 取建模对象为KUKA iiwa 7自由度机械臂[30]固连在一个漂浮基座上, 基座四周装有4个VSCMG, 基座的4个顶点安装3组推力器, 执行机构具体安装示意如图5和图6所示.
基座、机械臂以及控制力矩陀螺的惯量参数可以参考表3 ~ 表5.
表 3 基座惯量参数Table 3. Base inertia parametersParameter Value mass 27 kg inertia ${I_{xx}} = {I_{yy}} = {I_{zz}} = 0.405{\text{ kg}} \cdot {{\text{m}}^2}$ center of mass ${c_x} = {c_y} = {c_z} = 0$ dimension ${s_x} = {s_y} = {s_z} = 0.3{\text{ m}}$ 表 4 控制力矩陀螺惯量参数Table 4. VSCMG inertia parametersParameter Value mass 24.66 kg inertia $\begin{gathered} {I_{xx}} = 0.123\;3{\text{ kg}} \cdot {{\text{m}}^2} \\ {I_{yy}} = {I_{zz}} = 0.082\;21{\text{ kg}} \cdot {{\text{m}}^2} \\ \end{gathered} $ center of mass ${c_x} = {c_y} = {c_z} = 0$ 表 5 机械臂惯量参数Table 5. Arm inertia parametersMass/kg Inertia/(${\text{kg}} \cdot {{\text{m}}^2}$)
[${I_{xx}},{I_{yy}},{I_{zz}},{I_{yz}},{I_{xz}},{I_{xy}}$]Center of mass/m
[$x,y,z$]D-H parameters
$[{a_i},{d_i},{\alpha _i},{\theta _{0,i}}]$link 1 4 [0.1612, 0.1476, 0.0236, 0.0144, 0, 0] [0, −0.0300, 0.1200] [0 m,0.36 m,$ - {90^\circ }$,${0^\circ }$] link 2 4 [0.071, 0.0251, 0.0579, −0.0099, −5.04 × 10−5, −7.08 × 10−5] [3.0 × 10−4, 0.0590, 0.0420] [0 m,0 m,${90^\circ }$,${0^\circ }$] link 3 3 [0.1334, 0.1257, 0.0127, −0.0117, 0, 0] [0, 0.0300, 0.1300] [0 m,0.42 m,$ - {90^\circ }$,${0^\circ }$] link 4 2.7 [0.0452, 0.0131, 0.0411, −0.0062, 0, 0] [0, 0.0670, 0.0340] [0 m,0 m,${90^\circ }$,${0^\circ }$] link 5 1.7 [0.0306, 0.0278, 0.0057, −0.0027, −1.292 × 10−5, −3.57 × 10−6] [1.0 × 10−4, 0.0210, 0.0760] [0 m,0.4 m,$ - {90^\circ }$,${0^\circ }$] link 6 1.8 [0.0050, 0.0036, 0.0047, −4.320 × 10−7, 0, 0] [0, 6.0 × 10−4, 4.0 × 10−4] [0 m,0 m,${90^\circ }$,${0^\circ }$] link 7 0.3 [0.0011, 0.0011, 1.0 × 10−3, 0, 0, 0] [0, 0, 0.0200] [0 m,0.126 m,${0^\circ }$,${0^\circ }$] 5.1 动力学验证
本节通过算例验证动力学建模的正确性. 仿真模型仅包含基座和机械臂, 在基座上施加如下力和力矩
$$\left. \begin{split} & {{{\boldsymbol{F}}_b}(t) = \sin (t) \cdot {{[1,2,3]}^{\text{T}}}} \\ & {{{\boldsymbol{T}}_b}(t) = \sin (t) \cdot {{[1,2,3]}^{\text{T}}}} \end{split} \right\}$$ (50) 基座和机械臂的初始状态均为0. 使用第3节方法与Simscape计算作对比, 两种方法都使用ODE45进行计算, 仿真参数设置为: 最大步长0.001, 相对误差和绝对误差都设为1.0 × 10−8. 通过对比两种方法计算得到的机械臂末端位置进行误差对比, 结果如图7和图8所示.
可以看到两种方法计算结果十分接近. 同时本算例还进行了不同步长下的误差结果对比, 从图9可以看出当进一步缩小积分步长时, 两种方法误差可以进一步缩小. 由此可以验证对偶四元数方法的正确性.
为了进一步比较对偶四元数方法与传统方法的效率, 在本算例中还对比了对偶四元数方法与指数积方法所需的正向运动学递推基本运算次数. 指数积方法对正逆运动学问题的求解过程已有较多文献[31-32]进行过详细介绍和分析.
依据式(22), 每一次正向运动学递推所需执行的计算即求变量$ {{\boldsymbol{\hat \omega }}_k} = {\boldsymbol{\hat q}}_{kk + 1}^* \circ {{\boldsymbol{\hat \omega }}_{k + 1}} \circ {{\boldsymbol{\hat q}}_{kk + 1}} + {{\boldsymbol{\hat H}}_k}{\dot \theta _k} $的过程. 根据式(2)描述的对偶四元数乘法的矩阵表示法, 计算$ {\boldsymbol{\hat q}}_{kk + 1}^* \circ {{\boldsymbol{\hat \omega }}_{k + 1}} \circ {{\boldsymbol{\hat q}}_{kk + 1}} $需要进行72次乘法操作和56次加法操作, 而本算例中涉及的运动副均为旋转类型, $ {{\boldsymbol{\hat H}}_k}{\dot \theta _k} $中仅有单一元素非零, 因此每次递推过程需执行72次乘法和57次加法. 相比之下, 采用指数积递推公式时, 每次求解需执行96次乘法和66次加法[31]. 这一结果说明, 使用对偶四元数进行运算所需的算术操作数量显著少于采用指数积公式所需的运算数量, 凸显了其在计算资源优化方面的潜在优势.
5.2 闭环控制仿真
本节通过算例验证整体闭环控制的正确性. 控制的目标为机械臂末端沿着给定路径跟踪预定义的轨迹, 同时将卫星基座保持在初始姿态和位置. 卫星基座初始时刻静止在原点, 初始姿态均为零. 机械臂末端轨迹为三维双纽线, 对应的参数方程为
$$ {\boldsymbol{r}}(t) = \left\{ \begin{split} & {x = r[1 + \cos (\zeta t + {\theta _m})] + {x_0}} \\ & {y = r\sin (\zeta t + {\theta _m}) + {y_0}} \\ & {z = 2r\sin (\zeta t/2 + {\theta _m}/2) + {z_0}} \end{split} \right. $$ (51) 本节中参数选取如下: $r = 0.1\;{\mathrm{m}},$ $ \zeta = 1, $ $ {\theta _m} = 0 ,$ ${x_0} = 0.3,$ ${y_0} = 0$和${z_0} = 0.7.$ 给定机械臂末端轨迹后, 可以通过逆向运动学求解每个关节的目标位姿. 然后根据第4节内容分别控制机械臂和基座位姿.
仿真结果如下: 从图10和图11结果可以看到, 在推力器和控制力矩陀螺的作用下, 姿态稳定在初始状态附近, 基座的位置保持精度在mm级. 图12 ~ 图14为执行机构输出值, 在第9 s左右机械臂的运动产生了较大扰动, 控制系统能较快地响应使基座位姿仍然保持较高的控制精度.
6. 结 论
本文研究了使用对偶四元数进行多刚体系统建模和控制的问题. 使用对偶四元数可以将位置与姿态在统一的建模框架中表达, 有利于后续的控制系统设计. 通过铰链间物理关系递推, 完成了运动学和动力学递推方程的推导. 为了便于控制系统仿真计算, 完成了正向动力学方程的矩阵形式推导, 可以直接通过关节上的对偶主动力求解广义加速度. 对于实际航天应用中的漂浮基座和机械臂组合体的建模与控制, 本文分别讨论了基座和机械臂的控制率设计方式, 并使用控制力矩陀螺和推力器完成了7自由度机械臂的位姿跟踪和基座稳定控制.
未来的研究计划包括两方面. 一方面是继续完善此仿真框架的内容, 包括将刚体动力学模型拓展到柔性附件及柔性体建模等; 另一方面, 可基于此仿真框架, 研究使用对偶四元数完成轨迹规划的内容, 例如在最小干扰约束下, 指定机械臂末端轨迹时的机械臂的任务空间位姿规划.
-
表 1 部分对偶四元数运算法则
Table 1 Dual quaternion operation rules
Operation Formula expression multiplication ${\hat \sigma _1} \circ {\hat \sigma _2} = {{\boldsymbol{p}}_1} \circ {{\boldsymbol{p}}_2} + {{\varepsilon }}\left( {{{\boldsymbol{p}}_1} \circ {{\boldsymbol{q}}_2} + {{\boldsymbol{q}}_1} \circ {{\boldsymbol{p}}_2}} \right)$ conjugate ${\hat \sigma ^*} = {{\boldsymbol{p}}^*} + {{\varepsilon }}{{\boldsymbol{q}}^*}$ cross product ${\hat \sigma _1} \times {\hat \sigma _2} = {{\boldsymbol{p}}_1} \times {{\boldsymbol{p}}_2} + {{\varepsilon }}\left( {{{\boldsymbol{p}}_1} \times {{\boldsymbol{q}}_2} + {{\boldsymbol{p}}_2} \times {{\boldsymbol{q}}_1}} \right)$ frame transformation ${\boldsymbol{\hat r}}_b^b = {\boldsymbol{\hat q}}_{ba}^{a*} \circ {\boldsymbol{\hat r}}_a^a \circ {\boldsymbol{\hat q}}_{ba}^a$ logarithm $ \log \left( {{\boldsymbol{p}} + {{\varepsilon }}{\boldsymbol{q}}} \right) = \log {\boldsymbol{p}} + {{\varepsilon }}\left( {{{\boldsymbol{p}}^{ - 1}} \circ {\boldsymbol{q}}} \right) $ 表 2 常见铰链的映射矩阵和广义速度
Table 2 Mapping matrix and generalized velocity for common hinges
Hinge type Map matrix Generalized velocity revolute $\left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 3}}} \\ 1 \\ {{{{\boldsymbol{0}}}_{1 \times 4}}} \end{array}} \right)$ $ \dot \theta $ prismatic $\left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 7}}} \\ 1 \end{array}} \right)$ $\dot d$ spherical $\left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 3}}} \\ {{{\boldsymbol{I}}_3}} \\ {{{{\boldsymbol{0}}}_{4 \times 3}}} \end{array}} \right)$ $\left( {\begin{array}{*{20}{c}} {{\omega _x}} \\ {{\omega _y}} \\ {{\omega _z}} \end{array}} \right)$ 6-DoF joint $ \left( {\begin{array}{*{20}{c}} {{{{\boldsymbol{0}}}_{1 \times 3}}}&{{{{\boldsymbol{0}}}_{1 \times 3}}} \\ {{{\boldsymbol{I}}_3}}&{{{{\boldsymbol{0}}}_{3 \times 3}}} \\ {{{{\boldsymbol{0}}}_{1 \times 3}}}&{{{{\boldsymbol{0}}}_{1 \times 3}}} \\ {{{{\boldsymbol{0}}}_{3 \times 3}}}&{{{\boldsymbol{I}}_3}} \end{array}} \right) $ $\left( {\begin{array}{*{20}{c}} {{\omega _x}} \\ {{\omega _y}} \\ {{\omega _z}} \\ {{v_x}} \\ {{v_y}} \\ {{v_z}} \end{array}} \right)$ 表 3 基座惯量参数
Table 3 Base inertia parameters
Parameter Value mass 27 kg inertia ${I_{xx}} = {I_{yy}} = {I_{zz}} = 0.405{\text{ kg}} \cdot {{\text{m}}^2}$ center of mass ${c_x} = {c_y} = {c_z} = 0$ dimension ${s_x} = {s_y} = {s_z} = 0.3{\text{ m}}$ 表 4 控制力矩陀螺惯量参数
Table 4 VSCMG inertia parameters
Parameter Value mass 24.66 kg inertia $\begin{gathered} {I_{xx}} = 0.123\;3{\text{ kg}} \cdot {{\text{m}}^2} \\ {I_{yy}} = {I_{zz}} = 0.082\;21{\text{ kg}} \cdot {{\text{m}}^2} \\ \end{gathered} $ center of mass ${c_x} = {c_y} = {c_z} = 0$ 表 5 机械臂惯量参数
Table 5 Arm inertia parameters
Mass/kg Inertia/(${\text{kg}} \cdot {{\text{m}}^2}$)
[${I_{xx}},{I_{yy}},{I_{zz}},{I_{yz}},{I_{xz}},{I_{xy}}$]Center of mass/m
[$x,y,z$]D-H parameters
$[{a_i},{d_i},{\alpha _i},{\theta _{0,i}}]$link 1 4 [0.1612, 0.1476, 0.0236, 0.0144, 0, 0] [0, −0.0300, 0.1200] [0 m,0.36 m,$ - {90^\circ }$,${0^\circ }$] link 2 4 [0.071, 0.0251, 0.0579, −0.0099, −5.04 × 10−5, −7.08 × 10−5] [3.0 × 10−4, 0.0590, 0.0420] [0 m,0 m,${90^\circ }$,${0^\circ }$] link 3 3 [0.1334, 0.1257, 0.0127, −0.0117, 0, 0] [0, 0.0300, 0.1300] [0 m,0.42 m,$ - {90^\circ }$,${0^\circ }$] link 4 2.7 [0.0452, 0.0131, 0.0411, −0.0062, 0, 0] [0, 0.0670, 0.0340] [0 m,0 m,${90^\circ }$,${0^\circ }$] link 5 1.7 [0.0306, 0.0278, 0.0057, −0.0027, −1.292 × 10−5, −3.57 × 10−6] [1.0 × 10−4, 0.0210, 0.0760] [0 m,0.4 m,$ - {90^\circ }$,${0^\circ }$] link 6 1.8 [0.0050, 0.0036, 0.0047, −4.320 × 10−7, 0, 0] [0, 6.0 × 10−4, 4.0 × 10−4] [0 m,0 m,${90^\circ }$,${0^\circ }$] link 7 0.3 [0.0011, 0.0011, 1.0 × 10−3, 0, 0, 0] [0, 0, 0.0200] [0 m,0.126 m,${0^\circ }$,${0^\circ }$] -
[1] 薛智慧, 刘金国. 空间机械臂操控技术研究综述. 机器人, 2022, 44(1): 107-128 (Xue Zhihui, Liu Jinguo. Review of space manipulator control technologies. Robot, 2022, 44(1): 107-128 (in Chinese) Xue Zhihui, Liu Jinguo. Review of space manipulator control technologies. Robot, 2022, 44(1): 107-128 (in Chinese)
[2] 王明明, 罗建军, 袁建平等. 空间在轨装配技术综述. 航空学报, 2021, 42(1): 523913 (Wang Mingming, Luo Jianjun, Yuan Jianping, et al. In-orbit assembly technology: Review. Acta Aeronautica et Astronautica Sinica, 2021, 42(1): 523913 (in Chinese) Wang Mingming, Luo Jianjun, Yuan Jianping, et al. In-orbit assembly technology: Review. Acta Aeronautica et Astronautica Sinica, 2021, 42(1): 523913 (in Chinese)
[3] Abad AF, Ma O, Pham K, et al. A review of space robotics technologies for on-orbit servicing. Progress in Aerospace Science, 2014, 68(8): 1-26
[4] 黄文虎, 曹登庆, 韩增尧. 航天器动力学与控制的研究进展与展望. 力学进展, 2012, 42(4): 367-394 (Huang Wenhu, Cao Dengqing, Han Zengyao. Advances and trends in dynamics and control of spacecrafts. Advances in Mechanics, 2012, 42(4): 367-394 (in Chinese) Huang Wenhu, Cao Dengqing, Han Zengyao. Advances and trends in dynamics and control of spacecrafts. Advances in Mechanics, 2012, 42(4): 367-394 (in Chinese)
[5] 丰保民. 自由漂浮空间机器人轨迹规划与轨迹跟踪问题研究. [博士论文]. 哈尔滨: 哈尔滨工业大学, 2008 (Feng Baomin. Study on path planning and trajectory tracking control of free-floating space robot. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2008 (in Chinese) Feng Baomin. Study on path planning and trajectory tracking control of free-floating space robot. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2008 (in Chinese)
[6] 胡仄虹, 黄攀峰, 孟中杰等. 空间绳系机器人逼近过程的位姿一体化控制. 航空学报, 2013, 34(11): 2635-2644 (Hu Zehong, Huang Panfeng, Meng Zhongjie, et al. Integrated pose control of tethered space robot in approaching process. Acta Aeronautica et Astronautica Sinica, 2013, 34(11): 2635-2644 (in Chinese) Hu Zehong, Huang Panfeng, Meng Zhongjie, et al. Integrated pose control of tethered space robot in approaching process. Acta Aeronautica et Astronautica Sinica, 2013, 34(11): 2635-2644 (in Chinese)
[7] 王剑颖. 航天器姿轨一体化动力学建模、控制与导航方法研究[博士论文]. 哈尔滨: 哈尔滨工业大学, 2014 (Wang Jianying. Research on spacecraft integrated orbit and attitude dynamics, control and navigation. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2014 (in Chinese) Wang Jianying. Research on spacecraft integrated orbit and attitude dynamics, control and navigation. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2014 (in Chinese)
[8] Dooley J, McCarthy JM. Spatial rigid body dynamics using dual quaternion components//Proceedings of 1991 IEEE International Conference on Robotics and Automation, IEEE, 1991: 90-95
[9] Brodsky V, Shoham M. Dual numbers representation of rigid body dynamics. Mechanism and Machine Theory, 1999, 34(5): 693-718
[10] Chen L, Zielinska T, Wang J, et al. Solution of an inverse kinematics problem using dual quaternions. International Journal of Applied Mathematics and Computer Science, 2020, 30(2): 351-361
[11] 杨一岱, 荆武兴, 张召. 一种挠性航天器的对偶四元数姿轨耦合控制方法. 宇航学报, 2016, 37(8): 946-956 (Yang Yidai, JIing Wuxing, Zhang Zhao. A new method for orbit and attitude coupling control problem of flexible spacecraft based on dual quaternions. Journal of Astronautics, 2016, 37(8): 946-956 (in Chinese) Yang Yidai, JIing Wuxing, Zhang Zhao. A new method for orbit and attitude coupling control problem of flexible spacecraft based on dual quaternions. Journal of Astronautics, 2016, 37(8): 946-956 (in Chinese)
[12] 董宏洋. 基于对偶四元数的航天器位姿一体化控制方法研究. [博士论文]. 哈尔滨: 哈尔滨工业大学, 2019 (Dong Hongyang. Spacecraft integrated attitude-position control using dual quaternions. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2019 (in Chinese) Dong Hongyang. Spacecraft integrated attitude-position control using dual quaternions. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2019 (in Chinese)
[13] 张洪珠. 基于对偶四元数的航天器姿轨一体化动力学建模与控制. [硕士论文]. 哈尔滨: 哈尔滨工业大学, 2012 (Zhang Hongzhu. Integrated dynamics modeling and control for spacecraft based on dual quaternion. [Master Thesis]. Haerbin: Harbin Institute of Technology, 2012 (in Chinese) Zhang Hongzhu. Integrated dynamics modeling and control for spacecraft based on dual quaternion. [Master Thesis]. Haerbin: Harbin Institute of Technology, 2012 (in Chinese)
[14] Benic J, Kasac J, Situm Z, et al. System modeling and control of 2dof robotic manipulator based on dual quaternion approach//2022 International Conference on Electrical, Computer and Energy Technologies (ICECET), 2022
[15] Antonello A, Valverde A, Tsiotras P. Dynamics and control of spacecraft manipulators with thrusters and momentum exchange devices. Journal of Guidance, Control, and Dynamics, 2019, 42(1): 15-29
[16] Valverde A, Tsiotras P. Dual quaternion framework for modeling of spacecraft-mounted multibody robotic systems. Frontiers in Robotics and AI, 2018, 5: 128
[17] Qi L, Ling C, Yan H. Dual quaternions and dual quaternion vectors. Communications in Applied Mathematics and Computational Science, 2022, 4: 1494-1508
[18] Wang P, Wang T, Sun J, et al. Dual-quaternion-based dynamics modeling for a rigid–flexible coupling satellite. Journal of Guidance, Control, and Dynamics, 2023, 46(7): 1298-1313
[19] 彭璇. 基于对偶四元数建模的航天器位姿一体化控制方法研究. [博士论文]. 哈尔滨: 哈尔滨工业大学, 2021 (Peng Xuan. Research on spacecraft integrated position-attitude control methods based on dual quaternion modeling. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2021 (in Chinese) Peng Xuan. Research on spacecraft integrated position-attitude control methods based on dual quaternion modeling. [PhD Thesis]. Haerbin: Harbin Institute of Technology, 2021 (in Chinese)
[20] Nehrenz M, Sorgenfrei M. A comparison of thruster implementation strategies for a deep space nanosatellite//AIAA Guidance, Navigation, and Control Conference, AIAA Paper, 2015, 2015-0866
[21] Curti F, Romano M, Bevilacqua R. Lyapunov-based thrusters’ selection for spacecraft control: Analysis and experimentation. Journal of Guidance, Control, and Dynamics, 2010, 33(4): 1143-1160
[22] 黄志来, 李新圆, 金栋平. 单框架控制力矩陀螺输出特性分析. 力学学报, 2021, 53(2): 511-523 (Huang Zhilai, Li Xinyuan, Jin Dongping. Output characteristic analysis of single gimbal control moment gyroscope. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 511-523 (in Chinese) Huang Zhilai, Li Xinyuan, Jin Dongping. Output characteristic analysis of single gimbal control moment gyroscope. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 511-523 (in Chinese)
[23] 汤亮, 徐世杰. 采用变速控制力矩陀螺的航天器自适应姿态跟踪和稳定控制研究. 航空学报, 2006, 27(4): 663-669 (Tang Liang, Xu Shijie. Spacecraft adaptive attitude tracking and stable control with variable speed control moment gyroscopes. Acta Aeronautica et Astronautica Sinic, 2006, 27(4): 663-669 (in Chinese) Tang Liang, Xu Shijie. Spacecraft adaptive attitude tracking and stable control with variable speed control moment gyroscopes. Acta Aeronautica et Astronautica Sinic, 2006, 27(4): 663-669 (in Chinese)
[24] 黄志来, 李新圆, 金栋平. 单框架控制力矩陀螺输出特性分析. 力学学报, 2021, 532: 511-523 (Huang Zhilai, Li Xinyuan, Jin Dongping. Output characteristic analysis of single gimbal control moment gyroscope. Chinese Journal of Theoretical and Applied Mechanics, 2021, 532: 511-523 (in Chinese) doi: 10.6052/0459-1879-20-306 Huang Zhilai, Li Xinyuan, Jin Dongping. Output characteristic analysis of single gimbal control moment gyroscope. Chinese Journal of Theoretical and Applied Mechanics, 2021, 532: 511-523 (in Chinese) doi: 10.6052/0459-1879-20-306
[25] Papakonstantinou C, Lappas V, Kostopoulos V, et al. A gimballed control moment gyroscope cluster design for spacecraft attitude control. Aerospace, 2021, 8(9): 273
[26] Huang XH, Jia YH, Xu S, et al. A new steering approach for VSCMGs with high precision. Chinese Journal of Aeronautics, 2016, 29(6): 1673-1684
[27] Wie B. Singularity escape/avoidance steering logic for control moment gyro systems. Journal of Guidance, Control, and Dynamics, 2005, 28: 948-956
[28] 贾英宏, 徐世杰. 采用变速控制力矩陀螺的一种姿态/能量一体化控制研究. 宇航学报, 2003, 24(1): 32-37 (Jia Yinghong, Xu Shijie. Study on integrated attitude/power control using variable speed control moment gyros. Journal of Astronautics, 2003, 24(1): 32-37 (in Chinese) Jia Yinghong, Xu Shijie. Study on integrated attitude/power control using variable speed control moment gyros. Journal of Astronautics, 2003, 24(1): 32-37 (in Chinese)
[29] Yang T. Optimal thruster selection with robust estimation for formation flying applications. [PhD Thesis]. Cambridge: Massachusetts Inst. of Technology, 2003
[30] 陈家睿. 七轴机械臂的轨迹规划. [硕士论文]. 哈尔滨: 哈尔滨工业大学, 2022 (Chen Jiarui. Trajectory optimization for 7-dof kinematically redundant arms. [Master Thesis]. Haerbin: Harbin Institute of Technology, 2022 (in Chinese) Chen Jiarui. Trajectory optimization for 7-dof kinematically redundant arms. [Master Thesis]. Haerbin: Harbin Institute of Technology, 2022 (in Chinese)
[31] Sharkawy AN, Aspragathos N. A comparative study of two methods for forward kinematics and Jacobian matrix determination//2017 International Conference on Mechanical, System and Control Engineering (ICMSC). IEEE, 2017: 179-183
[32] Aspragathos NA, Dimitros JK. A comparative study of three methods for robot kinematics. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 1998, 28(2): 135-145 doi: 10.1109/3477.662755
-
期刊类型引用(3)
1. 王传达,李飞,章子健,周云,周文,彭海军. 直升机复杂外形桨叶的模块化多体建模方法. 力学学报. 2025(05): 1202-1215 . 本站查看
2. 付江维,范志文,吕书锋,宋晓娟. 充液航天器姿轨耦合预设性终端滑模控制. 航天控制. 2024(05): 45-52 . 百度学术
3. 李婧涵,周荻,李思远,杜润乐,刘佳琪. 非质心推进航天器近距离接近姿轨一体化控制. 宇航学报. 2024(12): 1959-1973 . 百度学术
其他类型引用(2)