Optimizing Multi-Depot Mixed Fleet Vehicle–Drone Routing Under a Carbon Trading Mechanism

Abstract

The global pursuit of carbon neutrality requires the reduction of carbon emissions in logistics and distribution. The integration of electric vehicles (EVs) and drones in a collaborative delivery model revolutionizes last-mile delivery by significantly reducing operating costs and enhancing delivery efficiency while supporting environmental objectives. This paper presents a cost-minimization model that addresses transportation, energy, and carbon trade costs within a cap-and-trade framework. We develop a multi-depot mixed fleet, including electric and fuel vehicles, and a drone collaborative delivery routing optimization model. This model incorporates key factors such as nonlinear EV charging times, time-dependent travel conditions, and energy consumption. We propose an adaptive large neighborhood search algorithm integrating spatiotemporal distance (ALNS-STD) to solve this complex model. This algorithm introduces five domain-specific operators and an adaptive adjustment mechanism to improve solution quality and efficiency. Our computational experiments demonstrate the effectiveness of the ALNS-STD, showing its ability to optimize routes by accounting for both spatial and temporal factors. Furthermore, we analyze the influence of charging station distribution and carbon trading mechanisms on overall delivery costs and route planning, underscoring the global significance of our findings.

1. Introduction

The logistics industry, a vital contributor to economic growth, faces considerable challenges due to its inherent complexity. With the rapid pace of urbanization, the demand for distribution services is increasing sharply. The global urban population, which was 55% in 2018, is expected to reach 68% by 2050 [1]. This urban expansion, combined with population growth, has driven a dramatic surge in parcel deliveries, increasing from 43 billion in 2014 to 131 billion in 2020 across 13 major markets, including the U.S., Brazil, and China [2]. As customer demand for delivery services grows and urban logistics activities intensify, the environmental impact of expanding logistics operations, both in terms of scale and cost, has become a pressing concern.

Despite global efforts toward carbon neutrality and zero-carbon transitions, CO₂ emissions continued to increase in 2022, particularly in countries such as China, India, and Canada. In recent years, emissions from key countries, including China, the U.S., and India, have remained persistently high [3]. In response to global goals for carbon neutrality and growing competition within the industry, the logistics and distribution sector must continuously seek ways to improve profitability while minimizing total costs, including those associated with carbon trading.

Effective vehicle route planning under the carbon trading mechanism is essential to ensure the sustainability of the logistics industry. Effective vehicle path planning under a carbon trading mechanism is crucial for promoting sustainability in the logistics industry. Guo et al. [4] developed a time window green path model for cold chain logistics to minimize transportation, refrigeration, carbon emissions, and labor costs, using a two-stage hybrid search algorithm to solve the problem. Yao et al. [5] explored the impacts of carbon tax and freshness costs, proposing a green vehicle path model for agricultural products and employing an improved ant colony optimization algorithm (IACO) to solve it. Liu et al. [6] constructed a low-carbon distribution route optimization model based on an improved genetic algorithm for community group purchasing, focusing on minimizing total costs and carbon emissions while incorporating soft time windows and emission parameters. Zhang et al. [7] developed a green and low-carbon distribution model for supermarket chains to minimize fixed transportation and carbon emission costs. Their model integrates carbon emission costs and proposes a heuristic algorithm to optimize the vehicle routing problem (VRP).

While these studies focus on optimizing vehicle paths and improving operational efficiency in single distribution center scenarios, the capacity-constrained multi-distribution center vehicle routing problem (MDCVRP) applies more to real-world logistics, supply chain management, and transportation scheduling. Jabir et al. [8] incorporated environmental impacts into a multi-site green vehicle routing problem, aiming to minimize both economic and emission reduction costs. They used LINGO for small-scale integer linear programming and developed ant colony optimization (ACO)-based algorithms to efficiently solve larger instances. Wang et al. [9] optimized the multi-site green vehicle routing problem by applying the CMEM microscopic fuel consumption model, which accounts for instantaneous speed and load, and solved it using a multi-objective particle swarm algorithm. Fan et al. [10] addressed the multi-site vehicle routing problem within a time-varying network, aiming to minimize total costs, including fuel expenses, by employing a time-varying particle swarm algorithm that optimized routes through spatiotemporal distance clustering. Xue [11] explored the impact of carbon emissions on multi-warehouse distribution networks, proposing a method to quantify the relationship between fuel consumption and carbon emissions and a two-stage heuristic framework for path optimization.

MDCVRP helps companies optimize resource allocation and reduce costs, while electric vehicle distribution, under carbon emission restrictions, offers an effective solution to promote environmental sustainability. The electric vehicle routing problem (EVRP), a key variant of the traditional vehicle routing problem, introduces additional constraints such as battery capacity and charging operations [12]. Froger et al. [13] modeled a nonlinear charging function using a concave linear approach, developing an EVRP optimization model that addresses real-world distribution challenges. Zhang et al. [14] incorporated uncertainties in service times, power consumption, and travel times across different regions and customer types, representing these uncertainties with fuzzy numbers and constructing a fuzzy optimization model based on plausibility theory. Lee [15] proposed an optimization method to simultaneously reduce EV path mileage and charging time, considering a nonlinear charging time function. Wang et al. [16] designed a multi-depot EV routing model with time windows and shared charging stations using an enhanced multi-objective genetic algorithm to optimize both path selection and charging station placement.

Under carbon reduction policies, the adoption of new energy technologies has accelerated. However, as many companies cannot fully transition their energy assets, hybrid fleet distribution systems have emerged as a practical solution for logistics companies. Ene et al. [17] introduced the green vehicle routing problem (GVRP) for heterogeneous fleets to reduce fuel consumption and CO₂ emissions. They developed a hybrid metaheuristic algorithm (HMA) to assess heterogeneous fleets’ impacts on GVRP variants. Islam et al. [18] designed a hybrid metaheuristic that combines ant colony optimization (ACO) with variable neighborhood search (VNS) to address distribution challenges involving a mixed fleet of conventional and green vehicles under a carbon emission cap. Al-Dalain and Celebi [19] proposed a hybrid fuel-electric fleet model for single-vehicle yard path planning, focusing on minimizing total operating costs while optimizing both distribution routes and fleet configurations. Yu et al. [20] explored optimizing the green path for hybrid fuel-electric fleets in a yard environment, incorporating actual energy consumption data and charging strategies, and solving the problem with an adaptive large neighborhood search algorithm. Amiri et al. [21] developed a bi-objective planning model to minimize total transportation costs and greenhouse gas emissions, tackling the path optimization problem for heavy-duty electric and conventional trucks and employing a hybrid adaptive neighborhood search algorithm integrating three multi-objective methods.

Developing a low-altitude economy has become a promising strategy for easing traffic congestion and reducing environmental pollution. As a key driver of this economy, drones are increasingly valuable in addressing the last-mile delivery challenge due to their flexibility, efficiency, and environmentally friendly nature. Drone-assisted vehicle delivery has become increasingly prevalent in last-mile logistics within the supply chain. Compared to traditional vehicle-only delivery solutions, incorporating drones into last-mile operations can significantly enhance efficiency, reducing customer waiting times by up to 60% [22]. The drone-assisted vehicle delivery routing problem (VRPD) extends the traditional vehicle routing problem by incorporating drones into delivery operations. Chiang et al. [23] developed a mixed integer green routing model for drones, showing that optimizing drone delivery routes is cost-effective and energy-efficient, significantly reducing CO₂ emissions. Han et al. [24] proposed a VRPTW model involving drones, aiming to minimize the weighted sum of total truck and drone energy and the total number of trucks, using an improved artificial swarm algorithm and introducing a scout bee strategy to enhance global search capability. Zhang et al. [25] introduced a multi-objective optimization model integrating drones and trucks to enhance delivery efficiency and reduce energy consumption. They utilized an extended non-dominated sorting genetic algorithm to solve the drone delivery route problem. Xiao and Gao [26] tackled the problem of drone delivery routing by developing a mixed-integer planning model that minimizes fixed transportation and carbon trading costs. They also proposed a two-stage heuristic algorithm to address carbon emissions in truck and drone transportation.

The VRPD significantly improves distribution efficiency, reduces costs, alleviates traffic congestion, and mitigates environmental impact. Electric vehicles (EVs) and drones offer environmentally friendly alternatives with low carbon emissions. Their coordinated use reduces dependence on traditional fuel vehicles, lowering the carbon footprint. Baek et al. [27] proposed a coordinated distribution system for EVs and drones focused on energy savings and emissions reductions, assuming that each drone serves only one customer at a time and battery replacement is impossible during operation. Zhu et al. [28] demonstrated that, considering the need for EVs to visit charging stations, the adaptive large neighborhood search algorithm outperforms the variable neighborhood search algorithm in EV-drone cooperative delivery scenarios. Kyriakakis et al. [29] developed a mixed-integer planning model that optimizes energy consumption for the problem of EV-drone co-delivery path, solved using a hybrid ant colony optimization algorithm. Building on this research, Mara et al. [30] proposed a comprehensive co-delivery model where an EV carries a drone for joint delivery, with the EV visiting charging stations to ensure route feasibility.

In constructing logistics distribution models, considering only the spatial distances between customers while neglecting temporal factors can lead to an incomplete assessment of the logistics system’s actual operation. Thus, developing an optimization model that integrates customer service times with spatial locations is crucial. Such a model provides a more accurate and comprehensive reflection of the system’s operational dynamics, enabling the formulation of more effective vehicle scheduling strategies, reducing delivery vehicle operating costs, and improving the efficiency and service quality of logistics distribution. For example, Wang et al. [31] proposed a multi-objective vehicle routing problem with time windows for perishables (MO-VRPTW-P), aiming to ensure the freshness of perishable goods and minimize value loss. They utilized a two-stage heuristic algorithm based on a Pareto variable neighborhood search genetic algorithm (STVNS-GA), incorporating spatiotemporal distances to solve the problem. Similarly, Hou et al. [32] addressed the heterogeneous vehicle routing problem in multi-depot networks, considering time-dependent road networks and soft time windows. They generated an initial solution by clustering customers based on spatiotemporal distances and employed a variable neighborhood search algorithm within a genetic algorithm for local optimization. Their approach was validated through examples, offering theoretical guidance for logistics enterprises in optimizing delivery plans. Wang and Li [33] tackled the vehicle routing problem with time windows and spatiotemporal distances (VRPTWTSD) by proposing an improved genetic algorithm (IGA). Their approach included a hybrid initialization method combining the temporal-spatial insertion heuristic (TSDIH) and the earliest ready time heuristic (ERH). Additionally, they developed two knowledge-based crossover operators to expand the search space during encoding. Statistical analyses validated the effectiveness of their algorithm, further advancing the theoretical and practical understanding of spatiotemporal integration in logistics optimization.

The existing literature provides valuable insights for this study; however, several critical gaps remain: (1) the development of a green distribution system that integrates mixed fleets and drones in the context of corporate energy asset transitions has not been fully addressed; (2) the multi-depot cooperative distribution model for mixed fleets and drones under carbon emission reduction targets remains unexplored; and (3) the comprehensive impact of spatial and temporal distances on transportation costs has yet to be thoroughly investigated.

This paper proposes a multi-depot collaborative distribution model that integrates a mixed fleet and drones within a carbon trading mechanism to bridge these gaps. The mixed fleet consists of traditional fuel vehicles and EVs, accounting for crucial factors such as resource allocation, customer time windows, EV charging strategies, energy consumption, time-dependent road networks, and the technical specifications of different transport modes. We apply a normalization process using weighted averages for temporal and spatial distances and employ an ALNS algorithm to optimize delivery routes. This approach enables logistics companies to improve cargo transportation efficiency while promoting low-carbon, green, and cost-effective operations. It also helps businesses reduce logistics costs, enhancing overall distribution efficiency.

The remainder of the paper is organized as follows. Section 2 introduces the model and defines the relevant parameters and variables; Section 3 details the design of the proposed ALNS-STD algorithm; Section 4 presents the algorithm’s validation using computational experiments; and Section 5 concludes with a summary of the findings.

2. Problem Description and Modeling Framework

2.1. Problem Description

In the time-dependent road network illustrated in Figure 1, there are multiple distribution centers, customer nodes, drones, EVs, and fuel vehicles, all with identical specifications. Both electric and fuel vehicles are equipped with drones to perform delivery tasks subject to capacity and time window constraints. While vehicles serve customers, drones can deliver to other customers, each being served only once. After completing their deliveries, the drones rejoin their designated vehicles at the next customer node. The mixed fleet, consisting of EVs and fuel vehicles, completes deliveries and returns to any distribution center. During delivery, EVs must visit charging stations to recharge and ensure sufficient power for continued operation. Two factors significantly impact the charging strategy and the efficiency of vehicle–drone coordination: (1) EV charging follows a nonlinear function influenced by variables such as battery performance and heat dissipation; and (2) during EV charging, the station can serve as a takeoff and landing platform for drones, regardless of customer demand or time window, acting as a special node. However, drones carried by fuel vehicles are not allowed to use charging stations as take-off or landing points.

This study defines the multi-depot mixed fleet–drone collaborative delivery routing optimization problem with time windows as MD-MFVRPDTW. The following assumptions are made for the proposed MD-MFVRPDTW model:

• All vehicles of the same type have identical specifications, and each vehicle is equipped with a drone.

• All delivery tasks begin and end at a designated distribution center.

• Drones are launched and recovered exclusively by their assigned vehicle, and these operations occur at specific nodes (distribution centers, customer nodes, or charging stations).

• Vehicles are fully charged or fueled at the start of the route, and energy consumption during customer service is not considered.

• Charging stations are standardized in equipment and operation, offering a full charge strategy but with limited availability for charging vehicles.

• Drones always begin flights fully charged and follow an instantaneous power change mode during operation.

• Each customer receives one delivery service from a vehicle or drone.

• The distribution capacity of each distribution center is sufficient to meet the total demand of all customers, and all centers are treated as homogeneous.

• There is no speed difference between fuel-powered and electric vehicles.

2.2. Notation

In this section, we present the notation for the MD-MFVRPDTW model within the context of a carbon trading framework. Table 1, Table 2, Table 3 and Table 4 define and detail the relevant sets, model parameters, decision variables, and non-decision variables, along with their corresponding descriptions. These components form the foundation of the MD-MFVRPDTW modeling framework, allowing a clear and comprehensive representation of the problem structure.

We use a tuple <i, j, k> to represent a drone’s flight route, where node i represents the drone’s launch node, node j represents the customer node served by the drone, and node k represents the drone’s recovery node.

2.3. Nonlinear Charging Time for EVs

The charging rate of electric vehicles (EVs) varies over time due to environmental factors such as temperature. To capture this more accurately, we adopt the nonlinear charging model proposed by Montoya et al. [34], offering a more realistic representation of the charging process. As illustrated in Figure 2, let ${\alpha }_k$ and $t_k$ represent the state of charge and the charging time at an inflection point, respectively. Let $B_t$ denote the full charge capacity, and $T_e$ the total time required for a full charge. Extracting the endpoint coordinates of the piecewise charging function from Montoya et al. [34] reveals the following: when the charge reaches 83.99% of $B_t$, the charging time is 59.9% of $T_e$. At 94.41% of $B_t$, the time increases to 73.05% of $T_e$. And, at 100% of $B_t$, the full charging time $T_e$ is required. By presetting the linear charging rate for the first segment and the full charge capacity $B_t$, the piecewise function describing the relationship between the state of charge and charging time can be defined. When an EV $f$ arrives at charging station $i$ with a remaining charge of $e_{if}^{(1)}(f\in F,i\in V_S)$, the charging time $t{c}_{if}$ for the EV can be calculated using Equation (1).

(1)$t{c}_{if}=\left\{\begin{array}{l}{T}_e-{\displaystyle \frac{e_{if}^{(1)}*t_1}{{\alpha }_1}},e_{if}^{(1)}\in [0,{\alpha }_1)\\ T_e-{\displaystyle \frac{(e_{if}^{(1)}-{\alpha }_1)*(t_2-t_1)}{{\alpha }_2-{\alpha }_1}}-t_1,e_{if}^{(1)}\in [{\alpha }_1,{\alpha }_2)\\ T_e-{\displaystyle \frac{(e_{if}^{(1)}-{\alpha }_2)*(t_3-t_2)}{{\alpha }_3-{\alpha }_2}}-t_2,e_{if}^{(1)}\in [{\alpha }_2,{\alpha }_3)\end{array}\right.$

2.4. Time-Dependent Vehicle Travel Time

Vehicle travel time varies throughout the day due to uncertainties in the vehicle’s speed under real-world conditions. This study assumes that the difference in speed between fuel and electric vehicles is negligible within an urban road network. Consequently, we model the time-dependent travel time for both vehicle types using a unified approach, following the method proposed by Ichoua et al. [35]. Specifically, travel time is represented as a continuous piecewise function that fluctuates based on the time of day.

2.5. Vehicle Energy Consumption

Our study uses the fuel consumption model proposed by Demir et al. [36] within the context of the pollution routing problem (PRP). Fuel consumption during transit is calculated based on the actual mechanical power of the fuel vehicle. The model assumes that the constant mechanical power $p_{ij}$ of a fuel vehicle traveling from node $i$ to node $j$ is a binary function of the load $l_{ij}$ and speed $v_{ij}$, which can be approximated by Equation (2). This approximation disregards the effects of road gradients and urban intersections.

(2)$P_{ij}=0.5A_f\rho c_d{v_{ij}}^3+(m_f+l_{ij})g{c}_r{v}_{ij}$

In Equation (1), $A_f$ is the windward area of the fuel vehicle, $\rho$ represents air density, $c_d$ is the air resistance coefficient, $g$ denotes the universal gravitational constant, $c_r$ is the rolling resistance coefficient, and $m_f$ is the empty weight of the fuel vehicle. The fuel consumption rate $G_{ij}$ and mechanical power $P_{ij}$ per unit time of fuel vehicle $f_2$ traveling between node $i$ and node $j$ are determined by Equation (3):

(3)$G_{ij}=\xi (\epsilon \varphi +P_{ij}/\varpi )/(\kappa \psi )$

where $\xi$ represents the fuel-to-air ratio, $\epsilon$ is the engine friction coefficient, $\varphi$ denotes the engine speed, $\varpi$ is the gasoline engine efficiency parameter, $\kappa$ is the calorific value of gasoline, and $\psi$ is the fuel rate conversion coefficient. The fuel consumption rate $G_{ij}$ per unit time is multiplied by the fuel vehicle’s travel time between nodes $i$ and $j$ to obtain the total fuel consumption, as shown in Equation (4).

(4)$C{F}_{ij}=G_{ij}{d}_{ij}/v_{ij}$

In a time-dependent road network, fuel consumption calculation must consider the time zone in which the fuel vehicle operates. Let $t_1$ represent the departure time of the fuel vehicle, $d_{ij}^t$ denote the distance between node $i$ and node $j$, and $t_2$ represent the arrival time. The $k$-th time zone is designated as $T_K$ with its upper and lower bounds represented as $\overline{t_k}$ and $t_k$, respectively. The travel speed is denoted as $v_{T_k}$. Based on these parameters, the fuel consumption $C{F}_{ij}$ between nodes $i$ and $j$ can be calculated as follows.

• (1). Scenario 1: Let $d_1=v_{T_k}(\overline{t_k}-t_1)$ represent the remaining distance the fuel vehicle can travel within the current time zone. If it is determined that $d_1\ge d_{ij}$ as illustrated in Figure 3, the calculation of the fuel consumption between nodes $i$ and $j$ is given by Equation (5).

  (5)$C{F}_{ij}=[\xi \epsilon \varphi /k\phi +0.5A_f\rho c_d\xi v_{T_k}^3/k\phi \varpi +g{c}_r\xi (m_f+l_{ij})v_{T_k}/k\phi \varpi ]d_{ij}/v_{T_k}$

• (2). Scenario 2: On the other hand, if $d_1\le d_{ij}$, it indicates a change in travel speed between the two nodes, as shown in Figure 4. In this scenario, the remaining distance $d_2=d_{ij}-d_1$, which the fuel vehicle can travel during the second phase starting at time $t_2$, must be calculated. The fuel consumption for the distance between nodes $i$ and $j$ is then expressed by Equation (6) below.

  (6)$\begin{array}{l}C{F}_{ij}=[\xi \epsilon \varphi /k\phi +0.5A_f\rho c_d\xi v_{T_k}^3/k\phi \varpi +g{c}_r\xi (m_f+l_{ij})v_{T_k}/k\phi \varpi ]d_1/v_{T_k}\\ +[\xi \epsilon \varphi /k\phi +0.5A_f\rho c_d\xi v_{T_k}^3/k\phi \varpi +g{c}_r\xi (m_f+l_{ij})v_{T_{k+1}}/k\phi \varpi ]d_2/v_{T_{k+1}}\end{array}$

2.6. Optimization Model

The objective function of the proposed MD-MFVRPDTW model is to minimize the overall delivery costs $Z$, which encompasses the transportation cost of EVs $Z_{transport}^{t1}$, the charging cost of EVs $Z_{charging}^{t1}$, the transportation cost of fuel vehicles $Z_{transport}^{t2}$, the transportation cost of drones $Z_{transport}^d$, the carbon trading cost $Z_{carbon}$, and the time window penalty cost $Z_{punishment}$. The calculation of these costs is detailed as follows.

(7)$Z_{transport}^{t1}=c_1{\displaystyle \sum _{f\in F_1}{\displaystyle \sum _{i\in N,i\ne j}{\displaystyle \sum _{j\in N}{d}_{ij}^t}}}{x}_{ij{f}_1}+c_{f_1}{\displaystyle \sum _{f\in F_1}{\displaystyle \sum _{i\in N_0}{\displaystyle \sum _{j\in N_V\cup N_0^{+}}{x}_{ij{f}_1}}}}$

(8) $Z_{ch\arg ing}^{t1}=c_2{\displaystyle \sum _{f\in F_1}{\displaystyle \sum _{i\in N_S}e{c}_{i{f}_1}}{y}_{i{f}_1}}$

(9) $Z_{transport}^{t2}=c_3{\displaystyle \sum _{f\in F_1}{\displaystyle \sum _{i\in N,i\ne j}{\displaystyle \sum _{j\in N}C{F}_{ij}}}}{x}_{ij{f}_2}+c_{f_2}{\displaystyle \sum _{f\in F_1}{\displaystyle \sum _{i\in N_0}{\displaystyle \sum _{j\in V_N\cup N_0^{+}}{x}_{ij{f}_2}}}}$

(10) $Z_{transport}^d=c_4{\displaystyle \sum _{f\in F_r}{\displaystyle \sum _{r\in R}{\displaystyle \sum _{i\in N_L}{\displaystyle \sum _{j\in N_R}{d}_{ij}^d}}}}{x}_{ijk{f}_r}$

(11) $Z_{carbon}=cp(EF{\displaystyle \sum _{f\in F_2}{\displaystyle \sum _{i\in N,i\ne j}{\displaystyle \sum _{j\in N}C{F}_{ij}{x}_{ij{f}_2}}}}-\overline{CE})$

(12) $Z_{punishment}=({\delta }_1{\displaystyle \sum _{i\in N_C}{\displaystyle \sum _{k\in R_r}\max \{0,E{T}_i-a_{i{f}_r}^t\}+{\delta }_2{\displaystyle \sum _{i\in N_C}{\displaystyle \sum _{k\in R_r}\max \{0,a_{i{f}_r}^t-ET\}}}}})$

Therefore, our MD-MFVRPDTW model can be formulated as follows.

Objective function:

(13)$Z=Z_{transport}^{t1}+Z_{charging}^{t1}+Z_{transport}^{t2}+Z_{transport}^d+Z_{carbon}+Z_{punishment}$

Constraints:

(14)${\displaystyle \sum _{i\in N_L\backslash \left\{j\right\}}{x}_{ij{f}_r}}+{\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{\displaystyle \sum _{k\in N_R\backslash \{i,j\}}{y}_{ijk{f}_r}}=1},\forall f_r\in F_r,r\in R,j\in N_C$

(15)${\displaystyle \sum _{j\in N_V\cup N_0^{+}}{x}_{ij{f}_r}}\le 1,\forall f_r\in F_r,r\in R,i\in N_0$

(16)${\displaystyle \sum _{i\in N_0\cup N_V}{x}_{ij{f}_r}}\le 1,\forall f_r\in F_r,r\in R,j\in N_0^{+}$

(17)${\displaystyle \sum _{i\in N_L\backslash \left\{j\right\}}{x}_{ij{f}_r}-}{\displaystyle \sum _{k\in N_R\backslash \left\{j\right\}}{x}_{jk{f}_r}}=0,\forall f_r\in F_r,r\in R,j\in N_V$

(18)$u_{j{f}_r}\le M{\displaystyle \sum _{i\in N_L\backslash \left\{j\right\}}{x}_{ij{f}_r}},\forall f_r\in F_r,r\in R,j\in N_R$

(19)$u_{i{f}_r}-u_{j{f}_r}\le M(1-x_{ij{f}_r})-1,\forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{j\right\},j\in N_R\backslash \left\{i\right\}$

(20)$u_{j{f}_r}-u_{i{f}_r}\le M{p}_{ij{f}_r},\forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{j\right\},j\in N_V\backslash \left\{i\right\}$

(21)$u_{j{f}_r}-u_{i{f}_r}\le M(p_{ij{f}_r}-1)+1,\forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{j\right\},j\in N_V\backslash \left\{i\right\}$

(22)${\displaystyle \sum _{j\in N_C}{\displaystyle \sum _{k\in N_R\backslash \left\{j\right\}}{q}_j{x}_{jk{f}_r}}+}{\displaystyle \sum _{j\in N_C}{\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{\displaystyle \sum _{k\in N_R\backslash \{i,j\}}{q}_j{y}_{ijk{f}_r}}}}\le Q_r,\forall f_r\in F_r,r\in R$

(23)$y_{ijk{f}_r}=0,\forall f_r\in F_r,r\in R,i\in N_L\backslash \{j,k\},j\in N_C\backslash \{N_D,i,k\},k\in N_R\backslash \{i,j\}$

(24)${\displaystyle \sum _{j\in N_D\backslash \{i,k\}}{\displaystyle \sum _{k\in N_R\backslash \{i,j\}}{y}_{ijk{f}_r}}}\le 1,\forall f_r\in F_r,r\in R,i\in N_L$

(25)${\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{\displaystyle \sum _{j\in N_D\backslash \{i,k\}}{y}_{ijk{f}_r}}}\le 1,\forall f_r\in F_r,r\in R,k\in N_R$

(26)$$\begin{gathered} {\displaystyle \sum _{n\in N_R\backslash \left\{i\right\}}{x}_{{\inf }_r}+}{\displaystyle \sum _{n\in N_L\backslash \left\{k\right\}}{x}_{nk{f}_r}}\ge 2y_{ijk{f}_r} \\ \forall f_r\in F_r,r\in R,i\in N_L\backslash \{j,k\},j\in N_C\backslash \{i,k\},k\in N_R\backslash \{i,j\} \end{gathered}$$

(27)$$\begin{gathered} a_{i{f}_r}^t+t_{ik}^t+s_k^t+lt{\displaystyle \sum _{j\in N_D\backslash \{i,n\}}{\displaystyle \sum _{n\in N_R\backslash \{i,j\}}{y}_{ijn{f}_r}}}+rt{\displaystyle \sum _{n\in N_L\backslash \{j,k\}}{\displaystyle \sum _{j\in N_D\backslash \{n,k\}}{y}_{njk{f}_r}}}-M*(1-x_{ik{f}_r})\le a_{k{f}_r}^t \\ \forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{k\right\},k\in N_R\backslash \{i,N_S\} \end{gathered}$$

(28)$$\begin{gathered} a_{i{f}_r}^t+t_{ik}^t+lt{\displaystyle \sum _{j\in N_D\backslash \{i,n\}}{\displaystyle \sum _{n\in N_R\backslash \{i,j\}}{y}_{ijn{f}_r}}}+rt{\displaystyle \sum _{n\in N_L\backslash \{j,k\}}{\displaystyle \sum _{j\in N_D\backslash \{n,k\}}{y}_{njk{f}_r}}}+t{c}_{k{f}_1}-M*(1-x_{ik{f}_r})\le a_{k{f}_r}^t \\ \forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{k\right\},k\in (N_R\cap N_S)\backslash \left\{i\right\} \end{gathered}$$

(29)$$\begin{gathered} a_{i{f}_r}^t+lt+t_{ij}^d+s_i^d-M(1-{\displaystyle \sum _{k\in N_R\backslash \{i,j\}}{y}_{ijk{f}_r}})\le a_{i{f}_r}^d \\ \forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{j\right\},j\in N_D\backslash \left\{i\right\} \end{gathered}$$

(30)$$\begin{gathered} a_{j{f}_r}^d+t_{jk}^d+rt+s_k^t-M(1-{\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{y}_{ijk{f}_r}})\le a_{k{f}_r}^d \\ \forall f_r\in F_r,r\in R,j\in N_D\backslash \left\{k\right\},k\in N_R\backslash \{j,N_S\} \end{gathered}$$

(31)$$\begin{gathered} a_{j{f}_r}^d+t_{jk}^d+rt+t{c}_{k{f}_1}-M(1-{\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{y}_{ijk{f}_r}})\le a_{k{f}_r}^d \\ \forall f_r\in F_r,r\in R,j\in V_D\backslash \left\{k\right\},k\in (V_R\cap V_S)\backslash \left\{j\right\} \end{gathered}$$

(32)$a_{i{f}_r}^t-M(1-{\displaystyle \sum _{j\in N_D\backslash \{i,k\}}{\displaystyle \sum _{k\in N_R\backslash \{i,j\}}{y}_{ijk{f}_r}}})\le a_{i{f}_r}^d,\forall f_r\in F_r,r\in R,i\in N_L$

(33)$a_{i{f}_r}^t+M(1-{\displaystyle \sum _{j\in N_D\backslash \{i,k\}}{\displaystyle \sum _{k\in N_R\backslash \{i,j\}}{y}_{ijk{f}_r}}})\le a_{i{f}_r}^d,\forall f_r\in F_r,r\in R,i\in N_L$

(34)$a_{k{f}_r}^t-M(1-{\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{\displaystyle \sum _{j\in N_D\backslash \{i,k\}}{y}_{ijk{f}_r}}})\le a_{k{f}_r}^d,\forall f_r\in F_r,r\in R,k\in N_R$

(35)$a_{k{f}_r}^t+M(1-{\displaystyle \sum _{i\in N_L\backslash \{j,k\}}{\displaystyle \sum _{j\in N_D\backslash \{i,k\}}{y}_{ijk{f}_r}}})\le a_{k{f}_r}^d,\forall f_r\in F_r,r\in R,k\in N_R$

(36)$$\begin{gathered} a_{k{f}_r}^d+M(3-p_{im{f}_r}-{\displaystyle \sum _{i\in N_C\backslash \{i,k\}}{y}_{ijk{f}_r}}-{\displaystyle \sum _{p\in N_C\backslash \{m,q\}}{\displaystyle \sum _{q\in N_R\backslash \{m,p\}}{y}_{mpq{f}_r}}})\le a_{m{f}_r}^d \\ \forall f_r\in F_r,r\in R,i\in N_L\backslash \{m,k\},m\in N_V\backslash \{i,k\},k\in N_R\backslash \{i,m\} \end{gathered}$$

(37)$E{T}_i\le a_{i{f}_r}^t\le L{T}_i,\forall f_r\in F_r,r\in R,i\in N_C-N_D$

(38)$E{T}_i\le a_{i{f}_r}^d\le L{T}_i,\forall f_r\in F_r,r\in R,i\in N_D$

(39)$l_{ij}-q_j-M(1-x_{ij{f}_2})\le l_{jk},\forall f_2\in F_2,i\in N_L\backslash \left\{j\right\},j\in N_C\backslash \left\{i\right\},k\in N_R\backslash \{i,j\}$

(40)$t{c}_{i{f}_1}=U_{i{f}_1}(e_{i{f}_1}^{(1)}),\forall f_1\in F_1,i\in N_S$

(41)$e{c}_{i{f}_1}=(B_{\mathrm{t}}-e_{i{f}_1}^{(1)})x_{i{f}_1},\forall f_1\in F_1,i\in N_S$

(42)$e_{i{f}_1}^{(2)}-h{d}_{ij}-M(1-x_{ij{f}_1})\le e_{j{f}_1}^{(1)},\forall f_1\in F_1,i\in N_L\backslash \left\{j\right\},j\in N_R\backslash \left\{i\right\}$

(43)$e_{i{f}_1}^{(1)}\ge 0,\forall f_1\in F_1,i\in N$

(44)$e_{i{f}_1}^{(2)}=B_{\mathrm{t}},\forall f_1\in F_1,i\in N_R$

(45)$e_{i{f}_1}^{(1)}=e_{i{f}_1}^{(2)},\forall f_1\in F_1,i\in N_C$

(46)$e_{i{f}_1}^{(2)}=B_{\mathrm{t}},\forall f_1\in F_1,i\in N_0$

(47)$a_{i{f}_r}^t=0,\forall f_r\in F_r,r\in R,i\in N_0$

(48)$a_{i{f}_r}^d=0,\forall f_r\in F_r,r\in R,i\in N_0$

(49)$$\begin{gathered} t_{ij}^d+s_i^d+t_{jk}^d-M*(1-y_{ijk{f}_r})\le B_d \\ \forall f_r\in F_r,r\in R,i\in N_L\backslash \{j,k\},j\in N_D\backslash \{i,k\},k\in N_R\backslash \{i,j\} \end{gathered}$$

(50)$x_{ij{f}_r}\in \{0,1\},\forall f_r\in F_r,r\in R,i,j\in A$

(51)$y_{ijk{f}_r}\in \{0,1\},\forall f_r\in F_r,r\in R,i\in N_L\backslash \{j,k\},j\in N_D\backslash \{i,k\},k\in N_R\backslash \{i,j\}$

(52)$p_{ij{f}_r}\in \{0,1\},\forall f_r\in F_r,r\in R,i\in N_L\backslash \left\{j\right\},j\in N_V\backslash \left\{i\right\}$

(53)$u_{i{f}_r}\in N,\forall f_r\in F_r,r\in R,i\in V$

Equation (13) represents the objective function that aims to minimize the total delivery costs. Equation (14) ensures that each customer is served exactly once, either by a vehicle or a drone. Constraints (15) and (16) require that each route start and end at any distribution center. Equation (17) enforces node flow conservation, while Constraints (18) and (19) prevent the formation of subtours. Constraints (20) and (21) impose restrictions on vehicle route sequencing to avoid backtracking. Constraint (22) sets limits on vehicle payload capacity, and Equation (23) stipulates that drones can only visit nodes that satisfy their delivery criteria. Constraints (24) and (25) ensure that the drone’s launch and recovery nodes are used only once, and Constraint (26) requires that the drone be launched and recovered at two distinct nodes visited by its assigned vehicle.

Constraints (27) and (28) ensure that the EV visits a charging station, while Constraints (29) and (30) regulate the EV’s arrival time at both charging and non-charging station nodes, ensuring that it exceeds the time at the preceding node. Constraint (31) requires that a drone visit a customer node only after the vehicle arrives at the launch point. Constraints (32) and (33) synchronize the arrival times of the vehicle and the drone at the launch node, while Constraints (34) and (35) enforce synchronization at the recovery node. Constraint (36) limits each vehicle to carrying only one drone, and Constraints (37) and (38) enforce time window restrictions for both the vehicle and the drone.

Constraint (39) calculates the fuel vehicle’s load between any two nodes, while Equation (40) determines the EV’s charging time. Equation (41) calculates the amount of charging required for the EV at a charging station, and Constraint (42) assesses the EV’s power levels between any two nodes. Constraint (43) ensures that the EV’s power remains above zero, while Equation (44) prescribes that the EV follows a full charging strategy. Equation (45) indicates that the EV does not consume power while serving a customer, and Equation (46) ensures that the EV departs from the distribution center fully charged. Equations (47) and (48) state that the mixed fleet and drones depart from the distribution center at time zero. Lastly, Constraint (49) enforces the drone’s range limitations, while Constraints (50)–(53) define the attributes of the decision variables.

3. Algorithm Design

3.1. Algorithm Design for ALNS-STD

In large-scale distribution problems, relying solely on the Euclidean distance often results in the separation of spatial and temporal distances between customers. This issue, when overlooked in the basic adaptive large neighborhood search (ALNS) framework, can lead to route fragmentation and the generation of inefficient vehicle paths. To overcome this, we propose an improved adaptive large neighborhood search that takes into account spatiotemporal distance (ALNS-STD) to optimize routing. To ensure the algorithm’s robustness, several strategies are employed. The integration of the Metropolis criterion from the simulated annealing algorithm helps prevent premature convergence to local optima by allowing occasional acceptance of suboptimal solutions, thereby maintaining diversity in the search process. Additionally, the implementation of an unmodified threshold strategy provides stability by setting consistent evaluation benchmarks, while a random scoring system introduces stochasticity to further diversify the search space and improve the algorithm’s ability to explore a broader solution landscape. The proposed algorithm integrates key factors, including nonlinear EV charging times, time-dependent vehicle travel characteristics, and EV mechanical power consumption. Algorithm 1 presents the pseudo-code for the ALNS-STD algorithm, specifically designed to address the unique characteristics of the MD-MFVRPDTW model within the carbon trading mechanism framework.

Algorithm 1 Adaptive Large Neighborhood Search Algorithm Integrating Temporospatial Distance
Input: Set of customer nodes $N_C=\{m+1,\dots ,m+n\}$, set of departure distribution centers $N_0=\{0,1,\dots ,m\}$, set of charging stations $N_S=\{m+n+1,\dots ,m+n+r\}$, set of operators $N=\{N_{a,}a=1,2,3,4,5,6\}$, customer location, customer time window.
Output: $S_{best}$
1:	Construct initial solution based on spatiotemporal distance $x_0$
2:	$x_c\leftarrow x_0$, $x_b\leftarrow x_0$
3:	Repeat
4:	Choose the operator $N_a\in N$ according to the probability, $\sigma =\{{\sigma }_a,1,2,3,4,5,6\}$
5:	$x_n\leftarrow X(x_c)$, Perform field operations on $x_0$
6:	if $\exists$$f(x_n)\le f(x_c)$ then
7:	$x_c\leftarrow x_n$
8:	if $\forall$$f(x_c)\le f(x_b)$ then
9:	$x_b\leftarrow x_c$
10:	Update the score of the corresponding operator $N_a$, ${\varpi }_a={\varpi }_a+{\varpi }_1$
11:	Else
12:	Update the score of the corresponding operator $N_a$$,{\varpi }_a={\varpi }_a+{\varpi }_2$
13:	end if
14:	Else
15:	if $p_n$ < $p_{\omega }$ then
16:	$x_c\leftarrow x_n$
17:	Update the score of the corresponding operator $N_a$, ${\varpi }_a={\varpi }_a+{\varpi }_3$
18:	end if
19:	end if
20:	Update the operator selection probability based on the operator score, ${\varpi }_a=\{{\varpi }_a,a=1,2,3\}$
21:	until Algorithmic stopping criterion satisfied
22:	return $S_{best}$

3.2. Spatial and Temporal Distances

Our study tackles the challenge of balancing EV power constraints with customer time windows. When an EV is fully charged but cannot meet the time window requirements at subsequent customer nodes, it can trigger a vehicle route decomposition operation, leading to inefficient routing. We introduce a normalization method for spatial and temporal distances between customers to improve vehicle resource utilization. By applying weighted averages, we derive a combined spatiotemporal distance, enabling the construction of a high-quality initial solution for more efficient route optimization.

The time windows for customer nodes $i$ and $j$ in our model are defined as $[E{T}_i,L{T}_i]$ and $[E{T}_j,L{T}_j]$, respectively, with the assumption that $E{T}_i<E{T}_j$. Let ${\tau }_{if}^t\in [E{T}_i,L{T}_i]$ represent the moment when vehicle $f$ arrives at customer node $i$, and let node $i$’s service time be $s_i^t$. Therefore, the moment when vehicle $f$ arrives at customer node $j$ can be expressed as ${\tau }_{jf}^t\in [E{T}_i+s_i^t+t_{ij}^t,L{T}_i+s_i^t+t_{ij}^t]$, where $E{T}_i+s_i^t+t_{ij}^t=a$ and $L{T}_i+s_i^t+t_{ij}^t=b$. At this point, there will always be four possible distinctive cases, as illustrated in Figure 5.

• Case 1: When $b<E{T}_j$, vehicle $f$ arrives at customer node $j$ before the lower bound of the time window. In this case, the time distance between customer nodes $i$ and $j$ is defined as the minimum possible waiting time, represented as ${\delta }_1(E{T}_j-b)$, where ${\delta }_1$ is the waiting time penalty coefficient.

• Case 2: When $a>L{T}_j$, vehicle $f$ arrives at customer point $j$ after the upper limit of the time window. In this scenario, the time distance between customer nodes $i$ and $j$ is defined as the minimum possible delay time, represented as ${\delta }_2(a-L{T}_j)$, where ${\delta }_2$ is the delay time penalty coefficient.

• Case 3: When $a\le E{T}_j<b$ or $a<L{T}_j<b$, there exists an overlap between the moment vehicle $f$ arrives at customer node $j$ and the time window of customer node $j$. In this case, the time distance between customer nodes $i$ and $j$ is provided as ${\delta }_3(s_i^t+t_{ij})$, where ${\delta }_3$ is the deviation coefficient.

• Case 4: When $E{T}_j<a<b<L{T}_j$, the moment vehicle $f$ arrives at customer node $j$ precisely aligns with node $j$’s time window. In this case, the time distance between customer nodes $i$ and $j$ is given by $s_i^t+t_{ij}$.

The time distance between any two customer nodes $i$ and $j$ in the road network, denoted as $d_{ij}^T$ or $d_{ji}^T$, but $d_{ij}^T\ne d_{ji}^T$, is treated as non-directional during the algorithm’s calculation. As a result, the greater value between $d_{ij}^T$ and $d_{ji}^T$ is used as the time distance. Equation (54) provides the method for calculating this time distance.

(54)$d_{ij}^T=\left\{\begin{array}{l}{\delta }_1(E{T}_j-b),b<E{T}_j\\ {\delta }_2(a-L{T}_j),a>L{T}_j\\ {\delta }_3(s_i^t+t_{ij}),a\le E{T}_j<b\vee a<L{T}_j<b\\ s_i^t+t_{ij},E{T}_j<a<b<L{T}_j\end{array}\right.$

The spatial distance $d_{ij}^S$ between any two customer nodes $i$ and $j$ in the road network is defined using the Euclidean distance. Depending on the specific problem context, different weights are assigned to the temporal and spatial distances to calculate the combined spatiotemporal distance $d_{ij}^{ST}$, as shown in Equation (55).

(55)$d_{ij}^{ST}=\alpha d_{ij}^S+\beta d_{ij}^T,\alpha +\beta =1$

3.3. Solution Representation

The MD-MFVRPDTW problem is NP-hard, representing a typical combinatorial optimization challenge. In our study, the solution is encoded using integer permutation encoding. Suppose the road network includes $n$ customer and charging station nodes and $m$ distribution centers, represented by natural numbers $0,1,2,...,m$. The initial solution is described using both upper and lower codes: the upper code Ⓤ identifies the vehicle route, while the lower binary code Ⓛ, of the same length, determines the service type at each customer node. For a non-negative integer $i\in \{2,3,\dots ,n+1\}$, customer points are classified based on the following criteria.

• If the $i$-th node in Ⓛ is encoded as 0, then the corresponding customer node in Ⓤ will be served by the vehicle.

• If the $i$-th node in Ⓛ is encoded as 1, and the $\left(i-1\right)$-th node is encoded as 0, then the corresponding customer node in Ⓤ will be served by the drone.

• If the $i$-th node in Ⓛ is encoded as 1, and the $\left(i-1\right)$-th node is also encoded as 1, then the corresponding $i$-th customer node in Ⓤ will be served by the vehicle, the $\left(i-1\right)$-th customer node of Ⓤ will be served by the vehicle.

The encoded sequence for Ⓛ (with a minimum length of 1) can also be generated based on the three scenarios described earlier. In this sequence, a node encoded as “1” represents a drone customer served, while the drone’s launch and recovery nodes correspond to the Ⓤ sequence, specifically the immediately preceding and following “0s” in the sequence, respectively. Figure 6 provides an example of a routing strategy involving 2 distribution centers, 16 customer nodes, and 2 charging stations, along with the corresponding route encoding method.

3.4. Construction of Initial Solutions

• Customers are classified as vehicle-serviced or drone-serviced based on the technical capabilities of the drones. An EV carrying a drone is dispatched from the distribution center to serve the assigned customers. When conditions allow, the vehicle prioritizes customer deliveries with the shortest spatiotemporal distance from the previous node. If the EV cannot meet the customer’s capacity or time window requirements, it returns to the distribution center, following proximity-based principles. This process generates a new route for the EV and continues until all vehicle-serviced deliveries are completed.

• Evaluate the feasibility of the EV’s remaining power along the current route. If the route is deemed infeasible, implement the greedy charging station insertion strategy proposed by Keskin et al. [37]. During this process, if the EV’s charge is insufficient to reach the charging station or does not satisfy the time window constraints for upcoming customer nodes before the charging station, a fuel vehicle will be dispatched to complete the delivery.

• To simplify the algorithm code, dispatch a fuel vehicle from the distribution center, treating it as an EV with unlimited battery capacity. Designate the unassigned customer node as the first stop and construct a new route starting from the distribution center, prioritizing the minimization of spatiotemporal distance. After incorporating all customer nodes that were infeasible for the EV into the fuel vehicle’s route, perform route optimization. Set a threshold to merge customer nodes, eliminating any routes that do not meet this criterion. Then, a greedy insertion strategy can be used to place customer nodes from inefficient routes on EV or fuel vehicle routes. Finally, the objective function values of the original and newly optimized routes will be compared, keeping the solution more efficient.

• In the complete set of vehicle routes, drone customers are inserted last. The initial solution, constructed based on spatiotemporal distance, is represented by the 3D graph in Figure 7. This route includes five vehicle customers, two drone customers, and one charging station. The horizontal plane represents a two-dimensional space, while the vertical axis represents time. Cylinders represent customer time windows, and their height indicates the tolerance of each window. At this point, the system includes an EV spatiotemporal path 1, $O_1\to A_1\to A_2\to C_1\to C_2\to D_1\to D_2\to H_1\to H_2\to O_4$, a drone spatiotemporal path 1, $O_1\to B_1\to B_2\to C_2$, a fuel vehicle spatiotemporal path 2, $O_2\to F_1\to F_2\to G_1\to G_2\to O_3$, and a drone spatiotemporal path 2, $O_2\to E_1\to E_2\to F_2$. Both the vehicles and drones meet the time window requirements of their respective customers, with the distribution center’s time window being wider than those of all the customer nodes.

3.5. Neighborhood Search Approach

Based on the characteristics of the proposed MD-MFVRPDTW model, we develop four disruption operators (removal operators) and two repair operators (insertion operators), integrating perturbation techniques to enhance solution diversity during the neighborhood search process.

3.5.1. Removal Operators

• (1). Random Removal Operator: Randomly select several customer nodes from the current vehicle route, including all associated drone customers, for removal, and eliminate any resulting invalid routes. As shown in Figure 8, vehicle customer nodes 5, 10, 12, and 15, along with their respective drone customer nodes, are randomly removed from the route.

• (2). Extreme Cost Removal Operator: The extreme cost removal strategy eliminates all drone customers and removes invalid routes by calculating the potential savings for both types of vehicles in the mixed fleet. Customer nodes with higher fuel consumption savings are removed from fuel vehicle routes, while those with higher electric consumption savings are removed from EV routes. As shown in Figure 9, fuel consumption savings for each fuel vehicle customer node and electric consumption savings for each EV customer node are calculated. These savings values are then multiplied by a random factor between 0 and 1, sorted and ranked accordingly. Nodes 15 and 17 of the EV route and nodes 5 and 10 of the fuel vehicle route are selected for removal based on descending savings. After all drone customer points are removed, any newly generated invalid routes (e.g., $D_0\to 18\to D_0$) are also eliminated.

• (3). Extreme Distance Removal Operator: To optimize the routes, select the closest customer node to the first and last distribution centers in the fuel vehicle route and the farthest customer node from these centers in the EV route for removal. After extracting all drone customers, any invalid routes (if present) should also be eliminated. As illustrated in Figure 10, consider EV Route 2 and fuel vehicle Route 3 from the current solution. The distances $d_1$ and $d_2$ between each customer node and the first and last distribution centers are calculated. The customer node with the shortest distance (node 14) in the fuel vehicle route and the node with the maximum product of $d_1\times d_2$ (node 7) in the EV route are removed, followed by the extraction of all drone customer nodes.

• (4). Congestion Avoidance Operator: A congestion avoidance strategy is designed to adjust the delivery plan, reducing delays and excess fuel consumption during peak hours. Fuel vehicles are less efficient when idle or moving at low speeds, leading to increased fuel consumption and carbon emissions due to the time-varying conditions of the road network. To mitigate this, customer nodes with time windows that overlap with peak hours are removed from fuel vehicle routes. Similarly, customer nodes with time windows that intersect low-peak hours are removed from EV routes. Afterward, all drone customers are extracted to eliminate invalid routes. As shown in Figure 11, all vehicles on the network can experience both low-peak and peak traffic periods. In this scenario, customer nodes 15, 8, and 7 on EV Routes 2 and 4 are removed due to their overlap with low-peak hours, while customer node 14 on fuel vehicle Routes 1 and 3 is removed due to its overlap with peak hours. Finally, the invalid routes (e.g., $D_1\to 19\to D_1$) are deleted after extracting all the drone customer nodes.

3.5.2. Insertion Operators

• (1). Random Insertion Operator: For the customer nodes removed by the previously mentioned removal operators, a random insertion point is selected to reintroduce the vehicle customer nodes into the mixed fleet route. The drone customer nodes are then reinserted in the same manner as described in Section 3.4, repeating this process until all customer nodes have been reintegrated into the route. As illustrated in the route encoding in Figure 12, the removed vehicle customer nodes (5, 10, 12, and 15) are reinserted after their new insertion locations are chosen at random. When reintroducing drone customer nodes into the vehicle route, adjustments are made to the drone’s launch and retrieval nodes, generating a new solution.

• (2). Greedy Insertion Operator: This operator processes the removed vehicle and original drone customers after merging by inserting those that meet the maximum carrying capacity and flight distance requirements of the drone into the vehicle route with the smallest increase in cost, ensuring that all restrictions are met. In Figure 13, we first assess whether the removed customer nodes (5, 10, 12, 15, 4, 9, 6, 11, 13, and 16) meet the conditions for drone insertion. Original drone customer nodes 6 and 13, which no longer satisfy these conditions, are swapped with original vehicle customer nodes 10 and 12, which can now accommodate the drone. Then, customer nodes 5, 6, 13, and 15 are selected for insertion based on the smallest cost difference before and after insertion (each multiplied by a random number between 0 and 1), while also ensuring that the vehicle’s load capacity and customer time windows are respected. These nodes are reinserted into the vehicle route using the previously described construction method, resulting in a feasible solution after the drone customers are reintroduced.

3.6. Adaptive Adjustment Mechanisms and Solution Acceptance Criteria

The adaptive adjustment mechanism and solution acceptance criteria incorporate a random scoring system to strengthen the algorithm’s robustness. This approach enables the algorithm to adapt more effectively to the complexities of distribution scenarios involving the collaboration of a mixed fleet and drones, ensuring improved solution quality in dynamic and challenging environments. This adaptation increases the overall efficiency of algorithm execution. Detailed specifications of the adaptive adjustment mechanism and solution acceptance criteria are outlined below.

3.6.1. Choice of Operator

The selection of destroy and repair operators is carried out using the roulette wheel method. Let $Op$ represent the set of all operators, ${\omega }_i(i\in Op)$ denote the weight assigned to each operator, ${\delta }_i(i\in Op)$ represent the score attributed to each operator, and $v_i$ track the frequency of operator usage. During each iteration, operators are assigned random scores within three predefined scoring ranges based on the quality of the solutions they generate.

• To achieve the global optimal solution, add 1.4 to 1.6 points.

• To obtain a suboptimal solution, add 1.1 to 1.3 points.

• To achieve a feasible solution, add 0.7 to 0.9 points.

• Set the update coefficient $\xi$ of the weights between 0 and 1, then calculate the operator weights using Equation (56) based on the corresponding operator scores.

  (56)${\omega }_i=\left\{\begin{array}{l}{\omega }_i,\mathrm{if}{v}_i=0\\ (1-\xi )\ast {\omega }_i+\xi \ast {\displaystyle \frac{{\delta }_i}{v_i}},\mathrm{if}{v}_i>0\end{array}\right.,\forall i\in Op$

3.6.2. Stopping and Acceptance Criteria

• (1). The simulated annealing algorithm uses the Metropolis criterion to accept inferior solutions. To avoid premature convergence to a local optimum, the algorithm must accept a neighboring solution $f{(x)}_{new}$, even if it is worse than the current solution $f{(x)}_{old}$, with a certain probability $p$. At the beginning of the simulated annealing process, the initial temperature is set to $T_0$, and a decay rate between 0 and 1 is applied. The temperature $T$ for each iteration is computed as the product of the previous iteration’s temperature and the decay rate. The probability $p$ is calculated as shown in Equation (57).

  (57)$p=\left\{\begin{array}{l}1,\mathrm{if}f{(x)}_{new}<f{(x)}_{old}\\ \exp \left({\displaystyle \frac{f{(x)}_{old}-f{(x)}_{new}}{T}}\right),\mathrm{if}f{(x)}_{new}\ge f{(x)}_{old}\end{array}\right.$

• (2). The algorithm’s iteration process is terminated using a threshold for consecutive unimproved solutions. If the number of consecutive unimproved solutions, denoted as $num\_con$, reaches the predefined threshold $num\_end$, the algorithm stops and outputs the final satisfactory solution. The values for $num\_end$ and the maximum number of iterations $num\_iterations$ are set on a case-by-case basis, allowing the number of adaptive neighborhood searches to be customized for the specific problem being addressed.

4. Experimental Study

In this section, we validate the proposed ALNS-STD algorithm through numerical examples related to the multi-depot mixed fleet–drone collaborative delivery routing problem within the carbon trading framework, and we analyze the experimental results. All experiments were conducted on a computer with an Intel Core i7-10875H 2.30 GHz CPU, 16 GB RAM, running Windows 10, using Python 3.8 as the programming environment.

4.1. Parameters for Numerical Experiments

Since no existing algorithm is tailored to solve the proposed MD-MFVRPDTW model, this article aims to validate the effectiveness of the ALNS-STD algorithm. To do so, we adapted the instance C101_21 from the electric vehicle routing problem with time windows (EVRPTW), originally introduced by Schneider et al. [38]. As there is no standard benchmark for the MD-MFVRPDTW model with charging stations, we customized the C101_21 instance accordingly. First, we applied a K-means clustering algorithm using silhouette coefficients to group the nodes and determine both the number of distribution centers and the radial range of all nodes. Next, a simulated annealing algorithm was used to optimize the location of the distribution centers. This process transformed the EVRPTW problem into an MD-MFVRPDTW instance, labeled C101_21M, comprising 4 distribution centers, 21 charging stations, and 100 customer nodes.

4.1.1. Parameter Configuration for EV’s Nonlinear Charging Function

The nonlinear charging function for EVs is derived from the study by Montoya et al. [34]. The linear charging rate in the initial stage is set at 0.29 kWh/min, with the full charge capacity denoted as $B_t=77.75$. Using these parameters, the charging time $t_k$ and state of charge ${\alpha }_k$ at each inflection point $k$ can be determined. Specifically, the charging time $t_k$ and charge state ${\alpha }_k$ corresponding to the inflection point $k$ in the function are set as follows: $t_1,{\alpha }_1=[226.7814,65.3047],t_2,{\alpha }_2=[276.6302,73.4023],t_3,{\alpha }_3=[378.6145,77.75]$. Additionally, the state of charge ${\alpha }_k$ upon an EV’s arrival at a charging station and the charging time $t_k$ required for a full charge can be determined using Equation (58).

(58)${\mathrm{t}}_k=\left\{\begin{array}{l}378.6145-3.47{\alpha }_k,0\le {\alpha }_k<65.30\\ 554.1145-6.16{\alpha }_k,65.30\le {\alpha }_k<73.40\\ 1823.9945-23.46{\alpha }_k,73.40\le {\alpha }_k\le 77.75\end{array}\right.$

4.1.2. Parameters for Time-Dependent Vehicle Speeds

To model the time-dependent variations in vehicle speed using a step function, we randomly generate seven data points, each with an average speed of 56.3 km/h. These points represent the vehicle’s speed across different time intervals. The detailed time zones and their corresponding speed values are provided in Table 5.

To examine the effects of speed fluctuations on distribution path planning under time-varying conditions, this study defines three distinct vehicle speed scenarios and performs a comparative analysis of their respective distribution outcomes. Specifically, Case 1 assumes a constant speed of 56.3 km/h. Case 2 reflects time-varying speeds, as outlined in Table 6. Case 3 simulates a highly congested scenario typical of extended holiday periods, in which speed values from Case 2 are reduced by half. The specific time zones and corresponding speed values for each scenario are presented in Table 6.

4.1.3. Parameters for Power Consumption

We introduce a mixed fleet of vehicles with comparable load capacities, including a commonly used cargo van such as the Changan Star 5, which has a payload capacity of 450 kg, serving as the fuel delivery vehicle. The parameters are summarized in Table 7.

4.1.4. Initial Values for Remaining Experiment-Related Parameters

The initial values for the remaining experiment parameters in this section are provided in Table 8.

4.2. Analysis of Experimental Results

The experiment tested three algorithms on the C101_21M case, conducting 10 runs for each. These algorithms included an improved ALNS algorithm incorporating an EV module (ALNS-1), an improved ALNS algorithm utilizing the neighborhood operator from this study without considering spatiotemporal distance (ALNS-2), and an improved ALNS algorithm that integrates spatiotemporal distance (ALNS-STD). Each algorithm was executed 10 times to generate solutions. In this analysis, $f_1$ and $f_2$ represent the number of dispatched EVs and fuel vehicles, respectively, $Z$ refers to the best solution obtained by the algorithm, while $T$ and $I$ indicate the solution time and the optimal number of iterations, respectively. The computational results are presented in Table 9.

The results in Table 9 demonstrate the performance of the three algorithms, each applied 10 times to the C101_21M instance. The ALNS-1 algorithm, which does not account for the equilibrium between total delivery costs and carbon trading costs, prioritizes dispatching fuel vehicles to minimize EV transportation costs. This approach achieves an optimal objective function value of 1548.98. In contrast, the ALNS-2 algorithm incorporates both fuel and electricity consumption in optimizing the neighborhood operator, making it more responsive to cost variations and better balancing the dispatch of fuel and electric vehicles. This results in an optimal objective function value of 1615.24. The ALNS-STD algorithm goes further by factoring in both spatial and temporal customer distributions to construct a better initial solution, achieving improved transportation costs through a balanced dispatch of fuel and electric vehicles. The ALNS-STD algorithm produces objective function values ranging from a maximum of 1628.22 to a minimum of 1562.86, with deviations from the median at 2.39% and −1.72% and from the mean at 2.07% and −2.03%. These deviations are within acceptable limits. Additionally, the minimum objective function value of ALNS-STD represents improvements of −0.9% and 3.24% compared to ALNS-1 and ALNS-2, respectively.

Although the ALNS-STD algorithm’s optimal objective function value is slightly higher than that of ALNS-1, its integration of low-carbon considerations makes it a more practical and theoretically sound solution for logistics and distribution companies pursuing low-carbon transitions. Regarding iterations, the ALNS-STD algorithm consistently converges to a satisfactory solution within 200 generations. The convergence curve in Figure 14 demonstrates the algorithm’s stability, coinciding with the known optimal solution of 1562.86 after 200 iterations. In summary, the proposed ALNS-STD algorithm effectively solves the MD-MFVRPDTW model by incorporating a carbon trading mechanism. Figure 15 further visualizes the solution’s satisfactory results, confirming the algorithm’s applicability and effectiveness.

4.3. Distribution System Composition Analysis of the Composition of the Distribution System

4.3.1. Impact of Drones on Delivery System Performance

This experiment establishes a mixed fleet delivery system without incorporating drone customer nodes by using the C101_21M case as the basis ($Q_d\leftarrow 1<q_i\Rightarrow {V_D}^{\prime }=0$) and examines the impact of drone integration on delivery performance by comparing the solution results of the MD-MFVRPDTW and MD-MFVRPTW models. The number of customers satisfying the drone technical features constraints is set to $V_D=25$, $V_D=50$, and $V_D=75$, while other parameters remain at their default values. In this context, ${V_D}^{\prime }$ denotes the actual number of customers served by drones, $V_T$ represents the number of customers served by vehicles, while $D_d$ and $D_t$ represent the total mileage covered by drones and vehicles, respectively. Furthermore, ERC indicates the electric consumption and charging costs of EVs, while FRC represents the fuel consumption costs of fuel vehicles. The ALNS-STD algorithm is applied to solve both models, with the results summarized in Table 10 and the comparison illustrated in Figure 16.

The experimental results demonstrate that the MD-MFVRPDTW model offers a more cost-effective solution by reducing total delivery costs through drones to decrease vehicle travel distances. As the number of customers meeting the technical constraints for drone operations increases, the disparity in the total vehicle travel distance between the two delivery models grows more pronounced. Notably, at ${V_D}^{\prime }=24$, the MD-MFVRPDTW model achieves lower EV energy consumption and charging costs compared to the MD-MFVRPTW model. Although the contribution of drones to reducing fuel vehicle consumption was less than expected, the overall delivery costs in the MD-MFVRPDTW model were substantially optimized due to the integration of the drone. The improvement in total delivery costs is primarily attributed to the fact that in the MD-MFVRPTW model, vehicles visit more customer nodes and perform route-splitting operations more frequently, leading to higher fixed vehicle dispatching costs. In conclusion, incorporating drones into logistics and distribution systems effectively reduces delivery costs, although their efficiency in mixed fleets largely depends on the proportion of customers they can serve.

4.3.2. Impact of a Single Vehicle Type on Delivery Costs

The experiments generate two distinct delivery systems: the MD-VRPDTW system, consisting solely of fuel vehicles, and the MD-EVRPDTW system, comprising only EVs. These systems are used to analyze the impact of fuel vehicles and EVs on delivery costs by comparing their solution results with those of the MD-MFVRPDTW system. The results of the algorithm runs are presented in Table 11, with the MD-MFVRPDTW data sourced from Table 10. A comparative analysis of the three models is illustrated in Figure 17.

The experimental results reveal that in terms of total vehicle distance traveled, the MD-EVRPDTW model consistently records higher values than the MD-MFVRPDTW model, which in turn exceeds those of the MD-VRPDTW model. This disparity arises because EVs are limited by battery capacity and must often take detours to visit charging stations, thereby increasing the total vehicle miles traveled. For drone mileage, the MD-VRPDTW model shows higher values than the MD-MFVRPDTW model, which in turn exceeds the MD-EVRPDTW model. However, the number of customer nodes visited by drones remains similar or identical across all three models. This outcome can be attributed to two factors: (i) the models exhibit a uniform tendency in selecting drones for deliveries; (ii) the use of charging stations as dual-purpose nodes for drone takeoff and landing reduces drone flight mileage while simultaneously expanding the solution space for the algorithm. These findings highlight the trade-offs and advantages inherent in different delivery systems.

The total delivery cost of the MD-MFVRPDTW model is higher than that of the MD-VRPDTW model. However, its fuel consumption cost is consistently lower. This finding offers valuable insights: under the carbon trading mechanism, an essential tool for advancing low-carbon initiatives, EVs cannot fully replace fuel vehicles in logistics distribution systems in the short term. Nevertheless, the energy-saving technology of EVs is poised to become a critical competitive advantage, driving the logistics market and playing a pivotal role in the transition toward greener commercial logistics.

Similarly, the total delivery cost of the MD-EVRPDTW model exceeds that of the MD-MFVRPDTW model. While the MD-EVRPDTW model reduces fuel consumption costs, the increasing number of EVs leads to a rapid rise in electricity consumption and charging costs. Additionally, fluctuations in the number of customers served by drones significantly influence EV electricity consumption, charging costs, and overall delivery costs. These results highlight the cost-reduction potential of drones and underscore that, given current technology, EVs alone cannot entirely replace fuel vehicles in logistics. A mixed fleet effectively reconciles the trade-off between economic efficiency and environmental sustainability in the logistics and distribution industry.

4.4. Sensitivity Analysis of the Mixed Fleet Parameters

4.4.1. Impact of EV Battery Capacity on Delivery System Performance

In this experiment, we analyze the effect of EV battery capacity on delivery route planning by solving the MD-MFVRPDTW problem across various levels of battery capacity. The EV battery capacity, $E_t$, varied from 60 kWh to 100 kWh in increments of 10 kWh, with all other parameters held constant. The algorithm performance results under these different battery capacities are summarized in Table 12.

The experimental results demonstrate that as the battery capacity of EVs in the hybrid fleet increases from 60 kWh to 100 kWh, both total delivery costs and electricity consumption, including charging costs, progressively increase across various delivery routes. Conversely, fuel consumption costs and total vehicle mileage steadily decline. This trend occurs because logistics EVs with smaller battery capacities frequently detour to charging stations within a given distribution range, thereby increasing total vehicle mileage. As battery capacity increases, the need for detours decreases, and vehicle mileage stabilizes.

However, in the short term, a larger battery capacity leads to a higher number of EVs being dispatched within the mixed fleet. The associated rise in fixed costs outweighs the savings from reduced detour mileage, resulting in an overall increase in total delivery costs. Specifically, the cost of electricity consumption and charging increases by 41.91%, while fuel consumption costs decrease by 86.10% and total vehicle mileage reduces by 6.97%. Despite these changes, the total delivery cost still rises by 4.69%.

From a fleet composition perspective, the fixed cost of EVs in this problem is higher than that of fuel vehicles. Reducing the ratio of EVs to fuel vehicles theoretically lowers total delivery costs, a trend supported by the experimental results. Conversely, increasing the EV-to-fuel vehicle ratio results in poorer cost optimization. Figure 18 illustrates that variations in battery capacity input parameters lead to differences in delivery routes and costs. This demonstrates the generality of the proposed model and its ability to maintain robust performance under changing input conditions.

In summary, variations in EV battery capacity within the mixed fleet have a substantial impact on delivery costs and distribution routing. Achieving more optimal distribution solutions will require advancements in battery technology, enhancements in supporting infrastructure, and reductions in EV fixed costs. This highlights the current challenges faced by logistics companies, emphasizing that EVs alone are not yet capable of fully replacing fuel vehicles to achieve a complete energy transition at this stage.

4.4.2. Impact of Fuel Consumption Rate on Delivery System Performance

The fuel consumption rate $G$ of a fuel vehicle during operation is influenced by the actual mechanical power $p$, which is, in turn, affected by the vehicle’s time-dependent travel speeds. To better understand the impact of varying fuel consumption rates on delivery route planning, this experiment defines three fuel consumption scenarios: Case 1, Case 2, and Case 3, while keeping all other parameters constant. Secondly, to evaluate the model’s robustness against inaccuracies in input parameters, the optimal delivery route derived from Case 1 is applied to the conditions of Case 2, creating a new scenario designated as Case 4. The performance results of the algorithm on these three fuel consumption levels are presented in Table 13.

The experimental results show that as the fuel consumption rate of fuel vehicles in the mixed fleet increases from Case 1 to Case 2, the electricity consumption and charging costs for EVs increase by 5.05%, the fuel consumption costs increase by 12.02%, and the total delivery cost increases by 1.73%. When the fuel consumption rate changes from Case 1 to Case 3, electricity consumption and charging costs for EVs increase by 17.42%, while fuel consumption costs decrease by 8.30%, resulting in a 4.04% increase in total delivery costs. The optimal cost in Case 4 exceeds that of Case 1, highlighting the model’s sensitivity to input parameter accuracy. This underscores the importance of precise parameter inputs, as inaccuracies can result in suboptimal solutions. These findings suggest that both the costs related to EVs, and overall delivery costs increase as the fuel consumption rate for fuel vehicles increases.

However, beyond Case 2, fuel consumption costs decrease as fewer fuel vehicles are dispatched. Additionally, the flight mileage of the drone in this experiment behaves differently from the previous experiment. In the earlier experiment, drone mileage remained relatively stable across mixed fleet configurations, influenced by EV battery capacity and charging station use. However, in this experiment, the fuel consumption rate primarily affects the number of EV dispatches, and fuel vehicles become more advantageous at lower consumption rates. As a result, the algorithm reduces the reliance on EVs and charging stations, leading to an increase in drone mileage across the delivery network.

In summary, variations in the fuel consumption rate of fuel vehicles within a mixed fleet have a significant influence on delivery costs and route planning. Furthermore, advancements in fuel-efficient and emission-reduction technologies for fuel vehicles could enhance the economic efficiency of mixed fleet delivery systems.

4.5. Analysis of the Geographical Distribution of Charging Stations

Schneider et al. [38] proposed a set of EVRPTW algorithms that include three types of instances classified based on the geographical distribution of nodes: cluster distribution (Cluster, C), random distribution (Random, R), and random-cluster distribution (Random-Cluster, RC). In the C and RC instances, the charging station coordinates are identical. In this experiment, we use the charging station data from the C and R instances of Schneider et al. [38] and independently design a new instance with a central-cluster distribution of charging stations (Central-Cluster, CC). In the CC instance, the customer nodes follow the same distribution pattern as in C101_21M. This approach allows us to analyze the impact of different geographical distributions of charging stations on overall distribution performance. The number of charging stations used in each case is indicated, and the path visualizations for the three charging station distributions are shown in Figure 19. The corresponding algorithm performance results are presented in Table 14.

The experimental results reveal significant differences in solution outcomes based on the geographic distribution of charging stations—clustered, random, and centralized clustering. These differences affect fleet composition, charging station utilization, and delivery costs in distinct ways:

• In clustered distribution, the ratio of EVs to fuel vehicles is the lowest, charging station utilization is minimal, and fuel vehicles contribute more significantly to fuel consumption costs. This occurs because the limited coverage area of clustered charging stations restricts EV flexibility in recharging, reducing their operational range and utility. Consequently, fewer EVs are dispatched, lowering fixed costs. EVs that do undertake long-distance deliveries require a full or nearly full charge to complete tasks, leading to extended charging times and higher energy levels per session. This increases the cost per charging session and overall energy costs, especially for high-load or long-distance deliveries. Fuel vehicles, unaffected by charging station locations, are used more frequently, further reducing the reliance on EVs.

• In random distribution, the randomly distributed charging stations allow EVs to effectively utilize the nearest available station, minimizing energy waste and empty mileage. This reduces energy consumption and improves efficiency. The distribution of charging demands enhances station utilization rates and lowers unit charging costs. This prompts the distribution center to dispatch more EVs, increasing the mixed fleet size and fixed costs.

• In centralized clustering distribution, charging stations are readily accessible, enabling EVs to recharge quickly after completing deliveries and minimizing idle time. This facilitates frequent operations and increases the ratio of EVs to fuel vehicles. Higher charging station utilization rates optimize fuel consumption costs and support EV delivery activities. While fixed costs increase due to greater reliance on EVs, the system achieves improved efficiency and reduced fuel vehicle dependency. However, when the sole optimization objective is to minimize total delivery costs, the clustered distribution in the original instance remains more economically favorable.

In short, each charging station distribution impacts fleet composition and costs differently. Clustered distributions favor fuel vehicles and lower fixed costs, random distributions balance costs, and fleet size, while centralized clustering supports higher EV utilization but with increased fixed costs. The choice of distribution should align with operational priorities, such as cost minimization or the transition to greener logistics.

4.6. Analysis of Carbon Trading Mechanisms

4.6.1. Impact of Carbon Pricing on Delivery Route Optimization

To evaluate the impact of carbon pricing within the carbon trading framework, the experiment gradually increases the carbon price from 0 CNY/kg to 1 CNY/kg in increments of 0.1 CNY/kg. The carbon allowance of the company is fixed at 150 kg, while all other parameters remain at their default values. The results of the solution generated by the proposed algorithm are presented in Table 15. In this context, $DRC$ represents drone transportation costs, $CTC$ denotes carbon trading costs, and $CE$ refers to carbon emissions.

The experimental results show that as the carbon price increases from 0 CNY/kg to 1 CNY/kg, the power consumption and charging costs of EVs increase steadily, while the fuel consumption costs of fuel vehicles decrease. Additionally, drone transportation costs also tend to decline. This can be explained as follows.

• Under the carbon trading mechanism, the rising carbon price incentivizes the logistics system to reduce the number of fuel vehicles dispatched, resulting in increased EV mileage.

• The increased use of EVs and their longer travel distances lead to higher utilization of charging stations and more frequent selection of drone recovery nodes, which in turn shortens drone flight distances and reduces transportation costs.

With regard to carbon emissions and carbon trading costs, the system’s carbon allowance exceeds the requirements for daily operations as a result of the reliance on fuel vehicles. As carbon emissions drop from 25.08 kg to 17.97 kg, the revenue from selling excess carbon credits increases. However, it is important to note that total delivery costs show an upward trend, rising from 1419.62 to 1718.62, which represents a 21.06% increase. This indicates that the high fixed costs of EVs, along with electricity expenses, significantly affect delivery route planning, outweighing the impact of carbon pricing. Only when the carbon price crosses a certain threshold can the logistics system’s route planning under the carbon trading mechanism achieve both carbon reduction and cost savings.

4.6.2. Impact of Carbon Allowance on Delivery Route Optimization

To examine the impact of carbon allowances on delivery route planning within the carbon trading mechanism, we adjusted the company carbon allowance ($\overline{CE}$) to increase in increments of 50 kg, ranging from 0 kg to 300 kg. Meanwhile, the carbon price ($cp$) was held constant at 0.5 CNY/kg, and all other parameters were kept at their default values. The solution results of the experiment are presented in Table 16.

The experimental results presented in Table 16 reveal that as the company carbon allowance increases from 0 kg to 300 kg and the carbon trading cost shifts from 10.80 to −139.19, there are only minor fluctuations in the power consumption and charging costs of EVs, fuel consumption costs of fuel vehicles, drone transportation costs, and the carbon emissions associated with the solutions. Notably, only the total delivery cost exhibits a downward trend, decreasing from 1695.77 to 1503.59—a reduction of 12.78%. This is because the relaxation of carbon quotas impacts only the deployment and utilization of fuel vehicles by companies, without influencing other cost components. In this logistics and delivery system, the company incurs additional expenses for purchasing carbon allowances only when the carbon allowance is less than 50 kg. Since the carbon allowance is directly related to the carbon trading costs in the objective function, variations in the carbon allowance theoretically have minimal impact on the optimization of the delivery routes.

5. Conclusions

This paper presents a multi-depot vehicle–drone collaborative delivery routing optimization problem that integrates mixed fleets and drones within the framework of a carbon trading mechanism. To minimize total costs—including those associated with electric vehicle (EV) transportation, EV charging, fuel vehicle transportation, drone operations, carbon trading, and time window penalties—we develop a mixed-integer programming model, MD-MFVRPDTW. This model captures critical factors such as nonlinear charging times, time-varying vehicle travel speeds, fuel consumption, and mechanical power dynamics.

To solve this complex problem, we propose an adaptive large neighborhood search algorithm that integrates spatiotemporal distance (ALNS-STD). The algorithm includes four types of destroy and two types of repair operators for targeted improvements, along with an adaptive weight adjustment mechanism, which enhances solution efficiency and adaptability. This approach not only generates high-quality initial solutions but also optimizes routes by accounting for spatial and temporal distances. We summarize the conclusions of this study as follows.

• (1). The proposed model and algorithm were validated through experimental studies, which highlight the significant influence of various factors, such as fleet composition, technical parameters, geographic distribution of charging stations, and the carbon trading mechanism, on economic efficiency and route optimization. These findings underscore the practical value of the model in enhancing low-carbon, cost-effective logistics operations, offering actionable insights for logistics companies seeking to balance environmental sustainability with operational efficiency.

• (2). Traditional VRP studies often rely solely on fuel-powered vehicles for delivery tasks, contributing significantly to air pollution. In contrast, the model developed in this study leverages drones, which, due to their technical advantages, can handle a portion of delivery tasks. This reduces dependence on fuel-powered vehicles, thereby lowering greenhouse gas emissions and supporting environmental protection. Furthermore, the integration of electric logistics vehicles enhances carbon emission reductions and contributes to improved urban air quality.

• (3). Amid the low-carbon transformation and the rise of the digital economy, our proposed model extends the traditional VRP problem to address the challenges of last-mile logistics. It effectively reduces delivery costs, enhances operational efficiency, and promotes the adoption of green, environmentally friendly delivery practices. This model offers a more efficient and sustainable solution for the e-commerce and online retail industries, while providing logistics companies with a powerful tool to strengthen their market competitiveness.

• (4). In addition to providing a theoretical foundation for logistics companies to make informed, profit-driven decisions, this research also offers guidance for policymakers in regulating and shaping market behavior. It contributes to the growing body of knowledge on vehicle routing problems and opens new avenues for future research.

There are several areas for improvement in the current research. Our model is developed under idealized conditions and does not fully account for the complexities of real-world scenarios, such as multiple constraints, dynamic environments, power energy management, and the effects of charging methods on battery capacity degradation. The simulation of real-world conditions remains limited. Future research could benefit from incorporating more advanced features, such as heterogeneous fleets and electric vehicle charging and swapping strategies, to enhance the model’s practicality and realism.

Furthermore, while the ALNS-STD algorithm improves the breadth and diversity of the search and enhances global optimization capabilities, it still faces challenges, including high uncertainty, significant computational costs, and slow convergence rates. To address these issues, it may be necessary to fine-tune the search strategy for specific applications or integrate the strengths of other algorithms. Developing more efficient solution methods tailored to the VRP and its extended problems could significantly improve the algorithm’s performance and applicability.

Author Contributions

Conceptualization, Y.P.; methodology, Y.P., Y.Z. (Yanlong Zhang), D.Z.Y., S.L. and Y.Z. (Yali Zhang); formal analysis, Y.P., Y.Z. (Yanlong Zhang) and D.Z.Y.; data curation, Y.P., Y.Z. (Yanlong Zhang), S.L. and Y.Z. (Yali Zhang); writing—original draft, Y.Z. (Yanlong Zhang) and S.L.; writing—review and editing, D.Z.Y. and Y.S.; visualization, Y.Z. (Yanlong Zhang), D.Z.Y. and Y.S.; supervision, D.Z.Y.; funding acquisition, Y.P. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data are available upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. An example of the multi-depot mixed fleet–drone collaborative delivery problem.

Figure 2. Nonlinear charging function.

Figure 3. Scenario 1 of the fuel vehicle’s fuel consumption calculation.

Figure 4. Scenario 2 of the fuel vehicle’s fuel consumption calculation.

Figure 5. Four cases of time window constraints.

Figure 6. An example routing plan of the MD-MFVRPDTW and its coding representation.

Figure 7. Illustration of a complete set of spatiotemporal routes.

Figure 8. Illustration of the random removal operator.

Figure 9. Illustration of the extreme cost removal operator.

Figure 10. Illustration of the extreme distance removal operator.

Figure 11. Illustration of the congestion avoidance operator.

Figure 12. Illustration of the random insertion operator.

Figure 13. Illustration of the greedy insertion operator.

Figure 14. Convergence curve of the algorithm iterations.

Figure 15. Route diagrams of satisfactory solutions for C101_21M.

Figure 16. Comparative analysis results of MD-MFVRPDTW and MD-MFVRPTW.

Figure 17. Comparative analysis results of MD-MFVRPDTW, MD-VRPDTW, and MD-EVRPDTW.

Figure 18. Delivery routes with different battery capacities.

Figure 19. Route diagrams for the three charging station distributions. (a) C distribution. (b) R distribution. (c) CC distribution.

Table 1

Set parameters.

Symbol	Definition
$m$	Number of distribution centers
$n$	Number of customer nodes
$r$	Number of charging stations
$N_0$	Set of departure distribution centers, $N_0=\{0,1,\dots ,m\}$
$N_C$	Set of customer nodes, $N_C=\{m+1,\dots ,m+n\}$
$N_S$	Set of charging stations, $N_S=\{m+n+1,\dots ,m+n+r\}$
$N_0^{+}$	Set of return distribution centers, $N_0^{+}=\{m+n+r+1,\dots ,2m+n+r+1\}$
$N$	Set of all nodes, $N=N_0\cup N_S\cup N_C\cup N_0^{+}$
$N_V$	Set of all nodes except distribution centers, $N_V=N_S\cup N_C$
$N_D$	Set of drone customer nodes, $N_D\subset N_C$
$N_L$	Set of drone launch nodes, $N_L=N_0\cup N_S\cup N_C$
$N_R$	Set of drone recovery nodes, $N_R=N_S\cup N_C\cup N_0^{+}$
$A$	Set of all arcs through which vehicles travel, $A=\{(i,j)\left|i,j\in N\right.\}$
$R_r$	Set of vehicle types: Let $r\in R_r$ represent any vehicle type. Specifically, let 1 denote an electric vehicle and 2 denote a fuel vehicle.
$F_r$	A fleet of $r$-type vehicles: Any $r$-type vehicle is associated with a drone. Let $f_r\in F_r$ denote the f-th $r$-type vehicle in the fleet.

Table 2

Model parameters.

Symbol	Definition
$d_{ij}^t$	Vehicle’s travel distance from node $i$ to node $j$
$d_{ij}^d$	Drone’s flight distance from node $i$ to node $j$
$t_{ij}^t$	Vehicle’s travel time from node $i$ to node $j$
$t_{ij}^d$	Drone’s flight time from node $i$ to node $j$
$s_i^t$	Vehicle service time at node $i$
$s_i^d$	Drone service time at node $i$
$[E{T}_i,L{T}_i]$	Service time window for customer node $i$
${\delta }_1$	Waiting time penalty cost coefficient
${\delta }_2$	Delay time penalty cost coefficient
$q_i$	Demand at node $i$
$l_{ij}$	Fuel vehicle’s load from node $i$ to node $j$
$Q_r$	$r$-Type vehicle’s payload capacity
$Q_d$	Drone payload capacity
$B_t$	Battery capacity of an electric vehicle
$B_d$	Maximum drone endurance
$v_d$	Drone flight speed
$h$	Energy consumption coefficient per kilometer for electric vehicles
$lt$	Drone launch operation time
$rt$	Drone recovery operation time
$\theta$	Fuel vehicle’s carbon emissions coefficient
$c_1$	Electric vehicle’s unit transportation cost
$c_2$	Charging station’s unit charging cost
$c_3$	Fuel vehicle’s unit fuel price
$c_4$	Drone’s unit transportation cost
$c_{f_r}$	$r$-Type vehicle’s fixed cost
$cp$	Carbon trading price
$\overline{CE}$	Corporate carbon allowances
$U_{i{f}_1}$	Electric vehicle $f_1$’s piecewise function charging time at charging station $i$
$M$	A large positive number

Table 3

Decision variables.

Symbol	Definition
$x_{ij{f}_r}$	Binary variable indicates whether vehicle $f_r$ travels from node $i$ to node $j$ or not
$x_{ijk{f}_r}$	Binary variable indicates whether drone $f_r$ takes route $<i,j,k>$ or not
$y_{i{f}_r}$	Binary variable indicates whether vehicle $f_r$ chooses to charge at charging station $i$ or not
$p_{ij{f}_r}$	Binary variable indicates whether node $i$ is visited by vehicle $f_r$ before node $j$ or not

Table 4

Non-decision variables.

Symbol	Definition
$a_{i{f}_r}^t$	Vehicle $f_r$’s arrival time at node $i$
$a_{i{f}_r}^d$	Drone $f_r$’s arrival time at node $i$
$e_{i{f}_1}^{(1)}$	EV $f_1$’s power level when arriving at node $i$
$e_{i{f}_1}^{(2)}$	EV $f_1$’s power level when leaving node $i$
$e{c}_{i{f}_1}$	EV $f_1$’s charging amount at node $i$
$t{c}_{i{f}_1}$	EV $f_1$’s charging time at node $i$
$u_{i{f}_r}$	Vehicle $f_r$’s load level when visiting node $i$

Table 5

Vehicle travel speeds across different time zones.

Serial Number	Time Zone	Speed (km/h)
1	7:00–9:30	43.4
2	9:30–13:30	64.7
3	13:30–15:30	68.1
4	15:30–17:30	52.4
5	17:30–20:30	33.8
6	20:30–24:00	56.3
7	24:00–31:00	75.4

Table 6

Vehicle speed variations in different scenarios.

Scenario	Speed (km/h)
	7:00–9:30	9:30–13:30	13:30–15:30	15:30–17:30	17:30–20:30	20:30–24:00	24:00–31:00
Case 1	56.3	56.3	56.3	56.3	56.3	56.3	56.3
Case 2	43.4	64.7	68.1	52.4	33.8	56.3	75.4
Case 3	21.7	32.25	34.05	26.2	16.9	28.15	37.7

Table 7

Parameters related to fuel vehicle power consumption.

Symbol	Definition	Value	Symbol	Definition	Value
$\epsilon$	Engine friction coefficient (kJ/r/L)	0.2	$\xi$	Fuel-to-air ratio	0.07
$\varphi$	Engine speed (r/s)	105	k	Calorific value of gasoline (kg/g)	43
$A_f$	Fuel vehicle windward area (m2)	2.7	φ	Fuel conversion factor (L/g)	737
$m_f$	Empty fuel vehicle weight (kg)	1125	$\rho$	Air density (kg/m3)	1.2041
$\varpi$	Engine efficiency parameter	0.35	$g$	Gravitational constant (m/s2)	9.81
$c_r$	Rolling resistance coefficient	0.01	$c_d$	Air resistance coefficient	0.7

Table 8

Initial parameter values related to experiments.

Parameter	Value	Parameter	Value
EV payload (kg)	200	Gasoline carbon emissions coefficient	2.356
Fuel vehicle payload (kg)	450	Unit carbon price (CNY)	0.5
EV battery capacity (kWh)	77.75	Gasoline price (CNY)	7.88
Drone payload (kg)	10	Corporate carbon allowances (kg)	150
Drone endurance (min)	55	Initial temperature of simulated annealing	100
Drone flight speed (km/h)	80	Decay rate of simulated annealing	0.94
Drone launch and recovery time (min)	1	Initial weights of operators	0
EV power consumption per km	1	Operator’s initial score	1
Electricity cost per km (CNY)	0.6	Operator weight update factor	0.3
EV’s unit charging cost (CNY)	0.4	Threshold number of unimproved solutions	200
EV’s unit fixed costs (CNY)	100	Maximum number of iterations	10,000
Fuel vehicle’s unit fixed cost (CNY)	90
Drone’s cost per unit distance (CNY)	0.3

Table 9

Computational results for the C101_21M instance.

No.	ALNS-1				ALNS-2				ALNS-STD
	${\mathit{f}}_1/{\mathit{f}}_2$	$\mathit{Z}$	$\mathit{T}$	$\mathit{I}$	${\mathit{f}}_1/{\mathit{f}}_2$	$\mathit{Z}$	$\mathit{T}$	$\mathit{I}$	${\mathit{f}}_1/{\mathit{f}}_2$	$\mathit{Z}$	$\mathit{T}$	$\mathit{I}$
1	6/8	1563.81	128.26	226	7/6	1631.60	120.07	208	6/6	1585.34	122.60	192
2	5/7	1579.45	131.78	230	7/7	1619.92	128.30	225	7/6	1605.71	121.64	188
3	6/6	1574.06	137.90	238	7/7	1625.25	122.79	215	7/6	1619.68	118.77	184
4	6/7	1581.37	127.92	217	6/7	1615.24	137.44	230	8/6	1628.22	115.99	179
5	5/8	1560.80	136.60	221	7/7	1627.71	120.01	212	6/6	1570.30	125.03	196
6	6/6	1565.33	150.85	247	7/6	1638.10	122.37	192	6/6	1585.60	122.87	190
7	5/8	1558.71	135.34	236	7/6	1635.47	118.64	198	6/6	1594.78	120.91	188
8	5/7	1562.07	137.29	229	6/6	1628.31	127.65	216	7/6	1618.97	117.81	181
9	4/9	1548.98	145.15	242	6/7	1630.06	129.70	195	5/7	1562.86	129.01	198
10	7/7	1574.10	125.91	220	6/7	1629.88	125.17	210	6/7	1581.25	122.30	195

Table 10

Solution results of MD-MFVRPDTW and MD-MFVRPTW examples.

Problem	${\mathit{V}}_{\mathit{D}}$	${\mathit{V}}_{\mathit{D}}^{’}$	${\mathit{V}}_{\mathit{T}}$	${\mathit{f}}_1/{\mathit{f}}_2$	${\mathit{D}}_{\mathit{d}}$	${\mathit{D}}_{\mathit{t}}$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{Z}$
MD-MFVRPDTW	25	23	77	6/6	247.73	840.29	397.55	76.30	1624.58
	50	24	76	6/6	264.85	782.07	359.61	74.62	1589.84
	75	27	73	7/6	325.25	760.17	226.33	66.98	1565.90
MD-MFVRPTW	25	0	100	7/8	0	851.91	389.08	74.31	1819.50
	50	0	100	8/7	0	848.60	402.75	61.07	1827.95
	75	0	100	7/8	0	837.51	387.42	68.38	1811.02

Table 11

Solution results of MD-MFVRPDTW, MD-VRPDTW, and MD-EVRPDTW examples.

Problem	${\mathit{V}}_{\mathit{D}}$	${\mathit{V}}_{\mathit{D}}^{’}$	${\mathit{V}}_{\mathit{T}}$	${\mathit{D}}_{\mathit{d}}$	${\mathit{D}}_{\mathit{t}}$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{Z}$
MD-MFVRPDTW	25	23	77	247.73	840.29	397.55	76.30	1624.58
	50	24	76	264.85	782.07	359.61	74.62	1589.84
	75	27	73	325.25	760.17	226.33	66.98	1565.90
MD-VRPDTW	25	23	77	252.05	812.02	0	151.42	1524.67
	50	24	76	297.42	767.98	0	124.35	1507.17
	75	27	73	354.64	731.73	0	106.45	1503.76
MD-EVRPDTW	25	23	77	229.62	893.48	417.41	0	1683.76
	50	24	76	257.74	873.84	385.96	0	1669.76
	75	27	73	291.65	830.25	306.78	0	1635.81

Table 12

Solutions results of MD-MFVRPDTW with various EV battery capacity levels.

No.	${\mathit{E}}_{\mathit{t}}$	${\mathit{f}}_1/{\mathit{f}}_2$	${\mathit{D}}_{\mathit{d}}$	${\mathit{D}}_{\mathit{t}}$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{Z}$
1	60	4/8	275.21	782.02	300.77	122.72	1512.19
2	70	5/7	273.88	748.37	328.13	92.72	1567.41
3	80	6/6	275.66	741.59	351.80	72.83	1569.40
4	90	7/5	271.31	730.45	364.55	40.42	1571.87
5	100	7/4	279.20	727.50	426.82	17.06	1583.22

Table 13

Solutions results of MD-MFVRPDTW with various fuel consumption levels.

	$\mathit{G}$	${\mathit{f}}_1/{\mathit{f}}_2$	${\mathit{D}}_{\mathit{d}}$	${\mathit{D}}_{\mathit{t}}$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{Z}$
Case 1	0.67	4/8	337.60	744.84	330.73	58.47	1544.22
Case 2	1.33	6/6	277.59	787.62	347.44	65.50	1571.01
Case 3	2.66	7/5	272.04	792.31	388.35	53.62	1606.60
Case 4	1.33	5/7	345.53	753.45	348.29	57.22	1553.16

Table 14

Solutions results of MD-MFVRPDTW with different charging station distributions.

Distribution	$\mathit{r}\mathit{s}$	${\mathit{f}}_1/{\mathit{f}}_2$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{Z}$
C	2	5/6	345.66	74.80	1585.37
R	4	8/5	308.08	64.38	1671.95
CC	5	9/4	329.83	38.20	1696.90

Table 15

Solutions results of MD-MFVRPDTW under various carbon prices.

$\mathit{c}\mathit{p}$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{D}\mathit{R}\mathit{C}$	$\mathit{C}\mathit{T}\mathit{C}$	$\mathit{C}\mathit{E}$	$\mathit{Z}$
0	214.62	83.89	130.77	0	25.08	1419.62
0.1	220.06	79.51	128.34	−12.62	23.77	1440.54
0.2	246.01	78.35	125.75	−25.31	23.42	1469.20
0.3	296.46	77.60	117.59	−38.04	23.20	1493.75
0.4	330.61	76.63	115.96	−50.84	22.91	1558.88
0.5	365.94	72.38	109.83	−64.18	21.64	1582.04
0.6	397.18	70.83	108.66	−77.29	21.18	1606.52
0.7	418.66	68.79	107.47	−90.60	20.57	1619.91
0.8	435.96	65.06	105.36	−104.44	19.45	1642.70
0.9	460.78	63.55	101.50	−117.90	19.00	1679.47
1	488.43	60.11	97.47	−132.03	17.97	1718.62

Table 16

Solutions results of MD-MFVRPDTW under various carbon allowances.

$\overline{\mathit{C}\mathit{E}}$	$\mathit{E}\mathit{R}\mathit{C}$	$\mathit{F}\mathit{R}\mathit{C}$	$\mathit{D}\mathit{R}\mathit{C}$	$\mathit{C}\mathit{T}\mathit{C}$	$\mathit{C}\mathit{E}$	$\mathit{Z}$
0	351.58	72.25	108.29	10.80	21.60	1695.77
50	355.73	72.16	106.02	−14.21	21.57	1652.86
100	349.49	71.92	107.33	−39.25	21.50	1611.33
150	365.94	72.38	109.83	−64.18	21.64	1582.04
200	352.34	72.85	104.63	−89.11	21.78	1545.73
250	364.21	73.06	103.47	−114.08	21.84	1523.70
300	359.30	72.34	102.96	−139.19	21.63	1503.59

References

1. Lehmann, J.; Winkenbach, M. A metaheuristic for the two-echelon multi-trip vehicle routing problem with mixed pickup and delivery demand and time windows. Transp. Res. Part C Emerg. Technol; 2024; 160, 104522. [DOI: https://dx.doi.org/10.1016/j.trc.2024.104522]

2. Yaiprasert, C.; Hidayanto, A.N. AI-powered ensemble machine learning to optimize cost strategies in logistics business. Int. J. Inf. Manag. Data Insights; 2024; 4, 100209. [DOI: https://dx.doi.org/10.1016/j.jjimei.2023.100209]

3. Ritchie, H.; Roser, M. CO₂ Emissions. OurWorldinDat. Available online: https://ourworldindata.org/co2-emissions (accessed on 12 October 2024).

4. Guo, X.; Zhang, W.; Liu, B. Low-carbon routing for cold-chain logistics considering the time-dependent effects of traffic congestion. Transp. Res. Part D Transp. Environ.; 2022; 113, 103502. [DOI: https://dx.doi.org/10.1016/j.trd.2022.103502]

5. Yao, Q.; Zhu, S.; Li, Y. Green vehicle-routing problem of fresh agricultural products considering carbon emission. Int. J. Environ. Res. Public Health; 2022; 19, 8675. [DOI: https://dx.doi.org/10.3390/ijerph19148675] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35886525]

6. Liu, Z.; Niu, Y.; Guo, C.; Shitong, J. A vehicle routing optimization model for community group buying considering carbon emissions and total distribution costs. Energies; 2023; 16, 931. [DOI: https://dx.doi.org/10.3390/en16020931]

7. Zhang, C.; Tang, L.; Zhang, J.; Gou, L. Optimizing distribution routes for chain supermarket considering carbon emission cost. Mathematics; 2023; 11, 2734. [DOI: https://dx.doi.org/10.3390/math11122734]

8. Jabir, E.; Panicker, V.; Sridharan, R. Design and development of a hybrid ant colony-variable neighborhood search algorithm for a multi-depot green vehicle routing problem. Transp. Res. Part D Transp. Environ.; 2017; 57, pp. 422-457. [DOI: https://dx.doi.org/10.1016/j.trd.2017.09.003]

9. Wang, Y.; Assogba, K.; Fan, J.; Xu, M.; Liu, Y.; Wang, H. Multi-depot green vehicle routing problem with shared transportation resource: Integration of time-dependent speed and piecewise penalty cost. J. Clean. Prod.; 2019; 232, pp. 12-29. [DOI: https://dx.doi.org/10.1016/j.jclepro.2019.05.344]

10. Fan, H.; Zhang, Y.; Tian, P.; LV, Y.; Fan, H. Time-dependent multi-depot green vehicle routing problem with time windows considering temporal-spatial distance. Comput. Oper. Res.; 2021; 129, 105211. [DOI: https://dx.doi.org/10.1016/j.cor.2021.105211]

11. Xue, G. A two-stage heuristic solution for multi-depot collaborative pickup and delivery network with transfers to reduce carbon emissions. J. Clean. Prod.; 2022; 373, 133839. [DOI: https://dx.doi.org/10.1016/j.jclepro.2022.133839]

12. Kucukoglu, I.; Dewil, R.; Cattrysse, D. The electric vehicle routing problem and its variations: A literature review. Comput. Ind. Eng.; 2021; 161, 107650. [DOI: https://dx.doi.org/10.1016/j.cie.2021.107650]

13. Froger, A.; Mendoza, J.E.; Jabali, O.; Laporte, G. Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions. Comput. Oper. Res.; 2019; 104, pp. 256-294. [DOI: https://dx.doi.org/10.1016/j.cor.2018.12.013]

14. Zhang, J.; Yan, J.; Liu, Y.; Zhang, H.; Lv, G. Daily electric vehicle charging load profiles considering demographics of vehicle users. Appl. Energy; 2020; 274, 115063. [DOI: https://dx.doi.org/10.1016/j.apenergy.2020.115063]

15. Lee, C. An exact algorithm for the electric-vehicle routing problem with nonlinear charging time. J. Oper. Res. Soc.; 2021; 72, pp. 1461-1485. [DOI: https://dx.doi.org/10.1080/01605682.2020.1730250]

16. Wang, Y.; Zhou, J.; Sun, Y.; Sun, Y. Collaborative multi-depot electric vehicle routing problem with time windows and shared charging stations. Expert Syst. Appl. Int. J.; 2023; 219, 119654. [DOI: https://dx.doi.org/10.1016/j.eswa.2023.119654]

17. Ene, S.; Küçükoğlu, İ.; Aksoy, A.; Aksoy, A.; Öztürk, N. A hybrid metaheuristic algorithm for the green vehicle routing problem with a heterogeneous fleet. Int. J. Veh. Des.; 2016; 71, pp. 75-102. [DOI: https://dx.doi.org/10.1504/IJVD.2016.078771]

18. Islam, M.A.; Gajpal, Y. Optimization of conventional and green vehicles composition under carbon emission cap. Sustainability; 2021; 13, 6940. [DOI: https://dx.doi.org/10.3390/su13126940]

19. Al-Dalain, R.; Celebi, D. Planning a mixed fleet of electric and conventional vehicles for urban freight with routing and replacement considerations. Sustain. Cities Soc.; 2021; 73, 103105. [DOI: https://dx.doi.org/10.1016/j.scs.2021.103105]

20. Yu, V.F.; Jodiawan, P.; Gunawan, A. An adaptive large neighborhood search for the green mixed fleet vehicle routing problem with realistic energy consumption and partial recharges. Appl. Soft Comput.; 2021; 105, 107251. [DOI: https://dx.doi.org/10.1016/j.asoc.2021.107251]

21. Amiri, A.; Amin, S.H.; Zolfagharinia, H. A bi-objective green vehicle routing problem with a mixed fleet of conventional and electric trucks: Considering charging power and density of stations. Expert Syst. Appl.; 2023; 213, 119228. [DOI: https://dx.doi.org/10.1016/j.eswa.2022.119228]

22. Moshref-Javadi, M.; Hemmati, A.; Winkenbach, M. A comparative analysis of synchronized truck-and-drone delivery models. Comput. Ind. Eng; 2021; 162, 107648. [DOI: https://dx.doi.org/10.1016/j.cie.2021.107648]

23. Chiang, W.C.; Li, Y.; Shang, J.; Urban, T.L. Impact of drone delivery on sustainability and cost: Realizing the UAV potential through vehicle routing optimization. Appl. Energy; 2019; 242, pp. 1164-1175. [DOI: https://dx.doi.org/10.1016/j.apenergy.2019.03.117]

24. Han, Y.; Li, J.; Liu, Z.; Liu, C.; Tian, J. Metaheuristic algorithm for solving the multi-objective vehicle routing problem with time window and drones. Int. J. Adv. Robot. Syst.; 2020; 17, 1729881420920031. [DOI: https://dx.doi.org/10.1177/1729881420920031]

25. Zhang, S.; Liu, S.; Xu, W.; Wang, W. A novel multi-objective optimization model for the vehicle routing problem with drone delivery and dynamic flight endurance. Comput. Ind. Eng.; 2022; 173, 108679. [DOI: https://dx.doi.org/10.1016/j.cie.2022.108679]

26. Xiao, Q.; Gao, J. The Multi-Visit Vehicle Routing Problem with Drones under Carbon Trading Mechanism. Sustainability; 2024; 16, 6145. [DOI: https://dx.doi.org/10.3390/su16146145]

27. Baek, D.; Chen, Y.; Chang, N. Energy-efficient coordinated electric truck-drone hybrid delivery service planning. Proceedings of the 2020 AEIT International Conference of Electrical and Electronic Technologies for Automotive (AEIT AUTOMOTIVE); Turin, Italy, 18–20 November 2020; pp. 1-6.

28. Zhu, T.; Boyles, S.D.; Unnikrishnan, A. Electric vehicle traveling salesman problem with drone with partial recharge policy. arXiv; 2022; arXiv: 2205.13735

29. Kyriakakis, N.A.; Stamadianos, T.; Marinaki, M. The electric vehicle routing problem with drones: An energy minimization approach for aerial deliveries. Clean. Logist. Supply Chain; 2022; 4, 100041. [DOI: https://dx.doi.org/10.1016/j.clscn.2022.100041]

30. Mara, S.T.W.; Sarker, R.; Essam, D. Solving electric vehicle–drone routing problem using memetic algorithm. Swarm Evol. Comput.; 2023; 79, 101295. [DOI: https://dx.doi.org/10.1016/j.swevo.2023.101295]

31. Wang, X.; Wang, M.; Ruan, J.; Zhan, H. The multi-objective optimization for perishable food distribution route considering temporal-spatial distance. Procedia Comput. Sci.; 2016; 96, pp. 1211-1220. [DOI: https://dx.doi.org/10.1016/j.procs.2016.08.165]

32. Hou, D.; Fan, H.; Ren, X.; Tian, P.; Lv, Y. Time-dependent multi-depot heterogeneous vehicle routing problem considering temporal–spatial distance. Sustainability; 2021; 13, 4674. [DOI: https://dx.doi.org/10.3390/su13094674]

33. Wang, J.; Li, J. An improved genetic algorithm for vehicle routing problem with time windows considering temporal-spatial distance. Proceedings of the International Conference on Intelligent Computing; Zhengzhou, China, 10–13 August 2023; Springer Nature Singapore: Singapore, 2023; pp. 397-409.

34. Montoya, A.; Guéret, C.; Mendoza, J.E.; Villegas, J.G. The electric vehicle routing problem with nonlinear charging function. Transp. Res. Part B Methodol.; 2017; 103, pp. 87-110. [DOI: https://dx.doi.org/10.1016/j.trb.2017.02.004]

35. Ichoua, S.; Gendreau, M.; Potvin, J.Y. Vehicle dispatching with time-dependent travel times. Eur. J. Oper. Res.; 2003; 2, pp. 379-396. [DOI: https://dx.doi.org/10.1016/S0377-2217(02)00147-9]

36. Demir, E.; Bektaş, T.; Laporte, G. An adaptive large neighborhood search heuristic for the pollution-routing problem. Eur. J. Oper. Res.; 2012; 223, pp. 346-359. [DOI: https://dx.doi.org/10.1016/j.ejor.2012.06.044]

37. Keskin, M.; Çatay, B. Partial recharge strategies for the electric vehicle routing problem with time windows. Transp. Res. Part C Emerg. Technol.; 2016; 65, pp. 111-127. [DOI: https://dx.doi.org/10.1016/j.trc.2016.01.013]

38. Schneider, M.; Stenger, A.; Goeke, D. The electric vehicle-routing problem with time windows and recharging stations. Transp. Sci.; 2014; 48, pp. 500-520. [DOI: https://dx.doi.org/10.1287/trsc.2013.0490]