Standard Monte Carlo Simulation Research Articles

Smooth Nested Simulation: Bridging Cubic and Square Root Convergence Rates in High Dimensions

Nested simulation concerns estimating functionals of a conditional expectation via simulation. In this paper, we propose a new method based on kernel ridge regression to exploit the smoothness of the conditional expectation as a function of the multidimensional conditioning variable. Asymptotic analysis shows that the proposed method can effectively alleviate the curse of dimensionality on the convergence rate as the simulation budget increases, provided that the conditional expectation is sufficiently smooth. The smoothness bridges the gap between the cubic root convergence rate (that is, the optimal rate for the standard nested simulation) and the square root convergence rate (that is, the canonical rate for the standard Monte Carlo simulation). We demonstrate the performance of the proposed method via numerical examples from portfolio risk management and input uncertainty quantification. This paper was accepted by Baris Ata, stochastic models and simulation. Funding: The authors acknowledge financial support from the National Natural Science Foundation of China [Grant NSFC 12101149] and the Hong Kong Research Grants Council [Grants GRF 17201520 and GRF 17206821]. Supplemental Material: The e-companion and data files are available at https://doi.org/10.1287/mnsc.2022.00204 .

Multi-staged ERO2.0 simulation of material erosion and deposition in recessed ITER mirror assemblies

In reactor class devices, with long pulse, high fluence plasma operation and with diagnostic systems that cannot easily be maintained, a critical issue is the deposition/erosion of wall materials and its impact on key optical diagnostic components — mirrors. A workflow for recessed volume simulation was developed for the three-dimensional Monte-Carlo code ERO2.0 (Romazanov, 2020), consisting of multiple simulation stages starting from first wall erosion to more focussed volumes around the key components. This complex workflow has been found to be mandatory to have meaningful statistics of interaction processes on the mirror because the physical fluxes of particles into recessed volumes are extremely low compared to first wall fluxes, leading to unsatisfying statistics in the case of standard Monte-Carlo simulations. The workflow is introduced with the example of a diagnostic mirror assembly in equatorial port 12 behind panels 14 and 15 in the diagnostic first wall in ITER, assuming a beryllium first wall of the main vessel, a stainless steel mirror unit and two molybdenum mirrors embedded in a port plug for neutron shielding. The simulation case corresponds to an H-mode burning plasma case from SOLPS modelling, for which the incoming flux distributions of deuterium neutrals into the recessed volume were simulated by the neutrals transport code EIRENE.

Interval-Probabilistic Electricity-Heat-Gas Flow Calculation by Dual-Level Surrogate Structure

The study of highly efficient methods for performing multiple energy flow (MEF) computation with epistemic and aleatory uncertainties, known as interval-probabilistic energy flow (IPEF) calculation, is a fast-growing area of research. However, existing approaches are easily caught in a quandary among accuracy, efficiency, and scope of applicability. This paper presents a dual-level surrogate structure (DLSS) for dealing with the IPEF calculation of an integrated energy system (IES) in an accurate and efficient manner. The structure accommodates electricity, heat, and gas, attaining a high-level balance among precision, speed, and range of applicability. At the lower level of DLSS, a sparse polynomial chaos expansion (sPCE) model is employed to directly map random input variables to output boundary values, saving massive repetition of deterministic calculation. At the upper level of DLSS, a heterogeneous-learning-based multiple energy flow (HL-MEF) model is proposed to accelerate the indispensable deterministic MEF calculation for online sampling, thereby easing the computational burden of sPCE training. Finally, numerical results on a real-world electricity-heat-gas system demonstrate that the proposed DLSS can solve the IPEF problem at nearly 400 times the speed of the standard double-layer Monte Carlo simulation (DLMCS) with guaranteed accuracy.

Learning the random variables in Monte Carlo simulations with stochastic gradient descent: Machine learning for parametric PDEs and financial derivative pricing

AbstractIn financial engineering, prices of financial products are computed approximately many times each trading day with (slightly) different parameters in each calculation. In many financial models such prices can be approximated by means of Monte Carlo (MC) simulations. To obtain a good approximation the MC sample size usually needs to be considerably large resulting in a long computing time to obtain a single approximation. A natural deep learning approach to reduce the computation time when new prices need to be calculated as quickly as possible would be to train an artificial neural network (ANN) to learn the function which maps parameters of the model and of the financial product to the price of the financial product. However, empirically it turns out that this approach leads to approximations with unacceptably high errors, in particular when the error is measured in the ‐norm, and it seems that ANNs are not capable to closely approximate prices of financial products in dependence on the model and product parameters in real life applications. This is not entirely surprising given the high‐dimensional nature of the problem and the fact that it has recently been proved for a large class of algorithms, including the deep learning approach outlined above, that such methods are in general not capable to overcome the curse of dimensionality for such approximation problems in the ‐norm. In this article we introduce a new numerical approximation strategy for parametric approximation problems including the parametric financial pricing problems described above and we illustrate by means of several numerical experiments that the introduced approximation strategy achieves a very high accuracy for a variety of high‐dimensional parametric approximation problems, even in the ‐norm. A central aspect of the approximation strategy proposed in this article is to combine MC algorithms with machine learning techniques to, roughly speaking, learn the random variables (LRV) in MC simulations. In other words, we employ stochastic gradient descent (SGD) optimization methods not to train parameters of standard ANNs but instead to learn random variables appearing in MC approximations. In that sense, the proposed LRV strategy has strong links to Quasi‐Monte Carlo (QMC) methods as well as to the field of algorithm learning. Our numerical simulations strongly indicate that the LRV strategy might indeed be capable to overcome the curse of dimensionality in the ‐norm in several cases where the standard deep learning approach has been proven not to be able to do so. This is not a contradiction to the established lower bounds mentioned above because this new LRV strategy is outside of the class of algorithms for which lower bounds have been established in the scientific literature. The proposed LRV strategy is of general nature and not only restricted to the parametric financial pricing problems described above, but applicable to a large class of approximation problems. In this article we numerically test the LRV strategy in the case of the pricing of European call options in the Black‐Scholes model with one underlying asset, in the case of the pricing of European worst‐of basket put options in the Black‐Scholes model with three underlying assets, in the case of the pricing of European average put options in the Black‐Scholes model with three underlying assets and knock‐in barriers, as well as in the case of stochastic Lorentz equations. For these examples the LRV strategy produces highly convincing numerical results when compared with standard MC simulations, QMC simulations using Sobol sequences, SGD‐trained shallow ANNs, and SGD‐trained deep ANNs.

Strong error analysis of Euler methods for overdamped generalized Langevin equations with fractional noise: Nonlinear case

This paper considers the strong error analysis of the Euler and fast Euler methods for nonlinear overdamped generalized Langevin equations driven by the fractional noise. The main difficulty lies in handling the interaction between the fractional Brownian motion and the singular kernel, which is overcome by means of the Malliavin calculus and fine estimates of several multiple singular integrals. Consequently, these two methods are proved to be strongly convergent with order nearly min{2(H + α − 1), α}, where H ∈ (1/2, 1) and α ∈ (1 − H, 1) respectively characterize the singularity levels of fractional noises and singular kernels in the underlying equation. This result improves the existing convergence order H + α − 1 of Euler methods for the nonlinear case, and gives a positive answer to the open problem raised in Fang and Li [ESAIM Math. Model. Numer. Anal. 54 (2020) 431–463]. As an application of the theoretical findings, we further investigate the complexity of the multilevel Monte Carlo simulation based on the fast Euler method, which turns out to behave better performance than the standard Monte Carlo simulation when computing the expectation of functionals of the considered equation. Finally, numerical experiments are carried out to support the theoretical results.

In Vitro Bioaccessibility and Health Risk of Heavy Metals from PM2.5/PM10 in Arid Areas—Hotan City, China

The impact of heavy metals in particulates plays an assignable role in human health—especially in the northwest region of China, which is affected by severe dust storms—and the bioaccessibility and health risks of heavy metals in particulate matter have not yet been quantified and evaluated. This study used Gamble’s solution and PBET (physiologically based extraction test) experiments to simulate the human respiratory and digestive systems, and analyzed the concentrations and bioaccessibility of Pb, Mn, Ni, Cd and As in PM2.5 and PM10 samples in the urban area of Hotan City during summer (July) and winter (January). The result shows that Mn and Pb are the most abundant elements in five metals. The bioaccessibility of Mn in gastric fluid was the highest in both summer (PM2.5: 64.1%, PM10: 52.0%) and winter (PM2.5: 88.0%, PM10: 85.1%). Meanwhile, in the respiratory system, the highest bioaccessibility of PM2.5 and PM10 in summer was Ni (53.3%), and Pb (47.9%), respectively. Although the concentration of Cd is low in winter, its bioaccessibility in lung fluid was the highest (PM2.5: 74.7%, PM10: 62.3%). The USEPA standard model and Monte Carlo simulation results show that the heavy metals in PM2.5 and PM10 would give rise to non-carcinogenic risk for both adults and children through the respiratory system in summer but had little risk in winter. However, the metal may have non-carcinogenic risk to children through intake. In addition, there is a cancer risk to adults through the respiratory system in winter (PM2.5: CR = 1.80 × 10−6, PM10: CR = 2.82 × 10−6), while there is a carcinogenic risk through the digestive system regardless of season and age.

Multifidelity conditional value-at-risk estimation by dimensionally decomposed generalized polynomial chaos-Kriging

We propose novel methods for Conditional Value-at-Risk (CVaR) estimation for nonlinear systems under high-dimensional dependent random inputs. We develop a novel DD-GPCE-Kriging surrogate that merges dimensionally decomposed generalized polynomial chaos expansion and Kriging to accurately approximate nonlinear and nonsmooth random outputs. We use DD-GPCE-Kriging (1) for Monte Carlo simulation (MCS) and (2) within multifidelity importance sampling (MFIS). The MCS-based method samples from DD-GPCE-Kriging, which is efficient and accurate for high-dimensional dependent random inputs, yet introduces bias. Thus, we propose an MFIS-based method where DD-GPCE-Kriging determines the biasing density, from which we draw a few high-fidelity samples to provide an unbiased CVaR estimate. To accelerate the biasing density construction, we compute DD-GPCE-Kriging using a cheap-to-evaluate low-fidelity model. Numerical results for mathematical functions show that the MFIS-based method is more accurate than the MCS-based method when the output is nonsmooth. The scalability of the proposed methods and their applicability to complex engineering problems are demonstrated on a two-dimensional composite laminate with 28 (partly dependent) random inputs and a three-dimensional composite T-joint with 20 (partly dependent) random inputs. In the former, the proposed MFIS-based method achieves 104x speedup compared to standard MCS using the high-fidelity model, while accurately estimating CVaR with 1.15% error.

Study on angle-only relative navigation for unmanned aerial vehicle formation in GPS-denied environment

To satisfy the autonomous relative navigation demand of unmanned aerial vehicle (UAV) formation in a GPS-denied environment, an angle-only relative navigation algorithm is proposed. First, the relative states in 3-dimension space are decoupled into relative motion along the horizontal plane and altitude by applying the relative azimuth between the north-east-down global frame and the body-fixed horizontal frame belonging to UAVs. Second, a relative dynamic and measurement model considering the measurement deviation of sensors are established. This solves the problem of colored noise in a majority of sensors used in practical engineering. Third, a consensus constraint model of relative states is established by introducing a geometric constraint among formation members. Additionally, a consensus-extended Kalman filter is designed. Finally, the angle-only relative navigation method is validated using standard Monte Carlo simulations. The simulation results show that all the states of the system are observable when the observability criterion is satisfied. The relative position accuracy achieved using gyros with a drift rate of 50°/h, accelerometers with 1 mg deviation, and optical cameras with 1° Gaussian white noise is within 10 m (3σ). Additionally, the method proposed in this study can be used to estimate the sensor measurement deviations effectively. Furthermore, it can be potentially applied to improve the absolute navigation accuracy of UAVs.

An efficient stratified sampling scheme for the simultaneous estimation of small failure probabilities in wind engineering applications

With the advent of performance-based engineering frameworks for extreme winds, greater efficiency in the estimation of failure probabilities associated with multiple limit states at ultimate load levels, including collapse, is required. To meet this demand, an optimal stratified sampling-based simulation scheme is proposed that shares several meritorious characteristics of standard Monte Carlo simulation (MCS), including the capability to deal with high-dimensional uncertain spaces and multiple limit state functions. The optimality of the scheme is governed by the allocation of MCS samples among the strata as well as the strata division itself. A constrained optimization problem is formulated to address efficient sample allocation which is an extension of the Neyman allocation. Details for the practical implementation of the scheme within wind engineering reliability assessment applications are discussed. It is observed that the proposed methodology can reduce the variance by several orders of magnitude compared to MCS, which is illustrated on a wind-excited steel structure for the limit states of member/system first yield and collapse. The applicability of Latin hypercube sampling within stratified sampling-based simulation and the resulting efficiency gains are also critically studied. The results show promise in the wider applicability of the proposed scheme in performance-based wind engineering.

Seismic collaborative reliability evaluation of slopes using subset simulation via support vector machine

The computational effort places a limit on the low failure probability estimate of geotechnical engineering exposed to seismic stimulation. In this research, a novel method is proposed for the seismic reliability assessment of complicated geotechnical engineering slope structures with low failure probability. First, an improved subset simulation (SS) approach is provided for the low-failure events with high-dimensional indirect cross-correlation parameters. This method expands the usefulness of SS by using the support vector machine (SVM) method. Then, in order to reduce computing effort, the seismic collaborative reliability evaluation (SCRE) framework for seismic reliability assessment based on SS and the response conditioning method (RCM) is constructed. Finally, the accuracy of SCRE is verified by comparison with the standard Monte Carlo simulation based on 2D slope. The results show that SCRE's computational effort is approximately 2% that of Monte Carlo simulation. Meanwhile, SVM does a good job of preserving the significant indirect cross-correlation that occurs in the subspace samples that were split by SS. The time-saving benefit of SCRE is further enhanced when the aforementioned situations are expanded to 3D, and the findings demonstrate that the failure probability and sliding surface of the slope are different in 2D and 3D.

Model two-particle kinetic equation for pairs of quasiparticles

Because the main process of interaction in the Boltzmann model of real gases is the binary collision of molecules, it is convenient to use the two-particle kinetic equation to describe the dynamics of a rarefied gas. This equation was written using the same physical assumptions as those used by Ludwig Boltzmann. The right-hand side of this equation contains the product of the linear scattering operator and chaos projector. The Boltzmann equation follows from this equation without any additional approximations after simple integration of the velocities and positions of the second particle. Using the divergence form of the scattering operator, this equation can be represented as the Liouville equation, which implies that real molecules can be replaced by quasiparticles whose distribution function is the same as that of real molecules but whose dynamics are completely different. Pairs of quasiparticles do not collide but move along continuous trajectories in the phase space. The relative velocities in pairs of quasiparticles slowly rotate with an angular velocity vector depending on the distribution function. We provide an explicit approximate expression for the angular velocity through the first few velocity moments, using a special covariant Grad expansion for the velocity distribution function, which reduces to the exact Bobylev–Kruk–Wu solution in the isotropic case. We simulated the relaxation of distribution function to equilibrium and compared results with the existing exact solutions. The described algorithm will be effective for modeling flow regions with low Knudsen numbers, where the standard Direct Monte Carlo Simulation (DSMC) method encounters significant difficulties.

Temperature Effects in Flexible Adsorption Processes for Amorphous Microporous Polymers.

A collection of atomistic molecular simulations is reported that illustrate the impact of adsorption temperature on species uptake and adsorbate-induced structural rearrangement for amorphous polymers of intrinsic microporosity. Temperature-sensitive structural rearrangement is evaluated by contrasting two methods: standard grand canonical Monte Carlo simulations using a rigid framework approximation and a combined Monte Carlo/molecular dynamics approach that fully incorporates framework flexibility. We report single-component gas phase adsorption isotherms for CH4, C2H4, C2H6, C3H6, C3H8, and CO2 across a temperature range of 250-400 K for models of an archetypal polymer of intrinsic microporosity, PIM-1. A quadratic model is presented that captures two main mechanisms of temperature-dependent adsorption-induced deformation of PIM-1 up to a relative swelling of 1.15: thermal expansion and an increased propensity to swell as a function of species uptake. Two case studies are reported that highlight the critical role of operating temperature in industrial storage and separation applications. The first study focuses on methane storage and delivery applications using a pressure-temperature swing adsorption application (PTSA). We demonstrate that larger working capacities are accompanied by increased volumetric strain between adsorption-desorption steps. The second case study considers PIM-1 as an adsorbent to separate an exemplar ternary syngas mixture at operating temperatures ranging 300-550 K. A temperature threshold of ∼400 K is identified, beyond which adsorption-induced PIM-1 swelling is negligible and the solubility selectivity-loading curve transitions to exhibiting a nearly linear relationship.

Pseudoparticle vertex solver for quantum impurity models

We present a quantum impurity solver based on a pseudo-particle framework, which combines diagrammatic resummations for a three-point vertex with diagrammatic Monte Carlo sampling of a four-point vertex. This recently proposed approach [A. J. Kim et al., arXiv:2112.15549] is generalized here to fermionic impurity problems and we discuss the technical details of the implementation, including the time-stepping approach, the Monte Carlo updates, and the routines for checking the two-particle irreducibility of the four-point vertex. We also explain how the vertex information can be efficiently stored using a Dubiner basis representation. The convergence properties of the algorithm are demonstrated with applications to exactly solvable impurity models and dynamical mean field theory simulations of the single-orbital Hubbard model. It is furthermore shown that the algorithm can handle a two-orbital problem with off-diagonal hybridizations, which would cause a severe sign problem in standard hybridization-expansion Monte Carlo simulations. Since the vertex-based algorithm successfully handles sign oscillating integrals in equilibrium and samples only connected diagrams, it may be a promising approach for real-time simulations.

A Position‐Free Path Integral for Homogeneous Slabs and Multiple Scattering on Smith Microfacets

AbstractWe consider the problem of multiple scattering on Smith microfacets. This problem is equivalent to computing volumetric light transport in a homogeneous slab. Although the symmetry of the slab allows for significant simplification, fully analytic solutions are scarce and not general enough for most applications. Standard Monte Carlo simulation, although general, is expensive and leads to variance that must be dealt with.We present the first unbiased, truly position‐free path integral for evaluating the BSDF of a homogeneous slab. We collapse the spatially‐1D path integral of previous works to a position‐free form using an analytical preintegration of collision distances. Evaluation of the resulting path integral, which now contains only directions, reduces to simple recursive manipulation of exponential distributions. Applying Monte Carlo to solve the reduced integration problem leads to lower variance.Our new algorithm allows us to render multiple scattering on Smith microfacets with less variance than prior work, and, in the case of conductors, to do so without any bias. Additionally, our algorithm can also be used to accelerate the rendering of BSDFs containing volumetrically scattering layers, at reduced variance compared to standard Monte Carlo integration.

Classically emulated digital quantum simulation of the Schwinger model with a topological term via adiabatic state preparation

We designed a protocol for digital quantum computation of a gauge theory with a topological term in Minkowski spacetime, which is practically inaccessible by standard lattice Monte Carlo simulations. We focus on $1+1$ dimensional quantum electrodynamics with the $\ensuremath{\theta}$ term known as the Schwinger model and test our protocol for this on an IBM simulator. We construct the true vacuum state of a lattice Schwinger model using adiabatic state preparation which, in turn, allows us to compute an expectation value of the fermion mass operator with respect to the vacuum. Upon taking a continuum limit we find that our result in the massless case agrees with the known exact result. In the massive case, we find an agreement with mass perturbation theory in the small-mass regime and deviations in the large-mass regime. We estimate computational costs required to take a reasonable continuum limit. Our results imply that digital quantum simulation appears a promising tool to explore nonperturbative aspects of gauge theories with real time and topological terms.

Thermodynamic calculations using reverse Monte Carlo: convergence aspects, sources of error and guidelines for improving accuracy

ABSTRACT The reverse Monte Carlo (RMC) method has been widely used to gain 3D structural ordering information from experimental scattering data. Recently, we have introduced a novel application of RMC, viz., calculating thermodynamic properties of crystalline materials, construction of phase diagram and rapid estimation of the local structural order [Agrahari and Chatterjee, Physical Review E, 104, 044129 (2021).]. The method has been shown to be accurate and orders-of-magnitude faster than standard Monte Carlo simulations – this makes the approach quite promising. However, the error in RMC has never been systematically quantified. The goal here is to perform a thorough investigation into the types of error in RMC-based thermodynamic calculations, assess the relative magnitude of these errors, and develop strategies to improve the accuracy. Some of the conclusions presented are also relevant to previous experiment-based RMC implementations.

The density–temperature range of exchange–correlation exciton existence by the fermionic path integral Monte Carlo method

A recently developed fermionic path integral Monte Carlo approach [Filinov et al., Phys. Rev. E 102, 033203 (2020) and Filinov et al., J. Phys. A 55, 035001 (2021)] has been applied for the estimation of the density–temperature range of exchange–correlation exciton existence in a strongly coupled degenerate uniform electron gas. The approach allows us to reduce the “fermionic sign problem” taking into account the interference effects of the Coulomb and exchange interaction of electrons in the basic Monte Carlo cell and its periodic images. Our results for radial distribution functions demonstrate the formation and decay of a short-range quantum ordering of electrons associated with exchange–correlation excitons in the literature. Such excitons have never been observed earlier in standard path integral Monte Carlo simulations.

Estimation of steady-state quantities of an HMM with some rarely generated emissions

Abstract We propose to apply the importance sampling and the antithetic variates statistical techniques to estimate steady-state quantities of an Hidden Markov chain (HMM) of which certain emissions are rarely generated. Compared to standard Monte Carlo simulation, the use of these techniques, allow a significant reduction in simulation time. Numerical Monte Carlo examples are studied to show the usefulness and efficiency of the proposed approach.

Stochastic finite element analysis of composite cycloidal propeller blade during crash-stop ship maneuver

This article quantifies uncertainty in the dynamic responses of marine cycloidal propeller (MCP) blade by Stochastic Finite Element Method (SFEM) during crash-stop ship maneuver. Stopping a ship in an emergency is known as crash-stop astern maneuver. A MCP contains a horizontal circular metallic disc on which six aerofoil blades are hung vertically. All blades are made up of carbon fiber-reinforced composite laminates. Uncertainty in ship speed and composite material properties are propagated to find out the statistical response of natural frequencies, displacement, and stresses on the blade. SFEM is a combination of Stochastic Response Surface Method (SRSM) and Finite Element Method (FEM). SRSM is a collocation-based non-intrusive polynomial chaos expansion to transform any random distribution into a stochastic response surface equation (surrogate model). During ship maneuvering motion, the whole unit rotates about a vertical axis at the centre. For each rotation of the disc, all blades undergo periodic oscillations. The SRSM for uncertainty propagation is computationally more efficient than the standard Monte Carlo Simulation (MCS) technique without compromising the accuracy of the results. It is observed that uncertainty in vibration responses of MCP blade is significant and worth study.

Bound States of the Exchange—Correlation Excitons in the Uniform Electron Gas by the Monte Carlo Simulations

The modified path integral representation of Wigner functions and the new Monte Carlo approach has been suggested to account for the impact of the interparticle interaction on the Pauli exclusion principle of fermions. This approach also allows to calculate the momentum distribution functions and to reduce the “sign problem” that is inaccessible to the standard path integral Monte Carlo methods. The obtained pair electron–electron distribution functions for the “uniform electron gas” demonstrate the short-range quantum ordering of electrons associated with exchange-correlation excitons. The exchange-correlation exciton is caused by the interaction of electrons with positively charged exchange holes and the excluded volume effect. The developed approach allows one to study the density–temperature range of the exciton arising, existence, and decay. Using the potential of the mean force and semiclassical Bohr–Sommerfeld quantization condition, we have demonstrated the existence of bound states disturbing the Maxwellian distribution and estimated their average energy levels. The exchange-correlation excitons have not been observed earlier in the standard path integral Monte Carlo simulations.