罚函数法，Lagrangian乘数法，递增Lagrangian乘数法*1， 内点法
0 本文的用意和写法
该系列博文面对对有限元解析有所了解，但未接触过接触解析的初学者。本文以实用性为第一目标，试图在概念准确的前提下追求直观描述，特别强调给出某些实用算法的细节，并附有相应程序。由于这些算法细节都经过本人实装测试，其有效性是有一定的保证的。希望读完此系列博文可以自己动手组装自己的解析程序，并希望能够指出其中问题或错误，共同提高。本文的解析对象是泛函的极值问题，这是本文的起点。对泛函的极值问题不了解的读者可以跳过泛函的变分运算步骤，着重于其后的物理概念和算法的讨论。
1 考虑约束条件的最优化算法
泛函的极值问题的求解从数值计算方法的角度来说是最优化算法，而接触问题为其约束条件。罚函数法，Lagrangian乘数法，递增Lagrangian乘数法， 内点法(Interior point method, barrier method或内部罚函数法)都是标准的考虑约束条件的最优化算法。这一节将以图1的问题为例解释这几种解法。
图1 一自由度的接触问题* 2
图1所示问题如没有接触约束，其泛函J为
(0.1) $J=1/2K x^2 -F x$
这里我们将m的位移记为x，外力F<0。该泛函的极值条件为其变分为零，即
(0.2) $\delta J=K x\delta x -F \delta x=0$
基于变分$\delta x$的任意性，要求
(0.3) $K x -F=0$
即得
(0.4) $ x=F/K $
在此我们假设F/K足够大使得$|x|>g_0$，但是接触条件要求物体m不能进入接触面, 即$g=g_0+x>=0$。如g<0, 可以用下述方法导入接触约束。
1.1 罚函数法
采用罚函数法考虑接触约束，其泛函J为(1.1) $J=1/2K x^2 -F x +1/2\varepsilon_N (x+g_0)^2$
这里$\varepsilon_N>0$为罚函数。该泛函的极值条件为其变分为零，即
(1.2) $\delta J=(K+\varepsilon_N) x\delta x -(F-\varepsilon_Ng_0) \delta x=0$
基于变分$\delta x$的任意性，要求
(1.3) $(K+\varepsilon_N) x -(F-\varepsilon_Ng_0)=0$
即得
(1.4) $ x=(F-\varepsilon_Ng_0)/(K+\varepsilon_N) $
当罚函数$\varepsilon_N$足够大时，从该公式可以得到逼近于精确解$x=-g_0$的近似解。
1.2 Lagrangian乘数法
采用Lagrangian乘数法考虑接触约束，其泛函J为
(2.1) $J=1/2K x^2 -F x +\lambda (x+g_0)$
这里$\lambda$为Lagrangian乘数,这是一个待定自变量。泛函J的极值条件为其变分为零，即
(2.2) $\delta J=(Kx+\lambda-F)\delta x +\delta \lambda (x+g_0)=0$
基于变分$\delta x$和$\delta \lambda$的任意性，要求
(2.3) $\begin{eqnarray}
Kx+\lambda=F \\
x=-g_0
\end{eqnarray}$
即得
(2.4) $ x=-g_0; \lambda=F+Kg_0 $
这是该问题的精确解。
(2.1) $J=1/2K x^2 -F x +\lambda (x+g_0)$
这里$\lambda$为Lagrangian乘数,这是一个待定自变量。泛函J的极值条件为其变分为零，即
(2.2) $\delta J=(Kx+\lambda-F)\delta x +\delta \lambda (x+g_0)=0$
基于变分$\delta x$和$\delta \lambda$的任意性，要求
(2.3) $\begin{eqnarray}
Kx+\lambda=F \\
x=-g_0
\end{eqnarray}$
即得
(2.4) $ x=-g_0; \lambda=F+Kg_0 $
这是该问题的精确解。
1.3 递增Lagrangian乘数法
采用递增Lagrangian乘数法考虑接触约束，其泛函J为*3
(3.1) $J=1/2K x^2 -F x +\lambda (x+g_0)+1/2\varepsilon_N (x+g_0)^2$
注意到该泛函既包括Lagrangian乘数$\lambda$也包括罚函数$\varepsilon_N>0$。但与前述Lagrangian乘数法和罚函数法不同的是，在该算法中
* Lagrangian乘数$\lambda$不是自变量
* 罚函数$\varepsilon_N$可以不是很大的数值
-----------------------------------------------------------------
注： 递增Lagrangian乘数法可以看成是罚函数法或Lagrangian乘数法的变种。如从Lagrangian乘数法出发，我们将（2.1）改写为
$J=1/2K x^2 -F x +\lambda (x+g_0)+(\lambda^*-\lambda) (x+g_0)$
这里，我们记$\lambda^*$为Lagrangian乘数的准确值。考虑 $\lambda^*-\lambda$为Lagrangian乘数的增分并将此式与（3.1）比较，可以看出
$\varepsilon_N (x+g_0)\approx \lambda^*-\lambda$
如此Lagrangian乘数的准确值可以如下逼近
$ \lambda^*\approx\lambda+\varepsilon_N (x+g_0)$
对应于下面的（3.5）式。
-----------------------------------------------------------------
泛函J的极值条件为其变分为零，即
(3.2) $\delta J=K x\delta x+ (\lambda+\varepsilon_Nx)\delta x -(F-\varepsilon_Ng_0) \delta x=0$
在此由于Lagrangian乘数$\lambda$不是自变量，我们不需要对其变分。基于变分$\delta x$的任意性，要求
(3.3) $(K+\varepsilon_N)x+\lambda-(F+\varepsilon_Ng_0)=0$
由此可得
(3.4) $ x=(F-\varepsilon_Ng_0-\lambda)/(K+\varepsilon_N) $
但是这里$\lambda$也是待定值，我们需要用迭代方法解此方程。
* 假设Lagrangian乘数的初值$\lambda_0=0$，n=0
* 带入上式求得$x_{n+1}$
* 判断接触条件$g_0+x>=0$是否满足。如不满足，则更新$\lambda$，
(3.5) $\lambda_{n+1}\approx \lambda_n+\varepsilon_N(g_0+x)$
令n->n+1; 并返回上一步重求x。
在上述计算过程中，Lagrangian乘数$\lambda$从零出发逐渐递增，这是我们称之为递增Lagrangian乘数法的原因。同时，x的值从$ x=(F-\varepsilon_Ng_0)/(K+\varepsilon_N) $出发逐渐递减，当$g_0+x$接近于0时计算结束。
可以看出，这里给出的x的初值与上述罚函数法给出的公式相同。因此，如果$\varepsilon_N$的值足够大，只需一步计算就可以使$g_0+x$趋于0，如此递增Lagrangian乘数法将退化为罚函数法。但是在递增Lagrangian乘数法中罚函数$\varepsilon_N$一般不是一个很大的数值，其原因将在下节讨论。
(4.1) $J=1/2K x^2 -F x -r log(g_0+x)=0$
这里$G(g)=log(x+g_0)$称为围墙（barrier）函数，因为该函数在未接触区域为一正值，但当其靠近接触面时急剧增大，趋近接触面时，其值趋向于正无穷大。它像围墙一样阻止迭代点越出约束边界。围墙函数可以有不同形式，如$G(g)=1/(g_0+x)$. 另外系数r>0的作用类似于前述罚函数。但在此它是个待定变数。在下面（4.3）式中可以看到，当r趋于零时，泛函J趋于极值。
泛函J的极值条件为其变分为零，即
(4.2) $\delta J=K x\delta x -(F+r/(g_0+x)) \delta x$
基于变分$\delta x$的任意性，要求
(4.3) $\begin{eqnarray}
Kx-\lambda-F=0 \\
\lambda（g_0+x）=r
\end{eqnarray}$
在这里我们导入了一个新变量$\lambda=r/(g+x)$. 可以看出该式与方程（2.3）非常类似，当r=0时两者等价。该方程组有三个变量x，r和$\lambda$, 需要用迭代方法求解。
* 取初值r〉0，递减系数$0<\gamma<1，\lambda_0=0, x_0=0$
* 带入(4.3)式求得x和$\lambda$. 注意到这是一个非线性方程组，需要将其线性化
(4.4) $\begin{eqnarray}
K\Delta x-\Delta \lambda=F-Kx+\lambda \\
\lambda \Delta x+(g_0+x)\Delta \lambda =r-\lambda(g_0+x)
\end{eqnarray}$
采用Newton迭代法求解。从此联立方程中可以直接求得x和$\lambda$的增量（primal-dual法），也可仅以x为自变量，靠迭代求解$\lambda$的增量（primal算法）。
* 令$r=r\gamma$, 回到上一步更新x和$\lambda$，直到r趋于0。
(3.1) $J=1/2K x^2 -F x +\lambda (x+g_0)+1/2\varepsilon_N (x+g_0)^2$
注意到该泛函既包括Lagrangian乘数$\lambda$也包括罚函数$\varepsilon_N>0$。但与前述Lagrangian乘数法和罚函数法不同的是，在该算法中
* Lagrangian乘数$\lambda$不是自变量
* 罚函数$\varepsilon_N$可以不是很大的数值
-----------------------------------------------------------------
注： 递增Lagrangian乘数法可以看成是罚函数法或Lagrangian乘数法的变种。如从Lagrangian乘数法出发，我们将（2.1）改写为
$J=1/2K x^2 -F x +\lambda (x+g_0)+(\lambda^*-\lambda) (x+g_0)$
这里，我们记$\lambda^*$为Lagrangian乘数的准确值。考虑 $\lambda^*-\lambda$为Lagrangian乘数的增分并将此式与（3.1）比较，可以看出
$\varepsilon_N (x+g_0)\approx \lambda^*-\lambda$
如此Lagrangian乘数的准确值可以如下逼近
$ \lambda^*\approx\lambda+\varepsilon_N (x+g_0)$
对应于下面的（3.5）式。
-----------------------------------------------------------------
泛函J的极值条件为其变分为零，即
(3.2) $\delta J=K x\delta x+ (\lambda+\varepsilon_Nx)\delta x -(F-\varepsilon_Ng_0) \delta x=0$
在此由于Lagrangian乘数$\lambda$不是自变量，我们不需要对其变分。基于变分$\delta x$的任意性，要求
(3.3) $(K+\varepsilon_N)x+\lambda-(F+\varepsilon_Ng_0)=0$
由此可得
(3.4) $ x=(F-\varepsilon_Ng_0-\lambda)/(K+\varepsilon_N) $
但是这里$\lambda$也是待定值，我们需要用迭代方法解此方程。
* 假设Lagrangian乘数的初值$\lambda_0=0$，n=0
* 带入上式求得$x_{n+1}$
* 判断接触条件$g_0+x>=0$是否满足。如不满足，则更新$\lambda$，
(3.5) $\lambda_{n+1}\approx \lambda_n+\varepsilon_N(g_0+x)$
令n->n+1; 并返回上一步重求x。
在上述计算过程中，Lagrangian乘数$\lambda$从零出发逐渐递增，这是我们称之为递增Lagrangian乘数法的原因。同时，x的值从$ x=(F-\varepsilon_Ng_0)/(K+\varepsilon_N) $出发逐渐递减，当$g_0+x$接近于0时计算结束。
可以看出，这里给出的x的初值与上述罚函数法给出的公式相同。因此，如果$\varepsilon_N$的值足够大，只需一步计算就可以使$g_0+x$趋于0，如此递增Lagrangian乘数法将退化为罚函数法。但是在递增Lagrangian乘数法中罚函数$\varepsilon_N$一般不是一个很大的数值，其原因将在下节讨论。
1.4 内点法
内点法是从可行域（即未接触区域g>0）内某一初始内点出发，在可行域内求极值的方法。其泛函J为(4.1) $J=1/2K x^2 -F x -r log(g_0+x)=0$
这里$G(g)=log(x+g_0)$称为围墙（barrier）函数，因为该函数在未接触区域为一正值，但当其靠近接触面时急剧增大，趋近接触面时，其值趋向于正无穷大。它像围墙一样阻止迭代点越出约束边界。围墙函数可以有不同形式，如$G(g)=1/(g_0+x)$. 另外系数r>0的作用类似于前述罚函数。但在此它是个待定变数。在下面（4.3）式中可以看到，当r趋于零时，泛函J趋于极值。
泛函J的极值条件为其变分为零，即
(4.2) $\delta J=K x\delta x -(F+r/(g_0+x)) \delta x$
基于变分$\delta x$的任意性，要求
(4.3) $\begin{eqnarray}
Kx-\lambda-F=0 \\
\lambda（g_0+x）=r
\end{eqnarray}$
在这里我们导入了一个新变量$\lambda=r/(g+x)$. 可以看出该式与方程（2.3）非常类似，当r=0时两者等价。该方程组有三个变量x，r和$\lambda$, 需要用迭代方法求解。
* 取初值r〉0，递减系数$0<\gamma<1，\lambda_0=0, x_0=0$
* 带入(4.3)式求得x和$\lambda$. 注意到这是一个非线性方程组，需要将其线性化
(4.4) $\begin{eqnarray}
K\Delta x-\Delta \lambda=F-Kx+\lambda \\
\lambda \Delta x+(g_0+x)\Delta \lambda =r-\lambda(g_0+x)
\end{eqnarray}$
采用Newton迭代法求解。从此联立方程中可以直接求得x和$\lambda$的增量（primal-dual法），也可仅以x为自变量，靠迭代求解$\lambda$的增量（primal算法）。
* 令$r=r\gamma$, 回到上一步更新x和$\lambda$，直到r趋于0。
2 算法评价
2.1 算法评价基准
评价一个数值算法的基本判据有
* 精度
* 收敛性
* 计算量
* 内存使用量
* 并行算法的相容性
对于有限元计算来说，其大部分计算时间消耗于联立方程的求解计算。求解联立方程大致可分为两类* 收敛性
* 计算量
* 内存使用量
* 并行算法的相容性
* 直接法
* 以CG法为代表的迭代法
其比较如下
直接法 迭代法
------------------------------------------------------------------- ---------
* 精度 好 可控
* 收敛性 不存在此问题 依赖于问题类型
* 计算量 较大 较小
* 内存使用量 较大 较小
* 并行算法的相容性 不好 好
从上表可以看出，对于大型计算而言，迭代法是较好的选择。在本文中将考虑各种算法对迭代法的收敛性的影响。如上所述，迭代法的收敛性依赖于问题类型，其具体的指标则是联立方程A的条件数。简单地，条件数可以理解为矩阵对角项的最大值与最小值的比
$cond(A) = max|A_{diag}|/min|A_{diag}|$
条件数增大则迭代法收敛性变差，这意味着迭代次数和计算时间的增加甚至计算发散。
另外，我们希望方程是正定的，主对角集中的，这些都可以提高迭代法求解器的收敛性。
2.2 算法评价
罚函数法:缺点： 近似解，要求极大罚函数值， 使方程组条件数变差。
优点： 计算量不大，实现最简单。
Lagrangian乘数法:
缺点： 增加了计算变量（计算量增加），方程性能变差（参见式（2.3）,可以看出其导入了零对角项，该方程变为非正定方程）。
优点： 精确解。
递增Lagrangian乘数法:
缺点： 需要迭代求解Lagrangian乘数（计算量增加）。
优点： 回避了计算变量的增加；由于选用较小的罚函数值，可以回避或减缓方程组条件数的恶化；对于接触问题解析来说，可以利用此方法把非对称的接触刚性矩阵对称化，大幅度节省内存和计算时间。
内点法:
缺点： 不能处理初始点已经位于接触面内的问题；采用primal-dual算法增加了计算变量；需要迭代求解。
优点： 内点法应用于接触问题解析主要是其可以回避解析过程中，接触点—〉非接触-〉接触，即所谓接触active set的变化，因为该方法对接触面附近的所有点都加上了约束。因此可以期待其提高计算的收敛性。
3 递增Lagrangian乘数法的实装方法（待续）
*1 Augmented Lagrangian的译法有多种。因为 augment的原意为 To make (something already developed or well under way) greater, as in size, extent, or quantity （http://www.thefreedictionary.com/augment），Make (something) greater by adding to it; increase (http://www.oxforddictionaries.com/definition/english/augment). 而Augmented Lagrangian法中，Lagrangian乘数的值正是逐渐增大的，因此提出该译法。
*2 接触条件一般在接触点为零点，接触面法线方向为正向的局部坐标系下表述。由于本例题只有一个自由度，本文简单地取局部坐标系为全局坐标。
*3 由于习惯的不同，(3.1)式也可以写为$J=1/2K x^2 -F x -\lambda (x+g_0)+1/2\varepsilon_N (x+g_0)^2$。这时得到的$\lambda$与(3.1)式相比差一个负号。
