懲罰函數也叫乘子法，求解帶約束的非線性規划問題時，常用KKT條件列出滿足條件的方程組，解方程組后即可得到最值點，當滿足KKT條件的方程組是一個非線性方程組，利用計算機求解很難給出通用算法。本篇介紹的懲罰函數可以將一個帶約束非線性問題轉化為無約束的非線性規划，而無約束線性規划可以用梯度法等實現求解，利用懲罰函數更方便我們制成計算機算法，在現代計算機算法中，凡涉及到求解最值，都會大量的運用懲罰函數或者借鑒懲罰函數思想。
懲罰函數可以分為外點法和內點法，其中外點法更通用，可解決約束為等式和不等式混合的情形，外點法對初始點也沒有要求，可以任意取定義域內任意一點，外點法是真正的將問題轉化為無約束問題；而內點法初始點必須為可行區內一點，在約束比較復雜時，這個選擇內點法的初始點是有難度的，並且內點法只能解決約束為不等式情形。
一、 外點法
首先引入一個帶等式約束的非線性規划問題：
min f(x)
s.t. hi(x)=0 i=1,2,3,...m
目標函數f(x)是非線性函數，約束條件為m個等式，上面問題中可以通過拉格朗日乘數法或KKT條件獲得最值，而懲罰函數的解決方案是引入一個新的函數:
⑴
上式中F(x,σ)為懲罰函數，σ稱為懲罰因子，
稱為懲罰項，F(x,σ)中參數x沒有限制，可以取任意值。當x的取值在可行區內時，即x滿足所有等式時，F(x,σ)=f(x);而x在可行區外時，由於懲罰項是平方和形式，這時F(x,σ)>f(x)。綜上所述，懲罰函數F(x,σ)函數值大於等於f(x),或者說F(x,σ)是目標函數f(x)的上界函數,兩者的關系可以通過下圖說明：
x'點是原問題最小值點，即x'也是可行點，滿足原問題所有等式條件，F(x,σ)與f(x)在該點函數值相等：F(x',σ)=f(x'),x'點是目標最值點。通過觀察可以發現，雖然F(x,σ)是f(x)的上界函數，但最值F(x*)是小於f(x')的，這是由於F(x*)是F(x,σ)是在無約束情形下的最小值，F(x,σ)參數x可取值范圍比原帶約束的問題大的多（原問題x必須滿足所有等式），極點x*是函數F(x,σ)全局最小值點,外點懲罰函數通過每輪增大參數σ，逐漸讓x*逼近x',通過不斷迭代計算，最終的x*便是近似的目標函數極值點x',這樣一來，就實現了利用求解無約束最值得到另一個有約束問題最值。
再來看約束為不等式的情形，如求以下函數f(x)的最小值：
min f(x)
s.t. gi(x)>=0 i=1,2,3,...n
約束為不等式時，懲罰函數是下面的形式：
②
當懲罰函數參數x在可行區內時,懲罰項
，此時F(x,σ)=f(x)；當x在可行區外時F(x,σ)>f(x)，與約束為等式時一樣，F(x,σ)的最值F(x*)在初始階段是小於f(x')的，可以用等值線來解釋這種特性：
g(x)>=0為不等式約束，F(x,σ)函數極值點x*初始時在可行區外，點x*往g(x)=0方向移動時，F(x*,σ)函數值逐漸增大。在可行區內，F(x,σ)與f(x)等值線是一樣的，代表兩個函數在可行區內函數值相同，約束問題極值點x'一定是在可行區內，即在上圖虛線上以及虛線右側的范圍內；在虛線左側時F(x,σ)極值F(x*)小於f(x')。懲罰函數解決不等式約束時，是通過調整參數σ，讓x*逐漸接近可行區， 從而得到約束極值點x'。
通過上面兩種情形的分析，可以發現外點法是從可行區外慢慢接近邊界，在接近的過程中計算每個階段極值點，一旦到達可行區范圍內，該極值點即為原帶約束非線性規划的極值點，接下來還有最后一個問題，為什么增大參數值σ可以讓極值點x*接近可行區，以公式（2）為例，做適當變換可得：
σ增大時，懲罰項將逐漸等於0，代表F(x,σ)等值線整體逐漸向可行區移動。這里也可以看出，σ之所以叫懲罰因子是對懲罰函數的極值點不在可行區時，對懲罰函數進行相應的懲罰，因此，調整懲罰因子將移動懲罰函數極值點，這個動態的過程可參考下圖，懲罰函數Φ(x,rk)不斷增大rk后,Φ(x,rk)的極值點x*不斷接近目標函數f(x)極值點：
綜合上面兩種情形，當約束條件為不等式和等式混合時，有以下通用的懲罰函數形式：
③
可以看到，不等式與等式約束公用一個懲罰因子σ，σ>0。根據上面的分析，接下來通過一段python代碼實現外點懲罰函數，下面是一個含有兩個不等式、一個等式約束的非線性規划問題：
示例最優解可用以通過圖形觀察出來，目標函數最小值是求原點到可行區內最短距離的平方，等式3x+y-7.5=0與不等式約束所規定區域相交於A、B兩點，其中在點B處(2.25,0.75)為最小值點：
import math
import sys
sys.setrecursionlimit(500)
errorRate=1e-4
stepLowerlimit=1e-5
#計算懲罰項
def punishfunc(sigma, x, y):
punishvalue=(3 * x + y-7.5) ** 2
if (2 * x - y < 0 ):
punishvalue=punishvalue+ (2 * x - y) ** 2
if ( x + y - 3 < 0):
punishvalue = punishvalue + (x + y - 3) ** 2
return punishvalue*sigma
#懲罰函數一階導數
def P_firstderivative(stepx,x, y,gx,gy, sigma):
d1=2*(gx *stepx-x) *gx+2*(gy*stepx-y)*gy+2*sigma*(-3*gx-gy)*(3*x-3*gx*stepx+y-gy*stepx-7.5)
if (2 * x - y < 0):
d1=d1+2*sigma*(-2*gx+gy)*(2*x-2*gx*stepx-y+gy*stepx)
if(x + y - 3 < 0):
d1 = d1 + 2*sigma*(-gx-gy)*(x-gx*stepx+y-gy*stepx-3)
return d1
#懲罰函數二階導數
def P_secondderivative( x,y,gx,gy, sigma):
d2=2*gx**2+2*gy**2+2*sigma*(-3*gx-gy)**2
if (2 * x - y < 0):
d2=d2+2*sigma*(-2*gx+gy)**2
if (x + y - 3 < 0):
d2 = d2 + 2*sigma*(-gx-gy)**2
return d2
#牛頓法實現一維搜索
def newton(x, y,gx,gy, sigma):
stepx,l,iterNum=0,1e-4,0
dx_1=P_firstderivative(stepx,x, y,gx,gy, sigma )
while abs( dx_1)>l:
dx_1 = P_firstderivative(stepx, x, y, gx, gy, sigma)
dx_2 = P_secondderivative(x,y,gx, gy, sigma)
stepx=stepx-dx_1/dx_2
iterNum = iterNum + 1
return stepx
#求每次sigma對應的最小值,sigma每次遞增,逐步逼近目標函數最值
def calmin(x, y, sigma, counter):
gradpart1_x, gradpart1_y, gradpart2_x, gradpart2_y = 0, 0, 0, 0
if(2 * x - y<0 ):
gradpart1_x = 2 * (2 * x - y) * 2 * sigma
gradpart1_y = 2 * (2 * x - y) * -1 * sigma
if ( x + y - 3 < 0):
gradpart2_x = 2 * (x + y - 3) * sigma
gradpart2_y = 2 * (x + y - 3) * sigma
gradpart3_x = 2 * (3 * x + y - 7.5) * 3 * sigma
gradpart3_y = 2 * (3 * x + y - 7.5) * 1 * sigma
gradx = 2*x + gradpart1_x + gradpart2_x+gradpart3_x
grady = 2*y + gradpart1_y + gradpart2_y+gradpart3_y
distance = math.sqrt(gradx ** 2 + grady ** 2)
gradx = gradx / distance
grady = grady / distance
#牛頓法一維搜索得到步長
stepx=newton(x,y,gradx,grady,sigma)
x = x - gradx * stepx
y = y - grady * stepx
punishvalue = punishfunc(sigma, x, y)
value_end = x ** 2 + y ** 2 + punishvalue
if (stepx <= stepLowerlimit ):
return x, y, value_end, punishvalue, sigma
else:
counter=counter+1
return calmin(x, y, sigma, counter)
def outpoint(x, y, sigma, c ):
_x, _y, _v, _p=x,y,0,1
while ( _p>errorRate ):
_x, _y, _v, _p, _sigma = calmin(x, y, sigma, 0)
x=_x
y=_y
sigma = sigma * c
pass
value_end = x ** 2 + y ** 2 + punishfunc(sigma, x, y)
print('最小值=%0.2f x=%0.2f y=%0.2f'%(value_end,x,y))
if __name__ == '__main__':
sigma = 1#懲罰因子
c = 10
print('從外點開始:')
outpoint(0, 0, sigma, c )# 外點測試
print('-'*100)
print('從內點開始:')
outpoint(3 ,3, sigma, c )#內點測試
可以看到外點懲罰函數的初始點既可在可行區外，也可在可行區內，這兩種情形得到的結果是一致的，運行結果如下：
余下文章請轉至鏈接：