SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

Ang SciPy (binibigkas na sai pie) ay isang numpy-based mathematics package na kasama rin ang mga C at Fortran na aklatan. Ginagawa ng SciPy ang iyong interactive na Python session sa isang kumpletong kapaligiran ng agham ng data tulad ng MATLAB, IDL, Octave, R, o SciLab.

Sa artikulong ito, titingnan natin ang mga pangunahing pamamaraan ng mathematical programming - paglutas ng mga problema sa conditional optimization para sa isang scalar function ng ilang variable gamit ang scipy.optimize package. Napag-usapan na ang mga unconstrained optimization algorithm huling artikulo. Ang mas detalyado at up-to-date na tulong sa mga scipy function ay palaging makukuha gamit ang help() command, Shift+Tab o sa opisyal na dokumentasyon.

Pagpapakilala

Ang isang karaniwang interface para sa paglutas ng parehong kondisyon at walang limitasyong mga problema sa pag-optimize sa scipy.optimize package ay ibinibigay ng function. minimize(). Gayunpaman, alam na walang unibersal na pamamaraan para sa paglutas ng lahat ng mga problema, kaya ang pagpili ng isang sapat na pamamaraan, gaya ng lagi, ay nahuhulog sa mga balikat ng mananaliksik.
Ang naaangkop na algorithm ng pag-optimize ay tinukoy gamit ang argumento ng function minimize(..., method="").
Para sa kondisyon na pag-optimize ng isang function ng ilang variable, ang mga pagpapatupad ng mga sumusunod na pamamaraan ay magagamit:

  • trust-constr — maghanap ng lokal na minimum sa rehiyon ng kumpiyansa. Artikulo ng Wiki, artikulo sa Habré;
  • SLSQP — sequential quadratic programming na may mga hadlang, Newtonian na pamamaraan para sa paglutas ng Lagrange system. Artikulo ng Wiki.
  • TNC - Pinutol na Newton Constrained, limitadong bilang ng mga pag-ulit, mabuti para sa mga nonlinear na function na may malaking bilang ng mga independent variable. Artikulo ng Wiki.
  • L-BFGS-B — isang paraan mula sa pangkat ng Broyden–Fletcher–Goldfarb–Shanno, na ipinatupad nang may pinababang pagkonsumo ng memorya dahil sa bahagyang pag-load ng mga vector mula sa Hessian matrix. Artikulo ng Wiki, artikulo sa Habré.
  • COBYLA — MARE Constrained Optimization Sa pamamagitan ng Linear Approximation, constrained optimization na may linear approximation (nang walang gradient na pagkalkula). Artikulo ng Wiki.

Depende sa napiling paraan, ang mga kondisyon at paghihigpit para sa paglutas ng problema ay itinakda nang iba:

  • bagay sa klase Bounds para sa mga pamamaraan L-BFGS-B, TNC, SLSQP, trust-constr;
  • ang listahan (min, max) para sa parehong mga pamamaraan L-BFGS-B, TNC, SLSQP, trust-constr;
  • isang bagay o isang listahan ng mga bagay LinearConstraint, NonlinearConstraint para sa COBYLA, SLSQP, mga pamamaraan ng trust-constr;
  • diksyunaryo o listahan ng mga diksyunaryo {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} para sa COBYLA, mga pamamaraan ng SLSQP.

Balangkas ng artikulo:
1) Isaalang-alang ang paggamit ng conditional optimization algorithm sa trust region (method=”trust-constr”) na may mga hadlang na tinukoy bilang mga object Bounds, LinearConstraint, NonlinearConstraint ;
2) Isaalang-alang ang sequential programming gamit ang least squares method (paraan = "SLSQP") na may mga paghihigpit na tinukoy sa anyo ng isang diksyunaryo {'type', 'fun', 'jac', 'args'};
3) Suriin ang isang halimbawa ng pag-optimize ng mga ginawang produkto gamit ang halimbawa ng isang web studio.

Conditional optimization method="trust-constr"

Pagpapatupad ng pamamaraan trust-constr batay sa EQSQP para sa mga problema sa mga hadlang ng anyo ng pagkakapantay-pantay at sa TRIP para sa mga problema sa mga hadlang sa anyo ng mga hindi pagkakapantay-pantay. Ang parehong mga pamamaraan ay ipinatupad ng mga algorithm para sa paghahanap ng lokal na minimum sa rehiyon ng kumpiyansa at angkop na angkop para sa malalaking problema.

Ang pagbabalangkas ng matematika ng problema ng paghahanap ng isang minimum sa pangkalahatang anyo:

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

Para sa mahigpit na mga hadlang sa pagkakapantay-pantay, ang lower bound ay itinakda na katumbas ng upper bound SciPy, pag-optimize na may mga kundisyon.
Para sa isang one-way na paghihigpit, ang itaas o mas mababang limitasyon ay itinakda np.inf na may kaukulang tanda.
Hayaang kailanganin upang mahanap ang minimum ng isang kilalang Rosenbrock function ng dalawang variable:

SciPy, pag-optimize na may mga kundisyon

Sa kasong ito, ang mga sumusunod na paghihigpit ay nakatakda sa domain ng kahulugan nito:

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

Sa aming kaso, mayroong isang natatanging solusyon sa punto SciPy, pag-optimize na may mga kundisyon, kung saan ang una at ikaapat na paghihigpit lamang ang may bisa.
Dumaan tayo sa mga paghihigpit mula sa ibaba hanggang sa itaas at tingnan kung paano natin ito maisusulat sa scipy.
Paghihigpit SciPy, pag-optimize na may mga kundisyon и SciPy, pag-optimize na may mga kundisyon tukuyin natin ito gamit ang Bounds object.

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

Paghihigpit SciPy, pag-optimize na may mga kundisyon и SciPy, pag-optimize na may mga kundisyon Isulat natin ito sa linear form:

SciPy, pag-optimize na may mga kundisyon

Tukuyin natin ang mga hadlang na ito bilang isang bagay na LinearConstraint:

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

At panghuli ang nonlinear constraint sa matrix form:

SciPy, pag-optimize na may mga kundisyon

Tinukoy namin ang Jacobian matrix para sa hadlang na ito at isang linear na kumbinasyon ng Hessian matrix na may arbitrary na vector SciPy, pag-optimize na may mga kundisyon:

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

Ngayon ay maaari nating tukuyin ang isang nonlinear na pagpilit bilang isang bagay 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)

Kung malaki ang sukat, maaari ding tukuyin ang mga matrice sa kalat-kalat na anyo:

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)

o bilang isang bagay 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)

Kapag kinakalkula ang Hessian matrix SciPy, pag-optimize na may mga kundisyon nangangailangan ng maraming pagsisikap, maaari kang gumamit ng isang klase HessianUpdateStrategy. Ang mga sumusunod na diskarte ay magagamit: BFGS и SR1.

from scipy.optimize import BFGS

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

Ang Hessian ay maaari ding kalkulahin gamit ang mga may hangganang pagkakaiba:

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

Ang Jacobian matrix para sa mga hadlang ay maaari ding kalkulahin gamit ang mga may hangganang pagkakaiba. Gayunpaman, sa kasong ito ang Hessian matrix ay hindi maaaring kalkulahin gamit ang mga may hangganang pagkakaiba. Dapat tukuyin ang Hessian bilang isang function o gamit ang klase ng HessianUpdateStrategy.

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

Ang solusyon sa problema sa pag-optimize ay ganito ang hitsura:

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]

Kung kinakailangan, ang function para sa pagkalkula ng Hessian ay maaaring tukuyin gamit ang LinearOperator class

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)

o ang produkto ng Hessian at isang di-makatwirang vector sa pamamagitan ng parameter 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)

Bilang kahalili, ang una at pangalawang derivative ng function na ino-optimize ay maaaring tantiyahin. Halimbawa, ang Hessian ay maaaring tantiyahin gamit ang function SR1 (quasi-Newtonian approximation). Ang gradient ay maaaring tinantya ng may hangganang pagkakaiba.

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)

Conditional optimization method="SLSQP"

Ang pamamaraan ng SLSQP ay idinisenyo upang malutas ang mga problema sa pagliit ng isang function sa anyo:

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

SciPy, pag-optimize na may mga kundisyon

saan SciPy, pag-optimize na may mga kundisyon и SciPy, pag-optimize na may mga kundisyon — mga hanay ng mga indeks ng mga expression na naglalarawan ng mga paghihigpit sa anyo ng mga pagkakapantay-pantay o hindi pagkakapantay-pantay. SciPy, pag-optimize na may mga kundisyon — set ng lower at upper bounds para sa domain ng definition ng function.

Ang mga linear at nonlinear na hadlang ay inilalarawan sa anyo ng mga diksyunaryo na may mga susi 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])
          }

Ang paghahanap para sa minimum ay isinasagawa tulad ng sumusunod:

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 ]

Halimbawa ng Optimization

Kaugnay ng paglipat sa ikalimang teknolohikal na istraktura, tingnan natin ang pag-optimize ng produksyon gamit ang halimbawa ng isang web studio, na nagdudulot sa atin ng maliit ngunit matatag na kita. Isipin natin ang ating sarili bilang direktor ng isang galley na gumagawa ng tatlong uri ng mga produkto:

  • x0 - nagbebenta ng mga landing page, mula 10 tr.
  • x1 - mga website ng kumpanya, mula 20 tr.
  • x2 - mga online na tindahan, mula 30 tr.

Kasama sa aming magiliw na pangkat sa pagtatrabaho ang apat na junior, dalawang middle at isang senior. Ang kanilang buwanang pondo sa oras ng pagtatrabaho:

  • Hunyo: 4 * 150 = 600 чел * час,
  • gitna: 2 * 150 = 300 чел * час,
  • senor: 150 чел * час.

Hayaang gumastos ang unang available na junior (0, 1, 2) na oras sa pagbuo at pag-deploy ng isang site ng uri (x10, x20, x30), gitna - (7, 15, 20), senior - (5, 10, 15 ) oras ng pinakamagandang oras ng iyong buhay.

Tulad ng anumang normal na direktor, gusto naming i-maximize ang buwanang kita. Ang unang hakbang sa tagumpay ay isulat ang layunin ng function value bilang ang halaga ng kita mula sa mga produktong ginawa bawat buwan:

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

Ito ay hindi isang error; kapag naghahanap para sa maximum, ang layunin ng function ay minimize sa kabaligtaran sign.

Ang susunod na hakbang ay upang ipagbawal ang aming mga empleyado sa labis na pagtatrabaho at ipakilala ang mga paghihigpit sa mga oras ng pagtatrabaho:

SciPy, pag-optimize na may mga kundisyon

Ano ang katumbas:

SciPy, pag-optimize na may mga kundisyon

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

Ang isang pormal na paghihigpit ay ang output ng produkto ay dapat na positibo lamang:

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

At sa wakas, ang pinaka-rosas na palagay ay dahil sa mababang presyo at mataas na kalidad, isang pila ng mga nasisiyahang customer ang patuloy na pumipila para sa amin. Maaari naming piliin ang mga buwanang dami ng produksyon sa aming sarili, batay sa paglutas sa napipigilan na problema sa pag-optimize 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]

Mag-ikot tayo nang maluwag sa mga buong numero at kalkulahin ang buwanang pagkarga ng mga rowers na may pinakamainam na pamamahagi ng mga produkto x = (8, 6, 3) :

  • Hunyo: 8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час;
  • gitna: 8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час;
  • senor: 8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.

Konklusyon: upang matanggap ng direktor ang kanyang karapat-dapat na maximum, pinakamainam na lumikha ng 8 landing page, 6 na medium-sized na site at 3 tindahan bawat buwan. Sa kasong ito, ang senior ay dapat mag-araro nang hindi tumitingin mula sa makina, ang load ng middles ay humigit-kumulang 2/3, ang juniors ay mas mababa sa kalahati.

Konklusyon

Binabalangkas ng artikulo ang mga pangunahing pamamaraan para sa pagtatrabaho sa package scipy.optimize, ginagamit upang malutas ang mga problema sa conditional minimization. Personal na ginagamit ko scipy para lamang sa mga layuning pang-akademiko, kaya naman ang ibinigay na halimbawa ay may likas na komiks.

Maraming teorya at virtual na mga halimbawa ang matatagpuan, halimbawa, sa aklat ni I.L. Akulich "Pagprograma ng matematika sa mga halimbawa at problema." Higit pang hardcore application scipy.optimize upang bumuo ng isang 3D na istraktura mula sa isang hanay ng mga larawan (artikulo sa Habré) ay maaaring matingnan sa scipy-cookbook.

Ang pangunahing mapagkukunan ng impormasyon ay docs.scipy.orgmga nagnanais na mag-ambag sa pagsasalin nito at sa iba pang mga seksyon scipy Maligayang pagdating sa GitHub.

salamat mephistophees para sa pakikilahok sa paghahanda ng publikasyon.

Pinagmulan: www.habr.com

Magdagdag ng komento