SciPy, optimeerimine tingimustega

SciPy (hääldatakse nagu sai pai) on numpy'le põhinev matemaatikapakett, mis sisaldab ka C- ja Fortran-raamatukogusid. SciPy muudab interaktiivse Python seansi täielikuks andmetöötluskeskkonnaks, nagu MATLAB, IDL, Octave, R või SciLab.

Selles artiklis vaatleme matemaatilise programmeerimise põhimeetodeid — tingimustega optimeerimise probleemide lahendamine mitme muutuja skalaarsest funktsioonist scipy.optimize paketi abil. Ilma tingimusteta optimeerimise algoritmid on juba käsitletud eelmisel artiklil. Üksikasjalikku ja ajakohast teavet funktsioonide kohta scipy's saab alati kasutada käsu help(), Shift+Tab või ametlikust dokumentatsioonist.

Sissejuhatus

Üks ühtne liides tingimustega ja ilma tingimusteta optimeerimise probleemide lahendamiseks scipy.optimize pakendis antakse funktsiooni kaudu minimize(). Siiski on teada, et universaalset viisi kõigi probleemide lahendamiseks ei eksisteeri, seega langeb sobiva meetodi valik alati uurija õlgadele.
Sobiv optimeerimisalgoritm määratakse funktsiooni argumendi kaudu minimize(..., method="").
Mitme muutuja tingimustega optimeerimise jaoks on saadaval järgmiste meetodite rakendused:

• trust-constr — lokaalse miinimumi otsimine usaldusväärsuspiirkonnas. Artikkel vikis, artikkel habras;
• SLSQP — järjestikune kvadratiivne programmeerimine piirangutega, Newtoni meetod Lagrange'i süsteemi lahendamiseks. Artikkel vikis.
• TNC — Truncated Newton Constrained, piiratud iteratsioonide arv, sobiv mitteeline funktsioonide jaoks, millel on palju sõltumatuid muutujaid. Artikkel vikis.
• L-BFGS-B — Broyden–Fletcher–Goldfarb–Shanno meetod, mis on rakendatud vähendatud mälu tarbimisega, osalise Hess'i maatriksi vektorite laadimise kaudu. Artikkel vikis, artikkel habras.
• COBYLA — KOBYLA Constrained Optimization By Linear Approximation, piiratud optimeerimine lineaarsest lähenemisest (ilma gradientide arvutamiseta). Artikkel vikis.

Sõltuvalt valitud meetodist määratakse tingimused ja piirangud probleemide lahendamiseks erinevalt:

• klassist objekti Bounds meetodite jaoks L-BFGS-B, TNC, SLSQP, trust-constr;
• loendi (min, max) nende sama meetodite L-BFGS-B, TNC, SLSQP, trust-constr jaoks;
• objekti või objektide loend LinearConstraint, NonlinearConstraint meetodite jaoks COBYLA, SLSQP, trust-constr;
• sõnastiku või sõnastike loendi {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} meetodite jaoks COBYLA, SLSQP.

Artikli plaan:
1) Uurida tingimuslikku optimeerimisalgoritmi usaldusväärses valdkonnas (method=»trust-constr») piirangutega, mis on määratletud objektidena. Bounds, LinearConstraint, NonlinearConstraint ;
2) Uurida järjestikust programmeerimist vähenenud ruutude meetodil (method=»SLSQP») piirangutega, mis on määratletud sõnastikuna. {'type', 'fun', 'jac', 'args'};
3) Analüüsida toodangu optimeerimise näidet veebistuudio näitel.

Tingimuslik optimeerimine method=»trust-constr»

Meetodi rakendamine trust-constr on põhinenud EQSQP võrdsete piirangutega probleemidele ja TRIP ebaühtsete piirangute probleemidele. Mõlemad meetodid on ellu viidud usaldusväärse piirkonna kohaliku miinimumi otsimise algoritmidega ja sobivad hästi suurte probleemide lahendamiseks.

Probleemi matemaatiline vorm üldiselt minimumi leidmiseks:

Tehete range võrdsuse puhul seotakse alumine piir ülemise külge. .
Ühepoolse piirangu korral seatakse ülemine või alumine piir np.inf vastava märgiga.
Oletame, et tuleb leida miinimum tuntud Rosenbrocki funktsioonist, mis sõltub kahest muutuja:

Sellel on määratud järgmised piirangud selle määramispiirkonnale:

Meie juhul on olemas ainus lahendus punktis , kus kehtivad ainult esimene ja neljas piirang.
Liigume piirangute kaudu alt üles ja uurime, kuidas neid scipy-s üles kirjutada.
Piirangud ja määra objekti Bounds abil.

from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])

Piirangud ja kirjutame lineaarse vormina:

Määrame need piirangud objekti LinearConstraint kujul:

import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])

Ja lõpuks mitte-lineaarne piirang maatriksivormis:

Määrame selle piirangu jaoks Jakobian ja Hesse maatriksi lineaarse kombinatsiooni suvalise vektoriga. :

Nüüd saame mitte-lineaarse piirangu määrata objekti 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)

Kui suurus on suur, saab maatriseid määrata ka hõredas vormis:

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)

või objekti LinearOperator:

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)

Kui Hessian-matriisi arvutamine nõuab suurt ressursikasutust, saab kasutada klassi HessianUpdateStrategy. Saadaval on järgmised strateegiad: BFGS ja SR1.

scipy.optimize import BFGS

nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())

Hessianti saab samuti arvutada lõplike erinevuste abil:

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')

Piirangute Jakobianit saab samuti arvutada lõplike erinevuste abil. Kuid sel juhul ei saa Hessiani enam lõplike erinevuste abil arvutada. Hessian peab olema määratletud funktsioonina või läbi klassi HessianUpdateStrategy.

nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())

Optimeerimisprobleemi lahendus on järgmine:

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` lõpetamise tingimus on täidetud.
Iteratsioonide arv: 12, funktsioonide hindamised: 8, CG iteratsioonid: 7, optimaalne: 2.99e-09, piirangute rikkumine: 1.11e-16, täitmise aeg: 0.033 s.
[0.41494531 0.17010937]

Vajadusel saab Hessiani arvutamise funktsiooni määrata klassi LinearOperator abil

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)

või Hessiani ja suvalise vektori korrutamine parametri 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)

Alternatiivselt saab optimeeritava funktsiooni esimesi ja teisi tuletusi ligikaudselt arvutada. Näiteks saab Hessiani approximiseerida funktsiooni SR1 (kvasi-Newtoni lähenemine). Gradienti saab ligikaudselt arvutada lõplike erinevustega.

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)

Konventsionaalne optimeerimine method=»SLSQP»

SLSQP meetod on mõeldud funktsiooni minimeerimise probleemide lahendamiseks kujul:

Kus ja — hulk indekseid, mis kirjeldavad piiranguid võrduste või ebavõrdsuste kujul. — hulk alumisi ja ülemisi piire funktsiooni määramispiirkonnale.

Lineaarseid ja mitte-lineaarseid piiranguid kirjeldatakse sõnaraamatutena, millel on võtmed type, lõbu ja 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])
          }

Minimeerimise protseduur toimub järgmiselt:

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)

Optimeerimine lõpetati edukalt.    (Väljamineku režiim 0)
            Praegune funktsiooni väärtus: 0.34271757499419825
            Iteratsioonid: 4
            Funktsiooni hindamised: 5
            Gradientide hindamised: 4
[0.41494475 0.1701105 ]

Optimeerimise näide

Seoses viienda tehnoloogilise tsükli üleminekuga vaatleme tootmise optimeerimist veebistuudio näitel, mis toob meile väikese, kuid stabiilse tulu. Kujutame ette, et oleme galerii direktor, kus toodetakse kolme tüüpi tooteid:

• x0 — müüvad landindid, alates 10 t. r.
• x1 — ettevõtte veebilehed, alates 20 t. r.
• x2 — internetipoed, alates 30 t. r.

Meie sõbralik töörühm koosneb neljast algajatest, kahest keskastme töötajast ja ühest kogenud spetsialistist. Nende tööaja fond kuus:

• algajad: 4 * 150 = 600 inimene * tund,
• keskastme töötajad: 2 * 150 = 300 inimene * tund,
• kogenud spetsialist: 150 inimene * tund.

Oletame, et ühe tüüpi veebilehe (x0, x1, x2) arendamiseks ja juurutamiseks peab esimene juba tekkiv algaja kulutama (10, 20, 30) tundi, keskastme töötaja — (7, 15, 20), kogenud spetsialist — (5, 10, 15) parima ajaga oma elust.

Nagu igale normaalsele direktorile, soovime maksimeerida igakuist kasumit. Esimene samm edusammudeks — kirjutame sihtfunktsiooni value kui kuu jooksul toodetud toodete tulu summana:

def value(x):
    return - 10*x[0] - 20*x[1] - 30*x[2]

See ei ole viga; sihtfunktsiooni maksimeerimisel minimeeritakse see vastupidise märgiga.

Järgmine samm — keelame oma töötajatel ületunde ja kehtestame piirangud tööaja fondile:

Mis on ekvivalentne:

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]])
            }

Formaalne piirang on see, et toodang peab olema ainult positiivne:

bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])

Ja lõpuks kõige roosilisem eeldus — madala hinna ja kõrge kvaliteedi tõttu seisavad meie juures pidevalt rahulolevate klientide read. Saame ise valida kuukohustused toodangute tootmisel, lähtudes tingimuslikust optimeerimisülesandest. 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]

Ümmardame ligikaudu täisarvudeni ja arvutame kuukoormuse optimiseeritud toodangu korral. x = (8, 6, 3) :

• algajad: 8 * 10 + 6 * 20 + 3 * 30 = 290 inimest * tund;
• keskastme töötajad: 8 * 7 + 6 * 15 + 3 * 20 = 206 inimest * tund;
• kogenud spetsialist: 8 * 5 + 6 * 10 + 3 * 15 = 145 inimest * tund.

Järeldus: et direktor saaks oma teenitud maksimumi, on optimaalne kuus teha 8 maandumist, 6 keskmist veebilehte ja 3 poodi. Seniors peab samal ajal tõeliselt vaeva nägema, loaderite koormus jääb umbes 2/3, juniorite vähem kui pool.

Kokkuvõte

Artiklis on esitatud põhitehnikad paketi scipy.optimize, mida kasutatakse tingimuslikku minimeerimise ülesannete lahendamiseks. Isiklikult kasutan ma scipy puhas akadeemiliste eesmärkide nimel, seega on antud näide naljakas.

Palju teooriat ja haruldasi näiteid saab leida näiteks I.L. Akulichi raamatust "Matemaatiline programmeerimine näidete ja ülesannetega". Kaheldamatum rakendus scipy.optimize 3D struktuuri loomine pildigruppide põhjal (artikkel habras) on nüüd võimalik näha scipy-cookbook.

Peamine teabe allikas on docs.scipy.org, kes soovivad panustada tõlkesesse ja teistesse osadesse, scipy teretulemast GitHub.

Aitäh mephistopheies osalemise eest väljaande ettevalmistamisel.

Allikas: habr.com