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:
In unserem Fall gibt es eine einzige Lösung an einem Punkt 
, fĂŒr den nur die erste und die vierte EinschrĂ€nkung gĂŒltig sind.
Lassen Sie uns die EinschrÀnkungen von unten nach oben durchgehen und sehen, wie wir sie in scipy formulieren können.
EinschrÀnkungen 
und 
wir definieren sie mit dem Bounds-Objekt.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])EinschrÀnkungen 
und 
wir schreiben in linearer Form:
Lassen Sie uns diese EinschrÀnkungen in Form des Objekts LinearConstraint definieren:
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 die lineare Kombination der Hessian-Matrix mit einem beliebigen Vektor. 
:
Nun können wir die 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 Dimension groà ist, können die Matrizen auch in spÀrlicher 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)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
