Dynamic modeling and parameter identification for a gantry-type automated fiber placement machine

Dynamic models play a critical role in the design of model-based controllers, which has a significant effect on the dynamic characteristics of motion equipment. This paper mainly focuses on the dynamic modeling and parameter identification for a gantry-type automated fiber placement (AFP) machine. First, a dynamic modeling process combining prismatic axes and revolute axes is conducted by Newton-Euler method, in which the effects of friction and hydraulic balance system are also considered. Then, as the convenience for parameter identification and the application in linearity control, the methods of dynamic model linearization and determination of minimum inertial parameters based on the multi-body system (MBS) theory are proposed, with which a dynamic model in the form of linearized minimum inertial parameters is consequently established. To identify the parameters in the model, key issues regarding excitation trajectory, filtering, and identification algorithm are discussed in detail. Finally, corresponding experiments are performed on the AFP machine, and experimental results show that there is a good agreement between the prediction of the model and the measurement in actuality. Data analysis shows that except for Z-axis, the relative error rates of the others are not greater than 5%, which proves the effectiveness of the established dynamic model and the identified parameters.


Introduction
Advanced carbon fiber reinforced composite (CFRP) materials have been extensively used in the aerospace field due to their excellent performance over traditional metal materials, and their usage has reached 50% in advanced civil aircraft such as Airbus A350 and Boeing B787 [1,2]. Automatic fiber placement (AFP) machines are the key equipment for the industrial manufacturing of composite components, which mainly consists of a set of motion equipment and a self-contained fiber placement head. Fiber placement head is the core part of the AFP machine, which can independently clamp, refeeding and cut several  composite tows to form a tape, and then heat and compact them on mould surface simultaneously [1,3]. It is thus flexible and suitable for the placement of complex composite components. Since the movement of the fiber placement head is provided by the motion equipment, there is no doubt that the dynamic characteristics of the motion equipment directly determine the efficiency and quality of the composite component manufacturing.
With the continuous pursuit of the fiber placement efficiency, the speed of the AFP machine is also increasing, which has been reached 50 m/min. Such high speed requires the driving force of the motion axis can be able to change rapidly in a large range. However, limited by the performance of traditional controllers such as PID control, problems of low control accuracy and slow response speed are inevitable to be introduced, making the machine hardly meets the requirement of high-speed and high-precision. Due to this, advanced controllers based on dynamic models have attracted widespread attention and become a focus in control field [4][5][6][7]. As the basis of these advanced controllers, the dynamic model of the motion equipment has to be studied first.
Of all the dynamic modeling methods, the two most frequently used are the Newton-Euler method and the Lagrangian method. Among them, the Newton-Euler method is based on the idea of iteration, which has a clear physical meaning and high computational efficiency, and the Lagrangian method is based on the idea of energy balance, which avoids the derivation of internal force of joints and reduces the difficulty of modeling. But the dynamic equations established by these two methods are equivalent. Hwang [8] developed a recursive formulation for flexible dynamic manufacturing analysis of openloop systems with the Newton-Euler method. Lee and Shah [9] presented the dynamic analysis of a 3-degree-of-freedom (DOF) in-parallel actuated manipulator using the Lagrangian approach. These two methods are also widely used in the dynamic modeling of robotics [10][11][12][13]. Besides, Mahmoud et al. [14] pointed out that the dynamic model established by the Newton-Euler method was more appropriate to be linearized. Indeed, with the advantage of clear physical meaning, the dynamic equation established by the Newton-Euler method can be converted into a linear form of multiplication of inertia parameters and kinematic parameters through appropriate derivation, which is more suitable for the application in parameter identification and linearity control.
In addition, inertia parameters in the dynamic model of the linear form above can be further linearly combined to generate the minimum inertial parameters. This is an essential property of the dynamic equation, which has great advantages in simplifying the dynamic calculation, inertia parameter identification, control system design and analysis, etc. [15]. There are many attempts have been made to determine the minimum inertial parameters in robot fields. At the end of the twentieth century, Gautier and Khalil [16][17][18] and Gautier and Khalil [19] and Khalil and Bennis [20] proposed recursive formulations to directly calculate the minimum inertial parameters of robots with D-H parameters. As time goes, professional toolkits such as Symoro [21,22] and Pybotics [23] are developed for robots, which bring great convenience to researchers in the calculation. Unfortunately, the calculation method of the minimum inertia parameters is highly associated with the definition of joint coordinate systems. Different from the D-H method for robots, the joint coordinate systems of the AFP machine are usually established by the multi-body system (MBS) theory, leading to the recursive formulations in the above literature are no longer applicable. After the dynamic model in the form of linearized minimum inertial parameters is established, how to acquire accurate inertia parameters is an essential problem that has a great impact on the accuracy of the dynamic model. Due to the existence of manufacture errors, assembly errors as well as friction during the motion, the inertia parameters are generally not obtained from the CAD model directly, but by means of parameter identification [24][25][26]. A standard parameter identification procedure consists of modeling, experiment design, data acquisition, signal processing, parameter estimation, and model validation. Atkeson et al. [27] selected the fifth-order polynomial in joint space as the excitation trajectory and used the least-square method for parameter estimation; Swevers et al. [28] proposed an excitation trajectory based on a finite Fourier series, which was periodic and could average multiple sampled data to improve the signal-to-noise ratio; Olsen and Petersen [29] adopted the maximum likelihood method to estimate dynamic parameters with measurement noise successfully. In this paper, the improved finite Fourier series is selected as the excitation trajectory, and the weighted least-squares (WLS) method is adopted as the parameter estimation algorithm.
Dynamic modeling and parameter identification have been widely studied in robot fields, and the methods that involve dynamic model linearization and determination of minimum inertial parameters are mature. However, in the fields of traditional machine tools, whose kinematic and dynamic modeling is generally utilizing the MBS theory, because of their low operating speed, the dynamic models of the machine tools seem less important, hence these methods based on MBS theory are rarely mentioned in past literature. In contrast, as far as the AFP machine is concerned, the dynamic model is indispensable for the high performance under high speed and high load. Therefore, issues about the dynamic modeling and parameter identification for a gantry-type AFP machine are studied, and remains are organized as follows: Sect. 2 is devoted to the dynamic modeling process of the AFP machine and analyzes the methods of dynamic model linearization and determination of minimum inertial parameters; Sect. 3 investigates the excitation trajectory and the parameter identification method; Sect. 4 develops the corresponding experiments and discusses the results, and Sect. 5 concludes the paper and proposes the direction of future.

3
axes named X-, Y-, and Z-axis, and they are perpendicular to each other. The last three are revolute axes called B-, A-, and C-axis, and they rotate around Y-axis, X-axis, and Z-axis respectively. It is worth mentioning that the Z-axis adopts a hydraulic balance system to reduce the influence of gravity on the drive motors, and the B-axis is designed as an arc guide rail structure to reduce the distance between the A-axis rotation center and the tail end of the fiber placement head, thus improving the flexibility of the machine.
Several coordinate systems are established by MBS theory to describe the motion of each axis, which are the For efficient placement of large size composite components, three prismatic axes of the gantry-type AFP machine are designed with large strokes and high speeds. According to the design indexes, the motion parameters of each axis are shown in Table 1.

Dynamic modeling process based on the Newton-Euler method
The Newton-Euler method is adopted to establish the dynamic model of the AFP machine. In order to simplify the calculation, the motion of the C-axis is ignored during the modeling process for the following reasons: (a) The structure of the head is symmetrical, and the distribution of the mass is uniform; therefore, the offset torque generated by the C-axis motion is small. (b) Due to the material properties of composite tows, a small turning radius of the head in placing direction may cause defects such as wrinkles; therefore, the C-axis cannot rotate rapidly. In a summary, the dynamic model of the AFP machine is reduced to five axes. Based on the coordinate systems established in Fig. 1, a dynamic modeling process of the AFP machine combining prismatic axes and revolute axes is introduced. At first, the axes of the AFP machine are numbered 1, 2, 3, 4, and 5 for X, Y, Z, B, and A, and then, forward kinematic recursion and reverse dynamic recursion are conducted.
The forward kinematic recursion mainly calculates the angular velocity, angular acceleration, and linear acceleration of each axis, in this stage, the prismatic axis and the revolute axis need to be distinguished.
For prismatic axis: For revolute axis: where i+1 i+1 , i+1̇i +1 , and i+1̇i +1 are the 3 × 1 vector of angular velocity, angular acceleration, and linear acceleration of axis i + 1 in the coordinate system i + 1 respectively; (1) Fig. 1 The gantry-type AFP machine and its coordinate systems i+1 i R and i P i+1 are the 3 × 3 rotation matrix and the 3 × 1 position vector of the coordinate system i to the coordinate system i + 1 ; ḋ i+1 and d i+1 are the 3 × 1 vector of translation velocity and acceleration of the prismatic axis; ̇i +1 and ̈i +1 are the 3 × 1 vector of rotation velocity and acceleration of the revolute axis.
According to the above two equations, the linear acceleration at the center of mass of each axis can be deduced: where i+1 P C i+1 is the 3 × 1 position vector of the center of mass of axis i + 1 in the coordinate system i + 1.
By substituting i+1̇C i+1 , i+1 i+1 , and i+1̇i +1 into the Newton-Euler equation, the inertial force i+1 F C i+1 and inertia torque i+1 N C i+1 at the center of mass of each axis are obtained: where m i+1 and C i+1 I i+1 are the mass and the symmetric inertia matrix in the coordinate system of the center of mass of axis i + 1.
Finally, transfers of the force and the torque of each axis are calculated by reverse dynamic recursion: where i f i and i n i are the resultant force and resultant torque of axis i . The component of i f i in the direction of prismatic axis is the driving force i f i , and the component of i n i in the direction of revolute axis is the driving torque i n i . According to Eqs. (1)-(5), the dynamic modeling process of the gantry-type AFP machine combining prismatic axes and revolute axes is shown in Fig. 2, where n is the number of axes, and g is the acceleration of gravity.
From the modeling process in Fig. 2, a dynamic model of the AFP machine can be summarized as a second-order differential equation: is the mass matrix, C(q,q) is the centrifugal force and Coriolis force, G(q) is the gravity component, and vector.

Considering friction and hydraulic balance system
When the actual situations of the AFP machine are considered, the dynamic equation established above needs some modifications, to be specific, the effect of friction and the hydraulic balance system. Friction should be treated as an additional factor since it is always regarded as the main error of motion. Although friction is a complex nonlinear phenomenon, its linear approximation model including Coulomb and viscous friction can guarantee accuracy, which is extensively used and can be expressed as [14]: where f c is the coefficient of Coulomb friction, f v is the coefficient of viscous friction, ̇q is the velocity of axis, and sign is a symbolic function.
Therefore, the dynamic model considering friction is revised as follows: where F f (̇q) is the component of friction.
In addition to friction, the hydraulic balance system also has a great impact on the Z-axis. Since the total weight of the Z-axis, B-axis, A-axis, C-axis, and the fiber placement head exceeds 10 t, the Z-axis servo drive system will always maintain a large unilateral force if there is no gravity balance method taken, and this unilateral force severely limits the Z-axis motion performance and shortens the life of the Z-axis motor. This problem can be effectively alleviated by the use of the hydraulic balance system, and the driving force of the Z-axis can be expressed as: where f s z is the driving force provided by the servo drive system, f h z is the driving force provided by the hydraulic balance system, P is the hydraulic pressure which can be obtained from the pressure sensor, and A is the effective area of the hydraulic cylinder piston. So far, the establishment of the dynamic model is fully completed.
This dynamic model characterizes the relationship between the motion and the driving force or torque of the AFP machine. However, it is a nonlinear equation where the inertial parameters and the kinematic parameters are highly coupled, making it less convenient in parameter identification and application in linearity control.

Dynamic model linearization and determination of minimum inertial parameters
With the help of the parallel axis theorem, the symmetric inertia matrix of each axis in the coordinate system of the center of mass can be transformed into a representation in the joint coordinate system, by this way, the dynamic model can be rewritten as a linear form of inertia parameters: this is an important property of the dynamic equation, where Y is an observation matrix composed of the kinematic parameters, including the position, velocity, acceleration of each axis, and the offset lengths between coordinate systems; P is an inertia vector composed of the inertial parameters, including the symmetric inertia matrix, mass, position of the center of mass, and friction coefficients of each axis. Due to the complicated process, the specific derivation of the linearization for the dynamic equation based on the MBS theory is presented in the Appendix 1.
Unfortunately, not every inertial parameter in P has an impact on the dynamic model [14], in other word, there are column vectors of 0 or linearly dependent columns in the observation matrix Y , which make the matrix Y not fullrank. Therefore, the inertial vector P needs to be simplified through recombining some inertial parameters by certain linear relationships or using simple closed-form rules to eliminate some inertial parameters to obtain minimum inertial parameters.
For this reason, recursive formulations of the determination of minimum inertial parameters based on the MBS theory are proposed. Similarly, a detailed derivation process is shown in the Appendix 2 for concision of the body of the paper. Finally, the dynamic equation of the AFP machine in the form of linearized minimum inertia parameters is obtained: this equation is completely equivalent to Eqs. (9) and (11), but it is a linear form, where redundant items are removed by linear combination; therefore, it not only dramatically reduces the calculation of the equation, but also has a huge advantage in the identification of P m .
Each axis has 10 inertia parameters and two friction coefficients, and the AFP machine therefore has a total of 60 dynamic parameters in vector P . After calculating the minimum inertial parameters, there are only 23 parameters left in the minimum parameter vector P m , 13 of which are the combined inertial parameters, and the remaining 10 are the friction coefficients.
By combining Eq. (12) with N sampled data on a trajectory, a statically indeterminate linear equation can be constructed: where Y N m is the measurement matrix of 5N × 23 , and N d is the generalized driving force vector of 5N × 1 . Since this equation is linear, the parameters in the vector P m can be easily calculated as long as the values of the Y N m and the N d are collected from the experiment.

Excitation trajectory and parameter identification
As the dynamic equation in the form of linearized minimum inertia parameters is established above, the determination of the inertial parameter values becomes a key issue, which directly affects the accuracy of the dynamic model. In actual situations, due to the additional weight of cables, actuators, guide rails, sliding blocks, drive motors and guard bars, as well as the influence of manufacture errors, friction, etc., the accuracy of the dynamic model established by the theoretical inertia parameters must be poor. However, using the method of parameter identification to obtain a set of comprehensive parameters that meet the accuracy of dynamic calculation is one of the effective ways to improve the accuracy of the dynamic model. The main process of parameter identification is as follows: (a) design an excitation trajectory to drive the machine; (b) collect necessary data while the machine is moving; (c) process the collected data; (d) substitute the processed data into the dynamic model and identify the inertial parameters with parameter estimation algorithm; and (e) verify the identified parameters.

Excitation trajectory
The result of parameter identification depends on the quality of the excitation trajectory. In order to identify all the unknown parameters of the model, the excitation trajectory should be able to excite the whole system; otherwise, some parameters may be unable to be identified or sensitive to measurement noise.
There are many forms of excitation trajectory, such as fifth-order polynomial, finite Fourier series, et al. But these excitation trajectories cannot meet the requirement that the velocity and acceleration of axes are 0 at the beginning and ending of the trajectory, which result in poor trajectory tracking accuracy and even jitter when the machine follows these excitation trajectories, affecting the identification accuracy of dynamic parameters. In this paper, an improved finite Fourier series is used as the excitation trajectory of each axis: where t is time, f is fundamental frequency, and n is the harmonic item number of Fourier series, usually 5 is adopted, a i,k , b i,k , and c i,k are the harmonic or polynomial coefficients, and the velocity and acceleration can be calculated by differentiating this equation.
Define the boundary conditions of the trajectory to guarantee the velocity and acceleration at the start and stop of the trajectory are 0: where q i1 and q i2 are the positions of the start and stop of the trajectory, respectively, and the coefficients c i,k can be eliminated by solving the above six equations. The improved finite Fourier series excitation trajectory adopts a polynomial to replace the constant term of traditional finite Fourier series to realize the constraint of boundary conditions, thus reducing the start-stop shock and improving the validity of the measurement data of the entire trajectory.
Unmodeled errors and measurement noise will affect the accuracy of parameter identification, and their influence can be reduced by optimizing the excitation trajectory coefficients. According to the statically indeterminate linear Eq. (13), the condition number of the measurement matrix Y N m reflects the anti-noise ability of the identification method and the convergence rate of parameter estimation in identification process, the smaller the value, the smaller the influence of the disturbance on the identification result. Therefore, its condition number is defined as the objective function of the excitation trajectory coefficient optimization: where q min , q max , ̇q min , ̇q max , q min , and q max are the constraints of position, velocity, and acceleration, see Table 1 for specific values. To solve this multi-objective optimization problem with inequality constraints, genetic algorithm (GA), particle swarm optimization (PSO), or "fmincon" function in MAT-LAB are common methods. Considering the randomness and experience of parameter selection of GA and PSO, we selected "fmincon" function as the final choice to calculate the optimal coefficients of the excitation trajectory.

Measurement data processing
After the excitation trajectory is optimized, the AFP machine should track it repeatedly. At the same time, the relevant data including the position of each axis, the torque of each axis drive motor, and the hydraulic pressure are sampled at a constant, user-specified frequency. Subsequently, as the sampled data are often polluted by noise, they need to be processed accordingly to improve their validity. Since the improved finite Fourier series is used as the excitation trajectory in this paper, the signal-to-noise ratio of the sampled data can be improved by the time-domain averaging method first, so that the contingency of the identification results is reduced: where M is the times of repetition and k represents the k th sampling.
Then, in order to eliminate or reduce the influence of high-frequency random noise and improve the smoothness of curves, the sampled data should be further smoothed. Common smoothing methods include average method, spline function method, and five-point triple smoothing method. The average method is relatively simple and the filtering effect is poor. The spline function method adopts spline interpolation to approximate sampling point, so its algorithm is varied, and the effect is better, but its calculation is relatively complicated. The five-point triple smoothing method uses polynomial least square approximation to smooth the sampling points, which is simple and efficient, and it can be expressed as: where y = y 1 y 2 ⋯ y m is the original sampled data and ȳ = ȳ 1ȳ2 ⋯ȳ m is the filtered data. To get enough smooth data, this method could to be repeated several times.
In addition, velocity and acceleration of each axis are also needed in parameter estimation process, but if they are calculated by numerical differential of position data directly, high-frequency noise will be introduced. In order to avoid transfer error caused by differential, an analytical approach [28] is employed to estimate the velocity and acceleration of each axis. The premise of this approach is that each axis of the machine can accurately track the excitation trajectory, and specific steps are as follows: (a) fit the position data to Fourier series via least-square method and (b) differentiate the fitting Fourier series to get the velocity and acceleration. This approach can effectively decrease noise.

Parameter identification
With the processed data obtained above, unknown parameters in the dynamic model can be finally identified. Consequently, by substituting the processed measurement data into Eq. (13), an overdetermined linear equation is constructed: where P m is the minimum parameter vector that needs to be identified, and represents the measurement error noise with 0 mean normal distribution. The classic least-square method is the universal method in solving overdetermined systems of linear equations and as follows: but this method assumes that measurement torque is polluted by white Gaussian noise and all standard deviation of actuator torque noise is the same. This assumption cannot be met, resulting in less accuracy of estimation of parameters. However, due to its simplicity, the least-square method is still widely used in many applications.
The improvement of the least-square method is the WLS method, whose error is weighted by a weighted matrix determined by diagonal matrix of torque noise variance, thereby the accurate and inaccurate data can be distinguished: where ∑ is the weighted matrix. Therefore, the WLS method is selected as the parameter estimation method in this paper.

Model verification
The effectiveness of the identification method and the accuracy of the identified parameters must be demonstrated in convincing ways. Generally, corresponding verification experiments need to be performed, as illustrated in Fig. 3.
The absolute prediction error of the dynamic model can be quantified by the residual root mean square (RMS), which is expressed as follows: where n is sampling times, m (i) is the i th measured value of the verification trajectory, and d (i) is the i th predictive value of the dynamic model. In addition, the relative prediction error of the dynamic model can be described by the relative error rate: where m is a vector of measurement data. The RMS and the relative error rate represent the uncertainty of the prediction error, the smaller their value, the better the model prediction.

Experiment and discussion
To identify the inertial parameters, as well as verify the effectiveness of the established dynamic model and the parameter identification method, corresponding experiments were carried out on the self-developed gantry-type AFP machine, whose control system is developed based on the TwinCAT3 of Beckhoff company. The experimental platform is shown in Fig. 4, including a gantry-type AFP machine and an operating console.
The first one performed was the identification experiment, in which the basic frequency of the excitation trajectory was set to f = 0.08 Hz, and the motion period was thereby 25 s. After optimization calculation, an optimal excitation trajectory with the condition number of the measurement matrix of 505.8 was obtained. The excitation trajectory was discretized every 0.01 s, to create an executable file, and then imported into FIFO (First Input First Output) of the control system to drive the AFP machine. FIFO is a stack area with an interval of 0.01 s, and is used for multi-axis motion control of the AFP machine. At the same time, the position of each axis, the torque of each axis drive motor, and the hydraulic pressure were sampled in real-time with a sampling interval of 0.01 s. The AFP machine was operated to repeat the excitation trajectory 10 times. Took the average value of the 10 groups of sampled data, and the measurement trajectory of each axis was compared with the excitation trajectory, as shown in Fig. 5. It can be found from the figures that the measurement trajectory is highly coincident with the excitation trajectory, and the condition number of the measurement trajectory is 569.0, which is very close to that of the excitation trajectory, indicating that the axes of the AFP machine track the excitation trajectory well. Subsequently, the velocity and the acceleration of each axis were calculated by the analytical approach mentioned in Sect. 3.
For the output torque in the sampled data, in addition to the average processing, they must be filtered and then multiplied by the reduction ratio of each axis respectively. Moreover, the torque of the prismatic axes also needed to divide by the gear index circle radius to get the driving force, and the servo force of the Z-axis should additionally add the averaged hydraulic force. Finally, the driving force of the three prismatic axes and the driving torque of the two revolute axes before and after filtering were obtained, as shown in Fig. 6.
With all the necessary data obtained above, the minimum parameters of the AFP machine were identified, as shown in Table 2, where the theoretical values of combined inertial parameters were also listed to facilitate comparison.
From Table 2, we can see that the identification values of the inertia parameter are slightly larger than the theoretical values. Since the theoretical values are obtained from the CAD model, which ignores the quality of servo motors, electrical hardware, cables, accumulators, oil circuit, etc., the identification result can be considered credible. However, it is found that there are big errors between the identified values and the theoretical values of the A-axis, it is because we ignore the influence of the C-axis, and the combined inertia parameters of the A-axis thus consist of inertia parameters of the A-axis, C-axis, and fiber placement head. While the fiber placement head is highly integrated, including 32 fiber/paper servo motors and drivers, 16 groups of control mechanisms, etc., which cannot be calculated by the CAD model well, hence it is reasonable that the identified values are significantly greater than the theoretical values.
What is more, it also can be found that the friction coefficients of X-and Z-axis are significantly larger than those of other axes, whose coefficients of viscous friction all exceed 20,000 N·s/m, which implies these two axes are under great friction during the movement and may need some lubrication.
The identified parameters were then substituted into the Eq. (12), and a complete dynamic model of the gantry-type AFP machine was established. To further verify the accuracy of the identified parameters, the concept of cross-validation was adopted, and two verification trajectories different from the above-mentioned excitation trajectory were designed to drive the AFP machine. Each verification trajectory was repeated 5 times.
After the sampled data were processed accordingly, the measured values of the driving force and driving torque of each axis were obtained. At the same time, the predictive values of the driving force and driving torque of each axis were also calculated through the dynamic model and the identified parameters. Subsequently, the measured values and the predictive values were compared, the comparison of the verification trajectory I is shown in Fig. 7, and the comparison of the verification trajectory II is shown in Fig. 8. Figures 7 and 8 illustrate that the predictive curves are highly consistent with the measured curves, which efficiently proves the accuracy of the dynamic model and the identified parameters. However, when observing the difference curve of each figure, we can discover that the prediction errors are small but there are some jump values. The reason of this phenomenon is that the static continuous linear friction   Particularly, it should be pointed out that there are apparent errors at the corner between the predictive curves and the measured curves of the Z-axis (red circles in the figures). Because the driving force of the Z-axis is provided by both the hydraulic balance system and the servo drive system, the hydraulic system responds slowly while the servo system responds fast, the difference between them may have significant effects on the prediction accuracy of the Z-axis, introducing a huge error when the direction of the driving force switches.
Finally, the RMS and the relative error rate of the excitation trajectory and two verification trajectories were calculated, and the results are shown in Table 3. It can be seen that, except for Z-axis, the RMS of each axis in three trajectories is relatively small, and the relative error rates are only within 5%, which further verify the good performance of the identification result. Although the relative error rates of the Z-axis are relatively large, they still do not exceed 10%.
Through the analysis of the above figures and tables, it can be concluded that in spite of the friction model used in this paper does not reflect the actual situation well, and there is a mismatch between the response speed of the hydraulic system and the servo system, the established dynamic model and the identified dynamic parameters still have high accuracy.

Conclusions
Dynamic modeling and parameter identification of a gantrytype AFP machine are studied in this paper, and the research results can provide a model foundation for the model-based controllers. The main contribution of this paper is proposing recursive formulations of the dynamic model linearization and the determination of minimum inertial parameters based on the MBS theory. The methods studied in this paper can also provide guidance for dynamic modeling and inertial parameter identification for other high-speed machine tools.
Firstly, the Newton-Euler method is adopted to establish the dynamic model of the gantry-type AFP machine, in which friction and the hydraulic balance system are also considered. Secondly, based on the MBS theory, the recursive formulations of the dynamic model linearization and the determination of minimum inertial parameters are deduced, which bring great convenience for subsequent parameter identification. Thirdly, discuss the problems of parameter identification including excitation trajectory, filtering, parameter identification, and verification. Finally, experiments have been carried out on the gantrytype AFP machine, and the dynamic parameters are identified and verified successfully. The experimental results show that, although the friction model used in this paper cannot reflect the actual situation well, except for Z-axis, the prediction relative error rates of others still do not exceed 5%, which proves the accuracy of the dynamic model and the identified parameters.
Satisfactory results indeed are obtained in this paper; however, due to the lack of research about the hydraulic balance system, the problem that large prediction errors of Z-axis driving force has not been solved yet, which will be the study focus in the future.

Appendix 1. Dynamic equation linearization
The key to linearization of the dynamic equation is the separation of kinematic parameters and inertial parameters. Observing the dynamic modeling process, we found that these two parts of the parameters are coupled in the Newton-Euler equation and the reverse torque recursion process. Therefore, as long as these equations are linearized, the entire dynamic equation can be linearized.
Firstly, define 10 inertia parameters for each axis: where I xx i ∼ I zz i are the parameters in the symmetric inertia matrix expressed in the coordinate system of the center of mass of axis i ; m i r cx i ∼ m i r cz i are the product of the axis mass m i and the center of mass coordinate vector i P C i , m i is the mass of the axis i.
For Newton's equation: the cross product in the above formulation can be expressed in the form of the cross-product matrix: , and the linearization of Newton's equation is completed. For Euler's equation and reverse torque recursion equation, since the existence of i P C i × i F C i in the equations, they become nonlinear. However, with the help of the parallel axis theorem, the symmetric inertia matrix in the coordinate system of the center of mass can be transformed into a representation in the joint coordinate system {i} . By this way, two items in Euler's equation can be expressed: By combining all the inertial parameters in the dynamic equation into a matrix, the inertia matrix P in equation d = YP can be formed: (26) i F C i = S ̇i m i i P C i + S i i S i i m i i P C i + m i i̇i = 0 3×6 S ̇i + S i i S i i i̇i P i = H i P i (27) x P T y P T z P T b P T a T With the reverse dynamic recursion, the force of each axis is reversed: and the torque of each axis is reversed: finally, the observation matrix Y in equation d = YP is: The dynamic model established by the Newton-Euler method therefore is written as a linear form of inertia parameters: The friction model of each axis is also expressed in linear form: by adding Y f = sign(q)q to the observation matrix and P f = f c f v T to the inertia matrix, the linear form of the dynamic equation considering friction is obtained.

Appendix 2. Minimum inertial parameters calculation based on the MBS theory
From equation d = YP , it is obvious that the inertial parameters of axis i only affect 1 … i , rather than i+1 … n . Therefore, the main idea of calculation the minimum parameters of axis i is: (a) Let i+1 f i+1 = i+1 n i+1 = 0 0 0 T and calculate the identifiable inertial parameters and the unidentifiable inertial parameters of i . (b) Integrate the unidentifiable inertial parameters into the inertial parameter of axis i − 1 by the reverse dynamic recursion. The calculation process of the minimum parameters for the prismatic axis and the revolute axis is different. (31) n a = 0 3×10 0 3×10 0 3×10 0 3×10 A a P = Y n a P n b = 0 3×10 0 3×10 0 3×10 A b 0 3×10 + b a RY n a + S b P a b a RY f a P = Y n b P n z = 0 3×10 0 3×10 A z 0 3×10 0 3×10 + a z RY n b + S a P z a z RY f b P = Y n z P n y = 0 3×10 A y 0 3×10 0 3×10 0 3×10 + z y RY n z + S z P y z y RY f z P = Y n y P n x = A x 0 3×10 0 3×10 0 3×10 0 3×10 + y x RY n y + S y P x y x RY f y P = Y n x P