Design of a Control System for a Maglev Planar Motor Based on Two-Dimension Linear Interpolation

Abstract: In order to realize the high speed and high-precision control of a maglev planar motor, a high-precision electromagnetic model is needed in the first place, which can also contribute to meeting the real-time running requirements. Traditionally, the electromagnetic model is based on analytical calculations. However, this neglects the model simplification and the manufacturing errors, which may bring certain errors to the model. Aiming to handle this inaccuracy, this paper proposes a novel design method for a maglev planar motor control system based on two-dimensional linear interpolation. First, the magnetic field is divided into several regions according to the symmetry of the Halbach magnetic array, and the uniform grid method is adopted to partition one of these regions. Second, targeting this region, it is possible to sample the electromagnetic forces and torques on each node of the grid and obtain the complete electromagnetic model in this region through the two-dimensional linear interpolation method. Third, the whole electromagnetic model of the maglev planar motor can be derived according to the symmetry of the magnetic field. Finally, the decoupling method and controller are designed according to this electromagnetic model, and thereafter, the control model can be established. The designed control system is demonstrated through simulations and experiments to feature better accuracy and meet the requirements of real-time control.


Introduction
Planar motion is required for various applications in the high-precision industry, such as pick and place machines, semiconductor lithography, and transport devices [1,2].Planar movement can be realized through super-accurate 2-D positioning stages.Traditionally, these stages are often constructed by stacking several one-degree-of-freedom linear drives supported by mechanical bearings [3].However, this type of positioning stage has a high moving mass and low stiffness, which restricts its dynamic performance.In recent years, the 2-D positioning stage has been driven directly by the maglev planar motor [4][5][6][7], which makes it an integrated structure featuring many advantages, such as direct driving, fast response, high sensitivity, good dynamic nature and simple structure [8,9].
A maglev planar motor is a multi-input and multi-output, non-linear and strongly coupled system.The inputs of a maglev planar motor are the currents imposed on each winding, and the outputs are the moving forces on the six degrees of freedom.Solving the corresponding current according to the required forces on the different degrees of freedom is a critical step for achieving six degrees movement of this motor.To realize the position control on the six degrees of freedom of this motor, decoupling the currents of multiple sets of windings and the forces and torques on six degrees of freedom is essential.To obtain this decoupling relation, an accurate electromagnetic model of a maglev planar motor is Energies 2017, 10, 1132 2 of 17 needed.Meanwhile, this electromagnetic model is supposed to have fast running speed to meet the requirements of high dynamic response of the control system.Therefore, on the basis of meeting the requirements of real-time running speed, a more precise electromagnetic model is indispensable for a high speed and high precision maglev planar motor control system.
It has been concluded in [10] that the surface charge model and harmonic model have high model accuracy, but they are not suitable for practical implementation because the computation time needed is too long, often within the scale of seconds.On the other hand, analytical models have relatively lower model accuracy (still deemed high enough for industrial applications) but the computation time needed (normally less than 2 µs) on a typical DSP control system is very small.Therefore, the analytical method has been widely used in the planar motor control implementation.
Hamers points out three analytical ways to solve magnetic force in his research report, namely the virtual work method, Maxwell stress method, and the Lorentz force equation method [11].Normally, the Lorentz force equation method can compute the electromagnetic force and torque of the coils in a maglev planar motor, which is the most convenient method and widely adopted in establishing electromagnetic models of the maglev planar motor [12][13][14].However, this method can not be directly used to calculate the electromagnetic force on the corner segment of the coil, therefore, the shape of the real coil needs to be simplified.In [15], the coils are simplified into four rectangles, but there are overlapping parts during the calculation process.In [16], the winding is regarded as four trapezoids, which reduces the overlapping parts.However, this simplified model still takes a right angle as a real circle corner.Reference [17] proposes a method that converts the electromagnetic force and torque modeling of the corner coils to the polar coordination system.This method can improve the accuracy of the electromagnetic force and torque of the corner coils.However, in all the methods above, in order to achieve the real-time computation, the Lorentz force equation method only takes the fundamental component of magnetic flux density into account.As a result, the ignored harmonic components will cause relatively significant errors in the electromagnetic model.
An improved semi-analytical method based on a single-point magnetostatic field simulation is proposed in [18] to get a more accurate electromagnetic model.This method does not assume the relative permeability to be 1 as done in many traditional models.Moreover, it does not neglect the difference between the attractive force and the repulsive force.However, the needed computation time is too long to make it suitable for real-time implications.In [19], numerical calculations and a measurement setup are used to model the force and torque generated by each coil.A novel single-deck six degree of freedom magnetic levitation stage has been introduced in [20] which also analyzes semi-empirical equations of the motor and furthermore builds a model based on the ratio of the measured force to the current.The forces in the direction of x/y/z are generated by different coil arrays, thus only one-dimensional fitting has been adopted due to the lack of electromagnetic coupling.
Interpolation is a mathematical method to get a continuous function based on discrete data points, and makes all the original data points located on the new continuous curve.The interpolation is an important method for the approximation of discrete functions, and it can estimate the approximate value of other points based on the original data points.Interpolation can get the complete model based on limited measured data points; moreover, the dimensions can be easily extended if needed.Therefore, it has been widely used in numerous application fields [21][22][23][24].
In order to minimize the errors arising from the model simplification and manufacturing process and to meet the requirements of real-time control simultaneously, this paper proposes a new method to obtain accurate electromagnetic force and torque based on the two-dimensional interpolation method.Based on the above method, the decoupling computation of the maglev planar motor can be realized and a high performance real-time control model can be designed.

Structure and Working Principle
Figure 1 shows the moving-coil ironless maglev permanent magnet planar motor.It consists of a mover, a stator, and the gap between them.The mover consists of the base plate and the 4×4 coil array, and the stator consists of the yoke plate and the permanent magnet array.The stator permanent magnet is arranged in a Halbach structure, which can increase the flux density near the coils.The coils on the mover are connected to the outer power supply.More specifically, each coil is supplied by one independent current amplifier, thus the decoupled current can be realized through the connected independent power supply.

Structure and Working Principle
Figure 1 shows the moving-coil ironless maglev permanent magnet planar motor.It consists of a mover, a stator, and the gap between them.The mover consists of the base plate and the 4×4 coil array, and the stator consists of the yoke plate and the permanent magnet array.The stator permanent magnet is arranged in a Halbach structure, which can increase the flux density near the coils.The coils on the mover are connected to the outer power supply.More specifically, each coil is supplied by one independent current amplifier, thus the decoupled current can be realized through the connected independent power supply.The working principle of the proposed planar motor is described as follows.The permanent magnet generates a magnetic field, in which the coils with current in them will generate electromagnetic force and electromagnetic torque.Since the magnets array remains still, the electromagnetic force and torque on the whole coils drive the mover move.As long as the proper current according to the relative position between the coils and the magnetic field are given, the control of electromagnetic force and torque on six degrees of freedom can be realized and the position control of the mover on six degrees of freedom can be achieved.

Coordinate Definition and Transformation
In Figure 2, the coil is divided into four straight segments.τn denotes the pole pitch, lc denotes the coil length of the straight coil segment, hc is the thickness of the coil, wc is the width of the coil, J is current density.Two different coordinate systems are defined in Figure 2.  The working principle of the proposed planar motor is described as follows.The permanent magnet generates a magnetic field, in which the coils with current in them will generate electromagnetic force and electromagnetic torque.Since the magnets array remains still, the electromagnetic force and torque on the whole coils drive the mover move.As long as the proper current according to the relative position between the coils and the magnetic field are given, the control of electromagnetic force and torque on six degrees of freedom can be realized and the position control of the mover on six degrees of freedom can be achieved.

Coordinate Definition and Transformation
In Figure 2, the coil is divided into four straight segments.τ n denotes the pole pitch, l c denotes the coil length of the straight coil segment, h c is the thickness of the coil, w c is the width of the coil, J is current density.Two different coordinate systems are defined in Figure 2.

Structure and Working Principle
Figure 1 shows the moving-coil ironless maglev permanent magnet planar motor.It consists of a mover, a stator, and the gap between them.The mover consists of the base plate and the 4×4 coil array, and the stator consists of the yoke plate and the permanent magnet array.The stator permanent magnet is arranged in a Halbach structure, which can increase the flux density near the coils.The coils on the mover are connected to the outer power supply.More specifically, each coil is supplied by one independent current amplifier, thus the decoupled current can be realized through the connected independent power supply.The working principle of the proposed planar motor is described as follows.The permanent magnet generates a magnetic field, in which the coils with current in them will generate electromagnetic force and electromagnetic torque.Since the magnets array remains still, the electromagnetic force and torque on the whole coils drive the mover move.As long as the proper current according to the relative position between the coils and the magnetic field are given, the control of electromagnetic force and torque on six degrees of freedom can be realized and the position control of the mover on six degrees of freedom can be achieved.

Coordinate Definition and Transformation
In Figure 2, the coil is divided into four straight segments.τn denotes the pole pitch, lc denotes the coil length of the straight coil segment, hc is the thickness of the coil, wc is the width of the coil, J is current density.Two different coordinate systems are defined in Figure 2.  The first one is the stator global coordinate system denoted with the superscript m, which is fixed on the top face of the stator permanent magnet.It can be expressed as: The second one is the mover local coordinate system denoted with the superscript c, which is fixed on the mover bottom surface center.It can be expressed as: The coordinate position of the k-th coil center in the mover local coordinate system can be expressed as: The position vector transformation between the mover local coordinate system and the stator global coordinate system is expressed as: The vector in the mover local coordinate system and that in the stator global coordinate system can be transformed into each other through the following [25,26]: The orientation transformation matrix is defined as: where, ϕ, θ and φ are the rotation angles about the c x-, c y-, and c z-axes, respectively.

Electromagnetic Force and Torque
According to the Lorenz's law, the electromagnetic force and torque on the coil can be derived as: where, B is the flux density distribution of the magnet array; J is the current density in the coil; V is the volume of the coil; r is the vector from the point about which the torque is calculated.
In the mover local coordinate system, the electromagnetic force and torque on the straight segment of the coil can be derived as: where, c B xy is the flux density distribution component on x-y plane.C z can be expressed as: Based on Equations ( 14) and ( 15), the electromagnetic force and torque on the k-th coil can be derived as: where, C 1 and C 2 are the coefficients related to the motor structure; C is harmonic compensation factor; c r z is the equivalent force arm in the c z-direction from the coil to the mover mass center.
Once the coil is put on the mover, the position deviation on the x-y plane can be seen as constant, therefore, the whole electromagnetic force and torque on the mover is the sum result of that on each winding, taking the position deviation into account.

Electromagnetic Model Based on Two-Dimensional Linear Interpolation
Analytical modeling seeks mathematical functions and equations that are obtained from the closed-form (exact or approximate) solution to the original physics-governing equations.The outcome of analytical modeling is behavioral models [27].For instance, Equations ( 17)-( 22) provide the analytical modeling for the electromagnetic force/torque of the maglev planar motor.The interpolation method, on the other hand, processes the discrete data with the numerical approach based on the known data, thus a complete mathematical model of the control target object can be derived.Normally, interpolation includes linear interpolation, polynomial interpolation, cubic spline interpolation, which extend from single dimension to multiple dimensions [28].To meet the real-time computation needs, Energies 2017, 10, 1132 6 of 17 both the accuracy and speed should be taken into account to choose a proper interpolation method; therefore, the two-dimension interpolation method has been chosen.

Interpolation Region Partition
Since the magnetic field of a Halbach array is strictly periodic symmetrical, it can be simplified to take points in a small range to perform the interpolation.The whole planar model can be derived from this symmetric characteristic.After taking these sampling points, the precious relation between position, current, electromagnetic force and torque can be obtained from two-dimension interpolation method, thus a precious real-time control model can established.
The interval 0-2τ n is partitioned into four regions according to the symmetry of the magnetic field, as shown in Figure 3. Thus, the corresponding relation between the current and electromagnetic force and torque in other regions can be deducted from the interpolation results in region I as shown in Figure 3.
outcome of analytical modeling is behavioral models [27].For instance, Equations ( 17)-( 22) provide the analytical modeling for the electromagnetic force/torque of the maglev planar motor.The interpolation method, on the other hand, processes the discrete data with the numerical approach based on the known data, thus a complete mathematical model of the control target object can be derived.Normally, interpolation includes linear interpolation, polynomial interpolation, cubic spline interpolation, which extend from single dimension to multiple dimensions [28].To meet the real-time computation needs, both the accuracy and speed should be taken into account to choose a proper interpolation method; therefore, the two-dimension interpolation method has been chosen.

Interpolation Region Partition
Since the magnetic field of a Halbach array is strictly periodic symmetrical, it can be simplified to take points in a small range to perform the interpolation.The whole planar model can be derived from this symmetric characteristic.After taking these sampling points, the precious relation between position, current, electromagnetic force and torque can be obtained from two-dimension interpolation method, thus a precious real-time control model can established.
The interval 0-2τn is partitioned into four regions according to the symmetry of the magnetic field, as shown in Figure 3. Thus, the corresponding relation between the current and electromagnetic force and torque in other regions can be deducted from the interpolation results in region I as shown in Figure 3.

Interpolation Calculation
Bilinear interpolation method is adopted to interpolate in region I of Figure 3 with.After the partition in region I, the interpolation points (xg, yh) are located within: Then the value of interpolation point s(xg, yh) can be determined by the value on (xi, yj), (xi+1, yj), (xi, yj+1), (xi+1, yj+1).
The value on the interpolation point can be calculated: , where:

Interpolation Calculation
Bilinear interpolation method is adopted to interpolate in region I of Figure 3 with.After the partition in region I, the interpolation points (x g , y h ) are located within: Then the value of interpolation point s(x g , y h ) can be determined by the value on (x i , y j ), (x i+1 , y j ), (x i , y j+1 ), (x i+1 , y j+1 ).
The value on the interpolation point can be calculated: where: When interpolating with Equation ( 27), the feedback position of maglev planar motor's mover under global coordination system is s(x g , y h ), thus, the corresponding electromagnetic model can be obtained.

Current Decoupling
The input of the maglev planar motor is the current in each coil and the output is the force and torque on the mover.It is a multi-input and multi-output coupling system.Thus, the decoupling current control of each coil is extra critical to realize the overall performance of the control system.In the proposed maglev planar motor, there are 16 sets of independent coils.Each coil can generate the forces and torques in multiple degrees of freedom.The total force and torque on the mover mass center can be seen as the sum of those on the 16 coils.The force and torque on the mover can be expressed as: The current in each coil of the mover and the electromagnetic force and torque of the mover can be expressed as: where, p is the position vector; K(p) is a coefficient matrix about the electromagnetic force and torque relevant to the position.Each element in K(p) is the relation of the corresponding current and force, this relation can be obtained from analytical method or sampling and interpolation.The Moore-Penrose generalized inverse matrix of K(p) is shown as: The current in each coil can be obtained as:

Dynamic Model
The coil and the back plate of the moving-coil maglev planar motor are deemed as one object, which forms the mover of this motor, and this can be regarded as rigid motion.Based on the above analysis, the rigid body dynamic of the mover can be described as: where, m is the mover mass; g is the gravitational acceleration; I x , I y and I z are rotary inertia.

Control System Structure
After the decoupling of the six degree of freedom maglev permanent magnetic synchronous planar motor, six independent single degree of freedom control system can be derived.Thus, an independent controller applied on each single degree of freedom can achieve the corresponding closed loop control.Thereafter, the multi-degree of freedom control of the six degrees maglev motor can be achieved.
In order to achieve the six dimension control of the maglev planar motor, a closed loop control is needed for each dimension.The current can be calculated according to the required electromagnetic force and torque and the relative position of the magnetic field.In Figure 4, p ref is the given position vector; W ref is the needed force and torque on each degree of freedom; i ref is the needed current in each coils; W is the generated forces and torques on the six degrees of freedom; and p is the displacement vector of the planar motor.0 0 0 0 0 0 0 0 0 0 where, m is the mover mass; g is the gravitational acceleration; Ix, Iy and Iz are rotary inertia.

Control System Structure
After the decoupling of the six degree of freedom maglev permanent magnetic synchronous planar motor, six independent single degree of freedom control system can be derived.Thus, an independent controller applied on each single degree of freedom can achieve the corresponding closed loop control.Thereafter, the multi-degree of freedom control of the six degrees maglev motor can be achieved.
In order to achieve the six dimension control of the maglev planar motor, a closed loop control is needed for each dimension.The current can be calculated according to the required electromagnetic force and torque and the relative position of the magnetic field.In Figure 4, pref is the given position vector; Wref is the needed force and torque on each degree of freedom; iref is the needed current in each coils; W is the generated forces and torques on the six degrees of freedom; and p is the displacement vector of the planar motor.

Controller Design
A PID controller has the advantages of simple structure and being easy to control, and it can get the desired control results under certain external conditions, therefore, it is widely used in the industrial control field [29][30][31][32][33]. Through the PID controller, the system can process the proportional, integral and differential calculation, and output the sum of the three calculations as the control signal for the controlled object.According to the characteristics of the controlled object and the required performance, the controller can be conducted as P control, PI control and PID control.

Control Law Analysis
A maglev system is a non-linear and open loop unstable system, the linearized model consists of displacement stiffness coefficient and current stiffness coefficient, but what differs from the general spring damp system is that the displacement stiffness coefficient of this system is negative without feedback control.The controlled object can only be suspended and remain stable with

Controller Design
A PID controller has the advantages of simple structure and being easy to control, and it can get the desired control results under certain external conditions, therefore, it is widely used in the industrial control field [29][30][31][32][33]. Through the PID controller, the system can process the proportional, integral and differential calculation, and output the sum of the three calculations as the control signal for the controlled object.According to the characteristics of the controlled object and the required performance, the controller can be conducted as P control, PI control and PID control.

Control Law Analysis
A maglev system is a non-linear and open loop unstable system, the linearized model consists of displacement stiffness coefficient and current stiffness coefficient, but what differs from the general spring damp system is that the displacement stiffness coefficient of this system is negative without feedback control.The controlled object can only be suspended and remain stable with certain disturbance when imposing enough current on the coil.The adjustment and choose of the maglev system can be explained through Figure 5.In Figure 5a, the kinematics equation of the maglev system is: where, k > 0 and d > 0.
In Figure 5b, the kinematics equation of the maglev system is: Make the two systems equivalent: Therefore, the control law of the current is chosen as: ( ) where: Equation ( 43) demonstrates that in order to make the two systems to be equivilent, the current of the controlled maglev system is supposed to consist of at least a proportional control part and differntial control part.From the perspective of mechanical engineering, the proportional part of the current is used to offset the restoring force caused by the negative stiffness and change the system into positive stiffness.The differential part of the current is used to make the system of positive damping, and thus of sufficient stability.

Controller Struture
From the analysis of the previous section, in order to realize the stable suspension of the controlled object, there must be proportional and differential law in this closed loop controller.Meanwhile, to neglect and minimize the stable error of this system, integral control law can be added into this control system.
The time domain expression of this PID controller is: In Figure 5a, the kinematics equation of the maglev system is: where, k > 0 and d > 0.
In Figure 5b, the kinematics equation of the maglev system is: Make the two systems equivalent: Therefore, the control law of the current is chosen as: where: Equation ( 43) demonstrates that in order to make the two systems to be equivilent, the current of the controlled maglev system is supposed to consist of at least a proportional control part and differntial control part.From the perspective of mechanical engineering, the proportional part of the current is used to offset the restoring force caused by the negative stiffness and change the system into positive stiffness.The differential part of the current is used to make the system of positive damping, and thus of sufficient stability.

Controller Struture
From the analysis of the previous section, in order to realize the stable suspension of the controlled object, there must be proportional and differential law in this closed loop controller.Meanwhile, to neglect and minimize the stable error of this system, integral control law can be added into this control system.
The time domain expression of this PID controller is: The general transfer function of the PID controller is: where, K p is the proportional gain; K i is the integral gain; K d is the differential gain.The structure of this PID controller is shown in Figure 6, where r is the system reference value, u is the control value, G(s) is the controlled object, p is the system output.The proportional component provides the overall control of the different error e, the integral component minimizes the stable error through low frequency integral compensation, and differential component improves the dynamic response through high frequency differential compensation.
The general transfer function of the PID controller is: where, Kp is the proportional gain; Ki is the integral gain; Kd is the differential gain.The structure of this PID controller is shown in Figure 6, where r is the system reference value, u is the control value, G(s) is the controlled object, p is the system output.The proportional component provides the overall control of the different error e, the integral component minimizes the stable error through low frequency integral compensation, and differential component improves the dynamic response through high frequency differential compensation.

Validation and Analysis
The validation and analysis consists of three parts.First, the precision of the model is analyzed through 3-D finite element method (FEM) simulation.Then the precision of the control system is analyzed through MATLAB/Simulink simulation.Finally, the real time ability of this system has been validated and analyzed.The parameters of the planar motor for the validation and analysis are shown in Table 1.

Validation and Analysis
The validation and analysis consists of three parts.First, the precision of the model is analyzed through 3-D finite element method (FEM) simulation.Then the precision of the control system is analyzed through MATLAB/Simulink simulation.Finally, the real time ability of this system has been validated and analyzed.The parameters of the planar motor for the validation and analysis are shown in Table 1.

Model Accuracy Analysis
The 3-D simulation is applied to validate the accuracy of the analytical model and the interpolation model respectively.
In Figures 7 and 8, the displacement in the direction of x is from 0 mm to 36 mm, the same displacement length and the same sampling method are applied to the direction of y.The parameters for the simulation are shown in Table 1, and a current of 3 A is imposed on the coil.The amplitude Energies 2017, 10, 1132 11 of 17 of F x , F y , F z , T x , T y , and T z calculated by FEM are 12.39 N, 12.39 N, 25.67 N, 1.68 N•m, 1.69 N•m, and 1.83 N•m, respectively.It can be seen from Figure 7 that the maximum absolute value errors are 0.83 N for F x , 0.91 N for F y , 2.38 N for F z , 0.154 N•m for T x , 0.166 N•m for T y , and 0.148 N•m for T z .Figure 8 shows the error between the interpolation model and 3-D simulation.It can be seen from Figure 8 that the maximum absolute value errors are 0.18 N for F x , 0.16 N for F y , 0.23 N for F z , 0.033 N• m for T x , 0.038 N• m for T y , and 0.032 N•m for T z .It has been found from the comparison above that interpolation method has a better model accuracy than the analytical method.

Model Accuracy Analysis
The 3-D simulation is applied to validate the accuracy of the analytical model and the interpolation model respectively.
In Figures 7 and 8, the displacement in the direction of x is from 0 mm to 36 mm, the same displacement length and the same sampling method are applied to the direction of y.The parameters for the simulation are shown in Table 1, and a current of 3 A is imposed on the coil.The amplitude of Fx, Fy, Fz, Tx, Ty, and Tz calculated by FEM are 12.39 N, 12.39 N, 25.67 N, 1.68 N•m, 1.69 N•m, and 1.83 N•m, respectively.It can be seen from Figure 7 that the maximum absolute value errors are 0.83 N for Fx, 0.91 N for Fy, 2.38 N for Fz, 0.154 N•m for Tx, 0.166 N•m for Ty, and 0.148 N•m for Tz. Figure 8 shows the error between the interpolation model and 3-D simulation.It can be seen from Figure 8 that the maximum absolute value errors are 0.18 N for Fx, 0.16 N for Fy, 0.23 N for Fz, 0.033 N•m for Tx, 0.038 N•m for Ty, and 0.032 N•m for Tz.It has been found from the comparison above that interpolation method has a better model accuracy than the analytical method.

Model Accuracy Analysis
The 3-D simulation is applied to validate the accuracy of the analytical model and the interpolation model respectively.
In Figures 7 and 8, the displacement in the direction of x is from 0 mm to 36 mm, the same displacement length and the same sampling method are applied to the direction of y.The parameters for the simulation are shown in Table 1, and a current of 3 A is imposed on the coil.The amplitude of Fx, Fy, Fz, Tx, Ty, and Tz calculated by FEM are 12.39 N, 12.39 N, 25.67 N, 1.68 N•m, 1.69 N•m, and 1.83 N•m, respectively.It can be seen from Figure 7 that the maximum absolute value errors are 0.83 N for Fx, 0.91 N for Fy, 2.38 N for Fz, 0.154 N•m for Tx, 0.166 N•m for Ty, and 0.148 N•m for Tz. Figure 8 shows the error between the interpolation model and 3-D simulation.It can be seen from Figure 8 that the maximum absolute value errors are 0.18 N for Fx, 0.16 N for Fy, 0.23 N for Fz, 0.033 N•m for Tx, 0.038 N•m for Ty, and 0.032 N•m for Tz.It has been found from the comparison above that interpolation method has a better model accuracy than the analytical method.After the decoupling, the originally strongly decoupled system becomes a weakly decoupled system.However, the decoupling performance depends on the modeling accuracy.The model decoupling and its control parts together determine the position loop accuracy of the whole control system.Under a given position loop accuracy, a high modeling accuracy can make the other parts in the control system, like the decoupling and noise filtering, relatively easier and thus more efficient.Therefore, it is easier to realize the desired position loop accuracy goal.

System Precision Analysis
Through MATLAB/Simulink software, the control system is established as shown in Figure 4, and it consists of six PID controllers to realize the control on the six degrees of freedom.The parameters of motor are as shown in Table 1.The simulation time is set as 0.1 s, and 1 mm displacement reference order on z-direction is given at 0.01 s, which is the working suspension height of the maglev planar mover.10 mm position step signal on y-direction is given at 0.03 s, and the reference order remains to be 0 on the rest directions to observe the decoupling performance.Upon the simulation above, the displacement curve and the force/toque curve from the control system based on the traditional analytical method is shown in Figures 9 and 11.The displacement curve and the force/torque curve from the control system based on the two-dimension linear interpolation method are shown in Figures 10 and 12.
After the decoupling, the originally strongly decoupled system becomes a weakly decoupled system.However, the decoupling performance depends on the modeling accuracy.The model decoupling and its control parts together determine the position loop accuracy of the whole control system.Under a given position loop accuracy, a high modeling accuracy can make the other parts in the control system, like the decoupling and noise filtering, relatively easier and thus more efficient.Therefore, it is easier to realize the desired position loop accuracy goal.

System Precision Analysis
Through MATLAB/Simulink software, the control system is established as shown in Figure 4, and it consists of six PID controllers to realize the control on the six degrees of freedom.The parameters of motor are as shown in Table 1.The simulation time is set as 0.1 s, and 1 mm displacement reference order on z-direction is given at 0.01 s, which is the working suspension height of the maglev planar mover.10 mm position step signal on y-direction is given at 0.03 s, and the reference order remains to be 0 on the rest directions to observe the decoupling performance.Upon the simulation above, the displacement curve and the force/toque curve from the control system based on the traditional analytical method is shown in Figures 9 and 11.The displacement curve and the force/torque curve from the control system based on the two-dimension linear interpolation method are shown in Figures 10 and 12   In Figure 9, the position disturbance on x-axis caused by y-axis and z-axis position variation is from −0.081 μm to 0.076 μm, the position disturbance on ϕ-axis is from −0.51 μrad to 0.73 μrad, the After the decoupling, the originally strongly decoupled system becomes a weakly decoupled system.However, the decoupling performance depends on the modeling accuracy.The model decoupling and its control parts together determine the position loop accuracy of the whole control system.Under a given position loop accuracy, a high modeling accuracy can make the other parts in the control system, like the decoupling and noise filtering, relatively easier and thus more efficient.Therefore, it is easier to realize the desired position loop accuracy goal.

System Precision Analysis
Through MATLAB/Simulink software, the control system is established as shown in Figure 4, and it consists of six PID controllers to realize the control on the six degrees of freedom.The parameters of motor are as shown in Table 1.The simulation time is set as 0.1 s, and 1 mm displacement reference order on z-direction is given at 0.01 s, which is the working suspension height of the maglev planar mover.10 mm position step signal on y-direction is given at 0.03 s, and the reference order remains to be 0 on the rest directions to observe the decoupling performance.Upon the simulation above, the displacement curve and the force/toque curve from the control system based on the traditional analytical method is shown in Figures 9 and 11.The displacement curve and the force/torque curve from the control system based on the two-dimension linear interpolation method are shown in Figures 10 and 12   In Figure 9, the position disturbance on x-axis caused by y-axis and z-axis position variation is from −0.081 μm to 0.076 μm, the position disturbance on ϕ-axis is from −0.51 μrad to 0.73 μrad, the In Figure 9, the position disturbance on x-axis caused by y-axis and z-axis position variation is from −0.081 µm to 0.076 µm, the position disturbance on ϕ-axis is from −0.51 µrad to 0.73 µrad, the position disturbance on θ-axis is from −0.23 µrad to 0.21 µrad, the position disturbance on φ-axis is from −0.38 µrad to 0.35 µrad.In Figure 10, because of the position variation of y-axis and z-axis, the position disturbance of x-axis is from −0.038 µm to 0.004 µm, the position disturbance of ϕ-axis is Energies 2017, 10, 1132 13 of 17 from −0.003 µrad to 0.12 µrad, the position disturbance of θ-axis is from -0.035 µrad to 0.013 µrad, the position disturbance of φ-axis is from −0.058 µrad to 0.078 µrad.
Figure 11 shows the force and torque disturbance caused by y-axis and z-axis position variation: for F x , it varies from −2.65 N to 2.48 N; For T x , it varies from −1.86 N•m to 0.79 N•m; For T y , it varies from −0.91 N•m to 0.16 N•m; For T z , it varies from −0.27 N•m to 0.95 N•m.Due to the saturation restriction of the controller, the forces on y-axis and z-axis are limited within ±1000 N. In Figure 12, it shows the force and torque disturbance caused by y-axis and z-axis position variation: for F x , it varies from −0.57N to 0.13 N; For T x , it varies from −0. position disturbance on θ-axis is from −0.23 μrad to 0.21 μrad, the position disturbance on φ-axis is from −0.38 μrad to 0.35 μrad.In Figure 10, because of the position variation of y-axis and z-axis, the position disturbance of x-axis is from −0.038 μm to 0.004 μm, the position disturbance of ϕ-axis is from −0.003 μrad to 0.12 μrad, the position disturbance of θ-axis is from -0.035 μrad to 0.013 μrad, the position disturbance of φ-axis is from −0.058 μrad to 0.078 μrad.
Figure 11 shows the force and torque disturbance caused by y-axis and z-axis position variation: for Fx, it varies from −2.65 N to 2.48 N; For Tx, it varies from −1.86 N•m to 0.79 N•m; For Ty, it varies from −0.91 N•m to 0.16 N•m; For Tz, it varies from −0.27 N•m to 0.95 N•m.Due to the saturation restriction of the controller, the forces on y-axis and z-axis are limited within ±1000 N. In Figure 12, it shows the force and torque disturbance caused by y-axis and z-axis position variation: for Fx, it varies from −0.57N to 0.13 N; For Tx, it varies from −0.  Figure 13 shows the difference of the forces and torques in FEM, analytical method and interpolation method.It has been found that the interpolation method is more close to the FEM than the analytical method.Through the comparison of the simulation results, it can be seen that the control system based on the interpolation method takes into account not only the model simplification, but also the error of materials manufacturing, processing and installation, which makes the established electromagnetic model reflect the real physical model better.Therefore, the established model is of higher position precision.position disturbance on θ-axis is from −0.23 μrad to 0.21 μrad, the position disturbance on φ-axis is from −0.38 μrad to 0.35 μrad.In Figure 10, because of the position variation of y-axis and z-axis, the position disturbance of x-axis is from −0.038 μm to 0.004 μm, the position disturbance of ϕ-axis is from −0.003 μrad to 0.12 μrad, the position disturbance of θ-axis is from -0.035 μrad to 0.013 μrad, the position disturbance of φ-axis is from −0.058 μrad to 0.078 μrad.
Figure 11 shows the force and torque disturbance caused by y-axis and z-axis position variation: for Fx, it varies from −2.65 N to 2.48 N; For Tx, it varies from −1.86 N•m to 0.79 N•m; For Ty, it varies from −0.91 N•m to 0.16 N•m; For Tz, it varies from −0.27 N•m to 0.95 N•m.Due to the saturation restriction of the controller, the forces on y-axis and z-axis are limited within ±1000 N. In Figure 12, it shows the force and torque disturbance caused by y-axis and z-axis position variation: for Fx, it varies from −0.57N to 0.13 N; For Tx, it varies from −0.    Figure 13 shows the difference of the forces and torques in FEM, analytical method and interpolation method.It has been found that the interpolation method is more close to the FEM than the analytical method.Through the comparison of the simulation results, it can be seen that the control system based on the interpolation method takes into account not only the model simplification, but also the error of materials manufacturing, processing and installation, which makes the established electromagnetic model reflect the real physical model better.Therefore, the established model is of higher position precision.Figure 13 shows the difference of the forces and torques in FEM, analytical method and interpolation method.It has been found that the interpolation method is more close to the FEM than the analytical method.Through the comparison of the simulation results, it can be seen that the control system based on the interpolation method takes into account not only the model simplification, but also the error of materials manufacturing, processing and installation, which makes the established electromagnetic model reflect the real physical model better.Therefore, the established model is of higher position precision.

Model Real-Time Verification
We programmed the proposed interpolation model and ran it on a NI PXI 1042Q real time platform to verify the real time ability.Figure 14a

Model Real-Time Verification
We programmed the proposed interpolation model and ran it on a NI PXI 1042Q real time platform to verify the real time ability.Figure 14a shows the real time control platform, including an industrial PC, embedded controller, plug-in I/O module, a PC and a display.PC connects the real time controller via a cable, and the PC is used to write LabVIEW program and download it to the controller.Controller needs to complete the real time computation of the program, read data from I/O board and upload data to the PC.
Figure 14a shows the measured results of the running time for each loop.From Figure 14b, it can be seen that the average running time of the proposed model based on two-dimensional interpolation is 200 µs and the maximum value is 232 µs.
We programmed the proposed interpolation model and ran it on a NI PXI 1042Q real time platform to verify the real time ability.Figure 14a  Table 2 shows the computation time of each modeling method.The model size of FEM is 0.2 m × 0.2 m × 0.1 m, which has been further divided into about 3 mm 2 units.The simulation time is about 1 h.Numerous points need to be simulated to get the force on a certain plane.FEM takes a large amount of time and thus it is suitable for theoretical validation.The computation times for the analytical method and interpolation method are quite similar.On the NI control platform it is about 200 µs.For real-time control, it is normally required that the computation time be less than 1 ms.Thus, both of the two above methods are suitable for the real-time control system in industrial applications.However, the interpolation method is more accurate than the analytical method.

Conclusions
This paper proposes an electromagnetic model establishing method based on two-dimensional linear interpolation, and designs a maglev planar motor control system accordingly.The proposed model can reduce the model simplification and the manufacturing error, and thus better represents the real physical model.The 3-D FEM simulation approach is used to validate the proposed model error and the analytical model error, respectively.The validation result shows that the proposed model is found to have smaller model error.Adopting MATLAB/Simulink to build the control system, comparing the control system position precision based on the two models, it shows the control system based on the two-dimensional linear interpolation method is of higher position precision.Meanwhile, NI-PXI platform is used to validate the real-time ability of the control model, the average running time is 200 µs.Thus, it can be concluded that the proposed control model is more suitable to be applied in the high speed and high precision maglev planar control system.

Figure 2 .
Figure 2. Top views and A-A view of the components of a maglev planar motor: a Halbach array and a coil.

Figure 2 .
Figure 2. Top views and A-A view of the components of a maglev planar motor: a Halbach array and a coil.

Figure 2 .
Figure 2. Top views and A-A view of the components of a maglev planar motor: a Halbach array and a coil.

Figure 5 .
Figure 5.Comparison of the mechanical system and the maglev system.(a)Mechanical system; (b) Maglev system.

Figure 5 .
Figure 5.Comparison of the mechanical system and the maglev system.(a) Mechanical system; (b) Maglev system.

Figure 7 .
Figure 7. Force and torque errors between analytical method and FEM.

Figure 8 .
Figure 8. Force and torque errors between interpolation method and FEM.

Figure 7 .
Figure 7. Force and torque errors between analytical method and FEM.

Figure 7 .
Figure 7. Force and torque errors between analytical method and FEM.

Figure 8 .
Figure 8. Force and torque errors between interpolation method and FEM.Figure 8. Force and torque errors between interpolation method and FEM.

Figure 8 .
Figure 8. Force and torque errors between interpolation method and FEM.Figure 8. Force and torque errors between interpolation method and FEM. . μm

Figure 9 .Figure 10 .
Figure 9. Displacement curve of the control system based on the analytical method.

Figure 9 .
Figure 9. Displacement curve of the control system based on the analytical method. . μm

Figure 9 .Figure 10 .
Figure 9. Displacement curve of the control system based on the analytical method.

Figure 10 .
Figure 10.Displacement curve of the control system based on the interpolation method.
13 N•m to 0.07 N•m; For T y , it varies from −0.005 N•m to 0.086 N•m; For T z , it varies from −0.089 N•m to 0.052 N•m.Energies 2017, 10, 1132 13 of 17

Figure 11 .
Figure 11.Force and torque curve of the control system based on the analytical method.

Figure 12 .
Figure 12.Force and torque curve of the control system based on the interpolation method.

Figure 11 .
Figure 11.Force and torque curve of the control system based on the analytical method.
13 N•m to 0.07 N•m; For Ty, it varies from −0.005 N•m to 0.086 N•m; For Tz, it varies from −0.089 N•m to 0.052 N•m.

Figure 11 .
Figure 11.Force and torque curve of the control system based on the analytical method.

Figure 12 .
Figure 12.Force and torque curve of the control system based on the interpolation method.

Figure 12 .
Figure 12.Force and torque curve of the control system based on the interpolation method.

Figure 13 .
Figure 13.Forces /torques of a coil in FEM, analytical method and interpolation method.

Figure 13 .
Figure13.Forces /torques of a coil in FEM, analytical method and interpolation method.

Figure 14 .
Figure 14.Program running time analysis.(a) Real-time control platform based on NI PXI-1042Q; (b) Results of the program running time.

Table 1 .
Main parameters of the planar motor.

Table 1 .
Main parameters of the planar motor.

Table 2 .
Computation time of FEM, interpolation method and analytical method.