Đề xuất phương pháp nhận dạng đa thông số của hệ thống được mô tả bởi phương trình đạo hàm riêng dựa trên phương pháp gradient phối ngẫu hiệu chỉnh

10 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
ĐỀ XUẤT PHƯƠNG PHÁP NHẬN DẠNG ĐA THÔNG SỐ CỦA HỆ 
THỐNG ĐƯỢC MÔ TẢ BỞI PHƯƠNG TRÌNH ĐẠO HÀM RIÊNG DỰA 
TRÊN PHƯƠNG PHÁP GRADIENT PHỐI NGẪU HIỆU CHỈNH 
PROPOSING MULTI-PARAMETER IDENTIFICATION OF THE SYSTEM 
DESCRIBED BY PARTIAL DIFFERENTIAL EQUATIONS BASED ON 
MODIFIED CONJUGATE GRADIENT METHOD 
Trần Thanh Phong1,*, Nguyễn Hoàng Phương1,2 
1 Trường Đại học Tiền Giang, Việt Nam 
2
 Trường Đại học Sư phạm Kỹ thuật Tp. Hồ Chí Minh, Việt Nam 
Ngày toà soạn nhận bài 23/12/2019, ngày phản biện đánh giá 11/2/2020, ngày chấp nhận đăng 19/2/2020. 
TÓM TẮT 
Bài báo giới thiệu phương pháp nhận dạng đồng thời nhiều thông số khiếm khuyết của hệ 
thống được mô tả bởi phương trình đạo hàm riêng bậc hai với phương trình truyền nhiệt là 
một ví dụ. Theo đó, một nguồn nhiệt di động trên bề mặt một tấm nhôm được xem xét với một 
nhóm cảm biến cố định đặt trên khu vực khảo sát để đo sự tiến triển của nhiệt độ theo thời 
gian và không gian. Việc giải bài toán ngược này đòi hỏi các dữ liệu đầu vào phải tối ưu để 
giảm thời gian tính toán bằng việc loại bỏ các giá trị đo đạc của các cảm biến không hữu 
dụng. Một giải thuật lựa chọn cảm biến được đề xuất kết hợp với phương pháp gradient phối 
ngẫu để giúp nhận dạng hiệu quả. Hơn nữa, giải thuật lặp cũng được hiệu chỉnh bằng việc đề 
xuất giải thuật lựa chọn cửa sổ chập và trượt linh hoạt để tối ưu thời gian nhận dạng. Nhằm 
giúp đánh giá hiệu quả của phương pháp đề xuất, giá trị nhiệt độ đo đạc với tác động bởi các 
nhiễu tuân theo hàm phân phối chuẩn Gauss. Kết quả cho thấy phương pháp được đề xuất có 
khả năng nhận dạng tốt hàm mật độ dòng nhiệt và quỹ đạo di chuyển với độ trễ thấp và sai số 
đáp ứng yêu cầu đặt ra. 
Từ khóa: Bài toán ngược; nguồn nhiệt; nhận dạng thông số; phương pháp gradient phối 
ngẫu; phương trình đạo hàm riêng. 
ABSTRACT 
The paper introduces the method of simultaneous identification of multi unknown 
parameters of the system described by the quadratic partial differential equation with the heat 
transfer equation as an example. Accordingly, a mobile heat source on the aluminium plate is 
considered with a fixed group of sensors located on the survey area to measure the evolution 
of temperature over time and space. Solving this inverse problem requires optimal input data 
to reduce computational time by eliminating the measurement values of useless sensors. A 
sensor selection algorithm is proposed in combination with the conjugate gradient method for 
effective identification. Moreover, the iterative algorithm is also corrected by proposing a 
flexible shutter and sliding window selection algorithm to optimize identification time. In 
order to help evaluate the effectiveness of the proposed method, the measured temperature 
value with the effect of disturbances follows the standard Gaussian distribution function. The 
results show that the proposed method has good ability to identify the function of heat flow 
density and trajectory with low latency and error to meet the requirements. 
Keywords: Heat source; identification; inverse problems; conjugate gradient method; partial 
differential equations. 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
11 
1. ĐẶT VẤN ĐỀ 
Trong thực tế, các hiện tượng vật lý 
thường gặp đều có thể mô hình hóa bằng 
phương trình toán học bởi một hệ phương 
trình vi phân từng phần bậc hai (parabol) 
hoặc bậc ba (hyperbol). Đặc biệt là các hiện 
tượng vật lý liên quan đến thời tiết, đám cháy 
rừng, vết dầu loang trên mặt nước biển, việc 
xả thải của các khu công nghiệp. Trong bài 
báo này, tác giả sử dụng phương trình truyền 
nhiệt tổng quát trong vật dẫn như một ví dụ. 
Nghiên cứu được đề xuất nhằm phát 
triển một phương pháp để nhận dạng đồng 
thời nhiều thông số bất định của hệ thống bất 
kỳ được mô tả bởi phương trình đạo hàm 
riêng dựa trên phương pháp gradient phối 
ngẫu với việc sử dụng nhóm cảm biến hữu 
dụng [1, 2]. Các cảm biến này được xác định 
nhờ vị trí của nhóm cảm biến cố định được 
thiết lập tại khu vực di chuyển của nguồn 
nhiệt bằng một giải thuật dựa trên mối quan 
hệ giữa vị trí nguồn nhiệt và cảm biến. Việc 
này nhằm giảm thời gian tính toán của hệ 
thống. 
Vấn đề giải bài toán ngược của hệ thống 
được mô tả bởi phương trình đạo hàm riêng, 
cụ thể là quá trình truyền nhiệt được đề xuất 
một cách không đầy đủ bởi Hadamard với 
việc giải bài toán bằng phương pháp lặp [3]. 
Phương pháp gradient phối ngẫu (CGM) 
được sử dụng để tối thiểu hóa sai lệch bằng 
kỹ thuật lặp và chứng tỏ đây là phương pháp 
ổn định để ước lượng các thông số [4]. Hạn 
chế của các nghiên cứu này là chỉ nhận dạng 
riêng lẻ từng thông số. Trong khi, việc truyền 
nhiệt trong không gian đa chiều rất phức tạp 
nên làm cho mô hình toán của hệ thống trở 
nên cồng kềnh và mất nhiều thời gian để xử 
lý [5, 6]. Để cải thiện hiệu năng của phương 
pháp, một giải thuật lặp có hiệu chỉnh được 
đề xuất để nhằm rút ngắn thời gian tính của 
quá trình nhận dạng các thông số bất định 
của hệ thống cần xem xét [4], [5], [7]. 
Theo đó, nghiên cứu này sẽ đề xuất một 
giải thuật để ước lượng thông số bất định của 
hệ thống dựa trên phương pháp gradient phối 
ngẫu để nhận dạng một cách đầy đủ theo 
hướng của Hadamrd kết hợp với giải thuật 
lựa chọn cảm biến tối ưu và cửa sổ trượt linh 
hoạt [6, 8, 9]. Cấu trúc của bài báo này bao 
gồm: mô tả hệ thống thực nghiệm, mô hình 
toán của vấn đề bài toán ngược dựa trên 
phương pháp CGM, phương pháp cửa sổ 
trượt kết hợp với giải thuật xác định cảm 
biến tối ưu sẽ được giới thiệu trong phần nội 
dung nghiên cứu. Cuối cùng là kết quả 
nghiên cứu, thảo luận và phần kết luận. 
2. NỘI DUNG NGHIÊN CỨU 
2.1 Mô tả hệ thống 
Việc nghiên cứu và mô phỏng quá trình 
truyền nhiệt tổng quát được thực hiện bằng 
việc xây dựng một mô hình thí nghiệm cho 
phương pháp nhận dạng đã đề xuất. Theo đó, 
giả thuyết rằng một nguồn nghiệt di động có 
hàm mật độ ( )t di chuyển theo quỹ đạo 
( )I t trên bề mặt một tấm nhôm hình vuông 
có kích thước cạnh bên L và độ dày là e . 
Giới hạn biên của miền làm việc được ký 
hiệu 2 R [7, 10, 11]. Nó được đốt nóng 
bởi nguồn nhiệt được giả định bằng đĩa đồng 
chất ( ),D I t r có tâm ( )I t và bán kính 
Ir . 
Hàm phân bố nhiệt độ của tấm kim loại 
 , ,x y t là hàm liên tục theo không gian và 
thời gian được tính bằng Kelvin. Biến số 
không gian của hệ thống ,x y  được 
tính bằng mét và biến số thời gian là t  
được tính bằng giây. Giả định rằng các giá trị 
của các thông số của hệ thống để xây dựng 
cho mô hình thí nghiệm đều biết trước và 
được liệt kê trong Bảng 1. 
Quỹ đạo di chuyển của nguồn nhiệt được 
mô tả như Hình 2. Đồng thời, hàm mật độ 
dòng nhiệt của các nguồn được cho bởi hàm 
có đồ thị được thể hiện như Hình 3. Biểu 
thức hàm mật độ công suất nhiệt tổng của 
nguồn , ,x y t để đốt nóng tấm nhôm 
thực nghiệm được diễn đạt như sau: 
( ) ,
, ,
0 ,
t khi x y D
x y t
khi x y D

 (1) 
12 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
Theo một cách khác, biểu thức , ,x y t 
còn có thể được biểu diễn một cách liên tục 
và khả vi dưới dạng hàm tổng hợp của các 
hàm mật độ thành phần theo biến thời gian 
và theo các tọa độ trong không gian như sau: 
( )
, , arccot ( , , )
t
x y t x y t r

  
 (2) 
với 
2 2
( , , ) ( ) ( )s sx y t x x t y y t  
Trong đó, hệ số  R được chọn 
nhằm mục đích rời rạc hóa hàm mật độ dòng 
nhiệt liên tục. Khoảng thời gian  có thể 
được chia ra thành n đoạn  
1
1
0
,
n
i i
i
T t t
với it i và bước chia rời rạc hóa được 
định nghĩa bởi T n . Để tránh làm mất 
tính tổng quát, phương trình quỹ đạo của tất 
cả các vị trí định vị bất kỳ của các nguồn 
nhiệt ( ), ( )s sI x t y t cũng được thành lập lại 
dưới dạng các hàm rời rạc một cách tuyến 
tính và được viết lại bằng cách sử dụng các 
hàm nón cơ bản ( )
is t với 0,1, ,i n : 
Hình 1. Biểu diễn hàm nón cơ bản 
Hàm mật độ dòng nhiệt được diễn đạt lại 
bởi ( ) ( )i it s t  và quỹ đạo di chuyển của 
nguồn nhiệt cũng được diễn đạt lại bởi 
( ) ( )i is sx t x s t , ( ) ( )
i i
s sy t y s t . 
2.2 Vấn đề thuận 
Nếu tất cả các thông số được biết trong 
bảng phụ lục, sự tiến triển của nhiệt độ trong 
không gian và thời gian là kết quả nghiệm 
của phương trình đạo hàm riêng bậc hai: 
0
( )
( ) ( )
( , ,0) ( , )
( )
0
c P
t
x y x y
n
 
     
 
 
 
 
  
  
 
 
 
 (3) 
với ( ) ( , , )x y t , 
( ) 2 ( )
( )
h
P
e
   

 và
2 2
2 2
( ) ( )
( )
   
 
 
 x y
 là toán tử Laplace 
của nhiệt độ theo không gian và thời gian. 
Trong phương trình trên, điều kiện đầu của 
phương trình vi phân được xem là nhiệt độ 
của môi trường xung quanh. 
Hình 2. Hàm mật độ dòng nhiệt 
Hình 3. Quỹ đạo di chuyển của nguồn nhiệt 
-0.4 -0.2 0 0.2 0.4
-0.2
-0.1
0
0.1
0.2
0.3
X [m]
Y
 [
m
]
Trajectory of source Sensor Ck
0 500 1000 1500
1
1.5
2
2.5
3
3.5
4
x 10
4
Time [seconde]
F
lu
x
 d
e
n
s
it
y
 [
W
/m
2
]
Real flux Flux to be estimated Initial flux
 τ 2τ 3τ τ (n-1)τ nτ 
1 
t [s] 
S0 S1 S2 Si-1 Si Si+1 Sn-1 Sn 
Si(t) 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
13 
Điều kiện biên phụ thuộc vào hàm mật độ 
dòng nhiệt và hệ số đối lưu nhiệt tự nhiên h. 
Nó rất khó đo đạc trực tiếp và thường được 
xác định thông qua kinh nghiệm hoặc mô hình 
và được chọn sao cho độ nhạy hàm phân bố 
nhiệt độ thay đổi ít nhất so với sự biến thiên 
của các thông số khác của hệ thống. 
Kết quả số thu được nhờ sử dụng phương 
pháp phần tử hữu hạn (FEM, Finite Element 
Method) thực hiện trong phần mềm COMSOL 
MultiphysicsTM được nhúng vào phần mềm 
Matlab® [12-18]. Quá trình phân bố nhiệt độ 
trên tấm nhôm theo thời gian được thể hiện tại 
các thời thời điểm như Hình 4. 
Để đánh giá độ tin cậy của mô hình toán 
đề xuất để đơn giản hóa phương trình truyền 
nhiệt trong không gian hai chiều, 25 cảm 
biến nhiệt Ck được đặt cố định trên tấm kim 
loại nhằm mục đích thu thập dữ liệu nhiệt độ 
của điểm đặt cảm biến trong suốt quá tình 
thực nghiệm. Hơn nữa, để đánh giá ảnh 
hưởng của các sai số trong quá trình đo đạc, 
giả định rằng nhiệt độ thu thập từ các cảm 
biến đã bị tác động bởi các nhiễu, tuân theo 
hàm phân phối xác suất Gauss   2( , )N 
với giá trị trung bình 0 và độ lệch 
chuẩn 1 . 
(a) Phân bố nhiệt độ tại t=300s (b) Phân bố nhiệt độ tại t=600s 
(c) Phân bố nhiệt độ tại t=900s (d) Phân bố nhiệt độ tại t=1500s 
Hình 4. Phân bố nhiệt độ trên tấm kim loại theo thời gian 
Hình 5. Phân bố nhiệt độ của các cảm biến theo thời gian 
0 250 500 750 1000 1250 1500
280
315
350
385
420
455
490
525
560
Time [s]
T
e
m
p
e
ra
tu
re
 [
K
]
Temperature evolution in time
14 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
2.3 Phương pháp gradient phối ngẫu 
2.3.1 Hàm mục tiêu 
Để nhận dạng hàm mật độ dòng nhiệt 
( )t , xét rằng “nhiệt độ đo” ˆ( , )kc t tại vị trí 
cảm biến 1,2,...,25kc , một vấn đề ngược có thể 
được thiết lập và giải nghiệm bằng việc tối 
thiểu hóa một tiêu chuẩn bậc hai: 
2
10
1 ˆ( , , ) ( , )
2
T n
k k
k
J c t c t dt 
  (4) 
Trong phần tiếp theo, việc thiết lập mô 
hình toán của các vấn để sẽ được giới thiệu 
nhằm tính toán các thông số trung gian của 
phương pháp nhận dạng thông số bất định hệ 
thống dựa trên phương pháp CGM. 
2.3.2 Vấn đề độ nhạy 
Xét rằng độ thay đổi của nhiệt độ 
( , , )x y t được sinh ra bởi sự thay đổi của 
mật độ dòng nhiệt tổng được cho bởi: 
 , , , , , ,x y t x ty xt y  . Nên: 
0
, , , ,
, , lim
x y t x y t
x y t

 


 (5) 
với , , , , , ,x y t x y t x y t   . 
Vấn đề độ nhạy của hệ thống được mô tả 
bởi hệ phương trình sau: 
0
( )
( ) ( )
( , ,0) ( , )
( )
0
c P
t
x y x y
n
 
      
 
 
 
 
  
  
 
 
 
 (6) 
với 
( ) 2 ( )
( )
h
P
e
   
 
 . 
Trong đó, sự thay đổi của hàm mật độ là: 
( ) ( )
( ) ( ) ( )
( ) ( )
t I t
t I t
  
   

 
 
 (7) 
Kết hợp với phương trình (2), ta có: 
( )
arccotan ( )
t
t r
t
  
   
 


(8) 
( )
( ) ( )( ( )) ( )
( )
s s s
s
x t H x x t x t
x t
 
   


 (9) 
( )
( ) ( )( ( )) ( )
( )
s s s
s
y t H y y t y t
y t
 
   


 (10) 
với 
2
( )
( )
( ) 1 ( ( ) )
t
H
r


    
Từ đó, suy ra: 
( )
( ) arccot ( )
(11)
 ( ) ( ( )) ( ) ( ( )) ( )s s s s
t
r
x x t x t y y t y t

    
  
  
Nghiệm của phương trình (6) cho phép 
xác định giá trị thay đổi của giá trị nhận dạng 
của các thông số đối với chu kỳ kế tiếp và 
được tính dựa trên sự thay đổi của hàm mật 
độ tổng tại vòng lặp thứ k, như sau: 
1 1
1
k k k
k d 
 (12) 
Trong đó, đại lượng 
1k
 
 được xác định 
nhờ vào việc tối thiểu hóa hàm mục tiêu với: 
2
1
10
1 ˆ, ( ) 
2 k k
T n
k
c c
k
J t t dt 
  (13) 
Đặt 1 ˆ( ) , ( )k k
k
c cd t t t  
 , từ đó suy ra: 
 1 1, , ,k k k
k k k
k
c c ct t t     
 (14) 
với ,k
k
c t là sự thay đổi của nhiệt độ 
sinh ra từ sự thay đổi của hàm mật độ dòng 
nhiệt trong chiều hướng thay đổi 
1k
d 
. Thay 
thế các giá trị vào phương trình (12) và khi 
đó, giá trị 
1k
 
 được xác định bởi: 
 11 Arg min , kk J

   
 (15) 
Đồng nghĩa với việc giải phương trình 
 1
1
0
,
k
k
J
 



, có thể được viết dưới dạng: 
2
1
10
ˆ( ) ( ) ( ) 0
k k
T n
k
c c
k
t t d t dt    
  (16) 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
15 
Giải phương trình này, ta nhận được độ lớn 
của hệ số thay đổi 
1k
 
 cho chu kỳ k+1: 
11 0
2
10
ˆ( ) ( ) 
ˆ ( ) 
k
k
T n
c
kk
T n
c
k
d t t dt
t dt
 


 
 
 (17) 
Việc xác định các giá trị cho quỹ đạo 
cũng được xác định một cách tương tự. Hơn 
nữa, việc giải vấn đề độ nhạy phụ thuộc vào 
vector chiều thay đổi cho vòng lặp kế tiếp và 
được xác định bởi vấn đề phụ trợ. 
2.3.3 Vấn đề phụ trợ 
Vector chiều thay đổi 
1k
d 
 của giá trị 
cần nhận dạng được xác định thông qua các 
gradient của hàm mục tiêu cho mỗi vòng lặp 
của giải thuật tối thiểu hóa được xác định bởi: 
; ;
i i i
s s
J J J
J
x y 
   
 
   
, 1, 2,...,i n (18) 
Để làm được điều này, công thức Lagrange 
( ( ), , )   được định nghĩa: 
0
( ( ), , ) ( ( ), ) ( )
ft
J A d dt     

  
(19) 
với 
( )
( ) ( ) ( )A c P
t
 
    


. 
Kết hợp với phương trình (3) và hàm số 
nhân ( )  được lựa chọn và xác định sao 
cho 
( ( ), , ( ))
( ) 0
( )
   
 
 
 

 và là 
nghiệm của vấn đề phụ trợ sau đây: 
( )
( ) ( )
( , ,0) 0 ( , )
( )
0
c E
t
x y x y
n
 
     

 
 
 
  
  
 
 
 
 (20) 
Với 
1
( ) ( ) ( , ) 2 ( ) /
k k
n
D c c
k
E d t x y h e    
 
và ( , ) ( ) ( )
k k k kD c c D c D c
x y x x y y   là 
hàm phân phối Dirac. Cuối cùng, ( )  là 
kết quả của vấn đề vừa nêu và khi đó ta có: 
 ( , , ), , ( , , ),x y t J x y t      (21) 
Từ phương trình (21) ta có: 
0
( ) ( )
arccot ( )
ft is t
J r d dt
e
 
  

  (22) 
0
( )
( ) ( ( )) ( )
f
i
ss
t
i
s xx
J H x x t s t d dt
e
 


  (23) 
0
( )
( ) ( ( )) ( )
f
i
ss
t
i
s yy
J H y y t s t d dt
e
 


  (24) 
Từ những công thức của gradient trên, 
chiều hướng tăng giảm sẽ được ước lượng cho 
mỗi vòng lặp mới bởi: 
2
1
2
1
( , )
( , )
( , )
k
k k
k
k
J
d J d
J
 
 
 

  

 (25) 
với , ,s sx y và  là module chuẩn 
Euclidean. Dựa các kết quả trên, giá trị dự 
báo của hàm mật độ dòng nhiệt và quỹ đạo 
của nguồn nhiệt di động được xác định bởi : 
1( ) ( ) k kk k dt t   
 với , ,s sx y (26) 
Từ mô hình toán học được đề xuất cho 
vấn đề nghịch như trên, tác giả đề xuất có thể 
nghiên cứu ảnh hưởng lẫn nhau về toán học 
giữa các biến số khác nhau của hệ thống 
bằng việc lấy đạo hàm riêng của hàm được 
mô hình hóa của hệ thống. Ngoài ra, cần phải 
tính đến tác dụng đối lưu và dẫn nhiệt của 
nhiệt trong miền khảo sát. Sau khi xác thực 
mô hình này, tác giả có thể khái quát hệ 
thống từ một số nguồn hoặc suy luận cho một 
số trường hợp cụ thể bằng cách loại các dữ 
liệu gây bất lợi (như là để xác định mật độ 
dòng nhiệt riêng, xác định quỹ đạo di chuyển 
của nguồn nhiệt hoặc xác định đồng thời cặp 
mật độ dòng nhiệt và quỹ đạo di chuyển của 
nguồn). 
Trong phần sau, áp dụng phương pháp 
nhận dạng trực tuyến dựa trên thuật toán lặp 
chính quy của gradient liên hợp sử dụng mạng 
16 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
cảm biến cố định để xác định các thông số của 
hệ thống thông qua các bài toán như sau: 
(1) Nhận dạng ( )t 
(2) Nhận dạng ( ( ), ( ))I x t y t 
(3) Nhận dạng ( ) & ( ( ), ( ))t I x t y t 
2.4 Nhận dạng thông số bằng cửa sổ trượt 
kết hợp giải thuật xác định cảm biến 
tối ưu 
Trong nghiên cứu [16,17], các tác giả đã 
nghiên cứu sự quan tâm việc thích ứng phương 
pháp gradient liên hợp gần đúng với nhận dạng 
thông số dựa trên phương pháp lặp. Thời gian 
tính toán phụ thuộc quan trọng vào độ phức tạp 
của mô hình và do đó lượng dữ liệu nhập vào 
thuật toán. Để có kết quả theo cách gần như 
trực tuyến, tác giả đề xuất việc chọn số lượng 
dữ liệu tối ưu, hay xác định khoảng thời gian 
tốt nhất cho mỗi quy trình nhận dạng là cần 
thiết (xem Hình 6). Thay vì chạy giải thuật 
nhận dạng trên toàn bộ khoảng thời gian thực 
nghiệm thì giải thuật lựa chọn khoảng thời gian 
tối ưu trong cửa sổ chập và độ lớn cửa sổ này 
có thể thay đổi một cách linh hoạt. Do đó, việc 
nhận dạng gần như chạy trực tuyến này dựa 
trên các cửa sổ thời gian chập và trượt. 
Hình 6. Giải thuật xác định cửa sổ chập 
Áp dụng thuật toán ước lượng giá trị của 
thông số ở mỗi khoảng thời gian chập được 
chọn để xác định cơ bản như sau: (i) mật độ 
của nguồn nhiệt, (ii) quỹ đạo của nguồn nhiệt. 
Thuật toán này lặp lại cho đến khi điều kiện 
dừng giải thuật được thỏa mãn. Lưu đồ sau 
đây cho thấy các bước của quy trình xác định 
cặp mật độ và quỹ đạo của một nguồn trong 
khoảng thời gian. 
Giải thuật nhận dạng trực tuyến CGM 
Bước 1 : Thiết lập thông số 
 Xác định của chập ,i i iT  
 . 
 Vector trạng thái: 0 ( )k t , 0 ( )kI t ; 0 0k 
 Mode ( ( ); ( ( ), ( )))j j jt I x t y t , 1k . 
Bước 2 : Thực thi giải thuật 
while 
  Vấn đề thuận và hàm mục tiêu 
 Nạp dữ liệu cảm biên ( , , ( ), )k knC t t I  
 Tính hàm mục tiêu ( ( ); ( ))k kJ t I t . 
If 
max( )|| || i istopJ J N tk T T 
Kết thúc giải thuật 
 Vấn đề phụ trợ và gradient 
 If ( 1)mode 
k kJ J  với 0
k
IJ , ( ) 0.I t 
Else 
k k
IJ J  với 0
kJ , ( ) 0.t 
 Tìm vector hướng 
1k k k kd J d  
 Vấn đề độ nhạy và độ thay đổi 
 Tính giá trị , ,x y t với hướng 1kd . 
 Tính 1 Argmin k J ( ).  
 Ước lượng giá trị 
 If ( 1)mode 
1 1 1k k k kd  
 và 0.mode 
Else 
1 1 1k k k k
II I d
 và 1.mode 
 1 0k và 1k k : về Bước 2.
Các giải pháp tính toán số được trình bày 
trong bài viết này được dựa trên phương pháp 
phần tử hữu hạn với phần mềm COMSOL 
MultiphysicsTM được nhúng trong phần mềm 
Matlab®. Kết quả của việc nhận dạng được 
trình bày trong phần sau. 
Bắt đầu 
Thử độ lớn của sổ chập 
Tính toán giá trị hàm mục tiêu 
Đọc dữ liệu nhiệt độ 
Kết thúc 
Cửa sổ chập 
Đúng 
Sai 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
17 
2.5 Kết quả và thảo luận 
Mục đích của nghiên cứu này là ước 
lượng giá trị của các thông số hệ thống được 
mô tả bởi phương trình đạo hàm riêng, ứng 
dụng cho phương trình truyền nhiệt trong 
không gian là ví dụ. Hàm mật độ dòng nhiệt 
và quỹ đạo của nguồn nhiệt di động được đề 
xuất để xây dựng mô phỏng thực nghiệm 
nhằm để chứng minh hiệu quả của phương 
pháp đề xuất. Dựa trên việc tối thiểu hóa sai 
số dữ liệu ngõ vào từ các cảm biến và giá trị 
từ mô hình với nhiễu được tuân theo phân 
phối chuẩn Gauss   2( , )N . 
Để kiểm nghiệm hiệu quả phương pháp 
nhận dạng thông số của hệ thống được mô tả 
bởi phương trình đạo hàm riêng (trường hợp 
nhận dạng hàm mật độ dòng nhiệt và quỹ đạo 
di chuyển). Tác giả tiến hành xây dựng mô 
hình thí nghiệm với một nguồn nhiệt có hàm 
mật độ ( )t và quỹ đạo di chuyển I(x(t), y(t)). 
Nguồn nhiệt di chuyển trên tấm nhôm 
hình vuông, kích thước các cạnh 1L m và 
độ dày 2e mm . Thời gian thực nghiệm 
1500ft s với thời gian rời rạc hóa 20s 
(hay suy ra, có 76x3 hệ số cần nhận dạng). 
Các thông số thiết lập hệ thống thí nghiệm 
được cho bởi Bảng 1 như sau. 
Bảng 1. Giá trị các thông số hệ thống 
Ký hiệu Đơn vị Giá trị 
c 3 1Jm K 2,34.10
4
h 2 1Wm K 15 
ft s
1500 
 1 1Wm K 160 
max 
2Wm 
 3.10
4
0 K 291 
L m 1 
e m 2.10
-3
r m 6.10
-2
n - - 
Để thiết lập điều kiện đầu cho giải thuật 
tại bước k=0, giả thuyết rằng giá trị rời rạc 
ban đầu của hàm mật độ dòng nhiệt 
ϕk=0=1,5.103(W/m2). Kết quả nhận dạng ( )t 
bằng phương pháp CGM sử dụng giải thuật 
lặp kết hợp với phương pháp cửa sổ trượt 
được mô tả trong Hình 7. 
Hình 7. Kết quả nhận dạng trên cửa sổ chập 
Kết quả nhận dạng ( )t trong cửa sổ 
chập  81, 222iT (đường màu đỏ là giá trị 
thật, đường màu xanh là giá trị nhận dạng cho 
bởi CGM). Theo đó, giá trị ( , 81)t t là kết 
quả đã nhận dạng và giá trị ( , 222)t t là 
giá trị dự báo cho các vòng lặp tiếp theo. 
Kết quả nhận dạng hàm mật độ dòng 
nhiệt và quỹ đạo di chuyển của nguồn nhiệt 
di động bằng phương pháp gradient phối 
ngẫu kết hợp giải thuật lựa chọn cảm biến tối 
ưu và xác định cửa sổ trượt linh hoạt được 
thể hiện lần lượt trong Hình 8 và Hình 9. 
Hình 8. Kết quả nhận dạng hàm mật độ 
Hình 9. Kết quả nhận dạng quỹ đạo 
0 500 1000 1500
1
1.5
2
2.5
3
3.5
4
x 10
4
Time [seconde]
F
lu
x 
de
ns
ity
 [W
/m
2 ]
Real flux Flux to be estimated
(t) [81:222] windows of 141 s
0 150 300 450 600 750 900 1050 1200 1350 1500
1
1.5
2
2.5
3
3.5
4
x 10
4
Time [s]
Fl
ux
 d
en
si
ty
 [W
/m
2 ]
Real flux
Estimated flux
-0.4 -0.2 0 0.2 0.4
-0.2
-0.1
0
0.1
0.2
0.3
X [m]
Y
 [m
]
-0.4 -0.2 0 0.2 0.4
-0.2
-0.1
0
0.1
0.2
0.3
X [m]
Y
 [
m
]
data1
data2
data3
data4
Real trajectory
Estimated trajectory
18 
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
Để đánh giá hiệu quả của phương pháp 
được đề xuất, tác giả sử dụng các tiêu chí 
như: thời gian tính là thời gian cần thiết để 
xác định tất cả 228 hệ số; trung bình sai số là 
trung bình độ lệch giữa nhiệt độ đo đạc được 
và nhiệt độ từ mô hình dựa trên giá trị nhận 
dạng; độ lệch chuẩn sai số là giá trị trung 
bình bình phương độ lệch giữa nhiệt độ đo 
đạc được và nhiệt độ từ mô hình, và độ trễ 
của quá trình nhận dạng. 
Trung bình sai số (K): 
1
1 ˆ ˆ(t) ( , t)
c
f
N
residus n
nc t
C dt
N
  
  (27) 
Độ lệch chuẩn (K): 
2
1
1 ˆ ˆ(t) ( , t)
c
f
N
residus n
nc t
C dt
N
  
  (28) 
Hiệu quả phương pháp CGM hiệu chỉnh 
được đề xuất cho việc nhận dạng thể hiện 
trong Bảng 2. 
Bảng 2. Giá trị các kết quả nhận dạng 
Thông số ( )t ( )I t ( ) & ( )t I t 
identift 1.741 2.212 2.340 
residus 0,501 0,4997 0,481 
residus 0,979 0,952 0,912 
delayt 309,1 327,5 350,9 
Số liệu thống kê ở Bảng 2 cho thấy rằng: 
giá trị của trung bình sai số và độ lệch chuẩn 
của các thí nghiệm nhận dạng đều tiệm cận 
với giá trị trung bình và độ lệch chuẩn của 
hàm phân phối xác suất Gauss. Điều này 
chứng minh phương pháp nhận dạng trực 
tuyến thông số hệ thống đề xuất dựa trên 
phương pháp đề xuất đáng tin cậy và hiệu quả. 
Hơn nữa, giải thuật lựa chọn cảm biến kết 
hợp với phương pháp CGM thể hiện ưu điểm 
trong nhận dạng thông số của hệ thống với 
thời gian tính toán là 1741s của nhận dạng 
hàm mật độ dòng nhiệt với độ trễ trung bình 
của quá trình nhận dạng là 309s. Nghĩa là, các 
kết quả nhận dạng hàm mật độ dòng nhiệt 
nhận được 309s sau khi kết thúc (Hình 10). 
Hình 10. Hàm trễ của giải thuật nhận dạng 
Bảng 2 cũng cho thấy rằng, số lượng 
biến cần nhận dạng (hay độ phức tạp tăng 
lên) thì thời gian tính toán cần thiết tăng lên 
là 2.212s và 2.340s tương ứng với nhận dạng 
quỹ đạo và nhận dạng đồng thời hàm mật độ, 
quỹ đạo. Theo đó, độ trễ của giải thuật cũng 
tăng theo tương ứng là 327,5s và 350,9s. 
3. KẾT LUẬN 
Bài báo đã tập trung giới thiệu một 
phương pháp nhận dạng đồng thời nhiều 
thông số khiếm khuyết của hệ thống được mô 
tả bởi phương trình đạo hàm riêng bậc hai 
với phương trình truyền nhiệt được đề cập 
đến như một ví dụ. Theo đó, phương pháp 
nhận dạng trực tuyến hàm mật độ dòng nhiệt, 
quỹ đạo và nhận dạng đồng thời hàm mật độ 
dòng nhiệt và quỹ đạo nguồn nhiệt di động 
được đề xuất dựa trên phương pháp gradient 
phối ngẫu hiệu chỉnh kết hợp với giải thuật 
lặp trên cửa sổ chập linh hoạt và giải thuật 
lựa chọn cảm biến hữu dụng. Kết quả cho 
thấy rằng phương pháp đề xuất mang tính 
hiệu quả và tin cậy cao cho quá trình nhận 
dạng. Kết quả này cho phép phát triển xây 
dựng thí nghiệm để nhận dạng thông số của 
hệ thống bất kỳ được mô tả bởi phương trình 
đạo hàm riêng bậc cao. 
TÀI LIỆU THAM KHẢO 
[1] Y. Jarny, MN. Ozisik, JP. Bardon, A general optimization method using adjoint equation 
for solving multidimensional inverse heat conduction, International Journal of Heat and 
Mass Transfer, 34-11, pp. 2911-2919, 1991. 
0 250 500 750 1000 1250 1500
0
100
200
300
400
500
600
Time [s]
D
el
ay
 [s
]
Delay time
Tạp Chí Khoa Học Giáo Dục Kỹ Thuật Số 58 (06/2020) 
Trường Đại Học Sư Phạm Kỹ Thuật TP. Hồ Chí Minh 
19 
[2] Keith A. Woodbury, Inverse Engineering Handbook: Handbook Series for Mechanical 
Engineering, Editor CRC Press, ISBN 9780849308611, pp. 480, 2002. 
[3] Oleg M. Alifanov, Inverse Heat Transfer Problems , International Series in Heat and 
Mass Transfer, Editor Springer Science & Business Media, ISBN9783642764363, pp. 
348, 2012. 
[4] M. Prud’homme, TH Nguyen, On the iterative regularization of inverse heat conduction 
problems by conjugate gradient method, International Communications in Heat and Mass 
Transfer, Vol. 25, 1998. 
[5] S. Beddiaf, et al., Time-dependent heat flux identification: Application to a three-dimensional 
inverse heat conduction problem, Proceedings of International Conference on Modelling, 
Identification and Control, Wuhan, Hubei, China, 2012, pp. 1242-1248, 2012. 
[6] S. Beddiaf, et al., Parametric identification of a heating mobile source in a three 
dimensional geometry, Inverse Problems in Science and Engineering, Vol. 23-1, pp. 
93-111, 2015. 
[7] Gillet, L.P., et al., Implementation of a conjugate gradient algorithm for thermal 
diffusivity identification in a moving boundaries system. Journal of Physics: Conference 
Series, 135(1), p.12082, 2008. 
[8] S Beddiaf, et al., Simultaneous determination of time-varying strength of the heat flux and 
location of a fixed source in a three-dimensional domain, Inverse Problems in Science 
and Engineering, 22-1-2, pp. 166-183, 2014. 
[9] C. Huang, W. Chen, A 3D inverse forced convection problem in estimating surface heat 
flux by conjugate gradient method, International Journal of Heat and Mass Transfer, 43, 
pp. 317-3181, 2000. 
[10] Lefèvre, F., & Le Niliot, C., Multiple transient point heat sources identification in heat 
diffusion: application to experimental 2D problems. International Journal of Heat and 
Mass Transfer, 45(9), 1951-1964, 2002. 
[11] Martin, T. J., & Dulikravich, G. S., Inverse Determination of Boundary Conditions and 
Sources in Steady heat conduction with heat generation. Journal of Heat Transfer, 118(3), 
546-554, 1996. 
[12] Coles, C., & Murio, D. A., Simultaneous space diffusivity and source term reconstruction 
in 2D IHCP. Computers & Mathematics with Applications, 42(12), 1549-1564, 2001. 
[13] Yi, Z., & Murio, D. A., Identification of source terms in 2-D IHCP. Computers & 
Mathematics with Applications, 47(10-11), 1517-1533, 2004. 
[14] DW. Pepper & JC Heinrich, The finite element method - basic concepts and applications, 
Taylor & Francis, Group, pp. 240, 1992. 
[15] L. Edsberg, Introduction to computation and modeling for differential equations, 
Wiley-Interscience, pp. 256, 2008. 
[16] WBJ Zimmerman, Multiphysics modeling with finite element methods, World Scientific 
Publishing, pp. 432, 2006. 
[17] RW Pryor, Multiphysics modeling using Comsol v.4 - A first principles approach, 
Mercury Learn Inform, 2012. 
[18] Lefèvre, F., & Le Niliot, C., A boundary element inverse formulation for multiple point 
heat sources estimation in a diffusive system: Application to a 2D experiment. Inverse 
Problems in Engineering, 10(6), 539-557, 2002. 
Tác giả chịu trách nhiệm bài viết: 
TS. Trần Thanh Phong 
Trường Đại học Tiền Giang 
Email: tranthanhphong@tgu.edu.vn