Simultaneous calibration method for magnetic localization and actuation systems

Abstract

The invention relates to a method of simultaneously calibrating magnetic actuation and sensing systems for a workspace, wherein the actuation system comprises a plurality of magnetic actuators and the sensing system comprises a plurality of magnetic sensors, wherein all the measured data is fed into a calibration model, wherein the calibration model is based on a sensor measurement model and a magnetic actuation model, and wherein a solution of the model parameters is found via a numerical solver order to calibrate both the actuation and sensing systems at the same time.

Claims

1. A method of simultaneously calibrating magnetic actuation and sensing systems for a workspace, wherein the actuation system comprises a plurality of magnetic actuators and the sensing system comprises a plurality of magnetic sensors, the method further comprising the following steps: generating a plurality of arbitrary magnetic fields at the workspace by applying a plurality of electric currents to the magnetic actuators; measuring each magnetic field with the sensing system in the workspace and measuring each respective electric current with a current sensor; measuring a magnetic field of a random known sample magnet with the sensing system at an arbitrary position and orientation in the workspace; and feeding all the measured data into a calibration model, wherein the calibration model is based on a sensor measurement model and a magnetic actuation model, and finding a solution of the model parameters via a numerical solver in order to calibrate both the actuation and sensing systems at the same time.

2. The method in accordance with claim 1, wherein the actuation and sensing systems each comprise an actuation system level and a sensing system level, respectively, wherein the workspace is sandwiched in between the actuation system level and the sensing system level.

3. The method in accordance with claim 1, wherein the magnetic actuators comprise electromagnets.

4. The method in accordance with claim 1, wherein the magnetic sensors can be on of mono-axial, triaxial and multi-axial sensors.

5. The method in accordance with claim 1, wherein each of the plurality of generated magnetic fields comprises a non-uniform shape throughout the workspace.

6. The method in accordance with claim 5, wherein the non-uniform shape is generated from a solenoidal electromagnet.

7. The method in accordance with claim 1, wherein the calibration model comprises a parameter set, which includes sensor parameters and actuator parameters, wherein the sensor parameters comprise parameters for position, orientation, axis distortion, axes distortion and gain of each sensor, and wherein the actuator parameters comprise parameters for position, orientation and magnetic moment of each actuator.

8. The method in accordance with claim 1, wherein the sensor measurement model includes a measurement of each generated magnetic field according to $\begin{array}{cc}{\stackrel{\_}{b}}_{\mathrm{ik}}\left(p_{\mathrm{si}},v_{\mathrm{ik}}\right)\propto G_i{v}_{\mathrm{ik}}+\left[\begin{array}{ccc}{v}_{\mathrm{ik}}^T& H_{x,i}& v_{\mathrm{ik}}\\ v_{\mathrm{ik}}^T& H_{y,i}& v_{\mathrm{ik}}\\ v_{\mathrm{ik}}^T& H_{z,i}& v_{\mathrm{ik}}\end{array}\right],& \left(1\right)\end{array}$ wherein p.sub.s,i is an array of parameters for G.sub.i, H.sub.x,i, H.sub.y,i and H.sub.z,i, and wherein H.sub.x,i, H.sub.y,i and H.sub.z,l are symmetric 3×3 matrices which involve rotations and quadratic terms together, G.sub.i is a linear mapping which represents effects of scaling, rotation and shear, and v.sub.ik is the magnetic senor reading.

9. The method in accordance with claim 1, wherein recovered B-field values are expressed in a stacked matrix form as $\begin{array}{cc}{\stackrel{\_}{𝕍}}_s\left(p_s,v\right)\propto \left[\begin{array}{ccc}{\stackrel{\_}{b}}_{11}\left(p_{s,1},v_{11}\right)& \mathrm{.Math.}& {\stackrel{\_}{b}}_{1L}\left(p_{s,1},v_{1L}\right)\\ \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ {\stackrel{\_}{b}}_{N1}\left(p_{s,N},v_{N1}\right)& \mathrm{.Math.}& {\stackrel{\_}{b}}_{\mathrm{NL}}\left(p_{s,N},v_{\mathrm{NL}}\right)\end{array}\right],& \left(2\right)\end{array}$ and wherein p.sub.s and v are stacked vectors for all p.sub.s,i and v.sub.ik, respectively.

10. The method in accordance with claim 1, wherein the magnetic actuation model is based on a dipole model, where a dipole magnetic field on the i-th sensor from the j-th actuator at the k-th measurement is expressed as $\begin{array}{cc}{b}_{\mathrm{ij},k}\propto B_{\mathrm{ij}}{m}_j{I}_{\mathrm{jk}};B_{\mathrm{ij}}\propto \frac{{\mu }_0}{4\pi {.Math.r_{\mathrm{ij}}.Math.}^3}\left(3{\hat{r}}_{\mathrm{ij}}{\hat{r}}_{\mathrm{ij}}^T-𝕝\right),r_{\mathrm{ij}}\propto r_i-r_j& \left(3\right)\end{array}$ where b.sub.ij,k is a magnetic field vector, B.sub.ij is a 3×3 matrix which only depends on the displacement vector r.sub.ij, is the identity matrix and μ.sub.0 is the permeability of free space.

11. The method in accordance with claim 1, wherein the magnetic fields at multiple measurement points are expressed as $\begin{array}{cc}{b}_k\mathrm{𝔹𝕄}{i}_k,\mathrm{where}𝔹\propto \left[\begin{array}{ccc}{B}_{11}& \mathrm{.Math.}& B_{1M}\\ \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ B_{N1}& \mathrm{.Math.}& B_{\mathrm{NM}}\end{array}\right],𝕄\propto \left[\begin{array}{ccc}{m}_1& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ 0& \mathrm{.Math.}& m_M\end{array}\right],& \left(4\right)\end{array}$ b.sub.k are modelled magnetic fields on all the sensors, is an actuation matrix mapping magnetic moments to the magnetic field on each sensor, is a magnetic moment package for the actuators and i.sub.k is a measured electric current for all the actuators in a stacked form at the k-th measurement.

12. The method in accordance with claim 1, wherein B-field values are expressed as a stacked matrix form as
.sub.a(p.sub.a,i)∝[b.sub.1. . . b.sub.L]∝[i.sub.1. . . i.sub.L]∝  (5) where p.sub.a is a vector including all actuation parameters i is a stacked vector for all i.sub.k, and is a stacked matrix for all current measurements.

13. The method in accordance with claim 1, wherein the calibration model comprises a cost function which minimizes errors between the B-field values, wherein the cost function is formulated as $\begin{array}{cc}\underset{p_s.Math.p_a}{\mathrm{minimize}}{.Math.{\stackrel{\_}{𝕍}}_s\left(p_s,v\right)-{𝕍}_a\left(p_a,i\right).Math.}_F^2& \left(6\right)\\ \begin{array}{cc}\mathrm{subject}\mathrm{to}& p_{s,lb}⪯p_s⪯p_{s,\mathrm{ub}}\\ & p_{a,lb}⪯p_a⪯p_{a,\mathrm{ub}}\\ & g\left(p_s\right)=0,\end{array}& \end{array}$ where ∥⋅∥.sub.F is the Frobenius norm, p.sub.s are sensor parameters p.sub.s,lb and p.sub.s,ub are sensor parameters' lower and upper boundaries, respectively, p.sub.a are actuator parameters, and P.sub.a,lb and p.sub.a,up are actuator parameters' lower and upper boundaries, respectively.

14. The method in accordance with claim 1, wherein a temporary magnetic scale is added to the calibration model which is derived from $\begin{array}{cc}\mathrm{mean}\left(\frac{{.Math.G_i.Math.}_F^2}{3}\right)-g^2=0,& \left(8\right)\end{array}$ where g is a gain value of the sensors.

15. The method in accordance with claim 1, wherein the calibration model comprises an optimization, which is formulated as $\begin{array}{cc}\underset{\alpha ,r_r,{\hat{n}}_r}{\mathrm{minimize}}{.Math.\alpha {\stackrel{\_}{𝕍}}_s\left(p_s,v_r\right)-{\left[B_{1r}\mathrm{.Math.}{B}_{N_r}\right]}^T{\hat{n}}_{r}m.Math.}_2^2,& \left(9\right)\end{array}$ where α is a magnetic scale correction factor, v.sub.r is the magnetic field measurement data when the sample magnet is in the workspace level, B.sub.ir is formulated using relative positions r.sub.i-r.sub.r, r.sub.r is the position of the sample magnet, {circumflex over (n)}.sub.r is an orientation vector of the sample magnet and m is a known magnetic moment of the sample magnet.

16. A magnetic actuation and sensing system, comprising a plurality of magnetic actuators and a plurality of magnetic sensors, a respective output of the plurality of magnetic actuators and of the plurality of magnetic sensors being made available at a common evaluation unit, the evaluation unit simultaneously calibrating the magnetic actuators and the magnetic sensors in accordance with a method of claim 1 using a calibration model.

17. The magnetic actuation and sensing system according to claim 16, wherein the plurality of magnetic actuators comprises 3 to 25 magnetic actuators.

18. The magnetic actuation and sensing system according to claim 16, wherein the plurality of magnetic sensors comprises between 9 and 256 magnetic sensors.

19. The magnetic actuation and sensing system according to claim 16, wherein the number of magnetic sensor parameters and the number of magnetic actuator parameters are 27 N and 6 M, respectively, which formulates an inequality which, in the case of a triaxial sensor according to claim 4, can be expressed as 3NL+1≥27N+6M, which results in $\begin{array}{cc}L\ge 9+\frac{6M-1}{3N}.& \left(7\right),\end{array}$ where L represents a total number of measurements.

Description

(1) Further embodiments of the invention are described in the following description of the Figures. The invention will be explained in the following in detail by means of embodiments and with reference to the drawing in which is shown:

(2) FIG. 1: an exemplary experimental setup for a sensing and actuating system;

(3) FIG. 2: a table explaining the different parameters;

(4) FIG. 3: an exemplary explanation of the geometric parameters of the sensing and actuating systems;

(5) FIGS. 4A to 4F: exemplary simulation results for the experimental setup of FIG. 1,

(6) FIG. 5: a table showing a specification for a real application of a magnetic actuation and sensing system;

(7) FIG. 6: a table showing the initial parameters and tolerances for the application of FIG. 5;

(8) FIGS. 7A to 7B: experimental calibration results; and

(9) FIGS. 8A to 8D: a comparison of the model fitting accuracy using different actuation models.

(10) FIG. 1 shows the targeted calibration system 10, which is composed of a magnetic sensing system 12 with a plurality of magnetic sensors 14, a workspace 16 for the objects (or patients) to be measured (not shown) as well as a magnetic actuation system 18 with a plurality of magnetic actuators 20.

(11) The magnetic sensing system 12 lies within its sensing system level 22 whereas the magnetic actuation system 18 lies within an actuation system level 24. Both levels 22, 24 are shown by an exemplary dashed line. The sensing system level 22 is spaced apart from the actuation system level 24, with the workspace being arranged between the two.

(12) The magnetic actuators 20 are composed of electromagnets which each comprise a coil 26 and core 28. A respective current sensor (not shown) can be attached to the magnetic actuators 20.

(13) The outputs of the plurality of magnetic actuators 20, the plurality of magnetic sensors 14 and the current sensors is made available at a common evaluation unit 30. The evaluation unit 30 is configured to simultaneously calibrate the magnetic actuators 20 and the magnetic sensors 14 using a calibration model.

(14) The magnetic actuation and sensing system 10 can be calibrated simultaneously by generating a plurality of arbitrary magnetic fields at the workspace 16 by applying a plurality of electric currents at the magnetic actuators 20.

(15) Each of the respectively generated magnetic fields is measured with the sensing system 12 in the workspace 16 and each respective electric current is measured with the current sensor at each of the magnetic actuators 20.

(16) A magnetic field of a random known sample magnet (not shown) is also measured with the sensing system 12. For this the sample magnet is placed at an arbitrary position and orientation in the workspace 16.

(17) The measured data is then fed into a calibration model, with the calibration model being based on a sensor measurement model in combination with a magnetic actuation model. Using the calibration model a solution of the model parameters is found via a numerical solver in order to calibrate both the actuation and sensing systems at the same time.

(18) Thus, the uncalibrated sensors and electromagnets work as calibration sources. Although the parameters of the system are inaccurate prior to the calibration, all the parameters for the sensor and actuation systems are calibrated simultaneously by the property of the calibration model. For this it is assumed that the current sensors are pre-calibrated before the calibration and the positions of the magnetic sensors are known.

(19) It should be noted that each of the plurality of generated magnetic fields is non-uniform throughout the workspace 16, which is generated from a solenoidal electromagnet.

(20) It should further be noted that the calibration model may comprise a parameter set, which includes sensor parameters and actuator parameters, wherein the sensor parameters comprise parameters for position, orientation, axis distortion, axes distortion and gain of each sensor, and wherein the actuator parameters comprise parameters for position, orientation and magnetic moment of each actuator.

(21) FIG. 2 shows a table explaining all the parameters which are necessary for the simultaneous calibration of the sensing and actuation system 10, which will be described later.

(22) FIG. 3 shows a geometric description of some of the parameters used for calibration. All the position vectors and rotation matrices are expressed in a global frame {G}. The local frame {L} is used to express gains of the sensors, m.sub.j and m′.sub.j show magnetic moments of j- and j′-th actuators 20, respectively. Black dots show the shifted center of the magnetic sources due to mutual inductions.

(23) FIGS. 4A to 4F show simulation results from a simultaneous calibration, wherein FIGS. 4A and 4B show rendered images of configurations of the OctoMag- and a large-scale OctoMag-system. The hollow cylinders stand for the coils 28, whereas the inner cylinders stand for the cores 26. The sphere in the middle of the Figures shows an exemplary workspace 16, whereas the sensing system 12 is shown via a plate 12 below the workspace 16.

(24) FIGS. 4C to 4D and FIGS. 4E to 4F show normalized error distributions of the sensing system and actuation system, respectively, after the calibration. The Figures show that the standard deviation (STD) of the normalized errors is less than 0.45% in all the cases. The sensing and the actuation systems' coefficients of determination (R.sup.2 values) are over 99.88% in all the cases.

(25) FIGS. 5 and 6 show further tables which reveal the parameters and tolerances of a real application of the simultaneous calibration of a sensing and actuation system.

(26) FIGS. 7A and 7B show the experimental calibration results, wherein FIG. 7A shows the sensing system calibration result which resulted from using the method according to the invention in contrast to results from a known state of the art individual calibration of a sensing system. The state of the art method showed a R.sup.2 value of 99.57% with 0.97% STD of the normalized error. The method according to the invention on the other hand shows a R.sup.2 value of 99.84% with 0.59% STD of the normalized error.

(27) In FIG. 7B one can see the result of an actuation system calibration by using the method according to invention in contrast to the result of a known state of the art individual calibration. The state of the art calibration shows a R.sup.2 value of 99.91% with 0.19% STD of the normalized errors, while the method according to the invention shows a R.sup.2 value of 99.45% with 0.47% STD of the normalized errors.

(28) FIGS. 8A to 8B show a comparison of model fitting accuracy using different actuation models. In FIG. 8A, one can see an experimental setup, where a coil with a core and an extra core are used. The workspace 16 is shown by the sphere with a radius of r.sub.w. The workspace 16 is placed above the closest coil 26 with a center-to-center distance d. The length of the cores 28 is d.sub.c.

(29) In this connection it should be noted that models to which the current invention is compared are respective models used for only calibrating either the sensors or the actuators, but neither are calibrated together at the same time using one model as is presently the case in the present invention.

(30) FIG. 8B shows the actuation model fitting results using different orders of actuation models by varying the workspace's range and distance. For each workspace configuration, three dipoles, three dipoles with three quadrupoles and three dipoles with three quadrupoles and three hexapoles are used. The model actuation model fitting results using a dipole model together with up to different orders of multi-pole models and the results using different numbers of dipole sources are shown in the FIGS. 8C and 8D, respectively.

(31) The simultaneous calibration method is described in more detail below on the basis of a specific example.

A. System Description

(32) The targeted calibration system is composed of a plurality of magnetic sensors (sensor array) and a set of electromagnets. The final application of the system is to control a magnetic robot using the actuation system (electromagnets) while the magnetic robot is localized by a localization algorithm using sensor data.

(33) The sensor array measures magnetic fields from a magnetic robot and the electromagnets. The magnetic sensor can be either mono-axial or multi-axial. In described example a triaxial sensor is used as that sensor is harder to calibrate. With the coupled magnetic fields from the magnetic robot and the actuators, the robot can be localized by a known algorithm.

(34) The actuation system generates desired magnetic fields and their spatial gradients at the center of the magnetic robot. The magnetic field is non-uniform throughout the workspace. Systems using non-uniform magnetic fields from solenoidal electromagnets, such as OctoMag or MiniMag, could be categorized as an applicable system. The example system (FIG. 1) has 64 sensors and 9 electromagnets.

B. Framework: Calibration as Bundle Adjustment

(35) Bundle adjustment (BA) has been employed in online simultaneous localization and mapping (SLAM) and offline structure from motion (SfM), which enables precise recovery of the parameters of both the sensors and landmarks. It is essentially an optimization-based method that can handle the simultaneous estimation of large numbers of parameters. It can be robust to external noises in the measurements. In the case where the measurement noise is zero-mean and normally distributed, BA solves a maximum likelihood estimation problem.

(36) In the method according to the invention, the measurements of the electric currents and the magnetic field can be treated equivalently as the feature points of images in BA. The positions, orientations and gains of the actuators and sensors are similar to the parameters of landmarks and camera(s). These parameters are constrained by the magnetic field dipole model under specific conditions.

(37) In order to solve the calibration problem in the BA framework, it must be considered parametrization, error modeling, building constraints, removing gauge freedom and choosing a proper optimization process.

C. Parametrization

(38) The proposed method involves a large set of parameters (the example system in the experiment has 1836 parameters to be calibrated). In such a large parameter space, an appropriate parametrization is required for a reliable optimization result. The parameter set is composed of sensor parameters and actuator parameters. The systems of the parameters are represented in FIGS. 2 and 3.

(39) 1) Sensor Parameters: A sensor has parameters for position, orientation, axis distortion, axes distortion and gain (either linear or non-linear). For a general sensor calibration, each sensor is regarded as having a non-linear response and a cross-coupling effect among the axes. The non-linearity is attributed by hard-iron or soft-iron effect from ferrous material on the circuit board. The distortion of axes is attributed by fabrication errors from the manufacturing process or thermal stress from the soldering process. In the case of an extremely sensitive triaxial magnetometer, such as a fluxgate magnetometer, the transverse field effect (TFE) becomes an additional source for the non-linearity and cross-coupling effects.

(40) A couple of models are presented to capture the non-linearity and the axes distortion. These models are for single triaxial magnetometers, so the axis rotation is neglected assuming that the axis rotation is pre-calibrated. However, in general, a system with multiple sensors might have a different orientation for each sensor. This is because the fabrication process might change the orientations or different orientations can be intended to increase the signal to noise ratio to a certain direction. Involving the non-linearity, cross-coupling and sensor rotation, a 3-D polynomial model could be used. The model can be easily expanded by increasing the order of the polynomials with the expense of computation.

(41) In this example a 3-D quadratic model is used. This model is not only able to capture the non-linearity of the sensors precisely, but is also computationally practical. The 3-D quadratic model with sensor orientations involves 27 parameters for each sensor. It is assumed that the positions of the sensors are known. The experimental setup has 64 sensors, thus 1782 parameters for sensors are presented.

(42) 2) Actuator Parameters: A magnetic actuator has a position, an orientation and a magnetic moment. A typical electromagnet has a coil and core and their axes are aligned together. In such a case, the orientation and the magnitude of the magnetic moment of the electromagnet can be combined into a 3-D vector. A system with more than two cores has a mutual induction which might change the magnetic center's position and the magnetic moment's orientation compared to a single electromagnet presented. The experimental setup has 9 actuators, thus 54 parameters for the actuators presented.

D. Calibration Model

(43) The calibration model is based on a sensor measurement model and magnetic actuation model. The magnetic fields from the actuators are measured by the magnetic sensors and the values are expressed in arbitrary units. Additionally, the electric currents are measured by the current sensors attached to the actuators. It is assumed that the current sensors are the only reliable sensors with zero-mean Gaussian noise corruption.

(44) 1) Sensor Measurement Model: The magnetic field is measured by the magnetic sensor array. A single sensor reads the magnetic flux density or the B-field through the voltage difference across the ends of a magneto-sensitive resistor. It is assumed that the sensor is non-linear, of which the axes are skewed and rotated, and the gains are different from the factory specifications. Using a quadratic model, the recovered B-field in the i-th sensor at the k-th measurement is

(45) ${\stackrel{\_}{b}}_{\mathrm{ik}}\left(p_{\mathrm{si}},v_{\mathrm{ik}}\right)=G_i{v}_{\mathrm{ik}}+\left[\begin{array}{ccc}{v}_{\mathrm{ik}}^T& H_{x,i}& v_{\mathrm{ik}}\\ v_{\mathrm{ik}}^T& H_{y,i}& v_{\mathrm{ik}}\\ v_{\mathrm{ik}}^t& H_{z,i}& v_{\mathrm{ik}}\end{array}\right],$
where p.sub.s,i is an array of parameters for G.sub.i, H.sub.y,i, and H.sub.z,l. In this expression G.sub.i, is a linear mapping with no additional constraints, meaning it represents the effects of scaling, rotation and shear. H.sub.x,i, H.sub.y,i, and H.sub.z,i are symmetric 3×3 matrices, which involve rotations and quadratic terms together. Note that the other symbols are explained in FIG. 2.

(46) As each sensor in the sensor array has multiple measurements, it is convenient to express the recovered B-field values from the whole measurements in a stacked matrix form as

(47) 0${\stackrel{\_}{𝕍}}_s\left(p_s,v\right)=\left[\begin{array}{ccc}{\stackrel{\_}{b}}_{11}\left(p_{s,1},v_{11}\right)& \mathrm{.Math.}& {\stackrel{\_}{b}}_{1L}\left(p_{s,1},v_{1L}\right)\\ \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ {\stackrel{\_}{b}}_{N1}\left(p_{s,N},v_{N1}\right)& \mathrm{.Math.}& {\stackrel{\_}{b}}_{\mathrm{NL}}\left(p_{s,N},v_{\mathrm{NL}}\right)\end{array}\right],$
where p.sub.s and v are stacked vectors for all p.sub.s,i and v.sub.ik, respectively.

(48) 2) B-field Actuation Model: The system described in FIG. 1 can be effectively modelled by the dipole magnetic field. This assumption is useful when a large-scale optimization is involved because the dipole model avoids intensive computations, but gives accurate estimation. Methods using the magnetic dipole model to fit the magnetic field for a workspace are already known. The known method uses a magnetic dipole to fit the magnetic field of a unit-current contribution. By increasing the number and the order of magnetic poles, more accurate model could be achieved. Due to the large search space in the optimization, it is preferable to have the least amount of parameters. In the current problem, a single dipole model for each unit-current contribution is chosen. A detailed analysis on the accuracy of the model with respect to the model complexity is present in the following.

(49) The dipole magnetic field on the i-th sensor form the j-th actuator at the k-th measurement is expressed as

(50) $b_{\mathrm{ij},k}=B_{\mathrm{ij}}{m}_j{I}_{\mathrm{jk}};$$B_{\mathrm{ij}}\equiv \frac{{\mu }_0}{4\pi {.Math.r_{\mathrm{ij}}.Math.}^3}\left(3{\hat{r}}_{\mathrm{ij}}{\hat{r}}_{\mathrm{ij}}^T-𝕀\right),r_{\mathrm{ij}}\equiv r_i-r_j,$
where b.sub.ij,k is the magnetic field vector in R.sup.3, B.sub.ij is a 3×3 matrix which only depends on the displacement vector r.sub.ij, .{circumflex over ( )}. is a unit vector normalization operator, is the identity matrix and μ.sub.0 is the permeability of free space.

(51) For magnetically actuated robots, multi-degrees-of-freedom (DOF) actuation is preferred as it can give more authority to control the robot and avoid singularities or reduce the energy consumption using the redundancy. Because the magnetic field is superimposed at a point, it can be considered the sum of the magnetic fields from each actuator as the resultant magnetic field.

(52) The magnetic field at multiple measurement points from multiple current inputs could be expressed in a matrix equation as b.sub.k=i.sub.k, where

(53) $𝔹=\left[\begin{array}{ccc}{B}_{11}& \mathrm{.Math.}& B_{1M}\\ \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ B_{N1}& \mathrm{.Math.}& B_{\mathrm{NM}}\end{array}\right]\mathrm{and}𝕄=\left[\begin{array}{ccc}{m}_1& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ 0& \mathrm{.Math.}& m_M\end{array}\right],$
b.sub.k is a modelled magnetic field on all the sensors in a stacked matrix form, is the actuation matrix from magnetic moments to the magnetic field on each sensor, is the magnetic moment package for the actuators (diagonal: magnetic moments, off-diagonal: zero-vectors) and i.sub.k is a measured electric current for all the actuators in a stacked form at the k-th measurement.

(54) Similar to the sensor measurement model, it is convenient to express the B-field values into a stacked matrix form as .sub.a(p.sub.a, i)=[b.sub.1 . . . b.sub.L]=[i.sub.1 . . . i.sub.L]=, where p.sub.a is a vector including all the actuation parameters, i is a stacked vector for all i.sub.k and is a stacked matrix for all the current measurements.

E. Formulation of a Cost Function

(55) The cost function can be formulated by minimizing the errors between the B-field values given by the sensor model and the actuation model. The problem can be stated in the least square manner,

(56) $\underset{p_s.Math.p_a}{\mathrm{minimize}}{.Math.{\stackrel{\_}{𝕍}}_s\left(p_s,v\right)-{\stackrel{\_}{𝕍}}_a\left(p_a,i\right).Math.}_F^2$$\begin{array}{cc}\mathrm{subject}\mathrm{to}& p_{s,lb}⪯p_s⪯p_{s,\mathrm{ub}}\\ & p_{a,lb}⪯p_a⪯p_{a,\mathrm{ub}}\\ & g\left(p_s\right)=0\end{array}$
where ∥⋅∥.sub.F is the Frobenius norm, p.sub.s are sensor parameters, p.sub.s,lb and p.sub.s,ub are sensor parameters' lower and upper boundaries, respectively, p.sub.a are actuator parameters and p.sub.a,lb and p.sub.a,up are actuator parameters' lower and upper boundaries, respectively; is an element-wise (general) inequality symbol. These boundaries come from the manufacturing tolerance or estimation tolerance for each parameter. g(⋅) is an additional constraint to resolve the gauge freedom, which is described in (8).

(57) 1) Minimum Number of Measurements: The total number of the equations should be larger than the total number of parameters to solve the optimization problem. The total number of the equations is a multiplication of the number of the sensor elements, 3N (in the case of a triaxial sensor), and the total number of measurements L. Additionally, there is an additional constraint for the optimization, g(⋅), which gives one more equation. The total number of parameters are the sum of the number of the sensor parameters, 27 N, and the number of actuator parameters, 6 M. This formulates an inequality, 3NL+1≥27N+6M, which results in

(58) $L\ge 9+\frac{6M-1}{3N}.$

(59) Although the minimum number of the measurements is defined, it is beneficial to collect measurement data more than the minimum number to improve calibration accuracy. The large number of data minimizes Gaussian error from the sensor system effectively through the optimization process.

(60) 2) Resolving Gauge Freedom: There is a gauge freedom in the magnetic scale. The measurements are invariant under a scaling of either the norm of one of the G.sub.i values. To resolve the gauge freedom, a physical reference needs to be introduced, which can give an absolute value for the magnetic scale.

(61) The physical reference is introduced using a permanent magnet whose magnetic moment is known. By measuring the magnetic field of the permanent magnet inside the workspace, the magnetic moments of the actuators and the gain for the sensors are related with the physical values.

(62) A “temporary” value for the gauge freedom parameter would be needed for a better calibration. This is because the measurement of the magnetic field of the permanent magnet is not redundant as the measurement of the magnetic fields form the electromagnets (number of the sample is much smaller and the noise is not compensated). The temporary value can be corrected after the first calibration. Although the measurement of the magnetic field form the reference magnet is a one-time measurement, it is preferred to expose the sensors to the magnetic field of the magnet for a long period of time and take a mean value of the measured data. This proves would minimize the effect of Gaussian noise during the measurement of the reference magnet.

(63) The temporary magnetic scale is added on the cost function. The magnetic scale can be fixed by constraining either the sensor gain or the magnetic moment of the actuator to a certain number. The sensor gain might be close to the real value as it is specified in the manufacturer's data sheet. Fixing the mean value of the gain of the sensors as the best pre-knowledge (magnetometer specification data) as

(64) $\mathrm{mean}\left(\frac{{.Math.G_i.Math.}_F^2}{3}\right)-g^2=0$
where g is the gain value of the sensors in the specification sheet. Adding this to the constraint in the scalar function g(⋅) of (6) finalizes the formulation for the optimization task.

(65) 3) Reassignment of Global Reference: The gauge freedom is fixed to the “temporary” value in the previous section. The value should be corrected according to the absolute physical value. This procedure should be performed after the first main optimization task.

(66) The reference reassignment proves involves additional optimization; however we use the result parameters from the first optimization. We can formulate the optimization task to find the magnetic scale of the system compared with the reference magnet. The optimization is formulated as

(67) $\underset{\alpha ,r_r,{\hat{n}}_r}{\mathrm{minimize}}{.Math.\alpha {\stackrel{\_}{𝕍}}_s\left(p_s,v_r\right)-{\left[B_{1r}\mathrm{.Math.}{B}_{N_{\tau }}\right]}^T{\hat{n}}_{r}m.Math.}_2^2,$
where α is the magnetic scale correction factor, v.sub.r is the magnetic field measurement data when the reference magnet is in the workspace, B.sub.ir is formulated using relative positions, r.sub.i-r.sub.r, r.sub.r is the position of the reference magnet, {circumflex over (n)}.sub.r is the orientation vector of the reference magnet and m is the known magnetic moment of the reference magnet.

(68) Using α from the optimization, the magnetic scales can be assigned as G*.sub.i=αG.sub.i, H*.sub.i=αH.sub.i, m*=αm, where the ‘*’ symbol means a corrected value.

F. Choosing a Numerical Solver

(69) The interior-point algorithm (IPA) is chosen to solve the bundle adjustment problem. Because the non-linear optimization problem stated in (6) has both equality and inequality constraints and a high-dimensional parameter space to optimize over, an optimization solver is needed which can handle the constraints well.

Simulation Results

(70) The proposed method was tested in a simulation environment, where all the system parameters are known as the ground-truth for the ease of comparison. As a general solenoid type magnetic actuation system, OctoMag and its large-scaled version is chosen as simulation models (configurations are shown in FIGS. 4A and 4B). Those systems are modelled in a commercialized finite-element-analysis (FEA) software (COMSOL™, Comsol group).

(71) The sensor system was simulated based on the data sheet of the sensor for the experiments. The sensor array has sixteen sensors in 2-D grid and the sensor array was placed tangent to the workspace. The placement of the sensor array is to simulate a system which uses the sensor array as a localization system, where the sensor array does not disturb the movement of the magnetic robot, but close to sense the magnetic field from it.

(72) Form the ground-truth configuration, random magnetic fields were generated with random electric current inputs. The random electric currents and magnetic fields were measured with virtual electric current sensors and magnetic sensors with Gaussian noises. The magnetic sensors have 900 nT standard deviation (STD) Gaussian noise, and the current sensors have 10 mA STD Gaussian noise. These measurements were fed into the optimization task. The parameter boundaries were 5% offset from the ground-truth values, except the magnetic moments as they are hard to predict before the calibration. From the given configuration, there were 675 parameters and 1 non-linear constraints in total. For the termination condition of the optimization, the first-order optimality tolerance was set to the noise level of the sensors.

(73) The validation data set is generated separately to check the accuracy of the model. The actuation model validation set is the B-field set inside the workspace, when coils are running with arbitrary inputs. The sensing model validation set is the B-field set at the place of the sensor array with different current inputs tan training datasets.

(74) FIGS. 4A-4F shows the simulation result. After the calibration, the sensor and actuation model had shown low normalized errors. STD of the normalized error of the sensing and actuation systems were 0.35%, 0.45% for OctoMag system, respectively, and 0.59% and 0.062% for the large-scale OctoMag system, respectively. The determinants of coefficient (R.sup.2 value) for all cases were over 99.88%.

Experimental Results

(75) A. Experimental setup and procedure: The proposed method was applied to a custom-designed magnetic localization and actuation system 10 (FIG. 1). The system consists of nine electromagnets 20 and sixty-four triaxial magnetic sensors 14 located on an 8×8 grid with 30 mm nodal distance. The sensor array was placed 250 mm above the central electromagnet in z-direction and their centers were aligned on the xy-plane. The system specifications are shown in FIG. 5.

(76) Manufacturing errors are present in both the localization and actuation systems. The sensing system involves errors from the orientation and the sensitivity of each sensor. The electromagnet contains a position shift due to the mutual inductance, and a model error from the magnetic moment mapping. All of these design parameters and manufacturing tolerances are shown in FIG. 6. Note that the initial parameters define the initial conditions, and tolerances defines the boundary conditions for the numerical solver.

(77) Before the calibration procedure, the DC magnetic field around the system and the current sensor offset are nullified. The sensor system's values are measured for a period of time (20 seconds) and their mean values are set as nullifying offsets. The magnetic field noise floor was 50 μT and current sensor offset was ±5 mA.

(78) In the experiment, the electromagnet generated five-hundred random magnetic fields (L=500) from current inputs ranging in [−3,3] A. The currents were measured by the current sensors while the magnetic field was measured by each magnetic sensor independently. After the measurements, the magnetic field of the reference magnet was measured inside the workspace in an arbitrary orientation and position. All the measurement data were fed into (6) and the optimization task was solved by the interior-point algorithm. The data measurements and the optimization took three minutes and fifteen minutes, respectively. After the optimization, the parameters were recorrected to fix the gauge freedom using (9) and (10).

(79) B. Magnetic Sensor System Calibration Result: The sensor measurement was compared with the ground-truth data. A triaxial Gauss meter was used for the ground-truth data measurement. To generate magnetic fields to be measure, a set of random current inputs ran in the coils. The random current inputs were generated using a uniform distribution within [−1, 1] A. The input currents should be within a proper range, such that the magnitude of the magnetic field is linear to the input currents and the magnetic field does not saturate the magnetic sensors. The magnetic field was measured first by the sensor array. After that, the Gauss meter measured the magnetic field at the same location of each sensor with the same current input. In total 1350 data points were measured for both the sensor array and the Gauss meter.

(80) As a reference method, we applied the calibration method from Plotkin et al. to the sensor system calibration. The method uses a permanent magnet to calibrate a magnetic sensor system. The method is modified to calibrate triaxial magnetic sensors, whereas the original method is for single axis magnetic sensors. ϕ 6×6 mm.sup.3 NdFeB N52 magnet was used as the calibration source. The magnet was placed in hundred different locations and its magnetic field was measured by the sensor array.

(81) FIG. 7A shows the sensor system calibration result. Using the proposed method, the sensor model had R.sup.2 value of 99.84% and the normalized STD of 0.59%. The calibration result from the reference method resulted in R.sup.2 value of 99.57% and a normalized error STD of 0.97%. The result shows that the proposed method is comparable with the previous sensing system calibration method.

(82) C. Actuation System Calibration Result: The actuation system was calibrated simultaneously by the proposed method. To validate the calibration results, the modelled magnetic field from (4) was compared with the measured magnetic field using the Gauss meter. The input currents for the electromagnet were generated randomly within [−1, 1] A using a uniform distribution. The errors where calculated by the difference between the measured magnetic field and the modelled magnetic field. The Gauss meter was placed in a 50×50×50 mm.sup.3 cubic grid which fit into the workspace. In total 1350 data points were measured and compared.

(83) The result of the proposed method is compared with the result from the method of Petruska et al. as a reference. To collect the data for the reference method, the Gauss meter measured B-field directly inside the workspace when each coil ran −3 A, 0 A and 3 A. These data were fit into spherical harmonics model with one B-coefficient for each magnetic source.

(84) FIG. 7B shows the experimental results. The proposed method resulted in the actuation model with R.sup.2 value of 99.45% and the normalized error STD of 0.47%. The known calibration method resulted in R.sup.2 value of 99.91% and a normalized error STD of 0.19%. The result shows that the proposed method is comparable with the previous actuation system calibration method.

Discussion

(85) A. Actuator model's complexity and accuracy: More complex models for actuators, such as multipole expansion models using spherical harmonics, could also be potentially used in the simultaneous calibration method. A magnetic dipole field model was used in this experiment because the model is computationally tractable and yields accurate results in the tested systems. For these systems, a more complex model may improve the actuator calibration accuracy but will increase the computational costs and may result in over-fitting.

(86) The model fitting accuracy was compared when using different actuator models, in terms of the number and order of the magnetic poles by varying the size of the workspace and the distance from magnetic sources to the workspace. A set of control experiments was conducted in a simulated environment using finite-element-analysis (FEA). FIG. 8A shows the simulated configuration where a coil 26 with a core 28 is placed near to another core 28. The workspace 16 dimensions, the bounding sphere radius r.sub.w and the distance between the center of the coil 26 and the center of the workspace 16, d, were varied and the magnetic field data from FEA was obtained when a unit-current flows through the coil 26. The magnetic field data were fitted with actuation models using already known methods.

(87) FIGS. 8A-D shows the result of the actuation model fitting. In FIG. 8B, the general trend shows all multi-pole models work better when increasing the distance between the workspace and the magnetic sources. With a fixed number of magnetic sources (three sources), even the dipole model yields an STD value less than 1% in all ranges. FIGS. 8C and 8D show the accuracy of the model when the size and position of the workspace is fixed (r.sub.w/d.sub.c=0.375 and d/d.sub.c=1.75). In FIG. 8C, with a fixed number of sources (three sources), by increasing the order of the poles, the accuracy of the model is improved. In FIG. 8D, the actuation model improves with an increased number of sources. In all cases, the simplest model (with a single dipole) fits the B-field in the workspace less than 1% STD of the normalized errors.

(88) B. Placement of sensor arrays: The sensor array placement might affect the accuracy of the calibration result. If the sensor is too far away from the workspace, the sensor array could not observe the magnetic field around the workspace. In such a case, the actuation model on the sensor array is accurate, but might not in the workspace. To achieve the best result, it is preferred to place the sensor array inside the workspace. Then the magnetic field observed in the sensor array closely matches with the magnetic field of the workspace.

(89) The sensor array could be placed near by the workspace as well. This is important for the applications where it is hard to place the sensor array inside the workspace due to mechanical constraints. Additionally, a system needs to be automated and calibrated regularly would be one of the applications. In such a case, it is recommended to place the sensor array as close as possible to the workspace, so that the B-fields on the sensors are close to the B-field of the workspace. Note that the simulations and experiments were conducted to simulate those cases where the sensor array is near to the workspace.

(90) C. Applications of the method: There are two major applications of the method according to the invention. First, the method could be employed for an automatic calibration with minimal human intervention. This would be beneficial for a system which needs recalibration often. For example, magnetic sensors require gain updated due to their time-variant soft-iron effect or a long-term gain-offset drift from the environmental change, such as the Earth's magnetic field change. Depending on the time scale of the gain change of the system or environment, the update can be scheduled regularly. A user could expose the reference magnet to the sensor system before the calibration, and the rest of the calibration could be done automatically.

(91) Second, the proposed method allows the automatic calibration for reconfigurable magnetic sensing and actuation systems. The reconfiguration may be beneficial for improving dexterity, energy efficiency or avoiding singularities of an actuation system. Especially, the reconfiguration of the system would be critical in medical robotics, because a treatment of a patient might have a preferred direction. In the case of the reconfiguration, the parameter space could be minimized including only the actuation parameters enabling in-Situ calibration.