Reduced order modelling using neural networks for predictive modelling of 3d-magneto-mechanical problems with application to magnetic resonance imaging scanners

The design of magnets for magnetic resonance imaging (MRI) scanners requires the numerical simulation of a coupled magneto-mechanical system where the effects that different material parameters and in-service loading conditions have on both imaging and MRI performance are key to aid with the design and the manufacturing process. To correctly capture the complex physics, and to obtain accurate solutions, ﬁnite element simulations with dense meshes and high order elements are needed. Reduced order model approaches, based on the established proper orthogonal decomposition (POD) approach, are attractive as they can rapidly predict the numerical simulations needed under changing parameters or conditions. However, the projected (PODP) approach has an invasive computational implementation, whilst the interpolated (PODI) approach presents challenges when the dimension of the space of parameters to be investigated becomes large. As an alternative, we investigate a POD technique based on using a neural network regression, which is not as invasive as PODP, but has superior approximation properties compared to PODI. We apply this to the coupled magneto-mechanical system to understand three pressing industrial problems: ﬁrstly, the accurate and rapid computation of the resonant frequencies associated with this coupled magneto-mechanical system, secondly, the effects of magnet motion on the Ohmic power and kinetic energy curves, and, thirdly, the prediction of the uncertainty in Ohmic power and kinetic energy curves as a function of exciting frequency for uncertain material parameters.


Introduction
The use of magnetic resonance imaging (MRI) [1] has been an essential tool in modern medical diagnosis due to the scanners' high in-built resolution when imaging fractures [2], joints [3], and soft tissues, such as damaged cartilage [2] or tumors [4].
Time varying magnetic fields generated as part of the imaging process by an MRI scanner interact with the conducting shields and generate unwanted eddy currents.The purpose of these shields is to isolate or block the strong static magnetic field generated by the main MRI magnet from the region outside the scanner.However, gradient coils, which generate weaker time varying magnetic fields, generate eddy currents and give rise to electromagnetic stresses causing deformations and vibrations within the MRI scanner [5][6][7].This, in turn, can lead to patient discomfort and ghosting.Still further, eddy currents and vibrations release energy, which can cause helium boil-off, an expensive coolant used to keep the magnets superconducting, and this may, ultimately, lead to a magnet quench.Vibrations are also caused by external machinery (e.g.lifts, cleaning equipment and building maintenance), which can lead to eddy currents being generated in conducting components of the scanner and may cause disturbance to the uniformity of the strong magnet field generated by the main magnet and also lead to ghosting effects.Manufacturers of MRI scanners want to predict and understand these vibrations to formulate a better understanding of the coupling phenomena and design new and better scanners, which exploit alternative materials.To do this, they need to carefully optimise the configuration of the model parameters in order to ensure that the new magnets meet their design specification.
To help design scanners, in collaboration with our industrial partner Siemens Healthineers, this paper addresses three pressing industrial problems: firstly, the accurate and rapid computation of the resonant frequencies associated with this coupled magneto-mechanical system, secondly, the effects of magnet motion on the Ohmic power and kinetic energy curves, and, thirdly, the prediction of the uncertainty in Ohmic power and kinetic energy curves as a function of exciting frequency for uncertain material parameters.Accurate full order model simulations require fine finite element discretisations using dense meshes and/or high order elements, and can take significant time to compute [8][9][10][11][12].To reduce the cost of solutions for new parameters, proper orthogonal decomposition (POD) reduced order model (ROM) schemes [13][14][15][16] have been explored that use either Galerkin projection or interpolation [17,18].However, the projected POD (PODP) approach has an invasive computational implementation, while the interpolated POD approach (PODI) presents challenges when the dimension of the space of parameters to be investigated becomes large.As an alternative, we investigate a POD technique based on using a neural network regression [19,20], which is not invasive but has superior approximation properties compared to PODI.Using multilayer perceptrons (MLP) neural networks in combination with POD has been employed by Hesthaven and Ubbiali [21] and this technique is called PODNN (or POD-NN) in the literature.The use of neural networks in the context of POD range from the prediction of wind pressure and wind induced responses for high-rise buildings [22] with a emphasis on back-propagating networks, to the use of recurrent neural networks (RNN) closure of parametric POD-Galerkin ROMs [23].In the latter, it is shown that a long short-term memory (LSTM) type of RNN can significantly improve accuracy and efficiency for non-linear problems, even beyond the time interval of training data.More recently, the use of other deep learning based reduced order models (DL-ROMs) in POD has been considered in [24,25] with a more extensive literature review given in [26] covering artificial neural networks (ANN), physics informed neural networks (PINNs) and feed forward neural networks and the differences between them.The application of this technology to the three aforementioned problems will demonstrate the advantages of using the neural network based POD over our previous POD approaches.
The structure of the paper is as follows: Sect. 2 briefly outlines the mathematical model of the problem of the coupled magneto-mechanical problem, in the form of a transmission problem, that we wish to tackle in this project.In this section, we also describe the finite element discretisation and the structure of the algebraic system to be solved.Then, in Sect.3, we describe the reduced order modelling approach based on POD and describe the offline and online stages.For the online stage, we focus on PODNN and show comparisons with other POD approaches.We also emphasize the importance of how the neural networks should be applied to POD to ensure accurate results efficiently.Sect. 4 presents a series of numerical examples where the PODNN technique is applied to the problems of interest in this work.This includes the prediction of resonant modes in the Ohmic power and kinetic energy curves using different parameter spaces, predicting the response in the same outputs of interest for different external vibrational loading conditions and predicting the response in outputs of interest when model parameters are uncertain.Finally, in Sect. 5 we conclude our findings and recall the work carried out in this paper.

Magneto-mechanical model
To computationally simulate the coupled electro-mechanical interactions in a MRI scanner we choose to follow [12] and formulate the equations in a Lagrangian setting.We recall that this mathematical model assumes that both the displacements and the velocities of the conducting components are small.The problem of interest is described by Fig. 1, where C denotes a conducting region, c C := R 3 \ C denotes the (unbounded) non conducting region and ∂ C represents the boundary of the conducting region.The unknown fields for this problem are the magnetic vector potential A and the mechanical displacement u, while J is a known external current source.
Continuing to follow [12], by linearising about the static (DC) solution, the transient (AC) problem becomes linear in time (t) dependent fields and the latter can be represented in a time-harmonic fashion so that Fig. 1 General representation of the magneto-mechanical problem, illustrating a conducting region c with magnetic permeability μ = μ * and electrical conductivity γ = γ * that is within an unbounded nonconducting region c C with μ = μ 0 and γ = 0.An excitation is applied by a current source J s that is prescribed in the coils with a similar decomposition for J, where ω = 2π f is the angular frequency of the driving excitation, f is the frequency in Hz and i := √ −1.A linearised transmission problem has already been derived in [12,27] for this problem leading to the strong form where γ , ρ and μ are the electric conductivity, mass density and magnetic permeability, respectively; n is a unit outward normal (pointing from the conducting to the non-conducting side); σ e (U U U AC ) is defined as the mechanical contribution to the Cauchy stress tensor; representing the nonconducting and conducting side of ∂ C , respectively.Once (2) has been solved, the complex amplitudes of the electric and magnetic AC fields in the Eulerian setting can be recovered as where B B B AC 0 = curl A A A AC and B DC 0 = curl A DC and the complete physical fields for both stages are In this paper, we assume that the DC fields are known.Our goal is to provide a new reduced order model procedure to allow the rapid prediction of approximations to A A A AC and U U U AC for new material parameters and loading conditions.In turn, this allows us to obtain associated outputs of interest as well as to understand the uncertainty in outputs of interest for uncertain material parameters.To this end, we truncate c C a finite distance from C , and create the truncated domain .On ∂ , we impose n × A A A AC = 0 as an approximation to the decay condition (2d).To circumvent the Coulomb type gauge in (2b), we add a small perturbation εA A A AC to the left hand side of (2a) in \ C and replace A A A AC with the regularised solution A A A AC ε .We then use the finite element method (FEM) to approximate full order model solutions to A A A AC ε and U U U AC over the .For further details we refer to [12].In the next section, we briefly recall the finite element (FE) discretisation employed.

FEM discretisation
FEM discretisation involves partitioning the domain into non-overlapping regions called elements and then approximating the solution over each of these elements by appropriate piecewise polynomials.We will use an unstructured grid of tetrahedral elements due to the automatic procedures that are available for generating these meshes around complex configurations and adopt a higher order (hp-version) FE framework allowing the combination of refining the mesh spacing (h) and increasing the order ( p) of these polynomials in order to improve the accuracy of the solution [8,10,11], as was applied in the case in our previous work [12].The use of such technology allows accurate solutions to be obtained and, in the case of p-refinement for smooth solutions, and combinations of h-and p-refinements, for nonsmooth solutions due to sharp edges and corners, is known to lead to exponential convergence of the solution [8].The mathematics and continuity requirements of the fields dictate that the mechanical displacement and magnetic vector potential should be approximated differently, for details of the weak form of (2) and its derivation we refer to [12].This then means that H(curl) conforming basis functions be used for A A A AC ε while H 1 conforming basis functions should be used for each component of U U U AC .The particular choice of basis functions employed correspond to those proposed by Schöberl and Zaglmayr [28,29].This means that fields A A A AC ε and U U U AC are approximated by where the subscripts hq and hp are used to indicate that these are discrete approximations on a mesh of spacing h with different orders q and p, additionally, P AC global and Q global refers to the number of degrees of freedom (DOF) in the AC system in the same respective manner.Still further, N g and L s are typical H(curl) and H 1 basis functions respectively [17].

Full order model
Using the FEM recalled above, the application of the linearised strategy in [12] leads to the need to solve parametrised linear systems of the form: Find A AC and U AC such that where K AC AA , C AC AA , K AC UA , K AC UU , K AC UU and M AC UU refer to the system matrices that are obtained by discretisation of the weak forms.The explicit definitions for these matrix blocks can be seen in [18] with K, M and C referring to stiffness, mass and damping type matrices, respectively.The R AC  A term drives the solution to the problem.The structure of (6) means that a computational efficient solution strategy is first employed to solve A(w) q(w) = r(w) , (7) for q(w) := A AC ε (w), where A(w) := K AC AA (w) + C AC AA (w) and r(w) := R AC A (w), we have emphasised how the solution may depend on parameters of interest w, and have introduced the simplified notation to help with later sections.Once (7) has been solved, U AC can be obtained from solving and, in the case of a non-zero Dirichlet boundary condition in (2c), to the right hand side in (8), an additional contribution corresponding to columns of K AC UU − ω 2 M AC UU that have been moved to the right hand side multiplied by the corresponding solution coefficients found from the Dirichlet condition [9], should be subtracted.Of the systems ( 7) and (8), in practice, the former is larger and more computationally expensive to solve since, typically, P AC global Q AC global as C is just a small part of and the latter only needs to be solved for those mechanical degrees of freedom in C .Henceforth, following the justification considered in [17], the former will be our focus in the following section.

Reduced order modelling
To aid with the development and design of a new MRI scanner, linear systems of the form (6) must be repeatability solved a large number of times for different parameters (such as frequency or conductivity of radiation shields), which can lead to exponential growth in computational cost.To address this, we consider a reduced order model (ROM) to optimise work flow and reduce the overall computational cost.In particular, we consider a form of proper orthogonal decomposition (POD) approach e.g.[14,30,31], which seeks to extract modal information from a small number of solution snapshots corresponding to sets of snapshot parameters and use the most important of these to construct a problem of reduced size.These ROMs have been successfully applied to a range of engineering problems and, in the context of the linearised magneto-mechanical problem [17], interpolated (PODI) and projected versions (PODP) of the POD for the system (7) have been considered.In this work, we investigate a POD technique based on using a neural network regression (PODNN).Our implementation is equivalent to that proposed by Hesthaven and Ubbiali [21], and is less invasive than PODP and is amenable to a large number of potential parameters w.We briefly outline the offline (shared by all POD approaches) and the differing online stages below.We also provide some practical insights in to how PODNN should be applied to achieve accurate results.

Offline stage
Computed solutions to (7) for different sets of (snapshot) parameters w = w j ∈ R P are arranged in D as: D = (q(w 1 ), q(w 2 ), q(w 3 ), . . ., q(w N s )) ∈ C η×N s , (9) with N s representing the number of computed snapshots and η = P AC global .To extract the modal information, a singular value decomposition (SVD) is applied so that D = H G H , where H ∈ C η×η and G ∈ C N s ×N s being unitary matrices, ∈ R η×N s being a diagonal matrix padded out with zeros and ii being the singular values, arranged in descending order decay towards zero and H denoting the matrix Hermitian.Then, by dropping those singular values corresponding to ii / 11 ≤ TOL, D can be well approximated by a truncated singular value decomposition (TSVD) where m is the number of singular values that are retained.Since m < N s << η, the truncated matrices of left and right singular vectors follow the form We would like m to be small whilst also maintaining a good approximation to D.

Online stage
The online stage involves approximating the solution of (7) for new sets of parameters at reduced computational cost compared to the offline stage.The aim is to calculate an approximation q(w) for new w that are different to the snapshot parameters w j .

Interpolated POD (PODI)
In the case of interpolated POD (PODI), the approximation is typically chosen as e.g [17,32] where I k (w, (G m ) H ) is an interpolation through the kth row of (G m ) H or equivalently the kth column of G m , with the overbar denoting the complex conjugate, for the parameter set of interest w.For further details of the implementation see [17].

Projected POD (PODP)
In the case of projected POD (PODP), the form of the approximation is chosen as e.g [15,17,33] q PODP (w) = H m p(w), (12) where p(w) ∈ C m is obtained from solving a small parameter dependent linear system obtained by Galerkin projection of ( 7) where For further details of the implementation see [17].

Neural network POD (PODNN)
In the case of a neural network POD (PODNN), one possible choice is analogous to (11) i.e.
where a Neural Network is constructed to provide a regression (instead of interpolation) through the kth row of (G m ) H or equivalently the kth column of G m ).An alternative is to use The version in (15) can be seen to be equivalent to the strategy proposed by Hesthaven and Ubbiali [21], which, in the above notation, is (10).In both the case of ( 14) and ( 15), the neural network regression involves creating a function defined by parameters p((G m ) H ) or p( m (G m ) H ), respectively, which enables the prediction of q PODNN at low cost for new parameters w.To establish the advantages that (15) offers over ( 14), we consider the ability of neural networks to capture oscillatory functions in the next section.

Capture of oscillatory functions by neural networks
As a simple example, we wish to consider the ability of a feed forward neural network to capture a vector y(x) ∈ R 100 of oscillatory functions with components where w k are amplitudes of the functions and x ∈ [0, 1].To do this, we set up a dataset with N samples in the form of ordered pairs D := (x (1) , y (1) ), (x (2) , y (2) ), . . ., (x (N ) , y (N ) ) , (17) where x ( j) = x ( j) is a scalar input feature in this case and y ( j)  is the corresponding vectorial output labels.Using a mean squared error (MSE) as a loss function, the training of a neural network for this problem corresponds to finding parameters p such that where R(x ( j) , p) describes the regression function provided by a particular Neural Network architecture.We start with a neural network with L = 1 layers and n = 10 neurons in each layer and a dataset with N = 101 entries.In Fig. 2, the results obtained for the cases where a single network is trained to approximate y(x) with w k = 1, separate networks are trained to approximate y k (x) with w k = 1 individually and a single network is applied to capture y(x) with w k = 1/k 2 (where the oscillatory functions have reducing amplitude) are compared.The MSE loss for the single network with In this figure, we observe that a single network is unable to capture the oscillatory functions for w k = 1, and there are clearly errors present even in the prediction of y 1 (x).On the other hand, training separate networks capture the behaviour of y(x) well with w k = 1.However, if the oscillatory functions are such that w k = 1/k 2 , a single network is able to capture the behaviour of the lowest modes accurately.
As L = 1 and n = 10 may not represent the best choice of hyper parameters, we now apply a Bayes optimisation parameter search [34] to determine the best choices of L opt and n opt in each case so that L( p) is minimised.Other alternative optimisation strategies could be used, but we expect them to lead to similar results.The resulting loss surfaces are shown in Fig. 3 and, in Fig. 4, we repeat the investigation shown in Fig. 2 for the new hyper-parameters obtained.Again, we can state the MSE loss for the optimised single network with  4 shows that a single network with optimised hyper-parameters is unable to capture well y(x) when w k = 1, but is able to capture the behaviour well when w k = 1/k 2 .When networks are trained to predict y k (x) individually, the prediction is accurate for the first few oscillatory functions but as they become more oscillatory, the loss is higher, as seen on Fig. 3 for the last oscillatory function, which hence means they are unable to capture y k (x) for higher k.
The computational costs of using a single neural network with w k = 1, multiple networks with w k = 1 and a single network with w k = 1/k 2 are compared in Fig. 5 for both a fixed choice of L and n and also those obtained from the Bayes optimisation L opt and n opt .These, together with the aforementioned plots, illustrate that a single network can perform well and is efficient provided that the oscillatory functions have reducing amplitudes.

Implications for PODNN
Returning again to ( 14) and ( 15), we see that the latter involves R k w, p m (G m ) H , where the rows of m (G m ) H represent functions that are both oscillatory and reducing in amplitude, whilst the former involves R k w, p (G m ) H , with the rows of (G m ) H being oscillatory but not reducing in amplitude.Thus, by Sect.3.2.4,the latter case is preferred as it requires less computational effort and requires only the training of a single network.In practice, the datasets we will use for PODNN can be setup in a similar way to D described above and, for later use, it is useful to introduce x ( j) = w j for linear spaced snapshots, log 10 (w j ) for logarithmically spaced snapshots, (19) where the log 10 (w j ) should be interpreted as log 10 (w 1 j ), log 10 (w 2 j ), . . ., log 10 (w P j ) for P parameters and use a similar approach for describing the continuous variable x(w).Analogous to (18), we also state the corresponding MSE loss function for this case as where the notation :, j implies all rows and the jth column.

Numerical examples
Numerical examples that demonstrate the advantages of the PODP scheme over a PODI scheme in the magnetomechanical problem described in Sect. 2 has been presented in [17] and so we focus on PODNN and the additional benefits this offers.We begin by describing the test-magnet geometry in Sect.4.1, which will be the focus of our comparisons presented in this work.Then, in Sect.4.2, we describe the software used to generate our numerical results followed by a description of the discretisation in Sect.4.3.We present predictions of resonant modes in the Ohmic power and kinetic energy curves using different parameter spaces for the test magnet geometry in Sect.4.4 followed by the predictions of the response in the same outputs of interest for different external vibrational loading conditions in Sect.4.5.Then, in Sect.4.6, the PODNN technique is then also used to predict the response in outputs of interest when model parameters are uncertain.

Test magnet geometry
Our focus will be a simplified MRI configuration called the test magnet problem, which is rotationally symmetric, and was previously considered in [12].For the simplified model, only the rotational symmetric z-gradient coil is considered and the non-rotationally symmetric x and y gradient coils are removed.If desired, this allows the problem to be reduced to a axisymmetric problem although we will model it as one   quarter of the full three-dimensional geometry. 1An illustration of this MRI configuration can be seen on Fig. 6, which shows a typical 2D cross-section and 3D view of the general setting.In this figure, the conducting shields (4 K and 77 K) and the outer vacuum chamber (OVC) are shown in different shades of blue and these will deform under the presence of electromagnetic stresses generated due to the presence of eddy current currents in the conducting components.The main coils and gradient coils are shown in shades of red.The dimensions, exciting currents and materials of this problem are commercially sensitive, but indicative values are provided in [35].Unless otherwise stated, we assume that ) to reflect the fact that the magnet geometry is fixed in position.Furthermore, the transmission problem described by ( 2) must be supplemented by the symmetry boundary conditions [17] that are to be imposed on the x = 0 and y = 0 planes.As outputs of interest for the MRI configurations involve the output (dissipated) power P 0 (ω, A A A AC ε , B DC 0 , U U U AC ) and the kinetic energy E k (ω, U U U AC ) which are defined as: as a function of ω.In practice, our approximation to these will be evaluated for a particular conducting shield C .At an intermediate stage, we will also consider plots of approximations to 1 Of course, if desired, a different fraction of the full geometry could be considered provided the correct symmetry boundary conditions are imposed. 2Reminder that a represents the complex conjugate of a vector field a.

Software
The magneto-mechanical model will be simulated using a simulation tool that has been extended from the previously developed higher order finite element and reduced order model tool developed by our group to simulate the problems described in [12,17,18].The latest additions include the extension to handle the PODI and PODNN reduced order models for large parameter space w, uncertainty quantification features and non-zero Dirichlet displacements of the shields.The neural networks were implemented with MATLAB's Deep Learning Toolbox [36], which has similar functionalities to popular python packages such as tensorflow [37], scikit-learn [38] or pytorch [39], since our simulation tool was also written in MATLAB [40].In particular, the feedforwardnet tool was used to build multilayer perceptron (MLP) neural networks.In the following, we focus on numerical results for the new reduced order models as the accuracy of the linearised strategy in (2) and the hp-finite element discretisation has already been established in [12].
For PODNN, the split between D train and D test is 85% and 15%, respectively, and is applied throughout.Prior to training using D train , values were normalised to the standard deviation and mean of ( m (G m ) H ) 1,: .Additionally, the form of the activation function followed a tan-sigmoid function.To solve the minimisation problems (20) the Bayesian regularisation solver with the option to use the Jacobian for calculations was employed with the maximum number of epochs set to 5000 and min grad= 10 −10 .Bayes optimisation, implemented with MATLAB's Statistics and Machine Learning Toolbox [41], is applied to determine the optimum number of layers L opt and neurons per layer n opt , in a similar way to Sect.3.2.4 with other hyper-parameters set as per [36].All computations were performed on a workstation that comprised of a 12-core Intel Xeon W-2265 Processor with 128GB (8 × 16GB) RAM and a NVIDIA Quadro RTX 4000 8GB graphical processing unit.

Discretisation
For this problem, c C is truncated to form a finite computational domain by excluding the region outside of a suitably sized cylinder.We discretise by an unstructured mesh of 33 805 tetrahedral elements and employ a discretisation with q = 3 and p = 3 order elements for A A A AC ε,hq and U U U AC hp , respectively.To dampen the resonant peaks in the model, we employ a fixed damping ratio of ξ = 2×10 −3 and to circumvent the Coulomb gauge we use a regularisation parameter of ε = 10 −4 .This choice of discretisation has been shown to provide accurate full order model solutions to the test magnet geometry over the parameters of interest [12].We first consider the effects of different choices of snapshots and the behaviour it has on the decay of ii .Considering the case of P = 1 and w = (ω) = 2π( f ) with 5 ≤ f ≤ 5000 Hz, we will explore five different selections of snapshot frequencies w i = 2π( f i ) and decide which will be most appropriate for our problem.From previous analysis in [17], we know that the test-magnet geometry is sensitive to lower frequencies and, hence, may require more frequencies within that lower range.Secondly, we consider the choice of snapshot parameters for case of P = 2 with w = (ω, γ OVC 0 ) where 0.1 ≤ γ OVC 0 ≤ 10 is a non-dimensional factor used to scale the conductivity in the OVC such that where γ OVC ref is a reference value for conductivity of the shield.
In [17] a piece-wise linear choice for the snapshot frequencies was made, but it is not clear that this is the best choice.Therefore, we compare the performance of selecting the snapshots f i in a linear, piece-wise linear and logarithmic fashion in the following.Linear The linear snapshot frequencies f i in Hz are chosen such that: for i = 1, . . ., N fs , where N fs is the total number of snapshot frequencies and, in the case of P = 1, N s = N fs .

Piece-wise linear
The piece-wise linear snapshot frequencies f i are chosen by first introducing the splitting N fs = N (1) fs , which are then used to calculate the spacings where the numerators are determined the range of frequencies that we intend to split by.The frequencies f i are then constructed as fs ) f (2) for N (1) fs , (27) with N (1) fs and N (2) fs chosen as in [17].Logarithmically spaced For logarithmically-spaced frequency snapshots, f i is chosen such that log 10 ( f i ) = log 10 (5) + log 10 (5000) − log 10 (5) with i = 1, . . ., N fs .
For N fs = 180, we compare the distribution of f i according to (25), ( 27) and (28) in Fig. 7, with the first being equally spaced and second and third showing a clustering of samples for smaller frequencies as expected.Making these different choices of snapshots to generate full order model solutions q(w j ), creating D in (9), and then applying the TSVD in (10) leads to the decay of ii shown in Fig. 8. Comparing the different choices, the log-spaced snapshots in (28) leads to fastest decay of the singular values.Setting T O L = 10 −6.5 (due to the decay in Fig. 8 indicating that the decay is no  27) and m = 20 for (28), hence, leading to the smallest reduced order model when the log-spaced frequencies are used.An additional benefit of ( 28) over ( 27) is that we do not need to define a priori where to split the piece-wise linear spacing of frequencies.
The complex amplitudes |( m (G m ) H ) i j | shown as lines for a fixed i = 1, . . ., 10 against f j , j = 1, . . ., N fs , are included on Fig. 9 for the different choices of snapshots.The general trend for linear, piece-wise linear and logarithmically spaced snapshots is that the lower modes (e.g.i = 1, 2) are simple smooth functions that are not very oscillatory, but, as i increases, the modes become increas-ingly more oscillatory.For higher modes, the oscillatory nature is very pronounced for small f j and Fig. 9d indicates there are additional benefits by showing |( m (G m ) H ) i j | against f j , j = 1, . . ., N fs on a log-log graph indicating that x ( j) = log 10 ω j = log 10 (2π f j ) may be a good choice in this case.Henceforth, we choose the frequency snapshots according to (28) and use x ( j) = log 10 ω j , but in the following we will not only consider N fs = 180, but also cases of N fs = 13, 23, 45, 90 to reduce computational effort needed for the off-line stage.The snapshots N fs = 45, 90 follow a similar decay and lead to models of size m = 20, but, the N fs = 13, 23 cases are seen to decay faster and henceforth a reduced model of size m = 18.However, as seen in Sect.4.4.2 for N fs = 13, it will not be a sufficient amount of snapshots to accurately capture the solution of the problem.
Having found that a logarithmic spacing f i is preferable, we now consider how best to choose the conductivity factor snapshots γ OVC 0,i .Given that γ and ω appear as a product in (2a) we argue that the best choice for the spacing of the snapshots γ OVC 0,i is also logarithmic.Hence, we propose log 10 (γ OVC 0,i ) = log 10 (0.1) + log 10 (10) − log 10 (0.1) with i = 1, . . ., N cs and N cs is the chosen number of snapshots for the conductivity factor.This means that for the case of P = 2 we have N s = N cs N fs snapshots in total.The decay of singular values, shown in Fig. 8 for P = 1, can also be found for P = 2 and shows a similar behaviour.In this case, however, setting T O L = 10 −6.5 corresponds to a reduced model of size m = 44 if N s = 23 × 23 and, hence, still offers a considerable saving.

Online stage
Predictions for w = (ω) = 2π( f ) with P = 1 In the online stage of either PODI or PODNN we can easily evaluate ( 23) and ( 22) for any frequency of interest.However, for the purpose of visualisation, we will show results for the case where the ROMs are evaluated for the following output frequencies in Hz for i = 1, . . ., N fo , with N fo = 500, unless otherwise stated.We use Bayes optimisation as in Sect.3.2.4 to find L opt and n opt using the search space L = 1, . . ., 5 and n = 1, . . ., 16.The resulting MSE loss L PODNN ( p) is shown on Fig. 10a for N fs = 180.As previously observed in Fig. 3, a similar L PODNN ( p) trend is observed with L PODNN ( p) < 10 −6 for n > 5 independent of L in this case.This investiga- (d) as (c) but with logarithmic axis Fig. 9 Test magnet problem with z (longitudinal) gradient coil; offline stage for P = 1, N fs = 180.Complex amplitudes |( m (G m ) H ) i j | as lines for a fixed i = 1, . . ., 10 against f j , j = 1, . . ., N fs , for a linearly spaced f j , b piece-wise linearly spaced f j , c log-spaced f j and d against log f j for log-spaced f j tion will be repeated for N fs = 13, 23, 45, 90 as each will give rise to a slightly different |( m (G m ) H ) i j |.We additionally show in Fig. 10b the ability of the optimised network to fit |( m (G m ) H ) i j | for a fixed i = 1, 9, 19 against f j , j = 1, . . ., N fs for N fs = 180.As expected, the higher modes are more oscillatory and, hence, the fitting becomes comparatively worse for these modes.For example, i = 9 is shown to have reasonable fitting at higher frequencies.However, for i = 19, the network is partially able to capture |( m (G m ) H ) i j | against f j , j = 1, . . ., N fs .Due to the fact that m 19,19 / m 1,1 ≈ 5.9 × 10 −7 , which is close enough to the TOL chosen, the partial fitting of this mode is well enough for a reasonable prediction as seen in Fig. 11.
To illustrate the interdependence between N fs and some of the hyper-parameters of the MLP, we show in Fig. 11

the behaviour of
as a function of ω = 2π f obtained by applying (15) and then (5a) for different N fs and different values of the hyper-parameters L = L opt and n = n opt found by the Bayes optimisation with their respective L PODNN ( p).The different lines in the figure correspond to the cases where C is considered to be the 4K, 77K shields and the OVC in turn, also shown are the corresponding solution snapshot solutions for A A A AC ε,hq (ω) L 2 ( C ) .For N fs = 23, 45, 90, 180 we get good agreement for all shields, but not in the case of the 4K shield with N fs = 13.
Having now established that the performance of PODNN is similar for N fs > 23 snapshots, we show, in Fig. 12 a com- parison between the accuracy and computational efficiencies of PODI, PODP and PODNN.As an approximate measure of accuracy, we introduce where are the evaluations of (23) with the reduced techniques and the full order model, respectively.A more precise calculation would involve the which does not average out possible spatial errors.However, we have high confidence in the PODNN solutions given the results in Fig. 11 and those presented later in this section.
To compare computational efficiencies, we consider the wall clock time of the online stages of the calculation of q POD for the POD techniques.We see that PODP is most accurate with independent of the N fs considered.PODI performs worse than PODP, especially for lower N fs by a few orders of magnitude but improves for larger N fs .PODNN does not perform as well as the PODP case but for small N fs , it is shown to perform better than PODI.However, the decrease in accuracy may be due to the tolerances chosen within the network but the accuracy achieved by PODNN is sufficient for this practical application.Additionally, we show the online computational expense for the , with PODNN showing the best in terms of performance.In Figure 10 and Figure 15 of [17], the total online and offline costs PODI and a full order as well as PODP with a full order were compared for the same problem (but using a different machine).In both cases, PODP and PODI offer similar computational costs, but just like in Fig. 12b, the accuracy of PODP was substantially better than PODI.Discounting training and optimisation time, the total costs for PODNN would be similar to PODP and PODI shown in [17] if the same machine was used.However, the cost of optimisation and training can be significant as Fig. 12c shows.Whilst the training and optimisation forms part of the offline stage in PODNN, this is significantly more expensive than the offline stages of PODP and PODI if included in the costs.
Using the most efficient and accurate PODNN strategy with N fs = 23, we show, in Fig. 13, the dissipated power P 0 C (ω) and kinetic energy E k C (ω) obtained by additionally solving (8) for each output frequency f i of interest using the reduced order model ( 15), field representations (5) and applying (22).The results are in good-agreement with the full order model and, importantly, accurately predict the resonant modes.
To emphasise that this is a 3D magneto-mechanical problem, we visualise the displacement and magnetic flux density fields obtained from solving the PODNN using the aforementioned network within a conducting domain.Figure 14 showcases both the variation in the magnitudes of the real displacement field and the real magnetic flux density field in the outer 4 K shield at output frequencies of f = 1000 Hz and f = 5000 Hz.Bands of higher and lower magnitudes can be seen in the displacement field at 5000 Hz, but at 1000 Hz, there is a higher displacement near the outer edges of the domain.A similar effect can be seen with the real magnetic flux density for 1000 Hz but at 5000 Hz, the domain shows almost uniform field of low displacement in the conducting domain.
for i = 1, . . ., N co , with N co = 4, unless otherwise stated.This means that we evaluate the ROM for a total of N fo N co sets of output parameters.
A similar investigation to that shown in Fig. 11 was repeated for the two-parameter case, again using the Bayes optimisation but with a search space L = 1, . . ., 10 and n = 1, . . ., 32.It should be noted that the optimisation time for this case takes considerably longer than for P = 1 and is more costly than the PODP method, with a time of 5.29 × 10 5 s, which is several days.It was found that, in this case, using a L opt = 3, n opt = 32 network for N fs = 23 and N cs = 23 produces accurate results provided that min grad=10 −10 and the number of singular values that is retained is m = 44 using the same TOL = 10 −6.5 as carried out for P = 1.As an example, we show the behaviour of for this separation is due to the log spaced nature of the conductivity snapshots, with a higher abundance of snapshots occurring less than unity than above, the idea here is to see how well the PODNN model with both those cases.In both conductivity factor regions, the PODNN model is seen to match reasonably well with the snapshot solutions at these specified conductivities.
Henceforth, we fix a PODNN using a L = 3, n = 32 network and snapshot data corresponding to N fs = 23 and N cs = 23 using this choice we show, in Figs.16 and 17, the dissipated power P 0 C (ω) and kinetic energy E k C (ω) obtained by additionally solving (8) for each output frequency f i and conductivity factor γ OVC 0,i of interest using the reduced order model (15), field representations (5) and applying (22), respectively.Alongside, a comparison with the snapshot solutions is evaluated and can be seen to match well with the PODNN solutions.We observe that the higher conducting factors generally lead lower values of P 0 C (ω) and E k C (ω) for each shield and frequency of interest, although the position of the resonant peaks is not substantially affected.

Prediction of the response to external vibrational loading
The consideration of external vibrational loading for the mathematical model described in Sect. 2 introduces a nonzero boundary condition , which reflects the amplitude of the external vibrational loading.Given the nature of the rotational symmetry of the problem considered in Sect.4.1, we consider only the application of non-zero amplitudes in the z direction since specifying other non-zero components would break the rotational symmetry of the problem and necessitate a full three-dimensional computation rather than simulating quarter geometry and using the symmetry boundary conditions in (21) on the x = 0 and y = 0 planes.

Offline stage
The offline stage follows the approach described in Sect.4.4.1 where a PODNN network with L = 3, n = 32 and N fs = 23 snapshots for w = (ω) = 2π( f ) with P = 1 is employed as described previously.

Online stage
Once (15) has been evaluated for a set of parameters of interest, the smaller dimensional mechanical problem (8) is solved.The magnitudes of the external vibrations are usually very small and so present results obtained by imposing U AC D z = 0, 10 −6 , 10 −5 , 10 −4 m, in turn, simultaneously on all the conducting shields.To obtain, we solve (8) for these frequencies and each U AC D z of interest.Following this, we obtain the dissipated power and kinetic energy using ( 5) and (22), respectively.
In Fig. 18, we show the resulting the dissipated power P 0 C (ω) and kinetic energy E k C (ω) curves, where the different lines correspond to different amplitudes of displacement.In these figures, we observe that the dissipated power in the 77K and OVC shields are not substantially affected by introducing non-zero U AC D z , however, there is a noticeable increase in the magnitude of the kinetic energy curves for these shields when using U AC D z = 10 −4 m, although, the frequencies of the resonant peaks remain unchanged.In the case of the 4K shield, we see the most noticeable changes, with the resonant peaks all but disappearing for the case of U AC D z = 10 −4 m and the kinetic energy curve.

Response when model parameters are uncertain
To understand the influence of an uncertain material parameter on the Ohmic power and kinetic energy curves, we may use the previously developed ROM.We focus on the situation where the conductivity γ OVC = γ OVC

Offline stage
The offline stage follows the approach described in Sect. ) with P = 2, as described previously.

Online stage
In this case, N co = 50 samples of γ OVC 0 ∼ N (m, s) are drawn and for each sample, N fo = 500 output frequencies considered.For each combination of f i and conductivity factor γ OVC 0,i , the dissipated power P 0 C (ω) and kinetic energy E k C (ω) are obtained by solving (8) using the reduced order model ( 15), field representations (5) and applying (22), respectively.Using this data, the mean values and 95% confidence intervals of P 0 C (ω) and E k C (ω), at each output frequency of interest, are evaluated and the results shown in Fig. 19.We see, in general, the position of the peaks are not significantly changed, but there is an noticeable difference in the amplitudes, particularly between the frequencies of 1000 ≤ f ≤ 3000 Hz for the 4 K and 77 K shields.The OVC shield does not seem to show as much variation, with the mean and 95% confidence intervals for P 0 C and E k C being almost indistinguishable on this scale.

Conclusion
This paper has presented a practical use of a POD technique based on using a neural network (PODNN), to aid in the design phase of new MRI scanner configurations.We compare the performance of PODNN to different POD techniques including projection (PODP) and interpolating (PODI) for this challenging example.The PODNN methodology employed is equivalent to that proposed by Hesthaven and Ubbiali [21].The PODP and PODI techniques build on our groups earlier work for the coupled magneto-mechanical problem of interest [12,17].The PODP approach produced relatively accurate results with a small numbers of snapshots (N fs and N cs ) but was heavily intrusive on the software and required recalculation of FEM matrices.The PODI method allowed the use of simple interpolants (such as a Lagrangian) and was non-intrusive, but had limitations when considering small snapshots.Instead, the PODNN method, does not have an invasive implementation and is well suited to this and produces accurate results but the optimisation time (using Bayes optimisation to obtain hyper-parameters) in the offline stage increases significantly when the parameter space is P = 2, making the PODP technique more computationally efficient if overall time is considered.
Results have been presented for the prediction of resonant modes in the Ohmic power and kinetic energy curves using different parameter spaces and the prediction of the response in the same outputs of interest for different external vibrational loading conditions.The PODNN technique has also been applied to predict the response in outputs of interest when model parameters are uncertain.These results have shown that an appropriate network architecture, can provide high-fidelity solutions to the magneto-mechanical problem for both the P = 1 (frequency) and P = 2 (frequency and conductivity) parameter cases and provide rapid online prediction of the outputs of interest for the design of MRI scanners.In further work, we also plan to make experimental comparisons with the fields induced in moving magnets and make improvements to the mathematical model presented in Sect. 2 that no longer assumes that the velocities of the conducting components are small.
then using multiple networks with w k = 1 the loss for the first oscillatory function is L( p) = 4.76 × 10 −13 and L( p) = 4.99 × 10 −1 for the last oscillatory function.Finally, the MSE loss for a single network with w k = 1/k 2 is L( p) = 5.13 × 10 −8 .
Single network with wk = 1/k 2

Fig. 2
Fig. 2 Oscillatory functions; feed forward network fit y k (x), k = 1, . . ., 10 with L = 1, n = 10 for a single network with w k = 1, b multiple networks with w k = 1, c single network with w k = 1/k 2 (a) Single network with w k = 1 (b) Multiple networks with w k = 1; first oscillatory function (c) Multiple networks with w k = 1; last oscillatory function (d) Single network with w k = 1/k 2

Fig. 4 Fig. 5
Fig. 4 Oscillatory functions; feed forward network fit y k (x), k = 1, . . ., 10 with optimum hyper-parameters for a single network with w k = 1, b multiple networks with w k = 1 and c single network with w k = 1/k 2

Fig. 6
Fig. 6 Test magnet problem with z (longitudinal) gradient coil; illustration of the components of the problem a in the axisymmetric meridian plane and b 3D view

Fig. 7 Fig. 8
Fig. 7 Test magnet problem with z (longitudinal) gradient coil; Offline stage for N fs = 180.A comparison of linear, piece-wise linear and log-spaced snapshots for 5 ≤ f i ≤ 5000 Hz

Fig. 11
Fig. 11 Test magnet problem with z (longitudinal) gradient coil; online stage for P = 1, comparison of A A A AC ε,hq (ω) L 2 ( ) .Results for a N fs = 13, b N fs = 23, c N fs = 45, d N fs = 90 and e N fs = 180

Fig. 12
Fig. 12 Test magnet problem with z (longitudinal) gradient coil; online stage for P = 1, comparisons of a ( A A A AC ε,hq L 2 ( C) , b online stage computational expense and c PODNN training and optimisation expense for N fs = [23, 45, 90, 180]

Fig. 13 Fig. 14 5 Fig. 15
Fig. 13 Test magnet problem with z (longitudinal) gradient coil; Online stage for P = 1, comparisons of Full order and PODNN using N fs = 23.Results for a dissipated power P 0 C and b kinetic energy E k C shield is uncertain.We consider a non-dimensional γ OVC 0 ∼ N (m, s) with population mean m = 1, population standard deviation s = log e (2) and γ OVC ref is a reference value with units of S/m and of order 10 6 .

5 Fig. 16 = 3 5 Fig. 17 = 3 Fig. 18 Fig. 19
Fig.16 Test magnet problem with z (longitudinal) gradient coil; online stage for P = 2 and N fs = 23, N cs = 23 comparing the dissipated power P 0 C for each of the three conducting shields.Results for a 4 K