Giải thuật lpsolve và ứng dụng quy hoạch tuyến tính

1. Cơ sở Toán học (Mathematical Foundations)

1.1. Định nghĩa bài toán

Quy hoạch tuyến tính (LP) là phương pháp tối ưu hóa một hàm mục tiêu tuyến tính, thỏa mãn các ràng buộc đẳng thức hoặc bất đẳng thức tuyến tính. Mô hình toán học tổng quát (Canonical Form):

MaximizecTxSubject toAxbx0\begin{aligned} & \text{Maximize} & & \mathbf{c}^\mathrm{T} \mathbf{x} \\ & \text{Subject to} & & A \mathbf{x} \leq \mathbf{b} \\ & & & \mathbf{x} \geq 0 \end{aligned}

Trong đó:

  • x\mathbf{x} : Vector biến quyết định (variables).
  • c\mathbf{c} : Vector hệ số hàm mục tiêu (objective coefficients).
  • AA : Ma trận hệ số ràng buộc (constraint matrix).
  • b\mathbf{b} : Vector giới hạn (bounds).

1.2. Phương pháp giải (Algorithms)

  1. Simplex Method (Đơn hình): Phát triển bởi George Dantzig (1947). Di chuyển dọc theo các cạnh của đa diện chấp nhận được (feasible polytope) để tìm cực trị tại các đỉnh.
  2. Interior Point Method (Điểm trong): Tiếp cận nghiệm tối ưu bằng cách di chuyển qua phần trong của vùng chấp nhận được (Karmarkar, 1984).
  3. Branch and Bound (Nhánh và Cận): Sử dụng cho quy hoạch nguyên (Integer Programming - MILP) khi các biến x\mathbf{x} phải là số nguyên.

2. Mô hình hóa 3 Bài toán Thực tế (Case Studies)

Bài toán 1: Quản lý Sản xuất (Production Mix)

Mục tiêu: Tối đa hóa lợi nhuận từ việc sản xuất 2 loại sản phẩm P1,P2P_1, P_2 .

  • Hàm mục tiêu: Max Z=30x1+40x2Z = 30x_1 + 40x_2
  • Ràng buộc nguyên liệu: 1x1+2x21001x_1 + 2x_2 \leq 100
  • Ràng buộc nhân công: 2x1+1x2802x_1 + 1x_2 \leq 80
  • x1,x20x_1, x_2 \geq 0

Bài toán 2: Vận chuyển (Transportation Problem)

Mục tiêu: Tối thiểu hóa chi phí vận chuyển từ 2 kho (S1, S2) đến 2 cửa hàng (D1, D2).

  • Cung (Supply): S1=50,S2=60S_1=50, S_2=60
  • Cầu (Demand): D1=40,D2=70D_1=40, D_2=70
  • Ma trận chi phí CC: (2431)\begin{pmatrix} 2 & 4 \\ 3 & 1 \end{pmatrix}
  • Hàm mục tiêu: Min Z=2x11+4x12+3x21+1x22Z = 2x_{11} + 4x_{12} + 3x_{21} + 1x_{22}

Bài toán 3: Quản lý Dự án (Project Management / CPM)

Mục tiêu: Tối thiểu hóa thời gian hoàn thành dự án gồm 3 công việc A, B, C.

  • A (3 ngày), B (2 ngày) làm song song; C (4 ngày) làm sau A và B.
  • Biến: t0,tA,tB,tCt_0, t_A, t_B, t_C (thời điểm bắt đầu/kết thúc).
  • Hàm mục tiêu: Min TendT_{end}
  • Ràng buộc: tAt03,tBt02,tCtA4,tCtB4t_A - t_0 \geq 3, t_B - t_0 \geq 2, t_C - t_A \geq 4, t_C - t_B \geq 4

3. Triển khai và Tính toán (Implementation)

3.1. Python (Thư viện SciPy)

Cài đặt: pip install scipy numpy matplotlib

Triển khai: Sử dụng scipy.optimize.linprog. Lưu ý: SciPy mặc định tìm Min, để tìm Max cần đảo dấu hàm mục tiêu.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog

# Bài toán 1
c_prod = [-30, -40]
A_prod = [[1, 2], [2, 1]]
b_prod = [100, 80]
res_prod = linprog(c_prod, A_ub=A_prod, b_ub=b_prod, bounds=(0, None), method='highs')

print(f"Lợi nhuận tối đa: {-res_prod.fun:.2f}")

3.2. R (Thư viện lpSolveAPI)

⚠️Lưu ý quan trọng về bộ nhớ trong R

Đối tượng mô hình lpSolve tạo bởi make.lp là con trỏ ngoài (external pointers) tới cấu trúc C lprec. R không thể tự động nhân bản chúng. Sử dụng rm() để giải phóng bộ nhớ.

library(lpSolveAPI)
lprec <- make.lp(0, 2)
set.objfn(lprec, c(30, 40))
lp.control(lprec, sense='max')
add.constraint(lprec, c(1, 2), "<=", 100)
add.constraint(lprec, c(2, 1), "<=", 80)
solve(lprec)
get.objective(lprec)
rm(lprec)

3.3. Octave (Sử dụng GLPK)

c = [30; 40];
A = [1, 2; 2, 1];
b = [100; 80];
[x_opt, z_max] = glpk(c, A, b, [0;0], [], "UU", "CC", -1);

3.4. Julia (Thư viện JuMP & GLPK)

using JuMP, GLPK
model = Model(GLPK.Optimizer)
@variable(model, x[1:2] >= 0)
@objective(model, Max, 30x[1] + 40x[2])
@constraint(model, c1, x[1] + 2x[2] <= 100)
@constraint(model, c2, 2x[1] + x[2] <= 80)
optimize!(model)

4. Phân tích độ nhạy (Sensitivity Analysis)

Trong thực tế, các tham số đầu vào thường xuyên biến động. Phân tích độ nhạy giúp trả lời câu hỏi: "Nghiệm tối ưu sẽ thay đổi thế nào nếu điều kiện biên thay đổi?"

4.1. Các khái niệm then chốt

  1. Giá mờ (Shadow Price / Dual Value):* Là giá trị tăng thêm của hàm mục tiêu khi ta nới lỏng một đơn vị ràng buộc.
  2. Biến bù/dư (Slack/Surplus):
    • Slack: Lượng tài nguyên chưa sử dụng hết (trong ràng buộc \leq).
    • Surplus: Lượng vượt mức tối thiểu (trong ràng buộc \geq).
  3. Giá giảm (Reduced Cost): Cho biết hàm mục tiêu sẽ thay đổi bao nhiêu nếu ta ép buộc một biến đang bằng 0 phải tham gia vào phương án.

4.2. Trích xuất dữ liệu độ nhạy

Python (SciPy):

# Shadow Prices của các ràng buộc
shadow_prices = res_prod.ineqlin.marginals
# Biến bù (Slack)
slacks = res_prod.slack

Julia (JuMP):

# Shadow Price của ràng buộc c1
sp_nguyen_lieu = shadow_price(c1)
# Reduced Cost của biến x1
rc_x1 = reduced_cost(x[1])

4.3. Diễn giải dưới góc độ Quản trị

📄abstract

Operational Insights Nếu Giá mờ* của "Nhân công" là 15.0:

  • Doanh nghiệp nên sẵn sàng thuê thêm giờ làm với giá < 15.0 để tăng lợi nhuận.
  • Nếu một tài nguyên có Biến bù (Slack) > 0, thì Giá mờ* của nó luôn bằng 0 (tài nguyên dư thừa thì việc có thêm nó không làm tăng lợi nhuận).

(*) Giá mờ: là một ước tính giá trị tiền tệ được gán cho các hàng hóa, dịch vụ hoặc nguồn lực không có giá thị trường rõ ràng, phản ánh giá trị thực hoặc chi phí cơ hội của chúng.


5. Biểu đồ hóa vùng chấp nhận (Visualization)

Biểu đồ thể hiện vùng nghiệm khả thi (Feasible Region) cho Bài toán 1 (Sản xuất).

import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 100, 200)
y1 = (100 - x) / 2
y2 = 80 - 2 * x

plt.figure(figsize=(8, 6))
plt.plot(x, y1, label=r'$x_1 + 2x_2 \leq 100$')
plt.plot(x, y2, label=r'$2x_1 + x_2 \leq 80$')
plt.fill_between(x, 0, np.minimum(y1, y2), where=(np.minimum(y1, y2)>=0), color='gray', alpha=0.3)
plt.xlim((0, 60))
plt.ylim((0, 60))
plt.legend()
plt.title('Vùng nghiệm chấp nhận được - Bài toán Sản xuất')
plt.grid(True)
plt.show()

6. Nguồn trích dẫn (References)

  1. Dantzig, G. B. (1963). Linear Programming and Extensions. Princeton University Press.
  2. Berkelaar, M., et al. (2024). lpSolveAPI: Interface to 'Lp_solve'.
  3. Virtanen, P., et al. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.
  4. Dunning, I., Huchette, J., & Lubin, M. (2017). JuMP: A Modeling Language for Mathematical Optimization.