Xem mẫu

  1. TỔNG HỢP ĐIỀU KHIỂN THÍCH NGHI HỆ THỐNG CHỐNG BÓ CỨNG BÁNH XE Ô TÔ KHI PHANH TRÊN CƠ SỞ MÔ HÌNH MỜ PGS.TS. LÊ HÙNG LÂN ThS. NGUYỄN VĂN TIỀM Bộ môn Điều khiển học, Khoa Điện - Điện tử Trường Đại học Giao thông Vận tải TS. LÊ CHUNG Khoa Kỹ thuật điều khiển Học viện Kỹ thuật Quân sự Tóm tắt: Hệ thống chống bó phanh (ABS – Anti-lock braking system) có vai trò rất quan trọng trong việc đảm bảo chất lượng khi phanh và tính dẫn hướng của ôtô. Đa số các bộ điều khiển ABS có bán ở trên thị trường đều dựa trên nguyên lý điều khiển on-off. Trên các xe ôtô hiện đại đều được trang bị ở mỗi bánh xe một bộ điều khiển ABS, mục đích là để điều khiển độ trượt tương đối giữa bánh xe và mặt đường khi phanh. Bài báo này đưa ra phương pháp tổng hợp hệ thống điều khiển thích nghi độ trượt này trên cơ sở logic mờ. Đánh giá hiệu quả của phương pháp thông qua các kết quả mô phỏng máy tính. Summary: The anti-lock braking system (ABS) is an important component of a complex steering system for the modern automobiles. Most of ABS controllers available on the market CT 2 are based on on-off controlling principle. All automobiles of latest type are fitted with an ABS controller, which aims to maintain a specified tire slip for each wheel during braking. This paper proposes a model of adaptive controller, based on fuzzy logic control to regulate the tire-slip. Simulation and test results are presented to form assessment of the method. I. ĐẶT VẤN ĐỀ 1.1. Chức năng, nhiệm vụ của bộ chống bó cứng bánh xe ô tô – ABS Fz Bộ ABS có nhiệm vụ là đảm bảo hiệu quả phanh tối ưu (quãng đường phanh ngắn nhất), trong khi đó vẫn đảm bảo tốt tính ổn v định hướng khi phanh và tính dẫn hướng ω ω của ô tô. Tb Fx Trong tính toán động lực học của quá trình phanh ô tô, người ta thường sử dụng Hình 1. Các mô men và lực tác động lên bánh xe. giá trị hệ số bám cho trong các bảng. Các hệ số này được xác định bằng thực nghiệm bằng phương pháp kéo bánh xe bị bó cứng hoàn toàn, nghĩa là khi bánh xe bị trượt lê 100%. Trong quá trình phanh ô tô thường xảy ra sự trượt bánh Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  2. μ(λ) xe tương đối với mặt đường, mà hệ số bám 1,0 của bánh xe với mặt đường lại phụ thuộc rất khô nhiều bởi độ trượt này, do đó làm ảnh hưởng 0,8 đến chất lượng phanh. Đồ thị thực nghiệm chỉ ướt sự thay đổi hệ số bám dọc của bánh xe với mặt 0,6 đường theo độ trượt λ giữa bánh xe và mặt tuyết 0,4 đường (hình 2). Theo [6] hệ số bám dọc bằng không khi lực phanh tiếp tuyến bằng không, 0,2 nghĩa là ứng với lúc chưa phanh. λ0 λ% 0 0,2 0,4 0,6 0,8 1 1.2. Yêu cầu của hệ thống điều khiển ABS Hình 2. Hệ số bám dọc theo λ khi phanh. Hệ thống điều khiển ABS phải đảm bảo độ trượt tương đối giữa bánh xe và mặt đường ở giá trị độ trượt tối ưu λ0 = 0,2 (20%) khi phanh. Khi điều kiện mặt đường thay đổi thì tính phi tuyến của ma sát giữa lốp xe và mặt đường cũng thay đổi theo. II. MÔ HÌNH ĐỘNG HỌC BÁNH XE Ô TÔ Các biểu thức của chuyển động của một trong 4 bánh xe ô tô như (1): Jω = rFx − Tb & (1) mv = −Fx & CT 2 trong đó: m là ¼ khối lượng xe; v là tốc độ của xe; ω là tốc độ bánh xe; Fz là lực pháp tuyến; Fx là lực ma sát. Tb là mô men phanh; r là bán kính bánh xe; J là mô men quán tính. v − ωr λ= Độ trượt của bánh xe được định nghĩa như sau [5]: , (2) v Fx = Fzμ(λ, μH, α, Fz, v), Lực ma sát bánh xe Fx được định nghĩa bởi: (3) ở đây μ(λ, μH, α, Fz, v) là hệ số ma sát giữa lốp xe và mặt đường, đây là một hàm phi tuyến với một kiểu phụ thuộc vào độ trượt như hình 2 [8], μH là hệ số ma sát lớn nhất và thay đổi theo điều kiện mặt đường, α là góc lái. Chúng ta chỉ xét trường hợp không có góc lái (α = 0). Thiết kế mô hình: Từ các biểu thức chuyển động, với quan niệm giá trị vận tốc của xe biến đổi chậm hơn rất nhiều so với sự thay đổi của các giá trị khác ở trên, động học của độ trượt bánh xe như sau: r 2 Fz r & λv = Tb − μ, (4) J J với tác động trễ điều khiển một thời gian T, mô hình ABS sẽ có cấu trúc như hình 3 và có thể λ(t )v = −βμ (λ (t )) + αu (t − T ) , & tổng hợp theo biểu thức sau: (5) Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  3. r 2 Fz r α = ;β = ở đây v là một hằng số nhưng không chắc chắn, trong đó: (6) J J β.k βμ(λ) λ 1 u Tb λ 1 Tb u e −sT α e −sT α vs vs Hình 4. Thiết kế mô hình đối tượng ABS_TT. Hình 3. Thiết kế mô hình đối tượng ABS. f (λ ) Theo [4] xấp xỉ khâu trễ bằng khâu quán tính bậc nhất với T = τ. u 1 1 Tb α Thành phần phi tuyến chưa biết sẽ Ts + 1 λ vs có dạng (7): d(t ) Ts + 1 f (λ ) = βμ (λ ) (7) Ts + 1 α β.k α Theo [8] đặc tính ma sát giữa lốp xe và mặt đường như hình 2. Hình 5. Mô hình đối tượng ABS_TT, cộng PT, nhiễu. Thiết kế cho trường hợp mặt đường nhựa khô. Xét phần tuyến tính ma sát: μ(λ ) ≈ k.λ (t ) Để tính toán hàm truyền tuyến tính của ABS ta biến đổi sơ đồ ở hình 4 tương đương với sơ CT 2 đồ hình 5. c WABS _ TT (s ) = (8) 2 + as + b s α T.β.k β.k trong đó: c = ;a= ; b= T.v T.v T.v Các tham số của xe ô tô [7]: ⎡N⎤ J = 1,0 [kg.m2]; m = 450 [kg]; r = 0,32 [m]; Fz = 4414 [N]; β = 451,584 ⎢ ⎥ ⎣ kg ⎦ ⎡1⎤ ⎡ Km ⎤ ⎡m⎤ ⎥ ; τ = T = 14 [ms] = 0,014 [s]; v = 126 ⎢ h ⎥ = 35⎢ s ⎥ ; α = 0,32 ⎢ ⎣ kg.m ⎦ ⎣ ⎦ ⎣⎦ với đặc tính μ(λ) như hình 2, lấy tuyến tính đoạn k = 4,5, thay các giá trị vào biểu thức tham số của đối tượng ta được: c = 0,6531; a = 129,4894; b = 4147,2. Có nhiều phương pháp để tính tham số PID cho đối tượng này [4]. Giả sử tham số bộ PID như sau: kP = 2,5808.103; kI = 1,8434.105, kD = 10. (9) Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  4. III. TỔNG HỢP BỘ ĐIỀU KHIỂN THÍCH NGHI TRÊN CƠ SỞ ĐIỀU KHIỂN MỜ Đối tượng điều khiển phi tuyến có dạng [3]: λ = −aλ − bλ + cu + c[f (λ ) + d(t )] , && & (10) trong đó: f (λ ) là hàm phi tuyến trơn không rõ và tín hiệu nhiễu d(t) không rõ có giới hạn trên cho trước, λ và u lần lượt là tín hiệu ra vào của hệ thống. Khi chưa xét đến thành phần phi tuyến thì đối tượng có dạng (8), thường được điều khiển bằng bộ PID kinh điển: (t ) = u 0 (t ) = k p e(t ) + k i ∫ e(t )dt + k d e(t ) (11) u & PID Có nhiều phương pháp để tổng hợp bộ điều khiển PID [3], [1]. 3.1. Xây dựng thuật toán điều khiển thích nghi mờ với bộ đánh giá TSK cho ABS Phương pháp tiến hành bao gồm hai bước: Thuật toán tổng quát với đối tượng phần tuyến tính có dạng khâu bậc hai, phần phi tuyến chưa biết đã được trình bày trong [2]. Sau đây chúng tôi sẽ áp dụng cho trường hợp điều khiển hệ thống ABS. Bước 1: Thiết kế bộ PID cho đối tượng danh định (không xét đến thành phần phi tuyến). Bước 2: Thiết kế mạch điều khiển bù phi tuyến, để đáp ứng được chất lượng điều khiển thì bộ điều khiển phải được cập nhật thay đổi cho phù hợp với sự ảnh hưởng của phi tuyến và nhiễu tác động nên hệ thống. CT 2 Mô hình mờ TSK sau là thích hợp cho việc mô tả hàm phi tuyến f (λ ) : R i : If λ is A i Then f i = t i λ + v i ; i = 1,2,..., M . Nếu sử dụng bộ mờ hoá singleton và bộ giải mờ trung bình trọng tâm thì hàm phi tuyến f (λ ) có thể xấp xỉ với độ chính xác bất kỳ bằng đánh giá: M ∑ μi fi ˆ (λ ) = i=1 f (12) M ∑μ i=1 i trong đó μ i là độ tin cậy của luật thứ i. f (λ ) = φ T (t ).θ(t ) ˆ Viết lại công thức (12) như sau: (13) [ ] 1 ψ T (t ) = [λ(t ) 1] trong đó: φ T (t ) = μ1ψ T (t ) μ 2 ψ T (t ) ... μ M ψ T (t ) ; (14) M ∑μ i =1 i Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  5. [ ] [ ] θ(t ) = θ1 (t ) θ T (t ) ... θ T (t ) : θi (t ) = t i T vi . và các tham số chưa biết: (15) M 2 3.2. Tổng hợp bộ điều khiển mờ thích nghi Với đối tượng (10) ta xây dựng luật điều khiển như sau: ( ) ( ) 1 1 && u (t ) = λ + aλ + bλ d + u PID (t ) − f (λ ) − f bu (t ) = d + u (t ) − fˆ (λ ) − f bu (t ) ˆ & (16) d d c1 PID c trong đó f (λ ) = φ T (t ).θ(t ) là bộ đánh giá mờ; e(t ) = λ d (t ) − λ (t ), λ d (t ) là độ trượt mong ˆ muốn; f bu (t ) là thành phần bù nhiễu. ˆ && & && & Thay (16) vào (10) ta có: λ = −aλ − bλ + λ + aλ + bλ d + u PID − cf − cf bu + cf + cd , d d ( )( ) ( )( ) ˆ && + k + a e + k + b e + k ∫ e = −c f − f − c d − f hay (17) e & D P I bu ⎡ ∫ e(t )⎤ Định nghĩa véc tơ sai số bám: E(t ) = ⎢ e(t ) ⎥ và véc tơ sai số đánh giá θ e (t ) = θ∗ − θ(t ) , ⎥ ⎢ ⎢ e(t ) ⎥ &⎦ ⎣ ∗ là véc tơ tham số tối ưu của bộ đánh giá mờ TSK. trong đó θ E(t ) = A e E + B f + B d − B c Fbu , (18) & Khi đó phương trình trên có thể viết lại như sau: CT 2 ⎡0 ⎤ ⎡0⎤ ⎡ ⎤ ⎡0⎤ 1 0 0 trong đó: ⎢ ⎥ ⎥ , B = ⎢ 0 ⎥, B = ⎢ 0 ⎥ , Bf = ⎢ Ae = ⎢ 0 0 1 0 ⎥ ⎥c⎢⎥ ⎢ ⎥d⎢ () ⎢− k ⎥ ˆ ⎢− c f − f ⎥ ⎢− cd ⎥ ⎢− c⎥ − b − kP − a − kD ⎦ ⎣⎦ ⎣ ⎦ ⎣ ⎦ ⎣I 1T 1T Chọn hàm Lyapunov xác định dương như sau: V(E, θ e ) = E PE + θe θe , (19) 2γ 2 ⎡ P ⎤ ⎡p p12 p13 ⎤ ⎢ 1 ⎥ ⎢ 11 ⎥ trong đó γ > 0 và P = ⎢P2 ⎥ = ⎢p 21 p 22 p 23 ⎥ là ma trận đối xứng xác định dương thoả ⎢ P ⎥ ⎢p ⎥ ⎣ 3 ⎦ ⎣ 31 p 32 p 33 ⎦ A T P + PA e = −Q mãn phương trình Lyapunov: (20) e với Q là ma trận đối xứng xác định dương chọn trước. Khi biết các tham số bộ điều khiển PID, ta có thể tìm được ma trận P từ (20). ( ) 1& &⎟ 1 Lấy đạo hàm V(E, θ e ) ta có: V E, θ E = ⎛ E T PE + E T PE ⎞ + θ T θ & & ⎜ ⎠ γee ⎝ 2 Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  6. ( ) 1 1 = − E T QE − c f − f P3 E − c.d.P3 E + θ T θ + cf bu P3 E ˆ & (21) ee γ 2 ( ) () () () ∗ T T ˆ ˆ ˆ ˆ Vì có thể viết (22) f = f opt − f opt − f = φ t θ − φ t θ e t nên biểu thức (21) có dạng: ) ( 1 1 V (E, θ e ) = − E T QE − c f − φ T θ∗ + φ T θ e P3 E − cdP3 E + θ T θ + cf bu P3 E & & (23) ee γ 2 ( ) Chọn luật thích nghi cho bộ đánh giá TSK: θ (t ) = −θ (t ) = γcφ (t )P3 E , γ > 0 , & & (24) e () ( ) f bu = − D u + ε .sign P3 E và thành phần bù nhiễu bất định f bu sau: (25) d ≤ Du , ˆ f − f opt ≤ ε trong đó: (26) D u là giá trị xác định, chính là giới hạn trên của nhiễu tác động vào hệ thống, còn ε là một hệ số, ý nghĩa của nó chính là sai số cho phép khi tính toán nhận dạng thành phần phi tuyến [ ()] Khi đó: V (E , θ e ) = − 1 E T QE + c ⎡− ε .sign(P E ) − ⎛ f − fˆ ⎤ ⎞ P E + c − D u sign P E − d P E (27) & ⎜ opt ⎟⎥ 3 ⎢ 3 3 3 ⎝ ⎠⎦ ⎣ 2 Từ (26),(27) ta xác định được: ⎧ 1T [ ] ⎡⎛ ⎞⎤ u ⎪ − E QE − c ⎢d + ⎜ f − fˆopt ⎟⎥ P3 E − c D + ε P3 E < 0 if P3 E > 0 ⎣⎝ ⎠⎦ 2 ⎪ CT 2 (28) V (E , θ e ) = ⎨0 & if P3 E = 0 ⎪ [ ] ⎤ ⎡ 1T ⎪ − E QE − c ⎢d + ⎛ f − fˆopt ⎞⎥ P3 E + c D + ε P3 E < 0 if P3 E < 0 u ⎜ ⎟ ⎣⎝ ⎠⎦ ⎩ 2 () Khi có sự đổi dấu qua bề mặt P3 E = 0 thì sign P3 E trong (25) để tính thành phần bù nhiễu bất định có thể thay bằng hàm sat như sau: ⎧ P3 E ⎪ −1 ≤ −1 if φb ⎛ ⎞⎪ E⎟ ⎪ E sat ⎜ P3 − 1 < P3 E =P ≤1 (29) if ⎜ φ ⎟ ⎨ 3 φb φb b⎠ ⎪ ⎝ P3 E ⎪1 >1 if φb ⎪ ⎩ trong đó φ b là độ mỏng của giá trị mờ hoá tại 0, khi đó (25) trở thành: ( ) ⎛ ⎞ f bu (t ) = m D u + ε sat⎜ P3 E ⎟ (30) ⎜ φb ⎟ ⎝ ⎠ Như vậy theo (19) thì V(.) là dương, V (.) là âm (28) do vậy thuật toán đã tổng hợp được là & ổn định. Sơ đồ cấu trúc của hệ thống điều khiển thích nghi ABS với đánh giá TSK như hình 6. Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  7. f (.) Đánh giá TSK f (.) ˆ Phần phi tuyến f (.) (− ) ˆ λd = λ0 e(t ) λ(t ) u (t ) Đối tượng phần 1 Bộ ĐK PID tuyến tính (− ) (−) c Tính d1 d (t ) Tính toán bù f bu nhiễu fbu Hình 6. Sơ đồ cấu trúc điều khiển ABS thích nghi trên cơ sở điều khiển mờ. IV. KẾT QUẢ MÔ PHỎNG KIỂM CHỨNG HIỆU QUẢ CỦA THUẬT TOÁN Hình 7 là kết quả khi phanh có ABS lúc xe λ0 chạy trên đường nhựa khô λ(t) v[m/s] với đối tượng ABS tuyến tính, bộ điều khiển PID ω.r[m/s] ban đầu (9). Tại thời điểm 0,5 s Hình 7. Các đáp ứng của hệ thống điều khiển ABS tuyến tính. bắt đầu phanh, λ0 là độ trượt tối ưu mong muốn, λ(t) là đáp ứng đầu ra của hệ thống, trên hình 7 ta thấy rằng giá trị này luôn bám sát theo λ0. Tốc độ của xe v[m/s] và tốc độ dài của bánh xe ω.r [m/s] luôn cách nhau CT 2 một khoảng bằng nhau để đạt được độ trượt ở giá trị tối ưu khi phanh. Khi phanh mà xe chạy trên điều kiện mặt đường thay đổi từ đường nhựa khô, sang đường tuyết và sau đó sang đường nhựa ướt. Kết quả như hình 8, từ kết quả đó ta thấy rằng nếu chỉ sử dụng tham số bộ PID ban đầu thì không đảm bảo được độ trượt tối ưu (λ(t) không bám theo λ0) và khi đó sẽ làm cho khả năng chệch hướng của xe ô tô. λ(t) v[m/s] μH-khô μH-ướt ω.r[m/s] λ0 μH-tuyết Hình 8. Các đáp ứng của hệ thống điều khiển ABS phi tuyến sử dụng PID. Các kết quả mô phỏng khi áp dụng thuật toán điều khiển ABS thích nghi mờ - Thiết kế bộ đánh giá mờ: bộ mờ dùng để nhận dạng hệ số ma sát mặt đường. Cấu trúc bộ mờ như hình 9.a; mờ hoá đầu vào như hình 9.b; giá trị rõ ở đầu ra bộ TSK như hình 9.c. Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  8. Hình 9.a. Bộ đánh giá mờ TSK. Hình 9.c. Giá trị rõ ở đầu ra: Dry = 1; Wet-Dry = 0,75; Wet = 0,5; Icy- Wet = 0,25; Icy = 0 Hình 9.b. Mờ hoá đầu vào. Luật mờ: IF (lamda is Icy) THEN (nguy(lamda) is Icy) IF (lamda is Icy-Wet) THEN (nguy(lamda) is Icy-Wet) IF (lamda is Wet) THEN (nguy(lamda) is Wet) IF (lamda is Wet-Dry) THEN (nguy(lamda) is Wet-Dry) IF (lamda is Dry) THEN (nguy(lamda) is Dry) - Các kết quả mô phỏng: Khi áp dụng thuật toán điều khiển thích nghi với bộ đánh giá TSK đã tổng hợp ở trên thì kết quả là đảm bảo được độ trượt tối ưu khi phanh, xem hình 10. v[m/s] CT 2 μH-khô λ(t) μH-ướt λ0 ω.r[m/s] μH-tuyết Hình 10. Các đáp ứng của hệ thống điều khiển mờ thích nghi ABS phi tuyến. Kết quả khi cho tốc độ của xe ô tô thay đổi trong quá trình phanh như hình 11. Nhìn vào kết quả đó ta thấy rằng với thuật toán thích nghi trên cơ λ(t) v[m/s] sở lôgíc mờ thì vẫn đảm bảo λ0 được độ trượt tối ưu khi phanh, có nghĩa là độ trượt ω.r[m/s] đầu ra λ(t) vẫn bám sát được giá trị trượt tối ưu λ0 = 0,2. Hình 11. Các đáp ứng của hệ thống điều khiển mờ thích nghi ABS Còn khi chỉ sử dụng bộ PID khi v thay đổi. ban đầu trong trường hợp này thì kết quả còn kém hơn rất nhiều (hình 12) so với trường hợp khi coi tốc độ của xe không đổi (hình 8). Thể hiện ở λ(t) không bám được giá trị trượt tối ưu λ0 = 0,2. Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
  9. VI. KẾT LUẬN Nhìn vào các đáp ứng của hệ thống điều khiển ABS chúng ta thấy v[m/s] λ0 rằng với phương pháp tổng hợp mà bài báo này λ(t) ω.r[m/s] đưa ra đạt được chất lượng điều khiển rất tốt. Hình 12. Các đáp ứng của hệ thống điều khiển ABS khi v thay đổi Đạt được kết quả này chỉ sử dụng PID. chính là do sử dụng bộ TSK để tự động nhận dạng ma sát mặt đường, trên cơ sở đánh giá ma sát này, thuật toán sẽ tự tính toán lượng điều khiển thích nghi. Tài liệu tham khảo [1]. Cao Tiến Huỳnh, Đào Tuấn, Trần Quang Oánh, Nguyễn Văn Tiềm (2004). “Điều khiển mờ thích nghi áp dụng cho đối tượng chuyển động”, Chuyên san Kỹ thuật điều khiển tự động, Tự động hoá ngày nay, Hội khoa học công nghệ tự động Việt Nam, tr. 16-22. [2]. Lê Hùng Lân, Nguyễn Văn Tiềm (2005). “Xây dựng thuật toán điều khiển mờ thích nghi áp dụng để CT 2 điều khiển đối tượng chuyển động trên cơ sở bộ đánh giá TSK”, Hội nghị khoa học kỹ thuật đo lường toàn quốc lần thứ IV, Tuyển tập báo cáo khoa học, NXB KHKT, tr. 666 –671. [3]. Lê Hùng Lân, Nguyễn Văn Tiềm, Trần Quang Oánh (2002), “Điều khiển thích nghi gián tiếp chuyển động trên cơ sở các bộ xấp xỉ mờ”, Tuyển tập các báo cáo khoa học, Hội nghị toàn quốc lần thứ 5 về Tự động hoá, tr. 289-294. [4]. Nguyễn Doãn Phước (2005), Lý thuyết điều khiển tuyến tính, NXB KHKT. [5]. CHIH-KENG CHEN, MING-CHANG SHIH, “PID-Type Fuzzy Control for Anti-Lock Brake Systems with Parameter Adaptation”, JSME International Journal, Series C, Vol. 47, No. 2, (2004), pp.675-685. [6]. FANGJUN JIANG, ZHIQIANG GAO, “An application of Nonlinear PID Control to a Class of Truck ABS Problems”, academic.csuohio.edu/aerl/papers/cdc01_abs.pdf . [7]. TOR A. JOHANSEN, IDAR PETERSEN, JENS KALKKUHL and JENS LÜDEMANN, “Gain- scheduled Wheel Slip Control in Automotive Brake Systems”, ieeexplore.ieee.org/iel5/87 /28090/01255656.pdf. [8]. WEI-EN TING and JUNG-SHAN LIN, “Nonlinear Control Design of Anti-lock Braking Systems Combined with Active Suspensions”,ieeexplore.ieee.org/iel5/9768/30803/01426017.pdf♦ Số 21 - 03/2008 Tạp chí KHOA HỌC GIAO THÔNG VẬN TẢI
nguon tai.lieu . vn