罚函数法又称外点法,其基本思想为将违背约束作为求最小值的一种惩罚,将约束带入目标函数得到一个辅助的无约束最优化问题。利用已有的无约束最优化方法求解。
考虑如下约束优化问题:
m i n f ( x ) s . t . g i ( x ) ≥ 0 ( i = 1 , . . , m ) h j ( x ) = 0 ( j = 1 , . . . l ) minf(x)\\ s.t. g_i(x) \geq 0 (i=1,..,m)\\ h_j(x)=0 (j=1,...l)\\minf(x)s.t.gi(x)≥0(i=1,..,m)hj(x)=0(j=1,...l)
构造如下辅助函数:
F ( x , λ ) = f ( x ) + λ P ( x ) F(x, \lambda)=f(x)+\lambda P(x)F(x,λ)=f(x)+λP(x)
P ( x ) P(x)P(x) 称为罚函数。其具有以下性质:
a. 当点x xx位于可行域外时,F FF取值很大,且离可行域越远F FF越大;
b. 当点x xx位于可行域内时,F = f F=fF=f;
因此,上述约束优化问题可转化为无约束优化问题:
m i n F ( x , λ ) = f ( x ) + λ P ( x ) min F(x,\lambda)=f(x)+\lambda P(x)minF(x,λ)=f(x)+λP(x)
P ( x ) P(x)P(x)的定义一般如下:
P ( x ) = ∑ i = 1 m Φ ( g i ( x ) ) + ∑ j = 1 l Ψ ( h j ( x ) ) P(x)=\sum_{i=1}^{m} \Phi(g_i(x))+\sum_{j=1}^{l}\Psi(h_j(x))P(x)=i=1∑mΦ(gi(x))+j=1∑lΨ(hj(x))
Φ , Ψ \Phi,\PsiΦ,Ψ是满足如下条件的连续函数:
当y ≥ 0 y\geq0y≥0,Φ ( y ) = 0 \Phi(y)=0Φ(y)=0;当y < 0 y<0y<0,Φ ( y ) > 0 \Phi(y)>0Φ(y)>0;
当y = 0 y=0y=0,Ψ ( y ) = 0 \Psi(y)=0Ψ(y)=0;当y ≠ 0 y \neq0y=0,Ψ ( y ) > 0 \Psi(y)>0Ψ(y)>0;
即,当点x xx在可行域内,F = f F=fF=f,在可行域外,F FF越来越大。
Φ , Ψ \Phi,\PsiΦ,Ψ一般定义如下:
Φ = [ m a x 0 , − g i ( x ) ] a , Ψ = ∣ h j ( x ) ∣ b \Phi=[max{0,-g_i(x)}]^a,\Psi=|h_j(x)|^bΦ=[max0,−gi(x)]a,Ψ=∣hj(x)∣b
a ≥ 1 , b ≥ 1 a\geq1,b\geq1a≥1,b≥1,a,b为常数,通常取a=b=2。
综上,罚函数法求解线性规划问题步骤如下:
- 设置初始点x ( 0 ) x^{(0)}x(0),初始罚因子λ ( 1 ) \lambda^{(1)}λ(1), 放大系数c>1,允许误差e>0,k=1;
- 以x ( k − 1 ) x^{(k-1)}x(k−1)作为搜索初始点,求解无约束规划问题:m i n f ( x ) + λ P ( x ) min f(x)+\lambda P(x)minf(x)+λP(x),令x ( k ) x^{(k)}x(k)为该无约束优化问题的极小值;
- 当λ P ( x ( k ) ) < e \lambda P(x^{(k)})<eλP(x(k))<e,停止计算,输出点x ( k ) x^{(k)}x(k);
- 否则,放大罚因子:令λ ( k + 1 ) = c λ ( k ) \lambda ^{(k+1)}=c\lambda ^{(k)}λ(k+1)=cλ(k)
例1:
m i n x s . t . 2 − x ≤ 0 min x\\ s.t. 2-x \leq 0minxs.t.2−x≤0
该问题最优解为:当λ → ∞ , x ∗ = 2 \lambda \to \infty, x^*=2λ→∞,x∗=2。
罚函数求解过程参考代码如下:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin,fminbound
lamda=2
def f(x):
return x
def p(x):
return (2-x)**2
def f_p(x):
return x+lamda*p(x)
#start
x0=np.array([1])
c=10
e=1e-6
k=1
while p(x0)>=e:
x0=fmin(f_p,x0)
print(p(x0),x0,lamda)
lamda*=c
输出结果:
Optimization terminated successfully.
Current function value: 1.875000
Iterations: 16
Function evaluations: 32
[0.0625] [1.75] 2
Optimization terminated successfully.
Current function value: 1.987500
Iterations: 13
Function evaluations: 26
[0.00062561] [1.97498779] 20
Optimization terminated successfully.
Current function value: 1.998750
Iterations: 11
Function evaluations: 22
[6.46615472e-06] [1.99745714] 200
Optimization terminated successfully.
Current function value: 1.999880
Iterations: 11
Function evaluations: 22
[4.08417449e-08] [1.99979791] 2000
参考内容:
https://wenku.baidu.com/view/72d423ada45177232e60a2c6.html