Search

When there is a fault, the protection system automatically removes one or more transmission lines on a fast time scale of less than one minute. The outaged lines form a pattern in the transmission network. We extract these patterns from utility outage data, determine some key statistics of these patterns, and then show how to generate new patterns consistent with these statistics. The generated patterns provide a new and easily feasible way to model the overall effect of the protection system at the scale of a large transmission system. This new generative modeling of protection is expected to contribute to simulations of disturbances in large grids so that they can better quantify the risk of blackouts. Analysis of the pattern sizes suggests an index that describes how much outages spread in the transmission network at the fast timescale.

Power grid expansion planning requires making large investment decisions in the present that will impact the future cost and reliability of a system exposed to wide-ranging uncertainties. Extreme temperatures can pose significant challenges to providing power by increasing demand and decreasing supply and have contributed to recent major power outages. We propose to address a modeling challenge of such high-impact, low-frequency events with a bi-objective stochastic integer optimization model that finds solutions with different trade-offs between efficiency in normal conditions and risk to extreme events. We propose a conditional sampling approach paired with a risk measure to address the inherent challenge in approximating the risk of low-frequency events within a sampling based approach. We present a model for spatially correlated, county-specific temperatures and a method to generate both unconditional and conditionally extreme temperature samples from this model efficiently. These models are investigated within an extensive case study with realistic data that demonstrates the effectiveness of the bi-objective approach and the conditional sampling technique. We find that spatial correlations in the temperature samples are essential to finding good solutions and that modeling generator temperature dependence is an important consideration for finding efficient, low-risk solutions., Comment: 43 pages, 10 figures

This paper explores two condensed-space interior-point methods to efficiently solve large-scale nonlinear programs on graphics processing units (GPUs). The interior-point method solves a sequence of symmetric indefinite linear systems, or Karush-Kuhn-Tucker (KKT) systems, which become increasingly ill-conditioned as we approach the solution. Solving a KKT system with traditional sparse factorization methods involve numerical pivoting, making parallelization difficult. A solution is to condense the KKT system into a symmetric positive-definite matrix and solve it with a Cholesky factorization, stable without pivoting. Although condensed KKT systems are more prone to ill-conditioning than the original ones, they exhibit structured ill-conditioning that mitigates the loss of accuracy. This paper compares the benefits of two recent condensed-space interior-point methods, HyKKT and LiftedKKT. We implement the two methods on GPUs using MadNLP.jl, an optimization solver interfaced with the NVIDIA sparse linear solver cuDSS and with the GPU-accelerated modeler ExaModels.jl. Our experiments on the PGLIB and the COPS benchmarks reveal that GPUs can attain up to a tenfold speed increase compared to CPUs when solving large-scale instances.

This paper demonstrates the scalability of open-source GPU-accelerated nonlinear programming (NLP) frameworks -- ExaModels.jl and MadNLP.jl -- for solving multi-period alternating current (AC) optimal power flow (OPF) problems on GPUs with high memory capacities (e.g., NVIDIA GH200 with 480 GB of unified memory). There has been a growing interest in solving multi-period AC OPF problems, as the increasingly fluctuating electricity market requires operation planning over multiple periods. These problems, formerly deemed intractable, are now becoming technologically feasible to solve thanks to the advent of high-memory GPU hardware and accelerated NLP tools. This study evaluates the capability of these tools to tackle previously unsolvable multi-period AC OPF instances. Our numerical experiments, run on an NVIDIA GH200, demonstrate that we can solve a multi-period OPF instance with more than 10 million variables up to $10^{-4}$ precision in less than 10 minutes. These results demonstrate the efficacy of the GPU-accelerated NLP frameworks for the solution of extreme-scale multi-period OPF. We provide ExaModelsPower.jl, an open-source modeling tool for multi-period AC OPF models for GPUs.

The theory of dual control was introduced more than seven decades ago. Although it has provided rich insights to the fields of control, estimation, and system identification, dual control is generally computationally prohibitive. In recent years, however, the use of Koopman operator theory for control applications has been emerging. This paper presents a new reformulation of the stochastic optimal control problem that, employing the Koopman operator, yields a standard LQR problem with the dual control as its solution. We provide a numerical example that demonstrates the effectiveness of the proposed approach compared with certainty equivalence control, when applied to systems with varying observability., Comment: 6 pages

In this paper, we address the challenge of solving large-scale graph-structured nonlinear programs (gsNLPs) in a scalable manner. GsNLPs are problems in which the objective and constraint functions are associated with nodes on a graph and depend on the variables of adjacent nodes. This graph-structured formulation encompasses various specific instances, such as dynamic optimization, PDE-constrained optimization, multistage stochastic optimization, and general network optimization. By leveraging the sequential quadratic programming (SQP) framework, we propose a globally convergent overlapping graph decomposition method to solve large-scale gsNLPs under standard mild regularity conditions on the graph topology. In each iteration, we perform an overlapping graph decomposition to compute an approximate Newton direction in a parallel environment. Then, we select a suitable stepsize and update the primal-dual iterate by performing a backtracking line search on an exact augmented Lagrangian merit function. Built on the exponential decay of sensitivity of gsNLPs, we show that the approximate Newton direction is a descent direction of the augmented Lagrangian, which leads to global convergence with a local linear convergence rate. In particular, global convergence is achieved for sufficiently large overlaps, and the local linear convergence rate improves exponentially in terms of the overlap size. Our results match existing state-of-the-art guarantees established for dynamic programs (which simply correspond to linear graphs). We validate the theory on a semilinear elliptic PDE-constrained problem., Comment: 3 figures

We propose a data-driven approach for propagating uncertainty in stochastic power grid simulations and apply it to the estimation of transmission line failure probabilities. A reduced-order equation governing the evolution of the observed line energy probability density function is derived from the Fokker--Planck equation of the full-order continuous Markov process. Our method consists of estimates produced by numerically integrating this reduced equation. Numerical experiments for scalar- and vector-valued energy functions are conducted using the classical multimachine model under spatiotemporally correlated noise perturbation. The method demonstrates a more sample-efficient approach for computing probabilities of tail events when compared with kernel density estimation. Moreover, it produces vastly more accurate estimates of joint event occurrence when compared with independent models., Comment: 12 pages, 6 figures, and 1 table

Accurate dynamic modeling of power systems is essential to assess the stability of electrical power systems when faced with disturbances, which can trigger cascading failures leading to blackouts. A continuum model proves to be effective in capturing Electromechanical Wave (EMW) propagation characteristics, including its velocity, arrival time, and deviations. Analyzing these characteristics enables the assessment of the impacts of EMW on the performance of the protection system. Prior research has often modeled nonlinear EMW propagation through Partial Differential Equations (PDEs) within a homogeneous and uniform frame structure, assuming constant bus voltages across the entire power system. However, this assumption can produce inaccurate results. In this paper, we relax this assumption by introducing a second-order nonlinear hyperbolic EMW propagation equation model that accounts for voltage variations. Additionally, we present numerical solutions for the EMW propagation equation using the Lax-Wendroff integration method. To validate our approach, we conduct simulations on two test systems: a two-bus one-machine system and the New England 39-bus 10-machine system. The simulation results demonstrate the effectiveness of our proposed model and emphasize the importance of including the bus voltage equations in the analysis.

This paper presents a novel centralized, variational data assimilation approach for calibrating transient dynamic models in electrical power systems, focusing on load model parameters. With the increasing importance of inverter-based resources, assessing power systems' dynamic performance under disturbances has become challenging, necessitating robust model calibration methods. The proposed approach expands on previous Bayesian frameworks by establishing a posterior distribution of parameters using an approximation around the maximum a posteriori value. We illustrate the efficacy of our method by generating events of varying intensity, highlighting its ability to capture the systems' evolution accurately and with associated uncertainty estimates. This research improves the precision of dynamic performance assessments in modern power systems, with potential applications in managing uncertainties and optimizing system operations., Comment: 9 pages, 8 figures, and 1 table

This paper introduces a framework for solving alternating current optimal power flow (ACOPF) problems using graphics processing units (GPUs). While GPUs have demonstrated remarkable performance in various computing domains, their application in ACOPF has been limited due to challenges associated with porting sparse automatic differentiation (AD) and sparse linear solver routines to GPUs. We address these issues with two key strategies. First, we utilize a single-instruction, multiple-data abstraction of nonlinear programs. This approach enables the specification of model equations while preserving their parallelizable structure and, in turn, facilitates the parallel AD implementation. Second, we employ a condensed-space interior-point method (IPM) with an inequality relaxation. This technique involves condensing the Karush--Kuhn--Tucker (KKT) system into a positive definite system. This strategy offers the key advantage of being able to factorize the KKT matrix without numerical pivoting, which has hampered the parallelization of the IPM algorithm. By combining these strategies, we can perform the majority of operations on GPUs while keeping the data residing in the device memory only. Comprehensive numerical benchmark results showcase the advantage of our approach. Remarkably, our implementations -- MadNLP.jl and ExaModels.jl -- running on NVIDIA GPUs achieve an order of magnitude speedup compared with state-of-the-art tools running on contemporary CPUs., Comment: Accepted for publication in PSCC 2024

Measures of power grid vulnerability are often assessed by the amount of damage an adversary can exact on the network. However, the cascading impact of such attacks is often overlooked, even though cascades are one of the primary causes of large-scale blackouts. This paper explores modifications of transmission line protection settings as candidates for adversarial attacks, which can remain undetectable as long as the network equilibrium state remains unaltered. This forms the basis of a black-box function in a Bayesian optimization procedure, where the objective is to find protection settings that maximize network degradation due to cascading. Notably, our proposed method is agnostic to the choice of the cascade simulator and its underlying assumptions. Numerical experiments reveal that, against conventional wisdom, maximally misconfiguring the protection settings of all network lines does not cause the most cascading. More surprisingly, even when the degree of misconfiguration is limited due to resource constraints, it is still possible to find settings that produce cascades comparable in severity to instances where there are no resource constraints., Comment: 13 pages, 5 figures

Physics-based covariance models provide a systematic way to construct covariance models that are consistent with the underlying physical laws in Gaussian process analysis. The unknown parameters in the covariance models can be estimated using maximum likelihood estimation, but direct construction of the covariance matrix and classical strategies of computing with it requires $n$ physical model runs, $n^2$ storage complexity, and $n^3$ computational complexity. To address such challenges, we propose to approximate the discretized covariance function using hierarchical matrices. By utilizing randomized range sketching for individual off-diagonal blocks, the construction process of the hierarchical covariance approximation requires $O(\log{n})$ physical model applications and the maximum likelihood computations require $O(n\log^2{n})$ effort per iteration. We propose a new approach to compute exactly the trace of products of hierarchical matrices which results in the expected Fischer information matrix being computable in $O(n\log^2{n})$ as well. The construction is totally matrix-free and the derivatives of the covariance matrix can then be approximated in the same hierarchical structure by differentiating the whole process. Numerical results are provided to demonstrate the effectiveness, accuracy, and efficiency of the proposed method for parameter estimations and uncertainty quantification., Comment: Accepted by SIAM/ASA Journal on Uncertainty Quantification

Modern computational methods, involving highly sophisticated mathematical formulations, enable several tasks like modeling complex physical phenomenon, predicting key properties and design optimization. The higher fidelity in these computer models makes it computationally intensive to query them hundreds of times for optimization and one usually relies on a simplified model albeit at the cost of losing predictive accuracy and precision. Towards this, data-driven surrogate modeling methods have shown a lot of promise in emulating the behavior of the expensive computer models. However, a major bottleneck in such methods is the inability to deal with high input dimensionality and the need for relatively large datasets. With such problems, the input and output quantity of interest are tensors of high dimensionality. Commonly used surrogate modeling methods for such problems, suffer from requirements like high number of computational evaluations that precludes one from performing other numerical tasks like uncertainty quantification and statistical analysis. In this work, we propose an end-to-end approach that maps a high-dimensional image like input to an output of high dimensionality or its key statistics. Our approach uses two main framework that perform three steps: a) reduce the input and output from a high-dimensional space to a reduced or low-dimensional space, b) model the input-output relationship in the low-dimensional space, and c) enable the incorporation of domain-specific physical constraints as masks. In order to accomplish the task of reducing input dimensionality we leverage principal component analysis, that is coupled with two surrogate modeling methods namely: a) Bayesian hybrid modeling, and b) DeepHyper's deep neural networks. We demonstrate the applicability of the approach on a problem of a linear elastic stress field data., Comment: 17 pages, 16 figures, IDETC Conference Submission

Existing methods to determine the stability of a power system to small perturbations are based on eigenvalue analysis and focus on the asymptotic (long-term) behavior of the power grid. During the preasymptotic (short-term) transient, however, the system can exhibit large growth that is not explained by eigenvalues alone. In this paper we propose a new framework to determine the maximum (optimal) preasymptotic growth using the singular value decomposition. The approach is tailored to the analysis of quantities of interest in power system dynamics, such as the set of rotor speed deviations. Matrix-free techniques are developed to avoid the explicit formation of dense matrices and enable the analysis of large-scale systems without reaching memory bounds. Extensive results carried out from small to very large-scale systems (e.g., 70k-bus system) verify the theoretical aspects of the technique.

Gaussian processes (GPs) are an attractive class of machine learning models because of their simplicity and flexibility as building blocks of more complex Bayesian models. Meanwhile, graph neural networks (GNNs) emerged recently as a promising class of models for graph-structured data in semi-supervised learning and beyond. Their competitive performance is often attributed to a proper capturing of the graph inductive bias. In this work, we introduce this inductive bias into GPs to improve their predictive performance for graph-structured data. We show that a prominent example of GNNs, the graph convolutional network, is equivalent to some GP when its layers are infinitely wide; and we analyze the kernel universality and the limiting behavior in depth. We further present a programmable procedure to compose covariance kernels inspired by this equivalence and derive example kernels corresponding to several interesting members of the GNN family. We also propose a computationally efficient approximation of the covariance matrix for scalable posterior inference with large-scale data. We demonstrate that these graph-based kernels lead to competitive classification and regression performance, as well as advantages in computation time, compared with the respective GNNs., Comment: ICLR 2023. Code is available at https://github.com/niuzehao/gnn-gp

We investigate how to port the standard interior-point method to new exascale architectures for block-structured nonlinear programs with state equations. Computationally, we decompose the interior-point algorithm into two successive operations: the evaluation of the derivatives and the solution of the associated Karush-Kuhn-Tucker (KKT) linear system. Our method accelerates both operations using two levels of parallelism. First, we distribute the computations on multiple processes using coarse parallelism. Second, each process uses a SIMD/GPU accelerator locally to accelerate the operations using fine-grained parallelism. The KKT system is reduced by eliminating the inequalities and the state variables from the corresponding equations, to a dense matrix encoding the sensitivities of the problem's degrees of freedom, drastically minimizing the memory exchange. We demonstrate the method's capability on the supercomputer Polaris, a testbed for the future exascale Aurora system. Each node is equipped with four GPUs, a setup amenable to our two-level approach. Our experiments on the stochastic optimal power flow problem show that the method can achieve a 50x speed-up compared to the state-of-the-art method., Comment: 23 pages, 8 figures

This article presents a performance analysis for stochastic predictive control (SPC) in linear systems with quadratic performance index and additive and multiplicative uncertainties. Under a finite support assumption, the problem can be cast as a finite-dimensional quadratic program, but the problem becomes quickly intractable as the problem size grows exponentially in the horizon length. SPC aims to compute approximate solutions by solving a sequence of problems with truncated prediction horizons and committing the solution in a receding-horizon fashion. While this approach is widely used in practice, its performance relative to the optimal solution is not well understood. This article reports for the first time a rigorous performance guarantee of SPC: Under stabilizability and detectability conditions, the dynamic regret of SPC is exponentially small in the prediction horizon length, allowing SPC to achieve near-optimal performance at a substantially reduced computational expense.

This article presents a constrained policy optimization approach for the optimal control of systems under nonstationary uncertainties. We introduce an assumption that we call Markov embeddability that allows us to cast the stochastic optimal control problem as a policy optimization problem over the augmented state space. Then, the infinite-dimensional policy optimization problem is approximated as a finite-dimensional nonlinear program by applying function approximation, deterministic sampling, and temporal truncation. The approximated problem is solved by using automatic differentiation and condensed-space interior-point methods. We formulate several conceptual and practical open questions regarding the asymptotic exactness of the approximation and the solution strategies for the approximated problem. As a proof of concept, we provide a numerical example demonstrating the performance of the control policy obtained by the proposed method.

We report numerical results on solving constrained linear-quadratic model predictive control (MPC) problems by exploiting graphics processing units (GPUs). The presented method reduces the MPC problem by eliminating the state variables and applies a condensed-space interior-point method to remove the inequality constraints in the KKT system. The final condensed matrix is positive definite and can be efficiently factorized in parallel on GPU/SIMD architectures. In addition, the size of the condensed matrix depends only on the number of controls in the problem, rendering the method particularly effective when the problem has many states but few inputs and moderate horizon length. Our numerical results for PDE-constrained problems show that the approach is an order of magnitude faster than a standard CPU implementation. We also provide an open-source Julia framework that facilitates modeling (DynamicNLPModels.jl) and solution (MadNLP.jl) of MPC problems on GPUs.

Nonstationary Gaussian process models can capture complex spatially varying dependence structures in spatial datasets. However, the large number of observations in modern datasets makes fitting such models computationally intractable with conventional dense linear algebra. In addition, derivative-free or even first-order optimization methods can be slow to converge when estimating many spatially varying parameters. We present here a computational framework that couples an algebraic block-diagonal plus low-rank covariance matrix approximation with stochastic trace estimation to facilitate the efficient use of second-order solvers for maximum likelihood estimation of Gaussian process models with many parameters. We demonstrate the effectiveness of these methods by simultaneously fitting 192 parameters in the popular nonstationary model of Paciorek and Schervish using 107,600 sea surface temperature anomaly measurements.

The classification of structural phase transitions as displacive or order-disorder in character is usually based on spectroscopic data above the transition. We use single crystal x-ray diffraction to investigate structural correlations in the quasi-skutterudites, (Ca$_x$Sr$_{1-x}$)$_3$Rh$_4$Sn$_{13}$, which have a quantum phase transition at $x\sim0.9$. Three-dimensional pair distribution functions show that the amplitudes of local atomic displacements are temperature-independent below the transition and persist to well above the transition, a signature of order-disorder behavior. The implications for the associated electronic transitions are discussed., Comment: 5 pages, 4 figures

This paper studies the trade-off between the degree of decentralization and the performance of a distributed controller in a linear-quadratic control setting. We study a system of interconnected agents over a graph and a distributed controller, called $\kappa$-distributed control, which lets the agents make control decisions based on the state information within distance $\kappa$ on the underlying graph. This controller can tune its degree of decentralization using the parameter $\kappa$ and thus allows a characterization of the relationship between decentralization and performance. We show that under mild assumptions, including stabilizability, detectability, and a subexponentially growing graph condition, the performance difference between $\kappa$-distributed control and centralized optimal control becomes exponentially small in $\kappa$. This result reveals that distributed control can achieve near-optimal performance with a moderate degree of decentralization, and thus it is an effective controller architecture for large-scale networked systems.

The interior-point method (IPM) has become the workhorse method for nonlinear programming. The performance of IPM is directly related to the linear solver employed to factorize the Karush--Kuhn--Tucker (KKT) system at each iteration of the algorithm. When solving large-scale nonlinear problems, state-of-the art IPM solvers rely on efficient sparse linear solvers to solve the KKT system. Instead, we propose a novel reduced-space IPM algorithm that condenses the KKT system into a dense matrix whose size is proportional to the number of degrees of freedom in the problem. Depending on where the reduction occurs we derive two variants of the reduced-space method: linearize-then-reduce and reduce-then-linearize. We adapt their workflow so that the vast majority of computations are accelerated on GPUs. We provide extensive numerical results on the optimal power flow problem, comparing our GPU-accelerated reduced space IPM with Knitro and a hybrid full space IPM algorithm. By evaluating the derivatives on the GPU and solving the KKT system on the CPU, the hybrid solution is already significantly faster than the CPU-only solutions. The two reduced-space algorithms go one step further by solving the KKT system entirely on the GPU. As expected, the performance of the two reduction algorithms depends intrinsically on the number of available degrees of freedom: their performance is poor when the problem has many degrees of freedom, but the two algorithms are up to 3 times faster than Knitro as soon as the relative number of degrees of freedom becomes smaller., Comment: 21 pages, 3 figures, 2 tables

This paper presents an efficient method for extracting the second-order sensitivities from a system of implicit nonlinear equations on upcoming graphical processing units (GPU) dominated computer systems. We design a custom automatic differentiation (AutoDiff) backend that targets highly parallel architectures by extracting the second-order information in batch. When the nonlinear equations are associated to a reduced space optimization problem, we leverage the parallel reverse-mode accumulation in a batched adjoint-adjoint algorithm to compute efficiently the reduced Hessian of the problem. We apply the method to extract the reduced Hessian associated to the balance equations of a power network, and show on the largest instances that a parallel GPU implementation is 30 times faster than a sequential CPU reference based on UMFPACK., Comment: SIAM-PP22

Motivated by the increasing availability of high-performance parallel computing, we design a distributed parallel algorithm for linearly-coupled block-structured nonconvex constrained optimization problems. Our algorithm performs Jacobi-type proximal updates of the augmented Lagrangian function, requiring only local solutions of separable block nonlinear programming (NLP) problems. We provide a cheap and explicitly computable Lyapunov function that allows us to establish global and local sublinear convergence of our algorithm, its iteration complexity, as well as simple, practical and theoretically convergent rules for automatically tuning its parameters. This in contrast to existing algorithms for nonconvex constrained optimization based on the alternating direction method of multipliers that rely on at least one of the following: Gauss-Seidel or sequential updates, global solutions of NLP problems, non-computable Lyapunov functions, and hand-tuning of parameters. Numerical experiments showcase its advantages for large-scale problems, including the multi-period optimization of a 9000-bus AC optimal power flow test case over 168 time periods, solved on the Summit supercomputer using an open-source Julia code.

We propose a novel feasible-path algorithm to solve the optimal power flow (OPF) problem for real-time use cases. The method augments the seminal work of Dommel and Tinney with second-order derivatives to work directly in the reduced space induced by the power flow equations. In the reduced space, the optimization problem includes only inequality constraints corresponding to the operational constraints. While the reduced formulation directly enforces the physical constraints, the operational constraints are softly enforced through Augmented Lagrangian penalty terms. In contrast to interior-point algorithms (state-of-the art for solving OPF), our algorithm maintains feasibility at each iteration, which makes it suitable for real-time application. By exploiting accelerator hardware (Graphic Processing Units) to compute the reduced Hessian, we show that the second-order method is numerically tractable and is effective to solve both static and real-time OPF problems., Comment: 7 pages, 3 figures, conference proceeding

We study nonlinear optimization problems with a stochastic objective and deterministic equality and inequality constraints, which emerge in numerous applications including finance, manufacturing, power systems and, recently, deep neural networks. We propose an active-set stochastic sequential quadratic programming (StoSQP) algorithm that utilizes a differentiable exact augmented Lagrangian as the merit function. The algorithm adaptively selects the penalty parameters of the augmented Lagrangian and performs a stochastic line search to decide the stepsize. The global convergence is established: for any initialization, the KKT residuals converge to zero almost surely. Our algorithm and analysis further develop the prior work of Na et al., (2022). Specifically, we allow nonlinear inequality constraints without requiring the strict complementary condition; refine some of the designs in Na et al., (2022) such as the feasibility error condition and the monotonically increasing sample size; strengthen the global convergence guarantee; and improve the sample complexity on the objective Hessian. We demonstrate the performance of the designed algorithm on a subset of nonlinear problems collected in CUTEst test set and on constrained logistic regression problems., Comment: 65 pages, 9 figures

We propose a fast temporal decomposition procedure for solving long-horizon nonlinear dynamic programs. The core of the procedure is sequential quadratic programming (SQP) that utilizes a differentiable exact augmented Lagrangian as the merit function. Within each SQP iteration, we approximately solve the Newton system using an overlapping temporal decomposition strategy. We show that the approximate search direction is still a descent direction of the augmented Lagrangian, provided the overlap size and penalty parameters are suitably chosen, which allows us to establish the global convergence. Moreover, we show that a unit stepsize is accepted locally for the approximate search direction, and further establish a uniform, local linear convergence over stages. This local convergence rate matches the rate of the recent Schwarz scheme by Na et al., 2022. However, the Schwarz scheme has to solve nonlinear subproblems to optimality in each iteration, while we only perform a single Newton step instead. Numerical experiments validate our theories and demonstrate the superiority of our method., Comment: 41 pages, 1 figure

In this work we present a novel technique, based on a trust-region optimization algorithm and second-order trajectory sensitivities, to compute the extreme trajectories of power system dynamic simulations given a bounded set that represents parametric uncertainty. We show how this method, while remaining computationally efficient compared with sampling-based techniques, overcomes the limitations of previous sensitivity-based techniques to approximate the bounds of the trajectories, when the local approximation loses validity because of the nonlinearity. In addition, we show how this method can be adapted to account for those cases in which the initial conditions depend on the uncertain parameter. To conclude, we present several numerical experiments that showcase the accuracy and scalability of the technique, including a demonstration on the IEEE New England test system.

We present the implementation of a trust-region Newton algorithm ExaTron for bound-constrained nonlinear programming problems, fully running on multiple GPUs. Without data transfers between CPU and GPU, our implementation has achieved the elimination of a major performance bottleneck under a memory-bound situation, particularly when solving many small problems in batch. We discuss the design principles and implementation details for our kernel function and core operations. Different design choices are justified by numerical experiments. By using the application of distributed control of alternating current optimal power flow, where a large problem is decomposed into many smaller nonlinear programs using a Lagrangian approach, we demonstrate computational performance of ExaTron on the Summit supercomputer at Oak RidgeNational Laboratory. Our numerical results show the linear scaling with respect to the batch size and the number of GPUs and more than 35 times speedup on 6 GPUs than on 40 CPUs available on a single node.

Randomized algorithms have propelled advances in artificial intelligence and represent a foundational research area in advancing AI for Science. Future advancements in DOE Office of Science priority areas such as climate science, astrophysics, fusion, advanced materials, combustion, and quantum computing all require randomized algorithms for surmounting challenges of complexity, robustness, and scalability. This report summarizes the outcomes of that workshop, "Randomized Algorithms for Scientific Computing (RASC)," held virtually across four days in December 2020 and January 2021.

We consider solving nonlinear optimization problems with a stochastic objective and deterministic equality constraints. We assume for the objective that its evaluation, gradient, and Hessian are inaccessible, while one can compute their stochastic estimates by, for example, subsampling. We propose a stochastic algorithm based on sequential quadratic programming (SQP) that uses a differentiable exact augmented Lagrangian as the merit function. To motivate our algorithm design, we first revisit and simplify an old SQP method \citep{Lucidi1990Recursive} developed for solving deterministic problems, which serves as the skeleton of our stochastic algorithm. Based on the simplified deterministic algorithm, we then propose a non-adaptive SQP for dealing with stochastic objective, where the gradient and Hessian are replaced by stochastic estimates but the stepsizes are deterministic and prespecified. Finally, we incorporate a recent stochastic line search procedure \citep{Paquette2020Stochastic} into the non-adaptive stochastic SQP to adaptively select the random stepsizes, which leads to an adaptive stochastic SQP. The global "almost sure" convergence for both non-adaptive and adaptive SQP methods is established. Numerical experiments on nonlinear problems in CUTEst test set demonstrate the superiority of the adaptive algorithm., Comment: 60 pages, 24 figures

For optimal power flow problems with chance constraints, a particularly effective method is based on a fixed point iteration applied to a sequence of deterministic power flow problems. However, a priori, the convergence of such an approach is not necessarily guaranteed. This article analyses the convergence conditions for this fixed point approach, and reports numerical experiments including for large IEEE networks.

We study solution sensitivity for nonlinear programs (NLPs) whose structures are induced by graphs. These NLPs arise in many applications such as dynamic optimization, stochastic optimization, optimization with partial differential equations, and network optimization. We show that for a given pair of nodes, the sensitivity of the primal-dual solution at one node against a data perturbation at the other node decays exponentially with respect to the distance between these two nodes on the graph. In other words, the solution sensitivity decays as one moves away from the perturbation point. This result, which we call exponential decay of sensitivity, holds under the strong second-order sufficiency condition and the linear independence constraint qualification. We also present conditions under which the decay rate remains uniformly bounded; this allows us to characterize the sensitivity behavior of NLPs defined over subgraphs of infinite graphs. The theoretical developments are illustrated with numerical examples.

Despite cascading failures being the central cause of blackouts in power transmission systems, existing operational and planning decisions are made largely by ignoring their underlying cascade potential. This paper posits a reliability-aware AC Optimal Power Flow formulation that seeks to design a dispatch point which has a low operator-specified likelihood of triggering a cascade starting from any single component outage. By exploiting a recently developed analytical model of the probability of component failure, our Failure Probability-constrained ACOPF (FP-ACOPF) utilizes the system's expected first failure time as a smoothly tunable and interpretable signature of cascade risk. We use techniques from bilevel optimization and numerical linear algebra to efficiently formulate and solve the FP-ACOPF using off-the-shelf solvers. Extensive simulations on the IEEE 118-bus case show that, when compared to the unconstrained and N-1 security-constrained ACOPF, our probability-constrained dispatch points can significantly lower the probabilities of long severe cascades and of large demand losses, while incurring only minor increases in total generation costs., Comment: 15 pages, 6 figures

Many physical datasets are generated by collections of instruments that make measurements at regular time intervals. For such regular monitoring data, we extend the framework of half-spectral covariance functions to the case of nonstationarity in space and time and demonstrate that this method provides a natural and tractable way to incorporate complex behaviors into a covariance model. Further, we use this method with fully time-domain computations to obtain bona fide maximum likelihood estimators---as opposed to using Whittle-type likelihood approximations, for example---that can still be computed efficiently. We apply this method to very high-frequency Doppler LIDAR vertical wind velocity measurements, demonstrating that the model can expressively capture the extreme nonstationarity of dynamics above and below the atmospheric boundary layer and, more importantly, the interaction of the process dynamics across it.

We develop a novel computational method for evaluating the extreme excursion probabilities arising from random initialization of nonlinear dynamical systems. The method uses excursion probability theory to formulate a sequence of Bayesian inverse problems that, when solved, yields the biasing distribution. Solving multiple Bayesian inverse problems can be expensive; more so in higher dimensions. To alleviate the computational cost, we build machine-learning-based surrogates to solve the Bayesian inverse problems that give rise to the biasing distribution. This biasing distribution can then be used in an importance sampling procedure to estimate the extreme excursion probabilities., Comment: arXiv admin note: text overlap with arXiv:2001.11904

We study the convergence properties of an overlapping Schwarz decomposition algorithm for solving nonlinear optimal control problems (OCPs). The algorithm decomposes the time domain into a set of overlapping subdomains, and solves all subproblems defined over subdomains in parallel. The convergence is attained by updating primal-dual information at the boundaries of overlapping subdomains. We show that the algorithm exhibits local linear convergence, and that the convergence rate improves exponentially with the overlap size. We also establish global convergence results for a general quadratic programming, which enables the application of the Schwarz scheme inside second-order optimization algorithms (e.g., sequential quadratic programming). The theoretical foundation of our convergence analysis is a sensitivity result of nonlinear OCPs, which we call "exponential decay of sensitivity" (EDS). Intuitively, EDS states that the impact of perturbations at domain boundaries (i.e. initial and terminal time) on the solution decays exponentially as one moves into the domain. Here, we expand a previous analysis available in the literature by showing that EDS holds for both primal and dual solutions of nonlinear OCPs, under uniform second-order sufficient condition, controllability condition, and boundedness condition. We conduct experiments with a quadrotor motion planning problem and a PDE control problem to validate our theory; and show that the approach is significantly more efficient than ADMM and as efficient as the centralized solver Ipopt., Comment: 16 pages

In this paper we focus on the parameter estimation of dynamic load models with stochastic terms, in particular, load models where protection settings are uncertain, such as in aggregated air conditioning units. We show how the uncertainty in the aggregated protection characteristics can be formulated as a stochastic differential equation with process noise. We cast the parameter inversion within a Bayesian parameter estimation framework, and we present methods to include process noise. We demonstrate the benefits of considering stochasticity in the parameter estimation and the risks of ignoring it., Comment: 7 pages, 9 figures, to be published in Electric Power Systems Research PSCC 2020

We present an overlapping Schwarz decomposition algorithm for constrained quadratic programs (QPs). Schwarz algorithms have been traditionally used to solve linear algebra systems arising from partial differential equations, but we have recently shown that they are also effective at solving structured optimization problems. In the proposed scheme, we consider QPs whose algebraic structure can be represented by graphs. The graph domain is partitioned into overlapping subdomains (yielding a set of coupled subproblems), solutions for the subproblems are computed in parallel, and convergence is enforced by updating primal-dual information in the overlapping regions. We show that convergence is guaranteed if the overlap is sufficiently large and that the convergence rate improves exponentially with the size of the overlap. Convergence results rely on a key property of graph-structured problems that is known as exponential decay of sensitivity. Here, we establish conditions under which this property holds for constrained QPs (as those found in network optimization and optimal control), thus extending existing work that addresses unconstrained QPs. The numerical behavior of the Schwarz scheme is demonstrated by using a DC optimal power flow problem defined over a network with 9,241 nodes.

We develop a novel computational method for evaluating the extreme excursion probabilities arising for random initialization of nonlinear dynamical systems. The method uses a Markov chain Monte Carlo or a Laplace approximation approach to construct a biasing distribution that in turn is used in an importance sampling procedure to estimate the extreme excursion probabilities. The prior and likelihood of the biasing distribution are obtained by using Rice's formula from excursion probability theory. We use Gaussian mixture biasing distributions and approximate the non-Gaussian initial excitation by the method of moments to circumvent the linearity and Gaussianity assumptions needed by excursion probability theory. We demonstrate the effectiveness of this computational framework for nonlinear dynamical systems of up to 100 dimensions.

We develop a one-Newton-step-per-horizon, online, lag-$L$, model predictive control (MPC) algorithm for solving discrete-time, equality-constrained, nonlinear dynamic programs. Based on recent sensitivity analysis results for the target problems class, we prove that the approach exhibits a behavior that we call superconvergence; that is, the tracking error with respect to the full horizon solution is not only stable for successive horizon shifts, but also decreases with increasing shift order to a minimum value that decays exponentially in the length of the receding horizon. The key analytical step is the decomposition of the one-step error recursion of our algorithm into algorithmic error and perturbation error. We show that the perturbation error decays exponentially with the lag between two consecutive receding horizons, while~the algorithmic error, determined by Newton's method, achieves quadratic convergence instead. Overall this approach induces our local exponential convergence result in terms of the receding horizon length for suitable values of $L$. Numerical experiments validate our theoretical findings., Comment: 21 pages

In this work, cascading transmission line failures are studied through a dynamical model of the power system operating under fixed conditions. The power grid is modeled as a stochastic dynamical system where first-principles electromechanical dynamics are excited by small Gaussian disturbances in demand and generation around a specified operating point. In this context, a single line failure is interpreted in a large deviation context as a first escape event across a surface in phase space defined by line security constraints. The resulting system of stochastic differential equations admits a transverse decomposition of the drift, which leads to considerable simplification in evaluating the quasipotential (rate function) and, consequently, computation of exit rates. Tractable expressions for the rate of transmission line failure in a restricted network are derived from large deviation theory arguments and validated against numerical simulations. Extensions to realistic settings are considered, and individual line failure models are aggregated into a Markov model of cascading failure inspired by chemical kinetics. Cascades are generated by traversing a graph composed of weighted edges representing transitions to degraded network topologies. Numerical results indicate that the Markov model can produce cascades with qualitative power-law properties similar to those observed in empirical cascades.

In this paper, we study the sensitivity of discrete-time dynamic programs with nonlinear dynamics and objective to perturbations in the initial conditions and reference parameters. Under uniform controllability and boundedness assumptions for the problem data, we prove that the directional derivative of the optimal state and control at time $k$, $\boldsymbol{x}^*_{k}$ and $\boldsymbol{u}^*_{k}$, with respect to the reference signal at time $i$, $\boldsymbol{d}_{i}$, will have exponential decay in terms of $|k-i|$ with a decay rate $\rho$ independent of the temporal horizon length. The key technical step is to prove that a version of the convexification approach proposed by Verschueren et al. can be applied to the KKT conditions and results in a convex quadratic program with uniformly bounded data. In turn, Riccati techniques can be further employed to obtain the sensitivity result, borne from the observation that the directional derivatives are solutions of quadratic programs with structure similar to the KKT conditions themselves. We validate our findings with numerical experiments on a small nonlinear, nonconvex, dynamic program., Comment: 28 pages