SciPy (uttalas sai pie) är ett numpy-baserat matematikpaket som även inkluderar C- och Fortran-bibliotek. SciPy förvandlar din interaktiva Python-session till en komplett datavetenskapsmiljö som MATLAB, IDL, Octave, R eller SciLab.
I den här artikeln kommer vi att titta på de grundläggande teknikerna för matematisk programmering - att lösa problem med villkorad optimering för en skalär funktion av flera variabler med hjälp av paketet scipy.optimize. Obegränsade optimeringsalgoritmer har redan diskuterats i . Mer detaljerad och uppdaterad hjälp om scipy-funktioner kan alltid erhållas med hjälp av kommandot help(), Shift+Tab eller i .
Inledning
Ett gemensamt gränssnitt för att lösa både villkorade och obegränsade optimeringsproblem i paketet scipy.optimize tillhandahålls av funktionen minimize(). Det är dock känt att det inte finns någon universell metod för att lösa alla problem, så valet av en adekvat metod faller som alltid på forskarens axlar.
Den lämpliga optimeringsalgoritmen specificeras med funktionsargumentet minimize(..., method="").
För villkorlig optimering av en funktion av flera variabler finns implementeringar av följande metoder tillgängliga:
trust-constr— Sök efter ett lokalt minimum i konfidensregionen. , ;SLSQP— Sekventiell kvadratisk programmering med begränsningar, Newtonsk metod för att lösa Lagrange-systemet. .TNC- Trunkerad Newton Begränsad, begränsat antal iterationer, bra för olinjära funktioner med ett stort antal oberoende variabler. .L-BFGS-B— en metod från Broyden–Fletcher–Goldfarb–Shanno-teamet, implementerad med minskad minnesförbrukning på grund av partiell laddning av vektorer från den hessiska matrisen. , .COBYLA— MARE begränsad optimering genom linjär approximation, begränsad optimering med linjär approximation (utan gradientberäkning). .
Beroende på den valda metoden ställs villkoren och begränsningarna för att lösa problemet på olika sätt:
- klassobjekt
Boundsför metoderna L-BFGS-B, TNC, SLSQP, trust-constr; - listan
(min, max)för samma metoder L-BFGS-B, TNC, SLSQP, trust-constr; - ett objekt eller en lista med objekt
LinearConstraint,NonlinearConstraintför COBYLA, SLSQP, trust-constr metoder; - ordbok eller lista över ordböcker
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}för COBYLA, SLSQP-metoder.
Artikelöversikt:
1) Överväg användningen av en villkorad optimeringsalgoritm i förtroenderegionen (metod=”trust-constr”) med begränsningar specificerade som objekt Bounds, LinearConstraint, NonlinearConstraint ;
2) Överväg sekventiell programmering med minsta kvadratmetoden (metod = "SLSQP") med begränsningar specificerade i form av en ordbok {'type', 'fun', 'jac', 'args'};
3) Analysera ett exempel på optimering av tillverkade produkter med hjälp av exemplet med en webbstudio.
Villkorlig optimeringsmetod = "trust-constr"
Implementering av metoden trust-constr baserat på för problem med begränsningar av formen av jämlikhet och på för problem med begränsningar i form av ojämlikheter. Båda metoderna implementeras av algoritmer för att hitta ett lokalt minimum i konfidensregionen och är väl lämpade för storskaliga problem.
Matematisk formulering av problemet med att hitta ett minimum i allmän form:
För strikta likhetsbegränsningar sätts den nedre gränsen lika med den övre gränsen .
För en enkelriktad begränsning är den övre eller nedre gränsen inställd np.inf med motsvarande tecken.
Låt det vara nödvändigt att hitta minimum av en känd Rosenbrock-funktion av två variabler:
I det här fallet sätts följande begränsningar på dess definitionsdomän:
I vårt fall finns det en unik lösning , för vilka endast den första och fjärde begränsningen är giltiga.
Låt oss gå igenom begränsningarna från botten till toppen och titta på hur vi kan skriva dem i scipy.
Begränsningar и låt oss definiera det med Bounds-objektet.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Begränsningar и Låt oss skriva det i linjär form:
Låt oss definiera dessa begränsningar som ett LinearConstraint-objekt:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Och slutligen den olinjära begränsningen i matrisform:
Vi definierar den jakobiska matrisen för denna begränsning och en linjär kombination av den hessiska matrisen med en godtycklig vektor :
Nu kan vi definiera en icke-linjär begränsning som ett 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)Om storleken är stor kan matriser också anges i gles 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 ett föremål 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)Vid beräkning av den hessiska matrisen kräver mycket ansträngning, du kan använda en klass . Följande strategier är tillgängliga: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian kan också beräknas med ändliga skillnader:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Den jakobiska matrisen för begränsningar kan också beräknas med ändliga skillnader. Men i detta fall kan den hessiska matrisen inte beräknas med ändliga skillnader. Hessian måste definieras som en funktion eller använda klassen HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Lösningen på optimeringsproblemet ser ut så här:
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]Vid behov kan funktionen för beräkning av hessian definieras med 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 produkten av hessian och en godtycklig vektor genom parametern 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örsta och andra derivatan av den funktion som optimeras approximeras. Till exempel kan hessian approximeras med funktionen SR1 (kvasi-newtonsk approximation). Gradienten kan approximeras med ändliga skillnader.
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)Villkorlig optimeringsmetod="SLSQP"
SLSQP-metoden är utformad för att lösa problem med att minimera en funktion i formen:
var и — Uppsättningar av index av uttryck som beskriver begränsningar i form av jämlikheter eller ojämlikheter. — uppsättningar av nedre och övre gränser för definitionsdomänen för funktionen.
Linjära och olinjära begränsningar beskrivs i form av ordböcker med nycklar 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ökandet efter minimum utförs på följande sätt:
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 ]Optimeringsexempel
I samband med övergången till den femte tekniska strukturen, låt oss titta på produktionsoptimering med exemplet med en webbstudio, vilket ger oss en liten men stabil inkomst. Låt oss föreställa oss oss själva som chef för ett kök som producerar tre typer av produkter:
- x0 - säljande målsidor, från 10 tr.
- x1 - företagswebbplatser, från 20 tr.
- x2 - nätbutiker, från 30 tr.
Vårt vänliga arbetslag består av fyra juniorer, två mitten och en senior. Deras månatliga arbetstidsfond:
- juni:
4 * 150 = 600 чел * час, - mitten:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Låt den första tillgängliga junioren spendera (0, 1, 2) timmar på utveckling och driftsättning av en plats av typen (x10, x20, x30), mellan - (7, 15, 20), senior - (5, 10, 15) ) timmar av den bästa tiden i ditt liv.
Som alla vanliga direktörer vill vi maximera månatliga vinster. Det första steget till framgång är att skriva ner den objektiva funktionen value som inkomstbeloppet från produkter som produceras per månad:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Detta är inte ett fel, när man söker efter maximum minimeras objektivfunktionen med motsatt tecken.
Nästa steg är att förbjuda våra anställda att överarbeta och införa restriktioner för arbetstiden:
Vad är likvärdigt:
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 formell begränsning är att produktutdata endast får vara positiv:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Och slutligen, det mest rosenröda antagandet är att på grund av det låga priset och den höga kvaliteten ställer en kö av nöjda kunder ständigt upp för oss. Vi kan välja månatliga produktionsvolymer själva, baserat på att lösa det begränsade optimeringsproblemet 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]Låt oss runda löst till heltal och beräkna den månatliga belastningen av roddare med en optimal fördelning av produkter x = (8, 6, 3) :
- juni:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - mitten:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Slutsats: för att regissören ska få sitt välförtjänta maximum är det optimalt att skapa 8 målsidor, 6 medelstora sajter och 3 butiker per månad. I det här fallet måste senioren plöja utan att titta upp från maskinen, belastningen på mitten blir cirka 2/3, juniorerna mindre än hälften.
Slutsats
Artikeln beskriver de grundläggande teknikerna för att arbeta med paketet scipy.optimize, används för att lösa problem med villkorad minimering. Personligen använder jag scipy rent akademiskt, varför exemplet som ges är av så komisk karaktär.
Många teorier och virtuella exempel finns till exempel i boken av I.L. Akulich "Matematisk programmering i exempel och problem." Mer hardcore applikation scipy.optimize att bygga en 3D-struktur från en uppsättning bilder () kan ses i .
Den huvudsakliga informationskällan är de som vill bidra till översättningen av detta och andra avsnitt scipy Välkommen till .
Tack för deltagande i utarbetandet av publikationen.
Källa: will.com