SciPy (uitgesproken als sai pie) is een op numpy gebaseerd wiskundepakket dat ook C- en Fortran-bibliotheken bevat. SciPy verandert uw interactieve Python-sessie in een complete data science-omgeving zoals MATLAB, IDL, Octave, R of SciLab.
In dit artikel zullen we kijken naar de basistechnieken van wiskundig programmeren - het oplossen van voorwaardelijke optimalisatieproblemen voor een scalaire functie van verschillende variabelen met behulp van het scipy.optimize-pakket. Onbeperkte optimalisatie-algoritmen zijn al besproken in . Meer gedetailleerde en actuele hulp over scipy-functies kan altijd worden verkregen met behulp van de opdracht help(), Shift+Tab of in .
Introductie
Een gemeenschappelijke interface voor het oplossen van zowel voorwaardelijke als onbeperkte optimalisatieproblemen in het scipy.optimize-pakket 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 voor een adequate methode valt, zoals altijd, op de schouders van de onderzoeker.
Het juiste optimalisatiealgoritme wordt gespecificeerd met behulp van het functieargument minimize(..., method="").
Voor voorwaardelijke optimalisatie van een functie van meerdere variabelen zijn implementaties van de volgende methoden beschikbaar:
trust-constr— zoeken naar een lokaal minimum in de vertrouwensregio. , ;SLSQP— sequentiële kwadratische programmering met beperkingen, Newtoniaanse methode voor het oplossen van het Lagrange-systeem. .TNC- Afgeknotte Newton Beperkt, beperkt aantal iteraties, goed voor niet-lineaire functies met een groot aantal onafhankelijke variabelen. .L-BFGS-B- een methode van het Broyden-Fletcher-Goldfarb-Shanno-team, geïmplementeerd met verminderd geheugengebruik als gevolg van het gedeeltelijk laden van vectoren uit de Hessische matrix. , .COBYLA— MARE beperkte optimalisatie door lineaire benadering, beperkte optimalisatie met lineaire benadering (zonder gradiëntberekening). .
Afhankelijk van de gekozen methode worden de voorwaarden en beperkingen voor het oplossen van het probleem anders gesteld:
- klasse object
Boundsvoor methoden L-BFGS-B, TNC, SLSQP, trust-constr; - de lijst
(min, max)voor dezelfde methoden L-BFGS-B, TNC, SLSQP, trust-constr; - een object of een lijst met objecten
LinearConstraint,NonlinearConstraintvoor COBYLA-, SLSQP-, trust-constr-methoden; - woordenboek of lijst met woordenboeken
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}voor COBYLA, SLSQP-methoden.
Artikeloverzicht:
1) Overweeg het gebruik van een voorwaardelijk optimalisatie-algoritme in de vertrouwensregio (method=”trust-constr”) met beperkingen gespecificeerd als objecten Bounds, LinearConstraint, NonlinearConstraint ;
2) Overweeg sequentiële programmering met behulp van de kleinste kwadratenmethode (methode = "SLSQP") met beperkingen gespecificeerd in de vorm van een woordenboek {'type', 'fun', 'jac', 'args'};
3) Analyseer een voorbeeld van optimalisatie van gefabriceerde producten met behulp van het voorbeeld van een webstudio.
Voorwaardelijke optimalisatiemethode = "trust-constr"
Implementatie van de methode trust-constr gebaseerd op voor problemen met beperkingen van de vorm van gelijkheid enzovoort voor problemen met beperkingen in de vorm van ongelijkheid. Beide methoden worden geïmplementeerd door algoritmen voor het vinden van een lokaal minimum in het vertrouwensgebied en zijn zeer geschikt voor grootschalige problemen.
Wiskundige formulering van het probleem van het vinden van een minimum in algemene vorm:
Voor strikte gelijkheidsbeperkingen wordt de ondergrens gelijk gesteld aan de bovengrens .
Voor eenrichtingsbeperking wordt de boven- of ondergrens ingesteld np.inf met het bijbehorende teken.
Laat het nodig zijn om het minimum van een bekende Rosenbrock-functie van twee variabelen te vinden:
In dit geval worden de volgende beperkingen gesteld aan het definitiedomein:
In ons geval is er op dit punt een unieke oplossing , waarvoor alleen de eerste en vierde beperking geldig zijn.
Laten we de beperkingen van onder naar boven doornemen en kijken hoe we ze in scipy kunnen schrijven.
Beperkingen и laten we het definiëren met behulp van het Bounds-object.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Beperkingen и Laten we het in lineaire vorm schrijven:
Laten we deze beperkingen definiëren als 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 de niet-lineaire beperking in matrixvorm:
We definiëren de Jacobiaanse matrix voor deze beperking en een lineaire combinatie van de Hessische matrix met een willekeurige vector :
Nu kunnen we een niet-lineaire beperking als een object definiëren 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 maat groot is, kunnen matrices ook in beperkte vorm worden gespecificeerd:
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 voorwerp 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)Bij het berekenen van de Hessische matrix vereist veel inspanning, je kunt een klas gebruiken . De volgende strategieën zijn beschikbaar: BFGS и 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 Jacobiaanse matrix voor beperkingen kan ook worden berekend met behulp van eindige verschillen. In dit geval kan de Hessische matrix echter niet worden berekend met behulp van eindige verschillen. De Hessian 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 voor het optimalisatieprobleem ziet er als volgt uit:
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]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 Hessische 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)Als alternatief kunnen de eerste en tweede afgeleide van de functie die wordt geoptimaliseerd, worden benaderd. De Hessiaan kan bijvoorbeeld worden benaderd met behulp van de functie SR1 (quasi-Newtoniaanse benadering). De gradiënt kan worden benaderd door 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)Voorwaardelijke optimalisatiemethode = "SLSQP"
De SLSQP-methode is ontworpen om problemen op te lossen bij het minimaliseren van een functie in de vorm:
Где и — reeksen indices van uitdrukkingen die beperkingen beschrijven in de vorm van gelijkheden of ongelijkheden. — sets van onder- en bovengrenzen voor het domein van de definitie van de functie.
Lineaire en niet-lineaire beperkingen worden beschreven in de vorm van woordenboeken met sleutels 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])
}Het zoeken naar het minimum wordt als volgt uitgevoerd:
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 ]Optimalisatie voorbeeld
Laten we, in verband met de overgang naar de vijfde technologische structuur, eens kijken naar productie-optimalisatie aan de hand van het voorbeeld van een webstudio, die ons een klein maar stabiel inkomen oplevert. Laten we ons voorstellen dat we de directeur zijn van een kombuis die drie soorten producten produceert:
- x0 - landingspagina's verkopen, vanaf 10 tr.
- x1 - bedrijfswebsites, vanaf 20 tr.
- x2 - online winkels, vanaf 30 tr.
Ons vriendelijke team bestaat uit vier junioren, twee middenklassers en één senior. Hun maandelijkse werktijdfonds:
- juni:
4 * 150 = 600 чел * час, - middens:
2 * 150 = 300 чел * час, - señor:
150 чел * час.
Laat de eerste beschikbare junior (0, 1, 2) uur besteden aan de ontwikkeling en implementatie van één site van het type (x10, x20, x30), midden - (7, 15, 20), senior - (5, 10, 15 ) uur van de beste tijd van je leven.
Zoals elke normale directeur willen wij de maandelijkse winst maximaliseren. De eerste stap naar succes is het opschrijven van de objectieve functie value als het bedrag aan inkomsten uit 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 doelfunctie geminimaliseerd met het tegenovergestelde teken.
De volgende stap is het verbieden van overwerk voor onze medewerkers en het invoeren van beperkingen op de werktijden:
Wat is gelijkwaardig:
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]])
}Een formele beperking is dat de productoutput alleen positief mag zijn:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])En tot slot is de meest rooskleurige veronderstelling dat er vanwege de lage prijs en hoge kwaliteit voortdurend een rij tevreden klanten voor ons in de rij staat. We kunnen zelf de maandelijkse productievolumes kiezen, op basis van het oplossen van het beperkte optimalisatieprobleem 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 losjes afronden op hele getallen en de maandelijkse belasting van roeiers berekenen met een optimale verdeling van producten x = (8, 6, 3) :
- juni:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - middens:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - señor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Conclusie: om de directeur zijn welverdiende maximum te laten ontvangen, is het optimaal om 8 landingspagina's, 6 middelgrote sites en 3 winkels per maand te creëren. In dit geval moet de senior ploegen zonder op te kijken van de machine, de belasting van de middenploegen zal ongeveer 2/3 zijn, de junioren minder dan de helft.
Conclusie
In het artikel worden de basistechnieken voor het werken met het pakket beschreven scipy.optimize, gebruikt om voorwaardelijke minimalisatieproblemen op te lossen. Persoonlijk gebruik ik scipy puur voor academische doeleinden, en daarom is het gegeven voorbeeld zo komisch van aard.
Veel theorie en virtuele voorbeelden zijn bijvoorbeeld te vinden in het boek van I.L. Akulich “Wiskundig programmeren in voorbeelden en problemen.” Meer hardcore applicatie scipy.optimize een 3D-structuur opbouwen op basis van een reeks afbeeldingen () kan worden bekeken .
De belangrijkste informatiebron is degenen die willen bijdragen aan de vertaling van deze en andere secties scipy Welkom bij .
Dank voor deelname aan de voorbereiding van de publicatie.
Bron: www.habr.com