Optimasi

~ Ujian Akhir Semester ~


Kontak : \(\downarrow\)
Email
Instagram https://www.instagram.com/boring.garr/
RPubs https://rpubs.com/garr/

*Library

library(lpSolve)
library(tidyquant)
library(plotly)
library(timetk)
library(DT)
library(tidyr)
library(forcats)
library(nloptr)
import numpy as np
import pandas as pd
import scipy
from scipy.optimize import linprog
from numpy import log
import plotly.graph_objects as go
import matplotlib.pyplot as plt
import pandas_datareader as web
import yfinance as yf
import datetime
from datetime import date
import plotly.express as px
import random

Buatlah contoh kasus Optimisasi Program Linear dalam Industri atau bisnis!

contoh Kasus

Sebuah bisnis kecil menjual dua produk, yaitu Produk 1 dan Produk 2. Setiap ton Produk 1 menghabiskan 30 jam kerja, dan setiap ton Produk 2 menghabiskan 20 jam kerja. Bisnis ini memiliki maksimum 2700 jam kerja untuk periode yang dipertimbangkan. Sedangkan untuk jam kerja mesin, setiap ton Produk 1 dan 2 masing-masing menghabiskan 5 dan 10 jam kerja mesin. Ada 850 jam mesin yang tersedia.

Setiap ton Produk 1 menghasilkan laba 20 dollar , sedangkan Produk 2 menghasilkan 60 dollar untuk setiap ton yang terjual. Untuk alasan teknis, perusahaan harus memproduksi minimal 95 ton secara total di antara kedua produk. Kita perlu mengetahui berapa ton Produk 1 dan 2 yang harus diproduksi untuk memaksimalkan total laba.

Variabel

maka diketahui untuk fungsi kendala: 1. x jumlah ton yang diproduksi dan dijual dari Produk 1 2. y jumlah ton yang diproduksi dan dijual dari produk 2

fungsi tujuan

Koefisien biaya dari variabel-variabel ini masing-masing adalah 20 dan 60. Oleh karena itu, fungsi tujuan didefinisikan dengan mengalikan setiap variabel dengan koefisien biaya yang sesuai.

\[\begin{align*} maxZ=20x+60y \end{align*}\]

fungsi kendala

Kendala Jam Kerja membuat bahwa jumlah total jam kerja yang digunakan pada Produk 1 dan Produk 2, yang sama dengan 30x+20y, adalah kurang atau sama dengan 2.700 jam.

Kendala yang sama, membuat jam kerja bahwa jumlah total jam kerja mesin 5x+10y kurang atau sama dengan 850.

Kendala yang membuat bahwa total unit yang diproduksi dan dijual x+y lebih besar atau sama dengan 95

perhitungan

Menyatukan semua ini, dan mempertimbangkan bahwa variabel-variabel keputusan adalah non negatif, LP yang memaksimalkan keuntungan adalah:

R

\[\begin{align*} x & \geq 0\\ y & \geq 0\\ 30x+20y & \leq 2700\\ 5x+10y & \leq 850\\ x+y & \geq 95 \end{align*}\]

tuj = c(20,60)
var = matrix(c(30, 20, 5, 10, 1, 1), nrow = 3, byrow = TRUE)
arah = c("<=", "<=", ">=")
bat = c(2700, 850, 95)

python

\[\begin{align*} x & \geq 0\\ y & \geq 0\\ 30x+20y & \leq 2700\\ 5x+10y & \leq 850\\ x+y & \geq 95 \end{align*}\]

z = [-20, -60]

lhs_ineq = [[30, 20],
[5, 10],
[-1,-1]]

rhs_ineq = [2700,
850,
-95]
bnd = [(0, float("inf")), 
(0, float("inf"))]
z = linprog(c=z, 
A_ub=lhs_ineq,
b_ub=rhs_ineq,
bounds=bnd,
method="highs")
print(z)
##            con: array([], dtype=float64)
##  crossover_nit: 0
##          eqlin:  marginals: array([], dtype=float64)
##   residual: array([], dtype=float64)
##            fun: -4900.0
##        ineqlin:  marginals: array([ -0.,  -8., -20.])
##   residual: array([600.,   0.,   0.])
##          lower:  marginals: array([0., 0.])
##   residual: array([20., 75.])
##        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
##            nit: 3
##          slack: array([600.,   0.,   0.])
##         status: 0
##        success: True
##          upper:  marginals: array([0., 0.])
##   residual: array([inf, inf])
##              x: array([20., 75.])

hasil

R

lp("max", tuj, var, arah, bat)
## Success: the objective function is 4900
lp("max", tuj, var, arah, bat)$solution
## [1] 20 75

maka nilai maksimum keuntungan yang dapat diperoleh sebesar 4900 dollar

dengan:

\[\begin{align*} x = 20\\ y = 75 \end{align*}\]

python

z.success
## True
z.fun
## -4900.0
z.x
## array([20., 75.])

maka nilai maksimum keuntungan yang dapat diperoleh sebesar 4900 dollar

dengan:

\[\begin{align*} x = 20\\ y = 75 \end{align*}\]

Buatlah contoh kasus Optimisasi Program NonLinear dalam Industri atau bisnis!

R

f <- function(x, y) {
        x^2 + y^2 + 3* x * y
}
n <- 30
xpts <- seq(-1.5, 1.5, len = n)
ypts <- seq(-1.5, 1.5, len = n)
gr <- expand.grid(x = xpts, y = ypts)
feval <- with(gr, matrix(f(x, y), nrow = n, ncol = n))
par(mar = c(5, 4, 1, 1))
contour(xpts, ypts, feval, nlevels = 20, xlab = "x", ylab = "y")
points(-1, -1, pch = 19, cex = 2)
abline(h = -1)

Python

objective_function = lambda x: (4*x[0]**2) - (4*x[0]**4) + (x[0]**(6/3) + (x[0]*x[1]) - (4*x[1]**2) + (4*x[1]**4))

X_1 = np.linspace(0,1, 100)
X_2 = np.linspace(0,1, 100)
X_1, X_2 = np.meshgrid(X_1, X_2)
Z = objective_function((X_1, X_2))

#Contour plot
fig, ax = plt.subplots()
fig.set_size_inches(14.7, 8.27)

cs = ax.contour(X_1, X_2, Z, 50, cmap='jet')
plt.clabel(cs, inline=1, fontsize=10) # plot objective function

plt.axvline(0.1, color='g', label=r'$x_1 \geq 0.1$') # constraint 1
plt.axhline(0.25, color='r', label=r'$x_2 \geq 0.25$') # constraint 2
plt.legend(bbox_to_anchor=(1, 1), loc=2, borderaxespad=0.1)

# tidy up
plt.xlabel(r'$x_1$', fontsize=16)
plt.ylabel(r'$x_2$', fontsize=16)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.show()

def geneate_starting_points(number_of_points):
    '''
    number_of_points: how many points we want to generate
    
    returns a list of starting points
    '''
    starting_points = []
    
    for point in range(number_of_points):
        
        starting_points.append((random.random(), random.random()))
        
    return starting_points
  
objective_function = lambda x: (4*x[0]**2) - (4*x[0]**4) + (x[0]**(6/3) + (x[0]*x[1]) - (4*x[1]**2) + (4*x[1]**4))

cons = [
    {'type': 'ineq', 'fun': lambda x: x[0] - 0.1}, # indoor seating >= 0.1
    {'type': 'ineq', 'fun': lambda x: x[1] - 0.25} # outdoor seating >= 0.25
]

boundaries = [(0,1), (0,1)]

from scipy.optimize import minimize

starting_points = geneate_starting_points(50)

first_iteration = True
for point in starting_points:
    res = minimize(
        objective_function,
        [point[0], point[1]],
        method='SLSQP',
        bounds=boundaries,
        constraints=cons
    )
    if first_iteration:
        better_solution_found = False
        best = res
    else:
        if res.success and res.fun < best.fun:
            better_solution_found = True
            best = res
            
if best.success:
    print(f"""Optimal solution found:
      -  Proportion of indoor seating to make available: {round(best.x[0], 3)}
      -  Proportion of outdoor seating to make available: {round(best.x[1], 3)}
      -  Risk index score: {round(best.fun, 3)}""")
else: 
    print("No solution found to problem")
## Optimal solution found:
##       -  Proportion of indoor seating to make available: 0.1
##       -  Proportion of outdoor seating to make available: 0.701
##       -  Risk index score: -0.88

Buatlah contoh kasus Optimisasi Portofolio dengan memilih 5 saham terbaik di Indonesia dengan cara melakukan analisis data keuangan terlebih dahulu pada saham LQ45!

data

R

tick <- c('BBCA.JK', 'GGRM.JK', 'JSMR.JK', 'INDF.JK', 'BBRI.JK')
price_data <- tq_get(tick,
                     from = '2020-01-01',
                     to = Sys.Date(),
                     get = 'stock.prices')
head(price_data)

Python

tick1 = ['BBCA.JK','GGRM.JK','JSMR.JK','INDF.JK','BBRI.JK']
price_data = yf.download(tick1, start='2020-01-01',
                               end= datetime.date.today())['Adj Close']
## 
[                       0%                       ]
[*******************   40%                       ]  2 of 5 completed
[**********************60%****                   ]  3 of 5 completed
[**********************80%*************          ]  4 of 5 completed
[*********************100%***********************]  5 of 5 completed
price_data.head
## <bound method NDFrame.head of                 BBCA.JK      BBRI.JK       GGRM.JK      INDF.JK      JSMR.JK
## Date                                                                        
## 2020-01-02  6322.238281  4009.359863  49895.953125  7022.275879  5155.947754
## 2020-01-03  6426.190918  4018.451172  50597.398438  7066.302734  5230.671387
## 2020-01-06  6364.764648  3972.994141  50854.593750  7044.289551  5155.947754
## 2020-01-07  6369.489746  4000.268555  52374.386719  7220.396484  5056.315918
## 2020-01-08  6312.788574  3982.085449  53683.750000  7242.409668  5056.315918
## ...                 ...          ...           ...          ...          ...
## 2022-12-16  8600.000000  4980.000000  18800.000000  6950.000000  2960.000000
## 2022-12-19  8650.000000  4970.000000  19050.000000  7000.000000  2990.000000
## 2022-12-20  8575.000000  4910.000000  18700.000000  6975.000000  2970.000000
## 2022-12-21  8675.000000  4890.000000  18250.000000  6900.000000  2960.000000
## 2022-12-22  8575.000000  4960.000000  18575.000000  6800.000000  2960.000000
## 
## [729 rows x 5 columns]>

pengembalian

log_ret_tidy <- price_data %>%
  group_by(symbol) %>%
  tq_transmute(select = adjusted,
               mutate_fun = periodReturn,
               period = 'daily',
               col_rename = 'ret',
               type = 'log')

datatable(log_ret_tidy)
log_ret_xts <- log_ret_tidy %>%
  spread(symbol, value = ret) %>%
  tk_xts()

datatable(log_ret_xts)

Rata-rata pengembalian

R

mean_ret <- colMeans(log_ret_xts)
print(round(mean_ret, 4))
## BBCA.JK BBRI.JK GGRM.JK INDF.JK JSMR.JK 
##  0.0004  0.0003 -0.0014  0.0000 -0.0008

Python

log_ret = np.log(price_data/price_data.shift(1))
log_ret
##              BBCA.JK   BBRI.JK   GGRM.JK   INDF.JK   JSMR.JK
## Date                                                        
## 2020-01-02       NaN       NaN       NaN       NaN       NaN
## 2020-01-03  0.016309  0.002265  0.013960  0.006250  0.014389
## 2020-01-06 -0.009605 -0.011377  0.005070 -0.003120 -0.014389
## 2020-01-07  0.000742  0.006841  0.029447  0.024693 -0.019513
## 2020-01-08 -0.008942 -0.004556  0.024693  0.003044  0.000000
## ...              ...       ...       ...       ...       ...
## 2022-12-16  0.011696  0.014156 -0.018446 -0.003591 -0.026668
## 2022-12-19  0.005797 -0.002010  0.013210  0.007168  0.010084
## 2022-12-20 -0.008708 -0.012146 -0.018544 -0.003578 -0.006711
## 2022-12-21  0.011594 -0.004082 -0.024358 -0.010811 -0.003373
## 2022-12-22 -0.011594  0.014213  0.017652 -0.014599  0.000000
## 
## [729 rows x 5 columns]

Matriks Kovariansi

R

cov_mat <- cov(log_ret_xts) * 252
print(cov_mat)
##            BBCA.JK    BBRI.JK    GGRM.JK    INDF.JK    JSMR.JK
## BBCA.JK 0.08235240 0.06324033 0.03482309 0.03257385 0.04708422
## BBRI.JK 0.06324033 0.14450743 0.04491092 0.04588698 0.06113619
## GGRM.JK 0.03482309 0.04491092 0.13178566 0.04787952 0.04365387
## INDF.JK 0.03257385 0.04588698 0.04787952 0.09740420 0.03875918
## JSMR.JK 0.04708422 0.06113619 0.04365387 0.03875918 0.16183444

Python

cov_mat = log_ret.cov() * 252

print (cov_mat)
##           BBCA.JK   BBRI.JK   GGRM.JK   INDF.JK   JSMR.JK
## BBCA.JK  0.082466  0.063327  0.034871  0.032619  0.047149
## BBRI.JK  0.063327  0.144706  0.044973  0.045950  0.061220
## GGRM.JK  0.034871  0.044973  0.131966  0.047945  0.043714
## INDF.JK  0.032619  0.045950  0.047945  0.097538  0.038812
## JSMR.JK  0.047149  0.061220  0.043714  0.038812  0.162057

Penerapan Metode Portofolio

R

wts = runif(n = length(tick))
wts = wts/sum(wts)

port_returns = (sum(wts * mean_ret) + 1)^252 - 1

port_risk = sqrt(t(wts) %*% (cov_mat %*% wts))

sharpe_ratio = port_returns/port_risk

num_port = 5000

all_wts = matrix(nrow = num_port, ncol = length(tick))

port_returns = vector('numeric', length = num_port)
 
port_risk = vector('numeric', length = num_port)

sharpe_ratio = vector('numeric', length = num_port)
for (i in seq_along(port_returns)) {
  wts <- runif(length(tick))
  wts <- wts/sum(wts)
  
  all_wts[i,] <- wts
  
  port_ret <- sum(wts * mean_ret)
  port_ret <- ((port_ret + 1)^252) - 1
  
  port_returns[i] <- port_ret
  
  port_sd <- sqrt(t(wts) %*% (cov_mat %*% wts))
  port_risk[i] <- port_sd
  
  sr <- port_ret/port_sd
  sharpe_ratio[i]<-sr
}
portfolio_values <- tibble(Return = port_returns,
                           Risk = port_risk,
                           SharpeRatio = sharpe_ratio)

all_wts <- tk_tbl(all_wts)

colnames(all_wts) <- colnames(log_ret_xts)

portfolio_values <- tk_tbl(cbind(all_wts, portfolio_values))

datatable(portfolio_values)

Python

num_port = 5000
all_wts = np.zeros((num_port, len(price_data.columns)))
port_returns = np.zeros((num_port))
port_risk = np.zeros((num_port))
sharpe_ratio = np.zeros((num_port))
for i in range(num_port):
 wts = np.random.uniform(size = len(price_data.columns))
 wts = wts/np.sum(wts)

 all_wts[i,:] = wts

 port_ret = np.sum(log_ret.mean() * wts) 
 port_ret = (port_ret + 1) ** 252 - 1

 port_returns[i] = port_ret

 port_sd = np.sqrt(np.dot(wts.T, np.dot (cov_mat, wts)))

 port_risk[i] = port_sd

 sr = port_ret/port_sd
 sharpe_ratio[i] = sr
names = price_data.columns
min_var = all_wts[port_risk.argmin()]
print(min_var)
## [0.42253581 0.04006174 0.17059677 0.28802686 0.07877882]
max_sr = all_wts[sharpe_ratio.argmax()]
print(max_sr)
## [0.80095553 0.03014462 0.00889255 0.15770316 0.00230415]
print(port_risk.min())
## 0.23633168168296384
print(sharpe_ratio.max())
## 0.32531171482990856

Variansi Minimum

R

min_var <- portfolio_values[which.min(portfolio_values$Risk),]

p <- min_var %>%
  gather(BBCA.JK:JSMR.JK, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Aset',
       y = 'Bobot',
       title = "Bobot Portofolio dengan Variansi Minimum") +
  scale_y_continuous(labels = scales::percent) +
  theme(legend.position="none")

ggplotly(p)

Python

min_var = pd.Series(min_var, index=names)
min_var = min_var.sort_values()

fig = px.bar(min_var, color=min_var.index, title="Portofolio dengan variansi umum")
fig.show()

Portofolio Tangensi

R

max_sr<- portfolio_values[which.max(portfolio_values$SharpeRatio), ]

p <- max_sr %>%
  gather(BBCA.JK:JSMR.JK, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Aset',
       y = 'Bobot',
       title = "Bobot Portofolio Tangensi (Maksimum Sharpe Ratio)") +
  scale_y_continuous(labels = scales::percent)+
  theme(legend.position="none")

ggplotly(p)

Python

max_sr = pd.Series(max_sr, index=names)
max_sr = max_sr.sort_values()

fig = px.bar(max_sr, color=max_sr.index, title="Portofolio dengan tangensi")
fig.show()

Batas Efisien Portofolio

R

p <- portfolio_values %>%
  ggplot(aes(x = Risk, y = Return, color = SharpeRatio)) +
  geom_point() +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = 'Risiko Tahunan',
       y = 'Pengembalian Tahunan',
       title = "Optimasi Portofolio & Perbatasan yang Efisien") +
  geom_point(aes(x = Risk, y = Return), data = min_var, color = 'black') +
  geom_point(aes(x = Risk, y = Return), data = max_sr, color = 'red') +
  annotate('text', x = 0.25, y = 0.10, label = "Portofolio Tangensi") +
  annotate('text', x = 0.255, y = -0.05, label = "Portofolio Varians minimum")
  

ggplotly(p)

Python

fig = plt.figure()
ax1 = fig.add_axes([0.1,0.1,0.8,0.8])
ax1.set_xlabel('risk')
ax1.set_ylabel("returns")
ax1.set_title("Portofolio optimization and efficient frontier")
plt.scatter(port_risk, port_returns)
plt.show()