Bài viết này sẽ nói về cách thiết lập và giải 1 biến thể bài toán người giao hàng (được "bịa" ra) dùng quy hoạch tuyến tính có số nguyên (MILP - mixed integer linear programming) qua solver SCIP & Julia.

Giới thiệu bài toán

Bài toán người giao hàng (Travelling Salesman)

  • Input: các thành phố và khoảng cách giữa chúng
  • Output: chu trình ngắn nhất thăm các thành phố đúng 1 lần

Biến thể trình bày ở bài viết này

Giống người giao hàng nhưng tìm số công nhỏ nhất. Cụ thể, các món hàng có cân nặng nhất định (hàng đến thành phố \(i\) có cân nặng \(c_i\)). Người giao hàng sẽ có cân nặng \(W + \sum_{i=1}^n c_i\) ban đầu với \(W\) là cân nặng của người đưa hàng (người đưa hàng phải mang tất cả hàng + cân nặng của anh ta). Khi đến thành phố \(i\) thì có thể bỏ lại món hằng nặng \(c_i\). Số công được tính bằng \(\sum\) (cân nặng \(\times\) khoảng cách). Trong thực tế thì không chỉ có \(1\) người giao hàng mà sẽ có \(m\) người.

  • Input: các thành phố, khoảng cách giữa chúng, cân nặng các món hàng ứng với từng thành phố
  • Output: chu trình với số "công" tối thiểu (người giao hàng khá lười, làm xong nhanh không quan trọng bằng làm tốn ít công sức nhất)

Phát biểu dưới dạng toán học

Tham khảo: Travelling_salesman_problem#/Miller-Tucker-Zemlin_formulation

Tương tự với bài toán TSP gốc, gọi

  • \(x_{i j} \in \{0,1\}\) là biến cho biết có chọn cạnh đi từ \(i\) tới \(j\) hay không
  • \(d_{i j} = d_{j i}\) là khoảng cách giữa thành phố \(i\)\(j\) (đầu vào)
  • \(w_{i j}\) là cân nặng mà người giao hàng phải mang qua cạnh \(ij\)
  • \(M\) là trọng lượng lớn nhất 1 người giao hàng có thể mang (đầu vào)
  • \(m, c, W\) được định nghĩa như ở trên (đầu vào)
  • Các thành phố được đánh số từ \(1\) đến \(n\) (mảng trong Julia bắt đầu từ \(1\) 😊 )

trong đó, coi \(x\)\(w\) là biến quyết định (decision variables) và \(m\) người giao hàng đang ở thành phố \(1\). Chúng ta có bài toán như sau:

Hàm tối ưu \(\min z = \sum_{i=1}^{n}\sum_{j=1}^{n}d_{i j} w_{ij} x_{ij}\)

Với điều kiện

  • \(\sum_{i=2}^n x_{1 i} = m\)\(m\) lần giao hàng đi ra khỏi thành phố \(1\)
  • \(\sum_{i=2}^n x_{i 1} = m\)\(m\) lần giao hàng đi vào lại thành phố \(1\)
  • \(x_{i i} = 0\qquad \forall i \) cái này hiển nhiên rồi
  • \(\sum_{j=1}^n x_{i j} = 1\qquad \forall i \in [2,n]\)\(1\) lần giao hàng đi ra khỏi thành phố \(i\)
  • \(\sum_{j=1}^n x_{j i} = 1\qquad \forall i \in [2,n]\)\(1\) lần giao hàng đi vào lại thành phố \(i\)
  • \(c_i = \sum x_{i j} w_{i j} - \sum x_{j i} w_{j i}\qquad \forall i \in [2,n]\) cân nặng khi đi ra bằng cân nặng khi đi vào - \(c_i\)
  • \((W + c_i) x_{i j} \le x_{i j} w_{i j} \le (W + M - c_j)x_{i j}\) cân nặng lúc nào cũng không dưới \(W\) và không trên \(W + M\)

Chúng ta dễ thấy công thức hiện tại không tuyến tính vì phụ thuộc vào tích của \(w_{ij} x_{ij}\) nhưng, may thay (cũng không hẳn là may mắn 😉), chỉ cần đặt \(y_{ij} = w_{ij} x_{ij}\) chúng ta sẽ có bài toán tối ưu tuyến tính (thường được giải nhanh và cách giải đơn giản hơn).

Hàm tối ưu \(\min z = \sum_{i=1}^{n}\sum_{j=1}^{n}d_{i j} y_{ij}\)

Với điều kiện

  • \(\sum_{i=2}^n x_{1 i} = m\)
  • \(\sum_{i=2}^n x_{i 1} = m\)
  • \(x_{i i} = 0\qquad \forall i \)
  • \(\sum_{j=1}^n x_{i j} = 1\qquad \forall i \in [2,n]\)
  • \(\sum_{j=1}^n x_{j i} = 1\qquad \forall i \in [2,n]\)
  • \(c_i = \sum y_{i j} - \sum y_{j i}\qquad \forall i \in [2,n]\)
  • \((W + c_i) x_{i j} \le y_{i j} \le (W + M - c_j)x_{i j}\)

Giải với Julia và SCIP

SCIP là 1 trình giải miễn phí các bài toán tối ưu có lẫn số nguyên khá nhanh và có interface dễ dùng cho Julia qua JuMP.

Sau khi cài đặt SCIPJulia.

Mở Julia console, gõ ] để vào phần pkg> rồi gõ add JuMP, add SCIP. Ta chỉ việc chuyển đổi công thức toán trên vào code

using JuMP, SCIP
... 
@variable(model, y[1:n, 1:n]) 
@variable(model, x[1:n, 1:n], Bin)

@objective(model, Min, sum(y[i, j] * d[i, j] for i = 1:n, j = range[1:n .!= i]))

@constraint(model, sum(x[1, i] for i = 2:n) == m)
@constraint(model, sum(x[i, 1] for i = 2:n) == m)

for i = 2:n
    @constraint(model, y[1, i] == w0 * x[1, i])
    @constraint(model, sum(x[i, j] for j = range[1:n .!= i]) == 1)
    @constraint(model, sum(x[j, i] for j = range[1:n .!= i]) == 1)
    @constraint(model, c[i] == sum(y[i, j] - y[j, i] for j = range[1:n .!= i]))
end

for i = 1:n
    for j = 1:n
        if i == j continue end
        @constraint(model, y[i, j] >= (w0 + c[i]) * x[i, j])
        @constraint(model, y[i, j] <= (w0 + cap - c[j]) * x[i, j])
    end
end
optimize!(model)

Chi tiết hơn github:quangio/tTSP (một số biến có tên khác)

Bài viết nháp gốc dạng pdf xanh sạch đẹp bằng tiếng anh tại đây

Note

Đừng bê code này vào ứng dụng của công ty bạn hoặc các thứ tương tự :v. SCIP và các solvers thông thường dùng branch-and-bound để giải mấy bài kiểu này. Theo quan sát của mình thì các thuật metaheuristic/local search/hybrid thường cho kết quả tốt nhanh hơn (VD: simulated annealing/tabu search + k-opt).