Xem mẫu

  1. 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
  2. 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 Wm2 đượ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
  3. 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 , ti1  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
  4. 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
  5. 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
  6. 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 Cn1,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 n1 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
  7. 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
  8. 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) ddt tf si (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
  9. 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 Ti1 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
  10. 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 n1  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 si (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 dk 1    0 cập nhật giá trị k  k  1 và lặp lại Bước 2. k 1 81
  11. 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 Jm3 K 1 2,34.104 h Hệ số truyền nhiệt tự nhiên Wm2 K 1 15 tf Thời gian thực nghiệm s 1500  Độ dẫn nhiệt Wm1K 1 160 max Mật độ dòng nhiệt cực đại Wm2 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
  12. 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
  13. 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
  14. 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
  15. 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