Pada tulisan sebelumnya, saya telah membahas introduction terhadap spiral optimization optimization (SOA) dan bagaimana cara kerjanya untuk menyelesaikan permasalahan pencarian minimasi, maksimasi, dan pencarian akar.
Sekarang saya akan mengangkat topik bagaimana menyelesaikan mixed integer programming (baik linear dan non linear) menggunakan SOA.
Salah satu trik yang bisa dilakukan agar SOA bisa menyelesaikan mixed integer programming adalah dengan mengubah constrained optimization problem menjadi unconstrained optimization problem kemudian memanfaatkan penalty constant.
Misal suatu permasalahan MILP atau MINLP bisa ditulis secara umum sebagai berikut:
\[\min_{x \in \mathbb{R}^n} f(x)\]
\[\text{subject to: } g_i(x) = 0, i = 1,2,..,M\]
\[\text{and } h_j(x) \leq 0,i = 1,2,..,N\]
\[x = (x_1,x_2,...,x_n)^T \in \mathbb{N}\]
Bentuk di atas bisa kita ubah menjadi:
\[F(x,\alpha,\beta) = f(x) + \sum_{i=1}^M \alpha_i g_i^2(x) + \sum_{j = 1}^N \beta_j (\max{(h_i(x),0)})^2\]
dimana \(\alpha,\beta\) merupakan penalty constant yang bisa dibuat sangat besar.
SOA relatif mudah untuk dituliskan dalam bentuk algoritma bahasa pemrograman manapun. Tapi ada satu hal yang bisa menjadi batu ganjalan dalam menuliskan algoritmanya. Apa itu? Yaitu pendefinisian matriks rotasi untuk masalah dengan n-dimensi.
Bentuk umum dari matriks rotasi adalah sebagai berikut:
\[R^{(n)} (\theta_{1,2},\theta_{1,3},..,\theta_{n,n-1}) = \prod_{i=1}^{n-1} \left ( \prod_{j=1}^i R^{(n)}_{n-i,n+1-j} (\theta_{n-i,n+1-j}) \right )\]
Perhatikan bahwa perkalian matriks rotasi yang dilakukan adalah cross product.
Alasan: Rotasi tidak mengubah norm suatu vektor.
Berikut adalah function yang digunakan untuk membuat rotation matrix:
buat_rot_mat = function(theta,n){
# buat template sebuah matriks identitas
temp_mat = matrix(0,ncol = n,nrow = n)
diag(temp_mat) = 1
# buat matriks identitas terlebih dahulu
mat_rot = temp_mat
for(i in 1:(n-1)){
for(j in 1:i){
temp = temp_mat
idx = n-i
idy = n+1-j
# print(paste0("Matriks rotasi untuk ",idx," - ",idy,": DONE"))
temp[idx,idx] = cos(theta)
temp[idx,idy] = -sin(theta)
temp[idy,idx] = sin(theta)
temp[idy,idy] = cos(theta)
# assign(paste0("M",idx,idy),temp)
mat_rot = mat_rot %*% temp
mat_rot = mat_rot
}
}
return(mat_rot)
}Oke, untuk contoh kasus pertama saya akan gunakan persoalan yang pernah saya tulis di blog sebelumnya.
Cari \((x,y)\) yang memaksimalkan \(7000 x + 12000 y\) dengan constraints sebagai berikut:
\[4x + 20y \leq 1960 \\ x + y \leq 250 \\ x,y \geq 0 \\ x,y \in \mathbb{Z}^+\]
Langkah pertama yang harus dilakukan adalah mengubah constrained problem menjadi unconstrained problem.
rm(list=ls())
f1 = function(x,y){-7000*x - 12000*y}
h1 = function(x,y){4*x + 20*y - 1960}
h2 = function(x,y){x + y - 250}
beta = 10^25
f = function(x,y){
el_1 = f1(x,y)
el_2 = beta * (max(h1(x,y),0))^2
el_3 = beta * (max(h2(x,y),0))^2
return(el_1 + el_2 + el_3)
}Fungsi \(f\) adalah bentuk unconstrained problem dari bentuk awalnya.
Berikut ini adalah function SOA yang dipakai untuk menyelesaikan masalah ini:
soa_mrf_ip = function(
N, # banyak titik
x1_d, # batas bawah x1
x1_u, # batas atas x1
x2_d, # batas bawah x2
x2_u, # batas atas x2
rot, # berapa banyak rotasi
k_max, # iterasi maks
r){ # berapa rate konstraksi
# N pasang titik random di selang [a,b] di R2
x1 = runif(N,x1_d,x1_u)
x2 = runif(N,x2_d,x2_u)
#x2 = 250 - x1
# hitung theta
theta = 2*pi / rot
# definisi matriks rotasi
A = matrix(c(cos(theta),-sin(theta),
sin(theta),cos(theta)),
ncol = 2,byrow = T)
# bikin data frame
temp =
data.frame(x1,x2) %>%
rowwise() %>%
mutate(f = f(round(x1,0),round(x2,0))) %>%
ungroup()
# proses iterasi
for(i in 1:k_max){
# mencari titik x* dengan max(f)
f_min =
temp %>%
# memastikan titik ada di D
filter(x1 >= x1_d & x1 <= x1_u) %>%
filter(x2 >= x2_d & x2 <= x2_u) %>%
# mencari titik min fungsi
filter(f == min(f))
# definisi pusat rotasi
pusat = c(f_min$x1[1],f_min$x2[1])
for(j in 1:N){
# kita akan ambil titiknya satu persatu
x0 = c(temp$x1[j],temp$x2[j])
# proses rotasi dan konstraksi terhadap pusat x*
xk = A %*% (x0-pusat) # diputar dengan x_bin sebagai pusat
xk = pusat + (r * xk)
# proses mengembalikan nilai ke temp
temp$x1[j] = xk[1]
temp$x2[j] = xk[2]
}
# hitung kembali nilai f(x1,x2)
temp =
temp %>%
rowwise() %>%
mutate(f = f(round(x1,0),round(x2,0))) %>%
ungroup()
}
# proses output hasil
output =
temp[N,] %>%
filter(f == min(f)) %>%
mutate(x1 = round(x1,0),x2 = round(x2,0),f = f(x1,x2))
return(output)
}Sekarang kita akan selesaikan dengan beberapa parameter berikut:
soa_mrf_ip(100, # banyak titik
0, # batas bawah x1
250, # batas atas x1
0, # batas bawah x2
250, # batas atas x2
20, # berapa banyak rotasi
900, # iterasi maks
0.8)## # A tibble: 1 x 3
## x1 x2 f
## <dbl> <dbl> <dbl>
## 1 190 60 -2050000Cari \(x_1,x_2,x_3\) yang memaksimalkan \(7x_1 + 3x_2 + x_3\), dengan constraints sebagai berikut:
\[6x_1 + 4x_2 + 5x_3 \leq 60 \\ 8x_1 + x_2 + 2x_3 \leq 80 \\ 9x_1 + x_2 + 7x_3 \leq 70 \\ x_3 \geq 0 \\ x_1,x_2 \in \mathbb{Z}^+\]
Masalah di atas termasuk ke dalam mixed integer linear programming karena ada dua variabel integer dan satu kontinu.
Berikut adalah modifikasi menjadi unconstrained optimization problem:
rm(list=ls())
f1 = function(x1,x2,x3){
7*x1 + 3*x2 + x3
}
h1 = function(x1,x2,x3){6*x1 + 4*x2 + 5*x3 - 60}
h2 = function(x1,x2,x3){8*x1 + x2 + 2*x3 - 80}
h3 = function(x1,x2,x3){9*x1 + x2 + 7*x3 - 70}
beta = 10^15
f = function(x1,x2,x3){
el_1 = - f1(x1,x2,x3)
el_2 = beta * (max(h1(x1,x2,x3),0))^2
el_3 = beta * (max(h2(x1,x2,x3),0))^2
el_4 = beta * (max(h3(x1,x2,x3),0))^2
return(el_1 + el_2 + el_3 + el_4)
}Berikut adalah function yang digunakan:
soa_mrf_ip_3_var = function(
N, # banyak titik
x1_d, # batas bawah x1
x1_u, # batas atas x1
x2_d, # batas bawah x2
x2_u, # batas atas x2
x3_d, # batas bawah x3
x3_u, # batas atas x3
rot, # berapa banyak rotasi
k_max, # iterasi maks
r){ # berapa rate konstraksi
# N pasang titik random di selang [a,b] di R3
x1 = runif(N,x1_d,x1_u)
x2 = runif(N,x2_d,x2_u)
x3 = runif(N,x3_d,x3_u)
# hitung theta
theta = 2*pi / rot
# definisi matriks rotasi
R12 = matrix(c(cos(theta),-sin(theta),0,
sin(theta),cos(theta),0,
0,0,1),
ncol = 3,byrow = T)
R13 = matrix(c(cos(theta),0,-sin(theta),
0,1,0,
sin(theta),0,cos(theta)),
ncol = 3,byrow = T)
R23 = matrix(c(1,0,0,
0,cos(theta),-sin(theta),
0,sin(theta),cos(theta)),
ncol = 3,byrow = T)
# bikin data frame
temp =
data.frame(x1,x2,x3) %>%
rowwise() %>%
mutate(f = f(round(x1,0),
round(x2,0),
x3)) %>%
ungroup()
# proses iterasi
for(i in 1:k_max){
# mencari titik x* dengan max(f)
f_min =
temp %>%
# memastikan titik ada di D
filter(x1 >= x1_d & x1 <= x1_u) %>%
filter(x2 >= x2_d & x2 <= x2_u) %>%
filter(x3 >= x3_d & x3 <= x3_u) %>%
# mencari titik max fungsi
filter(f == min(f))
# definisi pusat rotasi
pusat = c(f_min$x1[1],f_min$x2[1],f_min$x3[1])
for(j in 1:N){
# kita akan ambil titiknya satu persatu
x0 = c(temp$x1[j],temp$x2[j],temp$x3[j])
# proses rotasi dan konstraksi terhadap pusat x*
# diputar dengan x_bin sebagai pusat
xk = (R23 %*% (R13 %*% R12)) %*% (x0-pusat)
xk = pusat + (r * xk)
# proses mengembalikan nilai ke temp
temp$x1[j] = xk[1]
temp$x2[j] = xk[2]
temp$x3[j] = xk[3]
}
# hitung kembali nilai f(x1,x2,x3)
temp =
temp %>%
rowwise() %>%
mutate(f = f(round(x1,0),round(x2,0),x3)) %>%
ungroup()
}
# proses output hasil
output =
temp[N,] %>%
filter(f == max(f)) %>%
mutate(x1 = round(x1,0),x2 = round(x2,0),x3 = x3,
f = f1(x1,x2,x3))
return(output)
}Berikut solusinya:
soa_mrf_ip_3_var(
20, # banyak titik
0, # batas bawah x1
20, # batas atas x1
0, # batas bawah x2
20, # batas atas x2
0, # batas bawah x3
20, # batas atas x3
10, # berapa banyak rotasi
100, # iterasi maks
.9
)## # A tibble: 1 x 4
## x1 x2 x3 f
## <dbl> <dbl> <dbl> <dbl>
## 1 7 4 0.400 61.4