Requirement already satisfied: gurobipy in c:\dev\web\tlm\lectures\tlm_lecture_heuristics\intro_cop_heuristics_vrp\.venv\lib\site-packages (12.0.1)
Note: you may need to restart the kernel to use updated packages.Local search
Local search algorithms (a.k.a. greedy heuristics) aim at iteratively moving to higher quality solutions in the neighborhood until no further improvement can be made.
Neighborhood
A neighborhood function \(N\) is a mapping \(N : S \rightarrow 2^S\)[1] that assigns to each solution \(s \in S\) a set of neighbor solutions \(N(s) \subset S\).
Neighbor solutions \(s' \in N(s)\) are generated through a move operator \(m\) that performs a small pertubation to the solution \(s\).
[1] How many subsets does a set have?
Defining a neighborhood with the 2-opt operator
The neighborhood for the 2-opt operator (a.k.a. 2-exchange or subtour reversal operator) is represented by all the permutations obtained by removing and reconnecting all pairs of edges (except for the adjacent ones).
For example, for a solution (A, B, C, E, D, A) the following 2-opt moves can be performed:
| # | i | i+1 | j | j+1 | Neighbor |
|---|---|---|---|---|---|
| 1 | A | B | C | E | (A, C, B, E, D, A) |
| 2 | A | B | E | D | (A, E, C, B, D, A) |
| - | A | B | D | A | |
| 3 | B | C | E | D | (A, B, E, C, D, A) |
| 4 | B | C | D | A | (A, B, D, E, C, A) |
| 5 | C | E | D | A | (A, B, C, D, E, A) |
Notice that the exchange AB-DA (where i=A and j+1=A) is useless when considering a symmetric distance matrix: (A, B, C, E, D, A) and (D, E, C, B, A, D) have the same cost.
For a symmetric TSP with \(n\) nodes (i.e., \(n\) edges):
- Edge \(1\) can be exchanged with \(n-2\) edges (skipping adjacent) \(-1\) edge (useless exchange);
- Edge \(2\) can be exchanged with \(n-3\) edges;
- Edge \(k\) edge can be exchanged with \(n-k-1\) edges;
- Edge \(n-3\) can be exchanged with \(1\) edge.
Thus, the total number of exchanges is given by:
\[ (n-2) -1 + (n-3) + \dots + 1 = \frac{(n-1)(n-2)}{2} -1, \] which is equal to: \[ \frac{n(n-1)}{2}-n. \]
Hence, the neighborhood size for our example tour (A, B, C, E, D, A) comprised of \(n=5\) nodes is \(5\).
Example of 2-opt (sub-tour reversal) for the TSP (Talbi, 2009)
Iterative improvement algorithm
At iteration \(k\), the search moves from the current solution \(X_k\) to a new solution \(X_{k+1} \in N(X_k)\) (the neighborhood of \(X_k\)) only if the new point improves the value of the objective function \(F(X)\).
If no better \(X_{k+1}\) can be found in \(N(X_k)\) or if a user-specified number of iterations is reached, stop searching.
Local search behavior example (Talbi, 2009).
Immediate versus expanded neighborhood
Visiting the immediate neighborhood of \(X_k\) requires less local search computation at the expense of solution quality.
Whereas visiting the expanded neighborhood may lead to better solution quality at the expense of increased computation costs.
Study cases
Setting up the resources used throughout the document:
In the following, we present two study cases:
- Random city distribution: you can set the number of nodes that will be spread randomly in a [0,1]x[0,1] plane.
- 48 US capitals: a classic TSP instance comprised of the 48 US capitals with minimal tour length 33,523. You can find more classic TSP instances here.
Random city distribution
Set parameter Username
Set parameter LicenseID to value 2617855
Academic license - for non-commercial use only - expires 2026-02-04
Set parameter LazyConstraints to value 1
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
LazyConstraints 1
Optimize a model with 61 rows, 1830 columns and 3660 nonzeros
Model fingerprint: 0x7c2ee7ad
Variable types: 0 continuous, 1830 integer (1830 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e-02, 1e+00]
Bounds range [1e+00, 1e+00]
RHS range [2e+00, 2e+00]
Presolve time: 0.00s
Presolved: 61 rows, 1830 columns, 3660 nonzeros
Variable types: 0 continuous, 1830 integer (1830 binary)
Root relaxation: objective 5.942182e+00, 95 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 5.94218 0 6 - 5.94218 - - 0s
H 0 0 7.3932802 5.94218 19.6% - 0s
H 0 0 7.2963021 5.94218 18.6% - 0s
H 0 0 6.8939536 5.94218 13.8% - 0s
H 0 0 6.7775465 5.94218 12.3% - 0s
0 0 6.31651 0 14 6.77755 6.31651 6.80% - 0s
H 0 0 6.6922498 6.31651 5.61% - 0s
0 0 6.39880 0 - 6.69225 6.39880 4.38% - 0s
* 0 0 0 6.4067916 6.40679 0.00% - 0s
Cutting planes:
Zero half: 3
Lazy constraints: 18
Explored 1 nodes (133 simplex iterations) in 0.16 seconds (0.04 work units)
Thread count was 8 (of 8 available processors)
Solution count 6: 6.40679 6.69225 6.77755 ... 7.39328
Optimal solution found (tolerance 1.00e-04)
Best objective 6.406791568643e+00, best bound 6.406791568643e+00, gap 0.0000%
User-callback calls 208, time in user-callback 0.06 sec
2-opt implementation example
For the symmetric TSP a 2-opt local search only affects the tour length of the 4 cities involved in the move. Thus, by finding the greatest decrease in the tour length we can determine the best edges to swap.
Below, the complete 2-opt neighborhood for a random starting solution. The iterations where a improvement was made are marked in red.
Let us apply the 2-opt algorithm on better quality solution generated using the farthest addition (FA) construction heuristic:
Iterative improvement algorithm implementation example
We can use the 2-opt neighborhood to implement the iterative improvement algorithm:
Now, let us apply the iterative improvement algorithm on the solution constructed using the farthest addition (FA) construction heuristic:
48 US capitals
[ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47] 48
[(6734, 1453), (2233, 10), (5530, 1424), (401, 841), (3082, 1644), (7608, 4458), (7573, 3716), (7265, 1268), (6898, 1885), (1112, 2049), (5468, 2606), (5989, 2873), (4706, 2674), (4612, 2035), (6347, 2683), (6107, 669), (7611, 5184), (7462, 3590), (7732, 4723), (5900, 3561), (4483, 3369), (6101, 1110), (5199, 2182), (1633, 2809), (4307, 2322), (675, 1006), (7555, 4819), (7541, 3981), (3177, 756), (7352, 4506), (7545, 2801), (3245, 3305), (6426, 3173), (4608, 1198), (23, 2216), (7248, 3779), (7762, 4595), (7392, 2244), (3484, 2829), (6271, 2135), (4985, 140), (1916, 1569), (7280, 4899), (7509, 3239), (10, 2676), (6807, 2993), (5185, 3258), (3023, 1942)] 48
Set parameter LazyConstraints to value 1
Gurobi Optimizer version 12.0.1 build v12.0.1rc0 (win64 - Windows 11.0 (26100.2))
CPU model: 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
LazyConstraints 1
Optimize a model with 48 rows, 1128 columns and 2256 nonzeros
Model fingerprint: 0xdc564e0d
Variable types: 0 continuous, 1128 integer (1128 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e+02, 8e+03]
Bounds range [1e+00, 1e+00]
RHS range [2e+00, 2e+00]
Presolve time: 0.00s
Presolved: 48 rows, 1128 columns, 2256 nonzeros
Variable types: 0 continuous, 1128 integer (1128 binary)
Root relaxation: objective 3.166926e+04, 78 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 31669.2574 0 14 - 31669.2574 - - 0s
H 0 0 38335.094169 31669.2574 17.4% - 0s
H 0 0 37725.521880 31669.2574 16.1% - 0s
H 0 0 37029.139102 31669.2574 14.5% - 0s
H 0 0 36754.544644 31669.2574 13.8% - 0s
0 0 32320.8108 0 6 36754.5446 32320.8108 12.1% - 0s
H 0 0 35949.039512 32320.8108 10.1% - 0s
H 0 0 35606.532587 32320.8108 9.23% - 0s
H 0 0 35541.897685 32320.8108 9.06% - 0s
0 0 33021.6741 0 - 35541.8977 33021.6741 7.09% - 0s
0 0 33035.1856 0 - 35541.8977 33035.1856 7.05% - 0s
0 0 33043.1273 0 - 35541.8977 33043.1273 7.03% - 0s
0 0 33081.8621 0 - 35541.8977 33081.8621 6.92% - 0s
0 0 33255.2634 0 20 35541.8977 33255.2634 6.43% - 0s
0 0 33255.2634 0 14 35541.8977 33255.2634 6.43% - 0s
H 0 0 34966.443166 33255.2634 4.89% - 0s
0 0 33299.9040 0 12 34966.4432 33299.9040 4.77% - 0s
0 0 33299.9040 0 12 34966.4432 33299.9040 4.77% - 0s
H 0 0 34047.150934 33299.9040 2.19% - 0s
0 0 33434.9079 0 22 34047.1509 33434.9079 1.80% - 0s
0 0 33434.9079 0 12 34047.1509 33434.9079 1.80% - 0s
0 0 33484.2144 0 - 34047.1509 33484.2144 1.65% - 0s
0 0 33508.6532 0 33 34047.1509 33508.6532 1.58% - 0s
H 0 0 33706.845991 33508.6532 0.59% - 0s
0 0 33511.1608 0 19 33706.8460 33511.1608 0.58% - 0s
* 0 0 0 33523.708507 33523.7085 0.00% - 0s
Cutting planes:
Gomory: 6
MIR: 1
Zero half: 9
Lazy constraints: 2
Explored 1 nodes (442 simplex iterations) in 0.27 seconds (0.05 work units)
Thread count was 8 (of 8 available processors)
Solution count 10: 33523.7 33706.8 34047.2 ... 37725.5
Optimal solution found (tolerance 1.00e-04)
Best objective 3.352370850744e+04, best bound 3.352370850744e+04, gap 0.0000%
User-callback calls 292, time in user-callback 0.10 sec
Starting from a bad quality solution (e.g., random), many improvements can be made:
Limitation: Local optima entrapment
As shown in the examples above, at the end of the process, the search may be entrapped at local optima.
In order to escape entrapment, metaheuristics are designed to occasionally allow inferior moves, hoping to find better solutions by exploring neighborhoods that look unpromising (at first sight).