Nghiên cứu giải bài toán điều khiển tối ưu sử dụng phương pháp tựa theo dãy và giải thuật tính vi phân tự động ad cho hệ thống động học phi tuyến có ràng buộc - Ngô Phương Thanh

Ngô Phương Thanh và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 87(1): 9 - 13 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên 9 NGHIÊN CỨU GIẢI BÀI TOÁN ĐIỀU KHIỂN TỐI ƢU SỬ DỤNG PHƢƠNG PHÁP TỰA THEO DÃY VÀ GIẢI THUẬT TÍNH VI PHÂN TỰ ĐỘNG AD CHO HỆ THỐNG ĐỘNG HỌC PHI TUYẾN CÓ RÀNG BUỘC Ngô Phƣơng Thanh1*, Vũ Quốc Đông2 1Trường ĐH Kỹ thuật Công nghiệp – ĐHTN, 2Đại học Kỹ thuật Tổng hợp Ilmenau TÓM TẮT Bài báo đề cập đến bài toán tối ƣu động có ràng buộc ứng dụng phƣơng pháp Tựa theo dãy (Quasi- Sequential Approach) và phƣơng pháp Rời rạc hóa Collocation trực giao. Ƣu điểm chính của phƣơng pháp này dựa trên việc làm giảm kích thƣớc, hay chính là việc làm giảm số lƣợng biến trực tiếp của bài toán tối ƣu phi tuyến. Điều này cho phép thuật toán đƣợc ứng dụng với những bài toán lớn, phức tạp, đòi hỏi khối lƣợng và thời gian tính toán lớn. Mục tiêu của nghiên cứu này là phát triển một gói phần mềm đƣợc viết trên môi trƣờng C và C++ ứng dụng trong bài toán tối ƣu động có sử dụng giải thuật tính vi phân tự động AD (Automatic Differentiation). Từ khóa: tối ưu động, phương trình vi phân-đại số, phương pháp Tựa theo dãy, tối ưu phi tuyến, collocation trực giao, vi phân tự động.  ĐẶT VẤN ĐỀ Các bài toán tối ƣu động nói chung liên quan đến một hệ thống các phƣơng trình vi phân-đại số (DAEs) mô tả các hệ thống động học rất phổ biến trong các lĩnh vực cơ khí, cơ-điện tử, điện và điện tử cũng nhƣ công nghệ hóa học. Các mô hình đại diện cho các quá trình hóa học thƣờng bao gồm một hệ thống những phƣơng trình vi phân thƣờng (ODEs) mô tả những cân bằng khối lƣợng và năng lƣợng động của hệ thống, mà trong đó các phản ứng hóa học xảy ra, cùng với các phƣơng trình đại số (AEs) thể hiện các quan hệ cân bằng nhiệt động lực học, những giá trị ở chế độ làm việc xác lập, v.v... Các hệ thống mạch điện bao gồm những phần tử cơ bản nhƣ điện trở, tụ điện và điện cảm đƣợc mô tả bằng những hệ phƣơng trình vi phân mà đƣợc tổng hợp lại bằng các định luật Kirchhoff dƣới dạng các phƣơng trình đại số. Trong các hệ thống cơ khí, những hệ phƣơng trình vi phân thƣờng đƣợc dùng để mô tả các quá trình động học của những hệ thống con và các phƣơng trình đại số đƣợc dùng để tổng hợp các ràng buộc tại các khớp nối. Nhiệm vụ của bài toán tối ƣu động là thực hiện việc tìm kiếm một luật điều khiển cho một hệ thống cho trƣớc nhằm đạt đƣợc một tiêu chí tối ƣu nhất định. Bài toán nhƣ vậy bao gồm một hàm chi phí chứa các biến trạng thái (còn gọi là biến phụ thuộc) và các biến điều khiển (còn gọi là  Tel: 0915660599 biến độc lập) cùng với một tập hợp các phƣơng trình vi phân-đại số mô tả hệ thống động học. Việc giải bài toán tối ƣu động nói trên chính là việc đi tìm quỹ đạo của các biến điều khiển nhằm giảm thiểu giá trị hàm chi phí nhƣ: tìm đƣờng đi ngắn nhất, tìm thời gian xảy ra quá trình ngắn nhất, cực tiểu hóa chi phí, cực tiểu hóa thời gian tác động, giảm giá thành sản phẩm, v.v... Một bài toán tối ƣu với hệ thống động học phi tuyến đƣợc mô tả bằng hệ phƣơng trình DAEs dƣới dạng tổng quát nhƣ sau:       , 0 min , (1. ) . . , , 0 (1. ) (1. ) (1. ) 0 (1. ) z u L U L U z u a s t F z z u b z z z c u u u d z z e         Trong đó z(t) và u(t) lần lƣợt là các biến trạng thái (phục thuộc) và biến điều khiển (độc lập). Hệ phƣơng trình (1.b) là một hệ phƣơng trình DAEs mô tả hệ thống cần xét, (1.c) và (1.d) là những ràng buộc trên quỹ đạo của biến điều khiển và biến trạng thái, và (1.e) là giá trị khởi tạo của các biến trạng thái tại thời điểm ban đầu. SƠ LƢỢC NHỮNG PHƢƠNG PHÁP GIẢI BÀI TOÁN ĐIỀU KHIỂN TỐI ƢU Những phương pháp gián tiếp Bài toán điều khiển tối ƣu có thể đƣợc thực hiện bằng cách sử dụng nguyên lý cực đại Pontryagin (điều kiện cần), hoặc bằng cách giải phƣơng trình Hamilton-Jacobi-Bellman (điều kiện đủ). Những Ngô Phương Thanh và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 87(1): 9 - 13 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên 10 phƣơng pháp này đƣợc xếp vào nhóm phƣơng pháp gián tiếp. Phƣơng pháp này có thể đƣợc sử dụng để các bài toán biên trị (còn gọi là bài toán bờ) với những bƣớc nhảy và chuyển mạch trạng thái mà về mặt số học thƣờng rất khó giải do nó có bán kính hội tụ địa phƣơng rất nhỏ. Việc hình thành điều kiện tối ƣu cho bài toán thƣờng rất khó khăn và đòi hỏi những kiến thức nền tảng tốt về giải tích. Mặt khác, những phƣơng pháp giải kiểu này rất tốn thời gian. Hơn nữa, lời giải của rất nhiều bài toán tối ƣu động không thể tìm đƣợc theo cách giải tích, chẳng hạn khi tồn tại các ràng buộc về quỹ đạo với các biến điều khiển và/hoặc các biến trạng thái. Những phương pháp trực tiếp Ngƣợc lại, những phƣơng pháp trực tiếp, ví dụ nhƣ phƣơng pháp Collocation, trƣớc hết thực hiện việc tƣơng đối dễ tiến hành là rời rạc hóa bài toán tối ƣu động dựa trên phƣơng pháp phần tử hữu hạn, mà thƣờng tốn ít thời gian hơn và không yêu cầu ngƣời thực hiện phải có nền tảng tốt về giải tích; rồi sau đó thực hiện việc giải bài toán tối ƣu phi tuyến hữu hạn bằng các thuật toán quy hoạch phi tuyến sẵn có. Phƣơng pháp này thƣờng cho bán kính hội tụ địa phƣơng khá rộng. Với phƣơng pháp trực tiếp, bài toán tối ƣu động đƣợc rời rạc hóa và chuyển thành bài toán tối ƣu phi tuyến. Một số phƣơng pháp sử dụng các phƣơng án rời rạc hóa các phƣơng trình vi phân bằng những thuật toán nhƣ Euler, hình thang hay Runge-Kutta, bằng cách xác định một lƣới các điểm N trong khoảng thời gian 0[ , ]ft t , 0 1 2... N ft t t t t    . Theo những phƣơng pháp này, các phƣơng trình vi phân (1.b) trở thành các ràng buộc cân bằng của bài toán tối ƣu phi tuyến hữu hạn. Các phƣơng pháp trực tiếp khác thì sử dụng phƣơng pháp xấp xỉ các biến điều khiển và biến trạng thái sử dụng các hàm cơ sở, nhƣ hàm Splines hoặc đa thức Lagrange. Các phƣơng pháp trực tiếp mà chủ yếu là các phƣơng pháp số lại đƣợc chia thành hai nhóm, phương pháp giải Theo dãy (Sequential) và phương pháp giải Đồng thời (Simultaneous). Trong phƣơng pháp giải Theo dãy, chỉ có các biến điều khiển đƣợc rời rạc hóa thành một tập hữu hạn các biến điều khiển, và bài toán tối ƣu đƣợc thực hiện trên không gian với tập biến điều khiển này. Do đó các phƣơng pháp này thƣờng đƣợc đề cập đến nhƣ là các phƣơng pháp rời rạc hóa các biến điều khiển. Trong phƣơng pháp giải Đồng thời [2], bài toán điều khiển tối ƣu đƣợc chuyển thành bài toán tối ƣu phi tuyến hữu hạn bằng cách rời rạc hóa toàn bộ các tín hiệu điều khiển cũng nhƣ các các tín hiệu trạng thái; do đó, phƣơng pháp này đƣợc coi nhƣ phƣơng pháp rời rạc hóa đầy đủ. Do những ƣu điểm của nhóm phƣơng pháp này so với các pháp trƣớc đó, nhóm phƣơng pháp tối ƣu đồng thời trực tiếp đã đƣợc sử dụng rộng rãi trong các bài toán điều khiển tối ƣu hiện đại. Ngƣời ta cũng đã áp dụng những phƣơng pháp số đƣợc coi là khá hiệu quả để giải quyết bài toán quy hoạch phi tuyến với các ràng buộc, chẳng hạn nhƣ Quy hoạch toàn phƣơng liên tiếp với Tập tích cực (Active-Set SQP), phƣơng pháp Điểm trong (Interior-Point). Các phƣơng pháp trực tiếp sử dụng quy hoạch phi tuyến đã giải quyết khá tốt các bài toán liên quan đến các ràng buộc về quỹ đạo. Những độc giả có quan tâm có thể tìm hiểu thêm nội dung chi tiết trong tài liệu [4]. Phương pháp Tựa theo dãy - Quasi-Sequential Method (QSA) Phƣơng pháp Tựa theo dãy [3] kết hợp hai thuật toán đồng thời và tuần tự trong bài toán tối ƣu động nhằm kết hợp ƣu điểm của cả hai phƣơng pháp này bằng cách chia không gian các biến đƣợc rời rạc hóa thành hai không gian riêng rẽ: không gian các biến điều khiển và không gian các biến trạng thái. Phương pháp Collocation trực giao giải hệ phương trình vi phân-đại số DAEs Theo tài liệu [3], với việc áp dụng phƣơng pháp collocation trực giao, ở đây là phƣơng pháp Radau Collocation, một đa thức Lagrange trực giao với những điểm mút thời gian Radau, đƣợc sử dụng để xấp xỉ các biến trạng thái trong hệ phƣơng trình DAEs: , , , 0 0 0 , , ( ) ( ) NCNC NC l i l j l j l j j j i l j l i i j t t z t l t z z t t                  (2.a) Trong đó: l = 1, 2, , NL là số đoạn thời gian đƣợc chia trong khoảng thời gian cần tính toán, i,j = 1, 2, , NC là số điểm collocation trong một đoạn thời gian. Các biến trạng thái và các thành Ngô Phương Thanh và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 87(1): 9 - 13 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên 11 phần vi phân tại các điểm collocation của đoạn thời gian thứ l đƣợc tính nhƣ sau: , , , 0 ,, , 0 ( ) ( ) ( )( ) , 0,1,..., NC l l i j l i l j j NC j l il l i l j j z t l t z dl tdz t z dt dt i j NC          (2.b) Rời rạc hóa-Phương pháp Đồng thời Các biến điều khiển đƣợc rời rạc hóa thành từng đoạn hằng số u trong từng đoạn thời gian cùng với các biến trạng thái [2]. Bài toán (1) sau đó đƣợc chuyển thành bài toán tối ƣu phi tuyến thông thƣờng:     , min , (3. ) . . , 0 (3. ) (3. ) (3. ) m n mz u L U L U f z u a s t c z u b z z z c u u u d             Rời rạc hóa-Phương pháp Tựa theo dãy QSA Theo phƣơng pháp Tựa theo dãy, hệ thống các phƣơng trình (3.b) đƣợc loại bỏ, thay vào đó là một lớp bài toán con thực hiện việc mô phỏng, trong đó các biến trạng thái đƣợc tính toán dựa vào các biến điều khiển theo phƣơng trình mô tả hệ thống đã đƣợc rời rạc hóa thành hệ phƣơng trình đại số phi tuyến và giải theo phƣơng pháp Newton-Raphson, và bài toán (3) chuyển thành:     min ( ) : (4. ) . . ( ) : (4. ) (4. ) n m n u R n m m L U n m L U f z u ,u f R R a s t z z u z z u R R b u u u u R c                 Những chi tiết về phƣơng pháp Tựa theo dãy QSA đƣợc trình bày trong tài liệu [3]. Phương pháp tối ưu Điểm trong với bài toán quy hoạch phi tuyến Bài toán (4) chỉ chứa các biến điều khiển và các ràng buộc không cân bằng. Bằng cách thêm các biến phụ s để chuyển các ràng buộc không cân bằng (4.b) thành các ràng buộc cân bằng, ta có:      min : (5. ) . . 0 (5. ) (5. ) (5. ) n m n u R m L U n m L U f z u ,u f R R a s t z u s b z s z s R c u u u u R d                  Không làm mất tính tổng quát, ta có thể viết lại bài toán (5) dƣới dạng:     min : (6. ) . . 0 : (6. ) 0 (6. ) n n x R n m n f x f R R a s t h x h R R b x x R c       Trong đó: ( ) , ( ) ( ) (7)Tx s u h x z u s    Nhằm áp dụng phƣơng pháp tối ƣu điểm trong, bằng cách sử dụng hàm chặn Lôgarit, ta có:         1 min ln (8. ) . . 0 (8. ) n n i μ x R i φ x f x μ x a s t h x b       Bài toán (8) là dạng chuẩn tắc để sử dụng phƣơng pháp tối ƣu Điểm trong. Những mô tả chi tiết về phƣơng pháp giải bài toán tối ƣu phi tuyến này đƣợc thể hiện trong tài liệu [2]. Phương pháp Newton-Raphson giải hệ phương trình đại số phi tuyến Trong phƣơng pháp Tựa theo dãy [3], hệ phƣơng trình DAEs mô tả hệ thống (1.b) sau khi áp dụng phƣơng pháp Collocation trực giao đƣợc chuyển thành hệ phƣơng trình đại số phi tuyến tại mỗi đoạn collocation nhƣ hệ phƣơng trình (3.b). Tại vị trí collocation thứ th phƣơng trình hệ thống với các điểm collocation đƣợc viết gọn nhƣ sau:  ,0 , , 0 1,..., (9)c z u z NL       Trong đó ,0z là giá trị khởi tạo của phần tử z , NL là số khoảng thời gian đƣợc chia. Áp dụng phƣơng pháp Newton-Raphson giải hệ phƣơng trình đại số phi tuyến (9) với các biến z và tham số u , ta cần tính các giá trị đạo hàm với biến z : c z     (10) Tính toán độ nhạy (Sensitivities) Để có thể giải bài toán (5) đối với biến u theo phƣơng pháp Điểm trong, ta phải tính các giá trị đạo hàm (gradients) f u    và (jacobians) (Trong đó ( )z u thực chất là một hàm ẩn.) Ngoài ra, theo luật di truyền, ta còn cần tính thêm giá trị gradient ,0 ( )z u z      . Ngô Phương Thanh và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 87(1): 9 - 13 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên 12 Theo định lý về hàm ẩn, ta có thể tính các giá trị đạo hàm trên nhƣ sau: ,0 ,0 và (11) T T z zc c c c u z u z z z                                                  Cần nhắc lại (9) là một hệ phƣơng trình đại số phi tuyến nên việc tính giá trị các đạo hàm (10) và (11) theo phƣơng pháp giải tích hàm một cách chính xác là việc không đơn giản. Điều này đã làm nảy sinh yêu cầu tìm các giải thuật tính toán các giá trị đạo hàm gần đúng theo phƣơng pháp giải tích số dựa trên các phép vi phân số một cách tự động, ví dụ nhƣ phép sai phân tiến hoặc sai phân lùi [4]. Phép tính vi phân tự động-Tổng quan các công cụ phần mềm Theo tài liệu [4] thì các phƣơng pháp tính vi phân phổ biến hiện nay là: 1) Vi phân ký hiệu; 2)Vi phân hữu hạn; 3) Vi phân tự động. Hai phƣơng pháp 1) và 2) đều tồn tại nhiều nhƣợc điểm, đặc biệt khi tính toán các đạo hàm bậc cao. Phƣơng pháp 1) có tốc độ tính toán tƣơng đối chậm và khó khăn trong lập trình khi chuyển đổi một chƣơng trình máy tính thành một biểu thức. Phƣơng án 2) thì có nhƣợc điểm cơ bản là tạo ra những sai số trong quá trình rời rạc hóa và ƣớc lƣợc. Phƣơng án 3) thì không gặp phải những nhƣợc điểm trên. Do đó việc lựa chọn phƣơng án 3) trong tính toán các đạo hàm (10) và (11) là điều tất yếu, việc còn lại là lựa chọn công cụ phần mềm sẵn có thích hợp cho việc phát triển. Dựa theo [1], phần mềm mã nguồn mở cppAD [6] đƣợc chọn do nó có một số ƣu điểm nhất định so với các phần mềm còn lại nhƣ không yêu cầu nhiều bộ nhớ, cho tốc độ tính toán nhanh, dễ tìm hiểu lập trình. PHÁT TRIỂN GÓI PHẦN MỀM Phần mô phỏng hệ DAEs Phần mô phỏng hệ DAEs mà thực chất là giải hệ DAEs bằng phƣơng pháp collocation trực giao để tìm các giá trị biến trạng thái từ các giá trị của biến điều khiển đƣợc cung cấp từ bài toán tối ƣu phi tuyến. Hệ phƣơng trình đại số phi tuyến đƣợc thành lập theo phƣơng pháp Collocation sau đó sẽ đƣợc giải bằng phƣơng pháp Newton-Raphson với thông tin về đạo hàm đƣợc cung cấp từ các chƣơng trình con trong phần mềm cppAD viết trên C++. Giải bài toán tối ưu phi tuyến Bài toán tối ƣu phi tuyến với tập hợp biến là các biến điều khiển đƣợc giải bằng phƣơng pháp Điểm trong sử dụng phần mềm IPOPT viết trên C và C++. Tài liệu [3] trình bày những chi tiết về phƣơng pháp Điểm trong và phần mềm IPOPT. VÍ DỤ Để đánh giá kết quả của nghiên cứu này, một bài toán trong điều khiển tối ƣu đƣợc lấy làm ví dụ minh họa. Đó là bài toán điều khiển hệ dao động Val der Pol đã đƣợc trình bày trong tài liệu [5]:     3 2 1 2 1 2 2 1 2 2 2 3 1 2 1 min 5 (12. ) . . 1 (12. ) (12. ) (12. ) 0.3 ( ) 1.0 (12. ) 0.4 ( ) (12. ) (0) [0,1.0,0] (12. ) [0,5] (12. ) T J x a s t x x x x u b x x c x x x u d u t e x t f x g t h                    Bài toán (12) đƣợc giải trong hai trƣờng hợp có và không có ràng buộc quỹ đạo (12.f). Với trƣờng hợp không có ràng buộc (12.f), kết quả J=2.867391 có tốt hơn so với giá trị 2.86875 đƣợc thể hiện trong [5]. Với trƣờng hợp có ràng buộc (12.f), giá trị J=2.95428 cũng tốt hơn so với kết quả 2.95436 đƣợc thông báo trong [5]. Các đồ thị tín hiệu u(t) và x(t) đƣợc thể hiện trong các Hình 1 và 2 mà ta có thể thấy rõ một số cung suy biến do sự xuất hiện các ràng buộc về quỹ đạo đối với tín hiệu điều khiển u(t) và trạng thái x(t). Hình 1. Đồ thị tín hiệu điều khiển u(t) Ngô Phương Thanh và Đtg Tạp chí KHOA HỌC & CÔNG NGHỆ 87(1): 9 - 13 Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên 13 Hình 2. Đồ thị tín hiệu trạng thái x(t) KẾT LUẬN VÀ HƢỚNG PHÁT TRIỂN Bài báo đã đề cập đến phƣơng pháp giải bài toán tối ƣu với các hệ thống động học có ràng buộc theo phƣơng pháp Tựa theo dãy có sử dụng các giải thuật tính vi phân tự động AD. Một gói phần mềm cũng đã đƣợc phát triển nhằm cung cấp cho ngƣời dùng một môi trƣờng thân thiện dễ sử dụng nhờ việc tổng quát hóa bài toán tối ƣu động nhằm giảm thiểu công tác chuẩn bị và thiết lập bài toán. Đƣợc kết hợp với một phần mềm giải bài toán tối ƣu phi tuyến theo phƣơng pháp Điểm trong, gói phần mềm đƣợc viết đã chứng minh đƣợc tính hiệu quả của nó thông qua các ví dụ đã đƣợc trình bày trong bài báo. Thừa kế đƣợc những ƣu điểm của phƣơng pháp Tựa theo dãy, gói phần mềm này đƣợc trông đợi là có thể áp dụng đƣợc với những bài toán có cấu trúc phức tạp với khối lƣợng biến tƣơng đối lớn. Đây là vấn đề cần đƣợc tiến hành trong những nghiên cứu tiếp theo sâu hơn trong tƣơng lai cùng với những mở rộng về mặt lý thuyết, đặc biệt là với bài toán điều khiển tối ƣu với đối tƣợng có những quỹ đạo suy biến. TÀI LIỆU THAM KHẢO [1] Andersson, Joel; Houska, Boris; and Diehl, Moritz (2010). Towards a Computer Algebra System with Automatic Differentiation for use with Object-Oriented modelling languages. [2] A. Wächter and L. T. Biegler (2006), On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1), 25-57. [3] Hong, W.; P. Li; G. Wozny; L. T. Biegler and S. Wang (January 2006), A Quasi-Sequential Approach to Large-Scale Dynamic Optimization Problems, AIChE Journal, Volume 52, No.1, , 255-268. [4] Jorge Nocedal and Stephen J. Wright (2006), Numerical Optimization, Second Edition, Springer,. [5] V. S. Vassiliadis; R. W. H. Sargent; and C. C. Pantelides (1994). Solution of a Class of Multistage Dynamic Optimization Problems. 2. Problems with Path Constraints. Ind, Eng. Chem. Res., 33, 2123-2133. [6] A Package for Differentiation of C++ Algorithms. SUMMARY DYNAMIC OPTIMIZATION USING QUASI-SEQUENTIAL APPROACH AND AUTOMATIC DIFFERENTIATION ALGORITHM FOR CONSTRAINED NONLINEAR DYNAMICAL SYSTEMS Ngo Phuong Thanh 1 , Vu Quoc Dong 2 1College of Technology – TNU, 2Ilmenau University of Gerneral Technologgy This paper deals with dynamic optimization problems with path constraints based on Quasi-sequential Approach (QSA) and Automatic Differentiation algorithms embedded in orthogonal collocation for solving Differentiation- Algebraic Equation (DAEs) systems. This method takes advantage of reducing the dimension or decreasing the number of variables of the Nonlinear Programming (NLP) problems. This makes it possible for the QSA approach to handle time-consuming large-scale problems. The goal of this research is to develop a software environment written in C and C++ to solve dynamic optimization problems by using Automatic Differentiation algorithms. Key words: dynamic optimization, differential-algebraic equations, orthogonal collocation, quasi-sequential approach, nonlinear optimization, automatic differentiation.  Tel:0915660599