摩擦锥与软接触

球受到重力，并由平面法向力支撑。若法向加速度为零，则 f_n=mg=19.62\,\mathrm N。现在沿水平面施加外力 F_t。保持静止需要接触摩擦力与它相抵：

静摩擦力随平衡需求取值。当所需合力到达边界时，静止平衡达到摩擦能力的极限。发生滑动后，各向同性理想 Coulomb 摩擦取

本章用同一个 \mu 描述静、滑两类条件。若材料采用不同的静摩擦与滑动摩擦系数，可以分别用 \mu_s 与 \mu_k 代入对应关系。

质量为 m 的方块静止在倾角 \gamma 的平面上时，法向支撑为 mg\cos\gamma，沿斜面的摩擦需求为 mg\sin\gamma。两者的比值为 \tan\gamma；在 0\leq\gamma<\pi/2 时，静止条件就是 \tan\gamma\leq\mu。

在速度恰为零的瞬间，是否继续静止由后续动力学和平衡条件共同决定。速度条件、允许力范围与动力学方程合在一起，才构成完整的摩擦规律。

在三维接触面上，滑动可以指向两个独立切向的任意组合。先以法线和两个切向建立局部坐标，再把所有允许的接触力收集到力空间中；摩擦锥描述的就是这组允许力。

接触标架与各向同性圆锥

把接触轴排成矩阵的行：

设接触点两侧的世界速度分别为 u_1,u_2。接触速度与作用在几何 2 上的世界力为

几何 1 承受相反的力。v_n>0 表示分离，f_n\geq0 表示沿法线推开两物体。各向同性切向摩擦域为

固定 f_n 得到半径为 \mu f_n 的圆盘。法向力加倍，截面半径也加倍，因此全部截面组成圆锥。f_n=0 时截面收为原点。

若法向支撑为 20\,\mathrm N、滑动系数为 0.5，切向合力的大小最多为 10\,\mathrm N。两切向分量分别为 (6,8)\,\mathrm N 时合力恰好到达边界；分别把两个分量都取成 10\,\mathrm N 会使合力超出允许范围。圆盘把两个切向方向之间共享的摩擦容量表达出来。

椭圆摩擦锥与各向异性

接触面的两个切向方向可以具有不同的摩擦能力。例如，沿表面纹理与横跨纹理滑动时，允许的摩擦力尺度可能不同。用 \mu_1,\mu_2 分别描述两个主方向的滑动摩擦系数。在这个局部标架中，把圆形摩擦域沿两轴分别伸缩，得到椭圆摩擦锥（elliptic friction cone）：

固定一个正的法向力，分别除以 f_n^2，就得到熟悉的椭圆方程。它的两个半轴是 \mu_1f_n 和 \mu_2f_n。系数决定两个方向的相对尺度，法向力决定整个截面的大小。两系数相等时，椭圆成为圆，式（14.1a）便回到式（14.1）。

取 f_n=20\,\mathrm N、\mu_1=0.8,\mu_2=0.4，两个半轴分别为 16\,\mathrm N 和 8\,\mathrm N。切向力 (9.6,6.4)^{\mathsf T}\,\mathrm N 恰在边界上，因为 (9.6/16)^2+(6.4/8)^2=0.36+0.64=1。两个分量若分别取到各自半轴端点，平方和便等于二，已经超出允许域。两个方向共同满足一个椭圆条件。

若关心沿某一实际方向能够承受多大的力，可以把方向先固定。设该方向与第一切向轴的夹角为 \alpha，令 f_t=\ell(\cos\alpha,\sin\alpha)^{\mathsf T}。代入椭圆式得到

令 D=\operatorname{diag}(\mu_1,\mu_2)，同一摩擦域可以简记为

该写法还覆盖零摩擦方向。若 \mu_2=0，就有 f_{t2}=0，正法向截面退化为第一切向轴上的线段；若两系数都为零，只剩法向力射线。无论系数如何，f_n=0 的截面都只有零力。后面的凸性证明、最大耗散和软接触求解均以这组允许力为出发点。

法向支撑越大，允许的摩擦力范围越大。下图用竖向的法向力和水平的两个摩擦分量组成力空间，截面上的每个点都是一组可能的切向力；它的几何位置与物体在世界中的位置分别具有不同含义。

默认法向力为二十牛顿，截面半轴为十六牛顿和八牛顿。把法向力减半，两个半轴也同时减半；保持法向力而只减小第二个摩擦系数，截面只沿第二切向轴收缩。允许域确定后，具体选取哪个点还要由平衡、滑动功率或软接触目标决定。

摩擦锥为凸集：两组允许力的任意凸组合仍是允许力。这个性质使上一章的凸优化条件可以直接用于摩擦接触，含零摩擦方向的退化情形也保留凸性。

两物体受到互为反向的接触力。把两侧功率相加，并利用正交坐标变换，得到

切向部分 P_t=f_t^{\mathsf T}v_t 是相对滑动的摩擦功率。给定法向力和切向速度，理想 Coulomb 滑动摩擦在允许域内使这项功率最小。

在法向力固定、切向速度非零时，功率最小的圆盘边界点位于速度反方向，大小为 \mu f_n。因此最小功率为 -\mu f_n\|v_t\|。

两个切向摩擦尺度不同时，先用 f_t=Dy 把椭圆变成圆盘，再在圆盘中最小化功率，最后变回物理力。变换往返各乘一次尺度，因此物理力的每个分量带有摩擦系数的平方。当 Dv_t\ne0 时，得到 f_t^*=-f_nD^2v_t/\|Dv_t\|，其切向功率为 -f_n\|Dv_t\|。两个摩擦系数的平方改变分量比例，滑动力因而一般不与速度共线。

将矩阵式按两个分量展开，可以直接看出不同摩擦尺度怎样改变方向：

分母是两方向共用的归一化长度。分子中的 \mu_j^2 改变每个方向的相对分配；当两系数相同，公共因子约去，便回到速度反方向上的各向同性公式。

例如，取 f_n=20\,\mathrm N、(\mu_1,\mu_2)=(0.8,0.4)、v_t=(1,1)^{\mathsf T}\,\mathrm{m/s}，则

静止时先计算维持平衡所需的摩擦力，再检查它是否位于允许域。滑动时则固定非零速度，在允许域中寻找功率最小的点。下图沿用同一个摩擦截面，可比较这两种选取规则。

各向异性下，最小功率点通常偏离速度的反向延长线，但与速度的内积仍为负。判断耗散应检查内积；判断静止能否维持则检查所需平衡力是否可行。

MuJoCo 的 condim 决定接触包含哪些力与力矩。前三项是力，单位为牛顿；后三项是力矩，单位为牛顿米。数组顺序如下：

六维接触采用力在前、力矩在后的存储顺序。第五章的空间力按 [\tau;f] 排列；在两种接口之间传递时按各自顺序组装。固定朝向的平移球使用三个平移分量；扭转和滚动分量对应可转动刚体的一般接触块。

令 r=d-1，把法向以外的摩擦分量组成 \phi\in\mathbb R^r，摩擦尺度组成 D=\operatorname{diag}(\mu_1,\ldots,\mu_r)。前两个尺度是无量纲滑动系数，后三个把法向力换为力矩，量纲为长度。统一锥条件为

滑动、扭转与滚动的共同摩擦域

在 condim="6" 时，把五个摩擦尺度分别写成 \mu_1,\mu_2,\mu_s,\mu_{r1},\mu_{r2}。前两项作用于切向力，第三项作用于绕法线的扭矩，后两项作用于绕切向轴的滚动力矩。尺度均为正时，统一锥条件展开为

每个商的单位都是牛顿：滑动力除以无量纲系数，摩擦力矩除以具有长度单位的尺度。固定法向力后，五个摩擦分量共享同一个椭球域。固定全部力矩为零，得到第二节的切向力椭圆；给定非零力矩后，留给切向力的范围相应减小。condim="4" 保留前三项摩擦分量；condim="3" 保留两个切向力。

棱锥的非负系数表示

棱锥为每个摩擦分量设置正负两条棱，系数 \alpha_{j+},\alpha_{j-}\geq0。将棱系数解码为力：

对于正摩擦尺度，非负棱系数可表示的力恰好满足 \sum_j|\phi_j|/\mu_j\leq f_n。必要性来自每对系数的差不超过和；充分性可通过把尚未分配的法向力平均加入各对系数来构造。

对于两个相同的摩擦系数，固定法向截面的圆盘和菱形分别为

在两个摩擦系数相同的截面中，沿切向坐标轴，菱形与圆盘达到相同半径；沿对角线，菱形允许的力大小为圆盘的 1/\sqrt2。因此棱锥近似会随切向方向改变摩擦容量。

把圆形或椭圆截面换成棱锥后，某些方向上的允许力会缩小。下面同时画出两种截面，并将选定的物理力编码成非负棱系数；它们的和与成对差分别恢复法向和摩擦分量。

棱表示让每个变量只受非负限制，代价是允许力域具有方向相关的边界。编码与解码的数值核对还应区分完整的非负构造和源码采用的具体分配规则。

多接触系统中，一个切向力的改变会影响其他方向的加速度。保持法向力暂时固定，可以先在一个椭圆截面内选择最合适的切向力。目标函数仍保留这些方向之间的耦合，因此需要同时考虑代价曲线与椭圆边界。

第十三章的式（13.6）给出了完整二次目标。本节固定法向力和其他接触的力，只留下当前接触的切向分量作为未知量，因而得到一个较小的二次问题。

将接触块的法向力暂时固定，摩擦变量记作 x=\phi。其余接触力和固定法向力的作用并入一次项 b_c，切向主子块记为 H_c\succ0，得到

这两个局部系数可从完整目标直接分组得到。把所有固定变量记作 z，当前摩擦变量记作 x，矩阵相应分块为 H_{xx},H_{xz},H_{zz}。含有 x 的部分为

因此 H_c=H_{xx}，b_c=b_x+H_{xz}z。固定的法向力也包含在 z 中；改变它会改变切向一次项和允许椭圆的大小。

本节先取所有摩擦尺度为正、f_n>0。这是二次约束二次规划（Quadratically Constrained Quadratic Program，QCQP）。它选择当前二次动力学代价最小的摩擦力；其一次项由当前方程残差构成。

先求无约束极小点，若它落在允许椭圆内就已得到答案。若在外部，则引入非负乘子，把边界最优条件转为一族正定线性方程；只需调节这个乘子，使解的缩放长度达到 f_n。

对需要边界解的问题，\psi(0)>0，而乘子充分大时 \psi<0。正定性保证 \psi 严格递减，所以根唯一。这个单调关系既解释了增大乘子会缩短力向量，也为区间保护提供了依据。

取 H_c=\operatorname{diag}(1,4)、b_c=(-2,-2)^{\mathsf T}、摩擦尺度 (1,1/2) 和 f_n=1。缩放后 Hessian 恰好成为单位阵，故边界乘子为 \sqrt5-1，可完整手算下面的过程。

当二次目标的无约束极小点落在椭圆外，最优允许力位于椭圆边界。图中的等值线在最优点与边界相切；切换到乘子曲线，则能看到同一条件怎样变成唯一的一维零点。

耦合项改变等值线方向，因此一般边界解会离开无约束解的径向。逐次增大乘子使缩放力长度缩短，最终达到给定法向力允许的半径；这一步构成椭圆接触 PGS 的切向更新。

当某个 \mu_j=0 时，先将 x_j=0 固定，再在其余分量上求解。f_n=0 时整个可行域只有原点，直接返回零力；此时二次约束在原点的梯度为零，有限的该形式乘子可能无法表示驻点条件。将锥尖单独处理保留了清楚的数学含义。

第十三章中的完整力空间问题可以直接采用摩擦锥可行域：

式（13.6）确定目标函数，式（14.1a）及其高维形式确定允许力域。将这两部分组合后，法向力与切向力需要一起满足最优条件。

这里法向力和摩擦分量一起变化。固定法向 QCQP 是它在一个切片上的局部更新。另一方面，当当前广义加速度已给定时，MuJoCo 用伪加速度残差 j=Ja-a_{\rm ref} 计算局部软响应：

把局部接触代价与惯性代价相加，再对加速度求驻点，得到的正是 Newton 动力学方程。与此同时，局部允许力条件还给出完整力空间问题的最优性条件，所以两种求解变量描述同一个软接触系统。

式（13.7）已给出独立法向行的加速度目标。这里把各行单独的非负限制合成一个摩擦锥，局部代价由式（14.6）统一生成；惯性项及其动力学含义保持相同。

为得到解析的锥响应，MuJoCo 调整切向正则量，使每个正摩擦尺度满足 R_j\mu_j^2=\kappa。记法向正则量为 R_n，并假定各个参与摩擦方向的尺度均为正；零摩擦方向按前述规则固定为零、从缩放计算中移除。

先令 f_t=Dy，把允许域改写成 \|y\|\leq f_n。由于各方向的 R_j\mu_j^2 都等于同一个数 \kappa，切向正则项随之变成 \kappa\|y\|^2/2。再把切向一次项写成 u^{\mathsf T}y，所需的新量只有

u 是按摩擦尺度加权后的切向残差，T 是它的长度。给定 f_n 后，最优的 y 沿 -u；剩下的选择只关乎它的大小以及法向力。源码还使用 \bar\mu=\sqrt{\kappa/R_n} 归一化法向尺度。

在本地实现中，设 impratio 为 \eta>0，首先设置 R_1=R_n/\eta，从而 \kappa=R_n\mu_1^2/\eta、\bar\mu=\mu_1/\sqrt\eta。contact.friction 保存物理摩擦尺度，contact.mu 保存 \bar\mu。

局部二次目标的极小点可能在锥尖、锥内或锥面上。锥尖给出零力；整空间极小点落在锥内时直接采用线性软响应；余下情形沿锥面把问题化为法向力的一元二次函数。三个分支及条件会由同一个目标逐一得到。

法向力也可以直接按同一分区读取：顶部取零，底部取 -j_n/R_n，中间取 (T-j_n)/(R_n+\kappa)。摩擦分量在底部为 -R_t^{-1}j_t，在中间为 -f_nD^2j_t/T。这里 R_t 是由各摩擦方向正则量组成的对角块。

固定软约束参数后，改变伪加速度残差会使最优力先落在锥尖，再进入锥面或锥内。下面在残差平面上显示这些区域，并沿相同输入画法向力或代价，检查分区公式在边界处怎样相接。

例如选 R_n=1、\kappa=1 并沿 T=1 的水平线改变 j_n，两个分区边界是 j_n=1 与 j_n=-1。在二者之间，法向力为 (1-j_n)/2；越过左侧边界后变为 -j_n。两段在 j_n=-1 都给出力 1，斜率则由 -1/2 变为 -1，所以力连续而曲率可以改变。

分区边界控制采用哪一个力公式，力本身连续变化。曲率在边界处可以改变，所以 Newton 法需要在新试探点重新识别区域。物体是否正在滑动还要同时查看实际切向速度。

三个区域给出的代价与力在边界处连续。Newton 求解还需要区域内部的 Hessian；在锥面区，它由一个向量外积和一个非负的切向投影项组成，因此仍为半正定。

理想滑动摩擦通过 v_t 最小化瞬时功率；软锥响应通过 j_t 与正则项确定力。顶部、底部和中间是残差代价的分区名称。实际粘着、滑动和分离要结合接触速度、几何间隙与系统运动判断。

椭圆锥选项与接触对参数

MuJoCo 通过 <option cone="elliptic"/> 选择椭圆锥，通过 cone="pyramidal" 选择棱锥。这里的 elliptic 也包含两个切向系数相等时的圆锥。锥类型确定允许力域的形状，condim 确定参与的力与力矩分量。

几何体的 geom friction 有三个数，依次给出滑动、扭转、滚动摩擦尺度；滑动值和滚动值分别复制到两个方向。显式接触对的 pair friction 则有五个数，依次对应两个滑动方向、一个扭转方向和两个滚动方向，可以分别设置各向异性尺度。参数顺序见 MuJoCo 接触对说明。

椭圆接触示例模型包含平面与可自由运动的方块，使用 cone="elliptic"、condim="3"，五个摩擦参数为 0.8 0.4 0.02 0.005 0.005。当前维数使用前两个数，对应第二节的椭圆；将维数改为四或六，可继续启用扭转或滚动力矩。各向异性的方向由实际接触标架确定，世界方向的解释应同时检查 contact.frame。实际接触参数经过 mj_assignFriction 时，各摩擦尺度会被限制到至少 mjMINMU；完全无摩擦的接触采用 condim="1"。第二节的零尺度截面用于说明数学上的退化极限。

接触表示与块求解

以下对应本地 MuJoCo 3.10.0。mjContact 的 frame 按法向、两个切向存放；friction[5] 依次为两个滑动、扭转、两个滚动摩擦尺度。dim 保存接触维数，mu 保存正则化锥参数，H 保存需要时计算的锥代价 Hessian。

engine_util_misc.c 中 mju_decodePyramid 直接累加棱系数，并对每对系数取差：

force[0] = 0;
for (int i=0; i<2*(dim-1); i++) {
  force[0] += pyramid[i];
}
for (int i=0; i<dim-1; i++) {
  force[i+1] = (pyramid[2*i]-pyramid[2*i+1])*mu[i];
}

同文件的 mju_encodePyramid 选择每一对棱相同的法向份额 a=f_n/r，源码中的差值为

在 |\phi_j/\mu_j|\leq a 的条件下，每对系数均非负且能恢复该分量。对于其他输入，这段具体实现的截取与份额选择会改变结果。例如 f_n=1、\mu_1=\mu_2=1、\phi=(0.8,0) 本来位于菱形内，源码编码再解码得到 (0.5,0)；取 \phi=(-0.8,0) 时，第一条棱系数为 -0.15。第四节的剩余量构造则在整个棱锥域上提供非负表示。

solPGS 对椭圆接触先作法向更新或沿当前整组接触力的射线更新，再固定新的法向力、更新摩擦分量。solveQCQP 在 dim=3 时调用 mju_QCQP2，在 dim=4 时调用 mju_QCQP3，在 dim=6 时调用五变量的通用 mju_QCQP。

mju_QCQP 的通用实现先把椭圆缩放为球，再执行最多 20 轮 Newton 更新。下列片段对应第五节的线性求解和导数：

mju_cholSolve(res, Ala, b, n);
mju_scl(res, res, -1, n);
val = mju_dot(res, res, n) - r*r;
mju_cholSolve(tmp, Ala, res, n);
deriv = -2.0*mju_dot(res, tmp, n);
la += -val/deriv;

调用层在约束活跃时还将摩擦结果径向缩放回椭圆边界，用于清除近似求解留下的径向误差。教学程序使用独立的区间保护和相对残差检查，便于观察每步乘子。

为了与源码的缩放变量对照，令 N=\bar\mu j_n。正文中的顶部条件 j_n\geq T 等价于 N\geq\bar\mu T，底部条件 \kappa j_n+R_nT\leq0 等价于 \bar\mu N+T\leq0。两种写法采用相同边界，源码形式把常用尺度预先合并。

mj_constraintUpdate_impl 用 U[0]=jar[i]*con->mu、U[j]=jar[i+j]*friction[j-1] 构造第六节的缩放量；三个区域分别写入 SATISFIED、QUADRATIC、CONE。同文件的正则量调整满足

R[i+1] = R[i]/mju_max(mjMINVAL, m->opt.impratio);
d->contact[id].mu = friction[0]*mju_sqrt(R[i+1]/R[i]);
R[i+j+1] = R[i+1]*friction[0]*friction[0]
             / (friction[j]*friction[j]);

本地 mj_assignFriction 把输入摩擦尺度限制到至少 mjMINMU=1e-5。理想零摩擦由 condim=1 表示；本章图示和教学程序也计算椭圆尺度恰为零的数学极限。efc_D 表示 R^{-1}，与本章表示摩擦尺度的矩阵 D 按上下文区分。

理论背景可进一步对照 MuJoCo 摩擦锥说明、软约束模型。源码数组和分支以上述本地版本为准。

下载 friction_contact.py，与第十三章的 convex_contact.py 放在同一目录。教学模块包含各向异性最大耗散、非负棱系数、固定法向 QCQP 与软锥响应。Python 用二分搜索求乘子，与图示的区间保护 Newton 形成独立计算对照。

python3 examples/friction_contact.py

B = [[mu[i]*H[i][j]*mu[j] for j in active]
     for i in active]
c = [mu[i]*b[i] for i in active]
# Given a multiplier, solve in the scaled force coordinates.
y = solve(add_diagonal(B, multiplier), [-v for v in c])
value = sum(v*v for v in y) - normal*normal
# value > 0: increase multiplier; value < 0: decrease it.

零尺度先删除对应自由变量，锥尖直接返回零力。一般正尺度问题检查椭圆可行性、乘子非负性、驻点残差与互补乘积；软锥问题另检查 f=-\nabla S 和区域边界的连续性。