Bài tập dài môn điều khiển tự động

1 BÀI TẬP DÀI MÔN ĐIỀU KHIỂN TỰ ĐỘNG PHẦN I/ YÊU CẦU THIẾT KẾ 1) Yêu cầu thiết kế hệ thống điềukhiển tự động có: - Bộ điều khiển PID có hàm truyền : WPID(s) = KP*(1+ sTI * 1 +TD*s) - Đối tượng điều khiển là khâu trễ và quán tính bậc nhất có hàm truyền: WDT(s) = 1* + − sT e Ls Với các tham số: L/T = 0.9; T = 95; Hệ thống có sơ đồ như hình vẽ: Wpid(s) Wdt(s) 2) Tính toán các tham số Kp, Ti, Td đảm bảo tính ổn định và chất lượng của hệ thống (theo Ziegler – Nichols). 3) Xét ổn định của hệ thống. Tìm các điểm cực và điểm không. 4) Khảo sát chất lượng, chọn các tham số với các luật điều khiển P, PI, PID đảm bảo cho hệ thống có chất lượng tốt nhất (chỉnh định bằng tay). 5) Tính tham số tối ưu của bộ điều khiển PID dùng hàm least – squares (sai số bình phương nhỏ nhất) với các tham số L và T đã cho ở trên. PHẦN II/ TÍNH TOÁN VỚI MATLAB I/ Tính toán các tham số Kp, Ti, Td theo tiêu chuẩn Ziegler – Nichols Áp dụng bảng công thức thực nghiệm của Ziegler – Nichols: Luật điều khiển Kp Ti Td P T/L ∞ 0 PI 0.9T/L 10L/3 0 PID 1.2T/L 2L 0.5L Với các tham số đã cho: L/T = 0.9; T = 95 ⇒L=85.5 Thay vào bảng trên ta có: Luật điều khiển Kp Ti Td P 1.111 ∞ 0 PI 1 285 0 PID 1.333 171 42.75 2 II/ Phân tích hệ thống với các tham số chọn theo phương pháp Ziegler – Nichols và hiệu chỉnh bằng tay Để có thể phân tích hệ thống, ta cần khai triển Taylor khâu trễ e-Ls, lấy 3 số hạng đầu là đủ chính xác. Việc khai triển này đựơc thực hiện nhờ MatLab qua hàm xấp xỉ Pade: Chương trình MatLab như sau: >> [ts,ms]=pade(85.5,3) ts = -1.0000 0.1404 -0.0082 0.0002 ms = 1.0000 0.1404 0.0082 0.0002 >> wtre=tf(ts,ms) Transfer function: -s^3 + 0.1404 s^2 - 0.008208 s + 0.000192 --------------------------------------------------- s^3 + 0.1404 s^2 + 0.008208 s + 0.000192 1/ Khảo sát hệ thống với luật điều khiển tỉ lệ P: >> T=95;L=85.5; >> [ts,ms]=pade(L,3) ts = -1.0000 0.1404 -0.0082 0.0002 ms = 1.0000 0.1404 0.0082 0.0002 >> wdt=tf(ts,ms)*tf(1,[T 1]) Transfer function: -s^3 + 0.1404 s^2 - 0.008208 s + 0.000192 ------------------------------------------------------ 95 s^4 + 14.33 s^3 + 0.9201 s^2 + 0.02645 s + 0.000192 >> Kp=T/L Kp = 1.1111 >> wkinp=feedback(Kp*wdt,1) Transfer function: -1.111 s^3 + 0.1559 s^2 - 0.00912 s + 0.0002133 ------------------------------------------------------ 95 s^4 + 13.22 s^3 + 1.076 s^2 + 0.01733 s + 0.0004053 >> step(wkinp) >> nyquist(wkinp) >> pzmap(wkinp) >> [p,z]=pzmap(wkinp) p = -0.0628 + 0.0722i -0.0628 - 0.0722i -0.0068 + 0.0205i 3 -0.0068 - 0.0205i z = 0.0430 + 0.0410i 0.0430 - 0.0410i 0.0543 Đặc tính quá độ Đặc tính tần số 4 Đồ thị điểm không, điểm cực *Nhận xét: - Các điểm cực đều nằm bên trái trục ảo nên hệ thống ổn định. - Độ quá điển khiển maxδ = 39% - Thời gian tăng tốc (Rise time): 54.1s - Thời gian quá độ (Settling time): 573s. *Hiệu chỉnh để có đặc tính quá độ tốt hơn: cho Kp = 0.3 >> Kp = 0.3; >> wkinp = feedback(Kp*wdt,1) Transfer function: -0.3 s^3 + 0.04211 s^2 - 0.002462 s + 5.76e-005 ------------------------------------------------------- 95 s^4 + 14.03 s^3 + 0.9622 s^2 + 0.02398 s + 0.0002496 >> step(wkinp) - Sau khi hiệu chỉnh có: + Độ quá điều chỉnh maxδ = 1.24% + Thời gian quá độ: 239 (s) 5 + Thời gian tăng tốc: 112 (s) 2/ Khảo sát hệ thống với luật điều chỉnh tỉ lệ – tích phân PI Chương trình Matlab: >> T=95; L=85.5 >> Kp=0.9*T/L;Ti=10*L/3; >> [ts,ms]=pade(L,3) ts = -1.0000 0.1404 -0.0082 0.0002 ms = 1.0000 0.1404 0.0082 0.0002 >> wdt=tf(ts,ms)*tf(1,[T 1]) Transfer function: -s^3 + 0.1404 s^2 - 0.008208 s + 0.000192 ------------------------------------------------------ 95 s^4 + 14.33 s^3 + 0.9201 s^2 + 0.02645 s + 0.000192 >> wpi=tf([Kp*Ti Kp],[Ti 0]) Transfer function: 285 s + 1 --------- 285 s >> wkinpi=feedback(wdt*wpi,1) Transfer function: -285 s^4 + 39 s^3 - 2.199 s^2 + 0.04651 s + 0.000192 ------------------------------------------------------------------ 27075 s^5 + 3800 s^4 + 301.2 s^3 + 5.339 s^2 + 0.1012 s + 0.000192 >> step(wkinpi) >> nyquist(wdt*wpi) >> pzmap(wkinpi) >> [p,z]=pzmap(wkinpi) p = -0.0618 + 0.0694i -0.0618 - 0.0694i -0.0073 + 0.0184i -0.0073 - 0.0184i -0.0021 z = 0.0430 + 0.0410i 0.0430 - 0.0410i 0.0543 -0.0035 6 Đặc tính quá độ Đặc tính tần số Đồ thị điểm không - điểm cực 7 * Nhận xét: - Các điểm cực nằm bên trái trục ảo nên hệ thống ổn định. - Độ quá điều chỉnh: maxδ = 0% - Thời gian tăng tốc (Rise time): 675s - Thời gian quá độ (Settling time): 1500s * Hiệu chỉnh để đạt chất lượng tốt hơn: cho Kp = 0.65; Ti = 110; >> Kp=0.65; Ti=110; >> wpi=tf([Kp*Ti Kp],[Ti 0]) Transfer function: 71.5 s + 0.65 ------------- 110 s >> wkinpi=feedback(wdt*wpi,1) Transfer function: -71.5 s^4 + 9.385 s^3 - 0.4956 s^2 + 0.008392 s + 0.0001248 ------------------------------------------------------------------------------------- 10450 s^5 + 1505 s^4 + 110.6 s^3 + 2.414 s^2 + 0.02951 s + 0.0001248 >> step(wkinpi) Sau khi hiệu chỉnh ta có: - Độ quá điều chỉnh: maxδ = 1.94% - Thời gian quá độ: 286 (s) - Thời gian tăng tốc: 148 (s) 3/ Khảo sát hệ thống với luật điều chỉnh PID: Chương trình Matlab như sau: >> L=85.5;T=95;Kp=1.2*T/L;Ti=2*L;Td=0.5*L; >> [ts,ms]=pade(85.5,3) ts = -1.0000 0.1404 -0.0082 0.0002 ms = 8 1.0000 0.1404 0.0082 0.0002 >> wdt=tf(ts,ms)*tf(1,[T 1]) Transfer function: -s^3 + 0.1404 s^2 - 0.008208 s + 0.000192 ------------------------------------------------------ 95 s^4 + 14.33 s^3 + 0.9201 s^2 + 0.02645 s + 0.000192 >> wpid=tf([Kp*Ti*Td Kp*Ti Kp],[Ti 0]) Transfer function: 9747 s^2 + 228 s + 1.333 ------------------------ 171 s >> wkin=feedback(wpid*wdt,1) Transfer function: -9747 s^5 + 1140 s^4 - 49.33 s^3 + 0.1871 s^2 + 0.03283 s + 0.000256 ------------------------------------------------------------------------------------- 6498 s^5 + 3591 s^4 + 108 s^3 + 4.71 s^2 + 0.06566 s + 0.000256 >> step(wkin) >> nyquist(wkin) >> pzmap(wkin) >> [p,z]=pzmap(wkin) p = -0.5235 -0.0060 + 0.0322i -0.0060 - 0.0322i -0.0102 -0.0069 z = 0.0430 + 0.0410i 0.0430 - 0.0410i 0.0543 -0.0117 -0.0117 9 Đặc tính quá độ Đặc tính tần số Đồ thị điểm không - điểm cực 10 * Nhận xét: - Các điểm cực nằm bên trái trục ảo nên hệ thống ổn định, tuy nhiên do có 2 điểm cực phức liên hợp nằm sát trục ảo nên độ dự trữ ổn định nhỏ. - Độ quá điều chỉnh: maxδ =-250% - Thời gian tăng tốc: 108 (s) - Thời gian quá độ: 477 (s) *Hiệu chỉnh để có đặc tính tốt hơn: Kp=0.7; Ti=110; Td=10 Chương trình MatLab: >> Kp=0.7; Ti=110;Td=10; >> wpid=tf([Kp*Ti*Td Kp*Ti Kp],[Ti 0]) Transfer function: 770 s^2 + 77 s + 0.7 -------------------- 110 s >> wkin=feedback(wpid*wdt,1) Transfer function: -770 s^5 + 31.07 s^4 + 3.787 s^3 - 0.3859 s^2 + 0.009038 s + 0.0001344 ---------------------------------------------------------------------- 9680 s^5 + 1608 s^4 + 105 s^3 + 2.523 s^2 + 0.03016 s + 0.0001344 >> step(wkin) Đặc tính quá độ sau khi hiệu chỉnh - Độ quá điều chỉnh: maxδ =1.57% - Thời gian tăng tốc: 224 (s) - Thời gian quá độ: 284 (s) III/ Mô hình hoá hệ thống bằng Simulink: Bước 1/ Vào Simulink mô hình hoá hệ thống dưới dạng sơ đồ khối như sau: 11 Bước 2/ Nạp tham số cho các khâu: - Bộ điều khiển PID: Các tham số Kp, Ki, Kd - Khâu bão hoà Saturation: + Upper limit: 1.5 + Lower limit: -1.5 - Khâu hạn chế gia tốc Rate Limiter: + Rising slew rate: 1 + Faling slew rate: -1 12 - Khâu trễ Transport Delay: + Time Delay: 85.5 (ứng với L=85.5). - Khâu hàm truyền Transfer Fcn: + Numerator: [1] + Denominator: [95 1] - Ghi chương trình mô phỏng lại thành file “BaiTapDai” trong thư mục work của MatLab. Bước 2/ Trong cửa sổ Command Windows của MatLab vào File → New → M-file → sẽ hiện ra cửa sổ soạn thảo. Trong cửa sổ này ta viết các hàm tính toán tham số tối ưu cho bộ điều khiển: * Hàm thứ nhất: function e=yout_1(pid,model,tspan) opt=simset('solver','ode5'); assignpid(pid); [t,x,y]=sim(model,tspan,opt); e=y-1; 13 Ghi lại hàm này thành file: yout_1.m * Hàm thứ hai: function assignpid(pid) assignin('base','Kp',pid(1)); assignin('base','Ki',pid(2)); assignin('base','Kd',pid(3)); Ghi lại hàm này thành file: assignpid.m * Hàm thứ ba: function pid = lsqpid(pid0,model,tspan) options=foptions; options=[1 0.01 0.01]; pid=leastsq('yout_1',pid0,options,[],model,tspan); Ghi lại hàm này thành file: lsqpid.m Bước 3/ Trở về cửa sổ Command Windows viết chương trình: >> Kp=0.7; Ti=110;Td=10;Ki=Kp/Ti;Kd=Kp*Td; >> pid0=[Kp Ki Kd]; >> pid=lsqpid(pid0,'BaiTapDai',[0 5000]) Warning: The function LEASTSQ is obsolete and has been replaced by LSQNONLIN. LEASTSQ will be removed in a future release of the Optimization Toolbox. Update your code to call LSQNONLIN to suppress this warning message. See "Converting your code to Version 2 Syntax" in the Optimization Toolbox User's Guide for more information. > In C:\ChuongTrinhUngDung\Matlab\toolbox\optim\leastsq.m at line 53 In C:\ChuongTrinhUngDung\Matlab\work\lsqpid.m at line 4 f-COUNT RESID STEP-SIZE GRAD/SD LAMBDA 4 2.02681 1 -0.0534 10 2.00227 1 0.00408 1.03732 17 2.00085 1.32 2.44e-007 0.518662 24 2.0005 1.37 -9.51e-009 0.219114 30 2.00048 1.12 -1.22e-010 0.103279 Optimization Terminated Successfully pid = 0.7757 0.0071 12.8041 * Kết quả chương trình cho ta tham số tính toán tối ưu thu được: Kp = 0.7757; Ki = 0.0071; Kd = 12.8041 Bước 4/ Quay lại Simulink nháy vào Scope ta thu được đặc tính quá độ của hệ thống với tham số tính toán tối ưu tìm được: 14 - Độ quá điều khiển: maxδ =1.65% - Thời gian quá độ: 241(s) KẾT LUẬN Các tham số tối ưu của bộ điều khiển PID ứng với đối tượng đã cho là: Kp = 0.7757 Ki = 0.0071 Kd = 12.8041 Ứng với các tham số này, đặc tính quá độ thu được của hệ thống là tốt nhất.