👩🏽‍🔬 🤷🏻 ⏩ SciPy、条件付き最適化 🥌 🌏 🔶

SciPy（サイパイと発音）は、numpyベースの数学パッケージで、CおよびFortranライブラリも含まれています。 SciPyを使用すると、インタラクティブなPythonセッションは、MATLAB、IDL、Octave、R、SciLabなどの完全なデータ処理環境に変わります。

この記事では、scipy.optimizeパッケージを使用して、いくつかの変数のスカラー関数の条件付き最適化問題を解決する数学プログラミングの基本的な手法を検討します。無条件最適化アルゴリズムは、前回の記事ですでに検討されています。 scipy関数に関するより詳細で最新のヘルプは、常にhelp（）コマンド、Shift + Tab、または公式ドキュメントを使用して入手できます。

はじめに

scipy.optimizeパッケージの条件付き最適化と無条件最適化の両方の問題を解決するための一般的なインターフェースは、 minimize()関数によって提供されます。しかし、すべての問題を解決する普遍的な方法は存在しないことが知られているため、常に適切な方法を選択することは研究者次第です。

minimize(..., method="")関数の引数を使用して、適切な最適化アルゴリズムが指定されます。

いくつかの変数の関数の条件付き最適化には、次の方法が利用できます。

trust-constrトラストエリアのローカルミニマムを検索します。 Wiki 記事、 Habr記事、
SLSQP制約付き逐次二次計画法、ラグランジュ系を解くニュートン法。ウィキの記事。
TNCトランケートされたニュートン制約、反復回数の制限は、多数の独立変数を持つ非線形関数に適しています。ウィキの記事。
L-BFGS-B -4つのBroyden – Fletcher – Goldfarb – Shannoのメソッド。ヘッセ行列からのベクトルの部分的な読み込みによりメモリ消費を削減して実装されました。 wikiの記事、Habrの記事。
COBYLA - 馬線形近似による制約付き最適化、線形近似による制限付き最適化（勾配計算なし）。ウィキの記事。

選択した方法に応じて、問題を解決するための条件と制限が異なって設定されます。

メソッドL-BFGS-B、TNC、SLSQP、trust-constrのBoundsクラスのオブジェクト。
同じメソッドL-BFGS-B、TNC、SLSQP、trust-constrのリスト(min, max) 。
オブジェクトまたはオブジェクトのリストLinearConstraint 、メソッドCOBYLA、SLSQP、trust-constrのNonlinearConstraint 。
辞書または辞書リスト{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} 、COBYLA、SLSQPメソッド。

記事の概要：
1）オブジェクトBounds 、 LinearConstraint 、 NonlinearConstraint形式で指定された制限を持つ、信頼領域での条件付き最適化アルゴリズム（method = "trust-constr"）の適用を検討してください。
2）辞書{'type', 'fun', 'jac', 'args'}の形式で指定された制約を使用した連続最小二乗プログラミング（method = "SLSQP"）を検討します。
3）Webスタジオの例で製品の最適化の例を分解する。

条件付き最適化方法= "trust-constr"

trust-constrの実装は、等式の制約を持つ問題のEQSQPに基づいており、不等式の制約を持つ問題のTRIPに基づいています。どちらの方法も、信頼領域で極小値を見つけるためのアルゴリズムによって実装されており、大規模なタスクに適しています。

一般形式での最小探索問題の数学的定式化：

m i n_{x} f （ x ）

$\ min_x f（x）$

c^{l} l e q c （ x ） l e q c^{u} 、

$c ^ l \ leq c（x）\ leq c ^ u、$

x^{l} l e q x l e q x^{u}

$x ^ l \ leq x \ leq x ^ u$

厳密な等式制約の場合、下限は上限に等しく設定されます $c ^ l_j = c ^ u_j$ 。
一方向の制限の場合、上限または下限はnp.infによって対応する記号で設定されます。
2つの変数の既知のRosenbrock関数の最小値を見つける必要があるとします。

m i n_{x_{0} 、 x_{1}} 100 （ x_{0} - x_{1}^{2} ）^{2} + （ 1 - x_{0} ）^{2}

$\ min_ {x_0、x_1} 100（x_0-x_1 ^ 2）^ 2 +（1-x_0）^ 2$

この場合、次の制限が定義のドメインに設定されます。

x_{0}^{2} + x_{1} l e q 1

$x_0 ^ 2 + x_1 \ leq 1$

x_{0}^{2} - x_{1} l e q 1

$x_0 ^ 2-x_1 \ leq 1$

2 x_{0} + x_{1} = 1

$2x_0 + x_1 = 1$

x_{0} + 2 x_{1} l e q 1

$x_0 + 2x_1 \ leq 1$

0 l e q x_{0} l e q 1

$0 \ leq x_0 \ leq 1$

- 0.5 l e q x_{1} l e q 2.0

$-0.5 \ leq x_1 \ leq 2.0$

私たちの場合、その時点でユニークなソリューションがあります $[x_0、x_1] = [0.4149、0.1701]$ 最初と4番目の制限のみが有効です。
制限を下から順に見て、scipyでそれらをどのように書くことができるかを見てみましょう。
制限事項 $0 \ leq x_0 \ leq 1$ そして $-0.5 \ leq x_1 \ leq 2.0$ Boundsオブジェクトを使用して定義されます。

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

制限事項 $x_0 + 2 x_1 \ leq 1$ そして $2 x_0 + x_1 = 1$ 線形形式で記述します。

b e g i n b m a t r i x - i n f t y 1 e n d b m a t r i x l e q b e g i n b m a t r i x 1 ＆ 2 2 ＆ 1 e n d b m a t r i x b e g i n b m a t r i x x_{0} x_{1} e n d b m a t r i x l e q b e g i n b m a t r i x 1 1 e n d b m a t r i x

$\ begin {bmatrix}-\ infty \\ 1 \ end {bmatrix} \ leq \ begin {bmatrix} 1＆2 \\ 2＆1 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \ end {bmatrix } \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$

これらの制約をLinearConstraintオブジェクトとして定義します。

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

そして最後に、行列形式の非線形制約：

c （ x ） = b e g i n b m a t r i x x_{0}^{2} + x_{1} x_{0}^{2} - x_{1} e n d b m a t r i x l e q b e g i n b m a t r i x 1 1 e n d b m a t r i x

$c（x）= \ begin {bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2-x_1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$

この制約のヤコビ行列と、ヘッセ行列と任意のベクトルの線形結合を定義します $v$ ：

J （ x ） = b e g i n b m a t r i x 2 x_{0} ＆ 1 2 x_{0} ＆ - 1 e n d b m a t r i x

$J（x）= \ begin {bmatrix} 2x_0＆1 \\ 2x_0＆-1 \ end {bmatrix}$

H （ x 、 v ） = s u m l i m i t s_{i = 0}^{1} v_{i} n a b l a^{2} c_{i} （ x ） = v_{0} b e g i n b m a t r i x 2 ＆ 0 2 ＆ 0 e n d b m a t r i x + v_{1} b e g i n b m a t r i x 2 ＆ 0 2 ＆ 0 e n d b m a t r i x

$H（x、v）= \ sum \ limits_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i（x）= v_0 \ begin {bmatrix} 2＆0 \\ 2＆0 \ end {bmatrix} + v_1 \ begin {bmatrix} 2＆0 \\ 2＆0 \ end {bmatrix}$

これで、非線形制約を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)

ヘッセ行列の計算 $H（x、v）$ コストがかかるため、 HessianUpdateStrategyクラスを使用できます。次の戦略を使用できます： BFGSおよびSR1 。

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

ヘッセ行列は、有限差分を使用して計算することもできます。

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

制約のヤコビ行列は、有限差分を使用して計算することもできます。ただし、この場合、ヘッセ行列は有限差分を使用して計算できません。ヘッセ行列は、関数として定義するか、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クラスを使用して定義できます

 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)

または、ヘッセ行列とhesspパラメーターを介した任意のベクトルの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)

あるいは、最適化された関数の1次および2次導関数を近似的に計算できます。たとえば、ヘッシアンは関数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メソッドは、次の形式で関数を最小化する問題を解決するように設計されています。

m i n_{x} f （ x ）

$\ min_x f（x）$

c_{j} （ x ） = 0 、 j i n m a t h c a l E

$c_j（x）= 0、j \ in \ mathcal {E}$

c_{j} （ x ） g e q 0 、 j i n m a t h c a l I

$c_j（x）\ geq 0、j \ in \ mathcal {I}$

l b_{i} l e q x_{i} l e q u b_{i} 、 i = 1 、 . . . 、 N

$lb_i \ leq x_i \ leq ub_i、i = 1、...、N$

どこで $\ mathcal {E}$ そして $\ mathcal {I}$ -等式または不等式の形式で制約を記述する式インデックスのセット。 $[lb_i、ub_i]$ -関数定義ドメインの下限と上限のセット。

線形および非線形の制約は、キーの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 ]

最適化の例

5番目の技術モードへの移行に関連して、Webスタジオの例を使用して生産の最適化を検討してみましょう。 3つのタイプの製品を生産するギャラリーのディレクターとして自分自身を紹介しましょう。

x0-着陸の販売、10 trから
x1-企業サイト、20 trから
x2-オンラインストア、30 trから

私たちのフレンドリーな作業チームには、4人のジュン、2人のミドル、1人のシニアが含まれます。 1か月間の労働時間の基金：

6月： 4 * 150 = 600 * 、
ミドル： 2 * 150 = 300 * 、
シニア： 150 *

タイプ（x0、x1、x2）の1つのサイトの開発と展開に費やす最初の後輩が、最高の（10、20、30）時間、中期-（7、15、20）、上級-（5、10、15）時間を費やすと仮定します彼の人生の時間。

通常のディレクターと同様に、私たちは月間利益を最大化したいと考えています。成功への最初のステップは、目的関数valueを月次生産からの収入の合計として書き留めることです。

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

これは間違いではなく、最大値を検索するとき、目的関数は反対の符号で最小化されます。

次のステップは、従業員への処理を禁止し、労働時間資金に制限を導入することです。

b e g i n b m a t r i x 10 ＆ 20 ＆ 30 7 ＆ 15 ＆ 20 5 ＆ 10 ＆ 15 e n d b m a t r i x b e g i n b m a t r i x x_{0} x_{1} x_{2} e n d b m a t r i x l e q b e g i n b m a t r i x 600 300 150 e n d b m a t r i x

$\ begin {bmatrix} 10＆20＆30 \\ 7＆15＆20 \\ 5＆10＆15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ leq \ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}$

どちらが同等ですか：

b e g i n b m a t r i x 600 300 150 e n d b m a t r i x - b e g i n b m a t r i x 10 お よ び 20 お よ び 30 7 お よ び 15 お よ び 20 5 お よ び 10 お よ び 15 e n d b m a t r i x b e g i n b m a t r i x x_{0} x_{1} x_{2} e n d b m a t r i x g e q 0

$\ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}-\ begin {bmatrix} 10および20および30 \\ 7および15および20 \\ 5および10および15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ geq 0$

 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、6、3）の最適な分布でrowぎ手の毎月の負荷を計算します。

6月： 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で、6月の半分未満です。

おわりに

この記事では、条件付き最小化の問題を解決するために使用されるscipy.optimizeパッケージを使用するための基本的な手法の概要を説明します。個人的には、私は学術的な目的でのみscipyを使用しているため、上記の例はとてもscipy isです。

理論とwinrarの例の多くは、たとえばI. L. Akulichの本「例と問題における数学プログラミング」にあります。画像セット（ハブに関する記事）の3D構造を構築するためのscipy.optimizeよりハードコアなアプリケーションは、 scipy-cookbookにあります。

主な情報源はdocs.scipy.orgであり、 scipyこのセクションおよび他のセクションを翻訳に翻訳する場合は、 GitHubにようこそ。

この出版物に貢献してくれたmephistopheiesに感謝します。

SciPy、条件付き最適化

はじめに

条件付き最適化方法= "trust-constr"

条件付き最適化方法= "SLSQP"

最適化の例

おわりに

More articles: