SciPy (udtales sai pie) er en numpy-baseret matematikpakke, der også inkluderer C- og Fortran-biblioteker. SciPy gør din interaktive Python-session til et komplet datavidenskabsmiljø som MATLAB, IDL, Octave, R eller SciLab.
I denne artikel vil vi se på de grundlæggende teknikker for matematisk programmering - løsning af betingede optimeringsproblemer for en skalarfunktion af flere variable ved hjælp af scipy.optimize-pakken. Ubegrænsede optimeringsalgoritmer er allerede blevet diskuteret i . Mere detaljeret og opdateret hjælp til scipy-funktioner kan altid fås ved at bruge help()-kommandoen, Shift+Tab eller i .
Indledning
En fælles grænseflade til løsning af både betingede og ubegrænsede optimeringsproblemer i scipy.optimize-pakken leveres af funktionen minimize(). Det er dog kendt, at der ikke er nogen universel metode til at løse alle problemer, så valget af en passende metode falder som altid på forskerens skuldre.
Den passende optimeringsalgoritme er specificeret ved hjælp af funktionsargumentet minimize(..., method="").
Til betinget optimering af en funktion af flere variabler er implementeringer af følgende metoder tilgængelige:
trust-constr— søg efter et lokalt minimum i tillidsområdet. , ;SLSQP— sekventiel kvadratisk programmering med begrænsninger, Newtonsk metode til løsning af Lagrange-systemet. .TNC- Trunkeret Newton Begrænset, begrænset antal iterationer, god til ikke-lineære funktioner med et stort antal uafhængige variable. .L-BFGS-B— en metode fra Broyden–Fletcher–Goldfarb–Shanno-teamet, implementeret med reduceret hukommelsesforbrug på grund af delvis belastning af vektorer fra den hessiske matrix. , .COBYLA— MARE begrænset optimering ved lineær approksimation, begrænset optimering med lineær tilnærmelse (uden gradientberegning). .
Afhængigt af den valgte metode er betingelser og begrænsninger for at løse problemet indstillet forskelligt:
- klasseobjekt
Boundsfor metoderne L-BFGS-B, TNC, SLSQP, trust-constr; - listen
(min, max)for de samme metoder L-BFGS-B, TNC, SLSQP, trust-constr; - et objekt eller en liste over objekter
LinearConstraint,NonlinearConstrainttil COBYLA, SLSQP, tillid-konstr metoder; - ordbog eller liste over ordbøger
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}for COBYLA, SLSQP metoder.
Artikeloversigt:
1) Overvej brugen af en betinget optimeringsalgoritme i tillidsregionen (metode=”trust-constr”) med begrænsninger angivet som objekter Bounds, LinearConstraint, NonlinearConstraint ;
2) Overvej sekventiel programmering ved hjælp af mindste kvadraters metode (metode = "SLSQP") med begrænsninger specificeret i form af en ordbog {'type', 'fun', 'jac', 'args'};
3) Analyser et eksempel på optimering af fremstillede produkter vha. eksemplet med et webstudie.
Betinget optimeringsmetode = "trust-constr"
Implementering af metoden trust-constr baseret på for problemer med begrænsninger af form for lighed og videre for problemer med begrænsninger i form af uligheder. Begge metoder er implementeret af algoritmer til at finde et lokalt minimum i konfidensregionen og er velegnede til store problemer.
Matematisk formulering af problemet med at finde et minimum i generel form:
For strenge lighedsbegrænsninger sættes den nedre grænse lig med den øvre grænse .
For en envejsbegrænsning er den øvre eller nedre grænse sat np.inf med det tilhørende tegn.
Lad det være nødvendigt at finde minimum af en kendt Rosenbrock-funktion af to variable:
I dette tilfælde er følgende begrænsninger sat på dets definitionsdomæne:
I vores tilfælde er der en unik løsning på det punkt , for hvilke kun den første og fjerde begrænsning er gyldig.
Lad os gennemgå begrænsningerne fra bund til top og se på, hvordan vi kan skrive dem i scipy.
Begrænsninger и lad os definere det ved hjælp af Bounds-objektet.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Begrænsninger и Lad os skrive det i lineær form:
Lad os definere disse begrænsninger som 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 den ikke-lineære begrænsning i matrixform:
Vi definerer den jakobiske matrix for denne begrænsning og en lineær kombination af den hessiske matrix med en vilkårlig vektor :
Nu kan vi definere en ikke-lineær 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 matricer også angives i sparsom form:
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:
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)Ved beregning af den hessiske matrix kræver en stor indsats, kan du bruge en klasse . Følgende strategier er tilgængelige: BFGS и SR1.
from 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 endelige forskelle:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Den jakobiske matrix for begrænsninger kan også beregnes ved hjælp af endelige forskelle. Men i dette tilfælde kan den hessiske matrix ikke beregnes ved hjælp af endelige forskelle. Hessian skal defineres som en funktion eller ved at bruge klassen HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Løsningen på optimeringsproblemet ser sådan ud:
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]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 produktet 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 den første og anden afledede af den funktion, der optimeres, tilnærmes. For eksempel kan hessian tilnærmes ved hjælp af funktionen SR1 (kvasi-newtonsk tilnærmelse). Gradienten kan tilnærmes ved endelige forskelle.
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)Betinget optimeringsmetode="SLSQP"
SLSQP-metoden er designet til at løse problemer med at minimere en funktion i form:
hvor и — sæt af indekser af udtryk, der beskriver restriktioner i form af ligheder eller uligheder. — sæt af nedre og øvre grænser for funktionens definitionsdomæne.
Lineære og ikke-lineære begrænsninger er beskrevet i form af ordbøger med nøgler 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])
}Søgningen efter minimum 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)Optimization terminated successfully. (Exit mode 0)
Current function value: 0.34271757499419825
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
[0.41494475 0.1701105 ]Eksempel på optimering
Lad os i forbindelse med overgangen til den femte teknologiske struktur se på produktionsoptimering ved at bruge eksemplet med et webstudie, som giver os en lille, men stabil indkomst. Lad os forestille os os selv som direktør for en kabys, der producerer tre typer produkter:
- x0 - sælger landingssider, fra 10 tr.
- x1 - firmahjemmesider, fra 20 tr.
- x2 - netbutikker, fra 30 tr.
Vores venlige arbejdsteam omfatter fire juniorer, to mellem- og en senior. Deres månedlige arbejdstidsfond:
- juni:
4 * 150 = 600 чел * час, - midterste:
2 * 150 = 300 чел * час, - senior:
150 чел * час.
Lad den første ledige junior bruge (0, 1, 2) timer på udvikling og implementering af ét sted af typen (x10, x20, x30), mellem - (7, 15, 20), senior - (5, 10, 15) ) timer af den bedste tid i dit liv.
Som enhver normal direktør ønsker vi at maksimere den månedlige fortjeneste. Det første skridt til succes er at skrive den objektive funktion ned value som mængden af indkomst fra produkter produceret pr. måned:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Dette er ikke en fejl; når man søger efter maksimum, minimeres objektivfunktionen med det modsatte fortegn.
Næste skridt er at forbyde vores medarbejdere at overarbejde og indføre restriktioner for arbejdstiden:
Hvad 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]])
}En formel begrænsning er, at produktoutput kun må være positivt:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Og endelig er den mest rosenrøde antagelse, at der på grund af den lave pris og høje kvalitet konstant står en kø af tilfredse kunder i kø for os. Vi kan selv vælge de månedlige produktionsmængder, baseret på at løse det begrænsede optimeringsproblem 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 runde løst til hele tal og beregne den månedlige belastning af roere med en optimal fordeling af produkter x = (8, 6, 3) :
- juni:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - midterste:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Konklusion: For at direktøren kan modtage sit velfortjente maksimum, er det optimalt at oprette 8 landingssider, 6 mellemstore sider og 3 butikker om måneden. I dette tilfælde skal senior pløje uden at se op fra maskinen, belastningen af de midterste vil være cirka 2/3, juniorerne mindre end halvdelen.
Konklusion
Artiklen beskriver de grundlæggende teknikker til at arbejde med pakken scipy.optimize, bruges til at løse betingede minimeringsproblemer. Personligt bruger jeg scipy rent akademiske formål, hvorfor det angivne eksempel er af så komisk karakter.
En masse teori og virtuelle eksempler kan for eksempel findes i bogen af I.L. Akulich "Matematisk programmering i eksempler og problemer." Mere hardcore applikation scipy.optimize at bygge en 3D-struktur ud fra et sæt billeder () kan ses i .
Hovedkilden til information er dem, der ønsker at bidrage til oversættelsen af dette og andre afsnit scipy Velkommen til .
Tak for deltagelse i udarbejdelsen af publikationen.
Kilde: www.habr.com