SciPy, Optimierung mit Bedingungen

SciPy (ausgesprochen Sai Pie) ist ein numpy-basiertes Mathematikpaket, das auch C- und Fortran-Bibliotheken enthält. Mit SciPy wird eine 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 betrachtet - das Lösen bedingter Optimierungsprobleme für eine Skalarfunktion mehrerer Variablen mithilfe des scipy.optimize-Pakets. Unbedingte Optimierungsalgorithmen werden bereits im letzten Artikel berücksichtigt. Eine detailliertere und aktuellere Hilfe zu Scipy-Funktionen erhalten Sie immer mit dem Befehl help (), Umschalt + Tab oder in der offiziellen Dokumentation .

Einführung

Die allgemeine Schnittstelle zum Lösen von Problemen sowohl der bedingten als auch der bedingungslosen Optimierung im Paket scipy.optimize wird durch die Funktion minim minimize() bereitgestellt. Es ist jedoch bekannt, dass es keine universelle Methode zur Lösung aller Probleme gibt. Daher liegt die Wahl einer geeigneten Methode wie immer beim Forscher.

Ein geeigneter Optimierungsalgorithmus wird unter Verwendung des Arguments für die minimize(..., method="") .

Zur bedingten Optimierung der Funktion mehrerer Variablen stehen folgende Methoden zur Verfügung:

• trust-constr - Suche nach einem lokalen Minimum im Vertrauensbereich. Wiki- Artikel ; Habr-Artikel ;
• SLSQP - sequentielle quadratische Programmierung mit Einschränkungen, die Newtonsche Methode zur Lösung des Lagrange-Systems. Wiki-Artikel .
• TNC - Truncated Newton Constrained, eine begrenzte Anzahl von Iterationen, eignet sich für nichtlineare Funktionen mit einer großen Anzahl unabhängiger Variablen. Wiki-Artikel .
• L-BFGS-B - eine Methode aus den vier Broyden - Fletcher - Goldfarb - Shanno, die aufgrund des teilweisen Ladens von Vektoren aus der hessischen Matrix mit reduziertem Speicherverbrauch implementiert wurde. Wiki- Artikel, Habr-Artikel .
• COBYLA - Stuten Eingeschränkte Optimierung durch lineare Approximation, eingeschränkte Optimierung mit linearer Approximation (ohne Gradientenberechnung). Wiki-Artikel .

Abhängig von der gewählten Methode sind die Bedingungen und Einschränkungen für die Lösung des Problems unterschiedlich festgelegt:

• ein Objekt der Bounds Klasse für die Methoden L-BFGS-B, TNC, SLSQP, trust-constr;
• die Liste (min, max) für die gleichen Methoden L-BFGS-B, TNC, SLSQP, trust-constr;
• Objekt oder LinearConstraint , NonlinearConstraint für Methoden COBYLA, SLSQP, trust-constr;
• Wörterbuch- oder Wörterbuchliste {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} für die Methoden COBYLA, SLSQP.

Der Umriss des Artikels:
1) Betrachten Sie die Anwendung des bedingten Optimierungsalgorithmus im Konfidenzbereich (method = "trust-constr") mit den Einschränkungen, die in Form von Objekten angegeben sind. Bounds , LinearConstraint , NonlinearConstraint ;
2) Betrachten Sie die sequentielle Programmierung der kleinsten Quadrate (method = "SLSQP") mit den Einschränkungen, die in Form eines Wörterbuchs angegeben sind {'type', 'fun', 'jac', 'args'} ;
3) Zerlegen eines Beispiels für die Optimierung von Produkten am Beispiel eines Webstudios.

Bedingte Optimierungsmethode = "trust-constr"

Die Implementierung der trust-constr basiert auf EQSQP für Probleme mit Einschränkungen der Form der Gleichheit und auf TRIP für Probleme mit Einschränkungen in Form von Ungleichungen. Beide Methoden werden durch Algorithmen zum Finden eines lokalen Minimums im Konfidenzbereich implementiert und eignen sich gut für große Aufgaben.

Die mathematische Formulierung des minimalen Suchproblems in der allgemeinen Form:

$$

$$\ min_x f (x)$$

$$

$$c ^ l \ leq c (x) \ leq c ^ u,$$

$$

$$x ^ l \ leq x \ leq x ^ u$$

Bei strengen Gleichheitsbeschränkungen wird die Untergrenze gleich der Obergrenze gesetzt $$$c ^ l_j = c ^ u_j$.
Für eine np.inf wird eine Ober- oder Untergrenze von np.inf mit dem entsprechenden Vorzeichen festgelegt.
Es sei notwendig, das Minimum der bekannten Rosenbrock-Funktion zweier Variablen zu finden:

$$

$$\ min_ {x_0, x_1} 100 (x_0 - x_1 ^ 2) ^ 2 + (1-x_0) ^ 2$$

In diesem Fall gelten die folgenden Einschränkungen für den Definitionsbereich:

$$

$$x_0 ^ 2 + x_1 \ leq 1$$

$$

$$x_0 ^ 2 - x_1 \ leq 1$$

$$

$$2x_0 + x_1 = 1$$

$$

$$x_0 + 2x_1 \ leq 1$$

$$

$$0 \ leq x_0 \ leq 1$$

$$

$$-0,5 \ leq x_1 \ leq 2,0$$

In unserem Fall gibt es derzeit eine einzigartige Lösung $$$[x_0, x_1] = [0,4149, 0,1701]$für die nur die erste und vierte Einschränkung gelten.
Lassen Sie uns die Einschränkungen von unten nach oben durchgehen und überlegen, wie sie in scipy geschrieben werden können.
Einschränkungen $$$0 \ leq x_0 \ leq 1$und $$$-0,5 \ leq x_1 \ leq 2,0$definiert mit dem Objekt Bounds.

 from scipy.optimize import Bounds bounds = Bounds ([0, -0.5], [1.0, 2.0]) 

Einschränkungen $$$x_0 + 2 x_1 \ leq 1$und $$$2 x_0 + x_1 = 1$wir schreiben in linearer Form:

$$

$$\ begin {bmatrix} - \ infty \\ 1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 & 2 \\ 2 & 1 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \ end {bmatrix } \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Definieren Sie 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 eine nichtlineare Einschränkung in Matrixform:

$$

$$c (x) = \ begin {bmatrix} x_0 ^ 2 + x_1 \\ x_0 ^ 2 - x_1 \ end {bmatrix} \ leq \ begin {bmatrix} 1 \\ 1 \ end {bmatrix}$$

Wir definieren die Jacobi-Matrix für diese Einschränkung und die lineare Kombination der hessischen Matrix mit einem beliebigen Vektor $$$v$::

$$

$$J (x) = \ begin {bmatrix} 2x_0 & 1 \\ 2x_0 & -1 \ end {bmatrix}$$

$$

$$H (x, v) = \ sum \ limit_ {i = 0} ^ 1 v_i \ nabla ^ 2 c_i (x) = v_0 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix} + v_1 \ begin {bmatrix} 2 & 0 \\ 2 & 0 \ end {bmatrix}$$

Jetzt können wir eine nichtlineare Einschränkung als NonlinearConstraint Objekt definieren:

 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 einer spärlichen 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 LinearOperator Objekt:

 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 $$$H (x, v)$teuer können Sie die HessianUpdateStrategy Klasse verwenden. Folgende Strategien stehen zur Verfügung: BFGS und SR1 .

 from scipy.optimize import BFGS nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS()) 

Hessisch kann auch mit endlichen 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 unter Verwendung endlicher Differenzen berechnet werden. In diesem Fall kann die hessische Matrix jedoch nicht mit endlichen Differenzen berechnet werden. Hessisch muss als Funktion oder unter Verwendung der HessianUpdateStrategy-Klasse definiert werden.

 nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ()) 

Die Lösung für das Optimierungsproblem lautet wie folgt:

 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 hessische Berechnungsfunktion mit der LinearOperator-Klasse 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 Hessischen und eines beliebigen Vektors durch den hessp Parameter:

 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 die zweite Ableitung der optimierten Funktion ungefähr berechnet werden. Beispielsweise kann der Hessische Wert mit der Funktion SR1 (quasi-Newtonsche Näherung) approximiert werden. 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 wurde entwickelt, um das Problem der Minimierung der Funktion in Form von:

$$

$$\ min_x f (x)$$

$$

$$c_j (x) = 0, j \ in \ mathcal {E}$$

$$

$$c_j (x) \ geq 0, j \ in \ mathcal {I}$$

$$

$$lb_i \ leq x_i \ leq ub_i, i = 1, ..., N$$

Wo $$$\ mathcal {E}$und $$$\ mathcal {I}$- Sätze von Ausdrucksindizes, die Einschränkungen in Form von Gleichheiten oder Ungleichungen beschreiben. $$$[lb_i, ub_i]$- Sätze von Unter- und Obergrenzen für die Funktionsdefinitionsdomäne.

Lineare und nichtlineare Einschränkungen werden in Form von Wörterbüchern mit den Tasten type , fun 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 Mindestsuche wird wie folgt durchgeführt:

 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 zum fünften technologischen Modus die Produktionsoptimierung am Beispiel eines Webstudios, das uns ein kleines, aber stabiles Einkommen bringt. Stellen wir uns als Direktor einer Galerie vor, die drei Arten von Produkten herstellt:

• x0 - Verkauf von Landungen ab 10 tr
• x1 - Unternehmensstandorte, ab 20 tr
• x2 - Online-Shops, ab 30 tr

Unser freundliches Arbeitsteam besteht aus vier Juns, zwei Middle und einem Senior. Fund ihrer Arbeitszeit für einen Monat:

• Juni: 4 * 150 = 600 * ,
• Mitten: 2 * 150 = 300 * ,
• Senior: 150 *

Angenommen, der erste Junior, der für die Entwicklung und Bereitstellung eines Standorts vom Typ (x0, x1, x2) aufgewendet wird, sollte (10, 20, 30) Stunden, mittlere - (7, 15, 20), ältere - (5, 10, 15) Stunden der Besten verbringen Zeit seines Lebens.

Wie jeder normale Direktor möchten wir unseren monatlichen Gewinn maximieren. Der erste Schritt zum Erfolg besteht darin, den Zielfunktionswert als Summe der Einnahmen aus der monatlichen Produktion aufzuschreiben:

 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 dem entgegengesetzten Vorzeichen minimiert.

Der nächste Schritt besteht darin, die Verarbeitung für unsere Mitarbeiter zu verbieten und den Arbeitszeitfonds einzuschränken:

$$

$$\ begin {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ leq \ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix}$$

Welches ist gleichwertig:

$$

$$\ begin {bmatrix} 600 \\ 300 \\ 150 \ end {bmatrix} - \ begin {bmatrix} 10 & 20 & 30 \\ 7 & 15 & 20 \\ 5 & 10 & 15 \ end {bmatrix} \ begin {bmatrix} x_0 \\ x_1 \\ x_2 \ end {bmatrix} \ geq 0$$

 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]]) } 

Formale Einschränkung - die Ausgabe sollte nur positiv sein:

 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 steht ständig eine Reihe zufriedener Kunden für uns an. Wir können die monatlichen Produktionsmengen selbst auswählen, basierend auf der Lösung des Problems 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] 

Wir runden den Lutscher auf die nächste ganze Zahl und berechnen die monatliche Belastung der Ruderer mit einer optimalen Verteilung der Produktion x = (8, 6, 3) :

• Juni: 8 * 10 + 6 * 20 + 3 * 30 = 290 * ;
• Mitten: 8 * 7 + 6 * 15 + 3 * 20 = 206 * ;
• Senior: 8 * 5 + 6 * 10 + 3 * 15 = 145 * .

Fazit: Damit der Regisseur sein wohlverdientes Maximum erreicht, ist es optimal, 8 Landungen, 6 mittlere Standorte und 3 Geschäfte pro Monat durchzuführen. Gleichzeitig muss der Seigneur pflügen, ohne sich von der Maschine zu lösen. Die mittlere Ladung beträgt etwa 2/3, weniger als die Hälfte des Juni.

Fazit

Der Artikel beschreibt die grundlegenden Techniken für die Arbeit mit dem Paket scipy.optimize , mit denen Probleme der bedingten Minimierung gelöst werden. Persönlich scipy ich scipy nur für akademische Zwecke, daher ist dieses Beispiel so komisch.

Viele theoretische und winrare Beispiele finden sich beispielsweise im Buch von I. L. Akulich "Mathematische Programmierung in Beispielen und Problemen". Eine Hardcore-Anwendung von scipy.optimize zum Erstellen einer 3D-Struktur für eine Reihe von Bildern ( ein Artikel auf dem Hub ) finden Sie im scipy-cookbook .

Die Hauptinformationsquelle ist docs.scipy.org. Wenn Sie diesen und andere Abschnitte von scipy in eine Übersetzung übersetzen scipy , sind Sie bei GitHub willkommen.

Vielen Dank an mephistopheies für ihren Beitrag zu dieser Veröffentlichung.