R và thuật toán điểm trong cho quy hoạch tuyến tính

307 KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN Tóm tắt Ngôn ngữ mã nguồn mở R là môi trường tốt cho tính toán thống kê và phân tích dữ liệu. R cũng là công cụ cực tốt để tiếp cận Toán, và sử dụng Toán cho các ứng dụng của người dùng. Bài viết dùng R cho việc giải quy hoạch tuyến tính bằng thuật toán điểm trong. Khía cạnh được đưa ra trong bài viết là thực thi thuật toán bằng ngôn ngữ R, cả phía người dạy và người học. Từ khóa: Quy hoạch tuyến tính, thuật toán điểm trong, Affine scaling, barrier, R, dạy và học. 1. GIỚI THIỆU Đơn hình (simplex method, SM) là phương pháp phổ biến để giải các quy hoạch tuyến tính, đến khi phương pháp điểm trong phát triển. Với SM, việc đạt tới nghiệm tối ưu của quy hoạch tuyến tính thực hiện qua hành động chuyển từ đỉnh này tới đỉnh khác dọc theo các cạnh của miền ràng buộc đa diện, theo hướng thay đổi hàm mục tiêu. Điểm trong, như tên gọi, đi từ điểm nằm hẳn bên trong tập chấp nhận được, không nằm trên biên, theo hướng tốt dần để tới nghiệm tối ưu. Karmarkar (1984), được xem là người tiên phong của lĩnh vực phương pháp điểm trong. Phương pháp Projective scaling của Karmarkar có thể cạnh tranh với SM khi áp dụng vào các bài toán thực tế. Thuật toán điểm trong mượn những ý đơn giản từ R VÀ THUẬT TOÁN ĐIỂM TRONG CHO QUY HOẠCH TUYẾN TÍNH 29. ThS. Phạm Việt Huy Trường Đại học Tài chính - Marketing KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 308 tối ưu phi tuyến, và với quy hoạch tuyến tính thì chúng đủ cơ bản để phát triển. Với quy hoạch tuyến tính, phương pháp điểm trong đạt đến nghiệm tối ưu (giờ được biết là nằm trên biên của miền chấp nhận được) thông qua một dãy các điểm trong. Kết quả các phép lặp của phương pháp điểm trong không nằm trên biên, mà là các điểm trong thật sự của miền chấp nhận được. Có nhiều hướng đi đến các điểm trong kế tiếp từ điểm trong khởi đầu, tập trung vào hai nhóm chính. Một là thay đổi tỷ lệ để giữ các điểm hiện tại cách xa biên, và hạn chế bước nhảy để không chạm biên. Hướng này có trong phương pháp affine scaling của Dikin. Hai là thêm lực đẩy biên vào hàm mục tiêu để chặn di chuyển đến biên của miền ràng buộc. Hướng này được biết là phương pháp barrier. Chính ngay trong từng hướng, đã có sự cải tiến từng chút từng chút một, làm cho các phương pháp điểm trong càng ngày càng tinh tế. Chẳng hạn, Boah D.K và Twum S.B. (2019) đưa ra một thay đổi cho phương pháp affine scaling. Thay đổi rõ từ nghiên cứu này là số lượng các bước lặp ứng với từng lựa chọn của tỷ lệ. [2] Boah D.K và Twum sử dụng MatLab cho các tính toán của mình. Và ngay MatLab cũng có gói linprog để giải quy hoạch tuyến tính, mà lựa chọn điểm trong đã xuất hiện. Gói linprog đã xuất hiện trong R, và cùng với lpSolve, là công cụ đắc lực để giải quy hoạch tuyến tính. Bài viết tập trung vào hướng thực hiện thuật toán điểm trong gốc đối ngẫu affine scaling và thuật toán điểm trong barrier cho quy hoạch tuyến tính. Trước đó, chúng tôi tóm lại ý tưởng của các thuật toán và các kết quả liên quan, có tham khảo tài liệu [3], [4], [6]. Bài viết dùng ngôn ngữ R minh họa mã hóa đoạn chính thể hiện bước lặp chuyển sang điểm trong mới. Với hướng này, người dạy và người học đọc ra sự di chuyển của dãy điểm trong, có thể điều chỉnh theo ý muốn, và nắm lại ý tưởng của từng thuật toán. Hơn nữa, thực hiện R như một ngôn ngữ lập trình cũng là một điều phù hợp. Ở Việt Nam, thuật toán điểm trong đã được giới thiệu và nghiên cứu từ lâu. Riêng cho quy hoạch tuyến tính giải bằng thuật toán điểm trong, chúng tôi có tham khảo một sản phẩm từ tác giả Nguyễn Thị Hồng Lê và Trần Vũ Thiệu (2011) [7]. Hướng thực thi là phần để mở phía sau của tài liệu [7]. Thuật toán đơn hình, hay thuật toán điểm trong, ngay khi xuất hiện đã đi kèm luôn nội dung thực thi. Hoạt động thực thi được hỗ trợ cực mạnh từ nhiều nguồn lực. [5] cho một tổng quan về công cụ giải quyết quy hoạch tuyến tính bằng thuật toán điểm trong, cả nguồn mở và nguồn thương mại. 309 KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 2. THUẬT TOÁN ĐIỂM TRONG VÀ QUY HOẠCH TUYẾN TÍNH 2.1. Bài toán Bài toán được xem xét trong bài viết là: mincT Ax = b x ≥ 0, x là vector n thành phần. Đối ngẫu của bài toán này: maxbT AT π ≤ c. Điều kiện độ lệch bù phát biểu rằng x* và π* là tối ưu nếu: - x* là chấp nhận được của gốc: Ax* = b, x* ≥ 0; - π* là chấp nhận được của đối ngẫu: AT π* ≤ c; - Hoặc x j * = 0 hoặc a j T π j * = c j . Khi lấy σ ≥ 0 là biến bù từ ràng buộc AT π ≤ c, AT π + σ = c, thì a j T π j * = c j trở thành x j * = 0 hoặc σ j * = 0 hoặc cả hai, tức là x j *σ j * = 0. Viết lại điều kiện độ lệch bù: x j * và π* tối ưu nếu: Ax* = b, x* ≥ 0; AT π*+σ* = c, σ* ≥ 0; x j * π j * = 0, với mọi j = 1,...,n. Hệ x j σ j = 0 viết ở dạng ma trận là XΣ = 0, với X và Σ đều là các ma trận chéo, X = diag(x), [X] j,j = x j , và Σ = diag(σ), [Σ] j,j = σ j . Tối ưu gốc cho x và tối ưu đối ngẫu cho π giờ tương đương giải (1), là hệ sau: KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 310 Thuật toán affine scaling dùng (1). Thuật toán barrier điều chỉnh thành , để có hệ (2) sau đây: Điểm trong của hệ ràng buộc, , thỏa , . Điểm trong này dùng cho cả thuật toán affine scaling và thuật toán barrier. 2.2. Ý tưởng thuật toán affine scaling - Khi đã có điểm trong là nghiệm của hệ: đi tìm nghiệm và thỏa , mà và không là điểm trong, chỉ cần , và ít nhất một trong hai bằng 0, với mọi . - Từ , tìm bước nhảy để đến điểm tốt hơn, theo nghĩa gần hơn với việc thỏa điều kiện độ lệch bù. Để tìm bước nhảy này, thay vào hệ (1) cho thuật toán affine scaling: thì có hệ (3) sau đây: Hệ (3) có ở vế phải thì phi tuyến. Lựa chọn bỏ đi ở vế phải là hợp lý, và từ phương trình cuối của (3) ta có: 311 KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN Thay vào phương trình giữa, thì . Như vậy, hệ bước nhảy là (4): - Hy vọng dùng như nghiệm mới. Có thể có một thành phần nào đó của hoặc âm, dù chúng phải dương để giữ mình trong miền ràng buộc. Do đó, lấy nghiệm mới , với là phân số dương không lớn hơn 1. Vì và dương nên dương nếu chọn đủ nhỏ. Với , dẫn tới . Tương tự, với . phải thỏa , và . Lấy , rồi , thì có nghiệm mới. - Theo cách này, có thể không tới được nghiệm tối ưu, bởi các thành phần của vector và đều dương, mà điều kiện độ lệch bù buộc phải có một hoặc một bằng 0. Tuy nhiên, thành phần bù tiến gần đến 0 khi có đủ số bước lặp được thực hiện. 2.3. Ý tưởng thuật toán barrier - Thay đổi thuật toán affine scaling theo hướng giữ các phép lặp gần đường trung tâm gốc để đến nghiệm tối ưu. Một, điều kiện độ lệch bù giờ là , với . Hai, bắt đầu với giá trị lớn của , rồi giảm dần về 0. - Khi là , có sự chuyển dịch (1) thành (2), còn hệ tìm các bước nhảy thành phần chuyển từ (3) sang (5): KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 312 Cũng bỏ đi để sắp xếp lại (5) thành (6): - Log barrier đến từ việc buộc thành phần trong nghiệm tối ưu tránh 0, xem tài liệu tham khảo [4]. Kết quả: thỏa hệ điều kiện độ lệch bù nới lỏng với với mọi nếu và chỉ nếu là tối ưu của bài toán log barrier sau: Lựa chọn đến từ nới lỏng , điều chỉnh bằng phân số dương bé hơn 1. 2.4. Thuật toán điểm trong affine scaling - Chuẩn bị điểm trong xuất phát , và chọn tỷ lệ chấp nhận được gốc , cùng với mức sai lệch giới hạn . Điểm trong xuất phát : - Lặp tới khi với mọi : Tìm và từ (4): và Cập nhật 2.5. Thuật toán barrier - Chuẩn bị điểm trong xuất phát , chọn tỷ lệ chấp nhận được gốc , tỷ lệ chắn , và mức sai lệch giới hạn . 313 KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN Điểm trong xuất phát : 3. THỰC HIỆN Bài viết lấy các ví dụ chỉ để minh họa việc thực thi các thuật toán điểm trong cho quy hoạch tuyến tính. Một ví dụ cho affine scaling. Một ví dụ cho barrier. Kết quả từng đoạn lặp không thể hiện. Các kết quả đã được so sánh sơ bộ với kết quả từ MatLab, dùng gói lệnh linprog với tùy chọn optimoptions là interior point hoặc interior point legacy hoặc dual simplex đều được. MatLab thực hiện với cả ràng buộc đẳng thức, và các biến có cận trên cận dưới hữu hạn. Viết mã R để đưa ra kết quả từng đoạn lặp, hay việc giải quyết cận trên cận dưới của biến của quy hoạch tuyến tính, bằng R, thì không là khó khăn. Cho mỗi ví dụ, chúng tôi đã lựa chọn một điểm khởi đầu thuật toán, . Việc tìm điểm khởi đầu cũng là nội dung để người dạy và người học tiến hành với R. Chúng tôi cũng đã thử nghiệm với các điểm khởi đầu thuật toán khác, cùng với các giá trị , , khác. Có lựa chọn khởi đầu làm số bước lặp rất lớn, và cũng có lựa chọn làm số bước lặp nhỏ hơn. Đơn hình kết thúc ở chính xác điểm tối ưu. Điều này còn khả thi cho phương pháp điểm trong không? Phương pháp điểm trong sinh ra dãy điểm trong hội tụ tới tối ưu. “Nhận diện điểm tối ưu như thế nào từ kết quả sau các bước lặp của thuật toán điểm trong?”, “Dùng được tính toán hình thức cho các thuật toán điểm trong này hay không?”, ... không bàn đến ở phần này. KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 314 3.1. Ví dụ 1 Chọn điểm trong bắt đầu thuật toán là , Lấy , và sai lệch giới hạn . Đoạn lặp để chuyển đến điểm trong mới được viết bằng R như sau: repeat { Xb=diag(xb) Sigmab=diag(sigmab) Dxp=solve(rbind(cbind(-solve(Xb)%*%Sigmab,t(A)), cbind(A,matrix(0,nrow(A),ncol(t(A))))), rbind(sigmab,matrix(0,nrow(A),1))) Dx=matrix(1,nrow(xb),1) for (i in 1:length(xb)) { Dx[i,1]=Dxp[i] } Dp=matrix(1,nrow(pib),1) for (i in 1:length(pib)) { Dp[i,1]=Dxp[length(xb)+i] } Dsigma=-sigmab-solve(Xb)%*%Sigmab %*% Dx temp_x=xb/(-Dx) 315 KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN theta_x=min(temp_x[temp_x>0]) temp_s=sigmab/(-Dsigma) theta_sigma=min(temp_s[temp_s>0]) theta_v=c(1,alpha*theta_x,alpha*theta_sigma) theta=min(theta_v) xb=xb+theta*Dx pib=pib+theta*Dp sigmab=sigmab+theta*Dsigma tol=max(xb*sigmab) if (tol<epsilon) { break } } Kết quả sau 9 bước lặp: Bước lặp 9 có , đã nhỏ hơn chọn trước. 3.2. Ví dụ 2 Ví dụ 2 đến từ tài liệu tham khảo của Boah D.K. và Twum S.B. (2019) [2]. Ở [2], đây là quy hoạch hai biến, dễ cho việc hình dung hình học về phương pháp. Kết quả của ví dụ này không có thành phần 0. Đồng thời các tác giả này đưa ra đánh giá về số bước lặp ứng với các lựa chọn giá trị khác nhau, là giá trị chọn trước trong liên hệ và của [2]. Ví dụ 2: KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 316 Ví dụ 2 được xử lý bằng thuật toán điểm trong barrier. Đoạn lặp chính viết bằng ngôn ngữ R. Kết quả của [2] thực hiện từ MatLab. [2] có đưa ra hướng cải tiến so với thuật toán affine scaling ban đầu. Chọn điểm trong bắt đầu thuật toán barrier là , Sai lệch giới hạn như trước. Ví dụ 2, nếu vẫn dùng thuật toán affine scaling, với xuất phát vẫn như trên, , , kết quả có sau 8 bước lặp như sau Sai lệch giới hạn ở bước 8 là 3.17045e-07. Đoạn mã R cho phần lặp chính của thuật toán barrier như sau repeat { mub=beta*(t(sigmab)%*%xb)/length(xb) Xb=diag(xb) Sigmab=diag(sigmab) ev=matrix(1,length(xb)) Dxp=solve(rbind(cbind(-solve(Xb)%*%Sigmab,t(A)), cbind(A,matrix(0,nrow(A),ncol(t(A))))), rbind(sigmab-mub[1,1]*(solve(Xb)%*%ev),matrix(0,nrow(A),1))) 317 KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN Dx=matrix(1,nrow(xb),1) for (i in 1:length(xb)) { Dx[i,1]=Dxp[i] } Dp=matrix(1,nrow(pib),1) for (i in 1:length(pib)) { Dp[i,1]=Dxp[length(xb)+i] } Dsigma=-sigmab-solve(Xb)%*%(Sigmab%*%Dx-mub[1,1]*ev) temp_x=xb/(-Dx) theta_x=min(temp_x[temp_x>0]) temp_s=sigmab/(-Dsigma) theta_sigma=min(temp_s[temp_s>0]) theta_v=c(1,alpha*theta_x,alpha*theta_sigma) theta=min(theta_v) xb=xb+theta*Dx pib=pib+theta*Dp sigmab=sigmab+theta*Dsigma tol=max(xb*sigmab) if (tol<epsilon) { break } } Kết quả, sau 9 bước lặp, và còn sai lệch giới hạn ở bước 9 là 2.159964e-06. KỶ YẾU HỘI THẢO KHOA HỌC ĐÀO TẠO NGÀNH TOÁN KINH TẾ TRONG BỐI CẢNH HIỆN NAY VÀ CÁC VẤN ĐỀ LIÊN QUAN 318 2. KẾT LUẬN Bài viết minh họa việc thực hiện thuật toán điểm trong cho quy hoạch tuyến tính với công cụ là ngôn ngữ lập trình R. Hai thuật toán điểm trong được thực thi là thuật toán affine scaling và thuật toán barrier. Chính ngay trong R đã có gói lệnh để giải quy hoạch tuyến tính. Tiếp cận của bài viết là sử dụng R để tái hiện nét chính của các thuật toán khi dạy và học về thuật toán điểm trong. TÀI LIỆU THAM KHẢO 1. Karmarkar, N. (1984), Polynomial-Time Algorithm for Linear Programming. Combinatorica, 4, 373-395. 2. Boah, D.K. and Twum, S.B. (2019), An Improved Affine - Scaling Interior Point Algorithm. Journal of Applied Mathematics and Physics, 7, 2531-2536. 3. Eiselt, H.A. and Sandblom, C.L. (2007), Linear Programming and Its Applications. Springer-Verlag, Berlin Heidelberg, 280. 4. Robert J. Vanderbei (2014), Linear Programming Foundations and Extensions. International Series in Operations Research Management Science, Vol 196, Springer. 5. Handey Y. Benson (2011), Interior-Point Linear Programming Solvers. Encyclopedia of Operations Research and Management Science, Editor James J. Cochran, John Wiley Sons. 6. Singh, J.N. and Singh, D. (2002), Interior-Point Methods for Linear Programming: A Review. International Journal of Mathematical Education in Science and Technology, 33, 405-423. 7. Nguyễn Thị Hồng Lê, Trần Vũ Thiệu (2011), Phương pháp điểm trong giải quy hoạch tuyến tính. Đại học Thái Nguyên.