凸优化与 KKT 条件

先只考虑一个能够推开物体的法向力。给定当前自由加速度，如果它已经让物体离开接触目标，就可以取零力；如果物体仍加速压向表面，则需要正的支撑力。下面把这个选择写成一个带非负边界的最小化问题。

式（12.3）已求出一条软双边行的响应。接触再要求力非负，负的无约束试算力便要改成零；这一条简单边界将引出后面的最优性条件。

设不含接触力时的加速度为 a_u，参考法向伪加速度为 a_{\rm ref}。定义零力残差 b=\boldsymbol n^{\mathsf T}a_u-a_{\rm ref}。对固定平面和平移球，法线保持不变，因此这一行的伪加速度也等于普通相对法向加速度。

软约束引入正则量 R\geq0。令 h=1/m+R>0，把正则化残差记作 w：

将这条线性关系作为代价对力的导数，就得到非负半轴上的二次问题：

非负半轴上的极小点为 f_n^*=\max(0,-b/h)。零力残差为负时，增加正法向力可以降低目标；零力残差非负时，允许区间的起点已经最优。

这三个条件也可直接从允许的微小力变化得到。若 f_n>0，力可以增大也可以减小，所以最优点要求导数 w=0。若 f_n=0，允许变化只有 \Delta f_n\geq0；为使代价的一阶变化 w\Delta f_n 非负，必须有 w\geq0。

三式称为互补条件。正的接触力对应零正则化残差；正残差对应零接触力。两者同时为零时，接触位于两类状态的交界处。实际加速度偏差仍由 j=w-Rf_n 给出。

以 a_{\rm ref}=0、a_u=(0,0,-9.81)^{\mathsf T}\,\mathrm{m/s^2} 和 R=0.1\,\mathrm{kg^{-1}} 为例，h=0.6\,\mathrm{kg^{-1}}、b=-9.81\,\mathrm{m/s^2}，故 f_n^*=16.35\,\mathrm N、w^*=0、j^*=-1.635\,\mathrm{m/s^2}。有限正则量允许参考加速度存在有限偏差。

一个接触何时承载，可以从抛物线最低点是否落在非负半轴看出。切换为导数图后，同一判断变成：允许区间内能否到达零导数，或者在零力边界处导数已经非负。

曲线、导数与互补条件描述同一个选择。内部最优点允许力向两侧微调，所以导数为零；边界点只允许增大力，因此最优时的右导数非负。

\Phi 的单位为 \mathrm{kg\,m^2/s^4}，由力乘加速度得到。它是选取约束力的优化代价；与机械能对应的焦耳单位各有用途。

单个力的代价是一条抛物线，最小点可以从斜率为零找到。多个力一起变化时，代价成为曲面：梯度描述当前位置的上升方向，Hessian 描述各方向的弯曲程度。这两个量将决定求解方向，也用于判断找到的驻点是否为全局最小值。

多个接触力组成列向量 f。对一般变量 x\in\mathbb R^n，梯度收集每个分量的一阶偏导，Hessian 矩阵（二阶导数矩阵）收集二阶偏导：

梯度给出小位移上的一阶变化：F(x+\Delta x)=F(x)+\nabla F(x)^{\mathsf T}\Delta x+O(\|\Delta x\|^2)。对于对称矩阵 H 构成的二次函数，这个展开可以精确写到二阶：

凸集包含任意两点之间的整条线段。凸函数在每条线段上位于端点连弦之下。精确地说，对 0\leq\theta\leq1：

沿任意方向切开二次曲面，得到一元二次函数。它的二阶导数为 z^{\mathsf T}Hz；所有方向的曲率非负，恰好保证整张曲面位于任意端点连弦下方。

Hessian 的符号可以通过实际截线理解。下图使用 H=\operatorname{diag}(1,k)：旋转截线方向时，绿色曲线的曲率是 z^{\mathsf T}Hz；铜色连弦检验凸函数的线段定义。

当 k<0 时，总能找到向下弯曲的方向；当 k=0 时存在平坦方向；当 k>0 时每个方向都向上弯曲。这个逐方向判断将用于确认接触目标有没有唯一极小点。

在凸可行域上的凸函数中，一个局部极小点已经达到全局最小值。原因是任意更低的可行点都会通过连接线段，在当前点附近也产生更低的函数值。

取目标点 p=(2,1)^{\mathsf T}，在直线 x_1+x_2=s 上寻找离它最近的点。距离平方的一半为

直线上的允许微小位移满足 \nabla h^{\mathsf T}\Delta x=0。在最优点，目标沿正、反两个切向方向的导数都为零，因此目标梯度沿直线法向。用乘子 \beta 写成

等式乘子可取任意实数。其符号表示约束法向的哪一侧需要与目标梯度平衡。

在这条直线上，最优乘子为 \beta=(3-s)/2，最优点为 ((1+s)/2,(s-1)/2)^{\mathsf T}。当 s=1 时得到点 (1,0)^{\mathsf T}，目标梯度与乘子项正好相反。

约束乘子把边界施加的限制写成梯度平衡。下图先用直线等式观察两支相反的梯度箭头，再改为半空间：目标点进入可行内部时，可以直接到达它，乘子随之降为零。

同一边界在等式模型中必须始终满足，在不等式模型中则允许进入内部。两种情况的差别由乘子的符号和互补条件表达，下一节将把它推广到多条边界。

把直线改为半空间 g(x)=x_1+x_2-s\leq0。其梯度指向可行域外侧。最优点位于边界时，目标梯度的反方向正好被边界挡住，于是

乘子的非负号可以用可行方向核对。边界上的 \nabla g 指向外侧，因此 d=-\nabla g 指向内部。若 \nabla F=-\alpha\nabla g，沿该方向的变化为 \nabla F^{\mathsf T}d=\alpha\|\nabla g\|^2。边界最优点要求这项非负，故 \alpha\geq0。

例如上节取 s=1 时，边界点 (1,0) 的负梯度指向被约束挡住的外侧。增大到 s=4 后，目标点 (2,1) 已位于允许内部，距离目标可以降为零，乘子也取零。

若最优点位于内部，附近各方向都可移动，必须有 \nabla F=0，并可取 \alpha=0。边界与内部统一为互补松弛 \alpha g(x)=0。本例的乘子为 \alpha^*=\max(0,(3-s)/2)。

对于多条不等式 g_i(x)\leq0 和仿射等式 h_j(x)=0，Lagrange 函数与四组 KKT 条件为

等号成立的不等式叫作活跃约束。活跃约束的乘子可以为零；互补关系允许这种边界交会状态。

若目标和不等式约束可微且凸、等式约束仿射，一组满足 KKT 的点与乘子足以证明全局最优。证明把原始目标夹在同一个 Lagrange 函数的两侧，依次使用驻点、互补和可行性。

要从一个最优点推出乘子的存在，还需约束资格条件。一个常用的充分条件是 Slater 条件：在目标和约束的共同定义域相对内部，存在满足仿射等式、并严格满足各凸不等式的点。结合有限且可达的原始最优值，它保证强对偶及最优乘子的存在。因此在这些条件下，KKT 同时给出必要条件与充分条件。

法向力的可行条件可写成 -f_i\leq0。在这一符号约定下，不等式乘子恰好等于正则化残差 w_i，四组一般 KKT 条件便收敛为接触的三组互补条件。

一个候选解给出了已经达到的代价，还需要判断是否可能存在更好的解。如果另一个计算能够给出所有可行解代价的下界，那么候选值与下界之间的差距，就能衡量剩余的优化空间。对偶函数用约束乘子构造这个下界。

固定 \alpha\geq0 和任意 \beta，让原始变量遍历整个定义域，得到对偶函数

对于每个原始可行点，非负乘子乘上非正约束值，得到一个非正修正。因此

右侧给出可行方案的代价，左侧给出任何方案都无法低于的下界。对偶问题在允许乘子中最大化该下界。即使原问题存在非凸部分，这条弱对偶不等式仍成立。强对偶表示最佳下界与原始最优值相等；上一节的凸性和 Slater 条件提供了一个常用保证。

对正定二次目标，固定非负乘子后可通过配方消去力。得到的对偶函数只含乘子；任意原始可行力与任意非负乘子都给出一个可计算的上下界差。差值为零即可同时确认两者最优。

下图先取 h=1、b=-1、\alpha=0，对偶下界为 -1/2。原始力取 f=0 时，目标为零，间隙为 1/2；移到 f=1 时，目标也降到 -1/2，间隙归零。分别改变力与乘子，就能比较上界和下界怎样接近。

在这个例子中，即使原始力已经到达最优值，任取一个较差的乘子仍会留下正的间隙。要让上下界相等，原始力与乘子需要同时满足驻点和互补关系；下面将这组条件用于多个接触。

单变量中 \mathcal D(\alpha)=-(b-\alpha)^2/(2h)。若 b<0，最优乘子为零，f_n^*=-b/h；若 b\geq0，最优乘子为 b，f_n^*=0。两种情况的对偶间隙都为零。

设 M\succ0 为质量矩阵，J 按约束行排列，a_u 为自由广义加速度。接触力经 J^{\mathsf T} 回传后得到

约束响应矩阵的物理来源见式（11.6），正则力律见式（12.2）。本节将它们合成为同一个力空间目标，并继续保留每行非负的接触条件。

将力引起的广义加速度先投影，再减去目标，得到 j=Af+b。这里的三个系数分别为

A 给出惯性响应，b 给出零力时的目标偏差，R 给出软约束正则量。加上正则项后，残差为 w=Hf+b，恰好是以下二次目标的梯度：

前两章的伪加速度约定在此继续使用：参考项与 Ja 采用相同的约束行定义。若要把它换成某个几何间隙的普通二阶时间导数，还应按该间隙的定义加入曲率项。

当每条行的正则量均为正时，H=A+R 正定，即使几何约束存在冗余也有唯一最优力。约束行独立时，硬约束极限仍保证力的唯一性；存在冗余时，还需结合允许力域与零力残差判断。

所有变量均为法向力时，可行域是非负正交象限 f\geq0。取各分量均为正的向量即可严格满足这些线性不等式，因而本章的正定问题具备 KKT 和强对偶所需条件。下一章将把可行域替换成摩擦锥。

用加速度表示同一个软接触问题

求解还可以从物体加速度出发。对一个候选 a，先计算第 i 行低于参考值的程度 j_i=J_i a-a_{{\rm ref},i}，再按式（12.2）的软关系得到 f_i=\max(0,-j_i/R_i)。所有法向行的代价加到第七章的惯性代价上，构成

第一项是式（7.4）的质量度量，第二项让低于参考值的接触加速度付出代价。这称为以加速度为变量的原始问题。式（13.6）则以接触力为变量，消去了广义加速度，称为力空间对偶问题。两者的极小点通过物体动力学一一对应；第十五章的 PGS 使用后者，CG 与 Newton 使用前者。

原始问题与对偶问题的名称取决于所选起点。第五节从力的非负二次问题出发，使用边界乘子构造下界；这里从加速度目标出发，接触力本身成为对偶变量。两组关系可以分别通过各自的拉格朗日函数核对。

继续采用第十一章的叠放方块。两物体质量各为一千克，只沿竖直方向运动，重力加速度大小取每二次方秒十米。地面反力记为第一行，层间反力记为第二行。现在对上方块施加向上的外力，观察上方块何时离开下方块。

外力大小记为 P。先取参考加速度和正则化为零，得到一个硬接触的可手算模型。自由加速度、接触响应与零力残差为

这里 H 的单位为每千克，b 的单位为加速度。两行独立，且 u^{\mathsf T}Hu=(u_1-u_2)^2+u_2^2>0，所以这个硬模型仍有唯一极小点。由 w=Hf+b，两行残差分别为 w_1=f_1-f_2-10 与 w_2=-f_1+2f_2+P。

上拉力引起的接触分离

取 P=15\,\mathrm N。若先假定两处都承载，就令两行残差为零。由第一行 f_1=10+f_2，代入第二行得 f_2=-5\,\mathrm N，再得 f_1=5\,\mathrm N。层间负力违反非负条件，说明“两处都承载”的假设需要调整。

直接把负分量截成零，会留下 f=(5,0)^{\mathsf T}\,\mathrm N。此时地面行残差为 w_1=-5\,\mathrm{m/s^2}：下方块仍在向地面加速。层间力改为零后，地面力也必须重新计算。由 f_2=0,w_1=0 得到 f_1=10\,\mathrm N；另一行随后给出 w_2=5\,\mathrm{m/s^2}，上方块开始分离。

活动集方法把“哪些力固定为零”作为一个候选假设，再求其余力并核对该假设。对于两个接触，只有四种组合，可以完整列举。

用一般记号表示：令 \mathcal P 收集允许为正的力分量，\mathcal Z 收集固定为零的分量。在所选面内，先解

再检查 f_{\mathcal P}\geq0 和 w_{\mathcal Z}\geq0。构造时已有 w_{\mathcal P}=0，三者合在一起就是完整的非负力 KKT 条件。此处“活动”指力域中被压在零点的非负边界；物理上承载的接触则具有正力，图中直接写出承载行以免混淆。

取 P=0，两处接触都承载，解为 (20,10)^{\mathsf T}\,\mathrm N。逐渐增加上拉力，在 0\leq P\leq10 范围内得到 f=(20-P,10-P)^{\mathsf T}\,\mathrm N；超过 10\,\mathrm N 后，层间力保持为零，地面力保持为 10\,\mathrm N。因而接触分离的边界可以从物体受力和最优性条件两侧互相核对。

加入正则量 R 后，图中使用 H=A+R，候选组合的检查方式保持不变。对 n 个非负变量，完整枚举有 2^n 种组合，适合小规模推导与独立校验。第十五章沿同一模型逐步更新力，说明大型求解器怎样逐渐满足这些条件。

本地头文件 mujoco.h 的版本宏为 mjVERSION_HEADER 3010000。力空间二次函数中的量与源码对应如下：

engine_forward.c 在 mj_fwdConstraint 中先形成零力残差：

mj_mulJacVec(m, d, d->efc_b, d->qacc_smooth);
mju_subFrom(d->efc_b, d->efc_aref, nefc);

engine_core_constraint.c 的 mj_makeAR 在约束空间惯性矩阵对角线上加 efc_R。同文件的 mj_constraintUpdate_impl 先计算软响应，再处理非负分支：

force[i] = -D[i]*jar[i];
// non-negative contact row
if (jar[i] >= 0) {
  force[i] = 0;
  state[i] = mjCNSTRSTATE_SATISFIED;
}

这段分支是在给定当前加速度时计算 f_n=\max(0,-j/R)。求解器收敛后，j=Af+b 与该响应同时成立，便得到 (A+R)f+b=w 的互补条件。j 和 w 因此在代码中承担不同的计算角色。

engine_solver.c 的 residual 计算整行 efc_b + efc_AR * efc_force。PGS 对标量法向行的核心操作是

force[i] -= res[0]*ARinv[c];
if (force[i] < 0) {
  force[i] = 0;
}

这里 ARinv[c] 存放当前行对角元素的倒数。保持其他力不变时，这一步精确极小化该坐标上的抛物线。全局耦合由后续遍历持续更新。

概念说明可对照 MuJoCo 的约束模型文档；凸性、Slater 条件和对偶理论可对照 Boyd 与 Vandenberghe 的 Convex Optimization 第 2、3、5 章。

下载 convex_contact.py。程序用高斯消元求每个候选面上的极小点，再检查非负性与互补。输入矩阵按嵌套行数组给出，返回力、残差及候选记录。

python3 examples/convex_contact.py

在 examples 目录可以直接运行下面的方块算例。候选记录保留每组承载假设，因此输出能够逐行对应第七节的表格。

from convex_contact import solve_contact

H = [[1.0, -1.0], [-1.0, 2.0]]
b = [-10.0, 15.0]
result = solve_contact(H, b)
print(result["force"])     # [10.0, 0.0]
print(result["residual"])  # [0.0, 5.0]
for row in result["candidates"]:
    print(row["support"], row["force"], row["residual"], row["valid"])

def candidate(H, b, support):
    n = len(b)
    f = [0.0] * n
    block = [[H[i][j] for j in support] for i in support]
    values = solve(block, [-b[i] for i in support])
    for i, value in zip(support, values):
        f[i] = value
    w = [sum(H[i][j]*f[j] for j in range(n)) + b[i]
         for i in range(n)]
    return f, w

枚举先给出有限个可核对的候选；接触数增加后可以用迭代算法替换搜索过程。保留 f\geq0、w\geq0 和 f_iw_i\approx0 作为统一检查。