Calibration of a Hub Dynamic Balancing Machine Based on the Least Squares Method and Systematic Error Analysis

To resolve issues such as excessive residual vibrations and unsatisfactory balance effects in the balancing process, the particle swarm optimization (PSO)algorithm is combined with the least squares influence coefficient method of rotor dynamic balance to perform dynamic balance calibration based on the research of the least squares influence coefficient method of wheel dynamic balance. The influence coefficient generally has a large error due to the influence of the vibration measured value, thereby lowering the accuracy of the calibrated influence coefficient. Therefore, the maximum likelihood estimate (MLE) method is employed to address the influence coefficient error, and the result is compared with the calibration value of the influence coefficient (IC) method. The analysis results indicate that the residual value generated by the calibration of the influence coefficient through the maximum likelihood estimate (MLE) is 1.036 while the residual value obtained through the influence coefficient (IC) method is 1.513, suggesting that the former exhibits a smaller systematic error and is closer to the true value.


I. INTRODUCTION
The calibration of the dynamic balance of a hub is a reverse mathematical mapping process that estimates the unbalance of the rotor based on the measured vibration response of the support structure. The calibration model and model parameter accuracy are the key factors affecting the accuracy of dynamic balance measurements. In addition, unbalanced vibrations of the rotors are the main excitation sources in rotary machines. Currently, the main methods used to balance rotors are the modal balance method and the influence coefficient method.
Currently, the industrial use of high-speed rotating machines has increased. Therefore, the characteristics and performance of these rotating systems need to be improved [1]. Wheel hub dynamic balancing machines are The associate editor coordinating the review of this manuscript and approving it for publication was Inês Domingues . often used to realize the online monitoring of the dynamic balance quality of high-speed rotating wheel hubs. The monitored items mainly include the upper and lower unbalance, static and even unbalance, and phase angle. The calibration of the dynamic balance of the wheel hub is an inverse mathematical mapping process that estimates the imbalance of the rotor based on the measured vibration response of the support structure. Due to the high requirements for the equipment's detection accuracy, repeatability, stability, and detection cycle, the calibration of the wheel hub, that is, the accuracy of the calibration model and model parameters, is a key factor affecting the accuracy of dynamic balance measurements [2], [3]. The unbalanced vibrations of rotors are the main sources of excitation in rotating machinery. There are two methods for wheel calibration: the modal balance method and the influence coefficient method. The modal balance method takes the structural parameters of the rotor support system, the mass moment of inertia of the rotor, the front position, and the measured speed as model parameters, and the entire calculation process requires a long period of time and the assembly tolerance of the mechanical structure has high requirements. The accuracy of the model also needs to be verified [4], [5]. The influence coefficient method is widely used in the field of dynamic balance because the principle is simple and does not need to consider the complex information of the rotor system, and it is also easily assisted by computers. In practical applications, residual vibrations exist. The least squares influence coefficient method is proposed to solve unbalanced equations. However, the residual vibrations of some measuring points obtained by the least squares that influence the coefficient method are comparatively large. Thus, the method needs to be further optimized [6]. LUND and TONNESON extended this method by taking into account the measurement error when optimizing the influence coefficient matrix. However, in the current literature, field balancing and active balancing are the primary focuses of influence coefficient research, while few studies have focused on the influence coefficient estimation of balancing machines utilized on motor assembly lines [7], [8]. The two-plane influence coefficient method acts as the mainstream calibration method for dynamic balancing measurements and material-removing integrated machines for rigid rotors [9]. SW DYER proposed the application of a single-plane influence coefficient method to rotating machinery without considering the two-plane influence on the influence coefficient. Studies that considered unbalance did not involve even unbalance, and the research was relatively simple [10]. For influence coefficient estimation, GOOD-MAN et al. adopted the least squares method to improve the accuracy of the influence coefficient method [11], [12]. Considering different tire widths, YANG combined permanent calibration and the complex coefficient least squares influence coefficient method, which requires the mechanical structure to have high accuracy and does not regard the system as a black box. The system error has a great effect on unbalance measurement solutions [13], [14]. XU designed a computer program for the least squares iterative calculation of the influence coefficient method [15]. H TAPLAK proposed the application of a genetic algorithm to rotating machinery, although this method has not been applied to actual equipment [16], [17]. TIWARI proposed the application of the least squares method of residual unbalances in flexible rotors, which decay with time, and studied the residual imbalance of the rigid rotor dynamic balancing machine of a wheel hub dynamic balancing machine [19]. ZHANG proposed a Monte Carlo-based evaluation method for the dynamic uncertainty of a dynamic balance measurement system, which pushed the accuracy of the static state to the dynamic evaluation. However, although the determination of dynamic uncertainty improves the influence coefficient algorithm, it is difficult to realize in the running equipment [20].
Research Contributions of This Work: The main contribution of this work are summarized as follows.
The paper presents the influence of measurement errors on the quality and accuracy of wheel dynamic balance and discusses the methods of multiple linear regression to compare the residual values obtained by the maximum likelihood estimate (MLE) method and the influence coefficient (IC) method.
The paper discuss the influence of residual vibration, and describe residual vibration as the objective function, and use the particle swarm optimization (PSO) algorithm to optimize the residual vibration of each measurement point.
Finally, the paper will extend this method to calibrate wheels of various sizes in the following work.
The rest of this article is arranged as follows. The theoretical background of the proposed method is introduced in Section 2. The experimental verification of the proposed method is presented in Sections 3 and 4. Finally, the conclusions are drawn in Section 5.

II. CALIBRATION OF THE INFLUENCE COEFFICIENT
The measured value of the vibration response of the support structure is applied in the calibration process to estimate the influence coefficient instead of the true value. Therefore, the measurement error of the measured value of the vibration response is the main factor affecting the calibration accuracy. If the vibration response data containing the measurement errors are directly used according to the influence coefficient model in the ideal environment, the correct result will not be obtained when solving the unbalance U = f −1 (V ), thus making it necessary to address the measurement error and eliminate its influence on the previous test model. The measurement data obtained by repeated measurements are a sample of the measured values of the vibration response, and the main purpose of processing the measured data in the calibration process is to infer the overall information from the sample information of the measured values of the vibration response.
In the IC calibration process, the numbers of correction planes and vibration points are equal; the equations have a unique solution that does not consider the effect of system errors on the measured value. In the least squares method, one set of correction masses (the solution to the equations) that minimizes the quadratic sum of the residual vibration values at each equilibrium rotation speed of each vibration point is obtained. According to the principle of the IC method, the ideal relationship between the IC and the unbalance and original vibration is expressed by Eq. (1).
The IC method assumes that the rotor possesses linear characteristics. The linear relationship between the rotor unbalance and the vibration, which is caused by the mass unbalance, is represented by an IC. For two balance planes A and B, assuming that the rotor to be balanced has m test points in total and that the vibration value of the m measurement point is v AM AA0 after the test weight is added to the A plane, VOLUME 8, 2020 the IC is expressed by Eq. (2).
The double-plane IC method for the dynamic balance of the rigid rotor hub is applied to determine the unbalance on the condition of two balance planes, two vibration sensors, and a single-speed as well as the IC of a certain point. The unbalance of the rotor is distributed on two balance planes, the amplitude of the vibration response is obtained from the sensors installed on two balance planes, and the phase of the unbalance is estimated through a reference point.
Assuming that v 1 and v 2 are defined as the predicted values of the periodic vibration signals of the balance planes A and B, then the mathematical relationship between the sum of the unbalance and the amount of unbalance scattered on the balance planes A and B is expressed by Eq. (3): where A denotes the IC matrix and its elements are all complex variables. In matrix A, a ij represents the IC and refers to the influence of the unbalance of the planej on the vibration of the plane i; i = A, B; j = A, B.
Since the amplitude and phase of the vibration response measured by independent sensors are not correlated, separate estimations are allowed for the average and standard deviation of the amplitude and phase. Assuming that with the m (m = 1, 2 . . . l) test weight on the balance plane, k measurements of the vibration response of the support structure with sensor i(i = A, B) have been performed and the measurement error is ε AM , it can be concluded from the IC format that the amplitudes of the vibration response and unbalance have a linear relationship. Therefore, the linear regression method and the MLE method can be used to estimate the parameters. Taking a AM AA as an example for analysis, the following parameters can be set in Eq. (4): Then, in Eq. (5): The maximum likelihood estimation is given by Eq. (6): where L is the likelihood function.
Therefore, the maximum likelihood estimation of the IC a AM AA can be obtained with Eq. (7): The phase of the IC is the difference between the phase of the vibration response and the phase of the unbalance. Then, the value of IC is given by Eq. (8):

III. LEAST SQUARES METHOD BASED ON PARTICLE SWARM OPTIMIZATION
The least squares method (LSM) is a mathematical optimization technique that can find the best function match for the data through the minimized quadratic sum of errors. Goodman introduced the least squares method to the IC method [21]. The principle is to minimize the quadratic sum of the residual vibrations at each vibration point to determine the counterweight U . However, the least squares method is not the optimal solution in practical applications because it focuses on solving large residual vibrations. In addition, some large counterweight values cannot be easily achieved in practical counterweight installation. Therefore, the PSO algorithm has been applied to optimize the least squares IC method, obtaining a counterweight solution more in line with the actual requirements.

A. DESIGN OF THE LEAST SQUARES INFLUENCE COEFFICIENT METHOD
According to the principle of the least squares method, δ is the equation system for the residual vibration of each vibration point at this time; the matrix is used to represent Eq. (10) after simplification. A T is the conjugate transposed matrix of A, and the counterweight U can be obtained by solving the equation.
In actual operation, the amount of counterweight born by the high-speed rotor is limited. Since the least squares method emphasizes the residual vibration while ignoring the constraints of the counterweight range, a larger balance mass is the result; therefore, constraints are required. Additionally, the residual vibration cannot be too large in order to obtain satisfactory results. The mass of the counterweight is restricted to the required range [0, q lim ], and the quadratic sum of the residual vibrations in the least squares method is set as the fitness function, which is expressed in Eq. (11): The weighting factor needs to be repeatedly modified until the counterweight is less than u lim ; then, the residual vibration can be calculated. If the constraint conditions are met, then the calculation is finished; otherwise, the above steps need to be repeated.

B. PARTICLE SWARM OPTIMIZATION (PSO)
ZHU uses a genetic algorithm to determine the flexible rotor balance, but the genetic algorithm requires encoding and decoding and cross-mutation, and the solution parameters are more complicated [11]. WAN uses the IC method based on a genetic algorithm to address the unbalance of flexible rotors. However, there are no experimental demonstrations, and the application in actual engineering is more complicated [18]. WANG applies the least squares IC method of PSO to the theoretical analysis and experimental verification of the dynamic balance of flexible rotors [25]. Compared with the genetic algorithm, the PSO algorithm has a faster search speed, higher efficiency, simple algorithm implementation, and fast optimization speed and can more easily jump out of local extremes. PSO is an intelligent optimization method that simulates the foraging process of a flock of birds. This method has certain advantages because it only requires limited mathematical analysis, is independent of gradient information, and presents rapid convergence. Considering the industrial production environment, timeliness, and economy, PSO is widely used in industry [9], [10].
PSO algorithms are based on swarms and fitness. Each particle represents a potential solution, which is characterized by three parameters: position, velocity, and fitness. The velocity and position of the particles are updated according to their own best position (pBest) and the whole group's best position (gBest). The merit of the particle is measured by the fitness and represented by the objective function value (Fig. 1).
The dimension is D = 10, the number of iterations is M = 150, and the initial population is N = 10. In addition, the initial position and velocity of each particle are randomly generated, and the particles continue at the given speed in the environmental space to find the best update. The acceleration coefficients C1 = 2, C2 = 4, and r 1 and r 2 are random numbers in the range of [0,1]. The formula is provided Eq. (13).
The dynamic balance should achieve a satisfactory balance effect with the smallest possible balance quality. The least squares IC method may cause the balance quality to be too large for an application because the least squares method emphasizes that a smaller residual vibration does not impose restrictions on the range of the counterweight. Therefore, when the counterweight meets the requirements, the residual vibration must not be too large in order to obtain satisfactory results. To ensure that all particles are within the feasible range, the particle speed and position need to be restricted before each iteration. If the particle exceeds the defined range, then the extreme value of the range is taken. Therefore, the boundary range is set as follows: yp(i, j) > X max(j)|yp(i, j) < X min(j) yp(i, j) = rand * (X max(j)−X min(j))+X min(j) Although the standard PSO algorithm has simple principles and is easy to program, this method has inherent shortcomings. This article improves the standard PSO algorithm for specific situations. As the number of iterations increases, the particles move away from each other. The algorithm approaches the optimal solution, and the inertia coefficient of the particles should be gradually reduced; otherwise, the algorithm readily oscillates near the optimal solution. For this reason, this article applies an attenuation process to the inertia coefficient, ω = ω 0 (1−0.9 i−1 n−1 ), where ω 0 is the initial inertia coefficient and n is the total number of iterations. Due to the limited length of computer bytes, the numerical precision is limited. In the search process, there may be situations where multiple different particle positions correspond to the global optimal solution, and the unknown mean value of the particles corresponds to the optimal fitness of each generation VOLUME 8, 2020 in the program, which is determined by absorbing the position information of each corresponding particle.
The operation steps are described as follows: (1) Set the PSO population size to S, which is the number of balancing weights; (2) Initialize the velocity and position of each particle; (3) Calculate the fitness of each particle according to Eq. (14); (4) Compare the fitness of each particle, and record pBest; (5) Update the velocity and position of the particles; (6) Calculate each choice probability according to the fitness; (7) Apply the crossover and mutation operations, and generate new populations; and (8) Check whether the result meets the end condition.
If not, go to step (4); otherwise, output the optimal solution.

C. SOLUTION FOR THE SYSTEMATIC ERROR AND CALIBRATION OF THE INFLUENCE COEFFICIENT
According to the characteristics of system error and random error, multiple linear regression theory is used for the fitting relationship between the measurement error and the existing IC model to establish a multiple linear regression-based IC calibration model and simultaneously achieve accurate calibration of the IC and acquire a solution for the system errors. In addition, an updated IC model with error processing capabilities is constructed to effectively deal with measurement errors, thereby improving the accuracy of the dynamic balance test system and the IC calibration. After the influence of abnormal values is eliminated, the measured value of the vibration response is composed of only the true vibration value, systematic error, and random error. The establishment of the IC model can be obtained as shown in Eq. (16): where V 0 is a fixed value that represents the system error in the collection of the current set of data and ε denotes the random error following the normal distribution. According to the measured values of V , the IC can be calculated through parameter estimation, and a new set of values can be estimated through the recalibration and the collection of the next set of data. The system error can be corrected by the updated IC model. Additionally, the values of a and σ can be used as the basis for measurement performance evaluation and error source analysis. Assuming that the random error of the vibration response of the support structure on balance planes A and B is the variance in ε A and ε B , then V and a can be estimated through the least squares method: (18) where V denotes the systematic errors. The measured value can be corrected in the subsequent balancing process to bring the measured value closer to the true value. The value of σ 2 estimated by the matrix method is as follows:

IV. DISCUSSION
The structure of the dynamic tester is presented in Fig. 2; this tester is used for the inspection of wheels before they leave the factory.Equipment designed and developed according to China Academy of Machinery Science and Technology Group. The equipment brand is the MASS, and the model of the hub dynamic balancing equipment is AM32.
The hardware system of the tester includes the choice of the dynamic balancing spindle end motor (Beckhoff AM8553-1K20-0000), driver (Beckhoff AX5112-0000-0200), reducer (AE090-008), hard support installation rotor device and two force sensors that detect rotor vibration (installed at the fixed end). In addition, the dynamic balance of the hub is automatically measured. The rotor used in these experiments is in the motor of a wheel hub (illustrated in Fig. 2). The permissible unbalance of the rotor on one plane is 80 g, and the balancing rotating speed is 1800 r/min. An already balanced wheel hub is selected as the experimental target, and multiples of the permissible unbalance on one plane are chosen as the trial weights. Partially collected data are presented in Fig. 3. The wheel hub diameter is 22 inches; the centre hole diameter is 550-130 mm; and the efficiency is 18 seconds/piece, namely, 10 g, 20 g, 30 g, 40 g, 50 g, 60 g, 70 g, and 80 g. Each trial weight is positioned on only one balancing plane at the angular positions of 0 • and  180 • , and the vibration responses of the rotor are recorded. One group of experimental data is presented in Fig. 3.
In actual production, multiple additions of counterweights (i.e., multiple balances) are required to meet the standards during the whole dynamic balance process. Therefore, the effectiveness of the balancing method can be verified by the performance of a single balance. The sum of the added counterweights in the whole process of reaching the required dynamic balance is set as the true value of the unbalance of the main shaft. If the vibration response data containing the measurement errors are directly used according to the IC model in the ideal environment, then the correct result is not obtained when determining the unbalance (u = f (v)). Thus, it is necessary to address the measurement error and eliminate its influence in the previous test model. Moreover, different trial weights on the A plane and B plane were selected for repetitions of the above experiment, and the results are in line with the above analysis. One group of experimental data is presented in Table 1.
The calibration steps are described as follows.
(1) Drive the hub without any trial weight on the plane and measure the original vibration;    Tables 2 and 3) is an essential criterion for evaluating the rotor dynamic balance. The IC method and the MLE-based IC method are separately used to acquire coefficient estimates, and the effectiveness of the new method is verified through comparisons of the IC values obtained by the two methods.  First, Eqs. (1) and (2) can be used to determine the true value of the IC at different phase angles after a counterweight has been added. Second, the residual values obtained through the MLE-based IC method (Eqs. (4) and (7)) and the PSObased least squares method are 1.513 and 1.036, respectively, indicating that the MLE-based IC method exhibits a smaller system residual. The closer the calibrated value of the IC is to the true value, the closer the residual vibration of the balanced main shaft system is to 0. If the dynamic characteristics of the main shaft system change or the calibrated value of the IC is inaccurate, the system quality and balance accuracy decrease significantly and dynamic balance failure may even occur in several cases. The mean square value of the residual vibration with the rotor dynamic balance obtained by the MLE-based IC method is smaller; additionally, the balance effect is better and has practical value.
The calibration of the IC shows that the measured value of the vibration response of the support structure is adopted in the calibration process to estimate the IC rather than the true value. Therefore, the measurement error of the measured value of the vibration response is the main factor influencing the calibration accuracy. Considering that measurement errors are inevitable when the test system measures the vibration response, the measured value of the vibration response of the support structure is applied in the calibration process to estimate the impact instead of the true value, which causes a certain deviation in the measurement results and cannot accurately reflect the true value of the measured value. Therefore, the influence of the measurement error must be considered in the process of system calibration. Comparison of the calibration results obtained by the two methods and multiple linear regression analysis of the residual error indicates that the residual value calibrated through the MLE method is 0.011 and the value of the IC method is 0.02. Although both values are less than 0.05 and within the standard range, the residual error of the MLE method is smaller and has less system error (Figs. 6 and 7).

V. CONCLUSION
In PSO, in addition to considering the minimum residual vibration, the weight must be constrained (because the least squares method increases the weight of the distribution at certain moments and causes gross errors, the distribution weight needs to be restricted). PSO has the characteristics of selforganization, self-learning and self-adaptation and requires less coding. After the objective function and fitness function are determined, the information self-organizing search ability is obtained.
The amount of residual vibration imbalance is an important criterion for evaluating rotor dynamic balance. Experiments have proven that the residual vibration mean square value obtained after the rotor dynamic balance estimated based on the MLE influence coefficient is smaller, an improved balance effect can be achieved. In-depth analyses of the impact of the main errors of the dynamic balancing machine measurement system on calibration are still required.
With the error characteristics as the basis, maximum likelihood estimation theory is introduced based on the correctability of the system error and the normal distribution of the random error, and the IC calibration model, which can synchronously solve system errors, is proposed to realize the IC in the calibration process. Finally, the use of multiple linear regression analysis to compare the calibration results of the two methods indicates that the proposed method yields a much smaller residual mean square error than the traditional method and that the unbalanced phase change is small. YUAN CAO is currently pursuing the Ph.D. degree with the University of Science and Technology Beijing, Beijing, China. From 2018 to 2020, she designs precision instruments with the Machinery Science and Technology Group. From 2015 to 2018, she was a Teacher with the College. She is mainly engaged in the research of artificial intelligence, such as particle swarm optimization algorithm (PSO), genetic algorithm optimization (GA), neural network algorithm, and data analysis. TAO WANG is currently pursuing the master's degree with the School of Mechanical Engineering, University of Science and Technology Beijing. His tutor is Prof. Jianguo Cao. During his master's degree, he has participated in many topics of the tutor, mainly studying the mathematical model and process optimization of strip rolling, involving nuclear power alloy full-process integrated shape control technology. He is engaged in some artificial intelligence optimization, such as particle swarm optimization algorithm (PSO), genetic algorithm optimization (GA), and neural network algorithm.