Adaptive Routing Policies for Stochastic EVRPTW via Genetic Programming

1. Introduction

The vehicle routing problem (VRP) is a fundamental combinatorial optimization problem that has been extensively studied in the operations research literature [1]. Its practical relevance arises from a wide range of applications in transportation, logistics, and production systems [2]. Over time, numerous variants have been developed to better reflect operational constraints, including the vehicle routing problem with time windows (VRPTW) and pickup-and-delivery settings [3,4].

In the context of sustainable transportation, electric vehicle routing problems (EVRPs) have received increasing attention [5,6,7]. In these problems, routing decisions must consider not only customer service requirements, but also battery capacity, charging time, and charging infrastructure constraints. In particular, logistics systems using electric vehicles must explicitly address limited driving ranges and the need for recharging along routes, which is especially critical in mid- and long-haul operations [8].

An important development in EVRP research is the electric vehicle routing problem with time windows (EVRPTW) introduced i n [9], which extends the classical Solomon benchmark framework [10] and has become the standard reference model. Despite substantial progress in EVRPTW, most existing approaches rely on deterministic assumptions, where customer demand, travel times, and service durations are known in advance [11,12,13].

In practice, routing systems are inherently stochastic, with uncertainty in demand, travel time, and service duration, all of which can significantly affect routing decisions. Such sources of uncertainty have been recognized in the stochastic VRP literature [14,15], yet their integration into EV routing remains limited, with existing approaches typically relying on stochastic programming, simulation, or reinforcement learning frameworks [16,17,18]. Although expressive, these methods are computationally demanding and less suitable for fast or real-time decision-making. Recent work continues to explore how uncertainty can be modeled and incorporated into EV routing, highlighting that this remains an open and actively researched problem [19].

A complementary direction is provided by recent GP-based hyper-heuristic approaches for EVRPTW [20], which generate routing policies via a route generation scheme (RGS) instead of relying on iterative solution improvement, enabling fast and adaptive decision-making. However, such approaches have so far been developed only under deterministic assumptions, leaving a gap between stochastic modeling and policy-based routing.

Genetic programming, as a form of evolutionary computation, is directly inspired by biological evolution, where populations of candidate solutions evolve through selection, variation, and inheritance [21]. These mechanisms enable the emergence of adaptive behaviors without explicit design, reflecting how natural systems respond to changing environments. In the context of routing under uncertainty, this perspective is particularly relevant, as effective solutions must continuously adapt to evolving conditions rather than rely on fixed plans.

This paper addresses this gap by extending the GP-based RGS framework to stochastic EVRPTW, incorporating uncertainty in customer demand, travel time, and service duration while retaining the standard Schneider-based problem structure. The proposed approach introduces stochastic variants of the route generation scheme and extends the GP representation to support decision-making under uncertainty.

The main contributions of this work are as follows:

• Two extensions of the route generation scheme (RGS), namely Semi-parallel-B and Parallel-B, are introduced. These variants incorporate candidate-based vehicle selection to improve robustness under stochastic demand conditions.
• The GP-based priority function representation is extended with stochastic and global-state descriptors, as well as objective-specific terminals such as $CminV$ and $SlackSelf\left(i\right)$.
• A stochastic adaptation of the Schneider EVRPTW benchmark instances is developed, incorporating controlled variability in demand, service time, and travel time while preserving the original problem structure.
• A systematic experimental evaluation is conducted across multiple stochastic scenarios and objectives, analyzing the impact of uncertainty on routing performance and comparing different RGS variants.

The remainder of this paper is organized as follows. Section 2 reviews related work on electric vehicle routing and GP-based hyper-heuristics. Section 3 presents the proposed methodology, including the stochastic problem setting, the routing policy framework, and the GP-based approach for evolving priority functions. Section 4 describes the experimental setup and reports the computational results across different objectives and stochastic scenarios. Section 5 discusses the observed behavior of the routing strategies and their interaction with uncertainty. Finally, Section 6 concludes the paper and outlines directions for future research.

2. Related Work

The electric vehicle routing problem with time windows was introduced by [9] as an extension of the classical VRPTW benchmark from [10], incorporating electric vehicle constraints such as limited battery capacity and charging requirements, and has since become the standard reference in the field.b

Most existing work in this domain assumes deterministic problem parameters and relies on metaheuristic optimization methods, particularly adaptive large neighborhood search (ALNS) [11,22,23], iterative local search (ILS) [12], and variable neighborhood search (VNS) [13]. In addition, exact and hybrid methods have also been proposed, including column generation and branch-and-price-based approaches [24,25,26,27,28,29].

In contrast to deterministic formulations, stochastic vehicle routing problems explicitly model uncertainty in key parameters. Early work [30] established fundamental properties of routing with stochastic demands, including the distinction between full and split delivery and different modes of information revelation. A comprehensive overview of stochastic routing variants, including stochastic customers, stochastic demands, and stochastic travel times, is provided in [14]. Similarly, stochastic VRP is classified into three main classes: stochastic customers, stochastic demands, and stochastic times [15]. The authors also discuss solution approaches such as chance-constrained programming and stochastic programming with recourse. These frameworks typically follow a two-stage structure, where an initial solution is constructed before uncertainty is realized, followed by corrective actions once stochastic information becomes available.

Several works have addressed specific stochastic VRP variants. In [31] routing problems with stochastic travel and service times were studied, while in [32] the authors investigated stochastic travel times using branch-and-cut approaches combined with sampling. For stochastic demands, exact and decomposition-based methods were proposed in [33,34], and reinforcement learning approaches were explored in [35]. Although these studies provide important theoretical insights, they are generally based on exact optimization, decomposition, or scenario-based modeling, which limits their applicability in large-scale or real-time settings.

Compared to the classical VRP, stochastic extensions of EVRP and EVRPTW are relatively limited and typically focus on specific sources of uncertainty. For example, stochastic waiting times at charging stations were introduced [16], which modelled waiting time as a random variable revealed upon arrival. Other works consider stochastic customer demand and dynamic routing decisions. In [17] the authors studied EVRP with stochastic demands and dynamic remedial measures using a mixed-integer programming framework, while in [18] the authors addressed a dynamic stochastic EVRP using reinforcement learning. More recent contributions extend stochastic modeling to additional EV-specific factors: in [19] a stochastic electric vehicle location–routing problem incorporating demand uncertainty and environmental effects such as ambient temperature was considered, while in [36] stochastic energy consumption using a two-stage stochastic programming approach with recourse policies was modeled. Despite capturing important real-world aspects, these approaches rely on scenario-based optimization and computationally intensive methods, which limits their applicability in fast or real-time decision-making settings.

Overall, both deterministic metaheuristics and stochastic optimization approaches typically rely on computationally demanding procedures and fixed or scenario-based solution structures, making them less suitable for dynamic environments that require fast and adaptive decision-making. This limitation motivates the exploration of alternative solution paradigms.

Constructive heuristics represent an alternative class of methods that build solutions incrementally rather than improving an existing one. They are commonly used to generate initial solutions or to repair partial solutions within hybrid frameworks. Due to their relatively low computational cost, constructive methods are particularly suitable for dynamic and real-time environments. However, designing effective constructive heuristics is non-trivial, and such approaches are often less explored compared to improvement-based metaheuristics [37].

Genetic programming (GP) [21] has been increasingly used for automatic design of heuristics, particularly in scheduling [38,39] and graph-based optimization problems [40,41]. In the context of vehicle routing, however, its application remains limited. In [42] GP was for the first time used to generate heuristics for VRP. More specifically, the authors considered a dynamic variant of the problem in which customers arrive over time and can be accepted or rejected. The experimental results demonstrate that GP evolved heuristics achieve significantly better results compared to manually designed ones. A similar problem was also considered in [43], in which problems with different levels of dynamic customer arrivals were considered. The results again demonstrated that GP evolved heuristics performed better compared to the manually designed nearest neighbor heuristic. In [41] the authors considered a completely different problem variant, namely VRP with zone-based pricing. This variant models the situation in which profit needs to be maximized by determining which customers should be served. The results demonstrated that GP evolved heuristics achieved better results compared to a genetic algorithm, and were even competitive to a branch & price algorithm, thus outlining the suitability of GP for evolving heuristics for different problem variants.

One of the first notable applications of a GP-based hyper-heuristic framework for EVRPTW was done in [20]. This approach combines GP with a route generation scheme (RGS) to evolve routing policies instead of directly optimizing complete solutions. This shifts the focus from solution search to heuristic generation, enabling faster decision-making and improved adaptability in dynamic settings. This methodology has been further extended by using surrogate models to improve the performance of GP, as well as its convergence [44].

Overall, GP-based approaches provide an efficient and adaptive alternative to traditional EVRPTW solution methods, as they learn routing policies that can be applied quickly in dynamic settings. However, all the previous applications of GP on VRP or EVRP have been done on deterministic problem settings, completely disregarding the uncertainty that can be present in such problems. This highlights the need for extending this methodology further to also cover stochastic variants of the EVRP problem. This adaptive capability reflects the underlying evolutionary principles from which GP originates, making it particularly suitable for decision-making in dynamic and uncertain environments.

3. Materials and Methods

3.1. Electric Vehicle Routing Problem with Time Windows

EVRPTW extends the classical VRP by considering electric vehicles and time-constrained service at customers. The formulation introduced by [9] has become the standard reference model, and most existing EVRPTW studies are based on this framework. Each customer must be serviced within a predefined time window, while routing decisions are constrained by limited battery capacity, energy consumption, and recharging requirements. In contrast to conventional vehicles, battery recharging is significantly slower than refueling, which introduces additional temporal constraints and complicates route planning. As a result, routing decisions involve not only customer assignment and sequencing, but also battery feasibility, charging scheduling, and strict adherence to customer time windows.

The problem is defined on a fully connected symmetric graph $G=(N,L)$, where nodes represent the depot, customers, and charging stations. A single depot is assumed, serving as both the start and end location for all routes. Each link $(i,j)\in L$ is associated with a distance $dis{t}_{ij}=dis{t}_{ji}$, which determines the corresponding travel time, while energy consumption is assumed proportional to traveled distance.

To formally define the EVRPTW setting considered in this work, the following notation is used:

$N_c$: set of customer nodes;

$N_s$: set of charging-station nodes;

0: depot node;

$N=\left\{0\right\}\cup N_c\cup N_s$: set of all nodes;

$dis{t}_{ij}$: distance between nodes i and j;

$d_i$: nominal demand of customer i;

$s_i$: nominal service time of customer i;

$[a_i,b_i]$: service time window associated with customer i;

C: vehicle cargo capacity;

Q: vehicle battery capacity;

v: nominal vehicle speed;

r: energy-consumption rate;

g: charging rate;

$A_i$: arrival time at customer i;

$F_i$: completion time of service at customer i;

$T_i$: tardiness at customer i, defined as $T_i=max(F_i-b_i,0)$.

Each customer must be serviced exactly once by a single vehicle, prohibiting split deliveries. Service at customer $i\in N_c$ requires service duration $s_i$ and must be performed within the associated time window $[a_i,b_i]$, where $a_i$ denotes the earliest allowable service start time and $b_i$ denotes the due time. Vehicles arriving before $a_i$ must wait until service can begin. Depending on the formulation, time windows may be treated as either hard, where late arrivals are infeasible, or soft, where late service is permitted but penalized.

A homogeneous fleet of vehicles is available at the depot. Each vehicle has cargo capacity C, battery capacity Q, and travels at nominal speed v. Vehicles depart from the depot fully loaded and fully charged. Cargo capacity decreases as customers are served and cannot be replenished during route execution, while battery consumption is assumed proportional to traveled distance.

Charging stations allow vehicles to recharge their batteries during route execution. A full recharging strategy is assumed, meaning that each charging visit restores the battery to full capacity. Charging time is modeled linearly using charging rate g, such that recharging duration is proportional to the replenished energy amount. Charging stations are assumed to have unlimited capacity, allowing simultaneous charging of multiple vehicles.

Three optimization objectives are considered independently: minimization of fleet size, total energy consumption, and total tardiness, subject to customer-service, cargo-capacity, battery, charging, and time-window constraints.

The total tardiness objective is defined as

$$T_{tot}=\sum _{i\in N_c}{T}_i.$$

(1)

Total energy consumption is computed as

$$E_{tot}=\sum _{(i,j)\in L_{used}}r·dis{t}_{ij},$$

(2)

where $L_{used}$ denotes the set of traversed route links and r is the energy-consumption rate.

The fleet-size objective is defined as

$$V_{tot}=\left|V_{used}\right|,$$

(3)

where $V_{used}$ denotes the set of vehicles used to service all customers.

3.2. Stochastic Extension

In this work, a stochastic extension of the EVRPTW with soft time windows is considered. Uncertainty is introduced through multiplicative random factors applied independently to customer demand, service time, and vehicle speed, denoted by ${\lambda }^{\left(d\right)}$, ${\lambda }^{\left(s\right)}$, and ${\lambda }^{\left(v\right)}$, respectively. These sources of uncertainty follow established stochastic VRP formulations [14,15,31], where variability in demand, service duration, and travel time forms the standard modeling basis. Travel-time uncertainty is implemented indirectly via stochastic vehicle speed, preserving the network structure while inducing variability in traversal times.

The variability of each stochastic factor is controlled through a normalized variability parameter denoted by $CV$, while preserving the deterministic parameter as the expected value of the corresponding stochastic variable. Consequently, the deterministic instance represents the mean-case scenario, whereas increasing $CV$ introduces progressively stronger stochastic perturbations around the nominal values. In this study, values between $CV=0.1$ and $CV=0.3$ are used to represent moderate to pronounced uncertainty levels while maintaining stable and interpretable routing behavior.

Two distributions are considered in order to model qualitatively different uncertainty regimes. The first is a bounded symmetric uncertainty model implemented using a uniform distribution,

$$\lambda \sim \mathcal{U}[1-CV,1+CV],$$

(4)

which represents controlled fluctuations around the nominal value. The second is a positively skewed uncertainty model implemented using a lognormal distribution,

$$\sigma =\sqrt{ln(1+C{V}^2)},\mu =-{\displaystyle \frac{{\sigma }^2}{2}},\lambda \sim \mathrm{Lognormal}(\mu ,{\sigma }^2),$$

(5)

allowing strictly positive and asymmetric deviations with occasional larger realizations, better reflecting rare but substantial operational disruptions. Such parametric uncertainty models are commonly adopted in stochastic VRP and EVRP literature, where variability is typically represented through analytically defined probability distributions rather than real-world observational data [14,15,17].

The stochastic parameters are defined as

$${\tilde{d}}_i=d_i{\lambda }^{\left(d\right)},{\tilde{s}}_i=s_i{\lambda }^{\left(s\right)},\tilde{v}=v{\lambda }^{\left(v\right)},$$

(6)

where ${\lambda }^{\left(d\right)}$, ${\lambda }^{\left(s\right)}$, and ${\lambda }^{\left(v\right)}$ are independently sampled random variables representing demand, service-time, and speed variability, respectively. Correlated stochastic processes are not considered in this work, as the objective is to isolate the individual and combined effects of the principal uncertainty sources while maintaining a controlled and computationally tractable experimental setting.

Vehicle travel times are computed using the realized speed as

$${\tilde{t}}_{ij}={\displaystyle \frac{dis{t}_{ij}}{\tilde{v}}}={\displaystyle \frac{dis{t}_{ij}}{v{\lambda }^{\left(v\right)}}}.$$

(7)

A new realization of ${\lambda }^{\left(v\right)}$ is independently sampled for each vehicle movement, inducing stochastic travel times throughout route execution. In contrast, the realized customer demand ${\tilde{d}}_i$ and service time ${\tilde{s}}_i$ become known only upon arrival at customer i, reflecting the standard assumption of online information revelation in stochastic VRP models [33,34,35].

3.3. Routing Policy

Since this study considers a stochastic variant of the EVRPTW, a static precomputed routing plan is not an appropriate solution paradigm. In deterministic settings, routes can be constructed offline under the assumption that all input parameters are known in advance. In contrast, in the stochastic environment considered here, key parameters such as customer demand, service time, and vehicle speed are revealed only during route execution. Consequently, routes that appear feasible or near-optimal under nominal values may become suboptimal or infeasible once uncertainty is realized.

Instead of searching for a fixed routing solution, the objective is to develop a dynamically adaptive routing policy that makes sequential decisions online based on the current system state and realized stochastic information. This approach builds upon the routing policy framework introduced for the deterministic EVRPTW in [20], and extends it to the stochastic setting considered in this work. This perspective is consistent with evolutionary systems in nature, where adaptive behavior emerges from local decision rules interacting with a changing environment.

The proposed methodology consists of two components:

1.

Route Generation Scheme (RGS), which defines how routes are incrementally constructed;

2.

Priority Function (PF), which determines which customer should be selected next at each decision point.

Figure 1 provides an overview of the proposed framework, including the offline evolution of priority functions and the online adaptive route construction process under stochastic realizations. The route construction process defined by the RGS is detailed in the following section, while the representation and learning of the priority function are described in Section 3.5.

3.4. Route Generation Scheme

RGS defines the mechanism by which vehicle routes are constructed. In this study, five variants are considered: the Serial, Semi-parallel, and Parallel schemes adopted from [20], along with two extensions proposed in this work, namely the Semi-parallel-B and Parallel-B variants.

3.4.1. Common Feasibility and Operational Rules

Energy feasibility is verified before each movement decision. A vehicle may proceed to a destination only if it has sufficient energy to (i) reach the destination and (ii) subsequently reach the nearest charging station, preventing infeasible states without access to recharging. If this condition is not satisfied, the vehicle first travels to a feasible charging station, recharges, and then continues its route. Charging stations are selected by minimizing the total energy required to travel from the current location to the station and from the station to the intended destination.

Under stochastic demand, the realized demand of a selected customer may exceed the vehicle’s remaining capacity. In such cases, the service attempt is aborted, the vehicle returns to the depot (via a charging station if necessary), and the route is finalized. The customer remains unserved and is reconsidered later. This mechanism ensures feasibility under stochastic realizations and is applied implicitly in all RGS variants, although omitted from Algorithm 1 for clarity.

Time windows are modeled as soft. Under stochastic variations in service time and travel speed, enforcing hard time windows would frequently lead to infeasible routes, introducing additional feasibility-repair mechanisms that are orthogonal to the focus of this study. Moreover, strict hard time windows are often unrealistic in practical settings, where delays can occur and are typically tolerated to some extent. Therefore, late arrivals are allowed but penalized in the objective function.

In the Semi-parallel and Parallel variants, $LB$ denotes a capacity-based lower bound on the required number of vehicles, computed as

$$LB=⌈{\displaystyle \frac{{\sum }_{i\in N_c}{d}_i}{C}}⌉,$$

(8)

where $N_c$ is the set of customer nodes. This bound represents the minimum number of vehicles required to satisfy total nominal customer demand under cargo-capacity constraints.

3.4.2. Baseline RGS Variants

The Serial, Semi-parallel, and Parallel RGS variants adopted from [20] differ in how vehicles are activated and extended during route construction. In the Serial scheme, routes are constructed sequentially by activating one vehicle at a time. In contrast, the Semi-parallel and Parallel schemes extend multiple routes concurrently, initially activating $LB$ vehicles. The Semi-parallel variant selects the earliest available vehicle for extension and may switch to sequential construction if all active vehicles complete their routes. The Parallel variant maintains continuous parallelism by immediately activating a new vehicle whenever one completes its route. All variants follow the same route construction logic: routes are incrementally extended by selecting an unserved customer, verifying capacity and energy feasibility, and either continuing the route or returning the vehicle to the depot. The differences between variants lie solely in vehicle initialization and selection policies.

3.4.3. Extended RGS with Candidate Set Selection

The proposed Semi-parallel-B and Parallel-B variants extend their respective base schemes by modifying the vehicle selection mechanism. Instead of selecting only the earliest available vehicle, the set of the K earliest available vehicles is considered, and the vehicle with the largest remaining capacity is chosen. In this study, K is set to 3, providing a small but diverse candidate set that allows capacity-aware selection while keeping the decision process local and computationally efficient. This modification improves robustness under stochastic demand by reducing the likelihood of premature capacity saturation and unnecessary route termination.

3.4.4. Unified RGS Framework

All variants described above can be expressed within a single unified route construction framework. The variants differ only in vehicle initialization, selection, and update rules, while the underlying route construction logic remains identical. This unified procedure is presented in Algorithm 1, with variant-specific rules summarized in Table 1.

Algorithm 1 Unified RGS framework

Require:  variant $\in \{$Serial, Semi-parallel, Parallel, Semi-parallel-B, Parallel-B} 1: $UC\leftarrow N_c$▹ set of unserved customers 2: $Routes\leftarrow \varnothing$ 3: $V\leftarrow InitializeActiveVehicles(variant,LB)$ 4: $vehicle\leftarrow null$ 5: while $UC\ne \varnothing$ do 6:     if $vehicle=null$ then 7:         $vehicle\leftarrow SelectOrActivateVehicle(variant,V)$ 8:     end if 9:     $customer\leftarrow$ select a customer from $UC$ using the PF 10:     if $vehicle$ has sufficient remaining capacity for nominal demand of $customer$ then 11:         $destination\leftarrow customer$ 12:     else 13:         $destination\leftarrow 0$ 14:     end if 15:     if $vehicle$ does not have sufficient battery charge to reach $destination$ and then the nearest charging station in $N_s$ then 16:         $chargingStation\leftarrow$ select a charging station from $N_s$ 17:         move $vehicle$ to $chargingStation$ 18:         recharge $vehicle$ 19:     end if 20:     move $vehicle$ to $destination$ 21:     if $destination\in UC$ then 22:         serve $customer$ 23:         $UC\leftarrow UC\setminus \left\{customer\right\}$ 24:     else 25:         add the completed route of $vehicle$ to $Routes$ 26:         $V\leftarrow UpdateActiveVehicles(variant,V,vehicle)$ 27:         $vehicle\leftarrow null$ 28:     end if 29: end while 30: return $Routes$

3.5. Genetic Programming for Evolving Priority Functions

While the RGS defines how routes are constructed, the quality of the routing policy depends on the rule used to select the next customer at each decision point. In this work, this rule is represented by a priority function (PF), which assigns a numerical score to each unserved customer, and the customer with the highest score is selected.

The PF is automatically evolved using Genetic Programming (GP), which constructs decision rules by combining problem-specific features through mathematical operators [21,45]. Each GP individual represents a candidate PF encoded as an expression tree, where terminal nodes correspond to features of the current system state and internal nodes represent mathematical operators.

During route construction, the expression tree is evaluated for each candidate customer. Since some features depend on stochastic quantities (e.g., travel times or arrival times), a single evaluation may be sensitive to random realizations. To mitigate this effect, each candidate is evaluated multiple times under independent stochastic samples. In this work, the evaluation is repeated five times, providing a small number of samples that stabilizes the selection while keeping the computational cost manageable. The most frequently preferred customer is then selected. In the case of ties, preference is given to the customer selected by the smaller GP expression tree. This sampling-based decision mechanism improves robustness while maintaining computational efficiency.

The GP procedure follows a standard generational evolutionary scheme, where a population of candidate PFs is iteratively improved through selection, crossover, and mutation. The initial population consists of expression trees generated using the full method [21], with a depth of five. During evolution, the tree depth is limited to $2^8-1$, and individuals exceeding this limit are considered invalid and removed from the population. The evolutionary process runs for a fixed number of generations. The configuration largely follows standard GP settings, consistent with common practice [21] and with prior GP-based routing approaches [20]. The general workflow is shown in Algorithm 2, and the configuration parameters used in the experiments are summarized in Table 2. The implementation is based on the Jenetics1 evolutionary computation library (version 8.3.0).

Algorithm 2 Descriptive evolutionary workflow in Jenetics

1: Generate the initial population $P_0$ 2: Evaluate the fitness of all individuals in $P_0$ 3: while the termination criterion is not satisfied do 4:     Increment the generation counter: $g\leftarrow g+1$ 5:     Select the survivor population $S_g$ from the previous population $P_{g-1}$ 6:     Select the offspring population $O_g$ from the previous population $P_{g-1}$ 7:     Apply genetic alterations (crossover and mutation) to $O_g$ 8:     Remove invalid individuals from $S_g$ and $O_g$ 9:     Construct the new population $P_g$ by combining the filtered survivors and offspring 10:     Evaluate the fitness of all individuals in $P_g$ 11: end while

3.5.1. Terminal and Function Nodes

The primitive set builds upon the terminals and functions proposed in [20], which are retained to preserve the original EVRPTW routing policy framework. To address the stochastic setting considered in this work, the terminal set is extended with stochastic descriptors ($Var\_Dni$, $Var\_Tij$, $Var\_Sni$) intended to capture expected variability in customer demand, travel time, and service time, as well as global-state and coordination descriptors ($UC$, $DsumUC$, $CsumV$, $BestOtherET{A}_i$) intended to provide information about the remaining workload, available vehicle resources, and interactions between concurrently active vehicles. In addition, $Slack\_TW$ is included to provide information about the remaining time-window slack available for servicing a customer.

For objective-specific optimization, $CminV$ is introduced for the vehicle minimization objective to capture the minimum remaining capacity among active vehicles, while $SlackSelf\left(i\right)$ is introduced for tardiness minimization to capture the available temporal slack of the currently considered vehicle for customer i, computed as the difference between the customer due time and the estimated vehicle arrival time.

Due to the online nature of the stochastic setting, realized stochastic values are not available during customer selection. Consequently, all terminal nodes whose computation depends on stochastic quantities are evaluated using nominal parameter values. The complete set of terminal nodes is summarized in Table 3, while their relevance and usage patterns within evolved routing policies are further analyzed in Section 4.6.

The function nodes used in the GP representation are summarized in Table 4. The first row contains binary operators, and the second row contains unary operators. To ensure numerical stability during tree evaluation, the operators div, log, and sqr are implemented as safe functions. Specifically, div returns 0 when the denominator is close to zero, while log and sqr return 0 for non-positive arguments.

4. Results

4.1. Experimental Protocol

The experiments were conducted on a subset of Schneider’s benchmark instances [9] for the EVRPTW. In total, 48 instances were used, each consisting of 100 customers, 21 charging stations, and one depot. The focus on 100-customer instances was motivated by the stochastic setting considered in this work, as smaller instances tend to exhibit high variability and less stable routing behavior under random perturbations. According to the spatial distribution of customers, the instances are categorized into Clustered (C), Random (R), and Random-Clustered (RC) sets.

The dataset was split into mutually disjoint training and test sets. The training set contains 30 instances (10 C, 10 RC, and 10 R), while the test set contains the remaining 18 instances (6 C, 6 RC, and 6 R). The split preserves the balance across instance types and allocates a larger portion of instances to the test set to enable a more robust generalization evaluation. All results reported in this section are obtained on the test set, ensuring evaluation on previously unseen instances.

During training, each instance was evaluated twice per fitness computation in order to account for stochastic variability. All GP models were evolved under a fixed stochastic setting in which customer demand, service time, and vehicle speed follow a lognormal distribution with variability parameter $CV=0.2$. To account for both stochasticity and the randomness of GP initialization, the entire training procedure was repeated 10 times, resulting in 10 independently evolved PFs per optimization objective. This provides a representative sample of independently evolved PFs while keeping the computational cost manageable.

Each evolved PF was subsequently evaluated on a range of stochastic test scenarios. These scenarios systematically vary the distribution type and the level of variability for demand, service time, and vehicle speed. A test scenario is denoted as

$$\left\{\mathtt{distribution}\right\}-\left\{{\mathrm{CV}}_d\right\},\left\{{\mathrm{CV}}_s\right\},\left\{{\mathrm{CV}}_v\right\},$$

where distribution specifies the probability distribution (DET for deterministic, LN for lognormal, and U for uniform), and $C{V}_d$, $C{V}_s$, and $C{V}_v$ denote the variability parameters associated with customer demand, service time, and vehicle speed, respectively.

The deterministic scenario DET-0,0,0 corresponds to the standard EVRPTW without stochastic perturbations and serves as an internal baseline, enabling comparison with the original deterministic GP-RGS framework [20] under identical evaluation conditions and isolating the effect of incorporating stochastic information into the routing policy. This abbreviated notation is used consistently throughout the result tables.

The following evaluation scenarios were considered:

• DET-0,0,0
• LN-0.1,0,0, LN-0.2,0,0, LN-0.3,0,0
• LN-0,0.1,0, LN-0,0.2,0, LN-0,0.3,0
• LN-0,0,0.1, LN-0,0,0.2, LN-0,0,0.3
• LN-0.2,0.2,0, LN-0.2,0,0.2, LN-0,0.2,0.2
• LN-0.2,0.2,0.2, LN-0.3,0.3,0.3
• U-0.2,0.2,0.2, U-0.3,0.3,0.3

For each test scenario, every PF is evaluated on all 18 test instances (6 C, 6 RC, and 6 R). To account for stochastic variability, each instance is evaluated six times with independently sampled realizations of demand, service time, and travel time, resulting in 108 runs per scenario. For a given PF and scenario, the objective values are summed over all 108 runs, yielding a single score per PF. The reported results summarize the distribution of these scores across the 10 independently evolved PFs, both through descriptive statistics (min, avg, max) and through the full distributions shown in the violin plots. This evaluation procedure is applied consistently across all optimization objectives. The average online route-construction time on the considered 100-customer instances was approximately 10 ms per instance. To support reproducibility, the implementation used in this study is publicly available at: https://github.com/AlanDurdevic/stohastic-EVRP-GP.

4.2. Results for Vehicle Minimization Objective

Table 5 summarizes the results obtained when minimizing the number of vehicles. Figure 2 provides two complementary views of these results. The line plot (Figure 2) shows the average number of vehicles across all test scenarios and serves as the primary basis for comparing the methods. The results reveal a highly stable overall pattern. Semi-parallel and Semi-parallel-B consistently achieve the lowest average values, with only minor differences between them across scenarios. Serial consistently performs slightly worse than these two variants, but remains clearly better than the Parallel-based strategies. Among the latter, Parallel-B improves upon Parallel, although both remain substantially inferior to the Serial and Semi-parallel variants. This ordering is preserved across all scenarios, including the deterministic case (DET-0,0,0), indicating that stochasticity has little effect on the relative ranking of the RGS variants, although it does influence absolute performance levels.

The violin plots (Figure 2) complement this view by showing the full distribution of results across all runs. These distributions reinforce the separation between the methods: Semi-parallel and Semi-parallel-B are concentrated at the lowest values, Serial occupies a slightly higher but still compact range, while Parallel and Parallel-B are shifted to markedly higher values and exhibit much broader variability. This indicates that the Parallel-based strategies are not only worse on average, but also less stable under stochastic conditions. For completeness, detailed violin plots for each individual test scenario are provided in Appendix A.

The effect of the proposed B-modification varies across routing strategies. For the Parallel variant, Parallel-B consistently improves performance, leading to lower average values and reduced variability. This is reflected in both the average trends and the distributional view, where Parallel shows the highest values and widest spread, while Parallel-B shifts toward lower and more concentrated results. For the Semi-parallel strategy, the impact is less pronounced. Semi-parallel and Semi-parallel-B achieve very similar average results, with only minor differences across scenarios. However, the minimum values in Table 5 indicate that Semi-parallel-B more frequently attains the lowest minima, suggesting occasional improvements that are not consistently reflected in the average performance.

Regarding the influence of stochasticity, demand variability has the most pronounced effect on fleet size. As the variability parameter of customer demand increases, the average number of vehicles also increases across all methods, with the effect being particularly strong for the Parallel-based variants, which also exhibit increased variability. In contrast, variability in service time has little to no impact on fleet size, with both average values and distributions remaining largely unchanged across scenarios. Variability in vehicle speed has a moderate effect. While the overall ranking of methods remains unchanged, speed variability leads to a slight increase in dispersion and, in some cases, small shifts in average performance.

4.3. Results for Energy Consumption Minimization Objective

Table 6 summarizes the results obtained when minimizing energy consumption, while Figure 3 shows the corresponding average values across test scenarios and distributions. The average values (Figure 3) reveal a clear and stable separation for the best-performing method: the Serial variant consistently achieves the lowest energy consumption in all scenarios. Notably, this behavior is also observed in the deterministic scenario (DET-0,0,0), indicating that the superiority of the Serial strategy is not driven by stochastic effects, but rather reflects its inherently more energy-efficient routing structure. Beyond this, the ordering of the remaining methods is less consistent. Semi-parallel and Semi-parallel-B generally outperform the Parallel-based strategies, but their relative ranking varies across scenarios. Similarly, Parallel-B most often yields the highest energy consumption, although not in every instance. This indicates that, unlike the vehicle minimization objective, energy consumption is more sensitive to the interaction between routing strategy and stochastic variability, leading to scenario-dependent performance differences among the non-serial variants.

The influence of stochasticity depends on the source of uncertainty. Demand variability has only a minor effect, producing small changes in average values without altering the overall structure of the results. In contrast, variability in service time has a more structural impact: while its effect on average energy consumption is moderate, it can alter the relative performance of closely competing methods. In particular, the ordering between Semi-parallel and Semi-parallel-B is occasionally reversed under service-time variability. Variability in vehicle speed has a weaker influence, introducing modest shifts in average performance and increased dispersion without changing the overall ranking. Additional insight is provided by the minimum values reported in Table 6. Although the B-variants often exhibit worse average performance, they consistently achieve lower minimum values than their non-B counterparts, indicating that more energy-efficient solutions can occasionally be obtained.

The distributional view (Figure 3) supports these observations. Serial remains tightly concentrated at low values, indicating both strong performance and high stability. Semi-parallel and Semi-parallel-B exhibit a relatively compact spread, while the Parallel-based variants are shifted to higher values and display substantially greater dispersion. This is particularly evident for Parallel-B, which shows the widest spread and long upper tails, reflecting occasional highly inefficient solutions.

4.4. Results for Total Tardiness Minimization Objective

Table 7 summarizes the results for the tardiness objective, while Figure 4 presents the corresponding averages and distributions. In contrast to the energy objective, the performance hierarchy is clearly inverted. As shown in Figure 4, Parallel-based strategies achieve the lowest tardiness across all scenarios, with Parallel-B consistently achieving the lowest average values. The remaining methods exhibit more variability: Parallel generally performs well but is less stable, Semi-parallel variants occupy an intermediate range, while Serial performs worst by a substantial margin. This indicates that increased parallelism can be beneficial for handling time-window constraints, although its effectiveness depends on the specific routing strategy. This behavior is also observed in the deterministic scenario (DET-0,0,0), indicating that the advantage of Parallel-based strategies is not driven by stochastic effects but reflects their ability to better handle time-window constraints.

The distributional view (Figure 4) further supports these observations. Parallel-B achieves low average tardiness with relatively low variability, while Parallel exhibits comparable performance with somewhat higher dispersion. Parallel performs similarly, although with higher dispersion. In contrast, the Serial strategy shows both substantially higher average tardiness and a wide spread with long upper tails, indicating frequent occurrences of highly delayed solutions. Semi-parallel strategies occupy an intermediate position, with moderate dispersion and performance. Overall, these results suggest that higher degrees of parallelism tend to reduce tardiness and improve robustness, although differences between closely related variants are less pronounced.

The effect of the B-modification is most noticeable for the Parallel strategy, where Parallel-B generally achieves lower average tardiness and reduced variability. For the Semi-parallel strategy, the effect is less consistent, with Semi-parallel and Semi-parallel-B exhibiting similar performance across scenarios.

The influence of stochasticity does not follow a simple or monotonic pattern. Variations in individual factors often have limited or inconsistent effects on their own, while combinations of stochastic elements lead to abrupt changes in performance. In particular, changes in demand alone have little impact on tardiness, whereas variability in service time and vehicle speed can trigger sharp increases for certain strategies, most notably the Parallel variant. However, these effects are not cumulative in a predictable way: in scenarios where all three sources of uncertainty are present at moderate levels, the methods often stabilize and even show slight improvements in average performance.

4.5. Pairwise Comparison of Routing Strategies

While the previous analysis focuses on aggregated performance through averages and distributions, this section examines the consistency of pairwise differences between routing strategies across identical problem realizations. By comparing methods on a per-instance and per-run basis, the analysis isolates structural differences in decision behavior that may be obscured by aggregation.

Since all methods are evaluated on the same instances and stochastic realizations, observations are paired. Accordingly, statistical comparisons are performed using the Friedman test followed by pairwise Wilcoxon signed-rank tests with Holm correction. The Friedman test confirms the presence of statistically significant differences among routing strategies for all objectives (vehicles: ${\chi }^2=3199.65$, energy: ${\chi }^2=4744.43$, tardiness: ${\chi }^2=3321.88$; all $p<0.001$), justifying subsequent pairwise analysis.

The results of the pairwise comparisons are summarized in Table 8. Rows and columns correspond to routing variants, and each cell reports the outcome of the Wilcoxon test with Holm correction. A checkmark () indicates a statistically significant difference ($p<0.05$), while exact p-values are shown only for non-significant comparisons. For completeness, the median difference $\tilde{\Delta }$ (row−column) is also provided; negative values indicate that the row method tends to achieve lower objective values than the column method.

For the vehicle objective, all pairwise differences are statistically significant. However, several comparisons exhibit zero median difference, indicating that the observed differences, although highly consistent across realizations, are practically small. This highlights that statistical significance in this setting primarily reflects consistency rather than effect magnitude.

For the energy objective, all pairwise differences are also statistically significant, reflecting a clear and consistent ordering of the methods. In this case, statistical significance aligns closely with practical relevance, as the observed median differences are substantial and consistent with the separation already visible in the aggregated results. This confirms that the dominance of the Serial strategy and the inferior performance of Parallel-based variants are not only evident on average, but persist across individual realizations.

For the tardiness objective, the results exhibit a different pattern. Most pairwise differences are statistically significant; however, the comparison between Parallel and Parallel-B is not significant ($p=0.578$), indicating that these two strategies achieve comparable performance when evaluated on identical realizations. This contrasts with the distributional analysis, where differences in variability are apparent, and suggests that the distinction between these methods is driven more by dispersion and extreme outcomes than by consistent differences in central tendency. Overall, these findings reinforce that the effect of parallelism depends strongly on the objective, and that paired comparisons provide additional insight into the consistency of method behavior beyond aggregated performance measures.

4.6. Analysis of Terminal Node Usage

To gain insight into the decision mechanisms learned by the GP-based routing policies, we analyze terminal node usage across routing strategies and objectives. Specifically, we examine (i) the relative frequency of terminal nodes within evolved policies, grouped by routing strategy (RGS), and (ii) the association between terminal usage and objective values. Terminal frequencies are visualized using heatmaps (Figure 5), enabling identification of dominant features and structural differences between strategies. To assess functional relevance, we compute Spearman rank correlations between terminal usage and objective values (Figure 6). Only correlations with $p<0.1$ are reported, with stronger signals ($p<0.05$) highlighted, in order to retain consistent trends while accounting for the limited sample size per configuration.

For the vehicle minimization objective (Figure 5), Semi-parallel and Semi-parallel-B strategies are dominated by $D_{ni}$ and $E_{ni}$, each accounting for nearly 50% of node usage. Parallel variants rely on the same features but exhibit a more distributed pattern, while the Serial strategy shows the most dispersed usage without a clearly dominant feature.

The correlation analysis (Figure 6) is consistent with the observed usage patterns. The feature $D_{ni}$ exhibits a consistent negative correlation with the number of vehicles across semi-parallel and parallel strategies, indicating that prioritizing demand is associated with improved performance. In contrast, despite its high frequency, $E_{ni}$ shows weak correlation with the objective, suggesting a supporting rather than driving role. The Serial strategy follows a different pattern, with stronger correlations for timing- and state-related features (e.g., $R{T}_{ni}$, $E_{vk}$), while demand-related signals remain less influential.

The terminal node distribution for the energy minimization objective (Figure 5) is dominated by $E_{ni}$ across all strategies, particularly in semi-parallel variants. Beyond this common pattern, strategies differ in how additional features are used. The Serial strategy relies on a combination of energy-, demand-, and routing-related features, including $R{T}_{ni}$, $D{D}_{ni}$, and $EDe{p}_{ni}$, all of which appear with notable frequency alongside $E_{ni}$. Semi-parallel variants remain strongly concentrated on $E_{ni}$, with relatively limited use of other terminals. In contrast, the Parallel strategy exhibits a more distributed pattern across multiple terminals, including additional energy-related and stochastic descriptors such as $Va{r}_{Dni}$, $Va{r}_{Tij}$, and $Va{r}_{Sni}$. The Parallel-B variant follows a more constrained pattern, relying on fewer additional features compared to Parallel.

A clear pattern emerges in the correlation results (Figure 6): each routing strategy is primarily associated with a single dominant feature. The Serial strategy shows a strong negative correlation with $E{C}_{ni}$, while semi-parallel variants exhibit a positive correlation with the demand-related feature $D{D}_{ni}$. In contrast, the Parallel strategy shows significant correlations for multiple features, most notably $Va{r}_{Sni}$, indicating a broader set of influencing features. The Parallel-B variant exhibits a negative correlation with $EDe{p}_{ni}$, but without additional strongly correlated features.

For the tardiness minimization objective (Figure 5), $E_{ni}$ remains the most frequently used feature across all strategies, but is accompanied by a broader use of additional terminals, particularly in Parallel and Parallel-B. These strategies show increased utilization of customer- and timing-related features (e.g., $D{D}_{ni}$, $R{T}_{ni}$, $S{T}_{ni}$, $Tvk$, $Slac{k}_{TW}$), as well as system-level descriptors such as $UC$, indicating a more diverse decision basis combining feasibility, timing, and urgency information. In contrast, the Serial strategy distributes importance across multiple features without a clearly dominant signal, while the Semi-parallel strategy remains largely centered on $E_{ni}$ with limited contribution from timing- or urgency-related features.

The correlation results (Figure 6) show that $UC$, one of the newly introduced system-level descriptors, exhibits the strongest positive correlation with the objective in the Parallel strategy. In contrast, $E_{ni}$, although the most frequently used terminal, does not show consistently strong correlation. In Parallel-B, multiple features (including $R{T}_{ni}$, $E_{ni}$, and $Va{r}_{Tij}$) exhibit significant correlations, indicating a more distributed pattern. Other strategies are typically associated with a single dominant correlated feature (e.g., $UC$ or $E_{vk}$ in Serial, $Va{r}_{Sni}$ in Semi-parallel), although these signals are generally weaker or less consistent.

4.7. Comparison with Greedy Baselines

To provide additional context regarding the effectiveness of the proposed GP-based routing policies, we compare them against several simple constructive heuristics commonly used in routing settings:

• Nearest Neighbor (NN): selects the geographically closest feasible customer.
• Minimum Travel Energy (MTE): selects the customer with the lowest feasible travel-energy cost.
• Minimum Slack (MS): selects the customer with the smallest remaining time-window slack.
• Earliest Due Time (EDT): selects the customer with the earliest due time.

All greedy heuristics were implemented within the same route-construction framework as the proposed approach, replacing the GP-evolved priority function with the corresponding hand-designed rule. To isolate the effect of customer-selection policies, all greedy baselines use the Serial RGS variant. Since each heuristic corresponds to a single deterministic policy, only aggregated objective values are reported. The results are shown in Figure 7.

The nearest-neighbor heuristic remains competitive for fleet-size minimization in deterministic and low-variability scenarios, occasionally achieving slightly lower values than the GP-based Semi-parallel policy. However, its performance deteriorates more rapidly as demand variability increases, while the GP-based approaches remain more stable across stochastic settings.

For the energy objective, the GP-based Serial policy consistently achieves the lowest consumption across all tested scenarios. Although the minimum-travel-energy heuristic remains competitive, the remaining hand-designed heuristics perform substantially worse.

The largest performance differences are observed for tardiness, where the GP-based Parallel-B strategy consistently outperforms all greedy baselines, suggesting that the evolved policies better balance routing, scheduling, and charging decisions under stochastic operating conditions.

5. Discussion

5.1. Vehicle Minimization

The results largely follow the main conclusions reported in [20], where Serial and Semi-parallel strategies outperform the Parallel approach in terms of fleet size. However, under stochastic conditions the balance shifts more clearly toward limited parallelism, with Semi-parallel and Semi-parallel-B consistently achieving better results than Serial across all scenarios. Fully sequential construction delays the activation of additional vehicles, which becomes problematic when stochastic realizations disrupt routes, while fully parallel construction fragments the solution too early and reduces capacity utilization efficiency. Semi-parallel approaches maintain enough structure to avoid fragmentation while still reacting to emerging constraints. The B-modification further reinforces this behavior by prioritizing vehicles with larger remaining capacity and reducing the risk of premature capacity saturation. Among the considered stochastic factors, demand variability has the strongest influence on fleet size, directly affecting capacity feasibility and often forcing the use of additional vehicles, whereas service-time variability has relatively limited effect and travel-time variability mainly increases dispersion without changing the relative ordering of routing strategies. Although the absolute differences between methods remain relatively small, the statistical analysis confirms that these patterns are highly consistent across stochastic scenarios.

The node analysis further supports these observations. Vehicle minimization is primarily associated with the combined use of $D_{ni}$ and $E_{ni}$, particularly in the Semi-parallel and Semi-parallel-B variants. The negative correlation of $D_{ni}$ with the objective indicates that prioritizing demand contributes to reducing fleet size, while the high usage but weak correlation of $E_{ni}$ suggests a supporting feasibility-related role. In contrast, the Serial strategy relies more heavily on timing- and state-related features and exhibits weaker demand-related signals, consistent with its less effective handling of capacity utilization. Stochastic descriptors show limited influence in both usage and correlation analyses, suggesting that vehicle minimization is driven primarily by stable structural routing behavior rather than explicit modeling of variability.

5.2. Energy Minimization

For the energy minimization objective, the relative behavior of routing strategies changes substantially. Consistent with [20], the Serial strategy achieves the best overall performance. Unlike vehicle minimization, where limited parallelism is beneficial, energy minimization favors sequential route construction, which produces geographically compact routes. Because of this structure, the Serial strategy remains comparatively stable even when stochastic travel-time variability increases. In contrast, Parallel-based strategies activate multiple routes early, producing longer, less coherent routes with higher energy consumption. Under stronger stochastic variability, this leads to dispersed and occasional highly inefficient solutions, particularly for Parallel-B, where prioritizing vehicles with larger remaining capacity further reduces route compactness. Semi-parallel variants maintain a more balanced trade-off and therefore remain considerably more stable across stochastic scenarios.

The node analysis further clarifies these differences. In the Serial strategy, energy minimization is associated with a strong negative correlation of $E{C}_{ni}$, indicating that decisions are guided by a feature directly related to energy cost. In contrast, semi-parallel strategies mostly use $E_{ni}$ node but their correlation is driven by demand-related features such as $D{D}_{ni}$, suggesting a weaker alignment between dominant signals and the objective. Parallel variants do not exhibit a clearly dominant energy-related decision signal, consistent with their weaker performance for the energy objective.

5.3. Tardiness Minimization

Unlike energy minimization, which favors sequential route construction, tardiness minimization benefits from early and distributed route activation. As also observed in [20], Parallel-based strategies achieve the best overall performance by reducing the accumulation of delays within individual routes. In contrast, the Serial strategy delays vehicle activation, making recovery from stochastic disruptions more difficult. In highly parallel settings, the B-modification further improves average tardiness and solution stability by prioritizing vehicles with larger remaining capacity among the earliest available candidates. However, the lack of statistically significant differences between Parallel and Parallel-B suggests that the primary advantage lies in the degree of parallelism itself rather than in the specific vehicle-selection variant. Unlike the energy objective, the influence of stochasticity on tardiness is less structured and does not follow simple monotonic trends. Individual sources of variability often produce limited effects, while specific combinations can trigger abrupt performance changes and occasional spikes, indicating stronger interactions between stochastic components.

The node analysis suggests that each routing strategy tends to rely on a small number of dominant signals. In the Parallel strategy, the strongest correlation is observed for the system-level descriptor $UC$, suggesting that global information regarding the number of remaining unserved customers becomes particularly relevant in highly parallel routing settings, where coordination between concurrently active vehicles plays a larger role. Parallel-B exhibits a more distributed set of significant correlations, including $R{T}_{ni}$, $E_{ni}$, and $Va{r}_{Tij}$, while other strategies are typically associated with only one dominant correlated feature.

5.4. Overall Observations

Taken together, the results indicate that no single routing strategy is universally optimal across objectives. Instead, the effectiveness of route construction depends strongly on the interaction between the optimization objective, the degree of parallelism, and the underlying stochastic conditions. Vehicle minimization benefits from limited parallelism and balanced capacity utilization, energy minimization favors sequential and spatially coherent routing, while tardiness minimization benefits from early workload distribution across multiple vehicles. Across all objectives, stochasticity primarily amplifies these structural differences rather than fundamentally altering the relative behavior of routing strategies. The node analysis further suggests that effective routing policies emerge through different combinations of dominant decision signals depending on the objective and route-construction scheme.

6. Conclusions

This paper presented an extension of a GP-based RGS framework to the stochastic electric vehicle routing problem with time windows. The proposed approach incorporates uncertainty in customer demand, service time, and travel time through controlled stochastic perturbations, while extending the GP priority function with stochastic and global-state descriptors. In addition, two new route-generation variants, Semi-parallel-B and Parallel-B, were introduced to improve capacity-aware vehicle selection under uncertainty.

The results show that routing strategy performance depends strongly on the optimization objective. Semi-parallel strategies achieved the best fleet-size performance by balancing capacity utilization and routing flexibility, while the Serial strategy consistently achieved the lowest energy consumption by preserving spatial coherence and compact routes. In contrast, Parallel-based strategies achieved the best tardiness performance by distributing workload earlier across multiple vehicles and reducing delay accumulation. Across objectives, stochasticity primarily amplified the structural differences between routing strategies rather than altering their relative behavior. The node analysis further showed that effective routing policies rely on a relatively small number of objective-relevant decision signals, with the newly introduced system-level descriptor $UC$ exhibiting a particularly strong association with tardiness minimization in highly parallel routing settings.

Additional experiments with hand-designed constructive heuristics further showed that the proposed GP-based policies remain competitive across stochastic scenarios and substantially outperform simple greedy baselines for the tardiness objective. Overall, the results support the suitability of lightweight evolutionary hyper-heuristics for constructing adaptive routing policies in stochastic EVRPTW environments.

The present study focuses on policy-based constructive heuristics under controlled stochastic settings rather than on large-scale stochastic optimization approaches. While comparisons with stochastic programming or reinforcement learning methods would provide additional perspective, such approaches operate under substantially different modeling and computational assumptions. Furthermore, uncertainty was modeled using parametric distributions in a static stochastic setting, enabling systematic and reproducible evaluation but not fully capturing the complexity of real-world logistics systems.

In addition, charging-station locations were treated as fixed exogenous components of the benchmark instances, and the robustness of the proposed routing policies under varying charging-network configurations was not investigated. Future work will therefore focus on extending the proposed framework to dynamic and real-time routing environments, incorporating richer data-driven uncertainty models, evaluating the approach on real-world operational datasets, and analyzing routing-policy behavior under alternative charging-infrastructure layouts and configurations.