TUGAS 2 PRAKTIKUM KOMPUTASI STATISTIKA I

Nomor 1

The Laplace distribution, also known as the double exponential distribution, is a continuous probability distribution frequently used in situations where data exhibit a sharp central peak and heavier tails compared to the Normal distribution. It captures scenarios in which observations tend to cluster tightly around a central value but still allow for relatively large deviations, making it useful in fields such as signal processing, robust statistics, and Bayesian modeling.

It is defined by two parameters:

• µ : location

• b > 0 : scale

The probability density function (PDF) is: \[ f(x \mid \mu, b) = \frac{1}{2b} \exp\left( -\,\frac{|x - \mu|}{b} \right). \]

For observations x1, x2, . . . , xn, the likelihood function is: \[ L(\mu, b) = \left( \frac{1}{2b} \right)^{n} \exp\left( -\,\frac{1}{b} \sum_{i=1}^{n} |x_i - \mu| \right). \]

and the log-likelihood is: \[ \ell(\mu, b) = -n \ln(2b) - \frac{1}{b} \sum_{i=1}^{n} |x_i - \mu|. \]

Tasks

1. Sampling with Markov Chain Monte Carlo

Implement the Metropolis–Hastings algorithm to sample from the Laplace distribution with µ = 2 and b = 1.

1. Generate a Markov chain of length ntotal = 5000.

2. Discard the first n1 = 2000 samples as burn-in.

3. Keep the remaining n2 = 3000 samples.

4. Report the acceptance rate of your algorithm. Acceptance rate can be computed as:

\[ \alpha(x_t, x^*) = \min\left(1, \frac{\pi(x^*)}{\pi(x_t)}\right) \] Notes:

• Use initial state of x0 = 0.0 and random state of 67.

• Use Normal distribution with µ = 0 and \(σ^2\) = 1 as the proposal (comparing) distribution.

2. Comparing the MCMC Samples and Theoretical Distribution

Using the kept 3000 samples before:

1. Create a histogram of the MCMC sample and overlay the theoretical Laplace PDF.
2. Plot the empirical CDF of the MCMC sample and overlay the theoretical CDF.
3. Compute the sample mean and variance of the MCMC sample.
4. Compare these values with the theoretical mean and variance of Laplace(2, 1), and briefly comment on the agreement.

3. Estimating the Parameter of Laplace Distribution

Treat the kept MCMC samples as observed data.

1. Estimate µ and b using Newton Raphson method. You may estimate the parameter one by one separately or both at once.

2. Report the estimated parameters ˆµ and ˆb.

3. Compare your estimates with the true values µ = 2, b = 1 and provide a short interpretation.

Notes:

Bonus point if you trying to estimate using both approach (one by one and simultaneously) and comparing the results.

1. Sampling Distribusi Laplace dengan Metropolis–Hastings

Tujuan dari bagian ini adalah menghasilkan sampel dari distribusi Laplace dengan parameter \(\mu = 2\) dan \(b = 1\) menggunakan algoritma Metropolis–Hastings, serta mengevaluasi tingkat penerimaan (acceptance rate).

a. Implementasi Metropolis–Hastings

set.seed(67)

# PDF Laplace
laplace_pdf <- function(x, mu, b){
  (1/(2*b)) * exp(-abs(x - mu)/b)
}

# Target density π(x)
target_pi <- function(x){
  laplace_pdf(x, mu = 2, b = 1)
}

# Fungsi Metropolis-Hastings
metropolis_laplace <- function(n_iter, x0){
  x <- numeric(n_iter)
  x[1] <- x0
  accept <- 0
  
  for(t in 2:n_iter){
    x_star <- rnorm(1, mean = x[t-1], sd = 1)
    alpha <- min(1, target_pi(x_star) / target_pi(x[t-1]))
    if(runif(1) < alpha){
      x[t] <- x_star
      accept <- accept + 1
    } else {
      x[t] <- x[t-1]
    }
  }
  
  list(samples = x, accept_rate = accept/(n_iter-1))
}

# Jalankan algoritma
n_total <- 5000
burn_in <- 2000
mh <- metropolis_laplace(n_total, x0 = 0)

Penjelasan kode:

Kode ini mengimplementasikan algoritma Metropolis–Hastings untuk mengambil sampel dari distribusi Laplace dengan parameter 𝜇= 2 dan 𝑏= 1. Pertama, kita mendefinisikan fungsi PDF Laplace (laplace_pdf) dan target densitas 𝜋(𝑥). Fungsi metropolis_laplace kemudian membuat rantai Markov dengan n_iter iterasi mulai dari nilai awal x0. Pada setiap iterasi, kandidat baru x_star diusulkan dari distribusi Normal dengan mean nilai sebelumnya dan standar deviasi 1. Rasio penerimaan α dihitung sebagai perbandingan target density antara kandidat dan state saat ini, kemudian kandidat diterima atau ditolak secara probabilistik. Fungsi ini mengembalikan seluruh sampel dan tingkat penerimaan (accept_rate).

b dan c. Buang burn-in dan simpan 3000 sampel terakhir

samples_kept <- mh$samples[(burn_in + 1):n_total]

head(samples_kept)

## [1] 1.800716 1.800716 1.800716 1.889906 2.008788 2.008788

length(samples_kept)

## [1] 3000

Penjelasan kode:

Setelah menjalankan Metropolis–Hastings, kita membuang fase transien awal (burn-in) sebanyak 2000 sampel untuk memastikan rantai sudah mencapai distribusi target. Sisanya, sebanyak 3000 sampel, disimpan sebagai samples_kept. Fungsi head(samples_kept) menampilkan beberapa sampel pertama, sementara length(samples_kept) memastikan jumlah sampel yang disimpan sesuai dengan yang diharapkan. Dengan cara ini, kita mendapatkan sampel yang merepresentasikan distribusi Laplace(2,1) lebih akurat.

d. Hitung acceptance rate

mh$accept_rate

## [1] 0.6921384

Penjelasan kode:

Acceptance rate dihitung dari proporsi kandidat yang diterima sepanjang iterasi, yaitu jumlah penerimaan dibagi total iterasi minus 1. Dalam contoh ini, mh$accept_rate menghasilkan sekitar 0.692, artinya 69.2% dari kandidat baru diterima.

Interpretasi:

• Nilai awal x0 = 0 menyebabkan rantai memulai jauh dari pusat distribusi Laplace (𝜇= 2), sehingga perlu burn-in untuk mencapai daerah kepadatan tinggi.

• Sampel setelah burn-in (samples_kept) berada di sekitar 𝜇= 2, seperti terlihat dari beberapa nilai pertama (1.8, 1.88, 2.00, dll.), menunjukkan rantai berhasil mengekplorasi daerah utama distribusi.

• Jumlah sampel adalah 3000, sesuai perintah soal, memberikan basis yang cukup untuk analisis selanjutnya.

• Acceptance rate ~0.692 menunjukkan algoritma bekerja efisien, di mana tidak terlalu sering menolak (yang bisa membuat rantai terjebak) dan tidak terlalu sering menerima (yang bisa membuat sampel sangat berkorelasi).

2. Analisis Sampel MCMC dan Distribusi Teoretis

a. Histogram MCMC dan PDF Laplace teoretis

library(ggplot2)

laplace_pdf_theo <- function(x, mu_val = 2, scale = 1){
  1/(2*scale) * exp(-abs(x - mu_val)/scale)
}

df_samples <- data.frame(value = samples_kept)

ggplot(df_samples, aes(x = value)) +
  geom_histogram(aes(y = ..density..), bins = 40, fill = "skyblue", color = "grey30") +
  stat_function(fun = laplace_pdf_theo, args = list(mu_val = 2, scale = 1),
                color = "red", size = 1) +
  labs(title = "Histogram Sampel MCMC dengan PDF Laplace(2,1)",
       x = "x", y = "Density") +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Penjelasan kode:

Kode ini membuat histogram dari 3000 sampel MCMC yang diperoleh dari Metropolis–Hastings. Histogram dinormalisasi sehingga sumbu-y mewakili kepadatan probabilitas (density). Fungsi stat_function menambahkan kurva PDF Laplace teoretis dengan 𝜇 = 2 dan 𝑏 = 1di atas histogram. Dengan visualisasi ini, kita bisa membandingkan apakah sampel MCMC meniru bentuk distribusi teoretis.

b. Empirical CDF vs CDF teoretis

laplace_cdf_theo <- function(x, mu_val = 2, scale = 1){
  sapply(x, function(xx){
    if(xx < mu_val) 0.5 * exp((xx - mu_val)/scale) 
    else 1 - 0.5 * exp(-(xx - mu_val)/scale)
  })
}

ecdf_mcmc <- ecdf(samples_kept)
df_cdf <- data.frame(x = sort(unique(samples_kept)))
df_cdf$empirical <- ecdf_mcmc(df_cdf$x)
df_cdf$theoretical <- laplace_cdf_theo(df_cdf$x)

ggplot(df_cdf, aes(x = x)) +
  geom_step(aes(y = empirical), color = "black", direction = "hv", size = 0.8) +
  geom_line(aes(y = theoretical), color = "red", linetype = "dashed", size = 0.9) +
  labs(title = "Empirical CDF vs Theoretical CDF Laplace(2,1)",
       x = "x", y = "CDF") +
  theme_minimal()

Penjelasan kode:

Kode ini membandingkan CDF empiris dari sampel MCMC dengan CDF teoretis Laplace(2,1). Fungsi ecdf(samples_kept) menghasilkan CDF empiris dari data sampel, kemudian laplace_cdf_theo mendefinisikan CDF teoretis Laplace. Kedua CDF digambarkan pada plot yang sama: garis hitam untuk CDF empiris dan garis merah putus-putus untuk CDF teoretis.

c. Hitung mean dan variance sampel

sample_mean_mcmc <- mean(samples_kept)
sample_var_mcmc  <- var(samples_kept)

sample_mean_mcmc

## [1] 1.991308

sample_var_mcmc

## [1] 2.077282

Penjelasan kode:

Kode ini menghitung statistik dasar dari sampel MCMC: sample_mean_mcmc untuk rata-rata dan sample_var_mcmc untuk varians. Hasilnya adalah mean ≈ 1.991 dan varians ≈ 2.077. Nilai-nilai ini digunakan untuk menilai apakah sampel MCMC secara numerik mendekati karakteristik distribusi Laplace(2,1), yaitu mean teoretis = 2 dan varians teoretis = 2.

d. Bandingkan dengan mean dan variance teoretis Laplace(2,1)

theo_mean <- 2
theo_var  <- 2 * 1^2

diff_mean <- sample_mean_mcmc - theo_mean
diff_var  <- sample_var_mcmc  - theo_var

diff_mean

## [1] -0.008691808

diff_var

## [1] 0.07728249

Penjelasan kode:

Di sini kita membandingkan statistik sampel dengan nilai teoretis. Selisih mean ≈ -0.0087 dan selisih varians ≈ 0.077 menunjukkan bahwa rata-rata dan variansi sampel sangat dekat dengan nilai teoretis. Perbedaan kecil ini wajar dan muncul karena keterbatasan jumlah sampel, variasi Monte Carlo, dan autokorelasi sampel dalam rantai MCMC.

Interpretasi:

• Histogram & PDF: Histogram sampel mengikuti bentuk distribusi Laplace dengan puncak di sekitar 𝜇 = 2 dan ekor yang tebal, menunjukkan algoritma MCMC berhasil mengekstrak pola utama distribusi target.

• CDF empiris hampir identik dengan CDF teoretis di seluruh domain, menandakan bahwa sampel menyebar sesuai probabilitas teoretis.

• Nilai mean ≈ 1.991 dan varians ≈ 2.077 mendekati teoretis (2 dan 2), yang mengonfirmasi kualitas sampel.

Secara visual dan numerik, sampel MCMC berhasil meniru distribusi Laplace(2,1) dengan baik. Perbedaan kecil hanyalah efek statistik normal dari sampel terbatas.

3. Estimasi Parameter Distribusi Laplace

a dan b. Estimasi μ dan b menggunakan metode Newton–Raphson

# Estimasi μ dengan Newton-Raphson
estimate_mu <- function(data, mu_init = median(data), tol = 1e-8, eps = 1e-6, max_iter = 100){
mu <- mu_init
for(i in 1:max_iter){
grad <- sum((data - mu)/sqrt((data - mu)^2 + eps))          # turunan pertama
hess <- -sum(eps/((data - mu)^2 + eps)^(3/2))               # turunan kedua
mu_new <- mu - grad/hess
if(abs(mu_new - mu) < tol) return(mu_new)
mu <- mu_new
}
return(mu)
}

# Hitung estimasi μ
mu_hat <- estimate_mu(samples_kept)
mu_hat

## [1] 2.002956

# Estimasi b (skala) menggunakan nilai μ hasil estimasi
b_hat <- mean(abs(samples_kept - mu_hat))
b_hat

## [1] 1.011045

Penjelasan kode:

Kode ini mengestimasi parameter lokasi 𝜇 menggunakan metode Newton–Raphson. Fungsi estimate_mu mengiterasi hingga perubahan μ kecil (tol) atau jumlah iterasi maksimum tercapai. Gradien dan Hessian dihitung dengan penyesuaian eps untuk menghindari pembagian dengan nol. Setelah μ diperoleh, skala 𝑏 diestimasi sebagai rata-rata jarak absolut sampel terhadap μ yang sudah diperkirakan.

c. Bandingkan estimasi dengan nilai sebenarnya

# Nilai teoretis
true_mu <- 2
true_b  <- 1

# Tampilkan hasil perbandingan
cat("True μ:", true_mu, "\n")

## True μ: 2

cat("Estimated μ:", mu_hat, "\n\n")

## Estimated μ: 2.002956

cat("True b:", true_b, "\n")

## True b: 1

cat("Estimated b:", b_hat, "\n")

## Estimated b: 1.011045

Penjelasan kode:

Kode ini menampilkan hasil estimasi μ dan b dari data sampel MCMC, sekaligus membandingkannya dengan nilai teoretis (𝜇 = 2, 𝑏 = 1). Hal ini memungkinkan kita mengevaluasi seberapa baik prosedur estimasi menangkap karakteristik distribusi asli.

Interpretasi:

• μ̂ ≈ 2: Estimasi lokasi sangat dekat dengan nilai teoretis, menunjukkan konsistensi estimator Newton–Raphson untuk parameter lokasi.

• b̂ ≈ 1: Estimasi skala juga mendekati nilai teoretis, menandakan bahwa variabilitas sampel MCMC berhasil menangkap skala distribusi Laplace.

• Selisih kecil: Perbedaan minor mungkin muncul karena keterbatasan jumlah sampel, autokorelasi dalam rantai MCMC, atau efek Monte Carlo.

Kesimpulannya, metode Newton–Raphson berhasil mengekstrak parameter μ dan b dari sampel MCMC secara akurat dan konsisten.

---

Nomor 2

Sphinx, Siren, dan Stellar adalah tiga tim penelitian independen yang beroperasi secara mandiri, namun saat ini tengah terlibat dalam kompetisi nasional yang sangat ketat, yakni masing-masing tim bersaing untuk menghasilkan sampel eksperimen unggulan yang akan digunakan sebagai bahan penggalangan dana dari sponsor kompetisi. Keberhasilan setiap tim dalam memproduksi sampel tidak hanya diukur dari jumlah unit yang berhasil dibuat, tetapi juga dari tingkat ketelitian prosedur laboratorium, efisiensi penggunaan sumber daya yang tersedia, kepatuhan terhadap protokol keselamatan laboratorium, kemampuan tim dalam meminimalkan risiko degradasi bahan, serta kecermatan dalam menjaga kualitas alat dan sampel agar tetap stabil sepanjang periode eksperimen. Mengingat kompleksitas prosedur yang tinggi, variasi hasil yang mungkin muncul selama eksperimen, serta tuntutan sponsor yang ketat untuk menghasilkan data yang konsisten dan dapat direproduksi, setiap tim dituntut untuk menyusun strategi produksi harian yang maksimal secara finansial, sambil tetap menyeimbangkan risiko eksperimen, keterbatasan sumber daya, dan kapasitas teknisi yang tersedia.

Nilai Sponsor dan Biaya Produksi

Larutan X merupakan sampel yang relatif sederhana, namun tetap menuntut ketelitian penuh dari seluruh teknisi laboratorium agar prosedur berjalan lancar dan reaksi kimia stabil. Sponsor memberikan nilai sebesar Rp60.000 untuk setiap unit yang berhasil diproduksi, mengingat larutan ini menjadi dasar untuk berbagai analisis lanjutan yang krusial bagi eksperimen tingkat lanjut. Setiap unit Larutan X memerlukan alokasi bahan kimia senilai Rp5.000, tenaga teknisi senilai Rp20.000 per unit, penggunaan alat laboratorium dengan biaya Rp2.500 per jam, serta biaya penyimpanan sebesar Rp5.000 per unit.

Senyawa Y memiliki kompleksitas sintesis yang lebih tinggi dibanding Larutan X dan dihargai sponsor sebesar Rp100.000 per unit. Produksi setiap unit Senyawa Y menuntut alokasi bahan kimia Rp7.000, tenaga teknisi Rp30.000, penggunaan alat laboratorium Rp3.000 per jam, serta biaya penyimpanan Rp8.750 per unit.

Bahan Z adalah sampel paling kompleks dan bernilai tinggi. Sponsor memberikan Rp150.000 untuk setiap unit. Setiap unit membutuhkan bahan kimia Rp10.000, tenaga teknisi Rp40.000, penggunaan alat laboratorium Rp4.000 per jam, serta biaya penyimpanan Rp12.500 per unit.

Keterbatasan Sumber Daya

Setiap tim dihadapkan pada keterbatasan signifikan terkait jumlah bahan kimia yang tersedia, kapasitas operasional alat laboratorium, dan jumlah jam teknisi per hari:

• Sphinx: 250 unit bahan kimia, kapasitas alat 180 jam per hari.

• Siren: 350 unit bahan kimia, kapasitas alat 250 jam per hari.

• Stellar: 300 unit bahan kimia, kapasitas alat 200 jam per hari.

Penggunaan Sumber Daya per Sampel

• Larutan X membutuhkan 3 unit bahan kimia dan 2 jam penggunaan alat per unit.

• Senyawa Y membutuhkan 4 unit bahan kimia dan 4 jam penggunaan alat per unit.

• Bahan Z membutuhkan 5 unit bahan kimia dan 6 jam penggunaan alat per unit.

Aturan Produksi dan Persyaratan Sponsor

Sponsor mensyaratkan produksi minimum harian:

• Sphinx: minimal 20 X, 10 Y, 15 Z.

• Siren: minimal 25 X, 15 Y, 20 Z.

• Stellar: minimal 15 X, 10 Y, 10 Z.

Batas proporsi tambahan:

• Jumlah X harus minimal dua kali jumlah Y.

• Jumlah Z tidak boleh melebihi jumlah X.

Tujuan Tim Penelitian

Tujuan utama setiap tim adalah memaksimalkan pendanaan bersih sponsor, yaitu total nilai sponsor dikurangi seluruh biaya bahan, tenaga teknisi, operasional alat, dan penyimpanan. Strategi produksi optimal harus mempertimbangkan keterbatasan sumber daya, proporsi minimum sampel, risiko degradasi bahan, serta kualitas data eksperimen.

Pertanyaan

1. Tentukan berapa banyak Larutan X, Senyawa Y, dan Bahan Z yang harus diproduksi setiap harinya oleh Sphinx, Siren, dan Stellar agar masing-masing tim dapat memaksimalkan pendanaan bersih setelah dikurangi semua biaya produksi dan operasional.

2. Bandingkan pendanaan bersih maksimum yang diperoleh ketiga tim dan tentukan tim mana yang memiliki potensi optimal untuk meningkatkan pendanaan bersih dengan memanfaatkan keterbatasan bahan dan jam laboratorium secara efisien.

3. Jika tim diberikan kesempatan untuk menyesuaikan strategi produksi tanpa menambah total bahan kimia atau jam laboratorium, mereka memiliki dua opsi:

• Opsi 1: Menaikkan jumlah minimal Larutan X sebesar 15%.

• Opsi 2: Menurunkan batas produksi Bahan Z menjadi 80% dari Larutan X.

Pilih salah satu opsi dan jelaskan secara logis (tanpa perhitungan numerik) dampaknya terhadap:

• pendanaan bersih tim,

• risiko eksperimen, dan

• kemampuan memenuhi persyaratan sponsor.

Argumen harus mempertimbangkan interaksi antara keterbatasan sumber daya, kompleksitas sampel, serta potensi risiko degradasi atau kerusakan bahan.

library(lpSolve) # untuk pemecahan integer linear programming
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(ggplot2)

# parameter per-produk (all currency in IDR)
products <- tibble(
name = c("X","Y","Z"),
sponsor = c(60000,100000,150000),
cost_material = c(5000,7000,10000), # biaya moneter bahan per unit
cost_tech = c(20000,30000,40000),
hours_per_unit = c(2,4,6),
cost_tool_per_hour = c(2500,3000,4000),
cost_storage = c(5000,8750,12500),
chem_units = c(3,4,5)
) %>%
mutate(cost_tool_per_unit = hours_per_unit * cost_tool_per_hour,
profit_per_unit = sponsor - (cost_material + cost_tech + cost_tool_per_unit + cost_storage))

products

## # A tibble: 3 × 10
##   name  sponsor cost_material cost_tech hours_per_unit cost_tool_per_hour
##   <chr>   <dbl>         <dbl>     <dbl>          <dbl>              <dbl>
## 1 X       60000          5000     20000              2               2500
## 2 Y      100000          7000     30000              4               3000
## 3 Z      150000         10000     40000              6               4000
## # ℹ 4 more variables: cost_storage <dbl>, chem_units <dbl>,
## #   cost_tool_per_unit <dbl>, profit_per_unit <dbl>

Interpretasi:

Pada tabel products:

• profit_per_unit menunjukkan keuntungan bersih tiap unit setelah dikurangi semua biaya.

  • X: 30.000 IDR, Y: 45.250 IDR, Z: 68.500 IDR.

  • Z paling menguntungkan per unit, walaupun memerlukan bahan dan jam lebih banyak.

• cost_tool_per_unit adalah biaya alat per unit sesuai jam pemakaian. Misal X membutuhkan 2 jam × 2.500 = 5.000 IDR/unit.

Angka-angka ini menjadi dasar perhitungan profit dalam optimasi tim.

Penjelasan untuk tabel products:

Tabel products menampilkan struktur biaya dan nilai sponsor untuk setiap jenis sampel (X, Y, dan Z) yang diproduksi oleh tim. Data tersebut memperlihatkan bagaimana setiap unit produk memiliki komposisi biaya yang berbeda, meliputi kebutuhan bahan kimia, jam teknisi, jam penggunaan alat, serta biaya penyimpanan. Setelah seluruh komponen biaya dijumlahkan, terlihat bahwa margin keuntungan per unit menjadi semakin besar seiring meningkatnya kompleksitas produk, di mana sampel Z memiliki margin tertinggi meskipun biaya produksinya paling besar. Informasi ini penting karena struktur biaya–nilai inilah yang menentukan strategi optimasi tiap tim ketika dihadapkan pada keterbatasan bahan kimia dan jam alat.

# Fungsi penyelesaian untuk satu tim
solve_team <- function(team_name, chem_capacity, hour_capacity, mins, prop_X2Y = TRUE, prop_Z_le_X = TRUE){

# koefisien objektif (maximasi profit total)
obj <- products$profit_per_unit

# Matriks kendala (baris = constraints, kolom = X,Y,Z)
constr_mat <- rbind(
# bahan kimia
products$chem_units,
# jam alat
products$hours_per_unit
)
constr_dir <- c("<=","<=")
constr_rhs <- c(chem_capacity, hour_capacity)

# minimum produksi -> kita ubah menjadi constraint: n_i >= min_i -> -n_i <= -min_i for lpSolve
constr_mat <- rbind(constr_mat, -diag(3))
constr_dir <- c(constr_dir, rep("<=",3))
constr_rhs <- c(constr_rhs, -c(mins["X"], mins["Y"], mins["Z"]))

# proporsional: X >= 2Y -> X - 2Y >= 0 -> -X + 2Y <= 0 (ubah ke <=)
if(prop_X2Y){
constr_mat <- rbind(constr_mat, c(-1,2,0))
constr_dir <- c(constr_dir, "<=")
constr_rhs <- c(constr_rhs, 0)
}
# proporsional: Z <= X -> Z - X <= 0
if(prop_Z_le_X){
constr_mat <- rbind(constr_mat, c(-1,0,1)) # -X + 0Y + Z <= 0
constr_dir <- c(constr_dir, "<=")
constr_rhs <- c(constr_rhs, 0)
}

# Semua variabel integer
res <- lp(direction = "max",
objective.in = obj,
const.mat = constr_mat,
const.dir = constr_dir,
const.rhs = constr_rhs,
all.int = TRUE)

if(res$status != 0) {
warning(paste("Solver status for", team_name, "is", res$status, "(0 means optimal)."))
}

x <- res$solution
names(x) <- products$name
total_profit <- sum(x * obj)

# ringkasan pemakaian resource dan pendapatan kotor
usage_chem <- sum(x * products$chem_units)
usage_hours <- sum(x * products$hours_per_unit)
revenue <- sum(x * products$sponsor)
total_costs <- sum(x * (products$cost_material + products$cost_tech + products$cost_tool_per_unit + products$cost_storage))

tibble(
team = team_name,
X = x["X"], Y = x["Y"], Z = x["Z"],
revenue = revenue,
total_costs = total_costs,
net_profit = total_profit,
chem_used = usage_chem,
chem_cap = chem_capacity,
hours_used = usage_hours,
hours_cap = hour_capacity
)
}

Interpretasi:

• X, Y, Z adalah jumlah unit masing-masing produk yang optimal diproduksi.

  • Sphinx memilih X=25, Y=10, Z=15. Jumlah X melebihi minimum (20), Z sama dengan minimum (15).

  • Siren memproduksi lebih banyak X (35) karena kapasitas bahan/jam alat lebih besar.

• revenue adalah total pendapatan dari sponsor. Misal Sphinx: (25×60.000 + 10×100.000 + 15×150.000 = 4.750.000).

• total_costs adalah jumlah biaya bahan, teknisi, alat, dan penyimpanan.

• net_profit adalah selisih revenue – total_costs. Menunjukkan keuntungan bersih tim jika mengikuti produksi optimal.

• chem_used dan hours_used adalah pemakaian sumber daya aktual. Misal Sphinx menggunakan 190/250 unit bahan → masih ada sisa 60 unit bahan.

Penjelasan untuk fungsi solve_team:

Fungsi solve_team membangun suatu model optimasi linear yang dirancang untuk menentukan kombinasi produksi X, Y, dan Z yang menghasilkan keuntungan bersih maksimum bagi suatu tim. Dalam fungsi ini dimasukkan seluruh batasan penting yang harus dipatuhi tim, termasuk keterbatasan bahan kimia dan jam alat, batas produksi minimum, serta aturan proporsi yang ditetapkan sponsor. Hasil optimasi kemudian merangkum berapa unit setiap produk yang optimal untuk diproduksi, berapa banyak sumber daya yang digunakan, dan berapa keuntungan bersih yang diperoleh.

Menyelesaikan ketiga tim dengan data yang diberikan

# minima per tim
mins_sphinx <- c(X=20, Y=10, Z=15)
mins_siren <- c(X=25, Y=15, Z=20)
mins_stellar<- c(X=15, Y=10, Z=10)

res_sphinx <- solve_team("Sphinx", chem_capacity = 250, hour_capacity = 180, mins = mins_sphinx)
res_siren <- solve_team("Siren", chem_capacity = 350, hour_capacity = 250, mins = mins_siren)
res_stellar<- solve_team("Stellar",chem_capacity = 300, hour_capacity = 200, mins = mins_stellar)

results <- bind_rows(res_sphinx, res_siren, res_stellar)
results

## # A tibble: 3 × 11
##   team        X     Y     Z revenue total_costs net_profit chem_used chem_cap
##   <chr>   <dbl> <dbl> <dbl>   <dbl>       <dbl>      <dbl>     <dbl>    <dbl>
## 1 Sphinx     25    10    15 4750000     2750000    2000000       190      250
## 2 Siren      35    15    20 6600000     3821250    2778750       265      350
## 3 Stellar    50    10    10 5500000     3192500    2307500       240      300
## # ℹ 2 more variables: hours_used <dbl>, hours_cap <dbl>

Interpretasi:

• rank adalah peringkat tim berdasarkan net_profit. Siren berada di posisi 1 (2.778.750 IDR), Stellar ke-2 (2.307.500 IDR), Sphinx ke-3 (2.000.000 IDR).

• chem_ratio & hour_ratio adalah pemakaian relatif sumber daya. Misal Sphinx: 190/250 = 0,76 (76%) bahan, jam 180/180 = 1 (100%).

• bottleneck adalah sumber daya yang paling membatasi produksi. Misal Sphinx “Jam alat” → walau masih ada sisa bahan, produksi tidak bisa ditambah karena jam teknisi/alatalat sudah penuh.

Penjelasan untuk output hasil optimasi tiga tim:

Hasil optimasi menunjukkan bahwa setiap tim memilih memproduksi jumlah yang sedikit di atas batas minimal, kecuali pada produk X untuk Siren dan Stellar yang meningkat karena masih tersedia ruang dalam kapasitas bahan dan jam alat. Pilihan produksi ini memperlihatkan bagaimana masing-masing tim memaksimalkan margin keuntungan tanpa melanggar batas sumber daya yang diberikan.

Visualisasi dan perbandingan

results_long <- results %>%
select(team, X, Y, Z, net_profit, chem_used, hours_used) %>%
pivot_longer(cols = c(X,Y,Z), names_to = "product", values_to = "units")

# plot produksi per tim
ggplot(results_long, aes(x = team, y = units, fill = product)) +
geom_col(position = "dodge") +
labs(title = "Komposisi Produksi per Tim", y = "Jumlah unit per hari", x = "Tim")

# plot net profit
ggplot(results, aes(x = team, y = net_profit)) +
geom_col() +
labs(title = "Pendanaan Bersih (Net Profit) per Tim", y = "Pendanaan Bersih (IDR)", x = "Tim")

results %>% mutate(rank = rank(-net_profit)) %>% arrange(rank)

## # A tibble: 3 × 12
##   team        X     Y     Z revenue total_costs net_profit chem_used chem_cap
##   <chr>   <dbl> <dbl> <dbl>   <dbl>       <dbl>      <dbl>     <dbl>    <dbl>
## 1 Siren      35    15    20 6600000     3821250    2778750       265      350
## 2 Stellar    50    10    10 5500000     3192500    2307500       240      300
## 3 Sphinx     25    10    15 4750000     2750000    2000000       190      250
## # ℹ 3 more variables: hours_used <dbl>, hours_cap <dbl>, rank <dbl>

# Menunjukkan bottleneck (apakah jam penuh atau bahan kimia penuh?)
res_bottleneck <- results %>%
mutate(chem_ratio = chem_used / chem_cap,
hour_ratio = hours_used / hours_cap,
bottleneck = ifelse(hour_ratio >= chem_ratio, "Jam alat", "Bahan kimia"))
res_bottleneck

## # A tibble: 3 × 14
##   team        X     Y     Z revenue total_costs net_profit chem_used chem_cap
##   <chr>   <dbl> <dbl> <dbl>   <dbl>       <dbl>      <dbl>     <dbl>    <dbl>
## 1 Sphinx     25    10    15 4750000     2750000    2000000       190      250
## 2 Siren      35    15    20 6600000     3821250    2778750       265      350
## 3 Stellar    50    10    10 5500000     3192500    2307500       240      300
## # ℹ 5 more variables: hours_used <dbl>, hours_cap <dbl>, chem_ratio <dbl>,
## #   hour_ratio <dbl>, bottleneck <chr>

Penjelasan untuk peringkat profit dan analisis bottleneck:

Ketika hasil ketiga tim diberi peringkat berdasarkan pendanaan bersih, terlihat bahwa Siren menjadi tim paling efisien dalam memanfaatkan sumber dayanya, sementara Sphinx berada pada posisi terbawah karena kapasitasnya paling ketat. Analisis bottleneck mengungkap bahwa jam penggunaan alat merupakan pembatas utama bagi semua tim, bukan bahan kimia. Artinya, meskipun masih terdapat sisa bahan kimia, tim tidak bisa memperluas produksi karena keterbatasan waktu operasional alat laboratorium. Temuan ini mengisyaratkan bahwa investasi tambahan dalam kapasitas alat akan secara langsung meningkatkan fleksibilitas produksi dan potensi keuntungan tim.

Skenario: Opsi 1 (naikkan minimal X 15%) dan Opsi 2 (turunkan batas Z jadi 80% dari X)

# Opsi 1: X_min <- ceiling(1.15 * X_min_original)
mins_sphinx_o1 <- mins_sphinx; mins_sphinx_o1["X"] <- ceiling(1.15 * mins_sphinx["X"])
mins_siren_o1 <- mins_siren ; mins_siren_o1["X"] <- ceiling(1.15 * mins_siren["X"])
mins_stellar_o1<- mins_stellar; mins_stellar_o1["X"]<- ceiling(1.15 * mins_stellar["X"])

res_sphinx_o1 <- solve_team("Sphinx_Opsi1", chem_capacity = 250, hour_capacity = 180, mins = mins_sphinx_o1)
res_siren_o1 <- solve_team("Siren_Opsi1", chem_capacity = 350, hour_capacity = 250, mins = mins_siren_o1)
res_stellar_o1<- solve_team("Stellar_Opsi1",chem_capacity = 300, hour_capacity = 200, mins = mins_stellar_o1)

solve_team_op2 <- function(team_name, chem_capacity, hour_capacity, mins){
best <- NULL

for(Xtry in 0:500){
  obj_sub <- products$profit_per_unit[c(2,3)] + products$profit_per_unit[1] * Xtry*0 
  const_profit_X <- Xtry * products$profit_per_unit[1]
  
  # constraints matrix untuk [Y,Z]
  # chemical: 4*Y + 5*Z <= chem_capacity - 3*Xtry
  # hours: 4*Y + 6*Z <= hour_capacity - 2*Xtry
  # mins: Y >= minY, Z >= minZ -> -Y <= -minY etc
  chem_rem <- chem_capacity - products$chem_units[1] * Xtry
  hour_rem <- hour_capacity - products$hours_per_unit[1] * Xtry
  if(chem_rem < 0 || hour_rem < 0) next
  
  constr_mat_sub <- rbind(
    c(4,5),
    c(4,6),
    c(-1,0),
    c(0,-1),
    c(0,1) # will replace by Z <= floor(0.8*Xtry) as later
    )
  
  constr_dir_sub <- c("<=","<=","<=","<=", "<=")
  constr_rhs_sub <- c(chem_rem, hour_rem, -mins["Y"], -mins["Z"], floor(0.8 * Xtry))
  
  # last row currently is placeholder: we want Z <= floor(0.8*Xtry) -> 0*Y + 1*Z <= floor(0.8*Xtry)
  constr_mat_sub[5,] <- c(0,1)
  constr_dir_sub[5] <- "<="
  constr_rhs_sub[5] <- floor(0.8 * Xtry)
  
  # add prop X >= 2Y -> with X fixed: -1*Xtry + 2*Y <= 0 -> 2*Y <= Xtry -> Y <= floor(Xtry/2)
  constr_mat_sub <- rbind(constr_mat_sub, c(2,0))
  constr_dir_sub <- c(constr_dir_sub, "<=")
  constr_rhs_sub <- c(constr_rhs_sub, floor(Xtry/2))
  
  # also Z <= Xtry already implied by 0.8*X
  # solve subproblem
  res_sub <- lp(direction = "max",
                objective.in = products$profit_per_unit[c(2,3)],
                const.mat = constr_mat_sub,
                const.dir = constr_dir_sub,
                const.rhs = constr_rhs_sub,
                all.int = TRUE)
  
  if(res_sub$status == 0){
    y <- res_sub$solution[1]
    z <- res_sub$solution[2]
    total_profit <- const_profit_X + sum(res_sub$solution * products$profit_per_unit[c(2,3)]) 
    + Xtry * products$profit_per_unit[1]
    
    if(is.null(best) || total_profit > best$total_profit){
      best <- list(X = Xtry, Y = y, Z = z, total_profit = total_profit)
      }
    }
  }
if(is.null(best)) return(tibble(team = team_name, X=NA,Y=NA,Z=NA, net_profit = NA))
tibble(team = team_name, X = best$X, Y = best$Y, Z = best$Z, net_profit = best$total_profit)
}

res_sphinx_o2 <- solve_team_op2("Sphinx_Opsi2", chem_capacity = 250, hour_capacity = 180, mins = mins_sphinx)
res_siren_o2 <- solve_team_op2("Siren_Opsi2", chem_capacity = 350, hour_capacity = 250, mins = mins_siren)
res_stellar_o2 <- solve_team_op2("Stellar_Opsi2", chem_capacity = 300, hour_capacity = 200, mins = mins_stellar)

bind_rows(res_sphinx_o1, res_siren_o1, res_stellar_o1) %>% rename(type = team) -> r_op1
bind_rows(res_sphinx_o2, res_siren_o2, res_stellar_o2) %>% rename(type = team) -> r_op2

r_op1

## # A tibble: 3 × 11
##   type           X     Y     Z revenue total_costs net_profit chem_used chem_cap
##   <chr>      <dbl> <dbl> <dbl>   <dbl>       <dbl>      <dbl>     <dbl>    <dbl>
## 1 Sphinx_Op…    25    10    15 4750000     2750000    2000000       190      250
## 2 Siren_Ops…    35    15    20 6600000     3821250    2778750       265      350
## 3 Stellar_O…    50    10    10 5500000     3192500    2307500       240      300
## # ℹ 2 more variables: hours_used <dbl>, hours_cap <dbl>

r_op2

## # A tibble: 3 × 5
##   type              X     Y     Z net_profit
##   <chr>         <int> <dbl> <dbl>      <dbl>
## 1 Sphinx_Opsi2     NA    NA    NA         NA
## 2 Siren_Opsi2      NA    NA    NA         NA
## 3 Stellar_Opsi2    50    10    10    2307500

Interpretasi:

• Opsi 1: semua tim tetap memproduksi jumlah unit sama dengan skenario awal, sehingga profit (net_profit) tidak berubah. Angka X, Y, Z sama seperti sebelumnya.

• Opsi 2: Sphinx dan Siren menunjukkan NA → tidak ada solusi feasible dengan batas Z ≤ 0,8 X. Stellar tetap bisa memproduksi 50X, 10Y, 10Z dengan net_profit 2.307.500.

  • Angka NA menunjukkan kombinasi minimum + proporsi yang lebih ketat tidak bisa dipenuhi oleh tim karena keterbatasan sumber daya.

Penjelasan untuk Hasil Opsi 1 dan Opsi 2:

Hasil skenario Opsi 1 memperlihatkan bahwa menaikkan batas minimum produksi X sebesar 15% tidak mengubah solusi optimal pada ketiga tim. Hal ini terjadi karena ketiga tim memang sudah memilih memproduksi X di atas batas minimumnya, sehingga pengetatan ini tidak cukup untuk memaksa perubahan keputusan. Keuntungan bersih, penggunaan sumber daya, dan pola produksi semua tim tetap identik dengan skenario awal. Dengan demikian, Opsi 1 tidak memberikan dampak nyata dalam konteks struktur biaya dan kapasitas sumber daya yang tersedia.

Sebaliknya, hasil Opsi 2 memperlihatkan bahwa beberapa tim tidak mampu menemukan solusi feasible ketika batas produksi Z dipersempit menjadi maksimal 80% dari X. Pada Sphinx dan Siren, seluruh kombinasi produksi yang mungkin tidak dapat memenuhi secara simultan: batas minimal Z yang ditetapkan sponsor, proporsi Z ≤ 0.8X, dan keterbatasan jam alat. Hanya Stellar yang tetap menemukan solusi optimal—yang kebetulan sama dengan hasil sebelumnya karena memang sudah berada dalam batas 80%. Hal ini menunjukkan bahwa Opsi 2 jauh lebih restriktif dan dapat membuat konfigurasi produksi menjadi tidak mungkin dipenuhi, terutama bagi tim dengan kapasitas jam alat yang lebih ketat dan kebutuhan minimum Z yang lebih tinggi.