Adrian Sandu

Jostein Barry-Straume, Changmin Son, Adrian Sandu et al.

Engine Health Management (EHM) depends on reliable forecasting of Remaining Useful Life (RUL) and on tracking thermal indicators such as turbine gas temperature (TGT). In practice, real-world fleet data are heterogeneous and non-stationary, and point predictions alone are insufficient for risk-aware maintenance decisions. This paper presents a multi-task scientific machine learning framework for turbine prognostics that jointly predicts turbine gas temperature untrimmed (TGTU), Delta Turbine Gas Temperature (DTGT), and RUL, with quantified uncertainty in the form of prediction intervals whose empirical coverage is evaluated. A shared sequence encoder (convolutional front-end with residual bidirectional LSTM layers and attention pooling) feeds task-specific heads, including mean--variance estimation for probabilistic regression and, optionally, a survival head for threshold-based event modeling. The framework is designed to be tunable via a small set of practitioner-facing parameters (e.g., DTGT thresholding rules and RUL target construction) so that deployment can align with in-house policies and proprietary criteria. The predictive performance of the proposed framework is evaluated using both point and interval metrics, including mean absolute error (MAE), prediction interval coverage probability (PICP), mean prediction interval width (MPIW), and the coverage--width criterion (CWC). Results are reported both in aggregate and stratified by flight phase and maintenance segment to highlight operational-context effects and to support uncertainty-aware monitoring.

Elias D. Nino-Ruiz, Adrian Sandu, Jeffrey Anderson

We present a practical implementation of the ensemble Kalman (EnKF) filter based on an iterative Sherman-Morrison formula. The new direct method exploits the special structure of the ensemble-estimated error covariance matrices in order to efficiently solve the linear systems involved in the analysis step of the EnKF. The computational complexity of the proposed implementation is equivalent to that of the best EnKF implementations available in the literature when the number of observations is much larger than the number of ensemble members. Even when this conditions is not fulfilled, the proposed method is expected to perform well since it does not employ matrix decompositions. Computational experiments using the Lorenz 96 and the oceanic quasi-geostrophic models are performed in order to compare the proposed algorithm with EnKF implementations that use matrix decompositions. In terms of accuracy, the results of all implementations are similar. The proposed method is considerably faster than other EnKF variants, even when the number of observations is large relative to the number of ensemble members.

Elias D. Nino, Adrian Sandu, Xinwei Deng

This paper discusses an efficient parallel implementation of the ensemble Kalman filter based on the modified Cholesky decomposition. The proposed implementation starts with decomposing the domain into sub-domains. In each sub-domain a sparse estimation of the inverse background error covariance matrix is computed via a modified Cholesky decomposition; the estimates are computed concurrently on separate processors. The sparsity of this estimator is dictated by the conditional independence of model components for some radius of influence. Then, the assimilation step is carried out in parallel without the need of inter-processor communication. Once the local analysis states are computed, the analysis sub-domains are mapped back onto the global domain to obtain the analysis ensemble. Computational experiments are performed using the Atmospheric General Circulation Model (SPEEDY) with the T-63 resolution on the Blueridge cluster at Virginia Tech. The number of processors used in the experiments ranges from 96 to 2,048. The proposed implementation outperforms in terms of accuracy the well-known local ensemble transform Kalman filter (LETKF) for all the model variables. The computational time of the proposed implementation is similar to that of the parallel LETKF method (where no covariance estimation is performed). Finally, for the largest number of processors, the proposed parallel implementation is 400 times faster than the serial version of the proposed method.

Adrian Sandu

Differential equations arising in many practical applications are characterized by multiple time scales. Multirate time integration seeks to solve them efficiently by discretizing each scale with a different, appropriate time step, while ensuring the overall accuracy and stability of the numerical solution. In a seminal paper Knoth and Wolke (APNUM, 1998) proposed a hybrid solution approach: discretize the slow component with an explicit Runge-Kutta method, and advance the fast component via a modified fast differential equation. The idea led to the development of multirate infinitesimal step (MIS) methods by Wensch et al. (BIT, 2009.)Günther and Sandu (BIT, 2016) explained MIS schemes as a particular case of multirate General-structure Additive Runge-Kutta (MR-GARK) methods. The hybrid approach offers extreme flexibility in the choice of the numerical solution process for the fast component. This work constructs a family of multirate infinitesimal GARK schemes (MRI-GARK) that extends the hybrid dynamics approachin multiple ways. Order conditions theory and stability analyses are developed, and practical explicit and implicit methods of up to order four are constructed. Numerical results confirm the theoretical findings. We expect the new MRI-GARK family to be most useful for systems of equations with widely disparate time scales, where the fast process is dispersive, and where the influence of the fast component on the slow dynamics is weak.

Jostein Barry-Straume, Changmin Son, Adrian Sandu et al.

Effective prognostics and health management of modern engines relies on accurate turbine gas temperature predictions and robust uncertainty quantification to ensure reliability and safety. This paper investigates five major approaches for constructing prediction intervals -- namely the Delta method, Bayesian Monte Carlo Dropout, Bootstrap method, Lower-Upper Bound Estimation, and Mean-Variance Estimation -- as a means of capturing the uncertainty in neural network predictions of turbine gas temperature. Each approach is implemented within a unified experimental framework that employs cross-validation for hyperparameter selection, repeated train-test splits for performance robustness, and multiple metrics to evaluate both the accuracy and tightness of the intervals. In particular, Coverage Probability, Normalized Mean Prediction Interval Width, and the Coverage Width-based Criterion are measured to comprehensively assess each method's reliability and sharpness. Experiments conducted on a representative turbine gas temperature dataset reveal distinct trade-offs among the five methods in terms of interval coverage, width, and stability. These findings provide a practical guide for selecting and tuning prediction interval methods in engine health management and prognostics, ensuring both interpretability and precision in real-world applications.

Ahmed Attia, Razvan Stefanescu, Adrian Sandu

Hybrid Monte-Carlo (HMC) sampling smoother is a fully non-Gaussian four-dimensional data assimilation algorithm that works by directly sampling the posterior distribution formulated in the Bayesian framework. The smoother in its original formulation is computationally expensive due to the intrinsic requirement of running the forward and adjoint models repeatedly. Here we present computationally efficient versions of the HMC sampling smoother based on reduced-order approximations of the underlying model dynamics. The schemes developed herein are tested numerically using the shallow-water equations model on Cartesian coordinates. The results reveal that the reduced-order versions of the smoother are capable of accurately capturing the posterior probability density, while being significantly faster than the original full order formulation.

Azam Moosavi, Razvan Stefanescu, Adrian Sandu

This paper introduces multivariate input-output models to predict the errors and bases dimensions of local parametric Proper Orthogonal Decomposition reduced-order models. We refer to these multivariate mappings as the MP-LROM models. We employ Gaussian Processes and Artificial Neural Networks to construct approximations of these multivariate mappings. Numerical results with a viscous Burgers model illustrate the performance and potential of the machine learning based regression MP-LROM models to approximate the characteristics of parametric local reduced-order models. The predicted reduced-order models errors are compared against the multi-fidelity correction and reduced order model error surrogates methods predictions, whereas the predicted reduced-order dimensions are tested against the standard method based on the spectrum of snapshots matrix. Since the MP-LROM models incorporate more features and elements to construct the probabilistic mappings they achieve more accurate results. However, for high-dimensional parametric spaces, the MP-LROM models might suffer from the curse of dimensionality. Scalability challenges of MP-LROM models and the feasible ways of addressing them are also discussed in this study.

Arash Sarshar, Steven Roberts, Adrian Sandu

Multirate time integration methods apply different step sizes to resolve different components of the system based on the local activity levels. This local selection of step sizes allows increased computational efficiency while achieving the desired solution accuracy. While the multirate idea is elegant and has been around for decades, multirate methods are not yet widely used in applications. This is due, in part, to the difficulties raised by the construction of high order multirate schemes. Seeking to overcome these challenges, this work focuses on the design of practical high-order multirate methods using the theoretical framework of generalized additive Runge-Kutta (MrGARK) methods, which provides the generic order conditions and the linear and nonlinear stability analyses. A set of design criteria for practical multirate methods is defined herein: method coefficients should be generic in the step size ratio, but should not depend strongly on this ratio; unnecessary coupling between the fast and the slow components should be avoided; and the step size controllers should adjust both the micro- and the macro-steps. Using these criteria, we develop MrGARK schemes of up to order four that are explicit-explicit (both the fast and slow component are treated explicitly), implicit-explicit (implicit in the fast component and explicit in the slow one), and explicit-implicit (explicit in the fast component and implicit in the slow one). Numerical experiments illustrate the performance of these new schemes.

Mahesh Narayanamurthi, Paul Tranquilli, Adrian Sandu et al.

Exponential integrators are special time discretization methods where the traditional linear system solves used by implicit schemes are replaced with computing the action of matrix exponential-like functions on a vector. A very general formulation of exponential integrators is offered by the Exponential Propagation Iterative methods of Runge-Kutta type (EPIRK) family of schemes. The use of Jacobian approximations is an important strategy to drastically reduce the overall computational costs of implicit schemes while maintaining the quality of their solutions. This paper extends the EPIRK class to allow the use of inexact Jacobians as arguments of the matrix exponential-like functions. Specifically, we develop two new families of methods: EPIRK-W integrators that can accommodate any approximation of the Jacobian, and EPIRK-K integrators that rely on a specific Krylov-subspace projection of the exact Jacobian. Classical order conditions theories are constructed for these families. A practical EPIRK-W method of order three and an EPIRK-K method of order four are developed. Numerical experiments indicate that the methods proposed herein are computationally favorable when compared to existing exponential integrators.

Jostein Barry-Straume, Arash Sarshar, Andrey A. Popov et al.

A fundamental problem in science and engineering is designing optimal control policies that steer a given system towards a desired outcome. This work proposes Control Physics-Informed Neural Networks (Control PINNs) that simultaneously solve for a given system state, and for the optimal control signal, in a one-stage framework that conforms to the underlying physical laws. Prior approaches use a two-stage framework that first models and then controls a system in sequential order. In contrast, a Control PINN incorporates the required optimality conditions in its architecture and in its loss function. The success of Control PINNs is demonstrated by solving the following open-loop optimal control problems: (i) an analytical problem, (ii) a one-dimensional heat equation, and (iii) a two-dimensional predator-prey problem.

Azam Moosavi, Vishwas Rao, Adrian Sandu

Complex numerical weather prediction models incorporate a variety of physical processes, each described by multiple alternative physical schemes with specific parameters. The selection of the physical schemes and the choice of the corresponding physical parameters during model configuration can significantly impact the accuracy of model forecasts. There is no combination of physical schemes that works best for all times, at all locations, and under all conditions. It is therefore of considerable interest to understand the interplay between the choice of physics and the accuracy of the resulting forecasts under different conditions. This paper demonstrates the use of machine learning techniques to study the uncertainty in numerical weather prediction models due to the interaction of multiple physical processes. The first problem addressed herein is the estimation of systematic model errors in output quantities of interest at future times, and the use of this information to improve the model forecasts. The second problem considered is the identification of those specific physical processes that contribute most to the forecast uncertainty in the quantity of interest under specified meteorological conditions. The discrepancies between model results and observations at past times are used to learn the relationships between the choice of physical processes and the resulting forecast errors. Numerical experiments are carried out with the Weather Research and Forecasting (WRF) model. The output quantity of interest is the model precipitation, a variable that is both extremely important and very challenging to forecast. The physical processes under consideration include various micro-physics schemes, cumulus parameterizations, short wave, and long wave radiation schemes. The experiments demonstrate the strong potential of machine learning approaches to aid the study of model errors.

Andrey A. Popov, Arash Sarshar, Austin Chennault et al.

A rapidly growing area of research is the use of machine learning approaches such as autoencoders for dimensionality reduction of data and models in scientific applications. We show that the canonical formulation of autoencoders suffers from several deficiencies that can hinder their performance. Using a meta-learning approach, we reformulate the autoencoder problem as a bi-level optimization procedure that explicitly solves the dimensionality reduction task. We prove that the new formulation corrects the identified deficiencies with canonical autoencoders, provide a practical way to solve it, and showcase the strength of this formulation with a simple numerical illustration.

Ali Haisam Muhammad Rafid, Adrian Sandu

Deep neural networks (DNN) have found wide applicability in numerous fields due to their ability to accurately learn very complex input-output relations. Despite their accuracy and extensive use, DNNs are highly susceptible to adversarial attacks due to limited generalizability. For future progress in the field, it is essential to build DNNs that are robust to any kind of perturbations to the data points. In the past, many techniques have been proposed to robustify DNNs using first-order derivative information of the network. This paper proposes a new robustification approach based on control theory. A neural network architecture that incorporates feedback control, named Feedback Neural Networks, is proposed. The controller is itself a neural network, which is trained using regular and adversarial data such as to stabilize the system outputs. The novel adversarial training approach based on the feedback control architecture is called Feedback Looped Adversarial Training (FLAT). Numerical results on standard test problems empirically show that our FLAT method is more effective than the state-of-the-art to guard against adversarial attacks.

Paul Tranquilli, Adrian Sandu, Hong Zhang

This paper develops a new class of linearly implicit time integration schemes called Linearly-Implicit Runge-Kutta-W (LIRK-W) methods. These schemes are based on an implicit-explicit approach which does not require a splitting of the right hand side and allow for arbitrary, time dependent, and stage varying approximations of the linear systems appearing in the method. Several formulations of LIRK-W schemes, each designed for specific approximation types, and their associated order condition theories are presented.

Abhinab Bhattacharjee, Andrey A. Popov, Arash Sarshar et al.

The Adam optimizer, often used in Machine Learning for neural network training, corresponds to an underlying ordinary differential equation (ODE) in the limit of very small learning rates. This work shows that the classical Adam algorithm is a first-order implicit-explicit (IMEX) Euler discretization of the underlying ODE. Employing the time discretization point of view, we propose new extensions of the Adam scheme obtained by using higher-order IMEX methods to solve the ODE. Based on this approach, we derive a new optimization algorithm for neural network training that performs better than classical Adam on several regression and classification problems.

Jostein Barry-Straume, Adwait D. Verulkar, Arash Sarshar et al.

The objective of designing a control system is to steer a dynamical system with a control signal, guiding it to exhibit the desired behavior. The Hamilton-Jacobi-Bellman (HJB) partial differential equation offers a framework for optimal control system design. However, numerical solutions to this equation are computationally intensive, and analytical solutions are frequently unavailable. Knowledge-guided machine learning methodologies, such as physics-informed neural networks (PINNs), offer new alternative approaches that can alleviate the difficulties of solving the HJB equation numerically. This work presents a multistage ensemble framework to learn the optimal cost-to-go, and subsequently the corresponding optimal control signal, through the HJB equation. Prior PINN-based approaches rely on a stabilizing the HJB enforcement during training. Our framework does not use stabilizer terms and offers a means of controlling the nonlinear system, via either a singular learned control signal or an ensemble control signal policy. Success is demonstrated in closed-loop control, using both ensemble- and singular-control, of a steady-state time-invariant two-state continuous nonlinear system with an infinite time horizon, accounting of noisy, perturbed system states and varying initial conditions.

Ali Haisam Muhammad Rafid, Adrian Sandu

Regularization techniques such as $\mathcal{L}_1$ and $\mathcal{L}_2$ regularizers are effective in sparsifying neural networks (NNs). However, to remove a certain neuron or channel in NNs, all weight elements related to that neuron or channel need to be prunable, which is not guaranteed by traditional regularization. This paper proposes a simple new approach named "Guided Regularization" that prioritizes the weights of certain NN units more than others during training, which renders some of the units less important and thus, prunable. This is different from the scattered sparsification of $\mathcal{L}_1$ and $\mathcal{L}_2$ regularizers where the the components of a weight matrix that are zeroed out can be located anywhere. The proposed approach offers a natural reduction of NN in the sense that a model is being trained while also neutralizing unnecessary units. We empirically demonstrate that our proposed method is effective in pruning NNs while maintaining performance.

Austin Chennault, Andrey A. Popov, Amit N. Subrahmanya et al.

Data assimilation is the process of fusing information from imperfect computer simulations with noisy, sparse measurements of reality to obtain improved estimates of the state or parameters of a dynamical system of interest. The data assimilation procedures used in many geoscience applications, such as numerical weather forecasting, are variants of the our-dimensional variational (4D-Var) algorithm. The cost of solving the underlying 4D-Var optimization problem is dominated by the cost of repeated forward and adjoint model runs. This motivates substituting the evaluations of the physical model and its adjoint by fast, approximate surrogate models. Neural networks offer a promising approach for the data-driven creation of surrogate models. The accuracy of the surrogate 4D-Var solution depends on the accuracy with each the surrogate captures both the forward and the adjoint model dynamics. We formulate and analyze several approaches to incorporate adjoint information into the construction of neural network surrogates. The resulting networks are tested on unseen data and in a sequential data assimilation problem using the Lorenz-63 system. Surrogates constructed using adjoint information demonstrate superior performance on the 4D-Var data assimilation problem compared to a standard neural network surrogate that uses only forward dynamics information.

Rachel Cooper, Andrey A. Popov, Adrian Sandu

Reduced order modeling (ROM) is a field of techniques that approximates complex physics-based models of real-world processes by inexpensive surrogates that capture important dynamical characteristics with a smaller number of degrees of freedom. Traditional ROM techniques such as proper orthogonal decomposition (POD) focus on linear projections of the dynamics onto a set of spectral features. In this paper we explore the construction of ROM using autoencoders (AE) that perform nonlinear projections of the system dynamics onto a low dimensional manifold learned from data. The approach uses convolutional neural networks (CNN) to learn spatial features as opposed to spectral, and utilize a physics informed (PI) cost function in order to capture temporal features as well. Our investigation using the quasi-geostrophic equations reveals that while the PI cost function helps with spatial reconstruction, spatial features are less powerful than spectral features, and that construction of ROMs through machine learning-based methods requires significant investigation into novel non-standard methodologies.

Andrey A Popov, Adrian Sandu

Data assimilation is a Bayesian inference process that obtains an enhanced understanding of a physical system of interest by fusing information from an inexact physics-based model, and from noisy sparse observations of reality. The multifidelity ensemble Kalman filter (MFEnKF) recently developed by the authors combines a full-order physical model and a hierarchy of reduced order surrogate models in order to increase the computational efficiency of data assimilation. The standard MFEnKF uses linear couplings between models, and is statistically optimal in case of Gaussian probability densities. This work extends MFEnKF to work with non-linear couplings between the models. Optimal nonlinear projection and interpolation operators are obtained by appropriately trained physics-informed autoencoders, and this approach allows to construct reduced order surrogate models with less error than conventional linear methods. Numerical experiments with the canonical Lorenz '96 model illustrate that nonlinear surrogates perform better than linear projection-based ones in the context of multifidelity filtering.

Andrey A Popov, Adrian Sandu

Ever since its inception, the Ensemble Kalman Filter has elicited many heuristic methods that sought to correct it. One such method is localization---the thought that `nearby' variables should be highly correlated with `far away' variable not. Recognizing that correlation is a time-dependent property, adaptive localization is a natural extension to these heuristics. We propose a Bayesian approach to adaptive Schur-product localization for the DEnKF, and extend it to support multiple radii of influence. We test both the empirical validity of (multivariate) adaptive localization, and of our approach. We test a simple toy problem (Lorenz'96), extending it to a multivariate model, and a more realistic geophysical problem (1.5 Layer Quasi-Geostrophic). We show that the multivariate approach has great promise on the toy problem, and that the univariate approach leads to improved filter performance for the realistic geophysical problem.

Azam Moosavi, Ahmed Attia, Adrian Sandu

Data assimilation plays a key role in large-scale atmospheric weather forecasting, where the state of the physical system is estimated from model outputs and observations, and is then used as initial condition to produce accurate future forecasts. The Ensemble Kalman Filter (EnKF) provides a practical implementation of the statistical solution of the data assimilation problem and has gained wide popularity as. This success can be attributed to its simple formulation and ease of implementation. EnKF is a Monte-Carlo algorithm that solves the data assimilation problem by sampling the probability distributions involved in Bayes theorem. Because of this, all flavors of EnKF are fundamentally prone to sampling errors when the ensemble size is small. In typical weather forecasting applications, the model state space has dimension $10^{9}-10^{12}$, while the ensemble size typically ranges between $30-100$ members. Sampling errors manifest themselves as long-range spurious correlations and have been shown to cause filter divergence. To alleviate this effect covariance localization dampens spurious correlations between state variables located at a large distance in the physical space, via an empirical distance-dependent function. The quality of the resulting analysis and forecast is greatly influenced by the choice of the localization function parameters, e.g., the radius of influence. The localization radius is generally tuned empirically to yield desirable results.This work, proposes two adaptive algorithms for covariance localization in the EnKF framework, both based on a machine learning approach. The first algorithm adapts the localization radius in time, while the second algorithm tunes the localization radius in both time and space. Numerical experiments carried out with the Lorenz-96 model, and a quasi-geostrophic model, reveal the potential of the proposed machine learning approaches.

Ahmed Attia, Azam Moosavi, Adrian Sandu

This paper presents a fully non-Gaussian version of the Hamiltonian Monte Carlo (HMC) sampling filter. The Gaussian prior assumption in the original HMC filter is relaxed. Specifically, a clustering step is introduced after the forecast phase of the filter, and the prior density function is estimated by fitting a Gaussian Mixture Model (GMM) to the prior ensemble. Using the data likelihood function, the posterior density is then formulated as a mixture density, and is sampled using a HMC approach (or any other scheme capable of sampling multimodal densities in high-dimensional subspaces). The main filter developed herein is named "cluster HMC sampling filter" (ClHMC). A multi-chain version of the ClHMC filter, namely MC-ClHMC is also proposed to guarantee that samples are taken from the vicinities of all probability modes of the formulated posterior. The new methodologies are tested using a quasi-geostrophic (QG) model with double-gyre wind forcing and bi-harmonic friction. Numerical results demonstrate the usefulness of using GMMs to relax the Gaussian prior assumption in the HMC filtering paradigm.

Azam Moosavi, Razvan Stefanescu, Adrian Sandu

Reduced order models are computationally inexpensive approximations that capture the important dynamical characteristics of large, high-fidelity computer models of physical systems. This paper applies machine learning techniques to improve the design of parametric reduced order models. Specifically, machine learning is used to develop feasible regions in the parameter space where the admissible target accuracy is achieved with a predefined reduced order basis, to construct parametric maps, to chose the best two already existing bases for a new parameter configuration from accuracy point of view and to pre-select the optimal dimension of the reduced basis such as to meet the desired accuracy. By combining available information using bases concatenation and interpolation as well as high-fidelity solutions interpolation we are able to build accurate reduced order models associated with new parameter settings. Promising numerical results with a viscous Burgers model illustrate the potential of machine learning approaches to help design better reduced order models.

Angelamaria Cardone, Zdzislaw Jackiewicz, Hong Zhang et al.

For many systems of differential equations modeling problems in science and engineering, there are natural splittings of the right hand side into two parts, one non-stiff or mildly stiff, and the other one stiff. For such systems implicit-explicit (IMEX) integration combines an explicit scheme for the non-stiff part with an implicit scheme for the stiff part. In a recent series of papers two of the authors (Sandu and Zhang) have developed IMEX GLMs, a family of implicit-explicit schemes based on general linear methods. It has been shown that, due to their high stage order, IMEX GLMs require no additional coupling order conditions, and are not marred by order reduction. This work develops a new extrapolation-based approach to construct practical IMEX GLM pairs of high order. We look for methods with large absolute stability region, assuming that the implicit part of the method is A- or L-stable. We provide examples of IMEX GLMs with optimal stability properties. Their application to a two dimensional test problem confirms the theoretical findings.