Lê Thị Huyền Linh và Đtg<br />
<br />
Tạp chí KHOA HỌC & CÔNG NGHỆ<br />
<br />
86(10): 195 - 200<br />
<br />
MỘT PHƢƠNG PHÁP XÂY DỰNG ĐIỀU KHIỂN DỰ BÁO<br />
DỰA TRÊN MÔ HÌNH GAUSS<br />
Lê Thị Huyền Linh*, Nguyễn Thị Mai Hƣơng<br />
Trường Đại học Kỹ thuật Công nghiệp – ĐH Thái Nguyên<br />
<br />
TÓM TẮT<br />
Mô hình Gauss cung cấp một mô hình xác xuất không tham số cho việc nhận dạng mô hình hộp<br />
đen phi tuyến. Mô hình Gauss có thể xác định các khu vực không gian đầu vào có chất lượng dự<br />
báo kém (do thiếu dữ liệu hoặc quá phức tạp) bằng cách chỉ ra những phương sai cao hơn lân cận<br />
giá trị trung bình dự đoán. Khi sử dụng mô hì nh Gauss chỉ cần rất ít tham số để đạt được sự tối ưu .<br />
Bài báo này giới thiệu khả năng áp dụng của mô hình Gauss trong điều khiển dự báo . Phương pháp<br />
này được ứng dụng trong việc điều khiển quá trình trung hòa độ pH .<br />
Từ khóa: Mô hình Gauss, Mô hình điều khiển dự báo, Mô hình xác xuất không tham số.<br />
<br />
<br />
ĐẶT VẤN ĐỀ<br />
Model Predictive Control (MPC) là thuật toán<br />
điều khiển sử dụng mô hì nh tường minh để<br />
dự đoán đối tượng đáp ứng tương lai . Thuật<br />
toán này sẽ tối ưu các tham số điều khiển để<br />
có được đối tượng đáp ứng tương lai tối ưu .<br />
MPC được sử dụng trong bài toán xử lý<br />
những ràng buộc , tính ổn định thường gặp<br />
trong thực tế điều khiển mà các phương pháp<br />
điều khiển tối ưu thông thường không đáp<br />
ứng được. Trong điều khiển dự báo, chất<br />
lượng điều khiển phụ thuộc vào mô hình dự<br />
báo vì vậy một trong những phương pháp<br />
nâng cao chất lượng của hệ điều khiển dự báo<br />
là phương pháp xây dựng mô hình điều khiển.<br />
Nguyên tắc điều khiển mô hình dự báo phi<br />
tuyến (NMPC) với mô hì nh Gaus s thường<br />
được sử dụng để mô hì nh hóa các thông số<br />
phi tuyến tĩnh [1], [2]. Trong nghiên cứu này ,<br />
mô hình Gauss được phát triển để ứng dụng<br />
cho việc mô hình hóa hệ thống động.<br />
MÔ HÌNH HÓA HỆ THỐNG ĐỘNG VỚI<br />
GAUSS<br />
Mô hình Gauss là một mô hình xá<br />
c suất<br />
không tham số linh hoạt với những dự báo<br />
không xác định . Mô hình này là tập hợp các<br />
biến ngẫu nhiên được phân bố theo đường<br />
<br />
<br />
Tel: 0982 847826, Email: lethihuyenlinh@gmail.com<br />
<br />
Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br />
<br />
cong Gauss . Giả thiết hàm y = f(x) là mối<br />
quan hệ giữa đầu vào x và đầu ra y. Trong đó:<br />
y1,..yn ~ N(0,∑), ∑pq= Cov(yp, yq) = C(xp, xq)<br />
là hiệp phương sai giữa giá trị đầu ra và đầu<br />
vào tương ứng x p, xq. Do đó giá trị trung bình<br />
μ(x) ( luôn giả đị nh bằng<br />
0) và hàm hiệp<br />
phương sai C(xp, xq) sẽ xác đị nh mô hình<br />
Gauss một cách đầy đủ . Chú ý rằng hàm hiệp<br />
phương sai C(xp, xq) có thể là bất cứ hàm nào<br />
với đặc tí nh là nó có thể tạo ra một ma trận<br />
hiệp phương sai xác đị nh dương.<br />
Hàm hiệp phương sai thường được chọn là:<br />
C ( x p , xq ) = u1exp[-<br />
<br />
1 D<br />
å wd ( xdp - xqd )2 ]+u0<br />
2 d=1<br />
<br />
(1)<br />
<br />
T<br />
<br />
Trong đó : Q = [w1...wD u0 u1 ] là các tham<br />
số bậc cao của hàm hiệp phương sai và D là<br />
chiều đầu vào . Ở đây , các tham số có thể<br />
được xác đị nh dựa vào các dữ liệu đã có . Sau<br />
đó người ta có thể sử dụng tham số như là<br />
những chỉ số để xác định mức độ quan trọng<br />
của các đầu vào tương ứng: nếu d 0 hoặc<br />
gần bằng 0, điều đó có nghĩ a là đầu vào kí ch<br />
thước d chứa í t thông tin và có thể được bỏ<br />
qua. Hãy xét tập hợp các vector D chiều: đầu<br />
vào X = [x1, x2,. . ,xN] và đầu ra y =[y1, y2,. . .<br />
, yN]T. Dựa trên các dữ liệu (X, y), và một<br />
vector đầu vào mới x* đã cho, ta mong muốn<br />
tìm thấy những phân bố dự báo của đầu ra y *<br />
tương ứng. Với mô hình này, mục tiêu chính<br />
195<br />
<br />
http://www.lrc-tnu.edu.vn<br />
<br />
Lê Thị Huyền Linh và Đtg<br />
<br />
Tạp chí KHOA HỌC & CÔNG NGHỆ<br />
<br />
là điều chỉnh các tham số của hàm hiệp<br />
phương sai để đạt được sự tối ưu. Điều này<br />
được thực hiện bằng cách phát huy tối đa biểu<br />
đồ hợp lệ của tham số , yêu cầu của việc tính<br />
toán này chỉ là tương đối bởi vì nghịch đảo<br />
của ma trận hiệp phương sai dữ liệu (N x N)<br />
phải được tính theo mỗi bước lặp.<br />
Phương pháp này có thể dễ dàng sử dụng<br />
trong tí nh toán hồi quy. Trên cơ sở huấn luyện<br />
một lượng X, một ma trận hiệp phương sai K<br />
(N x N) được xác định. Như đã đề cập trước,<br />
mục đích là để tìm thấy sự phân bố của đầu ra<br />
y* tương ứng với một vector đầu vào x*:<br />
x* = [x1 (N + l), x2 (N + l) ,..., xD (N +1)]T<br />
Đối với một thử nghiệm đầu vào x*, việc phân<br />
bố dự đoán của đầu ra tương ứng y *|x*, (X,<br />
y) là đường cong Gauss, với giá trị trung bình<br />
và phương sai:<br />
( x* ) k ( x* )T K T y<br />
(2)<br />
<br />
2 ( x* ) k ( x* ) k ( x* )T K 1k ( x* ) v0<br />
<br />
(3)<br />
Trong đó: k x* [C x1 , x* ,. . . ,C x N , x* T<br />
là Nx1 vector của hiệp phương sai giữa<br />
trường hợp kiểm tra và huấn luyện, và<br />
k x* C x* , x* là hiệp phương sai giữa<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
các thử nghiệm đầu vào và chính nó.<br />
Trước khi dự đoán ta phải tính toán các giá trị<br />
bất định cho các dự báo trong tương lai. Điều<br />
đó sẽ cung cấp đầu vào cho việc ước lượng<br />
những giá trị trung bình và những giá trị bất<br />
định trong tương lai xa hơn.<br />
Nếu bây giờ chúng ta xét một đầu vào ngẫu<br />
nhiên mới x* ~ N ( x ,<br />
) và một khai<br />
x*<br />
<br />
<br />
<br />
triển Taylor là ( x ) và 2 ( x* ) lân cận x ,<br />
thì phân bố dự báo lại chính là mô hình Gauss<br />
với giá trị trung bình và phương sai.<br />
*<br />
<br />
86(10): 195 - 200<br />
<br />
m( x* ) Ex* [ ( x* )]<br />
k ( ( x* )T ) K 1 y<br />
<br />
(4)<br />
<br />
r ( x ) Ex* [ ( x )] varx* ( ( x ))<br />
*<br />
<br />
2<br />
<br />
*<br />
<br />
*<br />
<br />
2 ( ( x* ))<br />
1 2 2 ( x* )<br />
*<br />
*<br />
2 x*x*T x ( x )<br />
( x* )<br />
( x* )<br />
<br />
*T<br />
* ]}<br />
*<br />
*<br />
( x* ) x ( x ) ( x* ) x ( x )<br />
trace{ x*[<br />
<br />
(5)<br />
<br />
Phương trình (4) và (5) có thể được áp dụng<br />
cho tính toán trước dự báo với sự phân tán<br />
của các giá trị bất định. Mô hình Gauss giống<br />
như mạng Nơron được dùng cho các mô hình<br />
phi tuyến tĩnh và có thể được áp dụng cho mô<br />
hình hóa các hệ thống động, nếu đầu vào có trễ<br />
và tín hiệu đầu ra quay trở về đầu vào (tín hiệu<br />
hồi quy). Ưu điểm của mô hình Gauss là khi<br />
đạt được tối ưu xuất hiện ít tham số hơn việc<br />
sử dụng mô hình Nơron. Như vậy đầu ra hiện<br />
tại ở bước thứ k phụ thuộc vào đầu vào trước<br />
đó giống như việc điều khiển trước đầu vào.<br />
<br />
x(k ) [ yˆ (k -1), yˆ (k - 2),.... yˆ (k - L), u (k -1),<br />
u (k - 2),..., u (k - L)]T<br />
yˆ (k ) f ( x(k )) e<br />
<br />
(6)<br />
<br />
Ở đây x biểu thị vector trạng thái bao gồm các<br />
giá trị đầu ra y, giá trị đầu vào u có độ trễ L và<br />
e là nhiễu trắng.<br />
Dựa trên các mối quan hệ đã trình bày trên,<br />
mô hình được tạo ra sẽ mô tả các đặc tính<br />
động của hệ thống phi tuyến, đồng thời cung<br />
cấp thông tin về độ tin cậy của các dự báo đó<br />
ở cùng một thời điểm. Mô hình Gauss có thể<br />
xác định các khu vực không gian đầu vào có<br />
chất lượng dự báo kém (do thiếu dữ liệu hoặc<br />
quá phức tạp) bằng cách chỉ ra những phương<br />
sai cao hơn lân cận giá trị trung bình dự đoán.<br />
MÔ HÌNH ĐIỀU KHIỂN DỰ BÁO PHI TUYẾN<br />
<br />
Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br />
<br />
196<br />
<br />
http://www.lrc-tnu.edu.vn<br />
<br />
Lê Thị Huyền Linh và Đtg<br />
<br />
Tạp chí KHOA HỌC & CÔNG NGHỆ<br />
<br />
86(10): 195 - 200<br />
<br />
tín hiệu điều khiển tối ưu u (k + j|k), j = 0,. . .<br />
Nu-1 là được sử dụng.<br />
Trong lần lấy mẫu tiếp theo một mẫu đầu ra<br />
mới được xác định và toàn bộ quy trình sẽ<br />
được lặp lại. Nguyên lý này gọi là phạm vi<br />
chiến lược dịch dần (RHS). Hàm mục tiêu có<br />
dạng như sau:<br />
Hình 1. Sơ đồ khối mô hình hệ thống điều khiển<br />
dự báo<br />
<br />
Mô hình điều khiển dự báo phi tuyến được áp<br />
dụng với mô hình Gauss được mô tả bằng sơ<br />
đồ khối tổng quát như Hình 1. Mô hình được<br />
sử dụng là cố định, được nhận dạng off-line,<br />
có nghĩa là sử dụng thuật toán điều khiển<br />
không thích nghi. Vì vậy cấu trúc của toàn bộ<br />
vòng lặp điều khiển ít phức tạp như trong<br />
trường hợp mô hình biến thiên theo thời gian.<br />
Nội dung của phương pháp điều khiển dự báo<br />
được trình bày dưới đây:<br />
Dự báo tín hiệu đầu ra của hệ thống y(k+j)<br />
được tính cho mỗi lần lấy mẫu rời rạc k đối<br />
với một phạm vi rộng trong tương lai<br />
(j=N1,...,N2). Những dự báo này được biểu<br />
diễn là yˆ k j | k với j là bước dự đoán<br />
trước. N1 và N2 lần lượt là giới hạn dưới và<br />
trên của phạm vi dự báo. Giới hạn dưới và<br />
trên của miền tín hiệu đầu ra dự báo xác định<br />
trùng với phạm vi dự báo, trong đó có sự phối<br />
hợp giữa đầu ra và tín hiệu tham chiếu mong<br />
đợi. Đầu ra tín hiệu dự báo được tính toán từ<br />
mô hình đối tượng. Các dự báo này phụ thuộc<br />
vào khả năng điều khiển trong tương lai u (k<br />
+ j|k), j = 0,. . . , Nu-1, được áp dụng từ trước<br />
thời điểm k.<br />
Quỹ đạo tham chiếu được xác định r(k+j|k), j<br />
= N1,...,N2, xác định quá trình đáp ứng tham<br />
chiếu từ giá trị hiện tại y(k) tới quỹ đạo điểm<br />
đặt (k ) .<br />
Vector của tín hiệu điều khiển tương lai (u (k<br />
+ j|k), j = 0,. . . , Nu-1) được tính toán bằng<br />
cách tối thiểu hóa hàm mục tiêu sao cho sai<br />
lệch giữa đáp ứng tham chiếu r(k+j|k) và đáp<br />
ứng dự báo yˆ k j | k , j = N1….N2 là nhỏ<br />
nhất. Chỉ có giá trị đầu tiên u(k|k) của véctơ<br />
<br />
Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br />
<br />
min[r (k P) yˆ (k P)]2<br />
U (k )<br />
<br />
(7)<br />
<br />
Với điều kiện biên:<br />
var yˆ (k P) kv<br />
<br />
(8)<br />
<br />
U (k ) kih<br />
<br />
(9)<br />
<br />
U (k ) kir<br />
<br />
(10)<br />
<br />
x(k ) ksh<br />
<br />
(11)<br />
<br />
x (k ) ksr<br />
<br />
(12)<br />
<br />
Trong đó U(k)=[u(k) …u(k+P)] là tín hiệu<br />
đầu vào, P là điểm trùng hợp của giá trị đầu ra<br />
và giá trị tham chiếu, các bất đẳng thức từ (8)<br />
đến (12) biểu diễn sự ràng buộc của phương<br />
sai đầu ra kv, ràng buộc cứng của đầu vào kih,<br />
ràng buộc biến thiên của đầu vào kir, ràng<br />
buộc cứng của vectơ trạng thái ksh , ràng buộc<br />
biến thiên của vectơ trạng thái ksr. Mô hình<br />
như vậy được gọi là mô hình Gauss.<br />
Thuật toán tối ưu hóa thực chất là một<br />
chương trình để giải bài toán phi tuyến có<br />
ràng buộc, được tiến hành ở mỗi lần lấy mẫu<br />
trong phạm vi dự báo P, bao gồm hàng loạt<br />
các di chuyển mà mỗi di chuyển tương ứng<br />
với phạm vi điều khiển. Trong trường hợp<br />
này chỉ một phạm vi điều khiển được chọn để<br />
biểu diễn ràng buộc của phương sai, các ràng<br />
buộc còn lại không được tính đến. Tuy nhiên,<br />
điều này chỉ làm thay đổi về mặt số học chứ<br />
không thay đổi tính chất chung của giải pháp.<br />
Ngoài ra, NMPC với mô hình Gauss còn có<br />
thể được nhận biết như sau:<br />
- Hàm mục tiêu: Hàm mục tiêu được sử dụng<br />
ở công thức (7) chỉ là một lựa chọn. Sự lựa<br />
chọn hàm mục tiêu này có ảnh hưởng lớn đến<br />
khối lượng tính toán.<br />
- Vấn đề tối ưu hóa được nghiên cứu cho ΔU<br />
(k) thay vì U(k): Điều này không chỉ là một<br />
thay đổi về hình thức mà còn tạo ra các dạng<br />
197<br />
<br />
http://www.lrc-tnu.edu.vn<br />
<br />
Lê Thị Huyền Linh và Đtg<br />
<br />
Tạp chí KHOA HỌC & CÔNG NGHỆ<br />
<br />
NMPC khác. Chương trình điều khiển ma<br />
trận động (DMC) sử dụng nguyên tắc chung<br />
cùng với hàm mục tiêu tương ứng được đưa<br />
về dạng bình phương cực tiểu.<br />
<br />
hình toán học và các thông số cần thiết được<br />
trình bày trong [1], [3].<br />
<br />
- Hạn chế về phần mềm: Việc sử dụng thuật<br />
toán tối ưu có ràng buộc đòi hỏi rất nhiều<br />
công việc tính toán trong khi ràng buộc mềm,<br />
còn gọi là trọng số đối với các biến ràng buộc<br />
trong hàm mục tiêu lại có thể được sử dụng<br />
để giảm khối lượng tính toán.<br />
Có nhiều vấn đề quan tâm để áp dụng NMPC.<br />
Một trong số đó là giải pháp hiệu quả về mặt<br />
số học. Chương trình thuật toán phi tuyến tối<br />
ưu hóa yêu cầu khối lượng tính toán lớn. Các<br />
phương pháp ước lượng và các phương pháp<br />
khác chẳng hạn như phương pháp dự đoán kết<br />
quả được sử dụng để giảm khối lượng tính<br />
toán, nhất là trong các trường hợp đặc biệt,<br />
giống như các mô hình tuyến tính hoặc các<br />
hàm mục tiêu đặc biệt.<br />
Sự ổn định của hệ thống vòng lặp kín là vấn<br />
đề tiếp theo. Hiện nay chưa có điều kiện ổn<br />
định nào được suy ra từ mô hình Gauss như là<br />
sự biểu diễn của mô hình xác xuất không<br />
tham số.<br />
MÔ PHỎNG<br />
Quá trình pH<br />
Sơ đồ đơn giản hóa quá trình trung hòa độ pH<br />
[1], [2] được đưa ra trong Hình 2. Quá trình<br />
này bao gồm một dòng axit (Q1), dòng đệm<br />
(Q2) và dòng cơ sở (Q3) được trộn lẫn trong<br />
một bồn chứa T1. Trước khi trộn, đưa axit vào<br />
bồn T2 để tạo thêm động lực dòng chảy. Tỷ lệ<br />
các axit và dòng chảy cơ bản được điều khiển<br />
bằng các van điều khiển dòng, trong khi tỷ lệ<br />
dòng chảy đệm được điều khiển bằng tay với<br />
một lưu lượng kế. Dòng chảy ra pH là biến<br />
đo. Vì bộ cảm biến pH nằm ở hạ nguồn bồn<br />
T1, cho nên việc đo độ pH có một độ trễ (Td).<br />
Trong nghiên cứu này, độ pH được điều khiển<br />
bằng cách điều chỉnh tốc độ dòng chảy cơ<br />
bản. Mô tả chi tiết hơn của quá trình với mô<br />
<br />
Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br />
<br />
86(10): 195 - 200<br />
<br />
Hình 2. Hệ thống trung hòa độ pH<br />
<br />
Mô hình động của hệ thống trung hòa độ pH<br />
thể hiện trong Hình 2 được suy ra từ việc sử<br />
dụng các phương trình bảo toàn và các mối<br />
quan hệ cân bằng. Mô hình này cũng bao gồm<br />
các van và máy phát động lực cũng như các<br />
mối quan hệ thuỷ lực cho các dòng chảy thoát<br />
nước. Những giả định của mô hình bao gồm<br />
sự pha trộn một cách hoàn hảo nồng độ đồng<br />
đều, và sự hòa tan một cách hoàn toàn của các<br />
ion tham gia. Mô hình mô phỏng quá trình pH<br />
đã được sử dụng cho việc tạo ra các dữ liệu<br />
cần thiết, vì thế nó bao gồm rất nhiều các yếu<br />
tố phi tuyến cũng như hàm tính toán tính ẩn<br />
có giá trị đối với đường cong chuẩn độ phi<br />
tuyến cao.<br />
Nhận dạng mô hình<br />
Dựa trên các đáp ứng và chu trình cắt và thử<br />
lặp đi lặp lại một khoảng lấy mẫu thời gian 25<br />
giây đã được lựa chọn (thời gian lấy mẫu này<br />
lớn hơn thời gian chết). Tín hiệu nhận dạng<br />
được chọn của 400 lần lấy mẫu đã được tạo ra<br />
từ một phân bố ngẫu nhiên đồng nhất và tỷ lệ<br />
là 50 giây.<br />
Siêu tham số thu được trong mô hình Gauss ở<br />
lệnh thứ ba là:<br />
<br />
[1 , 2 , 3 , 4 , 5 , 6 , v0 , v1 ]<br />
=[ 6.0505, 2.0823, 0.4785, 5.3388,<br />
8.7080, 0.8754, 5.4164]<br />
(13)<br />
Trong đó những siêu tham số từ 1 đến 3<br />
biểu thị một trọng số cho mỗi đầu ra hồi quy, từ<br />
4 đến 6 biểu thị một trọng số cho mỗi đầu<br />
<br />
198<br />
<br />
http://www.lrc-tnu.edu.vn<br />
<br />
Lê Thị Huyền Linh và Đtg<br />
<br />
Tạp chí KHOA HỌC & CÔNG NGHỆ<br />
<br />
86(10): 195 - 200<br />
<br />
vào hồi quy, v0 là phương sai nhiễu được dự<br />
đoán và v1 là phương sai dọc được dự đoán.<br />
<br />
Hình 4. Đáp ứng của mô hình Gauss đối với tín<br />
hiệu hợp lệ<br />
Hình 3. Đáp ứng của mô hình Gauss với tín hiệu<br />
kích thích được sử dụng cho nhận dạng.<br />
<br />
Có thể thấy rằng mô hình này rất tương thích<br />
đối với việc nhận biết các tín hiệu đầu vào để<br />
sử dụng cho quá trình tối ưu hóa. Tuy nhiên,<br />
mô hình thu được chỉ chứa các thông tin chủ<br />
yếu trong vùng có pH dưới 7 như đã chỉ ra từ<br />
đáp ứng trong hình 3. Các tín hiệu hợp lệ có<br />
tần suất và quy mô thấp hơn với tín hiệu nhận<br />
dạng. Đó là, nếu như mô hình nhận dạng<br />
được kích thích bằng một tín hiệu cao hơn thì<br />
nó chắc chắn phải đáp ứng tốt với cả các tín<br />
hiệu có thành phần thấp hơn. Các tín hiệu hợp<br />
lệ đã được thu bằng máy phát điện của nhiễu<br />
ngẫu nhiên với sự phân phối đồng nhất và tỷ<br />
lệ 500 giây. Hình 4 mô tả đáp ứng của mô<br />
hình Gauss đối với tín hiệu hợp lệ so sánh với<br />
đáp ứng của quá trình. Sự phù hợp của đáp<br />
ứng đối với tín hiệu hợp lệ được biểu diễn<br />
bằng các thông số như sau:<br />
- Giá trị trung bình kiểm tra của sai số tuyệt<br />
đối: AE = 0.1276<br />
- Giá trị trung bình kiểm tra của sai số bình<br />
phương: SE = 0.0373<br />
<br />
Sau khi đạt được mô hình hợp lệ, mô hình<br />
này được dùng để thiết kế bộ điều khiển.<br />
<br />
Hình 5. Hình trên là đáp ứng điều khiển dựa trên<br />
mô hình Gauss, hình dưới là tín hiệu điều khiển<br />
trong trường hợp không có ràng buộc.<br />
<br />
Điều khiển<br />
Thuật toán điều khiển được mô tả ở trên dùng<br />
để mô phỏng điều khiển quá trình trung hoà<br />
độ pH. Quỹ đạo tham chiếu r được xác định<br />
sao cho nó tiệm cận với điểm đặt của hàm số<br />
mũ từ giá trị đầu ra hiện tại. Như đã nói ở trên<br />
điểm trùng hợp được chọn là 8 lần lấy mẫu và<br />
phạm vi điều khiển là một lần lấy mẫu.<br />
Những kết quả của điều khiển không ràng<br />
buộc được đưa ra ở hình 5 và 6.<br />
<br />
- Sai số mật độ: LD = 1.9889<br />
<br />
Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br />
<br />
199<br />
<br />
http://www.lrc-tnu.edu.vn<br />
<br />