## Kontrol Model Prediksi Robas Pada Waste Heat Boiler : Parametric Ellipsoidal Uncertainty
F. Mukhlish, E. Joelianto, D. Chaerani.
Program Studi Teknik Fisika Institut Teknologi Bandung
Received date : 21 Mei 2010 Accepted date : 29 Mei 2010

## Abstrak
Kontrol Model Prediksi adalah suatu strategi kontrol yang dirancang dengan menggunakan model dari proses suatu sistem. Model proses ini digunakan untuk menghitung himpunan prediksi ke depan dari proses  tersebut  berdasarkan  himpunan  sinyal  kontrol  yang  diberikan  pada  model.  Himpunan sinyal kontrol tersebut diperoleh dari suatu algoritma optimisasi yang meminimalkan fungsi kinerja MPC dan solusinya  memenuhi  kendala  sistem.  Oleh  karena  itu,  pemilihan  model  yang  tepat  ak an mempengaruhi kinerja dari pengontrol.
Optimisasi merupakan bagian terpenting pada perhitungan sinyal control ke depan pada MPC. Pada tugas  akhir  ini  digunakan  optimisasi  robas.  Optimisasi  robas  atau  Robust  Optimization  (RO)  adalah metoda optimisasi yang dapat menangani ketidakpastian data. Ketidakpastian dalam hal ini adalah data dari masalah optimisasi tidak diketahui secara pasti saat menentukan solusi optimal. PT.Petrokimia  Gresik  Pabrik  Pupuk  Nitrogen  memiliki  boiler  bertipe  Waste  Heat  Boiler  (WHB), pembangkit uap yang memanfaatkan gas panas sisa dari generator turbin gas. Sistem ini digunakan untuk  melihat  kinerja  MPC  dengan  optimisasi  robas  dalam  sistem  di  industri  dan  melakukan penghematan bahan bakar  pada sistem WHB.  Mekanisme pengontrolan dilakukan dengan keadaan awal  0  untuk  setiap  masukan  dan  keluaran,  lalu  diamati  kestabilan sistem  dan  variabel  prosesnya terhadap kendala sistem.
RO-MPC  dengan  ketidakpastian  elips  parametrik  dapat  mengkompensasi  gangguan  terukur  dari generator  turbin  gas.  Nilai  RMSE  dan  IAE  dari  keluaran  yang  mempengaruhi  kualitas steam  (uap), Tekanan  uap  RMSE  1  =  0.1327  dan  IAE  1  =  32.1710,  Flow  uap  RMSE  4  =  1.7132  dan  IAE  4  = 2.3603, Temperatur uap RMSE 5 = 8.9758 dan IAE5 = 2242.4.
Kata  Kunci:  Kontrol  Model  Prediksi,  Optimisasi  Robas,  Optimisasi  Konik,  Robust  Counterpart, receding horizon, waste heat boiler.

## 1 Pendahuluan
Kontrol Model Prediksi biasanya dikenal dengan nama Model Predictive Control (MPC). MPC merupakan  metode  pengontrolan  yang  menggunakan  model  proses  untuk  menentukan sinyal kontrol dengan cara meminimalkan fungsi kinerja.
Di  dalam  strategi  kontrol  ini  sangat  diperlukan  adanya  metode  optimisasi.  Hal  ini dikarenakan  optimisasi  merupakan  bagian  penting  untuk  mendapatkan  sinyal  kontrol optimal.  Metode  yang  digunakan  pada  tugas  akhir  ini  adalah  optimisasi  robas.  Metode optimisasi  ini  diterapkan  dalam  menentukan  besar  sinyal  kontrol  sehingga  meminimalkan fungsi objektif pengontrolan yakni fungsi kinerja MPC.  Keungulan optimisasi robas mampu menangani  ketidakpastian  data.  Ketidakpastian  dalam  hal  ini  berarti  bahwa  data  suatu masalah  optimisasi  tidak  diketahui  secara  pasti  pada  saat  solusi  optimalnya  akan ditentukan.
Untuk melihat implementasi MPC dengan optimisasi robas  di industri, maka digunakanlah sistem Waste Heat Boiler (WHB) sebagai sistem yang akan dikontrol.

## 2 Waste Heat Boiler

## 2.1 Sistem Waste Heat Boiler
PT.Petrokimia  Gresik  memiliki  boiler  bertipe  waste  heat  boiler  (WHB).  Waste  heat  boiler adalah pembangkit uap yang memanfaatkan gas panas sisa dari generator turbin gas.

## Gambar 1 Sistem waste heat boiler

## 2.2 Identifikasi Sistem WHB
Identifikasi  dilakukan  dengan  pendekatan  black  box  yaitu  sistem  dimodelkan  dengan menggunakan  data  masukan  dan  data  keluaran  dari  data  pengukuran,  tanpa memperhatikan proses didalamnya.
Tabel 1  Masukan dan keluaran WHB Masukan Keluaran u1 Fuel flow, N m3/hr y1  Steam pressure, kg/cm2G u2 Boiler feedwater to HP drum, ton/hr y2  Furnace temperature, 0C u3 Boiler feedwater to superheater, ton/hr y3  HP drum level, % of high Ganggugan Terukur y4  Steam flow, ton/hr d1 Exhaust gas temperature, 0C y5  Steam temperature, 0C
Model  yang  digunakan  untuk  identifikasi  sistem  pada  tugas  akhir  in  adalah  model  ARX. Struktur  model  ARX  merepresentasikan  hubungan  masukan  keluaran  dalam  bentuk persamaan diferensi sebagai berikut
) ( ) ( ) ( ) ( ) ( t e t u q B t y q A     1 1
atau dapat ditulis sebagai
) ( ) ( ) ... ( ) ( ) .... ( t e t u q b q b t y q a q a nb nb na na            1 1 1 1 1 Gambar 2 Diagram blok struktur model ARX
Pendekatan  pada  sistem  WHB  dilakukan  dengan  pendekatan  MJKJ  pada  model  ARX, dengan  memanfaatkan  perangkat  lunak  Matlab.  Struk tur  model  yang  dipilih  dalam identifikasi adalah model multivariable ARX, dengan orde-orde polynomial:
na = 1, merepresentasikan jumlah pole dari masukan nb = 1, merepresentasikan (jumlah zero + 1) dari masukan dengan  nilai  orde  di  atas,  sesuai  dengan  persamaan  ,  hasil  identifikasi  WHB  dengan menggunakan  perangkat  lunak  Matlab,  diperoleh  parameter  polynomial  A  dan  B  sebagai berikut:
               
               9932 0 0455 0 0074 0 0118 0 0728 0 0067 0 9446 0 0208 0 0063 0 0807 0 0261 0 0199 0 9953 0 0083 0 01017 0 0624 0 0103 0 0111 0 9936 0 4095 0 0076 0 0100 0 0013 0 0011 0 9520 0 . . . . . . . . . . . . . . . . . . . . . . . . . A                          0007 0 3415 0 0109 0 0000 0 0082 0 1093 0 0139 0 0012 0 0001 0 1721 0 0077 0 0001 0 0015 0 0453 0 0094 0 0001 0 0018 0 0073 0 0043 0 0002 0 . . . . . . . . . . . . . . . . . . . . B

## 3 Kontrol Model Prediksi

## 3.1 Dasar – dasar MPC
Model yang digunakan pada M PC adalah bentuk ruang keadaan diskrit
) ( ) ( ) ( ) ( k d B k u B k x A k x d p p     1
) ( ) ( k x C k y p  ) ( ) ( k x C k z z  d p p B B A , , merupakan  matriks  masukan, p C merupakan  matriks  keluaran  teramati  dari system,  dan z C merupakan  matriks  keluaran  terkontrol  dari  sistem .  Kita  dapat memprediksikan himpunan output dari proses berdasarkan  iterasi model berikut
                            

                                               
                ) | ( ˆ ) | ( ˆ ) | ( ˆ ) | ( ˆ ) | ( ˆ ) | ( ˆ ) ( ) | ( ˆ ) | ( ˆ ) | ( ˆ k H k d k k d k k d B C B A C B A C B C B A C B C k H k u k k u k k u B C B A C B A C B C B A C B C k x A C A C A C k H k z k k z k k z v v v v H v v v H v v v v v v v v v v v v v H v v v H v v v v v v v v v H v v z v v v v z v v v v v 1 1 0 0 1 1 0 0 2 1 2 1 2 1                 Prediksi  masukan ) | ( k i k u  dapat  diekspresikan  dalam  prediksi  perubahan  masukan sebagai berikut:
                                                
    
            ) | ( ˆ ) | ( ˆ ) | ( ˆ ) ( ) ( ) ( ) | ( ˆ ) | ( ˆ ) | ( ˆ k H k u k k u k k u I I I I I I k u k u k u k H k u k k u k k u v z 1 1 0 0 1 1 1 1 1      
 
Representasi persamaan ruang keadaan menjadi model  prediksi  yang hanya mengandung U  .
) ( ) ( ) ( ) ( ) ( ) ( k d B k u I B k u k x I B A k u k x d v v v                                     0 1 0 1
Persamaan ruang keadaan ini dapat direpresentasikan dengan
) ( ) ( ) ( k d B u B k A k m      
 1
) ( ) ( k C k y  
Maka persamaan prediksi menjadi
                                                                                             ) | ( ) | ( ) | ( ) } ( ˆ ) } ( ˆ ) } ( ) ( ) | ( ˆ ) | ( ˆ ) | ( ˆ k H k d k k d k k d
CB B CA CAB CB CAB CB k H k u k k u k k u B CA B CA B CA CB CAB CB k CA CA CA k H k y k k y k k y v m m H m m m m u H H H H H v v u B B B B 1 1 0 0 1 1 0 0 2 1 2 2 2 2        
        
Atau
       g da akan masa lampau masa U D L k Y tan ˆ ˆ ) ( ˆ     

## 

## 3.2 Fungsi kinerja MPC
Fungsi kinerja diberikan oleh :
2 2 0 2 1 ) ( ) ( ) | ( ˆ ) ( ) | ( ˆ i R M i i n M i MPC v v v k k u i k r k i k y J             
dengan matriks pembobot
           ) ( ) ( v u H Q H Q Q      0 0            ) ( ) ( u H R R R      0 0 0 Disederhanakan menjadi 2 2 R Q MPC U J Y
J ˆ ˆ    
Penjejakan
D L k k J k ˆ ) ( ) ( ) (     
Dimana
merupakan  himpunan  penjejakan  referensi
) | ( ˆ k i k r  dan  didefinisikan
sebagai
  ) ( ) ( ) ( ) | ( / k y k w e k w k k r Tref TS       1 1
Fungsi kinerja menjadi
2 2 R Q mpc U U J ˆ ˆ        U H U G U J T T mpc ˆ ˆ ˆ konstanta      
  Q G T 2 
R Q H T    
Sehingga  untuk  meminimumkan  fungsi  kerja  sebanding  dengan  meminimumkan  fungsi kuadratik G U U H U T T     ˆ ˆ .  Fungsi  kuadratik  tersebut  dapat  diubah  menjadi
U G U H U T T ˆ ˆ ˆ     , karena matriks H adalah matriks Hessian. Kendala sistem tersebut dapat dituliskan dalam persamaan matematis, yaitu:
max min ) | ( ˆ u k i k u u      
max min ) | ( ˆ u k i k u u    
max min ) | ( ˆ x k i k x x   
max min ) | ( ˆ y k i k y y   
Persamaan diatas dapat dituliskan dalam bentuk umum
) ( ) ( ) ( ˆ k Xid k M k Fu U         1
dan dapat disederhanakan menjadi
   U ˆ

## 4 Optimisasi Robas

## 4.1 Optimisasi dengan Ketidakpastian
Permasalahan optimisasi diberikan dalam  bentuk:
    m i x f x f i o x ,......, , ) , ( : , min
1 0    
Dimana    , x f o merupakan  fungsi  objektif  dari  permasalahan  optimisasi  dan    , x f i dengan m i ,......, 1  ,  adalah  fungsi  kendala.  Fungsi  ini  tergantung  kepada  vector x dan  , dimana n R x  adalah vector variable keputusan dan N R   adalah vector yang berisi data masalah. Sebagai contoh dalam permasalahan optimisasi linier pada umumnya
  0   x b Ax x C T , : min ,
 merupakan tiga variable (c,A,b) dimana N R c  adalah vector objektif, A adalah matriks kendala m x n dan m R b  adalah vector sisi kanan pada kendala yang diberikan.
Untuk  menangani  ketidakpastian  digunakan  metoda  robast  counterpart,  sehingga diasumsikan  merupakan  anggota  himpunan  ketidakpastian u .    Dengan  demikian permasalahan optimisasi dapat diekspresikan sebagai
    u m i x f x f i o
      }, ,......, , , : , min{ 1 0
Dari  persamaan  diatas,  ketidakpastian  dapat  dikeluarkan  dari  fungsi  objektif  dari persamaan optimisasi dengan cara menggantinya dengan masalah yang ekivalen
u m i x f t x f t i o
       } ,..., , ) , ( , ) , ( : { min 1 0

## Sehingga RC memiliki solusi dari dengan meminimumkan   diberikan oleh
} , ,..., , ) , ( , ) , ( : ) , {( u m i x f t x f x t i o         1 0
Maka RC dari permasalahan ketidakpastian yang meminimumkan   menjadi
} , ,..., , ) , ( , ) , ( : { min u m i x f t x f t i o         1 0

## 4.2 Optimisasi Konik
Bentuk umum dari permasalahan optimisasi konik adalah:
} : { min k b Ax x C T R x h   
Fungsi objektif dari persamaan diatas adalah x C T , dimana vector objektif n R C  . b Ax  merepresentasikan fungsi affine dari n R ke m R , k menunjukan convex cone pada m R dan matriks kendala A berukuran m x n .
Jika  diasumsikan  ketidakpastian u diberikan  sebagai  bidang  affine  dari  himpunan  batas n R z   } {  , dan z diberikan oleh fungsi kendala:
1. Kendala linier, yaitu system dengan pertidaksamaan linier
2. P P  
3. Kendala kuadratik konik, yaitu system dengan pertidaksamaan kuadratik konik
4. m i r x q p P i T i i i ,...., , 1     
Kendala  linier  tersebut  ekivalen  dengan  linier  cone,  begitu  pun  dengan  kendala  kuadratik konik  sebanding  dengan  second-order  cone.  Dengan  demikian  dapat  dibentuk  dalam second-order cone programming dengan fungsi kendala linier, yaitu:
x b  min
m i r x q c x A c x A t s i i i , , . ,  1 0 1          

## 4.3 Robust Counterpart pada MPC dengan Ketidakpastian
MPC  dengan  Optimisasi  Robas  dalam  mengatasi  ketidakpastian  memerlukan  metoda Robust  Counterpart  yang  diusulkan  oleh  Ben-Tal  dan  Nemirovskii.  masalah  pokok  dalam metoda  ini  adalah  perumusan  RO  dari  masalah  ketidakpastian  sebagai  sebuah  masalah optimisasi  dapat  dikerjakan  dengan  mudah  atau  sebaliknya.  Oleh  karena  itu,  Ben-Tal  dan Nemirovskii  mengusulkan  memodelkan  ketidakpastian  ke  dalam  bentuk  khusus  yang disebut dengan ketidakpastian berbentuk elips.  Pendekatan ini menghasilkan keuntungan, karena  pemodelan  RC  akan  menghasilkan  permasalahan  khusus  dari  CO  seperti permasalahan  optimisasi  linier  (LO),  optimisasi  kuadratik  konik  (CQO)  atau  optimisasi semidefinite  (SDO),  yang  dapat  dengan  mudah  dipecahkan  oleh  metoda  interior  point. Ketidakpastian  berbentuk  elips  dengan  variasi  parameter  dilakukan  dengan  maksud mencari  parameter  ukuran  elips  yang  dapat  mendekati  ketidakpastian  sebenarnya  dan memberikan RO dengan solusi yang lebih optimal.
    U G U U G U H U T T        , ˆ : ˆ ˆ ˆ min  0 U G U U G U H U t t s t T T         
 ˆ ˆ ˆ ˆ . min    0 1      r w w w rQ G U n , , , (                     
  U zr U zr U G U H U t t s t T T n T ˆ ˆ ˆ ˆ ˆ . min 0 0 1 1 H dapat  disusun  sebagai S S H T  ,  sehingga U S U H U T ˆ ˆ ˆ     dan  diasumsikan
   U S ˆ .  Sehingga  persamaan  permasalahan  RC-M PC  dapat  dirum uskan  sebagai berikut:
                         U zr U U S zr U G t t s t T T n ˆ ˆ ˆ ˆ . min 0 0 1 1

## 5 Mpc Dengan Optimisasi Robas Pada Sistem Whb

## 5.1 Pemodelan Sistem
Model  sistem  yang  digunakan  adalah  model  persamaan  ruang  keadaan  dari  Waste  Heat Boiler  (WHB).  Seperti  yang  telah  dibahas  pada  bab  2,  persamaan  ruang  keadaan  sistem diskrit dapat ditulis berdasarkan  struktur model ARX dengan waktu  cuplik 1 menit sebagai berikut:
) ( ) ( ) ( ) ( k d Bd k Bu k Ax k x     1
) ( ) ( ) ( k Du k Cx k y  
Dimana kinerja matriks A,B,C dan D adalah
                           9559 0 0057 0 0216 0 0048 0 0154 0 0024 0 9922 0 0125 0 0002 0 0007 0 0013 0 0024 0 9441 0 0021 0 0009 0 0189 0 0028 0 4013 0 0072 1 0374 0 0564 0 0104 0 1790 0 0038 0 9792 0 . . . . . . . . . . . . . . . . . . . . . . . . . A                           0 2405 0 0107 0 0011 0 0 0408 0 0102 0 0000 0 0 0488 0 0052 0 0001 0 0 1127 0 0114 0 0000 0 0 2594 0 0064 0 0006 0 . . . . . . . . . . . . . . . . B                     0073 0 0 0 0 0 0005 0 0 0 0 0 0014 0 0 0 0 0 0007 0 0 0 0 0 0046 0 0 0 0 0 . . . . . Bd                              3082 0 4545 0 2510 0 5322 0 6180 0 8458 0 1453 0 0329 0 0507 0 4140 0 3791 0 8187 0 0860 0 0044 0 4587 0 1846 0 2889 0 3751 0 8607 0 4027 0 0172 0 0703 0 9353 0 0549 0 1532 0 . . . . . . . . . . . . . . . . . . . . . . . . . C
Dan sistem ini merupakan sistem MJMK dengan gangguan terukur pada masukan.

## 5.2 Pengujian dan Analisis Pengontrol
Parameter untuk pengontrol MPC  dengan optimisasi robas diberikan sebagai berikut:
Horizon kontrol = 3 Horizon prediksi = 5 Horizon minimum  = 1 Matriks pembobot keluar kedepan
                                 30 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 30 0 0 0 0 0 0 0 70 0 0 0 0 0 0 0 20 0 0 0 0 0 0 0 20 0 0 0 0 0 0 0 120     Q Matriks pembobot sinyal keluaran                        005 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 005 0 0 0 0 0 0 0005 0 0 0 0 0 0 5 0 . . . . .   R Kendala sistem : 480 5 108 4 80 3 670 2 66 1 0 5 0 4 0 3 0 2 0 1 575 9 12 3 108 2 4643 1 349 0 3 0 2 0 1                   max max max max max min min min min min max max max max min min min min , , , , , . , , , , , y y y y y y y y y y d u u u d u u u Ketidakpastian pada optimisasi robas                        6 0 6 0 6 0 67 0 6 0 . . . . .   Q
Dengan  selang  parameter  dimensi  dari  ketidakpastian  adalah
2 0   r .  Parameter ketidakpastian  dipilih  dengan  uji  dan  coba.  Pemilihan  nilai  diagonal Q semakin  mengecil dan  diatas  rata-rata  nilai  mutlak  rata-rata  masukan.  Hal  ini  bertujuan  agar  perubahan masukan  yang  dihasilkan  akan  semakin  mengecil.  Parameter  dimensi  elips  akan  dirubah dengan rentang 2 0 1 0   . . Sedangkan parameter untuk LQR adalah:
Gambar 3 Gangguan terukur dari panas sisa generator turbin
Gambar 4 Masukan terkontrol
Gambar 5 Nilai masukan sistem WHB dengan pengontrol MPC dan optimisasi robas
Tabel 2 Perbandingan MPC optimisasi robas dan data pabrik karakteristik MPC dengan Optimisasi Robas Data pabrik RMSE 1 0.1327 0.5474 RMSE 2 1.3808 2.7610 RMSE 3 9.0878 2.9364 RMSE 4 1.7132 2.3603 RMSE 5 8.9758 9.9719 IAE 1 32.1710 119.3452 IAE 2 341.5004 573.5380 IAE 3 2.2713e+003 730.2716 IAE 4 427.1013 549.4128 IAE 5 2.2424e+003 .4868e+003
Gambar 6 Grafik batang dari nilai kasus terburuk untuk semua parameter
Gambar 7 Perbesaran grafik batang dari nilai kasus terburuk untuk semua parameter
Gambar 8 Nilai parameter   yang paling meminimumkan  pada setiap waktu cacah
Gambar 9 Nilai kasus terburuk  pada setiap waktu cacah untuk ketidakpastian elips parametrik
Jika  kita  bandingkan  nilai  kasus  terburuk  dari  ketidakpastian  elips  parametrik  dengan ketidakpastian  elips  dengan  nilai  pertubasi  sama  dan  parameter  dimensinya  bernilai  1, akan memiliki perbedaan seperti pada gambar berikut:
Gambar 10 Nilai kasus terburuk  pada setiap waktu cacah untuk ketidakpastian elips parametrik
Gambar  9  dan  Gambar  10  menunjukan  bahwa  nilai  kasus  terburuk  pada  ketidakpastian elips parametrik lebih minimum jika dibandingkan dengan ketidakpastian elips.

## 6 Kesimpulan Dan Saran

## 6.1 Kesimpulan
1. Optimisasi robas dengan ketidakpastian elips parametrik dapat digunakan pada MPC
2. Pengontrol  MPC  dengan  optimisasi  robas  dapat  menjaga  sistem  tetap  stabil  terhadap gangguan dan ketidakpastian.
3. Pengontrol  MPC  dengan  optimisasi  robas  lebih  baik  jika  dibandingkan  dengan pengontrol yang diaplikasikan pada plant saat ini.
4. Dengan  kondisi  awal  tunak  pada  seluruh  keluaran  yang  mempengaruhi  kualitas  uap didapat:
Tekanan uap RMSE 1 = 0.1327 IAE 1 = 32.1710 Flow uap RMSE 4 = 1.7132 IAE 4 = 2.3603 Temperatur uap RMSE 5 = 8.9758 IAE5 = 2242.4 Dan nilai  RMSE dan IAE  ini lebih baik jika dibandingkan dengan proses pada plant saat ini. Sehingga  dapat  disimpulkan  pengontrol  MPC  sengan  optimisasi  robas  lebih  optimal dibandingkan dengan pengo ntrol yang diaplikasikan saat ini

## 6.2 Saran
Berdasarkan  pengalaman  pada  pelaksanaan  tugas  akhir  ini  dapat  dikemukakan  saran untuk pengembangan pengontrol M PC kedepannya,  yaitu metoda pemilihan ketidakpastian pada optimisasi robas untuk MPC.

## 7 Daftar Pustaka
[1] R.  I.  Simangunsong,  Kontrol  Model  Prediksi  dengan  Optimisasi  Robas  :  Studi  Kasus Double Rotary Inverted Pendulum , Tugas Akhir Teknik Fisika, ITB, Oktober 2008.
[2] J. M. Maciejowski, Predictive Control with Constraints, England: Prentice Hall, 2000.
[3] K. V. Ling, "Introduction to Model Predictive Control," course notes, Bandung Institute of Technology, May 2008.
[4] K.  V.  Ling  and  B.F.  Wu,  "Embedded  Model  Predictive  Control,"  Proceedings  of  EEE  The AUN/SEED-net Fieldwise Seminar on Control Engineering, 2007, pp. EL-2 – EL-6.
[5] W.  S.  Lu,  Use  SeDuMi  to  Solve  LP,  SDP,  SCOP  Problems  :  Remarks  and  Examples, Departement of Electrical and Computer Engineering, University of Victoria, 2007.
[6] Ben-Tal  and  A.  Nemirovski,  “Robust  Optimization  –  Methodology  and  Applications”, Mathematical Programming, 2002, vol. 92, no. 3 Seri B, pp. 453-480.
[7] D.  Chaerani,  Modelling  Robust  Design  Problems  via  Conic  Optimization,  PhD.  Thesis, Technische Universiteit Delft, 2006.
[8] PT. Petrokimia Gresik, Standard Fasel-Lenjes Instrument List, 1993.