Model Sistem dan Metode Estimasi Parameter Rekursif

Bab 2 Model Sistem dan Metode Estimasi Parameter Rekursif

Tujuan Pengajaran Tujuan pengajaran pada bab ini adalah untuk memahami cara memodelkan sistem untuk sistem kendali swa-tala dengan metode kuadrat terkecil, extended least-squares, dan instrumental variable. Model sistem ini berhubungan dengan sistem fisik melalui proses linierisasi dan diskritisasi dengan menggunakan transformasi z. Tujuan yang kedua adalah untuk mempelajari teknik-teknik estimasi parameter model sistem. Hal ini dilengkapi dengan mengasumsikan bentuk waktu diskrit untuk model sistem dan kemudian menggunakan suatu algoritma estimasi rekursif untuk menghasilkan estimasi parameter model.

2.1 Model Sistem Untuk Estimasi Rekursif Untuk data identifikasi dengan menggunakan komputer sinyal keluaran proses dicuplik dan dikuantisasi oleh analog-digital converter (ADC). Hasil pengolahan ini berupa sinyal diskrit. Diasumsikan bahwa unit amplitudo dan waktu kuantisasi kecil maka sinyal diskrit tersebut dapat dianggap sebagai sinyal kontinu. Jika pencuplikan diaktifkan secara periodik dengan waktu pencuplikan h, maka serangkaian pulsa dibangkitkan setelah keluar dari ADC seperti terlihat pada gambar 2.1. Sinyal terdigitalisasi ini kemudian diproses oleh komputer digital dengan menggunakan algoritma terprogram, dan sinyal keluaran dihasilkan. Jika komputer digital digunakan untuk mengendalikan proses dan sinyal analog dibutuhkan oleh aktuator, maka sinyal keluaran terhitung oleh algoritma tadi disalurkan melewati digital-analog converter (DAC) dan diikuti oleh divais zeroorder-hold (zoh) untuk diubah menjadi sinyal analog. Waktu yang diperlukan untuk mengubah sinyal analog menjadi sinyal digital atau sebaliknya didalam struktur sistem kendali digital dapat diabaikan, karena nilainya jauh lebih kecil dibandingkan konstanta waktu aktuator, sensor, dan proses.

16

Bab 2: Metode Estimasi Parameter Rekursif

Diasumsikan bahwa sistem stabil adalah time invariant dan terlinierisasi, sehingga sistem tersebut dapat dinyatakan oleh persamaan difference y (t ) + a1 y (t − 1) + L + ana y (t − na ) = b0 u(t ) + b1u(t − 1) + L + abnb u(t − nb ) (2.1) Dengan mendefinisikan operator backward shift z-1 sebagai z −1 y (t ) = y (t − 1)

(2.2)

maka model persamaan (3.1) dapat diekspresikan didalam fungsi alih diskrit dalam bentuk Gp (z) =

y (t ) B( z −1 ) = u(t ) A( z −1 )

(2.3)

dengan A dan B adalah polinomial dalam operator z-1 A( z −1 ) = 1 + a1 z −1 + L + ana z − na B( z −1 ) = b0 + b1 z −1 + L + bnb z −nb

nh

(2.4)

nh

ud(t)

yd(t)

y(t)

u(t) Zero order hold

Control algorithm

Sampler & ADC

y(t) Process

Sampler & DAC

Gambar 2.1 Sistem kendali dengan komputer sebagai pengendali cuplik data. Sistem dinamik mempunyai respon yang tidak spontan terhadap sinyal masukan yang diberikan, sehingga umumnya fungsi alih pada persamaan (2.3) juga mengandung unsur satu waktu tunda Gp (z) =

y (t ) z −1B( z −1 ) = u (t ) A( z −1 )

(3.5)

Sinyal keluaran proses terukur y(t) diasumsikan terkontaminasi oleh gangguan n(t) seperti terlihat pada gambar 2.2 y ( t ) = y u ( t ) + n(t )

(2.6)


17

Sinyal gangguan tersebut diasumsikan dapat diwakili model autoregressive moving average (ARMA)

n(t ) + d1n(t − 1) + L + d p n(t − p ) = e(t ) + c1e(t − 1) + L + c p e(t − p )

(2.7)

dengan e(t) adalah derau yang tak terukur, terdistribusi normal, dan bebas linier yang karakterisasinya dinyatakan dalam

E {e(t )} = 0

cov[e(t ),τ ] = E {e(t )e(t + τ } = σ e2δ (t ) dengan σ e2 adalah variansi dan δ (τ ) adalah fungsi impuls diskrit ⎧1, ⎩0,

δ (τ ) = ⎨

τ =0 τ ≠0

e(t)

C(z-1) D(z-1)

u(t)

z-1B(z-1)

yu(t)

y(t)

A(z-1)

Gambar 2.2 Model proses dan derau

Fungsi alih filter derau ini adalah −1 −p n(t ) C( z −1 ) 1 + c1 z + L + c p z Ge ( z ) = = = e(t ) D( z −1 ) 1 + d1 z −1 + L + d p z − p

(2.8)

Persamaan (3.5) dan (3.8) menghasilkan kombinasi model proses dan derau y (t ) =

B( z −1 ) −1 C( z −1 ) z u ( t ) + e(t ) A( z −1 ) D( z −1 )

(2.9)

Tujuan identifikasi model proses adalah untuk mengestimasi parameterparameter proses dalam polinomial A(z-1) dan B(z-1) dan parameter derau C(z-1) dan D(z-1) dengan berbasiskan sinyal masukan u(t) dan keluaran y(t) terukur. Sebagai pengetahuan a priori diasumsikan orde model na, nb, dan p diketahui nilainya. Gangguan n(t) diasumsikan stasioner, yaitu akar-akar D(z-1) terletak

18


didalam unit lingkaran bidang-z. Jika polinomial D(z-1) sama dengan polinomial A(z-1) proses maka model ini disebut sebagai autoregressive moving average exogenous (ARMAX) y (t ) =

B( z −1 ) −1 C( z −1 ) z u ( t ) + e(t ) A( z −1 ) A( z −1 )

(2.10)

Model ARMAX mempunyai bentuk khusus untuk C(z-1) = 1 yang disebut sebagai model least-squares (LS) y (t ) =

1 B( z −1 ) −1 z u (t ) + e(t ) −1 A( z ) A( z −1 )

(2.11)

2.2 Metode Kuadrat-Terkecil Rekursif Identifikasi adalah penentuan karakteristik dinamik proses secara eksperimental berdasarkan sinyal-sinyal yang terukur kedalam suatu model matematik. Tujuan metode identifikasi model ini adalah meminimumkan selisih antara sinyal keluaran proses sebenarnya dengan keluaran model matematis. Didalam sistem kendali swa-tala proses identifikasi ini dilakukan secara on-line, dengan pengertian bahwa komputer beroperasi secara on-line dengan proses. Setiap penerimaan data baru langsung diolah oleh algoritma terprogram pada setiap periode pencuplikan, sehingga sistem kendali swa-tala termasuk kedalam katagori sistem waktu nyata (real time systems). Salah satu metode yang paling banyak dipakai dalam estimasi parameter adalah metode kuadrat terkecil (least-squares). Metode kuadrat terkecil pertama kali ditemukan oleh seorang matematikawan Jerman Gauss pada abad 18. Oleh Gauss metode ini digunakan untuk mengestimasi orbit planet-planet dari sistem tata surya. Inti dari metode ini adalah bahwa kecocokan antara model dengan sistem yang akan diidentifikasi diperoleh dengan meminimumkan selisih kuadrat antara keluaran model dengan keluaran sistem yang diidentifikasi untuk semua N data pengamatan N

N

i =1

i =1

J = ∑∈2 (t i ) = ∑ [y (t i ) − yˆ (t i )]

2

(2.12)

Gambar 2.3 menggambarkan ide metode ini. Terlihat bahwa antara sistem yang diidentifikasi dengan model yang digunakan terhubung secara pararel. Bentuk skema identifikasi ini tidak cocok dipakai untuk sistem dinamik, karena bentuk persamaan minimasi yang diperoleh menjadi non linier. Dari gambar 2.3 sinyal kesalahan estimasi mempunyai bentuk persamaan:


19

n(t) Proses u(t)

+

yu(t)

z-1B(z-1)

y(t)

A(z-1) Model ^

z-1B(z-1)

^ y(t)

-

kesalahan

+

^ -1 ) A(z Optimasi

Gambar 2.3 Skema identifikasi proses secara pararel

e(t ) = y (t ) −

Bˆ ( z −1 ) u (t ) Aˆ ( z −1 )

(2.13)

Hubungan antara parameter dari polinomial B(z-1) dan parameter dari polinomial A(z-1) adalah non linier, dan menyulitkan dalam solusi minimasi fungsi kriteria. Umumnya dibutuhkan metode yang rumit untuk menyelesaikannya, contohnya adalah seperti metode pemrograman dinamik. Oleh karena itu dipakai bentuk skema lain seperti yang terlihat pada gambar 2.4

n(t) u(t)

+

z-1B(z-1) A(z-1)

^ B(z-1)

-

+

y(t)

^ A(z-1)

∈(t)

Gambar 2.4 Skema identifikasi berdasarkan kesalahan prediksi

Bentuk persamaan kesalahannya adalah

20


∈ (t ) = Aˆ ( z −1 )y (t ) − z −1Bˆ ( z −1 )u(t ) = 1 + aˆ1 z −1 + L + aˆ na z −na y (t ) − z −1 bˆ0 + bˆ1 z −1 + L + bˆnb z −nb u(t ) = y (t ) − − aˆ1 y (t − 1) − L − aˆ na y (t − na ) + bˆ0 u(t − 1) + L + bˆnb u(t − nb − 1) 144444444444424444444444443

(

(

)

(

)

)

yˆ ( t |t −1)

(2.14) Diantara parameter-parameter yang akan diestimasi terlihat memiliki hubungan yang linier, sehingga memudahkan dalam mencari persamaan penyelesaiannya. Persamaan residu di atas merupakan persamaan kesalahan prediksi, yaitu antara keluaran proses yang baru dengan keluaran proses terprediksi satu langkah kedepan = ∈(t) kesalahan = prediksi

y(t) pengukuran baru

-

y(t|t-1) prediksi satu langkah kedepan (2.15)

Kesalahan ini dapat disebabkan karena keluaran sistem y(t) yang terkontaminasi derau dan juga akibat kesalahan estimasi parameter. Persamaan keluaran terprediksi yˆ (t | t − 1) =

− aˆ1 y (t − 1) − L − aˆ na y (t − na ) + bˆ0 u(t − 1) + L + bˆnb u(t − nb − 1) ⎡ aˆ1 ⎤ ⎢ M ⎥ ⎢ ⎥ ⎢aˆ ⎥ = [− y (t − 1) L − y (t − na ) − u(t − 1) L − u(t − nb − 1)]⎢ na ⎥ ˆ ⎢ b0 ⎥ ⎢ M ⎥ ⎢ ⎥ ⎣⎢bˆnb ⎦⎥

=

Ψ (t )θˆ(t − 1) T

(2.16) dengan

Ψ (t ) = [− y (t − 1) L − y (t − na ) u(t − 1) L u(t − nb − 1)] T

(2.17)

adalah vektor data dan

θˆ (t ) = [a1 L ana T

b0 L bnb ]

(2.18)

adalah vektor parameter. Sekarang diasumsikan sistem telah berjalan untuk t = 1, 2, ..., m+N, dengan


m = max{na,nb}

21

(2.19)

Maka terdapat N+1 persamaan yang mempunyai bentuk y (t ) = Ψ (t )θˆ(t − 1)+ ∈ (t ) T

yang dapat direpresentasikan sebagai persamaan vektor y (m + N ) = Ψ(m + N )θˆ(m + N − 1) + ∈(m + N )

(2.20)

dengan ⎡ y (m ) ⎤ ⎢ y (m + 1) ⎥ ⎥ y (m + N ) = ⎢ ⎥ ⎢ M ⎥ ⎢ ⎣ y (m + N )⎦

(2.21)

⎡ ψ (m ) ⎤ ⎢ ψ (m + 1) ⎥ ⎢ ⎥ ⎢ ⎥ M ⎥ ⎢ ⎢⎣ψ (m + N )⎥⎦ u(m − 1) L − y (1) u(m ) L − y ( 2)

Ψ( m + N ) =

u( m − 2) L u(1) ⎤ − y ( m − 2) ⎡ − y (m − 1) ⎢ u(m − 1) L u( 2) ⎥⎥ − y (m ) − y (m − 1) = ⎢ ⎢ M M O M M M O M ⎥ ⎥ ⎢ ⎣− y (m + N − 1) − y (m + N − 2) L − y (N ) u(m + N − 1) u(m + N − 2) L u(N )⎦ (2.22) ⎡ ∈ (m ) ⎤ ⎢ ∈ (m + 1) ⎥ ⎥ ∈(m + N ) = ⎢ ⎥ ⎢ M ⎥ ⎢ ⎣∈ (m + N )⎦

(2.23)

Fungsi kriteria sekarang merupakan perkalian vektor ∈ dengan vektor transposnya ∈T J = ∈ (m + N )∈(m + N ) = T

m +N

∑∈

2

(t )

(2.24)

t =m

Dengan menguraikan definisi vektor ∈ pada persamaan (2.20) kedalam persamaan (2.24) diperoleh persamaan

22


J

(y − Ψθ ) (y − Ψθ ) = T T T = y − θ Ψ (y − Ψθ ) T T T T T T = y y − y Ψθ − θ Ψ y + θ Ψ Ψθ T

(

)

(2.25)

Minimasi fungsi kriteria ∂J T T T T = −Ψ y − Ψ y + Ψ Ψθ + Ψ Ψθ = 0 ∂θ

(2.26)

menghasilkan bentuk persamaan estimasi parameter

θˆ(m + N − 1) = P (m + N )Ψ (m + N )y (m + N ) T

(2.27)

dengan P adalah matriks varians yang didefinisikan sebagai

[

]

P (m + N ) = Ψ (m + N )Ψ(m + N ) T

−1

(2.28)

Contoh 2.1 Sebuah proses dinamik yang tidak diketahui diidentifikasi dengan menggunakan model y (t ) = −aˆ1 y (t − 1) − aˆ 2 y (t − 2) + bˆ0 u(t − 2)

Data pengukuran dari sinyal masukan dan keluaran proses diberikan pada tabel berikuti ini: t u(t) y(t)

1 10 0

2 -10 0

3 10 1

4 -10 0.9

5 10 1.81

6 -10 1.629

Dengan menggunakan metode kuadrat terkecil, tentukanlah nilai parameter proses yang tidak diketahui! Solusi

Persamaan model ditulis kembali kedalam bentuk perkalian vektor ⎡ a1 ⎤ y (t ) = [− y (t − 1) − y (t − 2) u(t − 2)]⎢⎢a 2 ⎥⎥ ⎢⎣b0 ⎥⎦


23

Terlihat bahwa keluaran proses bergantung kepada nilai sebelumnya hingga sampai 2 kali pencuplikan. Dengan demikian isi vektor keluaran y baru dapat dimulai dari pencuplikan ketiga ⎡ y (3)⎤ ⎡ 1 ⎤ ⎢ y ( 4)⎥ ⎢ 0.9 ⎥ ⎥=⎢ ⎥ y=⎢ ⎢ y (5)⎥ ⎢ 1.81 ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ y (6)⎦ ⎣1.629⎦ dan matriks ψ berisi ⎡− y ( 2) − y (1) ⎢ − y (3) − y ( 2) Ψ=⎢ ⎢− y ( 4) − y (3) ⎢ ⎣ − y (5) − y ( 4)

u(1) ⎤ ⎡ 0 0 10 ⎤ ⎢ ⎥ u( 2)⎥ ⎢ − 1 0 − 10⎥⎥ = u(3)⎥ ⎢ − 0.9 10 ⎥ −1 ⎥ ⎥ ⎢ u( 4)⎦ ⎣− 1.81 − 0.9 − 10⎦

Dengan menggunakan persamaan (2.27) untuk menghitung nilai parameter yang dicari, maka ⎡ aˆ1 ⎤ ⎡ 1.9972 − 2.8472 − 0.1025⎤ θˆ = ⎢⎢aˆ 2 ⎥⎥ = ⎢⎢− 2.8472 4.6122 0.1475 ⎥⎥ ⎢⎣bˆ0 ⎥⎦ ⎢⎣ − 0.1025 0.1475 0.0078 ⎥⎦ =

−1

⎡− 1.9⎤ ⎢ 0.9 ⎥ ⎢ ⎥ ⎢⎣ 0.1 ⎥⎦

⎡ 1 ⎤ ⎡ 0 − 1 − 0.9 − 1.81⎤ ⎢ ⎥ ⎢0 ⎥ ⎢ 0.9 ⎥ − − 0 1 0 . 9 ⎢ ⎥ ⎢ 1.81 ⎥ ⎢⎣10 − 10 10 10 ⎥⎦ ⎢ ⎥ ⎣1.629⎦

2.2.1 Kriteria Evaluasi Estimator Disamping metode kuadrat-terkecil yang sering digunakan dalam identifikasi sistem, juga dikenal metode lain seperti algoritma extended least-squares, intrumental variable, dan maximum likelihood. Untuk mengetahui metode identifikasi terbaik maka diperlukan indikator kinerja sebagai acuan untuk mengevaluasi hasil kuaslitas estimasi. Di bawah ini adalah 3 definisi yang dijadikan bahan acuan untuk mengevaluasi estimator.

Definisi 2.1: Unbiased Pada beberapa kali estimasi akan diperoleh vektor parameter model proses hasil estimasi θˆ yang berbeda. Estimator disebut unbiased jika nilai rata-rata dari vektor parameter terestimasi sama dengan nilai vektor parameter sebenarnya θ. E{θˆ } = E{θ} (2.29)

24


Gambar 2.5 Kriteria estimasi yang unbiased Definisi 2.2: Effisien Estimator dinilai efisien, jika nilai kovarian kesalahan estimasi untuk sebuah model proses yang diberikan yang berdasarkan kepada N data pengukuran masukan dan keluaran proses adalah minimal, artinya tidak ada estimator lain yang sanggup memberikan nilai kovarian kesalahan estimasi lebih kecil lagi dibawah kondisi berikut ini

Cθ~θ~ = E {(θ − θˆ ).(θ − θˆ )T } ⇒ min

Gambar 2.6 Kriteria efisiensi estimator

(2.30)


25

Gambar 2.7 Kriteria konsistensi estimator Definisi 2.3: Konsisten Dengan bertambahnya data pengukuran maka kovarian kesalahan estimasi menjadi lebih kecil. Untuk N → ∞ nilai estimasi bergerak menuju nilai sebenarnya.

lim θˆ = θ ,

N →∞

~ lim Cθ~θ~ = 0 , dengan θ = θˆ − θ

N →∞

(2.31)

Untuk N Æ ∞ nilai estimasi θˆ mendekati nilai sebenarnya θ . Hal ini berarti bahwa peningkatan jumlah data untuk estimasi juga memperbaiki hasil estimasi.

Contoh 2.2 Diketahui suatu proses yang akan diidentifikasi mempunyai fungsi alih y (t ) =

0.5z −1 u ( t ) + v (t ) 1 − 0.97z −1

Simulasikan sistem ini dengan menggunakan program bantu MATLAB jika sistem diberikan masukan berupa sinyal pseudo random binary signal (PRBS) dengan jumlah pencuplikan sebanyak 1000. Sinyal gangguan v(t) adalah sinyal derau putih e(t) yang tertapis oleh lowpass filter: v (t ) =

1 e(t ) 1 − 0.9z −1

26


Atur varian sinyal gangguan v(t) sedemikian, sehingga sinyal keluaran mempunyai noise-to-ratio sebesar 1%, dan 10% var (v (t )) = 1%,10% var (y (t )) Program simulasi: %Create a PRBS signal N=1000; u = sign(randn(N,1)); %Simulate the process A0 = [1 -0.97]; B0 = [0 0.5]; y0 = filter(B0,A0,u); % Add 1% and 10% output disturbances v = filter(1,[1 -0.9],randn(N,1)); v01 = v/std(v)*sqrt(0.01)*std(y0); y01 = y0+v01; v10 = v/std(v)*sqrt(0.1)*std(y0); y10 = y0+v10; % Estimate 1st order ARX models thARX01 = arx([y01 u],[1 1 1]); thARX10 = arx([y10 u],[1 1 1]); % Simulation of the models yh01 = idsim(u,thARX01); yh10 = idsim(u,thARX10); % Calculate step responses Nstp=200; STP0 = filter(B0,A0,ones(Nstp,1)); STP01 = idsim(ones(Nstp,1),thARX01); STP10 = idsim(ones(Nstp,1),thARX10); % Plot model fit k=1:N; figure; subplot(211), plot(k,y01,'r--',k,yh01,'b-'); title('Fit of 1st order ARX model'); axis([0 1000 -8 8]); xlabel('Samples'); subplot(211), plot(k,y10,'r--',k,yh10,'b-'); title('Fit of 1st order ARX model'); axis([0 1000 -8 8]); xlabel('Samples'); % Plot step responses figure; k=1:Nstp; subplot(211), plot(k,STP0,'r--',k,STP01,'b-'); title('Step responses of the process (dashed) and of ARX model, 1% noise'); xlabel('Samples'); subplot(212), plot(k,STP0,'r--',k,STP10,'b-'); title('Step responses of the process (dashed) and of ARX model, 10% noise'); xlabel('Samples');


Gambar 2.8 Perbandingan kurva keluaran sistem dengan keluaran model

Gambar 2.9 Keluaran sistem dan model untuk masukan fungsi step

27

28


Dari hasil simulasi terlihat bahwa model yang diperoleh dari hasil identifikasi dengan algoritma kuadrat terkecil cukup baik meniru karakteristik sistem terutama untuk nilai variansi derau putih yang tidak besar.

2.2.2 Metode Kuadrat-Terkecil Rekursif Bentuk persamaan (2.27) adalah estimator parameter non rekursif, dimana parameter terestimasi diperoleh hanya setelah seluruh pengukuran data dilakukan. Agar estimator dapat digunakan dalam sistem kendali swa-tala, maka bentuk persamaannya harus diubah menjadi iteratif, sehingga parameter model proses dapat diperbaharui setiap periode pencuplikan untuk setiap perolehan data baru. Untuk estimator dengan menggunakan data dari waktu 1 ke t mempunyai bentuk persamaan

[

]

−1

θˆ(t ) = Ψ (t )Ψ(t ) Ψ (t )y (t ) T

T

(2.32)

Pada waktu t+1 diperoleh data pengukuran yang baru dari proses, sehingga dapat didefinisikan ⎡ ψ T (m ) ⎤ ⎥ ⎢ T ⎢ψ (m + 1)⎥ ⎡ Ψ(t ) ⎤ ⎥ ⎢ ⎢ M ⎥ ⎥=⎢ L ⎥ Ψ(t + 1) = ⎢ T ⎢ ψ (t ) ⎥ ⎢ T ⎥ ⎥ ⎣ψ (t + 1)⎦ ⎢ L ⎥ ⎢ T ⎢⎣ ψ (t + 1) ⎥⎦

(2.33)

⎡ y (m ) ⎤ ⎢ y (m + 1)⎥ ⎢ ⎥ ⎡ y (t ) ⎤ ⎢ ⎥ ⎢ M ⎥ y (t + 1) = ⎢ ⎥=⎢ L ⎥ ⎢ y (t ) ⎥ ⎢ y (t + 1)⎥ ⎦ ⎢ L ⎥ ⎣ ⎢ ⎥ ⎢⎣ y (t + 1) ⎥⎦

(2.34)

Persamaan estimasi parameter pada waktu t+1 sekarang mempunyai bentuk

[

]

−1

θ (t + 1) = Ψ (t + 1)Ψ(t + 1) Ψ (t + 1)y (t + 1) T

T

(2.35)

Suku yang berada pada bagian invers dapat diuraikan kembali menurut persamaan (2.36)


Ψ (t + 1)Ψ(t + 1) = T

29

⎡ Ψ(t ) ⎤ ⎢ ⎥ Ψ (t ) M ψ (t + 1) ⎢ L ⎥ ⎢ψ T (t + 1)⎥ ⎣ ⎦ T T Ψ (t )Ψ(t ) + ψ (t + 1)ψ (t + 1)

[

=

]

T

(2.36)

Dengan cara yang sama suku sisa dari persamaan (2.35) dapat ditulis kembali menjadi

Ψ (t + 1)y (t + 1) = T

=

⎡ y (t ) ⎤ Ψ (t ) M ψ (t + 1) ⎢⎢ L ⎥⎥ ⎢⎣ y (t + 1)⎥⎦ T Ψ (t )y (t ) + ψ (t + 1)y (t + 1)

[

]

T

(2.37)

Untuk menghindari waktu komputasi yang besar pada persamaan estimasi parameter akibat adanya suku invers, maka diperlukan cara untuk memperbaharui invers matriks ΨT(t+1)Ψ(t+1). Matriks kovarian pada waktu t berdasarkan persamaan (2.28) dapat ditulis kembali dalam bentuk

[

]

P (t ) = Ψ (t )Ψ(t ) T

−1

(2.38)

Invers matriks ini pada waktu t+1 adalah −1

−1

P (t + 1) = P (t ) + ψ (t + 1)ψ (t + 1) T

(2.39)

Dengan memanfaatkan teorema invers matriks Lemma

(A + BC D )−1 = A −1 − A −1 B (C −1 + D A −1 )

−1

DA

−1

dan menetapkan −1

A = P (t ),

C = 1,

B = ψ (t + 1)

D = ψ (t + 1) T

memberikan persamaan

(

)

−1 T T P (t + 1) = P (t )⎡I m − ψ (t + 1) 1 + ψ (t + 1)P (t )ψ (t + 1) ψ (t + 1)P (t )⎤ ⎢⎣ ⎥⎦

(2.40)

Untuk sistem dengan satu-masukan-satu-keluaran (SISO) suku invers pada persamaan (2.40) memberikan hasil skalar. Dengan demikian persamaan untuk memperbaharui matriks kovarian tidak lagi mengandung suku invers.

30


Langkah selanjutnya adalah membuat persamaan estimasi parameter bergantung kepada hasil estimasi sebelumnya. Persamaan (2.35) ditulis kembali menjadi

θˆ(t + 1) = P (t + 1)B(t + 1)

(2.41)

θˆ(t ) = P (t )B(t ) dengan B( t ) = Ψ ( t ) y ( t ) T

(2.42)

B(t + 1) = Ψ (t + 1)y (t + 1) = B(t ) + ψ (t + 1)y (t + 1) T

Dengan mendefinisikan e(t+1) sebagai kesalahan prediksi ∈ (t + 1) = y (t + 1) − ψ (t + 1)θˆ(t ) T

(2.43)

dan mensubstitusikan kedalam persamaan (2.41) memberikan persamaan T B(t + 1) = B(t ) + ψ (t + 1)ψ (t + 1)θˆ(t ) + ψ (t + 1) ∈ (t + 1)

(2.44)

Substitusi persamaan (2.44) kedalam persamaan (2.41) diperoleh persamaan

θˆ(t + 1) = θˆ(t ) + P (t + 1)ψ (t + 1)∈(t + 1)

(2.45)

1442443 γ ( t +1)

dengan γ(t+1) adalah vektor pengkoreksi. Dengan demikian persamaan estimasi dapat ditulis kembali menjadi

[

]

θˆ(t + 1) = θˆ(t ) + γ (t + 1) y (t + 1) − ψ (t + 1)θ (t ) T

(2.46)

Karena pencuplikan dimulai dari waktu t=1 maka persamaan (2.43), (2.40), dan (2.46) ditulis kembali menjadi T ∈ (t ) = y (t ) − ψ (t )θˆ(t − 1)

(

(2.47)

)

−1 T T P (t ) = P (t − 1)⎡I m − ψ (t ) 1 + ψ (t )P (t − 1)ψ (t ) ψ (t )P (t − 1)⎤ ⎢⎣ ⎥⎦

(2.48)


[

31

]

θˆ(t ) = θˆ(t − 1) + γ (t ) y (t ) − ψ (t )θ (t − 1) T

(2.49)

Untuk memulai algoritma rekursif dipilih harga awal vektor parameter dan matriks kovarian

θˆ(0) =

0 P ( 0) = α I

dengan α adalah konstanta yang bernilai tinggi. Nilai ekspektasi matriks P proporsional terhadap matriks kovarian parameter terestimasi E {P (t )} =

1

σ ∈2

cov[Δθ (t − 1)]

(2.50)

dengan

{ }

σ ∈2 = E ∈ ∈ T

Contoh 2.3 Untuk contoh 2.1 carilah nilai parameter model proses yang tidak diketahui dengan menggunakan metode kuadrat-terkecil rekursif dengan harga awal: P(0) = 1000.I dan θ(0) = 0!

Solusi Agar vektor data dapat diisi sempurna dengan data yang tersedia, maka algoritma rekursif baru dapat dijalankan pada pencuplikan ketiga. Nilai awal matriks kovarian dan vektor parameter tetap tidak berubah. Buatlah program simulasi dengan MATLAB untuk menampilkan estimasi parameter, nilai diagonal matriks kovarian, dan nilai keluaran prediksi satu-langkah-kedepan hingga t=200. Dari data yang diberikan dibentuk tabel data proses pada setiap pencuplikan: t 1 2 3 4 5 6

y(t) 0 0 1 0.9 1.81 1.629

y(t-1) 0 0 1 0.9 1.81

y(t-2) 0 0 1 0.9

u(t-2) 10 -10 10 -10

32


Pencuplikan t=3:

Vektor data ψ (3) = [0 0 10] . Menghitung kesalahan T

prediksi dengan menggunakan persamaan (2.47): ∈ (3) = y (3) − ψ (3)θ (2) ⎡0⎤ = 1 − [0 0 10]⎢⎢0⎥⎥ ⎢⎣0⎥⎦ = 1 T

Karena sebelum pencuplikan t=3 tidak dilakukan estimasi parameter, maka P(2) = P(0). Dengan demikian nilai baru matriks kovarian P(3) dengan menggunakan persamaan (2.48) adalah: T ⎡ ψ (3)ψ (3)P (2) ⎤ ⎥ P (3) = P (2)⎢I 3 − T 1 + ψ (3)P (2)ψ (3) ⎥⎦ ⎢⎣ 0 0 ⎤ ⎡1000 ⎢ 1000 0 ⎥⎥ = ⎢ 0 ⎢⎣ 0 0 0.01⎥⎦

Jadi nilai estimasi parameter berdasarkan persamaan (2.49) adalah

θˆ(3) = θˆ(2) + P (3)ψ (3) ∈ (3) T

=

⎡0⎤ ⎢0⎥ ⎢ ⎥ ⎢⎣0.1⎥⎦

Pencuplikan t=4: Vektor data ψ ( 4) = [− 1 0 − 10] . Dengan cara yang sama T

seperti pada pencuplikan ketiga, maka diperoleh:

∈ ( 4) = 1.9

0 − 0.1⎤ ⎡1.996 ⎢ P ( 4) = ⎢ 0 1000 0 ⎥⎥ ⎢⎣ − 0.1 0 0.01⎥⎦

Program simulasi: >> >> >> >> >> >>

v=0.0; u = sign(randn(200,1)); e = v*randn(200,1); th0 = idpoly([1 -1.9 0.9],[0 0 0.1]); y = sim(th0,[u e]); z = iddata(y,u);

⎡− 1.8924 ⎤ ⎥ θ ( 4) = ⎢⎢ 0 ⎥ ⎢⎣ 0.1 ⎥⎦


33

>> >> >> >> >>

alpha = 1e3; lambda = 1; [th,yh,p,phi]=rarx(z(1,:),[2 1 2],’ff’,lambda,zeros(3,1),alpha.*eye(3)); theta(1,:) = th; cov(:,:,1)=p; yhat(1)=yh; for k = 2:200 [th,yh,p,phi] = rarx(z(k,:),[2 1 2],'ff',lambda,th',p,phi); theta(k,:) = th; cov(:,:,k) = p; yhat(k) = yh; >> end

Gambar 2.5 dan 2.6 memperlihatkan hasil simulasi kurva estimasi parameter, elemen diagonal matriks kovarian, dan sinyal prediksi satu-langkah-kedepan. Dalam simulasi diasumsikan nilai sebenarnya dari proses adalah: a1 = -1.9, a2 = 0.9 dan b0 = 0.1. Dalam pengukuran tidak terdapat derau yang mengkontaminasi keluaran proses. Untuk kasus dimana proses adalah noise-free, maka algoritma kuadrat-terkecil rekursif memberikan hasil yang unbiased terhadap nilai parameter proses yang sebenarnya, seperti terlihat pada gambar (2.10) atas. Hasil kecocokan kurva antara kurva keluaran proses sebenarnya dengan keluaran prediksi satu-langkah-kedepan memperkuat pernyataan tersebut, seperti terlihat pada gambar 2.11.

Gambar 2.10 (atas) Kurva estimasi parameter, (bawah) diagonal matriks kovarian

34


Gambar 2.11 Kurva perbandingan antara sinyal keluaran proses dengan keluaran prediksi satu-langkah-kedepan

Contoh 2.4 Diketahui suatu proses yang akan diidentifikasi mempunyai fungsi alih y (t ) =

z −1 + 0.5z −1 u (t ) + v ( t ) 1 − 1.5z −1 + 0.7z −2

Simulasikan sistem ini dengan menggunakan program bantu MATLAB jika sistem diberikan masukan berupa sinyal pseudo random binary signal (PRBS) dengan jumlah pencuplikan sebanyak 200. Sinyal gangguan v(t) adalah sinyal derau putih e(t) yang tertapis oleh lowpass filter: v (t ) =

1 − 1.8z −1 + 0.81z −2 e(t ) 1 − 1.5z −1 + 0.7z −2

Atur varian sinyal gangguan v(t) sedemikian, sehingga sinyal keluaran mempunyai noise-to-ratio sebesar 1%, dan 10% var (v (t )) = 1%,10% var (y (t ))

Program simulasi: echo off, clear, clc niter = 3; orde=2;


N = 200; thau=0.03; NStp=70; % Create the PRBS signal and random signal u = idinput(N,'prbs'); randn('state',0); e = randn(N,1); %Simulate the process A0 = [1 -1.5 0.7]; B0 = [0 1 0.5]; C0 = [1 -1.8 0.81]; yp = filter(B0,A0,u); v = filter(C0,A0,e); ve = v/std(v)*sqrt(thau)*std(yp); y = yp+ve; z = iddata(y,u); % Estimate an ARX model alpha = 1e3; lambda = 1; [th,yh,p,phi] = rarx(z(1,:),[2 2 1],'ff',lambda,zeros(4,1),alpha.*eye(4)); theta(1,:) = th; cov(:,:,1)=p; yhat(1)=yh; for k = 2:N [th,yh,p,phi] = rarx(z(k,:),[2 2 1],'ff',lambda,th',p,phi); theta(k,:) = th; cov(:,:,k) = p; yhat(k) = yh; end k=0:N-1; figure; plot(k,theta(:,1),'r',k, ones(N,1)*-1.5,'r:'); hold on, plot(k,theta(:,2),'b',k, ones(N,1)*0.7,'b:'); plot(k,theta(:,3),'m',k, ones(N,1),'m:'); plot(k,theta(:,4),'g',k, ones(N,1)*0.5,'g:');

Gambar 2.12 Estimasi nilai parameter model proses setiap iterasi

35

36


Gambar 2.13 Kecenderungan penurunan nilai elemen diagonal matriks kovarian

2.2.3 Efek Harga Awal Untuk mengimplementasikan algoritma estimasi rekursif didalam identifikasi model proses, terdapat dua parameter penalaan estimator yang harus ditentukan pada saat inisialisasi, yaitu harga awal vektor parameter model θ(0) dan matriks kovarian P(0). Penentuan kedua parameter penalaan ini akan sangat terkait dengan knowledge yang dimiliki oleh seorang engineer terhadap proses yang akan diidentifkasi. •

•

Jika seorang engineer tidak memiliki a priori knowledge, yaitu pengetahuan sebelumnya akan proses yang akan dikendalikan baik berupa struktur model proses maupun model proses yang diturunkan berdasarkan pemodelan fisik berdasarkan first principle, maka harga awal dari kedua parameter penalaan tersebut dapat ditentukan sebagai θ(0) = 0 dan P(0) = αI, dengan α nilai cukup besar (α>1000). Nilai α yang besar mengindikasikan bahwa model awal estimator memiliki tingkat ketidakpastian yang tinggi disebabkan ketidaktersediaan pengetahuan sebelumnya akan proses. Jika pemodelan fisik sudah dilakukan sebelum mengimplementasikan algoritma identifikasi, dimana model linier awal diperoleh dengan menurunkan model nonlinier proses dan kemudian disederhanakan melalui proses linierisasi, maka model fisik ini dimasukkan kedalam estimator sebagai model awal proses. Nilai parameter α ditetapkan cukup kecil (α<10), untuk mengurangi perubahan estimasi parameter proses


37

yang drastis, disebabkan karena estimator sudah memiliki model awal yang baik. Gambar 2.14 menunjukkan hasil estimasi simulasi nilai parameter a dan b yang tidak diketahui terhadap harga awal matriks kovarian yang berbeda. Terlihat bahwa untuk nilai α yang besar untuk harga awal parameter θ(0) = 0 akan cepat konvergen untuk mendekati nilai parameter proses sebenarnya.

Gambar 2.14 Pengaruh harga awal P(0) terhadap konvergensi estimasi

2.3 Metode Estimasi Untuk Proses Time-Varying Pada beberapa proses tertentu nilai parameter model berubah terhadap waktu. Perubahan nilai parameter tersebut disebabkan oleh faktor internal atau eksternal dengan waktu. Biasanya penyebabnya diakibatkan linierisasi proses di sekitar titik kerja untuk perubahan sinyal yang kecil. Karena terjadinya perpindahan titik kerja menyebabkan karakteristik non linier mendominasi sifatsifat proses. Jika sifat ketidaklinieran tidak sangat dominan dan perubahan titik kerja lambat, maka karakteristik proses seperti ini dapat direpresentasi oleh persamaan difference dengan parameter yang berubah terhadap waktu.

2.3.1 Pembobotan Eksponensial Dengan Faktor Pelupa Konstan Jika pada metode kuadrat-terkecil standar semua data diperlakukan sama, artinya baik data baru maupun data yang sudah lama mempunyai bobot yang sama. Jika parameter proses berubah, maka data proses yang lama sebenarnya sudah tidak valid lagi, karena data tersebut tidak mewakili karakteristik proses

38


yang baru. Hal ini yang menjelaskan alasan kenapa algoritma kuadrat-terkecil standar tidak mampu menjejaki perubahan nilai parameter proses. Nilai matriks kovarian juga mengecil seiring dengan bertambahnya waktu, dimana tingkat kepercayaan akan model proses bertambah dan tingkat ketidakpastiannya berkurang. Ketika terjadi perubahan nilai parameter proses, nilai elemen matriks kovarian tidak cukup besar untuk memberikan faktor pengkoreksi yang memadai, sehingga tingkat adaptasinya sangat lambat. Prinsip dasarnya adalah dengan menerapkan konsep faktor pelupa, dimana data lama yang mewakili karakteristik sistem sebelum berubah mulai dilupakan, dan lebih menitikberatkan pada data baru yang diukur dari sistem setelah terjadi perubahan karakteristik. Mekanisme pelupa harus secara perlahan-lahan melupakan data lama. Fungsi kriteria pada persamaan (2.12) sekarang ditambah faktor bobot yang bernilai eksponensial untuk setiap data N

N

i =1

i =1

2 J = ∑ λN −i ∈2 (t i ) = ∑ λN −i [y (t i ) − yˆ (t i )]

(2.51)

dengan λ adalah faktor pelupa. Data baru akan mempunyai bobot yang paling tinggi. Nilai bobot itu akan berkurang secara eksponensial seiring dengan bertambahnya waktu. Data lama sedikit demi sedikit akan dilupakan, dan algoritma lebih mencurahkan perhatian kepada data baru. Algoritma estimasi rekursif dengan faktor pelupa ini adalah

[

]

θˆ(t ) = θˆ(t − 1) + γ (t ) y (t ) − ψ (t )θˆ(t − 1)

γ (t ) =

T

P (t − 1)ψ (t )

(2.52)

(2.53)

λ + ψ (t )P (t − 1)ψ (t ) T

P (t ) = [I − γ (t )ψ (t + 1)]P (t − 1)

1

λ

(2.54)

Pengaruh faktor pelupa λ dapat dikenali secara langsung dari persamaan invers matriks kovarian −1

−1

P (t ) = λ P (t − 1) + ψ (t )ψ (t ) T

(2.55)

P-1 proporsional terhadap matriks informasi L L=

1

σ

2 ∈

{

}

EΨ Ψ = T

1

σ

2 ∈

{ }

EP

−1

Faktor pelupa λ harus ditetapkan sebagai berikut:

(2.56)


39

λ kecil (misal λ=0.90), jika kecepatan perubahan parameter besar atau derau tidak terlalu besar. λ besar (misal λ=0.98), jika kecepatan perubahan parameter kecil atau derau besar.

Algoritma estimasi parameter dengan faktor pelupa konstan sangat cocok untuk proses dengan perubahan parameter kecil dan masukan cukup tereksitasi. Jika parameter proses knostan, hasil yang baik dapat diperoleh jika derau yang mempengaruhi proses tidak terlalu besar. Masalah baru timbul jika λ<1 dan masukan tidak cukup tereksitasi. Maka nilai P-1(t) berkurang karena ψ (t ) ≈ 0 , atau elemen P(t) bertambah secara kontinu dan tidak terkendali (blow up). Karena vektor pengkoreksi adalah

γ (t ) = P (t )ψ (t ) maka estimator menjadi sangat sensitif. Gangguan kecil atau kesalahan numerik menyebabkan perubahan estimasi parameter yang tiba-tiba. Estimator menjadi tidak stabil. Situasi ini dapat diamati dengan sistem kendali adaptif. Dengan demikian eksitasi masukan harus dimonitor atau faktor pelupa harus merupakan fungsi waktu.

Contoh 2.5 Diketahui suatu proses yang tidak diketahui mempunyai bentuk model leastsquares dengan bentuk persamaan difference adalah

(1 − 1.5z

−1

)

(

)

+ 0.7 z −2 y (t ) = 1 + 0.5z −1 u(t − 1) + e(t )

dengan variansi derau sebesar 0.2. Setelah pencuplikan ke 200 nilai parameter proses berubah. Bentuk proses yang baru adalah

(1 − 1.0z

−1

)

(

)

+ 0.4z −2 y (t ) = 1 + 0.5z −1 u(t − 1) + e(t )

Buatlah program dalam MATLAB untuk mensimulasikan respon keluaran proses dengan masukan berupa sinyal acak dan untuk menampilkan kurva hasil estimasi dengan menggunakan algoritma RLS dengan faktor pelupa λ=0.9 dan P(0) = 1000.I, dan θ(0) = 0. Solusi

Program simulasi: >> >> >> >> >>

echo off u = sign(randn(200,1)); e = 0.2*randn(200,1); th0 = idpoly([1 -1.5 0.7],[0 1 0.5]); th1 = idpoly([1 -1. 0.4],[0 1 0.5]); y1 = sim(th0,[u e]);

40


>> >> >> >> >> >> >>

y2 = sim(th1,[u e]); y = [y1;y2]; u = [u;u]; z = iddata(y,u); alpha = 1e3; lambda = .9; [th,yh,p,phi] = rarx(z(1,:),[2 2 1],'ff',lambda,zeros(4,1),alpha.*eye(4)); theta(1,:) = th; cov(:,:,1)=p; yhat(1)=yh; for k = 2:400 [th,yh,p,phi] = rarx(z(k,:),[2 2 1],'ff',lambda,th',p,phi); theta(k,:) = th; cov(:,:,k) = p; yhat(k) = yh; >> end

Dari hasil grafik estimasi pada gambar 2.15 terlihat bahwa algoritma RLS dengan faktor pelupa mampu menjejaki perubahan parameter proses a1 dan a2, dimana nilai rata-rata parameter proses yang baru mendekati nilai parameter proses yang sebenarnya. Untuk dapat menjejaki perubahan parameter sistem, maka estimator harus memiliki sifat adaptasi terhadap perubahan karakteristik sistem, dengan cara menetapkan nilai faktor pelupa λ lebih kecil dari 1. Gambar 2.16 menunjukkan perbandingan kurva keluaran proses dengan keluaran prediksi satu-langkah-kedepan. Penyimpangan terjadi pada t=200, namun setelah algoritma RLS dengan faktor pelupa mampu menghasilkan keluaran yang identik dengan keluaran proses.

Gambar 2.15 Hasil estimasi parameter untuk proses yang berubah


41

Gambar 2.16 Respon keluaran proses dengan keluaran prediksi satu-langkah-kedepan

Gambar 2.17 Pengaruh nilai faktor pelupa terhadap kemampuan menjejaki perubahan nilai parameter proses

2.4 Metode Extended Least-Squares (ELS) Permasalahan yang dijumpai dalam identifikasi sistem yang terkontaminasi derau dengan variansi derau yang cukup besar dengan menggunakan algoritma

42


kuadrat-terkecil adalah bahwa estimasi parameter yang dihasilkan bias dan tidak konsisten. Hal ini terutama disebabkan oleh fakta bahwa model ARX yang digunakan dalam algoritma kuadrat-terkecil mengasumsikan model derau yang sangat khusus yang sangat jarang dijumpai dalam sistem riil, seperti terlihat pada gambar 2.18 di bawah ini.

ARX model e(t) u(t)

1

B(z-1)

y(t)

A(z-1)

Real process v(t) u(t) B(z-1)

1

y(t)

A(z-1)

Gambar 2.18 Model ARX mengasumsikan model derau 1/A, sementara keluaran proses riil terkontaminasi oleh derau putih atau derau warna

Bias pada estimasi parameter mempunyai arti bahwa parameter hasil estimasi terdeviasi dari nilai optimal atau nilai sebenarnya. Ketidakkonsistenan algoritma kuadrat-terkecil mengancung arti bahwa nilai bias tersebut tidak mendekati nilai nol walaupun jumlah data yang digunakan mendekati tak hingga. Gambar 2.19 menggambarkan ketidakkonsisten metode kuadrat terkecil ketika jumlah data yang lebih besar tidak memperbaiki kualitas estimasi model. Karena konsistensi merupakan sifat estimator utama yang harus dimiliki, beberapa strategi telah dikembangkan untuk menghindari ketidakkonsistensi estimasi untuk model ARX. Ide dasar yang melatarbelakangi sebagian besar dari pendekatan tersebut adalah dengan menjaga sifat linieritas pada parameter model ARX, karena hal ini merupakan keuntungan terbesar yang dapat membuat solusi dari permasalahan identifikasi dapat berupa solusi analitik. Salah satu alternatif pemecahan permasalahan di atas adalah dengan memilih struktur model yang lebih umum seperti model ARMAX yang mempunyai hubungan tidak linier antar parameter, dan mengembangkan algoritma yang dapat mengestimasi parameter tidak linier dengan teknik metode kuadrat-terkecil linier bertingkat. Dalam estimasi dengan model least-squares diasumsikan proses diwakili oleh persamaan A( z −1 )y (t ) = B( z −1 )u(t − 1) + e(t )

(2.57)


43

Jika proses terkontaminasi oleh derau warna model yang tepat digunakan adalah model ARMAX A( z −1 )y (t ) = B( z −1 )u(t − 1) + C( z −1 )e(t )

(2.58)

dengan polinomial C(z-1) secara umum tidak dapat diasumsikan sama dengan satu. Jika proses dengan derau warna diestimasi dengan algoritma leastsquares standar, maka estimator akan memberikan hasil yang bias.

Gambar 2.19 Sifat unconsistent dari algoritma kuadrat terkecil dalam estimasi parameter model dengan keluaran proses terkontaminasi oleh derau

Untuk mengestimasi polinomial C(z-1) diperlukan data derau e(t-1), ..., e(t-nc). Akan tetapi {e(t)} adalah derau putih yang tidak dapat terukur. Hai yang terbaik yang dapat dilakukan adalah menggunakan kesalahan prediksi ∈(t) untuk menggantikan data derau yang tidak dapat diperoleh. Bentuk persamaan prediktor optimal untuk model ARMAX adalah yˆ (t | t − 1) =

⎛ A( z −1 ) ⎞ B( z −1 ) ⎜⎜1 − ⎟ y (t ) + u ( t ) −1 ⎟ C ( z −1 ) ⎝ C( z ) ⎠

(2.59)

Formula prediktor ARMAX akan selalu stabil walaupun polinom A(z-1) mempunyai akar yang tidak stabil. Hanya saja, polinom C(z-1) harus dijaga untuk selalu stabil. Struktur galat prediksi untuk model ARMAX ini seperti terlihat pada gambar 2.20, mempunyai persamaan

44


∈ (t ) =

A( z −1 ) B( z −1 ) y ( t ) − u (t ) C( z −1 ) C( z −1 )

(2.60)

process

e(t) CP(z-1) AP(z-1)

u(t)

AP

v(t)

+

BP(z-1)

y(t)

(z-1)

1

1

C(z-1)

C(z-1) B(z-1)

-

+

A(z-1)

ε(t)

Gambar 2.20 Model ARMAX dalam konfigurasi persamaan kesalahan

Perbedaan dengan skema estimasi pada model ARX adalah bahwa data masukan-keluaran terukur difilter terlebih dahulu dengan polinom 1/C(z-1). Kesalahan prediksi dapat diekspresikan dalam bentuk C( z −1 ) ∈ (t ) = ∈ (t )

A( z −1 )y (t ) − B( z −1 )u(t )

(

)

= A( z −1 )y (t ) − B( z −1 )u(t ) + 1 − C( z −1 ) ∈ (t )

(2.61)

Dan jika diubah kedalam domain waktu menghasilkan persamaan diferens seperti persamaan berikut ∈ (t ) = y (t ) + a1y (t − 1) + L + ana y (t − na ) + L − b0 u(t − 1) − L − bnb u(t − nb ) − L

(2.62)

− c1 ∈ (t − 1) − L − c nc ∈ (t − nc ) 144444 42444444 3 − yˆ ( t |t −1)

Sinyal ∈(t-i) bukan hasil pengukuran melainkan harus dikalkulasi dari residu sebelumnya dan estimasi sinyal derau e(t-i) yang tidak diketahui. Permasalahan identifikasi model ARMAX berbasiskan kepada persamaan residu 2.62. Ada dua pendekatan berbeda untuk solusi permasalahan tersebut. Pendekatan pertama berdasarkan kepada permasalahan teknik optimasi nonlinier multivariabel. Beberapa metode yang dikenal dalam teknik optimasi seperti metode Newton, Levenberg-Marquardt, conjugate direction, Quasi-Newton, dapat digunakan untuk memecahkan permasalahan tersebut. Hanya saja waktu komputasi yang dibutuhkan sangat membebankan perhitungan terutama untuk aplikasi real time seperti sistem kendali swa-tala.


45

Pendekatan kedua dengan memanfaatkan sifat yang dimiliki oleh algoritma kuadrat terkecil namun digunakan untuk mengestimasi parameter model nonlinier. Teknik ini dikenal sebagai metode kuadrat terkecil bertingkat. Dari aspek kecepatan konvergensi, pendekatan kuadrat terkecil bertingkat lebih cepat dibandingkan teknik optimasi nonlinier. Hanya saja, solusi yang diperoleh dapat terjebak kedalam permasalahan lokal minima tidak dapat dipecahkan. Langkah-langkah yang diperlukan penyelesaian identifikasi sistem dengan menggunakan algoritma extended least-squares adalah sebagai berikut: Langkah (i) Estimasi model ARX A(z-1)y(t) = B(z-1)u(t-1) + v(t) berdasarkan data pengukuran {u(t), y(t)} dengan menggunakan persamaan

(

)

−1

θ ARX = Φ Φ Φ y T

T

Langkah (ii) Kalkulasi kesalahan prediksi model ARX ∈ARX(t) = A(z-1)y(t) – B(z-1)u(t) dimana B(z-1) dan A(z-1) ditentukan dari θARX. Langkah (iii) Estimasi parameter model ARMAX ai, bi, dan ci dari persamaan 2.58 dengan algoritma kuadrat terkecil dengan pendekatan residu model ARMAX ∈(t-i) ≈∈ARX(t-i). Langkah (iv) Langkah (ii) dan (iii) dapat diiterasi hingga konvergensi estimasi tercapai. Dengan iterasi perbaikan kualitas estimasi parameter model, membuat algoritma ELS sangat cocok digunakan dalam identifikasi sistem dengan derau yang cukup besar. Kesalahan prediksi ARMAX akan mendekati karakteristik derau putih, yaitu uncorrelated dan rata-rata nol, seiring dengan semua informasi sistem yang diperlukan dapat dieksploitasi oleh model dengan baik.

Algoritma Recursive Extended Least-Squares (RELS) Algoritma rekursif ELS sama dengan algoritma least-squares standar yaitu persamaan (2.47), (2.48), dan (2.49), yang membedakan hanya bentuk vektor parameter θˆ dan vektor data ψT(t) T

[

θˆ = aˆ1 , L , aˆ na , bˆ1 , L , bˆnb , cˆ1 , L , cˆ nc dan

]

(2.63)

46


ψ (t ) = [− y (t − 1),L,− y (t − na ),u(t − 1),L, u(t − nb − 1),∈ (t − 1),L,∈ (t − nc )] (2.64) T

Nilai sinyal kesalahan prediksi ∈(t) yang juga disebut a priori error didalam vektor data ψT(t) dihitung rekursif dengan menggunakan persamaan (2.51). Oleh karena itu akar-akar dari polinomial C(z-1) harus terletak didalam unit lingkaran bidang z. Untuk mempercepat konvergensi estimasi parameter, umumnya perhitungan estimasi menggunakan a posteriori error, dimana kesalahan prediksi berbasiskan kepada nilai estimasi parameter saat sekarang, dibandingkan menggunakan a priori error yang berdasarkan nilai estimasi parameter sebelumnya. Sinyal a posteriori error mempunyai bentuk persamaan ~ (t ) = y (t ) − ϕ T (t )θˆ(t ) ∈

(2.64)

Dalam perhitungan algoritma RELS vektor data yang digunakan sekarang berubah menjadi ~ (t − 1),L,∈ ~ (t − nc )] ψ (t ) = [− y (t − 1),L,− y (t − na ),u(t − 1),L,u(t − nb − 1),∈ T

Dibandingkan metode RLS, konvergensi estimasi algoritma RELS lebih lambat, khususnya pada estimasi parameter model derau ci. Hal tersebut disebabkan oleh fakta bahwa selama prosedur estimasi parameter diperlukan pendekatan sinyal derau putih. Terutama pada awal iterasi, dimana pendekatan derau putih sangat buruk yang mengakibatkan buruknya kualitas estimasi parameter ai dan bi. Efeknya dapat menurunkan konvergensi estimasi. Langkah-langkah yang diperlukan algoritma RELS adalah, pada setiap waktu diskrit t: Langkah (i) Bentuk vektor data ϕT(t) menggunakan data baru terukur Langkah (ii) Hitung kesalahan prediksi ∈(t) menggunakan ∈ (t ) = y (t ) − ϕ (t )θˆ(t − 1) T

Langkah (iii) Perbaharui matriks kovarian P(t) menggunakan T ⎡ ϕ (t )ϕ (t )P (t − 1) ⎤ P (t ) = P (t − 1)⎢I − ⎥ T ⎢⎣ λ + ϕ (t ) P (t − 1)ϕ (t ) ⎥⎦

Langkah (iv) Perbaharui vektor parameter θ(t)

θˆ(t ) = θˆ(t − 1) + P (t )ϕ (t ) ∈ (t )


47

Langkah (v) Hitung sinyal a posteriori error menggunakan ~ (t ) = y (t ) − ϕ T (t )θˆ(t ) ∈ Langkah (vi) Tunggu periode waktu pencuplikan berikutnya dan ulangi langkah (i) Untuk contoh 2.4 dengan variansi derau σ ∈2 =1, hasil estimasi yang diperoleh algoritma ELS juga lebih baik dibandingkan dengan hasil dari algoritma kuadratterkecil standar seperti terlihat pada gambar 2.10. Nilai rata-rata estimasi parameter proses berada di sekitar nilai parameter proses yang sebenarnya. Program simulasi: >> >> >> >> >> >> >> >>

v=1; u = sign(randn(200,1)); e = v*randn(200,1); th0 = idpoly([1 -1.9 0.9],[0 0 0.1],[1 -1.2 0.6]); y = sim(th0,[u e]); Data = iddata(y,u); na=2; nb=1; nc=2; nk=2; alpha = 1e3; [th1,yh1,P1] = rarx([y u],[na nb nk],'ff',1,zeros(na+nb,1), ... alpha*eye(na+nb)); >> [th2,yh2,P2] = rarmax([y u],[na nb nc nk],'ff',1,zeros(na+nb+nc,1),... alpha*eye(na+nb+nc));

Gambar 2.21 Perbandingan respons step proses dengan model

48


Gambar 2.22 Hasil estimasi parameter proses dengan RLS dan RELS

2.5 Metode Instrumental Variable Pendekatan algoritma kuadrat-terkecil bertingkat seperti algoritma ELS memberikan keuntungan dalam estimasi paramater model nonlinier dari aspek waktu komputasi yang lebih singkat dibandingkan jika harus diselesaikan dengan teknik optimasi. Namun metode dikhususkan untuk estimasi model ARMAX yang mempunyai hubungan nonlinier diantara parameter model. Metode lain yang sangat populer dan mudah digunakan untuk mengatasi permasalahan konsistensi dari estimasi model ARX adalah metode instrumental variable (IV). Ide yang mendasari metode ini adalah dengan mengalikan vektor kesalahan prediksi dengan suatu matriks Z yang mempunyai dimensi yang sama dengan matriks regressor Φ. Setiap kolom matriks Z disebut sebagai variabel instrumen yang ditentukan oleh pemakai agar tidak berkorelasi dengan derau, sehingga jika semua informasi sistem dapat dieksploitasi model dengan baik, variabel instrumen juga tidak berkorelasi dengan ∈(t). Hal ini berarti bahwa setiap baris dari matriks Z orthogonal terhadap vektor kesalahan prediksi ∈, ZT∈=0. Dengan mengalikan persamaan vektor kesalahan prediksi dengan ZT menghasilkan 0 = ZTy - ZTΦθ

(2.65)


49

dan konsekuensinya ZTy = ZTΦθ

(2.66)

Jika matriks ZT non-singular, yang mewakili sifat masukan yang cukup tereksitasi dan pemilihan variabel instrumen yang representatif, algoritma estimasi instrumental variable dinyatakan dalam persamaan

(

)

−1

θˆ = Z Φ Z y T

T

(2.67)

Terlihat bahwa algoritma IV memiliki persamaan yang sederhana seperti yang dimiliki oleh algoritma kuadrat terkecil. Yang menjadi permasalahan dalam penggunaan metode IV adalah penentuan variabel instrumen agar dapat memperbaiki sifat konsistensi. Yang menyebabkan terjadinya bias estimasi pada algoritma kuadrat terkecil adalah penggunaan data keluaran terukur pada matriks regressor. Keluaran sistem sudah berkorelasi dengan derau. Hal ini akan membuat matriks regressor menjadi matriks singular. Salah satu cara untuk menentukan variabel instrumen adalah dengan menggunakan data keluaran model ARX yu(t) dalam skema output error dibandingkan data keluaran terukur y(t). Data keluaran model menggantikan data keluaran sistem terukur pada matriks regressor sementara data masukan proses u(t) dibuat tetap. Matriks baru ini yang disebut sebagai matriks Z. Algoritma lengkap dari metode instrumental variable dapat dilihat di bawah ini. Langkah (i) Estimasi model ARX berdasarkan data pengukuran {u(t),y(t)} menggunakan persamaan

(

)

−1

θˆ ARX = Φ Φ Φ y T

T

Langkah (ii) Simulasikan model ini untuk menghasilkan keluaran model yu(t) y u (t ) =

Bˆ ( z −1 ) u(t − 1) Aˆ ( z −1 )

dengan B(z-1) dan A(z-1) ditentukan oleh θARX Langkah (iii) Bentuk variabel instrumen ⎡ − y u (n ) ⎢ − y (n + 1) Z=⎢ u ⎢ M ⎢ ⎣− y u (N − 1)

− y u (1)

u(1) ⎤ L − y u ( 2) u(n + 1) L u(2) ⎥⎥ ⎥ O M M O M ⎥ L − y u (N − n ) u(N − 1) L u(N − n )⎦ L

u(n )

L

50


Langkah (iv) Estimasi parameter model dengan metode IV

(

)

−1

θˆ IV = Z Φ Z y T

T

Langkah (v) Iterasi langkah 2-4 algoritma IV hingga konvergensi tercapai

Gambar 2.23 Perbandingan kemiripan model dengan keluaran proses yang dihasilkan metode LS dan IV

Gambar 2.24. Sifat konsistensi metode IV terhadap bertambahnya jumlah data


51

Untuk contoh 2.4 dengan variansi derau σ ∈2 =1, hasil simulasi model dengan masukan berupa fungsi step yang diperoleh algoritma IV juga lebih baik dibandingkan dengan hasil dari algoritma kuadrat-terkecil standar seperti terlihat pada gambar 2.23. Hingga iterasi ke-5, metode IV dapat membuat kemiripan yang lebih baik antara model dengan keluaran proses. Sifat konsistensi metode IV juga lebih baik dibandingkan metode kuadrat-terkecil. Seperti terlihat pada gambar 2.24, metode IV dapat memperbaiki kemiripan model dengan proses untuk perolehan data estimasi yang lebih banyak.

Metode Recursive Instrumental Variable (RIV) Algoritma RLS secara umum menghasilkan estimasi parameter yang tidak konsisten dari nilai parameter proses sebenarnya seperti yang dihasilkan oleh algoritma kuadrat terkecil. Seperti yang sudah dijelaskan pada bagian sebelumnya, salah satu metode yang dapat mengatasi permasalahan ketidakkonsistenan adalah metode instrumental variable. Dengan adanya variabel instrumen untuk menghindari matriks regressor agar tidak menjadi singular, maka bentuk rekursif dari algoritma instrumental variable adalah

θˆ(t ) = γ (t ) = P (t ) =

θˆ(t ) + γ (t )e(t ) 1

ϕ (t )P (t − 1)z(t ) + λ T

P (t − 1)z(t )

(2.68)

(I − γ (t )ϕ (t ))P(t − 1) λ 1

T

Variabel instrumen ditentukan oleh keluaran model agar bebas linier terhadap derau z(t ) = [− y u (t − 1) L − y u (t − n ) u(t − 1) L u(t − n )]

(2.69)

dengan

[

]

y u (t ) = − aˆ1y u (t − 1) L − aˆ na y u (t − na ) + bˆ0 u(t − 1) L bˆnb u(t − nb ) (2.70)

Penggunaan model bebas IV dapat menyebabkan permasalahan kestabilan karena model ini berbasis kepada estimasi parameter yang merepresentasikan lup dalam lup. Permasalahan yang dijumpai ini menyebabkan aplikasi metode RIV tidak terlalu sukses dalam praktis.

52


Penggunaan model bebas ini juga membuat selama beberapa kali iterasi pertama estimasi algoritma RIV, parameter terestimasi ai dan bi tidak akurat yang disebabkan model bebas yang dimiliki masih memiliki unsur ketidakpastian yang tinggi. Untuk mengatasinya, algoritma RLS digunakan selama fasa awal estimasi. Algoritma RLS dijalankan selama beberapa periode pencuplikan, misalkan 4n, hingga estimasi parameter terlihat sudah stabil.

Gambar 2.26 Hasil perbandingan estimasi model dengan RLS dan RIV

Model Sistem dan Metode Estimasi Parameter Rekursif

Recommend Documents