🔙 Quay lại trang tải sách pdf ebook Matlab Lecture Ebooks Nhóm Zalo 9/16/2015 Lecture on Matlab Bài mở đầu Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology 1 Chương 1. Môi trường Matlab 1. Matlab là gì ? • MATLAB được viết tắt từ “MATrix LABoratory”, là một công cụ tính toán số và mô phỏng số, công cụ này được phát triển dựa trên các thư viện hàm tính toán số viết bằng ngôn ngữ lập trình FORTRAN. Matlab có giao diện thân thiện với người sử dụng. • Một định nghĩa khác: Matlab là một phần mềm đa năng tính toán số, thể hiện các số liệu bằng hình ảnh trực quan và là một ngôn ngữ đa năng cung cấp một môi trường linh hoạt cho việc tính toán kỹ thuật. • Matlab bao gồm nhiều công cụ để: - thu thập dữ liệu, phân tích và xử lý dữ liệu, - hiển thị hình ảnh trực quan và xử lý hình ảnh, - tạo mẫu và phát triển thuật toán - mô hình hóa và mô phỏng, lập trình và phát triển ứng dụng 3 Chương 1. Môi trường Matlab 1. Matlab là gì ? 2. Giao diện người sử dụng 3. Các phép tính số học cơ bản 4. Phép gán 5. Các định nghĩa toán học cơ bản 6. Số phức 7. Xử lý lỗi khi gõ lệnh 8. Kết thúc phiên làm việc với Matlab 2 Chương 1. Môi trường Matlab • Matlab có thể chạy trên hầu hết các hệ máy tính: máy tính sách tay - Labtop, máy tính cá nhân PC, đến các hệ thống siêu máy tính (super computer). Matlab được điều khiển bởi các tập lệnh, tác động qua bàn phím trên cửa sổ điều khiển. Nó cũng cho phép một khả năng lập trình với cú pháp thông dịch lệnh – còn gọi m file. Các lệnh hay bộ lệnh của Matlab lên đến con số hàng trăm và ngày càng được mở rộng nhờ các thư viện trợ giúp hay do người sử dụng tạo ra. 4 1 9/16/2015 Giao diện người sử dụng 1. Cửa sổ lệnh (Command 4 Window). Phép gán • Trong Matlab dấu bằng “=” được sử dụng cho phép gán. Mặc dù trong cách viết thông thường ta hiểu dấu bằng thể hiện một phương trình, nhưng trong Matlab dấu bằng được định nghĩa là phép gán. Để phân biệt giữa hai cách thể 2. Cửa sổ ghi lại các lệnh 5 đã thực hiện (Command History). 3 3. Cửa sổ cho biết thư mục hiện thời (Current Directory) và danh sách 1 các biến đang sử dụng 2 (Workspace). 4. Thanh chỉ đường dẫn của thư mục hiện thời, cho 6 phép chọn thư mục hiện thời. 5. Thanh Shortcut. 6. Nút Start. 5 hiện ta xét ví dụ sau. Nếu trong cửa sổ lệnh ta viết >> x + 18 = 120 • Matlab sẽ báo lỗi với dòng hiển thị mầu đỏ: ??? x+18=120 Error: The expression to the left of the equals sign is not a valid target for an assignment. • Matlab không hiểu cách bạn viết một phương trình như trên giấy, mà nó chỉ có thể thực hiện được các phép tính và gán giá trị tính được cho một biến nào đó. Chẳng hạn ta có thể gán giá trị của phép tính 120-18 cho biến x bằng cách viết >> x=120 - 18; • Các dấu “ ; ” ở cuối dòng lệnh báo cho Matlab biết là ta không muốn các giá trị của x hiện ra ở dòng dưới. >> x=120 - 18 6 Các phép tính số học cơ bản Các phép tính số học cơ bản • Đối với người mới bắt đầu sử dụng Matlab, trước hết ta cần làm quen với các phép tính số học (cộng +, trừ -, nhân * và chia /). Với hai số a và b trong Matlab ta có thể thực hiện được các phép tính số học như liệt kê trong bảng dưới đây. Phép tính Cách viết toán học Cách viết trong Matlab phép cộng c = a + b c = a + b phép trừ c = a – b c = a – b phép nhân c = ab c = a*b phép chia c = a : b c = a / b phép chia phải c = a : b c = a / b phép chia trái c = b : a c = a \ b lũy thừa c = ab c = a^b 7 >> a=3.5; b=7.5; % gan gia tri so cho hai bien a va b >> a+b % phep cong hai so ans = 11 >> a-b % phep tru hai so ans = -4 >> a*b % phep nhan hai so ans = 26.2500 >> a/b % phep chia hai so, a:b (= phep chia phai) ans = 0.4667 8 >> a/b % phep chia hai so, a:b (= phep chia phai) ans = 0.4667 >> a\b % phep chia trai (hieu la b chia cho a, b/a) ans = 2.1429 >> a^b % phep luy thua, a mu b ans = 1.2037e+004 2 9/16/2015 Một số dạng format • format long, short, rat, short e, bank Các định nghĩa và hàm toán học cơ bản • Để thuận tiện cho việc tính toán, trong Matlab người ta đã định nghĩa sẵn rất • Ví dụ >> format long >> co = cos(0.2), si = sin(0.2), tn = tan(0.2) co = 0.98006657784124 si = 0.19866933079506 tn = 0.20271003550867 >> format short >> co = cos(0.2), si = sin(0.2), tn = tan(0.2) co = 0.9801 si = 0.1987 tn = 0.2027 >> format short e >> x=5.125*3.16 x = 1.6195e+001 >> format long e >> x=5.125*3.16 x = 1.619500000000000e+001 >> format rat >> x=5.125*3.16 x = 3239/200 >> format bank >> luonggio=35.55 luonggio = 35.55 >> luongtuan=luonggio*40 luongtuan = 1422.00 9 Các định nghĩa và hàm toán học cơ bản nhiều đại lượng toán học và các hàm cơ bản. Ví dụ như số pi đã được định nghĩa sẵn >> R = 2; >> V = 4/3*pi*R^3 V = 33.5103 >> exp(1) ans = 2.7183 >> exp(2) ans = 7.3891 >> x = sqrt(9) x = 3 >> y = sqrt(11) y = 3.3166 >> log(3.2) % ln(3.2) ans = 1.1632 >> x = 5; log(x) ans = 1.6094 >> x = 3; log10(x) ans = 0.4771 10 Các định nghĩa và hàm toán học cơ bản Các hàm toán học sqrt(x) Tính căn của biến x exp(x) Hàm e mũ, Exponential log(x) Logarithm cơ số tự nhiên, cơ số e log1p(x) Tính giá trị log(1+x) log10(x) Logarithm cơ số 10 log2(x) Logarithm cơ số 2 pow2(x) Hàm mũ cơ số 2, pow2(x) = 2^x realpow(x,y) Tính x^y, với x, y là thực reallog(x) Logarithm cơ số tự nhiên của số thực nthroot(x, n) Căn bậc n của số thực nextpow2(n) Trả lại số p đầu tiên sao cho 2^p >= abs(n) 11 Các hàm lượng giác - Trigonometric functions sin Hàm sin >> sin(pi/3) cho 0.8660 sind Hàm sin với đối số là độ >> sind(60) cho 0.8660 sinh Hàm sin Hyperbolic asin Hàm sin ngược, tức arcsin >> asin(0.866) cho 1.0471 asind Hàm sin ngược cho kết quả là độ >> asind(0.866) cho 59.9971 asinh Hàm sin Hyperbolic ngược cos Hàm cos cosd Hàm cos với đối số là độ e e− x x − cosh Hàm cos Hyperbolic sinh( )x= acos Hàm cos ngược, tức arccos 2 acosd Hàm cos ngược cho kết quả là độ x xe e−+ acosh Hàm cos Hyperbolic ngược cosh( )2x= 12 3 9/16/2015 Số phức Chúng ta cũng có thể đưa số phức vào trong Matlab. im Nhớ lại rằng đơn vị phức được định nghĩa là căn của (-1), Trong Matlab số phức ký hiệu là i hoặc j • Một số phức có thể được viết dạng z = x + i b trong đó x là phần thực và y là phần ảo của z re • Ví dụ a = 3 + 3i, b = 1 – i 🡺 a + b = 4 + 2i Số phức • Một số hàm liên quan đến số phức được liệt kê trong bảng dưới đây Các hàm với số phức abs Cho độ lớn hay modul véctơ phức angle Cho góc pha hay argument của số phức, radian complex Tạo lập dữ liệu phức từ các phần thực vào ảo conj Cho số phức liên hợp imag Cho phần ảo real Cho phần thực isreal Trả lại giá trị đúng, 1, nếu đối số là thực cplxpair Sắp xếp các số phức thành cặp liên hợp >> a = 3 + 3i, b = 1 - i; >> c = a + b c = 4.0000 + 2.0000i 2 2 | | ( ) ( ) ( ) a abs a re a im a arg( ) arctan( ( ) / ( )) a im a re a | | , i a a ecos sin i e i >> abs(a) 🡪 4.2426 = 2 sqrt(3) >> angle(a) 🡪 0.7854 = pi/4 14 13 Xử lý lỗi khi gõ lệnh, Kết thúc phiên làm việc với Matlab • Như đã thấy ở trên, khi thực hiện các dòng lệnh bạn có thể làm không đúng và Matlab đã báo lỗi cho bạn. Error !. Nếu khi bạn kết thúc dòng lệnh bằng gõ phím ENTER và nhận ra rằng dòng lệnh trên có lỗi, bạn không nhất thiết phải gõ lại dòng lệnh đó, mà đơn giản bạn chỉ cần sử dụng các phím mũi tên lên hoặc xuống, ↓ ↑ , để hiển thị lại dòng lệnh cần sửa. Sau khi sửa lỗi, và đánh phím ENTER Matlab sẽ đưa ra kết quả. • Chúng ta vừa bắt đầu với một số lệnh cơ bản của Matlab, và bây giờ có thể bạn muốn ghi lại các công việc vừa thực hiện và thoát khỏi Matlab. Làm thế nào để kết thúc Matlab trên màn hình của bạn? Thật đơn giản bạn vào exit trong menu File, như khi sử dụng các chương trình khác trong Windows. Một cách khác là bạn gõ lệnh quit vào dấu nhắc trong cửa sổ lệnh và Matlab sẽ đóng lại. Kết thúc phiên làm việc. 15 16 4 Lecture 2 Véctơ, ma trận và Matlab Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology [email protected] 9/16/2015 Chương 2. Véctơ, ma trận và Matlab 2.1 Véctơ và các phép tính trên véctơ 2.2 Biểu diễn đa thức và các phép tính đa thức 2.3 Ma trận và các phép tính cơ bản trên ma trận 2.1 Véctơ và các phép tính trên véctơ Nhập véctơ Các phép tính cộng và trừ hai véctơ cùng cỡ Tạo một véctơ từ các véctơ khác Tạo véctơ hàng có các phần tử cách đều Các đặc trưng của véctơ Tích vô hướng và tích có hướng hai véctơ Tham chiếu đến các phần tử của véctơ Một số hàm định sẵn cho véctơ Véctơ và các phép tính trên véctơ Véctơ và các phép tính trên véctơ • Nhập véctơ >> a= [2; 1; 4] % vt cot >> v= [2 0 4] % vt hang 2 a y a 1 ' 2 1 4 Tạo một véctơ từ các véctơ khác Tạo véctơ hàng có các phần tử cách đều (toán tử hai chấm : ) x = [a1 : h : a2] 4 >> w= [1,1,9] % sd dau phay • Chuyển vị vectơ (cột thành hàng, hàng thành cột) >> a = [2; 1; 4]; >> y = a' • Các phép tính cộng và trừ hai véctơ cùng cỡ >> a = [1; 4; 5]; 1 2 a b 4 , 3 5 3 3 1 Ghép hai véc tơ cột >> a = [1; 4; 5]; >> b = [2; 3; 3]; >> d = [a; b] Ghép hai véc tơ hàng >> a = [1 4 5]; >> b = [2 3 3]; d = 1 4 5 2 3 3 x(i) = a1 + (i - 1)h, i = 1, 2, 3, . . . >> x = [0:2:10] x = 0 2 4 6 8 10 >> x = [0:0.1:1] x = 0 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 >> b = [2; 3; 3]; >> c1 = a + b >> c2 = a - b c a b c a b 7 , 1 1 2 8 2 >> d = [a b] >> help linspace d = 1 4 5 2 3 3 1 Véctơ và các phép tính trên véctơ • Cho các vectơ cùng cỡ v và u ta có các hàm tác động lên véc tơ Tên hàm Ý nghĩa length(v) Cho biết chiều dài, số phần tử của véctơ v sum(v) Trả lại tổng các phần tử của véctơ v prod(v) Trả lại tích các phần tử của véctơ v min(v) Trả lại phần tử nhỏ nhất của véctơ v max(v) Trả lại phần tử lớn nhất của véctơ v sort(v) Xếp các phần tử tăng dần find(v) Tìm các phần tử khác không trong v dot(u,v) Trả lại tích vô hướng hai véctơ cross(u,v) Trả lại tích có hướng hai véctơ cỡ 3*1 9/16/2015 Véctơ và các phép tính trên véctơ Ví dụ >> a = [8 4 4 1 7 11 2 0]; >> a = [2;3;3;4;5]; >> length(a) ans = 5 >> b = [1, 1, 2]; >> length(b) ans = 3 >> a = [8 4 4 1 7 11 2 0]; >> max(a) ans = 11 >> sum(a) ans = 37 >> prod(a) ans = 0 >> u = [1; 4; 5]; v = [2; 3; 3]; >> x = dot(u,v) x = 29 >> w = cross(u,v) >> min(a) ans = 0 w = -3 7 -5 Véctơ và các phép tính trên véctơ Biểu diễn đa thức và các phép tính đa thức Tham chiếu đến phần tử vectơ, Toán tử : >> a = [12; 17; –2; 0; 4; 27]; >> a(2) ans = 17 >> a(6) ans = 27 >> a(:) ans = 12 17 –2 0 4 27 >> v = a(1:3) v = 12 17 –2 2.2 Biểu diễn đa thức và các phép tính đa thức Nhập đa thức Các phép tính trên đa thức Phép nhân đa thức / Phép chia đa thức Phép cộng và trừ đa thức Không điểm hay nghiệm của phương trình đa thức Xây dựng đa thức từ các không điểm cho trước Giá trị của đa thức tại một điểm Đạo hàm đa thức 2 Biểu diễn đa thức và các phép tính đa thức • Nhập đa thức Trong Matlab một đa thức được đưa vào dưới dạng một véctơ hàng có các phần tử là các hệ số của đa thức. Khi đưa vào hay khai báo một đa thức ta cần biểu diễn dạng đầy đủ của nó, ví dụ xét đa thức 5 5 4 2 1 0 p s a s a s a s a s a ( ) = + + + + 5 4 2 với dạng đầy đủ như sau 5 4 3 2 1 0 5 5 4 2 1 0 p s a s a s s a s a s a s ( ) 0 = + + + + + Véc tơ ứng với đa thức bậc 5 trên là một vectơ hàng gồm sáu phần tử 5 5 4 2 1 0 p a a a a a = [ 0 ] Ví dụ 9/16/2015 Biểu diễn đa thức và các phép tính đa thức • Matlab sẽ xác định bậc của đa thức từ số phần tử của véctơ, trong ví dụ trên >> n = length(p)-1 n = 5 • Các hàm liên quan đến đa thức p = conv(p1,p2) Nhân hai đa thức [q,r]=deconv(p1,p2) Phép chia đa thức p1 cho p2, bậc của p1 lớn hơn của p2 roots(p) Tìm các nghiệm của phương trình đa thức, p(x) = 0 p = poly(nu) Xây dựng đa thức từ các không điểm cho trước, đưa racho ta một đa thức nhận các giá trị trong véctơ nu làm nghiệm. polyval(p,s0) Giá trị của đa thức tại một điểm polyder(p) Đạo hàm đa thức Phép cộng và trừ đa thức ?? >> p = [8 3 0 -4 6 0] p = 8 3 0 -4 6 0 5 p s s s s s ( ) 8 3 4 6 = + − + 5 4 2 Biểu diễn đa thức và các phép tính đa thức >> nu = [-2 -1 0 1 2] >> poly(nu) Biểu diễn đa thức và các phép tính đa thức • Chia đa thức - Phép chia đa thức p1 cho p2, bậc của p1 lớn hơn của p2 [ q, r ] = deconv(p1, p2) ans = 5 3 p s s s s ( ) 5 4 = − + 54 2 = − + s s s p s r s ( ) ( ) ( ) 1 0 -5 0 4 0 ( 5 4) 1 q s p s p s ( ) ( ) 2 2 • polyder(p) Đạo hàm đa thức n n k k 1 p s a s p s ka s ( ) ( ) k k k k 0 1 3 Ma trận và các phép tính cơ bản trên ma trận Ma trận và các phép tính cơ bản trên ma trận Nhập ma trận Nhân ma trận với một số, phép cộng và trừ hai ma trận Chuyển vị ma trận Phép nhân ma trận 9/16/2015 Ma trận và các phép tính cơ bản trên ma trận • Nhập ma trận Ma trận là một mảng số hai chiều, được sắp xếp theo hàng và cột. Để tạo một ma trận trong Matlab, chúng ta đặt từng hàng vào trong dấu ngoặc vuông [], các hàng được phân biệt với nhau bằng dấu chấm phẩy (;), và trong mỗi hàng các phần tử được phân biệt với nhau bằng dấu cách hoặc bằng dấu phẩy. Một số phép tính cơ bản khác • Ví dụ với hai ma trận 1 6 ⎡ ⎤ 2 0 1 Các ma trận đặc biệt Tham chiếu đến các phần tử của ma trận Ma trận khối ⎡ ⎤ − ⎢ ⎥ A B = = − ⎢ ⎥ ⎢ ⎥ , 1 7 4 7 11 ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ 3 0 1 ⎢ ⎥ ⎣ ⎦ Các phép nhân và chia phần tử hai ma trận cùng cỡ Tính định thức và giải hệ phương trình đại số tuyến tính Tìm hạng của ma trận Tìm ma trận nghịch đảo và ma trận tựa nghịch đảo (pseudo-inverse) Trị riêng và véctơ riêng của ma trận vuông Phân tích ma trận vuông A thành tích các ma trận • Ta đưa vào Matlab như sau: >> A = [-1, 6; 7, 11] A = -1 6 7 11 >> B = [2,0,1; -1,7,4; 3,0,1] B = 2 0 1 -1 7 4 3 0 1 Ma trận và các phép tính cơ bản trên ma trận Ma trận và các phép tính cơ bản trên ma trận Nhân ma trận với một số >> A = [-3, 6; 7, 10] ; Chuyển vị ma trận ⎡ ⎤ − 1 2 1 3 7 Phép nhân hai ma trận có cỡ thích hợp, C = A*B = = = { }, 1.. , 1.. , >> B = [5, 1; -1, 3] ; ⎢ ⎥ ⎡ ⎤ T A A = = ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦ A a i m j p ij >> C=2*A 3 5 , 2 5 4 = = = { }, 1.. , 1.. B b i p j n C = -6 12 7 4 ij p 14 20 Sử dụng dấu nháy đơn (’) C = AB = = = = ∑ { }, , 1.. , 1.. , c c a b i m j n ij ij ik kj k = 1 Cộng và trừ hai ma trận cùng cỡ >> A = [–1 2 0; 6 4 1] >> A = [2 1; 1 2]; B = [3 4; 5 6]; >> A = [1 4; 8 0; -1 3]; >> A+B A = –1 2 0 >> C = A*B >> B = [-1 7 4; 2 1 -2]; ans = 2 7 6 13 >> A-B ans = -8 5 8 7 6 4 1 >> B = A' B = –1 6 2 4 0 1 C = 11 14 13 16 >> C = A*B C = 7 11 -4 -8 56 32 7 -4 -10 4 Ma trận và các phép tính cơ bản trên ma trận Phép nhân phần tử hai ma trận .* { }, * , 1.. , 1.. ij ij ij ij C = A B = = = = c c a b i m j n Cần phân biệt phép nhân phần tử ( .* ) hai ma trận cũng cỡ C = A .* B 9/16/2015 Ma trận và các phép tính cơ bản trên ma trận • Các ma trận đặc biệt Khởi tạo các ma trận đặc biệt eye(m,m) eye(m) eye(m,n) tạo ma trận đơn vị cỡ mxmtạo ma trận đơn vị cỡ mxm tạo ma trận đơn vị mở rộng cỡ mxn zeros(m,n) tạo ma trận không cỡ mxn zeros(m,m) zeros(m) tạo ma trận không cỡ mxm, hay ma trận vuông ones(m,n) cho ta một ma trận 1 cỡ mxn ones(m,m) ones(m) cho ta một ma trận 1 cỡ mxm, hay ma trận vuông magic(n) tạo ma trận magic vuông cỡ nxn, n > 2 [ma phương] >> A = [2 1; 1 2]; >> B = [3 4; 5 6]; >> C = A .* B C = 6 4 5 12 >> A = [1 4; 8 0; -1 3]; B = [-1 7 4; 2 1 -2]; >> C = A * B C = 7 11 -4 -8 56 32 7 -4 -10 >> C = A .* B ??? Error using ==> times Matrix dimensions must agree. Ma trận và các phép tính cơ bản trên ma trận Các ma trận đặc biệt >> eye(3) ans = 1 0 0 0 1 0 0 0 1 >> eye(2,3) ans = 1 0 0 0 1 0 >> zeros(3,3) ans = 0 0 0 0 0 0 0 0 0 >> zeros(3,2) ans = 0 0 0 0 0 0 >> zeros(2) ans = 0 0 0 0 >> ones(2,2) or ones(2) ans = 1 1 1 1 >> magic(3) ans = 8 1 6 3 5 7 4 9 2 >> magic(5) ans = 17 24 1 8 15 23 5 7 14 16 4 6 13 20 22 10 12 19 21 3 11 18 25 2 9 Ma trận và các phép tính cơ bản trên ma trận Hàm diag • Lấy ra đường chéo của matrận bằng lệnh diag, và tạo ra ma trận đường chéo từ một vector bằng lệnh diag • c = diag(A) cho ta đường chéo của ma trận A • A = diag(c) cho ta ma trận đường chéo A, với A(i,i) = c(i) >> A1=[1 2 3; 3 4 5; 3 5 7] A1 = 1 2 3 3 4 5 3 5 7 >> Ad=diag(A1) Ad = 1 4 7 >> Add=diag(Ad) Add = 1 0 0 0 4 0 0 0 7 5 Ma trận và các phép tính cơ bản trên ma trận • Tham chiếu đến các phần tử của ma trận 9/16/2015 Ma trận và các phép tính cơ bản trên ma trận Ma trận khối Các phần tử, các cột hay các hàng của ma trận đều có thể được tác động đến nhờ cách đánh chỉ số của chúng. Chỉ số của các phần tử của ma trận là cặp số nguyên (i,j), i là chỉ số hàng và j là chỉ số cột. Các số này bắt đầu từ 1 (1,2,...), Matlab không sử dụng chỉ số 0 như một số ngôn ngữ lập trình khác. Ví dụ đối với ma trận • Cho các ma trận có kích cỡ thích hợp: A, B, C, M • Tạo ma trận cỡ lớn hơn từ các ma trận trên • G = [ A, B ; C , M ]; • Ví dụ, với A và B A B ⎡ ⎤ = ⎢ ⎥ G C M ⎢ ⎥ ⎣ ⎦ >> A = [1 2 3; 4 5 6; 7 8 9]; >> A(2,3) ans = 6 >> c2 = A(:,2) ans = 2 5 8 >> r3=A(3,:) r3 = 7 8 9 >> A(:,2:3) ans = 2 3 5 6 8 9 ⎡ ⎤ ⎡ ⎤ 2 4 4 2 A B = = ⎢ ⎥ ⎢ ⎥ , 7 5 5 7 ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ • Hãy xây dựng ma trận Q sau 0 E ⎡ ⎤ 2 2 2 2 × × Q A B BA = ⎢ ⎥ − 1 T ⎢ ⎥ −⎣ ⎦ >> A = [2 4; 7 5] >> B = [4 2; 5 7] >> Q = [ zeros(2) eye(2); -A\B B*A’ ] Q = 0 0 1 0 0 0 0 1 0 -1 16 38 -1 0 38 70 Ma trận và các phép tính cơ bản trên ma trận Định thức và giải hệ phương trình đại số tuyến tính • Định thức của ma trận vuông là một số. Ví dụ đối với ma trận vuông cỡ 2 × 2 Ma trận và các phép tính cơ bản trên ma trận Định thức và giải hệ phương trình đại số tuyến tính a a ⎡ ⎤ det( )a aa a a a >> A = [5 2 –9; –9 –3 2; 6 7 3]; A11 1211 22 12 21 11 12 = ⎢ ⎥ ⎣ ⎦ , a a 21 22 >> A = [1 3; 4 5]; >> det(A) ans = –7 A = = − a a 21 22 >> B = [3 –1 2 4; 0 2 1 8; –9 17 11 3; 1 2 3 –3]; >> det(B) ans = –533 >> b = [44;11;5]; >> det(A) ans = 368 • Xét hệ phương trình đại số • Viết lại dạng ma trận Ax = b >> x = A \ b tuyến tính sau ⎡ ⎤ − ⎡ ⎤ ⎡ ⎤ x = 5 2 9 x 44 –5.1250 ⎢ ⎥ = − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ x b = = ⎢ ⎥ ⎢ ⎥ , 11 5 2 9 44 x y z + − = − − + = 9 2 2 11 x y z + + = 6 7 3 44 x y z A 9 2 2 ⎢ ⎥ 6 7 3 ⎢ ⎥ ⎣ ⎦ y ⎢ ⎥ ⎢ ⎥ z 44 ⎣ ⎦ ⎢ ⎥ ⎣ ⎦ 7.6902 –6.0272 • Nếu det(A)≠0 pt Ax = b có nghiệm duy nhất x = A \ b 6 Ma trận và các phép tính cơ bản trên ma trận Tìm ma trận nghịch đảo • Nghịch đảo của ma trận vuông A là một ma trận, được ký hiệu là A-1 thỏa mãn A A-1 = A-1A = E, với E là ma trận đơn vị cùng cỡ. • Điều kiện để tồn tại ma trận nghịch đảo là det(A) ≠ 0. • Nếu tồn tại ma trận nghịch đảo của A, thì phương trình đại số tuyến tính A x = b có nghiệm duy nhất là x = A-1 b • Tính nghịch đảo ma trận vuông bằng lệnh inv(A) >> A = [2 3; 4 5] >> det(A) ans = -2 >> iA = inv(A) iA = -2.5000 1.5000 2.0000 -1.0000 >> S = [1 0 -1 2; 4 -2 -3 1; 0 2 -1 1; 0 0 9 8];>> det(S) ans = -108 >> iS=inv(S) iS = -0.9259 0.4815 0.4815 0.1111 -0.6296 0.1574 0.6574 0.0556 -0.5926 0.1481 0.1481 0.1111 0.6667 -0.1667 -0.1667 0 9/16/2015 Ma trận và các phép tính cơ bản trên ma trận • Để kiểm tra lại ta tính tích hai ma trận: >> iS*S ans = 1.0000 0 0 0 0 1.0000 0 0.0000 0 0 1.0000 0 0 0 0.0000 1.0000 >> S*iS ans = 1.0000 0 0 0 0.0000 1.0000 -0.0000 0 0.0000 -0.0000 1.0000 0 0 0 0 1.0000 • Ví dụ giải hệ A.x = b >> A = [3 –2; 6 –2]; b = [5; 2]; >> det(A) ans = 6 >> x=inv(A)*b x = -1.0000 -4.0000 >> e =A*x – b % thu lai e = 1.0e-015 * 0.0000 0.8882 Ma trận và các phép tính cơ bản trên ma trận Trị riêng và véctơ riêng của ma trận vuông • Cho A là một ma trận vuông cấp n, số λ được gọi là trị riêng và véctơ khác không x là véctơ riêng của A nếu chúng thoả mãn điều kiện Ax x = λ hay ( ) A E x 0 − = λvới E là ma trận đơn vị. • Cho A và B là hai ma trận vuông cấp n, nếu tồn tại số λ và véctơ x khác không thoả mãn điều kiện Ax Bx =λ hay ( ) 0 A B x − λ = thì số λ gọi là trị riêng suy rộng của hai ma trận A và B, véctơ x là véctơ riêng tương ứng. Ma trận và các phép tính cơ bản trên ma trận Trị riêng và véctơ riêng của ma trận vuông Lệnh eig tính toán trị riêng và véctơ riêng của ma trận vuông l = eig(A); cho một véctơ l chứa các trị riêng của ma trận vuông A. [V, D] = eig(A); cho ta ma trận chéo D chứa các trị riêng của A, ma trận V có các cột là các véctơ riêng của A thỏa mãn AV = VD. d = eig(A,B); cho ta véctơ chứa các trị riêng suy rộng của các ma trận vuông A và B. [V, D] =eig(A,B); cho ta ma trận chéo D chứa các trị riêng suy rộng của hai ma trận A và B, ma trận V có các cột là các véctơ riêng thỏa mãn AV = BVD. eig(A,B,'chol') tương tự như eig(A,B) đối với ma trận đối xứng A và ma trận đối xứng xác định dương B. Phương pháp tính ở đây dựa trên khai triển Cholesky ma trận B. 7 Ma trận và các phép tính cơ bản trên ma trận Trị riêng và véctơ riêng của ma trận vuông • Ví dụ >> C = [2 1 0; 1 4 -1; 0 -1 7]; >> eig(C) ans = 1.5573 4.1233 7.3194 >> [V, D]=eig(C) V = 0.9119 -0.4064 -0.0571 -0.4037 -0.8630 -0.3037 -0.0742 -0.3000 0.9510 D = 1.5573 0 0 0 4.1233 0 0 0 7.3194 9/16/2015 Ma trận và các phép tính cơ bản trên ma trận Phân tích LU và phân tích Cholesky • Phân tích LU. Với mỗi ma trận vuông A không suy biến, luôn phân tích được thành tích của hai ma trận tam giác A LU = với L là ma trận tam giác dưới, U là ma trận tam giác trên Ý nghĩa của phân tích LU có thể thấy trong việc giải hệ Ax = b. Sau phi phân tích LU ta nhận được Ax LUx b = = Đặt y Ux Ly b = ⇒ =ta giải hai hệLy b =vàUx y = Do L và U là các ma trận tam giác nên việc giải hai hệ này khá đơn giản. • Phân tích Cholesky. Ma trận vuông A đối xứng xác định dương luôn phân tích được thành tích của hai ma trận T A LL = với L là ma trận tam giác dưới Ma trận và các phép tính cơ bản trên ma trận các phép tính cơ bản trên ma trận Ví dụ phân tích LU >> A = [3 2 -9; -9 -5 2; 6 7 3]; b = [-65; 16; 5]; >> [L, U, P] = lu(A) • Ví dụ phân tích cholesky >> M = [2 3 4;3 5 6;4 6 9]; >> L = chol(M) L = .* .\ ./ .^ \ / Các phép tính ma trận, vector Cac hàm toán học tác động trên ma trận, vector Các phép tính trên từng phần tử Chia trái và chia phải L = -0.3333 0.0909 1.0000 1.0000 0 0 -0.6667 1.0000 0 U = -9.0000 -5.0000 2.0000 0 3.6667 4.3333 0 0 -8.7273 >> x = U\(L\b) x = 2.0000 -4.0000 7.0000 1.4142 2.1213 2.8284 0 0.7071 0.0000 0 0 1.0000 >> L‘ * L - M ans = 1.0e-015 * 0.4441 0 0 0 0 0 0 0 0 transpose(A), A’, A.’ ctranspose(A), A’, A.’ inv(A) det(A) linsolve(A, B) eig(A) rank(A) <[m, n]= > size(A <,i>) length(A) Chuyển vị (ma trận, vector) Chuyển vị, ma trận, vector liên hợp Nghịch đảo ma trận Tính định thức ma trận vuông Giải hệ AX=B, phân tích LU Trị riêng, vector riêng ma trận Hạng ma trận Cho biết cỡ ma trận, vector (i=1 số hàng, i=2 số cột) Giá trị lớn nhất của số hàng và số cột (=max(sizeA)) 8 9/16/2015 các phép tính cơ bản trên ma trận Các phép tính ma trận, vector Cac hàm toán học tác động trên ma trận, vector sum(v) prod(v) min(v) max(v) max(A(:)) sort(v) find(v) Tổng các phần tử của vector v Tích các phần tử của vector v Phần tử nhỏ nhất của vector v Phần tử lớn nhất của vector v Phần tử lớn nhất của ma trận A Xếp các phần tử tăng dần Tìm các phần tử khác không 9 Ma trận tựa nghịch đảo Ma trận tựa nghịch đảo [pseudo-inverse], help pinv 9/16/2015 Ma trận tựa nghịch đảo Nếu ma trận A suy biến, det(A) = 0, tức là trong ma trận A có những hàng là tổ Giải hệ A x = b; A là ma trận chữ nhật -Số ẩn > số phương trình -Số ẩn < số phương trình Ax b = hợp tuyến tính của các hàng khác. Do đó, ta có thể biến đổi hệ phương trình ban đầu về một hệ gồm m phương trình độc lập (m số phương trình: vô số nghiệm, do đó có thể đưa thêm điều kiện vào bài toán. Ví dụ cần tìm x sao cho độ lớn (chuẩn) của x nhỏ nhất Ax b A A × = ∈ ℜ > = , , , ( ) Hệ phương trình đại số tuyến tính Ax = b sẽ có vô số nghiệm. Điều này cho phép ta tìm được một nghiệm tối ưu theo một nghĩa nào đó từ tập vô số các nghiệm. m n = → x x 1 min n m rank m Đưa vào tiêu chuẩn tối ưu cho nghiệm cần tìm Tìm nghiệm của phương trình J T 2 Ax b J x x W x x 2 ( ) ( ) min T sao cho 1 o o Bài toán xuất hiện : -Robot dư dẫn động (số tọa độ khớp >6, robot kg; số tọa độ khớp >3, robot -Ma trận trọng số (xác định dương) W phẳng; [kinematics of redundant manipulator] -Tàu lặn mini (số cánh quạt >6). x 0 -Véc tơ tham chiếu (được chọn trước) 2 2 ( ) min T W E x x x J Ref. Mot so bai toan lien quan ma tran Ma trận tựa nghịch đảo x 0 x Wx J2 2 min T 1 o Ma trận tựa nghịch đảo [pseudo-inverse] x W A x 1 T Sử dụng PP nhân tử Lagrange để giải quyết bài toán tối ưu trên: W A AW A b Ax x [ ] ( ) λ1 1 [ ] ( ) T AW A b Axo ta sẽ tìm cực trị của hàm Lagrange như sau o 1 1 1 T T λ W A AW A b E W A AW A A x [ ] ( [ ] ) o o 2 ( , ) ( ) ( ) [ ] T T x x x W x x b Ax λ λ 1 1 1 1 1 1 1 o o m 1 λ T T T T o Đạo hàm theo biến x ta nhận được xλ λ ( , ) ( ) 0 T x W x x A Ma trận † 1 1 1 [ ] T T A W A AW A W o Giải tìm được x W A x λ 1 T o được gọi là ma trận tựa nghịch đảo có trọng số của ma trận A. Nếu W là ma trận đơn vị ta có là ma trận tựa nghịch đảo (phải) của A. Thay trở lại phương trình ban đầu b Ax AW A Ax λ 1 T † 1 [ ] T T A A AA† † x A b E A A x ( ) ; o o Do A có hạng đầy đủ nên Nếu A vuông, ko suy biến A A AA A A A [ ] [ ] T T T T khi x 0 x A b † o 1 ( ) det[ ] 0 T rank m A AW A † 1 1 A A A A λ 1 1 [ ] ( ) T AW A b Axo T T 1 1 10 Ma trận tựa nghịch đảo [pseudo-inverse] Hệ PT đại số tuyến tính với số PT nhiều hơn số ẩn. Thông thường PT không có nghiệm thỏa mãn. Ax b A x b , , , , m n n mm n 9/16/2015 Ma trận tựa nghịch đảo [pseudo-inverse] Đạo hàm J theo biến x và cho bằng 0, ta nhận được 1 1 rank n ( ) A ( ) 0 J T A W Ax b x Khi W là một ma trận đơn vị [ ] T T x A WA A Wb 1 PT không có nghiệm thỏa mãn, Tuy nhiên, ta có thể tìm được vectơ x sao J e and e e x A A A b A b n [ ] cho sai số 1 1 2 1 † T T T e Ax b min 2 2 i 1 i Một cách tổng quát ta có thể đặt lại bài toán như sau: tìm véctơ x sao cho hàm sau đây đạt cực tiểu 2 2 ( ) ( ) T T J e We Ax b W Ax b 1 1 W là ma trận trọng số đối xứng xác định dương † 1 [ ] T T A A A Ađược gọi là ma trận tựa nghịch đảo trái của A (xác định nghiệm hệ có số ẩn < số pt) † 1 [ ] T T A A AAđược gọi là ma trận tựa nghịch đảo phải của A (xác định nghiệm hệ có số pt < số ẩn) Ma trận tựa nghịch đảo [pseudo-inverse] Trường hợp ma trận A ở điều kiện yếu [ ill-condition ] Trong các trường hợp det(A) khác không nhưng nhỏ hơn chuẩn của nó nhiều lần, thì nghiệm tìm được rất nhạy đối với các phương pháp tính, đặc biệt một sự thay đổi rất nhỏ của A có thể làm cho nghiệm thay đổi rất nhiều, (A là ma trận ở điều kiện yếu) n n n Ax b A x b Ma trận tựa nghịch đảo [pseudo-inverse] Trường hợp ma trận A ở điều kiện yếu [ ill-condition ] f x Ax b Ax b x E x ( ) ( ) ( ) ( ) T T x A A E x b Ax x A b b b ( ) min T T T T T T det( ) A A , , , Đạo hàm theo biến véctơ x và cho kết quả bằng không, ta nhận được 1 ( ) 2( ) 2 0 T Khi giải các hệ ở điều kiện yếu, người ta hay sử dụng phương pháp bình phương tối thiểu có trọng số (damped least square inverse). Theo phương pháp này, nghiệm của phương trình được tìm với điều kiện bổ sung như f x T T A A E x A b x sau: 2 2 min( ) xAx b x f x Ax b Ax b x E x ( ) ( ) ( ) ( ) Chọn trọng số α > 0 Giải được 1 ( ) T T x A A E A b † 1 0 : [ ] T T x A b A A A b T T x A A E x b Ax x A b b b ( ) min T T T T T T Bài toán dạng này thường xuất hiện trong các bài toán động học ngược rôbốt công nghiệp khi mà rôbốt chuyển động gần hoặc qua các cấu hình kỳ dị, tại đó ma trận Jacôbi suy biến hay bị giảm hạng. 11 9/16/2015 Application of pseudo-inverse: Inverse kinematics Application of pseudo-inverse: Inverse kinematics Problems E E x r O tdetermined Given() n q t ∈ℜ ( ) [ , ] = ∈ m Method based on the pseudo-inverse of jacobian matrix ⎡ ⎤ ∂ ∂ ∂ ⎢ ⎥ f f f 1 1 1 ... ⎢ ⎥ ∂ ∂ ∂ ∂ = = ⎢ ⎥ q q q ∂= = f q q x J q q x f J qq 1 2 n f x q ( , ) 0 = nonlinear equations ( ) ( ) dor q ( ) ... ... ... ... ∂ ⎢ ⎥ more unknowns than equations ∂ dt ⎢ ⎥ ∂ ∂ ∂ f f f m m m ... ⎢ ⎥ ∂ ∂ ∂ ⎣ ⎦ q q q Methods: analytical and numerical Number of solutions system of m linear equations, with n unknown introduction the functional 1 2 n 12 min T J = → q Wq Solution W : positive weighting matrix q q q t dt Finite number of solutions Infinite number of solutions 1 1 1 ( )[ ( ) ( )] − − − T T q W J q J qW J q x = if W is identical matrix 1 ( )[ ( ) ( )] T T − + q J q J q J q x J x = = ( ) (0)t = + ∫0 1 ( )[ ( ) ( )] + − T T J J q J q J q = Nguyen Q. Hoang, Nguyen V. Khang 45 Application of pseudo-inverse: Inverse kinematics Blockdiagram of Inverse kinematics Nguyen Q. Hoang, Nguyen V. Khang 46 Application of pseudo-inverse: Inverse kinematics Method of Error feedback x† J q W ( ) 1/s q e Ke e x f q , ( ) e x J q q Ke ( ) e( ) 0 t Solving for q0 x(0) Drawback of this diagram : q(t) q(0) x q ( ), ( ) t t J q q x Ke ( ) q J x Ke ( ) Null space of jacobian matrixJq() o q J x Ke E J J z ( ) ( ) accumulated errors due to rounding and integral method makes don’t satisfy the constraint equations, so the end-effector moves out the given trajectory. f x q ( ( ), ( )) 0 t t ≠ Nguyen Q. Hoang, Nguyen V. Khang 47 By putting o zthe redundancy of manipulator can be exploited: - advoiding obstacles - advoiding impact with the joint limitations - advoiding singular configurations Nguyen Q. Hoang, Nguyen V. Khang 48 12 9/16/2015 Application of pseudo-inverse: Inverse kinematics Method of adjustment in error of joint variables Application of pseudo-inverse: Numerical experiments A 5-DOF planar manipulator Manipulator parameter x( )t ( ) x ( )tt q J x Ke ( ) E J J z * qx q y A q2 q3 BC C3 C4 q4 φ Link 1 2 3 4 5 m [kg] 10 7.5 5 4 1 Iz [kgm2] 0.11 0.10 0.03 0.02 0.02 l [m] 0.55 0.50 0.45 0.40 0.20 a [m] 0.15 0.16 0.19 0.20 0.10 x1s Adjusting ( ) o (0) x x o solving for qo q Determ. q q q1 O C2 C1 D g E x q5 2 1.5 Main idea of this method: find q closed to q*, such that q satisfy constraint equation f(x,q)=0 In this simulation, the end-effector (link 5) will be forced to move at the velocity of 2 m/s along a cirlular trajectory, while its orientation is constant, 1.0 rad. The trajectory has a center at point (0.8, 0.5) m and radius of R = 0.5 m. 1 0.5 0 0 2 4 6 8 10 t [s] x1 [m] x2 [m] x3 [rad] Nguyen Q. Hoang, Nguyen V. Khang 49 Application of pseudo-inverse: Numerical experiments A. Without consideration of joint limitations (planar manipulator) Nguyen Q. Hoang, Nguyen V. Khang 50 Application of pseudo-inverse: Numerical experiments / without post-adjusting q [rad] 3 2 1 1.5x 10-15 1 0.5 u [Nm] 30 20 10 0 q J x q J x Ke ( ) 0 -1 -2 0 -0.5 -1 0 2 4 6 8 10 -10 -20 -30 0 2 4 6 8 10 e1 [m] e2 [m] e3 [rad]0 2 4 6 8 10 t [s] q1 q2 q3 q4 q5 t [s] t [s] u1 u2 u3 u4 u5 B. Consideration of joint limitations (planar manipulator) q [rad] 2 1 0 1.5x 10-15 1 0.5 0 u [Nm] 30 20 10 0 [ ( ) ( )] q J x E J q J q z0 ( ) [ ( ) ( )] q J x Ke E J q J q z0 -1 -2 0 2 4 6 8 10 t [s] q1 q2 q3 q4 q5 -0.5 -1 0 2 4 6 8 10 t [s] e1 [m] e2 [m] e3 [rad] -10 -20 -30 0 2 4 6 8 10 t [s] u1 u2 u3 u4 u5 Nguyen Q. Hoang, Nguyen V. Khang 51 Nguyen Q. Hoang, Nguyen V. Khang 52 13 Application of pseudo-inverse: Numerical experiments / with post-adjusting q J x Ke ( ) q J x 9/16/2015 Application of pseudo-inverse:Hệ thống chân vịt đẩy Hệ thống chân vịt đóng một vai trò quan trọng trong việc điều khiển tàu lặn, đặc biệt trong các trường hợp yêu cầu độ chính xác định vị cao. Do đó cần thiết phải xây dựng mô hình chính xác nhất có thể cho hệ thống này. Trong phần này trình bày việc đưa ra mối quan hệ giữa lực điều khiển và các lực đẩy của chân vịt và mô hình động lực của hệ thống lực đẩy chân vịt. Việc điều khiển chuyển động tàu lặn phụ thuộc vào kiểu dáng của từng loại tàu và cách bố trí các động cơ – chân vịt đẩy. Phần này trình bày việc lựa chọn số lượng và cách bố trí hệ thống đẩy cho tàu lặn cỡ nhỏ dạng ROV. Các ưu nhược điểm của các phương pháp bố trí [ ( ) ( )] q J x E J q J q z0 ( ) [ ( ) ( )] q J x Ke E J q J q z0 Nguyen Q. Hoang, Nguyen V. Khang 53 được phân tích. Phần này cũng đề xuất một phương án phân bổ tối ưu lực điều khiển cho các động cơ. Ngoài ra, sự phân bố lại được thực hiện khi một vài động cơ ngừng hoạt động cũng được chú ý nhằm đảm bảo để tàu làm việc liên tục. Application of pseudo-inverse: Hệ thống chân vịt đẩy Application of pseudo-inverse: Hệ thống chân vịt đẩy u4 Quan hệ lực/mômen điều khiển và lực u1 Một số trường hợp bố trí hệ thống chân vịt đẩy chân vịt u2 T u u u 1 2 u [ , ,..., ] u3 O z y b u6 y b 1 1 u4 n Bu u8 τp φ 4 a a 4 5 u5 c [ , , , , , ]T F F F m m m x y z x y z τ y θ x a x a 3 x a b y uu6 3 2 u2 u1 5 ψ 2 u3 Để điều khiển được tất cả 6 dof của ROV x z np dim( ) dim( ) 6 & ( ) 6 u B rank a) 4 chân vịt, đk 4 dof b) 5 chân vịt, đk 4 dof c) 6 chân vịt, đk 6 dof Yêu cầu bố trí các chân vịt sao cho ma trận B có hạng đẩy đủ (full rank) 1 5 y y4 1 4 y8 n n p o tàu lặn thiếu hụt dẫn động 5 x 5 x 6 1 α α α α 4 x n n p o n n p o tàu lặn có dẫn động đầy đủ, và tàu lặn có hệ dẫn động dư (có tính chất dự phòng) 2 3 2 3 2 3 6 7 d) 5 chân vịt, đk 4 dof e) 6 chân vịt, đk 5 dof f) 8 chân vịt, đk 6 dof 14 Hệ thống chân vịt đẩy 9/16/2015 Application of pseudo-inverse: Hệ thống chân vịt đẩy z z Quan hệ lực điều khiển – lực đẩy Ma trận trạng thái động cơ c b0 0 0 0 c c c c τ = Bu { , , ..., } Ψ = diag a1 a2 an 0 0 0 0 s s s s Tiêu chuẩn tối ưu ai =1 - động cơ i không hỏng 1 1 ai = 0 - động cơ i hỏng (lỗi) 1 1 1 1 0 0 0 0 T J min = u u = u 2 || || y x B c c c c b b b b 0000 0000 2 Nghiệm tối ưu 2 2 Phân bổ lại lực điều khiển τ BΨ u B u 0 0 0 0 as as as as [ ]− = =T T # 1 = = ψ u B τ B BB τ [ ]− = T T 1 223 7 6 a 3 2 y 0.250[ ], 0.473[ ], 20o b c m a m Sơ đồ phân bổ lực Lực điều u ΨB BΨ B τ ψτ 4 1 8 5 x u B= khiển #Tàu lặn Phát hiện lỗi Application of pseudo-inverse: Hệ thống chân vịt đẩy Ví dụ B cos cos cos cos sin sin sin sin a a a a / 2 / 2 / 2 / 2 T τ [ , , ] X Y M 1 2 3 4 [ , , , ]T u u u u u Bốn động cơ-chân vịt làm việc bình thường a aX sin cos sin2 1 sin cos sin2 u B a aY 2 sin2 sin cos sin2 a a aM a a sin cos sin2 Khi động cơ số 4 lỗi, ta có ma trận phát hiện lỗi như sau Ψ diag([1,1,1,0]) 0 cos sin2 a 1 sin cos 0 a a B u B τ Y M 2sin a X Y 2cos 2sin sin2 sin 0 sin2 a a u X M 0 0 0 a 2cos0 15 9/16/2015 Lecture 3 Lập trình trong Matlab Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology 1 Các kiểu dữ liệu Chương 3. Lập trình trong Matlab 3.1 Các kiểu dữ liệu Kiểu véctơ và ma trận Chuỗi ký tự (ký tự, xâu - string) Kiểu ô, mảng (Cell-Array) Kiểu cấu trúc 3.2 Soạn thảo Script file trong Matlab Scripts Các hàm, MATLAB - Function 3.3 Các vòng lặp và rẽ nhánh Vòng lặp FOR Vòng lặp WHILE Lệnh if, cấu trúc if - else - end Cấu trúc switch-case 3.4 Các Mat-File 2 Các kiểu dữ liệu • Trong Matlab có nhiều kiểu dữ liệu, sau đây chỉ trình bày một số kiểu liên quan đến dạng thể hiện và lưu trữ. Để có thể nhận được sự trợ giúp trực tiếp từ Matlab, bạn cần tận dụng câu lệnh help. – Kiểu véctơ và ma trận – Chuỗi ký tự (ký tự, xâu - string) – Kiểu ô, mảng (Cell-Array) • Kiểu ô, mảng (Cell-Array) – Các dữ liệu của các kiểu khác nhau, ví dụ chuỗi ký tự, các ma trận có cỡ khác nhau, có thể được sử dụng như các ô của một mảng. Các phần tử ô được gọi đến thông qua chỉ số của nó. Các phần tử của một ô được đặt trong dấu ngoặc nhọn {…}. Ví dụ – Kiểu cấu trúc • Kiểu véctơ và ma trận Biểu diễn số các ma trận hai hay nhiều chiều cũng như các véctor trong chương 2 là một kiểu dữ liệu đặc biệt. Mỗi phần tử của véctơ hay ma trận theo chuẩn cần ô nhớ 8 Byte (class double) hoặc 4 Byte (class single). Các đại lượng phức cần một ô nhớ gấp đôi, phần thực và ảo được ghi riêng rẽ. • Chuỗi ký tự (ký tự, xâu - string) Chuỗi ký tự được đặt trong cặp dấu nháy đơn trên >> ’Day la chuoi ky tu’ Mỗi phần tử của chuỗi ký tự (Ma trận ký tự) cần ô nhớ 2 Byte. 3 >> A(1,1)={[1 2 3; 4 2 6; 1 7 9]}; >> A(1,2)={'Ma tran thu'}; >> A(2,1)={3+8i}; >> A(2,2)={0:pi/20:2*pi}; % goi den cac phan tu >> A(1,1) ans = [3x3 double] >> A{1,1} % o (1,1) ans = 1 2 3 4 2 6 1 7 9 >> A{1,2} % o (1,2) ans = Ma tran thu >> A(1,2) ans = 'Ma tran thu' 4 1 9/16/2015 Các kiểu dữ liệu • Kiểu cấu trúc Các kiểu dữ liệu • Kiểu cấu trúc (ví dụ) – Cấu trúc được sử dụng để làm việc với các kiểu dữ liệu khác nhau. Tên của một cấu trúc gồm hai phần, tên cấu trúc trước dấu chấm ( . ) và tên trường trong cấu trúc sau dấu chấm, (struct_name.field_name). Các phần tử của cấu trúc được gọi đến qua tên và chỉ số. – Cú pháp % tao mot cau truc structur=struct(’name_1’, value1, ’name_2’, value2,.. ) % tro den phan tu cua cau truc structur.name hoặc structur.name_1 = value_1; structur.name_2 = value_2; 5 >> A=[1 2 3; 4 2 6; 1 7 9]; >> my_structur = struct('data', A, 'dimension',[3 3]) my_structur = data: [3x3 double] dimension: [3 3] >> my_structur.dimension ans = 3 3 >> my_structur.data ans = 1 2 3 4 2 6 1 7 9 >> my_structur(2).data=inv(A) % mo rong truong data my_structur = 1x2 struct array with fields: data dimension >> my_structur.data % in ra ma tran A va inv(A) ans = 1 2 3 4 2 6 1 7 9 ans = 4.0000 -0.5000 -1.0000 5.0000 -1.0000 -1.0000 -4.3333 0.8333 1.0000 6 Soạn thảo Script file trong Matlab • Matlab Script files là một chuỗi các lệnh được gõ vào trong cửa sổ soạn thảo và được ghi vào tệp có phần mở rộng là m (filename.m, hay được gọi là m-file). Để tạo m-file, trong menu chính bạn chọn New ==> M-file hoặc bạn nháy chuột vào biểu tượng “tờ giấy trắng” trên thanh công cụ (hình 3-1). • Việc chạy m-file tương đương với việc đánh toàn bộ các dòng lệnh trên cửa sổ lệnh tại dấu nhắc ‘>>’ của Matlab. Các biến sử dụng trong m-file được đặt vào trong không gian làm việc của Matlab. Không gian làm việc này là trống khi khởi động Matlab, và nó sẽ chứa tất cả các biến được định nghĩa trong phiên làm việc. Muốn xóa tất cả các biến ta sử dụng lệnh clear all 7 Soạn thảo Script file trong Matlab • Ví dụ soạn m-file vẽ đường tròn là hàm của bán kính r.8 2 9/16/2015 Soạn thảo Script file trong Matlab Ví dụ về chương trình không có biến vào function V = volume(r) %r = input('Hay cho ban kinh cua hinh cau') vol = (4/3)*pi*r^3; disp('The tich hinh cau nay la:') disp(vol) Save với tên file volume.m, và gọi hàm trong command windows >> volume 9 Soạn thảo Script file trong Matlab Ví dụ về script file r1 = input('Hay cho ban kinh cua hinh cau') r2 = input('Hay cho ban kinh cua hinh cau') disp(‘Tong the tich hai hinh cau’) Tong_volume(r1, r2) Save với tên file TinhV.m, và gọi hàm trong command windows >> tinhV Soạn thảo Script file trong Matlab Ví dụ về chương trình không có biến vào function vol = volumef(r) vol = (4/3)*pi*r^3; Save với tên file volumef.m, và gọi hàm trong command windows >> volumef(2) function V = Tong_volume(r1, r2) V1 = volumef(r1); V2 = volumef(r2); V = V1 + V2;>> Tong_volume(2, 5) 10 Soạn thảo Script file trong Matlab • Các m-file dạng hàm function [output 1, output 2] = function_name(input1,input2) % Some comments that explain what the function does go here. global x y z persistent a b c MATLAB command 1; MATLAB command 2; MATLAB command 3; output 1 = ….; output 2 = ….; 11 12 3 9/16/2015 Soạn thảo Script file trong Matlab • Biến toàn cục (global variable) và biến cục bộ hay biến địa phương (local variable). Biến toàn cục được khai báo theo cú pháp global var_1 var_2 • Các biến toàn cục được xóa khỏi môi trường tính toán bằng lệnh clear [tên biến]. Biến cục bộ chỉ được sử dụng trong một hàm, biến này không được liệt kê trong cửa sổ workspace, và do đó không thể tác động đến được trong môi trường làm việc. • Biến persistent trong một hàm có thể được thống nhất với cách khai báo persistent var_1 var_2 • Trái với các đại lượng đã được khai báo bằng global, các đại lượng với khai báo persistent chỉ được biết đến trong hàm mà nó được khai báo. Do đó các hàm khác không thể tác động được đến biến này. Các biến persistent chỉ được xóa khi hàm chứa nó thoát ra khỏi bộ nhớ (clear function_name) hoặc khi hàm được thay đổi sau đó được ghi lại. 13 Các vòng lặp và rẽ nhánh Ví dụ về lệch for. Viết một hàm sử dụng vòng lặp for để tính tổng các Các vòng lặp và rẽ nhánh Lệnh for for var_= bien_chay các lệnh end Lệnh while while dieukien, các lệnh end Lệnh if (nếu .. thì ..) if dieukien ... ……. ……. end Lệnh switch (chuyển) switch .. case .. end Lệnh ngắt vòng lặp trong for và while break; 14 Các vòng lặp và rẽ nhánh Ví dụ về lệch while. Viết hàm tính tổng S = 1+1/2+1/3+ ... +1/100 phần tử của một véctơ v. Trên thanh công cụ, ta lựa chọn New –> m-file để mở một cửa sổ soạn thảo, trong cửa sổ này ta soạn nội dung sau function tong = tongvector(v) % ham tinh tong cac phan tu cua mot vector n=length(v); % tra lai so phan tu, chieu dai cua vector s = 0; % khoi gan cho tong for i = 1:n % hoac for i = 1 : 1 : n s = s + v(i); end; tong = s; % save with file name tongvector.m >> x=[1:2:50]; % su dung vong lap while function tongS n=100; S=0; i=1; while i<=n S = S+1/i; i = i+1; end tong=S; disp('Voi n = 100, tong la :') disp(tong) >> >> tong_S % su dung vong lap for function tong_Sn n=100; S=0; for i=1:n S=S+1/i; end tong=S; disp('Voi n = 100, tong la :') disp(tong) >> >> tong_S Sử dụng hàm >> s=tongvector(x) s = 625 15 Voi n = 100, tong la : 5.1874 Voi n = 100, tong la : 5.1874 16 4 9/16/2015 Các vòng lặp và rẽ nhánh Sử dụng lệng while giải bài toán. Ví dụ tìm số n lớn nhất sao cho tổng Các vòng lặp và rẽ nhánh Lệnh if, cấu trúc if - else - end sau nhỏ hơn 27 n 1 1 1 1 1 ... 27 = + + + + = < ∑ Lệnh điều kiện nếu - thì. Cú pháp IF expression Sn i = statements 1 2 3 i 1 ELSEIF expression function n = timn(S) % tim so n lon nhat thoa man % tongS(n) = 1+1/2+1/3+...+1/n < S % gia tri S > 1 vao tu ban phim T=1; while T disp('Vao so duong S') S = input('S = ') if (S>1 & S<10) T = 0; end end format long format long tong=0; i=0; while(tong> A A = 4 0 1 0 0 0 4 0 1 0 1 0 4 0 1 0 1 0 4 0 0 0 1 0 4 19 Các vòng lặp và rẽ nhánh Cấu trúc switch-case Lệnh chuyển trong nhiều trường hợp dựa trên biểu thức và có dạng tổng quát như sau SWITCH switch_expr CASE case_expr, statement, ..., statement CASE {case_expr1, case_expr2, case_expr3,...} statement, ..., statement ... OTHERWISE, statement, ..., statement END 20 5 9/16/2015 Các vòng lặp và rẽ nhánh diem = 1; % 2, 3, 4, 5, 6 switch diem case 1 disp('Diem gioi') case 2 disp('Diem tot') case 3 disp('Diem kha') case 4 disp('Diem trung binh') otherwise disp('Diem kem.') end 21 Bài tập Các Mat-File Mat-file là tệp có phần mở rộng là mat và do đó được gọi là tệp chấm mat (filename.mat). Mat-file là những tệp nhị phân được nén để lưu trữ kết quả số. Các file này được sử dụng để ghi các kết quả mà đã được tạo ra bởi chuỗi các lệnh, chỉ thị của Matlab. Chẳng hạn, để ghi giá trị của hai biến x1 và x2 trong file có tên dulieu.mat ta đánh dòng lệnh >> x1=5; >> x2=10; >> save dulieu.mat x1 x2 Muốn tải (đọc) một mat-file vào MATLAB ta gõ dòng lệnh >> load filename (or load filename.mat) 22 1. Viết một hàm Matlab hỏi người sử dụng bán kính và chiều cao của hình trụ, sau đó tính toán và đưa ra màn hình diện tích toàn phần của hình trụ và thể tích của hình trụ. 2. Viết một chương trình con có sử dụng lệnh while để tính tổng S(x,n) ( , ) 1 ... n S x n x x x = + + + + 2 Viết thuật toán (thuật giải) f x x a a x a ( ) 0, 0, = − = > = 2 k x kf x = = 1;, ( ) 1; ( ) x xf x 3. Viết một chương trình con có sử dụng lệnh while để tính căn của a>0 theo = −′ k k k + 1 pương pháp lặp Newton-Raphson, theo công thức [heron]? Với x(1) = 1, độ chính xác yêu cầu 2 x a x a ( ) k 2 2 x a x a x a − + = − = = + k k k + n n xx x là 0.0000001 k xx x x 2 2 2 2 n = = + + 12 2 2 x x− k k k n n 7 | | 10 n n +− < 1 1 ⎛ ⎞ a 4. Hãy sử dụng vòng lặp while để tìm thương và số dư khi chia hai số nguyên a = + ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ k xx 2 cho b, a > b. k 23 24 6 9/16/2015 Viết thuật toán (thuật giải) 1. Nhập a > 0; 2. Khởi động (khởi gán), k = 1, x(1) =1; ss = 1; eps = 1e-5; 3. While ss > eps x(k+1) = x(k)/2 + a/2/x(k); ss = abs(x(k+1) – x(k)); k = k + 1; end 4. x = x(k - 1) 25 7 9/16/2015 Lecture 4 Đồ họa trong Matlab Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology Nguyen Quang Hoang 1 Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D • Trước hết ta bắt đầu với việc vẽ đồ thị hàm một biến y = f(x). Công việc này trong Matlab bao gồm ba bước: – Định nghĩa hàm cần vẽ, y = f(x) – Xác định miền giá trị của biến, x thuộc miền [a,…,b] – Gọi hàm plot(x,y) của Matlab. Chương 4. Đồ họa trong Matlab 4.1 Đồ thị trong mặt phẳng - 2D Đặt màu và kiểu đường cho đồ thị Một số tùy chọn khi vẽ đồ thị 2D Vẽ nhiều đồ thị trên cùng một hệ trục Các lệnh axis Đặt giới hạn miền vẽ với lệnh axis Lệnh Subplots Vẽ các đồ thị xếp chồng và lệnh linspace Vẽ biểu đồ với lệnh contour – vẽ đường đồng mức Thêm chú thích trên đồ thị Vẽ đồ thị các hàm có điểm không xác định 4.2 Các lệnh vẽ trong không gian – 3D Nguyen Quang Hoang 2 Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D • Đặt màu và kiểu đường cho đồ thị – Khi muốn thể hiện nhiều đồ thị trong một hình người ta thường phân biệt các đồ thị bằng các màu, và kiểu đường khác nhau (nét liền, nét đứt) – Trong Matlab đã định nghĩa sẵn các lựa chọn • Ví dụ cần vẽ đồ thị hàm y = sin(x) trong khoảng từ 0 đến 10. – Chia khoảng cần vẽ [0,10] bằng các điểm chia cách đều, sử dụng toán tử (:) , chẳng hạn x = [a:h:b], với h là bước chia. – Tính giá trị của hàm f(x) tại các điểm chia tương ứng, y(i) = f(x(i)); kết quả được hai véc tơ x và y có cùng số các phần tử. – Ví dụ: • >> x=[0:0.01:10]; y1=sin(x); • >> y2=cos(x); plot(x,y1, x,y2,'--') 1 0.5 0 – Gọi hàm plot(x,y). Cụ thể: >> x = [0:0.1:10]; >> y = sin(x); >> plot(x,y) 1 0.5 0 -0.5 -1 0 2 4 6 8 10 -0.5 -1 0 2 4 6 8 10 Hãy thực hiện với h = 0.01; h = 0.5; h = 1. Nhận xét đồ thị. Nguyen Quang Hoang 3 Department of Applied Mechanics Nguyen Quang Hoang 4 Department of Applied Mechanics 1 9/16/2015 Đồ thị trong mặt phẳng – 2D Đồ thị trong mặt phẳng – 2D Kiểu đường - solid line nét liền -- dashed line nét đứt -. dashdot line nét chấm gạch : dotted line nét chấm (none) no line (không hiện) Kiểu đánh dấu (marker) . point điểm o circle vòng tròn x x-mark chữ x + plus dấu cộng * star dấu sao s square dấu vuông ….. …… Ký hiệu k Màu black – đen b blue – xanh da trời c cyan – lục lam g green – xanh lá cây r red – đỏ m magenta– đỏ tươi y yelow – vàng w white – trắng • lệnh legend cho phép ta điền các chú giải lên hình vẽ x [m] 3 • Ví dụ: 2 >> t = [0:0.1:10]; 1 0 >> x = sin(2*t); v = 2*cos(2*t); -1 >> plot(t,x,'k-', t,v,'k--'), -2 >> grid on, xlabel('t [s]') v [m/s] Sử dụng lệnh trợ giúp: >> help plot Nguyen Quang Hoang 5 >> legend('x [m]','v [m/s]') >> axis([0 10 -2.6 3.5]) 0 1 2 3 4 5 6 7 8 9 10 t [s] Nguyen Quang Hoang 6 Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D • Lệnh subplot - cho phép ta vẽ nhiều đồ thị trong nhiều hệ trục tọa độ vào cùng một hình vẽ. Lệnh này được gọi với cú pháp subplot(m, n, p), trong đó m và n cho biết số hàng và số cột trong hình vẽ như là một ma trận hay mảng. Số p cho biết đồ thị sẽ được vẽ vào ô thứ mấy trong mảng. Thứ tự các ô trước hết theo chỉ số hàng. Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D • Ví dụ: t = [0:0.01:8]; x = sin(3*t); v=3*cos(3*t); a=-9*sin(3*t); subplot(3,1,1), plot(t,x,'k-','Linewidth',1), grid on, ylabel('x [m]'), subplot(3,1,2), plot(t,v,'k-','Linewidth',1), grid on, ylabel('v [m/s]'), subplot(3,1,3), plot(t,a,'k-','Linewidth',2), grid on, ylabel('a [m/s^2]'), subplot(m,n,p) % p = 1,2,...mxn >> for i = 1:6 subplot(2,3,i), end Fig. 1 Fig. 2 Fig. 3 xlabel('t [s]‘) plot 1 plot 2 1 x [m] 0 -1 0 1 2 3 4 5 6 7 8 2 v [m/s] 0 -2 0 1 2 3 4 5 6 7 8 ] 2 10 Fig. 4 Fig. 5 Fig. 6 plot 3 a [m/s 0 Nguyen Quang Hoang 7 Department of Applied Mechanics -10 0 1 2 3 4 5 6 7 8 t [s] Nguyen Quang Hoang 8 Department of Applied Mechanics 2 9/16/2015 Vẽ hàm y=y(x), biết x = x(t), y=y(t) x a t Đồ thị trong mặt phẳng – 2D • Đồ thị theo tọa độ cực và Đồ thị với thang chia Logarith 0.25 0.2 sin( ) 90 Bien do 8 0.15 120 60 0.1 6 0.05 4 150 30 y a t y y x cos( ) ( ) 180 010-1 100 101 102 0 2 0 -1 -2 210 240 270 330 300 Goc pha -3 -4 10-1 100 101 102 log(ω) Lệnh polar(phi,r) a = 1; phi=[0:pi/90:2*pi]; r = a*phi; polar(phi, r, '-k'); Nguyen Quang Hoang 9 Department of Applied Mechanics Ngoài ra ta có thể vẽ các đồ thị với lệnh semilogx hoặc semilogy, với các lệnh này chỉ có một trục x hoặc y tương ứng sử dụng thang logarith còn trục kia sử dụng thang thập phân Nguyen Quang Hoang 10 Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D Các lệnh plot 2D plot(x,y) Vẽ đồ thị bằng cách nối thẳng các điểm [x(i),y(i)] comet(x,y) Vẽ hoạt hình một quĩ đạo stair(x), stairs(x,y,..) Vẽ đồ thị kiểu bậc thang stem(x,y) spy(matrix) hiển thị một ma trận thưa semilogx(x,y) Hiển thị y(x), sử dụng thang logarith cho trục x semilogy(x,y) Hiển thị y(x), sử dụng thang logarith cho trục y loglog(x,y) Hiển thị y(x), sử dụng thang logarith cho trục x và y fplot(ham,mien) Vẽ hàm số trong miền hold giữ đối tượng sẵn có, đóng lại Nguyen Quang Hoang 11 Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D Để vẽ nhiều đồ thị trong một cửa sổ đồ họa, ta sử dụng lệnh hold on ngay sau lệnh vẽ đồ thị đầu tiên. Lựa chọn này được dỡ bỏ nếu sau lệnh vẽ plot cuối cùng ta sử dụng lệnh hold off. Các giải thích trên về kiểu đường, màu sắc, đánh dấu (line style, colour, marker) được làm sáng tỏ trong thí dụ sau. figure(1) clf % xoa tat ca cac do thi hien co t = [0:pi/40:2*pi]; plot(t,sin(t),'-.r*'), % ve do thi 1 hold on plot(t,sin(t-pi/2),'mo') % ve do thi 2 plot(t,sin(t-pi),':bs') % ve do thi 3 axis([0, 2*pi, -1.2, 1.2]) hold off grid on Nguyen Quang Hoang 12 Department of Applied Mechanics 3 9/16/2015 Đồ thị trong mặt phẳng – 2D Khi vẽ đồ thị hàm phức tập f(x), ta có thể sử dụng các phép tính phần tử: .*, ./ , .^ Ví dụ: 0.40 2 sin 2 x y e x x Đồ thị trong mặt phẳng – 2D Các trục tọa độ và các chú thích trên đồ thị axis([xmin, xmax, ymin, ymax]) Điều chỉnh phạm vi các trục tọa độ axis([-inf, xmax, ymin, ymax]) Tự động chọn xmin grid Sử dụng chia lưới hay không chia lưới axis Năm lựa chọn đối với các trục xlabel(String) ylabel(String) Chú thích cho trục x Chú thích cho trục y title(String) text(String) legend(string_1, string_2,… <,position>) Hiển thị tên đồ thị Chèn chữ vào đồ thị Hiển thỉ các chú giải clc, clear all x=[0:0.05:12]; y=exp(-0.4*x).*x.^2+sin(2*x); plot(x,y,'k-','linewidth',2), grid on xlabel('x'), ylabel('y') 5 4 3 y 2 1 0 0 2 4 6 8 10 12 x Nguyen Quang Hoang 13 Department of Applied Mechanics Đồ thị trong mặt phẳng – 2D Lưu ý khi vẽ đồ thị của các hàm f(x) có điểm không xác định. Ví dụ 2 4(1 ) ( )(10 7 )(1 ) 2 − x =− − − 1 2 x x 0.7990, 1.3380 f xx x 2 2 Nguyen Quang Hoang 14 Department of Applied Mechanics Đồ thị trong không gian – 3D Các lệnh vẽ trong không gian – 3D plot3(x,y,z <,plotstil>) comet3(x,y,z <,comet length>) mesh(x,y,z <,color>) surf(x,y,z <,color>) surfc(x,y,z <,color>) surfl(x,y,z <,color>) patch(x,y <,z> ,color) waterfall(x,y,z..<,..>,..) contour3(x,y,z.. <,..>,..) contour(x,y,z.. <,..>,..) Ví dụ Cần thể hiện đồ thị hàm 2 2 z ye− − = x y trong miền Nguyen Quang Hoang 15 Department of Applied Mechanics − ≤ ≤ 2 , 2 x y Nguyen Quang Hoang 16 Department of Applied Mechanics 4 9/16/2015 Đồ thị trong không gian – 3D • [x,y] = meshgrid(-2:0.1:2); • z = y.*exp(-x.^2-y.^2); • Ví dụ: Đồ thị trong không gian – 3D 2 2 0.2 ( ) 2 2 ( , ) sin( ), 3 , 3 x y z f x y e x y x y − ⋅ + = = ⋅ + − ≤ ≤ • mesh(x,y,z), xlabel('x'), ylabel('y'), zlabel('z') mesh(x,y,z) surf(x,y,z) surfc(x,y,z) x = [-3 : 0.1 : 3]; y = [-3 : 0.1 : 3]; [X, Y] = meshgrid(x,y); f = exp(-0.2*(X.^2+Y.^2)).*sin(X.^2+Y.^2); mesh(X,Y,f); surf(X,Y,f) contour3(X,Y,f,30); Nguyen Quang Hoang 17 Department of Applied Mechanics Nguyen Quang Hoang 18 Department of Applied Mechanics 5 16/09/2015 Lecture 5 Nội suy và làm mịn đường cong Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology Nguyen Quang Hoang Department of Applied Mechanics Nêu bài toán • Trong kỹ thuật khi đo đạc ta nhận được các số liệu rời rạc, từ các số liệu đo này ta cần biểu diễn các đường để biết được đặc tính của hệ khảo sát và có thể suy ra được những giá trị lận cận các giá trị đã đo. Matlab đã cung cấp sẵn các hàm cho phép đưa ra được các đường biểu diễn phù hợp nhất với bộ số liệu đã có. Trong chương này ta sẽ trình bày các kỹ thuật đơn giản thực hiện công việc đó. • Mịn hóa bằng đa thức Lệnh polyfit(x,y,n) cho ta một đa thức bậc n dạng 1 2 1 1 ( ) ... n n p x p x p x p x p x p n n n n − = + + + + + − + 1 2 1 các hệ số của đa thức này được suy ra hai véctơ số liệu x và y, nhờ phương pháp sai số bình phương bé nhất (least mean square). 1 2 3 1 2 3 [ ... ], [ ... ] m m x y x x x x y y y y m n > Chương 5. Lập trình trong Matlab 5.1 Mịn hóa đường cong bằng đa thức 5.2 Nội suy đa thức 5.3 Bài tập thực hành 1 Nguyen Quang Hoang 2 Department of Applied Mechanics Phương pháp sai số bình phương bé nhất Trường hợp đơn giản nhất, n = 1, đa thức có dạng một phương trình đường thẳng y ax b Theo phương pháp sai số bình phương bé nhất hai hệ số a và b được tìm như sau. Sai số giữa giá trị của số liệu đo y(i) và giá trị tính theo đường thẳng vừa đưa ra được tính theo công thức ( ) , 1,2,..., i i i e ax b y i m = + − = Như thế tổng bình phương các sai số là m m = = + − ∑ ∑ 2 2 S a b e ax b y ( , ) ( ) i i i i i = = 1 1 Bài toán đặt ra ở đây là tìm a, b sao cho hàm S đạt cực tiểu. Nguyen Quang Hoang 3 Department of Applied Mechanics Nguyen Quang Hoang 4 Department of Applied Mechanics 1 16/09/2015 Phương pháp sai số bình phương bé nhất y ( ) i i n i e y P x = − i e ( , ) i i x y 1 n P x ax b = = + 1, ( ) y y ax b = + n y P x = ( ) x . Phương pháp sai số bình phương bé nhất Đạo hàm S theo các biến a và b ta nhận được m m ∂ ∂ S S ax b y x ax b y ∂ ∂ ∑ ∑ = + − = = + − = 2( ) 0, 2( ) 0 i i i i i a b = = i i 1 1 m m m m m = = = = = ∑ ∑ ∑ ∑ ∑ + − = + − = a x x b x x y a x mb y 0, 0 i i i i i i i i i i i i 1 1 1 1 1 Từ đây giải được m m m ∑ ∑ ∑ y x x x x y x y mx y − − i i i i i i i i i m m = = = 1 1 1 b a = = , ∑ ∑ 2 2 x x mx x x mx − − i i i i i i = = 1 1 i e m với các giá trị trung bình y x x ∑ ( ) ( , ) i i x y i i i e y ax b = − − = = ∑ ∑ 1 = i i − m m a b y xa = = − 1 1 x x y y , im x x x , i i ∑ ( ) x Nguyen Quang Hoang 5 Department of Applied Mechanics m m = = i i 1 1 Nguyen Quang Hoang Department of Applied Mechanics i i i = 1 − 6 Phương pháp sai số bình phương bé nhất Sử dụng ma trận tựa nghịch đảo (nghịch đảo trái) Quan hệ tuyến tính , 1,2,..., i i x a b y i m x y 1 1 1 x a y x yAp y 1 2 2 1 ... ... ... Phương pháp sai số bình phương bé nhất Sử dụng ma trận tựa nghịch đảo (nghịch đảo trái) Quan hệ bậc 2 1 2 3 , 1,2,..., i i i x p x p p y i m 2 x x y 2 1 1 1 1 x x y x x yAp y 1 p 2 1 ... ... ... ... 2 2 2 b 1 m m 1 p 2 p [ ] T T p A A A y Nguyen Quang Hoang 7 Department of Applied Mechanics 2 3 m m 1 m Nguyen Quang Hoang Department of Applied Mechanics [ ] T T p A A A y 8 2 16/09/2015 Phương pháp sai số bình phương bé nhất Ví dụ 1. Cho biết các số liệu đo được như sau x 6 8 10 12 14 16 18 20 22 24 y 3.94 3.8 4.1 3.87 4.45 4.33 4.12 4.43 4.6 4.5 Phương pháp sai số bình phương bé nhất Theo cách biểu diễn đa thức bằng một véctơ bắt đầu từ hệ số ứng với số hạng có bậc lũy thừa lớn nhất, ta có Nếu biểu diễn các số liệu trên bằng một đường thẳng, y = a x + b. Hãy tìm các giá trị của a, b và tìm giá trị của y khi x = 11. Trước hết ta cần nhập hai véctơ số liệu x và y vào Matlab. >> x = [6: 2 : 24] x = 6 8 10 12 14 16 18 20 22 24 >> y = [3.94 3.8 4.1 3.87 4.45 4.33 4.12 4.43 4.6 4.5] y = 3.94 3.80 4.10 3.87 4.45 4.33 4.12 4.43 4.60 4.50 Tiếp theo sử dụng lệnh polyfit để tìm các hệ số của đa thức >> p = polyfit(x,y,1) p = 0.0392 3.6267 >> a=p(1) a = 0.0392 >> b=p(2) b = 3.6267 Giá trị cần tìn tại x = 11 được tính bằng lệnh polyval(p,x) ans = 4.0574 >> polyval(p,11) 5 4 3 2 y 1 0 6 8 10 12 14 16 18 20 22 24 x Nguyen Quang Hoang 9 Department of Applied Mechanics Phương pháp sai số bình phương bé nhất Xét tổng bình phương của sai số là >> w = a*x+b; >> e = w-y; >> S1 = sum(e.*e) S1 = 0.2274 Nguyen Quang Hoang 10 Department of Applied Mechanics Mịn hóa bằng đa thức Ví dụ 2. Cho bộ số liệu sau x 0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 5.0 6.0 6.1 7.0 y 300 281 261 244 228 214 202 191 181 164 151 149 141 Nếu các điểm trên được làm mịn bằng đường bậc 2: >> n=2; >> p = polyfit(x,y,n); p = -0.0002 0.0443 3.5940 >> a=p(1); b=p(2); c = p(3); >> w = a*x.^2+b*x+c; >> e = w-y; >> S2= sum(e.*e) S2 = 0.2272 Ta thấy sai số có nhỏ hơn (S2> x = [0, 0.5, 1.0,1.5, 2, 2.5,3, 3.5,4, 5, 6, 6.1, 7]; >> y = [300,281,261,244, 228,214,202,191,181,164,151, 149,141]; >> plot(x,y,'ko', x,y, 'k-'), grid on, xlabel('x'), ylabel('y') Nguyen Quang Hoang Department of Applied Mechanics 11 Nguyen Quang Hoang Department of Applied Mechanics 12 3 16/09/2015 Mịn hóa bằng đa thức Mịn hóa bằng hàm e mũ >> n=2; % Xap xi bac 2 >> p = polyfit(x,y,n) >> a=p(1); b=p(2); >> c = p(3); >> w = a*x.^2+b*x+c; >> e = w-y; >> S2= sum(e.*e) S2 = 15.4079 300 250 200 y 150 100 0 1 2 3 4 5 6 7 x >> n=3; % Xap xi duong bac 3 >> p = polyfit(x,y,n) >> a=p(1); b=p(2); c = p(3); d = p(4); >> w = a*x.^3+b*x.^2+c*x+d; >> e = w-y; >> S3= sum(e.*e) S3 = 2.7933 350 300 250 y 200 150 100 0 1 2 3 4 5 6 7 x Trong nhiều trường hợp khi mịn hóa bằng đa thức không đáp ứng được yêu cầu đặt ra, ta cần phải tìm một phương án khác. Biểu diễn các số liệu thông qua hàm e mũ là một phương án nên được xem xét. Giả sử bộ số liệu (x,y) có thể lấy xấp xỉ bằng hàm ax y be = Bây giờ ta cần tìm hai hệ số a và b để các điểm dữ liệu nằm gần đường cong nhất. Nếu áp dụng trực tiếp phương pháp sai số bình phương bé nhất ta sẽ được một hệ phương trình phi tuyến đối với các ẩn a và b. Để tránh điều đó, ta lấy lôgarit cơ số tự nhiên hai vế được ln ln y ax b = + Đặt w = ln(y), z= x và p1 = a và p2 = ln(b) bài toán trở thành tìm các hệ số p1, p2 để xấp xỉ các số liệu (z, w) = (x, ln(y)) bằng đường thẳng w p z p = + 1 2 Nguyen Quang Hoang Department of Applied Mechanics Mịn hóa bằng hàm e mũ 13 Nguyen Quang Hoang Department of Applied Mechanics Mịn hóa bằng hàm e mũ 14 Ví dụ xét bộ số liệu sau Nhập số liệu vào Matlab x 1.2 2.8 4.3 5.4 6.8 7.9 y 7.5 16.1 38.9 67.0 146.6 266.2 Xấp xỉ bằng hàm e mũ >> p = polyfit(x,log(y),1) p = 0.5366 1.3321 300 250 200 >> x = [1.2 2.8 4.3 5.4 6.8 7.9]; >> y = [7.5 16.1 38.9 67.0 146.6 266.2]; 300 >> a = p(1); b=exp(p(2)); >> w = b*exp(a*x); y 150 >> plot(x,y,'ko', x,y, 'k-'), grid on, xlabel('x'), 250 200 >> err = w-y; >> S3= sum(err.*err) Se = 17.6259 100 50 ylabel('y‘) y 150 100 50 0 1 2 3 4 5 6 7 8 x >> x3 = [0:0.1:8]; >> y3 = b*exp(a*x3); >> plot(x,y,'ko', x3,y3, 'k-'), grid on, xlabel('x'),ylabel('y') 0 0 1 2 3 4 5 6 7 8 x Nguyen Quang Hoang Department of Applied Mechanics 15 Nguyen Quang Hoang Department of Applied Mechanics 16 4 16/09/2015 Nội suy bằng đa thức Lệnh polyfit(x,y,n) cho ta một đa thức bậc n mà đồ thị của nó là đường gần với đường tạo ra khi nối các điểm cho bởi bộ số liệu, còn nội suy đa thức cũng cho ta một đường cong nhưng đi qua các điểm số liệu đó. Có nhiều phương pháp để tạo ra đa thức đi qua các điểm cho như phương pháp nội suy theo công thức Lagrange, phương pháp Newton, phương pháp Neville, phương pháp nội suy với đường spline bậc 3, .... Chi tiết về các phương pháp này có thể xem trong các tài liệu về giải tích số và phương pháp số. Trong phần này chỉ trình bày việc sử dụng các lệnh interp của Matlab để thực hiện tìm các giá trị nội suy từ các số liệu sẵn có. Bài toán đặt ra là cho biết các giá trị tương ứng của hai đại lượng x và y như trong bảng số liệu x x1 x2 x3 … xn y y1 y2 y3 … yn Nội suy bằng đa thức Hãy tìm giá trị của biến y* ứng với giá trị nào đó của biến x* với x1<= x* <= xn. Bài toán này trong Matlab được giải quyết dễ dàng với lệnh interp. Cách gọi lệnh này như sau: y_star = interp1(x,y,x_star) Nếu x là véctơ gồm các phần tử nguyên từ 1 đến n (x = 1:n), với n = length(y), thì ta có thể gọi lệnh y_star = interp1(y,x_star) Để nói rõ hơn phương pháp nào được sử dụng trong việc nội suy này ta sử dụng cú pháp y_star = interp1(x, y, x_star, method) với method được chọn là một trong các phương pháp sau 'nearest' - nội suy theo lân cận gần nhất 'linear' - nội suy tuyến tính 'spline' - nội suy sử dụng đường spline bậc 3 'pchip' - nội suy bậc 3 từng đoạn 'cubic' - tương tự như 'pchip' Nguyen Quang Hoang Department of Applied Mechanics Nội suy bằng đa thức 17 Nguyen Quang Hoang Department of Applied Mechanics Nội suy bằng đa thức 18 Nếu giá trị x* nằm ngoài khoảng [x1, xn], khi đó ta gọi phép tính giá trị y* = y(x*) là phép tính ngoại suy và sử dụng lệnh với cú pháp: y_star = interp1(x, y, x_star, method, 'extrap') Sau đây xét ví dụ: >> x = 1:20; y = sin(x); >> x_star = 1.5; >> y_star = interp1(x,y, x_star) y_star = 0.8754 >> y_star = interp1(x,y, x_star, 'linear') y_star = 0.8754 >> y_star = interp1(x,y, x_star, 'cubic') y_star = 0.9008 >> y_star = interp1(x,y, x_star, 'nearest') y_star = 0.9093 >> y_star = interp1(x,y, x_star, 'spline') y_star = 1.0200 Như thế các phương pháp nội suy khác nhau sẽ cho ta các kết quả khác nhau. Để thấy rõ hơn kết quả của các phương pháp nội suy, ta xét ví dụ sau x = 0:8; y = sin(x); plot(x,y,'-ko'), hold on xi = 0:0.1:8; yi1 = interp1(x,y, xi, 'linear'); yi2 = interp1(x,y, xi, 'cubic'); yi3 = interp1(x,y, xi, 'spline'); plot(xi,yi1,'-k','linewidth',1.5) plot(xi,yi2,'-r','linewidth',1.5) plot(xi,yi3,'-b','linewidth',1.5), grid on xlabel('x'), ylabel('y'), axis([0 8, -1.1 1.1]) legend('points','linear','cubic','spline') Nguyen Quang Hoang Department of Applied Mechanics 19 Nguyen Quang Hoang Department of Applied Mechanics 20 5 16/09/2015 y Nội suy bằng đa thức 1 0.5 0 Bài tập 2. Tìm phương trình đường thẳng y = a x + b thể hiện quan hệ tuyến tính của bộ số liệu cho trong bảng x 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 y 3.076 2.810 2.588 2.297 1.981 1.912 1.653 1.478 1.399 1.018 0.794 4. Tỷ trọng tương đối của không khí thay đổi theo chiều cao đo được như -0.5 - 1 points linear cubic trong bảng. h (km) 0 1.525 3.050 4.575 6.10 7.625 9.150 ρ 1 0.8617 0.7385 0.6292 0.5328 0.4481 0.3741 0 1 2 3 4 5 6 7 8 spline x Tìm phương trình bậc hai thể hiện quan hệ giữa độ cao và tỷ trọng không khí, từ đó tính tỷ trọng của không khí ở độ cao h = 10.5 km. ρ = + + a bh ch 2 Nguyen Quang Hoang Department of Applied Mechanics Bài tập 21 Nguyen Quang Hoang Department of Applied Mechanics Bài tập 22 3. Bảng số liệu dưới đây cho biết khối lượng m (kg) và suất tiêu hao nhiên liệu y (km/lít) của một số xe ôtô sản suất bởi Ford và Honda năm 1999. Tìm phương trình đường thẳng y = a.m + b thể hiện quan hệ tuyến tính giữa khối lượng và suất 5. Độ nhớt động học μkcủa nước thay đổi theo nhiệt độ T được cho như trong bảng dưới đây tiêu hao nhiên liệu. Model m (kg) y (km/lít) Contour 1310 10.2 Crown 1810 8.1 Escort 1175 11.9 Expedition 2360 5.5 Explorer 1960 6.8 F-150 2020 6.8 Ranger 1755 7.7 Taurus 1595 8.9 Accord 1470 9.8 CR-V 1430 10.2 Civic 1110 13.2 Passport 1785 7.7 Xác định đường bậc ba thể hiện quan hệ này, từ đó tính độ nhớt động học tại các nhiệt độ T = 10◦, 30◦, 60◦, và 90 ◦C T (◦C) 0 21.1 37.8 54.4 71.1 87.8 100 μk (10−3 m2/s) 1.79 1.13 0.696 0.519 0.338 0.321 0.296 Nguyen Quang Hoang Department of Applied Mechanics 23 Nguyen Quang Hoang Department of Applied Mechanics 24 6 16/09/2015 Lecture 6 Giải phương trình vi phân trong Matlab Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology Nguyen Quang Hoang Department of Applied Mechanics Ví dụ giải phương trình vi phân thường Tìm hàm y(t) thỏa mãn phương trình vi phân Chương 6. Lập trình trong Matlab 6.1 Giải phương trình vi phân bậc nhất với ode23 và ode45 6.2 Giải các phương trình vi phân bậc hai 6.3 Bài tập thực hành 1 Nguyen Quang Hoang 2 Department of Applied Mechanics Giải số bằng lệnh ode23 Cú pháp lệnh ode23 cos( ) dyt dt= với điều kiện đầuy(0) 2 = [t, y] = ode23(‘funct_name’, [t0, t_end], y0) Với phương trình đơn giản này ta có thể viết ra được nghiệm chính xác (nghiệm giải tích) , y t t C ( ) sin( ) = + với C là hằng số tích phân phụ thuộc vào điều kiện đầu. Từ điều kiện đầu nhận được hằng số C: y C C (0) sin(0) 2 2 = + = ⇒ = Tên m-file chứa phương trình vi phân cần giải Khoảng thời gian cần giải t từ t0 đến t_end Kết quả tính được ghi vào hai véctơ y t t ( ) sin( ) 2 = + Nguyen Quang Hoang 3 t và y chứa các cặp [t(i), y(ti)] Nguyen Quang Hoang Giá trị điều kiện đầu 4 Department of Applied Mechanics Department of Applied Mechanics 1 16/09/2015 Giải số bằng lệnh ode23 Cú pháp lệnh ode23 [t, y] = ode23(‘funct_name’, [t0:h:t_end], y0) Tên m-file chứa phương trình vi phân cần giải Ví dụ giải phương trình vi phân thường Ví dụ: sử dụng ode23 giải phương trình vi phân cos( ) dyt dt=với điều kiện đầuy(0) 2 = Trong khoảng thời gian t = [0 2*pi] Trước hết tạo m-file chứa phương trình vi phân cần giải: Kết quả tính được ghi vào hai véctơ Khoảng thời gian cần giải t từ t0 đến t_end, với h là bước thời gian Giá trị điều kiện đầu function ydot = eq1(t,y) ydot = cos(t); end Sử dụng lệnh ode23: Ghi lại với tên file là tên hàm eq1.m t và y chứa các cặp [t(i), y(ti)] Nguyen Quang Hoang 5 >> [t, y] = ode23(‘eq1’, [0, 2*pi], 2) Nguyen Quang Hoang 6 Department of Applied Mechanics Ví dụ giải phương trình vi phân thường • Vẽ kết quả và so sánh với nghiệm chính xác >> f = sin(t) + 2; >> plot(t, f,'k-'), hold on >> axis([0 2*pi 0 4]) Department of Applied Mechanics Giải số bằng lệnh ode45 Cú pháp lệnh ode45 [t, y] = ode45(‘funct_name’, [t0, t_end], y0) >> plot(t,y,'ko') >> xlabel('t'), >> ylabel('y(t)') y(t) 4 3 2 1 Tên m-file chứa phương trình vi phân cần giải Khoảng thời gian cần giải t từ t0 đến t_end 0 0 1 2 3 4 5 6 t Nguyen Quang Hoang 7 Kết quả tính được ghi vào hai véctơ t và y chứa các cặp [t(i), y(ti)] Nguyen Quang Hoang Giá trị điều kiện đầu 8 Department of Applied Mechanics Department of Applied Mechanics 2 16/09/2015 Giải số bằng lệnh ode45 Ví dụ: sử dụng ode45 giải phương trình vi phân Giải số bằng lệnh ode23 >> [t, y2] = ode45(‘eq1’, [0:0.01:30], 0); 5 5 ( ) dy y F t với điều kiện đầu y(0) 0 = >> plot(t,y2) dt= − + / ( ) cos( ) t tc F t te t ω − = Trong khoảng thời gian t = [0 30] function ydot = eq2(t,y) tc = 0.1; w = 10; tc = = 0.1, 10 ω Ghi lại với tên Ft = t*exp(-t/tc)*cos(w*t); ydot = -5*y + 5*Ft; end >> [t, y] = ode45(‘eq2’, [0, 30], 0); >> % [t, y] = ode45(‘eq2’, [0:0.01:30], 0); >> plot(t,y) Nguyen Quang Hoang Department of Applied Mechanics Giải số bằng lệnh ode45 Ví dụ: sử dụng ode45 giải phương trình vi phân file là tên hàm eq2.m 9 y2(t) y(t) Nguyen Quang Hoang 10 Department of Applied Mechanics Hệ phương trình vi phân cấp 1 Hệ n phương trình vi phân cấp 1. 5 5 ( ) dy y F t dt= − + với điều kiện đầu y(0) 0 = dy 1 y f t y y = = ( , ,..., ) y y (0) / ( ) cos( ) t tc F t te t ω − = tc = = 0.1, 10 ω 1 1 1 dtdy 2 n Các điều kiện đầu = 1 10 y y (0) = y f t y y = = ( , ,..., ) 2 20 2 2 1 n ... dt Trong khoảng thời gian t = [0 30] ... dy function ydot = odevidu3(t,y, tc, w) n = y y (0) n n 0 y f t y y = = ( , ,..., ) n n n Ghi lại với tên dt 1 Ft = t*exp(-t/tc)*cos(w*t); ydot = -5*y + 5*Ft; end file là tên hàm odevidu3.m Hay viết ở dạng véc tơ y dt y f y y y = = = ( , ), (0) 0 >> [t,y] = ode45(@(t,y)odevidu3(t,y,0.1, 10), [0: 0.01:30], 0); >> plot(t,y) Nguyen Quang Hoang 11 Department of Applied Mechanics dty y y f f f [ , ,..., ] , [ , ,..., ] T T y f = = 1 2 1 2 n n Nguyen Quang Hoang 12 Department of Applied Mechanics 3 16/09/2015 Phương trình vi phân cấp cao Xét 1 phương trình vi phân cấp n. n x x (0) = Phương trình vi phân cấp cao Ví dụ hạ bậc phương trình vi phân cấp 2 sau − = =00 d x x f t x x x x ( ) ( 1) ( , , , ,..., ) n n n dt Đặt các biến mới và hạ bậc Các điều kiện đầu x x (0) = ... = ( 1) ( 1) (0) n n − − x x 0 0 x f t x x f t x x bx cx F t = = − − + Ω ( , , ), ( , , ) sin Đặt các biến mới và hạ bậc y x y x y = = = , y x y x y = = 1 1 2 y x y = , y x = = =1 2 1 = , y x y x f t y y 2 1 2 ( , , ) = = 2 2 3 2 , ... ... ( 1) n − ( 2) n y x y = = Cụ thể trong trường hợp trên y x = nn −− − n n − 1( ) 1( 1) = y x n n = = y x f t y y ( , ,..., ) n n 1 y y = 1 2 y f y = ( , ) t Hay viết gọn lại dạng véc tơ y dt y f y y y = = = ( , ), (0) 0 dty y y y y y f t [ , ,..., ] , [ , ,..., , ( , )] T T y f y = = 1 2 2 3 n n = − − + Ω y by cy F t 2 2 1 0 sin ⎡ ⎤ ⎡ ⎤ y y 1 2 y f y = = ⎢ ⎥ ⎢ ⎥ , ( , )sin t ⎢ ⎥ ⎢ ⎥ − − + Ω ⎣ ⎦ ⎣ ⎦ y by cy F t 2 2 1 0 Nguyen Quang Hoang 13 Department of Applied Mechanics function ydot = daodong1dof(t,y) Ví dụ hệ dao động 1dof Nguyen Quang Hoang Department of Applied Mechanics y x y 1 2 14 % t =1; y = [1; 2] m = 1; b = 0.0; c = 16; >> [t, y] = ode45('daodong1dof', 0 mx cx kx F t sin( ) y x F t cy ky 1sin( ) 2 0 2 1 m F0 = 10; Ome = 4.1; ydot = zeros(2,1); Ft = F0*sin(Ome*t); ydot(1) = y(2); ydot(2) = 1/m*(Ft - b*y(2) - c*y(1)); [0:0.01:100], [1, 2]) >> plot(t,y(:,1)) function ydot = vibr1dof(t,y, m,c,k,F0, Ome) ydot = zeros(2,1); Ft = F0*sin(Ome*t); ydot(1) = y(2); ydot(2) = 1/m*(Ft - c*y(2) - k*y(1)); end 10 5 0 -5 -10 0 10 20 30 40 50 60 70 80 90 100 Nguyen Quang Hoang 15 Department of Applied Mechanics end >> [t, y] = ode45(@(t,y)vibr1dof(t,y, m,c,k,F0, Ome), [0:0.01:100], [1, 2]) >> [t, y] = ode45(@(t,y)vibr1dof(t,y, 1, 0.1,16,10, 4.1), [0:0.01:100], [1, 2]) >> plot(t,y(:,1)) Nguyen Quang Hoang 16 Department of Applied Mechanics 4 16/09/2015 4 4 >> [t, y] = ode45(@(t,y)vibr1dof(t,y, 1, 0.0,16,1, 4.3), [0:0.01:70], [0,0]); 3 >> plot(t,y(:,1)), ylabel('y_1 [mm]'), xlabel('t[s]') 2 2 2 [mm/s] 4 1 0 2 [mm/s] 2 0 4 3 y -2 0 -1 -2 y 2 -4 0 10 20 30 40 50 60 70 -3 -4 0 10 20 30 40 50 60 70 -2 -4 0 10 20 30 40 50 60 70 1 0 -1 t[s] -1 -0.5 0 0.5 1 1 [mm] 1 0.5 0 t[s] -2 -3 -4 -1 -0.5 0 0.5 1 t[s] y -0.5 1 0.5 0 -0.5 -1 1 [mm] y -1 0 10 20 30 40 50 60 70 t[s] Nguyen Quang Hoang Nguyen Quang Hoang 17 Department of Applied Mechanics Giải hệ phương trình vi phân cấp 1 Trước hết xét việc giải hệ hai phương trình vi phân cấp một sau 18 Department of Applied Mechanics Giải hệ phương trình vi phân cấp 1 >> [t,x] = ode45('eqx',[0 10],[0,1]); dx dy điều kiện đầu Nghiệm cần tìm x được ghi lại bằng hai véctơ cột, cột thứ nhất x1 = x(:,1) = − + = − − x y (0) 0, (0) 1 = = 2, x y x xy dt dt Viết lại hệ dạng véc tơ như sau: dx dx 1 2 2 x x x y x x x x x 1 2 1 2 1 1 2 , , và cột thứ hai x2 = x(:,2). Với lệnh vẽ >> plot(t,x(:,1),t,x(:,2),'--'),xlabel('t'), >> axis([0 10 –1.12 1.12]) = = ⇒ = − + = − − dt dt Để giải hệ trên bằng ode45, ta viết một m-file thể hiện vế phải của hệ trên function xdot = eqx(t,x); xdot = zeros(2,1); xdot(1) = -x(1)^2 + x(2); xdot(2) = -x(1) - x(1)* x(2); end % save with file name eqx.m 1 0.5 0 -0.5 - 1 0 2 4 6 8 10 >> plot(x(:,1),x(:,2), 'k-'), 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Nguyen Quang Hoang Department of Applied Mechanics t 19 Nguyen Quang Hoang Department of Applied Mechanics 20 5 16/09/2015 Giải hệ phương trình vi phân cấp 1 Ví dụ dưới đây trình bày việc giải phương trình vi phân bậc hai. Hãy tìm nghiệm của phương trình vi phân sau điều kiện đầu y y t + = 16 sin(4.3 ) y y (0) 0, (0) 0 = = Bằng cách đặt ta nhận được hệ hai ptvp cấp một 1 2 x y x y = = , x x Giải hệ phương trình vi phân cấp 1 >> [t,x] = ode45('eqx2',[0 2*pi],[0,0]); >> plot(t,x(:,1),t,x(:,2),'--'), xlabel('t'), >> axis([0 2*pi –3 3]) >> plot(x(:,1), x(:,2),'k-'), >> xlabel('x_1'), ylabel('x_2‘) 3 = 1 2 = − x t x 2 1 sin(4.3 ) 16 x f x = ( , ) t 2 1 x1 x2 3 2 1 Vế phải của hệ trên được thể hiện trong một m-file như sau: - 3- 2- 10 2 function xdot = eqx2(t,x); x 0 xdot = zeros(2,1); xdot(1) = x(2); xdot(2) = sin(4.3*t)–16*x(1); end % save with file name eqx2.m 0 1 2 3 4 5 6 t -1 -2 -3 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 x1 Nguyen Quang Hoang Department of Applied Mechanics 21 Giải hệ phương trình vi phân cấp 1 Nguyen Quang Hoang 22 Department of Applied Mechanics Giải hệ phương trình vi phân cấp 1 >> [t,x] = ode45('eqx2',[4*pi 20*pi],[0,0]); >> plot(t,x(:,1)),xlabel('t') 1 0.5 0 -0.5 - 1 10 20 30 40 50 60 70 t Nguyen Quang Hoang 23 TT Lệnh solver Dạng bài toán Độ chính xác Trường hợp áp dụng 1 ode45 không cứng (nonstiff) trung bình Hay được sử dụng, là sự lựa chọn đầu tiên (sử dụng công thức R-K45) 2 ode23 không cứng (nonstiff) thấp Sử dụng khi không cần độ chính xác cao, hoặc khi giải bài toán cứng vừa phải. (sdụng công thức R-K23) 3 ode113 không cứng (nonstiff) thấp đến cao (sử dụng công thức Adams Bashforth-Moulton) 4 ode15s cứng (stiff) thấp đến trung bình sử dụng khi ode45 cho độ chính xác không cao do bài toán cứng Nguyen Quang Hoang Department of Applied Mechanics Department of Applied Mechanics 24 6 16/09/2015 Giải hệ phương trình vi phân cấp 1 TT Lệnhsolver Dạng bài toán Độ chính xác Trường hợp áp dụng 5 ode23s cứng (stiff) thấp Khi độ chính xác không của hệ cứng và khi ma trận khối lượng là hằng số 6 ode23t cứng vừa phải (moderately stiff) thấp Giải phương trình có độ cứng vừa phải, nghiệm không giảm dần, cóthể giải được hệ DAEs 7 ode23tb cứng (stiff) thấp giải hệ cứng với độ chính xác khôngcao Ví dụ về hệ cứng (stiffness system). Ref. Chapter 5. Solving Diff. Equations 0 y Ay y y + = = 0, (0) Ma trận A có các trị riêng phân biệt. y v t C t ( ) exp( ) i i i i Ví dụ y y y y , 0 1 y y y 1001 1000 0 1 2 1 1 y y y y y 1000 1001 , 1000 1001 2 2 2 1 2 0 1 | | 0 A I = A =1 Nguyen Quang Hoang Department of Applied Mechanics Bài tập 1000 1001 25 1000 1001 1 2 1, 1000 Nguyen Quang Hoang Department of Applied Mechanics Bài tập h < 2/ λ max 26 Kết hợp m-file và lệnh ode45 giải các phương trình vi phân sau: Phương trình vi phân dao động cưỡng bức của hệ tuyến tính n dof Mq Cq Kq f + + = ( ), t q q q q (0) , (0) = = eqbt1.m eqbt2.m dyy y t t dt= − = = Δ = 2.3 , (0) 0, 10, 0.01 end dx dx 1 2 2 = = − = = = 1 2 1 1 2 2 , , (0) (0) 0, 10 end Bằng cách hạ bậc ta nhận được − q v v M f Cv Kq = = − − 1 , ( ( ) ) t 0 0eqbt7.m ⎡ ⎤ ⎡ ⎤ q v y y = = ⎢ ⎥ ⎢ ⎥ − − ⎣ ⎦ ⎢ ⎥ ⎣ ⎦ ,( ( ) ) t − 1 v M f Cv Kq eqbt3.m eqbt4.m eqbt5.m x x x x x t dt dt dx dx 1 2 = = − = = = x x x x t 2 1 1 2 , , (0) (0) 1, 10 end dt dt dyty y t t dt= − + = = Δ = 1, (0) 1, 10, 0.01 end 2 d y dyy e y y t − − + = = = = 22 , (0) 2, (0) 0, 20 tend dt dt Hãy viết các m-file thể hiện các phương trình vi phân và thực hiện việc giải số tìm dao động của hệ, với ⎡ ⎤ ⎡ ⎤ − 12 0 0 10 10 0 ⎢ ⎥ ⎢ ⎥ = = − − ⎢ ⎥ ⎢ ⎥ M C 0 15 0 , 10 20 10 , ⎢ ⎥ ⎢ ⎥ − 0 0 20 0 10 10 ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎡ ⎤ − 2000 1000 0 5sin10 t ⎢ ⎥ ⎢ ⎥ = − − = ⎢ ⎥ ⎢ ⎥ K f 1000 2000 1000 , ( ) 0 t ⎢ ⎥ ⎢ ⎥ − 0 1000 1000 0 ⎣ ⎦ ⎣ ⎦ Nguyen Quang Hoang Department of Applied Mechanics 27 Nguyen Quang Hoang Department of Applied Mechanics 28 7 16/09/2015 Bài tập function ydot = daodong3dof(t,y); ydot = zeros(6,1); q = y(1:3); v = y(4:6); M = [12, 0, 0; 0, 15, 0; 0, 0, 20]; C = [10 -10 0; -10 20 -10; 0 -10 10]; K = [2, -1 0; -1 2 -1; 0 -1 1]*1000; ft = [5*sin(10*t); 0; 0]; qdot = v; vdot = inv(M)*(ft-C*v-K*q); ydot = [qdot; vdot]; Bài tập Phương trình vi phân chuyển động của tay máy 1 dof: 2 2 2 ( ) ( / ) cos( ) ( / ) O m m e a m a J J r K K r R br mgl K r R U + + + + = ϕ ϕ ϕ Hãy viết các m-file thể hiện các phương trình vi taymay1dof.m phân và thực hiện việc giải số tìm chuyển động của tay máy: Tay máy end % save with file name daodong3dof.m U >> [t,y] = ode45(‘daodong3dof',[0 20],[0,0,0,0,0,0]); >> plot(t,y(:,1)),xlabel('t') Nguyen Quang Hoang 29 Động cơ DC Hộp giảm tốc Nguyen Quang Hoang 30 Department of Applied Mechanics Bài tập Xét hệ dao động n dof, cho M, C, K, f(t): 1. Thiết lập phương trình mô tả hệ 2. Xét dao động tự do không cản Mq Cq Kq f + + = ( ), t Department of Applied Mechanics Bài tập Động lực học robot dạng chuỗi hở s M q q C q q q Dq F q g q u + + + + = t ( ) ( , ) sgn( ) ( ) ( ), 1. Thiết lập phương trình mô tả hệ • Tần số dao động riêng, các dạng dao động riêng • Vẽ các dạng dao động riêng của hệ • Kiểm tra tính trực giao của ma trận dạng riêng thông qua M và K 3. Xét dao động cưỡng bức có cản • Đồ thị biên độ tần số • Tính dao động cưỡng bức [thành phần dao động cùng tần số cưỡng bức] • Tính dao động của hệ theo thời gian q(t) • Vẽ quĩ đao pha của từng tọa độ suy rộng 4. Mô phỏng tìm q(t) với các kích động tùy ý [xung, bước nhảy, f(t) – tùy ý] 2. Cho q(t) tìm lực/mô men u(t) • Vẽ đồ thị q(t), q’(t), q’’(t) và u(t) 3. Cho lực/mô men u(t) tìm chuyển động q(t) • Vẽ đồ thị q(t), q’(t), và u(t) ( ) ( ) ( ) p d u g q K q q K q t = − − − 0 Nguyen Quang Hoang Department of Applied Mechanics 31 Nguyen Quang Hoang Department of Applied Mechanics 32 8 16/09/2015 Bài tập Động lực học robot song song C Bài 4. Cho một rôbốt song song phẳng ba q3 T M q q C q q q Dq g q u q t ( ) ( , ) ( ) ( ) ( ) , bậc tự do, chuyển động trong mặt phẳng C1 ngang. Hãy thiết lập phương trình vi phân + + + = + ( ) Φ λ chuyển động cho rôbốt với các tọa độ suy p3 C2 ϕ q 0 = a p q q q x = [ , , ] T T T T p2 rộng dư. Viết chương trình Matlab để mô phỏng động lực học cho rôbốt. Biết rằng B1 các tam giác ABC và (ABC)2là các tam φ B2 y 1. Thiết lập phương trình mô tả hệ 2. Động học thuận/ngược robot song song giác đều. 1 2 3 1 2 3 2 2 [ , , , , , , , , ]T A A q q q q p p p x y A B q1 A2 A1 p1 q2 x 3. Cho chuyển động của bàn máy [khâu cuối], x(t), tìm lực/mô men u(t) 4. Cho lực/mô men u(t) tìm chuyển động bàn máy x(t) Nguyen Quang Hoang 33 Department of Applied Mechanics Hình bài 4 Nguyen Quang Hoang 34 Department of Applied Mechanics 9 9/16/2015 Lecture 8 Giới thiệu về Simulink Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology 1 Khái niệm về simulink Chương 8. Giới thiệu về Simulink 8.1 Khái niệm về simulink 8.2 Nguyên lý hoạt động và việc thực hành trong simulink Khởi động simulink Các khối trong simulink 8.3 Một số ví dụ đơn giản a) Mô hình hóa một phương trình bằng sơ đồ khối b) Mô phỏng một quá trình động học c) Mô hình hóa một hệ động lực liên tục đơn giản d) Mô tả hệ dao động một bậc tự do e) Mô phỏng số tay máy một bậc tự do 8.4 Đơn giản sơ đồ simulink Đơn giản sơ đồ simulink bằng khối Fcn Đơn giản sơ đồ simulink bằng khối subsystem Kết hợp simulink và script file (m-file) 8.5 Xử lý kết quả mô phỏng 8.6 Bài tập thực hành 2 Biểu diễn sơ đồ khối của biểu thức • Công cụ simulink trong Matlab cho phép mô phỏng các hệ động lực với một giao diện đồ họa đặc biệt bằng các khối kết nối với nhau. Mục đích của chương này là giới thiệu về simulink thông qua các ví dụ. Sau khi thiết lập các phương trình vi phân mô tả hệ, một sơ đồ khối tương ứng được đưa ra và từ đó sử dụng các khối của simulink để mô phỏng hệ. Xét phương trình vi phân sau mv bv Ku t + = ( ) Biểu thức nghiệm tìm được dạng 1 ⇒ = − v Ku t bv [ ( ) ] m t t ∫ ∫ 1 v t v v d v Ku bv d ( ) (0) ( ) (0) [ ( ) ] = + = + − τ τ τ τ • Simulink là một phần mở rộng của MATLAB được phát triển bởi Mathworks Inc. Simulink là một gói chương trình để mô hình hóa, mô phỏng và phân tích các hệ động lực. Với simulink ta có thể mô hình hóa các hệ thống tuyến tính và phi tuyến liên tục theo thời gian, rời rạc theo thời gian, hoặc hỗn hợp cả liên tục và rời rạc theo thời gian. Các hệ cũng có thể là đa tốc độ khác nhau, nghĩa là các phần khác nhau có thể được lấy mẫu hoặc cập nhật số liệu ở các tốc độ khác nhau. 3 Biểu diễn dạng sơ đồ khối 0 0 v(0) m vt() 4 1 9/16/2015 • Khởi động simulink Khái niệm về simulink Các khối trong simulink – Kích đúp vào biểu tượng hoặc gõ lệnh >> simulink • Biểu tượng Các khối trong simulink của Simulink 5 6 Các khối trong simulink 7 8 2 9/16/2015 Các khối trong simulink Các khối trong simulink 9 10 Các khối trong simulink Các khối trong simulink 11 12 3 9/16/2015 Các khối trong simulink Một số ví dụ đơn giản Mô hình hóa một phương trình bằng sơ đồ khối Ví dụ ta cần chuyển đổi từ thang nhiệt độ C sang thang nhiệt độ Fahrenheit, phương trình mô tả quan hệ nhiệt độ giữa hai thang đo này như sau 932 T T F C = + 5 13 Một số ví dụ đơn giản Các khối cần được kéo vào model window Vị trí các khối trong thư viện Simulink Ramp Sources Constant Sources Gain Math Operation Add hoặc Sum Math Operation Scope Sinks 14 Một số ví dụ đơn giản Nút play để chạy chương trình 15 16 4 9/16/2015 Mô phỏng một quá trình động học Trong ví dụ này ta sẽ sử dụng Simulink để mô tả một phương trình động học saux t A t ( ) sin( ) = + ω ϕ Mô phỏng một quá trình động học Trong ví dụ này các thông số được nhận các giá trị A = = = 2, 4, /3 ω ϕ π Các khối cần được kéo vào model window Vị trí các khối trong thư viện Simulink Constant Sources Ramp Sources Gain Math Operation Sum Math Operation Product Math Operation Trigonometry Function Math Operation Mux Signal Routing Scope Sinks 17 Mô phỏng một quá trình động học x t A t ( ) sin( ) = + ω ϕ A = = = 2, 4, /3 ω ϕ π 18 Giải phương trình vi phân Xét phương trình vi phân mô tả một hệ động lực liên tục như sau u t() v m bv x t x t u t ( ) 2.5 ( ) ( ) = − + x t A t ( ) cos( ) = + ω ϕ A = = = 2, 4, /3 ω ϕ π 19 1 ⇒ = − v Ku t bv [ ( ) ] m 20 5 9/16/2015 Mô tả hệ dao động một bậc tự do Phương trình vi phân chuyển động x x Mô tả hệ dao động một bậc tự do Phương trình vi phân chuyển động m f t() mx t f t ( ) ( ) = Viết lại phương trình vi phân chuyển động trên dưới dạng sau k c m f t() mx t cx t kx t f t ( ) ( ) ( ) ( ) + + = Viết lại phương trình vi phân chuyển động trên dưới dạng sau = 1 Hình 8-12. 1 = − − x f t [ ( )] m x f t cx kx [ ( ) ] m x x x ∑1m f t() 1x x x f t() 1 1 m 1 s 1 s 21 s −c −k s 22 Mô tả hệ dao động một bậc tự do Mô tả hệ dao động một bậc tự do mx t cx t kx t f t ( ) ( ) ( ) ( ) + + = 23 24 6 9/16/2015 Mô tả hệ dao động một bậc tự do mx t cx t kx t f t ( ) ( ) ( ) ( ) + + = Mô hình hóa và mô phỏng hệ cơ điện Xét tay máy 1 dof được dẫn động bởi động cơ điện một chiều (DC motor) Các thông số của động cơ điện một chiều: • điện trở, Ra • cuộn cảm có hệ số tự cảm, La • hằng số phản sức điện động, Ke • hằng số mômen, Km • mô men quán tính của rô-to, Jm • Hệ số cản ổ trục, b1 Các thông số tay máy: • khối lượng m, khối tâm C, OC = L. U DC động cơ Tay máy Hộp giảm tốc • mômen quán tính khối đối với trục quay O, Jo; • hệ số cản tại ổ trục b2 >> load ketqua.mat t = ans.time(:,1) u = ans.data(:,1) x = ans.data(:,2) Hộp giảm tốc: • tỷ số truyền của r r r = = > ω ω 1 2 2 1 / / 1 Tương tác cơ – điện: = = θ , V K M K i emf e m plot(t,u), 25 Mô hình hóa và mô phỏng hệ cơ điện 26 Mô hình hóa và mô phỏng hệ cơ điện Ra La U Μ Tload r2 M1 r1 g M1 O C1 ϕ Đối với hộp giảm tốc (không khối lượng, không ma sát) ta có công suất vào bằng công suất ra T M load 1 θ ϕ = r r r r r 1 2 2 1 θ ϕ θ ϕ ϕ = ⇒ = = / Vemf Động cơ DC θ Tload T M M r load 1 1 = = ϕ θ/ / g Phương trình chuyển động của tay máy do tác C1 dụng của mômen Phương trình điện (từ phương trình cân bằng điện áp của một mạch vòng) 1 2 cos( ) J M b mgl Oϕ ϕ ϕ = − − M1 ϕ di L R i U K dt+ = − ω U U U V R L emf = + + ⇒ a a e Phương trình cơ (phương trình chuyển động quay của rôto) d L m F i q = e O Khử Tload & M1 dt= Σ ⇒ ( ) z z k 1, J K i b T m m load ω ω ω θ = − − = 27 28 7 9/16/2015 Mô hình hóa và mô phỏng hệ cơ điện Các PTVP mô tả 1, J K i b T m m load ω ω ω θ = − − = Mô hình hóa và mô phỏng hệ cơ điện Hệ 3 phương trình vi phân cấp 1 mô tả qua hệ đầu vào U và đầu ra ϕ di L U R i K rq = − − a a e dt 2 2 1 2 ( ) ( ) cos( ) J J r q K ri b r b q mgl di L R i U K i q dt+ = − = ω , + = − + − O m m ϕ a a e e ϕ = q 1 2 cos( ) J M b mgl Oϕ ϕ ϕ = − − Nếu bỏ qua Ladi/dt ≈ 0 2 ( ) J J rK K r R b r b + + Tay máy Các ràng buộc O m ϕ 2 2 ( / ) θ ω ϕ = = r Khử Tload & M1 m e a + + 1 2 ϕ + = mgl K r R U cos( ) ( / ) T M load 1 θ ϕ = d 1 1 / / T M M r loa= = ϕ θ 29 ϕ U ϕϕ, Hệ thống m a U DC động cơ Hộp giảm tốc 30 Mô hình hóa và mô phỏng hệ cơ điện Mô hình hóa và mô phỏng hệ cơ điện Các khối cần được kép vào model window Vị trí các khối trong thư viện Simulink Constant Sources Sum Math Operation Product Math Operation Divide Math Operation Gain Math Operation Trigonometric functions Math Operation Integrator Continuous Scope Sinks 31 32 8 9/16/2015 % Can nhap thong so he: % - DC motor: Ra, La, Ke, Km, Jm, b1 Mô hình hóa và mô phỏng hệ cơ điện Bài tập: Ghép 2 tay máy 1dof thành 1 tay máy phẳng 2dof. R K K L b J % - Tay may: m, lc, g, JO, b2 , , , , , 2a 2 2e 2a 21 2 Tay máy % - hop giam toc: ti so truyen r U , 1 Hệ thống ϕ ϕ , , ϕ ϕ 1 1 r 2 m m % - DC motor: Ra=1, La=0.0001, Ke=1, Km=1, U 2 , , 2 2 m J l b , , , 2 2 2 22 C C U Hộp giảm tốc Khâu 2 (Link2) Jm=0.001, b1=0.01; % - Tay may: m=5, lc=0.5, g=9.81, JO=0.2, b2=0.2; % -hop giam toc: ti so truyen U2 R K K L b J , , , , , 1a 1 1e 1a 11 1 Tay máy DC động cơ r = 2.0 r 1 m m U Hộp giảm tốc Khâu 1 (Link1) m J l b , , , 1 1 1 12 C C 33 U1 DC động cơ 34 Mô hình hóa và mô phỏng hệ cơ điện 2 ( ) J J rK K r R b r b + + O m ϕ 2 2 ( / ) m e a + + 1 2 ϕ + = mgl K r R U cos( ) ( / ) ϕ U ϕϕ, Hệ thống Điều khiển m a U DC động cơ Tay máy Hộp giảm tốc Xét hai bộ điều khiển 0 ϕ 1. ( / ) ( ) U K r R k k 1 m a p d U K r R k k mglPD-controller 0 2. ( / ) ( ) cos 1 m a p d 0 35 PD-controller + bù trọng lực 36 9 9/16/2015 Mô hình hóa và mô phỏng hệ cơ điện Bài tập 37 38 Bài tập Bài tập 6. Cho sơ đồ mạch điện như trên hình vẽ. Áp dụng định luật Kirchhoff ta nhận được phương trình vi phân của mạch điện như sau: di L R i R i i E t L 1 + + − = ( ) ( ) 1 1 2 1 2 E(t) i1 i1 R1 i2 R2 C dtdi q L R i i 2 2 − − + = ( ) 0 2 1 2 dt C dq i L i2 2 dt = 2 Cho biết các số liệu R R L C 1 2 = Ω = Ω = = 4 , 10 , 0.032H, 0.53F và điện áp đặt lên mạch có dạngHãy giải và vẽ ra đồ thị các dòng điện⎧ < < ⎪ 20 V khi 0 0.005 s t E t t ( ) 0, khi 0.005 s = ⎨ ≥ ⎪⎩1 2 i t i t ( ), ( ) t = 0...0.05 s 40 39 10 9/16/2015 Bài tập R1 − dt i1i Đơn giản hoá sơ đồ simulink 1. Đơn giản sơ đồ simulink bằng khối Fcn 2. Đơn giản sơ đồ simulink bằng khối subsystem d ∑1L 1 1 E t() s 3. Kết hợp simulink và script file (m-file) R2 - + dt i2 d 1 i - + di 1 L E t R i R i i = − − − ( ) ( ) 2 L 1 s 1 1 2 1 2 dt q2 dtdi q 2 2 L R i i = − − ( ) 2 1 2 dt C dqi d C−1 − 2 1 s q 2 dt = 2 41 Đơn giản hoá sơ đồ simulink 42 Đơn giản hoá sơ đồ simulink 1. Đơn giản sơ đồ simulink bằng khối Fcn 1. Đơn giản sơ đồ simulink bằng khối Fcn x t = + 5sin(4 /3) π 43 44 11 9/16/2015 Đơn giản hoá sơ đồ simulink 2. Đơn giản sơ đồ simulink bằng khối subsystem di U R i K rq 1 = − − ( ) Đơn giản hoá sơ đồ simulink di dt / dt L a a e dq K ri b r b q mgl = − + − ⎡ ⎤ 2 1( ) cos( ) ϕ + ⎣ ⎦ 2 1 2 dt J J r ( ) m ϕ dq O m dt = 45 M1 ϕ 46 Đơn giản hoá sơ đồ simulink 47 Đơn giản hoá sơ đồ simulink 48 12 9/16/2015 Đơn giản hoá sơ đồ simulink Đơn giản hoá sơ đồ simulink 3. Kết hợp simulink và script file (m-file) Xét trường hợp sử dụng simulink mô phỏng hệ có phương trình động lực phức tạp, như robot công nghiệp nhiều bậc tự do M q q C q q q g q u ( ) ( , ) ( ) + + = g M2 A C1 q2 C2 E q v = ( ) [ ( , ) ( )] − 1 M1 O v M q u C q v v g q = − − q1 49 Đơn giản hoá sơ đồ simulink 3. Kết hợp simulink và script file (m-file) u⎡ ⎤ M = ⎢ ⎥ ⎣ ⎦ 1 M 2 4. Xử lý kết quả 0 y f y u y y = = ( , ), (0) 50 Đơn giản hoá sơ đồ simulink Soạn một MATLAB function [m-file] function ydot = xvdot(in) ydot = …. ; end 51 load sim_results.mat t=qv(1,:); q1=qv(2,:); q2=qv(3,:); q1_dot=qv(4,:); q2_dot=qv(5,:); figure(1) plot(t,q1,'-k','linewidth',2), grid on xlabel ('t [s]'); ylabel ('q_1 [rad]'); figure(2) plot(t,q2,'-k','linewidth',2), grid on xlabel ('t [s]'); ylabel ('q_2 [rad]'); % save ve_q_t.m 52 13 9/16/2015 Đơn giản hoá sơ đồ simulink Đơn giản hoá sơ đồ simulink Ví dụ: kết hợp simulink và m-file Ví dụ: kết hợp simulink và m-file x ux O M q q C q q q g q Bu ( ) ( , ) ( ) + + = q v = ( ) [ ( , ) ( )] − 1 v M q u C q v v g q = − − mt l ul M q m m m l t p p cos Con lắc elliptic O x ux ul mt z θ mp l ( )cos 2 m l m l p p m lp C q q 0 sin ( , )0 0 g q 0 1 z θ mp 53 ( )sin m gl p B 0 54 Đơn giản hoá sơ đồ simulink function xdot = xdot_elliptic_pendul(ins) Đơn giản hoá sơ đồ simulink Kết hợp simulink-model với lệnh sim % cac tham so he m1 = 5; % kg m2 = 2; % kg l = 0.5; % m gr = 9.81; % m/s^2 d1 = 6.2; d2 = 0; % damping coefs % ins = [1 2 3 4 5]'; % for test % viet ra cac phuong trinh vi phan q = ins(1:2); % bien chua toa do suy rong qdot = ins(3:4); % bien chua van toc suy rong u = ins(5); % bien vao (bien dieu khien) M = [m1+m2, m2*l*cos(q(2)); m2*l*cos(q(2)), m2*l^2 ]; C = [0, -m2*l*qdot(2)*sin(q(2)); 0, 0]; D = diag([d1, d2]); % damping matrix g = [0; m2*l*gr*sin(q(2))]; B = [1; 0]; % he hut dan dong % B = [1, 0; 0, 1]; % he đủ dan dong vdot = inv(M)*(B*u- (C+D)*qdot - g); xdot = [qdot; vdot]; 55 • Thực hiện mô phỏng mà không cần mở simulink • Sử dụng được mô hình simulink trong một m-file • Thuận tiện cho các bài toán lặp như tìm tham số tối ưu cho bộ điều khiển bằng GA, PSO, Nelder-Mead,…, • Cú pháp: xem >> help sim 5. Mask and unmask ?? 56 14 9/16/2015 Đơn giản hoá sơ đồ simulink 57 Đơn giản hoá sơ đồ simulink % Can nhap thong so he: % - DC motor: Ra, La, Ke, Km, Jm, b1 % - Tay may: m, lc, g, JO, b2 % - hop giam toc: ty so truyen r Ra=1, La=0.0001, % - DC motor Ke=11, Km=1, Jm=0.001, b1=0.00000001; m=5, lc=0.5, g=9.81, JO=0.2, b2=2e-7; % - Tay may r = 2.0; % -hop giam toc: ti so truyen kp = 100; kd = 50; % -tham so bo dieu khien [t, y] = sim('taymay_1DOF_DCmotor',[0, 8]) plot(t,y(:,1), 'r-', 'linewidth',2), grid on, xlabel('t [s]'), ylabel('\phi (t) [rad]') % su dung neu dua mdl vao trong mot function options = simset('SrcWorkspace','current'); [t, y] = sim('modelname',[],options) 58 15 12/2/2015 Phân tích hệ động lực Điều khiển tối ưu LQR Linear Quadratic Regulators Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology Nguyen Quang Hoang 1 Department of Applied Mechanics Chương x. Xây dựng mô hình 1. Nêu bài toán điều khiển tối ưu 2. Cơ sở lý thuyết 3. Áp dụng điều khiển tối ưu hệ tuyến tính 4. Điều khiển tối ưu LQR 5. Ví dụ • Ref. Alok Sinha : Linear Systems- Optimal and Robust Control, CRC press, 2007. Nguyen Quang Hoang 2 Department of Applied Mechanics 1 12/2/2015 Nêu bài toán điều khiển tối ưu Để điều khiển hệ từ trạng thái x(0) về gốc tọa độ bằng điều khiển phản hồi u = u(x), có thể có rất nhiều thuật toán điều khiển khác nhau. Ngay cả đối với một thuật toán điều khiển, chẳng hạn u = -Kx, thì có rất nhiều khả năng chọn ma trận phản hồi K. Bài toán điều khiển tối ưu sẽ đưa ra phương thức lựa chọn bộ thông số của bộ điều khiển để một chỉ tiêu chất lượng nào đó được tối thiểu hóa. Trong phần này sẽ nêu bài toán điều khiển tối ưu, và phương pháp thiết kế tối ưu LQR. Nguyen Quang Hoang 3 Department of Applied Mechanics Nêu bài toán điều khiển tối ưu Xét hệ động lực x fxu x x  = = 0 ( , ), (0) y gx Điều khiển phản hồi trạng thái u ux = ( ) = ( ) Các tiêu chuẩn tối ưu cần được cực tiểu có thể là t t T T f f ò ò J dt dt e Qe x x Q x x ( )( ) = =- - (1) 1 0 0 t T d d ò f (2) J dt ( )( ) y y Qy y =- - 2 0 ò d d t T f (3) J dt = 3 0 ò t f u Ru 1 (4) J dt = 4 0 0, 0 T T QQ RR =³ => Nguyen Quang Hoang 4 Department of Applied Mechanics 2 12/2/2015 Nêu bài toán điều khiển tối ưu Để giải quyết tình trạng khó xử, chúng ta có thể thỏa hiệp giữa hai mục tiêu mâu thuẫn nhau bằng cách giảm thiểu các chỉ số thực hiện mà là một sự kết hợp lồi của J1 và J3, (1 ) ( ) ( ) ft T T d d J J J dt l l =- + = - - + é ù ê ú ò ë û x x Q x x u Ru (5) 1 3 0 Trong các ứng dụng nhất định, chúng ta có thể muốn trạng thái cuối cùng x(tf) là càng gần 0 càng tốt. Sau đó, một chỉ số thực hiện có thể được giảm thiểu là ( ) ( ), 0 T T f f x Fx F F t t = > (6) Chúng ta có thể kết hợp các biện pháp thực hiện (1), (3), và (6) khi mục tiêu điều khiển là để giữ cho trạng thái “nhỏ", sự điều khiển "không quá lớn," và trạng thái cuối cùng càng gần 0 càng tốt. Chỉ số thực hiện kết quả là () () ( ) ( ) ft T TT ff d d J t t dt = + - -+ é ù ê ú ò ë û x Fx x x Q x x u Ru (7) 0 Nguyen Quang Hoang 5 Department of Applied Mechanics Nêu bài toán điều khiển tối ưu Xét hệ động lực tuyến tính, với cặp (A,B) là điều khiển được x Ax Bu x x  0 =+ = , (0) y Cx = Tiêu chuẩn 2 trở thành 2 00 0, ff f tt t T TT T T J dt dt dt == = = ò ò ò y y x C Cx x Qx Q C C Bài toán tối ưu được đặt ra là: Tìm luật điều khiển u(t) để làm tối thiểu một trong các chỉ số thực hiện nêu trên, hoặc tối thiểu một chỉ tiêu thỏa hiệp trong các trường hợp: 1. Cố định thời gian cuối tf ; 2. Không cố định thời gian cuối tf ; 3. Không có ràng buộc đặt lên u(t); 4. Có ràng buộc đặt lên điều khiển, ui(t) có cận trên và cận dưới. Nguyen Quang Hoang 6 Department of Applied Mechanics 3 12/2/2015 Cơ sở lý thuyết điều khiển tối ưu Xét hệ động lực điều khiển được 0 x fxu x x  = = ( , ), (0) Và tiêu chuẩn tối ưu sau cần được cực tiểu (hoặc cực đại) [ ( )] ( , ) ft f J G t F dt = + ò x xu 0 Các ràng buộc đặt lên điều khiển ui(t), cận trên và cận dưới , , , 1,2,..., im i iM u uu i m ££ = dim( ) dim( ) , dim( ) xf u == = n m Nguyen Quang Hoang 7 Department of Applied Mechanics Cơ sở lý thuyết điều khiển tối ưu Các điều kiện cần cho tối ưu: Nếu điều khiển u0(t) là nghiệm của bài toán tối ưu và x0(t) là véc tơ trạng thái tương ứng: 0 00 x fx u  ( ) ( ( ), ( )), t tt = Khi đó ta sẽ đánh giá hàm chất lượng điều khiển (chỉ tiêu tối ưu), khi điều khiển thay đổi một lượng nhỏ δu(t), tức là 0 uu u ( ) ( ) ( ), ttt = + d Tiêu chuẩn tối ưu thay đổi một lượng là d d JJ t t J t = +- [ ( ) ( )] [ ( )] uu u 0 0 ¶ ¶¶ ¶ é ùé ù G FF G t dt F t t t () () () ft = + + ++ ê úê ú ¶ ¶¶ ¶ êë ûë û ò d dd d x xu f f f ff x xu x 0 (8) Nguyen Quang Hoang 8 Department of Applied Mechanics 4 12/2/2015 Cơ sở lý thuyết điều khiển tối ưu Tính gần đúng theo xấp xỉ bậc nhất nghiệm δx(t) ứng với δu(t): 0 00 x x fx x u u   ( ) ( ) ( ( ) ( ), ( ) ( )), t t t tt t += + + d dd == + ¶ ¶  f f x xx u d dd d ¶ ¶ ( ) ( ) ( ), d t tt dt x u Nhân hai vế phương trình trên với λT(t) và lấy tích phân từ 0 đến tf được éù é ù ¶ ¶ êú ê ú - += ¶ ¶ ëû ë û ò ò x ff x u ( ) () () () () 0 f f t t T T d t dt t t t dt d d d l l (9) dt 0 0 x u Cộng (8) với (9), thực hiện tích phân từng phần số hạng đầu của (9) được é ù éù ¶ ¶ G G J Ft t t t t T T = + + +- ê ú êú l l ( ) ( ) (0) (0) ( ) ( ) fx x d dd d f ff f f x x ¶ ¶ ë û ëû é ù ¶ ¶ T t l ( ) () () f H H dt t t dt + ++ ê ú ¶ ¶ ë û ò xu x dd d (10) x u 0 dt Nguyen Quang Hoang 9 Department of Applied Mechanics Cơ sở lý thuyết điều khiển tối ưu Trong (10) H được gọi là Hamilton, định nghĩa như sau ( , ) () ( , ) T HF t = + xu fxu l Do λ(t) là tùy ý, nên ta chọn nó thỏa mãn T dt H l l ( ) (, , ) ¶ = - ¶x u (11)(12) dt x Số hạng ngoài dấu tích phân của (10) là các số hạng liên quan điều kiện biên, chúng sẽ không xuất hiện trong các bài toán cụ thể: Nếu tf và x(0) xác định (cố định) thì ta có δtf = 0 và δx(0) = 0; số hạng thứ 3 ngoài dấu tích phân triệt tiêu nếu chọn çæ ö ¶ ÷ = ç ÷ T G t l (13) ( ) f t t ç ÷ ç ¶ ÷ è ø x = f Nguyen Quang Hoang 10 Department of Applied Mechanics 5 12/2/2015 Cơ sở lý thuyết điều khiển tối ưu Do đó, phương trình (10) được viết thành é¶ ù = ê ú ê ¶ ú ë û ò u ( ) ft H d d J t dt (14) u 0 Cực tiểu hàm J Nếu u0(t) là nghiệm tối ưu (J đạt cực tiểu) thì d d J t ³ " 0, ( ) u Trường hợp 1. Không có ràng buộc đặt trên u(t), thì , Ht ¶ = " ¶0 uĐây là điều kiện cần để tối ưu J. Trường hợp 2. Khi có ràng buộc đặt trên u(t), [bị chặn trên và dưới], thì ¶ ¶ =£ =³ ¶ ¶ H H ,0 , ,0 , , 0, & , 0 i iM i im uu uu u u i i Nguyen Quang Hoang 11 Department of Applied Mechanics Cơ sở lý thuyết điều khiển tối ưu Do đó, phương trình (10) được viết thành é¶ ù = ê ú ê ¶ ú ë û ò u ( ) ft H d d J t dt (14) u 0 Cực đại hàm J Nếu u0(t) là nghiệm tối ưu (J đạt cực đại) thì d d J t £ " 0, ( ) u Trường hợp 1. Không có ràng buộc đặt trên u(t), thì , Ht ¶ = " ¶0 uĐây là điều kiện cần để tối ưu J. Trường hợp 2. Khi có ràng buộc đặt trên u(t), [bị chặn trên và dưới], thì ¶ ¶ =³ =£ ¶ ¶ H H ,0 , ,0 , , 0, & , 0 i iM i im uu uu u u i i Nguyen Quang Hoang 12 Department of Applied Mechanics 6 12/2/2015 Cơ sở lý thuyết điều khiển tối ưu Các tính chất của Hamilton đối với hệ autonomous x fxu  = ( , ), l dH H d H d H d ¶¶¶ x u H H =  (, , ) x u l =++ ¶¶¶ l (15) dt dt dt dt x u T dt H lvà theo cách chọn( ) (, , ) l l Nhưng (, ) ¶H T = ¶ f x u dt nên dH H d ¶ = ¶u dt dt u ¶ = - ¶x u x (16) Trường hợp 1. Không có ràng buộc đặt trên u(t), thì 0 ¶H = ¶u Trường hợp 2. Khi có ràng buộc đặt trên u(t), thì u = umax hoặcu=umin nên 0 ddt = u dH Do đó, đối với điều khiển tối ưu ta có hay 0 H = const. dt = (16) Nguyen Quang Hoang 13 Department of Applied Mechanics Cơ sở lý thuyết điều khiển tối ưu Trường hợp đặc biệt: Thời điểm cuối không xác định δtf ≠ 0. é ¶ ¶ ù éù G G J Ft t t t t T T = + + +- ê ú êú d dd d l l ( ) ( ) (0) (0) ( ) ( ) fx x f ff f f x x ¶ ¶ ë û ëû é ù ¶ ¶ T t l H H dt t t dt ( ) () () f + ++ ê ú ¶ ¶ ë û ò dd d xu x 0 (10) x u Để δJ không phụ thuộc δtf, ta cho G Ft t ¶ + = ¶ f () () 0 f f x dt (17) Nếu điều kiện cuối về trạng thái không xác định δx(tf), điều kiện (13) được thỏa mãn và do đó T HF t = + ( , ) () ( , ) xu fxu l T Ht Ft t t () () ()() 0 =+ = l (18) f f f ff Nguyen Quang Hoang 14 Department of Applied Mechanics 7 12/2/2015 Cơ sở lý thuyết điều khiển tối ưu Trường hợp nếu trạng thái cuối không xác định, với xấp xỉ bậc nhất trong khai triển Taylor, ta có (19) 0 () () () () f f f ff d d x xx f t t t tt = - =- Với các biểu thức trên, ta có é¶ ù é ù =+ + + ê ú ê ú ë û ê ¶ ú ë û ò f xu ( ) ( ) ( ) (0) (0) ( ) ft T T H d dd d J F t t t t t dt l l f f ff (20) 0 u Để loại bỏ số hạng liên quan điều kiện biên khi tf không xác định () () ()() 0 T Ht Ft t t f f ff =+ = l f Vì H(t) = const, nên H(t) = H(tf)=0dọc theo quỹ đạo tối ưu đối với hệ autonomous, khi thời điểm cuối tf không xác định. Nguyen Quang Hoang 15 Department of Applied Mechanics Điều khiển tối ưu LQR hệ tuyến tính Xét hệ động lực tuyến tính, với cặp (A,B) là điều khiển được 0 x Ax Bu x x  =+ = , (0) Và tiêu chuẩn tối ưu sau cần được cực tiểu 1 1 () () 2 2ft T TT (21) f f J t t dt = ++ é ù ê ú ò ë û x Fx x Qx u Ru (22) 0 0, 0, 0 TT T F F QQ RR =³ =³ => , , , 1,2,..., im i iM u uu i n ££ = Hàm Hamilton trong bài toán này là { } ( ) 1 () () () 2TT T H tt t = ++ + x Qx u Ru Ax Bu l (23) Nguyen Quang Hoang 16 Department of Applied Mechanics 8 12/2/2015 Điều khiển tối ưu LQR hệ tuyến tính Lời giải bài toán tối ưu vòng mở: Từ điều kiện cần cho nghiệm tố ưu ¶= += ¶ Ht Ru B T l 0 () 0 u T t - 1  =- u RB l (24) ( ) Động học của λ(t) đưa ra từ (12) cùng với điều kiện cuối (13) dt =- - = Qx A Fx l l l d tttt (25) ( ) ( ), ( ) ( ) Tf f Kết hợp (21)-(24)-(25) ta được - =- =  1 T x Ax BR B x x 0 ( ) ( ), (0) t t l  T l ll ( ) ( ), ( ) ( ) t t tt =- - = Qx A Fx f f é ù é ù - - é ù éù ê ú = = ê ú ê ú êú ê ú ê ú êú ê ú - - ê ú êú ë û ë û ë û ëû 1 () () () T  t tt x A BR B x x M l Q A l l (26)  T ( ) () () t t t Nguyen Quang Hoang 17 Department of Applied Mechanics Điều khiển tối ưu LQR hệ tuyến tính Giải (26) tìm được nghiệm éù é ù êú ê ú = M x x ( ) (0) t te l l ( ) (0) t ëû ë û Xác định điều kiện đầu λ(0) ta viết é ù E E Mt t t () () = ê ú ê ú ê ú ë û 11 12 et t E E () () 21 22 Suy ra x Ex E ( ) ( ) (0) ( ) (0) tt t = + 11 12 l (27) l l ( ) ( ) (0) ( ) (0) tt t = + Ex E 21 22 Và từ điều kiện () () f f l t Ft = x += + é ù ê ú ë û ( ) (0) ( ) (0) ( ) (0) ( ) (0) tt tt E x E FE x E l l 21 22 11 12 ff ff - = - é ùé ù ê úê ú ë ûë û ( ) ( ) (0) ( ) ( ) (0) t t tt E FE FE E x l 22 12 11 21 f f ff 1 t t tt - = - - é ùé ù ê úê ú ë ûë û (0) ( ) ( ) ( ) ( ) (0) E FE FE E x l 22 12 11 21 f f ff Nguyen Quang Hoang 18 Department of Applied Mechanics 9 12/2/2015 Điều khiển tối ưu LQR hệ tuyến tính Sau khi xác định được λ(0), ta sẽ có λ(t) và u(t). Tuy nhiên, bài toán cần phải giải lại khi điều kiện đầu thay đổi. Ngoài ra, đầu vào điều khiện u(t) là khả thi chỉ trong biểu diễn, vì điều khiển ở dạng này không phải là hàm của các biến trạng thái (vòng điều khiển mở). Lời giải bài toán tối ưu vòng điều khiển kín (phản hồi trạng thái): Sử dụng phép biến đổi sau ( ) ( ) ( ), T l t tt = = Sx S S (28) Thay vào (25) ta được dt= + =- - Sx Sx Qx A Sx  l (29) ( ) , d t T Cùng với (21, 24) ta được ( ) 1 0 - T T S SA SBR B S Q A S x  + - ++ = (30) Nguyen Quang Hoang 19 Department of Applied Mechanics Điều khiển tối ưu LQR hệ tuyến tính Biểu thức (30) đúng với mọi x(t) nên 1 0 - T T S SA SBR B S Q A S  + - ++ = (31) Từ (13) và (28) cho ta æ ö ç¶ ÷ G t t t tt t T = = = = ç ÷ l l (32) ç ÷ ç ¶ ÷ è øFx S x S F ( ) ( ), ( ) ( ) ( ) ( ) f f f ff f x = t t f Phương trình (31) được gọi là phương trình Riccati. Phương trình vi phân phi tuyến này có thể được giải quyết bằng phương pháp số để nhận được ma trận S(t). Từ đó xác định được véc tơ điều khiển - 1 T u R BS x K x =- =- ( ) ( ) ( ) ( ), tt tt (33) - 1 T K R BS () () t t = Cấu trúc của (33) có dạng điều khiển phản hồi trạng thái với ma trận tối ưu K(t). Vì lời giải S(t) không phụ thuộc vào trạng thái hệ, nên điều khiển này là tối ưu cho tất cả các điều kiện ban đầu của véc tơ trạng thái, x(0). Nguyen Quang Hoang 20 Department of Applied Mechanics 10 12/2/2015 Điều khiển tối ưu LQR hệ tuyến tính Cho hệ tuyến tính Thuật giải x Ax Bu A B  = + , (,) là điều khiển được 1. Chọn các ma trận Q = QT ≥ 0, R = RT > 0 và F = FT ≥ 0 2. Giải phương trình Riccati (31), hệ ptvp thỏa mãn điều kiện cuối 1 0, ( ) T Tf t - S SA SBR B S Q A S S F  + - ++ = = 3. Xác định ma trận K(t) từ đó xác định véc tơ điều khiển 1 ( ) ( ), ( ) ( ) T t t tt - K R BS u K x = =- 4. Đánh giá (xác định) giá trị hàm mục tiêu: 1 1 () () 2 2ft T TT f f J t t dt = ++ é ù ê ú ò ë û x Fx x Qx u Ru 0 Nhận xét: gặp khó khăn khi giải phương trình Riccati. Hàm mục tiêu đôi khi không có ý nghĩa vật lý, nhưng với các Q và R khác nhau sẽ cho đáp ứng khác nhau. Nguyen Quang Hoang 21 Department of Applied Mechanics Điều khiển tối ưu LQR hệ tuyến tính Trường hợp quan trọng: tf → ∞ Xét trường hợp thời gian cuối tf tiến đến vô cùng (thời gian điều khiển tương đối lớn). Nghiệm S(t) hay K(t) tiến đến ma trận hằng. [xem tài liệu tk]. Trong trường hợp này, ta sẽ giải phương trình đại số Riccati (Algebraic Riccati Equation-ARE) thay vì giải phương trình vi phân Riccati. 1. Chọn các ma trận Q = QT ≥ 0, R = RT > 0 và F = FT ≥ 0 2. Xác định ma trận S từ phương trình đại số Riccati 1 0, 0, T T do - SA SBR B S Q A S S - ++ = = 3. Xác định ma trận K(t) từ đó xác định véc tơ điều khiển 1 , () T t - K R B S u Kx = =- 4. Kiểm tra tính ổn định của hệ kín x A BK x  = - ( ) 5. Đánh giá (xác định) giá trị hàm mục tiêu: 1 1 () () 2 2ft T TT f f J t t dt = ++ é ù ê ú ò ë û x Fx x Qx u Ru 0 Lệnh Matlab: >> K = lqr(A,B,Q,R) Nguyen Quang Hoang 22 Department of Applied Mechanics 11 12/2/2015 Điều khiển tối ưu LQR hệ tuyến tính Xét trường hợp có số hạng chéo trong hàm mục tiêu 1 1 () () 2 2 2ft T TT T f f J t t dt = + ++ é ù ê ú ò ë û x Fx x Qx u Ru x Nu 0 Khi đó sử dụng phép biến đổi T TT T T x Qx x Nu u Ru x Q x v Rv + += + 2 Với m - - T T 1 1 Q Q NR N v u R N x =- =+ m Khi đó hàm mục tiêu trở thành & 1 1 () () 2 2ft T TT ff m J t t dt = ++ é ù ê ú ò ë û x Fx x Q x v Rv 0 Hệ ban đầu được biến đổi tương ứng thành - x A x Bv A A BR N  = + =- 1 , with , T m m Nguyen Quang Hoang 23 Department of Applied Mechanics Điều khiển tối ưu LQR hệ tuyến tính Giả sử rằng Qm đối xứng, bán xác định dương, khi đó bài toán được biến đổi về dạng như đã trình bày. Lời giải của bài toán này có nghiệm như sau: Thuật giải 1. Chọn các ma trận Q = QT ≥ 0, R = RT > 0, F = FT ≥ 0, và ma trận N 2. Tính các ma trận Qm và Am (kiểm tra ma trận Qm ≥ 0, bán x.đ.d) - - Q Q NR N A A BR N =- =- 1 1 , T T m m 3. Giải phương trình Riccati (31), hệ ptvp thỏa mãn điều kiện cuối m m m m m m mm m f t - S S A S BR B S Q A S S F  + - ++ = = 1 0, ( ) T T 4. Xác định ma trận Km(t) từ đó xác định véc tơ điều khiển = + - é ù ê ú ë û = - 1 () () , T T K R BS N t t m m u Kx ()() Lệnh Matlab: m t t >> K = lqr(A,B,Q,R) Nguyen Quang Hoang 24 Department of Applied Mechanics 12 12/2/2015 Ví dụ 1 Xét hệ động lực bậc nhất và hàm mục tiêu như sau 5 1 ( ) (2 3 ) 2 2ft x xu  = + , 2 22 f J x t x u dt =+ + ò 0 VớiAbF Q R == = = = 1, 1, 5, 2, 3 Phương trình Riccati trở thành 1 2 2S 2, ( ) 5 3 f S S St F  =- + - = = Với tf cho trước, giải phương trình vi phân trên ta nhận được S(t), 0 ≤ t ≤ tf Luật điều khiển tối ưu sẽ là 1 ( ) ( ) ( ), 0 3 f ut Stxt t t =- £ £ Nguyen Quang Hoang 25 Department of Applied Mechanics Ví dụ 1 Nếu cho tf →∞ 5 1 , ( ) (2 3 ) 2 2 f x x u J x t x u dt ¥ 2 22 =+ = + + ò  0 VớiAbF Q R == = = = 1, 1, 5, 2, 3 Phương trình đại số Riccati trở thành 1 0 2S 2 3 15 2 =- + -  =  S S 3 Để hệ kín ổn định ta chọn 1,2 1 S S = =+3 15 3 15 Luật điều khiển tối ưu sẽ là u Kx = - , K + = 15 3 x x  = - 3 Nguyen Quang Hoang 26 Department of Applied Mechanics 13 12/2/2015 Ví dụ 2 Xét hệ tuyến tính sau: é ù é ùé ù é ù é ù ê ú ê úê ú ê ú ê ú = += é ù ê ú ë û ë û ë ûë û ë û ë û  xx x 01 0, 10 00 1 11 1  u y xx x 22 2 1 1 22 u k x r kx =- - - ( ) Tìm luật điều khiển để duy trì 1 2 x rx = = , 0 cực tiểu chỉ tiêu Đặt biến mới  x xr = - J x r u dt ¥ 1 (( ) ) 2 2 2 = -+ ò 1 0    é ù é ùé ù é ù ê ú ê úê ú ê ú = + x xu (tức cho tf →∞) 1 1 01 0 , 00 1 J x u dt ¥ 1 ( ) 2    2 2  x x = 2 2 1 1 x x 2 2 ë û ë ûë û ë û = + ò  1 0 é ù 2 0; 2 0 0 = = ê ú ê ú ê ú ë û Q R Phương trình đại số Riccati trở thành 1 0 - T T SA SBR B S Q A S - ++ = Nguyen Quang Hoang 27 Department of Applied Mechanics Ví dụ 2 Giải phương trình đại số Riccati é ùé ù é ùé ù é ùé ù é ù é ù é ù ê úê ú ê úê ú ê úê ú ê ú ê ú é ù +- + ê ú ê ú ê ú ë û ë û ë ûë û ë ûë û ë ûë û ë û ë û 00 01 0 20 1 0 1 SS SS SS SS 11 12 11 12 11 12 11 12 10 00 1 00 2 SS SS SS SS 12 22 12 22 12 22 12 22 é ù 0 0 = ê ú ê ú ê ú ë û 0 0 2 2 S SS S 12 12 22 22 - += - = - + = S S 11 12 2 0, 0, 2 0 2 22 Kết quả nhận được ma trận S và luật điều khiển phản hồi trạng thái như sau S é ù ê ú = ê ú ê ú ë û 22 2 2 22 2 ( )2 T u x x xr x - 1 R B Sx  =- =- - =- - - 121 2 Nguyen Quang Hoang 28 Department of Applied Mechanics 14 12/2/2015 Ví dụ 2 - Matlab Xét hệ tuyến tính sau: x Ax B  = + u , J x u dt x x r ¥ 1 ( ), 2 2 2 = + =- ò   1 11 0 Hệ và bài toán tối ưu mô tả bởi: é ù éù é ù 01 0 20 ==== ê ú êú ê ú , , ;2 00 1 00 ABQR ë û ëû ë û Phương trình đại số Riccati 1 0 - T T SA SBR B S Q A S - ++ = Lệnh lqr của matlab >> [K, S, e] = lqr(A,B,Q,R) Nguyen Quang Hoang 29 Department of Applied Mechanics Ví dụ 2 - Matlab Lệnh lqr của matlab A = [0 1; 0 0]; B = [0; 1]; C = [1 0]; Q = [2 0; 0 0]; R = [2]; [K, S, e] = lqr (A, B, Q, R) closed_Syst = ss(A - B*K, B*K(1,1), C, 0) t = 0:0.1:10; r = 2*ones (size(t)); [y, t, x] = lsim (closed_Syst, r, t, [5 0]); u = -K*x' + K(1,1)*r; plot (t, y, 'k', 'linewidth', 2), grid on, hold on plot(t,u, 'r', 'linewidth', 2) legend('x_1(t)', 'u(t)') K = 1.0000 1.4142 S = 2.8284 2.0000 2.0000 2.8284 e = -0.7071 + 0.7071i -0.7071 - 0.7071i 5 x1(t) u(t) 4 3 2 1 0 -1 -2 e là trị riêng của ma trận (A - BK) 0 2 4 6 8 10 -3 Nguyen Quang Hoang 30 Department of Applied Mechanics 15 12/2/2015 Bài tập Bài 1. Xét hệ tuyến tính sau cùng với hàm mục tiêu cần được cực tiểu hóa: é ù é ùé ù é ù ê ú ê úê ú ê ú = + - - ë û ë ûë û ë û  x xu 01 0 , 12 1 1 1  x x 2 2 5 2 22 2 1 2 20 (5) ( 2 ) , 5 2 2 f J x x x x x u dt t =´ + + + + = ò 1 1 2 12 0 o Thiết lập các phương trình vi phân Riccati với điều kiện biên thích hợp; o Giải các phương trình vi phân Riccati bằng các lệnh ode23 hoặc ode45 của Matlab; o Tìm ma trận S(t), ma trận K(t) và đáp ứng vòng kín. Ngoài ra, tìm giá trị tối ưu của J. Nguyen Quang Hoang 31 Department of Applied Mechanics Tối ưu năng lượng tiêu thụ động cơ DC Bài 2. Thiết lập bài toán điều khiển vị trí cho động cơ điện một chiều, với chỉ tiêu cực tiểu hóa năng lượng. Hãy đưa vào hàm cần tối ưu, từ đó tìm hàm tác động (điều khiển) u(t), 0 ≤ t ≤ tf, để tiêu hao năng lượng trong khoảng thời gian xác định (cho biết tf) bởi động cơ DC là nhỏ nhất. Biết điều kiện đầu và điều kiện cuối như sau Nguyen Quang Hoang 32 Department of Applied Mechanics 16 12/2/2015 Tối ưu năng lượng tiêu thụ động cơ DC Năng lượng tiêu thụ của động cơ DC trong một thời gian cố định (từ 0 đến tf) được cho bởi tích phân sau: Để thỏa mãn các yêu cầu đối với sự tồn tại nghiệm của bài toán điều khiển LQR, hàm mục tiêu được sửa đổi để được Nguyen Quang Hoang 33 Department of Applied Mechanics Tối ưu năng lượng tiêu thụ động cơ DC Hàm mục tiêu cần tối thiểu Đưa biểu thức trong dấu tích phân về dạng của bài toán điều khiển LQR 2 T TT T T x Qx x Nu u Ru x Q x v Rv + += + m ta có được (xem trường hợp có số hạng chéo, x’Nu) - 1 T Q Q NR N m T = - - 1 v u R Nx = + Nguyen Quang Hoang 34 Department of Applied Mechanics 17 12/2/2015 Tối ưu năng lượng tiêu thụ động cơ DC Để ma trận Qm ≥ 0, các tham số ρ>0 và γ>0 cần được chọn sao cho Để có điều kiện trạng thái cuối bằng 0, ta sử dụng phép biến đổi Khi đó phương trình trạng thái với biến mới sẽ là Do hàm mục tiêu không chứa các biến x1 và x2 = x2~, hàm mục tiêu được viết lại thành Nguyen Quang Hoang 35 Department of Applied Mechanics Tối ưu năng lượng tiêu thụ động cơ DC Trong điều khiển LQ, không có ràng buộc về điều kiện cuối cùng (trạng thái cuối). Do đó, hàm mục tiêu được chỉnh sửa để chứa các vector trạng thái cuối: Với ma trận Sf Bỏ qua ràng buộc về trạng thái cuối, bài toán LQ có thể được giải để tối thiểu hàm mục tiêu của hệ tuyến tính với điều kiện đầu khác không. Bằng cách chọn các số a và b đủ lớn, sẽ giảm ảnh hưởng của trạng thái cuối đến lời giải [giảm vai trò của của x~(tf)]. Nguyen Quang Hoang 36 Department of Applied Mechanics 18 02/12/2015 Phân tích hệ động lực Hệ LTI Ổn định của hệ tuyến tính bất biến theo thời gian Nguyen Q.Hoang Department of Applied Mechanics Hanoi University of Science and Technology Nguyen Quang Hoang 1 Department of Applied Mechanics Tổng quan về tính ổn định Về mặt tổng quát, rất khó có thể nhận được nghiệm dạng hiển của phương trình vi phân phi tuyến. Tuy nhiên, thông thường ta không cần biết nghiệm dạng hiển của hệ phương trình mô tả hệ. Nhưng chúng ta cần quan tâm đáp ứng của hệ khi thời gian tiến về vô cùng [sau khoảng thời gian đủ dài]. Hơn nữa, trong nhiều mô hình của hệ thống vật lý một sự thay đổi "nhỏ" các điều kiện đầu dẫn đến một sự thay đổi "nhỏ" của các nghiệm. Đặc biệt, nếu một hệ thống bị nhiễu từ trạng thái nghỉ của nó, hay vị trí cân bằng, sau đó nó bắt đầu di chuyển. Chúng ta gần như có thể nói rằng vị trí cân bằng là ổn định nếu hệ thống không đi xa khỏi vị trí này đối với nhiễu loạn nhỏ ban đầu. Tuy nhiên, có những hệ mà sự thay đổi nhỏ của các điều kiện đầu dẫn đến sự thay đổi lớn trong ứng xử của hệ. Nguyen Quang Hoang 2 Department of Applied Mechanics 1