Lecture Note 5.docx

HỌC VIỆN KỸ THUẬT QUÂN SỰ KHOA HÀNG KHÔNG VŨ TRỤ Bản số: 05

BÀI GIẢNG Môn học: Phân tích và tổ hợp hệ thống Cơ điện tử Bài: Tổng hợp bộ điều khiển PID Đối tượng: Cơ điện tử 50 Năm học: 2017-2018

Thiếu tá, ThS Hoàng Văn Tiến

HÀ NỘI, THÁNG 08 NĂM 2018

HỌC VIỆN KỸ THUẬT QUÂN SỰ KHOA HÀNG KHÔNG VŨ TRỤ

PHÊ DUYỆT Ngày tháng năm 2018 CHỦ NHIỆM KHOA

Đại tá, PGS, TS Đặng Ngọc Thanh

BÀI GIẢNG Môn học: Phân tích và tổ hợp hệ thống Cơ điện tử Bài: Tổng hợp bộ điều khiển PID Đối tượng: Cơ điện tử 50 Năm học: 2017-2018

Ngày tháng năm 2018 CHỦ NHIỆM BỘ MÔN

Thượng tá, TS Hoàng Quang Chính

HÀ NỘI, THÁNG 08 NĂM 2018

2

MỞ ĐẦU

Bài giảng này giúp sinh viên nắm được các phương pháp tổng hợp bộ điều khiển PID: - Khái niệm và vai trò bộ điều khiển PID - Tổng hợp bộ điều khiển điều khiển PID + Tổng hợp bộ điều khiển tỷ lệ P + Tổng hợp bộ điều khiển tỷ lệ - tích phân PI + Tổng hợp bộ điều khiển tỷ lệ - vi phân PD + Tổng hợp bộ điều khiển PID - Các phương pháp chỉnh định thông số bộ điều khiển PID + Phương pháp Ziegler-Nichols 1 + Phương pháp Ziegler-Nichols 2 + Phương pháp Cohen-Coon + Phương pháp phản hồi rơ-le

3

NỘI DUNG

Chương 5 TỔNG HỢP BỘ ĐIỀU KHIỂN PID

Bộ điều khiển PID là một trong những bộ điều khiển phổ biến nhất trong các lĩnh vực điều khiển công nghiệp vì chúng có cấu trúc đơn giản và làm việc tin cậy. Tuy nhiên việc xác định các thông số cho bộ điều khiển đòi hỏi phải có sự hiểu biết và chi phí một khoảng thời gian nhất định. Chính vì lý do đó trong hầu hết các thiết bị điều khiển công nghiệp thường được các hãng sản xuất tích hợp các chương trình tự động điều chỉnh các thông số. 2.1. CƠ SỞ CỦA ĐIỀU KHIỂN PID

Mặc dù những tiến bộ của lý thuyết và các phương pháp thiết kế đang liên tục phát triển trong lĩnh vực điều khiển tự động, các bộ điều khiển Tỉ lệ-Tích phân-Vi phân (PID - Proportional, Integral, Derivative control) vẫn được sử dụng rộng rãi nhất trong công nghiệp do ưu điểm về tỉ số giá thành và lợi nhuận của chúng đem lại. Thực tế, mặc dù chúng tương đối đơn giản trong việc sử dụng, nhưng vẫn có khả năng thỏa mãn về thực thi trong nhiều tác vụ điều khiển quá trình. Thực ra quá trình lịch sử lâu dài và các kiến thức được để lại qua nhiều năm đã làm cho cách sử dụng chúng như một bộ điều khiển phản hồi chuẩn mực. Tuy nhiên khả năng đáp ứng cao của các bộ vi xử lý, các công cụ phần mềm và sự tăng yêu cầu nâng cao chất lượng sản phẩm đồng thời giảm giá thành sản phẩm đã thúc đẩy các nhà nghiên cứu phát minh ra các phương pháp mới cho việc cải tiến khả năng hoạt động và sự đơn giản trong sử dụng các bộ điều khiển. Các bộ điều khiển PID đã được phát triển trong quá trình phát triển công nghệ, và ngày nay nó được hiện thực rất phổ biến dưới dạng số với các phần tử điện hoặc khí nén,... Nó có thể thấy trong hầu hết các loại thiết bị điều khiển như một bộ điều khiển đơn - độc lập hoặc dưới dạng các khối chức năng trong các bộ

4

điều khiển PLC (Programable Logic Controller) hoặc DCS (Distributed Control System). Các bộ điều khiển PID là công cụ chuẩn dùng cho tự động hóa công nghiệp. Sự linh hoạt của bộ điều khiển tạo cho nó khả năng sử dụng trong nhiều trường hợp. Các bộ điều khiển cũng có thể được sử dụng trong điều khiển tầng và các cấu hình bộ điều khiển khác. Nhiều vấn đề điều khiển đơn giản có thể được giải quyết rất tốt bởi điều khiển PID - với những yêu cầu chất lượng không quá cao. Thuật toán PID được đóng gói trong các bộ điều khiển tiêu chuẩn cho điều khiển quá trình và cũng là cơ sở cho nhiều hệ thống điều khiển đơn giản. Thuật toán PID được biểu diễn dạng công thức gồm 3 thành phần :

1t de u (t )  K c (e(t )   e( s)ds  Td ) Ti 0 dt

(2.1)

Trong đó: - u là biến điều khiển; - e = r – y là sai số; - r là giá trị đặt (giá trị mong muốn); - y là đầu ra của quá trình (biến đầu ra). Bộ điều khiển PID số có dạng như sau:

T t k u(tk )  K P e(tk )   e(ti )  d  e(tk )  e(tk 1 )  Ti i 1 t Thuật toán thực tế sử dụng gồm có nhiều sửa đổi. Nó thông thường là tiêu chuẩn cho phép thành phần vi phân chỉ tác dụng tới đầu ra quá trình. Có thể rất có lợi khi để thành phần tỉ lệ hoạt động chỉ như trong một phần của giá trị đặt. Tác dụng của thành phần vi phân được thay thế bởi một đại lượng xấp xỉ có thể giảm độ tăng tại tần số cao. Thành phần tích phân được điều chỉnh sao cho nó không được tiếp tục tích lũy khi bão hòa các biến điều khiển (anti-windup). Sự điều chỉnh này có hiệu lực, và sẽ hạn chế sự thay đổi nhanh chóng khi bộ điều 5

khiển chuyển từ chế độ bằng tay sang tự động hoặc khi các tham số được thay đổi. 2.1.1. Tổng hợp các bộ điều khiển PID Trong phần này trình bày việc tổng hợp các bộ điều khiển khác nhau như P, PI, PD và PID. Cấu trúc chung của bộ điều khiển PID được thể hiện trên hình 2.1 dưới đây:

Hình 2.1. Cấu trúc điển hình của bộ điều khiển PID a. Tổng hợp bộ điều khiển P Sơ đồ khối hệ điều khiển P được trình bày như sau:

r(t)

e(t)

u(t) Kp

§ èi t- î ng §K

y(t)

Hình 2.2. Sơ đồ khối bộ điều khiển P Từ hình vẽ 2.2 ta có: u (t )  Kp * e(t )

Để minh họa quá trình làm việc của bộ điều khiển, trong các phần tổng hợp bộ điều khiển, chúng ta tiến hành cho đối tượng điều khiển có hàm truyền đạt như sau: 6

G( s) 

1 ( s  1)3

Thực hiện tổng hợp bộ điều khiển P trên MATLAB với đoạn code sau: ---------------------------------------------------------------G =tf(1, [1 3 3 1]) % Nhập hàm truyền đạt đối tượng for Kp=0.1:0.2:2 % Vòng FOR thực hiện thay đổi các hệ số Kp G_c = feedback(Kp*G,1) % Tạo vòng điều khiển hệ kín step(G_c); % Lấy đáp ứng đơn vị ứng với từng Kp Hold on % Cho phép vẽ tất các các đáp ứng trên cùng 1 đồ thị end ---------------------------------------------------------------Từ đoạn code trên ta nhận được 1 họ đường đáp ứng đầu ra của hệ thống với bộ điều khiển tỉ lệ. Các đáp ửng này được biểu diễn trên hình 2.3. Từ hình vẽ 2.3 ta có những nhận xét sau: - Bộ điều khiển P là bộ điều khiển đơn giản - Khi thay đổi hệ số khuếch đại Kp ta thấy o Sai số tĩnh giảm o Tác động nhanh của hệ thống tăng hay thời gian tăng giảm. o Độ quá chỉnh tăng o Không có ảnh hưởng nhiều đến thời gian quá độ

7

Step Response 0.9

0.8

Kp=1.9 0.7

0.6

Amplitude

0.5

0.4

0.3

0.2

0.1

Kp=0.1 0

0

2

4

6

8

10

12

14

16

Time (sec)

Hình 2.3. Kết quả đáp ứng của hệ khi thay đổi hệ số tỷ lệ (Kp)  Chứng minh ảnh hưởng của thành phần tỷ lệ: • Sai số tĩnh giảm

Hình 1. Sơ đồ khối hệ thống điều khiển vòng kín sử dụng BĐK có hàm truyền C(s) - Biến đổi Laplace của sai số được xác định như sau:

E s  R s  Y s Y  s   E(s)C(s)G  s 

(Error! No text of specified style in document.1) (2)

8

Thế (2) vào (1) thu được:

E  s   R  s   E  s C s G s  E  s  1  C  s  G  s   R  s 

E s 

R s 1  C sG s

- Sai số tĩnh được xác định như sau:

sR  s  s 0 1  C  s  G  s 

ess  lim sE  s   lim s 0

Với tín hiệu đầu vào r(t) là tín hiệu bậc thang r  t   Rus  t  với us  t  là hàm bậc thang đơn vị, thì biến đổi Laplace của r(t) là R  s   R s , do đó:

sR s R R  lim  s 0 1  C  s  G  s  s 0 1  C  s  G  s  1  lim C  s  G  s 

ess  lim

s 0

Bộ điều khiển tỷ lệ có hàm truyền C  s   K p do đó:

ess 

R R  1  lim  K pG  s  1  K p lim G  s  s 0

(3)

s 0

Từ biểu thức (3) có thể thấy rằng khi tăng K p thì sai số tĩnh giảm và ngược lại. • Tác động nhanh của hệ thống tăng hay thời gian tăng giảm. • Độ quá chỉnh tăng • Không có ảnh hưởng nhiều đến thời gian quá độ

9

Hình 2. Ảnh hưởng của thành phần tỷ lệ tới chất lượng động học của hệ thống y(t) là tín hiệu đầu ra, e(t) là sai số giữa tín hiệu đặt và tín hiệu phản hồi, u(t) là tín hiệu đầu ra của bộ điều khiển tỷ lệ Để giải thích ảnh hưởng của thành phần tỷ lệ tới chất lượng động học của hệ thống, giả sử rằng hệ thống chứa động cơ điện một chiều với mô men tỷ lệ với tín hiệu u(t). Khi đó ảnh hưởng của thành phần tỷ lệ được giải thích như sau: 1. Trong khoảng thời gian 0  t  t1 : sai số e  t  dương. Mô men động cơ là dương và tăng nhanh, tạo ra gia tốc lớn, vận tốc tăng nhanh và đầu ra y(t) tăng nhanh. Do u  t   K p e  t  và mô men tỷ lệ thuận với u(t), nên K p càng lớn thì mô men càng lớn  gia tốc càng lớn, vận tốc tăng càng nhanh  y(t) tăng càng nhanh  thời gian tăng giảm và hệ thống đáp ứng càng nhanh. Ngoài ra, thì do quán tính nên khi y(t) tăng càng nhanh thì độ quá chỉnh lớn. 2. Trong khoảng thời gian t1  t  t3 : sai số e  t  âm, và tín hiệu u(t) và mô men cũng âm  gia tốc âm, vận tốc giảm và dần đổi chiều  đầu ra y(t) tăng đến giá trị cực đại tại t2 , sau đó giảm về giá trị đặt và tiếp tục giảm xuống dưới giá trị đặt do quán tính. 3. Trong khoảng thời gian t3  t  t5 : sai số e  t  lại dương  tín hiệu u(t) và mô men dương  gia tốc dương, vận tốc tăng trở lại  đầu ra y(t) giảm chậm dần đến thời điểm t4 , sau đó tăng trở lại. 4. Sau đó, y(t) tiếp tục dao động với biên độ nhỏ dần cho tới khi đạt trạng thái xác lập. Có thể thấy rằng: nếu K p càng lớn thì u(t) và mô men biến 10

đổi càng nhanh, và y(t) biến thiên nhanh hơn, nhưng số dao động cũng nhiều hơn. Vì vậy, K p không ảnh hưởng nhiều tới thời gian quá độ. b. Tổng hợp bộ điều khiển PI Sơ đồ khối hệ điều khiển PI được trình bày trên hình 2.4. Thực hiện tổng hợp bộ điều khiển PI trên MATLAB với đoạn code sau: ---------------------------------------------------------------G =tf(1, [1 3 3 1]) % Nhập hàm truyền đạt đối tượng Kp=1; s=tf('s'); for Ki=0.7:0.1:1.5

% Vòng FOR thực hiện thay đổi hệ số tích phân

Ki G_c = Kp+Ki/s % Bộ điều khiển PI G_s = feedback(G_c*G,1) %Tạo vòng điều khiển hệ kín step(G_s); % Lấy đáp ứng đơn vị ứng với từng Ti Hold on % Cho phép vẽ tất các các đáp ứng trên cùng 1 đồ thị end ---------------------------------------------------------------Từ đoạn code trên ta nhận được 1 họ đường đáp ứng đầu ra của hệ thống với bộ điều khiển PI. Các đáp ửng này được biểu diễn trên hình 2.5.

r(t)

e(t)

u(t)

Kp

Ki *

§ èi t- î ng §K

y(t)



Hình 2.4. Sơ đồ khối bộ điều khiển PI 11

Step Response 2

1.8

Ki=1.5

1.6

1.4

Amplitude

1.2

1

0.8

0.6

0.4

Ki=0.1

0.2

0

0

2

4

6

8

10

12

14

16

18

20

Time (sec)

Hình 2.5. Kết quả đáp ứng của hệ khi thay đổi hệ số Ki khi Kp = 1 Thêm thành phần tích phân vào bộ điều khiển và thay đổi hệ số khuếch đại Ki, Kp = const, chúng ta thấy đáp ứng của hệ thống có những đặc điểm sau (từ hình vẽ 2.5): o Thành phần tích phân có tác dụng triệt tiêu sai số tĩnh o Tác động nhanh của hệ thống giảm hay thời gian tăng tăng. o Độ quá chỉnh tăng o Thời gian quá độ tăng  Chứng minh ảnh hưởng của thành phần tích phân: • Triệt tiêu sai số tĩnh

Hình 1. Sơ đồ khối hệ thống điều khiển vòng kín sử dụng BĐK có hàm truyền C(s) - Biến đổi Laplace của sai số được xác định như sau: 12

E s  R s  Y s

(1)

Y  s   E(s)C(s)G  s 

(2)

Thế (2) vào (1) thu được:

E  s   R  s   E  s C s G s  E  s  1  C  s  G  s   R  s 

E s 

R s 1  C sG s

- Sai số tĩnh được xác định như sau:

sR  s  s 0 1  C  s  G  s 

ess  lim sE  s   lim s 0

Với tín hiệu đầu vào r(t) là tín hiệu bậc thang r  t   Rus  t  với us  t  là hàm bậc thang đơn vị, thì biến đổi Laplace của r(t) là R  s   R s , do đó:

sR s R R  lim  s 0 1  C  s  G  s  s 0 1  C  s  G  s  1  lim C  s  G  s 

ess  lim

s 0

Bộ điều khiển tỷ lệ có hàm truyền C  s  

ess 

K p s  Ki s

do đó:

R R   K s  Ki  K s  Ki   G  s  1  lim  p  G  s  1  lim  p s 0 s 0 s s  

(3)

Vì hàm truyền G  s  có dạng tổng quát như sau:

trong đó K và tất cả các hằng số T là thực, và j  0,1,2... . Do đó:

K s  K  G s  lim  p

s 0

i

s

và

ess 

R 0 1 

(3)

Như vậy, thành phần tích phân làm triệt tiêu sai số tĩnh. 13

 Tác động nhanh của hệ thống giảm hay thời gian tăng tăng.  Độ quá chỉnh tăng  Thời gian quá độ tăng c. Tổng hợp bộ điều khiển PD Sơ đồ khối hệ điều khiển PD được trình bày trên hình 2.6.

r(t)

e(t)

u(t)

Kp

Kd *

§ èi t- î ng §K

y(t)

d dt

Hình 2.6. Sơ đồ khối bộ điều khiển PD Thực hiện tổng hợp bộ điều khiển PI trên MATLAB với đoạn code sau: ---------------------------------------------------------------G =tf(1, [1 3 3 1]) % Nhập hàm truyền đạt đối tượng Kp=1; s=tf('s'); for Kd=0.1:0.2:2 % Vòng FOR thực hiện thay đổi hế số Kd G_c = Kp+Kd*s % Bộ điều khiển PD G_s = feedback(G_c*G,1) %Tạo vòng điều khiển hệ kín step(G_s); % Lấy đáp ứng đơn vị ứng với từng Td Hold on % Cho phép vẽ tất các các đáp ứng trên cùng 1 đồ thị End ---------------------------------------------------------------Từ đoạn code trên ta nhận được 1 họ đường đáp ứng đầu ra của hệ thống với bộ điều khiển PD. Các đáp ửng này được biểu diễn trên hình 2.7. 14

Step Response 0.6

Kd=2 0.5

0.4

Amplitude

Kd=0.1

0.3

0.2

0.1

0

0

1

2

3

4

5

6

7

8

9

10

Time (sec)

Hình 2.7. Kết quả đáp ứng của hệ khi thay đổi hệ số Kd khi Kp = 1 Thêm thành phần vi phân vào bộ điều khiển và thay đổi hệ số khuếch đại Kd, Kp=const chúng ta thấy đáp ứng của hệ thống có những đặc điểm sau (từ hình vẽ 2.7): o Không ảnh hưởng nhiều đến sai số tĩnh o Không ảnh hưởng nhiều đến tính tác động nhanh của hệ thống. o Độ quá chỉnh giảm o Thời gian quá độ giảm o Nhạy cảm với nhiễu  Chứng minh ảnh hưởng của thành phần vi phân: • Không ảnh hưởng tới sai số tĩnh Ở trạng thái tĩnh có sai số tính ess  constant , khi đó thành phần tín hiệu điều khiển tạo ra bởi thành phần vi phân như sau: 15

uKd  t   K d

dess  Kd  0  0 dt

(1)

Như vậy, thành phần vi phân không ảnh hưởng tới sai số tĩnh. • Nhạy cảm với nhiễu Vì nhiễu thường dao động với tần số lớn nên đạo hàm của nhiễu cũng là dao động với tần số lớn và biên độ tương ứng. Do đó, thành phần vi phân nhạy cảm với nhiễu và khuếch đại nhiễu. • Không ảnh hưởng nhiều đến tính tác động nhanh của hệ thống. • Độ quá chỉnh giảm Để giải thích ảnh hưởng của thành phần vi phân tới chất lượng động học của hệ thống, giả sử rằng hệ thống chứa động cơ điện một chiều với mô men tỷ lệ với tín hiệu u(t). Khi đó ảnh hưởng của thành phần tỷ lệ được giải thích như sau: 5. Trong khoảng thời gian 0  t  t1 : sai số e  t  dương và giảm, do đó: Thành phần vi phân tạo ra phần tín hiệu điều khiển:

uKd  t   Kd

de  t   0, dt

và thành phần tỷ lệ tạo ra tín hiệu: uKp  t   K p e  t   0 Do đó tín hiệu điều khiển u  t  tạo ra bơi BĐK PD là: u  t   uKp  t   uKi  t   0

Khi đó mô men động cơ là dương và tăng nhanh, tạo ra gia tốc lớn, vận tốc tăng nhanh và đầu ra y(t) tăng nhanh. Do đó, khi tăng K d  u  t  giảm (do

uKd  t  ngược dấu với uKp  t  )  mô men giảm  gia tốc giảm, vận tốc tăng chậm lại  y(t) tăng chậm hơn  thời gian tăng tăng và hệ thống đáp ứng càng nhanh, nhưng tác động này của thành phần vi phân không lớn .

16

Hình 2. Ảnh hưởng của thành phần vi phân tới chất lượng động học của hệ thống y(t) là tín hiệu đầu ra, e(t) là sai số giữa tín hiệu đặt và tín hiệu phản hồi, uKd  t  là tín hiệu đầu ra do thành phần vi phân gây ra 6. Trong khoảng thời gian t1  t  t2 : sai số e  t   0 âm, uKp  t   0 ;

uKd  t   0 . Do đó, trong khoảng thời gian này uKd  t  làm tăng mức âm của u  t   uKp  t   uKd  t  , tăng mô men cản  tăng gia tốc âm  vận tốc giảm và dần đổi chiều nhanh  giá trị cực đại của y(t) tại t2 giảm. Như vậy, cùng với tác động làm giảm thời gian tăng trong khoảng 0  t  t1 , tác động thành phần vi phân trong khoảng thời gian này làm giảm độ quá điều chỉnh và khi tăng hệ số K d thì độ quá chỉnh giảm.

7. Trong khoảng thời gian t2  t  t3 : sai số e  t   0 và tăng, khi đó: uKp  t   0 , uKd  t   0 . Do đó, trong khoảng thời gian này uKd  t  làm

giảm mức âm của u  t   uKp  t   uKd  t   giảm mô men cản  giảm gia tốc hãm  vận tốc giảm chậm về 0  y(t) giảm chậm về giá trị đặt, đồng thời độ vọt dưới điều chỉnh giảm. 8. Trong khoảng thời gian t3  t  t4 : sai số e  t   0 và tăng, khi đó: uKp  t   0 , uKd  t   0 . Do đó, trong khoảng thời gian này uKd  t  làm

tăng u  t   uKp  t   uKd  t   0  tăng mô men tăng tốc  tăng gia tốc dương  vận tốc tăng dần về 0  y(t) giảm nhanh về giá cực tiểu tại t4 . 17

9. Như vậy, thành phân vi phân tỷ lệ với tốc độ biến thiên của sai số nên nó làm giảm mức biến thiên của sai số và đầu ra. Do đó, khi tăng hệ số vi phân K d thì đáp ứng của hệ thống giảm dao động, và tăng giảm chấn, làm cho hệ thống nhanh chóng đạt đến trạng thái xác lập và thời gian xác lập giảm. d. Tổng hợp bộ điều khiển PID

Kp

r(t)

e(t)

Ki *



Kd *

d dt

u(t)

§ èi t- î ng §K

y(t)

Hình 2.6. Sơ đồ khối bộ điều khiển PID Thực hiện tổng hợp bộ điều khiển PI trên MATLAB với đoạn code sau: ---------------------------------------------------------------G =tf(1, [1 3 3 1]) Kp=1;Ki=1; s=tf('s');

% Nhập hàm truyền đạt đối tượng

% Vòng FOR thực hiện thay đổi hằng số vi % phân Td G_c = Kp+Ki/s+Kd*s) % Bộ điều khiển PD

for Kd=0.1:0.2:2

18

G_s = feedback(G_c*G,1) %Tạo vòng điều khiển hệ kín step(G_s); % Lấy đáp ứng đơn vị ứng với từng Td Hold on % Cho phép vẽ tất các các đáp ứng trên cùng 1 đồ thị end ---------------------------------------------------------------Từ đoạn code trên ta nhận được 1 họ đường đáp ứng đầu ra của hệ thống với bộ điều khiển PI. Các đáp ửng này được biểu diễn trên hình 2.9.

Step Response 1.6

1.4

Kd=0.1

1.2

Amplitude

1

0.8

Kd=2 0.6

0.4

0.2

0

0

2

4

6

8

10

12

14

16

18

20

Time (sec)

Hình 2.5. Kết quả đáp ứng của hệ khi thay đổi hệ số Kd khi Kp = 1 và Ki=1 Tóm lại: Từ các nghiên cứu trên chúng ta có thể kết luận sự ảnh hưởng của các thành phần đến chất lượng của hệ kín được trình bày trong bảng 2.1 sau: 19

Bảng 2.1. Tác dụng các tham số P, I, D tới chất lượng điều khiển Thời gian

Độ quá

Thời gian

tăng

điều chỉnh

quá độ

Kp

Giảm

Tăng

Thay đổi nhỏ

Giảm

Ki

Giảm

Tăng

Tăng

Loại trừ

Kd

Thay đổi nhỏ

Giảm

Giảm

Đáp ứng

Sai lệch tĩnh

không thay đổi

2.2. CÁC PHƯƠNG PHÁP CHỈNH ĐỊNH THÔNG SỐ BỘ ĐIỀU KHIỂN PID

Ngày nay với sự phát triển mạnh mẽ của kỹ thuật điều khiển số, nên hầu hết các bộ điều khiển PID đều được thực hiện dưới dạng mạch số cứng hoặc mềm. Từ đó cho phép các nhà thiết kế các bộ điều khiển PID thực hiện các thủ tục tự động chỉnh định thông số cho bộ điều khiển này. Tự động chỉnh định thông số bộ điều khiển PID là quá trình tự xác định các thống số cho bộ điều khiển PID theo một thuật toán nào đó dựa trên đáp ứng của đối tượng được điều khiển. Có nhiều phương pháp tự động điều chỉnh tham số PID đã được đề xuất. Phương pháp phổ biến nhất chính là thực hiện một thí nghiệm đơn giản với đối tượng điều khiển. Thí nghiệm có thể thực hiện trong vòng điều khiển hở hoặc vòng kín. Trong các thí nghiệm với vòng hở, đầu vào của quá trình được kích thích bởi một xung hoặc một cặp bước bậc thang. Một mô hình quá trình đơn giản, ví dụ như khâu bậc hai, sau đó được ước lượng bằng cách sử dụng lặp lại bằng phương pháp bình phương tối thiểu, hoặc vài phương pháp ước lượng khác. Nếu một mô hình quá trình khâu giữ mẫu bậc hai được ước lượng, sau đó bộ điều khiển PID có thể được sử dụng để gán các điểm cực. Tốc độ và sự suy giảm của hệ thống sau đó là do các tham số thiết kế. Một phương pháp thiết kế phổ biến là chọn các điểm không của bộ điều khiển, giống như khi bỏ qua hai

20

điểm cực của quá trình. Điều này cho phép các đáp ứng tốt với sự thay đổi điểm đặt, trong khi đáp ứng tới các nhiễu tải được xác định bởi tác động vòng hở. Phương pháp đáp ứng tức thời cho việc tự động điều chỉnh của các bộ điều khiển PID được sử dụng trong các bộ điều khiển của Yokogawa, Eurotherm(SSD), và Honeywell. Nó được dùng cho việc điều chỉnh trước trong các bộ điều khiển thích nghi của Leeds, Northrup và Turnbull Control. Các thí nghiệm điều chỉnh có thể được thực hiện trong vòng kín. Một ví dụ điển hình là phương pháp tự dao động của Ziegler -Nichols hoặc dạng biến đổi. Phương pháp tự chỉnh rơle dựa trên sự tự dao động được sử dụng trong các sản phẩm của SattControl và Fisher-Rousemount. Trong các bộ điều khiển này, sự điều chỉnh được bắt đầu đơn giản bằng cách bấm nút. Một ưu điểm của việc tạo các thí nghiệm trong vòng kín là đầu ra của quá trình có thể được giữ trong vùng giới hạn hợp lý, mà có thể khó cho các quá trình có thành phần tích phân nếu thí nghiệm được thực hiện trong vòng hở. Chức năng tự động điều chỉnh thường được tích hợp sẵn trong các bộ điều khiển PID. Tự động điều chỉnh cũng có thể được thực hiện bằng cách sử dụng thiết bị ngoài. Bộ điều chỉnh sẽ được nối với một quá trình và thực hiện một thí nghiệm thông thường với mạch hở. Bộ điều chỉnh sau đó sẽ gợi ý các tham số cài đặt, và chúng sẽ được cài vào bộ điều khiển PID một cách tự động hoặc bằng tay. Do yêu cầu bộ điều chỉnh bên ngoài phải có khả năng làm việc với các bộ PID của nhiều nhà sản xuất khác nhau, nên một điều quan trọng là nó phải có thông tin chi tiết về cách thực hiện của thuật toán PID trong từng trường hợp cụ thể. Một phương pháp khác cho việc tự động điều chỉnh là sử dụng một hệ chuyên gia để điều chỉnh bộ điều khiển. Việc này được thực hiện trong thời gian hoạt động bình thường của quá trình. Hệ chuyên gia sẽ đợi cho sự thay đổi của điểm đặt hoặc nhiễu lên tải và sau đó ước lượng sự hoạt động của hệ điều khiển vòng kín. Các thuộc tính như: độ suy giảm, chu kì dao động, và hệ số tăng tĩnh đều được ước lượng. Các tham số bộ điều khiển được thay đổi tùy thuộc tới quy 21

luật xây dựng bên trong, bắt chước hoạt động của một kĩ sư điều khiển có kinh nghiệm. Nhận dạng mẫu hoặc hệ chuyên gia được dùng trong các bộ điều khiển của các hãng như: Foxboro và Fenwal. 2.2.1. Phương pháp xác định thông số bộ điều khiển PID bằng thực nghiệm Sau đây là các bước thực hiện việc xác định các thông số của bộ điều khiển PID bằng tay: - Bước 1: Cho các hệ số Ki=0 và Kd=0 - Bước 2: Tăng Kp cho đến khi đáp ứng đầu ra dao động - Bước 3: Lấy giá trị Kp cho hệ thống gần bằng nửa giá trị nhận được trong bước 2. - Bước 4: Tăng Ki cho đến khi sai số tĩnh đủ nhỏ, tuy nhiên nếu Ki quá lớn gây cho hệ mất ổn định. - Bước 5: Muốn hệ thống đáp ứng tốt như đáp ứng nhanh đạt đến giá trị mong muốn, giảm độ quá chỉnh thì chúng ta cho thêm thành phần vi phân. Tuy nhiên hệ số khuếch đại này thường nhỏ để đảm bảo hệ ổn định vì nó rất nhạy với nhiễu.

2.2.2. Phương pháp xác định thông số bộ điều khiển PID theo ZieglerNichols Ziegler-Nichols đã đưa ra hai phương pháp thực nghiệm để xác định tham số bộ điều khiển PID. Phương pháp thứ nhất sử dụng mô hình xấp xỉ quán tính bậc nhất có trễ của một lớp đối tượng điều khiển. Phương pháp thứ hai có ưu điểm không cần sử dụng tới mô hình toán học của đối tượng. Tuy nhiên nó có hạn chế chỉ áp dụng với một lớp đối tượng nhất định. a. Phương pháp Ziegler-Nichols thứ nhất

22

Một cách đơn giản để xác định các tham số của bộ điều khiển PID dựa trên dữ liệu đáp ứng bậc thang đã được Ziegler-Nichols phát triển và công bố năm 1942. Phương pháp này được gọi là phương pháp đáp ứng bậc thang. Một hệ thống với đáp ứng bậc thang như hình 2.2 ta có thể xấp xỉ hóa bằng hàm truyền : ke  Ls G (s)  1  Ts

(2.6)

trong đó k là hệ số khuyếch đại tĩnh, L là thời gian trễ, và T là hằng số thời gian. Tham số a được tính như sau: ak

L T

(2.7)

Hình 2.2. Đáp ứng bậc thang và các xác định thông số cho bộ điều khiển PID bằng phương pháp Ziegler-Nichols 1 Các tham số của bộ điều khiển PID được Ziegler-Nichols xây đựng từ hai a và L trình bày trong bảng 2: Bảng 2.2. Các công thức xác định theo phương pháp Ziegler-Nichols 1 Bộ điều khiển

Kc

Ti

Td

P

1/a

.....

.....

PI

0.9/a

3L

..... 23

PID

1.2/a

2L

L/2

Quy luật điều chỉnh Ziegler-Nichols được phát triển bởi các mô phỏng kinh nghiệm trong nhiều hệ thống khác nhau. Luật điều chỉnh này có sự hạn chế cho các hệ thống điều chỉnh vòng kín mà sự hạn chế rung động là rất kém. Các hệ thống với sự chống dao động tốt hơn có thể đạt được bằng cách sửa đổi các giá trị số trong bảng trên. Bằng cách sử dụng các tham số phụ ta cũng có thể xác định khi nào quy tắc Ziegler-Nichols có thể áp dụng. Nếu hằng số thời gian T cũng được xác định trước, một quy tắc kinh nghiệm được xác lập khi công thức Ziegler-Nichols có thể ứng dụng nếu 0.1
A0 , vùng A1 tính từ dưới đáp ứng bước nhảy tới thời gian T + L k

được xác định trước, T được tính theo công thức : T = eA1/k ,trong đó e được tính theo logarit tự nhiên. Những hạn chế cơ bản của phương pháp này là có thể khó biết được kích cỡ của bước nhảy trong tín hiệu điều khiển và khó xác định khi nào đạt tới trạng thái ổn định. Bước nhảy có thể là lớn khiến cho đáp ứng được nhận thấy rõ ràng hơn nhiễu, tuy nhiên nó không quá lớn tới mức gây nhiễu cho hệ thống. Các nhiễu sẽ ảnh hưởng một cách đáng chú ý tới kết quả. Ngoài ra việc lấy đáp ứng tín hiệu bậc thang rất dễ bị ảnh hưởng của nhiễu và không áp dụng được cho quá trình dao động hoặc quá trình không ổn định (trừ khi chỉ chứa một khâu tích 24

phân bậc nhất). Với các quá trình có tính phi tuyến mạnh các số liệu đặc tính nhận được phụ thuộc rất nhiều vào biên độ và chiều thay đổi giá trị đặt. Phương pháp kẻ tiếp tuyến đồ thị dựa trên sự đánh giá trực quan nên thường cho số liệu L và a kém chính xác. Điểm đáng chú ý nữa là đặc tính đáp ứng của hệ kín với giá trị đặt thường hơi quá dao động (với hệ số tắt dần 0.25).

Hình 2.3. Phương pháp vùng cho việc xác định L và T Ví dụ: Sử dụng Ziegler-Nichols 1 để xác định các thông số của bộ điều khiển PID cho đối tượng có hàm truyền đạt sau:

Chúng ta thực hiện trên Matlab như sau: --------------------------------------------------------------------------------------------------s=tf('s'); G=10/(s+1)/(s+2)/(s+3)/(s+4) step(G); k=dcgain(G) --------------------------------------------------------------------------------------------------Chúng ta nhận được:

25

Step Response 0.45 0.4 0.35

Amplitude

0.3 0.25 0.2 0.15 0.1 0.05 0

0

1

2

3

4

5

6

7

Time (sec)

b. Phương pháp Ziegler-Nichols thứ 2 Phương pháp Ziegler-Nichols thứ hai hay còn gọi là phương pháp ZieglerNichols vòng kín. Phương pháp này cho phép việc xác định các tham số của bộ điều khiển PID dựa trên hệ số khuyếch đại tới hạn và chu kỳ tới hạn của hệ kín được xác định qua thực nghiệm. Hệ số khuyếch đại tới hạn của quá trình là giá trị khuyếch đại mà một bộ điều khiển P đưa hệ kín tới trạng thái dao động xác lập. Các thông số của bộ điều khiển được xác đinh theo bảng 3, những tham số này dùng cho một hệ kín với ít sự giảm dao động. Phương pháp Ziegler-Nichols 2 được thực hiện như sau: - Thay bộ điều khiển PID trong hệ kín hình 2.4a bằng bộ khuyếch đại tỷ lệ. Sau đó tăng hệ số khuyếch đại tới giá trị tới hạn Ku để hệ kín ở biên giới ổn định, đáp ứng đầu ra y(t) có dạng dao động điều hòa. Từ đó xác định chu kì dao động Tu. - Từ hai thông số Ku và Tu các thông số bộ điều khiển P, PI hay PID được xác định theo các công thức trong bảng 3. 26

w

e

ku

Đối tượng

u

-

y

(a)

(b) Hình 2.4. Phương pháp xác định các thông số PID bằng phương pháp Ziegler-Nichols 2 Bảng 2.3. Các công thức xác định theo phương pháp Ziegler-Nichols 2 Bộ điều khiển

Kc

Ti

Td

P

0.5Ku

.....

….

PI

0.4Ku

0.8Tu

….

PID

0.6Ku

0.5Tu

0.125Tu

Hai phương pháp của Ziegler-Nichols đều có những ưu nhược điểm nhất định: phương pháp thứ hai dễ thực hiện hơn trong các bộ tự chỉnh, đơn giản hơn vì không cần dùng mô hình ước lượng. Tuy nhiên nhược điểm của phương pháp thứ hai là quá trình thử nghiệm đặc tính dao động phải lặp đi lặp lại nhiều lần, tốn thời gian và dễ gây mất ổn định cho hệ thống. Phương pháp chỉnh định này

27

được đưa ra theo kinh nghiệm để đạt hệ số tắt dần khoảng 25% nên đáp ứng hệ kín sẽ hơi quá dao động và hơi kém bền vững với sai lệch mô hình. Phương pháp thứ nhất không xác định được rõ ràng cách đo các khoảng thời gian – do dựa vào trực quan đồ thị và tiếp tuyến. Kết quả xác định được bằng phương pháp thứ nhất trong nhiều trường hợp thường cao hơn 25%-40% so với phương pháp thứ hai, dẫn tới đáp ứng dao động nhiều hơn. Chất lượng của các bộ điều khiển được chỉnh định bằng phương pháp Ziegler-Nichols phụ thuộc rất nhiều vào giá trị thời gian trễ của đối tượng điều khiển, do đó thường có độ suy giảm kém và độ quá điều chỉnh khi đáp ứng với thay đổi điểm đặt, thay đổi tải là quá lớn. Phương pháp Ziegler-Nichols cho kết quả tốt hơn khi dùng chỉnh định với bộ điều khiển PID và kết quả kém hơn khi dùng với bộ điều khiển PI. 2.2.3. Phương pháp Cohen-Coon Phần lớn các đối tượng điều khiển trong công nghiệp có thể được xấp xỉ hóa bằng mô hình bậc nhất có trễ (FOPDT - First Order Plus Dead Time). Dựa trên đặc tính đó phương pháp chỉnh định Cohen-Coon được xây dựng cho một lớp đối tượng có hàm truyền đạt: Gp (s)  G1(s) 

K1e Tt 1s 1 Ts 1

(2.8)

Tiêu chuẩn cần đạt được chính là khử bỏ nhiễu phụ tải. Chỉ tiêu thiết kế chính là tỉ lệ suy giảm phần tư biên độ khi đáp ứng với nhiễu. Các chỉ tiêu thiết kế là cực đại hóa hệ số, cực tiểu hóa sai số trạng thái xác lập và tỉ lệ suy giảm phần tư biên độ cho các bộ điều khiển P và PD. Các tham số của bộ điều khiển PI được tính bằng cách cực tiểu hóa tiêu chuẩn tích phân sai lệch (IE - Integral Error) và yêu cầu đáp ứng tỉ lệ suy giảm phần tư biên độ. 

I IE   e(t )dt

(2.9)

0

Các tham số bộ điều khiển PID được tính toán với cùng chỉ tiêu như bộ điều khiển PI. Mối liên hệ giữa các tham số bộ điều khiển trong bảng 4 được cho 28

dưới dạng các tham số và được tính toán từ các tham số của mô hình FOPDT (2.8): 

K1Tt1 T1



Tt1 Tt1  T1

(2.10)

Theo nghiên cứu của K.J.Åstrõm và T.Hãgglund, phương pháp chỉnh định Cohen-Coon có thể áp dụng để chỉnh định bộ điều khiển PID on-line nếu đã biết các tham số của mô hình đối tượng điều khiển FOPDT (2.8)

Bảng 4: Luật chỉnh định PID theo phương pháp Cohen-Coon Bộ điều khiển

KP

TI

TD

P

1  0.35  1     1 

.....

.....

PI

0.9  0.92  1     1 

3.3  3 Tt1 1  1.2

.....

PD

1.24  0.13  1     1 

PID

1.35  0.18  1     1 

0.27  0.36 Tt1 1  0.87 2.5  2 Tt1 1  0.39

0.37  0.37 Tt1 1  0.81

Nếu thời gian trễ quá trình là nhỏ (trong giới hạn gần cận không) thì việc tăng lượng lớn các hệ số khuyếch đại sẽ là điều được báo trước. Do đó phương pháp này không thích hợp cho các hệ thống mà không có hoặc gần như không có thời gian trễ. 2.2.4. Phương pháp phản hồi rơle (Relay feedback) Hạn chế chính của phương pháp đáp ứng quá độ là rất nhạy cảm với nhiễu bởi vì dựa trên các thí nghiệm vòng hở. Các phương pháp dựa trên phản hồi rơle tránh được khó khăn này bởi vì nó yêu cầu các thí nghiệm được thực hiện trong vòng kín. 29

Tư tưởng cơ bản là sự nhận thấy nhiều quá trình có chu kì dao động giới hạn dưới tác dụng phản hồi rơle. Sơ đồ khối của hệ thống hệ thống như vậy được chỉ ra ở hình dưới: uc



+d

u -d

G(s)

y

-1 Hình 2.5. Sơ đồ khối của hệ thống phản rơle Phương pháp phản hồi rơle thực chất là một phương pháp nhận dạng trực tiếp trong vòng kín. Phương pháp này được K.J.Åstrõm và T.Hãgglund đưa ra dựa trên sự cải tiến phương pháp Ziegler-Nichols thứ hai. Điểm khác biệt là việc thay bộ điều khiển P bằng một khâu rơ le (khâu hai vị trí). Hệ tuyến tính với điều khiển rơ le. Nếu quá trình có đặc tính dao động tới hạn, hệ kín sẽ tự động đi tới trạng thái dao động tới hạn. Các tín hiệu đầu vào và đầu ra đạt được khi tín hiệu đặt uc là không và được chỉ ra ở đồ thị dưới:

Hình 2.6. Đặc tính vào ra của hệ phản hồi rơle Đồ thị trên cho thấy một chu kỳ dao động giới hạn được xác lập khá nhanh chóng. Ta có thể hiểu một cách trực quan hiện tượng xảy ra theo cách sau: đầu vào quá trình là một sóng vuông với tần số u. Bằng cách khai triển một chuỗi Fourier ta có thể thay thế đầu vào uc bằng tổng của các hàm sin với tần số u, 3u, ... Đầu ra được các hàm sin xấp xỉ hóa, có nghĩa là quá trình làm suy giảm ảnh hưởng của các hài bậc cao. Nếu coi biên độ của sóng vuông là d, 30

thành phần cơ bản có biên độ là 4d/. Tạo sự xấp xỉ, tất cả các hài bậc cao hơn có thể bỏ qua, ta tìm đầu ra quá trình là một hàm sin với tần số u và biên độ: a

4d



| G (iu ) |

(2.11)

Để có một dao động, đầu ra cũng phải đi qua không khi rơ le chuyển mạch. Hơn nữa, thành phần cơ bản của đầu vào và đầu ra phải ngược pha. Do đó ta có thể kết luận rằng tần số u phải giống như quá trình có sự trễ pha 180o. Các điều kiện cho sự dao động gồm có: argG(iu )   và | G (iu ) |

a 1  4d K u

(2.12)

trong đó Ku có thể được xem như tương đương với hệ số khuyếch đại của rơle cho quá trình truyền tín hiệu hình sin với biên độ a. Vì các lý do lịch sử, tham số này được gọi là hệ số khuyếch đại tới hạn. Đó là hệ số khuyếch đại đem một hệ thống với hàm truyền G(s) tới vùng biên giới ổn định dưới tác dụng điều khiển chỉ với khâu tỉ lệ. Chu kì Tu = 2/u cũng tương tự được gọi là chu kì tới hạn. Một thí nghiệm với phản hồi rơle do đó là một cách thích hợp để xác định chu kì tới hạn và hệ số khuyếch đại tới hạn. Chú ý rằng một tín hiệu đầu vào mà năng lượng của nó được tập trung tại tần số u là được sinh ra một cách tự động trong thí nghiệm.

PID A



u

T

Process

y

Relay -1

Hình 2.7. Sơ đồ khối của một bộ tự chỉnh rơle Khi việc điều chỉnh được yêu cầu, chuyển mạch được đặt tới T, có nghĩa là phản hồi rơle được kích hoạt và bộ điều chỉnh PID được ngắt ra. Khi chu kỳ 31

ổn định giới hạn được xác lập, các tham số PID được tính và bộ điều khiển PID sau đó được kết nối lại với đối tượng điều khiển. Toàn bộ quá trình thực hiện này được minh họa trên hình 2.8 y(t)

t u(t)

t

Hình 2.8. Kết quả minh họa việc sử dụng bộ tự chỉnh để xác định các thông số cho bộ điều khiển PID. Tuy nhiên, phương pháp này sẽ không làm việc cho tất cả các hệ thống. Thứ nhất, sẽ không có các dao động chu kỳ giới hạn duy nhất cho một hàm truyền tùy ý. Thứ hai, điều khiển PID không thích hợp cho tất cả các quá trình. Tự chỉnh rơle được tìm thấy một cách kinh nghiệm để làm việc tốt cho một phạm vi lớn các hệ thống trong điều khiển quá trình. Thuật toán thực hiện phương pháp rơle như sau: 1. Sử dụng đầu vào hằng không đổi uo, xác định đầu ra ổn định tương ứng yo. 2. Tiếp theo, thay đổi điểm đặt là yo. 3. Chuyển bộ điều khiển sang chế độ rơle và chọn một giá trị nhỏ h. Điều này sẽ thực hiện biểu thức sau:

4. Dưới chế độ này, đầu ra sẽ tiến đến một chu kỳ tới hạn.

32

Chu kỳ giới hạn xảy ra từ đây

Hình 2.9. Minh họa việc xác định chu kỳ tới hạn 5. Từ hình vẽ, đo biên độ của chu kỳ đầu ra a. Sử dụng cả a và h, hệ số khuyếch đại tới hạn có thể được xác định như sau: Ku = (4/π)(h/a) 6. Sử dụng Ku và chu kỳ tới hạn Pu (chu kỳ tới hạn được xác định như trên hình vẽ), áp dung bảng tính như phương pháp Ziegler-Nichols (Bảng 2.3) để xác định các thông số của bộ điều khiển PID. Phương pháp nhận dạng dựa trên phản hồi rơle đã được khẳng định và áp dụng rộng rãi trong thực tế. Ưu điểm của nó là việc tiến hành thử nghiệm đơn giản và phù hợp với đa số quá trình công nghiệp. Nếu quá trình có đặc tính dao động tới hạn thì hệ kín sẽ tự động tiến tới trạng thái dao động, ngay cả tín hiệu vào u cũng được tự động tạo ra. Ngoài ra nhờ khả năng tự do lựa chọn biên độ cho khâu rơle nên việc thử nghiệm với hệ kín đảm bảo quá trình nằm trong phạm vi an toàn mà quy trình công nghệ cho phép. Hơn nữa việc thử nghiệm trong vòng kín kết hợp với dạng tín hiệu vào – ra đặc biệt dẫn tới ảnh hưởng của nhiễu gần như được loại bỏ hoàn toàn.

33

Tuy nhiên hạn chế của phương pháp phản hồi rơle là dựa trên phương pháp Ziegler-Nichols thứ hai, phương pháp xác định các thông số bộ điều khiển PID hoàn toàn dựa trên kinh nghiệm, mặt khác công thức xác định các tham số đặc tính dao động tới hạn ở đây chỉ là xấp xỉ nên trong nhiều trường hợp, chất lượng điều khiển có thể không cao lắm. Mặc dù vậy phương pháp này cũng là mở đường cho hàng loạt phương pháp cải tiến khác được đề xuất và áp dụng. § ¸ p øng bËc thang TÝn hiÖu ®Çu ra

1.5

1

ZIEGLER-NICHOLS 2 RELAYFEEDBACK

0

TÝn hiÖu ®iÒu khiÓn

ZIEGLER-NICHOLS 1

0.5

0

5

20

25

30

20 15 [sec] Time Thêi gian [gi©y]

25

30

15

10

1.4 1.2 1 0.8

ZIEGLER-NICHOLS 1 ZIEGLER-NICHOLS 2

0.6

RELAYFEEDBACK 0.4

0

5

10

HƯỚNG DẪN NGHIÊN CỨU

Để hiểu và phục vụ cho tiết bài tập, học viên thực hiện ở nhà các nội dung sau: 1. Tổng hợp bộ điều khiển PID cho hệ có hàm truyền đạt mạch hở sau:: W

1 ( s  1)( s  2)( s  4)

34

Với các thông số yêu cầu như sau:   20%;

tqd  2s

2. Tổng hợp bộ điều khiển PID cho hệ có hàm truyền đạt mạch hở sau: W

1 ( s  1)( s  2)( s  4)

Với các thông số yêu cầu như sau:   20%;

tqd  2s

3. Tổng hợp bộ điều khiển PID cho hệ có hàm truyền đạt mạch hở sau: W

1 s( s  3)( s  7)

Với các thông số yêu cầu như sau:   15%;

tqd  2.3s

4. Tổng hợp bộ điều khiển PID cho hệ có hàm truyền đạt mạch hở sau: W

100 s(s  5)( s  9)

Với các thông số yêu cầu như sau:   15%;

tqd  1.5s

5. Tổng hợp bộ điều khiển PID cho hệ có hàm truyền đạt mạch hở sau: W

10 ( s  3)( s  7)( s  15)

Với các thông số yêu cầu như sau:   25%;

Ngày

tqd  2.5s

tháng

năm 2018

NGƯỜI BIÊN SOẠN

Thiếu tá, ThS Hoàng Văn Tiến

35