🥇SciPy, optimizare cu condiții

SciPy (se pronunță ca 'sai pai') este un pachet matematic bazat pe numpy, care include de asemenea biblioteci în C și Fortran. Cu SciPy, o sesiune interactivă Python devine un mediu complet de procesare a datelor, similar cu MATLAB, IDL, Octave, R sau SciLab.

În acest articol vom explora principalele tehnici de programare matematică — soluționarea problemelor de optimizare condiționată pentru o funcție scalară cu mai multe variabile, utilizând pachetul scipy.optimize. Algoritmii de optimizare necondiționată au fost deja examinați în articolul precedent. O documentație mai detaliată și actualizată despre funcțiile scipy este întotdeauna disponibilă folosind comanda help(), Shift+Tab sau în documentația oficială.

Introducere

Interfața generală pentru rezolvarea problemelor atât de optimizare condiționată, cât și necondiționată în pachetul scipy.optimize este oferită de funcția minimize(). Cu toate acestea, se știe că nu există o metodă universală pentru rezolvarea tuturor problemelor, așa că alegerea unei metode adecvate revine, ca de obicei, cercetătorului.
Algoritmul de optimizare potrivit se stabilește prin argumentul funcției minimize(..., method="").
Pentru optimizarea condiționată a unei funcții cu mai multe variabile sunt disponibile implementările următoarelor metode:

trust-constr — căutarea unui minim local în domeniul de încredere. Articol pe wiki, articol pe Habr;
SLSQP — programare cvadratică secvențială cu constrângeri, metoda lui Newton de soluționare a sistemului Lagrange. Articol pe wiki.
TNC — Truncated Newton Constrained, un număr limitat de iterații, bun pentru funcții neliniare cu un număr mare de variabile independente. Articol pe wiki.
L-BFGS-B — metoda din cvartetul Broyden–Fletcher–Goldfarb–Shanno, implementată cu un consum redus de memorie prin încărcarea parțială a vectorilor din matricea Hessian. Articol pe wiki, articol pe Habr.
COBYLA — Constrained Optimization By Linear Approximation, optimizare restricționată prin aproximarea liniară (fără calcularea gradientului). Articol pe wiki.

În funcție de metoda aleasă, condițiile și constrângerile pentru rezolvarea problemei sunt stabilite în mod diferit:

obiectul clasei Bounds pentru metodele L-BFGS-B, TNC, SLSQP, trust-constr;
lista (min, max) pentru aceleași metode L-BFGS-B, TNC, SLSQP, trust-constr;
obiectul sau lista obiectelor LinearConstraint, NonlinearConstraint pentru metodele COBYLA, SLSQP, trust-constr;
un dicționar sau o listă de dicționare {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} pentru metodele COBYLA, SLSQP.

Planul articolului:
1) Considerarea aplicării algoritmului de optimizare condiționată în domeniul de încredere (method=»trust-constr») cu restricții definite sub formă de obiecte Bounds, LinearConstraint, NonlinearConstraint ;
2) Considerarea programării secvențiale prin metoda celor mai mici pătrate (method=»SLSQP») cu restricții definite sub formă de dicționar {'type', 'fun', 'jac', 'args'};
3) Analiza unui exemplu de optimizare a producției utilizând un exemplu de studio web.

Optimizare condiționată method=»trust-constr»

Implementarea metodei trust-constr se bazează pe EQSQP pentru probleme cu restricții de egalitate și pe TRIP pentru probleme cu restricții de tip inegalitate. Ambele metode sunt implementate prin algoritmi de căutare a minimului local în domeniul de încredere și sunt bine adaptate pentru probleme de mari dimensiuni.

Formularea matematică a problemei de căutare a minimului în forma generală:

Pentru restricțiile de egalitate strictă, limita inferioară este stabilită egală cu limita superioară .
Pentru restricția unidirecțională, limita superioară sau inferioară este stabilită np.inf cu semnul corespunzător.
Să presupunem că trebuie să găsim minimul unei funcții cunoscute, funcția Rosenbrock, de două variabile:

În acest caz, sunt specificate următoarele restricții asupra domeniului său de definiție:

În cazul nostru există o singură soluție în punctul , pentru care sunt valabile doar prima și a patra restricție.
Să parcurgem restricțiile de la bază în sus și să vedem cum le putem scrie în scipy.
Limitări și definim cu ajutorul obiectului Bounds.

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

Limitări și să le scriem sub formă liniară:

Definim aceste restricții sub formă de obiect LinearConstraint:

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

Și în sfârșit restricția neliniară în formă matricială:

Definim matricea Jacobi pentru această restricție și combinația liniară a matricei Hessian cu un vector arbitrar :

Acum putem defini restricția neliniară ca un obiect 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)

Dacă dimensiunea este mare, matricele se pot defini și în formă rară:

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)

sau ca obiect LinearOperator:

de la 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)

Când calculul matricei Hessian necesită resurse mari, se poate utiliza clasa HessianUpdateStrategy . Următoarele strategii sunt disponibile:BFGS SR1 și de la scipy.optimize import BFGSnonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS()).

Hessianul poate fi, de asemenea, calculat prin diferențe finite:

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

Matricea Jacobian pentru constrângeri poate fi, de asemenea, calculată prin diferențe finite. Totuși, în acest caz, matricea Hessian nu mai poate fi calculată prin diferențe finite. Hessianul trebuie definit ca o funcție sau prin clasa HessianUpdateStrategy.

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

Soluția problemei de optimizare arată astfel:

de la scipy.optimize import minimize from scipy.optimize import rosen, rosen_der, rosen_hess, rosen_hess_prodx0 = 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` condiția de terminare este satisfăcută.
Numărul de iterații: 12, evaluări de funcție: 8, iterații CG: 7, optimalitate: 2.99e-09, încălcarea constrângerii: 1.11e-16, timpul de execuție: 0.033 s.
[0.41494531 0.17010937]

Dacă este necesar, funcția de calcul a Hessianului poate fi definită prin clasa 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)

sau produsul Hessianului și unui vector arbitrar prin parametrul

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):

Alternativ, primele și cele de-a doua derivate ale funcției optimizate pot fi calculate aproximativ. De exemplu, Hessianul poate fi aproximat printr-o funcție

(aproximarea quasi-Newton). Gradientul poate fi aproximat prin diferențe finite. de la scipy.optimize import BFGSnonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS()) de la 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)

Optimizarea condiționată method=»SLSQP»

Metoda SLSQP este destinată rezolvării problemelor de minimizare a funcției sub formă:

Метод SLSQP предназначен для решения задач минимизации функции в виде:

Unde și — mulțimi de indecși ai expresiilor care descriu constrângeri sub forma de egalități sau inegalități. — mulțimi de limite inferioare și superioare pentru domeniul de definiție al funcției.

Constrângerile liniare și neliniare sunt descrise sub formă de dicționare cu chei type, distracția și 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])
          }

Căutarea minimului se realizează astfel:

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)

Optimizarea s-a încheiat cu succes.    (Cod de ieșire 0)
            Valoarea funcției curente: 0.34271757499419825
            Iterații: 4
            Evaluări ale funcției: 5
            Evaluări ale gradientului: 4
[0.41494475 0.1701105 ]

Exemplu de optimizare

Având în vedere tranziția către a cincea spirală tehnologică, să analizăm optimizarea producției pe exemplul unei agenții web care ne aduce un venit mic, dar stabil. Să ne imaginăm că suntem directorul unei galerii care produce trei tipuri de produse:

x0 — landing page-uri vândute, de la 10 mii RON.
x1 — site-uri corporate, de la 20 mii RON.
x2 — magazine online, de la 30 mii RON.

Echipa noastră unită de muncă este formată din patru juniori, doi mid-level și un senior. Fondul de muncă pentru o lună este:

juniori: 4 * 150 = 600 ore de muncă,
mid-level: 2 * 150 = 300 ore de muncă,
senior: 150 ore de muncă.

Să presupunem că pentru dezvoltarea și implementarea unui site de tip (x0, x1, x2) primul junior întâlnit trebuie să aloce (10, 20, 30) ore, mid-level — (7, 15, 20), senior — (5, 10, 15) ore din cele mai bune zile ale vieții sale.

Ca orice director normal, dorim să maximizăm profitul lunar. Primul pas spre succes este să notăm funcția obiectivă valoare ca suma veniturilor din produsele realizate în luna respectivă:

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

Aceasta nu este o eroare, la căutarea maximului, funcția obiectivă este minimizată cu semnul invers.

Pasul următor — interzicem angajaților să lucreze ore suplimentare și introducem constrângeri asupra fondului de muncă:

Ceea ce este echivalent cu:

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

Limitarea formală - producția trebuie să fie doar pozitivă:

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

Și în cele din urmă, cea mai optimistă presupunere - datorită prețului mic și calității ridicate, avem mereu o coadă de clienți mulțumiți. Putem alege volumele lunare de producție, bazându-ne pe soluția unei probleme de optimizare condiționată cu 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]

Să rotunjim lax la întregi și să calculăm încărcarea lunară a vâslelor în scenariul optim al producției x = (8, 6, 3) :

juniori: 8 * 10 + 6 * 20 + 3 * 30 = 290 persoane * oră;
mid-level: 8 * 7 + 6 * 15 + 3 * 20 = 206 persoane * oră;
senior: 8 * 5 + 6 * 10 + 3 * 15 = 145 persoane * oră.

Concluzie: pentru ca directorul să primească maximul meritat, este optim să facă lunar 8 landing-uri, 6 site-uri medii și 3 magazine. Seniorul trebuie să lucreze neîncetat la stație, încărcarea mijlocie va fi de aproximativ 2/3, iar juniorii sub jumătate.

Concluzie

Articolul prezintă principalele tehnici de lucru cu pachetul scipy.optimize, utilizate pentru a rezolva probleme de minimizare condiționată. Personal folosesc scipy în scopuri pur academice, de aceea exemplul dat are un caracter glumeț.

Multe teorii și exemple raritate pot fi găsite, de exemplu, în cartea lui I.L. Akulich «Programarea matematică în exemple și probleme». O aplicație mai complexă scipy.optimize pentru construirea unei structuri 3D dintr-un set de imagini (articol pe Habr) poate fi consultată în scipy-cookbook.

Principalul sursă de informații - docs.scipy.org, doritorii de a contribui la traducerea acestei și altor secțiuni scipy sunt bineveniți pe GitHub.

Mulțumesc mephistopheies pentru contribuția adusă pregătirii publicației.

Sursa: habr.com