Cấu trúc dữ liệu và Giải thuật Bài 143: Thuật toán Network Simplex (Network Simplex Algorithm) - Tối ưu hóa luồng chi phí cực tiểu trên cây cơ sở với hiệu năng thực tế vượt trội.

Trong quy hoạch tuyến tính (Linear Programming - LP), thuật toán Đơn hình (Simplex Algorithm) của George Dantzig (1947) là một phương pháp nền tảng. Tuy nhiên, khi áp dụng thuật toán Đơn hình tổng quát lên bài toán Luồng chi phí cực tiểu (Min-Cost Flow), mỗi bước xoay (pivot) đòi hỏi nghịch đảo ma trận hoặc cập nhật phân rã ma trận cơ sở tốn $O(V^3)$ phép tính số học.

Năm 1969, George Dantzig cùng các cộng sự nhận thấy: Ma trận ràng buộc của bài toán mạng luồng có tính chất hoàn toàn đơn mô (Totally Unimodular). Mọi cơ sở hợp lệ (Basis) của bài toán quy hoạch tuyến tính mạng đều tương ứng song ánh 1-1 với một Cây khung cơ sở (Spanning Tree Basis) trên đồ thị. Phép nghịch đảo ma trận phức tạp biến đổi thành các thao tác duyệt cây, lần vết chu trình và cập nhật thế năng trên cây với chi phí trung bình cực nhỏ.

Đó chính là sự ra đời của Network Simplex (Đơn hình Mạng). Dù độ phức tạp lý thuyết trong trường hợp xấu nhất có thể không đa thức với một số luật chọn bước xoay ngây thơ, nhưng trong thực tế sản xuất và thi đấu, Network Simplex chạy nhanh hơn các thuật toán đa thức (như Cost-Scaling hay Push-Relabel) từ 5 đến 50 lần. Đây là thuật toán mặc định trong các bộ giải tối ưu hóa hàng đầu thế giới (như CPLEX, Gurobi, LEMON, GLPK).

1. Cấu trúc Cây cơ sở (Spanning Tree Basis)

Xét bài toán Luồng chi phí cực tiểu có nhu cầu tại đỉnh (Transshipment / Min-Cost Flow):

$$\min \sum_{(u, v) \in E} w(u, v) \cdot f(u, v)$$

thỏa mãn:

$$\sum_{v: (u, v) \in E} f(u, v) - \sum_{v: (v, u) \in E} f(v, u) = b(u) \quad (\forall u \in V)$$

$$0 \le f(u, v) \le c(u, v) \quad (\forall (u, v) \in E)$$

(với $b(u) > 0$ là đỉnh phát, $b(u) < 0$ là đỉnh thu, thỏa mãn $\sum b(u) = 0$).

Phân hoạch Cây cơ sở $(T, L, U)$

Một phương án cơ sở khả thi (Basic Feasible Solution) của Network Simplex phân hoạch tập cạnh $E$ thành 3 tập rời nhau:

1. Tập cạnh cơ sở $T$ (Tree Arcs): Tạo thành một cây khung gồm đúng $\vert{}V\vert{} - 1$ cạnh nối tất cả các đỉnh của đồ thị. Luồng trên các cạnh cơ sở có thể nhận giá trị tự do trong đoạn $[0, c(e)]$.

2. Tập phi cơ sở cận dưới $L$ (Lower-bound Non-tree Arcs): Các cạnh ngoài cây có luồng bằng đúng cận dưới: $f(e) = 0$.

3. Tập phi cơ sở cận trên $U$ (Upper-bound Non-tree Arcs): Các cạnh ngoài cây có luồng bão hòa bằng đúng cận trên: $f(e) = c(e)$.

                     (Tập cạnh cơ sở T - Cây khung)
                                  (1)
                                /     \
                              (2)     (3)
                              /         \
                            (4)         (5)

         Cạnh ngoài (L): f(e) = 0  --- - - - - ---> Cạnh ngoài (U): f(e) = c(e)

2. Thế năng đỉnh và Điều kiện tối ưu (Optimality Conditions)

Mỗi đỉnh $u \in V$ được gán một biến đối ngẫu gọi là thế năng đỉnh (Node Potential) $\pi(u)$.

Chi phí rút gọn (Reduced Cost)

Với mỗi cạnh $e = (u, v)$, chi phí rút gọn $w^\pi(u, v)$ được định nghĩa:

$$w^\pi(u, v) = w(u, v) - \pi(u) + \pi(v)$$

Tính toán Thế năng trên Cây cơ sở

Theo điều kiện bù (Complementary Slackness), mọi cạnh thuộc cây cơ sở $T$ bắt buộc phải có chi phí rút gọn bằng $0$:

$$\forall (u, v) \in T: \quad w^\pi(u, v) = 0 \iff \pi(v) = \pi(u) - w(u, v)$$

Gán thế năng cho một đỉnh gốc tùy ý $\pi(\text{root}) = 0$. Khi đó, toàn bộ thế năng $\pi(u)$ của tất cả các đỉnh còn lại được xác định duy nhất thông qua một lượt duyệt DFS/BFS trên cây $T$ trong thời gian $O(V)$.

Tiêu chuẩn tối ưu toàn cục

Phương án cơ sở hiện tại đạt tối ưu khi và chỉ khi mọi cạnh ngoài cây đều thỏa mãn:

1. $\forall e \in L: \quad w^\pi(e) \ge 0$ (Không thể giảm chi phí bằng cách tăng luồng).

2. $\forall e \in U: \quad w^\pi(e) \le 0$ (Không thể giảm chi phí bằng cách giảm luồng).

3. Bước xoay Đơn hình Mạng (Network Simplex Pivot)

Nếu tồn tại một cạnh vi phạm tiêu chuẩn tối ưu, ta thực hiện một bước xoay cơ sở (Pivot Step) gồm 4 pha:

Pha 1: Chọn Cạnh vào cơ sở (Entering Arc - $e_{\text{in}}$)

Tìm một cạnh $e_{\text{in}} = (u, v) \notin T$ vi phạm điều kiện tối ưu:

• $e_{\text{in}} \in L$ nhưng có $w^\pi(e_{\text{in}}) < 0$ (Cần đẩy thêm luồng theo chiều xuôi).

• $e_{\text{in}} \in U$ nhưng có $w^\pi(e_{\text{in}}) > 0$ (Cần rút bớt luồng theo chiều ngược).

Pha 2: Xác định Chu trình cơ sở (Fundamental Cycle)

Thêm cạnh $e_{\text{in}}$ vào cây $T$ sẽ tạo ra đúng một chu trình đơn duy nhất $C$. Chu trình này được tìm bằng cách lần ngược từ $u$ và $v$ về tổ tiên chung gần nhất (LCA) trên cây $T$.

Pha 3: Tăng/giảm luồng quanh chu trình và Chọn Cạnh ra (Leaving Arc - $e_{\text{out}}$)

Định hướng chu trình $C$ theo chiều tăng luồng của $e_{\text{in}}$:

• Với các cạnh trên $C$ có cùng chiều duyệt: dung lượng còn lại có thể đẩy thêm là $c(e) - f(e)$.

• Với các cạnh trên $C$ ngược chiều duyệt: lượng luồng có thể rút bớt là $f(e)$.

Lượng luồng điều chỉnh tối đa quanh chu trình là điểm nghẽn (bottleneck):

$$\theta = \min_{e \in C} \{ \text{khả năng điều chỉnh của } e \}$$

Đẩy lượng luồng $\theta$ quanh chu trình $C$. Cạnh đạt tới giới hạn (chạm $0$ hoặc chạm dung lượng $c$) được chọn làm Cạnh ra khỏi cơ sở ($e_{\text{out}}$). Cạnh này sẽ rời khỏi cây $T$ để gia nhập vào tập $L$ hoặc $U$.

Pha 4: Cập nhật Cây và Cập nhật Thế năng

1. Cắt bỏ cạnh $e_{\text{out}}$ khỏi cây $T$ và thêm cạnh $e_{\text{in}}$ vào cây $T$, hình thành cây khung cơ sở mới.

2. Việc thay đổi cây làm thay đổi thế năng của toàn bộ cây con bị ngắt kết nối. Thế năng của các đỉnh trong cây con này được cộng hoặc trừ một lượng cố định $\Delta = w^\pi(e_{\text{in}})$ trong thời gian $O(\text{size của subtree})$.

           Thêm e_in vào T:                       Cắt bỏ e_out:
                (LCA)                                  (LCA)
               /     \                                /     \
             (a)     (b)                            (a)     (b)
             /         \                            /
           (u) ======= (v)                        (u) ======= (v)
             \--e_in--/                             \--e_in--/
        (Tạo chu trình đơn C)                   (Cây khung mới hoàn chỉnh)

4. Khởi tạo Cơ sở ban đầu bằng Phương pháp Đỉnh nhân tạo (Big-M Method)

Để bắt đầu thuật toán Network Simplex, ta cần một cây khung cơ sở khả thi ban đầu. Phương pháp chuẩn mực nhất là thêm một Đỉnh siêu gốc nhân tạo (Artificial Root) $R$:

1. Với mỗi đỉnh $i \in V$:

   • Nếu $b(i) \ge 0$ (Đỉnh phát): Thêm cạnh $(R \to i)$ với dung lượng $\infty$, chi phí lớn $M$, và gán luồng ban đầu $f(R \to i) = b(i)$.

   • Nếu $b(i) < 0$ (Đỉnh thu): Thêm cạnh $(i \to R)$ với dung lượng $\infty$, chi phí lớn $M$, và gán luồng ban đầu $f(i \to R) = -b(i)$.

2. Tập hợp $\vert{}V\vert{}$ cạnh nối với đỉnh nhân tạo $R$ tạo thành một cây khung cơ sở ban đầu hoàn chỉnh với gốc tại $R$.

3. Mọi cạnh gốc của bài toán được đưa vào tập $L$ với $f(e) = 0$.

4. Chạy các bước xoay của Network Simplex. Vì chi phí $M$ rất lớn, thuật toán sẽ tự động đẩy toàn bộ luồng ra khỏi các cạnh nhân tạo và loại chúng khỏi cơ sở để thu được phương án tối ưu của đồ thị gốc.

5. Hiện thực hóa Network Simplex bằng Golang

Chương trình dưới đây cài đặt hoàn chỉnh cấu trúc dữ liệu Network Simplex, lưu trữ cây cơ sở qua các mảng quan hệ cha-con (parent, predEdge), và luật xoay ứng viên (First-Eligible Pivot Rule):

Go

package main

import (
	"fmt"
	"math"
)

const INF = math.MaxInt64 / 4
const BIG_M = 1_000_000_000 // Chi phí lớn cho các cạnh nhân tạo

type Edge struct {
	id       int
	u, v     int
	cap      int64
	cost     int64
	flow     int64
	inTree   bool
	isUpper  bool // True nếu e thuộc tập U (f = cap)
}

type NetworkSimplex struct {
	n        int
	edges    []*Edge
	supplies []int64

	// Cấu trúc cây cơ sở
	parent   []int // parent[u] trên cây cơ sở
	predEdge []*Edge
	depth    []int
	pi       []int64 // Thế năng đỉnh
	root     int
}

func NewNetworkSimplex(n int, supplies []int64) *NetworkSimplex {
	return &NetworkSimplex{
		n:        n,
		supplies: supplies,
		parent:   make([]int, n+1),
		predEdge: make([]*Edge, n+1),
		depth:    make([]int, n+1),
		pi:       make([]int64, n+1),
		root:     n, // Đỉnh nhân tạo làm gốc
	}
}

func (ns *NetworkSimplex) AddEdge(u, v int, cap, cost int64) {
	e := &Edge{
		id:      len(ns.edges),
		u:       u,
		v:       v,
		cap:     cap,
		cost:    cost,
		flow:    0,
		inTree:  false,
		isUpper: false,
	}
	ns.edges = append(ns.edges, e)
}

// Khởi tạo cây cơ sở ban đầu bằng Big-M method
func (ns *NetworkSimplex) initTree() {
	artRoot := ns.root

	for i := 0; i < ns.n; i++ {
		b := ns.supplies[i]
		if b >= 0 {
			// Cung nhân tạo từ artRoot -> i
			e := &Edge{
				id:     len(ns.edges),
				u:      artRoot,
				v:      i,
				cap:    INF,
				cost:   BIG_M,
				flow:   b,
				inTree: true,
			}
			ns.edges = append(ns.edges, e)
			ns.parent[i] = artRoot
			ns.predEdge[i] = e
			ns.depth[i] = 1
		} else {
			// Cung nhân tạo từ i -> artRoot
			e := &Edge{
				id:     len(ns.edges),
				u:      i,
				v:      artRoot,
				cap:    INF,
				cost:   BIG_M,
				flow:   -b,
				inTree: true,
			}
			ns.edges = append(ns.edges, e)
			ns.parent[i] = artRoot
			ns.predEdge[i] = e
			ns.depth[i] = 1
		}
	}

	ns.parent[artRoot] = -1
	ns.depth[artRoot] = 0
	ns.updatePotentials()
}

// Cập nhật thế năng đỉnh trên toàn cây cơ sở
func (ns *NetworkSimplex) updatePotentials() {
	// Dựng danh sách kề cho cây cơ sở
	treeAdj := make([][]int, ns.n+1)
	for i := 0; i < ns.n; i++ {
		p := ns.parent[i]
		if p != -1 {
			treeAdj[p] = append(treeAdj[p], i)
		}
	}

	ns.pi[ns.root] = 0
	queue := []int{ns.root}

	for len(queue) > 0 {
		u := queue[0]
		queue = queue[1:]

		for _, v := range treeAdj[u] {
			e := ns.predEdge[v]
			if e.u == u {
				// Cạnh hướng u -> v: w_pi = w - pi[u] + pi[v] = 0 => pi[v] = pi[u] - w
				ns.pi[v] = ns.pi[u] - e.cost
			} else {
				// Cạnh hướng v -> u: w_pi = w - pi[v] + pi[u] = 0 => pi[v] = pi[u] + w
				ns.pi[v] = ns.pi[u] + e.cost
			}
			ns.depth[v] = ns.depth[u] + 1
			queue = append(queue, v)
		}
	}
}

func (ns *NetworkSimplex) reducedCost(e *Edge) int64 {
	return e.cost - ns.pi[e.u] + ns.pi[e.v]
}

// Tìm LCA của hai đỉnh trên cây
func (ns *NetworkSimplex) findLCA(u, v int) int {
	for u != v {
		if ns.depth[u] > ns.depth[v] {
			u = ns.parent[u]
		} else {
			v = ns.parent[v]
		}
	}
	return u
}

// Thực hiện một bước xoay cơ sở
func (ns *NetworkSimplex) pivot(eIn *Edge) {
	u, v := eIn.u, eIn.v
	lca := ns.findLCA(u, v)

	// Xác định chiều tăng luồng trên chu trình tạo bởi eIn
	// Nếu eIn thuộc L (f=0): đẩy luồng theo chiều u -> v
	// Nếu eIn thuộc U (f=cap): đẩy luồng theo chiều v -> u
	forward := !eIn.isUpper

	type CycleEdge struct {
		edge      *Edge
		isForward bool
	}

	var cycle []CycleEdge

	// Lần vết từ u lên LCA
	curr := u
	for curr != lca {
		e := ns.predEdge[curr]
		isFwd := (e.v == curr) != forward
		cycle = append(cycle, CycleEdge{edge: e, isForward: isFwd})
		curr = ns.parent[curr]
	}

	// Lần vết từ v lên LCA
	curr = v
	for curr != lca {
		e := ns.predEdge[curr]
		isFwd := (e.u == curr) != forward
		cycle = append(cycle, CycleEdge{edge: e, isForward: isFwd})
		curr = ns.parent[curr]
	}

	// Xác định điểm nghẽn theta và Cạnh ra khỏi cơ sở (Leaving Arc)
	var theta int64 = INF
	if forward {
		theta = eIn.cap - eIn.flow
	} else {
		theta = eIn.flow
	}

	var eOut *Edge = eIn
	outHitsLower := !forward

	for _, ce := range cycle {
		e := ce.edge
		if ce.isForward {
			resCap := e.cap - e.flow
			if resCap < theta {
				theta = resCap
				eOut = e
				outHitsLower = false
			}
		} else {
			if e.flow < theta {
				theta = e.flow
				eOut = e
				outHitsLower = true
			}
		}
	}

	// Cập nhật luồng trên chu trình
	if forward {
		eIn.flow += theta
	} else {
		eIn.flow -= theta
	}

	for _, ce := range cycle {
		e := ce.edge
		if ce.isForward {
			e.flow += theta
		} else {
			e.flow -= theta
		}
	}

	// Nếu cạnh vào chính là cạnh ra (chỉ đổi từ L sang U hoặc ngược lại)
	if eIn == eOut {
		eIn.isUpper = (eIn.flow == eIn.cap)
		return
	}

	// Cập nhật cấu trúc cây: eIn vào cây, eOut rời cây
	eIn.inTree = true
	eOut.inTree = false
	eOut.isUpper = !outHitsLower

	// Tái cấu trúc liên kết parent và predEdge dọc theo đường đi
	var joinNode, splitNode int
	if ns.depth[eIn.u] > ns.depth[eIn.v] {
		joinNode = eIn.u
		splitNode = eIn.v
	} else {
		joinNode = eIn.v
		splitNode = eIn.u
	}

	// Đảo chiều cây con từ joinNode về eOut
	curr = joinNode
	pEdge := eIn
	pNode := splitNode

	for {
		nextCurr := ns.parent[curr]
		nextEdge := ns.predEdge[curr]

		ns.parent[curr] = pNode
		ns.predEdge[curr] = pEdge

		if curr == eOut.u || curr == eOut.v {
			break
		}

		pNode = curr
		pEdge = nextEdge
		curr = nextCurr
	}

	ns.updatePotentials()
}

// Solve giải bài toán Min-Cost Flow
func (ns *NetworkSimplex) Solve() (bool, int64) {
	ns.initTree()

	// Vòng lặp xoay Đơn hình
	for {
		var entering *Edge = nil

		// Chiến lược chọn cạnh vi phạm (First-Eligible Pivot Rule)
		for _, e := range ns.edges {
			if e.inTree {
				continue
			}

			rc := ns.reducedCost(e)
			if !e.isUpper && rc < 0 {
				entering = e
				break
			}
			if e.isUpper && rc > 0 {
				entering = e
				break
			}
		}

		if entering == nil {
			break // Không còn cạnh nào vi phạm -> Đạt tối ưu toàn cục!
		}

		ns.pivot(entering)
	}

	// Kiểm tra xem có luồng nào còn tồn đọng trên các cạnh nhân tạo không
	for _, e := range ns.edges {
		if (e.u == ns.root || e.v == ns.root) && e.flow > 0 {
			return false, 0 // Không thể thỏa mãn bài toán luồng (Infeasible)
		}
	}

	var totalCost int64 = 0
	for _, e := range ns.edges {
		if e.u != ns.root && e.v != ns.root {
			totalCost += e.flow * e.cost
		}
	}

	return true, totalCost
}

func main() {
	// Mạng 4 đỉnh:
	// Đỉnh 0: Nguồn phát (+4)
	// Đỉnh 3: Đích tiêu thụ (-4)
	// Đỉnh 1, 2: Trung chuyển (0)
	supplies := []int64{4, 0, 0, -4}
	ns := NewNetworkSimplex(4, supplies)

	// Các cung: (u, v, cap, cost)
	ns.AddEdge(0, 1, 3, 2)
	ns.AddEdge(0, 2, 3, 5)
	ns.AddEdge(1, 2, 2, 1)
	ns.AddEdge(1, 3, 4, 4)
	ns.AddEdge(2, 3, 3, 1)

	feasible, minCost := ns.Solve()
	if !feasible {
		fmt.Println("Bài toán vô nghiệm!")
		return
	}

	fmt.Println("=== KẾT QUẢ THUẬT TOÁN NETWORK SIMPLEX ===")
	fmt.Printf("Chi phí tối ưu toàn mạng: %d\n\n", minCost)
	fmt.Println("Phân bổ luồng trên các cạnh gốc:")
	for _, e := range ns.edges {
		if e.u != ns.root && e.v != ns.root && e.flow > 0 {
			fmt.Printf("- Cung %d -> %d: Luồng = %d/%d | Đơn giá = %d | Chi phí = %d\n",
				e.u, e.v, e.flow, e.cap, e.cost, e.flow*e.cost)
		}
	}
}

6. Chiến lược chọn Cạnh vào (Pivot Rules) và Kỹ thuật Chống lặp chu trình (Anti-Cycling)

Trong thực tế, tốc độ hội tụ của Network Simplex phụ thuộc tới $80\%$ vào Luật chọn Cạnh vào cơ sở (Pivot Rule):

Suy thoái (Degeneracy) và Luật Cơ sở mạnh của Cunningham (Strong Basis Rule)

Một bước xoay được gọi là suy thoái (degenerate) nếu lượng luồng tăng quanh chu trình $\theta = 0$. Khi xảy ra suy thoái liên tiếp, thuật toán có thể rơi vào hiện tượng lặp vòng vô hạn (cycling).

• Quy tắc Cunningham (1976): Cây cơ sở luôn được duy trì sao cho mọi cạnh có luồng bằng $0$ đều hướng về phía gốc, và mọi cạnh bão hòa ($f = c$) đều hướng ra xa gốc. Quy tắc này đảm bảo rằng hiện tượng suy thoái không bao giờ tạo thành chu trình lặp, chứng minh thuật toán luôn kết thúc sau một số hữu hạn bước.

7. So sánh Toàn diện Hiệu năng Thực tế

8. Ứng dụng Thực tế

• Bài toán Khoảng cách người xúc đất (Earth Mover's Distance - EMD / Wasserstein Metric): Trong thị giác máy tính và học máy, tính toán khoảng cách phân phối xác suất giữa hai bức ảnh hoặc đám mây điểm 3D được mô hình hóa thành bài toán vận tải cỡ lớn ($N \approx 10^4$). Network Simplex là thuật toán cốt lõi tính toán EMD nhanh nhất hiện nay.

• Tối ưu hóa Lưới điện thông minh & Điều độ Kinh tế (Optimal Power Flow Dispatch): Giải bài toán dòng công suất mạng lưới truyền tải điện quốc gia có tổn thất đường dây theo thời gian thực.

• Khớp lệnh thị trường giao dịch Logistics (Ride-Hailing & Freight Matching Engines): Hệ thống điều phối hàng chục nghìn tài xế và chuyến hàng theo từng khung giờ cố định.