Combination of Measurement Data and Domain Knowledge for Simulation of Halbach Arrays With Bayesian Inference

Accelerator magnets made from blocks of permanent magnets (PMs) in a zero-clearance configuration are known as Halbach arrays. The objective of this article is the fusion of knowledge from the magnetic field and material measurements and domain knowledge (magnetostatics) to obtain an updated magnet model of a Halbach array. From Helmholtz coil measurements of the magnetized blocks, a prior distribution of the magnetization is estimated. Measurements of the magnetic flux density are used to derive a posterior distribution by means of Bayesian inference. The method is validated on simulated data and applied to measurements of a dipole of the FASER detector. The updated magnet model of the FASER dipole describes the magnetic flux density one order of magnitude better than the prior magnet model.


I. INTRODUCTION
A SPECIAL type of accelerator magnet, referred to as a Halbach array [1], is made of circularly arranged blocks of permanent magnets (PMs).The magnet system of the FASER detector [2] consists of three 0.57 T Halbach dipole magnets.In [3], it is shown that imperfections of the PM block magnetizations affect the field quality by generating higher-order multipole field errors.
Quality assurance of the FASER dipole consisted of two measurement campaigns [3]: 1) the measurement of the magnetization of each PM block before the magnet assembly using a Helmholtz coil system and 2) the measurement of the magnetic flux density of the built magnet with a 3-D Hall probe mapper.In both cases, uncertainty arises from sensor calibration, stage misalignment, magnetization errors in the PM blocks, and manufacturing errors in the magnet assembly.
In this article, Bayesian inference is used to combine the prior magnetization data obtained from the Helmholtz coil measurements with the magnetic flux density measurements in the magnet, taken with the Hall probe mapper, and the domain knowledge imposed by the magnetostatic problem.This method is validated on simulated data and applied to measured data of a FASER dipole.Inserting the posterior magnetization in the simulation leads to an updated model that predicts the magnetic flux density in the homogeneous field  Insertion of the magnetic vector potential A defined by curl A = B and the constitutive equation in the magnetostatic problem yields the well-known curl-curl formulation of the problem, which is discretized and solved with the finiteelement (FE) software getDP [9].
The FASER dipoles are composed of a series of 12-18 rings with 16 PM blocks per ring (Fig. 1, bottom).We use the index i to indicate the PM block number in a ring and index j to indicate the ring number.The nominal magnetization of the PM blocks For all j = 1, . . ., 12, with nominal orientation 1, top left) and nominal magnetic moments m i = 330 A m 2 .Due to manufacturing tolerances and production errors, the magnetizations M i j of the procured PM blocks differ from the nominal (specified) values.This can be expressed as adding the deviation M i, j to the nominal magnetization The magnetization of each PM block was measured by means of a Helmholtz coil before the assembly of the Halbach array [3].We interpret M i, j and M i, j as random vectors that can be deduced from the measurements.We define the parameter vector p of the PM block magnetization in the 3-D case by The 2-D case is analogous for only the x and y components.The probability space reflects the random occurrence of production errors of the PM blocks.Having measurements for multiple PM blocks of the same magnetization type, corresponding to different rings of the Halbach dipole, a prior distribution π 0 of p can be derived by computing the sample mean µ and sample covariance C 0 , leading to p ∼ π 0 := N (µ 0 , C 0 ), where N is the Gaussian distribution.The Anderson-Darling test justifies the usage of the Gaussian distribution.
The simulation model H : p → q maps the parameter vector p to a field-related measurable quantity q such as the magnetic flux density B in the bore of the magnet.The two definitions of q that are used in this article are where each B ⋆ depends on measurement positions x, y, z, and where A k (z) and B k (z) are the Fourier coefficients of the magnetic flux density component B r (z) on a circle with radius r 0 = 75 mm, centered in the magnet bore.The advantage of working with the Fourier coefficients is the averaging effect on random noise.However, especially in the 2-D case, considering only the Fourier coefficient leads to an under-determined inverse problem.
If the reluctivity in D is assumed to be constant, for example, the iron ring is neglected D iron = ∅, the model is linear and we can write Hp = q by abuse of notation.
To select positions x, y, at which the magnetic flux density B is sensitive to variations of p, a sensitivity analysis is performed.Similar to [10], it can be shown that the Gateaux derivative B ′ of the mapping can be determined by finding the weak solution and computing B ′ = curl A ′ .For illustration, Fig. 1 shows the absolute value of the Gateaux derivative It increases toward the segment whose magnetization is varied.Thus, equally distributed measurement positions on a centered circle inside the inner air region are chosen.

III. BAYESIAN INFERENCE
Due to uncertainties of the magnetic flux density measurements, the measurable quantity q = H(p) is not directly observable, but q obs := H(p) + ε with ε ∼ N (0, ).Thus, the distribution of q obs given p is π(q obs |p) = N (H(p), ).Following [11], Bayes' rule: can be applied to determine the posterior distribution of the magnetization parameter vector p, given an observation q obs of the measurable quantity.To generate samples from the posterior distribution, we use different approaches depending on whether the model H is linear or nonlinear.
In the linear case, the posterior distribution is proportional to a Gaussian distribution [12] π(p|q obs ) ∝ N (µ 1 , C 1 ).
The expected value µ 1 and the covariance matrix C 1 of this distribution are given by In the nonlinear case, the Metropolis-Hastings (MH) algorithm [11] is applied.The principle of this algorithm is the construction of a Markov chain, whose states are samples of π(p|q obs ).Therefore, the transition kernel K , for which π(p|q obs ) is an invariant distribution, is defined by with the proposal density π ( p|p) and a( p, p) := min 1, π( p|q obs ) π (p| p) π(p|q obs ) π ( p|p) .
The preconditioned Crank Nicolson proposal [13] π ( p|p Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. with step size s = 1/80 is chosen, which incorporates the prior knowledge on the covariance C 0 .The sample {p k } k≥0 is then generated by conducting the following steps of the MH sampling [11]. 1) Initialize p 0 := µ 0 .Set k = 1.
2) Generate the proposed sample p ∼ π ( p|p k−1 ).Accept p by setting p k := p with probability a( p, p k−1 ).Else, set Inserting the Gaussian distributions of the prior, the measurement uncertainty, and the preconditioned Crank Nicolson proposal in (13) simplifies the acceptance probability to

IV. VALIDATION
For validation of the algorithm, observation data is generated based on a ground-truth parameter vector p true ∼ π 0 by setting For q obs B based on magnetic flux densities, we set σ = 10 −4 T, for the observation of Fourier coefficients q obs F , we set σ = 10 −6 T, which are reasonable values for modern-day measurement equipment.
The simulation models in 2-D and 3-D are implemented using the FEM software getDP [9], and the input files to the numerical models are available at [14].
In the linear case, the algorithms are validated on the 3-D simulation model.The ground-truth p true is compared to the prior distribution π 0 and the two posterior distributions that are based on the different observation vectors.In Fig. 2, the ground-truth of the magnetization deviation p true − M restricted to the fifth ring of the Halbach array is shown together with the expected values and the variances.The maximal difference between µ B 1 and p true is decreased by 70% compared to the maximal difference between µ 0 and p true .The variances are also smaller.
In the nonlinear case, the method is applied to the 2-D simulation model and simulated observations q obs B of magnetic flux densities.To obtain the sample {p k } k , the MH sampling is applied for k = 18 000 steps.In Fig. 3, the sample mean and the sample variance of the prior and posterior distributions are plotted with the ground-truth p true .The posterior sample mean follows the ground-truth better than the prior sample mean, the maximal difference is decreased by 50%.
Applying the MH algorithm in combination with the 3-D nonlinear simulation model is computationally prohibitive, because the simulation is more time-consuming, and a higher number k of samples is required because of the increasing dimension of the parameter vector in 3-D.

V. APPLICATION
The algorithm in the linear case (neglecting the iron ring) is applied to observations q obs F of K = 8 Fourier coefficients on a centered circle with r 0 = 75 mm radius in dim(z) = 156 positions along the full range of the magnet, including the fringe field.Measurements of the magnetic flux density are taken with a Hall probe mapper for each coordinate of z in 60 equally distributed points on the circle and q obs F is obtained by Fourier analysis.
For estimating the noise covariance matrix positioning uncertainties that result from small misalignment of the various coordinate systems (magnet, mapper, and simulation), as well as small oscillations of the Hall probe mapper system have to be taken into account.In the fringe field region, where the magnetic flux density changes the most, positioning errors lead to larger measurement errors than in regions where the magnetic flux density is almost constant.Therefore, = σ (z) 2 I is chosen position-dependent with σ (z) = 5 × 10 −5 and 5 × 10 −3 T in the homogeneous field region, and σ (z) = 5 × 10 −5 and 5 × 10 −3 T in the fringe field.Fig. 4 shows the relative error Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.between the measured and simulated B x component along the z-axis for the simulation with the prior parameter vector µ 0 and the updated posterior parameter vector µ F 1 .The measured flux density data used in this comparison was not part of the training set used for updating the model.The relative error E rel (z, µ F 1 ) of the update is one order of magnitude smaller than the relative error E rel (z, µ 0 ) of the prior, almost everywhere in the homogeneous field region.

VI. CONCLUSION
Bayesian inference was applied to combine domain knowledge and observations from material and magnetic flux density measurements for deriving a posterior magnetization distribution of the PM blocks of a Halbach array.The method is not only validated on simulated data, but also applied on measurements of a Halbach dipole of the FASER experiment.Updating the magnetization of the PM blocks with the expected value of the posterior distribution decreases the relative error of the simulated magnetic flux density compared to the measured magnetic flux density inside the aperture by one order of magnitude.A more detailed analysis of the systematic measurement errors will be investigated in future research.Furthermore, surrogate models of the nonlinear 3-D model shall be investigated to reduce the computational cost of the MH algorithm.

ACKNOWLEDGMENT
This work was supported in by the Wolfgang Gentner Programme of the German Federal Ministry of Edu-cation and Research under Grant 13E18CHA and in part the Graduate School Computational Engineering (CE) within the Centre for Computational Engineering, Technische Universität The authors would like to thank Erik Schnaubelt for helping with the getDP models and Olaf Dunkel and Mariano Pentella for the PM block measurements.

Fig. 1 .
Fig. 1.Top left: Domain D with the cross section of the FASER Halbach dipole with PM blocks D i and their nominal magnetization M i and nominal orientation α i .Top right: Absolute value of the Gateaux derivative B ′ of the mapping (7).Bottom: 3-D geometry of the FASER Halbach dipole with D iron = ∅.region one order of magnitude better than the prior simulation model.II.PROBLEM STATEMENT The magnetostatic field problem of a Halbach array is defined on a domain D = D air ∪ D iron ∪ D m , consisting of the air region D air , iron region D iron , and permanent magnet region D m (see Fig. 1).The fields are defined by curl H = 0 in D, divB = 0 in D, and B•n = 0 on ∂ D, with magnetic flux density B, magnetic field strength H, and outward pointing unit normal vector n.The corresponding constitutive equation is ν(∥B∥)B = H + M, where M denotes the magnetization with supp(M) = D m = 16 i=1 D i and D i denotes the PM blocks (Fig. 1 top, left).The reluctivity function ν : R + 0 → R + 0 is defined by ν(s) = 1/µ 0 in D air and ν(s) = 1/(µ 0 µ r ) in D m

Fig. 2 .
Fig. 2. Validation of posterior derivation algorithm in the linear case on the 3-D simulation model.Comparison of ground-truth (red), prior N (µ 0 , C 0 ) (black), posterior based on magnetic flux density observation N (µ B 1 , C B 1 ) (blue) and posterior based on Fourier coefficient observation N (µ F 1 , C F 1 ) (green).

Fig. 3 .
Fig. 3. Validation of posterior derivation algorithm in the nonlinear case on the 2-D simulation model.Comparison of ground-truth (red), prior N (µ 0 , C 0 ) (black) and posterior (blue) sample mean and variance.

Fig. 4 .
Fig. 4. Relative error of prior and posterior simulation model compared to magnetic flux density measurements of the first FASER dipole.Area of the fringe field marked in gray.