SciPy (wymawia się jak "saj paj") to oparty na numpy pakiet matematyczny, który zawiera również biblioteki w C i Fortranie. Dzięki SciPy interaktywny sesja Python staje się równie pełnoprawnym środowiskiem przetwarzania danych, jak MATLAB, IDL, Octave, R czy SciLab.
W tym artykule omówimy podstawowe metody programowania matematycznego – rozwiązywanie zadań warunkowej optymalizacji dla skalarnej funkcji wielu zmiennych za pomocą pakietu scipy.optimize. Algorytmy bezwarunkowej optymalizacji zostały już omówione w . Szczegółowy i aktualny poradnik dotyczący funkcji scipy zawsze można uzyskać za pomocą polecenia help(), Shift+Tab lub w .
Wprowadzenie
Ogólny interfejs do rozwiązywania zadań zarówno warunkowej, jak i bezwarunkowej optymalizacji w pakiecie scipy.optimize jest zapewniany przez funkcję minimize(). Jednak wiadomo, że nie ma uniwersalnej metody do rozwiązania wszystkich zadań, dlatego wyboru odpowiedniej metody jak zawsze dokonuje badacz.
Odpowiedni algorytm optymalizacji określa się za pomocą argumentu funkcji minimize(..., method="").
Dla warunkowej optymalizacji funkcji wielu zmiennych dostępne są implementacje następujących metod:
trust-constr— poszukiwanie lokalnego minimum w obszarze zaufania. , ;SLSQP— sekwencyjne programowanie kwadratowe z ograniczeniami, metoda Newtona do rozwiązania systemu Lagrange'a. .TNC— Truncated Newton Constrained, ograniczona liczba iteracji, dobra dla nieliniowych funkcji z dużą liczbą zmiennych niezależnych. .L-BFGS-B— metoda z czwórki Broyden–Fletcher–Goldfarb–Shanno, zrealizowana przy zmniejszonym zużyciu pamięci dzięki częściowemu ładowaniu wektorów z macierzy Hessego. , .COBYLA— KOBYLA Constrained Optimization By Linear Approximation, ograniczona optymalizacja z liniową aproksymacją (bez obliczania gradientu). .
W zależności od wybranej metody, warunki i ograniczenia do rozwiązania zadania określane są w różny sposób:
- obiektem klasy
Boundsdla metod L-BFGS-B, TNC, SLSQP, trust-constr; - listą
(min, max)dla tych samych metod L-BFGS-B, TNC, SLSQP, trust-constr; - obiektem lub listą obiektów
LinearConstraint,NonlinearConstraintdla metod COBYLA, SLSQP, trust-constr; - słownikiem lub listą słowników
{'type':str, 'fun':callable, 'jac':callable,opt, 'args':sequence,opt}dla metod COBYLA, SLSQP.
Plan artykułu:
1) Rozważyć zastosowanie algorytmu optymalizacji warunkowej w zakresie zaufania (method=»trust-constr») z ograniczeniami zdefiniowanymi w postaci obiektów Bounds, LinearConstraint, NonlinearConstraint ;
2) Rozważyć sekwencyjne programowanie metodą najmniejszych kwadratów (method=»SLSQP») z ograniczeniami zdefiniowanymi w postaci słownika {'type', 'fun', 'jac', 'args'};
3) Rozpatrzyć przykład optymalizacji produkcji na przykładzie studia webowego.
Optymalizacja warunkowa method=»trust-constr»
Realizacja metody trust-constr opiera się na dla zadań z ograniczeniami w postaci równości oraz na dla zadań z ograniczeniami w postaci nierówności. Oba metody są wdrażane algorytmami poszukiwania lokalnego minimum w zakresie zaufania i dobrze nadają się do problemów o dużej skali.
Matematyczne sformułowanie problemu poszukiwania minimum w ogólnej postaci:
Dla ograniczeń ścisłej równości dolna granica ustalana jest równa górnej 
.
Dla jednostronnego ograniczenia górna lub dolna granica ustalana jest np.inf z odpowiednim znakiem.
Załóżmy, że konieczne jest znalezienie minimum znanej funkcji Rosenbrocka dwóch zmiennych:
Przy tym zadane są następujące ograniczenia na jej dziedzinę:
W naszym przypadku istnieje jedno rozwiązanie w punkcie 
, dla którego spełnione są tylko pierwsze i czwarte ograniczenia.
Przejrzymy ograniczenia od dołu do góry i rozważymy, jak można je zapisać w scipy.
Ograniczenia 
i 
zdefiniujemy za pomocą obiektu Bounds.
from scipy.optimize import Bounds
bounds = Bounds ([0, -0.5], [1.0, 2.0])Ograniczenia 
i 
zapiszemy w postaci liniowej:
Zdefiniujmy te ograniczenia w postaci obiektu LinearConstraint:
import numpy as np
from scipy.optimize import LinearConstraint
linear_constraint = LinearConstraint ([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])A na koniec nieliniowe ograniczenie w postaci macierzy:
Zdefiniujmy macierz Jacobiego dla tego ograniczenia i liniową kombinację macierzy Hessego z dowolnym wektorem 
:
Teraz nieliniowe ograniczenie możemy zdefiniować 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 określać również w formie 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 obiekt LinearOperator:
z 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)Kiedy obliczanie macierzy Hessiana 
wymaga dużych zasobów, można wykorzystać klasę . Dostępne są następujące strategie: BFGS i SR1.
z scipy.optimize import BFGS
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())Hessian 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 Jacobiego dla ograniczeń można również obliczyć za pomocą różnic skończonych. W tym przypadku jednak macierzy Hessiana nie można obliczyć za pomocą różnic skończonych. Hessian musi być zdefiniowany w postaci funkcji lub za pomocą klasy HessianUpdateStrategy.
nonlinear_constraint = NonlinearConstraint (cons_f, -np.inf, 1, jac = '2-point', hess = BFGS ())Rozwiązanie problemu optymalizacji wygląda następująco:
z scipy.optimize import minimize
z 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` warunek zakończenia jest spełniony.
Liczba iteracji: 12, ewaluacje funkcji: 8, iteracje CG: 7, optimalność: 2.99e-09, naruszenie ograniczenia: 1.11e-16, czas wykonania: 0.033 s.
[0.41494531 0.17010937]W razie potrzeby funkcję obliczania Hessiana 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 Hessiana i dowolnego wektora przez 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, pierwsze i drugie pochodne funkcji optymalizowanej mogą być obliczone w sposób przybliżony. Na przykład, Hessian może być przybliżony za pomocą funkcji SR1 (przybliżenia quasi-Newtona). Gradient może być przybliżony za pomocą różnic skończonych.
z 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)Optymalizacja warunkowa method=»SLSQP»
Metoda SLSQP jest przeznaczona do rozwiązywania problemów minimizacji funkcji w postaci:
Gdzie 
i 
— zbiory indeksów wyrażeń opisujących ograniczenia w postaci równości lub nierówności. 
— zbiory dolnych i górnych granic dla dziedziny określenia funkcji.
Ograniczenia liniowe i nieliniowe są opisywane w postaci słowników z kluczami type, fun 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])
}Szukając minimum, postępujemy 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)Optymalizacja zakończona sukcesem. (Tryb wyjścia 0)
Aktualna wartość funkcji: 0.34271757499419825
Iteracje: 4
Oceny funkcji: 5
Oceny gradientu: 4
[0.41494475 0.1701105 ]Przykład optymalizacji
W związku z przejściem do piątego układu technologicznego, przyjrzymy się optymalizacji produkcji na przykładzie studia internetowego, które przynosi nam niewielki, ale stabilny dochód. Wyobraźmy sobie, że jesteśmy dyrektorem galerii, w której produkowane są trzy rodzaje produktów:
- x0 — lądowania sprzedażowe, od 10 tys. zł.
- x1 — strony korporacyjne, od 20 tys. zł.
- x2 — sklepy internetowe, od 30 tys. zł.
Nasz zgrany zespół roboczy składa się z czterech juniorów, dwóch midów i jednego seniora. Fundusz ich czasu pracy na miesiąc:
- juniorzy:
4 * 150 = 600 godz. * człowiek, - midzi:
2 * 150 = 300 godz. * człowiek, - senior:
150 godz. * człowiek.
Niech na opracowanie i wdrożenie jednej strony typu (x0, x1, x2) pierwszy z brzegu junior musi poświęcić (10, 20, 30) godzin, mid — (7, 15, 20), senior — (5, 10, 15) godzin najlepszego czasu swojego życia.
Jak każdy normalny dyrektor, chcemy zmaksymalizować miesięczny zysk. Pierwszym krokiem do sukcesu jest zapisanie funkcji celu value jako sumy dochodów z wyprodukowanych w ciągu miesiąca produktów:
def value(x):
return - 10*x[0] - 20*x[1] - 30*x[2]To nie jest błąd, podczas poszukiwania maksimum funkcja celu jest minimalizowana z odwrotnym znakiem.
Kolejny krok — zabrania się przetwarzania naszym pracownikom i wprowadza się ograniczenia na fundusz czasu pracy:
Co odpowiada:
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 – produkcja musi być wyłącznie pozytywna:
bnds = Bounds ([0, 0, 0], [np.inf, np.inf, np.inf])A na koniec najbarwniejsze założenie – z powodu niskiej ceny i wysokiej jakości nieustannie ustawia się kolejka z zadowolonych klientów. Możemy sami decydować o miesięcznych wolumenach produkcji, na podstawie rozwiązania zadania warunkowej optymalizacji z 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ąglimy umiarkowanie do całości i obliczymy miesięczne obciążenie wioślarzy przy optymalnym rozkładzie produkcji x = (8, 6, 3) :
- juniorzy:
8 * 10 + 6 * 20 + 3 * 30 = 290 osób * godzina; - midzi:
8 * 7 + 6 * 15 + 3 * 20 = 206 osób * godzina; - senior:
8 * 5 + 6 * 10 + 3 * 15 = 145 osób * godzina.
Wniosek: aby dyrektor otrzymywał swój zasłużony maksymalny zysk, optymalnie jest produkować miesięcznie 8 landingów, 6 średnich stron i 3 sklepy. Senior przy tym musi pracować bez przerwy, obciążenie midów wyniesie około 2/3, a juniorów mniej niż połowę.
Podsumowanie
Artykuł przedstawia podstawowe metody pracy z pakietem scipy.optimize, stosowane do rozwiązywania zadań warunkowej minimalizacji. Osobiście używam scipy wyłącznie w celach akademickich, dlatego podany przykład ma charakter żartobliwy.
Wiele teorii i przykładowych zadań można znaleźć na przykład w książce I.L. Akułicza „Programowanie matematyczne w przykładach i zadaniach”. Bardziej hardcore'owe zastosowanie scipy.optimize do budowy struktury 3D na podstawie zestawu zdjęć () można zobaczyć w .
Głównym źródłem informacji jest , chętni do współpracy w tłumaczeniu tej i innych sekcji scipy są mile widziani na .
Dziękuję za udział w przygotowaniu publikacji.
Źródło: habr.com
