罚函数法又称外点法,其基本思想为将违背约束作为求最小值的一种惩罚,将约束带入目标函数得到一个辅助的无约束最优化问题。利用已有的无约束最优化方法求解。
考虑如下约束优化问题:
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
x
x位于可行域外时,
F
F
F取值很大,且离可行域越远
F
F
F越大;
b. 当点
x
x
x位于可行域内时,
F
=
f
F=f
F=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\geq0
y≥0,
Φ
(
y
)
=
0
\Phi(y)=0
Φ(y)=0;当
y
<
0
y<0
y<0,
Φ
(
y
)
>
0
\Phi(y)>0
Φ(y)>0;
当
y
=
0
y=0
y=0,
Ψ
(
y
)
=
0
\Psi(y)=0
Ψ(y)=0;当
y
≠
0
y \neq0
y=0,
Ψ
(
y
)
>
0
\Psi(y)>0
Ψ(y)>0;
即,当点
x
x
x在可行域内,
F
=
f
F=f
F=f,在可行域外,
F
F
F越来越大。
Φ
,
Ψ
\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\geq1
a≥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 0
minxs.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