线性规划求解之罚函数法(外点法)

罚函数法又称外点法,其基本思想为将违背约束作为求最小值的一种惩罚,将约束带入目标函数得到一个辅助的无约束最优化问题。利用已有的无约束最优化方法求解。
考虑如下约束优化问题:
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=1mΦ(gi(x))+j=1lΨ(hj(x))
Φ , Ψ \Phi,\PsiΦ,Ψ是满足如下条件的连续函数:
y ≥ 0 y\geq0y0,Φ ( 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\geq1a1,b1,a,b为常数,通常取a=b=2。
综上,罚函数法求解线性规划问题步骤如下:

  1. 设置初始点x ( 0 ) x^{(0)}x(0),初始罚因子λ ( 1 ) \lambda^{(1)}λ(1), 放大系数c>1,允许误差e>0,k=1;
  2. x ( k − 1 ) x^{(k-1)}x(k1)作为搜索初始点,求解无约束规划问题:m i n f ( x ) + λ P ( x ) min f(x)+\lambda P(x)minf(x)+λP(x),令x ( k ) x^{(k)}x(k)为该无约束优化问题的极小值;
  3. λ P ( x ( k ) ) < e \lambda P(x^{(k)})<eλP(x(k))<e,停止计算,输出点x ( k ) x^{(k)}x(k);
  4. 否则,放大罚因子:令λ ( 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.2x0
该问题最优解为:当λ → ∞ , 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


版权声明:本文为qq_35356741原创文章,遵循CC 4.0 BY-SA版权协议,转载请附上原文出处链接和本声明。