SciPy ، والتحسين مع الظروف

SciPy (وضوحا sai pie) عبارة عن حزمة رياضيات مبهمة تتضمن أيضًا مكتبات C و Fortran. باستخدام SciPy ، تتحول جلسة Python التفاعلية إلى بيئة كاملة لمعالجة البيانات مثل MATLAB أو IDL أو Octave أو R أو SciLab.

في هذه المقالة ، سننظر في التقنيات الأساسية للبرمجة الرياضية - حل مشكلات التحسين الشرطي لدالة عددية من المتغيرات المتعددة باستخدام الحزمة scipy.optimize. تعتبر خوارزميات التحسين غير المشروطة بالفعل في المقالة الأخيرة . يمكن دائمًا الحصول على مساعدة أكثر تفصيلًا وحديثة حول وظائف scipy باستخدام الأمر help () أو Shift + Tab أو في الوثائق الرسمية .

مقدمة

يتم توفير الواجهة العامة لحل مشاكل التحسين الشرطي وغير المشروط في الحزمة. يتم تحسين الحزمة بواسطة الدالة التصغير minimize() . ومع ذلك ، من المعروف أنه لا توجد طريقة عالمية لحل جميع المشكلات ، لذلك ، كما هو الحال دائمًا ، فإن اختيار الطريقة المناسبة يقع على عاتق الباحث.

تم تحديد خوارزمية تحسين مناسبة باستخدام الوسيطة للدالة minimize(..., method="") .

من أجل التحسين الشرطي لوظيفة العديد من المتغيرات ، تتوفر الطرق التالية:

• trust-constr - ابحث عن الحد الأدنى المحلي في منطقة الثقة. مقالة ويكي ؛ مقالة هبر .
• SLSQP - البرمجة التربيعية المتسلسلة مع القيود ، الطريقة النيوتونية لحل نظام لاغرانج. مقالة ويكي .
• TNC - اقتطاع Newton Restrained ، وهو عدد محدود من التكرارات ، مفيد للوظائف غير الخطية بعدد كبير من المتغيرات المستقلة. مقالة ويكي .
• L-BFGS-B - طريقة من أربعة Broyden - Fletcher - Goldfarb - Shanno ، تم تنفيذها مع تقليل استهلاك الذاكرة بسبب التحميل الجزئي للمتجهات من المصفوفة Hessian. مقالة ويكي ، مقالة هبر .
• COBYLA - MARE التحسين المقيد بواسطة التقريب الخطي ، التحسين الأمثل مع التقريب الخطي (بدون حساب التدرج اللوني). مقالة ويكي .

اعتمادًا على الطريقة المختارة ، يتم تعيين الشروط والقيود الخاصة بحل المشكلة بشكل مختلف:

• كائن من فئة Bounds لأساليب L-BFGS-B و TNC و SLSQP و trust-constr ؛
• قائمة (min, max) لنفس الأساليب L-BFGS-B ، TNC ، SLSQP ، trust-constr ؛
• كائن أو قائمة الكائنات LinearConstraint ، NonlinearConstraint للأساليب COBYLA، SLSQP، trust-constr؛
• القاموس أو القواميس {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} لأساليب COBYLA، SLSQP.

الخطوط العريضة للمقال:
1) النظر في تطبيق خوارزمية التحسين الشرطي في منطقة الثقة (الطريقة = "trust-constr") مع القيود المحددة في شكل Bounds ، LinearConstraint ، NonlinearConstraint الكائنات ؛
2) ضع في اعتبارك برمجة المربعات الصغرى المتسلسلة (method = "SLSQP") مع القيود المحددة في شكل قاموس {'type', 'fun', 'jac', 'args'} ؛
3) لتفكيك مثال لتحسين المنتجات على مثال استوديو الويب.

طريقة التحسين الشرطي = "trust-constr"

يعتمد تطبيق trust-constr على EQSQP للمشاكل المتعلقة بقيود شكل المساواة وعلى TRIP للمشاكل ذات القيود في شكل عدم المساواة. يتم تنفيذ كلتا الطريقتين بواسطة خوارزميات لإيجاد حد أدنى محلي في منطقة الثقة ومناسبة تمامًا للمهام واسعة النطاق.

الصياغة الرياضية لمشكلة البحث الدنيا في النموذج العام:

$$

$$\ min_x f (x)$$

$$

$$c ^ l \ leq c (x) \ leq c ^ u،$$

$$

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

لقيود المساواة الصارمة ، يتم تعيين الحد الأدنى مساوٍ للجزء العلوي $$$c ^ l_j = c ^ u_j$ .
بالنسبة للقيود np.inf يتم تعيين الحد العلوي أو السفلي بواسطة np.inf مع العلامة المقابلة.
دع من الضروري إيجاد الحد الأدنى من وظيفة Rosenbrock المعروفة لمتغيرين:

$$

$$\ min_ {x_0، x_1} 100 (x_0 - x_1 ^ 2) ^ 2 + (1-x_0) ^ 2$$

في هذه الحالة ، يتم تعيين القيود التالية على مجال التعريف الخاص بها:

$$

$$x_0 ^ 2 + x_1 \ leq 1$$

$$

$$x_0 ^ 2 - x_1 \ leq 1$$

$$

$$2x_0 + x_1 = 1 دولا$$

$$

$$x_0 + 2x_1 \ leq 1$$

$$

$$0 \ leq x_0 \ leq 1$$

$$

$$-0.5 \ leq x_1 \ leq 2.0 دولا$$

في حالتنا ، هناك حل فريد في هذه المرحلة $$$[x_0، x_1] = [0.4149 ، 0.1701]$ التي القيود الأولى والرابعة فقط صالحة.
دعنا نذهب من خلال القيود من أسفل إلى أعلى ونفكر في كيفية كتابتها بسرعة.
قيود $$$0 \ leq x_0 \ leq 1$ و $$$-0.5 \ leq x_1 \ leq 2.0 دولا$ المعرفة باستخدام كائن الحدود.

 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 دولا$ نكتب في شكل خطي:

$$

$$\ تبدأ {bmatrix} - \ infty \\ 1 \ end {bmatrix} \ leq \ start {bmatrix} 1 & 2 \\ 2 & 1 \ end {bmatrix} \ start {bmatrix} x_0 \\ x_1 \ end {bmatrix } \ leq \ تبدأ {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) = \ تبدأ {bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2 - x_1 \ end {bmatrix} \ leq \ تبدأ {bmatrix} 1 \\ 1 \ end {bmatrix}$$

نحدد مصفوفة Jacobi لهذا التقييد والمزيج الخطي لمصفوفة Hessian مع متجه تعسفي $$$v$ :

$$

$$J (x) = \ تبدأ {bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \ end {bmatrix}$$

$$

$$H (x، v) = \ sum \ limit_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i (x) = v_0 \ تبدأ {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix} + v_1 \ تبدأ {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') 

يمكن أيضًا حساب مصفوفة Jacobi للقيود باستخدام اختلافات محدودة. ومع ذلك ، في هذه الحالة ، لا يمكن حساب مصفوفة هيس باستخدام اختلافات محدودة. يجب تعريف 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] 

إذا لزم الأمر ، يمكن تعريف وظيفة حساب Hessian باستخدام فئة 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) 

أو منتج 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) 

بدلاً من ذلك ، يمكن حساب المشتقات الأولى والثانية من الوظيفة المُحسّنة تقريبًا. على سبيل المثال ، يمكن تقريب هسي باستخدام الدالة 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 لحل مشكلة تقليل الوظيفة إلى الحد الأدنى في شكل:

$$

$$\ min_x f (x)$$

$$

$$c_j (x) = 0 ، j \ in \ mathcal {E}$$

$$

$$c_j (x) \ geq 0 ، j \ in \ mathcal {I}$$

$$

$$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 ] 

مثال التحسين

فيما يتعلق بالانتقال إلى الوضع التكنولوجي الخامس ، دعونا نفكر في تحسين الإنتاج باستخدام مثال استوديو الويب ، الذي يوفر لنا دخلًا صغيرًا ولكنه ثابت. دعونا نقدم أنفسنا كمدير لمعرض ، والذي ينتج ثلاثة أنواع من المنتجات:

• X0 - بيع الهبوط ، من 10 طن
• X1 - مواقع الشركات ، من 20 طن
• X2 - المخازن على الانترنت ، من 30 طن

يضم فريق العمل الودي لدينا أربعة Juns ، واثنين من كبار الأوسط و Senior. تمويل وقت العمل لمدة شهر:

• يونيو: 4 * 150 = 600 * ،
• المتوسطات: 2 * 150 = 300 * ،
• كبار: 150 *

افترض أن أول مبتدئ ينفق على تطوير ونشر موقع واحد من النوع (x0 ، x1 ، x2) يجب أن يقضي (10 ، 20 ، 30) ساعة ، وسط - (7 ، 15 ، 20) ، أقدم - (5 ، 10 ، 15) ساعة وقت حياته.

مثل أي مدير عادي ، نريد زيادة أرباحنا الشهرية. الخطوة الأولى للنجاح هي تدوين value الوظيفة الموضوعية كمجموع الدخل من الإنتاج الشهري:

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

هذا ليس خطأ ؛ عند البحث عن الحد الأقصى ، يتم تقليل الوظيفة الهدف مع الإشارة المعاكسة.

والخطوة التالية هي حظر المعالجة لموظفينا وفرض قيود على صندوق وقت العمل:

$$

$$\ تبدأ {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ start {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ leq \ تبدأ {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}$$

وهو ما يعادل:

$$

$$\ تبدأ {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix} - \ start {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ start {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 * 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 للأغراض الأكاديمية فقط ، لذلك هذا المثال هزلي للغاية.

يمكن العثور على الكثير من الأمثلة النظرية و winrar ، على سبيل المثال ، في كتاب I. L. Akulich "البرمجة الرياضية في الأمثلة والمشاكل". يمكن العثور على تطبيق أكثر scipy.optimize من scipy.optimize لإنشاء بنية ثلاثية الأبعاد لمجموعة من الصور ( مقالة على المحور ) في كتاب الطبخ ذي النمط .

المصدر الرئيسي للمعلومات هو docs.scipy.org ، إذا كنت ترغب في ترجمة هذا القسم وغيره من أقسام scipy إلى ترجمة ، فمرحبا بكم في GitHub .

شكرا ل mephistopheies للمساهمة في هذا المنشور.