SciPy (ausgesprochen "sai pai") ist ein auf NumPy basierendes Mathematikpaket, das auch Bibliotheken in C und Fortran umfasst. Mit SciPy wird eine interaktive Python-Session zu einer vollwertigen Datenverarbeitungsumgebung wie MATLAB, IDL, Octave, R oder SciLab.
In diesem Artikel werfen wir einen Blick auf die grundlegenden Techniken der mathematischen Programmierung – insbesondere die Lösung von Problemen der bedingten Optimierung für eine skalarwertige Funktion mehrerer Variablen mithilfe des Pakets scipy.optimize. Die Algorithmen der unbedingten Optimierung wurden bereits in . Eine detaillierte und aktuelle Dokumentation der Funktionen von SciPy ist jederzeit über den Befehl help(), Shift+Tab oder in .
Einführung
erhältlich. Eine allgemeine Schnittstelle zur Lösung sowohl der bedingten als auch der unbedingten Optimierungsprobleme im Paket scipy.optimize wird durch die Funktion minimize()bereitgestellt. Es ist jedoch bekannt, dass es keinen universellen Ansatz zur Lösung aller Probleme gibt, weshalb die Wahl der geeigneten Methode wie immer in der Verantwortung des Forschers liegt.
Der geeignete Optimierungsalgorithmus wird über das Argument der Funktion minimize(..., method="").
festgelegt. Für die bedingte Optimierung von Funktionen mehrerer Variablen stehen die folgenden Methoden zur Verfügung:
trust-constr— Suche nach dem lokalen Minimum im vertrauenswürdigen Bereich. , ;SLSQP— sequenzielle quadratische Programmierung mit Einschränkungen, Newtons Methode zur Lösung des Lagrange-Systems. .TNC— Truncated Newton Constrained, eine beschränkte Anzahl von Iterationen, geeignet für nichtlineare Funktionen mit vielen unabhängigen Variablen. .L-BFGS-B— Methode der Broyden–Fletcher–Goldfarb–Shanno-Vierers, implementiert mit reduziertem Speicherbedarf durch teilweise Speicherung der Vektoren aus der Hesse-Matrix. , .COBYLA— COBYLA, die begrenzte Optimierung durch lineare Approximation (ohne Berechnung des Gradienten). .
Je nach gewählter Methode werden die Bedingungen und Einschränkungen zur Lösung des Problems unterschiedlich definiert:
- Objekt der Klasse
Boundsfür die Methoden L-BFGS-B, TNC, SLSQP, trust-constr; - Liste
(min, max)für dieselben Methoden L-BFGS-B, TNC, SLSQP, trust-constr; - Objekt oder Liste von Objekten
LinearConstraint,NonlinearConstraintfür die Methoden COBYLA, SLSQP, trust-constr; - Wortbuch oder Liste von Wörterbüchern
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}für die Methoden COBYLA, SLSQP.
Artikelübersicht:
1) Betrachtung der Anwendung des bedingten Optimierungsalgorithmus im Vertrauensbereich (method=»trust-constr») mit durch Objekte definierten Restriktionen. Bounds, LinearConstraint, NonlinearConstraint ;
2) Analyse der sequenziellen Programmierung mit der Methode der kleinsten Quadrate (method=»SLSQP») unter Verwendung von Dictionary-definierten Einschränkungen. {'type', 'fun', 'jac', 'args'};
3) Analyse eines Beispiels zur Optimierung der Produktionsausgabe anhand eines Webstudios.
Bedingte Optimierung mit der Methode=»trust-constr»
Implementierung der Methode trust-constr basiert auf für Probleme mit Gleichheitsrestriktionen und auf für Probleme mit Ungleichheitsrestriktionen. Beide Methoden nutzen Algorithmen zur lokalen Minimalwertsuche im Vertrauensbereich und eignen sich gut für groß angelegte Probleme.
Mathematische Formulierung des Problems der Minimierung im Allgemeinen:
Für Gleichheitsrestriktionen wird die untere Grenze auf die obere Grenze gesetzt. 
.
Für einseitige Einschränkungen wird die obere oder untere Grenze gesetzt auf np.inf mit dem entsprechenden Vorzeichen.
Angenommen, es muss das Minimum der bekannten Rosenbrock-Funktion in Bezug auf zwei Variablen gefunden werden:
Dabei sind folgende Einschränkungen für ihren Definitionsbereich gegeben:
В нашем случае имеется единственное решение в точке 
, для которой справедливы только первое и четвертое ограничения.
Пройдемся по ограничениям снизу вверх и рассмотрим, как можно их записать в scipy.
Einschränkungen 
und 
определим с помощью объекта Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Einschränkungen 
und 
запишем в линейной форме:
Определим эти ограничения в виде объекта 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)Wenn die Berechnung der Hesse-Matrix 
hohe Kosten verursacht, kann die Klasse verwendet werden. Folgende Strategien sind verfügbar: BFGS und SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Die Hesse-Matrix kann auch mittels endlicher Differenzen berechnet werden:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Die Jacobi-Matrix für die Nebenbedingungen kann ebenfalls mit endlichen Differenzen berechnet werden. In diesem Fall kann die Hesse-Matrix jedoch nicht mit endlichen Differenzen berechnet werden. Der Hesse muss als Funktion oder mit Hilfe der Klasse HessianUpdateStrategy definiert werden.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Die Lösung des Optimierungsproblems sieht folgendermaßen aus:
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` Abbruchbedingung ist erfüllt.
Anzahl der Iterationen: 12, Funktionsauswertungen: 8, CG-Iterationen: 7, Optimalität: 2.99e-09, Verletzung der Einschränkung: 1.11e-16, Ausführungszeit: 0.033 s.
[0.41494531 0.17010937]Falls erforderlich, kann die Funktion zur Berechnung des Hessian mit der Klasse LinearOperator definiert werden.
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)oder das Produkt des Hessian und eines beliebigen Vektors über das Parameterfeld. 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)Alternativ können die ersten und zweiten Ableitungen der zu optimierenden Funktion approximiert werden. Zum Beispiel kann der Hessian durch eine Funktion approximiert werden. SR1 (quasi-Newton-Anäherung). Der Gradient kann mit Hilfe von Finite-Differenzen approximiert werden.
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)Bedingte Optimierung method="SLSQP"
Die SLSQP-Methode dient zur Lösung von Minimierungsproblemen in der Form:
Dabei ist 
und 
— eine Menge von Indexausdrücken, die Einschränkungen in Form von Gleichungen oder Ungleichungen beschreiben. 
— eine Menge von unteren und oberen Grenzen für den Definitionsbereich der Funktion.
Lineare und nichtlineare Einschränkungen werden in Form von Dictionaries mit den Schlüsseln beschrieben: type, . Bei und 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])
}Die Minimierung erfolgt wie folgt:
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)Optimierung erfolgreich abgeschlossen. (Exit-Modus 0)
Aktueller Funktionswert: 0.34271757499419825
Iterationen: 4
Funktionsauswertungen: 5
Gradientenbewertungen: 4
[0.41494475 0.1701105 ]Beispiel für Optimierung
Im Zuge des Übergangs zur fünften technologischen Ära betrachten wir die Produktionsoptimierung anhand eines Beispiels aus einer Web-Agentur, die uns ein kleines, aber stabiles Einkommen sichert. Stellen wir uns als Direktor einer Galerie vor, in der drei Produktarten hergestellt werden:
- x0 — Verkaufs-Landingpages, ab 10.000 €
- x1 — Unternehmenswebseiten, ab 20.000 €
- x2 — Online-Shops, ab 30.000 €
Unser eingespieltes Team besteht aus vier Junior Entwicklern, zwei Mid-Level und einem Senior Entwickler. Der Arbeitszeitfonds für einen Monat lautet:
- Junior Entwickler:
4 * 150 = 600 Person * Stunden, - Mid-Level:
2 * 150 = 300 Person * Stunden, - Senior Entwickler:
150 Person * Stunden.
Nehmen wir an, dass ein beliebiger Junior für die Entwicklung und das Deployment einer Webseite vom Typ (x0, x1, x2) (10, 20, 30) Stunden, ein Mid-Level (7, 15, 20) Stunden und ein Senior (5, 10, 15) Stunden seine wertvollste Zeit benötigt.
Wie ein normaler Direktor wünschen wir uns, den monatlichen Gewinn zu maximieren. Der erste Schritt zum Erfolg ist, die Zielfunktion aufzustellen value als Summe der Einnahmen aus der im Monat produzierten Ware:
def value(x):
return -10*x[0] - 20*x[1] - 30*x[2]Dies ist kein Fehler; bei der Maximierung wird die Zielfunktion mit umgekehrtem Vorzeichen minimiert.
Der nächste Schritt besteht darin, Überstunden für unsere Mitarbeiter zu verbieten und Arbeitszeitkontingente einzuführen:
Was entspricht:
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]])
}Die formale Einschränkung besagt, dass die Produktion nur positiv sein darf:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Und schließlich die optimistischste Annahme: Aufgrund des niedrigen Preises und der hohen Qualität bilden sich ständig Warteschlangen von zufriedenen Kunden. Wir können die monatlichen Produktionsmengen selbst wählen, basierend auf der Lösung des Optimierungsproblems mit 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]Runden wir auf und berechnen die monatliche Auslastung der Arbeiter in der optimalen Produktionsverteilung x = (8, 6, 3) :
- Junior Entwickler:
8 * 10 + 6 * 20 + 3 * 30 = 290 Personenstunden; - Mid-Level:
8 * 7 + 6 * 15 + 3 * 20 = 206 Personenstunden; - Senior Entwickler:
8 * 5 + 6 * 10 + 3 * 15 = 145 Personenstunden.
Fazit: Damit der Direktor sein verdientes Maximum erhält, sollte er monatlich optimalerweise 8 Landingpages, 6 mittelgroße Webseiten und 3 Onlineshops erstellen. Der Senior sollte dabei durchgehend arbeiten, die Auslastung der Mid-Level-Techniker wird etwa 2/3 betragen, während die Junior-Techniker weniger als die Hälfte ausmachen werden.
Fazit
In diesem Artikel sind die grundlegenden Methoden zur Arbeit mit dem Paket beschrieben scipy.optimize, die zur Lösung von Aufgaben der bedingten Minimierung verwendet werden. Persönlich benutze ich scipy ausschließlich zu akademischen Zwecken, daher hat das gezeigte Beispiel einen humorvollen Charakter.
Viele theoretische Konzepte und Beispiele in WinRAR-Formaten sind beispielsweise im Buch von I.L. Akulich "Mathematische Programmierung in Beispielen und Aufgaben" zu finden. Eine anspruchsvollere Anwendung scipy.optimize zur Erstellung einer 3D-Struktur anhand eines Bilderstapels () kann im .
gefunden werden. Die Hauptquelle für Informationen ist , Interessierte, die zur Übersetzung dieses und anderer Abschnitte beitragen möchten, scipy sind herzlich willkommen auf .
Danke für die Teilnahme an der Veröffentlichungsvorbereitung.
Quelle: habr.com
