Magnetic Position System Design Method Applied to Three-Axis Joystick Motion Tracking

This manuscript discusses the difficulties with magnetic position and orientation (MPO) system design and proposes a general method for finding optimal layouts. The formalism introduces a system quality measure through state separation and reduces the question “How to design an MPO system?” to a global optimization problem. The latter is then solved by combining differential evolution algorithms with magnet shape variation based on analytical computations of the field. The proposed formalism is then applied to study possible realizations of continuous three-axis joystick motion tracking, realized with just a single magnet and a single 3D magnetic field sensor. The computations show that this is possible when a specific design condition is fulfilled and that large state separations as high as 1mT/∘ can be achieved under realistic conditions. Finally, a comparison to state-of-the-art design methods is drawn, computation accuracy is reviewed critically, and an experimental validation is presented.


Introduction
Magnetic position and orientation (MPO) detection systems determine the relative motion between permanent magnets and magnetic field sensors by measuring the modulation of the magnetic field. Such systems offer many advantages like robustness against dirt and temperature, long lifetimes ensured by contactless operation, as well as high resolutions at low cost and low power operation [1,2]. Multiple applications are treated in the literature, including proximity detection, linear motion systems [3,4], angle and rotation sensing [5,6], encoders [7], and more complex forms like motion tracking of six degrees of freedom (DoFs) [8,9]. Nowadays, more than one hundred MPO system applications exist in the automotive sector alone, including gas and brake pedals, gear shifts, indicator levers, side mirror position, wheel speed, and anti-lock braking system (ABS) sensors [10].
A position system of specific interest is the three-axis joystick, which combines regular 2D joystick motion with rotation of the lever about its own axis. This concept is commonly employed for control elements in the automotive [11,12], nautical [13,14], medical [15], and aerospace [16] domains as well as for consumer electronics applications like arcade sticks [17] and closed-circuit television (CCTV) steering [18]. State-of-the-art magnetic implementations combine two MPO systems, namely a 2D joystick and a rotation sensor [11,13]. The latter must be integrated into the movable joystick shaft, thereby making it difficult to mechanically manufacture this type of implementation. Recent proposals show that special cases of three-axis joystick motion tracking can be realized with a limited number of DoFs in the form of an MPO system with only a single magnet and a single 3D sensor [19][20][21]. Thanks to their simplified construction, such devices are highly cost-efficient, which is a critical aspect for industrial applications, as is also shown by recent sensor and magnetic material developments [22,23].
Readout of an MPO system requires reconstruction of the magnet position from the sensor outputs. This is closely related to the magnetostatic inverse problem that is often mathematically ill-posed and computationally demanding [24,25]. One of the biggest challenges is to solve the inverse problem within milliseconds in order to enable real-time measurement and interaction. State-of-the-art systems address this problem by approximating the field with simple harmonics that can be easily inverted [3,4]. Several field-shaping proposals improve this approach by engineering fields with specific desired forms through multi-magnet arrangements [26], shimming techniques, or complex magnet shapes with inhomogeneous magnetization [27,28]. When more than two DoFs and more complex motions are involved, a direct numerical inversion of 3D approximations of the magnetic field seems to be the only viable solution. Such approximations are based on analytical solutions of permanent magnet problems [8], dipole approximations [9], pre-computed solutions, or simple look-up tables [21] to achieve the necessary computation times for fast inversion.
While sophisticated techniques exist to solve the inverse problem, there is no patented procedure on how to layout an MPO system, i.e., how to generally arrange magnets and sensors in order to realize the desired motion parameters of interest in the best possible way. For this, state-of-the-art implementations rely mostly on experience and educated guesses combined with point-wise finite element (FE) simulations for layout testing and optimization, an approach which is limited by the intrinsically long computation times involved. It is the aim of this paper to overcome these limitations by using computationally efficient analytical methods to find possible MPO system realizations and optimal layouts. The paper structure is as follows: In Section 2.1, a general formalism is introduced which describes the conditions for a feasible and optimal MPO system layout. This formalism is compared to state-of-the-art field shaping methods in Section 2.2. A description of the three-axis-joystick is then given in Section 2.3. Magnetic field computations based on analytical solutions are discussed in Section 2.4, and differential evolution is suggested in Section 2.5 as a suitable global optimization method. The formalism is then applied to the three-axis-joystick problem, demonstrating under which conditions a continuous three-axis motion can be realized (Section 3.1) and how to find optimal layouts (Section 3.2). Two realistic and optimized layouts are proposed in Section 3.3. Finally, a comparison with an experiment is performed in Section 3.4.

General Formalism
This section introduces a general formalism for the design and optimization of MPO sensor systems, starting with the following relevant quantities: • the observables of interest α ∈ P, where P is the parameter space of interest. Typical observables can be the lever position in a gear shift or the tilt angle of a joystick. The parameter space P then simply corresponds to the allowed range of mechanical motion.

•
the system design parameters s ∈ S describe a specific MPO system implementation attempting to realize α. The system parameters include all quantities that can be varied within an allowed system parameter space S in a design process. They include magnet and sensor choice and placement within the system, component geometries, or material parameters. The parameter range S can be for instance the result of limited construction space. • the sensor outputs B(α, s) ∈ B space that correspond to single or multiple components of the magnetic field at the sensor positions. Although linear sensing technology is considered in this work, the formalism can be easily extended to include arbitrary sensor transfer functions. Hereafter, the sensor output and the magnetic field will be treated the same, and therefore, B space will simply denote the set of all possible sensor outputs. • a set of constraints C(s) that must be fulfilled in the design process. They can for example describe weighted sensing resolutions, maximal cost limitations, or the influence of system fabrication tolerances and external stray fields.
The central goal of MPO system design and optimization is to understand how well the observables of interest α can be determined from the sensor outputs B for a given system s and how to optimally design such a system when subjected to a limited system parameter range S and constraints C. It is critical to understand that the proposed formalism holds generally for any MPO system as the fundamental position detection limitation is given only by the field(s) at the sensor position(s).
Any function composition f(B) including sensor transfer functions or combinations of different field components, like the arctan 2 (B eve , B odd ) scheme in linear position systems [3], cannot improve system performance that is fundamentally limited by the magnetic field state density with respect to the natural sensor noise. In the case of the linear position system, for example, the goal is only to ease interpretation of the sensor output signal: the field ratio reduces the strong airgap dependence and the arctan 2 function naturally stitches discontinuities and linearizes the output.
The fundamental requirement for a system s to realize the observable parameter space P is that each state α ∈ P is associated with a unique sensor output B s ∈ B space . In other words, the magnetic field at sensor B(α, s) must be a bijective map between P and B space . Invertibility of the magnetic field on P is guaranteed when B is smooth and locally invertible for every α ∈ P and when B space is simply connected. The smoothness of B can be assumed from the structure of the problem (smooth input variables and smooth fields of permanent magnets). Simple connectedness of B space , however, must always be checked [29]. A brief discussion of this is presented in Appendix A. In agreement with the inverse function theorem, local invertibility of a multivariate vector function is ensured when the Jacobian matrix has full rank. The Jacobian J is defined in the usual way: B is thus invertible on P for a specific s if under the assumption that α and B have the same dimension. In the more general case, when the dimension of B is larger than the one of α, meaning there are more sensor outputs than observables of interest, then B space is simply a manifold with similar dimension as α, embedded in a higher dimensional space, and the requirement (2) becomes det J(α, s) T J(α, s) = 0 ∀α ∈ P.
A system s is feasible (can in principle be realized) if Equation (3) can be fulfilled. Moreover, by means of Jacobian analysis, the density of states (DOS) denoted by D can be directly obtained as The inverse of the DOS corresponds to the state separation ∆p = D −1 . Large state separations make it easier to detect different system states and, as such, ∆p(α, s) provides an estimation of the quality of a specific MPO system implementation s. Hence, a quality factor Q can be introduced: as the minimal state separation (weakest link) for a given system s. Systems with large Q-factors are preferable from a technical point of view because their individual states are better separated-and therefore easier to detect-and distortions from external influences typically play a lesser role. For a given system parameter range S and set of constraints C, the best possible system s opt is thus obtained as a result of the following optimization problem: Equations (2)-(6) are the core instruments to answer the following questions: "How can an MPO system be realized?" and "What is the best possible realization?". While state-of-the-art MPO design relies on educated guesses on which magnet-sensor arrangement can realize the desired observables of interest, Equation (6) simply reformulates this task as a global optimization problem.
To better understand the formalism, Figure 1 shows a sketch of the state separation for three different system implementations s 1 , s 2 , and s 3 together with the parameter space P of interest. Figure 1a represents a system where Equation (3) is violated. Figure 1b shows an implementation which is in principle possible (i.e., (3) is satisfied) but where state separations are bad (Equation (6) not fulfilled). Finally, Figure 1c represents an optimal solution satisfying Equation (6) with large state separations inside P.  3 represents an optimal implementation that not only is feasible but also exhibits a large state separation.

Field Shaping and Shape Variation
The optimization problem (6) is reminiscent of inverse magnet design for field shaping [24,26]. There is however a crucial difference. While field shaping attempts to give the magnetic field a specific, favorable form (e.g., a linear component [26]) which can then be easily processed for readout, Equation (6) only aims to find the configuration with maximal state separation and relies as such on a more complex form of direct inversion for readout [8]. Field shaping thus requests preliminary knowledge of a target field, whereas the proposed MPO design method circumvents this requirement by simply looking at all possible solutions.
To solve the optimization problem (6), it seems nevertheless reasonable to extend field shaping techniques to MPO design. In this context, topology optimization using the adjoint method is able to address thousands of DoFs and to find optimal magnet forms ab initio [24,[30][31][32]. However, remarkable results can already be achieved by variation of simple magnet shapes, as proposed in [26]. Comparison and discussion of the two methods are given in Appendix B.
In this paper, an extended version of the shape variation approach is chosen instead of complex procedures like the adjoint method for several reasons. First, the simple magnet forms are commercially available and cheap, while the slightly better performing but very complex magnet shapes that result from the adjoint method are hard to obtain and expensive in fabrication. In addition, the development effort required by shape variation is relatively small, whereas it can be demanding to formulate the complete MPO problem so that it can be treated with the adjoint method, which requires a priori calculation of an analytical variational derivation of the cost function. Complications may arise, for instance, when non-differentiable functions like min or max are involved.
Finally, it is worth noting that field shaping and the proposed MPO design should not be viewed as opposing strategies but rather as synergetic approaches. Indeed, MPO design by shape variation is limited to fewer DoFs and focuses on answering the question of how to realize a motion, which could, in turn, hint at potential target magnetic field shapes suitable for more sophisticated field shaping methods dealing with more DoFs and complex magnet forms.

The Three-Axis Joystick System
Here, the general formalism developed in Section 2.1 is applied to describe a specific MPO system of interest: a three-axis joystick in which the lever can be continuously tilted in two directions and rotated about its own axis. The objective is to resolve this three-axis motion by using only a single magnet and a single 3D magnetic field sensor. Previous works [19][20][21] focused on the extreme cost-efficiency of such an implementation in comparison to state-of-the-art solutions [11,13]. However, they only concentrated on a much simpler implementation with discrete tilt directions and tilt angles.
A convenient representation of the observables of interest describing the three-axis motion is given by the three angles ψ, θ, and ϕ sketched in Figure 2a: ψ is the azimuth angle corresponding to the lever tilt direction, θ is the polar angle indicating the amplitude of tilt, and ϕ is the rotation angle tracking the lever rotation about its own axis. According to the general formalism, the observables of interest are thus α = (ψ, θ, ϕ). The parameter space P is defined by the corresponding allowed angle ranges, which-consistent with typical three-axis joystick motion-are chosen as ψ ∈ [0, 360] • , θ ∈ [0, θ max ], and ϕ ∈ [ϕ min , ϕ max ].
The MPO system is realized by fixing a permanent magnet at the bottom of the lever and by mounting a 3D magnetic field sensor directly below. This configuration implies that α and the sensor output B have the same dimension, which reduces the feasibility study to a simpler sign analysis through the use of Equation (2).
The magnet is defined in a local coordinate system (LCS) fixed to the lever which is denoted by barred symbols (x,ȳ,z) and coincides with the global coordinates (x, y, z) when α = 0. A set of critical system parameters s is given by the following geometrical and physical quantities: • the position of a 3D magnetic field sensor r s = (x s , y s , z s ). The sensor output is the magnetic field vector B. These definitions naturally introduce an additional pair of critical parameters characteristic for such MPO systems, i.e., the airgap g = z s −z m between magnet and sensor as well as the magnet distance from the center of tilt d CoT =z m − c/2. Figure 2a,b shows a schematic of the magnetic joystick with all the corresponding relevant system parameters. The lever, the sensor, and all the other system component materials are chosen to be nonmagnetic (stainless steel, plastics, and silicon). The influences of sensor noise, possible magnetic shielding, external stray fields, or imperfect magnetization are neglected in this study.
The observables of interest α = (ψ, θ, ϕ) allow to infer the magnet position r m and orientation e m i in the global coordinate system in terms of a rotation R(α) of the magnet positionr m and orientations e m i in the LCS, A derivation of the rotation matrix R is reported in Appendix C. Assuming that the magnetic field can be computed in the LCS asB(r), it is possible to determine the sensor output as: Equation (9) expresses the sensor output B in terms of both the observables of interest α and the system parameters s, which forms the basis for a further system analysis starting with Equation (1).

Magnetic Field Computation
Equation (9) requires the magnetic fieldB(r) to be computed in the LCS. To this end, several viable options exist, the most commonly employed one being the FE method, as FE environments are readily available from multiple commercial and noncommercial sources [33][34][35][36], though the long computation times involved make such numerical approaches unpractical for dealing with the global multivariate optimization problem (6) that lies at the core of this paper.
As discussed in Section 2.2, the use of analytical solutions of permanent magnet problems [37,38] is envisaged in this work. Specifically, the field of cuboid magnets can be brought to a closed form [39,40], which enables a specifically fast computation of the sensor output. Moreover, cuboid-shaped magnets are commercially available off-the-shelf and their position and orientation can be defined very precisely when they are integrated into a mechanical setup, as opposed to spherical or cylindrical forms.
For implementation of the analytical formulas, the Magpylib Python package [41] is used, which is specially designed for dealing with MPO systems by integrating complex motions like those expressed by Equations (7)- (9). The analytical formulas are fully tested, vectorized, and optimized to ensure computational efficiency, achieving sub-microsecond computation times for calculation of the sensor output on standard x86 CPUs. In addition, a minimal development effort is required, since Magpylib enables system implementation with only a few lines of code, as demonstrated in Appendix D.
The analytical solution provides an excellent approximation of the magnetic field of realistic modern magnets despite neglecting demagnetization effects. A detailed discussion on the validity of this approximation is provided in Appendix E. In general, the error is less than 1% when µ r < 1.05 (realistic for high-grade NdFeB, SmCo, or ferrite materials) and when the distance between sensor and magnet is of the same order or larger than the size of the magnet.
Finally, it must be noted that many MPO systems use soft magnetic materials as magnetic shields, flux guides, and concentrators or shimming elements [1,42]. However, not only can soft magnetic materials not be simulated analytically but also they bring a decisive disadvantage to MPO systems: they generate position-dependent external stray fields that cannot be compensated by differential measurement. The current trend in MPO systems is towards increased stray-field stability [6,43,44], and for this reason, soft magnetic components are generally avoided.

Optimization Algorithm
Equation (6) describes an optimization problem leading to the best possible MPO system layout s opt that is subject to the constraints C and that realizes the parameters of interest α. The number of layout parameters and, therefore, the difficulty to solve the problem depend strongly on the complexity of the system. A single cuboid magnet already features 12 DoFs through its position, orientation, magnetization, and dimensions alone. For the typical 10-50 critical DoFs in MPO systems (see also the discussion in Appendix B), the differential evolution (DE) algorithm [45] is an excellent choice to solve such an optimization problem.
DE is a population-based evolutionary algorithm for the optimization of continuous variables in multidimensional spaces. Similar to genetic algorithms, it relies on an iterative process where a population is evolved by mutation, crossover, and selection to improve each generation. In contrast to genetic algorithms, however, DE avoids the harmful effect of mutation by carrying it out before the selection process [46].
The DE algorithm is especially well-suited to the MPO optimization problem (6) for several reasons. Firstly, the objective function is complex and its derivatives cannot be easily calculated, thereby favoring the use of such a black box optimization. Secondly, there can be multiple local optima which require the application of a global treatment. Finally, constraints can be easily included by nature of the algorithm, which adds to every new generation-only-allowed solutions. DE for magnet shape variation relies heavily on the fast computation times provided by the analytical solutions proposed in Sections 2.2 and 2.4 due to the large populations required for multivariate global optimization. In terms of computational efficiency, this approach can make optimal use of the computational resources on x86 type processors, as different field evaluations are completely independent of each other. The multiple field evaluations necessary to compute one objective function solution (α-grid spanning the space P) can be performed on the single instruction multiple data (SIMD) modules using the vectorized code from Magpylib, while the multiple objective function solutions (different values of s) required for constructing a population can be generated in parallel on separate cores. The SciPy [47] implementation of the DE algorithm provides automated multiprocessing and is used in this paper to perform the optimizations.
A similar computational setup was successfully used in the past for calibration of an MPO system [21]. The proposed implementation allows for reproducible and convergent parameter variations with up to several tens of parameters without resorting to extreme computation resources or distributed computation. A variation involving 20 system parameters and 460 field evaluations in the objective function was performed within 56 minutes and 44 seconds on an Intel R Xeon R Scalable Processor "Skylake" Gold 6126 (2.60 GHz, 12-Core Socket 3647, 19.25MB L3 Cache) running on 12 cores and converging within 2365 generations with population sizes of 2000 (n pop = 100 in the algorithm). This corresponds to ∼639 field evaluations per millisecond, not counting the algorithm effort.

Feasibility Analysis
A system s is defined as feasible if it can theoretically solve a given task, i.e., if there is a one-to-one correspondence between the mechanical states of interest α ∈ P and the sensor output B. Mathematically, this is expressed by Equation (2) for implementation of the three-axis joystick proposed in Section 2.3, where α and B are of similar dimensions.
A dipole moment is used as a magnetic field source for this feasibility analysis instead of a cuboid magnet. Since the dipole moment is the fundamental entity in magnetism and any magnetization distribution can be constructed from it by superposition, this approach provides a basic physical insight and is easily extended to finite-sized magnets. The magnetic field at the sensor location r s generated by a dipole moment m placed at the position r m is where r(α, s) = r s (s) − r m (α, s) is the distance between the sensor and the magnetic dipole. The position r m and the orientation of the moment m are obtained for each mechanical state α through Equations (7) and (8), while the sensor position r s is fixed completely by the system parameters s. For a generic implementation, it is always possible to find a connected parameter space P that satisfies (2). One must only make sure that P does not cross a breakdown hypersurface, i.e., surfaces where ∆p = 0, as sketched in Figure 1 for two observables. Figure 3a displays such a breakdown surface for a system withr m = (0, 0, 0) mm, r s = (3, 0, 0) mm and and a magnetic moment µ 0m = (0, 1.25 · 10 5 , 0) mT · mm 3 that corresponds to the magnetic moment of a cube having 5-mm-long sides and a 1000 mT remanence field. Any feasible connected parameter space P must lie completely either above or below the breakdown surface. The limited choices for sensing regions for this specific implementation are immediately apparent: it is, for example, possible to realize a system addressing all tilt directions ψ ∈ [0, 360] • , but only when restricting tilt angle θ and rotation angle ϕ. The possibility to include all tilt directions can be understood through a projection of the breakdown surface along the ψ-direction: the 3D breakdown surface then becomes a 2D breakdown region which is shown in Figure 3b.
For a feasibility analysis including all angles, ψ, ϕ ∈ [0, 360] • , magnet and sensor positions can be reduced tor m = (x m , 0,z m ) and r s = (x s , 0, z s ) withx m , x s ≥ 0 without loss of generality due to the rotation symmetry. Fundamentally different behaviors are observed when the magnet displacement exceeds the sensor displacement (x m > x s ) or when the opposite is the case (x m < x s ). This is demonstrated in Figure 4, where four variations are displayed: • system type s 1 : sensor in the center, magnet displaced (x s = 0,x m > 0), • system type s 2 : sensor and magnet displaced, magnet further out (x s <x m ), • system type s 3 : sensor and magnet displaced, sensor further out (x s >x m ), • system type s 4 : magnet in center, sensor displaced (x s > 0,x m = 0).
Systems s 1 and s 4 are special cases of s 2 and s 3 , respectively, with specific technical relevance.   Figure 4 shows that the shape of the sensing region strongly depends on the implementation s. Specifically, the sensing region becomes maximal in size when the magnetic moment is perpendicular to the plane spanned by the sensor position and the lever axis for α = 0 (blue sensing regions): m ⊥ (r s × e z ). (11) It is interesting to observe that full rotation ϕ ∈ [0, 360] • can only be realized when the magnet displacement exceeds the sensor displacement (x m > x s ) and while the perpendicularity condition (11) is simultaneously fulfilled. At the transition fromx m > x s tox m < x s (see Figure 4f,g), the large connected sensing region splits up into two disconnected lobes with opposite signs of the Jacobian determinant.
To better understand the breakdown for different implementations, the magnetic field of the two systems s 1 and s 4 is displayed in Figure 5  This forms the basis of the Mini-Drive implementation [19,21], where only 4 discrete tilt directions are considered.
As a result of a systematic investigation, a general empirical rule is proposed to determine the limits of the sensing region for s 1 and s 4 under the assumption that the perpendicularity condition is fulfilled. The maximum tilt angle θ max can be estimated as The values of θ max predicted by means of Equation (12) are reported in Figure 4e,h.

Quality Analysis
The feasibility analysis specifies which systems are in principle possible but provides no information about the quality, which is expressed in (4) by the DOS and the state separation ∆p. A specific implementation s is considered to be of high quality when all its mechanical states α ∈ P are well-separated in magnetic space, as expressed by Equation (5), which ensures easy state identification by the sensor.
For a further analysis of possible implementations, the state separations of the system types s 1 , s 2 , s 3 , s 4 defined in Section 3.1 are computed, assuming that the perpendicularity condition (11) is valid. The surprising results are shown in Figure 6a-d, where the sensing region is delimited in blue (in accordance with Figure 4) and the shading corresponds to the state separation. All systems exhibit low state separations for very small tilt angles. This is the result of the chosen coordinate representation, where the mechanical state density tends to infinity when θ → 0. However, the singularity is typically avoided by the dead-band θ dead of several degrees, which typically limits position computation to θ > θ dead in most applications [13].
When the sensor lies in the center and the magnet is displaced (system s 1 ), the rotation symmetry makes the DOS independent of the rotation angle ϕ, which is optimal for the realization of 360 • rotation. However, the state separation decays quickly for increasing tilt angles, as the distance between magnet and sensor increases, which makes it difficult to exploit the large θ max provided by the feasibility study above.
For small sensor displacement, the system makes a transition to s 2 . The state separation then becomes inhomogeneous in the rotation direction. Rotation angles about ϕ = 0 • can be resolved better at the expense of a lower state separation at ϕ = 180 • . This is opposite to the variation of θ max , which becomes smaller around ϕ = 0 • but increases at ϕ = 180 • .
When the sensor and magnet are at the same position, there is only a single feasibility lobe with maximum at ϕ = 180 • and complete breakdown at ϕ = 0 • . Then, by displacing the sensor beyond the magnet, system s 3 is realized. The large central lobe reduces in size and a second small lobe appears at ϕ = 0 • . Unfortunately, the large lobe-which is good, from a feasibility point of view, for realizing large ranges of the parameters of interest-features a low state separation, while the small lobe exhibits a large one. Both lobes become of equal size with similar state separation when the magnet is located in the center, thereby realizing the system s 4 .
In summary, an optimal system requires large state separation in the parameter region of interest. Small and inhomogeneous state separations reduce the quality of sensing regions. As shown by the above study, the low quality of the large lobe in s 3 and the inhomogeneity of the state separation in system s 2 make these implementations inferior in a technical sense to the implementations s 1 and s 4 .

Optimized Systems with Cuboid Magnets
The two most relevant system types s 1 and s 4 , introduced in Section 3.1, are here optimized in accordance with Equation (6) using the realistic parameters and finite-sized magnets outlined in Section 2.3. The goal of the optimization is to determine the best possible set of realistic system parameters s opt for given observable ranges of interest.
In both systems, the perpendicularity condition (11)  As expected, for the two optimized systems s 1,opt and s 4,opt the optimization procedure yields the minimum allowed value for the gap g = 2 mm and the maximum value of the distance of the magnet from the center of tilt d CoT = 4.97 mm: these correspond to small magnet-sensor distances and maximal mechanical state separation. The remaining displacements lead to a maximal state separation whenx m = 1.52 mm and x s = 3.26 mm for s 1,opt and s 4,opt , respectively. Figure 7 displays the state separation ∆p for the optimal systems. Panels (a) and (b) show that ∆p(s 1,opt ) is generally larger than ∆p(s 4,opt ) for these implementations with quality factors, computed via Equation (5), of Q(s 1,opt ) = 0.64 mT and Q(s 4,opt ) = 0.23 mT, respectively. The lower quality of s 4,opt is a result of the much larger lateral displacement which leads to a greater distance between magnet and sensor.

Experimental Results
In this section, the theoretical predictions are tested in an experiment realizing a system of type s 4 with observables of interest θ ∈ [2,15] • and ϕ ∈ [− 30,30] • in accordance with the nautical device presented in [13]. The s 4 configuration was chosen because, in contrast to s 1 systems, the symmetry allows for integration of two sensors in the same plane (i.e., on a single printed circuit board (PCB)), each one in a different sensing lobe. The second sensor can be used both for redundancy reasons and for the need to compensate external magnetic stray fields by evaluating a differential signal.
The chosen system parameters are the following: a cubical magnet with side length 5 mm, an airgap of g = 2 mm, a magnet distance to the center of tilt of d CoT = 6.83 mm, and a sensor displacement of x s = ±5.52 mm. These parameters are selected as a result of optimization, as outlined in Section 3.3. This optimization includes potential large system fabrication tolerances (δx m = ±0.5 mm, δȳ m = ±0.5 mm, δz m = ±0.5 mm, δx s = ±0.25 mm, δy s = ±0.25 mm, and δg = ±0.5 mm) to ensure reliability and stability of the mechanical system.
The experimental setup is outlined in Figure 8. A custom three-axis joystick is coupled to a robot arm to realize precise mechanical states. The chosen four-axis pick-and-place robot is commonly used to examine linear position and angle sensors. The robot arm is parallel to the z-axis and can move linearly to any point in x-y-z-space with a precision of 30 µm. A rotation of the arm around the z-axis is possible with a precision of 0.02 • . The joystick lever is connected to the robot arm with the help of homokinetic coupling with an integrated length compensation element. This configuration overcomes several difficulties that arise from coupling the linear robot motion to the spherical joystick motion, including a necessary robot arm length compensation, a nonlinear connection between joystick rotation and robot rotation as well as a nonuniform transmission of torque which is not ideal to account for clearances in the setup. Further details on this kind of couplings are reported in [48]. The expected positioning error in this implementation is 0.01 • in tilt and 0.02 • in rotation when clearances are neglected.
The joystick itself is realized through a center ball (half sphere) which is connected to the rod and pressed into a spherical cavity by two springs and a ball-bearing in order to minimize further clearances in the setup. Figure 8a shows a sketch of the setup, whereas a picture of the actual mechanical realization is displayed in Figure 8b. The magnet is fixed to the center ball, while two 3D sensors are mounted on a PCB and integrated into the system along the positive and negative directions of the x-axis; see Figure 8c,d. To account for the low state separation, a high-precision 3D Hall sensor with a resolution of 8 µT is chosen. The sensor is configured in such way that it provides raw values for the three components of the magnetic flux density.
In the experiment 16,200 mechanical positions are measured in an angle grid with step sizes ∆ψ = 8 • , ∆θ = 1 • and ∆ϕ = 2.5 • . A comparison between experimental data and theoretical predictions (simulated as discussed in Section 2.4) is shown in Figure 9 for the two sensors at x s = 5.5 mm and x s = −5.5 mm.
The theoretical values show a mean deviation from experimental measurements by 0.1 mT for both the sensors in x s = −5.5 mm and x s = 5.5 mm. The high consistency between theory and experiment is achieved by fitting the theoretical predictions onto the experimental data by variation of 27 tolerances that include sensor position, orientation, gains and offsets, magnet position, orientation, dimensions, magnetization, and experimental angle offsets. All computed tolerances lie within reasonable ranges. The resulting mean angles errors are e θ = 0.067 • , e ϕ = 0.126 • and e ψ = 0.491 • and at the 99 percentile, the maximum errors are e θ,max = 0.481 • , e ϕ,max = 0.658 • and e ψ,max = 2.278 • for single sensor evaluation. Such errors are at the same level with literature values of 3-DoF and 6-DoF motion tracking [8] where, however, multiple sensors are used for read-out. Figure 9. Comparison between experimental measurements (red dots) and analytical calculation (black lines) in an s 4 -type system with two sensors located at x s = −5.5 mm (a) and x s = 5.5 mm (b), respectively.

Conclusions
In this work, the difficulties related to magnet position and orientation (MPO) system design are discussed. A new method for the computation of MPO system layouts is proposed, aimed at maximizing the state separation through a global optimization procedure that is enabled by computationally fast analytical solutions of permanent magnet problems. A comparison to sophisticated topology optimization shows that the proposed computationally inexpensive ansatz can achieve excellent results.
The formalism is applied to study the three-axis-joystick problem. It is shown for the first time that continuous three-axis motion tracking can be realized with only a single magnet and a single 3D magnetic field sensor when fulfilling a specific design criterion that relies on perpendicularity between rotation axis, magnet or sensor displacement, and magnetization direction. For realistic cubical magnets with 5 mm sides and 1000 mT remanence field, a full 360 • rotation and tilts up to 12 • can be realized at 2 mm airgap with state separation close to 1 mT/ • . Much larger tilt angles beyond 60 • can also be achieved at the expense of state separation. The proposed systems allow for three-axis-joystick motion tracking at unparalleled cost-efficiency. The computations are confirmed by an experimental study where the measured fields are consistent with the theoretical predictions, and resulting mean angle errors are below 0.5 • .
This study demonstrates the potential of the proposed formalism in designing novel MPO systems. Future work is dedicated to include system fabrication tolerances in the design process as well as efficient calibration, inversion strategies, and stray field compensation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Connectedness Requirement
The requirement for B space to be simply connected suppresses a potential periodicity of the sensor output that would render the one-to-one correspondence between mechanical states P and sensor outputs B space invalid. This is best demonstrated with an example. Consider a quadrupole magnet disc that is used in an end-of-shaft application [6,43]. The shaft angle ϕ ∈ [0, 360] • should be determined by a 2D magnetic field sensor located above the disc. The sensor outputs are then B = B 0 (cos(2ϕ), sin(2ϕ)). In this case, J T J = 4B 2 0 is computed so that Equation (3) always holds true. The sensor output is also always smooth; however, it is periodic with 180 • so that a one-to-one correspondence is not possible for all rotations. Such a periodicity, which makes the chosen implementation invalid, becomes immediately visible by checking the iso-surfaces of B space (as outlined in Section 3.1).

Appendix B. Field Shaping with Topology Optimization (Adjoint Method) and Shape Variation
Magnetic field shaping relates to optimization problems for designing magnet systems that are able to generate given target fields B target . Therefore, the error d(B, B target ) is minimized with respect to a suitable chosen deviation measure d(·, ·). Using the formalism introduced in Section 2.1, this can be expressed through the quality factor Q(s) := −d(B, B target ), which has to be maximized with respect to the system design parameters s ∈ S. In this appendix, two different fundamental approaches and height of the new magnet pair). The result of such an optimization is displayed in Figure A2c,d for two and three layers of magnets, respectively. x y z (b) Figure A1. Illustration of the dimensions and the initial degrees of freedom (DoFs) for (a) topology optimization (magnetic material separated in cubic cells with edge length 0.2 within the red target area) and (b) the shape variation method (length of magnets and distance in-between, marked with green and yellow arrows, respectively). Deviations of the optimization results with respect to the optimal solution are reported in Table A1. Table A1. Deviation of the optimum quality factor Q(s) defined in Equation (A3) for the shape variation with various number of magnet layers, corresponding to Figure A2b-d: values are reported in % compared to the solution of the adjoint method in Figure A2a. It can be seen that shape variation performs surprisingly well, reaching more than 99% of the possible quality in the presented optimization problem, despite the crude geometric approximation when compared to the optimal solution. This discrepancy is expected to become even smaller by increasing the distance from the magnet target region. For many applications, such a small gain might not justify the additional effort to realize the complex optimal magnet shape. While the presented problem is reminiscent of MPO systems, it is also of a very simple nature and it must be pointed out that no general conclusion should be drawn. A reasonable summary of the findings of this section is given in Table A2. Table A2. Overview of the comparison between optimization with topology optimization and shape variation.

Topology Optimization Shape Variation
number DoF large small optimization algorithm gradient based (local) function evaluation based (global) quality factor requires derivation requires fast evaluation optimum (magnet shape) very accurate, but possibly local inaccurate, depending on choice of DoF optimum (magnetic field) very accurate, but possibly local possibly very accurate, depending on DoF choice application case theoretical optimum solutions good practical solutions

Appendix C. Computation of the Rotation Matrix
Here, the rotation matrices used to track the motion of the magnet are given explicitly. To this aim, it is convenient to split the three-axis motion in two independent (non-commutative) rotation operations. First, the lever-and with it, the magnet-is rotated about the z-axis by the angle ϕ expressed via the rotation matrix R 1 (ϕ): Then, for the application of the tilt, the lever is rotated by an angle θ about an axis (cos(ψ), sin(ψ), 0) that passes through the center of tilt in the origin of the global coordinate system. Such an angle-axis rotation is mathematically formulated by means of the quaternion notation. The unitary quaternion u representing the rotation of an angle θ about the axis u i can be written as u = cos θ 2 , u 1 sin θ 2 , u 2 sin θ 2 , u 3 sin θ 2 . The components of the rotation axis (u 1 , u 2 , u 3 ) are calculated by means of the spherical coordinates transformation: u i = (sin(π/2) cos ψ, sin(π/2) sin ψ, cos(π/2)), where ψ is the azimuth angle and π 2 is the polar angle. Hence, the explicit expression of the unitary quaternion takes form: u = cos θ 2 , cos(ψ) sin θ 2 , sin(ψ) sin θ 2 , 0 and the corresponding rotation matrix is Finally, the full three-axis motion is described via R(α) = R 2 (θ, ψ)R 1 (ϕ). Such a matrix is used in Equations (7), (8), and (9) to compute the magnet position r m (α, s), orientation e m i (α, s) and the magnetic field B(α, s) of a given implementation s in the LCS. Eventually, R(α) is given in full by in principle invalid. In the following, accuracy of the analytical formulas is tested when comparing with magnets having a realistic material response. The simplest material law for permanent magnets is given by M(H) = M 0 + χ r H, which corresponds to a linearized hysteresis curve about an imprinted remanence magnetization M 0 . The remanent susceptibility χ r describes a linear material response to the field H.
Modern high-grade SmCo, NdFeB, or ferrite materials feature extremely large coercive fields [53] and small demagnetization slopes as low as χ r > 0.05, as found in the product portfolio of any permanent magnet manufacturer. As a consequence, Equation (A8) is very accurate and demagnetization effects are expected to be small. A method of moments code [54] is implemented to determine the demagnetization effects of a cubic magnet with M 0 e z . For a chosen discretization of the cube into 6859 (19 3 ) cells, the relative error of the magnetic field of the MOM computation easily undercuts 10 −4 in this study, as was determined by an additional FE simulation (ANSYS Maxwell).
The simulation results are displayed in Figure A3. The field is sampled along several representative lines (blue): perpendicular to the magnet surfaces (solid lines), perpendicular to the edges (dashed lines), and outwards from the corner (dotted line) in the first octant defined by the cube. The relative error along these lines is shown (in percent of the field amplitude) in Figure A3b for two different susceptibility values χ r = 0.1 (red) and χ r = 0.05 (yellow). The seven different lines are represented by multiple markers at each distance.
Without additional compensation, demagnetization effects result in an error up to a few percent of the field amplitude (circles). However, this error is mostly of a quantitative nature and can be analytically compensated simply by correcting the total dipole moment. The correction factors are M/M 0 ≈ 0.9679 for χ r = 0.1 and M/M 0 ≈ 0.9837 for χ r = 0.05. An approximation of the correction factors can also be computed analytically [55]. The relative error of the compensated analytical solution (triangles) decays quickly with the distance from the magnet as the cube field starts to approach the form of a dipole field. Variation problems that include the magnetization amplitude (see, e.g., the fitting in Section 3.4) automatically draw he high accuracy of the compensated analytical solution, but even the few percent errors associated with the uncompensated solution are typically acceptable for magnet system design.