problemreductions/rules/
travelingsalesman_ilp.rs

1//! Reduction from TravelingSalesman to ILP (Integer Linear Programming).
2//!
3//! Uses position-based variables x_{v,k} with McCormick linearization.
4//! - Variables: x_{v,k} for vertex v at position k (binary), plus auxiliary y variables
5//! - Constraints: assignment, non-edge consecutive, McCormick
6//! - Objective: minimize total edge weight of the tour
7
8use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
9use crate::models::graph::TravelingSalesman;
10use crate::reduction;
11use crate::rules::traits::{ReduceTo, ReductionResult};
12use crate::topology::{Graph, SimpleGraph};
13
14/// Result of reducing TravelingSalesman to ILP.
15#[derive(Debug, Clone)]
16pub struct ReductionTSPToILP {
17    target: ILP<bool>,
18    /// Number of vertices in the source graph.
19    num_vertices: usize,
20    /// Edges of the source graph (for solution extraction).
21    source_edges: Vec<(usize, usize)>,
22}
23
24impl ReductionTSPToILP {
25    /// Variable index for x_{v,k}: vertex v at position k.
26    fn x_index(&self, v: usize, k: usize) -> usize {
27        v * self.num_vertices + k
28    }
29}
30
31impl ReductionResult for ReductionTSPToILP {
32    type Source = TravelingSalesman<SimpleGraph, i32>;
33    type Target = ILP<bool>;
34
35    fn target_problem(&self) -> &ILP<bool> {
36        &self.target
37    }
38
39    /// Extract solution: read tour permutation from x variables,
40    /// then map to edge selection for the source problem.
41    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
42        let n = self.num_vertices;
43
44        // Read tour: for each position k, find vertex v with x_{v,k} = 1
45        let mut tour = vec![0usize; n];
46        for k in 0..n {
47            for v in 0..n {
48                if target_solution[self.x_index(v, k)] == 1 {
49                    tour[k] = v;
50                    break;
51                }
52            }
53        }
54
55        // Map tour to edge selection
56        let mut edge_selection = vec![0usize; self.source_edges.len()];
57        for k in 0..n {
58            let u = tour[k];
59            let v = tour[(k + 1) % n];
60            // Find the edge index for (u, v) or (v, u)
61            for (idx, &(a, b)) in self.source_edges.iter().enumerate() {
62                if (a == u && b == v) || (a == v && b == u) {
63                    edge_selection[idx] = 1;
64                    break;
65                }
66            }
67        }
68
69        edge_selection
70    }
71}
72
73#[reduction(
74    overhead = {
75        num_vars = "num_vertices^2 + 2 * num_vertices * num_edges",
76        num_constraints = "num_vertices^3 + -1 * num_vertices^2 + 2 * num_vertices + 4 * num_vertices * num_edges",
77    }
78)]
79impl ReduceTo<ILP<bool>> for TravelingSalesman<SimpleGraph, i32> {
80    type Result = ReductionTSPToILP;
81
82    fn reduce_to(&self) -> Self::Result {
83        let n = self.graph().num_vertices();
84        let graph = self.graph();
85        let edges_with_weights = self.edges();
86        let source_edges: Vec<(usize, usize)> =
87            edges_with_weights.iter().map(|&(u, v, _)| (u, v)).collect();
88        let edge_weights: Vec<f64> = edges_with_weights
89            .iter()
90            .map(|&(_, _, w)| w as f64)
91            .collect();
92        let m = source_edges.len();
93
94        // Variable layout:
95        // [0, n²): x_{v,k} = vertex v at position k
96        // [n², n² + 2mn): auxiliary y variables for McCormick linearization
97        //   For edge_idx e and position k:
98        //     y_{forward} at n² + e * 2n + 2k      (x_{u,k} * x_{v,(k+1)%n})
99        //     y_{reverse} at n² + e * 2n + 2k + 1  (x_{v,k} * x_{u,(k+1)%n})
100        let num_x = n * n;
101        let num_y = 2 * m * n;
102        let num_vars = num_x + num_y;
103
104        let x_idx = |v: usize, k: usize| -> usize { v * n + k };
105        let y_idx =
106            |edge: usize, k: usize, dir: usize| -> usize { num_x + edge * 2 * n + 2 * k + dir };
107
108        let mut constraints = Vec::new();
109
110        // Constraint 1: Each vertex has exactly one position
111        for v in 0..n {
112            let terms: Vec<(usize, f64)> = (0..n).map(|k| (x_idx(v, k), 1.0)).collect();
113            constraints.push(LinearConstraint::eq(terms, 1.0));
114        }
115
116        // Constraint 2: Each position has exactly one vertex
117        for k in 0..n {
118            let terms: Vec<(usize, f64)> = (0..n).map(|v| (x_idx(v, k), 1.0)).collect();
119            constraints.push(LinearConstraint::eq(terms, 1.0));
120        }
121
122        // Constraint 3: Non-edge consecutive prohibition
123        // For each ordered pair (v, w) where {v, w} ∉ E and v ≠ w:
124        //   x_{v,k} + x_{w,(k+1) mod n} <= 1 for all k
125        for v in 0..n {
126            for w in 0..n {
127                if v == w {
128                    continue;
129                }
130                if graph.has_edge(v, w) {
131                    continue;
132                }
133                for k in 0..n {
134                    constraints.push(LinearConstraint::le(
135                        vec![(x_idx(v, k), 1.0), (x_idx(w, (k + 1) % n), 1.0)],
136                        1.0,
137                    ));
138                }
139            }
140        }
141
142        // Constraint 4: McCormick linearization for auxiliary variables
143        // For each edge (u, v) at index e:
144        //   Forward (dir=0): y = x_{u,k} * x_{v,(k+1)%n}
145        //   Reverse (dir=1): y = x_{v,k} * x_{u,(k+1)%n}
146        for (e, &(u, v)) in source_edges.iter().enumerate() {
147            for k in 0..n {
148                let k_next = (k + 1) % n;
149
150                // Forward: y_{e,k,0} = x_{u,k} * x_{v,k_next}
151                let y_fwd = y_idx(e, k, 0);
152                let xu = x_idx(u, k);
153                let xv_next = x_idx(v, k_next);
154                constraints.push(LinearConstraint::le(vec![(y_fwd, 1.0), (xu, -1.0)], 0.0));
155                constraints.push(LinearConstraint::le(
156                    vec![(y_fwd, 1.0), (xv_next, -1.0)],
157                    0.0,
158                ));
159                constraints.push(LinearConstraint::ge(
160                    vec![(y_fwd, 1.0), (xu, -1.0), (xv_next, -1.0)],
161                    -1.0,
162                ));
163
164                // Reverse: y_{e,k,1} = x_{v,k} * x_{u,k_next}
165                let y_rev = y_idx(e, k, 1);
166                let xv = x_idx(v, k);
167                let xu_next = x_idx(u, k_next);
168                constraints.push(LinearConstraint::le(vec![(y_rev, 1.0), (xv, -1.0)], 0.0));
169                constraints.push(LinearConstraint::le(
170                    vec![(y_rev, 1.0), (xu_next, -1.0)],
171                    0.0,
172                ));
173                constraints.push(LinearConstraint::ge(
174                    vec![(y_rev, 1.0), (xv, -1.0), (xu_next, -1.0)],
175                    -1.0,
176                ));
177            }
178        }
179
180        // Objective: minimize Σ_{e=(u,v)} w_e * Σ_k (y_{e,k,0} + y_{e,k,1})
181        let mut objective: Vec<(usize, f64)> = Vec::new();
182        for (e, &w) in edge_weights.iter().enumerate() {
183            for k in 0..n {
184                objective.push((y_idx(e, k, 0), w));
185                objective.push((y_idx(e, k, 1), w));
186            }
187        }
188
189        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize);
190
191        ReductionTSPToILP {
192            target,
193            num_vertices: n,
194            source_edges,
195        }
196    }
197}
198
199#[cfg(test)]
200#[path = "../unit_tests/rules/travelingsalesman_ilp.rs"]
201mod tests;
problemreductions/rules/travelingsalesman_ilp.rs

problemreductions/rules/
travelingsalesman_ilp.rs
