SciPy (prononcita sai torto) estas numpy-bazita matematikpakaĵo kiu ankaŭ inkludas C kaj Fortran-bibliotekojn. SciPy igas vian interagan Python-sesion en kompletan datuman sciencan medion kiel MATLAB, IDL, Octave, R aŭ SciLab.
En ĉi tiu artikolo, ni rigardos la bazajn teknikojn de matematika programado - solvante kondiĉajn optimumigajn problemojn por skalara funkcio de pluraj variabloj uzante la scipy.optimize-pakaĵon. Nelimigitaj optimizaj algoritmoj jam estis diskutitaj en . Pli detala kaj ĝisdata helpo pri scipy-funkcioj ĉiam povas esti akirita per la help() komando, Shift+Tab aŭ en .
Enkonduko
Ofta interfaco por solvi kaj kondiĉajn kaj nelimigitajn optimumigajn problemojn en la scipy.optimize-pakaĵo estas disponigita per la funkcio minimize(). Tamen oni scias, ke ne ekzistas universala metodo por solvi ĉiujn problemojn, do la elekto de taŭga metodo, kiel ĉiam, falas sur la ŝultrojn de la esploristo.
La taŭga optimumiga algoritmo estas specifita uzante la funkcioargumenton minimize(..., method="").
Por kondiĉa optimumigo de funkcio de pluraj variabloj, efektivigoj de la sekvaj metodoj estas haveblaj:
trust-constr— serĉu lokan minimumon en la konfida regiono. , ;SLSQP— sinsekva kvadrata programado kun limoj, Newtoniana metodo por solvi la Lagrange-sistemon. .TNC- Detranĉita Newton Limigita, limigita nombro da ripetoj, bona por neliniaj funkcioj kun granda nombro da sendependaj variabloj. .L-BFGS-B— metodo de la Broyden–Fletcher–Goldfarb–Shanno-teamo, efektivigita kun reduktita memorkonsumo pro parta ŝarĝo de vektoroj de la hesia matrico. , .COBYLA— MARE Constrained Optimization By Linear Approximation, limigita optimumigo kun lineara aproksimado (sen gradienta kalkulo). .
Depende de la elektita metodo, kondiĉoj kaj limigoj por solvi la problemon estas agordita malsame:
- klasobjekto
Boundspor metodoj L-BFGS-B, TNC, SLSQP, trust-constr; - la listo
(min, max)por la samaj metodoj L-BFGS-B, TNC, SLSQP, trust-constr; - objekto aŭ listo de objektoj
LinearConstraint,NonlinearConstraintpor COBYLA, SLSQP, trust-constr-metodoj; - vortaro aŭ listo de vortaroj
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}por COBYLA, SLSQP-metodoj.
Artikola skizo:
1) Konsideru la uzon de kondiĉa optimumiga algoritmo en la fidregiono (metodo = "trust-constr") kun limoj specifitaj kiel objektoj Bounds, LinearConstraint, NonlinearConstraint ;
2) Konsideru sinsekvan programadon uzante la metodon de malplej kvadrataj (metodo = "SLSQP") kun limigoj specifitaj en la formo de vortaro {'type', 'fun', 'jac', 'args'};
3) Analizu ekzemplon de optimumigo de fabrikitaj produktoj uzante la ekzemplon de retejo-studio.
Kondiĉa optimumiga metodo="trust-constr"
Efektivigo de la metodo trust-constr bazita sur por problemoj kun limoj de la formo de egaleco kaj plu por problemoj kun limoj en formo de neegalaĵoj. Ambaŭ metodoj estas efektivigitaj per algoritmoj por trovado de loka minimumo en la fidregiono kaj estas bone konvenitaj por grandskalaj problemoj.
Matematika formuliĝo de la problemo de trovado de minimumo en ĝenerala formo:
Por striktaj egaleclimoj, la malsupra limo estas metita egala al la supra limo .
Por unudirekta limo, la supra aŭ malsupra limo estas fiksita np.inf kun la responda signo.
Estu necese trovi la minimumon de konata Rosenbrock-funkcio de du variabloj:
En ĉi tiu kazo, la sekvaj restriktoj estas fiksitaj sur ĝia domajno de difino:
En nia kazo, estas unika solvo ĉe la punkto , por kiuj validas nur la unua kaj kvara limigoj.
Ni trarigardu la limigojn de malsupre al supro kaj rigardu kiel ni povas skribi ilin en scipy.
Restriktoj и ni difinu ĝin uzante la objekton Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Restriktoj и Ni skribu ĝin en linia formo:
Ni difinu ĉi tiujn limojn kiel LinearConstraint-objekton:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])Kaj finfine la nelinia limo en matrica formo:
Ni difinas la jakobian matricon por ĉi tiu limo kaj lineara kombinaĵo de la hesia matrico kun arbitra vektoro :
Nun ni povas difini nelinian limon kiel objekton 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)Se la grandeco estas granda, matricoj ankaŭ povas esti precizigitaj en malabunda formo:
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)aŭ kiel objekto 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)Kiam oni kalkulas la hesan matricon postulas multan penon, vi povas uzi klason . La sekvaj strategioj haveblas: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())La Hesio ankaŭ povas esti kalkulita uzante finhavajn diferencojn:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')La jakobia matrico por limoj ankaŭ povas esti kalkulita uzante finhavajn diferencojn. Tamen, en ĉi tiu kazo la Hessa matrico ne povas esti kalkulita uzante finhavajn diferencojn. La Hessian devas esti difinita kiel funkcio aŭ uzante la HessianUpdateStrategy klason.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())La solvo al la optimumiga problemo aspektas jene:
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]Se necese, la funkcio por kalkuli la Hessian povas esti difinita uzante la klason 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)aŭ la produkto de la Hessian kaj arbitra vektoro tra la parametro 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)Alternative, la unua kaj dua derivaĵoj de la funkcio estanta optimumigita povas esti proksimumataj. Ekzemple, la Hessian povas esti proksimuma uzante la funkcion SR1 (kvazaŭ-newtona aproksimado). La gradiento povas esti proksimuma per finhavaj diferencoj.
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)Kondiĉa optimumiga metodo="SLSQP"
La SLSQP-metodo estas dizajnita por solvi problemojn de minimumigado de funkcio en la formo:
Kie и — aroj de indeksoj de esprimoj priskribantaj restriktojn en formo de egalecoj aŭ neegalaĵoj. — aroj de malsuperaj kaj superaj baroj por la domajno de difino de la funkcio.
Liniaj kaj neliniaj limoj estas priskribitaj en la formo de vortaroj kun klavoj 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])
}La serĉo de la minimumo estas farita jene:
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 ]Optimumigo Ekzemplo
Lige kun la transiro al la kvina teknologia strukturo, ni rigardu produktadon-optimumigo uzante la ekzemplon de retejo-studio, kiu alportas al ni malgrandan sed stabilan enspezon. Ni imagu nin kiel la direktoro de galero kiu produktas tri specojn de produktoj:
- x0 - vendado de landpaĝoj, de 10 tr.
- x1 - kompaniaj retejoj, ekde 20 tr.
- x2 - interretaj vendejoj, de 30 tr.
Nia amika laborteamo inkluzivas kvar junulojn, du mezajn kaj unu maljunulojn. Ilia monata labortempa fonduso:
- junioj:
4 * 150 = 600 чел * час, - mezoj:
2 * 150 = 300 чел * час, - sinjoro:
150 чел * час.
Lasu la unuajn disponeblajn junulojn pasigi (0, 1, 2) horojn por la disvolviĝo kaj disvastigo de unu loko de tipo (x10, x20, x30), meza - (7, 15, 20), altranga - (5, 10, 15). ) horoj de la plej bona tempo de via vivo.
Kiel ĉiu normala direktoro, ni volas maksimumigi monatajn profitojn. La unua paŝo al sukceso estas noti la objektivan funkcion value kiel la kvanto de enspezo de produktoj produktitaj monate:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Ĉi tio ne estas eraro; serĉante la maksimumon, la objektiva funkcio estas minimumigita kun la kontraŭa signo.
La sekva paŝo estas malpermesi al niaj dungitoj trolabori kaj enkonduki limigojn pri laborhoroj:
Kio estas ekvivalenta:
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]])
}Formala restrikto estas ke produktoproduktaĵo devas nur esti pozitiva:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])Kaj finfine, la plej roza supozo estas, ke pro la malalta prezo kaj alta kvalito, vico da kontentaj klientoj konstante viciĝas por ni. Ni mem povas elekti la monatajn produktadajn volumojn, surbaze de solvado de la limigita optimumiga problemo 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]Ni rondiru malfikse al tutaj nombroj kaj kalkulu la monatan ŝarĝon de remistoj kun optimuma distribuado de produktoj. x = (8, 6, 3) :
- junioj:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - mezoj:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - sinjoro:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Konkludo: por ke la direktoro ricevu sian merititan maksimumon, estas optimume krei 8 landpaĝojn, 6 mezgrandajn retejojn kaj 3 vendejojn monate. En ĉi tiu kazo, la maljunulo devas plugi sen rigardi supren de la maŝino, la ŝarĝo de la mezoj estos proksimume 2/3, la junuloj malpli ol duono.
konkludo
La artikolo skizas la bazajn teknikojn por labori kun la pakaĵo scipy.optimize, uzata por solvi kondiĉajn minimumigproblemojn. Persone mi uzas scipy pure por akademiaj celoj, tial la ekzemplo donita estas de tia komika naturo.
Multaj teorioj kaj virtualaj ekzemploj troveblas, ekzemple, en la libro de I.L. Akulich "Matematika programado en ekzemploj kaj problemoj". Pli hardcore aplikaĵo scipy.optimize konstrui 3D strukturon el aro da bildoj () videblas en .
La ĉefa fonto de informo estas tiuj, kiuj volas kontribui al la traduko de ĉi tiu kaj aliaj sekcioj scipy Bonvenon al .
Спасибо por partopreno en la preparado de la eldonaĵo.
fonto: www.habr.com