Finite Element Modelling of a Field-Sensed Magnetic Suspended System for Accurate Proximity Measurement Based on a Sensor Fusion Algorithm with Unscented Kalman Filter

The presented paper describes accurate distance measurement for a field-sensed magnetic suspension system. The proximity measurement is based on a Hall effect sensor. The proximity sensor is installed directly on the lower surface of the electro-magnet, which means that it is very sensitive to external magnetic influences and disturbances. External disturbances interfere with the information signal and reduce the usability and reliability of the proximity measurements and, consequently, the whole application operation. A sensor fusion algorithm is deployed for the aforementioned reasons. The sensor fusion algorithm is based on the Unscented Kalman Filter, where a nonlinear dynamic model was derived with the Finite Element Modelling approach. The advantage of such modelling is a more accurate dynamic model parameter estimation, especially in the case when the real structure, materials and dimensions of the real-time application are known. The novelty of the paper is the design of a compact electro-magnetic actuator with a built-in low cost proximity sensor for accurate proximity measurement of the magnetic object. The paper successively presents a modelling procedure with the finite element method, design and parameter settings of a sensor fusion algorithm with Unscented Kalman Filter and, finally, the implementation procedure and results of real-time operation.


Introduction
Linear proximity sensors (LSPs) with mid-and low-range measurement capabilities are devices that are used widely in many industrial and non-industrial applications. They are mostly used to determine the displacement, direction of movement, orientation, speed, etc., of a measured object. The LSP exploits different physical principles of operation, where capacitive, inductive, ultrasonic, optical and magnetic phenomena are the most commonly deployed in a sensing operation. Many of these physical phenomena, especially accurate optical technology and the ultrasonic principle, require complex pre-and post-processing operations which, unintentionally, result in a high price and relatively large dimensions of the measuring unit. The dimensional obstacles and the high price of the sensor devices often prevent installation of precise sensing technology on small/miniature and low cost devices. In the time of high expansion and pervasive sensing technologies, especially in the field of miniature sensors, as well as a highly efficient processing unit of relatively small dimensions and price, they offer many applicable solutions which can effectively replace many complex and expensive Figure 1 shows electromagnet-1 with driven current i, linear Hall Effect sensor-2 and actuated body with permanent magnet-3. Variable d is the distance between the actuated body and the electromagnet. The following research will not be relayed directly on the specified structure (a) or (b) of the EMAwS depicted in Figure 1, but on the modelling of the magnetic force F e for the proximity measurement d. The modelling of the magnetic force applied to the magnetic body is principally the same for systems (a) and (b) in Figure 1. For the electromagnetic force modelling and for further data fusion algorithm the state space description of the system is taken.
Sensors 2016, 16,1504 4 of 23 (a) or (b) of the EMAwS depicted in Figure 1, but on the modelling of the magnetic force for the proximity measurement . The modelling of the magnetic force applied to the magnetic body is principally the same for systems (a) and (b) in Figure 1. For the electromagnetic force modelling and for further data fusion algorithm the state space description of the system is taken. The state space representation of the system is: ( ) ∈ ℜ and ( ) ∈ ℜ represent smooth nonlinear uncertainty functions and they are assumed to be bounded. With regard to the law of motion the model structure is: and the coil dynamic equation is: where ( ) , ( ) , , , ( ), , , ( ) , ( , , ) are the applied voltage, coil current, winding resistance, winding inductance, induction constant, distance, body mass, opposite force and magnetic force, respectively. The meaning of opposite force depends on the system structure in Figure 1 and can be treated as gravity force = for system (a) and the friction force = ℎ for the system (b). The parameter is the friction coefficient. Calculation of the magnetic force between two magnets is a complex task and can be done with the Gilbert and Ampere model [29]. With regard to the mentioned models the magnetic force between the actuated body and the electromagnet can be determined with inspection of the magnetic field as a function of the separation distance presented by Naumović and Hajjaji [22,23,29,30]. It can be done with the integral method, where the solenoid is modelled as a volume current density and the permanent magnet as a surface current density around its circumference presented by Furlani and Robertson [31,32]. The derived magnetic force models are mostly highly complex functions, which can be used only for offline analysis, parameter estimation and structure optimization [32]. Such models have very limited use in real-time systems, especially The state space representation of the system is: y (n) (t) = g (x (t)) + f (x (t)) u (t) (1) where the terms y (t) = [y 1 (t) , y 2 (t) , . . . , y m (t)] ∈ m , x (t) = [x 1 (t) , x 2 (t) , . . . , x m (t)] ∈ mn and u (t) = [u 1 (t) , u 2 (t) , . . . , u m (t)] ∈ m denote the system output vector, states vector and system input vector, respectively. The functions g (t) ∈ m and f (t) ∈ m represent smooth nonlinear uncertainty functions and they are assumed to be bounded. With regard to the law of motion the model structure is: and the coil dynamic equation is: where U (t), i (t), R, L, k m (x) , x, m, F f g (t), F e (x, i, t) are the applied voltage, coil current, winding resistance, winding inductance, induction constant, distance, body mass, opposite force and magnetic force, respectively. The meaning of opposite force depends on the system structure in Figure 1 and can be treated as gravity force F g = mg for system (a) and the friction force F f = k . h for the system (b). The parameter k is the friction coefficient. Calculation of the magnetic force between two magnets is a complex task and can be done with the Gilbert and Ampere model [29]. With regard to the mentioned models the magnetic force F e between the actuated body and the electromagnet can be determined with inspection of the magnetic field as a function of the separation distance presented by Naumović and Hajjaji [22,23,29,30]. It can be done with the integral method, where the solenoid is modelled as a volume current density and the permanent magnet as a surface current density  [31,32]. The derived magnetic force models are mostly highly complex functions, which can be used only for offline analysis, parameter estimation and structure optimization [32]. Such models have very limited use in real-time systems, especially in the systems with fast execution demand and the systems with modest computational power. In the presented research the magnetic force model will be used in the UKF algorithm in the systems with relativity short execution time demand and feedback controller. Algorithms, UKF and feedback controller need to be executed in each sampling period iteratively. The main focus of the modelling is obtaining an accurate model dynamic with an adequately simple structure.

Magnetic Force Modelling with the Finite Element Method
Partial nonlinear equations arise in mathematical modelling in many diverse areas, such as material science, fluid dynamics, electromagnetism, economics, etc. [33]. In most cases the equations of the described system are so complicated that their solution in purely analytical form is impossible or impractical. Due to the complexity, many times the circumstances have compelled users to search for approximate solutions to the unknown analytical solutions. The FEM is one of the numerical techniques used for finding approximate solutions of partial differential equations with known boundary values [34]. In the presented paper the FEM technique is used for modelling and analysis of the electromagnetic force F e and EMAwS system. The open source 2D Finite Element Method Magnetics (FEMM) created by Meeker [35] was used. FEMM offers broad coupling with different external simulation and analysis software and self-created scripts in LUA-language. Such external coupling and self-created scripts' capability allows many options in simulation, analysis and measurement. Each created model's parameters, dimensions, electrical circuit properties, etc. in FEMM can be controlled by coupled external software or the self-created scripts in LUA. We used the coupling possibility with MATLAB software, where MATLAB scripts are controlling the multiphysics model in FEMM, similar to the work presented by Benamimour [36].
The analysis of the electromagnetic force F e was started by drawing the EMAwS system in FEMM. The model drawing is an important task, where all material characteristics, geometric and circuit parameters, measurement units, mesh polygonal angle (grid generation) are determined. The system (a) in Figure 1 will be examined for further analysis. The system parameters are recorded in Table 1 and in Figure 2. All the materials, wire characteristic, number of coil turns, driven current, PM characteristic were assigned in the FEMM software as depicted in Figure 3. The picture in Figure 3 represents a 2D-EMAwS model in central cross-section view, with computational boundary.  In Figure 3. with the selection of solver precision (1 ) and minimal angle value setting (30˚) in a millimeter scale for the inside computational boundary surface, we get a computational mesh with 17,257 nodes. The selection of the solver parameters (solver precision and mesh angle) is crucial in the accuracy of the simulation results. In our approach the solver parameters were selected experimentally in a manner to acquire accurate and reliable results in a reasonable time. By increasing the solver precision and angle value of the FEM procedure, the computational effort was greatly increased, wherein the accuracy of the simulation was insignificant compared to the computation time used. In regard to the drawn model and solver parameters' selection, the simulation lasted approximately 4 min on a Windows 10-based PC, with a i7-3770@3.4GHz CPU and 8GB-RAM memory. The computational mesh in the selected computational boundary surface and calculated quantities with FEM are depicted in Figure 4.   In Figure 3. with the selection of solver precision (1 ) and minimal angle value setting (30˚) in a millimeter scale for the inside computational boundary surface, we get a computational mesh with 17,257 nodes. The selection of the solver parameters (solver precision and mesh angle) is crucial in the accuracy of the simulation results. In our approach the solver parameters were selected experimentally in a manner to acquire accurate and reliable results in a reasonable time. By increasing the solver precision and angle value of the FEM procedure, the computational effort was greatly increased, wherein the accuracy of the simulation was insignificant compared to the computation time used. In regard to the drawn model and solver parameters' selection, the simulation lasted approximately 4 min on a Windows 10-based PC, with a i7-3770@3.4GHz CPU and 8GB-RAM memory. The computational mesh in the selected computational boundary surface and calculated quantities with FEM are depicted in Figure 4. In Figure 3. with the selection of solver precision (1e −8 ) and minimal angle value setting (30) in a millimeter scale for the inside computational boundary surface, we get a computational mesh with 17, 257 nodes. The selection of the solver parameters (solver precision and mesh angle) is crucial in the accuracy of the simulation results. In our approach the solver parameters were selected experimentally in a manner to acquire accurate and reliable results in a reasonable time. By increasing the solver precision and angle value of the FEM procedure, the computational effort was greatly increased, wherein the accuracy of the simulation was insignificant compared to the computation time used. In regard to the drawn model and solver parameters' selection, the simulation lasted approximately 4 min on a Windows 10-based PC, with a i7-3770@3.4GHz CPU and 8GB-RAM memory. The computational mesh in the selected computational boundary surface and calculated quantities with FEM are depicted in Figure 4.
The FEMM software offers calculation of many different quantities and values of the designed electromechanical system. The electromagnetic force F e for the given system in Figure 4 was calculated via a weighted stress tensor, where the force value components were mapped in a 2D plane.  The FEMM software offers calculation of many different quantities and values of the designed electromechanical system. The electromagnetic force for the given system in Figure 4 was calculated via a weighted stress tensor, where the force value components were mapped in a 2D plane.

Coupling FEMM with Matlab and Data Analysis
The coupling procedure is attained after designing a model in FEMM software with all the material characteristics. There were many choices for how to achieve and how to handle the objects inside the FEMM from outer programmes or self-defined scripts. As we mentioned before, the Matlab software was used to handle model objects externally in FEMM. Matlab-FEMM coupling has many possibilities for managing model objects and FEMM properties. The whole design of the model and simulation procedure inside the Matlab software was presented by Benamimour [36]. In our case, we used the possibility of simulation commissioning and model objects handling in an already existing FEMM model. The coupled Matlab software was used to control and acquire simulated data from FEMM. Matlab software scheduled the simulation parameters, where the driven coil current and vertical position of the PM were adjusted. The whole experiment was based on the measurement of the electromagnetic force, which was acquired from the FEMM at the exact specified value of the current and PM position. The current value was changed on interval [0 A ÷ 1.4 A] by steps of 0.2 A, where position was changed on interval [15 mm ÷ 45 mm] by steps of 0.5 mm. The simulation schedule is depicted in Figure 5.
The simulated values are further used for model parameters estimation given in Equation (2). The obtained results from the simulation are presented in Figures 6 and 7.

Coupling FEMM with Matlab and Data Analysis
The coupling procedure is attained after designing a model in FEMM software with all the material characteristics. There were many choices for how to achieve and how to handle the objects inside the FEMM from outer programmes or self-defined scripts. As we mentioned before, the Matlab software was used to handle model objects externally in FEMM. Matlab-FEMM coupling has many possibilities for managing model objects and FEMM properties. The whole design of the model and simulation procedure inside the Matlab software was presented by Benamimour [36]. In our case, we used the possibility of simulation commissioning and model objects handling in an already existing FEMM model. The coupled Matlab software was used to control and acquire simulated data from FEMM. Matlab software scheduled the simulation parameters, where the driven coil current and vertical position of the PM were adjusted. The whole experiment was based on the measurement of the electromagnetic force, which was acquired from the FEMM at the exact specified value of the current and PM position. The current value was changed on interval [0 A ÷ 1.4 A] by steps of 0.2 A, where position was changed on interval [15 mm ÷ 45 mm] by steps of 0.5 mm. The simulation schedule is depicted in Figure 5.
The simulated values are further used for model parameters estimation given in Equation (2). The obtained results from the simulation are presented in Figures 6 and 7.

Modelling of Electromagnetic Force F e
The modelling of the electromagnetic force will be presented from the obtained simulation results from Matlab and FEMM software. The modelling of the electromagnetic F e has been based on the function fitting procedure to the set of the simulated result data's with preselected model structure. Model estimation and data fitting was performed with the quadratic programming (QP) method. The QP technique is a well-known and efficient approach for convex optimization with introduced, boundaries, inequalities and equalities [37]. Before the objective function for optimization will be derived, the same optimization relaxation will be assigned. The modelling of the F e will be considered only on the system with fixed structure and geometrical form. The current and distance dependence by given fixed preselected system will be considered. From the simulated data presented in Figure 7 it is straightforward to recognise that the current dependence is completely linear to the force F e . This assumption can be proven with the calculated difference between each adjacent current's characteristic in Figure 7. Each current characteristic has been labelled as I k , where k represents corresponding current value from 0 A to 1.4 A by steps of 0.2 A. The linearity test L p of current to the F e is given with: where p represents the number of pairs of adjacent difference characteristics. In the presented case with 8 current characteristics k ∈ [0 ÷ 8] in Figure 7, we have seven difference pairs; p = 7. With deviation and mean calculation of the L p we have proved the linear dependence of the current to the force F e .
The calculated values are presented in Table 2.
The distance dependence needs to be determined in regard to the linear correlation of current and force. The distance and force dependence was estimated with a curve fitting algorithm and the QP optimization technique. The QP optimization technique allows finding an optimal solution in the given interval of search parameters with selected equality and inequality constraints. In the given case, only positive solutions are allowed. The optimization procedure was divided into two stages. In the first stage the general model structure of the model was specified, where the polynomial parameters stay uncertain. The corresponding preselected model was fitted for each current characteristic in Figure 7. After the first stage of the optimization, we get a set of different models with fixed structure and different polynomial coefficients. In regard to Figure 7 and the simulation schedule we derived eight different models. The second stage was intended to estimate polynomial coefficients, where each coefficient in the polynomial has eight values from a fitted model. Each coefficient is also approximated with a new coefficient polynomial function. The coefficient polynomial functions are inserted further in the preselected model structure from the previous stage one.

Model Structure Selection F e -Stage One
The first stage begins with selection of the simple model structure. The best two model candidates were selected after a few iterations and with the assumption of the current-force linearity. The best model candidates for F e are: where coefficients a, b, c are estimated with QP optimization, for each separated current characteristic. The variable d is PM proximity. Selected residual functions for QP optimization are: where n is a number of simulation data from FEMM. The parameters a, b, c are solution of the QP programming technique with residual Equations (7) and (8). The QP solutions with corresponding matrixes are: Results of model fitting F 1e and F 2e over current characteristic from 0 A to 1.4 A, are presented in Table 3 and Figures 8 and 9. Table 3. Parameter estimation of the model F 1e and F 2e , by current characteristics 0 A-1.4 A.

Current
Model    The model fitting on current characteristics is presented in Figure 9, where Figures A and B show the fitting properties with models and , respectively. Only three current characteristics 0 A, 0.6 A and 1.2 A are presented for better clarity of the results. It is obvious from Table 3, Figures 8 and  9, that model has better data matching to the current characteristics than model . From Table 3 it can be seen that the residual values indicate that model structure has quite accurate fitting properties over the whole area, where residual values remain constant with neglect deviation in regard to the residual values of the model . The model will be considered for further analysis.  The model fitting on current characteristics is presented in Figure 9, where Figure 9a,b show the fitting properties with models F 1e and F 2e , respectively. Only three current characteristics 0 A, 0.6 A and 1.2 A are presented for better clarity of the results. It is obvious from Table 3, Figures 8 and 9, that model F 1e has better data matching to the current characteristics than model F 2e . From Table 3 it can be seen that the residual values indicate that model structure F 1e has quite accurate fitting properties over the whole area, where residual values remain constant with neglect deviation in regard to the residual values of the model F 2e . The model F 1e will be considered for further analysis.    The model fitting on current characteristics is presented in Figure 9, where Figures A and B show the fitting properties with models and , respectively. Only three current characteristics 0 A, 0.6 A and 1.2 A are presented for better clarity of the results. It is obvious from Table 3, Figures 8 and  9, that model has better data matching to the current characteristics than model . From Table 3 it can be seen that the residual values indicate that model structure has quite accurate fitting properties over the whole area, where residual values remain constant with neglect deviation in regard to the residual values of the model . The model will be considered for further analysis.

The F 1e Coefficient Estimation and F e Model Derivation-Stage Two
After model structure derivation F 1e the next step of the optimization procedure is estimation of coefficient functions. The coefficient function describes the coefficient change in model F 1e , where it is apparent from Table 3 that, the estimated model F 1e has variable coefficients. Figure 10 shows the parameter change of model F 1e with regard to the current characteristics.

The Coefficient Estimation and Model Derivation-Stage Two
After model structure derivation the next step of the optimization procedure is estimation of coefficient functions. The coefficient function describes the coefficient change in model , where it is apparent from Table 3 that, the estimated model has variable coefficients. Figure 10 shows the parameter change of model with regard to the current characteristics. With regard to Figure 10 only parameter needs to be estimated, where parameter has a constant value with neglected deviation. The mean value of parameter is 5.996 × 10 . Parameter can be approximate with linear function. Such function indicates linear dependence between force and current, which was confirmed before with Table 2. The selected linear function of parameter is: where ( ) is a parameter function of current and unknown parameter . The same optimization procedure was used as for model Equations (5) and (6), where residual function is: The obtained results are presented in Table 4 and in Figure 10. Fitting solutions of the parameters and are presented in Figure 10. After derivation of the parameter , function ( ) and model structure in Equation (5), we get an accurate model of the electromagnetic force . The electromagnetic force model is:     (13), consideration of objective function given in Equation (7) and simulated FEMM results. With regard to Figure 10 only parameter b needs to be estimated, where parameter a has a constant value with neglected deviation. The mean value of parameter a is 5.996 × 10 −9 . Parameter b can be approximate with linear function. Such function indicates linear dependence between force and current, which was confirmed before with Table 2. The selected linear function of parameter b is: where f b (i) is a b parameter function of current and unknown parameter m. The same optimization procedure was used as for model Equations (5) and (6), where residual function is: The obtained results are presented in Table 4 and in Figure 10. Fitting solutions of the parameters a and b are presented in Figure 10. After derivation of the parameter a, function f b (i) and model structure in Equation (5), we get an accurate model of the electromagnetic force F e . The electromagnetic force model F e is:  (13), consideration of objective function given in Equation (7) and simulated FEMM results. The graphical model fitting results with derived model F e given in Equation (13) and FEMM was already presented in Figure 9a.

Distance Measurement and Sensor Fusion Algorithm
After derivation of a proper mathematical model, the next step was to design an algorithm for accurate position measurement. There are many approaches and estimation algorithms which offer many promising and useful results. The basic problem of distance measurement with a Hall sensor is the noisy and biasing output voltage, which is caused by structural and external factors. In our application, the disturbances on the Hall sensor can be generally defined as: external magnetic field influence, from the driven coil, the temperature dependence and structural imperfections of the Hall element. All influences have a complex nonlinear dependence, where the magnetic field from a driven coil can be prior estimated with static characteristics. With determined static characteristic the disturbances can be mainly supressed but the accuracy of the measured distance can still remain incorrect due to the other unknown disturbances, voltage drift, sensor noise, imperfect static characteristics etc. The measured static characteristics for the system in Figure 1a with an Allegro MicroSystems A1321 Hall sensor (Allegro MicroSystems, LLC, Worcester, MA, USA) are depicted in Figure 11. Figure 11a presents the slack linear static characteristic between the coil and Hall voltage, where the Hall voltage and PM distance have a nonlinear relation. Both characteristics in Figure 11 are used for direct proximity measurement with a Hall sensor; conversion from a Hall voltage to PM distance. The characteristics Hall-coil voltage and Hall voltage and PM position are estimated in the same fashion as the characteristic for electrical force F e given in Equations (7) and (8).
The equations for direct distance measurement with static characteristics compensation are: where CH (v), MH (v), v, v hall , d hall represents the Hall-coil voltage characteristic, magnet-Hall characteristic, coil voltage, Hall voltage and PM proximity, respectively. Direct distance measurement is composed with coil voltage compensation CH (v), after which the PM proximity is obtained from the MH (v) characteristic.

Unscented Kalman Filter and Sensor Fusion Algorithm
The Kalman filter-KF is a broadly used estimation, prediction and sensor fusion algorithm [19,38,39]. There exist many variants of the KF algorithm, where the Extended Kalman Filter-EKF and Unscented Kalman Filter-UKF are particularly used to deal with nonlinear systems. The central operation of the linear KF is the propagation of the Gaussian random variable through the system dynamics, where covariance of the estimation error needs to be minimized. The EKF used the same procedure as the linear KF, where the Gaussian random variable is approximated analytically with Jacobian or Hessian matrices (Taylor series-TS approximation-linearization). Such approximations can, in some cases, introduce large errors in the true posterior mean and covariance of the Gaussian random variable, which may lead to poor performance of the filter [40]. The UKF addresses the issues with approximation of the Gaussian random variables with the Taylor series. Similar to the TS approximation the Unscented Transformation (UT) can be used for forming a Gaussian approximation of the joint distribution of random variables, which are defined in the nonlinear dynamic system. UT transformation deterministically chooses a relatively small amount of the fixed number of sigma points, which capture the true mean and covariance of the Gaussian random state variable [41]. The UT transformation is a method for calculation of the statistics of a random variable which undergoes a nonlinear system. In regard to the TS approximation in the EKF algorithm, the UT transformation is better at capturing the higher order of moments caused by the non-linear transform and is less error prone in regard to the calculation of the Jacobian and Hessian matrices-TS [41][42][43]. The basic framework of the UKF involves the estimation of the state of the discrete time nonlinear system. The discrete time non-linear system is: where represents the unobservable states of the system, is a system input and is the measured output signal. The and are process and measurement noise, respectively, with corresponding process noise and measurement noise covariance matrices. The non-linear state and output functions are known. The UKF algorithm is just a straightforward extension of the UT transformation to the recursive estimation of standard state update equation: where is Kalman gain and is estimated state vector-filter output. The UT transformation sigma points are applied in the new augmented sigma point matrix , which is generalized obtained from a value of the vector and the state covariance matrix [41]. The UKF algorithm equations and calculation procedure are given below.

Unscented Kalman Filter and Sensor Fusion Algorithm
The Kalman filter-KF is a broadly used estimation, prediction and sensor fusion algorithm [19,38,39]. There exist many variants of the KF algorithm, where the Extended Kalman Filter-EKF and Unscented Kalman Filter-UKF are particularly used to deal with nonlinear systems. The central operation of the linear KF is the propagation of the Gaussian random variable through the system dynamics, where covariance of the estimation error needs to be minimized. The EKF used the same procedure as the linear KF, where the Gaussian random variable is approximated analytically with Jacobian or Hessian matrices (Taylor series-TS approximation-linearization). Such approximations can, in some cases, introduce large errors in the true posterior mean and covariance of the Gaussian random variable, which may lead to poor performance of the filter [40]. The UKF addresses the issues with approximation of the Gaussian random variables with the Taylor series. Similar to the TS approximation the Unscented Transformation (UT) can be used for forming a Gaussian approximation of the joint distribution of random variables, which are defined in the nonlinear dynamic system. UT transformation deterministically chooses a relatively small amount of the fixed number of sigma points, which capture the true mean and covariance of the Gaussian random state variable [41]. The UT transformation is a method for calculation of the statistics of a random variable which undergoes a nonlinear system. In regard to the TS approximation in the EKF algorithm, the UT transformation is better at capturing the higher order of moments caused by the non-linear transform and is less error prone in regard to the calculation of the Jacobian and Hessian matrices-TS [41][42][43]. The basic framework of the UKF involves the estimation of the state of the discrete time nonlinear system. The discrete time non-linear system is: where x k represents the unobservable states of the system, u k is a system input and y k is the measured output signal. The w k and n k are process and measurement noise, respectively, with corresponding process noise Q and measurement noise R covariance matrices. The non-linear state F and output H functions are known. The UKF algorithm is just a straightforward extension of the UT transformation to the recursive estimation of standard state update equation: where κ k is Kalman gain andx k is estimated state vector-filter output. The UT transformation sigma points are applied in the new augmented sigma point matrix χ a k−1 , which is generalized obtained from a value of the vectorx k−1 and the state covariance matrix P k−1 [41]. The UKF algorithm equations and calculation procedure are given below.

•
Initialization of the UKF-filter:x Variablex 0 represents initial states of nonlinear system given in Equation (17) and P 0 is an initial covariance matrix of the state variable x 0 . • Prediction: UT-transformation, calculation of 2L + 1 sigma points: and associated weights: The variable L is the number of system states and λ is the scaling parameter; λ = α 2 (L + k i ) − L.
Parameter α determinates the spreads of the sigma points, k i is a secondary scaling parameter and β is used to incorporate prior knowledge of the distribution χ a k . The weights W m and W c represent mean weighting factor and estimation error covariance weighting factor respectively [42].
• Measurement Update Equation x k|k =x k|k−1 + κ k y k −ŷ k|k−1 , where variables χ a k|k−1 , Υ Υ Υ k|k−1 ,x k|k−1 ,ŷ k|k−1 , Q k , R k , P k|k−1 , P y k ,y k , P x k ,y k , κ k ,x k|k , P k|k are prior sigma points states, prior output sigma points obtained from a nonlinear model given in Equation (17), prior estimated state vector, prior estimated output vector, process noise covariance matrix, measured noise covariance matrix, prior state covariance matrix, prior output covariance matrix, cross covariance matrix, Kalman gain, posterior state vector-UKF output and a posteriori state covariance matrix, respectively. The UKF algorithm starts with the a priori selected weights W m , W c and scaling parameters λ, k i . The UKF filter is used for the sensor fusion procedure in a manner to improve accuracy of the distance measurement with the Hall sensor. The direct measured value from the Hall sensor with Equation (14) represented variable y k , whereŷ k is the estimated variable from the derived nonlinear model of Equations (2), (3) and (13).

Deployment of the UKF Filter
The UKF filter, as a sensor fusion algorithm, is based on two main baselines; the real direct proximity measurement from the Hall sensor with nonlinear characteristic compensation given in Equation (14) and the proximity estimation from the selected non-linear dynamic model in Equations (2), (3) and (13). Both the proximity values are used in the UKF algorithm. The state space representation of the derived nonlinear model is: .
is a system state vector x (t) R 3 and represents; d-PM vertical position, v-PM velocity and i-coil current. The size of the state vector is L = 3, which means that the UKF operated with 2L + 1 sigma points. For the further use on a real-time system, the continuous-time Equation (32) were discretized using the Euler integration scheme. The discrete form of the continuous system in Equation (32) is: where t s is a sampling time with preselected value of, t s = 0.5 ms. Model system output is the predicted proximity d (k). For better transparency of the paper text, we assumed that the d (k) = d k and it applies for all other variables v (k) , u (k) , etc. The state variablex k of the UKF algorithm in regard to the Equations (22)-(31) is equal tox k = [d k ,v k ,î k ] T , where input vector u k is a driven coil voltage. An important part of the algorithm is the state estimation update in Equation (30), which is based on error calculation between the direct proximity measurement d hall k obtained from Equation (14) and UKF's output stated k|k−1 given in Equations (23) and (25). The state update equation is,x k|k =x k|k−1 + κ k (d hall k −d k|k−1 ). The derived model in Equation (33) was used for sigma state estimation χ a k|k−1 in Equation (22) and sigma output estimation Υ Υ Υ k|k−1 in Equation (25). The discrete process noise Q k and measured noise R k were set to: Other fixed parameters of the UKF have been selected arbitrary with values: α = 0.002, k i = 0, β = 2 [40]. The selected weighting matrices W m and W c are:   Figure 12 represents deploying UKF as a sensor fusion algorithm for accurate proximity measurement from a Hall sensor.

Experimental Results
Testing the sensor fusion algorithm with a UKF filter was done predominantly in a closed loop system, with nonlinear and linear controllers. The nonlinear controller was obtained on the basis of backstepping controller design. The backstepping method is a recursive design technique, which stabilises the origin of the system in strict feedback form [44,45]. The nonlinear controller was synthesised from the derived model in Equation (32). A simplified second order nonlinear model in Equation (32) was used for the simplification of the Backstepping controller design. In the design procedure we assumed that the dynamic of the electrical coil was much higher than the mechanical part, where we got

Experimental Results
Testing the sensor fusion algorithm with a UKF filter was done predominantly in a closed loop system, with nonlinear and linear controllers. The nonlinear controller was obtained on the basis of backstepping controller design. The backstepping method is a recursive design technique, which stabilises the origin of the system in strict feedback form [44,45]. The nonlinear controller was synthesised from the derived model in Equation (32). A simplified second order nonlinear model in Equation (32) was used for the simplification of the Backstepping controller design. In the design procedure we assumed that the dynamic of the electrical coil was much higher than the mechanical part, where we got di dt = 0. The static current value was (t) = 61.61 165.6 u (t). The static current value i (t) was placed to the velocity equation dv dt of the model in Equation (32), wherein the second order model was obtained. For better reference tracking capability of the feedback system, the supplementary error state was introduced; e (t) = re f (t) − d (t). The re f (t) is a reference tracking value. The Backstepping controller procedure has had two 'back steps', from the proximity-error sate e (t) over v (t) to the driven voltage as a system input and controller output. The nonlinear controller structure was: where u k , re f k , re f d k , re f dd k are controller output-coil voltage, reference value, first derivative and second derivatives of the reference value, respectively. For the first and second derivatives a cut-off filter was used with frequency at the 150 Hz and 125 Hz. The velocity of the suspended object v k was not measured directly with separated velocity sensor, but was estimated from a sensor fusion algorithm with UKF. For the second test a classic linear PID structure with output clamping algorithm was used to validate the system efficiency with UKF. The PID controller was tuned with the linear approximation of the model given in Equation (32) and Integra of Square Error (ISE) performance index. The controller was designed with Control System Design and Optimization toolboxes in a MATLAB 2016a environment. The clamping algorithm was used to prevent proper operation of the system and prevent integrator windup. The PID controller structure with clamping output was: where u k , e k are controller output and current controller error, respectively. The sensor UKF fusion algorithm and nonlinear controller are running on an ARM Cortex-M4 STM32F407VGT6 microcontroller (STMicroelectronics, Geneva, Switzerland) with floating point unit-FPU and 12 bit AD conversion for accurate Hall voltage measurement. Figure 13 represents a real-time system configuration with a Hall proximity sensor and ARM microcontroller. The schematic in Figure 14 represents the structure of the sensor fusion algorithm with the feedback controller in Figure 14a and the algorithm procedure flow chart in Figure 14b. algorithm with UKF.
For the second test a classic linear PID structure with output clamping algorithm was used to validate the system efficiency with UKF. The PID controller was tuned with the linear approximation of the model given in Equation (32) and Integra of Square Error (ISE) performance index. The controller was designed with Control System Design and Optimization toolboxes in a MATLAB 2016a environment. The clamping algorithm was used to prevent proper operation of the system and prevent integrator windup. The PID controller structure with clamping output was: where , are controller output and current controller error, respectively. The sensor UKF fusion algorithm and nonlinear controller are running on an ARM Cortex-M4 STM32F407VGT6 microcontroller (STMicroelectronics, Geneva, Switzerland) with floating point unit-FPU and 12 bit AD conversion for accurate Hall voltage measurement. Figure 13 represents a real-time system configuration with a Hall proximity sensor and ARM microcontroller. The schematic in Figure 14 represents the structure of the sensor fusion algorithm with the feedback controller in Figure 14a and the algorithm procedure flow chart in Figure 14b.    (14)). The reference value of the feedback system represents the true proximity of the PM obtained from the external independent ruler. The accuracy of the external ruler was around 0.2 mm. The experiments were tested with   (14)). The reference value of the feedback system represents the true proximity of the PM obtained from the external independent ruler. The accuracy of the external ruler was around 0.2 mm. The experiments were tested with feedback controllers (backstepping and PID), where feedback signals (position, velocity) are taken from the UKF sensor fusion algorithm. in the presented results, the reference value was changed in the span of 22 mm to the 30 mm from the driven coil.   (14)). The reference value of the feedback system represents the true proximity of the PM obtained from the external independent ruler. The accuracy of the external ruler was around 0.2 mm. The experiments were tested with feedback controllers (backstepping and PID), where feedback signals (position, velocity) are taken from the UKF sensor fusion algorithm. in the presented results, the reference value was changed in the span of 22 mm to the 30 mm from the driven coil.           where d k is a true value andd k is the measurement from direct approximation or sensor fusion algorithm with UKF. Table 6. RMS values of the signals.

RMS Backstepping Controller PID Controller
Direct measurement 12.32 × 10 −5 15.71 × 10 −5 Sensor fusion with UKF 1.52 × 10 −5 1.68 × 10 −5 Figure 19 presents a comparison between the UKF and EKF sensor fusion algorithms with the backstepping controller given in Equation (34). The EKF algorithm uses a linear approximation (Jacobian matrix) of the model in Equation (32). For better comparison of both algorithms the span of reference value was selected between 22 mm and 25 mm.   Figure 19 presents a comparison between the UKF and EKF sensor fusion algorithms with the backstepping controller given in Equation (34). The EKF algorithm uses a linear approximation (Jacobian matrix) of the model in Equation (32). For better comparison of both algorithms the span of reference value was selected between 22 mm and 25 mm. The results presented in Figures 15-18 show the effectiveness and reliability of the sensor fusion algorithm with UKF. The UKF, in regard to the direct measurement application improves proximity measurements significantly, suppresses Hall drift and lowers noise influence on the proximity information signal. The effectiveness of the noise suppression is confirmed with the Frequency spectrum plot in Figures 17 and 18 and with Table 6, where the RMS values of the signals are calculated. From Figure 19 it can be seen clearly that the UKF sensor fusion algorithm outperforms the EKF. The advantage of UKF can be noticed also in its reference tracking and noise suppression capability. The resolution of the measurement in regard to the feedback stability region between 18 mm and 32 mm is estimated at 0.05 mm with regard to Equations (14)-(16), Hall sensitivity value and 12 bit AD resolution of the ARM microcontroller. The accuracy of the proximity measurement in regard to the external ruler is estimated at 0.2 mm. At the end of the experiment it needs to be The results presented in Figures 15-18 show the effectiveness and reliability of the sensor fusion algorithm with UKF. The UKF, in regard to the direct measurement application improves proximity measurements significantly, suppresses Hall drift and lowers noise influence on the proximity information signal. The effectiveness of the noise suppression is confirmed with the Frequency spectrum plot in Figures 17 and 18 and with Table 6, where the RMS values of the signals are calculated.
From Figure 19 it can be seen clearly that the UKF sensor fusion algorithm outperforms the EKF. The advantage of UKF can be noticed also in its reference tracking and noise suppression capability. The resolution of the measurement in regard to the feedback stability region between 18 mm and 32 mm is estimated at 0.05 mm with regard to Equations (14)- (16), Hall sensitivity value and 12 bit AD resolution of the ARM microcontroller. The accuracy of the proximity measurement in regard to the external ruler is estimated at 0.2 mm. At the end of the experiment it needs to be mentioned that the close-loop system in both configurations with backstepping and PID controllers are unstable if direct measurement from the Hall sensors is used in controller operation.

Conclusions
The main contribution of the presented paper is accurate distance measurement with low cost sensing devices in the presence of a magnetic field. In the present case a low cost linear Hall sensor was used. The rough data from the sensor are relatively noisy and contain exogenous disturbance effects, which need to be removed or suppressed efficiently. The paper findings show great distance measurement results with regard to the system open loop instability and feedback controller sensitivity to the sensor noise and signal drift. The sensor fusion algorithm with UKF improved the accuracy of distance measurement and system states estimation drastically, which are also used in feedback control. The efficiency of the sensor fusion algorithm originated from the accuracy of the system model derived with the Finite Element Method. Both approaches, modelling of the system dynamics with the Finite Element Method and UKF filter settings, are crucial in the proximity measuring with a Hall sensor. Such an approach can be used in many different electro-mechanical applications where a relatively uncertain sensor is used and system behavior is known. The approach offers a great potential to acquiring the quantities, which are not directly measured with separated sensors but are estimated with the model and sensor fusion algorithm (measuring the coil current and the velocity of the levitating object).