Kategorie
coding

Powrót do szkoły czyli optymalizacja liniowa

Tworząc oprogramowanie sporą część naszego czasu zajmują raczej mało interesujące i trywialne zadania, np. kolejny endpoint z użyciem Springa. Ta monotonia, co jakiś czas jest przeplatana czymś bardziej interesującym, czymś co pozwoli pobudzić nasze szare komórki. Jest wiele powodów, dla których programiści wybierają swój zawód. Jednym z nich, dla mnie przynajmniej, jest spora satysfakcja z rozwiązywania ciekawych problemów, a poszukiwanie optymalnego rozwiązania, jak najbardziej się do nich zalicza.

Nasz problem

Zdefiniujmy ramy problemu, który chcemy rozwiązać. Mamy szereg inwestorów, i każdy z nich chce zainwestować swój kapitał w akcje na giełdzie. Dobrym nawykiem jest nieinwestowanie wszystkich swoich pieniędzy w jedną firmę, więc zakładamy arbitralnie, że maksymalnie w jedną firmę możemy zainwestować 25% kapitału. Same akcje i firmy za nimi stojące mają ocenę jakości (losowo przypisaną) podobną do tej jaką oferuje S&P, choć patrząc na film Big Short, nasze oceny mogą być równie trafne. Co bardziej zorientowani inwestorzy mają pewne preferencje i chcą mieć ekspozycję na konkretny sektor, np. energetykę. To daje nam już pewien zarys zadania.

Krystalizuje się więc powoli nasz model, który będzie stanowić bazę dla naszego problemu:

...
// Peter has 100k $ to invest and wants to invest at least 0.3 of his capital in energy, finance and it
new Investor("peter", 100000, new Strategy("it", 0.3, "energy", 0.3, "finance", 0.3))
...
// OXY is a ticker for a company from energy sector, quality AAA and price 22.39 $
new Stock("OXY", "Occidental Petroleum Corp", AAA, "energy", 22.39)

Problem specyficzny i niszowy, podobnie jak poszukiwanie szybkości w systemach, jednak bardzo szybko można zapędzić się w kozi róg, jeśli podejdzie się do problemu naiwnie. Możemy zmierzyć się z zadaniem w brutalny sposób, ale brute-force ma tę cechę, że próbuje sprawdzić całą możliwą gamę danych, które niekoniecznie spełniają kryteria zadania. Konstrukcja takiego „prostego” podejścia często wymaga tworzenia wielu dodatkowych i małoczytelnych fragmentów, które utrzymywałyby kryteria wewnątrz zadanych przedziałów. Możliwe, że nawet uda stworzyć poprawny algorytm, ale jeśli będziemy chcięli wziąć pod uwagę kolejne czynniki czy też zwiększyć próbkę danych, rozszerzając rozmiar poszukiwanych wyników sto razy, skomplikowanie problemu może nasz algorytm zabić na miejscu. Jeśli jakimś cudem przeżyje, dowiemy się o tym z ogromnym opóźnieniem, bo przejście przez całą powiększoną przestrzeń rozwiązań zajmie odpowiednio długo. Nie mamy na to czasu…

Nowy młotek w skrzynce z narzędziami

Skoro problem nie jest zwykłym gwoździem, trzeba poszukać równie niezwykłego młotka. Patrząc pod odpowiednim kątem, dostrzerzemy, że problem, który mamy rozwiązać to jest nic innego jak zestaw równań liniowych jakimi zajmowaliśmy się w liceum (albo i wcześniej). Jedyną różnicą jest to, że zamiast trzech niewiadomych możemy mieć ich znacznie więcej. Przełożeniem takich zadań na problem programistyczny jest programowanie liniowe, które jest niczym innym jak zaprezentowaniem powyższych równań w akceptowalnej formie dla naszego języka programowania oraz oznaczenie jednego z tych równań specjalną rolą, mianowicie aby pełniło rolę funkcji kosztu. To ta funkcja będzie odpowiadać za proces optymalizacji, bo wedle niej, wszelkie inne parametry będą dopasowywane, aby maksymalizować bądź też minimalizować końcowy wynik, na którym nam zależy. Znając języki jak choćby java, nie do końca intuicyjnie jest przejście pokroju

Ax + b = 0

do

public int linear(a, x, b);

już nie mówiąc jak implementacja takiej funkcji powinna wyglądać. Musi ona zapewniać mechanizm znany jako solver, czyli oprogramowanie obliczeniowe potrafiące rozwiązywać równania. Nie jest to trywialne, ale jak z wieloma aspektami IT, ktoś już stanął przed tym problemem przed Tobą i aktualnie jest dostępnych wiele rozwiązań i bibliotek oferujących taką funkcjonalność. Każde, oferuje pewien format oraz konwencję, jak równania powinny być opisane, ale rezultat co do zasady jest podobny. Różnice oddzielające lepsze solvery od gorszych jest prędkość wykonywania obliczeń.

Jedną z istniejących już bibliotek jest ortools, stworzona przez Google i oferuje dokładnie to czego potrzebujemy, czyli liniowy solver. W takim wypadku nie pozostaje nam nic innego, jak przedstawić nasz problem w formie zestawu równań, a później wyrazić je w konwencji biblioteki.

Szybki skok z teorii do praktyki

Na początek potrzebujemy zdefiniować dla każdego inwestora zbiór danych reprezentujący ilość akcji każdej firmy jaka składa się na jego portfolio.

kapitał\ inwestora = ilość_{1}*cena_1 + ilość_{2}*cena_2 + ilość_{2}*cena_2 + \cdots +  ilość_{n}*cena_n

Jeśli powyższe zrobimy dla każdego inwerstora i uszeregujemy je kolejno pod sobą, to okaże się, że umieszczając cenę akcji jako nagłówek kolumny otrzymujemy doskonała reprezentację naszego problemu w postaci macierzy:

\begin{matrix} 
& 
\begin{matrix} c_1 \ & c_2 \  & \cdots \ & c_n \ \end{matrix}
\\ 
\begin{matrix} ki_1 \\ ki_2 \\ \vdots \\ ki_m\end{matrix} 
& 
\begin{bmatrix}
i_{1,1} & i_{1,2} & \cdots & i_{1,n} \\
i_{2,1} & i_{2,2} & \cdots & i_{2,n} \\
\vdots  & \vdots  & \ddots & \vdots  \\
i_{m,1} & i_{m,2} & \cdots & i_{m,n} \end{bmatrix} 
\end{matrix} 
\\ \ 
\\ ki \Rightarrow kapitał \ inwestora, c \Rightarrow cena \  akcji, i \Rightarrow ilość \ akcji \ danej \ firmy \\

Powyższa macierz natomiat daje nam doskonałe możliwości, nie tylko do zdefiniowania bazowego problemu, ale także pozwala na doprecyzowanie kolejnych wymagań, które stawiamy przed naszym problemem. Pierwszym, które zakładaliśmy jest umieszczenie do 25% kapitału w akcjach jednej firmy, czyli:

ilość_{m,n} \le 0.25 * \frac{kapitał_m}{cena\ akcji_n}

Mając tak opisany problem, mamy już na tyle dużo informacji, aby zacząć tworzyć rozwiązanie. Struktura programu nie będzie wyjątkowo skomplikowana i skupia się na aplikowaniu kolejnych kryteriów wobec naszych zmiennych. Dla potencjalnych inwestorów investors i dostępnych akcji stocks, pobierzemy dostępny solver z ortools, a następnie stworzymy nasze rozwiązanie Solution, które będzie zamykać w sobie wszelkie potencjalne zmienne dla pary inwestor:akcja.

    ...
    public static void main(String[] args) {
        Investor[] investors = new Investor[]{ ... };
        Stock[] stocks = new Stock[]{ ... };

        Loader.loadNativeLibraries();
        MPSolver solver = MPSolver.createSolver("GLOP");
        Solution solution = new Solution(investors, stocks);
        createAmountVariablesForEachStockPerInvestor(solver, solution);
        ...
        // next constraints and requirements
    }
    ...

Zmienne wraz z pierwszym ograniczeniem nadmiernego inwestowania w jedno aktywo zapisujemy do tablicy zmiennych, która przypomina macierz wartości opisaną nieco wyżej. Przyjmujemy, że każdy inwestor może nie posiadać akcji danej firmy, więc najmniejszą przydzieloną wartością jest 0.

   ...
   private static void createAmountVariablesForEachStockPerInvestor(MPSolver solver, Solution solution) {
        Investor[] investors = solution.getInvestors();
        Stock[] stocks = solution.getStocks();
        for (int i = 0; i < investors.length; i++) {
            for (int j = 0; j < stocks.length; j++) {
                Investor investor = investors[i];
                Stock stock = stocks[j];
                double amount = investor.getAvailableFunds() / stock.getPrice() * 0.25;
                solution.getVars()[i][j] = solver.makeIntVar(0.0, amount, investor.getName() + "-" + stock.getName());
            }
        }
    }
    ...

Następnie chcemy opisać główną cechę naszego zadania, tę, która również posłużyła za pierwowzór macierzy, czyli sumę wartości poszczególnych inwestycji w akcje dla każdego inwestora.

ki_m \ge \sum_{n=1}^N i_{m,n}*c_{m,n}
\\ \ 
\\ ki_m \Rightarrow kapitał \ inwestora_m, c \Rightarrow cena \  akcji, i \Rightarrow ilość \ akcji \ danej \ firmy \\

Jest to kolejne ograniczenie, a w bibliotece jakiej używamy wyraża się to poprzez stworzenie obiektu MPConstraint wraz ze zdefiniowanych współczynnikiem (dla nas jest to cena danej akcji), bo zmienną, której poszukujemy jest ilość akcji.

    ...
    private static void applyConstraintToLimitPortfolioToAvailableFunds(MPSolver solver, Solution solution) {
        Investor[] investors = solution.getInvestors();
        Stock[] stocks = solution.getStocks();
        for (int i = 0; i < investors.length; i++) {
            Investor investor = investors[i];
            MPConstraint constraint = solver.makeConstraint(-INFINITY, investor.getAvailableFunds(), investor.getName() + "-available-funds");
            for (int j = 0; j < stocks.length; j++) {
                Stock stock = stocks[j];
                MPVariable variable = solution.getVars()[i][j];
                constraint.setCoefficient(variable, stock.getPrice());
            }
        }
    }
    ...

Następnie, musimy wziąć pod uwagę bardziej zorientowanych inwestorów, którzy chcą mieć akcje z pewnego, wybranego przez siebie sektora. W naszej implementacji przykładowa wartość 0.3 poniżej, informuje, że taka część całego kapitału powinna trafić do tego sektora, a więc taką część wszystkich dostępnych pieniędzy przeznaczamy na akcje w danym sektorze.

new Strategy("it", 0.3, "energy", 0.3, "finance", 0.3)

Nie każdy inwestor ma takie wymaganie, ale ci, którzy mają, będą mieli dodane warunki do swojego „wiersza” zmiennych w macierzy rozwiązania dla każdego zdefiniowanego sektora którym się interesują.

    ...
    private static void applyConstraintToLimitMarketExposure(MPSolver solver, Solution solution) {
        Investor[] investors = solution.getInvestors();
        Stock[] stocks = solution.getStocks();
        for (int i = 0; i < investors.length; i++) {
            Investor investor = investors[i];
            Map<String, Double> sectorToStockShare = investor.getStrategy().getSectorToStockShare();
            for (Map.Entry<String, Double> entry : sectorToStockShare.entrySet()) {
                ...
                double limit = entry.getValue();
                double partOfFunds = investor.getAvailableFunds() * limit;
                MPConstraint constraint = solver.makeConstraint(partOfFunds, INFINITY, name);
                for (int j = 0; j < stocks.length; j++) {
                    Stock stock = stocks[j];
                    if (Objects.equals(sector, stock.getSector())) {
                        MPVariable variable = solution.getVars()[i][j];
                        constraint.setCoefficient(variable, stock.getPrice());
                    }
                }
            }
        }
    }
    ...

Mając tak przygotowane zadanie, pozostaje nam już tylko podać faktyczny cel optymalizacji, który bierze pod uwagę wszystkie obostrzenia, które dodaliśmy powyżej. W ortools służy do tego MPObjective, który definiuje się bardzo podobnie jak ograniczenia, ale nie podaje wartości ograniczającej sumę składowych. Określa się jedynie poszukiwane ekstremum i definiuje czy jest to maksymalna czy minimalna wartość. Reszta należy już do biblioteki.

max\left( \  
\sum_{n=1}^{N_{AAA}} AAA_{współ}*i_{n}*c_{n} \  + 
\sum_{m=1}^{M_{AA}} AA_{współ}*i_{m}*c_{m} \  + 
\cdots +
\sum_{o=1}^{O_{B}} B_{współ}*i_{o}*c_{o} \ 
\right)

\\ \ 
\\ AAA_{współ},...,B_{współ} \Rightarrow przypisana \ waga , c \Rightarrow cena \  akcji , i \Rightarrow ilość \ akcji \ danej \ jakości \\

W naszym przypadku zależy nam, aby wybierane akcje do portfolio każdego inwestora były możliwie najwyższej jakości, biorąc pod uwagę wszystkie dostępne aktywa i jednocześnie wymagania jakie zdefiniowaliśmy. Pan Belfort pewnie z nami się nie zgodzi, ale akcjom o najwyższej jakości AAA, przyznamy pełną wartość zakupu jako punkty. Im gorszej jakości akcje, tym mniej tych punków zostanie przyznanych i w ten sposób zdefiniujemy naszą preferencję.

    ...
    private static void applyObjectiveToMaximizeUsageOfHighestQualityAsset(MPSolver solver, Solution solution) {
        Investor[] investors = solution.getInvestors();
        Stock[] stocks = solution.getStocks();
        MPObjective objective = solver.objective();
        for (int i = 0; i < investors.length; i++) {
            for (int j = 0; j < stocks.length; j++) {
                Stock stock = stocks[j];
                MPVariable variable = solution.getVars()[i][j];
                double gainValue = stock.getRating().getFactor() * stock.getPrice();
                objective.setCoefficient(variable, gainValue);
            }
        }
        objective.setMaximization();
    }
    ...

Uruchomienie rozwiązania,

    ...
    MPSolver.ResultStatus result = solver.solve();
    System.out.println("The solver came with " + result.name() + " solution. It took "
            + solver.iterations() + " iterations within "
            + solver.wallTime() + " milliseconds."
    );
    solution.summary();
    ...

spowoduje wypisanie rezultatu i zarazem sugestii dla naszych inwestorów, co zrobić z dostępnym kapitałem.

The solver came with OPTIMAL solution. It took 22 iterations within 33 milliseconds.
            | 'OXY',r:AAA|  'CXO',r:AA|   'APA',r:B|  'COP',r:AA|  'NOV',r:BB|  'AFL',r:AA| 'AIG',r:AAA|  'AIZ',r:AA| 'AJG',r:AAA|   'ALL',r:A| 'ACN',r:AAA|'ADBE',r:AAA|  'AMD',r:BB|'AKAM',r:AAA| 'APH',r:BBB
            |     p:22.39|     p:65.60|     p:17.27|     p:45.12|     p:14.75|     p:46.48|     p:41.35|    p:140.17|    p:116.11|    p:108.40|    p:253.65|    p:458.08|     p:88.21|    p:106.45|    p:131.74
            |      energy|      energy|      energy|      energy|      energy|     finance|     finance|     finance|     finance|     finance|          it|          it|          it|          it|          it
------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------
       peter|     1116.57|       76.22|        0.00|        0.00|        0.00|        0.00|      120.92|        0.00|      215.31|        0.00|        0.00|       32.75|        0.00|      234.85|        0.00
         bob|        0.00|        0.00|        0.00|        0.00|        0.00|        0.00|      362.76|        0.00|      215.31|        0.00|       39.42|       54.58|        0.00|      234.85|        0.00
        john|     1116.57|        0.00|        0.00|      554.08|        0.00|        0.00|      604.59|        0.00|      215.31|        0.00|        0.00|        0.00|        0.00|        0.00|        0.00
        fred|        0.00|        0.00|        0.00|        0.00|        0.00|        0.00|        0.00|        0.00|      215.31|        0.00|       98.56|       54.58|        0.00|      234.85|        0.00

Wygląda to nieźle, bo nasze wymagania pozwalały na znalezienie optymalnego rozwiązania i zgodnie z oczekiwaniami, wszystkie kryteria są zachowane. Jest to istotne, bo wystarczy, że nieco „przegniemy” ze strategiami danego inwestora i problem może okazać się nierozwiązywalny. Podobnie z wyborem akcji danej firmy, ponieważ ograniczamy ich użycie tylko na poziomie jednego inwestora, a nie wszystkich inwestorów. Możemy zaobserwować tą sytuację w podobnych przydziałach akcji w niektórych miejscach (peter i john mają tyle samo akcji OXY). Jedynym mankamentem jaki możemy zauważyć jest pewna trudność z ewentualnym zakupem 1116.57 akcji i wypadałoby to naprawić. Winowajcą jest tutaj użyty na początku solver GLOP, działający na liczbach rzeczywistych, a wyniki jakich oczekujemy lepiej, aby były wyrażone w liczbach całkowitych. Użycie solvera SCIP (czyli zaledwie zamiana wartości w pierwszych liniach kodu) sprawia, że otrzymujemy dokładnie to czego potrzebujemy. Jest to okupione niestety dłuższymi obliczeniami, ale rozwiązanie dla nas jest dużo bardziej użyteczne. Różnice natomiast pomiędzy tymi dwoma solverami, i co za tym idzie czasem obliczeń, najlepiej wyjaśni ten opis.

The solver came with OPTIMAL solution. It took 5289233 iterations within 540606 milliseconds.
            | 'OXY',r:AAA|  'CXO',r:AA|   'APA',r:B|  'COP',r:AA|  'NOV',r:BB|  'AFL',r:AA| 'AIG',r:AAA|  'AIZ',r:AA| 'AJG',r:AAA|   'ALL',r:A| 'ACN',r:AAA|'ADBE',r:AAA|  'AMD',r:BB|'AKAM',r:AAA| 'APH',r:BBB
            |     p:22.39|     p:65.60|     p:17.27|     p:45.12|     p:14.75|     p:46.48|     p:41.35|    p:140.17|    p:116.11|    p:108.40|    p:253.65|    p:458.08|     p:88.21|    p:106.45|    p:131.74
            |      energy|      energy|      energy|      energy|      energy|     finance|     finance|     finance|     finance|     finance|          it|          it|          it|          it|          it
------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------|------------
       peter|     1116.00|       73.00|        0.00|        5.00|        0.00|        0.00|      604.00|        0.00|      126.00|        0.00|       98.00|        0.00|        0.00|       52.00|        0.00
         bob|      264.00|        0.00|        0.00|        0.00|        0.00|        0.00|      361.00|        0.00|      165.00|        0.00|       98.00|       50.00|        0.00|      115.00|        0.00
        john|     1116.00|      264.00|        4.00|      169.00|        0.00|        2.00|      600.00|        2.00|      210.00|        4.00|        0.00|        0.00|        0.00|        0.00|        0.00
        fred|      156.00|        0.00|        0.00|        0.00|        0.00|        0.00|      604.00|        0.00|      215.00|        0.00|       98.00|       46.00|        0.00|        6.00|        0.00
Słowo na koniec

Jak widać (mam nadzieję, że faktycznie jest to widoczne), pomimo mało urzekającego pierwszego wrażenia, takie podejście do problemu pokazuje spory potencjał przy procesie optymalizacji. Trochę bardziej wyrafinowany sposób niż zwykłe przeszukanie całej dostępnej przestrzeni rozwiązań, wiążące się z pewnym dodatkowym nakładem pracy, ale w momencie wzrostu skomplikowania przestrzeni rozwiązań czy też znaczącego wzrostu możliwych wyników szybko pokaże swoją przewagę. Od razu też zaznaczę, iż kod, który stworzyłem, ma za zadanie pokazać użycie biblioteki, a nie podążać za najnowszymi wytycznymi clean code, więc purystów proszę o odrobinę wyrozumiałości, tym razem 😉 .

Swoją drogą, ciekawe jak wiele osób na codzień mierzy się z tego typu problemami, czy jest to tylko mój wymysł czy jednak narzędzia podobne do ortools faktycznie mają miejsce w naszym ulubionym miejscu… na produkcji. Jeśli masz okazję rozwiązywać tego typu problemy, z chęcią usłyszę jak do tego podchodzisz w swoich projektach.

Cały projekt, który posłużył jako baza do tego wpisu możesz znaleźć tutaj: https://bitbucket.org/PawelRebisz/ortools-usage/src/master/ i uruchomić własnoręcznie. Aby w pełni z tego skorzystać będzie potrzebna instalacja biblioteki ortools, która jest doskonale opisana tutaj.

Dodaj komentarz

Twój adres e-mail nie zostanie opublikowany. Wymagane pola są oznaczone *