Industrial mechanical arm joint rigidity online identification method based on optimal excitation trajectory

Abstract

The invention discloses an industrial mechanical arm joint rigidity online identification method based on an optimal excitation track, and belongs to the field of mechanical arm parameter identification. The method comprises the following steps: firstly, establishing a mechanical arm rigid-flexible coupling dynamic model, and constructing a linear regression equation taking joint rigidity and viscous damping as parameters; a Fisher information matrix is derived based on the regression model to quantify the amount of information of the excitation trajectory. And further, parameterizing the excitation trajectory by using finite term Fourier series, and solving the optimal excitation trajectory by using a scalar index of a Fisher information matrix and a sustainable excitation condition as optimization targets under the condition of meeting joint motion constraints. And finally, the mechanical arm is driven to execute the track, motion data are collected and input into parameter identification algorithms such as unscented Kalman filtering and recursive least square, and online and high-precision estimation of the time-varying joint stiffness is achieved. According to the method, the excitation track is designed directly from the information theory, and the precision, robustness and real-time performance of joint rigidity identification are effectively improved.

Description

Technical Field

[0001] This invention belongs to the field of robotic arm parameter identification technology, specifically a method for online identification of joint stiffness of industrial robotic arms based on optimal excitation trajectory. It involves excitation trajectory design and online identification of joint stiffness in parameter identification, aiming to improve the accuracy and robustness of online identification of joint stiffness of robotic arms by designing the optimal excitation trajectory. Background Technology

[0002] In robotic systems, harmonic reducers and RV reducers, with their advantages of high transmission ratio and large torque, have become core components of joint transmission, and their performance directly determines the stability and accuracy of robot motion. Ideally, the reducer possesses constant torsional stiffness, but under the influence of dynamic loads, machining errors, assembly errors, and periodic changes in gear meshing rates, it exhibits significant stiffness fluctuations. This time-varying stiffness, as the main excitation source within traditional systems, is easily coupled with the system's natural frequency, inducing complex torsional vibrations and harmonic resonance risks, directly leading to problems such as chattering and decreased positioning accuracy at the robot's end effector. Therefore, identifying and studying joint stiffness as a time-varying parameter, rather than simplifying it to a constant value, is a necessary prerequisite for improving robot accuracy and stability. Accurately identifying joint dynamic stiffness is the foundation for designing control strategies and is of great significance for applications such as precision assembly and high-end manufacturing.

[0003] Existing technologies have significant limitations in the design of excitation trajectories for parameter identification. Most current methods focus on minimizing the condition number of the observation matrix, which, while suppressing noise to some extent, fails to fundamentally assess the amount of information contained in the trajectory and cannot ensure that the dynamic characteristics of the robotic arm are fully excited. Meanwhile, in joint stiffness identification, mainstream methods typically apply a load to the end effector of the robotic arm and calculate the Cartesian space stiffness by measuring the end-effector displacement. This is then mapped to joint space using the Jacobian matrix to estimate the joint stiffness. However, the end-effector displacement measured by this method is essentially a combined result of the joint and link flexibility deformations. Therefore, the joint stiffness obtained through this mapping is actually an equivalent joint stiffness coupled with link flexibility, rather than the true stiffness of the joint transmission components themselves, leading to inherent biases in the identification results. Furthermore, this method is offline identification and cannot capture the dynamic changes in joint stiffness during actual operation. Summary of the Invention

[0004] In view of this, the purpose of this invention is to provide an online identification method for the joint stiffness of industrial robotic arms based on the optimal excitation trajectory, which can effectively improve the identification accuracy and numerical stability, and effectively reveal the inherent laws of complex time-varying stiffness.

[0005] To achieve the above objectives, the present invention provides the following technical solution: An online method for identifying the joint stiffness of an industrial robotic arm based on an optimal excitation trajectory includes the following steps: Step 1: For a robotic arm system containing flexible joints, establish its rigid-flexible coupling dynamic model; Step 2: Based on the aforementioned dynamic model, parameterize the motor-side dynamic equations to construct a model with joint stiffness. With joint viscous damping coefficient For the linear regression model of the parameters to be identified; Step 3: Set the observation noise to independent and identically distributed Gaussian white noise, construct the conditional probability density function of the observations based on the linear regression model, and derive the parameters to be identified. Fisher Information Matrix The parsing expression; Step 4: Design the excitation trajectory And parameterize it as a finite-term Fourier series; based on the motion capability of the robotic arm joints, set the physical constraints of the angle, angular velocity and angular acceleration of the excitation trajectory; Step 5: Construct the Fisher information matrix. The scalar properties and sustainable excitation conditions are used as the objective of the optimization function. Under the condition of satisfying the physical constraints described in step five, the Fourier coefficients are solved by the optimization algorithm to generate the optimal excitation trajectory. Step Six: Drive the robotic arm to execute the optimal excitation trajectory and simultaneously acquire the motor side rotation angle. , connecting rod side angle and motor drive torque Time series data; Step 7: Input the collected data into the parameter identification algorithm constructed based on the dynamic model described in Step 1 to estimate and update the joint stiffness online in real time. The value of .

[0006] Furthermore, in step one, the robotic arm system containing flexible joints is a single-degree-of-freedom flexible joint robotic arm, and its rigid-flexible coupling dynamic model is as follows: in: , , , , , and These are the mass of the robotic arm, the length of the link, the Coulomb coefficient of friction, the viscous coefficient of friction, the joint stiffness, the joint viscous damping coefficient, and the motor drive torque, respectively. It is the acceleration due to gravity; and These are the connecting rod side rotation angle and the motor side rotation angle, respectively. The moment of inertia is the rotational inertia on the motor side. is a sign function; the superscript indicates differentiation with respect to time.

[0007] Furthermore, in step two, the motor-side dynamic equations in the rigid-flexible coupling dynamic model are rewritten in linear form: in: This represents the observed value of the motor drive torque; Indicates the moment of inertia on the motor side; Indicates the angle of twist; Indicates the speed of torsion; The linear regression model is then expressed as: in: These are the observed values; Denotes the regression matrix, and ; Let be the parameter vector, and ; The observed noise is represented by a mean of 0 and a variance of . The noise is Gaussian distributed and the noise at different times is independent of each other.

[0008] Furthermore, in step three, the Fisher information matrix The parsing expression is: in: The number of observations; For the incentive trajectory; These are intermediate parameters.

[0009] Furthermore, in step four, the excitation trajectory The finite-term Fourier series is expressed as: in: Indicates the offset; Indicates the fundamental frequency; Indicates the harmonic order; Represents the coefficient of the cosine term in the velocity expression; Represents the coefficient of the sine term in the velocity expression; The physical constraints on the angle, angular velocity, and angular acceleration of the excitation trajectory are expressed as follows: in: and These are the minimum and maximum values ​​of the excitation trajectory angle, respectively; and These are the minimum and maximum values ​​of the angular velocity of the excitation trajectory, respectively. and These are the minimum and maximum values ​​of the angular acceleration of the excitation trajectory, respectively.

[0010] Furthermore, in step five, the optimization function is: For optimization function First item Used to maximize the Fisher information matrix The determinant of the matrix; the second term Used to maximize the Fisher information matrix The trace; the third item Used to minimize the Fisher information matrix The condition number; the fourth term Used to guarantee the matrix The minimum eigenvalue is large enough to satisfy the sustainable incentivization condition, and: in: This represents the regression matrix.

[0011] Furthermore, in step seven, the parameter identification algorithm employs an unscented Kalman filter algorithm, which includes the following steps: Based on the rigid-flexible coupling dynamics model, the state space vector is obtained. and input vector Using the Euler method and discrete-time equations, the relationship between the state variables at the next time step and the state variables and inputs at the previous time step can be obtained: The observation equation is: in: Indicates the angle of twist Observation noise; In each forecast period The prediction phase is as follows: in: The dimension representing the state variable; This represents the scaling factor, which directly determines the diffusion distance of the sigma point relative to the mean. Used to encode prior distribution knowledge of state variables; The selection of the sampling points controls the distribution of the sampling points; Used for higher-order moment correction and ensuring numerical stability; The process covariance matrix is ​​a diagonal matrix. and Indicates the first Sigma points, ; This represents the estimated current state. Represents the error covariance matrix; and Indicates the first The mean weight of each Sigma point ; and Indicates the first Covariance weights of Sigma points ; express The predicted state Sigma point at any given time; express The mean of the predicted state at any given time; express The prediction covariance matrix at time step; Indicates the input quantity; Indicates process noise; In each forecast period The update phase is as follows: in: It is the Kalman gain matrix; It is the error covariance matrix; It is the prediction observation error covariance matrix; It is the state-observation cross-covariance matrix; express The error covariance matrix at time t; Represents the observation noise covariance matrix; These are actual observed values; express The observation Sigma point for time prediction; express The observed mean of the time prediction; express The predicted torsion angle Sigma point at any given time.

[0012] Furthermore, in step seven, the parameter identification algorithm employs a recursive least squares algorithm with a forgetting factor, and its recursive formula is as follows: in: This represents the gain vector, which determines the weight of the contribution of new data to parameter updates; The covariance matrix represents the uncertainty of parameter estimation and converges gradually with iteration. Represents the parameter estimation vector; For each iteration step; This represents the forgetting factor, used to adjust the weights of historical data.

[0013] The beneficial effects of this invention are as follows: This invention presents an online identification method for the joint stiffness of industrial robotic arms based on optimal excitation trajectories. Based on optimal excitation trajectory design, it achieves high-precision online identification of the joint stiffness of industrial robotic arms and has the following technical advantages: (1) Improved the accuracy and robustness of parameter identification: Traditional excitation trajectory design only focuses on the condition number of the observation matrix, while this invention uses the determinant and trace of the Fisher information matrix as optimization targets, ensuring that the excitation trajectory contains the maximum amount of parameter information from the information theory level, thus showing higher identification accuracy in algorithms such as UKF and RLS, and remaining stable in noisy environments. (2) The true stiffness of the joint is directly identified online: By establishing a rigid-flexible coupling dynamic model and parameterizing the motor side equation, the joint stiffness K and damping D can be directly identified, avoiding the mapping error introduced by the flexible coupling of the link, and the stiffness can be tracked online. (3) Enhanced the engineering applicability and reliability of the method: By parameterizing the trajectory with Fourier series and applying physical constraints, the feasibility of the excitation trajectory is guaranteed; the continuous excitation condition is introduced into the optimization function to make the identification process continuously converge.

[0014] In summary, the method of this invention provides accurate joint stiffness information for high-precision robot control, and has important application value in fields such as precision assembly and high-end manufacturing. Attached Figure Description

[0015] To make the objectives, technical solutions, and beneficial effects of this invention clearer, the following figures are provided for illustration: Figure 1 This is a flowchart of the online identification method for joint stiffness of industrial robotic arms based on optimal excitation trajectories according to the present invention. Figure 2 Schematic diagram of a single-degree-of-freedom flexible joint robotic arm; Figure 3The graphs show a comparison between the traditional excitation trajectory and the excitation trajectory of the present invention; where (a), (b), and (c) are comparison graphs of trajectory angle, angular velocity, and angular acceleration, respectively. Figure 4 The results are from the Unscented Kalman Filter (UKF) algorithm. Figure 5 The identification result is for the Recursive Least Squares (RLS) algorithm; Figure 6 Comparison of root mean square errors for different trajectories under two identification algorithms; Figure 7 This invention demonstrates the improvement in root mean square error compared to traditional methods under two identification algorithms. Detailed Implementation

[0016] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0017] This invention, based on information theory, proposes an online method for identifying the joint stiffness of industrial robotic arms based on optimal excitation trajectories. First, starting from the foundations of information theory, an analytical expression for the Fisher information matrix is ​​derived, and some properties of the Fisher information matrix are established as core criteria for excitation trajectory design, thereby effectively quantifying the information contribution of the excitation trajectory to the parameters to be identified. Second, to overcome the problem of limited identification performance due to data saturation, this invention introduces a sustainable excitation condition to ensure that the identification results are continuously updated and converged, significantly improving the accuracy and numerical stability of the parameter identification algorithm. Compared with traditional excitation trajectory design criteria, the method proposed in this invention exhibits higher identification accuracy and numerical stability in both identification algorithms, effectively revealing the inherent laws of complex time-varying stiffness. Furthermore, this design criterion possesses good universality and scalability, and can be widely applied to other fields requiring high-precision parameter identification.

[0018] Specifically, such as Figure 1 As shown in the figure, the online identification method for joint stiffness of industrial robotic arms based on optimal excitation trajectory in this embodiment includes the following steps.

[0019] Step 1: For a robotic arm system containing flexible joints, establish its rigid-flexible coupling dynamic model.

[0020] This embodiment takes a single-degree-of-freedom flexible joint robotic arm as the research object, and its simplified schematic diagram is as follows: Figure 2 As shown, the effectiveness of the proposed excitation trajectory design algorithm is verified. The rigid-flexible coupling dynamic model of the single-degree-of-freedom flexible joint manipulator is shown in the following equation: in: , , , , , and These are the mass of the robotic arm, the length of the link, the Coulomb coefficient of friction, the viscous coefficient of friction, the joint stiffness, the joint viscous damping coefficient, and the motor drive torque, respectively. It is the acceleration due to gravity; and These are the connecting rod side rotation angle and the motor side rotation angle, respectively. The moment of inertia is the rotational inertia on the motor side. is a sign function; the superscript indicates differentiation with respect to time.

[0021] Step 2: Based on the aforementioned dynamic model, and assuming the motor rotor inertia is known, the motor-side dynamic equations are parameterized to construct a model with joint stiffness. With joint viscous damping coefficient For the linear regression model of the parameters to be identified.

[0022] The motor-side dynamic equations in the rigid-flexible coupling dynamic model are rewritten in linear form: in: This represents the observed value of the motor drive torque; Indicates the moment of inertia on the motor side; Indicates the angle of twist; Indicates the speed of torsion.

[0023] The linear regression model is then expressed as: in: These are the observed values; Denotes the regression matrix, and ; Let be the parameter vector, and ; The observed noise is represented by a mean of 0 and a variance of . The noise is Gaussian distributed and the noise at different times is independent of each other.

[0024] Step 3: Set the observation noise to independent and identically distributed Gaussian white noise, construct the conditional probability density function of the observations based on the linear regression model, and derive the parameters to be identified. Fisher Information Matrix The parsing expression.

[0025] (1) Conditional probability density function: Assuming that the torque measurement noise in the regression equation is independent and identically distributed Gaussian white noise, the conditional probability density function of the torque measurement value is constructed based on this.

[0026] For a single moment Observed values The conditional probability density function is: Because the noise is independent, The joint probability density of the observations is the product of the individual probabilities, which further leads to: Taking a logarithmic transformation to convert the product into a summation, and extracting the constant term, we get: Among them: the first item It is the sum of squared residuals, and the parameter and Related to; the second item Is with parameters and An irrelevant constant that depends only on the variance of the noise.

[0027] (2) Based on the probability density function, the analytical expression of the Fisher information matrix of the parameter to be identified is calculated; then, the torque term in the regression equation is expressed as a function of the excitation trajectory parameter, thereby obtaining the Fisher information matrix completely represented by the excitation trajectory.

[0028] The Fisher information matrix is ​​introduced below to design the excitation trajectory. This matrix is ​​used to measure the effect of the observed data on the parameters. and The information content, whose elements are defined as the negative expectation of the second derivative of the log-likelihood function: in: The first derivative with respect to the parameter is as follows: Furthermore, the second derivative is obtained: Since the second derivative does not contain Its expected value depends on and noise variance The initial expression for the Fisher information matrix is ​​obtained as follows: The following solution is obtained from equation (1) using the connecting rod side dynamic equation: For ease of writing, let: Under conditions where the system's frictional characteristics are unknown, to reduce the impact of unmodeled dynamics on identification accuracy, they are uniformly modeled as observation noise. This allows us to obtain the final expression for the Fisher information matrix. Only with incentive trajectories Related, The value is a positive constant and does not affect the solution to the optimization problem. Therefore, the Fisher information matrix... The parsing expression is: in: The number of observations; For the incentive trajectory; These are intermediate parameters.

[0029] Step 4: Design the excitation trajectory And parameterize it as a finite-term Fourier series; based on the motion capability of the robotic arm joints, set the physical constraints of the angle, angular velocity and angular acceleration of the excitation trajectory.

[0030] A finite-term Fourier series is used as the parameterized form of the excitation trajectory, and physical constraints are set for the angle, angular velocity, and angular acceleration of the robotic arm joints. The Fourier series is used as the excitation trajectory, as shown in (13), and constraint boundaries for joint angle, velocity, and acceleration are set to ensure that the excitation trajectory meets the physical constraints: in: Indicates the offset; Indicates the fundamental frequency; Indicates the harmonic order; Represents the coefficient of the cosine term in the velocity expression; Represents the coefficient of the sine term in the velocity expression; and These are the minimum and maximum values ​​of the excitation trajectory angle, respectively; and These are the minimum and maximum values ​​of the angular velocity of the excitation trajectory, respectively. and These are the minimum and maximum values ​​of the angular acceleration of the excitation trajectory, respectively.

[0031] Step 5: Construct the Fisher information matrix. The scalar properties and sustainable incentive conditions are used as the objective of the optimization function. Under the condition of satisfying the physical constraints described in step five, the Fourier coefficients are solved by the optimization algorithm. The optimal excitation trajectory is generated to achieve high-precision online identification of joint stiffness.

[0032] To ensure the excitation trajectory provides sufficient information, the Fisher information matrix must not be a singular matrix. To guarantee the continuous updating of the parameters to be identified, a sustainable excitation condition is introduced: the existence of a regression matrix. This makes equation (14) true: The constraint function and optimization objective of the excitation trajectory are shown in Equation (15). Since the numerical differences between the indicators in the optimization objective are significant, in order to better iteratively optimize, a method is adopted. and By scaling the function and combining it with optimization algorithms, the coefficients of the Fourier series can be obtained: For optimization function First item Used to maximize the Fisher information matrix The determinant of the matrix; the second term Used to maximize the Fisher information matrix The trace; the third item Used to minimize the Fisher information matrix The condition number; the fourth term Used to guarantee the matrix The minimum eigenvalue is large enough to satisfy the sustainable incentive condition.

[0033] Specifically, the optimization function used in traditional excitation trajectory design is: To objectively evaluate the performance of this invention, a comparison scheme was designed in this embodiment: the traditional excitation trajectory and the optimal excitation trajectory described in this invention were compared under two identification algorithms: unscented Kalman filtering (UKF) and recursive least squares (RLS). The time-varying stiffness of this single-degree-of-freedom flexible joint manipulator is given by equation (17), with units of Nm / rad: Step Six: Drive the robotic arm to execute the optimal excitation trajectory and simultaneously acquire the motor side rotation angle. , connecting rod side angle and motor drive torque Time series data.

[0034] Step 7: Input the collected data into the parameter identification algorithm constructed based on the dynamic model described in Step 1 to estimate and update the joint stiffness online in real time. The value of .

[0035] Specifically, based on the rigid-flexible coupling dynamic model established in step one, unscented Kalman filter (UKF) and recursive least squares (RLS) identification algorithms are constructed respectively. The system data collected by the traditional excitation trajectory and the optimal trajectory in this embodiment are substituted into the above identification algorithms, and the identification accuracy of the two trajectories is comprehensively evaluated and compared using the root mean square error (RMSE) as an indicator.

[0036] To verify the effectiveness of the excitation trajectory design method proposed in this invention, this embodiment constructs two identification algorithms: Unscented Kalman Filter (UKF) and Recursive Least Squares (RLS), and compares and analyzes the traditional trajectory design method with the method of this invention.

[0037] (1) Unscented Kalman Filter (UKF) Algorithm The parameter identification algorithm employs an unscented Kalman filter algorithm, including: Based on the rigid-flexible coupling dynamics model, the state space vector is obtained. and input vector Furthermore, we obtain: in: The process noise represents the state variable, assuming it follows a mean of 0 and a variance of . The Gaussian distribution. Using the Euler method, the discrete-time equations can be used to obtain the relationship between the state variables at the next time step and the state variables and inputs at the previous time step, as shown in equation (19): in: The observation equation is given by the following formula: in: Indicates the angle of twist The observation noise is assumed to have a mean of 0 and a variance of . The algorithm uses a Gaussian distribution. The entire process is mainly divided into two stages: prediction and update.

[0038] In each forecast period The prediction phase is as follows: in: The dimension representing the state variable; This represents the scaling factor, which directly determines the diffusion distance of the sigma point relative to the mean. Used to encode prior distribution knowledge of state variables; The selection of the sampling points controls the distribution of the sampling points; Used for higher-order moment correction and ensuring numerical stability; The process covariance matrix is ​​a diagonal matrix. and Indicates the first Sigma points, ; This represents the estimated current state. Represents the error covariance matrix; and Indicates the first The mean weight of each Sigma point ; and Indicates the first Covariance weights of Sigma points ; express The predicted state Sigma point at any given time; express The mean of the predicted state at any given time; express The prediction covariance matrix at time step; Indicates the input quantity; This indicates process noise.

[0039] In each forecast period The update phase is as follows: in: It is the Kalman gain matrix; It is the error covariance matrix; It is the prediction observation error covariance matrix; It is the state-observation cross-covariance matrix; express The error covariance matrix at time t; Represents the observation noise covariance matrix; These are actual observed values; express The observation Sigma point for time prediction; express The observed mean of the time prediction; express The predicted torsion angle Sigma point at any given time.

[0040] Updated state estimate Covariance Using the initial value for the next moment, the prediction and update steps are repeated to achieve recursive filtering, so that the identified value continuously approaches the true value.

[0041] (2) Recursive Least Squares (RLS) Algorithm The parameter identification algorithm employs a recursive least squares (RLS) algorithm with a forgetting factor. Specifically, the RLS algorithm achieves online real-time calculation of parameter estimates by constructing a recursive update mechanism. The recursive formula for the RLS algorithm is: in: This represents the gain vector, which determines the weight of the contribution of new data to parameter updates; The covariance matrix represents the uncertainty of parameter estimation and converges gradually with iteration. Represents the parameter estimation vector; For each iteration step; This represents the forgetting factor, used to adjust the weights of historical data. The smaller the value, the faster the old data is forgotten, making it suitable for fast time-varying systems.

[0042] Based on the aforementioned physical constraints and optimization objective function, the traditional excitation trajectory and the optimal excitation trajectory generated by this invention were obtained using MATLAB's fmincon optimization solver. The results are comparable to, for example... Figure 3 As shown in (a), (b) and (c) in the figure.

[0043] To verify the effectiveness and robustness of the excitation trajectory proposed in this invention, the following comparative experiment was set up: Under two typical working conditions, namely no noise and Gaussian white noise, the traditional excitation trajectory and the optimal excitation trajectory of this invention were used as command signals, and two typical algorithms, namely unscented Kalman filtering (UKF) and recursive least squares (RLS), were applied simultaneously for parameter identification. Due to the time-varying stiffness, the influence of the viscous damping coefficient on the dynamics is easily masked by the time-varying stiffness, resulting in inaccurate identification by both the traditional algorithm and the algorithm of this invention under different identification algorithms, so it is not shown. Figure 4 and Figure 5 The results of UKF and RLS algorithms for identifying two different trajectories are shown in the comparison.

[0044] To quantify the recognition performance of traditional trajectories and the trajectory of this invention, a comparative analysis was conducted using the root mean square error (RMSE) metric. Figure 6-7 As shown in the figure. The results show that the trajectory generated by this invention has improved the parameter identification accuracy to varying degrees compared with the traditional trajectory, regardless of whether there is noise or noise. At the same time, it exhibits stronger robustness in noisy environments, demonstrating its significant advantages in handling uncertainty.

[0045] The above-described embodiments are merely preferred embodiments provided to fully illustrate the present invention, and the scope of protection of the present invention is not limited thereto. Equivalent substitutions or modifications made by those skilled in the art based on the present invention are all within the scope of protection of the present invention. The scope of protection of the present invention is defined by the claims.

Claims

1. An optimal excitation trajectory based online identification method for joint stiffness of industrial robot arms, characterized in that: The method comprises the following steps: Step one: for a flexible joint containing mechanical arm system, a rigid-flexible coupling dynamics model is established; Step two: based on the kinetic model, parameterize the motor side kinetic equation, construct a linear regression model with joint stiffness and joint viscous damping coefficient as to-be-identified parameters; Step three: set the observation noise as independent and identically distributed Gaussian white noise, construct the conditional probability density function of the observation value based on the linear regression model, and deduce the analytical expression of the Fisher information matrix of the to-be-identified parameter ​​ Step four: design the excitation trajectory and parameterize it as a finite Fourier series; set physical constraints on the angle, angular velocity and angular acceleration of the excitation trajectory according to the movement ability of the joints of the robot arm; Step five: construct an optimization function targeting the scalar properties of the Fisher information matrix of step five, generate the optimal excitation trajectory by solving the Fourier coefficients with an optimization algorithm under the condition of satisfying the physical constraints of step five; Step six: drive the mechanical arm to execute the optimal excitation trajectory, and synchronously collect the time sequence data of motor side rotation angle , connecting rod side rotation angle and motor driving torque ​ Step seven: input the collected data into the parameter identification algorithm based on the dynamic model constructed in step one, estimate and update the value of joint stiffness online in real time. Step seven: input the collected data into the parameter identification algorithm based on the dynamic model constructed in step one, estimate and update the value of joint stiffness online in real time.

2. The method of claim 1, wherein the method further comprises: In the step one, the flexible joint containing mechanical arm system is a single degree of freedom flexible joint mechanical arm, and the rigid-flexible coupling dynamics model is: wherein: , , , , , and are the mass of the robot arm, the link length, the Coulomb friction coefficient, the viscous friction coefficient, the joint stiffness, the joint viscous damping coefficient, the motor driving torque, respectively; is the gravitational acceleration; and are the link side angle and the motor side angle, respectively; is the motor side moment of inertia; is the sign function; the superscript dot denotes the derivation with respect to time.

3. The method of claim 1, wherein the method further comprises: In the step two, the motor side dynamics equation in the rigid-flexible coupling dynamics model is rewritten into a linear form: wherein: represents a motor driving torque observation value; represents a motor-side rotational inertia; represents a twist angle; represents a twist speed; The linear regression model is represented as: where are observations; denotes a regression matrix, and ; is a parameter vector, and ; denotes observation noise, which is subject to a Gaussian distribution with mean 0 and variance , and the noises at different time instants are mutually independent.

4. The method of claim 1, wherein: In the third step, the Fisher information matrix The analytical expression of the Fisher information matrix is wherein: is the number of observations; is the excitation trajectory; is the intermediate parameter.

5. The method of claim 1, wherein: In step four, the excitation trajectory is represented in a finite-term Fourier series form as wherein: represents an offset amount; represents a fundamental frequency; represents a harmonic number; represents a coefficient of a cosine term in a velocity expression; represents a coefficient of a sine term in a velocity expression; The physical constraints of the angle, angular velocity and angular acceleration of the excitation trajectory are represented as: wherein: and are the minimum and maximum values of the angle of the excitation trajectory, respectively; and are the minimum and maximum values of the angular velocity of the excitation trajectory, respectively; and are the minimum and maximum values of the angular acceleration of the excitation trajectory, respectively.

6. The method of claim 1, wherein: In the step five, the optimization function is: For the optimization function : the first term is for maximizing the determinant of the Fisher information matrix ; the second term is for maximizing the trace of the Fisher information matrix ; the third term is for minimizing the condition number of the Fisher information matrix ; the fourth term is for ensuring that the minimum eigenvalue of the matrix is large enough to satisfy the sustainable incentive condition: wherein: represents a regression matrix.

7. The method of claim 1, wherein: In the step seven, the parameter identification algorithm adopts an unscented Kalman filter algorithm, comprising the following steps: According to the rigid-flexible coupling dynamics model, a state space vector and an input vector The Euler method is used to discretize the time equation to obtain the relationship between the state variable at the next time and the state variable and input at the last time: The observation equation is: wherein: represents the torsion angle an observation noise; In each prediction cycle The prediction phase of each prediction cycle is: where denotes the dimension of the state variable; denotes the scaling factor, which directly determines the spread of the sigma points relative to the mean; the prior distribution knowledge of the state variable; the selection of controls the distribution of the sampling points; for high-order moment correction and numerical stability; denotes the process covariance matrix, which is a diagonal matrix; and denotes the th Sigma point, ; denotes the current state estimate value; denotes the error covariance matrix; and denotes the th Sigma point mean weight, ; and denotes the th Sigma point covariance weight, ; denotes the predicted state Sigma point at time ; denotes the predicted state mean at time ; denotes the predicted covariance matrix at time ; denotes the input quantity; denotes the process noise; At each prediction cycle The update phase of each prediction cycle is: where: is the Kalman gain matrix; is the error covariance matrix; is the predicted observation error covariance matrix; is the state-observation cross-covariance matrix; denotes the error covariance matrix at time instant; denotes the observation noise covariance matrix; is the actual observation value; denotes the predicted observation Sigma point at time instant; denotes the predicted observation mean at time instant; denotes the predicted twist angle Sigma point at time instant.

8. The method of claim 1, wherein: In the step seven, the parameter identification algorithm adopts a recursive least square algorithm with a forgetting factor, and the recursive formula is: wherein: represents a gain vector, which determines the contribution weight of new data to parameter update; represents a covariance matrix, which reflects the uncertainty of parameter estimation and gradually converges with iteration; represents a parameter estimation vector; is an iteration step; represents a forgetting factor, which is used to adjust the weight of historical data.

Images