Ước lượng khối lượng tải trọng và quan sát nhiễu của hệ thống cầu trục

Kỹ thuật điều khiển & Tự độnghóa N. T. Thành, , N. T. T. Hằng, “Ước lượng khối lượng tải trọng hệ thống cầu trục.” 98 ƯỚC LƯỢNG KHỐI LƯỢNG TẢI TRỌNG VÀ QUAN SÁT NHIỄU CỦA HỆ THỐNG CẦU TRỤC Nguyễn Trung Thành1*, Trần Ngọc Quý2, Nguyễn Thanh Tiên3, Nguyễn Thị Thu Hằng1 Tóm tắt: Cầu trục là một hệ thống thiếu cơ cấu chấp hành và có trễ, quá trình làm việc thường chịu ảnh hưởng của nhiễu như ma sát hay nhiễu gió. Hơn thế khối lượng của tải trọng không phải lúc nào cũng dễ dàng xác định được, đặc biệt đối với những loại tải lớn. Đây là một trong những nguyên nhân làm giảm hiệu quả chống lắc tải trọng cầu trục. Tuy nhiên, nghiên cứu này đã tìm ra hướng giải quyết. Để góp phần nâng cao hiệu suất điều khiển, nhóm nghiên cứu thiết kế bộ quan sát nhiễu ma sát và gió để có thể loại bỏ ảnh hưởng của chúng. Một phương pháp ước lượng khối lượng tải trọng cũng đã được giới thiệu để lấy số liệu cho bài toán điều khiển. Ngoài ra nghiên cứu có đề xuất mô hình chống lắc sử dụng kết hợp các bộ điều khiển PID để có được kết quả tốt. Từ khóa: Điều khiển cầu trục; Bộ quan sát trạng thái; Góc lắc tải trọng; nhiễu gió; Nhiễu ma sát. 1. ĐẶT VẤN ĐỀ Cầu trục là một loại máy nâng vận chuyển được dùng nhiều trong các nhà máy, bến cảng, công trình xây dựng. Đặc điểm của cầu trục là sử dụng phối hợp một số chuyển động để nâng và di chuyển tải trọng lớn đến các vị trí cần thiết. Do kết cấu cơ khí cồng kềnh, nâng hạ tải bằng dây mềm, lại thêm môi trường làm việc ngoài trời nên cầu trục thường chịu ảnh hưởng của các loại nhiễu, đặc biệt là nhiễu ma sát và gió. Nếu không được kiểm soát thì nhiễu sẽ ảnh hưởng rất lớn đến chất lượng điều khiển làm giảm hiệu suất của cầu trục, đôi khi còn gây mất an toàn cho người và thiết bị tại hiện trường làm việc. Để giải quyết vấn đề này, người ta có thể lắp đặt cảm biến vận tốc góc và ước lượng góc lắc bằng bộ quan sát [2, 3]. Tuy nhiên phương pháp này rất khó đảm bảo ước lượng chính xác do sai lệch tham số. Đã có những nghiên cứu trước đây [1, 4] thực hiện ước lượng nhiễu, bù nhiễu và bù trễ trong luật điều khiển để cải thiện chất lượng hệ thống. Tuy nhiên góc lắc vẫn còn dao động với biên độ khá lớn, thậm chí chưa dập tắt được hoàn toàn khi tải trọng đã về đến đích. Vì vậy nghiên cứu này tập trung vào vấn đề quan sát, ước lượng tham số nhiễu, từ đó làm tiền đề cho việc thiết kế bộ điều khiển dập dao động tải trọng hiệu quả đồng thời ước lượng khối lượng tải trọng thông qua dòng điện của động cơ truyền động trong quá trình nâng hạ tải. Kết quả mô phỏng đã chứng tỏ hiệu quả của giải pháp đề xuất. Phương pháp ước lượng tải trọng được trình bày trong mục 2 của bài báo, mục 3 đề xuất phương pháp ước lượng nhiễu gió và nhiễu ma sát cho cầu trục 2-D và mục 4 là kết luận. 2. ƯỚC LƯỢNG KHỐI LƯỢNG TẢI TRỌNG CẦU TRỤC Khối lượng của tải trọng có quan hệ trực tiếp với công suất động cơ nâng hạ tải. Nên căn cứ vào dòng điện, điện áp động cơ nâng hạ ta có thể dự đoán tương đối chính xác khối lượng tải. Bằng cách ước lượng khối lượng tải trọng trong quá trình nâng của cầu trục, các tham số để vận chuyển có thể được xác định. Hình 1 mô tả hoạt động nâng tải. Động cơ servo quay với vận tốc góc g (rad/s) với Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san TĐH, 04 - 2019 99 momen y (Nm). Giá trị này được chuyển đổi sang tốc độ nâng ( )y t (m/s) và tín hiệu tuyến tính vào cơ cấu nâng ref y F (N), dòng điện cấp cho động cơ nâng hạ i(t). Hình 1. Sơ đồ phân tích lực cơ cấu nâng hạ tải trọng. Quãng đường nâng được xác định là ( )y t (m) và giả thiết khối lượng hàng hóa chưa biết rõ là m(kg). Ở đây ta coi hệ thống là vận hành trong môi trường thực tế nên khối lượng m sẽ là tác nhân trong chuyển động cuốn cáp nâng hàng mà hệ thống có thể nhận được. Xét một loại cầu trục lớn vận chuyển container, khối lượng của pully không đổi và được giả thiết là rất nhỏ so với tải trọng nên có thể bỏ qua. Lực ma sát giữa pully và dây kéo không tính đến. Để thuận tiện cho việc biểu diễn, gia tốc trọng trường g được coi là một đầu vào, thì phương trình chuyển động được suy từ phương trình cân bằng lực tại vị trí pully như sau: ( ) ref y my t F mg  (1) Khi tính đến truyền động của động cơ chấp hành (với giả thiết là động cơ điện một chiều kích từ độc lập) thì: ref y r p r F Ki K i r   trong đó: r là tỷ số truyền của động cơ tới pulley, rp là bán kính hiệu dụng tang trống, K là hệ số momen và ( : ) r p K r r K , i là dòng điện cuộn dây phần ứng. suy ra ( ) 1r id K y t gdt m               (2) Phương trình (1) được sử dụng để tạo ra bộ ước lượng tải trọng. Ban đầu một bộ quan sát để lấy tốc độ ước lượng ˆ( )y t được tạo ra như công thức (3). Trường hợp khối lượng tải không biết trước, khối lượng tham chiếu mn (kg) biết trước: ( )ˆ ˆ( ) ( ( ) ( )) 1 ( ) r m n i td K y t K y t y t g tdt m                     (3) Kỹ thuật điều khiển & Tự độnghóa N. T. Thành, , N. T. T. Hằng, “Ước lượng khối lượng tải trọng hệ thống cầu trục.” 100 trong đó Km là hệ số hiệu chỉnh tốc độ hội tụ của bộ quan sát, ( )y t là gia tốc nâng tải, ˆ( )y t là gia tốc nâng tải quan sát. Toán tử hóa các phương trình (2), (3) sau đó thay (2) vào (3) rút ra ˆ( )Y s : 2 1ˆ( ) . ( )r n m m n K ms m K Y s I s g s K mm s s      (4) Gọi sai lệch giữa tốc độ nâng ước lượng và tốc độ nâng thực tế là nhiễu dm thì: ˆ( ) ( ( ) ( )) m m d s K Y s Y s   (5) Thay (4) và (2) (sau khi đã toán tử hóa) vào biểu thức (5) suy ra: ( ) . ( ) ( ). ( ) 1 1 n nr m m n n m m m m mK d s I s D s I s m m m ms K      (6) Lượng dm xuất hiện là do sự chênh lệch giữa khối lượng tham chiếu mn và khối lượng thực m. Khối lượng thực m có thể được suy ra từ lượng chênh lệch này. Thành phần ( )mD s là trễ bậc 1 ở hệ số khuếch đại Kr nên sẽ hội tụ tại Kr, và do đó m d cuối cùng hội tụ theo (7). ( ) ( ) . n m r n m m d t K i t m m   (7) Từ (7) suy ra khối lượng ước lượng mˆ theo phương trình (8). Mẫu số biểu thức (8) được giả thiết là khác 0 trong hệ thống thực và hệ thống thí nghiệm được xét với tham số bảng 1, trong đó các tham số bộ điều khiển vị trí và bộ điều khiển vận tốc được xác định theo phương pháp Ziegler–Nichols 2 [6] cho đối tượng không cần biết rõ mô tả toán học. suy ra ( ) ˆ ( ) ( ) . ( ) r n r n m K m i t m t K i t m d t   (8) Bảng 1. Tham số cơ cấu nâng. Mô tả Ký hiệu Giá trị Hệ số tỷ lệ bộ điều khiển vị trí Krp 35 Hệ số tỷ lệ bộ điều khiển vận tốc Klp 27 Hệ số tích phân bộ điều khiển vận tốc Kli 15 Hệ số quan sát Km 100 Tỷ số truyền r 1/20 Bán kính tang trống rp 0.03 (m) Hệ số mô men K 1.8 Khối lượng tham chiếu mn 0.8 (kg) Khối lượng thực m 0.87 (kg) Phương trình (8) khẳng định rằng khối lượng tải trọng có thể được ước lượng từ hoạt động quan sát với khối lượng tham chiếu. Hệ thống điều khiển nâng hàng bao gồm bộ ước lượng khối lượng tải trọng có cấu trúc như hình 2. Bộ ước lượng tải trọng được thiết kế sử dụng ( )i t và y thu được từ hệ thống điều khiển nâng tải. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san TĐH, 04 - 2019 101 Các tham số được đặt như trong bảng 1, trong đó các bộ điều khiển PID được thiết kế theo phương pháp Ziegler–Nichols 2 [6] cho đối tượng không cần biết rõ hàm truyền đạt nhằm điều khiển vị trí và giảm dao động tải trọng cầu trục. Các thông số cầu trục được lấy từ một mô hình sử dụng cho phòng thí nghiệm dựa trên thiết bị thực được xây dựng với tỷ lệ mô hình là 1/500. Hình 3 biểu diễn kết quả của mô phỏng trên phần mềm Matlab và Simulink. Thời gian nâng được đặt tới 2 giây và tham số điều khiển Km trong phương trình (8) là để điều chỉnh thời gian hội tụ ước lượng được đặt tới 100. Hằng số thời gian là đủ nhỏ so với thời gian nâng hàng. Quan sát kết quả mô phỏng ta thấy 1 giây trước khi kết thúc hoạt động nâng tải thì khối lượng ước lượng là ˆ 0.8705m kg , 1 giây sau khi kết thúc hoạt động nâng tải thì khối lượng ước lượng là ˆ 0.8695m kg so với giá trị thực tế 0.87m kg và ổn định hơn so với trước khi kết thúc nâng tải. Do đó ước lượng của khối lượng thu được khi kết thúc quá trình nâng sai khác không đáng kể so với giá trị thực có thể được sử dụng làm dữ liệu cho việc thiết kế điều khiển sau này. Mô phỏng với các tải trọng có số liệu khác nhau ta thu được các thông số ước lượng khá gần với giá trị thực. Các số liệu so sánh được thống kê trong bảng 2. Bảng 2. Kết quả tham số ước lượng với các tải trọng khác nhau. m(kg) 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 mˆ(kg) 0.2001 0.4005 0.6010 0.7985 0.9980 1.1980 y y efry efry  yˆ   md mˆ 1 lp liK K s  rpK mK i Hình 2. Sơ đồ hệ thống điều khiển động cơ nâng hạ tải. Hình 3. Kết quả mô phỏng khối lượng ước lượng, vị trí tải trọng, dòng điện động cơ truyền động nâng tải. 0 1 2 3 4 5 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 thoigian(s) v it ri (m ) yref y 0 1 2 3 4 5 0 0.1 0.2 0.3 0.4 thoigian(s) d o n g d ie n (A ) 0 1 2 3 4 5 0.8 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.9 thoigian(s) k h o ilu o n g u o c lu o n g (k g ) 0.8705 kg 0.8695 kg Kỹ thuật điều khiển & Tự độnghóa N. T. Thành, , N. T. T. Hằng, “Ước lượng khối lượng tải trọng hệ thống cầu trục.” 102 3. XÂY DỰNG BỘ QUAN SÁT NHIỄU VÀ MÔ HÌNH CHỐNG LẮC Hình 4. Mô hình vật lý cầu trục 2-D. Hình 4 mô tả loại cầu trục 2D (dạng cổng trục). Động cơ di chuyển với momen m  và tốc độ góc m  và dịch chuyển dọc theo dầm cầu. Gọi efr F là lực tác động làm dịch chuyển xe với tốc độ x trong suốt quá trình chuyển động tuyến tính, x là khoảng cách dịch chuyển. Gọi M là khối lượng xe con, L là chiều dài dây cáp sau khi đã nâng tải,  là góc lắc tải so với phương thẳng đứng, d F là nhiễu ma sát giữa đường ray và xe con, w F là nhiễu gió. Xét phương trình chuyển động Lagrange: efrac ac d d L L F F dt x x          (9) w ac ac d L L F dt           (10) Sử dụng phương pháp mô tả toán học theo tài liệu [5] thu được hệ tuyến tính sau: ef w 2 0 0 1 1 ( ) ( ) ( ) ( ) ( ) 0 0 1 ( ) 1 ( ) ( ) ( ) ( ) 0 0 p t r t t d t p t p t t p t m g mx t x t m m L F t F td t t m m F tdt t m m g t m L m m L m L                                                                       (11) Trường hợp tổng quát, ảnh hưởng của nhiễu gió tác động vào tải thuộc không gian đa chiều, nhưng hoạt động chống lắc của cầu trục 2-D chỉ thực hiện theo một chiều dọc dầm cầu. Do vậy trong nghiên cứu này, tác giả giới hạn ảnh hưởng của nhiễu gió chỉ theo một chiều trên mặt phẳng chuyển động. Hình 5 biểu diễn sơ đồ khối nguyên lý làm việc và chống lắc của cầu trục. Vị trí đặt efrx được đặt theo quỹ đạo tối ưu thiết kế dạng hình thang. Bằng việc căn chỉnh thời gian tăng tốc, giảm tốc theo chu kỳ lắc của tải trọng, góc lắc có thể được giữ ở cuối hành trình di chuyển. Hệ điều khiển chống lắc được xây dựng để tạo ra một góc lắc đặt efr dựa trên vị trí đặt efrx . Vì thực tế điều khiển thấy có trễ nên giá trị thu được từ cảm biến không được sử dụng trực tiếp để điều khiển. Ở nghiên cứu này tác giả chọn một giải pháp là lấy góc lắc ước lượng từ giá trị của bộ quan sát. x 0 mp mt Fref x(t) (t) L Fd Fw Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san TĐH, 04 - 2019 103 x efrx 1 vp viK K s  efr efrF 1 pp piK K s  efrx x Hình 5. Hệ thống điều khiển cầu trục hồi tiếp cơ bản. Phù hợp với công thức (7) sau khi chuyển đổi sang mô hình chuẩn quan sát được thu được một bộ quan sát với góc lắc ước lượng ˆ từ tốc độ xe con x và lực tác động vào cầu trục F bao gồm cả nhiễu, đã được thực hiện theo phương pháp của Gopinath [7] theo công thức (12, 13). ˆ ( ) q v t là giá trị trạng thái quan sát của bộ quan sát, ˆ ˆˆ ˆ ˆ, , , , q q q q q A B J C D là các ma trận hệ số sử dụng trong tính toán quan sát. ˆ ˆ ˆ ˆˆ( ) ( ) ( ) ( ) q q q q q t t x t t  v A v B J F  (12) ˆ ˆˆ ˆ( ) ( ) ( ) q q q q t t x t z C v D  (13) trong đó: 1 2 ˆ ˆ ˆ( ) ( ) ( ) T q q q t t t   v v v , ˆˆˆ ( ) ( ) ( ) T q t t t         z  , ef w( ) ( ) ( ) ( ) T r dt F t F t F t      F Góc lắc quan sát có thể được ước lượng kém chính xác bởi sai lệch tham số do trễ. Để đối phó với vấn đề này, cần phải xét góc lắc thu được x  với thời gian trễ d T . Cũng vì xét trong chu kỳ điều khiển với bộ quan sát ước lượng trạng thái nên góc lắc được ước lượng cũng theo chu kỳ, do đó tác giả đề xuất “bộ quan sát điều chỉnh thời gian” trong cùng một thời gian trễ d T . Bằng việc sử dụng điều chỉnh trạng thái góc lắc với thời gian trễ có thể thu được một góc lắc ước lượng chính xác hơn. Để giải quyết bài toán điều khiển (12) trước hết cần ước lượng nhiễu gió và nhiễu ma sát. Trong một số nghiên cứu trước đây, các bộ quan sát nhiễu gió và nhiễu ma sát thường được xét riêng hoặc chỉ xét một loại [1], nên độ chính xác ước lượng khó được đảm bảo do thực tế nhiều khi chúng cùng tác động. Nếu kể đến tất cả các tín hiệu vào thì đầu ra của bộ quan sát nhiễu có thể được biểu diễn như sau: W ˆ ( ) ( ) ( ) ( ) ( ) ( ) d d L F s H s F s H s H s F s     (14) W W ˆ ( ) ( ) ( ) ( ) ( ) d F H s H s F s H s F s      (15) trong đó: 2 2 2 ( ) ( 2 ) n n n H s s s        ; 2 2 ( )L s H s Ls g   ; 2 1 1 02 1 02 2 2 2 2 1 0 41 2 ( ) n n s s mL mL H s L s s                         (16) Kỹ thuật điều khiển & Tự độnghóa N. T. Thành, , N. T. T. Hằng, “Ước lượng khối lượng tải trọng hệ thống cầu trục.” 104 1 0 , ,H    và n  có thể được chọ tùy ý nhưng cần thỏa mãn 1 1 / n   và 2 1 1 0 4 n      . Kết quả là mô hình đơn giản theo (16) có thể được sử dụng bởi sự xấp xỉ: ( )H s L  (17) Ở đây ước lượng nhiễu không bao gồm nhiễu ma sát được thiết kế là ( ) m F s và được tìm ra theo (18). Công thức (18) thu được bằng ước lượng thành phần nhiễu ma sát ( ) d F s có trong công thức (14). Ảnh hưởng của nhiễu ma sát là không bao gồm trong ( ) m F s như thể hiện bởi phương trình này, và do đó một ngoại lực ước lượng chỉ chịu ảnh hưởng của nhiễu gió W ( )F s thu được: W W2 1 1ˆ ˆ ˆ( ) ( ) ( ) ( ) ( ) m d g F s F s F s H s F s L L Ls g      (18) ( ) m F s trong phương trình (18) có thể được tính toán từ ước lượng ˆ ( ) d F s và W ˆ ( )F s của hai bộ quan sát nhiễu khi nhiễu gió và nhiễu ma sát cùng ảnh hưởng. Tuy nhiên để tìm thấy thành phần nhiễu gió W ( )F s trực tiếp từ ( ) m F s trong một hàm truyền thì rất khó thực hiện. Mặt khác, bởi vì điểm cực dao động /s j g L được bao gồm trong mẫu số của hàm truyền đạt từ công thức (18), sự ước lượng được thực hiện với đặc tính dao động có thêm nhiễu gió là khá khó. Do vậy, cần thêm vào bộ lọc ( ) f G s theo công thức (19), các cực dao động có mặt trong mẫu số của (18) có thể được thực hiện loại bỏ điểm cực, điểm không và giá trị cố định của nhiễu ước lượng (giá trị cuối cùng của nhiễu step) có thể thu được ổn định: 2 2 2 2 ( ) 2 n f n n L Ls g G s g s s        , (19) ' 2 W ˆ( ) ( ). ( ) ( )m f mF s G s F s H s F  (20) Do vậy thành phần nhiễu ' ( ) m F s được cho bởi (28) tỷ lệ với nhiễu gió cố định. Kết quả là giá trị ước lượng ' W W 1ˆ ˆ ˆ( ) ( ) ( ) ( )f dF s G s F s F s L       (21) Một khâu trễ được đưa thêm vào ước lượng nhiễu, nhưng vì trong nghiên cứu này nhiễu gió được giả thiết có dạng hàm step, nên sự hội tụ được giả sử là xảy ra ở cuối quá trình. Hơn nữa, mặc dù nhiễu gây bởi thời gian lấy mẫu để phát hiện góc lắc có được xét, độ chính xác của ước lượng góc lắc có thể được cải thiện bằng cách coi nhiễu của gió tác động một chiều. Hình 6 biểu diễn một sơ đồ khối của quá trình ước lượng nhiễu và hình 7, hình 8 là các kết quả mô phỏng. Các kết quả ước lượng khá chính xác trong cả hai trường hợp trước và sau khi có sự tham gia của cả hai loại nhiễu. Nghiên cứu khoa học công nghệ Tạp chí Nghiên cứu KH&CN quân sự, Số Đặc san TĐH, 04 - 2019 105 x xefrx ' wFˆ 1 vp viK K s  1 L fG wFˆ dFˆ ef d ˆrF F w w 1 p iK K s  efr ˆ  a efrF tF dF wF 1 pp piK K s  efrx Hình 6. Sơ đồ hệ thống ước lượng nhiễu và điều khiển chống lắc. (a) (b) Hình 7. Kết quả ước lượng nhiễu gió trước (a) và sau (b) khi kết hợp nhiễu ma sát. (c) (d) Hình 8. Kết quả ước lượng nhiễu ma sát trước (c) và sau (d) khi kết hợp nhiễu gió. 4. KẾT LUẬN Ở nghiên cứu này nhóm tác giả tập trung giải quyết được hai vấn đề. Một là ước lượng được khối lượng tải trọng nâng hạ cầu trục trực tiếp từ thông số hệ truyền động, cụ thể là từ dòng điện phần ứng động cơ nâng hạ tải trọng. Hai là xây dựng được bộ quan sát nhiễu là ma sát và gió cho hệ cầu trục 2-D, đây là các nhiễu chính tác động vào hệ thống mà trước đây vốn chỉ được nghiên cứu độc lập. Các kết quả thu được từ thuật toán và mô phỏng đủ tin cậy để có thể sử dụng vào việc giải quyết tiếp bài toán điều khiển chống lắc tải trọng một cách thuận lợi hơn. TÀI LIỆU THAM KHẢO [1]. Ngô Quang Hiếu, “Điều khiển chống lắc hệ cần cẩu container có bù ma sát”, Tạp chí Khoa học trường Đại học cần thơ (2013); tr. 8-14 Kỹ thuật điều khiển & Tự độnghóa N. T. Thành, , N. T. T. Hằng, “Ước lượng khối lượng tải trọng hệ thống cầu trục.” 106 [2]. Kim YS et al, “Anti-sway control of container cranes: Inclinometer, observer, and state feedback”, Int J Control Automat Syst (2004);2(4): pp.435–449. [3]. Honda F. et al, “Sway angle measurement in an overhead crane using an observer with an angular speed sensor”, Proc JSME (2006);72-C(722). [4]. Kaneko K, Mine H, Ohishi K. “Crane anti-sway control using a sway angle observer to correct sensor delay”, Trans IEICE (2009);129-D(6): pp.555– 562. [5]. Nguyen Trung Thanh et al, “Dynamic modeling and analysis of a three – dimentional overhead crane system with the variation of load mass and hoisting/lowering force”, Tạp chí KH&CN, Đại học Thái Nguyên số 3(2017). [6]. Nguyễn Doãn Phước, “Lý thuyết điều khiển tuyến tính”, NXB Khoa học và Kỹ thuật, (2009); tr.172-174 [7]. Mita T. et al, “Fundamentals of digital control”, Corona Press; 1988. ABSTRACT ESTIMATING THE MASS AND ANTI-SWAY OF PAYLOAD OF CRANE SYSTEM Crane is a un-actuator and delay system, the working process is often affected by disturbances such as friction or wind. Moreover, the payload is not always easy to determine, especially for super loads. This is one of the reasons to reduce the effectiveness of load anti-sway. However, this study has found a solution. To improve control performance, we designed a combined observer of friction and wind disturbance to cancel their effects. A method of estimating the payload is used to obtain data for the control problem. In addition, this study used a combination model including of PID controllers to obtain better results. Keywords: Crane control; States observer; Payload sway; Wind disturbance; Friction disturbance. Nhận bài ngày 09 tháng 01 năm 2019 Hoàn thiện ngày 08 tháng 3 năm 2019 Chấp nhận đăng ngày 15 tháng 3 năm 2019 Địa chỉ: 1 Trường Đại học SPKT Hưng Yên; 2 Viện KH-CNQS; 3 Học viện KTQS. * Email: [email protected].