Giáo trình Matlab cơ bản

   ‐0.7814;0.7814  0]; 
c = [1.9691  6.4493]; 
 506
x0 = [1 ; 0] 
sys = ss(a, [], c, []); 
initial(sys, x0) 
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 ctimpulse.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 ctstep2.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) 
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 ctlsim.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) 
 507
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ỳ. 
11. Đá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. 
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 
ctfreqs.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 ctfreqplot.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 
 508
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 ctinvfreqz.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 
ctfreqsplot.m): 
ts = [10]; 
ms = [1 4.5 9 10]; 
w = [1:0.01:3]; 
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 
ctbode.m) như sau: 
 509
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) 
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 
ctmatgin6_32.m) như sau: 
clc 
sys = zpk([], [‐1 ‐1 ‐1], 4) 
margin(sys) 
 510
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 
ctnewstabil6_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 ctnyquist6_34.m) như sau: 
clc 
H = tf([2 5 1], [1 2 3]) 
nyquist(H) 
12. 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]; 
 511
roots(a) 
ans = 
   ‐4.0000 
   ‐3.0000 
   ‐2.0000 
   ‐1.0000  
Như vậy hệ ổn định. 
13. Độ 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: 
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à: 
K
)1s(
b
+ 
h
R(s)  C(s) 
- 
Bộ bù Thiết bị
Sensor
 512
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 ctsensibility.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ʹ); 
Độ 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. 
14. Sai số xác lập: Khảo sát hệ như hình vẽ: 
Hàm truyền của hệ kín là: 
)s(G)s(H1
)s(G
)s(R
)s(C
+=  
H(s)
G(s)R(s) C(s)
-
 513
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. 
15. 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: 
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 
cG 
k
0G 
MG 
u y 
-
 514
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ệnh MATLAB  sau  (lưu  trong 
ctrlocus.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.