接触是机器人与外部世界交互的根本方式。无论是人形机器人的行走、机械臂的抓取，还是四足机器人的奔跑，都离不开对接触力的精确建模与控制。本章深入探讨接触力学的核心理论，从库仑摩擦模型到互补性问题求解，再到接触稳定性分析，为后续的接触规划和全身控制奠定坚实的理论基础。我们将特别关注如何在计算效率和物理真实性之间取得平衡，这是实时机器人控制的关键挑战。
在机器人学中，我们通常将接触简化为点接触模型。这种简化虽然忽略了接触面积和压力分布，但在大多数控制应用中提供了良好的精度-复杂度平衡。考虑机器人某个链接与环境表面的接触点，定义接触坐标系，其中z轴沿接触法向量指向机器人，x-y平面为切平面。
接触力 $\mathbf{f} = [f_x, f_y, f_z]^T$ 可分解为：
单边接触约束（Signorini条件）确保接触力只能是压力： \(f_n \geq 0\)
这个简单的不等式约束是接触问题非凸性的根源，也是使接触动力学问题在数学上具有挑战性的核心原因。当 $f_n = 0$ 时，接触断开；当 $f_n > 0$ 时，接触激活。这种离散的开关行为导致了系统的混合动力学特性。
接触运动学约束： 当接触激活时，接触点的法向速度必须为零（假设刚性接触）： \(v_n = 0 \quad \text{if} \quad f_n > 0\)
这与单边约束一起构成了完整的法向接触条件。
库仑摩擦定律是1785年由Charles-Augustin de Coulomb提出的经验定律，尽管其简单性掩盖了摩擦的微观复杂性，但它在工程应用中极其有效。该定律描述了切向力与法向力的关系，将可行接触力限制在一个锥形区域内。
静摩擦情况： 当接触点无相对滑动时，切向力必须满足： \(\|\mathbf{f}_t\| \leq \mu_s f_n\)
其中 $\mu_s$ 是静摩擦系数，通常在0.3-1.5范围内（取决于材料对）。这定义了一个二阶锥（Second-Order Cone）： \(\mathcal{K} = \{\mathbf{f} \in \mathbb{R}^3 : \sqrt{f_x^2 + f_y^2} \leq \mu_s f_z, f_z \geq 0\}\)
几何上，这是一个顶点在原点、半锥角为 $\theta = \arctan(\mu_s)$ 的圆锥：
z (法向)
^
|
| f_n
____|____
/ | \
/ | \ 摩擦锥
/ |θ \ tan(θ) = μ_s
/ | \
--------|---------> y
| /
| /
| /
| /
x (切向)
动摩擦情况： 当发生滑动时（$|\mathbf{v}_t| > 0$），切向力具有确定的大小和方向： \(\mathbf{f}_t = -\mu_d f_n \frac{\mathbf{v}_t}{\|\mathbf{v}_t\|}\)
其中 $\mu_d$ 是动摩擦系数（通常 $\mu_d < \mu_s$），$\mathbf{v}_t$ 是接触点的相对切向速度。这个关系表明：
Stribeck效应： 实际摩擦在静-动转换时表现出复杂行为： \(\mu(\|\mathbf{v}_t\|) = \mu_d + (\mu_s - \mu_d) \exp(-\|\mathbf{v}_t\|/v_{stribeck})\)
其中 $v_{stribeck}$ 是Stribeck速度，典型值约0.001-0.01 m/s。这种速度依赖性在精密控制中很重要，但在大多数机器人应用中可以忽略。
二阶锥约束虽然精确，但在优化问题中处理较为复杂，需要专门的锥规划求解器。实践中，特别是在实时控制中，常用多面体锥近似将问题转化为线性规划或二次规划。
多面体近似原理： 将圆形摩擦锥用 $m$ 个半平面的交集近似，每个半平面对应一个摩擦锥的生成方向。将摩擦锥离散化为 $m$ 个均匀分布的方向： \(\mathbf{d}_i = \frac{1}{\sqrt{1 + \mu^2}}[\mu\cos(2\pi i/m), \mu\sin(2\pi i/m), 1]^T, \quad i = 1, ..., m\)
这些向量是摩擦锥边缘的单位方向。接触力可表示为这些方向的非负线性组合： \(\mathbf{f} = \sum_{i=1}^m \lambda_i \mathbf{d}_i, \quad \lambda_i \geq 0\)
近似精度分析：
近似误差定义为： \(\epsilon = \max_{\mathbf{f} \in \mathcal{K}_{poly}} \frac{d(\mathbf{f}, \mathcal{K})}{\|\mathbf{f}\|}\)
其中 $\mathcal{K}_{poly}$ 是多面体锥，$d(\cdot, \cdot)$ 是点到集合的距离。
计算优势： 线性化后，接触力优化问题变为： \(\min_{\boldsymbol{\lambda}} f(\boldsymbol{\lambda}) \quad \text{s.t.} \quad \boldsymbol{\lambda} \geq \mathbf{0}\)
这是标准的非负最小二乘问题，可用高效的QP求解器（如qpOASES、OSQP）在微秒级时间内求解。
当机器人与环境有多个接触点时（如四足机器人的四只脚、人形机器人的双脚、或多指抓取），需要考虑所有接触力的联合可行域。这个高维空间的结构对理解系统的力能力至关重要。
多接触力的笛卡尔积： 考虑 $n$ 个接触点，每个接触点 $i$ 有其摩擦锥 $\mathcal{K}_i$。整体接触力的可行域是各接触锥的笛卡尔积： \(\mathcal{F} = \mathcal{K}_1 \times \mathcal{K}_2 \times ... \times \mathcal{K}_n \subset \mathbb{R}^{3n}\)
这意味着每个接触力独立地受其摩擦锥约束。
可行域的几何性质： 这个可行域具有以下重要性质：
非凸性：由于单边约束 $f_{ni} \geq 0$，当某个接触断开时（$f_{ni} = 0$），可行域发生不连续变化。这导致优化问题可能有多个局部最优。
多面体结构：在摩擦锥线性化后，$\mathcal{F}$ 成为多面体锥： \(\mathcal{F}_{poly} = \{\mathbf{f} : \mathbf{A}\mathbf{f} \leq \mathbf{0}, \mathbf{f} \geq \mathbf{0}\}\) 其中 $\mathbf{A}$ 编码了所有摩擦锥约束。
维度爆炸：随接触点数量增长，可行域维度呈线性增长（$3n$ 维），但约束数量可能呈超线性增长（每个锥 $m$ 个约束）。
稀疏结构：约束矩阵 $\mathbf{A}$ 具有块对角结构，每个块对应一个接触点，这可被高效求解器利用。
扳手空间的可行域： 通过抓取矩阵 $\mathbf{G} \in \mathbb{R}^{6 \times 3n}$ 将接触力映射到机器人质心的扳手（力和力矩）： \(\mathbf{w} = \mathbf{G}\mathbf{f} = \sum_{i=1}^n \mathbf{G}_i \mathbf{f}_i\)
扳手可行域： \(\mathcal{W} = \mathbf{G}(\mathcal{F}) = \{\mathbf{w} : \exists \mathbf{f} \in \mathcal{F}, \mathbf{w} = \mathbf{G}\mathbf{f}\}\)
这是 $\mathcal{F}$ 在 $\mathbf{G}$ 下的像，通常是6维空间中的非凸集。
计算复杂度考虑：
接触的本质特征——非贯穿、单边性、和摩擦——可用互补性条件精确而优雅地描述。互补性理论源于数学规划，完美捕捉了接触的”非此即彼”特性：要么有间隙无力，要么有力无间隙。
法向互补性（Signorini条件，1933）： 对于每个潜在接触点，法向力 $f_n$ 和法向间隙 $\delta_n$ 满足： \(0 \leq f_n \perp \delta_n \geq 0\)
其中 $\perp$ 表示互补性，即 $f_n \cdot \delta_n = 0$。这个简洁的数学表达包含三个约束：
物理含义的完整解释：
切向互补性（摩擦条件）： 摩擦的互补性更为复杂，涉及力-速度关系。定义滑动速度 $\mathbf{v}t$ 和最大静摩擦力 $f{max} = \mu f_n$：
\[\begin{cases} \|\mathbf{v}_t\| = 0 \Rightarrow \|\mathbf{f}_t\| \leq f_{max} & \text{（粘滞/静摩擦）} \\ \|\mathbf{v}_t\| > 0 \Rightarrow \mathbf{f}_t = -f_{max} \frac{\mathbf{v}_t}{\|\mathbf{v}_t\|} & \text{（滑动/动摩擦）} \end{cases}\]这可以用最大耗散原理统一表述： \(\mathbf{f}_t \in \arg\min_{\|\mathbf{f}'_t\| \leq \mu f_n} \mathbf{v}_t^T \mathbf{f}'_t\)
速度级互补性： 在速度级别（常用于仿真），互补性条件变为： \(0 \leq f_n \perp v_n + \epsilon \delta_n \geq 0\)
其中 $\epsilon > 0$ 是稳定化参数（Baumgarte稳定化），防止约束漂移。
能量一致性： 互补性条件确保接触不产生能量： \(\mathbf{f}^T \mathbf{v} = f_n v_n + \mathbf{f}_t^T \mathbf{v}_t \leq 0\)
等号仅在完全非弹性碰撞时成立。
线性互补性问题（Linear Complementarity Problem）提供了求解接触动力学的数学框架。通过时间离散化和适当的假设，复杂的接触动力学可以转化为标准的LCP形式。
时间离散化： 考虑离散时间步长 $h$，从时刻 $t$ 到 $t+h$。使用半隐式Euler积分：
动力学方程（Newton-Euler方程的离散形式）： \(\mathbf{M}\mathbf{v}^{+} = \mathbf{M}\mathbf{v}^{-} + h\mathbf{f}_{ext} + \mathbf{J}^T\boldsymbol{\lambda}\)
其中：
接触速度约束： 接触点的速度由广义速度通过雅可比映射得到： \(\mathbf{v}_c = \mathbf{J}\mathbf{v}^{+}\)
LCP推导： 将动力学方程解出 $\mathbf{v}^{+}$： \(\mathbf{v}^{+} = \mathbf{v}^{-} + h\mathbf{M}^{-1}\mathbf{f}_{ext} + \mathbf{M}^{-1}\mathbf{J}^T\boldsymbol{\lambda}\)
代入接触速度约束： \(\mathbf{v}_c = \mathbf{J}(\mathbf{v}^{-} + h\mathbf{M}^{-1}\mathbf{f}_{ext}) + \mathbf{J}\mathbf{M}^{-1}\mathbf{J}^T\boldsymbol{\lambda}\)
这给出标准LCP形式： \(\mathbf{v}_c = \mathbf{A}\boldsymbol{\lambda} + \mathbf{b}\) \(0 \leq \boldsymbol{\lambda} \perp \mathbf{v}_c \geq 0\)
其中：
Delassus矩阵的性质：
存在性与唯一性：
求解LCP是接触仿真的计算瓶颈。不同算法在精度、速度和鲁棒性之间有不同权衡。
Lemke算法（主元法）： Lemke算法是求解LCP的经典直接法，基于线性规划的单纯形法思想。
算法特点：
实现要点：
投影Gauss-Seidel（PGS）： PGS是接触仿真中最流行的迭代算法，因其简单性和鲁棒性。
迭代公式（分块形式）： \(\lambda_i^{k+1} = \text{proj}_{[0,\infty)}\left(\lambda_i^k - \omega \frac{1}{A_{ii}}\left(\sum_{j<i} A_{ij}\lambda_j^{k+1} + \sum_{j>i} A_{ij}\lambda_j^k + b_i\right)\right)\)
其中 $\omega \in (0, 2)$ 是超松弛因子。
收敛性分析：
实践优化：
加速方法：
算法选择指南：
当考虑库仑摩擦时，问题变为NCP。摩擦锥约束是非线性的：
最大耗散原理： 接触力应最小化功率耗散： \(\boldsymbol{\lambda} = \arg\min_{\boldsymbol{\lambda}' \in \mathcal{K}} \mathbf{v}_c^T \boldsymbol{\lambda}'\)
锥互补性问题（CCP）： \(\mathcal{K} \ni \boldsymbol{\lambda} \perp \mathbf{v}_c \in \mathcal{K}^*\)
其中 $\mathcal{K}^*$ 是对偶锥。
求解策略：
隐式Euler方法： 稳定但数值阻尼大： \(\mathbf{q}^{t+h} = \mathbf{q}^t + h\mathbf{v}^{t+h}\) \(\mathbf{M}\mathbf{v}^{t+h} = \mathbf{M}\mathbf{v}^t + h(\mathbf{f}_{ext} + \mathbf{f}_c^{t+h})\)
半隐式方法： 平衡稳定性和精度： \(\mathbf{v}^{t+h} = \mathbf{v}^t + h\mathbf{M}^{-1}(\mathbf{f}_{ext} + \mathbf{f}_c^{t+h})\) \(\mathbf{q}^{t+h} = \mathbf{q}^t + h\mathbf{v}^{t+h}\)
变分积分器： 保持能量和动量守恒： \(\delta \int_0^T L(\mathbf{q}, \dot{\mathbf{q}}) dt = 0\)
力封闭（Force Closure）： 接触配置能够抵抗任意外力/力矩： \(\exists \boldsymbol{\lambda} \in \mathcal{K} : \mathbf{G}\boldsymbol{\lambda} = -\mathbf{w}_{ext}\)
其中 $\mathbf{G}$ 是抓取矩阵，$\mathbf{w}_{ext}$ 是外部扳手。
检验条件： \(\text{rank}(\mathbf{G}) = 6 \text{ 且 } \mathbf{0} \in \text{int}(\text{ConvexHull}(\mathbf{G}\mathcal{K}))\)
形封闭（Form Closure）： 更强的条件，纯几何约束阻止所有运动： \(\mathbf{G}\boldsymbol{\lambda} = \mathbf{0}, \boldsymbol{\lambda} \in \mathcal{K} \Rightarrow \boldsymbol{\lambda} = \mathbf{0}\)
对于点接触，至少需要7个接触点才能实现形封闭。
实际接触并非刚性，存在局部变形。Hertz接触模型：
法向刚度： \(f_n = k_n \delta_n^{3/2}\)
其中 $k_n$ 取决于材料属性和接触几何。
切向刚度（预滑动）： \(\mathbf{f}_t = \mathbf{K}_t \boldsymbol{\delta}_t\)
刚度矩阵： \(\mathbf{K}_c = \begin{bmatrix} k_t & 0 & 0 \\ 0 & k_t & 0 \\ 0 & 0 & k_n \end{bmatrix}\)
典型值：$k_n \sim 10^6$ N/m，$k_t \sim 0.5k_n$
Hunt-Crossley模型包含非线性阻尼： \(f_n = k_n \delta_n^{3/2}(1 + \frac{3(1-e^2)}{4}\dot{\delta}_n)\)
其中 $e$ 是恢复系数。
线性化阻尼模型： \(\mathbf{f}_c = -\mathbf{K}_c\boldsymbol{\delta} - \mathbf{D}_c\dot{\boldsymbol{\delta}}\)
阻尼比：$\zeta = \frac{d}{2\sqrt{km}}$，典型取 $\zeta \in [0.1, 0.3]$。
考虑接触系统的稳定性，定义Lyapunov函数：
机械能函数： \(V = \frac{1}{2}\mathbf{v}^T\mathbf{M}\mathbf{v} + U(\mathbf{q})\)
其中 $U(\mathbf{q})$ 是势能（包括重力势和弹性势）。
耗散不等式： \(\dot{V} \leq -\alpha V - \beta \|\mathbf{v}\|^2\)
其中 $\alpha, \beta > 0$ 由摩擦和阻尼决定。
接触Lyapunov函数： 包含接触惩罚项： \(V_c = V + \frac{1}{2}\sum_i k_i \delta_i^2\)
稳定性条件：
参数不确定性： 摩擦系数、接触位置、表面法向的不确定性： \(\mu \in [\mu_{min}, \mu_{max}]\) \(\mathbf{p}_c = \mathbf{p}_{nom} + \Delta\mathbf{p}, \|\Delta\mathbf{p}\| \leq \epsilon\)
鲁棒力封闭： 在所有不确定性下保持力封闭： \(\forall \mu \in [\mu_{min}, \mu_{max}] : \mathbf{0} \in \text{int}(\mathcal{W})\)
稳定裕度： 定义最小抗扰动能力： \(\sigma = \min_{\|\mathbf{w}\|=1} \max_{\boldsymbol{\lambda} \in \mathcal{K}} \mathbf{w}^T\mathbf{G}\boldsymbol{\lambda}\)
$\sigma > 0$ 表示鲁棒力封闭。
人形机器人脚掌通常建模为矩形，四个角点接触：
脚掌俯视图
p1 -------- p2
| |
| C | C: 压力中心
| |
p3 -------- p4
x: 前进方向
y: 侧向
每个接触点 $i$ 的力：$\mathbf{f}i = [f{xi}, f_{yi}, f_{zi}]^T$
总接触扳手： \(\mathbf{w} = \sum_{i=1}^4 \begin{bmatrix} \mathbf{f}_i \\ \mathbf{p}_i \times \mathbf{f}_i \end{bmatrix}\)
压力中心： 所有接触力的等效作用点： \(\mathbf{p}_{CoP} = \frac{\sum_i f_{zi} \mathbf{p}_i}{\sum_i f_{zi}}\)
CoP必须在支撑多边形内： \(\mathbf{p}_{CoP} \in \text{ConvexHull}(\{\mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3, \mathbf{p}_4\})\)
零力矩点： 地面反力的等效作用点，水平力矩为零： \(\mathbf{p}_{ZMP} = \frac{\mathbf{n} \times (\mathbf{M}_O + \mathbf{p}_O \times \mathbf{F})}{\mathbf{n} \cdot \mathbf{F}}\)
其中 $\mathbf{n}$ 是地面法向量。
稳定性条件：$\mathbf{p}_{ZMP} \in$ 支撑多边形
给定期望的总扳手 $\mathbf{w}_{des}$，求解各接触点的力：
优化问题： \(\begin{align} \min_{\mathbf{f}_1, ..., \mathbf{f}_4} & \quad \sum_{i=1}^4 \|\mathbf{f}_i\|^2 \\ \text{s.t.} & \quad \sum_{i=1}^4 \mathbf{G}_i \mathbf{f}_i = \mathbf{w}_{des} \\ & \quad \mathbf{f}_i \in \mathcal{K}_i, \quad i = 1, ..., 4 \end{align}\)
其中 $\mathbf{G}_i$ 是第 $i$ 个接触点的抓取矩阵。
解析解（无摩擦情况）： 使用伪逆： \(\mathbf{f} = \mathbf{G}^+ \mathbf{w}_{des} + (\mathbf{I} - \mathbf{G}^+\mathbf{G})\mathbf{f}_0\)
其中 $\mathbf{f}_0$ 是内力（不产生净扳手）。
行走时的接触模式切换：
双支撑相：
单支撑相：
切换条件：
起脚条件：摆动脚法向力降至阈值 \(\sum_{i \in swing} f_{zi} < \epsilon_{\text{lift}}\)
落脚条件：摆动脚检测到接触 \(\delta_n < \epsilon_{\text{touch}} \text{ 且 } \dot{\delta}_n < 0\)
平滑过渡： 使用时变权重： \(\mathbf{f} = \alpha(t)\mathbf{f}_{\text{double}} + (1-\alpha(t))\mathbf{f}_{\text{single}}\)
其中 $\alpha(t)$ 是平滑过渡函数。
简化模型：
快速求解器：
传感器融合：
典型控制频率：1kHz（全身控制），100Hz（规划层）
本章系统地介绍了接触力学的核心概念和计算方法：
关键概念：
核心公式：
计算方法：
接触力学是机器人控制的基石，本章的理论将在后续的接触规划、MPC和全身控制中反复应用。
习题4.1 考虑一个质量为 $m = 10$ kg的箱子放在水平面上，静摩擦系数 $\mu_s = 0.5$。如果施加水平力 $F = 40$ N，箱子是否会滑动？计算接触力。
提示：比较施加力与最大静摩擦力。
习题4.2 给定4个接触点的摩擦锥，摩擦系数 $\mu = 0.7$。使用4面体近似，写出每个基向量 $\mathbf{d}_i$。
提示：考虑东南西北四个方向。
习题4.3 一个机器人手指与物体接触，接触点位于 $\mathbf{p} = [1, 0, 0]^T$ m，接触力 $\mathbf{f} = [0, 0, 10]^T$ N。计算该接触力产生的扳手。
提示：扳手包含力和力矩。
习题4.4 证明对于平面中的刚体，3个不共线的点接触（无摩擦）可以实现力封闭。需要多少个点才能实现形封闭？
提示：考虑平面刚体的自由度和约束数量。
习题4.5 推导PGS算法的收敛条件。证明当Delassus矩阵 $\mathbf{A}$ 是对称正定时，算法必然收敛。
提示：将PGS视为不动点迭代，分析谱半径。
习题4.6 设计一个算法，给定人形机器人质心位置和期望的地面反力，计算每个脚掌四个角点的接触力分配。考虑摩擦锥约束和ZMP稳定性约束。
提示：构造二次规划问题。
习题4.7（开放题）讨论软接触模型（如Hertz模型）与刚性接触模型在机器人控制中的优缺点。什么情况下必须使用软接触模型？
提示：考虑计算效率、物理真实性、数值稳定性。
陷阱：过度简化的摩擦锥（如4面体）可能导致不可行解。
示例：
# 错误：4面体可能过于保守
friction_cone_4 = [[1,0,μ], [0,1,μ], [-1,0,μ], [0,-1,μ]]
# 正确：使用8面体或更多
friction_cone_8 = []
for i in range(8):
angle = 2*π*i/8
friction_cone_8.append([μ*cos(angle), μ*sin(angle), 1])
解决：根据精度需求选择合适的离散化程度。
陷阱：严格的互补性 $f_n \cdot \delta_n = 0$ 在数值上难以满足。
正确做法： \(f_n \cdot \delta_n \leq \epsilon_{comp}\)
其中 $\epsilon_{comp} \sim 10^{-6}$ 是容差。
常见原因：
调试技巧：
陷阱：仅在离散时刻检测接触可能遗漏快速碰撞。
解决：连续碰撞检测（CCD）或自适应时间步长。
问题：多接触时存在无穷多解。
错误做法：
# 可能产生巨大内力
f = pinv(G) @ w_des # 伪逆解
正确做法： 添加正则化或物理意义的目标函数（如最小化接触力）。
陷阱：接触模式突变导致控制不连续。
解决：
误区：追求高精度物理仿真而牺牲实时性。
实践经验：
记住：在机器人控制中，一个及时的近似解往往胜过延迟的精确解。