Daftar Isi Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
TAKSIRAN INTERVAL KEPERCAYAAN PADA PERHITUNGAN k eff MONTE CARLO Entin Hartini*
ABSTRAK TAKSIRAN INTERVAL KEPERCAYAAN PADA PERHITUNGAN keff MONTE CARLO . Hasil Monte Carlo menggambarkan rata-rata kontribusi dari beberapa history sampel sepanjang perjalanan neutron . Satu kuantitas penting dalam Monte Carlo adalah kekeliruan statistik dihubungkan dengan hasil. Pentingnya kekeliruan dan jumlah history karena user tidak hanya memperoleh kualitas dari hasil juga menentukan hasil secara statistik baik, sehingga mencerminkan interval konfidensi benar. Tujuan estimasi ini adalah agar user yang mempelajari MCNP memahami tentang interpretasi dari estimasi rata-rata, interval kepercayaan untuk mencapai nilai yang diinginkan. Makalah ini berisi kalkulasi Monte Carlo untuk interval kepercayaan estimasi k eff pada MCNP. Estimator keff adalah kombinasi least square atau kombinasi maksimum likelihood dengan distribusi multivariate normal. Kombinasi dari tiga estimator keff pada MCNP ditunjukan, secara teori dan empirik dengan study statistik untuk memberikan estimator k eff terbaik. Kata-kata kunci: Monte Carlo, Interval Konfidensi, Kritikalitas, ENDF/B-VI MCNP5
ABSTRACT CONVIDENCE INTERVAL ESTIMATE ON MONTE CARLO keff CALCULATION. Monte Carlo's result present an average of the contribution from many history samp led during the course of neutron transport. One important thing in Monte Carlo is statistical error associated with result. The importance for error and history's amount because user not only get quality of result also determine result statistically good, so reflects confidence interval is right. To the effect this estimation is that user that studies MCNP to understand about interpretation of estimation average, confidence interval for reach desirable point. This therefore contains Monte Carlo's calculation to estimations confidence interval k eff on MCNP. Estimator keff is least square's combine or likelihood's maximum combine with multivariate's distribution normal. Combine of three estimator k eff on MCNP at indication, in theory and empirical with study statistic for given the best estimator k eff. Keywords: Monte Carlo, Confidence Interval, Criticality, ENDF/B-VI MCNP5
PENDAHULUAN Metodologi kritikalitas MCNP dan beberapa dasar statistik termasuk Interval kepercayaan sangat diperlukan dalam presentasi dari hasil monte carlo. Kombinasi dari tiga estimator keff pada MCNP ditunjukkan, secara teori dan empirik dengan *
Pusat Pengembangan Informatika Nuklir - BATAN
305
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
studi statistik untuk estimator keff terbaik. Metoda kombinasi estimator didasarkan pada teori Gaus markov, untuk metoda least squares . MCNP tidak menghasilkan estimasi titik dari keff, melainkan memberikan nilai rentang, dengan beberapa konfidensi yang ditetapkan sehingga akan mendapatkan harga keff yang benar. Untuk pertambahan peluang rentang keff , history yang dibutuhkan harus ditambah. Rentang yang merupakan interval kepercayaan sejalan dengan hasil Monte Carlo. Pada makalah ini akan dibahas bagaimana membentuk interval kepercayaan dari tiga estimasi keff terpisah: tumbukan, serapan dan panjang lintasan serta estimasi untuk kombinasi tiga estimator. Estimator keff dan estimator kombinasinya dihitung dengan kombinasi least square atau kombinasi maksimum likelihood dengan distribusi normal variabel ganda untuk kovarian matriks diketahui. Contoh kasus adalah menghitung interval kepercayaan untuk eksperimen kritikalitas benchmark bola plutonium dengan reflektor uranium dan energi kontinyu pada temperatur ruang. Jumlah total partikel yang disimulasikan 1500 neutron persiklus dengan total siklus 200 melewati 25 siklus yang pertama untuk menghindari problem konvergensi sumber.
METODOLOGI Teori Central Limit dan Interval Kepercayaan Monte Carlo Untuk menentukan presisi interval kepercayaan dari Monte Carlo, sentral limit teori dan teori peluang ditetapkan sebagai: [1]
σ σ < x < E(x) + β lim Pr E(x) + α = N N N →∞
1
β
∫e 2π α
−t 2 / 2
dt
(1)
di mana α dan β nilai yang berubah-ubah dan Pr(Z) peluang dari Z . Estimasi deviasi standar dari x adalah S x . Untuk N besar ditulis sebagai:
x − E ( x) < βS x = Pr αS x < σ/ N
1 2π
β
−t ∫e
2
/2
dt
(2)
α
306
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
Interval Kepercayaan Pada MCNP MCNP Secara umum adalah code Monte Carlo untuk transport netron dan radiasi untuk aplikasi kritikalitas nuklir. Estimasi keff menggunakan estimasi standar deviasi untuk membangun interval konfidensi keff. Bilamana dari n sample estimasi x dari variabel acak dengan simpangan baku s dengan interval konfidensi berada dalam jangkauan x ± s , menurut teori central limit distribusi dari estimasi mean mendekati distribusi normal dengan n mendekati ∞. Untuk jumlah sampel tidak terbatas, distribusi didekati oleh distribusi student t, simetrik sekitar nol dan mendekati distribusi normal dengan n →∞. Ini digunakan untuk mendeskripsikan variable random t, dimana [2]
t=
t=
x−µ s (x ) x−µ
s( x ) / n
(3) (4)
Distribusi student t mempunyai distribusi yang berbeda untuk masing-masing n, ditulis sebagai t n-1, dimana n-1 adalah jumlah derajat kebebasan, yaitu jumlah pengukuran independent yang mungkin. Titik dari absis pada graf dari distribusi student t adalah “ persentil dari distribusi” dan ditulis sebagai tn-1, 1-α/2 , dimana indek bawah kedua adalah tingkat kepercayaan.
Estimator keff Pada MCNP Variable bebas x, dengan distribusi dan mean µ tidak diketahui, dengan n sample acak independent (x1, x2, …,xn) dipilih dari distribusi peluang densitas x. Sebuah estimator, X (x1, x2, …,xn) adalah fungsi dari sample acak mewakili secara statistik dengan rata-rata tidak diketahui. Dalam Monte Carlo estimator adalah ratarata dari sample acak. Kualitas yang diinginkan dari estimator adalah : tidak bias, konsisten dan efisien.Estimator x dikatakan tidak bias, jika nilai ekspektasi sama dengan nilai µ dimana E(x) = µ untuk semua µ. Konsep tidak bias telah mendasari pengulangan eksperimen, konsistensi yang lain diterapkan untuk eksperimen tunggal dimana ukuran sample n menjadi besar. Estimator x adalah konsisten jika pendekatan, probabilistik nilai µ dengan n besar. Pada MCNP diasumsikan estimator keff tidak bias, bila varian dari estimasi berkurang 1/n , dimana n jumlah siklus maka kombinasi tiga estimator akan efisien, dimana mempunyai varians terkecil diantara estimator linier. MCNP mempunyai tiga estimator berbeda digunakan untuk mengestimasi keff, kritikalitas dari multy system, tumbukan , serapan dan panjang lintasan estimator. Estimator serapan mempunyai dua bentuk analog serapan atau implicit serapan Ekspresi matematik dari masing-masing estimator pada MCNP dan masingmasing estimasi keff pada siklus khusus. Rata-rata semua siklus memberikan estimasi 307
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
ahir dari keff untuk masing- masing dari tiga estimator,. Indeks atas C,AA,IA dan TL indikasi tumbukan, analog serapan , implicit serapan dan panjang lintasan.[2] Untuk tumbukan, [2] C = keff
∑ f kν kσ Fk 1 Wi k ∑ N i ∑ k f kσ Tk
(5)
di mana
i k
= = =
tumbukan dimana fisi mungkin terjadi isotop yang dikandung pada tumbukan ke i total penampang lintang mikroskopik untuk nuklida k
σ Tk σ Fk = νk =
penampang lintang mikroskopik fisi untuk nuklida k Jumlah rata-rata dari total netron yang dihasilkan perfisi dengan tumbukan isotop pada energi tertentu fraksi atom untuk nuklid k ukuran nominal sumber untuk siklus bobot partikel masuk dalam tumbukan
= = =
fk N Wi
Untuk analog serapan [2]
k effAA =
1 N
∑ν ∑ σ
∑Wi i
k
k
k
Ak
σ Fk
+ σ Fk
(6)
di mana
σ Ak
= penampang lintang mikroskopik serapan, tidak termasuk fisi Untuk Implicit serapan, [2]
k effIA =
1 N
σ Ak + σ Fk σ Tk
∑Wi i
ν k σ Fk σ Ak + σ Fk
(7)
di mana i = dijumlah untuk semua tumbukan dimana terjadi fisi
σ Ak + σ Fk = σ Tk
frekuensi dari analog capture pada masing masing tumbukan, kualitas dengan pembobotan adalah disesuaikan pada masing-masing tumbukan.
308
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
Untuk panjang lintasan, [2] TL k eff =
1 N
∑Wiρd ∑ f υ σ k
i
k
Fk
(8)
k
di mana d = jarak lintasan neutron i = lintasan semua neutron ρ = densitas atom dalam sel
Kombinasi Estimator Distribusi Multivariat yang diperoleh dari sampel multivariat X :{x1,1, ...,x1,n; x2,1,...,x2,n; ...,xk,n,} dimana n adalah jumlah siklus aktif keff dan k adalah jumlah dari estimator, dimana dalam MCNP k= 3. rata-rata n siklus aktif menghasilkan vektor dari k estimator, X = ( x1,..., x k ), dimana masing-masing elemen varians yang berbeda dan berkorelasi satu sama lain. Nilai ekpektasi dari masing-masing x adalah nilai dari keff yaitu µ , varian dan kovarian membentuk kovarian matriks Σ . maka sampel variabel ganda diasumsikan berdistribusi normal dengan vektor rata-rata µ e dan varian Σ, X ~ N( µ e, Σ) dimana vektor e adalah vektor satu dan X ~ N( µ e, Σ) berdistribuasi normal dengan rata-rata µ e dan varian Σ. Penyelesaian Persamaan maximum likelihood variabel ganda memberikan µ , σ dan estimassi invers kovarian matriks dengan meminimumkan jumlah kuadrat kekeliruan.
Kombinasi Dua Estimator Pada MCNP Transformasi dari x pada z untuk 2 estimasi berbeda diberikan dengan: [2]
1 0 x1 x1 = z = Ax = 1 - 1 x 2 x1 − x 2
(9)
Regresi untuk z1 pada z2, yang mana nilai ekspektasinya nol.
309
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
Untuk dua kasus estimator 2 1 0 σ 11 2 ∑ x = 1 - 1 σ 11
σ 122 1 1 σ 222 0 - 1
σ 112 σ 112 - σ 122 = 2 2 2 2 2 σ - σ σ σ 2 σ + − 11 22 12 11 12
(10)
Di mana varian σ ij diestimasi dengan σˆ ij : [3]
σˆ 112 =
s ij2
n -1 1 n σˆ ij2 = ∑ xil − xi )( x jl − x j ) n − 1 l =1 n l 1 n = x x − x ∑ il jl ∑ il ∑ x j l / n n − 1 l =1 l =1 1
(11)
σˆ ij2 adalah kovarian populasi sampel antara n nilai dari xil dan x j l , dimana i=j , σˆ ij2 estimasi varian populasi dari n. Sedangkan sij2 adalah jumlah kuadrat dari deviasi dan rata-rata dibagi dengan derajat kebebasan (n-1) memberikan estimasi tak bias untuk σˆ ij2 . Σ
ˆ berisi σˆ 2 sedangkan matrik S berisi berisi elemen σ ij2 dan ∑ ij
ˆ = S /(n − 1) elemen sij2 .dan jumlah kuadrat memberikan hubungan ∑ Untuk dua estimasi matriks skalar partisi adalah: [3]
S x12 = s112 − s122
(12)
dan [3] 1 2 S x−22 = ( s112 + s22 − 2s122 ) −1
=
1
(13)
2 2 (s11 + s22 − 2s122 )
Untuk dua estimasi persamaan 17 , 20 dan 21 menjadi [2]
310
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
µˆ = x1 − = =
( s112 + s122 )( x1 − x 2 ) 2 − 2 s122 s112 + s 22
2 − σˆ 122 ) x1 + (σˆ 112 − σˆ 122 ) x 2 (σˆ 22 σˆ 112 + σˆ 222 − 2σˆ 122
(14)
(σˆ − σˆ ) (σˆ − σˆ ) + x x2 1 σˆ 112 + σˆ 222 − 2σˆ 122 σˆ 112 + σˆ 222 − 2σˆ 122 2 22
2 12
2 11
2 12
= w 1 x1 + w2 x 2 Varian kombinasi dari dua estimator adalah:
σˆ µ2 =
1 s2 − s2 )2 (( x − x ) 2 1 2 s11 − 2 11 2 12 2 −1 + 2 1 2 2 2 n−2 s11 + s 22 − 2 s12 ) n s11 + s 22 − 2s12 )
(15)
Kombinasi Tiga Estimator Pada MCNP Varian matriks untuk 3tiga estimator berbeda menggunakan transformasi vektor z, maka [3]
1 ∑ x = 1 1
0
0 -1 0 0 - 1
σ 121 σ 122 2 σ 13
σ 122 σ 222 σ
2 23
σ 132 2 σ 23 σ 332
∑ x12 = σ 112 − σ 122 − σ 112 − σ 132
∑ −x122 =
2 2 σ 112 + σ 332 − 2σ 132 − (σ 112 + σ 23 − σ 212 − σ 132 ) 1 2 2 g − (σ 112 + σ 23 − σ 122 − σ 132 ) σ 112 + σ 221 − 2σ 122
(16)
(17)
(18)
di mana :
g = (σ 112 + σ 222 − 2σ 122 )(σ 112 + σ 332 − 2σ 132 ) − (σ 112 + σ 232 − σ 122 − σ 132 ) 2
(19)
rata-rata tiga estimatornya [2]
311
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
µˆ = x1 − −
(σˆ 112 − σˆ 122 ) 2 (σˆ 112 + σˆ 332 − 2σˆ 132 )( x1 − x 2 ) − (σˆ 112 + σˆ 23 − σˆ 122 − σˆ 132 )( x1 − x3 ) g
[
(σˆ 112 − σˆ 132 ) 2 2 − (σˆ 112 + σˆ 231 − σˆ 122 − σˆ 132 )( x1 − x 2 ) + (σˆ 112 + σˆ 22 − 2σˆ 122 )( x1 − x3 ) g
[
]
Bila σˆ ij2 = σˆ 2ji Maka [2] 3
µˆ =
∑fx l =1 3
l
∑f l =1
l
(20)
l
Di mana
f l = σˆ 2jj (σˆ kk2 − σˆ ik2 ) − σˆ kk2 σˆ ij2 + σˆ 2jk (σˆ ij2 + σˆ ik2 − σˆ 2jk )
(21)
di mana l adalah permutasi parsial dari i, j, k seperti berikut ini l 1 2 3
i 1 2 3
j 2 3 1
k 3 1 2
Estimasi varian [2]
σˆ µ2ˆ =
S − 2S 3 ) S1 (1 + n( 2 gn gs
(22)
di mana gs = (n-1)2 g
(23)
3
S1 = ∑ flσˆ 12l l =1 3
S 2 = ∑ ( s 2jj + s kk2 − 2 s 2jk ) xl2
(24)
l =1 3
S 3 = ∑ s kk2 + s ij2 − s 2jk − s ik2 ) xl x j l =1
312
]
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
HASIL DAN PEMBAHASAN KASUS Perhitungan MCNP untuk eksperimen kritikalitas terdiri atas 6,06±0,03 kg teras bola plutonium (densitas=15,53 g/cm3 ) yang dikelilingi lapisan sferis berbatasan konsentris dari uranium dengan ketebalan 7,72 inci dan densitas 19,0 g/cm3 . Eksperimen diasumsikan pada temperatur ruang (2930 K). Radius teras benchmark berbentuk bola Pu =4,5332 cm. Bola direfleksikan oleh uranium dengan tebal 19,6088 cm, sehingga radius luarnya menjadi 24,1420 cm. Densitas atom teras Pu dan reflektor uranium yang dihitung dari komposisi isotopik diberikan pada Tabel 2. Perhitungan MCNP untuk teras bola plutoniumdengan reflector uranium dikerjakan dengan opsi “KCODE” dan ”KSRC” dengan tampang lintang ENDF/B-V energi kontinu pada temperatur ruang. Jumlah total partikel yang disimulasikan adalah 1500 neutron per siklus dengan total siklus 200 melewati 25 siklus.
Gambar 1. Model Geometri MCNP Tabel 1. Komposisi Bahan Isotop/ elemen 239 Pu 240 Pu 241 Pu Ga
Fraksi Atom % 93,80 4,80 0,30 1,10 313
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
Tabel 2. Densitas Atom Isotop Teras 239 Pu 240 Pu 241 Pu Ga Reflektor 234 U 235 U 238 U
Densitas atom 3,6697x10-2 1,8700x10-3 1,1639x10-4 1,4755x10-3 2,6438x10-6 3,4610x10-4 4,7721x10-2
OUTPUT MCNP Estimasi rata-rata keff pada simpangan baku 68, 95, dan 99 % interval kepercayaan pada tabel 3 Tabel 3. Interval Kepercayaan keff Pada Taraf 68%, 95% dan 99% Estimator keff Tumbukan Serapan Panjang lintasan Tumbukan/ serapan Serapan/ pj lintasan Tumbukan/ pj lintasan Tumbukan/ serapan/ pj lintasan
keff
Kepercayaan 68%
Kepercayaan 95%
1,00181 1,00170 1,00138
Simpang an Baku 0,00155 0,00157 0,00142
Kepercayaan 99%
1,00026-1,00337 1,00013-1,00326 0,99996-1,00281
0,99872-1,00491 0,99857-1,00482 0,99855-1,00422
0,99771-1,00592 0,99755-1,00584 0,99762-1,00515
1,00179
0,00156
1,00023-1,00335
0,99869-1,00490
0,9967-1,00591
1,00150
0,00136
1,00014-1,00286
0,99879-1,00420
0,99790-1,00508
1,00154
0,00136
1,00019-1,00290
0,99884-1,00424
0,99796-1,00513
1,00153
0,00136
1,00019-1,00290
0,99882-1,00424
0,99794-1,00512
Estimasi akhir dari kombinasi tumbukan /serapan /panjang lintasan untuk keff = 1.00153 dengan simpangan baku 0.00136 dan estimasi interval kepercayaan keff untuk 68, 95, & 99 persen adalah 1.00017 sampai 1.00289, 0.99882 sampai 1.00424, dan 0.99794 sampai 1.00512
314
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
KESIMPULAN MCNP memberikan interval kepercayaan untuk tiga estimator keff , rata-rata sampel dari masing-masing estimasi dan rata-rata sampel dari gabungan ketiga estimasi . rata-rata terbaik didapat pada kombinasi tiga estimasi dimana mempunyai simpangan baku terkecil.
DAFTAR PUSTAKA 1. RONALD E. WALPOLE AND RAYMOND h. MYERS, “Pobability and Statistics for Engineers and Scie ntists”, Macmillan and Publishing Company, 1989. 2. URBATSCH T.J , “Estimation and Interpretation of keff Confidence Intervals in MCNP”, Los Alamos: LA-12658, 1995. 3. JOHN O. RAWLINGS, “ Applied Regression Analysis”, California , 1988.
Pacific Grove
DISKUSI
GUNARWAN PRAYITNO Hasil estimasi keff terbaik, dasarnya apa ? Dan dimana kontribusi hasil yang anda hitung pada program 5 pilar BATAN. ENTIN HARTINI Dibuat estimasi keff untuk masing-masing dan kombinasinya dari tumbukan, serapan dan panjang lintasan. Estimasi terbaik adalah estimasi dengan varian terkecil. Berhubungan dengan penelitian 2008. JAN SETIAWAN 1.
Taksiran interval kepercayaan merupakan output MCNP5 atau perhitungan terpisah ?
315
Risalah Lokakarya Komputasi dalam Sains dan Teknologi Nuklir: 6-7 Agustus 2008(305-316)
2.
Ketidakpastian pada output MCNP hasil perhitungan atau output interval dari MCNP
ENTIN HARTINI 1. 2.
Modul perhitungan MCNP dikupas dan dibuat modul tersendiri Dibandingkan keff dari MCNP terhadap keff = 1
DAFTAR RIWAYAT HIDUP
1.
Nama
: Dra. Entin Hartini
2.
Tempat/Tanggal Lahir
: Majalengka, 19 Februari 1962
3.
Instansi
: BATAN
4.
Pekerjaan / Jabatan
: Staf Peneliti PPIN - BATAN
5.
Riwayat Pendidikan
: S1 Statistik Universitas Padjadjaran
6.
Pengalaman Kerja
: PPIN - BATAN
Daftar Isi
316