JAREX: A N ACQUISITION F UNCTION FOR M ULTI -O BJECTIVE A LGORITHMIC P ROCESS C HARACTERIZATION A P REPRINT
arXiv:2609.24954v1 [stat.ML] 21 Sep 2026
Xinyang Li* 1 , Kevin Stone1 , and Ajit Vikram†1 1 Pharmaceutical Analysis & Digital Technologies, Merck & Co., Inc., Rahway, New Jersey, 07065, United States
A BSTRACT Pharmaceutical process characterization is central to Quality by Design because it defines how variations in process parameters affect the ability to meet product quality specifications, thereby supporting proven acceptable ranges and robust manufacturing. In practice, however, characterization still relies largely on factorial design of experiments (DOE) approaches, which are inefficient for resolving multivariate pass/fail boundaries in higher-dimensional spaces. While Bayesian optimization has transformed process optimization, adaptive methods for multi-objective process characterization remain lacking. Here, we introduce JAREX (Joint Acceptable Region EXploration), a Bayesian active-learning acquisition function for multi-objective process characterization. JAREX formulates characterization as a joint boundary-learning problem and adaptively selects experiments to recover the joint pass region defined by simultaneous satisfaction of threshold criteria across multiple objectives. JAREX combines an optimistic joint-feasibility mask with a multi-objective extension of randomized straddle, focusing sampling on the joint edge of failure. Our benchmark study suggests that JAREX provides more accurate and sample-efficient recovery of the joint pass region than factorial DOE, space-filling designs, and greedy objective-wise strategies over the full experimental budget range. For batched experimentation, it reduces the number of iterative process characterization experiments by more than half while preserving high accuracy for the boundary-identification task. Implemented in the open-source obsidian package, JAREX provides a modular framework for adaptive, data-efficient multi-objective algorithmic process characterization, supporting sampleefficient range finding in high-dimensional spaces.
1
Introduction
Pharmaceutical manufacturing under the Quality-by-Design (QbD) paradigm depends on a quantitative understanding of how process parameters (e.g., temperature, pressure, time, solvent ratio, etc.) impact in-process attributes and critical quality attributes (CQAs) such as identity, purity, and yield. [1–8] Two complementary approaches build this understanding: process optimization and process characterization. Process optimization identifies optimal operating conditions that maximize desirable outcomes while minimizing undesirable ones. Process characterization, however, is not about identifying an optimum but about measuring the operating ranges of the process: how CQAs respond across the parameter space and identifying where the edge of failure lies within relevant operational windows, beyond which specifications can no longer be met. In current practice, process characterization most commonly results in establishing proven acceptable ranges (PARs), informed by univariate windows varying one factor at a time (OFAT). [1, 3, 9] PARs alone, however, do not constitute a design space, the multidimensional region demonstrated to deliver quality assurance against joint parameter excursions rather than OFAT variations. [1, 2, 10] Exploring a design space therefore requires multivariate characterization of the parameter space, which is costly and data-intensive in practice, and this is the problem this work addresses. Traditionally, process characterization workflows rely heavily on factorial design of experiments (DOE) or other statistical screening designs. [8, 9, 11, 12] These methods provide systematic but non-adaptive sampling of the parameter ∗ †
[email protected] [email protected]
JAREX for Multi-Objective Process Characterization
A P REPRINT
space and scale poorly with dimensionality. Alternatively, adaptive, data-driven strategies that iteratively select informative experiments based on prior observations may be useful to explore design space efficiently. A useful precedent comes from process optimization, where Bayesian optimization (BO) has emerged as a transformative tool. [13, 14] BO combines probabilistic surrogate models, typically Gaussian process (GP) models, with acquisition functions such as Expected Improvement, [15] Upper Confidence Bound, [16, 17] and Probability of Improvement [18] to balance exploration and exploitation. Open-source frameworks such as BoTorch [19] have made these methods broadly accessible. An analogous adaptive paradigm is needed for process characterization, but the gap between theory and practice remains substantial. While mathematical frameworks for this type of problem do exist, [20–23] their translation into practical, application-ready tools for pharmaceutical process characterization is still limited. Moreover, to our knowledge, no acquisition function exists for multi-objective process characterization (MOPC), despite nearly all pharmaceutical processes involving multiple CQAs simultaneously. Here, we present JAREX (Joint Acceptable Region EXploration), a new acquisition function designed for characterizing high-dimensional processes with multiple objectives. JAREX adaptively targets the edge of failure and accounts for correlations among objectives, while standard outputs such as PARs and multi-factor interaction plots are recovered from the fitted surrogate via post-processing. On a simulated kinetic model, JAREX substantially outperforms other methods in sample efficiency and boundary recovery. The method is available in the open-source obsidian package [24] and brings to pharmaceutical process characterization the kind of adaptive data-driven guidance that Bayesian optimization has brought to process optimization.
2
Methodology
Process characterization is fundamentally a global learning task: the goal is to understand how quality attributes respond across an entire parameter space, with particular emphasis on where they cross their acceptable quality limits. This differs sharply from optimization, which seeks a single optimal point (or, in multi-objective settings, a Pareto surface). Consider a process with parameter space X ⊆ Rd , where d is the number of process parameters, and let x ∈ X denote a specific set of operating conditions. For each of m quality attributes, we define an objective Oi (x) for i = 1, . . . , m, together with a threshold hi that encodes the pass/fail specification. Without loss of generality, we take “pass” to mean Oi (x) ≥ hi , so the pass region for attribute i is Pi = {x ∈ X : Oi (x) ≥ hi }.
(1)
Our characterization task is to recover these pass regions from a surrogate model fitted on a finite sample budget. There are two natural ways to frame this task algorithmically, and they motivate distinct acquisition strategies. The first is a feasibility perspective. At every candidate x, query the likelihood of the objective passing, and sample where this likelihood is most informative. Ierapetritou and coworkers [23] took this route by adapting Expected Improvement into a feasibility-oriented acquisition function (feasibility EI) that favors candidates near the threshold with high model uncertainty. The full definition and interpretation are given in Section S1.1 of the Supplementary Information. The second is a boundary perspective: rather than reasoning pointwise about feasibility, treat the surface {x : Oi (x) = hi } itself, the “edge of failure” in the ICH Q8(R2) guideline [1], as the object to recover. This is precisely the level set estimation (LSE) problem studied in applied mathematics. Two main algorithm families have emerged for LSE: stepwise uncertainty reduction, [25, 26] which reduces global uncertainty about the level set, and straddle-based methods, [20, 21, 27] which target points that are simultaneously near the threshold and uncertain. Our method builds on the boundary perspective via the straddle family, for reasons that will become clear once we introduce the multi-objective setting. 2.1
Joint Pass Region: the Target of Multi-Objective Characterization
Nearly all real pharmaceutical processes involve multiple quality attributes, and characterizing each objective independently ignores the correlations that ultimately determine whether a set of conditions is acceptable. The relevant target is the joint pass region, the intersection of all individual pass regions, Pjoint =
m \
Pi = {x ∈ X : Oi (x) ≥ hi , ∀i}.
(2)
i=1
Outside Pjoint , at least one attribute fails its specification, and the identity of the failing objective is often the critical piece of information for process understanding. Pjoint makes precise the ICH Q8 concepts already introduced in Section 1: its boundary is exactly the edge of failure, beyond which CQAs cannot be met. In principle, a fully converged estimate of Pjoint would also describe the largest 2
JAREX for Multi-Objective Process Characterization
A P REPRINT
possible design space, since every point inside satisfies every quality specification. In practice, however, a regulatorapproved design space is almost never drawn to coincide with our estimate of Pjoint . Real processes are subject to parameter variability and measurement noise, and any estimate of Pjoint obtained from finite data carries residual model uncertainty. The design space is therefore deliberately chosen as a conservative subset strictly inside Pjoint , with margins that reflect both process variability and the confidence level of the underlying characterization. ICH Q8 explicitly notes that determining the edge of failure is not essential for establishing a design space. [1] However, characterizing it efficiently is what makes better-justified design space choices possible, and this is exactly what JAREX is built to do. We therefore adopt Pjoint as the mathematical target: it is the largest region consistent with all specifications, the upper bound against which any practical design space is measured. 2.2
Joint Acceptable Region EXploration (JAREX)
With the target Pjoint established, we now construct an acquisition function that directly learns its boundary, the edge of failure. We build up from two well-understood single-objective components, upper confidence bound (UCB), and the randomized straddle algorithm, then extend to multiple objectives. The classic UCB acquisition function [16, 17] selects candidates by maximizing UCB(x) = µ(x) + β 1/2 · σ(x),
(3)
where µ(x) and σ(x) are the predicted mean and standard deviation from the surrogate, and β 1/2 tunes the balance between exploitation (high mean) and exploration (high uncertainty). UCB targets extrema, not boundaries, but its exploitation-exploration trade-off structure is the scaffold for what follows. The straddle algorithm of Bryan et al. [20], later formalized by Gotovos et al. [21], reshapes UCB for level set estimation by optimizing against the threshold instead of toward an extremum, STR(x) = −|µ(x) − h| + β 1/2 · σ(x).
(4)
The exploitation term µ(x) in UCB is replaced by −|µ(x) − h|, which penalizes predictions far from the threshold. Maximizing STR(x) therefore prefers candidates that are either close to the boundary (small |µ(x) − h|) or highly uncertain (large σ(x)), exactly the two kinds of points that are most informative for learning the boundary. The main practical difficulty with Eq. 4 is choosing β 1/2 and the conventional value β 1/2 = 1.96 limits efficient boundary exploration in practice. To address this, Inatsu et al. [27] proposed a randomizing straddle variant that samples β at every iteration, STR∗ (x) = max [STR(x), 0] ,
(5)
with β drawn fresh each step from a χ22 distribution, β ∼ χ22 (ξ) = 21 e−ξ/2 ,
ξ ≥ 0,
(6)
which adaptively balances aggressive boundary refinement with occasional broader exploration. More details are provided in Section S1.2 of the Supplementary Information. The randomized straddle algorithm STR∗ gives us a well-tuned way to learn a single boundary. MOPC, however, requires more than a per-objective treatment. Running STR∗i independently per objective i (a greedy strategy) ignores correlations and spends effort refining each individual boundary rather than the joint boundary that defines Pjoint . Joint Acceptable Region EXploration (JAREX) addresses both issues at once by restricting the search to regions that are plausibly within Pjoint and then choosing the objective whose boundary is most informative there. The first ingredient is an optimistic view of each individual pass region. With surrogate models providing µi (x) and σi (x), we cannot apply Eq. 1 directly: a point near the boundary with large σi is exactly where sampling is most informative and must not be prematurely excluded. Borrowing directly from UCB, we define the per-objective 1/2 optimistic pass region using the Upper Confidence Bound UCBi (x) = µi (x) + βi σi (x) (Eq. 3) as P̂i = {x ∈ X : UCBi (x) ≥ hi },
(7)
and the joint optimistic pass region as their intersection, P̂joint =
m \
P̂i = {x : UCBi (x) ≥ hi , ∀i}.
i=1
3
(8)
JAREX for Multi-Objective Process Characterization
A P REPRINT
P̂joint contains every candidate that is either confidently inside Pjoint or too uncertain to rule out. It acts as a mask: only candidates x̂ ∈ P̂joint are eligible for selection. Within this region, we compute the randomized-straddle score per objective, h i 1/2 STR∗i (x̂) = max −|µi (x̂) − hi | + βi · σi (x̂), 0 ,
(9)
which quantifies how informative sampling x̂ would be for the i-th boundary. We then combine the per-objective scores with a softmin-weighted sum, Ajarex (x̂) =
m X
softmini [STR∗i (x̂)] · STR∗i (x̂)
(10)
i=1 ∗
e−STRi (x̂)/τ softmini = P −STR∗ (x̂)/τ , j je with τ = 0.5 by default. The softmin-weighted sum smoothly interpolates between the mean (τ → ∞) and the min (τ → 0). It rewards points where multiple objectives are simultaneously informative for their boundaries and suppresses points where only one objective is critical while the others are confidently far from their thresholds. This AND-aggregation aligns with the joint pass region being an intersection (Eq. 8): the joint boundary is, by definition, where several objectives co-approach their thresholds. A detailed comparison with max, mean, fixed weighted-sum, and softmax combinations is given in Section S1.3 of the Supplementary Information. The next experiment is then selected via the familiar x∗ = arg max [Ajarex (x̂)] , x̂
(11)
the surrogate is updated with the new observation, and the procedure repeats. One implementation detail is worth addressing here, as it materially affects performance. The joint optimistic pass region P̂joint is an irregular, generally non-convex subset of X and cannot be written as a simple closed-form constraint. A hard indicator mask creates a sharp cliff at its boundary with a flat zero plain outside, which blocks the gradientbased multi-start search, a strategy widely used by common optimizers to refine candidate points. [19] Section S1.4 of the Supplementary Information gives the full discussion. Instead, we use a differentiable soft mask. For each x, define the worst-case slack across objectives, d(x) = min[UCBi (x) − hi ], i
(12)
so that d ≥ 0 iff x ∈ P̂joint . The soft mask is then M (x) =
1 d≥0 , sech(kd) d < 0
(13)
where k = 1 by default and sech(x) = 2/(ex +e−x ). M is identically 1 throughout P̂joint and decays smoothly outside, preserving the intended acceptance region while remaining differentiable everywhere. The softmin combination and the soft mask are the key components for extending classical randomized straddle to JAREX. These two modifications break a specific step in the convergence proof of Ref. 27. As a result, the single-objective guarantee does not carry through to the joint misclassification loss. A much weaker statement can be recovered, but only by assuming that every objective is queried sufficiently often. This assumption contradicts JAREX’s explicit deprioritization of regions confidently outside P̂joint . We therefore do not claim a formal bound here and defer a full analysis to future work. Section S1.5 of the Supplementary Information identifies exactly which step of the argument fails and why. The empirical results in Section 3 demonstrate that this is an acceptable trade-off for the multi-objective setting. Algorithm 1 summarizes the full procedure as implemented in our open-source package obsidian, [24] a library for algorithmic process design for pharmaceutical applications built on BoTorch [19] and PyTorch. [28] The JAREX object inherits from BoTorch’s MCAcquisitionFunction class, so it can be readily reused in other packages of similar design.
3
Results and Discussion
In this section, we demonstrate the effectiveness of the randomized straddle and JAREX algorithms for process characterization. We begin by establishing benchmark metrics and illustrating the boundary-focused selection strategy of 4
JAREX for Multi-Objective Process Characterization
A P REPRINT
Algorithm 1 JAREX: Joint Acceptable Region EXploration Require: Initial dataset D0 = {(xj , yi,j )} for m objectives Require: Thresholds {hi }m i=1 , temperature τ , decay k 1: for t = 1, 2, . . . , T do 2: Train surrogate models: 3: for each objective i = 1, . . . , m do 4: Fit GP model to obtain µi (x) and σi (x) 5: end for 6: Sample βi ∼ χ22 for i = 1, . . . , m 7: Compute JAREX acquisition for all candidates x: 8: for each objective i = 1, . . . , m do 1/2 9: STR∗i (x) = max[−|µi (x) − hi | + βi σi (x), 0] 10: end for exp(−STR∗ i /τ ) 11: softmini = P exp(−STR ∗ /τ ) j Pjm 12: Ajarex (x) = i=1 softmini [STR∗i (x)] · STR∗i (x) 13: Compute a soft mask for all candidates x: 14: d(x) = mini [UCBi (x) − hi ] 1 d≥0 15: M (x) = sech(kd) d < 0 16: Select next point: 17: x∗ = arg maxx [M (x) · Ajarex (x)] 18: Evaluate yi,t = Oi (x∗ ) for all objectives 19: Update Dt = Dt−1 ∪ {(x∗ , {yi,t })} 20: end for
▷ Eq. 6
▷ Eq. 9
▷ Eq. 10 ▷ Eq. 12 ▷ Eq. 13 ▷ Eq. 11
randomized straddle. We then present benchmark results for single-objective characterization using randomized straddle and multi-objective characterization using JAREX on a simulated chemical kinetics problem. Our results show that JAREX achieves superior joint characterization performance compared to traditional methods (factorial DOE and space-filling) and greedy multi-objective approaches, while maintaining robust performance across different transformation strategies. Finally, we demonstrate that surrogate models trained with JAREX can recover traditional process characterization outputs, including proven acceptable ranges and multi-factor interaction plots.
3.1
Benchmark Details
We evaluate method performance by comparing the identified pass regions with the true pass regions (known ground truth) using the Jaccard index, [29] J(Ppred , Ptrue ) =
|Ppred ∩ Ptrue | , |Ppred ∪ Ptrue |
(14)
where Ppred and Ptrue are the predicted and true pass regions, respectively. The Jaccard index ranges from 0 (no overlap) to 1 (perfect match), providing a set-based metric focused on boundary classification rather than prediction error magnitude. The Jaccard index is an appropriate metric for process characterization because correct boundary identification inherently requires accurate predictions near the threshold, where it matters most (SI Section S2 includes an extended discussion of this index). For all benchmark studies presented here, we use response surface methodology (RSM) designs and space-filling (SF) sampling as baseline approaches. The RSM designs comprise a full two-level factorial and a central composite design (CCD). Each test function is evaluated with 10 independent campaigns using fixed random seeds across methods. All methods share the same initial observations generated via Latin hypercube sampling (LHS), with campaign budgets and implementation details provided in the Supplementary Information. Gaussian process surrogate models with Matérn kernels (BoTorch implementation) are used throughout, with hyperparameters optimized via multi-start optimization (see SI Section S2 for optimization details). Jaccard indices are calculated on Sobol-discretized parameter spaces, with evaluation point density adapted to surrogate uncertainty. 5
JAREX for Multi-Objective Process Characterization
95% CI
True Function
Predicted Mean
Observed Points
Acquisition
1.0
0 (a) Step 5 New trial: y = −0.49 Threshold: h = −0.31
−1 −2
A P REPRINT
0.5
−4
0.0 1.0
0 (b) Step 8 New trial: y = −0.30 Threshold: h = −0.31
y
−1 −2
0.5
−3 −4
0.0 1.0
0 (c) Step 15 New trial: y = −0.31 Threshold: h = −0.31
−1 −2
Acquisition Function (a.u.)
−3
0.5
−3 −4
0.0 −3
−2
−1
0
1
2
3
Jaccard Index
x 1.0
(d)
0.5 0.0 0.0
2.5
5.0
7.5
10.0
12.5
15.0
17.5
20.0
Characterization Step
Fig. 1: Next experiment selected by the randomized straddle algorithm at steps 5 (panel a), 8 (panel b), and 15 (panel c) of a one-dimensional characterization campaign. Blue curve: true function; grey dashed line: threshold; purple curve and shaded area: GP mean and uncertainty; red dots: existing observations; black square: next selected experiment. The black curve shows the acquisition function value (arbitrary scale), with higher values indicating stronger sampling preference. Panel d: characterization history across iterations.
3.2
Selection Rules in Randomized Straddle
To select samples efficiently in an iterative setting, the randomized straddle algorithm (Eq. 5) prefers candidates that are either near the threshold or have high uncertainty. Fig. 1 illustrates this behavior on a one-dimensional doublewell potential test case (see SI Section S3.1.1 for function definition), showing the next experiments at different campaign steps (panels a-c) and the overall characterization history (panel d). The algorithm selects candidates with high uncertainty early in the campaign (e.g., step 5) and focuses on the proximity of the threshold as uncertainty decreases (e.g., step 8). By step 15, the algorithm continues to sample exactly at the threshold to refine the boundary, which causes the Jaccard index to plateau as further improvements become marginal. Importantly, the algorithm ignores regions far from the boundary even when the surrogate mean deviates from the true function, as these regions do not affect pass/fail classification. This boundary-focused sampling strategy exemplifies the fundamental difference between characterization (exploring pass/fail boundaries) and optimization (finding optimal points). 3.3
Single-Objective Characterization
Fig. 2 shows the benchmark results of applying the above-mentioned randomized straddle algorithm to a singleobjective case study. We use noise-free response functions here (results with 5% Gaussian noise are provided in the Supplementary Information) so that the comparison isolates each method’s intrinsic sampling strategy at modest budgets, rather than its tolerance to observation noise. This choice does not amplify any single method. The results reveal several key insights. First, UCB, which is designed for optimization tasks rather than characterization ones, performs poorly beyond 2D despite high exploration (β 1/2 = 6). This confirms that boundary exploration requires fundamentally different acquisition strategies than optimum-seeking. Second, the RSM designs and SF provide competitive performance in low-dimensional settings but fall off sharply as the dimensionality and complexity increases, highlighting the curse of dimensionality inherent in uninformed sampling. Third, among characterizationspecific methods, randomized straddle shows comparable performance to feasibility EI in simple cases but demon6
JAREX for Multi-Objective Process Characterization
(a) Six-Hump Camel (2D)
2-level (5)
1.00
Straddle
0.75
CCD (9)
UCB (β = 6) Feasibility EI
Space Filling UCB (β = 1.96)
A P REPRINT
0.60
0.50 0.25
5
29
41
53
65
100
130
160
2-level (17)
(b) Hartmann 4D
0.75 0.50 0.25
0.00 0.01
0.00 10
70
(c) Hartmann 6D
2-level (65)
1.00
40
0.75
CCD (77)
Jaccard Index
1.00
17
CCD (25)
0.01
0.00
0.50 0.25
0.25
0.12
0.00 15
75
135
195
255
Number of Observations
315
DOE
Fig. 2: Benchmark results for single-objective characterization on noise-free analytical test functions: six-hump camel (2D), Hartmann 4D (4D), and Hartmann 6D (6D). Jaccard index convergence (mean and standard deviation averaged across 10 campaigns) for SF, UCB (β 1/2 = 1.96 and β 1/2 = 6), feasibility EI, and randomized straddle. The side panel (DOE) reports the RSM designs, a full two-level factorial and a central composite design (CCD), as singledesign Jaccard values, with run counts in parentheses. strates superior convergence speed and numerical stability in high-dimensional settings (6D), where it substantially outperforms all other approaches. 3.4
Multi-Objective Characterization
Next, we benchmark the developed acquisition function, JAREX, on a multi-objective task using a simulated reaction kinetic model. Specifically, we compare how JAREX performs relative to RSM designs, SF, and a greedy multiobjective extension of randomized straddle (which is designed for single-objective problems). The greedy approach cycles through objectives one at a time. At each step, randomized straddle selects the next experiment based on a single objective, while all objectives are measured and used to train the surrogate models. The joint pass region is obtained by intersecting individual pass regions. While this approach seems plausible, it does not explicitly target the joint boundary where multiple objectives simultaneously approach their thresholds. We developed a four-dimensional kinetic model simulating competing reaction pathways in a chemical process. The reaction scheme consists of a main catalytic pathway producing the desired product P, alongside a competing reversible side reaction forming impurity X: kmain (T, Cat)
Y + Z −−−−−−−→ P kside,fwd (Cat)
−− 2Y + Z ↽ −− −− −− −− −− −⇀ −X kside,rev (T )
(15) (16)
where Z is the reactant, Y is the co-reactant, P is the desired product, X is an undesired impurity, and Cat represents the catalyst. The parameter space is defined by four input parameters: initial reactant concentration ([Z0 ]), initial co-reactant concentration ([Y0 ]), temperature (T ), and catalyst loading ([Cat0 ]). Together, these parameters create coupled tradeoffs among reactant Z conversion, product purity, and cost per unit product (Fig. 3). Detailed reaction mechanisms and mechanistic interpretation are provided in Section S4.1 of the Supplementary Information. 7
JAREX for Multi-Objective Process Characterization
Experimental Cost
Reactant Conversion
100
200
500 1000
0.6
0.8
A P REPRINT
Cost per Unit Product 100
1000
10000 100000
1.0
0.8
0.6
0.4 0.4
1.0 0.4 0.6 Product Purity
0.8
1.0
Fig. 3: Trade-offs among kinetic model objectives. Reactant Z conversion vs. product purity colored by raw experimental cost in panel (a) and cost per unit product in panel (b). We evaluate three competing objectives: conversion (maximize), product purity (maximize), and cost per unit product (minimize). Conversion is the fraction of initial reactant Z0 consumed. Purity is the molar ratio [P ]/([P ] + [X]) at the final time point. Cost is defined per mole of desired product P , accounting for reactant Z, co-reactant Y , catalyst, and solvent consumption (temperature affects cost indirectly through kinetics and final yield). These objectives exhibit inherent trade-offs, as shown in Fig. 3. For instance, high conversion requires high [Z0 ] and [Y0 ], while high purity requires suppression of the side reaction, which often reduces conversion. Conversely, the high-conversion, high-purity regime (upper right in Fig. 3) incurs the highest raw material cost yet yields the lowest cost per unit product due to superior process efficiency. The objectives are physically coupled through the underlying kinetics with trade-offs and reinforcements that shift across the parameter space. This non-trivial interplay makes the kinetic model a meaningful benchmark for multi-objective characterization methods and more representative of real industrial problems. Table 1 lists the parameter space and objective thresholds. See Supplementary Information for complete kinetic equations, rate expressions, and integration details. Table 1: Parameter space and objective thresholds for the kinetic model. Parameter space
Objectives
Variable
Range
Objective
Threshold
[Z0 ] [Y0 ] [Cat0 ] T
[0.8, 3.0] M [0.8, 6.0] M [0.005, 0.2] M [273, 333] K
Z conversion Purity Cost per unit product
≥ 90% ≥ 99.3% ≤ $570
Fig. 4 presents benchmark results for the kinetic model. All methods use a standard transformation of the response function (see Section S4.2 for log-transformed results and transformation analysis). Fig. 4a shows joint Jaccard convergence on this benchmark. JAREX (black curve) delivers the strongest performance across the entire budget range: the joint Jaccard rises rapidly within the first 30–40 iterations, continues to improve smoothly thereafter with tight error bars across the 10 campaigns, and settles well above every competing approach. This confirms that explicitly targeting the joint edge of failure (Eq. 8) yields efficient characterization of the joint pass region even in this four-dimensional, three-objective, noisy setting. The non-adaptive baseline falls well short of JAREX. The RSM designs perform poorly despite being conventional goto methods for process characterization, illustrating that grid-based designs scale poorly to realistic parameter spaces. SF improves on the RSM designs through better coverage but still plateaus far below JAREX, because uninformed sampling cannot redirect effort toward the joint boundary as observations accrue. The greedy multi-objective extension of randomized straddle, although adaptive, suffers from two critical flaws. First, performance becomes highly sensitive to objective ordering, producing the erratic convergence histories shown for different orderings; the circular dependency is fundamental, because choosing the “right” objective to prioritize requires prior knowledge of the response surface, but characterization exists to acquire that knowledge. Second, greedy methods refine each objective boundary separately without considering their interplay, learning individual pass regions efficiently but failing to characterize the joint pass region that actually defines acceptable operating space. JAREX avoids 8
JAREX for Multi-Objective Process Characterization
0.8 0.6
0.35
0.4 SF JAREX Greedy: Cv-Pu-UC
0.2 0.0 8
Jaccard Index
CCD (27)
2-level (19)
Jaccard Index
1.0 (a)
A P REPRINT
32
56
Greedy: Pu-Cv-UC Greedy: UC-Cv-Pu 0.00
80
104
128
DOE Number of Observations 1.0 (b) JAREX (c) Greedy: Cv-Pu-UC Purity Conversion Cost Joint
0.5
0.0 0
25
50
75
100
0
25
50
75
100
Jaccard Index
Number of Observations 1.0 (d) 0.5
Batch size=1 Batch size=4 Batch size=8
0.0 0
24
48
72
96
120
Number of Iterations
Fig. 4: Multi-objective characterization benchmarks for reactant Z conversion, product purity, and cost per unit product with 5% Gaussian noise added to the simulation. (a) Joint Jaccard index convergence for SF, greedy randomized straddle (multiple objective orderings), and JAREX; the side panel (DOE) reports the RSM designs—a full two-level factorial and a CCD—as single-design Jaccard values, with run counts in parentheses. (b-c) Individual objective and joint Jaccard indices for JAREX and greedy, respectively. (d) JAREX batching comparison: serial (q = 1), moderate (q = 4), and high (q = 8) batch sizes. Curves show average Jaccard indices across 10 campaigns. Error bars indicate ±1σ.
both pitfalls by explicitly targeting the joint boundary (the edge of failure). This yields stable, order-independent convergence that continues improving even after greedy and SF methods plateau. The individual-versus-joint performance trade-off is illustrated in Fig. 4b-c, where panel b shows JAREX and panel c shows the greedy results. JAREX trades off individual objective performance to achieve superior joint characterization performance by deprioritizing high-uncertainty regions known to fail for at least one objective (Eq. 8). This strategy focuses the experimental budget on the joint pass region, which is what matters for process understanding and decisionmaking. The greedy method, by contrast, achieves high Jaccard indices for purity and conversion but performs notably worse for cost. Interestingly, the joint Jaccard index exceeds the worst individual objective (cost). This occurs because the three objectives are physically coupled through product yield Pf . The true joint pass region is 2.4 times larger than expected under independence assumptions. The cost pass set includes a large region failing purity or conversion. While these points inflate the cost surrogate’s classification errors, the joint classifier correctly rejects them, improving joint precision. A detailed quantitative breakdown of these effects is provided in the Supplementary Information (Section S4.3). In practical experimental settings, real-time campaign duration is sometimes more constraining than total sample budget especially when development timelines are more precious than material cost. Batching, where multiple experiments are run in parallel, is a critical strategy for reducing calendar time, especially when individual experiments require hours or days. However, batching naturally poses a challenge for the efficiency of BO-based methods. Without the ability to update the surrogate between experiments within a batch, the algorithm must commit to multiple points simultaneously, potentially leading to redundant sampling or suboptimal exploration. Fig. 4d evaluates JAREX under three batch scenarios: serial (q = 1, 120 iterations), moderate batching (q = 4, 80 iterations/320 total samples), and high batching (q = 8, 50 iterations/400 total samples). The results reveal a favorable trade-off. While batched campaigns consume more total samples, they achieve superior accuracy with dramatically fewer iterations. Specifically, q = 4 exceeds serial Jaccard after only 50 iterations and q = 8 after only 40, each cutting real-time campaign duration by more than half, while improving final Jaccard index through the larger total sample budget. The larger total sample budget further improves the final Jaccard index. This demonstrates that JAREX maintains strong performance 9
JAREX for Multi-Objective Process Characterization
A P REPRINT
under batching constraints, making it practical for time-sensitive experimental workflows where reducing the number of "wait-and-decide" cycles is as important as minimizing total experimental burden. 3.5
Application to Standard Process Characterization Metrics
Traditional process characterization studies are interpreted through metrics such as proven acceptable ranges (PARs) and multi-factor interactions. Having demonstrated the effectiveness of JAREX for efficient characterization of the joint pass region, we next evaluated its utility for recovering example process characterization outputs. Although randomized straddle and JAREX are not formulated to directly target these conventional metrics, they are designed to learn the underlying pass/fail boundary. The resulting surrogate therefore provides a high-dimensional representation of process robustness, from which traditional characterization metrics can be recovered as lower-dimensional projections through post-processing. Because this is a broader learning problem than dedicated univariate or bivariate characterization, accuracy for any single derived metric at a fixed experimental budget may be somewhat lower than that of approaches optimized specifically for that metric. We ran a JAREX campaign on the kinetic model starting from 8 initial observations. JAREX then iteratively designed 32 additional experiments by maximizing its acquisition function (Eq. 10), with the GP surrogate retrained after each new observation. The resulting trained surrogate is the input to all PAR analyses below. The JAREX-derived PARs are obtained directly from the trained posterior, evaluated along univariate scans through the centroid of the identified joint pass region. Ground Truth
JAREX
Threshold
1.00 Product 0.75 Purity
0.50 1.00 Reactant 0.75 Conversion
0.50 900 Cost per 300 Unit Product
100
Joint Acceptable Range
1.2 1.8 2.4 [Z0] (mol/L)
1.5 3 4.5 [Y0] (mol/L)
285 300 315 330 0.05 0.1 0.15 T (K) [Cat0] (mol/L)
Fig. 5: Comparison of proven acceptable ranges (PAR) for reactant Z conversion, product purity, and cost per unit product between JAREX and traditional linear search. Top three rows: individual objective pass regions (green shaded areas). Bottom row: joint PAR shown as horizontal bars, with the ground truth in blue and the JAREX prediction in orange. The brackets mark the 70% confidence bounds. JAREX predictions closely track ground truth within and near the joint pass region, while showing larger deviations outside, reflecting its strategy of focusing experimental budget on the joint pass region. Fig. 5 compares JAREX-derived PARs with traditional linear-search PARs for all three objectives and for the joint region. In this case, purity is the binding constraint for the joint pass region, while conversion and cost have broader feasible ranges. JAREX tracks the ground truth closely within and near the joint boundary across all dimensions. Outside that region, individual-objective uncertainty increases, which is expected because JAREX intentionally deprioritizes areas already likely to fail at least one objective. This allocation concentrates experimental effort where it most affects joint characterization. Even so, the 70% confidence bounds remain consistent with a coherent joint evaluation across all dimensions. Two-dimensional interaction effects between pairs of process parameters are commonly used in process characterization to understand the interplay between parameters and identify potential synergies or trade-offs. Traditionally, these interaction maps require dedicated two-factor campaigns at multiple fixed settings of the remaining parameters, which adds experimental cost and time to the process characterization study. With a trained JAREX-based model, they can be obtained directly by evaluating the posterior on two-parameter grids while fixing the remaining parameters, with no additional experiments. Fig. 6 shows the interaction effects between the initial reactant concentration [Z0 ] and initial catalyst concentration [Cat0 ] obtained from JAREX with other parameters fixed at various values. From top to bottom, the joint pass region transitions from nearly all passing to almost all failing. This variation is primarily driven by the purity constraint, 10
T=300 T=305
[Cat0]
0.15 0.1 0.05 0.15 0.1 0.05
T=310
JAREX for Multi-Objective Process Characterization
0.15 0.1 0.05 1.5 2 2.5 [Y0]=2.0
Fail
Mean pass
1.5 2 2.5 [Y0]=3.0
1.5 2 2.5 [Y0]=2.0
70% confidence
1.5 2 2.5 [Y0]=3.0
1.5 2 2.5 [Y0]=2.0
1.5 2 2.5 [Y0]=3.0
A P REPRINT
95% confidence
1.5 2 2.5 [Y0]=2.0
1.5 2 2.5 [Y0]=3.0
[Z0] Product Purity
Reactant Conversion
Cost per Unit Product
Joint Evaluation
Fig. 6: Two-dimensional interaction effects for reactant Z conversion, product purity, and cost per unit product between initial reactant concentration [Z0 ] and catalyst concentration [Cat0 ] at various fixed values of other parameters, obtained from a JAREX model trained with 40 samples. Confidence levels: mean pass (yellow), 70% (light green), 95% (dark green), fail (red). Last column: joint evaluation with confidence determined by intersection of individual pass regions.
which responds sensitively to the other parameters (temperature T and initial co-reactant concentration [Y0 ]) held at different fixed values. Such sensitivity highlights the importance of exploring interactions at multiple operating conditions rather than assuming a single interaction map is representative. As shown in Fig. 6, JAREX focuses sampling on the joint boundary, where experiments are most informative for defining the feasible operating region. This targeted allocation of the experimental budget yields high-confidence characterization of the joint pass region while maintaining consistent boundary structure as the budget increases (Supplementary Information, Section S4.4). These analyses and visualizations are also available through the modules in the obsidian package. [24] Moreover, JAREX’s iterative framework naturally supports targeted follow-up experiments, selecting new conditions that maximize the information gain and efficiently refine any remaining uncertainty where it matters most.
Conclusions In this paper, we introduce JAREX (Joint Acceptable Region EXploration), an acquisition function specifically designed for multi-objective process characterization using a Bayesian active-learning framework. Methodologically, JAREX formulates process characterization as a joint boundary-learning problem. Rather than optimizing a single response or characterizing each objective independently, it directly targets the joint pass region defined by simultaneous satisfaction of all quality specifications. This aligns the acquisition rule with the actual objective of Quality-by-Designdriven characterization, namely, understanding the multidimensional edge of failure that governs acceptable process operation. By combining an optimistic joint-feasibility mask with a randomized-straddle strategy, JAREX focuses on proposing experiments on the most informative parts of the joint boundary while accounting for cross-objective interactions. Notably, we report JAREX as a modular open-source framework, implemented as part of a Python-based package (obsidian). Moreover, our benchmark studies demonstrate that the JAREX approach improves learning efficiency for the process characterization task. In the single-objective formulation, the underlying randomized-straddle strategy outperforms alternative approaches, including factorial DOE, space-filling methods, and acquisition functions such as Upper Confidence Bound and feasibility Expected Improvement. On a four-dimensional, three-objective kinetic model benchmark, JAREX outperforms the greedy multi-objective extension of randomized straddle across the entire budget range. This indicates that adaptive sampling directed at the joint edge of failure is more resource efficient (fewer experiments needed) than uniform coverage or sequential refinement of separate boundaries. Under batched iterative experimentation (i.e., multiple experiments per iteration), JAREX further halves the number of iterations needed while matching the accuracy of serial runs (i.e., one experiment per iteration). A practical strength of the framework presented here is its compatibility with the expected outputs from the traditional process characterization workflow. Although JAREX is designed to learn the full multidimensional joint boundary, the trained surrogate can be leveraged to recover metrics such as proven acceptable ranges and critical multi-factor process interactions. This makes the method directly relevant to existing pharmaceutical development workflows while providing a more data-efficient route to process understanding. More broadly, JAREX and the accompanying open-source framework aim to bring to process characterization the same shift that Bayesian optimization enabled for process optimization: replacing largely static, uninformed experi11
JAREX for Multi-Objective Process Characterization
A P REPRINT
mentation with adaptive, model-driven decision-making. In that sense, this work establishes a practical foundation for data-efficient, multivariate characterization of pharmaceutical processes. Future efforts should focus on broader experimental validation across diverse unit operations, principled handling of correlated and heteroscedastic measurement noise, and integration with closed-loop robotic platforms to enable fully autonomous characterization workflows.
Author contributions XL derived and implemented the method and performed the benchmarking experiments. AV and KS proposed and provided guidance on the project. All authors contributed to manuscript revision.
Conflicts of interest The authors declare that they have no conflict of interest.
Data availability The single- and multi-objective characterization functions are available as part of the open-source obsidian package. The benchmarking scripts and all supporting data (benchmark campaigns, kinetic model scan, and process-characterization data, as JSON/CSV with figure-reproduction scripts) are openly available at Zenodo, DOI: 10.5281/zenodo.21923038.
Acknowledgements XL acknowledges valuable discussions with Eugene Zakharov and Yingjie Chen at Merck & Co., Inc., Rahway, NJ, USA.
12
JAREX for Multi-Objective Process Characterization
A P REPRINT
Supplementary Information JAREX: An Acquisition Function for Multi-Objective Algorithmic Process Characterization
S1
Methodology Details
This section collects technical details that support the Methodology section in the main text but are not essential to the narrative. S1.1
Feasibility Expected Improvement
The main text mentions the Feasibility Expected Improvement (feasibility EI) acquisition function of Ierapetritou and coworkers [23] only conceptually. Its explicit form is µ(x) − h EIfeas (x) = σ(x) · ϕ , (S1) σ(x) where µ(x) and σ(x) are the predicted mean and standard deviation at x, h is the threshold, and ϕ(·) is the standard normal probability density function. Conceptually, feasibility EI measures how likely a candidate is to lie near the threshold boundary: it is maximized when µ(x) is close to h or when σ(x) is large. Sampling where feasibility EI is large therefore concentrates experiments both near the feasibility boundary and in regions of high uncertainty, efficiently classifying the parameter space into feasible (µ(x) > h) and fail (µ(x) < h) regions. Although the original authors do not use the term “characterization”, the goal of classifying the parameter space into pass/fail regions is the same. S1.2
The χ22 Distribution
The central difficulty with the deterministic straddle acquisition function (Eq. 4 in the main text) is choosing a single value for β 1/2 . If β 1/2 is too small (too aggressive), the search concentrates too narrowly around the threshold and may miss other informative regions; if it is too large (too conservative), it wastes samples in regions that are not critical for boundary identification. The conventional value β 1/2 = 1.96 turns out to be far too conservative in practice, which is what motivates randomizing β rather than fixing it. The randomized straddle acquisition function (Eq. 5 √ in the main text) instead draws β from a χ22 distribution at every 1/2 iteration. The expectation of the multiplier β is 2π/2 ≈ 1.25, noticeably more aggressive than the traditional fixed choice β 1/2 = 1.96. Smaller β 1/2 values concentrate the search near the threshold for focused boundary refinement, while occasional larger values drawn from the long exponential tail drive broader exploration and prevent over-concentration on a single already-known boundary fragment. This scheduling thus spends most iterations aggressively near the boundary while occasionally allowing broader exploration via the distribution’s long tail. Because χ22 has a simple closed-form CDF, β can be sampled efficiently via inverse transform sampling at negligible cost per iteration. Crucially, the schedule does not depend on the current sample count, so its adaptive character comes entirely from the distribution itself rather than from a hand-tuned annealing schedule. Inatsu et al. provided the full theoretical analysis and empirical comparison with fixed-β variants in Ref. 27. S1.3
Aggregating Per-Objective Scores in JAREX
Once the per-objective randomized-straddle scores STR∗i (x̂) ≥ 0 are available, JAREX must reduce them to a single acquisition value at each x̂. The main text uses a softmin-weighted sum (Eq. 10 in the main text). Several simpler alternatives are natural candidates: ∗ • Maximum: Amax jarex (x̂) = maxi STRi (x̂). OR-style; ignores objectives other than the locally most informative one and is essentially what the greedy baseline does. P ∗ 1 • Arithmetic mean: Aavg jarex (x̂) = m i STRi (x̂). Treats all objectives symmetrically and dilutes the joint signal with contributions from objectives that are confidently far from their thresholds. P ∗ • Fixed weighted sum: Aw jarex (x̂) = i wi STRi (x̂) with user-chosen wi . Requires prior knowledge about which objective deserves the most effort, which is precisely what characterization is trying to determine.
13
JAREX for Multi-Objective Process Characterization
A P REPRINT
The softmin combination interpolates between these: as τ → ∞ it reduces to the arithmetic mean, and as τ → 0 it concentrates on mini STR∗i (x̂). This AND-aggregation aligns with Pjoint being an intersection of per-objective pass regions (Eq. 8 in the main text): JAREX is large only when several objectives co-approach their thresholds. Flipping the sign of the exponent in Eq. 10 of the main text gives the corresponding softmax-weighted sum, P e+STR∗i /τ ∗ ∗ /τ STRi (x̂), which is the smooth counterpart of the maximum aggregator above. For STRi on the same i P +STR∗ j j e
scale as τ , the exponential weighting collapses almost all mass onto whichever objective is locally largest, so the aggregator becomes effectively a hard maximum: small noise-driven fluctuations in STR∗i flip which objective dominates from one iteration to the next, leading to numerically unstable, abruptly-switching acquisition values. Empirically this manifested as substantially worse joint-Jaccard convergence than softmin on the kinetic benchmark. The default τ = 0.5 is robust to moderate perturbations. S1.4
Hard vs. Soft Mask and BoTorch Integration
The joint optimistic pass region P̂joint (Eq. 8 in the main text) defines the set of candidates eligible for selection by JAREX. A naive implementation uses P̂joint as a hard indicator mask, multiplying the acquisition value by 1 inside P̂joint and 0 outside. This has two consequences that interact poorly with BoTorch’s default candidate optimization. First, when proposing new candidates, BoTorch [19] launches multi-start gradient-based searches from the best of a large set of random samples to refine candidates to a numerical local optimum of the acquisition function. A hard indicator mask is non-differentiable at the boundary of P̂joint and identically zero outside it, so the gradient is undefined on the boundary and vanishes in the exterior. Any starting point that falls outside P̂joint therefore remains stuck at zero acquisition with no gradient signal to guide it back into the feasible region. Second, P̂joint is generally irregular and non-convex; it cannot be expressed as a simple constraint such as a box or a polytope that BoTorch’s constrained optimizers can handle natively. In combination, these two properties reduce the candidate search to random sampling, discarding the numerical refinement that is one of BoTorch’s principal advantages. The soft mask in Eq. 13 of the main text avoids both problems. It preserves P̂joint exactly as the region where the mask equals 1, but outside this region the mask decays smoothly via sech(k d(x)), providing a non-zero gradient that pulls the optimizer back toward the feasible region from any starting point. The smoothness parameter k controls how steeply the mask decays; we use k = 1 by default and have observed robustness to reasonable variations. The trade-off is a thin band immediately outside P̂joint where the mask is non-zero but less than one, introducing a small amount of acquisition value in nominally infeasible territory. In practice this band is narrow enough that selected candidates remain essentially within P̂joint , while the differentiability restores BoTorch’s numerical candidate refinement. S1.5
A Note on Mathematical Rigor
The single-objective randomized straddle algorithm of Inatsu et al. [27] builds on the level-set classification framework [20, 21] and comes with a clean theoretical guarantee. Under standard Gaussian-process assumptions, its expected misclassification loss converges to zero at a sublinear rate. This section explains exactly where that argument breaks down for JAREX and what could in principle be recovered under a restrictive additional assumption. The argument is included because a reader interested in the theoretical status of the method deserves to see the concrete step that fails, not just a statement that “the proof does not apply”. Throughout this section we take the joint misclassification loss to be the sum of per-objective losses, ℓjoint (x) = t P m i=1 ℓt,i (x); equivalently, joint misclassification at x implies misclassification on at least one objective, so the joint loss is set-theoretically dominated by the sum of per-objective losses. The Inatsu proof combines four ingredients. 1/2
(i) A Gaussian-process concentration bound that places f (x) inside [µt−1 (x) ± βδ σt−1 (x)] with high probability, where βδ = 2 log(1/δ). (ii) A pointwise inequality ℓt (x) ≤ at−1,δ (x) between the misclassification loss and the straddle score; (iii) The distributional identity that, when δ ∼ U (0, 1), the random variable 2 log(1/δ) is exactly χ22 -distributed, which lets one swap the expectation over δ for an expectation over βt ∼ χ22 and, combined with the selection rule xt = arg maxx at−1 (x), upgrades the pointwise bound in (ii) to an expected bound evaluated at the selected point xt ; 14
JAREX for Multi-Objective Process Characterization
(iv) A Cauchy–Schwarz step combined with the maximum-information-gain bound converts the per-iteration bound into a sublinear cumulative rate.
A P REPRINT
2 t σt−1 (xt ) ≤ Cγt , which
P
The critical step for JAREX is (iii). In the single-objective case the selection rule maximizes the same scalar quantity at−1 that upper-bounds the loss, so the inequality can be evaluated at the selected point and then summed over iterations. In JAREX the selection rule maximizes M (x) · Ajarex (x), a data-dependent softmin-weighted combination of m per-objective straddle scores multiplied by a smoothly masked indicator of P̂joint . Neither the softmin-weighted combination nor the masked acquisition equals the single-objective straddle score of any individual objective, and in general the selected point x∗ is not the argmax of ai,t−1 for any fixed i. As a result, the per-objective inequality in step (ii) cannot be pushed to the selected point through step (iii), and the Cauchy–Schwarz / information-gain argument in step (iv) can no longer be applied per objective at x∗ to yield a sublinear rate. Per-objective ingredients (i) and (ii) themselves remain valid, and the joint misclassification loss is still bounded pointwise by the sum of per-objective losses by the set-inclusion argument above; what is missing is the link between these pointwise bounds and the specific point JAREX chooses to sample. A conditional, much weaker statement can be recovered only by restoring that link by force. If we additionally assumed that every objective is queried sufficiently often, for example by appending a round-robin fallback that guarantees each objective receives a non-vanishing share of the iteration budget, then each per-objective surrogate would satisfy Inatsu’s bound on its own subsequence of iterations. The joint misclassification loss, bounded above by the sum of per-objective losses, would then inherit a sublinear rate via a union bound, with the constants degraded by at most a factor of m and by the fraction of iterations assigned to each objective. This is a genuine but essentially trivial guarantee: it only says that if we force uniform coverage over objectives, no objective can be arbitrarily neglected and each individual boundary is eventually characterized. It does not capture any of the joint-boundary behavior that motivates JAREX in the first place, and it contradicts the intended behavior of the algorithm, which is to deprioritize regions confidently outside P̂joint and objectives whose boundary is locally uninformative. Imposing uniform perobjective coverage would recover a formal bound only by disabling the feature of JAREX that actually delivers the empirical gains in Section 3 of the main text, so we do not pursue this route. A tighter analysis that accommodates the softmin weighting and the soft mask, and that targets the joint misclassification loss directly rather than reducing it to per-objective losses, is what a full theoretical treatment of JAREX would require. We view this as an interesting but substantial direction for future work: the pointwise loss control step (ii) remains a promising starting point, but a new argument is needed to handle a selection rule that depends jointly on all m surrogates and on a data-dependent feasibility mask.
S2
Additional Benchmark Details
This section collects technical details on implementation choices referenced in the main text that apply to both singleobjective and multi-objective characterization benchmarks. The Jaccard index is a set-based similarity metric that scores only whether each parameter set is assigned the correct pass/fail label, making it particularly appropriate for process characterization, where the primary goal is to correctly identify the pass/fail boundary rather than to minimize prediction error uniformly across the parameter space. Traditional metrics such as mean absolute error (MAE) or root mean squared error (RMSE) weight all predictions equally and are therefore dominated by contributions from regions far from the boundary, where accuracy is less critical; a boundary-focused method may report a higher overall MAE while being substantially more accurate at the pass/fail threshold, which is the only region that affects characterization outcomes. A high Jaccard score also implicitly rewards global coherence of the response surface, since it is unlikely that a model produces accurate boundary predictions while being grossly wrong elsewhere. For all test cases, Gaussian process surrogate models were implemented using BoTorch [19] with the default Matérn kernel. Model hyperparameters were optimized using scipy.optimize with BFGS-family methods (L-BFGS-B for box-constrained optimization) with 10 multi-start restarts to ensure robust convergence. Initial hyperparameter values for each restart were sampled randomly within physically reasonable ranges to avoid local minima. Jaccard indices were calculated by discretizing the parameter space using quasi-random Sobol sequences, with the number of evaluation points determined adaptively based on the overall uncertainty of the surrogate model to ensure accurate boundary representation; the adaptive strategy increases the point density when the GP posterior variance is high near decision boundaries. 15
JAREX for Multi-Objective Process Characterization
S3
A P REPRINT
Single-Objective Characterization Details
This section provides detailed descriptions of the benchmark test cases used to evaluate the performance of the randomized straddle algorithm for single-objective characterization. We describe the analytical test functions used for benchmarking, starting with the one-dimensional double-well potential illustrated in Figure 1 of the main text. S3.1 S3.1.1
Single-Objective Analytical Functions Double-Well Potential (1D)
By design (Eq. 5 in the main text), the randomized straddle algorithm selects candidates that are either near the threshold or have high uncertainty. This preference can easily be visualized using a classical one-dimensional doublewell potential, which has two local minima separated by a barrier. The double-well potential is defined as, f (x) = ax4 − bx2 + c.
(S2)
2
b for the test function. −f (x) is shown as the blue curve in Figure 1 of the main text. We chose a = 1, b = 21 , c = 4a We flipped the sign purely for visualization purposes. The threshold h = − 2c is shown as the horizontal grey dashed line. 2 initial observations are added at the beginning of the campaign.
To validate performance across higher-dimensional spaces, we use three additional standard analytical test functions: the six-hump camel function (2D), Hartmann 4D function (4D), and Hartmann 6D function (6D). These functions are commonly used benchmarks in the optimization literature and provide increasing levels of complexity for testing boundary identification algorithms. Each function is described below. S3.1.2
Six-Hump Camel Function
The six-hump camel function is a two-dimensional multimodal test function. The standard form is defined as, x41 2 x21 + x1 x2 + (−4 + 4x22 )x22 , f (x1 , x2 ) = 4 − 2.1x1 + 3
(S3)
where x1 ∈ [−2, 2] and x2 ∈ [−1, 1]. For our characterization benchmarking, we apply a transformation to convert this minimization problem to a maximization problem, y(x1 , x2 ) = 11.0316 − f (x1 , x2 ). We use a threshold value of h = 10.6316 and classify the parameter space into regions where y(x1 , x2 ) ≥ h (pass) and y(x1 , x2 ) < h (fail). Initial sampling consists of 5 observations generated using Latin hypercube sampling (LHS). [30] S3.1.3
Hartmann 4D Function
The Hartmann 4D function is a four-dimensional multimodal test function defined as, 4 4 X X f (x) = − αi exp − Aij (xj − Pij )2 , i=1
(S4)
j=1
where x = (x1 , x2 , x3 , x4 ) ∈ [0, 1]4 , and the parameters are given by, α = (1.0, 1.2, 3.0, 3.2)T 10 3 17 3.5 17 0.1 0.05 10 A= 3 3.5 1.7 10 17 8 0.05 10 1312 1696 5569 −4 2329 4135 8307 P = 10 2348 1451 3522 4047 8828 8732
124 3736 2883 5743
The function has four local minima, with a global maximum at approximately f (x∗ ) ≈ 3.73. For characterization, we normalize the function output to y = (1.1 − f (x))/0.839 and use a threshold of h = 1.5. Initial sampling consists of 10 observations generated using LHS. [31] 16
JAREX for Multi-Objective Process Characterization
S3.1.4
A P REPRINT
Hartmann 6D Function
The Hartmann 6D function is a six-dimensional multimodal test function defined analogously to the 4D version,
f (x) = −
4 X
αi exp −
i=1
6 X
Aij (xj − Pij )2 ,
(S5)
j=1
where x = (x1 , x2 , x3 , x4 , x5 , x6 ) ∈ [0, 1]6 , with parameters,
α = (1.0, 1.2, 3.0, 3.2)T 10 3 17 3.5 1.7 8 17 0.1 8 14 0.05 10 A= 3 3.5 1.7 10 17 8 17 8 0.05 10 0.1 14 1312 1696 5569 124 8283 5886 2329 4135 8307 3736 1004 9991 P = 10−4 2348 1451 3522 2883 3047 6650 4047 8828 8732 5743 1091 381
The function has multiple local minima. For characterization, we normalize the function output to y = −(2.58 + f (x))/1.94 and use a threshold of h = 1.5. Initial sampling consists of 15 observations generated using LHS. [32]
S3.2
Single-Objective Benchmarks with Noise
To evaluate the robustness of the randomized straddle algorithm under realistic experimental conditions, we also benchmarked all methods with 5% Gaussian noise added to the analytical function outputs. The noise is modeled as N (0, σ 2 ) where σ = 0.05 × |y|, representing a common level of experimental noise in real-world applications. Fig. S1 shows the benchmark results with noise for all analytical functions. 17
JAREX for Multi-Objective Process Characterization
(a) Six-Hump Camel (2D)
2-level (7)
1.00
Straddle
0.75
CCD (11)
UCB (β = 6) Feasibility EI
Space Filling UCB (β = 1.96)
A P REPRINT
0.44
0.50 0.25
0.07
0.00 29
41
53
65
(b) Hartmann 4D
2-level (19)
0.75 0.50 0.25
0.00 0.00
0.00 10
40
70
100
130
160
(c) Hartmann 6D
2-level (67)
1.00 0.75
CCD (79)
Jaccard Index
17
CCD (27)
5 1.00
0.50 0.18
0.25 0.01
0.00 15
75
135
195
255
Number of Observations
315
DOE
Fig. S1: Benchmark results for single-objective characterization with 5% Gaussian noise tested on various analytical functions: six-hump camel (2D), Hartmann 4D (4D), and Hartmann 6D (6D). A Gaussian noise N (0, σ 2 ) with σ = 0.05 × |y| is added to all function outputs. Each method is run with 10 independent campaigns with different random seeds, and the mean and standard deviation of the Jaccard index are plotted. All methods show some performance degradation compared to the noise-free case (Fig. 2 in main text), but the relative performance ranking remains consistent. The randomized straddle algorithm continues to demonstrate superior efficiency and stability, particularly for high-dimensional problems. While convergence to near-perfect characterization (approaching 100% Jaccard index) is achievable with sufficiently large sample budgets, economical sample sizes for practical method development are used here. The side panel (DOE) in each transformation reports the RSM designs, a full two-level factorial and a central composite design (CCD), as single-design Jaccard values, with run counts in parentheses.
S4
Multi-Objective Characterization Details
S4.1
Kinetic Model
The kinetic model consists of a main catalytic pathway with a competing reversible side reaction: kmain (T, Cat)
Y + Z −−−−−−−→ P kside,fwd (Cat)
− 2Y + Z − ↽ −− −− −− −− −− −⇀ −X kside,rev (T )
(S6) (S7)
where Z is the reactant, Y is the co-reactant, P is the desired product, X is an undesired impurity, and Cat represents the catalyst. The main reaction exhibits both temperature and catalyst dependence (high activation energy), while the side forward reaction depends only on catalyst (no temperature dependence) and the side reverse reaction depends only on temperature (no catalyst dependence). 18
JAREX for Multi-Objective Process Characterization
S4.1.1
A P REPRINT
Rate Expressions
The reaction kinetics is described by the following system of ordinary differential equations (ODEs), d[Z] = −rmain − rsf + rsr (S8) dt d[Y] = −rmain − 2 · rsf + 2 · rsr (S9) dt d[P] = rmain (S10) dt d[X] = rsf − rsr (S11) dt d[Cat] = −kdecay · [Cat]. (S12) dt The stoichiometry reflects that the main reaction consumes both Y and Z, while the side reaction involves 2 moles of Y per mole of X. The main reaction rate (Y + Z → P) follows a Langmuir-Hinshelwood dual-site mechanism with substrate inhibition: rmain = kmain (T, [Cat]) ·
[Y] · [Z] , (1 + [Y]/KY,main + [Z]/KZ,main )2
where the rate constant exhibits Arrhenius temperature dependence and catalyst substrate inhibition: Ea,main [Cat] kmain (T, [Cat]) = Amain exp − · 2 , RT Km,main + [Cat] + K[Cat] i,main
(S13)
(S14)
with Amain = 6 × 106 L/(mol·min) as the pre-exponential factor, Ea,main = 40 kJ/mol as the activation energy (high, providing temperature-dependent purity control), Km,main = 0.01 mol/L as the Michaelis constant, Ki,main = 0.15 mol/L as the substrate inhibition constant, KY,main = 0.5 mol/L and KZ,main = 1.2 mol/L as the competitive adsorption constants for Y and Z, respectively. The squared denominator reflects a dual-site mechanism where both Y and Z compete for catalyst sites. R = 8.314 J/(mol·K) is the gas constant. The forward side reaction rate (2Y + Z → X) follows a Langmuir-Hinshelwood mechanism with fractional order in Y and catalyst substrate inhibition (no temperature dependence): rsf = ksf ([Cat]) ·
[Y]nY,sf · [Z] , (1 + [Y]/KY,sf + [Z]/KZ,sf )2
(S15)
where the rate constant depends only on catalyst loading (no temperature dependence): ksf ([Cat]) = ksfmax ·
[Cat] 2
Km,sf + [Cat] + [Cat] Ki,sf
,
(S16)
with ksfmax = 4.5 L2 /(mol2 ·min) as the maximum rate at saturating catalyst, Km,sf = 0.015 mol/L as the Michaelis constant, Ki,sf = 0.35 mol/L as the substrate inhibition constant, KY,sf = 4.0 mol/L and KZ,sf = 1.5 mol/L as the competitive adsorption constants, and nY,sf = 2.0 as the fractional order in Y. The squared denominator reflects dual-site competition. The absence of temperature dependence makes this pathway relatively more favorable at low temperatures. The reverse side reaction rate (X → 2Y + Z) depends only on temperature (no catalyst dependence) with product inhibition: [X] rsr = ksr (T ) · , (S17) 1 + [Y]/KY,rev + [Z]/KZ,rev where the rate constant follows Arrhenius temperature dependence: Ea,sr ksr (T ) = Asr · exp − , (S18) RT with Asr = 5×106 min−1 as the pre-exponential factor, Ea,sr = 32 kJ/mol as the activation energy, KY,rev = 1.0 mol/L and KZ,rev = 0.8 mol/L as product inhibition constants. The relatively low activation energy (compared to the main 19
JAREX for Multi-Objective Process Characterization
A P REPRINT
reaction) and absence of catalyst dependence makes this pathway relatively more favorable at higher temperatures, enabling impurity recycling. The following table summarizes all kinetic parameters used in the model. Parameter
S4.1.2
Value
Description
Main reaction: Y + Z → P Amain 6 × 106 L/(mol·min) Ea,main 40 kJ/mol Km,main 0.01 mol/L Ki,main 0.15 mol/L KY,main 0.5 mol/L KZ,main 1.2 mol/L
Pre-exponential factor Activation energy Michaelis constant Substrate inhibition constant Y competitive adsorption constant Z competitive adsorption constant
Side forward reaction: 2Y + Z → X ksfmax 4.5 L2 /(mol2 ·min) Km,sf 0.015 mol/L Ki,sf 0.35 mol/L KY,sf 4.0 mol/L KZ,sf 1.5 mol/L nY,sf 2.0
Maximum rate at saturating catalyst Michaelis constant Substrate inhibition constant Y competitive adsorption constant Z competitive adsorption constant Fractional order in Y
Side reverse reaction: X → 2Y + Z Asr 5 × 106 min−1 Ea,sr 32 kJ/mol KY,rev 1.0 mol/L KZ,rev 0.8 mol/L
Pre-exponential factor Activation energy Y product inhibition constant Z product inhibition constant
Catalyst decay kdecay 2 × 10−4 min−1
First-order decay constant
Objectives
The three objectives are calculated from the final concentrations as follows: [Z]0 − [Z]final [maximize] [Z]0 [P]final Product P purity = [maximize] [P]final + [X]final uZ [Z]0 + uY [Y]0 + uCat [Cat]0 + usol Unit product cost = [P]final Z conversion rate =
(S19) (S20) [minimize]
(S21)
The unit cost represents the total cost per mole of product P formed ($/mol), accounting for: • Material costs: reactant Z (uZ = 10 $/mol), co-reactant Y (uY = 45 $/mol), and catalyst (uCat = 6, 500 $/mol) • Solvent cost: fixed cost per batch (usol = 1 $/batch) All costs are normalized by the amount of product P formed. Note that temperature T does not directly enter the cost function but affects the kinetics and thus the final product yield. For multi-objective characterization benchmarking, we define the following target specifications: • Purity: Product purity ≥ 99.3% • Conversion: Reactant Z conversion ≥ 90% • Cost: Unit product cost ≤ 570 $/mol P 20
JAREX for Multi-Objective Process Characterization
A P REPRINT
These thresholds define a pass/fail classification of the four-dimensional parameter space. The goal of characterization is to identify the boundary between the pass region (where all three constraints are satisfied) and the fail region (where at least one constraint is violated). S4.1.3
Integration Details
The system of ordinary differential equations is integrated numerically using the BDF (Backward Differentiation Formula) method from scipy.integrate.solve_ivp. The integration time span is t ∈ [0, tfinal ] minutes, where tfinal is typically 120 min, with early termination if 99.9% of substrate Z is consumed. Initial conditions are specified by the decision variables: • Reactant Z initial concentration [Z]0 ∈ [0.8, 3.0] mol/L • Co-reactant Y initial concentration [Y]0 ∈ [0.8, 6.0] mol/L • Reaction temperature T ∈ [273, 333] K • Catalyst loading [Cat]0 ∈ [0.005, 0.2] mol/L with fixed initial values [P]0 = 0 mol/L and [X]0 = 0 mol/L. The solver uses a relative tolerance of 10−6 and an absolute tolerance of [Z]0 × 10−7 to ensure accurate integration. The output consists of species concentrations at the final time point. S4.2
Transformation Strategies
The kinetic model objectives (purity, conversion, cost) exhibit highly skewed distributions, making standardization alone less effective for Gaussian process modeling. This section provides detailed analysis of transformation strategies and their impact on method performance. Two transformation strategies were evaluated: • Standard transformation: Simple standardization (zero mean, unit variance) applied to each objective independently. • Log transformation: Logit transformation for purity and conversion (bounded objectives) combined with log transformation for cost (positive objective), followed by standardization. The log-based transformations help stabilize the surrogate model by mapping the highly skewed objective distributions to more Gaussian-like distributions, which better match the GP modeling assumptions. Full benchmark results with both transformation strategies are shown in Fig. S2 for serial characterization and Fig. S3 for batched characterization. The log transformation significantly improves the performance of the RSM designs, SF, and greedy randomized straddle methods. However, JAREX demonstrates remarkable numerical stability, achieving strong performance with standard transformation and only modest improvement with log transformation. This robustness is a desirable property for real-world applications where the optimal transformation strategy may not be known beforehand. S4.3
Quantitative Analysis of Multi-Objective Trade-offs
As discussed in the main text, the joint Jaccard index can exceed the worst individual Jaccard index, which appears counter-intuitive. This section provides a quantitative breakdown of the underlying mechanisms. The Jaccard index is defined as J = TP/(TP + FP + FN). Because the denominator counts errors over the specific region being evaluated, the joint Jaccard index is not the product of the individual ones, Y J(Â1 ∩ Â2 ∩ Â3 , A1 ∩ A2 ∩ A3 ) ̸= J(Âi , Ai ), (S22) i
even when the objectives are fully independent. In the kinetic model benchmark, the product Jpurity × Jconv × Jcost = 0.846 × 0.859 × 0.621 = 0.451, which is far below the observed Jjoint = 0.718. Fig. S4 shows the Venn diagram of the true pass sets evaluated on a grid of N = 10,000 points. The joint pass region covers 16.4% of the parameter space, roughly 2.4 times larger than the 7.8% expected if the three objectives were independent (0.407 × 0.438 × 0.436 = 0.078). This strong positive correlation arises because all three objectives are physically coupled through the product yield Pf . High conversion and high purity both drive Pf up, which in 21
JAREX for Multi-Objective Process Characterization
0.8
CCD (27)
(a) Transform: Standard
2-level (19)
1.0
A P REPRINT
0.6 0.35
0.2
1.0
(b) Transform: Log
0.8 0.6 SF Greedy: Cv-Pu-UC Greedy: Pu-Cv-UC Greedy: UC-Cv-Pu JAREX
0.4 0.2 0.0 8
32
56
80
104
Number of Observations
128
CCD (27)
0.00
0.0
2-level (19)
Jaccard Index
0.4
0.47
0.05
DOE
Fig. S2: Full benchmark results for sequential multi-objective characterization showing both standard (panel a) and log (panel b) transformations. JAREX (black) maintains robust performance across both transformation strategies, while other methods show stronger dependence on transformation choice. The side panel (DOE) in each transformation reports the RSM designs, a full two-level factorial and a central composite design (CCD), as single-design Jaccard values, with run counts in parentheses.
turn drives cost per unit product down. As a result, conditions satisfying purity and conversion have a chance to automatically satisfy the cost criterion,
P (cost passes | purity ∩ conv pass) = 22
1843 = 70.5%. 2614
(S23)
JAREX for Multi-Objective Process Characterization
A P REPRINT
1.0
Jaccard Index
0.8
0.6
0.4
Batch size=1, Transform: Log Batch size=4, Transform: Log Batch size=8, Transform: Log Batch size=1, Transform: Standard Batch size=4, Transform: Standard Batch size=8, Transform: Standard
0.2
0.0 0
24
48
72
96
120
Number of Iterations
Fig. S3: Full batched characterization results comparing serial (q = 1), moderate (q = 4), and high (q = 8) batch sizes with both standard (panel a) and log (panel b) transformations. Batching provides consistent benefits across transformation strategies.
Product purity ≥ 99.3% Cost per unit product ≤ 570 Reactant conversion ≥ 90% Product purity
Cost per unit product
40.79%
24.31%
40.67%
16.43% 24.13%
23.28%
40.83% Reactant conversion
Fig. S4: Venn diagram of the true pass sets for the three objectives in the kinetic model, evaluated on a grid of 10,000 points. Each circle represents the fraction of parameter space satisfying the corresponding threshold. The triple intersection (joint pass region, 16.4%) is 2.4× larger than expected under independence (7.8%), reflecting the strong positive correlation among objectives driven by the shared dependence on product yield Pf . The cost pass set contains a large “extra” region, points that pass the cost threshold but fail purity or conversion, comprising 25.1% of the parameter space (|C \ (A ∩ B ∩ C)| = 2513 points). These points lie near the cost boundary and are intrinsically difficult for the cost surrogate to classify correctly, inflating both its false positive (FPcost = 1352) and false negative (FNcost = 812) counts. However, because these points fail at least one other objective, the joint classifier correctly rejects them on purity or conversion grounds; they become true negatives for the joint evaluation rather than contributing to joint FP or FN. This error attribution effect is visible in the confusion matrix comparison obtained from the GP surrogate models of one sequential characterization campaign. 23
JAREX for Multi-Objective Process Characterization
Recall Precision FP count
A P REPRINT
Cost alone
Joint
3544/4356 = 81.4% 3544/4896 = 72.4% 1352
1494/1843 = 81.1% 1494/1731 = 86.3% 237
0.621
0.718
Jaccard
The recall rates are nearly identical (∼81%), meaning both classifiers find roughly the same fraction of their respective true pass regions. The decisive difference is precision. The joint classifier achieves 86% compared to 72% for cost alone, because the intersection operation eliminates most false positives from the cost extra region. This precision improvement is the direct mechanism by which Jjoint exceeds Jcost . In summary, the observation that the joint Jaccard index exceeds the worst individual Jaccard index is not a paradox but a natural consequence of (i) objective correlation concentrating the joint pass set in a well-characterized core, and (ii) the intersection operation filtering out classification errors that occur in the periphery of individual pass sets. S4.4
Two-Dimensional Interaction Plots
T=300
0.15 0.1 0.05
T=305
0.15 0.1 0.05
T=310
[Cat0]
Fig. S5 shows the same two-dimensional interaction effects as Fig. 6 in the main text, but with the JAREX model trained on 80 samples (8 initial and 72 iterative) instead of 40. Doubling the sample budget leads to noticeably sharper boundaries and expanded high-confidence regions. In particular, the 95% confidence pass regions (dark green) grow substantially, indicating that the model is considerably more certain about its predictions. Meanwhile, the narrow yellow band (mean pass only) that borders the fail region in the 40-sample model shrinks or disappears in many panels, confirming that the additional data has resolved much of the remaining boundary uncertainty. The overall topology of the pass/fail regions remains consistent between the two models, demonstrating that even the 40-sample model captures the correct qualitative structure of the parameter space. However, the wider high-confidence regions in the 80-sample model would provide stronger statistical support for defining proven acceptable ranges in a regulatory context.
0.15 0.1 0.05 1.5 2 2.5 [Y0]=2.0
Fail
Mean pass
1.5 2 2.5 [Y0]=3.0
1.5 2 2.5 [Y0]=2.0
70% confidence
1.5 2 2.5 [Y0]=3.0
1.5 2 2.5 [Y0]=2.0
1.5 2 2.5 [Y0]=3.0
95% confidence
1.5 2 2.5 [Y0]=2.0
1.5 2 2.5 [Y0]=3.0
[Z0] Product Purity
Reactant Conversion
Cost per Unit Product
Joint Evaluation
Fig. S5: Two-dimensional interaction effects between the initial reactant concentration [Z0 ] and initial catalyst concentration [Cat0 ] obtained from JAREX with other parameters fixed at various values. 80 samples are used in the campaign to train the models. Each panel shows the predicted pass region at three confidence levels derived from the Gaussian process posterior: mean pass (predicted mean exceeds threshold) in yellow, 70% confidence in light green, and 95% confidence in dark green, with the remaining area (red) classified as fail based on the predicted mean values. The last column shows the joint evaluation, where confidence levels are determined by the intersection of individual objective pass regions at each confidence level. Compared with Fig. 6 in the main text (40 samples), the 95% confidence regions are substantially wider, indicating improved model certainty with the larger training set.
References [1] ICH Q8(R2): Pharmaceutical Development. Technical report, International Council for Harmonisation, 2009. 24
JAREX for Multi-Objective Process Characterization
A P REPRINT
[2] ICH Q11: Development and Manufacture of Drug Substances. Technical report, International Council for Harmonisation, 2012. [3] Lawrence X. Yu, Gordon Amidon, Mansoor A. Khan, Stephen W. Hoag, James Polli, G. K. Raju, and Janet Woodcock. Understanding Pharmaceutical Quality by Design. AAPS J., 16(4):771–783, 2014. doi: 10.1208/ s12248-014-9598-3. [4] Anurag S. Rathore and Helen Winkle. Quality by design for biopharmaceuticals. Nat. Biotechnol., 27(1):26–34, 2009. doi: 10.1038/nbt0109-26. [5] Robert A. Lionberger, Sau L. Lee, LaiMing Lee, Andre Raw, and Lawrence X. Yu. Quality by Design: Concepts for ANDAs. AAPS J., 10(2):268–276, 2008. doi: 10.1208/s12248-008-9026-7. [6] Lawrence X. Yu and Michael Kopcha. The future of pharmaceutical quality and the path to get there. Int. J. Pharm., 528(1–2):354–359, 2017. doi: 10.1016/j.ijpharm.2017.06.039. [7] Helena B. Grangeia, Catarina Silva, Sergio P. Simões, and Marco S. Reis. Quality by design in pharmaceutical manufacturing: A systematic review of current status, challenges and future perspectives. Eur. J. Pharm. Biopharm., 147:19–37, 2020. doi: 10.1016/j.ejpb.2019.12.007. [8] Dave am Ende, Kimber S. Bronk, Jason Mustakis, Greg O’Connor, Charles L. Santa Maria, Richard Nosal, and Timothy J. N. Watson. API Quality by Design Example from the Torcetrapib Manufacturing Process. J. Pharm. Innov., 2(3):71–86, 2007. doi: 10.1007/s12247-007-9015-x. [9] Pankaj D. Rege, Andreas Schuster, Jens Lamerz, Christian Moessner, Wolfgang Göhring, Pirmin Hidber, Helmut Stahr, Oana Mihaela Andrei, Janine Burren, Alexandre Moesching, Daniel Coleman, and Stefan Hildbrand. QbD Approach to Process Characterization and Quantitative Criticality Assessment of Process Parameters. Org. Process Res. Dev., 28(4):1003–1017, 2024. doi: 10.1021/acs.oprd.3c00356. [10] John J. Peterson and Kevin Lief. The ICH Q8 Definition of Design Space: A Comparison of the Overlapping Means and the Bayesian Predictive Approaches. Stat. Biopharm. Res., 2(2):249–259, 2010. doi: 10.1198/sbr. 2009.08065. [11] Douglas C. Montgomery. Design and Analysis of Experiments. Wiley, 9 edition, 2017. [12] Steven A. Weissman and Neal G. Anderson. Design of Experiments (DoE) and Process Optimization. A Review of Recent Publications. Org. Process Res. Dev., 19(11):1605–1633, 2015. doi: 10.1021/op500169m. [13] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando de Freitas. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proc. IEEE, 104(1):148–175, 2016. doi: 10.1109/JPROC.2015. 2494218. [14] Peter I. Frazier. A Tutorial on Bayesian Optimization. arXiv preprint arXiv:1807.02811, 2018. [15] Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient Global Optimization of Expensive BlackBox Functions. J. Glob. Optim., 13(4):455–492, 1998. doi: 10.1023/A:1008306431147. [16] Peter Auer, Nicolo Cesa-Bianchi, and Paul Fischer. Finite-time Analysis of the Multiarmed Bandit Problem. Mach. Learn., 47:235–256, 2002. doi: 10.1023/A:1013689704352. [17] Niranjan Srinivas, Andreas Krause, Sham M. Kakade, and Matthias Seeger. Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. In Proc. 27th Int. Conf. Mach. Learn., pages 1015–1022, 2010. [18] Harold J. Kushner. A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. J. Basic Eng., 86(1):97–106, 1964. doi: 10.1115/1.3653121. [19] Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, Andrew Gordon Wilson, and Eytan Bakshy. BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. In Adv. Neural Inf. Process. Syst., volume 33, pages 21524–21538, 2020. [20] Brent Bryan, Robert C. Nichol, Christopher R. Genovese, Jeff Schneider, Christopher J. Miller, and Larry Wasserman. Active Learning For Identifying Function Threshold Boundaries. In Y. Weiss, B. Schölkopf, and J. Platt, editors, Adv. Neural Inf. Process. Syst., volume 18. MIT Press, 2005. URL https://proceedings. neurips.cc/paper_files/paper/2005/file/8e930496927757aac0dbd2438cb3f4f6-Paper.pdf. [21] Alkis Gotovos, Nathalie Casati, Gregory Hitz, and Andreas Krause. Active Learning for Level Set Estimation. In Proc. 23rd Int. Joint Conf. Artif. Intell., pages 1344–1350, 2013. [22] Fani Boukouvala and Marianthi G. Ierapetritou. Feasibility analysis of black-box processes using an adaptive sampling Kriging-based method. Comput. Chem. Eng., 36:358–368, 2012. doi: 10.1016/j.compchemeng.2011. 06.005. 25
JAREX for Multi-Objective Process Characterization
A P REPRINT
[23] Nirupaplava Metta, Rohit Ramachandran, and Marianthi Ierapetritou. A novel adaptive sampling based methodology for feasible region identification of compute intensive models using artificial neural network. AIChE J., 67 (2):e17095, 2021. doi: 10.1002/aic.17095. [24] MSD. obsidian: Bayesian Optimization for Pharmaceutical Development, 2024. URL https://github.com/ MSDLLCPapers/obsidian. [25] Julien Bect, David Ginsbourger, Ling Li, Victor Picheny, and Emmanuel Vazquez. Sequential design of computer experiments for the estimation of a probability of failure. Stat. Comput., 22(3):773–793, 2012. doi: 10.1007/ s11222-011-9241-4. [26] Clément Chevalier, Julien Bect, David Ginsbourger, Emmanuel Vazquez, Victor Picheny, and Yann Richet. Fast parallel kriging-based stepwise uncertainty reduction with application to the identification of an excursion set. Technometrics, 56(4):455–465, 2014. doi: 10.1080/00401706.2013.860918. [27] Yu Inatsu, Shion Takeno, Kentaro Kutsukake, and Ichiro Takeuchi. Active Learning for Level Set Estimation Using Randomized Straddle Algorithms. arXiv preprint arXiv:2408.03144, 2024. [28] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, High-Performance Deep Learning Library. In Adv. Neural Inf. Process. Syst., volume 32, pages 8026–8037, 2019. [29] Paul Jaccard. Étude comparative de la distribution florale dans une portion des Alpes et du Jura. Bull. Soc. Vaud. Sci. Nat., 37:547–579, 1901. [30] Sonja Surjanovic and Derek Bingham. Six-Hump Camel Function. Virtual Library of Simulation Experiments: Test Functions and Datasets, 2013. URL https://www.sfu.ca/~ssurjano/camel6.html. Retrieved March 10, 2026. [31] Sonja Surjanovic and Derek Bingham. Hartmann 4-D Function. Virtual Library of Simulation Experiments: Test Functions and Datasets, 2013. URL https://www.sfu.ca/~ssurjano/hart4.html. Retrieved March 10, 2026. [32] Sonja Surjanovic and Derek Bingham. Hartmann 6-D Function. Virtual Library of Simulation Experiments: Test Functions and Datasets, 2013. URL https://www.sfu.ca/~ssurjano/hart6.html. Retrieved March 10, 2026.
26