9/16/2015

Chương 1. Môi trường Matlab

1. Matlab là gì ?

Lecture on Matlab

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

Bài mở đầu

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

Nguyen Q.Hoang

Department of Applied Mechanics

8. Kết thúc phiên làm việc với Matlab

Hanoi University of Science and Technology

1

2

Chương 1. Môi trường Matlab

Chương 1. Môi trường Matlab

1. Matlab là gì ?

• Matlab có thể chạy trên hầu hết các hệ máy tính: máy tính sách

• MATLAB được viết tắt từ “MATrix LABoratory”, là một công cụ tính tay - Labtop, máy tính cá nhân PC, đến các hệ thống siêu máy toán số và mô phỏng số, công cụ này được phát triển dựa trên các tính (super computer). Matlab được điều khiển bởi các tập lệnh, thư viện hàm tính toán số viết bằng ngôn ngữ lập trình FORTRAN.

tác động qua bàn phím trên cửa sổ điều khiển. Nó cũng cho phép Matlab có giao diện thân thiện với người sử dụng.

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à

• Một định nghĩa khác: Matlab là một phần mềm đa năng tính toán số, ngày càng được mở rộng nhờ các thư viện trợ giúp hay do người thể hiện các số liệu bằng hình ảnh trực quan và là một ngôn ngữ đa sử dụng tạo ra.

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

4

1

9/16/2015

Giao diện người sử dụng

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 1. Cửa sổ lệnh (Command

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 4

Window).

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

hiện ta xét ví dụ sau. Nếu trong cửa sổ lệnh ta viết 5

đã thực hiện (Command

>> x + 18 = 120

History).

•

Matlab sẽ báo lỗi với dòng hiển thị mầu đỏ:

3

3. Cửa sổ cho biết thư mục

??? x+18=120

hiện thời (Current

Directory) và danh sách

Error: The expression to the left of the equals sign is not a valid target for an 1

các biến đang sử dụng

assignment.

(Workspace).

2

•

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ó 4. Thanh chỉ đường dẫn

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 đó.

của thư mục hiện thời, cho

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 phép chọn thư mục hiện

6

>> x=120 - 18;

thời.

•

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ị

5. Thanh Shortcut.

của x hiện ra ở dòng dưới.

6. Nút Start.

>> x=120 - 18

5

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

>> a=3.5; b=7.5; % gan gia tri so cho

>> a/b % phep chia hai so, a:b (=

phép tính số học (cộng +, trừ

hai bien a va b

phep chia phai)

-, 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

>> a+b % phep cong hai so

ans =

dưới đây.

ans =

0.4667

11

>> a\b % phep chia trai (hieu la b

Phép tính

Cách viết toán học

Cách viết trong

>> a-b % phep tru hai so

chia cho a, b/a)

Matlab

ans =

ans =

phép cộng

c = a + b

c = a + b

-4

2.1429

>> a*b % phep nhan hai so

>> a^b % phep luy thua, a mu b

phép trừ

c = a – b

c = a – b

ans =

ans =

phép nhân

c = ab

c = a* b

26.2500

1.2037e+004

phép chia

c = a : b

c = a / b

>> a/b

% phep chia hai so, a:b (= phep chia

phép chia phải

c = a : b

c = a / b

phai)

phép chia trái

c = b : a

c = a \ b

ans =

lũy thừa

c = ab

c = a^ b

0.4667

7

8

2

9/16/2015

Một số dạng format

Các định nghĩa và hàm toán học cơ bản

•

format long, short, rat, short e, bank

•

Để 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ụ

nhiều đại lượng toán học và các hàm cơ bản. Ví dụ như số pi đã được định

>> format long

>> format long e

nghĩa sẵn

>> co = cos(0.2),

>> x=5.125*3.16

>> R = 2;

si = sin(0.2), tn = tan(0.2)

x = 1.619500000000000e+001

>> V = 4/3*pi*R^3

co = 0.98006657784124

si = 0.19866933079506

>> format rat

V = 33.5103

tn = 0.20271003550867

>> x=5.125*3.16

>> exp(1)

>> log(3.2) % ln(3.2)

x = 3239/200

ans = 2.7183

ans = 1.1632

>> format short

>> format bank

>> co = cos(0.2), si = sin(0.2),

>> luonggio=35.55

>> exp(2)

>> x = 5; log(x)

tn = tan(0.2)

luonggio = 35.55

ans = 7.3891

ans = 1.6094

co = 0.9801 si = 0.1987 tn = 0.2027

>> luongtuan=luonggio*40

>> x = sqrt(9)

>> x = 3; log10(x)

luongtuan = 1422.00

x = 3

ans = 0.4771

>> format short e

>> y = sqrt(11)

>> x=5.125*3.16

y = 3.3166

x = 1.6195e+001

9

10

Các định nghĩa và hàm toán học cơ bản

Các định nghĩa và hàm toán học cơ bản

Các hàm lượng giác - Trigonometric functions

Các hàm toán học

sin

Hàm sin

>> sin(pi/3) cho 0.8660

sqrt(x)

Tính căn của biến x

sind

Hàm sin với đối số là độ

>> sind(60) cho 0.8660

exp(x)

Hàm e mũ, Exponential

sinh

Hàm sin Hyperbolic

log(x)

Logarithm cơ số tự nhiên, cơ số e

asin

Hàm sin ngược, tức arcsin

>> asin(0.866) cho 1.0471

log1p(x)

Tính giá trị log(1+x)

asind

Hàm sin ngược cho kết quả là độ

>> asind(0.866) cho 59.9971

log10(x)

Logarithm cơ số 10

asinh

Hàm sin Hyperbolic ngược

log2(x)

Logarithm cơ số 2

cos

Hàm cos

pow2(x)

Hàm mũ cơ số 2, pow2(x) = 2^x

cosd

Hàm cos với đối số là độ



realpow(x,y)

Tính x^y, với x, y là thực

x

x

e  e

cosh

Hàm cos Hyperbolic

sinh( x) 

reallog(x)

Logarithm cơ số tự nhiên của số thực

2

nthroot(x, n)

Căn bậc n của số thực

acos

Hàm cos ngược, tức arccos

x

 x

nextpow2(n)

Trả lại số p đầu tiên sao cho 2^p >= abs(n)

acosd

Hàm cos ngược cho kết quả là độ

e  e

acosh

Hàm cos Hyperbolic ngược

cosh( x) 

2

11

12

3

9/16/2015

Số phức

Số phức

Chúng ta cũng có thể đưa số phức vào trong Matlab.

Nhớ lại rằng đơn vị phức được định nghĩa là căn của im

•

Một số hàm liên quan đến số phức được liệt kê trong bảng dưới đây (-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

Các hàm với số phức

trong đó x là phần thực và y là phần ảo của z

abs

Cho độ lớn hay modul véctơ phức

•

Ví dụ

re

angle

Cho góc pha hay argument của số phức, radian

a = 3 + 3i,

b = 1 – i  a + b = 4 + 2i

complex

Tạo lập dữ liệu phức từ các phần thực vào ảo

2

2

| a | abs( a)

r (

e a)

im( a)

conj

Cho số phức liên hợp

arg( a)

arctan( im( a) / r (

e a))

imag

Cho phần ảo

>> a = 3 + 3i, b = 1 - i;

a | a | i

e ,

i

e

cos

i sin

real

Cho phần thực

>> c = a + b

isreal

Trả lại giá trị đúng, 1, nếu đối số là thực

c = 4.0000 + 2.0000i

>> abs(a)  4.2426 = 2 sqrt(3) cplxpair

Sắp xếp các số phức thành cặp liên hợp

>> angle(a)  0.7854 = pi/4

13

14

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

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ơ

Lecture 2

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

Véctơ, ma trận và Matlab

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

Nguyen Q.Hoang

Tạo véctơ hàng có các phần tử cách đều

Department of Applied Mechanics

Hanoi University of Science and Technology

Các đặc trưng của véctơ

hoang.nguyenquang@hust.edu.vn

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ơ

Tạo một véctơ từ các véctơ

Tạo véctơ hàng có các phần tử

2

>> a= [2; 1; 4] % vt cot

khác

cách đều (toán tử hai chấm : )

a

1

y

a '

2 1 4

>> v= [2 0 4] % vt hang

x = [a1 : h : a2]

>> w= [1,1,9] % sd dau phay

4

Ghép hai véc tơ cột

x(i) = a1 + (i - 1)h, i = 1, 2, 3, . . .

d =

• Chuyển vị vectơ (cột thành hàng,

>> a = [1; 4; 5];

hàng thành cộ

>> b = [2; 3; 3];

t)

1

>> d = [a; b]

>> x = [0:2:10]

>> a = [2; 1; 4];

1

2

4

x = 0 2 4 6 8 10

>> y = a'

a

4 ,

b

3

5

2

• Các phép tính cộng và trừ hai

5

3

Ghép hai véc tơ hàng

>> x = [0:0.1:1]

3

véctơ cùng cỡ

3

x = 0 0.10 0.20 0.30 0.40 0.50 0.60

>> a = [1 4 5];

>> a = [1; 4; 5];

3

1

0.70 0.80 0.90 1.00

>> b = [2 3 3];

>> b = [2; 3; 3];

c

a

b

7 , c

a

b

1

>> help linspace

>> d = [a b]

d = 1 4 5 2 3 3

>> c1 = a + b

1

2

8

2

>> c2 = a - b

1

9/16/2015

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ơ

• Cho các vectơ cùng cỡ v và u ta có các hàm tác động lên véc tơ

Ví dụ

Tên hàm

Ý nghĩa

>> a = [8 4 4 1 7 11 2 0];

length(v)

Cho biết chiều dài, số phần tử của véctơ v

>> a = [2;3;3;4;5];

>> sum(a)

>> length(a)

ans = 37

sum(v)

Trả lại tổng các phần tử của véctơ v

ans = 5

>> prod(a)

prod(v)

Trả lại tích các phần tử của véctơ v

>> b = [1, 1, 2];

ans = 0

min(v)

Trả lại phần tử nhỏ nhất của véctơ v

>> length(b)

max(v)

Trả lại phần tử lớn nhất của véctơ v

ans = 3

>> u = [1; 4; 5]; v = [2; 3; 3];

sort(v)

Xếp các phần tử tăng dần

>> x = dot(u,v)

>> a = [8 4 4 1 7 11 2 0];

find(v)

Tìm các phần tử khác không trong v

x = 29

>> max(a)

dot(u,v)

Trả lại tích vô hướng hai véctơ

ans = 11

>> w = cross(u,v)

cross(u,v)

Trả lại tích có hướng hai véctơ cỡ 3*1

>> min(a)

w =

ans = 0

-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ử :

>> v = a(1:3)

v =

2.2 Biểu diễn đa thức và các phép tính đa thức

>> a = [12; 17; –2; 0; 4; 27];

12

Nhập đa thức

>> a(2)

17

ans = 17

–2

Các phép tính trên đa thức

>> a(6)

Phép nhân đa thức / Phép chia đa thức

ans = 27

>> a(:)

Phép cộng và trừ đa thức

ans =

Không điểm hay nghiệm của phương trình đa thức 12

17

Xây dựng đa thức từ các không điểm cho trước

–2

0

Giá trị của đa thức tại một điểm

4

Đạo hàm đa thức

27

2

9/16/2015

Biểu diễn đa thức và các phép tính đa thức Biểu diễn đa thức và các phép tính đa thức

• Nhập đ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 Trong Matlab một đa thức được đưa vào dưới dạng một véctơ hàng có các

>> n = length(p)-1

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 n = 5

cần biểu diễn dạng đầy đủ của nó, ví dụ xét đa thức

•

Các hàm liên quan đến đa thức

5

4

2

p ( s)  a s  a s  a s  a s  a 5

5

4

2

1

0

p = conv(p1,p2)

Nhân hai đa thức

với dạng đầy đủ như sau

[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

5

4

3

2

1

0

p ( s)  a s  a s  0 s  a s  a s  a s roots(p)

Tìm các nghiệm của phương trình đa thức, p(x) = 0

5

5

4

2

1

0

p = poly(nu)

Xây dựng đa thức từ các không điểm cho trước, đưa ra 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ử

cho ta một đa thức nhận các giá trị trong véctơ nu làm nghiệm.

p  [ a a 0 a a a ]

5

5

4

2

1

0

polyval(p,s0)

Giá trị của đa thức tại một điểm

Ví dụ

polyder(p)

Đạo hàm đa thức

>> p = [8 3 0 -4 6 0]

Phép cộng và trừ đa thức ??

5

4

2

p ( s)  8 s  3 s  4 s  6 s 5

p = 8 3 0 -4 6 0

Biểu diễn đa thức và các phép tính đa thức 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

>> nu = [-2 -1 0 1 2]

>> poly(nu)

[ q, r ] = deconv(p1, p2)

ans =

5

3

p ( s)  s  5 s  4 s 5

4

2

 s( s  5 s  4)

p ( s)

r( s)

1 0 -5 0 4 0

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

0

k

1

3

9/16/2015

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 Ma trận và các phép tính cơ bản trên ma trận

•

Nhập 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 Nhân ma trận với một số, phép cộng và trừ hai ma trận ma trận trong Matlab, chúng ta đặt từng hàng vào trong dấu ngoặc vuông [], Chuyển vị ma trận

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.

Phép nhân ma trận

•

Ví dụ với hai ma trận

Một số phép tính cơ bản khác





Các ma trận đặc biệt

2

0 1

 1



6 





Tham chiếu đến các phần tử của ma trận

A  

,

B   1



7 4

7

11









Ma trận khối

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

•

Ta đưa vào Matlab như sau:

>> B = [2,0,1; -1,7,4; 3,0,1]

Tìm hạng của ma trận

>> A = [-1, 6; 7, 11]

B =

Tìm ma trận nghịch đảo và ma trận tựa nghịch đảo (pseudo-inverse) A =

2 0 1

Trị riêng và véctơ riêng của ma trận vuông

-1 6

-1 7 4

Phân tích ma trận vuông A thành tích các ma trận 7 11

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ố

Chuyển vị ma trận

Phép nhân hai ma trận có cỡ thích hợp, C = A*B

>> A = [-3, 6; 7, 10] ;

1

2

 

A  { a },

i  1.. ,

m j  1.. ,

p

>> B = [5, 1; -1, 3] ;





 1 3 7

ij

A  3

5  ,

T

A 

>> C=2*A





B  { b },

i  1.. ,

p j  1.. n

2



5 4









ij

C = -6 12

7

4











p

C = AB 

 





14 20

{ c },

c

a b , i

1.. ,

m j

1.. ,

n

ij

ij

ik kj

Sử dụng dấu nháy đơn (’)

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 =

>> C = A*B

>> B = [-1 7 4; 2 1 -2];

ans =

–1 2 0

C =

>> C = A*B

2 7

6 4 1

11 14

C =

6 13

>> B = A'

>> A-B

B =

13 16

7 11 -4

–1 6

ans =

-8 56 32

2 4

-8 5

7 -4 -10

0 1

8 7

4

9/16/2015

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 Phép nhân phần tử hai ma trận

•

Các ma trận đặc biệt

Khởi tạo các ma trận đặc biệt

C = .

A * B  { c },

c  a * b ,

i  1.. ,

m j  1.. n

ij

ij

ij

ij

eye(m,m)

tạo ma trận đơn vị cỡ mxm

Cần phân biệt phép nhân phần tử ( .* ) hai ma trận cũng cỡ C = A .* B

eye(m)

tạo ma trận đơn vị cỡ mxm

eye(m,n)

tạo ma trận đơn vị mở rộng cỡ mxn

>> A = [2 1; 1 2];

>> A = [1 4; 8 0; -1 3];

zeros(m,n)

tạo ma trận không cỡ mxn

>> B = [3 4; 5 6];

B = [-1 7 4; 2 1 -2];

>> C = A .* B

>> C = A * B

zeros(m,m)

tạo ma trận không cỡ mxm, hay ma trận vuông

zeros(m)

C =

C =

6 4

7 11 -4

ones(m,n)

cho ta một ma trận 1 cỡ mxn

5 12

-8 56 32

ones(m,m)

cho ta một ma trận 1 cỡ mxm, hay ma trận

7 -4 -10

ones(m)

vuông

>> C = A .* B

magic(n)

tạo ma trận magic vuông cỡ nxn, n > 2

??? Error using ==> times

[ma phương]

Matrix dimensions must agree.

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 Các ma trận đặc biệt

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

>> eye(3)

>> zeros(3,2)

>> magic(3)

chéo từ một vector bằng lệnh diag

ans =

ans =

ans =

•

c = diag(A) cho ta đường chéo của ma trận A

1 0 0

0 0

8 1 6

•

A = diag(c) cho ta ma trận đường chéo A, với A(i,i) = c(i) 0 1 0

0 0

3 5 7

0 0 1

0 0

4 9 2

>> eye(2,3)

>> zeros(2)

>> A1=[1 2 3; 3 4 5; 3 5 7]

>> Add=diag(Ad)

ans =

ans =

>> magic(5)

A1 =

Add =

1 0 0

0 0

ans =

1 2 3

1 0 0

0 1 0

0 0

17 24 1 8 15

3 4 5

0 4 0

>> ones(2,2) or ones(2)

23 5 7 14 16

3 5 7

0 0 7

>> zeros(3,3)

4 6 13 20 22

>> Ad=diag(A1)

ans =

ans =

10 12 19 21 3

Ad =

0 0 0

1 1

11 18 25 2 9

1

0 0 0

1 1

4

0 0 0

7

5

9/16/2015

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

• Tham chiếu đến các phần tử của 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

•

Cho các ma trận có kích cỡ thích hợp: A, B, C, M

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

A B 

số

•

Tạo ma trận cỡ lớn hơn từ các ma trận trên

G  



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

C M

(1,2,...), Matlab không sử dụng chỉ số 0 như một số ngôn ngữ lập trình khác.

•

G = [ A, B ; C , M ];





Ví dụ đối với ma trận

•

Ví dụ, với A và B









>> A = [2 4; 7 5]

>> A = [1 2 3;

>> r3=A(3,:)

2 4

4 2

A  

,

B  



>> B = [4 2; 5 7]

4 5 6;

r3 =

7 5

5 7









>> Q = [ zeros(2) eye(2); -A\B B*A’ ]

7 8 9];

7 8 9

Q =

>> A(:,2:3)

•

Hãy xây dựng ma trận Q sau

0 0 1 0

>> A(2,3)

>> c2 = A(:,2)

ans =

 0

E

2 2



2 2

 

0 0 0 1

ans = 6

ans =

2 3

Q  



1



T

0 -1 16 38

2

5 6

A B BA 





-1 0 38 70

5

8 9

8

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 thức và giải hệ phương trình đại số tuyến tính Đị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

a

a



a

a

11

12 

11

12

>> A = [5 2 –9; –9 –3 2; 6 7 3];

A 

,

a



a 

det( )

A 

 a a  a a

11 22

12 21

a

a

>> b = [44;11;5];

21

22





21

22

>> A = [1 3; 4 5];

>> B = [3 –1 2 4; 0 2 1 8; –9 17 11 3; 1 2 3 –3];

>> det(A)

>> det(A)

>> det(B)

ans = 368

ans = –7

ans = –533

•

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

5 x  2 y  9 z  44

A   9



2



2 

x  y ,

 

b  11





 

 

7.6902

9

 x  2 y  2 z  11

6

7

3





z





44

 

 

 

–6.0272

6 x  7 y  3 z  44

•

Nếu det(A)0 pt Ax = b có nghiệm duy nhất

x = A \ b

6

9/16/2015

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 Tìm ma trận nghịch đảo

• Để kiểm tra lại ta tính tích hai ma trận:

•

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

>> iS*S

>> S*iS

A A-1 = A-1A = E, với E là ma trận đơn vị cùng cỡ.

ans =

ans =

•

Điều kiện để tồn tại ma trận nghịch đảo là det(A)  0.

1.0000 0 0 0

1.0000 0 0 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 0 1.0000 0 0.0000

0.0000 1.0000 -0.0000 0

A x = b có nghiệm duy nhất là x = A-1 b 0 0 1.0000 0

0.0000 -0.0000 1.0000 0

•

Tính nghịch đảo ma trận vuông bằng lệnh inv(A) 0 0 0.0000 1.0000

0 0 0 1.0000

>> A = [2 3; 4 5]

>> S = [1 0 -1 2; 4 -2 -3 1; 0 2 -1 1; 0 0 9 8];

• Ví dụ giải hệ A.x = b

>> det(A)

>> det(S)

>> A = [3 –2; 6 –2]; b = [5; 2];

>> e =A*x – b % thu lai

ans = -2

ans = -108

>> iS=inv(S)

>> det(A)

e =

>> iA = inv(A)

iS =

ans = 6

1.0e-015 *

iA =

-0.9259 0.4815 0.4815 0.1111

>> x=inv(A)*b

0.0000

-2.5000 1.5000

-0.6296 0.1574 0.6574 0.0556

x =

0.8882

2.0000 -1.0000

-0.5926 0.1481 0.1481 0.1111

-1.0000

0.6667 -0.1667 -0.1667 0

-4.0000

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 Trị riêng và véctơ riêng của ma trận vuông

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 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 l = eig(A);

cho một véctơ l chứa các trị riêng của ma trận vuông A.

Ax  x hay (A   )

E x  0 với E là ma trận đơn vị.

cho ta ma trận chéo D chứa các trị riêng của A, ma trận V có

[V, D] = eig(A);

•

các cột là các véctơ riêng của

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 A thỏa mãn AV = VD.

thoả mãn điều kiện

cho ta véctơ chứa các trị riêng suy rộng của các ma trận d = eig(A,B);

Ax  Bx hay ( A   )

B x  0

vuông A và B.

thì số 

cho ta ma trận chéo D chứa các trị riêng suy rộng của hai ma 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

[V, D] =eig(A,B); 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.

tương tự như eig(A,B) đối với ma trận đối xứng A và ma trận eig(A,B,'chol')

đố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

9/16/2015

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 Trị riêng và véctơ riêng của ma trận vuông

Phân tích LU và phân tích Cholesky

•

Ví dụ

•

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

>> C = [2 1 0;

>> [V, D]=eig(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 1 4 -1;

V =

Ý 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 0 -1 7];

0.9119 -0.4064 -0.0571

tích LU ta nhận được

>> eig(C)

-0.4037 -0.8630 -0.3037

Ax  LUx  b

ans =

-0.0742 -0.3000 0.9510

Đặt y  Ux  Ly  b

ta giải hai hệ Ly  b

và Ux  y

1.5573

D =

4.1233

1.5573 0 0

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.

7.3194

0 4.1233 0

0 0 7.3194

•

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

•

Ví dụ phân tích cholesky

Các phép tính ma trận, vector

>> A = [3 2 -9; -9 -5 2; 6 7 3];

>> M = [2 3 4;3 5 6;4 6 9];

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ử

b = [-65; 16; 5];

>> L = chol(M)

.* .\ ./ .

Chia trái và chia phải

>> [L, U, P] = lu(A)

L =

\ /

L =

transpose(A), A’, A.’

Chuyển vị (ma trận, vector)

1.4142 2.1213 2.8284

-0.3333 0.0909 1.0000

ctranspose(A), A’, A.’

Chuyển vị, ma trận, vector liên hợp

0 0.7071 0.0000

1.0000 0 0

Nghịch đảo ma trận

0 0 1.0000

inv(A)

-0.6667 1.0000 0

>> L‘ * L

Tính định thức ma trận vuông

- M

det(A)

U =

Giải hệ AX=B, phân tích LU

ans =

linsolve(A, B)

-9.0000 -5.0000 2.0000

eig(A)

Trị riêng, vector riêng ma trận

0 3.6667 4.3333

1.0e-015 *

rank(A)

Hạng ma trận

0 0 -8.7273

0.4441 0 0

Cho biết cỡ ma trận, vector

>> x = U\(L\b)

0 0 0

<[m, n]= > size(A <,i>)

(i=1 số hàng, i=2 số cột)

x = 2.0000

0 0 0

length(A)

Giá trị lớn nhất của số hàng và số cột

-4.0000

(=max(sizeA))

7.0000

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)

Tổng các phần tử của vector v

prod(v)

Tích các phần tử của vector v

min(v)

Phần tử nhỏ nhất của vector v

max(v)

Phần tử lớn nhất của vector v

max(A(:))

Phần tử lớn nhất của ma trận A

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

9

9/16/2015

Ma trận tựa nghịch đảo

Ma trận tựa nghịch đảo

Ma trận tựa nghịch đảo [pseudo-inverse], help pinv 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ổ

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 Giải hệ A x = b; A là ma trận chữ nhật

Ax  b

ban đầu về một hệ gồm m phương trình độc lập (m<n). Trường hợp này sẽ

-Số ẩn > số phương trình

được xem xét trong phần sau.

-Số ẩn < số phương trình

Ma trận A chữ nhật cỡ mxn, (m<n), hạng bằng m, [full rank]

Hệ phương trình đại số tuyến tính Ax = b sẽ có vô số nghiệm. Điều này cho

-Số ẩn > 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

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 x sao cho độ lớn (chuẩn) của x nhỏ nhất nghiệm.

Ax

,

m n

b A





 

,

n  ,

m

rank( )

A  m

Đưa vào tiêu chuẩn tối ưu cho nghiệm cần tìm

1

Tìm nghiệm của phương trình

T

J  x x  min

2

1

Ax

b

sao cho J

(x

x ) T

(

W x

x )

min

2

o

o

Bài toán xuất hiện :

W

-Ma trận trọng số (xác định dương)

-Robot dư dẫn động (số tọa độ khớp >6, robot kg; số tọa độ khớp >3, robot phẳng; [kinematics of redundant manipulator]

x0

-Véc tơ tham chiếu (được chọn trước)

-Tàu lặn mini (số cánh quạt >6).

1

T

x

0

J

x Wx

min

T

2

W

2E

J

x x

( x )

min

o

2

2

Ref. Mot so bai toan lien quan ma tran

Ma trận tựa nghịch đảo

Ma trận tựa nghịch đảo [pseudo-inverse]

Sử dụng PP nhân tử Lagrange để giải quyết bài toán tối ưu trên: 1

T

x

W A λ

x

1

T

1

λ

[AW A ] (b

Ax )

o

o

ta sẽ tìm cực trị của hàm Lagrange như sau

1

T

1

T

1

W A [AW A ] (b

Ax )

x

o

o

1

( ,

x )

λ

(x

x ) T

(

W x

x )

T

λ [b

Ax]

m 1

λ

2

o

o

1

T

1

T

1

1

T

1

T

1

W A [AW A ] b

(E

W A [AW A ]

)

A xo

Đạo hàm theo biến x ta nhận được

Ma trận †

1

T

1

T

1

A

W A [AW A ]

( ,

x λ)

(

W x

x )

T

A λ

0

W

o

x

được gọi là ma trận tựa nghịch đảo có trọng số của ma trận A.

Giải tìm được

1

T

x

W A λ

x

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.

o

Thay trở lại phương trình ban đầu

†

T

T

1

A

A [AA ]

†

†

1

T

x

A b

(E

A )

A x ;

b

Ax

AW A λ

Ax

o

o

Nếu A vuông, ko suy biến

Do A có hạng đầy đủ nên

†

khi x

0

x

A b

o

1

†

T

T

1

T

T

1

( )

A

det[

T

rank

m

AW A ]

0

A

A [AA ]

A [A A ]

T

T

1

1

1

T

1

A A A

A

λ

[AW A ] (b

Ax )

o

10

9/16/2015

Ma trận tựa nghịch đảo [pseudo-inverse]

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 Đạo hàm J theo biến x và cho bằng 0, ta nhận được có nghiệm thỏa mãn.

m n

n 1

m 1

Ax

,

b

A

, x

, b

, m

n

J

T

A

(

W Ax

)

b

0

T

1

[

]

T

x

A WA A Wb

x

rank( )

A

n

Khi W là một ma trận đơn vị

PT không có nghiệm thỏa mãn, Tuy nhiên, ta có thể tìm được vectơ x sao cho sai số

n

1

T

1

2

T

1

T

†

J

e e

e

and x

[A A] A b

A b

e

Ax

b

min

2

2

i

i 1

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

†

T

1

[

]

T

A

A A A

được gọi là ma trận tựa nghịch đảo trái của A hàm sau đây đạt cực tiểu

(xác định nghiệm hệ có số ẩn < số pt)

1

T

1 (

) T

J

e We

Ax

b

(

W Ax

)

b

2

2

W là ma trận trọng số đối xứng xác định dương

†

T

T

1

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]

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 ]

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 f ( )

x

(Ax

) T

b (Ax

)

T

b

x (

)

E x

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 T

x ( T

A A

)

T

T

T

T

E x

b Ax

x A b

b b

min

trận ở điều kiện yếu)

det( )

A

A

n n

n 1

Ax

,

b

A

, ,

x b

Đạo hàm theo biến véctơ x và cho kết quả bằng không, ta nhận được 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 f ( )

x

2( T

A A

)

E x

2 T

A b

0

phương tối thiểu có trọng số ( damped least square inverse). Theo phương x

pháp này, nghiệm của phương trình được tìm với điều kiện bổ sung như

sau:

Giải được

2

2

min( Ax

b

x )

x

Chọn trọng

T

1

(

)

T

x

A A

E

A b

†

T

1

0 :

[

]

T

x

A b

A A A b

số  > 0

f ( )

x

(Ax

) T

b (Ax

)

T

b

x (

)

E x

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 T

x ( T

A A

)

T

T

T

T

E x

b Ax

x A b

b b

min

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 Method based on the pseudo-inverse of jacobian matrix Problems

Given (

x t)  [r , O ]  m determined (

q )

n

t  





E

E

f



f



f



1

1

1



...



 f ( q) dq

 q



q



q

 

1

2

n

 f

nonlinear equations

J( q) 

  ...

...

...

... 

f( ,

x )

q  0

 x or

J( q) q  x

 q dt

 q





more unknowns than equations

 f



f



f

 

m

m

...

m

 q

q

q 







 1

2

n 

system of m linear equations, with n unknown Methods: analytical and numerical

Number of solutions

introduction the functional

1

T

J  q Wq

 min

2

W : positive weighting matrix

Solution

1



T

1



T

1

q W J ( q)[ J( q) W J ( q)]



x

t

q( t)  q(0)   d

q t

0

if W is identical matrix

T

T

1

q

J ( q)[ J( q) J ( q)]





x  J x



T

T

1

J

J ( q)[ J( q) J ( q)]



Finite number of solutions

Infinite number of solutions

Nguyen Q. Hoang, Nguyen V. Khang

45

Nguyen Q. Hoang, Nguyen V. Khang

46

Application of pseudo-inverse: Inverse kinematics Application of pseudo-inverse: Inverse kinematics Blockdiagram of Inverse kinematics

Method of Error feedback

e

K ,

e

e

x

f ( q)

e t

( )

0

q

e

x

J( q) q

Ke

†

J ( q)

W

1/ s

x

q( t)

J( q) q

x

Ke

q

J ( x

Ke)

q(0)

x(0)

Solving for q0

Null space of jacobian matrix J( q) Drawback of this diagram :

q

J ( x

Ke) ( E

J J) zo

accumulated errors due to rounding and integral method makes

x(

t), q( t)

z

don’t satisfy the constraint equations, so the end-effector moves out the By putting

o

the redundancy of manipulator can be exploited:

given trajectory.

- advoiding obstacles

- advoiding impact with the joint limitations

f ( x( t), q( t))  0

- advoiding singular configurations

Nguyen Q. Hoang, Nguyen V. Khang

47

Nguyen Q. Hoang, Nguyen V. Khang

48

12

9/16/2015

Application of pseudo-inverse: Inverse kinematics Application of pseudo-inverse: Numerical experiments Method of adjustment in error of joint variables A 5-DOF planar manipulator

Manipulator parameter

Link

1

2

3

4

5

q 3

B

C

y

m [kg]

10

7.5

5

4

1

q

q 4

x( t)

2

x



C

q

3

x( t)

I [kgm2]

q

J ( x

Ke)

A

C

z

0.11

0.10

0.03

0.02 0.02

4

1

C2

x( t)

s

Adjusting

( E

J J) z

D

o

*

q

q

l [m]

q

0.55

0.50

0.45

0.40 0.20

q

1

5

g

E

C1

a [m]

0.15

0.16

0.19

0.20 0.10

solving for

x

x(0)

O

o

qo

q

x

2

Determ. q

1.5

1

In this simulation, the end-effector (link 5)

Main idea of this method:

will be forced to move at the velocity of 2

0.5

find q closed to q*, such that q satisfy constraint equation f(x, q)=0

m/s along a cirlular trajectory, while its

orientation is constant, 1.0 rad. The

00

2

4

6

8

10

trajectory has a center at point (0.8, 0.5) m

t [s]

and radius of R = 0.5 m.

x [m]

x [m]

x [rad]

1

2

3

Nguyen Q. Hoang, Nguyen V. Khang

49

Nguyen Q. Hoang, Nguyen V. Khang

50

Application of pseudo-inverse:

Application of pseudo-inverse: Numerical experiments

A. Without consideration of joint limitations (planar manipulator) Numerical experiments / without post-adjusting

-15

x 10

3

1.5

30

20

2

1

q

J x

q

J ( x

Ke)

10

]

1

0.5

]

d

m

0

ar[

N

[

q

0

0

u -10

-20

-1

-0.5

-30

-2

-1

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

t [s]

t [s]

t [s]

q

q

q

q

q

e [m]

e [m]

e [rad]

u

u

u

u

u

1

2

3

4

5

1

2

3

1

2

3

4

5

B. Consideration of joint limitations (planar manipulator)

-15

q

J x

[ E

J ( q) J( q)] z

q

J ( x

Ke) [ E

J ( q) J( q)] z

x 10

0

0

2

30

1.5

20

1

1

10

]

]

d

0.5

a

m

r

0

[

N

0

[

q

0

u -10

-1

-0.5

-20

-2

-1

-30

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

t [s]

t [s]

t [s]

q

q

q

q

q

e [m]

e [m]

e [rad]

u

u

u

u

u

1

2

3

4

5

1

2

3

1

2

3

4

5

Nguyen Q. Hoang, Nguyen V. Khang

51

Nguyen Q. Hoang, Nguyen V. Khang

52

13

9/16/2015

Application of pseudo-inverse:

Application of pseudo-inverse:Hệ thống chân vịt đẩy Numerical experiments / with post-adjusting

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

q

J x

q

J ( x

Ke)

đó 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)] z

q

J ( x

Ke) [ E

J ( q) J( q)] z

được phân tích. Phần này cũng đề xuất một phương án phân bổ tối ưu 0

0

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.

Nguyen Q. Hoang, Nguyen V. Khang

53

Application of pseudo-inverse: Hệ thống chân vịt đẩy Application of pseudo-inverse: Hệ thống chân vịt đẩy Quan hệ lực/mômen điều khiển và lực

u4

Một số trường hợp bố trí hệ thống chân vịt

u

đẩy chân vịt

1

u3

u2

y

z

b

u

6

T

u

u

[ , u ,..., u ]

y

b

1

O

1

2

n

1

u

4



Bu

p

u



4

5

4

a

u

c

5

y

8



a

[ F , F , F , m , m , m ] T



x

b

x

y

z

x

y

z

a

x

u

a

x

a

y

6

u

3

u

u

2

1

5

3



2

Để điều khiển được tất cả 6 dof củ

2

a ROV

u

3

z

x

n

dim( u)

dim( )

6 &

(

rank

)

B

6

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

p

Yêu cầu bố trí các chân vịt sao cho ma trận B có hạng đẩy đủ (full rank) y

y

y

5

8

1

4

1

4

n

n

tàu lặn thiếu hụt dẫn động

1

4

5

5

6





p

o

x

x





x

n

n

2

3

p

o

tàu lặn có dẫn động đầy đủ, và

2

3

2

3

6

7

n

n

p

o

tàu lặn có hệ dẫn động dư (có tính chất dự phòng) 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

9/16/2015

Hệ thống chân vịt đẩy

Application of pseudo-inverse: Hệ thống chân vịt đẩy Quan hệ lực điều khiển – lực đẩy

Ma trận trạng thái động cơ

z

z

τ  Bu

Ψ  diag{ a , a , ..., a }

1

2

n

a 

0

0

0

0

c

c

c

c

1

Tiêu chuẩn tối ưu

i

- động cơ i không hỏng

b

c

0

0

0

0

s

s

s

s

1 T

1

a  0

2

i

- động cơ i hỏng (lỗi)

J 

u u 

|| u ||

min

1

1

1

1

0

0

0

0

2

B

2

2

Phân bổ lại lực điều khiển

y

x

c

c

c

c

0

0

0

0

Nghiệm tối ưu

b b

b

b

0

0

0

0

τ  BΨ u  B u



T

T 

0

0

0

0

as

as

as

as

u  B# τ  B [ BB

1

] τ

T

T 

u  ΨB [ BΨ B

1

] τ

223

0.250[ ],

0.473[ ],

20 o

b

c

m a

m

7

6

3

2

a

Sơ đồ phân bổ lực

#

Tàu lặn

y

Lực điều

u  

B 

khiển

Phát hiện lỗi

4

1

8

5

x

Application of pseudo-inverse: Hệ thống chân vịt đẩy Ví dụ

T

cos

cos

cos

cos



X

[ Y

, , M ]

B

sin

sin

sin

sin

a / 2

a / 2

a / 2

a / 2

u

[ , , , ] T

u u u u

1

2

3

4

Bốn động cơ-chân vịt làm việc bình thường

a sin

a cos

sin 2

X

1

a sin

a cos

sin 2

u

B

Y

a

2 sin 2

a sin

a cos

sin 2

M

a sin

a cos

sin 2

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])

Y

M

0

a cos

sin 2

2 sin

a

1

a sin

a cos

0

X

Y

B

u

B 

u

a sin 2 a sin

0

sin 2

2 cos

2 sin

X

M

0

0

0

2 cos

a

0

15

9/16/2015

Chương 3. Lập trình trong Matlab

3.1 Các kiểu dữ liệu

Lecture 3

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

Lập trình trong Matlab

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

Nguyen Q.Hoang

Vòng lặp FOR

Department of Applied Mechanics

Vòng lặp WHILE

Hanoi University of Science and Technology

Lệnh if, cấu trúc if - else - end

Cấu trúc switch-case

3.4 Các Mat-File

1

2

Các kiểu dữ liệu

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

• Kiểu ô, mảng (Cell-Array)

đế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ệ

– 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 nh help.

có cỡ khác nhau, có thể được sử dụng như các ô của một mả

–

ng.

Kiểu véctơ và ma trận

Các phần tử ô được gọi đến thông qua chỉ số của nó. Các phần tử

– Chuỗi ký tự (ký tự, xâu - string)

của một ô được đặt trong dấu ngoặc nhọn {…}. Ví dụ

– Kiểu ô, mảng (Cell-Array)

– Kiểu cấu trúc

>> A(1,1)={[1 2 3; 4 2 6; 1 7 9]};

>> A{1,1}

% o (1,1)

•

Kiểu véctơ và ma trận

>> A(1,2)={'Ma tran thu'};

ans =

Biểu diễn số các ma trận hai hay nhiều chiều cũng như các véctor trong

>> A(2,1)={3+8i};

1 2 3

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

>> A(2,2)={0:pi/20:2*pi};

4 2 6

đại lượng phức cần một ô nhớ gấp đôi, phần thực và ảo được ghi riêng rẽ.

1 7 9

•

Chuỗi ký tự (ký tự, xâu - string)

% goi den cac phan tu

>> A{1,2} % o (1,2)

Chuỗi ký tự được đặt trong cặp dấu nháy đơn trên

>> A(1,1)

ans = Ma tran thu

>> ’Day la chuoi ky tu’

ans = [3x3 double]

>> A(1,2)

Mỗi phần tử của chuỗi ký tự (Ma trận ký tự) cần ô nhớ 2 Byte.

ans = 'Ma tran thu'

3

4

1

9/16/2015

Các kiểu dữ liệu

Các kiểu dữ liệu

• Kiểu cấu trúc

• 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

>> A=[1 2 3; 4 2 6; 1 7 9];

>> my_structur.data % in ra

nhau. Tên của một cấu trúc gồm hai phần, tên cấu trúc trước ma tran A va inv(A)

>> my_structur = struct('data', A, 'dimension',[3 3]) dấu chấm ( . ) và tên trường trong cấu trúc sau dấu chấm, my_structur =

data: [3x3 double]

ans =

(struct_name.field_name). Các phần tử của cấu trúc được gọi

dimension: [3 3]

1 2 3

đến qua tên và chỉ số.

>> my_structur.dimension

4 2 6

– Cú phá

1 7 9

p % tao mot cau truc

ans = 3 3

>> my_structur.data

ans =

structur=struct(’name_1’, value1, ’name_2’, value2,.. ) ans = 1 2 3

4.0000 -0.5000 -1.0000

% tro den phan tu cua cau truc

4 2 6

5.0000 -1.0000 -1.0000

structur.name

1 7 9

-4.3333 0.8333 1.0000

hoặc

>> my_structur(2).data=inv(A) % mo rong truong data my_structur =

structur.name_1 = value_1;

1x2 struct array with fields:

structur.name_2 = value_2;

data

dimension

5

6

Soạn thảo Script file trong Matlab

Soạn thảo Script file trong Matlab

•

Matlab Script files là một chuỗi các lệnh

• Ví dụ soạn m-file vẽ đường tròn là hàm của bán kính r.

đượ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

8

2

9/16/2015

Soạn thảo Script file trong Matlab

Soạn thảo Script file trong Matlab

Ví dụ về chương trình không có biến vào

Ví dụ về chương trình không có biến vào

function V = volume(r)

function vol = volumef(r)

vol = (4/3)*pi*r^3;

%r = input('Hay cho ban kinh cua hinh cau')

vol = (4/3)*pi*r^3;

Save với tên file volumef.m, và gọi hàm trong command windows

disp('The tich hinh cau nay la:')

>> volumef(2)

disp(vol)

function V = Tong_volume(r1, r2)

Save với tên file volume.m, và gọi hàm trong command V1 = volumef(r1);

windows

V2 = volumef(r2);

>> volume

V = V1 + V2;

>> Tong_volume(2, 5)

9

10

Soạn thảo Script file trong Matlab

Soạn thảo Script file trong Matlab

Ví dụ về script file

• Các m-file dạng hàm

function [output 1, output 2] = function_name(input1,input2) r1 = input('Hay cho ban kinh cua hinh cau')

% Some comments that explain what the function does go here.

r2 = input('Hay cho ban kinh cua hinh cau')

global x y z

disp(‘Tong the tich hai hinh cau’)

persistent a b c

Tong_volume(r1, r2)

MATLAB command 1;

MATLAB command 2;

MATLAB command 3;

output 1 = ….;

Save với tên file TinhV.m, và gọi hàm trong command output 2 = ….;

windows

>> tinhV

11

12

3

9/16/2015

Soạn thảo Script file trong Matlab

Các vòng lặp và rẽ nhánh

Lệnh for

for var_= bien_chay

• Biến toàn cục ( global variable) và biến cục bộ hay biến địa phương các lệnh

( local variable). Biến toàn cục được khai báo theo cú pháp end

global var_1 var_2

Lệnh while

while dieukien,

• 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 các lệnh

[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 end

được liệt kê trong cửa sổ workspace, và do đó không thể tác động đến Lệnh

được trong môi trường làm việc.

if (nếu .. thì ..)

if dieukien

...

• Biến persistent trong một hàm có thể được thống nhất với cách khai

<elseif> ……. <else> …….

báo

end

persistent var_1 var_2

Lệnh switch (chuyển)

switch ..

• Trái với các đại lượng đã được khai báo bằng global, các đại lượng với case ..

khai báo persistent chỉ được biết đến trong hàm mà nó được khai báo.

<otherwise ..>

Do đó các hàm khác không thể tác động được đến biến này. Các biến end

persistent chỉ được xóa khi hàm chứa nó thoát ra khỏi bộ nhớ (clear Lệnh

function_name) hoặc khi hàm được thay đổi sau đó được ghi lại.

ngắt vòng lặp trong for break;

và while

13

14

Các vòng lặp và rẽ nhánh

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 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

% su dung vong lap while

% su dung vong lap for

function tongS

function tong_Sn

function tong = tongvector(v)

n=100;

% ham tinh tong cac phan tu cua mot vector

S=0; i=1;

n=100;

n=length(v); % tra lai so phan tu, chieu dai cua vector while i<=n

S=0;

s = 0; % khoi gan cho tong

S = S+1/i;

for i=1:n

for i = 1:n % hoac for i = 1 : 1 : n

i = i+1;

S=S+1/i;

s = s + v(i);

end

end

end;

tong=S;

tong=S;

tong = s;

disp('Voi n = 100, tong la :')

disp('Voi n = 100, tong la :')

% save with file name tongvector.m

disp(tong)

disp(tong)

>> x=[1:2:50];

>> >> tong_S

>> >> tong_S

Sử dụng hàm

>> s=tongvector(x)

Voi n = 100, tong la :

Voi n = 100, tong la :

s = 625

5.1874

5.1874

15

16

4

9/16/2015

Các vòng lặp và rẽ nhánh

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 Lệnh if, cấu trúc if - else - end

sau nhỏ hơn

Lệnh điều kiện nếu - thì. Cú pháp

27

1

1

1

n 1

S  1 

  ... 

   27

IF expression

2

3

n

statements 1

i 1 i

ELSEIF expression

function n = timn(S)

format long

statements 2

% tim so n lon nhat thoa man

tong=0; i=0;

ELSE

% tongS(n) = 1+1/2+1/3+...+1/n

while(tong<S)

statements 3

< S

i = i+1;

END

% gia tri S > 1 vao tu ban phim

tong = tong+1/i;

T=1;

end

while T

n=i-1;

disp('Vao so duong S')

tong=tong-1/i;

S = input('S = ')

disp('So n tim duoc la: ')

if (S>1 & S<10) T = 0; end

disp(n)

end

disp('Gia tri cua tong S(n) la: ')

format long

disp(tong)

17

18

Các vòng lặp và rẽ nhánh

Các vòng lặp và rẽ nhánh

k=5;

% ket qua chuong trinh

Cấu trúc switch-case

for i=1:k

>> A

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 for j=1:k

A =

như sau

if i==j

4 0 1 0 0

SWITCH switch_expr

A(i,j)=4;

0 4 0 1 0

CASE case_expr,

elseif(abs(j-i)==2)

1 0 4 0 1

statement, ..., statement

A(i,j)=1;

0 1 0 4 0

CASE {case_expr1, case_expr2, case_expr3,...}

else

0 0 1 0 4

statement, ..., statement

A(i,j)=0;

...

end

OTHERWISE,

end

statement, ..., statement

end

END

19

20

5

9/16/2015

Các vòng lặp và rẽ nhánh

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 diem = 1; % 2, 3, 4, 5, 6

chấm mat (filename.mat). Mat-file là những tệp nhị phân được nén switch diem

để lưu trữ kết quả số. Các file này được sử dụng để ghi các kết quả

case 1

mà đã được tạo ra bởi chuỗi các lệnh, chỉ thị của Matlab. Chẳng disp('Diem gioi')

hạn, để ghi giá trị của hai biến x1 và x2 trong file có tên dulieu.mat case 2

ta đánh dòng lệnh

disp('Diem tot')

>> x1=5;

case 3

>> x2=10;

disp('Diem kha')

>> save dulieu.mat x1 x2

case 4

disp('Diem trung binh')

Muốn tải (đọc) một mat-file vào MATLAB ta gõ dòng lệnh otherwise

disp('Diem kem.')

>> load filename (or load filename.mat)

end

21

22

Bài tập

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ụ, Viết thuật toán (thuật giải)

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 2

của hình trụ

f ( x)  x  a  0, a  0, x  a

.

2. Viết một chương trình con có sử dụng lệnh while để tính tổng S(x,n) k  1;, x( k)  1;

2

( , )  1  

 ...

n

S x n

x x

 x

f ( 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

x

 x 

pương pháp lặp Newton-Raphson, theo công thức [heron]?

k 1

k

f (

 x )

k

2

x  a

x

Với x(1) = 1, độ chính xác yêu cầu

a

2

2





n

n

x







là 0.0000001

x

a

x

a

x

a

k

k

k

n 1











2 x

2

2 x

x

7



k

n

n

| x

 x | 10

2 x

2 x

2

2 x

n 1

n

k

k

k





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 1

a

  

cho b, a > b.

x



2

k



x 



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

Chương 4. Đồ họa trong Matlab

Lecture 4

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

Đồ họa trong Matlab

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ị

Nguyen Q.Hoang

Vẽ đồ thị các hàm có điểm không xác định

Department of Applied Mechanics

Hanoi University of Science and Technology

4.2 Các lệnh vẽ trong không gian – 3D

Nguyen Quang Hoang

1

Nguyen Quang Hoang

2

Department of Applied Mechanics

Department of Applied Mechanics

Đồ thị trong mặt phẳng – 2D

Đồ thị trong mặt phẳng – 2D

• Trướ

• Đặt màu và kiểu đường cho đồ thị

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:

– Khi muốn thể hiện nhiều đồ thị trong một hình người ta thường

– Định nghĩa hàm cần vẽ, y = f(x)

phân biệt các đồ thị bằng các màu, và kiểu đường khác nhau (nét

– Xá

liền, nét đứt)

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.

– 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.

– Ví dụ:

1

– Chia khoảng cần vẽ [0,10] bằng các điểm chia cách đều, sử dụng toán tử

• >> x=[0:0.01:10]; y1=sin(x);

(:) , chẳng hạn x = [a:h:b], với h là bước chia.

• >> y2=cos(x); plot(x,y1, x,y2,'--') 0.5

– 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ử.

0

– Gọ

1

i hàm plot(x,y).

-0.5

Cụ thể:

0.5

>> x = [0:0.1:10];

0

-10

2

4

6

8

10

>> y = sin(x);

- 0.5

>> plot(x,y)

- 10

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

Nguyen Quang Hoang

4

Department of Applied Mechanics

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

Ký

Màu

• lệnh legend cho phép ta điền các chú giải lên hình vẽ

hiệu

- solid line

nét liền

k black – đen

x [m]

-- dashed line

nét đứt

•

3

b blue – xanh da trời

Ví dụ:

v [m/s]

2

-. dashdot line nét chấm gạch

c cyan – lục lam

>> t = [0:0.1:10];

1

: dotted line

nét chấm

g green – xanh lá cây

>> x = sin(2*t); v = 2*cos(2*t);

0

(none) no line

(không hiện)

>> plot(t,x,'k-', t,v,'k--'),

r red – đỏ

-1

>> grid on, xlabel('t [s]')

-2

m magenta– đỏ tươi

Kiểu đánh dấu (marker)

0

1

2

3

4

5

6

7

8

9

10

>> legend('x [m]','v [m/s]')

y yelow – vàng

t [s]

. point

điểm

>> axis([0 10 -2.6 3.5])

o circle

vòng tròn

w white – trắng

x x-mark

chữ x

+ plus

dấu cộng

* star

dấu sao

s square

dấu vuông

Sử dụng lệnh trợ giúp:

…..

……

>> help plot

Nguyen Quang Hoang

5

Nguyen Quang Hoang

6

Department of Applied Mechanics

Department of Applied Mechanics

Đồ thị trong mặt phẳng – 2D

Đồ 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í dụ:

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), t = [0:0.01:8]; x = sin(3*t); v=3*cos(3*t); a=-9*sin(3*t); trong đó m và n cho biết số hàng và số cột trong hình vẽ như là một subplot(3,1,1), plot(t,x,'k-','Linewidth',1), grid on, ylabel('x [m]'), 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.

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

xlabel('t [s]‘)

1

]m[ 0

plot 1 x

-1

0

1

2

3

4

5

6

7

8

>> for i = 1:6 subplot(2,3,i), end

2

]s/

plot 2 m 0

Fig. 1

Fig. 2

Fig. 3

[ v -2

0

1

2

3

4

5

6

7

8

10

]2s

Fig. 4

Fig. 5

Fig. 6

/

0

m

plot 3 [ a

-100

1

2

3

4

5

6

7

8

t [s]

Nguyen Quang Hoang

7

Nguyen Quang Hoang

8

Department of Applied Mechanics

Department of Applied Mechanics

2

9/16/2015

Đồ thị trong mặt phẳng – 2D

Vẽ hàm y=y(x), biết x = x(t), y=y(t)

• Đồ thị theo tọa độ cực và Đồ thị với thang chia Logarith 0.25

x

a sin( t)

90 8

0.2

120

60

o 0.15

d

6

nei 0.1

B

4

150

30

0.05

0

y

a cos( t)

y

y( x)

2

-1

0

1

2

10

10

10

10

180

0

0

-1

ahp -2

210

330

co

G

-3

240

300

-4 -1

0

1

2

270

10

10

10

10

log()

Lệnh polar(phi,r)

Ngoài ra ta có thể vẽ các đồ thị với

a = 1;

lệnh semilogx hoặc semilogy, với phi=[0:pi/90:2*pi];

các lệnh này chỉ có một trục x hoặc y

r = a*phi;

tương ứng sử dụng thang logarith còn

polar(phi, r, '-k');

trục kia sử dụng thang thập phân

Nguyen Quang Hoang

9

Nguyen Quang Hoang

10

Department of Applied Mechanics

Department of Applied Mechanics

Đồ thị trong mặt phẳng – 2D

Đồ thị trong mặt phẳng – 2D

Các lệnh plot 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 plot(x,y)

Vẽ đồ thị bằng cách nối thẳng các điểm [x(i),y(i)]

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 comet(x,y)

Vẽ hoạt hình một quĩ đạo

đường, màu sắc, đánh dấu (line style, colour, marker) được làm sáng tỏ

stair(x),

trong thí dụ sau.

Vẽ đồ thị kiểu bậc thang

stairs(x,y,..)

figure(1)

stem(x,y)

clf % xoa tat ca cac do thi hien co

spy(matrix)

hiển thị một ma trận thưa

t = [0:pi/40:2*pi];

semilogx(x,y)

Hiển thị y(x), sử dụng thang logarith cho trục x plot(t,sin(t),'-.r*'), % ve do thi 1

hold on

semilogy(x,y)

Hiển thị y(x), sử dụng thang logarith cho trục y plot(t,sin(t-pi/2),'mo')

% ve do thi 2

loglog(x,y)

Hiển thị y(x), sử dụng thang logarith cho trục x và y plot(t,sin(t-pi),':bs')

% ve do thi 3

fplot(ham,mien) Vẽ hàm số trong miền

axis([0, 2*pi, -1.2, 1.2])

hold <on| off> giữ đối tượng sẵn có, đóng lại hold off

grid on

Nguyen Quang Hoang

11

Nguyen Quang Hoang

12

Department of Applied Mechanics

Department of Applied Mechanics

3

9/16/2015

Đồ thị trong mặt phẳng – 2D

Đồ 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ử: Các trục tọa độ và các chú thích trên đồ thị

.*, ./ , .

Điều chỉnh phạm vi các trục tọa

axis([xmin, xmax, ymin, ymax])

Ví dụ:

độ

Tự động chọn xmin

axis([-inf, xmax, ymin, ymax])

0.40 x 2

y

e

x

sin2 x

Sử dụng chia lưới hay không

grid <on | off >

chia lưới

5

clc,

clear all

axis <off | on | auto| equal | square>

Năm lựa chọn đối với các trục

4

x=[0:0.05:12];

3

xlabel(String)

Chú thích cho trục x

y=exp(-0.4*x).*x.^2+sin(2*x); y

ylabel(String)

Chú thích cho trục y

plot(x,y,'k-','linewidth',2), 2

Hiển

grid on

title(String)

thị tên đồ thị

1

xlabel('x'), ylabel('y')

text(String)

Chèn chữ vào đồ thị

legend(string_1, string_2,… <,position>)

0

Hiển thỉ các chú giải

0

2

4

6

8

10

12

x

Nguyen Quang Hoang

13

Nguyen Quang Hoang

14

Department of Applied Mechanics

Department of Applied Mechanics

Đồ thị trong mặt phẳng – 2D

Đồ thị trong không gian – 3D

Lưu ý khi vẽ đồ thị của các hàm f(x) có điểm không xác định. Ví dụ

Các lệnh vẽ trong không gian – 3D

2

4(1  x )

f ( x) 

x

0.7990,

x

1.3380

2

2

(10  7 x )(1  x )  2

1

2

plot3(x,y,z <,plotstil>)

surfl(x,y,z <,color>)

comet3(x,y,z <,comet length>) patch(x,y <,z> ,color) mesh(x,y,z <,color>)

waterfall(x,y,z..<,..>,..)

surf(x,y,z <,color>)

contour3(x,y,z.. <,..>,..)

surfc(x,y,z <,color>)

contour(x,y,z.. <,..>,..)

Ví dụ

2

2

x

y

Cần thể hiện đồ thị hàm

z

ye 



trong miền

2

  x, y  2

Nguyen Quang Hoang

15

Nguyen Quang Hoang

16

Department of Applied Mechanics

Department of Applied Mechanics

4

9/16/2015

Đồ thị trong không gian – 3D

Đồ thị trong không gian – 3D

• [x,y] = meshgrid(-2:0.1:2);

• Ví dụ:

2

2

0

 .2 ( x  y )

2

2

z  f ( x, y)  e

 sin( x  y ),

3

  x, y  3

• z = y.*exp(-x.^2-y.^2);

• mesh(x,y,z), xlabel('x'), ylabel('y'), zlabel('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,z) surf(x,y,z) surfc(x,y,z) mesh(X,Y,f); surf(X,Y,f) contour3(X,Y,f,30); Nguyen Quang Hoang

17

Nguyen Quang Hoang

18

Department of Applied Mechanics

Department of Applied Mechanics

5

16/09/2015

Chương 5. Lập trình trong Matlab

5.1 Mịn hóa đường cong bằng đa thức

Lecture 5

5.2 Nội suy đa thức

Nội suy và làm mịn đường cong

5.3 Bài tập thực hành

Nguyen Q.Hoang

Department of Applied Mechanics

Hanoi University of Science and Technology

Nguyen Quang Hoang

1

Nguyen Quang Hoang

2

Department of Applied Mechanics

Department of Applied Mechanics

Nêu bài toán

Phương pháp sai số bình phương bé nhất

• Trong kỹ thuật khi đo đạc ta nhận được các số liệu rời rạc, từ các số

Trường hợp đơn giản nhất, n = 1, đa thức có dạng một phương trình đường liệu đo này ta cần biểu diễn các đường để biết được đặc tính của hệ

thẳng

khảo sát và có thể suy ra được những giá trị lận cận các giá trị đã đo.

y

ax

b

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ẽ

Theo phương pháp sai số bình phương bé nhất hai hệ số a và b được tìm như

trình bày các kỹ thuật đơn giản thực hiện công việc đó.

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

• 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 e  ( ax  b)  y , i  1,2,..., m

i

i

i

n

n 1



2

1

p ( x)  p x  p x

 ...  p x  p x  p

n

1

2

n 1



n

n 1



Như thế tổng bình phương các sai số là

các hệ số của đa thức này được suy ra hai véctơ số liệu x và y, nhờ

m

m

2

2

S( ,

a b)   e  ( ax  b  y ) phương pháp sai số bình phương bé nhất (least mean square).

i

i

i

i 1



i 1



x

[ x x x ... x ], y

[ y y y ... y ] m  n Bài toán đặt ra ở đây là tìm a, b sao cho hàm S đạt cực tiểu.

1

2

3

m

1

2

3

m

Nguyen Quang Hoang

3

Nguyen Quang Hoang

4

Department of Applied Mechanics

Department of Applied Mechanics

1

16/09/2015

Phương pháp sai số bình phương bé nhất

Phương pháp sai số bình phương bé nhất

y

e  y  P ( x )

Đạo hàm S theo các biến a và b ta nhận được

i

i

n

i

y  P ( x)

m

m

S



S



n

 2( ax  b  y ) x  0,

 2( ax  b  y )  0

e

i

i

i

i

i

a





b



i

i 1

i 1



m

m

m

m

m

( x , y )

a x x  b x   x y  0, a x  mb   y  0

i i

i

i i

i

i

i

i

i 1



i 1



i 1



i 1



i 1



Từ đây giải đượ

x

c

n  1, P ( x)  ax  b 1

m

m

m

y

y  x x  x  x y

 x y  mx y

i i

i i

i i

i 1



i 1



i 1

b







m

,

a

m

y  ax  b

2

2

.

 x x  mx

 x x  mx

i i

i i

e

i 1



i 1



i

m

với các giá trị trung bình

 y ( x  x)

i

i

( x , y )

e  y  ax  b

1 m

1 m

i 

i

i

i

i

i



 

x 

 x ,

y 



1

y

a

m

,

b

y

xa

i

i

m

 x ( x  x)

i 1



m i 1



i

i

x

i 1



Nguyen Quang Hoang

5

Nguyen Quang Hoang

6

Department of Applied Mechanics

Department of Applied Mechanics

Phương pháp sai số bình phương bé nhất

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) Sử dụng ma trận tựa nghịch đảo (nghịch đảo trái) Quan hệ tuyến tính

Quan hệ bậc 2

x a

b

y ,

i

1,2,..., m

2

i

i

x p

x p

p

y ,

i

1,2,..., m

i

1

i 2

3

i

x

1

y

1

1

2

x

x

1

y

x

1 a

y

1

1

1

p

2

2

Ap

y

2

1

x

x

1

y

... ... b

...

2

2

2

p

Ap

y

2

...

... ...

...

x

1

y

p

m

m

T

1

[

]

T

p

A A A y

2

3

x

x

1

y

T

1

T

m

m

m

p

[A A] A y

Nguyen Quang Hoang

7

Nguyen Quang Hoang

8

Department of Applied Mechanics

Department of Applied Mechanics

2

16/09/2015

Phương pháp sai số bình phương bé nhất

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

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 x

6

8

10

12

14

16

18

20

22 24

có bậc lũy thừa lớn nhất, ta có

>> a=p(1)

y 3.94 3.8 4.1 3.87 4.45 4.33 4.12 4.43 4.6 4.5

a = 0.0392

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

>> b=p(2)

5

tìm các giá trị của a, b và tìm giá trị của y khi x = 11.

b = 3.6267

Trước hết ta cần nhập hai véctơ số liệu x và y vào Matlab.

4

Giá trị cần tìn tại x =

>> x = [6: 2 : 24]

11 được tính bằng

x =

lệnh polyval(p,x)

3

6 8 10 12 14 16 18 20 22 24

>> polyval(p,11)

y

>> y = [3.94 3.8 4.1 3.87 4.45 4.33 4.12 4.43 4.6 4.5]

ans = 4.0574

2

y =

3.94 3.80 4.10 3.87 4.45 4.33 4.12 4.43 4.60 4.50

1

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)

0

p = 0.0392 3.6267

6

8

10

12

14

16

18

20

22

24

x

Nguyen Quang Hoang

9

Nguyen Quang Hoang

10

Department of Applied Mechanics

Department of Applied Mechanics

Phương pháp sai số bình phương bé nhất

Mịn hóa bằng đa thức

Xét tổng bình phương của sai số là

Ví dụ 2. Cho bộ số liệu sau

>> w = a*x+b;

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

>> e = w-y;

y

300 281 261 244 228 214 202 191 181 164 151 149 141

>> S1 = sum(e.*e)

S1 = 0.2274

Nếu các điểm trên được làm mịn bằng đường bậc 2: Hãy tìm một đa thức bậc 3 phù hợp với số liệu cho. So sánh sai số so với

>> n=2;

trường hợp nếu sử dụng đa thực bậc hai.

>> p = polyfit(x,y,n);

Trước hết nhập số liệu vào hai véctơ và nối các điểm ta được p =

>> x = [0, 0.5, 1.0,1.5, 2, 2.5,3, 3.5,4, 5, 6,

-0.0002 0.0443 3.5940

6.1, 7];

>> a=p(1); b=p(2); c = p(3);

>> y = [300,281,261,244, 228,214,202,191,181,164,151,

>> w = a*x.^2+b*x+c;

149,141];

>> e = w-y;

>> plot(x,y,'ko', x,y, 'k-'), grid on, xlabel('x'),

>> S2= sum(e.*e)

ylabel('y')

S2 = 0.2272

Ta thấy sai số có nhỏ hơn (S2<S1) nhưng không đáng kể.

Nguyen Quang Hoang

11

Nguyen Quang Hoang

12

Department of Applied Mechanics

Department of Applied Mechanics

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

>> n=3; % Xap xi duong bac 3

Trong nhiều trường hợp khi mịn hóa bằng đa thức không đáp ứng được

>> p = polyfit(x,y,n)

>> p = polyfit(x,y,n)

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

>> a=p(1); b=p(2);

>> a=p(1); b=p(2);

thông qua hàm e mũ là một phương án nên được xem xét. Giả sử bộ số

>> c = p(3);

c = p(3); d = p(4);

liệu (x,y) có thể lấy xấp xỉ bằng hàm

>> w = a*x.^2+b*x+c;

>> w = a*x.^3+b*x.^2+c*x+d;

ax



>> e = w-y;

y

be

>> e = w-y;

>> S2= sum(e.*e)

>> S3= sum(e.*e)

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 S2 = 15.4079

S3 = 2.7933

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. Để

350

300

tránh điều đó, ta lấy lôgarit cơ số tự nhiên hai vế được 300

250

ln y  ax  ln b

250

y

y

200

200

Đặ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ố

150

p1, p2 để xấp xỉ các số liệu (z, w) = (x, ln(y)) bằng đường thẳng 150

100





0

1

2

3

4

5

6

7

w

p z

p

1000

1

2

3

4

5

6

7

x

1

2

x

Nguyen Quang Hoang

13

Nguyen Quang Hoang

14

Department of Applied Mechanics

Department of Applied Mechanics

Mịn hóa bằng hàm e mũ

Mịn hóa bằng hàm e mũ

Ví dụ xét bộ số liệu sau

x

1.2

2.8

4.3

5.4

6.8

7.9

Xấp xỉ bằng hàm e mũ

300

y

7.5

16.1

38.9

67.0

146.6

266.2

>> p = polyfit(x,log(y),1)

p =

0.5366 1.3321

250

Nhập số liệu vào Matlab

>> x = [1.2 2.8 4.3 5.4 6.8 7.9]; 200

>> a = p(1); b=exp(p(2));

>> y = [7.5 16.1 38.9 67.0 146.6 266.2];

>> w = b*exp(a*x);

y

150

300

>> plot(x,y,'ko', x,y,

>> err = w-y;

'k-'),

250

>> S3= sum(err.*err)

100

grid on, xlabel('x'),

Se = 17.6259

200

ylabel('y‘)

>> x3 = [0:0.1:8];

50

y

150

>> y3 = b*exp(a*x3);

0

100

>> plot(x,y,'ko', x3,y3, 'k-'),

0

1

2

3

4

5

6

7

8

grid on, xlabel('x'),ylabel('y')

x

50

01

2

3

4

5

6

7

8

x

Nguyen Quang Hoang

15

Nguyen Quang Hoang

16

Department of Applied Mechanics

Department of Applied Mechanics

4

16/09/2015

Nội suy bằng đa thức

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 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*

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

<= xn. Bài toán này trong Matlab được giải quyết dễ dàng với lệnh thức cũng cho ta một đường cong nhưng đi qua các điểm số liệu đó. Có interp. Cách gọi lệnh này như sau:

nhiều phương pháp để tạo ra đa thức đi qua các điểm cho như phương y_star = interp1(x,y,x_star)

pháp nội suy theo công thức Lagrange, phương pháp Newton, phương 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 =

pháp Neville, phương pháp nội suy với đường spline bậc 3, .... Chi tiết length(y), thì ta có thể gọi lệnh

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à y_star = interp1(y,x_star)

phương pháp số.

Để nói rõ hơn phương pháp nào được sử dụng trong việc nội suy này ta Trong phần này chỉ trình bày việc sử dụng các lệnh interp của Matlab sử dụng cú pháp

để 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 y_star = interp1(x, y, x_star, method)

là cho biết các giá trị tương ứng của hai đại lượng x và y như trong bảng với method được chọn là một trong các phương pháp sau số liệu

'nearest' - nội suy theo lân cận gần nhất

'linear' - nội suy tuyến tính

x

x1

x2

x3

…

xn

'spline' - nội suy sử dụng đường spline bậc 3

'pchip' - nội suy bậc 3 từng đoạn

y

y1

y2

y3

…

yn

'cubic' - tương tự như 'pchip'

Nguyen Quang Hoang

17

Nguyen Quang Hoang

18

Department of Applied Mechanics

Department of Applied Mechanics

Nội suy bằng đa thức

Nội suy bằng đa thức

Nếu giá trị x* nằm ngoài khoảng [x1, xn], khi đó ta gọi phép tính giá trị y*

Như thế các phương pháp nội suy khác nhau sẽ cho ta các kết quả khác

= y(x*) là phép tính ngoại suy và sử dụng lệnh với cú pháp: 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ụ

y_star = interp1(x, y, x_star, method, 'extrap') sau

x = 0:8; y = sin(x); plot(x,y,'-ko'), hold on Sau đây xét ví dụ:

xi = 0:0.1:8;

>> x = 1:20; y = sin(x);

yi1 = interp1(x,y, xi, 'linear');

>> x_star = 1.5;

yi2 = interp1(x,y, xi, 'cubic');

>> y_star = interp1(x,y, x_star)

yi3 = interp1(x,y, xi, 'spline');

y_star = 0.8754

>> y_star = interp1(x,y, x_star, 'linear') plot(xi,yi1,'-k','linewidth',1.5)

y_star = 0.8754

plot(xi,yi2,'-r','linewidth',1.5)

>> y_star = interp1(x,y, x_star, 'cubic')

plot(xi,yi3,'-b','linewidth',1.5), grid on

y_star = 0.9008

>> y_star = interp1(x,y, x_star, 'nearest') xlabel('x'), ylabel('y'), axis([0 8, -1.1 1.1])

y_star = 0.9093

legend('points','linear','cubic','spline')

>> y_star = interp1(x,y, x_star, 'spline') y_star = 1.0200

Nguyen Quang Hoang

19

Nguyen Quang Hoang

20

Department of Applied Mechanics

Department of Applied Mechanics

5

16/09/2015

Nội suy bằng đa thức

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

1

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

0.5

y

0

4. Tỷ trọng tương đối của không khí thay đổi theo chiều cao đo được như

trong bảng.

- 0.5

points

linear

h (km)

0

1 . 525

3 . 050

4 . 575

6 . 10

7 . 625

9 . 150

cubic

ρ

1

0 . 8617 0 . 7385 0 . 6292 0 . 5328

0 . 4481 0 . 3741

- 1

spline

0

1

2

3

4

5

6

7

8

Tìm phương trình bậc hai thể hiện quan hệ giữa độ cao và tỷ trọng không x

khí, từ đó tính tỷ trọng của không khí ở độ cao h = 10.5 km.

2

  a  bh  ch

Nguyen Quang Hoang

21

Nguyen Quang Hoang

22

Department of Applied Mechanics

Department of Applied Mechanics

Bài tập

Bài tập

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 5. Độ nhớt động học μ của nước thay đổi theo nhiệt độ T được cho như trong bảng k

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 dưới đây

tiêu hao nhiên liệu.

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 Model

m (kg)

y (km/lít)

nhiệt độ T = 10◦, 30◦, 60◦, và 90 ◦C

Contour

1310

10.2

Crown

1810

8.1

T (◦C)

0

21 . 1

37 . 8

54 . 4

71 . 1

87 . 8

100

Escort

1175

11.9

μ

1 . 79

1 . 13

0 . 696

0 . 519

0 . 338

0 . 321

0 . 296

Expedition

2360

5.5

k

Explorer

1960

6.8

(10−3 m2 / s)

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

Nguyen Quang Hoang

23

Nguyen Quang Hoang

24

Department of Applied Mechanics

Department of Applied Mechanics

6

16/09/2015

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

Lecture 6

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

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

1

Nguyen Quang Hoang

2

Department of Applied Mechanics

Department of Applied Mechanics

Ví dụ giải phương trình vi phân thường

Giải số bằng lệnh ode23

Tìm hàm y(t) thỏa mãn phương trình vi phân

Cú pháp lệnh ode23

dy  cos( t)

với điều kiện đầu

(

y 0)  2

dt

[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) ,

Tên m-file chứa phương trình vi

(

y t)  sin( t)  C

phân cần giải

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 Khoảng thời gian cần giải

đầu nhận được hằng số C:

t từ t0 đến t_end

(

y 0)  sin(0)  C  2



C  2

Kết quả tính được ghi vào hai véctơ

Giá trị điều kiện đầu

(

y t)  sin( t)  2

t và y chứa các cặp [t(i), y(ti)]

Nguyen Quang Hoang

3

Nguyen Quang Hoang

4

Department of Applied Mechanics

Department of Applied Mechanics

1

16/09/2015

Giải số bằng lệnh ode23

Ví dụ giải phương trình vi phân thường

Cú pháp lệnh ode23

Ví dụ: sử dụng ode23 giải phương trình vi phân dy  cos( t)

với điều kiện đầu

(

y 0)  2

[t, y] = ode23(‘funct_name’, [t0:h:t_end], y0) dt

Trong khoảng thời gian t = [0 2*pi]

Tên m-file chứa phương trình vi

phân cần giải

Trước hết tạo m-file chứa phương trình vi phân cần giải: function ydot = eq1(t,y)

Khoảng thời gian cần giải

ydot = cos(t);

Ghi lại với tên file là tên

t từ t0 đến t_end, với h là

hàm eq1.m

bước thời gian

end

Kết quả tính được ghi vào hai véctơ

Sử dụng lệnh ode23:

Giá trị điều kiện đầu

t và y chứa các cặp [t(i), y(ti)]

>> [t, y] = ode23(‘eq1’, [0, 2*pi], 2)

Nguyen Quang Hoang

5

Nguyen Quang Hoang

6

Department of Applied Mechanics

Department of Applied Mechanics

Ví dụ giải phương trình vi phân thường

Giải số bằng lệnh ode45

• Vẽ kết quả và so sánh với nghiệm chính xác Cú pháp lệnh ode45

>> f = sin(t) + 2;

>> plot(t, f,'k-'), hold on

[t, y] = ode45(‘funct_name’, [t0, t_end], y0)

>> axis([0 2*pi 0 4])

4

>> plot(t,y,'ko')

Tên m-file chứa phương trình vi

3

phân cần giải

>> xlabel('t'),

>> ylabel('y(t)') )t 2

y(

Khoảng thời gian cần giải

t từ t0 đến t_end

1

Kết quả tính được ghi vào hai véctơ

Giá trị điều kiện đầu

0

t và y chứa các cặp [t(i), y(ti)]

0

1

2

3

4

5

6

t

Nguyen Quang Hoang

7

Nguyen Quang Hoang

8

Department of Applied Mechanics

Department of Applied Mechanics

2

16/09/2015

Giải số bằng lệnh ode45

Giải số bằng lệnh ode23

Ví dụ: sử dụng ode45 giải phương trình vi phân

>> [t, y2] = ode45(‘eq1’, [0:0.01:30], 0); dy

với điều kiện đầu

>> plot(t,y2)



(

y 0)  0

5

 y  5 F( t)

dt

t

 /

( )

tc

F t  te

cos( t

 )

tc  0.1,

  10

Trong khoảng thời gian t = [0 30]

function ydot = eq2(t,y)

tc = 0.1; w = 10;

Ghi lại với tên

Ft = t*exp(-t/tc)*cos(w*t);

file là tên

ydot = -5*y + 5*Ft;

hàm eq2.m

end

>> [t, y] = ode45(‘eq2’, [0, 30], 0);

>> % [t, y] = ode45(‘eq2’, [0:0.01:30], 0); y2(t)

y(t)

>> plot(t,y)

Nguyen Quang Hoang

9

Nguyen Quang Hoang

10

Department of Applied Mechanics

Department of Applied Mechanics

Giải số bằng lệnh ode45

Hệ phương trình vi phân cấp 1

Ví dụ: sử dụng ode45 giải phương trình vi phân Hệ n phương trình vi phân cấp 1.

dy

với điều kiện đầu



(

y 0)  0

5

 y  5 F( t)

dt

dy 1

y 

 f ( t, y ,..., y ) y (0)  y

t

 /

( )

tc

F t  te

cos( t

 )

tc  0.1,

  10

1

1

1

n

dt

1

10

dy

y (0)  y

2

2

20

y 

 f ( t, y ,..., y ) Các điều kiện đầu

...

Trong khoảng thời gian t = [

2

2

1

n

0 30]

dt

...

y (0)  y

n

n 0

function ydot = odevidu3(t,y, tc, w)

dyn

y 

 f ( t, y ,..., y ) n

n

1

n

dt

Ghi lại với tên

Ft = t*exp(-t/tc)*cos(w*t);

file là tên hàm

Hay viết ở dạng véc tơ

ydot = -5*y + 5*Ft;

odevidu3.m

end

d

 y

y

 f( t, )

y ,

(

y 0)  y0

dt

y  [ y , y ,..., y ] T , f  [ f , f ,..., f ] T

>> [t,y] = ode45(@(t,y)odevidu3(t,y,0.1, 10), [0: 0.01:30], 0); 1

2

n

1

2

n

>> plot(t,y)

Nguyen Quang Hoang

11

Nguyen Quang Hoang

12

Department of Applied Mechanics

Department of Applied Mechanics

3

16/09/2015

Phương trình vi phân cấp cao

Phương trình vi phân cấp cao

Xét 1 phương trình vi phân cấp n.

Ví dụ hạ bậc phương trình

x(0)  x

vi phân cấp 2 sau

n

d x

0

( n)

( n 1

 )

x



 f ( t, x, x, x,..., x

)

x(0)  x

Các điều kiện đầu

0

x  f ( t, x, x), f ( t, x, x)  b

 x  cx  F sin t



n

dt

...

0

( n 1

 )

( n 1

 )

x

(0)  x 0

Đặt các biến mới và hạ bậc

Đặt các biến mới và hạ bậc

y  x,

y  x  y

1

1

2

y  x,

y  x  y

y  x,

y  x  y

1

1

2

2

2

3

y  x,

y  x  f ( t, y , y )

...

...

2

2

1

2

( n 2)

y

 x

( n 1

 )

y

 x

 y

Cụ thể trong trường hợp trên

n 1



n 1



n

( n 1

 )

y  x

( n )

y  x

 f ( t, y ,..., y ) n

n

1

n

y  y

1

2

y  f( t, )

y

Hay viết gọn lại dạng véc tơ

y  b

 y  cy  F sin t



2

2

1

0

d

 y

y

 f( t, )

y ,

(

y 0)  y

y

 



y



0

dt

1

2

y   ,

f( t, )

y  



y  [ y , y ,..., y ] T , f  [ y , y ,..., y , f ( t, ) y ] T

y

b

 y  cy  F sin t



1

2

n

2

3

n

 2

 



2

1

0





Nguyen Quang Hoang

13

Nguyen Quang Hoang

14

Department of Applied Mechanics

Department of Applied Mechanics

y

x

y

function ydot = daodong1dof(t,y)

1

2

Ví dụ hệ dao động 1dof

mx

cx

kx

F sin( t)

1

% t =1; y = [1; 2]

0

y

x

F sin( t) cy

ky

2

0

2

1

m = 1; b = 0.0; c = 16;

m

>> [t, y] = ode45('daodong1dof',

F0 = 10; Ome = 4.1;

[0:0.01:100], [1, 2])

function ydot = vibr1dof(t,y, m,c,k,F0, Ome)

>> plot(t,y(:,1))

ydot = zeros(2,1);

ydot = zeros(2,1);

Ft = F0*sin(Ome*t);

Ft = F0*sin(Ome*t);

ydot(1) = y(2);

ydot(1) = y(2);

ydot(2) = 1/m*(Ft - b*y(2) - c*y(1));

ydot(2) = 1/m*(Ft - c*y(2) - k*y(1));

10

end

end

5

0

>> [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])

-5

>> plot(t,y(:,1))

-100

10

20

30

40

50

60

70

80

90

100

Nguyen Quang Hoang

15

Nguyen Quang Hoang

16

Department of Applied Mechanics

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

4

]

1

s/m 0

0

m[ 2

2

y

-1

4

]s

-2

/

-2

m

3

0

m[

-3

2y

2

-4

-4

0

10

20

30

40

50

60

70

-1

-0.5

0

0.5

1

70

-2

1

t[s]

0

60

-4

-1

0

10

20

30

40

50

60

70

50

t[s]

1

-2

40

-3

]s[t

0.5

-4

30

-1

-0.5

0

0.5

1

]mm 0

20

[ 1y

10

-0.5

0

1

5.

0

5.

-1

0

0-

1

]

m

m

[

y

-10

10

20

30

40

50

60

70

t[s]

Nguyen Quang Hoang

17

Nguyen Quang Hoang

18

Department of Applied Mechanics

Department of Applied Mechanics

Giải hệ phương trình vi phân cấp 1

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

>> [t,x] = ode45('eqx',[0 10],[0,1]);

dx

dy

2

 x

  ,

y

 x

  xy điều kiện đầu x(0)  0, (

y 0)  1

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) dt

dt

và cột thứ hai x2 = x(:,2). Với lệnh vẽ

Viết lại hệ dạng véc tơ như sau:

>> plot(t,x(:,1),t,x(:,2),'--'),xlabel('t'), dx

dx

>> axis([0 10 –1.12 1.12])

1

2

2

x  x, x  y 

 x

  x ,

 x

  x x

1

2

1

2

1

1 2

>> plot(x(:,1),x(:,2), 'k-'),

dt

dt

1

1

Để 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 0.5

0.8

0.6

function xdot = eqx(t,x);

0

0.4

xdot = zeros(2,1);

0.2

xdot(1) = -x(1)^2 + x(2);

- 0.5

0

xdot(2) = -x(1) - x(1)* x(2);

- 0.2

end % save with file name eqx.m

- 1

- 0.4 -0.6

- 0.4

- 0.2

0

0.2

0.4

0.6

0

2

4

6

8

10

t

Nguyen Quang Hoang

19

Nguyen Quang Hoang

20

Department of Applied Mechanics

Department of Applied Mechanics

5

16/09/2015

Giải hệ phương trình vi phân cấp 1

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

>> [t,x] = ode45('eqx2',[0 2*pi],[0,0]); nghiệm của phương trình vi phân sau

>> plot(t,x(:,1),t,x(:,2),'--'), xlabel('t'), y  16 y  sin(4.3 t)

>> axis([0 2*pi –3 3])

điều kiện đầu

(

y 0)  0, (

y 0)  0

Bằng cách đặt x  ,

y x  y

>> plot(x(:,1), x(:,2),'k-'),

1

2

ta nhận được hệ hai ptvp cấp một

>> xlabel('x_1'), ylabel('x_2‘)

x  x

1

2

x  f( t, )

x

3

x

x

3

x  sin(4.3 t)  16 x

1

2

2

1

2

2

Vế phải của hệ trên được thể hiện trong mộ

1

t m-file như sau:

1

0

function xdot = eqx2(t,x);

0

x 2

-1

xdot = zeros(2,1);

- 1

xdot(1) = x(2);

-2

- 2

xdot(2) = sin(4.3*t)–16*x(1);

-3

end

% save with file name eqx2.m

0

1

2

3

4

5

6

- 3 -0.6 -0.4 -0.2

0

0.2

0.4

0.6

t

x1

Nguyen Quang Hoang

21

Nguyen Quang Hoang

22

Department of Applied Mechanics

Department of Applied Mechanics

Giải hệ phương trình vi phân cấp 1

Giải hệ phương trình vi phân cấp 1

>> [t,x] = ode45('eqx2',[4*pi 20*pi],[0,0]); TT

Lệnh

Dạng bài

Độ chính

Trường hợp áp dụng

>> plot(t,x(:,1)),xlabel('t')

solver

toán

xác

1

không cứng

1

ode45

trung bình Hay được sử dụng, là sự lựa chọn (nonstiff)

đầu tiên (sử dụng công thức R-K45)

0.5

không cứng

Sử dụng khi không cần độ chính xác

0

2

ode23

thấp

cao, hoặc khi giải bài toán cứng vừa

(nonstiff)

phải. (sdụng công thức R-K23)

- 0.5

không cứng

3

ode113

thấp

đến (sử dụng công thức Adams-

- 1

(nonstiff)

cao

Bashforth-Moulton)

10

20

30

40

50

60

70

t

cứng

4

ode15s

thấp

đến sử dụng khi ode45 cho độ chính xác

(stiff)

trung bình không cao do bài toán cứng

Nguyen Quang Hoang

23

Nguyen Quang Hoang

24

Department of Applied Mechanics

Department of Applied Mechanics

6

16/09/2015

Giải hệ phương trình vi phân cấp 1

Ví dụ về hệ cứng (stiffness system). Ref. Chapter 5. Solving Diff. Equations TT

Lệnh

y  Ay  0,

(

y 0)  y

solver

Dạng bài toán Độ chính

xác

Trường hợp áp dụng

0

cứng

Khi độ chính xác không của hệ cứng

Ma trận A có các trị riêng phân biệt.

5

ode23s

thấp

và khi ma trận khối lượng là hằng

(stiff)

số

(

y t)

C v exp(

t)

i i

i

cứng vừa phải

Giải phương trình có độ cứng vừa

i

6

ode23t

(moderately

thấp

phải, nghiệm không giảm dần, có

Ví dụ

y

y ,

1

2

y

0

1 y

stiff)

thể giải được hệ DAEs

1

1

y

1001 y

1000 y

0

y

1000 1001

1000

1001 ,

y

y

y

y

2

2

cứng

2

1

2

7

ode23tb

thấp

giải hệ cứng với độ chính xác không

(stiff)

cao

0

1

1

A 

| A

I |

0

1000 1001

1000 1001

1,

1000





1

2

h

2 / max

Nguyen Quang Hoang

25

Nguyen Quang Hoang

26

Department of Applied Mechanics

Department of Applied Mechanics

Bài tập

Bài tập

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),





dy

(

q 0)

q ,

(

q 0)

q

eqbt7.m

0

0

eqbt1.m

 2

 .3 ,

y

(

y 0)  0,

t

 10,

t

  0.01

end

dt

Bằng cách hạ bậc ta nhận được

q



v



1



q  v v  M f

 Cv 

y   , y 

dx

dx

,

( ( t)

K )

q



1



v

M (f( t)  Cv  Kq)

1

2

2

 





eqbt2.m

 2 x x ,

 x

 ,

x (0)  x (0)  0,

t

 10



1 2

1

1

2

end

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

dx

dx











1

2

eqbt3.m

 x ,

 x

 ,

x (0)  x (0)  1,

t

 10

12

0

0

10

10

0

2

1

1

2

end

dt

dt









M  0 15

0 ,





C 

1

 0 20

1

 0 ,













dy

0

0

20

0

10

10









eqbt4.m

 t

 y  1, y(0)  1,

t

 10,

t

  0.01

end

dt

 2000

1

 000

0 

5 sin10 t









2

d y

dy

K 

1

 000 2000

1

 000 , f( t)





0





eqbt5.m

2

t

y

e



 

,

y(0)  2, y(0)  0,

t

 20









2

end

dt

dt

0

1

 000 1000

0









Nguyen Quang Hoang

27

Nguyen Quang Hoang

28

Department of Applied Mechanics

Department of Applied Mechanics

7

16/09/2015

Bài tập

Bài tập

function ydot = daodong3dof(t,y);

Phương trình vi phân chuyển động của tay máy 1 dof:

ydot = zeros(6,1);

2

2

2



 



 



( J

J r )

( K K r / R

br )

mgl cos( )  ( K r / R ) U

q = y(1:3); v = y(4:6);

O

m

m

e

a

m

a

Hãy viết các m-file thể hiện các phương trình vi

M = [12, 0, 0; 0, 15, 0; 0, 0, 20];

taymay1dof.m

phân và thực hiện việc giải số tìm chuyển động của

C = [10 -10 0; -10 20 -10; 0 -10 10];

tay máy:

K = [2, -1 0; -1 2 -1; 0 -1 1]*1000;

ft = [5*sin(10*t); 0; 0];

Tay máy

qdot = v;

vdot = inv(M)*(ft-C*v-K*q);

ydot = [qdot; vdot];

end

% save with file name daodong3dof.m

U

>> [t,y] = ode45(‘daodong3dof',[0 20],[0,0,0,0,0,0]); Hộp giảm tốc

>> plot(t,y(:,1)),xlabel('t')

Động cơ DC

Nguyen Quang Hoang

29

Nguyen Quang Hoang

30

Department of Applied Mechanics

Department of Applied Mechanics

Bài tập

Bài tập

Mq  Cq  Kq  f( t), Động lực học robot dạng chuỗi hở

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ệ

(

M )

q q  (

C ,

q )

q q  Dq  F sgn( )

q  (

g )

q  (

u t),

s

2. Xét dao động tự do không cản

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ệ

2. Cho q(t) tìm lực/mô men u(t)

•

Kiểm tra tính trực giao của ma trận dạng riêng thông qua M và K

•

Vẽ đồ thị q(t), q’(t), q’’(t) và u(t)

3. Xét dao động cưỡng bức có cản

3. Cho lực/mô men u(t) tìm chuyển động q(t)

•

Đồ thị biên độ tần số

•

Vẽ đồ thị q(t), q’(t), và u(t)

•

Tính dao động cưỡng bức [thành phần dao động cùng tần số

cưỡng bức]

(

u t)  (

g )

q  K (q  q )  K q p

0

d

•

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 ý]

Nguyen Quang Hoang

31

Nguyen Quang Hoang

32

Department of Applied Mechanics

Department of Applied Mechanics

8

16/09/2015

Bài tập

Động lực học robot song song

Bài 4. Cho một rôbốt song song phẳng ba C

q

3

bậc tự do, chuyển động trong mặt phẳng C

(

M q)q  (

C ,

q q)q  Dq  (

g q)  (

u )

T

t   (q),

1

ngang. Hãy thiết lập phương trình vi phân

(

 q)  0

chuyển động cho rôbốt với các tọa độ suy

p

3

C

2

p

2

T

rộng

q 

dư. Viết chương trình Matlab để mô

[ T

q , T

q , T

x ]

a

p

phỏng động lực học cho rôbốt. Biết rằng B

1

các tam giác ABC và (ABC) là các tam

2

 B 2

1. Thiết lập phương trình mô tả hệ

y

giác đều.

q

2

q

A2

1

A

p

1

q

[ q , q , q , p , p , p , x , y , ] T

2. Động học thuận/ngược robot song song

1

2

3

1

2

3

2

A

2

A

x

B

A

1

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) Hình bài 4

4. Cho lực/mô men u(t) tìm chuyển động bàn máy x(t) Nguyen Quang Hoang

33

Nguyen Quang Hoang

34

Department of Applied Mechanics

Department of Applied Mechanics

9

9/16/2015

Chương 8. Giới thiệu về Simulink

8.1 Khái niệm về simulink

Lecture 8

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 Giới thiệu về Simulink

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

Nguyen Q.Hoang

8.4 Đơn giản sơ đồ simulink

Department of Applied Mechanics

Đơn giản sơ đồ simulink bằng khối Fcn

Hanoi University of Science and Technology

Đơ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

1

2

Khái niệm về simulink

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 Xét phương trình vi phân sau

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 mv  bv  K (

u t)

1

 v 

[ Ku( t)  bv]

khi thiết lập các phương trình vi phân mô tả hệ, một sơ đồ khối tương m

ứng được đưa ra và từ đó sử dụng các khối của simulink để mô phỏng Biểu thức nghiệm tìm được dạng

hệ.

t

t 1

(

v t)  (

v 0) 

(

v ) d  (

v 0) 

[ K (

u )  bv] d





• Simulink là một phần mở rộng của MATLAB được phát triển bởi m

0

0

Mathworks Inc. Simulink là một gói chương trình để mô hình hóa, mô Biểu diễn dạng sơ đồ khối

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 v( t)

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 v(0)

có thể được lấy mẫu hoặc cập nhật số liệu ở các tốc độ khác nhau.

3

4

1

9/16/2015

Khái niệm về simulink

Các khối trong simulink

• Khởi động simulink

– Kích đúp vào biểu tượng

hoặc gõ lệnh

>> simulink

•

Biểu tượng

của Simulink

5

6

Các khối trong simulink

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 9

thang nhiệt độ Fahrenheit, phương trình mô tả

T  T  32

F

5 C

quan hệ nhiệt độ giữa hai thang đo này như sau

Các khối cần được kéo vào model

Vị trí các khối trong thư

window

viện Simulink

Ramp

Sources

Constant

Sources

Gain

Math Operation

Add hoặc Sum

Math Operation

Scope

Sinks

13

14

Một số ví dụ đơn giản

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

Mô phỏng một quá trình động học

Trong ví dụ này ta sẽ sử dụng Simulink để

x( t)  A sin( t

  )

mô tả một phương trình động học sau

Trong ví dụ này các thông số được nhận

A  2,   4,    / 3

các giá trị

Các khối cần được kéo vào

Vị trí các khối trong

model window

thư viện Simulink

Constant

Sources

Ramp

Sources

Gain

Math Operation

Sum

Math Operation

Product

Math Operation





Trigonometry

Math Operation

x( t) A sin( t



)

Function

A  2,   4,    / 3

Mux

Signal Routing

Scope

Sinks

17

18

Mô phỏng một quá trình động học

Giải phương trình vi phân

Xét phương trình vi phân mô tả một hệ

1

động lực liên tục như sau

 v 

[ Ku( t)  bv]

m

u( t)

v

m

bv

x( t)  2

 .5 x( t)  (

u t)

x( t)  A cos( t

  )

A  2,   4,    / 3

19

20

5

9/16/2015

Mô tả hệ dao động một bậc tự do

Mô tả hệ dao động một bậc tự do

x

Phương trình vi phân chuyển động

x

Phương trình vi phân chuyển động

mx( t)

f ( t)

 f( t)

c

mx( t)  cx( t)  kx( t)  f ( t) f ( t)

m

Viết lại phương trình vi phân chuyển

m

Viết lại phương trình vi phân chuyển

động trên dưới dạng sau

k

động trên dưới dạng sau

1

Hình 8-12.

1

x 

[ f ( t)]

x 

[ f ( t)  cx  kx ]

m

m

x

1

x

1

x

f ( t)

1

x

1

x

1

x

f ( t)



1

m

s

s

m

s

s

c



k



21

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

Mô hình hóa và mô phỏng hệ cơ điện

mx( t)  cx( t)  kx( t)  f ( t) 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ở, R

Tay máy

a

•

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

U

Các thông số tay máy:

Hộp giảm tốc

•

khối lượng m, khối tâm C, OC = L.

DC động 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 b

2

Hộp giảm

>> load ketqua.mat

tốc:

•

tỷ số truyền của

Tương tác cơ – điện:

t = ans.time(:,1)

u = ans.data(:,1)

r   /  r / r  1

1

2

2

1

V

 K ,

M  K i

emf

e

m

x = ans.data(:,2)

plot(t,u),

25

26

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

Đối với hộp giảm tốc (không khối lượng, không g

ma sát) ta có công suất vào bằng công suất ra

C1

T   M 

lo d

a

1

M 1 M 1



Ra

La

r   r 

   r

 / r  r

O

1

2

2

1

 T

r

load

2

T

 M  /  M / r

U

V

r

lo d

a

1

1

emf

1

g



T

load

Động cơ DC

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) J   M  b   mgl cos() M

1



di

O

1

2

U  U  U  V



L

 R i  U  K 

i  q

O

R

L

emf

a

a

e

dt

e

Phương trình cơ (phương trình chuyển động quay của rôto) Khử T

& M

load

1

d L  m

 ( F )

 J   K i  b   T ,

  

m

m

1

load

z

z

k

dt

27

28

7

9/16/2015

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 PTVP mô tả

Hệ 3 phương trình vi phân cấp 1 mô tả qua hệ đầu vào U và đầu ra 

J   K i  b   T ,

  

m

m

1

load

di

L

 U  R i  K rq

a

a

e

dt

di

L

 R i  U  K ,



i  q

2

2

( J  J r ) q  K ri  ( b r  b ) q  mgl cos() a

a

e

e

dt

O

m

m

1

2

  q

J   M  b   mgl cos() O

1

2

Nếu bỏ qua L di/dt  0

a

2

Tay máy

Các ràng buộc

( J  J r ) 

O

m

2

2









( K K r / R

b r

b )

   r

Khử T

& M

load

1

m

e

a

1

2

m

 g c

l os()  ( K r / R ) U

m

a

T   M 

lo d

a

1

U

U

,

 

Hộp giảm tốc

T

 M  /   M / r

Hệ thống

d

loa

1

1

DC động cơ

29

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

Vị trí các khối trong thư viện

window

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

Mô hình hóa và mô phỏng hệ cơ điện

% Can nhap thong so he:

Bài tập: Ghép 2 tay máy 1dof thành 1 tay

máy phẳng 2dof.

% - DC motor: Ra, La, Ke, Km, Jm, b1

% - Tay may: m, lc, g, JO, b2

R , K , K , L , b , J

 , ,

2a

2 m

2e

2a

21

2 T

m ay máy

% - hop giam toc: ti so truyen r

U ,

1

Hệ thống

1

1

r

2

U

 , ,

2

2

2

m , J , l , b

2

C 2 C 2

22

% - DC motor:

Ra=1, La=0.0001,

Khâu 2

Ke=1, Km=1,

(Link2)

Jm=0.001, b1=0.01;

UU 2

Hộp giảm tốc

% - Tay may:

Tay máy

DC động cơ

m=5, lc=0.5, g=9.81,

JO=0.2, b2=0.2;

% -hop giam toc: ti so truyen

R , K , K , L , b , J

1a

1 m

1e

1a

11

1 m

r = 2.0

r 1

Khâu 1

m , J , l , b

U U

(Link1)

1

C 1 C 1

12

1

Hộp giảm tốc

DC động cơ

33

34

Mô hình hóa và mô phỏng hệ cơ điện

2

( J  J r ) 

O

m

2

2

( K K r / R  b r  b )

m

e

a

1

2

m

 g c

l os()  ( K r / R ) U

Tay máy

m

a

U

,

 

Hệ thống

U

Điều khiển

Hộp giảm tốc

DC động cơ



Xét hai bộ điều khiển

0

1

1. U

( K r / R )

k (

) k

PD-controller

m

a

p

0

d

1

2. U

( K r / R )

k (

) k

mgl cos

PD-controller + bù

m

a

p

0

d

trọng lực

35

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 1

L

 R i  R ( i  i )  E( t) 1 1

2

1

2

L

dt

R

i

i

1

2

1

di

q

2

2

E(t)

L

 R ( i  i ) 

 0

2

1

2

R

2

dt

C

C

i

1

dq

i

2

L

2

 i 2

dt

Cho biết các số liệu

R  4 ,

 R  10 ,

 L  0.032 H, C  0.53 F

1

2

và điện áp đặt lên mạch có dạng

2

 0 V khi 0  t  0.005 s Hãy giải và vẽ ra đồ thị các dòng điện



E( t)  0,

khi t  0.005 s

i ( t), i ( t)

t  0...0.05 s



1

2

39

40

10

9/16/2015

Bài tập

Đơn giản hoá sơ đồ simulink

R



1. Đơn giản sơ đồ simulink bằng khối Fcn

1

d i

i

2. Đơn giản sơ đồ simulink bằng khối subsystem

dt 1

1

1

E( t)



1

L

s

3. Kết hợp simulink và script file (m-file)

-

R 2

+

d i

dt 2

i

-

1

1

2

+

L

s

di

L 1  E t

( )  R i  R i

(  i )

dt

1 1

2

1

2

d q

1

q

di

q

dt

2

2

L 2  R i

(  i

2

) 

dt

2

1

2

C

s

dq 2  i

C 1





dt

2

41

42

Đơn giản hoá sơ đồ simulink

Đơ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  5sin(4 t   /3)

43

44

11

9/16/2015

Đơn giản hoá sơ đồ simulink

Đơn giản hoá sơ đồ simulink

2. Đơn giản sơ đồ simulink bằng khối subsystem

di / dt

di

1



 U  R i  K rq

a

e



dt

La

dq

1

2



K



ri  ( b r  b ) q  mgl cos()

2

m

1

2

dt

( J

J r ) 





O

m

d  q

dt

M 1



45

46

Đơn giản hoá sơ đồ simulink

Đơn giản hoá sơ đồ simulink

47

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

E

M

2

g

q

C

2

2

A

q  v

C

1

1

v 

(

M q) [u  (

C ,

q )

v v  (

g q)]

M

1

q

1

O

y  f( ,

y )

u ,

(

y 0)  y0

 M 

1

u  M





2





49

50

Đơn giản hoá sơ đồ simulink

Đơn giản hoá sơ đồ simulink

3. Kết hợp simulink và script file (m-file)

4. Xử lý kết quả

load sim_results.mat

t=qv(1,:);

q1=qv(2,:);

q2=qv(3,:);

q1_dot=qv(4,:);

q2_dot=qv(5,:);

Soạn một MATLAB function [m-file]

figure(1)

plot(t,q1,'-k','linewidth',2), grid on

function ydot = xvdot(in)

xlabel ('t [s]');

ylabel ('q_1 [rad]');

figure(2)

plot(t,q2,'-k','linewidth',2), grid on

xlabel ('t [s]');

ylabel ('q_2 [rad]');

ydot = …. ;

% save ve_q_t.m

end

51

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

(

M )

q q  (

C ,

q )

q q  (

g )

q  Bu

x

q  v

ux

1

v 

(

M q) [u  (

C ,

q )

v v  (

g q)]

O

mt

l

ul

m

m

m l cos

(

M q)

t

p

p

2

x

z



m l cos

m l

p

p

Con lắc elliptic

m

u

p

x

0

m l sin

p

O

(

C ,

q q)

0

0

mt

ul

l

0

(

g q)

1

z



m gl sin

B

p

0

mp

53

54

Đơn giản hoá sơ đồ simulink

Đơn giản hoá sơ đồ simulink

function xdot = xdot_elliptic_pendul(ins)

Kết hợp simulink-model với lệnh sim

•

Thực hiện mô phỏng mà không cần mở simulink

% cac tham so he

M = [m1+m2, m2*l*cos(q(2));

•

Sử dụng được mô hình simulink trong một m-file

m1 = 5; % kg

m2*l*cos(q(2)), m2*l^2 ];

•

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 m2 = 2; % kg

l = 0.5; % m

C = [0, -m2*l*qdot(2)*sin(q(2));

khiển bằng GA, PSO, Nelder-Mead,…,

gr = 9.81; % m/s^2

0, 0];

•

Cú pháp: xem >> help sim

d1 = 6.2; d2 = 0; % damping coefs

D = diag([d1, d2]); % damping matrix

% ins = [1 2 3 4 5]'; % for test

% viet ra cac phuong trinh vi phan

g = [0; m2*l*gr*sin(q(2))];

q = ins(1:2);

% bien chua toa do suy rong

B = [1; 0]; % he hut dan dong

qdot = ins(3:4);

% B = [1, 0; 0, 1]; % he đủ dan dong

% bien chua van toc suy rong

u = ins(5);

vdot = inv(M)*(B*u- (C+D)*qdot - g);

% bien vao (bien dieu khien)

5. Mask and unmask ??

xdot = [qdot; vdot];

55

56

14

9/16/2015

Đơn giản hoá sơ đồ simulink

Đơ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)

57

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

Điều khiển phản hồi trạng thái

 x = f ( x, u), x(0) = x

u = u( x) 0

y = g( x) Các tiêu chuẩn tối ưu cần được cực tiểu có thể là t

t

f

f

T

J =

e Qedt =

( x - x ) TQ( x - x ) dt ò

ò

(1)

1

0

0

d

d

tf

J =

( y - y ) TQ( y - y ) dt ò

(2)

2

0

d

d

tf

T

J =

u Rudt

ò

(3)

3

0

tf

J =

1 dt

ò

(4)

4

0

T

Q

Q

0,

T

=

³

R = R > 0

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,

t

J

(1 l)

f

J

J

l

(é

) T (

)

T

ù

=

-

+

=

-

-

+

dt

ò x x Q x x

u Ru

(5)

1

3

êë

ú

0

d

d

û

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à T

x ( t ) Fx( t ), T

F = F > 0

(6)

f

f

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à t

T

J

x ( t ) Fx( t ) f

(é x

x ) TQ( x

x )

T

ù

=

+

-

-

+ u Ru dt

ò

(7)

f

f

êë

ú

0

d

d

û

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(0) = x 0

y = Cx

Tiêu chuẩn 2 trở thành

t

t

t

f

f

f

T

T

T

T

=

y y

=

x C Cx

=

x Qx ,

T

J

dt

dt

dt Q =

ò

ò

ò

C C

2

0

0

0

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

 x = f ( x, u), x(0) = x 0

Và tiêu chuẩn tối ưu sau cần được cực tiểu (hoặc cực đại) t

J = G[ x( t )]

f

+

F( x, u) dt ò

f

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 u

£ u £ u ,

i = 1,2,..., m

i, m

i

i, M

dim( x) = dim( f ) = , n

dim( u) = 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:

 x ( t) = f ( x ( t), u ( t)), 0

0

0

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à

u( t) = u ( t) + du( t), 0

Tiêu chuẩn tối ưu thay đổi một lượng là

J

d = J[ u ( t) + du( t)]- J[ u ( t)]

0

0

t

G

é F

F

ù

é

G

ù

¶

¶

¶

¶

=

dx( t )

f

+

ê

dx +

du ú dt + F

ê ( t ) +

f ( t )ú t

d

ò

f

0

f

f

f

ê

ú

ê

ú

¶ x

¶ x

¶ u

¶

ë

û

ë

x

û

(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):

 x ( t) + d  x( t) = f ( x ( t) + dx( t), u ( t) + du( t)), 0

0

0

d

¶ f

¶ f

d  x( t) =

dx =

dx( t) +

du( t),

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 t é

ù

t

d

é

ù

x

¶ f

¶

f

d

f

T

( )

ê l ( t)

f

T

ú dt -

l ( t)ê

dx( t) +

du( t)ú dt = 0

ò

ò

(9)

0

ê

ú

0

dt

ê

ú

¶ 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

d = F

ê ( t ) +

f ( t )

T

ú t

d + l (0) dx(0) T

+ ê

- l ( t )ú dx( t ) f

f

f

f

f

ê

ú

ê

ú

¶ x

¶

ë

û

ë x

û

é

T

t

ù

¶

¶

l

f

H

H

d

( t)

ê

dx( t)

du( t)

d ú

+

+

+

x dt

ò ê

ú

0

(10)

¶ x

¶ u

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

= ( x, u)

T

H

F

+ l ( t) f ( x, u) (11)

Do (t) là tùy ý, nên ta chọn nó thỏa mãn

T

dl ( t)

H

¶

= -

( x, u, l) (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

æ G ö

¶

T

l ( t )

ç

÷

= ç

÷

f

ç

÷

ç

(13)

è ¶ x ÷ø t t= 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

t é

ù

¶

f

H

J

d =

ê

d ( t)ú dt

ò

u

(14)

0

ê

ú

¶

ë u

û

Cực tiểu hàm J

Nếu u

J

d ³ 0 ,

" du( t)

0(t) là nghiệm tối ưu (J đạt cực tiểu) thì

Trường hợp 1. Không có ràng buộc đặt trên u(t), thì H

¶

= 0,

" t

Đây là điều kiện cần để tối ưu J.

¶ u

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

¶

u = u ,

£ 0,

&

u = u ,

³ 0

i,0

i, M

i,0

i, m

¶ 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

t é

ù

¶

f

H

J

d =

ê

d ( t)ú dt

ò

u

(14)

0

ê

ú

¶

ë u

û

Cực đại hàm J

Nếu u

J

d £ 0 ,

" du( t)

0(t) là nghiệm tối ưu (J đạt cực đại) thì

Trường hợp 1. Không có ràng buộc đặt trên u(t), thì H

¶

= 0,

" t

Đây là điều kiện cần để tối ưu J.

¶ u

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

¶

u = u ,

³ 0,

&

u = u ,

£ 0

i,0

i, M

i,0

i, m

¶ 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 = f ( x, u), dH

H

¶ dx

H

¶ du

¶ H dl

H = H( x, u, l) 

=

+

+

(15)

dt

¶ x dt

¶ u dt

¶ l dt

H

¶

T

dl ( t)

¶ H

Nhưng

T

= f ( x, u) và theo cách chọn

= -

( x, u, l)

¶ l

dt

¶ x

dH

¶ H du

nên

=

dt

¶ u dt

(16)

Trường hợp 1. Không có ràng

Trường hợp 2. Khi có ràng buộc đặt trên

buộc đặt trên u(t), thì

u(t), thì u = umax hoặc u = umin nên

¶ H = 0

du =

¶

0

u

dt

dH

Do đó, đối với điều khiển tối ưu ta có

= 0 hay H = const.

(16)

dt

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

d = F

ê ( t ) +

f ( t )

T

ú t

d + l (0) dx(0) T

+ ê

- l ( t )ú dx( t ) f

f

f

f

f

ê

ú

ê

ú

¶ x

¶

ë

û

ë x

û

é

T

t

ù

¶

¶

l

f

H

H

d

( t)

ê

dx( t)

du( t)

d ú

+

+

+

x dt

ò ê

ú

0

¶ x

¶ u

dt

(10)

êë

úû

Để J không phụ thuộc tf, ta cho

G

¶

F( t ) +

f ( t ) = 0

(17)

f

f

¶ x

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 đó

H = F( x, u) T

+ l ( t) f ( x, u) H( t ) = F( t )

T

+ l ( t ) f ( t ) = 0

(18)

f

f

f

f

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ó

dx( t ) = x( t ) - x ( t ) = - f ( t ) t d

(19)

f

f

0

f

f

f

Với các biểu thức trên, ta có

t é

H

ù

¶

J

d = F

é ( t )

T

+ l ( t ) f ( t )ù T

t

d + l (0) dx(0) f

+

ê

du( t)ú dt (20)

ê

ò

f

f

f

ë

ú f

û

0

ê

ú

¶

ë u

û

Để loại bỏ số hạng liên quan điều kiện biên khi tf không xác định H( t ) = F( t )

T

+ l ( t ) f ( t ) = 0

f

f

f

f

Vì H(t) = const, nên H(t) = H(tf) = 0 dọ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

 x = Ax + Bu,

x(0) = x 0

(21)

Và tiêu chuẩn tối ưu sau cần được cực tiểu

1

1

t

T

J

x ( t ) Fx( t ) f é T

T

ù

=

+

x Qx + u Ru dt

ò

(22)

f

f

êë

ú

0

2

2

û

T

F

F

0,

T

Q

Q

0,

T

=

³

=

³

R = R > 0

u

£ u £ u ,

i = 1,2,..., n

i, m

i

i, M

Hàm Hamilton trong bài toán này là

1

= { T

T

x Qx + u Ru}

T

H

+ l ( t)( Ax( t) + Bu( t )

)

(23)

2

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 H

¶

= 0

T



Ru + B

(

l t) = 0

¶ u

-1

T



u = - R B

(

l t)

(24)

Động học của (t) đưa ra từ (12) cùng với điều kiện cuối (13) d (

l t)

T

= - Qx - A l( t),

l( t ) = Fx( t ) (25)

f

f

dt

Kết hợp (21)-(24)-(25) ta được

-1

 x( t)

T

= Ax - BR B (

l t), x(0) = x 0



l( t)

T

= - Qx - A l( t),

l( t ) = Fx( t ) f

f

é  ù

é

-1

x( t)

T ù é

A

BR B

x( t)ù

é x( t)ù

-

ê

ú

ê

ú ê

ú

ê

ú

(26)

=

=

ê 

M

(

l t)ú

ê

T

ú ê

Q

A

(

l t)ú

ê (

l 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

Xác định điều kiện đầu (0) ta viết

é x( t)ù

é x(0)ù

é

ù

ê

ú

E ( t) E ( t) t

eM ê

ú

=

ê

ú

ê

Mt

11

12

e

=

l( t)ú

ê (

l 0)ú

ê

ú

êë

úû

êë

úû

E ( t) E ( t) ê 21

22

ë

úû

Suy ra

x( t) = E ( t) x(0) + E ( t) (

l 0)

11

12

(27)

(

l t) = E ( t) x(0) + E ( t) (

l 0)

21

22

Và từ điều kiện

(

l t ) = Fx( t ) f

f

E ( t ) x(0) E ( t ) (

l 0)

é

F E ( t ) x(0) E ( t ) (

l 0)ù



+

=

+

21

f

22

f

ê 11 f

12

f

ë

úû

é E ( t ) FE ( t )ù (

l 0)

é FE ( t ) E ( t )ù



-

=

-

x(0)

ê 22 f

12

f

ë

úû

ê

11

f

21

f

ë

úû

-1

(0

l )

é E ( t ) FE ( t )ù é FE ( t ) E ( t )ù



=

-

-

x(0)

ê 22 f

12

f

ë

úû ê

11

f

21

f

ë

úû

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

l( ) = S( ) x( ), T

t

t

t

S = S

(28)

Thay vào (25) ta được

d (

l t)

T

= 

Sx + 

Sx = - Qx - A Sx, (29)

dt

Cùng với (21, 24) ta được

(

-1

T

T



S + SA - SBR B S + Q + A S) x = 0

(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

T

T



S + SA - SBR B S + Q + A S = 0

(31)

Từ (13) và (28) cho ta

æ G ö

¶

T

l ( t )

ç

÷

= ç

÷

= Fx( t ), l( t ) = S( t ) x( t )  S( t ) = F

(32)

f

ç

÷

f

f

f

f

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 B S( t) x( t) = - K( t) x( t), (33)

1

K( t)

-

T

= R B S( 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

x = Ax + B ,

u

( ,

A )

B

là điều khiển được

Thuật giải

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

-

T

T

S + SA - SBR B S + Q + A S = 0, (

S t ) = F

f

3. Xác định ma trận K(t) từ đó xác định véc tơ điều khiển 1

(

K )

-

T

t = R B (

S t),

u = - (

K t) (

x t)

4. Đánh giá (xác định) giá trị hàm mục tiêu:

1

1

t

T

J

x ( t )F (

x t )

f é T

T

ù

=

+

x Qx + u Ru dt

ò

f

f

êë

ú

0

2

2

û

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: t  

f

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

-

T

T

SA - SBR B S + Q + A S = 0, do

S = 0,

3. Xác định ma trận K(t) từ đó xác định véc tơ điều khiển 1

-

T

K = R B ,

S

u = -K (

x t)

4. Kiểm tra tính ổn định của hệ kín

x = (A - B )

K x

5. Đánh giá (xác định) giá trị hàm mục tiêu:

Lệnh Matlab:

1

1

t

T

J

x ( t )F (

x t )

f é T

T

ù

=

+

x Qx + u Ru dt

ò

>> K = lqr(A,B,Q,R)

f

f

êë

ú

0

2

2

û

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

t

T

J

x ( t )F (

x t )

f é T

T

x Qx

u Ru

2 T

ù

=

+

+

+ x Nu dt

ò

f

f

êë

ú

0

2

2

û

Khi đó sử dụng phép biến đổi

T

x Qx + 2 T

T

T

T

x Nu + u Ru = x Q x + v Rv m

Với

1

-

T

-1

Q = Q - NR N

&

T

v = u + R N x

m

Khi đó hàm mục tiêu trở thành

1

1

t

T

J

x ( t )F (

x t )

f é T

T

ù

=

+

x Q x + v Rv dt

ò

f

f

êë

ú

0

2

2

m

û

Hệ ban đầu được biến đổi tương ứng thành

-1

x = A x + B ,

v

with

T

A = A - BR N ,

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 Q



m và Am (kiểm tra ma trận Qm

0, bán x.đ.d)

1

-

T

1

Q = Q - NR N ,

-

T

A = A - BR N

m

m

3. Giải phương trình Riccati (31), hệ ptvp thỏa mãn điều kiện cuối

-1 T

T

S + S A - S BR B S + Q + A S = 0, S ( t ) = F

m

m

m

m

m

m

m m

m

f

4. Xác định ma trận Km(t) từ đó xác định véc tơ điều khiển 1

K ( t)

- é T

R

B S ( t)

T ù

=

+ N ,

Lệnh Matlab:

m

ê

m

ë

úû

u = -K ( t) (

x t)

>> K = lqr(A,B,Q,R)

m

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

t

x = x + ,

u

2

f

2

2

J = x( t ) +

(2 x + 3 u ) dt

ò

f

0

2

2

Với

A = 1,

b = 1, F = 5,

Q = 2,

R = 3

Phương trình Riccati trở thành

1 2

S = -2S + S - 2,

S( t ) = F = 5

3

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

u( t) = - S( t) x( t), 0 £ t £ t 3

f

Nguyen Quang Hoang

25

Department of Applied Mechanics

Ví dụ 1

Nếu cho t 

f

5

1

¥

2

2

2

x = x + ,

u

J = x( t ) +

(2 x + 3 u ) dt

ò

f

0

2

2

Với

A = 1,

b = 1, F = 5,

Q = 2,

R = 3

Phương trình đại số Riccati trở thành

1 2

0 = -2S + S - 2  S

= 3  15

1,2

3

Để hệ kín ổn định ta chọn

+

S = S = 3 + 15

3

15

K =

1

3

Luật điều khiển tối ưu sẽ là

u = K

- x,

15

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:

é x ù

é0 1ù é x ù é0ù

é x ù

ê 1 ú

ê

ú ê 1 ú

ê ú

é

ù ê 1

,

u

y

1 0

ú

=

+

=

x

ê ú

ê



0 0ú ê x ú

ê1ú

êë

úû x ê ú

ê 2

ë úû

êë

úû ê 2

ë úû

êë úû

ê 2

ë úû

Tìm luật điều khiển

u

=

-

k

( x

-

r

)

-

k

x để duy trì x = r, x = 0

1

1

2 2

1

2

cực tiểu chỉ tiêu

1

¥

2

2

J =

(( x - r) + u ) dt

ò

(tức cho t )

1

f

0

2

Đặt biến mới

x = x - r

é x ù

é0 1ù é x ù é0ù

1

¥

1

1

ê 1 ú

ê

ú ê 1 ú

ê ú

=

+

,

u

2

2

J =

( x + u ) dt

ò 

x = x

x

ê  ú

ê



0 0ú ê x ú

ê



1ú

1

0

2

2

2

ê 2

ë úû

êë

úû ê 2

ë úû

êë úû

é2 0ù

ê

ú

Q =

; R = 2

ê0 0ú

Phương trình đại số Riccati trở thành

êë

úû

1

-

T

T

SA - SBR B S + Q + A S = 0

Nguyen Quang Hoang

27

Department of Applied Mechanics

Ví dụ 2

Giải phương trình đại số Riccati

é0 0ù S

é

S ù

é S

S ù é0 1ù

é S

S ù é0ù é1ù

S

é

S ù

é2 0ù

ê

ú ê 11

12 ú

ê 11

12 ú ê

ú

ê 11

12 ú ê ú

+

-

ê ú é

ù ê 11

12

0 1

ú

ê

ú

+

ê1 0ú S

ê

S ú

S

ê

S ú ê0 0ú

S

ê

S ú ê1ú ê2ú êë

úû S ê

S ú

ê0 0ú

êë

úû ê 12

22

ë

úû ê 12

22

ë

úû êë

úû ê 12

22

ë

úû êë úû ë û

ê 12

22

ë

úû êë

úû

é0 0ù

ê

ú

= ê0 0ú

êë

úû

2

2

S

S S

S

12

12 22

22

-

+ 2 = 0,

S -

= 0, -

+ 2 S = 0

11

12

2

2

2

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

2

é 2 2 ù

ê

ú

S = ê

ú

ê 2

2 2ú

ë

û

-1 T

u = -R B Sx = x

- - 2 x = (

- x - r) - 2 x

1

2

1

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

1

¥

2

2

J =

( x + u ) dt,

x = x - r

ò

1

1

1

0

2

Hệ và bài toán tối ưu mô tả bởi:

é0 1ù

é0ù

é2 0ù

ê

ú

A

,

ê ú

B

,

ê

ú

=

=

Q =

; R = 2

ê0 0ú

ê1ú

ê0 0ú

êë

úû

êë úû

êë

úû

Phương trình đại số Riccati

1

-

T

T

SA - SBR B S + Q + A S = 0

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

closed_Syst = ss(A - B*K, B*K(1,1), C, 0)

A = [0 1; 0 0];

t = 0:0.1:10;

B = [0; 1];

r = 2*ones (size(t));

C = [1 0];

[y, t, x] = lsim (closed_Syst, r, t, [5 0]);

Q = [2 0; 0 0];

u = -K*x' + K(1,1)*r;

R = [2];

plot (t, y, 'k', 'linewidth', 2), grid on, hold on

[K, S, e] = lqr (A, B, Q, R)

plot(t,u, 'r', 'linewidth', 2)

legend('x_1(t)', 'u(t)')

5

K =

x (t)

4

1

1.0000 1.4142

u(t)

3

S =

2

2.8284 2.0000

1

2.0000 2.8284

e =

0

-0.7071 + 0.7071i

-1

-0.7071 - 0.7071i

-2

-3

e là trị riêng của ma trận (A - BK)

0

2

4

6

8

10

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 ù

é 0

1 ù é x ù

é0ù

ê 1 ú

ê

ú ê 1 ú

ê ú

=

+

,

u

x

ê ú

ê



1

2ú x

ê ú

ê1ú

-

-

ê 2

ë úû

êë

úû ê 2

ë úû

êë úû

5

1

2

2

2

2

2

J = ´20 x (5) +

( x + 2 x + x x + u ) dt, t = 5

ò

1

1

2

1 2

0

2

2

f

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

T

x Qx + 2 T

T

T

T

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

1

-

T

v = u + R N x

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 Q 

m

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

02/12/2015

Tổng quan về tính ổn định

ỔN ĐỊNH

 Đặc điểm kỹ thuật hết sức quan trọng trong thiết kế điều khiển

 Hệ thống không ổn định phải được ổn định bởi thông tin phản hồi.

 Hệ thống vòng kín không ổn định là vô ích.

 Điều gì xảy ra nếu một hệ thống không ổn định?

• có thể đạt các giới hạn cơ khí / hay bão hòa điện

• có thể phá vỡ hoặc gây cháy máy, ...

Nguyen Quang Hoang

3

Department of Applied Mechanics

Tổng quan về tính ổn định

Xét quả cầu nhỏ đặt tại bốn vị trí như trên hình vẽ. Nếu bỏ qua lực ma sát, một sự thay đổi nhỏ từ trạng thái cân bằng dẫn đến những điều sau đây: 1. Chuyển động dao động quanh trạng thái cân bằng không phụ thuộc vào độ lớn của nhiễu loạn ban đầu. Trong trường hợp này, vị trí cân bằng ổn định toàn cục.

2. Vị trí mới cũng là một trạng thái cân bằng, và quả bóng nhỏ sẽ ở lại trong nó.

Đây là vị trí cân bằng ổn định (phiếm định).

3. Quả bóng di chuyển đi mà không trở về trạng thái cân bằng. Vị trí cân bằng này là không ổn định.

4. Chuyển động dao động về trạng thái cân bằng, trừ khi các nhiễu kích động ban đầu quá lớn làm quả cầu buộc phải chuyển đến vị trí cân bằng mới. Vị trí cân bằng này là ổn định cục bộ.

(1)

(2)

(3)

(4)

Nguyen Quang Hoang

4

Department of Applied Mechanics

2

02/12/2015

Tổng quan về tính ổn định

Một tiêu chuẩn quan trọng để đánh giá một hệ động lực là đáp ứng ổn định của nó. Tổng quan về tính ổn định được thể hiện như sau: Một hệ tuyến tính là ổn định nếu đáp ứng đầu ra bị chặn khi tác động đầu là bị chặn (ổn định BIBO), hoặc nó trở về trạng thái gốc ban đầu (ổn định tiệm cận). Nếu xuất hiện đáp ứng khác với hai điều nêu trên thì hệ tuyến tính là không ổn định.

Kích động ngoài bao gồm cả điều kiện đầu và tác động điều khiển, cũng như nhiễu tác động.

Nguyen Quang Hoang

5

Department of Applied Mechanics

Ổn định BIBO

u

y

Hệ thống

Hệ tuyến tính LTI

D

x = Ax + B ,

u

u

x

x

y

y = Cx + D ,

u

1

B

C

s

n

x Î Â ,

m

u Î Â ,

l

y Î Â

A

Hệ phi tuyến

y

y = (

g ,

x u)

x = f( ,

x u),

y = (

g ,

x u),

x

1

x

u

n

x Î Â ,

m

u Î Â ,

l

y Î Â

x = f( ,

x u)

s

Nguyen Quang Hoang

6

Department of Applied Mechanics

3

02/12/2015

Ổn định BIBO

Định nghĩa ổn định BIBO: Một hệ thống được gọi là ổn định BIBO

(Bounded Input – Bounded Output], nếu khi kích động đầu vào u(t) bị

chặn thì đáp ứng của hệ (đầu ra) y(t) cũng bị chặn.

Nếu với mọi số dương N, M. Mọi đầu vào bị chặn sẽ tạo ra đầu ra tương ứng cũng bị chặn.

(

u t) £ N < ¥  (

y t) £ M < ¥

ICs =0

u( t)

u

Hệ ổn

y

y( t)

định BIBO

Nguyen Quang Hoang

7

Department of Applied Mechanics

Ổn định BIBO

Đối với hệ tuyến tính (SISO)

x = Ax + b ,

n

u

x Î Â , u Î Â

y =

,

cx

y Î Â

Tính được hàm truyền

(

p s)

1

G( s)

(

c sI

)-

=

=

- A

b

q( s)

BIBO Stability ⇔ Tất cả các điểm cực của hàm truyền nằm bên trái trục ảo trong mặt phẳng phức.

Đưa ra phương trình đặc trưng, sau đó giải tìm nghiệm của phương trình đặc trưng.

q( s)  0

Nghiệm của phương trình đặc trưng chính là các điểm cực của hàm truyền G(s).

Nguyen Quang Hoang

8

Department of Applied Mechanics

4

02/12/2015

Ổn định tiệm cận hệ tuyến tính

Ổn định tiệm cận [Asymptotic stablility]

ICs

u( t) = 0

Hệ ổn định

y

Khi u(t) = 0, hệ trở thành

tiệm cận

x = A ,

x

y = c ,

x

y( t)

(

x t)  0 khi t 

,

¥ hay y( t)  0

Đối với hệ tuyến tính :

x = Ax + B ,

u

t

y = C ,

x

• Hệ ổn định tiệm cận ⇔ tất cả các trị riêng của ma trận A có phần thực âm (nằm bên trái trục ảo)

(

P s)

C

I -

-

adj[ s

A]B

1

(

G s) =

= (

C sI - )

A

B =

q( s)

sI - A

sI - A = 0

Giải tìm các trị riêng của A.

Lưu ý: Ổn định tiệm cận không phụ thuộc ma trận B và C

Nguyen Quang Hoang

9

Department of Applied Mechanics

Ổn định tiệm cận hệ tuyến tính

Ổn định tiệm cận từ việc tách mô hình

x = Ax + B ,

u

Giả sử A có tất cả các trị riêng phân biệt y = C ,

x

(chéo hóa được A).

l , v là cặp trị riêng và véc tơ riêng của A Av = l v , i = 1 ,  ,n i

i

i

i i

Xây dựng ma trận chuyển đổi

T = [v , v , ,

 v ]

1

2

n

Biến đổi tọa độ

x = T x  x = T 



x

x = Ax + Bu

-1

1

-

 x = T AT x + T Bu x = Ax + Bu

y = Cx

y = CT x

y = Cx + Du

x = Tz

 x = A x + B , u

1

-

-1

z = T ATz + T Bu

-1

A = T A ,

T

y = C x

y = CTz + Du

1

-

B = T

,

B

C = CT

Nguyen Quang Hoang

10

Department of Applied Mechanics

5

02/12/2015

Ổn định tiệm cận hệ tuyến tính

Lưu ý rằng

l

é

0

0 0 ù

ê 1

ú

ê 0 l

0

0 ú

-1

ê

2

ú

A = T AT = ê  0

ú

ê



 ú

ê 0 0  l ú

ê

n ú

ë

û

Vậy nên, khi u = 0, ta có



é

 x = l x , i = 1,2,..., n x ù

l

é

0

0 0 ù é x ù

i

i i

ê 1 ú

ê 1

ú ê 1 ú

ê x ú ê 0 l

0

0 ú ê x ú

 x ( t) = x (0)exp( l t)



i

i

i

x = A x

ê 2 ú

ê

2

ú ê 2 ú



=

ê ú

ê





0

ú ê ú

ê ú

ê



 ú ê  ú

Với điều kiện đầu

ê x ú ê 0 0  l ú ê x ú ê

1

(0

x )

-

= T

(

x 0)

n ú

ê

n ú ê n ú

ë û

ë

û ë û

Tìm được x(t)

l t

l t

l t

1

2

(

x t) = T (

x t) = v x (0) e + v x (0) e +  + v x (0) n e ,

1 1

2 2

n n

Như vậy,

hệ ổn định tiệm cận ⇔ tất cả các trị riêng của ma trận A nằm bên trái trục ảo Nguyen Quang Hoang

11

Department of Applied Mechanics

Hệ tuyến tính: Ổn định tiệm cận và ổn định BIBO

Trong trường hợp không có sự loại trừ điểm không và điểm cực (hay không có điểm không nào trùng với các điểm cực), các điểm cực của hàm truyền chính là các trị riêng của ma trận A. Khi đó, ổn định BIBO là tương đương với ổn định tiệm cận.

Kết luận: Nếu hệ vừa điều khiển được vừa quan sát được, thì ổn định BIBO tương đương với ổn định tiệm cận.

Conclusion: If the system is both controllable and observable, then BIBO Stability ⇔ Asymptotical Stability

Một số phương pháp kiểm tra tính ổn định hệ tuyến tính

•

Ổn định tiệm cận [Asymptotically stable]

o Nếu tất cả các trị riêng của ma trận A nằm bên trái trục ảo [All the eigenvalues of A lie in the LHP]

•

Ôn định BIBO

o Tiêu chuẩn Routh-Hurwitz [Routh-Hurwitz criterion]

o Phương pháp quĩ tích nghiệm [Root locus method]

o Tiêu chuẩn Nyquist, [Nyquis criterion], …...

Nguyen Quang Hoang

12

Department of Applied Mechanics

6

02/12/2015

Ổn định tiệm cận hệ tuyến tính

Nhắc lại một số thuật ngữ

Tính được hàm truyền

Đối với hệ tuyến tính (SISO)

x = Ax + b ,

n

u

x Î Â , u Î Â

n( s)

-1

G( s) =

= (

c sI - )

A

b

y =

,

cx

y Î Â

d( s)

Điểm không: là nghiệm của n(s)=0

Điều kiện ổn định trong miền phức (s)

Điểm cực: là nghiệm của d(s)=0

Với một hệ được thể hiện bằng hàm

Đa thức đặc trưng: d(s)=0

truyền G(s),

Hệ là ổn định BIBO

Phương trình đặc trưng: d(s)=0

n( s)

( s - 1)( s + 1)

Tất cả các điểm cực của G(s) nằm ở nửa

G( s) =

=

2

d( s)

( s + 2)( s + 1)

trái của mặt phẳng phức

Điểm không của G(s):  1

Điểm cực của G(s): 2,  j

Hệ là ổn định tiệm cận

Nguyen Quang Hoang

13

Department of Applied Mechanics

Ổn định của hệ tuyến tính

Ý tưởng của điều kiện ổn định

Xét ví dụ:

y + y

a = u( t),

y(0) = y 0

1

 sY( s) - y(0) + Y

a ( s) = U( s)  Y( s) =

( U( s)+ y(0 ))

s + a

Xét ổn định tiệm cận

1

 Y( s) =

y(0)

u(t)=0 hay U(s)=0

s + a

1

 y( t)

-

=  [ Y( s)] = y(0) - t e a  0  Re( )

a > 0

1

Xét ổn định BIBO

 Y ( s) =

U( s) = G( s) U( s) y(0)=0

s + a

1

-

1

 y( t) =  [ Y( s)]

-

=  [ G( s) U( s)]

t

t

=

g( t) u( t - t) dt =

e at

-

u( t - t) dt

ò

ò

0

0

t

t

y( t) £

e at

-

⋅ u( t - ) d £

e a

- t

t

t

dt ⋅ u

ò

ò

max

0

0

Bị chặn nếu Re( a) > 0

Nguyen Quang Hoang

14

Department of Applied Mechanics

7

02/12/2015

Ổn định tiệm cận hệ tuyến tính

Xét hệ bậc hai với kích động xung va chạm

2

2

y + 2 zw y + w y = w (

d t),

y(0) = 0,

y(0) = 0

n

n

n

2

w

2

2

 s Y ( s) + 2 zw Y ( s) + w Y ( s) = U( s)  Y ( s) n

=

U( s)

n

n

2

2

s + 2 zw s + w

n

n

1 (

d t)

2

w

y( t)

z > 1

Cản lớn

n

2

2

s + 2 zw s + w

z = 1

Cản tới hạn

n

n

1 > z > 0 Cản nhỏ

z = 0

Không cản

1 > z > 0

z = 0

2

p = zw

-

 jw 1 - z

1,2

n

n

p =  jw

1,2

n

w

h( t)

n

zw

-

t

n

=

e

sin(

2

1 - z w t

h( t) = C sin( w t) + C cos( w t) 2

n )

1 - z

1

n

n

Nguyen Quang Hoang

15

Department of Applied Mechanics

Ổn định tiệm cận hệ tuyến tính

Xét hệ bậc hai với kích động xung va chạm

1 (

d t)

2

w

y( t)

n

2

2

s + 2 zw s + w

n

n

Cản nhỏ

Không cản

1 > z > 0

z = 0

2

p = zw

-

 jw 1 - z

1,2

n

n

p =  jw

1,2

n

w

h( t)

n

zw

-

t

n

=

e

sin(

2

1 - z w t

h( t) = C sin( w t) + C cos( w t) 2

n )

1 - z

1

n

n

lim h( t) = lim

- w t

lim h( t) = lim C sin( w t + a) = ???

n

C e

sin( w t + a) = 0

n

1

d

t ¥

t¥

t ¥

t¥

h( t) £ C

Ổn định tiệm cận

Ổn định biên

Nguyen Quang Hoang

16

Department of Applied Mechanics

8

02/12/2015

Ổn định tiệm cận hệ tuyến tính

Xét hệ bậc hai với kích động xung va chạm

1 (

d t)

2

w

y( t)

n

2

2

s + 2 zw s + w

n

n

Cản tới hạn

Cản lớn

z = 1

p = w

-

= l

x > 1

1,2

n

1,2

h( t)

- w t

p = w

z

-  z -

<

n (

2

1

0

1,2

)

n

= e

( C + tC

1

2 )

p t

p t

1

2

- w t

h( t) = C e

+ C e

n

1

2

lim h( t) = lim e

( C + tC = 0

1

2 )

t ¥

t ¥

lim h( t) = lim ( p t

p t

1

2

C e

+ C e

= 0

1

2

)

t ¥

t ¥

Ổn định tiệm cận

Ổn định tiệm cận

Nguyen Quang Hoang

17

Department of Applied Mechanics

Ổn định tiệm cận hệ tuyến tính

w

h( t)

n

zw

-

t

w = 3;

n

=

e

sin(

2

1 - z w t

n

2

n )

1 - z

z = [0, 0.3, 0.7, 0.98];

3

2

1

0

mm]

y [

3

-1

2

-2

1

-30

2

4

6

8

10

0

t [s]

3

-1

2

-2

1

-3

-2

-1

0

1

2

0

mm]

y [ -1

-2

-30

2

4

6

8

10

Xem: he_bac_2.m

t [s]

Nguyen Quang Hoang

18

Department of Applied Mechanics

9

02/12/2015

Ổn định của hệ tuyến tính

w

h( t)

n

zw

-

t

w = 3;

n

=

e

sin(

2

1 - z w t

n

2

n )

1 - z

z = [0, 0.3, 0.7, 0.98];

3

2

1

0

mm]

y [ -1

3

2

-2

1

-30

2

4

6

8

10

0

t [s]

3

-1

2

-2

1

-3

-2

-1

0

1

2

0

mm]

y [ -1

-2

-30

2

4

6

8

10

t [s]

Nguyen Quang Hoang

19

Department of Applied Mechanics

Ổn định của hệ tuyến tính

1 (

d t)

2

w

y( t)

n

2

p = zw

-

 jw 1 - z

2

2

s + 2 zw s + w

1,2

n

n

n

n

ổn định tiệm cận

ổn định tiệm cận

ổn định biên

Nguyen Quang Hoang

20

Department of Applied Mechanics

10

02/12/2015

Ổn định của hệ tuyến tính

không ổn định

không ổn định

không ổn định

Nguyen Quang Hoang

21

Department of Applied Mechanics

Ổn định của hệ tuyến tính

Một số lưu ý

 Đối với một hệ tổng quát (phi tuyến), điều kiện ổn định BIBO và điều kiện ổn định tiệm cận là khác nhau.

 Đối với hệ thống tuyến tính bất biến thời gian (LTI) (mà chúng ta có thể sử

dụng biến đổi Laplace và nhận được hàm truyền), các điều kiện ổn định như

nhau.

 Trong khóa học này, chúng ta chỉ quan tâm hệ thống LTI, khái niệm "ổn định"

sử dụng cho cả hai: ổn định BIBO và ổn định tiệm cận.

 Ổn định biên nếu (hệ dao động):

o

G(s) không có điểm cực nằm trong nửa trái mặt phẳng phức, và o

G(s) có một điểm cực đơn nằm trên trục ảo, và

o

G(s) không có điểm cực bội (nghiệm bội) nằm trên trục ảo.

1

1

G( s) =

G( s) =

2

s( s + 4)( s + 1)

2

2

s( s + 4) ( s + 1)

Ổn định biên

Không ổn định biên

 Unstable if a system is neither stable nor marginally stable.

Nguyen Quang Hoang

22

Department of Applied Mechanics

11

02/12/2015

Ổn định tiệm cận hệ tuyến tính

 Hệ có điểm cực kép (nghiệm kép)

é

ù

é 2

2

ù

-

-

s

w

-

2 s

w

1 ê

ú



= t sin t

w

1 ê

ú



= t cos t

w

ê 2

2 2

( s

w ) ú

+

ê 2

2 2 ú

ë

û

( s + w )

ë

û

é

ù

é

ù

-

...

-

....

1

2

ê

ú



= t sin t

w

1

2

ê

ú



= t cos t

w

ê 2

2 3

( s

w ) ú

+

ê 2

2 2 ú

ë

û

( s + w )

ë

û

Ổn định biên có dẫn đến ổn định BIBO ?

2 s

 Xét ví dụ

G( s) =

 Với đầu vào tác động

2

( s + 1)

1

u( t) = sin t  U( s) =

 Đầu ra của hệ

2

( s + 1)

2 s

1

Y ( s) = G( s) U( s) =

Thực chất đây là trường hợp hệ ổn

2

2

( s + 1) ( s + 1)

định biên (hệ dao động), nhưng do

é

ù

kích động điều hòa có tần số kích

-

2 s

1

y( t)

ê

ú



= 

= t sin t

ê

động trùng tần số dao động riêng.

2

2

( s

1) ú

+

ë

û

Xảy ra hiện tượng công hưởng.

Nguyen Quang Hoang

23

Department of Applied Mechanics

Ổn định của hệ tuyến tính

Tóm lại

jw

Với s

Vùng ổn

Vùng

i là điểm cực của hàm truyền G(s), hệ

định

không ổn

G(s) là:

định

 ổn định BIBO, ổn định tiệm cận nếu

Re(s

0

s

i) < 0, với mọi i.

 ổn định biên nếu

Vùng ổn

Vùng

o

Re(si)  0 với mọi i, và

định

không ổn

o

có một nghiệm đơn Re(s

định

i) = 0

 không ổn định nếu hệ không nằm trong

hai trường hợp trên (ổn định và ổn định

biên).

Nguyen Quang Hoang

24

Department of Applied Mechanics

12

02/12/2015

Ổn định của hệ tuyến tính

Xét ổn định của hệ bậc 2 (hệ cơ)

x

x

k

f(t)

X( s)

f(t)

X( s)

m=1

G( s) =

G( s) =

F( s)

m=1

F( s)

1

1

=

=

2

s

2

s + k

Các điểm cực =

Các điểm cực =

Ổn định?

Ổn định?

k

s

Ngoại trừ

F( s) =

, F( s) =

2

2

s + k

s + k

x

x

k

f(t)

X( s)

f(t)

X( s)

m=1

G( s) =

G( s) =

F( s)

m=1

F( s)

c

1

c

1

=

=

2

s + cs

2

s + cs + k

Các điểm cực =

Các điểm cực =

Ổn định?

Ổn định?

Nguyen Quang Hoang

25

Department of Applied Mechanics

Ổn định của hệ tuyến tính

Tóm lại

 Ổn định của hệ LTI

o ổn định (BIBO và tiệm cận), ổn định biên, và không ổn định.

o Ổn định của hệ có hàm truyền G(s) xác định bởi các điểm cực của G(s).

 Tiêu chuẩn ổn định Routh-Hurwitz xác định tính chất ổn định mà không cần tính các điểm cực của hệ.

Nguyen Quang Hoang

26

Department of Applied Mechanics

13

02/12/2015

Ổn định của hệ tuyến tính

Bài tập

G( s)

Các điểm không z

Ổn định/ổn định

và điểm cực p

biên/không ổn định

5( s + 2)

2

( s + 1)( s + s + 1)

5( s

- + 2)

2

( s + 1)( s + s + 1)

5

2

( s - 2)( s + 3)

2

( s + 3)

2

( s + 1)( s - s + 1)

5( s

- + 2)

2

2

( s + 1)( s + 1)

1

2

( s - 1)( s + 1)

Nguyen Quang Hoang

27

Department of Applied Mechanics

Phương pháp Lyapunov

Ổn định của hệ tuyến tính (LTI): Phương

pháp Lyapunov

Tiêu chuẩn ổn định Lyapunov xét tính ổn định của hệ trực tiếp trong không gian trạng thái, thích hợp đối với các hệ được mô tả bởi mô hình trạng thái.

phương pháp này áp dụng được cho cả hệ tuyến tính và phi tuyến.

x  Ax  Bu

x  f( ,

x u)

Nguyen Quang Hoang

28

Department of Applied Mechanics

14

02/12/2015

Phương pháp Lyapunov

Định lý: Hệ động lực

x  Ax  Bu

y  Cx  Du

ổn định BIBO khi và chỉ khi nó ổn định tiệm cận Lyapunov, tức là khi và chỉ khi các quỹ đạo trạng thái tự do (u =0) có hướng tiến về gốc tọa độ và kết thúc tại đó.

Chứng minh: Hệ động lực trên ổn định BIBO, tức là tất cả các trị riêng của ma trận hệ thống A có phần thực âm.

Quĩ đạo trạng thái tự do (u = 0), xác định bởi

Do tất cả các trị riêng của ma trận hệ thống A có phần thực âm, nên lim ( t)  lim t

eA

x

(

x 0)  0

t 

t 

Do đó, để kiểm tra tính ổn định tiệm cận Lyapunov (cũng là tính ổn định BIBO), ta chỉ cần kiểm tra xem quĩ đạo trạng thái của hệ trong quá trình tự do có hướng về

gốc tọa độ và kết thúc tại đó không.

Nguyen Quang Hoang

29

Department of Applied Mechanics

Phương pháp Lyapunov

Tư tưởng của phương pháp Lyapunov

Dựa trên phương diện năng lượng: Năng lượng còn tồn tại bên trong hệ do tác động tức thời bên ngoài đưa vào được mô bởi một hàm không âm. Hệ sẽ ổn định (tiệm cận) ở trạng thái cân bằng của nó nếu trong lân cận điểm cân bằng đó hàm năng lượng này của hệ luôn có xu hướng giảm dần về 0.

x(0)

x( t)

Giả sử bao quanh điểm cân bằng (x = 0) có họ

các đường cong khép kín v. Các đường cong

này được xem như các biên lân cận điểm gốc

x

0. Để xem quĩ đạo trạng thái x(t) tự do (u=0) 0

có tiến về gốc tọa độ hay không, ta chỉ cần

xem quĩ đạo x(t) có cắt tất cả các đường cong

thuộc họ v từ bên ngoài vào trong hay không

x( t)

và nếu điều đó xảy ra thì chắc chắn x(t) phải có hướng tiến về gốc tọa độ và kết thúc tại đó.

Nguyen Quang Hoang

30

Department of Applied Mechanics

15

02/12/2015

Phương pháp Lyapunov

Phương pháp Lyapunov gồm hai bước:

1. Xây dựng họ các đường cong v khép kín

x(0)

x( t)

 v

chứa điểm cân bằng (gốc tọa độ) bên trong.

2. Kiểm tra xem quĩ đạo trạng thái x(t) tự do (u=0) có cắt tất cả các đường cong thuộc họ



v từ bên ngoài vào trong hay không. Hiển

x

nhiên, để x(t) cắt đường cong thuộc họ v từ

0

ngoài vào trong là tại điểm cắt đó, tiếp tuyến

với quĩ đạo tại điểm cắt đó phải tạo với véc tơ

v vuông góc với đường cong theo hướng từ

x( t)

trong ra ngoài một góc lớn hơn 90.

Định lý Lyapunov (về ổn định tiệm cận): Nếu tồn tại hàm V(x) thỏa mãn điều kiện:

 Khả vi, xác định dương, tức V ( )

x  0

x

  0 & V ( )

x  0  x  0

dV ( )

x

V

 dx

 và





0 (đạo hàm V(x) dọc theo quĩ đạo trạng thái tự do âm).

dt

x dt

Thì hệ là ổn định tiệm cận Lyapunov tại gốc (x = 0) hay ổn định BIBO.

Hàm V(x) khi đó được gọi là hàm Lyapunov.

Nguyen Quang Hoang

31

Department of Applied Mechanics

Phương pháp Lyapunov

Chứng minh Định lý Lyapunov (về ổn định tiệm cận): Hàm khả vi xác định dương V(x) có tính chất: khi ta cắt nó bằng một mặt phẳng V = Ck song song với đáy (không gian trạng thái) va chiếu xuống đáy thì ta sẽ

được một đường cong kín chứa gốc tọa độ. Đường cong vk ứng với Ck nhỏ hơn thì nằm bên trong đường cong ứng với Ck lớn hơn.

C  C  v

nằm bên trong

v

1

2

1

2

V

x 2

V(x)= C 3

V(x)= C 2

x 1

V(x)= C 1

x



 v

0

x 2

0 < C1 < C2 < C3

x( t)

x 1

x(0)

Nguyen Quang Hoang

32

Department of Applied Mechanics

16

02/12/2015

Phương pháp Lyapunov

Véc tơ vuông góc với tiếp tuyến của đường cong vk và hướng theo chiều tăng Ck là

T

T

 V

 

 V



V



V

 

  grad V  

  

,

,,



v

x







 x



x



x







1

2

n 

Sẽ hướng từ trong ra ngoài đường cong vk.

V

x 2

V(x)= C 3

V(x)= C 2

x 1

V(x)= C 1

x



0

x 2

 v

x 1

0 < C1 < C2 < C3

x( t)

Minh họa hàm xác định dương đối với hai biến x

x(0)

1 và x2

Nguyen Quang Hoang

33

Department of Applied Mechanics

Phương pháp Lyapunov

Do có

T

dV

V

 dx

 V



V



V

  dx



 

,

,,



dt

x

 dt

 x



x



x





 dt

1

2

n 

T dx

dx

 grad V 

 grad V

cos



 dt

dt

0 < C1 < C2 < C3

mà x  dx / dt chính là tiếp tuyến V(x)= C 3

của quĩ đạo trạng thái x(t), nên với

điều kiện dV/dt < 0, góc  tạo bởi

gradV và dx/dt phải là góc tù (lớn

V(x)= C 2

hớn 90), tức là quĩ đạo trạng thái

x(t) sẽ cắt các đường cong vk theo

V(x)= C 1

x

hướng từ ngoài vào trong.



Lưu ý:

•

Định lý này áp dụng cho cả hệ phi tuyến.

 v

x(0)

•

Nếu không tìm được hàm Lyapunov thì không thể

x( t)

khẳng định về tính ổn định của hệ. Định lý cho ta điều kiện đủ để xác định tính ổn định của hệ.

Nguyen Quang Hoang

34

Department of Applied Mechanics

17

02/12/2015

Phương pháp Lyapunov

Ví dụ minh họa. Xét hệ động lực mô tả bởi

 4



2

0   x 

1 0

1

dx

 u 



   

 1

 Ax  Bu 

2



5



1

x

 1 2  



  2  

 u

dt

2

 0

1



2

  x

  0 1  



3



   



Chọn hàm khả vi, xác định dương

2

2

2

V ( )

x  x  x  x  0

x

  0

1

2

3

Đạo hàm V(x) dọc quĩ đạo x(t) của quá trình tự do (u = 0), ta được dV

¶ V

=

x u=0

dt

¶x

=[ 2 x

2 x

2 x ]Ax

1

2

3

2

2

2

= - 8 x - 10 x - 4 x < 0

" x ¹ 0

1

2

3

Theo định lý Lyapunov về ổn định tiệm cận, ta kết luận hệ ổn định.

Nguyen Quang Hoang

35

Department of Applied Mechanics

Ổn định của hệ tuyến tính (LTI) : PP Lyapunov Khi áp dụng phương pháp trực tiếp Lyapunov đối với hệ tuyến tính x = Ax Người ta hay sử dụng hàm trơn xác định dương dạng det A ¹ 0

( )

T

x

x P ,

T

V

=

x

P = P > 0

P là ma trận đx xác định dương

Đạo hàm V(x) theo thời gian dọc quỹ đạo ta có V( )

T

T

x = x Px + x Px = (A ) T

T

x Px + x Px

T

T

T

T

= x A Px + x PAx = x ( T

A P + PA)

T

x = -x Qx

Như thế, nếu tồn tại ma trận Q đối xứng x.đ.d. thỏa mãn ( T

A P + PA) = -Q

thì V là hàm Lyapunov và điểm cân bằng x = 0 là ổn định toàn cục. Phương trình trên được gọi là phương trình ma trận Lyapunov.

Nguyen Quang Hoang

36

Department of Applied Mechanics

18

02/12/2015

Ổn định của hệ tuyến tính

Xét hệ động lực tuyến tính mô tả bởi phương trình x = A ,

x

(

x 0) = x

det A ¹ 0

0

Hệ trên là ổn định tiệm cận khi và chỉ khi đối với mỗi ma trận Q thực đối xứng xác định dương nghiệm P của phương trình ma trận Lyapunov sau T

A P + PA = -Q

là (thực đối xứng) xác định dương.

Chứng minh (điều kiện cần)

1. Nếu hệ ổn định tiệm cận thì với mỗi Q thực đối xứng x.đ.d tồn tại P thực đối xứng xác định dương thỏa mãn pt ma trận Lyapunov.

Sử dụng phương pháp phản chứng: giả sử rằng với mọi x  0, ta có T

x Px

0,

&

T

>

-x Qx < 0

" x ¹ 0

Mà hệ không ổn định tiệm cận hay ma trận hệ A không ổn định tiệm cận. Điều này có nghĩa là với điều kiều đầu x(0)  0, nghiệm ( t t )

0

( t)

e -

= A

x

(

x 0) giới nội và không hội tụ về gốc tọa độ.

Nguyen Quang Hoang

37

Department of Applied Mechanics

Ổn định của hệ tuyến tính

Chứng minh (điều kiện đủ)

2. Nếu tồn tại Q thực đối xứng xác định dương & P x.đ.d. thỏa mãn phương trình ma trận Lyapunov thì hệ là ổn định tiệm cận. Thực vậy, ta chọn hàm Lyapunov và tính đạo hàm của nó

V ( )

T

x = x P ,

x

.

p d.

Kết luận: hệ

V( )

T

T

T

T

T

x = x Px + x Px = x PAx + x A Px ổn định tiệm

cận

T

= x (

T

PA + A

)

T

P x = -x Qx

n. d.

Định lý về sự tồn tại P thỏa mãn phương trình ma trận Lyapunov. Nếu ma trận A là Hurwitz (ma trận bền) thì với mọi Q đối xứng x.đ.d sẽ tồn tại P đối xứng x.đ.d thỏa mãn phương trình ma trận Lyapunov.

Định nghĩa: Ma trận A được gọi là Hurwitz nếu các trị riêng của nó có phần thực âm, Re(i) < 0, i = 1,2,…,n.

Nguyen Quang Hoang

38

Department of Applied Mechanics

19

02/12/2015

Ổn định của hệ tuyến tính

Chứng minh (We make a constructive proof)

Với ma trận Q thực đối xứng x.đ.d, ta xét nghiệm P

¥

T

T

T

T

t

t

=

e

e dt

ò A

A

P

Q

Từ Q = Q > 0  P = P > 0

0

Chú ý rằng, tích phân trên chỉ hội tụ khi ma trận A là Hurwitz (các trị riêng của A có phần thực âm).

Ta sẽ chỉ ra, ma trận P thỏa mãn phương trình ma trận Lyapunov (

T

PA + A

)

P = -Q

Ta có

¥

T

T

T

T

T

é T A t

A t

A t

A t

PA

A P

A e Q e

e

e

ù

+

=

+

Q

A dt

ò ê

ú

ë

û

0

¥ d

¥

T

T

T

T

é A t

A t ù

A t

A t

=

e

e

Q

dt = e

e

Q

= -

ò

Q

dt ê

ú

ë

û

0

0

Nguyen Quang Hoang

39

Department of Applied Mechanics

Ổn định của hệ tuyến tính (LTI): PP Lyapunov

Kết luận về tính ổn định tiệm cận từ trị riêng của ma trận A Xét điểm cân bằng x = 0 của hệ LTI với ma trận A chéo hóa được bằng ma trận Q

x = Ax

x = Q ,

y x = Qy

-1

1

x = Qy = AQy  y = Q AQ , y

(

y 0)

-

= Q

(

x 0)

Với A chéo hóa được, tức tồn tại ma trận Q : 1

-

L = Q AQ = diag( l , l ,..., l ) 1

2

n

y = l y ,

y (0) = y , i = 1,2,..., n

i

i i

i

i,0

y ( t)

l ti

= e y (0), i = 1,2,..., n

i

i

Như thế, hệ tuyến tính với A chéo

l t

l t

l t

hóa được là ổn định tiệm cận nếu

1

2

(

y t) = diag( e , e ,..., n e ) (

y 0)

Re( l ) < 0, " i = 1,2,..., n L t

1

(

x t)

Q (

y t)

e

-

=

= Q

Q

(

x 0)

i

Nếu có Re(i) = 0, thì hệ không ổn

L t

l t

l t

l t

1

2

e

= diag( e , e ,..., n e ,)

định tiệm cận, nhưng ổn định biên.

Nguyen Quang Hoang

40

Department of Applied Mechanics

20

Table of Contents

ức