最优化方法Python计算：求解约束优化问题的拉格朗日乘子算法

从仅有等式约束的问题入手。设优化问题(7.8)
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t.}\mathbf{h}(\mathbf{x})=\mathbf{o}\end{array}(1)$$
中$f(\mathbf{x})$，$\mathbf{h}(\mathbf{x})$均在${\text{R}}^n$上二阶连续可微，且$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)$是问题(1)满足二阶充分条件的KKT点。即$\nabla f({\mathbf{x}}_0)-\nabla \mathbf{h}({\mathbf{x}}_0){\mathbf{\lambda }}_0=0$且 ∀ d ∈ { d ∣ d ≠ o , ∇ h ( x ) ⊤ d = o } \forall\boldsymbol{d}\in\{\boldsymbol{d}|\boldsymbol{d}\not=\boldsymbol{o},\nabla\boldsymbol{h}(\boldsymbol{x})^\top\boldsymbol{d}=\boldsymbol{o}\} ∀d∈{d∣d=o,∇h(x)⊤d=o}，${\mathbf{d}}^{\top }{\nabla }_{\mathbf{xx}}L({\mathbf{x}}_0,{\mathbf{\lambda }}_0)\mathbf{d}>0$。其中，$L(\mathbf{x},\mathbf{\lambda })=f(\mathbf{x})-{\mathbf{\lambda }}^{\top }\mathbf{h}(\mathbf{x})$为问题(1)的拉格朗日函数，其中的$\mathbf{\lambda }\in {\text{R}}^l$为拉格朗日乘子。
定义问题(7.8)的{\heiti{增广拉格朗日函数}}
$$L(\mathbf{x},\mathbf{\lambda },\sigma )=f(\mathbf{x})-{\mathbf{\lambda }}^{\top }\mathbf{h}(\mathbf{x})+\frac{\sigma }{2}\mathbf{h}(\mathbf{x}{)}^{\top }\mathbf{h}(\mathbf{x})(2)$$
其中，罚因子$\sigma >0$。与问题(1)的拉格朗日函数$L(\mathbf{x},\mathbf{\lambda })$相比，增广拉格朗日函数$L(\mathbf{x},\mathbf{\lambda },\sigma )$多了个罚项$\frac{\sigma }{2}\mathbf{h}(\mathbf{x}{)}^{\top }\mathbf{h}(\mathbf{x})$。对问题(1)的增广拉格朗日函数(2)，不难验证$\forall \mathbf{x}\in \Omega$，$L(\mathbf{x},\mathbf{\lambda },\sigma )=f(\mathbf{x})$。可以证明以下命题：
定理1：$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)$是问题(1)满足二阶充分条件的KKT点，则存在${\sigma }^{\prime }\ge 0$，使得对任意的$\sigma >{\sigma }^{\prime }$ $${\mathbf{x}}_0=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }L(\mathbf{x},{\mathbf{\lambda }}_0,\sigma ).$$
反之，若${\mathbf{x}}_0=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }L(\mathbf{x},{\mathbf{\lambda }}_0,\sigma )$，且$\mathbf{h}({\mathbf{x}}_0)=\mathbf{o}$，则${\mathbf{x}}_0$是问题(1)的局部最优解，即${\mathbf{x}}_0=\arg \underset{\mathbf{x}\in \Omega }{\min }f(\mathbf{x})$。
定理1将存在最优解满足二阶充分条件KKT点$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)$的等式约束优化问题(1)求解，转化为对其增广拉格朗日函数$L(\mathbf{x},{\mathbf{\lambda }}_0,\sigma )$作为目标函数的辅助无约束优化问题的求解，其中$\sigma >0$是一个充分大的实数。然而，一般情况下，拉格朗日乘子${\mathbf{\lambda }}_0$未必已知。需设法用一系列 { λ k } \{\boldsymbol{\lambda}_k\} {λk​}去逼近${\mathbf{\lambda }}_0$。
利用定理1，我们得到由初始的乘子${\mathbf{\lambda }}_1$和罚因子${\sigma }_1>0$出发，计算问题(1)的满足二阶充分条件的KKT点$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)$近似值的迭代方法：计算
$$\{\begin{array}{l}{\mathbf{x}}_{k+1}=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }L(\mathbf{x},{\mathbf{\lambda }}_k,{\sigma }_k)\\ {\mathbf{\lambda }}_{k+1}={\mathbf{\lambda }}_k-{\sigma }_k\mathbf{h}({\mathbf{x}}_k)\end{array}.$$
对给定的容错误差$\epsilon >0$，若$\Vert \mathbf{h}({\mathbf{x}}_{k+1})\Vert <\epsilon$，则停止迭代。$\left(\begin{array}{c}{\mathbf{x}}_{k+1}\\ {\mathbf{\lambda }}_{k+1}\end{array}\right)$即为$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\end{array}\right)$近似值。\par
否则，即$\Vert \mathbf{h}({\mathbf{x}}_{k+1})\Vert \ge \epsilon$，用给定的放大系数$c>1$计算
$${\sigma }_{k+1}=\{\begin{array}{ll}{\sigma }_k& \frac{\Vert {\mathbf{h}}_{k+1}\Vert }{\Vert {\mathbf{h}}_k\Vert }<1\\ c\times {\sigma }_k& \frac{\Vert {\mathbf{h}}_{k+1}\Vert }{\Vert {\mathbf{h}}_k\Vert }\ge 1\end{array},$$
进入下一轮迭代。
以上讨论的对等式约束优化问题的乘子算法是由 Powell和Hestenes于1969年提出的，常称为PH算法。
对仅有不等式约束
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t. }\mathbf{g}(\mathbf{x})\ge \mathbf{o}\end{array}$$
其中，$f:{\text{R}}^n\to \text{R}$，$\mathbf{g}:{\text{R}}^n\to {\text{R}}^m$在${\text{R}}^n$上二阶连续可微。可行域 Ω = { x ∣ g ( x ) ≥ o } \Omega=\{\boldsymbol{x}|\boldsymbol{g}(\boldsymbol{x})\geq\boldsymbol{o}\} Ω={x∣g(x)≥o}。添加松弛变量$\mathbf{s}\ge \mathbf{o}$，构造等价问题
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t. }\mathbf{g}(\mathbf{x})-\mathbf{s}=\mathbf{o}\\ \mathbf{s}\ge \mathbf{o}\end{array}$$
通过变换，上述问题可转换成等价的
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t.}{\mathbf{h}}_1(\mathbf{x})=\mathbf{o}\end{array}.$$
其中，${\mathbf{h}}_1(\mathbf{x})=\min \left(\frac{1}{\sigma }\mathbf{\mu },\mathbf{g}(\mathbf{x})\right)$。
相仿地，对一般的优化问题
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t. }\mathbf{h}(\mathbf{x})=\mathbf{o}\\ \mathbf{g}(\mathbf{x})\ge \mathbf{o}\end{array}(3)$$
其中$f(\mathbf{x})$、$\mathbf{h}(\mathbf{x})$和$\mathbf{g}(\mathbf{x})$在${\text{R}}^n$上二阶连续可微，可行域 Ω = { x ∣ h ( x ) = o , g ( x ≥ o ) } \Omega=\{\boldsymbol{x}|\boldsymbol{h}(\boldsymbol{x})=\boldsymbol{o},\boldsymbol{g}(\boldsymbol{x}\geq\boldsymbol{o})\} Ω={x∣h(x)=o,g(x≥o)}。为求解其满足二阶充分条件的KKT点$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\\ {\mathbf{\mu }}_0\end{array}\right)$，将其转换为与之等价的问题
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t. }\mathbf{h}(\mathbf{x})=\mathbf{o}\\ {\mathbf{h}}_1(\mathbf{x})=\mathbf{o}\end{array}.(4)$$
其增广拉格朗日函数为
$$L(\mathbf{x},\mathbf{\lambda },\mathbf{\mu },\sigma )=f(\mathbf{x})-{\mathbf{\lambda }}^{\top }\mathbf{h}(x)-{\mathbf{\mu }}^{\top }{\mathbf{h}}_1(\mathbf{x})+\frac{\sigma }{2}(\mathbf{h}(\mathbf{x}{)}^{\top }\mathbf{h}(\mathbf{x})+{\mathbf{h}}_1(\mathbf{x}{)}^{\top }{\mathbf{h}}_1(\mathbf{x}))$$
对初始乘子${\mathbf{\lambda }}_1$，${\mathbf{\mu }}_1$，罚因子${\sigma }_1$，用拉格朗日乘子法计算问题(7.10)满足二阶充分条件的KKT点$\left(\begin{array}{c}{\mathbf{x}}_0\\ {\mathbf{\lambda }}_0\\ {\mathbf{\mu }}_0\end{array}\right)$的近似值，迭代式为
$$\{\begin{array}{l}{\mathbf{x}}_{k+1}=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }L(\mathbf{x},{\mathbf{\lambda }}_k,{\mathbf{\mu }}_k,{\sigma }_k)\\ {\mathbf{\lambda }}_{k+1}={\mathbf{\lambda }}_k-{\sigma }_k\mathbf{h}({\mathbf{x}}_k)\\ {\mathbf{\mu }}_{k+1}=\max (\mathbf{o},{\mathbf{\mu }}_k-{\sigma }_k\mathbf{g}({\mathbf{x}}_k))\end{array},$$
对给定的容错误差$\epsilon >0$，记$\beta =\Vert \mathbf{h}({\mathbf{x}}_k)\Vert +\Vert {\mathbf{h}}_1({\mathbf{x}}_k)\Vert$，则迭代终止条件为
$$\beta <\epsilon ,$$
罚因子扩大条件为
$$\frac{{\beta }_{new}}{{\beta }_{old}}=\frac{\Vert \mathbf{h}({\mathbf{x}}_{k+1})\Vert +\Vert {\mathbf{h}}_1({\mathbf{x}}_{k+1})\Vert }{\Vert \mathbf{h}({\mathbf{x}}_k)\Vert +\Vert {\mathbf{h}}_1({\mathbf{x}}_k)\Vert }\ge \eta .$$
其中$\eta \in (0,1)$。这一算法是Rockfellar于1973年，在PH算法的基础上提出的，常称为PHR算法。下列代码实现PHR算法。

import numpy as np											#导入numpy
from scipy.optimize import minimize, OptimizeResult			#导入minimize和OptimizeResult
def phr(f, x1, h = None, g = None, eps = 1e-5):
    k = 0
    sigmak = 10
    c = 2.5
    eta = 0.8
    if callable(h):											#有等式约束
        l = h(x1).size
        lamdak = 0.1 * np.ones(l)
    else:													#无等式约束
        lamdak = np.zeros(1)
        h = lambda x: np.zeros(1)
    if callable(g):											#有不等式约束
        m = g(x1).size
        muk = 0.1 * np.ones(m)
        h1 = lambda x: np.minimum(muk / sigmak, g(x))
    else:													#无不等式约束
        muk = np.zeros(1)
        h1 = g = lambda x: np.zeros(1)
        m = 1
    zero = np.zeros(m)										#零向量
    P = lambda x: np.dot(h(x), h(x)) + np.dot(h1(x), h1(x))	#罚函数
    L = lambda x: f(x) - np.dot(lamdak, h(x)) - \			#增广拉格朗日函数
                  np.dot(muk, h1(x)) + 0.5 * sigmak * P(x)
    xk = x1													#迭代点
    bnew = bold = np.linalg.norm(h(xk)) + np.linalg.norm(h1(xk))#新、旧beta值
    while bnew >= eps:
        k += 1
        xk = minimize(L, xk).x								#更新迭代点
        lamdak = lamdak - sigmak * h(xk)					#更新乘子lambda
        muk = np.maximum(zero, muk - sigmak * g(xk))		#更新乘子mu
        if bnew / bold >= eta:								#须扩大罚因子
            sigmak *= c
        bold = bnew											#保存beta
        bnew = np.linalg.norm(h(xk)) + np.linalg.norm(h1(xk))#更新beta
    return OptimizeResult(fun = f(xk), x = xk, nit = k)

程序的第3~37行定义的函数phr实现图7.4所示的乘子算法PHR。该函数的5个参数f，x1，h，g和eps分别表示优化问题(3)的目标函数$f(\mathbf{x})$，初始迭代点${\mathbf{x}}_1$，等式约束函数$\mathbf{h}(\mathbf{x})$，不等式约束函数$\mathbf{h}(\mathbf{x})$和容错误差$\epsilon$。其中，h，g和eps的缺省值分别为None、None和$10^{-5}$。
函数体内，第4~7行分别初始化迭代次数$k$，罚因子${\sigma }_k$，罚因子放大系数$c$和罚因子放大判别值$\eta$。
第8~13行的if-else语句根据是否有等式约束，初始化乘子$\mathbf{\lambda }$：若是，初始化为含有$l$个0.1的向量，否则置为0。相仿地，第14~21行的if-else语句设置乘子$\mathbf{\mu }$函数及${\mathbf{h}}_1(\mathbf{x})=\min \left\{\frac{\mathbf{\mu }}{{\sigma }_k},\mathbf{g}(\mathbf{x})\right\}$。
第23~25行定义罚函数和增广拉格朗日函数
$$\begin{array}{l}P(\mathbf{x})=\mathbf{h}(\mathbf{x}{)}^{\top }\mathbf{h}(\mathbf{x})+{\mathbf{h}}_1(\mathbf{x}{)}^{\top }{\mathbf{h}}_1(\mathbf{x})\\ L(\mathbf{x})=f(\mathbf{x})-{\mathbf{\lambda }}^{\top }\mathbf{h}(x)-{\mathbf{\mu }}^{\top }{\mathbf{h}}_1(\mathbf{x})+\frac{{\sigma }_k}{2}P(\mathbf{x})\end{array}.$$
第26行初始化迭代点xk，第27行将${\beta }_{new}$和${\beta }_{old}$初始化为$\Vert \mathbf{h}({\mathbf{x}}_k)\Vert +\Vert {\mathbf{h}}_1({\mathbf{x}}_k)\Vert$。
第28~36行的while循环执行算法的迭代。其中，第30~32行按
$$\{\begin{array}{l}{\mathbf{x}}_{k+1}=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }L(\mathbf{x},{\mathbf{\lambda }}_k,{\mathbf{\mu }}_k,{\sigma }_k)\\ {\mathbf{\lambda }}_{k+1}={\mathbf{\lambda }}_k-{\sigma }_k\mathbf{h}({\mathbf{x}}_k)\\ {\mathbf{\mu }}_{k+1}=\max (\mathbf{o},{\mathbf{\mu }}_k-{\sigma }_k\mathbf{g}({\mathbf{x}}_k))\end{array},$$
更新迭代点xk，乘子lamdk和muk。第33~34行根据是否
$$\frac{{\beta }_{new}}{{\beta }_{old}}\ge \eta$$
决定扩大罚因子${\sigma }_k$。第35、36行更新${\beta }_{old}$和${\beta }_{new}$，以备进入下一轮迭代。一旦测得
$${\beta }_{new}<\epsilon$$
迭代结束，第37行返回最优解近似值xk，最优值近似值f(xk)和迭代次数k。
例1：用程序7.11定义的phr函数求解优化问题
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})=(x_1-2{)}^2+(x_2-1{)}^2\\ \text{s.t. }{x}_1-2x_2+1=0\\ \frac{1}{4}{x}_1^2+x_2^2\le 1\end{array},$$
给定初始迭代点${\mathbf{x}}_1=\left(\begin{array}{c}3\\ 3\end{array}\right)$。
解：本题中目标函数$f(\mathbf{x})=(x_1-2{)}^2+(x_2-1{)}^2$，等式约束函数$h(\mathbf{x})=x_1-2x_2+1$，不等式约束函数$g(\mathbf{x})=1-\frac{1}{4}{x}_1^2-x_2^2$。下列代码完成计算。

import numpy as np								#导入numpy
from utility import phr							#导入phr
np.set_printoptions(precision=4)				#设置输出精度
f = lambda x: (x[0] - 2) ** 2 + (x[1] - 1) ** 2	#目标函数
h = lambda x: x[0] - 2 * x[1] + 1				#等式约束函数
g = lambda x: 1 - 0.25 * x[0] ** 2 - x[1] ** 2	#不等式约束函数
x1 = np.array([3, 3])							#初始迭代点
print(phr(f, x1, h, g))

利用代码内注释信息不难理解程序。运行程序，输出

 fun: 1.3934640206427056
 nit: 5
   x: array([0.8229, 0.9114])

意为phr用5次迭代，以$10^{-5}$的容错误差，算得最优解近似值${\mathbf{x}}_0\approx \left(\begin{array}{c}0.8229\\ 0.9114\end{array}\right)$，最优值的近似值$f({\mathbf{x}}_0)\approx 1.3934$。
例2：用phr函数求解求解线性规划
$$\{\begin{array}{l}\text{minimize}{x}_1-x_2\\ \text{s.t. }-x_1+2x_2+x_3\le 2\\ -4x_1+4x_2-x_3=4\\ x_1-x_3=0\\ x_1,x_2,x_3\ge 0\end{array},$$
给定初始迭代点${\mathbf{x}}_1=\mathbf{o}$。\par
解：这个线性规划的数据矩阵为
$$\mathbf{c}=\left(\begin{array}{c}1\\ -1\\ 0\end{array}\right),{\mathbf{A}}_{eq}=\left(\begin{array}{ccc}-4& 4& -1\\ 1& 0& -1\end{array}\right),{\mathbf{b}}_{eq}=\left(\begin{array}{c}4\\ 0\end{array}\right),{\mathbf{A}}_{iq}=\left(\begin{array}{ccc}1& -2& -1\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),{\mathbf{b}}_{iq}=\left(\begin{array}{c}-2\\ 0\\ 0\\ 0\end{array}\right).$$
于是
$$f(\mathbf{x})={\mathbf{c}}^{\top }\mathbf{x},\mathbf{h}(\mathbf{x})={\mathbf{A}}_{eq}-{\mathbf{b}}_{eq},\mathbf{g}(\mathbf{x})={\mathbf{A}}_{iq}-{\mathbf{b}}_{ip}.$$，下列程序完成计算。

import numpy as np					#导入numpy
c = np.array([1, -1, 0])			#向量c
Ae = np.array([[-4, 4, -1],			#矩阵Aeq
               [1, 0, -1]])
be = np.array([4, 0])				#向量beq
Ai = np.array([[1, -2, -1],			#矩阵Aiq
               [1, 0, 0],
               [0, 1, 0],
               [0, 0, 1]])
bi = np.array([-2, 0, 0, 0])		#向量biq
f = lambda x: np.dot(c, x)			#目标函数
h = lambda x: np.matmul(Ae, x) - be	#等式约束函数
g = lambda x: np.matmul(Ai, x) - bi	#不等式约束函数
x1 = np.zeros(3)					#初始迭代点
print(phr(f, x1, h, g))

运行程序，输出

 fun: -0.9999999986126238
 nit: 2
   x: array([-5.4602e-08,  1.0000e+00, -8.5667e-09])

意为phr只用了2次迭代，即得到最优解${\mathbf{x}}_0=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$，最优值$f({\mathbf{x}}_0)=-1$。而对同一问题，sumt（见博文《最优化方法Python计算：求解约束优化问题的罚函数算法》）却用了8次迭代获得同样的结果。
例3：用phr函数求解二次规划
$$\{\begin{array}{l}\text{minimize}{x}_1^2+x_1{x}_2+2x_2^2+x_3^2-6x_1-2x_2-12x_3\\ \text{s.t. }{x}_1+x_2+x_3-2=0\\ x_1-2x_2+3\ge 0\\ x_1,x_2,x_3\ge 0\end{array},$$
给定初始可行解${\mathbf{x}}_1=\left(\begin{array}{c}1\\ 1\\ 0\end{array}\right)$。
解：本二次规划的数据矩阵为
$$\begin{array}{l}\mathbf{H}=\left(\begin{array}{ccc}2& 1& 0\\ 1& 4& 0\\ 0& 0& 2\end{array}\right),\mathbf{c}=\left(\begin{array}{c}-6\\ -2\\ -12\end{array}\right)\\ {\mathbf{A}}_{eq}=\left(\begin{array}{ccc}1& 1& 1\end{array}\right),{\mathbf{b}}_{eq}=(2)\\ {\mathbf{A}}_{iq}=\left(\begin{array}{ccc}1& -2& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right),{\mathbf{b}}_{iq}=\left(\begin{array}{c}-3\\ 0\\ 0\\ 0\end{array}\right)\end{array}$$
下列代码完成计算。

import numpy as np									#导入numpy
H = np.array([[2, 1, 0],							#矩阵H
              [1, 4, 0],
              [0, 0, 2]])
c = np.array([-6, -2, -12])							#向量c
Ae = np.array([[1, 1, 1]])							#矩阵Aeq
be = np.array([2])									#向量beq
Ai = np.array([[-1, -2, 0],							#矩阵Aiq
               [1, 0, 0],
               [0, 1, 0],
               [0, 0, 1]])
bi = np.array([-3, 0, 0, 0])						#向量biq
x1 = np.array([1, 1, 0])							#初始迭代点
f = lambda x: 0.5 * np.matmul(x, np.matmul(H, x)) + np.matmul(c, x)	#目标函数
h = lambda x: np.matmul(Ae, x) - be					#等式约束函数
g = lambda x: np.matmul(Ai, x) - bi					#不等式约束函数
print(phr(f, x1, h, g))
\end{lstlisting}\par

运行程序，输出

 fun: -20.000001376228557
 nit: 7
   x: array([-9.6384e-08, -1.2540e-07,  2.0000e+00])

即phr函数用7次迭代，算得该二次规划的最优解${\mathbf{x}}_0=\left(\begin{array}{c}0\\ 0\\ 2\end{array}\right)$，最优值$f({\mathbf{x}}_0)=-20$。相较于sumt函数（见博文《最优化方法Python计算：求解约束优化问题的罚函数算法》）用16次迭代获得同样结果，效率提高是明显的。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！