Convex optimization of measurement allocation for magnetic tracking systems

Magnetic tracking is a popular technique that exploits static and low-frequency magnetic fields for positioning of quasi-stationary objects. One important system design aspect, which substantially influences the performance of the tracking system, is how to collect as much information as possible with a given number of measurements. In this work, we optimize the allocation of measurements given a large number of possible measurements of a generic magnetic tracking system that exploits time-division multiplexing. We exploit performance metrics based on the Fisher information matrix. In particular, the performance metrics measure worst-case or average performance in a measurement domain, i.e. the domain where the tracking is to be performed. An optimization problem with integer variables is formulated. By relaxing the constraint that the variables should be integer, a convex optimization problem is obtained. The two performance metrics are compared for several realistic measurement scenarios with planar transmitter constellations. The results show that the worst performance is obtained in the most distant parts of the measurement domain. Furthermore, measurement allocations optimized for worst-case performance require measurements in a larger area than measurement allocations optimized for average performance.


Introduction
Magnetic tracking systems are designed to estimate the position and/or orientation of a specially designed object by means of its interaction with static or lowfrequency magnetic fields. Given that the human body is transparent to magnetic fields at these frequencies, magnetic tracking systems are popular within the biomedical engineering community. For example, magnetic tracking has been used for eye tracking to diagnose Ménière's disease , positioning of wireless capsule endoscopes within the gastro-intestinal tract (Yang et al. 2009), real-time organ-positioning during radiotherapy of cancer tumors (Iustin et al. 2008), catheter tracking (Krueger et al. 2005; Biosense Webster 2011), monitoring of heart valve prostheses (Baldoni and Yellen 2007), tongue movement tracking (Gilbert et al. 2010;Wang et al. 2013), tracking of lung segment movements (Leira et al. 2012), and positioning of bone-embedded implants (Sherman et al. 2007). Examples of non-medical applications of magnetic tracking include head tracking for helmet-mounted sights in military aircraft (Raab et al. 1979), underground drilling guidance (Ripka et al. 2012), augmented and virtual reality (Liu et al. 2004), and tracking of the ball during an American football game (Arumugam et al. 2011).
In general, the performance of a measurement system is improved if the number of measurements increases because more information is collected and noise tends to be averaged. Nevertheless, this comes at the cost of more expensive hardware, lengthier measurement time and longer post-processing of the collected data. Therefore, a key issue in the design of measurement systems is how to maximize the information gained per measurement.
This question is fundamental to the theory on the (optimal) design of experiments that has been extensively applied to geo-spatial sciences for problems in agriculture, geology, meteorology etc. The reader is referred to Walter and Pronzato (1997), Uciński (2005), Pukelsheim (2006), Atkinson et al. (2007), and Pronzato and Pázman (2013) for an introduction to the subject. Joshi and Boyd (2009) studied sensor selection by means of convex optimization without a specific application in mind. Examples of electromagnetic applications include optimization of measurement setups for antenna measurements in the near-field (Nordebo and Gustafsson 2006), tracking of human tongue movements (Wang et al. 2013), estimation of current densities in magnetic resonance imaging magnets (Begot et al. 2002), and reconstruction of AC electric currents flowing in massive parallel conductors (Di Rienzo and Zhang 2010).
Within the magnetic tracking community, the impact of the number of measurements has been studied by Schlageter et al. (2001) and Plotkin and Paperno (2003). Schlageter et al. (2001) found that the accuracy of their magnetic tracking system was improved when the number of transmitters, and thus the number of measurements, was doubled. Plotkin and Paperno (2003) found that using more transmitters reduces the number of local minima present in the inverse problem. In contrast, how to obtain as much information as possible from a given number of measurements has received little attention. A rare example of such a study is the work by Shafrir et al. (2010) in which the positions of a fixed number of transmitters are optimized using a two-step evolutionary algorithm. However, their approach is devoted to a specific estimator and it requires that a positioning algorithm is executed a large number of times to build statistics.
In this work, we consider magnetic tracking systems that exploit time-division multiplexing and study how to allocate measurement efforts in an optimal way given a large number of possible measurements. We exploit the theory on the optimal design of experiments and formulate performance metrics based on the Fisher information matrix. The optimization of measurement allocation yields an optimization problem with integer variables. We approximate the integer variables by real variables, which gives us a convex optimization problem. In contrast to the method presented by Shafrir et al. (2010), the proposed method is valid for all unbiased estimators and it does not require a massive amount of computation. Furthermore, the convex nature of the proposed method is very attractive because it removes two difficulties commonly encountered in design of experiments optimization problems, namely, high dimensionality and presence of several local minima that are not globally optimal. Also, the convexity of the method proposed in this work makes it feasible to treat large scale problems.
In this work, we optimize for a measurement domain of arbitrary shape by formulating two cost functions that improve (i) the worst-case performance (minimax approach) and (ii) the expected performance for an assumed prior distribution of the position and orientation of the object we wish to track (average approach). The two approaches are compared for several test cases. Furthermore, we investigate optimal measurement allocation for a realistic measurement scenario. Finally, we study the impact of restrictions on the transmitter positions, which are commonly encountered in practice.
The paper is organized as follows. The modeling of a generic magnetic tracking system is presented in Sect. 2. Section 3 presents performance metrics and the proposed solution methods. The results are then presented in Sect. 4 and discussed in Sect. 5. Finally, the work is concluded in Sect. 6.

Modeling of the measurement system
Consider a quasi-magnetostatic tracking system operating at a single frequency. The tracking system consists of (i) one receiving coil with unknown positionx r ;ỹ r ;z r ð Þ and unknown orientationm r ¼ ðm r x ; m r y ; m r z Þ, and (ii) N t identical transmitting coils (also referred to as transmitters) with known positionsx k t ;ỹ k t ;z k t ð Þ and known orientationsm t k . Here and in the following, a vector a ¼ aâ is represented by the magnitude a and the unit vectorâ. The tracking system exploits time-division multiplexing to separate the signals from the different transmitters, i.e., the transmitters are operated in sequence such that only one transmitter is transmitting at any given time instant.
The aim of the tracking system is to estimate the position and orientation of the receiving coil, i.e. to estimate p ¼x r ;ỹ r ;z r ;m r x ;m r y ;m r We assume that the physical properties of the sensor are known, which is normally the case in practice. Thus,m r is known and we can use the unit vectorm r ¼m r =m r directly instead ofm r , without loss of generality.
To obtain entries with identical units in the vector that we wish to estimate, the spatial coordinates are normalized with the distance d, which yields Thus, the p ¼ 5 degrees of freedom that are to be estimated are described by Let R k ¼ r r À r t k denote the distance vector of length R k from the transmitting coil k to the receiving coil. By modeling the transmitting and receiving coils as magnetic dipoles and exploiting Faraday's law, the scaled induced voltage in the receiving coil generated by transmitting coil k is given by Jackson (1998) where x is the angular frequency, l 0 is the permeability of free space, and V 0 is a reference voltage that renders V k unit-less and thereby independent of the unit of measurement. The parameter a k is assumed to be known and it describes the diameter, number of turns, and the excitation current for each transmitting coil k. We use xa k =V 0 ¼ xa=V 0 ¼ 4:33 Â 10 6 Am/Vs for all k throughout this work, which implies that all transmitting coils are identical. The gradient of the scaled induced voltage in the receiver generated by transmitting coil k with respect to the position of the receiver r r is given by and the gradient of the scaled voltage with respect to the magnetic dipole moment of the receiver m r is given by Notice that the gradient with respect to the position in (6) scales as R À4 k whereas the gradient with respect to m r in (7) scales as R À3 k . Consider a measurement scenario where the receiver can be assumed to be stationary in both position and orientation during a time DT. The time to perform one measurement is Dt, which corresponds to the time required to record and process the signal generated by one of the transmitters. In this article, we are focused on the positioning of a mechanical object that is quasi-stationary on time scales that are many order of magnitudes larger than the time scale associated with the electrical system that performs the measurement. Therefore, we make the assumption that DT is many orders of magnitude larger than Dt, which reflects many real-life situations. Thus, the maximum number of measurements that can be collected during the time DT with stationary conditions is limited by the large number N meas ¼ DT=Dt, which follows from that the measurement system is based on time-division multiplexing. For convenience, we assume that N meas is an integer in the following given the nature of an actual measurement system, i.e. measurements are collected and processed as single units by standard off-the-shelf measurement instruments.

Optimization problem
In this work, we seek to improve the performance of the tracking system by allocating the N meas measurements among N t candidate transmitters in an optimal way. Let w k 2 N be the number of measurements performed with transmitter k. Clearly, it is advantageous to perform as many measurements as possible during the stationary time-interval DT and, thus, we have P k w k ¼ N meas . Let w ¼ w 1 ; w 2 ; . . .; w N t ½ T . Thus, we want to minimize w Jðp; wÞ subject to w k 2 0; 1; 2; . . .; N meas where J is a cost function quantifying the system's performance and X p is the measurement domain for which we want to optimize the tracking system. This is a combinatorial optimization problem and an exhaustive search requires ðN t Þ N meas cost function evaluations, which is prohibitive from a computational perspective. This is Convex optimization of measurement allocation for magnetic... 853 particularly true for more complicated measurement scenarios that may require parameter studies that involve the solution of many optimization problems. Now, let k k ¼ w k =N meas denote the fraction of the total number of measurements that are performed with transmitter k. Thus, k k 2 0; 1=N meas ; 2=N meas ; . . .; 1 f g . We use the approximation k k 2 0; 1 ½ because N meas is large as discussed in Sect. 2 above. By using the notation K ¼ k 1 ; k 2 ; . . .; k N t ½ T , we obtain the relaxed optimization problem minimize K Jðp; KÞ subject to k k 2 0; 1 which is a good approximation to the problem in (8). The feasible domain dictated by the constraints is convex. Thus, if the cost function J is convex with respect to K, the entire optimization problem is convex and can be readily solved.
In the following subsections, we introduce a performance metric, present the cost functions that are used, and present the method to solve the corresponding optimization problems.

Performance metric
Let V meas k ðp 0 Þ denote the measured signal generated by transmitting coil k for an arbitrary receiver position and orientation p 0 2 R 3 Â S 3 in the parameter space, where R 3 is the position in three dimensional space and S 3 is all possible directions of orientation on the unit sphere in R 3 . Noise that is caused by, for example, thermal noise in amplifiers can degrade the performance of the positioning system. Therefore, we model the measured signal as the true signal V k ðp 0 Þ corrupted with additive Gaussian noise as where the noise terms n k $ N ð0; r 2 Þ are independent and identically distributed and N ðl; r 2 Þ denotes the Gaussian distribution with mean l and variance r 2 . Below, we denote the gradient of V k ðpÞ with respect to the parameters in p at the point p 0 by r p V k ðp 0 Þ. A metric for the performance of the parameter estimation is provided by the Fisher information matrix M 2 R pÂp (Kay 1993) given by because the Cramér-Rao inequality (Walter and Pronzato 1997) yields a lower bound for the covariance of the estimatep for all unbiased estimators.
Here, A # B signifies that the matrix A À B is positive semi-definite. Furthermore, the bound can be attained, for example, asymptotically by the maximum-likelihood estimator. Therefore, the performance of the measurement system is expected to improve by maximizing M (in some sense). However, to find an optimal Fisher information matrix M Ã that fulfills M Ã # M; 8M 6 ¼ M Ã is, in general, not possible (Uciński 2005) and a real-valued function JðMÞ is often optimized instead.
Here, we use that yields a so-called D-optimal (Determinant-optimal) solution. If the model V k ðpÞ is linear in p, the D-optimal solution minimizes the volume of the lower bound for the confidence ellipsoid described by M À1 in (12). The volume of the b-confidence ellipsoid is given by Pronzato and Pázman (2013) where F v 2 p is the cumulative distribution function for the v 2 -distribution with p degrees of freedom and C denotes the Gamma function. The geometric mean of the lengths of the confidence ellipsoid's semi-axes, which we refer to as the mean confidence-radius, is given by Joshi and Boyd (2009) An attractive feature of the D-optimality criterion is that it is invariant to scaling of the parameters in p (Uciński 2005). By using the cost function from (11) and (13) as well as X p ¼ p 0 in (9), we obtain the relaxed local design problem which is a convex optimization problem as shown by Boyd and Vandenberghe (2004, Section 7.5).

Local and non-local designs
The Cramér-Rao inequality in (12) yields a lower bound for the covariance of the estimated parameters. Given that the covariance is a measure of the linear relationship between the estimated parameters, it does not capture their true relationship when the functions V k ðpÞ are non-linear in p. In addition, M À1 is a function of p 0 because of this non-linearity. This is the reason why an optimal experiment design based on (11) and (12) for V k ðpÞ non-linear in p is referred to as a local design (Walter and Pronzato 1997). The region of validity of a local design depends on the size of the region where the linearization V k ðpÞ ffi V k ðp 0 Þ þ r p V k ðp 0 Þðp À p 0 Þ is a good approximation to the true non-linear V k ðpÞ.
In contrast to local designs, it is often desired to optimize the performance of the measurement system not just for one point in the parameter space R p but rather for a measurement domain X p R p . In this work, we optimize the measurement performance in X p for (i) average optimality and (ii) minimax optimality. To this aim, we exploit a discrete set of linearization points X lin ¼ fp i g N lin i¼1 X p that constitutes a sufficiently dense discretization of X p .
Average optimality In this case, we assign a prior probability distribution p p ðpÞ for the parameters that are to be estimated. We then find the so-called ELD-optimal (Expectation of Log Determinant-optimal) experiment design (Walter and Pronzato 1997) by minimizing the cost function where E denotes the expectation with respect to p p ðpÞ. In our case, we assume a uniform prior probability density and so the cost function can be written as We evaluate this integral by quadrature at the linearization points p i 2 X lin with weights q i as with the quadrature scheme described in Quadrature. This quadrature scheme features weights q i that are positive, which is important to preserve the convexity of the optimization problem. In addition, it preserves symmetries with respect to the m r x m r zand m r y m r z -planes. Furthermore, the value of det M is unaffected ifm r is multiplied by -1, which is exploited by performing the quadrature for the half sphere m r z ! 0 only. Notice that minimization of the continuous cost function in (17) and its discretized version in (20) are equivalent within the accuracy of the quadrature scheme. Also, notice that there are other more efficient quadrature schemes for evaluating the expectation in (17) than the one presented in Quadrature. An example of such a scheme is the one given by Gotwalt et al. (2009) and Gotwalt (2010) that, however, includes negative weights in situations where more than 7 parameters are to be estimated.
Thus, we obtain the relaxed average optimality problem that is also a convex problem because the cost function is a sum of convex functions with positive weights.
Minimax optimality In many applications, it is often desired to guarantee a certain accuracy of the measurements. In these cases, the worst-case performance is optimized instead of the average performance. This leads to a so-called minimax problem, where we seek to minimize the MMLD (MiniMax Log of Determinant) cost function The computation of J MMLD ðKÞ involves solving a separate optimization problem defined by the right hand side of (22). We solve this optimization problem by computing J MMLD ðKÞ at N lin linearization points p i 2 X lin and taking the maximum value, i.e.
J MMLD ðKÞ % max which yields the relaxed minimax optimality problem This is a convex problem because the pointwise maximum of convex functions is convex (Boyd and Vandenberghe 2004).
Convex optimization of measurement allocation for magnetic... 857

Solution method
In this work, we seek to allocate a limited number of measurements given a large number of possible candidate measurements in an optimal way. Apart from the information on the measurement allocation, the solutions to our optimization problems also inform us of the number of transmitters to use and their positions. A lower bound on the number of transmitters is given by the number of parameters p that are to be estimated. The number of measurements to use is limited by the constraint P N t k¼1 k k ¼ 1 in (21) and (24) and this constraint shows strong similarity with penalty terms encountered in compressed sensing and related problems (Bruckstein et al. 2009). Such a penalty term typically involves the L 1 -norm of the solution vector and it is added with a weight to the cost function that should be minimized. In the context of compressed sensing and related problems, the penalty term favors a sparse solution with only a few non-zero entries, should such a solution be consistent with the rest of the problem statement. Here, we find that the optimized measurement allocation vectors K Ã computed from (21) and (24) feature only a few non-zero entries in comparison to the number of transmitter candidates, which is confirmed by the results presented in this article.

Thresholding and clustering of weights
The weights k k that are obtained in the solutions of the convex problems (21) and (24) above can, for the examples we have studied in this work, be grouped as follows: (i) a handful of the weights are large ( [ 10 À3 ); (ii) many are zero; (iii) several are nearly zero (\10 À9 ). That the weights of the last group are not zero is due to the finite precision arithmetic and termination criteria tolerances of the exploited numerical solver. In addition, the weights of the last group are several orders of magnitude smaller than the weights of the first group. Therefore, we use the threshold k th ¼ 10 À6 and set all weights k k \k th equal to zero. We refer to weights k k ! k th as non-zero.
Furthermore, the finite resolution of a Cartesian grid of transmitter candidates may cause several neighboring transmitters to obtain non-zero weights k k . We replace such a cluster X cl of non-zero weights k k ; k 2 X cl with only one weight k k cl placed at r k cl according to if all the non-zero weights k k in the cluster are vertices of the same cell in the Cartesian grid. If this is not the case, e.g. there are five non-zero k k in the cluster or the cluster consists of three non-zero k k on a straight line, we do not perform the clustering. Instead, the problem should be solved again for a denser grid of transmitter candidates.

Evaluation of derivatives
To solve the optimization problems (21) and (24), the Fisher information matrix for a given receiver position and orientation p i 2 X lin must be computed. Thus, the derivatives with respect to the two degrees of freedom given by are needed. However, the gradient r ðm r Þ V in (7) includes the derivatives with respect to the Cartesian components of the receiver's magnetic dipole moment m r . If all three of these components are included in the Fisher information matrix, the constraint in (26) makes the Fisher information matrix rank-deficient. We therefore introduce a local Cartesian coordinate system ðû i ;v i ;ŵ i Þ withŵ i ¼m r i to express r ðm r i Þ V. The u i -and v i -components of r ðm r i Þ V are then used in the computation of the Fisher information matrix. (The w i -component of r ðm r i Þ V is always zero because of the constraint in (26).) The cost functions that are exploited in this work are unaffected by a rotation ofû i andv i aroundŵ i because determinants are invariant to rotations.

Solver
The relaxed average optimality problem in (21) and the relaxed minimax optimality problem in (24) are solved directly with the routine SNOPT (Gill et al. 2005) provided in the TOMLAB (Tomlab Optimization AB 2012) package of optimization algorithms. The SNOPT-routine is an implementation of the sequential quadratic programming algorithm. All gradients that are needed are computed analytically by SNOPT.

Results
Planar transmitter constellations have become increasingly popular, see for example (Iustin et al. 2008;Plotkin et al. 2010). In this work, we therefore consider only planar constellations of transmitters. More specifically, we consider constellations where all transmitters lie in the plane z ¼ 0 with dipole moments oriented along the z-axis, i.e. z t k ¼ 0 andm t k ¼ẑ for all k. (It should be noted that the proposed method can handle any geometry of the transmitter constellation. Furthermore, transmitters with different orientations can also be considered with the method, should this be desired.) In particular, we consider two types of planar transmitter constellations based on (i) a Cartesian grid of transmitter candidates and (ii) a polar grid of transmitter candidates. These transmitter constellations are referred to as Cartesian arrays and polar arrays, respectively, in the following.
Examples of the two transmitter array types are shown in Fig. 1. The transmitters in a Cartesian array are placed on a Cartesian grid with jx t k j x max , jy t k j y max and an inter-transmitter distance h in both the x-and y-directions. The transmitters in a polar array are placed on a polar grid with N arms transmitters per circle, 0 ðx t k Þ 2 þ ðy t k Þ 2 r 2 max and a radial distance h between neighboring circles. On each circle, the transmitters are placed at the polar angles w l ¼ l 2p N arms where l ¼ 1; . . .; N arms . Below, we compare the average and minimax cost functions in Sect. 4.1. Then, we study a realistic measurement scenario in Sect. 4.2. Finally, we investigate the impact of restrictions on the permissible size and position of the transmitter array (Sect. 4.3). Note that we use r 2 ¼ 1 in the following tests without loss of generality. All the results in this section are presented in terms of the continuous weights K, which we find useful and informative in an engineering setting. The continuous weight k k can directly be interpreted as the fraction of measurements that are to be collected based on transmitter candidate k, which is useful since k k is not explicitly dependent on N meas . In other words, the solution k k describes a variety of measurement systems that feature different values of N meas , which may involve widely different hardware implementations. In addition, the continuous weights may be used as a good starting point for the combinatorial optimization problem (8), which may be approached in a number of different ways depending on the application at hand and computational resources available. For sufficiently large values of N meas , the weights k k can be rounded to an integer multiple of 1=N meas without any significant change in the performance of the measurement system, i.e. the objective function in (8) is basically unaltered given the real-world measurement situation. Should N meas not be sufficiently large for the application at hand, the approach presented by Joshi and Boyd (2009) can be used to pursue the solution of the combinatorial optimization problem in (8), where the relaxed solution may be used as a starting guess. However, we are focused on positioning of a mechanical system that is quasi-stationary on time scales that are many orders of magnitude larger than the time scale associated with the electrical system that performs the measurement, which implies that N meas is indeed very large for any practical purposes.

Cost function comparison
In order to illustrate the differences between average and minimax optimality, we consider a simple problem with cylindrical symmetry. The measurement domain is given by where N lin ¼ 1000 linearization points for z r 2 ½0:1; 1 are exploited by the quadrature scheme described in Quadrature. We consider a polar array defined by r max ¼ 1:1, h ¼ 0:0025 and N arms ¼ 3 with N t ¼ 1321 candidate transmitters. Next, we constrain the transmitters on each circle of constant radius in the array to have equal weights, which is motivated by the symmetry of the problem. (It should be noted that the symmetry is broken by the transmitter array. However, we obtain identical results for N arms ¼ 3; 4; . . .; 7. Furthermore, circles centered at the origin are formed by the transmitters with non-zero weights obtained by solving the problem with a Cartesian array of transmitter candidates, where no additional constraints on the weights are incorporated.) Figure 2 shows the non-zero weights of the solution to the relaxed average optimality problem (21) and the relaxed minimax optimality problem (24). For this measurement scenario, the clustering procedure described in Sect. 3.2.1 is modified such that all radially adjacent weights k k ! k th are clustered into one single weight k k cl placed at r k cl according to (25). The corresponding radii and total weights of the circles with non-zero weights are given in Table 1. The non-zero weights for minimax optimality are fewer and constitute a larger constellation than the non-zero weights for average optimality. Furthermore, optimizing for minimax optimality yields the same result as optimizing only for the sensor position that is furthest away from the transmitter plane, i.e. z r ¼ 1, cf. (Talcoth and Rylander 2013). Figure 3 shows the pointwise cost J D ðp i ; K Ã Þ as a function of z r by dashed and solid curves for the solutions to the average and minimax optimality problems,  respectively. As can be seen in Fig. 3, the performance of the system degrades as the distance between the sensor and the transmitter plane increases. This can also be seen by combining (6), (7) and (11) where the distances R to the sensor are assumed to scale in the same way for all contributing transmitters. This also explains why the optimum of the minimax optimality problem is identical to the optimum for pointwise D-optimality at z r ¼ 1 (Talcoth and Rylander 2013) and the larger constellation size as compared to the clustered non-zero weights that correspond to average optimality. As can be seen from the curves in Fig. 3, optimizing for minimax optimality gives a slight improvement in worst-case performance as compared to optimizing for average optimality because qðK Minimax Þ=qðK Average Þ % 0:81 at z r ¼ 1. Here, q is the mean confidence-radius from (15). Further, K Minimax and K Average denote the measurement allocations optimized for minimax and average optimality, respectively, and their clustered non-zero weights are shown in Fig. 2. However, the improvement in worst-case performance comes at the expense of a large  degradation in performance close to the transmitter plane, e.g. qðK Minimax Þ=qðK Average Þ % 125 at z r ¼ 0:1. We also examine the impact of linearization point density by varying the number of linearization points in the measurement domain. For average optimality, at least 30 linearization points are needed to obtain the same constellation of clustered nonzero weights as described above. In contrast, only one linearization point at z r ¼ 1 is needed for minimax optimality because the worst performance is governed by the point furthest away from the transmitter plane.

A realistic measurement scenario
We investigate a realistic measurement scenario and quantify the potential for improvement of measurement allocation optimization as compared to an ad-hoc measurement allocation procedure.
We solve the relaxed average optimality problem in (21) and the relaxed minimax optimality problem in (24) with a transmitter candidate array of Cartesian type defined by x max ¼ 1:44, y max ¼ 1:44, and h ¼ 0:09 with N t ¼ 1089 candidate transmitters. Thresholding and clustering is applied as described in Sect. 3.2.1. Furthermore, we introduce an ad-hoc measurement allocation procedure consisting of a Cartesian array defined by x max ¼ 1, y max ¼ 1, and h ¼ 0:5 with 25 equally weighted transmitters, i.e. k k ¼ 1=25 for all k. This ad-hoc measurement allocation constitutes a natural choice for collecting measurements, should an optimization algorithm for measurement allocation not be available. (c) Minimax Fig. 4 Ad-hoc measurement allocation (left) and measurement allocations optimized for average optimality (middle) and minimax optimality (right). The clustered non-zero weights are represented by circular markers whose size is proportional to the weight k kcl . Transmitter candidate array boundaries are indicated with dashed lines Convex optimization of measurement allocation for magnetic... 863 Figure 4 shows the ad-hoc measurement allocation as well as the measurement allocations optimized for average and minimax optimality. The optimized measurement allocations are symmetric with respect to the x-and y-axes. Notice that these symmetries are not imposed during the solution of the optimization problems but are due to the symmetries present in the optimization problems. Also notice that all three measurement allocations require 25 transmitters. Similar to the results in Sect. 4.1, minimax optimality requires measurements over a larger area than average optimality. This is because of the scaling with respect to the distance R for the derivatives (6) and (7) and the cost function (28). Figure 5 shows J D ðMðp i Þ; K Ã Þ for all linearization points p i 2 X lin as a function of the linearization point index. Notice that these indices have been sorted individually for each case in non-decreasing order of the cost. All curves in Fig. 5 show four different levels that correspond to the different z r -values of the linearization points. Larger cost and, thus, worse performance is obtained for the most distant linearization points.
The cost function values for average optimality J ELD and minimax optimality J MMLD are given in Table 2 for the different measurement allocations. The best performance in terms of average optimality is shown by the measurement allocation optimized for average optimality, as expected. The mean confidence-radius of the minimax-optimal measurement allocation is 15% larger than the mean confidenceradius of the average-optimal measurement allocation. Similarly, the mean confidence-radius of the ad-hoc measurement allocation is 73% larger than the mean confidence-radius of the average-optimal measurement allocation. For minimax optimality, the increase in mean confidence-radius is 23% for the ad- Ad−hoc Minimax Average Fig. 5 Pointwise cost J D ðMðp i Þ; K Ã Þ as a function of linearization point index i for the ad-hoc measurement allocation as well as for the measurement allocations optimized for average and minimax optimality. Note that the indices i are sorted individually for each case to yield non-decreasing curves hoc measurement allocation and 35% for the average-optimal measurement allocation as compared to the minimax-optimal measurement allocation. Thus, we have shown an example where measurement allocation optimization provides substantial improvement of a measurement system as compared to an ad-hoc measurement allocation procedure. In this example, the improvement is especially large when average optimality is considered. Next, we consider the possible effect of rounding the elements in the optimized measurement allocation vector K Ã to multiples of 1=N meas , where we use a simplified analysis. Given an optimized measurement allocation vector K Ã , we consider the corresponding perturbed vectorK ¼ nK Ã . Here, all weights k k are scaled by the multiplicative factor n ¼ 1 þ dn, where dn is small in comparison to unity. For a perturbation dn, the relative perturbation in the mean confidence-radius (15) isqðbÞ=qðbÞ ¼ ð1 þ dnÞ À1=2 % 1 À dn=2. If k k N meas [ 20 for all non-zero weights k k , a pessimistic estimate of the relative change in the mean confidenceradius could be coarsely approximated by rounding all weights downwards to an integer multiple of 1=N meas . If we assume that the rounding (in the worst-case scenario) would correspond to roughly dn ¼ À0:05, we would have a relative perturbation in the mean confidence-radius ofqðbÞ=qðbÞ ¼ 1:026, i.e. a degradation of about 2.5%. This is a rather small degradation in the performance of the measurement system in relation to the improvements achieved by the relaxed solution, when the relaxed solution is compared to the ad-hoc measurement allocation. Should the combinatorial problem be solved, it is rather likely that the degradation in mean confidence-radius is much smaller than 2.5% for such a situation. Given the vast difference in time-scales of the quasi-static object and the measurement of the electrical system, we find that such improvements are in many cases of minor importance but could be pursued by, e.g., the technique presented by Joshi and Boyd (2009).

Impact of restrictions on transmitter candidate array size and position
In some measurement situations, it may be impossible to perform measurements underneath the sensor, i.e. the receiver is located such that its orthogonal projection onto the transmitter plane is located outside the region occupied by transmitter candidates. Also, only a part of the transmitter plane may be available for measurements. We study this scenario by considering the measurement domain Convex optimization of measurement allocation for magnetic... 865 where d corresponds to the distance to the point r r center ¼ ð0; 0; 1Þ above the center of the transmitter candidate array and ½v y ; v z ¼ ½3; 1= ffiffiffiffiffi 10 p . The relaxed minimax optimality problem in (24) is solved for different values of the distance d with N lin ¼ 605 linearization points on half of the unit sphere as described in Quadrature and a Cartesian transmitter candidate array defined by x max ¼ 1:2, y max ¼ 1:2, and h ¼ 0:04 with N t ¼ 3721 candidate transmitters. The thresholding and clustering procedure from Sect. 3.2.1 is exploited. The receiver is above the edge of transmitter candidate array for d ¼ d edge ¼ y max =v y % 1:26. Furthermore, the receiver is above the transmitter candidate array for 0 d\d edge and outside the transmitter candidate array for d [ d edge . Figure 6 shows the cost function and some optimized measurement allocations for different values of d. For small d, the receiver is above the transmitter candidate array and close to r r center . For these receiver positions, non-zero weights are found in all parts of the transmitter candidate array without any effect of its limited size and the performance of the measurement system is almost constant. For increased values of d, the receiver is found further away from r r center either above or outside the transmitter candidate array. For receiver positions in this region, the limited size of the transmitter candidate array strongly influences the measurement allocation and non-zeros weights are primarily obtained for positive y t -coordinates. Thus, weights with a shorter distance to the receiver are preferred to weights with a longer distance to the receiver. The measurement system performance decreases moderately with Examples of optimized measurement allocations are shown as inlaid plots with clustered non-zero weights (circles) and transmitter candidate array boundaries (dashed rectangle). The size of the circular markers is proportional to the corresponding weights k kcl increasing d above the transmitter candidate array and substantially outside the transmitter candidate array. In contrast, for receiver positions far outside the transmitter candidate array corresponding to large values of d, measurements are also allocated to weights with a negative y t -coordinate far away from the receiver. This suggests that more information can be gained by diversifying the measurements than what is lost by the increased distance to the receiver. For d [ 4, the cost function scales approximately as 49 log d instead of 36 log d as indicated by (28). The increase in the distance scaling factor is likely due to that measurements can only be allocated within a region that does not scale with d.

Discussion
Many of the difficulties and approximations in this work are related to the model being non-linear in the parameters that we wish to estimate. For example, the Fisher information matrix approximates the confidence volume with an ellipsoid. If the model is linear in the parameters, the confidence volume is indeed an ellipsoid. However, our model is non-linear in the parameters and, then, the confidence volume can take other shapes and does not even have to form a connected set. Therefore, the mean confidence-radius should only be considered as a qualitative metric because it is based on this approximation.
Local designs are based on the assumption that the parameter values that we wish to estimate are known. However, if the parameters are known, we do not need to estimate them. In this work, we have addressed this issue by optimizing for a range of possible sensor positions and orientations, where we have considered minimax and average optimality. An alternative approach is to exploit so-called sequential designs that updates the measurement procedure depending on already measured data. Plotkin and Paperno (2003) constructed a magnetic tracking system based on this idea (without using the design of experiments-terminology), where a subset of the transmitters in a 8 by 8 transmitter array is activated as a function of the most recently estimated sensor position.
To solve the average and minimax optimality problems, we perform quadrature in the measurement domain at a finite set of linearization points. As shown by the results in Sect. 4.1, few linearization points are needed when optimizing for minimax optimality if they are positioned at the most distant part of the measurement domain, i.e. where the worst performance is obtained due to the considerable distance scaling of the cost function. In contrast, more linearization points are needed for average optimality.
The optimization method for measurement allocation presented in this work can also be useful in other situations. For example, the measurement allocation result could be exploited as a starting guess for an optimization method that considers integer variables, a more elaborate physical model, or the impact of non-linearities and the choice of estimation procedure. Moreover, the convex nature of the method is advantageous. In particular, it permits large scale problems to be addressed and extensive parameter studies to be performed.
We have limited this study to planar transmitter constellations with known transmitter orientations due to their importance in practice. However, the proposed method can handle any type of transmitter constellation geometry. That is, transmitter candidates can take any position and orientation. Thus, transmitter constellations that occupy curved surfaces, several disjoint surfaces, volumes, etc., can be handled.
We have studied the situation where N meas is large. Should a situation where N meas is not large be encountered, the proposed method could be exploited as a first step. The obtained weights would then have to be ensured to be multiples of 1=N meas . This could for example be achieved by a local optimization procedure similar to the ones proposed by Joshi and Boyd (2009).
In this work, we have optimized measurement allocations to yield large changes in the measured signals for a change in the parameters that are to be estimated. We have not considered in full the characteristics of the estimation problem that is obtained with the optimized measurement allocation during the optimization procedure; for example, if the parameters can be uniquely determined everywhere in the measurement domain and if there are local minima present in the estimation problem. This is related to the concepts of identifiability and estimability and the reader is referred to Pronzato and Pázman (2013) for further information.

Conclusion
Magnetic tracking is a popular technique that exploits static and low-frequency magnetic fields for positioning of quasi-stationary objects. In this work, we have proposed a method for optimizing the allocation of measurements given a large number of candidate transmitters of a generic magnetic tracking system that exploits time-division multiplexing. The sensor and the transmitters are modeled as magnetic dipoles in free space. Performance metrics based on the Fisher information matrix are exploited to quantify the worst-case performance (minimax optimality) and the expected performance with respect to a prior distribution of the sensor's position and orientation (average optimality). Optimization problems with integer variables are formulated. By means of a convex relaxation, the integer variables are approximated with real variables and convex optimization problems are obtained. The proposed method is valid for all unbiased estimators and it avoids two commonly encountered problems, namely, high dimensionality and the presence of local minima that are not globally optimal.
The two performance metrics are compared for several realistic measurement scenarios where planar transmitter constellations are considered. Given the strong distance dependence of the measured signal, the worst-case performance is obtained in the most distant regions of the measurement domain. Consequently, measurement allocations optimized for minimax optimality requires measurements over a larger area than measurement allocations optimized for average optimality.
The optimized measurement allocations that are the result of solving the convex optimization problems can be used directly or as a starting guess for the solution of more detailed optimization problems.
In conclusion, the proposed method works well for optimization of measurement allocation for magnetic tracking systems that exploit time-division multiplexing and it provides useful and informative designs of such experiments.
Acknowledgements Oskar Talcoth has been supported by the Swedish research council VINNOVA via a project within the VINN Excellence center CHASE, and by the Swedish National Graduate School in Scientific Computing. The computations were performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Quadrature
Standard trapezoidal quadrature is exploited in the three spatial dimensions ðx r ; y r ; z r Þ whereas the quadrature on the surface of the three-dimensional unit sphere (corresponding tom r ) uses a subdivision of the curved surface in triangular elements as described below. The linearization points in five dimensions are obtained by combining each spatial quadrature point with all of the quadrature points on the unit sphere. Finally, the prior probability density is incorporated in the weights q i by normalizing their sum to one.
For quadrature on the entire or a part of the unit sphere, the surface of the considered part of the sphere is approximated by a mesh M with triangular elements and N nodes nodes. (The positions of the nodes on the sphere are symmetric with respect to the m r x m r z -and m r y m r z -planes. Thus, possible symmetries of the integrand with respect to these planes are preserved.) A piece-wise linear basis function v i ðm r Þ is associated with each node m r i in the mesh where v i ðm r j Þ ¼ The integrand f ðm r Þ is approximated by Integrating this over the unit sphere in R 3 and changing the order of integration and summation gives where T j is mesh element j and X i is the set of indices for the mesh elements that include node i. Note that the weights (the expression within curly brackets on the last line of (33)) are positive, which is also the case for trapezoidal quadrature. Thus, the weights q i exploited in the five-dimensional quadrature are positive.