Python SciPy optimize.shgo用法及代碼示例

本文簡要介紹 python 語言中 scipy.optimize.shgo 的用法。

用法:

scipy.optimize.shgo(func, bounds, args=(), constraints=None, n=100, iters=1, callback=None, minimizer_kwargs=None, options=None, sampling_method='simplicial', *, workers=1)#

使用 SHG 優化查找函數的全局最小值。

SHGO 代表“單純同調全局優化”。

參數 ：：

func： 可調用的

要最小化的目標函數。必須采用 f(x, *args) 的形式，其中 x 是一維數組形式的參數，而 args 是完全指定函數所需的任何其他固定參數的元組。

bounds： 序列或 Bounds

變量的界限。有兩種方法可以指定邊界：

1. Bounds 類的實例。

2. 序列(min, max)中每個元素的對x.

args： 元組，可選

完全指定目標函數所需的任何其他固定參數。

constraints： {Constraint, dict} 或 {Constraint, dict} 列表，可選

約束定義。僅適用於 COBYLA、SLSQP 和trust-constr。有關指定約束的更多詳細信息，請參閱教程 [5]。

注意

當前僅 COBYLA、SLSQP 和 trust-constr 局部最小化方法支持約束參數。如果局部優化問題中使用的constraints序列未在minimizer_kwargs中定義，並且使用約束方法，則將使用全局constraints。 (在minimizer_kwargs中定義constraints序列意味著不會添加constraints，因此如果需要添加等式約束等，則constraints中的不等式函數也需要添加到minimizer_kwargs中)。 COBYLA 僅支持不平等約束。

n： 整數，可選

構造單純複形時使用的采樣點數量。對於默認的simplicial采樣方式 生成2**dim + 1個采樣點，而不是默認的n=100。對於所有其他指定值n生成采樣點。為了sobol,halton和其他任意sampling_methods n=100或者生成其他指定數量的采樣點。

iters： 整數，可選

用於構造單純複形的迭代次數。默認值為 1。

callback： 可調用的，可選的

每次迭代後調用，如 callback(xk) ，其中 xk 是當前參數向量。

minimizer_kwargs： 字典，可選

要傳遞給最小化器的額外關鍵字參數 scipy.optimize.minimize 一些重要的選項可能是：

• method str

  The minimization method. If not given, chosen to be one of BFGS, L-BFGS-B, SLSQP, depending on whether or not the problem has constraints or bounds.

• args tuple

  Extra arguments passed to the objective function (func) and its derivatives (Jacobian, Hessian).

• options dict, optional

  Note that by default the tolerance is specified as {ftol: 1e-12}

options： 字典，可選

求解器選項字典。為全局例程指定的許多選項也傳遞給scipy.optimize.minimize 例程。也傳遞給本地例程的選項標有“(L)”。

停止條件，如果滿足任何指定條件，算法將終止。但是，默認算法不需要指定任何內容：

• maxfev int (L)

  可行域中函數評估的最大數量。 (請注意，隻有支持此選項的方法才會在精確指定的值處終止例程。否則，標準將僅在全局迭代期間終止)

• f_min

  指定最小目標函數值(如果已知)。

• f_tol float

  停止準則中 f 值的精度目標。請注意，如果全局例程中的采樣點在此容差範圍內，全局例程也將終止。

• maxiter int

  要執行的最大迭代次數。

• maxev int

  要執行的最大抽樣評估次數(包括在不可行點中搜索)。

• maxtime float

  允許的最大處理運行時間

• minhgrd int

  最小同源組等級差異。在每次迭代期間(近似地)計算目標函數的同調群。該組的排名與目標函數中局部凸子域的數量一一對應(在足夠的采樣點之後，每個子域都包含唯一的全局最小值)。如果 maxhgrd 指定迭代之間的 hgr 差異為 0，則算法將終止。

目標函數知識：

• symmetry list or bool

  指定目標函數是否包含對稱變量。在完全對稱的情況下，搜索空間(以及因此的性能)最多減少 O(n!) 倍。如果指定 True，則所有變量將設置為與第一個變量對稱。默認設置為 False。

  例如： f(x) = (x_1 + x_2 + x_3) + (x_4)**2 + (x_5)**2 + (x_6)**2

  在此方程中，x_2 和 x_3 與 x_1 對稱，而 x_5 和 x_6 與 x_4 對稱，這可以向求解器指定為：

  對稱性 = [0, # 變量 1

  0, # 與變量 1 對稱 0, # 與變量 1 對稱 3, # 變量 4 3, # 與變量 4 對稱 3, # 與變量 4 對稱 ]

• jac bool or callable, optional

  目標函數的雅可比(梯度)。僅適用於 CG、BFGS、Newton-CG、L-BFGS-B、TNC、SLSQP、dogleg、trust-ncg。如果jac是一個布爾值並且是 True，fun假設與目標函數一起返回梯度。如果為 False，梯度將以數值方式估計。jac也可以是返回目標梯度的可調用對象。在這種情況下，它必須接受與fun. (傳遞給scipy.optimize.minmize自動地)

• hess, hessp callable, optional

  目標函數的 Hessian(二階導數矩陣)或目標函數的 Hessian 乘以任意向量 p。僅適用於Newton-CG、狗腿、trust-ncg。僅其中之一hessp或者hess需要給予。如果hess提供，然後hessp將被忽略。如果兩者都沒有hess也不hessp提供，則 Hessian 乘積將使用有限差分來近似jac.hessp必須計算 Hessian 乘以任意向量。 (傳遞給scipy.optimize.minmize自動地)

算法設置：

• minimize_every_iter bool

  如果為 True，那麽有希望的全局采樣點將在每次迭代時傳遞給局部最小化例程。如果為 True，則僅運行最終的最小化器池。默認為 True。

• local_iter int

  每次迭代隻評估幾個最佳的最小化池候選者。如果為 False，則所有潛在點都將傳遞給局部最小化例程。

• infty_constraints bool

  如果為真，那麽在可行域之外生成的任何采樣點都將被保存，並賦予目標函數值 inf 。如果為 False，那麽這些點將被丟棄。在找到全局最小值之前，使用此函數可能會導致函數評估的性能更高，指定 False 將使用更少的內存，但性能會略有下降。默認為真。

回饋：

• disp bool (L)

  設置為 True 以打印收斂消息。

sampling_method： str 或函數，可選

當前內置的采樣方法選項是halton,sobol和simplicial.默認simplicial提供了在有限時間內收斂到全局最小值的理論保證。halton和sobol方法在采樣點生成方麵更快，但以損失保證收斂為代價。它更適合大多數收斂速度相對較快的“easier” 問題。用戶定義的采樣函數必須接受兩個參數n維度采樣點dim每次調用並輸出具有形狀的采樣點數組n x dim.

workers： int 或 map-like 可調用，可選

並行采樣並運行本地串行最小化。提供 -1 以使用所有可用的 CPU 內核，或提供 int 以使用多個進程(使用 multiprocessing.Pool )。

或者提供 map-like 可調用，例如多處理.Pool.map用於並行評估。本次評價按如下方式進行workers(func, iterable).要求函數可以 pickle 。

返回 ：：

res： OptimizeResult

優化結果表示為 OptimizeResult 對象。重要的屬性是：x 對應於全局最小值的解數組，fun 全局解處的函數輸出，xl 局部最小值解的有序列表，funl 對應局部解處的函數輸出， success 一個布爾標誌，指示優化器是否成功退出，message 說明終止的原因，nfev 包括采樣調用在內的目標函數評估總數，nlfev 目標函數評估總數從所有本地搜索優化中，nit 全局例程執行的迭代次數。

注意：

使用單純同調全局優化的全局優化[1]。適用於解決通用 NLP 和黑盒優化問題到全局最優(低維問題)。

一般來說，優化問題的形式如下：

minimize f(x) subject to

g_i(x) >= 0,  i = 1,...,m
h_j(x)  = 0,  j = 1,...,p

其中 x 是一個或多個變量的向量。 f(x)是目標函數R^n -> R，g_i(x)是不等式約束，h_j(x)是等式約束。

或者，也可以使用 bounds 參數指定 x 中每個元素的下限和上限。

雖然 SHGO 的大多數理論優勢僅在 f(x) 是 Lipschitz 平滑函數時得到證明，但該算法也被證明在 f(x) 是非連續、非凸的更一般情況下收斂到全局最優如果使用默認采樣方法 [1]，則不平滑。

可以使用傳遞給 scipy.optimize.minimize 的 minimizer_kwargs 參數指定本地搜索方法。默認情況下，使用SLSQP 方法。通常，如果為問題定義了不等式約束，則建議使用SLSQP 或COBYLA 局部最小化，因為其他方法不使用約束。

halton 和 sobol 方法點是使用 scipy.stats.qmc 生成的。可以使用任何其他 QMC 方法。

參考：

[1] (1,2)

Endres, SC, Sandrock, C, Focke, WW (2018) “一種簡單的lipschitz 優化同源算法”，全局優化雜誌。

[2]

Joe, SW 和 Kuo, FY (2008) “使用更好的二維投影構建 Sobol 序列”，SIAM J. Sci。計算。 30, 2635-2654。

[3] (1,2)

Hock, W 和 Schittkowski, K (1981)“非線性編程代碼的測試示例”，經濟學和數學係統講義，187。Springer-Verlag，紐約。http://www.ai7.uni-bayreuth.de/test_problem_coll.pdf

[4]

威爾士，DJ (2015) “透視：從勢能景觀中洞察反應坐標和動力學”，化學物理學雜誌，142(13)，2015。

[5]

https://docs.scipy.org/doc/scipy/tutorial/optimize.html#constrained-minimization-of-multivariate-scalar-functions-minimize

例子：

首先考慮最小化 Rosenbrock 函數的問題， rosen ：

>>> from scipy.optimize import rosen, shgo
>>> bounds = [(0,2), (0, 2), (0, 2), (0, 2), (0, 2)]
>>> result = shgo(rosen, bounds)
>>> result.x, result.fun
(array([1., 1., 1., 1., 1.]), 2.920392374190081e-18)

請注意，邊界確定目標函數的維度，因此是必需的輸入，但是您可以使用 None 或將轉換為大浮點數的對象(如 np.inf)指定空邊界。

>>> bounds = [(None, None), ]*4
>>> result = shgo(rosen, bounds)
>>> result.x
array([0.99999851, 0.99999704, 0.99999411, 0.9999882 ])

接下來，我們考慮 Eggholder 函數，這是一個具有多個局部最小值和一個全局最小值的問題。我們將演示參數的使用和 shgo 的函數。 (https://en.wikipedia.org/wiki/Test_functions_for_optimization)

>>> import numpy as np
>>> def eggholder(x):
...     return (-(x[1] + 47.0)
...             * np.sin(np.sqrt(abs(x[0]/2.0 + (x[1] + 47.0))))
...             - x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47.0))))
...             )
...
>>> bounds = [(-512, 512), (-512, 512)]

shgo具有內置的低差異采樣序列。首先，我們將輸入 64 個初始采樣點索博爾順序：

>>> result = shgo(eggholder, bounds, n=64, sampling_method='sobol')
>>> result.x, result.fun
(array([512.        , 404.23180824]), -959.6406627208397)

shgo 還返回找到的任何其他局部最小值，可以使用以下命令調用：

>>> result.xl
array([[ 512.        ,  404.23180824],
       [ 283.0759062 , -487.12565635],
       [-294.66820039, -462.01964031],
       [-105.87688911,  423.15323845],
       [-242.97926   ,  274.38030925],
       [-506.25823477,    6.3131022 ],
       [-408.71980731, -156.10116949],
       [ 150.23207937,  301.31376595],
       [  91.00920901, -391.283763  ],
       [ 202.89662724, -269.38043241],
       [ 361.66623976, -106.96493868],
       [-219.40612786, -244.06020508]])

>>> result.funl
array([-959.64066272, -718.16745962, -704.80659592, -565.99778097,
       -559.78685655, -557.36868733, -507.87385942, -493.9605115 ,
       -426.48799655, -421.15571437, -419.31194957, -410.98477763])

這些結果在存在許多全局最小值並且需要其他全局最小值的值或局部最小值可以提供對係統的洞察力的應用中很有用(例如物理化學中的形態[4])。

如果我們想找到更多的局部最小值，我們可以增加采樣點的數量或迭代次數。我們將采樣點的數量增加到 64 個，並將迭代次數從默認值 1 增加到 3。使用 simplicial 這將為我們提供 64 x 3 = 192 個初始采樣點。

>>> result_2 = shgo(eggholder,
...                 bounds, n=64, iters=3, sampling_method='sobol')
>>> len(result.xl), len(result_2.xl)
(12, 23)

請注意，例如 n=192, iters=1 和 n=64, iters=3 之間的區別。在第一種情況下，minimiser 池中包含的有希望的點隻被處理一次。在後一種情況下，每 64 個采樣點處理一次，總共處理 3 次。

為了演示解決具有非線性約束的問題，請考慮以下來自 Hock 和 Schittkowski 問題 73 (cattle-feed) [3] 的示例：

minimize: f = 24.55 * x_1 + 26.75 * x_2 + 39 * x_3 + 40.50 * x_4

subject to: 2.3 * x_1 + 5.6 * x_2 + 11.1 * x_3 + 1.3 * x_4 - 5    >= 0,

            12 * x_1 + 11.9 * x_2 + 41.8 * x_3 + 52.1 * x_4 - 21
                -1.645 * sqrt(0.28 * x_1**2 + 0.19 * x_2**2 +
                              20.5 * x_3**2 + 0.62 * x_4**2)      >= 0,

            x_1 + x_2 + x_3 + x_4 - 1                             == 0,

            1 >= x_i >= 0 for all i

[3] 中給出的近似答案是：

f([0.6355216, -0.12e-11, 0.3127019, 0.05177655]) = 29.894378

>>> def f(x):  # (cattle-feed)
...     return 24.55*x[0] + 26.75*x[1] + 39*x[2] + 40.50*x[3]
...
>>> def g1(x):
...     return 2.3*x[0] + 5.6*x[1] + 11.1*x[2] + 1.3*x[3] - 5  # >=0
...
>>> def g2(x):
...     return (12*x[0] + 11.9*x[1] +41.8*x[2] + 52.1*x[3] - 21
...             - 1.645 * np.sqrt(0.28*x[0]**2 + 0.19*x[1]**2
...                             + 20.5*x[2]**2 + 0.62*x[3]**2)
...             ) # >=0
...
>>> def h1(x):
...     return x[0] + x[1] + x[2] + x[3] - 1  # == 0
...
>>> cons = ({'type': 'ineq', 'fun': g1},
...         {'type': 'ineq', 'fun': g2},
...         {'type': 'eq', 'fun': h1})
>>> bounds = [(0, 1.0),]*4
>>> res = shgo(f, bounds, n=150, constraints=cons)
>>> res
 message: Optimization terminated successfully.
 success: True
     fun: 29.894378159142136
    funl: [ 2.989e+01]
       x: [ 6.355e-01  1.137e-13  3.127e-01  5.178e-02]
      xl: [[ 6.355e-01  1.137e-13  3.127e-01  5.178e-02]]
     nit: 1
    nfev: 142
   nlfev: 35
   nljev: 5
   nlhev: 0

>>> g1(res.x), g2(res.x), h1(res.x)
(-5.062616992290714e-14, -2.9594104944408173e-12, 0.0)