Incorporating Large Language Models into the Bayesian Optimization Process for Engineering Design

Bayesian optimization (BO) is a framework typically applied to optimizing expensive black-box functions with minimal evaluations (Jones et al. 1998; Shahriari et al. 2016; Frazier 2018). It maintains a surrogate model approximating the function of interest, typically via a Gaussian process (GP), and uses an acquisition function to iteratively select points for evaluation. This approach is widely applied within engineering design cycles, where each function evaluation may require costly simulations or physical experiments.

In practice, engineering optimization rarely starts from zero knowledge. Domain experts provide prior knowledge about the design space, physical constraints that prevent infeasible configurations, and information about parameter interactions. This expertise is traditionally incorporated through manual design-of-experiments, hand-tuned search bounds, or custom kernel functions (Snoek et al. 2012; Shahriari et al. 2016), all of which require human effort and close coordination between domain experts and the optimization workflow. This coordination is often costly, as domain experts may not be available at each iteration, formalizing prior knowledge into mathematical constraints is labor-intensive, and the process must be repeated for each new problem.

Large language models (LLMs) present an opportunity to bridge this gap. Motivated by human-in-the-loop optimization platforms such as CAMEO (Meyur et al. 2025; Purohit et al. 2025), where domain experts guide automated experimental campaigns, we investigate whether LLMs can use natural language domain knowledge to generate quantitative signals useful for BO. As LLMs are trained on a broad corpus of information, including scientific literature, engineering textbooks, and technical discussions, they are able to reason over generalized domain knowledge (Phan et al. 2025). LLMs have consequently been used as both candidate generators and surrogates within BO for hyperparameter tuning (T. Liu et al. 2024), making them promising candidates for this task.

In this paper, we investigate how LLMs can incorporate natural language domain knowledge into the BO process for optimizing NACA four-digit airfoil profiles for maximum lift-to-drag ratio. We study two injection points: (1) informed initialization (Section 3.1), where LLMs replace standard space-filling schemes with domain-aware point selection, and (2) informed surrogates (Section 3.2), where LLMs are used to augment a GP surrogate with few-shot in-context predictions that leverage domain knowledge.

Our contributions are threefold. First, we compare LLM-based initialization strategies, including multiple agentic reasoning topologies, against standard sampling methods. Second, we evaluate three LLM augmented surrogate strategies that vary the degree to which LLMs influence the BO loop. Third, we use controlled comparisons to assess how domain context and model choice affect optimization performance. The evaluation is limited to one discrete NACA airfoil case study; therefore, transfer to other engineering domains, higher-dimensional parameterizations, and continuous design spaces remains to be demonstrated.

The remainder of this paper is organized as follows: Section 2 reviews relevant background and related work, Section 3 describes our LLM-augmented BO framework, Section 4 presents experimental results, and Section 5 discusses implications and limitations.

2.1 Related Work

The use of prior knowledge to warm-start BO has been studied through transfer learning and meta-learning (Feurer et al. 2018), where surrogate models are pre-trained on related tasks. LLMs offer a complementary modality requiring only natural language specification of the domain rather than structured historical data. More broadly, LLMs have been applied to molecular design (Jablonka et al. 2024) and parametric shape optimization for airfoils (Zhang et al. 2024), though these efforts adopt evolutionary strategies rather than integrating the LLM into an iterative BO loop.

LLAMBO (T. Liu et al. 2024) was among the first to demonstrate that, via few-shot prompting over past observations, LLMs can serve as effective surrogates and candidate generators within Bayesian optimization for machine learning hyperparameter tuning. Subsequent extensions include LLINBO (Chang et al. 2025), which adds formal regret bounds and trust-controlled LLM proposals, and BOLT (Zeng et al. 2025), which trains LLMs on BO trajectories for multi-task transfer. Beyond hyperparameter tuning, LLM-guided BO has been applied to database configuration (Lao et al. 2024) and analog circuit design (Chen et al. 2024). Our work extends this direction to engineering design, where physics-grounded domain knowledge plays a prominent role. We additionally evaluate a cyclic surrogate strategy and a running hypothesis mechanism that, to the best of our knowledge, have not been reported in prior LLM-BO work. The proposed agentic initialization topologies (Diamond, Tournament, Debate) draw on work in multi-agent LLM collaboration (Z. Liu et al. 2024), where specialist agents with complementary roles are coordinated to improve LLM-based reasoning.

2.2 Bayesian Optimization

Bayesian optimization aims to identify 𝒙^∗ that maximizes a function 𝑓 over a design space 𝒳, where . The process alternates between two steps: (1) fitting a probabilistic surrogate model to the observed data and (2) selecting the next evaluation point by maximizing an acquisition function . A commonly used surrogate is a Gaussian process (GP) (Rasmussen and Williams 2006), which provides both a mean prediction 𝜇(𝒙) and an uncertainty estimate 𝜎(𝒙) at any point in the design space. In this work, we use a GP surrogate with Expected Improvement (EI) (Jones et al. 1998) as the acquisition function. EI balances exploitation of regions with high predicted mean against exploration of regions with high uncertainty. The overall procedure is summarized in Algorithm 1.

Algorithm 1: Bayesian optimization (BO) overview

The BO loop requires an initial set of observations 𝒟₀ (line 2) in order to fit the first surrogate. Standard approaches use random sampling, Latin hypercube sampling (LHS), or Sobol sequences to provide space-filling coverage of 𝒳. These schemes are agnostic to any prior knowledge about the system.

2.3 Problem Domain: NACA Airfoil Optimization

The NACA four-digit airfoil family (Jacobs et al. 1933) parameterizes airfoil geometry with three shape parameters, all expressed relative to the chord (the distance from the leading to the trailing edge): maximum camber 𝑚 (first digit, as % of chord), location of maximum camber 𝑝 (second digit, in tenths of chord), and maximum thickness 𝑡 (third and fourth digits, as % of chord). Together with the angle of attack 𝛼, these four parameters form the design vector 𝒙 = (𝑚, 𝑝, 𝑡, 𝛼) we optimize over.

We evaluate performance using the lift-to-drag ratio (𝐶_L/𝐶_D), a standard aerodynamic efficiency metric that measures how much lift an airfoil generates per unit of drag, and is a primary driver of range, endurance, and fuel burn in fixed-wing aircraft design. 𝐶_L/𝐶_D is computed via the XFOIL panel-method solver (Drela 1989) at Reynolds number 𝑅𝑒 = 10⁶ with critical amplification ratio 𝑁_crit = 5. We construct a discrete design pool over 𝒙 as defined in Table 1. The full factorial enumeration yields 3,234 combinations; after filtering for XFOIL convergence failures, 3,206 feasible designs remain with a true optimum of 𝐶_L/𝐶_D = 147.6 at NACA 6409, 𝛼 = 5.0∘. Figure 1 illustrates this airfoil with its shape parameters annotated. The design space exhibits non-trivial parameter interactions, with multiple local optima (Figure 2).

Figure 1: NACA 6409, the optimal lift-to-drag ratio airfoil in the design pool, with the dashed camber line and
its three shape-defining parameters annotated

Table 1: Design space parameters for the NACA four-digit airfoil optimization problem.

Figure 2: Objective landscape (𝐶_L/𝐶_D) as a function of thickness and angle of attack. (a) At the globally optimal camber (𝑚 = 6, 𝑝 = 4), (b) At a local optimum (𝑚 = 4, 𝑝 = 3), (c) Maximum 𝐶_L/𝐶_D over all (𝑚, 𝑝) at each (𝑡, 𝛼). The optimal configuration is denoted by a red star.

We investigate two points where domain knowledge can enter the BO loop. Both operate within the same discrete BO framework: at each iteration, candidates are drawn from the unobserved portion of the design space, scored by the surrogate(s), ranked by the acquisition function, and the top candidate is evaluated.

3.1 LLM-Informed Initialization

In standard BO, the initial observation set 𝒟₀ is typically generated by a space-filling scheme (random, LHS, or Sobol). Following the initialization approach introduced by LLAMBO (T. Liu et al. 2024) for hyperparameter tuning, we adapt LLM-based candidate generation to aerodynamic design, replacing the space-filling scheme with LLM-generated initial points that leverage domain knowledge. All informed LLM calls in this work, both initialization and surrogate, are prefixed with a shared domain context block Ω (Figure 3) that specifies the task, objective metric, simulator, parameter definitions drawn from the NACA four-digit convention (Abbott and Doenhoff 1959), and a physical-reasoning guide grounded in standard low-speed aerodynamics references (Anderson 2017; Abbott and Doenhoff 1959). The initialization prompt (Figure 4) prepends Ω and requests 𝑛_init = 5 diverse initial designs suitable for BO exploration. The physical-reasoning portion of Ω covers:

• Parameter-level mechanisms: how each variable generally influences aerodynamic performance (higher camber raises both lift and drag; thicker sections delay stall but add form drag; camber position shifts the pressure peak).
• Coupling and non-separability: the best value of any one parameter depends strongly on
  the others, so independent reasoning breaks down.
• Low-Re regime caveats: at 𝑅𝑒 = 10⁶, 𝑁_crit = 5, textbook heuristics can invert; e.g., higher
  camber can still raise L/D despite added pressure drag.

The LLM returns parameter combinations which are mapped to the nearest valid design in the pool. Initialization is a zero-shot task, with the LLM receiving only Ω and no prior evaluation data. The parameter definitions and physical-reasoning guide are listed in full in Appendix A; both are stripped in the uninformed ablation (Section 3.3) to isolate the effect of domain knowledge.

Figure 3: Domain context prepended to every LLM call, inclusive of: task, metric, and simulator descriptions, parameter definitions, and the physical-reasoning guide (Appendix A).

Figure 4: Initialization prompt card: red domain context (Figure 3) followed by task instructions.

3.1.1 Agentic Topologies

We evaluate four agentic topologies, each expressed as a directed acyclic graph (DAG) of LLM calls, for generating initial points (Figure 5). We use agentic in the sense of role-differentiated LLM calls communicating across a DAG, as in DyLAN (Z. Liu et al. 2024) and related multi-agent frameworks. The topology determines what each node receives as input, with all nodes receiving the initialization prompt (Figure 4) and non-root nodes additionally receiving their parents’ proposed points. Each node’s role is assigned based on its position within the DAG by a single LLM planning call per topology, following the team-optimization concept of Z. Liu et al. (2024), and prepended to that node’s prompt.

The four topologies are: the Single topology, a one-node DAG that issues a single LLM call and returns its proposal directly; Diamond, which splits the task across two proposer agents whose outputs are merged into a single output; Debate, in which three proposers generate candidates and iteratively refine them over two rounds before being merged; and Tournament, which brackets four proposer agents through two pairwise merging rounds.

Figure 5: The four agentic topologies evaluated for LLM-informed initialization.

3.2 LLM-Informed Surrogate

Beyond initialization, LLMs can inform the iterative BO loop itself by augmenting existing surrogate models. We evaluate three strategies that vary the degree to which an LLM augments a GP surrogate.

3.2.1 Surrogate Iteration

All three surrogate strategies follow the common template in Algorithm 2. Each iteration:

• Updates the GP posterior, the running LLM generated natural-language hypothesis (𝐻_t), and the current LLM defined region of interest (𝑅_t).
• Constructs a shortlist of candidates selected by GP expected improvement and/or randomly sampled from the current LLM defined region using the budgets in Table 2.
• Ranks the shortlist using blended GP-LLM predictions.

The strategies differ only in the composition of the shortlist, and apply the common template as follows:

• GP-Led: All shortlisted candidates are selected based on expected improvement from the GP over the unobserved pool.
• LLM-Led: The shortlist combines equal numbers of candidates selected by GP expected improvement and candidates sampled from the current LLM-defined region of interest.
• Cycle-k: The LLM-led shortlist strategy is retained; however, every k-th BO iteration restricts the shortlist to candidates sampled only from the current LLM-defined region of interest.

Algorithm 2 incorporates two distinct LLM roles: LLM_update jointly updates the hypothesis and region (Figure 6), while LLM_eval predicts the objective values of shortlisted designs (Figure 7). Both are prefixed with the same domain context (Ω) used at initialization (Figure 3), so that LLM contributions to the surrogate loop combine domain context with the observations 𝒟_t.

Algorithm 2: LLM-Augmented BO Iteration, all experiments use 𝜔=0.5

Table 2: Candidate budgets used by each surrogate strategy. The reported experiments use K=16

3.2.2 Hypothesis Mechanism

The LLM maintains a running hypothesis (𝐻_t), a natural-language summary of which parameter regions are believed to produce high 𝐶_L/𝐶_D, together with a region of interest (𝑅_t) defined by bounds on the design variables and the permitted angle values. At each iteration, the LLM is prompted with the domain context (Ω), the full observation history (𝒟_t), along with the previous hypothesis (𝐻_t-1 ), and region (𝑅_t-1 ), to produce updated versions reflecting the observed trends (Figure 6). The region is then used to filter the discrete design pool, from which LLM-region candidates are sampled. A representative hypothesis trajectory is shown in Appendix C.

Figure 6: Joint hypothesis and region update prompt card (𝐿𝐿𝑀_update)

3.2.3 LLM Prediction and Prompt-Order Variation

Following T. Liu et al. (2024), we obtain LLM predictions by prompting the model with the observation history as few-shot examples and requesting a predicted objective value for a query design (Figure 7). To approximate uncertainty, we exploit LLM sensitivity to example ordering: for each query point 𝒙 we perform three calls with uniquely shuffled observation orderings and calculate 𝜇_LLM(𝒙) and 𝜎_LLM(𝒙) as the sample mean and standard deviation of the predictions. This standard deviation provides an empirical proxy for LLM predictive uncertainty and is blended with the GP uncertainty, but it should not be interpreted as calibrated or converged predictive uncertainty.

Figure 7: LLM evaluation prompt card for predicting shortlisted objective values (𝐿𝐿𝑀_eval)

3.3 Experimental Setup

All BO runs use 𝑛_init = 5 initial points and 𝑛_iter = 10 iterations with a Matérn-5/2 GP (Rasmussen and Williams 2006) and expected improvement (EI) acquisition function. Claude Sonnet 4 (Anthropic 2025b) is used as the baseline LLM, with Claude Opus 4.6 (Anthropic 2025a) used to assess the effect of increased model capability during informed initialization. Capability metrics for both models are reported in Appendix D. For each initialization method, 50 independent sets of initial points are generated, each used as the warm start for 20 independent BO runs (varying only in GP hyperparameter optimization stochasticity), yielding 1,000 total runs per method. All surrogate methods share the same random initial points (seeded for reproducibility) across 30 independent repeats, enabling paired statistical comparisons.

Both experiments include an ablation where the domain context (Ω, Figure 3) is stripped, and is referred to as uninformed. The parameter names (𝑚, 𝑝, 𝑡, 𝛼) are replaced by (𝑋1, 𝑋2, 𝑋3, 𝑋4), the objective (𝐶_L/𝐶_D) is replaced by (𝑌), and references to NACA airfoils, aerodynamics, XFOIL, Reynolds number, 𝑁crit, and aerodynamic expertise are removed. In the surrogate experiment, the LLM receives anonymous input-output observations accumulated during BO. This condition tests whether the LLM-augmented process provides value without explicit domain knowledge, when applied to the same underlying optimization landscape.

The experiments are structured to compare the main design choices of the proposed frameworks. For initialization, we compare the effects of domain context, model capability, and agent topology. For the modified surrogate, we compare the overall LLM-augmented process with standard GP-only BO and examine how changing the candidate-selection strategy and domain context affects performance. Appendix B summarizes the interpretation of each comparison.

We report a range-normalized score 𝑆 that places performance on a common [0,100] scale, with 𝑓_min = 57.95 and 𝑓_max = 147.61 the worst and best feasible 𝐶_L/𝐶_D values in the design pool. Each BO iteration produces a set 𝒟 of evaluated designs, with the quantity of interest being the best design found, .

(2)

4.1 Informed Initialization Results

Table 3 summarizes initialization results for the three statistical baselines (Random, LHS, Sobol) and six LLM-guided variants, spanning four agentic topologies and two LLM models. The Informed column marks whether the full domain context is used or stripped of the parameter definitions and physical-reasoning guide. We report the mean and standard deviation of 𝑆^* at initialization (the best of the 5 initial points at iteration 0) and after 10 BO iterations, aggregated across 50 initialization rounds each with 20 BO reruns. Opus 4.6 achieves the highest observed mean final 𝑆^⋆ of 89.3, followed by Sonnet 4 at 85.2. All three baselines (Random, LHS, Sobol) cluster around 80 after 10 BO iterations.

Table 3: Informed initialization results.

Figure 8 shows the convergence curves (left) and final iteration 𝑆^* distributions (right) for all methods. The informed Single configurations generally maintain higher mean performance than the three non-LLM space-filling baselines throughout BO, while the multi-agent configurations show more variable performance. The uninformed Single configuration starts from a substantially lower initial score and remains below the informed configurations. Formal statistical comparisons are reported in Section 4.3.

Analysis of parameter selection frequencies (Figure 9) reveals that LLM-based methods develop non-uniform preferences over the design space. Informed LLMs under the single configuration are biased towards thin airfoils (𝑡 < 13), with the stronger Opus 4.6 model increasing this preference in comparison to Sonnet 4. All multi-agent topologies (Diamond, Tournament, Debate) result in noticeable clusters forming at the boundaries of 𝑡, 𝑚, and 𝑝.

Figure 8: Informed Initialization. Left: convergence curves; each line shows the mean of S^* across n = 1,000 runs at every BO iteration. Right: distribution of per-run S^* at the final iteration.

Figure 9: Parameter selection frequency for initial points (n = 250 points per method)

4.2 Modified Surrogate Results

Table 4 summarizes surrogate experiments for the GP-only baseline and LLM-augmented blending strategies. The Informed column indicates whether the full domain context is provided to LLM calls. We report the mean and standard deviation of the final 𝑆^* after 10 BO iterations across 𝑛 = 30 paired repeats, two reliability readouts, denoting 𝑃(𝑆^* ≥ 85) and 𝑃(𝑆^* ≥ 90) , and the mean final 𝑆^* based on the quality of the initial points (low, mid, high; 𝑛 = 10 per tertile bin). All methods share the same random initial points with a mean 𝑆^* = 65.3.

The LLM-augmented configurations show higher observed final means than GP-only, with informed Cycle-3 achieving the highest observed mean (89.1). Removing explicit aerodynamic and NACA context reduces the Cycle-3 mean to 84.9, but the uninformed configuration remains above GP-only (78.8). The largest separation occurs for low-quality initial sets, where informed Cycle-3 reaches 84.6 compared with 67.2 for GP-only. These subgroup values are descriptive as each bin contains only n=10 repeats. Figure 10 shows the convergence curves (left) and final iteration 𝑆^* distributions (right) for all surrogate strategies. All LLM-augmented methods achieve higher observed means than GP-only.

Table 4: Modified surrogate results. Reliability columns give the fraction of the 𝑛 = 30 repeats that clear the 𝑆^* threshold. The low/mid/high columns report the mean final 𝑆^* conditional on the quality of the initial set.

Figure 10: Modified Surrogate. Left: convergence curves; each line shows the mean of S^* across n = 30 paired repeats at every BO iteration. Right: distribution of per-run S^* at the final iteration.

To examine how the methods traverse the design space, Figure 11 reports the normalized distance between consecutive BO selections. Each design parameter is normalized to [0,1], and the step size is defined as the Euclidean distance between consecutive selected designs. The LLM-augmented methods take smaller steps on average than GP-only, consistent with dampening the long-range jumps characteristic of GP-EI and promoting more localized movement through the design space.

Figure 11: Distribution of step sizes between consecutive BO selections, by surrogate method.

4.3 Statistical Significance

We conduct two statistical tests that assess distinct questions. Paired Wilcoxon signed-rank tests evaluate whether final 𝑆^* values differ systematically between two methods, while Levene tests evaluate whether their final-score spreads differ. We apply the Benjamini-Hochberg correction separately across the eleven Wilcoxon and eleven Levene tests. The complete set of raw and adjusted values is reported in Table 5, with adjusted p-values below 0.05 denoted in bold.

Table 5: Pairwise tests for key comparisons. ΔS^* is the mean paired difference; 𝑝_Wilcoxon is the raw Wilcoxon p-value, 𝑝_{Wilcoxon (BH)} and 𝑝_Levene(BH) are global Benjamini-Hochberg-adjusted Wilcoxon and Levene p-values, respectively. Wilcoxon and Levene tests are corrected as separate families. Experiment 1 pairs by warm start seed (n = 50); Experiment 2 pairs by shared initial points (n = 30).

4.4 Computational Cost

The LLM-augmented methods introduce substantial inference overhead relative to GP-only. Across the reported 10-iteration surrogate configurations, each run used upwards of 490 LLM calls at a cost of $1.96-$2.34. Initialization required approximately $0.02-$0.18 per generated warm start. Across all reported configurations and repeats, the aggregate LLM cost was $353.20. Timing validation gave median runtimes of approximately 323 seconds for Cycle-3 and 0.84 seconds for GP-only. Detailed token counts and rates are provided in Appendix D. We note that minimal effort was placed into optimizing for runtime, and large portions of the proposed LLM modified surrogate algorithms could be run in parallel.

In this work, we examined two ways of incorporating LLM guidance into BO. First, we investigated whether LLM-informed initialization can provide useful starting designs compared with standard space-filling methods. Second, we assessed whether the proposed LLM-augmented surrogate strategies can improve iterative optimization relative to the GP-only implementation. Here we consider the initialization and surrogate results separately, before addressing attribution and generalization.

Table 3 shows that domain context and model capability can provide useful guidance when selecting starting designs. Claude Opus 4.6 produces the highest observed initialization and final mean 𝑆^*, with its improvement over Random remaining significant after Benjamini-Hochberg correction (p=0.0003). Opus 4.6 also produces a narrower final-score distribution than Random, with the corresponding Levene comparison remaining significant after correction (p=0.0034). This spread difference is not observed consistently across the other comparisons. While the informed Claude Sonnet 4 configuration also has a higher observed mean compared to the space-filling baseline, this comparison does not remain significant after correction (p=0.0964). Although Opus 4.6 achieves a higher observed initialization mean than Sonnet 4, their paired comparison does not remain significant after correction (p=0.0537). Removing the domain context from the Single configuration substantially reduces performance (p=0.0006), indicating that the supplied context can provide useful guidance at initialization.

Multi-agent topologies (Diamond, Tournament, Debate) underperform the Single topology across the reported initialization metrics. Inspecting the proposed initial points, we observe that the topologies involving multiple rounds of message passing tend to cluster their selections at the boundaries of the design space (Tournament concentrating on low camber and forward position, Debate clustering at extreme values, as described in Figure 9). We attribute this to automatic role assignments offloading a complex decomposition task onto the LLM itself. Without curated roles, the agents appear to anchor on extreme exemplars rather than converging on the region of interest. We hypothesize that curated roles that better decompose the four-parameter problem, or higher-dimensional and multi-objective problems where decomposition is more natural, would give multi-agent approaches more room to contribute, though this remains to be demonstrated.

For the surrogate experiment, Table 4 shows that the LLM augmented methods have higher observed final performance than the GP-only implementation, including when aerodynamic context is removed. Informed Cycle-3 improves upon GP-only after Benjamini-Hochberg correction (p=0.0049), while the uninformed Cycle-3 variant also improves upon the GP-only baseline (p=0.0444). This indicates that the cyclic mechanism retains value without explicit aerodynamic or NACA context. The informed Cycle-3 configuration has a 4.2-point higher observed mean than the uninformed variant, but the direct comparison falls slightly above the significance threshold after correction (p=0.0537). The GP-led and LLM-led configurations also have higher observed means than GP-only, but their comparisons do not remain significant after BH correction.

One observable difference between the methods is how they traverse the design space. Figure 11 shows that GP-only takes larger normalized steps than the LLM-augmented methods, consistent with more localized movement through the design space. This pattern may help explain how the LLM-augmented methods more readily exploit promising regions than GP based BO, but this does not establish that smaller steps cause the observed performance gains. Table 4 shows that the methods separate most strongly when the initial set is poor, with GP-only achieving a mean final 𝑆^* of 67.2 in the lowest initialization bin, while informed Cycle-3 reaches 84.6. Across all initializations, the cyclic configurations achieve the highest observed means and reliability, with informed Cycle-3 reaching a mean final 𝑆^* of 89.1. However, Cycle-1 and Cycle-3 are not statistically distinguishable (p=0.5034). No difference in spread among the surrogate configurations remains significant after Benjamini-Hochberg correction, as shown by the adjusted Levene tests in Table 5.

These results are specific to the tested airfoil problem and do not establish whether the benefit arises from general aerodynamic knowledge acquired during training or possible exposure to related benchmark data during LLM training. Recall of the optimum NACA airfoil may be less likely because the problem combines a particular XFOIL configuration, including the selected Reynolds number, critical amplification ratio, angle-of-attack range, and discrete parameter bounds. Furthermore, the initialization results do not show consistent identification of the optimum from the prompt alone, with LLM-generated starting sets varying considerably. This pattern is more consistent with a useful yet imperfect prior as opposed to reliable recall of a specific optimum. We note however, that this interpretation cannot be verified because the training data for Claude Opus 4.6 and Claude Sonnet 4 is not fully observable.

Future work could vary the source and structure of the injected information, including literature, simulator documentation, design-history data, or expert-curated constraints. The natural-language hypothesis also provides a possible interface for future human-in-the-loop studies in which experts inspect or edit the LLM’s working hypothesis during optimization. We evaluate only a single discrete NACA airfoil problem, so generalization to continuous spaces, higher-dimensional problems, and other engineering domains remains to be demonstrated and should be prioritized. Future studies should also conduct component-level ablations of hypothesis updating, prompt-order repetition, blending weights, and shortlist size to determine how these mechanism choices contribute independently to the observed performance gains.

We investigated two modes in which LLMs can support Bayesian optimization: LLM-guided initialization, where space-filling schemes are replaced with domain-aware initial points, and LLM-augmented surrogates, which blend a Gaussian process with few-shot LLM predictions. The evaluation used lift-to-drag ratio maximization on a discrete NACA four-digit airfoil benchmark.

The LLM-guided initialization strategy using Claude Opus 4.6 improves final 𝑆^* by 8.2 points over a randomized space filling scheme after Benjamini-Hochberg correction (p<0.001). This result shows that LLM guidance can improve the selection of starting designs for BO on the specific NACA airfoil problem. The LLM-augmented surrogate strategy, Cycle-3, using Claude Sonnet 4 improves final 𝑆^* by 10.2 points over the standard GP-only implementation (p=0.0049). The same surrogate strategy without explicit aerodynamic or NACA context improves by 6.1 points (p=0.0444), indicating that the proposed LLM based surrogate retains value without explicit domain knowledge on this benchmark. The result supports the complete cyclic strategy rather than attributing its improvement solely to explicit aerodynamic context.

Together, these findings show that LLM guidance can contribute both to BO initialization and to iterative surrogate-guided optimization, with the natural-language hypothesis providing a possible interface for future expert inspection or intervention. Within the scope of this discrete NACA airfoil benchmark, the results support LLMs as a flexible source of guidance for both BO initialization and iterative optimization, motivating evaluation on higher-dimensional, continuous, and other engineering problems.

Abbott, Ira H., and Albert E. von Doenhoff. 1959. Theory of Wing Sections: Including a Summary of Airfoil Data. Dover Publications.

Anderson, John D., Jr. 2017. Fundamentals of Aerodynamics. 6th ed. McGraw-Hill Education.

Anthropic. 2025a. Claude Opus 4.6 System Card. Anthropic.

Anthropic. 2025b. Claude Sonnet 4 System Card. Anthropic.

Artificial Analysis. 2026. “Artificial Analysis Intelligence Index.” Accessed July 2026. https://artificialanalysis.ai/.

Chang, Chih-Yu, Milad Azvar, Chinedum Okwudire, and Raed Al Kontar. 2025. “LLINBO: Trustworthy LLM-in-the-Loop Bayesian Optimization.” arXiv preprint arXiv:2505.14756.

Chen, Guojin, Keren Zhu, Seunggeun Kim, Hanqing Zhu, Yao Lai, Bei Yu, and David Z. Pan. 2024. “LLM-Enhanced Bayesian Optimization for Efficient Analog Layout Constraint Generation.” arXiv preprint arXiv:2406.05250.

Drela, Mark. 1989. “XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils.” In Low Reynolds Number Aerodynamics, edited by T. J. Mueller, vol. 54. Lecture Notes in Engineering. Springer. https://doi.org/10.1007/978-3-642-84010-4_1.

Feurer, Matthias, Benjamin Letham, Frank Hutter, and Eytan Bakshy. 2018. “Practical Transfer Learning for Bayesian Optimization.” ICML 2018 AutoML Workshop. https://arxiv.org/abs/1802.02219.

Frazier, Peter I. 2018. “A Tutorial on Bayesian Optimization.” arXiv preprint arXiv:1807.02811.

Jablonka, Kevin Maik, Philippe Schwaller, Andres Ortega-Guerrero, and Berend Smit. 2024. “Leveraging Large Language Models for Predictive Chemistry.” Nature Machine Intelligence 6 (2): 161-69. https://doi.org/10.1038/s42256-023-00788-1.

Jacobs, Eastman N., Kenneth E. Ward, and Robert M. Pinkerton. 1933. The Characteristics of 78 Related Airfoil Sections from Tests in the Variable-Density Wind Tunnel. NACA-TR-460. National Advisory Committee for Aeronautics (NACA).

Jones, Donald R., Matthias Schonlau, and William J. Welch. 1998. “Efficient Global Optimization of Expensive Black-Box Functions.” Journal of Global Optimization 13 (4): 455-92. https://doi.org/10.1023/A:1008306431147.

Lao, Jiale, Yibo Wang, Yufei Li, Jianping Wang, Yunjia Zhang, Zhiyuan Cheng, Wanghu Chen, Mingjie Tang, and Jianguo Wang. 2024. “GPTuner: A Manual-Reading Database Tuning System via GPT-Guided Bayesian Optimization.” Proceedings of the VLDB Endowment 17 (8): 1939-52. https://doi.org/10.14778/3659437.3659449.

Liu, Tennison, Nicolás Astorga, Nabeel Seedat, and Mihaela van der Schaar. 2024. “Large Language Models to Enhance Bayesian Optimization.” International Conference on Learning Representations (ICLR). https://arxiv.org/abs/2402.03921.

Liu, Zijun, Yanzhe Zhang, Peng Li, Yang Liu, and Diyi Yang. 2024. “A Dynamic LLM-Powered Agent Network for Task-Oriented Agent Collaboration.” Conference on Language Modeling (COLM). https://arxiv.org/abs/2310.02170.

Meyur, Rounak, Sam Donald, Tonya Martin, Thiagarajan Ramachandran, and Sumit Purohit. 2025. “CAMEO: A Co-Design Architecture for Multi-Objective Energy System Optimization.” Proceedings of the SC ’25 Workshops of the International Conference for High Performance Computing, Networking, Storage and Analysis (New York, NY, USA), SC workshops ’25, 2201-12. https://doi.org/10.1145/3731599.3767576.

Phan, Long, Alice Gatti, Ziwen Han, Nathaniel Li, Josephina Hu, et al. 2025. “Humanity’s Last Exam.” arXiv preprint arXiv:2501.14249.

Purohit, Sumit, Rounak Meyur, Sam A. Donald, Tonya J. Martin, and Thiagarajan Ramachandran. 2025. CAMEO: A Co-Design Architecture for Multi-Objective Energy System Optimization (Project Report). PNNL-SA-3015290. Pacific Northwest National Laboratory (PNNL), Richland, WA. https://doi.org/10.2172/3015290.

Rasmussen, Carl Edward, and Christopher K. I. Williams. 2006. Gaussian Processes for Machine Learning. MIT Press.

Shahriari, Bobak, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando de Freitas. 2016. “Taking the Human Out of the Loop: A Review of Bayesian Optimization.” Proceedings of the IEEE 104 (1): 148-75. https://doi.org/10.1109/JPROC.2015.2494218.

Snoek, Jasper, Hugo Larochelle, and Ryan P. Adams. 2012. “Practical Bayesian Optimization of Machine Learning Algorithms.” Advances in Neural Information Processing Systems 25 (NIPS 2012), 2960-68.

Zeng, Yimeng, Natalie Maus, Haydn Thomas Jones, Jeffrey Tao, Fangping Wan, Marcelo Der Torossian Torres, Cesar de la Fuente-Nunez, Ryan Marcus, Osbert Bastani, and Jacob R. Gardner. 2025. “Large Scale Multi-Task Bayesian Optimization with Large Language Models.” arXiv preprint arXiv:2503.08131.

Zhang, Xinxin, Zhuoqun Xu, Guangpu Zhu, Chien Ming Jonathan Tay, Yongdong Cui, Boo Cheong Khoo, and Lailai Zhu. 2024. “Using Large Language Models for Parametric Shape Optimization.” arXiv preprint arXiv:2412.08072.

Appendix A: Parameter Definitions and Domain Context

The two prompt blocks that carry the domain knowledge injected at initialization, parameter definitions (the model card) and the domain-context block, are reproduced verbatim below as they appear in the prompt for the 4D C_L/C_D objective. Their content is grounded in the NACA four-digit convention (Abbott and Doenhoff 1959) and standard references on low-speed airfoil behaviour (Anderson 2017). In the uninformed ablation (3.3), both blocks are stripped, leaving only the parameter ranges and the task instruction.

Parameter Definitions:

Each design is a NACA MPTT airfoil defined by four parameters (4D optimisation):

m (max camber, % of chord): 1 to 6 (valid: [1, 2, 3, 4, 5, 6])

– m=0 means symmetric (zero camber) — no lift at zero angle of attack
– Higher camber generates more lift (𝐶_L) at any given angle of attack
– Higher camber also increases pressure drag (CD)
– The NET effect on C_L/C_D depends on how much C_L gains outweigh C_D gains, which varies with thickness, camber position, and operating alpha
– The optimal camber is non-obvious and depends on the other parameters

p (chordwise position of max camber, in tenths of chord): 1 to 7 (valid: [1, 2, 3, 4, 5, 6, 7])

– Controls where the camber peak (and thus the pressure peak) occurs along the chord
– Lower p: camber peak near the leading edge — steeper pressure recovery
– Higher p: camber peak further aft — gentler pressure gradients
– The combination of m and p determines the camber line shape and thus the pressure distribution, which directly affects stall behavior and L/D

t (max thickness, % of chord): 8 to 18 (valid: [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18])

– Controls the airfoil thickness-to-chord ratio
– Thinner airfoils have less form drag but may stall earlier
– Thicker airfoils have more form drag but can delay flow separation
– The effect on C_L/C_D depends on the operating alpha and camber

alpha_center (operating angle of attack, degrees): 5.0 to 8.0 (valid: [5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0])

– This is the EXACT angle of attack at which C_L/C_D is evaluated
– Each airfoil has a different alpha where L/D peaks, depending on its camber, thickness, and camber position
– Beyond the peak alpha, drag rises due to flow separation and incipient stall, causing C_L/C_D to decline — the rate of decline varies by airfoil shape
– The optimal alpha depends on the airfoil: thicker, higher-camber designs tend to tolerate higher alphas before separation, while thinner, lower-camber designs tend to peak at lower alphas
– You must jointly choose the airfoil shape AND the operating angle

Domain Context (Aerodynamic Reasoning Guide):

You are maximizing C_L/C_D at ONE specific alpha. Key considerations:

1. L/D = C_L/C_D. To maximise this ratio, you want high lift AND low drag. These often conflict: changes that increase lift may also increase drag. The NET effect on the ratio is what matters.
2. Each airfoil shape has a different alpha where L/D peaks. The peak shifts depending on camber, thickness, and camber position.
3. The four parameters (m, p, t, alpha) interact strongly. The optimal value of any one parameter depends on the values of all the others. Do not assume that a parameter value that is good in one combination is good in all combinations.
4. Consider the full range of each parameter. The best design may be at parameter values you would not initially expect based on simple aerodynamic heuristics.

IMPORTANT: These polars come from XFOIL at Re=1,000,000 with ncrit=5 (relatively low turbulence). The C_L/C_D landscape in this regime may differ significantly from textbook intuitions or rules of thumb. In particular: – The relationship between camber level and drag is not straightforward at this Reynolds number – Higher camber does increase pressure drag, but the C_L gains may MORE than compensate, leading to higher C_L/C_D ratios than expected – Your aerodynamic intuition about which parameter ranges are “good” or “bad” may not apply here. – spread your picks across the full range of each parameter to avoid anchoring to potentially incorrect priors.

Appendix B: Experimental Comparisons

Table B1: Comparisons used to examine design choices of the initialization and surrogate frameworks, with an asterisk denoting uninformed condition.

Appendix C: Hypothesis Trajectory

We provide representative hypothesis snapshots from Iterations 1, 5, and 10 of a Cycle-3 run. Iterations (3), (6), and (9) are region-enforcement steps, at which candidate selection is restricted to the LLM-defined region while blended GP-LLM scoring is retained. All text is recorded as produced by the model.

Iteration 1 (best = 115.2)

High C_L/C_D appears strongly correlated with moderate to high camber (𝑚 = 2 − 4) and lower angles of attack (𝛼 = 5.0 − 5.5). The best design uses maximum tested camber (𝑚 = 4) with mid-chord camber position (𝑝 = 4) and high thickness (𝑡 = 16) at minimum angle (𝛼 = 5.0). Higher angles of attack (7.0 − 7.5∘) consistently produce lower C_L/C_D ratios, suggesting these airfoils are most efficient near their design point rather than at high lift coefficients. The combination of significant camber for lift generation and conservative angle of attack for drag minimization appears optimal.

Iteration 5 (best = 128.4)

The data confirms that maximum C_L/C_D occurs with high camber (𝑚 = 4 − 5), moderate camber position (𝑝 = 4 − 5), thick airfoils (𝑡 = 16 − 18), and low angles of attack (𝛼 = 5.0). The top performers all use 𝛼 = 5.0, indicating these highly cambered airfoils operate optimally near their design point where the favorable pressure gradient from camber maximizes lift while thick sections maintain attached flow. The combination of 𝑚 = 5, 𝑝 = 4, 𝑡 = 16 at 𝛼 = 5.0 achieving the highest C_L/C_D of 121.2 suggests this represents near-optimal balance between camber-induced circulation and thickness-enabled pressure recovery. Higher angles of attack dramatically reduce performance due to increased induced drag and potential flow separation on these cambered profiles.

Iteration 10 (best = 147.6)

The data strongly confirms that maximum C_L/C_D requires high camber (𝑚 = 6) with moderate-forward camber position (𝑝 = 4) and thin-to-moderate thickness (𝑡 = 11 − 12) at low angle of attack (𝛼 = 5.0). The top three performers all share 𝑚 = 6, 𝑝 = 4, 𝛼 = 5.0 with only thickness varying, showing a clear optimum at 𝑡 = 11 (C_L/C_D = 142.3) with performance degrading as thickness increases beyond 𝑡 = 12. The consistent superiority of 𝑝 = 4 over 𝑝 = 3 (108.7 vs 142.3 at same 𝑚 = 6, 𝑡 = 16, 𝛼 = 5.0) indicates that moderate-forward camber positioning provides optimal pressure recovery without premature flow separation. Lower camber values (𝑚 ≤ 5) show significantly reduced performance even with optimal thickness, confirming that maximum circulation from high camber is essential for peak efficiency.

Appendix D: Computational Cost and Model Comparison

Table D1 reports LLM calls, median input and output tokens, and estimated dollar cost per run. The surrogate configurations used 448-490 calls and approximately 442,000-532,000 total tokens per 10-iteration run, corresponding to estimated costs of approximately $1.96-$2.34 under the stated rates in Table D2. Initialization required 1-7 calls and an estimated $0.02-$0.18 per generated warm start. Across the reported configurations and repeats, the estimated aggregate LLM cost was $326.70 for the surrogate experiment and $26.50 for initialization ($353.20 total). These figures are based on the listed cost as of July 2026 within Table D2. In addition, model capability is contextualized using the Intelligence Index (Artificial Analysis 2026). We note that this is included only as external model context and is not a measure of performance on the NACA optimization task.

Timing was measured for runs initialized with five observations and followed by ten BO iterations. Across five Cycle-3 runs, median runtime was 323.2 seconds, with median usage of 484 LLM calls, 469,800 input tokens, and 21,067 output tokens. For comparison, ten local GP-only runs required a median 0.836 seconds. The variation in LLM calls reflects occasional candidate-list deduplication.

Table D1: Median LLM usage and cost per reported run

Table D2: Cost and intelligence score breakdown per model