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 [1]–[5].

The recent integration of Large Language Models (LLMs) into evolutionary search has dramatically advanced AHD [6]–[17]. For example, [6] propose Evolution of Heuristics (EoH), which uses LLMs as evolutionary operators to iteratively generate, recombine, and refine heuristic functions for combinatorial optimization. [9] 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 [11] and PartEvo [12] 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 [14], [18]. This robustness gap is particularly problematic in real-world deployment, where problem instances are drawn from unknown or shifting distributions [4], [18], [19].

Recent work has attempted to address this limitation. [15] 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. [14] 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 [20] and [21] 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_{h}\; \min_{s' \in \mathcal{B}_\varepsilon(S)}\; \mathrm{eval}(h,\, s'),\] where \(\mathcal{B}_\varepsilon(S)\) denotes the ball of radius \(\varepsilon\) 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 \(\varepsilon\)-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 worst-case 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\times\) under distribution shift and RAISE remains competitive on nominal instances.

2 Methodology↩︎

2.1 Problem Formulation: Robust Automated Heuristic Design (RAHD)↩︎

Let \(\mathcal{H}\) denote the space of heuristics. For a given design task, each heuristic \(h \in \mathcal{H}\) is evaluated on an instance \(s\) via a black-box oracle \(\mathrm{eval}(h, s) \in \mathbb{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_{h \in \mathcal{H}}\, \mathbb{E}_{s \sim P_0}[\mathrm{eval}(h, s)]\] under a fixed nominal distribution \(P_0\) (usually the distribution is not available, so it is a training instance set). This objective is brittle under distributional shift: a heuristic tuned to \(P_0\) may degrade sharply when test instances are drawn from a nearby but a different distribution [14], [15], [18].

To address this limitation, we adopt a distributionally robust formulation. Given a nominal instance set \(S = \{s_1, \ldots, s_m\} \sim P_0\) and a robustness radius \(\varepsilon > 0\), we define the uncertainty set \[\mathcal{B}_\varepsilon(P_0) = \bigl\{P : d(P, P_0) \leq \varepsilon\bigr\},\] where \(d\) is a distance measure between distributions. In the distribution-level formulation, \(d\) captures the discrepancy between \(P\) and \(P_0\) directly [22]. There are different distance metrics (e.g., optimal transport metrics) [22] can be employed. In practice, however, the true distribution \(P_0\) 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_{h \in \mathcal{H}}\;\min_{s' \in \mathcal{B}_\varepsilon(S)}\; \mathrm{eval}(h,\, s'), \label{eq:robust95ahd}\tag{1}\] where \[\mathcal{B}_\varepsilon(S) = \bigl\{s : \exists\, s_i \in S,\; d(s, s_i) \leq \varepsilon\bigr\}\] is the instance-level uncertainty set, where \(d(s, s_i)\) is a distance metric between instances (e.g., any \(\ell_p\) norm or problem-specific similarity measure). In this work, we simply adopt the normalized \(\ell_1\) distance, \[d(s_{a}, s_{b}) = \frac{1}{n}\sum_{j=1}^{n} \lvert s_{a,j} - s_{b,j} \rvert,\] 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 \(\mathcal{B}_\varepsilon(P_0)\). 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 \(\varepsilon\) continuously interpolates between two limiting regimes:

• When \(\varepsilon \to 0\), the uncertainty set \(\mathcal{B}_\varepsilon(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 [6] and ReEvo [9]. 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 \(\varepsilon \to \infty\), 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.

In practice, \(\varepsilon\) is set to a finite intermediate value, providing a good balance of robustness performance. Larger values of \(\varepsilon\) 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 \(\tau\) generations, run the inner loop (adversarial instance search) to search for worst-case instances within the instance-level uncertainty set \(\mathcal{B}_\varepsilon(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 2 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 \(\mathcal{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 7.

After new heuristics are generated, each candidate \(h\) is evaluated on the current adversarial instance set \(\mathcal{S} = \{s'_1, \ldots, s'_K\} \subseteq \mathcal{B}_\varepsilon(S)\) through the robust score \[\hat{r}(h) = \mathcal{A}\bigl(\mathrm{eval}(h, s'_1), \ldots, \mathrm{eval}(h, s'_K)\bigr),\] where \(\mathcal{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 \(\mathcal{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 \(\mathcal{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 \(\tau\) outer-loop generations, RAISE runs an LLM-free evolutionary search to find worst-case instances within \(\mathcal{B}_\varepsilon(S)\) for the current best heuristic \(h \in \mathcal{P}\). Each candidate adversarial instance is parameterized by an 18-dimensional gene vector \(\mathbf{g} \in [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: \[p_{\mathrm{adv}} = \mathrm{clip}\!\left(\bar{p} + \sum_{i=1}^{9} w_i \cdot \phi_i(\mathbf{g}),\; 0, 1\right), \label{eq:adv95dist}\tag{2}\] where \(\bar{p}\) is the mean nominal distribution, \(w_i = g_i - 0.5\) are signed mixture weights, and \(\phi_i\) are nine parametric basis distributions: Uniform, Small, Large, Center, Bimodal, Gaussian, Periodic, Poisson-like, and Peak (Appendix 7.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 \(p_{\mathrm{adv}}\) via 2 , we project it onto the \(\varepsilon\)-ball of the nearest nominal instance to enforce feasibility with respect to \(\mathcal{B}_\varepsilon(S)\): \[p_{\mathrm{proj}} = b + \frac{\varepsilon}{\max(\varepsilon,\, d(p_{\mathrm{adv}}, b))} \cdot (p_{\mathrm{adv}} - b),\] where \(b\) is the closest nominal instance and \(d(\cdot,\cdot)\) is the mean absolute difference defined in Section 2. When \(d(p_{\mathrm{adv}}, b)\) is larger than \(\varepsilon\), this projection places candidates near the boundary of \(\mathcal{B}_\varepsilon(S)\).

Inner evolutionary algorithm. The inner search applies uniform crossover and Gaussian mutation over \(G_{\mathrm{in}}\) generations on a population of \(P_{\mathrm{in}}\) candidate gene vectors, with each vector evaluated by its ability to degrade the current best heuristic \(h\). The worst-scoring feasible instances \(s^*\), minimizing \(\mathrm{eval}(h, s')\) within \(\mathcal{B}_\varepsilon(S)\), is appended to instance set \(\mathcal{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 \(\mathcal{S}\).

3 Experimental Studies↩︎

3.1 Experimental Setup↩︎

3.1.0.1 Tasks:

1) Online Bin Packing (OBP): \(n\) items arrive sequentially and are assigned to bins with a capacity \(C\). Waste ratio \(= (\text{bins used} - \text{lower bound}) / \text{lower bound}\) (lower is better) is used as the objective. Trained on 5 Weibull instances [6], [23] and Tested with \(C \in \{100, 200, 300, 400\}\) and \(n \in \{1000, 5000, 10000\}\). Each size combination has 5 different distributions (Uniform, Normal, Lognormal, Exponential, and Triangular) of items sizes, which results in \(4 \times 3 \times 5 = 60\) datasets for OBP. 2) Online Job Shop Scheduling (JSP): Jobs with 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 \in \{10,20\}\) and \(j \in \{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 \in \{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 8.

3.1.0.2 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 [6], ReEvo [9], and PartEvo [12], each optimised on the nominal training set. Robustness-aware LLM-based AHD: EoH-S [15] (top heuristic trained on 128 diverse instances) and MoH [14] (evaluated with the heuristic reported in the original paper).

3.1.0.3 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 \(\tau = 5\), with inner population \(P_{\mathrm{in}} = 8\) and \(G_{\mathrm{in}} = 4\) inner generations. We ablate four values \(\varepsilon \in \{0.001, 0.002, 0.005, 0.010\}\) on OBP and use \(\varepsilon = 0.002\) for OJSP and OVRP. All hyperparameter and prompt details are in Appendix 7.

4pt

3.2 Results Across Different Distributions↩︎

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

3.2.0.1 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.

3.2.0.2 RAISE achieves consistent OOD robustness.

RAISE\(^{\varepsilon=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\(^{\varepsilon=0.005}\) achieves \(0.273\%\) versus EoH’s \(5.303\%\), a \(19.4\times\) 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.

3.2.0.3 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 = 1\text{k}\) owing to the lower difficulty of small-instance packing, but RAISE decisively surpasses it at \(n \geq 5\text{k}\). 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 (\(\Delta = 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 robustness-aware search generalises to structurally different online combinatorial problems beyond bin packing.

7pt

3.4 Convergence process↩︎

Figure 4 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; (ii) w/o Distribution Constraint removes the \(\varepsilon\)-ball constraint-checking but keep the adversary search, which is similar to [20] and [21]; (iii) w/o Base Distributions replaces the nine-basis parameterisation with flat-vector crossover and Gaussian mutation; (iv) w/o \(\varepsilon\)-Mapping retains the nine-basis parameterisation but skips the boundary projection.

6pt

Table 4 shows that removing any component consistently degrades OOD robustness. Three findings stand out. (i) The \(\varepsilon\)-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\% \to 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 \(\varepsilon\)-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 \(\varepsilon\)-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.

5.0.0.1 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.

Appendix↩︎

Table of Contents↩︎

6 Related Work↩︎

6.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 [6] 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 [9] 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 competitive results with fewer LLM calls. HSEvo [10] proposes a harmony search framework to balance exploitation and exploration. MCTS-AHD [11] employs Monte Carlo tree search to organize heuristics in a tree structure, enabling more effective exploration of the search space. PartEvo [12] incorporates feature-assisted niche construction within abstract search spaces, enabling the seamless integration of niche-based search strategies from evolutionary computation.

6.1.0.1 Generalizable heuristic design

EoH-S [15] 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 [14] 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.

6.1.0.2 Adversarial instance search.

[20] employ an LLM to identify regions of the input space where a heuristic underperforms and to suggest targeted improvements. [21] 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.

6.2 Hyper-Heuristics and Automated Algorithm Configuration↩︎

[24] 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 [25], ParamILS [26], and irace [27] 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.

6.3 Adversarial Training and Minimax Learning↩︎

DRO [22] 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 \(\mathbb{E}_{P_0}[f(x)]\) with the worst-case expected objective \(\sup_{P \in \mathcal{U}} \mathbb{E}_P[f(x)]\) over an ambiguity set \(\mathcal{U}\).

Adversarial training [28] 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 [29] study the theoretical conditions under which a model trained on a source distribution generalizes to a target distribution, providing H-divergence-based 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 [30] apply group DRO to improve worst-case 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.

7 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.

7.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 \(\tau = 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 \(\tau\) 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 \(\mathrm{eval}(h^*_t, s')\) subject to \(d(s', s_i) \leq \varepsilon\) for some \(s_i \in 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 \(\mathcal{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.

7.2 Inner Adversarial Instance Search↩︎

The inner evolutionary algorithm is a compact evolutionary search. Its goal is to find, within the feasible uncertainty set \(\mathcal{B}_\varepsilon(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 \(G_{\mathrm{in}}\) generations using selection, crossover, and mutation. Algorithm 5 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 \(P_{\mathrm{in}} = 8\) gene vectors, each drawn independently from \(\mathrm{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.

Step 2: Decoding and distribution construction (line 4). Each gene vector \(\mathbf{g}_i\) is decoded into a candidate adversarial item-size distribution \(p_{\mathrm{adv}}\) via the basis-mixture model (Eq. 2 in the main paper). The first nine genes \(g_{1:9}\) define signed mixture weights \(w_i = g_i - 0.5 \in [-0.5, 0.5]\) over nine parametric basis distributions, while genes \(g_{10:17}\) control the shape hyperparameters of those bases and \(g_{18}\) adjusts the instance length within \(\pm 10\%\) of the nominal. The resulting \(p_{\mathrm{adv}}\) 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 \(p_{\mathrm{adv}}\) produced by decoding may fall anywhere in distribution space and need not satisfy the feasibility constraint \(d(p_{\mathrm{adv}}, s_j) \leq \varepsilon\) for any \(s_j \in S\). To enforce feasibility, we first identify the nearest nominal instance \(b = \arg\min_{s_j \in S} d(p_{\mathrm{adv}}, s_j)\) under the mean absolute difference \(d\), and then project \(p_{\mathrm{adv}}\) radially onto the \(\varepsilon\)-boundary of the ball centered at \(b\): \[p_{\mathrm{proj}} = b + \frac{\varepsilon}{\max(\varepsilon,\, d(p_{\mathrm{adv}}, b))} \cdot (p_{\mathrm{adv}} - b).\] When \(d(p_{\mathrm{adv}}, b) > \varepsilon\), this formula rescales the displacement vector \((p_{\mathrm{adv}} - b)\) so that the projected point lies exactly on the boundary at distance \(\varepsilon\) from \(b\); when \(d(p_{\mathrm{adv}}, b) \leq \varepsilon\), the point is already feasible and is returned unchanged (since \(\max(\varepsilon, d) = \varepsilon\) 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 \(p_{\mathrm{proj}}\), a concrete problem instance \(s'_i\) is generated. The instance is then passed to the black-box evaluator to compute \(\mathrm{eval}(h^*, s'_i)\), which measures how well \(h^*\) performs on \(s'_i\). The adversarial score is defined as the negation, \(\mathrm{score}_i = -\mathrm{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 \(P_{\mathrm{in}}\) candidates in a generation, the population is ranked by adversarial score in descending order and the top half (i.e., \(\lfloor P_{\mathrm{in}}/2 \rfloor = 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 \((\mathbf{g}_i, \mathbf{g}_j)\) produces one offspring via uniform crossover with probability \(p_{\mathrm{cross}} = 0.5\): for each of the 18 gene positions independently, the offspring inherits the value from \(\mathbf{g}_i\) with probability 0.5 and from \(\mathbf{g}_j\) 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 \(\mathcal{N}(0, 0.12^2)\) with probability \(p_{\mathrm{mut}} = 0.35\), and the result is clipped to \([0,1]^{18}\) to maintain valid gene values. The mutation standard deviation \(\sigma = 0.12\) is chosen to produce moderate perturbations relative to the unit gene range: on average, a mutated gene shifts by approximately \(0.12 \times 0.35 \approx 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 \(G_{\mathrm{in}}\) 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 \(\mathcal{S}\), after which all heuristics in the outer population are re-scored on the updated \(\mathcal{S}\) before evolution continues.

Gene decoding detail. The 18-dimensional gene vector \(\mathbf{g}\) is partitioned as:

• \(g_{1:9}\): basis function weights \(w_i = g_i - 0.5 \in [-0.5, 0.5]\)

• \(g_{10: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)

• \(g_{18}\): length gene (fractional adjustment to item count within \(\pm 10\%\) of nominal)

The nine basis distributions \(\phi_i\) are:

1. Uniform: \(q \sim \mathrm{Uniform}(0,1)\)

2. Small: \(q \sim \mathrm{Beta}(1, p)\) with shape \(p = 1 + 4 \cdot g_{10}\) (favors small items)

3. Large: \(q \sim \mathrm{Beta}(p, 1)\) with shape \(p = 1 + 4 \cdot g_{11}\) (favors large items)

4. Center: \(q \sim \mathrm{Triangular}(0, 0.5, 1)\) (concentrated around 0.5)

5. Bimodal: mixture of \(\mathrm{Uniform}(0, 0.3)\) and \(\mathrm{Uniform}(0.7, 1)\)

6. Gaussian: \(q \sim \mathrm{Normal}(\mu, \sigma^2)\) clipped to \([0,1]\), with \(\mu = g_{12}\), \(\sigma = 0.1 + 0.3 \cdot g_{13}\)

7. Periodic: sinusoidal weight with frequency controlled by \(g_{14}\)

8. Poisson-like: discretized Poisson CDF with rate \(\lambda = 1 + 9 \cdot g_{15}\)

9. Peak: spike at position \(g_{16}\), spreading controlled by \(g_{17}\)

7.3 LLM Evolutionary Operators↩︎

The outer loop employs five LLM operators from [6] with distinct roles, as described in Section 2: Create, 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.

The operators divide into two families. The exploration operators (Create, E1, E2) are designed to introduce structural novelty: Create 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 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.

8 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.

8.1 Benchmark Tasks and Instance Generation↩︎

8.1.1 Online Bin Packing (OBP)↩︎

Problem definition. In Online Bin Packing, items of integer size \(a_i \in \{1, \ldots, C\}\) arrive one at a time from an online sequence \(\sigma = (a_1, a_2, \ldots, a_n)\). 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 \(L = \lceil \sum_i a_i / C \rceil\) the lower bound (optimal offline solution). The waste ratio is: \[\mathrm{WasteRatio} = \frac{B - L}{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 [6], [23], each with \(n = 5{,}000\) items and bin capacity \(C = 200\). These nominal instances constitute the initial adversarial instance set \(\mathcal{S}_0\), 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 (\(\lfloor \cdot \rceil\)), and clipped to the valid range \([1, C]\). Table 6 gives the exact parameterizations used in the generation script.

6pt

Test suite structure. The OBP test suite spans all five distributions, three problem sizes (\(n \in \{1{,}000,\;5{,}000,\;10{,}000\}\)), and four bin capacities (\(C \in \{100,\;200,\;300,\;400\}\)), yielding \(5 \times 3 \times 4 = 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.

6pt

8.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 \(p_{jm} \in \{1, \ldots, 99\}\). The machine order within each job is a random permutation of all \(M\) machines. Jobs carry a release time \(r_j\) (before which they cannot be started) and a due date \(d_j\). 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 \(C_{\max} = \max_j C_j\), where \(C_j\) denotes the completion time of job \(j\).

The normalized makespan divides \(C_{\max}\) by the theoretical lower bound \(\mathrm{LB} = \max\!\bigl(\max_j \sum_m p_{jm},\;\lceil\sum_{j,m} p_{jm} / M\rceil\bigr)\); 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 \(p_{jm}\) 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.

3. Planning horizon. \(H = 2 \times \max(\max_j \sum_m p_{jm},\; \lceil\sum_{j,m} p_{jm} / M\rceil)\).

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.

5. Due dates. \(d_j = r_j + p_j + \delta_j\), where \(p_j = \sum_m p_{jm}\), and \(\delta_j\) combines a multiplicative slack (\(\mathrm{Uniform}(1.35,\;1.75) \times p_j\)), an additive congestion allowance (\(0.35 \times \bar{p}\), with \(\bar{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.

6pt

Test suite structure. The JSP test suite spans all five distributions, two machine counts (\(M \in \{10, 20\}\)), and two job counts (\(J \in \{20, 50\}\)), yielding \(5 \times 2 \times 2 = 20\) dataset files, each with \(I = 5\) instances. The complete test suite covers 100 JSP test instances. Table 9 summarizes the structure.

6pt

8.1.3 Online Vehicle Routing (VRP)↩︎

Problem definition. In Online Vehicle Routing (OVRP), a fleet of \(V\) homogeneous vehicles serves \(N_c = 50\) customers from a central depot. Each customer \(i\) has an integer demand \(d_i \in \{1,\ldots,10\}\) and a location \((x_i, y_i) \in [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 \(N_c = 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 \in \{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 \(d_i\) are drawn from the specified distribution (Table 10), with \(d_i \in [1,\;10]\).

3. Vehicle capacity. Fixed at \(Q = 40\) for all instances with \(N_c = 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\,D_{\max}\) vs.\(0.25\,D_{\max}\) for OBP/JSP) and a reduced scale (\(\sigma_\ell = 0.55\) vs.\(0.60\)), reflecting the narrower integer demand range.

6pt

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

8.2 Baseline Algorithm Descriptions↩︎

8.2.1 Classical Baselines for OBP↩︎

8.2.1.1 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.

8.2.1.2 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 \(\approx 1.7\), but typically performs slightly worse than BestFit because it does not account for remaining capacity when selecting the bin.

8.2.2 Classical Baselines for JSP↩︎

8.2.2.1 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.

8.2.2.2 Minimum Slack (MinSlack).

Dispatching priority is assigned to the job with the smallest slack value \(s_j = d_j - r_j - \sum_m p_{jm}\), where \(d_j\) is the due date and \(r_j\) the release time. MinSlack targets due-date adherence but can produce suboptimal makespan when slack values cluster closely.

8.2.2.3 Earliest Due Date (EDD).

Operations are scheduled by ascending due date \(d_j\). EDD is optimal for minimizing maximum lateness in the single-machine case but offers no direct makespan guarantee on multi-machine instances.

8.2.3 Classical Baselines for VRP↩︎

8.2.3.1 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.

8.2.3.2 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.

8.2.3.3 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.

8.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 (\(N_{\max} = 1{,}000\) samples). Training uses the nominal instance sets described in Section 8.1.

8.2.4.1 EoH [6].

Evolution of Heuristics maintains a population of heuristic Python functions and applies four LLM operator types (i1: create, i2: improve, e1: crossover, e2: 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.

8.2.4.2 MoH [14].

MoH (Meta-objective Heuristics) employs an LLM to iteratively refine a meta-optimizer that autonomously constructs and diversifies heuristic-optimizers via self-invocation, 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.

8.2.4.3 ReEvo [9].

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.

8.2.4.4 PartEvo [12].

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.

8.2.4.5 EoH-S [15].

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.

8.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.

9 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\(^{\varepsilon=0.002}\) throughout unless noted otherwise.

9.1 Online Bin Packing (OBP)↩︎

The metric is shifted waste ratio (%), i.e., \((1 - \text{filled capacity}/\text{bin capacity})\times 100\); lower is better.

9.1.0.1 Per-size breakdown.

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

5.5pt

5.5pt

5.5pt

9.2 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.

4pt

4pt

9.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 \(\dagger\). Tables 17–19 report results for each fleet size separately. Column abbreviations: U=Uniform, N=Normal, LN=Lognormal, Exp=Exponential, Tri=Triangular.

5pt

5pt

5pt

10 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.

10.1 Online Bin Packing↩︎

10.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).

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

    eps = 1e-12
    NEG_INF = -1e9

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

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

    # Feasibility mask
    feasible = bins >= item - 1e-12
    scores = np.full(n, NEG_INF, dtype=float)
    if not np.any(feasible):
        return scores

    b_all = bins.copy()
    # Rank preference: favor packing into smaller bins first (normalized 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)

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

    # Residual after packing
    r = np.maximum(b - item, 0.0)

    # Conservative expected-slot size for future items
    mu = max(0.6 * item, 1e-12)
    n_slots = np.floor(r / mu).astype(int)

    # Geometric-series robustness: diminishing 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)

    # Tightness: how well the item fills the bin now
    tightness = item / (b + eps)

    # Mid-sliver penalty: relative leftover in [0,1]
    frac_left = r / (b + eps)
    mid_sliver_penalty = np.exp(-8.0 * frac_left)

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

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

    combined = (
        w_rank * rank_norm * (1.0 + 0.5 * robustness)
        + w_tight * tightness
        + w_rob * (robustness / (1.0 + robustness))
        + w_pen * mid_sliver_penalty
        + w_exact * exact_fit
    )
    combined = combined - 0.002 * r   # small tie-breaker

    # Normalize to [0, 1] across feasible 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)

    scores[feasible] = norm
    return scores

10.1.2 Robustness Analysis↩︎

Five design choices in Listing [lst:obp] jointly explain why RAISE generalises across Uniform, Normal, Lognormal, Exponential, and Triangular test distributions.

10.1.2.1 (1) Relative rather than absolute metrics.

Tightness is \(\texttt{item}/(b{+}\epsilon)\) 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.

10.1.2.2 (2) Conservative geometric-series robustness term.

The term \(\mu = 0.6\,{\times}\,\texttt{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 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.

10.1.2.3 (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 near-optimal in the average case for any continuous item distribution. Combining it multiplicatively with the robustness term (\(1{+}0.5\,{\times}\,\texttt{robustness}\)) selectively boosts bins that are both tight and capable of accommodating future items.

10.1.2.4 (4) Strong exact-fit bonus.

The \(w_\text{exact}{=}2.5\) bonus saturates the score whenever \(r{\le}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.

10.1.2.5 (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.

10.2 Online Job Shop Scheduling↩︎

10.2.1 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.

import numpy as np
def priority(current_time: int,
             candidate_operations: np.ndarray) -> np.ndarray:
    ops = np.asarray(candidate_operations, dtype=float)
    if ops.size == 0:
        return np.array([], dtype=float)

    duration           = ops[:, 1]
    release_time       = ops[:, 2]
    due_date           = ops[:, 3]
    remaining_work     = ops[:, 4]
    remaining_operations = ops[:, 5]

    start_time      = np.maximum(current_time, release_time)
    completion_time = start_time + duration
    slack           = due_date - completion_time

    # Urgency: small/negative slack gets high priority
    urgency = np.where(
        slack < 0,
        1.0 + (-slack) / (1.0 + duration),   # overdue: grows with tardiness
        1.0 / (1.0 + slack)                   # on-time: diminishing as slack grows
    )

    # SRPT-like finishing: favour small remaining work
    finishing = 1.0 / (1.0 + remaining_work)

    # Downstream risk: more remaining ops => higher priority
    ops_density = remaining_operations / (1.0 + remaining_operations)

    # Duration penalty: penalise long operations unless urgent
    duration_penalty = duration / (1.0 + remaining_work)

    w_urgency = 3.0
    w_finish  = 2.0
    w_ops     = 1.5
    w_dur     = 0.75

    score = (w_urgency * urgency
             + w_finish  * finishing
             + w_ops     * ops_density
             - w_dur     * duration_penalty)

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

10.2.2 Robustness Analysis↩︎

10.2.2.1 (1) Dynamic, state-derived urgency.

Slack \(= \texttt{due\_date} - (\max(\texttt{current\_time},\texttt{release\_time}) + \texttt{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.

10.2.2.2 (2) Asymmetric urgency for overdue vs.on-time jobs.

When slack \({\ge}0\), urgency \(= 1/(1+\texttt{slack})\) decreases gently as slack grows, de-prioritising comfortable jobs. When slack \(<0\), urgency \(= 1 + |\texttt{slack}|/(1+\texttt{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.

10.2.2.3 (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{+}\texttt{remaining\_work})\) approximates SRPT for the multi-machine online setting, providing a distribution-independent backbone.

10.2.2.4 (4) Downstream cascade risk.

\(\texttt{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.

10.2.2.5 (5) Duration penalty modulated by remaining work.

\(\texttt{duration\_penalty} = d/(1{+}\texttt{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.

10.3 Online Vehicle Routing↩︎

10.3.1 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, 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.

import numpy as np
def priority(customer_features: np.ndarray,
             candidate_insertions: np.ndarray) -> np.ndarray:
    import numpy as np

    if candidate_insertions is None or candidate_insertions.size == 0:
        return np.array([], dtype=float)

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

    delta_distance    = cand[:, 2]   # incremental distance cost
    residual_capacity = cand[:, 3]   # remaining capacity after insertion
    current_load      = cand[:, 4]   # route utilisation in [0, 1]
    route_size        = cand[:, 5]   # normalised number of stops
    customer_dist     = cand[:, 6]   # normalised distance to insertion point
    projected_route   = cand[:, 7]   # normalised projected route length

    demand   = float(customer_features[2])
    progress = float(np.clip(customer_features[4], 0.0, 1.0))

    # Adaptive weights: shift from capacity/balance to travel cost as routes fill
    w_travel    = 0.3 + 1.0 * progress       # increases with progress
    w_projected = 0.1 + 0.6 * progress       # increases with progress
    w_prox      = 0.5                         # fixed proximity weight
    w_capacity  = 0.2 + 1.0 * (1.0 - progress) # decreases with progress
    w_balance   = 0.1 + 0.8 * (1.0 - progress) # decreases with progress
    w_size      = 0.4                          # fixed route-size penalty

    # Load-balance score: peaks at 50    balance_score = 1.0 - 2.0 * np.abs(current_load - 0.5)

    # Demand scaling: larger demands amplify capacity importance
    demand_scale = 1.0 + 0.5 * max(0.0, demand)

    score = (
        w_travel    * (-delta_distance)
        + w_projected * (-projected_route)
        + w_prox      * (-customer_dist)
        + w_capacity  * residual_capacity * demand_scale
        + w_balance   * balance_score
        - w_size      * route_size
    )

    return np.asarray(score, dtype=float)

10.3.2 Robustness Analysis↩︎

10.3.2.1 (1) Progress-adaptive weighting.

All six weights in Listing [lst:vrp] are functions of progress \(\in[0,1]\), which measures how far route construction has advanced (fraction of customers already assigned). At \(\texttt{progress}{=}0\): \(w_\text{travel}{=}0.3\), \(w_\text{capacity}{=}1.2\), \(w_\text{balance}{=}0.9\). At \(\texttt{progress}{=}1\): \(w_\text{travel}{=}1.3\), \(w_\text{capacity}{=}0.2\), \(w_\text{balance}{=}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.

10.3.2.2 (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 demand distribution the demand-scaled capacity term \(w_\text{capacity}{\times}\texttt{residual}{\times} \texttt{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.

10.3.2.3 (3) Demand-adaptive capacity scaling.

\(\texttt{demand\_scale} = 1 + 0.5\,{\times}\,\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.

10.3.2.4 (4) Balance score peaks at 50% utilisation.

The term \(1 - 2|\texttt{current\_load} - 0.5|\) equals 1 at perfect half-load and \(-1\) at completely empty or fully loaded routes. The 50% target is distribution-agnostic: 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 distribution-specific optimal fill level degrade whenever the actual distribution shifts.

10.3.2.5 (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.

11 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).

11.1 RAISE Cost Breakdown↩︎

Table 20 lists the cost of a single RAISE run at the standard budget of \(N_{\max} = 1000\) LLM samples, using a GPT-5-mini class model used in our experiments.

11.2 Comparison with Existing Frameworks↩︎

Table 21 compares RAISE against EoH, ReEvo, PartEvo, MoH, and EoH-S at \(N_{\max}{=}1000\) on online bin packing problem.

5pt

11.2.0.1 Token overhead.

All methods share the same 1,000-call budget. RAISE’s total token cost (\(\approx\)​2.82M) is comparable to EoH (\(\approx\)​2.80M). ReEvo and PartEvo are substantially more expensive in tokens (\(\approx\)​4.30M and \(\approx\)​4.51M respectively) due to their multi-stage reflection and critique prompting strategies, which usually generate longer and more elaborate heuristic programs.

11.2.0.2 EoH-S incurs a hidden evaluation overhead.

EoH-S has a token cost comparable to RAISE, but its instance evaluation count (\(\approx\)​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.

11.2.0.3 MoH incurs dual-task evaluation overhead.

MoH’s \(\approx\)​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.

11.2.0.4 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.

References↩︎