1 Introduction↩︎

Many materials and product design tasks are governed by specification satisfaction rather than single-property optimization or constraint. In these settings, a design is valid so long as its properties lie within prescribed ranges [1], [2]. We refer to these acceptable property intervals as “target ranges.” For example, polymers used in flexible electronics and aerospace composites are often designed around application-specific glass-transition windows [3], [4]; photovoltaic absorbers perform best when their band gaps are near values suited for spectral matching [5]–[7]; and lubricants are viable within suitable viscosity ranges [8], [9]. This range-based perspective is also common in broader design practice, including Ashby-style materials selection, where property ranges guide candidate screening before objective-based ranking [10]; pharmaceutical quality-by-design, where critical quality attributes are specified by acceptable ranges [11]; and product development, where needs are translated into target specifications with ideal and marginally acceptable values [12]. Because several candidates may satisfy the same prescribed range yet differ in synthesis, cost, robustness, scalability, or processability, the utility of design campaigns can be enhanced should they return a portfolio of valid options rather than a single solution [13]–[15]. Thus, the relevant search problem is often not to optimize a single property but to recover, under limited evaluation budgets, distinct candidates whose properties fall within prescribed ranges.

Bayesian optimization (BO) is well suited to materials design because it enables sample-efficient exploration through probabilistic surrogate modeling [16]–[18]. It has been applied successfully to reaction screening [19], [20], alloy formulation [21], [22], and polymer or protein design [23]–[26]. Standard BO acquisition functions, including Expected Improvement (EI) [17], Probability of Improvement (PI) [27], and confidence-bound methods [28], are primarily designed to identify extrema rather than valid regions defined by specifications. Multi-objective BO methods [29]–[31] extend this paradigm to Pareto-front discovery but still prioritize frontier points rather than candidates that lie within acceptable property ranges. Even when Pareto points happen to fall within a target range, interior valid designs that are dominated by frontier points are not preferentially recovered, although they may be equally desirable in practice. Consequently, conventional BO formulations are often misaligned with specification-driven design.

Recent work has extended BO toward goal-seeking and feasibility-aware tasks. Level-set estimation [32] and contour-finding methods [33], [34], including the Expected Feasibility Function (EFF) [35], prioritize threshold crossings or boundary structure rather than discovery of distinct candidates within prescribed ranges. Bayesian Algorithm Execution (BAX) and its multi-target extension MultiBAX [1], [36] encode general target-seeking utilities through posterior sampling but can incur substantial computational cost in continuous or high-dimensional domains and may propose points outside the target range under high predictive uncertainty [2]. Target-constrained approaches such as t-EGO [2] focus sampling around a desired scalar target value but do not naturally address vector-valued specifications or multiple independent target ranges. Constrained BO formulations [37], [38] treat feasibility as a multiplier on a separate objective acquisition, subordinating range-satisfaction probability to an external optimization target. Hierarchical scalarization methods such as Chimera [39] let the user rank the objectives and assign each a tolerance, then fold them into a single objective function. Because the search is driven by this one aggregate score, such methods pursue the objectives in priority order rather than tracking whether each criterion lies within its target range. Together, these limitations motivate a BO framework in which range satisfaction is the primary search signal.

In this work, we develop a range-aware BO framework for specification-driven materials design, built around two range-aware, sampling free acquisition functions: the Tolerance Ball (TB) and the Heaviside (HV). Unlike point-target methods, which preferentially sample near a specified value, or constrained BO approaches, which treat feasibility as a modifier of a separate objective, our framework uses the posterior probability of satisfying a target range as the acquisition objective itself. Because many design campaigns must pursue multiple specifications under a shared experimental or computational budget [23], [26], [40]–[44], including different product grades or operating windows within the same design space, our framework naturally extends to parallel target search across multiple target ranges. We compare the proposed acquisition functions against existing goal-seeking and target-driven methods based on their ability to recover diverse valid designs across target ranges of varying width and difficulty. Finally, beyond idealized benchmarks, we demonstrate the framework in two practical design settings. These include polymerization reaction condition design using kinetic Monte Carlo simulations and sequence-defined oligomer design for prescribed optical absorption bands supported by quantum chemical calculations.

2 Methods↩︎

2.1 Problem Formulation↩︎

2.1.1 Single-Target Inverse Design↩︎

We consider the problem of identifying designs \(\mathbf{x} \in \mathcal{X} \subset \mathbb{R}^M\) whose \(K\)-dimensional property vector lies within a prescribed Euclidean tolerance of the target \(\mathbf{y}^{\mathrm{tgt}} \in \mathbb{R}^K\). Here, \(\mathcal{X}\) denotes the design space, \(\mathbb{R}^M\) is the \(M\)-dimensional real-valued space of possible designs, and \(f \colon \mathcal{X} \to \mathbb{R}^K\) denotes the property map from a design to its predicted properties. A design \(\mathbf{x}\) is considered valid if \[\|f(\mathbf{x}) - \mathbf{y}^{\mathrm{tgt}}\|_2^2 \le \varepsilon^2, \label{eq:valid95solution}\tag{1}\] where \(\|\cdot\|_2\) denotes the Euclidean norm and \(\varepsilon > 0\) is the tolerance radius. When target properties differ in scale or admissible tolerance radius, each output dimension is normalized before applying Eq. 1 . We use the Euclidean norm in the normalized output space because, after normalization, it treats deviations across output dimensions symmetrically and provides a simple measure of overall mismatch from the target. For a fixed target, we denote by \(\mathcal{S} = \{\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N\}\) the set of valid designs identified during the search, where \(N = |\mathcal{S}|\) is the cardinality of \(\mathcal{S}\) (i.e., the number of members of the set). We denote by \(n_{\mathrm{eval}}\) the total number of function evaluations allocated to that target.

2.1.2 Parallel Search for Multiple Targets↩︎

We extend the single-target formulation to a collection of \(T\) targets, \(\mathcal{Y}^{\mathrm{tgt}} = \{\mathbf{y}_1^{\mathrm{tgt}}, \dots, \mathbf{y}_T^{\mathrm{tgt}}\}\). For each target \(t \in \{1,\dots,T\}\), we then define the corresponding valid set by \[\mathcal{S}_t = \left\{ \mathbf{x} \in \mathcal{X} : \|f(\mathbf{x}) - \mathbf{y}_t^{\mathrm{tgt}}\|_2^2 \le \varepsilon^2 \right\},\] and let \(N_t = |\mathcal{S}_t|\) denote the number of valid designs identified for target \(t\). Throughout and for convenience, we use a common tolerance \(\varepsilon\) for all targets, although the formulation readily accommodates target-specific tolerances.

2.2 Performance Metrics↩︎

We quantify search performance by the diversity of valid designs discovered per function evaluation. Diversity is used as an evaluation metric only, not as an explicit term in the acquisition functions [45], [46].

2.2.1 Continuous Design Spaces↩︎

For continuous design spaces, diversity is measured using \(\delta\)-uniqueness. A subset \(\mathcal{S}_{\delta} \subseteq \mathcal{S}\) is \(\delta\)-unique if \[\|\mathbf{x}_i - \mathbf{x}_j\|_2^2 > \delta^2 \label{eq:delta95unique}\tag{2}\] for all distinct \(\mathbf{x}_i, \mathbf{x}_j \in \mathcal{S}_{\delta}\).

Let \(N_{\mathrm{u}}(\delta)\) denote the maximum cardinality of a \(\delta\)-unique subset of \(\mathcal{S}\). We define the raw diversity score as the area under the uniqueness curve, \[A_{\mathrm{u}} = \int_0^{\delta_{\max}} N_{\mathrm{u}}(\delta)\, d\delta, \label{eq:auc95metric}\tag{3}\] where \(\delta_{\max}\) is fixed at \(0.1\phi\), with \(\phi\) denoting the diameter of the design space. For the normalized design space in which each of the \(M\) design variables ranges from 0 to 1, \(\phi = \sqrt{M}\). The normalized diversity score for continuous design spaces is then defined as \[D_\mathrm{c} = \frac{A_{\mathrm{u}}}{n_{\mathrm{eval}}\,\delta_{\max}}.\] When performance is aggregated across multiple target tolerance ranges, we report the geometric mean of the corresponding \(D_\mathrm{c}\) values.

2.2.2 Discrete Design Spaces↩︎

For discrete design spaces, the raw diversity is the number of distinct valid designs found, \(N = |\mathcal{S}|\). Let \(N_{\mathrm{valid}}^{\mathrm{all}}\) denote the total number of valid designs in the full discrete candidate set. Because each function evaluation can identify at most one previously unseen valid design, the maximum attainable number of valid discoveries after \(n_{\mathrm{eval}}\) evaluations is \[N_{\max} = \min \left( n_{\mathrm{eval}},\, N_{\mathrm{valid}}^{\mathrm{all}} \right).\] If \(N_{\mathrm{valid}}^{\mathrm{all}}\) is unknown, we set \(N_{\max} = n_{\mathrm{eval}}\). We then define the normalized diversity score for discrete design spaces as \[D_\mathrm{d} = \frac{N}{N_{\max}}.\]

2.3 Gaussian Process Surrogate Model↩︎

We model the mapping from candidates to properties using Gaussian process (GP) regression [47]. For a \(K\)-dimensional property vector, we fit one GP independently to each output dimension: \[f_k(\mathbf{x}) \sim \mathcal{GP}\bigl(m_k(\mathbf{x}),\, \kappa_k(\mathbf{x},\mathbf{x}')\bigr), \qquad k = 1,\dots,K,\] using a zero prior mean, \(m_k(\mathbf{x}) \equiv 0\). Each GP uses a Matérn 5/2 kernel with signal variance \(\sigma_k^2\) and length scale \(l_k\): \[\kappa_k(\mathbf{x},\mathbf{x}') = \sigma_k^2 \left( 1 + \frac{\sqrt{5}\|\mathbf{x}-\mathbf{x}'\|_2}{l_k} + \frac{5\|\mathbf{x}-\mathbf{x}'\|_2^2}{3l_k^2} \right) \exp\left( -\frac{\sqrt{5}\|\mathbf{x}-\mathbf{x}'\|_2}{l_k} \right).\] The kernel hyperparameters are re-estimated after each BO iteration by maximizing the marginal likelihood. The Matérn 5/2 kernel is a standard and robust default for Bayesian optimization [18]; its twice-differentiable sample paths impose weaker smoothness than the squared-exponential (RBF) kernel [48] and thereby better accommodate the moderately rough structure–property relationships common in molecular data [49].

Given observations \(\{(\mathbf{x}_i, y_{ik})\}_{i=1}^{n}\) for output \(k\), the posterior predictive mean and variance at a new point \(\mathbf{x}_*\) are \[\begin{align} \mu_k(\mathbf{x}_*) &= \mathbf{k}_{k,*}^{\top}\left(\mathbf{K}_k + \sigma_{n,k}^2 \mathbf{I}\right)^{-1}\mathbf{y}_k, \\ s_k^2(\mathbf{x}_*) &= \kappa_k(\mathbf{x}_*, \mathbf{x}_*) - \mathbf{k}_{k,*}^{\top}\left(\mathbf{K}_k + \sigma_{n,k}^2 \mathbf{I}\right)^{-1}\mathbf{k}_{k,*}, \end{align}\] where \(\mathbf{y}_k = [y_{1k},\dots,y_{nk}]^\top\), \(\mathbf{k}_{k,*} \in \mathbb{R}^{n}\) has entries \[(\mathbf{k}_{k,*})_i = \kappa_k(\mathbf{x}_*, \mathbf{x}_i),\] and \(\mathbf{K}_k \in \mathbb{R}^{n \times n}\) has entries \[(\mathbf{K}_k)_{ij} = \kappa_k(\mathbf{x}_i,\mathbf{x}_j).\]

2.4 Bayesian Optimization Procedure↩︎

Bayesian optimization is used to sequentially select candidates for evaluation. For a target \(\mathbf{y}^{\mathrm{tgt}}\), we define the squared discrepancy \[d(\mathbf{x}) = \|f(\mathbf{x}) - \mathbf{y}^{\mathrm{tgt}}\|_2^2\] and the posterior-mean discrepancy \[\Delta^2(\mathbf{x}) = \|\boldsymbol{\mu}(\mathbf{x}) - \mathbf{y}^{\mathrm{tgt}}\|_2^2.\]

For a single target, the BO loop proceeds as follows. First, the surrogate model is initialized with \(n_{\mathrm{init}}\) candidates generated by Latin hypercube sampling (LHS) [50]. The GP surrogate is then fit to the observed data, the acquisition function is maximized over \(\mathcal{X}\), the selected candidates is evaluated, and the surrogate is updated. This process is repeated until the evaluation budget is exhausted. Unless otherwise noted, we use a total budget of 50 evaluations per run, and duplicate evaluations are not allowed.

For parallel target search, all targets are modeled by a single GP surrogate and a shared observation set. Search proceeds as a synchronous batch procedure. At the beginning of each outer iteration, the GP is fit once to all observed data. Each target then proposes one candidate by maximizing a target-specific acquisition function under the shared posterior. The resulting batch of proposed candidates is evaluated, and the GP is updated only after all evaluations in the batch are complete. The posterior therefore remains fixed while the batch is constructed. If two targets propose the same unevaluated candidate, that candidate is evaluated only once and assigned to the earlier target under a fixed ordering. The later target then selects the highest-ranked unevaluated alternative under the same posterior. This duplicate-removal step prevents redundant evaluation while preserving a batch size of \(T\) whenever sufficiently many unevaluated candidates remain.

2.5 Acquisition Functions↩︎

We compare six acquisition strategies: Expected Improvement (EI), Lower Confidence Bound (LCB), the Tolerance Ball (TB) and Heaviside (HV) acquisitions, a posterior-mean BAX baseline (BAX), and random sampling (RS).

2.5.1 Expected Improvement↩︎

Expected Improvement (EI) [17] favors selecting candidates that are expected to improve upon the smallest discrepancy observed so far, with consideration also of the magnitude of that improvement. To obtain a closed-form expression for vector-valued targets, we adopt the isotropic predictive-variance approximation introduced by Uhrenholt et al. [51], treating the predictive distribution of \(\mathbf{f}(\mathbf{x})\) given the current \(n\) observations \(\mathcal{D}_n = \{(\mathbf{x}_i, \mathbf{y}_i)\}_{i=1}^n\) as: \[\mathbf{f}(\mathbf{x}) \mid \mathcal{D}_n \approx \mathcal{N}\left( \boldsymbol{\mu}(\mathbf{x}), \eta^2(\mathbf{x}) \mathbf{I}_K \right),\] where \[\eta^2(\mathbf{x}) = \frac{1}{K}\sum_{k=1}^{K} s_k^2(\mathbf{x}).\] Under this approximation, \[\frac{d(\mathbf{x})}{\eta^2(\mathbf{x})} \sim \chi^2_K\bigl(\lambda(\mathbf{x})\bigr),\] with noncentrality parameter \[\lambda(\mathbf{x}) = \frac{\Delta^2(\mathbf{x})}{\eta^2(\mathbf{x})},\] where \(d(\mathbf{x})\) is the squared discrepancy and \(\Delta^2(\mathbf{x})\) is the posterior-mean squared discrepancy defined above. We then let \[d_{\min} = \min_{i=1,\dots,n} d(\mathbf{x}_i)\] denote the smallest observed discrepancy, and define \[r = \frac{d_{\min}}{\eta^2(\mathbf{x})}.\] The EI acquisition is then \[A_{\mathrm{EI}}(\mathbf{x}) = d_{\min}\, F_{K,\lambda}(r) - \eta^2(\mathbf{x}) \left[ K\,F_{K+2,\lambda}(r) + \lambda(\mathbf{x})\,F_{K+4,\lambda}(r) \right],\] where \(F_{K,\lambda}(\cdot)\) denotes the cumulative distribution function of the noncentral \(\chi^2\) distribution with \(K\) degrees of freedom and noncentrality parameter \(\lambda(\mathbf{x})\).

2.5.2 Lower Confidence Bound↩︎

As a complementary baseline, we use the closed-form lower confidence bound (LCB) acquisition [51], \[A_{\mathrm{LCB}}(\mathbf{x};\beta) = -\bigl(\alpha - \beta \rho\bigr)^{1/\ell} \bigl(K + \lambda(\mathbf{x})\bigr)\eta^2(\mathbf{x}),\] where \(\beta = 1\) controls the exploration–exploitation trade-off. The auxiliary quantities are \[\begin{align} \ell &= 1 - \frac{r_1 r_3}{3 r_2^2}, \\ r_s &= 2^{s-1}(s-1)!\,\bigl(K + s\,\lambda(\mathbf{x})\bigr), \qquad s = 1,2,3, \\ \alpha &= 1 + \ell(\ell-1)\left[ \frac{r_2}{2 r_1^2} - (2-\ell)(1-3\ell)\frac{r_2^2}{8 r_1^4} \right], \\ \rho &= \frac{\ell r_2^2}{r_1} \left[ 1 - \frac{\ell(1-3\ell)}{4 r_1^2} r_2 \right]. \end{align}\]

2.5.3 Tolerance-Aware Acquisitions↩︎

We consider two tolerance-aware, sampling-free acquisitions computed directly from the GP posterior. The first, Tolerance Ball (TB), is the posterior probability that a candidate satisfies the output-space target range. Using the same isotropic predictive-variance approximation as in the EI derivation, we obtain \[A_{\mathrm{TB}}(\mathbf{x}) = \mathbb{P}\left[d(\mathbf{x}) \le \varepsilon^2\right] = F_{K,\lambda(\mathbf{x})} \left( \frac{\varepsilon^2}{\eta^2(\mathbf{x})} \right).\] For \(K=1\), this reduces to the probability that the scalar posterior lies within a tolerance interval around the target; for \(K>1\), it is the corresponding probability that the posterior vector lies inside a Euclidean ball. This is analogous to probability-of-feasibility criteria widely used in constrained Bayesian optimization [37], [38], [52], with feasibility here defined as membership in the prescribed output-space tolerance region.

The second, Heaviside (HV), prioritizes candidates whose posterior mean lies inside the tolerance range while retaining a smooth transition near the boundary. We define \[A_{\mathrm{HV}}(\mathbf{x}) = \bigl(1 - w(\mathbf{x})\bigr) + w(\mathbf{x})\,A_{\mathrm{TB}}(\mathbf{x}),\] where \[w(\mathbf{x}) = \frac{1}{2} \left[ 1 + \tanh\left( \frac{\Delta^2(\mathbf{x}) - \varepsilon^2}{\tau} \right) \right]\] and \(\tau = 10^{-3}\varepsilon^2\) is a smoothing parameter that controls the width of the transition region around the tolerance boundary. Thus, \(A_{\mathrm{HV}}(\mathbf{x}) \approx 1\) for candidates whose posterior mean lies well inside the tolerance range \(\bigl(\Delta^2(\mathbf{x}) \ll \varepsilon^2\bigr)\), whereas \(A_{\mathrm{HV}}(\mathbf{x}) \approx A_{\mathrm{TB}}(\mathbf{x})\) for candidates whose posterior mean lies well outside it \(\bigl(\Delta^2(\mathbf{x}) \gg \varepsilon^2\bigr)\). At the boundary, \(\Delta^2(\mathbf{x}) = \varepsilon^2\), so \(w(\mathbf{x}) = 1/2\) and \(A_{\mathrm{HV}}(\mathbf{x})= \frac{1}{2}(1+A_{\mathrm{TB}}(\mathbf{x}))\).

Notably, TB and HV are pointwise feasibility acquisitions; they do not explicitly penalize similarity to previously discovered valid candidates. Any diversity observed in the discovered set therefore arises implicitly from posterior uncertainty, the geometry of the feasible regions, and duplicate-evaluation exclusion. The approaches would thus be complementary to other strategies to promote solution diversity.

2.5.4 Posterior-Mean BAX Baseline↩︎

As a baseline for tolerance-constrained design discovery, we include a feasible-set exploration strategy based on a posterior-mean variant of BAX [1], [36]. Given that the corresponding feasible set of valid candidates is \[\mathcal{T} = \left\{ \mathbf{x} \in \mathcal{X} : d(\mathbf{x}) \le \varepsilon^2 \right\},\] we approximate [1], [36] this set by thresholding the GP posterior mean, resulting in \[\bar{\mathcal{T}} = \left\{ \mathbf{x} \in \mathcal{X} : \Delta^2(\mathbf{x}) \le \varepsilon^2 \right\}.\]

For continuous design spaces, \(\bar{\mathcal{T}}\) is approximated on a set of 1,000 candidate points sampled uniformly from \(\mathcal{X}\) at each BO iteration. For discrete design spaces, the feasible-set approximation is evaluated over the full candidate set. The resulting BAX-style acquisition is defined as \[A_{\mathrm{BAX}}(\mathbf{x}) = \begin{cases} \dfrac{1}{K}\sum_{k=1}^{K} s_k(\mathbf{x}), & \text{if } \mathbf{x} \in \bar{\mathcal{T}} \setminus \mathcal{X}_{\mathrm{obs}}, \\ 0, & \text{otherwise} \end{cases}\] where \(\mathcal{X}_{\mathrm{obs}}\) is the set of previously evaluated candidates. In the first logic branch, ‘\(\setminus\)’ indicates the set-difference operator, such that the expression is evaluated over the set of candidates in \(\bar{\mathcal{T}}\) that have not already been evaluated. If \(\bar{\mathcal{T}} \setminus \mathcal{X}_{\mathrm{obs}} = \varnothing\) (i.e., there are no unevaluated candidates in the target set), the acquisition defaults to global uncertainty sampling: \[A_{\mathrm{US}}(\mathbf{x}) = \frac{1}{K}\sum_{k=1}^{K} s_k(\mathbf{x}).\]

2.5.5 Random Sampling↩︎

Random sampling selects the next unevaluated candidate at random from the available candidate set with uniform probability. For continuous design spaces, the candidate set consists of 1,000 points drawn with uniform probabilities from \(\mathcal{X}\) at each iteration. For discrete design spaces, the candidate set consists of all remaining unevaluated candidates.

2.6 Datasets↩︎

We evaluated the BO methods on four classes of inverse-design tasks. In each, the goal is to find candidates whose properties fall within a prescribed tolerance of a target specification: (i) synthetic benchmark functions (Branin, Hartmann-3, Ackley-5, Layeb-6) on continuous domains, for which the target is a scalar function value; (ii) pool-based datasets of pre-characterized candidates from nanoparticle synthesis, intrinsically disordered proteins, polymer topology, and small-molecule libraries, for which targets are measured or computed physicochemical properties; (iii) kinetic Monte Carlo (KMC) styrene polymerization, for which the design variables are reaction conditions and the target is the shape of the polymer molecular weight distribution; and (iv) sequence-defined conjugated oligomers, for which the design variables are molecular building-block choices and the target is the shape of the UV–vis absorption band. In the following, we describe these design tasks and the underlying functions in more detail.

2.6.1 Synthetic Functions↩︎

We considered four analytic benchmark functions, Branin (\(\mathbb{R}^2 \to \mathbb{R}\)), Hartmann-3 (\(\mathbb{R}^3 \to \mathbb{R}\)), Ackley-5 (\(\mathbb{R}^5 \to \mathbb{R}\)), and Layeb-6 (\(\mathbb{R}^6 \to \mathbb{R}\)). The analytic expressions and domains of all four functions are given in the SI Section S1. Fig. 1 shows representative two-dimensional slices with the target levels overlaid. These benchmark functions span a range of input dimensionalities and landscape characteristics, including differing degrees of multimodality and roughness, and they induce inverse-design tasks with multiple valid solutions that may lie in distinct regions of the design space [53].

2.6.2 Pool-Based Datasets↩︎

For pool-based experiments, the design space is a fixed finite candidate set. At each BO iteration, one previously unevaluated candidate is selected from the pool and evaluated. The datasets used were nanoparticle synthesis (NS), intrinsically disordered proteins (IDP), TopoRg, BACE, ESOL, FreeSolv, Lipophilicity (Lipo), and QM9.

SI Table S1 summarizes the number of candidates, the design-space representation, and the target property or properties for each dataset. The NS dataset comprises 1,997 candidate synthesis conditions, each defined by four normalized design variables, Ti(TEOA)\(_2\) concentration, TEOAH\(_3\) concentration, pH, and temperature [54]. The target outputs are nanoparticle radius and polydispersity index. The IDP dataset includes 2,031 polypeptide sequences of length 20–50 residues. Each sequence is characterized by three normalized phase-behavior properties, radius of gyration (\(\tilde{R}_g\)), second virial coefficient (\(\tilde{B}_2\)), and expenditure density (\(\rho^*\)) [55]. The TopoRg dataset consists of 1,340 candidates represented by eight-dimensional latent vectors obtained from a variational autoencoder trained on graph representations of polymer architectures. The target property is the mean single-chain radius of gyration [56]. The remaining datasets (BACE, ESOL, FreeSolv, Lipophilicity, and QM9) are molecular libraries with experimentally measured or computationally derived physicochemical properties [57]–[61]. The specific target property used for each dataset is listed in Table S1.

2.6.3 Kinetic Monte Carlo Styrene Polymerization↩︎

We considered a task in which the objective is to identify reaction conditions that reproduce a prescribed molecular weight distribution (MWD). To this end, we modeled styrene polymerization using a simplified two-stage radical polymerization process inspired by the sequential controlled/conventional strategy reported by Lenzi et al. [62]. The first stage is a nitroxide-mediated polymerization (NMP) step, and the second stage is a conventional free-radical polymerization (FRP) step. This model is intended as a stylized batch, bulk, isothermal benchmark that captures the mechanistic contrast between controlled radical growth in stage one and conventional radical propagation in stage two, rather than as an exact reproduction of the heterogeneous reactor configurations used experimentally by Lenzi et al. [62], [63].

The complete reaction networks and kinetic parameters for both stages are listed in SI Tables S2–S4, using literature-based approximations for styrene NMP and FRP kinetics. The KMC scheme follows the stochastic simulation algorithm originally developed by Gillespie [64], [65], with the polymerization formulation reported by Romero et al. [66]. Each candidate is represented by a five-dimensional design vector \[\mathbf{x} = \left( X_{\mathrm{NMP}}, T_{\mathrm{NMP}}, C_{\mathrm{init,NMP}}, T_{\mathrm{FRP}}, C_{\mathrm{init,FRP}} \right),\] where \(X_{\mathrm{NMP}}\) denotes the prescribed monomer conversion at the end of the NMP stage, \(T_{\mathrm{NMP}}\) and \(T_{\mathrm{FRP}}\) are the reaction temperatures in the two stages, and \(C_{\mathrm{init,NMP}}\) and \(C_{\mathrm{init,FRP}}\) are the initiator concentrations in the NMP and FRP stages, respectively. Temperatures are reported in units of Kelvin (K), and initiator concentrations are reported in molar units (mol L\(^{-1}\)). The FRP stage was terminated at an overall conversion of 95%.

To generate the polymerization design space, we enumerated a full-factorial grid over these five variables. The prescribed NMP conversion \(X_{\mathrm{NMP}}\) was varied from 0.4 to 0.8 using 8 uniformly spaced values. The NMP temperature \(T_{\mathrm{NMP}}\) was varied from 353 to 403 K using 10 uniformly spaced values, and the FRP temperature \(T_{\mathrm{FRP}}\) was varied from 313 to 363 K using 10 uniformly spaced values. The initial monomer concentration was 8.56 mol L\(^{-1}\), and the KMC simulation was initialized with a control volume of \(8.86\times10^{-16}\) L. The initiator concentrations \(C_{\mathrm{init,NMP}}\) and \(C_{\mathrm{init,FRP}}\) were determined from logarithmically spaced monomer-to-initiator ratios of 50 to 200 for NMP and 500 to 2000 for FRP, respectively, yielding 8 values for each stage. This corresponds to \(C_{\mathrm{init,NMP}}\) values from \(4.28\times10^{-2}\) to \(1.712\times10^{-1}\) mol L\(^{-1}\) and \(C_{\mathrm{init,FRP}}\) values from \(4.28\times10^{-3}\) to \(1.712\times10^{-2}\) mol L\(^{-1}\). The resulting full-factorial grid comprised 51,200 unique reaction conditions.

For a given candidate \(\mathbf{x}\), the KMC simulation returns the final MWD in weight-fraction form. We represent this output as a 100-dimensional vector \(\mathbf{y}\), where each component \(y_b\) is the polymer weight fraction within the \(b\)-th bin. The molecular weight of a chain with degree of polymerization \(n\) is \(M = M_\mathrm{M}\,n\), where \(M_\mathrm{M} = 104.15\) g mol\(^{-1}\) is the styrene monomer molar mass (SI Table S2). The bins are defined on a logarithmic molecular-weight axis spanning \(M_\mathrm{M}\) to \(10^{5}M_\mathrm{M}\) (approximately \(10^{2}\) to \(10^{7}\) g mol\(^{-1}\)), with 100 uniformly spaced intervals in log-10 space. Specifically, bin \(b\) corresponds to the interval \[\log_{10}(M) \in [\xi_b,\xi_{b+1}),\] where \(\{\xi_b\}_{b=1}^{101}\) are uniformly spaced between \(\log_{10}(M_\mathrm{M})\) and \(5 + \log_{10}(M_\mathrm{M})\). The output vector is normalized such that \[\sum_{b=1}^{100} y_b = 1.\] This binned weight-fraction output serves as the property vector used for inverse design. We treat the 100-bin vector as a fixed numerical descriptor of the MWD rather than as a full probabilistic representation.

2.6.4 Sequence-Defined Conjugated Oligomers↩︎

We considered a task defined over a library of sequence-defined conjugated oligomers. The library was constructed from 21 synthetically accessible oligomer cores and 100 commercially available pinacol boronic esters that serve as terminal caps. Combining each core with each cap yields 2,100 candidate structures; after removing duplicate SMILES, which arise when distinct core–cap pairings are structurally equivalent, 1,980 unique oligomers remain and constitute the optimization library (Fig. 6a and SI Fig. S15).

The objective is the shape of the simulated UV–vis absorption band. For each oligomer we computed a broadened absorption spectrum over 240–500 nm by time-dependent density functional theory (TD-DFT) [67], [68] and represented it by a compact, physically interpretable parametrization. Working in photon-energy space (eV), where the electronic bands are well described as Gaussian, we fit each spectrum with \(K=5\) Gaussian components, \[S(E) = \sum_{k=1}^{K} A_k \exp\!\left[-\frac{(E-\mu_k)^2}{2\sigma^2}\right], \label{eq:gaussian95spectrum}\tag{4}\] where \(A_k\) and \(\mu_k\) are the amplitude and center energy of band \(k\), and the width \(\sigma = 0.12\) eV is held fixed for all bands and all molecules, matching the broadening used to construct the spectra. Spectra contain between one and five resolved bands, which we encode in a fixed-length representation by ordering the components by descending amplitude. When a spectrum has fewer than five bands, each unused component is assigned zero amplitude (\(A_k = 0\)) together with a single canonical center energy common to all molecules (\(3.35\) eV). Because this canonical center is constant across the library, it has zero variance and is therefore uninformative to the surrogate, and it cancels exactly in the tolerance distance whenever both a target and a candidate lack the corresponding band. Moreover, \(S(E)\) is independent of \(\mu_k\) when \(A_k = 0\), such that the center of an absent band has no effect on the absorption profile being matched; \(\mu_k\) influences the objective only when band \(k\) is physically present (\(A_k > 0\)).

Each spectrum is peak-normalized such that its most intense band has unit amplitude (\(A_1 \equiv 1\)). Because \(A_1\) is constant across the library, it is discarded. The resulting nine-dimensional optimization target \((\mu_1, A_2, \mu_2, A_3, \mu_3, A_4, \mu_4, A_5, \mu_5)\) comprises the position of the dominant band together with the relative amplitudes and positions of the four secondary bands. This representation captures both where a molecule absorbs and the overall shape of its absorption envelope, allowing specifications to be posed over an entire spectral band rather than over an isolated \(\lambda_{\max}\). Before optimization the target vector is scaled group-wise to \([0,1]\), with all band amplitudes sharing one scaling and all band energies sharing another, such that the dimensionless amplitudes and the eV-scale energies contribute comparably to the tolerance distance rather than being normalized per dimension. The five resulting targets \(T_1\) through \(T_5\) are listed in SI Table S6 and displayed in Fig. S16.

2.6.4.1 Time-Dependent Density Functional Theory Calculations.

For each oligomer, long alkyl side chains were truncated to methyl groups to reduce computational cost [69], explicit hydrogen atoms were added, and 100 three-dimensional conformers were generated in RDKit [70] using the ETKDGv3 embedding procedure [71], each subsequently optimized with the MMFF94 force field [72]; the lowest-energy conformer was retained. To promote planarity of the conjugated backbone while preserving side-chain flexibility, dihedral-angle constraints (\(0^\circ\) or \(180^\circ\), chosen from the force-field-optimized conformer) were applied only to the single bonds linking distinct ring or fused-ring systems [73]. Constrained geometry optimizations were then performed in ORCA [74] using GFN2-xTB [75] with ALPB implicit solvation [76] in dichloromethane. Vertical excitation energies and oscillator strengths for the 15 lowest singlet excited states were computed for the optimized geometries by TD-DFT at the \(\omega\)B97X-3c/CPCM(dichloromethane) level [77], [78] in ORCA. To obtain a continuous absorption profile rather than a discrete stick spectrum, only transitions above 240 nm were broadened (each by a Gaussian of width \(\sigma = 0.12\) eV) and summed on a common grid spanning 240–500 nm, such that intense deep-UV bands below the cutoff did not contaminate the profile. The resulting broadened spectrum is the quantity subsequently parametrized by the five-Gaussian fit of Eq. 4 .

2.6.5 Input Representation and Preprocessing↩︎

All input variables were mapped to the unit hypercube \([0,1]^M\) using min–max normalization performed separately for each dataset. For the molecular datasets and the oligomer case study, each candidate was represented by the full set of molecular descriptors computed with Mordred from its chemical structure [79], without manual feature selection. Descriptor columns with zero variance were removed, and principal component analysis (PCA) [80] was then performed separately for each dataset, retaining the number of components needed to explain 95% of the total variance; the surrogate therefore operates on these principal components rather than on individual descriptors. For pool-based datasets, output properties were standardized to zero mean and unit variance separately for each property before target selection and tolerance definition. This standardization makes tolerance values comparable across datasets and across property dimensions. For analytic benchmark functions and case studies without a predefined output scale, no output normalization was applied.

2.6.6 Target and Tolerance Selection↩︎

For each dataset, we selected five representative targets that span the output space. For pool-based datasets, we applied \(k\)-medoids clustering [81] with \(k=5\) to the standardized property vectors of all candidates and used the resulting medoids as targets. For synthetic functions, we first sampled 20,000 inputs uniformly from the domain, evaluated the corresponding function values, and then applied \(k\)-means clustering [82] with \(k=5\) in output space. The resulting cluster centers were used as targets. We defined a dataset-specific base tolerance \(\varepsilon_0\) as the geometric mean of the pairwise Euclidean distances between the five targets. Performance was then evaluated at three tolerance levels, obtained by scaling this base tolerance by a tolerance ratio \(r\), \[\varepsilon = r\,\varepsilon_0, \qquad r \in \{0.2,\,0.4,\,0.6\},\] where a larger tolerance ratio \(r\) yields a wider acceptance range and less stringent specification.

3 Results↩︎

3.1 Tolerance–Ball Acquisition Achieves the Best Overall Discovery Performance↩︎

We first compare the acquisition functions on the benchmark suite, comprising the four analytical functions and the eight pool-based datasets, in terms of how effectively they recover multiple distinct valid designs under a fixed evaluation budget. This is the central criterion in our setting because specification-driven design often admits more than one acceptable solution, and a useful method should uncover a broad set of valid candidates rather than repeatedly refine a single promising region. Accordingly, performance here refers to the normalized diversity score (\(D_\mathrm{c}\) for continuous tasks, \(D_\mathrm{d}\) for discrete tasks) of valid discoveries per evaluation, aggregated across tasks and tolerance ratios \(r\).

Figure 2 shows that TB exhibits the strongest overall discovery performance across the benchmark suite. It achieves the lowest average rank at all three tolerance ratios (Fig. 2a) and the highest fraction of best-performing tasks (Fig. 2b), ranging from 0.42 at \(r=0.4\) to 0.67 at \(r=0.6\). LCB and HV form a competitive second tier but neither matches TB on either aggregate metric. In contrast, RS performs worst overall and accounts for the majority of worst-performing outcomes (Fig. 2c), with EI and BAX also producing occasional last-place outcomes among the non-random methods.

The performance gap between TB and HV highlights the importance of retaining probabilistic resolution within and near the tolerance range. Both acquisitions are designed for targeting ranges of valid discoveries, but HV trails TB at every tolerance ratio. One likely reason is that HV assigns similar acquisition values to many candidates whose posterior means lie inside the range, even when those candidates have different posterior probabilities of satisfying the range. This behavior can reduce its ability to prioritize the candidates with the highest probability of range satisfaction. This loss of resolution parallels the classical trade-off between Probability of Improvement (PI) and Expected Improvement (EI), in which PI saturates once a candidate is likely to improve and therefore disregards the magnitude of improvement that EI captures [83]. TB avoids this loss of resolution by directly scoring the posterior probability of range satisfaction, which provides a finer ranking of candidates according to their likelihood of meeting the specification.

The stronger performance of LCB relative to EI further indicates that extremum-seeking improvement is poorly aligned with range discovery. Although both methods use target discrepancy, EI evaluates improvement relative to the best discrepancy observed so far. This criterion is well suited to refining a single incumbent solution, but it can concentrate search near an already promising region after an acceptable candidate has been found. LCB combines posterior mean and uncertainty, allowing several plausible target-consistent regions to remain competitive during the search. This broader exploration behavior is more compatible with the goal of discovering multiple acceptable candidates, which likely explains why LCB consistently outperforms EI across tolerance ratios.

BAX underperforms relative to TB and LCB despite its set-oriented objective. This likely reflects its reliance on a posterior-mean approximation of the feasible region. When that approximation is incomplete, the acquisition may favor uncertainty reduction over candidates with high probability of satisfying the target range. Consequently, evaluations can be allocated to refining the feasible-set estimate rather than to recovering valid designs. This behavior is consistent with the weaker aggregate performance of BAX in Fig. 2.

3.2 Tolerance–Ball Acquisition Maintains Balanced Discovery Across Targets↩︎

Figure 3 examines the per-target structure of discovery performance using the Branin benchmark and the QM9 molecular dataset as representative examples from the continuous and discrete settings, respectively; radar plots for the remaining benchmark tasks are provided in SI Section S4. The aggregate advantage of TB arises from its ability to sustain strong discovery across multiple targets rather than concentrating performance on only a subset. At an intermediate tolerance ratio (\(r = 0.4\)), TB achieves a high mean normalized diversity score \(D_\mathrm{c}\) of 0.42 across the five targets in the Branin task (Fig. 3a), with comparable coverage across all five targets. HV and LCB also perform well in this low-dimensional setting, with mean scores of 0.41 and 0.42, whereas EI achieves a mean of 0.20 and shows pronounced concentration on \(T_1\). BAX achieves an intermediate mean score of 0.36. A similar pattern appears in the QM9 inverse-design task (Fig. 3b), where TB again combines the highest mean score (0.15) with comparatively even performance across the full target set, while EI, BAX, and RS achieve lower scores and concentrate on a subset of targets. This balance is reflected in the across-target variance of \(D_\mathrm{c}\) and \(D_\mathrm{d}\) (SI Fig. S1). RS and EI show the largest across-target variation at moderate to high tolerance ratios, consistent with discovery concentrated on a subset of targets, whereas TB falls below both at \(r = 0.4\) and \(r = 0.6\). At \(r = 0.2\), the methods are comparable. These results indicate that TB not only improves overall discovery but also tends to distribute discovery outcomes more evenly among targets.

The discovered-point plots further indicate that TB recovers candidates that are both valid and well distributed across the tolerance regions. In the Branin example (Fig. 3c), TB identifies candidates across multiple disjoint acceptable regions while keeping the discovered points closely confined to the prescribed tolerance ranges. HV also finds candidates in multiple regions, although its discoveries are less consistently confined to the admissible ranges. LCB selects many points near the target contours, but several fall outside the tolerance regions, consistent with a search bias toward the target contour rather than the full range-satisfaction criterion. BAX yields an even more diffuse distribution of discovered points, with a larger fraction of points outside the admissible regions.

The QM9 property-space visualization reveals the same qualitative pattern in a discrete molecular dataset. In Fig. 3d, TB produces a denser and more localized occupancy around the target locations, whereas LCB and HV generate more diffuse point clouds and BAX is the most dispersed. These visualizations support the interpretation that TB improves discovery by prioritizing candidates with high posterior probability of satisfying the tolerance criterion, rather than by broadly sampling near the target region. Additional radar plots in SI Section S4 show consistent behavior across the remaining benchmark tasks.

3.3 Off-Target Hit Dynamics Reveal Target Fidelity During Parallel Target Search↩︎

In principle, algorithms for parallel target search should yield not only efficient discovery of valid candidates but also faithful allocation of those discoveries to their intended targets. In other words, the efficacy of an algorithm should not necessarily hinge on “accidental” discoveries where suitable candidates for one objective originate from the search for another. Figure 4 tracks this property through the cross-target hit ratio. For each candidate proposed for a given target, we recorded whether the resulting property vector fell inside the tolerance range of another target. Because the target ranges considered here are disjoint, such an event indicates that the acquisition selected a valid design for the wrong target. The evolution of this ratio over BO iterations therefore provides a measure of target fidelity. Methods with low off-target hit ratios more consistently preserve the intended target, whereas methods with high ratios more frequently recover valid designs from other target ranges.

Target fidelity reflects both the acquisition function and the structure of the function landscape. In the Branin task (Fig. 4a), TB, HV, LCB, and BAX maintain near-zero off-target hit ratios across all tolerance ratios, whereas RS rises rapidly to a high plateau, which is expected since its search strategy is not actually targeted, and EI increases more gradually. The QM9 task (Fig. 4b) shows the same qualitative ordering between RS and the model-guided methods but the separation is smaller, and BAX and EI also drift upward at larger tolerance ratios. This behavior suggests that the design regions associated with different targets are less distinctly separated in the chosen representation. The nanoparticle synthesis (NS) task (Fig. 4c) similarly preserves the advantage of TB, HV, and LCB, although BAX and EI here climb closer to RS. In the Lipophilicity (Lipo) task (Fig. 4d), all methods show more rapid growth in off-target hit ratio, indicating that the property landscape or candidate distribution makes target-specific discrimination intrinsically more difficult. Even in this more challenging setting, TB and LCB attain the lowest off-target hit ratios.

Higher off-target hit ratios are associated with weaker target-specific discovery in these experiments. RS produces the highest off-target hit ratios and also shows the weakest overall discovery performance, whereas TB combines the lowest off-target hit ratios with the strongest aggregate discovery (Fig. 2). This relationship suggests that the advantage of TB is unlikely to arise from opportunistically recovering valid designs from any target range but instead from selecting candidates with high posterior probability of satisfying the intended target range. The off-target analysis therefore provides additional evidence that explicit range-satisfaction scoring can improve both discovery efficiency and target fidelity in parallel target search.

3.4 Case Study: Polymer Molecular Weight Distribution Design↩︎

Whereas the preceding benchmark tasks are drawn from existing datasets, we now turn to case studies designed to more closely reflect real-world considerations. The KMC polymerization task addresses inverse polymer design for a prescribed full molecular weight distribution (MWD). Rather than optimizing a scalar property such as number-average molecular weight or dispersity, the objective here is to identify reaction conditions that reproduce a target MWD. Designing the full MWD is well motivated, because its breadth, skew, and modality influence polymer properties beyond what summary statistics capture [84], [85]. This goal has been pursued through controlled-synthesis and flow strategies [86]–[89] and more recently through machine-learning approaches, including reinforcement learning and KMC-based active learning and multi-objective reverse engineering [90]–[94]. However, recovering reaction conditions for a target MWD under an explicit, range-aware tolerance criterion remains comparatively less explored. Here, each specification is a 100-bin weight-fraction distribution, and a design counts as valid when its simulated distribution falls within the prescribed Euclidean tolerance of the target. Because the design space is substantially larger (51,200 candidates) and valid regions are sparser than in the pool-based benchmarks, we extended the evaluation budget to 200 evaluations for this case study. The five target distributions span a range of shapes, including unimodal profiles centered at different molecular weights and bimodal profiles with a dominant peak and a weaker high-molecular-weight shoulder (Fig. 5a). Figure 5b also shows that the residuals between target and recovered distributions remain small in magnitude across the full molecular-weight range, indicating that the Euclidean tolerance criterion captures the principal features of the distribution in this fixed representation.

TB provides the strongest and most balanced recovery of valid polymerization conditions. In the radar plot of normalized diversity score (Fig. 5c), TB achieves the highest mean \(D_\mathrm{d}\) of 0.93 across the five targets, followed by HV at 0.85, LCB at 0.69, and BAX at 0.64. EI and RS perform poorly, with mean scores of 0.10 and 0.04 respectively, reflecting that valid reaction conditions occupy only a small region of the five-dimensional design space and are unlikely to be recovered by extremum-seeking or unguided search. The target-resolved discovery trajectories further show that TB identifies valid designs more rapidly and consistently than the other methods. Figure 5c reports the cumulative number of valid designs found for each target over 200 BO iterations. TB rises steeply during the early stage of the search and reaches the largest, or near-largest, number of valid designs for most targets by the end of the evaluation budget. LCB and HV also accumulate valid designs at competitive rates, whereas BAX trails the leading methods. EI and RS remain near zero throughout the search.

The valid designs discovered by TB are also diverse in reaction-condition space, even when they produce similar MWDs. For each target, a representative pair of valid designs that yield nearly identical distributions is reported in SI Fig. S14; the two designs differ across all five reaction variables (SI Table S5), confirming that distinct combinations of conditions can satisfy the same product-level specification. The UMAP [95] projection in Fig. 5d illustrates this, with the valid designs for every target occupying multiple separated regions of the embedding. This degeneracy is valuable for process design because it reveals multiple reaction pathways to the same product specification and provides flexibility for subsequent optimization of cost, robustness, scalability, or experimental feasibility.

3.5 Case Study: Sequence-Defined Conjugated Oligomers Discovery↩︎

The sequence-defined oligomer task tests range-aware discovery in a finite molecular library, where discrete changes in core sequence or terminal-cap chemistry reshape the entire optical absorption band rather than merely shifting its peak. Accordingly, each target is specified not by a single wavelength but as a region of the spectral-shape space introduced above, that is, the nine-dimensional target vector \((\mu_1, A_2, \mu_2, A_3, \mu_3, A_4, \mu_4, A_5, \mu_5)\) comprising the dominant-band energy together with the relative amplitudes and positions of the four secondary bands; a representative spectrum and its five-Gaussian decomposition are shown in Fig. 6b. The residual between the simulated spectrum and this five-Gaussian fit (Fig. 6b, bottom) remains small and largely structureless, confirming that the compact parametrization reproduces the absorption band. The five targets \(T_1\)–\(T_5\) were obtained by \(k\)-medoids clustering of the Gaussian spectral vectors and span the dominant-band variation observed in the library (see SI Table S6 and Fig. S16).

Performance differences among acquisition functions are less pronounced in this discrete molecular setting than in the polymerization case study (Fig. 6c). TB attains the highest mean normalized diversity score (\(D_\mathrm{d} = 0.32\)) and remains among the strongest methods across the five target bands, but the separation among methods is modest and random sampling is competitive for several targets. The target-resolved discovery trajectories in SI Fig. S17 corroborate this view. From this data, our observation is that no single acquisition function dominates all bands across the full evaluation budget, although TB and HV sustain higher discovery rates.

The TB discoveries nevertheless reveal a chemically important form of degeneracy, where structurally diverse oligomers can lead to comparable absorption profiles. Figures 6d and 6e present three TB-discovered oligomers for targets \(T_1\) and \(T_3\), respectively. Within each panel the molecules differ in conjugated-core sequence and terminal-cap chemistry, yet their simulated spectra are similar, sharing an almost identical dominant-band position to within the prescribed tolerance. This many-to-one mapping from molecular structure to optical response is valuable from a design perspective because it provides multiple synthetic routes to a single absorption specification, leaving room for subsequent selection according to criteria outside the present objective, such as synthetic accessibility or stability.

4 Conclusion↩︎

Many materials and product design problems are guided by specification satisfaction. A candidate is valuable when its properties fall within prescribed tolerance ranges, rather than when a single objective reaches an optimum. We reformulated this challenge as a parallel target-range discovery problem and derived a closed-form Tolerance Ball (TB) acquisition that directly scores the posterior probability of range satisfaction through a noncentral \(\chi^2\) approximation. Because the acquisition is sampling-free and analytically tractable, it offers a computationally efficient alternative to extremum-seeking strategies and sampling-based set-estimation methods.

Across synthetic functions, molecular libraries, and case studies, TB consistently achieved the strongest and most balanced recovery of valid designs, outperforming standard Bayesian optimization baselines and recent goal-oriented approaches at multiple tolerance ratios. This advantage likely arises because TB explicitly scores range satisfaction for the intended target, yielding low cross-target hit ratios and high target fidelity. Improvement-based or uncertainty-driven methods, in contrast, can drift toward off-target ranges or concentrate on refining a single incumbent solution. In the large, sparse design space of the kinetic Monte Carlo polymerization task, TB recovered distinct reaction conditions that accurately reproduced prescribed molecular weight distributions (MWDs). This distributional target class has practical relevance for optimizing mechanical, rheological, and processing behavior. A similar structure–property degeneracy emerged in the sequence-defined oligomer library, where chemically distinct molecules produced nearly identical absorption bands, providing multiple molecular candidates for a single optical specification.

Although TB is not always the top-ranked method for every individual target or dataset, it achieves strong performance across the full range of settings considered, making it a reliable default choice for specification-driven search when the structure of the feasible region is not known in advance. The benefit of surrogate-guided search also depends on the structure of the design space: when many valid molecules fall within the target ranges in a discrete library, the performance gap between TB and random sampling narrows, indicating that model guidance adds the most value when valid regions are sparse or disconnected. Future work could incorporate explicit diversity incentives into the acquisition, relax the isotropic variance approximation through full multi-output Gaussian process models, and introduce adaptive tolerances that reflect evolving design specifications. Integrating range-aware optimization with automated synthesis and characterization pipelines could enable closed-loop, specification-driven materials discovery.

5 Acknowledgments↩︎

S.J. and M.A.W. acknowledge support from the National Science Foundation under Grant No. 2118861. Simulations and analyses were performed using resources from Princeton Research Computing at Princeton University, which is a consortium led by the Princeton Institute for Computational Science and Engineering (PICSciE) and Office of Information Technology’s Research Computing. All code required to recreate the results above can be obtained at https://github.com/webbtheosim/target-range-optimization.

References↩︎