Paper deep dive
PyVRP$^+$: LLM-Driven Metacognitive Heuristic Evolution for Hybrid Genetic Search in Vehicle Routing Problems
Manuj Malik, Jianan Zhou, Shashank Reddy Chirra, Zhiguang Cao
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 97%
Last extracted: 4/10/2026, 4:12:03 AM
Summary
The paper introduces Metacognitive Evolutionary Programming (MEP), a framework that enhances Large Language Models (LLMs) to act as strategic discovery agents for evolving metaheuristics. By implementing a structured 'Reason-Act-Reflect' cycle and domain-aware initialization, MEP improves the Hybrid Genetic Search (HGS) algorithm for various Vehicle Routing Problem (VRP) variants, achieving up to 2.70% better solution quality and 45% faster runtimes.
Entities (5)
Relation Signals (3)
MEP → utilizes → LLM
confidence 100% · elevates the LLM to a strategic discovery agent
MEP → optimizes → HGS
confidence 95% · By applying MEP to evolve core components of the state-of-the-art Hybrid Genetic Search (HGS) algorithm
PyVRP → implements → HGS
confidence 90% · PyVRP is a state-of-the-art, open-source implementation of HGS specific to VRPs
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Designing high-performing metaheuristics for NP-hard combinatorial optimization problems, such as the Vehicle Routing Problem (VRP), remains a significant challenge, often requiring extensive domain expertise and manual tuning. Recent advances have demonstrated the potential of large language models (LLMs) to automate this process through evolutionary search. However, existing methods are largely reactive, relying on immediate performance feedback to guide what are essentially black-box code mutations. Our work departs from this paradigm by introducing Metacognitive Evolutionary Programming (MEP), a framework that elevates the LLM to a strategic discovery agent. Instead of merely reacting to performance scores, MEP compels the LLM to engage in a structured Reason-Act-Reflect cycle, forcing it to explicitly diagnose failures, formulate design hypotheses, and implement solutions grounded in pre-supplied domain knowledge. By applying MEP to evolve core components of the state-of-the-art Hybrid Genetic Search (HGS) algorithm, we discover novel heuristics that significantly outperform the original baseline. By steering the LLM to reason strategically about the exploration-exploitation trade-off, our approach discovers more effective and efficient heuristics applicable across a wide spectrum of VRP variants. Our results show that MEP discovers heuristics that yield significant performance gains over the original HGS baseline, improving solution quality by up to 2.70\% and reducing runtime by over 45\% on challenging VRP variants.
Tags
Links
- Source: https://arxiv.org/abs/2604.07872v1
- Canonical: https://arxiv.org/abs/2604.07872v1
Trouble viewing inline? Open PDF directly →
Full Text
86,900 characters extracted from source content.
Expand or collapse full text
PyVRP + : LLM-Driven Metacognitive Heuristic Evolution for Hybrid Genetic Search in Vehicle Routing Problems Manuj Malik* Singapore Management University Singapore, Singapore manujm@smu.edu.sg Jianan Zhou* Nanyang Technological University Singapore, Singapore jianan004@e.ntu.edu.sg Shashank Reddy Chirra University of Oxford Oxford, England shashank@robots.ox.ac.uk Zhiguang Cao Singapore Management University Singapore, Singapore zgcao@smu.edu.sg ABSTRACT Designing high-performing metaheuristics for NP-hard combinato- rial optimization problems, such as the Vehicle Routing Problem (VRP), remains a significant challenge, often requiring extensive domain expertise and manual tuning. Recent advances have demon- strated the potential of large language models (LLMs) to automate this process through evolutionary search. However, existing meth- ods are largely reactive, relying on immediate performance feedback to guide what are essentially black-box code mutations. Our work departs from this paradigm by introducing Metacognitive Evolu- tionary Programming (MEP), a framework that elevates the LLM to a strategic discovery agent. Instead of merely reacting to per- formance scores, MEP compels the LLM to engage in a structured Reason-Act-Reflect cycle, forcing it to explicitly diagnose failures, formulate design hypotheses, and implement solutions grounded in pre-supplied domain knowledge. By applying MEP to evolve core components of the state-of-the-art Hybrid Genetic Search (HGS) al- gorithm, we discover novel heuristics that significantly outperform the original baseline. By steering the LLM to reason strategically about the exploration-exploitation trade-off, our approach discov- ers more effective and efficient heuristics applicable across a wide spectrum of VRP variants. Our results show that MEP discovers heuristics that yield significant performance gains over the origi- nal HGS baseline, improving solution quality by up to 2.70% and reducing runtime by over 45% on challenging VRP variants. KEYWORDS Metaheuristic, Vehicle Routing Problem, Combinatorial Optimiza- tion, Large Language Model, Evolutionary Algorithm ACM Reference Format: Manuj Malik*, Jianan Zhou*, Shashank Reddy Chirra, and Zhiguang Cao. 2026. PyVRP + : LLM-Driven Metacognitive Heuristic Evolution for Hybrid Genetic Search in Vehicle Routing Problems. In Proc. of the 25th International Conference on Autonomous Agents and Multiagent Systems (AAMAS 2026), ∗ Corresponding authors. This work is licensed under a Creative Commons Attribution Inter- national 4.0 License. Proc. of the 25th International Conference on Autonomous Agents and Multiagent Systems (AAMAS 2026), C. Amato, L. Dennis, V. Mascardi, J. Thangarajah (eds.), May 25 – 29, 2026, Paphos, Cyprus.© 2026 International Foundation for Autonomous Agents and Multiagent Systems (w.ifaamas.org). https://doi.org/10.65109/ Paphos, Cyprus, May 25 – 29, 2026, IFAAMAS, 18 pages. https://doi.org/10. 65109/ 1 INTRODUCTION Vehicle routing problems (VRPs) represent the foundation of combi- natorial optimization (CO) with profound implications for logistics, supply chain management, and transportation. Given their NP-hard nature, exact methods are intractable for large-scale, real-world instances, leading to the widespread adoption of metaheuristics. These high-level algorithmic frameworks are designed to efficiently explore the solution space for near-optimal solutions. Among these, Hybrid Genetic Search (HGS) has emerged as a State-Of-The-Art (SOTA) approach, demonstrating exceptional performance across numerous VRP variants [24]. Despite its effectiveness, HGS remains heavily dependent on human-crafted algorithmic components. Like many metaheuristics, its performance hinges on carefully tuned internal mechanisms, such as parent selection, survivor selection, and penalty adjustment, which are designed through tedious trial and error and require significant domain expertise. Recently, Large Language Models (LLMs) have emerged as a promising tool for automatic heuristic design and scientific discov- ery [18]. By leveraging their ability to reason over code, data, and mathematical structure, LLMs have the potential to surpass human- engineered solutions across a wide range of domains. For instance, FunSearch [19] integrated LLMs into an evolutionary framework to discover novel solutions to mathematical problems, while Al- phaEvolve [18] extended this paradigm to generate strategies that outperform human-designed algorithms in domains such as ma- trix multiplication, kernel engineering, and beyond. Beyond raw performance, LLM-driven evolution condenses the costly manual design cycle—often requiring weeks of domain expert effort—into a manageable automated process. Similar progress has also been made in the domain of CO. Evo- lution of Heuristics (EoH) [13] demonstrated that LLMs can ef- fectively automate heuristic design within an evolutionary loop by co-evolving solutions at two levels of abstraction: high-level thoughts expressed in natural language and their corresponding executable code. ReEvo [32] further leveraged “verbal gradients” derived from performance feedback to iteratively refine heuristics. While promising, these approaches mainly focus on simple heuris- tics for standard CO problems such as the traveling salesman prob- lem (TSP), and advancing traditional solvers beyond human capabil- ity remains a challenge. While these pioneering works successfully arXiv:2604.07872v1 [cs.NE] 9 Apr 2026 demonstrate the potential of LLM-driven evolution, they often treat the LLM as a reactive code generator, relying predominantly on immediate performance feedback to steer the search. This paper ar- gues that such an approach underutilizes the reasoning capabilities of modern LLMs. Our core contribution is therefore not only the use of an LLM for evolution but also the introduction of Metacogni- tive Evolutionary Programming (MEP), a framework that shifts the paradigm from reactive mutation to proactive, structured discovery. MEP elevates the LLM from a simple code mutator to a strategic dis- covery agent by embedding two key innovations: a Domain-Aware Initialization phase that grounds the model in established strategic principles, and a mandatory Reason-Act-Reflect cycle that compels the LLM to explicitly diagnose failures, formulate testable design hypotheses, and critically self-assess its own output. This metacog- nitive scaffolding is the primary driver of our results, enabling the discovery of heuristics that not only perform high but also advance SOTA in VRP solving with minimal human intervention. By framing the evolutionary search as a hypothesis-driven dis- covery process, our aim is to show that heuristic evolution can sys- tematically uncover novel heuristics that outperform their human- designed counterparts. Our work builds on the framework intro- duced by [1], which models human-like metacognitive reasoning. We structure our methodology in two key phases: •Domain-Aware Initialization: We equip the LLM with do- main knowledge via three critical components: common pit- falls (퐾 푝 ) that undermine the target component, proven mit- igation strategies (퐾 푠 ) to address these pitfalls, and problem- specific traps (퐾 푡 ) unique to VRP solving. This knowledge foundation primes the LLM with the necessary domain un- derstanding before the evolutionary search begins. •Reason-Act-Reflect Cycle: Each generation in the evolu- tionary loop of MEP follows a structured cognitive cycle. The LLM is prompted to: 1) Reason: Explicitly diagnose the shortcomings of parent heuristics. 2) Act: Formulate a clear design hypothesis and implement it as a new code-based heuristic. 3) Reflect: Assess its own generated code by artic- ulating its rationale and potential limitations directly within the code’s documentation. We apply MEP to the challenging task of discovering high- performance components for the HGS algorithm, using the open- source PyVRP library [26] as our testbed. The discovered collection of high-performing heuristics, which we term PyVRP + , demon- strates that a structured, hypothesis-driven evolution can produce novel HGS variants that advance the state-of-the-art in VRP solving. The results show that MEP is capable of evolving heuristics that deliver performance gains of up to 2.70%, while simultaneously reducing runtime by over 45% on challenging VRP variants. These improvements underscore MEP’s potential as a compelling alterna- tive to hand-crafted HGS components and represent an encouraging step toward automated algorithm discovery in CO. 2 RELATED WORK Recent research has increasingly adopted LLMs to automate the design of heuristics for CO. Pioneering this direction, FunSearch [19] demonstrated that LLMs could effectively generate core heuris- tic components, circumventing the need for manual redesign of entire algorithms. Similarly, Liu et al. [12]employed LLMs to gen- erate key evolutionary algorithm components, such as crossover and mutation operators. Further advancing this line of work, the Evolution of Heuristics (EoH) framework [13] introduced a method for iteratively refining heuristic seed functions through evolution- ary computation. Most of these approaches have adopted evolu- tionary frameworks, integrating advanced techniques to enhance their efficacy. These techniques include reflective feedback mecha- nisms [32], predictive modeling for performance estimation [27], and strategies that encourage heuristic diversity [8]. This para- digm has successfully extended beyond single-objective scenar- ios, addressing multi-objective optimization problems [31] and di- verse CO domains including SAT [21], MILP [33], and dynamic Job Shop Scheduling Problems (JSSP) [10]. Moving beyond single- agent approaches, recent studies [21,29] have explored multi-agent frameworks, achieving enhanced performance through collabora- tive heuristic generation. To overcome limitations associated with population-based methods, Zheng et al. [35]introduced a Monte Carlo Tree Search (MCTS) strategy, integrating LLM-generated heuristics to boost exploration efficiency. More recent advance- ments focus on reducing manual intervention. Thach et al. [23] proposed a novel end-to-end framework based on problem reduc- tion, which enables heuristic methods to function independently of predefined general algorithmic templates. Similarly, Shi et al. [20]developed a meta-optimization framework that systematically discovers robust optimizers for a range of heuristic design tasks, pro- moting broader exploration and operating at a higher abstraction level. Furthermore, Šurina et al. [22]suggested integrating evolu- tionary heuristic generation with reinforcement learning-based fine-tuning, shifting LLM usage from static heuristic generation towards dynamic refinement of both heuristics and the underlying models. However, these approaches primarily focus on evolving simple heuristics or individual components thereof, rather than tackling sophisticated and well-established traditional solvers. Most existing work targets constructive heuristics for basic problem vari- ants or designs standalone neural policies that bypass traditional optimization entirely [2,5]. Concurrently, LLMs have also been explored as end-to-end CO solvers [11], for refining HGS via RL- finetuned LLMs [37], and within agentic frameworks for complex VRPs [34]. In this paper, we bridge this gap by demonstrating that LLM-guided evolution can enhance the performance of a SOTA solver like HGS across multiple challenging VRP variants, moving beyond proof-of-concept demonstrations to practical algorithmic improve- ments. For a comprehensive review of the use of LLMs in CO, we direct readers to [7]. LLM-guided scientific discovery is rapidly expanding beyond combinatorial optimization. For instance, Chen et al. [6]apply LLMs to Neural Architecture Search, Goldie et al. [9]use them to discover novel loss functions and optimizers for reinforcement learning, and Nadimpalli et al. [16]explore their use in evolving activation func- tions. While most of these methods focus on evolving function-level programs, Novikov et al. [18]scale this approach to entire files with hundreds of lines of code, enabling the discovery of novel algo- rithms for matrix multiplication, solutions to open mathematical problems, kernel optimization, and more. Although these methods often outperform existing human-designed baselines within their domains, more ambitious efforts, such as AI Scientist [14,28], aim for fully autonomous end-to-end scientific discovery with minimal human intervention. While such systems demonstrate impressive capabilities, they still fall short of matching human-level perfor- mance (e.g., producing research papers that are accepted at top-tier conferences). To support progress in this direction, MLGym [17] introduces a suite of benchmark environments specifically designed for training and evaluating models on scientific discovery tasks, signaling a promising future for this emerging field. 3 PRELIMINARIES This section presents the necessary background on Vehicle Routing Problems and the Hybrid Genetic Search algorithm. 3.1 Vehicle Routing Problem The VRP is a class of combinatorial optimization problems focused on finding optimal sets of routes for a fleet of vehicles to serve a set of customers. A VRP is typically defined on a graph퐺=(푉,퐴), where푉is a set of nodes and퐴is a set of arcs. The node set푉 is partitioned into a set of depots and a set of customers. Each arc(푖, 푗) ∈ 퐴is associated with a travel cost푐 푖,푗 . The objective is to determine a set of vehicle routes, starting and ending at a depot, that visit all customers while minimizing a specific objective function, usually the total travel cost. To systematically handle the diverse landscape of VRP variants, we adopt the composable attribute framework used in [3,4,36]. This allows any variant to be defined by a combination of fundamental properties, facilitating structured analysis and evaluation. These attributes are categorized as follows: •Node-specific attributes, which define local requirements for each customer or depot. The primary examples are: –Demand: The quantity of goods for pickup or delivery at a node푖, denoted as푞 푖 . This directly influences vehicle load throughout a route. –Time Windows (TW): The required service interval[푎 푖 ,푏 푖 ] for a customer, which imposes critical temporal constraints on the route schedule. – Service Time: The duration푠 푖 required to complete service at a node, which consumes time and affects the feasibility of subsequent stops. •Global attributes, which apply to entire routes or the problem as a whole. Key global constraints include: – Capacity (C): The maximum load푄 푘 that a vehicle can carry, limiting the total demand of customers on a single route. –Duration Limits (L): A maximum allowable travel distance or time for a route, ensuring operational constraints are met. –Backhauling (B/MB): Rules governing routes with both de- liveries (linehauls) and pickups (backhauls), often enforcing precedence constraints. –Open Routes (O): A structural variant where vehicles are not required to return to the depot after their final delivery. –Multi-Depot (MD): Scenarios involving multiple starting and ending points for vehicles, adding a layer of assignment com- plexity. This unified taxonomy provides a clear and systematic method for representing complex VRPs, enabling the evaluation of algorith- mic components across a wide spectrum of problem types. 3.2 VRP Variants Considered In this work, we evaluate our approach on a diverse set of VRPs to test the generalization of the discovered heuristics. We consider the following 7 VRP variants: •Traveling Salesman Problem (TSP): The TSP represents the simplest case of the VRP family, involving a single vehicle (or salesman) that must visit every customer exactly once and return to the starting depot. Despite its conceptual simplicity, the TSP is NP-hard and serves as a building block for more complex routing problems. The objective is to minimize the total travel distance or cost of the tour. •Capacitated VRP (CVRP): The CVRP extends the basic VRP by introducing vehicle capacity constraints. Each ve- hicle푘has a finite capacity푄 푘 , and each customer푖has a demand푞 푖 that must be satisfied. The constraint ensures that the total demand of all customers served by a single route cannot exceed the vehicle’s capacity: Í 푖∈푅 푘 푞 푖 ≤ 푄 푘 , where푅 푘 represents the set of customers served by vehicle푘. This variant is fundamental in logistics applications where vehicles have physical limitations. •Generalized Vehicle Routing Problem (GVRP): In this variant, graph nodes are grouped into mutually exclusive clusters, and serving any single node within a cluster satisfies the service requirement for the entire cluster. This creates a combinatorial challenge where the algorithm must select both which clusters to visit and which specific node within each cluster to serve. The generalized structure is applicable to scenarios where customers have multiple service locations or where service points can substitute for each other. •Multi-Depot VRP with Time Windows (MDVRPTW): This variant combines the complexity of multiple depot lo- cations with time window constraints. Each vehicle is based at one of several depots and must serve customers within their specified time windows[푎 푖 ,푏 푖 ]. The challenge lies in optimally assigning customers to depots and vehicles while respecting both spatial proximity and temporal feasibility. Each depot typically has its own fleet of vehicles, and routes must start and end at the same depot. This variant is partic- ularly relevant for logistics companies with multiple distri- bution centers serving geographically dispersed customers with delivery time requirements. •Prize-Collecting VRP with Time Windows (PCVRPTW): Unlike traditional VRP variants where all customers must be visited, the PCVRPTW allows for selective customer service. Each customer푖has an associated prize휋 푖 collected when visited, and the objective is to maximize the total collected prize minus the total travel cost, subject to a total cost bud- get퐵. Customers have time windows, and the challenge is to identify the most profitable subset of customers to serve within the available resources. •VRP with Backhauls (VRPB): The VRPB addresses scenar- ios where customers are divided into two categories: linehaul customers (requiring delivery from the depot) and backhaul customers (requiring pickup to the depot). A critical con- straint is that all linehaul customers must be served before any backhaul customers on the same route. This precedence Phase I: Domain-Aware InitializationPhase I: LLM Generation Large Language Model (LLM) REASON diagnosis, hypothesis ACT code gen, valid. check REFLECT rationale, offspring Phase I: Solution Selection Population Update Generation Evolutionary Operators (Improved) Generation Fitness Evaluation Structured Prompt Repeat Inputs HGS-Parent Selection Base heuristics: Evaluation instances: Metacognitive Mandate Define optimization context Selection Strategy Best-of- Selection Figure 1: An overview of MEP, which outlines a process for evolving high-performing heuristics. The process continues until convergence criteria or a maximum number of generations is reached, ultimately returning the best-performing heuristic. constraint reflects practical considerations in logistics, where vehicles must deliver goods before they have space to collect returns or recyclables. • VRP with Time Windows (VRPTW): Each customer푖 must be visited within a specified time window[푎 푖 ,푏 푖 ], where 푎 푖 is the earliest service time and푏 푖 is the latest service time. If a vehicle arrives before푎 푖 , it must wait until the service can begin. Arrival after푏 푖 renders the solution infeasible. This variant is particularly relevant in applications such as package delivery, where customers have preferred delivery times, or in service industries with appointment scheduling. A formal mathematical definition of the base VRP is provided in Appendix A. 3.3 Hybrid Genetic Search HGS is a powerful metaheuristic framework that combines the global search capabilities of a genetic algorithm with the fine-tuning power of a highly effective local search. The general workflow of HGS is as follows: •Initialization: A diverse initial population is generated, often using problem-specific constructive heuristics. •Parent Selection: Two parent solutions (푃 1 ,푃 2 ) are selected us- ing a fitness-biased mechanism that balances solution quality (elitism) and genetic diversity, typically measured by solution distance metrics. •Crossover: An offspring solution (퐶) is generated using prob- lem specific recombination operators, such as Selective Route- EXchange (SREX) for VRPs, that effectively preserve mean- ingful structural characteristics of the parent solutions. •Education (Local Search): The offspring undergoes intensive improvement via a problem-tailored local search procedure. For VRPs, this typically involves neighborhood moves like relocate, swap, and 2-opt operators. This hybridization with local search is fundamental to HGS’s performance. • Survivor Selection: The educated offspring is considered for population inclusion. If population size limits are exceeded, removal prioritizes either low-quality solutions or those ex- hibiting high similarity to others, maintaining both quality and diversity through a distance-based diversity manage- ment mechanism. •Diversification: When stagnation is detected, the search restarts through partial population reinitialization or injection of new diverse solutions, preventing premature convergence. PyVRP is a state-of-the-art, open-source implementation of HGS specific to VRPs, written in Python with performance-critical com- ponents in C++. Its modular architecture enables isolated replace- ment of core components (e.g., parent selection, local search opera- tors) with custom implementations, creating an ideal experimental platform for heuristic evolution via LLM-generated functions. 4 METACOGNITIVE EVOLUTIONARY PROGRAMMING We propose Metacognitive Evolutionary Programming (MEP), an evolutionary framework inspired by Bai et al. [1]that emulates human metacognitive reasoning by guiding the LLM through a structured and strategic design process. As illustrated in Figure 1, MEP enables the LLM to evolve heuristics through a human-like reason–act–reflect cycle, augmented with domain knowledge and guided by a planning strategy that fosters strategic improvement. The central thesis of MEP is that engaging the LLM in a struc- tured, hypothesis-driven process makes code evolution more effi- cient and more likely to yield novel, high-performing heuristics. The framework is composed of two main phases: a one-time Domain- Aware Initialization, and an iterative Reason–Act–Reflect cycle. 4.1 Phase 1: Domain-Aware Initialization Before the evolutionary search begins, MEP equips the LLM with comprehensive knowledge about the target component. This plan- ning phase provides three critical knowledge types: • Common Pitfalls (퐾 푝 ): Known failure modes specific to the component being evolved. •Mitigation Strategies (퐾 푠 ): Proven techniques to address iden- tified pitfalls. • Problem-Specific Traps (퐾 푡 ): Domain-specific challenges (e.g., in VRP solving). This knowledge foundation serves as domain expertise that guides the LLM’s reasoning throughout the evolutionary process, ensuring that generated heuristics are informed by established algorithmic principles. 4.2 Phase 2: Reason-Act-Reflect Cycle Each generation of a new heuristic in MEP follows a structured three-step cognitive process, called the Reason-Act-Reflect cycle. Step 1: REASON (Diagnosis and Hypothesis). At the start of each evolutionary timestep, two parent heuristics are selected from the current population. Before generating any code, the LLM is prompted to perform a structured analysis of the parents, encour- aging deliberate reasoning over random mutation. • Diagnosis: The LLM begins by analyzing the provided parent heuristics, i.e.,code 1 ,code 2 , along with their performance scores and any available feedback. It is explicitly prompted to diagnose their weaknesses, guided by the principles intro- duced in Phase 1. •Design Hypothesis: Building on the identified weaknesses, the LLM must formulate a concise, one-sentence Design Hy- pothesis that proposes a specific improvement for the new heuristic (code 3 ). This step is crucial, as it ensures the LLM explicitly reasons about how to address the limitations of its predecessors, rather than simply combining elements of existing code. Step 2: ACT (Code Implementation). Once the hypothesis is formu- lated, the LLM is tasked with implementing it. This step is governed by strict guardrails including a fixed function signature, approved libraries, and a required output format, to ensure the generated code is syntactically valid and fully compatible with the PyVRP evaluation sandbox. Step 3: REFLECT (Rationale and Self-Critique). After generat- ing the code, the LLM engages in a final metacognitive step: self- reflection. In this phase, it critically evaluates its prior reasoning, hypothesis formulation, and implementation choices. This self- assessment is recorded and carried forward to inform the gener- ation of future offspring, enabling improvement over subsequent evolutionary cycles. The Reason-Act-Reflect cycle enforces a principled generation process, ensuring that each new heuristic emerges from structured reasoning rather than arbitrary mutation. Figure 1 outlines the process for evolving high-performing heuristics using the pro- posed method. A high-level prompt structure is presented in Figure 2, breaking-down the LLM interaction into three distinct phase- aligned segments. Planning provides high-level strategic context, explicitly injecting domain knowledge to ground the LLM’s initial reasoning in established principles. Reasoning directs the LLM to perform a structured diagnosis of parent heuristics’ weaknesses and mandates the formulation of a concise, testable Design Hypoth- esis before any code generation, ensuring deliberate improvement intent. Reflection requires the LLM to self-critique its generated Role: AI Optimization Researcher (Vehicle Routing Problem) 1) Planning Phase • Pitfalls (퐾 푝 ): List potential failures in VRP parent selection... • Strategies (퐾 푠 ): Propose countermeasures for each pitfall... • Hidden Traps (퐾 푡 ): Identify misleading problem instance features... Overall Objective: Write anovelPython functionselect_parentsfor a Hybrid Genetic Search algorithm, guided by the planning insights. Context: Background on parent selection in genetic algorithms for VRP, including standard approaches (e.g., elitist, tournament). ... Common Existing Approaches: Standard parent selection methods in genetic algorithms include elitist selection (top-ranked individuals), tournament selection (best from random subset), fitness-proportionate selection, etc. ... 2) Reasoning Phase •Analyze Patterns: Assess exploration/exploitation bias in provided selec- tors. • Assess Shortcomings: Evaluate selectors against the 퐾 푝 /퐾 푡 list. • Design Sketch: Outline the proposed new selector, sel3, in prose. Your Task in This Interaction: (1) Analyze sel1 and sel2. (2) Design and implement sel3 = select_parents(...). (3) Implementation Requirements: • Inputs: population, rng, cost_evaluator, k=2 • Output: tuple[Solution, Solution] • Behavior: Balance feasibility, cost, and diversity. Reference Python Classes: class CostEvaluator: ... class Solution: ... class ProblemData: ... # ... (other PyVRP class definitions) 3) Reflection Phase (1) How does the final code address the initial pitfalls (퐾 푝 )? (2) Did new challenges emerge during coding? (3) Propose two concrete revisions for further improvement. Response Format: A single Python code block containing only the required function and imports. Dynamic Inputs for a Given Task: • Selector 1: $code1 • Performance: Cost: $score1, Feedback: $feedback1 • Selector 2: $code2 • Performance: Cost: $score2, Feedback: $feedback2 Figure 2: The prompt template used to guide the LLM. It consists of distinct instructions for planning, reasoning, and reflection, providing both high-level guidance and detailed technical specifications. code and the preceding reasoning/hypothesis within the code doc- umentation, embedding metacognitive assessment directly into the output. Each segment combines high-level instructions (e.g., “Diag- nose the weaknesses...”) with detailed technical specifications (e.g., fixed function signatures, approved libraries, required output for- mat) to enforce both strategic coherence and syntactic correctness. The full prompt is provided in Appendix B. 5 MEP CASE STUDY ON HGS We apply MEP to evolve three critical components within the HGS algorithm. This includesselect_parentsfor fitness-diversity bal- anced selection,select_survivorsfor population management, andupdate_penaltiesfor constraint handling of infeasible solu- tions. These operators directly govern the exploration-exploitation trade-off of the genetic algorithm. In this context, exploitation refers Table 1: Average cost comparison of the original PyVRP components and the best components evolved by MEP on 100 TSP instances. Lower cost values indicate better performance. Component Evolved Baseline Cost Best Evolved Cost Improvement select_parents7 942 16677648152.23 % select_survivors7 942 16678868520.69 % update_penalties7 942 16679329730.12 % to selecting high-quality (low-cost) parents to refine the best-known solutions, while exploration involves selecting diverse parents to in- troduce new genetic material and prevent premature convergence to local optima. A better heuristic must balance these two com- peting objectives. To this end, we use the PyVRP library as our experimental sandbox. Thanks to its modular design, we can isolate and replace core HGS components (such as the parent selection operator) with LLM-generated functions. This approach creates PyVRP + , a collection of better-performing heuristics. Evolutionary Process: It is outlined as follows: • Base Population: The base population consists of default implementations of selected modules from PyVRP. •Evaluation: Each heuristic in the population is evaluated by running HGS on a set of VRP instances, with its fitness score defined as the negative average cost of the solution achieved. •Generation: In each generation, two parent heuristics and their performance scores are selected. These, together with domain knowledge, are provided to the MEP prompt to pro- duce a new offspring heuristic. • Population Update: The new heuristic is evaluated and considered for inclusion in the population. •Termination: The process runs for a fixed number of gen- erations (N), and the best heuristic discovered is reported. The evolutionary process runs for 10 generations. In each genera- tion, 10 offspring heuristics are produced, and the top 5 performers, i.e., selected from both the new offspring and the existing popula- tion, are retained for the next generation. To ensure the robustness of our findings, the entire evolutionary process for each component was conducted five times using different random seeds, and the best-performing heuristic from these runs is reported in our results. The LLM interactions were performed using the GPT-4.1 model. We set the temperature parameter to 1.0 to encourage creative and diverse heuristic designs, while all other parameters, such as top-p and max_tokens, were left at their default values. Benchmark Suite: We evaluate our approach on two datasets. First, we use the publicly available TSP100 dataset [30], a well- established benchmark in combinatorial optimization, to assess LLM-generated heuristics during the evolution process. Second, we evaluate the discovered HGS variants on a comprehensive bench- mark suite spanning six VRP variants. Each variant, sourced from the PyVRP library 1 , includes a dedicated folder containing 60 in- stances, ensuring statistically robust evaluation. The benchmark spans a range of problem complexities, i.e., from the standard Ca- pacitated VRP (CVRP) to more challenging formulations like the Multi-Depot VRP with Time Windows (MDVRPTW). The diversity of these variants, each with distinct constraint structures and ob- jectives, provides a comprehensive testing ground for evaluating the robustness and generalizability of evolved heuristics. 6 EXPERIMENTAL RESULT AND ANALYSIS Our primary goal is to demonstrate that MEP can discover better heuristic components for a SOTA VRP solver. Before presenting the results, it is crucial to distinguish between two different time costs: the one-time design cost of the evolutionary process and the recur- ring execution cost of the final heuristic. The performance tables in our analysis report the execution cost, which is the wall-clock time required for the final, evolved heuristic to solve a given VRP instance. This metric allows for a direct comparison against the baseline solver’s runtime. In contrast, the design cost is the upfront, one-time computational investment required to discover the heuris- tic. This involves the entire MEP process, including all LLM API calls and the evaluation of candidate solutions over multiple genera- tions. We consider this design cost analogous to the manual effort a human researcher would invest over weeks or months to develop a novel algorithm. For this work, the total design cost for discovering each final component was approximately 3–4 hours on a single NVIDIA H100 GPU. This involved approximately 100 LLM API calls per component at a total cost of $15–20 USD, substantially lower than the weeks of manual effort typically required. We posit that this one-time investment is a practical trade-off for producing novel heuristics that significantly improve solution quality and reduce subsequent execution costs. We use GPT-4.1 in the evolutionary search, as mentioned earlier, chosen for its strong performance and cost-effectiveness. The following standard metrics are used to evaluate solution quality: •Cost/Obj.: The total cost of a solution, typically representing the sum of travel distances or travel times across all routes. Lower values indicate better solutions. •Improvement: The relative percentage by which the method’s objective value is better than the best-known solution (BKS). It is computed as: Improvement= Obj BKS − Obj method Obj BKS × 100%(1) A positive value indicates that the method outperforms the BKS, while a negative value means it performs worse. 1 https://github.com/PyVRP/Instances/ Table 2: Performance comparison of the baseline PyVRP solver versus the PyVRP + equipped with the MEP-evolved operators across various VRP variants (average values are presented). Lower cost and shorter times indicate better performance. VariantBaseline Cost MEP-Evolved Cost Improvement Baseline Time (s) MEP-Evolved Time (s) CVRP64 450.0364236.300.33 %9.559.11 GVRP7 686 993.647532007.092.01 %0.871.89 MDVRPTW8955.118825.461.45 %12.439.74 PCVRPTW22 975.8522406.422.48 %7.924.19 VRPB83 696.2183146.060.66 %14.0713.85 VRPTW33 424.4832521.752.70 %17.0314.66 6.1 Hardware Settings All experiments are conducted on a high-performance computing cluster with the following specifications: • GPU: NVIDIA H100 GPU with dedicated allocation • CPU: 32-core processor allocation per job • Memory: 120 GB RAM per compute node •Cluster Configuration: Single-node allocation with GPU acceleration 6.2 Evolving Effective HGS Components Table 1 summarizes the average cost of the best-evolved heuristic for each target component compared to the PyVRP baseline on the TSP100 benchmark. In all cases, MEP successfully discovered an op- erator with a better cost. The most significant improvements were observed in theselect_parentsandselect_survivorscompo- nents, which saw a 2.23% and 0.69% improvement, respectively. This highlights MEP’s ability to innovate on complex, high-impact components of the metaheuristic. The cost-generation curve illus- trating the iterative improvement of the evolutionary process is shown in Figure 3. 02468 Generations 7770000 7780000 7790000 7800000 7810000 7820000 7830000 Cost TSP Cost Evolution Over Generations Figure 3: Evolution of average cost values over 10 generations for TSP optimization using MEP for parent selection, with averages computed across 100 instances. A key outcome of the MEP process is the emergence of sophisti- cated and novel heuristics that go beyond simple variations of exist- ing strategies. For example, the top-performingselect_parents operator, described by the LLM as a “Hybrid Adaptive Stratified Diversity-Pressure Selection”, divides the population into cost- based strata and applies tailored diversity pressures to each sub- group. The complete evolved code is provided in Appendix C. These results highlight that MEP’s hypothesis-driven framework enables the discovery of complex, well-motivated algorithmic innovations, far beyond simple parameter tuning or superficial modifications. To evaluate the synergistic potential of these evolved heuris- tics, we created an integrated HGS solver that combines the best- performing LLM-generated modules for parent selection, survivor selection, and penalty adjustment. This integrated solver repre- sents a more holistic test of MEP’s capabilities, assessing whether components evolved in isolation can cooperate effectively. As we will demonstrate in the following section, this integrated approach not only matches but often exceeds the performance of individ- ually evolved modules, achieving consistent cost reductions and significant runtime improvements on complex VRP variants. This finding validates that MEP can produce a suite of compatible, high- performing heuristics that collectively advance the state-of-the-art. 6.3 Integrated Solver Performance and Generalization To assess the robustness and overall performance of our approach, we evaluated the fully integrated HGS solver, i.e., equipped with the best evolved heuristics for parent selection, survivor selection, and penalty updates against the PyVRP baseline on a diverse set of VRP variants. The results, shown in Table 2, demonstrate a clear and consistent advantage for the MEP-evolved solver. The integrated solver achieves notable cost reductions across all variants, with particularly strong improvements on complex in- stances like VRPTW (2.70%), PCVRPTW (2.48% improvement) and GVRP (2.01% improvement). Furthermore, the evolved heuristics lead to significant computational efficiencies, reducing runtime by over 45% on PCVRPTW and 10% on MDVRPTW. While some vari- ants show a modest increase in runtime, the superior solution qual- ity underscores the effectiveness of the discovered strategies. This strong generalization performance confirms that MEP is capable of producing a synergistic set of heuristics that are not over-fitted to a single problem type and can enhance a state-of-the-art solver across a wide range of operational scenarios. 6.4 Generalization Across Evolved Components We also evaluated the generalization capabilities of the other best- evolved. Two of these experimental results are presented in Table 4, which also demonstrate positive, albeit more modest, improvements Table 3: Ablation study results for theselect_parentsop- erator, showing percentage improvement over the baseline. Average cost values are presented in the table. VariantBaseline Cost Full MEP Impr. MEP-noInit Impr. Reactive Evol. Impr. CVRP64 450.030.33 %0.28 %−0.20 % GVRP7 686 993.641.82 %0.98 %1.30 % MDVRPTW8955.111.49 %0.90 %−2.05 % PCVRPTW22 975.852.43 %2.36 %1.46 % VRPB83 696.210.57 %0.53 %−0.39 % VRPTW33 424.482.25 %2.24 %0.10 % Table 4: Generalization performance of individual versus combined MEP-evolved operators. The table shows the per- centage improvement over the HGS baseline, with the final column (PyVRP + ) representing a solver that integrates all individually evolved components. Variantselect_parents select_survivors update_penalties PyVRP + CVRP0.33 %0.05 %0.01 %0.33 % GVRP1.82 %0.08 %0.13 %2.01 % MDVRPTW1.49 %0.00 %0.16 %1.45 % PCVRPTW2.43 %0.15 %0.06 %2.48 % VRPB0.57 %0.04 %0.05 %0.66 % VRPTW2.25 %0.29 %0.35 %2.70 % over the baseline across various VRP instances. For example, the evolved ‘select_survivors‘ operator improves VRPTW performance by 0.29%, while the evolved ‘update_penalties‘ operator improves it by 0.35%. This indicates that the MEP framework is capable of discovering generally beneficial heuristics for multiple components of the HGS algorithm, even if the primary gains are concentrated in the highest-impact operators, such as parent selection. Notably, none of the evolved HGS components underperformed their baseline counterparts on any of the six VRP variants, despite not being explicitly evaluated on these during the evolution process. This highlights the robustness and generalizability of the discovered algorithmic strategies. 6.5 Ablation Studies To assess the individual contributions of key components in our framework, we performed a series of ablation studies focused on the select_parents module, chosen for its strong generalization across tasks. These experiments aim to isolate the effects of the Domain- Aware Initialization and the structured Reason–Act–Reflect cycle. Our ablation compares full MEP framework against two variants: •MEP-noInit: This version removes the Domain-Aware Initial- ization from the prompt. The LLM still performs the Reason- Act-Reflect cycle, but its reasoning is purely local, based only on the two parent heuristics provided in each generation. • Reactive Evolution: This variant removes the structured rea- soning components. The prompt is simplified, providing the parent heuristics and their scores and asking the LLM to gen- erate an improved version. This setup mimics prior reactive evolution methods such as EoH [13] and ReEvo [32]. The results in Table 3 reveal a clear performance hierarchy. The full MEP framework consistently achieves the best results. The MEP-noInit variant performs reasonably well but is worse than the full MEP, highlighting the benefit of the strategic guidance provided by the Domain-Aware Initialization. For example, on GVRP, the improvement drops from 1.82% with the full MEP to 0.98% with- out the Initialization. The Reactive Evolution variant performs the worst, often yielding only marginal gains or even performance degradation (e.g., -0.20% on CVRP), underscoring the importance of the structured reason-act-reflect cycle. 7 CONCLUSION In this work, we introduced Metacognitive Evolutionary Program- ming (MEP), a new paradigm for automated algorithm design that enhances LLM-driven evolution with strategic planning and struc- tured reasoning for the Vehicle Routing Problem (VRP). By ground- ing the LLM in domain knowledge and guiding it through a struc- tured Reason–Act–Reflect cycle, we elevate it from a basic code mutator to a strategic agent capable of deliberate heuristic dis- covery. Our application of MEP to evolve core components of the state-of-the-art HGS solver demonstrates a promising path toward discovering heuristics that are not only high-performing but also novel and strategically sound. Looking ahead, this work opens up several avenues for future research. Our work demonstrates that MEP can successfully enhance individual components of a state-of-the-art solver. More importantly, we have shown that these independently evolved heuristics can be successfully integrated, creating a holistically improved solver that generalizes across mul- tiple challenging VRP variants. Future work could explore the joint, simultaneous evolution of all components to uncover even deeper synergistic interactions, though this remains a computationally intensive challenge. While this presents significant challenges with the current LLM context length and computational costs, we be- lieve it is crucial for uncovering synergistic interactions between components. Our work demonstrates that a manageable one-time design cost can yield heuristics with lasting performance benefits, providing a strong incentive to tackle these future challenges. Fur- thermore, to rigorously test the generalizability of our discovered heuristics, we plan to evaluate their performance across the full set of 48 VRP variants covered in the original PyVRP. Overall, this work lays a foundation for more autonomous and intelligent algo- rithmic systems capable of scientific discovery, transcending the constraints imposed by manually designed algorithms. The code is available at https://github.com/ra-MANUJ-an/pyvrp-code. ACKNOWLEDGMENTS This research is supported by the National Research Foundation, Singapore under the AI Singapore Programme (AISG Award No: AISG3-RPGV-2025-017), and the Singapore Ministry of Education (MOE) Academic Research Fund (AcRF) Tier 1 grant. REFERENCES [1]Tian Bai, Yongwang Cao, Yan Ge, and Haitao Yu. 2025. MP: Endowing Large Language Models with Lateral Thinking. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39. 23460–23468. [2]Yoshua Bengio, Andrea Lodi, and Antoine Prouvost. 2021. Machine learning for combinatorial optimization: a methodological tour d’horizon. European Journal of Operational Research 290, 2 (2021), 405–421. [3]Federico Berto, Chuanbo Hua, Junyoung Park, Laurin Luttmann, Yining Ma, Fanchen Bu, Jiarui Wang, Haoran Ye, Minsu Kim, Sanghyeok Choi, Nayeli Gast Zepeda, André Hottung, Jianan Zhou, Jieyi Bi, Yu Hu, Fei Liu, Hyeonah Kim, Jiwoo Son, Haeyeon Kim, Davide Angioni, Wouter Kool, Zhiguang Cao, Qingfu Zhang, Joungho Kim, Jie Zhang, Kijung Shin, Cathy Wu, Sungsoo Ahn, Guojie Song, Changhyun Kwon, Kevin Tierney, Lin Xie, and Jinkyoo Park. 2025. RL4CO: an Extensive Reinforcement Learning for Combinatorial Optimization Benchmark. In SIGKDD Conference on Knowledge Discovery and Data Mining. [4]Federico Berto, Chuanbo Hua, Nayeli Gast Zepeda, André Hottung, Niels A Wouda, Leon Lan, Junyoung Park, Kevin Tierney, and Jinkyoo Park. 2025. RouteFinder: Towards Foundation Models for Vehicle Routing Problems. Trans- actions on Machine Learning Research 2025 (2025). [5]Quentin Cappart, Didier Chételat, Elias B Khalil, Andrea Lodi, Christopher Morris, and Petar Veličković. 2023. Combinatorial optimization and reasoning with graph neural networks. Journal of Machine Learning Research 24, 130 (2023), 1–61. [6]Angelica Chen, David Dohan, and David So. 2023. EvoPrompting: Language Models for Code-Level Neural Architecture Search. In Thirty-seventh Conference on Neural Information Processing Systems. https://openreview.net/forum?id= ifbF4WdT8f [7] Francesca Da Ros, Michael Soprano, Luca Di Gaspero, and Kevin Roitero. 2025. Large Language Models for Combinatorial Optimization: A Systematic Review. arXiv preprint arXiv:2507.03637 (2025). [8] Pham Vu Tuan Dat, Long Doan, and Huynh Thi Thanh Binh. 2025. Hsevo: Elevating automatic heuristic design with diversity-driven harmony search and genetic algorithm using llms. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39. 26931–26938. [9] Alexander David Goldie, Zilin Wang, Jakob Nicolaus Foerster, and Shimon Whiteson. 2025. How Should We Meta-Learn Reinforcement Learning Algo- rithms?. In Reinforcement Learning Conference. https://openreview.net/forum? id=jKzQ6af2DU [10]Jin Huang, Qihao Liu, Xinyu Li, Liang Gao, and Yue Teng. 2026. Automatic pro- gramming via large language models with population self-evolution for dynamic fuzzy job shop scheduling problem. IEEE Transactions on Fuzzy Systems (2026). [11]Xia Jiang, Yaoxin Wu, Minshuo Li, Zhiguang Cao, and Yingqian Zhang. 2025. Large language models as end-to-end combinatorial optimization solvers. In Advances in Neural Information Processing Systems. [12]Fei Liu, Xialiang Tong, Mingxuan Yuan, and Qingfu Zhang. 2023. Algorithm evolution using large language model. arXiv preprint arXiv:2311.15249 (2023). [13]Fei Liu, Tong Xialiang, Mingxuan Yuan, Xi Lin, Fu Luo, Zhenkun Wang, Zhichao Lu, and Qingfu Zhang. 2024. Evolution of Heuristics: Towards Efficient Automatic Algorithm Design Using Large Language Model. In International Conference on Machine Learning. [14] Chris Lu, Cong Lu, Robert Tjarko Lange, Jakob Foerster, Jeff Clune, and David Ha. 2024. The AI Scientist: Towards Fully Automated Open-Ended Scientific Discovery. arXiv preprint arXiv:2408.06292 (2024). arXiv:2408.06292 [cs.AI] https://arxiv.org/abs/2408.06292 [15] Aman Madaan, Niket Tandon, Prakhar Gupta, Skyler Hallinan, Luyu Gao, Sarah Wiegreffe, Uri Alon, Nouha Dziri, Shikha Prabhumoye, Yiming Yang, et al.2023. Self-Refine: Iterative Refinement with Self-Feedback. In Advances in Neural Infor- mation Processing Systems. [16]Kalyan Varma Nadimpalli, Shashank Reddy Chirra, Pradeep Varakantham, and Stefan Bauer. 2025. Evolving RL: Discovering New Activation Functions using LLMs. In Towards Agentic AI for Science: Hypothesis Generation, Comprehension, Quantification, and Validation. https://openreview.net/forum?id=H2x9juCuJg [17]Deepak Nathani, Lovish Madaan, Nicholas Roberts, Nikolay Bashlykov, Ajay Menon, Vincent Moens, Mikhail Plekhanov, Amar Budhiraja, Despoina Magka, Vladislav Vorotilov, et al.[n.d.]. MLGym: A New Framework and Benchmark for Advancing AI Research Agents. In Second Conference on Language Modeling. [18]Alexander Novikov, Ngân V ̃ u, Marvin Eisenberger, Emilien Dupont, Po-Sen Huang, Adam Zsolt Wagner, Sergey Shirobokov, Borislav Kozlovskii, Francisco JR Ruiz, Abbas Mehrabian, et al.2025. AlphaEvolve: A coding agent for scientific and algorithmic discovery. arXiv preprint arXiv:2506.13131 (2025). [19]Bernardino Romera-Paredes, Mohammadamin Barekatain, Alexander Novikov, Matej Balog, M Pawan Kumar, Emilien Dupont, Francisco JR Ruiz, Jordan S Ellenberg, Pengming Wang, Omar Fawzi, et al.2024. Mathematical discoveries from program search with large language models. Nature 625, 7995 (2024), 468–475. [20]Yiding Shi, Jianan Zhou, Wen Song, Jieyi Bi, Yaoxin Wu, and Jie Zhang. 2025. Generalizable Heuristic Generation Through Large Language Models with Meta- Optimization. arXiv preprint arXiv:2505.20881 (2025). [21]Yiwen Sun, Xianyin Zhang, Shiyu Huang, Shaowei Cai, BingZhen Zhang, and Ke Wei. 2024. AutoSAT: Automatically optimize sat solvers via large language models. arXiv preprint arXiv:2402.10705 (2024). [22]Anja Šurina, Amin Mansouri, Lars CPM Quaedvlieg, Amal Seddas, Maryna Via- zovska, Emmanuel Abbe, and Caglar Gulcehre. [n.d.]. Algorithm Discovery With LLMs: Evolutionary Search Meets Reinforcement Learning. In Second Conference on Language Modeling. [23]Nguyen Thach, Aida Riahifar, Nathan Huynh, and Hau Chan. 2025. RedAHD: Reduction-Based End-to-End Automatic Heuristic Design with Large Language Models. arXiv preprint arXiv:2505.20242 (2025). [24]Thibaut Vidal. 2022. Hybrid genetic search for the CVRP: Open-source implemen- tation and SWAP* neighborhood. Computers & Operations Research 140 (2022), 105643. [25] Xuezhi Wang, Jason Wei, Dale Schuurmans, Quoc Le, Ed Chi, Sharan Narang, Aakanksha Chowdhery, and Denny Zhou. 2023. Self-Consistency Improves Chain of Thought Reasoning in Language Models. In International Conference on Learning Representations. [26] Niels A Wouda, Leon Lan, and Wouter Kool. 2024. PyVRP: A high-performance VRP solver package. INFORMS Journal on Computing 36, 4 (2024), 943–955. [27] Xuan Wu, Di Wang, Chunguo Wu, Lijie Wen, Chunyan Miao, Yubin Xiao, and You Zhou. 2025. Efficient Heuristics Generation for Solving Combinatorial Optimization Problems Using Large Language Models. In SIGKDD Conference on Knowledge Discovery and Data Mining. [28]Yutaro Yamada, Robert Tjarko Lange, Cong Lu, Shengran Hu, Chris Lu, Jakob Foerster, Jeff Clune, and David Ha. 2025. The AI Scientist-v2: Workshop- Level Automated Scientific Discovery via Agentic Tree Search. arXiv preprint arXiv:2504.08066 (2025). arXiv:2504.08066 [cs.AI] https://arxiv.org/abs/2504.08066 [29]Xianliang Yang, Lei Song, Yapu Zhang, and Jiang Bian. 2025. HeurAgenix: A Multi-Agent LLM-Based Paradigm for Adaptive Heuristic Evolution and Selection in Combinatorial Optimization. https://openreview.net/forum?id=xxSK3ZNAhh [30]Shunyu Yao, Xi Lin, Jiashu Wang, Qingfu Zhang, and Zhenkun Wang. 2024. Rethinking supervised learning based neural combinatorial optimization for routing problem. ACM Transactions on Evolutionary Learning (2024). [31] Shunyu Yao, Fei Liu, Xi Lin, Zhichao Lu, Zhenkun Wang, and Qingfu Zhang. 2025. Multi-objective evolution of heuristic using large language model. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39. 27144–27152. [32]Haoran Ye, Jiarui Wang, Zhiguang Cao, Federico Berto, Chuanbo Hua, Haeyeon Kim, Jinkyoo Park, and Guojie Song. 2024. ReEvo: Large Language Models as Hyper-Heuristics with Reflective Evolution. In Advances in Neural Information Processing Systems. [33] Huigen Ye, Hua Xu, An Yan, and Yaoyang Cheng. 2025. Large Language Model- driven Large Neighborhood Search for Large-Scale MILP Problems. In Interna- tional Conference on Machine Learning. [34]Ni Zhang, Zhiguang Cao, Jianan Zhou, Cong Zhang, and Yew-Soon Ong. 2025. An Agentic Framework with LLMs for Solving Complex Vehicle Routing Problems. In International Conference on Learning Representations. [35]Zhi Zheng, Zhuoliang Xie, Zhenkun Wang, and Bryan Hooi. 2025. Monte Carlo Tree Search for Comprehensive Exploration in LLM-Based Automatic Heuristic Design. In International Conference on Machine Learning. [36]Jianan Zhou, Zhiguang Cao, Yaoxin Wu, Wen Song, Yining Ma, Jie Zhang, and Xu Chi. 2024. MVMoE: Multi-Task Vehicle Routing Solver with Mixture-of-Experts. In International Conference on Machine Learning. [37]Rongjie Zhu, Cong Zhang, and Zhiguang Cao. 2025. Refining Hybrid Genetic Search for CVRP via Reinforcement Learning-Finetuned LLM. In International Conference on Learning Representations. A VRP ATTRIBUTE TAXONOMY To handle the vast landscape of VRP variants in a structured manner, we describe VRP instances via composable attributes. This allows us to represent any variant as a combination of fundamental prop- erties, enabling systematic generation and evaluation of evolved components across a multitude of problem types. The attributes are categorized into node-specific and global properties, each with precise mathematical formulations. Let퐺= (푉,퐴)be a directed graph where푉= 0,1,2, . . .,푛is the set of nodes with node 0 representing the depot and1,2, . . .,푛representing customers, and퐴=(푖, 푗):푖, 푗 ∈ 푉,푖≠ 푗is the set of arcs. Let푚denote the number of available vehicles. Node attributes: Properties specific to each customer or depot node푖 ∈ 푉. These attributes define the local requirements and characteristics that must be satisfied during service. •Demand: The demand attribute represents the quantity of goods to be delivered to or picked up from each customer. Let 푞 푖 ∈R denote the demand at node 푖, where: 푞 푖 = > 0 if customer 푖 requires delivery (linehaul) < 0 if customer 푖 requires pickup (backhaul) = 0if node 푖 is the depot (2) • Time Window (TW): Each customer푖has a designated time window[푎 푖 ,푏 푖 ]within which service must commence. Let푡 푖 denote the service start time at customer푖. The time window constraints are: 푎 푖 ≤ 푡 푖 ≤ 푏 푖 , ∀푖 ∈ 푉 \0(3) 푡 푖 +푠 푖 +휏 푖푗 ≤ 푡 푗 + 푀(1− 푥 푖푗푘 ),(4) ∀푖, 푗 ∈ 푉,푘 ∈ 1, . . .,푚(5) where푠 푖 ≥0 is the service time at customer푖,휏 푖푗 ≥0 is the travel time from푖to푗,푀is a sufficiently large constant, and 푥 푖푗푘 ∈ 0, 1 indicates whether vehicle 푘 traverses arc(푖, 푗). The depot time window[푎 0 ,푏 0 ]constrains the total route duration: 푎 0 ≤ 푡 start 0 ≤ 푡 end 0 ≤ 푏 0 (6) where푡 start 0 and푡 end 0 represent the departure and return times at the depot, respectively. •Service Time: Each customer푖requires a non-negative service time푠 푖 ≥0 during which the vehicle is occupied. This affects the temporal feasibility of routes and is incorporated into time window constraints through equation (5). The total service time for vehicle 푘 is: 푆 푘 = ∑︁ 푖∈푉 푠 푖 푦 푖푘 , ∀푘 ∈ 1, . . .,푚(7) where푦 푖푘 ∈ 0,1indicates whether customer푖is served by vehicle 푘 . Global attributes: Constraints or properties that apply to entire routes or the complete problem instance, affecting the optimization objective and feasibility conditions. • Capacity (C): Each vehicle푘has a finite capacity푄 푘 . The capacity constraint ensures that the total demand of all cus- tomers served by a single route cannot exceed the vehicle’s capacity: ∑︁ 푖∈푉 푞 푖 푦 푖푘 ≤ 푄 푘 , ∀푘 ∈ 1, . . .,푚(8) •Duration Limit (L): This constraint imposes a maximum allowable cost or duration퐷 푘 for each vehicle route푘. The constraint is formulated as: ∑︁ (푖,푗)∈퐴 푐 푖푗 푥 푖푗푘 + ∑︁ 푖∈푉 푠 푖 푦 푖푘 ≤ 퐷 푘 ,(9) ∀푘 ∈ 1, . . .,푚(10) where 푐 푖푗 represents the travel cost between nodes 푖 and 푗 . For homogeneous duration limits across all vehicles,퐷 푘 = 퐷 for all 푘 . • Backhaul (B): In backhaul problems, customers are parti- tioned into linehaul customers퐿= 푖 ∈ 푉:푞 푖 >0and backhaul customers퐵=푖 ∈ 푉:푞 푖 <0. The precedence constraint ensures all linehauls are served before backhauls: ∑︁ 푗∈퐵 ∑︁ (푖,푗)∈퐴 푥 푖푗푘 = 0,(11) ∀푖 ∈ 퐿,푘 ∈ 1, . . .,푚(12) Alternatively, using auxiliary variables푢 푖푘 ∈R + represent- ing the position of customer 푖 in vehicle 푘 ’s route: 푢 푖푘 + 1≤ 푢 푗푘 +푛(1− 푥 푖푗푘 ),(13) ∀푖 ∈ 퐿, 푗 ∈ 퐵,푘 ∈ 1, . . .,푚(14) where 푛=|푉| is the total number of nodes. •Open Route (O): In open route problems, vehicles are not re- quired to return to the depot after serving their last customer. This modifies the route structure by eliminating mandatory return arcs. The flow conservation constraints become: ∑︁ 푗∈푉\0 푥 0푗푘 = 푧 푘 , ∀푘 ∈ 1, . . .,푚(15) ∑︁ 푗∈푉 푥 푗푖푘 − ∑︁ 푗∈푉 푥 푖푗푘 = 0,(16) ∀푖 ∈ 푉 \0,푘 ∈ 1, . . .,푚(17) ∑︁ 푗∈푉\0 푥 푗0푘 ≤ 푧 푘 , ∀푘 ∈ 1, . . .,푚(18) where 푧 푘 ∈ 0, 1 indicates whether vehicle 푘 is used. • Multi-Depot (MD): Let퐷=푑 1 ,푑 2 , . . .,푑 푝 ⊂ 푉denote the set of depots. Each vehicle푘is assigned to exactly one depot depot(푘) ∈ 퐷 : ∑︁ 푑∈퐷 푤 푘푑 = 1, ∀푘 ∈ 1, . . .,푚(19) where푤 푘푑 ∈ 0,1indicates whether vehicle푘is assigned to depot 푑 . The route must start and end at the assigned depot: ∑︁ 푗∈푉\퐷 푥 푑푗푘 =푤 푘푑 , ∀푑 ∈ 퐷,푘 ∈ 1, . . .,푚(20) ∑︁ 푖∈푉\퐷 푥 푖푑푘 =푤 푘푑 ,(21) ∀푑 ∈ 퐷,푘 ∈ 1, . . .,푚 (if not open routes)(22) Each customer must be served by exactly one vehicle: 푚 ∑︁ 푘=1 푦 푖푘 = 1, ∀푖 ∈ 푉 \ 퐷(23) •Mixed Backhaul (MB): Mixed backhaul relaxes the strict precedence constraint of standard backhaul problems, al- lowing pickups and deliveries to be intermixed on the same route. Letℓ 푖푘 ≥0 represent the load of vehicle푘after serving customer 푖. The constraints are: 0≤ ℓ 푖푘 ≤ 푄 푘 , ∀푖 ∈ 푉,푘 ∈ 1, . . .,푚(24) ℓ 푗푘 = ℓ 푖푘 +푞 푗 − 푀(1− 푥 푖푗푘 ),(25) ∀(푖, 푗) ∈ 퐴,푘 ∈ 1, . . .,푚(26) ℓ 0푘 = 0, ∀푘 ∈ 1, . . .,푚(27) Using this unified format, any VRP variant can be systematically represented as: VRP Variant= Base VRP∪A node ∪A global (28) whereA node ⊆ TW, Service TimeandA global ⊆ C, L, B, O, MD, MB. B FULL LLM PROMPT The complete prompt template used in our experiments for parents selection is provided below. This template dynamically incorporates the two existing parent selectors ($code1and$code2) along with their performance metrics to guide the generation of improved solutions. We’l be providing more prompt templates alongside the release of the code. The prompt template for our experiments is structured into the following sequential sections: (1) Role and Planning (2) Context and Existing Approaches (3) Reasoning Phase and Task Description (4) Implementation Instructions (5) Implementation Requirements (6) Reference Classes (7) Reflection Phase and Response Format (8) Existing Parent Selectors Complete Prompt Template - Role and Planning Role: AI Optimization Researcher (Vehicle Routing Problem) 1) Planning Phase Before you inspect any code, pause and map out your thinking. (1)Pitfalls (퐾 푝 ): List three ways a parent-selector can go wrong in VRP (e.g., over-exploitation, diversity collapse, cost-feasibility trade-offs). (2)Strategies (퐾 푠 ): For each pitfall, propose one high-level counter- measure (e.g., partial elitism + random injection). (3)Hidden Traps (퐾 푡 ): What tricky instance features (e.g., clustered deliveries, time-windows) could mislead selection? Overall Objective: Only write a Python functionselect_parentsimple- menting a novel parent selection strategy for a genetic algorithm (Hybrid Genetic Search) solving the Vehicle Routing Problem (VRP). Unlike standard strategies, your goal is to propose a new method—do not replicate existing implementations exactly. The function will be used within the PyVRP li- brary. Use your planning insights to guide design. Do not redefine any of the provided classes; just import them if needed. Lower cost is better. Context and Existing Approaches Context: In Hybrid Genetic Search for VRP, parents are selected from a population of solutions to undergo crossover. Parent selection directly influences exploration (diversity) and exploitation (quality) in the search. Standard selection strategies include: • Fitness-based selection • Tournament selection • Random selection with diversity bias Common Existing Approaches: • Elitist Selection: Always select top-ranked (lowest-cost) individ- uals. • Tournament Selection: Randomly sample a subset of the popu- lation and pick the best. •Fitness Proportionate Selection: Choose solutions with proba- bility proportional to inverse cost. •Diversity-Aware Selection: Combine one elite solution with one structurally different (e.g., high Hamming distance) solution. •Random Selection: Uniform random sampling from feasible in- dividuals. Reasoning Phase and Task Description 2) Reasoning Phase Now dig into sel1 and sel2—but think aloud. •Identify Patterns: What exploration/exploitation bias does each show? • Assess Shortcomings: Where does each stumble relative to your 퐾 푝 /퐾 푡 list? •Design Sketch: Using your퐾 푠 strategies, outline in prose how sel3 will avoid those failures. Your Task in This Interaction: You will be presented with two existing implementations of parent selection logic which are in Python,sel1and sel2, along with observed performance metrics. Your goal is to propose a new Python implementation,sel3, that performs better (produces better children in crossover) and less cost compared to the given ones. Implementation Instructions Instructions: (1) Analyze sel1 and sel2: • Understand how it balances exploration and exploitation. • Evaluate whether it favors elite, random, or diverse solu- tions—and if it filters infeasible ones. •Identify reasons why it may not generalize or may produce low-quality parents. (2)Designsel3 = select_parents(population, rng, cost_evaluator, k=2): • Propose a novel selection mechanism that improves the ex- pected quality and diversity of selected parents. • Consider hybrid strategies (e.g., partial elitism, feasibility filters, diversity-aware tournaments). • Avoid trivial combinations of existing methods. Implementation Requirements (3) Implementation Requirements: • Input: – population: list[Solution] – rng: RandomNumberGenerator – cost_evaluator: CostEvaluator – k: int (number of parents, always 2) • Output: tuple[Solution, Solution] • Function Signature: Must be exactly: 1 def select_parents( 2 population: list[Solution], 3 rng: RandomNumberGenerator , 4 cost_evaluator: CostEvaluator , 5 k: int = 2 6 ) -> tuple[Solution , Solution ]: • Behavior: –Prioritize feasible solutions but allow occasional infeasible ones if beneficial. –Use cost_evaluator.penalised_cost(solution) to assess solu- tion quality. – Ensure diversity between selected parents. – Avoid modifying the population in-place. • Performance Hints: – Cache penalized costs if reused multiple times. – Minimize sorting or scanning the entire population. – Userng.randint(len(population))to sample indices ef- ficiently. • Constraints: – Use only Python built-ins and PyVRP classes. Existing Parent Selectors Pair of Existing Parent Selectors: Selector 1: 1 $code1 Current Cost: $score1 Cost to beat (Target Cost/Current best cost): $baseline_score Feedback: $feedback1 Selector 2: 1 $code2 Current Cost: $score2 Cost to beat (Target Cost/Current best cost): $baseline_score Feedback: $feedback2 Reference Classes Reference Python Classes Implementation: The following class definitions are already available in another file and are provided here for reference only. Please do not redefine them. They can be imported using:from pyvrp._pyvrp import CostEvaluator, RandomNumberGenerator, Solution, etc. 1 class CostEvaluator: 2 def __init__( 3 self , 4 load_penalties: list[float], 5 tw_penalty: float , 6 dist_penalty: float , 7 ) -> None: ... 8 def load_penalty( 9 self , load: int , capacity: int , dimension: int 10 ) -> int: ... 11 def tw_penalty(self , time_warp: int) -> int: ... 12 def dist_penalty(self , distance: int , max_distance: int) -> int: ... 13 def penalised_cost(self , solution: Solution) -> int: ... 14 def cost(self , solution: Solution) -> int: ... 15 16 class Solution: 17 def __init__( 18 self , 19 data: ProblemData , 20 routes: list[Route] | list[list[int]], 21 ) -> None: ... 22 @classmethod 23 def make_random( 24 cls , data: ProblemData , rng: RandomNumberGenerator 25 ) -> Solution: ... 26 def is_feasible(self) -> bool: ... 27 # ... (additional methods omitted for brevity) 28 29 class RandomNumberGenerator: 30 @overload 31 def __init__(self , seed: int) -> None: ... 32 @overload 33 def __init__(self , state: list[int]) -> None: ... 34 def rand(self) -> float: ... 35 def randint(self , high: int) -> int: ... 36 # ... (additional methods omitted for brevity) 37 38 # ... (other PyVRP class definitions omitted for brevity) Reflection Phase and Response Format 3) Reflection Phase After you draft sel3, self-evaluate before returning code. (1)Which of your original퐾 푝 pitfalls does the final code best address? Which remain? (2) Did any new traps surface during coding? How would you tweak sel3 to handle them? (3)Propose two concrete revisions to improve quality/diversity trade- off further. Response Format: Enclose your answer in a single Python code block (python) with only the functionselect_parentsand necessary imports. The code must be directly extractable. Response Format Code: 1 from pyvrp._pyvrp import CostEvaluator , RandomNumberGenerator , Solution 2 3 def select_parents( 4 population: list[Solution], 5 rng: RandomNumberGenerator , 6 cost_evaluator: CostEvaluator , 7 k: int = 2 8 ) -> tuple[Solution , Solution ]: 9 """ 10 [Brief description of your novel parent selection strategy] 11 """ 12 # [Your implementation here] 13 parent1 = ... 14 parent2 = ... 15 return parent1 , parent2 C PARENT SELECTION CODE This section presents the complete implementation, from Figures 4 to 7, of the parent selection algorithm that was discovered. The code is divided to ensure related functionality stays together while keeping each chunk at a reasonable size for page layout. D SUPPLEMENTARY DISCUSSION D.1 Comparison with Self-Refine and Self-Consistency While approaches such as Self-Refine [15] and Self-Consistency [25] have demonstrated the value of iterative LLM self-improvement, they operate at fundamentally different levels of abstraction com- pared to MEP. Self-Refine enables an LLM to critique and polish a generated output through successive rounds of feedback, function- ing as an output-level refinement loop. Self-Consistency samples multiple reasoning paths and selects the most frequent answer, im- proving reliability through redundancy. Neither approach enforces the explicit diagnostic reasoning that MEP mandates. MEP draws specific inspiration from the metacognitive frame- work of Bai et al. [1], which models human-like planning and rea- soning before acting. The critical distinction lies in the structured Reason-Act-Reflect cycle: before generating any code, the LLM must (1) explicitly diagnose the weaknesses of parent heuristics against known failure modes (퐾 푝 ,퐾 푡 ), (2) formulate a concise, testable de- sign hypothesis, and (3) reflect on its own implementation post-hoc. This proactive, hypothesis-driven process is qualitatively different 1 from pyvrp._pyvrp import CostEvaluator , RandomNumberGenerator , Solution 2 3 def select_parents(population: list[Solution], rng: RandomNumberGenerator , 4 cost_evaluator: CostEvaluator , k: int =2) -> tuple[Solution , Solution ]: 5 """ 6 Hybrid Adaptive Stratified Diversity -Pressure Parent Selector 7 ----------------------------------------------------------------------------------- 8 9 Summary: 10 This selector provides scalable parent selection that maintains structural diversity awareness 11 for manageable populations while gracefully degrading to efficient cost -based selection 12 for large -scale problems. It employs population -size adaptive strategies to balance quality 13 and computational efficiency , preventing timeouts in large populations while preserving 14 diversity -driven selection for smaller populations. 15 16 Step -by -step: 17 --------------- 18 1. ** Precompute Costs & Stratification **: 19 - Cache penalized costs and feasibility flags for all solutions upfront. 20 - Sort population by cost and partition into three strata: 21 * Elite: top 1/6 solutions (highest quality) 22 * Mid: solutions from 1/6 to 2/3 (moderate quality) 23 * Tail: bottom 1/3 solutions (lower quality) 24 25 2. ** Parent 1 Selection **: 26 - Select from elite stratum using feasibility -cost biased tournament. 27 - Tournament emphasizes best feasible solutions with occasional infeasible survivors 28 (20% probability) for diversity maintenance. 29 30 3. ** Adaptive Parent 2 Strategy **: 31 - Determine complementary strata based on Parent 1 characteristics: 32 * If Parent 1 is elite & feasible : select from mid+tail (exploration) 33 * If Parent 1 is infeasible or non -elite : select from elite+mid (intensification) 34 35 4. ** Population -Size Adaptive Selection **: 36 - ** Large populations ( >500) **: Skip structural analysis entirely , use simplified scoring: 37 * Cost advantage (60% weight) 38 * Feasibility bonus (30% weight) 39 * Basic diversity bonus (10% weight) 40 - ** Small populations (= <500) **: Full structural analysis with lazy encoding: 41 * Structural dissimilarity via client -set overlap (55% weight) 42 * Cost advantage (30% weight) 43 * Feasibility bonus (15% weight) 44 45 5. **Lazy Encoding & Sampling **: 46 - Only encode solutions that are actually evaluated as candidates. 47 - Adaptive sampling limits: 20-30 candidates per stratum depending on population size. 48 - Add small randomness to scores to avoid local minima. 49 50 ----------------------------------------------------------------------------------- 51 This approach explicitly counters: 52 - Computational timeouts in large populations through adaptive complexity reduction. 53 - Over -exploitation by maintaining structural diversity awareness when computationally feasible. 54 - Population stratification collapse through complementary strata selection strategy. 55 - Unnecessary computational overhead by employing lazy encoding and adaptive sampling. 56 - Pure cost or pure diversity bias through weighted multi -criteria scoring. 57 58 Implementation is efficient: 59 Costs and feasibility are cached once; structural encoding is lazy and only applied 60 when needed; population -size thresholds automatically switch between full analysis 61 and simplified selection; adaptive sampling prevents excessive candidate evaluation. 62 """ 63 64 n = len(population) 65 assert n >= 2,'Population must have at least two solutions.' 66 assert k == 2,'Only k=2 supported.' 67 68 # Cache costs and feasibility upfront for efficient repeated access 69 costs = [cost_evaluator.penalised_cost(sol) for sol in population] 70 feasible = [sol.is_feasible () for sol in population] 71 sorted_indices = sorted(range(n), key=lambda i: costs[i]) 72 73 # Partition population into three strata by cost quantiles 74 # Elite: top 1/6, Mid: 1/6 to 2/3, Tail: bottom 1/3 75 elite_cut = max(1, n // 6) 76 mid_cut = max(elite_cut + 1, n * 2 // 3) 77 elite_indices = sorted_indices [: elite_cut] 78 mid_indices = sorted_indices[elite_cut:mid_cut] 79 tail_indices = sorted_indices[mid_cut :] Figure 4: Hybrid Adaptive Stratified Diversity-Pressure Parent Selector. 1 def encode_solution_simple(solution: Solution): 2 """ 3 Lightweight client -set based structural encoding for performance. 4 Captures global client distribution and per -route client assignments 5 without expensive edge -pair analysis from original version. 6 """ 7 all_clients = set() 8 route_clients = [] 9 for route in solution.routes (): 10 clients_in_route = set(v for v in route if v > 0) 11 route_clients.append(clients_in_route) 12 all_clients.update(clients_in_route) 13 return (all_clients , route_clients) 14 15 def simple_structural_distance(enc1 , enc2): 16 """ 17 Simplified structural dissimilarity based on client overlap patterns. 18 Combines global client set similarity with route -level client distribution , 19 replacing the complex harmonic mean of node/edge Jaccard distances. 20 """ 21 all_clients1 , route_clients1 = enc1 22 all_clients2 , route_clients2 = enc2 23 24 # Global client set similarity (replacement for node -set similarity) 25 if not all_clients1 and not all_clients2: 26 return 0.0 27 28 intersection = len(all_clients1 & all_clients2) 29 union = len(all_clients1 | all_clients2) 30 31 if union == 0: 32 return 0.0 33 34 global_similarity = intersection / union 35 36 # Route structure similarity (simplified replacement for edge -set analysis) 37 min_routes = min(len(route_clients1), len(route_clients2)) 38 if min_routes == 0: 39 return 1.0 - global_similarity 40 41 route_overlaps = [] 42 for i in range(min_routes): 43 if not route_clients1[i] and not route_clients2[i]: 44 continue 45 r_intersection = len(route_clients1[i] & route_clients2[i]) 46 r_union = len(route_clients1[i] | route_clients2[i]) 47 if r_union > 0: 48 route_overlaps.append(r_intersection / r_union) 49 50 if route_overlaps: 51 avg_route_similarity = sum(route_overlaps) / len(route_overlaps) 52 else: 53 avg_route_similarity = 0.0 54 55 # Combine global and route -level similarity (70% global , 30% route -level weighting) 56 combined_similarity = 0.7 * global_similarity + 0.3 * avg_route_similarity 57 return 1.0 - combined_similarity Figure 5: Hybrid Adaptive Stratified Diversity-Pressure Parent Selector. from reactive refinement. In Self-Refine, the model generates first and critiques second; in MEP, diagnosis and hypothesis formulation necessarily precede implementation. This ordering is critical for complex algorithmic design tasks, where generating code without a clear strategic intent leads to superficial mutations rather than prin- cipled innovations. In the classical hyper-heuristic taxonomy, MEP thus functions as a metacognitive construction hyper-heuristic, operating in the space of heuristic generation rather than heuris- tic selection, with the LLM serving as both the generator and the strategic reasoner. D.2 Illustrative Failure-Reflection Case To concretize the Reason-Act-Reflect cycle, we describe a represen- tative case observed during the evolution of theselect_parents operator. In generation 3, the LLM was presented with two parent heuris- tics: one that selected parents purely by cost rank (achieving low average cost but exhibiting premature convergence on clustered instances), and another that used uniform random selection (main- taining diversity but producing poor-quality offspring). Reason. The LLM diagnosed that the first parent suffered from diversity collapse—a pitfall identified in 퐾 푝 —noting that “selecting exclusively from the top cost stratum causes the population to converge to a narrow region of the solution space, particularly on instances with clustered customer distributions.” It further identified that the second parent’s lack of quality pressure wasted computational budget on unpromising crossover pairs. Act. The LLM formulated the design hypothesis: “A stratified selection mechanism that draws the first parent from the elite stra- tum and the second parent from complementary strata, weighted by structural dissimilarity, will balance exploitation of high-quality solutions with exploration of diverse genetic material.” This led to the 1 def tournament_select(indices , t_size=7, allow_infeasible_prob =0.2): 2 """ 3 Combined feasibility -cost biased tournament emphasizing best feasible solutions 4 with occasional infeasible survivors for diversity maintenance. 5 """ 6 if len(indices) <= t_size: 7 candidates = indices [:] 8 else: 9 candidates = [indices[rng.randint(len(indices))] for _ in range(t_size)] 10 11 # Filter candidates by feasibility with probabilistic infeasible inclusion 12 filtered = [] 13 for idx in candidates: 14 if feasible[idx]: 15 filtered.append(idx) 16 elif rng.rand() < allow_infeasible_prob: 17 filtered.append(idx) 18 19 if not filtered: 20 filtered = candidates 21 22 best_idx = min(filtered , key=lambda i: costs[i]) 23 return best_idx 24 25 # Select parent1 from elite stratum using feasibility -cost biased tournament 26 parent1_idx = tournament_select(elite_indices) 27 parent1 = population[parent1_idx] 28 parent1_feas = feasible[parent1_idx] 29 p1_cost_rank = sorted_indices.index(parent1_idx) 30 31 # Determine complementary parent2 strata based on parent1 characteristics 32 if parent1_feas and p1_cost_rank < elite_cut: 33 parent2_strata = mid_indices + tail_indices 34 else: 35 parent2_strata = elite_indices + mid_indices 36 37 if not parent2_strata: 38 parent2_strata = sorted_indices 39 40 # Adaptive sampling limits to balance quality vs. performance 41 # Reduce sample size for large populations to prevent timeout 42 max_sample_size = min(20, len(parent2_strata)) if n > 500 else min(30, len(parent2_strata)) 43 44 # Sample parent2 candidates with collision avoidance 45 sampled_parent2_indices = [] 46 seen = set() 47 attempts = 0 48 while len(sampled_parent2_indices) < max_sample_size and attempts < max_sample_size * 2: 49 candidate = parent2_strata[rng.randint(len(parent2_strata))] 50 if candidate != parent1_idx and candidate not in seen: 51 sampled_parent2_indices.append(candidate) 52 seen.add(candidate) 53 attempts += 1 54 55 # Fallback candidate selection if sampling fails 56 if not sampled_parent2_indices: 57 for idx in parent2_strata: 58 if idx != parent1_idx: 59 sampled_parent2_indices.append(idx) 60 break 61 else: 62 # Last resort fallback 63 sampled_parent2_indices.append(parent2_strata [0]) Figure 6: Hybrid Adaptive Stratified Diversity-Pressure Parent Selector. implementation of cost-based stratification with adaptive diversity pressure. Reflect. The LLM self-critiqued that its implementation addressed the diversity collapse pitfall (퐾 푝 ) but noted a potential vulnerability: “On very small populations (푛<10), the stratification boundaries may yield empty strata, requiring fallback logic.” This reflection was car- ried forward and informed the next generation, which introduced explicit handling of degenerate population sizes. This case illustrates how the structured cycle produces targeted improvements grounded in diagnosed weaknesses, rather than ar- bitrary code mutations. D.3 Comparison with Neural Solvers To contextualize the performance of MEP-evolved heuristics rel- ative to neural combinatorial optimization approaches, we eval- uated POMO-MTL, a representative deep reinforcement learning (DRL)-based solver, on the standard CVRP test set. The results are presented in Table 5. While MEP improves upon the SOTA HGS baseline, the DRL- based solver lags behind by a substantial margin of 11.06%. This comparison underscores a key motivation of our work: rather than replacing traditional solvers with end-to-end neural policies—which currently cannot match the performance of well-engineered meta- heuristics on realistic VRP instances—MEP enhances the existing 1 # Population -size adaptive parent2 selection strategy 2 if n > 500: 3 # Large population: Simplified selection focusing on cost and feasibility 4 # Skip structural analysis entirely to prevent timeouts 5 best_score = -float('inf') 6 parent2_idx = None 7 8 p2_costs = [costs[i] for i in sampled_parent2_indices] 9 p2_min_cost = min(p2_costs) 10 p2_max_cost = max(p2_costs) if max(p2_costs) > p2_min_cost else p2_min_cost + 1 11 12 for idx in sampled_parent2_indices: 13 if idx == parent1_idx: 14 continue 15 16 cand_cost = costs[idx] 17 cand_feas = feasible[idx] 18 19 # Simplified three -component scoring without structural analysis 20 cost_score = 1.0 - (cand_cost - p2_min_cost) / (p2_max_cost - p2_min_cost) 21 feas_bonus = 1.0 if cand_feas else 0.0 22 diversity_bonus = 0.5 if abs(costs[idx] - costs[parent1_idx ]) > (p2_max_cost - p2_min_cost) * 0.1 else 0.0 23 24 # Weighted combination emphasizing cost with feasibility and basic diversity 25 score = 0.6 * cost_score + 0.3 * feas_bonus + 0.1 * diversity_bonus 26 score += (rng.rand() - 0.5) * 0.02 # Small randomness to avoid local minima 27 28 if score > best_score: 29 best_score = score 30 parent2_idx = idx 31 32 else: 33 # Small population: Full structural analysis with lazy encoding 34 # Only encode parent1 and actual candidates (not entire population) 35 p1_enc = encode_solution_simple(population[parent1_idx ]) 36 37 p2_costs = [costs[i] for i in sampled_parent2_indices] 38 p2_min_cost = min(p2_costs) 39 p2_max_cost = max(p2_costs) if max(p2_costs) > p2_min_cost else p2_min_cost + 1 40 41 # Weighted combination of three normalized signals 42 w_struct = 0.55 # Structural dissimilarity weight 43 w_cost = 0.3 # Cost advantage weight 44 w_feas = 0.15 # Feasibility bonus weight 45 best_score = -float('inf') 46 parent2_idx = None 47 48 for idx in sampled_parent2_indices: 49 if idx == parent1_idx: 50 continue 51 52 cand_cost = costs[idx] 53 cand_feas = feasible[idx] 54 cand_enc = encode_solution_simple(population[idx]) # Lazy encoding 55 56 # Three -component scoring with structural dissimilarity 57 struct_dist = simple_structural_distance(p1_enc , cand_enc) 58 cost_score = 1.0 - (cand_cost - p2_min_cost) / (p2_max_cost - p2_min_cost) 59 feas_bonus = 1.0 if cand_feas else 0.0 60 61 # Weighted combination avoiding pure cost or pure diversity bias 62 score = w_struct * struct_dist + w_cost * cost_score + w_feas * feas_bonus 63 score += (rng.rand() - 0.5) * 0.02 # Small randomness to avoid local minima 64 65 if score > best_score: 66 best_score = score 67 parent2_idx = idx 68 69 # Final fallback if no suitable parent2 found 70 if parent2_idx is None: 71 parent2_idx = tournament_select(parent2_strata) 72 73 parent2 = population[parent2_idx] 74 return (parent1 , parent2) Figure 7: Hybrid Adaptive Stratified Diversity-Pressure Parent Selector. SOTA solver from within, leveraging the LLM’s reasoning capabil- ities to discover better algorithmic components while preserving the proven search framework of HGS. D.4 Detailed Comparison with ReEvo and EoH The “Reactive Evolution” baseline in Table 3 is explicitly designed to simulate the core mechanics of prior feedback-driven LLM-based hyper-heuristics, including ReEvo [32] and Evolution of Heuristics (EoH) [13]. In this variant, the prompt provides only the parent heuristics and their performance scores, asking the LLM to gen- erate an improved version based on direct feedback, without any structured diagnosis, hypothesis formulation, or self-reflection. Table 5: Performance comparison of MEP-Evolved, PyVRP Baseline, and POMO-MTL (DRL Solver) on the CVRP benchmark. Lower cost indicates better performance. MethodAverage Cost Gap vs. HGS Baseline MEP-Evolved (Ours)64 236.30 −0.33 % (Improvement) PyVRP Baseline (HGS)64 450.030.00 % POMO-MTL (DRL Solver)71 575.70+11.06 % (Worse) The results reveal a stark contrast. On VRPTW, Reactive Evo- lution achieves only a 0.10% improvement, while the full MEP framework achieves 2.25% under identical conditions—a 22×differ- ence in effectiveness. More critically, Reactive Evolution produces performance degradation on three of six variants (−0.20% on CVRP, −0.39% on VRPB,−2.05% on MDVRPTW), indicating that feedback- driven mutation without structured reasoning can actively harm the solver when applied to complex algorithmic components. In contrast, the full MEP framework achieves positive improvements on all six variants without exception. These results suggest that while the reactive paradigm of ReEvo and EoH is effective for evolving simpler constructive heuristics for basic problem variants, it is insufficient for the more challeng- ing task of improving components within a sophisticated, well- optimized SOTA solver. The metacognitive scaffolding of MEP— specifically the mandate to diagnose failures and formulate hypothe- ses before code generation—is the critical ingredient that enables effective heuristic discovery at this level of algorithmic complexity. D.5 Limitations We acknowledge several limitations of the current work. Proprietary Model Dependence. Our framework relies on GPT-4.1 for heuristic generation. While our ablation study (Table 3) demon- strates that it is the MEP methodology—not the specific LLM—that drives performance gains (as evidenced by the clear hierarchy be- tween Full MEP, MEP-noInit, and Reactive Evolution under the same model), the dependence on a proprietary API limits repro- ducibility and accessibility. In preliminary experiments, smaller open-source models (e.g., LLaMA-70B) struggled with the complex reasoning and strict syntactic requirements needed for valid Python code generation in this context. As open-source models continue to mature in code generation and structured reasoning capabilities, adapting MEP to leverage these models is a priority for future work. Evolutionary Budget. Budget constraints limited the evolutionary process to푁=10 generations. While Figure 3 demonstrates rapid convergence within this limit, a more extensive search with larger generation budgets could potentially discover even more effective heuristics. The current results nonetheless provide strong evidence that MEP efficiently locates high-quality regions of the search space with a manageable one-time investment. Generalization Scope. While the evolved heuristics generalize well across the six VRP variants evaluated—despite being evolved ex- clusively on TSP100 instances—further validation across a broader set of problem distributions, larger instance scales, and the full set of 48 VRP variants supported by PyVRP is needed to fully characterize the boundaries of this generalization.