Bài giảng Các mô hình và phần mềm tối ưu hoá và ứng dụng trong nông nghiệp

Trường Đại học Nông nghiệp I Hà Nội Khoa Công nghệ thông tin PGS.TS. NGUYỄN HẢI THANH CÁC MÔ HÌNH VÀ PHẦN MỀM tối ưu hoá và ứng dụng trong nông nghiệp (bài giảng điện tử trong khuôn khổ dự án CNTT) HÀ NỘI, THÁNG 10 NĂM 2007 1 MỤC LỤC CHƯƠNG I. ỨNG DỤNG MỘT SỐ MÔ HÌNH TỐI ƯU TRONG NÔNG NGHIỆP 3 1. MÔ HÌNH QUY HOẠCH ĐƠN MỤC TIÊU 3 1.1. Mô hình tối ưu một mục tiêu (tuyến tính và phi tuyến) 3 1.2. Các ví dụ minh hoạ bài toán tối ưu một mục tiêu 4 2. MÔ HÌNH QUY HOẠCH ĐA MỤC TIÊU TUYẾN TÍNH VÀ PHI TUYẾN 8 2.1. Giới thiệu bài toán quy hoạch đa mục tiêu 8 2.2. Các khái niệm cơ bản của bài toán tối ưu đa mục tiêu 9 2.3. Các ví dụ minh hoạ bài toán quy hoạch đa mục tiêu 10 CHƯƠNG II. GIẢI BÀI TOÁN QUY HOẠCH TUYẾN TÍNH BẰNG PHƯƠNG PHÁP ĐƠN HÌNH 23 1. PHƯƠNG PHÁP ĐƠN HÌNH GIẢI BTQHTT DẠNG CHÍNH TẮC 23 1.1. Ví dụ 23 1.2. Thuật toán đơn hình cho BTQHTT dạng chính tắc 27 2. PHƯƠNG PHÁP ĐƠN HÌNH HAI PHA GIẢI BTQHTT TỔNG QUÁT 28 2.1. Ví dụ 28 2.2. Thuật toán đơn hình hai pha giải BTQHTT dạng tổng quát 30 3. GIẢI CÁC BÀI TOÁN TỐI ƯU TRÊN MICROSOFT EXCEL 32 3.1. Giải BTQHTT 32 3.2. Giải bài toán quy hoạch phi tuyến có ràng buộc tuyến tính 34 3.3. Một số ví dụ khác 36 4. GIẢI BTQHTT TRONG LINGO 36 5. GIẢI BTQHTT BẰNG PHẦN MỀM QHTT 38 CHƯƠNG III. BÀI TOÁN QUY HOẠCH PHI TUYẾN 40 1. PHƯƠNG PHÁP RST2ANU GIẢI BÀI TOÁN TỐI ƯU PHI TUYẾN TOÀN CỤC HỖN HỢP NGUYÊN 40 1.1. Đặt vấn đề 40 1.2. Thuật giải tìm kiếm ngẫu nhiên có kiểm soát RST2ANU 41 1. 3. Một số nhận xét về phiên bản nâng cấp của phần mềm 43 2. MỘT SỐ VÍ DỤ ÁP DỤNG RST2ANU 44 2.1. Bài 1: Bài toán xác định tham số sàng phân loại 44 2.2. Bài 2: Bài toán xác định cơ cấu đầu tư chăn nuôi cá 46 2 3. TÍCH HỢP RST2ANU VỚI MATLAB 48 3.1. Nhập từ bàn phím 48 3.2. Nhập từ tệp 50 CHƯƠNG IV. GIẢI BÀI TOÁN QUY HOẠCH TUYẾN TÍNH ĐA MỤC TIÊU BẰNG PHƯƠNG PHÁP THOẢ DỤNG MỜ 52 1. CÁC KHÁI NIỆM CƠ BẢN 52 1.1. Phát biểu mô hình 52 1.2. Phương án tối ưu Pareto 53 1.3. Phương pháp thoả dụng mờ giải BTQHTT đa mục tiêu 54 2. GIẢI BTQHTT ĐA MỤC TIÊU BẰNG CHƯƠNG TRÌNH MÁY TÍNH MULTIOPT 59 2.1. Ví dụ 59 2.2. Bài toán quy hoạch đất xã Nhân Chính 63 2.3. Bài toán quy hoạch đất xã Trâu Quỳ 70 CHƯƠNG V. MÔ HÌNH VÀ PHẦN MỀM TỐI ƯU PHI TUYẾN ĐA MỤC TIÊU 71 1. BÀI TOÁN TỐI ƯU PHI TUYẾN TRONG MÔI TRƯỜNG MỜ / NGẪU NHIÊN 71 1.1. Phát biểu bài toán và phương pháp mức ưu tiên 71 1.2. Xử lý các ràng buộc 72 1.3. Xử lý các mục tiêu 74 1.4. Sử dụng thông tin pay-off để đoán nhận k e , jd ~ 77 1.5. Mô hình tất định tương đương của bài toán 79 1.6. Khái niệm tối ưu hoá PL-Pareto 79 2. THUẬT GIẢI TƯƠNG TÁC LẶP PRELIME VÀ MỘT SỐ ỨNG DỤNG 80 2.1. Phát biểu thuật giải 80 2.2. Bài toán Chakraborty 81 2.3. Bài toán xác định cơ cấu đầu tư cho các hộ chăn nuôi cá 87 2.4. Bài toán quy hoạch sử dụng đất trên địa bàn huyện Trùng Khánh 88 CHƯƠNG VI. KẾT LUẬN VÀ ĐỀ XUẤT MỘT SỐ HƯỚNG NGHIÊN CỨU 92 1. ÁP DỤNG CÁC MÔ HÌNH TỐI ƯU TRONG NÔNG NGHIỆP 92 2. NGHIÊN CỨU ÁP DỤNG VÀ ĐỀ XUẤT CÁC PHƯƠNG PHÁP TỐI ƯU 92 3. XÂY DỰNG CÁC PHẦN MỀM TỐI ƯU 93 4. XÂY DỰNG HỆ HỖ TRỢ RA QUYẾT ĐỊNH CÀI ĐẶT TRÊN MẠNG MÁY TÍNH 94 DANH MỤC TÀI LIỆU THAM KHẢO 96 3 Chương I ỨNG DỤNG MỘT SỐ MÔ HÌNH TỐI ƯU TRONG NÔNG NGHIỆP Tối ưu hoá là một trong những lĩnh vực kinh điển của toán học có ảnh hưởng đến hầu hết các lĩnh vực, trong đó có nông nghiệp. Trong thực tế, việc tìm ra giải pháp tối ưu cho một vấn đề nào đó chiếm một vai trò hết sức quan trọng. Phương án tối ưu là những phương án tốt nhất, tiết kiệm chi phí, tài nguyên, sức lực mà lại cho hiệu quả cao. Có thể phát biểu mô hình (bài toán) tối ưu tổng quát như sau: F(X) ÆMax (Min) với X ∈D được gọi là miền ràng buộc. F ở đây có thể là một hàm vô hướng hay hàm véc tơ, tuyến tính hay phi tuyến. Trong trường hợp F là hàm vô hướng thì ta có mô hình quy hoạch (tối ưu) đơn mục tiêu, còn nếu F là véc tơ thì có mô hình quy hoạch (tối ưu) đa mục tiêu. X có thể là một biến đơn lẻ hay một tập hợp nhiều biến tạo thành một vectơ hay thậm chí là một hàm của nhiều biến khác. Biến có thể nhận các giá trị liên tục hay rời rạc. D là miền ràng buộc của X, thường được biểu diễn bởi các đẳng thức, bất đẳng thức, và được gọi là miền phương án khả thi hay phương án chấp nhận được. 1. MÔ HÌNH QUY HOẠCH ĐƠN MỤC TIÊU 1.1. Mô hình tối ưu một mục tiêu (tuyến tính và phi tuyến) Dạng chính tắc của bài toán tối ưu toàn cục một mục tiêu được biểu diễn như sau: Max (Min) f(X) X = (x1, x2, , xn) với: (i) gj(X) ≤ 0, j = 1, 2, , k, (ii) gj(X) = 0, j = k+1, k+2, , m, Trong các bài toán thực tế có thể bổ sung các ràng buộc dạng: (iii) ai ≤ xi ≤ bi, i = 1, 2, , n. Hàm mục tiêu f(X) và các hàm ràng buộc gj(X) với j=1,2, ,m có thể là tuyến tính hay phi tuyến.Véc tơ X có thể bao gồm các thành phần rời rạc hay liên tục hoặc là sự kết hợp giữa các thành phần rời rạc và các thành phần liên tục. Các dạng khác của bài toán tối ưu một mục tiêu đều có thể đưa về dạng chính tắc theo những quy tắc nhất định. Nếu ký hiệu D là miền các phương án (miền ràng buộc) cho bởi các ràng buộc (i), (ii) và/hoặc (iii) thì bài toán trên đây có thể viết gọn hơn như sau: f(X) → Max (Min) với X ∈ D. Lúc này, X* ∈ D được gọi là phương án tối ưu toàn cục nếu ∀ X∈D ta luôn có: f(X*) ≤ f(X). Trong trường hợp f(X*) ≤ f(X) chỉ đúng với ∀X∈D trong một lân cận của X* thì X* được gọi là phương án tối ưu địa phương. Trong trường hợp có ít nhất một trong các hàm mục tiêu hay ràng buộc là hàm phi tuyến, chúng ta có bài toán quy hoạch phi tuyến. Trong các bài toán tối ưu phi tuyến ứng dụng nói chung, và trong nông nghiệp nói riêng, lời giải tối ưu toàn cục có một ý 4 nghĩa quan trọng. Chẳng hạn trong thiết kế máy nông nghiệp, sau khi dùng phương pháp phân tích hồi qui nhiều chiều, ta thường thu được hàm mục tiêu f(X) có dạng phi tuyến. Các bài toán tối ưu toàn cục cũng có thể nảy sinh trong quy hoạch kinh tế - sinh thái vùng, hay xác định cơ cấu đất canh tác - cây trồng. Bài toán đặt ra là phải tìm được phương án tối ưu toàn cục. Có rất nhiều phương pháp giải các lớp bài toán tối ưu phi tuyến, nhưng chưa có phương pháp nào tỏ ra hữu hiệu cho mọi bài toán tối ưu phi tuyến, đặc biệt là các bài toán có các biến nhận các giá trị liên tục cũng như nguyên. Trong trường hợp tất cả các hàm mục tiêu cũng như ràng buộc đều là các hàm tuyến tính, chúng ta có BTQHTT. Trái với bài toán quy hoạch phi tuyến, BTQHTT có thể giải bằng một số phương pháp tối ưu quen biết (như phương pháp đơn hình cải biên, phương pháp hai pha, phương pháp điểm trong v.v) và được sử dụng rộng rãi trong quy hoạch sử dụng đất cũng như nhiều lĩnh vực của kinh tế và quản trị kinh doanh nông nghiệp. Đặc biệt, khi các ràng buộc đều cho ở dạng bất đẳng thức với dấu ≤ thì ta có mô hình tối ưu (quy hoạch tuyến tính) một mục tiêu sau: Min CX với ràng buộc DX ∈ , trong đó: C là véc tơ ∈ Rn D = { ∈X Rn : AX ≤ B, X ≥ 0 } với A là ma trận cấp m × n và ∈B Rm 1.2. Các ví dụ minh hoạ bài toán tối ưu một mục tiêu Bài toán quy hoạch sử dụng đất (Mô hình tối ưu tuyến tính một mục tiêu giải bài toán quy hoạch sử dụng đất trên địa bàn xã Đông Dư, huyện Gia Lâm, tỉnh Hà Nội) Chúng ta xét mô hình tối ưu một mục tiêu với mục tiêu cần cực đại hoá là hiệu quả kinh tế. Để thiết lập mô hình, trước hết chọn các biến quyết định. Dựa vào kết quả các dữ liệu đã thu được, ta chọn các biến quyết định như sau: xj với j = 1, 2, , 18 là diện tích các loại cây trồng ( theo thứ tự là: lúa xuân, lúa mùa, ngô xuân, ngô đông, ngô bao tử đông, lạc xuân, đậu xanh xuân, đậu tương đông đất chuyên màu, đậu tương đông đất ba vụ, dưa chuột xuân, dưa chuột bao tử, mướp đắng xuân, rau mùi tàu, rau gia vị, đậu cô ve đông, ớt xuân, cà chua xuân, cà chua đông), x19 là diện tích ao hồ thả cáao cá, xj với j = 20, , 24 là số đầu vật nuôi trong năm (trâu, bò, lợn, gia cầm). x24 là số công lao động thuê ngoài, x25 là lượng tiền vốn vay ngân hàng. Lúc đó chúng ta có bài toán tối ưu tuyến tính một mục tiêu với 33 ràng buộc (chưa kể điều kiện không âm của các biến) như sau: Hiệu quả kinh tế: f(X) = 4306,14 x1 + 4168,73 x2 + 3115,21 x3 + 3013,11 x4 + 4158,68 x5+ 4860,91 x6 + 4295,31 x7 + 3706,11 x8 + 3788,25 x9+ 12747,31 x10+ 12752,96 x11 + 12064,81 x12 + 79228,88 x13 + 35961,31 x14 + 10823,91 x15+ 7950,16 x16 + 7928,06 x17 +5738, 46 x18 + 11129,50 x19 + 429,00 x20 + 674,00 x21+ 219,50 x22+ 11,10 x23– 15,50 x24 – 0,12 x25 → Max Với các ràng buộc sau đây : 5 x1 ≤ 80,88; x2 ≤ 75,78; x3 ≤ 64,89; x4 ≤ 64,89; x5 ≤ 10,50; x6 ≤ 64,89; x7 ≤ 64,89; x8 ≤ 16,50; x9 ≤ 45,30; x10 ≤ 5,50; x11 ≤ 8,5; x12 ≤ 6,80; x13 ≤ 13,70; x14 ≤ 14,50; x15 ≤ 4,80; x16 ≤ 4,50; x17 ≤ 4,20; x18 ≤ 10,20; x19 ≤ 33,11; x20 ≤ 40,00; x21 ≤ 180,00; x22 ≤ 4280; x23 ≤ 18800; x5 + x9 + x11 + x13 + x18 ≤ 45,30; x3 + x6 + x7 + x10 + x 12 + x16 + x17 ≤ 64,89; x4 + x8 + x14+ x15 ≤ 64,89; x1 + x13 ≤ 80,88; x2 + x13 ≤ 75,88; 205,5x1 + 150x3 + 75,75x4 + 75x5 + 225,5x6 + 221,5x7 + 102,7x8+ 100,75x9 +360 x10 +140x11 + 385x 12 + 1833,6x13 + 1446,3x14+210,25 x15 + 410,5x16 +360,5 x17 + 176x18 + 67x19 +20x20 + 16x21 + 9x22 + 0,3x23 - x24 ≤ 226149,00; 201,5x2 + 150x3 + 75,25x4 + 102,7x8+ 100,75x9 + 140x11 + 2475,4x13 + 1446,3x14+ 210,25x15 + 176x18 + 58x19 + 16x20 + 12x21 + 7x22 + 0,2x23 - x24 ≤ 152190,00; 2871,89x1 +2691,89 x2 + 2243,62x3 + 2243,66x4 + 3630,89x5 + 4780,06x6 + 2229,11x7 + 2401,41x8+ 2326,88x9 + 16440,61x10 + 16058,39x11 +15960,61x 12 + 68494,59x13 + 23146,11x14+ 13676,26x15 +6061,76x16 + 11083,11x17 + 10391,89x18 + 18058x19 + 1223x20 + 1098,5x21 + 624,5x22 + 12x23 - x24 ≤ 3881500; 3,5x5 + 8x6 + 3,5x7 +4,1x8+ 3,5x9 + 4,16x10 + 3,5x11 + 4x 12 + 12,1x13 + 14,4x14+ 3,42x15 + 11,58x16 + 8x17 + 7,5x18 -3 x20 –2x21 - 0,95x22 – 0,0052x23 ≤ 0; 5,1x1 + 4,96x2 + 3,85x3 + 3,8x4 ≥ 921,25; Điều kiện không âm của các biến: ∀xj ≥ 0 ( j = 1, 2, , 25). Bằng phần mềm thương phẩm thích hợp có sẵn Lingo hay sử dụng Solver của Excel (xem chương II) có thể tìm được phương án tối ưu của bài toán trên như sau: x1=67,18, x2=62,08, x3=25,32, x4=45,59, x5=10,50, x6=3,37, x9=2,40, x10=6,50, x11=8,50, x12=6,50, x13=13,70, x14=14,50, x15=4,80, x16=4,50, x17=4,20, x18=10,20, x19=33,11, x20=40,00, x21=180, x22=4280, x23=18800, x25=230701010,78. Hiệu quả kinh tế cực đại đạt được là 4270,36. Bài toán tối ưu hoá giá trị sản xuất (Mô hình tối ưu phi tuyến một mục tiêu giải bài toán tối ưu hoá giá trị sản xuất trên một héc ta nuôi cá tại huyện Văn Giang, tỉnh Hưng Yên) Sử dụng số liệu điều tra 112 hộ nuôi cá vùng đồng trong đê thuộc 4 xã Văn Giang, Hưng Yên, chúng tôi chạy mô hình hồi quy tương quan trong Excel và nhận đư- ợc kết quả hàm sản xuất Cobb – Douglas như sau (cần cực đại hoá): Z = f(X) = 19,375 x10,236 x20,104 x30,096 x40,056 x50,056 e0,168 x6 e0,066 x7 →Max trong đó: z : Giá trị sản xuất bình quân triệu/ ha/năm (GO), x1 : Chi phí giống bình quân 1 ha 1 năm (tr/ ha), x2 : Chi phí thức ăn bình quân 1 ha 1 năm (tr/ ha), 6 x3 : Chi phí lao động bình quân 1 ha 1 năm (tr/ ha), x4 : Chi phí khấu hao và thuê đất bình quân 1 ha 1 năm (tr/ ha), x5 : Các chi phí khác bình quân 1 ha 1 năm (tr/ ha), x6 , x7: Biến giả định về hình thức nuôi, x6 = 1 đối với nuôi chuyên canh, x6 = 0 đối với nuôi tổng hợp, x7 = 1 với hình thức nuôi 1 loại cá chính kết hợp với các loại cá khác, x7 = 0 với hình thức nuôi 2 loại cá chính kết hợp với các loại cá khác. Với từng mức đầu tư/ tổng chi phí TC ta có các ràng buộc: - Với mức đầu tư dưới 40 tr đ/ ha: TC < 40 - Với mức đầu tư 40 - 50 tr đ/ ha: 40 ≤ TC < 50 - Với mức đầu tư 50 – 60 tr đ/ ha: 50 ≤ TC < 60 - Với mức đầu tư 60 - 70 tr đ/ ha: 60 ≤ TC < 70 - Với mức đầu tư trên 70 tr đ/ ha: TC ≥ 70 trong đó: x1 + x2 + x3 + x4 + x5 = TC . Với hình thức nuôi ta có: x6+ x7 = 1 (x6 , x7 chỉ nhận các giá trị 0 hoặc 1). Trên đây là bài toán tối ưu phi tuyến, với 5 biến liên tục và 2 biến nguyên. Sử dụng phần mềm RST2ANU (xemchương III) để giải bài toán tối ưu phi tuyến toàn cục hỗn hợp nguyên đã thiết lập trên đây ta có kết quả trong bảng I.1. Bảng I.1. Kết quả cơ cấu đầu tư tối ưu vùng đồng Đầu tư (tr/ha) 70 x1 35 – 45% 40 – 45% 40 – 45% 35 – 45% 35 – 40% x2 15 – 20% 17 – 25% 17 – 23% 15 – 20% 18 – 25% x3 15 – 20% 15 – 20% 15 – 20% 16- 19% 17 – 23% x4 10 – 15% 7 – 15% 8 – 15% 9 – 13% 10 – 15% X5 10 – 15% 10 – 15% 10 - 15% 9 - 15% 10 - 15% GO (tr/ ha) 106 NI (tr/ ha) - 38,1-38,3 38,3-37,5 37,5-36 - Việc thực hiện cơ cấu đầu tư tối ưu làm giá trị sản xuất (GO) cũng như thu nhập ròng (NI = GO - TC) ở từng mức đầu tư tăng lên rõ rệt so với thực tế sản xuất tại địa phương. Đặc biệt, mức đầu tư 50 tr/ha cho ta thu nhập hỗn hợp cao nhất 38,3 tr/ha, lớn hơn 8 tr/ha so với hiện tại không áp dụng cơ cấu đầu tư tối ưu cũng như hình thức nuôi thích hợp. Tại mức đầu tư này, cơ cấu đầu tư tối ưu là x1 từ 19,6 – 21,5 triệu (39,2 – 42,2%); x2 từ 8,6 - 9,8 triệu (17,2 – 19,6%); x3 từ 8,6 – 9,9 triệu ( 17,2 – 19,8%); x4 từ 4,7 – 6,4 triệu (9,4 – 12,8%); x5 từ 4,9 – 6,3 triệu (9,8 –12,6%) với hình thức nuôi chuyên canh (x6=1). 7 Kết quả áp dụng phần mềm RST2ANU (xem chương III) tại mức đầu tư 50 triệu đồng/ha cho phương án tối ưu sau: zmax=88,360733với x1=21,498072, x2=9,528987, x3=8,758034, x4=5,138906, x5=5,076000, x6=1,000000, x7=0,000000. Bài toán tối ưu thông số sàng phân loại (Mô hình tối ưu phi tuyến một mục tiêu giải quyết vấn đề tính toán một số thông số hình học và động học của cơ cấu sàng phân loại dao động) Trong ví dụ này chúng tôi chỉ xin nêu vắn tắt một ứng dụng của mô hình tối ưu phi tuyến một mục tiêu trong việc tìm nghiệm của hệ phương trình phi tuyến sau phát sinh trong việc tính toán một số thông số hình học và động học của cơ cấu sàng phân loại dao động (cần chú ý rằng nhiều phương pháp tính toán thông dụng khác của giải tích số đã tỏ ra không hiệu quả): r cosϕ1 + l cosϕ2 + l’’3 cosϕ3 + l4 cosϕ4 – xC1 = 0; r sinϕ1 + l sinϕ2 + l’’3 sinϕ3 + l4 sinϕ4 – yC1 = 0; r cosϕ1 + l cosϕ2 + l’3 cos(ϕ3 - α)+ l5 cosϕ5 – xD1 = 0; r sinϕ1 + l sinϕ2 + l’3 sin(ϕ3 - α)+ l5 sinϕ5 – yD1 = 0. Trong hệ phi tuyến trên các thông số đã biết là: r = 0,05m; l=0,30m; l’’3 = 0,15m; l’3 = 1,075m; l3 = 1,025m; l4 = 0,50m; l5 = 0,40m; xC1 = 0,365m; yC1 = 0,635m; xD1 = 1,365m; yD1 = 0,635m; α = π/8. Để sử dụng phần mềm tính toán tối ưu phi tuyến RST2ANU giải hệ phương trình phi tuyến cho ϕ = kπ/8 (k=0,, 9), trước hết chúng ta cần thiết lập cực tiểu hoá hàm mục tiêu sau: z = (r cosϕ1 + l cosϕ2 + l’’3 cosϕ3 + l4 cosϕ4 – xC1)2 + (r sinϕ1 + l sinϕ2 + l’’3 sinϕ3 + l4 sinϕ4 – yC1)2+ (r cosϕ1 + l cosϕ2 + l’3cos(ϕ3 - α)+ l5 cosϕ5 – xD1)2 + (r sinϕ1 + l sinϕ2 + l’3sin(ϕ3 - α)+ l5sinϕ5 – yD1)2 Æ Min Kết quả được cho trong bảng I.2 với zmin = 0. Bảng I.2. Kết quả tính toán giá trị các thông số của sàng phân loại ϕ1 ∈[0,2π] ϕ2 ∈[0,π] ϕ3∈[0,π] ϕ4∈[0,π] ϕ5∈[0,π] 0 0,226128 0,551311 1,783873 1,666775 π/18 0,199269 0,550518 1,784628 1,670250 2π/18 0,170835 0,550590 1,782751 1,668853 3π/18 0,143343 0,550490 1,778826 1,663697 4π/18 0,112669 0,552073 1,770032 1,652171 5π/18 0,090986 0,551991 1,759350 1,639575 6π/18 0,066036 0,553576 1,745374 1,622823 7π/18 0,051284 0,554296 1,730174 1,602970 8π/18 0,039053 0,555262 1,713242 1,581813 9π/18 0,033773 0,556277 1,695605 1,560720 8 2 MÔ HÌNH QUY HOẠCH ĐA MỤC TIÊU TUYẾN TÍNH VÀ PHI TUYẾN 2.1. Giới thiệu bài toán quy hoạch đa mục tiêu Trong các bài toán kỹ thuật, công nghệ, quản lý kinh tế, nông nghiệp v.v... nảy sinh từ thực tế, chúng ta thường phải xem xét đồng thời một lúc nhiều mục tiêu. Các mục tiêu này thường là khác về thứ nguyên, tức là chúng được đo bởi các đơn vị khác nhau. Những tình huống như vậy tạo ra các bài toán đa mục tiêu. Người kỹ sư / người ra quyết định lúc này cần phải tối ưu hoá (cực đại hoá hoặc cực tiểu hoá tuỳ theo tình huống thực tế) không phải là chỉ một mục tiêu nào đó, mà là đồng thời tất cả các mục tiêu đã đặt ra. Tuy nhiên, các mục tiêu này thường đối chọi cạnh tranh với nhau. Việc làm tốt hơn mục tiêu này thường dẫn tới việc làm xấu đi một số mục tiêu khác. Vì vậy việc giải các bài toán tối ưu đa mục tiêu, tức là tìm ra một phương án khả thi tốt nhất theo một nghĩa nào đó, thực chất chính là một bài toán ra quyết định. Có thể thấy lại ở đây một lần nữa khẳng định " Tối ưu hoá chính là một công cụ định lượng chủ yếu nhất của quá trình ra quyết định". Hiện tại các tài liệu, sách chuyên khảo, tạp chí cập nhật về các lĩnh vực liên ngành giữa Toán, Vận trù học, Khoa học Quản lý, Tin học, Công nghệ, Kinh tế, Nông nghiệp... đề cập rất nhiều tới bài toán tối ưu đa mục tiêu. Vấn đề nghiên cứu cơ sở lý thuyết, thuật toán, lập mô hình, xây dựng hệ máy tính trợ giúp quyết định, và áp dụng các mô hình tối ưu đa mục tiêu cho các quá trình công nghệ, quản lý... là một vấn đề liên ngành được rất nhiều nhà nghiên cứu khoa học và kĩ sư thực hành quan tâm. Mô hình tối ưu đa mục tiêu có dạng sau đây: Min fj(X), X = (x1, x2, , xn) j = 1, 2, , p (p ≥2) với: (i) gj(X) ≤ 0, j = 1, 2, , k, (ii) gj(X) = 0, j = k+1, k+2, , m, Trong các bài toán thực tế có thể bổ sung các ràng buộc đạng: (iii) ai ≤ xi ≤ bi, i = 1, 2, , n. Trong mô hình này, ta có p mục tiêu cần tối ưu hoá, các hệ số của các hàm mục tiêu và ràng buộc nói chung được giả sử là các giá trị thực xác định (cũng gọi là giá trị rõ). Trong trường hợp có ít nhất một trong các hàm mục tiêu hay các hàm ràng buộc là hàm phi tuyến, chúng ta có bài toán quy hoạch đa mục tiêu phi tuyến. Còn nếu tất cả các hàm mục tiêu và các hàm ràng buộc đều là hàm tuyến tính, chúng ta có mô hình quy hoạch tuyến tính đa mục tiêu với dạng chính tắc như sau: Min CX với ràng buộc DX ∈ , trong đó: C là ma trận cấp p×m D = { ∈X Rn : AX ≤ B, X ≥ 0 } với A là ma trận cấp m x n và ∈B Rm 9 2.2. Các khái niệm cơ bản của bài toán tối ưu đa mục tiêu Khái niệm then chốt trong tối ưu hoá đa mục tiêu là khái niệm phương án tối ưu Pareto. Định nghĩa 1: Một phương án tối ưu Pareto X* có tính chất sau đây: i) Trước hết nó phải thuộc vào miền các phương án khả thi của bài toán, tức là phải thoả mãn tất cả các ràng buộc: X* ∈ D. ii) Với mọi phương án khả thi khác X ∈ D mà có một mục tiêu nào đó tốt hơn (tồn tại chỉ số i sao cho fi(X) tốt hơn fi(X*)) thì cũng phải có ít nhất một mục tiêu khác xấu hơn (tồn tại j ≠ i sao cho fj(X) xấu hơn fj (X*)). Nói một cách khác, không tồn tại một phương án khả thi nào X ∈ D có thể trội hơn X* trên tổng thể. Định nghĩa 2: Giải bài toán tối ưu toàn cục đa mục tiêu là chọn ra từ tập hợp P các phương án tối ưu Pareto của bài toán một (một số) phương án tốt nhất theo một nghĩa nào đó dựa trên cơ cấu ưu tiên của người ra quyết định. Các phương án như vậy còn được gọi là phương án thoả dụng. Cách 1: Bằng một phương pháp tối ưu toán học thích hợp tìm ra tập hợp P tất cả các phương án tối ưu Pareto. Người ra quyết định sẽ đề ra cơ cấu ưu tiên của mình đối với tập P. Lúc này các phương pháp toán chẳng hạn như giải tích phân loại, các phương pháp lọc v.v được áp dụng để tìm ra phương án tối ưu cho bài toán đa mục tiêu ban đầu. Cách 2: Việc tìm tập hợp P trong trường hợp các bài toán tối ưu phi tuyến là khá khó, nếu không nói là không thể tìm được. Vì vậy, so với cách 1, cách 2 sẽ tiến hành theo trình tự ngược lại. Trước hết người ra quyết định sẽ đề ra cơ cấu ưu tiên của mình. Dựa vào cơ cấu ưu tiên đó, các mục tiêu sẽ được tổ hợp vào một mục tiêu duy nhất, tiêu biểu cho hàm tổng tiện ích của bài toán. Bài toán tối ưu với hàm mục tiêu tổ hợp này sẽ được giải bằng một phương pháp tối ưu toán học thích hợp, để tìm ra một (hoặc một số) phương án tối ưu Pareto. Lúc này, người ra quyết định sẽ chọn ra trong số các phương án tối ưu Pareto đó một phương án tốt nhất. Chúng ta sẽ tiếp tục phân tích cách thứ 2. Rõ ràng, người ra quyết định không thể đề ra cơ cấu ưu tiên của mình một cách chính xác ngay từ đầu. Trong quá trình giải bài toán, trong mỗi bước lặp, sau khi xem xét lại cơ cấu ưu tiên đã đề ra, cũng như phương án tối ưu trung gian, người ra quyết định có thể dựa vào các thông tin đó để thay đổi lại cơ cấu ưu tiên của mình. Sau đó, quá trình giải lại được tiếp tục, cho tới khi một phương án tối ưu cuối cùng được đưa ra. Định nghĩa 3: Phương pháp giải bài toán tối ưu dựa trên sự trợ giúp của hệ máy tính, nhằm giúp người ra quyết định từng bước thay đổi các quyết định trung gian một cách thích hợp để đi tới một phương án tối ưu Pareto thoả mãn nhất, được gọi là phương pháp tương tác người - máy tính. Cho tới thời điểm hiện nay, hàng chục phương pháp giải tương tác bài toán tối ưu đa mục tiêu đã được đề cập tới trong các tạp chí chuyên ngành, và đa số chúng đều 10 có những ứng dụng rất thành công trong nhiều lĩnh vực. Một trong các lớp phương pháp quan trọng và khá thuận tiện cho người sử dụng là phương pháp tương tác người - máy tính giải bài toán tối ưu đa mục tiêu với các yếu tố cấu thành sau: - Cơ cấu ưu tiên của người ra quyết định và hàm tổ hợp tương ứng. - Kiểu tương tác người - máy tính: các thông tin nào máy tính phải đưa ra lại trong các bước lặp trung gian, và cách thay đổi các thông số của cơ cấu ưu tiên từ phía người ra quyết định. - Kỹ thuật tối ưu toán học được xây dựng dựa trên lý thuyết tối ưu hoá nhằm tìm ra các phương án tối ưu Pareto cho các bài toán cần giải trong các bước lặp trung gian. 2.3. Các ví dụ minh hoạ bài toán quy hoạch đa mục tiêu Bài toán xác định cơ cấu cây trồng (Mô hình quy hoạch tuyến tính đa mục tiêu giải quyết vấn đề đánh giá hiệu quả sử dụng đất và xác định cơ cấu cây trồng xã Nhân Chính, huyện Lý Nhân, tỉnh Hà Nam) Trong nhiều mô hình quy hoạch đất, cũng như xây dựng cơ cấu cây trồng hợp lý các mô hình tối ưu được sử dụng rộng rãi. Ở nước ngoài, cũng như ở Việt Nam các phương pháp mô hình toán như vậy đã đem lại hiệu quả cao trong việc quản lý sử dụng đất cũng như việc bố trí một cơ cấu cây trồng thích hợp. Tuy nhiên đa số các mô hình tối ưu đó là mô hình một mục tiêu. Ngày nay, mô hình toán tối ưu nhiều mục tiêu được áp dụng trong nhiều lĩnh vực, đặc biệt trong lĩnh vực quản lý kinh tế, thiết kế chế tạo. Ví dụ này trình bày một cách áp dụng mô hình tối ưu nhiều mục tiêu cho bài toán xác định cơ cấu cây trồng - sử dụng đất trên địa bàn một xã vùng đồng bằng sông Hồng. Bài toán tối ưu nhiều mục tiêu nhằm đạt tới việc tối ưu hoá (cực đại hoá ) các mục tiêu sau : - Hiệu quả kinh tế của việc bố trí cơ cấu cây trồng và sử dụng đất. - Độ thích hợp tổng thể của cơ cấu cây trồng. - Hiệu quả môi trường. Bài toán có những điều kiện ràng buộc như sau về mặt sử dụng đất và cơ cấu cây trồng: Qua nghiên cứu về đất, toàn bộ diện tích canh tác được chia ra làm M loại. Còn kết quả nghiên cứu về trồng trọt cho phép áp dụng trên các loại đất đó N công thức luân canh, được đánh giá ở hai mức thích hợp. Độ thích hợp là 1 nếu công thức luân canh là thích hợp và là 2 nếu công thức luân canh ít thích hợp. Ký hiệu diện tích đất loại k là bk với k = 1, 2, ... M ; còn xijk (các biến quyết định) là diện tích luân canh công thức i (i = 1, 2, ..., N) với độ thích hợp j (j = 1, 2) trên đất loại k. Đặt aịjk là hệ số của xijk trong ràng buộc về diện tích đất loại k, ta sẽ có aijk = 0 nếu ta không áp dụng công thức luân canh i cho đất loại k; và aijk = 1 nếu công thức luân canh i được áp dụng cho đất loại k. Vậy các điều kiện ràng buộc của bài toán là : 11 ⎪⎩ ⎪⎨ ⎧ ∀∀∀≥ ==∑∑ == kjiX Mkbxa kji j kkjikji N i ,,0 )...,,2,1( 2 11 Ký hiệu cịjk là hệ số lợi nhuận trên một đơn vị diện tích của công thức luân canh i trên đất loại k với độ thích hợp j, thì mục tiêu về hiệu quả kinh tế được viết như sau : ∑∑∑ === → 2 111 j kjikjikji N i M k Maxxac (Mục tiêu 1 : Cực đại hoá lợi nhuận tổng) Để cực đại hoá độ thích hợp tổng thể, cần cực đại hoá tổng diện tích công thức luân canh có độ thích hợp loại 1. Vậy ta có mục tiêu sau : Maxxa kiki N i M k →∑∑ == 11 11 (Mục tiêu 2 : Cực đại hoá độ thích hợp tổng thể) Một mục tiêu nữa cần xem xét, như trên đã nói, là hiệu quả môi trường của cơ cấu cây trồng - sử dụng đất. Đây là một vấn đề thực tiễn, tuy nhiên các đánh giá về môi trường rất khó định lượng. Để xác định hệ số môi trường cho các công thức luân canh, chúng tôi dùng phương pháp tổng hợp ý kiến chuyên gia. Mỗi chuyên gia sẽ đưa ra đánh giá của mình về hiệu quả môi trường cho các công thức luân canh ở các mức : Tốt, khá, trung bình, xấu. Sau đó các mức đó sẽ được định lượng bởi các số 100, 75, 50 và 25. Tỉ lệ % các ý kiến cho từng mức đánh giá sẽ được coi là xác suất thực nghiệm của các giá trị trên. Và do đó mỗi công thức luân canh thứ i sẽ ứng với một cặp số mi (kỳ vọng) và σi (độ lệch tiêu chuẩn) của phân phối thực nghiệm thu được. Thay cho các phân phối xác suất thực nghiệm, ta sẽ xem xét hệ số mờ im~ = ( mi - 3σi, mi , mi + 3σ) của hiệu quả môi trường cho công thức luân canh thứ i. Do đó mục tiêu hiệu quả môi trường là mục tiêu mờ được viết như sau: ∑∑∑ === → 2 111 ~ j kjikjii N i M k Maxxam (Mục tiêu 3 : Cực đại hoá hiệu quả môi trường) Tóm lại ta có mô hình tối ưu (mờ) ba mục tiêu sau đây : ∑∑∑ === → 2 111 j kjikjikji N i M k Maxxac Maxxa kiki N i M k →∑∑ == 11 11 ∑∑∑ === → 2 111 ~ j kjikjii N i M k Maxxam 12 Với các ràng buộc : ⎪⎩ ⎪⎨ ⎧ ∀∀∀≥ ==∑∑ == kjiX Mkbxa kji j kkjikji N i ,,0 )...,,2,1( 2 11 Việc xử lý mục tiêu mờ (mục tiêu cực đại hoá hiệu quả môi trường) được tiến hành bằng cách cực đaị hoá hiệu quả môi trường với hệ số rõ mi và có thể kèm thêm một ràng buộc bổ sung : ∑∑∑ === ≥− 2 111 )( j kjikjiii N i M k exam ε Ở đây εi được tuỳ ý chọn thoả mãn điều kiện : 0 < εi < 3σi . Thông thường ta chọn εi = 90% × 3σi hay εi = 2,7 σi. Còn số e là một số được lựa chọn thích hợp từ những thông tin có thể khai thác từ bài toán trên (chẳng hạn từ thông tin của bảng pay- off), thường gọi là ngưỡng tối thiểu có thể chấp nhận của mục tiêu. Mô hình tối ưu (mờ) ba mục tiêu chính là một mô hình ra quyết định nhiều mục tiêu, và được giải bằng phương pháp đối thoại người ra quyết định - máy tính, nhằm giúp người ra quyết định từng bước tìm hiểu và thích nghi với các thông tin nội tại của mô hình, để cuối cùng đi tới một lời thoả mãn nhất. Các mục tiêu của mô hình được chuyển sang các mục tiêu mờ phản ánh độ thoả dụng của người ra quyết định. Đây là một cách làm rất hợp lý, bởi vì các đơn vị đo khác nhau của các mục tiêu được chuyển vào đơn vị thống nhất đo độ thoả dụng của người ra quyết định. Đây cũng là một lợi thế của tối ưu mờ (Fuzzy Optimization) bên cạnh rất nhiều ưu điểm khác. Số liệu thực tế qua khảo sát cơ cấu cây trồng, sử dụng đất ở địa phương Qua khảo sát thực tế các số liệu về các mặt sau được thu thập: - Hiệu quả kinh tế của một số công thức luân canh (bảng I.3). - Bố trí cơ cấu cây trồng theo các đơn vị đất (bảng I.4). - Số liệu về 15 loại đơn vị đất với các đặc tính đã được khảo sát (bảng I.5). - Tổng hợp các phiếu đánh giá hiệu quả môi trường một số công thức luân canh áp dụng cho vùng đồng bằng sông Hồng (bảng I.6). Bảng I.3. Hiệu quả kinh tế một số công thức trồng trọt luân canh Công thức trồng trọt Giá trị SX (ngànđ/ha) Năng suất (tạ/ha) Chi phí trung gian (ngàn đ) Yêu cầu lao động (công/ha) Giá trị tăng thêm (ngàn/ha) I.Bí xanh Cà chua Su hào 36634 39000 16064 457,92 65,00 200,80 9431 5879 4109 730 554 517 27203 33121 11955 Tổng 91698 19419 1801 72279 II.Bí xanh Cà chua Bắp cải 36634 39000 15050 457,92 65,00 215,00 9431 5879 4173 730 554 462 27203 33121 10877 Tổng 90684 19483 1746 71201 III.Bí xanh Bí xanh 36634 31948 457,92 399,35 9431 10056 730 758 27203 21892 13 Bắp cải 15050 215,00 4173 462 10877 Tổng 83632 23660 1950 59972 IV. L.xuân TQ L.mùa sớm Su hào Su hào 11223 8796 16060 14898 66,00 43,98 146,00 186,23 3557 3601 4048 4837 310 277 376 356 7666 5195 12012 10061 Tổng 50987 16043 1319 34934 V. L.xuân TQ L.mùa TQ Hành Tây 11223 11700 25178 66,00 65,00 209,82 3557 3468 7201 310 310 695 7666 8232 17977 Tổng 48101 14226 1315 33875 VI. L.xuân TQ Lúa mùa Dưa chuột 11223 11700 10139 66,00 65,00 207,79 3557 3468 7737 310 310 450 7666 8232 2402 Tổng 33062 14762 1070 18300 VII. L.xuân TQ L.mùa TQ Khoai tây 11223 11700 25401 66,00 65,00 169,34 3557 3468 6213 310 310 425 7666 8232 19188 Tổng 48324 13238 1045 35086 VIII. Cà chua Lúa mùa Bắp cải 10500 11700 15050 105,00 65,00 215,00 5879 3468 4173 554 310 462 4621 8232 10877 Tổng 37250 13520 1326 23730 IX. Lúa xuân Lúa mùa 11223 11700 66,00 65,00 3557 3468 310 310 7666 8232 Tổng 22923 7025 620 15898 X. Bí xanh Lúa mùa Bắp cải 36634 11700 15050 66,00 65.00 215,00 9431 3468 4173 310 310 462 27203 8232 10877 Tổng 63384 17072 1082 46312 XI. Bí xanh Lúa mùa Khoai tây 36634 11700 25401 457,92 65,00 169,34 9431 3468 6213 310 310 425 27203 8232 19188 Tổng 73735 19112 1465 54623 Bảng I.4. Bố trí cơ cấu cây trồng theo các đơn vị đất 1. (92,87 ha) VII. Lúa xuân - Lúa mùa - Khoai tây (2) VIII. Cà chua - Lúa mùa - Bắp cải (1) X. Bí xanh - Lúa mùa - Bắp cải (1) 2 (27,62) X. Bí xanh - Lúa mùa - Bắp cải (2) I. Bí xanh - Cà chua - Su hào (1) III. Bí xanh - Bí xanh - Bắp cải (1) 3, 4, 13, 14 (41,52) II. Bí xanh - Cà chua - Bắp cải (1) I. Bí xanh - Cà chua - Su hào (1) XI. Bí xanh - Lúa mùa - Khoai tây (2) 14 5, 6, 7 (163,15 ha) IX. Lúa xuân - Lúa mùa - ải (2) 8, 15 (16,49) VII. Lúa xuân - Lúa mùa - Khoai tây (2) V. Lúa xuân - Lúa mùa - Hành tây (2) X. Bí xanh - Lúa mùa - Bắp cải (1) 9, 10, 11 (85,02ha) VII. Lúa xuân - Lúa mùa - Khoai tây (1) V. Lúa xuân - Lùa mùa - Hành tây (1) IV. Lúa xuân-Lúa mùa-Su hào-Su hào (2) 12 (48,49) III. Bí xanh - Bí xanh - Bắp cải (1) X. Bí xanh - Lúa mùa - Bắp cải (2) Ghi chú : (1) Thích hợp (2) Kém thích hợp Bảng I.5. Diện tích đơn vị đất đai 1 92,87 ha 9 43,88 ha 2 27,62 10 8,85 3 5,42 11 32,29 4 12,02 12 48,49 5 99,89 13 14,03 6 27,96 14 10,05 7 35,30 15 8,41 8 8,08 Bảng I.6. Tổng hợp phiếu đánh giá hiệu quả môi trường của một số công thức luân canh áp dụng cho vùng đồng bằng sông Hồng Công thức Công thức Mức độ luân canh Tốt (%) Khá (%) T.B (%) Xấu (%) 1 Bí xanh-Cà chua-Su hào 45,4 27,3 27,3 2 Bí xanh-Cà chua-Bắp cải 45,4 18,2 36,4 3 Bí xanh-Bí xanh-Bắp cải 36,4 36,4 27,2 4 LX-LM sớm-Su hào-Su hào 72,7 18,2 9,1 5 Lúa xuân-LM sớm-Hành tây 72,7 18,2 9,1 6 Lúa xuân-LM sớm-Dưa chuột 54,5 27,3 18,2 7 Lúa xuân-Lúa mùa-Khoai tây 81,8 18,2 8 Cà chua-Lúa mùa sớm-Bắp cải 27,3 45,4 27,3 9 Lúa xuân-Lúa mùa 45,4 36,4 18,2 10 Bí xanh-Lúa mùa sớm-Bắp cải 27,3 54,5 18,2 11 Bí xanh-LM sớm-Khoai tây 27,3 45,4 27,3 Thiết lập mô hình đa mục tiêu xác định cơ cấu cây trồng Để có thể chọn những công thức trồng trọt phù hợp với điều kiện đất đai, đồng thời đảm bảo đạt hiệu quả mong muốn, cần xét ba mục tiêu sau: i) Hiệu quả kinh tế, ii) Độ thích hợp đất đai, iii) Hiệu quả môi trường. Tiến hành thiết lập mô hình, trước hết chọn các biến quyết định. Dựa vào kết quả các dữ liệu đã thu được, đặt tên các biến như sau: x1 (x112): diện tích trồng các loại cây công thức 1, ở độ thích hợp đất (1) trên khu đất (2) x2 (x113): diện tích trồng các loại cây công thức 1, ở độ thích hợp đất (1) trên khu đất (3) x3 (x213): diện tích trồng các loại cây công thức 2, ở độ thích hợp đất (1) trên khu đất (3) 15 x4 (x312): diện tích trồng các loại cây công thức 3, ở độ thích hợp đất (1) trên khu đất (2) x5 (x317): diện tích trồng các loại cây công thức 3, ở độ thích hợp đất (1) trên khu đất (7) x6 (x516): diện tích trồng các loại cây công thức 5, ở độ thích hợp đất (1) trên khu đất (6) x7 (x716): diện tích trồng các loại cây công thức 7, ở độ thích hợp đất (1) trên khu đất (6) x8 (x811): diện tích trồng các loại cây công thức 8, ở độ thích hợp đất (1) trên khu đất (1) x9 (x914): diện tích trồng các loại cây công thức 9, ở độ thích hợp đất (1) trên khu đất (4) x10 (x10,11): diện tích trồng các loại cây công thức 10, ở độ thích hợp đất (1) trên khu đất (1) x11 (x10,15): diện tích trồng các loại cây công thức 10, ở độ thích hợp đất (1) trên khu đất (5) x12 (x426): diện tích trồng các loại cây công thức 4, ở độ thích hợp đất (2) trên khu đất (6) x13 (x525): diện tích trồng các loại cây công thức 5, ở độ thích hợp đất (2) trên khu đất (5) x14 (x721): diện tích trồng các loại cây công thức 7, ở độ thích hợp đất (2) trên khu đất (1) x15 (x725): diện tích trồng các loại cây công thức 7, ở độ thích hợp đất (2) trên khu đất (5) x16 (x10,22): diện tích trồng các loại cây công thức 10, ở độ thích hợp đất (2) trên khu đất (2) x17 (x10,27): diện tích trồng các loại cây công thức 10, ở độ thích hợp đất (2) trên khu đất (7) x18 (x11,23): diện tích trồng các loại cây công thức 1, ở độ thích hợp đất (2) trên khu đất (3). Với các số liệu thực tế thu được qua khảo sát, chúng tôi có mô hình tối ưu (mờ) sau đây : Mục tiêu 1: z1 = 72279(x112 + x113) + 71201x213 + 59972 (x312 + x317) + 34934 x426 + 33875 (x516 + x525 ) + 35086 (x716 + x721 + x725) + 23730 x811 + 15898 x914 + 46312 (x10, 11 + x10, 15 + x10, 22 + x10, 27) + 54623 x11, 23 → Max Mục tiêu 2 : z2 = x112 + x113 + x213 + x312 + x317 + x516 + x716 + x811 + x914 + x10, 11 + x10, 15 → Max Mục tiêu 3 : z3 = 1~m (x112 + x113) + 2~m x213 + 3~m (x312 + x317) + 4~m x426 + 5~m (x516 + x525 ) + 7~m (x716 + x721 + x725) + 8~m x811 + 9~m x914 + 10~m (x10, 11 + x10, 15 + x10, 22 + x10, 27) + 11~m x11, 23 → Max Với các ràng buộc sau đây : x721 + x811 + x10, 11 = 92,87 (ha); x112 + x312 + x10, 22 = 27,62; x113 + x213 + x11, 23 = 41,52; x914 = 163,15; x525 + x725 + x10, 25 = 16,49; x426 + x516 + x716 = 85,02; x317 + x10, 27 = 48,49; xịk ≥ 0 ∀i, ∀j , ∀k. Ở mô hình trên 1~m được xác định như sau : Trước hết ta tính m1 = 27,3% × 25 + 27,3% × 50 + 45,4% × 75 σ1 = 21222 m75%4,4550%3,2725%3,27 −×+×+× Sau đó có : 1~m = (m1 - 3σ1 , m1 , m1 + 3σ1 ). Tương tự có thể tính được các hệ số mờ khác. Do đó mục tiêu 3 được được thay bởi: Mục tiêu 3 : z3 = 54,525(x112 + x113) + 52,25x213 + 52,3(x312 + x317) + 90,9 x426 + 84,075(x516 + x525 ) + 95,375(x716 + x721 + x725) + 75x811 + 81,8x914 + 77,275(x10, 11 + x10, 15 + x10, 22 + x10, 27) + 77,275x11, 23 → Max 16 Các bước giải bài toán tối ưu cơ cấu cây trồng Tiến hành giải bài toán trên bằng phần mềm MULTIOPT (xem chương IV) theo các bước sau: Trước hết, nhằm giúp người ra quyết định xác định hàm thoả dụng, bài toán tối ưu cho từng mục tiêu riêng rẽ (như vậy có ba bài toán tối ưu một mục tiêu) sẽ được giải quyết. Riêng mục tiêu về hiệu quả môi trường sẽ lấy mi thay thế cho im~ đối với mọi công thức canh tác i. Màn hình máy tính sẽ thông báo về ba lời giải tối ưu cho ba bài toán, ký hiệu là X1, X2, và X3. X1 : x112 = 27,62; x113 = 41,52 ; x312 = 48,49 ; x716 = 85,02 ; x914 = 163,15; x10, 11 = 92,87 và x10, 15 = 16,49 (các biến khác có giá trị bằng 0). X2 : x112 = 27,62 ; x113 = 41,52 ; x317 = 48,49 ; x516 = 85,02 ; x811 = 92,87; x914 = 163,15; x10, 15 = 16,49 (các biến khác có giá trị bằng 0). X3 : x716 = 85,02 ; x721 = 92,87 ; x725 = 16,49 ; x914 = 163,15 ; x10, 22 = 27,62; x10, 23 = 48,49 và x11, 23 = 41,52 (các biến khác có giá trị bằng 0). Lúc này ta có thông tin pay- off cho ở bảng I.7. Bảng I.7. Thông tin pay- off Phương án Mục tiêu 1 Mục tiêu 2 Mục tiêu 3 X1 18544625,75 475,16 36218,401 X2 16344476,19 475,16 35039,98 X3 15206528,66 248,17 40985,02 Dựa trên các thông tin pay-off trên, máy tính sẽ giúp người ra quyết định xây dựng được các hàm thoả dụng cho các mục tiêu đã đặt ra. Các hàm đó còn gọi là các hàm liên thuộc mờ (fuzzy membership functions). Chẳng hạn hàm thoả dụng cho mục tiêu hiệu quả kinh tế là : μ1 (z1) = w 1 B 1 w 11 ZZ ZZ − − → Max Với z1B (tốt nhất) = 18544625,75 z1w (xấu nhất) = 15206528,66. Sau khi có các hàm thoả dụng, hàm liên hợp (aggregation function) được thiết lập cho các hàm thoả dụng đó để có mục tiêu sau : w1 μ1(z1) + w2 μ2 (z2) + w3 μ3 (z3) → Max Với w1 + w2 + w3 = 1, 0 ≤ w1 , w2 , w3 ≤ 1 Các hệ số w1 , w2 , w3 được gọi là các trọng số phản ánh tầm quan trọng của từng hàm thoả dụng thành phần trong hàm liên hợp. Bằng cách thay đổi giá trị w1 , w2 , w3, người ra quyết định có thể tìm ra được các phương án thoả dụng thích hợp (có tính chất tối ưu Pareto yếu). Xem xét các phương án đó cùng với độ thoả dụng đạt được cho từng mục tiêu (có thể dựa vào phương pháp ra quyết định tập thể - group decision making) người ra quyết định có thể đi tới một quyết định hợp lý về cơ cấu cây trồng - sử dụng đất như trong bảng I.8. 17 Bảng I.8. Kết quả lựa chọn phương án tối ưu Đơn vị đất Khu vực Diện tích (ha) Phương án 1 Phương án 2 Phương án 3 Phương án chọn 1 1 92.,8 7 Bí xanh (xuân) – LM – bắp cải Cà chua (đông) – LM – bắp cải Lúa xuân (TQ) - LM (TQ) – khoai tây Bí xanh (xuân) – LM – bắp cải 2 2 27,6 2 Bí xanh (xuân) – cà chua (hè thu)- su hào Bí xanh (xuân) – cà chua (hè thu)- su hào Bí xanh (xuân) – LM – bắp cải Bí xanh (xuân) – cà chua (hè thu)- su hào 3; 4; 14; 13 3 41,5 2 Bí xanh (xuân) – cà chua (hè thu)- su hào Bí xanh (xuân) – cà chua (hè thu)- su hào Bí xanh (xuân) – LM – khoai tây Bí xanh (xuân) – cà chua (hè thu)- su hào 5; 6; 7 4 163, 15 Lúa xuân – lúa mùa Lúa xuân – lúa mùa Lúa xuân – lúa mùa Lúa xuân – lúa mùa 8; 15 5 16,4 9 Bí xanh (xuân) – LM – bắp cải Bí xanh (xuân) – LM – bắp cải Lúa xuân (TQ) - LM (TQ)– khoai tây Bí xanh (xuân) – LM – bắp cải (6.51234 ha) LX TQ) – LM (TQ)– khoai tây (9.97766 ha) 9; 10; 11 6 85,0 2 Lúa xuân (TQ) - LM (TQ) – khoai tây Bí xanh (xuân) – LM – bắp cải Lúa xuân (TQ) - LM (TQ)– khoai tây Lúa xuân (TQ) - LM (TQ)– khoai tây 12 7 48,4 9 Bí xanh (xuân)- bí xanh (mùa)- bắp cải Bí xanh (xuân) - bí xanh (mùa)- bắp cải Bí xanh (xuân) – LM – bắp cải Bí xanh (xuân)- bí xanh mùa)- bắp cải Trong bảng I.8, chúng ta có: Phương án 1: Các công thức trồng trọt cho hiệu quả kinh tế cao nhất. Phương án 2: Các công thức trồng trọt cho hiệu quả cao nhất về mức độ thích nghi đất đai. Phương án 3: Các công thức trồng trọt cho hiệu quả môi trường cao nhất . Phương án cho ở cột phương án chọn là phương án thoả dụng ứng bộ giá trị các trọng số: w1 = 0.4, w2 = 0.2, w3 = 0.4 khi sử dụng phần mềm MULTIOPT. Bài toán đánh giá hiệu quả sử dụng đất và lao động (Mô hình quy hoạch tuyến tính đa mục tiêu giải quyết vấn đề đánh giá hiệu quả sử dụng đất và lao động trên địa bàn xã Trâu Quỳ, huyện Gia Lâm, tỉnh Hà Nội) 18 Để quy hoạch sử dụng đất, đồng thời đảm bảo đạt hiệu quả môi trường, cần xem xét năm mục tiêu sau: i) Tổng lợi nhuận, ii) Hiệu quả sử dụng vốn, iii) Giá trị ngày công lao động, iv) Số công lao động, vi) Hiệu quả môi trường. Để tiến hành giải bài toán, trước hết phải chọn các biến quyết định. Dựa vào cơ cấu cây trồng xã năm 1999, các biến sau được lựa chọn: x1 : Diện tích trồng lúa xuân (ha), x2 : Diện tích trồng lúa mùa (ha), x3 : Diện tích trồng ngô (ha), x4 : Diện tích trồng đậu tương (ha), x5 : Diện tích trồng khoai tây (ha), x6 : Diện tích trồng rau (ha), x7 : Diện tích trồng mùi (ha), x8 : Diện tích trồng táo (ha), x9 : Diện tích trồng nhãn (ha), x10 : Diện tích trồng xoài (ha). Các mục tiêu cần cực đại hoá là: z1 = 4,48 x1 + 4,2 x2 + 2,59 x3 + 0,98 x4 + 5,8 x5 + 15,61 x6 + 29,67 x7 + 39,21 x8 + 116,58 x9 + 105,13 x10 z2 = 0,6205 x1 + 0,5915 x2 +0,465 x3 + 0,1583 x4 + 0,7065 x5 + 0,5864 x6 + 1,2996 x7 + 1,2735 x8 + 1,1726 x9 + 1,756 x10 z3 = 0,0217 x1 + 0,0206 x2 + 0,0154 x3 + 0,0045 x4 + 0,0248 x5 + 0,0109 x6 + 0,0241 x7 + 0,0349 x8 + 0,09 x9 + 0,0811 x10 z4 = 206 x1 + 204 x2 + 168 x3 + 216 x4 + 234 x5 + 1428 x6 + 1232 x7 + 1124 x8 + 1296 x9 + 1296 x10 z5 = 0,7 x1 + 0,778 x2 + 1,273 x3 + 1,75 x4 + x5 + 0,368 x6 + 0,875 x7 + 3 x8 + 3 x9 + 3 x10 Với các ràng buộc sau (về cơ cấu diện tích đất canh tác): x1≤ 189,6407; x2 ≤ 189,6407; x3 ≤ 17,4931; x4 ≤ 17,4931; x5 ≤17,4931; x6 ≤189,6407; x7 ≤ 17,4931; x8 ≤18; x9 ≤18; x10 ≤ 18. Diện tích trồng rau trên cả đất ba vụ và đất chuyên màu: x6 ≥ 26,4 Diện tích đất trồng cây vụ đông trên đất ba vụ: x3 + x4 + x5 + x6 + x7 = 43,8931 Diện tích đất trồng các cây ăn quả trên đất trồng cây hàng năm khác: x8 + x9 + x10 = 18 Điều kiện để có lợi nhuận là: 4,48 x1 + 4,2 x2 + 2,59 x3 + 0,98 x4 + 5,8 x5 + 15,61 x6 + 29,67 x7 + 39,21 x8 + 116,58 x9 + 105,13 x10 > 0 Điều kiện về sản lượng lương thực: Các cây lương thực của xã gồm lúa xuân, lúa mùa và ngô: 5,14 x1 + 4,98 x2 + 3,77 x3 ≥ 1700,5 Điều kiện không âm của bài toán: xi (i = 1, 2, , 10) ≥ 0. Bảng I.9 cho phép so sánh kết quả sử dụng đất canh tác xã Trâu quỳ năm 1999 với các phương án tối ưu thu được khi áp dụng phần mềm MULTIOPT (xem chương IV ). 19 Bảng I.9. So sánh kết quả sử dụng đất canh tác STT Các chỉ tiêu Kết quả thực tế Kết quả mô hình 1 Diện tích các cây trồng Lúa xuân Lúa mùa Ngô Đậu tương Khoai tây Rau Mùi Táo Nhãn Xoài 189,64 189,64 2,53 0,82 0,91 27,73 12,4 9,27 4,85 3,88 189,64( x1) 189,64(x2) 0,00 0,00 0,00 26,40(x6) 17,49(x7) 0,00 0,00 18,00(x10) 2 3 4 5 6 7 8 Tổng thu nhập Tổng chi phí Tổng lợi nhuận Hiệu quả sử dụng vốn Tổng số công lao động Giá trị ngày công lao động Môi trường 8538,62 4750,00 3788,47 288,19 154462 9,74 360,72 9364,92 4895,38 4469,54(z1) 299,67(z2) 160331(z3) 10,19(z4) 359,31(z5) ( Giá trị các trọng số: w1= 0,1; w2 = 0,1; w3 = 0,2; w4 = 0,2; w5 = 0,4 ). Tối ưu hoá kết quả và hiệu quả kinh tế chăn nuôi cá(Mô hình quy hoạch phi tuyến đa mục tiêu giải quyết vấn đề tối ưu hoá kết quả và hiệu quả kinh tế chăn nuôi cá của các hộ nông dân nuôi cá tại huyện Văn Giang, tỉnh Hưng Yên) Việc các hộ nuôi cá quyết định nên sản xuất như thế nào, mức đầu tư bao nhiêu phụ thuộc chủ yếu vào những lợi thế và tiềm lực kinh tế của từng nông hộ. Những nông hộ có nhiều lao động dư thừa chú trọng đến thu nhập hỗn hợp hơn thu nhập ròng, còn những hộ có sẵn nguồn thức ăn tận dụng trong gia đình mà không mất tiền mua thì lại xem trọng giá trị sản xuất hơn. Chính vì vậy, cần không chỉ quan tâm đến giá trị sản xuất mà còn cần quan tâm cả đến thu nhập hỗn hợp và thu nhập ròng. Do đó, phải giải quyết bài toán tối ưu ba mục tiêu sau: i) Tối đa hoá giá trị sản xuất cá trên một ha (z1 → Max), 20 ii) Tối đa hoá thu nhập ròng trên một ha (z2 →Max), iii) Tối đa hoá thu nhập hỗn hợp trên một ha(z3 → Max). z1 = 19,387 x10,236 x20,103 x30,096 x40,056 x5 0,056 e0,168x6 e0,066x7 → Max z2 = z1 - TC → Max z3 = z1 - TC + x 3 → Max Trong bài toán đa mục tiêu trên, các ràng buộc là: - Với điều kiện thực tế tại địa phương và kết quả tính toán cho thấy mức đầu tư trên 75 triệu/ha sẽ không đem lại hiệu quả kinh tế cao. Do vậy tổng chi phí giới hạn ở mức 75 triệu/ha hay: 0 < TC ≤ 75. - Việc phân tích các kết quả tính toán cho thấy, với mức đầu tư tối ưu thì chi phí giống x1 không thể vượt quá 40% tổng chi phí, tương tự chi phí thức ăn x2, chi phí lao động x3, chi phí khấu hao x4 và chi khác x5 không thể vượt quá 30%, 20%, 15% và 15% so với tổng chi phí. Vì vậy, các ràng buộc cần xem xét là: - 0 < TC = x1 + x2 + x3 + x4 +x5 ≤ 75 - x1≤ 40; x2≤ 30; x3≤ 20; x4≤ 15; x5 ≤ 15 - x6 + x7 ≤ 1 (x6 , x7 ∈ {0, 1} ) Khi giải mô hình trên bằng phần mềm PRELIME ta thu được các kết quả cho trong bảng I.10. Bảng I.10. Kết quả tối ưu hoá từng mục tiêu đơn lẻ vùng đồng Tỉ lệ đầu tư(triệu/ha) Biến giả Giá trị Max (triêụ/ha) Phương án x1 x2 x3 x4 x5 Tổng chi x6 x7 GO NI MI X1 30,0 15,7 12,6 8,5 8,1 75,0 1 0 110,2 - - X2 20,0 8,8 8,2 4,8 4,8 46,5 1 0 - 38,4 - X3 23,7 10,3 20,0 5,6 5,6 65,2 1 0 - - 54,5 Trước hết, giá trị cực đại của từng mục tiêu đơn lẻ với mức đầu tư tương ứng đ- ược thể hiện trong bảng trên. Biến hình thức nuôi x6 = 1, nghĩa là ở vùng đồng KQ và HQKT đạt lớn nhất khi các hộ chuyển sang nuôi chuyên canh và ở mức đầu tư 75 triệu/ha, với cơ cấu đầu tư x1 : x2 : x3 : x4 : x5 = 30,0 : 15,7 : 12,6 : 8,5 : 8,1 sẽ cho ta giá trị sản xuất lớn nhất (đạt 110,15 triệu); còn ở mức đầu tư 46,5 triệu, với cơ cấu x1 : x2 : x3 : x4 : x5 = 20,0 : 8,8 : 8,2 : 4,8 : 4,8 thì thu được thu nhập ròng lớn nhất (đạt 38,42 triệu); các hộ cần cực đại hoá thu nhập hỗn hợp thì nên đầu tư ở mức 65,19 triệu, với cơ cấu đầu tư là x1 : x2 : x3 : x4 : x5 = 23,7 : 10,3 : 20,0 : 5,6 : 5,6. Như vậy, mức và cơ cấu đầu tư khác nhau làm cho KQ và HQKT thu được của các hộ cũng khác nhau. Trên tập {X1 , X2 , X3} ta thu được giá trị tốt nhất của GO, NI và MI lần lượt là: GOMax = 110,2, NIMax = 38,4, MIMax = 54,5 và giá trị thấp nhất của GO là 71,8, NI là 18,1, MI là 30,8. Từ đó ta xác định được các giải giá trị thích hợp cho GO, NI và MI tương ứng là: [71,8 - 110,2], [18,1 - 38,4] và [30,8 - 54,5]. Với các thông tin thu được, phần mềm PRELIME (xem bài báo của C. Mohan và Nguyen Hai Thanh,“An interactive satisficing 21 method for solving multiobjective mixed fuzzy-stochastic programming problems”, International Journal for Fuzzy Sets and Systems, Vol. 117, No.1, pp. 61-79, 2001) sẽ cho biết: Đối với từng mục tiêu z1 (GO), z2 (NI) và z3 (MI) các hàm thoả dụng lần lượt là μ1, μ2 và μ3 với các mức đạt được được thể hiện trên các đồ thị. - Đối với hàm thoả dụng về giá trị sản xuất: 0 nếu z1 ≤ 71,8 = a1 μ1 = 0,5 nếu z1 = 80,0 = b1 1 nếu z1 ≥110,2 = c1 - Đối với hàm thoả dụng về thu nhập ròng: 0 nếu z2 ≤ 15,1 = a2 μ2 = 0,5 nếu z2 = 34,0 = b2 1 nếu z2 ≥ 38,4 = c2 - Đối với hàm thoả dụng về thu nhập hỗn hợp: 0 nếu z3 ≤ 38,8 = a3 μ3 = 0,5 nếu z3 = 45,0 = b3 1 nếu z3 ≥ 54,5 = c3 Chúng ta có thể hiểu ý nghĩa của các hàm thoả dụng trên như sau: Đối với hàm z1 (GO) thì một phương án X = (x1, x2, x3, x4, x5, x6, x7) với giá trị z1 ≤ 110,2 triệu/ha = c1 cho ta độ thoả dụng μ1 = 1; với giá trị z1 ≥ 71.8 triệu/ha = a1 thì độ thoả dụng μ1 = 0; với z1 = 80,0 triệu/ha lại cho ta độ thoả dụng μ1 = 0,5. Các giá trị khác của z1 cho các độ thoả dụng được nội suy căn cứ vào các đồ thị trên. Giá trị z1 = 80,0 triệu/ha = b1 (tương ứng với độ thoả dụng 0,5) được gọi là giá trị chốt và được xác định cho bước lặp đầu tiên trong quá trình giải bài toán tối ưu đa mục tiêu. ở các bước sau giá trị chốt b1 có thể được thay đổi căn cứ vào thông tin nội tại của bài toán do máy tính đưa ra. Tương tự ta có thể giải thích ý nghĩa của các hàm thoả dụng μ2 và μ3 ứng với các mục tiêu z2 (NI) và z3 (MI). Một phương án X như vậy cho ta độ thoả dụng ứng với ba mục tiêu z1 (GO), z2 (NI), z3 (MI) là μ1, μ2 và μ3 . Từ ba hàm thoả dụng này ta thiết lập được hàm thoả dụng tổng hợp (Aggregation Utility Function) cho đồng thời ba mục tiêu như sau: μ = Max { Min [μ1 , μ2 , μ3 ] } Trong quá trình giải bài toán qua các bước lặp nếu thu được phương án X với μ > 0,5 thì điều đó có nghĩa là tất cả ba hàm thoả dụng μ1, μ2 và μ3 lớn hơn 0,5. Từ đó ta thấy các giá trị chốt b1, b2, b3 đã đặt ra là thấp nên ta có thể tăng (ít nhất) một trong ba giá trị chốt này, việc ưu tiên mục tiêu nào (z1, z2 hay z3) tuỳ thuộc vào sự quan tâm của nông hộ nhằm sử dụng có hiệu quả nhất nguồn lực của mình. Ngược lại, nếu μ < 0,5, thì ta giảm ít nhất một trong ba giá trị ưu tiên b1, b2 hay b3. Như vậy, trong quá trình giải các giá trị ưu tiên này sẽ được thương lượng với nhau để điều chỉnh sao cho thu được phương án X* thoả mãn nhất với độ thoả dụng tổng hợp μ gần với 0,5. 22 Kết quả tính toán bài toán đa mục tiêu ở vùng đồng khi sử dụng phần mềm PRELIME được thể hiện trong bảng I.11. Có thể thấy rằng, tuỳ vào mục đích của người ra quyết định mà có sự lựa chọn khác nhau về GO, NI hay MI. Đối với vùng đồng, các hộ đều lựa chọn hình thức nuôi chuyên canh (x6 = 1), còn mức và cơ cấu đầu tư tối ưu lại phụ thuộc vào mục đích của hộ sẽ chọn GO, NI hay MI yếu tố nào là quan trọng. Chẳng hạn, người thứ nhất quan tâm nhiều đến GO thì cần đầu tư ở mức 62,7 triệu/ha, tương ứng với cơ cấu đầu tư x1 : x2 : x3 : x4 : x5 = 25,0 : 11,2 : 14,3 : 6,3 : 6,0. Nhìn chung, giá trị lớn nhất đồng thời của ba mục tiêu không có sự thay đổi lớn qua các sự lựa chọn; GOMax đạt từ 95,4 đến 99,6 triệu/ha, NIMax đạt 35,6 đến 36,9 triệu/ha và MIMax đạt từ 52,1 đến 53,4 triệu/ha. Tổng mức đầu tư cũng không có sự biến động lớn, chúng chỉ dao động từ 58,8 đến 62,6 triệu/ha, ứng với các tỉ lệ cơ cấu cũng không có sự biến đổi nhiều. Bảng I.11. Kết quả bài toán đa mục tiêu vùng đồng Tỉ lệ đầu tư (triệu/ha) Biến giả Giá trị max (triệu/ha) Bước lặp thứ Độ thoả dụng (%) x1 x2 x3 x4 x5 Tổn g chi x6 x7 GO NI MI 1 0,82 25, 0 11, 2 14, 3 6,3 6,0 62,7 1 0 99,6 36,9 51, 2 2 0,74 22, 6 10, 0 15, 5 5,4 5,4 58,8 1 0 95,4 36,6 52, 1 3 0,49 23, 1 10, 3 17, 9 5,6 5,7 62,6 1 0 98,1 35,6 53, 4 4 0,46 22, 7 10, 0 16, 6 5,4 5,3 60,0 1 0 96,2 36,2 52, 7 5 0,50 22, 2 9,9 15, 0 5,3 5,3 57,7 1 0 94,5 36,8 51, 8 23 Chương II GIẢI BÀI TOÁN QUY HOẠCH TUYẾN TÍNH BẰNG PHƯƠNG PHÁP ĐƠN HÌNH 1. PHƯƠNG PHÁP ĐƠN HÌNH GIẢI BTQHTT DẠNG CHÍNH TẮC 1.1. Ví dụ . Xét BTQHTT: Max z = 8x1 + 6x2, với các ràng buộc 4x1 + 2x2 ≤ 60 2x1 + 4x2 ≤ 48 x1, x2 ≥ 0 Đưa BTQHTT dạng chuẩn tắc trên về dạng chính tắc bằng các biến bù không âm x3 và x4 như sau: Max z = 8x1 + 6x2 + 0x3 + 0x4 với các ràng buộc 4x1 + 2x2 + x3 = 60 2x1 + 4x2 + x4 = 48 x1, x2, x3, x4 ≥ 0. Chú ý. BTQHTT có dạng chính tắc là BTQHTT với các biến không âm, các ràng buộc có dấu “=”, hệ số vế phải của các ràng buộc không âm. Ngoài ra, mỗi phương trình bắt buộc phải có một biến đứng độc lập với hệ số +1. Cách lập và biến đổi các bảng đơn hình Để giải BTQHTT dạng chính tắc trên đây, cần lập một số bảng đơn hình như trong bảng II.1. Trước hết, cần điền số liệu của bài toán đã cho vào bảng đơn hình bước 1: – Cột 1 là cột hệ số hàm mục tiêu ứng với các biến cơ sở đã chọn. Phương án xuất phát có thể chọn là x1 = x2 = 0 (đây chính là điểm gốc toạ độ O(0, 0) trên hình II.1), do đó x3 = 60, x4 = 48. Như vậy tại bước này chúng ta chưa bước vào sản xuất, nên trong phương án chưa có đơn vị sản phẩm loại I hay loại II nào được sản xuất ra (chỉ “sản xuất” ra các lượng nguyên liệu dư thừa, ta cũng nói là các “sản phẩm” loại III và IV), và giá trị hàm mục tiêu z tạm thời bằng 0. Các biến bù có giá trị lớn hơn 0 có nghĩa là các nguyên liệu loại tương ứng chưa được sử dụng hết. Ta gọi các biến x3 và x4 là các biến cơ sở vì chúng có giá trị lớn hơn 0 còn x1 và x2 là các biến ngoài cơ sở vì chúng có giá trị bằng 0. Với bài toán có hai ràng buộc, tại mỗi bước chỉ có hai biến cơ sở. – Cột 2 là cột các biến cơ sở. Trong cột 3 (cột phương án) cần ghi các giá trị của các biến cơ sở đã chọn. – Các cột tiếp theo là các cột hệ số trong các điều kiện ràng buộc tương ứng với các biến x1, x2, x3 và x4 của bài toán đã cho. 24 Bảng II.1. Các bảng đơn hình giải BTQHTT c1 = 8 c2 = 6 c3 = 0 c4 = 0 Hệ số hàm mục tiêu cj Biến cơ sở Phương án x1 x2 x3 x4 Bảng đơn hình bước 1 0 0 x3 x4 60 48 4 2 2 4 1 0 0 1 Hàng z z0 = 0 z1 = 0 z2 = 0 z3 = 0 z4 = 0 Hàng Δj = cj – zj Δ1 = 8 Δ2 = 6 Δ3 = 0 Δ4 = 0 Bảng đơn hình bước 2 8 0 x1 x4 15 18 1 0 1/2 3 1/4 –1/2 0 1 Hàng z z0 = 120 z1 = 8 z2 = 4 z3 = 2 z4 = 0 Hàng Δj = cj – zj Δ1 = 0 Δ2 = 2 Δ3 = –2 Δ4 = 0 Bảng đơn hình bước 3 8 6 x1 x2 12 6 1 0 0 1 1/3 –1/6 –1/6 1/3 Hàng z z0 = 132 8 6 5/3 2/3 Hàng Δj = cj – zj 0 0 –5/3 –2/3 Phân tích bảng đơn hình bước 1 – Hệ số ứng với biến x1 trên hàng thứ nhất là a11 = 4 có nghĩa là tỷ lệ thay thế riêng giữa một đơn vị sản phẩm loại I và một đơn vị sản phẩm loại III là 4 (giải thích: xét phương trình (hay ràng buộc) thứ nhất 4x1 + 2x2 + x3 = 60, x1 tăng một đơn vị thì x3 phải giảm bốn đơn vị nếu giữ nguyên x2). Tương tự ta có thể giải thích được ý nghĩa của các hệ số aij khác cho trên hàng 1 và hàng 2 trong bảng đơn hình bước 1. – Chúng ta xét hàng z của bảng đơn hình. Để tính z1, cần áp dụng công thức z1 = (cột hệ số của hàm mục tiêu)× (cột hệ số của biến x1) = 0×4 + 0×2 = (giá một đơn vị sản phẩm loại III)×(tỷ lệ thay thế riêng loại I / loại III) + (giá một đơn vị sản phẩm loại IV)×(tỷ lệ thay thế riêng loại I / loại IV) = tổng chi phí phải bỏ ra khi đưa thêm một đơn vị sản phẩm loại I vào phương án sản xuất mới = 0. Các giá trị zj, với j = 1, 2, 3, 4, được tính tương tự và chính là các chi phí khi đưa một thêm một đơn vị sản phẩm loại xj vào phương án sản xuất mới. Còn z0 là giá trị của hàm mục tiêu đạt được tại phương án đang xét: z0 = (cột hệ số của hàm mục tiêu)× (cột phương án) = 0×60 + 0 × 48 = 0. – Trên hàng Δj cần ghi các giá trị Δj , j = 1, 2, 3, 4, tính theo công thức Δj = cj – zj = lợi nhuận / đơn vị sản phẩm – chi phí / đơn vị sản phẩm. Vậy Δj là "lãi biên" / một đơn vị sản phẩm khi đưa một thêm một đơn vị sản phẩm loại xj vào phương án sản xuất mới. Nếu Δj > 0 thì hàm mục tiêu còn tăng được khi ta đưa thêm các sản phẩm loại j vào phương án sản xuất mới. Có thể chứng minh được Δj chính là đạo hàm riêng jz / x∂ ∂ của hàm mục tiêu z theo biến xj . Như vậy, x1 tăng lên 1 thì z tăng lên 8 còn x2 tăng lên 1 thì z tăng lên 6 . Do Δ1 và Δ2 đều lớn hơn 0 nên vẫn còn khả năng cải thiện hàm mục tiêu khi chuyển sang (hay “xoay sang”) một phương án cực biên kề tốt hơn. 25 Thủ tục xoay (pivotal procedure) Bước 1: Chọn cột xoay là cột bất kỳ có Δj > 0. Lúc đó biến xj tương ứng với cột xoay được chọn làm biến cơ sở mới do xj tăng kéo theo hàm mục tiêu tăng. Ở đây ta chọn đưa x1 vào làm biến cơ sở mới. Bước 2: Chọn hàng xoay để xác định đưa biến nào ra khỏi tập các biến cơ sở (vì tại mỗi bước số biến cơ sở là không thay đổi). Để chọn hàng xoay, ta thực hiện quy tắc “tỷ số dương bé nhất” bằng cách lấy cột phương án (60, 48)T chia tương ứng cho cột xoay (4, 2)T để chọn tỷ số bé nhất. Một điều cần chú ý là ta chỉ xét các tỷ số có mẫu số dương. Vì Min {60/4, 48/2} = 60/4 đạt được tại hàng đầu, nên hàng xoay là hàng đầu (hàng tương ứng với biến x3). Do đó cần đưa x3 ra khỏi tập các biến cơ sở. Bước 3: Chọn phần tử xoay nằm trên giao của hàng xoay và cột xoay. Bước 4: Xoay sang bảng đơn hình mới, xác định các biến cơ sở mới để điền vào cột biến cơ sở, đồng thời thay các giá trị trong cột hệ số hàm mục tiêu. Sau đó, tính lại các phần tử của hàng xoay bằng cách lấy hàng xoay cũ chia cho phần tử xoay để có hàng mới tương ứng. Bước 5: Các phần tử còn lại của bảng đơn hình mới tính theo quy tắc “hình chữ nhật”: (1)mới = (1)cũ – (2)cũ× (4)cũ/(3)cũ, trong đó (3) là đỉnh tương ứng với phần tử xoay . Giải thích. Các bước xoay trên đây chỉ là phép biến đổi tương đương hệ phương trình 4x1 + 2x2 + x3 = 60 (2.1) 2x1 + 4x2 + x4 = 48 (2.2) để có hệ x1 + (1/2)x2 + (1/4)x3 = 15 (2.1’) 0x1 + 3x2 – (1/2)x3 + x4 = 18 (2.2’) bằng cách lấy phương trình (2,1) chia cho 4 (phần tử xoay) để có (2,1’), rồi lấy (2,2) trừ bớt 2× (2.1)/4 để có (2,2’). Đây chính là nội dung của bước 4 và bước 5. Còn việc thực hiện bước 3 sẽ đảm bảo rằng giá trị của các biến cơ sở mới không âm (x1 = 15, x4 = 18). Áp dụng thủ tục xoay cho các phần tử nằm trên hàng 1 và 2 của bảng đơn hình bước 1, sau đó tính các giá trị trên hàng zj và Δj tương tự như khi lập bảng đơn hình bước 1, chúng ta sẽ nhận được bảng đơn hình bước 2. (1) (2) (3) (4) Chẳng hạn: nếu (1)cũ = 4,(2)cũ = 2, (3)cũ = phần tử xoay = 4, (4)cũ = 2 thì (1)mới = 4 – 2×2/4 =3 Quy tắc hình chữ nhật 26 Phân tích bảng đơn hình bước 2 Bảng bước 2 có thể được phân tích tương tự như bảng bước 1. Lúc này giá trị của hàm mục tiêu là z0 = 120 đã được cải thiện hơn so với bước 1. Ta thấy Δ2 = 2 > 0 nên còn có thể cải thiện hàm mục tiêu bằng cách đưa biến x2 vào làm biến cơ sở mới. Thực hiện các bước xoay sang phương án cực biên kề tốt hơn, chúng ta sẽ có bảng đơn hình bước 3. Phân tích bảng đơn hình bước 3 Tại bảng đơn hình bước 3 ta thấy điều kiện tối ưu đã được thoả mãn (Δj ≤ 0, ∀j =1,4 ) nên không còn khả năng cải thiện phương án. Phương án tối ưu đã đạt được tại x1 = 12, x2 = 6, x3 = 0, x4 = 0, tức là tại điểm cực biên B(12, 6) với giá trị zmax = 132. Một số chú ý – Điều kiện tối ưu cho các BTQHTT dạng Max là Δj ≤ 0, ∀j . – Đối với các BTQHTT cần cực tiểu hoá hàm mục tiêu thì điều kiện tối ưu (hay tiêu chuẩn dừng) là Δj ≥ 0, ∀j (nếu ∃ j* sao cho j*Δ < 0 thì cần tiếp tục cải thiện hàm mục tiêu bằng cách chọn cột j* làm cột xoay). – Trong thực tiễn giải các BTQHTT dạng tổng quát có thể xảy ra trường hợp không tìm được phương án xuất phát (tức là không có phương án khả thi). Lúc này có thể kết luận mô hình đã thiết lập có các điều kiện ràng buộc quá chặt chẽ, cần xem xét nới lỏng các điều kiện này. – Trong trường hợp ta tìm được cột xoay mà không tìm được hàng xoay thì kết luận hàm mục tiêu không bị chặn trên (đối với các BTQHTT dạng Max) hoặc không bị chặn dưới (đối với các BTQHTT dạng Min). Trong các trường hợp trên cũng phải dừng lại và kết luận mô hình quy hoạch tuyến tính đã thiết lập không phù hợp với thực tế. Khung thuật toán đơn hình Sau đây là khung thuật toán của phương pháp đơn hình được phát biểu cho BTQHTT cực đại hóa dạng chính tắc. Bước khởi tạo – Tìm một phương án cực biên ban đầu (nếu không tìm được thì dừng và in ra thông báo “BTQHTT không có phương án”). – Tính Δj = cj – zj, ∀j = 1,n , trong đó n là số biến của bài toán đang xét. Các bước lặp Bước 1: Kiểm tra điều kiện tối ưu. Nếu điều kiện tối ưu Δj = cj – zj ≤ 0, ∀j = 1,n đã được thoả mãn thì in / lưu trữ kết quả của bài toán và chuyển sang bước kết thúc. Bước 2: Nếu tồn tại một chỉ số j sao cho Δj > 0 thì tiến hành thủ tục xoay gồm năm bước đã biết, tính lại các Δj, ∀j = 1,n và quay lại bước 1 (Chú ý: Trong trường hợp 27 ta tìm được cột xoay mà không tìm được hàng xoay thì kết luận hàm mục tiêu không bị chặn, in / lưu trữ kết quả của bài toán và chuyển sang bước kết thúc). Bước kết thúc. Dừng. 1.2. Thuật toán đơn hình cho BTQHTT dạng chính tắc Xét BTQHTT sau (xem thêm giáo trình của Nguyễn Hải Thanh, Tối ưu hoá, Nxb Bách Khoa, 2007): z = f(x) = c1x1 + c2x2 + .... + cnxn → Max (Min), với các điều kiện ràng buộc a11x1 + a12x2 + ... + a1nxn ≤ b1 a21x1 + a22x2 + ... + a2nxn ≤ b2 ... am1x1 + am2x2 + ... + amnxn ≤ bm x1, x2, ..., xn ≥ 0 (điều kiện không âm). Đưa BTQHTT trên về dạng chính tắc: z = f(x) = c1x1 + c2x2 + .... + cnxn + 0xn+1 + + 0xn+m → Max (Min), với các điều kiện ràng buộc a11x1 + a12x2 + ... + a1nxn + xn+1 = b1 a21x1 + a22x2 + ... + a2nxn + xn+2 = b2 ... am1x1 + am2x2 + ... + amnxn + xn+m = bm x1, x2, ..., xn+m ≥ 0 (điều kiện không âm). Bước khởi tạo – Nhập các hệ số hàm mục tiêu c, ma trận ràng buộc A và các hệ số vế phải b. – Đặt d1 = cn+1= 0, ..., dm = cn+m= 0 , tức là cB = (d1, ..., dm)T. – Đặt chỉ số các biến cơ sở: r(1) = n + 1, ..., r(m) = n + m. – Gán xr(i) = bi , i = 1,m . Đặt flag = 2. Các bước lặp Bước 1: – Tính cTx = z = d1xr(1) + ... + dmxr(m). – Tính zj = m pj p p 1 a d = ∑ , ∀j = 1, n m+ . – Tìm Δ = [ΔN , ΔB] = [ TNc – TBc B–1N, TBc – TBc B–1B], trong đó ΔB = 0. Như vậy Δj = cj – zj, với zj = m pj p p 1 a d = ∑ , ∀j ∈ N = {1, 2, , n + m} \ {r(1), , r(m)} và Δj = cj – zj = 0, ∀j ∈ B = {r(1), , r(m)}, (tức là zN = TBc B–1N và zB = TBc B–1B). 28 Bước 2: Nếu tồn tại chỉ số j ∈ N sao cho Δj > 0 thì thực hiện thủ tục xoay. – Xác định cột xoay: chọn cột xoay s ứng với một chỉ số j có tính chất Δj > 0. Thông thường chọn j ứng với Δj > 0 lớn nhất, hoặc chọn ngẫu nhiên. – Xác định hàng xoay q theo quy tắc tỷ số dương bé nhất: r (q ) r ( i ) is qs is x x Min , a 0 a a ⎧ ⎫= ∀ >⎨ ⎬⎩ ⎭ . Trong trường hợp không tồn tại ais > 0, đặt flag = 0 và chuyển sang bước kết thúc. – Xác định phần tử xoay aqs. – Tính lại (để chuyển sang bảng đơn hình mới): bq := bq/aqs, aqj := aqj/aqs, ∀j . ∀ i ≠ q tính lại bi := bi – bq*ais và aij := aij – aqj*ais, ∀j. – Đặt lại chỉ số các biến cơ sở: r(q) := s, dq := cs, và xr(i) = bi , i = 1,m . Quay về bước 1. Bước 3: Nếu Δj ≤ 0, ∀j ∈ N thì đặt flag = 1 và chuyển sang bước kết thúc. Bước kết thúc: Ghi lại dữ liệu đầu vào của BTQHTT và kết quả cuối cùng. Nếu flag = 0 thì kết luận BTQHTT có hàm mục tiêu không bị chặn trên. Còn nếu flag = 1 thì kết luận BTQHTT có phương án tối ưu đã tìm được. Dừng. 2. PHƯƠNG PHÁP ĐƠN HÌNH HAI PHA GIẢI BTQHTT TỔNG QUÁT Từ trước tới nay, chúng ta luôn giả sử rằng BTQHTT được xem xét luôn có phương án và có thể biết được một phương án (cực biên) ban đầu của nó để khởi tạo quá trình giải. Trong mục này chúng ta sẽ đi xét các trường hợp khi chưa biết BTQHTT có phương án hay không, cũng như chưa biết được phương án cực biên ban đầu. Đối với những trường hợp này có thể sử dụng phương pháp đơn hình hai pha. Chúng ta sẽ trình bày phương pháp đơn hình hai pha thông qua ví dụ sau. 2.1. Ví dụ . Max z = 8x1 + 6x2, với các ràng buộc 1 2 1 2 1 2 4x 2x 60 2x 4x 48 x ,x 0. + ≤⎧⎪ + ≥⎨⎪ ≥⎩ hay: Max z = 8x1 + 6x2 + 0x3 + 0x4, với các ràng buộc 1 2 3 1 2 4 1 2 3 4 4x 2x x 60 2x 4x x 48 x ,x ,x ,x 0. + + =⎧⎪ + − =⎨⎪ ≥⎩ Trước hết cần trả lời câu hỏi BTQHTT dạng chuẩn tắc trên đây có phương án hay không, nếu có thì cần tìm một phương án cực biên xuất phát của nó. Pha 1. Tìm một phương án cực biên xuất phát bằng cách xét BTQHTT sau đây: Min ω = x5, với các ràng buộc 29 1 2 3 1 2 4 5 1 2 3 4 5 4x 2x x 60 2x 4x x x 48 x ,x ,x ,x ,x 0. + + =⎧⎪ + − + =⎨⎪ ≥⎩ (2.3) Bảng II.2. Các bảng đơn hình giải bài toán pha 1 0 0 0 0 1 Hệ số hàm mục tiêu Biến cơ sở Phương án x1 x2 x3 x4 x5 0 1 x3 x5 60 48 4 2 2 4 1 0 0 –1 0 +1 Hàng ω ω0 = 48 ω1 = 2 ω2 = 4 ω3 = 0 ω4= – 1 ω5 = 1 Hàng Δj Δ1 = –2 Δ2 = – 4 Δ3 = 0 Δ4 = 1 Δ5 = 0 0 0 x3 x2 36 12 3 1/2 0 1 1 0 1/2 –1/4 –1/2 1/4 Hàng ω ω0 = 0 0 0 0 0 0 Hàng Δj 0 0 0 0 1 Mục đích của pha 1 là để giải BTQHTT với các ràng buộc (2.3) hay còn gọi là bài toán ω. Nếu tìm được phương án tối ưu của bài toán ω với các biến giả đều nhận giá trị bằng 0 thì điều này chứng tỏ BTQHTT ban đầu có phương án. Trong trường hợp đó dễ dàng tìm được một phương án cực biên của nó (xem bảng II.2). Tại bảng đơn hình cuối cùng, ta thấy Δj ≤ 0, ∀j, nên phương án tối ưu đã đạt được với x2 = 12, x3 = 36, x1 = x4 = x5 = 0 và ωmin = 0. Do đó chúng ta đưa ra kết luận là BTQHTT ban đầu có phương án x1 = 0, x2 = 12, x3 = 36, x4 = 0. Một cách tổng quát, có thể khẳng định được ngay rằng, nếu bài toán ω có phương án tối ưu với giá trị hàm mục tiêu là 0 thì BTQHTT ban đầu có phương án, trong trường hợp trái lại thì nó không có phương án. Pha 2. Giải BTQHTT ban đầu căn cứ phương án cực biên vừa tìm được ở pha 1 (xem bảng II.3): Max z = 8x1 + 6x2 +0x3 + 0x4, với các ràng buộc 1 2 3 1 2 4 1 2 3 4 4x 2x x 60 2x 4x x 48 x ,x ,x ,x 0. + + =⎧⎪ + − =⎨⎪ ≥⎩ Nhận xét. Kết quả giải ví dụ trên bằng phương pháp đơn hình hai pha cũng giống với kết quả đạt được khi giải bằng phương pháp đơn hình mở rộng. Tuy nhiên, khi sử dụng phương pháp đơn hình hai pha, chúng ta tránh được sự phiền phức trong việc khai báo giá trị dương đủ lớn của tham số M như trong phương pháp đơn hình mở rộng. 30 Bảng II.3. Các bảng đơn hình giải bài toán pha 2 8 6 0 0 Hệ số hàm mục tiêu Biến cơ sở Phương án x1 x2 x3 x4 0 6 x3 x2 36 12 3 1/2 0 1 1 0 1/2 –1/4 Hàng z z0 = 72 z1 = 3 z2 = 6 z3 = 0 z4 =–3/2 Hàng Δj Δ1 = 5 Δ2 = 0 Δ3 = 0 Δ4 = 3/2 0 6 x4 x2 72 30 6 2 0 1 1 1/2 1 0 Hàng z 180 12 6 3 0 Hàng Δj –4 0 –3 0 Tại bảng đơn hình cuối cùng, ta thấy Δj ≤ 0, ∀j, nên phương án tối ưu đã đạt được với x2 = 30, x4 = 72, x1 = x3 = 0 và zmax = 180. 2.2. Thuật toán đơn hình hai pha giải BTQHTT dạng tổng quát Xét bài toán gốc: z = n j 1= ∑ cj xj → Min/ Max n ij j i j 1 n ij j i j 1 n ij j i j 1 j i 1, m 1 i m 1, m m1 1 2 i m m 1, m m m1 2 1 2 3 j 1, n m m m1 2 3 a x b ( ) a x b ( ) a x b ( ) x 0 ( ) = = = = = + + = + + + + = + + + ∑ ∑ ∑ ⎧ ≤⎪⎪ ≥⎪⎪⎨⎪ =⎪⎪ ≥⎪⎩ Bước 1: - Nhập dạng bài toán Min, Max. - Nhập tổng số ràng buộc m bao gồm các ràng buộc mang dấu: ≤ (m1 ràng buộc) , ≥ (m2 ràng buộc) và = (m3 ràng buộc). - Nhập số biến: n biến. - Nhập véc tơ hệ số hàm mục tiêu: C = [ c1, c2, . . ., cn ]. - Nhập véc tơ hệ số vế phải: b = [ b1, b2, . . ., bm ]. - Nhập ma trận hệ số ràng buộc: A = [ai j ]m x n. Bước 2: Đưa bài toán về dạng chính tắc: dạng Max đưa về dạng Min. Đưa thêm biến bù thiếu: m1 biến xn+i, 1i 1,m= , Biến bù thừa: m2 biến xn+ m1+p, 2p 1, m= , Biến giả: m2 + m3 biến xn+m1+m2+q, 2 3q 1,m m= + . Nếu m2 + m3 = 0, chuyển sang bước 4. Nếu m2 + m3 ≠ 0, giải bài toán theo hai pha bằng cách chuyển sang bước 3. 31 Sơ đồ thuật giải (cho BTQHTT gốc là bài toán Max) Bước 3: Pha thứ nhất. Xây dựng và giải bài toán phụ: ω = m2 m3 q +∑ xm1 + m2+q+n → Min Giải pha I Tính Δk Giải pha II Khởi tạo, đổi dấu, thêm biến bù, biến giả Bắt đầu ∀k, Δk≤0 ? Y N Kết thúc Xuất kết quả Tính Max Δk, Δk> 0 Tìm cột Hàm mục tiêu không bị chặn, BT vô nghiệm Tìm phần tử Chuyển đổi cơ sở N Y Tìm hàng xoay Tính f( ) BT có nghiệm Tính Δk Y N Có biến giả ? ∃k,Δk<0? Y N Tính Max kΔ , Δk Tìm cột xoay Tìm hàng xoay và phần tử Chuyển đổi cơ sở BT vô nghiệm ∃biến giả≠ 0? Y N Xóa biến giả Pha I Pha II Sơ đồ thuật giải đơn hình hai pha 32 1 2 3m 2*m m ij j i 1 2 3 j 1 j 1 2 3 a x b (i 1, m m m ) x 0 j 1, (n m 2* m m ). + + = ⎧⎪ = = + +⎪⎨⎪ ≥ ∀ = + + +⎪⎩ ∑ Kết thúc pha 1: Xảy ra ba trường hợp sau: – Phương án tối ưu không có biến giả, lấy đó làm phương án xuất phát, thay lại hệ số hàm mục tiêu, loại trừ các cột biến giả và sang bước 4. – Phương án tối ưu có biến giả khác 0 thì bài toán tối ưu không có phương án. Dừng. – Phương án tối ưu có chứa biến giả nhưng biến giả bằng 0, xoá các dòng chứa các biến giả này, thay lại hệ số hàm mục tiêu, loại trừ các cột biến giả và sang bước 4. Bước 4: Pha thứ 2. Giải bài toán gốc với phương án xuất phát tìm được bằng phương pháp đơn hình. Bước 5: In kết quả . 3. GIẢI CÁC BÀI TOÁN TỐI ƯU TRÊN MICROSOFT EXCEL Microsoft Excel 2000, 2003 có các công cụ toán học rất mạnh để giải các bài toán tối ưu và thống kê toán học. Excel có thể giải được các loại bài toán tối ưu: BTQHTT tổng quát, các biến có thể có ràng buộc hai phía, ràng buộc cũng có thể viết ở dạng hai phía; bài toán vận tải có hai chỉ số; bài toán quy hoạch nguyên (các biến có điều kiện nguyên hay boolean); bài toán quy hoạch phi tuyến. Số biến của BTQHTT hay nguyên có thể lên tới 200 biến. Excel còn có thể giải các bài toán hồi quy trong thống kê toán học: hồi quy đơn, hồi quy bội, hồi quy mũ. Dùng Solver ta có thể tìm cực đại hay cực tiểu của một hàm số đặt trong một ô gọi là ô đích. Solver chỉnh sửa một nhóm các ô (gọi là các ô có thể chỉnh sửa) có liên quan trực tiếp hay gián tiếp đến công thức nằm trong ô đích để tạo ra kết quả. Ta có thể thêm vào các ràng buộc để hạn chế các giá trị mà Solver có thể dùng. Đối với BTQHTT Solver dùng phương pháp đơn hình, đối với quy hoạch phi tuyến Solver dùng phương pháp tụt gradient để tìm một cực trị địa phương. 3.1. Giải BTQHTT Xét bài toán quy hoạch c1x1 + c2 x2 + + cnxn = f(x) → max / min a11x1 + a12 x2 + + a1nxn Q b1 a21x1 + a22 x2 + " + a2nxn Q b2 am1x1 + am2 x2 + + amnxn Q bm xj ≥ 0, j = 1, . . . , n nguyên hoặc nhị phân 0–1. 33 Trong đó Q là một trong các phép toán quan hệ ≥ , ≤ hoặc = , thứ tự các phép toán quan hệ trong các ràng buộc là tuỳ ý. Như vậy bài toán trên có thể là BTQHTT thông thường, quy hoạch tuyến tính nguyên hay quy hoạch 0–1. Cách bố trí dữ liệu cho trên bảng tính: c[1] c[2] . . . . . . c[n] Σ c[j] x[j] a[1,1] a[1,2] . . . . . . a[1,n] Σ a[1,j] x[j] b[1] a[2,1] a[2,2] . . . . . . a[2,n] Σ a[2,j] x[j] b[2] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . a[m,1] a[m,2] . . . . . . a[m,n] Σ a[m,j] x[j] b[m] x[1] x[2] . . . . . . x[n] Hàng cuối cùng là các giá trị ban đầu của các biến để các công thức của Excel hoạt động, có thể lấy giá trị của tất cả các biến bằng 1. Xét bài toán: x1 +4x2 + x3 → Min , với các ràng buộc: 2x1 +3x2 +4x3 ≥ 20 5x1 – x2 +2x3 ≥12 x1 +2x2 – x3 ≤ 2 x1 +4x2 –2x3 ≤1 x1, x2, x3 ≥ 0 Các bước thực hiện để giải bài toán: Bước 1. Nhập dữ liệu bài toán vào bảng tính dưới dạng sau: Phương án ban đầu X = (1, 1, 1), nó có thể không chấp nhận được. Bước 2. Tính giá trị hàm mục tiêu tại ô E2 bằng công thức = SUMPRODOCT($B$7 : $D$7, B2 : D2) Hàm Sumproduct cho tích vô hướng của hai dãy ô. Copy công thức từ ô E2 sang dãy các ô E3 : E6 nhằm tính giá trị vế trái của bốn ràng buộc bài toán (1). Bước 3. Dùng lệnh Tools / Solver, xuất hiện hộp thoại Solver Parameters. 34 Mục Set Target Cell: chọn ô đích (chứa giá trị hàm mục tiêu), có thể nháy vào biểu tượng của Excel bên phải hộp văn bản để xác định ô, trong ví dụ chọn ô E2. Mục Equal To: chọn Max nếu cực đại hàm mục tiêu, chọn Min nếu cực tiểu hàm mục tiêu, chọn Value of và nhập giá trị nếu muốn ô đích bằng một giá trị nhất định, trong ví dụ chọn Min. Mục By Changing cells: chọn các ô chứa các biến của bài toán, ta chọn khối ô B7:D7. Nháy nút Add để nhập tất cả các ràng buộc vào khung Subject to the Constraints (dòng đầu trong khung ứng với ràng buộc không âm trên các biến, dòng thứ hai ứng với hai ràng buộc đầu bài toán (2), dòng cuối ứng với 2 ràng buộc cuối). Khi nháy nút Add, hiện hộp thoại Hộp văn bản Cell Reference để chọn các ô cần đặt ràng buộc lên chúng, hộp văn bản ở giữa để chọn loại ràng buộc (>= = <= interger, binary), hộp văn bản Constraint để chọn giá trị ràng buộc (có thể là số hay giá trị trong các ô). Sau khi nhập xong các ràng buộc, nháy vào nút Options, hiện hộp thoại Solver Options, đánh dấu kiểm vào mục Assume Linear Model (khảng định mô hình của ta là tuyến tính). Bước 4. Trong hộp thoại Solver Parameters nháy vào nút Solve để bắt đầu giải bài toán tối ưu. Giải xong bài toán xuất hiện hộp thoại Solver Results, chọn mục Keep Solver Solution (giữ lại lời giải), nháy OK, kết quả giải bài toán nằm ở các ô B7 : D7. Kết quả ta được phương án tối ưu là X = (0.5 , 0 , 4.75), giá trị cực tiểu hàm mục tiêu là 5.25 ở ô E2. 35 3.2. Giải bài toán quy hoạch phi tuyến có ràng buộc tuyến tính Xét bài toán quy hoạch phi tuyến Min{f(x)| gi (x) = 0, i = 1, 2, , m, x∈Rn}. Để giải bài toán quy hoạch phi tuyến bằng Solver ta cần xác định khối ô để chứa các biến (x[1], x[2], . . . , x[n]), một ô chứa giá trị hàm mục tiêu f(x), khối m ô chứa giá trị các hàm gi (x) . Ví dụ. Giải bài toán quy hoạch toàn phương: – x1 –2x2 +0,5x1 2 +0,5x2 2 →Min 2x1 +3x2 + x3 =6 x1 +4x2 + x4 =5 x1, x2, x3, x4 ≥ 0 Bảng tính để giải bài toán này như sau: Phương án trong khối ô B2:E2 (phương án ban đầu cho mọi phần tử bằng 0), hàm mục tiêu trong ô F2 xác định bởi công thức = - b2 - 2*c2 + 0.5*b2^2 + 0.5*c2^2. Ô F3 tính theo công thức = sumproduct ($b$2: $e$2, b3 : e3), công thức này chép sang ô F4. Các ràng buộc bài toán B2 : E2 >= 0, và F3:F4 = G3:G4. Khi giải các bài toán quy hoạch phi tuyến ta phải bỏ chọn mục Assume Linear Model trong hộp thoại Solver Options. Kết quả dùng Solver giải bài toán: trị tối ưu -2.0294, phương án tối ưu (0.7647, 1.0588, 1.294, 0). Tóm lại Solver có thể giải được nhiều bài toán tối ưu, số lượng biến tối đa của bài toán là 200 biến. Tuy nhiên cũng có nhiều bài toán nó không giải được, khi đó nó đưa thường đưa ra các thông báo: – Solver could not find a feasible solution: bài toán không có lời giải chấp nhận được. Hoặc có thể do các giá trị khởi đầu của những ô chứa biến khác quá xa các giá trị tối ưu, hãy thay đổi các giá trị khởi đầu và giải lại. – The maximum iteration limit was reached, continue anyway ? số bước lặp đã đến số cực đại. Ta có thể tăng số bước lặp ngầm định nhờ lệnh Tools/ Solver, chọn Options, nhập giá trị mới vào hộp Iterations. – The maximum time limit was reached, continue anyway ? thời gian chạy vượt quá thời gian tối đa ngầm định. Ta có thể sửa giá trị trong mục Max Time trong gộp thoại Solver Options. Chú ý, nếu các lệnh Solver và Data Analysis không có trong menu Tools ta phải cài đặt bổ sung từ đĩa CD: dùng lệnh Tools / Add-Ins, hiện hộp thoại, chọn mục Solver Add in và Analysis ToolPak. 36 3.3. Một số ví dụ khác Bài 1. Giải BTQHTT nguyên bộ phận: z=5x1 +x2 +x3 +2x4 +3x5 →min –x2 +5x3 –x4 –2x5 ≤ 2 5x1 –x2 +x5 ≥7 x1 +x2 +6x3 +x4 ≥4 xj ≥ 0 j=1, 2, 3, 4, 5 xj = interger, j = 1, 2, 3. Đáp số: Giá trị tối ưu là 12, phương án tối ưu (2, 2, 0, 0, 0). Bài 2. Giải BTQHTT 0–1 (bài toán cái túi) sau: 30x1 +19x2 +13x3 +38x4 +20x5 +6x6 +8x7 +19x8 +10x9 +11x10 →max 15x1 +12x2 +9x3 +27x4 +15x5 +5x6 +8x7 +20x8 +12x9 +15x10 → 62 xj ∈{0, 1}, j=1, 2,",10 Đáp số: Giá trị tối ưu là 95, phương án tối ưu là ( 1, 1, 0, 1, 0, 0, 1, 0, 0, 0). Bài 3. Giải bài toán quy hoạch lõm (có thể có nhiều cực tiểu địa phương) –x1 2 +2x1 –x2 2 +4x2 –x3 2 +8x3 –x4 2 +14x4 –x5 2 +18x5 –180 →Min –x1 –2x2 +x3 +2x4 +3x5 ≤ 85 –7x1 +9x2 –5x3 +33x4 –11x5 ≤500 2x1 –x2 +2x3 –x4 +2x5 ≤ 150 1.3x1 +x2 +x3 +x4 +x5 ≤ 300 x1 +x2 +x3 +x4 +x5 ≤ 300 x1, x2, x3, x4, x5 ≥0 Đáp số. Với phương án ban đầu X = (50, 50, 50, 50, 50) dùng Solver có phương án tối ưu là X = (0, 190, 0, 0, 110) và trị tối ưu hàm mục tiêu là - 45640. 4. GIẢI BTQHTT TRONG LINGO LINGO cho phép giải rất nhiều loại toán tối ưu, trong đó có BTQHTT (biến liên tục cũng như biến nguyên). Để giải bài toán này, chúng ta cần cài đặt Lingo vào trong máy tính. Nhấn vào biểu tượng Lingo trên màn hình để vào cửa sổ Lingo. Sau đó thực hiện các lệnh Lingo: Menu > New > và gõ vào các dữ liệu của bài toán. Nhập bài toán max = 8*x1+6*x2; 4*x1+2*x2<=76; 37 2*x1+5*x2<=52; @gin(x1); @gin(x2); Hai điều kiện sau cùng là các điều kiện biến nguyên. Kết quả chạy bài toán khi các biến đều liên tục Rows= 3 Vars= 2 No. integer vars= 0 ( all are linear) Nonzeros= 8 Constraint nonz= 4( 0 are +- 1) Density=0.889 Smallest and largest elements in absolute value= 2.00000 76.0000 No. : 0, Obj=MAX, GUBs <= 1 Single cols= 0 Optimal solution found at step: 0 Objective value: 159.0000 Variable Value Reduced Cost X1 17.25000 0.0000000E+00 X2 3.500000 0.0000000E+00 Row Slack or Surplus Dual Price 1 159.0000 1.000000 2 0.0000000E+00 1.750000 3 0.0000000E+00 0.5000000 Kết quả chạy bài toán khi biến x1 nguyên Rows= 3 Vars= 2 No. integer vars= 1 ( all are linear) Nonzeros= 8 Constraint nonz= 4( 0 are +- 1) Density=0.889 Smallest and largest elements in absolute value= 2.00000 76.0000 No. : 0, Obj=MAX, GUBs <= 1 Single cols= 0 Optimal solution found at step: 5 Objective value: 157.6000 Branch count: 1 Variable Value Reduced Cost X1 17.00000 -5.600000 X2 3.600000 0.0000000E+00 Row Slack or Surplus Dual Price 1 157.6000 1.000000 2 0.8000000 0.0000000E+00 3 0.0000000E+00 1.200000 38 Kết quả chạy bài toán khi các biến đều nguyên Rows= 3 Vars= 2 No. integer vars= 2 ( all are linear) Nonzeros= 8 Constraint nonz= 4( 0 are +- 1) Density=0.889 Smallest and largest elements in absolute value= 2.00000 76.0000 No. : 0, Obj=MAX, GUBs <= 1 Single cols= 0 Optimal solution found at step: 7 Objective value: 156.0000 Branch count: 2 Variable Value Reduced Cost X1 18.00000 -8.000000 X2 2.000000 -6.000000 Row Slack or Surplus Dual Price 1 156.0000 1.000000 2 0.0000000E+00 0.0000000E+00 3 6.000000 0.0000000E+00 5. GIẢI BTQHTT BẰNG PHẦN MỀM QHTT Sử dụng phần mềm QHTT trên mạng giáo dục edu.net.vn, dễ dàng nhập được dữ liệu BTQHTT và có đáp số với toàn bộ các bảng trung gian. Tuy nhiên phần mềm này chỉ áp dụng cho các biến liên tục và vẫn còn sai sót. Bài toán dạng chính tắc: F(x) = 8x1 + 6x2 => MAX Các ràng buộc: 4x1 + 2x2 + x3 = 60 2x1 + 4x2 + x4 = 48 Trong đó: 39 x3, x4 là biến bù x1 >=0, x2 >=0, x3 >=0, x4 >=0 Ci Xi Yi X1 X2 X3 X4 Lamda 0 X3 60 4 2 1 0 15 0 X4 48 2 4 0 1 24 F(x) 0 -8 -6 0 0 Ci Xi Yi X1 X2 X3 X4 Lamda 8 X1 15 1 1/2 1/4 0 30 0 X4 18 0 3 -1/2 1 6 F(x) 120 0 -2 2 0 D Ci Xi Yi X1 X2 X3 X4 Lamda 8 X1 12 1 0 1/3 -1/6 - 6 X2 6 0 1 -1/6 1/3 - F(x) 132 0 0 5/3 2/3 Phương án tối ưu của bài toán là : (12,6,0,0) Giá trị tối ưu của hàm mục tiêu là : F(x) = 132 40 Chương III BÀI TOÁN QUY HOẠCH PHI TUYẾN 1. PHƯƠNG PHÁP RST2ANU GIẢI BÀI TOÁN TỐI ƯU PHI TUYẾN TOÀN CỤC HỖN HỢP NGUYÊN 1.1. Đặt vấn đề Dạng chính tắc của bài toán tối ưu một mục tiêu được biểu diễn như sau: Min (Max) f(X) , X = (x1, x2, , xn)∈ Rn, với các điều kiện ràng buộc (i) gj(X) ≤ 0, j = 1, 2, , k, (ii) gj(X) = 0, j = k+1, k+2, , .m. Trong các bài toán thực tế có thể bổ sung thêm các ràng buộc (iii) ai ≤ xi ≤ bi, i = 1, 2, , n. Trong trường hợp hàm mục tiêu f(X) hay có ít nhất một trong các hàm ràng buộc gj(X), j = 1, 2, , m, là hàm phi tuyến, chúng ta có bài toán tối ưu phi tuyến. Khi tất cả các toạ độ xi đều bắt buộc nhận các giá trị nguyên, i = 1, 2, , n, thì ta có bài toán tối ưu nguyên. Còn nếu chỉ có một số toạ độ (nhưng không phải tất cả các toạ độ) bắt buộc nhận giá trị nguyên thì ta có bài toán tối ưu hỗn hợp nguyên. Ký hiệu D là miền các phương án (miền ràng buộc) cho bởi các ràng buộc (i), (ii) và / hoặc (iii) thì bài toán tối ưu trên đây có thể viết gọn hơn như sau: f(X) → Min (Max) với X ∈ D. Lúc này, đối với bài toán cực tiểu hoá, X* ∈ D được gọi là phương án tối ưu toàn cục nếu ∀ X∈D ta luôn có: f(X*) ≤ f(X). Trong trường hợp f(X*) ≤ f(X) chỉ đúng với ∀X∈D trong một lân cận nào đó của X* thì X* được gọi là phương án tối ưu địa phương. Một cách tương tự, ta có thể định nghĩa khái niệm phương án tối ưu toàn cục / địa phương cho bài toán cực đại hoá. Nếu chúng ta chỉ quan tâm tới việc tìm kiếm phương án tối ưu toàn cục thì ta có bài toán tối ưu toàn cục. Các phương pháp giải bài toán tối ưu toàn cục phi tuyến đơn mục tiêu được phân ra thành hai lớp: phương pháp tất định và phương pháp ngẫu nhiên (deterministic and stochastic methods). Phương pháp tất định sử dụng các tính chất giải tích của hàm mục tiêu và các hàm ràng buộc. Một số dạng bài toán tối ưu toàn cục với những tính chất giải tích nhất định của hàm mục tiêu và các hàm ràng buộc có thể giải được bằng các phương pháp tất định thích hợp, chẳng hạn như phương pháp quy hoạch toàn phương, quy hoạch tách, quy hoạch lồi, quy hoạch d.c Trong các trường hợp đó phương án tối ưu toàn cục có thể tìm được sau một số hữu hạn bước tính toán với độ chính xác chọn trước. Tuy nhiên, đối với nhiều lớp bài toán tối ưu toàn cục phương pháp tất định tỏ ra không có hiệu quả. Trong khi đó, các phương pháp ngẫu nhiên như: phương pháp đa khởi tạo (multistart), mô phỏng tôi (simulated annealing), thuật giải di truyền (genetic 41 algorithm) có thể áp dụng để giải các bài toán tối ưu toàn cục dạng bất kỳ, không đòi hỏi các tính chất đặc biệt của hàm mục tiêu hay các hàm ràng buộc. Các phương pháp ngẫu nhiên đặc biệt tỏ ra có hiệu quả đối với các bài toán tối ưu phi tuyến nguyên và hỗn hợp nguyên. Tuy nhiên, các phương pháp này thường chỉ cho phương án “gần” tối ưu khá tốt sau một số hữu hạn bước mà không kiểm soát được độ chính xác của phương án tìm được. Như vậy, hiện tại có nhiều phương pháp tối ưu toàn cục được đề xuất. Tuy nhiên chưa có một phương pháp nào tỏ ra hữu hiệu cho mọi bài toán tối ưu, đặc biệt là các bài toán tối ưu với biến nguyên hay hỗn hợp nguyên. Hơn nữa, các phương pháp tối ưu cần phải được lập trình để đóng gói thành các phần mềm thân thiện đối với người sử dụng. Đây là một đòi hỏi rất thực tế của các kĩ sư, các nhà khoa học, các doanh nghiệp trong nhiều lĩnh vực công nghiệp cũng như nông nghiệp. Trong bài báo này, chúng tôi trình bày một phần mềm tính toán khoa học (RST2ANU) có thể đáp ứng được phần nào các đòi hỏi nêu trên đối với người sử dụng để giải các bài toán tối ưu phi tuyến toàn cục với các biến liên tục, nguyên hoặc hỗn hợp nguyên. Phần mềm này được xây dựng dựa trên phương pháp tìm kiếm ngẫu nhiên có kiểm soát (cùng tên gọi RST2ANU) do Mohan và Nguyễn Hải Thanh đề xuất (Xem “A controlled random search technique incorporating the simulated annealing concept for solving integer and mixed integer global optimization problems”, Computational Optimization and Applications, Vol. 14, pp. 103-132, 1999). Đây là một phương pháp tối ưu đã được chạy kiểm thử trên hàng trăm bài toán mẫu và nhiều bài toán thực tế với độ tin cậy rất cao và tốc độ tính toán chấp nhận được. 1.2. Thuật giải tìm kiếm ngẫu nhiên có kiểm soát RST2ANU Thuật giải RST2ANU là thuật giải lặp, bao gồm hai pha, pha địa phương và pha toàn cục. Trong pha toàn cục, một số lượng thích hợp đủ lớn các phương án chấp nhận được được phát sinh ra một cách ngẫu nhiên và lưu trữ trong mảng có tên A. Đánh dấu hai điểm có giá trị hàm mục tiêu lớn nhất và nhỏ nhất tương ứng là M và L. Trong pha địa phương, các phương án được xử lý nhằm thu được giá trị ước lượng tốt hơn của hàm mục tiêu. Trong pha này, thuật giải xác định X* là điểm được nội suy bậc hai dựa trên phương án L và hai phương án khác được chọn ngẫu nhiên trong mảng A. Nếu như X* chấp nhận được thì với f(X*) ≤ f(M), M sẽ được thay thế bởi X* trong mảng A, còn với f(X*)>f(M) M sẽ được thay thế bởi X* với xác suất p= exp(-β(f(X*)-f(M))/(f(X*)-f(L))) , trong đó β >0 là tham số được lựa chọn thích hợp. Nếu X* không phải là phương án chấp nhận được, bỏ qua X* và chọn hai phương án khác trong A một cách ngẫu nhiên rồi cùng với L tiếp tục sinh ra phương án mới. Quá trình cứ thế tiếp diễn như vậy cho tới khi tập hợp các phương án trong A sẽ có xu hướng co cụm lại xung quanh một phương án tối ưu toàn cục. Lưu đồ thuật giải RST2AN được thể hiện trên hình minh hoạ, trong đó: • n, f(X), g(j), ai, bi, là các đầu vào. 42 • A = RandomNSolution (N) phát sinh N phương án ngẫu nhiên chấp nhận được, đồng thời tính giá trị của hàm mục tiêu và trả về kết quả cho mảng A. Như vậy, mảng A chứa luôn cả giá trị hàm mục tiêu tương ứng với từng phương án. Lưu đồ thuật giải RST2ANU • Arrange(A) sắp xếp mảng A theo thứ tự tăng dần của hàm mục tiêu. r = random(0,1) p=exp(-beta(f(X*)-f(M))/(f(X*)-f(L))) N 43 • Max(A), Min(A) trả về phương án có giá trị hàm mục tiêu lớn nhất và nhỏ nhất trong A. • Clustered(A, eps1, eps2) cho biết mảng A đã hội tụ theo hàm mục tiêu hay chưa. • Nếu (f(M) – f(L))/FM) < eps1 thì mảng A hội tụ, ngược lại chưa hội tụ, với FM = f(M) nếu f(M) > eps2, ngược lại FM = 1. • NewSolution() trả về một phương án mới được suy ra từ 3 điểm: L và hai điểm được chọn ngẫu nhiên khác trong mảng A theo phương pháp nội suy. • Feas(X) nhận giá trị TRUE nếu X chấp nhận được, ngược lại nhận giá trị FALSE • Random(0,1) trả về giá trị ngẫu nhiên nằm trong khoảng (0,1). • Replace(A, M, X*) thay thế M trong A bởi X* kèm theo cả giá trị hàm mục tiêu sao cho không cần phải sắp xếp lại mảng A mà vẫn đảm bảo các điểm được sắp xếp theo thứ tự giá trị hàm mục tiêu tăng dần. • Kết thúc 1: Số lần tìm kiếm liên tiếp mà không cải thiện được giá trị hàm mục tiêu vượt quá số lần cho phép. Thuật giải dừng với giá trị tốt nhất của hàm mục tiêu tìm được là FL tương ứng với phương án L. • Kết thúc 2: Phương án tối ưu toàn cục đã đạt được là L với giá trị hàm mục tiêu là FL. • Kết thúc 3: Số lần nội suy liên tiếp mà không tìm được phương án thay thế M trong A vượt quá số lần cho phép. Thuật giải dừng với giá trị tốt nhất của hàm mục tiêu tìm được là FL tương ứng với phương án L. • Kết thúc 4: Số lần lặp vượt quá số lần cho phép. Thuật giải dừng với giá trị tốt nhất của hàm mục tiêu tìm được là FL tương ứng với phương án L. 1.3. Một số nhận xét về phiên bản nâng cấp của phần mềm Phần mềm tính toán khoa học RST2ANU đã được thiết kế và xây dựng có thể sử dụng để giải quyết nhiều mô hình tối ưu phát sinh trong lĩnh vực nông nghiệp, hỗ trợ cho giảng dạy và nghiên cứu khoa học nông nghiệp cũng như trong các lĩnh vực khác. Phần mềm này đã được nâng cấp có tính thân thiện với người sử dụng và tránh được sao chép, có thể được phổ cập có bản quyền một cách rộng rãi. Việc tạo ra các giao diện thân thiện cho phép dễ dàng nhập các hàm mục tiêu và ràng buộc của nhiều dạng bài toán tối ưu phi tuyến là một vấn đề khá phức tạp đã được giải quyết thành công trong phần mềm này. Trong tình hình hiện tại, khi các phần mềm tối ưu phi tuyến không có sẵn trên thị trường trong và ngoài nước, phần mềm RST2ANU nên được triển khai sử dụng để giải quyết các bài toán tối ưu trong các lĩnh vực khác nhau, bao gồm cả các bài toán nguyên và hỗn hợp nguyên. Các nghiên cứu cần tiếp tục được triển khai và được hỗ trợ 44 về mặt tài chính để tích hợp RST2ANU vào các gói phần mềm trong điều khiển tự động hóa hay trong các hệ hỗ trợ ra quyết định. 2. MỘT SỐ VÍ DỤ ÁP DỤNG RST2ANU 2.1. Bài 1: Bài toán xác định tham số sàng phân loại Sau đây là cách viết chương trình con tính giá trị hàm mục tiêu và kiểm tra tính khả thi của phương án đã được phát sinh khi thực hiện chương trình máy tính RST2ANU. float f() { float fv, g; g= 0.05*cos(0+3.1417/18)+0.3*cos(x[0])+0.15*cos(x[1])+0.5*cos(x[2])- 0.365; g=1000*g*g; fv=g; g= 0.05*sin(0+3.1417/18)+0.3*sin(x[0])+0.15*sin(x[1])+0.5*sin(x[2])-0.635; g=1000*g*g; fv=fv+g; g=0.05*cos(0+3.1417/18)+0.3*cos(x[0])+1.075*cos(x[1]- 3.1417/8)+0.4*cos(x[3])-1.365; g=1000*g*g; fv=fv+g; g=0.05*sin(0+3.1417/18)+0.3*sin(x[0])+1.075*sin(x[1]- 3.1417/8)+0.4*sin(x[3])-0.635; g=1000*g*g; fv=fv+g; return(fv); } int feas() {int flag=1; return(1); } File vào v1.in 1 4 0 4 0.000001 0.000001 10000 2000 0 2000 0 0 1 0 0 0 0 0 1 2 3 0 3.1417 0 3.1417 0 3.1417 0 3.1417 123 234 345 456 567 45 3 3 0.01 0 0 File ra v1.out PROBLEM No 1 **crs2.c n=4,nc=0,nint=4,epsilon=0.000001,eps1=0.000001,iterlast=10000 ifailast=2000,imlast=0,islast=2000,iprint=0,id=0,imp=1,ipat=0 noninteger variables are: x[0] x[1] x[2] x[3] guess=0,nguess=0,lppatt=0 lower and upper bounds of coordinates: xmin[0]=0.000000,xmax[0]=3.141700 xmin[1]=0.000000,xmax[1]=3.141700 xmin[2]=0.000000,xmax[2]=3.141700 xmin[3]=0.000000,xmax[3]=3.141700 itnlast=3,istnlast=3,eps2=0.010000 ***CRS SOLUTION FOR NLPP by type 1 as usually **case 1 na=50 seed[0]=*, s*,f=0.000000,fm=0.000001,iter=763,ifun=1664,t1=0 *, s*,f=0.000000,fm=0.000001,iter=1069,ifun=1619,t1=0 *, s*,f=0.000000,fm=0.000001,iter=1026,ifun=1557,t1=0 *,f*,f=4.300038,fm=4.300157,iter=1941,ifun=25813,t1=0 t*, f= 0.000000,iter2=4799 ifun2=30653,t2=0 fopt=0.000000 0.198735,0.550661,1.784750,1.670850, seed[1]=*,f*,f=4.300035,fm=4.300137,iter=1891,ifun=32299,t1=0 *,f*,f=4.300043,fm=4.300156,iter=1705,ifun=31407,t1=0 *,f*,f=4.300050,fm=4.300186,iter=967,ifun=21132,t1=1 *, s*,f=0.000000,fm=0.000001,iter=938,ifun=1481,t1=0 t*, f= 0.000000,iter2=5501 ifun2=20783,t2=1 fopt=0.000000 0.198752,0.550653,1.784751,1.670846, seed[2]=*,f*,f=4.300051,fm=4.300158,iter=1168,ifun=-27879,t1=0 *, s*,f=0.000000,fm=0.000001,iter=933,ifun=1398,t1=0 *,f*,f=4.300040,fm=4.300148,iter=1295,ifun=-32336,t1=0 *, s*,f=0.000000,fm=0.000001,iter=1063,ifun=1643,t1=0 t*, f= 0.000000,iter2=4459 ifun2=8362,t2=0 fopt=0.000000 0.198742,0.550658,1.784749,1.670847, seed[3]=*, s*,f=0.000000,fm=0.000001,iter=943,ifun=1430,t1=0 46 *, s*,f=0.000000,fm=0.000001,iter=868,ifun=1345,t1=0 *, s*,f=0.000000,fm=0.000001,iter=1048,ifun=1634,t1=0 *,f*,f=4.300054,fm=4.300150,iter=2354,ifun=-30714,t1=0 t*, f= 0.000000,iter2=5213 ifun2=-26305,t2=0 fopt=0.000000 0.198741,0.550658,1.784751,1.670850, seed[4]=*, s*,f=0.000000,fm=0.000001,iter=1216,ifun=1905,t1=0 *, s*,f=0.000000,fm=0.000001,iter=1191,ifun=1859,t1=0 *,f*,f=4.300034,fm=4.300148,iter=1659,ifun=27170,t1=0 *, s*,f=0.000000,fm=0.000001,iter=801,ifun=1234,t1=0 t*, f= 0.000000,iter2=4867 ifun2=32168,t2=0 fopt=0.000000 0.198742,0.550658,1.784751,1.670848, ***b1=1196,b2=1883,*** ***a1=4967,a2=25,*** 2.2. Bài 2: Bài toán xác định cơ cấu đầu tư chăn nuôi cá float f() { float fv; fv=19.375*pow(x[0],0.236)*pow(x[1],0.104)*pow(x[2],0.096)*pow(x[3],0.056); fv=fv*pow(x[4],0.056)*exp(0.168*x[5])*exp(0.066*x[6]); return(fv); } int feas() { int flag=1; float g; g=x[0]+x[1]+x[2]+x[3]+x[4]; if(g>50) { flag=0; goto LAST; } g=x[5]+x[6]; if(g>1) { flag=0; goto LAST; } LAST: if(flag==0) {return(0);} 47 else {return(1);} File vào CUONG1.IN 1 7 2 5 0.001 0.001 10000 2000 0 2000 0 0 0 0 0 0 0 0 1 2 3 4 0 50 0 50 0 50 0 50 0 50 0 1 0 1 123 234 345 456 567 0 0 0.01 0 0 File ra CUONG1.OUT PROBLEM No 1 **crs2.c n=7,nc=2,nint=5,epsilon=0.001000,eps1=0.001000,iterlast=10000 ifailast=2000,imlast=0,islast=2000,iprint=0,id=0,imp=0,ipat=0 noninteger variables are: x[0] x[1] x[2] x[3] x[4] guess=0,nguess=0,lppatt=0 lower and upper bounds of coordinates: xmin[0]=0.000000,xmax[0]=50.000000 xmin[1]=0.000000,xmax[1]=50.000000 xmin[2]=0.000000,xmax[2]=50.000000 xmin[3]=0.000000,xmax[3]=50.000000 xmin[4]=0.000000,xmax[4]=50.000000 xmin[5]=0.000000,xmax[5]=1.000000 xmin[6]=0.000000,xmax[6]=1.000000 itnlast=0,istnlast=0,eps2=0.010000 ***CRS *** SOLUTION FOR NLPP by type 1 as usually **case 1 na=16 seed[0]=0*, s*,f=-88.334068,fm=-88.248016,iter=81,ifun=15061,t1=0 t*, f= -88.334068,iter2=81 ifun2=15061,t2=0 fopt=-88.334068 20.970308,9.272722,9.084186,5.446885,5.225899,1.000000,0.000000, seed[1]=0*, s*,f=-88.360970,fm=-88.274643,iter=688,ifun=15912,t1=0 t*, f= -88.360970,iter2=688 ifun2=15912,t2=0 fopt=-88.360970 21.514578,9.515970,8.769965,5.103930,5.095558,1.000000,0.000000, seed[2]=0*, s*,f=-88.360855,fm=-88.272774,iter=588,ifun=14737,t1=1 t*, f= -88.360855,iter2=588 ifun2=14737,t2=1 fopt=-88.360855 48 21.570715,9.461304,8.736217,5.097508,5.134254,1.000000,0.000000, seed[3]=0*, s*,f=-88.359039,fm=-88.270882,iter=691,ifun=17550,t1=0 t*, f= -88.359039,iter2=691 ifun2=17550,t2=0 fopt=-88.359039 21.375309,9.651163,8.772411,5.121580,5.079536,1.000000,0.000000, seed[4]=0*, s*,f=-88.360451,fm=-88.273888,iter=691,ifun=16477,t1=0 t*, f= -88.360451,iter2=691 ifun2=16477,t2=0 fopt=-88.360451 21.518408,9.449259,8.761295,5.180042,5.090996,1.000000,0.000000, ***t3=0/1*** ***b1=547,b2=2840,*** ***a1=547,a2=2840,*** 3. TÍCH HỢP RST2ANU VỚI MATLAB Trong Matlab không có lệnh nhảy vô điều kiện goto. Để xử lí vấn đề này cần lập trình theo từng khối tuần tự, kết hợp với cách sử dụng các lời gọi hàm thông qua các file chương trình, khối lệnh khác nhau. Chương trình RST2ANU gồm nhiều file, mỗi file có thể là một chương trình con, cũng có thể là một khối lệnh. Chương trình hoàn chỉnh được chứa trong thư mục RST2ANU, file CRS là file chính và cũng là lệnh gọi chương trình từ Matlab. Nhập / lưu bài toán và các dữ kiện Chạy chương trình bằng cách dùng lệnh crs từ dòng lệnh của Matlab. Khi chạy sẽ thấy xuất hiện: 1. Enter the problem. 2. Load problem from file. 3. Exit. Your choice: Nhập dữ liệu cho bài toán có thể nhập từ file dữ liệu có sẵn (chọn 2), hoặc có thể nhập lần lượt các tham số đầu vào (chọn 1), nếu chọn 3 là thoát khỏi chương trình. 3.1. Nhập từ bàn phím Ví dụ. Xét bài toán tối ưu phi tuyến toàn cục hỗn hợp nguyên z = x10,6 + x20,6 + x30,4 + 2x4 + 5x5 − 4x3 – x6, → Min với các ràng buộc: x2 − 3x1 − 3x4 = 0; x3 − 2x2 − 2x5 = 0; 4x4 – x6 = 0; x1 + 2x4 ≤ 4; x2 + x5 ≤ 4; x3 + x6 ≤ 6; x1 ≤ 3; x2 ≤ 4; x3 ≤ 4; x4 ≤ 1; x5 ≤ 2; x6 ≤ 6; x1, x2, x3, x4, x5, x6 ≥ 0; x4, x5, x6 là các biến nguyên. Your choice: 1 49 Sau đó lần lượt xuất hiện các dòng nhắc nhập dữ liệu đầu vào, ta lần lượt nhập như sau: +++++ THE PROBLEM +++++ Number of variables (N):6 Array to describe integer variables:[0 0 0 1 1 1] Goal function f(x)=(X(1)^0.6)+(X(2)^0.6)+(X(3)^0.4)+2*X(4)+5*X(5)-4*X(3)- X(6) Feasible conditions: (Leave blank for none) 1. X(1)=(X(2)-3*X(4))/3 2. X(3)=2*(X(2)+X(5)) 3. X(6)=4*X(4) 4. X(2)-3*X(1)-3*X(4)==0 5. X(3)-2*X(2)-2*X(5)==0 6. 4*X(4)-X(6)==0 7. X(1)+2*X(4)<=4 8. X(2)+X(5)<=4 9. X(3)+X(6)<=6 10. Number of feasible conditions entered: 9 Array to describe min values of X: [0 0 0 0 0 0] Array to describe max values of X: [3 4 4 1 2 6] +++++ CONDITIONS FOR ALGORITHM +++++ Limit Random Trying to find one feasible solution to (1000): 1000 Number of solutions in array A (20): 30 Limit Number of Interactive to(1000): 1000 Limit Number of Consecutive searches to (500): 200 Limit Number of Improve Failures to (500): 200 Epsilon 1 (1e-6): 1e-006 Epsilon 2(1e-1): 0.1 Beta (1): 1 Enter file name to save this problem (*.mat, leave blank for unsave):demo_01.mat Kết quả nhận được như sau: ===== Here is the result ===== Stop with number of failure on improve goal value exceed. X* = 0.6667 2.0000 4.0000 0 0 0 f(X*) = -11.9591 ± 1.1959e-005 50 Total time: 6.875 second(s) Source: demo_01.mat Results saved to res.txt End. 3.2. Nhập từ tệp Bài toán xác định cơ cấu đầu tư 1. Enter the problem. 2. Load problem from file. 3. Exit Your choice:2 Enter file name (*.mat):p01_50.mat The problem is: MIN f(X)=- 19.3753*X(1)^0.236*X(2)^0.104*X(3)^0.096*X(4)^0.056*X(5)^0.056*exp(0.168*X(6 ))*exp(0.066*X(7)) 1. X(5)=50-(X(1)+X(2)+X(3)+X(4)); 2. X(6)+X(7)==1 3. X(1)+X(2)+X(3)+X(4)+X(5)==50 0 <= X(1) <= 50 real 0 <= X(2) <= 50 real 0 <= X(3) <= 50 real 0 <= X(4) <= 50 real 0 <= X(5) <= 50 real 0 <= X(6) <= 1 integer 0 <= X(7) <= 1 integer Epsilon approximately: 0.0001% Searching. Be patient... Searching completed. ===== Here is the result ===== Stop with epsilon approximately. X* = Columns 1 through 6 21.0902 9.4856 9.2223 5.0966 5.1052 1.0000 Column 7 0 51 f(X*)= -88.3465 ± 8.8346e-005 Total time: 5.217 second(s) Source: p01_50.mat Results saved to res.txt End. Chú ý: – Để xoá màn hình dùng lệnh: clc – Để đọc tệp *.mat dùng lệnh: open(‘*.mat) 52 Chương IV GIẢI BÀI TOÁN QUY HOẠCH TUYẾN TÍNH ĐA MỤC TIÊU BẰNG PHƯƠNG PHÁP THOẢ DỤNG MỜ 1. CÁC KHÁI NIỆM CƠ BẢN 1.1. Phát biểu mô hình Trong các bài toán kĩ thuật, công nghệ, quản lý, kinh tế nông nghiệp v.v... nảy sinh từ thực tế, chúng ta thường phải xem xét để tối ưu hoá đồng thời một lúc nhiều mục tiêu. Các mục tiêu này thường là khác nhau về thứ nguyên, tức là chúng được đo bởi các đơn vị khác nhau. Những tình huống như vậy tạo ra các bài toán tối ưu đa mục tiêu. Như vậy, chúng ta cần phải tối ưu hoá (cực đại hoá hoặc cực tiểu hoá tuỳ theo tình huống thực tế) không phải là chỉ một mục tiêu nào đó, mà là đồng thời tất cả các mục tiêu đã đặt ra. Bài toán tối ưu đa mục tiêu mà trong đó miền ràng buộc D là tập lồi đa diện và các mục tiêu zi = zi(x), với i = 1, 2,, p, là các hàm tuyến tính xác định trên D, được gọi là BTQHTT đa mục tiêu. Với mục đích tìm hiểu bước đầu, BTQHTT đa mục tiêu (BTQHTT đa mục tiêu) được phát biểu như sau (Bài toán 1): Max Cx với ràng buộc x ∈ D, trong đó: C là ma trận cấp p × n và D = {x ∈ Rn: Ax ≤ b} với A là ma trận cấp m × n và b ∈ Rm. Các hàng của ma trận C là các véc tơ gradient c1, c2, , cp của các hàm mục tiêu z1 = c1Tx , z2 = c2Tx , , zp = cpTx. Ví dụ: Giải BTQHTT hai mục tiêu. z1 = 8x1+ 6x2 → Max z2 = x1 + 3x2 → Max với các ràng buộc: Ta có thể viết bài toán này dưới dạng ma trận như sau: Max Cx với ràng buộc x ∈ D = {x∈ R2 : Ax ≤ b}, trong đó x = (x1, x2)T, b = (60, 48)T, còn C = 8 1 ⎡⎢⎣ 6 3 ⎤⎥⎦ , A = 4 2 2 4 ⎡ ⎤⎢ ⎥⎣ ⎦ . Có thể nói, BTQHTT đa mục tiêu là BTQHTT mà trong đó chúng ta phải tối ưu hoá cùng một lúc nhiều mục tiêu. Tuy nhiên, các mục tiêu này thường đối chọi cạnh tranh với nhau. Việc làm tốt hơn mục tiêu này thường dẫn tới việc làm xấu đi 4x1 + 2x2 ≤ 60 2x1 + 4x2 ≤ 48 x1, x2 ≥ 0. (D) 53 một số mục tiêu khác. Vì vậy việc giải các bài toán tối ưu đa mục tiêu, tức là tìm ra một phương án khả thi tốt nhất theo một nghĩa nào đó, thực chất chính là một bài toán ra quyết định. 1.2. Phương án tối ưu Pareto Khái niệm then chốt trong tối ưu hoá đa mục tiêu là khái niệm phương án tối ưu Pareto. Xét Bài toán 1 , chúng ta cần biết các định nghĩa và định lý sau đây. Định nghĩa 1 Một phương án tối ưu Pareto x* có tính chất sau đây: − Trước hết nó phải thuộc vào miền các phương án khả thi của bài toán, tức là phải thoả mãn tất cả các ràng buộc: x* ∈ D. − Xét phương án khả thi x ∈ D, x ≠ x*. Nếu tồn tại tại một chỉ số i ∈ {1, 2, , p} sao cho zi(x) > zi(x*) thì tồn tại j ∈ {1, 2, , p}, j ≠ i, sao cho zj(x) < zj(x*). Nói một cách khác, không tồn tại một phương án khả thi nào x ∈ D có thể trội hơn x* trên tổng thể tất cả các mục tiêu. Định nghĩa 2 Một phương án tối ưu Pareto yếu x* có tính chất sau đây: − Trước hết nó phải thuộc vào miền các phương án khả thi của bài toán, tức là phải thoả mãn tất cả các ràng buộc: x* ∈ D. − Xét một phương án khả thi x ∈ D, x ≠ x*. Nếu tồn tại tại một chỉ số i ∈ {1, 2, , p} sao cho zi(x) > zi(x*) thì tồn tại j ∈ {1, 2, , p}, j ≠ i, sao cho zj(x) ≤ zj(x*). Để nhận biết tập phương án tối ưu Pareto chúng ta cần tới các định nghĩa sau. Định nghĩa 3 Nón cảm sinh bởi các véc tơ gradient c1, c2, , cp của các hàm mục tiêu được gọi là nón tiêu chuẩn (criterion cone). Để tìm tập các phương án tối ưu Pareto chúng ta có thể sử dụng tập các điểm trội. Định nghĩa 4 Cho x ∈ D. Tập điểm trội tại x là tập xD = { x }⊕ C≥ , với C≥ = {x = (x1, x2) ∈ R2 : Cx ≥ 0, Cx ≠ 0}là nón đối cực nửa dương (semi-positive polar cone). Định lý 1: Xét Bài toán 1. Lúc đó: x ∈ D là phương án tối ưu Pareto khi và chỉ khi xD ∩ D = { x }. Chứng minh. Giả sử x là phương án tối ưu Pareto và xD ∩ D ≠ { x }. Lúc đó tồn tại xˆ ∈ xD ∩ D sao cho xˆ ≠ x và xˆ = x + x với x ∈ C≥ . Do Cx ≥ 0, Cx ≠ 0 nên C xˆ ≥ C x và C xˆ ≠ C x . Điều này vô lí do x là phương án tối ưu Pareto. Ngược lại, giả sử xD ∩ D = { x }. Lúc này nếu tồn tại xˆ ≠ x sao cho C xˆ ≥ 54 C x và C xˆ ≠ C x thì xˆ ∉ D. Vậy x là phương án tối ưu Pareto. „ Để minh hoạ các định nghĩa 1, 3 và 4, chúng ta xét lại ví dụ đã biết. Ví dụ: Giải BTQHTT hai mục tiêu. z1 = 8x1+ 6x2 → Max z2 = x1 + 3x2 → Max với các ràng buộc: Miền các phương án khả thi D (miền giới hạn bởi tứ giác ABCD) được biểu thị trên hình I,