A Two-Stage Clustering–Evolutionary Framework with Adaptive Search Control for the Vehicle Routing Problem with Time Windows

As global e-commerce develops rapidly, logistics distribution, the critical link between supply and demand, has become a key factor of enterprise competitiveness. According to recent statistics, the total social logistics volume in China had exceeded 331.2 trillion Chinese Yuan (CNY) by November 2025, with a significant increase in express delivery demand. However, several challenges still hinder the quality development of the logistics industry, such as “last-mile” delivery, mismatched allocation of transportation resources, and the high terminal distribution cost.

The vehicle routing problem, particularly the multi-depot vehicle routing problem, is a classical combinatorial optimization problem in logistics [1]. Its fundamental objective is to identify an optimal set of delivery routes, including vehicle capacity, travel distance, and time window restrictions, minimizing total operational cost and satisfying multiple constraints. The vehicle routing problem with time windows has attracted a lot of attention because it is closely related to real-world logistics. Since 2020, research on the vehicle routing problem with time windows has changed toward multi-strategy integration in complex operation. In terms of algorithmic model innovation, Tamke and Buscher [2] proposed branch-and-cut-based routing optimization in drone-assisted delivery systems. Similarly, Weng et al. [3] investigated multi-depot vehicle routing incorporating drone collaboration for emergency logistics, extending the application scope in complex scenarios. Beyond traditional ground logistics, the integration of autonomous and unmanned technologies has introduced new dimensions to routing optimization. In the context of urban delivery, Morim et al. [4] investigated drone-assisted vehicle routing incorporated with robot stations, and Lupi et al. [5] proposed a novel optimization system based on platoons of automated vehicles for last-mile urban freight. To solve the dynamic vehicle routing problem with time windows in these environments, Li et al. [6] proposed an elastic strategy-based adaptive genetic algorithm. This method effectively enhances the global search capability of the algorithm and significantly improves its solution accuracy and adaptability in dynamic and complex logistics environments. From the perspective of heuristic algorithm enhancement, Niazy et al. [7] proposed hybrid optimization combining the chicken swarm optimization with genetic algorithms. This greatly improved the solution accuracy for large-scale capacitated vehicle routing. In terms of heuristic optimization, Rodrigues [8] designed a novel genetic algorithm to effectively solve capacitated vehicle routing problems (VRPs). Building upon such advanced optimization strategies, Thammasang and Arunyanart [9] addressed the unique constraints of ice distribution, while Moghaddasi et al. [10] and Jiang et al. [11] adapted these heuristic methods for cold chain logistics and fresh fruit harvest systems, respectively. Furthermore, regarding dynamic and multi-objective routing models in specialized fields, Shen et al. [12] developed a two-stage scheduling algorithm for medical waste collection. In the aviation sector, Bao et al. [13] proposed a proactive real-time scheduling method for apron service vehicles, and subsequently Bao et al. [14] optimized the scheduling of electric aircraft tow tractors under operator cooperation.

The concept of spatial dimensionality reduction has been widely adopted in preprocessing and collaborative optimization. Wang et al. [15] proposed a clustering-based extended genetic algorithm, addressing the multi-depot vehicle routing problem with three-dimensional loading constraints. In retail logistics distribution, the modeling and optimization of the multi-depot logistics network conducted by Mogale et al. [16] enhanced distribution efficiency, providing valuable insights into spatial decomposition strategies for large-scale logistics systems. In addition, Yan et al. [17] investigated a split-delivery vehicle routing problem with time window and three-dimensional loading constraints, with tailored optimization addressing the associated complexities.

To address the growing complexity of vehicle routing problems (VRP) in dynamic and diverse logistics environments. Mardešić et al. [18], researchers have heavily focused on advanced clustering techniques and heuristic algorithms. Silva et al. [19] proposed a cluster branching technique to enhance decision-making efficiency. Lau et al. [20] extensively investigated the vehicle routing problem with time windows and a limited number of vehicles. By optimizing the allocation of constrained vehicle resources, their study provides an effective routing solution for complex logistics scenarios with limited fleet sizes.

These contributions have further enriched the research on the vehicle routing problem with multiple constraints.

However, limitations remain in the application of conventional genetic algorithms to the vehicle routing problem with time windows. First, fixed crossover and mutation probabilities often cause premature convergence during the later evolution stages, leading to entrapment in local optima. Second, research on optimized facility location and routing with multi-depot and large-scale demand remains fragmented, and systematic dynamic validation is often lacking. Although adaptive strategies have been utilized, the effective integration of clustering analysis with dynamic simulation for end-to-end optimization (from facility location to route planning) has not been fully explored.

In summary, hybrid optimization integrating k-means clustering and an adaptive genetic algorithm is proposed for logistics distribution routing, with empirical validation using the AnyLogic platform. In the first-stage solution, k-means clustering is utilized for scientific partitioning of customer nodes and distribution center allocation, thereby reducing the computational complexity of route optimization. In the second-stage solution, an improved adaptive genetic algorithm is designed for the vehicle routing problem model with time windows, with adaptive operators dynamically adjusting crossover and mutation probabilities based on individual fitness values. This enhances the global search capability. In the real-world case study based on operational data from Company S, a dynamic simulation model is developed using AnyLogic. Comparative analyses show the changes of key indicators, including total operational cost, vehicle load factor, and customer satisfaction, before and after optimization. The results confirm that the proposed approach is feasible and superior in complex logistics.

Logistics distribution routing optimization is a two-stage optimization framework composed of distribution center allocation and vehicle routing optimization. In the first stage, improved k-means clustering addresses the spatial hierarchical structure between sorting centers and large-scale customer nodes. The aim is to minimize the within-cluster sum of squared errors through iterative computation using customer coordinates and demand distributions. Then the optimal locations of secondary distribution centers are determined within a complex geographic network. This stage not only optimizes the hierarchical configuration of the logistics network but also achieves spatial decomposition. This can reduce the computational dimensionality of subsequent routing optimization. In the second stage, the vehicle routing problem with time windows is optimized based on the distribution center locations. Delivery vehicles are dispatched from their respective distribution centers to satisfy operational constraints, including vehicle capacity limits, maximum travel distance per trip, and soft time window requirements concerning customer service deadlines. The optimal routing should balance service efficiency and operational cost. The overall objective of the system is to minimize total operational cost by optimizing location and routing decisions. The cost function incorporates fixed vehicle activation costs, distance-dependent transportation costs, and penalty costs concerning time window violations.

The following assumptions were proposed for easy computation and consistency with practical logistics:

Assumption 1: The locations of distribution centers are fixed, and all delivery vehicles should return to their respective origin distribution centers after completing service tasks.

Assumption 2: The demand at each customer node is known in advance and does not exceed the rated capacity of a single vehicle.

Assumption 3: Each customer node is served exactly once by a single vehicle.

Assumption 4: Vehicles travel at a constant speed, without considering stochastic factors such as traffic congestion and other unforeseen disruptions.

Assumption 5: All delivery vehicles have identical specifications, capacity limits, and cost parameters.

Table 1 shows the model variables and corresponding definitions.

The optimization objective is to minimize the total distribution cost ($Z$), including fixed vehicle activation cost (Z$_1$), transportation cost concerning travel distance (Z$_2$), and time window penalty cost (Z$_3$) using a linear penalty function.

$ Z=\min \left(Z_1+Z_2+Z_3\right) $

The three kinds of costs are expressed as follows:

$ \begin{aligned} &Z_1=\sum_{k \in V} \sum_{j \in N} x_{0 j k} \cdot c_2\\ &Z_2=\sum_{\mathrm{k} \in V} \sum_{i \in N \cup\{0\}} \sum_{j \in N \cup\{0\}} x_{i j k} \cdot d_{i j} \cdot c_1\\ &Z_3=\sum_{k \in V} \sum_{i \in N}\left[w_1 \cdot \max \left(E T_i-T_{i k}, 0\right)+w_2 \cdot \max \left(T_{i k}-L T_i, 0\right)\right] \end{aligned} $

In terms of the capacity constraint, the total load assigned to each vehicle must not exceed its rated capacity. In terms of the service uniqueness constraint, each customer node must be visited exactly once and served by a single vehicle. The constraints are expressed as follows:

$ \begin{aligned} & \sum_{i \in N} q_i y_{i k} \leq Q, \forall k \in V \\ & \sum_{k \in V} y_{i k}=1, \forall i \in N \end{aligned} $

As for the flow conservation constraint, if a node is entered for each vehicle, it must subsequently be exited. As for the distribution center constraint, each vehicle must depart from its assigned distribution center and return to the same distribution center. The constraints are expressed as follows:

$ \begin{aligned} &\sum_{i \in N \cup\{0\}} x_{i j k}=\sum_{i \in N \cup\{0\}} x_{i j k}=y_{j k}, \forall j \in N, \forall k \in V\\ &\sum_{j \in N} x_{0 j k}=\sum_{i \in N} x_{i 0 k} \leq 1, \forall k \in V \end{aligned} $

The time continuity constraint and variable constraint are expressed as follows:

$ \begin{aligned} &\mathrm{x}_{ijk}=1 \Rightarrow T_{i k}+s_i+t_{i j}=T_{j k}, \forall i, j \in N, \forall k \in V\\ &x_{i j k} \in\{1,0\}, y_{j k} \in\{0,1\} \end{aligned} $

A two-stage solution is developed to optimize large-scale logistics distribution under time window constraints. k-means clustering is employed in the first stage for scientific distribution region partitioning. An adaptive genetic algorithm is used in the second stage to optimize the routes within each subregion.

To reduce the computational complexity concerning the large-scale vehicle routing problem, customer nodes are initially clustered based on their geographical coordinates. The number of clusters ($K$) is determined using the silhouette coefficient as an evaluation metric. The silhouette coefficient $S$($i$) integrating intra-cluster cohesion and inter-cluster separation ranges from -1 to 1. A value closer to 1 indicates more appropriate clustering. The optimal value of K is identified by traversing candidate values and selecting the one that maximizes the average silhouette coefficient. During the clustering process, a set of K initial cluster centroids is randomly selected. The centroids are iteratively updated until the total Euclidean distance between each customer node and its assigned cluster centroid is minimized:

$ J=\sum_{j=1}^K \sum_{i \in C_j}\left\|x_i-\mu_j\right\|^2 $

where, $x_i$ denotes the customer coordinates, and $\mu_j$ represents the centroid of cluster $j$. Figure 1 shows the flowchart of the k-means clustering algorithm.

The procedural framework of the adaptive genetic algorithm is described in Figure 2.

Step 1: Encoding

To enable processing within the genetic algorithm framework, candidate solutions are encoded into chromosome representations. This transformation between problem solutions and chromosome structures is referred to as encoding and decoding.

Step 2: Population Initialization

The initial population generation starts from the random sampling of a predefined number of individuals from the solution space. A chromosome represents an individual. The chromosomes have various forms, such as binary strings, real-valued vectors, or other suitable structures.

Step 3: Fitness Evaluation

The fitness function constructed aims to evaluate the population quality at the individual level. The fitness value determines the probability of selecting an individual as a parent for reproduction and genetic propagation. Individuals with higher fitness values exhibit superior performance in solving the optimization problem.

Step 4: Selection

Appropriate selection strategies are used to select individuals based on their fitness values. These selections are based on rank or involve a roulette wheel or tournament. Certain individuals are chosen for subsequent crossover and mutation.

Step 5: Crossover and Mutation

After two parent individuals are selected, new offspring are generated via crossover. Then the application of mutation to the offspring introduces new genetic characteristics, enhancing population diversity.

Step 6: Parameter Adaptation

The crossover and mutation probabilities are dynamically adjusted according to the fitness distribution and diversity of the current population.

Step 7: Iteration

Steps 2 through 5 are iteratively repeated for a predefined number of generations. When reaching the maximum number of iterations, achieving a predefined fitness threshold, or obtaining stable fitness values over successive generations, the iterative process is terminated.

Step 8: Output

In the final generation, the individual with the highest fitness value is considered the optimal solution. The results are analyzed in the aspects of solution quality, convergence speed, and computational complexity.

In conventional genetic algorithms, the crossover probability ($P_c$ ) and mutation probability ($P_m$) are typically fixed, which may lead to premature convergence or reduced search efficiency. To address these limitations, an adaptive mechanism is introduced for dynamic genetic operator adjustments. The integer encoding scheme adopted is based on permutation. Specifically, the integers 1, 2, …, $n$ are used to denote customer nodes, with 0 representing the distribution center. With a chromosome as a sequence, multiple vehicle routes in the sequence are encoded by inserting the delimiter 0 to separate different routes. A typical chromosome is expressed as $0 \rightarrow i_1 \rightarrow i_2 \rightarrow \ldots \rightarrow 0 \rightarrow i_j \rightarrow \ldots \rightarrow 0$.

The optimization objective is to minimize the total distribution cost Z. Therefore, the fitness function f is the reciprocal of the total cost:

$ f=\frac{1}{Z} $

The total cost Z consists of transportation cost associated with travel distance, fixed vehicle activation cost, and time window penalty cost.

The vehicle routing problem with time windows is non-deterministic polynomial-time hard, involving multiple coupled constraints (e.g., time window, vehicle capacity, and travel distance constraints). The adaptive genetic algorithm, through its population-based search mechanism, handles multiple constraints and identifies feasible solutions that approximate the global optimum. Furthermore, premature convergence is mitigated by maintaining population diversity. This is very critical for complex optimization problems because diversity enhances the solution space exploration. In addition, the adaptive genetic algorithm integrates with other optimization techniques, such as local search methods and simulated annealing, due to its flexibility and extensibility. This facilitates the hybrid optimization development, further improving solution quality and computational efficiency. Based on these considerations, the algorithmic design is formulated below.

(a) Encoding Scheme

The encoding scheme is critical in the genetic algorithm, as it directly influences algorithmic performance and solution efficiency. The choice of encoding determines the representation of individuals and the application of genetic operators (e.g., crossover and mutation) to individuals. In the vehicle routing problem, commonly adopted encoding schemes include binary encoding, floating-point encoding, and integer-based encoding. Although binary encoding is straightforward to implement in routing problems, it often leads to a lot of infeasible solutions. Floating-point encoding offers higher numerical precision, but it often has lower efficiency in local search operations. Consequently, permutation-based integer encoding is more suitable for the vehicle routing problem.

The permutation-based encoding scheme enables direct mapping between chromosomes and vehicle routes without an explicit decoding process. As a result, the simplified representation and computational procedures can improve the computational efficiency. Under this scheme, customer nodes are encoded as a sequence (1, 2, 3, …, $N$), and vehicle indices are encoded as a sequence (1, 2, 3, …, $k$) for subsequent crossover and mutation. As for the encoding scheme of a distribution system (one depot, eight customer nodes, and two delivery vehicles), integers 1 and 2 represent the two vehicles, integers 3 to 10 correspond to the eight customers, and the chromosome length is set to 10. The structure is illustrated in Figure 3. This encoding approach completely represents the distribution routes, enhances the efficiency of the genetic algorithm, and reduces computational complexity.

With this encoding scheme, the genetic algorithm operates more efficiently in solving the vehicle routing problem. This reduces computational complexity and improves solution quality. The chromosome corresponds to a routing solution in which vehicle 1 serves customer nodes in the sequence $4 \rightarrow 7 \rightarrow 10 \rightarrow 3 \rightarrow 9$, while vehicle 2 serves customer nodes in the sequence $6 \rightarrow 5 \rightarrow 8$.

(b) Population Initialization

As for the random initialization strategy adopted, an initial population consisting of 100 chromosomes is constructed. The sequence of elements within each chromosome is randomly permuted to produce new individuals or solutions. The fitness value of each chromosome is then evaluated, and the population is structured accordingly for subsequent operations. This ensures that chromosomes corresponding to different distribution nodes are processed rationally. The primary advantage of this approach lies in its simplicity and effectiveness. Random initialization within an appropriate range increases the likelihood of identifying globally optimal or near-optimal solutions. Moreover, it reduces the risk of the algorithm becoming trapped in local optima. This initialization method is broadly applicable and widely covers the solution space, thereby enhancing population diversity.

(c) Fitness Function

In the genetic algorithm, the fitness function is a key criterion for evaluating the individual quality. A higher fitness value indicates superior performance of the corresponding individual. The fitness value is computed as the reciprocal of the objective function. Based on the formulated objective function, the fitness function is expressed as follows:

$${ fitness }(i)=\frac{1}{Z}\left(i=1,2, \ldots, N_i\right)$$

where, $ \mathit{fitness}(i) $ and $Z$ denote the fitness value and the objective function of chromosome $i$; and $N_i$ denotes the population size.

(d) Selection Operator

Tournament selection is adopted in the genetic algorithm. From the generated chromosome population, two individuals are randomly sampled to form a tournament group, and their fitness values are compared. The individual with the highest fitness is selected as the “winner” and carried forward to the next generation. This method is computationally efficient and straightforward to implement. Therefore, it is very suitable for large-scale populations. The stochastic selection preserves individuals with diverse characteristics, thereby reducing the risk of premature convergence. The probability $P$ of an individual being selected as the winner can be expressed as follows:

$ P_{\text {(iwnins })}=\frac{{ fitness }(i)}{\sum_{\mathrm{j}=1}^N { fitness }(i)} $

where, fitness (i) denotes the fitness value of each individual, and N denotes the tournament size, i.e., the number of individuals randomly selected in each tournament. A schematic representation of the minimum tournament selection tree is illustrated in Figure 4.

(e) Crossover Operator

Various crossover strategies have been proposed in genetic algorithms. Uniform crossover, although simple, does not preserve the relative order of genes. Single-point and two-point crossover methods can partially maintain gene sequence continuity; however, they often disrupt gene integrity. The core principle of the order crossover method is to preserve both the relative order and uniqueness of genes, thereby ensuring the feasibility of offspring solutions.

In the order crossover process, two parent chromosomes are first selected, and two crossover points are randomly determined. A subsequence of genes from the first parent is copied directly into the corresponding positions of the offspring chromosome. Subsequently, the remaining genes are filled by scanning the second parent in sequence, excluding genes that have already been inherited, and inserting the remaining genes into the unassigned positions in order. For example, a gene segment {5, 3, 4} may be selected from the first parent and placed into the offspring at the corresponding positions. The remaining genes are then taken sequentially from the second parent. Each gene is checked for duplication; if it is not already present in the offspring, it is inserted into the next available position. Upon reaching the end of the second parent, the process continues from the beginning of the chromosome. Genes already present in the offspring (e.g., gene 5) are skipped, and the procedure proceeds to the next gene (e.g., gene 1). This process is repeated until the results in Figure 5 can be obtained. This strategy effectively avoids the occurrence of duplicate or missing genes, which are common issues in conventional single-point or two-point crossover methods, thereby ensuring the validity of the generated routing solutions. The crossover procedure is illustrated in Figure 5.

(f) Mutation Operator

During the swap mutation operation process, each gene within a chromosome is sequentially examined, and mutation is performed based on a predefined mutation probability. When a gene is selected for mutation, another gene is randomly chosen from the same chromosome, and their positions are exchanged. The swap mutation process is illustrated in Figure 6.

(g) Adaptive Crossover and Mutation Probabilities

In the adaptive genetic algorithm, prior to the execution of crossover and mutation operations, the average fitness of the population is computed. This average fitness value is then compared with the fitness of each individual, based on which the crossover probability and mutation probability are assigned to each individual.

The parameters $P_c$ and $P_m$ are dynamically adjusted according to the evolutionary state of the population. When the fitness values of individuals tend to converge or the population approaches a local optimum, the probabilities are increased to enhance population diversity. Conversely, when the population exhibits high diversity, the probabilities are reduced to preserve high-quality genetic information. It can be expressed as follows:

$P_{\mathrm{c}}=\left\{\begin{array}{l}k_1 \frac{f_{\max }-f^{\prime}}{f_{\max }-f_{a v g}}, f^{\prime} \geq f_{a v g} \\ k_2, f^{\prime} \lt f_{a v g}\end{array}\right.$

$P_{\mathrm{m}}=\left\{\begin{array}{l}k_3 \frac{f_{\max }-f}{f_{\max }-f_{a v g}}, f \geq f_{a v g} \\ k_4, f \lt f_{a v g}\end{array}\right.$

where, $f_{\max}$ denotes the maximum fitness value within the population, $f_{\text {avg}}$ represents the average fitness value of the population in each generation, $f^{\prime}$ denotes the larger fitness value between the two individuals selected for crossover, $f$ represents the fitness value of the individual subject to mutation, and $k_1, k_2, k_3$, and $k_4$ are constants within the interval (0, 1). Table 2 shows the main parameters of the genetic algorithm.

A real-world last-mile logistics distribution scenario from Company S in Shanghai was selected as the case study. The distribution system consisted of two candidate distribution centers (DC$_1$ and DC$_2$) and 60 customer nodes randomly distributed across the service area. The geographical coordinates (latitude and longitude), demand quantities ($q_i$), and service time window requirements [$E T_i, L T_j$] were obtained from the company’s operational database. Table 3 and Table 4 show the parameters of the vehicle and the improved genetic algorithm, respectively. Figure 7 shows the geographic distribution and demand of customer nodes.

Based on the solution obtained using the adaptive genetic algorithm, four and eleven delivery vehicles were allocated to the distribution centers DC$_1$ and DC$_2$, respectively. The optimized routing configuration is highly compact, and route intersections are effectively avoided.

The best fitness value within the population is recorded at each generation. By comparing the optimal fitness value of the current generation with that of the previous generation, the algorithm optimization can be evaluated. An increase in the best fitness value indicates continued improvement in solution quality, whereas stagnation suggests the possibility of adjusting parameters (e.g., crossover and mutation probabilities) or terminating the criteria. This feedback mechanism monitors algorithmic performance in real time and ensures effectiveness and convergence stability. The evolving optimal solutions for the two distribution centers across iterations are illustrated in Figure 8.

As illustrated in Figure 8, the fitness values of both distribution centers show a decreasing trend over successive iterations. This indicates that the adaptive genetic algorithm effectively improves solution quality.

Based on the constructed function framework, the output results of the model were determined. A dedicated procedure was implemented to extract the optimal chromosome, obtaining key performance indicators, including fitness value, time window penalty cost, transportation cost, fixed vehicle activation cost, and chromosome representation. Five simulation runs were performed because the algorithm is stochastic, and the best-performing solution was selected. The corresponding results are presented in Table 5.

Table 5 shows a significant reduction in total distribution cost. The distribution costs in the two service regions are reduced by 28.12% and 20.62%, respectively, because the adaptive genetic algorithm has a strong global optimization capability in handling time window constraints. In addition, the vehicle utilization efficiency is greatly improved. The load factor of most delivery vehicles exceeds 90% after optimization, with seven vehicles in distribution center DC$_2$ reaching full capacity utilization (100%). This indicates a marked enhancement in utilizing transportation resources.

The AnyLogic simulation model was executed. The vehicle distribution process is presented in Figure 9.

Table 6 and Table 7 show the vehicle routing schemes and operational details for the two distribution centers.

The tables show that four delivery vehicles are allocated to distribution center DC$_1$, with load factors exceeding 90% for all vehicles. Among them, Vehicle 1 achieves full capacity utilization (100%). For distribution center DC$_2$, seven vehicles achieve full capacity utilization (100%) among all ten delivery vehicles, while only two vehicles (83% and 88%) exhibit load factors below 90%. These results indicate that vehicle load factors are significantly improved after optimization. Furthermore, an analysis of distribution costs shows substantial cost reductions in both centers (28.12% and 20.62%). A detailed comparative analysis is presented in Table 8.

The conventional genetic algorithm exhibits notable limitations in addressing large-scale, multi-customer distribution problems. Specifically, its solutions tend to produce longer routing distances and lower vehicle load factors without effectively considering customer time window constraints. This leads to higher operational costs and reduced customer satisfaction. In contrast, the adaptive genetic algorithm obtains significant improvements. The introduction of adaptive crossover and mutation functions substantially enhances the convergence performance and the quality of the optimal solutions. These findings indicate that the proposed approach effectively overcomes the limitations of the conventional genetic algorithm. The optimized routing solutions exhibit shorter travel distances, higher vehicle load factors, and lower time window penalty costs. Consequently, this improves customer service levels, significantly reduces overall distribution costs, and reduces overall logistics system costs.

An integrated optimization framework combining k-means clustering and an adaptive genetic algorithm was proposed to address the vehicle routing problem with time windows. The effectiveness of the proposed approach was further validated through simulation experiments conducted on the AnyLogic platform. Based on the case study of the logistics distribution system of Company S, the following conclusions were drawn:

Algorithmic performance: The proposed adaptive genetic algorithm, incorporating dynamic adjustment mechanisms for crossover and mutation probabilities, effectively mitigates premature convergence and entrapment in local optima commonly observed in the conventional genetic algorithm. Compared with the conventional genetic algorithm, the superior global search capability and faster convergence of the adaptive genetic algorithm can identify high-quality routing solutions more efficiently.

•Solution effectiveness: The two-stage optimization strategy—comprising k-means-based regional partitioning and adaptive genetic algorithm-based routing optimization—successfully reduces the dimensionality for large-scale distribution problems. Simulation results from the AnyLogic platform confirm that the optimized routing solutions are feasible under dynamic delivery conditions. Furthermore, significant improvements in vehicle load factor are observed while maintaining compliance with customer time window constraints.

•Economic benefits: The empirical analysis based on Company S demonstrates the practical applicability and economic value of the proposed model and algorithm. The operational costs of distribution centers DC$_1$ and DC$_2$ are reduced by 28.12% and 20.62%, respectively. These results indicate that the proposed approach not only effectively reduces logistics costs but also provides a scientifically grounded decision-support framework for addressing the “last-mile” delivery challenge in urban logistics systems.

Conceptualization, L.L.; methodology, L.L.; software, J.R.W.; validation, Y.J.L.; formal analysis, J.R.W; investigation, J.R.W.; resources, L.L.; data curation, L.L.; writing—original draft preparation, Y.J.L.; writing—review and editing, J.R.W; visualization, L.L.; supervision, L.L.; project administration, L.L.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript