TUGAS BESAR KONTROL HIBRID ROBUST “Pengontrolan dengan Metoda Linear Quadratic Regulator (LQR) dan Metoda Algebraic Riccati Equation (ARE) untuk Mengurangi Resonansi Ketika Despin pada Dual-Spin Spacecraft”
Disusun Oleh: Dita Ayu Banjarnahor/ 13312032 Qanun Miladial Hikmah/ 13312088
Dosen Pengajar: Dr. Ir. Endra Joelianto
PRODI TEKNIK FISIKA FAKULTAS TEKNOLOGI INDUSTRI INSTITUT TEKNOLOGI BANDUNG 2015
I. Deskripsi Sistem DUAL-SPIN SPACECRAFT Dual Spin Spacecraft merupakan satelit yang di gunakan untuk menangkap citra luar angkasa. Pada kondisi ruang angkasa banyak terdapat gangguan yang tak
terduga
yang
menyebabkan
keseimbangan dari satelit terganggu, salah satunya ketika Preccesion Phase Lock (PPL). Pada saat penangkapan PPL, gerakan dari ketidakseimbangan motor
Gambar 1.1. Satelit
disinkronkan dengan presisi beban inersia SIC. Ketika despin dari motor konstan dan relatif kecil terhadap ketidakseimbangan, maka PPL menyebabkan deviasi sekuler yang besar pada sumbu poros dari sumbu yang diinginkan. Untuk mengatasi hal tersebut, diterapkan konstan
torque
pada
kapasitas
maksimum
motor
serta
dikembangkan alternatif solusi berdasarkan non-linier closed loop feedback control dari torsi motor despin. SISTEM DUAL-SPIN SPACECRAFT Dual spin spacecraft terdiri dari 2 buah badan yang terhubung pada satu poros yang menyebabkan terjadinya rotasi relatif. Badan pertama disebut sebagai rotor yang berputar di sekitar sumbu axis relatif rotasi. Rotor dilengkapi dengan stiffness gyroscope yang akan menjaga gerak dari rotor ketika berputar. Badan ke dua disebut sebagai platform, yang tidak mengalami rotasi. Hal ini berfungsi sebagai acuan arah dari muatan satelit. Gambar 1.2 Dual Spin spacecraft NON-LINIER CONTROLLER Pada proses PPL, telah diketahui bahwa satelit dapat mengalami ketidakseimbangan sehingga diperlukan pengontrol untuk mengatur mengurangi resonasi pada sistem. Untuk
memudahkan pengembangan sistem kontrol pada dual-spin spacecraft, digunakan pendekatan sistem mekanik
yang
digunakan spring-mass
lebih
pendekatan oscillator
sederhana. sistem dengan
Sehingga
menggunakan embedded
unbalance rotor atau yang lebih dikenal dengan sistem TORA (Translation Oscillations with an Gambar 1.3 Sistem TORA
Eccentric Rotational Proof Mass Actuator). TORA
adalah troli dengan beban yang dihubungkan pada pegas dengan konstanta pegas k. Troli memiliki gerak satu dimensi kearah horizontal. F merupakan gangguan (gaya luar) yang diberikan kepada sistem. Pada troli terdapat massa yang berotasi dengan momen inersia I dan jarak e dari posisi putarnya dan pusat massa beban. Pengontrol Torsi N diberikan pada beban. Ө menyatakan sudut, ketika Ө=00 maka bebannya tegak lurus dari posisi troli sedangkan Ө=900 maka searah dengan q nya. Persamaan untuk sistem di atas adalah :
√
X(t) = ( x x(dot) Ө Ө(dot) )
Kondisi awal yang diketahui Tabel 1.1. Parameter sistem dan nilainya masing-masing No
Besaran
Parameter
Nilai
Satuan
1
Cart Mass
M
1.3
Kg
2
Arm Mass
m
0.1
Kg
3
Arm Eccentricity
e
0.06
M
4
Arm Inertia
I
0.0002
Kg m2
5
Spring Stiffness
K
186.3
N/m
6
epsilon
Ԑ
0.2
7
Torque
N
1.3
mN
8
Simpangan
q
0.1
M
Kondisi awal yang diketahui dimasukan ke dalam persamaan dalam matriks sehingga didapat persamaan berikut (dihitung dengan menggunakan Matlab) %kondisi awal yang diketahui q= 0.1; M=1.3; m=0.1; e=0.06; I=0.0002; k=186.3; epsilon= 0.2; N=1.3; pi=3.14; u=((M+m)/(k*(I+m*e^2)))*N;
x(1) = square((M+m)/(I+m*e^2))*q; x(2)=0; x(3)=-20/(180*pi); x(4)=-epsilon*x(2)*cos(x(3))+u; x=[ x(1); x(2); x(3) ; x(4) ]
%matriks A= matriks f(x(t)) A = [0 1 0 0; -1/(1-epsilon^2*cos(x(3))^2) 0 0 (epsilon*x(4)*sin(x(3)))/(1epsilon^2*cos(x(3))^2); 0 0 0 1 (epsilon*cos(x(3)))/(1-epsilon^2*cos(x(3))^2) 0 0 (x(4)*epsilon*sin(x(3)))/(1-epsilon^2*cos(x(3))^2)] %matriks B= matriks g(x(t)) B = [ 0 ; -epsilon*cos(x(3))/(1-epsilon^2*cos(x(3))^2) ; 0 ; 1/(epsilon^2*cos(x(3)^2))] C = [0 0 1 0] D = [0]
Didapat hasil matriks berikut:
Sedangkan fungsi transfer sistem TORA adalah -3.053e-016 s^3 + 25 s^2 - 1.055e-015 s + 26 G(s) = -----------------------------------------------------------s^4 - 0.1286 s^3 + 1.042 s^2 - 0.1072 s
Hasil plot respon sistem pada domain waktu terlihat pada gambar berikut.
Gambar 1.4 Respon sistem awal terhadap waktu Terlihat bahwa sistem tidak terkendali dan menunjukan nilai yang tak hingga. Oleh karena itu, diperlukan pengontrol untuk membuatnya stabil.
Sedangkan pada domain frekuensi, respon sistem ditunjukkan oleh gambar berikut.
Gambar 1.5 Respon sistem awal dalam domain frekuensi Terlihat bahwa sistem hanya bekerja pada frekuensi rendah. Akibatnya, sisem cenderung dipengaruhi oleh disturbance. Analisis kestabilan pada bode plot dengan teorema Small Gain Theorem menyebutkan bahwa pada saat fasa dari sistem berada pada -180 dan frekuensi = 1 Hz maka sistem akan stabil ketika magnitudanya kurang dari 1 dB. Berdasarkan plot frekuensi sistem, terlihat bahwa ketika fasanya -180 pada saat frekuensi = 1 Hz, besar magnituda dari sistem memiliki harga lebih besar dari 1 dB. Oleh karena itu, terlihat bahwa sistem tidak stabil. Sedangkan, dengan meninjau kestabilan berdasarkan posisi dan bentuk root locus sistem, diperoleh plot root locus sebagai berikut.
Gambar 1.6 Plot root locus Analisis kestabilan berdasarkan pole dari sistem menyebutkan bahwa sistem akan stabil ketika polenya berada pada sumbu kiri; artinya besar pole < 0. Jika dilihat pada gambar di atas, ditunjukkan bahwa sistem sangat tidak stabil karena seluruh polenya berada pada sumbu kanan atau pole > 0. Maka, untuk membuat sistem menjadi stabil, diperlukan langkah-langkah pengontrolan. Sebelum menentukan jenis pengontrolan yang akan dilakukan, terlebih dahulu dilakukan pengecekan terhadap keterkontrolan dan keteramatan sistem.
II. Cek Controllability dan Observability Controllability merupakan parameter untuk mengetahui apakah sistem yang dirancang dapat dikontrol atau tidak. Artinya, seberapa besar pengaruh pemberian input terhadap perubahan variabel keadaan. Besar controllability akan ditentukan oleh matriks A dan B. Observability diperlukan untuk mengetahui apakah sistem dapat ditinjau hasilnya. Besar observability ditentukan oleh matriks A dan C. Nilai observability dan controllability dihitung dengan menggunakan program Matlab berikut: %cek CTRB CB=ctrb(A,B) rank (CB) [Mcb,Ncb]=size(CB); OB=obsv(A,C) rank (OB) [Mob,Nob]=size(OB); if rank(CB)==Mcb disp('SISTEM CONTROLLABLE'); else disp('SISTEM NON CONTROLLABLE'); end if rank (OB)==Mob disp ('SISTEM OBSERVABLE'); else disp ('SISTEM TIDAK OBSERVABLE') end
Sistem akan controllable ketika Rank (CB) sama dengan size (CB) atau ukuran matriksnya. Sedangkan sistem observable ketika Rank(OB) sama dengan size (OB) atau ukuran matriksnya. Berdasarkan matriks yang kita punya, besar rank(CB) dan size(CB) sama yaitu 4
sehingga sistem CONTROLLABLE. Sedangkan rank(OB) dan size (OB) juga sama yaitu 4 sehingga system OBSERVABLE. Karena sistem terkontrol dan teramati, selanjutnya ditentukan metoda pengontrolan yang sesuai.
III. Metoda Linear Quadratic Regulator (LQR) Untuk menentukan kontroler yang optimal, dengan Matlab dihitung dengan menyelesaikan continuous-time, quadratic regulator problem, dan the Associated Riccati Equation. Persamaan Riccati sendiri diberikan oleh persamaan berikut.
Dari perintah tersebut akan muncul nilai K, P, dan E. Nilai K merupakan optimal feedback gain. Nilai P adalah matriks positive-definite dari solusi persamaan Riccati. Sedangkan nilai E adalah eigen value dari A-B*K. K= -1.0178 -0.9743
1.0000
1.0414
P= 16.9479 -0.0254 -0.9743 -0.0409 -0.0254 15.3366
0.9867
0.0887
-0.9743
0.9867
1.1621
0.0482
-0.0409
0.0887
0.0482
0.0424
E= -24.9808 -1.0008 -0.0639 + 1.0177i -0.0639 - 1.0177i
Untuk melihat respon dari sistem yang dirancang, ditinjau plot dari masing-masing variabel yang dicari. Berikut program Matlab untuk plotting %pembuatan plot K = [Khat(1) Khat(2) Khat(3) Khat(4)]; P; AA = [ A-B*K B*Ki ; -C 0 ]; BB = [ 0;0;0;0;1 ]; CC = [ C 0 ]; DD = [D]; t=0:0.02:15; [x,y,t]= step(AA,BB,CC,DD,1,t); plot (t,y) grid legend('x1','x2','x3','x4','input')
Dari plot di atas, besar Q dan R sistem pada proses penentuan KPE diubah-ubah untuk mengetahui pengaruh dari perubahan Q dan R. Q merupakan matriks identitas yang dapat diubah nilainya untuk memperkecil nilai x. Sedangkan R adalah quadratic performance index. Dilakukan manipulasi harga Q dan R sehingga didapatkan plot berikut.
Gambar 3.1 Respon sistem terhadap waktu dengan variasi nilai R pada Q = 1 Jika ditinjau dari masing-masing grafik, terlihat bahwa ketika nilai R kecil, grafik cenderung menunjukkan respon yang menyimpang, di mana pada R=0.0001 magnitudanya
bernilai negatif. Hingga nilai R=10000 grafik menjadi sangat linier. Grafik tidak berosilasi terlebih dahulu tetapi langsung mencapai nilai steadynya, dengan settling time di sekiar t=7 detik. Tetapi jika kita tinjau dari nilai integrator gainnya, didapat bahwa besar integral gainnya berbanding terbalik dengan nilai R. Semakin besar R, maka nilai Ki semakin kecil. Tetapi untuk R=1, nilai Ki=1; dan ketika R diperkecil, nilai Ki menjadi sangat besar. Hubungan Ki dan R sendiri ditunjukkan oleh persamaan berikut Ki=1/square(R).
Gambar 3.2 Respon sistem terhadap waktu dengan variasi nilai Q pada R = 1 Ketika nilai Q diubah-ubah, terlihat perbedaan pada grafik yang didapat. Saat nilai Q=0.0001 grafik menunjukan nilai yang langsung steady sama seperti ketika saat R=10000, sedangkan ketika nilai Q besar, grafik cenderung sama dengan ketika R kecil. Berdasarkan teori, besar Q berbanding terbalik dengan R dan hal ini terbukti pada grafik yang didapat. Maka dari itu, dipilihlah pengontrol dengan harga Q yang kecil dan R yang besar sehingga didapatkan harga K, P dan E adalah K= 0.0166 -0.0013
0.0001
0.0112
0.0079 -0.0013
6.6230
P= 54.4645
0.0079 52.4754 -0.0079 -0.0660 -0.0013 -0.0079
0.0050
0.0399
6.6230 -0.0660
0.0399
4.4609
E= -0.0127 + 1.0192i -0.0127 - 1.0192i -0.1001 -0.0250 Sehingga didapatkan respon dan pengontrol sistem sebagai berikut.
Gambar 3.3 Respon frekuensi sistem pengontrol
Gambar 3.4 Plot pengontrol LQR Berdasarkan hasil plot gambar di atas, terlihat bahwa sistem tidak mengalami osilasi dan cenderung dapat terkontrol. Hasil ini dirasa belum memuaskan, sebab terlihat bahwa cukup lama waktu yang dibutuhkan untuk mencapai steady state. Namun, dari sisi kestabilan, analisis kestabilan sistem pada domain frekuensi menunjukkan bahwa pada saat fasanya 180 dan frekensi 1 Hz, magnituda sistem kurang dari 1 dB, yang artinya sistem telah stabil. Dengan meninjau posisi dan bentuk root locus, diperoleh plot berikut.
Gambar 3.5 Plot root locus
Berdasarkan analisis pole pada sistem di atas, setelah diberikan pengontrol LQR, ditunjukkan bahwa pole berada pada sumbu sebelah kiri, artinya pole < 0, yang artinya sistem stabil.
IV. Analisis Respon Frekuensi
Gambar 4.1 Sensitivity bode plot
Gambar 4.2 Open loop bode plot
Gambar 4.3 Tracking bode plot Grafik sensitivity menunjukkan karakter yang mirip dengan High Pass Filter. Artinya, sinyal-sinyal pada frekuensi rendah akan difilter, sedangkan sinyal pada frekuensi lebih tinggi akan diteruskan sebagaimana adanya. Sehingga, disturbance yang biasa terdapat pada frekuensi rendah, dapat diminamilisir. Dapat dikatakan, grafik sensitivity ini menunjukkan respon sistem yang cukup baik. Grafik tracking pada frekuensi rendah meneruskan sinyal yang diterimanya. Pada frekuensi 1rad/s ke atas, mulai menunjukkan karakter memperkecil amplitude sinyal yang diterimanya. Sehingga, noise yang biasanya ada pada frekuensi tinggi, dapat diminimalisir. Grafik open loop gain menunjukkan profil yang cukup curam. Hal ini sesuai dengan yang diharapkan, di mana pada frekuensi yang semakin tinggi, noise akan diperkecil nilainya sampai menuju nol. Sedangkan pada frekuensi rendah, sinyal aslinya akan tetap dipertahankan, bahkan cenderung diamplifikasi.
V. Metoda PID pada LQR
Pada PID LQR ini, digunakan harga Q yang dimanipulasi agar didapat hasil yang optimal. Dipilih Q= 5 eye(5) serta R=1. Berdasarkan harga tersebut, didapatkan parameter KPE sebagai berikut. Kpid = -1.8195 -2.4948
2.2361
2.7663 12.0147
Ppid = 85.4165
0.6119 -5.5784 -0.9963 -1.8195
0.6119 77.5503
4.1276 -0.4804 -2.4948
-5.5784
4.1276
6.5737
1.1090
2.2361
-0.9963 -0.4804
1.1090
1.2090
2.7663
-1.8195 -2.4948
2.2361
2.7663 12.0147
Epid = -5.3789 + 5.1899i -5.3789 - 5.1899i -1.0006 -0.0639 + 1.0177i -0.0639 - 1.0177i Dari harga K yang didapat, dilakukan plot pada pengontrol dan inputnya. Dari gambar di bawah, terlihat bahwa pengontrol mengalami osilasi terlebih dahulu sebelum mencapai kestabilan.
Gambar 5.1 Plot pengontrol PID LQR terhadap waktu Pada metoda PID LQR, agar didapatkan 3 parameter K, maka matriks input diaugmentasi terlebih dahulu sehingga didapatkan matriks orde 5. K hasil LQR yang berjumlah 5 akan dikalikan dengan pinv(b) dimana b = [C 0; C*A C*B; C*A*A C*A*B] Hal inilah yang membedakan metoda PID LQR dengan metoda LQR. Sehingga didapatkan parameter PID sebesar Kd = 0.4799
Kp = 2.7046
Ki =2.2361
Nilai-nilai parameter PID tersebut kemudian dimasukkan ke program simulasi Simulink berikut.
Gambar 5.1 Simulasi PID Maka didapatkan hasil seperti pada gambar berikut.
Gambar 5.2 Simulasi Simulink respon sistem terhadap waktu dengan metoda PID LQR Terlihat bahwa sebelum diberi pengontrol, sistem menuju nilai tak hingga sehingga tidak stabil. Tetapi ketika diberikan pengontrol, sistem cenderung lebih stabil. Berdasarkan gambar, terlihat bahwa respon sitem mengalami overshoot sebesar 24% tetapi cepat mencapai keadaan stabilnya. Nilai overshoot yang cukup besar ini menunjukkan hasil yang belum meuaskan. Namun, besar settling time hanya sekitar 1.7 s sehingga sistem relatif cepat stabil. Ini menunjukkan respon sistem yang cukup terkontrol. PERBANDINGAN METODA LQR DAN PID LQR
Gambar 5.3 Perbandingan respon sistem dengan metoda LQR dan PID LQR untuk input step Respon sistem pada gambar di atas menunjukkan perbedaan yang cukup besar antara metoda pengontrol LQR dengan PID LQR. Pengontrol PID LQR menunjukkan respon sistem yang lebih cepat mencapai set point, namun berosilasi terlebih dulu sebelum mencapai keadaan tunak. Osilasi yang dihasilkan memiliki maximum peak yang sangat rendah, jauh di bawah overshoot 10%. Artinya, nilai tersebut masih dapat ditolerir pada plant secara umum. Sedangkan pengontrol LQR memberikan respon sistem yang langsung mencapai keadaan tunak tanpa mengalami overshoot. Namun, waktu yang dibutuhkan untuk mencapai set point lebih lama dibandingkan metoda PID LQR. Hasil yang ditunjukkan oleh metoda PID LQR lebih dapat diterima karena lebih cepat mencapai set point dengan osilasi yang dapat diabaikan. Hal ini
dikarenakan, pada metoda ini telah diperhitungkan parameter-parameter pengontrol PID, yaitu Kp, Ki, dan Kd yang akan memberikan hasil yang lebih optimal.
VI.
Metoda PID H-∞ (ARE) ARE (algebraic Riccati equation) adalah persamaan non linier pada sistem H tak hingga
untuk control optimal pada sistem kontinu atau diskrit. Persamaan ARE untuk sistem kontinu:
Di mana nilai B tergantung pada gamma dengan persamaan berikut:
B1 ditentukan sendiri dan B2 sesuai dengan B hasil perhitungan. Gamma diharapkan memiliki nilai yang sekecil mungkin. Berdasarkan uji coba yang dilakukan, nilai gamma yang digunakan adalah sebesar 1.605. Nilai A dan B yang digunakan juga merupakan hasil augmentasi agar didapatkan 5 parameter K sehingga dapat dikalikan dengan pinv(b) seperti pada kasus LQR. Sehingga didapatkan parameter PID Kd =0.4125
Kp=1.0967
Ki =1.9054
Hasil PID juga diplotkan dengan simulasi Simulink pada Gambar 5.1 sehingga didapatkan hasil sebagai berikut
Gambar 6.1 Simulasi Simulink respon sistem terhadap waktu dengan metoda PID ARE Berdasarkan hasil di atas, besar overshoot dengan metode ARE mencapai 20% dengan settling time sebesar 3 s yang lebih lama dari metoda sebelumnya. Oleh karena itu, sistem PID dengan LQR menghasilkan sistem yang cenderung lebih stabil dari pada metode H-∞. PERBANDINGAN PID PADA LQR DAN ARE
Gambar 6.2 Perbandingan respon sistem dengan metoda PID LQR dan PID ARE
Warna biru memperlihatkan sistem yang dikontrol dengan PID LQR, sedangkan yang berwarna ungu menunjukkan sistem yang dikontrol dengan PID ARE. Meskipun dengan overshoot yang lebih besar 4%, pengontrol PID LQR memberikan respon sistem yang lebih cepat mencapai keadaan tunak. Karena itu, untuk sistem TORA ini, metoda PID LQR memberikan hasil yang lebih baik dibanding metoda PID ARE.
VII.
GMDS – Plant mass damper spring GMDS menyatakan hubungan output-input pada sistem mass-damper-spring. Representasi
ruang keadaan ditunjukkan oleh:
di mana:
Sistem TORA dalam hal ini didekati sebagai sistem mass-damper-spring. Pada sistem ini, telah diinisiasi nilai parameter-parameter sistem, yaitu: No
Besaran
Parameter
Nilai
Satuan
1
Cart Mass
M
1.3
Kg
2
Arm Mass
m
0.1
Kg
3
Arm Eccentricity
e
0.06
M
4
Arm Inertia
I
0.0002
Kg m2
5
Spring Stiffness
K
186.3
N/m
6
epsilon
Ԑ
0.2
-
7
Torque
N
1.3
mN
8
Simpangan
q
0.1
M
Pada keadaan riil, parameter-parameter tersebut tidak dapat diketahui berapa nilai persisnya. Namun, nilai yang tercantum di atas dapat diperkirakan dengan disertai nilai ketidakpastian ( ). Ketidakpastian ini dinyatakan dalam rentang, sehingga menghasilkan nilai parameter baru dengan persamaan berikut:
Dengan mengambil nilai parameter baru (dengan
= 0.4, dan
= 0.3, dan -1 <
,
< -1, diperoleh nilai
):
̅ ̅ Dengan memasukkan nilai-nilai parameter di atas pada program di bawah dengan perangkat lunak MATLAB: %gmds m = 2.167; c = 0; k = 266.142; pm = 0.4; pc = 0.2; pk = 0.3; A = [ 0 1 -k/m -c/m]; B1 = [ 0 0 0 -pm -pc/m -pk/m]; B2 = [ 0 1/m]; C1 = [-k/m -c/m 0 c k 0]; C2 = [ 1 0 ]; D11 = [-pm -pc/m -pk/m 0 0 0 0 0 0]; D12 = [1/m 0 0 ]; D21 = [0 0 0];
D22 = 0; G = pck(A,[B1,B2],[C1;C2],[D11 D12;D21 D22]);
diperoleh matriks GMDS sebagai berikut:
Kemudian, GMDS yang diperoleh tersebut dianalisis respon frekuensinya. Dengan program MATLAB berikut: % Frequency responses of the perturbed plants % gmds omega = logspace(-1,1,100); [delta1,delta2,delta3] = ndgrid([-1 0 1],[-1 0 1], ... [-1 0 1]); for j = 1:27 delta = diag([delta1(j),delta2(j),delta3(j)]); olp = starp(delta,G); olp_ic = sel(olp,1,1); olp_g = frsp(olp_ic,omega); figure(1) vplot('bode',olp_g,'c-') subplot(2,1,1) hold on subplot(2,1,2) hold on end subplot(2,1,1) olp_ic = sel(G,4,4); olp_g = frsp(olp_ic,omega); vplot('bode',olp_g,'r--') subplot(2,1,1) title('BODE PLOTS OF PERTURBED PLANTS') hold off subplot(2,1,2) hold off
diperoleh bode plot dari sistem TORA yang melibatkan ketidakpastian sebagai berikut:
Gambar 7.1 Respon frekuensi dari plant TORA yang didekati sebagai mass damper spring dengan ketidakpastian Grafik magnituda terhadap frekuensi di atas menunjukkan bahwa pada frekuensi 1 rad/s dan fasa 180°, magnituda tidak melebihi 1 dB. Hal ini sesuai dengan Small Gain Theorem yang menunjukkan bahwa sistem stabil. Namun, respon frekuensi ini juga menunjukkan hasil yang belum sesuai dengan yang diinginkan. Profil closed-loop gain yang ditunjukkan diharapkan sesuai seperti pada Gambar 4.3, di mana noise pada frekuensi tinggi diminmalisir, dan sinyal asli pada frekuensi rendah dipertahankan. Maka, untuk memeroleh hasil yang lebih optimal, diberikan variasi pada nilai Dari variasi nilai yang dilakukan, dipilih sebagai berikut. Matriks GMDS:
= 110 dan
dan
.
= 100 sehingga diperoleh hasil
Dan respon frekuensi sebagai berikut.
Gambar 7.2 Respon frekuensi GMDS dengan
=110 dan
=100
Gambar 7.2 menunjukkan respon sistem yang lebih optimal. Terlihat pada frekuensi 1 rad/s, sistem mulai berubah fasanya dan menunjukkan karakter yang mirip dengan low-pass filter, meneruskan sinyal asli pada frekuensi rendah dan meredam noise pada frekuensi tinggi. Respon ini juga memenuhi Small Gain Theorem, yang menunjukkan bahwa sistem ini stabil.
VIII. Kesimpulan 1. Sistem TORA ini Controllable dan Observable 2. Nilai Q dan R berbanding terbalik terhadap hasil respon sistem. Ketika Q diperkecil/ R diperbesar, sistem menjadi tidak stabil dan nilainya jauh dari input 3. Berdasarkan analisis frekuensi yang ditunjukkan oleh grafik sensitivity dan tracking, sistem menunjukkan respon yang baik (meminimalisir noise dan disturbance) 4. Sistem PID dengan LQR menghasilkan sistem yang cenderung lebih stabil dari pada metode H-∞.
IX.
Daftar Pustaka
D.-W. Gu, P. Hr. Petkov and M. M. Konstantinov. Robust Control Design with MATLAB. 2005. Springer: Glasgow, Scotland, U.K. Ogata, Katsuhiko. Modern Control Engineering. 2010. Prentice Hall: United States of America R. J. Kinsey, D. L. Mingori, and R. H. Rand, “Nonlinear Controller to Reduce Resonance Effects of A Dual-Spin Spacecraft Through Precession Phase Lock,” Proc. Of the 31st Conference on Decision and Control, pp.3025-3030, 1992. Tavakoli Mahdi, Hamid D. Taghirad and Mehdu Abrishamchian, “Identification and Robust H∞ Control of The Rotational/ Translational Actuator System.” International Journal of Control, Automation and System, vol 3 no 3 pp.387-396, 2005 Zhou, Kemin and Doyle, Jhon C. Essentials of Robust Control. 1999. Prentice Hall: United States of America
X.
Lampiran
Program MATLAB: %kondisi awal yang diketahui q= 0.1; M=1.3; m=0.1; e=0.06; I=0.0002; k=186.3; epsilon= 0.2; N=1.3; pi=3.14; u=((M+m)/(k*(I+m*e^2)))*N; x(1) = square((M+m)/(I+m*e^2))*q; x(2)=0; x(3)=-20/(180*pi); x(4)=-epsilon*x(2)*cos(x(3))+u; x=[ x(1); x(2); x(3) ; x(4) ] %matriks A= matriks f(x(t)) A= [0 1 0 0; -1/(1-epsilon^2*cos(x(3))^2) 0 0 (epsilon*x(4)*sin(x(3)))/(1epsilon^2*cos(x(3))^2); 0 0 0 1 (epsilon*cos(x(3)))/(1-epsilon^2*cos(x(3))^2) 0 0 (x(4)*epsilon*sin(x(3)))/(1-epsilon^2*cos(x(3))^2)]; %matriks B= matriks g(x(t)) B= [ 0 ; -epsilon*cos(x(3))/(1-epsilon^2*cos(x(3))^2) ; 0 ; 1/(epsilon^2*cos(x(3)^2))]; C= [0 0 1 0]; D = [0]; % matriks A B dan C diaugmentasi untuk penghitungan parameter PID pada LQR dan ARE Aaug=[ A B ; 0 0 0 0 0] Baug= [ 0 ; 0 ; 0 ; 0 ; 1] Caug=[ C 0 ] %------------------------------------------------------------------------------------------------------------------%Transcfer function system [Num,Den]=ss2tf(A,B,C,D) G=tf(Num,Den) %------------------------------------------------------------------------------------------------------------------%Grafik respon sistem awal SS=ss(A,B,C,D); figure(1) step(SS) %respon sistem jika diberi input step title('Respon Waktu Sistem Awal') figure(2) bode(SS) title('Respon Frekuensi Sistem Awal') figure(3) rlocus(SS)
title ('Respon root locus') %------------------------------------------------------------------------------------------------------------------%cek CTRB dari sistem TORA CB=ctrb(A,B); rank (CB); [Mcb,Ncb]=size(CB); OB=obsv(A,C); rank (OB); [Mob,Nob]=size(OB); if rank(CB)==Mcb disp('SISTEM CONTROLABLE'); else disp('SISTEM NON CONTROLABLE'); end if rank (OB)==Mob disp ('SISTEM OBSERVABLE'); else disp ('SISTEM TIDAK OBSERVABLE') end %------------------------------------------------------------------------------------------------------------------%METODE LQR %harga R dan Q ddiubah2 agar didapat nilai respon sistem yg paling baik Q3=eye(size(A)); % memakai R R1=0.0001; R2=0.01; R3=1; R4=100; R5=10000; ts=0:1:50; figure(4) h1=subplot(5,1,1) [K,P,E]=lqr(A,B,Q3,R1); A1=A-B*K; B1=B*K(1); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=1 dan R=0.0001'); hold all h2=subplot(5,1,2); [K,P,E]=lqr(A,B,Q3,R2); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=1 dan R=0.01'); hold all h3=subplot(5,1,3);
[K,P,E]=lqr(A,B,Q3,R3); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=1 dan R=1'); hold all h4=subplot(5,1,4); [K,P,E]=lqr(A,B,Q3,R4); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=1 dan R=100'); hold all h4=subplot(5,1,5); [K,P,E]=lqr(A,B,Q3,R5); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=1 dan R=10000'); hold all Q1=0.0001*eye(size(A)); Q2=0.01*eye(size(A)); Q4=100*eye(size(A)); Q5=10000*eye(size(A)); figure(5) h1=subplot(5,1,1) [K,P,E]=lqr(A,B,Q1,R3); A1=A-B*K; B1=B*K(1); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=0.0001 dan R=1'); hold all h2=subplot(5,1,2); [K,P,E]=lqr(A,B,Q2,R3); A1=A-B*K; B1=B*K(4); C1=C; D1=D;
sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=0.01dan R=1'); hold all h3=subplot(5,1,3); [K,P,E]=lqr(A,B,Q3,R3); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=1 dan R=1'); hold all h4=subplot(5,1,4); [K,P,E]=lqr(A,B,Q4,R3); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=100 dan R=1'); hold all h4=subplot(5,1,5); [K,P,E]=lqr(A,B,Q5,R3); A1=A-B*K; B1=B*K(4); C1=C; D1=D; sys_lqr=ss(A1,B1,C1,D1); [y,ts]=step(sys_lqr,ts); plot(ts,y) title('Q=10000 dan R=1'); hold all %------------------------------------------------------------------------------------------------------------------%berdasarkan hasil manipulasi harga Q dan R didapatkan bahwa sistem lebih %stabil ketika harga Q sangat kecil atau ketika harga R sangat besar %seingga digunakan kombinasi harga keduanya yaitu Q=0.0001 dan R=10000 %penentuan K,P,E [K,P,E]=lqr(A,B,Q1,R5) AA=A-B*K; BB=B*K(4); CC=C; DD=D; t=0:1:50; figure(6) [x,y,t]=step(AA,BB,CC,DD,1,t); %plot respon sistem plot(t,y);grid legend('x1','x2','x3','x4','input') title('plot pengontrolan LQR') figure(7) SS1=ss(AA,BB,CC,DD);
rlocus(SS1) %cek posisi pole sistem untuk menunjukkan kestabilan sistem figure(8) bode(SS1) title('respon frekuensi sistem berpengontrol')% cek sistem pada domain frekuensi %plot T S L [num, den]=ss2tf(AA, BB, CC, DD); TF2=tf(num,den); Sensitivity=feedback(1,TF2); CSensitivity=feedback(TF2,1); %plot figure(9) h1=plot(1); bodemag(TF2); %hanya magnituda yang diplotkan title('open loop bode plot'); %menunjukan sistem kalang terbuka hold all figure(10) h2=plot(2); bodemag(Sensitivity); title('Sensitivity'); %meneunjukan sensitivity figure(11) hold all h3=plot(3); bodemag(CSensitivity); title('Tracking'); %menunjukan tracking sistem hold all %----------------------------------------------------------------------------------------------------------%PENENTUAN PID %PID dengan LQR Q=5*eye(5) R=1; [Kpid,Ppid,Epid]=lqr (Aaug,Baug,Q,R) b = [C 0; C*A C*B; C*A*A C*A*B]; K_hat = Kpid*pinv(b); % Parameter PID Ki_hat = K_hat(1,3); Kd = Ki_hat/(1+Ki_hat*C*B) % Parameter Derivatif K2_hat = K_hat(1,2); Kp = K2_hat*(1-Ki_hat*C*B) % Parameter Proporsional K1_hat = K_hat(1,1); Ki = K1_hat*(1-Ki_hat*C*B) % Parameter Integral %pembuatan plot PID LQR K = [Kpid(1) Kpid(2) Kpid(3) Kpid(4)]; Apid = Aaug-Baug*Kpid; Bpid = Baug*Kpid(5); Cpid = Caug; Dpid = 0; t=0:1:50; figure(12) [x2,y2,t]= step(Apid,Bpid,Cpid,Dpid,1,t); SS2=ss(Apid,Bpid,Cpid,Dpid); plot (t,y2) grid
legend('x1','x2','x3','x4','input') title('plot pengontrolan PID LQR') %perbandingan LQR dan PID LQR figure(13) step(SS1,SS2) legend('LQR','PID LQR') %------------------------------------------------------------------------------------------------------------------% H tak hingga PID B1=[0 0 0 0 1]; B2=Baug; gama=1.05; %dimanipulasi agar sistem dapat dikontrol, ketika menggunakan gama yang kecil sitem tidak terkontrol Aare=Aaug; Bare=(B2*B2')-(gama^-2*B1*B1'); Care=Caug'*Caug; Dare=[0]; Xare=ARE(Aare,Bare,Care); Kare=-Baug'*Xare; Kare_hat = Kare*pinv(b); % Parameter PID Kiare_hat = Kare_hat(1,3); Kdare = Kiare_hat/(1+Kiare_hat*C*B) % Parameter Derivatif K2are_hat = Kare_hat(1,2); Kpare = K2are_hat*(1-Kiare_hat*C*B) % Parameter Proporsional K1are_hat = Kare_hat(1,1); Kiare = K1are_hat*(1-Kiare_hat*C*B) % Parameter Integral %------------------------------------------------------------------------------------------------------------------%gmds m = 2.167; c = 0; k = 266.142; pm = 110; pc = 0.2; pk = 100; A = [ 0 1 -k/m -c/m]; B1 = [ 0 0 0 -pm -pc/m -pk/m]; B2 = [ 0 1/m]; C1 = [-k/m -c/m 0 c k 0]; C2 = [ 1 0 ]; D11 = [-pm -pc/m -pk/m 0 0 0 0 0 0]; D12 = [1/m
0 0 ]; D21 = [0 0 0]; D22 = 0; G = pck(A,[B1,B2],[C1;C2],[D11 D12;D21 D22]); %------------------------------------------------------------------------------------------------------------------% Frequency responses of the perturbed plants % gmds omega = logspace(-1,1,100); [delta1,delta2,delta3] = ndgrid([-1 0 1],[-1 0 1], ... [-1 0 1]); for j = 1:27 delta = diag([delta1(j),delta2(j),delta3(j)]); olp = starp(delta,G); olp_ic = sel(olp,1,1); olp_g = frsp(olp_ic,omega); figure(1) vplot('bode',olp_g,'c-') subplot(2,1,1) hold on subplot(2,1,2) hold on end subplot(2,1,1) olp_ic = sel(G,4,4); olp_g = frsp(olp_ic,omega); vplot('bode',olp_g,'r--') subplot(2,1,1) title('BODE PLOTS OF PERTURBED PLANTS') hold off subplot(2,1,2) hold off
