最优化方法Python计算：求解约束优化问题的罚函数算法

设等式约束优化问题
$$\{\begin{array}{l}\text{minimize}f(\mathbf{x})\\ \text{s.t. }\mathbf{h}(\mathbf{x})=\mathbf{o}\end{array}(1)$$
其中$f:{\text{R}}^n\to \text{R}$，$\mathbf{h}:{\text{R}}^n\to {\text{R}}^l$，在${\text{R}}^n$上连续。可行域 Ω = { x ∣ h ( x ) = o } \Omega=\{\boldsymbol{x}|\boldsymbol{h}(\boldsymbol{x})=\boldsymbol{o}\} Ω={x∣h(x)=o}。构造罚函数
$$P(\mathbf{x})=\frac{1}{2}\mathbf{h}(\mathbf{x}{)}^{\top }\mathbf{h}(\mathbf{x})$$
称${\text{R}}^n\to \text{R}$函数
$$F(\mathbf{x},\sigma )=f(\mathbf{x})+\sigma P(\mathbf{x})$$
为问题(1)的{\heiti{增广目标函数}}，其中$\sigma >0$，称为罚因子。注意到$F(\mathbf{x},\sigma )=f(\mathbf{x})$，当且仅当$\mathbf{x}\in \Omega$。考虑无约束优化问题
$$\{\begin{array}{l}\text{minimize}F(\mathbf{x},\sigma )\\ \text{s.t. }\mathbf{x}\in {\text{R}}^n\end{array},(2)$$
称为(1)的子问题。给定罚因子$\sigma >0$，设子问题(2)最优解${\mathbf{x}}_{\sigma }=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }F(\mathbf{x},\sigma )$。若${\mathbf{x}}_{\sigma }\in \Omega$，则$\sigma P({\mathbf{x}}_{\sigma })=0$，且$F({\mathbf{x}}_{\sigma },\sigma )=f({\mathbf{x}}_{\sigma })$。故${\mathbf{x}}_{\sigma }$可以视为$f(\mathbf{x})$最优解的近似值。否则，即${\mathbf{x}}_{\sigma }\notin \Omega$，则$\sigma P({\mathbf{x}}_{\sigma })$是一个正数，其存在是对${\mathbf{x}}_{\sigma }$脱离$\Omega$的一种“惩罚”。此时，若加大${\sigma }^{\prime }>\sigma$，再次尝试求解子问题${\mathbf{x}}_{{\sigma }^{\prime }}=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }F(\mathbf{x},{\sigma }^{\prime })$，意欲迫使${\mathbf{x}}_{{\sigma }^{\prime }}$向$\Omega$靠拢。可以证明以下命题
定理1：若约束优化问题(1)存在唯一最优解${\mathbf{x}}_0$，子问题(2)对任意$\sigma >0$，均有最优解${\mathbf{x}}_{\sigma }$。给定序列$0<{\sigma }_1<{\sigma }_2<\cdots <{\sigma }_k<\cdots$，且$\underset{k\to \infty }{\lim }{\sigma }_k=+\infty$。按迭代式
$${\mathbf{x}}_k=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }F(\mathbf{x},{\sigma }_k)$$
算得序列 { x k } \{\boldsymbol{x}_k\} {xk​}，则必有
$$\{\begin{array}{l}\underset{k\to \infty }{\lim }{\mathbf{x}}_k={\mathbf{x}}_0\\ \underset{k\to \infty }{\lim }{\sigma }_{k}P({\mathbf{x}}_k)=0\end{array}.$$
对一般的有约束优化问题
$$\{\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:{\text{R}}^n\to \text{R}$，$\mathbf{h}:{\text{R}}^n\to {\text{R}}^l$，$\mathbf{g}:{\text{R}}^n\to {\text{R}}^m$。$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}。定义罚函数
$$P(\mathbf{x})=\mathbf{h}(\mathbf{x}{)}^{\top }\mathbf{h}(\mathbf{x})+{\mathbf{g}}_1(\mathbf{x}{)}^{\top }{\mathbf{g}}_1(\mathbf{x}).$$
其中，${\mathbf{g}}_1(\mathbf{x})=\min (\mathbf{o},\mathbf{g}(\mathbf{x}))=\left(\begin{array}{c}\min (0,g_1(\mathbf{x}))\\ ⋮\\ \min (0,g_m(\mathbf{x}))\end{array}\right)$。不难证实
$$\{\begin{array}{ll}P(\mathbf{x})=0& \mathbf{x}\in \Omega \\ P(\mathbf{x})>0& \mathbf{x}\notin \Omega \end{array},$$
且在${\text{R}}^n$上连续。取罚因子$\sigma >0$，定义增广目标函数：
$$F(\mathbf{x},\sigma )=f(\mathbf{x})+\sigma P(\mathbf{x}).$$
及子问题
$$\{\begin{array}{l}\text{minimize}F(\mathbf{x},\sigma )\\ \text{s.t.}\mathbf{x}\in {\text{R}}^n\end{array}(4)$$
则关于等式约束优化问题(1)及其子问题(2)的定理1的结论对问题(3)及其子问题(4)均成立。根据定理1，对问题(4)给定罚因子初始值${\sigma }_1>0$（典型值为10）、放大系数$c>1$（典型值为2.5）及容错误差$\epsilon >0$，按
$${\mathbf{x}}_{k+1}=\arg \underset{\mathbf{x}\in {\text{R}}^n}{\min }F(\mathbf{x},{\sigma }_k)$$
做迭代，直至${\sigma }_{k}P({\mathbf{x}}_k)<\epsilon$为止。所得当前${\mathbf{x}}_k$即为问题(3)的最优解近似值。由于罚函数是关于约束函数的二次式，故此算法常称为二次罚函数算法。该算法以一系列子问题的最优解来逼近问题的最优解，也被称为序列无约束极小化方法简记为SUMT方法。
下列代码实现罚函数算法。

import numpy as np													#导入numpy
from scipy.optimize import minimize,OptimizeResult					#导入minimize,OptimizeResult
def sumt(f, x1, h = None, g = None, eps = 1e-5):
    if not callable(h):												#没有等式约束
        h = lambda x: np.zeros(1)									#置零
    if callable(g):
        m = g(x1).size												#g的维数
        zero = np.zeros(m)
        g1 = lambda x: np.minimum(zero, g(x))						#转换不等式约束
    else:															#没有不等式约束
        g1 = lambda x: np.zeros(1)									#置零
    sigmak = 10.0													#初始罚因子
    c = 2.5															#罚因子放大系数
    sPk = 10.0														#初始罚项
    k = 0															#初始迭代数
    xk = x1															#初始迭代解向量
    P = lambda x: 0.5 * (np.dot(h(x), h(x)) + np.dot(g1(x), g1(x)))	#罚函数
    F = lambda x: f(x) + sigmak * P(x)								#子问题目标函数
    while sPk >= eps:												#迭代
        k += 1														#迭代次数自增1
        result = minimize(F,xk)										#求解子问题
        xk = result.x												#子问题最优解
        sPk = result.fun - f(xk)									#罚项
        sigmak *= c													#放大罚因子
    return OptimizeResult(fun = f(xk), x = xk, nit = k)

程序的第3~25行定义的sumt函数实现求解一般约束优化问题的二次罚函数算法（序列无约束极小化SUMT方法）。该函数的参数f，x1表示约束优化问题的目标函数$f(\mathbf{x})$和初始解向量${\mathbf{x}}_1$。h表示等式约束函数$\mathbf{h}(\mathbf{x})$缺省值为None，表示优化问题无等式约束。g表示不等式约束函数$\mathbf{g}(\mathbf{x})$，缺省值为None，表示问题无不等式约束。参数eps表示容错误差$\epsilon$。
函数体内第4~5行对无等式约束的情况将h置为零函数，以便统一代码。第6~11行的if-else语句根据是否有不等式约束，设置表示${\mathbf{g}}_1(\mathbf{x})=\min (\mathbf{g}(\mathbf{x}),\mathbf{o})$便于统一代码的形式。第12~16行分别对罚因子${\sigma }_k$、罚因子放大系数$c$、罚项${\sigma }_{k}P({\mathbf{x}}_k)$、迭代次数$k$及解向量${\mathbf{x}}_k$作初始化操作。第19~24行的while语句完成算法的迭代。其中，第21行调用scipy.optimize模块中的minimize函数求解子问题，返回值赋予result。第22行读取迭代解向量xk，第23行计算罚项${\sigma }_{k}P({\mathbf{x}}_k)=F({\mathbf{x}}_k,{\sigma }_k)-f({\mathbf{x}}_k)$，赋予sPk。第24行放大罚因子${\sigma }_k=c\times {\sigma }_k$。迭代结束，第25行返回问题的最优解近似值${\mathbf{x}}_k$、最优值近似值$f({\mathbf{x}}_k)$和迭代次数$k$。
例1：用函数sumt求解求解线性规划
$$\{\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}$。
解：这个线性规划的数据矩阵为
$$\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(sumt(f, x1, h, g))

借助代码内的注释信息，不难理解本程序。第15行调用sumt函数，计算问题最优解近似值并输出。运行程序，输出

 fun: -1.0000094540269446
 nit: 8
   x: array([-6.39954354e-06,  1.00000305e+00,  6.29523189e-06])

意味着经过8次迭代，算得本问题的最优解为${\mathbf{x}}_0=\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right)$，最优解处函数值为$f({\mathbf{x}}_0)=-1$。
例2：用sumt函数求解二次规划
$$\{\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(sumt(f, x1, h, g))

对照代码内的注释信息，不难理解程序。运行程序，输出

 fun: -20.000011166902542
 nit: 16
   x: array([-2.22198421e-07, -6.51694981e-07,  2.00000173e+00])

这意味着经过16次迭代，算得本问题的最优解${\mathbf{x}}_0=\left(\begin{array}{c}0\\ 0\\ 2\end{array}\right)$，最优值$f({\mathbf{x}}_0)=-20$。
例3：用sumt函数求解优化问题
$$\{\begin{array}{l}\text{minimize}{x}_1^2+x_2^2-16x_1-10x_2\\ \text{s.t. }-x_1^2+6x_1-4x_2+11\ge 0\\ x_1{x}_2-3x_2-e^{x_1-3}+1\ge 0\\ x_1,x_2\ge 0\end{array},$$
给定初始迭代点${\mathbf{x}}_1=\mathbf{o}$。
解：下列代码完成计算：

import numpy as np												#导入numpy
f = lambda x: x[0] ** 2 + x[1] ** 2 - 16 * x[0] - 10 * x[1]		#目标函数
g = lambda x: np.array([-x[0] ** 2 + 6 * x[0] - 4 * x[1] + 11,	#不等式约束函数
                        x[0] * x[1] - 3 * x[1] - np.exp(x[0] - 3) + 1,
                        x[0],
                        x[1]])
x1 =  np.zeros(2)												#初始迭代点
print(sumt(f, x1, g = g))

借助代码内注释信息，不难理解本程序。需要注意，本问题只有不等式约束g，第8行调用sumt时应表示等式约束的参数h使用缺省值，在其后的不等式约束参数g须使用赋值形式传递。运行程序，输出

 fun: -79.80782889081777
 nit: 11
   x: array([5.23961012, 3.74603874])

意味着经过11次迭代，算得线性规划问题最优解近似值${\mathbf{x}}_0=\left(\begin{array}{c}5.3296\\ 3.7460\end{array}\right)$，最优值$f({\mathbf{x}}_0)\approx -79.8078$。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！