BEM-Based Magnetic Field Reconstruction by Ensemble Kálmán Filtering

2.1 Discretization

In boundary-element methods (BEM), the boundary data is approximated by locally supported basis functions on a boundary mesh. The algorithms presented in this paper have been implemented in C++ using the Boundary Element Based Engineering Library (BEMBEL) [7]. This enables a discretization of ∂ ⁡ Ω by means of non-uniform rational basis spline (NURBS) parameterizations as proposed in [9]. NURBS parameterizations are commonly used in computer-aided design and enable the exact representation of a large class of geometries. For an iso-geometric approach the boundary data ν is approximated in the finite-dimensional space 𝒱 h ⊂ H 1 2 ⁢ ( ∂ ⁡ Ω ) , which is spanned by patchwise tensor product basis-splines of degree p (see for example [3]). This yields the approximated boundary data ν h and surface current density 𝒔 h in terms of the N degrees of freedom ν k . A tensor product basis-spline of degree p = 3 is shown in Figure 1 (right), together with the resulting surface current density for some given ν over sixteen elements. In the following, the vector of degrees of freedom for the discretized double-layer potential will be denoted as the state vector 𝝂 = ( ν 1 , … , ν N ) T ∈ ℝ N .

The normal component of the magnetic flux density 𝒏 ⋅ 𝑩 is proportional to the Neumann trace, i.e., 𝒏 ⋅ 𝑩 = - μ 0 ⁢ 𝒏 ⋅ grad ⁡ ϕ m , where μ 0 = 4 ⁢ π ⋅ 10 - 7 ⁢ Vs / Am is the vacuum permeability. For Lipschitz domains this holds almost everywhere on ∂ ⁡ Ω .

The component 𝒏 ⋅ 𝑩 provides the Neumann data for the computation of the state vector 𝝂 . To this end, we apply the Neumann trace operator to the double layer potential in (2.3), which results in the hypersingular integral operator, denoted as D : H 1 2 ⁢ ( ∂ ⁡ Ω ) → H - 1 2 ⁢ ( ∂ ⁡ Ω ) (see [35, p. 6]). Following a Galerkin discretization scheme for D, we derive the linear equation system

where

The k-th tensor product basis-splines is denoted by φ k . The explicit definition of these functions in the iso-geometric framework is given in [3]. Practically we use the representation

where

is the vectorial surface curl and φ ~ is extending φ locally into the neighborhood of ∂ ⁡ Ω (see [35, Section 1.1]).

The regularization term 〈 φ k , φ l 〉 ∂ ⁡ Ω is introduced as a cure to the null space of D (see [35]) such that [ 𝑫 ] k , l is positive definite.

Figure 1

Left: A domain Ω, its boundary ∂ ⁡ Ω and the outward directed unit vector 𝒏 in the air gap of an electromagnet. Right: A tensor product basis spline of degree p = 3 , defined over sixteen square elements. The vector field illustrates the resulting surface current density 𝒔 .

We can derive an equation for the evaluation of 𝑩 at the position 𝒓 , by taking the curl of (2.4)

where φ n is the n-th tensor product basis spline on the boundary mesh.

The Aubin–Nitsche formalism allows to derive point-wise error estimates for the evaluation of the discretized boundary integral equations (see [45, Theorem 12.8], [38, Example 4.2.15]). As it will be relevant for the upcoming discussions, we provide the result for the indirect Neumann problem in Theorem 1.

The constant c ⁢ ( 𝒓 ) decreases rapidly with increasing distance of 𝒓 to the boundary. This is known as the smoothing property in boundary element methods.

2.2 A Local Expansion for Particle Tracking

In particle beam dynamics the motion of a relativistic particle is usually expressed in so called phase-space coordinates. The momentum variables in phase-space coordinates are more abstract and differ from the ones used in classical mechanics, but this has the advantage that certain invariants of motion are conserved for a passage of the particle through a static magnetic field (see [47, Chapter 2]).

From the Euler–Lagrange formalism, the equations of motion result in Hamilton’s equations

where 𝒛 := ( 𝒒 T , 𝒑 T ) T expresses the phase space coordinates, with 𝒒 denoting the particles coordinate vector and 𝒑 its canonical momentum. In the absence an electric field the Hamiltonian H is defined as

where c denotes the velocity of light, q the particle charge and m its mass. Hamilton’s equations (2.11) may be formulated as an ordinary differential equation system of the type

Its solution by numerical integration for a given initial condition is the main objective of single particle tracking.^[1]

Due to the partial derivatives ∂ ⁡ H ∂ ⁡ 𝒒 , the spatial derivatives of the magnetic vector potential ∂ i ⁡ A j for i , j ∈ x , y , z are involved on the right-hand side of (2.13). A useful field model must therefore provide 𝑨 , as well as its spatial derivatives of first order. Approaches for numerical integration of (2.13) may use explicit or implicit higher order Runge–Kutta schemes [19] or also Lie-algebraic methods [48].

In principle one could apply the representation formula (2.4) directly in the numerical integration scheme. The spatial derivatives can be computed by differentiating the Green’s function u * analytically and evaluating the resulting integral equation numerically. This, however, requires the repeated integration over the whole boundary ∂ ⁡ Ω for each particle at each position, which is not known a-priori.

More efficient approaches have been proposed, using local expansions of the magnetic vector potential around a reference trajectory by either Taylor approximations [11] or the fitting of 𝑨 by solid harmonic functions within intercepting spheres [2].

Using Taylor expansions makes it necessary to impose (2.1) to recover a magneto-static solution. This is implied by using the solid harmonic functions. We therefore follow the approach developed in [2], but in contrast to fitting 𝑨 on the surfaces of spheres, we make use of a local expansion of the Green’s function according to Theorem 2 for the computation of the local harmonic expansion. This allows us to derive analytical bounds for the truncation errors and also the fast calculation of the local expansion for an ensemble of state vectors.

The functions r l ⁢ Y l m and Y l m / r ′ ⁣ l + 1 are the regular and irregular solid harmonic functions, which are eigensolutions of the Laplace equation.

By Theorem 2, the magnetic vector potential may be expanded at the position 𝒓 j , according to

This expansion is valid inside the largest sphere with center 𝒓 j , which does not enclose any part of the domain boundary.

The evaluation position 𝒓 - 𝒓 j is expressed in spherical coordinates ( r , θ , φ ) and the moments 𝒎 l m ⁢ ( 𝒓 j ) are given by

where 𝒓 ′ - 𝒓 j = ( r ′ , θ ′ , φ ′ ) . It is possible to assemble the integral (2.16) for the basis functions used for the approximation of ν h . In this way, the moments can be evaluated rapidly by matrix vector products for an ensemble of state vectors.

The moments 𝒎 l m ⁢ ( 𝒓 j ) can be computed for J points 𝒓 j along a reference trajectory. The particle tracking can then be performed on (2.15) and it requires only the evaluations of the solid harmonics r l ⁢ Y l m ⁢ ( θ , φ ) . As the solid harmonics are eigensolutions of the Laplace equation, evaluating the truncated sum in (2.15) yields an exact magneto-static solution.

The truncation error for the evaluation of 𝑨 in Ω can be controlled by means of the truncation parameter L as well as the number of expansion positions J. This is conceptually equivalent to p and h refinement in boundary and finite element methods.

Figure 2

Truncation error for the evaluation of the magnetic vector potential over the ratio r / R . Dashed: Theoretical bound. Solid: Observed error in a numerical experiment.

Figure 2 shows the theoretical bounds for the truncation error as a function of the ratio r R and also the truncation errors observed in a numerical experiment. The rates of the theoretical bounds meet our observations, but the absolute values are pessimistic. However, this is uncritical, since the truncation error | 𝑨 ⁢ ( 𝒓 ) - 𝑨 L ⁢ ( 𝒓 ) | can be computed for a given L. In this way L can be adjusted for an adaptive expansion.

3.1 The Observation Operator

The voltages measured by Hall sensors can have complicated dependencies on the three components of the magnetic flux density. These stem from planar Hall effects and geometric imperfections, but also from non-linearity errors of the sensor transfer function. Calibration techniques allow for the characterization of planar and axial Hall effects by rotating the sensor in a reference field. This gives rise to the sensitivity function

with the complex calibration coefficients c i l , m ∈ ℂ , where c i l , m = - ( c i l , - m ) * , and the spherical harmonic functions Y l m . The angles θ ⁢ ( 𝑩 ) and φ ⁢ ( 𝑩 ) denote the polar and azimuth angles of the 𝑩 field in a sensor coordinate system (see Figure 3 (a)). We have introduced the index i ∈ { 0 , 1 , 2 } to distinguish the three Hall sensors involved in a three component measurement. The sensors do not need to be aligned with any of the magnet coordinate axes. Their misalignment is encoded in the coefficients c i l , m . Moreover, the three sensors cannot be mounted at exactly the same spot. For this reason, the positions of the three Hall sensors are denoted by 𝒓 m + 𝒐 i , where 𝒓 m is a measurement position (usually the position measured by the linear encoders) and the vectors 𝒐 i denote the displacements of the Hall sensors with respect to this position. Typically, the displacements between the sensors are in the mm range. The absolute sensor position and orientation can be identified by measurements in calibration magnets and is therefore assumed to be given.

To establish the observation operator, we define the three Hall voltages by

and sort the measurement vector according to the scheme 𝒚 ~ = ( U 0 ⁢ ( 𝒓 1 ) , … , U 0 ⁢ ( 𝒓 M ) , U 1 ⁢ ( 𝒓 1 ) , … , U 2 ⁢ ( 𝒓 1 ) , … ) T , considering M measurement positions. The m-th row of the observation operator 𝑯 : ℝ N → ℝ 3 ⁢ M × N is computed by

for i ∈ { 0 , 1 , 2 } and 𝑩 ⁢ ( 𝒓 , 𝝂 ) according to (2.9).

3.2 The Noise Model

For the uncertainty quantification, a noise model is needed which represents the real measurement data. We therefore follow the guide to the expression of uncertainty in measurement (GUM) [49]. It proposes to introduce input quantities to the measurement equation (3.3), in order to model the relevant physical phenomena, which are perturbing the measurements. Then, the uncertainty in the measurands is quantified by means of a first order Taylor approximation of the measurement equation with respect to the input quantities. This procedure is justified as long as the values of the input quantities are small deviations from zero.

In case of the Hall probe mapper system, the measurements are mostly affected by perturbations of the sensor position and orientation, due to the vibration of the mapper arm. We therefore include five degrees of freedom for each measurement position as input quantities. Three displacements, denoted as w x , w y and w z , and two rotations φ x and φ y ^[2]. The input quantities are illustrated in Figure 3. All input quantities are collected into the 5 ⁢ M -dimensional vector of mechanical perturbations 𝒅 .

Figure 3

Left: The Hall probe mapper system. It uses the stages of a coordinate measuring machine (C) to manipulate the position of a three axes Hall sensor (A) within an accelerator magnet. The sensor is mounted on an arm (B). The stage position is measured by linear encoders (D). Right: (a) A Hall sensor and a local coordinate system ( u , v , w ) . The angles θ and φ are the polar and azimuth angles in the ( u , v , w ) coordinates. (b) Input quantities for modeling the mechanical perturbations of the sensor position. These perturbations occur due to position errors and arm vibrations.

The mechanical perturbations are caused by random perturbations of the stage motion and are therefore not deterministic. However, the vibrations of the mapper arm lead to correlations between the measuring positions. These correlations have been determined by optical measurements in the machine setup. This yields a statistical model for the mechanical perturbations, which is of the type 𝒅 ∼ 𝒩 ⁢ ( 𝟎 , 𝚺 𝒅 ) , with a covariance matrix 𝚺 𝒅 .

As proposed in the GUM, the impact on the measurement vector 𝒚 is determined by means of a first order Taylor approximation. To this end, we introduce the rotation matrix 𝑴 ⁢ ( φ x , φ y ) ∈ ℝ 3 × 3 , which describes the rotations around the transversal axis^[3]. We then substitute 𝑩 ↦ 𝑴 ⁢ ( φ x , φ y ) ⋅ 𝑩 and 𝒐 i ↦ 𝒐 i + ( w x , w y , w z ) T in (3.3) and derive the resulting expressions with respect to the input parameters w x , w y , w z , φ x and φ y . This yields the five sensitivity matrices ∂ ⁡ 𝑯 j ∈ ℝ M × N for j ∈ { w x , w y , w z , φ x , φ y } . Collecting the vectors ∂ ⁡ 𝑯 j ⁢ 𝝂 in the matrix ∂ ⁡ 𝑯 ⁢ ( 𝝂 ) = ( diag ⁢ ( ∂ ⁡ 𝑯 w x ⁢ 𝝂 ) , … , diag ⁢ ( ∂ ⁡ 𝑯 φ y ⁢ 𝝂 ) ) ∈ ℝ M × 5 ⁢ M , the Taylor approximation can be denoted by

The matrices ∂ ⁡ 𝑯 j are fully populated and might occupy a large memory. For the estimation of the measurement noise, the computation of the products ∂ ⁡ 𝑯 j ⁢ 𝝂 with an accuracy in the percent range is tolerable. To this end, the matrix vector products can be computed on low rank approximation saving a considerable amount of memory. In this context fast multipole methods are particularly effective, since all five matrix vector products are using the same upward and downward pass operations [17].

5.1 The Kálmán Filter

The update step of a linear Kálmán filter is derived from the posterior under the assumption of a linear observation operator 𝑯 and normal distributions for likelihood and prior. Denoting the prior pdf by p ⁢ ( 𝝂 ) ∼ 𝒩 ⁢ ( 𝝂 0 , 𝑸 0 ) , the posterior takes on the form p ( 𝝂 | 𝒚 ) ∼ 𝒩 ( 𝝂 1 , 𝑸 1 ) , with the posterior covariance matrix

and the expected value

A well-established alternative representation for 𝝂 1 and 𝑸 1 is

where 𝑰 N × N ∈ ℝ N × N is an identity matrix and 𝑲 ∈ ℝ N × N is the Kálmán gain

This follows from the Sherman–Morrison–Woodbury formula [40].

Both formulations require a matrix inversion to compute the posterior covariance matrix 𝑸 1 . There are circumstances where this inversion is problematic. This for instance is the case in a high-dimensional setting, where 𝑯 can only be stored in a compressed format, such as a fast multipole matrix. Other compression possibilities would involve H/H2 matrices [9] or tensor decompositions [22]. Moreover, in some circumstances, the observation operator 𝑯 might require the inversion of a matrix itself (see Remark 7), although this is not the case for the indirect double-layer potential formulation according to (2.3). For these reasons it can be of advantage to use Theorem 4 to compute an ensemble of K state vectors 𝝂 k and estimate the covariance matrix from these samples.

Approximating the posterior covariance matrix by an ensemble of K state vectors with K ≪ N can be considered as a compression technique to reduce the computational effort significantly. This advantage is further enhanced by the use of fast algorithms for the computations of the products 𝑯 ⁢ 𝝂 , e.g., the fast multipole method.

Figure 6

Elements of the posterior covariance matrix in logarithmic scale. See Section 4 for details about the experiment. Here M = 4093 measurement positions were collected. Left: The covariance matrix was computed from direct inversion of the matrix 𝑸 1 - 1 using the LU-solver implemented in Eigen [18]. Center and right: The covariance matrix was approximated from K = 7000 and 1000 samples 𝝂 k . The approximation is visible by sampling errors, i.e. spurious covariances.

Figure 7

Uncertainty in the field reconstruction along a reference trajectory s, measured in three standard deviations. The reconstruction is based on the covariance matrices given in Figure 6. Top: (blue) The uncertainty is evaluated by the normal transformation 𝑯 ⁢ 𝑸 1 ⁢ 𝑯 T , where 𝑸 1 was computed by direct inversion and 𝑯 is the observation operator for B y . (orange and green): The covariance matrix is approximated from K = 7000 and 1000 samples. It is found that the approximation of the covariance matrix according to Figure 6 is not harmful for the uncertainty quantification at the reference trajectory. This is due to the smoothing property of the boundary integral equation. Bottom: The red line shows the mean value of B y ⁢ ( s ) . The red contour illustrates the uncertainty, whereas, for better visualization, the standard deviation was multiplied by a factor 300.

In Figure 6, we illustrate the posterior covariance matrix approximated for different ensemble sizes K and also the case of a direct inversion of 𝑸 1 - 1 . In this case M = 4093 measurement positions were taken with the Hall probe mapper system. The finite sampling size introduces sampling errors, which appear as spurious correlations. This is a well-known phenomenon in ensemble methods (see [14]). However it is not harmful for the field evaluation inside the domain, particularly in the vicinity of the particle beam, because details in the matrix 𝑸 1 are smoothed out when evaluating the boundary integral equation. This is shown in Figure 7, where the uncertainty in the evaluation of B y evaluated along a reference trajectory is shown.

5.2 The Ensemble Kálmán Filter

The equations for the Kálmán filter require the inverse of the prior covariance matrix 𝑸 0 to be accessible. In a high-dimensional setting they are therefore only applicable for simple prior covariance models, such as identity or sparse matrices, but also the covariance matrix which will be presented in Section 5.3.

In the case of the ensemble Kálmán filter, the prior and posterior covariance matrices are represented by ensembles of state vectors. In this way a significant reduction of computational complexity can be achieved.

Similarly to (5.7) the ensemble Kálmán filter is based on a randomized linear equation system. The proof for the following statement is completely analogous to the one provided for Theorem 4 and it is found in [13].

The following equations summarize the analysis step of the ensemble Kálmán filter as they are well known in the field of data assimilation (see for instance [4], [13] or [14]). The ensemble Kálmán filter follows from Proposition 2 after estimating the mean and covariance of the prior 𝝂 0 according to

The update step can then be performed directly on the ensemble 𝑽 0 := ( 𝝂 0 1 , … , 𝝂 0 K ) ∈ ℝ N × K and the data matrix 𝒀 := ( 𝒚 1 , … , 𝒚 K ) ∈ ℝ 3 ⁢ M × K , via the equations

The matrix 𝑽 1 = ( 𝝂 1 1 , … , 𝝂 1 K ) stores the posterior ensemble. The statistics of the posterior distribution can be computed from the sampled mean and covariance of the posterior ensemble (see (5.9)).

In the case of a non-linear observation operator 𝑯 , or when the observation operator is available in compressed form only, a matrix-free implementation is possible. One therefore defines the matrices 𝑯 ⁢ 𝑽 0 and 𝑯 ⁢ 𝑼 according to

where

In this work the ensemble Kálmán filter is used to update the prior ensemble 𝑽 0 in regions of the magnet where the remaining uncertainty is large. Compared to the dimension of the state vector N, less measurements are used for such local updates. The 𝑷 matrix is therefore of a size M ≪ N . For high-dimensional updates M ≫ N , it is shown in [30] how to reformulate 𝑷 to obtain a linear complexity in M.

5.3 The Prior

Due to the smoothing property of the boundary integral equation, the inverse problem is ill-conditioned. To ensure a numerically stable solution we derive a regularization term which is based on a prior field solution. It is now shown how mean 𝝂 0 and covariance 𝑸 0 of a Gaussian prior p ⁢ ( 𝝂 ) ∼ 𝒩 ⁢ ( 𝝂 0 , 𝑸 0 ) can be computed from a numerical field simulation.

From (2.5) we derive the randomized indirect Neumann problem:

The process noise ϵ p is used to impose uncertainty in the field simulation. From the linearity of the expected value function it follows

due to the symmetry of 𝑫 .

Equation (5.15) provides a Gaussian prior p ⁢ ( 𝝂 ) , with covariance matrix 𝑸 0 = δ - 1 ⁢ 𝑫 - 1 ⁢ 𝑫 - 1 . In this context, the discretized hypersingular operator 𝑫 already provides the square root of matrix 𝑸 0 - 1 . It therefore ties in nicely with equation (5.7). As an alternative to (5.7) it is possible to solve (5.14) for K realizations of ϵ p k to generate an initial ensemble and use it for the ensemble Kálmán filter.

The regularization parameter δ controls the impact of the prior. Since 𝑸 0 - 1 = δ ⁢ 𝑫 ⁢ 𝑫 , all terms related to the prior in (5.7) vanish for δ → 0 . It then holds a proportionality relation between the likelihood and the posterior. This prior is not informative and it may be considered as choosing a sufficiently flat pdf such that p ⁢ ( 𝝂 ) ≈ constant in Bayes rule (5.1).

On the other hand, the prior state vector 𝝂 0 is recovered in the limit δ → ∞ . This can be seen when dividing (5.7) by δ and substituting 𝑸 0 - 1 = δ ⁢ 𝑫 ⁢ 𝑫 such that

The parameter δ has the same role as the ridge parameter in the Tikhonov regularization. Different approaches have been developed for its selection (see for instance [1, Chapter 1]). In the following we make use of a set of validation measurements, providing us with the vector 𝒚 val ∈ ℝ 3 ⁢ M val and we select δ according to

where 𝝂 δ * is computed from (5.16) using δ = δ * , 𝑯 val is the observation operator for the validation set and E ⁢ ( ⋅ ) is the expected value function. By this definition ϵ is the root-mean-squared (RMS) error between the predicted mean and the validation measurements.

In Figure 8 (right), the validation error ϵ ⁢ ( δ * ) is shown in a logarithmic scale, for two scenarios. The results have been obtained from the application discussed in Section 4. Whereas only M = 4093 measurements have been collected in the first scenario, in the second one, we have taken measurements at M = 16 324 positions. The validation errors saturate at ( 3 ⁢ M ) - 1 2 ⁢ ∥ 𝒚 val - 𝑯 val ⁢ ( 𝝂 0 ) ∥ 2 for δ → ∞ in both cases. This is the RMS error of the prior field solution. On the other hand, the rank of the matrix 𝑯 T ⁢ 𝑹 - 1 ⁢ 𝑯 determines the validation errors for δ → 0 . For the case M = 4093 this matrix is rank-deficient, while for the second scenario, more data is available and ϵ ⁢ ( δ * ) saturates at the maximum likelihood solution. The values for δ * minimizing the validation errors are found in Table 2. A-priori estimates for δ that avoid the validation set are discussed in [1, Chapter 1].

Figure 8

Left: Prediction of | 𝑩 | along the x-axis in the magnet center, for different values of the regularization parameters, compared to a validation measurement, here M = 4093 . Right: The validation error over δ * in logarithmic scale for both scenarios.

5.4 Design of Experiment

Design of experiment is about determining the minimum amount of measurements required to meet the accuracy requirements. As we have seen in the previous section, the number of measurements, their spatial distribution and also the prior E ⁢ ( 𝝂 0 ) determine the matrix 𝑯 T ⁢ 𝑹 ⁢ ( E ⁢ ( 𝝂 0 ) ) - 1 ⁢ 𝑯 , which is equivalent to the inverse of the posterior covariance matrix for δ = 0 .

Using the procedure developed in the previous sections, we are capable to propagate the measurement uncertainty down to beam related quantities by particle tracking. In the following, we consider the trajectory of a proton with momentum P 0 = 100   MeV / c , for a single passage through the magnet as the quantity of interest.

This trajectory is illustrated in Figure 9, which shows the solution of (2.13) for a reference particle using an explicit, fourth-order Runge–Kutta scheme for the numerical integration (see [20, Chapter 2]).

We have calculated the trajectory for K = 1000 ensembles of the posterior, based on the ensemble generated from the scenario M = 16324 . The uncertainty in the fringe field region results in 3 standard deviations of 51 ⁢ μ ⁢ m for the prediction of the x-coordinate at the exit of the magnet. For comparison, a constant deviation in the vertical flux density B y of one unit in 10 000 with respect to the maximum field, gives a total deviation in the trajectory radius of 108 ⁢ μ ⁢ m . The uncertainty of 51 ⁢ μ ⁢ m is therefore acceptable for many applications.

In contrast to the posterior covariance matrix, the validation error cannot be determined a-priori. It mainly origins from locations in the fringe fields, where the spatial resolution of 6 ⁢ mm is not sufficient to properly resolve the differences between prior and posterior. As a mitigation measure, update measurements are taken at these positions. In order to avoid moving the stages back and forth between the left and the right-hand side of the magnet, both sides of the magnet are treated individually.

The initial ensemble is updated move-by-move, using the equations of the ensemble Kálmán filter with K = 1000 . We denote the posterior ensemble of update j by 𝑽 j . We have performed a total of J = 400 update measurements equally distributed between left and right-hand side. These updates provide us with the posterior ensembles 𝑽 200 and 𝑽 400 , which are the results after treating the left and right-hand side of the magnet.

Table 3 gives the resulting validation errors. By means of the ensemble Kálmán updates, the validation error was reduced by a factor 2.5. As a side effect, more measurement data increases the precision in the prediction of the particle trajectory. This is shown in Figure 9, where we have added the trajectories evaluated from the posterior ensembles 𝑽 200 and 𝑽 400 .

Figure 9

Prediction of the particle trajectory for different stages of the ensemble Kálmán filter. The deviations from the mean trajectory were amplified by a factor of 1000 for better visualization. The deviations are representing the uncertainty in the particle trajectory due to the measurement uncertainty in the fringe fields. First M = 16324 measurements were taken with the Hall probe mapper system and the trajectory was evaluated from K = 1000 samples from the posterior ( 𝑽 1 ). Update measurements in the left and right-hand side of the magnet (gray) were performed to reduce the validation error and also the precision in the particle trajectory. This yields the ensembles 𝑽 200 and 𝑽 400 . A proton with momentum P 0 = 100 MeV / c is considered. This corresponds to the injection energy of the ELENA ring.