SciPy (wymawiane sai pie) to pakiet matematyczny oparty na liczbach, który zawiera także biblioteki C i Fortran. SciPy zamienia interaktywną sesję Pythona w kompletne środowisko do nauki o danych, takie jak MATLAB, IDL, Octave, R lub SciLab.
W tym artykule przyjrzymy się podstawowym technikom programowania matematycznego - rozwiązywaniu problemów optymalizacji warunkowej dla funkcji skalarnej kilku zmiennych przy użyciu pakietu scipy.optimize. Algorytmy optymalizacji nieograniczonej zostały już omówione w . Bardziej szczegółową i aktualną pomoc dotyczącą funkcji scipy można zawsze uzyskać za pomocą polecenia help(), Shift+Tab lub w .
Wprowadzenie
Wspólny interfejs do rozwiązywania zarówno warunkowych, jak i nieograniczonych problemów optymalizacyjnych w pakiecie scipy.optimize zapewnia funkcja minimize(). Wiadomo jednak, że nie ma uniwersalnej metody rozwiązania wszystkich problemów, zatem wybór odpowiedniej metody jak zawsze spada na barki badacza.
Odpowiedni algorytm optymalizacji określa się za pomocą argumentu funkcji minimize(..., method="").
Do warunkowej optymalizacji funkcji kilku zmiennych dostępne są implementacje następujących metod:
trust-constr— poszukiwanie minimum lokalnego w obszarze ufności. , ;SLSQP— sekwencyjne programowanie kwadratowe z ograniczeniami, metoda newtonowska rozwiązywania układu Lagrange’a. .TNC- Obcięty Newton Ograniczona, ograniczona liczba iteracji, dobra dla funkcji nieliniowych z dużą liczbą zmiennych niezależnych. .L-BFGS-B— metoda zespołu Broydena–Fletchera–Goldfarba–Shanno, zaimplementowana przy zmniejszonym zużyciu pamięci w wyniku częściowego ładowania wektorów z macierzy Hessego. , .COBYLA— Optymalizacja z ograniczeniami MARE metodą aproksymacji liniowej, optymalizacja z ograniczeniami za pomocą aproksymacji liniowej (bez obliczania gradientu). .
W zależności od wybranej metody warunki i ograniczenia rozwiązania problemu są ustalane inaczej:
- obiekt klasy
Boundsdla metod L-BFGS-B, TNC, SLSQP, trust-constr; - Lista
(min, max)dla tych samych metod L-BFGS-B, TNC, SLSQP, trust-constr; - obiekt lub listę obiektów
LinearConstraint,NonlinearConstraintdla metod COBYLA, SLSQP, trust-constr; - słownik lub listę słowników
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}dla metod COBYLA, SLSQP.
Zarys artykułu:
1) Rozważ zastosowanie algorytmu optymalizacji warunkowej w obszarze zaufania (method=”trust-constr”) z ograniczeniami określonymi jako obiekty Bounds, LinearConstraint, NonlinearConstraint ;
2) Rozważ programowanie sekwencyjne metodą najmniejszych kwadratów (metoda = „SLSQP”) z ograniczeniami określonymi w formie słownika {'type', 'fun', 'jac', 'args'};
3) Przeanalizuj przykład optymalizacji wytwarzanych produktów na przykładzie studia internetowego.
Metoda optymalizacji warunkowej="trust-constr"
Implementacja metody trust-constr oparte na dla problemów z ograniczeniami postaci równości i dalej dla problemów z ograniczeniami w postaci nierówności. Obie metody są implementowane przez algorytmy znajdowania minimum lokalnego w obszarze ufności i dobrze nadają się do problemów o dużej skali.
Matematyczne sformułowanie problemu znalezienia minimum w postaci ogólnej:
W przypadku ścisłych ograniczeń równości dolna granica jest równa górnej granicy .
W przypadku ograniczenia jednokierunkowego ustawiana jest górna lub dolna granica np.inf z odpowiednim znakiem.
Niech będzie konieczne znalezienie minimum znanej funkcji Rosenbrocka dwóch zmiennych:
W tym przypadku na jego dziedzinę definicji nakładane są następujące ograniczenia:
W naszym przypadku mamy do czynienia z unikalnym rozwiązaniem , dla których obowiązują tylko pierwsze i czwarte ograniczenie.
Przejrzyjmy ograniczenia od dołu do góry i zobaczmy, jak możemy je zapisać w scipy.
Ograniczenia и zdefiniujmy to za pomocą obiektu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Ograniczenia и Zapiszmy to w postaci liniowej:
Zdefiniujmy te ograniczenia jako obiekt LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])I wreszcie ograniczenie nieliniowe w postaci macierzowej:
Definiujemy macierz Jakobiana dla tego ograniczenia oraz liniową kombinację macierzy Hessego z dowolnym wektorem :
Teraz możemy zdefiniować wiązanie nieliniowe jako obiekt 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)Jeśli rozmiar jest duży, macierze można również określić w postaci rzadkiej:
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)lub jako przedmiot 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)Przy obliczaniu macierzy Hessego wymaga dużo wysiłku, możesz skorzystać z zajęć . Dostępne są następujące strategie: BFGS и SR1.
from scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hesjan można również obliczyć za pomocą różnic skończonych:
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = cons_J, hess = '2-point')Macierz Jakobianu dla ograniczeń można również obliczyć przy użyciu różnic skończonych. Jednak w tym przypadku macierzy Hessego nie można obliczyć przy użyciu różnic skończonych. Wartość Hessian należy zdefiniować jako funkcję lub przy użyciu klasy HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Rozwiązanie problemu optymalizacyjnego wygląda następująco:
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]W razie potrzeby funkcję obliczania hesja można zdefiniować za pomocą klasy 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)lub iloczyn Hesja i dowolnego wektora poprzez parametr 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)Alternatywnie można aproksymować pierwszą i drugą pochodną optymalizowanej funkcji. Na przykład Hesjan można aproksymować za pomocą tej funkcji SR1 (przybliżenie quasi-newtonowskie). Gradient można aproksymować za pomocą różnic skończonych.
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)Metoda optymalizacji warunkowej="SLSQP"
Metoda SLSQP przeznaczona jest do rozwiązywania problemów minimalizacji funkcji w postaci:
Где и — zestawy wskaźników wyrażeń opisujących ograniczenia w postaci równości lub nierówności. — zbiory dolnych i górnych granic dziedziny definicji funkcji.
Więzy liniowe i nieliniowe opisane są w formie słowników z kluczami 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])
}Poszukiwanie minimum odbywa się w następujący sposób:
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 ]Przykład optymalizacji
W związku z przejściem na piątą strukturę technologiczną spójrzmy na optymalizację produkcji na przykładzie studia internetowego, które przynosi nam niewielki, ale stabilny dochód. Wyobraźmy sobie siebie jako dyrektora kuchni, która produkuje trzy rodzaje produktów:
- x0 - sprzedaż stron docelowych, od 10 tr.
- x1 - strony korporacyjne, od 20 tr.
- x2 - sklepy internetowe, od 30 tr.
Nasz przyjazny zespół roboczy składa się z czterech juniorów, dwóch środkowych i jednego seniora. Ich miesięczny fundusz czasu pracy:
- Czerwiec:
4 * 150 = 600 чел * час, - środki:
2 * 150 = 300 чел * час, - senor:
150 чел * час.
Niech pierwszy dostępny junior spędzi (0, 1, 2) godzin na opracowywaniu i wdrażaniu jednej witryny typu (x10, x20, x30), środkowej - (7, 15, 20), starszej - (5, 10, 15 ) godziny najlepszego okresu w Twoim życiu.
Jak każdy normalny dyrektor chcemy maksymalizować miesięczne zyski. Pierwszym krokiem do sukcesu jest zapisanie funkcji celu value jako wysokość dochodu z wytworzonych produktów miesięcznie:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]Nie jest to błąd, szukając maksimum, minimalizuje się funkcję celu ze znakiem przeciwnym.
Kolejnym krokiem jest zakazanie naszym pracownikom przepracowania i wprowadzenie ograniczeń w zakresie godzin pracy:
Co jest równoważne:
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]])
}Formalne ograniczenie polega na tym, że wynik produktu może być tylko dodatni:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])I na koniec najbardziej różowe założenie jest takie, że ze względu na niską cenę i wysoką jakość, nieustannie ustawia się do nas kolejka zadowolonych klientów. Sami możemy wybrać miesięczne wolumeny produkcji w oparciu o rozwiązanie problemu optymalizacji z ograniczeniami 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]Zaokrąglijmy luźno do liczb całkowitych i obliczmy miesięczne obciążenie wioślarzy przy optymalnym rozmieszczeniu produktów x = (8, 6, 3) :
- Czerwiec:
8 * 10 + 6 * 20 + 3 * 30 = 290 чел * час; - środki:
8 * 7 + 6 * 15 + 3 * 20 = 206 чел * час; - senor:
8 * 5 + 6 * 10 + 3 * 15 = 145 чел * час.
Wniosek: aby reżyser otrzymał zasłużone maksimum, optymalnie jest stworzyć miesięcznie 8 stron docelowych, 6 witryn średniej wielkości i 3 sklepy. W takim przypadku senior musi orać, nie odrywając wzroku od maszyny, obciążenie środkowych wyniesie około 2/3, juniorów mniej niż połowę.
wniosek
W artykule omówiono podstawowe techniki pracy z pakietem scipy.optimize, używany do rozwiązywania problemów minimalizacji warunkowej. Osobiście używam scipy wyłącznie do celów akademickich i dlatego podany przykład ma tak komiczny charakter.
Dużo teorii i wirtualnych przykładów można znaleźć chociażby w książce I.L. Akulicha „Programowanie matematyczne w przykładach i problemach”. Bardziej hardcorowa aplikacja scipy.optimize zbudować strukturę 3D z zestawu obrazów () można obejrzeć w .
Głównym źródłem informacji jest tych, którzy chcą przyczynić się do tłumaczenia tej i innych sekcji scipy Witamy w .
Dzięki za udział w przygotowaniu publikacji.
Źródło: www.habr.com