Industry-oriented method for dynamic force identification in peripheral milling based on FSC-LSQR using acceleration signals

Online measurement of milling force is very important for machining process monitoring and control. In practice, it is difficult to measure the milling force directly during the milling process. This paper develops a method for milling force identification called least square QR-factorization with the fast stopping criterion (FSC-LSQR) method, and the queue buffer structure (QBS) is employed for the online identification of milling force using acceleration signals. The convolution integral of milling force and acceleration signals is discretized, which turns the problem of milling force identification into a linear discrete ill-posed problem. The FSC-LSQR algorithm is adopted for milling force identification because of its high efficiency and accuracy, which can effectively handle the linear discrete ill-posed problem. The online identification of milling force can be realized using the acceleration signal enqueue and the milling force dequeue operations of the QBS. Finally, the effectiveness of the method is verified by milling tests under different milling parameters.


Introduction
Milling has been of great importance in the field of machining processes and the impact cannot be overemphasized. It is widely used in national defense, aviation, aerospace, and other fields. Since the goal of the milling process is to achieve higher productivity and better surface quality, process disturbances like chatter, collision, and tool wear need to be monitored and controlled [1,2]. All studies on these problems are closely related to milling forces. Therefore, becoming vital for the intelligent spindle to sense the milling force.
The piezoelectric dynamometer is the most commonly used instrument for online measurement of milling force [3,4]. Dynamometers for milling forces still have many limitations in industrial applications: (1) the dynamometer is too expensive; (2) the workpiece size is greatly limited; (3) inconvenient installation of the dynamometer. In order to overcome these shortcomings, Rizal et al. [5] mounted a force sensor between the spindle and the tool to measure the milling force, and verified its effectiveness through experiments. Kistler corporation has developed the wireless rotary milling force measurement instrument. The rotary dynamometer mentioned above is extremely expensive, complicated in structure, and changes the dynamic characteristics of the spindle, which means that the rotary dynamometer can just be used in laboratories rather than in the actual milling process. In summary, there is an urgent need for a low-cost, easy-to-install, fast, and accurate milling force measurement system.
In recent years, many experts have thought of using current, displacement, and other signals to identify milling forces. Machine tool feed motor contains a large amount of processing information. Jeong and Cho [6] and Aslan and Altintas [7] identified the milling force and extend its identification bandwidth using the current signal from the feed motor; however, this method requires accurate system parameters. Kim and Jeon [8], Chen et al. [9], and Song et al. [10] choose the spindle motor current signal with higher 1 3 static quasi-state sensitivity to monitor the milling force. The frequency bandwidth of the milling force identification using the current signal is relatively narrow [11]. Kim et al. [12] and Kim [13] proposed a method of using a cylindrical capacitive displacement sensor to quantitatively estimate the milling force. Albrecht et al. [14] used the Kalman filter to reduce the influence of process disturbance and measurement error on the identification of milling force by displacement signal. Multi-sensor fusion technology is used to improve the accuracy of milling force identification, such as displacement and acceleration signals (Zhou et al. [15] and Salehi et al. [16]), and current and acceleration signals (Hamid et al. [17]). Although the method above can identify the milling force during the high-speed milling process, the increase in cost and the decrease in reliability are inevitable.
Compared with the displacement sensor, the accelerometer has the advantages of easier installation and lower cost, which makes the acceleration signal more suitable for identifying the milling force during the milling process. Spiewak [18] first proposed a method to estimate the force indirectly from accelerometers mounted on the spindle box. In the actual industrial environment, it is difficult to measure the frequency response function of the spindle under rotating conditions, even if this modal characterization is more accurate. So the impact test is used to obtain the frequency response function of the spindle. Both the acceleration signal and the frequency response function can be measured. It seems that the milling force can be identified using the inverse filtering method [19]. However, Chae and Park [20] pointed out that the inverse filtering method is sensitive to noise at certain frequencies, and even small noises will have large errors. In actual milling conditions, the transfer matrix between acceleration signals and milling force is usually illconditioned and cannot be directly inverted to identify the milling force. This kind of problem is called a linear discrete ill-posed problem in the field of mathematics. Regularization is the main method for solving linear discrete ill-posed problems, including direct regularization and iterative regularization algorithms. The most common direct regularization algorithms (Tikhonov [21]) and truncated singular value decomposition (TSVD) [22]) are introduced to identify the milling force by measuring acceleration signals. Compared with direct regularization algorithms, iterative regularization algorithms can reduce the amount of calculation and quickly obtain the approximate solution when solving linear discrete ill-posed problems. Especially for solving largescale problems, iterative regularization algorithms will have more advantages. Wang et al. [23] introduced a conjugate gradient least square (CGLS) algorithm to achieve identification of milling force by using filtered acceleration signals. Morigi et al. [24] proved that the LSQR algorithm has better numerical stability than the CGLS algorithm. The iteration stopping criterion of the LSQR algorithm often requires reference value or high time cost, which makes the traditional LSQR algorithm only identify the milling force in theory. Up to now, the indirect identification methods of milling forces have been implemented under offline conditions, and no one has given the process of how to use the proposed method to realize the online identification of milling forces during the milling process.
The purpose of this paper is to realize online identification of dynamic milling force by using acceleration signals and to provide an industry-oriented milling force measurement method. The applied FSC-LSQR algorithm overcomes the excessive number of iterations and inaccurate results caused by the semi-convergence of the LSQR algorithm solution. Based on the FSC-LSQR algorithm, the online identification of milling force during the milling process is realized by using the QBS, which gives the idea of online identification of milling force. Finally, the milling tests validate the effectiveness and generalization of the proposed method. In addition, the CGLS algorithm [23] is also implemented for comparison, which shows that the FSC-LSQR algorithm has higher identification speed and accuracy.

Problem statement of milling force identification
The milling process of the three-flute milling cutter is shown in Fig. 1. Take the tool tip position as the origin, the milling feed direction as the x axis, and the direction perpendicular to the x axis as the y axis to establish a two-degree-of-freedom rectangular coordinate system.
The two-degree-of-freedom equation of motion for the milling process is given as follows: where , and are the 2 × 2 modal mass, damping and stiffness matrices, respectively. (t) = x(t) y(t) T , x(t) and y(t) are tool displacements in the feed (x) and cross-feed (y) directions, and F y (t) are milling forces in the x and y directions.
Since the x(y) direction excitation has much less influence on the y(x) direction response of the spindle system than the same magnitude excitation in the y(x) direction, the cross-talk effects between the x and y directions of the spindle system is not considered, that is, Taking the x direction as an example, the response x(t) under the milling force F x (t) can be calculated by Duhamel integration [25] as follows: where x 0 is initial value, m x is modal mass, x is damping ratio, x is natural frequency, and h xx (t) is transfer function.
According to [23], Eq. (2) can be omitted as follows: The x(t) in Eq. (3) has a general definition: it can be displacement, velocity, or acceleration, and the same as h xx (t) also corresponds to the displacement, velocity, or acceleration transfer functions. The accelerometer has been the most widely used in industrial environments. Therefore, the acceleration signal is used to identify the milling force during the milling process.
In practical applications, Eq. (3) should be discretized and the resulting discrete linear equation as follows: where Δt = t∕n is time interval; n is sampling length.
Equation (4) can be written in a compact form as follows: is the measured acceleration signal; the matrix H xx represents the dynamic characteristics of the system. It can be seen from Eq. (1) that the milling force includes two directions: x direction and y direction. According to Eq. (5), the following equation can be obtained.
where y, H yy , and F y are similar to the derivation of x, H xx , and F x in Eq. (5).
Equation (6) is simplified to the following form: can be seen as the combinatorial vector of different excitation and response points, respectively. From Eq. (7), it can be seen that the acceleration signal can be inverted by the matrix to identify the milling force, and the equation is as follows: From Eq. (8), it seems easy to identify the milling force by using acceleration signals. However, the unknown disturbance is inevitable during the actual milling process, and the obtained matrix is an ill-conditioned matrix that usually satisfies the following three conditions: 1. The singular value gradually decays to zero; 2. The number of conditions is large; 3. Discrete Picard condition is satisfied.
According to the characteristics of the matrix , the identification of milling force is a linear discrete ill-posed problem. We cannot solve Eq. (8) directly. Instead, a feasible method is used to find the approximate solution of as the numerical solution of Eq. (8).
The core problems to be solved are summarized as follows: Problem How to efficiently and accurately solve the approximate solution of the milling force by using the measured acceleration signal, so as to meet the requirements of milling force measurement during the actual milling process.

Online identification method of milling force
During the actual milling process, the accelerometers in the feed and cross-feed directions transmit the acceleration signals to the computer through the driver and DAQ-board. The computer uses the QBS to realize the online identification of the milling force based on the proposed FSC-LSQR algorithm. The specific process is shown in Fig. 2.

Theory of fast stopping criterion for least square QR-factorization algorithm
In order to solve the core problem, the identification method of milling force is required to have a fast solution speed, that is, to efficiently and accurately solve the approximate solution of the milling force in Eq. (7). We refer to the discussion and analysis of a series of related iterative algorithms by Meurant [26] and choose one of the most commonly used and most well-known algorithms, the LSQR algorithm. An LSQR algorithm with a fast stopping criterion that satisfies practical industrial application is given. which can fast and accurate identification of milling force. The LSQR algorithm is described in [27]. When the LSQR algorithm is used to identify the milling force, the initial value of the approximate solution of the milling force is the zero vector 0 = . The approximate solution k for the milling force satisfies the following properties.
It can be seen from Eq. (9) that the LSQR algorithm is essential to minimizing the residuals. At the k-th step, k is searched to minimize the euclidean norm of the residual k = k − . According to the properties of the LSQR algorithm, the residual satisfies the following inequality.
The inequality Eq. (10) is strictly established. According to inequality Eq. (10), it can only guarantee that the residual norm of the obtained approximate solution decreases monotonically as the number of iteration steps k increases. We cannot guarantee that the error norm ‖ ‖ k − ‖ ‖ will also continue to decrease. Normally, the error norm does gradually decrease in the initial few iteration steps, but after it is reduced to a certain level, the value of the error norm begins to increase as the number of iteration steps k further increases. This property exhibited by the iterative solution k is called semi-convergence [28]. The principle of the LSQR algorithm for the approximate solution is analyzed above. The stop iteration criterion of the LSQR algorithm is only to find the minimum residual k , the approximate solution obtained k is inaccurate, and the calculation time will increase as the number of iteration steps k increases. The LSQR algorithm cannot meet the requirements of online identification of milling force in the industry. The objective is to make the LSQR algorithm stop the iteration in time The process of online identification of milling force by using acceleration signals when its iterative solution error norm ‖ ‖ k − ‖ ‖ is as small as possible.
The fast stopping criterion for the LSQR iterative algorithm is considered. The properties of the iterative solutions obtained by two iterative algorithms (i.e., LSQR and Craig) are used to determine the stopping criterion of the LSQR algorithm and the approximate solution F k of the milling force. The specific method is to determine the stopping criterion and the corresponding iterative solution by comparing the residual k of the iterative solution obtained by the LSQR algorithm with the residual ̂k of the iterative solution obtained by the Craig algorithm. The property satisfied by the iterative solution of the Craig algorithm is different from the property of the LSQR algorithm.
where ̂k represents the approximate solution of the milling force obtained by the Craig algorithm.
The purpose of the Craig algorithm is not to minimize residuals, but to minimize errors, which are not easy to describe and measure in practical application. Although Craig algorithm will not simply be used to solve ill-posed problems, the residual ̂k =̂k − obtained by the Craig algorithm is very helpful to better select the stopping criterion of the LSQR algorithm.
An industry-oriented FSC-LSQR algorithm is given to realize the online identification of milling force by using acceleration signals. The FSC-LSQR algorithm is based on the Golub-Kahan bidiagonalization, and can run LSQR and the Craig algorithms at the same time. As the iteration proceeds, we can get the residual k of LSQR algorithm and the residual ̂k of Craig algorithm, and then use ‖ ‖ k ‖ ‖ and ‖ ‖̂k ‖ ‖ to determine when to stop the iterative of the LSQR algorithm.

Remark 1
The specific value of in the FSC-LSQR algorithm needs to be determined according to different industrial environments. The FSC-LSQR algorithm is based on the LSQR algorithm with the introduction of the Craig algorithm. The fast stopping criterion is given by the residual ratio of the Craig algorithm residual ‖ ‖̂k ‖ ‖ to the LSQR algorithm residual ‖ ‖ k ‖ ‖ , which makes the milling force identification more accurate and with fewer iterations. Both the LSQR algorithm and the Craig algorithm are based on the Golub-Kahan bidiagonalization process [24]. The FSC-LSQR algorithm only needs to do one vector product of matrices and T in each iteration of the algorithm. Therefore, the calculation amount required for each iteration of the FSC-LSQR algorithm is essentially the same as that of the LSQR algorithm.

Online identification strategy of milling force
We use the QBS to demonstrate the whole process of the online identification of milling force during the actual milling process based on the FSC-LSQR algorithm, which is not available in the previous paper. In addition, the minimum requirement for the speed of the identification method is given. Only when the identification speed is higher than this requirement, then can the online identification of milling force be realized.
In order to realize the online identification of the milling force, a storage structure that can output the identified milling force while inputting the acceleration signal is needed. The QBS in the data structure has the characteristics of first in first out (FIFO) [29], and then it is applied to online identification of milling force. Taking the QBS containing twenty cutter tooth passing cycles (CTPC) as an example, the strategy for online identification of milling force during the actual milling process is given as shown in Fig. 3. Where queue head pointer front and queue tail pointer rear are used to output milling force (dequeue) and input acceleration signal (enqueue), respectively. id j and j are the identified milling force and acceleration signal in the j-th CTPC, j = 1, ..., n.
It can be seen from Fig. 3 that the queue with a length of twenty CTPC is given, where the first ten CTPC are the identified milling forces, and the last ten CTPC are the acceleration signals. The milling force is dequeued in real-time according to the sampling frequency of the acceleration signal, and the acceleration signal is enqueued at the same time. When the milling force of the tenth CTPC is dequeued, the acceleration signals of the eleventh to the twentieth CTPC have been identified as milling forces. In this way, the online identification of milling force can be realized. This strategy has a requirement, that is, when the queue length is twenty CTPC, the identification time of ten CTPC must be less than the output time of ten CTPC. In other words, if the proposed method is not fast enough for milling force identification, it cannot meet the requirements of online identification of milling force. In the later tests, it is verified that the proposed FSC-LSQR algorithm can output the milling force at a frequency of up to 10,240 Hz in real-time, which fully satisfies the requirements of industrial applications. Note that no matter what identification method of the milling force is used, it takes time. The online identification strategy of the milling force will still have a time delay, which is unavoidable.

Experimental verification
The impact test and milling test are used to obtain the required frequency response functions for the method and to verify the validity of the proposed method, respectively.
As seen in Fig. 4, the accelerometer is mounted along the x and y directions on a box near the front bearing of the spindle. The dynamometer used in the test is the Kistler 9129A dynamometer, and the workpiece material is 6061 aluminum alloy. PCB-086C01 impact hammer (sensitivity is 2.25 mv/N) is used to generate and measure impact force, MI-7008 data acquisition system (sampling frequency is 10240Hz) is used to record impact force and the response. A series of downmilling tests are carried out in a three-axis CNC machining center (VMC-V5) by using a three-flute milling cutter (high-speed steel material, diameter 10 mm, flute length 40 mm, overhang length 60 mm, helix angle 45 • ). Table 1 lists the milling parameters of the conducted milling tests. The parameters in Table 1 have been set to be the values that can guarantee stable milling, and the change in spindle speed and milling width strengthens the persuasiveness of the results. Figure 5 shows the mounting points for accelerometers. To ensure the same direction of the Kistler dynamometer and the accelerometer, two screw holes (x and y directions) are machined in the spindle for mounting the accelerometer.
The accelerometer is mounted in a position that minimizes the transmission path of milling forces, minimizes energy loss, and has higher sensitivity.

Frequency response function measurement
The impact tests are carried out to obtain the frequency response function between the tool tip and the accelerometer. The frequency response functions (FRF-xx, FRF-yy) are measured five times by impact tests along the x and y directions, and the averaged frequency response functions are calculated for milling force identification. This method of averaging multiple impact tests can reduce the influence of random errors and make the obtained frequency response functions more accurate. The frequency response functions obtained through the impact tests are shown in Fig. 6.

Dynamic compensation of Kistler dynamometer
The output value of the Kistler dynamometer is compensated, and then the compensated value is used as the reference standard for identifying milling force. The reason for the distortion of the dynamometer measurement results is explained. The direct measurement system for milling forces (Fig. 7) consists of the workpiece, the bolts, and the dynamometer. The workpiece and the dynamometer are connected by a total of ten bolts in the front and rear rows. The dynamometer is also fixed on the workbench by bolts. The milling forces generated during the milling of the workpiece are measured using a Kistler dynamometer.
The input and output of the milling force measurement system are mi and mo , respectively. The relationship between mi and mo is given as follows: where f ( ) denotes the frequency response function between measured force and milling force, mo ( ) denotes the frequency domain representation of the measured force mo , and mi ( ) denotes the frequency domain representation of the milling force mi .
The frequency response function between the measured force and the milling force obtained by the impact test is shown in Fig. 8. Figure 8a, b correspond to the frequency response functions in the feed (x) and cross-feed (y) directions, respectively. In ideal situations, the measurement  result of the dynamometer is not affected by the dynamic characteristics of the milling force measurement system, that is, the amplitude of the frequency response function is 1 in the entire frequency bandwidth. As shown in Fig. 8, when the FRF−xx and FRF−yy are in the range of 0-500 Hz and 0-550 Hz, respectively, the amplitude of the frequency response function is close to 1; when the frequency is higher, the amplitude of the frequency response function is much greater than 1. Distortion of milling force measurement is not a unique phenomenon of the dynamometer. In fact, the influence of the dynamic characteristics of the milling force measurement system is widespread. Therefore, it is necessary to compensate for the measured milling forces.
Because the focus of this paper is not the dynamic compensation of the Kistler dynamometer, the improved inverse filter proposed by Wan et al. [30] is used to compensate for the measured milling force. The compensated milling force is regarded as the actual milling force , which serves as a reference standard for verifying the accuracy of the identified milling force.

Comparison of milling force identification results
The frequency spectrum components of milling force are mainly distributed on the tooth passing frequency f c and its multipliers. Where f c = NΩ∕60 , Ω is the spindle speed, and N is the number of milling cutter teeth. Inspired by [23], Butterworth band-pass filtering is chosen to remove low-and high-frequency noise for more accurate identification results. The frequency response functions in the x and y directions are obtained by impact tests. The acceleration signal measured in the milling test is passed through Butterworth band-pass filtering as input. The measured milling force after inverse filter compensation is used as the reference standard to verify the effectiveness of the proposed method. A fixed number of iterations is given to verify the effectiveness of the FSC-LSQR algorithm. Figure 9 shows the error norm of the milling force identified by the LSQR algorithm and the measured milling force after compensation under 100 iterations. Figure 10 shows the residual norm of the FSC-LSQR algorithm under 100 iterations.
As can be seen in Figs. 9 and 10, if only the LSQR algorithm is used, the error norm first decreases rapidly and then slowly increases as the number of iterations increases, which verifies that the error norm ( ‖F x k − F x ‖ and ‖F y k − F y ‖ ) of the approximate solution mentioned above does not necessarily satisfy monotonically decreasing. The residual norm ( ‖r x k ‖ = ‖H xx F x k − x‖ and ‖r y k ‖ = ‖H yy F y k − y‖ ) satisfies the inequality Eq. (10) that ‖r x k ‖ ≥ ‖r x k+1 ‖ and ‖r y k ‖ ≥ ‖r y k+1 ‖ . During the milling process, it is impossible to equip each machine tool with the reference value of the milling force provided by the Kistler dynamometer to obtain the optimal number of iterations. Morigi et al. [24] seemed to obtain a satisfactory approximate solution and its stopping criterion by finding the local minimum points of ‖ k+1 − k ‖ . However, it is pointed out in the [31] that the effect of this algorithm in numerical experiments is not very well.
As mentioned in Sect. 3.1, both the LSQR algorithm and the Craig algorithm are based on the Golub-Kahan bidiagonalization process, so the solution speed of the FSC-LSQR algorithm is almost the same as that of the LSQR algorithm in each iteration. In general, the ‖r x k ‖ and ‖r y k ‖ obtained by the Craig algorithm will be greater than the ‖r x k ‖ and ‖r y k ‖ obtained by the LSQR algorithm. Hanke [28] has given some numerical verifications on this property. We determine the value of parameter in the FSC-LSQR algorithm through milling tests.

Remark 2
The advantages of the FSC-LSQR algorithm in theory and practical applications are summarized as follows: (1) In theory, the FSC-LSQR algorithm gives the iterative stopping criterion of the optimal solution, so that the milling force can be identified quickly and accurately. (2) In  practical applications, the FSC-LSQR algorithm does not require the use of dynamometer results as the reference to give the iterative stopping criterion, and the optimal solution of the FSC-LSQR algorithm can be obtained by using the accelerometer.
The results obtained by using the FSC-LSQR and the CGLS algorithms to identify the milling force are shown in Figs. 11,12,13,14,15,and 16. From Figs. 11,12,13,14,15,and 16, it is obvious that the milling force identified by the FSC-LSQR algorithm is in good agreement with the measured milling force after compensation. There is a large error between the milling force identified by the CGLS algorithm and the measured milling force after compensation. From the magnified view, there is still a certain difference between the milling force identified by the FSC-LSQR algorithm and the measured milling force after compensation, which is mainly affected by the complicated path of the acceleration signal from the tool tip to the spindle box and the nonlinearity of the system.

Remark 3
Wang et al. [23] have demonstrated that the CGLS algorithm identifies milling forces with higher accuracy and speed than the Tikhonov regularization algorithm [32], so the proposed method is only compared with the CGLS algorithm. Figure 17 shows the absolute errors between the measured milling force after compensation and the identified milling force in the x and y directions, which illustrates that the milling force identified by the FSC-LSQR algorithm has less absolute error than the milling force identified by the CGLS algorithm.
Quantitative analysis of milling force identification based on CGLS and FSC-LSQR algorithms is carried out for milling condition I. The most common evaluation indexes (peak-to-peak value and root mean square value) are used to give the evaluation standard for the accuracy of milling force identification. where N the number of times the CTPC, P cp i and P id i are the peak-to-peak value of the measured milling force after compensation and the identified milling force for the i-th CTPC, and R cp i and R id i are the root mean square value of the measured milling force after compensation and the identified milling force for the i-th CTPC.
The CGLS and FSC-LSQR algorithms are run in MAT-LAB R2018a on the same computer to identify the milling force. CPU is AMD Ryzen 7 3700X 8-Core Processor 3.59GHz, RAM is 16G, and the operating system is Windows 10. Note that the quantitative analysis of ten CTPC is very meaningful, the online identification strategy of milling force given by the QBS in Sect. 3.2 is to take the acceleration signal of ten CTPC as an example. In addition, only the milling force of a complete CTPC can obtain additional information, such as tool breakage, work path deviation, etc. The results of the quantitative analysis are shown in Fig. 18.
It can be seen from Fig. 18 that the FSC-LSQR algorithm is superior to the CGLS algorithm in terms of both the peak-to-peak evaluation index ( V pp ) and the root mean square evaluation index ( V rms ). The most important thing is that the running time of the FSC-LSQR algorithm is much lower than that of the CGLS algorithm, which provides the possibility of online identification of milling force.

Discussions
According to the online identification strategy of milling force given in Sect. 3.2, if the running time of the algorithm to identify milling force for ten CTPC is less than the time of milling force dequeue, the online identification of milling force can be realized. The following is an example of milling condition I. According to the sampling frequency of 10,240 Hz and spindle speed of 4000 r/min, the output time of 512 points collected in ten CTPC is 0.05 s. Fig. 18 shows the running time of the milling force identification for ten CTPC using FSC-LSQR and CGLS algorithms. Among them, the time of the milling force identification using the FSC-LSQR algorithm is less than 0.05 s. In other words, the proposed FSC-LSQR algorithm can realize the online identification of milling force with the sampling frequency of 10,240 Hz. In addition, it is found that the use of acceleration signals for milling force identification gives an accurate AC component, which is the most important component in milling force identification. The DC component is compensated by the conclusions obtained from the existing milling force model.