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- Truncated Newton Constrained, מספר מוגבל של איטרציות, טוב לפונקציות לא ליניאריות עם מספר רב של משתנים בלתי תלויים. .L-BFGS-B- שיטה מצוות Broyden–Fletcher–Goldfarb–Shanno, המיושמת עם צריכת זיכרון מופחתת עקב טעינה חלקית של וקטורים מהמטריצה ההסיאנית. , .COBYLA— MARE מוגבל אופטימיזציה על ידי קירוב ליניארי, אופטימיזציה מוגבלת עם קירוב ליניארי (ללא חישוב שיפוע). .
בהתאם לשיטה שנבחרה, התנאים וההגבלות לפתרון הבעיה נקבעים באופן שונה:
- אובייקט כיתה
Boundsעבור שיטות L-BFGS-B, TNC, SLSQP, trust-constr; - הרשימה
(min, max)עבור אותן שיטות L-BFGS-B, TNC, SLSQP, trust-constr; - חפץ או רשימה של חפצים
LinearConstraint,NonlinearConstraintעבור COBYLA, SLSQP, אמון-constr שיטות; - מילון או רשימת מילונים
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}עבור שיטות COBYLA, SLSQP.
מתווה מאמר:
1) שקול את השימוש באלגוריתם אופטימיזציה מותנה באזור האמון (שיטה="trust-constr") עם אילוצים שצוינו כאובייקטים Bounds, LinearConstraint, NonlinearConstraint ;
2) שקול תכנות רציף בשיטת הריבועים הקטנים ביותר (שיטה = "SLSQP") עם הגבלות שצוינו בצורה של מילון {'type', 'fun', 'jac', 'args'};
3) נתח דוגמה לאופטימיזציה של מוצרים מיוצרים באמצעות הדוגמה של סטודיו אינטרנט.
שיטת אופטימיזציה מותנית = "trust-constr"
יישום השיטה trust-constr מבוסס על לבעיות עם אילוצים של צורת השוויון והלאה לבעיות עם אילוצים בצורה של אי שוויון. שתי השיטות מיושמות על ידי אלגוריתמים למציאת מינימום מקומי באזור הביטחון ומתאימות היטב לבעיות בקנה מידה גדול.
ניסוח מתמטי של הבעיה של מציאת מינימום בצורה כללית:
עבור אילוצי שוויון קפדניים, הגבול התחתון מוגדר שווה לגבול העליון 
.
עבור אילוץ חד כיווני, הגבול העליון או התחתון מוגדר np.inf עם השלט המתאים.
יהיה צורך למצוא את המינימום של פונקציית רוזנברוק ידועה של שני משתנים:
במקרה זה, ההגבלות הבאות נקבעות על תחום ההגדרה שלו:
במקרה שלנו, יש פתרון ייחודי בנקודה 
, שרק ההגבלה הראשונה והרביעית תקפות.
הבה נעבור על ההגבלות מלמטה למעלה ונראה כיצד נוכל לכתוב אותן ב-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])ולבסוף האילוץ הלא ליניארי בצורת מטריצה:
אנו מגדירים את המטריצה היעקוביאנית עבור אילוץ זה ושילוב ליניארי של המטריצה ההסיאנית עם וקטור שרירותי 
:
כעת נוכל להגדיר אילוץ לא ליניארי כאובייקט 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)בעת חישוב המטריצה ההסיאנית 
דורש הרבה מאמץ, אתה יכול להשתמש בכיתה . האסטרטגיות הבאות זמינות: 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')ניתן לחשב את המטריצה היעקוביאנית עבור אילוצים גם באמצעות הבדלים סופיים. עם זאת, במקרה זה לא ניתן לחשב את המטריצה ההסיאנית באמצעות הבדלים סופיים. יש להגדיר את ה-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]במידת הצורך, ניתן להגדיר את הפונקציה לחישוב ה-Hssian באמצעות המחלקה 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:
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 נועדה לפתור בעיות של מזעור פונקציה בצורה:
שם 
и 
- קבוצות של מדדים של ביטויים המתארים הגבלות בצורה של שוויון או אי שוויון. 
- קבוצות של גבול תחתון ועליון עבור תחום ההגדרה של הפונקציה.
אילוצים ליניאריים ולא ליניאריים מתוארים בצורה של מילונים עם מפתחות 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 לבנות מבנה תלת מימדי מקבוצת תמונות () ניתן לצפות ב .
מקור המידע העיקרי הוא אלה המעוניינים לתרום לתרגום של סעיפים זה ואחרים scipy ברוך הבא ל .
תודה להשתתפות בהכנת הפרסום.
מקור: www.habr.com
