SciPy (udtales som 'sai pai') er et matematisk bibliotek baseret på numpy, som også inkluderer C- og Fortran-biblioteker. Med SciPy bliver en interaktiv Python-session en lige så fuldgyldig databehandlingsmiljø som MATLAB, IDL, Octave, R eller SciLab.
I denne artikel vil vi se på de grundlæggende teknikker inden for matematisk programmering — løsningen af betingede optimeringsopgaver for en skalar funktion af flere variabler ved hjælp af pakken scipy.optimize. De ubetingede optimeringsalgoritmer er allerede blevet beskrevet i . Mere detaljeret og opdateret information om scipy-funktioner kan altid fås ved hjælp af kommandoen help(), Shift+Tab eller i .
Introduktion
Det overordnede interface til løsning af både betingede og ubetingede optimeringsopgaver i pakken scipy.optimize leveres af funktionen minimize(). Det er dog kendt, at der ikke findes én universel metode til at løse alle opgaver, så valget af en passende metode hviler altid på forskerens skuldre.
Den passende optimeringsalgoritme angives med argumentet til funktionen minimize(..., method="").
For betinget optimering af en funktion af flere variabler er følgende metoder tilgængelige:
trust-constr— lokalisering af et minimum i et tillidsområde. , ;SLSQP— sekventiel kvadratisk programmering med begrænsninger, Newtons metode til løsning af Lagrangesystemet. .TNC— Truncated Newton Constrained, et begrænset antal iterationer, godt til ikke-lineære funktioner med mange uafhængige variabler. .L-BFGS-B— metoden fra Broyden–Fletcher–Goldfarb–Shanno familien, implementeret med reduceret hukommelsesforbrug ved delvis indlæsning af vektorer fra Hessians matrix. , .COBYLA— COBYLA Constrained Optimization By Linear Approximation, begrænset optimering med lineær tilnærmelse (uden beregning af gradient). .
Afhængigt af den valgte metode angives betingelser og begrænsninger på forskellige måder:
- ved hjælp af objektet af klassen
Boundsfor metoderne L-BFGS-B, TNC, SLSQP, trust-constr; - som en liste
(min, max)for de samme metoder L-BFGS-B, TNC, SLSQP, trust-constr; - ved hjælp af objektet eller listen af objekter
LinearConstraint,NonlinearConstraintfor metoderne COBYLA, SLSQP, trust-constr; - ved hjælp af et dictionary eller en liste af dictionaries
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}for metoderne COBYLA, SLSQP.
Artikelens plan:
1) Overvej anvendelsen af den betingede optimeringsalgoritme i tillidsområdet (method=»trust-constr») med begrænsninger angivet som objekter. Bounds, LinearConstraint, NonlinearConstraint ;
2) Overvej sekventiel programmering ved hjælp af mindste kvadraters metode (method=»SLSQP») med begrænsninger angivet i form af en ordbog {'type', 'fun', 'jac', 'args'};
3) Gennemgå et eksempel på optimering af produktionen med et webbureau som eksempel.
Betinget optimering method=»trust-constr»
Implementeringen af metoden trust-constr er baseret på til opgaver med lighedsbegrænsninger og på til opgaver med ulighedsbegrænsninger. Begge metoder er implementeret med algoritmer til lokal minima søgning i et tillidsområde og er velegnede til storskala opgaver.
Matematisk formulering af problemet med at finde minimum i generel form:
For begrænsninger med streng lighed sættes den nedre grænse lig med den øvre 
.
For ensidige begrænsninger sættes den øvre eller nedre grænse np.inf med den entsprechende tegn.
Lad os sige, at vi skal finde minimum for den kendte Rosenbrok-funktion med to variable:
Her er der angivet følgende begrænsninger for dens definitionsområde:
I vores tilfælde er der en entydig løsning i punktet 
, hvor kun den første og fjerde begrænsning er gældende.
Lad os gennemgå begrænsningerne fra bunden og se, hvordan vi kan skrive dem i scipy.
Begrænsninger 
og 
defineret med et Bounds-objekt.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Begrænsninger 
og 
lad os skrive det i lineær form:
Definer disse begrænsninger i form af et LinearConstraint-objekt:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Og endelig et ikke-lineært begrænsning i matrixform:
Lad os definere Jacobimatricen for denne begrænsning og en lineær kombination af Hessiamatricen med en vilkårlig vektor 
:
Nu kan vi definere den ikke-lineære begrænsning som et objekt 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)Hvis størrelsen er stor, kan matricerne også angives i sparsomt format:
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)eller som et objekt LinearOperator:
fra 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)Når beregningen af Hessian-matrixen 
kræver betydelige ressourcer, kan klassen . Følgende strategier er tilgængelige: BFGS og SR1.
fra scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian kan også beregnes ved hjælp af finite differencer:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Jacobian-matricen for restriktionerne kan også beregnes ved hjælp af finite differencer. Men i dette tilfælde kan Hessian-matrixen ikke beregnes via finite differencer. Den skal defineres som en funktion eller ved hjælp af klassen HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Løsningen af optimeringsproblemet ser således ud:
fra scipy.optimize import minimize
fra 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 betingelse er tilfredsstillet.
Antal iterationer: 12, funktionsevalueringer: 8, CG-iterationer: 7, optimalitet: 2.99e-09, overtrædelse af restriktioner: 1.11e-16, udførelsestid: 0.033 s.
[0.41494531 0.17010937]Om nødvendigt kan funktionen til beregning af Hessian defineres ved hjælp af klassen 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)eller produktionen af Hessian og en vilkårlig vektor gennem parameteren 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)Alternativt kan de første og anden afledte af den optimerede funktion beregnes approximativt. For eksempel kan Hessian blive approksimeret ved hjælp af funktionen SR1 (quasi-Newton-approximation). Gradientet kan approximativt beregnes ved hjælp af finite differencer.
fra 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)Betinget optimering method="SLSQP"
SLSQP-metoden er designet til at løse optimeringsopgaver i form af:
Where 
og 
— mængder af indeksudtryk, der beskriver begrænsninger i form af lighed eller ulighed. 
— mængder af nedre og øvre grænser for definitionsområdet af funktionen.
Lineære og ikke-lineære begrænsninger beskrives i form af ordbøger med nøgler. type, sjov og 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])
}Minimale søgning udføres som følger:
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)Optimering afsluttet med succes. (Exit mode 0)
Nuværende funktionsværdi: 0.34271757499419825
Iterationer: 4
Funktionsvurderinger: 5
Gradientvurderinger: 4
[0.41494475 0.1701105 ]Eksempel på optimering
I takt med overgangen til den femte teknologiske bølge, lad os se på optimering af produktionen ved hjælp af et eksempel fra et webbureau, der giver os en lille, men stabil indkomst. Forestil dig, at vi er direktør for et galleri, hvor der produceres tre typer produkter:
- x0 — salgssider, fra 10 t.kr.
- x1 — virksomhedshjemmesider, fra 20 t.kr.
- x2 — internetbutikker, fra 30 t.kr.
Vores venlige arbejdsgruppe består af fire juniorer, to mellemledere og en senior. Deres arbejdskraftfond for en måned:
- juniorer:
4 * 150 = 600 mandetimer, - mellemledere:
2 * 150 = 300 mandetimer, - senior:
150 mandetimer.
Lad os sige, at den første junior skal bruge (10, 20, 30) timer på udvikling og implementering af en hjemmeside af typen (x0, x1, x2), mellemlederen — (7, 15, 20), senioren — (5, 10, 15) timer af deres bedste tid.
Som enhver normal direktør vil vi gerne maksimere den månedlige profit. Det første skridt mod succes er at notere den målrettede funktion value som summen af indtægterne fra de producerede produkter i løbet af en måned:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Dette er ikke en fejl; ved søgningen efter maksimumet minimeres den målrettede funktion med omvendt fortegn.
Næste skridt er at forbyde overarbejde for vores ansatte og indføre begrænsninger på arbeidskraftfondet:
Hvilket svarer til:
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]])
}Det formelle begrænsning er, at produktionen kun skal være positiv:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Og endelig, den mest optimistiske antagelse – på grund af den lave pris og høje kvalitet danner der konstant en kø af tilfredse kunder. Vi kan selv vælge de månedlige produktionsmængder baseret på løsningen af en betinget optimeringsopgave med 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]Lad os afrunde til hele tal og beregne den månedlige belastning af roere ved den optimale produktionsfordeling x = (8, 6, 3) :
- juniorer:
8 * 10 + 6 * 20 + 3 * 30 = 290 personer * timer; - mellemledere:
8 * 7 + 6 * 15 + 3 * 20 = 206 personer * timer; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 personer * timer.
Konklusion: For at direktøren kan få sit velfortjente maksimum, er det optimalt at lave 8 landing pages, 6 mellemstore websteder og 3 butikker om måneden. Senior skal arbejde koncentreret, mens belastningen for mid-level medarbejdere vil være cirka 2/3, og juniorer mindre end halvdelen.
Konklusion
Artiklen behandler de grundlæggende metoder til arbejd med pakken scipy.optimize, der anvendes til at løse opgaver om betinget minimalisering. Personligt bruger jeg scipy udelukkende til akademiske formål, hvorfor det givne eksempel har en sjov karakter.
Meget teori og sjove eksempler kan findes, for eksempel i bogen af I.L. Akulich "Matematisk programmering i eksempler og opgaver". En mere hardcore anvendelse scipy.optimize til at opbygge en 3D-struktur på baggrund af et sæt billeder () kan ses i .
Hovedkilden til information er , de, der ønsker at bidrage til oversættelsen af denne og andre sektioner scipy er velkomne til .
Tak for deltagelse i forberedelsen af publikationen.
Kilde: habr.com
