本文主要针对了解数值优化、但主要集中在调库，不理解算法 原理的读者。（比如日常工作是定义一个带约束的优化问题进行求解，求解过程则全部交给优化器。）
本系列讲解了等式/不等式约束、罚函数、对数障碍，拉格朗日乘子法，约束最优性方程，KKT条件，对偶变量等概念，并以casadi和ipopt为例，讲解它如何把不等式约束转为等式约束并进行求解的过程。因此，如果想深入了解，可以阅读本文。
但是，读者需要能留出至少数个小时来阅读本文，否则心急的情况下确实是很难很快理解消化本文的内容。

本系列内容分为两章：
第一节见：【机器人感知运控笔记三】约束优化基础：从简单问题到拉格朗日乘子与 KKT 条件
包含机器人 案例引言及第 1—7 节：牛顿法回顾、消元、平方罚函数、对数障碍、乘子与互补条件。

第二节是本章：
理解 IPOPT：原始—对偶内点法与二阶信息
包含第 8—14 节：松弛变量、最优性方程、联合牛顿步骤、障碍参数更新、Hessian 接口与导数检查。

如果读完觉得不错，欢迎关注我的专栏：
机器人感知运控笔记-微鉴道长的专栏

注意：以下内容在开始前，必须先阅读前一章内容！！ CSDN字数的限制使我无法把这两节合并为一节，拆分成两节是有阅读困难的。。

8. 把障碍子问题写成一组方程

第 5 节已经给出了一种直接的求解方式。对贯穿全文的一维例子，固定 $\mu >0$ 后，只需要求解：

$$\underset{x<1}{\min }(x-2{)}^2-\mu \log (1-x).$$

这个例子可以直接求导，也容易找到满足 $x<1$ 的初值。增加变量和乘子会增加计算步骤，一维例子采用上面的写法已经很方便。

当问题包含很多变量和约束时，求解器还需要处理几个实际问题：初值可能违反部分约束；受到约束限制的最优点，目标梯度可能仍然非零；每次更新变量时，需要同时考虑目标和各条约束的影响。

本章介绍的写法为这些计算提供了一套统一的组织方式：用松弛变量单独表示约束余量，用乘子表达约束在梯度平衡中的作用，再把两者与原变量一起求解。 辅助量会增加未知数，也会提供明确的方程结构，供求解器计算每一轮的联合修正。

下面仍用 $f(x)=(x-2{)}^2$、$g(x)=x-1\le 0$ 这个例子，先说明每个辅助量的用途，再展开推导，最后介绍这套结构怎样用于多变量、多约束问题。

8.1 把余量列为变量，并用等式把它与 x 联系起来

原约束要求 $x\le 1$。当 $x=0.8$ 时，约束函数值为 $g(x)=-0.2$，到边界还剩 $0.2$。因此，本例的余量等于 $1-x=-g(x)$。

第 5 节用 $t=1-x$ 简写这个余量。现在把 $t$ 单独列为一个待求解的变量，同时用等式规定它与 $x$ 的关系：

$$x-1+t=0,t\ge 0.$$

这个 $t$ 称为松弛变量。等式要求 $t=1-x$，非负要求又保证 $x\le 1$，所以改写前后允许的 $x$ 完全相同。例如，$x=0.8$ 可以配上 $t=0.2$；$x=1.2$ 则要求 $t=-0.2$，仍然不符合条件。

用 $g$ 表示，新增的等式就是 $g(x)+t=0$。这里的 $g$ 仍是原来的约束函数，$t$ 是额外引入的变量。对一般不等式 $g_i(x)\le 0$，也可以这样写：

$$g_i(x)+t_i=0,t_i\ge 0.$$

$t_i$ 统一称为约束余量。只有在本例这样的特定表达式中，它才恰好等于几何距离。

接着，用第 5 节的对数障碍处理 $t\ge 0$。固定 $\mu >0$ 后，得到：

$$\underset{x,t}{\min }(x-2{)}^2-\mu \log t,x-1+t=0,t>0.$$

等式成立时，把 $t=1-x$ 代回目标，就恢复了第 5 节的 $(x-2{)}^2-\mu \log (1-x)$。保留 $x,t$ 两个变量后，模型通过等式规定两者的关联；求解器在迭代过程中可以暂时允许等式有误差，再逐步修正。

这种处理允许算法从违反原约束的初值开始迭代。 假设初值为 $x=1.2$，原目标和约束函数都能够计算：$f(1.2)=0.64$、$g(1.2)=0.2$。直接使用障碍项时，需要计算 $\log (1-1.2)=\log (-0.2)$，问题出在这个对数的定义域。

引入 $t$ 后，求解器使用的目标变为 $(x-2{)}^2-\mu \log t$。保留 $x=1.2$，给新增变量选一个正初值，例如 $t=0.1$，就可以计算这个目标。这里的 0.1 是演示用的初值；只考虑对数的定义域，0.2 或 1 也可以。 初值的数值尺度会影响后续迭代。这个 $t$ 初值无需与 $\mu$ 相等：前者是待修正的变量，后者是障碍参数。当前关联等式的误差为：

$$x-1+t=1.2-1+0.1=\mathrm{0.3.}$$

求解器可以一边保持 $t>0$，一边调整 $x,t$，让这项误差逐步接近零。当 $x-1+t=0$ 成立时，$x=1-t<1$，原不等式也随之满足。这里，初始的 $t=0.1$ 是待修正的变量值；等式满足后，它才准确表示当前 $x$ 到边界的余量。

对一般约束 $g_i(x)\le 0$，这个计算安排也相同：先选正的 $t_i$，使 $\log t_i$ 可以计算，再修正 $g_i(x)+t_i=0$。复杂的约束关系保留在等式中，对数的定义域要求统一变成简单的 $t_i>0$。最终解仍然需要满足约束关系；增加 $t_i$ 保留了原问题允许的 $x$ 范围。

IPOPT 在初始化阶段就按确定的规则安排松弛变量。 它参考当前约束值，并根据内部的边界调整规则把松弛变量放到允许范围内部。之后通过步长控制保持其合法性。讲义中的 $t=0.1$ 用来说明这种安排，实际数值由求解器的规则与配置决定。这套过程在正式迭代前就已安排，初始化也只保证相关表达式可计算，问题能否达到可行解、能否收敛仍需继续求解。

实际建模时，应把原不等式及其边界交给 IPOPT，由它组织内部的障碍处理。如果自行把 $\log (1-x)$ 写进目标，求解器会按这个表达式求值；在 $x=1.2$ 处仍可能发生求值错误。模型中 $f,g$ 自身的定义域也需要满足。讲义用 $g(x)+t=0$、$t>0$ 统一说明数学关系，求解器内部可以使用不同的松弛变量符号和边界表示。

8.2 用乘子表达约束的作用，建立联合求解条件

有了余量变量，我们能够检查约束关系是否满足。接下来还需要判断目标是否已经达到约束允许的最优位置。第 7 节用乘子表达梯度平衡，现在把同样的方法用于固定 $\mu$ 的障碍子问题。

这里会同时求原始变量 $x,t$ 与相应的乘子。$t$ 的加入没有改变乘子的含义：它们仍然表达约束梯度怎样参与最优性条件。

现在需要处理等式 $x-1+t=0$。沿用第 7 节的方法，为它引入一个乘子，暂时记为 $y$。等式本身不规定乘子的正负。障碍子问题的拉格朗日函数为：

$${\mathcal{L}}_{\mu }(x,t,y)=(x-2{)}^2-\mu \log t+y(x-1+t).$$

此时 $\mu$ 是已给定的数，$x,t,y$ 是需要一起求出的未知量。我们用两条梯度平衡方程，加上关联等式，确定它们需要满足的关系。

先对 $x$ 求偏导。 在多变量函数中，对 $x$ 求偏导时把 $t,y$ 固定：$(x-2{)}^2$ 的导数为 $2(x-2)$，$-\mu \log t$ 的导数为零，$y(x-1+t)$ 的导数为 $y$。因此：

$$\frac{\partial {\mathcal{L}}_{\mu }}{\partial x}=2(x-2)+y=0.$$

这里保留了 $x,t$ 两个独立的变量位置；两者需要满足的关系会由 $x-1+t=0$ 同时保证。若先代入 $t=1-x$ 消去 $t$，就回到了第 5 节直接对单变量障碍目标求导的做法。

再对 $t$ 求偏导。 这次把 $x,y$ 固定：$(x-2{)}^2$ 的导数为零，$-\mu \log t$ 的导数为 $-\mu /t$，$y(x-1+t)$ 的导数为 $y$。因此：

$$\frac{\partial {\mathcal{L}}_{\mu }}{\partial t}=-\frac{\mu }{t}+y=0⟹ty=\mu .$$

因为 $\mu >0$、$t>0$，满足这条方程时必有 $y=\mu /t>0$。$y$ 的正号由关于 $t$ 的方程推导出来；引入这个等式乘子时，我们没有预设它的符号。

最后，把原有的关联等式也放进来：

$$\begin{array}{rl}2(x-2)+y& =0,\\ x-1+t& =0,\\ ty& =\mu .\end{array}$$

这三条方程分别要求：目标与约束的梯度平衡，位置与余量对应，余量与乘子的乘积等于当前的 $\mu$。到这里，障碍子问题已经变成了一组需要同时满足的方程。

这一乘积关系适用于所有采用标准对数障碍的问题：固定 $\mu >0$、满足一阶最优性条件时，每条不等式的余量与对应乘子的乘积都等于 $\mu$。该关系来自对数障碍项对余量变量的求导，与目标函数和约束函数的具体形式无关；迭代中的当前值可以尚未满足它。

8.3 对照原问题，明确 y 与 λ 的对应关系

上一步引入了 $y$，第 7 节使用的则是 $\lambda$。要理解两者的关系，可以比较它们在关于 $x$ 的方程中分别乘在哪一项前面。

第 7 节原问题的拉格朗日函数为：

$$L(x,\lambda )=f(x)+\lambda g(x).$$

对 $x$ 求偏导，得到：

$$f^{\prime }(x)+\lambda g^{\prime }(x)=0.$$

其中，$\lambda$ 是原约束梯度 $g^{\prime }(x)$ 前面的系数。现在把障碍子问题的拉格朗日函数展开：

$$\begin{array}{rl}{\mathcal{L}}_{\mu }(x,t,y)& =f(x)-\mu \log t+y[g(x)+t]\\ & =f(x)+yg(x)+yt-\mu \log t.\end{array}$$

对 $x$ 求偏导时，$t,y$ 固定，最后两项 $yt-\mu \log t$ 的导数为零，所以得到：

$$f^{\prime }(x)+y{g}^{\prime }(x)=0.$$

把两个结果放在一起：

$$\begin{array}{rl}{f}^{\prime }(x)+\lambda g^{\prime }(x)& =0,\\ f^{\prime }(x)+y{g}^{\prime }(x)& =0.\end{array}$$

$\lambda$ 和 $y$ 都是同一个约束梯度 $g^{\prime }(x)$ 前面的系数。 这就是两者在关于 $x$ 的方程中承担相同作用的具体含义。在本例中，$g^{\prime }(x)=1$，两行分别就是 $2(x-2)+\lambda =0$ 和 $2(x-2)+y=0$。

除了梯度平衡，还要比较它们与余量的关系。原问题的互补条件为 $\lambda g(x)=0$；代入 $g(x)=-t$，就得到 $t\lambda =0$。障碍子问题给出的则是 $ty=\mu$，且 $y>0$。随着 $\mu$ 减小，乘积要求逐步接近原来的零互补条件。

因此，后文把障碍子问题中的 $y$ 也记为 $\lambda$，保持它与原不等式的对应关系。这是对同一个未知量换用后文的统一名称。有限的 $\mu$ 下，障碍子问题的解和乘子数值可以与原问题的最优解和乘子不同。

采用这个记号后，三条方程写成：

$$\begin{array}{rl}2(x-2)+\lambda & =0,\\ x-1+t& =0,\\ t\lambda & =\mu ,\end{array}t>0,\lambda >0.$$

第三条称为扰动互补条件：“扰动”指原来的乘积要求 $t\lambda =0$ 被改为 $t\lambda =\mu$。它由上一步对 $t$ 求偏导得到。

第 7.5 节的 $\nu$ 继续用于原问题自带的等式 $h(x)=0$。$\nu$ 与 $\lambda$ 这两个名字按约束来源区分；这里写成等式的 $g(x)+t=0$，仍对应原来的那条不等式。

8.4 代入当前数值，计算每条方程的误差

方程已经列出，接下来需要求出满足它们的数值。求解通常从一组初值开始，因此首先要检查：当前数值代进去以后，每条方程还差多少。

对一条方程，计算“左边减右边”，这个差值就叫残差。例如，$t\lambda =\mu$ 的残差是 $t\lambda -\mu$。残差为零时，这条方程成立；残差的绝对值反映当前左右两边相差多少。

给三条方程的残差分别命名：

要满足的方程	残差的名称与表达式	检查的内容
$2(x-2)+\lambda =0$	$r_{\mathrm{dual}}=2(x-2)+\lambda$	梯度是否平衡
$x-1+t=0$	$r_{\mathrm{ineq}}=x-1+t$	余量是否与位置对应
$t\lambda =\mu$	$r_{\mathrm{cent}}=t\lambda -\mu$	乘积是否等于当前的 $\mu$

这些 $r$ 都由当前的 $x,t,\lambda ,\mu$ 计算得到。此处先按右侧的具体含义理解它们，后面再说明下标的来源。

例如，固定 $\mu =0.1$，当前取 $x=0.8$、$t=0.2$、$\lambda =0.4$。代入得到：

$$\begin{array}{rl}{r}_{\mathrm{dual}}& =2(0.8-2)+0.4=-2,\\ r_{\mathrm{ineq}}& =0.8-1+0.2=0,\\ r_{\mathrm{cent}}& =0.2\times 0.4-0.1=-\mathrm{0.02.}\end{array}$$

第二行已经为零，说明余量与位置对应。第一行仍为 $-2$，说明梯度还没有平衡。第三行的 $-0.02$ 表示：当前乘积为 $0.08$，比要求的 $0.1$ 小了 $0.02$。

如果只看第三行，把 $\lambda$ 从 $0.4$ 改为 $0.5$ 就能使它归零。但此时第一行等于 $-1.9$，仍未满足梯度平衡。这说明，只修正一条方程还不够；需要联合调整 $x,t,\lambda$，让三条方程同时成立。

本例可以直接算出三条方程同时成立时的结果。第二行给出 $x=1-t$，代入第一行得到 $\lambda =2(1+t)$，再代入第三行：

$$2t(1+t)=\mathrm{0.1.}$$

取正根并依次求出 $x,\lambda$，得到：

$$t\approx 0.047722558,x\approx 0.952277442,\lambda \approx \mathrm{2.095445115.}$$

这一组数值使三项残差都接近零，其中的 $x$ 与第 5.2 节直接求障碍目标最小值得到的结果一致。这里 $\mu =0.1$ 对应的解仍留有正余量；原问题的最优点 $x=1$ 则位于边界上。

8.5 多条不等式各自对应一项乘积残差

前面只有一条不等式，因此 $r_{\mathrm{cent}}$ 是一个数。如果有 $m$ 条不等式，第 $i$ 条有自己的余量 $t_i$ 和乘子 ${\lambda }_i$，各自需要满足 $t_i{\lambda }_i=\mu$。

把每条方程的残差按顺序放在一起，就得到向量：

$$r_{\mathrm{cent}}=\left[\begin{array}{c}{t}_1{\lambda }_1-\mu \\ t_2{\lambda }_2-\mu \\ ⋮\\ t_m{\lambda }_m-\mu \end{array}\right].$$

这个向量有 $m$ 个分量，每个分量计算同一条约束的“余量 × 乘子 − μ”。 例如，两条约束当前有 $t=(0.2,0.1{)}^T$、$\lambda =(0.4,1.5{)}^T$，且 $\mu =0.1$，则：

$$r_{\mathrm{cent}}=\left[\begin{array}{c}0.2\times 0.4-0.1\\ 0.1\times 1.5-0.1\end{array}\right]=\left[\begin{array}{c}-0.02\\ 0.05\end{array}\right].$$

第一个乘积比要求小 $0.02$，第二个乘积比要求大 $0.05$。向量保留了每条约束各自的误差。

为了简写，定义 $T=\mathrm{diag}(t_1,\dots ,t_m)$：对角线上放各个 $t_i$，其余位置为零。再用 $1$ 表示长度为 $m$ 的全 1 向量。于是：

$$T\lambda =\left[\begin{array}{c}{t}_1{\lambda }_1\\ ⋮\\ t_m{\lambda }_m\end{array}\right],\mu 1=\left[\begin{array}{c}\mu \\ ⋮\\ \mu \end{array}\right].$$

两者相减，得到同一个残差向量的紧凑写法：

$$r_{\mathrm{cent}}=T\lambda -\mu 1.$$

下标 cent 来自 centrality，表示中心性。第 5.2 节把不同 $\mu$ 下的障碍子问题解连成中心路径；路径上的解满足梯度平衡、约束关系和 $t_i{\lambda }_i=\mu$。$r_{\mathrm{cent}}$ 专门检查最后一项。它为零时，仍需检查其余方程，才能判断当前值是否满足整组条件。

8.6 加入原始等式，得到一般问题的四组残差

现在把一维例子推广到一般问题：$x$ 可以包含多个变量，约束中既有不等式 $g(x)\le 0$，也有原始等式 $h(x)=0$。不等式继续引入余量 $t$ 和乘子 $\lambda$，原始等式的乘子记为 $\nu$。障碍子问题的拉格朗日函数为：

$${\mathcal{L}}_{\mu }(x,t,\nu ,\lambda )=f(x)-\mu \sum _i\log t_i+{\nu }^{T}h(x)+{\lambda }^T[g(x)+t].$$

前面关于 $x$ 的方程只包含一个导数。现在需要对 $x$ 的各个分量求偏导，并把所有约束的贡献加起来。用 $J_h,J_g$ 表示约束的 Jacobian，每一行是一条约束对各个变量的导数，则：

$$\begin{array}{rl}{J}_h(x{)}^T\nu & =\sum _j{\nu }_j\nabla h_j(x),\\ J_g(x{)}^T\lambda & =\sum _i{\lambda }_i\nabla g_i(x).\end{array}$$

这些矩阵乘法表示：每条约束的梯度乘上自己的乘子，再求和。因此，梯度平衡残差也沿用前面的结构；原始等式另外增加一组残差 $h(x)$。四组残差合起来为：

$$\begin{array}{rl}{r}_{\mathrm{dual}}& =\nabla f(x)+J_h(x{)}^T\nu +J_g(x{)}^T\lambda ,\\ r_{\mathrm{eq}}& =h(x),\\ r_{\mathrm{ineq}}& =g(x)+t,\\ r_{\mathrm{cent}}& =T\lambda -\mu 1.\end{array}$$

它们依次检查：梯度平衡、原始等式、余量与约束值的关联、余量与乘子的乘积。与一维例子相比，新增了原始等式这一组，其余三组的计算含义保持相同。

这里可以补全下标的含义：eq 对应原始等式；ineq 对应原不等式的来源，具体计算的是改写后的等式误差 $g(x)+t$；dual 是梯度平衡残差的惯用名称。对应第 7.6 节的术语，扩展模型中的 $x,t$ 是原始变量，乘子 $\nu ,\lambda$ 是对偶变量。

8.7 根据残差修正变量，再逐步调整 μ

回到第 8.4 节的数值例子：算法从一组 $x,t,\lambda$ 出发，算出了 $(-2,0,-0.02)$ 这三个残差。它已经知道哪些方程还没满足，接下来需要确定变量应该各自改变多少。

牛顿法可以把这些方程在当前点附近线性化，再求解一个线性方程组，得到联合修正方向：

$$\Delta x,\Delta t,\Delta \lambda .$$

$\Delta x$ 表示原变量的修正方向，另外两项分别表示余量与乘子的修正方向。三个方向由同一组方程共同确定，因此约束关系、梯度平衡和乘积要求都会参与计算 $\Delta x$。求解器选取步长并更新这些量，随后重新计算残差。对一般问题，还需要一起计算原始等式乘子的修正方向 $\Delta \nu$。

保留这些辅助未知量，也便于利用大型问题中的稀疏结构。 如果一条约束只涉及部分变量，那么它对其余变量的偏导为零。每个 $t_i$ 又只出现在自己的关联等式和乘积方程中。把这些关系组织成线性方程组后，可以用稀疏矩阵保存非零项，并利用稀疏线性求解器计算修正方向。

这解释了增加未知量的计算价值：求解器获得了约束与变量之间的明确关联，可以利用这种结构组织存储和计算。代价是方程组变大，实际速度取决于问题规模、稀疏性和实现方式。一维例子适合直接代入求导；大型问题则可以利用联合方程的结构处理多条约束。

IPOPT 采用原始—对偶内点法，把障碍处理与上述联合更新结合起来，并保持松弛变量和对应的不等式乘子为正。第 8.1 节中 $x=1.2,t=0.1$ 的例子说明了这种迭代允许的中间状态：对数可以计算，原约束和关联等式的误差还需要继续修正。

固定正的 $\mu$ 时，四组残差都接近零，表示接近当前障碍子问题的一阶最优性条件。当前子问题要求 $t_i{\lambda }_i=\mu$。

原问题的互补条件写为 ${\lambda }_i{g}_i(x)=0$。引入余量后，关联等式 $g_i(x)+t_i=0$ 成立时，有 $g_i(x)=-t_i$，因此可以等价改写为：

$${\lambda }_i{g}_i(x)=0⟺-{\lambda }_i{t}_i=0⟺t_i{\lambda }_i=0.$$

这里用到了引入松弛变量时规定的关联关系。迭代中若 $g_i(x)+t_i$ 还有误差，$t_i$ 与 $-g_i(x)$ 就尚未完全对应。因此，求解器逐步减小 $\mu$，并继续调整变量和乘子，在减小关联等式误差的同时，让 $t_i{\lambda }_i$ 逐渐接近零。

进入下一节前，可以按作用整理已经出现的量：

量	表达的内容	求解过程中怎样变化
$x$	原模型需要的结果	与辅助变量一起迭代更新
$t$	每条不等式的余量	保持为正，并逐步与约束值对应
$\lambda ,\nu$	约束梯度参与平衡的系数	随最优性方程一起求出
$\mu$	当前障碍子问题的参数	由障碍参数策略调整
$\alpha$	一次迭代采用的步长系数	根据合法范围和迭代进展选择

后面两章分别处理两个计算环节：第 9 节讨论给定 $\mu$ 时怎样求联合修正方向；第 10 节讨论何时更新 $\mu$、选取多大的值。实际求解中，$\mu$ 的更新可以与变量迭代交替进行。

9. 用牛顿法求解约束最优性方程

上一节已经明确了各个未知量的含义，以及它们需要满足的方程。现在可以使用熟悉的牛顿法计算修正量，只是需要一起更新的量增加了。

固定 $\mu$，用 $z$ 收集原始变量与乘子，用 $R(z)$ 收集相应的残差。希望求出的状态满足 $R(z)=0$。在当前点附近做一阶展开：

$$R(z+\Delta z)\approx R(z)+J_R(z)\Delta z.$$

这里的 $z$ 表示当前迭代点，我们希望新位置的真实残差 $R(z+\Delta z)$ 为零。牛顿法先用右侧的线性展开预测这个残差，令预测值为零，求出本轮修正量，再到新位置重新检查真实残差。由此得到需要求解的线性方程：

$$J_R(z)\Delta z=-R(z).$$

这里的 $J_R$ 是残差对全部未知量的 Jacobian。它包含每条方程对变量和乘子的导数。求出 $\Delta z$ 后，还需要选择步长、检查新点，再决定本轮采用多大的修正。下面先用已经熟悉的一维例子写出这个矩阵，再推广到一般问题，最后解释其中的二阶信息。

9.1 从一维例子的三条残差写出联合方程

沿用第 8 节的一维例子。固定 $\mu$，需要同时满足的三条条件及其残差为：

$$\begin{array}{rl}{r}_{\mathrm{dual}}& =2(x-2)+\lambda ,\\ r_{\mathrm{ineq}}& =x-1+t,\\ r_{\mathrm{cent}}& =t\lambda -\mu .\end{array}$$

当前的 $x,t,\lambda$ 已经有数值，右侧残差也能够计算。本轮要决定的是它们的修正量 $\Delta x,\Delta t,\Delta \lambda$。分别查看修正后每条残差会怎样变化。

梯度平衡这一行对 $x$ 的导数为 2，对 $\lambda$ 的导数为 1，因此修正后的值为：

$$2(x+\Delta x-2)+(\lambda +\Delta \lambda )=r_{\mathrm{dual}}+2\Delta x+\Delta \lambda .$$

让它为零，得到 $2\Delta x+\Delta \lambda =-r_{\mathrm{dual}}$。

余量关联这一行同样是线性的。修正后的值为：

$$(x+\Delta x)-1+(t+\Delta t)=r_{\mathrm{ineq}}+\Delta x+\Delta t.$$

让它为零，得到 $\Delta x+\Delta t=-r_{\mathrm{ineq}}$。

乘积条件这一行包含两个未知量的乘积。展开后得到：

$$\begin{array}{rl}(t+\Delta t)(\lambda +\Delta \lambda )-\mu & =r_{\mathrm{cent}}+t\Delta \lambda +\lambda \Delta t\\ & +\Delta t\Delta \lambda .\end{array}$$

牛顿法的一阶近似省略最后的二阶项 $\Delta t\Delta \lambda$。令线性预测值为零，得到 $t\Delta \lambda +\lambda \Delta t=-r_{\mathrm{cent}}$。这也说明，走到新位置后仍需要重新检查真实残差。

现在把三行按同一个变量顺序排列。列依次对应 $\Delta x,\Delta \lambda ,\Delta t$，便得到联合线性方程组：

$$\left[\begin{array}{ccc}2& 1& 0\\ 1& 0& 1\\ 0& t& \lambda \end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \lambda \\ \Delta t\end{array}\right]=-\left[\begin{array}{c}{r}_{\mathrm{dual}}\\ r_{\mathrm{ineq}}\\ r_{\mathrm{cent}}\end{array}\right].$$

矩阵的每一行对应一条条件，每一列对应一个量的修正。 例如，第一行第三列为零，因为梯度平衡式不直接含 $t$；第三行第二列为 $t$，因为乘积 $t\lambda$ 对 $\lambda$ 的偏导为 $t$。

求解这个方程组，就能同时得到三种修正量。左上角的 2 来自目标函数的二阶导数；本例的约束是线性的，其二阶导数为零。一般问题会把这个位置扩展成矩阵，下面先看它在完整方程中的位置。

9.2 多变量、多约束的联合矩阵

上一节已经把三条标量方程排成矩阵。多个变量与约束沿用同样的过程：先写出每条修正方程，再按统一的未知量顺序收集系数。 其中一些系数会组成一个矩阵块。下面先把这些块展开成具体数字，再看一般写法。

本节所有系数和残差都在当前迭代点计算。在一次线性求解中，它们是已知量，要解的是各个 $\Delta$ 修正量。更新到新点后，再重新计算系数和残差。

先用两个变量、一条等式写出完整矩阵。 复用第 3 节的例子，把变量记为 $x_1,x_2$：

$$\underset{x_1,x_2}{\min }(x_1-2{)}^2+x_2^2,h(x)=x_1+x_2-1=0.$$

给等式配上乘子 $\nu$，拉格朗日函数为：

$$L(x_1,x_2,\nu )=(x_1-2{)}^2+x_2^2+\nu (x_1+x_2-1).$$

对 $x_1,x_2$ 分别求偏导，得到两条梯度平衡方程，再加上等式本身。一共有三项残差：

$$\begin{array}{rl}{r}_{\mathrm{dual},1}& =2(x_1-2)+\nu ,\\ r_{\mathrm{dual},2}& =2x_2+\nu ,\\ r_{\mathrm{eq}}& =x_1+x_2-1.\end{array}$$

现在分别修正 $x_1,x_2,\nu$。第一条残差随 $x_1$ 改变的系数为 2，随 $\nu$ 改变的系数为 1；其余两行同理。让修正后的线性预测值为零，就得到：

$$\begin{array}{rl}2\Delta x_1+\Delta \nu & =-r_{\mathrm{dual},1},\\ 2\Delta x_2+\Delta \nu & =-r_{\mathrm{dual},2},\\ \Delta x_1+\Delta x_2& =-r_{\mathrm{eq}}.\end{array}$$

列顺序固定为 $\Delta x_1,\Delta x_2,\Delta \nu$。第一行的系数依次是 $2,0,1$，第二行是 $0,2,1$，第三行是 $1,1,0$，于是：

$$\left[\begin{array}{ccc}2& 0& 1\\ 0& 2& 1\\ 1& 1& 0\end{array}\right]\left[\begin{array}{c}\Delta x_1\\ \Delta x_2\\ \Delta \nu \end{array}\right]=-\left[\begin{array}{c}{r}_{\mathrm{dual},1}\\ r_{\mathrm{dual},2}\\ r_{\mathrm{eq}}\end{array}\right].$$

例如，从 $x_1=x_2=0,\nu =0$ 出发，右侧为 $(4,0,1{)}^T$。求解得到 $\Delta x_1=1.5$、$\Delta x_2=-0.5$、$\Delta \nu =1$。本例的残差展开没有被省略的高阶项，因此采用这组完整修正后，三条方程恰好同时成立。

把同类变量放在一起，就得到分块写法。 将 $\Delta x_1,\Delta x_2$ 收集成向量 $\Delta x$，把前两项残差收集成 $r_{\mathrm{dual}}$。上面矩阵中的系数可以分成四块：

$$W=\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right],A=J_h(x)=\left[\begin{array}{cc}1& 1\end{array}\right],A^T=\left[\begin{array}{c}1\\ 1\end{array}\right].$$

因此，同一个 $3\times 3$ 方程组可以写为：

$$\left[\begin{array}{cc}W& A^T\\ A& 0\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \nu \end{array}\right]=-\left[\begin{array}{c}{r}_{\mathrm{dual}}\\ r_{\mathrm{eq}}\end{array}\right].$$

这里第一块行包含两条梯度平衡方程，第二块行包含一条等式方程。$W$ 是 $2\times 2$ 的矩阵，$A^T$ 是 $2\times 1$ 的列向量，$A$ 是 $1\times 2$ 的行向量，右下角的零在本例中只有一个元素。

$A$ 与 $A^T$ 的位置可以直接从前面的三条方程核对。变量变化使等式误差改变 $\Delta x_1+\Delta x_2=A\Delta x$。乘子变化则分别给两条梯度平衡式增加一个 $\Delta \nu$，合起来就是 $A^T\Delta \nu$。等式 $h(x)$ 本身不依赖 $\nu$，所以右下角为零。

一般问题也具有这个结构。$A=J_h(x)$ 的每一行是一条等式对各变量的导数，因此 $A\Delta x$ 计算约束值的线性变化；梯度平衡中各个等式乘子的贡献为 $A^T\nu$，因此改变乘子时会出现 $A^T\Delta \nu$。$W$ 则记录固定当前乘子时，梯度平衡对 $x$ 的导数，具体组成在第 9.3 节展开。

加入不等式后，先写出四组修正方程。 现在原问题同时包含 $h(x)=0$ 和 $g(x)\le 0$，需要一起修正 $x,\nu ,\lambda ,t$。沿用第 8 节的残差，并记 $A=J_h(x)$、$G=J_g(x)$。

第一组是梯度平衡。改变 $x$ 产生的线性变化记为 $W\Delta x$；改变两类乘子，分别产生 $A^T\Delta \nu$ 和 $G^T\Delta \lambda$。所以这一组方程为：

$$W\Delta x+A^T\Delta \nu +G^T\Delta \lambda =-r_{\mathrm{dual}}.$$

这里 $W$ 已经包含约束梯度随 $x$ 变化的影响。$G^T$ 的来源与前面 $A^T$ 相同：每条不等式的乘子乘在自己的约束梯度前面。

第二组是原始等式。它只依赖 $x$，因此：

$$A\Delta x=-r_{\mathrm{eq}}.$$

第三组是余量关联等式 $g(x)+t=0$。改变 $x$ 给出 $G\Delta x$，改变 $t$ 则直接增加 $\Delta t$，因此：

$$G\Delta x+\Delta t=-r_{\mathrm{ineq}}.$$

第四组是乘积条件。第 $i$ 条的线性变化为 $t_i\Delta {\lambda }_i+{\lambda }_i\Delta t_i$。定义 $T=\mathrm{diag}(t)$、$\Lambda =\mathrm{diag}(\lambda )$，把这些逐项乘积排列起来，就得到：

$$T\Delta \lambda +\Lambda \Delta t=-r_{\mathrm{cent}}.$$

四组方程已经分别写出，接下来只需收集系数。把列顺序统一为 $\Delta x,\Delta \nu ,\Delta \lambda ,\Delta t$：

对应的方程	$\Delta x$ 前的系数	$\Delta \nu$ 前的系数	$\Delta \lambda$ 前的系数	$\Delta t$ 前的系数
梯度平衡	$W$	$A^T$	$G^T$	$0$
原始等式	$A$	$0$	$0$	$0$
余量关联	$G$	$0$	$0$	$I$
乘积条件	$0$	$0$	$T$	$\Lambda$

表中的零块表示该行没有直接依赖相应的修正量。第三行最后一格为 $I$，因为单位矩阵满足 $I\Delta t=\Delta t$。第四行只使用 $\Delta \lambda ,\Delta t$，因为我们把 $t$ 作为独立变量，各项 $t_i{\lambda }_i$ 对 $x$ 的直接偏导为零；$x$ 与 $t$ 的联系由第三行负责。

将系数表与各行右侧的负残差放在一起，得到完整矩阵形式：

$$\left[\begin{array}{cccc}W& A^T& G^T& 0\\ A& 0& 0& 0\\ G& 0& 0& I\\ 0& 0& T& \Lambda \end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta \nu \\ \Delta \lambda \\ \Delta t\end{array}\right]=-\left[\begin{array}{c}{r}_{\mathrm{dual}}\\ r_{\mathrm{eq}}\\ r_{\mathrm{ineq}}\\ r_{\mathrm{cent}}\end{array}\right].$$

读任意一行时，把这一行的系数块分别乘上对应的修正量，再相加，就恢复了刚才逐组写出的方程。例如，第三行就是 $G\Delta x+I\Delta t=-r_{\mathrm{ineq}}$。

图中的四行表示四组方程，每一组可以包含多条。若 $x$ 有两个分量，另有一条等式和一条不等式，那么未知修正量是 $\Delta x_1,\Delta x_2,\Delta \nu ,\Delta \lambda ,\Delta t$，完整矩阵实际为 $5\times 5$。第一组包含两条梯度平衡方程，另外三组各有一条。

求解这个线性方程组，就能得到各组未知量的联合修正。下一节展开左上角的 $W$，说明目标和非线性约束的二阶导数怎样进入这一块。

9.3 联合矩阵中的拉格朗日 Hessian

第 9.2 节中的梯度平衡残差为：

$$r_{\mathrm{dual}}={\nabla }_{x}L=\nabla f+\sum _j{\nu }_j\nabla h_j+\sum _i{\lambda }_i\nabla g_i.$$

对它关于 $x$ 求导，就得到联合矩阵左上角的 $W$：

$$W={\nabla }_{xx}^{2}L={\nabla }^{2}f+\sum _j{\nu }_j{\nabla }^2{h}_j+\sum _i{\lambda }_i{\nabla }^2{g}_i.$$

这就是拉格朗日 Hessian。它回答的是：在当前乘子取值下，改变 $x$ 会怎样改变梯度平衡。因为平衡式包含目标梯度和约束梯度，对它求导就会同时出现目标与约束的二阶项。

计算这一块时，求的是关于 $x$ 的偏导，所以乘子保持不变。乘子自身也会更新，它对梯度平衡的影响由 $A^T\Delta \nu$、$G^T\Delta \lambda$ 表达，对应上一节完整矩阵第一行的另外两个非零块。

在一维示例中，$g=x-1$ 是线性的，$g^{\prime \prime }=0$，所以 $W=f^{\prime \prime }=2$。换成弯曲边界，约束就会产生自己的二阶项。

用一个圆形边界，可以具体算出非线性约束新增的二阶项。考虑在单位圆盘内寻找最接近 $(2,0)$ 的点：

$$\underset{x,y}{\min }f(x,y)=(x-2{)}^2+y^2,g(x,y)=x^2+y^2-1\le 0.$$

这里的 $x,y$ 是两个坐标，下面的梯度和 Hessian 都按 $(x,y)$ 的顺序排列。目标 $f$ 是当前位置到 $(2,0)$ 的距离平方，约束要求位置在单位圆盘内。我们分四步算出这个例子的 $W$。

第一步：确定原问题的最优点。 圆盘约束 $x^2+y^2\le 1$ 给出 $x\le 1$，因此 $2-x\ge 1$。对任意可行点，都有：

$$f(x,y)=(2-x{)}^2+y^2\ge 1.$$

点 $(1,0)$ 满足圆盘约束，且 $f(1,0)=1$，恰好达到这个下界。所以原问题的最优点是 $(1,0)$。这里先手算出最优点，用它检查后面的乘子和矩阵。

第二步：由梯度平衡求出乘子。 先分别对目标和约束求一阶偏导：

$$\nabla f(x,y)=\left[\begin{array}{c}2(x-2)\\ 2y\end{array}\right],\nabla g(x,y)=\left[\begin{array}{c}2x\\ 2y\end{array}\right].$$

本例只有一条不等式，因此只有一个对应的乘子 $\lambda$。把 $(x,y)=(1,0)$ 代入梯度平衡条件 $\nabla f+\lambda \nabla g=0$：

$$\left[\begin{array}{c}-2\\ 0\end{array}\right]+\lambda \left[\begin{array}{c}2\\ 0\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$

第一行要求 $-2+2\lambda =0$，所以 ${\lambda }^{\ast }=1$；第二行已经是 $0=0$。再检查其余条件：$g(1,0)=0$，${\lambda }^{\ast }=1\ge 0$，并且 ${\lambda }^{\ast }g(1,0)=0$。因此，这组点与乘子满足原问题的 KKT 条件。

第三步：继续求导，得到两个 Hessian。 对目标梯度的第一项 $2(x-2)$ 求导：关于 $x$ 的导数为 2，关于 $y$ 的导数为 0。对第二项 $2y$ 求导：关于 $x$ 的导数为 0，关于 $y$ 的导数为 2。按行排列，得到：

$${\nabla }^{2}f=\left[\begin{array}{cc}\frac{{\partial }^{2}f}{\partial x^2}& \frac{{\partial }^{2}f}{\partial x\partial y}\\ \frac{{\partial }^{2}f}{\partial y\partial x}& \frac{{\partial }^{2}f}{\partial y^2}\end{array}\right]=\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]=2I.$$

约束梯度的两项是 $2x$ 和 $2y$，继续求导同样得到：

$${\nabla }^{2}g=\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]=2I.$$

这里的 $I$ 是 $2\times 2$ 单位矩阵。目标梯度中的 $2(x-2)$ 与约束梯度中的 $2x$ 相差一个常数，对 $x$ 再求一次导数都得到 2，因此两个 Hessian 恰好相同。

第四步：按乘子加权，组成拉格朗日 Hessian。 本例的拉格朗日函数是 $L=f+\lambda g$。计算关于 $(x,y)$ 的二阶偏导时固定 $\lambda$，于是：

$$W={\nabla }^{2}f+\lambda {\nabla }^{2}g=2I+\lambda 2I=2(1+\lambda )I.$$

在刚才求出的原问题最优点处，乘子为 ${\lambda }^{\ast }=1$，所以：

$$W{|}_{\lambda =1}=2I+2I=4I=\left[\begin{array}{cc}4& 0\\ 0& 4\end{array}\right].$$

也可以直接展开来核对。固定 $\lambda =1$ 时：

$$\begin{array}{rl}L(x,y,1)& =(x-2{)}^2+y^2+(x^2+y^2-1)\\ & =2x^2+2y^2-4x+3.\end{array}$$

这个表达式关于 $x,y$ 的二阶偏导分别为 4，混合偏导为 0，因此同样得到 $W=4I$。

这里目标贡献了 $2I$，非线性约束在乘子 $\lambda =1$ 时又贡献了 $2I$。如果只提供目标函数的 Hessian，就会缺少后面这一部分。迭代过程中，求解器按当前乘子计算 $W=2(1+\lambda )I$；$4I$ 是本例在 $\lambda =1$ 时的具体值。

9.4 消去余量和乘子的修正量，联系障碍曲率

回看一维例子，拉格朗日 Hessian 是 $W=2$，而第 5 节算出的障碍目标二阶导数是 $2+\mu /(1-x{)}^2$。两种写法记录曲率的位置不同：直接消去余量后，障碍曲率包含在目标的二阶导里；保留余量和乘子后，这部分信息还分布在乘积方程及其导数中。本节通过消元把这两种写法联系起来。

先看第 9.2 节完整矩阵的第三行：

$$G\Delta x+\Delta t=-r_{\mathrm{ineq}}⟹\Delta t=-r_{\mathrm{ineq}}-G\Delta x.$$

它把 $\Delta t$ 表达为 $\Delta x$ 与当前残差的组合。代入第四行 $T\Delta \lambda +\Lambda \Delta t=-r_{\mathrm{cent}}$，得到：

$$\Delta \lambda =T^{-1}\Lambda G\Delta x+T^{-1}(\Lambda r_{\mathrm{ineq}}-r_{\mathrm{cent}}).$$

这里各个 $t_i>0$，所以对角矩阵 $T$ 可以求逆。再把这个 $\Delta \lambda$ 代入第一行，收集所有乘在 $\Delta x$ 前面的项，便得到：

$$W+G^T{T}^{-1}\Lambda G.$$

其他由当前残差组成的项可以移到方程右侧。这个结果说明，消去余量与乘子的修正量后，原变量方向上的系数同时包含 $W$，以及由约束 Jacobian、余量和乘子共同构成的附加项。

接下来要把消元结果与障碍目标的二阶导数对应起来，需要先回顾第 5.2 节的中心路径。固定一个 $\mu >0$，求出对应障碍子问题的理想解；让 $\mu$ 连续变化，再把这些解连起来，就得到本讲义例子中的中心路径。它描述不同障碍参数对应的解的位置。

例如，第 5.2 节的一维问题有以下结果：

固定的 $\mu$	对应子问题的最优点 $x_{\mu }$
$0.1$	约 $0.952277$
$0.01$	约 $0.995025$
$0.001$	约 $0.999500$

每个 $x_{\mu }$ 还对应自己的余量 $t_{\mu }$ 和乘子 ${\lambda }_{\mu }$，所以这里也可以把中心路径看成随 $\mu$ 变化的一系列 $(x_{\mu },t_{\mu },{\lambda }_{\mu })$。这些理想解满足相应子问题的完整一阶条件；实际迭代点通常还带有求解误差，可以偏离这条路径。

下面在其中一个固定 $\mu$ 的解处比较两种曲率写法。中心路径上的解同时满足余量关联和扰动互补条件：

$$g(x)+t=0,t\lambda =\mu .$$

第一条给出 $t=-g(x)$，第二条给出 $\lambda =\mu /t$。把两者连起来，就得到后面要用的关系：

$$\lambda =\frac{\mu }{t}=-\frac{\mu }{g(x)}.$$

因为 $g(x)<0$，所以右侧仍然是正数，符合当前障碍子问题对乘子的要求。

在一维例子中，$G=1$、$W=2$，消元后的系数为 $2+\lambda /t$。此时 $g(x)=x-1$，余量为 $t=1-x$。代入刚才的乘子关系，得到：

$$2+\frac{\lambda }{t}=2+\frac{\mu }{t^2}=2+\frac{\mu }{(1-x{)}^2}.$$

这就恢复了第 5 节直接对障碍目标求出的二阶导数。类似地，对于只有目标 $f$ 和一条不等式 $g(x)\le 0$ 的问题，在 $g(x)<0$ 时直接对障碍目标求二阶导数，得到：

$${\nabla }^2\left[f-\mu \log (-g)\right]={\nabla }^{2}f-\frac{\mu }{g}{\nabla }^{2}g+\frac{\mu }{g^2}\nabla g\nabla g^T.$$

这里的前两项可以与 $W$ 直接对应。本段只包含目标 $f$ 和一条不等式 $g$，拉格朗日函数为 $L=f+\lambda g$。按照固定当前乘子、对原变量求二阶偏导的定义：

$$W={\nabla }^{2}f+\lambda {\nabla }^{2}g.$$

在中心路径上的同一个点，代入 $\lambda =-\mu /g$，就有：

$$\begin{array}{rl}{\nabla }^{2}f-\frac{\mu }{g}{\nabla }^{2}g& ={\nabla }^{2}f+\left(-\frac{\mu }{g}\right){\nabla }^{2}g\\ & ={\nabla }^{2}f+\lambda {\nabla }^{2}g\\ & =W.\end{array}$$

这就是“前两项等于 $W$”的具体代入过程。计算 $W$ 时乘子仍作为常数，上面的代入用于比较同一点上的矩阵数值。因此，沿中心路径，障碍目标的完整 Hessian 可以写为：

$${\nabla }^2\left[f-\mu \log (-g)\right]=W+\frac{\mu }{g^2}\nabla g\nabla g^T.$$

其中最后一项对应消去余量和乘子后出现的附加曲率。对于这里的一条不等式，$G=\nabla g^T$，且 $\lambda /t=\mu /t^2=\mu /g^2$，所以它恰好与前面的 $G^T{T}^{-1}\Lambda G$ 对应。

IPOPT 在原始—对偶方程中结合这些信息计算方向：模型的二阶导数接口提供 $W$，求解器结合约束 Jacobian、余量和乘子处理其余部分。第 13 节的 hess_lag 对应其中的 $W$。实际求解还会结合边界处理、缩放、正则化、惯性修正和线搜索，完成本轮更新。

10. 障碍参数的更新策略：mu_strategy

第 8—9 节讨论了固定 $\mu$ 时的求解：计算残差和导数，形成线性方程组，求出变量与乘子的方向，再选择步长并更新状态。第 5.5 节还指出，为了接近原问题的解，需要继续调整 $\mu$。

这就产生了一个独立的算法选择：当前子问题求到什么程度时更新 $\mu$，以及下一次使用多大的 $\mu$。IPOPT 用 mu_strategy 选择这套更新规则。$\mu$ 的变化会调整乘积要求 $t_i{\lambda }_i=\mu$，从而影响接下来的修正方向；每一步实际移动多远仍由步长控制。

选择这个选项，相当于为求解器选择安排障碍子问题的规则。使用者提供原目标和约束，求解器根据规则与迭代表现更新 $\mu$。下面比较两种安排。

10.1 按阶段减小 μ：monotone

monotone 通常维持当前 $\mu$，让当前障碍子问题达到与该 $\mu$ 相配的精度，然后减小 $\mu$，进入下一阶段。概念流程是：

$${\mu }_0⟶{\mu }_1⟶{\mu }_2⟶\cdots ,{\mu }_{j+1}<{\mu }_j.$$

每个阶段内部可以包含多次牛顿迭代。这里的单调性指障碍参数的更新，并不要求原始目标函数在每次迭代都下降。

IPOPT 的相关选项包括 mu_init 和 barrier_tol_factor。后者参与控制当前障碍子问题需要达到的精度。忽略停止阈值等细节，减小规则的核心形式为：

$${\mu }_{\mathrm{new}}=\min \left({\kappa }_{\mu }\mu ,{\mu }^{{\theta }_{\mu }}\right),0<{\kappa }_{\mu }<1,{\theta }_{\mu }>1.$$

10.2 随迭代表现调整 μ：adaptive

我们还可以利用当前迭代的残差与进展，为下一次迭代选择 $\mu$。adaptive 的自由模式采用这种做法；“自由”指参数可以根据当前情况调整，序列可以增大或减小。这些值由算法规则根据当前状态计算。

默认的 mu_oracle="quality-function" 使用质量函数评价候选参数，综合考虑对偶残差、原始残差、互补性与中心性等信息。

进展不足时，IPOPT 会根据 adaptive_mu_globalization 的规则切换到单调模式。默认规则 obj-constr-filter 观察目标与约束违反的进展。这个切换有助于维持收敛行为；全局化机制处理从较远初值出发时的进展，无法保证非凸问题的全局最优解。

10.3 两种策略的迭代特点与比较方法

两种策略都利用障碍子问题逼近原问题，它们安排参数更新的方式不同。比较时，可以先观察参数更新时机，再检查这种安排对求解过程的影响。

问题	monotone	adaptive
何时更新 μ	当前阶段达到相应精度后	自由模式下随迭代表现选择
μ 的变化	阶段间减小	自由模式下可增可减
可能收益	阶段逻辑清晰，轨迹较有规律	某些问题可以减少阶段等待和迭代次数
可能代价	某些阶段可能耗费较多迭代	参数选择与搜索轨迹更复杂，某些问题可能变慢或失败

mu_strategy 会影响优化轨迹，在非凸问题中还可能影响最终到达的局部解。比较它们时，应固定目标、约束、初值、导数和停止条件，记录成功率、迭代数与总时间。

mu_init 主要用于 monotone。使用 adaptive 时，应按照 IPOPT 选项文档理解其独立的参数选择机制。

11. 最优性误差与求解停止条件

前面已经说明了变量、乘子和 $\mu$ 的更新过程。为了判断何时结束，还需要检查当前点距离原问题的最优性条件有多远。第 7 节列出的条件对应三类误差：可行性误差、梯度平衡误差和互补误差。直观上可以观察：

$$\begin{array}{rl}\text{primal error}& :h(x),\max (g(x),0),\\ \text{dual error}& :{\nabla }_{x}L(x,\nu ,\lambda ),\\ \text{complementarity}& :t_i{\lambda }_i.\end{array}$$

约束误差很小，表示当前结果已经接近可接受范围；梯度平衡误差很小，表示目标与约束的局部作用接近平衡；乘积接近零，表示互补关系接近原问题的要求。检查这三类误差，才能同时了解可行性与最优性的进展。

第 8 节固定 $\mu$ 时使用的残差为 $t_i{\lambda }_i-\mu$；这里判断原问题时关心的是 $t_i{\lambda }_i$ 接近零。即使当前障碍子问题已经解得很准，较大的 $\mu$ 仍可能让结果离原问题的最优性条件有距离。

IPOPT 的实际判断还包含变量边界、误差缩放和具体容差。tol 管理整体最优性误差，constr_viol_tol、dual_inf_tol、compl_inf_tol 分别约束相应分量。acceptable_* 配置可以允许在较宽松的条件连续满足时结束。

max_iter、CPU 或墙钟时间限制、外部回调中止也会结束求解。求解器可以在预算耗尽前收敛；预算耗尽本身无法说明已经得到满足要求的解。

IPOPT 日志中常用的字段为：

字段	阅读方式
objective	当前原始目标值
inf_pr	当前原始约束违反的度量；默认显示原问题约束违反
inf_du	内部缩放问题的梯度平衡误差；也称对偶不可行性度量
lg(mu)	障碍参数的十进制对数，例如 −3 对应 0.001
alpha_pr / alpha_du	原始变量 / 乘子采用的步长系数
ls	当前迭代的回溯线搜索次数

如果总时间增加，可以分开检查迭代次数、每次导数求值成本、线性方程求解成本和线搜索次数。仅凭一次求解时间，很难定位是哪一部分变化。

12. 二阶信息的来源：hessian_approximation

到这里，我们已经把求解方向、障碍参数更新和停止条件串了起来。实际实现中还要解决一个计算问题：第 9 节使用的拉格朗日 Hessian 从哪里获得。直接求二阶导可以利用模型曲率，也可能花费较多时间，因此 IPOPT 允许选择二阶信息的来源。

IPOPT 通过非线性规划（NLP）接口获取模型值与导数。选项 hessian_approximation 决定其中的二阶信息怎样获得，提供两种主要方式：

选项值	二阶信息来源	常见考虑
"exact"	使用 NLP 接口提供的拉格朗日 Hessian	模型能够高效计算可靠二阶导时，可以利用其曲率结构
"limited-memory"	根据近期迭代信息构造有限内存拟牛顿近似	二阶导难以获得或成本较高时，可以减少显式 Hessian 求值

在 CasADi 中可以明确写为：

options = {"ipopt.hessian_approximation": "exact"}
# 或者选择另一种二阶信息来源：
options = {"ipopt.hessian_approximation": "limited-memory"}

有限内存方法利用变量与梯度等信息的变化，近似拉格朗日函数的曲率，具体实现属于 L-BFGS 类方法。它仍然需要目标梯度和约束 Jacobian。两种选项使用相同的原目标、约束和最优性条件，变化的是计算修正方向时采用的曲率信息。

exact 模式调用接口给出的矩阵。当应用程序自定义了一个不完整的 Hessian 时，求解器会使用这个近似矩阵，接口名称不会自动补齐缺失项。

这两个选项管理不同的事情：mu_strategy 决定障碍参数如何更新；hessian_approximation 决定方向系统从哪里获得拉格朗日曲率。因此可以交叉组合并做对照实验。

13. 自定义 Hessian 的接口：hess_lag

上一节的 exact 模式需要 NLP 接口提供一个数值矩阵。对于纯符号模型，CasADi 通常可以自动生成相应的导数函数；需要接入专用导数或明确构造近似矩阵时，也可以通过 hess_lag 自定义它。下面把第 9 节的公式对应到这个接口。

13.1 接口的输入与输出

在每次需要二阶信息时，IPOPT 已经有当前的变量与乘子数值。CasADi 的 hess_lag 接口接收这些数值，返回对应的拉格朗日 Hessian 上三角部分，供求解器组装方向方程。

在 CasADi 的统一约束向量写法中：

$$\underset{x}{\min }f(x,p),l_g\le g(x,p)\le u_g.$$

接口要求的矩阵为：

$$H(x,p,\sigma ,\lambda )={\nabla }_{xx}^2\left[\sigma f(x,p)+{\lambda }^{T}g(x,p)\right].$$

这对应第 9 节的拉格朗日曲率。用户提供这个矩阵后，IPOPT 会组装松弛变量与互补条件产生的其他项。

输入	含义
$x$	当前决策变量
$p$	求解期间固定的模型参数；没有参数时长度为 0
$\sigma$，接口名 lam_f	接口传入的目标函数系数
$\lambda$，接口名 lam_g	与约束向量逐项对应的乘子

lam_f 与障碍参数 $\mu$ 各自独立。lam_g 由求解器更新，目标函数中的手工权重已经包含在 $f$ 内。实现这个接口时，我们提供计算规则；每次调用所用的乘子数值由求解器传入。

前文采用 $g\le 0$ 的写法，因此乘子非负。CasADi 支持上下两侧边界，合并后的 lam_g 可以带正负号。读取其数值时，应结合该约束使用的是上界、下界还是等式。

13.2 圆盘问题的符号实现

现在把第 9.3 节已经手算过的圆盘问题写成代码。这样可以同时核对目标、约束、乘子和矩阵。下面的 Python 写法与 C++ 中的 MX::hessian、MX::triu 对应：

import casadi as ca

x = ca.MX.sym("x", 2)
p = ca.MX.sym("p", 0)
sigma = ca.MX.sym("lam_f")
lam = ca.MX.sym("lam_g", 1)

f = (x[0] - 2)**2 + x[1]**2
g = ca.vertcat(ca.dot(x, x) - 1)

L = sigma * f + ca.dot(lam, g)
H = ca.hessian(L, x)[0]
hess_lag = ca.Function(
    "my_hess_lag", [x, p, sigma, lam], [ca.triu(H)]
)

solver = ca.nlpsol("circle", "ipopt", {"x": x, "f": f, "g": g}, {
    "ipopt.mu_strategy": "adaptive",
    "ipopt.hessian_approximation": "exact",
    "hess_lag": hess_lag,
})
result = solver(x0=[0, 0], lbg=[-ca.inf], ubg=[0])

逐步对应数学含义：

1. L 建立 $\sigma f+{\lambda }^{T}g$ 的符号表达式。
2. ca.hessian(L,x)[0] 对决策变量求二阶导。Python 接口同时返回 Hessian 与梯度，[0] 取 Hessian。
3. ca.triu(H) 保留上三角。对称矩阵只需传递一半元素，CasADi 的这个接口采用上三角约定。
4. ca.Function 保存计算规则。每次求值都会使用当时的 $x,p,\sigma ,\lambda$ 生成数值矩阵。
5. hess_lag 属于 CasADi 选项；带 ipopt. 前缀的两个选项传递给 IPOPT。

在传给 nlpsol 的问题字典中，x 是决策变量；sigma 和 lam 则是 hess_lag 的符号输入。IPOPT 在需要矩阵时把当前系数传进来。整个调用关系是：CasADi 保存求导规则，IPOPT 提供当前数值，函数返回矩阵，IPOPT 再用矩阵计算联合修正方向。

13.3 符号约束与外部约束的二阶信息

圆盘例子的目标和约束都由 CasADi 符号表达式给出，因此可以直接求出全部二阶导。实际模型有时还会调用外部几何库，外部接口能够提供的导数阶数决定了这部分曲率能否直接计算。把约束向量按来源分成两部分：

$$g=\left[\begin{array}{c}{g}_{\mathrm{sym}}\\ g_{\mathrm{ext}}\end{array}\right],\lambda =\left[\begin{array}{c}{\lambda }_{\mathrm{sym}}\\ {\lambda }_{\mathrm{ext}}\end{array}\right],$$

完整拉格朗日 Hessian 为：

$$W=\sigma {\nabla }^{2}f+\sum _i{\lambda }_{\mathrm{sym},i}{\nabla }^2{g}_{\mathrm{sym},i}+\sum _j{\lambda }_{\mathrm{ext},j}{\nabla }^2{g}_{\mathrm{ext},j}.$$

如果外部接口只提供约束值和一阶导，可以显式构造以下近似：

$$\tilde{W}=\sigma {\nabla }^{2}f+\sum _i{\lambda }_{\mathrm{sym},i}{\nabla }^2{g}_{\mathrm{sym},i}.$$

此时外部约束仍保留在完整的 $g$ 与 Jacobian 中，求解器继续使用它们判断可行性并形成搜索方向；二阶矩阵中省略了外部约束的曲率项。

原代码里的 Slice(0, m_sym) 选择前一段乘子，dot(...) 形成 ${\lambda }_{\mathrm{sym}}^T{g}_{\mathrm{sym}}$，随后求 Hessian。切片顺序必须与完整约束向量的排列一致。

提供更多二阶信息可能改善局部模型，也会增加导数求值成本。整体收益取决于迭代数、每次求值和线性求解的总成本，需要通过同条件实验判断。

14. 实际求解中的导数检查与结果验证

14.1 模型与初值

明确目标、约束和结果验收条件，检查变量尺度，并保证目标与约束在初值处能够计算。

14.2 导数一致性：derivative_test

IPOPT 用梯度、Jacobian 和 Hessian 计算修正方向。如果函数值与提供的导数不一致，方向计算就可能受到影响。接入外部 Callback 或自定义 hess_lag 后，可以先用 derivative_test 核对这些接口。

一阶与二阶检查各自比较什么。 检查器在正式优化迭代前，用有限差分估计导数，再与接口返回的导数逐项比较。

derivative_test 的取值	检查内容
"none"	关闭检查，默认值
"first-order"	检查目标梯度 $\nabla f$ 和约束 Jacobian $J_g$
"second-order"	先检查一阶导，再检查目标与各约束的 Hessian
"only-second-order"	只检查二阶导，适合一阶导已核对的情况

以一维函数为例，一阶检查利用函数值的变化：

$$f^{\prime }(x)\approx \frac{f(x+h)-f(x)}{h}.$$

二阶检查利用一阶导数的变化：

$$f^{\prime \prime }(x)\approx \frac{f^{\prime }(x+h)-f^{\prime }(x)}{h}.$$

这里 $h$ 是一个小扰动。多变量时，检查器逐个扰动变量，比较对应的梯度、Jacobian 和 Hessian 元素。先核对一阶导，再看二阶导，可以避免错误的一阶导继续影响二阶差分结果。

二阶检查通过拉格朗日 Hessian 接口分别核对目标和各约束的二阶项，因此要检查变量顺序、乘子顺序及目标系数 lam_f。使用 limited-memory 时，通常只启用一阶检查。

在 CasADi 中启用检查。 可以向创建求解器的选项字典中加入：

options = {
    "ipopt.derivative_test": "first-order",
    "ipopt.derivative_test_print_all": "yes",
    "ipopt.derivative_test_perturbation": 1e-8,
    "ipopt.derivative_test_tol": 1e-4,
    "ipopt.point_perturbation_radius": 0.0,
}

一阶检查通过后，将 "first-order" 改为 "second-order"。derivative_test 控制检查过程，hessian_approximation 选择优化迭代的二阶信息来源。

• derivative_test_print_all="yes" 输出全部比较项，便于查看正确项与异常项；默认 "no" 只输出可疑项。
• derivative_test_perturbation 设置相对扰动尺度，默认 $10^{-8}$。扰动太大会增加差分近似误差，太小则容易受到浮点舍入或外部函数噪声影响。
• derivative_test_tol 设置标记差异的阈值，默认 $10^{-4}$。它控制导数检查的判定，与优化停止容差 tol 各自独立。
• point_perturbation_radius=0.0 关闭检查中心点默认允许的随机扰动，便于围绕指定初值复查；有限差分本身的小扰动仍然保留。

阅读输出并定位差异。 grad_f、jac_g、obj_hess、constr_hess 分别对应目标梯度、约束 Jacobian、目标 Hessian 和约束 Hessian。每行给出接口值、有限差分估计值和偏差；行首的 * 标记超过阈值的项。Jacobian 或 Hessian 条目后的 v 表示该位置已经列入接口声明的稀疏结构。

发现差异后，先检查变量与约束的排列、导数符号及系数，再结合扰动尺度复查。最近距离、接触切换等函数可能在几何分支切换处不光滑，可以选择分支稳定的位置检查局部导数。

第 13.3 节主动省略外部约束曲率时，二阶检查会在对应项上报告差异。此时应核对差异是否与省略项一致，同时确认保留的二阶项以及完整的一阶导正确。检查通过说明测试位置的函数与导数一致；模型是否表达了任务要求，仍需要结合任务检查。

导数检查会增加函数调用次数，适合调试阶段使用。检查完成后恢复 derivative_test="none"，再测量正常求解耗时。

14.3 结果与耗时

结合求解状态、约束误差和任务要求判断结果。比较参数时固定模型与初值，一次改变一个因素，同时记录迭代次数与总耗时。