Ứ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

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: 72 Ứ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 TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN 73 1. Đặt vấn đề Quá trình truyền nhiệt của nguồn nhiệt di động trong bối cảnh của các vấn đề phi tuyến là một ví dụ có liên quan đến các hiện tượng vật lý có thể mô hình hóa bởi phương trình vi phân đạo hàm riêng. Phương pháp này được đề xuất nhằm kiểm chứng hiệu quả của phương pháp ước lượng dựa trên phương pháp gradient phối ngẫu (CGM) [1], [2] với việc sử dụng một nhóm cảm biến hữu dụng. Nhóm cảm biến này được xác định dựa trên vị trí của nhóm cảm biến cố định được thiết lập tại không gian di chuyển của nguồn nhiệt thông qua một giải thuật dựa trên mối quan hệ giữa vị trí của nguồn nhiệt và vị trí của cảm biến. Việc này nhằm giúp loại bỏ vị trí của các cảm biến không hữu ích để giảm thời gian tính toán của hệ thống. Trong một vài nghiên cứu trước đây, 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 đạ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 cũng đượ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]. Tuy nhiên, nó chỉ có thể nhận dạng riêng lẻ một thông số của nguồn nhiệt. Hơn nữa, phương trình truyền nhiệt tổng quá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à gây 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ũng đượ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]. Nghiên cứu này sẽ đề xuất một giải thuật mới của phương pháp nhận dạng hệ thống một cách đầy đủ theo hướng của Hadamard kết hợp phương pháp cửa sổ trượt trực tuyến và giải thuật lựa chọn cảm biến. Giải thuật nhận dạng này được xây dựng dựa trên phương pháp gradient phối ngẫu (conjugate gradient method) bao gồm ba vấn đề chính cần giải quyết: vấn đề trực tiếp (direct problem), vấn đề bổ trợ (adjoint problem) và vấn đề độ nhạy (sensitivity problem) [6], [8], [9] bằng mô hình toán của 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. 2. Nghiên cứu 2.1. Mô tả hệ thống Để xây dựng mô hình thí nghiệm cho phương pháp nhận dạng thông số của hệ thống được đề xuất trong bài viết này, giả sử có một nguồn nhiệt S di chuyển trên bề mặt của một tấm kim loại bằng nhôm hình vuông 3 có kích thước cạnh bên L và độ dày e được mô tả như trong Hình 1. Giới hạn biên của miền làm việc được ký hiệu 2 [7], [10], [11]. Biến số không gian của hệ thống  , , ,2 2 2 2 L L L Lx y               được tính bằng mét và biến số thời gian là 0, ft t     được tính bằng giây. Tấm kim loại này được đốt nóng bởi hai nguồn nhiệt lần lượt có các hàm mật độ thông lượng nhiệt tương ứng là  t được tính bằng 2Wm được giả định bằng một đĩa đồng chất di động  ( ), ID I t r có tâm  (t), ( )I II x y 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. SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020) 74 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 số  0Ψ , , , , , , ( ), ( ),c h t I t     của hệ thống sử dụng để xây dựng cho mô hình thí nghiệm là biết trước và được liệt kê trong 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ệ thống đo lường quốc tế SI (tiếng Pháp: Système International d'unités). Trong đó, tấm kim loại sẽ được đốt nóng bởi hai nguồn nhiệt di chuyển trên bề mặt (mặt phẳng xOy ) của tấm kim loại để giúp ta khảo sát quá trình truyền nhiệt trên 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 mô tả như trong Hình 2 (a). Đồng thời, các hàm mật độ dòng nhiệt của các nguồn được cho bởi hàm số có đồ thị được thể hiện như trong Hình 2 (b). Biểu thức của hàm mật độ công suất nhiệt tổng của cả hai nguồn trên  , ,x y t được sử dụng để đốt nóng tấm kim loại thực nghiệm được diễn đạt như sau:      ( ) , ( ), , , 0 otherwise It if x y D I t r x y t       (1) 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:       2 2 ( ) , , arccotan ( ) ( )I I I t x y t x x t y y t r              (2) Với, hệ số   đượ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 tN đoạn   1 1 0 , tN 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 f tt N  . Để tránh 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  ( ), ( )I II 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, , ti N : 1 1 1 / ( ) 1 / 0 othe rwise i i i i i t i if t t t s t t i if t t t                 (3) Hàm mật độ dòng nhiệt được diễn đạt lại như sau ( ) ( ) ( ) ( )i i trt s t s t    và phương trình quỹ đạo di chuyển của Nguồn nhiệt Quỹ đạo nguồn nhiệt     I x t y t, x y z h e L L h x y r TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN 75 các nguồn nhiệt cũng được diễn đạt lại như sau ( ) ( ) ( ) ( ) i i tr I I Ix t x s t x s t  , ( ) ( ) ( ) ( )i i trI I Iy t y s t y s t  . Với “tr” là ký hiệu ma trận chuyển vị. -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 X [m] Y [ m ] Trajectory of source Sensor Ci 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 (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) 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 ( , , ) ( , , ) 2 ( , , ) ( , , ) ( , , ) ( , ,0) ( , ) ( , , ) 0 ( , , ) x y t x y t h x y t c x y t x y t t e x y x y x y t x y t n                              (4) a) Phân bố nhiệt độ tại t=300s (b) Phân bố nhiệt độ tại t=600s SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020) 76 (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 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 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 TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN 77 được 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  . 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” ˆ( , )nC t tại vị trí cảm biến 1,2,...,25nC  , 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 ft N n n n J x y t C t x y t C t dt        (5) 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 x tt y    .       0 , , , , , , lim x y t x y t x y t               Vấn đề độ nhạy của hệ thống được mô tả bởi hệ thống sau: ( , , ) ( , , ) 2 ( , , ) ( , , ) ( , , ) ( , ,0) 0 ( , ) ( , , ) 0 ( , , ) x y t x y t h x y t c x y t x y t t e x y x y x y t x y t n                              (6) Trong đó, sự thay đổi của là:     2 2 ( , , ) ( ) ( , , ) ( ) arccotan ( ) ( ) ( )I I I x y t t x y t t t x x t y y t r                    Nghiệm của vấn đề độ nhạy ( , , )x y t là rất hữu ích trong viêc tính toán độ tăng giảm 1k  cho mỗi vòng lặp: 1 1 1 k k k k d      Độ tăng giảm 1k  làm tối thiểu tiêu chuẩn  1, kJ   :  11 Arg m n ,ik kJ       hoặc  1 1 1 0. , k k k k dJ            Điều này suy ra độ tăng giảm 1k  được tính toán cho mỗi vòng lặp là:           11 0 2 10 1 , , , , ˆ, f f c n n c t k n N k C C nk t k n n N C t dt C t d t t t                (7) Để tính toán vấn đề độ nhạy, chiều SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020) 78 hướng tăng giảm 1k d  6 tN phải là được biết thông qua việc tính toán với hàm mục tiêu, cụ thể là gradient của hàm mục tiêu phải được tính toán thông qua vấn đề phối ngẫu sau đây. 2.3.3. Vấn đề phối ngẫu (adjoint problem) Nhằm mục đích tính toán gradienr i J J          1,2,..., ti N  cho mỗi vòng lặp, một công thức Lagrange ( ( , , ), , )x y t  được giới thiệu như sau:  0 0 ( ( , , ), , ) ( ( , , ), ) 2 ( , , )( , , ) ( , , ) ( . ) , , ft x y t J x y t h x y tx y t c x y t x y t d dt t e                           Nếu ( , , )x y t là nghiệm của phương trình (4) do đó: ( ( , , ), , ) ( ( , , ), )x y t J x y t     và ( ( , , ), , ) ( ( , , ), )x y t J x y t       . Độ biến thiên Lagrange có thể được viết lại: ( , , ) ( , , ) ( , , ) ( , , ) ( , , ) ( , , ) ( , , ) ( , , ) ( , , ). ( , , ) x y t x y t x y t x y t x y t x y t                             Với ( , , )x y t được cố định để ( , , ) ( , , ) 0 ( , , ) x y t x y t         , và khi các phương trình phức hợp được kiểm chứng, nên: ( , , ) ( , , ) ( ( , , ), )J x y t                 Từ   ( ( , , ), , ) ( , , ), ( ) ( ) x y t J x y t t t            và với 1 ˆ( ) ( , ) ( ) n n k C Cd t t t       , độ biến thiên Lagrange có thể được viết lại như sau:       10 0 0 0 0 ( , , ) ( , , ), , ( ) ( , , ) ( , , ) 2 ( , , ) ( , , ) ( , , ) ( , , ) ( , , ) . f f c f f f t tN D n n t t t x y t x y t d t x y t C d dt c x y t d dt t h x y t x y t x y t d dt x y t d dt x y t d dt e e                                               Biết rằng     1 ( , , ) ( ) c n n N D C D C n E x y t d t x x y y       với ( ) ( )n nD C D Cx x y y   là hàm phân phối Dirac liên quan đến cảm biến ( , ) n nn C C C x y và ( ) ( , , )x y t  , nên:       0 0 0 ( ) 2 ( ) ( ), , ( ) ( ) ( ) ( ) , , ( ) ( ) ( . ) f f f t t t f h E c d dt t e c x y t t d dt d dt e                                                         TRẦN THANH PHONG - NGUYỄN HOÀNG PHƯƠNG TẠP CHÍ KHOA HỌC ĐẠI HỌC SÀI GÒN 79 Từ        , , , , ( , , ) 0 , , , , x y t x y t x y t x y t           , sau đó nhân tố ( , , ) 0x y t  , nên 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: ( , , ) 2 ( , , ) ( , , ) ( , , ) ( , , ) ( , , ) 0 ( , ) ( , , ) 0 ( , , ) f x y t h x y t c x y t E x y t x y t t e x t t x y x y t x y t n                            (8) Cuối cùng,  , ,x y t là kết quả của vấn đề vừa nêu và khi đó ta có:       0 ( , , ), , ( , , ) ( , , ), , ( , , ), ftx y t x y t x y t d dt J x y t 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:     2 2 0 ( ) ( , , ) arccotan ( ) ( ) f i t i I I I s t x y t J x x t y y t r d dt e                     (9) Từ những công thức của gradient như trên, chiều hướng tăng giảm sẽ được ước lượng cho mỗi vòng lặp mới 1k  (với  là module chuẩn Euclidean): 2 1 2 1 ( , ) ( , ) ( , ) k k k k k J d J d J              (10) Trong phần tiếp theo, việc áp dụng phương pháp nhận dạng trực tuyến dựa trên phương pháp lập hoá lặp đi lặp lại (CGM) xem xét một mạng lưới các cảm biến cố định được thực hiện theo số để xác định mật độ thông lượng của nguồn nhiệt. 2.4. 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 2.4.1. Phương pháp cửa sổ trượt Trong các nghiên cứu mới đây [16], các tác giả đều quan tâm của việc thích ứng phương pháp chuyển đổi liên hợp để đạt được kết quả nhận dạng mong muốn. Thời gian tính toán phụ thuộc vào độ phức tạp của mô hình và số lượng dữ liệu sử dụng làm đầu vào thuật toán. Để có kết quả nhận dạng trực tuyến, cần phải chọn số lượng dữ liệu cần thiết có liên quan, tức là xác định khoảng thời gian tốt nhất cho mỗi thủ tục xác định bắt buộc (xem Hình 4). Vì vậy, tác giả đề xuất một thuật toán tính toán độ dài của khoảng thời gian cửa sổ trượt được xác định dựa trên các giá trị của hàm mục tiêu và tiêu chí dừng của giải thuật. SCIENTIFIC JOURNAL OF SAIGON UNIVERSITY No. 71 (05/2020) 80 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 ,i i i      T với 1i i     Bước 2: Tính toán giá trị hàm mục tiêu trên iT 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 1iT        2 0 1 ˆ; , , ; , , 2 c i N n n n J I C t I C t dt        T Bước 3: Kiểm tra điều kiện If  ; , stopJ I J   (với  là hệ số để hiệu chỉnh hệ thống) 1i i     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 ,i i i      T Kết quả thực hiện giải thuật lựa chọn cửa sổ trượt như trên và áp dụng giải thuật này cho phương pháp nhận dạng thông số hệ thống được trình bày trong phần kết quả và thảo luận. 2.4.2. Giải thuật xác định cảm biến Như trình bày ở phần đặt vấn đề, bài báo đề xuất một giải thuật nhằm loại bỏ các cảm biến không hữu dụng cho quá trình nhận dạng thông số hệ thống trên cửa sổ trượt nhằm giảm bớt số lượng dữ liệu đầu vào của giải thuật. Từ đó, giúp làm giảm đáng kể độ phức tạp của giải thuật và đồng nghĩa với giảm thời gian tính toán. Việc này được thực hiện dựa trên đề xuất một nhóm gồm 25 cảm biến cố định trên không gian làm việc và dựa vào khoảng cách từ nguồn nhiệt đến các vị trí theo thời gian :         2 22 , , ( ) ( ) ( ) ( ) id i I i I i L I x y C x y x t x t y t y t      (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 ,i i i      T  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 iT    2 , ,