SciPy (ausgesprochen wie "sai pai") ist ein auf numpy basierendes mathematisches Paket, das auch Bibliotheken in C und Fortran umfasst. Mit SciPy wird die interaktive Python-Sitzung zu einer vollstÀndigen Datenverarbeitungsumgebung wie MATLAB, IDL, Octave, R oder SciLab.
In diesem Artikel werden die grundlegenden Techniken der mathematischen Programmierung behandelt â die Lösung von Problemen der bedingten Optimierung fĂŒr eine skalare Funktion mehrerer Variablen mithilfe des Pakets scipy.optimize. Die Algorithmen der unbedingten Optimierung wurden bereits in . Eine detaillierte und aktuelle Dokumentation zu den Funktionen von scipy kann jederzeit mit dem Befehl help(), Shift+Tab oder in .
EinfĂŒhrung
Die allgemeine Schnittstelle zur Lösung sowohl bedingter als auch unbedingter 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 des geeigneten Verfahrens wie immer in der Verantwortung des Forschers liegt.
Der geeignete Optimierungsalgorithmus wird durch das Argument der Funktion minimize(..., method="").
fĂŒr die bedingte Optimierung einer Funktion mehrerer Variablen stehen Umsetzungen der folgenden Methoden zur VerfĂŒgung:
trust-constrâ Suche nach einem lokalen Minimum im Vertrauenbereich. , ;SLSQPâ sequenzielle quadratische Programmierung mit EinschrĂ€nkungen, Newtons Methode zur Lösung des Lagrangesystems. .TNCâ Truncated Newton Constrained, begrenzte Anzahl von Iterationen, gut fĂŒr nichtlineare Funktionen mit einer groĂen Anzahl unabhĂ€ngiger Variablen. .L-BFGS-Bâ ein Verfahren von der Familie BroydenâFletcherâGoldfarbâShanno, das mit reduziertem Speicherbedarf durch partielle Beladung der Vektoren aus der Hessischen Matrix umgesetzt wird. , .COBYLAâ COBYLA Constrained Optimization By Linear Approximation, eingeschrĂ€nkte Optimierung mit linearer Approximation (ohne Berechnung des Gradienten). .
Je nach gewĂ€hlter Methode werden die Bedingungen und EinschrĂ€nkungen fĂŒr die Problemlösung unterschiedlich festgelegt:
- ein Objekt der Klasse
BoundsfĂŒr die Methoden L-BFGS-B, TNC, SLSQP, trust-constr; - eine Liste
(min, max)fĂŒr dieselben Methoden L-BFGS-B, TNC, SLSQP, trust-constr; - ein Objekt oder eine Liste von Objekten
LinearConstraint,NonlinearConstraintfĂŒr die Methoden COBYLA, SLSQP, trust-constr; - ein Wörterbuch oder eine Liste von WörterbĂŒchern
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}fĂŒr die Methoden COBYLA, SLSQP.
ArtikelĂŒbersicht:
1) Die Anwendung des Algorithmus zur bedingten Optimierung im Vertrauenbereich (method=»trust-constr») mit EinschrÀnkungen, die als Objekte definiert sind, in Betracht ziehen. Bounds, LinearConstraint, NonlinearConstraint ;
2) Betrachten Sie sequenzielle Programmierung nach der Methode der kleinsten Quadrate (method=»SLSQP») mit den in Form eines Dictionaries angegebenen EinschrÀnkungen {'type', 'fun', 'jac', 'args'};
3) Untersuchen Sie ein Beispiel zur Optimierung der produzierten Waren am Beispiel eines Webstudios.
Bedingte Optimierung method=»trust-constr»
Implementierung der Methode trust-constr basiert auf fĂŒr Aufgaben mit GleichheitsbeschrĂ€nkungen und auf fĂŒr Aufgaben mit UngleichheitsbeschrĂ€nkungen. Beide Methoden sind durch Algorithmen zur Suche nach lokalen Minima im Vertrauensbereich implementiert und eignen sich gut fĂŒr groĂangelegte Aufgaben.
Die mathematische Formulierung des Problems der Minimalwertsuche im Allgemeinen lautet:
FĂŒr strikte GleichheitsbeschrĂ€nkungen wird die untere Grenze auf die obere Grenze gesetzt 
.
FĂŒr einseitige EinschrĂ€nkungen wird die obere oder untere Grenze festgelegt np.inf mit dem entsprechenden Zeichen.
Angenommen, es muss das Minimum der bekannten Rosenbrock-Funktion von zwei Variablen gefunden werden:
Dabei sind die folgenden EinschrĂ€nkungen fĂŒr ihren Definitionsbereich gegeben:
In unserem Fall gibt es eine eindeutige Lösung im Punkt 
, fĂŒr den nur die erste und die vierte EinschrĂ€nkung gelten.
Lassen Sie uns die EinschrÀnkungen von unten nach oben durchgehen und betrachten, wie wir sie in scipy formulieren können.
EinschrÀnkungen 
und 
Definieren mit Hilfe des Objekts Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])EinschrÀnkungen 
und 
Schreiben wir in linearer Form:
Definieren wir diese EinschrÀnkungen in Form eines Objekts LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Und schlieĂlich eine nichtlineare EinschrĂ€nkung in matrizinĂ€rer Form:
Definieren wir die Jacobi-Matrix fĂŒr diese EinschrĂ€nkung und die lineare Kombination der Hessischen Matrix mit einem beliebigen Vektor 
:
Jetzt 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 GröĂe 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())Der Hesse kann auch mit Hilfe von endlichen Differenzen berechnet werden:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Die Jacobimatrix fĂŒr EinschrĂ€nkungen kann ebenfalls mit Hilfe von endlichen Differenzen berechnet werden. In diesem Fall kann die Hesse-Matrix jedoch nicht mehr mit Hilfe von endlichen Differenzen berechnet werden. Der Hesse muss als Funktion oder ĂŒber die Klasse HessianUpdateStrategy definiert werden.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ( ))Die Lösung des Optimierungsproblems sieht wie folgt 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 wurde erfĂŒllt.
Anzahl der Iterationen: 12, Funktionsauswertungen: 8, CG-Iteration: 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 Hesse mit Hilfe 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 der Hesse-Matrix und eines beliebigen Vektors ĂŒber 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 ersten und zweiten Ableitungen der zu optimierenden Funktion approximiert werden. Beispielsweise kann der Hesse mit Hilfe der Funktion SR1 (quasi-Newton-Approximation) approximiert werden. Der Gradient kann mit endlichen 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 Methode SLSQP dient zur Lösung von Minimierungsproblemen der Funktion in Form von:
Wo 
und 
â Mengen von Indexen von AusdrĂŒcken, die EinschrĂ€nkungen in Form von Gleichungen oder Ungleichungen beschreiben. 
â Mengen von unteren und oberen Grenzen fĂŒr den Definitionsbereich der Funktion.
Lineare und nichtlineare EinschrĂ€nkungen werden als WörterbĂŒcher mit SchlĂŒsseln beschrieben type, SpaĂ 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 beendet. (Austrittsmodus 0)
Aktueller Funktionswert: 0.34271757499419825
Iterationen: 4
Funktionsauswertungen: 5
Gradientenbewertungen: 4
[0.41494475 0.1701105 ]Beispiel fĂŒr Optimierung
Angesichts des Wechsels zum fĂŒnften technologischen Paradigma betrachten wir die Optimierung der Produktion am Beispiel eines Webstudios, das uns einen kleinen, aber stabilen Gewinn beschert. Stellen wir uns als Direktor einer Galerie vor, die drei Arten von Produkten herstellt:
- x0 â Verkaufslandings, ab 10 Tsd. RUB
- x1 â Unternehmenswebsites, ab 20 Tsd. RUB
- x2 â Onlineshops, ab 30 Tsd. RUB
Unser freundliches Arbeitsteam besteht aus vier Junioren, zwei Medioren und einem Senior. Der Stundenansatz fĂŒr ihre Arbeitszeit pro Monat betrĂ€gt:
- Junioren:
4 * 150 = 600 Person-Stunden, - Medioren:
2 * 150 = 300 Person-Stunden, - Senior:
150 Person-Stunden.
Angenommen, ein Junior muss fĂŒr die Entwicklung und das Deployment einer Website vom Typ (x0, x1, x2) (10, 20, 30) Stunden aufwenden, ein Medior (7, 15, 20) und ein Senior (5, 10, 15) Stunden das Beste aus seiner Lebenszeit.
Wie jedem normalen Direktor wollen wir den monatlichen Gewinn maximieren. Der erste Schritt zum Erfolg besteht darin, die Zielfunktion aufzuschreiben, value als Summe der Einnahmen aus der im Monat produzierten Ware:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Das ist kein Fehler, bei der Suche nach dem Maximum wird die Zielfunktion mit umgekehrtem Vorzeichen minimiert.
Der nĂ€chste Schritt besteht darin, unseren Mitarbeitern Ăberstunden zu untersagen und eine Begrenzung der Arbeitszeit festzulegen:
Was Àquivalent ist zu:
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 â die Produktion muss nur positiv sein:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Und schlieĂlich die erfreulichste Annahme - aufgrund des niedrigen Preises und der hohen QualitĂ€t steht stĂ€ndig eine Schlange von zufriedenen Kunden vor unserer TĂŒr. Wir können die monatlichen Produktionsmengen selbst wĂ€hlen, basierend auf der Lösung der Aufgabe der bedingten Optimierung 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 grob auf ganze Zahlen und berechnen wir die monatliche Auslastung der Ruderer bei optimaler Produktverteilung x = (8, 6, 3) :
- Junioren:
8 * 10 + 6 * 20 + 3 * 30 = 290 Personenstunden; - Medioren:
8 * 7 + 6 * 15 + 3 * 20 = 206 Personenstunden; - Senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 Personenstunden.
Fazit: Damit der Direktor seinen verdienten Maximalwert erhĂ€lt, sollten monatlich optimal 8 Landing Pages, 6 mittelgroĂe Websites und 3 Shops erstellt werden. Der Senior sollte dabei ununterbrochen arbeiten, die Auslastung der Mid-Level wird etwa 2/3 betragen, die der Junior weniger als die HĂ€lfte.
Fazit
In dem Artikel sind die grundlegenden Methoden zur Arbeit mit dem Paket scipy.optimize, die zur Lösung von Aufgaben der bedingten Minimierung verwendet werden, enthalten. Persönlich benutze ich scipy ausschlieĂlich zu akademischen Zwecken, weshalb das angefĂŒhrte Beispiel einen humorvollen Charakter hat.
Viel Theorie und ausgefallene Beispiele finden Sie beispielsweise im Buch von I.L. Akulich âMathematische Programmierung in Beispielen und Aufgabenâ. Eine anspruchsvollere Anwendung scipy.optimize zum Aufbau einer 3D-Struktur aus Bildern () kann man in .
Die Hauptquelle fĂŒr Informationen ist , Interessierte, die zur Ăbersetzung dieses und anderer Abschnitte beitragen möchten, scipy sind herzlich eingeladen auf .
Danke fĂŒr die Mitwirkung an der Veröffentlichung.
Quelle: habr.com
