arXiv CSopen access

RAISE: LLM-based Automated Heuristic Design with Robust Adversary Instance Search

RAISE: LLM-based Automated Heuristic Design with Robust Adversary Instance Search

arXiv:2606.31801v1 [cs.AI] 30 Jun 2026

Fei Liu1,2

Alessio Figalli2 Patrick Owen1 Nicola Serra1 1 University of Zurich 2 ETH Zurich

Abstract Automated Heuristic Design (AHD) with Large Language Models (LLMs) has shown remarkable progress in discovering high-quality heuristics. However, existing LLM-based AHD methods optimize heuristics for a fixed training instance set and may fail catastrophically when deployed under real-world distributional shifts. We propose Robust Adversary Instance Search (RAISE), a framework that integrates constrained worst-case instance search within a principled neighborhood of the training distribution into the LLM-based evolutionary search loop. RAISE treats robust AHD as a constrained adversarial instance search problem: the outer loop evolves heuristics via LLM operators, while an LLM-free inner loop efficiently identifies hard instances within an ε-ball around the training instance set using a basis distribution parameterization with boundary projection. Comprehensive experiments on Online Bin Packing (OBP), Online Job Shop Scheduling (OJSP), and Online Vehicle Routing (OVRP) across five distribution families demonstrate that existing LLM-based AHD methods degrade by up to 19× under distribution shift, while RAISE consistently maintains strong performance across all tested distributions and problem scales.

1

Introduction

Designing effective heuristic algorithms for optimization problems has traditionally required extensive domain expertise and manual engineering. Automated Heuristic Design (AHD) seeks to reduce this burden by automatically discovering high-quality heuristic programs, enabling rapid adaptation to new problem domains without much human intervention [Pillay and Qu, 2018, Stützle and López-Ibáñez, 2018, Qu et al., 2020, Zhou et al., 2024, Ma et al., 2025]. The recent integration of Large Language Models (LLMs) into evolutionary search has dramatically advanced AHD [Liu et al., 2024, 2023, Van Stein and Bäck, 2024, Ye et al., 2024, Dat et al., 2025, Zheng et al., 2025, Hu and Zhang, 2025, Ma et al., 2026, Shi et al., 2026, Liu et al., 2026, Novikov et al., 2025, Sun et al., 2026]. For example, Liu et al. [2024] propose Evolution of Heuristics (EoH), which uses LLMs as evolutionary operators to iteratively generate, recombine, and refine heuristic functions for combinatorial optimization. Ye et al. [2024] extend this paradigm with ReEvo, incorporating reflective short-term and long-term memory mechanisms that guide the LLM toward progressively better heuristics. Additionally, MCTS-AHD Zheng et al. [2025] and PartEvo Hu and Zhang [2025] integrate Monte Carlo tree search and feature-assisted niche construction, respectively, to enhance the efficiency of evolutionary search. The distribution shift problem. Despite their strong nominal performance (performance on training instance set), existing LLM-based AHD methods share a critical limitation: they optimize heuristics for a fixed instance distribution encountered during training, producing highly specialized solutions that may fail catastrophically when deployed under different conditions [Wang et al., 2026, Shi et al., 2026]. This robustness gap is particularly problematic in real-world deployment, where problem Preprint.

instances are drawn from unknown or shifting distributions [Zhou et al., 2024, Goh et al., 2025, Wang et al., 2026]. Recent work has attempted to address this limitation. Liu et al. [2026] propose EoH-S, which learns a portfolio of diverse heuristics to improve cross-distribution performance. However, it relies on a predefined set of diverse training instances. Shi et al. [2026] introduce MoH, a meta-optimization method for cross-distribution generalization, but still requires explicit multi-task training across manually-defined distributions. In both cases, the dependence on a predefined set of distributions limits practicality in dynamic or unpredictable environments. Furthermore, while Karimi et al. [2025] and Duan et al. [2025] have studied adversarial instance generation, they do not impose constraints on distributions and rely on LLMs to generate instances, which lacks controllability and introduces additional cost. This work. We propose RAISE (Robust Adversary Instance Search for LLM-based Heuristic Design), which incorporates adversarial instance search into the LLM evolutionary search loop. RAISE formulates robust AHD as a constrained minimax optimization problem to find the best robust heuristic h∗ : h∗ = arg max ′ min eval(h, s′ ), h

s ∈Bε (S)

where Bε (S) denotes the ball of radius ε around the nominal instance set S. The outer maximization is performed by the LLM evolutionary search, while the inner minimization is handled by an adversarial instance search that discovers the hardest instances within the uncertainty set. Contributions. This paper makes the following contributions: • Problem formulation. We formalize an instance-level robust AHD objective as a constrained minimax problem over an ε-ball uncertainty set, providing an operationalizable approximation to Distributionally Robust Optimization (DRO) that does not require distributional assumptions, to the best of our knowledge, the first constrained instance-level formulation in the LLM-based AHD literature. • Adversarial instance search. We propose RAISE, a bi-level evolutionary framework with an LLM-driven outer loop for heuristic evolution and an LLM-free inner loop for worstcase instance discovery. The inner loop employs a basis distribution parameterization with boundary projection to efficiently explore the uncertainty set without additional LLM cost. • Empirical validation. Experiments on three online combinatorial optimization tasks—OBP, OJSP, and OVRP—across five distribution families and 95 datasets with different distributions show that RAISE significantly improves the performance on shifted distributions among all learned methods, while existing LLM-based AHD methods degrade by up to 19× under distribution shift and RAISE remains competitive on nominal instances.

2

Methodology

2.1

Problem Formulation: Robust Automated Heuristic Design (RAHD)

Let H denote the space of heuristics. For a given design task, each heuristic h ∈ H is evaluated on an instance s via a black-box oracle eval(h, s) ∈ R. Without loss of generality, we define eval(h, s) so that larger values indicate better performance. For minimization objectives in our experiments, we negate the raw metric before evaluation. Standard AHD optimizes h∗ = arg max Es∼P0 [eval(h, s)] h∈H

under a fixed nominal distribution P0 (usually the distribution is not available, so it is a training instance set). This objective is brittle under distributional shift: a heuristic tuned to P0 may degrade sharply when test instances are drawn from a nearby but a different distribution [Wang et al., 2026, Shi et al., 2026, Liu et al., 2026]. To address this limitation, we adopt a distributionally robust formulation. Given a nominal instance set S = {s1 , . . . , sm } ∼ P0 and a robustness radius ε > 0, we define the uncertainty set  Bε (P0 ) = P : d(P, P0 ) ≤ ε , 2

Outer Loop: LLM-driven Heuristic Evolution Heuristic Population

ℎ$

</>

ℎ#

New Population

LLM-based Operator Exploration (Generate diversity)

ℎ$%

Exploitation (Refine & combine)

</>

ℎ#%

</>

⋮ ℎ!

Inner Loop: Adversarial Instance Search Input: current best heuristic h

Heuristic 𝒉

% ℎ!

Robust Evaluation Instance Set S (from nominal + augmented)

𝐴(ℎ") 𝐴(ℎ#)

Better

New Instance Set S

Projection onto ε-ball of nominal instances

⋯

Augment

Evolutionary search

⋮

𝐴(ℎ$ )

Adversarial distribution

ε

⋯

</>

Aggregation Score

Evaluate performance of each heuristic

Mixture of 9 basis distributions

g ∈ [0,1]18

</>

⋮ </>

Gene vector

S ← S + s*

Re-score population

Worst-case instance s*

Worse

Selection Keep top-N heuristics

Information flow

Feedback

Outer loop (heuristic evolution)

Inner loop (adversarial search)

Figure 1: Overview of the RAISE framework: an LLM-driven outer evolutionary loop for heuristic design coupled with an LLM-free inner adversarial instance search that discovers worst-case instances within an ε-ball uncertainty set around the nominal training distribution.

where d is a distance measure between distributions. In the distribution-level formulation, d captures the discrepancy between P and P0 directly [Rahimian and Mehrotra, 2019]. There are different distance metrics (e.g., optimal transport metrics) Rahimian and Mehrotra [2019] can be employed. In practice, however, the true distribution P0 is unavailable; only a finite set of nominal instances is observed. We therefore adopt an instance-level approximation, defining the practical, robust AHD objective as h∗ = arg max ′ min eval(h, s′ ), (1) h∈H s ∈Bε (S)

where  Bε (S) = s : ∃ si ∈ S, d(s, si ) ≤ ε is the instance-level uncertainty set, where d(s, si ) is a distance metric between instances (e.g., any ℓp norm or problem-specific similarity measure). In this work, we simply adopt the normalized ℓ1 distance, n 1X |sa,j − sb,j |, d(sa , sb ) = n j=1 where n is the dimension of the instance feature vector. Equation (1) is an instance-level approximation of distributionally robust optimization, rather than optimizing directly over all distributions in Bε (P0 ). We emphasize that this instance-level formulation is an engineering approximation that operationalizes robustness intuitions, rather than a formal solution to distribution-level DRO. The robustness radius ε continuously interpolates between two limiting regimes: • When ε → 0, the uncertainty set Bε (S) collapses to the nominal instance set S itself, and the robust objective (1) reduces to standard AHD—finding the heuristic with the best aggregated performance on the nominal instance set S, as in EoH [Liu et al., 2024] and ReEvo [Ye et al., 2024]. We note a subtle distinction: existing works typically aggregate performance via the mean over S, whereas our formulation uses a worst-case (min) aggregation; nonetheless, the instance set itself remains unchanged. • When ε → ∞, the uncertainty set expands to cover all feasible instances, and the objective becomes fully distribution-agnostic: finding the heuristic with the best worst-case performance across any conceivable instance. 3

Algorithm 1 RAISE: Robust Adversary Instance Search for LLM-Based Heuristic Design Require: Task description T , nominal instances S0 , robustness radius ε, max samples Nmax , max pop size Pmax , refresh interval τ Ensure: Best heuristic h∗ 1: Initialize population P ← ∅, instance set S ← {S0 } 2: Sample Pmax heuristics, evaluate on S, initialize P 3: while nsamples < Nmax do 4: for each LLM operator do ▷ Outer loop 5: Sample parent(s) via rank-weighted selection from P 6: Generate candidate  h via the selected LLM  operator  K 7: Evaluate r̂(h) ← A eval(h, s′k ) k=1 8: Update P via survival selection 9: end for 10: if generation mod τ = 0 then ▷ Inner loop 11: Run inner adversarial search on current best h ∈ P and get worst instance s∗ ; 12: Update S ← S + s∗ 13: Re-score all heuristics in P on updated S 14: end if 15: end while 16: return h∗ ← arg maxh∈P r̂(h)

In practice, ε is set to a finite intermediate value, providing a good balance of robustness performance. Larger values of ε yield more conservative but shift-resilient heuristics, while smaller values recover the higher nominal performance of standard AHD. 2.2

RAISE Framework Overview

RAISE solves the minimax problem (1) through a bi-level evolutionary search with two coupled loops. The overall workflow is: 1. Start from a nominal instance set S and initialize a population of candidate heuristics. 2. Evolve heuristics using LLMs in the outer loop (LLM-based heuristic evolution) under a fixed set of adversarial evaluation instances. 3. Every τ generations, run the inner loop (adversarial instance search) to search for worstcase instances within the instance-level uncertainty set Bε (S). 4. Augment the instance set with the discovered hard instances, re-score the population, and continue the outer search. This alternating procedure approximates the minimax objective in (1): the inner loop exposes current failure modes of the heuristic population, and the outer loop adapts that population to those harder conditions. Algorithm 1 summarizes the interaction between the two loops. While formal convergence guarantees are not available, this alternating design has a clear operational interpretation: the method repeatedly identifies weaknesses of the current heuristic population and then selects against those weaknesses in the next outer-loop phase. 2.3

Outer Loop: LLM-Based Evolutionary Heuristic Search

The outer loop maintains a population P of heuristic programs and improves it using LLM operators (e.g., synthesize a new heuristic from two high-performing heuristics or perturb one heuristic to explore its variants). The specific prompts are introduced in Appendix B. After new heuristics are generated, each candidate h is evaluated on the current adversarial instance set S = {s′1 , . . . , s′K } ⊆ Bε (S) through the robust score
r̂(h) = A eval(h, s′1 ), . . . , eval(h, s′K ) , 4

where A is an aggregation function, K is the number of instances in current instance set. This finite-sample score serves as a practical surrogate for the worst-case objective in (1): as the inner loop enlarges S with progressively harder instances, high-scoring heuristics are precisely those that remain robust across the most challenging distributional shifts encountered so far. After evaluation, candidate heuristics are admitted into the population P via survival selection, retaining only the top-N heuristics. In principle, the aggregation A should be the minimum, in direct correspondence with the worst-case objective in (1). In practice, directly optimizing the minimum over a finite and evolving adversarial set can produce unstable selection dynamics dominated by individual outlier instances. We therefore use mean aggregation as a more stable surrogate objective during search while retaining adversarial instance augmentation to encourage robustness. 2.4

Inner Loop: LLM-Free Adversarial Instance Search

Every τ outer-loop generations, RAISE runs an LLM-free evolutionary search to find worst-case instances within Bε (S) for the current best heuristic h ∈ P. Each candidate adversarial instance is parameterized by an 18-dimensional gene vector g ∈ [0, 1]18 . Basis distribution encoding. The distribution of a candidate instance s′ is constructed as a mixture of nine basis distributions centered on the nominal: ! 9 X padv = clip p̄ + wi · ϕi (g), 0, 1 , (2) i=1

where p̄ is the mean nominal distribution, wi = gi − 0.5 are signed mixture weights, and ϕi are nine parametric basis distributions: Uniform, Small, Large, Center, Bimodal, Gaussian, Periodic, Poisson-like, and Peak (Appendix B.2 lists the detailed settings). The 18 genes split into two groups: 9 genes define the mixture weights and the remaining 9 control the shape hyperparameters and support length. This parameterization is expressive enough to represent diverse distributional shifts while remaining low-dimensional enough for efficient inner-loop search. Moreover, it is not anchored to any predefined distribution family. Epsilon-ball projection. After constructing padv via (2), we project it onto the ε-ball of the nearest nominal instance to enforce feasibility with respect to Bε (S): pproj = b +

ε · (padv − b), max(ε, d(padv , b))

where b is the closest nominal instance and d(·, ·) is the mean absolute difference defined in Section 2. When d(padv , b) is larger than ε, this projection places candidates near the boundary of Bε (S). Inner evolutionary algorithm. The inner search applies uniform crossover and Gaussian mutation over Gin generations on a population of Pin candidate gene vectors, with each vector evaluated by its ability to degrade the current best heuristic h. The worst-scoring feasible instances s∗ , minimizing eval(h, s′ ) within Bε (S), is appended to instance set S for the outer loop. For example, if the current best heuristic is specialized to small-item instances on OBP, the inner loop may discover a large-item-heavy distribution. The outer loop then reevaluates all heuristics on the new worst instance and favors candidates that remain effective under the new set S.

3

Experimental Studies

3.1

Experimental Setup

Tasks: 1) Online Bin Packing (OBP): n items arrive sequentially and are assigned to bins with a capacity C. Waste ratio = (bins used − lower bound)/lower bound (lower is better) is used as the objective. Trained on 5 Weibull instances [Romera-Paredes et al., 2024, Liu et al., 2024] and Tested with C ∈ {100, 200, 300, 400} and n ∈ {1000, 5000, 10000}. Each size combination has 5 different distributions (Uniform, Normal, Lognormal, Exponential, and Triangular) of items sizes, which results in 4 × 3 × 5 = 60 datasets for OBP. 2) Online Job Shop Scheduling (JSP): Jobs with 5

EoH MoH ReEvo

Waste Ratio (%)

10

OBP Performance and Robustness Under Distribution Shift (b) Mean ± Std Across Distributions

PartEvo EoH-S RAISE (ours)

EoH

14 MoH

8

5.77±2.79

ReEvo

RAISE max = 2.85% EoH/MoH/ReEvo: 5.9 8.6%

6

12

5.08±2.28

PartEvo

4.04±1.89

4 EoH-S

2.13±1.05

RAISE (ours)

1.95±0.87

Normal

Lognormal

Distribution

Exponential

Triangular

0

2

4

10 8 6 4

2 0 Uniform

(c) Avg. Waste over Shifted Distributions by Problem Size 16

5.45±1.89

Waste Ratio (%)

(a) Waste Ratio Across Distributions (Avg. over n {1k,5k,10k}, C {100 400})

2 6

8

Waste Ratio (%)

10

12

0

n = 1, 000

n = 5, 000

n = 10, 000

Problem Size (Avg. over 5 Distributions, 4 Capacities)

Figure 2: a) OBP waste ratio (the lower the better) across five distributions. b) Average results averaged on all 60 intance datasets. c) Averaged results across different sizes. Table 1: OBP waste ratio (%) on n=5000, C=200 instances. Best per distribution in bold. underline: second best. Lower is better. Method

Uniform

Normal

Lognormal

Exponential

Triangular

Avg.

BestFit FirstFit

1.717 2.928

2.803 4.096

1.914 2.049

0.300 0.450

3.006 3.696

1.948 2.644

EoH Liu et al. [2024] ReEvo [Ye et al., 2024] PartEvo [Hu and Zhang, 2025] EoH-S [Liu et al., 2026] MoH [Shi et al., 2026]

5.182 5.466 3.853 2.318 5.277

6.308 6.539 5.062 3.010 6.531

1.321 0.714 0.835 1.334 0.917

5.303 3.011 2.236 0.436 1.443

5.765 3.218 4.619 2.043 5.113

4.776 3.790 3.321 1.828 3.856

RAISEε=0.001 RAISEε=0.002 RAISEε=0.005 RAISEε=0.010

3.814 1.954 2.785 1.670

4.225 2.914 3.402 2.763

0.822 1.253 1.011 2.062

3.435 0.314 0.273 0.354

1.997 1.919 2.636 2.803

2.859 1.671 2.021 1.930

stochastic processing times are assigned to machines online. The objective is to minimize normalized makespan. Trained on 5 uniform distribution instances with 10 machines (m) and 20 jobs (j). Tested on problem sizes: m ∈ {10, 20} and j ∈ {20, 50}. With 5 distributions for each combination, there is in total 20 datasets. 3) Online Vehicle Routing (VRP): Customers with stochastic demands arrive online. The objective is to minimize normalized route length ratio. Trained on 5 instances with 10 vehicles (v) and 50 customers. Tested on problem sizes: v ∈ {5, 10, 15}. With 5 distribution for each combination, there is in total 15 datasets. The details of the three tasks and the five different distributions families are introduced in Appendix C. Baselines. Classical: BestFit and FirstFit (OBP); SPT, MinSlack, EDD (JSP); Nearest Feasible Insertion, Slack-Preserving Insertion, and Urgency-Weighted Insertion (VRP). LLM-based AHD: EoH [Liu et al., 2024], ReEvo [Ye et al., 2024], and PartEvo [Hu and Zhang, 2025], each optimised on the nominal training set. Robustness-aware LLM-based AHD: EoH-S [Liu et al., 2026] (top heuristic trained on 128 diverse instances) and MoH [Shi et al., 2026] (evaluated with the heuristic reported in the original paper). Settings. All LLM-driven methods use 1,000 samples and a population size of 10, with GPT-5-mini as the backbone. All methods train on 5 nominal instances except EoH-S, which uses 128 diverse instances as in the original paper. For RAISE, the refresh interval is τ = 5, with inner population Pin = 8 and Gin = 4 inner generations. We ablate four values ε ∈ {0.001, 0.002, 0.005, 0.010} on OBP and use ε = 0.002 for OJSP and OVRP. All hyperparameter and prompt details are in Appendix B. 3.2

Results Across Different Distributions

Figure 2 summarizes the results on OBP. Table 1 and Table 2 report average results over different distribution families and all twelve size configurations (n ∈ {1k, 5k, 10k}, C ∈ {100, 200, 300, 400}). 6

Table 2: Average waste ratio (%) over five out-of-distribution test sets (Uniform, Normal, Lognormal, Exponential, Triangular) on Online Bin Packing across twelve size configurations (n ∈ {1k, 5k, 10k}, C ∈ {100, 200, 300, 400}). Bold: best per column; underline: second best. Lower is better. C = 100

Method 1k

5k

C = 200 10k

1k

5k

C = 300 10k

B EST F IT F IRST F IT

2.645 1.882 1.608 2.636 1.948 1.675 3.570 2.527 2.147 3.629 2.644 2.287

EOH R E E VO PART E VO E O H-S MOH

6.272 5.457 4.921 2.766 6.217

3.799 2.591 2.453 1.741 2.632

3.037 1.827 2.155 1.428 1.937

8.137 7.880 5.093 2.969 9.310

4.776 3.790 3.321 1.828 3.856

1k

10k

2.656 1.984 1.725 3.633 2.691 2.341

3.829 8.511 2.878 9.322 2.885 6.001 1.513 3.091 2.867 11.355

RAISE (ours) 2.742 1.641 1.307 2.685 1.671 1.377

5k

C = 400

4.957 4.630 3.922 1.865 4.685

1k

5k

10k

2.644 2.006 1.732 3.647 2.714 2.367

4.092 8.761 5.093 4.158 3.486 10.187 5.053 3.801 3.208 6.613 4.434 3.527 1.531 3.203 1.899 1.563 3.550 13.016 5.640 4.199

2.629 1.731 1.478

2.735 1.841 1.566

LLM-based methods collapse under distribution shift. EoH reaches a waste ratio of 5.182% on the Uniform out-of-distribution (OOD) setting and 6.308% on Normal, while ReEvo reaches 5.466% and 6.539% respectively—both far above classical BestFit (1.717% / 2.803%). MoH, despite meta-level training, similarly collapses to 5.277% on Uniform and 6.531% on Normal, because it still relies on explicit multi-task training across predefined distributions. These results confirm that heuristics optimised on a fixed nominal distribution are highly specialised and brittle: even modest distributional shifts cause severe performance collapse. RAISE achieves consistent OOD robustness. RAISEε=0.002 attains an average shifted waste of 1.671%, outperforming all baselines including the most competitive robustness-aware method EoH-S (1.828%). The gap is most pronounced on the Exponential distribution, where RAISEε=0.005 achieves 0.273% versus EoH’s 5.303%, a 19.4× improvement. Crucially, EoH-S trains on 128 diverse instances, which is substantially more data than RAISE’s 5 nominal instances and is often impractical in real deployments where diverse training instances may be unavailable. Robustness generalises across problem scales. Table 2 shows that RAISE’s advantage is consistent across all twelve size configurations. Among all learned methods, RAISE achieves the best average OOD performance on 8 of 12 configurations and is the second best on the remaining four. Classical BestFit remains strong at n = 1k owing to the lower difficulty of small-instance packing, but RAISE decisively surpasses it at n ≥ 5k. These trends indicate that the robustness gains scale with instance size, aligning with the expectation that adversarial instance search becomes more informative as problem scale grows. 3.3

Results on Online JSP and VRP

Table 3 reports OOD performance on JSP and VRP. RAISE achieves the lowest normalised makespan on both machine scales (1.2266 at 10-machine, 1.2390 at 20-machine), outperforming EoH (1.2284 / 1.2539) and ReEvo (1.2422 / 1.2500). The advantage widens at the 20-machine scale (∆ = 0.011 over EoH), where the richer combinatorial structure amplifies the benefit of adversarial worst-case exposure during search. All LLM-based methods substantially outperform classical scheduling rules, confirming the overall effectiveness of LLM-driven heuristic design. RAISE attains the best normalised route-length ratio on the 5-vehicle (0.8931) and 10-vehicle (0.9263) configurations, outperforming EoH (0.9017 / 0.9273) and ReEvo (0.8966 / 0.9287). On the 15-vehicle configuration, RAISE (0.9570) is marginally behind EoH (0.9561) while outperforming all classical insertion heuristics by a clear margin. Performance differences among LLM-based methods are smaller on OVRP than on OBP, suggesting that OVRP’s combinatorial structure provides some inherent distributional diversity across the tested OOD families. Nonetheless, RAISE consistently matches or exceeds competing methods across all configurations, demonstrating that its robustnessaware search generalises to structurally different online combinatorial problems beyond bin packing. 7

Table 3: Performance on online JSP (normalized makespan, ↓) and online VRP (route-length ratio ↓). OJSP columns report the mean over 10 configurations (5 processing-time distributions × 2 job counts: 20j/50j); VRP columns report the mean over 5 demand distributions. Classical baselines are task-specific (— = not applicable). Bold: best per column; underline: second best. Online JSP (makespan ↓)

Online VRP (route length, ↓)

Method

10-machine

20-machine

5-vehicle

10-vehicle

15-vehicle

SPT MinSlack EDD NF-Insert SP-Insert UW-Insert

1.3561 1.3812 1.4233 — — —

1.3956 1.4346 1.4886 — — —

— — — 0.9347 0.8971 0.9000

— — — 0.9747 0.9360 0.9345

— — — 1.0028 0.9665 0.9670

EoH ReEvo

1.2284 1.2422

1.2539 1.2500

0.9017 0.8966

0.9273 0.9287

0.9561 0.9613

RAISE (ours)

1.2266

1.2390

0.8931

0.9263

0.9570

Score convergence

a -0.0120

Shifted distributions

b

Early

Middle

6

End

4

Density

Heuristic score

Best

-0.0160

Nominal Early Middle End

6 5697 Nominal

-0.0200

Mean

2 -0.0240 Shaded band: population variability

20

40

60

80

100

Population generation

0

0.08

0.16

0.24

0.32

Item size (normalised capacity)

Figure 3: Robust Adversary Search Convergence in Online Bin Packing: Left, Population Score Evolution Across Generations; Right, Adversarial Distribution Shift Over Search Stages

3.4

Convergence process

Figure 3 illustrates the characteristic convergence dynamics of RAISE on OBP under iterative robust adversary search. It shows the heuristic score (which should be larger the better, as a negative sign is used for illustration) as a function of the number of populations. The left panel presents the evolution of both the best heuristic score and the population mean over generations, revealing clear stage-wise behavior. Following an initial exploration phase, the population rapidly adapts to increasingly challenging evaluation instances, indicating progressive consolidation of robust design principles rather than isolated lucky discoveries. The ramp-up stage corresponds to the emergence of a new robust heuristic, while the decreasing stage reflects the exploration of new worst-case instances, during which the current heuristic score declines. The right panel further reveals that this convergence is accompanied by a systematic shift in the adversarial item-size distribution away from the nominal regime toward harder distributions that expose the weakness of conventional packing rules. Taken together, these trends show that RAISE does not simply optimize for nominal performance, but converges by continually identifying, attacking, and repairing distribution-sensitive failure modes, ultimately yielding heuristics that remain effective under shifted distributions.

4

Ablation Study

We ablate four design choices of RAISE by independently removing each from the full system: (i) w/o Robust Search removes the adversarial inner loop entirely, reducing RAISE to standard AHD; 8

Table 4: Ablation study on the Online Bin Packing (OBP) problem. We report the shifted waste ratio (%) under five out-of-distribution test settings, averaged over all problem sizes (n ∈ {1000, 5000, 10000}) and bin capacities (C ∈ {100, 200, 300, 400}). Lower is better; bold marks the best result per column. Waste Ratio (%) by Distribution Method Avg. Uniform Normal Lognormal Exponential Triangular w/o Robust Search w/o Dis. Constraints w/o Base Distributions w/o ε-Mapping

4.94 2.84 4.49 7.16

5.32 2.93 4.84 8.40

1.40 1.17 2.77 1.67

1.21 0.48 2.53 2.33

5.11 2.60 5.28 10.23

3.60 2.01 3.98 5.96

RAISE (full)

2.54

2.85

1.45

0.48

2.43

1.95

(ii) w/o Distribution Constraint removes the ε-ball constraint-checking but keep the adversary search, which is similar to Karimi et al. [2025] and Duan et al. [2025]; (iii) w/o Base Distributions replaces the nine-basis parameterisation with flat-vector crossover and Gaussian mutation; (iv) w/o ε-Mapping retains the nine-basis parameterisation but skips the boundary projection. Table 4 shows that removing any component consistently degrades OOD robustness. Three findings stand out. (i) The ε-mapping is the most critical component, its removal raising the average waste ratio by +4.01 (+205%) and causing catastrophic degradation on Triangular (2.43% → 10.23%). Without boundary projection, the inner loop identifies adversarial instances far outside the nominal regime, causing the outer loop to over-optimise for distributional extremes at the expense of typical cases. (ii) The nine-basis parameterisation provides essential distributional diversity. Replacing it with flat-vector mutation increases average waste by +2.03 (+104%), as plain perturbations fail to expose structurally different distributions such as skewed or bimodal item-size regimes. (iii) The distribution constraint adds stability. Removing the ε-ball projection while retaining the basis parameterisation yields 2.01% average waste, marginally above the full model (1.95%), confirming that constraining divergence from the nominal set prevents the search from concentrating on deployment-unlikely extremes. (iv) The robust search procedure is indispensable: ablating it entirely degrades average performance by +1.65 (+85%), confirming that adversarial worst-case evaluation during search is necessary for consistent cross-distribution reliability.

5

Conclusion

We presented RAISE, a bi-level evolutionary framework that incorporates adversarial instance search under a constrained instance-level uncertainty set into LLM-based AHD. By formalizing robust AHD as a minimax problem over an instance-level ε-ball uncertainty set, RAISE decouples the robustness mechanism from any predefined distribution family: an LLM-driven outer loop evolves heuristics while an LLM-free inner loop continuously surfaces worst-case instances within a principled neighborhood of the training distribution, allowing RAISE to train on as few as five nominal instances yet generalize across substantially shifted distributions at test time. Comprehensive experiments on OBP, Online JSP, and Online VRP demonstrate that existing LLM-based AHD methods degrade by up to 19× under distribution shift, while RAISE consistently achieves the best average out-of-distribution performance among all learned methods across 95 test datasets. Limitations and future work. RAISE currently models single-dimensional distributional shifts and is evaluated only on online combinatorial optimization; extending to multi-dimensional correlated shifts, offline and mixed settings, and broader problem classes are important next steps.

9

References Shai Ben-David, John Blitzer, Koby Crammer, Alex Kulesza, Fernando Pereira, and Jennifer Wortman Vaughan. A theory of learning from different domains. Machine learning, 79(1):151–175, 2010. Edmund K Burke, Michel Gendreau, Matthew Hyde, Graham Kendall, Gabriela Ochoa, Ender Özcan, and Rong Qu. Hyper-heuristics: A survey of the state of the art. Journal of the Operational Research Society, 64(12):1695–1724, 2013. Pham Vu Tuan Dat, Long Doan, and Huynh Thi Thanh Binh. Hsevo: Elevating automatic heuristic design with diversity-driven harmony search and genetic algorithm using llms. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 39, pages 26931–26938, 2025. Ruibo Duan, Yuxin Liu, Xinyao Dong, and Chenglin Fan. Ealg: Evolutionary adversarial generation of language model-guided generators for combinatorial optimization. arXiv preprint arXiv:2506.02594, 2025. Yong Liang Goh, Zhiguang Cao, Yining Ma, Jianan Zhou, Mohammed Haroon Dupty, and Wee Sun Lee. SHIELD: Multi-task multi-distribution vehicle routing solver with sparsity and hierarchy. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview. net/forum?id=6DJEaz1cCj. Ian J Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in neural information processing systems, 27, 2014. Qinglong Hu and Qingfu Zhang. Partition to evolve: Niching-enhanced evolution with llms for automated algorithm discovery. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. Frank Hutter, Holger H Hoos, Kevin Leyton-Brown, and Thomas Stützle. Paramils: an automatic algorithm configuration framework. Journal of artificial intelligence research, 36:267–306, 2009. Frank Hutter, Holger H Hoos, and Kevin Leyton-Brown. Sequential model-based optimization for general algorithm configuration. In International conference on learning and intelligent optimization, pages 507–523. Springer, 2011. Pantea Karimi, Dany Rouhana, Pooria Namyar, Siva Kesava Reddy Kakarla, Venkat Arun, and Behnaz Arzani. Robust heuristic algorithm design with llms. arXiv preprint arXiv:2510.08755, 2025. Fei Liu, Xialiang Tong, Mingxuan Yuan, and Qingfu Zhang. Algorithm evolution using large language model. arXiv preprint arXiv:2311.15249, 2023. Fei Liu, Xialiang Tong, Mingxuan Yuan, Xi Lin, Fu Luo, Zhenkun Wang, Zhichao Lu, and Qingfu Zhang. Evolution of heuristics: Towards efficient automatic algorithm design using large language model. In Proceedings of the 41st International Conference on Machine Learning, ICML 2024, 2024. Fei Liu, Yilu Liu, Qingfu Zhang, Tong Xialiang, and Mingxuan Yuan. Eoh-s: Evolution of heuristic set using llms for automated heuristic design. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 40, pages 37090–37098, 2026. Manuel López-Ibáñez, Jérémie Dubois-Lacoste, Leslie Pérez Cáceres, Mauro Birattari, and Thomas Stützle. The irace package: Iterated racing for automatic algorithm configuration. Operations Research Perspectives, 3:43–58, 2016. Zeyuan Ma, Hongshu Guo, Yue-Jiao Gong, Jun Zhang, and Kay Chen Tan. Toward automated algorithm design: A survey and practical guide to meta-black-box-optimization. IEEE Transactions on Evolutionary Computation, 2025. Zeyuan Ma, Yue-Jiao Gong, Hongshu Guo, Jiacheng Chen, Yining Ma, Zhiguang Cao, and Jun Zhang. Llamoco: Instruction tuning of large language models for optimization code generation. IEEE Transactions on Evolutionary Computation, 2026. Alexander Novikov, Ngân Vũ, Marvin Eisenberger, Emilien Dupont, Po-Sen Huang, Adam Zsolt Wagner, Sergey Shirobokov, Borislav Kozlovskii, Francisco JR Ruiz, Abbas Mehrabian, et al. Alphaevolve: A coding agent for scientific and algorithmic discovery. arXiv preprint arXiv:2506.13131, 2025. Nelishia Pillay and Rong Qu. Hyper-heuristics: theory and applications. Springer, 2018. 10

Rong Qu, Graham Kendall, and Nelishia Pillay. The general combinatorial optimization problem: Towards automated algorithm design. IEEE Computational Intelligence Magazine, 15(2):14–23, 2020. Hamed Rahimian and Sanjay Mehrotra. Distributionally robust optimization: A review. arXiv preprint arXiv:1908.05659, 2019. 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. Mathematical discoveries from program search with large language models. Nature, 625(7995):468–475, 2024. Shiori Sagawa*, Pang Wei Koh*, Tatsunori B. Hashimoto, and Percy Liang. Distributionally robust neural networks. In International Conference on Learning Representations, 2020. URL https://openreview.net/forum?id=ryxGuJrFvS. Yiding Shi, Jianan Zhou, Wen Song, Jieyi Bi, Yaoxin Wu, Zhiguang Cao, and Jie Zhang. Generalizable heuristic generation through LLMs with meta-optimization. In The Fourteenth International Conference on Learning Representations, 2026. URL https://openreview.net/forum?id= tIQZ7pVN6S. Thomas Stützle and Manuel López-Ibáñez. Automated design of metaheuristic algorithms. In Handbook of metaheuristics, pages 541–579. Springer, 2018. Weiwei Sun, Shengyu Feng, Shanda Li, and Yiming Yang. Co-bench: Benchmarking language model agents in algorithm search for combinatorial optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 40, pages 33126–33134, 2026. Niki Van Stein and Thomas Bäck. Llamea: A large language model evolutionary algorithm for automatically generating metaheuristics. IEEE Transactions on Evolutionary Computation, 29(2): 331–345, 2024. Rongzheng Wang, Yihong Huang, Muquan Li, Jiakai Li, Di Liang, Bob Simons, Pei Ke, Shuang Liang, and Ke Qin. Rethinking llm-driven heuristic design: Generating efficient and specialized solvers via dynamics-aware optimization. arXiv preprint arXiv:2601.20868, 2026. Haoran Ye, Jiarui Wang, Zhiguang Cao, Federico Berto, Chuanbo Hua, Haeyeon Kim, Jinkyoo Park, and Guojie Song. Reevo: Large language models as hyper-heuristics with reflective evolution. Advances in neural information processing systems, 37:43571–43608, 2024. Zhi Zheng, Zhuoliang Xie, Zhenkun Wang, and Bryan Hooi. Monte carlo tree search for comprehensive exploration in LLM-based automatic heuristic design. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=Do1OdZzYHr. Jianan Zhou, Yaoxin Wu, Zhiguang Cao, Wen Song, Jie Zhang, and Zhiqi Shen. Collaboration! towards robust neural methods for routing problems. Advances in Neural Information Processing Systems, 37:121731–121764, 2024.

11

Appendix Table of Contents A Related Work

12

A.1 LLM-Based Automated Heuristic Design . . . . . . . . . . . . . . . . . . . . . .

12

A.2 Hyper-Heuristics and Automated Algorithm Configuration . . . . . . . . . . . . .

13

A.3 Adversarial Training and Minimax Learning . . . . . . . . . . . . . . . . . . . . .

13

B Detailed Method Description

14

B.1 Bi-Level Search Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

B.2 Inner Adversarial Instance Search . . . . . . . . . . . . . . . . . . . . . . . . . .

14

B.3 LLM Evolutionary Operators . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

C Detailed Experimental Design and Settings

18

C.1 Benchmark Tasks and Instance Generation . . . . . . . . . . . . . . . . . . . . . .

18

C.2 Baseline Algorithm Descriptions . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

C.3 RAISE Hyperparameter Settings . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

D Comprehensive Experimental Results

22

D.1 Online Bin Packing (OBP) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

D.2 Online Job Shop Scheduling (JSP) . . . . . . . . . . . . . . . . . . . . . . . . . .

23

D.3 Online Vehicle Routing (VRP) . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

E RAISE Heuristic Functions: Code and Robustness Analysis

25

E.1 Online Bin Packing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

E.2 Online Job Shop Scheduling . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

E.3 Online Vehicle Routing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

F Computational Cost Analysis

30

F.1

RAISE Cost Breakdown . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

F.2

Comparison with Existing Frameworks . . . . . . . . . . . . . . . . . . . . . . .

30

A

Related Work

A.1

LLM-Based Automated Heuristic Design

The explosion of capable LLMs has unlocked a new paradigm for AHD in which the LLM serves as a black-box mutation and crossover operator over algorithm space. EoH [Liu et al., 2024] is one of the preliminary works on LLM-based AHD. It maintains a population of heuristic Python functions, each scored on a fixed set of problem instances, and uses four LLM prompt types to create, crossover, and mutate programs. EoH achieves state-of-the-art performance on online bin packing and traveling salesman problems under nominal distribution conditions. ReEvo [Ye et al., 2024] extends EoH with a dual-memory architecture: a short-term memory of recent high-scoring programs and a long-term memory of reflective insights generated by the LLM analyzing performance gaps. ReEvo shows that guided memory injection substantially improves sample efficiency, achieving 12

competitive results with fewer LLM calls. HSEvo [Dat et al., 2025] proposes a harmony search framework to balance exploitation and exploration. MCTS-AHD [Zheng et al., 2025] employs Monte Carlo tree search to organize heuristics in a tree structure, enabling more effective exploration of the search space. PartEvo [Hu and Zhang, 2025] incorporates feature-assisted niche construction within abstract search spaces, enabling the seamless integration of niche-based search strategies from evolutionary computation. Generalizable heuristic design EoH-S [Liu et al., 2026] introduces portfolio-based AHD: instead of seeking a single best heuristic, it discovers a diverse set of complementary heuristics with varying strengths. EoH-S is particularly relevant as a baseline in our distribution-shift experiments, as maintaining a diverse portfolio implicitly confers some cross-distribution robustness. MoH [Shi et al., 2026] leverages LLMs to iteratively refine a meta-optimizer that autonomously constructs diverse heuristic-optimizers through self-invocation, thereby eliminating reliance on a predefined evolutionary computation heuristic-optimizer. These constructed heuristic-optimizers subsequently evolve heuristics for downstream tasks, enabling broader heuristic exploration. MoH further employs a multi-task training scheme to promote generalization. A key limitation shared by all of the above methods is that they optimize for a fixed nominal instance distribution. While EoH-S and MoH propose different strategies to enhance generalization, both still rely on predefined instance sets and offer no guarantee of robust performance under unknown distributions. RAISE explicitly addresses this limitation via distributional robustness through adversarial instance search. Adversarial instance search. Karimi et al. [2025] employ an LLM to identify regions of the input space where a heuristic underperforms and to suggest targeted improvements. Duan et al. [2025] propose a co-evolutionary framework that jointly evolves both the heuristic algorithm and an instance generator using LLMs, with dynamically generated instances serving to enhance heuristic robustness throughout the evolutionary process. While both approaches incorporate adversarial search into LLM-driven AHD, they impose no explicit constraints on the instance distribution and rely on LLMs to generate adversarial instances directly, which limits controllability and incurs additional inference cost. In contrast, RAISE employs a constrained adversarial search that bounds the divergence between adversarial and nominal instance distributions, striking a principled balance between worst-case robustness and performance on the nominal distribution. Crucially, the worst-case instance search loop in RAISE is LLM-free, making it lightweight and ready to work with a broad range of LLM-driven AHD frameworks. A.2

Hyper-Heuristics and Automated Algorithm Configuration

Burke et al. [2013] surveys hyper-heuristics across two dimensions: (1) heuristic selection vs. heuristic generation, and (2) constructive vs. perturbative. RAISE belongs to the generative, constructive quadrant—it synthesizes new heuristic functions from scratch rather than selecting or composing existing low-level heuristics. Traditional hyper-heuristics rely on hand-crafted low-level heuristic libraries and meta-level selectors (e.g., simulated annealing, reinforcement learning) to adaptively combine them. This approach is interpretable and efficient but fundamentally limited by the quality of the low-level library. Algorithm configuration methods such as SMAC [Hutter et al., 2011], ParamILS [Hutter et al., 2009], and irace [López-Ibáñez et al., 2016] find optimal hyperparameter settings for fixed algorithm templates using model-based search. They operate in continuous or mixed-integer parameter spaces, enabling principled exploration, but they cannot design new algorithmic logic. A.3

Adversarial Training and Minimax Learning

DRO [Rahimian and Mehrotra, 2019] is a classical framework in operations research and machine learning for decision-making under distributional uncertainty. The key idea is to replace the expected objective EP0 [f (x)] with the worst-case expected objective supP ∈U EP [f (x)] over an ambiguity set U. 13

Adversarial training [Goodfellow et al., 2014] and minimax learning originate in the GAN framework, where a generator and discriminator are co-trained via a minimax objective. The connection to distributional robustness was made explicit in subsequent work: training on adversarially perturbed data provably improves worst-case generalization under the assumed perturbation model. Domain adaptation methods [Ben-David et al., 2010] study the theoretical conditions under which a model trained on a source distribution generalizes to a target distribution, providing H-divergencebased bounds on transfer error. These bounds motivate training procedures that minimize distribution shift, aligning with RAISE’s adversarial instance generation objective. Distributionally robust neural networks [Sagawa* et al., 2020] apply group DRO to improve worstcase accuracy across subgroups in neural network training, demonstrating that explicit worst-case group weighting outperforms uniform ERM under distribution shift. RAISE is conceptually analogous but operates in the program-synthesis domain: instead of gradient-based adversarial perturbation of network weights, RAISE uses an evolutionary inner loop over a parameterized distribution space, since heuristic programs are non-differentiable. The key distinction between RAISE and all gradient-based adversarial/DRO methods is the black-box nature of heuristic evaluation: there is no gradient of the performance metric with respect to instance parameters, necessitating derivative-free optimization.

B

Detailed Method Description

This appendix provides comprehensive documentation of RAISE’s components beyond what is possible in the main paper. We describe each component in detail, provide complete pseudocode for sub-routines, document all LLM prompt templates, give a full hyperparameter table, and discuss implementation details and theoretical properties. B.1

Bi-Level Search Architecture

RAISE’s bi-level structure separates two timescales: the outer loop runs continuously (every generation), while the inner loop fires periodically (every τ = 5 generations). This separation is motivated by computational efficiency: the inner adversarial search involves re-evaluating an inner population of 8 candidate distributions over the current best heuristic; amortizing this cost over τ outer generations keeps the overhead below 20% of total compute. At each inner loop invocation, the best heuristic h∗t from the current outer population is used as the fixed “victim” for the adversarial search. The search finds instances that minimize eval(h∗t , s′ ) subject to d(s′ , si ) ≤ ε for some si ∈ S, where d denotes the mean absolute difference defined in Section 2. After the inner search, the worst-case instance is appended to the current instance set S, and the entire outer population is re-scored under the updated instances. This re-scoring ensures consistent ranking of all heuristics under the same adversarial context. B.2

Inner Adversarial Instance Search

The inner evolutionary algorithm is a compact evolutionary search. Its goal is to find, within the feasible uncertainty set Bε (S), the instances that most severely degrade the performance of the current best heuristic h∗ . Rather than searching directly in the space of problem instances—which is high-dimensional and unstructured—the inner loop operates on a compact 18-dimensional gene space that parameterizes a rich family of item-size distributions via basis function mixtures. Each gene vector is decoded into a candidate adversarial distribution, projected onto the boundary of the uncertainty set to ensure feasibility, and evaluated by running h∗ on a sampled instance. The population is then improved over Gin generations using selection, crossover, and mutation. Algorithm 2 provides the complete pseudocode; the remainder of this section describes each step in detail. Step 1: Population initialization (line 1). The inner search begins by randomly initializing Pin = 8 gene vectors, each drawn independently from Uniform([0, 1]18 ). This uniform initialization ensures broad coverage of the gene space at the start of the search, giving the evolutionary process diverse starting points without imposing any prior bias toward a particular distributional shift. 14

Algorithm 2 Inner Adversarial Instance Search Require: Current best heuristic h∗ , nominal instances S = {s1 , . . . , sk }, radius ε, inner population size Pin , inner generations Gin Ensure: Worst-case feasible instance 1: Initialize gene population G = {g1 , . . . , gPin } with gi ∼ Uniform([0, 1]18 ) 2: for g = 1 to Gin do 3: for each gi ∈ G do ▷ Eq. (2) in main paper 4: padv ← decode and mix(gi ) 5: b ← arg minsj ∈S d(padv , sj ) ▷ Nearest nominal instance ε · (padv − b) ▷ Epsilon-boundary projection 6: pproj ← b + max(ε, d(p adv ,b)) ′ 7: si ← sample instance(pproj ) 8: scorei ← −eval(h∗ , s′i ) ▷ Adversarial score: lower eval = higher adversarial score 9: end for 10: Sort G by score descending; keep top ⌊Pin /2⌋ 11: Crossover: for each pair of selected parents (gi , gj ), create offspring via uniform crossover with pcross = 0.5 12: Mutation: apply Gaussian noise N (0, 0.122 ) to each gene with pmut = 0.35; clip to [0, 1]18 13: end for 14: return The worst instance by adversarial score from final generation

Step 2: Decoding and distribution construction (line 4). Each gene vector gi is decoded into a candidate adversarial item-size distribution padv via the basis-mixture model (Eq. (2) in the main paper). The first nine genes g1:9 define signed mixture weights wi = gi − 0.5 ∈ [−0.5, 0.5] over nine parametric basis distributions, while genes g10:17 control the shape hyperparameters of those bases and g18 adjusts the instance length within ±10% of the nominal. The resulting padv is a flexible, smooth distribution that can represent a wide range of shifts—heavy small-item tails, large-item concentrations, bimodal splits, periodic patterns, and others—while remaining low-dimensional enough for efficient evolutionary search. The nine basis distributions and their gene encodings are described in detail at the end of this section. Step 3: Nearest nominal instance and epsilon-ball projection (lines 5–6). The unconstrained distribution padv produced by decoding may fall anywhere in distribution space and need not satisfy the feasibility constraint d(padv , sj ) ≤ ε for any sj ∈ S. To enforce feasibility, we first identify the nearest nominal instance b = arg minsj ∈S d(padv , sj ) under the mean absolute difference d, and then project padv radially onto the ε-boundary of the ball centered at b: pproj = b +

ε · (padv − b). max(ε, d(padv , b))

When d(padv , b) > ε, this formula rescales the displacement vector (padv − b) so that the projected point lies exactly on the boundary at distance ε from b; when d(padv , b) ≤ ε, the point is already feasible and is returned unchanged (since max(ε, d) = ε and the scaling factor equals 1). Projecting onto the boundary rather than the interior is deliberate: the hardest adversarial distributions concentrate near the boundary, where the maximum allowed shift from the nominal is achieved. Interior points represent smaller, less extreme shifts and are therefore typically less damaging to the heuristic. Step 4: Instance sampling and adversarial scoring (lines 7–8). From the projected distribution pproj , a concrete problem instance s′i is generated. The instance is then passed to the black-box evaluator to compute eval(h∗ , s′i ), which measures how well h∗ performs on s′i . The adversarial score is defined as the negation, scorei = −eval(h∗ , s′i ), so that maximizing the adversarial score is equivalent to finding instances on which h∗ performs worst. Step 5: Selection (line 9). After scoring all Pin candidates in a generation, the population is ranked by adversarial score in descending order and the top half (i.e., ⌊Pin /2⌋ = 4 individuals) is retained as parents for the next generation. This truncation-selection scheme is computationally simple and well-suited to the short inner search horizon. The selected parents define the gene pool from which the next generation is produced via crossover and mutation. Step 6: Crossover (line 10). Each pair of selected parents (gi , gj ) produces one offspring via uniform crossover with probability pcross = 0.5: for each of the 18 gene positions independently, the 15

offspring inherits the value from gi with probability 0.5 and from gj with probability 0.5. Uniform crossover is chosen over single-point or two-point crossover because it treats all gene positions symmetrically, which is appropriate here given that the 18 genes encode heterogeneous quantities (mixture weights, shape parameters, and a length parameter) with no meaningful positional ordering. Step 7: Mutation (line 11). Each gene of every offspring is independently perturbed by additive Gaussian noise N (0, 0.122 ) with probability pmut = 0.35, and the result is clipped to [0, 1]18 to maintain valid gene values. The mutation standard deviation σ = 0.12 is chosen to produce moderate perturbations relative to the unit gene range: on average, a mutated gene shifts by approximately 0.12 × 0.35 ≈ 0.04 units, which translates to a meaningful but not catastrophic change in the decoded distribution. Step 8: Return worst-case instances (line 13). After Gin generations, the worst candidate instance with the highest adversarial scores across the final population are returned to the outer loop. It is appended to the instance set S, after which all heuristics in the outer population are re-scored on the updated S before evolution continues. Gene decoding detail. The 18-dimensional gene vector g is partitioned as: • g1:9 : basis function weights wi = gi − 0.5 ∈ [−0.5, 0.5] • g10:17 : shape hyperparameters for basis distributions (e.g., power p for Small/Large, mean/std for Gaussian, period for Periodic, peak position for Peak, center width for Center) • g18 : length gene (fractional adjustment to item count within ±10% of nominal) The nine basis distributions ϕi are: 1. Uniform: q ∼ Uniform(0, 1) 2. Small: q ∼ Beta(1, p) with shape p = 1 + 4 · g10 (favors small items) 3. Large: q ∼ Beta(p, 1) with shape p = 1 + 4 · g11 (favors large items) 4. Center: q ∼ Triangular(0, 0.5, 1) (concentrated around 0.5) 5. Bimodal: mixture of Uniform(0, 0.3) and Uniform(0.7, 1) 6. Gaussian: q ∼ Normal(µ, σ 2 ) clipped to [0, 1], with µ = g12 , σ = 0.1 + 0.3 · g13 7. Periodic: sinusoidal weight with frequency controlled by g14 8. Poisson-like: discretized Poisson CDF with rate λ = 1 + 9 · g15 9. Peak: spike at position g16 , spreading controlled by g17 B.3

LLM Evolutionary Operators

The outer loop employs five LLM operators from [Liu et al., 2024] with distinct roles, as described in Section 2: C REATE, E1, E2, M1, and M2. Table 5 summarizes their inputs, outputs, and invocation characteristics. All operators require the model to first describe the proposed heuristic in one sentence (enclosed in a boxed environment) before implementing the Python function, ensuring interpretable outputs alongside executable code. Table 5: Summary of RAISE’s five LLM evolutionary operators. Operator

Inputs

Purpose

C REATE E1 E2 M1 M2

Task description only k parents + task k parents + task One parent + task One parent + task

Generate novel heuristic from scratch Generate heuristic with a totally different form from all parents Extract backbone idea from parents; synthesize a new motivated variant Produce a modified version of the parent heuristic Perturb numerical parameters of the parent’s scoring function

Parents

The operators divide into two families. The exploration operators (C REATE, E1, E2) are designed to introduce structural novelty: C REATE requires no parent at all; E1 explicitly instructs the model to produce a form totally different from all provided parents; and E2 first asks the model to identify 16

0 k k 1 1

the common backbone idea across parents before synthesizing a new heuristic motivated by—but structurally distinct from—that backbone. The exploitation operators (M1, M2) make incremental changes to a single parent. M1 allows free modification of the heuristic’s logic, while M2 targets only the numerical parameter settings of the existing scoring function, enabling fine-grained local search without altering algorithmic structure. C REATE Prompt Template {task description} Create a new heuristic. 1. First, describe your heuristic in one sentence inside boxed {}. 2. Next, implement the following Python function: {function signature} Do not give additional explanations. E1 Prompt Template {task description} I have {k} existing heuristics with their codes as follows: {indivs} Please help me create a new heuristic that has a totally different form from the given ones. 1. First, describe your heuristic in one sentence. The description must be inside boxed {}. 2. Next, implement the following Python function: {function signature} Do not give additional explanations. E2 Prompt Template {task description} I have {k} existing heuristics with their codes as follows: {indivs} Please help me create a new heuristic that has a totally different form from the given ones but can be motivated from them. 1. Firstly, identify the common backbone idea in the provided heuristics. 2. Secondly, based on the backbone idea describe your new heuristic in one sentence. The description must be inside boxed {}. 3. Thirdly, implement the following Python function: {function signature} Do not give additional explanations. M1 Prompt Template {task description} I have one heuristic with its code as follows: {indivs} Please create a new heuristic that can be a modified version of the given one. 1. First, describe your heuristic in one sentence. The description must be inside boxed {}. 2. Next, implement the following Python function: {function signature} Do not give additional explanations.

17

M2 Prompt Template {task description} I have one heuristic with its code as follows: {indivs} Please identify the main heuristic parameters and assist me in creating a new heuristic that has different parameter settings of the score function provided. 1. First, describe your heuristic in one sentence. The description must be inside boxed {}. 2. Next, implement the following Python function: {function signature} Do not give additional explanations.

C

Detailed Experimental Design and Settings

This appendix provides complete documentation of the experimental setup, including benchmark task definitions, dataset construction procedures with exact parameterizations derived from the generation scripts, baseline algorithm configurations, RAISE hyperparameter settings, and the full evaluation protocol. C.1 C.1.1

Benchmark Tasks and Instance Generation Online Bin Packing (OBP)

Problem definition. In Online Bin Packing, items of integer size ai ∈ {1, . . . , C} arrive one at a time from an online sequence σ = (a1 , a2 , . . . , an ). Upon arrival, each item must be irrevocably assigned to an open bin (without exceeding the bin capacity C) or to a new bin. The objective is to minimize the total number of bins used. Formally, let B denote the number of bins used and P L = ⌈ i ai /C⌉ the lower bound (optimal offline solution). The waste ratio is: B−L WasteRatio = . L Lower values indicate better performance; a waste ratio of 0 means the online algorithm matches the offline lower bound. Training instances. RAISE’s evolutionary search is seeded with a nominal OBP instance set comprising I = 5 independent instances in Weibull distribution [Romera-Paredes et al., 2024, Liu et al., 2024], each with n = 5,000 items and bin capacity C = 200. These nominal instances constitute the initial adversarial instance set S0 , which RAISE augments dynamically during the search. All LLM-based baselines (EoH, MoH, ReEvo, PartEvo, EoH-S) are trained on this same nominal set to ensure a fair comparison. Test distribution parameterization. Testing evaluates cross-distribution generalization across five standard families. Item sizes are drawn from the continuous distribution, rounded to the nearest integer (⌊·⌉), and clipped to the valid range [1, C]. Table 6 gives the exact parameterizations used in the generation script. Table 6: OBP item-size distribution parameterizations (capacity C). All samples are rounded to integers and clipped to [1, C]. Distribution

Continuous draw

Parameters

Support

Uniform Normal Lognormal Exponential Triangular

ai ∼ U(1, C) ai ∼ N (µ, σ 2 ) ai ∼ LogNormal(µℓ , σℓ2 ) ai ∼ Exp(λ) + 1 ai ∼ Tri(ℓ, m, r)

— µ = 0.5C, σ = 0.2C µℓ = ln(0.25C), σℓ = 0.6 λ = 0.3C ℓ = 1, m = 0.35C, r = C

[1, C] clip to [1, C] clip to [1, C] clip to [1, C] [1, C]

Test suite structure. The OBP test suite spans all five distributions, three problem sizes (n ∈ {1,000, 5,000, 10,000}), and four bin capacities (C ∈ {100, 200, 300, 400}), yielding 5×3×4 = 18

60 dataset files, each containing I = 5 independently generated instances. The complete test suite thus covers 300 OBP test instances. Table 7 summarizes the combinatorial structure. Table 7: OBP test suite summary. Each (distribution, n, C) combination produces one dataset file with I = 5 instances. Problem size n

Distributions

Capacities C

Files

Instances

5 5 5

{100, 200, 300, 400} {100, 200, 300, 400} {100, 200, 300, 400}

20 20 20

100 100 100

60

300

1,000 5,000 10,000 Total

C.1.2

Online Job Shop Scheduling (JSP)

Problem definition. In Online Job Shop Scheduling (OJSP), J jobs must be processed on M machines. Each job j consists of exactly M operations, where operation (j, m) must be processed on machine m for an integer processing time pjm ∈ {1, . . . , 99}. The machine order within each job is a random permutation of all M machines. Jobs carry a release time rj (before which they cannot be started) and a due date dj . Jobs arrive online at their release time; upon arrival, the dispatching heuristic irrevocably assigns all operations of the job to machines. The objective is to minimize the makespan Cmax = maxj Cj , where Cj denotes the completion time of job j. The normalized makespan divides
Cmax by the theoretical lower bound LB = P P max maxj m pjm , ⌈ j,m pjm /M ⌉ ; lower values indicate better performance, with 1.0 denoting an optimal solution. Instance structure. Each JSP instance is generated as follows. 1. Operation durations. For each job j and machine m, the processing time pjm is drawn from the specified distribution (Table 8), with integer values in [1, 99]. 2. Machine order. Each job’s operation sequence is a uniformly random permutation of the M machine indices. P P 3. Planning horizon. H = 2 × max(maxj m pjm , ⌈ j,m pjm /M ⌉). 4. Release times. Release times are drawn with random gaps (probability 0.35 of a zero gap, otherwise a random positive increment), then rescaled to fit within a release window of 0.45 H; the resulting times are randomly permuted across jobs. P 5. Due dates. dj = rj + pj + δj , where pj = m pjm , and δj combines a multiplicative slack (Uniform(1.35, 1.75) × pj ), an additive congestion allowance (0.35 × p̄, with p̄ the mean job workload), and a random noise term bounded by a congestion buffer. Training instances. The nominal training set comprises I = 5 instances with M = 10 machines and J = 20 jobs. The operation durations are uniformly sampled. Test distribution parameterization. Table 8 specifies the operation-duration distributions for all five test families. The parameterizations mirror those used for OBP (normalized to the same [1, 99] range), enabling direct comparison of distributional effects across benchmarks. Table 8: JSP operation-duration distribution parameterizations (Dmax = 99). All samples are rounded to integers and clipped to [1, 99]. Distribution

Continuous draw

Parameters

Support

Uniform Normal Lognormal Exponential Triangular

pjm ∼ U(1, Dmax ) pjm ∼ N (µ, σ 2 ) pjm ∼ LogNormal(µℓ , σℓ2 ) pjm ∼ Exp(λ) + 1 pjm ∼ Tri(ℓ, m, r)

— µ = 0.5Dmax , σ = 0.2Dmax µℓ = ln(0.25Dmax ), σℓ = 0.6 λ = 0.3Dmax ℓ = 1, m = 0.35Dmax , r = Dmax

[1, 99] clip to [1, 99] clip to [1, 99] clip to [1, 99] [1, 99]

19

Test suite structure. The JSP test suite spans all five distributions, two machine counts (M ∈ {10, 20}), and two job counts (J ∈ {20, 50}), yielding 5 × 2 × 2 = 20 dataset files, each with I = 5 instances. The complete test suite covers 100 JSP test instances. Table 9 summarizes the structure. Table 9: JSP test suite: number of instances per (M, J, distribution) cell. Each cell contains I = 5 instances. Subtotal counts instances per distribution. Distribution

M = 10, J = 20

M = 10, J = 50

M = 20, J = 20

M = 20, J = 50

Uniform Normal Lognormal Exponential Triangular

5 5 5 5 5

5 5 5 5 5

5 5 5 5 5

5 5 5 5 5

Total

25

25

25

25

C.1.3

Online Vehicle Routing (VRP)

Problem definition. In Online Vehicle Routing (OVRP), a fleet of V homogeneous vehicles serves Nc = 50 customers from a central depot. Each customer i has an integer demand di ∈ {1, . . . , 10} and a location (xi , yi ) ∈ [0, 1]2 sampled from a clustered spatial process. Customers arrive online; upon arrival, the routing heuristic irrevocably assigns each customer to a vehicle’s route. All vehicles have a fixed capacity of Q = 40 units (for Nc = 50). The depot is located at the centroid (0.5, 0.5) of the unit square. The evaluation metric is 1/r, where r is the route-length ratio defined as the total served demand normalized by the route length relative to a reference; lower 1/r is better, with smaller values indicating more efficient routing. Instance structure. Each VRP instance is generated as follows. 1. Customer locations. Locations are drawn from a Gaussian mixture with K ∈ {2, 3} clusters, where each cluster center is sampled uniformly from [0.15, 0.85]2 , cluster standard deviation is drawn uniformly from [0.08, 0.18], and cluster weights follow a symmetric Dirichlet prior. All coordinates are clipped to [0, 1]2 . 2. Customer demands. Integer demands di are drawn from the specified distribution (Table 10), with di ∈ [1, 10]. 3. Vehicle capacity. Fixed at Q = 40 for all instances with Nc = 50 customers. Training instances. The nominal training set comprises I = 5 instances with 50 customers. The demands are uniformly sampled. Test distribution parameterization. Table 10 gives the customer-demand distributions for all five test families. Note that the Lognormal parameterization for VRP uses a slightly different location parameter (0.28 Dmax vs. 0.25 Dmax for OBP/JSP) and a reduced scale (σℓ = 0.55 vs. 0.60), reflecting the narrower integer demand range. Table 10: VRP customer-demand distribution parameterizations (Dmax = 10). All samples are rounded to integers and clipped to [1, 10]. Distribution

Continuous draw

Parameters

Support

Uniform Normal Lognormal Exponential Triangular

di ∼ U(1, Dmax ) di ∼ N (µ, σ 2 ) di ∼ LogNormal(µℓ , σℓ2 ) di ∼ Exp(λ) + 1 di ∼ Tri(ℓ, m, r)

— µ = 0.5Dmax , σ = 0.2Dmax µℓ = ln(0.28Dmax ), σℓ = 0.55 λ = 0.3Dmax ℓ = 1, m = 0.35Dmax , r = Dmax

[1, 10] clip to [1, 10] clip to [1, 10] clip to [1, 10] [1, 10]

Test suite structure. The VRP test suite spans all five distributions and three fleet sizes (V ∈ {5, 10, 15}), yielding 5 × 3 = 15 dataset files, each with I = 5 instances. The complete test suite covers 75 VRP test instances, all with Nc = 50 customers. 20

C.2

Baseline Algorithm Descriptions

C.2.1

Classical Baselines for OBP

BestFit. Each arriving item is placed in the open bin with the least remaining capacity that still accommodates the item (tightest-fit first). If no open bin can accommodate the item, a new bin is opened. BestFit is a deterministic algorithm with no hyperparameters and achieves a theoretical competitive ratio of approximately 1.7 in the worst case. FirstFit. Each arriving item is placed in the first open bin (in order of creation) that can accommodate it, opening a new bin only if no existing bin suffices. FirstFit is also deterministic with competitive ratio ≈ 1.7, but typically performs slightly worse than BestFit because it does not account for remaining capacity when selecting the bin. C.2.2

Classical Baselines for JSP

Shortest Processing Time (SPT). Upon job arrival, operations are scheduled to prioritize the job with the shortest total remaining processing time among all ready operations. SPT tends to minimize average completion time but can sacrifice makespan on large instances. Minimum Slack (MinSlack). Dispatching priority is assigned to the job with the smallest slack P value sj = dj − rj − m pjm , where dj is the due date and rj the release time. MinSlack targets due-date adherence but can produce suboptimal makespan when slack values cluster closely. Earliest Due Date (EDD). Operations are scheduled by ascending due date dj . EDD is optimal for minimizing maximum lateness in the single-machine case but offers no direct makespan guarantee on multi-machine instances. C.2.3

Classical Baselines for VRP

Nearest-Feasible Insertion (NF-Ins.). Each arriving customer is inserted at the cheapest feasible position (minimum route-length increase) in the vehicle route with the highest remaining capacity that can serve the customer’s demand. If no feasible vehicle exists, a new route is initialized. Slack-Preserving Insertion (SP-Ins.). Insertion cost is adjusted to penalize routes that leave little residual capacity, thereby preserving flexibility for future customer arrivals. The customer is assigned to the vehicle whose insertion minimizes the normalized insertion cost relative to remaining capacity slack. Urgency-Weighted Insertion (UW-Ins.). Insertion cost is weighted by a demand urgency factor, which increases the priority of high-demand customers. The customer is routed to the vehicle that minimizes the urgency-weighted insertion distance. C.2.4

LLM-Based Baselines

All LLM-based methods are initialized with the same task description and function signature, and receive the same total LLM sampling budget (Nmax = 1,000 samples). Training uses the nominal instance sets described in Section C.1. EoH [Liu et al., 2024]. Evolution of Heuristics maintains a population of heuristic Python functions and applies four LLM operator types (I 1: create, I 2: improve, E 1: crossover, E 2: evolve) with equal probability. Population size is 10; no memory component. The heuristic trained on nominal instances is evaluated directly on all shifted test distributions without any re-training. MoH [Shi et al., 2026]. MoH (Meta-objective Heuristics) employs an LLM to iteratively refine a meta-optimizer that autonomously constructs and diversifies heuristic-optimizers via selfinvocation, eliminating reliance on a fixed evolutionary computation template. The constructed heuristic-optimizers then evolve downstream heuristics. MoH further applies a multi-task training scheme across multiple problem instances to promote cross-instance generalization. 21

ReEvo [Ye et al., 2024]. Reflective Evolution augments EoH with a dual-memory architecture: a short-term memory of recent high-scoring programs and a long-term reflective memory generated by the LLM analyzing performance gaps. Memory injection guides subsequent LLM calls toward progressively better regions of heuristic space. Configuration: same population size and budget as EoH, plus memory components enabled. PartEvo [Hu and Zhang, 2025]. Partition-based Evolution incorporates feature-assisted niche construction within an abstract search space, enabling seamless integration of niche-based strategies from evolutionary computation. Niching encourages maintenance of a diverse heuristic population, mitigating premature convergence. EoH-S [Liu et al., 2026]. EoH-Stochastic introduces portfolio-based AHD: rather than seeking a single best heuristic, it discovers a diverse set of complementary heuristics forming a Pareto-optimal portfolio. EoH-S is particularly relevant as a distribution-shift baseline because maintaining a diverse portfolio implicitly confers some cross-distribution robustness. Results in the main paper report the best single heuristic from the EoH-S portfolio, evaluated identically to other methods. C.3

RAISE Hyperparameter Settings

Table 11 provides a complete listing of all RAISE hyperparameters, including the values used in experiments, the range explored during preliminary tuning, and a sensitivity classification. Table 11: Complete hyperparameter settings.

D

Parameter

Symbol

Value Used

Outer loop parameters Population size Max samples (budget) Refresh interval

Pmax Nmax τ

10 1000 5

Robustness parameters Robustness radius (OBP) Robustness radius (JSP) Robustness radius (VRP) Robust aggregation

ε ε ε A

{0.001, 0.002, 0.005, 0.010} {0.002} {0.002} mean

Inner search parameters Inner population size Inner generations Gene dimension Mutation std. dev. Mutation probability Crossover probability

Pin Gin — σ pmut pcross

8 4 18 0.12 0.35 0.5

Comprehensive Experimental Results

This appendix reports the complete numerical results for all three benchmark problems across every evaluated distribution and configuration. All values are computed directly from the raw evaluation CSV files. Bold marks the best result in each column. RAISE denotes RAISEε=0.002 throughout unless noted otherwise. D.1

Online Bin Packing (OBP)

The metric is shifted waste ratio (%), i.e., (1 − filled capacity/bin capacity) × 100; lower is better. Per-size breakdown. capacities.

Tables 12–14 detail results at each problem size, averaged over the four bin

22

Table 12: OBP shifted waste ratio (%) at n = 1000, averaged over C ∈ {100, 200, 300, 400}. Lower is better. Method

Uniform

Normal

Lognormal

Exponential

Triangular

Avg

EoH MoH ReEvo PartEvo EoH-S

9.97 12.11 10.89 7.72 5.21

9.82 13.75 11.82 8.35 4.50

3.27 4.23 3.23 2.12 1.43

6.36 6.02 5.86 2.67 1.37

10.18 13.77 9.26 7.42 2.94

7.92 9.97 8.21 5.66 3.09

RAISE

3.96

3.74

1.69

0.86

3.24

2.70

Table 13: OBP shifted waste ratio (%) at n = 5000, averaged over C ∈ {100, 200, 300, 400}. Lower is better.

D.2

Method

Uniform

Normal

Lognormal

Exponential

Triangular

Avg

EoH MoH ReEvo PartEvo EoH-S

5.02 5.72 5.80 4.17 2.56

5.88 7.01 6.76 5.33 3.11

1.32 1.06 0.89 1.04 1.17

5.25 1.68 3.17 2.40 0.45

5.81 5.54 3.46 4.73 1.81

4.66 4.20 4.02 3.53 1.82

RAISE

1.94

2.86

1.38

0.35

2.08

1.72

Online Job Shop Scheduling (JSP)

The metric is normalized makespan (lower is better). Tables 15 and 16 report results for 10-machine and 20-machine instances respectively. Columns are grouped by distribution; within each group the two job-count variants (20-job and 50-job) are shown side by side. The Avg column averages across all ten (distribution, job-count) cells within the table. D.3

Online Vehicle Routing (VRP)

The metric is 1/r where r is the achieved route-length ratio (total served demand divided by vehicle capacity, weighted by route length) relative to a reference; lower 1/r is better. In two instances (Uniform and Triangular with 15 vehicles) the Nearest-Feasible baseline achieves r < 1.0, giving 1/r > 1.0; these are marked with †. Tables 17–19 report results for each fleet size separately. Column abbreviations: U=Uniform, N=Normal, LN=Lognormal, Exp=Exponential, Tri=Triangular.

Table 14: OBP shifted waste ratio (%) at n = 10000, averaged over C ∈ {100, 200, 300, 400}. Lower is better. Method

Uniform

Normal

Lognormal

Exponential

Triangular

Avg

EoH MoH ReEvo PartEvo EoH-S

3.94 4.53 4.75 3.32 2.15

4.02 5.03 4.95 3.87 2.12

0.79 0.52 0.44 0.61 1.08

5.00 0.86 2.54 1.90 0.28

5.15 4.75 2.31 5.02 1.71

3.78 3.14 3.00 2.94 1.47

RAISE

1.72

1.94

1.30

0.25

1.96

1.43

23

Table 15: JSP normalized makespan on 10-machine instances (lower is better). Two job counts (20j / 50j) per distribution. Avg is the mean over all ten columns. Uniform Method

20j

50j

SPT 1.3980 1.2690 MinSlack 1.4158 1.2612 EDD 1.5291 1.2904

Normal 20j

Lognormal

Exponential

Triangular

20j

50j

20j

50j

20j

50j

Avg

1.5293 1.4965 1.5347

1.2586 1.2961 1.3229

1.3885 1.4965 1.5383

1.2723 1.3561 1.3736

1.4360 1.4203 1.4936

1.2570 1.2902 1.2870

1.3561 1.3812 1.4233

50j

1.5062 1.2456 1.5065 1.2727 1.5372 1.3260

EoH ReEvo

1.2456 1.1513 1.2805 1.1918 1.2952 1.1957 1.2905 1.2618 1.1486 1.2940 1.2238 1.3326 1.2205 1.2750

1.1689 1.2912 1.1730 1.2284 1.2006 1.2692 1.1954 1.2422

RAISE

1.2378 1.1696

1.1809

1.2819 1.2135

1.2748

1.2143

1.2795

1.2448

1.1693

1.2266

Table 16: JSP normalized makespan on 20-machine instances (lower is better). Two job counts (20j / 50j) per distribution. Avg is the mean over all ten columns. Uniform Method

20j

50j

SPT 1.3117 1.4749 MinSlack 1.3051 1.4961 EDD 1.3278 1.5939

Normal 20j

Lognormal

Exponential

Triangular

20j

50j

20j

50j

20j

50j

Avg

1.2630 1.3671 1.4686

1.4332 1.4633 1.5498

1.2948 1.3735 1.4024

1.4193 1.5327 1.5830

1.3122 1.3988 1.4006

1.5199 1.5009 1.5844

1.3956 1.4346 1.4886

50j

1.3784 1.5482 1.3877 1.5204 1.3804 1.5948

EoH ReEvo

1.2109 1.2611 1.2557 1.3148 1.2486 1.1840 1.3035 1.2412 1.3241 1.2250

1.2950 1.2047 1.2815 1.1813

1.2517 1.2072 1.2850 1.2539 1.2394 1.2215 1.2989 1.2500

RAISE

1.1752 1.2755 1.2210 1.3128

1.2677

1.2676

1.2019

1.1936

1.1896

1.2850

1.2390

Table 17: VRP inverse route-length ratio 1/r, fleet size 5 vehicles. Lower is better; bold marks the best per column. Method

U

N

LN

Exp

Tri

Avg

NF-Ins. SP-Ins. UW-Ins.

0.9423 0.9009 0.9060

0.9262 0.8930 0.8946

0.9249 0.8935 0.8946

0.9359 0.8956 0.8993

0.9444 0.9024 0.9057

0.9347 0.8971 0.9000

EoH ReEvo

0.9049 0.9034

0.8990 0.8901

0.8987 0.8921

0.9004 0.8952

0.9056 0.9022

0.9017 0.8966

RAISE

0.8974

0.8901

0.8890

0.8919

0.8970

0.8931

Table 18: VRP inverse route-length ratio 1/r, fleet size 10 vehicles. Lower is better; bold marks the best per column. Method

U

N

LN

Exp

Tri

Avg

NF-Ins. SP-Ins. UW-Ins.

0.9895 0.9398 0.9440

0.9677 0.9285 0.9285

0.9638 0.9417 0.9283

0.9668 0.9309 0.9291

0.9859 0.9389 0.9428

0.9747 0.9360 0.9345

EoH ReEvo

0.9297 0.9321

0.9235 0.9198

0.9230 0.9341

0.9267 0.9262

0.9337 0.9312

0.9273 0.9287

RAISE

0.9345

0.9228

0.9215

0.9211

0.9316

0.9263

24

Table 19: VRP inverse route-length ratio 1/r, fleet size 15 vehicles. Lower is better; bold marks the best per column. † marks entries where the raw route-length ratio r < 1. Method

E

U

N

LN

NF-Ins. SP-Ins. UW-Ins.

†

1.0180 0.9745 0.9726

0.9931 0.9553 0.9537

0.9877 0.9732 0.9713

EoH ReEvo

0.9637 0.9668

0.9529 0.9631

RAISE

0.9709

0.9508

Exp

Tri

Avg

0.9904 0.9562 0.9639

†

1.0249 0.9734 0.9736

1.0028 0.9665 0.9670

0.9496 0.9618

0.9494 0.9516

0.9651 0.9634

0.9561 0.9613

0.9499

0.9487

0.9646

0.9570

RAISE Heuristic Functions: Code and Robustness Analysis

This appendix presents the exact Python heuristic functions synthesized by RAIS across three combinatorial optimization tasks: Online Bin Packing (OBP), Online Job Shop Scheduling (JSP), and Online Vehicle Routing (VRP). For each task we list the code produced by the LLM evolutionary loop and then provide a mechanistic analysis explaining why the learned design choices confer distribution-shift robustness. E.1

Online Bin Packing

E.1.1

Heuristic Code

The heuristic is a priority(item, bins) function: given an incoming item of size item and a vector bins of residual capacities of open bins, it returns a score vector whose argmax selects the bin to use (negative infinity for infeasible bins). 1 2 3

import numpy as np def priority ( item : float , bins : np . ndarray ) -> np . ndarray : import numpy as np

4 5 6

eps = 1e -12 NEG_INF = -1 e9

7 8 9

bins = np . asarray ( bins , dtype = float ) n = bins . size

10 11 12 13 14

if n == 0: return np . array ([] , dtype = float ) if item <= eps : return np . full (n , NEG_INF , dtype = float )

15 16 17 18 19 20

# F e a s i b i l i t y mask feasible = bins >= item - 1e -12 scores = np . full (n , NEG_INF , dtype = float ) if not np . any ( feasible ) : return scores

21 22 23 24 25 26 27 28 29 30

b_all = bins . copy () # Rank p r e f e r e n c e : favor packing into smaller bins first ( n o r m a l i z e d rank ) order = np . argsort ( b_all ) ranks = np . empty_like ( order ) ranks [ order ] = np . arange ( n ) if n > 1: rank_norm_all = 1.0 - ( ranks . astype ( float ) / ( n - 1) ) else : rank_norm_all = np . ones (n , dtype = float )

31 32 33

b = b_all [ feasible ] rank_norm = rank_norm_all [ feasible ]

34

25

35 36

# R e s i d u a l after packing r = np . maximum ( b - item , 0.0)

37 38 39 40

# C o n s e r v a t i v e expected - slot size for future items mu = max (0.6 * item , 1e -12) n_slots = np . floor ( r / mu ) . astype ( int )

41 42 43 44 45

# Geometric - series r o b u s t n e s s : d i m i n i s h i n g returns per extra slot q = 0.65 robustness = (1.0 - np . power (q , n_slots + 1) ) / (1.0 - q ) - 1.0 robustness = np . maximum ( robustness , 0.0)

46 47 48

# T i g h t n e s s : how well the item fills the bin now tightness = item / ( b + eps )

49 50 51 52

# Mid - sliver penalty : r e l a t i v e l e f t o v e r in [0 ,1] frac_left = r / ( b + eps ) m i d _ s l i v e r _ p e n a l t y = np . exp ( -8.0 * frac_left )

53 54 55

# Exact - fit bonus exact_fit = ( r <= 1e -9) . astype ( float )

56 57 58 59 60 61

w_rank = 0.72 w_rob = 0.85 w_tight = 0.95 w_pen = 0.55 w_exact = 2.5

62 63 64 65 66 67 68 69 70

combined = ( w_rank * rank_norm * (1.0 + 0.5 * robustness ) + w_tight * tightness + w_rob * ( robustness / (1.0 + robustness ) ) + w_pen * m i d _ s l i v e r _ p e n a l t y + w_exact * exact_fit ) combined = combined - 0.002 * r # small tie - breaker

71 72 73 74 75 76 77 78

# N o r m a l i z e to [0 , 1] across f e a s i b l e bins c_min = float ( np . min ( combined ) ) c_max = float ( np . max ( combined ) ) if c_max - c_min < 1e -12: norm = np . full_like ( combined , 0.5 , dtype = float ) else : norm = ( combined - c_min ) / ( c_max - c_min )

79 80 81

scores [ feasible ] = norm return scores

Listing 1: RAISE heuristic for Online Bin Packing. E.1.2

Robustness Analysis

Five design choices in Listing 1 jointly explain why RAISE generalises across Uniform, Normal, Lognormal, Exponential, and Triangular test distributions. (1) Relative rather than absolute metrics. Tightness is item/(b+ϵ) and the mid-sliver penalty is applied to the fractional leftover r/b, not to absolute residual sizes. Because both signals are scale-normalised by the current bin capacity, they behave consistently whether items are drawn from a concentrated distribution (Normal, Lognormal) or a heavy-tailed one (Exponential, Lognormal). An absolute-residual formulation would fire differently depending on whether C=100 or C=400; the relative formulation is invariant. (2) Conservative geometric-series robustness term. The term µ = 0.6 × item defines a conservative lower bound for future items and counts how many such items could still fit in the leftover space r, then weighs them with geometric decay q=0.65. This models future packing utility without assuming 26

any particular distribution: it uses only the current item size as an anchor. The diminishing-return geometry caps the bonus for large leftovers, preventing over-favouring empty bins when many small items are expected—a pattern that penalises heuristics trained solely on Uniform distributions when tested on heavy-tailed Lognormal or Exponential instances. (3) Rank normalisation as a distribution-agnostic prior. Preferring bins with smaller residuals (rank-normalised to [0, 1] with 1 for the smallest) embodies Best-Fit logic, which is known to be nearoptimal in the average case for any continuous item distribution. Combining it multiplicatively with the robustness term (1+0.5 × robustness) selectively boosts bins that are both tight and capable of accommodating future items. (4) Strong exact-fit bonus. The wexact =2.5 bonus saturates the score whenever r≤10−9 , ensuring that exact-fit opportunities—which are relatively common under Uniform and Triangular distributions where item sizes span the full range—are always exploited. This prevents waste on instances where exact fits are available but a pure tightness heuristic might slightly prefer a near-fit. (5) Per-call score normalisation. Normalising all feasible bin scores to [0, 1] before returning prevents the absolute magnitude of the score from depending on the number of open bins or on the distribution of residual capacities. In late-game situations (few bins, all nearly full) the score landscape stays numerically stable, whereas un-normalised scores can exhibit vanishing gradients that lead to nearly-random selection. E.2 E.2.1

Online Job Shop Scheduling Heuristic Code

The heuristic is a priority(current time, candidate operations) function: given the current simulation clock and a matrix of eligible operations (columns: machine id, duration, release time, due date, remaining work, remaining operations), it returns a score vector whose argmax selects the operation to dispatch next. 1 2 3 4 5 6

import numpy as np def priority ( current_time : int , c a n d i d a t e _ o p e r a t i o n s : np . ndarray ) -> np . ndarray : ops = np . asarray ( candidate_operations , dtype = float ) if ops . size == 0: return np . array ([] , dtype = float )

7 8 9 10 11 12

duration = release_time = due_date = rem aining_ work = remaining_operations

ops [: , 1] ops [: , 2] ops [: , 3] ops [: , 4] = ops [: , 5]

13 14 15 16

start_time = np . maximum ( current_time , release_time ) co mp l et io n_ t im e = start_time + duration slack = due_date - c om pl e ti on _t i me

17 18 19 20 21 22 23

# Urgency : small / n e g a t i v e slack gets high p r i o r i t y urgency = np . where ( slack < 0 , 1.0 + ( - slack ) / (1.0 + duration ) , # overdue : grows with tardiness 1.0 / (1.0 + slack ) # on - time : d i m i n i s h i n g as slack grows )

24 25 26

# SRPT - like f i n i s h i n g : favour small r e m a i n i n g work finishing = 1.0 / (1.0 + r emaining _work )

27 28 29

# D o w n s t r e a m risk : more r e m a i n i n g ops = > higher p r i o r i t y ops_density = r e m a i n i n g _ o p e r a t i o n s / (1.0 + r e m a i n i n g _ o p e r a t i o n s )

30 31

# D u r a t i o n penalty : p e n a l i s e long o p e r a t i o n s unless urgent

27

32

d u r a t i o n _ p e n a l t y = duration / (1.0 + remai ning_wo rk )

33 34 35 36 37

w_urgency w_finish w_ops w_dur

= = = =

3.0 2.0 1.5 0.75

38 39 40 41 42

score = ( w_urgency * + w_finish + w_ops - w_dur

urgency * finishing * ops_density * d u r a t i o n _p e n a l t y )

43 44 45 46 47

score = np . nan_to_num ( score , nan =0.0 , posinf = np . finfo ( float ) . max , neginf = np . finfo ( float ) . min ) return score

Listing 2: RAISE heuristic for Online Job Shop Scheduling.

E.2.2

Robustness Analysis

(1) Dynamic, state-derived urgency. Slack = due date − (max(current time, release time) + duration) is recomputed at every scheduling step from the actual system clock and operation data. No distributional assumption enters: whether job arrival times follow a Poisson process, a bursty heavy-tail, or a periodic pattern, the urgency signal correctly reflects the current tightness. This is in contrast to heuristics that encode timing expectations calibrated to the training distribution, which degrade when inter-arrival statistics change. (2) Asymmetric urgency for overdue vs. on-time jobs. When slack ≥0, urgency = 1/(1 + slack) decreases gently as slack grows, de-prioritising comfortable jobs. When slack < 0, urgency = 1 + |slack|/(1 + duration) increases linearly with tardiness while scaling inversely with duration— so a short overdue operation gets more additional priority than a long one that is barely overdue. This asymmetry is correct regardless of the due-date distribution: it never lets an already-late job be permanently deprioritised by a deeply-backlogged queue. (3) SRPT-like finishing component. Shortest Remaining Processing Time (SRPT) is a provably optimal policy for minimising mean completion time in single-machine scheduling and performs robustly across exponential, uniform, and heavy-tailed service-time distributions. The term 1/(1+remaining work) approximates SRPT for the multi-machine online setting, providing a distribution-independent backbone. (4) Downstream cascade risk. ops density = R/(1+R) where R is the number of remaining operations monotonically increases with R but saturates below 1. Prioritising jobs with more remaining operations reduces the risk of blocking downstream machines—a property whose value is independent of whether job sizes are drawn from a uniform or heavy-tailed distribution. (5) Duration penalty modulated by remaining work. duration penalty = d/(1+remaining work) penalises long operations only when there is little remaining work left in the job. A job with large remaining work is less penalised even if its current operation is long, because blocking that machine is justified by the downstream workload. This modulation prevents the heuristic from systematically disfavouring heavy jobs, which are more prevalent under Lognormal or Pareto arrival distributions. E.3 E.3.1

Online Vehicle Routing Heuristic Code

The heuristic is a priority(customer features, candidate insertions) function. customer features encodes the arriving customer (demand, progress, etc.); candidate insertions is a matrix of feasible insertion positions with columns: route id, 28

position, delta distance, residual capacity, current load, route size, customer dist, projected route. The function returns a score vector whose argmax selects the best insertion. 1 2 3 4

import numpy as np def priority ( c u s t o m e r _ f e a t u r e s : np . ndarray , c a n d i d a t e _ i n s e r t i o n s : np . ndarray ) -> np . ndarray : import numpy as np

5 6 7

if c a n d i d a t e _ i n s e r t i o n s is None or c a n d i d a t e _ i n s e r t i o n s . size == 0: return np . array ([] , dtype = float )

8 9

cand = np . asarray ( candidate_insertions , dtype = float )

10 11 12 13 14 15 16

del ta_dist ance residual_capacity current_load route_size customer_dist point pr oj e ct ed _r o ut e

= = = = =

cand [: , cand [: , cand [: , cand [: , cand [: ,

2] 3] 4] 5] 6]

= cand [: , 7]

# # # # #

i n c r e m e n t a l d i s t a n c e cost r e m a i n i n g c a p a c i t y after i n s e r t i o n route u t i l i s a t i o n in [0 , 1] n o r m a l i s e d number of stops n o r m a l i s e d d i s t a n c e to i n s e r t i o n

# n o r m a l i s e d p r o j e c t e d route length

17 18 19

demand = float ( c u s t o m e r _ f e a t u r e s [2]) progress = float ( np . clip ( c u s t o m e r _ f e a t u r e s [4] , 0.0 , 1.0) )

20 21 22 23 24 25 26 27

# A d a p t i v e weights : routes fill w_travel = 0.3 + w_projected = 0.1 + w_prox = 0.5 w_capacity = 0.2 + w_balance = 0.1 + w_size = 0.4

shift from c a p a c i t y / balance to travel cost as 1.0 * progress 0.6 * progress

# i n c r e a s e s with p r o g r e s s # i n c r e a s e s with p r o g r e s s # fixed p r o x i m i t y weight 1.0 * (1.0 - progress ) # d e c r e a s e s with p r o g r e s s 0.8 * (1.0 - progress ) # d e c r e a s e s with p r o g r e s s # fixed route - size penalty

28 29 30

# Load - balance score : peaks at 50% u t i l i s a t i o n balance_score = 1.0 - 2.0 * np . abs ( current_load - 0.5)

31 32 33

# Demand scaling : larger demands amplify c a p a c i t y i m p o r t a n c e demand_scale = 1.0 + 0.5 * max (0.0 , demand )

34 35 36 37 38 39 40 41 42

score = ( w_travel * + w_projected + w_prox + w_capacity + w_balance - w_size )

( - del ta_dist ance ) * ( - p r oj ec t ed _r ou t e ) * ( - customer_dist ) * r e s i d u a l _ c a p a c i t y * demand_scale * balance_score * route_size

43 44

return np . asarray ( score , dtype = float )

Listing 3: RAISE heuristic for Online Vehicle Routing. E.3.2

Robustness Analysis

(1) Progress-adaptive weighting. All six weights in Listing 3 are functions of progress ∈ [0, 1], which measures how far route construction has advanced (fraction of customers already assigned). At progress=0: wtravel =0.3, wcapacity =1.2, wbalance =0.9. At progress=1: wtravel =1.3, wcapacity =0.2, wbalance =0.1. This phase transition is triggered by the actual construction state, not by any distributional parameter. It correctly reflects the decision-making reality: when routes are nearly empty, placing high-demand customers in capacity-rich routes prevents infeasibility regardless of the demand distribution; when routes are nearly full, minimising incremental travel cost is paramount. (2) Multi-objective coverage. The score balances five complementary signals—immediate travel cost, projected route length, proximity, residual capacity, load balance—plus a route-size penalty. No single distribution shift can simultaneously mislead all signals. For example, under a heavy-tailed 29

demand distribution the demand-scaled capacity term wcapacity ×residual×demand scale strongly steers large-demand customers away from nearly-full routes, preventing constraint violations; under a spatially clustered instance the proximity and projected-route terms guide insertion toward the geographically nearest feasible route. (3) Demand-adaptive capacity scaling. demand scale = 1 + 0.5 × max(0, d) amplifies the importance of capacity whenever the arriving customer carries a large demand. Rather than relying on the expected demand from a training distribution, this adapts to the observed demand value at run time. Under heavy-tailed demand distributions, where occasional very-large customers appear, this prevents catastrophic insertions that leave a route with insufficient capacity for subsequent customers. (4) Balance score peaks at 50% utilisation. The term 1 − 2|current load − 0.5| equals 1 at perfect half-load and −1 at completely empty or fully loaded routes. The 50% target is distributionagnostic: it maximises the expected remaining capacity buffer regardless of whether future demands are drawn from a uniform, Gaussian, or exponential distribution. Heuristics that use a distributionspecific optimal fill level degrade whenever the actual distribution shifts. (5) Route-size penalty encourages spreading. Penalising large route size discourages concentrating many stops on a single vehicle. This implicit load-balancing across routes is beneficial across all distributions because it maintains spare capacity on multiple routes simultaneously, reducing the chance of total infeasibility when a burst of high-demand customers arrives—a concern particularly acute under heavy-tailed or bursty demand distributions.

F

Computational Cost Analysis

This appendix quantifies the LLM API request budget, token consumption, and evaluation overhead of RAISE, and benchmarks these costs against five competing LLM-based AHD frameworks. All LLM-call counts follow the unified sample budget convention used across the comparison (each budget unit = one call to the LLM API resulting in one candidate heuristic program). F.1

RAISE Cost Breakdown

Table 20 lists the cost of a single RAISE run at the standard budget of Nmax = 1000 LLM samples, using a GPT-5-mini class model used in our experiments. F.2

Comparison with Existing Frameworks

Table 21 compares RAISE against EoH, ReEvo, PartEvo, MoH, and EoH-S at Nmax =1000 on online bin packing problem. Token overhead. All methods share the same 1,000-call budget. RAISE’s total token cost (≈2.82M) is comparable to EoH (≈2.80M). ReEvo and PartEvo are substantially more expensive in tokens (≈4.30M and ≈4.51M respectively) due to their multi-stage reflection and critique prompting strategies, which usually generate longer and more elaborate heuristic programs. EoH-S incurs a hidden evaluation overhead. EoH-S has a token cost comparable to RAISE, but its instance evaluation count (≈128,000) exceeds all other methods by over an order of magnitude. This arises from its post-hoc portfolio selection step, in which the full heuristic population is evaluated across a diverse set of 128 training instances spanning multiple distributions to select the top-K performers. Beyond the evaluation cost, this step requires prior knowledge of deployment distributions—often unavailable in practice. RAISE requires only nominal training-distribution instances, and delivers a single deployable heuristic rather than a portfolio. MoH incurs dual-task evaluation overhead. MoH’s ≈10,000 instance evaluations (double that of EoH) stem from its simultaneous training on two task variants—e.g., OBP under two item distributions—requiring every candidate heuristic to be scored on both variants each generation. This broadens coverage across the two training distributions but provides no adversarial robustness mechanism. 30

Table 20: RAISE per-run cost breakdown (Nmax =1000, OBP primary configuration). Token counts are computed from the system design (Appendix B); API costs use published list prices (GPT-5-mini: $0.25/$2.00 per 1M in/out tokens). Inner adversarial search uses CPU only (no LLM calls). Component

Description

Quantity

LLM API Calls (Outer Evolutionary Loop) Total budget Nmax Hard limit on LLM calls — C REATE New heuristic — E1 Crossover — E2 Backbone extraction + Crossover — M1 Mutation — M2 Mutation

1,000 ≈20 ≈245 ≈245 ≈245 ≈245

Token Consumption per LLM Call (approximate averages) C REATE Task desc. + template E1 Task + 2 parents + instruction + template E2 Task + 2 parent + instruction + template M1 Task + parent + instruction + template M2 Task + parent + instruction + template Weighted average

1,200 in / 500 out 2,850 in / 650 out 3,050 in / 650 out 1,850 in / 600 out 1,650 in / 600 out 2,200 in / 620 out

Total Token Budget (1,000 calls) Input tokens Output tokens

1, 000 × 2, 200 1, 000 × 620

≈2.2M ≈0.62M

Estimated API Cost GPT-5-mini

$0.25/$2.00 per 1M in/out

≈$ 1.79

Inner Adversarial Search Overhead (CPU only; no LLM) instance refreshes Every τ =5 gens; 200 gens total Inner evaluations per refresh Pin =8 pop × Gin =4 gens Re-scoring after refresh Pmax = 10 programs Outer evaluations per generation Pmax = 10 programs × 15 instances Total evaluations 40×(32+40) + 10×15×100 = Wall-Clock Time (OBP, 20 parallel samplers/evaluators) Measured runtime

40 32 40 100 ≈ 17,880 ≈2.0 h

Table 21: Per-run LLM cost comparison across AHD frameworks at Nmax =1000 samples. † Estimated from published architecture. API cost uses GPT-5-mini: $0.25/$2.00 per 1M in/out tokens Method

LLM Calls

Total Tokens

Est. Cost (mini)

#Instance Evals

EoH [Liu et al., 2024] ReEvo [Ye et al., 2024] PartEvo [Hu and Zhang, 2025] MoH [Shi et al., 2026]† EoH-S [Liu et al., 2026]

1,000 1,000 1,000 1,000 1,000

≈2.80M ≈4.30M ≈4.51M ≈3.10M ≈2.86M

≈$ 1.78 ≈$ 2.74 ≈$ 2.89 ≈$ 1.95 ≈$ 1.81

≈5,000 ≈5,000 ≈5,000 ≈10,000 ≈128,000

RAISE (ours)

1,000

≈2.82M

≈$ 1.79

≈17,880

31

Inner adversarial search is essentially free at the API level. RAISE’s inner worst instance search accounts for additional heuristic evaluations in CPU time only, which is independent of the number of deployment distributions, in contrast to EoH-S’s evaluation overhead, which scales linearly with the size of its 128-instance training set.

32