约束优化方法--惩罚函数法

基本介绍

考虑一般的约束优化问题
$$\begin{array}{rl}& \min f(\mathbf{x})\\ & s.t.g_i(\mathbf{x})\le 0,i=1,2,...,m,\\ & h_j(\mathbf{x})=0,j=1,2,...,p,\\ & \mathbf{x}\in X\subset {\mathbf{R}}^{\mathbf{n}},\mathbf{x}=(x_1,x_2,...,x_n{)}^T.\end{array}$$
其中$f,g_i,h_j\in C^1$.

惩罚函数法的基本思想是：把约束问题转化为一个或一系列无约束问题求解, 所以也称为序列无约束极小化技术(Sequential Unconstrained Minimization Technique), 简称SUMT法.

具体思路：根据约束的特点，构造某种惩罚函数加入到目标函数中建立无约束问题，利用惩罚策略，对于无约束问题求解过程中企图违反约束的迭代点给予很大的目标函数值，迫使无约束问题的极小点无限靠近可行域或保持在可行域内移动，直至迭代点列收敛到原约束问题极小点。

构造方法

下面以仅含等式约束的优化问题为例给出相应的构造方法：
min ⁡ f ( x ) , x ∈ R n s . t c i ( x ) = 0 , i ∈ E = { 1 , . . . , l } . \begin{align*} & \min f(\mathbf{x}),\mathbf{x}\in \mathbf{R^n}\\ & s.t \quad c_i(\mathbf{x})=0,i\in E=\{1,...,l\}. \end{align*} ​minf(x),x∈Rns.tci​(x)=0,i∈E={1,...,l}.​
构造形如
$$F(\mathbf{x},M)=f(\mathbf{x})+M\sum _{i=1}^l|c_i(x){|}^{\beta },\beta \ge 1$$
的函数, 其中$M>0$为参数, 称为罚因子.

若记原约束问题的可行域为$D$, 对于惩罚项$p(\mathbf{x})={\sum }_{i=1}^l|c_i(x){|}^{\beta },\beta \ge 1$,

• 当$\mathbf{x}$为可行解, 即$\mathbf{x}\in D$时, $c_i(x)=0,p(x)=0,F(\mathbf{x},M)=f(x)$, 不受惩罚；
• 当$\mathbf{x}$不是可行解时, 即$\mathbf{x}\notin $D时, $c_i(x)\ne 0,p(x)>0,F(\mathbf{x},M)=f(x)+p(x)>0$, 即约束条件被破坏是一种惩罚, 且$M$越大, 惩罚越重.

此时, 为使$F(\mathbf{x},M)$取得极小值, $p(x)$应充分小, 即$F(\mathbf{x},M)$的极小点应充分逼近可行域, 当然也希望它能够逼近最优解.

从上述过程中, 我们可以体会到惩罚项及罚函数的意义, 即惩罚破坏约束的迭代点从而使目标函数的极小点逼近可行域或在可行域中移动.

此外, 我们还可以得出惩罚项构造要满足的条件：

1. $p(\mathbf{x})$连续；
2. $\forall \mathbf{x}\in {\mathbf{R}}^{\mathbf{n}},p(\mathbf{x})\ge 0$;
3. $\forall \mathbf{x}\in S,p(\mathbf{x})=0$.

下面考虑不等式约束的优化问题：
min ⁡ f ( x ) , x ∈ R n s . t . c i ( x ) ≥ 0 , i ∈ I = { 1 , 2 , . . . , m } . \begin{align*} & \min f(\mathbf{x}), \mathbf{x}\in \mathbf{R^n}\\ & s.t. \quad c_i(\mathbf{x})\ge 0, i\in I=\{1,2,...,m\}. \end{align*} ​minf(x),x∈Rns.t.ci​(x)≥0,i∈I={1,2,...,m}.​
类似地构造如下函数
$$F(\mathbf{x},M)=f(\mathbf{x})+Mp(\mathbf{x}),M>0,$$
其中
$$p(\mathbf{x})=\{\begin{array}{ll}& 0,c_i(\mathbf{x})\ge 0,\\ & {\sum }_{i=1}^m|c_i(\mathbf{x}){|}^{\alpha },\alpha \ge 1,c_i(\mathbf{x})<0.\end{array}$$
也可以写成
$$p(\mathbf{x})=\sum _{i=1}^m|\min (0,c_i(x)){|}^{\alpha }=\sum _{i=1}^m(\frac{|c_i(x)|-c_i(x)}{2}{)}^{\alpha },\alpha \ge 1.$$
同样地, 惩罚项仍符合前述的惩罚策略.

最后, 对于一般的约束优化问题
$$\begin{array}{rl}& \min f(\mathbf{x})\\ & s.t.c_i(\mathbf{x})=0,i\in E=1,2,...,l,\\ & c_i(\mathbf{x})\ge 0,i\in I=l+1,...,m.\end{array}$$
其可行域为 D = { x ∈ R n ∣ c i ( x ) = 0 , i ∈ E ; c i ( x ) ≥ 0 , i ∈ I } D=\{\mathbf{x}\in \mathbf{R^n}|c_i(x)=0,i\in E; c_i(x)\ge 0,i\in I\} D={x∈Rn∣ci​(x)=0,i∈E;ci​(x)≥0,i∈I}.

构造如下函数
$$F(\mathbf{x},M)=f(\mathbf{x})+Mp(\mathbf{x}),M>0,(1)$$
其中
$$p(\mathbf{x})=\sum _{i=1}^l|c_i(\mathbf{x}){|}^{\beta }+\sum _{j=l+1}^m|\min (0,c_j(\mathbf{x})){|}^{\alpha },\alpha \ge 1,\beta \ge 1.(2)$$
显然, 当$\mathbf{x}\in D$时, $p(\mathbf{x})=0$; 当$\mathbf{x}\notin D$时, $p(\mathbf{x})>0$.

式$(1)$称为约束问题的增广目标函数, $(2)$是一般约束问题的罚函数, $M>0$是罚因子, 通常取 $\alpha =\beta =2$.

也就是说, 一般约束问题的求解转化为求增广目标函数的无约束极小, 即 min ⁡ F ( x , M k ) , { M k } \min F(\mathbf{x},M_k),\{M_k\} minF(x,Mk​),{Mk​}为正的数列, 且 { M k } → ∞ \{M_k\}\rightarrow \infty {Mk​}→∞. 此外, 还要注意的是, 无约束问题最优解的极限是原问题的最优解.

应用举例

例1

求解约束问题
$$\begin{array}{rl}& \min f(x)=(x_1-3{)}^2+(x_2-2{)}^2\\ & s.t.h(x)=x_1+x_2-4=0.\end{array}$$
解：
$$F(\mathbf{x},M)=(x_1-3{)}^2+(x_2-2{)}^2+M(x_1+x_2-4{)}^2.$$
则有
$$\begin{array}{rl}& \frac{\partial F}{\partial x_1}=2(x_1-3)+2M(x_1+x_2-4),\\ & \frac{\partial F}{\partial x_2}=2(x_2-2)+2M(x_1+x_2-4).\end{array}$$
令$\frac{\partial F}{\partial x_1}=0,\frac{\partial F}{\partial x_2}=0$, 得
$$x_1=\frac{5M+3}{2M+1},x_2=\frac{3M+2}{2M+1}.$$
又由$\frac{{\partial }^{2}F}{\partial x_1^2}=2(M+1),\frac{{\partial }^{2}F}{\partial x_2^2}=2(M+1),\frac{{\partial }^{2}F}{\partial x_1\partial x_2}=2M$，可得
$${\nabla }^{2}F=\left[\begin{array}{cc}2(M+1)& 2M\\ 2M& 2(M+1)\end{array}\right].$$
$M>0\Rightarrow {\nabla }^{2}F$正定. 因此, $F(\mathbf{x},M)$在$(\frac{5M+3}{2M+1},\frac{3M+2}{2M+1}{)}^T$取得极小值. 令$M\to +\infty$, 则原约束问题的最优解为
$$x^{\ast }=\underset{M\to +\infty }{\lim }\mathbf{x}(M)=(\frac{5}{2},\frac{3}{2}{)}^T.$$

例2

求解约束问题
$$\begin{array}{rl}& \min f(x)=x_1^2+x_2^2\\ & s.t.x_1+1\le 0.\end{array}$$
解：

首先要将约束条件改写成$-x_1-1\ge 0$！
$$F(\mathbf{x},M)=x_1^2+x_2^2+M[\min (0,-x_1-1){]}^2=\{\begin{array}{l}{x}_1^2+x_2^2,x_1+1\le 0,\\ x_1^2+x_2^2+M(x_1+1{)}^2,x_1+1>0.\end{array}$$
故
$$\frac{\partial F}{\partial x_1}=\{\begin{array}{l}2x_1,x_1<-1,\\ 2x_1+2M(x_1+1),x_1>-1.\end{array}\frac{\partial F}{\partial x_2}=2x_2.$$
令$\frac{\partial F}{\partial x_1}=0,\frac{\partial F}{\partial x_2}=0$, 得
$$x_1=\frac{-M}{M+1},x_2=0.$$
类似地, 有$M>0\Rightarrow {\nabla }^{2}F$正定. 因此, $F(\mathbf{x},M)$在$(\frac{-M}{M+1},0{)}^T$取得极小值. 令$M\to +\infty$, 则原约束问题的最优解为
$$x^{\ast }=\underset{M\to +\infty }{\lim }\mathbf{x}(M)=(-1,0{)}^T.$$
且
$$F(\mathbf{x},M)=(-\frac{M}{M+1}{)}^2+M(\frac{1}{M+1}{)}^2=\frac{M}{M+1}\to f(x^{\ast })=1.$$
Remark: 虽然$\mathbf{x}(M)\to x^{\ast }(M\to +\infty )$, 但$\mathbf{x}(M)$往往不满足约束条件, 如在例1中$x_1(M)+x_2(M)=\frac{4M}{2M+1}\ne 1$, 在例2中$x_1(M)=-\frac{M}{M+1}>-1$, 且$\mathbf{x}(M)$都从可行域外部趋向于$x^{\ast }$, 因此惩罚函数法也称为外罚函数法.

算法流程

已知约束问题, 取控制误差$\epsilon >0$和罚因子放大系数$c>1(e.g.\epsilon =10^{-4},c=10)$.

步骤1 给定初始点$x_0$(可以不在可行域内)和初始惩罚因子$M_1=1$, 令$k=1$.

步骤2 以$x_{k-1}$为初始点求无约束问题
$$\min F(\mathbf{x},M_k)=f(\mathbf{x})+M_{k}p(x)$$
得最优解${\mathbf{x}}_k=\mathbf{x}(M_k)$.

步骤3 若$M_{k}p({\mathbf{x}}_{\mathbf{k}})<\epsilon$, 则以${\mathbf{x}}_k$为近似最优解, 停止; 否则令$M_{k+1}=c{M}_k,k=k+1$, 转步骤2.

计算实例

求解约束问题
$$\begin{array}{rl}& \min f(\mathbf{x})=(x_1-2{)}^4+(x_2-2x_1{)}^2\\ & s.t.x_1^2-x_2=0.\end{array}$$
易知该问题的最优解为$(2,4{)}^T$.

取${\mathbf{x}}_0=(2,1{)}^T,M_1=0.1,c=10,\epsilon =0.01,\beta =2,$求解
$$\min F(\mathbf{x},M)=(x_1-2{)}^4+(x_2-2x_1{)}^2+M_k(x_1^2-x_2{)}^2.$$
经过两次迭代, $M_{2}p({\mathbf{x}}_2)=1.805074\times 10^{-10}<\epsilon$, 算法停止. 近似最优解为$x^{\ast }=(1.999998,4.000007{)}^T$.

Code(简易版)

% 主程序
function SUMT
global M
x0=[2,1]';
M=0.1;
c=10;
eps=0.01;
k=1;
while M*p(x0)>=eps
    x0=fminsearch('Aim',x0);
    M=c*M;
    k=k+1;
end
format long
disp('最优解:')
disp(x0)
disp('迭代次数')
disp(k)
disp('终止条件')
disp(M*p(x0))
end

function r=p(x)
% 罚函数
r=(x(1)^2-x(2))^2;
end

function r=Aim(x)
global M
r=(x(1)-2)^4+(x(2)-2*x(1))^2+M*p(x);
end

算法缺点及改进

1. 每个近似最优解$\mathbf{x}(M_k)$不是可行解. 改进: 内罚函数法.
2. 由收敛性知, $M_k$越大越好, 而这会造成$F(\mathbf{x},M)$的Hessen矩阵条件数越大, 趋向于病态, 给无约束问题求解增加困难, 甚至于无法求解. 改进: 乘子法.

Reference

@book{唐焕文2000实用最优化方法,
  title={实用最优化方法.第2版},
  author={唐焕文 and 秦学志},
  publisher={实用最优化方法.第2版},
  year={2000},
}
@book{肖柳青2000实用最优化方法,
  title={实用最优化方法},
  author={肖柳青 and 周石鹏},
  publisher={实用最优化方法},
  year={2000},
}