SciPy (uitgesproken als 'sai pai') is een wiskundige bibliotheek, gebaseerd op numpy, die ook bibliotheken in C en Fortran omvat. Met SciPy transformeert een interactieve Python-sessie zich in een volwaardige data-analyseomgeving, zoals MATLAB, IDL, Octave, R of SciLab.
In dit artikel bekijken we de basisprincipes van wiskundige programmering — het oplossen van optimalisatieproblemen met voorwaarden voor een scalair functie van meerdere variabelen met behulp van de scipy.optimize bibliotheek. Algoritmen voor onvoorwaardelijke optimalisatie zijn al besproken in . Verdere gedetailleerde en actuele informatie over de functies van scipy kan altijd worden verkregen met behulp van de help() functie, Shift+Tab of in .
Inleiding
De algemene interface voor het oplossen van zowel voorwaardelijke als onvoorwaardelijke optimalisatieproblemen in de scipy.optimize bibliotheek wordt geleverd door de functie minimize(). Het is echter bekend dat er geen universele methode bestaat om alle problemen op te lossen, dus de keuze van de juiste methode ligt altijd bij de onderzoeker.
De geschikte optimalisatiealgoritme wordt ingesteld met behulp van het argument van de functie minimize(..., method="").
Voor voorwaardelijke optimalisatie van functies met meerdere variabelen zijn de volgende methoden beschikbaar:
trust-constr— het zoeken naar een lokaal minimum binnen een vertrouwensgebied. , ;SLSQP— sequentieel kwadratisch programmeren met beperkingen, de Newton-methode voor het oplossen van het Lagrange-systeem. .TNC— Truncated Newton Constrained, met een beperkt aantal iteraties, goed voor niet-lineaire functies met veel onafhankelijke variabelen. .L-BFGS-B— de methode van de viertien Broyden–Fletcher–Goldfarb–Shanno, geïmplementeerd met een verminderd geheugenverbruik door gedeeltelijke laadtijd van vectoren uit de Hess-matrix. , .COBYLA— COBYLA Constrained Optimization By Linear Approximation, beperkte optimalisatie met lineaire benadering (zonder het berekenen van de gradient). .
Afhankelijk van de gekozen methode worden verschillende voorwaarden en beperkingen ingesteld voor het oplossen van het probleem:
- met een object van de klasse
Boundsvoor de methoden L-BFGS-B, TNC, SLSQP, trust-constr; - met een lijst van
(min, max)voor dezelfde methoden L-BFGS-B, TNC, SLSQP, trust-constr; - met een object of een lijst van objecten
LinearConstraint,NonlinearConstraintvoor de methoden COBYLA, SLSQP, trust-constr; - met een dictionary of een lijst van dictionaries
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}voor de methoden COBYLA, SLSQP.
Artikelplan:
1) Toepassing van de voorwaardelijke optimalisatie-algoritme in een vertrouwensgebied (method=»trust-constr») met beperkingen, gedefinieerd als objecten Bounds, LinearConstraint, NonlinearConstraint ;
2) Overweeg sequentiële programmering met de methode van de kleinste kwadraten (method="SLSQP") met beperkingen die zijn gedefinieerd in de vorm van een woordenboek {'type', 'fun', 'jac', 'args'};
3) Bespreek een voorbeeld van optimalisatie van de geproduceerde output aan de hand van een webstudio.
Voorwaardelijke optimalisatie method="trust-constr"
Implementatie van de methode trust-constr is gebaseerd op voor problemen met gelijkheidsbeperkingen en op voor problemen met ongelijkheidsbeperkingen. Beide methoden zijn geïmplementeerd met algoritmes voor het zoeken naar een lokaal minimum binnen een betrouwbaarheidsgebied en zijn goed toepasbaar op grootschalige problemen.
Wiskundige formulering van de minimumzoekopdracht in algemene vorm:
Voor strikte gelijkheidsbeperkingen wordt de ondergrens gelijkgesteld aan de bovenste grens 
.
Voor eenzijdige beperkingen wordt de boven- of ondergrens vastgesteld np.inf met het bijbehorende teken.
Stel dat we het minimum van de bekende Rosenbrock-functie van twee variabelen moeten vinden:
Er zijn de volgende beperkingen op het domein:
In ons geval is er één enkele oplossing in het punt 
, waarvoor alleen de eerste en vierde beperkingen waar zijn.
Laten we de beperkingen van onder naar boven doornemen en bekijken hoe we deze in scipy kunnen vastleggen.
Beperkingen 
en 
wordt gedefinieerd met behulp van een Bounds-object.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Beperkingen 
en 
schrijf het in lineaire vorm:
Laten we deze beperkingen definiëren in de vorm van een LinearConstraint-object:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])En tenslotte een niet-lineaire beperking in matrixvorm:
Laten we de Jacobimatrix voor deze beperking en een lineaire combinatie van de Hessian-matrix met een willekeurige vector definiëren 
:
Nu kunnen we de niet-lineaire beperking definiëren als een object 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)Als de grootte groot is, kunnen matrices ook in spaarzame vorm worden gedefinieerd:
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)of als een object 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)Wanneer de berekening van de Hessiaan 
veel middelen vereist, kan de klasse . De volgende strategieën zijn beschikbaar: BFGS en SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())De Hessiaan kan ook worden berekend met behulp van eindige verschillen:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')De Jacobimatrix voor de beperkingen kan ook worden berekend met behulp van eindige verschillen. In dit geval kan de Hessiaan echter niet meer met eindige verschillen worden berekend. De Hessiaan moet worden gedefinieerd als een functie of met behulp van de klasse HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())De oplossing van het optimalisatieprobleem is als volgt:
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` stopcriterium is voldaan.
Aantal iteraties: 12, functieevaluaties: 8, CG-iteraties: 7, optimaliteit: 2.99e-09, schending van de beperking: 1.11e-16, uitvoeringstijd: 0.033 s.
[0.41494531 0.17010937]Indien nodig kan de functie voor het berekenen van de Hessiaan worden gedefinieerd met behulp van de klasse LinearOperator
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)of het product van de Hessiaan en een willekeurige vector via de 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)Alternatief kunnen de eerste en tweede afgeleiden van de te optimaliseren functie benaderend worden berekend. Bijvoorbeeld, de Hessiaan kan worden benaderd met behulp van de functie SR1 (quasi-Newton benadering). De gradient kan worden benaderd met eindige verschillen.
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)Gecodeerde optimalisatie method=»SLSQP»
De SLSQP-methode is bedoeld voor het oplossen van optimalisatieproblemen in de vorm van:
Waarbij 
en 
— verzamelingen van indexen van uitdrukkingen die beperkingen in de vorm van gelijkheden of ongelijkheden beschrijven. 
— verzamelingen van onder- en bovengrenzen voor het domein van de functie.
Lineaire en niet-lineaire beperkingen worden beschreven in de vorm van woordenboeken met sleutels type, fun en 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])
}De zoekactie naar de minimumwaarde gebeurt als volgt:
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)Optimalisatie succesvol beëindigd. (Exit mode 0)
Huidige functiewaarde: 0.34271757499419825
Iteraties: 4
Functiewaarden: 5
Gradiëntwaarden: 4
[0.41494475 0.1701105 ]Voorbeeld van optimalisatie
Met de overgang naar de vijfde technologische structuur, beschouwen we de optimalisatie van productie aan de hand van een webstudio, die ons een klein maar stabiel inkomen oplevert. Stel, we zijn de directeur van een galerie waar drie soorten producten worden vervaardigd:
- x0 — verkoopbare landingspagina's, vanaf 10k.
- x1 — bedrijfswebsites, vanaf 20k.
- x2 — webwinkels, vanaf 30k.
Onze gezellige werkgroep bestaat uit vier junioren, twee mediors en één senior. De beschikbare werkuren per maand zijn:
- junioren:
4 * 150 = 600 persoon * uur, - mediors:
2 * 150 = 300 persoon * uur, - senior:
150 persoon * uur.
Laten we zeggen dat een junior (10, 20, 30) uren nodig heeft voor de ontwikkeling en inzet van een website van het type (x0, x1, x2), een medior (7, 15, 20) en een senior (5, 10, 15) uren van hun beste tijd.
Als elke normale directeur, willen we de maandelijkse winst maximaliseren. De eerste stap naar succes is het opschrijven van de doelstelling. value Als de som van de inkomsten van de geproduceerde producten per maand:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Dit is geen fout, bij het zoeken naar het maximum wordt de doelstelling geminimaliseerd met een omgekeerd teken.
De volgende stap — we verbieden het overwerken van onze werknemers en stellen beperkingen op het beschikbare aantal werkuren in:
Wat equivalent is aan:
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]])
}De formele beperking is dat de productie alleen positief moet zijn:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])En tenslotte, het meest optimistische aanname is dat, vanwege de lage prijs en hoge kwaliteit, we voortdurend een rij tevreden klanten hebben. We kunnen zelf de maandelijkse productievolumes kiezen op basis van de oplossing van de voorwaardelijke optimalisatie taak met 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]Laten we dit globaal afronden naar hele getallen en de maandelijkse belasting van de roeiers berekenen bij de optimale productverdeling x = (8, 6, 3) :
- junioren:
8 * 10 + 6 * 20 + 3 * 30 = 290 personen * uur; - mediors:
8 * 7 + 6 * 15 + 3 * 20 = 206 personen * uur; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 personen * uur.
Conclusie: om ervoor te zorgen dat de directeur zijn verdiende maximum krijgt, is het optimaal om per maand 8 landingspagina's, 6 gemiddelde websites en 3 winkels te maken. De senior moet daarbij onafgebroken werken, de belasting voor de midden-vaardigen zal ongeveer 2/3 bedragen, en voor de juniors minder dan de helft.
Conclusie
In het artikel worden de belangrijkste technieken gepresenteerd voor het werken met het pakket scipy.optimize, gebruikt voor het oplossen van taken van voorwaardelijke minimalisatie. Persoonlijk gebruik ik scipy uit puur academische doeleinden, daarom heeft het gegeven voorbeeld een enigszins humoristisch karakter.
Veel theorie en zeldzame voorbeelden zijn te vinden, bijvoorbeeld in het boek van I.L. Akuliche "Wiskundige programmering in voorbeelden en taken". Voor een hardere toepassing scipy.optimize voor het construeren van een 3D-structuur op basis van een set afbeeldingen () kan men kijken in .
De belangrijkste informatiebron is , degenen die willen bijdragen aan de vertaling van deze en andere secties scipy zijn welkom op .
Dank voor deelname aan de publicatievoorbereiding.
Bron: habr.com
