A hybrid twin based on machine learning enhanced reduced order model for real-time simulation of magnetic bearings

The simulation of magnetic bearings involves highly non-linear physics, with high dependency on the input variation. Moreover, such a simulation is time consuming and can’t run, within realistic computation time for control purposes, when using classical computation methods. On the other hand, classical model reduction techniques fail to achieve the required precision within the allowed computation window. To address this complexity, this work proposes a combination of physics-based computing methods, model reduction techniques and machine learning algorithms, to tackle the requirements. The physical model used to represent the magnetic bearing is the classical Cauer Ladder Network method, while the model reduction technique is applied on the error of the physical model’s solution. Later on, in the latent space a machine learning algorithm is used to predict the evolution of the correction in the latent space. The results show an improvement of the solution without scarifying the computation time. The solution is computed in almost real-time (few milliseconds), and compared to the ﬁnite element reference solution.


Introduction
Recent development and democratization of the available toolbox in machine learning techniques made their use increasingly popular.Thus, they become increasingly popular in engineering applications [1][2][3][4].However, in multiple engineering systems and processes, data is expensive and hard to generate.Therefore, the ability to solely rely on data-driven techniques to generate representative models is hindered.Thus, many recent efforts aim to combine the basic principles and physical knowledge with data-driven methods [5][6][7][8][9].Some of these works aim to use thermodynamics principles [10][11][12][13], while others tend to improve or correct a simulation by infusing additional information coming from data [14][15][16].The result is the so called hybrid twin, or physics-aware digital twin, a combination of scientific knowledge, enhanced with machine learning techniques [17][18][19][20].Gaussian processes are also used to model data and include experimental variability into physical models [21,22].
In this work, we propose a novel digital twin building technique, based on the combination of a spectral approximation, similar to the one used in the Proper Orthogonal Decomposition method, and a time-dependent integrator based on machine learning algorithms.Similar methods already tried to combine the use of spectral decomposition with machine learning.For example, in [23], the authors built a machine learning method in the spectral domain, by aiming to change the eigenvalues and eigenvectors of the transfer operator in the direct space.This approach manipulates a matrix operator, in its eigenvalues/eigenvectors space; to adapt the solution and fit it to the measurements.In [24], the author uses a spectral decomposition of the input image before running a convolution neural network, to achieve a reduction in the memory requirements.The authors in [25,26] propose a deep learning algorithm to map or transform the reduced basis, allowing them to adapt their reduced basis to novel use cases.POD-Galerkin PINN ROM was proposed in [27] to tackle the Navier-Stokes problem.The method proposed here is different from the aforementioned ones available in the literature.The proposed method uses a weight prediction based on the basis decomposition [28].This work however operates on the error of the quantities of interest and not the full fields.The spectral neural network takes as input the errors of the physical quantities of interest, as well as the inputs of the model, to built a correction, based on the re-weighing of the correction's reduced basis, built a priori.The correction basis is build using the POD method [29].The aim of this work is to propose a hybrid modeling enhancing the available state of the art industrial simulation while conserving the same computation time.
The selected showcased industrial application in this work corresponds to the case of magnetic bearings.Magnetic bearings are used since several years to levitate heavy rotating shafts in critical machines.These bearings have the ability to perform quasi indefinitely since no contamination, no lubrication, no contact, and no friction are exhibited [30].Moreover, losses by friction are totally suppressed, while rolling resistance is tremendously reduced, leading to potentially very high rotation speed [31].Despite all these advantages, active magnetic bearings are only used on heavy equipment as they require external energy source to generate the levitation through magnetic fields.Moreover, the control of active magnetic bearings requires a feedback every few microseconds to achieve effective shaft positioning and active suppression of vibrations.The possibility of enhancing this control through simulation-based control is currently being investigated.However, such a control would require a simulation running within few microseconds.Multiple previous efforts aim to model the magnetic bearings using either analytical methods, numerical methods or data driven techniques [32][33][34].Thus, this work proposes a novel technique coupling model reduction, physics based methods and advanced machine learning algorithms to tackle the simulation time requirement.
When it comes to the simulation of magnetic fields, especially in the context of rotors and stators, model reduction was used on several cases to create fast simulations [35].However, to the best knowledge of the authors, there has been no hybrid modeling, combining scientific knowledge to data-driven techniques, in the context of creating a digital twin for active magnetic bearings simulation in real-time.
The work starts in Sect.'Review of the used physical modeling' by reviewing the used physical modeling of the process based on the Cauer Ladder Network method [36][37][38].The Cauer Ladder Network method replaces a full physical-space finite element problem with a sequence of resistors and self inductance, transforming the simulation into a reduced order model problem, based on an equivalent electric circuit model to solve.Section 'Hybrid twin construction' introduces the data-driven method and the coupling between this datadriven technique and the physical model.Section 'Application to the magnetic bearings modelling's quantities of interest' shows the selected application, the results, and discusses the strong and weak points of the method.The work uses the finite element solution as the ground truth because of missing experimental measurements, the same algorithm can be used to built hybrid modeling with experimental results as ground truth when available.The work ends with some conclusions in Sect.'Discussion and conclusions'.

Review of the used physical modeling
In this section we review the physical modeling of a magnetic bearing.The studied problem represents a single coil toy-model used to showcase the numerical technique.The schematic drawing of the coil is shown in Fig. 1.We can see in this figure the coil, which is used to create a magnetic field, leading to the levitation of the steel shaft on top of the coil.A small air gap is to be modeled with high precision to achieve a high fidelity solution.

Maxwell's equations
The problem to solve is the nonlinear magneto-quasistatic problem, a realistic assumption in the framework of magnetic bearings for industrial application, which simplifies much the full order model computation.In these settings, we consider the Maxwell's equations while neglecting the displacement current, i.e.: where B is the magnetic flux density, H the magnetic field, J ind the eddy current density, N a unit current density vector and i is the current flowing through the coil.The constitutive laws give us that: Note that the electric conductivity σ = 0 only in the steel shaft, because we assume that the stator is laminated.
To solve this problem, we use a A * -formulation, then a magnetic vector potential A is introduced such that: In order to impose a voltage in the electrical coil winding, we consider the circuit equation: where R is the resistance of the coil and φ the magnetic flux.The complete problem is then written as: U being the imposed electric tension of the magnetic bearing's coil.As we are addressing a 2D in-plane problem for the sake of simplicity, the potential A is reduced to a scalar representing only the out of plane component of A. At the boundary of the domain, we impose a homogeneous Dirichlet boundary condition.
The electromagnetic force is computed using the Maxwell's tensor Eq. [39], i.e.: where the line domain Force (included in the air) encloses the moving part, n is the normal to this line.Finally, μ 0 is the magnetic permeability of the empty space.

Constructing a reference solution using finite element
In the first part of this work, we use the finite element method (FE) to solve the Eq. ( 9), which will lead us to a reference solution.The reference solution is used later on as our target solution.We can see in Fig. 2 the adaptive mesh used for the FE discretization, which is refined near the airgap and in the conductive part to take into account the eddy currents.
The FE solution is computed in the 2D domain in a square of dimensions 700 × 700 mm, with an air box of 250 mm around the simulated magnetic bearing.
The transient finite element solution is computed for 14 different time variations of the input electric potential U , using 14 different control use cases.We know that PWM exacerbates the eddy current, then to have some example with many eddy currents we implement a simple p-controller with saturation to control the voltage and get as a set point the magnetic flux.
The available use-cases are illustrated in table 1.One of the finite elements solutions is shown in Fig. 3, for illustrative purposes.In 3(a) we illustrate the out of plane component of the magnetic vector potential A z , which, for sake of simplicity, will be referred as A for the rest of this work.Moreover, in 3(b) we illustrate the magnitude of the magnetic flux density B .

A fast physics-based modeled solution
Recent developments in the simulation of magnetic fields generated by coils have lead to the creation of a novel simulation technique based on the use of spectral, mode shape based, representation.This method is currently known by the name of Cauer Ladder Network Method [36][37][38].The method consists in assuming that A can be written in a separated form using: c q a q (x, y) (11) where n q is the number of selected layers in the ladder Network representation, c q is a constant to compute, equal to the current in the equivalent circuit, and a q are the orthogonal basis representation vectors.In this work, the presented results are shown for n q = 8, as beyond n q = 8 it is observed that the results do not improve anymore in the selected use case.To compute a q and c q we use the matrix formulation method presented in [38].The method produces an electric circuit of several (R, L) branches, where R is a resistance and L denotes an inductance.The current intensity in each branch of the ladder Network is noted c j and will be used later as input for the correction of the results.An example of the equivalent Cauer ladder network is shown in Fig. 4. We note: where A REF corresponds to the solution of a magnetostatic problem with i = 15 (that corresponds to the first equation of ( 9) with σ = 0).In fact, this term is used to construct the Cauer Ladder Network, that is linear, near a saturation operating point, in order to obtain a linear model as close as possible to the exact nonlinear model.Algorithm 1 describes how he Cauer network is built.
Algorithm 1: Algorithm for determining the coefficients of the Cauer network Result: The Cauer Ladder Network solution of the selected use case is shown in Fig. 5 when using eight ladder circuits (i.e.n q = 8).Is it noted that further increase of the number of Cauer Ladder circuits n q would no longer improve the solution considerably.We can clearly note a similar trend as for the finite element solution shown in Fig. 3, especially for the A field, but the B field seems to show large differences, especially when it comes to the in-depth penetration of the magnetic field.
The Cauer Ladder Network Method is an excellent method to simulate linear problems and linear dependencies.However, in highly nonlinear problems exhibiting large eddy currents and magnetic saturation, like the one showcased in this work, its prediction power degrades dramatically with respect to the reference solution.

Hybrid twin construction
This section aims at building-up a hybrid twin capable of improving the fast physics-based model presented in Sect.'A fast physics-based modeled solution'.The main driving reason to consider an improvement of the Cauer Ladder Network Method is its ability to capture a large part of the solution within a fraction of a second per time step, using a standard PC.
To start the hybrid modeling of the magnetic levitation, we first compute a reduced basis of the error, using only 10 out of the 14 available simulations in the database.To compute the reduced basis, we use the traditional Proper Orthogonal Decomposition technique.Thus, the solution errors in the physical space (x, y), are placed, at every time-step, as a vectors in a large matrix of error snapshots X.For instance: with A FEM j the finite element method's solution at time step j and A Cauer j the Cauer Ladder Network's one.Later on, a singular value decomposition (SVD) is used to select the most relevant singular vectors of the solution.For instance, we can note: V is the matrix of space singular vectors, the diagonal of the matrix contains the sorted decreasing singular values λ k , while T contains the time vector basis.The selected reduced basis vectors V k of space variation of the solution are selected using all k satisfying: where is a tolerance value.In this use case, the tolerance is set to = 10 −4 .This tolerance has lead to the use of approximately 120 vectors V i in the correction's reduced space.The correction term C is now supposed to have the form: with n = 120 the number of vectors in the reduced basis, and α i are the reduced coordinates, or the weights of the reduced basis vectors V k .
To efficiently compute α i within a very limited time, and possibly in a non-incremental manner, α i are supposed to be the output of a surrogate model H built using neural networks.To take into consideration the memory effect of the solution, depending on previous inputs/outputs of the system, the surrogate model H uses a memory layer.For instance, the time dependency of the solution on previous inputs/outputs combination is accounted for using a Long-short Term Memory layer (LSTM layer) [40].The use of LSTM requires the input to have several time-steps, and not only the current one.The inputs of the surrogate model are in fact the imposed electrical tension U , the computed Cauer Ladder Network currents in the eight branches c j for j = 1, • • • , 8, the Cauer computed force F , the magnetic flux φ, and the time t.The force F and the magnetic flux φ are our global quantities of interest computed from a post-processing of the Cauer Ladder Network's solution A Cauer .All these inputs are given in the form of the last κ time-steps, such as at every time step i we can write: Using this notation, we accept that the first κ time steps will not be corrected, and the correction will only start after κ time steps are computed.In the selected surrogate model H, a single dense fully connected layer follows the LSTM layer and outputs the n required values in the α vector, at every time step.One should note that the inputs of the surrogate model depend only on the Cauer Ladder Network solution and not on the correction itself.Thus, all α i for all time steps can be computed simultaneously.
The surrogate model H is named a spectral network.In fact, as previously explained, it uses the spectral decomposition of the field, through a singular value decomposition, and reconstructs the solution in a form of a reduced order model.The weights are predicted as a function of the current and previous time inputs, while fixing the reduced basis of the approximation V.

Remark 1
The corrected solution is finally computed such as: Fig. 6 Corrected hybrid model solution for the same selected use case shown in Fig. 3, illustrated at t = 5s however, since α k (t) is nonlinear as a function of global quantities computed using A Cauer , the correction is nonlinear as a function of the Cauer Ladder Network's solution.
To train the network, the weights and biases of the two selected layers can be computed using a gradient descent algorithm and the following loss function L: where W are the weights to be optimized in the neural network and β a regularization coefficient.The optimization is performed using ADAM optimization algorithm with a custom adaptive learning rate.The gradients are computed using automatic differentiation in Keras−Tensorflow toolbox, using tensorflow's GradientTape functionality.The training is only performed on the selected 10 data sets, out of the 14 available simulations.The solution of the corrected fields A corrected for κ = 5 is shown in Fig. 6 for the same use case illustrated in Fig. 3.The solution shows a good improvement with respect to the Cauer Ladder Network solution.One should note that the solution shown in Fig. 6 belongs to a training set.Another solution belonging to the testing set is illustrated in Fig. 7, which shows simultaneously the finite element solution in 7a, the Cauer Ladder Network solution in Fig. 7b, and the corrected solution in Fig. 7c.The testing solution shows the amplitude of the magnetic field B , and also showcases a large improvement of the solution after correction.This improvement is more tangible when exploring global quantities of interest, as illustrated in sect.'Application to the magnetic bearings modeling's quantities of interest'.
One should note that the computation of the solution's weights and correction, for all time steps, was performed within 1 to 1.1 s, for any of the 14 available use cases.
The relative error maps on the testing set are shown in Fig. 8.The error for the Cauer Ladder Network reaches more than 50% on the A fields.The error of the corrected Cauer by the proposed Spectral Net does not exceed 11%.Moreover, the error is localized away from the region of interest, which is mainly around the coil and shaft separation region.

Application to the magnetic bearings modeling's quantities of interest
In this section we illustrate the global quantities for all the training and testing sets, as a function of time.The global quantities of interest are the force F applied on the steel shaft and the magnetic flux φ.   (20), and E max is the maximum relative error in percentage, as defined in Eq. ( 21): Z being the quantity of interest and K the number of time steps available in the sequence.The testing set solutions, post processed to compute the global quantities of interest, are illustrated in Figs.11 and 12 for the force and the flux respectively.The testing set shows an excellent agreement of the hybrid twin solution with the finite element one, and a tremendous improvement when compared to the Cauer Ladder Network results.An effort to model the global quantities was taken by the authors using fully data-driven techniques, using many advanced machine learning techniques.When not infusing physical information in the model, the solution remains unstable and couldn't be generalized to different scenarios or use cases.The authors tested an LSTM method, a Resnet with embedded LSTM to compute the derivatives, as well as a NeuralODE, based on fully datadriven approximation of the derivatives.However, none of the tested techniques was able to produce stable results for thousands of time steps, within a fraction of a second.The presented results prove the robustness of the proposed approach.Many of the testing sets were performed at a different maximum tension U , different scenario schemes, and different frequencies, some never seen in the training zone.The algorithm can perform in extrapolation with very reasonable outputs, while computing the results of about 14,000 time steps within about 1 s on a usual PC.The results show excellent accuracy while the simulation performs in almost real-time.When it comes to control, the inference of a single time step is currently about 200 ms using Python libraries, limited by the Cauer ladder Network computation requirements.Nevertheless, the simulation of few thousand time steps takes about 1 to 1.1 s, leveraging the vector and tensor computation possibilities.The only calculation performed online is the prediction of the basis weights α i using the neural network at every time step i.
The required computation time in the electronic control is less a feedback frequency of 2000 Hz, that is around 0.5ms per feedback.The resultant algorithm is not implemented on the electronic material yet.A speedup of 100 times between python computing and PCB computing is realistic and achieved in multiple previous applications.The algorithm achieves the required computation time on PC.The code optimization and the use of dedicated PCB material achieves the required computation time for control purposes.

Discussion and conclusions
This work proposes a hybrid modeling, enhancing the available state of the art simulation of magnetic field propagation, without impacting computation time.The results show a relevant enhancement of the solution throughout all the simulated time domain.Moreover, the solution is stable and reliable, as depicted in the different testing cases.The Cauer + Specktral Net solution enhanced the Cauer Ladder Network solution, while preserving stability and computation time.
In this work, we illustrated a novel machine learning hybrid twin architecture, which relies first on the use of spectral decomposition of a solution, in a form of a reduced order basis.Later on, the method finds the solution weights using a machine learning algorithm.
The proposed spectral network improves a physics-based simulation by computing a correction, using a reduced basis of the error with respect to a ground truth, high fidelity, finite element solution.This work uses the finite element solution as the ground truth values, since experimental measurements are missing.The same algorithm can be used to build hybrid modeling corrections fitting the experimental measurements when available.
Moreover, the algorithm uses explainable physics, and its limitations can be understood, and controlled to a certain extent, using the basic principles of electromagnetic.

Fig. 1 A
Fig. 1 A schematic drawing of the magnetic levitation bearing in use

Fig. 2
Fig.2The adopted mesh considered for the high fidelity finite element simulation of the magnetic levitation bearing

Fig. 3 A
Fig. 3 A finite element solution for a selected use case solved in a time domain t ∈ [0; 5]s, and illustrated at t = 5s in this figure

Fig. 4 A 3 Fig. 5 A
Fig. 4 A equivalent Cauer ladder network with n q = 3

Fig. 7 AFig. 8
Fig.7A solution from the testing set at the same time step.We remind the reader to mind the axis

Figures 9
Figures 9 and 10 illustrate the global quantities for the training sets, for the forces and the magnetic flux respectively.Figs. 9 and 10 also compare the finite element solution, the Cauer Ladder Network's solution, and the corrected hybrid twin solution.in all the figures E m stands for the mean relative error in percentage point-wise, as defined in Eq.(20), and E max is the maximum relative error in percentage, as defined in Eq. (21):

Fig. 9
Fig. 9 Post processing of the force F in the training set

Fig. 10
Fig. 10 Post processing of the magnetic flux φ in the training set

Fig. 11 Fig. 12
Fig. 11 Post processing of the force F in the testing set

Table 1 A
review of the available use cases in this work