Paper deep dive

Accounting for Hysteresis and Eddy Currents in Finite Element Simulations of Ferromagnetic Laminated Cores using a Recurrent Neural Network

Florent Purnode, Louis Denis, François Henrotte, Gilles Louppe, Christophe Geuzaine

Year: 2026Venue: arXiv preprintArea: cs.CEType: PreprintEmbeddings: 106

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%

Last extracted: 7/17/2026, 5:16:45 AM

Summary

This paper introduces a recurrent neural network (RNN) surrogate model to efficiently simulate hysteresis and eddy currents in ferromagnetic laminated cores. By training the RNN on artificially generated magnetic field sequences derived from a reference energy-based model and 1D magneto-quasistatic simulations, the approach achieves high accuracy while limiting computational costs to approximately twice that of anhysteretic simulations. The model is designed for seamless integration into 2D magnetodynamic finite element frameworks and is publicly available.

Entities (10)

Eddy Currents · phenomenon · 99%Finite Element Simulation · method · 99%Hysteresis · phenomenon · 99%Recurrent Neural Network · model-architecture · 98%University of Liege · organization · 98%Energy-Based Model · model · 97%Laminated Core · component · 97%Gated Recurrent Unit · model-architecture · 96%Machine Learning Surrogate Modeling · technique · 94%M235-35A Steel · material · 92%

Relation Signals (7)

Proposed Approach → integratedinto → Finite Element Simulations

confidence 98% · is integrated into realistic two-dimensional magnetodynamic finite element simulations based on a magnetic vector potential formulation

Energy-Based Model → accountsfor → Hysteresis

confidence 97% · This hysteresis model is based on energy considerations and is intrinsically vectorial

1D FEM Formulation → accountsfor → Eddy Currents

confidence 96% · accounts for the eddy currents and the associated skin effect in the thickness direction via the resistivity term

Gated Recurrent Unit → istypeof → Recurrent Neural Network

confidence 96% · The GRU is adopted here, as it achieves accuracy comparable to LSTMs with fewer parameters

Recurrent Neural Network → trainedas → Surrogate Model

confidence 95% · a recurrent neural network is trained as a surrogate of a laminated-core material model

Surrogate Model → limits → Computational Cost

confidence 94% · limiting the computational cost to about twice that of an anhysteretic simulation

Training Dataset → consistsof → Artificially Generated Magnetic Field Sequences

confidence 93% · training the recurrent neural network on a sufficiently diverse set of artificially generated magnetic field sequences

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Incorporating hysteresis and eddy currents into finite element simulations of laminated-core electrical machines is computationally challenging. Resolving the fields inside the laminations at each integration point and at every nonlinear iteration leads to computational costs several orders of magnitude higher than anhysteretic simulations, making such approaches impractical for design applications. Conversely, simplified models accounting only for magnetic saturation are becoming increasingly inadequate as electrical machine topologies and operating conditions grow in complexity. In this context, machine learning surrogate modeling has emerged as a promising alternative, offering efficient and accurate approximations of complex electromagnetic behaviors. In this paper, a recurrent neural network is trained as a surrogate of a laminated-core material model for an isotropic laminated core, and is integrated into realistic two-dimensional magnetodynamic finite element simulations based on a magnetic vector potential formulation. The proposed approach achieves excellent agreement with the reference laminated-core model while limiting the computational cost to about twice that of an anhysteretic simulation. By training the recurrent neural network on a sufficiently diverse set of artificially generated magnetic field sequences designed to mimic those encountered in electrical machine simulations, the proposed approach can be readily applied across a wide range of finite element simulations. Furthermore, the trained surrogate model is provided as a standalone component that can be easily incorporated into existing computational frameworks. It is publicly available at this https URL.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

104,891 characters extracted from source content.

Expand or collapse full text

Accounting for Hysteresis and Eddy Currents in Finite Element Simulations of Ferromagnetic Laminated Cores using a Recurrent Neural Network Florent Purnode florent.purnode@uliege.be Louis Denis François Henrotte Gilles Louppe Christophe Geuzaine cgeuzaine@uliege.be Abstract Incorporating hysteresis and eddy currents into finite element simulations of laminated-core electrical machines is computationally challenging. Resolving the fields inside the laminations at each integration point and at every nonlinear iteration leads to computational costs several orders of magnitude higher than anhysteretic simulations, making such approaches impractical for design applications. Conversely, simplified models accounting only for magnetic saturation are becoming increasingly inadequate as electrical machine topologies and operating conditions grow in complexity. In this context, machine learning surrogate modeling has emerged as a promising alternative, offering efficient and accurate approximations of complex electromagnetic behaviors. In this paper, a recurrent neural network is trained as a surrogate of a laminated-core material model for an isotropic laminated core, and is integrated into realistic two-dimensional magnetodynamic finite element simulations based on a magnetic vector potential formulation. The proposed approach achieves excellent agreement with the reference laminated-core model while limiting the computational cost to about twice that of an anhysteretic simulation. By training the recurrent neural network on a sufficiently diverse set of artificially generated magnetic field sequences designed to mimic those encountered in electrical machine simulations, the proposed approach can be readily applied across a wide range of finite element simulations. Furthermore, the trained surrogate model is provided as a standalone component that can be easily incorporated into existing computational frameworks. It is publicly available at https://gitlab.onelab.info/getdp/lamnet. keywords: Finite Element Simulation , Homogenization , Hysteresis , Iron Losses , Laminated Core , Neural Network †journal: Journal of Computational Physics [1] organization=Department of Electrical Engineering and Computer Science, University of Liege, country=Belgium 1 Introduction Transformers and electric machines rely on laminated ferromagnetic cores to guide and concentrate magnetic flux. Accurate modeling of electromagnetic fields and magnetic losses within these laminated regions is crucial for the reliable prediction of induced voltages and overall device efficiency [1, 2]. However, the electromagnetic behavior of laminated cores is particularly challenging to model due to the inhomogeneous distribution of eddy currents, the nonlinear and history-dependent nature of magnetic hysteresis, and the coupling between these two phenomena. Modeling hysteresis alone is already challenging. The Preisach [3] and Jiles-Atherton [4] models have long been considered standard hysteresis models, and have been widely implemented in two-dimensional (2D) finite element (FE) simulations [5, 6, 7]. However, both models were originally formulated for scalar hysteresis, and their extensions to vector hysteresis lack a clear physical foundation which limits their applicability to vector field problems. In addition, they may suffer from limitations related to accuracy, computational cost, or numerical convergence [8, 9]. More recently, the energy-based (EB) model [10, 11] has emerged as a promising alternative. This hysteresis model is based on energy considerations and is intrinsically vectorial [12]. The forward EB model takes the magnetic field H as input and returns the history-dependent magnetic flux density B, while its inverse formulation [13, 14], computes H from a given B. The forward model can be evaluated either by solving a minimization problem [15] or a nonlinear vector function [11]. Besides these two formulations, a simplified explicit variant, akin to a vector play model [10, 11], also exists and provides a computationally efficient update at the cost of a slightly reduced accuracy [16]. The inverse formulation, by contrast, requires solving a minimization problem [11, 15, 17], whose robustness has been recently significantly improved [13, 14]. Eddy currents further add to the complexity of modeling ferromagnetic behavior. They form closed current loops whose distribution is governed by the spatial and temporal variations of the magnetic flux density. In laminated cores, eddy currents are largely confined to the plane of each lamination, but three-dimensional (3D) effects arise near edges and geometric discontinuities and cannot be fully captured by the 2D formulations commonly employed in electrical machine design. However, direct 3D modeling of laminated cores, in which the fields are resolved within each individual lamination, is several orders of magnitude more computationally expensive. Although feasible on high-performance computing platforms, such simulations remain impractical for engineering design and optimization. The periodicity of the laminated structure can, however, be exploited through homogenization techniques, yielding 2D formulations that accurately capture the bulk eddy-current behavior, although their accuracy deteriorates near edges and geometric discontinuities. In multiscale methods like the heterogeneous multiscale method (HMM) [18, 19] or the multilevel finite element method (FE2) [20], the effective material response unavailable at the macroscopic scale is computed during the simulation by solving a mesoscopic-scale problem on a representative volume element and homogenizing its solution. In the present context, the mesoscopic-scale model can be formulated as a magneto-quasistatic FE simulation solved on the cross-section of a lamination and incorporating the EB hysteresis model. However, this approach remains too computationally demanding for routine machine design for three main reasons: i) The mesoscopic-scale model must be queried in each macroscopic finite element, at each time step and at each nonlinear iteration. i) The mesoscopic-scale model must additionally provide the corresponding Jacobian, ∂/∂ / or ∂/∂ / , to construct the macroscopic tangent stiffness matrix required for a Newton-Raphson scheme. If obtained by finite differences, the mesoscopic-scale problem must be solved for the nominal state, and for additional perturbed states, further increasing the number of mesoscopic-scale model queries. i) Each mesoscopic-scale query itself requires evaluating the EB model at every element of the mesoscopic mesh. Advances in machine learning have led to a rapid expansion of neural network (N)-based methods dedicated to magnetic modeling. Among other applications, NNs are used to predict core losses [21], to predict hysteresis model parameters [22, 23], or to predict B from H and inversely [24, 25, 26, 27]. This enthusiasm is further evidenced by the emergence of dedicated international competitions in the power electronics community, such as the MagNet challenges, which have been specifically designed to stimulate innovation in the data-driven modeling of core losses [28] and transient magnetic behavior [29]. Yet, despite an intense research activity, and proven effectiveness in related fields such as heterogeneous elasto-plastic materials [30], integrating an N-based surrogate model to account for vector hysteresis and eddy currents in 2D FE simulations of realistic electrical machines remains an open problem. In a previous contribution [31], the relation required for magnetic vector-potential formulations was approximated by a parametric homogenized constitutive law ~(,˙,pk) H(B, B,p_k). The parameters pkp_k of this law were identified element-wise through a three-step procedure: i) A 2D macroscopic-scale anhysteretic FE simulation was performed to compute characteristic periodic (t)H(t) sequences. i) Corresponding (t)B(t) responses were obtained by feeding the resulting (t)H(t) sequences to the aforementioned mesoscopic-scale model using the explicit variant of the forward EB model, and homogenizing its solution. i) The parameters pkp_k were then identified by fitting the homogenized constitutive law ~(,˙,pk) H(B, B,p_k) to the computed ((t),(t))(B(t),H(t)) pairs. This identification was carried out through N training, which minimizes the fitting error over all sequences simultaneously and yields a smooth spatial distribution of the parameters pkp_k across the macroscopic domain. The resulting macroscopic 2D FE simulations efficiently captured the combined effects of hysteresis and eddy currents through the parametric constitutive law ~(,˙,pk) H(B, B,p_k). However, this approach is restricted to periodic steady-state operating conditions. Furthermore, the fixed analytical structure of the parametric law inherently limits the achievable accuracy, and the resulting macroscopic FE simulations may converge slowly. In the present work, rather than assuming a fixed parametric structure, an effective constitutive relation ~() H(B) is directly learned from precomputed data generated by the aforementioned mesoscopic-scale model. A recurrent neural network (RNN) is used to naturally handle the history-dependent nature of hysteresis, lifting the periodicity restriction of the previous approach, while the Jacobian ∂~/∂ H/ is obtained exactly through automatic differentiation [32]. The improved approach is evaluated by modeling two real-world applications. It demonstrates excellent accuracy and robust convergence properties, while maintaining computational times comparable to those of conventional anhysteretic simulations. In addition, the proposed N-based model provides an estimation of its error, which constitutes valuable information for assessing the reliability of the computed solution. The paper is organized as follows. Section 2 first presents the mesoscopic-scale model used to generate the training datasets and to assess the N-based model accuracy. Section 3 then details the model architecture, training dataset and training procedure, while Section 4 presents the conducted 2D FE simulations. Finally, Section 5 provides conclusions and perspectives for future work. 2 Lamination Model This section presents the detailed lamination model used as reference throughout this paper. It provides a homogenized description of the fields at the mesoscopic scale, i.e., at the scale of the lamination thickness. The model combines the explicit variant of the EB hysteresis model [10] with a one-dimensional (1D) magnetodynamic FE formulation accounting for eddy currents. The EB model is first introduced, followed by the associated 1D FE simulation. 2.1 Quasi-Static Mesoscopic Model: The Energy-Based Hysteresis Model The EB hysteresis model is grounded in thermodynamic principles; its theoretical derivation is presented in detail in [11, 15]. Without expanding on this derivation, the model captures the distributed nature of pinning strengths through an ensemble of N elementary cells. Each cell n decomposes the mesoscopic magnetic field h into a reversible component revnh_rev^n and an irreversible component irrnh_irr^n, whose magnitude is constrained within a threshold κnκ^n. As long as |irrn|≤κn|h_irr^n|≤κ^n, the cell behaves reversibly. Once the threshold is exceeded, the cell transitions into irreversible behavior. Among the update laws proposed to describe this transition, we adopt the explicit variant, illustrated on Fig. 1 and expressed as revn=rev,pnif |−rev,pn|≤κn,−κn−rev,pn|−rev,pn|otherwise,h_rev^n= casesh_rev,p^n&if |h-h_rev,p^n|≤κ^n,\\ h-κ^n h-h_rev,p^n|h-h_rev,p^n|&otherwise, cases (1) where rev,pnh_rev,p^n denotes the previous reversible component revnh_rev^n. From the weighted sum of the reversible fields rev=∑n=1Nwnrevnh_rev= _n=1^Nw^nh_rev^n, with the wnw^n parameters summing up to 1, the mesoscopic magnetic flux density is expressed as =μ0+[ℒ(|rev|ha)Ja|rev|+ℒ(|rev|hb)Jb|rev|]rev,b= _0h+ [ L ( |h_rev|h_a ) J_a|h_rev|+ L ( |h_rev|h_b ) J_b|h_rev| ]h_rev, (2) where μ0 _0 is the magnetic permeability of vacuum, ℒ(x)=coth⁡(x)−1/x L\,(x)= \,(x)-1/x denotes the Langevin function, and hah_a, hbh_b, JaJ_a and JbJ_b are parameters of the EB model. These parameters, along with the κnκ^n and wnw^n parameters, are determined from measurements, and have to be updated whenever the material changes. This work employs N=11N=11 elementary cells for the modeling of M235-35A non-oriented electrical steel sheets, with the identified parameter values reported in A. Figure 1: Illustration of (1). (a) h remains in the circle of radius κnκ^n centered in rev,pnh_rev,p^n, so that revnh_rev^n remains unchanged. (b) h exits the circle, so that revnh_rev^n moves along the line linking rev,pnh_rev,p^n and h, and stop at a distance κnκ^n from h. 2.2 Dynamic Mesoscopic Model: Laminations with Eddy Current and Hysteresis The aforementioned EB model solely accounts for quasi-static hysteresis. To also account for dynamic effects, it is included into the FE simulation of a stack of ferromagnetic laminations. Considering that all the laminations see the same imposed macroscopic magnetic field H, and neglecting edge effects, simulating a single lamination is sufficient to describe the whole stack. Furthermore, the mesoscopic magnetic field h and the mesoscopic magnetic flux density b inside a lamination can reasonably be assumed parallel to the lamination plane and to vary exclusively within the lamination depth. Eventually, the symmetry with respect to the median plane of a lamination is exploited to simulate solely half of a lamination thickness d (see Fig. 2). Figure 2: Definition of the domain of analysis of the 1D lamination model. Let t denote time, ρ the electrical resistivity, and (x,y,z)(x,y,z) a Cartesian coordinate system, with x and y lying in the plane of the laminations and z in the stacking direction. Furthermore, let z=0z=0 correspond to the lamination surface and z=d/2z=d/2 to its center. The 1D magneto-quasistatic formulation then reads: Find the mesoscopic magnetic field (z,t)=(hx(z,t),hy(z,t))h(z,t)=(h_x(z,t),h_y(z,t)) such that ∂t(bx(),hx′)[0,d/2]+(ρ∂zhx,∂zhx′)[0,d/2]=0∂t(by(),hy′)[0,d/2]+(ρ∂zhy,∂zhy′)[0,d/2]=0 cases _t (b_x (h ),h_x )_ [0,d/2 ]+ (ρ\, _zh_x, _zh_x )_ [0,d/2 ]&=0\\ _t (b_y (h ),h_y )_ [0,d/2 ]+ (ρ\, _zh_y, _zh_y )_ [0,d/2 ]&=0 cases (3) holds for all test functions hx′h_x , hy′h_y , where the notation (⋅,⋅)D(·,·)_D denotes the inner product of its arguments over bounded domain D. The macroscopic excitation (t)H(t) is imposed at the lamination surface, i.e. (0,t)=(t)h(0,t)=H(t), and the zero-current condition at the mid-plane leads to the Neumann boundary condition ∂z(d/2,t)=0 _zh(d/2,t)=0. The macroscopic magnetic flux density (t)B(t) is retrieved by averaging the mesoscopic magnetic flux density ()=(bx(),by())b(h)=(b_x(h),b_y(h)) over the half lamination thickness. This formulation accounts for the eddy currents and the associated skin effect in the thickness direction via the resistivity term ρ∂zρ\, _zh, and accounts for hysteresis via the ()b(h) term, solution of the EB model. In practice, the model is implemented as a self-contained C language code, solving the linear system, whose matrix is sparse, using an efficient LU sparse algorithm. The code is applied to a stack of M235-35A non-oriented electrical steel laminations, whose domain of analysis z∈[0,d/2]z∈ [0,d/2 ] is discretized with 51 equidistant nodes. Given that the maximum differential permeability is μmax=limh→0+(μ0+ℒ(h/ha)Ja/h+ℒ(h/hb)Jb/h)≈25.5 _max= _h→ 0^+( _0+ L(h/h_a)\,J_a/h+ L(h/h_b)\,J_b/h)≈ 25.5 mH/m, obtained from (2), the resulting mesh size d/100=3.5μd/100=3.5~ ensures that the skin depth δ=ρ/(πfμmax)δ= ρ/(π f _max) is resolved by at least five nodes up to approximately 44 kHz. Thus, skin-effect phenomena are accurately captured over the frequency range considered in this work (see Section 3.2). The temporal discretization is chosen sufficiently fine to avoid nonlinear iterations, and the time derivatives are discretized with the implicit Euler scheme. Finally, the cells of the EB model are initialized assuming an initially demagnetized material subjected to a slowly increasing excitation field until reaching (t=0)H(t=0), such that no dynamic effects are initially present and ((t=0))B(H(t=0)) lies on the virgin magnetization curve. 3 Architecture and Training of the N-Based Surrogate Model The choice of the N architecture is of predominant importance. Since hysteresis inherently depends on the material history, the N model must operate on time-sequence data. Some authors process the entire sequence at once using, for example, neural operators [33]. In FE simulations, however, the input sequence length evolves continuously during time integration, which makes whole-sequence architectures poorly suited for direct integration into incremental solvers. Another approach consists in using a fixed-length sliding window that captures only the most recent history [34, 35, 36]. Nevertheless, such a fixed-length window may omit parts of the history that significantly influence the response, thereby reducing accuracy. Alternatively, attention-based architectures [37] have shown high accuracy on sequence-to-sequence tasks [38], but the computational overhead becomes significant for long sequences. In this contribution, an RNN architecture is preferred. Indeed, RNNs are well suited to handle hysteresis in FE simulations: (i) they maintain and update state variables similarly as the EB model maintains and updates the reversible fields revnh_rev^n; (i) they are fast to infer a single time step, making them efficient within incremental solvers; and (i) they achieve high accuracy [39, 40, 41]. This section discusses the implementation of the N-based model that will be integrated into the FE simulations of Section 4, and details its training process. At first, Section 3.1 presents the general model architecture and Section 3.2 details the datasets used for training. Section 3.3 later demonstrates that the model cannot be trained naively to minimize the mean squared error, since saturated fields will dominate the training to the detriment of lower fields where the magnetic behavior is more intricate. Section 3.3 therefore introduces a scaled loss, which is further modified in Section 3.4 to assess the model confidence. In Section 3.5, the model is upgraded to ensure correct prediction even in extreme saturation, which lies out of the training range. Section 3.6 later describes how the training sequences, initially defined on constant time steps, are resampled to include variable time steps, thereby enabling adaptive time-stepping strategies in the FE simulations. Finally, Section 3.7 discusses the hyperparameter tuning, while Section 3.8 assesses the model performance, prior to any FE coupling. 3.1 Baseline Architecture To process temporal sequences, RNNs ingest a single time step at a time. At each time step k, the network updates a hidden state SkS^k based on the current input XkX^k and the previous hidden state Sk−1S^k-1. In its simplest form, this update is given by Sk=ϕ(WxXk+WsSk−1+b)S^k=φ(W_x\,X^k+W_s\,S^k-1+b), where ϕ(⋅)φ(·) is a nonlinear activation function, WxW_x and WsW_s are weight matrices, b is a bias term, and the output YkY^k is computed directly from the hidden state SkS^k. This basic formulation however suffers from the vanishing- and exploding-gradient problems, which impair its ability to learn long-range dependencies. To address this issue, advanced recurrent units have been developed, the most prominent being the Long Short-Term Memory (LSTM) [42] and the Gated Recurrent Unit (GRU) [43]. The GRU is adopted here, as it achieves accuracy comparable to LSTMs with fewer parameters and a simpler gating structure, making it both easier to train and more efficient to deploy. Without going into the details of the GRU cell, Fig. 3 presents the basic model architecture, which is further enhanced in Sections 3.4 and 3.5. At time step k, the model takes as input the x- and y-components of the magnetic flux density kB^k, together with the time increment Δtk t^k. The latter may vary between successive iterations to enable temporal refinement during FE simulations. Both kB^k and Δtk t^k are scaled, concatenated, and passed through an encoder block implemented as a multilayer perceptron, which progressively increases the feature dimensionality to match the GRU hidden-state size. The GRU cell is followed by a decoder block, also implemented as a multilayer perceptron, which maps the hidden state back to the x- and y-components of ~k H^k, that is the predicted magnetic field, intended to approximate the target magnetic field kH^k. Rather than initializing the hidden state to zero, it is initialized by a trainable vector, plus the output of a shallow network taking 0B^0 and Δt0 t^0 as inputs. This offloads the initial magnetization response, which is history-independent and lies on the virgin curve (see Section 2.2), to a dedicated subnetwork. The GRU cell therefore starts with a meaningful hidden state, and can devote its full capacity to learning history-dependent evolution, reducing the number of parameters that would otherwise be required. Figure 3: Baseline model architecture, enhanced in Sections 3.4 and 3.5. 3.2 Generation of the Training Datasets In order to propose an N-based surrogate model readily usable in a wide range of applications, the training dataset should be populated with a broad variety of field profiles. Since using data from simulations or measurements on existing devices may bias the model toward specific applications, the training H sequences are generated artificially to span the diversity of waveforms encountered in electrical machines. The corresponding B sequences are later obtained by feeding them to the lamination model of Section 2. As an initial strategy, the training H sequences can be expressed as a sum of harmonics in which each component is defined by two orthogonal contributions with random amplitudes and phases: (t)=∑h∈ℋαhH⋆h[cos⁡(βh)sin⁡(2πfht+φhx)^+sin⁡(βh)sin⁡(2πfht+φhy)^].H(t)= _h _hH h [ ( _h) (2π fht+ _h^x)\, x+ ( _h) (2π fht+ _h^y)\, y ]. (4) In this expression, the set ℋH is generated by first randomly selecting the maximum harmonic hmax∈[1,50]h_max∈[1,50], later randomly selecting the number of harmonics |ℋ|∈[1,hmax]|H|∈[1,h_max], and finally randomly selecting the remaining harmonic indices without replacement in 1,…,hmax−1\1,…,h_max-1\. The characteristic amplitude H⋆H is drawn from a log-uniform distribution over [10A/m,100kA/m][10~A/m,100~kA/m], while the fundamental frequency f is uniformly sampled between 0.1 Hz and 1 kHz, with harmonic h having frequency hfh\,f. The selected ranges for H⋆H and hmaxh_max are justified in Sections 3.5 and 3.6, respectively, where the associated anhysteretic extension and resampling strategy are introduced. The coefficients αh∈[0,1] _h∈[0,1] are also uniformly sampled and, together with the decay factor 1/h1/ h, modulate the relative amplitude of each harmonic, thereby controlling the harmonic distribution and enabling the generation of both smooth and highly distorted waveforms. The parameters βh∈[0,2π] _h∈[0,2π] in turn govern how each harmonic contribution is shared between the orthogonal directions x and y, thus defining its orientation in the plane. Eventually, the phase terms φhx _h^x and φhy∈[0,2π] _h^y∈[0,2π] independently shift the temporal evolution of each directional contribution, while their difference determines the trajectory of each harmonic component in the plane: identical phases lead to unidirectional oscillations, whereas a non-zero phase difference produces a rotating component, resulting in circular or elliptical trajectories. Up to now, (4) describes a broad class of excitations, encompassing alternating, rotating, and distorted magnetic field trajectories. However, it poorly accounts for several features commonly encountered in electrical machines: i) purely unidirectional fields, i) DC-biased fields, i) pulse-like fields, iv) and ramp-up fields. To improve the representation of these excitation types, several independent modifications are introduced. Each modification is activated with a probability of 50 %, avoiding both under-representation of the corresponding waveform characteristics and excessive bias toward them. i) Unidirectional Fields Without specific constraints on the parameters βh _h, φhx _h^x, and φhy _h^y, (4) almost never generates purely unidirectional fields. Such fields arise only when all βh _h are identical and when φhx _h^x and φhy _h^y are either equal or differ by an integer multiple of π for all h∈ℋh . However, unidirectional fields are commonly encountered in electrical machines. To ensure their adequate representation in the training dataset, all βh _h are assigned the same value with a probability of 50 %. Independently, the phase terms are generated such that φhx=φhy _h^x= _h^y for all h∈ℋh , again with a probability of 50 %. Since these two conditions are imposed independently, approximately 25 % of the generated sequences are purely unidirectional. i) DC-biased Fields Among other environments, DC-biased fields are encountered in the rotor of synchronous machines, in the stator of DC machines and in the presence of permanent magnets. With a probability of 50 %, (4) is therefore augmented by a DC component dc=Hdc(cos⁡(βdc)^+sin⁡(βdc)^)H_dc=H_dc\, ( ( _dc )\, x+ ( _dc )\, y ) (5) where both the amplitude Hdc∈[0,H⋆]H_dc∈[0,H ] and direction βdc∈[0,2π] _dc∈[0,2π] are uniformly sampled. i) Pulse-Like Fields In rotating electrical machines, the slotted geometry of the stator and rotor often causes the air-gap reluctance to vary periodically. As a result, the field near the air-gap takes a pulse-like shape: high when rotor and stator teeth are facing each other, and significantly reduced when a tooth faces a slot. In principle, (4) can build such pulse-like fields, but it is unlikely. To ensure adequate coverage of this field shape in the dataset, a pulse-like modulation is applied with a probability of 50 %. When applied, the x and y components of both (4) and dcH_dc are scaled by αpx+(1−αpx)p(t) _p^x+(1- _p^x)\,p(t) and αpy+(1−αpy)p(t) _p^y+(1- _p^y)\,p(t), respectively, where αpx _p^x and αpy∈[0,1] _p^y∈[0,1] are uniformly sampled, and p(t)=sin(2πfhpt+φp)epp(t)= (2π fh_pt+ _p)^e_p (6) is a pulse-like function. The number of pulses per period hp∈[1,20]h_p∈[1,20] and the phase φp∈[0,2π] _p∈[0,2π] are both drawn from uniform distributions. The upper bound on hph_p was chosen to ensure the resulting pulses are sufficiently resolved, while still allowing variability within each pulse owing to the harmonic content of (4). On the other hand, the exponent epe_p, sampled as a random integer between 3 and 20, controls the sharpness and symmetry of the pulses. When epe_p is odd, p(t)p(t) preserves the sign of the sine term, resulting in alternating positive and negative pulses, as can be observed in the stator teeth. Conversely, when epe_p is even, the function remains non-negative, producing pulses of constant sign, as is typical in the rotor near the air gap, a behavior that will be illustrated in Fig. 17 at the end of Section 4.2. iv) Ramp-up Transients The magnetic excitation can progressively build up from an initially demagnetized state. To account for such excitations, the generated field is multiplied, with a probability of 50 %, by the ramp-up function r(t)=1−exp⁡(−γπft)r(t)=1- (-γπ ft) (7) where γ is drawn from a log-uniform distribution over [1,10][1,10]. This parameter controls the rise speed of the excitation, ranging from smooth progressive growth to rapid establishment of the underlying waveform. From this procedure, 100 000100\,000 H sequences are generated and fed to the lamination model to form the training dataset. As will be shown in Section 3.8, this dataset size provides sufficient accuracy. Smaller datasets, on the other hand, lead to a degradation in model performance. While the lamination model is evaluated with a fine temporal discretization, as discussed in Section 2.2, the resulting sequences are downsampled to 500 uniform time steps spanning half of the fundamental period T=1/fT=1/f. Although the second half is not strictly redundant — owing to phase differences, pulse modulation, and ramp-up transients — it contains largely similar information. Furthermore, restricting the sequence length speeds up the RNN training while reducing the risk of vanishing or exploding gradients. Validation and test datasets are generated in the same manner, each comprising 10 00010\,000 sequences. Overall, the dataset generation process requires approximately four days on a single CPU core of an HPE ProLiant XL225n compute node from the Lucia cluster. The workload is, however, embarrassingly parallel and can therefore be completed in a much shorter wall-clock time when distributed across multiple CPU cores. 3.3 Permeability-Based Error Scaling A common loss for regression tasks is the Mean Squared Error (MSE). Restraining its definition to the time step k of a single sequence, it is given by ℒMSE=|k−~k|2.L_MSE= |H^k- H^k |^2. (8) Poor results are however obtained using this loss as is. Indeed, in the saturated regime, the magnetic response can vary over tens of kA/m or more, so even a small relative error produces a large absolute contribution to the loss. Conversely, low-field regions (where the material behavior is most sensitive and non-linear) contribute comparatively little to the total error, leading to poor accuracy in this regime. Common strategies to improve accuracy across a wide range of output values include: (i) using a relative loss, e.g. |k−~k|2/(|k|2+ϵ)|H^k- H^k|^2/(|H^k|^2+ ), with an ϵ added in the denominator term to control the trade-off between accuracy at low and high fields; (i) defining the loss on the time derivative, although it may suffer the same issues as the MSE, and may cause the predicted trajectories to drift progressively over time; and (i) applying a logarithmic scaling, which improves overall accuracy but does not fully resolve the imbalance between low- and high-field regions. To overcome these limitations, this paper introduces a loss that conjointly scales the predicted field ~k H^k and the target field kH^k by the anhysteretic magnetic permeability μanh(|k|) _anh(|H^k|). This scaling essentially expresses the error in terms of a magnetic flux density, naturally reducing the dominance of high saturated fields. The resulting loss, referred to as the Mean Squared Scaled Error (MSSE), is specifically designed for applications involving magnetic saturation and writes ℒMSSE=(μanh(|k|)|k−~k|)2,L_MSSE= ( _anh (|H^k| ) |H^k- H^k | )^2, (9) with μanh(|k|)=μ0+ℒ(|k|ha)Ja|k|+ℒ(|k|hb)Jb|k|, _anh (|H^k| )= _0+ L ( |H^k|h_a ) J_a|H^k|+ L ( |H^k|h_b ) J_b|H^k|, (10) which directly stems from (2). Fig. 4 demonstrates the suitability of (9), by comparing the predictions made on two simple sequences after having trained the model of Fig. 3 on the aforementioned losses. Figure 4: Given a low-field (left) and a high-field (right) sinusoidal H curves oscillating at 200 Hz with 500 time steps per period, and the corresponding B sequences returned by the lamination model, comparison between the ~() H(B) sequences returned by models trained to minimize different losses: 1) the MSE, 2) a relative loss, 3) a loss on the time derivative, 4) a loss using a logarithmic scaling, and 5) the proposed MSSE. All models were trained on the same datasets with uniformly spaced time steps and using identical numbers of layers, hidden units, and activation functions, as adopted for the final architecture in Section 3.7. A grid search over the learning rate and batch size was performed for each loss individually. Among the losses, the proposed MSSE loss shows the best performance in both low- and high-field regions. 3.4 Learning the Prediction Confidence Beyond predicting the material response itself, it is valuable for the N to also indicate how much confidence it places in each of its own predictions. Nix and Weigend [44] proposed in 1994 a loss function that enables an N to learn both the mean and standard deviation of a normal probability distribution, originally to capture the uncertainty inherent to data collection. Although this source of uncertainty is absent here – since the lamination model used to generate the training dataset is deterministic (see Sections 2 and 3.2) – a similar loss can be used to make the N predict an expected error. Evaluating the N previously-trained with the MSSE loss on the whole validation dataset shows that the error μanh(|k|)|k−~k| _anh(|H^k|)\,|H^k- H^k| is exponentially distributed. Following the steps presented in [45, Section 5.2], the loss is therefore derived by taking the negative log-likelihood of an exponential distribution with scale parameter ε . This yields the exponential negative log-likelihood (ENLL) loss, adopted as the final training loss and expressed here for a single time step k: ℒENLL=log⁡(εk)+μanh(|k|)|k−~k|εk.L_ENLL= ( ^k )+ _anh(|H^k|)\,|H^k- H^k| ^k. (11) The two terms in the ENLL loss play complementary roles: the logarithmic term penalizes overly large values of εk ^k, whereas the second term penalizes values that underestimate the prediction error. Their balance encourages εk ^k to match the actual prediction error. Consequently, the scale parameter εk ^k acts as a predicted error, expressed in tesla, which varies from sequence to sequence and from time step to time step. At evaluation, εk ^k is expected to grow when the model is given sequences it poorly handled during training. In practice, εk ^k is also expected to grow on sequences lying outside of the training domain, but this behavior remains a tendency rather than a certainty, as the N response to unseen inputs cannot itself be guaranteed. To ensure stable training, εk ^k must remain strictly positive. In practice, a minimum value of εmin=1 _min=1 mT is enforced using a softplus-like activation function. εk=1ϰln⁡[1+exp⁡(ϰ(ςok+1−εmin))]+εmin, ^k= 1 [1+ ( ( o^k+1- _min ) ) ]+ _min, (12) where oko^k denotes the raw output of the decoder, ς=0.1 =0.1 controls the slope of the activation function, and ϰ=5 =5 its sharpness. 3.5 Robustness at High Fields The model is trained on H sequences with maximum amplitudes of roughly 100100 kA/m (see Section 3.2), which, via (10), correspond to anhysteretic amplitudes ||=μanh(||)|||B|= _anh(|H|)\,|H| of approximately 2.052.05 T. The model is therefore expected to provide reliable predictions within this range. Predictions for B amplitudes exceeding 2.052.05 T, however, lie outside the training domain and their accuracy cannot be guaranteed. While exceeding 2.052.05 T is uncommon in practice, excursions to higher field values may occur during the Newton–Raphson iterations of the 2D macroscopic-scale FE simulations (see Section 4). Therefore, the model must remain well-behaved even when |||B| exceeds 2.05 T. Beyond 100100 kA/m, ferromagnetic materials are fully saturated. In this regime, the material response follows a simple anhysteretic law, such as (10). A practical strategy is thus to detect when the model operates outside its training domain, and to replace it with an anhysteretic law. However, this approach may introduce issues near the switching point, as it may create a discontinuity. The proposed solution relies on two points: i) having the neural network predict only the deviation from an anhysteretic curve; and i) clipping the amplitude of the input kB^k, ensuring that the model never receives data outside its training domain. The resulting enhanced model architecture is illustrated in Fig. 5. Figure 5: Enhanced model architecture. i) Residual formulation The predicted field ~k H^k is now written as the sum of an anhysteretic contribution anhkH_anh^k and a hysteretic deviation ~devk H_dev^k, output of the decoder. The term anhkH_anh^k is computed as νanh(|k|)k _anh(|B^k|)\,B^k, where νanh(|k|)=1/μanh(|k|) _anh(|B^k|)=1/ _anh(|B^k|) expresses a magnetic reluctivity, and is implemented as a piecewise cubic Hermite interpolating polynomial (PCHIP). By inverting (10), it is defined on a uniform partition of ||2|B|^2, each interval having a width ΔB2 B^2. The index of the appropriate spline segment is therefore directly given by ⌊(Bx2+By2)/ΔB2⌋ (B_x^2+B_y^2)/ B^2 . If the PCHIP was defined on ΔB B intervals, an additional square root would be required. Alternatively, if it was defined on ΔH H intervals, finding the corresponding spline index would require a lookup operation. For |k|≥3|B^k|≥ 3 T, the PCHIP is replaced by the analytical expression [(|k|−Bs)/μ0+Hs]/|k|[(|B^k|-B_s)/ _0+H_s]/|B^k|, where the parameters BsB_s and HsH_s are chosen to guarantee νanh(|k|) _anh(|B^k|) is continuous at 33 T. i) Input Clipping When the material is fully saturated, k(k)H^k(B^k) follows the anhysteretic law, so that anhk(k)H_anh^k(B^k) alone captures the behavior of ~k(k) H^k(B^k) and ~devk(k) H_dev^k(B^k) remains constant. In this situation, the magnitude of the input kB^k therefore has no impact on the prediction and can be clipped without impacting the model accuracy. Equivalently, the N can be given an approximate magnetic polarization ~k=[1−μ0νanh(|k|)]k J^k=[1- _0\, _anh(|B^k|)]\,B^k as input, since it is naturally bounded by Jmax=Ja+Jb≈1.95J_max=J_a+J_b≈ 1.95 T (see Section 2.1). With a training dataset populated with B sequences whose maximum amplitude reaches Bmax≈2.05B_max≈ 2.05 T, the N is trained on magnetic polarizations up to [1−μ0νanh(Bmax)]Bmax≈1.92[1- _0\, _anh(B_max)]\,B_max≈ 1.92 T, which is slightly smaller, yet sufficiently close to Jmax=1.95J_max=1.95 T. Note that extending the training dataset beyond 100100 kA/m is therefore unnecessary. It would also be counterproductive, as it would increase the cost of data generation and training, despite the saturated regime being fully described by anhkH_anh^k. In addition to ensuring robustness at high fields, the residual formulation offers further advantages related to the learning problem itself. By constraining the network to predict only the hysteretic deviation ~devk H_dev^k from the known anhysteretic baseline anhkH_anh^k, the target function is significantly smoother and smaller in magnitude than the full material response. This simplification reduces the complexity of the function the N must approximate, allowing good accuracy to be achieved with a more compact architecture. Furthermore, since the dominant, well-understood part of the behavior is handled analytically, the N is implicitly guided toward physically meaningful predictions, which tends to reduce the risk of spurious oscillations. 3.6 Sequence Resampling for Adaptive Time-Stepping At each training iteration, a batch of sequences is drawn at random from the 100 000100\,000 sequences composing the training dataset. To encourage accurate predictions across different sequence lengths, these sequences are randomly truncated, resulting in sub-sequences of lengths between 10 and 500 uniformly-spaced time steps. These sub-sequences are then resampled to expose the model to non-uniform temporal discretizations, making the trained model compatible with adaptive time-stepping schemes. The complete resampling algorithm is presented in B. In brief, the algorithm first selects a random subset of points to retain from the source sub-sequence, while enforcing a maximum gap g¯=5 g=5 between consecutive retained points. Since the highest considered harmonic (50th50^th harmonic) is initially sampled by 20 points per period (see Section 3.2), enforcing g¯=5 g=5 guarantees that at least 4 points per period are retained after resampling, keeping it sufficiently resolved. The algorithm then inserts up to ninsn_ins linearly interpolated points according to a skewed random distribution, allowing the interpolated points to be either smoothly spread across the source sub-sequence or concentrated within localized regions. In the present work, ninsn_ins is chosen equal to the length of the source sub-sequence, allowing substantial variation in the time steps while keeping the resampled sequence length comparable to that of the source sub-sequence. This limits training time and mitigates vanishing and exploding gradient issues. A minimum normalized spacing δmin=0.05 _min=0.05 is also enforced between consecutive output points within each segment, corresponding to 5 % of the source-sequence time step Δt t. This minimum spacing prevents the insertion of nearly coincident points that could cause numerical instabilities during training. With this spacing, up to 1/δmin−1=191/ _min-1=19 points can be inserted within a segment. However, the algorithm limits the number of inserted points per segment to m¯=10 m=10, in order to leave sufficient room for the remaining points to be distributed with varying, non-degenerate time steps. Moreover, during training, the loss is evaluated only at the retained source points to prevent biasing the model toward the interpolation scheme. Note that Neural Ordinary Differential Equations (Neural ODEs) could inherently handle variable time steps, as Neural ODEs learn the system dynamics and integrates it over the requested time interval, rather than directly predicting the system output at a given time. In the context of constitutive surrogate models for FE simulations, this approach may however be inappropriate, since it may require multiple N evaluations to perform a single nonlinear iteration, thereby increasing the computational cost. 3.7 Hyperparameter Tuning and Training Fig. 6 compares the training, validation, and test metrics obtained for different architectural variants, and justifies the hereunder selected architecture. Preliminary experiments revealed that using a single GRU cell is sufficient to achieve high quality predictions. Slightly better results are observed by stacking two GRU cells but, since it doubles the evaluation time, a single GRU cell is preferred. This GRU cell uses a hidden size of 300, chosen as a compromise between accuracy and computational cost. The encoder and decoder use respectively two and three layers, as deeper architectures did not improve performance significantly, whereas shallower configurations were found less accurate. The decoder’s final layer uses a linear activation function to produce continuous and unbounded ~devk H_dev^k values, whereas the activation functions used in the encoder and in the first two decoder layers were determined through a grid search. To guarantee a smooth Jacobian ∂/∂ / and thus well-behaved Newton–Raphson iterations, the activation function search was restricted to continuously differentiable candidates: tanh, softplus, and CELU. Although different activation function combinations produced similar final results, some exhibited instabilities during training. The final architecture uses softplus activations in both encoder layers, and CELU and softplus activations in the first and second decoder layers, respectively. Figure 6: Training, validation, and test metrics for different architectural variants, each differing from the selected architecture a) by a single parameter. All models were trained using the Adam optimizer, a learning rate of 2⋅10−42· 10^-4 and a batch size of 128. While the ENLL loss is used for training, the error μanh(||)|−~| _anh(|H|)\,|H- H| is used on the test dataset. For readability, a moving-average operator was applied on the training and validation curves, using a window of 30. The Adam optimizer was employed throughout the optimization procedure, while the learning rate and batch size were progressively adjusted, resulting in a final learning rate of 2⋅10−42· 10^-4 for a batch size of 128. Overall, the N-based model contains approximately 670 000 parameters, corresponding to a storage size of 2.8 MB. It was trained two days on a single NVIDIA A100 GPU from the Lucia cluster using Pytorch [46]. Although training was performed in single precision (float32), inference is carried out in double precision (float64), as the round-off errors inherent to single-precision arithmetic were observed to prevent the FE simulations of Section 4 from reaching sufficiently small residuals. Importantly, no additional training or fine-tuning was performed; all subsequent results reported in this work were obtained using this single trained model. 3.8 Standalone Performance To assess the model performance, the error μanh(||)|−~| _anh(|H|)\,|H- H| is evaluated on the test dataset described in Section 3.2, yielding a low average value of 8.78 mT. This error is defined as the norm of the difference between the target H and predicted ~ H magnetic field vectors, and expressed in T through the scaling factor μanh(||) _anh(|H|), allowing a direct comparison with the predicted error ε . Assuming the error follows an exponential distribution, the expected fractions of predictions whose observed error lies within the 1ε1 , 2ε2 , and 3ε3 predicted error intervals are 63.2 %, 86.5 % and 95 %, respectively. In practice, the observed fractions on the test dataset are 63.5 %, 93.6 %, and 99.2 %. The coverage curve shown on the left of Fig. 7 extends this comparison to all confidence levels by plotting the observed proportion of predictions lying within the corresponding predicted error interval. A perfectly calibrated model would yield a coverage curve lying exactly on the diagonal, where the observed coverage systematically equals the nominal confidence level. Curves below the diagonal indicate overconfidence, whereas curves above the diagonal correspond to conservative error estimates. As shown in Fig. 7, the model is slightly overconfident for small confidence levels, where the predicted error intervals are somewhat too narrow, and slightly conservative at higher confidence levels. The right panel of Fig. 7 also depicts how the mean observed error and mean predicted error vary with the magnetic flux density magnitude |||B|. In particular, both decrease as saturation is reached, reflecting the fact that the magnetic response becomes essentially anhysteretic, so that anhH_anh alone provides an accurate prediction. In the strongly saturated regime, the observed error drops below the imposed lower bound εmin _min, making the predicted error conservative. Figure 7: Comparison between the observed error μanh(||)|−~| _anh(|H|)\,|H- H| and the predicted error ε on the test dataset. Left: Coverage curve comparing the observed and ideal coverages, showing good overall calibration with slight overconfidence at low confidence levels, and slight underconfidence at high confidence levels. Right: Mean observed error and mean predicted error in bins of the magnetic flux density magnitude |||B|. Fig. 8 compares three target H sequences from the test dataset, representative of low-, medium-, and high-field regimes, with the corresponding predicted ~ H sequences. It also reports the observed error |−~||H- H| and the scaled predicted error ε/μanh(||) / _anh(|H|). The figure demonstrates close agreement between the target and predicted sequences across the full range of field amplitudes, including regions where small variations in B produce significantly amplified responses in H. Moreover, ε/μanh(||) / _anh(|H|) provides a good approximation of the error |−~||H- H|. Figure 8: Comparison between three target H sequences from the test dataset, and the corresponding predicted ~ H sequences. The mismatch |−~||H- H| and the predicted error ε/μanh(||) / _anh(|H|) show a good correlation. Similarly, Fig. 9 compares the ~ B field, obtained by inverting the N-based model via Newton-Raphson for a given H field, against the reference lamination solution B. The Jacobian ∂~/∂~∂ B/∂ H required by the scheme is obtained via automatic differentiation of the N-based model, yielding exact derivatives at low computational cost. This formulation also naturally enables adaptive time-stepping strategies. To illustrate this capability, the imposed H field is designed to switch from strongly saturated regions to low-field regions, and an adaptive time-stepping strategy, controlled by the predicted error ε , is employed: the time step Δt t is doubled when ε<25 <25 mT and halved when ε>40 >40 mT, with minimum and maximum time steps of 1.251.25 μ and 4040 μ , respectively. When fully-saturated, the anhysteretic contribution anhH_anh alone governs the model response, resulting in low predicted error ε and allowing the adaptive scheme to increase the time step. Conversely, at lower field levels, the irreversible contribution ~irr H_irr becomes significant and the predicted error ε increases, prompting the scheme to reduce the time step. This adaptive strategy proves effective: the predicted ~ B is in close agreement with the reference lamination solution B. Fig. 9 also demonstrates the ability of the model to extend out of its training range in the saturation regime, as |||B| reaches 3.2 T. Figure 9: Inversion of the N-based surrogate model via Newton-Raphson for a test H field spanning low- and high-field regimes, using an adaptive time-stepping strategy. 4 Macroscopic-Scale Simulation The magnetic vector potential A formulation is used for 2D FE simulations of electrical machines because it reduces to a single scalar unknown (the out of plane component of the vector potential) discretized with nodal Lagrange elements (rather than edge elements), makes div=0div\,B=0 exact with flux lines as its contours, lets prescribed winding currents enter trivially as a source while giving flux linkage and torque directly, and handles the dominant non-conducting regions without any topological cuts [47]. Let Ω denote the computational domain, Ωcore⊂Ω _core⊂ the laminated ferromagnetic core region, Ωs⊂Ω _s⊂ the stranded conductor region and sj_s the imposed current density, it seeks A such that (μ0−1curl,curl′)Ω\Ωcore+(μ−1curl,curl′)Ωcore−(s,′)Ωs=0 ( _0^-1curl\,A,curl\,A )_ _core+\; (μ^-1curl\,A,curl\,A )_ _core-\; (j_s,A )_ _s=0 (13) holds for all test function ′A in the appropriate function space. The coupling with the N-based model is introduced by replacing the constitutive relation μ−1curl=μ−1μ^-1curl\,A=μ^-1B in Ωcore _core by the N-based prediction ~() H(B). Owing to the nonlinearity of ~() H(B), the resulting problem is solved using a Newton-Raphson scheme. At iteration i, the constitutive law is linearized around i−1B_i-1 as ~(i)≈~(i−1)+∂~∂(i−1)(i−i−1), H(B_i)≈ H(B_i-1)+ ∂ H (B_i-1) (B_i-B_i-1 ), (14) where ~(i−1) H(B_i-1) is evaluated in a forward pass of the N, and the Jacobian ∂~/∂(i−1)∂ H/ (B_i-1) is computed by reverse-mode automatic differentiation. Note that the GRU hidden state is kept fixed during the Newton-Raphson iterations and updated only upon convergence. Consequently, once convergence has been reached, the N-based model must be queried one additional time to update its hidden state. Alternatively, two hidden states can be stored per element in Ωcore _core: the input hidden state Sk−1S^k-1, which remains fixed during the Newton-Raphson iterations, and the output hidden state SikS^k_i, returned by the GRU at each iteration, which becomes the new input hidden state SkS^k once convergence is reached. The latter approach is adopted in this work, requiring the storage of 2⋅300=6002· 300=600 double-precision numbers (4.8 kB) per element in Ωcore _core. The proposed methodology is evaluated on two representative case studies, both making use of the same N-based model trained for stacks of M235-35A steel laminations. Subsection 4.1 first considers a three-phase transformer, and subsection 4.2 then considers a synchronous reluctance motor. In both applications, the models are discretized using the free and open-source software Gmsh [48], while the magnetic vector potential field A is discretized with first order elements using the free and open-source FE solver GetDP [49]. All simulations reported in this work also perform element-wise integration using a single Gauss quadrature point, and consider relative and absolute residual tolerances of 10−710^-7 and 10−1010^-10, respectively. The convergence behavior of the Newton-Raphson scheme is discussed in in subsection 4.3, while subsection 4.4 analyses the computational cost. 4.1 Three Phase Transformer For each phase p∈A,B,Cp∈\A,B,C\ of the three-phase transformer, the current density in the primary and secondary windings is imposed as jp(t)=105(1−exp⁡(−2πft))sin⁡(2πft−ϕp)j_p(t)=10^5\;(1- (-2π ft))\; (2π ft- _p), where f=50f=50 Hz is the frequency, ϕA=0∘ _A=0 , ϕB=120∘ _B=120 and ϕC=240∘ _C=240 are phase lags, and the factor (1−exp⁡(−2πft))(1- (-2π ft)) introduces a gradual startup. A homogeneous Dirichlet boundary condition is also imposed on the boundary of Ω . Fig. 10 depicts the mesh of the analyzed three-phase transformer, together with the curves of the imposed current density in each phase and the corresponding By(t)B_y(t) and Hy(t)H_y(t) curves at points a, b, and c located at the center of each leg. These curves result either from an anhysteretic simulation (dotted lines) or from a simulation using the N-based model (solid lines), both running over 800 constant time steps for a single period. The impact of hysteresis is visible on the ByB_y curves which, compared to the anhysteretic case, are delayed due to remanent magnetization. Figure 10: Left: Mesh of the considered three-phase transformer, with circa 5300 nodes and 10 500 elements, of which 4600 belong to Ωcore _core. Right: Curves of the current densities jAj_A, jBj_B and jCj_C imposed in phase A, B and C, and curves of the magnetic flux density B and magnetic field H along the y direction in the center of each leg, resulting from a simulation with the proposed N-based model (solid lines) or with an anhysteretic law (dotted lines). To evaluate the effect of time-step refinement on simulation accuracy, two additional simulations are performed over a single period T=1/fT=1/f using the same trained N-based model. The first considers half less (400) constant time steps, and the second uses an adaptive time step Δtk t^k. At each time step k, the converged solution is accepted if the maximum predicted error maxΩ⁡εk _ ^k is lower than 100100 mT, it is otherwise rejected and the time step is halved. In case maxΩ⁡εk _ ^k is lower than 7070 mT, the time step is doubled back, with a maximum time step of T/100=200μT/100=200~ , and a minimum time step of T/1600=12.5μT/1600=12.5~ , resulting in a simulation of 210 accepted time steps. Note that this adaptive strategy is guided by the predicted error εk ^k rather than by an estimate of the local truncation error, and should therefore be regarded as a heuristic approach. Nevertheless, classical error-controlled time-stepping schemes remain applicable. To measure the error introduced by the N-based model, the H sequences (one per element in Ωcore _core) obtained during the simulation with 800 time steps are injected into the lamination model. The resulting sequences, noted ¯ B, serve as target and can be compared to the B sequences obtained directly during the three simulations using the N-based model. They are used to compute the mismatch |−¯||B- B|, which is a direct measure of the simulation accuracy. If |−¯||B- B| was zero at each Gauss point and at each time step, it would mean the N-based model is a perfect substitute to the lamination model. The same results could therefore be obtained using the lamination model as material law, using sufficiently fine temporal refinement and at a computational cost and complexity largely superior. Figure 11: Spatial distributions of the maximum-in-time mismatch (first column), time-averaged mismatch (second column), time-averaged predicted error (third column), and exceedance (last column), for the three N-based simulations. The first two columns of Fig. 11 are used to assess the simulation accuracy. The first column displays the spatial distribution of the maximum-in-time mismatch maxt⁡|−¯| _t|B- B|, while the second displays the spatial distribution of the time-averaged mismatch ⟨|−¯|⟩t |B- B| _t, both computed element-wise over the ferromagnetic core region Ωcore _core, and for the three aforementioned N-based simulations. For all three simulations, the maximum-in-time mismatch remains below 66 mT throughout Ωcore _core, while the mean mismatch ⟨|−¯|⟩t,Ω |B- B| _t, , averaged over time and space, is about 6 mT. When compared locally with the norm of B, these mismatches correspond to relative values ⟨|−¯|⟩t/⟨||⟩t |B- B| _t/ |B| _t that remain below 5.4 % throughout Ωcore _core, with a spatial average ⟨|−¯|⟩t/⟨||⟩t⟩Ω |B- B| _t/ |B| _t _ of about 0.7 %. Without access to the lamination model, the mismatch |−¯||B- B| cannot be computed, and one should rely exclusively on the predicted error ε . The third column of Fig. 11 shows the spatial distribution of its time average ⟨ε⟩t _t. It is observed that, except in the vicinity of point c, ⟨ε⟩t _t remains larger than the time-averaged mismatch ⟨|−¯|⟩t |B- B| _t. Furthermore, locations where the time-averaged mismatch is larger correspond to locations where the time-averaged predicted error is larger, showing the ability of the N to detect challenging regions. The time-averaged mismatch ⟨|−¯|⟩t |B- B| _t and time-averaged predicted error ⟨ε⟩t _t however provide only aggregate measures over time and do not indicate whether the predicted error consistently bounds the actual mismatch, or if the mismatch frequently exceeds the predicted error. The last column of Fig. 11 reports the percentage of time the mismatch is larger than the predicted error, hereafter referred as the exceedance ⟨(|−¯|>ε)⟩t 1(|B- B|> ) _t. The mean exceedance ⟨(|−¯|>ε)⟩t,Ω 1(|B- B|> ) _t, is obtained by further averaging over space. It reaches 16.4 %, 13.6 % and 7.8 % for respectively the 800-time-step, 400-time-step and variable-time-step simulations, showing the predicted error indeed remains larger than the mismatch most of the time. The mean exceedance decreases to 4.9 %, 4 % and 2.2 % when considering twice the predicted error, and to 2 %, 1.9 % and 0.9 % when considering three times the predicted error. Figure 12: Comparison of the B, ε , |−¯|/ε|B- B|/ , and ()B(H) curves obtained during the N-based simulations with 800 constant time steps, 400 constant time steps, and 210 adaptive time steps, at the points c (time-averaged mismatch greater than time-averaged predicted error), d (largest maximum-in-time mismatch), and e (largest time-averaged predicted error) identified in Fig. 10 and 11. The target ¯ B curves are given by the lamination model, evaluated on the H curves obtained during the 800-time-step N-based simulation. At the top of the figure, the Δt t curves are also depicted, together with the ByB_y curves obtained at the center of each leg of the three-phase transformer (points a, b, and c). While the previous analysis focuses on global and time-averaged quantities, it is also instructive to examine the behavior at specific points. Fig. 12 shows the time evolution of B, ε , and |−¯|/ε|B- B|/ , together with the corresponding ()B(H) curves, at points c, d, and e, identified in Fig. 10 and 11. Point c stands in the middle of the region where the time-averaged mismatch is greater than the time-averaged predicted error, while point d exhibits the worst maximum-in-time mismatch and point e showcases the largest time-averaged predicted error. Primarily focusing on the B and ()B(H) curves, it is observed that the N curves closely match the target curves. The ByB_y curve of point c also shows that the higher time-average mismatch does not originate from a discrepancy in shape, but rather results from a slight amplitude offset between the target and prediction curves. This curve explains why the time-averaged mismatch is slightly larger than the time-averaged predicted error in the vicinity of point c. Despite this deviation, the predicted By(Hy)B_y(H_y) loop still accurately reproduces the reference curve. By comparing the ε and B curves, it appears that ε is minimum when |||B| is sufficiently large to reach saturation, in which case the model prediction mainly relies on the anhysteretic contribution anhH_anh. Fig. 12 also depicts the evolution of the time step for the adaptive-time-step simulation, and the ByB_y curves at the center of each of the transformer leg. The time step is decreased whenever a transformer leg switches from nearly-zero field to saturation or inversely, reflecting a rise in the predicted error during these transitions. Eventually, the |−¯|/ε|B- B|/ curves reveal that the previously reported exceedance translates locally into ratios that stay close to unity for most of the time horizon. 4.2 Synchronous Reluctance Motor The considered synchronous reluctance motor is supplied with three-phase sinusoidal currents of 217.31ARMS217.31~A_RMS, delivering a power of 220220 kW. It has two pole pairs and operates at a frequency of 6565 Hz, resulting in a nominal speed of θn˙=1950 \, _n=1950 rpm. The flux barriers in the rotor result from a shape optimization process aimed at minimizing torque ripple [50]. Fig. 13 presents the considered mesh — which relies on a moving band technique to handle the relative motion between stator and rotor — together with a map of the time-averaged magnetic-flux-density amplitude ⟨||⟩t |B| _t. It also introduces two coordinate systems: a fixed xyxy frame and a rotating x′y′x y frame attached to the rotor. When evaluating the N-based model, the rotor sequences must be expressed according to the x′y′x y frame, in order to eliminate artificial time variations due to rotation. Figure 13: Left: Mesh of the considered synchronous reluctance motor, with circa 14700 nodes and 29600 elements, of which 17000 belong to Ωcore _core. Right: Norm of B, averaged over time. As with the transformer, the target ¯ B sequences (one per element in Ωcore _core) are obtained by first running an N-based simulation, and subsequently feeding the resulting H sequences to the lamination model. This simulation uses 4800 constant time steps to model a single period, a deliberately large number to ensure a sufficiently fine temporal resolution and accurately capture rapid variations due to rotor–stator tooth alignment. Simulations with a coarser temporal resolution may lack precision in the predicted H sequences and, although the impact on subsequent N predictions remains limited, it could drive the EB hysteresis model cells into a spurious state from which they cannot automatically recover, ultimately compromising the accuracy of the target ¯ B sequences. At first, these ¯ B sequences are compared against the B sequences obtained from an N-based simulation with 1200 time steps for a single period. Unlike the simulations of the transformer, the three-phase currents are imposed directly as sinusoidal waveforms at their nominal amplitude, without any ramp-up phase. As a result, the first time step is more challenging to converge, since the solution must evolve from an initially completely demagnetized state along the initial magnetization curve. To facilitate convergence, a relaxation factor of 0.8 is applied for this initial time step, and no strict limit is imposed on its number of iterations. Fig. 14 presents the spatial distributions of the Figure 14: Spatial distributions of the maximum-in-time mismatch (first column), time-averaged mismatch (second column), time-averaged predicted error (third column), and exceedance (last column), for the 1200-time-step N-based simulation. maximum-in-time mismatch maxt⁡|−¯| _t|B- B|, time-averaged mismatch ⟨|−¯|⟩t |B- B| _t, time-averaged predicted error ⟨ε⟩t _t, and exceedance ⟨(|−¯|>ε)⟩t 1(|B- B|> ) _t. The peak mismatch maxt,Ω⁡|−¯| _t, |B- B| reaches 73 mT at point c, while the time-averaged mismatch remains below 27.5 mT throughout Ωcore _core, reaching its maximum at point e. The time-averaged predicted error attains its maximum at the rotor air-gap interface (point f), where the field undergoes the most abrupt variations, ranging from strong saturation when aligned with a stator tooth, to nearly zero otherwise. The predicted error also remains strictly above 1 mT, as discussed in Section 3.4, and, except for a region near point a, consistently bounds the time-averaged mismatch. The mean relative mismatch is here computed as ⟨|−¯|⟩t/(⟨||⟩t+ϵ)⟩Ω |B- B| _t/( |B| _t+ε) _ , where ϵ=1ε=1 mT is introduced to avoid near-zero divisions in regions where ⟨||⟩t |B| _t vanishes (see Fig. 13). The resulting mean relative mismatch is 0.5 %. In Fig. 15 and 16, the B, ε , |−¯|/ε|B- B|/ , and ()B(H) curves are displayed at points a, b and c in the stator, and at points d, e and f in the rotor. Similarly to point c in the transformer, the B and ¯ B curves at point a exhibit a slight amplitude shift, which explains the localized larger error near point a. For all the other curves, B and ¯ B are in close agreement. Figure 15: Comparison of the B, ε , |−¯|/ε|B- B|/ , and ()B(H) curves obtained during the N-based simulation with 1200 time steps, at the stator points a, b and c, identified in Fig. 13 and 14. The target ¯ B curves are given by the lamination model, evaluated on the H curves obtained during a 4800-time-step N-based simulation. Figure 16: Comparison of the B, ε , |−¯|/ε|B- B|/ , and ()B(H) curves obtained during the N-based simulation with 1200 time steps, at the rotor points d, e and f, identified in Fig. 13 and 14. The target ¯ B curves are given by the lamination model, evaluated on the H curves obtained during a 4800-time-step N-based simulation. To assess the sensitivity of the method to reduced temporal resolution, three additional N-based simulations are performed with 600, 300, and 150 nominal time steps per period. Note that, in practice, an adaptive time-stepping strategy is employed: except for the first time step, the time step is halved whenever convergence is not achieved within 20 iterations, and doubled when convergence is reached in fewer than 10 iterations. Table 1 reports the peak mismatch maxt,Ω⁡|−¯| _t, |B- B|, mean mismatch ⟨|−¯|⟩t,Ω |B- B| _t, , and mean predicted error ⟨ε⟩t,Ω _t, at the stator and at the rotor for the aforementioned N-based simulations. It is first worth noting that the simulations maintain good accuracy on average, as the mean stator and rotor mismatches increase only modestly, from 6.1 to 6.9 mT and 1.8 to 2.6 mT, respectively, when reducing the number of nominal time steps from 1200 to 150. Moreover, the mean and peak mismatches remain nearly unchanged between 1200 and 300 nominal time steps, suggesting that a resolution of 300 time steps is already sufficient. This observation is further supported by the behavior of the mean predicted error, which increases only slightly from 1200 to 300 time steps, but rises more markedly from 300 to 150 time steps. Table 1: Peak mismatch, mean mismatch, and mean predicted error at the stator and rotor for the N-based simulations with 1200, 600, 300, and 150 nominal time steps per period. Nominal number Stator [mT] Rotor [mT] of time steps maxt,Ω⁡|−¯| _t, - B ⟨|−¯|⟩t,Ω - B _t, ⟨ε⟩t,Ω _t, maxt,Ω⁡|−¯| _t, - B ⟨|−¯|⟩t,Ω - B _t, ⟨ε⟩t,Ω _t, 1200 73.0 6.1 5.8 47.3 1.8 5.5 600 77.5 6.1 7.0 52.5 1.9 7.0 300 82.6 6.2 11.3 204.9 2.1 11.2 150 320.2 6.9 26.0 701.3 2.6 20.7 As the temporal resolution is reduced, the peak error increases more markedly at the rotor than at the stator, which can be explained by the coarser temporal discretization of the rotor sequences. Indeed, with 48 stator teeth and two pole pairs, the 24th24^th harmonic dominates at the rotor and is sampled by only about 6 time steps in the simulation with 150 nominal time steps per period. Nevertheless, Fig. 17, which shows the B and H fields at point f in the rotor for the 1200, 600, 300, and 150 nominal-time-step simulations, demonstrates the ability of the N-based model to handle such poorly resolved temporal variations. This robustness is particularly valuable during the design phase of electrical machines, where coarse time discretizations are often preferred to reduce computational costs. In this regime, the model strongly benefits from the anhysteretic contribution anhH_anh, which provides a physically consistent base to the N predictions, regardless of the temporal discretization. Fig. 18 further compares these N-based simulations by depicting the ()B(H) curves at point b in the stator. Owing to the contribution of anhH_anh, the N-based model provides perfect predictions in the strongly saturated regime, even for coarse temporal discretizations. At lower field levels, although some noise is observed, even the 150-nominal-time-step simulation remains capable of capturing the minor hysteresis loops. Figure 17: Comparison of the B and H curves obtained during the N-based simulations with 1200, 600, 300 and 150 nominal time steps, at the rotor point f. Figure 18: Comparison of the ()B(H) curves obtained during the N-based simulations with 1200, 600, 300 and 150 nominal time steps, at the stator point b. 4.3 Convergence of the Newton-Raphson scheme Fig. 19 presents the convergence behavior of the Newton-Raphson scheme for the synchronous reluctance motor simulations performed with 1200, 600, 300, and 150 nominal time steps, using either the N-based surrogate model or the anhysteretic law. Except for a limited number of instances in the 150-time-step case, where the N-based simulation reduces its time step after reaching the maximum of 20 Newton iterations, both models exhibit nearly identical convergence behavior. This close superposition indicates that the N-based model preserves the numerical behavior of the anhysteretic model. The initial time step also requires more Newton iterations to converge since, as previously mentioned, the three-phase currents are imposed directly as sinusoidal waveforms at their nominal amplitude. Consequently, the solution must evolve from a completely demagnetized state and stabilize along the first magnetization curve. Figure 19: Convergence history of the nonlinear Newton-Raphson iterations, with the anhysteretic material law and with the N-based model, considering the 1200, 600, 300 and 150 nominal-time-step simulations of the synchronous reluctance machine. 4.4 Computational Cost All simulations were carried out on a MacBook Pro laptop equipped with an Apple M1 Max chip, using GetDP for the FE computations [49] and Python to evaluate either the anhysteretic law or the N-based surrogate model. In practice, the transformer simulations run in 12.7, 6.8, and 4.3 minutes for the 800-, 400-, and adaptive-time-step cases respectively, while the synchronous reluctance motor simulations run in 55.4, 34.2, 20.9 and 14.8 minutes for 1200, 600, 300, and 150 nominal time steps respectively. Tables 2 and 3 provide a detailed breakdown of these timings for the transformer and synchronous reluctance motor simulations, decomposed into three contributions: evaluation of the material law in Python, assembly and solution of the system in GetDP, and input/output (I/O) operations corresponding to data exchange between Python and GetDP. This exchange is currently implemented via text files written to and read from disk, which, while functional, accounts for approximately 20 % of the total simulation time. In future work, this file-based communication will be replaced by shared-memory communication to remove this overhead. Table 2: Comparison of computation times and number of system solves for the transformer simulations. “Total time steps” denotes all attempted time steps (accepted and rejected), whereas “effective time steps” counts only the accepted ones; the two coincide for fixed time-step simulations. Time steps System solves Time [s] Time/system solve [s] Total Effective Total /Step Python GetDP I/O Total Python GetDP N-b. 800 800 3631 4.54 292 311 158 762 0.08 0.09 400 400 1944 4.86 157 167 83 409 0.08 0.09 240 210 1243 5.18 93 109 53 256 0.08 0.09 Anh. 800 800 3321 4.15 2 282 138 424 / 0.09 400 400 1751 4.38 1 147 72 221 / 0.08 Table 3: Comparison of computation times and number of system solves for the synchronous-reluctance-motor simulations. “Nominal time steps” refers to the target number of time steps per period, whereas “effective time steps” denotes the actual number of accepted time steps after adaptive refinement; the two differ only when additional intermediate time steps are introduced following a failed convergence attempt. Time steps System solves Time [s] Time/system solve [s] Nominal Effective Total /Step Python GetDP I/O Total Python GetDP N-based 1200 1200 4860 4.05 1376 1180 766 3323 0.28 0.24 600 600 3007 5.01 844 738 467 2050 0.28 0.25 300 300 1829 6.10 517 452 282 1252 0.28 0.25 150 153 1290 8.43 365 323 197 886 0.28 0.25 Anhyst. 1200 1200 4912 4.09 6 1171 739 1917 / 0.24 600 600 3001 5.00 3 726 452 1182 / 0.24 300 300 1849 6.16 2 450 277 730 / 0.24 150 150 1259 8.39 1 310 189 501 / 0.25 For both applications, the time required to evaluate the N-based model in Python is approximately that needed to assemble and solve the system in GetDP. Compared to the anhysteretic simulations, for which the time spent in Python is negligible, and assuming that the I/O overhead would be removed through improved communication, a simulation using the N-based model is therefore limited to roughly twice the cost of an anhysteretic simulation. Beyond these absolute timings, the computational cost scales linearly with problem size. The number of elements in Ω increases by a factor of 2.8 from the transformer to the synchronous reluctance motor (10500 to 29600), and the GetDP assembly and solve time per system solve grows by a consistent factor of 2.8 (0.09 to 0.25 s). The same holds for the Python evaluation, where a 3.7-fold increase in Ωcore _core elements (4600 to 17000) yields a 3.7-fold increase in evaluation time (0.08 to 0.28 s). 5 Conclusions and Perspectives This paper has discussed the 2D FE modeling of electrical machines with laminated ferromagnetic cores. The results show that homogenization of the laminated ferromagnetic core behavior via a surrogate N-based model offers key implementation and numerical advantages over direct computational homogenization: i) The automatic differentiation capabilities of NNs allow an exact evaluation of the N-based ferromagnetic model Jacobian, eliminating the need for finite differences. Furthermore, since the proposed N is smooth by construction, owing to its differentiable activation functions, the resulting Jacobian matrices are smooth, yielding Newton-Raphson convergence comparable to that obtained with conventional anhysteretic constitutive laws. In contrast, Jacobians obtained using direct homogenization would inherit the numerical noise of finite differences and could be affected by the non-smoothness of the EB hysteresis model. i) The N allows coupling an H-conform mesoscopic-scale model (to use the fast vector-play EB approximation) with an efficient B-conform macroscopic machine model. i) The N-based model inherits the accuracy of the fine temporal discretization used during dataset generation and largely preserves this accuracy at inference, even when coarser time steps are used at the macroscopic scale. Furthermore, after training on variable time steps, the N-based model seamlessly handles variable time increments, so that macroscopic-scale adaptive time-stepping strategies are naturally supported. Achieving the same with a direct homogenization approach would require maintaining a fine temporal discretization at the mesoscopic scale throughout the simulation, thereby incurring the full computational cost of fine mesoscopic-scale solves [51] iv) From the perspective of the FE solver, the proposed constitutive model behaves as a single function call. The magnetic history required to describe hysteresis is entirely encoded in the N hidden state, eliminating the need to explicitly manage the internal variables of an underlying mesoscopic-scale model. v) Once trained for a given lamination (material grade and sheet thickness), the proposed surrogate enables computationally efficient FE simulations. Although extending the approach to a new lamination requires generating a new dataset (approximately 4 days on a single CPU core, but embarrassingly parallel) and retraining the network (approximately 2 days), the model can subsequently be used in a wide range of FE simulations, with computational times roughly twice those of conventional anhysteretic simulations. Moreover, the predicted error returned by the N-based model provides an immediate accuracy measure of the FE simulations. Several directions are left for future work. The present study considers sinusoidal current excitations, whereas motors are often operated using pulse-width modulation (PWM); assessing the performance of the surrogate model under such conditions is a natural next step. The N architecture could also be extended to predict ferromagnetic core losses and to account for temperature effects, following an approach similar to [41], where N-based models are used for AC loss prediction in superconducting conductors. On a broader scale, the method currently assumes isotropic magnetic properties; extending it to anisotropic laminations, such as grain-oriented electrical steels, would require adapting both the dataset generation and the network architecture. A 3D extension would also be needed to accurately capture eddy-current paths near edges and geometric discontinuities. Appendix A Parameters of the Energy-Based Model ha=18.18h_a=18.18 A/m hb=3905.7h_b=3905.7 A/m Ja=1.387J_a=1.387 T Jb=0.559J_b=0.559 T ρ=690nΩmρ=690~n \,m d=350μd=350~ N=11N=11 elementary cells: κn[A/m]κ^n~[A/m] 0.0 7.3 18.8 32.1 45.5 55.8 66.9 80.6 99.1 143.0 213.5 wn[−]w^n~[-] 0.075 0.103 0.106 0.342 0.119 0.105 0.053 0.043 0.028 0.019 0.006 Appendix B Resampling algorithm The resampling algorithm is designed to produce sequences with irregular time steps of varying density. It first thins the source sub-sequence by retaining only a random subset of points, then re-enriches it by inserting linearly interpolated points, whose distribution can range from spatially uniform to highly localized. Inputs: srcB_src: source sub-sequence nsrcn_src: length of the source sub-sequence ninsn_ins: number of points to insert g¯ g: maximum gap between two consecutive kept source points m¯ m: maximum number of points inserted per segment δmin _min: minimum spacing between consecutive output points Output: upB_up: upsampled sequence ⊳ Step 1: Determine which source points are kept idxkeep[i]←False,∀i=1,…,nsrc−1idx_keep[i] ,\;∀ i=1,…,n_src-1 i←1i← 1 while i<nsrci<n_src do⊳ point nsrcn_src is always kept, see Step 3 idxkeep[i]←Trueidx_keep[i] Draw g∼Uniform1,…,g¯g \1,…, g\ i←i+gi← i+g end while ⊳ Step 2: Determine how many points are inserted between two source points Draw s∼Uniform[0,0.99]s [0,0.99] ⊳ skewness parameter Draw r[j]∼Uniform(0,1),∀j=1,…,nins+1r[j] (0,1),\;∀ j=1,…,n_ins+1 r[j]←ϵ∀j s.t. r[j]<sr[j]←ε\;∀\,j s.t. r[j]<s ⊳ large s forces many r[j]r[j] near zero (ϵ⪆0ε 0) r[j]←(∑l=1jr[l])/(∑l=1nins+1r[l]),∀j=1,…,ninsr[j]←( _l=1^jr[l])/( _l=1^n_ins+1r[l]),\;∀\,j=1,…,n_ins ⊳ low s: r grows smoothly; high s: r grows in steps idxins[j]←⌊r[j]⋅(nsrc−1)⌋+1,∀j=1,…,ninsidx_ins[j]← r[j]·(n_src-1) +1,\;∀\,j=1,…,n_ins mins[i]←#j:idxins[j]=i,∀i=1,…,nsrc−1m_ins[i]←\#\j:idx_ins[j]=i\,\;∀\,i=1,…,n_src-1 ⊳ number of inserted points in segment i mins[i]←m¯∀i s.t. mins[i]>m¯m_ins[i]← m\;∀\,i s.t. m_ins[i]> m⊳ cap insertions per segment ⊳ Step 3: Build upsampled (output) sequence k←1k← 1 for i=1,…,nsrc−1i=1,…,n_src-1 do if idxkeep[i]idx_keep[i] then up[k]←src[i]B_up[k] _src[i], k←k+1k← k+1 ⊳ append source point i if kept Draw δ[j]∼Uniform(0,1),∀j=1,…,mins[i]+1δ[j] (0,1),\;∀ j=1,…,m_ins[i]+1 δ[j]←∑l=1jδ[l],∀j=1,…,mins[i]+1δ[j]← _l=1^jδ[l],\;∀ j=1,…,m_ins[i]+1 δ[j]←δ[j]⋅1−(mins[i]+1)δminδ[mins[i]+1]+j⋅δmin,∀j=1,…,mins[i]δ[j]←δ[j]· 1-(m_ins[i]+1) _minδ[m_ins[i]+1]+j· _min,\;∀ j=1,…,m_ins[i] ⊳ enforce minimum spacing up[k+j−1]←(1−δ[j])⋅src[i]+δ[j]⋅src[i+1],∀j=1,…,mins[i]B_up[k+j-1]←(1-δ[j])·B_src[i]+δ[j]·B_src[i+1],\;∀ j=1,…,m_ins[i] k←k+mins[i]k← k+m_ins[i] end for up[k]←src[nsrc]B_up[k] _src[n_src]⊳ last source point is always kept Acknowledgments Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant no2.5020.11 and by the Walloon Region. The present research also benefited from computational resources made available on Lucia, the Tier-1 supercomputer of the Walloon Region, infrastructure funded by the Walloon Region under the grant agreement no1910247. Funding sources This work was supported in part by the Walloon Region of Belgium under grant 8139 M&SSCoT and 9001 ELECTRON. The work of Louis Denis was supported by the F.R.S-FNRS. Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work, the authors used ChatGPT in order to revise portions of the manuscript and improve its language and readability. After using this service, the authors reviewed and edited the content as needed and take full responsibility for the content of the published article. References [1] D. Ceylan, R. Zeinali, B. Daniels, K. O. Boynov, and E. A. Lomonova, “Magnetodynamic finite element analysis coupled with a vector hysteresis model applied to a variable flux reluctance machine,” IEEE Transactions on Industry Applications, vol. 60, no. 1, p. 80–92, 2023. doi: 10.1109/TIA.2023.3292806 [2] M. M. Billah, F. Martin and A. Belahcen, “A computationally effective method for iron loss estimation in a synchronous machine from a static field solution,” in International Conference on Electrical Machines, p. 751–757, 2020. doi: 10.1109/ICEM49940.2020.9271020 [3] I. D. Mayergoyz and G. Friedman, “Generalized Preisach model of hysteresis,” IEEE Transactions on Magnetics, vol. 24, no. 1, p. 212–217, 1988. doi: 10.1109/20.43892 [4] D.C. Jiles and D.L. Atherton, “Theory of ferromagnetic hysteresis,” Journal of Magnetism and Magnetic Materials, vol. 61, no. 1-2, p. 48–60, 1986. doi: 10.1016/0304-8853(86)90066-1 [5] A. Benabou, S. Clénet and F. Piriou, “Comparison of Preisach and Jiles–Atherton models to take into account hysteresis phenomenon for finite element analysis,” Journal of Magnetism and Magnetic Materials, vol. 261, no. 1-2, p. 139–160, 2003. doi: 10.1016/S0304-8853(02)01463-4 [6] W. Li and C. S. Koh, “Investigation of the vector Jiles–Atherton model and the fixed point method combined technique for time-periodic magnetic problems,” IEEE Transactions on Magnetics, vol. 51, no. 4, p. 1–6, 2015. doi: 10.1109/TMAG.2014.2360150 [7] M. Nierla, M. Loeffler, M. Kaltenbacher and S. J. Rupitsch, “Comparison of different vector Preisach models for the simulation of ferromagnetic materials,” COMPEL, vol. 38, no. 5, p. 1696–1710, 2019. doi: 10.1108/COMPEL-12-2018-0494 [8] A. H. S. Atyia and A. M. Ghanim, “Limitations of Jiles–Atherton models to study the effect of hysteresis in electrical steels under different excitation regimes,” COMPEL, vol. 43, no. 1, p. 66–79, 2024. doi: 10.1108/COMPEL-02-2023-0061 [9] S. Hussain and D. A. Lowther, “An efficient implementation of the classical preisach model,” IEEE Transactions on Magnetics, vol. 54, no. 3, p. 1–4, 2018. doi: 10.1109/TMAG.2017.2748100 [10] F. Henrotte, A. Nicolet and K. Hameyer, “An energy‐based vector hysteresis model for ferromagnetic materials,” COMPEL, vol. 25, no. 1, p. 71–80, 2006. doi: 10.1108/03321640610634344 [11] K. Jacques, “Energy-based magnetic hysteresis models - Theoretical development and finite element formulations,” Ph.D. dissertation, University of Liège, Liège, Belgium, 2018. [12] A. Bergqvist, “Magnetic vector hysteresis model with dry friction-like pinning,” Physica B: Condensed Matter, vol. 233, no. 4, p. 342–347, 1997. doi: 10.1016/S0921-4526(97)00319-0 [13] H. Egger, F. Engertsberger, L. Domenig, K. Roppert and M. Kaltenbacher, “On forward and inverse energy-based magnetic vector hysteresis operators,” IEEE Transactions on Magnetics, vol. 61, no. 4, p. 1–7, 2025. doi: 10.1109/TMAG.2025.3544507 [14] H. Egger, F. Engertsberger and A. Schafelner, “Efficient evaluation of forward and inverse energy-based magnetic hysteresis operators,” IEEE Transactions on Magnetics, 2025. doi: 10.1109/TMAG.2025.3621134 [15] V. François-Lavet, F. Henrotte, L. Stainier, L. Noels and C. Geuzaine, “An energy-based variational model of ferromagnetic hysteresis for finite element computations,” Journal of Computational and Applied Mathematics, vol. 246, p. 243–250, 2013. doi: 10.1016/j.cam.2012.06.007 [16] K. Roppert, M. Kaltenbacher, L. Domenig and L. Daniel, “Magnetoelastic vector hysteresis modeling for electromagnetic devices: A combination of a multiscale model with the energy-based hysteresis framework,” IEEE Transactions on Magnetics, vol. 61, no. 8, p. 1–14, 2025. doi: 10.1109/TMAG.2025.3584819 [17] K. Jacques, R. V. Sabariego, C. Geuzaine and J. Gyselinck, “Inclusion of a direct and inverse energy-consistent hysteresis model in dual magnetostatic finite-element formulations,” IEEE Transactions on Magnetics, vol. 52, no. 3, p. 1–4, 2016. doi: 10.1109/TMAG.2015.2490578 [18] W. E and B. Engquist, “The heterogeneous multi-scale method for homogenization problems,” in Multiscale Methods in Science and Engineering, p. 89–110, 2005. doi: 10.1007/3-540-26444-2_4 [19] I. Niyonzima, R. V. Sabariego, P. Dular, F. Henrotte and C. Geuzaine, “Computational homogenization for laminated ferromagnetic cores in magnetodynamics,” IEEE Transactions on Magnetics, vol. 49, no. 5, p. 2049–2052, 2013. doi: 10.1109/TMAG.2012.2237546 [20] J. Schröder, “A numerical two-scale homogenization scheme: The FE2-method,” in Plasticity and beyond: Microstructures, crystal-plasticity and phase transitions, J. Schröder and K. Hackl, Eds. Springer, 2014, p. 1–64. doi: 10.1007/978-3-7091-1625-8_1 [21] L. Solimene, C. S. Ragusa, A. Giuffrida, N. Lombardo, F. Marmello, S. Morra and M. Pasquale, “A hybrid data-driven approach in magnetic core loss modeling for power electronics applications,” in Mediterranean Electrotechnical Conference, p. 1356–1361, 2024. doi: 10.1109/MELECON56669.2024.10608787 [22] M. d’Aquino, S. Perna, C. Serpico and C. Visone, “Hysteresis compensation by deep learning algorithms,” Physica B: Condensed Matter, vol. 675, art. 415596, 2024. doi: 10.1016/j.physb.2023.415596 [23] H. Zhang, Q. Yang, C. Zhang, Y. Li and Y. Chen, “Magnetic properties simulation of electrical steel sheet based on recurrent neural network,” IEEE Transactions on Magnetics, vol. 60, no. 3, p. 1–4, 2024. doi: 10.1109/TMAG.2023.3316753 [24] Z. Guo, B. H. Nguyen, H. Hamzehbahmani and R. V. Sabariego, “Dynamic hysteresis model of grain-oriented ferromagnetic material using neural operators,” IEEE Transactions on Magnetics, vol. 61, no. 10, p. 1–7, 2025. doi: 10.1109/TMAG.2025.3600089 [25] V. Bertolini, M. Quercio, L. Sabino, M. Stella, A. Faba, F. R. Fulginei, F. Crescimbini and E. Cardelli, “Development of a neural network approach to model magnetic core power loss in high-frequency transformers,” IEEE Transactions on Industry Applications, 2026. doi: 10.1109/TIA.2026.3675154 [26] D. Serrano, H. Li, T. Guillod, S. Wang, M. Luo, C. R. Sullivan and M. Chen, “Neural network as datasheet: Modeling B-H loops of power magnetics with sequence-to-sequence LSTM encoder-decoder architecture,” in Workshop on Control and Modeling for Power Electronics, p. 1–8, 2022. doi: 10.1109/COMPEL53829.2022.9829998 [27] H. Li, D. Serrano, S. Wang, T. Guillod, M. Luo and M. Chen, “Predicting the B-H loops of power magnetics with transformer-based encoder-projector-decoder neural network architecture,” in Applied Power Electronics Conference and Exposition, p. 1543–1550, 2023. doi: 10.1109/APEC43580.2023.10131497 [28] M. Chen, H. Li, S. Wang, T. Guillod, D. Serrano, N. Förster, et al., “MagNet challenge for data-driven power magnetics modeling,” IEEE Open Journal of Power Electronics, vol. 6, p. 883–898, 2024. doi: 10.1109/OJPEL.2024.3469916 [29] S. Wang, H. Kwon, H. Li, Y. Elasser, G. -G. Kang, D. Zhou, D. Grigoryan and M. Chen, “MagNetX: Foundation neural network models for simulating power magnetics in transient,” in Applied Power Electronics Conference and Exposition, p. 2438–2445, 2025. doi: 10.1109/APEC48143.2025.10977420 [30] L. Wu, V. D. Nguyen, N. G. Kilingar and L. Noels, “A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths,” Computer Methods in Applied Mechanics and Engineering, vol. 369, art. 113234, 2020. doi: 10.1016/j.cma.2020.113234 [31] F. Purnode, F. Henrotte, G. Louppe and C. Geuzaine, “Neural network-based simulation of fields and losses in electrical machines with ferromagnetic laminated cores,” International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 37, no. 2, art. e3226, 2024. doi: 10.1002/jnm.3226 [32] A. G. Baydin, B. A. Pearlmutter, A. A. Radul and J. M. Siskind, “Automatic differentiation in machine learning: a survey,” Journal of Machine Learning Research, vol. 18, no. 153, p. 1–43, 2018. url: https://w.jmlr.org/papers/v18/17-468.html [33] A. Chandra, B. Daniels, M. Curti, K. Tiels and E. A. Lomonova, “Magnetic hysteresis modeling with neural operators,” IEEE Transactions on Magnetics, vol. 61, no. 1, p. 1–11, 2025. doi: 10.1109/TMAG.2024.3496695 [34] S. Quondam-Antonio, F. Riganti-Fulginei, A. Laudani, G. M. Lozito and R. Scorretti, “Deep neural networks for the efficient simulation of macro-scale hysteresis processes with generic excitation waveforms,” Engineering Applications of Artificial Intelligence, vol. 121, art. 105940, 2023. doi: 10.1016/j.engappai.2023.105940 [35] M. Stella, A. Faba, V. Bertolini, F. R. Fulginei, L. Sabino, H. Tiismus, A. Kallaste and E. Cardelli, “Modelling of magnetization processes of 3D-printed Fe-Si components by means of an artificial neural network implemented in a FEM scheme,” in International Conference on Electrical Machines, p. 1–7, 2024. doi: 10.1109/ICEM60801.2024.10700301 [36] N. Vuokila, C. Cunning, J. Zhang, N. Akel, A. Khan and D. A. Lowther, “The application of neural networks to the modeling of magnetic hysteresis,” IEEE Transactions on Magnetics, vol. 60, no. 3, p. 1–4, 2024. doi: 10.1109/TMAG.2023.3326708 [37] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser and I. Polosukhin, “Attention is all you need,” Advances in Neural Information Processing Systems, vol. 30, 2017. doi: 10.48550/arXiv.1706.03762 [38] A. Khan, S. Ceesay, Y. P. Teng, R. Wang, H. Yue and D. Lowther, “Generalizable deep neural network based multi-material hysteresis modeling,” in Conference on Electromagnetic Field Computation, p. 1–4, 2022. doi: 10.1109/CEFC-LONG57069.2022.10107520 [39] H. Li, D. Serrano, T. Guillod, S. Wang, E. Dogariu, A. Nadler, M. Luo, V. Bansal, N. K. Jha, Y. Chen, C. R. Sullivan and M. Chen, “How MagNet: Machine learning framework for modeling power magnetic material characteristics,” IEEE Transactions on Power Electronics, vol. 38, no. 12, p. 15829–15853, 2023. doi: 10.1109/TPEL.2023.3309232 [40] H. Li, B. Zhang, Y. Shen and K. Liu, “Research on prediction of nanocrystalline alloy hysteresis properties based on long short-term memory network,” Scientific Reports, vol. 15, no. 1, 2025. doi: 10.1038/s41598-025-91138-1 [41] L. Denis, J. Dular, V. Nuttens, M. Wozniak, B. Vanderheyden and C. Geuzaine, “Acceleration of multi-scale LTS magnet simulations with neural network surrogate models,” IEEE Transactions on Applied Superconductivity, vol. 36, no. 5, p. 1–6, 2026. doi: 10.1109/TASC.2026.3660620 [42] S. Hochreiter and J. Schmidhuber, “Long Short-Term Memory,” Neural Computation, vol. 9, no. 8, p. 1735–1780, 1997. doi: 10.1162/neco.1997.9.8.1735 [43] K. Cho, B. van Merriënboer, C. Gülçehre, D. Bahdanau, F. Bougares, H. Schwenk, and Y. Bengio, “Learning phrase representations using RNN encoder-decoder for statistical machine translation,” in Conference on Empirical Methods in Natural Language Processing, p. 1724–1734, 2014. doi: 10.3115/v1/D14-1179 [44] D. A. Nix and A. S. Weigend, “Estimating the mean and variance of the target probability distribution,” in International Conference on Neural Networks, p. 55–60, 1994. doi: 10.1109/ICNN.1994.374138 [45] S. J. Prince, Understanding deep learning, MIT Press, 2023. [46] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, et al., “PyTorch: An imperative style, high-performance deep learning library,” in Advances in Neural Information Processing Systems, vol. 32, 2019. [47] G. Meunier (ed.), The finite element method for electromagnetic modeling, 2010 [48] C. Geuzaine and J. F. Remacle, “Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, p. 1309–1331, 2009. doi: 10.1002/nme.2579 [49] P. Dular, C. Geuzaine, F. Henrotte and W. Legros, “A general environment for the treatment of discrete problems and its application to the finite element method,” IEEE Transactions on Magnetics, vol. 34, no. 5, p. 3395–3398, 1998. doi: 10.1109/20.717799 [50] C. De Greef, V. Kluyskens, F. Henrotte, C. Versele, C. Geuzaine and B. Dehez, “Time-efficient multi-physics optimization approaches for the design of synchronous reluctance motors,” in Energy Conversion Congress and Exposition, p. 4197–4204, 2021. doi: 10.1109/ECCE47101.2021.9595977 [51] I. Niyonzima and C. Geuzaine and S. Schöps, “Waveform relaxation for the computational homogenization of multiscale magnetoquasistatic problems,” Journal of Computational Physics, vol. 327, p. 416–433, 2016. doi: 10.1016/j.jcp.2016.09.011

Instant research discovery