JOURNAL TITLE
1
Bilevel Late Acceptance Hill Climbing for the Electric Capacitated Vehicle Routing Problem
arXiv:2604.13013v1 [cs.AI] 14 Apr 2026
Yinghao Qin, Mosab Bazargani, Edmund K. Burke, Carlos A. Coello Coello, Zhongmin Song, and Jun Chen
Abstract—This paper tackles the Electric Capacitated Vehicle Routing Problem (E-CVRP) through a bilevel optimization framework that handles routing and charging decisions separately or jointly depending on the search stage. By analyzing their interaction, we introduce a surrogate objective at the upper level to guide the search and accelerate convergence. A bilevel Late Acceptance Hill Climbing algorithm (b-LAHC) is introduced that operates through three phases: greedy descent, neighborhood exploration, and final solution refinement. b-LAHC operates with fixed parameters, eliminating the need for complex adaptation while remaining lightweight and effective. Extensive experiments on the IEEE WCCI-2020 benchmark show that b-LAHC achieves superior or competitive performance against eight state-of-theart algorithms. Under a fixed evaluation budget, it attains nearoptimal solutions on small-scale instances and sets 9/10 new best-known results on large-scale benchmarks, improving existing records by an average of 1.07%. Moreover, the strong correlation (though not universal) observed between the surrogate objective and the complete cost justifies the use of the surrogate objective while still necessitating a joint solution of both levels, thereby validating the effectiveness of the proposed bilevel framework and highlighting its potential for efficiently solving large-scale routing problems with a hierarchical structure. Index Terms—Metaheuristics, bilevel optimization, late acceptance hill climbing, E-CVRP.
I. I NTRODUCTION S the global shift toward sustainable transportation accelerates, electric vehicles (EVs) have become central to decarbonizing road-based logistics. In this context, the Electric Vehicle Routing Problem (EVRP) has emerged as a key research frontier. Its core variant, the Electric Capacitated Vehicle Routing Problem (E-CVRP) [1], forms the foundation of EV logistics. Extensions incorporate features such as non-linear charging, heterogeneous fleets, and time-window constraints. Methodological advances in solving E-CVRP are therefore crucial, underpinning applications from urban lastmile delivery to long-haul freight.
A
Y. Qin and J. Chen are with the Centre for Intelligent Transport, QMUL, London, U.K. (e-mail: [email protected]; [email protected]). M. Bazargani and E. K. Burke are with the School of Computer Science and Engineering, Bangor University, Bangor, U.K. (e-mail: [email protected]; [email protected]). C. A. Coello Coello is with the Department of Computer Science, Evolutionary Computation Group, CINVESTAV-IPN, Mexico City, Mexico (e-mail: [email protected]). He is also (as part of a sabbatical leave) with the Basque Center for Applied Mathematics & Ikerbasque, Spain. Z. Song is with the School of Future Technology, Shanghai University, Shanghai, China (e-mail: [email protected]). The first author is supported by the China Scholarship Council (Grant No. 202209110001). Manuscript received April 19, 2005; revised August 26, 2015. (Corresponding author: Jun Chen.)
9
9 9
9
9 9
Õ
9
9 9
9
Fig. 1. Illustrative solution to the E-CVRP, featuring one depot, four charging stations, and ten customer locations.
Given a fleet of homogeneous EVs, the objective of ECVRP is to determine the most efficient route for each vehicle that respects battery charge and vehicle cargo capacity constraints. Each vehicle must begin and end at a depot while serving a set of customers. Fig. 1 illustrates a schematic example of a feasible E-CVRP solution. In this example, the first route (solid lines) requires no charging stops. The second route (dashed lines) visits the same charging station twice. The third route (dotted lines) makes consecutive stops at two distinct charging stations before reaching the next customer, representing an extreme case of ultra-long-distance delivery. The E-CVRP can be viewed as a combination of the classical Capacitated Vehicle Routing Problem (CVRP) [2] and the Green Vehicle Routing Problem (GVRP) [3]. The former considers cargo capacity constraints, whereas the latter incorporates limited driving range under the assumption of unlimited cargo capacity. The E-CVRP simultaneously accounts for both vehicle capacity and driving range constraints. It is also a member of the family of Vehicle Routing Problems with Intermediate Stops (VRPIS) [4], where optional en-route stops are introduced to support service operations. VRPIS generally arise in three application contexts: replenishment and disposal operations, rest or idling periods, and refueling. The E-CVRP belongs to the refueling category. Current methodologies for solving the E-CVRP face two major challenges. First, at the modeling level, existing formulations either produce intractably large search spaces [1] or rely on aggressive pruning strategies [5] that risk discarding high-quality solutions. They often fail to scale effectively to large-scale instances. Second, from an algorithmic perspective, most existing approaches rely on highly customized,
JOURNAL TITLE
2
problem-specific heuristics, often embedded within vague con- space. At the upper level, a CVRP is solved, where customers ceptual frameworks [6], [7], [8], [9]. While effective, such are partitioned into feasible subsets and the visiting sequence is heuristics are difficult to reproduce, tune, and adapt. These lim- determined for each subset. At the lower level, a fixed-route itations motivate a modeling approach that reduces computa- vehicle charging problem (FRVCP) [11] optimizes charging tional complexity, together with a transparent, general-purpose station insertions along the given path from the upper level. solution approach capable of scaling to large instances. This bilevel structure substantially reduces the computational At the root, these difficulties stem from domain charac- complexity and provides a theoretical framework for integrated teristics specific to electric fleets. Compared to the classical routing and charging decision making. CVRP, the main additional challenge in E-CVRP arises from We further develop a bilevel Late Acceptance Hill Climbing the limited cruising range imposed by EV battery capac- (b-LAHC) algorithm for the E-CVRP. LAHC is a single-point ity. Unlike internal combustion engine vehicles, which can metaheuristic [12] that has shown success in various schedulbe refueled within minutes, EVs require significantly longer ing and combinatorial optimization problems [13], [14]. Its charging times, making detours for recharging costly in terms core mechanism, the “late acceptance” rule [12], is simple yet of operational time and scheduling. In real-world scenarios, effective: a candidate solution is accepted if it improves upon EVs may even need to interrupt their routes for multiple the current solution or outperforms the one encountered several recharges within a single trip. From the charging infrastructure iterations earlier. This mechanism introduces a well-balanced perspective, the limited number of stations implies that each trade-off between exploration and exploitation. of them may be visited once, multiple times, or not at all, The main contributions of this paper are the following: possibly by the same or different vehicles. These interdependencies between vehicle routing, battery management, and 1) A bilevel optimization model is explicitly adopted for the charging infrastructure create unique modeling challenges and E-CVRP for the first time. The model captures both routsignificantly increase the computational complexity in solving ing and charging decisions without duplicating nodes for E-CVRP at a large scale. charging stations. The solution and objective spaces are Most existing formulations model the E-CVRP as a singlefurther analyzed to guide efficient search strategies. Notably, level mixed-integer linear programming (MILP). A reprethe model enables partial evaluation of solutions without sentative example is provided by Mavrovouniotis et al. [1], solving the lower level, which serves as a surrogate for full where the problem is defined on a simple graph and binary evaluation and accelerates the search. decision variables are associated with each arc to indicate 2) A bilevel metaheuristic, b-LAHC, is proposed based on whether it is traversed by an EV. In order to accommodate late acceptance hill climbing. It features three phases: (i) multiple visits to charging stations within a route, charging initialization, where a greedy descent rapidly drives the nodes must be replicated. In some models, the number of upper-level routing solution to a local optimum; (ii) exploreplications is bounded by twice the number of customers [10], ration, a neighborhood search guided by the late accepcorresponding to a worst-case scenario in which a vehicle tance mechanism at the upper level, with efficient charging visits a charging station both before and after serving each decisions conditionally triggered at the lower level; and customer. Although mathematically rigorous, such a node (iii) refinement, where exhaustive charging optimization is replication approach dramatically enlarges the search space applied to the best-found solution at termination. This threeand makes the problem more difficult to solve. To overcome phase design balances computational efficiency and solution these drawbacks, Froger et al. [5] propose an alternative quality, while inheriting the convergence properties of the formulation based on pre-enumerating feasible charging staoriginal LAHC and adapting them to the bilevel structure. tion paths (CSPs) between each pair of non-charging nodes 3) An extensive study is conducted on the IEEE WCCI 2020 (customers or depots). The E-CVRP is then represented as E-CVRP benchmark [1], comparing b-LAHC with seven a multigraph, where multiple arcs between a pair of nonstate-of-the-art algorithms across 17 instances under two charging nodes correspond to different feasible CSPs. The termination criteria: a fixed number of evaluations and a model then selects appropriate CSPs to connect non-charging fixed runtime. Overall, b-LAHC achieves 10 new bestnodes and thereby construct complete EV routes. However, known solutions (BKS), matches 2 and remains competitive the number of feasible CSPs grows rapidly with instance on the rest. Its strength is particularly pronounced on largesize. Even with strong dominance pruning, hundreds of paths scale instances, improving the current BKS by an average are generated for 10-customer instances and over a thousand of 1.07% across 10 cases, and reducing solution standard for 20-customer instances on average [5]. Consequently, the deviation by 46.2%, relative to the second-best method. resulting MILP can be solved using a commercial solver only for small-scale instances, while most 20-customer cases The remainder of this paper is organized as follows. Secalready become computationally intractable. tion II reviews the literature on bilevel optimization in routing In light of these limitations of single-level formulations, problems and summarizes the current approaches for EVRPs. we adopt a bilevel modeling perspective for the E-CVRP. Section III presents the proposed bilevel mathematical model The proposed formulation decomposes the problem into two of the E-CVRP. Section IV introduces the b-LAHC algorithm. interdependent components: routing and charging. We analyze Section V reports our computational experiments and analyzes their interactions and leverage these insights to guide algorithm the results in detail. Section VI concludes the paper and design, focusing on the most promising regions of the search outlines potential directions for future research.
JOURNAL TITLE
3
TABLE I S UMMARY OF NOTATION USED IN THE BILEVEL E-CVRP MODEL .
General
Upper Level
Symbol
Description
G = (V, A) V = {d} ∪ Vc ∪ Vf A = {(i, j) | i ̸= j, i, j ∈ V } δi dij h Qc Qb σi M
Complete graph representing the E-CVRP instance. Node set consisting of the depot d, customer set Vc , and charging station set Vf . Set of arcs representing all possible connections between distinct nodes. Demand of customer i ∈ Vc . (δi > 0) Distance between nodes i and j, for all (i, j) ∈ A. (dij ≥ 0) Battery consumption rate per unit distance. (h > 0) Maximum cargo capacity of the EV. (Qc > 0) Maximum battery capacity of the EV. (Qb > 0) State-of-charge (SoC) of the EV upon departure from node i ∈ V . (σi ≥ 0) The available number of identical EVs. (M > 0)
Rv = {nv0 , . . . , nvLv +1 }
Ordered sequence of nodes visited by vehicle v, starting and ending at the depot d (nv0 = nvLv +1 = d), where Lv denotes the number of customers served. Complete route of vehicle v obtained by integrating charging decisions Sv into the customer route Rv . Upper-level decision variable representing a set of M customer-serving routes, where some routes Rv may be empty, indicating that vehicle v is not used. Feasible set of all upper-level routing solutions x. Total objective value of a solution with routing plan x and its optimal charging y ∗ (x).
R̄v = Rv ⊕ Sv x = {R1 , . . . , RM } X F (x, y ∗ (x)) Sv = {sv0 , . . . , svLv }
Lower Level
y(x) = {S1 , . . . , SM } Y(x) f (x, y)
The charging decision along route Rv , where each svℓ represents all possible charging station visit configurations between the paired nodes (nvℓ , nvℓ+1 ). v • If at most one charging station is allowed: sℓ ∈ Vf ∪ {NIL}. v • If up to two charging stations are allowed: sℓ ∈ Vf ∪ {(u,w) | u,w ∈ Vf ,u ̸= w} ∪ {NIL}, where (u, w) denotes an ordered pair of distinct charging stations inserted between two nodes in a route. Here, NIL indicates that no charging station is inserted. Lower-level decision variable representing the charging decisions for all M routes in the routing plan x. Feasible set of all lower-level charging solutions y given a routing plan x. Extra travel cost caused by charging detours for routing plan x under decision y.
II. R ELATED W ORK A. Solution Approaches for EVRPs Depending on the scale and realism of modern EVRPs, solution approaches broadly fall into two families: exact and heuristic. Each offers complementary strengths and trade-offs. Exact methods are highly effective for small-scale instances but scale poorly due to exponential complexity. Most rely on MILP formulations and commercial solvers or specialized branch-and-bound frameworks. Zuo et al. [15] employed secant-based piecewise linearization to approximate the concave nonlinear charging, solving instances with up to 25 customers and 4 stations using CPLEX. Kancharla and Ramadurai [16] extended this with nonlinear charging and loaddependent discharging, solving instances with up to 30 customers and 4 stations with Gurobi. Tahami et al. [17] applied branch-and-cut to E-CVRP instances with up to 30 customers and 21 charging stations. More recently, Lam et al. [18] proposed a branch-and-cut-and-price for an Electric Vehicle Routing Problem with Time Windows (EVRPTW) that incorporates piecewise-linear recharging and capacitated charging stations, tested on instances with up to 100 customers. Nafstad et al. [19] developed a branch-price-and-cut for EVRPTW with heterogeneous charging technologies and nonlinear functions, optimally solving instances with up to 100 customers and 21 stations within one hour. These studies adopt highly detailed modeling; when combined with exact algorithms, such models further compound the computational burden, making
it particularly difficult to scale to larger, real-world instances. In contrast to exact algorithms, heuristic methods dominate for medium- and large-scale EVRPs, where exact solvers become impractical. Their efficiency and robustness lead to high-quality approximate solutions within practical time, a key consideration for many real-world applications. A wide range of heuristic and metaheuristic frameworks have been developed for the EVRP and its variants, spanning both population-based and single-point algorithms. In the population-based category, Caillard and Ben Chabane [20] introduced a hybrid approach that combines genetic mechanisms with ACO for the dynamic EVRP with time windows. Yang et al. [21] developed a partheno-genetic algorithm to minimize total distribution costs under time-of-use electricity pricing, considering both fast and regular charging options. Zhen et al. [22] proposed a Particle Swarm Optimization algorithm for a variant of EVRP with different driving modes. In addition, several single-point metaheuristics have been proposed to EVRPs. Felipe et al. [23] used SA for a Green VRP variant that permits partial recharging and accommodates heterogeneous charging technologies. Montoya et al. [11] developed an Iterated Local Search algorithm to optimize nonlinear charging decisions along predetermined vehicle routes. Goeke [24] proposed a granular Tabu Search for the EV Pickup and Delivery Problem with Time Windows. Furthermore, a number of hybrid approaches have also been developed for EVRPs. Seyfi et al.. [25] studied the MultiMode Hybrid EVRP using a matheuristic that combines Vari-
JOURNAL TITLE
4
able Neighborhood Search and mathematical programming. Keskin and Çatay [26] addressed an EVRP with time windows, multiple charging speeds and partial recharges, solving it with an Adaptive Large Neighborhood Search (ALNS) and a MILP metaheuristic. Nolz et al. [27] introduced the consistent EVRP with backhauls and charging-slot management, and developed a hybrid framework that integrates ALNS with constraint programming for charging slots and quadratic programming for trip scheduling. Overall, existing approaches rely either on computationally intensive formulations or heavily customized heuristics, limiting scalability and transparency. B. Bilevel Optimization in Vehicle Routing Problems Many real-world routing problems exhibit hierarchical decision structures in which decisions at one level depend on and influence those at another. While such problems can be formulated as single-level models, doing so often results in inefficiencies and reduced interpretability. Bilevel optimization offers an alternative modeling paradigm for hierarchical routing problems. It reflects a divide-and-conquer idea, decomposing a complex problem into interdependent subproblems that are solved in a coordinated manner. Furthermore, when the resulting subproblems have been well studied in the literature, existing solution techniques can be directly leveraged rather than developing entirely new algorithms from scratch. This allows methodological efforts to focus more on how to decompose the problem itself and how to exploit the interaction and coordination between the two levels. Nearly two decades ago, Marinakis et al. [28] proposed one of the earliest bilevel formulations for the CVRP, combining a Set Partitioning Problem at the upper level and a Traveling Salesman Problem at the lower level, solved by a Genetic Algorithm (GA). Tu et al. [29] later proposed a bilevel Simulated Annealing (SA) for the large-scale MultiDepot Vehicle Routing Problem, where the upper level assigns customers to depots via Voronoi cells (“cluster-first, routesecond” [30]) and the lower level refines routes through reassignment and local search. Zhou et al. [31] developed a bilevel memetic algorithm for the Soft-Clustered VRP, with cluster assignments at the upper level and customer routing at the lower level, integrating group matching-based crossover, bilevel neighborhood search, and tabu-based reconstruction. Jia et al. [6], [7] introduced a bilevel ant colony optimization (ACO) approach for the E-CVRP, combining an “order-first, split-second” method [32] with local search at the upper level and feasible recharging schedules at the lower level. Qin and Chen [33] further proposed a memetic algorithm with adaptive selection, achieving state-of-the-art performance with 6 new best-known solutions on the 17-instance IEEE WCCI-2020 benchmark [1]. Although these studies decompose the ECVRP into routing and charging subproblems, the interaction between the two levels has not been systematically analyzed. III. M ATHEMATICAL M ODELING The E-CVRP is defined on a directed complete graph G = (V, A) with symmetric distances dij = dji , where V = {d} ∪ Vc ∪ Vf and A = {(i, j) | i ̸= j, i, j ∈ V }. Here, d denotes
the depot, Vc the set of customers, and Vf the set of charging stations. Each customer i ∈ Vc has demand δi > 0. Traveling from i to j consumes h dij units of energy, where h is the per-distance energy rate. Each EV is characterized by cargo capacity Qc and battery capacity Qb . We let σi denote the battery level upon departure from node i and M the number of available identical EVs. A. Bilevel Modeling The E-CVRP is formulated using a bilevel model involving two hierarchical decision-making levels. The upper-level decision maker, referred to as the leader, controls customerserving routes. The lower-level decision maker, the follower, determines feasible charging strategies following the routes specified by the leader. The complete bilevel optimization model is formulated as follows, with all relevant notations defined in Table I. Upper-level (routing decision): min F (x, y ∗ (x)),
x∈X
where F (x, y ∗ (x)) =
M X
cost(R̄v ) (1)
v=1
Subject to: M [
(Rv \ {d}) = Vc
(2)
v=1
(Ru \ {d}) ∩ (Rv \ {d}) = ∅, X
δi ≤ Qc ,
∀ u, v ∈ {1, . . . , M }, u ̸= v (3)
∀ v = 1, . . . , M
(4)
i∈Rv
Lower-level (charging decisions for a given x): y ∗ (x) = arg min f (x, y), y∈Y(x)
where f (x, y) =
M X v=1
cost(Sv ) (5)
Subject to: σi = Qb ,
∀i ∈ Vf ∪ {d}
σj = σi − h · dij , 0 ≤ σi ≤ Qb ,
∀(i, j) ∈ R̄v
∀i ∈ R̄v
(6) (7) (8)
The upper-level objective (1) seeks a routing plan x that minimizes the distance traveled by all vehicles, under the assumption that the follower always computes the best feasible charging strategy y ∗ (x) in response. Constraints (2)-(4) ensure that all customers are exclusively assigned to EVs and that the demand on each route does not exceed the vehicle’s capacity. Given a routing plan x, the follower aims to determine the charging strategy y ∗ (x) that minimizes the additional distance traveled when making a detour to charging stations (5). Constraints (6)-(8) ensure battery feasibility: EV departs from the depot or charging station with a full battery, maintains sufficient charge between consecutive nodes, and keeps the battery level within valid bounds throughout the route. The resulting optimal response y ∗ (x) is passed back to the leader, who updates the upper-level objective accordingly.
JOURNAL TITLE
5
s01
0
insertion points (between successive nodes), and assuming that up to two stations may be inserted in each gap, the number of insertion options per insertion point is:
s11
4
s02
0
s12
1 + |Vf | + |Vf |(|Vf | − 1)
s22
Each term of the above corresponds to no insertion, one station, or a pair of distinct stations. Therefore, the total number of feasible insertions for route Rv is:
s32
L +1
0
3
Depot
1
Customer
2
0
Charging Insertion
Fig. 2. An illustrative example of a solution representation for the E-CVRP.
|Yv | = (1 + |Vf | + |Vf |(|Vf | − 1)) v and the full lower-level space is: |Y(x)| =
M Y
|Yv |
v=1
This interaction continues iteratively until termination, yielding an approximate solution pair (x∗ , y ∗ ). Notably, to cover all feasible configurations in the follower’s charging decisions, it is necessary to bound the number of charging stations inserted between any two successive nodes. Thus, we introduce the following assumption: For any two successive nodes in a route, the EV requires no more than two charging station visits to traverse the segment. This assumption aligns with practical applications, as EVs are generally used for short- to medium-range travel, while long-distance routes that may require more than two charging stops are typically served by fuel-based vehicles.
The total solution space of the bilevel E-CVRP is combinatorially large. It grows super-exponentially with the number of customers (due to partition–permutation coupling at the upper level), and for a fixed route, exponentially with the number of insertion points. Moreover, for each insertion point, the number of possible charging station options scales quadratically with the number of charging stations. This makes exact methods intractable even for moderately sized instances; hence, heuristic methods are not only preferable but arguably the only practical option for real-world problems.
B. Solution Space Analysis
In most bilevel programming problems, the leader’s objective value can only be evaluated after solving the follower’s problem. However, in the bilevel model of the E-CVRP, each upper-level decision x naturally defines a preliminary objective ϕ(x), representing the routing cost ignoring battery-related constraints, which is the objective of the classical CVRP subproblem. This structure allows partial evaluation of solution quality without solving the lower-level problem, providing valuable guidance for designing efficient search strategies. When a routing plan x is given, the follower responds with the optimal charging strategy y ∗ (x), resulting in the complete solution (x, y ∗ (x)). Intuitively, ϕ(x) should be generally positively correlated with the complete objective F (x, y ∗ (x)), since y ∗ (x) aims to minimize additional detour distances imposed based on the given routing decision x. This correlation allows ϕ(x) to serve as a cheap surrogate for F (x, y ∗ (x)), especially valuable at early stages of the search process, as evaluating complete solutions is computationally expensive and not necessary. Therefore, when designing the algorithm, more search effort can be allocated to regions with lower ϕ(x) values to accelerate the search. For example, during early-stage convergence, only the routing solution x is evaluated via ϕ(x); full evaluations of (x, y ∗ (x)) are deferred until ϕ(x) drops below a predefined threshold. However, it is worth noting that the correlation between ϕ(x) and F (x, y ∗ (x)) is not strictly monotonic, indicating a potential objective misalignment phenomenon that needs both levels to be jointly optimized in a sophisticated way.
Fig. 2 illustrates a sample solution for the E-CVRP instance with depot 0, customers 1-4, and charging stations 5 and 6. The upper-level solution x consists of two routes: R1 : 0 → 4 → 0 and R2 : 0 → 3 → 1 → 2 → 0. In the corresponding lowerlevel decision y(x), each route Rv is associated with a charging decision sequence Sv , which together form the complete route R̄v . Specifically, R̄1 is given by 0 → s10 → 4 → s11 → 0, and R̄2 is given by 0 → s20 → 3 → s21 → 1 → s22 → 2 → s23 → 0. Each svℓ denotes a charging decision after leaving node ℓ on route Rv , where svℓ can be selected from the set {NIL, 5, 6, (5, 6), (6, 5)}. Upper-Level Space X : The upper level determines a set of routes x = {R1 , . . . , RM }, where each route Rv is an ordered sequence of customer nodes. This decision involves both the partitioning of the customer set and the permutation of nodes within each route. The number of possible partitions increases super-exponentially (Bell number), and the number of permutations per route grows factorially. Therefore, the size of X can be approximated as: |X | ≈ Partition(Vc ) ×
M Y
Lv !
v=1
where Lv is the number of customers served by vehicle v. Lower-Level Space Y(x): Given a routing plan x, the lower level determines where to insert charging stations along each route to ensure battery feasibility. For a route Rv with Lv + 1
C. Objective Space Analysis
JOURNAL TITLE
IV. B ILEVEL L ATE ACCEPTANCE H ILL C LIMBING A. Algorithm Overview The proposed b-LAHC algorithm operates through three key phases. In the initialization phase, a greedy descent procedure (G) first drives the upper-level routing solution x to a local optimum. The exploration phase then performs a systematic neighborhood search (N ) guided by a late acceptance mechanism, where the efficient follower optimizer (FSE ) is conditionally activated to update charging configurations. Throughout this phase, the algorithm maintains a fixed-length history list (P) and monitors search progress through the move acceptance ratio (ρns ). Finally, in the refinement phase, upon termination, the exhaustive enumeration follower optimizer (F) is invoked to polish the best-found solution. This threephase approach effectively balances computational efficiency during the search with solution quality guarantees in the final output, while inheriting the theoretical convergence properties of the original LAHC framework that have been adapted to the bilevel problem structure. Algorithm 1 details b-LAHC procedure. An initial solution x is generated via random customer permutation and the split method [32] (Line 2), followed by greedy descent G to reach a local optimum (Line 3). The efficient follower optimizer FSE is applied to compute the optimal lower-level decision ∗ ySE (x) using the simple enumeration (SE) method (Line 4; see Section IV-E), and the best-known solution (x∗ , y ∗ ) and the corresponding objective F ∗ are updated if improved (Lines 5– 7). Next, the algorithm initializes the best-known surrogate cost ϕ∗ , the iteration counter I, the idle iteration counter Iidle , the number of successful moves n (i.e., the moves within the current history list that accept the candidate solution), and the corresponding success ratio ρn (Line 8). From among these, ρn serves as an additional indicator of search progress, added to the original LAHC scheme. It is updated periodically and enhances convergence detection. Subsequently, a history list P of length Lh is created to store past surrogate costs for the late acceptance criterion. Each entry is initialized by applying noise to ϕ∗ , multiplying it by a random factor between αlb (slightly below 1.0) and αub (slightly above 1.0). This balanced perturbation avoids overly strict or lenient acceptance thresholds during the search (Lines 9–11). At each iteration, N (x, ϕvi , ηmax ) is invoked, yielding a candidate solution x′ and a flag isMoved that records whether the move is accepted (Line 13). If ϕ(x′ ) improves upon ϕ(x), the idle counter Iidle is reset to 0 and ϕ∗ is updated; otherwise, Iidle is incremented by 1 (Lines 14–18). The virtual index vi in the history list is then updated (Line 19). Whenever the history list P completes a cycle, ρn is recomputed and n is reset (Lines 20–22). If the flag isMoved is true, the update process is triggered (Lines 23–31): x is replaced by x′ , n is incremented, and the historical value ϕvi is updated (Line 24). To avoid unnecessary lower-level evaluations and to concentrate the search efforts on promising regions, as discussed in Section III-C, the follower is invoked only when ϕ(x) is close enough to the best surrogate cost ϕ∗ found so far (Line 25). Once the follower is triggered, the optimizer ∗ FSE computes the optimal lower-level decision ySE (x), and the
6
Algorithm 1 Bilevel Late Acceptance Hill Climbing Input: Greedy descent G, Neighborhood exploration N , Follower optimizers FSE and F , History list length Lh , Noise bounds αlb , αub , Max attempts ηmax , Follower activation threshold γ Output: Best bilevel solution (x∗ , y ∗ ) found 1: F ∗ ← +∞, x∗ ← ∅, y ∗ ← ∅ 2: Produce an initial upper-level solution x 3: Run greedy descent until local optimum x ← G(x) ∗ ← FSE (x) 4: ySE ∗ 5: if F (x, ySE ) < F ∗ then ∗ ∗ ∗ ∗ 6: (x , y ) ← (x, ySE ), F ∗ ← F (x, ySE ) 7: end if 8: Set ϕ∗ ← ϕ(x), I ← 0, Iidle ← 0, n ← Lh , ρn ← 1.0 9: for vi ← 0 to Lh − 1 do 10: P[vi] ← ϕ∗ · U NIFORM R ANDOM(αlb , αub ) 11: end for ▷ ϕvi corresponds to P[vi] 12: repeat 13: x′ , isMoved ← N (x, ϕvi , ηmax ) 14: if ϕ(x′ ) < ϕ(x) then 15: Iidle ← 0, ϕ∗ ← min(ϕ(x′ ), ϕ∗ ) 16: else 17: Iidle ← Iidle + 1 18: end if 19: vi ← I mod Lh ▷ Virtual index in history list 20: if vi = 0 then 21: ρn ← n/Lh , n ← 0 22: end if 23: if isMoved then 24: x ← x′ , n ← n + 1, ϕvi ← min(ϕ(x), ϕvi ) 25: if ϕ(x) < γ · ϕ∗ then ∗ 26: ySE ← FSE (x) ∗ 27: if F (x, ySE ) < F ∗ then ∗ ∗ ∗ ∗ 28: (x , y ) ← (x, ySE ), F ∗ ← F (x, ySE ) 29: end if 30: end if 31: end if 32: I ←I +1 33: hasConverged ← (I ≥ 105 and Iidle ≥ 0.02 · I) or ρn ≤ 0.001 34: budgetExceeded ← B UDGET C HECK() 35: until hasConverged or budgetExceeded 36: if hasConverged and ¬budgetExceeded then 37: Restart the algorithm from Line 2 38: end if 39: y ∗ ← F(x∗ ), F ∗ ← F (x∗ , y ∗ ) ▷ Final refinement ∗ ∗ 40: return (x , y )
best-known solution (x∗ , y ∗ ) together with its objective F ∗ are updated whenever an improvement is found (Lines 26–29). At the end of each iteration, I is incremented (Line 32), and two termination flags are updated: • hasConverged: set to true when either (i) excessive idle iterations are detected (I ≥ 105 and Iidle ≥ 0.02I) or (ii) the search success ratio deteriorates (ρn ≤ 0.001). • budgetExceeded: set to true when the computational budget is exhausted.
JOURNAL TITLE
7
M1
M3
M2
a
a
b
a
or
a a
b
b
M7
M6
a
α
α
b
b
b
b
M5
a
M4
β
β
a
α
M8
a
α
b
β
b
β
a
α
a
α
b
β
b
β
a d
d
d
a
d
Fig. 3. The eight move operators illustrated by simplified route segments.
The current loop terminates when either condition is met (Line 35). If convergence is detected before the budget is exhausted, the algorithm restarts from Line 2 and continues until the next termination check (Lines 36–38). Finally, before returning the solution, a final refinement step is performed: the exhaustive follower F is invoked to compute the optimal charging configuration y ∗ (x∗ ) for the best-found routing plan x∗ (Line 39), ensuring that the returned bilevel solution (x∗ , y ∗ ) represents a verified optimum (Line 40). B. Move Operators Move operators define how a current solution x is transformed into a neighboring solution x′ by modifying part of x. Both greedy descent (G) and neighborhood exploration (N ) are built upon these operators. Inspired by Prins [34], eight move operators are adopted to modify the upper-level routing solution x. Let a and b be two distinct customer nodes, which may belong to the same route or different routes. Let α and β denote the immediate successors of a and b, respectively. The eight operators are illustrated in Fig. 3 and summarized below: M1. Relocate a within the same route by inserting it before or after b. M2. Relocate a from its current route and insert it after b in another route. M3. Swap a and b in the same route. M4. Swap a and b between two different routes. M5. In the same route, replace arcs (a, α) and (b, β) with (a, b) and (α, β). M6. In two different routes, replace arcs (a, α) and (b, β) with (a, b) and (α, β). M7. In two different routes, replace arcs (a, α) and (b, β) with (a, β) and (b, α). M8. Remove a from its current route and insert it into an empty route. To apply an operator, a target T is first selected from the current solution x, where T can be either a single route or a pair of routes. From T , two distinct customer nodes (a, b) are chosen. The selected operator op is subsequently applied,
yielding a neighboring solution x′ . This process is formally expressed as x′ = op(x, T, a, b). These operators can be categorized into three groups: (i) intra-route moves (M1, M3, M5), (ii) inter-route moves (M2, M4, M6, M7), and (iii) inter-route moves, involving empty routes (M8). M8 differs from the others in that it moves a node to an empty route. Only M8 can increase the total number of routes in the solution x, whereas M1–M7 either reduce or preserve the current route count. This operator is designed to handle special cases where using more vehicles may result in a lower total cost. M8 facilitates the exploration of more diverse regions in the solution space and helps the search to escape from local optima. C. Greedy Descent The greedy descent procedure G(x) performs iterative improvement on x until reaching a local optimum xopt , i.e., G(x): x 7→ xopt . This computationally efficient exploitation rapidly reduces ϕ(x), directing the search toward promising regions of the search space. Algorithm 2 performs a multi-level best-improvement descent over an operator set I (excluding M8). Each outer iteration resets the any-improvement flag aImp and shuffles I to avoid ordering bias (Line 2). For each operator op ∈ I, the operator-level flag oImp is cleared (Line 4). Then, for every target T (a route or route pair) in the current solution x, the algorithm repeatedly scans pairs (a, b) ∈ T and applies op to generate x′ (Lines 5–18). If ϕ(x′ ) < ϕ(x), it updates x and sets both mDone (move-applied flag) and oImp; otherwise it continues (Lines 8–15). At the end of each attempt, information about the affected route or route pair is updated to reflect the latest structural changes in x (Line 16). This step is specifically for inter-route moves as they may create empty routes and thus change the subsequent targets in x (affecting Line 5). The per-target loop stops when no further improving move is found. After all targets for the current operator are processed, aImp is updated (Line 19). The outer loop repeats while any improvement exists. Once no operator yields further
JOURNAL TITLE
Algorithm 2 Greedy Descent G Input: Current solution x, set of move operators I Output: A locally optimal solution xopt 1: repeat 2: aImp ← False, S HUFFLE(I) 3: for each operator op ∈ I do 4: oImp ← False 5: for each target T ∈ x do ▷ Enumerate 6: repeat 7: mDone ← False 8: for each (a, b) ∈ T do ▷ Enumerate 9: x′ ← op(x, T, a, b) 10: if ϕ(x′ ) < ϕ(x) then 11: x ← x′ 12: mDone ← True, oImp ← True 13: break 14: end if 15: end for 16: U PDATE TARGET I NFO(x) 17: until mDone = False 18: end for 19: aImp ← aImp ∨ oImp 20: end for 21: until aImp = False 22: return xopt ← x
improvement, the procedure terminates and returns the locally optimal solution xopt . D. Neighborhood Exploration The neighborhood exploration component drives the upperlevel solution x through iterative transitions across its solution space, guided by the late acceptance mechanism. Algorithm 3 details this procedure. The process begins by initializing the candidate solution x′ and the move flag (Line 1). The attempt counter i is initialized (Line 2), and a move operator op is selected uniformly at random (u.a.r.) from the operator set M (Line 3). M includes all operators described in Section IV-B. At each attempt, a target T is u.a.r. chosen from the current solution x (Line 5). From T , a customer node a is selected (Line 6), and the candidate positions b are enumerated based on the chosen operator op (Line 7). For operator M8, b refers to the depot node of an empty route. For each feasible position b (Line 8), a candidate neighbour solution x′′ is generated by applying the operator (Line 9). The late acceptance condition then checks if x′′ is better than either ϕ(x) or the historical value ϕvi (Line 10). Once satisfied, x′′ is accepted as the new neighbour x′ , isMoved is set to true (Line 11). The process terminates early once the condition is satisfied or after reaching the maximum number of attempts ηmax (Line 16). E. Lower-level Decision The lower-level charging decision involves optimally inserting charging stations along a fixed route. This problem
8
Algorithm 3 Neighborhood Exploration N Input: Current solution x, history solution cost ϕvi , maximum attempts ηmax , move operators M Output: A neighbour solution x′ , a flag isMoved 1: x′ ← x, isMoved ← false 2: i ← 0 ▷ Attempt counter 3: Randomly select a move operator op ∈ M 4: repeat 5: Randomly select a target T ∈ x 6: Randomly select a node a ∈ T 7: Enumerate candidate positions b based on op 8: for each feasible b do 9: Generate x′′ ← op(x, T, a, b) 10: if ϕ(x′′ ) < ϕvi or ϕ(x′′ ) < ϕ(x) then 11: x′ ← x′′ , isMoved ← true 12: break 13: end if 14: end for 15: i←i+1 16: until isMoved or i ≥ ηmax 17: return x′ , isMoved
is referred to as the Fixed Route Vehicle Charging Problem (FRVCP), which has been proven to be NP-hard [11]. The lower-level decision space Y(x) is fully determined by the upper-level routing solution x, as introduced in Section III-B. Within Y(x), all potential charging insertion points, i.e., the gaps between consecutive nodes, along each route are considered simultaneously. Moreover, multiple possible charging configurations may exist for each insertion point. To reduce the complexity of the lower-level decision, we bound the number of charging station visits on each route Rv between ℓbv and ℓbv + 1, where ℓbv denotes the minimum required visits for Rv . It is calculated by dividing the total travel distance of the route ϕ(Rv ) by the maximum cruising range of an EV at full charge Qb /h: ϕ(Rv ) ℓbv = Qb /h Under this restriction, we exhaustively enumerate all feasible configurations with up to two consecutive charging station visits at each insertion point. The resulting optimal charging configuration is denoted as y ∗ , and the exhaustive enumeration procedure is referred to as F. However, since F is computationally expensive, it would be impractical to include it throughout the search. Thus, an efficient lower-level optimization alternative is essential. Specifically, we adopt the Simple Enumeration (SE) method proposed in [7], which significantly reduces computational complexity by restricting the search to configurations with at most one charging station per insertion point and precomputing the “best charging station” for every pair of customers (or depot-customer pairs). This best charging station is selected to minimize the additional travel distance incurred by the detour, as defined below: θij = arg min (diθ + dθj ) , θ∈Vf
∀i, j ∈ Vc ∪ {d}, i ̸= j
JOURNAL TITLE
9
TABLE II D ETAILS OF THE IEEE WCCI-2020 BENCHMARK SET
TABLE IV M AIN PARAMETERS OF THE B -LAHC ALGORITHM
name
|Vc |
|Vf |
M
Qc
Qb
h
UB
Parameter
Description
Value
E22 E23 E30 E33 E51 E76 E101
21 22 29 32 50 75 100
8 9 6 6 5 7 9
4 3 4 4 5 7 8
6000 4500 4500 8000 160 220 200
94 190 178 209 105 98 103
1.2 1.2 1.2 1.2 1.2 1.2 1.2
384.67 573.13 511.25 869.89 570.17 723.36 899.88
Lh ηmax γ
History length Maximum attempts Follower activation threshold
5723 60 1.01
X143 X214 X351 X459 X573 X685 X749 X819 X916 X1001
142 213 350 458 572 684 748 818 915 1000
4 9 35 20 6 25 30 25 32 9
7 11 40 26 30 75 98 171 207 43
1190 944 436 1106 210 408 396 358 33 131
2243 987 649 929 1691 911 790 926 1591 1684
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
– – – – – – – – – –
TABLE III OVERVIEW OF M ETHODS FOR THE WCCI-2020 EVRP B ENCHMARK Methods
Stop Criteria Language/Tool Code Available
MILP [1] VNS [35] SA GA [36] HHASA-TS [8] BACO [6] CBACO-I [7] TAMLS [9] CBMA [33] b-LAHC
Max Evals Max Evals Max Evals Max Evals Max Evals Max Time Max Time Max Time Max Time Both
Gurobi C++ C++ C++ MATLAB C++ C++ C++ C++ C++
No Yes No Yes Yes Yes Yes No Yes Yes
Under SE, the number of configurations per route is determined by selecting either ℓbv or ℓbv + 1 insertion points from the Lv + 1 available gaps, with the optimal charging station at each chosen point fixed in advance: Lv + 1 Lv + 1 + ℓbv ℓbv + 1 where Lv + 1 is the number of insertion points. The total complexity is: M Y Lv + 1 Lv + 1 |YSE (x)| = + ℓbv ℓbv + 1 v=1 which is polynomial in practice, since ℓbv is typically small (in most cases, an EV requires no more than three charging stops per trip), and Lv grows linearly with the size of the problem. SE thus reduces the theoretically exponential complexity of Y(x) to a tractable scale even for large-scale instances. We denote the optimal charging configuration returned by SE as ∗ ySE , and the SE process as FSE . In Algorithm 1, FSE serves as the default efficient optimizer, whereas F is reserved for final solution refinement to ensure optimality, thereby achieving a balance between computational efficiency and solution quality. V. E XPERIMENTS A. Benchmark Suite and Reference Methods We evaluate b-LAHC on the IEEE WCCI-2020 EVRP benchmark [1], covering 17 instances (7 small E-set, 10 large X-set) under two commonly used termination criteria (Max Evals and Max Time). Table II summarizes the key
TABLE V R ANK CORRELATION AND T OP -k% SET OVERLAP BETWEEN ϕ(x) AND F (x, y ∗ (x)). Instance #Samples E22 E23 E30 E33 E51 E76 E101 X143 X214 X351 X459 X573 X685 X749 X819 X916 X1001
τb
60,163 0.9245 100,717 0.9646 182,830 0.9750 115,460 0.9529 391,533 0.9802 712,352 0.9755 854,090 0.9750 1,354,989 0.9812 3,752,720 0.9785 12,244,328 0.9877 8,990,056 0.9903 12,317,882 0.9811 31,464,639 0.9869 46,890,428 0.9867 61,717,308 0.9661 78,677,293 0.9666 27,402,846 0.9889
Recall@1% Recall@5% Recall@10% Recall@20% 0.7425 0.8978 0.9043 0.8823 0.7439 0.7650 0.7685 0.9728 0.5141 0.8522 0.9291 0.0278 0.4730 0.5526 0.3169 0.4708 0.6183
0.8136 0.9255 0.9354 0.9531 0.9279 0.8990 0.9201 0.9341 0.8491 0.9558 0.9639 0.8562 0.9087 0.9034 0.6893 0.6854 0.9256
0.8764 0.9387 0.9567 0.9671 0.9631 0.9157 0.9478 0.9542 0.9423 0.9702 0.9902 0.9429 0.9630 0.9653 0.8790 0.8890 0.9671
0.9253 0.9646 0.9745 0.9792 0.9814 0.9676 0.9743 0.9792 0.9857 0.9876 0.9924 0.9846 0.9879 0.9848 0.9616 0.9620 0.9904
characteristics of all instances. These instances are divided into two categories: • Each instance in E-set contains no more than 100 customers. For these instances, upper bounds (UB) are provided as baseline results, some of which are optimal. • Each instance in X-set contains between 142 and 1,000 customers, and no known upper bounds are provided. These instances present different characteristics in customers’ spatial distribution and their demands, significantly affecting the structure of feasible solutions. Two different computational budgets are available. The first one is defined as a multiple of the problem size pz, where pz = |{d}|+|Vc |+|Vf |. The maximum number of evaluations is given by Max Evals = 25,000 × pz The evaluation of a solution requires O(pz) time. In the official evaluator, each access to an arc weight dij consumes a fraction 1/pz of the budget, meaning that even neighborhood exploration contributes to the overall evaluation count. The second budget is defined in [6], [7] as: |Vc | + |Vf | (hours) 100 where the parameter ω is set to 1, 2, and 3 for instance groups E22–E101, X143–X916, and X1001, respectively. Max Evals facilitates a fair comparison across algorithms. In contrast, Max Time is subject to differences in hardware, programming languages, compiler settings, and parallelization strategies. These factors make it difficult to quantify the computational effort and may lead to inconsistent performance assessments. Both criteria are adopted in this study. Table III summarizes the state-of-the-art methods for the benchmark. Among them, the VNS, SA, and GA implementations correspond to the top three algorithms submitted to the WCCI-2020 EVRP competition. Max Time = ω ·
JOURNAL TITLE
10
TABLE VI O BJECTIVE VALUE COMPARISON WITH STATE - OF - THE - ART ALGORITHMS USING THE M AX E VALS STOP CRITERION Instance
Index
BKS
VNS
SA
GA
HHASA-TS
BACO
CBACO-I
CBMA
b-LAHC
Gap
E22
best mean std.
384.67 384.67 –
384.67 384.67 2.11
384.67 384.67 0.00
384.67 384.67 0.00
384.67 384.67 0.00
384.67 384.67 0.00
384.67 384.67 0.00
384.67 384.67 0.00
384.67 385.25 0.28
– 0.15% –
E23
best mean std.
571.94 571.94 –
571.94 571.94 0.00
571.94 571.94 0.00
571.94 571.94 0.00
571.94 571.94 0.00
571.94 571.94 0.00
571.94 571.94 0.00
571.94 571.94 0.00
571.94 571.94 0.00
– 0.00% –
E30
best mean std.
509.47 509.47 –
509.47 509.47 0.00
509.47 509.47 0.00
509.47 509.47 0.00
509.47 509.47 0.00
509.47 509.47 0.00
509.47 509.47 0.00
509.47 509.47 0.00
509.47 509.47 0.00
– 0.00% –
E33
best mean std.
840.14 840.43 –
840.14 840.43 1.18
840.57 854.07 12.80
844.25 845.62 0.92
840.14 840.70 1.40
846.02 847.20 1.21
846.96 847.74 0.56
840.57 840.57 0.00
840.14 840.14 0.00
– -0.03% –
E51
best mean std.
529.90 529.90 –
529.90 543.26 3.52
533.66 533.66 0.00
529.90 542.08 8.57
529.90 536.98 7.27
546.85 555.48 4.47
529.90 535.39 4.72
529.90 529.90 0.00
529.90 531.19 3.59
– 0.24% –
E76
best mean std.
692.64 694.61 –
692.64 697.89 3.09
701.03 712.17 5.78
697.27 717.30 9.58
692.74 694.96 1.63
743.20 757.08 5.79
726.02 735.34 4.85
692.64 694.61 2.10
696.27 699.54 2.91
– 0.71% –
E101
best mean std.
837.10 843.10 –
839.29 853.34 4.73
845.84 852.48 3.44
852.69 872.69 9.58
837.10 843.10 3.90
905.09 905.09 0.00
905.09 905.09 0.00
841.80 844.96 2.63
840.55 844.01 3.16
– 0.11% –
X143
best mean std.
15910.86 16214.37 –
16028.05 16459.31 242.59
16610.37 17188.90 170.44
16488.60 16911.50 282.30
15910.86 16214.37 215.77
17407.82 17407.82 0.00
17407.82 17407.82 0.00
16288.88 16489.90 142.22
15916.19 16103.45 85.93
– -0.68% –
X214
best mean std.
11090.28 11206.60 –
11323.56 11482.20 76.14
11404.44 11680.35 116.47
11762.07 12007.06 156.69
11090.28 11206.60 84.58
12267.88 12267.88 0.00
12267.88 12267.88 0.00
11461.70 11608.98 90.60
11114.46 11245.67 87.27
– 0.35% –
X351
best mean std.
26622.42 26750.60 –
27064.88 27217.77 86.20
27222.96 27498.03 155.62
28008.09 28336.07 205.29
26622.42 26750.60 102.55
28639.35 28932.86 132.28
28394.79 28571.17 89.85
27204.23 27447.89 175.23
26370.42 26448.67 58.10
– -1.13% –
X459
best mean std.
24794.35 25041.10 –
25370.80 25582.27 106.89
25464.84 25809.47 157.97
26048.21 26345.12 185.14
24794.35 25041.10 237.58
27371.06 27371.06 0.00
27371.06 27371.06 0.00
25615.60 25942.22 140.26
24465.89 24600.78 79.95
– -1.76% –
X573
best mean std.
51436.90 51776.70 –
52181.51 52548.09 278.85
51929.24 52793.66 577.24
54189.62 55327.62 548.05
51436.90 51776.70 166.86
56884.09 57342.53 240.17
57018.02 57417.98 248.69
52458.56 52748.36 253.97
51234.65 51379.40 107.62
– -0.77% –
X685
best mean std.
69955.95 70401.25 –
71345.40 71770.57 197.08
72549.90 73124.98 320.07
73925.56 74508.03 409.43
69955.95 70401.25 218.98
78045.96 78045.96 0.00
78045.96 78045.96 0.00
73975.35 74412.32 253.08
69415.64 69610.37 123.60
– -1.12% –
X749
best mean std.
79779.87 80135.67 –
81002.01 81327.39 176.19
81392.78 81848.13 275.26
84034.73 84759.79 376.10
79779.87 80135.67 219.50
85781.59 85781.59 0.00
85640.27 85767.46 44.69
81939.07 82657.87 368.48
78809.51 78937.48 90.58
– -1.50% –
X819
best mean std.
161924.79 162530.67 –
164289.95 164926.41 318.62
165069.77 165895.78 403.70
170965.68 172410.12 568.58
161924.79 162530.67 289.41
174375.10 174840.49 230.07
173145.29 173596.85 276.31
165577.34 166464.86 507.42
159864.38 160001.36 175.36
– -1.56% –
X916
best mean std.
336717.71 337641.92 –
341649.91 342460.70 510.66
342796.88 343533.85 556.98
357391.57 360269.94 1192.57
336717.71 337641.92 461.47
362477.43 363118.31 489.37
360138.18 360745.35 423.35
343838.22 345689.82 1147.38
332309.16 332558.20 174.34
– -1.51% –
X1001
best mean std.
75469.29 75931.28 –
77476.36 77920.52 234.73
78053.86 78593.50 306.27
78832.90 79163.34 229.19
75469.29 75931.28 304.10
80523.17 80523.17 0.00
80523.17 80523.17 0.00
78075.73 78897.32 345.21
74893.30 75130.68 135.43
– -1.05% –
w/t/l
vs b-LAHC
2/2/13
1/2/14
1/2/14
4/2/11
1/2/14
1/2/14
3/2/12
3.53
4.59
5.59
2.59
7.00
6.59
3.94
rank
B. Experimental Setup Table IV summarizes the control parameters used in bLAHC. Default settings were obtained via automatic configuration using IRACE [37] on the full benchmark set. All experiments were conducted on the Sulis HPC platform using dedicated AMD EPYC 7742 nodes (2.25 GHz, 64-core processors), restricted to single-thread execution and 1 GB RAM per run. Unless otherwise specified, the budget is based on Max Evals to ensure a fair comparison across algorithms. All available algorithms listed in Table III were evaluated under strictly identical experimental protocols, including hardware, programming language (C++), compiler (GCC 13.3.0), flag(-O3), and 10 independent runs per instance. The only exception was HHASA-TS, which was executed using its
2.18
native MATLAB implementation due to unavailable sourcelevel porting. For SA, results are taken directly from the official competition report, while TAMLS results are reproduced directly from its original paper. For statistical significance testing, we employed the nonparametric Friedman test (α = 0.05), followed by a Holm post-hoc correction to account for multiple comparisons.
C. Correlation vs. Misalignment: Validating ϕ(x) as a Surrogate Objective To empirically validate that the surrogate routing cost ϕ(x) is positively correlated with the complete objective F (x, y ∗ (x)), we collected unique (ϕ(x), F (x, y ∗ (x))) pairs by
JOURNAL TITLE
11
VNS
X143
GA
BACO
X214
2.40 2.25
4.2
1.8
3.9
CBMA
b-LAHC
X459
X573
4.0
7.0
3.6
3.6
1.6
1.95
HHASA-TS
X351
2.0
2.10
Objective value (×104)
CBACO-I
3.3
3.2
3.0
2.8
6.5 6.0
1.4
1.80
1.2
1.65 1e−06
1e−04
1e−02
2.7 1e−06
X685
1e−04
1e−02
1
1e−06
X749
1e−02
1
1e−06
1e−04
1e−02
1
20
9.6
8.0
8.8
7.2
8.0 1e−06
1e−04
1e−02
1
1e−02
1
1e−02
1
X1001
11
37.5
10
36.0
18
1e−04
12
19 8.8
1e−06
X916 39.0
21
10.4
9.6
1e−04 X819
11.2
10.4
5.5
2.4
9 34.5
17
8 16 1e−06
1e−04
1e−02
1
1e−06
1e−04
1e−02
1
33.0
1e−06
1e−04
1e−02
1
1e−06
1e−04
Evaluations (×107)
Fig. 4. Convergence trends on the large-scale instances.
recording solutions encountered during the search process of b-LHAC. Table V summarizes the results in terms of Kendall’s τb rank correlation and Top-k% set recall between ϕ(x) and F (x, y ∗ (x)). Across all 17 instances, the Kendall’s τb values are on average above 0.97, with τb > 0.95 in 15 out of 17 instances. This confirms a very strong positive correlation between ϕ(x) and F (x, y ∗ (x)). However, the Top-k% set recall reveals a more nuanced picture. While Recall@10% and Recall@20% remain above 0.95 in most instances, Recall@1% exhibits noticeable drops in several large-scale X-instances (e.g., 0.0278 for X573 and 0.3169 for X819), even when τb is high (0.9811 and 0.9661, respectively). This discrepancy indicates a phenomenon of objective misalignment. Although ϕ(x) and F (x, y ∗ (x)) are largely well-aligned, they may disagree on the precise ranking of top-performing solutions. In other words, solutions that appear optimal under ϕ(x) are not necessarily optimal under F (x, y ∗ (x)). Nevertheless, the optimal or near-optimal complete solutions under F (x, y ∗ (x)) are highly likely to be found within the broader top-5% or top-10% regions defined by ϕ(x). In summary, using ϕ(x) as a cheap yet reliable surrogate objective during the search, significantly reduces lower-level evaluations while preserving solution quality. However, the observed objective misalignment highlights the necessity of a unified bilevel optimization framework, as a naive twostage approach may yield solutions far from the true optimal complete solution.
D. Comparison with State-of-the-art Algorithms 1) Solution Quality under Max Evals: Table VI reports a head-to-head comparison of b-LAHC against state-of-the-art algorithms. “BKS” denotes the best-known solution. “Gap” is computed as the percentage difference between the b-LAHC mean and the BKS (negative means b-LAHC outperforms BKS). “w/t/l” counts the wins, ties, and losses of a compared algorithm against b-LAHC across all instances. “rank” shows
the average rank over all instances (lower is better) based on the Friedman test. Small-scale (E-set, 7 instances). b-LAHC beats BKS on E33, ties BKS on E23 and E30, and is very close on the remaining four, with the largest deviation +0.71% on E76. • Large-scale (X-set, 10 instances). b-LAHC surpasses BKS on 9 out of 10 instances and underperforms only on X214 (+0.35%). These results confirm that b-LAHC scales particularly well with instance size. • Statistical significance and overall ranking. A Friedman test over all 17 instances and 8 algorithms yields χ2 = 74.67 with p = 1.67 × 10−13 < 0.05, rejecting the null hypothesis and indicating statistically significant performance differences. The overall ranking is: b-LAHC ≻ HHASA-TS ≻ VNS ≻ CBMA ≻ SA ≻ GA ≻ CBACO-I ≻ BACO. • Head-to-head against the second-best method (HHASATS). b-LAHC achieves 11 wins / 2 ties / 4 losses at the mean level. On the winning subset, b-LAHC’s average gap to BKS is -0.99%; on the losing subset it is +0.33%. b-LAHC shows lower standard deviation than HHASA-TS on 12 out of 17 instances overall (3 higher, 2 ties), including 9 out of 10 large X-instances. Averaged over non-zero HHASATS cases, b-LAHC reduces standard deviation by 39.5% overall (22.8% on E-instances; 46.2% on X-instances). This indicates that b-LAHC is not only more effective on most large X-instances but also more robust. • Takeaway. Under the Max Evals budget, b-LAHC delivers near-BKS or better results on small E-instances and shows a clear advantage on large X-instances, achieving the top overall Friedman ranking with strong statistical significance. •
2) Convergence Analysis: Figure 4 shows convergence on large-scale X-instances with a log-scaled X-axis, which highlights early convergence. b-LAHC shows no advantage in the early phase: performing worse than CBMA (X143, X214, X459), VNS (X685–X916), and CBMA, CBACO-I, and BACO (X1001). This suggests that these algorithms benefit from stronger initialization, whereas b-LAHC starts from a randomly generated solution and relies on the greedy descent
JOURNAL TITLE
12
Objective Value: 545.01
Instance Index Baseline
Remove G
Remove F
γ=0
Obj
Gap
Obj
Gap
Obj
Gap
best mean std.
384.67 385.25 0.28
384.67 384.96 0.35
– -0.08% –
384.67 385.25 0.30
– 0.00% –
385.39 394.81 13.61
– 2.48% –
E23
best mean std.
571.94 571.94 0.00
571.94 571.94 0.00
– 0.00% –
571.94 571.94 0.00
– 0.00% –
571.94 579.39 18.94
– 1.30% –
E30
best mean std.
509.47 509.47 0.00
509.47 511.26 1.98
– 0.35% –
509.47 509.47 0.00
– 0.00% –
509.47 530.06 20.49
– 4.04% –
E33
best mean std.
840.14 840.14 0.00
840.14 840.47 0.22
– 0.04% –
840.14 840.14 0.00
– 0.00% –
845.57 870.16 14.54
– 3.57% –
E51
best mean std.
529.90 531.19 3.59
535.44 538.47 2.47
– 1.37% –
529.90 531.19 3.79
– 0.00% –
529.90 568.05 19.51
– 6.94% –
E76
best mean std.
696.27 699.54 2.91
723.91 727.54 2.31
– 4.00% –
696.27 699.54 3.07
– 0.00% –
722.96 733.92 11.07
– 4.91% –
E101
best mean std.
840.55 844.01 3.16
880.66 882.26 1.07
– 4.53% –
840.55 844.32 3.67
– 0.04% –
865.80 894.58 16.10
– 5.99% –
X143
best 15916.19 16794.26 mean 16103.45 16922.98 std. 85.93 66.97
– 5.09% –
15916.19 – 17276.88 – 16110.37 0.04% 17794.60 10.50% 81.94 – 321.14 –
X214
best 11114.46 11202.04 mean 11245.67 11307.11 std. 87.27 52.07
– 0.55% –
11136.01 – 11823.46 – 11252.84 0.06% 12104.47 7.64% 88.41 – 133.49 –
X351
best 26370.42 26314.51 mean 26448.67 26448.35 std. 58.10 77.10
– 0.00% –
26383.02 – 28110.86 – 26488.96 0.15% 28310.11 7.04% 56.99 – 134.67 –
X459
best 24465.89 25510.22 mean 24600.78 26566.49 std. 79.95 492.94
– 7.99% –
24486.95 – 26309.21 – 24612.92 0.05% 26697.30 8.52% 92.52 – 192.83 –
X573
best 51234.65 67566.80 – 51308.77 – 52868.26 – mean 51379.40 73086.69 42.25% 51540.36 0.31% 53563.39 4.25% std. 107.62 3506.93 – 149.05 – 395.22 –
X685
best 69415.64 85787.12 – 70429.18 – 73197.36 – mean 69610.37 94291.09 35.46% 70634.31 1.47% 73491.89 5.58% std. 123.60 4974.22 – 150.31 – 193.54 –
X749
best 78809.51 93402.93 – 79139.04 – 82928.79 – mean 78937.48 97470.40 23.48% 79282.97 0.44% 83497.56 5.78% std. 90.58 2183.00 – 88.20 – 308.47 –
X819
best 159864.38 160299.64 – 161950.84 – 164463.38 – mean 160001.36 161290.69 0.81% 162238.74 1.40% 164893.76 3.06% std. 175.36 669.36 – 166.57 – 203.64 –
X916
best 332309.16 339313.98 – 335162.52 – 339887.86 – mean 332558.20 341834.68 2.79% 335479.71 0.88% 340409.34 2.36% std. 174.34 1583.63 – 261.13 – 478.17 –
E22
best 74893.30 326525.74 – 75312.64 – 80837.17 – X1001 mean 75130.68 360120.73 379.33% 75594.82 0.62% 82010.74 9.16% std. 135.43 20775.16 – 160.11 – 497.33 –
(G) phase to reach a local optimum. In the later phase, however, b-LAHC clearly dominates on all instances except for X214. This strongly indicates that the success of b-LAHC lies in the strength of its bilevel framework design, despite its weaker initialization. 3) Solution Quality under Max Time: Under the Max Time setting, a Friedman test conducted across all 17 instances and 5 algorithms yields χ2 = 40.78 with a p-value of 2.99 × 10−8 < 0.05, indicating statistically significant performance differences. The average ranking is b-LAHC ≻ TAMLS ≻ CBMA ≻ CBACO-I ≻ BACO, where b-LAHC achieves the lowest average rank and consistently outperforms competing methods, particularly on large X-instances. Detailed experimental results are provided in the Supplementary Material.
Y Coordinate
TABLE VII A BLATION STUDY ON THE EFFECTS OF G REEDY D ESCENT (G), C HARGING R EFINEMENT (F ), AND B ILEVEL FRAMEWORK
Objective Value: 509.47
410
410
400
400
390
390
380
380
370
370
360
360
350
350
340
340
330
330 100
120
140
160
180
200
220
100
120
140
160
180
200
220
X Coordinate Depot
Customer
Charging Station
Fig. 5. Comparative solution with and without M8 operator on E30.
E. Ablation Study on Key Components Table VII summarizes the results of the ablation study, evaluating the contributions of the Greedy Descent (G) and Charging Refinement (F) components. In addition, we also examine the case γ = 0, where b-LAHC degenerates into a naive two-stage pipeline with upper- and lower-level decisions solved separately. The baseline corresponds to b-LAHC executed with the parameters settings given in Table IV. Removing G causes a substantial degradation in performance on 15 out of 17 instances, with an average mean gap of +29.88% relative to the baseline, and a median gap of +2.79%. The deterioration is especially severe on large-scale instances such as X1001 and X573, while the impact on small E-instances is negligible. Removing F has a moderate effect but still worsens performance on 13 out of 17 instances, with an average mean gap of +0.32% and a median gap of only +0.05%. The influence of F is negligible on small E-instances but becomes critical on large-scale cases where charging refinement is needed. For instance, excluding F results in noticeable degradations on X685, X819, and X916. It is also worth noting that incorporating F incurs only a modest computational overhead. Enabling F adds less than 0.1 s of runtime for all instances. This negligible overhead justifies the inclusion of F given its significant improvement in solution quality. Setting γ = 0 leads to a consistent deterioration across nearly all instances, highlighting the importance of the bilevel framework. The average mean gap increases to +5.48% with a median gap of +5.58%, markedly higher than the baseline. The degradation is already visible on small E-instances, and it becomes even more pronounced on large-scale cases such as X143, X351, and X1001. Compared with removing F, the impact of disabling the bilevel framework is significantly larger, underscoring the necessity of jointly optimizing upper- and lower-level decisions rather than treating them in a decoupled two-stage manner. Overall, these findings indicate that G and F play complementary roles, while the bilevel framework itself is indispensable. Their synergy enables the full b-LAHC algorithm to achieve superior solution quality, scalability, and stability across diverse problem scales. The neighborhood exploration initially comprised M1–M7, all of which preserve or reduce the number of routes. They led to unstable behavior in certain instances. In particular, on E30
JOURNAL TITLE
(Fig. 5), some runs converged to inferior three-route solutions (objective value 545.01), whereas the best-known outcomes are four-route solutions (509.47). This counter-intuitive phenomenon arises from the EV setting, where charging detours and route geometry can favor additional routes. To address this, we introduced M8, the only move able to increase the route count. M8 explores regions of the solution space that are unreachable by M1–M7 and provides an escape from local optima caused by under-utilized fleets. Empirically, enabling M8 stabilizes convergence on E30 (Fig. 5), with runs consistently reaching the superior four-route solution, and we observed analogous improvements in other cases. F. Sensitivity Analysis of Hyperparameters The optimal choice of Lh depends strongly on problem scale, consistent with Lobo et al. [38]. Very small values lead to premature convergence, whereas excessively large values slow progress under a fixed budget. A moderate setting provides the best trade-off, and Lh = 5723 is adopted as it yields stable performance across instances. ηmax controls neighborhood exploration effort. The algorithm shows strong robustness to this parameter; overly small values reduce exploration on large instances, while overly large values waste evaluations on small ones. The default setting ηmax = 60 achieves a balanced performance. γ regulates how frequently the lower-level optimizer is invoked. Experiments indicate stable performance when γ ∈ [1.01, 1.05], which balances upperlevel and lower-level decisions. We therefore use γ = 1.01 as the default setting. Detailed sensitivity analyses for these parameters are provided in the Supplementary Material. VI. C ONCLUSIONS AND F UTURE W ORK This paper introduced a bilevel late acceptance hill climbing algorithm (b-LAHC) for the bilevel formulation of the ECVRP. By leveraging a surrogate objective at the upper level with an efficient charging optimizer at the lower level, bLAHC is lightweight and interpretable. Despite using fixed parameters, it delivers efficient search across instances. Compared with state-of-the-art methods, b-LAHC follows a distinct design philosophy. HHASA-TS [8] integrates SA with reinforcement learning to dynamically adjust the search, improving adaptiveness but incurring computational overhead and higher variability on large instances. CBMA [33] and CBACO-I [7] employ adaptive diversity control in population-based frameworks, enabling effective exploration on small/medium instances but showing slower convergence and reduced robustness on larger benchmarks. Notably, bLAHC attains comparable or superior results as a lightweight single-point method without relying on complex adaptive mechanisms. These findings suggest two directions: incorporating adaptive learning into b-LAHC or incorporating promising features in b-LAHC into a population-based framework. Interestingly, M8 operator plays a more critical role than expected. In some instances (e.g., E30), optimal solutions require more vehicles due to charging detours and battery limits. Without M8, the search converges prematurely to suboptimal solutions featuring fewer vehicles. This highlights
13
the importance of route-creation moves for handling complex EV charging behaviors. The experiments provide strong support for the proposed bilevel framework, which leverages the dominance of upperlevel routing decisions in the E-CVRP, while offering a principled mechanism to determine when to introduce joint optimization, thereby ensuring greater generality and adaptability across diverse problem instances. Future research can proceed in three main directions. First, we plan to extend the current b-LAHC algorithm by incorporating learning mechanisms to enable adaptive behavior during the search. Second, we aim to incorporate useful features in the single-point framework into a population-based metaheuristic capable of multimodal search, allowing a systematic comparison between single-point and population-based strategies. Third, the proposed bilevel optimization model can be applied to other routing problems with hierarchical decision structures, offering a unified perspective for tackling bilevel routing challenges beyond the E-CVRP. ACKNOWLEDGMENTS The authors thank the IT support team at QMUL. Computations were performed on the Sulis Tier-2 HPC platform, funded by EPSRC (EP/T022108/1) and HPC Midlands+. R EFERENCES [1] M. Mavrovouniotis, C. Menelaou, S. Timotheou, C. Panayiotou, G. Ellinas, and M. Polycarpou, “Benchmark set for the ieee wcci-2020 competition on evolutionary computation for the electric vehicle routing problem,” KIOS CoE, University of Cyprus, Tech. Rep., 2020, technical Report. [Online]. Available: https://mavrovouniotis.github.io/EVRPcompetition2020/ [2] P. Toth and D. Vigo, Vehicle routing: problems, methods, and applications. SIAM, 2014. [3] S. Erdoğan and E. Miller-Hooks, “A green vehicle routing problem,” Transportation research part E: logistics and transportation review, vol. 48, no. 1, pp. 100–114, 2012. [4] M. Schiffer, M. Schneider, G. Walther, and G. Laporte, “Vehicle routing and location routing with intermediate stops: A review,” Transportation Science, vol. 53, no. 2, pp. 319–343, 2019. [5] A. Froger, J. E. Mendoza, O. Jabali, and G. Laporte, “Improved formulations and algorithmic components for the electric vehicle routing problem with nonlinear charging functions,” Computers & Operations Research, vol. 104, pp. 256–294, 2019. [6] Y.-H. Jia, Y. Mei, and M. Zhang, “A bilevel ant colony optimization algorithm for capacitated electric vehicle routing problem,” IEEE transactions on cybernetics, vol. 52, no. 10, pp. 10 855–10 868, 2021. [7] ——, “Confidence-based ant colony optimization for capacitated electric vehicle routing problem with comparison of different encoding schemes,” IEEE Transactions on Evolutionary Computation, vol. 26, no. 6, pp. 1394–1408, 2022. [8] E. Rodrı́guez-Esparza, A. D. Masegosa, D. Oliva, and E. Onieva, “A new hyper-heuristic based on adaptive simulated annealing and reinforcement learning for the capacitated electric vehicle routing problem,” Expert Systems with Applications, vol. 252, p. 124197, 2024. [9] Y. Chen, J. Xue, Y. Zhou, and Q. Wu, “An efficient threshold acceptancebased multi-layer search algorithm for capacitated electric vehicle routing problem,” IEEE Transactions on Intelligent Transportation Systems, vol. 25, no. 6, pp. 5867–5879, 2024. [10] M. Mavrovouniotis, C. Menelaou, S. Timotheou, G. Ellinas, C. Panayiotou, and M. Polycarpou, “A benchmark test suite for the electric capacitated vehicle routing problem,” in 2020 IEEE Congress on evolutionary computation (CEC). IEEE, 2020, pp. 1–8. [11] A. Montoya, C. Guéret, J. E. Mendoza, and J. G. Villegas, “The electric vehicle routing problem with nonlinear charging function,” Transportation Research Part B: Methodological, vol. 103, pp. 87–110, 2017.
JOURNAL TITLE
[12] E. K. Burke and Y. Bykov, “The late acceptance hill-climbing heuristic,” European Journal of Operational Research, vol. 258, no. 1, pp. 70–78, 2017. [13] G. H. Fonseca, H. G. Santos, and E. G. Carrano, “Late acceptance hillclimbing for high school timetabling,” Journal of Scheduling, vol. 19, no. 4, pp. 453–465, 2016. [14] A. L. Bolaji, A. F. Bamigbola, and P. B. Shola, “Late acceptance hill climbing algorithm for solving patient admission scheduling problem,” Knowledge-Based Systems, vol. 145, pp. 197–206, 2018. [15] X. Zuo, Y. Xiao, M. You, I. Kaku, and Y. Xu, “A new formulation of the electric vehicle routing problem with time windows considering concave nonlinear charging function,” Journal of Cleaner Production, vol. 236, p. 117687, 2019. [16] S. R. Kancharla and G. Ramadurai, “Electric vehicle routing problem with non-linear charging and load-dependent discharging,” Expert Systems with Applications, vol. 160, p. 113714, 2020. [17] H. Tahami, G. Rabadi, and M. Haouari, “Exact approaches for routing capacitated electric vehicles,” Transportation Research Part E: Logistics and Transportation Review, vol. 144, p. 102126, 2020. [18] E. Lam, G. Desaulniers, and P. J. Stuckey, “Branch-and-cut-and-price for the electric vehicle routing problem with time windows, piecewiselinear recharging and capacitated recharging stations,” Computers & Operations Research, vol. 145, p. 105870, 2022. [19] G. M. Nafstad, G. Desaulniers, and M. Stålhane, “Branch-price-and-cut for the electric vehicle routing problem with heterogeneous recharging technologies and nonlinear recharging functions,” Transportation Science, vol. 59, no. 3, pp. 628–646, 2025. [20] S. Caillard and R. B. Chabane, “Evolutionary-based ant system algorithm to solve the dynamic electric vehicle routing problem.” in ICORES, 2024, pp. 285–293. [21] H. Yang, S. Yang, Y. Xu, E. Cao, M. Lai, and Z. Dong, “Electric vehicle route optimization considering time-of-use electricity price by learnable partheno-genetic algorithm,” IEEE Transactions on smart grid, vol. 6, no. 2, pp. 657–666, 2015. [22] L. Zhen, Z. Xu, C. Ma, and L. Xiao, “Hybrid electric vehicle routing problem with mode selection,” International Journal of Production Research, vol. 58, no. 2, pp. 562–576, 2020. [23] Á. Felipe, M. T. Ortuño, G. Righini, and G. Tirado, “A heuristic approach for the green vehicle routing problem with multiple technologies and partial recharges,” Transportation Research Part E: Logistics and Transportation Review, vol. 71, pp. 111–128, 2014. [24] D. Goeke, “Granular tabu search for the pickup and delivery problem with time windows and electric vehicles,” European Journal of Operational Research, vol. 278, no. 3, pp. 821–836, 2019. [25] M. Seyfi, M. Alinaghian, E. Ghorbani, B. Çatay, and M. S. Sabbagh, “Multi-mode hybrid electric vehicle routing problem,” Transportation Research Part E: Logistics and Transportation Review, vol. 166, p. 102882, 2022. [26] M. Keskin and B. Çatay, “A matheuristic method for the electric vehicle routing problem with time windows and fast chargers,” Computers & operations research, vol. 100, pp. 172–188, 2018. [27] P. C. Nolz, N. Absi, D. Feillet, and C. Seragiotto, “The consistent electric-vehicle routing problem with backhauls and charging management,” European Journal of Operational Research, vol. 302, no. 2, pp. 700–716, 2022. [28] Y. Marinakis, A. Migdalas, and P. M. Pardalos, “A new bilevel formulation for the vehicle routing problem and a solution method using a genetic algorithm,” Journal of Global Optimization, vol. 38, pp. 555– 580, 2007. [29] W. Tu, Z. Fang, Q. Li, S.-L. Shaw, and B. Chen, “A bi-level voronoi diagram-based metaheuristic for a large-scale multi-depot vehicle routing problem,” Transportation Research Part E: Logistics and Transportation Review, vol. 61, pp. 84–97, 2014. [30] B. E. Gillett and L. R. Miller, “A heuristic algorithm for the vehicledispatch problem,” Operations research, vol. 22, no. 2, pp. 340–349, 1974. [31] Y. Zhou, Y. Kou, and M. Zhou, “Bilevel memetic search approach to the soft-clustered vehicle routing problem,” Transportation Science, vol. 57, no. 3, pp. 701–716, 2023. [32] C. Prins, P. Lacomme, and C. Prodhon, “Order-first split-second methods for vehicle routing problems: A review,” Transportation Research Part C: Emerging Technologies, vol. 40, pp. 179–200, 2014. [33] Y. Qin and J. Chen, “A confidence-based bilevel memetic algorithm with adaptive selection scheme for capacitated electric vehicle routing problem,” in 2024 IEEE Congress on Evolutionary Computation (CEC). IEEE, 2024, pp. 1–10.
14
[34] C. Prins, “A simple and effective evolutionary algorithm for the vehicle routing problem,” Computers & operations research, vol. 31, no. 12, pp. 1985–2002, 2004. [35] D. Woller, V. Kozák, M. Kulich, and L. Přeučil, “Variable neighborhood search for the electric vehicle routing problem,” arXiv preprint arXiv:2511.09570, 2025. [36] V. Q. Hien, T. C. Dao, and H. T. T. Binh, “A greedy search based evolutionary algorithm for electric vehicle routing problem,” Applied Intelligence, vol. 53, no. 3, pp. 2908–2922, 2023. [37] M. López-Ibáñez, J. Dubois-Lacoste, L. P. Cáceres, M. Birattari, and T. Stützle, “The irace package: Iterated racing for automatic algorithm configuration,” Operations Research Perspectives, vol. 3, pp. 43–58, 2016. [38] F. G. Lobo, M. Bazargani, and E. K. Burke, “A cutoff time strategy based on the coupon collector’s problem,” European Journal of Operational Research, vol. 286, no. 1, pp. 101–114, 2020.