Hybrid Multi-Domain Analytical and Data-Driven Modeling for Feed Systems in Machine Tools

: Position error-compensation control in the servo system of computerized numerical control (CNC) machine tools relies on accurate prediction of dynamic tracking errors of the machine tool feed system. In this paper, in order to accurately predict dynamic tracking errors, a hybrid modeling method is proposed and a dynamic model of the ball screw feed system is developed. Firstly, according to the law of conservation of energy, a complete multi-domain system analytical model of a ball screw feed system was established based on energy ﬂow. In order to overcome the uncertainties of the analytical model, then the data-driven model based on the back propagation (BP) neural network was established and trained using experimental data. Finally, the data-driven model was coupled with the multi-domain analytical model and the hybrid model was developed. The model was veriﬁed by experiment at di ﬀ erent velocities and the results show that the prediction accuracy of the hybrid model reaches high levels. The hybrid modeling method combines the advantages of analytical modeling and data-driven modeling methods, and can signiﬁcantly improve the feed system’s modeling accuracy. The research results of this paper are of great signiﬁcance to improve the compensation control accuracy of CNC machine tools.


Introduction
Feed systems are used to position the machine tool components carrying the cutting tool and workpiece to the desired location [1]. They are one of the most important subsystems of computer numerical controlled (CNC) machine tools, as their positioning accuracy and speed determine the quality and productivity of machine tools [2]. Feed drive systems are either powered by linear motors directly, or by rotary motors via a ball screw and nut. Among them, ball screw feed drive systems are widely used in CNC machine tools for their high stiffness, reliable operation and ability to mitigate the impact of inertial and cutting force variations [3,4]. Machining accuracy significantly depends on the tracking performance of the feed systems for a given trajectory. Therefore, the ball screw should exhibit good transient and steady tracking performance to meet the demands of high productivity and high precision of machine tools [5,6].
To improve the tracking accuracy of the ball screw feed system, position error compensation control is commonly used. If the tracking error can be estimated accurately in advance, the control problem of the feed system can be eased. Therefore, a model of the ball screw feed system that can accurately predicted tracking error is essential to improve tracking performance. There are many factors affecting the tracking performance of the feed system, such as nonlinear friction, backlash, vibration, estimating the energy consumption of a milling machine, while Gaussian process (GP) has been used to develop the data-driven energy prediction model with uncertainty for a milling machine [36]. Jinkyoo Park et al. [42] have discussed the data-driven modeling methods that are based on Bayesian statistic inference and illustrate their potential applications for performance characterization, condition diagnostic and control optimization of machine tool systems, for which analytical modeling of such systems can be difficult. Data-driven models are widely used and have achieved good results, but do not have interpretability and expressiveness. They need training data, but the data can hardly cover all working conditions. So, the fully data-based models may suffer from inaccuracies as well when the working condition of the feed system varies in a large range. In order to overcome the shortcomings of purely data-driven modeling methods, a hybrid analytical and data-driven modeling method was proposed. Relatively little related literature has appeared so far. René Felix Reinhart et al. [43,44] have used the hybrid analytical and data-driven modeling method to establish the kinematics model for soft robots and the inverse dynamic model for rigid robots. The authors only consider the hybrid modeling of robots at the component level.
This paper argues that an analytical model of a feed system should be established as a whole, incorporating mechanical, electrical, control, and other subsystems. Considered at system level, this paper develops a hybrid multi-domain analytical and data-driven model for the ball screw feed system to predict the tracking error. According to the law of conservation of energy, based on energy flow and symmetry transformation, a complete multi-domain system analytical model of the ball screw feed system is established. In order to compensate for the uncertainty of the model structure and parameters, a data-driven model based on a neural network was coupled to the multi-domain analytical model in parallel. The modeling method can implement seamless integration of complex systems in multiple domains and improve modeling accuracy. The hybrid model of the ball screw feed system was validated by experiment results.
The remainder of the paper is organized as follows. Section 2 introduces the multi-domain analytical modeling process for the ball screw feed system. Section 3 introduces the data-driven model. Section 4 introduces the hybrid modeling method for the ball screw feed system. Section 5 introduces the experimental verification of the hybrid model. Section 6 summarizes the full text.

Multi-Domain Modeling Method for the Ball Screw Feed System
It is generally known that a ball screw feed system is a multi-domain integrated system which incorporates mechanical, control, electrical and other subsystems. This paper argues that the multi-domain coupling characteristics of the feed system should be taken into account when establishing an analytical model. Therefore, a modular non-causal modeling method based on Modelica language is used to establish the multi-domain analytical model of the feed system. Modelica is an object-oriented physical system modeling language, which supports the modeling of component models based on equations and the modeling of complex system based on component non-causal connection. In Modelica, the interface of the component model is called the connector, and the coupling relationship established on the component connector is called the connection. If the connection expresses a causal coupling relationship, it is called a causal connection. If the connection expresses a non-causal coupling relationship, it is called a non-causal connection. The component model connection mechanism is shown in Figure 1.
The key to establishing a multi-domain knowledge-driven model of ball screw feed system is how to divide the system and define the interfaces of component models in each domain. In this paper, the ball screw feed system is divided according to its real physical structure, then component models are established and its interfaces are defined, and finally the multi-domain integrated model of feed system is established by connecting component models. The ball screw feed system studied in this Symmetry 2019, 11, 1156 4 of 20 paper was driven by permanent magnet synchronous motor, and the module division of the system is shown in Figure 2. Each module can also be further divided into smaller modules. an analytical model. Therefore, a modular non-causal modeling method based on Modelica language is used to establish the multi-domain analytical model of the feed system. Modelica is an objectoriented physical system modeling language, which supports the modeling of component models based on equations and the modeling of complex system based on component non-causal connection. In Modelica, the interface of the component model is called the connector, and the coupling relationship established on the component connector is called the connection. If the connection expresses a causal coupling relationship, it is called a causal connection. If the connection expresses a non-causal coupling relationship, it is called a non-causal connection. The component model connection mechanism is shown in Figure 1.  The key to establishing a multi-domain knowledge-driven model of ball screw feed system is how to divide the system and define the interfaces of component models in each domain. In this paper, the ball screw feed system is divided according to its real physical structure, then component models are established and its interfaces are defined, and finally the multi-domain integrated model of feed system is established by connecting component models. The ball screw feed system studied in this paper was driven by permanent magnet synchronous motor, and the module division of the system is shown in Figure 2. Each module can also be further divided into smaller modules. In order that each component model can be connected, the uniform interface in each domain should be defined. In non-causal connection, the interface in each domain usually contains two types of variables: flow variables and potential variables. When the interface is connected, it satisfies the generalized Kirchhoff theorem, that is, the sum of flow variables is zero and the potential variables are equal. In causal connection, the connector includes input variables and output variables. In this paper, mainly the mechanical, control and electrical fields are considered, in which the connection in the mechanical and electrical fields is a non-causal connection, and the connection in the control field is a causal connection. The connectors defined are shown in Table 1: The multi-domain analytical model of the ball screw feed system was developed base on Modelica language for its advantages and capabilities of component-based non-causal modeling. And the model mainly included the mechanical subsystem model, servo control subsystem model, the permanent magnet synchronous motor (PMSM) model, sensor models and other models. The structure of the model is shown in Figure 3. In order that each component model can be connected, the uniform interface in each domain should be defined. In non-causal connection, the interface in each domain usually contains two types of variables: flow variables and potential variables. When the interface is connected, it satisfies the generalized Kirchhoff theorem, that is, the sum of flow variables is zero and the potential variables are equal. In causal connection, the connector includes input variables and output variables. In this paper, mainly the mechanical, control and electrical fields are considered, in which the connection in the mechanical and electrical fields is a non-causal connection, and the connection in the control field is a causal connection. The connectors defined are shown in Table 1: The key to establishing a multi-domain knowledge-driven model of ball screw feed system is how to divide the system and define the interfaces of component models in each domain. In this paper, the ball screw feed system is divided according to its real physical structure, then component models are established and its interfaces are defined, and finally the multi-domain integrated model of feed system is established by connecting component models. The ball screw feed system studied in this paper was driven by permanent magnet synchronous motor, and the module division of the system is shown in Figure 2. Each module can also be further divided into smaller modules. In order that each component model can be connected, the uniform interface in each domain should be defined. In non-causal connection, the interface in each domain usually contains two types of variables: flow variables and potential variables. When the interface is connected, it satisfies the generalized Kirchhoff theorem, that is, the sum of flow variables is zero and the potential variables are equal. In causal connection, the connector includes input variables and output variables. In this paper, mainly the mechanical, control and electrical fields are considered, in which the connection in the mechanical and electrical fields is a non-causal connection, and the connection in the control field is a causal connection. The connectors defined are shown in Table 1: The multi-domain analytical model of the ball screw feed system was developed base on Modelica language for its advantages and capabilities of component-based non-causal modeling. And the model mainly included the mechanical subsystem model, servo control subsystem model, the permanent magnet synchronous motor (PMSM) model, sensor models and other models. The structure of the model is shown in Figure 3. The key to establishing a multi-domain knowledge-driven model of ball screw feed system is how to divide the system and define the interfaces of component models in each domain. In this paper, the ball screw feed system is divided according to its real physical structure, then component models are established and its interfaces are defined, and finally the multi-domain integrated model of feed system is established by connecting component models. The ball screw feed system studied in this paper was driven by permanent magnet synchronous motor, and the module division of the system is shown in Figure 2. Each module can also be further divided into smaller modules. In order that each component model can be connected, the uniform interface in each domain should be defined. In non-causal connection, the interface in each domain usually contains two types of variables: flow variables and potential variables. When the interface is connected, it satisfies the generalized Kirchhoff theorem, that is, the sum of flow variables is zero and the potential variables are equal. In causal connection, the connector includes input variables and output variables. In this paper, mainly the mechanical, control and electrical fields are considered, in which the connection in the mechanical and electrical fields is a non-causal connection, and the connection in the control field is a causal connection. The connectors defined are shown in Table 1: The multi-domain analytical model of the ball screw feed system was developed base on Modelica language for its advantages and capabilities of component-based non-causal modeling. And the model mainly included the mechanical subsystem model, servo control subsystem model, the permanent magnet synchronous motor (PMSM) model, sensor models and other models. The structure of the model is shown in Figure 3.

Flow variable
Torque τ Flow variable Current i Input variable u The key to establishing a multi-domain knowledge-driven model of ball screw feed system is how to divide the system and define the interfaces of component models in each domain. In this paper, the ball screw feed system is divided according to its real physical structure, then component models are established and its interfaces are defined, and finally the multi-domain integrated model of feed system is established by connecting component models. The ball screw feed system studied in this paper was driven by permanent magnet synchronous motor, and the module division of the system is shown in Figure 2. Each module can also be further divided into smaller modules. In order that each component model can be connected, the uniform interface in each domain should be defined. In non-causal connection, the interface in each domain usually contains two types of variables: flow variables and potential variables. When the interface is connected, it satisfies the generalized Kirchhoff theorem, that is, the sum of flow variables is zero and the potential variables are equal. In causal connection, the connector includes input variables and output variables. In this paper, mainly the mechanical, control and electrical fields are considered, in which the connection in the mechanical and electrical fields is a non-causal connection, and the connection in the control field is a causal connection. The connectors defined are shown in Table 1: The multi-domain analytical model of the ball screw feed system was developed base on Modelica language for its advantages and capabilities of component-based non-causal modeling. And the model mainly included the mechanical subsystem model, servo control subsystem model, the permanent magnet synchronous motor (PMSM) model, sensor models and other models. The structure of the model is shown in Figure 3.

Potential variable
Angle ϕ

Potential variable
Voltage v Output variable y The key to establishing a multi-domain knowledge-driven model of ball screw feed system is how to divide the system and define the interfaces of component models in each domain. In this paper, the ball screw feed system is divided according to its real physical structure, then component models are established and its interfaces are defined, and finally the multi-domain integrated model of feed system is established by connecting component models. The ball screw feed system studied in this paper was driven by permanent magnet synchronous motor, and the module division of the system is shown in Figure 2. Each module can also be further divided into smaller modules. In order that each component model can be connected, the uniform interface in each domain should be defined. In non-causal connection, the interface in each domain usually contains two types of variables: flow variables and potential variables. When the interface is connected, it satisfies the generalized Kirchhoff theorem, that is, the sum of flow variables is zero and the potential variables are equal. In causal connection, the connector includes input variables and output variables. In this paper, mainly the mechanical, control and electrical fields are considered, in which the connection in the mechanical and electrical fields is a non-causal connection, and the connection in the control field is a causal connection. The connectors defined are shown in Table 1: Table 1. Connector defined in each subject area.

Mechanical Connectors Electrical Connectors Control Connectors Icon
Icon The multi-domain analytical model of the ball screw feed system was developed base on Modelica language for its advantages and capabilities of component-based non-causal modeling. And the model mainly included the mechanical subsystem model, servo control subsystem model, the permanent magnet synchronous motor (PMSM) model, sensor models and other models. The structure of the model is shown in Figure 3.
The multi-domain analytical model of the ball screw feed system was developed base on Modelica language for its advantages and capabilities of component-based non-causal modeling. And the model mainly included the mechanical subsystem model, servo control subsystem model, the permanent magnet synchronous motor (PMSM) model, sensor models and other models. The structure of the model is shown in Figure 3.

Modeling and Analysis of the Mechanical Subsystem
The mechanical part of a ball screw feed system is usually composed of sub-components, such as couplings, ball screws, supporting bearings at both ends, linear guide rails, screw nuts, the workbench, and so on. When establishing the mechanical subsystem model, the mechanical connector interface should be defined first, and then the component models are established based on the equation. The modeling schematic diagram is shown in Figure 4. The output torque of servo motors is represented as and the rotational and translational motion of mechanical parts are both driven by servo motors. The torsional stiffness of the coupling, the torsional stiffness and axial stiffness of the ball screw, the axial stiffness of the bearing and the screw nut, the backlash, and the nonlinear friction are all considered in the model.  For the servo motor, according to Euler's law of motion, the torque balance equation on its output shaft is where is the rotor inertia of the motor, is the angular displacement of the output shaft of the motor, is the viscous damping coefficient and is the load torque on the motor shaft. Considering the flexibility of each mechanical component, then where Θ is the angular displacement of the screw shaft. is the total torsional stiffness converted from the stiffness of each component.
can be calculated by the following equation, found in [27]:

Modeling and Analysis of the Mechanical Subsystem
The mechanical part of a ball screw feed system is usually composed of sub-components, such as couplings, ball screws, supporting bearings at both ends, linear guide rails, screw nuts, the workbench, and so on. When establishing the mechanical subsystem model, the mechanical connector interface should be defined first, and then the component models are established based on the equation. The modeling schematic diagram is shown in Figure 4. The output torque of servo motors is represented as T e and the rotational and translational motion of mechanical parts are both driven by servo motors. The torsional stiffness of the coupling, the torsional stiffness and axial stiffness of the ball screw, the axial stiffness of the bearing and the screw nut, the backlash, and the nonlinear friction are all considered in the model.

Modeling and Analysis of the Mechanical Subsystem
The mechanical part of a ball screw feed system is usually composed of sub-components, such as couplings, ball screws, supporting bearings at both ends, linear guide rails, screw nuts, the workbench, and so on. When establishing the mechanical subsystem model, the mechanical connector interface should be defined first, and then the component models are established based on the equation. The modeling schematic diagram is shown in Figure 4. The output torque of servo motors is represented as and the rotational and translational motion of mechanical parts are both driven by servo motors. The torsional stiffness of the coupling, the torsional stiffness and axial stiffness of the ball screw, the axial stiffness of the bearing and the screw nut, the backlash, and the nonlinear friction are all considered in the model. For the servo motor, according to Euler's law of motion, the torque balance equation on its output shaft is where is the rotor inertia of the motor, is the angular displacement of the output shaft of the motor, is the viscous damping coefficient and is the load torque on the motor shaft. Considering the flexibility of each mechanical component, then where Θ is the angular displacement of the screw shaft. is the total torsional stiffness converted from the stiffness of each component.
can be calculated by the following equation, found in [27]: For the servo motor, according to Euler's law of motion, the torque balance equation on its output shaft is where J m is the rotor inertia of the motor, θ m is the angular displacement of the output shaft of the motor, B m is the viscous damping coefficient and T L is the load torque on the motor shaft.
Considering the flexibility of each mechanical component, then where Θ L is the angular displacement of the screw shaft. K L is the total torsional stiffness converted from the stiffness of each component. K L can be calculated by the following equation, found in [27]: where K T is the total equivalent torsional stiffness of the feed system, K Z is the total equivalent axial tension compression stiffness of feed system, and η is efficiency of the driving mechanism. K T and K Z can be calculated as follows [27]: In the above equation, K tc is torsional stiffness of coupling, K ts is the torsional stiffness of the lead screw, K as is the axial stiffness of the lead screw, K an is the axial stiffness of the nut, K ab1 and K ab2 are the axial stiffness support bearing. The formulas for calculating the above stiffness parameters can be found in previous studies [10,27]: where E bs is the elastic modulus of the lead screw material, A s is the sectional area of the lead screw, L sg is the distance between two supporting bearings, x is the distance between the nut and one end of the bearing, F p is the nut preload level, K n0 is the contact stiffness of nut, and C n0 is double nut pre-loading dynamic load rating. For the lead screw, it is obvious that the moment balance equation on the shaft is where J c and J s are the moment of inertia of the coupling and the screw shaft respectively, T t is the load torque of the screw shaft. There is backlash between the ball screw and the nut. The backlash has non-linear characteristics, which has a great influence on the positioning accuracy of the feed drive system. In this paper, a nonlinear elastic method [13] is used to establish the backlash model. The characteristics of the backlash are described by an abstract spring in the method. When the screw is in the gap zone, the spring stiffness is zero, and when the screw passes through the gap zone, the spring becomes stiff. The torque-angular displacement curve is shown in Figure 5.
where is the total equivalent torsional stiffness of the feed system, is the total equivalent axial tension compression stiffness of feed system, and is efficiency of the driving mechanism. and can be calculated as follows [27]: In the above equation, is torsional stiffness of coupling, is the torsional stiffness of the lead screw, is the axial stiffness of the lead screw, is the axial stiffness of the nut, and are the axial stiffness support bearing. The formulas for calculating the above stiffness parameters can be found in previous studies [10,27]: where is the elastic modulus of the lead screw material, is the sectional area of the lead screw, is the distance between two supporting bearings, x is the distance between the nut and one end of the bearing, is the nut preload level, is the contact stiffness of nut, and is double nut preloading dynamic load rating.
For the lead screw, it is obvious that the moment balance equation on the shaft is where and are the moment of inertia of the coupling and the screw shaft respectively, is the load torque of the screw shaft.
There is backlash between the ball screw and the nut. The backlash has non-linear characteristics, which has a great influence on the positioning accuracy of the feed drive system. In this paper, a nonlinear elastic method [13] is used to establish the backlash model. The characteristics of the backlash are described by an abstract spring in the method. When the screw is in the gap zone, the spring stiffness is zero, and when the screw passes through the gap zone, the spring becomes stiff. The torque-angular displacement curve is shown in Figure 5. For the ball screw feed system, when the backlash is considered, the load torque on the screw shaft can be calculated as follows:  For the ball screw feed system, when the backlash is considered, the load torque on the screw shaft can be calculated as follows: In Equation (9), K c and d are the equivalent spring constant and equivalent damping of backlash respectively, b is the total backlash, and ∆θ = θ L − θ t . θ t is related to the displacement of the workbench. Regardless of vibration, it can be seen in [18] that the displacement of the workbench S is: S is the additional displacement of the worktable caused by the axial and torsional deformation of the lead screw, and where F is the driving force of the workbench, which is generated by the torque T t , For the workbench, according to Newton's theorem, the force equilibrium equation is where F x is the axis component of cutting force. According to Stribeck model [10], the friction force F f can be calculated as follows: where the µ is the Coulomb friction coefficient, M is load mass, µ v viscous friction coefficient, f s is the Stribeck effect, b 0 represents the exponential decay. By synthesizing the above formulas, the dynamic model of the mechanical part of the ball screw feed system can be established. The model is encapsulated and the interface for data exchange with other models is set up. The encapsulated model of the mechanical part can be seen in Figure 4.

Modeling of the Servo Drive System
The feed system studied in this paper was powered by a permanent magnet synchronous motor, and the space vector control is used in servo drive system, which includes position control loop, speed control loop and current control loop. The motor was fed by a pulse width modulation (PWM) voltage source thyristor inverter, and the voltage modulation mode was space vector pulse width modulation (SVPWM) [45]. The servo drive system mainly included proportional integral (PI) controllers, coordinate transformation modules, the SVPWM module and inverter. In this paper, based on the analysis of the theory of space vector control, the modular modeling method described above is applied to establish the servo drive system model. First of all, the module models containing input and output interfaces are established. Digital PI controllers are usually used in actual servo control systems. In order to keep the model consistent with reality, the model of digital PI controllers is established based on discrete control theory. The theory of digital PI control algorithm is as follows: The coordinate transformation module includes Park transformation, Clark transformation and corresponding inverse transformation. These module models are established according to the transformation equation. The details of the Park and Clark transformation can be seen in [21].
The implementation of the SVPWM algorithm in vector control can be summarized as consisting of three steps: determination of the sector, calculation of the basic space vectors duration and calculation of the PWM comparison threshold, as shown in Figure 6.
The coordinate transformation module includes Park transformation, Clark transformation and corresponding inverse transformation. These module models are established according to the transformation equation. The details of the Park and Clark transformation can be seen in [21].
The implementation of the SVPWM algorithm in vector control can be summarized as consisting of three steps: determination of the sector, calculation of the basic space vectors duration and calculation of the PWM comparison threshold, as shown in Figure 6. According to the above analysis process, each module model was established based on the Modelica language, and then the SVPWM module model was obtained by connection. The outputs of the SVPWM model are the control signal of the inverter switch tube. The topology circuit of the inverter is shown in Figure 6: VT1~VT6 are power switch tubes, and the on-off of each switch is controlled by the drive signal S1~S6, 、 and is the output voltage of the inverter. In this paper, the electrical interface model, power switch transistor and diode model are established first, and then the inverter model is established by connecting the component models according to the topology circuit structure of the inverter. Similar to the model of mechanical parts, the servo control subsystem model was encapsulated and several interfaces set up, as shown in Figure 6.

Modeling of the Permanent Magnet Synchronous Motor
The electromagnetic relations of the PMSM are very complicated. In order to simplify the analysis, some assumptions are made as follows: (1) The conductivity of permanent magnet material is zero; (2) There is no damping winding on the rotor; (3) Stator winding current produces only sine distribution of magnetic potential in the air gap, ignoring the high-order harmonic of magnetic field; (4) During the steady-state operation, the waveform of induction electromotive force in phase winding is sinusoidal.
The physical quantities of AC motor windings, such as voltage, current and flux, vary with time, and time phasors are often used to represent them in analysis. However, they can also be defined as spatial vectors if considering the space position of the windings. The synthetic vector formed by the According to the above analysis process, each module model was established based on the Modelica language, and then the SVPWM module model was obtained by connection. The outputs of the SVPWM model are the control signal of the inverter switch tube. The topology circuit of the inverter is shown in Figure 6: VT1~VT6 are power switch tubes, and the on-off of each switch is controlled by the drive signal S1~S6, U A , U B and U C is the output voltage of the inverter. In this paper, the electrical interface model, power switch transistor and diode model are established first, and then the inverter model is established by connecting the component models according to the topology circuit structure of the inverter. Similar to the model of mechanical parts, the servo control subsystem model was encapsulated and several interfaces set up, as shown in Figure 6.

Modeling of the Permanent Magnet Synchronous Motor
The electromagnetic relations of the PMSM are very complicated. In order to simplify the analysis, some assumptions are made as follows: (1) The conductivity of permanent magnet material is zero; (2) There is no damping winding on the rotor; (3) Stator winding current produces only sine distribution of magnetic potential in the air gap, ignoring the high-order harmonic of magnetic field; (4) During the steady-state operation, the waveform of induction electromotive force in phase winding is sinusoidal.
The physical quantities of AC motor windings, such as voltage, current and flux, vary with time, and time phasors are often used to represent them in analysis. However, they can also be defined as spatial vectors if considering the space position of the windings. The synthetic vector formed by the same physical quantity in three-phase windings is the synthetic space vector of this physical quantity. The symmetrical three-phase current, voltage and flux linkage in three-phase symmetrical stator windings of permanent magnet synchronous motor can be regarded as the space vector. Through coordinate transformation of the space vector, the mathematical equation of permanent magnet synchronous motor in a two-phase rotating coordinate system can be derived [45].
The voltage equation of d and q axis is shown as follows: where ω is the electrical rotor angular speed, R s is equivalent stator resistance, u d and u q are the d-and-q-axis voltage, i d and i q are the d-and-q-axis current, ψ d and ψ q are the d-and-q-axis flux linkage. The flux linkage equations of d and q axes are as follows: where L sσ is the damper stray inductance, L md is the stator main field inductance per phase in d-axis, L mq is stator main field inductance per phase in q-axis, i f is the equivalent excitation current of the rotor permanent magnet. i f is a parameter and the calculation of it can be seen in [46]. The electromagnetic moment equation acting on the motor shaft is as follows: where T em is the electromagnetic torque and p is the number of pole pairs. Based on the above PMSM mathematical model, the equivalent circuit of PMSM in d and q axis coordinates can be obtained, as shown in Figure 7: The symmetrical three-phase current, voltage and flux linkage in three-phase symmetrical stator windings of permanent magnet synchronous motor can be regarded as the space vector. Through coordinate transformation of the space vector, the mathematical equation of permanent magnet synchronous motor in a two-phase rotating coordinate system can be derived [45]. The voltage equation of d and q axis is shown as follows: where is the electrical rotor angular speed, is equivalent stator resistance, and are the d-and-q-axis voltage, and are the d-and-q-axis current, and are the d-and-q-axis flux linkage.
The flux linkage equations of d and q axes are as follows: where is the damper stray inductance, is the stator main field inductance per phase in daxis, is stator main field inductance per phase in q-axis, is the equivalent excitation current of the rotor permanent magnet.
is a parameter and the calculation of it can be seen in [46]. The electromagnetic moment equation acting on the motor shaft is as follows: where is the electromagnetic torque and is the number of pole pairs. Based on the above PMSM mathematical model, the equivalent circuit of PMSM in d and q axis coordinates can be obtained, as shown in Figure 7: Based on the analysis and the equivalent circuit mentioned above, and considering the loss models of copper, iron and friction, a permanent magnet synchronous motor model is established in this paper. The encapsulated model of the PMSM can be seen in Figure 7. Completely accurate PMSM modelling is difficult to obtain because that PMSM is a high-order, non-linear, strong coupling system. So, some assumptions have been made to analyze the PMSM model. The error produced by these assumptions is hard to avoid, but the error can be reduced by modifying the model parameters through experiments. What's more, the error can be compensated for to a certain extent by the datadriven model. Based on the analysis and the equivalent circuit mentioned above, and considering the loss models of copper, iron and friction, a permanent magnet synchronous motor model is established in this paper. The encapsulated model of the PMSM can be seen in Figure 7. Completely accurate PMSM modelling is difficult to obtain because that PMSM is a high-order, non-linear, strong coupling system. So, some assumptions have been made to analyze the PMSM model. The error produced by these assumptions is hard to avoid, but the error can be reduced by modifying the model parameters through experiments. What's more, the error can be compensated for to a certain extent by the data-driven model.

The Learned Data-Driven Model
The accuracy of the analytical model can hardly reach a high level because the model can hardly capture all aspects of the feed system's intrinsic properties, and there remain unmodeled dynamics. Therefore, a data-driven model was developed as a supplement to the analytical model.
In this thesis, the data-driven model was used to fit the deviation between the output of the analytical model and the measured value. A back propagation neural network (BPNN) was chosen to build the data-driven model because the data-driven model is small in scale and the BPNN structure is simple, the local extremum points are few, the fitting accuracy is high, and it is easy to transplant. The structure of BPNNs can be divided into input layer, hidden layer and output layer. In the hybrid model, the data-driven model was used to predict the deviation between theoretical model output and measured value. The scale of the data-driven model was relatively small, so a three-layer BPNN was chosen in which the position command signal and the working table displacement output from analytical model were the input layer neurons. There was only one output neuron which was the deviation between the working table displacement measured by experiment and that simulated by the analytical model. The relevant theory of BPNN can be referred to [47]. Combining with the empirical formula, the number of neurons in the hidden layer was determined by the trial and error method. The number of neurons in the hidden layer is 10 and the empirical formula is shown as follows In Equation (21), m, n and l are the number of nodes in the input layer, output layer and hidden layer respectively, and a is a constant between 0 and 10.
BP neural network represents the mapping relationship between n independent variables and m dependent variables. Before making a prediction, network training is needed to give it the ability to perform association and decision-making. If ω ij represents the connection weights between input layer i and hidden layer j, v jk represents connection weights between hidden layer j and output layer k, a j and b k represent the thresholds of hidden layer j and output layer k, respectively, then the output of the hidden layer is: ω ij x i − a j ), j = 1, 2, . . . , 10, n = 2 (22) In Equation (22), f 1 is the activation function of the hidden layer. The sigmoid function shown in Equation (23) was chosen in this paper.
The output of the BP network is shown as follows: In Equation (24), f 2 is the activation function of the output layer and that the function f 2 (x) = x was chosen in this paper.
If Y 1 represents the actual value of the deviation between measured working table displacement and that simulated from the analytical model, then the error function of the BP network is shown as follows.
The weights and thresholds of the BPNN are updated by gradient descent method of error function. In order to train the BPNN, some input and output label data are collected by experiment, as shown in Figure 8. The input working table command displacement shown in Figure 8a is the sinusoidal signals with periods of 12π seconds in the first 30 s and periods of 2π seconds in the rest of the 30 s, and that the amplitude of the signal is 0.1 m. The second input of the BPNN shown in Figure 8b is the simulation working table displacement of the analytical model. The output label of the BPNN shown in Figure 8c is the deviation between measured working table displacement and that simulated from analytical model. Use of the steepest gradient descent method to train the BPNN data-driven model and the results are shown in Figure 9. In Figure 9, the Loss is the variance of the deviation between the label value and the output of the BP network model. It can be seen that after about 2000 training sessions, the loss value tended to be stable. The data-driven model can be established by using the DBN [34], ANN [37], GP [36] and so on. In this paper, the simple BP network was chosen to avoid over-fitting because of the datadriven model is small in scale. After a great quantity of training, the weights and thresholds of the data-driven model can be learned.  The input working table command displacement shown in Figure 8a is the sinusoidal signals with periods of 12π seconds in the first 30 s and periods of 2π seconds in the rest of the 30 s, and that the amplitude of the signal is 0.1 m. The second input of the BPNN shown in Figure 8b is the simulation working table displacement of the analytical model. The output label of the BPNN shown in Figure 8c is the deviation between measured working table displacement and that simulated from analytical model. Use of the steepest gradient descent method to train the BPNN data-driven model and the results are shown in Figure 9. The input working table command displacement shown in Figure 8a is the sinusoidal signals with periods of 12π seconds in the first 30 s and periods of 2π seconds in the rest of the 30 s, and that the amplitude of the signal is 0.1 m. The second input of the BPNN shown in Figure 8b is the simulation working table displacement of the analytical model. The output label of the BPNN shown in Figure 8c is the deviation between measured working table displacement and that simulated from analytical model. Use of the steepest gradient descent method to train the BPNN data-driven model and the results are shown in Figure 9. In Figure 9, the Loss is the variance of the deviation between the label value and the output of the BP network model. It can be seen that after about 2000 training sessions, the loss value tended to be stable. The data-driven model can be established by using the DBN [34], ANN [37], GP [36] and so on. In this paper, the simple BP network was chosen to avoid over-fitting because of the datadriven model is small in scale. After a great quantity of training, the weights and thresholds of the data-driven model can be learned. In Figure 9, the Loss is the variance of the deviation between the label value and the output of the BP network model. It can be seen that after about 2000 training sessions, the loss value tended to be Symmetry 2019, 11,1156 12 of 20 stable. The data-driven model can be established by using the DBN [34], ANN [37], GP [36] and so on. In this paper, the simple BP network was chosen to avoid over-fitting because of the data-driven model is small in scale. After a great quantity of training, the weights and thresholds of the data-driven model can be learned.

Hybrid Model of the Ball Screw Feed System
The multi-domain hybrid modeling method proposed in this paper combines the advantages of analytical and data-driven modeling methods. The principle of this hybrid modeling method is shown in Figure 10.
In the formula, ∈ and ∈ are input and output vectors, is a function for the hybrid model relation between and , and and are the functions associated with analytical model and data-driven model. ∈ is the parameter vector of the function , and and are the parameter vectors associated with the function and , respectively. The symbol "⊕" represents a coupling operation between the analytical model and data-driven model. Based on experimental and simulation data and parameter identification algorithm, the parameters can be identified. The parameters can be obtained by the training algorithm. In general, a model that can learn from data without using any domain knowledge is called a data-driven model, for example, artificial neural networks, support vector machines, fuzzy methods, generalized linear models, and so on. In this thesis, the data-driven model is used to fit the deviation between the output of the analytical model and the measured value. The mathematical expressions of the coupling between knowledge-driven model and the data-driven model are given by the following: The data-driven model can be learned after training and then the data-driven model is coupled with the analytical model and the hybrid model of the ball screw feed system is established as shown in Figure 11. The hybrid model consists of two parts, which are the multi-domain analytical model and data-driven model. The model can be expressed by the following mathematical formula.
In the formula, x ∈ R m andŷ ∈ R n are input and output vectors, f is a function for the hybrid model relation between x andŷ, and f k and f d are the functions associated with analytical model and data-driven model. p ∈ R p is the parameter vector of the function f, and p k and p d are the parameter vectors associated with the function f k and f d , respectively. The symbol "⊕" represents a coupling operation between the analytical model and data-driven model. Based on experimental and simulation data and parameter identification algorithm, the parameters p k can be identified. The parameters p d can be obtained by the training algorithm.
In general, a model that can learn from data without using any domain knowledge is called a data-driven model, for example, artificial neural networks, support vector machines, fuzzy methods, generalized linear models, and so on. In this thesis, the data-driven model is used to fit the deviation between the output of the analytical model and the measured value. The mathematical expressions of the coupling between knowledge-driven model and the data-driven model are given by the following: The data-driven model can be learned after training and then the data-driven model is coupled with the analytical model and the hybrid model of the ball screw feed system is established as shown in Figure 11.

Layer1
Layer2 Layer3 Figure 11. The hybrid model of the ball screw feed system.

Simulation and Experiments
In order to validate the multi-domain hybrid system model, a ball screw feed system test platform was set up as shown in Figure 12. The position command signal of the working table was generated by the numerical control device. The angular displacement of the servo motor could be measured by the photoelectric encoder and the position of the workbench could be measured by the linear encoder. Given the different types of position command signals, the actual position, velocity and acceleration of the working table, and the output torque of the servo motor, and so on, could be measured by the corresponding experimental instruments.  The main parameters of the ball screw feed system are listed in Table 2. When these parameters were imported into the multi-domain hybrid model of feed system, the tracking error could be predicted by the model and the comparison between the measured results and the simulation results could be performed. Table 2. The parameters of the ball screw feed system.

Simulation and Experiments
In order to validate the multi-domain hybrid system model, a ball screw feed system test platform was set up as shown in Figure 12. The position command signal of the working table was generated by the numerical control device. The angular displacement of the servo motor could be measured by the photoelectric encoder and the position of the workbench could be measured by the linear encoder. Given the different types of position command signals, the actual position, velocity and acceleration of the working table, and the output torque of the servo motor, and so on, could be measured by the corresponding experimental instruments.

Layer1
Layer2 Layer3 Figure 11. The hybrid model of the ball screw feed system.

Simulation and Experiments
In order to validate the multi-domain hybrid system model, a ball screw feed system test platform was set up as shown in Figure 12. The position command signal of the working table was generated by the numerical control device. The angular displacement of the servo motor could be measured by the photoelectric encoder and the position of the workbench could be measured by the linear encoder. Given the different types of position command signals, the actual position, velocity and acceleration of the working table, and the output torque of the servo motor, and so on, could be measured by the corresponding experimental instruments. The main parameters of the ball screw feed system are listed in Table 2. When these parameters were imported into the multi-domain hybrid model of feed system, the tracking error could be predicted by the model and the comparison between the measured results and the simulation results could be performed.  The main parameters of the ball screw feed system are listed in Table 2. When these parameters were imported into the multi-domain hybrid model of feed system, the tracking error could be predicted by the model and the comparison between the measured results and the simulation results could be performed.

Single Axis Ball Screw Feed System Hybrid Model
First, the single-axis ball screw feed system hybrid model was verified under the sinusoidal position command signal (the given position command signal varies with time according to the sinusoidal law). The G code was loaded to make the displacement of the working table change with time according to sinusoidal law. The maximum feed speed was set to 1000 mm/min, 2000 mm/min, 3000 mm/min, 4000 mm/min, 5000 mm/min, and 6000 mm/min respectively, and the stroke was set to 100 mm. The given position command signals of the working table are shown in Figure 13.

Single Axis Ball Screw Feed System Hybrid Model
First, the single-axis ball screw feed system hybrid model was verified under the sinusoidal position command signal (the given position command signal varies with time according to the sinusoidal law). The G code was loaded to make the displacement of the working table change with time according to sinusoidal law. The maximum feed speed was set to 1000 mm/min, 2000 mm/min, 3000 mm/min, 4000 mm/min, 5000 mm/min, and 6000 mm/min respectively, and the stroke was set to 100 mm. The given position command signals of the working table are shown in Figure 13. The sinusoidal position command signal with the maximum feed speed of 1000mm/min and 6000 mm/min were used to train the data-driven model, and the rest of the signals were used to verify the hybrid model. The same G code was imported into the hybrid model and then the simulation results were compared with the measured results.
If represents the displacement of the workbench measured by the experiment, represents the displacement predicted by the pure analytical model, represents the displacement predicted by the hybrid model, and is the command displacement, then the dynamic tracking errors can be defined as follows: The sinusoidal position command signal with the maximum feed speed of 1000mm/min and 6000 mm/min were used to train the data-driven model, and the rest of the signals were used to verify the hybrid model. The same G code was imported into the hybrid model and then the simulation results were compared with the measured results.
If S m represents the displacement of the workbench measured by the experiment, S Ap represents the displacement predicted by the pure analytical model, S Hp represents the displacement predicted by the hybrid model, and S is the command displacement, then the dynamic tracking errors can be defined as follows: where E m is the tracking error measured by the experiment, E Ap is the tracking error predicted by the pure analytical model and E Hp is the tracking error predicted by the hybrid model. The hybrid model prediction error (E Hp − E m ) and the analytical model prediction error (E Ap − E m ) were compared, as shown in Figure 14.
where is the tracking error measured by the experiment, is the tracking error predicted by the pure analytical model and is the tracking error predicted by the hybrid model. The hybrid model prediction error ( − ) and the analytical model prediction error ( − ) were compared, as shown in Figure 14. It can be clearly seen in Table 3 that the prediction accuracy of the hybrid model was higher than that of the pure analytical model. To evaluate the difference between the simulation results and experiment results quantitatively, the maximum absolute error (MAE), root mean square error (RMSE) and relative error (RE, refers to the ratio of absolute error to measured true value, expressed as percentage) were evaluated. The simulation accuracy of the hybrid model and the analytical model can be seen in Table 3. The results show that the maximum deviation of the tracking error between the measured and predicted by the hybrid model were all under 6.79 μm. The RMSE of the prediction errors were all under 2.63 μm and the RE decreased gradually as the feed rate increased. The comparison of results showed that, compared with the pure analytical model, the hybrid model can improve the accuracy of tracking error prediction, and the accuracy of the hybrid model proposed in this paper can reach a relatively high level.  It can be clearly seen in Table 3 that the prediction accuracy of the hybrid model was higher than that of the pure analytical model. To evaluate the difference between the simulation results and experiment results quantitatively, the maximum absolute error (MAE), root mean square error (RMSE) and relative error (RE, refers to the ratio of absolute error to measured true value, expressed as percentage) were evaluated. The simulation accuracy of the hybrid model and the analytical model can be seen in Table 3. The results show that the maximum deviation of the tracking error between the measured and predicted by the hybrid model were all under 6.79 µm. The RMSE of the prediction errors were all under 2.63 µm and the RE decreased gradually as the feed rate increased. The comparison of results showed that, compared with the pure analytical model, the hybrid model can improve the accuracy of tracking error prediction, and the accuracy of the hybrid model proposed in this paper can reach a relatively high level.

Double Axis Ball Screw Feed System Hybrid Model
The double axis ball screw feed system consisted of the x-axis feeding system and the y-axis feeding system. The system model was established with the use of the hybrid multi-domain analytical and data-driven modeling method proposed in this paper, and it can be used to predict contour error. The accuracy of the model was verified by experiment under different speeds of circle drawing. The experimental velocities included 60 rad/min, 80 rad/min, 100 rad/min and 120 rad/min. The simulation and experiment results are shown in Figure 15. (The circle contour error is magnified 500 times in the figure).

Double Axis Ball Screw Feed System Hybrid Model
The double axis ball screw feed system consisted of the x-axis feeding system and the y-axis feeding system. The system model was established with the use of the hybrid multi-domain analytical and data-driven modeling method proposed in this paper, and it can be used to predict contour error. The accuracy of the model was verified by experiment under different speeds of circle drawing. The experimental velocities included 60 rad/min, 80 rad/min, 100 rad/min and 120 rad/min. The simulation and experiment results are shown in Figure 15. (The circle contour error is magnified 500 times in the figure). There are four circles in Figure 15, they are the standard circle, experiment circle, analytical model-predicted circle and the hybrid model-predicted circle. The standard circle is the trajectory command given by the numerical control device and the experiment circle is the trajectory measured by the linear encoder. The predicted circles of the pure analytical and the hybrid model are compared in the figure. As can be seen in the figure, the circle predicted by the hybrid model was much closer to the experimental circle than the circle predicted by the pure analytical model. Especially at high speed, the advantage of the hybrid model was more obvious. As the angular velocity increased, the prediction effect of the pure analytical model was getting worse. The contour error curves predicted by the pure analytical model and the hybrid model are show in Figure 16. It can be seen in the figure that the contour error predicted by the hybrid model was much smaller than that the analytical model predicted. Also, the contour error predicted by the analytical model increased as the angular velocity There are four circles in Figure 15, they are the standard circle, experiment circle, analytical model-predicted circle and the hybrid model-predicted circle. The standard circle is the trajectory command given by the numerical control device and the experiment circle is the trajectory measured by the linear encoder. The predicted circles of the pure analytical and the hybrid model are compared in the figure. As can be seen in the figure, the circle predicted by the hybrid model was much closer to the experimental circle than the circle predicted by the pure analytical model. Especially at high speed, the advantage of the hybrid model was more obvious. As the angular velocity increased, the prediction effect of the pure analytical model was getting worse. The contour error curves predicted by the pure analytical model and the hybrid model are show in Figure 16. It can be seen in the figure that the contour error predicted by the hybrid model was much smaller than that the analytical model predicted. Also, the contour error predicted by the analytical model increased as the angular velocity increased, whereas the contour error predicted by the hybrid model did not change much. Comparison of the maximum contour error predicted by the hybrid model and the pure analytical model are shown in Table 4. The results show that the maximum contour error predicted by the hybrid model can reach 6.96 μm and the prediction ability is better than the pure analytical model.

Conclusions
This paper proposed a multi-domain hybrid analytical and data-driven modeling method for the ball screw feed system in machine tools, and a hybrid model of a single-axis and double-axis ball screw feed system were developed by the proposed modeling method. The hybrid model comprised an approximate multi-domain analytical model and a learned, data-driven error model. Experiments were performed to validate the accuracy of the hybrid model and the results show that the model can significantly improve the feed system's modeling accuracy. The main conclusions of this article are as follows: 1. A hybrid multi-domain analytical and data-driven modeling method was proposed in this paper and a hybrid model of a ball screw feed drive system was established accurately using the hybrid modeling method. 2. In contrast to the traditional causal modeling method based on signal flow, the multi-domain Comparison of the maximum contour error predicted by the hybrid model and the pure analytical model are shown in Table 4. The results show that the maximum contour error predicted by the hybrid model can reach 6.96 µm and the prediction ability is better than the pure analytical model.

Conclusions
This paper proposed a multi-domain hybrid analytical and data-driven modeling method for the ball screw feed system in machine tools, and a hybrid model of a single-axis and double-axis ball screw feed system were developed by the proposed modeling method. The hybrid model comprised an approximate multi-domain analytical model and a learned, data-driven error model. Experiments were performed to validate the accuracy of the hybrid model and the results show that the model can significantly improve the feed system's modeling accuracy. The main conclusions of this article are as follows:

1.
A hybrid multi-domain analytical and data-driven modeling method was proposed in this paper and a hybrid model of a ball screw feed drive system was established accurately using the hybrid modeling method.

2.
In contrast to the traditional causal modeling method based on signal flow, the multi-domain integrated analytical model of the ball screw feed system was established using the non-causal modeling method based on energy flow. The analytical model of a feed system developed in this paper realized seamless integrated modeling of a complicated multi-domain system.

3.
A data-driven error model based on a BP neural network was established and the model was trained using experimental data. Then the learned data-driven error model was coupled with the analytical model of the ball screw feed system and the hybrid model was obtained. 4.
The hybrid model was validated using experimental data at different speeds, and the results show that, whether for the tracking error of a single-axis feeding system or the contour error of a double-axis feeding system, the prediction effect of the hybrid model is better than that of a pure analytical model, and the prediction accuracy of the hybrid model reaches a higher level.
The hybrid modeling method combined the advantages of the analytical modeling and the data-driven modeling methods. Complex system models often have numerous parameters, and some parameters are difficult to determine accurately. In addition, pure analytical models generally have assumptions, so pure analytical models inevitably have uncertainties and unmodeled dynamics. Coupling a data-driven model to the multi-domain integrated analytical model of feed systems, the uncertainty of the model can be reduced to a certain extent, and the accuracy of the model can be improved. The research results of this paper could be applied to error compensation control of CNC machine tools to improve their control accuracy.