Paper deep dive
Tools to Explain Neural Networks for Power System Dynamics
Petros Ellinas, Johanna Vorwerk, Spyros Chatzivasileiadis
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%
Last extracted: 8/11/2026, 5:03:24 AM
Summary
This paper introduces analytical tools to explain the training performance of machine learning surrogate models for power system dynamics. It utilizes the Neural Tangent Kernel (NTK) method to provide a modal interpretation of learning performance, linking physical stiffness in power system Differential-Algebraic Equations (DAEs) to optimization stiffness during Neural Network (NN) training. The authors develop adaptive loss-weighting strategies and demonstrate that structure-aware architectures like ActNet outperform vanilla NNs in handling stiff, multi-timescale dynamics. The approach is validated on physics-informed ML surrogates for Synchronous Machines and power electronic converters.
Entities (7)
Relation Signals (5)
Neural Tangent Kernel → usedforexplaining → training performance
confidence 95% · this paper presents... analytical tools to explain the training performance of machine learning surrogate models... uses the Neural Tangent Kernel (NTK) method.
Physics-Informed Machine Learning → incorporates → physical structure
confidence 93% · Physics-Informed Machine Learning (PIML) approaches have emerged that improve reliability and data efficiency by incorporating the physical structure directly into training
ActNet → outperforms → vanilla NNs
confidence 92% · structure-aware neural architectures, such as ActNet, perform better than vanilla NNs.
Physical Stiffness → manifestsas → optimization stiffness
confidence 90% · This connection explains how physical stiffness and timescale separation in power system dynamic models appear as optimization stiffness during Neural Network (NN) training.
Adaptive Loss-Weighting → improves → Convergence
confidence 88% · guides the use of adaptive loss-weighting strategies that rebalance modes and improve convergence in mildly stiff regimes.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:This paper presents, for the first time in power systems literature to our knowledge, analytical tools to explain the training performance of machine learning surrogate models for power system dynamics. Power system simulations are increasingly challenged by stiff and multi-timescale dynamics arising from converter-interfaced resources and fast control loops. Machine learning surrogates emerge as promising tools to handle this complexity and accelerate dynamic simulations. However, their performance remains difficult to interpret, which limits their adoption. Building on the small-signal eigenvalue analysis in power systems, this paper uses the Neural Tangent Kernel (NTK) method. NTK delivers a modal interpretation of the learning performance, identifying error modes that decay rapidly versus others that converge slowly. This connection explains how physical stiffness and timescale separation in power system dynamic models appear as optimization stiffness during Neural Network (NN) training. Based on this analysis, we develop adaptive loss-weighting strategies to improve and explain why structure-aware neural architectures, such as ActNet, perform better than vanilla NNs. We assess the proposed approach on physics-informed machine learning surrogate models of \acp{SM} and power electronic converters. The methods introduced in this paper can deliver the necessary analytical tools to interpret and improve the performance of machine learning surrogates, paving the way for the systematic, physics-aware design of NN architectures and training strategies. By moving beyond trial-and-error development, these tools reveal training dynamics and failure modes, support more reliable design decisions, and strengthen confidence in machine-learning surrogates for engineering applications.
Tags
Links
- Source: https://arxiv.org/abs/2608.08048v1
- Canonical: https://arxiv.org/abs/2608.08048v1
Trouble viewing inline? Open PDF directly →
Full Text
60,350 characters extracted from source content.
Expand or collapse full text
Tools to Explain Neural Networks for Power System Dynamics Petros Ellinas, Johanna Vorwerk, and Spyros Chatzivasileiadis This work is supported by the European Research Council (ERC) Starting Grant VeriPhIED, Grant Agreement No. 949899. Abstract This paper presents, for the first time in power systems literature to our knowledge, analytical tools to explain the training performance of machine learning surrogate models for power system dynamics. Power system simulations are increasingly challenged by stiff and multi-timescale dynamics arising from converter-interfaced resources and fast control loops. Machine learning surrogates emerge as promising tools to handle this complexity and accelerate dynamic simulations. However, their performance remains difficult to interpret, which limits their adoption. Building on the small-signal eigenvalue analysis in power systems, this paper uses the Neural Tangent Kernel (NTK) method. NTK delivers a modal interpretation of the learning performance, identifying error modes that decay rapidly versus others that converge slowly. This connection explains how physical stiffness and timescale separation in power system dynamic models appear as optimization stiffness during Neural Network (N) training. Based on this analysis, we develop adaptive loss-weighting strategies to improve and explain why structure-aware neural architectures, such as ActNet, perform better than vanilla NNs. We assess the proposed approach on physics-informed machine learning surrogate models of Synchronous Machines and power electronic converters. The methods introduced in this paper can deliver the necessary analytical tools to interpret and improve the performance of machine learning surrogates, paving the way for the systematic, physics-aware design of N architectures and training strategies. By moving beyond trial-and-error development, these tools reveal training dynamics and failure modes, support more reliable design decisions, and strengthen confidence in machine-learning surrogates for engineering applications. Index Terms: Neural Tangent Kernel, Physics-Informed Neural Networks, Scientific Machine Learning, Power System Dynamics, Stiff and Multi-Timescale Systems I Introduction The rapid growth of renewable generation and the electrification of heating and transport have increased the need for fast and reliable dynamic analysis of large-scale nonlinear power systems. Operators must assess system security across many contingencies and operating points. Each time-domain simulation, however, is computationally demanding because the governing Differential-Algebraic Equations (DAEs) are nonlinear and multi-timescale: fast electromagnetic and control modes coexist with slow electromechanical modes. While all have been present in power systems dominated by SMs, converter-interfaced generation convolutes timescale separation and requires considering multiscale DAEs in dynamic simulations [7, 10]. This need to consider a wide variety of time scales enhances stiffness, which forces explicit solvers to take very small time steps and causes implicit solvers to face increasingly ill-conditioned algebraic systems. When these simulations must be repeated thousands of times to cover the wide range of credible scenarios, classical methods quickly become too slow for real-time or near-real-time decision support. Many acceleration strategies have been proposed, including model reduction, network equivalencing, high-performance computing implementations, and specialized techniques to accelerate electromagnetic transient (EMT) simulations [14]. More recently, Machine Learning (ML) approaches have been explored to approximate the dynamics of specific power-system components, either by directly approximating system trajectories [13, 2] or by learning surrogate models of the underlying dynamics [1]. These efforts broadly follow two complementary directions: (i) stand-alone surrogates that directly map operating conditions and time to system trajectories for rapid open-loop screening [13, 2], and (i) learned surrogates embedded within a numerical DAE solver [1]. Both approaches exploit Neural Network (N) interpolation and generalization capabilities across operating conditions, enabling extremely fast evaluation after training. Along this line of research, Physics-Informed Machine Learning (PIML) approaches have emerged that improve reliability and data efficiency by incorporating the physical structure directly into training and penalizing violations of the governing DAEs and initial conditions at collocation points [11]. The ML model maps initial states, control inputs, and time to system trajectories, while minimizing the DAE residual to constrain the surrogate to the physically admissible solution manifold. When used in an open-loop setting, e.g., in a single-machine infinite bus, this enables the evaluation of several time points in parallel, rather than sequential time stepping, which drastically increases computation speed. Here, it is interesting to highlight the analogy between physics-based power-system dynamic models formulated as DAEs and surrogate models based on NNs. In the ML-based surrogate modeling of dynamical systems, N training determines the learnable parameters, i.e. weights and biases, by solving an optimization problem so that the model reproduces the time evolution of the states exactly as described by a physics-based system of DAEs. In this sense, the role of the explicit time integration of the DAE system is replaced by an optimization-based training procedure. As a result, classical notions of numerical stability of the solver reappear as questions of convergence, conditioning, and robustness of the underlying optimization process. Despite their advantages, ML surrogates are still treated with caution in power-system applications because their internal behavior is often opaque, and it is unclear how design choices and training objectives affect accuracy and reliability. In safety-critical contexts, this opacity is also increasingly at odds with emerging expectations for transparency and auditability [12]. In this work, we view transparency not only as the interpretability of model internals but also as the ability to predict and explain reliability and failure modes, as we usually do from the structure of the underlying DAEs. This motivates a structural question: How does physical stiffness in the governing DAEs affect the numerical stability and reliability of the N training optimization problem, and how can its failure modes be predicted? In this paper, we study ML surrogates for power-system components and address four central questions: (i) How do different neural architectures learn stiff and multi-scale dynamics? (i) How do modeling choices affect accuracy and reliability? (i) (i) how reliability and failure modes, meaning when and why training becomes inaccurate or fails to converge, can be predicted and explained? and (iv) how these choices influence generalization beyond the training conditions? We treat the surrogate architecture, the choice of activation functions, and the design of the loss function as deliberate modeling choices meant to adapt to the mathematical structure of the governing DAEs, rather than heuristic design decisions. To answer these questions, we focus on individual components, both because much of the recent literature studies component-level ML surrogates [2] and because low-order models keep the NTK analysis interpretable, allowing the underlying mechanisms and failure modes to be isolated. The proposed analysis is not restricted to low-order systems and can also be applied to ML surrogates of higher-order dynamic models. To study how N architectures interact with stiffness and the multi-timescale behavior of power system dynamics, we use the Neural Tangent Kernel (NTK) method [5] to analyze the training dynamics. More generally, NTK provides a tool for analyzing the convergence of individual objectives in multi-objective N training, where different loss components may drive parameter updates with different strengths and converge at different rates. In the PIML surrogate setting considered here, these loss components correspond to the initial-condition and physics-residual terms. Analogous to eigenvalues in small-signal stability analysis, the NTK provides a spectral decomposition of training dynamics: NTK eigenvalues determine how different loss-function error modes decay, just as Jacobian eigenvalues characterize the post-disturbance evolution and decay of physical modes. This connects physical stiffness in the DAEs with optimization stiffness in training and guides the use of adaptive loss-weighting strategies that rebalance modes and improve convergence in mildly stiff regimes. The main contributions of this work include: • A diagnostic tool for multi-objective learning. Using the NTK method, we propose a comprehensive analysis tool for multi-objective training, including PIML models. We present an analogy to small-signal stability analysis to provide an intuitive interpretation of the training dynamics. It identifies when training is well-posed or ill-conditioned by explaining how different loss components and embedded physical properties affect learnability and convergence. Furthermore, the NTK analysis guides the use of adaptive loss weighting to promote balanced convergence across loss terms. • We compare N architectures for stiff and multi-timescale dynamics, including ActNet. Based on observed performance and our NTK analysis, we provide practical guidelines for selecting architectures and activation functions according to the structural properties of the approximated dynamics. • We evaluate the proposed methods on SM and inverter models (ranging from 2nd order to 11th order models), highlighting both their strengths and limitations across dynamical regimes. Sections I–V present the DAE formulation, the PIML surrogate setup, the structure-aware neural architectures, the NTK training analysis, and the case-study results; the paper closes with conclusions on trustworthy surrogate modeling for power-system dynamics. I Differential–Algebraic Systems for Power Systems DAEs model power-system dynamics, which combine differential equations with algebraic constraints: d(t)dt−f((t),(t),t)=0,(t0)=0.M\, dx(t)dt-f (x(t),u(t),t )=0, (t_0)=x_0. (1) Here (t)∈ℝnx(t) ^n collects both dynamic states, e.g., generator rotor angles, speeds, and controller states, and algebraic variables, e.g., bus voltages. The input (t)u(t) represents disturbances or control actions. The mass matrix M contains ones on differential rows and zeros on algebraic rows. Rows with Mii=0M_i=0 impose instantaneous algebraic constraints, such as network power balance, which are enforced at each time step using a root-finding method, typically Newton–Raphson. To simulate the DAEs, time-stepping schemes advance the solution either explicitly or implicitly. Explicit methods compute k+1x_k+1 directly from known quantities and are inexpensive per step, but require small time steps in the presence of fast dynamics. Implicit methods allow larger time steps and are therefore common in transient-stability programs, but require solving a nonlinear system F(k,k+1)=0F(x_k,x_k+1)=0 [4]. Power-system models are stiff because fast electrical or control states and slower electromechanical states are coupled within the same DAE system, even though their time scales may differ by several orders of magnitude. As a result, explicit solvers must use time steps small enough to remain stable for the fastest dynamics, even when the main quantities of interest evolve much more slowly. Implicit solvers alleviate this restriction, but at the cost of solving more challenging nonlinear systems. A more detailed discussion of stiffness and numerical-solver stability is provided in the companion document [3]. These challenges motivate surrogate approaches that approximate system evolution without a full nonlinear solve at every step. I PIML for Differential–Algebraic Systems This section introduces the neural solution operator, its physics-informed training objective, and the structure-aware architectures considered in this work. I-A Neural Networks NNs are a fundamental class of ML models composed of L layers of interconnected neurons. Formally, a N defines a parametric mapping N:ℝn0→ℝnLN_ w:R^n_0 ^n_L, which transforms an input vector ∈ℝn0z ^n_0 into an output ∈ℝnLy ^n_L through a sequence of layer-wise transformations l _l: N()=L−1(L−2(…0()))N_ w(z)= _L-1\! ( _L-2\! (… _0(z) ) ) (2) Each layer applies an affine transformation followed by a nonlinear activation. Because the training data only constrain the model at sampled points, the architecture affects how the N behaves between those points and which solution it learns during training. Concretely, each layer l _l is defined as l(l−1)=ϕl(ll−1+l), _l(z_l-1)= _l\! ( W_l\,z_l-1+b_l ), where l W_l is the weight matrix, lb_l is the bias vector, l−1z_l-1 is the input to layer l, and ϕl _l is a nonlinear basis function, usually called an activation function. By stacking such layers, NNs construct hierarchical representations and efficiently approximate complex functions. For compact notation, we collect all trainable weights and biases in :=l,ll=0L−1 w:=\ W_l,b_l\_l=0^L-1. Ultimately, training a N consists of determining, iteratively, the parameters w that minimize a task-dependent objective or loss function ℒ(N),L(N_ w), thereby selecting a function from the set of functions representable by the chosen architecture. In the context of DAE systems, this objective must reflect the governing physical constraints. I-B Physics-Informed Machine Learning Building on the N representation above, we now formulate a physics-informed solution operator for dynamic power-system components. The solution operator, also called the flow map, maps a consistent initial condition and prescribed disturbances or control inputs to the system state at a future time. To approximate the solution operator (0,(⋅),t0,t)↦(t0+t)(x_0,u(·),t_0,t) (t_0+t) without explicit time stepping at inference, we use a N: N(0,(⋅),t0,t)≈(t0+t)N_ w(x_0,u(·),t_0,t) (t_0+t) (3) Rather than relying on labeled trajectories, we train NN_ w by enforcing the DAE residual and the initial condition as soft constraints. Let tii=1Nt\t_i\_i=1^N_t denote the set of NtN_t temporal collocation points, and let 0,j,j(⋅)j=1Nic\x_0,j,u_j(·)\_j=1^N_ic denote the set of NicN_ic training scenarios, where 0,jx_0,j is the initial state and j(⋅)u_j(·) is the corresponding input trajectory of scenario j. The network parameters w are obtained by minimizing ℒ =ρdae1NtNic∑i=1Nt∑j=1Nic‖∂tNNi,j−fi,j‖22 = _dae 1N_tN_ic _i=1^N_t _j=1^N_ic \|M _tN_i,j-f_i,j \|_2^2 +ρic1Nic∑j=1Nic‖N0,j−0,j‖22 + _ic 1N_ic _j=1^N_ic \|N_0,j-x_0,j \|_2^2 (4) where NNi,j _i,j :=N(0,j,j(⋅),0,ti), :=N_ w (x_0,j,u_j(·),0,t_i ), N0,j _0,j :=N(0,j,j(⋅),0,0), :=N_ w (x_0,j,u_j(·),0,0 ), fi,j f_i,j :=f(NNi,j,j(ti),ti), :=f (N_i,j,u_j(t_i),t_i ), ∂tNNi,j _tN_i,j :=∂tN(0,j,j(⋅),0,t)|t=ti := . ∂ tN_ w (x_0,j,u_j(·),0,t ) |_t=t_i The coefficients ρdae _dae and ρic _ic weight the DAE-residual and initial-condition losses, respectively, while ∥⋅∥2\|·\|_2 denotes the Euclidean norm. Since N(0,(⋅),t0,t)N_ w(x_0,u(·),t_0,t) depends on the continuous variables 0x_0, (⋅)u(·), and t, the problem formulation becomes ∂N(0,(⋅),t0,t)∂t−f(N(0,(⋅),t0,t),(t),t)=0,M _ w(x_0,u(·),t_0,t)∂ t-f\! (N_ w(x_0,u(·),t_0,t),u(t),t )=0, (5) which defines a partial differential equation over the full input domain. Training therefore requires representative sampling of (0,,t0,t)(x_0,u,t_0,t) so that the learned operator produces consistent trajectories across different initial conditions and inputs. When trajectory-simulated samples are available, a supervised data loss can be added to anchor the solution and ease the physics-informed optimization. In this work, we restrict attention to operators that take a single instantaneous input value (t)u(t) and set t0=0t_0=0 without loss of generality. Although the formulation is developed for general index-1 DAE systems, PIML methods have predominantly been applied to Ordinary Differential Equations (ODEs), since DAEs introduce algebraic constraints that tightly couple state variables and complicate both training and numerical stability. Geometric Interpretation of Integration General DAE systems define a flow map whose trajectories lie on a manifold. Time-stepping algorithms approximate the continuous evolution, starting from an initial condition, by moving along this manifold through a sequence of linear steps, according to the local dynamics at the current point, as shown in Fig. 1. images/flowmap Figure 1: Flowmap view of a continuous DAE. The surface represents the flowmap manifold. The curve shows the system’s continuous trajectory. Discrete solver steps move along this curve, approximating the underlying flow map. In contrast, the PIML approaches aim to approximate the flowmap itself as a global function, replacing step-by-step advancement with direct evaluation of the learned operator. Note that these PIML operators are often understood to be reduced order models. However, they do not project onto a lower-dimensional state. Instead, they can be understood as noisy replicas of the true system. I-C Embedding Physical Knowledge into ML Models Physics-informed losses enforce the governing equations during training, but they do not, on their own, determine how the ML models interpolate between the training points to represent the underlying dynamics; that also depends on the architecture. In power-system applications, where oscillatory and multi-timescale behavior is often known a priori, it is therefore natural to embed such structure directly into the model design [16]. A representative example is the Kolmogorov-Arnold Neural Network (KAN) [2, 9, 8], which replaces fixed activation functions such as tanh with learnable univariate spline functions. This makes the model more adaptable during training and can help it capture the geometry of the flow-map manifold [2]. However, spline evaluation and basis expansion make KANs significantly more expensive at inference time. In power-system studies, where a surrogate may need to be evaluated repeatedly for many operating points or contingencies, this computational burden can become limiting. The choice of activation function is central to this trade-off. Smooth monotone activations such as tanh train stably, whereas sinusoidal activations better match oscillatory behavior at the cost of more sensitive optimization. Although other periodic activations, e.g. [18], may also be useful for power-system dynamics, their systematic evaluation is outside the scope of this work. These considerations are directly relevant in power systems, where rotor-angle states are periodic and many transient responses contain oscillatory electromechanical or fast controller-driven modes. As a result, the activation choice affects not only training stability but also whether the surrogate reliably captures the system’s dominant physical behavior. These considerations motivate the use of architectures that retain the expressive benefits of learnable activation functions while remaining computationally lightweight. I-C1 ActNet Layers [width=0.7]images/actnet2 Figure 2: ActNet layer structure. Each input variable is passed through trainable sinusoidal channels and mixed through coefficients cijc_ij, allowing the layer to represent timescales within a feedforward architecture. Figure 2 illustrates the main idea of ActNet. Each input is processed by several trainable sinusoidal channels, which can represent different oscillation frequencies. The resulting channel responses are weighted and combined to form the layer output, allowing slow and fast oscillatory components to be represented simultaneously. Besides KANs, ActNet provides a lightweight alternative that embeds oscillatory structure in its activation functions. It represents the solution using trainable sinusoidal channels that can capture slow and fast oscillatory components within the same architecture. Here, a channel denotes one trainable sinusoidal basis function. For each channel i=1,…,mi=1,…,m, ActNet defines the univariate adaptive warp ϕi(t)=sin(ωit+pi)−μiσi+ε, _i(t)= ( _it+p_i)- _i _i+ , (6) where ωi _i and pip_i are trainable frequency and phase parameters, and μi,σi _i, _i normalize the output for numerical stability. Since ϕi _i depends on trainable frequency and phase parameters, ActNet can be interpreted as using learnable sinusoidal activation functions. Given an input vector =(x1,…,xn)x=(x_1,…,x_n), ActNet forms the mixed response of channel i as si()=∑j=1ncijϕi(xj),s_i(x)= _j=1^nc_ij\, _i(x_j), (7) and aggregates the channel responses as N()=∑i=1msi().N_ w(x)= _i=1^ms_i(x). (8) Here, cijc_ij determines how strongly input variable xjx_j contributes to the i-th sinusoidal channel. The parameter ωi _i determines the frequency of the i-th sinusoidal channel, while cijc_ij weights the contribution of input xjx_j to that channel. For ActNet, :=ωi,pii=1m∪ciji=1,j=1m,n. w:=\ _i,p_i\_i=1^m∪\c_ij\_i=1,j=1^m,n. This structure aligns the architecture with the multi-timescale nature of power-system dynamics. Channels with smaller ωi _i can represent slowly varying components, while channels with larger ωi _i can represent faster oscillatory components. This is useful for component models in which electromechanical, electrical, and fast-control interactions shape the dynamic response. Compared to spline-based KAN architectures, ActNet remains computationally lighter while retaining a structure-aware representation for stiff and multi-timescale dynamics. The resulting layer structure is shown in Fig. 2. IV Neural Tangent Kernel and Training Dynamics The input–output mapping of a general N is highly nonlinear and nonconvex in its parameters. Therefore, understanding how the trainable parameters w evolve and why training succeeds or fails requires tools beyond classical approximation theory. In this work, we adopt the NTK as a principled framework for analyzing how NNs adapt their parameters during training. The NTK for NNs is analogous to small-signal stability analysis of dynamic models. When a nonlinear dynamical system ˙=F() x=F(x) is linearized around an operating point ∗x^*, the resulting state-space matrix =∂F/∂|∗A= .∂ F/ |_x^* describes the local behavior of small perturbations around that operating point. Its eigenvalues determine whether modes decay, oscillate, or grow, and at what rates. This modal interpretation provides the conceptual basis for comparing physical-system modes with training-error modes. NTK theory provides a similar modal view of N training. During training, gradient descent updates the trainable parameters w, which changes the N output and reduces the training error. The NTK describes how sensitive the N outputs at different training samples are to the parameter updates of w, showing whether errors at different inputs or loss function terms reduce together or evolve independently during gradient descent. Under the infinite-width NTK approximation, the N can be linearized around its initialization and the NTK is treated as approximately fixed during training. In this regime, the NTK plays the role of a training-domain state-space matrix: its eigenvectors define independent training-error, i.e, loss function, modes, and its eigenvalues determine their decay rates. This allows us to interpret N training using the same modal intuition used for power-system small-signal dynamics. During training by gradient descent on a loss ℒ(N)L(N_ w), the parameters follow the continuous gradient-flow dynamics ddτ=−∇ℒ(N), d wdτ=-\, _ wL(N_ w), (9) where ∇ℒ _ wL denotes the gradient of the loss with respect to the trainable parameter vector w, and τ denotes continuous training time, i.e., the continuous-time analogue of training iterations or epochs. Consider the mean-squared error loss ℒ(N)=12∑i=1N(N(i)−i)2,L(N_ w)= 12 _i=1^N (N_ w(x_i)-y_i )^2, where iy_i denotes the target output at the training point ix_i. The induced learning dynamics of the network output satisfy dτN()=−∑i=1N ddτN_ w(x)=- _i=1^N (N(i)−i) (N_ w(x_i)-y_i ) (10) ×⟨∇N(i),∇N()⟩. × _ wN_ w(x_i), _ wN_ w(x) . The inner product (i,j)=⟨∇N(i),∇N(j)⟩ (x_i,x_j)= _ wN_ w(x_i), _ wN_ w(x_j) (11) is the NTK. Collecting the predictions on the training set into the vector (τ)=(N(τ)(1),…,N(τ)(N))⊤,U(τ)= (N_ w(τ)(x_1),…,N_ w(τ)(x_N) ) , the training error is (τ)=(τ)−e(τ)=U(τ)-y. Each eigenvector of the NTK describes a distinct way in which prediction errors are distributed across the training samples, and the total training error can be understood as a weighted sum of these error modes. The weights of this weighted sum indicate how strongly each mode contributes to the current error, while the eigenvalue aia_i associated with each mode i determines how quickly that contribution is reduced during training: large eigenvalues indicate fast learning, whereas small eigenvalues indicate slowly corrected or nearly unchanged errors. This is analogous to expressing a physical disturbance as a combination of dynamic modes in small-signal analysis. In these modal coordinates, the training dynamics decouple into independent scalar equations daidτ=−λiai da_idτ=- _i\,a_i (12) so each error component evolves independently as ai(τ)=ai(0)e−λiτa_i(τ)=a_i(0)e^- _iτ. Thus, the eigenvectors define independent error directions, while the eigenvalues λi _i directly determine the decay rate of each training-loss error mode: large eigenvalues lead to fast convergence, and small eigenvalues produce slow decay of the corresponding residual errors. Under the infinite-width NTK approximation, the NTK is treated as approximately constant during training, (i,j)≈0(i,j) (x_i,x_j)≈ _0(x_i,x_j) [5, 17]. This means that, although the N output changes, the sensitivity of the output to the trainable parameters changes only weakly. Under the fixed-NTK approximation, the output dynamics reduce to: ddτ=−0(−), dUdτ=-\, _0(U-y), (13) so each training mode decays at a rate determined by the eigenvalues of the initial NTK. Please note that this is an approximation. In particular, when the parameters change substantially throughout training, the NTK values may also vary. []images/SmallSignalStability []images/NTKAnalogy Figure 3: Analogy between small-signal stability analysis and NTK training dynamics. Top: eigenvalues of the state-space matrix determine the behavior of physical perturbation modes. Bottom: NTK eigenvalues determine the decay rates of training-error modes. Large NTK eigenvalues correspond to fast error decay, while small eigenvalues correspond to slow or nearly frozen modes. Fig. 3 summarizes the resulting parallel between power-system dynamic modes and training-error modes. In small-signal analysis, the eigenvalues of the linearized state-space matrix determine how dynamic modes behave after a system perturbation. The location of an eigenmode determines whether an oscillation occurs, at what rate, and whether it decays or grows. Similarly, the eigenvalues of the NTK determine how different training-error directions over the sampled points decay during the training process. An uneven NTK spectrum indicates optimization stiffness: some loss function error modes are corrected rapidly, while others remain nearly frozen. These slow NTK modes can limit the accuracy of the surrogate, especially for fast or oscillatory parts of the physical response. In contrast, an even NTK spectrum indicates that all errors decay at a similar rate and are learned with comparable speed and accuracy. Consequently, an ideal NTK profile would mean that NTK eigenmodes are similar across error components. IV-A NTK for PIML Models PIML training involves multiple loss terms that enforce data, initial conditions, and physics constraints. The coupling between the initial-conditions loss ℒb(ℒ(N(τ)))L_b(L(N_ w(τ))) and physics loss ℒp(N(τ)))L_p(N_ w(τ))) from (I-B) can be analyzed through the block NTK matrix K, so that their interaction is approximated by [ℒ˙b(N(τ)))ℒ˙p(N(τ)))]=−[bbpbpp]⏟[ℒb(N(τ)))ℒp(N(τ)))], bmatrix L_b(N_ w(τ)))\\ L_p(N_ w(τ))) bmatrix=- bmatrixK_b&K_bp\\ K_pb&K_p bmatrix_K bmatrixL_b(N_ w(τ)))\\ L_p(N_ w(τ))) bmatrix, (14) where each block is a Gram matrix of inner products between parameter gradients of the corresponding loss terms. For example, if i indexes samples in the initial conditions loss and j indexes samples in the physics loss: (bp)ij=⟨∂ℒb(N(τ)(xib))∂,∂ℒp(N(τ)(xjp))∂⟩.(K_bp)_ij= _b(N_ w(τ)(x_i^b))∂ w, _p(N_ w(τ)(x_j^p))∂ w . Here, K extends the NTK interpretation to the multi-objective setting. Its diagonal blocks quantify how strongly each loss function term drives w updates, while the off-diagonal blocks quantify the interaction between different loss terms. For each diagonal block, large leading eigenvalues correspond to training-error directions that respond strongly to changes in w and therefore decay rapidly under gradient descent. A long tail of small eigenvalues indicates error directions that respond weakly to parameter updates, so they decay slowly and may limit convergence. Since the trace of the matrix, Tr(Kii)Tr(K_i) equals the sum of the eigenvalues of the i-th block, it measures the overall strength with which the loss term i affects the training update. Large differences in Tr(Kii)Tr(K_i) therefore indicate imbalanced training, where some objectives dominate while others evolve slowly. IV-B NTK Eigenvalues for Adaptive Weights During Training Since the NTK eigenmodes provide insights into how fast each loss component is learned during training, they can be used to update the training weights during gradient descent. Guided by the NTK analysis, we formulate an adaptive weighting scheme that promotes uniform convergence across the individual loss objectives. Each loss term is assigned a weight ρi _i, chosen inversely to the trace of its corresponding self-block: ρi=∑jTr(Kjj)Tr(Kii), _i= _jTr(K_j)Tr(K_i), (15) where KiiK_i is the diagonal block associated with the i-th loss. A large trace indicates that the loss strongly drives the parameter update, while a small trace indicates a weak contribution. The inverse-trace weighting acts as a training preconditioner: dominant objectives are down-weighted and weaker objectives are up-weighted, promoting more balanced convergence across loss terms. V Case Studies In this section, we evaluate vanilla PINNs, KANs, and ActNet on power system time-domain simulation tasks with increasing stiffness and timescale separation. We first report standard accuracy and runtime metrics, and then use the NTK analysis to interpret the results in terms of training-error modes and loss-term imbalance. The goal is to connect physical stiffness in the underlying power-system dynamic model with optimization stiffness in PIML training, and to assess how different architectures represent oscillatory, stiff, and multi-timescale dynamics. We analyze three sets of power system test cases to illustrate three qualitatively different regimes: 1. Mildly Stiff Dynamics (Synchronous Generator Models): We first use classic SM models as baseline test cases to compare vanilla PINNs, KANs, and ActNet in terms of accuracy and runtime in a regime where stiffness is limited. 2. Stiff System with Moderate Time-Scale Separation (Modified Single Machine Infinite Bus (SMIB)): Next, we modify a classic SMIB system to introduce a fast electrical mode and examine whether the architectures can simultaneously learn slow and fast dynamics. The performed NTK analysis provides insights into the performance of the NNs. 3. Strongly Multiscale and Stiff Dynamics (Inverter Model): Finally, we attempt to learn an 11th-order inverter model with nested control loops that induce both stiffness and pronounced timescale separation. Again, we demonstrate how the NTK analysis offers valuable insights into the model performance. Note that to the best of our knowledge, there are no reports in the literature of accurately learning such a stiff DAE system with NNs. The test system, implementation details, training configuration, and resulting performance are analyzed per scenario below, emphasizing how stiffness and multi-timescale structure influence learnability across architectures. Code and supplementary material are available at [3]. V-A Regime 1: Classic Synchronous-Machine Models V-A1 Test Case We first compare the accuracy and runtime of the three architectures on the 2nd (2D), 4th (4D), and 6th-order (6D) SM models formulated in [3]. These models provide a controlled setting in which the dynamical complexity increases with model order. To keep the training problem tractable, the experiments were carried out on a restricted set of initial conditions around the nominal operating region, including in the training 200 simulated trajectories, including an extra term in the loss function [2]. V-A2 Implementation and Training We compare three physics-informed architectures: (i) a vanilla PINN, (i) a KAN, and (i) ActNet. For each SM order (2D, 4D, and 6D), one network is trained for each architecture, resulting in nine trained models in total. The same architecture settings are used across all SM implementations per N type to isolate the effect of model order and architecture type. Specifically, the vanilla PINN uses 3 hidden layers with 64 neurons per layer. ActNet uses the same depth and width as the vanilla PINN, with m=4m=4 trainable sinusoidal channels. The KAN uses a single hidden layer with 20 neurons, providing a compact spline-based alternative. All models are trained with the same physics-informed objective for a fair comparison. Training uses the SOAP optimizer [15] with learning rate 10−410^-4. As numerical baselines, we report runtimes from a CPU-based Runge-Kutta 5 (RK5) solver implemented in Julia and from a GPU-based solver implemented through DiffEqGPU.jl. All GPU-runtime measurements are reported on the same NVIDIA T4 GPU. Each runtime corresponds to a full trajectory calculation over a horizon of T=1sT=1\,s for 1000 scenarios. images/smib_all Figure 4: Accuracy and runtime comparison for the different architectures on the three SM test cases. Runtime is compared for evaluating 1000 trajectories and benchmarked against classic DAE GPU and CPU solvers. V-A3 Results Figure 4 summarizes the MSE and runtime of the three models on the 2D, 4D, and 6D SM test cases, while also detailing the runtimes of both a GPU-based numerical solver and a CPU-based numerical solver. KAN and ActNet reduce the trajectory error relative to the vanilla PINN in the 2D and 4D test cases, with the largest reduction observed in the 4D case. At the same time, all three architectures achieve errors of the same order in the 6D case. Note that ActNet achieves these gains with a parameter count comparable to the vanilla PINN, suggesting that the improvement comes from the design of the activation functions rather than from a simple increase in model size. The main difference between KAN and ActNet appears in inference cost. On the NVIDIA T4 GPU, the vanilla PINN is the fastest model, followed by ActNet, whereas KAN is substantially slower. Across all three test cases, ActNet is about 17×17× faster than KAN. Comparing runtimes, ActNet is slightly faster than the GPU-based solver in the 2D case and up to about 9×9× faster in the 6D case, while being nearly 30×30× faster than the CPU-based solver benchmark. KAN, on the other hand, is slower than the GPU solver. This slowdown of KAN is structural: spline evaluation and basis expansion reduce the degree of vectorization available to KAN, while ActNet remains a fully feedforward architecture. Overall, this baseline comparison shows that architectures with adaptive activation structures, namely KAN and ActNet, improve accuracy over the vanilla PINN. Among them, ActNet provides the best accuracy-runtime trade-off as it offers a computationally lightweight feedforward structure while using learnable sinusoidal activations. V-B Regime 2: Capturing Slow and Fast Dynamics We next introduce stiffness into a classic 2nd-order SMIB model to examine if a single N architecture can capture both slow and fast modes. We therefore augment the 2D2D SMIB model with terminal-voltage dynamics, obtaining a 3D3D system with an additional fast electrical mode and the resulting stiffness. The GPU solver evaluates in approximately 12ms12\,ms for 1000 scenarios, making it a promising candidate for PIML-based speedup. V-B1 The Modified SMIB System To introduce one fast electrical mode into the 2nd-order SMIB dynamics, the terminal voltage V(t)V(t) is introduced as a dynamic state, resulting in the following 3rd-order ODE system: [1000H000τ]dt[δωV]=[ωPm−EVsin(δ)+PL(t)−DωVref−V+ksin(δ)] bmatrix1&0&0\\ 0&H&0\\ 0&0&τ bmatrix ddt bmatrixδ\\ ω\\ V bmatrix= bmatrixω\\ P_m-EV (δ)+P_L(t)-Dω\\ V_ref-V+k (δ) bmatrix Here, δ is the rotor angle relative to the synchronous reference, ω is the speed deviation, and V is the terminal-voltage magnitude; k controls the voltage-oscillation amplitude, while δ sets its phase. Note that all symbol definitions and parameter values are provided in the supplementary material. A small time constant τ introduces a fast electrical mode that creates stiffness in the ODE system and introduces a clear timescale separation: Linearizing around (δ,ω,V)=(0,0,1)(δ,ω,V)=(0,0,1) highlights that the voltage state V contributes the fast electrical eigenvalue, while the rotor states δ and ω yield the slow electromechanical eigenvalue −D/(2H)-D/(2H). From these modal time constants, the corresponding stiffness ratio is: r≈2H/(Dτ)≈2×103r≈ 2H/(Dτ)≈ 2× 10^3, indicating a difference of three orders of magnitude between the fast and slow dynamics. Importantly, this stiffness is essentially one-dimensional, arising from a single dominant fast electrical mode. V-B2 Implementation and Training We train three models: a vanilla PINN with tanh activations, a KAN, and an ActNet, from a single initial condition to isolate architectural and training effects from generalization across initial conditions. All models use the physics-informed loss in (I-B) and are trained for 50005000 epochs using Adam [6]. We observe that the proposed adaptive weighting in (15) accelerates convergence and yields a lower final loss than unweighted training. The vanilla PINN has two hidden layers of 32 neurons, the KAN has one hidden layer of 9 splines with grid size G=7G=7 and spline order k=3k=3, and the ActNet has two hidden layers of 5 neurons with m=2m=2 trainable sinusoidal channels. These architectures were selected through a small grid search and then fixed for a fair comparison. TABLE I: Final loss values for ActNet, KAN, and the vanilla PINN on the modified SMIB test case. Loss Term ActNet KAN PINN LδL_δ !25 4.27×10−54.27× 10^-5 !25 9.91×10−59.91× 10^-5 3.93×10−23.93× 10^-2 LωL_ω !25 1.01×10−51.01× 10^-5 2.95×10−22.95× 10^-2 !252.25×10−12.25× 10^-1 LVL_V 3.89×10−33.89× 10^-3 4.44×10−24.44× 10^-2 !25 4.05×10−14.05× 10^-1 LicL_ic !25 8.54×10−88.54× 10^-8 !253.24×10−53.24× 10^-5 4.49×10−24.49× 10^-2 Ltot L_tot 3.94×10−33.94× 10^-3 7.41×10−27.41× 10^-2 !257.14×10−17.14× 10^-1 images/high_freq_trajectories Figure 5: Predicted trajectories for the modified SMIB test case, comparing the vanilla PINN, KAN, and ActNet. images/eigenvalues_MLP Figure 6: Unweighted NTK eigenvalue spectra of the LδL_δ, LωL_ω, LVL_V, and LicL_ic loss blocks for the three architectures. The horizontal axis gives the eigenvalue rank after sorting by decreasing magnitude. images/eigenvalues_scaled Figure 7: NTK eigenvalue spectra after applying the adaptive weights ρi _i. The horizontal axis gives the eigenvalue rank after sorting by decreasing magnitude. images/eigenvalues_iterations Figure 8: ActNet NTK eigenvalue spectra at selected training checkpoints. Each curve represents the complete spectrum at the indicated iteration, with the horizontal axis denoting the eigenvalue rank. V-B3 Results Figure 5 showcases one time-domain trajectory of the SM frequency for the classic time-domain solver and the three N architectures. While the PINN fails to capture the faster oscillation around the main trajectory, both ActNet and PI-KAN successfully reproduce it. Note that the KAN appears more accurate. These findings are supported by the final physics-loss components across all evaluation trajectories presented in Table I. Here, LδL_δ, LωL_ω, and LVL_V denote the mean-squared residuals of the δ-, ω-, and V-equations, respectively; LicL_ic denotes the initial-condition residual, and LtotL_tot the total loss. ActNet depicts the smallest total loss, more than one order of magnitude below KAN and nearly three below the vanilla PINN, with the largest gain in the residual associated with the oscillatory speed response. We next use the NTK spectra to analyze these differences, which occur despite applying the NTK-weighted training scheme. Generally, the NTK spectrum indicates how different training-error modes respond to parameter updates: large leading eigenvalues correspond to modes that are learned quickly, while small tail eigenvalues correspond to modes that converge slowly. We therefore look for two features: how rapidly the spectrum decays, and whether the dominant eigenvalues are balanced across the residual losses. A sharp decay suggests that only a few modes are learned efficiently. In contrast, a milder decay with larger tail eigenvalues indicates that a broader range of modes remains trainable, including modes associated with faster or more oscillatory response components. Figure 8 depicts the NTK spectra computed on 100 time-collocation points for the residual terms LδL_δ, LωL_ω, and LVL_V at the start of the training. The nearly flat LicL_ic spectrum indicates similar correction rates because the same initial condition is enforced for all collocation samples at t=0t=0, without temporal-derivative terms. The vanilla PINN displays a very heavy tail, with eigenvalues decaying from about 101010^10 to values close to zero. In contrast, KAN and ActNet show a milder decay with larger tail eigenvalues. The dominant eigenvalues are also more balanced across residual terms. In the vanilla PINN, LδL_δ dominates near 101010^10, while LωL_ω and LVL_V remain near 10110^1 to 10210^2. For KAN and ActNet, the largest eigenvalues cluster around 10810^8 to 10910^9, with LVL_V slightly largest, producing a more even distribution of gradient contributions. This more balanced distribution, together with the larger tail eigenvalues, is consistent with the time-domain results: KAN and ActNet learn the fast oscillatory component more effectively, whereas the vanilla PINN fails to capture it accurately. Figure 8 shows the NTK spectra, at epoch 0, after applying the adaptive weighting rule formulated in (15), which scales each loss term according to the inverse trace of its NTK block. This reduces disparities in the effective contribution of the residual terms and promotes more similar convergence rates across the state-equation residuals. Figure 8 compares the complete ActNet NTK spectra at initialization and after 1000, 3000, and 5000 training iterations. The horizontal axis gives the eigenvalue rank in decreasing order, while the legend identifies the training checkpoint. Their strong overlap suggests that the NTK remains approximately fixed during training. V-C Regime 3: Inverter Model Benchmark The final case study considers an 11-state PQ-controlled, grid-following inverter model. The GPU solver evaluates 1000 scenarios in approximately 27ms27\,ms, making this benchmark relevant for acceleration. Although its stiffness ratio, rinv≈1.7×102r_inv≈ 1.7× 10^2, is smaller than that of the modified SMIB, rSMIB≈2×103r_SMIB≈ 2× 10^3, the inverter contains several interacting current-control, power-control, PLL, and filtering modes with different timescales [3]. We therefore examine whether this distributed multi-timescale structure creates optimization difficulties that are not captured by the stiffness ratio alone. V-C1 Inverter Model and Training The inverter model includes a PLL, outer PQ control, inner current control, current dynamics, and measurement filters, as detailed in [3]. We train a vanilla PINN, ActNet, and KAN from a single initial condition for 10,00010,000 epochs with a learning rate of 10−410^-4, using the physics-informed objective in (I-B). The PINN and ActNet use three hidden layers, while the KAN uses two hidden layers of 20 neurons each. For the NTK comparison in Fig. 9, we additionally evaluate PINNs with widths 64, 128, 256, and 512. The PLL and outer-control states evolve substantially more slowly than the current dynamics, producing large differences in residual and gradient scales. We denote the residual loss of state z by LzL_z and use a prediction horizon of 0.1s0.1\,s, unless stated otherwise. V-C2 NTK-Based Diagnosis We diagnose the training difficulty at three levels: the per-residual NTK spectra identify slowly learned error modes, the residual gradients reveal imbalance between loss terms, and a short-horizon ablation separates optimization difficulty from the ability of the N to fit the trajectory. images/inverter_ntk_pic_new Figure 9: NTK eigenvalue spectra of the diagonal loss blocks Ki,iK_i,i for selected inverter residuals and different architectures. The residuals represent slow PLL dynamics, fast current dynamics, and intermediate filtering dynamics. The horizontal axis denotes the eigenvalue rank after sorting by decreasing magnitude. Figure 9 compares the slow PLL residual LξpllL_ _pll, the fast current residuals LidL_i_d and LiqL_i_q, and the intermediate filter residual LQfiltL_Q_filt. All architectures exhibit long tails of small eigenvalues, meaning that many error modes are only weakly corrected by parameter updates. Increasing the PINN width raises some eigenvalues but does not remove these slowly learned modes, while ActNet and KAN show the same qualitative limitation. Thus, increasing model size or changing among the tested architectures does not by itself resolve the training difficulty. images/gradient_imbalance Figure 10: Maximum residual gradients with respect to the N final-layer weights and the corresponding NTK-based weights from (15). The adaptive weights reduce the dominance of residuals with large gradients. Figure 10 examines the imbalance between residual losses. The current residuals LidL_i_d and LiqL_i_q produce the largest gradients and therefore receive the smallest adaptive weights, while weaker residuals receive larger weights. This weighting rebalances the influence of the complete loss terms, but each weight scales all modes of a loss together and therefore does not remove its small-eigenvalue tail. V-C3 Short-Horizon Ablation To determine whether the remaining error arises from limited trajectory representation or from the physics-informed optimization, we train a simplified one-layer N over 10−4s10^-4\,s using five configurations: data loss only; physics loss with and without NTK weighting; and combined data and physics losses with and without NTK weighting. Figure 11: Short-horizon ablation comparing supervised and physics-informed training. The data-only model provides a supervised reference, while NTK weighting improves the combined data-and-physics objective relative to its unweighted counterpart. Figure 11 shows that the data-only model achieves a much lower trajectory error, indicating that, over this short horizon, fitting direct state observations is easier than satisfying the derivative-based physics residuals. The N can therefore represent the local trajectory, while the physics-informed objective creates the main optimization difficulty. NTK weighting improves the combined data-and-physics case relative to unweighted training, although it remains less accurate than direct supervised fitting. Nevertheless, physics-informed training provides an explicit and physically interpretable objective and avoids reliance on labeled trajectories generated by numerical simulation; the result therefore motivates improving its optimization rather than replacing it with supervised training. Overall, the inverter is difficult to learn because its stiffness is distributed across several interacting fast and slow modes. This produces both slowly learned modes within individual residuals and strong imbalance between residual losses, making the inverter harder to learn than the modified SMIB despite its smaller stiffness ratio. These findings suggest future work on architectures with separate representations for fast current dynamics and slower control and filtering states, together with state-wise normalization and mode-aware weighting. Developing and validating such NTK-guided architectures is left for future work, with the goal of replacing trial-and-error model selection by systematic, physics-aware architecture and training design. VI Conclusion This paper applies the Neural Tangent Kernel (NTK) method to assess, interpret, and improve the training of neural surrogate models for power-system dynamics. As physics-informed machine learning is increasingly used to model aggregate power-system behavior, accelerate simulations of complex components such as inverters, and support black-box converter surrogates, understanding how architecture and hyperparameters affect performance becomes essential. We show that classical notions of numerical stability in differential-algebraic systems reappear as convergence, conditioning, and robustness issues in N training. Through its analogy with small-signal eigenvalue analysis, NTK provides a modal interpretation of learning dynamics, showing how states with widely different timescales are learned and guiding adaptive loss weighting. Focusing on individual power-system components, we studied physics-informed learning of stiff dynamics using vanilla PINNs, KANs, and ActNets. On SM models, KANs and ActNets captured fast and slow states more evenly than vanilla PINNs, with ActNet providing the best accuracy–runtime trade-off. For the 11-state inverter model, the results show that learnability depends less on the overall stiffness range than on the distribution of stiff directions. Although the inverter has a smaller stiffness ratio than the modified SMIB benchmark, its electrical, control, and synchronization dynamics create a highly anisotropic eigenvalue distribution that the NTK inherits. This makes training accurate physics-informed surrogates difficult, and neither NTK-based adaptive weighting nor the tested N architectures fully resolve the issue. Future work should use NTK-guided architecture design, loss weighting, and optimization strategies to address stiff multi-timescale power-system dynamics and produce fast, accurate neural surrogates. References [1] M. Bossart, J. D. Lara, C. Roberts, R. Henriquez-Auba, D. S. Callaway, and B. S. Hodge (2025) Acceleration of power system dynamic simulations using a deep equilibrium layer and neural ode surrogate. IEEE Transactions on Energy Conversion 40 (4), p. 2710–2722. External Links: Document Cited by: §I. [2] P. Ellinas, I. Karampinis, I. V. Nadal, R. Nellikkath, J. Vorwerk, and S. Chatzivasileiadis (2025) Physics-informed machine learning for power system dynamics: a framework incorporating trustworthiness. Sustainable Energy, Grids and Networks 43, p. 101818. External Links: Document Cited by: §I, §I, §I-C, §V-A1. [3] () GitHub - elpetros99/tools_explaining_power: Appendix for the paper ”Tools to Explain Neural Networks for Power System Dynamics” — github.com. Note: https://github.com/elpetros99/tools\_explaining\_power.git[Accessed 10-04-2026] Cited by: §I, §V-A1, §V-C1, §V-C, §V. [4] E. Hairer and G. Wanner (1996) Solving ordinary differential equations i: stiff and differential-algebraic problems. 2 edition, Springer Series in Computational Mathematics, Vol. 14, Springer, Berlin, Heidelberg. External Links: Document, ISBN 978-3-540-60452-5 Cited by: §I. [5] A. Jacot, F. Gabriel, and C. Hongler (2020) Neural tangent kernel: convergence and generalization in neural networks. External Links: 1806.07572, Link Cited by: §I, §IV. [6] D. P. Kingma and J. Ba (2015) Adam: a method for stochastic optimization. In 3rd International Conference on Learning Representations (ICLR), Cited by: §V-B2. [7] J. D. Lara, R. Henriquez-Auba, D. Ramasubramanian, S. Dhople, D. S. Callaway, and S. Sanders (2024) Revisiting power systems time-domain simulation methods and models. IEEE Transactions on Power Systems 39 (2), p. 2421–2437. External Links: Document Cited by: §I. [8] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljačić, T. Y. Hou, and M. Tegmark (2025) KAN: kolmogorov-arnold networks. In The Thirteenth International Conference on Learning Representations, External Links: Link Cited by: §I-C. [9] Z. Liu, Y. Wang, S. Vaidya, F. Ruehle, J. Halverson, M. Soljacic, T. Y. Hou, and M. Tegmark (2025) KAN: kolmogorov–arnold networks. In International Conference on Learning Representations (ICLR), External Links: Link Cited by: §I-C. [10] U. Markovic, O. Stanojev, P. Aristidou, E. Vrettos, D. Callaway, and G. Hug (2021) Understanding small-signal stability of low-inertia systems. IEEE Transactions on Power Systems 36 (5), p. 3997–4017. External Links: Document Cited by: §I. [11] G. S. Misyris, A. Venzke, and S. Chatzivasileiadis (2020) Physics-informed neural networks for power systems. In 2020 IEEE Power & Energy Society General Meeting (PESGM), Vol. , p. 1–5. External Links: Document Cited by: §I. [12] () Regulation - EU - 2024/1689 - EN - EUR-Lex — eur-lex.europa.eu. Note: https://eur-lex.europa.eu/eli/reg/2024/1689/oj/eng[Accessed 24-01-2025] Cited by: §I. [13] J. Stiasny and S. Chatzivasileiadis (2023) Physics-informed neural networks for time-domain simulations: accuracy, computational cost, and flexibility. Electric Power Systems Research 224, p. 109748. External Links: ISSN 0378-7796, Document, Link Cited by: §I. [14] S. Subedi, M. Rauniyar, S. Ishaq, T. M. Hansen, R. Tonkoski, M. Shirazi, R. Wies, and P. Cicilio (2021) Review of methods to accelerate electromagnetic transient simulation of power systems. IEEE Access 9 (), p. 89714–89731. External Links: Document Cited by: §I. [15] N. Vyas, D. Morwani, R. Zhao, M. Kwun, I. Shapira, D. Brandfonbrener, L. Janson, and S. Kakade (2025) SOAP: improving and stabilizing shampoo using adam. External Links: 2409.11321, Link Cited by: §V-A2. [16] S. Wang, Y. Teng, and P. Perdikaris (2020) Understanding and mitigating gradient pathologies in physics-informed neural networks. External Links: 2001.04536, Link Cited by: §I-C. [17] S. Wang, X. Yu, and P. Perdikaris (2022) When and why PINNs fail to train: a neural tangent kernel perspective. Journal of Computational Physics 449, p. 110768. External Links: Document Cited by: §IV. [18] L. Ziyin, T. Hartwig, and M. Ueda (2020) Neural networks fail to learn periodic functions and how to fix it. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33, p. 1583–1594. External Links: Link Cited by: §I-C.