## 2 Materials and Methods
### MILP Formulation  
Let \(N=\{1,2,\ldots,n\}\) be the set of all customers and \(N_{0}=N\cup\{0\}\) be the set of all customers and the depot. For each customer \(j\in N\), \(d_{j}\) and \(p_{j}\) denote, respectively, its delivery and the pickup amounts. The fleet of vehicles is denoted by \(V=\{1,2,\ldots,m\}\). Note that for the sake of modeling simplicity, the number of vehicles available of each type is already given in this formulation. This assumption, however, limits the decision on the fleet composition, and commonly, a sufficiently high number of vehicles has to be assumed for each vehicle type such that optimality can be guaranteed. For each vehicle \(k\in V\), \(q^{k}\), \(f^{k}\), and \(g^{k}\) denote, respectively, its capacity, fixed cost, and variable cost. In addition, for each pair of nodes \(i,j\in N_{0}\), \(c_{ij}\) represents the distance between the nodes.  
We aim to determine the fleet composition and vehicle routes that minimize the total routing cost. Let \(x_{ij}^{k}\) be the binary decision variable indicating whether vehicle \(k\) traverses are _ij_ or not. The continuous variable \(y_{ij}^{k}\) denotes the total load picked up by vehicle \(k\) while traversing arc _ij_. Conversely, the continuous variable \(z_{ij}^{k}\) denotes the total load to be delivered by vehicle \(k\) while traversing arc _ij_. The optimization problem can be formulated as follows.  
\[\min\sum_{k\in V}\sum_{j\in N}f^{k}x_{0j}^{k}+\sum_{k\in V}\sum_{i \in N_{0}}\sum_{j\in N_{0}}c_{ij}g^{k}x_{ij}^{k}\] (1) s.t. \[\sum_{k\in V}\sum_{i\in N_{0}}x_{ij}^{k}=1 \forall j\in N \tag{2}\] \[\sum_{i\in N_{0}}x_{iv}^{k}-\sum_{j\in N_{0}}x_{vj}^{k}=0 \forall v\in N_{0},k\in V\] (3) \[\sum_{j\in N}x_{0j}^{k}\leq 1 \forall k\in V\] (4) \[y_{0j}^{k}=0 \forall j\in N,k\in V\] (5) \[z_{0}^{k}=0 \forall i\in N,k\in V\] (6) \[\sum_{j\in N_{0}}\sum_{k\in V}y_{vj}^{k}-\sum_{i\in N_{0}}\sum_{k \in V}y_{iv}^{k}=p_{v} \forall v\in N\] (7) \[\sum_{i\in N_{0}}\sum_{k\in V}z_{iv}^{k}-\sum_{j\in N_{0}}\sum_{k \in V}z_{vj}^{k}=d_{v} \forall v\in N\] (8) \[\sum_{i\in N}\sum_{k\in V}y_{i0}^{k}=\sum_{v\in N}p_{v}\] (9) \[\sum_{j\in N}\sum_{k\in V}z_{0j}^{k}=\sum_{v\in N}d_{v}\] (10) \[y_{ij}^{k}+z_{ij}^{k}\leq d^{k}x_{ij}^{k} \forall i\in N_{0},j\in N_{0},k\in V\] (11) \[x_{ij}^{k}\in\{0,1\},y_{ij}^{k}\geq 0,z_{ij}^{k}\geq 0 \forall i\in N_{0},j\in N_{0},k\in V \tag{12}\]  
The objective function (1) minimizes the total fixed costs of vehicles and variable transportation costs. Constraints (2) guarantee that each customer is visited exactly once. Constraints (3) impose the condition that each customer is served by the same vehicle. Constraints (4) ensure that each vehicle can serve one route at most. Constraints (5) guarantee that there is no demand picked up in the beginning of routes. Conversely, constraints (6) impose the condition that there is no demand for delivery in the end of routes. Constraints (7) and (8) are flow equations for pickup and delivery loads, respectively. Constraints (9) and (10), respectively guarantee that the sum of the inflow to the origin is equal to the total pickup demand and that the sum of the outflow from the origin is equal to the total delivery demand. Constraints (11) ensure that the capacity of the vehicles is respected throughout the routes. Finally, constraints (12) specify the domain of variables. The MILP model (1)-(12) contains \(O(mn^{2})\) binary variables and \(O(mn)\) constraints.