SciPy, optimalisering med betingelser

SciPy (произносится как сай пай) — это основанный на numpy математический пакет, включающий в себя также библиотеки на C и Fortran. С SciPy интерактивный сеанс Python превращается в такую же полноценную среду обработки данных, как MATLAB, IDL, Octave, R или SciLab.

I denne artikkelen skal vi se på de grunnleggende teknikkene for matematisk programmering - å løse betingede optimaliseringsproblemer for en skalarfunksjon av flere variabler ved å bruke scipy.optimize-pakken. Ubegrensede optimaliseringsalgoritmer har allerede blitt diskutert i siste artikkel. Mer detaljert og oppdatert hjelp om scipy-funksjoner kan alltid fås ved å bruke help()-kommandoen, Shift+Tab eller i offisiell dokumentasjon.

Innledning

Et felles grensesnitt for å løse både betingede og ubegrensede optimaliseringsproblemer i scipy.optimize-pakken leveres av funksjonen minimize(). Imidlertid er det kjent at det ikke er noen universell metode for å løse alle problemer, så valget av en adekvat metode, som alltid, faller på forskerens skuldre.
Подходящий алгоритм оптимизации задается с помощью аргумента функции minimize(..., method="").
For betinget optimalisering av en funksjon av flere variabler, er implementeringer av følgende metoder tilgjengelige:

• trust-constr — поиск локального минимума в доверительной области. Статья на wiki, artikkel om Habré;
• SLSQP — последовательное квадратичное программирование с ограничениями, ньютоновский метод решения системы Лагранжа. Wiki-artikkel.
• TNC — Truncated Newton Constrained, ограниченное число итераций, хорош для нелинейных функций с большим числом независимых переменных. Статья на wiki.
• L-BFGS-B — en metode fra Broyden–Fletcher–Goldfarb–Shanno-teamet, implementert med redusert minneforbruk på grunn av delvis lasting av vektorer fra den hessiske matrisen. Статья на wiki, artikkel om Habré.
• COBYLA — КОБЫЛА Constrained Optimization By Linear Approximation, ограниченная оптимизация с линейной аппроксимацией (без вычисления градиента). Статья на wiki.

Avhengig av den valgte metoden, er betingelser og restriksjoner for å løse problemet satt annerledes:

• klasseobjekt Bounds для методов L-BFGS-B, TNC, SLSQP, trust-constr;
• listen (min, max) for de samme metodene L-BFGS-B, TNC, SLSQP, trust-constr;
• et objekt eller en liste over objekter LinearConstraint, NonlinearConstraint для методов COBYLA, SLSQP, trust-constr;
• ordbok eller liste over ordbøker {'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt} for COBYLA, SLSQP metoder.

Artikkeloversikt:
1) Рассмотреть применение алгоритма условной оптимизации в доверительной области (method=»trust-constr») с ограничениями, заданными в виде объектов Bounds, LinearConstraint, NonlinearConstraint ;
2) Рассмотреть последовательное программирование методом наименьших квадратов (method=»SLSQP») с ограничениями, заданными в виде словаря {'type', 'fun', 'jac', 'args'};
3) Analyser et eksempel på optimalisering av produserte produkter ved å bruke eksemplet med et webstudio.

Betinget optimaliseringsmetode = "trust-constr"

Implementering av metoden trust-constr basert på EQSQP для задач с ограничениями вида равенства и на TUR for problemer med begrensninger i form av ulikheter. Begge metodene er implementert av algoritmer for å finne et lokalt minimum i konfidensregionen og er godt egnet for store problemer.

Математическая постановка задачи поиска минимума в общем виде:

Для ограничений строгого равенства нижняя граница устанавливается равной верхней .
Для одностороннего ограничения верхняя или нижняя граница устанавливается np.inf с соответствующим знаком.
La det være nødvendig å finne minimum av en kjent Rosenbrock-funksjon av to variabler:

I dette tilfellet er følgende begrensninger satt på definisjonsdomenet:

В нашем случае имеется единственное решение в точке , для которой справедливы только первое и четвертое ограничения.
Пройдемся по ограничениям снизу вверх и рассмотрим, как можно их записать в scipy.
Restriksjoner и определим с помощью объекта Bounds.

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

Restriksjoner и запишем в линейной форме:

Определим эти ограничения в виде объекта LinearConstraint:

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

Og til slutt den ikke-lineære begrensningen i matriseform:

Vi definerer den jakobiske matrisen for denne begrensningen og en lineær kombinasjon av den hessiske matrisen med en vilkårlig vektor :

Теперь нелинейное ограничение можем определить как объект 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)

Если размер велик, матрицы можно задавать и в разреженном виде:

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)

или как объект 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)

Ved beregning av den hessiske matrisen требует больших затрат, можно использовать класс HessianUpdateStrategy. Følgende strategier er tilgjengelige: BFGS и SR1.

from scipy.optimize import BFGS

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

Гессиан также может быть вычислен с помощью конечных разностей:

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

Матрицу Якоби для ограничений также можно вычислить с помощью конечных разностей. Однако, в этом случае матрицу Гессе с помощью конечных разностей уже не вычислить. Гессиан должен быть определен в виде функции или с помощью класса HessianUpdateStrategy.

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

Løsningen på optimaliseringsproblemet ser slik ut:

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]

При необходимости, функцию вычисления гессиана можно определить с помощью класса 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 produktet av hessian og en vilkårlig vektor gjennom parameteren 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)

Альтернативно, первая и вторая производные оптимизируемой функции могут быть вычислены приближенно. Например, гессиан может быть аппроксимирован с помощью функции SR1 (kvasi-newtonsk tilnærming). Gradienten kan tilnærmes ved endelige forskjeller.

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)

Условная оптимизация method=»SLSQP»

SLSQP-metoden er designet for å løse problemer med å minimere en funksjon i formen:

hvor и — sett med indekser av uttrykk som beskriver restriksjoner i form av likheter eller ulikheter. — множества нижних и верхних границ для области определения функции.

Lineære og ikke-lineære begrensninger er beskrevet i form av ordbøker med nøkler 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])
          }

Поиск минимума осуществляется следующим образом:

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 ]

Optimaliseringseksempel

В связи с переходом к пятому технологическому укладу, рассмотрим оптимизацию производства на примере веб-студии, которая приносит нам небольшой, но стабильный доход. Представим себя директором галеры, на которой производится три вида продукции:

• x0 - selge landingssider, fra 10 tr.
• x1 - bedriftsnettsteder, fra 20 tr.
• x2 — интернет магазины, от 30 т.р.

Наш дружный рабочий коллектив включает в себя четырех джунов, двух мидлов и одного сеньора. Фонд их рабочего времени на месяц:

• juni: 4 * 150 = 600 чел * час,
• мидлы: 2 * 150 = 300 чел * час,
• сеньор: 150 чел * час.

La den første tilgjengelige junior bruke (0, 1, 2) timer på utvikling og distribusjon av ett nettsted av typen (x10, x20, x30), mellom - (7, 15, 20), senior - (5, 10, 15) ) timer av den beste tiden i livet ditt.

Som enhver vanlig direktør ønsker vi å maksimere månedlig fortjeneste. Det første trinnet til suksess er å skrive ned den objektive funksjonen value как сумму доходов от произведенной за месяц продукции:

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

Это не ошибка, при поиске максимума целевая функция минимизируется с обратным знаком.

Det neste trinnet er å forby våre ansatte å overarbeide og innføre restriksjoner på arbeidstid:

Что эквивалентно:

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

Формальное ограничение — выпуск продукции должен быть только положительным:

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

Og til slutt, den mest rosenrøde antagelsen er at på grunn av den lave prisen og den høye kvaliteten, står en kø av fornøyde kunder stadig i kø for oss. Vi kan velge de månedlige produksjonsvolumene selv, basert på å løse det begrensede optimaliseringsproblemet 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]

Нестрого округлим до целых и посчитаем месячную загрузку гребцов при оптимальном раскладе продукции x = (8, 6, 3) :

• juni: 8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час;
• мидлы: 8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час;
• сеньор: 8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.

Konklusjon: for at regissøren skal få sitt velfortjente maksimum, er det optimalt å lage 8 landingssider, 6 mellomstore nettsteder og 3 butikker per måned. I dette tilfellet må senioren pløye uten å se opp fra maskinen, belastningen på de midterste vil være ca 2/3, juniorene mindre enn halvparten.

Konklusjon

Artikkelen skisserer de grunnleggende teknikkene for å jobbe med pakken scipy.optimize, используемые для решения задач условной минимизации. Лично я использую scipy rent for akademiske formål, og det er derfor eksemplet som er gitt er av en så komisk karakter.

Mange teorier og virtuelle eksempler finnes for eksempel i boken til I.L. Akulich "Matematisk programmering i eksempler og problemer." Mer hardcore-applikasjon scipy.optimize å bygge en 3D-struktur fra et sett med bilder (artikkel om Habré) можно посмотреть в scipy-cookbook.

Hovedkilden til informasjon er docs.scipy.orgde som ønsker å bidra til oversettelsen av denne og andre seksjoner scipy Velkommen til GitHub.

Takk mephistophees за участие в подготовке публикации.

Kilde: www.habr.com