Phương pháp nhánh cận giải bài toán quy hoạch nguyên

MỤC LỤC • MỤC LỤC............................................................................................................................i • Lời cảm ơn...........................................................................................................................ii • Lời nói đầu .........................................................................................................................iii 1 QUY HOẠCH NGUYÊN 1.1 Giới thiệu........................................................................................................................1 1.2 Quy hoạch nguyên trong thực tế......................................................................................6 2 PHƯƠNG PHÁP NHÁNH CẬN 2.1 Giới thiệu thuật toán trên ví dụ.....................................................................................30 2.2 Chia nhánh, đặt cận, chạm đáy......................................................................................34 2.3 Thuật toán nhánh - và - cận...........................................................................................36 3 GIẢI BÀI TOÁN NGUYÊN TRÊN MATLAB 3.1 Giới thiệu Matlab..........................................................................................................39 3.2 Lập trình thuật toán trên bài toán nguyên.....................................................................40 3.3 Giải bài toán nguyên trên Matlab..................................................................................44 TÀI LIỆU THAMKHẢO...............................................................................................49 i Lời cảm ơn Tôi chân thành cám ơn công lao to lớn của PGS. TS Trần Huệ nương - Giảng viên hướng dẫn của tôi, người giới thiệu tôi đến với đề tài này và giúp tôi đạt được những kết quả tốt. Cô đã hỗ trợ tôi trong cả khi thành công và thất bại, động viên tôi tập trung bất cứ khi nào tôi nản chí; nếu không có những sự giúp đỡ ấy của Cô thì tôi nghĩ luận văn này không bao giờ hoàn thành. Viết luận văn Thạc sỹ là một quá trình vất vả, mất nhiều tháng liên tục. Tuy nhiên, tôi đã nhận được sự giúp đỡ, khuyến khích của những bạn bè, bạn học, Giáo sư. Với tất cả niềm vui tôi xin cám ơn tất cả. Điều quan trọng nhất, tôi không thể làm được luận văn này nếu không có sự giúp đỡ của gia đình. Luận văn không thể hoàn thành mà không có tình yêu, sự nhẫn nại, sự hỗ trợ giúp đỡ. Tôi chân thành cám ơn. TP. Hồ Chí Minh, ngày 01 tháng 05 năm 2009. Nguyễn Hữu Thương ii Lời nói đầu Chuyên ngành quy hoạch tuyến tính được các nhà Toán học trên thế giới nghiên cứu và phát triển kể từ sau đại chiến thế giới lần thứ hai. Người tiên phong trong lĩnh vực này là G.B.Dantzig. Có nhiều bài toán thực tế thuộc các lĩnh vực khác nhau có thể mô tả toán học bằng quy hoạch tuyến tính. Quy hoạch nguyên (IP - Integer Programming) rất phổ biến trong thực tế. Mảng bài toán này có vẻ đơn giản nhất mà cũng quan trọng nhất trong các bài toán chọn quyết định (có liên quan nhau nhưng đều ở dạng có hoặc không). Một số ý tưởng để giải IP là phương pháp mặt phẳng cắt của Gomory công bố năm 1958, phương pháp nhánh - và - cận (branch - and - bound method) xuất hiện trong [Land - Doig 1960]. Tác giả đã lựa chọn phương pháp nhánh - và - cận với những ưu điểm rõ rệt để trình bày trong luận văn này. Luận văn gồm ba chương: Quy hoạch nguyên, phương pháp nhánh - và cận và Một số ứng dụng của phương pháp nhánh cận giải quy hoạch nguyên. Tuy nhiên, luận văn cũng không tránh khỏi thiếu sót. Rất mong nhận được sự đóng góp của Quý Thầy Cô và các bạn đồng nghiệp. iii Chương 1 Quy hoạch nguyên 1.1 Giới thiệu Quy hoạch nguyên (Integer programming - IP) là bài toán quy hoạch trong đó tất cả hoặc một phần các biến bị ràng buộc chỉ lấy giá trị nguyên. Tương ứng tên cụ thể cho từng trường hợp là quy hoạch nguyên hoàn toàn (pure integer programming) và quy hoạch nguyên bộ phận (mixed integer programming - MIP). Tuy vậy, thuật ngữ "Quy hoạch nguyên" thực tế chỉ dành cho lĩnh vực toán nghiên cứu lớp bài toán quy hoạch tuyến tính nguyên. Vì vậy hai thuật ngữ được hiểu như nhau. Lớp bài toán này rất phổ biến trong thực tế. Mảng bài toán này có vẻ đơn giản nhất mà cũng quan trọng nhất trong lớp này là các bài toán chọn quyết định (có liên quan nhau nhưng đều ở dạng "có" hoặc "không" (yes - or - no decision). Chẳng hạn bài toán bổ nhiệm: Giả sử có tập S gồm m người và tập D gồm m việc hoặc chức vụ. Cước phí của việc bổ nhiệm người i ∈ S vào việc j ∈ D là cij. Bài toán đặt ra là tìm cách chia cho mỗi người đúng một việc sao cho tổng cước phí là nhỏ nhất. Đặt biến là xij = { 1 nếu người i được việc j, 0 trường hợp khác thì bài toán trở thành min ∑ i∈S ∑ j∈D cijxij, 1 ∑ i∈S xij = 1, i ∈ S (mỗi người chỉ được đúng một việc), ∑ j∈D xij = 1, j ∈ D (mỗi việc chỉ dành cho đúng một người). Nhiều bài toán khác quyết định có - không như vậy, chẳng hạn đầu tư không? có mở đại lý ở thành phố A không? có nhận thực hiện dự án không?.... Vì các biến (quyết định) chỉ nhận hai giá trị (biến 0-1) nên bài toán có tên là quy hoạch nguyên nhị phân (binary integer programming - BIP). Nếu bài toán có n biến và các ràng buộc cho thấy miền chấp nhận được là giới nội trong Rn, thì bài toán "chỉ có" hữu hạn nghiệm chấp nhận được, không như quy hoạch tuyến tính (tổng quát) có miền này là cả một tập lồi đa diện. Tuy vậy, sẽ rất sai lầm nếu nghĩ rằng quy hoạch nguyên dễ giải hơn, hoặc hơn nữa nghĩ rằng có thể điểm diện (enumeration) các phương án (tức là nghiệm) chấp nhận được để chọn phương án tối ưu. Chẳng hạn BIP với n biến có 2n phương án chấp nhận được. Nhưng đây là số khổng lồ nếu n lớn. Trong lịch sử Trung Hoa cổ đại đã có minh hoạ rất hay là nhà vua không đủ thóc, dù vét hết kho trong nước, để thưởng 264 hạt thóc theo ước muốn (được vua ban) của người nghĩ ra cờ tướng làm vua thích thú (bàn cờ đó có 64 ô). Nguyên nhân vì miền chấp nhận được, tuy hữu hạn, của quy hoạch nguyên không có các đỉnh để thực hiện phương pháp đơn hình, đã biết là rất hiệu quả. Chính vì vậy quy hoạch tuyến tính dễ giải hơn. Người ta có thể giải mọi quy hoạch tuyến tính với nhiều ngàn biến và ràng buộc bằng máy tính trong vài giờ, nhưng đến nay vẫn không đảm bảo giải được với thời gian tương tự bất kỳ dù chỉ là BIP với 100 biến. Người ta thấy rằng ở quy hoạch tuyến tính thì độ phức tạp tính toán tăng rõ theo số ràng buộc nhưng với quy hoạch nguyên thì số biến ảnh hương quyết định. Một số ý tưởng tự nhiên để giải IP (quy hoạch nguyên) là giải quy hoạch tuyến tính tương đương, tức là IP nhưng bỏ đi ràng buộc giá trị biến phải nguyên (quy hoạch tuyến tính này gọi là LP - nới lỏng (LP relaxation) của IP) rồi làm tròn nghiệm tối ưu thành nghiệm nguyên. Nó sẽ xấp xỉ nghiệm tối ưu của IP. Cách này có thể áp dụng trong thực tế, nhưng phải cẩn thận hai nguy cơ sau đây. Một là nghiệm làm tròn có thể thậm chí 2 không chấp nhận được đối với IP như ví dụ sau Ví dụ 1.1.1 Xét quy hoạch nguyên max z = x2, −x1 + x2 6 1 2 , x1 + x2 6 3 1 2 , x1 > 0, x2 > 0, x1, x2 nguyên. Hình 1.1: Nghiệm tối ưu làm tròn không chấp nhận được 3 H.1.1 cho thấy điểm A(3/2, 2) nới lỏng, B(1, 2) và C(2, 2) là hai nghiệm làm tròn từ A nhưng không chấp nhận được. Nghiệm tối ưu của IP là D(1, 1) và E(2, 1). Nguy cơ thứ hai là nghiệm làm tròn có thể chấp nhận được nhưng có thể có giá trị mục tiêu rất xa mục tiêu tối ưu của IP như ví dụ sau đây cho thấy. Ví dụ 1.1.2 Xét quy hoạch nguyên max z = x1 + 5x2, x1 + 10x2 6 20, x1 6 2, x1 > 0, x2 > 0, x1, x2nguyên. H.1.2 cho thấy điểm A(2, 7/4) là nghiệm tối ưu của LP - nới lỏng (với mục tiêu z = 11). B(2, 1) là nghiệm làm tròn cho IP, chấp nhận được và ứng mục tiêu z = 7, xa với mục tiêu tối ưu của IP là z = 10 đạt tại điểm C(0, 2). Bài tập: "Chứng minh rằng nếu quy hoạch nguyên tuyến tính có ít nhất một ràng buộc đẳng thức thì nghiệm làm tròn từ LP nới lỏng không thể là nghiệm chấp nhận được của quy hoạch nguyên" cho thấy trường hợp nghiệm làm tròn không chấp nhận được là rất phổ biến. Có thể nghĩ rằng dữ liệu của bài toán đã là xấp xỉ nên tính không chấp nhận được này có thể "tha thứ". Nhưng nhiều bài toán IP thực tế có dữ liệu chính xác. Chẳng hạn ràng buộc kiểu lựa chọn đúng một trong n quyết định "có" còn lại là "không" 4 Hình 1.2: Nghiệm tối ưu làm tròn xa nghiệm tối ưu của IP là ràng buộc với hệ số chính xác x1 + x2 + . . .+ xn = 1, xj = 0 hoặc 1. Hơn nữa ở bài toán nhị phân với biến có ý nghĩa là "có" hay "không" thì giá trị không nguyên là vô nghĩa rồi, nên không thể làm tròn đi nữa. Vì quy hoạch nguyên hay gặp như thế, lại không thể nhận được từ làm tròn nghiệm nhận được từ LP - nới lỏng của nó, nhu cầu tìm các thuật toán riêng là rất lớn, vẫn đang rất được quan tâm nghiên cứu hiện nay. Do giải quy hoạch nguyên khó và tính toán phức tạp nên với các bài toán cỡ lớn cách tiếp cận tốt nhất là dùng các thuật toán phán đoán (heuristic algorithm), tuy không đảm bảo đi đến nghiệm tối ưu trong mọi trường hợp, nhưng thường rất hiệu quả. Với các bài toán không lớn quá, đã có nhiều thuật toán giải chính xác. Cách tiếp cận chung là kết hợp việc giải các LP - nới lỏng với điểm diện ẩn (implicit emmeration), tức là điểm diện các nghiệm chấp nhận được nhưng ở mỗi bước tìm cách loại bỏ số đáng kể các nghiệm không thể là tối ưu không cần điểm diện. Trước khi trình bày phương pháp chính là phương pháp nhánh - và - cận ở mục sau, Mục ... giới 5 thiệu .... Để kết thúc mục này, ta đã thấy nhiều trường hợp may mắn là giải quy hoạch tuyến tính (không nguyên) bằng thuật toán đơn hình ta nhận được nghiệm tối ưu là nguyên. Đó là bài toán dòng trên mạng với các dữ liệu nguyên, mà trường hợp riêng là các bài toán sau với dữ liệu nguyên: bài toán dòng cực đại, bài toán đường ngắn nhất, bài toán vận tải. Nếu khi giải LP - nới lỏng của quy hoạch nguyên mà gặp may mắn như vậy thì ta được luôn nghiệm tối ưu của quy hoạch nguyên. 1.2 Quy hoạch nguyên trong thực tế Ý tưởng giải một bài toán nguyên trong thực tế: Bước 1. Giải bài toán tuyến tính gốc sau khi bỏ đi điều kiện biến nguyên, Bước 2. Giả sử được nghiệm tối ưu xk = t. Nếu nghiệm tối ưu đạt được thoả ràng buộc nguyên thì kết thúc bước lặp. Ngược lại đến bước 3, Bước 3. Chia bài toán thành hai phần: LP1: xk ≤ [t] max z = cx, Ax ≤ b, xk ≤ [t], x ≥ 0. LP2: xk ≥ [t] + 1 6 max z = cx, Ax ≤ b, xk ≥ [t] + 1, x ≥ 0. Bước 4. Giải LP1 và LP2 một cách riêng rẽ, Bước 5. Nếu với bất kỳ bài toán con, nghiệm tối ưu nguyên đạt được thì bài toán không chia nhánh nữa. Mặt khác quay trở lại Bước 3. Ví dụ 1.2.1 (Bài toán Balo - knapsack problem) Một nhà thám hiểm chỉ được mang theo một ba lô trọng lương không quá b, Có n loại vật dụng nên mang theo. Loại vật j có trọng lượng mỗi vật là aj và giá trị sử dụng mỗi vật là cj. Hỏi ông phải mang thế nào để có giá trị sử dụng lớn nhất? Gọi xj là số lượng vật j mang theo thì mô hình toán học của bài toán balo này là quy hoạch nguyên sau: max z = n∑ j=1 cjxj, n∑ j=1 ajxj ≤ b, xj ≥ 0 và nguyên, j = 1, . . . , n. 7 Ta xét và giải bài toán balo hai biến như sau: Một học sinh mang theo balo trọng lượng không quá 12kg, mang theo hai vật, trọng lượng vật một là 3kg, trọng lượng vật hai tương ứng là 2kg, giá trị sử dụng của cả hai vật là 1. Gọi x1, x2 là số lượng vật một và hai mang theo, sao cho vật hai mang không quá hai cái. Hỏi phải mang theo bao nhiêu để giá trị sử dụng lớn nhất? Mô hình bài toán: max z = x1 + x2, 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1, x2 ≥ 0 và nguyên. Đầu tiên, ta giải bài toán trên bằng phương pháp đơn giản (bỏ đi điều kiện nguyên), Sau khi thêm biến bù, ta có: 3x1 + 2x2 + w3 = 12, x2 + w4 = 2, với w3, w4 là những biến bù. Nghiệm cơ bản chấp nhận được là: x1 = x2 = 0, z = 0, Giải bài toán đơn hình dạng bảng: 8 Ta có nghiệm tối ưu: x1 = 8/3, x2 = 2 Vì nghiệm tối ưu đạt được không phải nghiệm nguyên, ta chọn x1 = 8/3(2 + 2/3) mà có giá trị phân số và chia bài toán thành hai bài toán con (LP2 và LP3) bằng cách thêm hai ràng buộc mới x1 ≤ 2, x1 ≥ 3, bởi vì [x1] = [8/3] = 2. LP2: max z = x1 + x2, 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1 ≤ 2, x1, x2 ≥ 0. Thêm biến bù: 3x1 + 2x2 + w5 = 12, x2 + w6 = 2, x1 + w7 = 2, 9 với w5, w6, w7 là những biến bù. Giải bài toán đơn hình dạng bảng: Nghiệm tối ưu của LP2: x1 = x2 = 2. Vì tất cả các biến trong LP2 có giá trị nguyên, nên không rẽ nhánh LP2 thêm nữa. LP3: max z = x1 + x2, 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1 ≥ 3, x1, x2 ≥ 0. Thêm biến bù: 3x1 + 2x2 + w8 = 12, x2 + w9 = 2, 10 x1 − w10 + A1 = 3, với w8, w9 là những biến bù, w10 là biến phụ, A1 là biến giả. Giải bài toán đơn hình dạng bảng: Nghiệm tối ưu của LP3: x1 = 3, x2 = 3/2. Nghiệm tối ưu đạt được không nguyên nên không là nghiệm tối ưu của bài toán. Vì x2 = [1 + 1/2] = 1, ta có thể chia nhánh LP3 thành hai bài toán con (LP4 và LP5) bằng cách thêm các ràng buộc x2 ≤ 1, x2 ≥ 2. LP4: max z = x1 + x2, 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1 ≥ 3, x2 ≤ 1, x1, x2 ≥ 0. 11 Thêm biến bù: 3x1 + 2x2 + w11 = 12, x2 + w12 = 2, x1 − w13 + A1 = 3, x2 + w14 = 1, với w11, w12, w14 là những biến bù, w13 là biến phụ, A1 là biến giả. Giải bài toán đơn hình dạng bảng: Nghiệm tối ưu của LP4: x1 = 10/3, x2 = 1. Nghiệm tối ưu đạt được không nguyên nên không là nghiệm tối ưu của bài toán. Vì x1 = [3 + 1/3] = 3, ta có thể chia nhánh LP4 thành hai bài toán con (LP6 và LP7) bằng cách thêm các ràng buộc x1 ≤ 3, x1 ≥ 4. LP5: max z = x1 + x2, 12 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1 ≥ 3, x2 ≥ 2, x1, x2 ≥ 0. Nghiệm của LP5 không chấp nhận được nên không xét LP5 nữa. LP6: max z = x1 + x2, 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1 ≥ 3, x2 ≤ 1, x1 ≤ 3, x1, x2 ≥ 0. 13 Thêm biến bù, biến phụ, biến giả: 3x1 + 2x2 + w15 = 12, x2 + w16 = 2, x1 − w17 + A1 = 3, x2 + w18 = 1, x1 + w19 = 3, với w15, w16, w18, w19 là những biến bù, w17 là biến phụ, A1 là biến giả. Giải bài toán đơn hình dạng bảng, ta có: LP6 có nhiều nghiệm tối ưu, khi x2 = 0 thì trùng nghiệm cơ bản. Nghiệm tối ưu của LP6: x1 = 3, x2 = 1. Vì tất cả các biến trong LP6 có giá trị nguyên, nên không rẽ nhánh LP6 thêm nữa. LP7: max z = x1 + x2, 3x1 + 2x2 ≤ 12, 14 x2 ≤ 2, x1 ≥ 3, x2 ≤ 1, x1 ≥ 4, x1, x2 ≥ 0. Thêm biến bù: 3x1 + 2x2 + w20 = 12, x2 + w21 = 2, x1 − w22 + A1 = 3, x2 + w23 = 1, x1 − w24 + A2 = 4, với w20, w21, w23 là những biến bù, w22, w24 là biến phụ, A1, A2 là biến giả. Giải bài toán đơn hình dạng bảng, ta có: 15 Nghiệm tối ưu của LP7: x1 = 4, x2 = 0, so sánh với giá trị mục tiêu của bài toán gốc và điều kiện ràng buộc biến nguyên thì nghiệm của LP7 cũng là nghiệm tối ưu của bài toán. Vì tất cả các biến trong LP7 có giá trị nguyên, nên không rẽ nhánh LP7 thêm nữa. Bài toán đạt nghiệm nguyên chấp nhận được tại LP2, LP6, LP7 đều có giá trị mục tiêu là z = 4 nhưng tổng trọng lượng LP2 mang được là 10, LP6 là 11 và LP7 là 12, vậy nghiệm tối ưu của cả bài toán là LP7 với tối ưu giá trị sử dụng z = 4 và tối đa trọng lượng mang theo. Tập nghiệm tối ưu cho bài toán gốc là: x∗ = (4, 0), z∗ = 4 Sơ đồ cây cho quá trình giải bài toán bằng phương pháp nhánh - và - cận có dạng như sau: 16 Ví dụ 1.2.2 (Capital Budgeting problem) Giả sử chúng ta ước đầu tư một giá trị tối đa b và có n mục tiêu đầu tư. Đầu tư mục tiêu j yêu cầu aj sẽ có giá trị "present" (giá trị chiết khấu kỳ hạn) là cj. Chúng ta nên đầu tư như thế nào để tổng giá trị "present" đạt maximum. Đây là dạng bài toán đầu tư "có - không" nên ta sẽ đưa về quy hoạch nguyên nhị phân: max z = n∑ j=1 cjxj, n∑ j=1 ajxj ≤ b, xj = { 0 không đầu tư 1 đầu tư xj là nguyên,j = 1, . . . , n. 17 Ta xét bài toán cụ thể với khoản đầu tư mơ ước tối đa b = 14, 000 USD, có bốn mục tiêu đầu tư, mục tiêu đầu tư một a1 = 5, 000 USD có giá trị c1 = 8, 000 USD, mục tiêu đầu tư hai a2 = 7, 000 USD có giá trị c2 = 11, 000 USD, mục tiêu đầu tư ba a3 = 4, 000 USD có giá trị c1 = 6, 000 USD, mục tiêu đầu tư bốn a4 = 3, 000 USD có giá trị c1 = 4, 000 USD. Hỏi phải đầu tư như thế nào để tổng giá trị "present" đạt max? Mô hình bài toán: max z = 8x1 + 11x2 + 6x3 + 4x4, 5x1 + 7x2 + 4x3 + 3x4 ≤ 14, xj ∈ {0, 1}, j = 1, . . . , 4. Theo đó, ta giải bài toán trên bằng phương pháp đơn giản (bỏ đi điều kiện nguyên), bài toán có dạng: Sau khi thêm biến, ta có: max z = 8x1 + 11x2 + 6x3 + 4x4, 5x1 + 7x2 + 4x3 + 3x4 + w5 = 14, x1 + w6 = 1, x2 + w7 = 1, x3 + w8 = 1, 18 x4 + w9 = 1, với w5, w6, w7, w8, w9 la biến bù. Phương pháp đơn hình dạng bảng: Chuyển hàm mục tiêu từ max thành min, ta giải: max z = −min−z = −8x1 − 11x2 − 6x3 − 4x4,  x1 x2 x3 x4 w5 w6 w7 w8 w9 −8 −11 −6 −4 0 0 0 0 0 0 5 7 4 3 1 0 0 0 0 14 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 1  Chọn cột xoay là cột 2, hàng xoay là hàng 3, ta có: x1 x2 x3 x4 w5 w6 w7 w8 w9 −8 0 −6 −4 0 0 11 0 0 11 5 0 4 3 1 0 −7 0 0 7 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 1  Chọn cột xoay là cột 1, hàng xoay là hàng 2, ta có: 19  x1 x2 x3 x4 w5 w6 w7 w8 w9 0 0 −6 −4 0 8 11 0 0 19 0 0 4 3 1 −5 −7 0 0 2 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 1 1  Chọn cột xoay là cột 3, hàng xoay là hàng 1, ta có: x1 x2 x3 x4 w5 w6 w7 w8 w9 0 0 0 1/2 3/2 1/2 1/2 0 0 22 0 0 1 3/4 1/4 −5/4 −7/4 0 0 1/2 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 1 0 0 0 −3/4 −1/4 5/4 7/4 1 0 1/2 0 0 0 1 0 0 0 0 1 1  Hệ số trên hàng đầu không âm, nghiệm tối ưu là : x1 = x2 = 1, x3 = 0.5, x4 = 0, z = 22 Chúng ta thấy nghiệm không nguyên, nên chúng ta phải đảm bảo x3 nguyên và rẽ nhánh x3 = 0 hoặc x3 = 1 và chia thành hai bài toán con LP2 và LP3, LP2: x3 = 0 Sau khi thêm biến, ta có: max z = 8x1 + 11x2 + 4x4, 20 5x1 + 7x2 + 3x4 + w10 = 14, x1 + w11 = 1, x2 + w12 = 1, x4 + w13 = 1, với w10, w11, w12, w13 la biến bù. Chuyển hàm mục tiêu từ max thành min, ta giải: max z = −min−z = −8x1 − 11x2 − 4x4,  x1 x2 x4 w10 w11 w12 w13 −8 −11 −4 0 0 0 0 0 5 7 3 1 0 0 0 14 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 1  Chọn cột xoay là cột 2, hàng xoay là hàng 3, ta có: x1 x2 x4 w10 w11 w12 w13 −8 0 −4 0 0 11 0 11 5 0 3 1 0 −7 0 7 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 1  21 Chọn cột xoay là cột 1, hàng xoay là hàng 2, ta có: x1 x2 x4 w10 w11 w12 w13 0 0 −4 0 8 11 0 19 0 0 3 1 −5 −7 0 2 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 1  Chọn cột xoay là cột 3, hàng xoay là hàng 1, ta có: x1 x2 x4 w10 w11 w12 w13 0 0 0 4/3 4/3 5/3 0 21.667 0 0 1 1/3 −5/3 −7/3 0 2/3 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 0 −1/3 5/3 7/3 1 1/3  Hệ số trên hàng đầu không âm, nghiệm tối ưu là : x1 = x2 = 1, x3 = 0, x4 = 0.667, z = 21.667 Nghiệm đạt được không nguyên nên không là nghiệm tối ưu của bài toán ban đầu. LP3: x3 = 1 Sau khi thêm biến, ta có: max z = 8x1 + 11x2 + 6 + 4x4, 5x1 + 7x2 + 3x4 + w14 = 10, 22 x1 + w15 = 1, x2 + w16 = 1, x4 + w17 = 1, với w14, w15, w16, w17 la biến bù. Chuyển hàm mục tiêu từ max thành min, ta giải: max z = −min−z = −8x1 − 11x2 − 6− 4x4,  x1 x2 x4 w14 w15 w16 w17 −8 −11 −4 0 0 0 0 6 5 7 3 1 0 0 0 10 1 0 0 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 1 0 0 0 1 1  Chọn cột xoay là cột 2, hàng xoay là hàng 1 (do nếu chọn hàng xoay là 2 bài toán trở lại ban đầu), ta có: x1 x2 x4 w14 w15 w16 w17 −1/7 0 5/7 11/7 0 0 0 152/7 5/7 1 3/7 1/7 0 0 0 10/7 1 0 0 0 1 0 0 1 −5/7 0 −3/7 −1/7 0 1 0 −3/7 0 0 1 0 0 0 1 1  23 Chọn cột xoay là cột 1, hàng xoay là hàng 2, ta có: x1 x2 x4 w14 w15 w16 w17 0 0 5/7 11/7 1 0 0 153/7 0 1 3/7 1/7 1 0 0 5/7 1 0 0 0 1 0 0 1 0 0 2/7 4/7 5/7 1 0 2/7 0 0 1 0 0 0 1 1  Hệ số trên hàng đầu không âm, nghiệm tối ưu là : x1 = x3 = 1, x4 = 0, x2 = 0.714, z = 21.85 Nghiệm đạt được không nguyên nên không là nghiệm tối ưu của bài toán ban đầu. Giải LP3 và LP2 vẫn chưa đạt nghiệm nguyên, vì vậy chúng ta sẽ chọn bài toán con và rẽ nhánh 1 trong các biến theo tiêu chuẩn sau: • Chọn bài toán con mà chưa từng chọn trước đây, • Chọn bài toán con với giá trị mục tiêu cao nhất. Trong trường hợp này, ta rẽ nhánh bài toán LP3 với x2 = 0 hoặc x2 = 1 thành hai bài toán con LP4 và LP5. LP4: x3 = 1, x2 = 0 Sau khi thêm biến, ta có: max z = 8x1 + 6 + 4x4, 5x1 + 3x4 + w18 = 10, 24 x1 + w19 = 1, x4 + w20 = 1, với w18, w19, w20 la biến bù. Chuyển hàm mục tiêu từ max thành min, ta giải: max z = −min−z = −8x1 − 6− 4x4,  x1 x4 x18 x19 x20 −8 −4 1 0 0 6 5 3 1 0 0 10 1 0 0 1 0 1 0 1 0 0 1 1  Chọn cột xoay là cột 1, hàng xoay là hàng 2, ta có: x1 x4 w18 w19 w20 0 −4 0 8 0 14 0 3 1 −5 0 5 1 0 0 1 0 1 0 1 0 0 1 1  Chọn cột xoay là cột 2, hàng xoay là hàng 3, ta có: 25  x1 x4 w18 w19 w20 0 0 0 8 4 18 0 0 1 −5 3 2 1 0 0 1 0 1 0 1 0 0 1 1  Hệ số trên hàng đầu không âm, nghiệm tối ưu là : x1 = x3 = x4 = 1, x2 = 0, z = 18 Chúng ta bây giờ đã có một nghiệm chấp nhận được với giá trị mục tiêu là 18. LP5: x3 = 1, x2 = 1 Sau khi thêm biến, ta có: max z = 8x1 + 17 + 4x4, 5x1 + 3x4 + w21 = 3, x1 + w22 = 1, x4 + w23 = 1, với w21, w22, w23 la biến bù. Chuyển hàm mục tiêu từ max thành min, ta giải: max z = −min−z = −8x1 − 17− 4x4, 26  x1 x4 w21 w22 w23 −8 −4 0 0 0 17 5 3 1 0 0 3 1 0 0 1 0 1 0 1 0 0 1 1  Chọn cột xoay là cột 1, hàng xoay là hàng 1 (do nếu chọn hàng xoay là 2 bài toán trở lại ban đầu), ta có:  x1 x4 w21 w22 w23 0 4/5 13/5 0 0 21.8 1 3/5 1/5 0 0 3/5 0 −3/5 −1/5 1 0 2/5 0 1 0 0 1 1  Hệ số trên hàng đầu không âm, nghiệm tối ưu là : x1 = 0.6, x3 = 1, x4 = 0, x2 = 1, z = 21 Nghiệm đạt được không nguyên nên không là nghiệm tối ưu của bài toán ban đầu. Chúng ta thấy nghiệm không nguyên, nên chúng ta phải đảm bảo x1 nguyên và rẽ nhánh x1 = 0 hoặc x1 = 1 và chia thành hai bài toán con LP6 và LP7, LP6: x1 = 0, x2 = x3 = 1 Nghiệm tối ưu: x1 = 0, x2 = x3 = x4 = 1, z = 21 LP7: x1 = 1, x2 = x3 = 1 27 Nghiệm tối ưu: không chấp nhận được Bài toán có hai nghiệm tối ưu thoả ràng buộc biến nguyên ở LP4 và LP6 nhưng so sánh với giá trị mục tiêu thì ta thấy nghiệm tối ưu của LP6 cũng là nghiệm tối ưu của bài toán do giá trị mục tiêu z = 21 là tốt nhất. Tập nghiệm tối ưu: x∗ = (0, 1, 1, 1)z∗ = 21 Sơ đồ cây cho quá trình giải bài toán bằng phương pháp nhánh - và - cận có dạng như sau: 28 29 Chương 2 Phương pháp nhánh - và - cận 2.1 Giới thiệu thuật toán trên ví dụ Phương pháp nhánh - và - cận (branch - and - bound method) xuất hiện đầu tiên trong [Land - Doig 1960] và nhất là dạng hoàn thiện trong [Dakin 1965], nó trở nên có ưu thế rõ rệt trong việc giải bài toán nguyên. Ta sẽ trình bày phương pháp thông qua ví dụụ Hình 2.1: Ví dụ 2.1.1 Xét quy hoạch nguyên (IP) 30 max z = 5x1 + 4x2, x1 + x2 ≤ 5, 10x1 + 6x2 ≤ 45, x1, x2 ∈ Z+. Gọi LP0 là quy hoạch tuyến tính tương ứng ở bước lặp 0, tức là thay x1, x2 ∈ Z − + bằng x1, x2 ≥ 0 có nghiệm tối ưu là x1 = 3.75, x2 = 1.25 và z = 23.75, xem H 2.1, ở đây A là điểm tối ưu của LP0 còn các dấu • là miền chấp nhận được (cnđ) của IP. Lấy bất kỳ một biến mà ở A có giá trị không nguyên, chẳng hạn x1 = 3.75. Vì phần miền chấp nhận được của LP0 thoả 3 < x1 < 4 không thể chứa nghiệm chấp nhận được của IP nên ta có thể bỏ qua không xét phần này và chia nhánh (branching) ra được hai miền và tương ứng là hai bài toán quy hoạch tuyến tính con ở bước 1 tiếp theo: 1. LP1 với miền chấp nhận được là miền chấp nhận được LP0, với x1 ≤ 3; 2. LP2 với miền chấp nhận được là miền chấp nhận được LP0, với x1 ≥ 4. x1 được gọi là biến chia nhánh (branching variable) H 2.2 (a) biểu diễn việc chia miền và (b) là sơ đồ cây của thuật toán ở bước đầu và đã được nghiệm tối ưu của IP Đến đây ta biết nghiệm tối ưu của IP, phải nằm trong LP (hoặc LP2 (tức là trong miền xác định của ít nhất một trong hay quy hoạch tuyến tính này). Ta đành giải bất kỳ, chẳng hạn LP1, là LP0 thêm ràng buộc x1 ≤ 3: max z = 5x1 + 4x2, 31 Hình 2.2: x1 + x2 ≤ 5, 10x1 + 6x2 ≤ 45, x1 ≤ 3, x1, x2 ∈ Z+. LP1 có thể giải bằng cách đưa về dạng chính tắc và dùng thuật toán đơn hình tiêu chuẩn. Cũng có thể giải bằng thuật toán đơn hình áp dụng cho quy hoạch tuyến tính tổng quát, có chứa cận dạng lj ≤ xj ≤ uj, xem [Khánh - Nương 2000], Chương 9. Nghiệm tối ưu của LP1 là x1 = 3, x2 = 2, z = 23, vô tình đã là nguyên, Đây là nghiệm tốt nhất tạm thời cho IP, ta gọi là một kỷ lục (incumbent). Do đó LP1 không thể suy ra nghiệm chấp nhận được tốt hơn cho IP, chỉ còn chờ so sánh với sự nghiêm cứu LP2. Ta nói đã làm chạm đáy (fathoming) LP1. Việc "chờ so sánh" này cũng quna trọng không kém việc (nếu như) LP1 còn được thao tác tiếp (như LP2, sau đây sẽ thấy) vì nó cho một cận dưới của mục tiêu tối ưu của IP do z = 23 đã đạt được (bởi nghiệm nguyên), z tối ưu không thể nhỏ hơn. Vậy ở đây ta đã đặt cận (bounding) cho thuật toán. Tiếp theo, nói chung ta phải giải LP2 như sẽ làm dưới đây. Nhưng trong trường hợp may mắn này ta thấy cận dưới 32 đã là 23 = [23.75] ([.] ký hiệu phần nguyên), với 23.75 là mục tiêu tối ưu của LP0 và cũng là của IP, đều nguyên, nên mục tiêu tối ưu của IP, là số nguyên không thể tốt hơn [23.75] = 23. Vậy ta đã may mắn được một nghiệm tối ưu của IP ngay sau một bước. Nếu ta không chọn (bất kỳ) LP1 như trên mà chọn LP2, thì dù đã chọn (bất kỳ và may mắn một phần) biến chia nhánh là x1, thuật toán kéo dài như sau. Giải LP2 (là quy hoạch có được khi thay ràng buộc x1 ≤ 3 bởi x1 ≥ 4 ở LP1) bằng thuật toán đơn hình ta được x1 = 4, x2 = 0.8333, z = 23, 3333. Vì x2 = 0.8333 chưa nguyên, ta phải nhánh LP2 để được LP3, LP4 với ràng buộc bổ sung tương ứng x2 ≤ 0 và x2 ≥ 1, tức là: 1. Miền chấp nhận được LP3 là miền chấp nhận được LP0 thêm ràng buộc x1 ≥ 4, x2 ≤ 0; 2. Miền chấp nhận được LP4 là miền chấp nhận được LP0 thêm ràng buộc x1 ≥ 4, x2 ≥ 1. Đến đây ta thấy hợp cả ba miền chấp nhận được của LP1, LP3 và LP4 chứa miền chấp nhận được của IP và còn phải xét cả ba. Có thể chọn bất kỳ để làm tiếp. Một trong các cách chọn hay dùng là chọn các bài toán con được lập sau cùng, ở đây là LP3 và LP4. Như vậy khi giải bằng phương pháp đơn hình có thể tận dụng kết quả ở lời giải LP2 ngay trước đó, vì so với LP2, ở LP3 và LP4 chỉ có thêm một ràng buộc. Do đó có thể áp dụng nghiên cứu hậu tối ưu lên nghiệm tối ưu của LP2 mà không cần giải lại (xem [Khánh - Nương 2000], Chương 8). Giả sử chọn ngẫu nhiên LP4 và giải ta thấy LP4 không có nghiệm chấp nhận được. Tất nhiên sau này chẳng cần nhìn đến LP4 nữa, ta cũng nói nó đã chạm đáy. Tiếp theo còn LP3 và LP1, ta lại chọn bất kỳ, nhưng theo nguyên tắc mới lập hơn thì ta lấy LP3 và giải được nghiệm tối ưu của LP3 là x1 = 4.5, x2 = 0, z = 22.5. Vì x1 chưa nguyên ta phải chia nhánh để được bài toán con: 1. Miền chấp nhận được LP5 là miền chấp nhận được LP0 thêm ràng buộc x1 ≥ 4, x2 ≤ 0, x1 ≤ 4; 2. Miền chấp nhận được LP4 là miền chấp nhận được LP0 thêm ràng buộc x1 ≥ 4, x2 ≤ 0, x1 ≥ 5. 33 Ta còn bài toán con chưa xét LP1, LP5 và LP6, chọn bất kỳ theo nguyên tắc mới hơn xét LP6 ta thấy nó không chấp nhận được, tức là đã chạm đáy. Đến lượt giải LP5, mới hơn LP1, ta được nghiệm tối ưu là x1 = 4, x2 = 0, z = 20 thoả cả ràng buộc về tính nguyên. Đến đây ta được một cận dưới của mục tiêu tối ưu của IP là z = 20 và do đó LP5 cũng chạm đáy, cung cấp một cận dưới cũng là kỷ lục là z = 20. Còn lại duy nhất LP1 phải xét (xem có được kỷ lục mới không). Theo trên đã thấy giải LP1 ta được nghiệm tối ưu x1 = 3, x2 = 2, z = 23, cũng thoả tính nguyên nên ta được một cận dưới và LP1 đã chạm đáy, cung cấp một kỷ lục mới. Vì không còn bài toán con nào để xét, kỷ lục mới này là nghiệm tối ưu của IP. Như vậy cũng với việc chọn x1 là biến phân nhánh nhưng chọn LP2 để xét thay cho LP1 ta đã phải làm dài nhất trong các khả năng có thể (nếu dùng quy tắc chọn trong các bài toán con mới lập) như H 2.3 minh hoạ sơ đồ cây thực hiện thuật toán. 2.2 Chia nhánh, đặt cận, chạm đáy Mọi thuật toán nhánh - và - cận đều gồm ba thủ tục chia nhánh, đặt cận, chạm đáy. Mỗi thủ tục được linh động trong cách thực hiện nên có các dạng thuật toán khác nhau. "Chia nhánh" là chọn một trong các bài toán con còn phải xét để chia tiếp ra. "Chọn" có thể lấy bất kỳ. Một quy tắc phổ biến là chọn trong các bài toán con mới lập. Vì việc đặt cận cũng có thể khác nhau, có nhiều cách định nghĩa cận (không chỉ là kỷ lục như định nghĩa trên), nên cũng sinh ra cách chọn khác là chọn trong các bài toán con có cận tốt nhất, với hy vọng sớm đi đến kỷ lục mới. "Chia" có nghĩa phổ biến là lấy một biến chia nhánh xj để chia ra hai bài toán con với ràng buộc thêm xj ≤ [a] và xj ≥ [a] + 1 tương ứng. Với bài toán nhị phân thì hai ràng buộc thêm sẽ là xj = 0 và xj = 1. Tuy vậy, một số thuật toán chi tiết hơn có thể dùng các chiến thuật nào đó để lấy biến chia nhánh. "Đặt cận" thường thực hiện bằng giải bài toán nới lỏng. Như phương pháp trên đây thì phải dùng LP - nới lỏng để đặt cận. Với các thuật toán nhánh - và - cận khác, còn có các kiểu bài toán nới lỏng khác. Chẳng hạn thuật toán nhánh - và - cận cho quy hoạch 34 Hình 2.3: nhị phân có thể áp dụng Lagrange nới lỏng (Lagrange relaxation). Định nghĩa tổng quát Lagrangian nới lỏng như sau. Giả sử cần xét quy hoạch nguyên max z = cTx, Ax < b, xj ∈ Z+, j = 1, . . . , n. Ở đây x ∈ Rn, A cấp mxn. Lấy một λ ≥ 0, λ ∈ Rm cố định. Khi đó Lagrangian nới lỏng được định nghĩa là bài toán max x≥0 L = cTx− λT (Ax− b). 35 Rõ ràng ta có quan hệ hai mục tiêu tối ưu z∗ ≤ L∗, tức là L∗ là một cận trên (có sát z∗ hay không còn phụ thuộc λ đã chọn) "Chạm đáy" có ba tiêu chuẫn như sau (cho bài toán tìm max): • Tiêu chuẩn 1: Nếu cận của bài toán con ≤ mục tiêu kỷ lục, hoặc nói chung hơn là nếu mọi mục tiêu chấp nhận được của bài toán con ≤ mục tiêu kỷ lục, thì bài toán này chạm đáy. • Tiêu chuẩn 2: Nếu bài toán con không chấp nhận được thì nó chạm đáy. • Tiêu chuẩn 3: Nghiệm tối ưu của LP - nới lỏng của bài toán con là nguyên thì nó chạm đáy. Rõ ràng nếu bài toán là tìm min thì ở Tiêu chuẩn 1 chỉ việc đổi dấu thành ≥, còn hai tiêu chuẩn kia giữ nguyên dạng. 2.3 Thuật toán nhánh - và - cận Xét bài toán quy hoạch nguyên tìm max. Gọi z là cận dưới - kỷ lục. Bước khởi động đặt i = 0 và đặt z = −∞. Ở một bước lập tổng quát có hai đoạn như sau. Đoạn 1 (Chạm đáy/đặt cận). Chọn LPi là bài toán con sẽ xét. Giải và xét nó có chạm đáy hay không. (a) Nếu LPi chạm đáy theo tiêu chuẩn 3 thì cập nhật z để được cận dưới mới. Trái lại, theo tiêu chuẩn khác, thì chọn bài toán con i mới và lập như trên. Nếu không còn bài toán con chưa xét thì dừng. Nghiệm tối ưu của Ip là tương ứng với cận dưới - kỷ lục cuối cùng, nếu có nghiệm như vậy. (b) Nếu LPi không chạm đáy, chuyển sang Đoạn 2 để chia nhánh LPi. 36 Đoạn 2 (Chia nhánh). Chọn một trong các biến xj mà x ∗ j ở nghiệm tối ưu LPi chưa là số nguyên. Lập hai bài toán con LP - nới lỏng bằng cách thêm ràng buộc tương ứng xj ≤ [x∗j ] và xj ≥ [x∗j ] + 1 Chuyển sang Đoạn 1. NHẬN XÉT. Nếu quy hoạch là tim min, thì ở thuật toán chỉ việc thay cận dưới z bởi cận trên z và ở bước khởi động đặt z = +∞. Các chi tiết khác không đổi. Thuật toán trên cũng áp dụng cho quy hoạch nguyên bộ phận. Nếu biến nào không bị ràng buộc lấy giá trị nguyên thì không được chọn là biến phân nhánh. Giả sử ở ví dụ 2.1, bây giờ x2 là biến liên tục. Nếu ta chọn LP1 để giải được nghiệm tối ưu x1 = 3, x2 = 2, z = 23 thì đây là nghiệm chấp nhận được của quy hoạch nguyên bộ phận MIP. Tuy nhiên ta không thể lam chạm đáy LP2 để được nghiệm vừa có là tối ưu của MIP bằng lý luận như ở Mục 2.1 vì z của MIP có thể không nguyên. Ta phải giải LP2 để được nghiệm như ở H 2.3. Nghiệm này chấp nhận được cho MIP nên LP2 chạm đáy theo Tiêu chuẩn 3. z = 23.3333 của LP2 là kỷ lục mới, và là nghiệm tối ưu của MIP vì không còn bài toán con nào để xét. Sơ đồ khối thuật toán: 37 38 Chương 3 Giải bài toán nguyên trên Mathlab 3.1 Giới thiệu Mathlab Matlab là chương trình phần mềm trợ giúp cho việc tính toán và hiển thị. Matlab có thể chạy trên hầu hết các hệ máy tính từ máy tính cá nhân đến máy tinh khổng lồ - super computer. Matlab được điều khiển bởi tập các bộ lệnh, tương tác bằng bàn phím trên cửa sổ điều khiển, đồng thời Matlab còn cho phép khả năng lập trình với cú pháp thông dịch lệnh hay còn gọi là scrift file. Các lệnh, bộ lệnh của Matlab lên đến con số hàng trăm và ngày càng được mở rộng bởi các toolbox trợ giúp hay các hàm ứng dụng tạo ra bởi người sử dụng. Các lệnh của Matlab rất mạnh và hiệu quả cho phép giải các loại hình bài toán khác nhau và đặc biệt hiệu quả cho hệ phương trình tuyến tính cũng như các thao tác trên các bài toán ma trận. Không những thế Matlab còn rất hữu hiệu trong việc trợ giúp thao tác và truy xuất đồ hoạ trong không gian 2D, 3D cũng như khả năng tạo hoạt cảnh minh hoạ cho việc mô tả bài toán một cách sinh động. Cùng với trên 25 toolbox (thư viện trợ giúp) khác nhau Matlab đưa đến cho các bạn sự lựa chọn hoàn chỉnh và phong phú về các công cụ trợ giúp đắc lực trong các lĩnh vực khác nhau trên con đường nghiên cứu mà bạn đã chọn. 39 Đây là một số lĩnh vực mà Matlab đang giải quyết: • Nghiên cứu và phát triển trong lĩnh vực công nghiệp • Giảng dạy, nghiên cứu lập các chương trình ứng dụng trong giảng dạy cho các môn như toán, lý, hoá, ... trong các trường phổ thông nhằm nâng cao khả năng tiếp thu cũng như ý sáng tạo trong học sinh • Giảng dạy, nghiên cứu lập các chương trình về toán đặc biệt là các loại hình nguyên lý cơ bản và các phương trình tuyến tính cho sinh viên cũng như học sinh các trường kỹ thuật • Giảng dạy và nghiên cứu trong lĩnh vực kỹ thuật và khoa học như: điện tử, lý thuyết điều khiển, vật lý, đồ hoạ, xử lý ảnh, vật liệu, ... • Giảng dạy và nghiên cứu trên mọi lĩnh vực có xuất hiện tính toán bao gồm toán kinh tế, hoá, cơ học, sinh học, ... 3.2 Lập trình thuật toán giải bài toán nguyên Thuật toán cần áp dụng cùng với Optimazition toolbox và thuật toán đơn hình gốc (linprog.m) trong Matlab: function[x, val, status] = IP1(f, A, b, Aeq, beq, lb, ub,M, e) options = optimset(′display′,′ off ′); bound = inf ; Khởi tạo chặn trên của biến [x0, val0] = linprog(f, A, b, Aeq, beq, lb, ub, [], options); [x, val, status, b] = rec(f, A, b, Aeq, beq, lb, ub, x0, val0,M, e, bound); 40 Hàm đệ quy truy xuất cây nhánh cận function[xx, val, status, bb] = rec(f, A, b, Aeq, beq, lb, ub, x, v,M, e, bound) options = optimset(′display′,′ off ′); x là nghiệm tối ưu khởi tạo và v là giá trị tối ưu tương ứng Giải mô hình bài toán tuyến tính sau khi bỏ ràng buộc nguyên [x0, val0, status0] = linprog(f, A, b, Aeq, beq, lb, ub, [], options); Nếu nghiệm không chấp nhận được hay giá trị mục tiêu cao hơn giá trị chặn, trả về nghiệm khởi tạo ifstatus0 bound xx = x; val = v; status = status0; bb = bound; return; end Nếu biến ràng buộc nguyên là nguyên ind = find(abs(x0(M)− round(x0(M))) > e); ifisempty(ind) status = 1; ifval0 < bound nghiệm này tốt hơn nghiệm hiện tại thì thay thế x0(M) = round(x0(M)); xx = x0; val = val0; bb = val0; else 41 xx = x; trả về nghiệm khởi tạo val = v; bb = bound; end return end Nếu nghiệm của LP relaxation là chấp nhận được và đưa một giá trị nhỏ hơn chặn hiện tại nhưng vài biến ràng buộc nguyên là không nguyên. Vì vậy chúng ta lọc biến không nguyên và từ hay bài toán LP và giải chúng đệ quy bằng cách gọi hàm tương tự (rẽ nhánh) Bài toán LP1 với ràng buộc Xi <= floor(Xi) , i=ind(1) brvar = M(ind(1)); brvalue = x(brvar); ifisempty(A) [rc] = size(Aeq); else [rc] = size(A); end A1 = [A; zeros(1, c)]; A1(end, brvar) = 1; b1 = [b; floor(brvalue)]; Bài toán LP với ràng buộc Xi >= ceil(Xi) , i=ind(1) A2 = [A; zeros(1, c)]; A2(end, brvar) = −1; b2 = [b;−ceil(brvalue)]; 42 giải LP1 [x1, val1, status1, bound1] = rec(f, A1, b1, Aeq, beq, lb, ub, x0, val0,M, e, bound); status = status1; if status1 > 0 và bound1 < bound nếu nghiệm tối ưu thì đưa vào chặn tốt hơn xx = x1; val = val1; bound = bound1; bb = bound1; else xx = x0; val = val0; bb = bound; end Giải LP2 [x2, val2, status2, bound2] = rec(f, A2, b2, Aeq, beq, lb, ub, x0, val0,M, e, bound); if status2 > 0 và bound2 < bound nếu nghiệm tối ưu thì đưa vào chặn tốt hơn status = status2; xx = x2; val = val2; bb = bound2; end 43 3.3 Giải bài toán nguyên trên Mathlab Ví dụ 3.3.1 Xét bài toán "knapsack problem" ở phần trên bằng phương pháp nhánh cận trên Matlab. Mô hình bài toán: max z = x1 + x2, 3x1 + 2x2 ≤ 12, x2 ≤ 2, x1, x2 ≥ 0 và nguyên. Chuyển thành ngôn ngữ Matlab, ta có:, f = [−1,−1]; lấy phủ định của bài toán max A = [3 2; 0 1]; B = [12; 2]; lb = [0 0]; ub = [inf inf ]; M = [1, 2]; vecto chỉ số cho biến mà có ràng buộc nguyên 44 45 Nghiệm tối ưu của bài toán: x1 = 4, x2 = 0 với z = 4 Ví dụ 3.3.2 Xét bài toán "Capital Budgeting" ở phần trên bằng phương pháp nhánh cận trên Matlab . Mô hình bài toán: max z = 8x1 + 11x2 + 6x3 + 4x4, 5x1 + 7x2 + 4x3 + 3x4 ≤ 14, xj ∈ 0, 1, j = 1, . . . , 4. 46 Chuyển thành ngôn ngữ Matlab, ta có:, f = [−8,−11,−6,−4]; lấy phủ định của bài toán max A = [5 7 4 3; 1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]; B = [14; 1; 1; 1; 1]; lb = [0 0 0 0]; ub = [1 1 1 1]; M = [1, 2, 3, 4]; vecto chỉ số cho biến mà có ràng buộc nguyên 47 Nghiệm tối ưu của bài toán: x1 = 0, x2 = x3 = x4 = 1 với z = 21 48 TÀI LIỆU THAM KHẢO 1. [Beasley 1996] (J. E. Beasley, editor), Advances in Linear and Integer Programminh, Clarendon Press, Oxford. 2. Belvaux, G. and Wolsey, L. A., 2001, Modelling Practical Lot-Sizing Problems as Mixed-Integer Programs. Management Science 47, 993-1007. 3. [Bertsimas - Tsitsiklis 1997] (D. Bertsimas, J. N. Tsitsiklis), Introduction to Linear Optimization, Athena Scientific, Belmont, Massachusetts. 4. [Bixby - Boyd - Rios Mercado 1998] (R. E. Bixby, E. A. Boyd, R. Z. Rios Mercado) Integer Programming and combinatorial Optimization, Springer, New York - Berlin. 5. [Danh 2003] (T. N. Danh), Toán rời rạc, Nhà xuất bản Đại học Quốc gia TP. Hồ Chí Minh. 6. [Gomory 1958] (R. E. Gomory), Outline of an algorithm for integer solutions to linear programs, Bulletin of the AMS, Vol. 64, 275-278. 7. [Hanselman - Littlefield] (Duane Hanselman, Bruce Littlefield), Mastering Matlab 7, University of Maine, New York. 8. [Hillier - Lieberman 2001] (F. S. Hillier, G. J. Lieberman), Introduction to Operation Reasearch, McGraw-Hill, New York. 9. [Khánh - Nương 2000] (P. Q. Khánh, T. H. Nương), Quy hoạch tuyến tính, Nhà xuất bản Giáo dục, Hà Nội. 10. [Khánh 2000] (P. Q. Khánh), Toán chuyên đề, Nhà xuất bản Đại học Quốc gia TP. Hồ Chí Minh. 11. [Khánh 2002] (P. Q. Khánh), Vận trù học, Nhà xuất bản Giáo dục, Hà Nội. 12. João P. Pedroso. An evolutionary solver for linear integer programming. BSIS Tehnial Report 98-7, Riken Brain Siene Institute, Wako-shi, Saitama, Japan, 1998. 13. [Nghĩa 1997] (N. Đ. Nghĩa), Tối ưu hoá, Nhà xuất bản Giáo dục, Hà Nội. 49 14. [Nghĩa - Thành 1997] (N. Đ. Nghĩa, N. T. Thành), Toán rời rạc, Nhà xuất bản Khoa học Kỹ thuật. 15. Robert E. Bixby, Sebastian an Ceria, CassandraM. M Zeal, and MartinW. P. Savels- bergh. An updated mixed integer programming library. Tehnial report, Rie Univer- sity, 1998. TR98-03. 16. [Schrijver 1999] (A. Schrijver), Theory of Linear and Integer Programming, John Wiley and Sons, New York. 17. S. Sen. Algorithms for Stochastic Mixed-Integer Programming Models. Technical Report, University of AZ. 2004. 18. [Taha 1975] (H. Taha), Integer Programming, Theory, Applications and Computa- tions, Academic Press, New York. 19. [Taha 1997] (H. Taha), Operations Research, An Introduction, McMillan, New York. 50