## 2 Materials and Methods
### Kernel Search
Linking two routes, one with the arc \((i,0)\) and another with the arc \((0,j)\), results in the cost of  
\[s_{ij}=c_{i0}+c_{ij}-c_{ij}. \tag{9}\]  
The CWS algorithm iteratively merges pairs of routes maximizing cost savings \(s_{ij}\) subject to the merged route feasibility (respecting the capacity restriction). There are two versions of CWS: sequential and parallel. In the sequential version, routes are expanded one at a time, until they are no longer feasible, whereas in the parallel version, several routes are built at the same time. According to [2], the parallel version usually offers better results than the sequential one. Due to its greedy nature, the CWS heuristics is not the most accurate. In order to achieve solutions with better qualities, the savings Formula (9) has been modified by adding new terms and parameters to it, as follows:  
\[s_{ij}=c_{0j}+c_{i0}-\lambda c_{ij}+\mu|c_{0i}+c_{0}|+\
u\frac{d_{i}+d_{j}}{\bar{ d}}. \tag{10}\]  
Here, \(\lambda\) is a route shape parameter [27; 28] that controls the relative significance of the distance between two customers, \(\mu\) is a weight parameter [29] that takes into account the asymmetry between two merged customers in terms of their distances from the depot, and \(\
u\) is a customer demand parameter [25] that prefers customers with larger demands relative to the average demand \(\bar{d}\), i.e., \(\bar{d}=\frac{\sum_{i=1}^{n}d_{i}}{n}\). To adjust the specific values of these three independent parameters, the authors evaluated 8820 different parameter vectors (\(\lambda,\mu,\
u\)) by varying \(\lambda\in\langle 0.1,2\rangle\), and \(\mu,\
u\in\langle 0,2\rangle\), respectively, with a step size of 0.1. Later, a genetic algorithm [30], and empirically adjusted greedy heuristics [31] were used as the parameter tuning procedures, which required much shorter computing times. In our implementation, the Taguchi method [32] was applied to set the parameter values for the parameterized CWS algorithm. The initialization phase of KS-CVRP works as follows: a parameterized CWS method is used to find the initial feasible solution of CVRP, and its value is saved as the initial upper bound \(z^{*}\). All the variables \(x_{ij}\) corresponding to the arcs \((i,j)\) between customers \(i\) and \(j\) that belong to the found CWS solution are taken to create the initial kernel \(K\). Variables \(x_{ij}\), \(i,j\in V\), \(i\
eq j\) not present in the kernel are ranked according to a suitable criterion reflecting the variable likelihood to be included in the optimal integer solution. Subsequently, a sequence of buckets with these such ordered variables is created. Buckets \(\{B_{i}\}_{i=1}^{N_{b}}\) are iteratively explored during the _improvement phase_. At each iteration a modified MILP sub-problem is solved on the union of the current kernel and the processed bucket \(B_{i}\). Reminder: we are only seeking such a solution that is better than the current best one, and that uses at least one new item from the current bucket. Figure 3 illustrates the evolution of solution improvement of problem instance No. 8 with 40 nodes, two vehicles, and 11 buckets. Because the original KS (KS-O) did not find an initial feasible  
Figure 2: An optimal solution of the LP relaxation of CVRP (the depot is red). solution, the initial value is infinity (in the picture represented by value of 618). This method reached the first feasible solution only in the third iteration. KS-CVRP starts from the initial solution with a value of 464. While KS-O achieved the optimal solution (value of 458) after exploring all buckets, KS-CVRP was enough to explore the first bucket. Regarding the latter constraint, in case of large-sized instances, a time limit must be set to solve each restricted problem. Therefore, some of these problems may not be solved optimally within the specified time. To facilitate computing, we devised an iterative improvement procedure that utilizes the strategy of the Exact Iterative Method for CVRP (EIM-CVRP) [33].