Thiết kế bộ điều khiển thích nghi theo mô hình mẫu cho hệ truyền động qua bánh răng

Lê Th ị Thu Hà và Đtg Tạp chí KHOA H ỌC & CÔNG NGH Ệ 112(12)/2: 63 - 68 THI ẾT K Ế B Ộ ĐIỀU KHI ỂN THÍCH NGHI THEO MÔ HÌNH M ẪU CHO H Ệ TRUY ỀN ĐỘ NG QUA BÁNH R ĂNG Lê Th ị Thu Hà*, Tr ần Th ị Thanh Th ảo Tr ường Đại h ọc K ỹ thu ật Công nghi ệp – ĐH Thái Nguyên TÓM T ẮT Bài báo trình bày ph ươ ng pháp thi ết k ế b ộ điều khi ển thích nghi theo mô hình m ẫu cho h ệ truy ền động qua bánh r ăng trên c ơ s ở s ử d ụng mô hình tr ạng thái c ủa h ệ. Kh ả n ăng bám ti ệm c ận t ốt theo mô hình m ẫu của h ệ có ch ứa đầ y đủ các thành ph ần b ất đị nh sinh ra t ừ hi ệu ứng khe h ở, ma sát, độ dẻo bánh r ăng đã được ch ứng minh c ả v ề lý thuy ết và mô ph ỏng. Từ khóa: điều khi ển thích nghi, mô hình m ẫu, h ệ th ống bánh r ăng, khe h ở, mômen ma sát. * ĐẶT V ẤN ĐỀ M, M − f1 f 2 là các moment ma sát c ủa hai Điều khi ển bám ổn đị nh h ệ truy ền độ ng qua bánh r ăng 1 và 2. bánh r ăng mang đầy đủ trong nó các y ếu t ố bất đị nh nh ư khe h ở, độ không c ứng v ững c ủa − c là ch ỉ s ố đo độ c ứng c ủa v ật li ệu làm vật li ệu làm bánh r ăng luôn gi ữ vai trò trung bánh r ăng. Nó chính là đại l ượng đánh giá độ tâm trong l ớp các bài toán điều khi ển h ệ cứng v ững c ủa h ệ truy ền độ ng. truy ền độ ng. α − L là góc kh ớp hai r ăng. Đây là ch ỉ s ố đo độ khe h ở gi ữa các bánh r ăng. V ới hai r ăng ăn M α = ° d kh ớp chính xác tuy ệt đố i thì L 20 . Các 1 α > ° cặp r ăng có khe h ở luôn có L 20 Mms 1 Md = + − J1 Jd JJ 1, 2 là các moment quán tính của c ặp bánh r ăng 1,2 và c ủa độ ng c ơ d ẫn động. Mc P = −1 M − i12, i 21 i 12 ms 2 là t ỷ s ố truy ền c ủa hai bánh răng. 2 Mc − Mc là moment c ản (t ải), được xem nh ư Hình 1: Minh h ọa h ệ truy ền động qua bánh r ăng nhi ễu tác độ ng vào h ệ. Theo [3] thì h ệ truy ền độ ng qua bánh r ăng, có ϕ, ϕ ɺ sơ đồ c ấu trúc minh h ọa ở hình 1, không có − 2 2 là v ị trí và t ốc độ c ủa bánh r ăng th ụ ϕ kho ảng ch ết gi ữa các bánh r ăng, s ẽ mô t ả động và 2 sẽ được xem là tín hi ệu ra c ủa hệ. được b ởi mô hình Euler-Lagrange: Nếu gi ữa hai bánh r ăng có các khe h ở thì khi  ϕ+2 2 αϕ +=− ϕ J111ɺɺ crcos ( 1122 i ) Md M f 1 ở ch ế độ khe h ở, moment d ẫn độ ng ở đầ u vào  không có tác d ụng thay đổ i t ốc độ c ủa bánh Jϕɺɺ − cr2cos 2 αϕ ( + i ϕ ) =−− M M  222 2211c f 2 (5) răng b ị độ ng và bánh r ăng b ị độ ng lúc đó ch ỉ trong đó còn ch ạy theo quán tính. Nói cách khác, ở giai đoạn khe h ở, h ệ s ẽ có mô hình [3]: r, r − 1 2 là bán kính vòng ngoài c ủa hai bánh  ϕ = − J1ɺɺ 1 Md M f 1 răng.  Jϕɺɺ = − M + M  2 2c f 2 (6) * Tel: 0977008928; Email: hahien1977@gmail.com 63 Lê Th ị Thu Hà và Đtg Tạp chí KHOA H ỌC & CÔNG NGH Ệ 112(12)/2: 63 - 68 THI ẾT K Ế B Ộ ĐIỀU KHI ỂN ϕ ϕ tuy ến tính v ới v ận t ốc và d11( ,), td 22 ( ,) t Xây d ựng mô hình cho h ệ truy ền độ ng qua là nh ững thành ph ần moment ma sát ph ụ bánh r ăng ở c ả hai ch ế độ làm vi ệc thu ộc gia t ốc, moment t ải. V ới nh ững ký hi ệu Ứng v ới t ừng lo ại mô hình (5) và (6) mô t ả cho trong (8) này, mô hình Euler-Lagrange hai ch ế độ làm vi ệc khác nhau c ủa h ệ, ng ười (7) được vi ết l ại thành: ta th ường áp dung các ph ươ ng pháp điều / Jϕθϕɺɺ+( + iM ϕ ) =−− θϕ ɺ d khi ển khác nhau. Th ường dùng nh ất là s ử 11 11122d 11 1  − − Jϕθϕɺɺ−1( + i 1 ϕ ) =− θϕ/ ɺ − d dụng các công c ụ nh ận d ạng ho ặc x ấp x ỉ khe  22 2 2 121 22 2 (9) hở để t ừ đó s ử d ụng nguyên lý điều khi ển bù nh ằm giúp loại b ỏ được mô hình (6) trong quá Để chuy ển (9) về d ạng mô hình tr ạng thái, tr ước trình thi ết k ế b ộ điều khi ển. tiên, t ừ ph ươ ng trình th ứ hai trong (9) ta có: ϕ= θϕθϕɺɺ ++−/ ɺ ϕ  Tuy nhiên, n ếu xem khe h ở c ũng là m ột thành 1i 12222( J 22 d 2) 2  ph ần b ất đị nh trong h ệ, gi ống nh ư các tham s ố =θϕɺɺ + θϕ ɺ −i ϕ + d c đo độ c ứng c ủa v ật li ệu làm bánh r ăng, 32 42 122 3 (10) M, M với: f1 f 2 mô t ả các thành ph ần ma sát hay α did=θ, θθ = J , θ = i θθ / góc kh ớp hai r ăng L , thì ta có th ể ghép hai 3 1222 3 22 4 1222 mô hình (5) và (6) chung l ại v ới nhau thành là các thành ph ần b ất đị nh h ằng s ố và hàm s ố một mô hình t ổng quát: tươ ng ứng. ϕ  ϕ+⌢ 2 2 αϕ +=− ϕ 1 J11ɺɺ crLL 1cos ( 1122 i ) M df M 1 Đạo hàm theo th ời gian hai v ế c ủa cho  ϕ−⌢ 2 2 αϕ + ϕ =−− trong công th ức (10), ta có: J22ɺɺ crLL 2cos ( 2211 i ) M cf M 2 (7) ϕɺ= θϕ ɺɺɺ + θϕ ɺɺ −i ϕ ɺ + d ⌢ 1 32 42 122 4 (11) c trong đó tham s ố được đị nh ngh ĩa là: = ɺ trong đó d4 d 3 là thành ph ần hàm b ất đị nh, ⌢ c ở ch ế độ ăn kh ớp c =  được gi ả thi ết c ũng b ị ch ặn. T ừ đây ta suy ra 0 ở ch ế độ khe h ở ϕɺɺ= θϕ(4) + θϕ ɺɺɺ −i ϕ ɺɺ + d Nh ư v ậy mô hình (7) này s ẽ ch ứa trong nó t ất 1 32 42122 5 (12) cả các y ếu t ố b ất đị nh c ủa h ệ. Đây là nh ững d= d ɺ tham s ố ho ặc các hàm r ất khó, ho ặc không th ể với 5 4 . xác định được m ột cách đủ chính xác. Có th ể Thay (11), (12) vào ph ươ ng trình th ứ nh ất c ủa kể đế n đó là độ không c ứng v ững c của v ật mô hình (9), ta được: α li ệu, góc ăn kh ớp L gi ữa hai bánh r ăng, θϕ(4) + θϕ − ϕ ++ J13( 2 42ɺɺɺ i 122 ɺɺ d 5 ) M, M moment ma sát f1 f 2 trên các tr ục +θθϕθϕ +−++ ϕ ϕ = 1() 32ɺɺ 42 ɺ i 122 d 3 i 122 M truy ền độ ng, moment t ải c , khe h ở. =−Mθθϕ/ ()ɺɺɺ + θϕ ɺɺ − idd ϕ ɺ +− d 132 42 122 4 1 Ti ếp theo, để đơn gi ản hóa trong trình bày, ta θ và điều này d ẫn đế n: sẽ s ử d ụng các ký hi ệu k cho h ằng s ố và dk =θϕ(4) ++( θ θθϕ/ )ɺɺɺ + cho hàm s ố b ất đị nh nh ư sau: Md J132 J 14 132 θ= ⌢ 2 2 α +θθ/ +− θθ ϕ + 1cr L 1 co s L ( 14 13J 112 i ) ɺɺ 2 θ− 1= ⌢ 2 2 α 2cr L 2 co s L +θθ − θ/ ϕ + ( 1 4 1i 12) ɺ 2 =θ/ ϕ + ϕ Mf 1 11ɺ d 11( , t ) +() +θ + θ / + − =−θ/ ϕɺ − ϕ Jd15 13 d 14 dd 1 MMc f 2 22 dt 22( , ) (8) Sử d ụng ký hi ệu vector c ủa tham s ố h ằng b ất θ/ θ / 1, 2 θθθ θ với là hai h ằng s ố b ất đị nh đo thành f, g ph ần moment ma sát độ ng được gi ả thi ết là định với: 64 Lê Th ị Thu Hà và Đtg Tạp chí KHOA H ỌC & CÔNG NGH Ệ 112(12)/2: 63 - 68 =− − − θ = 1 u v ax01 ax 12 ax 23 (19) g J θ 1 3 (13) với ba h ằng s ố a0, a 1 , a 2 tùy ch ọn cho h ệ (18), θ θ− θ / 14 112i  1 ta s ẽ thu được h ệ kín d ạng tuy ến tính chu ẩn θθθ =−θθ/ + θθ −  f 14 13J 112 i điều khi ển v ới mô hình tr ạng thái: J θ   1 3 J θ+ θ/ θ  14 13 (14) 0 1 0  xɺ =0 0 1  x +b()θθθT x + v d(x , t )   f và hàm s ố b ất đị nh : − − − a0 a 1 a 2  1 d=−( Jd +θ d + θ / dd + ) A J θ 15 13 14 1 m (20) 1 3 (15) trong đó cũng nh ư t ừ th ực t ế là ta ch ỉ c ần quan tâm t ới det(sI− A ) =++ aasass2 + 3 ϕɺ 3m 012 tốc độ 2 , t ức là ch ỉ c ần quan tâm t ới ba bi ến Điều này d ẫn ta t ới ý t ưởng r ằng có th ể ch ọn tr ạng thái: các h ằng s ố a0, a 1 , a 2 theo ph ươ ng pháp gán x  ϕɺ  1 2 các điểm c ực s1, s 2 , s 3 tươ ng ứng v ới ch ất x =x  = ϕɺɺ  2   2  lượng mong mu ốn đặ t tr ước. Ch ẳng h ạn để h ệ x  ϕɺɺɺ  ổn đị nh, không có dao độ ng trong quá trình 3 2 (16) s1, s 2 , s 3 thì v ới ký hi ệu c ủa tín hi ệu đầ u vào: quá độ, ta ch ọn ba h ằng s ố th ực âm , rồi tính: = u M d − − −= (ss1 )( ss 2 )( ss 3 ) ta có d ạng mô hình tr ạng thái t ươ ng đươ ng =−++s3( s s ss )( 2 + ss + ss + ssssss ) − của mô hình Euler-Lagrange (7): 1 2 3 12 23 13 123 =− = ++ a0 sss 123 , a 1( ss 12 ss 23 ss 13 )  = xɺ1 x 2 a=−( sss + + )  ⇔ 2 123 (21)  = xɺ2 x 3  Nói cách khác, h ệ thu được (20) nh ờ b ộ điều =θθθT + + θ xɺ3 fx dt( x , ) g u khi ển vòng trong (19), tr ừ thành ph ần b ất định (17) θθθ f Thi ết k ế b ộ điều khi ển thích nghi theo mô cho b ởi (14), đã có đầy đủ t ất c ả các ch ất hình m ẫu lượng mong mu ốn đặ t tr ước. B ởi v ậy nhi ệm v ụ điều khi ển ti ếp theo bây gi ờ ch ỉ còn là lo ại b ỏ Gi ả thi ết r ằng h ệ truy ền độ ng có ma sát không θθθ ph ụ thu ộc gia t ốc, t ức là có gia t ốc r ất nh ỏ. sự ảnh h ưởng c ủa f trong h ệ kín (20). Khi đó ta có th ể b ỏ qua thành ph ần d(x , t ) Để làm được điều này ta s ẽ áp d ụng nguyên trong (17). Ngoài ra, n ếu nh ư ta có th ể xấp x ỉ tắc thích nghi theo mô hình m ẫu, t ức là ta s ẽ θ thi ết k ế thêm b ộ điều khi ển vòng ngoài để h ệ g được là h ằng s ố xác đị nh thì không m ất (20) bám ti ệm c ận theo được mô hình m ẫu, θ = 1 tính t ổng quát ta có th ể cho r ằng g . Khi suy ra t ừ (20) sau khi lo ại b ỏ đi s ự ảnh h ưởng θθθ đó (17) tr ở thành: của thành ph ần b ất đị nh f nh ư sau:  = + xɺ = x xɺmA m x m b w  1 2 (22)  = xɺ2 x 3 Hình 2 minh h ọa nguyên t ắc điều khi ển thích  T nghi theo mô hình m ẫu cho đố i t ượng (5) ở xɺ =θθθ x + u  3 f (18) ch ế độ ch ạy g ần đề u ( để b ỏ qua được các ma Dễ th ấy được r ằng khi s ử d ụng b ộ điều khi ển sát ph ụ thu ộc gia t ốc), g ồm hai vòng điều vòng trong: khi ển trong và ngoài. 65 Lê Th ị Thu Hà và Đtg Tạp chí KHOA H ỌC & CÔNG NGH Ệ 112(12)/2: 63 - 68 w v Lyapunov (26) với m ọi ma tr ận đố i x ứng xác Điều khi ển vòng u Đối t ượng x ∈ 3× 3 trong (19) (18) định d ươ ng Q R luôn có nghi ệm 3× 3 P ∈R cũng đố i x ứng xác đị nh d ươ ng. Từ công th ức đạ o hàm (25) của hàm xác định z Điều khi ển vòng e ngoài (27) dươ ng (24) thì theo lý thuy ết Lyapunov II, với b ộ ch ỉnh đị nh thích nghi tham s ố p(t ) : x Mô hình m ẫu m  −1 T (22) pɺ = E x b P e  z = xT p Hình 2. Sơ đồ điều khi ển thích nghi theo mô hình  (27) mẫu cho h ệ truy ền độ ng bánh r ăng sẽ có Để thi ết k ế b ộ điều khi ển vòng ngoài v ới ɺ =−T < ∀≠ nhi ệm v ụ là cho h ệ (20) bám ti ệm c ận theo V eQ e0, e 0 (28) được mô hình m ẫu (22), tr ước tiên ta c ần đế n Đó là điều ki ện đủ để được ph ươ ng trình mô t ả sai l ệch mô hình. lime (t ) = 0 < ∞ = − t→∞ e(t ) Ký hi ệu e x x m là sai l ệch mô hình. Khi và đó v ới các phép gán: tức là s ẽ có được tính bám ti ệm c ận c ủa (20) T theo mô hình m ẫu (22). Tuy nhiên, do v ới z = x p v= w − z và công th ức (28) thì Vɺ ch ỉ xác đị nh âm theo sai e trong đó p(t ) là vector tham s ố b ộ điều khi ển lệch , nói cách khác nó ch ỉ bán xác đị nh âm θθθ − vòng ngoài c ần ph ải xác đị nh, ta có theo e và p , nên c ũng ch ỉ đảm b ảo có e → 0 eɺ =A ex + b T (θθθ − p ) được tính ti ệm c ận c ủa , ch ứ ch ưa m (23) →θθθ kh ẳng đị nh được c ũng s ẽ có p , nên c ơ Sử d ụng hàm tr ơn xác định d ươ ng cấu ch ỉnh đị nh (27) không thay th ế được c ơ V (eee ) =T P +−()()θ pT E θ − p cấu nh ận d ạng tham s ố b ất đị nh trong mô hình. (24) Tổng k ết l ại, b ộ điều khi ển thích nghi theo ∈ 3× 3 với P, E R là đối x ứng xác đị nh d ươ ng mô hình m ẫu cho h ệ truy ền độ ng qua bánh tùy ch ọn, ta s ẽ có v ới (23) răng làm vi ệc ở ch ế độ có moment ma sát không ph ụ thu ộc gia t ốc, xây d ựng trên n ền ɺ =T + T −−(θθθ )T V eeeeɺP P ɺ2 p E p ɺ mô hình tr ạng thái (18) của h ệ, s ẽ được t ổng =T T + − e()Am P PA m e hợp qua các b ước nh ư sau: T T s, s , s − 2 ()θθθ −pE pxɺ − () Pb e  1. Ch ọn các điểm c ực 1 2 3 nằm bên trái   tr ục ảo, ứng v ới ch ất l ượng ổn đị nh mong T T =−eeT Q −−2()θθθ p E pxɺ − () Pb e  mu ốn c ủa h ệ kín r ồi tính các tham s ố cho b ộ   (25) điều khi ển vòng trong a0, a 1 , a 2 theo công trong đó th ức (21). Để h ệ kín không nh ững ổn đị nh mà T + = − = ϕɺ APm PA m Q (26) ở ch ế độ còn có tín hi ệu đầ u ra x1 2 bám w( t ) Rõ ràng, v ới vi ệc ch ọn các tham s ố a0, a 1 , a 2 ti ệm c ận theo tín hi ệu m ẫu , ta c ần ch ọn = − của b ộ điều khi ển vòng trong theo (21) có các chúng th ỏa mãn thêm s1 s 2 s 3 1 điểm cực s1, s 2 , s 3 ch ọn tr ước n ằm bên trái ∈ 3× 3 2. Ch ọn Q R đối x ứng xác đị nh d ươ ng A tr ục ảo thì ma tr ận m là ma tr ận b ền. Điều 3× 3 và tìm nghi ệm P ∈R cũng đố i x ứng xác này đảm b ảo ch ắc ch ắn r ằng ph ươ ng trình định d ươ ng c ủa ph ươ ng trình Lyapunov (26). 66 Lê Th ị Thu Hà và Đtg Tạp chí KHOA H ỌC & CÔNG NGH Ệ 112(12)/2: 63 - 68 Tli¶egom ki¶¶§¶ 3× 3 0.15 Q ∈R Q ∞ Ma tr ận được ch ọn có càng 0.1 lớn, t ốc độ bám c ủa (20) theo mô hình m ẫu (22) càng cao. Chú ý khi đó ph ải tr ả giá là độ 0.05 quá điều ch ỉnh càng l ớn. 0 3× 3 3. Ch ọn E ∈R đối x ứng xác đị nh d ươ ng. -0.05 E -0.1 Nếu ch ọn E có ∞ càng nh ỏ, t ốc độ ch ỉnh -0.15 định p(t ) càng cao, do đó quá trình quá độ -0.2 của h ệ càng ng ắn. 0 10 20 30 40 50 60 70 4. Xây d ựng b ộ điều khi ển vòng trong theo Hình 4. Sai l ệch bám (19), b ộ điều khi ển vòng ngoài theo (27) và 0.4 0.2 mô hình m ẫu theo (22) 0 KẾT QU Ả MÔ PH ỎNG -0.2 -0.4 Xét h ệ truy ền độ ng có mô hình (18). Ch ọn -0.6 = = =− s1 s 1 s 3 1 ta s ẽ có v ới (21): -0.8 -1 = = = -1.2 a01 , a 1 3 , a 2 3 = = -1.4 Q E10 I 3 -1.6 Ch ọn các ma tr ận ta có đồ th ị 0 10 20 30 40 50 60 70 80 90 100 = ϕɺ qu ỹ đạ o x1 2 của h ệ và tín hi ệu m ẫu w( t ) Hình 5. Tham s ố b ộ điều khi ển vòng ngoài cho ở hình 3. Nó cho ta th ấy tr ực quan được KẾT LU ẬN kh ả n ăng bám t ốt c ủa tín hi ệu đầ u ra c ủa h ệ Bộ điều khi ển c ủa bài báo được thi ết k ế trên theo tín hi ệu m ẫu. nền thích nghi gi ả đị nh rõ. e Hình 4 mô t ả sai l ệch . Nó xác nh ận tính Bằng k ết qu ả mô ph ỏng, bài báo còn ch ỉ ra t ừ bám ti ệm c ận theo mô hình m ẫu c ủa h ệ kín. hình 5 r ằng b ộ điều khi ển gi ới thi ệu ở đây còn đảm b ảo ch ất l ượng bám ngay c ả khi các Ngoài ra, các hình 5 còn cho th ấy m ặc dù các θθθ f tham s ố p(t ) của b ộ điều khi ển vòng ngoài tham s ố b ất đị nh của h ệ truy ền độ ng không nh ất thi ết ph ải bám theo giá tr ị b ất đị nh không ph ải là h ằng s ố, m ặc dù ở ph ần ch ứng θθθ minh ta ph ải gi ả thi ết nó ch ỉ là tham s ố h ằng f , song h ệ v ẫn có được ch ất l ượng bám ổn bất đị nh để có được s ự bi ến đổ i t ừ công th ức định r ất t ốt. Đặ c bi ệt n ữa là ở mô ph ỏng này ta (24) thành (25). θθθ (t ) còn có f là hàm thay đổi theo th ời gian. Theo lý thuy ết, vi ệc v ẫn có được tính bám ti ệm c ận t ốt được ngay c ả khi có hàm bất đị nh 25 θθθ f (t ) 20 có th ể không ph ải là s ự ng ẫu nhiên mà 15 vẫn đúng cho m ọi tr ường h ợp, không ch ỉ 10 riêng ở ph ần mô ph ỏng này. Suy ngh ĩ đó là 5 0 hợp lý vì th ực ch ất ở đây, để đưa ra được -5 công th ức (27) cho c ơ c ấu ch ỉnh đị nh, ta đã s ử -10 dụng lý thuy ết Lyapunov, v ốn ch ỉ là m ột điều -15 ki ện đủ . -20 -25 Bởi v ậy bài toán ch ứng minh ch ặt ch ẽ tính 0 10 20 30 40 50 60 70 80 θθθ (t ) bám ti ệm c ận c ủa h ệ v ẫn th ỏa mãn khi f Hình 3. Qu ỹ đạ o tín hi ệu ra so sánh là hàm bất đị nh, s ẽ là bài toán m ở ti ếp theo với tín hi ệu đặ t của nhóm tác gi ả. 67 Lê Th ị Thu Hà và Đtg Tạp chí KHOA H ỌC & CÔNG NGH Ệ 112(12)/2: 63 - 68 TÀI LI ỆU THAM KH ẢO 3. Ha,L.T.T. (2012), Modelling of transmission two-weel gearing System . Reaserch report, TNUT. 1. Couwder, R. (2006), Electric Drivers and 4. Menon, K. and Krishnamurty (1999), Control Electromechanical Systems . Elservier GB. of low friction and gear backlash in machine tool 2. Eutebach, T. and Pacas, J.M. (1999), feed drive systems. . Mechatronics 9, pp.33-52. Damping of torsional vibration in high dynamic 5. Thosen,S abd Fuchs,F.W. (2009), Speed drivers . 8. European Conference on Power control of torsional driver systems with backlash . European Conference on Power Electronics and Electronics and Applications EPE 99. Applications EPE 09. SUMMARY DESIGN OF MODEL REFERENCE ADAPTIVE CONTROLLER FOR GEARING TRANSMISSION SYSTEM Le Thi Thu Ha*, Tran Thi Thanh Thao College of Technology - TNU This paper presents the design method of the model reference adaptive controller for gearing transmission systems based on using this model. The asymptotic tracking behavior of the system in the presence of all uncertainties caused by effect of backlash, friction or cogwheel elasticity is proved theoretically and experimentally. Key words: Adaptive tracking, model reference, gearing systems, backlash, torsional moment. Ph ản bi ện khoa h ọc: PGS.TS. Lại Kh ắc Lãi – Đại h ọc Thái Nguyên * Tel: 0977008928; Email: hahien1977@gmail.com 68