Python Tutorials

Part 1. 程式基礎與數學篇

Programming Basics and Math Python

 

Section 3. 極值問題

 

Section 3. 極值問題

最佳化演算(Optimization) 求極值


前面的主題是解滿足一階導函數的根,因此,我們必須先微分求出導函數。但是,最簡單的方法是直接對目標函數求解:同時解出目標極值和臨界值。
Python完成這件事有sympy和scipy,筆者覺得sympy需要宣告的參數太多,尤其在帶限制式時。所以我們使用模組scipy內的函數optimize()和minimize(),Python程式碼的步驟解說如下。

3-1 單變數

我們先看本範例,也就是主題一的簡單範例

第1步:定義函數
from scipy import optimize
def f(x, sign=-1):
        return sign*(2*x**3+3*x**2-12*x-7)

兩行就OK,相當簡易。待會我們再解釋sign的意義。皆下來執行求極值:

第2步:求解與結果
Result1 = optimize. minimize_scalar(f)
Result1.x
Result1.fun
f(Result1.x)

optimize. minimize_scalar()是求解函數。Result1內有許多物件,主要有三個:

我們看看列印在螢幕的結果
Result1.x
Out[2]: 1.0

Result1.fun
Out[3]: -14.0

f(Result1.x)
Out[4]: -14.0

我們回去看本範圍範例的此題,可能的臨界值有兩個,從圖形看的出來,我們解出的只是極小值 (1, -14)。那另一極大值的解呢?根據scipy說明文件,須把函數取負值 ,這也是我們為什麼寫函數時,要增加一個參數 sign,因為這樣比較方便,判斷極大值時,可以如下這樣處理
第1步:定義函數
def f(x, sign = -1):
        return sign*(2*x**3+3*x**2-12*x-7)

第2步:求解與結果
Result2 = optimize.minimize_scalar(f)

螢幕的結果如下(需注意極值須加負號)
Result2.x
Out[5]: -1.999999999777818

-Result2.fun
Out[6]: 13.0

-f(Result2.x)
Out[7]: 13.0

本範圍範例的函數f(x)=x4-x3相對極小值
第1步:定義函數
def f(x,sign=1):
        return sign*(x**4-x**3)

第2步:求解與結果
Result3 = optimize.minimize_scalar(f)

螢幕的結果如下
Result3.x
Out[9]: 0.7500000000447832

Result3.fun
Out[10]: -0.10546875

f(Result3.x)
Out[11]: -0.10546875

 

3-2 多變數
再來就是不帶限制條件的多變數函數,本範例的函數相對極小值,則

第1步:定義函數
import numpy as np
from scipy.optimize import minimize
def f(x, sign=1):
        x1 = x[0]
        x2 = x[1]
        return sign*(x1**3-4*x1*x2 +2*x2**2)
第2步:求解與結果
x0=[1,1]
Result4 = minimize(f, x0)

        雙變數以上的數值求解演算比較複雜,我們使用的函數是minimize(),如果用上面的optimize.minimize_scalar()執行會失敗。x0=[1,1]是起始值(initial values)

螢幕的結果如下
Result4.x
Out[13]: array([1.33333404, 1.3333353 ])

Result4.fun
Out[14]: -1.185185185181036

f(Result4.x)
Out[15]: -1.185185185181036

 

因為是數值結果,手解的臨界值是4/3,電腦則算出1.3333。相對極小值則是 -1.185。

此題還有一解(0,0)是鞍點,如果設定x0=[0,0],就會帶出0的極值。因為鞍點判斷的程式作法需要賦予更多的條件,不是本書涵蓋,也不是商學院微積分主題,我們大致知道目前學習的狀況即可。如果想挑戰Python 程式,把目前所學過的方法串起來,可以循以下步驟:
步驟 1. 使用diff()函數解一階偏微分,求取可能的臨界值
步驟 2. 求二階偏微分
步驟3. 參考主題一的二元一次方程式求解判斷的準則(Delta),定義Hessian 來判斷誰是鞍點,誰無解。
最後. 把可能的臨界值設為起始值,求解。

        這樣的四步驟其實是一個程式訓練的典型基礎,有程式興趣的同學可以用在本範例,因為避免微積分學習花太多時間講這個,我們就到此為止。

3-3 多變數帶限制式

        最後就是帶限制條件的極值,我們以本範圍範例來說明


第1步:定義目標函數
import numpy as np
from scipy.optimize import minimize
def f(x, sign=1):
        x1 = x[0]
        x2 = x[1]
        return sign*(x1+ 2*x2)

第2步:定義限制條件
def constraint1(x, sign=1):
        return sign*(x[0]*x[1]- 5000)

第3步:設定求解參數
        x0=[10,10]    # 起始值,
        b1 = (0, np.inf)  # 參數條件x>0,上界給予正無限大,np.inf就是¥
        b2 = (0, np.inf)  # 參數條件y>0,上界給予正無限大,np.inf就是¥
        bnds= (b1,b2)  # 邊界條件向量
        con1 = {'type': 'ineq', 'fun': constraint1} #把限制集定義成字典
        cons = [con1]  #把con1做成串列 (萬一有多個條件時,可以包在一起)

第4步:求解與結果
Result5 = minimize(f, x0,bounds=bnds,constraints=cons)
螢幕的結果如下

Result5.x
Out[17]: array([100.00691556,  49.99654246])

Result5.fun
Out[18]: 200.0000004780455

f(Result5.x)
Out[19]: 200.0000004780455

 

可對照數值的結果和代數的結果。

接下來我們看「三變數,兩條限制式」,如本範例。這個範例足以做很多的推廣


 

第1步:定義目標函數
import numpy as np
from scipy.optimize import minimize
def f(x, sign=1):
        x1 = x[0]
        x2 = x[1]
        x3 = x[2]
        return sign*(x1**2+ x2**2+x3**2)

第2步:定義限制條件
def constraint1(x, sign=1):
        return sign*(x[0]+x[1]-3)
        def constraint2(x, sign=1):
        return sign*(x[0]+x[2]- 5)

第3步:設定求解參數 (參數沒有bounds,故不需定義如前提的b1,b2)
x0=[1,1,1]
con1 = {'type': 'ineq', 'fun': constraint1}
con2 = {'type': 'ineq', 'fun': constraint2}
cons = [con1, con2]

第4步:求解與結果
Result6 = minimize(f,x0,constraints=cons)

螢幕的結果如下
Result6 = minimize(f,x0,constraints=cons)

Result6.x
Out[21]: array([2.6666667, 0.3333333, 2.3333333])

Result6.fun
Out[22]: 12.666666666666735

f(Result6.x)
Out[23]: 12.666666666666735

可以確認數值結果是一樣的。這樣本範圍所有的問題,我們都可以處理了