## 3 Problem Description and Mathematical Modeling
### Mathematical Formulation  
The MIP model based on node clones is described as follows:  
\[Min c_{f}\sum_{i\in DD}\sum_{j\in V^{\prime}}x_{ij}+c_{cs}\sum_{i\in C^{\prime }}\theta_{i}+c_{bs}\sum_{i\in S^{\prime}}w_{i}+c_{v}\sum_{i\in V^{\prime}_{DD} }\sum_{j\in V^{\prime}_{AD}}d_{ij}x_{ij} \tag{1}\] \[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\] (2) \[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\] (3) \[\sum_{j\in V^{\prime}}x_{ij}\leq 1,\forall i\in DD\] (4) \[\sum_{i\in V^{\prime}}x_{ij}\leq 1,\forall j\in AD\] (5) \[\sum_{i\in DD}\sum_{j\in V^{\prime}}x_{ij}=\sum_{j\in V^{\prime}}\sum_{i\in AD }x_{ji}\] (6) \[\sum_{i\in V^{\prime}_{DD}}x_{ij}=\sum_{i\in V^{\prime}_{AD}}x_{ji}, \forall j\in V^{\prime},i\
eq j\] (7) \[0\leq u_{j}\leq u_{i}-q_{i}x_{ij}+q_{v}\big{(}1-x_{ij}\big{)},\ \forall i\in V^{\prime}_{DD},\ \forall j\in V^{\prime}_{AD},i\
eq j\] (8) \[0\leq u_{i}\leq q_{v},\forall i\in DD\] (9) \[\tau_{i}\leq T_{max},\ \forall i\in AD\] (10) \[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]0 (11) \[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]1 (12) \[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]2 (13) \[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]3\[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]4  
\[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]5  
\[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]6  
\[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]7  
\[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]8  
\[s.t. \sum_{j\in V^{\prime}_{AD}}x_{ij}=1,\forall i\in V,i\
eq j\]9  
\[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\]0  
\[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\]1  
\[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\]2  
\[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\]3  
\[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\]4  
\[\sum_{j\in V_{AD}}x_{ij}\leq 1,\forall i\in F\]5  
Objective Function (1) minimizes the total distribution cost, including vehicle utilization, energy replenishment, and traveling distance expenses. Constraint (2) ensures each customer node is visited exactly once, while Constraint (3) restrict each cloned charging pump and BSI to at most one visit. Constraints (4) and (5) specify that vehicles depart from the depot to a single destination and return from a unique node, respectively. Constraint (6) balances the number of vehicles departing and returning to the depot. Constraint (7) enforces flow conservation to maintain route continuity. Constraints (8) and (9) manage vehicle load capacity, ensuring that loads at subsequent customer points do not exceed previous ones and that total loads do not surpass vehicle capacity. Constraint (10) ensures vehicles return within the driver's shift time. Constraints (10)-(13) track service durations at customer points, charging pumps, and BSI, respectively. Constraints (14)-(16) define the remaining battery level upon arrival at the next location after departing from customer points, charging pumps, and BSI, respectively. Constraints (17) and (18) specify the remaining battery level when vehicles arrive at and depart from charging pumps and BSI. Constraint (19) ensures vehicles depart the depot with fully charged batteries. Constraint (20) regulates energy replenishment at cloned charge pumps. Constraints (21)-(23) impose capacity limitations on CL. Constraint (21) ensures that the arrival time at a subsequent clone of the same charging pump is at least after the charging duration of the previous clone at the previous arrival time. Constraint (22) enforces the use of cloned charging pumps in ascending numerical order to avoid equivalent solutions. Constraint (23) prioritizes the use of charging pumps with smaller serial numbers within the same CI. Finally, Constraints (24)-(26) define the nature of the decision variables.