Giáo trình Matlab - Chương 6: Matlab và điều khiển tự động (Phần 2)

[],c,[]); 
initial(sys,x0) 
134
b. Đáp ứng xung Dirac: Ta tìm đáp ứng của hệ thống với xung nhờ hàm 
impulse. Các lệnh MATLAB (lưu trong ct6_24.m)như sau: 
clc 
a = [‐0.5572 ‐0.7814;0.7814  0]; 
b = [1 ‐1;0 2]; 
c = [1.9691  6.4493]; 
sys = ss(a,b,c,0); 
impulse(sys) 
Hình bên trái là đáp ứng của kênh thứ nhất và hình bên phải là đáp ứng của 
kênh thứ 2.  
c. Đáp ứng đối với hàm bước nhảy: Để tìm đáp ứng của hệ thống đối với 
hàm bước nhảy  ta dùng hàm  step. Các  lệnh MATLAB  (lưu  trong  ct6_25.m) 
như sau:  
clc 
a = [‐0.5572   ‐0.7814;0.7814  0]; 
b = [1 ‐1;0 2]; 
c = [1.9691  6.4493]; 
sys = ss(a,b,c,0); 
step(sys) 
Hình bên trái là đáp ứng của kênh thứ nhất và hình bên phải là đáp ứng của 
kênh thứ 2. 
d. Đáp ứng với tín hiệu bất kỳ: Để tìm đáp ứng của hệ thống đối với hàm 
bất kì ta dùng hàm lsim. Các lệnh MATLAB (lưu trong ct6_26.m) như sau:  
clc 
[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) 
Ta dùng hàm gensig để tạo một xung hình vuông, trong 4 chu kỳ và lấy mẫu 
sau 0.1s trong 10 chu kỳ. 
3. Đáp ứng trong miền tần số của hệ thống: Cho một hàm truyền của một hệ 
thống,thay  s bằng  jω  ta  có hàm  truyền  đạt  tần  số của hệ  thống  đó. Độ  rộng 
băng của hệ  thống ωB  là  tần số mà  tại đó biên độ của g giảm đi 1/√2. Tần số 
ứng với giá trị max của G(ω) gọi là ωr và có trị số là: 
2
nr 21 ζ−ω=ω  
Để vẽ đặc tính tần biên‐pha của một hệ thống ta dùng lệnh freqs. 
135
Ví dụ: Cho hàm truyền của một hệ thống là: 
4s2s
4)s(G 2 ++=  
Tìm  đặc  tính  tần  biên‐pha  của  hệ  thống  bằng  các  lệnh MATLAB(lưu  trong 
ct6_27.m): 
w = 0:0.01:3; 
ms = [1 2 4]; 
ts = [4]; 
freqs(ts,ms,w); 
Ta cũng có thể tạo đồ thị như sau(lưu trong ct6_28.m): 
ts = [4]; 
ms = [1 2 4]; 
w = 0:0.01:3; 
g = freqs(ts,ms,w); 
mag = abs(g); 
pha = angle(g); 
subplot(2,1,1); 
loglog(w,mag); 
grid on; 
subplot(2,1,2); 
semilogx(w,pha); 
grid on 
Ngược lại khi có đặc tính tần biên ‐ pha ta có thể tìm lại được hàm truyền bằng 
lệnh invfreqs.  
Ví dụ: Tìm hàm truyền của hệ thống(lưu trong ct6_29.m):  
ts = [1 2 3 2 1 4];  
ms = [1 2 3 2 3]; 
[h,w] = freqz(b,a,64); 
[tsm,msm] = invfreqz(h,w,4,5) 
Ta cũng có thể xây dựng đặc tính tần thực‐ảo  
Ví dụ: Cho hàm truyền : 
10s9s5.4s
10)s(G 23 +++=  
Tìm đặc tính tần thực ‐ ảo của hệ bằng các lệnh MATLAB (lưu trong ct6_30.m): 
ts = [10]; 
ms = [1 4.5 9 10]; 
w = [1:0.01:3]; 
136
h = freqs(ts,ms,w); 
t = real(h); 
a = imag(h); 
subplot(2,1,1); 
plot(w,t) 
subplot(2,1,2); 
plot(w,a) 
  Để vẽ đồ thị Bode của hệ thống ta dùng hàm bode. Đồ thị thứ nhất nhất 
là đặc tính biên‐tần logarit, được chia theo dB. Đồ thị thứ hai là đặc tính pha‐ 
tần logarit chia theo độ.  
  Các dạng của lệnh bode gồm: 
  bode(sys) 
  bode(sys,w) 
  [bien, pha, w] = bode(sys) 
Để  vẽ  đồ  thị Bode  của một  hệ  thống  ta dùng  các  lệnh MATLAB(lưu  trong 
ct6_31.m) như sau: 
clc 
g = tf([1 0.1 7.5],[1 0.12 9 0 0]); 
figure(1) 
bode(g) 
figure(2) 
bode(g,{0.1 , 100}) 
gd = c2d(g,0.5) 
figure(3) 
bode(g,ʹrʹ,gd,ʹb‐‐ʹ) 
Hàm margin cho biết dự  trữ ổn định của hệ  thống. Dự  trữ biên gm  là hệ số 
khuyếch đại Fr mà nếu ta thêm vào hàm truyền đạt của hệ hở thì hệ kín vừa 
đạt được giới hạn ổn định. Dự trữ pha pm được định nghĩa là khoảng cách góc 
pha ϕr tới ‐180°. Hàm cho biết gm tại tần số đảo pha wcg và pm tại tần số cắt 
pha wcp. Hàm allmargin có tác dụng rộng hơn hàm margin. Các kết quả trả về 
của allmargin gồm: 
GMFrequency: giá trị tần số mà tại đó đồ thị pha cắt đường thẳng nằm 
ngang ‐180° 
GainMargin:  dự  trữ  biên  ‐  giá  trị  đảo  của  biên  độ  tại  tần  số 
GMFrequency 
PMFrequency: giá trị tần số mà tại đó đồ thị biên cắt đường thẳng nằm 
ngang 0 dB(ứng với hệ số khuyếch đại 1) 
137
PhaseMargin: dự trữ pha ‐ khoảng cách góc (> 0) từ vị trí PMFrequency 
đến ‐180°. 
          DelayMargin: dự trữ thời gian trễ ‐ giá trị thời gian trễ mà nếu vượt quá, 
hệ thống sẽ mất ổn định. 
DMFrequency: giá trị tần số ứng với DelayMargin.   
Stable: =1 khi mach vòng kín ổn định; bằng 0 trong các trường hợp khác. 
Các đại lượng này có thể đọc được từ đồ thị tạo bởi margin. Để xác định 
dự trữ ổn định của một hệ thống cụ thể ta dùng các lệnh MATLAB(lưu trong 
ct6_32.m) như sau: 
clc 
sys = zpk([],[‐1 ‐1 ‐1],4) 
margin(sys) 
allmargin(sys) 
Kết quả hệ thống ổn định. Nó có DelayMargin = 0.3s. Bây giờ ta gán cho sys 
một khoảng thời gian trễ  là stabil.DelayMargin + 0.01, nghĩa là vượt quá thời 
gian trễ ổn định 0.01s. Kết quả tính toan mới của allmargin sẽ thông báo tính 
không ổn định của hệ thống. Các lệnh MATLAB (lưu trong ct6_33.m) như sau: 
clc 
sys = zpk([],[‐1 ‐1 ‐1],4) 
margin(sys) 
stabil = allmargin(sys) 
sys.ioDelay = stabil.DelayMargin + 0.01; 
newstabil = allmargin(sys) 
Một khả năng khác để mô tả đặc tính tần số là đồ thị Nyquist. Nó biểu 
diễn các giá trị thực và ảo thuộc hàm truyền đạt phức của mạch vòng hở F0(jω) 
trong dải tần số ω = 0 ÷ ∞ trên hệ toạ độ phức. Đường cong do các điểm tạo 
thành  được gọi  là quỹ đạo biên  ‐ pha F0(jω). Trên  cơ  sở  tiêu  chuẩn  ổn  định 
Nyquist ta có thể rút ra kết luận về tính ổn định của hệ kín(có phản hồi đơn vị 
âm) từ đồ thị Nyquist. Để vẽ đồ thị Nyquist ta dùng hàm Nyquist. Ta có các 
lệnh MATLAB(lưu trong ct6_34.m) như sau: 
clc 
H = tf([2 5 1],[1 2 3]) 
nyquist(H) 
138
§3. ĐẶC TÍNH CỦA HỆ THỐNG ĐIỀU KHIỂN 
1. Tính ổn định: Tiêu chuẩn ổn định nói rằng hệ sẽ ổn định nếu các nghiệm 
của phương trình đặc tính có phần thực âm. Phương trình đặc tính là đa thức 
mẫu số của hàm truyền. Do vậy chỉ cần tính nghiệm của đa thức đặc tính bằng 
lệnh roots là ta có thể xác dịnh hệ ổn định hay không. 
Ví dụ: Xét tính ổn định của hệ có phương trình đặc tính là: 
  s4 + 10 s3 + 35s2 + 50s + 24 
Các lệnh MATLAB là: 
a = [1 10 35 50 24]; 
roots(a) 
ans = 
   ‐4.0000 
   ‐3.0000 
   ‐2.0000 
   ‐1.0000  
Như vậy hệ ổn định. 
2. Độ nhạy: Độ nhạy của hệ thống được đo bằng tỉ số phần trăm sự thay đổi 
của hàm truyền theo sự thay đổi phần trăm của thông số b. Ví dụ độ nhạy của 
hàm truyền T(s) theo b được xác định bằng: 
b
)s(T
b
)s(T
b/b
)s(T/)s(TSTb ∆
∆=∆
∆=  
Khi ∆b gần đến 0 ta có: 
)s(T
b
b
)s(TSTb ∂
∂=  
Độ nhạy tĩnh là giá trị của S khi t→0. Độ nhạy động được tính bằng cách thay s 
bằng jω và vẽ đường S theo ω. Biên độ của S(jω) đo sai số của hệ thống. 
Ví dụ: Khảo sát hệ điều khiển như hình vẽ sau: 
K
Thiết bịBộ bù
)1s(
b
+
Sensor
h
- 
R(s) 
C(s)  
139
Trong đó b có trị định mức là 4 và h có trị định mức là 0,5. Tìm độ nhạy T(s) 
theo b, vẽ modul hàm độ nhạy theo ω với hai giá trị bù là K = 2 và K = 0.5. Tìm 
độ nhạy T(s) theo h, vẽ modul của hàm độ nhạy theo h với K = 2 và K = 0.5. 
Hàm truyền của hệ thống là: 
Kbh1s
Kb)Ts( 2 ++=  
Với b = 4 và h = 0.5 ta có ωB = 1 + 2K. 
Độ nhạy của T(s) theo b khi b = 4 và h = 0.5 là: 
K21s
1s
Kbh1s
1s
)s(T
b
b
)s(TSTb ++
+=++
+=∂
∂=  
K21s
K2
Kbh1s
Kbh
)s(T
h
b
)s(TSTh ++
−=++
−=∂
∂=  
Các lệnh MATLAB (lưu trong ct6_35.m) như sau: 
k1 = 1; 
k2 = 0.5; 
ts = [1 1]; 
ms1 = [1 1+2*k1]; 
ms2 = [1 1+2*k2]; 
w = 0:0.01:15; 
stb1 = abs(freqs(ts,ms1,w)); 
stb2 = abs(freqs(ts,ms2,w)); 
subplot(2,1,1); 
plot(w,stb1,w,stb2); 
title(ʹDo nhay cua T theo bʹ); 
ts1 = ‐2*k1; 
ts2 = ‐2*k2; 
stb1 = abs(freqs(ts1,ms1,w)); 
stb2 = abs(freqs(ts2,ms2,w)); 
subplot(212); 
plot(w,stb1,w,stb2); 
title(ʹDo nhay cua T theo hʹ); 
Hình trên đường có ứng với K = 0.5 và hình dưới đường cao ứng với K = 2. 
  Từ hình vẽ  ta  thấy  rằng  độ nhạy  của hệ  thống  theo b giảm khi hệ  số 
khuếch đại của vòng hở K tăng trong khi độ nhạy theo h tăng khi K tăng. Rõ 
ràng là độ nhạy theo b tăng nhanh bên ngoài ωB. 
3. Sai số xác lập: Khảo sát hệ như hình vẽ: 
140
 R(s)
H(s)
G(s)
-
  C(s)
Hàm truyền của hệ kín là: 
)s(G)s(H1
)s(G
)s(R
)s(C
+=  
Sai số của hệ kín là: 
    E(s) = R(s) – H(s)C(s) = 
)s(G)s(H1
)s(R
+  
Sử dụng định lí giá trị cuối ta có: 
)s(H)s(G1
)s(sRlime
sss += ∞→  
Đầu vào bước nhảy đơn vị: 
ps
ss K1
1
)s(H)s(Glim1
1e +=+=
∞→
Đầu vào tăng tuyến tính đơn vị: 
vs
ss K
1
)s(H)s(sGlim1
1e =+=
→∞
Đầu vào parabol đơn vị: 
a
2
s
ss K
1
)s(H)s(Gslim1
1e =+=
→∞
Ta có thể dùng Symbolic Math để tính các giới hạn trên. 
§4. PHÂN TÍCH VÀ THIẾT KẾ QUỸ ĐẠO NGHIỆM 
 Phương pháp kinh điển để tham số hoá khâu điều khiển của vòng điều 
hỉnh là phương pháp quỹ đạo nghiệm. Quỹ đạo nghiệm là quỹ đạo điểm cực, 
hợp  thành bởi  các  điểu  cực  của hệ  thống, phụ  thuộc vào hệ  số khuyếch  đại 
phản hồi k va được biểu diễ trên mặt phẳng phức với phần thưc Re(λ) = σ trên 
trục  hoành  x  và  phần  ảo  Im(λ)  = ω  trên  trục  tung  y.  Để  vẽ  được  quỹ  đạo 
nghiệm của hệ thống ta dung hàm rlocus. Ta xét hệ thống sau: 
cG
k
0G
MG
-
u  y
141
Cú pháp của rlocus là 
  rlocus(sys[,k]) 
  [r, k] = rlocus(sys) 
  r = rlocus(sys, k) 
  Mô hình sys trong lệnh trên là hàm truyền đạt của hệ thống hở GoGcGM 
được xác định bằng lệnh MATLAB: 
  sys = sysM*sysO*sysC 
mà chưa có hệ số khuyếch đại phản hồi k,  là  tham số tuỳ chọn sẽ được khai 
báo riêng. Điều đó có nghĩa là sys được ghép nối bởi các mô hình riêng lẻ. Khi 
gọi  rlocus(sys[, k]) mà không yêu  trả biến về  ta nhận  được  đồ  thị quỹ  đạo 
nghiệm  của  sys. Nếu  ta không khai báo  các hệ  số khuyêch  đại  trong vec  tơ 
tham số tuỳ chọn k, MATLAB sẽ tự động quyết định giá trị thích hợp. Sau khi 
dùng rlocus vẽ quỹ đạo điểm cực ta tìm các giá trị liên quan đến điểm cực bất 
kì năm tên quỹ đạo bằng cách nhấp chuột vào một điểm trên quỹ đạo. Lúc đó 
lệnh  rlocusfind  được  thưc  hiện.  Ta  dùng  các  lện MATLAB  sau  (lưu  trong 
ct6_36.m)để vẽ quỹ đạo nghiệm của một hệ thống: 
clc 
sys = zpk([],[‐0.1 ‐1‐j ‐1+j ], 1) 
rlocus(sys) 
[r, k] = rlocus(sys) 
sgrid 
  Để  trực  quan  ta  có  thể  dùng  công  cụ  thiết  kế  bằng  cách  nhâp  lệnh 
sisotool vào cửa sổ lệnh MATLAB. 
142