🥇SciPy，條件最佳化 |專業主機

SciPy（發音為 sai pie）是一個基於 numpy 的數學包，還包括 C 和 Fortran 庫。 SciPy 將您的互動式 Python 會話轉變為完整的資料科學環境，例如 MATLAB、IDL、Octave、R 或 SciLab。

在本文中，我們將了解數學規劃的基本技術 - 使用 scipy.optimize 套件解決多個變數的標量函數的條件最佳化問題。無約束優化演算法已經在上一篇文章。有關 scipy 函數的更詳細和最新的幫助始終可以使用 help() 命令、Shift+Tab 或官方文檔.

介紹

函數提供了一個通用接口，用於解決 scipy.optimize 包中的條件和無約束優化問題 minimize()。然而，眾所周知，沒有解決所有問題的通用方法，因此選擇適當的方法一如既往地落在研究人員的肩上。
使用函數參數指定適當的最佳化演算法 minimize(..., method="").
對於多個變數的函數的條件最佳化，可以使用以下方法的實作：

trust-constr — 在置信區域中搜尋局部最小值。維基文章, 關於哈布雷的文章;
SLSQP — 約束的順序二次規劃，求解拉格朗日系統的牛頓法。維基文章.
TNC - 截斷牛頓約束，迭代次數有限，適合具有大量自變數的非線性函數。維基文章.
L-BFGS-B — Broyden–Fletcher–Goldfarb–Shanno 團隊的一種方法，由於從 Hessian 矩陣中部分載入向量，因此減少了記憶體消耗。維基文章, 關於哈布雷的文章.
COBYLA — MARE Constrained Optimization By Linear Approximation，使用線性近似的約束優化（無梯度計算）。維基文章.

根據所選的方法，解決問題的條件和限制設定不同：

類對象 Bounds 對於方法 L-BFGS-B、TNC、SLSQP、trust-constr；
列表 (min, max) 對於相同的方法 L-BFGS-B、TNC、SLSQP、trust-constr；
一個物件或物件列表 LinearConstraint, NonlinearConstraint 對於 COBYLA、SLSQP、信任構造方法；
字典或字典列表 {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} 對於 COBYLA、SLSQP 方法。

文章大綱：
1) 考慮在信任區域（method=”trust-constr”）中使用條件最佳化演算法，並將約束指定為對象 Bounds, LinearConstraint, NonlinearConstraint ;
2) 考慮使用最小平方法（method =“SLSQP”）進行順序編程，並以字典的形式指定限制 {'type', 'fun', 'jac', 'args'};
3) 使用網路工作室的範例來分析製造產品最佳化的範例。

條件最佳化方法=“trust-constr”

方法的實施 trust-constr 基於 EQSQP 對於具有平等形式和約束條件的問題 TRIP 針對不平等形式的約束問題。這兩種方法都是透過在置信區域中尋找局部最小值的演算法來實現的，並且非常適合大規模問題。

尋找最小值問題的一般形式的數學表述：

對於嚴格的相等約束，下限設定為等於上限 .
對於單向約束，設定上限或下限 np.inf 並帶有相應的標誌。
需要找到兩個變數的已知 Rosenbrock 函數的最小值：

在這種情況下，對其定義域設定以下限制：

在我們的例子中，此時有一個獨特的解決方案，僅第一和第四限制有效。
讓我們從下到上看一下這些限制，看看如何在 scipy 中編寫它們。
限制 и 讓我們使用 Bounds 物件來定義它。

from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])

限制 и 我們把它寫成線性形式：

讓我們將這些約束定義為 LinearConstraint 物件：

import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

最後是矩陣形式的非線性約束：

我們為此約束定義雅可比矩陣以及 Hessian 矩陣與任意向量的線性組合 :

現在我們可以將非線性約束定義為對象 NonlinearConstraint:

from scipy.optimize import NonlinearConstraint

def cons_f(x):
     return [x[0]**2 + x[1], x[0]**2 - x[1]]

def cons_J(x):
     return [[2*x[0], 1], [2*x[0], -1]]

def cons_H(x, v):
     return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)

如果尺寸很大，矩陣也可以以稀疏形式指定：

from scipy.sparse import csc_matrix

def cons_H_sparse(x, v):
     return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
                                            jac=cons_J, hess=cons_H_sparse)

或作為一個對象 LinearOperator:

from scipy.sparse.linalg import LinearOperator

def cons_H_linear_operator(x, v):
    def matvec(p):
        return np.array([p[0]*2*(v[0]+v[1]), 0])
    return LinearOperator((2, 2), matvec=matvec)

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
                                jac=cons_J, hess=cons_H_linear_operator)

計算Hessian矩陣時需要很多努力，你可以使用一個類 HessianUpdateStrategy。可以採用以下策略： BFGS и SR1.

from scipy.optimize import BFGS

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())

Hessian 矩陣也可以使用有限差分來計算：

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')

約束的雅可比矩陣也可以使用有限差分來計算。然而，在這種情況下，Hessian 矩陣無法使用有限差分來計算。 Hessian 必須定義為函數或使用 HessianUpdateStrategy 類別。

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())

最佳化問題的解決方案如下所示：

from scipy.optimize import minimize
from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prod

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.033 s.
[0.41494531 0.17010937]

如果需要，可以使用 LinearOperator 類別定義計算 Hessian 的函數

def rosen_hess_linop(x):
    def matvec(p):
        return rosen_hess_prod(x, p)
    return LinearOperator((2, 2), matvec=matvec)

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
                 constraints=[linear_constraint, nonlinear_constraint],
                 options={'verbose': 1}, bounds=bounds)

print(res.x)

或 Hessian 矩陣與通過參數的任意向量的乘積 hessp:

res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_prod,
                constraints=[linear_constraint, nonlinear_constraint],
                options={'verbose': 1}, bounds=bounds)
print(res.x)

或者，可以近似正在最佳化的函數的一階和二階導數。例如，Hessian 矩陣可以使用下列函數進行近似 SR1 （準牛頓近似）。梯度可以用有限差分來近似。

from scipy.optimize import SR1
res = minimize(rosen, x0, method='trust-constr',  jac="2-point", hess=SR1(),
               constraints=[linear_constraint, nonlinear_constraint],
               options={'verbose': 1}, bounds=bounds)
print(res.x)

條件最佳化方法=“SLSQP”

SLSQP 方法旨在解決以下形式的函數最小化問題：

在哪裡 и ——以等式或不等式形式描述限制條件的表達式索引集。 — 函數定義域的下限和上限集合。

線性和非線性約束以帶有鍵的字典的形式描述 type, fun и jac.

ineq_cons = {'type': 'ineq',
             'fun': lambda x: np.array ([1 - x [0] - 2 * x [1],
                                          1 - x [0] ** 2 - x [1],
                                          1 - x [0] ** 2 + x [1]]),
             'jac': lambda x: np.array ([[- 1.0, -2.0],
                                          [-2 * x [0], -1.0],
                                          [-2 * x [0], 1.0]])
            }

eq_cons = {'type': 'eq',
           'fun': lambda x: np.array ([2 * x [0] + x [1] - 1]),
           'jac': lambda x: np.array ([2.0, 1.0])
          }

搜尋最小值的過程如下：

x0 = np.array([0.5, 0])
res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
               constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
               bounds=bounds)

print(res.x)

Optimization terminated successfully.    (Exit mode 0)
            Current function value: 0.34271757499419825
            Iterations: 4
            Function evaluations: 5
            Gradient evaluations: 4
[0.41494475 0.1701105 ]

優化示例

關於向第五種技術結構的過渡，讓我們以網路工作室為例來看看生產優化，它為我們帶來了少量但穩定的收入。讓我們想像自己是一家生產三種產品的廚房的主管：

x0 - 銷售登陸頁，從 10 tr 起。
x1 - 公司網站，20 tr 起。
x2 - 線上商店，30 tr 起。

我們友善的工作團隊包括四名初級人員、兩名中級人員和一名高級人員。他們每個月的工作時間基金：

六月： 4 * 150 = 600 чел * час,
中音： 2 * 150 = 300 чел * час,
先生： 150 чел * час.

讓第一個可用的初級人員花費(0, 1, 2) 小時來開發和部署類型為(x10, x20, x30) 的一個站點，中級人員- (7, 15, 20)，高級人員- (5 , 10, 15) ) 一生中最美好的時光。

像任何普通董事一樣，我們希望每月利潤最大化。成功的第一步是寫下目標函數 value 每月生產產品的收入金額：

def value(x):
    return - 10*x[0] - 20*x[1] - 30*x[2]

這不是錯誤；當搜尋最大值時，目標函數會以相反的符號最小化。

下一步是禁止我們的員工過度工作並限制工作時間：

什麼是等價的：

ineq_cons = {'type': 'ineq',
             'fun': lambda x: np.array ([600 - 10 * x [0] - 20 * x [1] - 30 * x[2],
                                         300 - 7  * x [0] - 15 * x [1] - 20 * x[2],
                                         150 - 5  * x [0] - 10 * x [1] - 15 * x[2]])
            }

正式的限制是產品輸出只能為正：

bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])

最後，最樂觀的假設是，由於價格低廉，品質高，滿意的客戶不斷為我們排隊。我們可以在解決約束優化問題的基礎上自行選擇每月的產量 scipy.optimize:

x0 = np.array([10, 10, 10])
res = minimize(value, x0, method='SLSQP', constraints=ineq_cons, bounds=bnds)
print(res.x)

[7.85714286 5.71428571 3.57142857]

讓我們鬆散地四捨五入為整數，並透過產品的最佳分佈來計算划船者的每月負荷 x = (8, 6, 3) :

六月： 8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час;
中音： 8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час;
先生： 8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.

結論：為了讓總監獲得他應得的最大獎勵，最好每月創建 8 個登陸頁面、6 個中型網站和 3 個商店。在這種情況下，高級人員必須在犁地時抬頭不離開機器，中級人員的負載約為2/3，初級人員不到一半。

結論

本文概述了使用該套件的基本技術 scipy.optimize，用於解決條件最小化問題。我個人使用 scipy 純粹出於學術目的，這就是為什麼給出的例子具有如此滑稽的性質。

許多理論和虛擬範例可以在 I.L. Akulich 所著的《範例與問題中的數學程式設計》一書中找到。更多硬派應用 scipy.optimize 從一組影像建立 3D 結構（關於哈布雷的文章) 可以在以下位置查看 scipy 食譜.

主要資訊來源是 docs.scipy.org希望為本節和其他部分的翻譯做出貢獻的人 scipy 歡迎來到 GitHub上.

謝謝梅菲斯托菲斯參與出版物的編寫。

來源： www.habr.com

SciPy，條件最佳化

介紹

條件最佳化方法=“trust-constr”

條件最佳化方法=“SLSQP”

優化示例

結論

添加評論 取消回复

添加評論取消回复