Email:
RPubs: https://rpubs.com/akun-rpubs-anda/


Packages

library(dplyr)      #data manipulation
library(ggplot2)    #data visualization
library(DT)         #beautify the data table
library(plotly)     #make a pie chart
library(Metrics)    #find error value
library(tidyverse)
library(yardstick)
library(data.table)
library(paletti)
library(GGally)
library(plotly)
library(rsample)      # Initial Split
library(partykit) 
library(rpart)        # Decision Tree
library(rpart.plot)   # Decision Tree
library(caret)        # Confussion Matrix
library(randomForest) # random forest

library(devtools)
#devtools::install_github('skinner927/reprtree')
#library(reprtree)

library(e1071)        # naive bayes
library(nnet)         # multinomial logistic regression
library(gridExtra)
library(grid)
library(knitr)
library(kableExtra)
library(cowplot)
library(formattable)

#COLORS
library(ggthemes)
library(paletti)
# WARNA
mycolorfill = c(
  
  light_blue="#2f4b7c", 
  smooth_blue ="#4B87CB",
  light_purple ="#665191",
  dark_pink="#a05195", 
  light_pink="#d45087", 
  light_red="#f95d6a", 
  dark_orange="#ff6347",
  semi_orange="#e79658",
  orange="#dda15a",
  cream="#b59378",
  dark_cream="#A57F5F",
  choc="#85664B",
  dark_choc="#6b5340",
  light_orange="#ff7c43"
)

#viz_palette(mycolorfill)
mycolor_fill  <- get_scale_fill(get_pal(mycolorfill))
mycolor_color <- get_scale_color(get_pal(mycolorfill))
mycolor_hex <- get_hex(mycolorfill)

1 Soal 1

Please figure out how to optimize the cost effectively to blend your sausages.

Ingredient      <- c("Chicken","Wheat","Starch")
Cost            <- c(4.32,2.46,1.86)
Availability    <- c(30,20,17)
Sausage         <- data.frame(Ingredient,Cost,Availability)
Sausage
##   Ingredient Cost Availability
## 1    Chicken 4.32           30
## 2      Wheat 2.46           20
## 3     Starch 1.86           17

So we need to make the following variables above in mathematical way,

\(es\) = Economy Sausage

\(ps\) = Premium Sausage

\(c\) = Chicken

\(w\) = Wheat

\(s\) = Starch

The function if we want to minimize the cost is:

\[MIN(Cost)=4.32(c_{es}+c_{ps})+2.46(w_{es}+w_{ps})+1.86(s_{es}+s_{ps})\]

And the constraints is:

~Wheat = \(w_{es}+w_{ps}≤20\)

~Chicken = $c_{es} + c{ps} ≤ 53 $

~Starch = \(s_{es}+s_{ps}≤17\)

~Ingredient of Chicken in Economy Sausage = \(c_{es}≥0.4(c_{es}+w_{es}+s_{es})\)

~Ingredient of Chicken in Premium Sausage = \(c_{ps}≥0.6(c_{ps}+w_{ps}+s_{ps})\)

~Ingredient of Starch in Economy Sausage = \(s_{es}≥0.25(c_{es}+w_{es}+s_{es})\)

~Ingredient of Starch in Economy Sausage = \(s_{ps}≥0.25(c_{ps}+w_{ps}+s_{ps})\)

~The Profit of Economy Sausage = \(0.05×350\)

~The Profit of Economy Sausage = \(0.05×500\)

2 Soal 2

Please find out the optimization problem of using Multinomial Linear Regression and Nonlinier Regresion at Rpubs or Kaggle (translate and submit it to your Rpubs!).

2.1 Multinomial Regression

2.1.1 Pengenalan

Ini adalah latihan regresi linier menggunakan dataset Fish Market. Dataset ini berisi beberapa pengukuran untuk tujuh spesies ikan yang berbeda. Dalam analisis ini, kita akan melakukan model regresi linier berganda dan juga memeriksa bagaimana data yang berbeda dapat mempengaruhi model regresi, misalnya dengan menghitung nilai jarak juru masak dan \(R^2\) dalam skenario yang berbeda.

Hal pertama yang pertama. Saya akan dengan mudah memuat sebagian besar paket yang akan saya gunakan dan membaca dataset:

suppressPackageStartupMessages(library(rmarkdown))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(broom))
suppressPackageStartupMessages(library(car))
suppressPackageStartupMessages(library(ggExtra))
suppressPackageStartupMessages(library(plotly))

fish <- read_csv("Fish.csv")
## 
## -- Column specification --------------------------------------------------------
## cols(
##   Species = col_character(),
##   Weight = col_double(),
##   Length1 = col_double(),
##   Length2 = col_double(),
##   Length3 = col_double(),
##   Height = col_double(),
##   Width = col_double()
## )

Semua variabel terlihat dalam format yang benar. Tidak perlu mempermasalahkan data ini, karena sudah rapi.

Sekarang akan memeriksa bingkai data:

head(fish)
## # A tibble: 6 x 7
##   Species Weight Length1 Length2 Length3 Height Width
##   <chr>    <dbl>   <dbl>   <dbl>   <dbl>  <dbl> <dbl>
## 1 Bream      242    23.2    25.4    30     11.5  4.02
## 2 Bream      290    24      26.3    31.2   12.5  4.31
## 3 Bream      340    23.9    26.5    31.1   12.4  4.70
## 4 Bream      363    26.3    29      33.5   12.7  4.46
## 5 Bream      430    26.5    29      34     12.4  5.13
## 6 Bream      450    26.8    29.7    34.7   13.6  4.93

Nama kolom tidak terlalu intuitif, jadi di bawah ini saya menjelaskan masing-masing (yang juga dijelaskan di halaman dataset asli):

Species: spesies yang diukur, total ada tujuh.

Weight: Berat ikan dalam gram.

Length1: panjang vertikal dalam cm.

Length2: panjang diagonal dalam cm.

Length3: panjang silang dalam cm.

Height: Tinggi dalam cm.

Width: lebar diagonal dalam cm.

2.1.2 Menguji variabel penjelas dan memeriksa kebisingan

Kita akan menggunakan berat sebagai variabel prediktor kita. Untuk memulai analisis kita, kita akan memplot korelasi antara variabel prediktor kita dan semua variabel numerik yang tersisa dan menghitung \(R^2\) untuk setiap kasus.

Untuk memplot semua informasi ini, kita perlu sedikit mengubah struktur bingkai data:

options(repr.plot.width = 20, repr.plot.height = 10)

fish %>% # changing the data frame structure
  pivot_longer(cols=names(fish)[3:length(names(fish))],names_to = "explanatory") %>% 
  #lotting
  ggplot(aes(y=Weight,x=value)) + geom_point() + geom_smooth(method="lm",linetype=2) +
  theme_bw() +
  theme(text = element_text(size = 25)) +
  ggtitle("Weight vs. height, length measurements and width") +
  ylab("Weight (g) - predictor variable") +
  xlab("Explanatory variable (free scale)") +
  facet_wrap(~explanatory, scales="free")
## `geom_smooth()` using formula 'y ~ x'

Dari slope, dapat dilihat bahwa korelasi antara bobot dan variabel yang tersisa adalah positif. Plot ini, bukanlah cara yangkita akan menghitung nilai-nilai ini dalam kerangka data, dan mengurutkannya dari nilai yang lebih tinggi. Caranya adalah:

rsquared_mod <- function(x) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = fish)
    )[["r.squared"]]
  )
}

# explanatory variables for which I will apply the function
explanatory <- names(fish)[3:length(names(fish))]

# data frame with r squared values in decreasing order

tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod))) %>% 
    arrange(desc(r_squared))
## # A tibble: 5 x 2
##   variables r_squared
##   <chr>         <dbl>
## 1 Length3       0.852
## 2 Length2       0.844
## 3 Length1       0.839
## 4 Width         0.786
## 5 Height        0.525

Nilai-nilai di atas mengacu pada \(R^2\) untuk nilai korelasi di setiap plot, yang pada dasarnya merupakan bobot nilai terhadap semua variabel lainnya, secara terpisah. Tiga ukuran panjang adalah ukuran yang paling tepat menggambarkan berat ikan. Tingginya kurang dapat diandalkan dan lebarnya berdiri di antaranya. Nilai \(R^2\) yang tinggi untuk ukuran panjang, diatasi dengan fakta bahwa ini pada dasarnya merupakan hal yang serupa, dan menunjukkan kemungkinan kolinearitas di antara kedua variabel tersebut. Ini berarti bahwa, dalam regresi multivariat, variabel-variabel ini menjadi mubazir. Kita akan memeriksanya nanti.

Untuk saat ini, saya akan melanjutkan dengan regresi linier sederhana, menggunakan variabel Length3 sebagai penjelas, karena ia memiliki nilai \(R^2\) highwest.

Di bawah ini, koefisien model:

lm(Weight~Length3,data=fish)
## 
## Call:
## lm(formula = Weight ~ Length3, data = fish)
## 
## Coefficients:
## (Intercept)      Length3  
##     -490.40        28.46

Menurut model ini, ada peningkatan sekitar 28 gram pada berat ikan tertentu seiring bertambahnya panjang silang dalam satuan (yaitu 1 cm).

Untuk memeriksanya secara visual, kita akan memplot regresi, lalu memetakan warna titik-titik ke spesies ikan yang dimaksud. Caranya adalah:

library(paletteer) # r package with color palettes

fish %>% ggplot(aes(y=Weight,x=Length3)) + 
  geom_point(alpha=0.7,aes(color=Species),size=10) + 
  geom_smooth(method="lm", linetype=2) +
  theme_bw() +
  ylab("Weight (g)") +
  xlab("Cross length (cm)") +
  labs(color='Species') +
  ggtitle(" Weight vs. cross length,\n each color refers to a species.") +
  scale_color_paletteer_d("pals::alphabet") +
  theme(text = element_text(size = 25))
## `geom_smooth()` using formula 'y ~ x'

Faktor menarik yang perlu diperhatikan adalah kemungkinan pengaruh yang tinggi dari pengukuran terkait dengan spesies Pike (berwarna ungu), pada kemiringan dan kebisingan regresi.

Hal ini mungkin terjadi karena tingginya nilai pengukuran berat dan panjang beberapa spesies ini, yang mungkin berbeda dari keseluruhan sampel. Spesies ini kemungkinan besar akan memiliki nilai rata-rata tertinggi untuk jarak juru masak yang terkait dengan pengamatannya, yang akan menunjukkan pengaruhnya terhadap model regresi. Untuk memeriksanya, nilai-nilai ini dapat diringkas untuk setiap spesies:

mod <- lm(Weight~Length3,data=fish)

augment(mod,data=fish) %>% 
  group_by(Species) %>% summarise(mean_cooksd=mean(.cooksd)) %>% 
  arrange(desc(mean_cooksd))
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 7 x 2
##   Species   mean_cooksd
##   <chr>           <dbl>
## 1 Pike         0.0428  
## 2 Smelt        0.0116  
## 3 Whitefish    0.00550 
## 4 Perch        0.00528 
## 5 Bream        0.00325 
## 6 Roach        0.00150 
## 7 Parkki       0.000504

Seperti yang diharapkan, pengamatan Pike memiliki rata-rata nilai jarak juru masak tertinggi. Spesies ini diikuti oleh Smelt (hijau tua di plot), yang memiliki konsentrasi nilai mendekati nol tinggi. Poin-poin ini, untuk kedua spesies, cukup jauh dari garis regresi.

Untuk memeriksa secara lebih teliti pengaruh spesies ini kepada regresinya, kita akan memodelkan variabel yang sama, dan menghitung \(R^2\) yang terkait dengan model baru:

fish %>% filter(!(Species %in% c("Pike","Smelt"))) %>% 
  lm(Weight~Length3,data=.) %>% summary() %>% {.$r.squared} %>% round(digits=2)
## [1] 0.91

Mengingat fakta bahwa nilai \(R^2\) untuk model sebelumnya adalah 0,85, ini dapat dianggap sebagai peningkatan yang bagus. Selain memerika nilai koefisien, kita dapat melihat peningkatan nilai kemiringan juga:

fish %>% filter(!(Species %in% c("Pike","Smelt"))) %>% 
  lm(Weight~Length3,data=.)
## 
## Call:
## lm(formula = Weight ~ Length3, data = .)
## 
## Coefficients:
## (Intercept)      Length3  
##     -656.48        34.14

Terakhir, pengaruh ini dapat dilihat secara visual dengan memplot setiap garis regresi secara bersamaan:

df_alt <- fish %>% filter(!(Species %in% c("Pike","Smelt")))
mod2 <- lm(formula = Weight ~ Length3, data = df_alt)


df_alt %>% ggplot(aes(y=Weight,x=Length3)) + geom_point(alpha=0) +
  geom_abline(intercept = -656.48, slope = 34.14, linetype=2,color="blue") +
  geom_abline(intercept = -490.40, slope = 28.46, linetype=2,color="red") +
  ylab("Weight (g)") +
  xlab("Cross length (cm)") +
  ggtitle("Regression lines for the models containing all species (red) and Pike and Smelt removed (blue).") +
  theme_bw() +
  theme(text = element_text(size = 25))

Dengan memvisualisasikan garis regresi, perubahan perilaku menjadi lebih jelas. Perubahan kemiringan ini menunjukkan korelasi yang lebih kuat antara variabel tanpa spesies. Fakta yang diperkuat dengan menghitung indeks korelasi pearson untuk setiap kasus dicari dengan:

tibble(df=c("df_original","df_sp_removed"),
       pearson_index=c(
         cor(fish$Weight,fish$Length3),
         cor(df_alt$Weight,df_alt$Length3)
       ))
## # A tibble: 2 x 2
##   df            pearson_index
##   <chr>                 <dbl>
## 1 df_original           0.923
## 2 df_sp_removed         0.954

Akhirnya, cara bagus lainnya untuk memeriksa bagaimana setiap spesies yang mungkin mempengaruhi regresinya adalah dengan memplot setiap baris secara terpisah untuk setiap spesies:

fish %>% ggplot(aes(y=Weight,x=Length3)) + geom_point(size=3) + 
  geom_smooth(method="lm", linetype=2) +
  ylab("Weight (g)") +
  xlab("Cross length (cm)") +
  coord_cartesian(ylim = c(0,1500)) + 
  ggtitle("Weight vs. crossed length, per species.") +
  theme_bw() +
  theme(text = element_text(size = 25)) +
  facet_wrap(~Species)
## `geom_smooth()` using formula 'y ~ x'

Pike adalah satu-satunya spesies dengan nilai bobot lebih tinggi dari 1500. Di sisi lain, Smelt hanya memiliki nilai yang mendekati nol, berarti memperkuat pengaruhnya kepada garis regresi. Informasi ini agak berlebihan dengan plot tempat semua spesies diambil sampelnya, tetapi ini adalah cara alternatif untuk melihat semua fakta-fakta ini.

Jenis wawasan ini menimbulkan pertanyaan, yaitu: spesies mana yang harus dipilih agar sesuai dengan model keseluruhan terbaik dengan data yang cukup? Pendekatan yang dapat membantu adalah dengan menghitung nilai \(R^2\) untuk model yang sama untuk setiap spesies, sementara juga memperhitungkan ukuran sampel dalam setiap kasus:

# previously used function with new data argument
rsquared_mod <- function(x,data) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = data)
    )[["r.squared"]]
  )
}

# explanatory variables
explanatory <- names(fish)[3:length(names(fish))]

# function to compute r squared per model per species

r_sp_x <- function(sp){
  data <- fish %>% filter(Species==sp)
  return( 
  tabela_r <- tibble(Species=rep(sp,5),
                     variables=explanatory,
                     r_squared=unlist(map(explanatory,rsquared_mod,data=data)))
  )
}

# list with data frames with r squared values per variable per species

list_r_sp_x <- map(unique(fish$Species),r_sp_x)

# above list as data frame
df_r_sp_x <- do.call("rbind",list_r_sp_x) 

# joining df with summarised original data to included number of observations per species

left_join( 
df_r_sp_x %>% filter(variables=="Length3") %>% arrange(desc(r_squared)),
fish %>% group_by(Species) %>% summarise(n=n()),
by="Species"
)
## `summarise()` ungrouping output (override with `.groups` argument)
## # A tibble: 7 x 4
##   Species   variables r_squared     n
##   <chr>     <chr>         <dbl> <int>
## 1 Whitefish Length3       0.993     6
## 2 Pike      Length3       0.958    17
## 3 Parkki    Length3       0.943    11
## 4 Perch     Length3       0.921    56
## 5 Smelt     Length3       0.898    14
## 6 Bream     Length3       0.897    35
## 7 Roach     Length3       0.842    20

Ketika dianalisis secara terpisah, Pike (yang memberikan nilai jarak juru masak yang tinggi untuk pengamatannya), memiliki nilai \(R^2\) yang tinggi untuk model tersebut. Namun penting untuk diperhatikan, bahwa jumlah observasi untuk sebagian besar spesies ini cukup rendah (di bawah 20). Hanya spesies Perch dan Beam yang memiliki lebih dari 20 pengamatan. Ikan bandeng yang memiliki nilai \(R^2\) tertinggi pada modelnya hanya memiliki 6 pengamatan. Jika idenya adalah untuk membuat model yang lebih umum, yang dilindungi oleh kemungkinan pencilan, mungkin masuk akal untuk memelihara spesies dengan jumlah pengamatan yang lebih tinggi atau mengambil sampel beberapa spesies lebih banyak.

2.1.3 Kolinearitas variabel

Seperti yang dinyatakan sebelumnya, semua pengukuran panjang menghasilkan nilai R2 yang tinggi disaat dikorelasikan dengan berat.

Mengingat bingkai data:

rsquared_mod <- function(x) {
  return(
    summary( 
    lm(paste0("Weight ~", x), data = fish)
    )[["r.squared"]]
  )
}


explanatory <- names(fish)[3:length(names(fish))]

tabela_r <- tibble(variables=explanatory,r_squared=unlist(map(explanatory,rsquared_mod)))

tabela_r %>% arrange(desc(r_squared))
## # A tibble: 5 x 2
##   variables r_squared
##   <chr>         <dbl>
## 1 Length3       0.852
## 2 Length2       0.844
## 3 Length1       0.839
## 4 Width         0.786
## 5 Height        0.525

Perilaku ini mungkin merupakan sebuah indikasi kolinearitas. Ini dapat diperiksa dengan menghitung faktor inflasi varians (vif) untuk setiap variabel, dalam model yang memperhitungkan semuanya:

fish %>% 
  select(-Species) %>%
  lm(Weight~.,data=.) %>%
  vif() %>% as.data.frame()
##                  .
## Length1 1681.49649
## Length2 2084.25783
## Length3  422.46825
## Height    14.57009
## Width     12.27536

Nilai vif yang tinggi berhubungan dengan pengukuran panjang, memiliki arti yaitu menunjukkan bahwa ini sebenarnya berlebihan untuk sebuah model, lalu model ini juga tidak membuat perbedaan besar yang mana digunakan untuk menjelaskan bobot nilai terebut.

Perilaku ini dapat divisualisasikan dengan memplot bobot dalam fungsi dua ukuran panjang (tidak peduli yang mana):

# model
mod3 <- lm(Weight~Length3+Length1, data=fish)

# 3d scatterplot
p <- plot_ly(data = fish, z = ~Weight, x = ~Length3, y = ~Length1, 
     opacity = 0.6, colorbar = list(title = "Predicted weight")) %>%
      add_markers()

# model coefficients

cf.mod <- coef(mod3)

# computing plane size

x1.seq <- seq(min(fish$Length3),max(fish$Length3),length.out=1000)
x2.seq <- seq(min(fish$Length1),max(fish$Length1),length.out=1000)

# computing plane

z <- t(outer(x1.seq, x2.seq, function(x,y) cf.mod[1]+cf.mod[2]*x+cf.mod[3]*y))


k <- p %>% # 3d scatterplot with regression plane
  add_surface(x = ~x1.seq, y = ~x2.seq, z = ~z, 
              showscale = TRUE,colorscale="Viridis") %>% 
  layout(scene = list(xaxis = list(title = "Crossed length (cm)"), 
                      yaxis = list(title = "Vertical length (cm)"),
                      zaxis = list(title = "Weight (g)")))

             
# need this to plot html widgets             

#library(IRdisplay)
#htmlwidgets::saveWidget(k, "m.html")
#display_html('<iframe src="m.html" width=100% height=450></iframe>')

Perilaku linier dekat antara kedua pengukuran panjang memperkuat kolinearitasnya: mengarahkan mouse pada setiap titik data akan mengungkapkan bahwa nilai x dan y di plot selalu memiliki nilai yang sama (dalam hal ini variabel respons kami berada di sumbu z). Juga, bidang yang diproyeksikan untuk model ini memprediksi nilai bobot yang tinggi untuk semua nilai panjang yang tinggi, sebesar apapun nilainya tersebut.

2.1.4 Kesimpulan

Secara keseluruhan, dapat disimpulkan bahwa nilai bobot dapat dijelaskan secara wajar dengan pengukuran panjang yang bersifat kolinear. Lebih dalam lagi, spesies dalam dataset memberikan kontribusi yang berbeda terhadap kualitas model regresi, dengan beberapa menghasilkan kebisingan dengan memiliki sedikit pengamatan atau sangat berpengaruh.

Sumber : https://www.kaggle.com/moraessaur/fish-market-simple-linear-regression

2.2 Non linear Regression

2.2.1 Required R packages

Package tidyverse untuk manipulasi dan visualisasi data yang mudah. Package caret untuk alur kerja machine learning yang mudah.

library(tidyverse)
library(caret)
theme_set(theme_classic())

2.2.2 Preparing the data

data <- structure(list(S = c(3.6, 1.8, 0.9, 0.45, 
                           0.225, 0.1125, 3.6, 
                           1.8, 0.9, 0.45, 0.225, 
                           0.1125, 3.6, 1.8, 0.9, 
                           0.45, 0.225, 0.1125, 0), 
                     v = c(0.004407692, 0.004192308, 
                           0.003553846, 0.002576923,
                           0.001661538, 0.001064286, 
                           0.004835714, 0.004671429, 
                           0.0039, 0.002857143,
                           0.00175, 0.001057143, 
                           0.004907143, 0.004521429, 
                           0.00375, 0.002764286,
                           0.001857143, 0.001121429, 
                           0)), 
                .Names = c("S", "v"), 
                class = "data.frame", 
                row.names = c(NA, -19L))
data
##         S           v
## 1  3.6000 0.004407692
## 2  1.8000 0.004192308
## 3  0.9000 0.003553846
## 4  0.4500 0.002576923
## 5  0.2250 0.001661538
## 6  0.1125 0.001064286
## 7  3.6000 0.004835714
## 8  1.8000 0.004671429
## 9  0.9000 0.003900000
## 10 0.4500 0.002857143
## 11 0.2250 0.001750000
## 12 0.1125 0.001057143
## 13 3.6000 0.004907143
## 14 1.8000 0.004521429
## 15 0.9000 0.003750000
## 16 0.4500 0.002764286
## 17 0.2250 0.001857143
## 18 0.1125 0.001121429
## 19 0.0000 0.000000000

Setelah kami memiliki data kita, kita dapat menggunakan paket drc untuk menyesuaikannya dengan kurva. Kita juga dapat menggunakan paket ggplot2 untuk memplot data juga.

library(drc)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:formattable':
## 
##     area
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 'drc' has been loaded.
## Please cite R and 'drc' if used for a publication,
## for references type 'citation()' and 'citation('drc')'.
## 
## Attaching package: 'drc'
## The following objects are masked from 'package:stats':
## 
##     gaussian, getInitial
library(ggplot2)
MM.model <- drm(v~S, data=data, fct=MM.2())
mmdf <- data.frame(S=seq(0,max(data$S),length.out=100))
mmdf$v <- predict(MM.model, newdata=mmdf)
ggplot(data, aes(x = S, y = v)) +
  theme_bw() +
  xlab("Concentration [mM]") +
  ylab("Speed [dE/min]") +
  ggtitle("Techvidvan Michaelis-Menten kinetics") +
  geom_point(alpha = 0.5) +
  geom_line(data = mmdf, 
            aes(x = S, y = v), 
            colour = "green")

Kita bisa melihat ringkasan model dengan menggunakan fungsi summary (). Ayo kita coba!

summary(MM.model)
## 
## Model fitted: Michaelis-Menten (2 parms)
## 
## Parameter estimates:
## 
##                 Estimate Std. Error t-value   p-value    
## d:(Intercept) 0.00542523 0.00011758  46.141 < 2.2e-16 ***
## e:(Intercept) 0.43980191 0.03051108  14.415 5.814e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error:
## 
##  0.0001727146 (17 degrees of freedom)

Kita juga bisa melakukan regresi dan memplotnya menggunakan fungsi nls (). Mari kita coba.

mm.model.nls <- nls(v~Vm*S/(K+S), data=data, 
                    start = list(K=max(data$v)/2, 
                                 Vm=max(data$v)))
summary(mm.model.nls)
## 
## Formula: v ~ Vm * S/(K + S)
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## K  0.4398016  0.0311612   14.11  8.1e-11 ***
## Vm 0.0054252  0.0001193   45.47  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0001727 on 17 degrees of freedom
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 4.666e-07

Dalam bab ini, kita belajar tentang regresi non-linier di R. Kita mempelajari apa itu regresi non-linier dan jenis model regresi apa yang dianggap non-linier. Kemudian kita melihat metode estimasi kemungkinan maksimum. Kita mempelajari lebih lanjut tentang regresi logistik, regresi Michaelis-Menten, dan model aditif umum. Terakhir, Kita juga mempelajari cara mengubah model non-linier menjadi model linier dan mengapa kami ingin melakukannya. Akhirnya, kita belajar bagaimana menerapkan model regresi non-linier di R.

3 Soal 3

Please read this link RPubs - Multi-Objective Optimization with NSGA-II (translate and submit to your Rpubs).

3.1 Optimisasi Muklti Objektif

3.1.1 Pembukaan

3.1.1.1 Tentang

Pengoptimalan itu sangat penting diberbagai bidang, termasuk dalam ilmu data. Sedangkan di bidang manufaktur, di mana setiap keputusan sangat penting untuk suatu proses dan keuntungan dari suatu organisasi, sering dilakukan optimasi, mulai dari jumlah setiap produk yang diproduksi, bagaimana suatu unit dijadwalkan untuk diproduksi, mendapatkan parameter proses yang terbaik atau optimal, dan juga perutean. Dalam ilmu data, kita sudah familiar dengan model tuning, dimana kita akan menyesuaikan model kita untuk meningkatkan performa dari sautu model. Algoritma pengoptimalan dapat membantu kita untuk mendapatkan performa model yang lebih baik.

Beberapa masalah memiliki tujuan yang bertentangan, atau trade-off, seperti penyeimbangan antara memaksimalkan pengembalian dengan meminimumkan resiko, atau memaksimalkan keuntungan sambil meminimumkan biaya produksi. Masalah seperti ini termasuk dalam Multi-Objective Optimization Problem (MOOP). Memecahkan MOOP lebih sulit dan membutuhkan pendekatan yang berbeda dari optimasi tujuan tunggal. Salah satu contoh algoritma untuk mengatasi masalah tersebut adalah Algoritma Evolusi Multi-Objektif, yang merupakan perpanjangan dari Algoritma Genetika untuk MOOP. Jika Anda masih awam dengan Algoritma Genetika, Anda bisa mengecek postingan saya sebelumnya. Salah satu dari Multi-Objective Evolutionary Algorithm yang paling populer adalah NSGA-II2. Banyak algoritma baru yang diturunkan dari metode ini. Kita akan belajar bagaimana menerapkannya untuk menyelesaikan MOOP dalam suatu masalah bisnis.

Posting ini didedikasikan untuk mempelajari caranya menyelesaikan masalah pengoptimalan multi-objektif dan membangun algoritme pengoptimalan dengan R.

3.1.1.2 Tujuan Pembelajaran

~Mempelajari cara menyelesaikan masalah pengoptimalan multi-objektif

~Mempelajari cara kerja NSGA-II

~Mempelajari implementasi NSGA-II

3.1.1.3 Library dan Setup

library(ranger)
## 
## Attaching package: 'ranger'
## The following object is masked from 'package:randomForest':
## 
##     importance
library(scales)
## 
## Attaching package: 'scales'
## The following objects are masked from 'package:formattable':
## 
##     comma, percent, scientific
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
library(plotly)
library(tidyverse)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:cowplot':
## 
##     stamp
## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
options(scipen = 999)
library(earth)
## Loading required package: Formula
## Loading required package: plotmo
## Loading required package: plotrix
## 
## Attaching package: 'plotrix'
## The following object is masked from 'package:scales':
## 
##     rescale
## Loading required package: TeachingDemos
## 
## Attaching package: 'TeachingDemos'
## The following object is masked from 'package:formattable':
## 
##     digits
## The following object is masked from 'package:plotly':
## 
##     subplot
library(kernlab)
## 
## Attaching package: 'kernlab'
## The following object is masked from 'package:scales':
## 
##     alpha
## The following object is masked from 'package:purrr':
## 
##     cross
## The following object is masked from 'package:ggplot2':
## 
##     alpha
library(nnet)
library(mice)
## 
## Attaching package: 'mice'
## The following object is masked from 'package:stats':
## 
##     filter
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
library(tidyverse)
library(kableExtra)
library(AppliedPredictiveModeling)
library(lpSolve)
library(tidyr)
library(dplyr)
library(ggplot2)
library(caret)
library(purrr)
library(rlang)
## 
## Attaching package: 'rlang'
## The following object is masked from 'package:data.table':
## 
##     :=
## The following objects are masked from 'package:purrr':
## 
##     %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
##     flatten_lgl, flatten_raw, invoke, list_along, modify, prepend,
##     splice
## The following object is masked from 'package:Metrics':
## 
##     ll
library(hrbrthemes)
## NOTE: Either Arial Narrow or Roboto Condensed fonts are required to use these themes.
##       Please use hrbrthemes::import_roboto_condensed() to install Roboto Condensed and
##       if Arial Narrow is not on your system, please see https://bit.ly/arialnarrow
library(GGally)
library(mlbench)

3.1.2 Optimasi Multi-Objektif: Konsep

3.1.2.1 Tujuan Ganda

Optimasi adalah proses untuk menjadi lebih baik. Optimasi bertujuan untuk menemukan solusi yang dapat memaksimalkan atau meminimumkan suatu fungsi tujuan tunggal yang sesuai dengan suatu kendala yang ada. Misalnya, kita ingin memaksimalkan sautu keuntungan sekaligus meminimumkan suatu biaya produksi. Jika masalah memiliki lebih dari satu tujuan, masalah dapat dimodifikasi menjadi masalah objektif tunggal, seperti menggabungkan nilai mereka dan memberi bobot kepada setiap fungsi tujuan. Namun, metode ini memiliki beberapa kelemahan:

Metode ini membutuhkan pengetahuan tentang suatu prioritas atau bobot dari setiap tujuan.

~ Fungsi agregat hanya memberikan solusi tunggal.

~ Pertukaran antara tujuan tidak dapat diketahui dan dievaluasi.

Solusi saya tidak akan tercapai jika ruang pencarian berbentuk cembung.

Sementara itu, dalam masalah optimasi multi objektif (MOOP), trade-off dapat diidentifikasikan karena algoritma MOOP memberikan banyak solusi di seluruh pareto front. Berbagai solusi juga memberikan pembuat keputusan memiliki beberapa alternatif untuk dipilih daripada hany satu solusi.

3.1.2.2 Pareto Front

Karena MOOP memiliki lebih dari satu fungsi objektif, kita tidak bisa langsung mengatakan apakah solusi A lebih baik daripada solusi B. Kita memiliki sesuatu yang disebut pareto-dominance.

Diberikan:

Sebuah vektor fungsi tujuan \(\bar f \bar {(x)} = [f_1 \bar{(x)}, f_2 \bar {(x)}, ..., f_k \bar {(x)}]\)

~ Dan solusi yang layak dari Ω

Permasalahan Optimasi Multi Objektif adalah mencari vektor $  bar x ϵ Ω $ yang mengoptimalkan fungsi vektor $  bar f  bar {(x)} $.

Dalam masalah minimisasi, vektor $  bar x $ dikatakan mendominasi $ bar y $ (dilambangkan  bar x ≺  bar y) jika:

~ $ f_i \ bar x ≤ f_i \ bar y $ untuk semua$ i $ fungsidi $ f_i $


~ Setidaknya ada satu $i $ sehingga $ f_i \ bar x <f_i \ bar y $

Dengan kata sederhana, jika kita memiliki solusi yang disebut solusi x, solusinya dikatakan mendominasi (lebih baik dari) solusi yang lain, yang disebut solusi y, jika semua fungsi tujuan dari solusi x kurang atau sama dari semua fungsi tujuan dari solusi y dan setidaknya ada satu fungsi tujuan dari solusi x yang lebih baik (memiliki nilai yang lebih rendah dalam minimumkan) daripada fungsi tujuan dari solusi y.

Suatu solusi disebut dominan jika memiliki setidaknya satu nilai fungsi tujuan yang lebih baik dari yang lain sedangkan memiliki nilai yang lebih baik atau sama pada fungsi tujuan yang lain. Lihatlah grafik di bawah ini.

Misalkan kita memiliki suatu solusi (disebut saja solusi A, B, C, dan D) dengan nilai fungsi tujuan pertama sebagai sumbu x (f1) dan fungsi tujuan kedua (f2) sebagai sumbu y. Kedua fungsi objektif tersebut merupakan masalah minimisasi, sehingga semakin dekat posisinya dengan kiri bawah maka semakin baik solusinya. Dari gambar di atas, kita dapat menyimpulkan beberapa hal:

Solusi B dikatakan tidak didominasi oleh solusi A, karena f1 dari A lebih kecil dari B tetapi f2 lebih besar dari B.

~ Karena solusi A dan B tidak didominasi masing-masing lain, mereka termasuk dalam grup / front yang sama.

~ Solusi C didominasi oleh solusi A dan B, karena f1 dari C lebih rendah dari A dan B, dengan f2 dari C lebih rendah dari B atau sama dengan A.

~ Karena solusi C didominasi oleh A dan B, itu milik bagian depan bawah (depan 2).

~ Solusi D didominasi oleh solusi A, B, dan C karena kedua nilai obyektif tersebut lebih buruk daripada keduanya.

~ Karena solusi D didominasi oleh A, B, dan C, solusi ini termasuk dalam front terendah (front 3).

~ True Pareto Front adalah seperangkat solusi optimal yang ingin kami capai dalam MOOP.

3.1.2.3 Algoritma Multi-Objectif Evolusionaritas

Multi-Objective Evolutionary Algorithm (MOEA) adalah aplikasi dari algoritma evolusioner, seperti algoritma genetik untuk menangani MOOP. Sudah ada banyak sekali penelitian di bidang ini, yang menghasilkan banyak variasi dari suatu algoritma. NSGA-II merupakan salah satu pilihan yang populer dan menjadi suatu inspirasi untuk penelitian selanjutnya.

Sebelum kita melanjutkan, mari kita tinjau cara kerja Algoritma Genetika.

Akan ada populasi kromosom atau larutan yang dipilih secara acak.

Nilai kebugaran atau nilai fungsi tujuan dari setiap solusi dihitung.

Dari populasi tersebut, dua solusi akan dipilih sebagai solusi induk, baik melalui pemilihan turnamen atau metode pemilihan lainnya.

Solusi yang dipilih akan disilangkan untuk menciptakan solusi baru yang disebut solusi anak.

Solusi anak dapat berubah karena mutasi acak (yang memiliki kemungkinan sangat rendah untuk terjadi).

Di akhir iterasi, populasi baru akan dipilih dari solusi induk atau populasi awal dan solusi anak berdasarkan nilai fitness.

Selama kriteria penghentian tidak terpenuhi, biasanya dalam hal jumlah iterasi, algoritme akan terus melakukan iterasi.

Solusi terbaik atau optimal adalah solusi dengan nilai kebugaran optimal.

Untuk masalah multi-objektif, NSGA-II bekerja dengan alur kerja yang sama dengan proses di atas. Bedanya, pada setiap akhir iterasi, dipilih populasi baru untuk iterasi berikutnya dengan menggunakan metode sortir yang tidak didominasi. Solusinya dibagi menjadi beberapa front berdasarkan status dominasinya terhadap satu sama lain. Jika kita memiliki ukuran populasi awal sebesar 10, maka 10 individu teratas dari populasi induk dan anak dari depan yang lebih baik akan dipilih untuk iterasi berikutnya. Juga, untuk mengontrol keragaman solusi, jarak crowding, yang sesuai dengan jarak Manhattan dari dua solusi yang bertetangga untuk dua tujuan, dari setiap solusi dihitung3. Burung dengan jarak kokok terbesar akhirnya dipilih.

3.1.3 Studi Kasus

3.1.3.1 Optimisasi Portofolio dengan 3 Objektif

Masalah ini direplikasi dari studi Anagnostopolous dan Mamanis 4, yang melibatkan 3 fungsi objektif. Studi ini menggunakan MOEA untuk pemilihan portofolio dan pengoptimalan dalam manajemen investasi.

Masalah pengoptimalan portofolio berkaitan dengan pengelolaan portofolio aset yang meminimalkan tujuan risiko yang menjadi kendala untuk menjamin tingkat pengembalian tertentu. Salah satu prinsip fundamental dari suatu investasi keuangan adalah diversifikasi, dimana investor mendiversifikasi investasinya menjadi berbagai jenis aset. Diversifikasi portofolio meminimalkan eksposur investor terhadap risiko, dan memaksimalkan pengembalian portofolio. Dalam beberapa kasus, sangat sulit bagi pengambil keputusan untuk mengetahui sebelumnya jumlah optimal sekuritas yang harus dimasukkan ke dalam portofolionya tanpa memeriksa semua tradeoff antara risiko, pengembalian dan kardinalitas portofolio. Bagaimanapun, seorang investor mungkin kehilangan solusi penting dengan tradeoff besar antara tujuan ketika dia dipaksa untuk memperbaiki jumlah aset dalam portofolio sebelumnya.

Fungsi tujuan terdiri dari:

  1. Memaksimalkan Pengembalian

Hasil

\[ max μ (x) = ∑ _ {i = 1} ^ n x_i μ_i \] \(mu(x)\): rata-rata total pengembalian \(x_i\): berat aset i \(mu_i\): rata-rata pengembalian aset i

  1. Meminimalkan Risiko

Risiko portofolio diwakili oleh varians pengembalian.

\[mindp (x) = ∑_ {i_1} ^ N ∑_ {j_1} ^ N x_i x_j σ_ {ij} \] \(p(x)\): risiko total portofolio

\(σ_{ij}\): kovariansi antara aset \(i\)dan \(j\)

  1. Meminimalkan Jumlah Aset

Jumlah aset adalah diminimalkan. Jumlah aset adalah jumlah dari semua aset yang berbobot > 0.

\[min kartu(x)= ∑_{i_1}^N 1_{x> 0}\] \(kartu (x)\): jumlah aset

Subjek Ke

\[∑_{i_1}^N x_i=1\] \[0≤x_i≤1\]

3.1.3.1.1 Import Data

Untuk mendapatkan nama perusahaan yang lebih jelas, ayo mengimport data Ticker Symbol dan Security.

#library(dplyr)
#library(tidyr)

#securities <- read.csv("securities.csv", sep = ",", header = TRUE)

#securities <- securities %>% 
#select(Ticker.symbol, Security)

#securities1<- rename(securities,
#       symbol = Ticker.symbol,
#       name = Security)

Data diperoleh dari Bursa Efek New York di Kaggle (https://www.kaggle.com/dgawlik/nyse). Kami hanya akan menggunakan data dari Januari hingga Maret 2015 untuk ilustrasi.

#library(tidymodels)
#library(ranger)
#library(scales)
#library(plotly)
#library(tidyverse)
#library(lubridate)

#nyse <- read.csv("prices.csv")

#nyse <- nyse %>% 
#  mutate(date = ymd(date)) %>% 
#  filter(year(date) == 2015,
#         month(date) %in% c(1:3))

#head(nyse)

Gabungkan mereka dengan data price

#nyse <- nyse %>% 
#   left_join(securities1, by = c("symbol")) %>% 
#   select(date, symbol, name, everything())

#head(nyse)

Biarkan limit dari jumlah stok yang tersedia hanya dipilih secara acak menjadi 30 stok. Terlalu banyak pilihan akan membuat waktu berjalan semakin lama. Sejak ini hanyalah sebuah pengenalan atau permulaan, Saya akan menguranginya.

#set.seed(13)
#selected_stock <- sample(nyse$symbol, 30)

#nyse <- nyse %>% 
#  filter(symbol %in% selected_stock)
#head(nyse)
3.1.3.1.2 Perhitungan Pengembalian
#nyse <- nyse %>%   
#  rename(price = close) %>% 
#  select(date, symbol, name, price) %>% 
#  group_by(symbol, name) %>%
#  mutate(price_prev = lag(price),
#         returns = (price - price_prev)/price_prev) %>% 
#  slice(-1) %>% 
#  ungroup()

#head(nyse)

Mari menghitung rata-rata pengembalian dari masing-masing stok

#mean_stock <- nyse %>% 
#  group_by(symbol) %>% 
#  summarise(mean = mean(returns)) %>% 
#  arrange(desc(mean))

#head(mean_stock)

Nilai \(R_f\) diperoleh dari tingkat bunga terbaru pada tagihan Treasury AS tiga bulan. Karena datanya dari tahun 2016, kita akan menggunakan data dari tahun 2015 (Gunakan data dari 27 Maret 2015), yaitu 0,04%. Tarif diperoleh dari https://ycharts.com/indicators/3_month_t_bill.

rf <- 0.04/100
3.1.3.1.3 Matriks Kovariansi diantara Portofolio

Menghitung nilai matriks kovarians diantara portofolio. Pertama, kita perlu untuk memisahkan pengembalian dari masing-masing portofolio kepada kolom yang ada untuk menyebarkan mereka.

#nyse_wide <- nyse %>%
#  pivot_wider(id_cols = c(date, symbol), names_from = symbol, values_from = returns) %>% 
#  select(-date)

# Create Excess Return
#for (i in 1:n_distinct(nyse$symbol)) {
#  nyse_wide[,i]<- nyse_wide[,i] - as.numeric(mean_stock[i,2])
#}
  
#head(nyse_wide)

Membuat matriks kovarians.

#nyse_cov <- cov(x = nyse_wide)
3.1.3.1.4 Membuat Fungsi Fitness

Mari membuat fungsi fitnes.

#fitness <- function(x){
  # Assign weight for each stocks
#  weight_stock <- x
  
 # Calculate the total returns
# f1 <- numeric()
 #for (i in 1:n_distinct(nyse$symbol)) {
#   f1[i] <- weight_stock[i]*mean_stock$mean[i]
# }
# mean_return <- sum(f1) - 1e9 * (round(sum(weight_stock),10)-1)^2 
 
 # Calculate the total risk
# f2 <- numeric()
# for (i in 1:n_distinct(nyse$symbol)) {
#   f3 <- numeric()
   
#   for (j in 1:n_distinct(nyse$symbol)) {
#    f3[j] <- weight_stock[i]*weight_stock[j]*nyse_cov[i,j]
#   }
# f2[i] <- sum(f3)
# }
# risk <- sum(f2) + 1e9 * (round(sum(weight_stock),10)-1)^2
 
 # Calculate the number of assets
# card <- length(weight_stock[weight_stock > 0]) 
 
# return(c(-mean_return, risk, card))
#}
3.1.3.1.5 Menjalankan Algoritma

Algoritma NSGA-II bisa dijalankan dengan fungsi nsga() dari package nsga2R.

#library(nsga2R)
Berikut adalah beberapa parameter untuk fungsi:
fn:fungsi kebugaran untuk diminimalkan

~ varNo:Jumlah variabel keputusan.

~ objDim:Jumlah fungsi tujuan

lowerBounds:Batas bawah masing-masing variabel keputusan.

~upperBounds:Batas atas masing-masing variabel keputusan.

~ popSize:ukuran populasi.

~ generations:Jumlah generasi.

Kita akan menjalankan algoritma untuk 1000 iterasi, dengan ukuran populasi sebesar 200. Probabilitas mutasi adalah 0,2 sedangkan probabilitas crossovernya sebesar 0.8.

#set.seed(123)
#n_asset <- n_distinct(nyse$symbol)
#finance_optim <- nsga2R(fn = fitness, varNo = n_asset, objDim = 3, generations = 1000,
#                        mprob = 0.2, popSize = 200, cprob = 0.8,
#                        lowerBounds = rep(0, n_asset), upperBounds = rep(1, n_asset))

Untuk mengecek solusi depan pareto, Anda bisa melihat dari nilai objectives dan nilai subset hanya dengan pareto front rank dari 1. Mari melihat berapa banyak perbedaan dari solusi optimal yang dibuat.

#finance_optim$objectives[finance_optim$paretoFrontRank == 1, ] %>% 
#   matrix(ncol = 3) %>% 
#   as.data.frame() %>% 
#   distinct() %>% 
#   mutate(V1 = round(-V1, 5),
 #         V2 = round(V2, 5)) %>% 
#   rename(`Total Return` = V1,
#          Risk = V2,
#          `Number of Assets` = V3)

Hanya ada 1 solusi.

Mari mengecek total beratnya

#sum(finance_optim$parameters[1, ])

Jumlah beratnya 1, sehingga kita dapat yakin bahwa solusi kita tidak melanggar suatu kendala.

Jika Anda ingin mendapatkan bobot untuk setiap aset, Anda dapat menggunakan parameter objekdari finance_optim. Mari kita lihat parameter untuk solusi pertama.

#data.frame(symbol = unique(nyse$symbol),
#           name = unique(nyse$name),
#           weight = finance_optim$parameters[1, ]) %>% 
#   mutate(return = round(weight * mean_stock$mean, 5),
#          weight = round(weight,5))

3.1.3.2 Nilai Kelayakan Kendaraan

Masalah kelayakan kecelakaan kendaraan adalah masalah multi-tujuan di mana tingkat keselamatan kecelakaan kendaraan yang dapat dioptimalkan. Tingkat keselamatan yang lebih tinggi berarti seberapa baik kendaraan dapat melindungi penghuninya dari efek kecelakaan frontal.

Studi ini mempertimbangkan uji tabrak frontal lebar penuh dan uji tabrak offset-frontal. Dampak lebar penuh biasanya menghasilkan perlambatan yang tinggi. Diketahui bahwa perlambatan yang tinggi dapat menyebabkan cedera biomekanik yang serius bagi penghuninya sementara uji offset menuntut integritas struktural kendaraan. Untuk meningkatkan kelayakan kecelakaan secara keseluruhan dari struktur depan kendaraan, kedua skenario tabrakan ini harus dipertimbangkan secara bersamaan. Ada lima variabel keputusan yang merepresentasikan ketebalan member tulangan di sekeliling bagian depan mobil, seperti terlihat pada gambar di bawah ini.

Berikut ini adalah 3 tujuan yang perlu dicapai agar crashworthiness optimal:

\[minimize f_i (x) = 1640.2823+2.3573285 x_1 + 2.3220035 x_2 + 4.5688768 x_3 +7.7213633 x_4 + 4.4559504 x_5\] \[ minimize f_2(x) = 6.5856+1.15 x_1 − 1.0427 x_2 + 0.9738 x_3 + 0.8364 x_4 − 0.3695 x_1 x_4 + 0.0861 x_1 x_5 + 0.3628 x_2 x_4 − 0.1106 x1^2 − 0.3437 x_3^2 + 0.1764 x_4^2\]

\[minimize f_3 (x) = −0.0551 + 0.0181 x_1 + 0.1024 x_2 + 0.0421 x_3 − 0.0073 x_1 x_2 + 0.024 x_2 x_3 − 0.0118 x_2 x_4 − 0.0204 x_3 x_4 − 0.008 x_3 x_5 − 0.0241 x_2^2 + 0.0109 x_4^2\]

Dibatasi dengan \[ 1≤ x_i ≤3\]

\(f_1\) : Nilai massa dari kendaraan.

\(f_2\) : Integrasi akselerasi tabrakan dalam tabrakan frontal penuh

\(f_3\) : Intrusi toe-board dalam kecelakaan offset-frontal 40%

3.1.3.2.1 Menghitung Fungsi Fitness

MAri menghitung fungsi fitnessnya.

fitness <- function(x){
  # Calculate f1(x)
  f1 <- 1640.2823 + 2.3573285*x[1] + 2.3220035*x[2] +4.5688768*x[3] +
    7.7213633*x[4] + 4.4559504*x[5]
  
  # Calculate f2(x)
  f2 <- 6.5856 + 1.15*x[1] - 1.0427*x[2] + 0.9738*x[3] + 0.8364*x[4] - 
    0.3695*x[1]*x[4] + 0.0861*x[1]*x[5] + 0.3628*x[2]*x[4] - 
    0.1106*x[1]^2 - 0.3437*x[3]^2 + 0.1764*x[4]^2
  
  # Calculate f3(x)
  f3 <- -0.0551 + 0.0181*x[1] + 0.1024*x[2] + 0.0421*x[3] - 
    0.0073*x[1]*x[2] + 0.024*x[2]*x[3] - 0.0118*x[2]*x[4] -
    0.0204*x[3]*x[4] - 0.008*x[3]*x[5] - 0.0241*x[2]^2 + 0.0109*x[4]^2
  
  return(c(f1,f2,f3))
}
3.1.3.2.2 Jalankan Algoritmanya

Kita akan menjalankan NSGA-II dengan 5000 generasi dan besar populasinya 100.

#set.seed(123)
#car_optim <- nsga2R(fn = fitness, varNo = 5, objDim = 3, generations = 5000, popSize = 100, 
                        #lowerBounds = rep(1, 5), upperBounds = rep(3, 5))

Mari melihat solusi optimal: semua yang termasuk dengan front pertama.

#car_result <- car_optim$objectives %>% 
# as.data.frame() %>% 
 # bind_cols(as.data.frame(car_optim$parameters)) %>% 
 # filter(car_optim$paretoFrontRank == 1) %>% 
 # mutate_all(.funs = function(x){round(x,3)}) %>%
 #distinct() %>%  
  #rename("f1" = V1, "f2" = V2, "f3" = V3,
         #"x_1" = V11, "x_2" = V21, "x_3" = V31, "x_4" = V4, "x_5" = V5)

#car_result

Seperti yang dapat kita lihat, kita dapat mencapai lebih dari 1 nilai obyektif yang serupa, yang menunjukkan bahwa semua solusi tersebut tidak didominasi satu sama lain: tidak ada solusi yang lebih baik dari yang lain. Mari kita periksa 3 solusi pertama.

#car_result[1:3, ]

Apakah solusi pertama lebih baik dari yang kedua? Ingatlah bahwa semuanya adalah masalah minimisasi.

Nilai f1 dari solusi pertama lebih tinggi dari solusi kedua

Nilai f2 dari solusi pertama lebih rendah dari solusi kedua

Nilai f3 dari solusi pertama lebih tinggi dari solusi kedua

Solusi kedua dapat mendominasi atau lebih baik dari solusi pertama jika nilai f2 lebih rendah dari solusi pertama. Karena f2 solusi kedua lebih tinggi, solusi pertama dan kedua tidak saling mendominasi, sehingga keduanya sama-sama baik.

Bagaimana dengan solusi pertama dengan yang ketiga?

Nilai f1 dari solusi pertama lebih rendah dari solusi ketiga

Nilai f2 dari solusi pertama lebih rendah dari solusi ketiga

Nilai f3 dari solusi pertama lebih tinggi dari solusi ketiga

Solusi pertama dapat mendominasi solusi ketiga jika nilai f3 lebih rendah dari solusi ketiga. Karena tidak, maka solusi pertama dan solusi ketiga tidak saling mendominasi.

Konsekuensinya adalah bahwa pembuat keputusan dapat memutuskan solusi mana yang akan dipilih dan dilihat trade-off antara tujuan. Jika pengambil keputusan memprioritaskan massa yang lebih rendah (f1), mereka harus memilih solusi yang memiliki nilaif1 terendah.

#car_result %>% 
#  filter(f1 %in% c(min(f1), max(f1)))

Jika Anda ingin memilih solusi dengan massa terendah, Anda dapat memilih baris pertama dari tabel di atas. Jika Anda ingin memilih solusi dengan massa terbesar, Anda dapat memilih baris kedua. Berdasarkan tabel di atas, jika kita memilih massa yang lebih rendah, itu akan menghasilkan intrusi f3 yang lebih tinggi. Pertukaran ini harus dipertimbangkan oleh pembuat keputusan.

Kita dapat memvisualisasikan nilai obyektif.

#plot_ly(car_result, x = ~f1, y = ~f2, z = ~f3) %>%
# add_markers() %>%
# layout(scene = list(xaxis = list(title = 'Mass'),
#                  yaxis = list(title = 'Collision Acceleration'),
#                zaxis = list(title = 'Intrusion')))

3.1.4 Konklusi

Multi-Objective Optimization Problem (MOOP) harus diselesaikan secara berbeda dari masalah tujuan tunggal. MOEA, seperti NSGA-II, dapat digunakan untuk menyelesaikan MOOP. Masih banyak algoritma yang dapat bekerja lebih baik dan meningkatkan beberapa aspek NSGA-II6. Kurangnya pengembangan paket untuk MOEA di R menjadi suatu perhatian. Kita mungkin perlu membangun fungsi kita sendiri jika kita tidak puas dengan NSGA-II.

3.1.5 Referensi

  1. http://rpubs.com/Argaadya/550805

2.http://www.dmi.unict.it/mpavone/nc-cs/materiale/NSGA-II.pdf

3.https://www.springer.com/gp/book/9783319521558

4.https://www.sciencedirect.com/science/article/abs/pii/S0305054809002275

5.https://www.researchgate.net/profile/Qing_Li24/publication/225123240_Multiobjective_optimization_for_crash_safety_design_of_vehicles_using_stepwise_regression_model/links/02e7e52cae2e11338d000000.pdf

6.http://i2pc.es/coss/Docencia/SignalProcessingReviews/Zhou2011.pdf