›
›

Giáo trình Matlab trong điều khiển tự động

pdf 58 trang phuongnguyen 40
Download
Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình Matlab trong điều khiển tự động", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

  • pdfgiao_trinh_matlab_trong_dieu_khien_tu_dong.pdf

Nội dung text: Giáo trình Matlab trong điều khiển tự động

  1.   GIÁO TRÌNH MATLAB TRONG ĐIỀU KHIỂN TỰ ĐỘNG NHĨM LỆNH VỀ ĐÁP ỨNG TẦN SỐ
  2. MATLAB trong điều khiển tự động Trang 159 NHÓM LỆNH VỀ ĐÁP ỨNG TẦN SỐ (Frequency Response) 1. Lệnh BODE a) Công dụng: Tìm và vẽ đáp ứng tần số giản đồ Bode. b) Cú pháp: [mag,phase,w] = bode(a,b,c,d) [mag,phase,w] = bode(a,b,c,d,iu) [mag,phase,w] = bode(a,b,c,d,iu,w) [mag,phase,w] = bode(num,den) [mag,phase,w] = bode(num,den,w) c) Giải thích: Lệnh bode tìm đáp ứng tần số biên độ và pha của hệ liên tục LTI. Giản đồ Bode dùng để phân tích đặc điểm của hệ thống bao gồm: biên dự trữ, pha dự trữ, độ lợi DC, băng thông, khả năng miễn nhiễu và tính ổn định. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh bode sẽ vẽ ra giản đồ Bode trên màn hình. bode(a,b,c,d) vẽ ra chuỗi giản đồ Bode, mỗi giản đồ tương ứng với một ngõ vào của hệ không gian trạng thái liên tục: . x= Ax + Bu y = Cx + Du với trục tần số được xác định tự động. Nếu đáp ứng thay đổi nhanh thì cần phải xác định nhiều điểm hơn. bode(a,b,c,d,iu) vẽ ra giản đồ Bode từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng giản đồ Bode. bode(num,den) vẽ ra giản đồ Bode của hàm truyền đa thức hệ liên tục G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. bode(a,b,c,d,iu,w) hay bode(num,den,w) vẽ ra giản đồ Bode với vector tần số w do người sử dụng xác định. Vector w chỉ ra các điểm tần số (tính bằng rad/s) mà tại đó đáp ứng tần số giản đồ Bode được tính. Nếu vẫn giữ lại các đối số ở vế trái của dòng lệnh thì: [mag,phase,w] = bode(a,b,c,d) [mag,phase,w] = bode(a,b,c,d,iu) [mag,phase,w] = bode(a,b,c,d,iu,w) [mag,phase,w] = bode(num,den) [mag,phase,w] = bode(num,den,w)
  3. MATLAB trong điều khiển tự động Trang 160 Sẽ không vẽ ra giản đồ Bode mà tạo ra các ma trận đáp ứng tần số mag, phase và w của hệ thống. Ma trận mag và phase có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector w. G(s) = C(sI –A)-1B + D mag(ω) = ⏐G(jω)⏐ phase(ω) = ∠G(jω) Góc pha được tính bằng độ. Giá trị biên độ có thể chuyển thành decibel theo biểu thức: magdB = 20*log10(mag) Chúng ta có thể dùng lệnh fbode thay cho lệnh bode đối với các hệ thống có thể chéo nhau. Nó sử dụng các thuật giải nhanh hơn dựa trên sự chéo hóa của ma trận hệ thống A. d) Ví dụ: Vẽ đáp ứng biên độ và pha của hệ bậc 2 với tần số tự nhiên ωn= 1 và hệ số tắt dần ζ = 0.2 [a,b,c,d] = ord2(1,0.2); bode(a,b,c,d) grid on và ta được giản đồ Bode đáp ứng tần số của hệ thống như sau: Bode Diagrams 0 -10 -20 -30 -40 0 -50 Phase (deg); Magnitude (dB) Magnitude (deg); Phase -100 -150 -1 0 1 10 10 10 Frequency (rad/sec) 2. Lệnh FBODE a) Công dụng: Vẽ đáp ứng tần số giản đồ Bode cho hệ tuyến tính liên tục. b) Cú pháp: [mag,phase,w] = fbode(a,b,c,d) [mag,phase,w] = fbode(a,b,c,d,iu)
  4. MATLAB trong điều khiển tự động Trang 161 [mag,phase,w] = fbode(a,b,c,d,iu,w) [mag,phase,w] = fbode(num,den) [mag,phase,w] = fbode(num,den,w) c) Giải thích: Lệnh fbode tìm nhanh đáp ứng tần số biên độ và pha của hệ liên tục LTI. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh fbode sẽ vẽ ra giản đồ Bode trên màn hình. fbode(a,b,c,d) vẽ ra chuỗi giản đồ Bode, mỗi giản đồ tương ứng với một ngõ vào của hệ không gian trạng thái liên tục: . x= Ax + Bu y = Cx + Du với trục tần số được xác định tự động. Nếu đáp ứng thay đổi nhanh thì cần phải xác định nhiều điểm hơn. fbode(a,b,c,d,iu) vẽ ra giản đồ Bode từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng giản đồ Bode. fbode nhanh hơn nhưng kém chính xác hơn bode. fbode(num,den) vẽ ra giản đồ Bode của hàm truyền đa thức hệ liên tục G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. fbode(a,b,c,d,iu,w) hay fbode(num,den,w) vẽ ra giản đồ Bode với vector tần số w do người sử dụng xác định. Vector w chỉ ra các điểm tần số (tính bằng rad/s) mà tại đó đáp ứng tần số giản đồ Bode được tính. Nếu vẫn giữ lại các đối số ở vế trái của dòng lệnh thì: [mag,phase,w] = fbode(a,b,c,d) [mag,phase,w] = fbode(a,b,c,d,iu) [mag,phase,w] = fbode(a,b,c,d,iu,w) [mag,phase,w] = fbode(num,den) [mag,phase,w] = fbode(num,den,w) sẽ không vẽ ra giản đồ Bode mà tạo ra các ma trận đáp ứng tần số mag, phase và w của hệ thống. Ma trận mag và phase có số cột bằng số ngõ ra và có số hàng là length(w). d) Ví dụ: Vẽ đáp ứng biên độ và pha của hệ bậc 2 với tần số tự nhiên ωn= 1 và hệ số tắt dần ζ = 0.2 [a,b,c,d] = ord2(1,0.2); fbode(a,b,c,d); grid on và ta được đáp ứng như sau:
  5. MATLAB trong điều khiển tự động Trang 162 Bode Diagrams 0 -10 -20 -30 -40 0 -50 Phase (deg); Magnitude (dB) Magnitude (deg); Phase -100 -150 -1 0 1 10 10 10 Frequency (rad/sec) 3. Lệnh DBODE a) Công dụng: Tìm và vẽ đáp ứng tần số giản đồ Bode của hệ gián đoạn. b) Cú pháp: [mag,phase,w] = dbode(a,b,c,d,Ts) [mag,phase,w] = bode(a,b,c,d,Ts,iu) [mag,phase,w] = bode(a,b,c,d,Ts,iu,w) [mag,phase,w] = bode(num,den,Ts) [mag,phase,w] = bode(num,den,Ts,w) c) Giải thích: Lệnh dbode tìm đáp ứng tần số biên độ và pha của hệ liên tục LTI. Lệnh dbode khác với lệnh freqz mà trong đó đáp ứng tần số đạt được với tần số chưa chuẩn hóa. Đáp ứng có được từ dbode có thể được so sánh trực tiếp với đáp ứng lệnh bode của hệ thống liên tục tương ứng. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh dbode sẽ vẽ ra giản đồ Bode trên màn hình. dbode(a,b,c,d,Ts) vẽ ra chuỗi giản đồ Bode, mỗi giản đồ tương ứng với một ngõ vào của hệ không gian trạng thái liên tục: x[n+] = Ax[n] + Bu{n] y[n] = Cx[n] + Du[n] với trục tần số được xác định tự động. Các điểm tần số được chọn trong khoảng từ π/Ts (rad/sec), trong đó π/Ts (rad/sec) tương ứng với nửa tần số lấy mẫu (tần số Nyquist). Nếu đáp ứng thay đổi nhanh thì cần phải xác định nhiều điểm hơn. Ts là thời gian lấy mẫu.
  6. MATLAB trong điều khiển tự động Trang 163 dbode(a,b,c,d,Ts,iu) vẽ ra giản đồ Bode từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng giản đồ Bode. dbode(num,den,Ts) vẽ ra giản đồ Bode của hàm truyền đa thức hệ liên tục gián đoạn. G(z) = num(z)/den(z) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. dbode(a,b,c,d,Ts,iu,w) hay dbode(num,den,Ts,w) vẽ ra giản đồ Bode với vector tần số w do người sử dụng xác định. Vector w chỉ ra các điểm tần số (tính bằng rad/s) mà tại đó đáp ứng tần số giản đồ Bode được tính. Hiện tượng trùng phổ xảy ra tại tần số lớn hơn tần số Nyquist. Nếu vẫn giữ lại các đối số ở vế trái của dòng lệnh thì: [mag,phase,w] = dbode(a,b,c,d,Ts) [mag,phase,w] = dbode(a,b,c,d,Ts,iu) [mag,phase,w] = bode(a,b,c,d,Ts,iu,w) [mag,phase,w] = bode(num,den,Ts) [mag,phase,w] = bode(num,den,Ts,w) sẽ không vẽ ra giản đồ Bode mà tạo ra các ma trận đáp ứng tần số mag, phase và w của hệ thống được tính tại các giá trị tần số w. Ma trận mag và phase có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector w. G(z) = C(zI –A)-1B + D mag(ω) = ⏐G(ejωT)⏐ phase(ω) = ∠G(ejωT) trong đó T là thời gian lấy mẫu. Góc pha được tính bằng độ. Giá trị biên độ có thể chuyển thành decibel theo biểu thức: magdB = 20*log10(mag) d) Ví dụ: Vẽ đáp ứng giản đồ Bode của hệ thống có hàm truyền như sau: 2z2 − 3.4 z + 1.5 H() z = z2 −1.6 s + 0.8 với thời gian lấy mẫu Ts = 0.1 num = [2 -3.4 1.5]; den = [1 -1.6 0.8]; dbode(num,den,0.1); grid on và ta được đáp ứng tần số giản đồ Bode của hệ gián đoạn như sau:
  7. MATLAB trong điều khiển tự động Trang 164 Bode Diagrams 20 10 0 -10 100 50 Phase (deg); Magnitude (dB) Magnitude (deg); Phase 0 -50 -1 0 1 2 10 10 10 10 Frequency (rad/sec) 4. Lệnh FREQS a) Công dụng: Tìm đáp ứng tần số của phép biến đổi Laplace. b) Cú pháp: h = freqs(b,a,w) [h,w] = freqs(b,a) [h,w] = freqs(b,a,n) freqs(b,a) c) Giải thích: Lệnh freqs trở thành đáp ứng tần số H(jω) của bộ lọc analog. B() s b(1) snb + b (2) s nb−1 + +b ( nb + 1) H() s = = A() s a(1) sna + a (2) s na−1 + +a ( na + 1) trong đó vector b và a chứa các hệ số của tử số và mẫu số. h = freqs(b,a,w) tạo ra vector đáp ứng tần số phức của bộ lọc analog được chỉ định bởi các hệ số trong vector b và a. Lệnh freqs tìm đáp ứng tần số trong mặt phẳng phức tại các thời điểm tần số được hcỉ định trong vector w. [h,w] = freqs(b,a) tự động chọn 200 điểm tần số trong vector w để tính vector đáp ứng tần số h. [h,w] = freqs(b,a,n) chọn ra n điểm tần số để tìm vector đáp ứng tần số h. Nếu bỏ qua các đối số ngõ ra ở vế trái thì lệnh freqs sẽ vẽ ra đáp ứng biên độ và pha trên màn hình. freqs chỉ dùng cho các hệ thống có ngõ vào thực và tần số dương. d) Ví dụ:
  8. MATLAB trong điều khiển tự động Trang 165 Tìm và vẽ đáp ứng tần số của hệ thống có hàm truyền: 0.2s2 + 0.3 s + 1 H() s = s2 +0.4 s + 1 % Khai báo hàm truyền: a = [1 0.4 1]; b = [0.2 0.3 1]; % Xác định trục tần số: w = logspace(-1,1); % Thực hiện vẽ đồ thị: freqs(b,a,w) 1 10 0 10 Magnitude -1 10 -1 0 1 10 10 10 Frequency (radians) 0 -50 -100 Phase (degrees) Phase -150 -1 0 1 10 10 10 Frequency (radians) 5. Lệnh FREQZ a) Công dụng: Tìm đáp ứng tần số của bộ lọc số. b) Cú pháp: [h,w] = freqz(b,a,n) [h,f] = freqz(b,a,n,Fs) [h,w] = freqz(b,a,n,‘whole’) [h,f] = freqz(b,a,n,‘whole’,Fs) h = freqz(b,a,w) h = freqz(b,a,f,Fs) freqz(b,a) c) Giải thích:
  9. MATLAB trong điều khiển tự động Trang 166 Lệnh freqz tìm đáp ứng tần số H(ejωT) của bộ lọc số từ các hệ số tử số và mẫu số trong vector b và a. [h,w] = freqz(b,a,n) tìm đáp ứng tần số của bộ lọc số với n điểm B() z b(1)+ b (2) z−1 + + b ( nb + 1) z −nb H() z = = A() z a(1)+ a (2) z−1 + + a ( na + 1) z −na từ các hệ số trong vector b và a. freqz tạo ra vector đáp ứng tần số hồi tiếp và vector w chứa n điểm tần số. freqz xác định đáp ứng tần số tại n điểm nằm đều nhau quanh nửa vòng tròn đơn vị, vì vậy w chứa n điểm giữa 0 và π. [h,f] = freqz(b,a,n,Fs) chỉ ra tần số lấy mẫu dương Fs (tính bằng Hz). Nó tạo ra vector f chứa các điểm tần số thực giữa 0 và Fs/2 mà tại đó lệng sẽ tính đáp ứng tần số. [h,w] = freqz(b,a,n,‘whole’) và [h,f] = freqz(b,a,n,‘whole’,Fs) sử dụng nđiểm quanh vòng tròn đơn vị (từ 0 tới 2π hoặc từ 0 tới Fs) h = freqz(b,a,w) tạo ra đáp ứng tần số tại các điểm tần số được chỉ trong vector w. Các điểm tần số này phải nằm trong khoảng (0 ÷2π). h = freqz(b,a,f,Fs) tạo ra đáp ứng tần số tại các điểm tần số được chỉ trong vector f. Các điểm tần số này phải nằm trong khoảng (0 ÷ Fs). Nếu bỏ qua các đối số ngõ ra thì lệnh freqz vẽ ra các đáp ứng biên độ và pha trên màn hình. Lệnh freqz dùng cho các hệ thống có ngõ vào thực hoặc phức. d) Ví dụ: Vẽ đáp ứng biên độ và pha của bộ lọc Butter. [b,a] = butter(5,0.2); freqz(b,a,128) và ta được đồ thị đáp ứng: 100 0 -100 -200 -300 Magnitude Response (dB) Response Magnitude 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Normalized frequency (Nyquist == 1) 0 -100 -200 -300 -400 Phase (degrees) Phase -500 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Normalized frequency (Nyquist == 1)
  10. MATLAB trong điều khiển tự động Trang 167 6. Lệnh NYQUIST a) Công dụng: Vẽ biểu đồ đáp ứng tần số Nyquist. b) Cú pháp: [re,im,w] = nyquist(a,b,c,d) [re,im,w] = nyquist(a,b,c,d,iu) [re,im,w] = nyquist(a,b,c,d,iu,w) [re,im,w] = nyquist(num,den) [re,im,w] = nyquist(num,den,w) c) Giải thích: Lệnh nyquist tìm đáp ừng tần số Nyquist của hệ liên tục LTI. Biểu đồ Nyquist dùng để phân tích đặc điểm của hệ thống bao gồm: biên dự trữ, pha dự trữ và tính ổn định. Nều bỏ qua các đối số ở vế trái của dòng lệnh thì nyquist sẽ vẽ ra biểu đồ Nyquist trên màn hình. Lệnh nyquist có thể xác định tính ổn định của hệ thống hồi tiếp đơn vị. Cho biểu đồ Nyquist của hàm truyền vòng hở G(s), hàm truyền vòng kín: G() s Gcl (s) = 1+ G ( s ) là ổn định khi biểu đồ Nyquist bao quanh điểm –1+j0 P lần theo chiều kim đồng hồ, trong đó P là số cực vòng hở không ổn định. nyquist(a,b,c,d) vẽ ra chuỗi biểu đồ Nyquist, mỗi đồ thị ứng vời mối quan hệ giữa một ngõ vào và một ngõ ra của hệ không gian trạng thái liên tục: . x= Ax + Bu y = Cx + Du với trục tần số được xác định tự động. Nếu đáp ứng thay đổi càng nhanh thì cần phải xác định càng nhiều điểm trên trục tần số. nyquist(a,b,c,d,iu) vẽ ra biểu đồ Nyquist từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng Nyquist. nyquist(num,den) vẽ ra biểu đồ Nyquist của hàm truyền đa thức hệ liên tục G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. nyquist(a,b,c,d,iu,w) hoặc nyquist(num,den,w) vẽ ra biểu đồ Nyquist với vector tần số w do người sử dụng xác định. Vector w chỉ ra các điểm tần số (tính bằng rad/s) mà tại đó đáp ứng Nyquist được tính. Nếu vẫn giữ lại các đối số ở vế trái của dòng lệnh thì: [re,im,w] = nyquist(a,b,c,d) [re,im,w] = nyquist(a,b,c,d,iu) [re,im,w] = nyquist(a,b,c,d,iu,w) [re,im,w] = nyquist(num,den) [re,im,w] = nyquist(num,den,w)
  11. MATLAB trong điều khiển tự động Trang 168 không vẽ ra biểu đồ Nyquist mà tạo ra đáp ứng tần số của hệ thống dưới dạng các ma trận re, im và w. Các ma trận re và im có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector w. d) Ví dụ: Vẽ biểu đồ Nyquist của hệ thống có hàm truyền: 2s2 + 5 s + 1 H() s = s2 +2 s + 3 num = [2 5 1]; den = [1 2 3]; nyquist(num,den); title(‘Bieu do Nyquist’) và ta được biểu đồ Nyquist như hình vẽ: 7. Lệnh DNYQUIST a) Công dụng: Vẽ biểu đồ đáp ứng tần số Nyquist của hệ gián đoạn. b) Cú pháp: [re,im,w] = dnyquist(a,b,c,d,Ts) [re,im,w] = dnyquist(a,b,c,d,Ts,iu) [re,im,w] = dnyquist(a,b,c,d,Ts,iu,w) [re,im,w] = dnyquist(num,den,Ts) [re,im,w] = dnyquist(num,den,Ts,w) c) Giải thích:
  12. MATLAB trong điều khiển tự động Trang 169 Lệnh dnyquist tìm đáp ừng tần số Nyquist của hệ gián đoạn LTI. Biểu đồ Nyquist dùng để phân tích đặc điểm của hệ thống bao gồm: biên dự trữ, pha dự trữ và tính ổn định. Đáp ứng tần số dùng lệnh dnyquist có thể so sánh trực tiếp với đáp ứng nyquist của hệ liên tục tương ứng. Nều bỏ qua các đối số ở vế trái của dòng lệnh thì dnyquist sẽ vẽ ra biểu đồ Nyquist trên màn hình. Lệnh dnyquist có thể xác định tính ổn định của hệ thống hồi tiếp đơn vị. Cho biểu đồ Nyquist của hàm truyền vòng hở G(s), hàm truyền vòng kín: G() z Gcl (z) = 1+G ( z ) là ổn định khi biểu đồ Nyquist bao quanh điểm –1+j0 P lần theo chiều kim đồng hồ, trong đó P là số cực vòng hở không ổn định. dnyquist(a,b,c,d,Ts) vẽ ra chuỗi biểu đồ Nyquist, mỗi đồ thị ứng vời mối quan hệ giữa một ngõ vào và một ngõ ra của hệ không gian trạng thái gián đoạn: x[n+] = Ax[n] + Bu{n] y[n] = Cx[n] + Du[n] với trục tần số được xác định tự động. Các điểm tần số được chọn trong khoảng từ 0 đến π/Ts radians tương ứng với nửa tần số lấy mẫu (tần số Nyquist). Nếu đáp ứng thay đổi càng nhanh thì cần phải xác định càng nhiều điểm trên trục tần số. Tần số là thời gian lấy mẫu. dnyquist(a,b,c,d,Ts,iu) vẽ ra biểu đồ Nyquist từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng Nyquist. dnyquist(num,den,Ts) vẽ ra biểu đồ Nyquist của hàm truyền đa thức hệ gián đoạn: G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. dnyquist(a,b,c,d,Ts,iu,w) hoặc dnyquist(num,den,w) vẽ ra biểu đồ Nyquist với vector tần số w do người sử dụng xác định. Vector w chỉ ra các điểm tần số (tính bằng rad/s) mà tại đó đáp ứng Nyquist được tính. Hiện tượng trùng phổ xảy ra tại tần số lớn hơn tần số Nyquist (π/Ts rad/s). Để tạo ra trục tần số với các khoảng tần số bằng nhau theo logarit ta dùng lệnh logspace. Nếu vẫn giữ lại các đối số ở vế trái của dòng lệnh thì: [re,im,w] = dnyquist(a,b,c,d,Ts) [re,im,w] = dnyquist(a,b,c,d,Ts,iu) [re,im,w] = dnyquist(a,b,c,d,Ts,iu,w) [re,im,w] = dnyquist(num,den,Ts) [re,im,w] = dnyquist(num,den,Ts,w) không vẽ ra biểu đồ Nyquist mà tạo ra đáp ứng tần số của hệ thống dưới dạng các ma trận re, im và w. Các ma trận re và im chứa các phần thực và phần ảo của đáp ứng tần số của hệ thống được tính tại các giá trị tần số w, re và im có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector w. d) Ví dụ: Vẽ biểu đồ Nyquist của hệ gián đoạn có hàm truyền: 2z2 − 3.4 z + 1.5 H() z = z2 −1.6 z + 0.8 với thời gian lấy mẫu Ts = 0.1
  13. MATLAB trong điều khiển tự động Trang 170 % Xác định hàm truyền: num = [2 -3.4 1.5]; den = [1 -1.6 0.8]; % Vẽ biểu đồ Nyquist: dnyquist(num,den,0.1) title(‘Bieu do Nyquist he gian doan’) và ta được biểu đồ Nyquist hệ gián đoạn như sau: 8. Lệnh NICHOLS a) Công dụng: Vẽ biểu đồ đáp ứng tần số Nichols. b) Cú pháp: [mag,phase,w] = nichols(a,b,c,d) [mag,phase,w] = nichols(a,b,c,d,iu) [mag,phase,w] = nichols(a,b,c,d,iu,w) [mag,phase,w] = nichols(num,den) [mag,phase,w] = nichols(num,den,w) c) Giải thích: Lệnh nichols tìm đáp ứng tần số Nichols của hệ liên tục LTI. Biểu đồ Nichols được dùng để phân tích đặc điểm của hệ vòng hở và hệ vòng kín. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh nichols sẽ vẽ ra biểu đồ Nichols trên màn hình.
  14. MATLAB trong điều khiển tự động Trang 171 nichols(a,b,c,d) vẽ ra chuỗi biểu đồ Nichols, mỗi đồ thị tương ứng với mối quan hệ giữa một ngõ vào và một ngõ ra của hệ không gian trạng thái liên tục: . x= Ax + Bu y = Cx + Du với trục tần số được xác định tự động. Nếu đáp ứng thay đổi nhanh thì cần phải xác định càng nhiều điểm trên trục tần số. nichols(a,b,c,d,iu) vẽ ra biểu đồ Nichols từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng Nichols. nichols(num,den) vẽ ra biểu đồ Nichols của hàm truyền đa thức hệ liên tục G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. nichols(a,b,c,d,iu,w) hay nichols(num,den,w) vẽ ra biểu đồ Nichols với vector tần số w do người sử dụng xác định. Vector w chỉ định những điểm tần số (tính bằng rad/s) mà tại đó đáp ứng Nichols được tính. Để tạo ra trục tần số với các khoảng tần số bằng nhau theo logarit ta dùng lệnh logspace. Nếu giữ lại các đối số ở vế trái của dòng lệnh thì: [mag,phase,w] = nichols(a,b,c,d) [mag,phase,w] = nichols(a,b,c,d,iu) [mag,phase,w] = nichols(a,b,c,d,iu,w) [mag,phase,w] = nichols(num,den) [mag,phase,w] = nichols(num,den,w) sẽ không vẽ ra biểu đồ Nichols mà tạo ra đáp ứng tần số của hệ thống dưới dạng các ma trận mag, phase và w. Các ma trận mag và phase chứa đáp ứng biên độ và pha của hệ thống được xác định tại những điểm tần số w. Ma trận mag và phase có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector w. G(s) = C(sI –A)-1B + D mag(ω) = ⏐G(jω)⏐ phase(ω) = ∠G(jω) Góc pha được tính bằng độ và nằm trong khoảng –3600 tới 00. Giá trị biên độ có thể chuyển về đơn vị decibel theo công thức: magdB = 20*log10(mag) Để vẽ lưới biểu đồ Nichols ta dùng lệnh ngrid. d) Ví dụ: Trích trang 11-150 sách ‘Control System Toolbox’ Vẽ đáp ứng Nichols của hệ thống có hàm truyền: −4s4 + 48 s3 − 18 s2 + 250 s + 600 H() s = s4 +30 s3 + 282 s2 +525 s + 60 num = [-4 48 -18 250 600]; den = [1 30 282 525 60]; nichols(num,den) title(‘Bieu do Nichols’) ngrid(‘new’)
  15. MATLAB trong điều khiển tự động Trang 172 và ta được biểu đồ Nichols như hình vẽ: 9. Lệnh DNICHOLS a) Công dụng: Vẽ biểu đồ đáp ứng tần số Nichols của hệ gián đoạn. b) Cú pháp: [mag,phase,w] = dnichols(a,b,c,d,Ts) [mag,phase,w] = dnichols(a,b,c,d,Ts,iu) [mag,phase,w] = dnichols(a,b,c,d,Ts,iu,w) [mag,phase,w] = dnichols(num,den,Ts) [mag,phase,w] = dnichols(num,den,Ts,w) c) Giải thích: Lệnh dnichols tìm đáp ứng tần số Nichols của hệ gián đoạn LTI. Biểu đồ Nichols được dùng để phân tích đặc điểm của hệ vòng hở và hệ vòng kín. Đáp ứng từ lệnh dnichols có thể so sánh trực tiếp với đáp ứng từ lệnh nichols của hệ liên tục tương ứng. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh dnichols sẽ vẽ ra biểu đồ Nichols trên màn hình. dnichols(a,b,c,d,Ts) vẽ ra chuỗi biểu đồ Nichols, mỗi đồ thị tương ứng với mối quan hệ giữa một ngõ vào và một ngõ ra của hệ không gian trạng thái gián đoạn: x[n+] = Ax[n] + Bu{n] y[n] = Cx[n] + Du[n] với trục tần số được xác định tự động. Các điểm tần số được chọn trong khoảng từ 0 tới π/Ts radians. Nếu đáp ứng thay đổi nhanh thì cần phải xác định càng nhiều điểm trên trục tần số.
  16. MATLAB trong điều khiển tự động Trang 173 dnichols(a,b,c,d,Ts,iu) vẽ ra biểu đồ Nichols trên màn hình từ ngõ vào duy nhất iu tới tất cả các ngõ ra của hệ thống với trục tần số được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được sử dụng cho đáp ứng Nichols. dnichols(num,den,Ts) vẽ ra biểu đồ Nichols của hàm truyền đa thức hệ gián đoạn G(z) = num(z)/den(z) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. dnichols(a,b,c,d,Ts,iu,w) hay dnichols(num,den,Ts,w) vẽ ra biểu đồ Nichols với vector tần số w do người sử dụng xác định. Vector w chỉ định những điểm tần số (tính bằng rad/s) mà tại đó đáp ứng Nichols được tính. Hiện tượng trùng phổ xảy ra tại tần số lớn hơn tần số Nyquist (π/Ts rad/s). Để tạo ra trục tần số với các khoảng tần số bằng nhau theo logarit ta dùng lệnh logspace. Nếu giữ lại các đối số ở vế trái của dòng lệnh thì: [mag,phase,w] = dnichols(a,b,c,d,Ts) [mag,phase,w] = dnichols(a,b,c,d,Ts,iu) [mag,phase,w] = dnichols(a,b,c,d,Ts,iu,w) [mag,phase,w] = dnichols(num,den,Ts) [mag,phase,w] = dnichols(num,den,Ts,w) không vẽ ra biểu đồ Nichols mà tạo ra đáp ứng tần số của hệ thống dưới dạng các ma trận mag, phase và w. Các ma trận mag và phase chứa đáp ứng biên độ và pha của hệ thống được xác định tại những điểm tần số w. Ma trận mag và phase có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector w. G(z) = C(zI –A)-1B + D mag(ω) = ⏐G(ejωT)⏐ phase(ω) = ∠G(ejωT) trong đó T là thời gian lấy mẫu. Góc pha được tính bằng độ và nằm trong khoảng –3600 tới 00. Giá trị biên độ có thể chuyển về đơn vị decibel theo công thức: magdB = 20*log10(mag) Để vẽ lưới biểu đồ Nichols ta dùng lệnh ngrid. d) Ví dụ: Vẽ đáp ứng Nichols của hệ thống có hàm truyền: 5.1 H() z = z4 +1.1 z3 + 1.36 z2 + 0.88 z + 0.31 num = 1.5; den = [1 1.1 1.36 0.88 0.31]; ngrid(‘new’) dnichols(num,den,0.05) title(‘Bieu do Nichols gian doan’) và ta được biểu đồ Nichols của hệ gián đoạn:
  17. MATLAB trong điều khiển tự động Trang 174 10. Lệnh NGRID a) Công dụng: Tạo lưới cho đồ thị Nichols. b) Cú pháp: ngrid ngrid(‘new’) c) Giải thích: Lệnh grid tạo lưới cho đồ thị Nichols. Đồ thị này có liên hệ với số phức H/(1+H), trong đó H là một số phức bất kỳ. Nếu H là một điểm trên đáp ứng tần số vòng hở của hệ SISO thì H/(1+H) là giá trị tương ứng trên đáp ứng tần số vòng kín của hệ thống. 0 0 ngrid tạo ra lưới trong vùng có biên độ từ –40 dB tới 40 dB và góc pha từ -360 tới 0 với các đường hằng số mag(H/(1+H)) và angle(H/(1+H)) được vẽ. ngrid vẽ lưới đồ thị Nichols ngoài biểu đồ Nichols đã có như biểu đồ được tạo ra bởi lệnh nichols hoặc dnichols. ngrid(‘new’) xóa màn hình đồ họa trước khi vẽ lưới và thiết lập trạng thái giữ để đáp ứng Nichols có thể được vẽ bằng cách dùng lệnh: ngrid(‘new’) nichols(num,den) hay nichols(a,b,c,d,iu) d) Ví dụ: Vẽ lưới trên biểu đồ Nichols của hệ thống: −4s4 + 48 s3 − 18 s2 + 250 s + 600 H() s = s4 +30 s3 + 282 s2 +525 s + 60 num = [-4 48 -18 250 600];
  18. MATLAB trong điều khiển tự động Trang 175 den = [1 30 282 525 60]; nichols(num,den) title(‘Bieu do Nichols’) ngrid(‘new’) và ta được đồ thị đáp ứng như sau: 11. Lệnh MARGIN a) Công dụng: Tính biên dự trữ và pha dự trữ. b) Cú pháp: [Gm,Pm,Wcp,Wcg] = margin(mag,phase,w) [Gm,Pm,Wcp,Wcg] = margin(num,den) [Gm,Pm,Wcp,Wcg] = margin(a,b,c,d) c) Giải thích: Lệnh margin tính biên dự trữ (gain margin), pha dự trữ (phase margin) và tần số cắt (crossover frequency) từ dữ liệu đáp ứng tần số. Biên dự trữ và pha dự trữ dựa trên hệ thống vòng hở SISO và cho biết tính ổn định tương đối của hệ thống khi hệ thống là hệ thống vòng kín. Nếu bỏ qua các đối số ở vế trái dòng lệnh thì giản đồ Bode với biên dự trữ và pha dự trữ sẽ được vẽ trên màn hình. Biên dự trữ là độ lợi cần tăng thêm để tạo ra độ lợi vòng đơn vị tại tần số mà góc pha bằng –1800. Nói cách khác, biên dự trữ là 1/g nếu g là độ lợi tại tần sồ góc pha –1800. Tương tự, pha dự trữ là sự khác biệt giữa góc pha đáp ứng và –1800 khi độ lợi là 1. Tần số mà tại đó biên độ là 1 được gọi là tần số độ lợi đơn vị (unity-gain frequency) hoặc tần số cắt.
  19. MATLAB trong điều khiển tự động Trang 176 margin(num,den) tính biên dự trữ và pha dự trữ của hàm truyền liên tục: G(s) = num/den Tương tự, margin(a,b,c,d) tính độ dự trữ của hệ không gian trạng thái (a,b,c,d). Với cách này, lệnh margin chỉ sử dụng cho hệ liên tục. Đối với hệ gián đoạn, ta sử dụng lệnh dbode để tìm đáp ứng tần số rồi gọi margin. [mag,phase,w] = dbode(a,b,c,d,Ts) margin(mag,phase,w) [Gm,Pm,Wcp,Wcg] = margin(mag,phase,w) sẽ không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận biên dự trữ Gm, pha dự trữ Pm, tần số kết hợp Wcp, Wcg được cho bởi các vector biên độ mag, phase và tần số w của hệ thống. Các giá trị chính xác được tìm ra bằng cách dùng phép nội suy giữa các điểm tần số. Góc pha được tính bằng độ. d) Ví dụ: Tìm biên dự trữ, pha dự trữ và vẽ giản đồ Bode của hệ bậc 2 có ωn = 1 và ζ = 0.2 [a,b,c,d] = ord(1,0.2); bode(a,b,c,d) margin(a,b,c,d) [Gm,Pm,Wcp,Wcg] = margin(a,b,c,d) và ta được kết quả: Gm = lnf(∞) Pm = 32.8599 độ Wcg = NaN (không xác định) Wcp = 1.3565 Giản đồ Bode của hệ: 12. Lệnh SIGMA a) Công dụng:
  20. MATLAB trong điều khiển tự động Trang 177 Tìm giản đồ Bode giá trị suy biến của hệ không gian trạng thái. b) Cú pháp: [sv,w] = sigma(a,b,c,d) [sv,w] = sigma(a,b,c,d,‘inv’) [sv,w] = sigma(a,b,c,d,w) [sv,w] = sigma(a,b,c,d,w,‘inv’) c) Giải thích: Lệnh sigma tính các giá trị suy biến của ma trận phức C(jωI-A)-1B+D theo hàm của tần số ω. Các giá trị suy biến là mở rộng của đáp ứng biên độ giản đồ Bode của hệ MIMO. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì sigma sẽ vẽ ra giản đồ Bode của giá trị suy biến trên màn hình. [sv,w] = sigma(a,b,c,d) vẽ ra giản đồ suy biến của ma trận phức: G(w) = C(jωI-A)-1B+D theo hàm của tần số. Trục tần số được chọn tự động và phối hợp nhiều điểm nếu đồ thị thay điểm nhanh. Đối với các ma trận vuông, sigma(a,b,c,d,‘inv’) vẽ đồ thị các giá trị suy biến của ma trận phức đảo: G-1(w) = [C(jωI-A)-1B+D]-1 sigma(a,b,c,d,w) hoặc sigma(a,b,c,d,w,‘inv’) vẽ đồ thị các giá trị suy biến với vector tần số do người sử dụng xác định. Vector w chỉ ra những tần số (tính bằng rad/s) mà tại đó đáp ứng các giá trị suy biến được tính. Nếu giữ lại các đối số ở vế trái dòng lệnh thì: [sv,w] = sigma(a,b,c,d) [sv,w] = sigma(a,b,c,d,‘inv’) [sv,w] = sigma(a,b,c,d,w) [sv,w] = sigma(a,b,c,d,w,‘inv’) không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận suy biến theo chiều giảm dần của bậc tương ứng với các điểm tần số trong vector w. Đối với phép phân tích rắn chắc, các giá trị suy biến của ma trận hàm truyền đặc biệt được phân tích. Về thực hiện các lệnh để đạt được ma trận hàm truyền mong muốn của một số khối được trình bày trong bảng sau: Ma trận hàm truyền Sơ đồ khối Lệnh G(jω) G(s) sigma(a,b,c,d) G-1(s) G-1(jω) sigma(a,b,c,d,‘inv’)
  21. G(s) MATLAB trong điều khiển tự động Trang 178 [a,b,c,d] = parallel(a,b,c,d,[ ],[ ],[ ],eye(d)) G(s) sigma(a,b,c,d) 1+G(jω) [a,b,c,d] = feedback([ ],[ ],[ ],eye(d),a,b,c,d) G-1(s) sigma(a,b,c,d,‘inv’) [a,b,c,d] = feedback(a,b,c,d,[ ],[ ],[ ],eye(d)) 1+G-1(jω) sigma(a,b,c,d) Đáp ứng giá trị suy biến của hệ SISO tương đương với đáp ứng biên độ giản đồ Bode của hệ đó. d) Ví dụ: Xét hệ bậc 2 có ωn = 1 và ζ = 0.2. Vẽ đồ thị giá trị suy biến của hệ thống. [a,b,c,d] = ord(1,0.2); margin(a,b,c,d) title(‘Gia tri suy bien’) và ta được đáp ứng như hình vẽ: 13. Lệnh DSIGMA a) Công dụng: Tìm giản đồ Bode giá trị suy biến của hệ không gian trạng thái. b) Cú pháp: [sv,w]= dsigma(a,b,c,d,Ts)
  22. MATLAB trong điều khiển tự động Trang 179 [sv,w]= dsigma(a,b,c,d,Ts,‘inv’) [sv,w]= dsigma(a,b,c,d,Ts,w) [sv,w]= dsigma(a,b,c,d,Ts,w,'inv') c) Giải thích: Lệnh dsigma tính các giá trị suy biến của ma trận phức C(ejωTI-A)-1+B+D theo hàm của tần số ω. Các gia trị suy biến là mở rộng của đáp ứng biên độ giản đồ Bode của hệ MIMO và có thể được dùng để xác định độ rắn chắc của hệ thống. Nếu bỏ qua các đối số ở vế trái dòng lệnh thì dsigma sẽ vẽ ra giản đồ Bode của giá trị suy biến trên màn hình. dsigma(a,b,c,d,Ts) vẽ giản đồ suy biến của ma trận phức : G(w) = C(ejωTI-A)-1+B+D theo hàm của tần số. Các điểm tần số được chọn tự động trong khoảng từ 0 tới π/Ts rad/sec trong đó π/Ts rad/sec tương ứng với nửa tần số lấy mẫu (tần số Nyquist). Nếu đồ thị thay đổi nhanh thì cần chọn nhiều điểm tần số hơn. Đối với các hệ thống có ma trận vuông, dsigma(a,b,c,d,Ts,‘inv’) vẽ đồ thị các giá trị suy biến của ma trận phức đảo : G-1(w) = [C(ejωTI-A)-1B+D]-1 dsigma(a,b,c,d,Ts,w) hoặc dsigma(a,b,c,d,Ts,‘inv’) vẽ đồ thị các giá trị suy biến với vector tần số do người sử dụng xác định. Vector w chỉ ra những tần số (tính bằng rad/sec) mà tại đó đáp ứng các giá trị suy biến được tính. Hiện tượng trùng phổ xảy ra tại tần số lớn hơn tần số Nyquist (π/Ts rad/sec). Để tạo ra vector tần số được chia đều theo logarit tần số ta dùng lệnh logspace. Nếu giữ lại các đối số ở vế trái dòng lệnh thì : [sv,w]= dsigma(a,b,c,d,Ts) [sv,w]= dsigma(a,b,c,d,Ts,‘inv’) [sv,w]= dsigma(a,b,c,d,Ts,w) [sv,w]= dsigma(a,b,c,d,Ts,w,‘inv’) không vẽ ra các đồ thị đáp ứng mà tạo ra các giá trị suy biến trong sv và các điểm tần số w. Mỗi hàng của ma trận sv chứa các giá trị suy biến theo chiều giảm dần của bậc tương ứng với các điểm tần số trong vector w. Đối với phép phân tích rắn chắc, các giá trị suy biến của ma trận hàm truyền đặc biệt được phân tích. Việc thực hiện các lệnh để đạt được ma trận hàm truyền mong muốn của một số khối được trình bày trong bảng sau : Ma trận hàm truyền Sơ đồ khối Lệnh G(s) G(jω) dsigma(a,b,c,d) G-1(s) G-1(jω) dsigma(a,b,c,d, ‘inv’) [a,b,c,d]= parallel(a,b,c,d,[ ],[ ],[ ],eye(d)) G(s) dsigma(a,b,c,d) G(s)
  23. MATLAB trong điều khiển tự động Trang 180 1+ G(jω) [a,b,c,d]=feedback([ ],[ ],[ ],eye(d),a,b,c,d) dsigma(a,b,c,d,‘inv’) [a,b,c,d]= feedback(a,b,c,d,[ ],[ ],[ ],eye(d)) G-1(s) 1+G-1(jω) dsigma(a,b,c,d) Đáp ứng giá trị suy biến của hệ SISO tương đương với đáp ứng biên độ giản đồ Bode của hệ đó. d) Ví dụ: Xét hệ bậc 2 có ωn = 1 và ζ = 0.2. Vẽ đồ thị giá trị suy biến của hệ thống với thời gian lấy mẫu Ts = 0.1 [a,b,c,d]= ord2(1,0.2); bode(a,b,c,d) dsigma(a,b,c,d,0.1) title('Gia tri suy bien gian doan') và ta có giản đồ Bode giá trị suy biến : 14. Lệnh LTIFR a) Công dụng: Đáp ứng tần số của hệ tuyến tính bất biến. b) Cú pháp: ltifr(a,b,s)
  24. MATLAB trong điều khiển tự động Trang 181 c) Giải thích: Lệnh ltifr dùng để mở rộng đáp ứng tần số của hệ không gian trạng thái tuyến tính bất biến. G = Ltifr(a,b,s) tìm đáp ứng tần số của hệ thống với một ngõ vào duy nhất : G(s) = (sI – A)-1B Vector s chỉ ra số phức mà tại đó đáp ứng tần số được xác định. Đối với đáp ứng giản đồ Bode hệ liên tục, s nằm trên trục ảo. Đối với đáp ứng giản đồ Bode hệ gián đoạn, s nhận các giá trị quanh vòng tròn đơn vị. ltifr tạo ra đáp ứng tần số dưới dạng ma trận phức G với số cột bằng số trạng thái hay số hàng của ma trận A và có số hàng là length(s). CÁC BÀI TẬP VỀ ĐÁP ỨNG TẦN SỐ Bài 1: hàm margin (bài tập này trích từ trang 11-138 sách ‘Control System Toollbox’ » hd=tf([0.04798 0.0464],[1 -1.81 0.9048],0.1)
  25. MATLAB trong điều khiển tự động Trang 182 Transfer function: 0.04798 z + 0.0464 z^2 - 1.81 z + 0.9048 Sampling time: 0.1 ; Thời gian lấy mẫu: 0,1 » [Gm,Pm,Wcg,Wcp]=margin(hd); » [Gm,Pm,Wcg,Wcp] ans = 2.0517 13.5712 5.4374 4.3544 » margin(hd) Kết quả: Bode Diagrams Gm=6.2424 dB (at 5.4374 rad/sec), Pm=13.571 deg. (at 4.3544 rad/sec) 20 0 -20 -40 -60 -80 0 -100 Phase Magnitude (dB) (deg); -200 -300 101 Frequency (rad/sec) Bài 2: lệnh modred (bài tập này trích từ trang 11-142 sách ‘Control System Toollbox’ s3+ 11 s 2+ 36 s + 26 h() s = s4+ 14,6 s 3 + 74,96 s 2 + 153,7 s + 99,65 » h=tf([1 11 36 26],[1 14.6 74.96 153.7 99.65]) Transfer function: s^3 + 11 s^2 + 36 s + 26
  26. MATLAB trong điều khiển tự động Trang 183 s^4 + 14.6 s^3 + 74.96 s^2 + 153.7 s + 99.65 » [hb,g]=balreal(h) a = x1 x2 x3 x4 x1 -3.6014 -0.82121 -0.61634 -0.058315 x2 0.82121 -0.59297 -1.0273 -0.090334 x3 -0.61634 1.0273 -5.9138 -1.1272 x4 0.058315 -0.090334 1.1272 -4.4918 b = u1 x1 1.002 x2 -0.10641 x3 0.086124 x4 -0.0081117 c = x1 x2 x3 x4 y1 1.002 0.10641 0.086124 0.0081117 d = u1 y1 0 Continuous-time model. g = 0.1394 0.0095 0.0006 0.0000 » g' ans = 0.1394 0.0095 0.0006 0.0000 » hmdc=modred(hb,2:4,'mdc') a = x1 x1 -4.6552 b = u1 x1 1.1392 c = x1 y1 1.1392
  27. MATLAB trong điều khiển tự động Trang 184 d = u1 y1 -0.017857 Continuous-time model. » hdel=modred(hb,2:4,'del') a = x1 x1 -3.6014 b = u1 x1 1.002 c = x1 y1 1.002 d = u1 y1 0 Continuous-time model. » bode(h,'-',hmdc,'x',hdel,'*') Kết quả: Bode Diagrams From: U(1) 0 -20 -40 -60 -80 0 -50 Phase Magnitude (dB) (deg); -100 To: Y (1) -150 -200 10-1 100 101 102 103 Frequency (rad/sec) Bài 3: (Trang 11-16 sách ‘Control System Toollbox’) Xem zero-pole-gain (zero-cực-độ lợi) của hệ thống sau: » sys=zpk([-10 -20.01],[-5 -9.9 -20.1],1) Zero/pole/gain: (s+10) (s+20.01)
  28. MATLAB trong điều khiển tự động Trang 185 (s+5) (s+9.9) (s+20.1) » » [sys,g]=balreal(sys) a = x1 x2 x3 x1 -4.9697 0.2399 -0.22617 x2 -0.2399 -4.2756 9.4671 x3 -0.22617 -9.4671 -25.755 b = u1 x1 1 x2 0.024121 x3 0.022758 c = x1 x2 x3 y1 1 -0.024121 0.022758 d = u1 y1 0 Continuous-time model. g = 0.1006 0.0001 0.0000 » g' ans = 0.1006 0.0001 0.0000 » sysr=modred(sys,[2 3],'del') a = x1 x1 -4.9697 b = u1 x1 1 c = x1 y1 1 d = u1 y1 0
  29. MATLAB trong điều khiển tự động Trang 186 Continuous-time model. » zpk(sysr) Zero/pole/gain: 1.0001 (s+4.97) » bode(sys,'-',sysr,'x') Bode Diagrams From: U(1) -10 -20 -30 -40 -50 0 -20 -40 Phase Magnitude (dB) (deg); -60 To: Y (1) -80 -100 100 101 102 Frequency (rad/sec) Bài 4: Trích từ trang 55 sách ‘Hướng dẫn sử dụng MATLAB’ tác giả Nguyễn Văn Giáp. Vẽ biểu đồ nyquist của hệ thống: H(s) = (s+4)/(s2 + 3s – 8) » num=[1 4]; » den=[1 3 -8]; » nyquist(num,den);
  30. MATLAB trong điều khiển tự động Trang 187 Nyquist Diagrams From: U(1) 0.3 0.2 0.1 0 To: Y (1) -0.1 Imaginary AxisImaginary -0.2 -0.3 -0.4 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 Real Axis BaØi 5: Trích trang 11-147 sách ‘Control System Toolbox’ Vẽ đáp ứng Nichols của hệ thống có hàm truyền: −4s4 + 48 s3 − 18 s2 + 250 s + 600 H() s = s4 +30 s3 + 282 s2 +525 s + 60 » H=tf([-4 48 -18 250 600],[1 30 282 525 60]) Transfer function: -4 s^4 + 48 s^3 - 18 s^2 + 250 s + 600 s^4 + 30 s^3 + 282 s^2 + 525 s + 60 Nichols(H) ngrid
  31. MATLAB trong điều khiển tự động Trang 188 Nichols Charts From: U(1) 20 15 10 5 To: Y (1) 0 Open-Loop Gain (dB) Open-Loop -5 -10 -15 -600 -500 -400 -300 -200 -100 0 Open-Loop Phase (deg) Bài 6: Trang 131 sách ‘Ứng dụng MATLAB trong điều khiển tự động’ tác giả Nguyễn Văn Giáp. Trên giản đồ Nichols vẽ đường cong logarit biên độ – pha của hàm truyền hệ thống k H(s) = S3+52s2+100s » k=438; » num=k; » den=[1 52 100 0]; » w=.1:.1:10; » [mag,phase]=bode(num,den,w); » ngrid, Kết quả:
  32. MATLAB trong điều khiển tự động Trang 189 40 0 dB 30 0.25 dB 0.5 dB 20 1 dB -1 dB 10 3 dB -3 dB 6 dB 0 -6 dB -10 -12 dB Open-Loop Gain (dB) Open-Loop -20 -20 dB -30 -40 -40 dB -350 -300 -250 -200 -150 -100 -50 0 Open-Loop Phase (deg)
  33. MATLAB trong điều khiển tự động Trang 190 NHÓM LỆNH VỀ ĐÁP ỨNG THỜI GIAN (Time Response) 1. Lệnh IMPULSE a) Công dụng: Tìm đáp ứng xung đơn vị. b) Cú pháp: [y,x,t] = impulse(a,b,c,d) [y,x,t] = impulse(a,b,c,d,iu) [y,x,t] = impulse(a,b,c,d,iu,t) [y,x,t] = impulse(num,den) [y,x,t] = impulse(num,den,t) c) Giải thích: Lệnh impulse tìm đáp ứng xung đơn vị của hệ tuyến tính. Nếu bỏ qua các đối số bên trái thì lệnh impulse sẽ vẽ ra đáp ứng xung trên màn hình. impulse(a,b,c,d) tạo ra chuỗi đồ thị đáp ứng xung, mỗi đồ thị ứng với một mối quan hệ vào ra của hệ liên tục LTI: . x = Ax + Bu y = Cx + Du với vector thời gian được xác định tự động. impulse(a,b,c,d,iu) tạo ra đáp ứng xung từ ngõ vào duy nhất iu tới toàn bộ các ngõ ra của hệ thống với vector thời gian được xác định tự động. iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được dùng cho đáp ứng xung. impulse(num,den) tạo ra đồ thị đáp ứng xung của đa thức hàm truyền: G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. impulse(a,b,c,d,iu,t) hay impulse(num,den,t) dùng vector thời gian t do người sử dụng quy định. Vector t chỉ định những thời điểm mà đáp ứng xung được tính và vector t phải được chỉ chia thành các khoảng đều nhau. Nếu giữ các đối số bên trái: [y,x,t] = impulse(a,b,c,d) [y,x,t] = impulse(a,b,c,d,iu) [y,x,t] = impulse(a,b,c,d,iu,t) [y,x,t] = impulse(num,den) [y,x,t] = impulse(num,den,t) không vẽ ra các đồ thị mà tạo ra các ma trận đáp ứng trạng thái và đáp ứng ngõ ra của hệ thống và vector thời gian t. Ma trận y và x chứa các đáp ứng trạng thái và đáp ứng ngõ ra của hệ thống được xác định tại những thời điểm t. Ma trận y có số cột là số ngõ ra và mỗi hàng ứng với một thành phần trong vector t. Ma trận x có số cột là số trạng thái và mỗi hàng ứng với một thành phần trong vector t. d) Ví dụ: (Trích từ trang 11-95 sách ‘control System Toolbox’)
  34. MATLAB trong điều khiển tự động Trang 191 Vẽ đáp ứng xung của hệ không gian trạng thái bậc 2 sau: ⎡ . ⎤ x1 ⎡−0.5 − 0.8⎤ ⎡x1 ⎤ ⎡1⎤ ⎢ . ⎥ = ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ u ⎢ ⎥ 0.8 0 x 0 ⎣x2 ⎦ ⎣ ⎦ ⎣ 2 ⎦ ⎣ ⎦ ⎡x1 ⎤ y= [][ 1.9 6.5 ⎢ ⎥ + 0] u ⎣x 2 ⎦ % Khai báo hệ thống: a = [-0.5 -0.8 ; 0.8 0]; b = [1 ; 0]; c = [1.9 6.5]; d = [0]; % Vẽ đáp ứng xung: impulse(a,b,c,d); title(‘Dap ung xung’) (đặt tiêu đề cho đồ thị) và cuối cùng ta nhận được đồ thị đáp ứng xung như sau: 2. Lệnh DIMPULSE a) Công dụng: Tìm đáp ứng xung đơn vị của hệ gián đoạn. b) Cú pháp: [y,x] = dimpulse(a,b,c,d) [y,x] = dimpulse(a,b,c,d,iu) [y,x] = dimpulse(a,b,c,d,iu,n)
  35. MATLAB trong điều khiển tự động Trang 192 [y,x] = dimpulse(num,den) [y,x] = dimpulse(num,den,n) c) Giải thích: Lệnh dimpulse tìm đáp ứng xung đơn vị của hệ tuyến tính gián đoạn. Nếu bỏ qua các đối số bên trái thì thì lệnh dimpulse sẽ vẽ ra đáp ứng xung trên màn hình. dimpulse(a,b,c,d) tạo ra chuỗi đồ thị đáp ứng xung, mỗi đồ thị đáp ứng với một mối quan hệ vào ra của hệ gián đoạn LTI: x[n + 1] = Ax[n] + Bu[n] y[n] = Cx[n] + Du[n] với số điểm lấy mẫu được xác định tự động. dimpulse(a,b,c,d,iu) tạo ra đồ thị đáp ứng xung từ ngõ vào duy nhất iu tới toàn bộ các ngõ ra của hệ thống với số điểm lấy mẫu được xác định tự động. iu là chỉ số ngõ vào của hệ thống và chỉ ra ngõ vào nào được dùng cho đáp ứng xung. dimpulse(num,den) tạo ra đồ thị đáp ứng xung của đa thức hàm truyền: G(z) = num(z)/den(z) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của z. dimpulse(num,den,n) hay dimpulse(a,b,c,d,iu,n) dùng số điểm lấy mẫu n do người sử dụng chỉ định. Nếu giữ các đối số bên trái: [y,x] = dimpulse(a,b,c,d) [y,x] = dimpulse(a,b,c,d,iu) [y,x] = dimpulse(a,b,c,d,iu,n) [y,x] = dimpulse(num,den) [y,x] = dimpulse(num,den,n) không vẽ ra các đồ thị mà tạo ra các ma trận đáp ứng ngõ ra và đáp ứng trạng thái của hệ thống. Ma trận y và x chứa các đáp ứng trạng thái và ngõ ra của hệ thống được xác định tại những điểm lấy mẫu. Ma trận y có số cột là số ngõ ra. Ma trận x có số cột là số trạng thái. d) Ví dụ: Vẽ đáp ứng xung của hệ gián đoạn có hàm truyền sau: 2z2 + 3.4 z + 1.5 H() z = z 2 −1.6 + 0.8 num = [2 -3.4 1.5]; den = [1 -1.6 0.8]; dimpulse(num,den); title(‘Dap ung xung he gian doan’) và cuối cùng ta được đồ thị đáp ứng xung hệ gián đoạn như sau:
  36. MATLAB trong điều khiển tự động Trang 193 3. Lệnh INITIAL a) Công dụng: Tìm đáp ứng điều kiện ban đầu. b) Cú pháp: [y,x,t] = initial(a,b,c,d,x0) [y,x,t] = initial(a,b,c,d,x0,t) c) Giải thích: Lệnh initial dùng để tìm đáp ứng của hệ tuyến tính liên tục ứng với điều kiện ban đầu của các trạng thái. Nếu bỏ qua các đối số ở bên trái thì lệnh initial sẽ vẽ ra đáp ứng điều kiện ban đầu trên màn hình. initial(a,b,c,d,x0) vẽ ra đồ thị đáp ứng điều kiện ban đầu của tất cả các ngõ ra của hệ liên tục LTI: . x = Ax + Bu y = Cx + Du với vector thời gian được xác định tự động. x0 là vector trạng thái ban đầu. initial(a,b,c,d,x0,t) vẽ ra đồ thị đáp ứng ban đầu với vector thời gian t do người sử dụng xác định. Vector t sẽ chỉ ra những thời điểm mà tại đó đáp ứng điều kiện ban đầu được tính. Nếu sử dụng các đối số ở vế trái của dòng lệnh thì: [y,x,t] = initial(a,b,c,d,x0) [y,x,t] = initial(a,b,c,d,x0,t) sẽ không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận đáp ứng trạng thái x, đáp ứng ngõ ra y và vector thời gian t của hệ thống đối với điều kiện ban đầu x0. Ma trận y và x chứa các
  37. MATLAB trong điều khiển tự động Trang 194 đáp ứng ngõ ra và đáp ứng trạng thái của hệ thống được tính tại thời điểm t. Ma trận y có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector t. Ma trận x có số cột bằng số trạng thái và mỗi cột ứng với một thành phần trong vector t. d) Ví dụ: Vẽ đáp ứng điều kiện ban đầu của hệ không gian trạng thái bậc 2 sau: . ⎡ ⎤ ⎡−0.5572 − 0.7814⎤ ⎡x ⎤ ⎡1⎤ ⎢x1 ⎥ = 1 + u . ⎢ 0.7814 0 ⎥ ⎢x ⎥ ⎢0⎥ ⎣⎢x2 ⎦⎥ ⎣ ⎦ ⎣ 2 ⎦ ⎣ ⎦ ⎡x1 ⎤ y = []1.9691 6.4493 ⎢ ⎥ + [0]u ⎣x2 ⎦ với điều kiện ban đầu x0 = [1 0] % Khai báo hệ thống, điều kiện ban đầu và trục thời gian: a = [-0.5572 -0.7814 ; 0.7814 0]; b = [1 ; 0]; c = [1.9691 6.4493]; d = [0]; x0 = [1 0]; t = 0:0.1:20; % Vẽ đáp ứng: initial(a,b,c,d,x0,t) title(‘Dap ung dieu kien ban dau’)
  38. MATLAB trong điều khiển tự động Trang 195 4. Lệnh DINITIAL a) Công dụng: Tìm đáp ứng điều kiện ban đầu của hệ gián đoạn. b) Cú pháp: [y,x] = dinitial(a,b,c,d,x0) [y,x] = dinitial(a,b,c,d,x0,n) c) Giải thích: Lệnh dinitial dùng để tìm đáp ứng của hệ tuyến tính gián đoạn ứng với điều kiện ban đầu của các trạng thái. Nếu bỏ qua các đối số ở bên trái thì lệnh dinitial sẽ vẽ ra đáp ứng điều kiện ban đầu trên màn hình. dinitial(a,b,c,d,x0) vẽ ra đồ thị đáp ứng điều kiện ban đầu của tất cả các ngõ ra của hệ gián đoạn LTI: x[n + 1] = Ax[n] + Bu[n] y[n] = Cx[n] + Du[n] với số điểm lấy mẫu được xác định tự động. x0 là vector trạng thái ban đầu. dinitial(a,b,c,d,x0,n) vẽ ra đồ thị đáp ứng ban đầu với số điểm lầy mẫu n do người sử dụng xác định. Nếu sử dụng các đối số ở vế trái của dòng lệnh thì: [y,x] = dinitial(a,b,c,d,x0) [y,x] = dinitial(a,b,c,d,x0,n) sẽ không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận đáp ứng trạng thái x, đáp ứng ngõ ra y của hệ thống đối với điều kiện ban đầu x0. Ma trận y có số cột bằng số ngõ ra và ma trận x có số cột bằng số trạng thái. d) Ví dụ: Vẽ đáp ứng điều kiện ban đầu của hệ không gian trạng thái bậc 2: x [n+ 1] ⎡ 1 ⎤ ⎡−0.7497 − 0.2027⎤ ⎡x1 [n]⎤ ⎡−4.1841⎤ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ u ⎣⎢x2 [n+ 1]⎦⎥ ⎣ 1 0 ⎦ ⎣x2 [n]⎦ ⎣−6.5049⎦ ⎡x1[] n ⎤ y = []3.9321 0 ⎢ ⎥ ⎣x2 [] n ⎦ với điều kiện ban đầu x0 = [1 0] a = [-0.7497 -0.2027 ; 1 0]; b = [-4.1841 ; -6.5049]; c = [3.9321 0]; d = [0]; dinitial(a,b,c,d,[1 0]); title(‘Dap ung dieu kien ban dau cua he gian doan’)
  39. MATLAB trong điều khiển tự động Trang 196 5. Lệnh LSIM a) Công dụng: Mô phỏng hệ thống liên tục với các ngõ vào tùy ý. b) Cú pháp: [y,c] = lsim(a,b,c,d,u,t) [y,c] = lsim(a,b,c,d,u,t,x0) [y,c] = lsim(num,den,u,t) c) Giải thích: Lệnh lsim dùng để mô phỏng hệ tuyến tính liên tục với các ngõ vào tùy ý. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh lsim vẽ ra ra đồ thị trên màn hình. Cho hệ không gian trạng thái LTI: . x = Ax + Bu y = Cx + Du lsim(a,b,c,d,u,t) vẽ ra đồ thị đáp ứng thời gian của hệ thống với ngõ vào thời gian ban đầu nằm trong ma trận u. Ma trận u phải có số cột bằng số ngõ vào u. Mỗi hàng của ma trận u tương ứng với một thời gian mới và ma trận u phải có số hàng là length(t). Vector t chỉ ra trục thời gian cho quá trình mô phỏng và phải chia thành các đoạn bằng nhau. Nếu dùng thêm đối số x0 ở vế phải thì lệnh lsim(a,b,c,d,u,t,x0) sẽ chỉ ra điều kiện ban đầu của các trạng thái. lsim(num,den,u,t) vẽ ra đáp ứng thời gian của hàm truyền đa thức: G(s) = num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s.
  40. MATLAB trong điều khiển tự động Trang 197 Nếu giữ lại các đối số ở vế trái thì: [y,c] = lsim(a,b,c,d,u,t) [y,c] = lsim(a,b,c,d,u,t,x0) [y,c] = lsim(num,den,u,t) sẽ không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận y và x, trong đó ma trận y là đáp ứng ngõ ra và ma trận x là đáp ứng trạng thái của hệ thống. Ma trận y có số cột bằng số ngõ ra y và mỗi hàng ứng với một hàng của ma trận u. Ma trận x có số cột bằng số trạng thái x và mỗi hàng ứng với một hàng của ma trận u. d) Ví dụ: (Trích từ trang 11-127 sách ‘Control System Toolbox’) Mô phỏng và vẽ đồ thị đáp ứng của hệ thống có hàm truyền: 2s2 + 5 s + 1 H() s = s2 +2 s + 3 với chu kỳ sóng vuông là 10s. num = [2 5 1]; den = [1 2 3]; t = 0:.1:10; period = 4; u = (rem(t,period)) >= period./2); lsim(num,den,u,t); title(‘Dap ung cua song vuong’) và ta được đồ thị đáp ứng của hệ như sau: 6. Lệnh DLSIM a) Công dụng:
  41. MATLAB trong điều khiển tự động Trang 198 Mô phỏng hệ thống gián đoạn với các ngõ vào tùy ý. b) Cú pháp: [y,c] = dlsim(a,b,c,d,u,t) [y,c] = dlsim(a,b,c,d,u,x0) [y,c] = dlsim(num,den,u) c) Giải thích: Lệnh lsim dùng để mô phỏng hệ tuyến tính gián đoạn với các ngõ vào tùy ý. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh dlsim vẽ ra ra đồ thị trên màn hình. Cho hệ không gian trạng thái LTI: x[n + 1] = Ax[n] + Bu[n] y[n] = Cx[n] + Du[n] dlsim(a,b,c,d,u) vẽ ra đồ thị đáp ứng thời gian của hệ thống với ngõ vào thời gian ban đầu nằm trong ma trận u. Ma trận u phải có số cột bằng số ngõ vào u. Mỗi hàng của ma trận u tương ứng với một thời điểm mới. Nếu dùng thêm đối số x0 ở vế phải thì lệnh lsim(a,b,c,d,u,x0) sẽ chỉ ra điều kiện ban đầu của các trạng thái. lsim(num,den,u) vẽ ra đáp ứng thời gian của hàm truyền đa thức: G(z) = num(z)/den(z) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. Nếu giữ lại các đối số ở vế trái thì: [y,c] = dlsim(a,b,c,d,u) [y,c] = dlsim(a,b,c,d,u,x0) [y,c] = dlsim(num,den,u) sẽ không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận y và x, trong đó ma trận y là đáp ứng ngõ ra và ma trận x là đáp ứng trạng thái của hệ thống. Ma trận y có số cột bằng số ngõ ra y và mỗi hàng ứng với một hàng của ma trận u. Ma trận x có số cột bằng số trạng thái x và mỗi hàng ứng với một hàng của ma trận u. d) Ví dụ: Mô phỏng đáp ứng của hệ thống gián đoạn có hàm truyền: 2z2 − 3.4 z + 1.5 H() z = z2 −1.6 z + 0.8 với 100 mẫu của nhiễu ngẫu nhiên. num = [2 -3.4 1.5]; den = [1 -1.6 0.8]; rand(‘nomal’) u = rand(100,1); dlsim(num,den,u) title(‘Dap ung nhieu’) và ta được đồ thị đáp ứng của hệ như sau:
  42. MATLAB trong điều khiển tự động Trang 199 7. Lệnh STEP a) Công dụng: Tìm đáp ứng nấc đơn vị. b) Cú pháp: [y,x,t] = step(a,b,c,d) [y,x,t] = step(a,b,c,d,iu) [y,x,t] = step(a,b,c,d,iu,t) [y,x,t] = step(num,den) [y,x,t] = step(num,den,t) c) Giải thích: Lệnh step tìm đáp ứng nấc đơn vị của hệ tuyến tính liên tục. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh step vẽ ra đáp ứng nấc trên màn hình. step(a,b,c,d) vẽ ra chuỗi đồ thị đáp ứng nấc, mỗi đồ thị tương ứng với mối quan hệ giữa một ngõ vào và một ngõ ra của hệ liên tục LTI: . x = Ax + Bu y = Cx + Du với vector thời gian được xác định tự động. step(a,b,c,d,iu) vẽ ra đồ thị đáp ứng nấc từ một ngõ vầo duy nhất tới tất cả các ngõ ra của hệ thống với vector thời gian được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và nó chỉ ra ngõ vào nào được sử dụng cho đáp ứng xung. step(num,den) vẽ ra đồ thị đáp ứng nấc của hàm truyền đa thức:
  43. MATLAB trong điều khiển tự động Trang 200 G(s) =num(s)/den(s) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. step(a,b,c,d,iu,t) hay step(num,den,t) cũng vẽ ra đáp ứng nấc của hệ không gian trạng thái hay hàm truyền với vector thời gian t do người sử dụng xác định. Vector t chỉ ra những thời điểm mà tại đó đáp ứng nấc được tính và vector t phải được chia thành những đoạn đều nhau. Nếu giữ lại các đối số ở vế trái của dòng lệnh thì: [y,x,t] = step(a,b,c,d) [y,x,t] = step(a,b,c,d,iu) [y,x,t] = step(a,b,c,d,iu,t) [y,x,t] = step(num,den) [y,x,t] = step(num,den,t) không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận đáp ứng ngõ ra y và ma trận đáp ứng trạng thái x củahệ thống được xác định tại những thời điểm t. Ma trận y có số cột bằng số ngõ ra và mỗi hàng ứng với một thành phần trong vector t. Ma trận x có số cột bằng số trạng thái và mỗi hàng ứng với một thành phần trong vector t. d) Ví dụ: Vẽ đồ thị đáp ứng nấc của hệ không gian trạng thái bậc 2 sau: . ⎡ ⎤ ⎡−0.5572 − 0.7814⎤ ⎡x ⎤ ⎡1⎤ ⎢x1 ⎥ = 1 + u . ⎢ 0.7814 0 ⎥ ⎢x ⎥ ⎢0⎥ ⎣⎢x2 ⎦⎥ ⎣ ⎦ ⎣ 2 ⎦ ⎣ ⎦ ⎡x1 ⎤ y = []1.9691 6.4493 ⎢ ⎥ + [0]u ⎣x2 ⎦ a = [-0.5572 -0.7814 ; 0.7814 0]; b = [1 ; 0]; c = [1.9691 6.4493]; d = [0]; step(a,b,c,d); title(‘Dap ung nac’) và ta được đồ thị đáp ứng nấc của hệ thống như sau:
  44. MATLAB trong điều khiển tự động Trang 201 8. Lệnh DSTEP a) Công dụng: Tìm đáp ứng nấc đơn vị của hệ gián đoạn. b) Cú pháp: [y,x] = dstep(a,b,c,d) [y,x] = dstep(a,b,c,d,iu) [y,x] = dstep(a,b,c,d,iu,n) [y,x] = dstep(num,den) [y,x] = dstep(num,den,n) c) Giải thích: Lệnh dstep tìm đáp ứng nấc đơn vị của hệ tuyến tính gián đoạn. Nếu bỏ qua các đối số ở vế trái của dòng lệnh thì lệnh dstep vẽ ra đáp ứng nấc trên màn hình. dstep(a,b,c,d) vẽ ra chuỗi đồ thị đáp ứng nấc, mỗi đồ thị tương ứng với mối quan hệ giữa một ngõ vào và một ngõ ra của hệ gián đoạn LTI: x[n + 1] = Ax[n] + Bu[n] y[n] = Cx[n] + Du[n] với số điểm lấy mẫu được xác định tự động. dstep(a,b,c,d,iu) vẽ ra đồ thị đáp ứng nấc từ một ngõ vầo duy nhất tới tất cả các ngõ ra của hệ thống với số điểm lấy mẫu được xác định tự động. Đại lượng vô hướng iu là chỉ số ngõ vào của hệ thống và nó chỉ ra ngõ vào nào được sử dụng cho đáp ứng xung. dstep(num,den) vẽ ra đồ thị đáp ứng nấc của hàm truyền đa thức:
  45. MATLAB trong điều khiển tự động Trang 202 G(z) =num(z)/den(z) trong đó num và den chứa các hệ số đa thức theo chiều giảm dần số mũ của s. dstep(a,b,c,d,iu,n) hay dstep(num,den,n) cũng vẽ ra đáp ứng nấc của hệ không gian trạng thái hay hàm truyền với số điểm lấy mẫu do người sử dụng xác định. Nếu giữ lại các đối số ở vế trái của dòng lệnh thì: [y,x] = dstep(a,b,c,d) [y,x] = dstep(a,b,c,d,iu) [y,x] = dstep(a,b,c,d,iu) [y,x] = dstep(num,den) [y,x] = dstep(num,den,n) không vẽ ra các đồ thị đáp ứng mà tạo ra các ma trận đáp ứng ngõ ra y và ma trận đáp ứng trạng thái x củahệ thống. Ma trận y có số cột bằng số ngõ ra. Ma trận x có số cột bằng số trạng thái. d) Ví dụ: Vẽ đáp ứng nấc của hệ gián đoạn của hệ có hàm truyền như sau: 2z2 − 3.4 z + 1.5 H() z = z2 −1.6 z + 0.8 num = [2 -3.4 1.5]; den = [1 -1.6 0.8]; dstep(num,den) title(‘Dap ung nac he gian doan’) và ta được đồ thị đáp ứng nấc của hệ như hình bên:
  46. MATLAB trong điều khiển tự động Trang 203 9. Lệnh LTITR a) Công dụng: Tìm đáp ứng thời gian của hệ tuyến tính bất biến. b) Cú pháp: ltitr(a,b,u) ltitr(a,b,u,x0) c) Giải thích: Lệnh ltitr dùng để mở rộng đáp ứng thời gian của hệ tuyến tính bất biến. Nó mô phỏng cho hệ không gian trạng thái gián đoạn: x = ltitr(a,b,u) mở rộng đáp ứng của hệ gián đoạn: x[n + 1] = Ax[n] + Bu[n] đối với ngõ vào u. Ma trận u phải có số cột bằng số ngõ vào u. Mỗi hàng của ma trận u tương ứng với một điểm thời gian mới. ltitr tạo ra ma trận x với số cột bằng số trạng thái x và có số hàng là length(u). Nếu thêm vào vế phải dòng lệnh tham số x0 thì điều kiện ban đầu sẽ được thiết lập với lệnh x = ltitr(a,b,u,x0) 10. Lệnh FILTER a) Công dụng: Lọc dữ liệu với đáp ứng xung không xác định hay đáp ứng xung xác định. b) Cú pháp: y = filter(b,a,X) [y,zf] = filter(b,a,X) [y,zf] = filter(b,a,X,zi) y = filter(b,a,X,zi,dim) [ ] = filter(b,a,X,[ ],dim) c) Giải thích: Lệnh fiter lọc dữ liệu tuần tự sử dụng bộ lọc số cho các ngõ vào thực và phức. y = filter(b,a,X) lọc dữ liệu trong vector X với bộ lọc được mô tả bởi vector hệ số tử số b và vector hệ số mẫu số a. Nếu a(1) không bằng 1, bộ lọc sẽ chuẩn hóa hệ số lọc bởi a(1). Nếu a(1) bằng 0 thì sẽ báo lỗi. Nếu X là một ma trận, bộ lọc sẽ thực hiện trên các cột của X. Nếu X là một mảng đa chiều, bộ lọc sẽ thực hiện theo chiều duy nhất. [y,zf] = filter(b,a,X) tạo ma trận điều kiện cuối cùng zf của bộ trễ. Ngõ ra zf là một vector của max(size(a),size(b)) hoặc một tập hợp các vector với mỗi vector là một cột của X. [y,zf] = filter(b,a,X,zi) chấp nhận điều kiện ban đầu zi và tạo ra điều kiện cuối cùng cuối cùng zf của bộ lọc trễ. Ngõ vào zi là một vector có kích thước length(a),length(b)) – 1. y = filter(b,a,X,zi,dim) và [ ] = filter(b,a,X,[ ],dim) thực hiện lọc theo chiều dim.
  47. MATLAB trong điều khiển tự động Trang 204 CÁC BÀI TẬP VỀ ĐÁP ỨNG THỜI GIAN Bài1: Lệnh pade: Tính toán sắp xỉ Bài này trích từ trang 11-66 sách ‘Control System Toollbox’ » pade(0.1,3) Step response of 3rd-order Pade approximation 1.5 1 0.5 0 Amplitude -0.5 -1 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Time (secs) Phase response 0 -200 -400 -600 Phase (deg.) -800 -1000 1 2 3 10 10 10 Frequency (rad/s) Bài 2: Trích từ trang 11-24 sách ‘Control System Toollbox’ s -1 H(s) = s2 + 4s +5 » H=tf([1 -1],[1 4 5],'inputdelay',035) Transfer function: s - 1 exp(-35*s) * s^2 + 4 s + 5 » Hd=c2d(H,0.1,'foh') Transfer function:
  48. MATLAB trong điều khiển tự động Trang 205 0.04226 z^2 - 0.01093 z - 0.03954 z^(-350) * z^2 - 1.629 z + 0.6703 Sampling time: 0.1 » step(H,'-',Hd,' ') Step Response From: U(1) 0.15 0.1 0.05 0 -0.05 To: Y (1) Amplitude -0.1 -0.15 -0.2 -0.25 0 5 10 15 20 25 30 35 40 Time (sec.) 2s2 + 5s + 1 s2 + 2s + 3 Bài 3: Trang 11-127, H(s) = s - 1 s2+s+5 » [u,t]=gensig('square',4,10,0.1); » H=[tf([2 5 1],[1 2 3]);tf([1 -1],[1 1 5])]; » lsim(H,u,t) Kết quả: Bài tập này được trích từ trang 11-127 sách ‘Control System Toolbox’
  49. MATLAB trong điều khiển tự động Trang 206 Linear Simulation Results 3 2 1 0 To: Y (1) -1 -2 0.4 Amplitude 0.2 0 -0.2 To: Y (2) -0.4 -0.6 0 1 2 3 4 5 6 7 8 9 10 Time (sec.) Bài 4: Dùng lệnh lsim, trích từ trang 11-130 sách ‘Control Systen Toollbox’ Dịch đề: Vẽ đáp ứng khâu bậc 2 của hàm truyền sau: ω2 h() s = s2 + 2s + ω2 ω = 62,83 » w2=62.83^2 w2 = 3.9476e+003 » h=tf(w2,[1 2 w2]); » t=0:0.1:5; %vector of time sample: » u=(rem(t,1)>=0.5); %square ware value : » lsim(h,u,t) Kết quả:
  50. MATLAB trong điều khiển tự động Trang 207 Linear Simulation Results 2 1.5 1 0.5 To: Y (1) Amplitude 0 -0.5 -1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (sec.) Bài 5: Trang 11-131 sách ‘Control Systen Toollbox’ Ta lấy số liệu bài 24 nhưng thời gian mẫu là 0,1. Chương trình: » w2=62.83^2; » hd=c2d(h,0.1); » t=0:0.1:5; %vector of time sample: » u=(rem(t,1)>=0.5); %square ware value : » lsim(hd,u,t)
  51. MATLAB trong điều khiển tự động Trang 208 Linear Simulation Results 0.7 0.6 0.5 0.4 To: Y (1) 0.3 Amplitude 0.2 0.1 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (sec.) Bài 6: Trang 11-132 sách ‘Control Systen Toollbox’ Cũng lấy số liệu 2 bài trên. » w2=62.83^2; » h=tf(w2,[1 2 w2]); » t=0:0.1:5; %vector of time sample: » u=(rem(t,1)>=0.5); %square ware value : » hd=c2d(h,0.1); » lsim(h,'b ',hd,'r-',u,t) %
  52. MATLAB trong điều khiển tự động Trang 209 Linear Simulation Results 2 1.5 1 0.5 To: Y (1) Amplitude 0 -0.5 -1 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time (sec.) Bài 7: Trích từ trang 46 sách ‘ứng dụng matlab trong điều khiển tự động’ Phương trình biến trang thái của hệ thống tuyến tính bất biến theo thời gian là: Chương trình được viết trong file.m: %function [yout,x] = lsim(A, B, C, D, U, t, x0) %Phuong trinh bien trang thai cua mot he thong tuyen tinh % bat bien theo thoi gian la: % . % x1 % . 0 1 0 x1 1 % {x2} = { 0 0 1 } { x2 } + {1} r(t) % . -6 -11 -6 x3 1 % x3 % 1 % y=[1 1 0]x, x(0)= 0.5 % -0.5 % Xac dinh x(t),y(t) khi r(t) la ham bac don vi hold on grid on A=[0 1 0;0 0 1;-6 -13 -6]; B=[1;1;1];%xac dinh vi ban dau va hinh dang cua do thi x1,y,x2,x3 C=[1 1 0]; D=0; x0=[1 .5 -.5]; %vecto hang dieu kien ban dau
  53. MATLAB trong điều khiển tự động Trang 210 t=0:.05:8; %buoc nhay U=ones(1,length(t));%tao vecto hang u(t) [x,y]=lsim(A,B,C,D,U,t,x0); plot(t,x,t,y) title('BAI GIAI BT15') xlabel('Thoi gian-giay') text(3.8,1.8,'y'),text(3.8,2.6,'x1');%Canh vi tri cua y va x1 tren do thi text(3.8,-0.6,'x2'),text(3.8,-1.4,'x3')%Canh vi tri cua x2 va x3 tren do thi BAI GIAI BT15 4 Bài 9: 3 trích x1 từ trang 2 48 y sách tác 1 giả Nguy ễn 0 Văn x2 Giáp. -1 x3 Cũng với -2 yêu 0 1 2 3 4 5 6 7 8 cầu Thoi gian-giay như bài 28, nhưng r(t)=sin(2Πt). Chương trình soạn trong file.m: %function [yout,x] = lsim(A, B, C, D, U, t, x0) %BT16:Ve do thi y(t),x(t) cua bai BT15 neu r(t)=sin(2pit) A=[0 1 0;0 0 1;-6 -11 -6]; B=[1;1;1];C=[1 1 0];D=0; x0=[1 .5 -.5]; %vecto hang dieu kien ban dau t=0:.05:4; %buoc nhay r=sin(2*pi*t); [y,x]=lsim(A,B,C,D,r,t,x0); plot(t,x,t,y) title('BAI GIAI BT16') xlabel('Thoi gian-giay') text(3.8, 1.8,'y'),text(3.8, 2.6,'x1')
  54. MATLAB trong điều khiển tự động Trang 211 text(3.8, -8,'x2'),text(3.8, -1.4,'x3') BAI GIAI BT16 2 y 1.5 1 0.5 0 -0.5 -1 x3 -1.5 -2 0 0.5 1 1.5 2 2.5 3 3.5 4 Thoi gian-giay ài Bài 10: Xét hàm truyền sau: s+ 10 G(s)= s2 + 8s + 25 Để tính đáp ứng bước của hệ thống này ta dùng cấu trúc như sau: [out,state,tt]=step([1 10],[1 8 25]) Giả sử ta muốn phân tích một đáp ứng bước của hệ thống thay đổi, với zero của hàm truyền thay đổi nhưng độ lợi dc (dc gain) của hệ thống không đổi, để giữ lại cho hệ thống cùng mẫu và thay đổi hệ số của số hạng đầu trong đa thức của tử,tức là hệ số của s, vì vậy mà dc gain là hằng số và zero thay đổi. Ví dụ : hệ thống như ví dụ trên nhưng số hạng ban đầu của đa thức ở tử số thay đổi thành (-4,- 2,-1,0,1,2,4) Ta thưc hiện trong cửa sổ lệnh của matlab như sau: » coef=[-4 -2 -1 0 1 2 4]; » den=[1 8 25]; » [y,x,t]=step([coef' 10*ones(length(coef),1)],den); » mesh(coef,t,y) Kết quả như hình:
  55. MATLAB trong điều khiển tự động Trang 212 Hình 3.7: So sánh giữa các đáp ứng step Bài 11: đáp ứng xung (impulse) Ví dụ hệ thống có hàm truyền sau: s+ 10 G(s)= s2 + 2s + 25 Vẽ đáp ứng xung của hệ thống: impulse([1 10],[1 2 25]) Giả sử ta muốn phân tích đáp ứng xung thay đổi như thế nào khi zero của hàm truyền thay đổi, không thay đổi dc gain của hệ thống. giống như ví dụ ở phần trước ta có : » coef=[-4 -2 -1 0 1 2 4]; » den=[1 2 25]; » impulse([coef' 10*ones(length(coef),1)],den);
  56. MATLAB trong điều khiển tự động Trang 213 Kết quả như hình sau: Bài 12: Trích từ trang 716 sách ‘The Student Edition of MATLAB’ Dịch đề: Thiết kế 1 khâu gồm 10 bộ lọc của dãi băng truyền ngang có tần số từ 100 đến 200 Hz và vẽ đáp ứng xung của nó: » n=5;wn=[100 200]/500; » [b,a]=butter(n,wn); » [y,t]=impz(b,a,101); » stem(t,y)
  57. MATLAB trong điều khiển tự động Trang 214 0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 0 10 20 30 40 50 60 70 80 90 100 Bài 13: Đáp ứng từng ngõ vào Một vấn đề tổng quát hơn là ta có thể tính được tín hiệu ngõ ra của hệ thống LTI với một tính hiệu ngõ vào không đồng nhất. Ví dụ như hệ thống bậc nhất sau: . x= − x + u y= x Hệ thống này bị tác động với một tín hiệu ngõ vào hình sin có tần số là 1Hz, tín hiệu ngõ ra thu được bởi cấu trúc: >> freq=1; t=0:0.05:10; >> u=sin(2*pi*freq*t); lsim(-1,1,1,0,u,t)
  58. MATLAB trong điều khiển tự động Trang 215 Kết quả là hình sau: Hình : Đáp ứng từng ngõ vào