- Trang Chủ
- Vật lý
- Ứng dụng phương pháp gradient phối ngẫu hiệu chỉnh để nhận dạng thông số của hệ thống
Xem mẫu
- TRƯỜNG ĐẠI HỌC SÀI GÒN SAIGON UNIVERSITY
TẠP CHÍ KHOA HỌC SCIENTIFIC JOURNAL
ĐẠI HỌC SÀI GÒN OF SAIGON UNIVERSITY
Số 71 (05/2020) No. 71 (05/2020)
Email: tcdhsg@sgu.edu.vn ; Website: http://sj.sgu.edu.vn/
ỨNG DỤNG PHƯƠNG PHÁP GRADIENT PHỐI NGẪU
HIỆU CHỈNH ĐỂ NHẬN DẠNG THÔNG SỐ CỦA HỆ THỐNG
Applying modified conjugate gradient method for identifiying parameter of system
TS. Trần Thanh Phong(1), ThS. Nguyễn Hoàng Phương(2)
(1),(2)Trường Đại học Tiền Giang
TÓM TẮT
Bài báo này nghiên cứu vấn đề giải bài toán ngược liên quan đến việc nhận dạng giá trị các thông số bất
định của hệ thống được mô tả bởi phương trình đạo hàm riêng. Phương pháp nhận dạng trực tuyến hàm
mật độ nhiệt lượng của 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 phương pháp cửa sổ trượt có dự báo. Theo đó, một nhóm cảm biến cố định được đặ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. Giải thuật lựa chọn
cảm biến được xây dựng để tìm ra vị trí của ba cảm biến hữu ích nhất cho việc nhận dạng trong khoảng
thời gian làm việc của cửa sổ trượt nhằm giảm thời gian tính toán. Giá trị nhiệt độ đo đạc bị tác động
bởi các nhiễu tuân theo hàm phân phối chuẩn Gauss. Kết quả nghiên cứu cho thấy, phương pháp được
đề xuất có khả năng nhận dạng tốt hàm mật độ nhiệt lượng với độ trễ thấp với 241s và sai số là 0,5K
đáp ứng yêu cầu đặt ra. Hiệu quả của phương pháp này được chứng minh bằng các kết quả số thông qua
việc mô phỏng.
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
This paper focuses on an ill-posed inverse problem with unknown parameters of a system, which can be
described by partial differential equations (PDE). An online identification method of the density of a
mobile heating source is proposed. This approach was based on the algorithm of the predictive online
conjugate gradient method (PCGM) with automatically adaptive sliding window size, which is based on
the iterative regularization methods for the general case. Assuming that the working area is limited, a set
of fixed sensors was placed on the domain to acquire the evolutions of the temperatures in a spatial-
temporal manner. A method was built on the order to select three sensors that are the most effective for
the identification procedure by decreasing the computational time. The measurements are disturbed
according to a realistic Gaussian noise. The results of this research show that the proposed method can
be able to identify the flux density of heat source with the low delay time about 241 seconds and the
good response error about 0.5K. The effectiveness of the proposed method were illustrated by the
numerical results.
Keywords: conjugate gradient method, heat source, identification, inverse problems, partial differential
equations
Email: tranthanhphong@tgu.edu.vn
72
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
1. Đặt vấn đề thuật mới của phương pháp nhận dạng hệ
Quá trình truyền nhiệt của nguồn nhiệt thống một cách đầy đủ theo hướng của
di động trong bối cảnh của các vấn đề phi Hadamard kết hợp phương pháp cửa sổ
tuyến là một ví dụ có liên quan đến các trượt trực tuyến và giải thuật lựa chọn cảm
hiện tượng vật lý có thể mô hình hóa bởi biến. Giải thuật nhận dạng này được xây
phương trình vi phân đạo hàm riêng. dựng dựa trên phương pháp gradient phối
Phương pháp này được đề xuất nhằm kiểm ngẫu (conjugate gradient method) bao gồm
chứng hiệu quả của phương pháp ước ba vấn đề chính cần giải quyết: vấn đề trực
lượng dựa trên phương pháp gradient phối tiếp (direct problem), vấn đề bổ trợ (adjoint
ngẫu (CGM) [1], [2] với việc sử dụng một problem) và vấn đề độ nhạy (sensitivity
nhóm cảm biến hữu dụng. Nhóm cảm biến problem) [6], [8], [9] bằng mô hình toán
này được xác định dựa trên vị trí của nhóm của bài toán ngược dựa trên phương pháp
cảm biến cố định được thiết lập tại không CGM, phương pháp cửa sổ trượt kết hợp
gian di chuyển của nguồn nhiệt thông qua với giải thuật xác định cảm biến tối ưu.
một giải thuật dựa trên mối quan hệ giữa vị 2. Nghiên cứu
trí của nguồn nhiệt và vị trí của cảm biến.
2.1. Mô tả hệ thống
Việc này nhằm giúp loại bỏ vị trí của các
Để xây dựng mô hình thí nghiệm cho
cảm biến không hữu ích để giảm thời gian
phương pháp nhận dạng thông số của hệ
tính toán của hệ thống.
thống được đề xuất trong bài viết này, giả
Trong một vài nghiên cứu trước đây,
sử có một nguồn nhiệt S di chuyển trên bề
vấn đề giải bài toán nghịch của các hiện
tượng vật lý được mô tả bởi phương trình mặt của một tấm kim loại bằng nhôm hình
đạo hàm riêng, cụ thể là quá trình truyền vuông 3 có kích thước cạnh bên L
nhiệt được đề xuất một cách không đầy đủ và độ dày e được mô tả như trong Hình 1.
bởi Hadamard với việc giải bài toán bằng Giới hạn biên của miền làm việc được ký
phương pháp lặp [3]. Phương pháp hiệu 2 [7], [10], [11].
gradient phối ngẫu cũng được sử dụng để Biến số không gian của hệ thống
tối thiểu hóa sai lệch bằng kỹ thuật lặp và
x, y L , L L , L được
2 2 2 2
chứng tỏ đây là phương pháp ổn định để
ước lượng các thông số [4]. Tuy nhiên, nó tính bằng mét và biến số thời gian là
chỉ có thể nhận dạng riêng lẻ một thông số t 0, t f được tính bằng giây. Tấm
của nguồn nhiệt. Hơn nữa, phương trình kim loại này được đốt nóng bởi hai nguồn
truyền nhiệt tổng quát trong không gian đa nhiệt lần lượt có các hàm mật độ thông
chiều rất phức tạp nên làm cho mô hình lượng nhiệt tương ứng là t được tính
toán của hệ thống trở nên cồng kềnh và gây
mất nhiều thời gian để xử lý [5], [6]. Để cải bằng Wm2 được giả định bằng một đĩa
thiện hiệu năng của phương pháp, một giải đồng chất di động D I (t ), rI có tâm
thuật lặp có hiệu chỉnh cũng được đề xuất I xI (t), yI (t ) và bán kính rI . Hàm phân
để 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 bố nhiệt độ của tấm kim loại x, y, t là
hệ thống cần xem xét [4], [5], [7]. hàm liên tục theo không gian và thời gian
Nghiên cứu này sẽ đề xuất một giải được tính bằng Kelvin.
73
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
L
y Quỹ đạo nguồn nhiệt
x L h
Nguồn nhiệt
z
r
y
h I x t , y t
e
x
Hình 1. Mô tả hệ thống thí nghiệm
Giả định rằng các giá trị của các thông biến thời gian và theo các tọa độ trong
số Ψ , , , c, , h, (t ), I (t ),0 của hệ không gian như sau:
(t )
thống sử dụng để xây dựng cho mô hình thí x, y , t arccotan
nghiệm là biết trước và được liệt kê trong (2)
Bảng 1 với đơn vị đo của các đại lượng
được tuân theo đơn vị đo lường trong hệ
2 2
x xI (t ) y y I (t ) rI
thống đo lường quốc tế SI (tiếng Pháp: Với, hệ số được chọn nhằm
Système International d'unités). mục đích rời rạc hóa hàm mật độ dòng
Trong đó, tấm kim loại sẽ được đốt nhiệt liên tục. Khoảng thời gian có thể
nóng bởi hai nguồn nhiệt di chuyển trên bề được chia ra thành N t đoạn
mặt (mặt phẳng xOy ) của tấm kim loại để
Nt 1
giúp ta khảo sát quá trình truyền nhiệt trên T
i 0
ti , ti1 với ti i và bước chia
bề mặt và bên trong tấm kim loại này. Quỹ
đạo di chuyển của hai nguồn nhiệt được rời rạc hóa được định nghĩa bởi
mô tả như trong Hình 2 (a). Đồng thời, các tf Nt . Để tránh mất tính tổng quát,
hàm mật độ dòng nhiệt của các nguồn được phương trình quỹ đạo của tất cả các vị trí
cho bởi hàm số có đồ thị được thể hiện như định vị bất kỳ của các nguồn nhiệt
trong Hình 2 (b). Biểu thức của hàm mật I xI (t ), yI (t ) cũng được thành lập lại
độ công suất nhiệt tổng của cả hai nguồn
dưới dạng các hàm rời rạc một cách tuyến
trên x, y, t được sử dụng để đốt nóng
tính và được viết lại bằng cách sử dụng các
tấm kim loại thực nghiệm được diễn đạt
hàm nón cơ bản s i (t ) với i 0, , Nt :
như sau:
(t ) if x, y D I (t ), rI (1)
1 t / i if ti 1 t ti
x, y , t
s (t ) 1 t / i
i
if ti t ti 1 (3)
0 otherwise
Theo một cách khác, biểu thức 0 otherwise
x, y, t còn có thể được biểu diễn một Hàm mật độ dòng nhiệt được diễn
cách liên tục và khả vi dưới dạng hàm tổng đạt lại như sau (t ) i si (t ) ( )tr s(t )
hợp của các hàm mật độ thành phần theo và phương trình quỹ đạo di chuyển của
74
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
các nguồn nhiệt cũng được diễn đạt lại yI (t ) yIi si (t ) ( yI )tr s(t ) . Với “tr” là ký
như sau xI (t ) xI s (t ) ( xI ) s(t ) ,
i i tr
hiệu ma trận chuyển vị.
4
x 10
Trajectory of source Sensor Ci 4
0.25 Real flux Flux to be estimated Initial flux
0.2
3.5
0.15
0.1
3
Flux density [W/m2]
0.05
Y [m]
0
2.5
-0.05
-0.1 2
-0.15
-0.2 1.5
-0.25
1
-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0 500 1000 1500
X [m] Time [seconde]
(a) Quỹ đạo di chuyển của nguồn nhiệt (b) Mật độ dòng nhiệt
Hình 2. Hàm mật độ dòng nhệt và quỹ đạo di chuyển của các nguồn nhiệt
2.2. Vấn đề thuận (direct problem) trong không gian và thời gian là kết quả
Nếu tất cả các thông số được biết trong nghiệm của phương trình đạo hàm riêng
bảng phụ lục, sự tiến triển của nhiệt độ bậc hai:
( x, y, t ) ( x, y, t ) 2h ( x, y, t )
c ( x , y , t ) ( x, y, t )
t e
( x, y,0) 0 ( x, y ) (4)
( x, y, t )
0 ( x, y, t )
n
a) Phân bố nhiệt độ tại t=300s (b) Phân bố nhiệt độ tại t=600s
75
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
(c) Phân bố nhiệt độ tại t=900s (d) Phân bố nhiệt độ tại t=1500s
Hình 3. Phân bố nhiệt độ trên tấm kim loại theo thời gian
Temperature evolution in time
560
525
490
455
Temperature [K]
420
385
350
315
280
0 250 500 750 1000 1250 1500
Time [s]
Hình 4. Hàm mật độ dòng nhiệt và quỹ đạo di chuyển của các nguồn nhiệt
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ố của 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 t=(300, 600, 900, 1500) 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 Cn được đặt cố định trên tấm kim loại như trong Hình
2(a) 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. Những nhiễu này
76
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
được tuân theo hàm phân phối xác suất Gauss N , 2 với giá trị trung bình 0 và độ
lệch chuẩn 1 .
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” ˆ(Cn , t ) tại vị trí
cảm biến Cn1,2,...,25 , 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:
tf
J , ( x, y, t )
1 N
Cn , t, ( x, y, t ) ˆ Cn , t dt
(5)
2
2 0 n1
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. Độ nhạy (sensitivity problem)
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, y, t x, y, t .
x, y , t x , y , t
x, y, t lim
0
Vấn đề độ nhạy của hệ thống được mô tả bởi hệ thống sau:
( x, y, t ) ( x, y, t ) 2h ( x, y, t )
c t
( x, y, t )
e
( x, y, t )
( x, y,0) 0 ( x, y )
( x, y, t ) (6)
0 ( x, y, t )
n
Trong đó, sự thay đổi của là:
( x, y, t )
( x, y, t )
(t )
(t )
arccotan
k 1 Arg min J ,
k 1
0.
(t ) k k 1
J , k 1 d
x x (t ) y y (t )
I
2
I
2
rI
hoặc
k 1
Nghiệm của vấn đề độ nhạy ( x, y, t)
Điều này suy ra độ tăng giảm
k 1
là rất hữu ích trong viêc tính toán độ tăng
k 1 k 1 được tính toán cho mỗi vòng lặp là:
giảm k 1 cho mỗi vòng lặp: k 1 d
k
t, ˆ t C , t, dt
tf N
c
k 1 k
Độ tăng giảm
k 1
làm tối thiểu tiêu Cn Cn n (7)
0 n 1
k 1
: C , t , dt
tf N
k 1 2
chuẩn J ,
c
k
n
0 n 1
Để tính toán vấn đề độ nhạy, chiều
77
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
k 1 problem)
hướng tăng giảm d
6N
t
phải là được
Nhằm mục đích tính toán gradienr
biết thông qua việc tính toán với hàm mục
J i 1,2,..., Nt cho mỗi vòng
tiêu, cụ thể là gradient của hàm mục tiêu J i
phải được tính toán thông qua vấn đề phối
ngẫu sau đây. lặp, một công thức Lagrange ( ( x, y, t ), , )
2.3.3. Vấn đề phối ngẫu (adjoint được giới thiệu như sau:
( ( x, y, t ), , ) J ( ( x, y, t ), )
( x, y, t ) 2h ( x, y, t ) 0
t
f
c ( x, y, t ) ( x, y , t ) d dt.
0
t e
Nếu ( x, y, t ) là nghiệm của phương ( , , )
( x, y, t ) 0 , và khi các
trình (4) do đó: ( x, y, t )
( ( x, y, t ), , ) J ( ( x, y, t ), ) phương trình phức hợp được kiểm chứng,
và ( ( x, y, t ), , ) J ( ( x, y, t ), ) . nên:
Độ biến thiên Lagrange có thể được ( , , )
viết lại:
( , , )
( , , ) J ( ( x, y, t ), )
( , , ) ( x, y, t )
( x, y , t )
( , , ) Từ J ( x, y, t ), ( ( x, y, t ), , ) (t )
( x, y, t ) (t )
( x, y , t )
k 1
( , , ) và với d (t ) Cn (t , ) ˆCn (t ) , độ
( x, y , t ).
( x, y , t ) biến thiên Lagrange có thể được viết lại
Với ( x, y, t ) được cố định để như sau:
tf tf
Nc
( x, y, t )
( x, y, t ), , d (t ) ( x, y, t ) D Cn d dt c ( x, y, t )d dt
0 n 1 0
t
2h ( x, y, t )
tf tf t
f
( x, y, t ) ( x, y, t )d dt ( x, y, t )d dt ( x, y, t )d dt.
0 0
e 0
e
với
Nc
Biết rằng E ( x, y, t ) d (t ) D x xCn D y yCn D ( x xCn ) D ( y yCn )
n 1
là hàm phân phối Dirac liên quan đến cảm biến Cn ( xCn , yCn ) và ( ) ( x, y, t ) , nên:
tf
( ) 2h ( )
( ), , E ( ) c ( ) ( ) d dt
0
t e
tf tf
( )
c x, y, t t f ( )d ( ) ( ) dt d dt.
0 0
e
78
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
x, y, t , ,
Từ ( x, y, t ) 0 x, y, t , sau đó nhân tố ( x, y, t ) 0 , nên
x, y, t
nó cần thiết rằng hàm phối ngẫu x, y, t là nghiệm của hệ phương trình sau:
( x, y, t ) 2h ( x, y, t )
c ( x , y , t ) E ( x , y , t ) ( x, y, t )
t e
( x, t , t f ) 0 ( x, y ) (8)
( x, y, t )
0 ( x, y, t )
n
Cuối cùng, x, y, t là kết quả của vấn đề vừa nêu và khi đó ta có:
( x, y, t ), ,
t
f
( x, y , t )
( x, y, t ), , d dt J ( x, y, t ),
0
e
Sau khi một số sự phát triển toán học không được trình bày do số lượng trang hạn chế
của bài viết này, các gradient chức năng của hàm mục tiêu có thể được xây dựng như sau:
x x (t) y y (t) r (xe, y, t) ddt
tf
si (t )
J i arccotan
2 2
0
I I I
(9)
Từ những công thức của gradient như 2.4.1. Phương pháp cửa sổ trượt
trên, chiều hướng tăng giảm sẽ được ước Trong các nghiên cứu mới đây [16],
lượng cho mỗi vòng lặp mới k 1 (với các tác giả đều quan tâm của việc thích ứng
là module chuẩn Euclidean): phương pháp chuyển đổi liên hợp để đạt
2 được kết quả nhận dạng mong muốn. Thời
k 1 J ( , k ) k gian tính toán phụ thuộc vào độ phức tạp
d J ( , k ) d (10)
2 của mô hình và số lượng dữ liệu sử dụng
J ( , k 1 )
làm đầu vào thuật toán. Để có kết quả nhận
Trong phần tiếp theo, việc áp dụng dạng trực tuyến, cần phải chọn số lượng dữ
phương pháp nhận dạng trực tuyến dựa liệu cần thiết có liên quan, tức là xác định
trên phương pháp lập hoá lặp đi lặp lại khoảng thời gian tốt nhất cho mỗi thủ tục
(CGM) xem xét một mạng lưới các cảm xác định bắt buộc (xem Hình 4). Vì vậy,
biến cố định được thực hiện theo số để xác tác giả đề xuất một thuật toán tính toán độ
định mật độ thông lượng của nguồn nhiệt. dài của khoảng thời gian cửa sổ trượt được
2.4. Phương pháp cửa sổ trượt kết hợp xác định dựa trên các giá trị của hàm mục
với giải thuật xác định cảm biến tối ưu tiêu và tiêu chí dừng của giải thuật.
79
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
Giải thuật xác định khoảng thời gian cửa sổ trượt
Bước 1: Thiết lập thông số
Độ rộng cơ bản của cửa sổ là 1 giây Ti i , i với i i 1
Bước 2: Tính toán giá trị hàm mục tiêu trên Ti tích hợp với việc xác định giá trị dự báo
của hàm mật độ trong khoảng Ti1
J ; , I
1 Nc
Cn , t ; , I ˆ Cn , t dt
2
2 Ti n 0
Bước 3: Kiểm tra điều kiện
If J ; , I J stop (với là hệ số để hiệu chỉnh hệ thống)
i i 1 và trở về Bước 2.
Else
Thoát khỏi giải thuật và thực hiện việc nhận dạng trên Ti i , i
Kết quả thực hiện giải thuật lựa chọn nhận dạng thông số hệ thống trên cửa sổ
cửa sổ trượt như trên và áp dụng giải thuật trượt nhằm giảm bớt số lượng dữ liệu đầu
này cho phương pháp nhận dạng thông số vào của giải thuật. Từ đó, giúp làm giảm
hệ thống được trình bày trong phần kết quả đáng kể độ phức tạp của giải thuật và đồng
và thảo luận. nghĩa với giảm thời gian tính toán. Việc
2.4.2. Giải thuật xác định cảm biến này được thực hiện dựa trên đề xuất một
Như trình bày ở phần đặt vấn đề, bài nhóm gồm 25 cảm biến cố định trên không
báo đề xuất một giải thuật nhằm loại bỏ các gian làm việc và dựa vào khoảng cách từ
cảm biến không hữu dụng cho quá trình nguồn nhiệt đến các vị trí theo thời gian :
L2di I x, y Ci x, y xI (t ) xi (t ) yI (t ) yi (t )
2 2
(11)
Giải thuật xác định cảm biến
Bước 1: Thiết lập thông số
Độ rộng cửa sổ trượt Ti i , i
Danh sách cảm biến tiềm năng Cn, lựa chọn Cs.
Bước 2: Tính toán giá trị khoảng cách từ nguồn đến từng cảm biết trên Ti
L2di I x, y Ci x, y
Bước 3: Kiểm tra điều kiện
If k 3 && L2dk min L2di
80
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
Cs=Cs&Ck và Cn=Cn\Ck
→ lặp lại Bước 2.
Else
Kết thúc giải thuật.
Hai giải thuật nêu trên được vận dụng để xây dựng giải thuật CGM để nhận dạng
thông số trực tuyến, được mô tả như sau:
Giải thuật CGM cho việc nhận dạng thông số trực tuyến
Bước 1 : Thiết lập
Chọn vector trạng thái đầu: (t )
k 0
Xác định giá trị cửa sổ trượt Ti [ i , i ] và vị trí cảm biến hữu dụng Cn
Bước 2 : Thực thi giải thuật lặp dựa trên phương pháp CGM
Giải vấn đề thuận và tính giá trị hàm mục tiêu
Cập nhật giá trị đo (Cn , t , k (t )) của cảm biến nth.
tf
Tính toán giá trị hàm mục tiêu
1 N
Cn , t ˆ Cn , t dt
2
J , ( x, y , t )
2 0 n1
If ( J ( k (t )) J stop || k Nmax || tTi Ti )
→ Kết thúc giải thuật lặp trên cửa sổ Ti
Giải vấn đề phối ngẫu và tính toán giá trị gradient
Giải phương trình (7) để tìm giá trị ( x, y, t )
Tính giá trị gradient của hàm mật độ công suất
tf
si (t ) ( x, y, t )
J i arccotan x xI (t ) y yI (t ) rI d dt
2 2
0
e
k 1 k
Tính hướng tăng giảm d J ( , k ) k d
Giải vấn đề độ nhạy để tính độ tăng giảm
Giải phương trình (6) để tìm giá trị x, y, t trong hướng d k 1 .
Tính độ tăng giảm
k 1
Argmin J (
).
Tính giá trị mới của thông số cho vòng lặp mới
Giá trị mới của mật độ : k 1 k k 1 dk 1
0 cập nhật giá trị k k 1 và lặp lại Bước 2.
k 1
81
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
2.5. Kết quả và thảo luận này di chuyển trên một tấm nhôm hình
Để kiểm nghiệm hiệu quả của phương vuông 2 , kích thước các cạnh
pháp nhận dạng thông số của hệ thống L 1m và độ dày e 2mm . Thời gian
được mô tả bởi phương trình đạo hàm thực nghiệm t f 1500s với thời gian rời
riêng (trường hợp nhận dạng hàm mật độ
dòng nhiệt). Tác giả tiến hành xây dựng rạc hóa 20s (hay suy ra, có 76 hệ số
mô hình thí nghiệm với một nguồn nhiệt có cần nhận dạng). Các thông số thiết lập
hàm mật độ ϕ(t) và quỹ đạo di chuyển hệ thống thí nghiệm được cho bởi Bảng 1
I(x(t),y(t)) (Hình 1 và Hình 2). Nguồn nhiệt như sau.
Bảng 1. Giá trị các thông số thiết lập cho hệ thống
Ký hiệu Định nghĩa tên gọi Đơn vị Giá trị
c Nhiệt lượng riêng Jm3 K 1 2,34.104
h Hệ số truyền nhiệt tự nhiên Wm2 K 1 15
tf Thời gian thực nghiệm s 1500
Độ dẫn nhiệt Wm1K 1 160
max Mật độ dòng nhiệt cực đại Wm2 3.104
0 Nhiệt độ ban đầu K 291
L Chiều dài/chiều rộng m 1
e Độ dày tấm kim loại m 2.10-3
r Bán kính nguồn nhiệt m 6.10-2
n Vector đơn vị (pháp tuyến hướng ra ngoài ) - -
Để thiết lập điều kiện đầu cho giải thuật được mô tả trong Hình 5. Kết quả nhận dạng
tại bước k=0, giả thuyết rằng giá trị rời rạc trong cửa sổ Ti [81, 222] s với đường màu
ban đầu của hàm mật độ dòng nhiệt xanh là giá trị thật, đường màu đỏ là giá trị
ϕk=0=1,5.103W/m2. Kết quả nhận dạng sử nhận dạng cho bởi CGM (giá trị t222s).
lặp kết hợp với phương pháp cửa sổ trượt
82
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
Hình 5. Phân bố nhiệt độ trên tấm kim loại theo thời gian
Để đánh giá kết quả của phương pháp trung bình bình phương độ lệch giữa nhiệt
đề xuất, tác giả sử dụng các tiêu chí chính độ đo đạc được và nhiệt độ từ mô hình, và
như thời gian tính toán tidentification là thời độ trễ tdelay của quá trình nhận dạng hệ
gian cần thiết để xác định được tất cả 76 thống. Hiệu quả của phương pháp CGM
hệ số; trung bình sai số residus là trung kết hợp với giải thuật lựa chọn cảm biến
bình độ lệch giữa nhiệt độ đo đạc được hữu dụng so với không ứng dụng giải thuật
và nhiệt độ từ mô hình dựa trên giá trị lựa chọn cảm biến được thể hiện trong
nhận dạng; residus độ lệch chuẩn sai số là Bảng 2 như sau.
Bảng 2. Bảng tổng hợp kết quả nhận dạng hệ thống
Giá trị
Chỉ tiêu Ký hiệu, công thức (*)
CGM1 CGM2(*)
Thời gian tính (s) tidentification 1741 5072
ˆ(t) ˆ(Cn , t) dt
Nc
1
Trung bình sai số (K) residus 0,5012 0,5115
Nc n 1 t f
ˆ(t) ˆ(Cn , t) dt
Nc
1 2
Độ lệch chuẩn (K) residus 0,9798 1,1054
Nc n 1 t f
Độ trễ (s) tdelay 309,14 2640,05
(*)
Chú thích: CGM1 có giải thuật lựa chọn cảm biến, CGM2 không giải thuật lựa chọn cảm biến.
83
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
Từ số liệu thống kê trong Bảng 2 cho báo đáng tin cậy và hiệu quả. Hơn nữa, giải
thấy, giá trị của trung bình sai số và độ lệch thuật lựa chọn cảm biến kết hợp với
chuẩn của hai giải thuật nhận dạng đều tiệm phương pháp CGM thể hiện ưu điểm trong
cận với giá trị trung bình và độ lệch chuẩn nhận dạng thông số của hệ thống với thời
hàm phân phối xác suất Gauss. Điều này gian tính toán là 1741s so với 5072s (giảm
chứng minh phương pháp nhận dạng trực 2331s, tỷ lệ giảm 155,4%) và độ trễ trung
tuyến thông số hệ thống đề xuất dựa trên bình của quá trình nhận dạng là 309s so với
phương pháp gradient phối ngẫu hiệu chỉnh 2640s. Kết quả nhận dạng hàm mật độ dòng
kết hợp với phương pháp cửa sổ trượt có dự nhiệt được thể hiện trong Hình 6.
4 (t) [1342:1500] window of 158 s
x 10
4
Real density
Estimated density
3.5
3
Flux density [W/m2]
2.5
2
1.5
1
0 250 500 750 1000 1250 1500
Time [s]
Hình 6. Hàm mật độ dòng nhệt và quỹ đạo di chuyển của các nguồn nhiệt
Thời gian nhận dạng hàm mật độ dòng kết quả nhận dạng nhận được 241 giây sau
nhiệt của phương pháp gradient phối ngẫu khi quá trình thực nghiệm kết thúc. Giá trị
chỉnh sửa có sử dụng giải thuật là hàm trễ trong quá trình nhận dạng thông số
tidentification 1741s , trong khi thời gian thực của hàm mật độ dòng nhiệt như Hình 7 với
giá trị lớn nhất là 551s và giá trị nhỏ nhất
nghiệm là t f 1500s . Hay nói cách khác, là 96s.
84
- TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN
600
Delay time
500
400
Delay [s]
300
200
100
0
0 250 500 750 1000 1250 1500
Time [s]
Hình 7. Hàm mật độ dòng nhệt và quỹ đạo di chuyển của các nguồn nhiệt
3. Kết luận và hướng phát triển và giải thuật lựa chọn các cảm biến hữu
Trong bài báo này, vấn đề giải bài toán dụng. Kết quả nghiên cứu cho thấy phương
ngược liên quan đến việc nhận dạng giá trị pháp được đề xuất có khả năng nhận dạng
các thông số bất định của hệ thống được mô rất tốt hàm mật độ nhiệt lượng với độ trễ
tả bởi phương trình đạo hàm riêng được thấp là 241s và sai số là 0,5K đáp ứng yêu
giới thiệu. Theo đó, phương pháp dạng trực cầu đặt ra. Kết quả này cho phép phát triển
tuyến hàm mật độ dòng nhiệt của nguồn việc xây dựng thí nghiệm để nhận dạng
nhiệt di động được đề xuất dựa trên phương thông số của hệ thống bất kỳ được mô tả
pháp gradient phối ngẫu hiệu chỉnh kết hợp bởi phương trình đạo hàm riêng bậc hai
với phương pháp lặp dựa trên cửa sổ trượt hoặc các hệ thống bậc cao hơn.
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), 2911-2919, 1991.
[2] Keith A. Woodbury, Inverse Engineering Handbook: Handbook Series for
Mechanical Engineering, Editor CRC Press, 2002.
[3] Oleg M. Alifanov, Inverse Heat Transfer Problems, Springer Science & Business
Media, 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, 25, 999-1008, 1998.
85
- SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020)
[5] S. Beddiaf, L. A, L. P, J. C, “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, 1242-
1248, 2012.
[6] S. Beddiaf, L. A, L. P, J. C, “Parametric identification of a heating mobile source in a
three dimensional geometry”, Inverse Problems in Science and Engineering, 23, 93-
111, 2015.
[7] Gillet. L.P, L. P, “Implementation of a conjugate gradient algorithm for thermal
diffusivity identification in a moving boundaries system”, Journal of Physics:
Conference Series, 135(1), 012082, 1-8, 2008.
[8] S. Beddiaf, L. A, L. P, J. C, “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), 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, 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, 1992.
[15] L. Edsberg, Introduction to computation and modeling for differential equations,
Wiley-Interscience, 2008.
[16] WBJ Zimmerman, Multiphysics modeling with finite element methods, World
Scientific Publishing, 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.
Ngày nhận bài: 24/10/2019 Biên tập xong: 15/5/2020 Duyệt đăng: 20/5/2020
86
nguon tai.lieu . vn