Identiﬁcation of Robot Joint Torsional Stiffness Based on the Amplitude of the Frequency Response of Asynchronous Data

: The difﬁculty of adding external excitation and the asynchronous data collection from the industrial robot operation limited the online parameter identiﬁcation of industrial robots. In this regard, this study proposes an identiﬁcation method that only uses the amplitude of the frequency response function (FRF) of the system to identify robot joint torsional stiffness and dynamic parameters. The error criterion function shows that this method is feasible and comparable to applying the complete frequency response for identiﬁcation. The Levenberg–Marquardt (L-M) algorithm is used to ﬁnd the global optimal value of the error criterion function. In addition, an operational excitation method is proposed to excite the system. The speed proﬁle is set as a triangle wave to excite the system using rectangular wave electromagnetic torques. The simulation results show that using the amplitude of the FRF to identify parameters applies to asynchronous data. The experiments on a single-degree-of-freedom articulated arm test bench show that the motion excitation method is effective, and both stiffness and inertia are identiﬁable.


Introduction
Servo systems have become wildly used with the development of industry and technology, particularly in high-end medical equipment, new energy, robots, CNC machine tools, and other fields [1,2]. Due to the fatigue of the mechanical part of the actual system, the stiffness of the shaft system declines, and the performance of the transmission mechanism reduces as the running time increases. Thus, the introduction of greater motion error in the system is inevitable [3,4]. Therefore, changes in the shaft system stiffness can reflect the operating status of the system and are essential indicators of the performance evaluation and prediction of the transmission system [5][6][7]. Identifying the stiffness of the shaft system during operation is crucial to predict an early failure of an industrial robot [8].
Parameter identification is typically used to estimate the parameters of a mathematical model based on the input and output information of a particular system [9]. The parameter identification in the servo drive system includes two aspects: internal parameter identification and external parameter identification [10,11]. For the dual inertia system with nonlinear characteristics, Beineke [12] proposed two identification methods for external mechanical parameters: standard minimum root mean square and four-step instrumental variable method. They extracted the speed response of the motor at different stages for the dual inertia system with nonlinear characteristics and identified the parameters corresponding to each stage. The identification results were compared with the traditional neural network, self-organizing mapping, and classic frequency domain analysis methods. Villwock and Pacas [13,14] proposed a non-parametric identification method to obtain the resonant frequency of the dual inertia elastic system directly. In the identification process, In this study, only one axis motion of the industrial robot is considered. The ty industrial robot joint is composed of a permanent magnet synchronous motor (PMS speed reducer, and an arm, as shown in Figure 1. In the figure, e T is electromag torque, m  is motor speed, l T is the load torque, and l  is the rotating speed o articulated arm. The joint structure of the robot is complex. Moreover, identifying each structura rameter is challenging. However, a robot joint can be equivalent to an identifiable inertia model. Equations (1) and (2) can be used to calculate the equivalent stiffness equivalent load inertia of the joint.   The joint structure of the robot is complex. Moreover, identifying each structural parameter is challenging. However, a robot joint can be equivalent to an identifiable dual inertia model. Equations (1) and (2) can be used to calculate the equivalent stiffness and equivalent load inertia of the joint.
where K 1 , K 2 , · · ·, K n are the stiffness of each part, or a series of parts, of the shaft. K is the equivalent joint stiffness. J L and J l are the actual load inertia and equivalent load inertia, respectively. i is the reduction ratio. The double inertia model of the system can be obtained, as shown in Figure 2.
In this study, only one axis motion of the industrial robot is considered. The typ industrial robot joint is composed of a permanent magnet synchronous motor (PMSM speed reducer, and an arm, as shown in Figure 1. In the figure, e T is electromagn torque, m  is motor speed, l T is the load torque, and l  is the rotating speed of articulated arm. The joint structure of the robot is complex. Moreover, identifying each structural rameter is challenging. However, a robot joint can be equivalent to an identifiable d inertia model. Equations (1) and (2) can be used to calculate the equivalent stiffness equivalent load inertia of the joint.    The mechanical dynamics of the system can be written as follows: where J m is the moment of inertia of the servo motor. b m and b l are the damping of the motor and the swing, respectively. θ m and θ l are the position angle of the motor shaft and articulated arm, respectively. T s is the shaft torque, C is the damping of the shaft, F is the gravity of the swing arm, and L is the distance from the shaft to the line of action of the load.
According to Equation (3), the transfer function from the electromagnetic torque to motor speed can be obtained as follows: Equation (4) shows that the motor speed feedback introduces conjugate zeros and conjugate poles to the closed-loop system, which will lead to mechanical resonance. Through the conjugate zeros and the conjugate pole presented in Equation (4), the anti-resonance frequency (ARF) ω ARF and the natural torsional frequency (NTF) ω NTF can be derived. If b m = 0, b l = 0, and C = 0 are assumed, the ω ARF and ω NTF can be expressed as follows: Since the value of inertia changes by a very small amount in the long-time operation, the stiffness determines the change of the resonance frequencies of the system.

The Excitation Signal of the Swing Movement
The direct closed-loop identification setup shown in Figure 3 is considered, where C(s) is the linear transfer function of the controller, P(s) is the transfer function of the drive system, u(t) is the system input, and y(t) is the system output. ω r is the reference rotating speed. The novelty of the identification setup is that it does not use external excitation signal because it uses operational excitation.  The identification accuracy depends heavily on the excitation signal, which should cover the frequency components contained in the system. The most used excitation signals are swept-frequency signals and PRBS [20,21]. The excitation signal is superimposed on the torque reference obtained from the speed controller as the input signal to excite the system. However, if the frequency band of the torque itself is wide enough, and the signal strength in the frequency band is large enough, the system can also be excited. Therefore, The identification accuracy depends heavily on the excitation signal, which should cover the frequency components contained in the system. The most used excitation signals are swept-frequency signals and PRBS [20,21]. The excitation signal is superimposed on the torque reference obtained from the speed controller as the input signal to excite the system. However, if the frequency band of the torque itself is wide enough, and the signal strength in the frequency band is large enough, the system can also be excited. Therefore, we considered the electromagnetic torque generated by the reciprocating motion of industrial robots as system's excitation.
The commonly used speed profile of reciprocating motion of robot joints is the trapezoidal curve. The trapezoidal curve divides the entire movement process into three stages: constant acceleration, constant speed, and constant deceleration [22]. When the motion period is short and the maximum speed is not reached, the trapezoidal curve becomes a triangular wave. For this case, the electromagnetic torque corresponds to a rectangular wave. The frequency spectrum of a rectangular wave is wideband, which can achieve a similar effect to the external excitation.

Calculation of the Electromagnetic Torque
Torque sensors cannot be implemented to measure torque in industrial robot joints. Thus, the electromagnetic torque is typically calculated from the line current. In this regard, it is necessary to perform a coordinate transformation in the alternating current (AC) motor physical model to obtain the relationship between the electromagnetic torque and the line current [23].
Coordinate transformation follows the principle of constant power before and after transformation, mainly including the following two forms: (1) Clarke transformation can be written as follows: (2) Park transformation can be written as follows: where i A , i B , i C are three-phase alternating currents. i α and i β are two-phase alternating current in the two-phase stationary coordinate system αβ. i d and i q are two direct currents in the two-phase rotating coordinate system dq. C 3s/2s is the Clarke transformation matrix, while C 2s/2r is the Park transformation matrix. θ represents the electrical angle. After the coordinate transformation, the decoupling of torque components is performed, and the PMSM mathematical model in the dq coordinate system is obtained. The torque equation can be written as follows: If i d is controlled to be zero, PMSM is decoupled, and the torque of PMSM can be written as follows [24,25]: where ϕ d and ϕ q are the d and q axis components of the stator flux linkage, p is the number of pole pairs, and k t is the torque constant.
For i d = 0, i q is equal to the peak value of the line current [26]. As the load torque varies with the position angle, the peak of the line current varies accordingly. As the line current is an analytical signal, its peak curve can be detected as the envelope using Hilbert transform (HT). The modulus of the analyzed signal is the instantaneous amplitude or envelope of the signal [27]. Assuming that the measured line current is i a (t), and its HT is denoted as h[i a (t)], then i q (t) and T e (t) can be obtained as follows:

Error Criterion Function of the Identification Method of the Amplitude of the FRF
In the operation of the industrial robot, it is difficult for conventional signal acquisition cards to meet the sampling rate requirements when analog sampling is used due to the high accuracy of the encoder. A better way to obtain the speed is to read the encoder data from the driver. This type of signal is digital, which is different from the signal type of current. Thus, the input and the output data are asynchronous, causing deviations in the phasefrequency characteristics of the system, which makes it impossible to identify the system accurately by frequency response. However, the amplitude-frequency characteristics are not affected; therefore, this study proposes a method for parameter identification using the amplitude of the FRF. The following analyzes the identification performance of this method.
In direct closed-loop identification setups, the transfer function in Equation (4) can also be written as follows: where The corresponding FRF is When using the FRF for identification, the error between the frequency response estimate ∧ H(ω) and the actual frequency response H(ω) is The generalized error criterion function can be written as follows: The function can be written further as Equation (18)  minimizes the error criterion function, the partial derivative of the error criterion function with respect to b 1 is zero.
When only the amplitude of the FRF is used for identification, the amplitude of the FRF estimate ∧ H(ω) and the actual amplitude of the FRF |H(ω)| should be equal in the absence of error.
Squaring both sides simultaneously The root of the error function is removed by considering the error. The generalized error criterion function can be written as follows: The partial derivative of the error criterion function E a with respect to b 1 is Comparing Equation (26) with Equation (20), the results obtained by the two equations differ. However, by transforming Equation (20), the following equation can be obtained: The difference between Equations (26) and (27) is According to Equation (14), when the damping C, b 1 , and b 2 are zero, L is zero. The damping in the actual system is difficult to identify accurately through parameter identification. Typically, the damping parameters are ignored under various complex working conditions [28]. Therefore, when the damping can be neglected, parameter identification using the amplitude of the FRF is equivalent to that using the frequency response.

Calculation of the Frequency Response Function
Generally, the FRF of the system input signal x(t) and output signal y(t) is given by the ratio between the spectral density function G xy (ω) and G xx (ω) as follows: where H(ω) is the FRF of the system, G xy (ω) is the cross-power spectral density of the input x(t) and output y(t), and G xx (ω) is the self-power spectral density of the input signal x(t).

Levenberg-Marquardt Algorithm
Section 3.3 introduced the error criterion function, and this section uses the L-M algorithm to iterate on it to find the optimal value.
Error function: For Newton's law, the iterative formula is: where e k (x) is the error, n is the total number of points in the identification frequency domain window, and x are the parameters to be identified. x k is the input vector at the k-th iteration, and x k+1 is the input vector at the (k+1)-th iteration. ∇E(x) and ∇ 2 E(x) are the gradient and the Hessian matrix of the error criterion function E(x) respectively.
In the above two formulas, S(x) = n ∑ k=1 e k (x)∇ 2 e k (x) , J is the Jacobian matrix, namely: The calculation method of the Gauss-Newton method is: The L-M algorithm is an improved form, namely: where I is the identity matrix. µ is the damping coefficient, and µ needs to be adjusted during each iteration: When L decreases quickly, a smaller µ value is used. At this time, the method is similar to the Gauss-Newton method. When L decreases slowly, the µ value is increased, and the method is similar to the gradient method [29]. Figure 4 shows the flowchart of the parameter identification process, which can be summarized as follows: (1) The dynamic model of a single robot joint is established.
(2) The speed command signal is set as a triangle wave so that the electromagnetic torque follows a rectangular window to excite the system. (3) The current signal and the motor encoder signal are collected to calculate the electromagnetic torque and the motor speed. (4) The ratio between the cross-power spectrum of input and output signals and the self-power spectrum of input signals is calculated to obtain the amplitude of the FRF. (2) The speed command signal is set as a triangle wave so that the electromagn torque follows a rectangular window to excite the system. (3) The current signal and the motor encoder signal are collected to calculate the elec magnetic torque and the motor speed. (4) The ratio between the cross-power spectrum of input and output signals and the s power spectrum of input signals is calculated to obtain the amplitude of the FRF (5) The L-M optimization algorithm is used to fit the experimental and theoretical a plitude of the FRF to minimize the error.

Simulations
According to the dynamic equations of the single joint degree of freedom system troduced in Section 1, a simulation model of the articulated arm motion was built. parameters used in the simulation correspond to the actual parameters of the test ben the structure of the test bench and its parameters are described in detail in Section 5 The joint stiffness is the equivalent stiffness of the coupling stiffness and the reducer s ness in series, and the value is 891N m / rad .The inertia of the motor side is the sum of inertia of the motor and the inertia of the reducer, and the value is 2 0.003027kg m . equivalent inertia on the load side is 2 0.00748kg m . The speed command was set as a triangular wave to excite the system using rect gular wave electromagnetic torques. The top speed was 1000 r/min, and the movem period was 0.4 s, and the total simulation time was set to 2 s.
The amplitude of the FRF and the phase of the FRF of the system are performed us synchronous data and asynchronous data, respectively, as shown in the Figure 5. Acco ing to Equations (5) and (6)

 =
Hz. The two frequencies are c sistent with the resonance frequency obtained by the simulation model in Figure 5, wh shows the accuracy of the simulation model. The amplitude frequency characteristics tained by asynchronous data and synchronous data are the same. However, the asynch nous data has a phase shift.

Simulations
According to the dynamic equations of the single joint degree of freedom system introduced in Section 1, a simulation model of the articulated arm motion was built. The parameters used in the simulation correspond to the actual parameters of the test bench, the structure of the test bench and its parameters are described in detail in Section 5.1. The joint stiffness is the equivalent stiffness of the coupling stiffness and the reducer stiffness in series, and the value is 891 N·m/rad.The inertia of the motor side is the sum of the inertia of the motor and the inertia of the reducer, and the value is 0.003027 kg·m 2 . The equivalent inertia on the load side is 0.00748 kg·m 2 .
The speed command was set as a triangular wave to excite the system using rectangular wave electromagnetic torques. The top speed was 1000 r/min, and the movement period was 0.4 s, and the total simulation time was set to 2 s.
The amplitude of the FRF and the phase of the FRF of the system are performed using synchronous data and asynchronous data, respectively, as shown in the Figure 5. According to Equations (5) and (6), the theoretical resonance frequencies ARF and ATF of the system are calculated as ω ARF = 54.96 Hz, ω NTF = 102.39 Hz. The two frequencies are consistent with the resonance frequency obtained by the simulation model in Figure 5, which shows the accuracy of the simulation model. The amplitude frequency characteristics obtained by asynchronous data and synchronous data are the same. However, the asynchronous data has a phase shift. Data in the range of 0-500 Hz in the FRF and the amplitude of the FRF, respectively, were used to identify the stiffness of the shaft system using the L-M optimization algorithm. Firstly, the conventional identification method based on the FRF was used to identify the synchronous data. In this case, the identification result is 866.13 N m / rad and the error is 2.79%, the identification accuracy is satisfactory. However, for asynchronous data, identification result of the identification method based on the FRF is 654.67 N m / rad and the error reaches 26.52%. The result is unacceptable, it shows that the FRF is not applicable for asynchronous data.
However, the identification method of the amplitude of the FRF obtains accurate results for asynchronous data. The stiffness identification result is 859. 13 N m / rad and the error is 3.49%. Table 1 shows the stiffness identification results of the three cases. It proves that the proposed method can avoid phase error, and it is feasible and comparable to applying the complete frequency response for identification.

Experimental System
A single-degree-of-freedom articulated arm test bench was built to simulate the motion of an industrial robot single joint. The entire setup is shown in Figure 6. The transmission system comprises a speed reducer and a shaft. Moreover, the elastic coupling is connected in series with the shaft and the speed reducer. The PLC model is a SIEMENS 1212c connected to the driver to control the motor motion. An oscilloscope is used to obtain the encoder signal from the driver. The current sensor is placed on one of the power lines, and the NI9234 acquisition card is used to collect the line current data. The input and output data collected are asynchronous. The specific parameters of the test bench are listed in Table 2. Data in the range of 0-500 Hz in the FRF and the amplitude of the FRF, respectively, were used to identify the stiffness of the shaft system using the L-M optimization algorithm. Firstly, the conventional identification method based on the FRF was used to identify the synchronous data. In this case, the identification result is 866.13 N·m/rad and the error is 2.79%, the identification accuracy is satisfactory. However, for asynchronous data, identification result of the identification method based on the FRF is 654.67 N·m/rad and the error reaches 26.52%. The result is unacceptable, it shows that the FRF is not applicable for asynchronous data.
However, the identification method of the amplitude of the FRF obtains accurate results for asynchronous data. The stiffness identification result is 859.13 N·m/rad and the error is 3.49%. Table 1 shows the stiffness identification results of the three cases. It proves that the proposed method can avoid phase error, and it is feasible and comparable to applying the complete frequency response for identification.

Experimental System
A single-degree-of-freedom articulated arm test bench was built to simulate the motion of an industrial robot single joint. The entire setup is shown in Figure 6. The transmission system comprises a speed reducer and a shaft. Moreover, the elastic coupling is connected in series with the shaft and the speed reducer. The PLC model is a SIEMENS 1212c connected to the driver to control the motor motion. An oscilloscope is used to obtain the encoder signal from the driver. The current sensor is placed on one of the power lines, and the NI9234 acquisition card is used to collect the line current data. The input and output data collected are asynchronous. The specific parameters of the test bench are listed in Table 2.

900N m / rad
During the experiment, the triangular wave was set as the command speed to produce the electromagnetic torque as a rectangular wave to excite the system. Two motion conditions were set, and the acceleration time and deceleration time of the first group were both set to 0.5 s. Furthermore, the acceleration time and deceleration time of the second group were both set to 0.2 s. In these two conditions, the joint arm swings back and forth 120 degrees in each motion cycle.

The Impact of the Acceleration Time on Excitation
The electromagnetic torque is obtained by multiplying the envelope of the current and the motor torque constant. The motor encoder signal is processed to obtain the motor speed. Because the current signal and the motor encoder signal are collected separately and have different sampling rates, it is necessary to resample the sampling rate of the two signals through interpolation so that they are consistent. Figure 7 shows the input signals and the output signals of the test bench for an acceleration time of 0.5 s. Figure 8 shows the frequency response of the system when the acceleration time is 0.5 s. No evident resonance peaks are found on the amplitude and phase of the FRF, indicating that the system is not sufficiently excited, and it is impossible to identify the parameters.  During the experiment, the triangular wave was set as the command speed to produce the electromagnetic torque as a rectangular wave to excite the system. Two motion conditions were set, and the acceleration time and deceleration time of the first group were both set to 0.5 s. Furthermore, the acceleration time and deceleration time of the second group were both set to 0.2 s. In these two conditions, the joint arm swings back and forth 120 degrees in each motion cycle.

The Impact of the Acceleration Time on Excitation
The electromagnetic torque is obtained by multiplying the envelope of the current and the motor torque constant. The motor encoder signal is processed to obtain the motor speed. Because the current signal and the motor encoder signal are collected separately and have different sampling rates, it is necessary to resample the sampling rate of the two signals through interpolation so that they are consistent. Figure 7 shows the input signals and the output signals of the test bench for an acceleration time of 0.5 s. Figure 8 shows the frequency response of the system when the acceleration time is 0.5 s. No evident resonance peaks are found on the amplitude and phase of the FRF, indicating that the system is not sufficiently excited, and it is impossible to identify the parameters. Machines 2021, 9, x FOR PEER REVIEW 13 of 17  Then, the acceleration time is set to 0.2 s, and the torque and speed obtained are shown in Figure 9. The amplitude and the phase of the FRF of the system are obtained, as shown in Figures 10 and 11. Evident resonance can be observed in the two Figures.  Then, the acceleration time is set to 0.2 s, and the torque and speed obtained are shown in Figure 9. The amplitude and the phase of the FRF of the system are obtained, as shown in Figures 10 and 11. Evident resonance can be observed in the two Figures. Then, the acceleration time is set to 0.2 s, and the torque and speed obtained are shown in Figure 9. The amplitude and the phase of the FRF of the system are obtained, as shown in Figures 10 and 11. Evident resonance can be observed in the two Figures.  Machines 2021, 9, x FOR PEER REVIEW 14 of 17    Figure 12 shows the power spectrum of the electromagnetic torque at 0.2 s and 0.5 s acceleration time. The power spectrum of the two types of electromagnetic torque has a wide enough frequency band. However, the intensity decreases with increasing frequency. At the position of resonance frequency, the intensity of the two types of excitation signals is not significant. The amplitude is approximately −170 dB when the acceleration time is 0.5 s, and the amplitude is approximately −100dB when the acceleration time is 0.2 s. The success of the latter excites the system, indicating that the shorter the reciprocating time, the stronger the excitation signal, and the easier it is to excite the system.    Figure 12 shows the power spectrum of the electromagnetic torque at 0.2 s and 0.5 s acceleration time. The power spectrum of the two types of electromagnetic torque has a wide enough frequency band. However, the intensity decreases with increasing frequency. At the position of resonance frequency, the intensity of the two types of excitation signals is not significant. The amplitude is approximately −170 dB when the acceleration time is 0.5 s, and the amplitude is approximately −100dB when the acceleration time is 0.2 s. The success of the latter excites the system, indicating that the shorter the reciprocating time, the stronger the excitation signal, and the easier it is to excite the system.    Figure 12 shows the power spectrum of the electromagnetic torque at 0.2 s and 0.5 s acceleration time. The power spectrum of the two types of electromagnetic torque has a wide enough frequency band. However, the intensity decreases with increasing frequency. At the position of resonance frequency, the intensity of the two types of excitation signals is not significant. The amplitude is approximately −170 dB when the acceleration time is 0.5 s, and the amplitude is approximately −100dB when the acceleration time is 0.2 s. The success of the latter excites the system, indicating that the shorter the reciprocating time, the stronger the excitation signal, and the easier it is to excite the system.  Figure 12 shows the power spectrum of the electromagnetic torque at 0.2 s and 0.5 s acceleration time. The power spectrum of the two types of electromagnetic torque has a wide enough frequency band. However, the intensity decreases with increasing frequency. At the position of resonance frequency, the intensity of the two types of excitation signals is not significant. The amplitude is approximately −170 dB when the acceleration time is 0.5 s, and the amplitude is approximately −100 dB when the acceleration time is 0.2 s. The success of the latter excites the system, indicating that the shorter the reciprocating time, the stronger the excitation signal, and the easier it is to excite the system.

Parameter Identification
The FRF and the amplitude of the FRF of the system with an acceleration time of 0.2 s were used for identification respectively, including the frequency range of 0-500 Hz. Because the torsional stiffness was our focus, we only identified stiffness first and considered that inertias were known. The identification results of the two methods are shown in Table 3. Since the data collected in the experiment are asynchronous, there is phase deviation in the FRF. Therefore, when the traditional identification method based on FRF is used, the identification error reaches 37.93%. The proposed identification method is based on the amplitude of the FRF, which avoids phase error and obtains accurate identification results. The identified value is 851.9N m / rad , and the error is 4.49%. The effectiveness of the method was tested for multiple unknowns. In this regard, the two inertia parameters were regarded as unknowns and identified along with torsional stiffness. The results are listed in Table 4. As seen in the table, the identification accuracy of torsional stiffness declines when the unknown parameters increase compared with that of the simulation; the error of the torsional stiffness reaches 9.39%. Thus, simultaneously identifying multiple parameters increases the calculation error. In contrast, the identification errors of the three parameters are all below 10%, indicating that not only the stiffness can be identified. However, motor inertia and load inertia can also be identified.  Figure 10 shows the spectrum of the amplitude of the FRF generated by the experimental and identified data. As seen in the Figure, the identified resonance frequencies agree well with the experimental resonance frequencies when only the locations of the anti-resonance and resonance frequencies are considered. However, note that a certain difference exists between the identified amplitude and the experimental amplitude. The

Parameter Identification
The FRF and the amplitude of the FRF of the system with an acceleration time of 0.2 s were used for identification respectively, including the frequency range of 0-500 Hz. Because the torsional stiffness was our focus, we only identified stiffness first and considered that inertias were known. The identification results of the two methods are shown in Table 3. Since the data collected in the experiment are asynchronous, there is phase deviation in the FRF. Therefore, when the traditional identification method based on FRF is used, the identification error reaches 37.93%. The proposed identification method is based on the amplitude of the FRF, which avoids phase error and obtains accurate identification results. The identified value is 851.9 N·m/rad, and the error is 4.49%. The effectiveness of the method was tested for multiple unknowns. In this regard, the two inertia parameters were regarded as unknowns and identified along with torsional stiffness. The results are listed in Table 4. As seen in the table, the identification accuracy of torsional stiffness declines when the unknown parameters increase compared with that of the simulation; the error of the torsional stiffness reaches 9.39%. Thus, simultaneously identifying multiple parameters increases the calculation error. In contrast, the identification errors of the three parameters are all below 10%, indicating that not only the stiffness can be identified. However, motor inertia and load inertia can also be identified.  Figure 10 shows the spectrum of the amplitude of the FRF generated by the experimental and identified data. As seen in the Figure, the identified resonance frequencies agree well with the experimental resonance frequencies when only the locations of the anti-resonance and resonance frequencies are considered. However, note that a certain difference exists between the identified amplitude and the experimental amplitude. The reason is the damping effect that is neglected in the identification, and there is noise interference in the experiment. Figure 11 shows the spectrum of the phase of the FRF generated by the experimental and identified data. When the input and output data of the experiment are asynchronous, although the resonance frequency can be identified, the phase of the FRF of the experimental data and the identification data have a substantial deviation. Therefore, the phase of the FRF leads to identification errors when the input and output data are asynchronous. Nevertheless, the identification method of the amplitude of the FRF avoids the phase error, which can replace the identification method based on the FRF.

Conclusions
In this study, aiming at identifying the torsional stiffness of industrial robot joints, we proposed an operational identification method. The novelty of this method relies on using motion excitation to excite the system and the amplitude of the FRF for parameter identification. The proposed method can function without affecting the normal working state of the robot. The experiments and simulations showed the effectiveness of the proposed method. More importantly, the stiffness, motor inertia, and load inertia could be identified. The future work will focus on further enriching the identifiable parameters, including the identification of electrical parameters.