A Metamodel-Based General-Purpose Autocalibration Tool for Simulation Models

Abstract

Background. Simulation calibration is the process of configuring a simulation model’s parameters to improve the agreement between the model output and the desired calibration targets (e.g., observed historical data). For most realistic simulation models, this calibration process can be quite computationally expensive, as it requires running the simulation model for each parameter combination. To alleviate this problem, metamodels offer a tradeoff between accuracy and computational efficiency for extensive simulative analysis with a highly complex parameter space. Method. In this study, we examine 4 simulation calibration approaches. Randomly-Simulate (RS) is a simulation-based benchmark widely used in the literature. Optimally-Predict (OP) is an optimization-based approach that we adapt to the simulation calibration setting. Building on these 2 baselines, we introduce 2 hybrid strategies, Predict-then-Simulate (PtS) and Simulate-then-Predict (StP), which combine simulation runs and metamodel-based optimization in complementary orders. We compare all 4 methods in terms of calibration accuracy and computational cost. Results. While the metamodel-based OP approach substantially reduced computational cost relative to RS and identified parameter combinations near the optimal configuration, the hybrid strategies delivered superior calibration performance. In particular, the PtS approach, which combines metamodel-based optimization with targeted simulation refinement, achieved on average a 46% reduction in total actual error compared with the RS benchmark, while maintaining computational efficiency. Conclusions. The study introduces a novel metamodel-based optimization approach to simulation calibration and illustrates its potential benefits for computationally expensive studies. While developed for deterministic targets, the method provides a foundation for future extensions to settings involving stochastic simulation outputs and other forms of model uncertainty. An open-access Python implementation of the proposed framework is provided to facilitate adoption and reproducibility.

Highlights

A metamodel-based optimization approach is proposed for calibrating simulation model parameters, which can offer computational advantages particularly in cases in which direct calibration is expensive due to complex or high-dimensional simulation models.

A hybrid approach that narrows down the search space via metamodel-based optimization and then uses simulation runs for fine-tuning offers a sweet spot in the tradeoff between computational efficiency and accuracy.

Keywords

simulation calibration metamodeling optimization over trained neural networks

Simulation models are widely used in health care, finance, and engineering due to their ability to replicate complex real-world systems with high fidelity.¹ Ensuring model validity requires calibrating input parameters so that simulated outputs align with observed data. The general calibration process is illustrated in Appendix A. Because simulations often involve many parameters and each run can take from several minutes to several days,² calibration can be computationally expensive. Traditional approaches such as grid search, random sampling, and heuristic optimization methods (e.g., simulated annealing or genetic algorithms) often require a large number of simulation runs and cannot guarantee identification of optimal parameter combinations.^3–5 Another widely used framework for calibration is Bayesian calibration, which represents input parameters as random variables with prior distributions and infers their posterior distribution given observed data.^6,7 While powerful and flexible, such approaches typically rely on sampling-based methods (e.g., Markov Chain Monte Carlo) to approximate the posterior and can be computationally demanding, particularly in high-dimensional settings or when simulation evaluations are costly.^8–11

To alleviate the computational cost of calibration, researchers have proposed using metamodels (surrogate models) to approximate the relationship between simulation inputs and outputs.² Once trained, a metamodel can be evaluated at a fraction of the cost of the original simulation, allowing broader exploration of the parameter space.¹² However, the number of candidate parameter combinations still grows exponentially with the number of inputs, and the “curse of dimensionality” remains a challenge even in metamodel-based calibration beyond a moderate number of parameters (p > 20). Common metamodeling approaches include polynomial regression, multivariate adaptive regression splines, Gaussian process models (kriging), radial basis function networks, and neural networks (NNs).¹³ Among these, NNs have gained popularity due to their flexibility and ability to approximate complex nonlinear relationships, with applications spanning healthcare, production, and energy systems.^14–18 In practice, metamodel-based calibration is typically implemented via sampling-based strategies, where candidate input configurations are evaluated using the metamodel and the best-performing solution is selected.¹⁹ While computationally efficient, such approaches do not guarantee identification of the optimal solution with respect to the metamodel-predicted error. A formal definition of the metamodel-based calibration problem is provided in Appendix B.

Our goal in this study is to propose an efficient general-purpose computational framework for the parameter calibration of simulation models. Central to this framework is a hybrid calibration strategy, Predict-then-Simulate (PtS), which first identifies the input vector that minimizes the calibration error predicted by a trained metamodel and subsequently conducts a targeted local search using simulation runs in the neighborhood of the identified solution. To support this strategy, we formulate the calibration task in deterministic simulations as a mixed-integer optimization (MIO) problem by embedding trained NNs into mathematical optimization models. This formulation directly addresses a core challenge in calibration, the combinatorial explosion of potential parameter combinations, and enables a structured and exact search over the feasible space with respect to the metamodel-predicted error. Furthermore, the framework allows the integration of logical constraints or domain-specific prior knowledge, making it well suited for complex calibration tasks.²⁰ We also discuss how the proposed framework can be extended to settings with stochastic simulation outputs.

This study contributes to the simulation calibration literature in several ways. First, we introduce hybrid calibration strategies that combine metamodel-based optimization with targeted simulation. Second, we adapt MIO formulations for trained NNs to the context of simulation calibration, enabling exact optimization over the metamodel with respect to predicted calibration error. Third, we demonstrate how contextual (nontunable) inputs can be incorporated into the metamodel to allow generalization across different calibration settings. Finally, we evaluate the proposed framework using both a real-world epidemic simulation model and a stylized numerical benchmark to assess performance across diverse scenarios. These contributions aim to provide a general-purpose calibration tool that is adaptable, computationally efficient, and broadly applicable across domains.

Methods

In this section, we describe how a trained NN metamodel can be embedded within an optimization framework to identify input parameter combinations that minimize metamodel-predicted calibration error and how this formulation supports a family of calibration strategies under a limited simulation budget.

Rather than relying on sampling-based approaches, we formulate metamodel-based calibration as an optimization problem in which the trained NN is embedded directly into a mixed-integer linear programming model. This enables a systematic search over the feasible parameter space and identifies the input configuration that minimizes the metamodel-predicted error. The formulation builds on prior work showing that feed-forward NNs with ReLU activation functions can be represented exactly within MIO models.²⁰ Full details of the mixed-integer formulation are provided in Appendix C. The solution to this problem provides a candidate parameter vector that is optimal with respect to the metamodel-predicted error; when the metamodel approximates the underlying simulation model well, this solution is expected to lie near the true optimal calibration.

In practice, calibration is often performed under a limited budget of simulation runs, which motivates the need to allocate these costly evaluations efficiently. A baseline strategy is to randomly sample parameter combinations, evaluate them using the simulation model, and select the best-performing solution. We refer to this strategy as Randomly-Simulate (RS). An alternative is to rely entirely on the metamodel-based optimization described in Appendix C, which we refer to as Optimally-Predict (OP). The OP strategy assumes the availability of a pretrained metamodel and therefore does not account for the simulation runs used during training. Building on these 2 baseline strategies, we propose hybrid calibration strategies that combine metamodel-based optimization with targeted simulation to improve calibration performance.

Hybrid Calibration Strategies

We now describe the hybrid calibration strategies that combine metamodel-based optimization with simulation. The first hybrid strategy is called Simulate-then-Predict (StP) (“Randomly” and “Optimally” prefixes are dropped for brevity). As the name suggests, we first apply the RS strategy (i.e., use our budget of r simulation runs in this stage), pick the point with the smallest actual error as an anchor, and then run the OP strategy within the proximity of that anchor point (by imposing lower and upper bounds on the input parameters). Finally, we choose the point that OP finds as the calibrated input parameter combination.

The second hybrid strategy is called Predict-then-Simulate (PtS), which combines the 2 base strategies in the reverse order. Namely, we first apply OP on the whole feasible input region, pick the point with the smallest predicted error as an anchor, then run RS within the proximity of that anchor point (i.e., use our budget of r simulation runs in this narrow region), and choose the point with the smallest actual error among the r simulated points.

Figure 1 graphically illustrates the 4 calibration strategies described above. Red dots represent the actual simulation runs, and the red quad is the anchor point with the smallest error among them. Blue-shaded regions represent the feasible space within which OP is run, and the blue quad represents the “optimal” point found by the corresponding OP run (i.e., the anchor point). The size of the proximity region around an anchor point can be defined by the user. In our computational results, we show the impact of choosing different sizes for proximity regions on the accuracy of the calibrated parameters.

Figure 1

Illustration of 4 calibration strategies. Red dots indicate simulation runs, and the red cross marks the best-performing point among them. Blue crosses and shaded regions depict the metamodel-based search (Optimally-Predict; OP), with hybrid strategies (Simulate-then-Predict [StP] and Predict-then-Simulate [PtS]) combining simulation and metamodel runs within user-defined proximity boxes.

Case Study: FluTE Simulation Model

To evaluate the proposed calibration strategies, we used FluTE, a publicly available stochastic influenza epidemic simulation model developed by Chao et al.²¹ The model simulates the spread of influenza in large metropolitan populations and has been widely used to study intervention strategies and epidemic dynamics. The full FluTE model includes 42 input parameters (e.g., transmissibility, vaccination coverage, and response delay) and produces multiple epidemiological outputs such as symptomatic attack rates across population groups. A single simulation run can take up to 2 h for a population of 10 million individuals, making large-scale calibration computationally expensive. For this study, we use a reduced version of the model to enable training of the metamodel. Specifically, we simulated a population of 500,000 individuals and restrict attention to 6 tunable input parameters and 9 output variables that exhibit substantial variability. With this setup, each simulation run takes approximately 5 to 6 min. The selected inputs and outputs, along with additional details of the FluTE model and simulation setup, are provided in Appendix D. We use an NN as the metamodel due to its compatibility with MIO formulations for calibration. Details of the network architecture, training procedure, and normalization are provided in Appendix E.

Proof-of-Concept Numerical Experiment: Rosenbrock Function

To complement the FluTE case study with a controlled benchmark setting, we introduce a stylized calibration experiment using a well-known benchmark function from the optimization literature: the Rosenbrock function.²² This function provides a controlled, analytically defined surface with a known global minimum, allowing us to rigorously test the behavior of our calibration algorithms independent of domain-specific assumptions. By treating the Rosenbrock function as a stand in for a simulation model, we can more directly evaluate the performance and consistency of different calibration strategies (RS, OP, PtS, and StP) in identifying input parameters that reproduce a given target output. The 2-dimensional Rosenbrock function is defined as

f (x, y; a, b) = (a - x)^{2} + b (y - x^{2})^{2}

where $x$ and $y$ are the inputs to be calibrated and a and b are fixed context parameters that define the shape and curvature of the function. Importantly, a and b are not part of the calibration problem; instead, they simulate variation in the underlying environment or scenario in which calibration takes place.

Results

In this section, we compare the performance of 4 calibration techniques introduced in the “Methods” section, namely, 1) RS, 2) OP, 3) StP, and 4) PtS. We use a budget of $r = 100$ simulation runs for each method (except for OP, which does not require any simulation runs). The main metric we use for comparing different calibration methods is the total actual error ( $Δ = \sum_{j = 1}^{P} | S_{j} (x) - {\hat{z}}_{j} |$ ), although, when needed, we also report the total predicted error ( $δ = \sum_{j = 1}^{P} | M_{j} (x) - {\hat{z}}_{j} |$ ). We first give the main results with the case study FluTE simulation model and then extend the study to different contexts using the Rosenbrock function.

Comparing Base Strategies

In the first analysis, we compare the performance of the 2 base calibration strategies, RS and OP. To do so, we conduct 50 independent experiments, each defined by a different randomly selected target output vector ${\hat{z}}^{i} = ({\hat{z}}_{1}^{i}, {\hat{z}}_{2}^{i}, \dots, {\hat{z}}_{P}^{i}$ ), representing distinct calibration scenarios. These target output vectors were generated by uniformly sampling input parameter vectors from the feasible region and running the simulation model to obtain the corresponding outputs. The 50 calibration targets were selected from a hold-out set not used for metamodel training, and because each target was generated from a known input vector, the true calibrated inputs are known for evaluation. For each of these experiments, we evaluate the performance of both RS and OP methods in identifying input parameters that produce output values close to the given target. Specifically, we let $x_{i}^{RS}$ and $x_{i}^{OP}$ denote the input parameter vectors calibrated using RS and OP, respectively. The actual errors for each method in experiment i are denoted by $Δ_{i}^{RS} = \sum_{j = 1}^{P} | S_{j} (x_{i}^{RS}) - {\hat{z}}_{j}^{i} |$ and $Δ_{i}^{OP} = \sum_{j = 1}^{P} | S_{j} (x_{i}^{OP}) - {\hat{z}}_{j}^{i} |$ , where $S_{j} (.)$ represents the simulation output for output dimension j. For comparison purposes, we also compute the corresponding metamodel-predicted (proxy) errors for both methods, denoted by $δ_{i}^{RS} = \sum_{j = 1}^{P} | M_{j} (x_{i}^{RS}) - {\hat{z}}_{j}^{i} |$ and $δ_{i}^{OP} = \sum_{j = 1}^{P} | M_{j} (x_{i}^{OP}) - {\hat{z}}_{j}^{i} |$ , where $δ_{i}^{RS}$ is evaluated post hoc at the input selected by RS.

For comparison purposes, we also compute the corresponding metamodel-predicted (proxy) errors for both methods, denoted by δ_i^RS and δ_i^OP, where δ_i^RS is evaluated post hoc at the input selected by RS.

Figure 2 summarizes the results from experiments with 50 random target outputs. It provides the boxplot of predicted errors from 50 experiments for RS (red) and OP (blue), respectively, in the first 2 boxplots. On the other hand, last 2 boxplots are provided for the actual errors. Figure 2 serves as a confirmation that the OP model, indeed, selects better (in fact, the best) input parameters with respect to the proxy/predicted error (i.e., $δ_{i}^{OP} \leq δ_{i}^{RS}$ holds for each individual observation i). Thus, we can conclude that the OP method works as intended, optimizing the proxy metric. Whether optimizing for a proxy metric translates into optimizing for an actual metric is the question we next turn to. To that end, actual error boxplots provide a stark contrast to predicted error boxplots. This time we see that the RS method (which selects the point with the smallest actual error among $r = 100$ randomly sampled input combinations) is clearly superior to OP. This contrast points out an important insight that calibration strategies like OP, which rely solely on a metamodel-based optimization (of a proxy error metric), should be handled with caution. Even when the metamodel approximates the simulation model with high fidelity, as in our case, where test sample R² values exceed 0.97 for all outputs (Table F.2), the selected input parameters may still lead to suboptimal calibration. The gap between predicted and actual errors reflects the metamodel’s inherent approximation error. In multioutput settings, this error accumulates across outputs, and the total discrepancy aligns closely with the sum of individual output Mean Absolute Error (MAEs). As shown in Table F.2, these small per-output errors can compound into a meaningful total error, despite the high R². Therefore, we next turn our attention to the hybrid strategies described in the “Methods” section.

Figure 2

Comparison of predicted error (PE) and actual error (AE) for Optimally-Predict (OP) and Randomly-Simulate (RS). Boxplots show the PE and AE distributions for 50 calibration targets. OP achieves lower predicted error than RS does, as expected, but RS outperforms OP in terms of actual error.

Hybrid (Fixed-Budget) Calibration Strategies

Comparison of the RS and OP base strategies suggests that if one would like more accurately calibrated input parameters, the better way of obtaining that is the RS strategy. It must be noted, however, that the RS strategy is significantly more costly compared with the OP method (100 × 5 min v. 30 s) since it requires actual simulation runs. Therefore, the following question becomes relevant: “If one has a budget of $r = 100$ simulation runs for a calibration task, how should the budget be best used?” Assuming we have a trained metamodel readily available and a fixed budget of $r = 100$ , we will next analyze the performance of 2 hybrid calibration strategies described earlier: 3) StP and 4) PtS.

Both hybrid strategies require an initial anchor point and a proximity range around that anchor point. StP uses its budget of $r$ simulation runs upfront with the RS strategy to produce a better anchor point and then searches for improved solutions within the vicinity of that point using the OP strategy. PtS, on the other hand, first finds an anchor point via the OP strategy (not necessarily the best one as the results in the previous subsection showed) and then uses its budget of $r$ simulation runs to find a better point within the vicinity of that anchor point via the RS strategy.

For both strategies, finding an anchor point is straightforward (it just requires running either RS or OP, respectively). Defining the proximity range of that anchor point for further search, however, requires user input. The wider the proximity range is, the less important the anchor point becomes. Figure 3 illustrates the search spaces for various proximity ranges in a 2-dimensional parameter space. The first image represents the optimal point found by the OP method. The last image represents the randomly sampled points in the RS phase. The point depicted with the anchor sign represents the point with the smallest actual error among $r$ samples. The other images show the vicinity of the anchor point widening from left to right by increasing the proximity range ( $R \in {\frac{1}{50}, \frac{1}{20}, \frac{1}{10}, \frac{1}{5}}$ ). For the next set of experiments, we use $R = \frac{1}{20}$ ; that is, we set the vicinity of the anchor point in each input dimension:

x_{i} \in AnchorPoin t_{i} \pm \frac{1}{20} \times (u b_{i} - l b_{i})

Figure 4 illustrates the actual error distributions from 50 experiments with randomly selected target outputs for all 4 strategies in boxplots. The median error values were calculated as 0.1408, 0.2059, 0.1803, and 0.0773. It becomes evident that the PtS strategy significantly outperforms the other 3 strategies. The performance of all methods can be ordered as follows: 4) PtS > 1) RS > 3) StP > 2) OP. These results highlight the added value of a smart use of limited computational budgets. Specifically, although solely relying on a metamodel-based calibration strategy, OP provides the worst performance; using OP to narrow down the search space for a more computationally intense search (i.e., random simulation runs), instead of blindly sampling the whole input space, provides a significant gain in terms of calibration accuracy without much added computational cost. Finally, the superior performance of PtS over StP provides an important insight that the order of execution of the base strategies does matter. Using first OP to narrow down the search space followed by an RS to pinpoint the best point performs significantly better compared with the other way around. On average, the PtS strategy achieved a 46% relative reduction in total actual error compared with the RS strategy, indicating substantial improvement in calibration accuracy.

Figure 3

Visualization of proximity ranges in anchor-based calibration strategies. Panels illustrate the effect of varying the proximity range R∈{1/50,1/20,1/10,1/5} around an anchor point (blue cross) found by Optimally-Predict (OP). Red dots indicate simulation evaluations, and red crosses mark the best-performing samples. As the proximity region expands, the influence of the anchor point diminishes. The final panel on the right shows the Randomly-Simulate (RS) strategy without any proximity restriction.

Figure 4

Actual error distributions for 4 calibration strategies. Boxplots show the total actual error across 50 calibration experiments for each method: Randomly-Simulate (RS), Optimally-Predict (OP), Simulate-then-Predict (StP), and Predict-then-Simulate (PtS). The PtS strategy achieves the lowest median error, outperforming all others and demonstrating the benefit of combining metamodel-based narrowing with targeted simulation.

Effect of Proximity Range on Calibration Accuracy

In the previous experiment, we fixed the proximity range around the anchor point to 1/10 of the original range of each input parameter. This was an arbitrary choice, so in the final set of experiments, we look at the effect of choosing different values for the proximity range. We showed before, in Figure 3, the rectangular borders drawn by taking the distance, in each parameter dimension, from the anchor point as much as the original parameter range multiplied by $R \in {\frac{1}{5}, \frac{1}{10}, \frac{1}{20}}$ . Thus, the search space was narrowed down by 5⁶, 10⁶, and 20⁶ times, respectively, since we have 6 input parameters. We will now look at the performance (i.e., the actual error distributions) of the PtS strategy for different values of proximity ranges $R \in {\frac{1}{1}, \frac{1}{5}, \frac{1}{10}, \frac{1}{20}, \frac{1}{50}, \frac{1}{1000}}$ . Note that $R = \frac{1}{1}$ corresponds to the full input domain and $R = \frac{1}{1000}$ corresponds (approximately) to the anchor point itself. Consequently, PtS with $R = \frac{1}{1}$ corresponds to RS, and PtS with $R = \frac{1}{1000}$ corresponds to the OP strategy, respectively.

Figure 5 shows the actual error distributions for the cases of $R =$ 1/1000 (OP), 1/50, 1/20, 1/10, 1/5, and 1/1 (RS) from left to right. The median error values were calculated as 0.2059, 0.0852, 0.0773, 0.0869, 0.0929, and 0.1408. As can be seen from the boxplots and the median error values, the performance of the PtS strategy improves as the proximity range around the anchor point expands until around the $R = \frac{1}{20}$ and then starts to deteriorate with further expansions of the proximity range.

Figure 5

Effect of proximity range on calibration accuracy in the predict-then-simulate (PtS) strategy. Boxplots depict actual error distributions for PtS under varying proximity ranges R∈{1/1000,1/50,1/20,1/10,1/5,1/1}, where R = 1/1000 approximates Optimally-Predict (OP) and R = 1/1 corresponds to Randomly-Simulate (RS). Performance improves as the search region expands from the anchor point up to R = 1/20, after which accuracy deteriorates.

These results confirm that the proximity range is a tunable hyperparameter that can significantly affect calibration performance. While we identified a “sweet spot” around R = 1/20 in our context, the optimal value is likely to vary depending on the problem structure, dimensionality, and metamodel accuracy. Practitioners are therefore encouraged to perform sensitivity analyses or validation-guided tuning when applying PtS or similar strategies in other settings.

Generalization across Contextual Parameter Settings

We performed calibration experiments with the Rosenbrock function following the same procedure used in our main analysis (section Effect of Proximity Range on Calibration Accuracy) starting by training the NN metamodel with 4500 (=9 × 500) randomly sampled instances. In our experiments with the Rosenbrock function, we vary $a \in {5, 10, 15}$ and $b \in {10, 20, 30}$ to generate a range of functional surfaces, reflecting different levels of nonlinearity and steepness. For each (a, b) combination, we generate 500 training data by sampling inputs $x ~ U (a - 1, a + 1)$ and $y ~ U (a^{2} - 5, a^{2} + 5)$ , and we compute the corresponding output $f (x, y; a, b)$ . The resulting dataset, after normalization, is used to train a feedforward NN metamodel with 3 hidden layers of 256 nodes each. The trained metamodel achieved a test sample R² value of 0.9005, indicating a reasonably accurate approximation of the underlying Rosenbrock surfaces across varying contexts. Then, for 50 randomly selected target outputs, each generated from a known (x, y ) pair under a fixed (a, b) context, we compared the performance of RS, OP, PtS, and StP in recovering the corresponding inputs. This stylized setting isolates calibration behavior from the complexities of domain-specific simulation models and validates that the observed patterns in our primary case study are not artifacts of model-specific dynamics. As shown in Figure 6, the comparison of the 4 strategies with the Rosenbrock function yields a similar pattern to our earlier results.

Figure 6

Calibration performance of 4 strategies on the Rosenbrock function. Boxplots show the absolute error in normalized outputs across 50 calibration targets for Randomly-Simulate (RS), Optimally-Predict (OP), Simulate-then-Predict (StP), and Predict-then-Simulate (PtS). Experiments were conducted using a neural network metamodel trained on Rosenbrock function data with varying nonlinearity (via parameters a∈{5,10,15}, b∈{10,20,30}). Results replicate the patterns observed in our main case study, confirming the robustness of PtS in consistently achieving the lowest calibration error.

Furthermore, because the NN metamodel is trained across multiple (a, b) contexts, this setup allows us to assess how well the metamodel-guided calibration generalizes to unseen contexts as well, that is, new (a, b) values not included in the training data. This feature mirrors real-world applications in which calibration must often operate in new or evolving environments. The graphs in Figure 7 confirm that while in-sample results are tighter, the model still captures the general behavior of the function in unseen (a, b) contexts. This also highlights an advantage over Gaussian process models, which typically struggle with extrapolation beyond the training domain.

Figure 7

In-sample versus out-of-sample calibration accuracy across varying Rosenbrock contexts. Scatter plots compare true versus optimized values for x and y under both in-sample (a∈{5,15,25}, b∈{10,20,30}) and out-of-sample (a∈{10,20}, b∈{15,25}) Rosenbrock parameter settings. In-sample predictions (top row) show tight alignment with the ground truth, whereas out-of-sample predictions (bottom row) exhibit slightly more variance but still preserve the overall functional structure. These results demonstrate that the neural network metamodel generalizes well across unseen calibration contexts.

While this approach strengthens the internal validity of our findings, we acknowledge that further work involving a broader set of synthetic and real-world functions is necessary for more comprehensive generalization. Nonetheless, the Rosenbrock experiments provide useful evidence that our proposed PtS method consistently outperforms alternative strategies in a controlled and analytically grounded setting, even when applied across varied simulation contexts.

Extension to Calibration with Stochastic Simulation Outputs

We also evaluated our approach under uncertainty by adding noise to the Rosenbrock function and training the metamodel on multiple noisy realizations. Predictive uncertainty was estimated using Monte Carlo dropout as an extension. A detailed description and additional results are provided in Appendix F.

Discussion

Computational results of this study present 3 key findings. First, Figure 4 shows that the PtS hybrid strategy has a better performance compared with the base strategy OP. This result shows that relying solely on metamodel-based calibration (like many studies in the literature do) can be problematic, even if the underlying metamodel has a high accuracy. On the other hand, OP-like metamodel-based strategies seem to be able to place us at the vicinity of the optimally calibrated point. Using this insight, a limited computational budget for simulation runs can be used in a smart way to get as close as possible to the optimally calibrated point. A similar idea underlies the StP strategy, in which simulation runs are used first to find the best actual output and then the metamodel is used to interpolate nearby points. In that sense, PtS and StP represent 2 ends of a single-step adaptive learning framework,²³ in which the emulator is used either before or after the simulation runs to guide selection. Future work could extend this idea into multistep adaptive calibration, dynamically alternating between emulator updates and targeted simulations based on evolving uncertainty or acquisition criteria.

Second, Figure 5 shows the change of calibration performance with respect to the change in the proximity range. We observed that there is a sweet spot of proximity range to look for. Although this sweet spot may change from application to application, it is worth keeping this in mind for the calibration procedures designed based on the PtS strategy. We anticipate that as the metamodel accuracy improves, the smaller proximity ranges will provide better results, and at the limit (when metamodel perfectly approximates the underlying simulation model), the anchor point found by the OP will be guaranteed to be the optimally calibrated point.

Third, the additional Rosenbrock experiments suggest that PtS is not only effective in structured case studies but also robust in more abstract and analytically challenging settings. Importantly, its consistent advantage in both in-sample and out-of-sample contexts indicates that the method adapts well to new parameter regimes, a property particularly valuable in real-world calibration tasks where model environments evolve over time. It is important to note, however, that the out-of-sample evaluations primarily reflect interpolation within context ranges represented in the training data rather than unrestricted extrapolation beyond the training domain. As with any surrogate-based approach, performance may deteriorate when applied far outside regions covered during metamodel training. In such cases, the hybrid PtS framework mitigates this risk by validating candidate solutions through additional simulation runs. These findings reinforce the practical promise of PtS as a reliable calibration strategy across varied applications.

Our approach also differs meaningfully from other NN-based calibration methods in the literature, such as BayCANN,¹⁷ which use NNs as emulators within a Bayesian framework and rely on sampling-based inference (e.g., Markov Chain Monte Carlo) to characterize posterior distributions of model parameters. While such approaches provide rich information about parameter uncertainty that is valuable for decision making under uncertainty, they can be computationally intensive, particularly in high-dimensional settings. In contrast, we treat the NN as a deterministic mapping and frame the calibration task as a mixed-integer linear program, allowing us to find a solution that minimizes the metamodel-predicted error. These approaches are therefore best viewed as complementary: our method is well suited to settings with deterministic or point-target calibration and tight computational budgets, whereas Bayesian approaches are more appropriate when full characterization of parameter uncertainty is required.

The focus of this study is on calibrating simulation models with deterministic outputs to fixed, known target values. That is, we assume the calibration targets are observed without uncertainty, and the input parameters are free to vary across their feasible ranges, without being constrained by prior probability distributions. This is effectively equivalent, in a decision-analytic sense, to assuming uniform (possibly improper) priors on the input parameters and point estimates for calibration targets. While this assumption is consistent with many applications in operations research and engineering, it contrasts with probabilistic decision-analytic models often used in health policy, where input parameters are specified by joint probability distributions and calibration targets may also be uncertain. In this context, “probabilistic” refers to uncertainty in inputs (and sometimes targets), not stochastic simulation noise in outputs. For a broader discussion of these distinctions, see Menzies et al.¹⁰ and Pineda-Antunez et al.¹¹ Extending our metamodel-based framework to accommodate parameter priors and uncertain targets, such as through integration with Bayesian calibration techniques,^8,9 is an important area for future research. We believe our work provides a complementary tool for deterministic or point-target calibration problems, particularly in settings requiring frequent, fast recalibration under changing scenarios or limited computational budgets.

A limitation of our study is that our NN-based metamodel does not propagate prediction uncertainty or explicitly model output stochasticity. While this is appropriate for the FluTE case study, in which outputs are treated as expected values and Monte Carlo variability is minimized, it may not generalize to settings where stochastic behavior persists. In such cases, calibration may require matching distributions rather than point estimates. Emulators that support uncertainty quantification, such as Gaussian processes,^24–26 may be more suitable, although they often face scalability challenges in high-dimensional or data-rich environments where deterministic ANNs remain a practical alternative. To assess the impact of stochastic simulation outputs, we conducted an additional experiment using a stochastic Rosenbrock function (see Appendix F). The results suggest that while central calibration accuracy is preserved, increased output variability can reduce calibration stability, highlighting the need for uncertainty-aware emulators and motivating future work in this direction. Another limitation concerns nonidentifiability in calibration, in which multiple parameter combinations may yield indistinguishably similar model outputs, leading to alternative optima in the OP formulation. To explore this, we examined multiple high-quality solutions returned by the optimization model across 100 calibration targets. In most cases, solutions with similar objective values exhibited small parameter differences, indicating stable calibration; however, in some cases (particularly near parameter boundaries), substantially different solutions with equivalent objective values were observed. This highlights the need for caution when interpreting calibrated parameters, especially in settings where identifiability and parameter stability are critical (see Appendix G for details). Finally, we do not explicitly account for the simulation runs required to train the metamodel when evaluating computational cost. In particular, methods such as OP and PtS assume a pretrained metamodel, which may not hold in all settings. While this assumption is reasonable in applications involving repeated or periodic calibration, it may be less appropriate in time-sensitive contexts. In these settings, metamodels can be trained offline to enable rapid calibration, but generating sufficiently rich training data can still be computationally demanding. Future work could explore more efficient data generation or incremental updating strategies.

For practitioners working with computationally expensive simulation models, our findings suggest a structured workflow: 1) train a metamodel using a diverse set of simulation runs, 2) identify a candidate input vector by minimizing the metamodel-predicted calibration error, and 3) allocate a limited number of additional simulation runs in a local neighborhood around this candidate solution, as in the PtS strategy. This approach avoids exhaustive parameter sweeps while maintaining fidelity to the true simulation outputs. The method is particularly advantageous in moderately high-dimensional settings or when repeated recalibration is required under tight computational budgets. Implementation details and example code are publicly available (https://github.com/elifsena000/SimulationCalibration).

Conclusions

Simulation metamodeling has become increasingly prominent as computational models grow more complex and expensive to evaluate. NNs, in particular, offer flexible surrogate representations capable of approximating highly nonlinear input–output relationships. However, calibrating simulation inputs to match target outputs remains computationally demanding, especially in high-dimensional settings.

In this study, we proposed a metamodel-based calibration framework that embeds trained NNs within an MIO model to identify the input parameter combination that minimizes metamodel-predicted calibration error. Building on this foundation, we introduced a hybrid strategy (PtS) that combines optimization over the metamodel with targeted simulation refinement. Computational experiments demonstrate that PtS substantially improves calibration accuracy relative to commonly used random search approaches while maintaining computational efficiency. The proposed framework makes no case-specific assumptions and is applicable to a wide range of deterministic simulation models. We believe this hybrid calibration strategy provides a practical and computationally efficient tool for simulation studies requiring repeated or resource-constrained calibration.

Supplemental Material

sj-docx-1-mdm-10.1177_0272989X261452258 – Supplemental material for A Metamodel-Based General-Purpose Autocalibration Tool for Simulation Models

Supplemental material, sj-docx-1-mdm-10.1177_0272989X261452258 for A Metamodel-Based General-Purpose Autocalibration Tool for Simulation Models by Taghi Khaniyev, Elif Sena Işık, Jagpreet Chhatwal, İsmail Fatih Yıldırım and Turgay Ayer in Medical Decision Making

Footnotes

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Taghi Khaniyev and Elif Sena Işık were funded by The Scientific and Technological Research Council of Türkiye (TUBITAK) under the International Fellowship Program (BIDEB-2232) for Early Stage Researchers (project No. 121C155) for their research including the preparation of this article. However, all responsibility regarding the publication belongs to the authors of the publication. Financial support received from TUBITAK does not mean that the content of the publication is scientifically approved by TUBITAK. Financial support for this study was provided in part by a grant from The Scientific and Technological Research Council of Türkiye (TUBITAK). The funding agreement ensured the authors’ independence in designing the study, interpreting the data, writing, and publishing the report.

Ethical Considerations

No ethical considerations. We confirm that the article uses only simulated data and does not use any data that require ethical consideration.

ORCID iDs

Taghi Khaniyev

Jagpreet Chhatwal

Data Availability

Research data and materials will be made available upon request to the corresponding author, Taghi Khaniyev.

References

Van Gelder

Das

Janssen

Roels

. Comparative study of metamodelling techniques in building energy simulation: guidelines for practitioners. Simulation Modelling Practice and Theory. 2014;49:245–57. DOI: 10.1016/j.simpat.2014.10.004

Patwary

Huang

HK.

Metamodel-based calibration of large-scale multimodal microscopic traffic simulation. Transp Res Part C: Emerg Technol. 2021;124:102859. DOI: 10.1016/j.trc.2020.102859

Malleson

Calibration of Simulation Models. New York: Springer; 2014. p 243–52. DOI: 10.1007/978-1-4614-5690-2_688

Santner

Williams

Notz

WI.

The Design and Analysis of Computer Experiments. New York: Springer; 2003.

Reddy

TA.

Literature review on calibration of building energy simulation programs: uses, problems, procedure, uncertainty, and tools. ASHRAE Trans. 2006;112(1):226–40.

Lim

Zhai

ZJ.

Comprehensive evaluation of the influence of meta-models on Bayesian calibration. Energy Build. 2017;155:66–75. DOI: 10.1016/j.enbuild.2017.09.009

Chen

Gao

Zeng

Liu

A meta-model-based optimization approach for fast and reliable calibration of building energy models. Energy. 2019;188:116046. DOI: 10.1016/j.energy.2019.116046

Kennedy

O’Hagan

Bayesian calibration of computer models. J Royal Stat Soc B: Stat Methodol. 2001;63(3):425–64. DOI: 10.1111/1467-9868.00294

Bayarri

Berger

Paulo

, et al. A framework for validation of computer models. Technometrics. 2007;49(2):138–54. DOI: 10.1198/004017007000000092

10.

Menzies

Soeteman

Pandya

Kim

JJ.

Bayesian methods for calibrating health policy models: a tutorial. Pharmacoeconomics. 2017;35(6):613–24. DOI: 10.1007/s40273-017-0507-4

11.

Pineda-Antunez

Levy

Owens

Goldhaber-Fiebert

JD.

Efficient surrogate modeling methods for Bayesian calibration of complex simulation models. Med Decis Making. 2023;43(3):355–68. DOI: 10.1177/0272989X231153018

12.

Myren

Lawrence

A comparison of Gaussian processes and neural networks for computer model emulation and calibration. Stat Anal Data Min ASA Data Sci J. 2021;14(6):606–23.

13.

Barton

Meckesheimer

Metamodel-based simulation optimization. In: Henderson

Nelson

, eds. Simulation, Handbooks in Operations Research and Management Science. Vol. 13. Amsterdam: Elsevier; 2006. p 535–74. DOI: 10.1016/S0927-0507(06)13018-2

14.

Fonseca

Navaresse

Moynihan

Simulation metamodeling through neural networks. Eng Appl Artif Intell. 2003;16:177–83. DOI: 10.1016/s0952-1976(03)00043-5

15.

Dunke

Nickel

Neural networks for the metamodeling of simulation models with online decision making. Simul Model Pract Theory. 2020;99:102016. DOI: 10.1016/j.simpat.2019.102016

16.

Fonseca

Navaresse

Neural networks for job shop simulation. Adv Eng Inform. 2002;16(4):241–46. DOI: 10.1016/S1474-0346(03)00005-3

17.

Jalal

Trikalinos

Alarid-Escudero

Baycann: streamlining Bayesian calibration with neural network metamodeling. Front Physiol. 2021;12:662314. DOI: 10.3389/fphys.2021.662314

18.

Laguna

Marti

Neural network prediction in a system for optimizing simulations. IIE Trans. 2002;34:273–82. DOI: 10.1023/A:1012485416856

19.

Vahdat

Alagoz

Chen

Saoud

Borah

Limburg

PJ.

Calibration and validation of the Colorectal Cancer and Adenoma Incidence and Mortality (CRC-AIM) microsimulation model using deep neural networks. Med Decis Making. 2023;43(6):719–36. DOI: 10.1177/0272989X231184175

20.

Anderson

Huchette

Tjandraatmadja

Vielma

JP.

Strong mixed-integer programming formulations for trained neural networks. Math Program 2020; 183(1): 3–9.

21.

Chao

Halloran

Obenchain

Longini

Jr.

FluTE, a publicly available stochastic influenza epidemic simulation model. PLoS Comput Biol. 2010;6:e1000656. DOI: 10.1371/journal.pcbi.1000656

22.

Shang

Y-W

Qiu

Y-H.

A note on the extended Rosenbrock function. Evol Comput. 2006;14(1):119–26. DOI: 10.1162/106365606776022733

23.

Ellis

Iskandar

Schmid

Wong

Trikalinos

TA.

Active learning for efficiently training emulators of computationally expensive mathematical models. Stat Med. 2020;39(18):3521–48. DOI: 10.1002/sim.8679

24.

Jones

DR.

Efficient global optimization of expensive black-box functions. J Glob Optim. 1998;13(4):455–92. DOI: 10.1023/A:1008306431147

25.

Oakley

O’Hagan

Probabilistic sensitivity analysis of complex models: a Bayesian approach. J R Stat Soc B Stat Methodol. 2002;64(3):751–69. DOI: 10.1111/1467-9868.00346

26.

Gramacy

RB.

Surrogates: Gaussian Process Modeling, Design, and Optimization for the Applied Sciences. Boca Raton (FL): CRC Press; 2020.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.50 MB

0.00 MB