SciPy (ausgesprochen „Sai Pie“) ist ein Numpy-basiertes Mathematikpaket, das auch C- und Fortran-Bibliotheken enthält. SciPy verwandelt Ihre interaktive Python-Sitzung in eine vollständige Data-Science-Umgebung wie MATLAB, IDL, Octave, R oder SciLab.
In diesem Artikel betrachten wir die grundlegenden Techniken der mathematischen Programmierung – das Lösen bedingter Optimierungsprobleme für eine Skalarfunktion mehrerer Variablen mithilfe des Pakets scipy.optimize. Unbeschränkte Optimierungsalgorithmen wurden bereits in besprochen . Ausführlichere und aktuellere Hilfe zu Scipy-Funktionen erhalten Sie jederzeit mit dem Befehl help(), Umschalt+Tab oder in .
Einführung
Die Funktion stellt eine gemeinsame Schnittstelle zum Lösen sowohl bedingter als auch uneingeschränkter Optimierungsprobleme im Paket scipy.optimize bereit minimize(). Es ist jedoch bekannt, dass es keine universelle Methode zur Lösung aller Probleme gibt, sodass die Wahl einer geeigneten Methode wie immer auf den Schultern des Forschers liegt.
Der entsprechende Optimierungsalgorithmus wird mit dem Funktionsargument angegeben minimize(..., method="").
Zur bedingten Optimierung einer Funktion mehrerer Variablen stehen Implementierungen der folgenden Methoden zur Verfügung:
trust-constr– Suche nach einem lokalen Minimum im Konfidenzbereich. , ;SLSQP— Sequentielle quadratische Programmierung mit Einschränkungen, Newtonsche Methode zur Lösung des Lagrange-Systems. .TNC- Truncated Newton Constrained, begrenzte Anzahl von Iterationen, gut für nichtlineare Funktionen mit einer großen Anzahl unabhängiger Variablen. .L-BFGS-B– eine Methode des Broyden-Fletcher-Goldfarb-Shanno-Teams, implementiert mit reduziertem Speicherverbrauch aufgrund des teilweisen Ladens von Vektoren aus der Hessischen Matrix. , .COBYLA— MARE Constrained Optimization By Linear Approximation, eingeschränkte Optimierung mit linearer Approximation (ohne Gradientenberechnung). .
Abhängig von der gewählten Methode werden Bedingungen und Einschränkungen zur Lösung des Problems unterschiedlich festgelegt:
- Klasse Objekt
Boundsfür Methoden L-BFGS-B, TNC, SLSQP, trust-constr; - Die Liste
(min, max)für die gleichen Methoden L-BFGS-B, TNC, SLSQP, trust-constr; - ein Objekt oder eine Liste von Objekten
LinearConstraint,NonlinearConstraintfür COBYLA-, SLSQP- und Trust-Constr-Methoden; - Wörterbuch oder Liste von Wörterbüchern
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}für COBYLA, SLSQP-Methoden.
Artikelübersicht:
1) Erwägen Sie die Verwendung eines bedingten Optimierungsalgorithmus in der Vertrauensregion (method=“trust-constr“) mit als Objekten angegebenen Einschränkungen Bounds, LinearConstraint, NonlinearConstraint ;
2) Erwägen Sie die sequentielle Programmierung unter Verwendung der Methode der kleinsten Quadrate (Methode = „SLSQP“) mit in Form eines Wörterbuchs angegebenen Einschränkungen {'type', 'fun', 'jac', 'args'};
3) Analysieren Sie ein Beispiel für die Optimierung hergestellter Produkte am Beispiel eines Webstudios.
Bedingte Optimierungsmethode = „trust-constr“
Implementierung der Methode trust-constr bezogen auf für Probleme mit Einschränkungen der Form der Gleichheit und so weiter für Probleme mit Einschränkungen in Form von Ungleichungen. Beide Methoden werden durch Algorithmen zum Finden eines lokalen Minimums im Vertrauensbereich implementiert und eignen sich gut für großräumige Probleme.
Mathematische Formulierung des Problems der Suche nach einem Minimum in allgemeiner Form:
Für strikte Gleichheitsbeschränkungen wird die Untergrenze gleich der Obergrenze gesetzt .
Für eine unidirektionale Einschränkung wird die Ober- oder Untergrenze festgelegt np.inf mit entsprechendem Schild.
Es sei notwendig, das Minimum einer bekannten Rosenbrock-Funktion zweier Variablen zu finden:
In diesem Fall gelten für seinen Definitionsbereich folgende Einschränkungen:
In unserem Fall gibt es an dieser Stelle eine einzigartige Lösung , für die nur die erste und vierte Einschränkung gelten.
Gehen wir die Einschränkungen von unten nach oben durch und schauen uns an, wie wir sie in Scipy schreiben können.
Einschränkungen и Definieren wir es mit dem Bounds-Objekt.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Einschränkungen и Schreiben wir es in linearer Form:
Definieren wir diese Einschränkungen als LinearConstraint-Objekt:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Und schließlich die nichtlineare Einschränkung in Matrixform:
Wir definieren die Jacobi-Matrix für diese Einschränkung und eine lineare Kombination der Hesse-Matrix mit einem beliebigen Vektor :
Jetzt können wir eine nichtlineare Einschränkung als Objekt definieren 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)Wenn die Größe groß ist, können Matrizen auch in dünner Form angegeben werden:
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)oder als Objekt 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)Bei der Berechnung der Hessischen Matrix erfordert viel Aufwand, Sie können eine Klasse verwenden . Folgende Strategien stehen zur Verfügung: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Der Hessesche Wert kann auch mit Hilfe endlicher Differenzen berechnet werden:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Die Jacobi-Matrix für Einschränkungen kann auch mithilfe endlicher Differenzen berechnet werden. In diesem Fall kann die Hesse-Matrix jedoch nicht mithilfe endlicher Differenzen berechnet werden. Der Hessian muss als Funktion oder mithilfe der HessianUpdateStrategy-Klasse 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` 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]Bei Bedarf kann die Funktion zur Berechnung des Hesse-Werts mithilfe 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 Hesseschen und eines beliebigen Vektors durch den Parameter 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 erste und zweite Ableitung der zu optimierenden Funktion angenähert werden. Beispielsweise kann der Hesse-Wert mit der Funktion angenähert werden SR1 (Quasi-Newtonsche Näherung). Der Gradient kann durch endliche Differenzen angenähert 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 Optimierungsmethode = „SLSQP“
Die SLSQP-Methode soll Probleme der Minimierung einer Funktion in der Form lösen:
Wo и – Sätze von Indizes von Ausdrücken, die Einschränkungen in Form von Gleichheiten oder Ungleichungen beschreiben. – Sätze von Unter- und Obergrenzen für den Definitionsbereich der Funktion.
Lineare und nichtlineare Einschränkungen werden in Form von Wörterbüchern mit Schlüsseln beschrieben 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])
}Die Suche nach dem Minimum 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)Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]Optimierungsbeispiel
Betrachten wir im Zusammenhang mit dem Übergang zur fünften technologischen Struktur die Produktionsoptimierung am Beispiel eines Webstudios, das uns ein kleines, aber stabiles Einkommen beschert. Stellen wir uns vor, wir wären der Direktor einer Galeere, die drei Arten von Produkten herstellt:
- x0 – Verkauf von Landingpages, ab 10 Tr.
- x1 - Unternehmenswebsites, ab 20 Tr.
- x2 - Online-Shops, ab 30 tr.
Zu unserem freundlichen Arbeitsteam gehören vier Junioren, zwei Mittelspieler und ein Senior. Ihr monatlicher Arbeitszeitfonds:
- Juni:
4 * 150 = 600 чел * час, - Mitten:
2 * 150 = 300 чел * час, - Herr:
150 чел * час.
Lassen Sie den ersten verfügbaren Junior (0, 1, 2) Stunden für die Entwicklung und Bereitstellung einer Site vom Typ (x10, x20, x30), Middle – (7, 15, 20), Senior – (5, 10, 15) aufwenden ) Stunden der schönsten Zeit Ihres Lebens.
Wie jeder normale Direktor wollen wir den monatlichen Gewinn maximieren. Der erste Schritt zum Erfolg besteht darin, die Zielfunktion aufzuschreiben value als Höhe des Einkommens aus produzierten Produkten pro Monat:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Dies ist kein Fehler; bei der Suche nach dem Maximum wird die Zielfunktion mit umgekehrtem Vorzeichen minimiert.
Der nächste Schritt besteht darin, unseren Mitarbeitern Überlastung zu verbieten und Arbeitszeitbeschränkungen einzuführen:
Was ist gleichwertig:
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]])
}Eine formale Einschränkung besteht darin, dass der Produktoutput nur positiv sein darf:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Und schließlich ist die rosigste Annahme, dass sich aufgrund des niedrigen Preises und der hohen Qualität ständig eine Schlange zufriedener Kunden für uns anstellt. Wir können die monatlichen Produktionsmengen selbst wählen, basierend auf der Lösung des eingeschränkten 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]Lassen Sie uns locker auf ganze Zahlen runden und die monatliche Belastung der Ruderer bei optimaler Produktverteilung berechnen x = (8, 6, 3) :
- Juni:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - Mitten:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - Herr:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Fazit: Damit der Regisseur sein wohlverdientes Maximum erhält, ist es optimal, 8 Landingpages, 6 mittelgroße Seiten und 3 Shops pro Monat zu erstellen. In diesem Fall muss der Senior pflügen, ohne von der Maschine aufzuschauen, die mittlere Belastung beträgt etwa 2/3, die Junioren weniger als die Hälfte.
Abschluss
Der Artikel beschreibt die grundlegenden Techniken für die Arbeit mit dem Paket scipy.optimize, wird zur Lösung bedingter Minimierungsprobleme verwendet. Persönlich benutze ich scipy rein wissenschaftlicher Natur, weshalb das angeführte Beispiel so komischen Charakter hat.
Viel Theorie und virtuelle Beispiele finden sich beispielsweise im Buch von I.L. Akulich „Mathematische Programmierung in Beispielen und Problemen“. Mehr Hardcore-Anwendung scipy.optimize um eine 3D-Struktur aus einer Reihe von Bildern zu erstellen () kann in eingesehen werden .
Die Hauptinformationsquelle ist diejenigen, die zur Übersetzung dieses und anderer Abschnitte beitragen möchten scipy Willkommen zu .
Danke für die Mitarbeit bei der Erstellung der Publikation.
Source: habr.com