Mathematical models and nonlinear dynamics of a linear electromagnetic motor

The aim of this paper was to construct a mathematical model of a drive system with ironless permanent magnet synchronous linear motor and to animalize the influence of the magnetic field distribution function on the model’s accuracy and ease of simulation computation. The studied motor employs an U-shaped stationary guideway with permanent magnets placed perpendicularly to the motor’s direction of motion and a forcer with three sets of rectangular coils subjected to alternating external electrical voltage. The system’s parameters are both mechanical (number of magnets and coils, size of magnets, distances between magnets, size of coils) and electromagnetic (auxiliary magnetic field, permeability, coil’s resistance). Lorentz force allows for the transition from electromagnetic parameters to mechanical force, and Faraday’s law of induction creates a feedback between the forcer’s speed and coils voltage. An Ampere’s model of permanent magnet is used to determine the function of auxiliary magnetic field distribution throughout the stator. Two simplified distribution functions are introduced and studied. During validation, external current function is applied to each coil (serving as excitation), while the displacement of forcer in time is the output function. Model parameters are found via genetic algorithms such that the numerical solution of the model best fits experimental data. Several cases of motor operation are compared against simulation results showing good coincidence between computation and experiment.


Introduction
It is well known and documented that linear machines, including motors, generators, and actuators, generate motion/translation directly without rotation and transmission conversion devices, and due to their compact, simple and relatively cheep structure with a simultaneous high dynamic efficiency and performance, they are extensively applied in various branches of industry. They are employed in manufacturing [1], transportation [2], aerospace industry [3], automation, and robotic industries [4,5], bioengineering, and other.
However, the proper design of the linear machines producing large forces and operating at high speeds requires advanced mathematical modelling in order to achieve precise control to carry out required performance. On the other hand, a proper mathematical modelling and derivation of the reliable governing dynamic equations requires a comprehensive and self-contained knowledge that span from mechanics, applied mathematics, electrical, engineering, and mechatronics with a particular emphasis on nonlinear effects.
In what follows, we briefly describe the state of the art of the problem by exhibition of the derived and employed simple mathematical models based on fundamental physical/mechanical lows and supported by experimental validation.
Karnopp [6] studied linear electromagnetic motors composed of coils of a copper wire interacting with permanent magnets serving as mechanical dampers. In particular, novel coreless design for use with high energy magnets was proposed.
Compter and Frissen [7] presented an innovative planar motor levitating and propelling platform. It was composed of a simple structure with four three-phase amplifiers. Full 6-DoF single-input-single-output control was developed, and a stable controlled electrodynamical planar motion with a long stroke was obtained.
Holmes et al. [8] designed, fabricated, and tested an axial-flux permanent magnet electromagnetic generator. The latter consisted of a polymer rotor with embedded permanent magnets sandwiched between two silicon stators. The effects on performance of design parameters (number of layers in the stator coils, the rotor stator gap, soft magnetic pole pieces on the stators) were investigated using the finite element simulations. The employed microfabricated axial-flow microturbine was used to produce a compact power conversion device for power generation and flow sensing operations.
Lemarquand and Lemarquand [9] considered the structures of permanent magnet pickups for an electric guitar with an emphasis put on the nonlinear dynamics generated by the induced electromotive force yielded by the string, magnet, and coil interactions. The study was supported by the analytical computations using the Coulombian model of magnets.
Yatchev and Ritchie [10] compared two approaches for simulation of dynamics of a permanent magnet linear actuator, i.e. the so-called full coupled model (standard) and the decoupled model (proposed by the authors). It was illustrated that the decoupled model of a linear actuator with moving permanent magnet employing bicubic spline approximations yielded only the required accuracy but demonstrated more flexibility in comparison with the standard approach.
Hamzehbahmani [11] derived a mathematical model for a single-side short-stator linear induction motor with the end effects. MATLAB/SIMULINK-based simulations were compared with the experimental results. The latter included measurement of the equiv-alent circuit parameters, both free acceleration and locked primary experiments.
The voice coil motor consisted of a wounded coil in the centre covered with Halbach magnet array with iron yoke was mathematically modelled with Lorentz force principle by Jeong et al. [12]. Relation between coil-passing area and suppressing parasitic motion was derived and the motor equation was expressed by geometric dimensions of coil and magnets, which opened a possibility for its design.
Yao et al. [13] designed based on the mathematical modelling a novel tubular linear machine with 3D hybrid permanent magnet arrays and multiple movers. The governing equations were derived using the source free property and Maxwell equations. The magnetic field distribution was analytically formulated using Bessel functions and harmonic approximation to the magnetization vector. The obtained analytical solutions were validated by numerical simulations. In addition, it was demonstrated that the variation of magnetic flux density in the linear machine was consistent with the magnet patterns.
Mathematical modelling and dynamic analysis of a novel moving magnet linear actuator were carried out by Hassan et al. [14]. The stator was composed of two reversely wound coils being electrically excited with single-phase AC power. The time-dependent current model of the stator winding was proposed, and the stroke, velocity, and acceleration of the armature were derived. The spring, damping, inertial and magnetic forces were estimated. The system dynamics were analytically approximated and experimentally validated.
As it has been already pointed out, the linear motors are taking on bigger and bigger shares of the market for precise positioning systems. They provide a costly but dynamically superior alternative to standard drives such as feed screw conveyors. In order to correctly project a desired motion profile, every positioning system has to equip a specific closed-loop controller between the motor and mains. The quality of such system, and therefore the quality of motion projection, depends on, amongst other things, the precision of the system elements' mathematical models used on the programing stage of building the controller. Most industrial controllers today use simplified models. This paper focuses on construction of a novel model for an exemplary positioning system which is different from standard available in the existing literature and has higher level of complexity. This paper is organized in the following way. Section 3 deals with basic electromagnetic interactions between elements of the motor. It presents a base model of a single winding in a vicinity of a C-shaped magnet and reports the values of Lorentz force acting on it as well as the value of electromagnetic induction. It then describes the magnetic distribution function of a single magnet with use of Ampere's model and Biot-Savart law. In Sect. 4, a model of a three-phase motor, highlighting the impact of choosing a proper magnetic field distribution function for a guideway with numerous magnets, is proposed. Two functions and therefore two models are presented. The first on simplifies the distribution to a sine function, while the second attempts for a more true to life distribution. Flaws of both are named, and for the simplified function, a scope of use is listed and explained. Section 5 presents the method of models' validation introducing final changes to models' equations. A resistive force of mechanical guideways is defined, and friction model is chosen to best suit the free-wheeling operations. Then, both models' parameters are identified using genetic algorithm and verified against experiment. Section 6 contains concluding remarks.

Validation platform
A base for model construction was a laboratory stand, with HIWIN's coreless linear motor and Copley Controls' servo drive. The same platform is also used for model validation. The stand's forcer (inductor), shown as "1" in Fig. 1, allows the load to move at speeds of up to 5 meters per second and generate a continuous and peak force of 45N and 180N, respectively. The length of three magnetic U-shaped guideway ("2" of Fig. 1) allows for a theoretical maximum stroke of 830 mm. For safety measures, the actual maximum stroke is limited to 780 mm. Two high-quality linear guideways (denoted by "5" in Fig. 1) provide for a reliable and quiet motion. The position and velocity measurement is done by a Renishaws analogue optical linear encoder (see "3" in Fig. 1). It is able to read its position with a resolution of 0.1 µm, and the entire system's positioning error is not trailing far behind. The stand program allows it to work both automatically based on a signal from external devices and via manual input. This stand has been designed and originally built [15].
The servo drive allows for a multitude of options to fine tune the motor's control loops. Amongst them is the "scope" function, which gives the option to subject the motor and the entire system to several functions of current, velocity, and position including sinusoidal function, square function and microstepping. It also measures the actual time graphs of said quantities and allows exporting them to a csv file. In this paper, the "scope" function is used to obtain data employed for model validation.

Base model
In the most basic model of a linear motor, the inductor can be represented as a single, perfectly rectangular conductor loop. This loop allows for a certain voltage function U g to be applied to it and cause a current of value I to flow through it. The exact means of this application is unimportant for the workings of the model and therefore will be omitted.
The magnetic guideway provides a source of magnetic field. A typical U-shaped guideway is composed of several pairs of strong magnets. Each pair's magnetic field distribution can be modelled as that of a single Cshaped magnet in a way depicted in Fig. 2 Assuming that the loop cannot deform or rotate and can only move along the direction of y-axis, it shall always remain a rectangle with a centre placed on yaxis and with shorter sides parallel to the direction of motion.
The described model is a base model for the paper. It should be noted that the actual linear motor bears many similarities to the base model and the only differences are number and displacement of magnets, number of conductor loops wound together within a single winding and number and displacement of windings.
This section focuses on the interactions within the base model, noting the differences and appropriate changes in the equations where the entire motors model is considered.

Lorentz force
Every charged particle moving through a magnetic field with a given speed is subjected to a Lorentz force, which works perpendicular to both the velocity vector of the Fig. 1 The laboratory stand used for model verification (1-motor's inductor, 2-magnetic guideway, 3-linear encoder, 4cable chain, 5-mechanical guideways)

Fig. 2
Model of a single winding particle and magnetic distribution vector. Same force works on a conductor placed within a magnetic field. In such a case, the force is perpendicular to the vector of the current flowing through the conductor [16].
For the base model, only the y-component of the Lorentz force is relevant (where y is the axis of motion). It can be obtained by calculating a closed-loop integral along the conductor path. For a specific case of a rectangular shape, the closed-loop integral can be rewritten as a sum of two definite integrals in the following form where B(x, y) is the function of the magnetic field's zcomponent distribution, l x , l y are the dimensions of the rectangular conductor loop in x an y axes, respectively, and y 0 is the position of the loop relative to centre of the magnet [17]. A single motor winding consists of n of such conductor loops placed a top each other along z-axis. Assuming that B distribution has the same values for each of the loop, the force for an entire single winding can be calculated as For a motor model, each coil's magnetic field also acts on adjacent windings; however, the sum of this forces for the entire motor is equal to zero. The stiffness of the inductor disallows movement of windings relative to each other, and therefore, these forces are omitted in the model.

Electromagnetic induction
When a closed conductor loop moves through dense magnetic field, the magnetic flux inside it changes. This causes an EMF (electromagnetic force) to induce inside the loop. The value of EMF for the base model can be calculated as follows where i is the EMF and l x , l y are the dimensions of the rectangular conductor loop in x an y axes, respectively. Similarly for a rectangular winding with n loops it can be recast to the following form The electric current I flowing through the rectangular coil with n conductor loops is the result of the external voltage function U G (t), the electromagnetic force i induced in the loop due to external magnetic flux change and the electromagnetic force caused by self inductance. It can be calculated by solving the following ODEİ where Re is the electric resistance of the coil and L is its leakage inductance. For a model with a set of coils in close proximity additional EMF resulting from windings' mutual inductance must be taken into consideration.

Magnetic field distribution
In accordance with reference [18], the value of magnetic field of a permanent magnet can be approximated  by using Amperes model, that is, by assuming that a magnets magnetic field is the same as that of a perfect, tightly wound solenoid, with current I S flowing through it. Then, Biot-Savart law can be applied to calculate the exact value of magnetic field at any point P in the vicinity of the magnet via the following formula where and r is the distance between point P and the centre of the magnet, l is the distance between the magnets centre and the infinitely small length of the solenoid dl, while μ is the magnetic permeability of the magnets environment.
In case of a C-shaped magnet, the field can be calculated as a resultant fields of two bar magnets placed perpendicular to each other with opposing poles facing each other. For a single bar magnet, the field can be then calculated by solving the following integral equation where and x p , y p , z p are the coordinates of point P, σ x , σ y , σ z are the bar magnet dimensions in their respective axes,x,ŷ,ẑ are the unit vectors of axes x, y and z, and i is the current density over the sides of the coil. Figure 3 shows the C-shaped magnets magnetic fields z-component distribution in a in y-axis (for z = 0) calculated with Eqs. (8)-(12) for a sample magnet.
For the ease of use and because of the good coincidence with actual values, the C-shaped magnet's magnetic fields distribution on a xy-plane (coordinate system set as in Fig. 2) will be approximated with 2D Gaussian function in the following form where B 0 is the magnet's given constant.
where y,ẏ,ÿ are, respectively, the position, velocity, and acceleration of the pusher, R is the resistive force of the guideways, F M j are the force functions of each of the motor's winding governed by Eqs. (2) and (3), y j is the distance between two adjacent phases, and U G j are the voltages on each winding. As mentioned, the stands magnetic guideway consists of several dozen C-shaped neodymium magnets. They are placed in such a way that each two neighbouring magnets have their poles placed oppositely and that the distance between magnets is the same for each pair and equal to χ (see Fig. 4).
For such a guideway, the magnetic field distribution function B(x, y) in Eq. (2) is a superposition of each individual magnet's field and can be calculated as a sum of magnetic distribution functions of every magnet. It should be note that the physical pusher, due to mechanical limitations, can only move in such a way that all the windings are enclosed inside the guideway. Therefore, the model will not loose any of its accuracy if an infinitely long guideway model is used, bearing in mind that the actual stands stroke is still limited and every responses with a pushers displacement greater than the stroke of the stand is physically unachievable in real life.
The infinitely long magnetic guideway field distribution can be calculated as a sum of distribution from infinitely many magnets, that is However Eq. (15) is cumbersome and provides significant complications for use in the stands ODE, making it highly nonlinear and difficult for numerical simulations. It can, however, be casted and simplified to a more computation friendly form. This paper will present two twin motor models, built based on Eq. (15).

Motor with sinusoidal field distribution
For different values of σ y and χ , the shape of the curve described by Eq. (15) changes significantly along the y-axis. The x-axis distribution is only dependent on the σ x and B 0 parameters. Let the distribution function be fixed in x = 0. In such case, the one argument function of distribution is given as follows Let there be a correlation between σ y and χ given in a form of a dimensionless parameter k such that In such case, the shape of the magnetic field distribution depends solely on the value of k. Figure 5 presents the magnetic field distribution in y-axis for varying k. The obtained plot affirms that B is periodic for every studied value of k (where k = 0.5 means that there is no free space between adjacent magnets and k = 0 means that magnets are infinitely apart.) For lower k values, magnetic distribution from individual magnets is more prominent and typical Gaussian distribution peaks can be observed with large ranges of y for which B ≈ 0 in between. Larger k produces a distribution where individual magnetic fields are distorted by adjacent magnets producing a net distribution similar to a sine wave.
It can easily be seen that for a certain value of k the distribution of B can be approximated with satisfactory precision by a cosine function in the form of The inaccuracy of such an approximation, at a given point along the guideway, is the square of difference between values of magnetic distributions given by Eqs. (16) and (18) and is therefore a function of y and k. As both functions rely on the magnitude factor B 0 in a similar fashion, the relative error in the form of can be used to measure the quality of cosine approximation in any given point. The absolute error of the cosine approximation is a function dependent solely on k coefficient. It can be obtained by calculating the sum of errors from each point along the guideway. For continuous function R and a theoretical infinitely long guideway, this sum becomes an improper integral in the form of In this paper, we assume that the scope of use for cosine approximation is reserved only to the motor model instances where A achieves a minimum value. The absolute error is theoretically a function of both k and χ ; however, in case of the latter it can only achieve minimum when χ approaches infinity, i.e. for a single infinitely long magnet. Because of that, only the impact of k on the error function is analysed.
The A function takes minimal values for k = k 0 , where k 0 is a value for which the first derivative of A is equal to 0 and the second derivative is greater then 0. Using the Leibniz integral rule [19,Chapter 8], the sum rule for integration and the value for the square of sum, the value of the first derivative can be obtained in the form of where Fig. 14 Comparison between simulation and experimental data during model identification It can be noted that for nonzero χ , as is the case with physical linear motors, the roots of Δ I do not depend on its value. The sums in Eq. (21) do not converge to a finite value; however, the result of Eq. (21) can be approximated numerically after substituting the infinite summations to a finite ones with boundaries of −N and N , where N is a natural number. The larger the N, the longer the computation time needed to approximate and the more accurate the approximation. For every N there exists a k = k N meeting the following criteria With the results of several numerically computed k N for different N , an interpolation functions k I (n), having all positive real values in its domain, can be constructed. The interpolated function has to have the property that N ∈N k I (N ) ≈ k N and lim n→∞ k I (n) = k I 0 ≈ k 0 . Figure 6 presents a contour plot for a Δ(n, k) function (with the values for n / ∈ N interpolated from two adjacent point for the sake of graph clarity). The red line on the graph joins points where Δ reaches zero, i.e. points with coordinates (k N , N ). As best fitting the results, the k I (n) function should take an exponential form of where the values of parameters a 0 , a 1 , a 2 , a 3 can be obtained with the least squares method. They are equal to: a 0 = 0.4594, a 1 = 0.1547, a 2 = 0.4118, a 3 = 0.4077. The limit of the exponential interpolation function is equal to lim n→∞ k I (n) = a 0 ; therefore, the scope of use for sine model is motors with the correlation factor k equal to 0.4594. The same can be proven true for a 2D function of magnetic distribution given by Eq. (15) and a 2D cosine approximation function in the form of The absolute error of the 2D approximation, equal to after integrating over x and dividing over the constants, reduces to the form given by Eq. (20). Therefore, the approximation is valid for the same values of k. For motors matching that criteria, the integrals in Eqs. (2) and (3), i.e. values of the Lorentz force and EMF, can be analytically calculated as where K f is the winding's force constant, equal to By introducing Eq. (28) to the stand's model's motion's equation (Eq. (14)) and calculating the values of currents on each coil with Eqs. (5) and (29), the following system of equations for the laboratory stand is obtained In Eq. (31), U G j is the external voltage function and K V , K e , K U , ω g and α pj are the motor's back EMF, electrical time and voltage constants and guideway density, windings displacement factors, respectively, defined as where L and R E are the leakage inductance and electrical resistance of each winding.

Motor with elliptic theta function field distribution
Another method of using Eq. (15) in laboratory stands ODE is to insert it directly, but after enclosing the infinite sum of Gaussian distributions inside a single, more computational friendly function. This can be done by first ridding Eq. (15) of the (−1) i and instead splitting the sum into the difference of two other sums, each of them representing the resultant field distribution of same faced magnets in the guideway. The resultant sums can be recast to a form in which a single Jacobi theta functions of the third kind is used for each of them. This paper follows the Whittaker and Watson [20, Chapter 22] notation of elliptic functions; therefore, the ϑ 3 symbol will be defined as Although the elliptic theta is a complex function, while both the argument z and nome q are complex numbers (with q restrained only to the upper half-plane), it becomes real if the argument and nome are real and when q ∈ (−1, 1). Figure 7 presents the plot of ϑ 3 for real z values and for different values of q.
Equation (15) can then be presented in the form where It can be seen that the nome is the same for both theta functions and is a function of the correlation parameter k [as defined in Eq. (17)]. The z 0 parameter, which affects the frequency of the Jacobi functions, is also equal for those functions. The π 2 component acts as a phase parameter shifting the negatively signed function by exactly half of its period. An auxiliary function Θ can be introduced to the equation for greater clarity. This function is defined as follows Similarly, a Θ 0 (z, q) function, used in stand's model construction to incorporate the size of coils, can be defined as where z L = l y z 0 /2 is a winding length factor for a motor with windings' sizes of l x , l y . By introducing Eqs. (34) and (37) to laboratory stand's equation of motion (Eq. (14)) and the currents' equation (Eq. (5)) as well as defining the motor's theta force constant, theta back emf constant, and windings displacement factor in the form of the governing ODE take the form of R (ẏ, y) ,

Model identification and validation
Following their construction, both models were validated against the physical motor. However, the stand's servo drive disallows subjecting motors to voltage functions directly and instead provides current functions as the source for external excitation. Because of that, the identification and validation was not done using a known voltage function (such as a sinusoidal or square wave), but instead by measuring both the currents in each phase [I j in Eqs. (31) and (39)] and the displacement of motor y. The models were then supplied with the measured current functions, and their parameters were found, using genetic methods, such that the difference between the values of displacement for the model and physical motor was minimal. With parameters found, several numerical simulations and experimental measurements were taken to further validate the models. It has to be noted here that the models need additional shift parameter to include the difference between the position y = 0 for the servo drive and for magnetic field distribution functions Eqs. (34) and (18).

Resistive force
Prior to identifying the model's parameters, the nature of system's resistive force should be discerned. In order to do so, the free-wheeling method was employed. The stand's motor was cut from any external excitation voltage and windings were kept electrically separated so to avoid the induction of EMF. The pusher was then manually accelerated to a fixed speed and released allowing it to move solely on inertia. By incorporating the stands servo drive in the measurement process, a precise values of motor's free-wheeling times as well as speed and displacement functions during that time were obtained. To further enhance the accuracy of collected data, the same process was repeated for different initial speeds and with additional weights added to the pusher. Figures 9, and 10 show a sample velocity and position graphs for three free-wheeling operations with different pusher weights. It can be noted that unlike the position graph the velocity profile is significantly jagged and overshot appears to manifest near the end of motion (especially noticeable during free-wheeling with 1 kg added mass). Obviously, those phenomena should not be present during free-wheeling with no closed-loop control. Instead, they are the result of the servo drive's measurement system specifics and the fact that there exists no method of accurately measuring the linear motor's velocity. The results shown in Fig. 9 are actually numerically calculated derivatives of motor's displacement (Fig. 10), as such they carry substantial noise and cannot be used in identification or validation processes. They were nevertheless included for the sake of completion.
As observed on the graphs, the behaviour of the laboratory stand during free-wheeling is typical of systems with high dry friction component. The experiments proved that the guideway was homogeneous and therefore that the resistive function R was independent form displacement y. Following [21] the LuGre dry friction model in the form of was incorporated as a resistive force model [22]. Here σ 0 , σ 1 , σ 2 represent, respectively, the stiffness, microdamping and viscous aspect of the resistive force and w is the internal friction state [the symbol w was chosen instead of more widely used z to avoid confusion with parameters in Eqs. (33)-(39)]. The value of w can be obtained by solving a differential equation in the form oḟ where g(ẏ) is a function of velocity used to incorporate the Stribeck effect and Coulomb friction into the model. It is defined as The values of F S and F C represent the values of stiction force and Coulomb friction force, while parameter v s determines how quickly does the value of function switch from F S to F C . According to references [23][24][25], the values of α should be in the range between 0.5 and 2. In this paper, α = 2 was chosen and used for both identification and validation.
The parameters in Eq. (43) were identified using the following free-wheeling model for the laboratory stanḋ where m is the forcer's mass and m A is the known value of the mass added during free-wheeling. As with all other models presented in this paper, the parameters were found using genetic algorithms. Table 1 shows the values of those parameters, and Fig. 11 shows sample validation of the free-wheeling model.

Sine model
By incorporating the obtained resistive force function in Eq. (31) and adding the shift parameter in the form of α S , the sine model is obtained. It is governed by equationṡ During identification, the excitation functions I 1 (t), I 2 (t) and I 3 (t) were constructed by setting the amplifier work mode to "Current sine function generator" and measuring actual values of current (as shown in Fig. 12) along with values of displacement and velocity. These values were then converted to a continuous function using the spline interpolation [26]. Figure 13 shows all three excitation functions. The model was then identified in similar manner as with the resistive force function by using a genetic algorithm. The comparison between best acquired simulation result and experimental data during identification is displayed in Fig. 15. Table 2 presents identified parameters value, and Fig. 14 shows examples of comparison between simulation and experimental data for four model cases of motor operation. It has to be noted that parameters of resistive force [Eq. (40)] were identified separately for the sine model and their values slightly differ from those listed in Table 1.
As can be seen in Figs. 14 and 15, the model with identified parameters fits closely to experimental data for most cases of closed-loop working conditions and with slowly changing sine current functions. However, during sudden and rapid changes of current value, as was the case with square wave current excitation function, the model exhibited strong oscillations of velocity and displacement not present in the physical counterpart.
The identified force constant parameter has negative value, which means that the phase order in mathematical model was inverted in comparison with the actual motor. It should also be noted that despite a slight change of resistive function parameters the obtained displacement and velocity functions for free-wheeling operations still very closely resemble the experiment (Fig. 16).

Elliptic theta model
The second tested model was identified in a very similar pattern. Firstly, a shift parameter z s and resistive force were incorporated into model's ODE (Eq. (39)) resulting in the following final model's equationṡ During identification, the current functions I j (t) were chosen and defined in the same manner as in Sect. 5.2 and are represented in Figs. 12 and 13. As in previous cases, genetic algorithm was used during the process. Figure 14 presents the obtained results, while values of model parameters are shown in Table 3. Same as with sine model, Fig. 17 presents examples of comparison between four experiment and simulation data for four cases of motor operation (experiment data and excitation functions are exactly the same as with sine model presented in Fig. 15).
The simulation results provide close resemblance to experiments for all studied cases of closed and openloop working conditions. In case of square current excitation, both the physical motor and simulation exhibit similar square wave oscillations; however, the experiment shows a slight drift of displacement values in negative direction, while the simulation continues oscillates around a single point. One possible explanation of this is that the stand tends to exhibit slightly lower friction for high negative values of velocity. In case of sine wave position excitation, the model exhibited an unexpected drop and peak in velocity during first moments of operation; it was, however, brief and did not greatly affect the desired position function shape and values. As with sine model, the dry friction parameters were changed compared to previous model; however, this did not adversely affect the model during free-wheeling operation. The negative value of z L demonstrates that, as with sine model, the phases for motor were chosen inversely.

Conclusions
The proposed model of coreless linear three-phase motor was constructed with respect to Lorentz force and electromagnetic induction. The functions of magnetic induction were assumed as a sum of infinitely many single magnets, while the magnets individual fields distribution was calculated using Amperes model of permanent magnets and Biot-Savart law. The numerically calculated distribution between two bar magnets was then superseded with a single Gaussian function. Two methods of calculating the sum from many magnets were proposed. First of them uses the sine approximation, and second provides exact values and incorporates Jacobi elliptic theta functions. For certain cases, with appropriate relationship between magnets size, distances between magnets and number of magnets (as shown by the red line in Fig. 6) the two solutions provide nearly identical values of field distribution. In accordance with these two methods, two mathematical models of the motor were created. Preceding their identification and validation, a resistive force function was chosen as a dry friction force and model with LuGre method. Both models were correctly identified and validated. In most studied cases, the results of simulations matched closely to experiment for both models; however, during operations with sudden and abrupt changes to current value the sine model produced significantly worse results. On the other hand, while the elliptic theta function model provided a more accurate results, the sine model's results were smoother and resembled the experiment better during sine wave position excitation. The value of theta function nome q 0 calculated during identification suggests that motor was indeed constructed in a manner in which sine approximation nets similar results as actual magnetic field values.
A huge differences between models were observed in both computation time and ease of finding the necessary parameters. On average, a single simulation of the model with elliptic theta field distribution took ten to twenty times longer than sine model. At the same time, finding satisfactory convergent solution took over ten thousand passes of genetic algorithm (so-called generations) for sine model, while elliptic theta model required only two thousand.
Both models can be used depending on required task. The elliptic theta model requires more computation time, but provides more accurate solution for all cases. The sine model works well for usual motor applications. It provides results faster and requires less parameters to identify, but falls short during atypical operations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.