Harmonic reducer dynamic transmission error distribution characteristic optimization method

Abstract

This invention discloses a method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer, comprising: establishing a static transmission error probability model to obtain the probability distribution of the overall static transmission error; constructing a dynamic transmission error mathematical model considering static transmission error and dynamic parameters; constructing a high-precision surrogate model of dynamic transmission error including the probability distribution of static transmission error and the range of dynamic parameters; obtaining the probability distribution of dynamic transmission error; and dynamically adjusting the range of dynamic parameters based on a particle swarm optimization strategy to find the parameter range that optimizes the dynamic transmission error distribution. By introducing static transmission error into the dynamic model of the harmonic reducer, considering the influence of the processing and assembly of the harmonic reducer on the transmission error, a probability-range hybrid uncertainty model is used to mathematically describe the probability distribution characteristics of static transmission error and the range of dynamic parameters, and a Chebyshev polynomial method is used to construct an approximate model of dynamic transmission error.

Description

Technical Field

[0001] This invention relates to a method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer, belonging to the technical field of harmonic reducers. Background Technology

[0002] In high-precision transmission fields such as industrial robots and humanoid robots, harmonic reducers, with their advantages of large transmission ratio, small size, high load-bearing capacity, and high transmission accuracy, have become core components of joint drive mechanisms. Their transmission performance directly determines the overall motion accuracy and stability of the robot. However, under complex and variable working conditions and during long-term operation, harmonic reducers often face the problem of increasing dynamic transmission errors. The generation and evolution mechanism of this error is complex, and accurate calculation and effective control are key to improving the reliability and accuracy of the transmission system.

[0003] Currently, research on dynamic transmission errors in harmonic reducers largely relies on deterministic dynamic models for analysis. These methods typically assume that system parameters (such as stiffness and damping) are fixed values ​​and often neglect the influence of static transmission errors during modeling. However, in actual manufacturing, assembly, and operation, the geometric transmission errors characterizing static accuracy, as well as the dynamic parameters reflecting the system's dynamic characteristics, are not constant but are affected by various factors such as machining tolerances, wear, load fluctuations, and temperature changes, exhibiting multi-source uncertainties that fluctuate within a certain range. Existing deterministic analysis methods that ignore parameter uncertainties and static errors cannot accurately reflect the actual distribution characteristics of dynamic transmission errors under complex operating conditions, leading to deviations between predicted results and engineering realities. Therefore, they cannot provide a reliable basis for subsequent accuracy control and design optimization. Summary of the Invention

[0004] Purpose of the invention: To address the shortcomings of existing technologies, this invention provides a method for optimizing the dynamic transmission error distribution characteristics of harmonic reducers. By fully considering the impact of static transmission errors (including machining and assembly) and uncertainties in dynamic parameters such as system stiffness and damping on dynamic transmission errors, this invention reflects the probability distribution characteristics of dynamic transmission errors of harmonic reducers under real complex working conditions to the greatest extent. Furthermore, based on a particle swarm optimization strategy, it dynamically adjusts the interval positions of uncertain dynamic parameters, thereby optimizing the probability distribution of dynamic transmission errors.

[0005] Technical solution: A method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer, comprising the following steps:

[0006] S1. Based on the measured statistical parameters of the error sources in the processing and assembly of the harmonic reducer, a static transmission error probability model is established through Monte Carlo (MC) simulation to obtain the probability distribution of the static transmission error of the whole machine.

[0007] S2. Based on the probability distribution of the static transmission error of the whole machine in S1, a dynamic transmission error mathematical model considering static transmission error and dynamic parameters is constructed using nonlinear differential dynamic equations.

[0008] S3. Using the Chebyshev polynomial expansion method, a high-precision surrogate model of dynamic transmission error is constructed based on the mathematical model of S2, which includes the probability distribution of static transmission error and the range of dynamic parameters.

[0009] S4. Using Monte Carlo sampling, the probability distribution of static transmission error and the range of dynamic parameters in the high-precision surrogate model of S3 are sampled and substituted into the model calculation in S3 to obtain the probability distribution of dynamic transmission error.

[0010] S5. Based on the particle swarm optimization strategy, the dynamic parameter range in S3 is dynamically adjusted to find the parameter range that makes the dynamic transmission error distribution optimal.

[0011] In a preferred embodiment, S1 specifically comprises:

[0012] The manufacturing and assembly errors of the components of the harmonic reducer include the cumulative error ΔF of the flexspline tooth pitch. p1 Flex gear tooth tangential comprehensive error Δf i1 The clearance ΔE between the flexible wheel and the mounting hole 13 Radial runout ΔE between the flexible wheel and the mounting hole 14 Cumulative error ΔF of the tooth pitch of the rigid wheel p1 Δf, the combined tangential error of the rigid gear teeth i2 Radial runout ΔE between the rigid wheel and the mounting hole 23 The clearance ΔE between the rigid wheel and the mounting hole 24 Radial runout ΔE of wave generator 31 The gap ΔE between the wave generator and the input shaft 32 Wave generator profile error ΔE 33 The gap ΔE between the wave generator and the flexible bearing 34 The clearance ΔE between the flexible bearing and the flexure wheel 41 ;

[0013] The above errors are categorized into four main types based on their period:

[0014] Fixed eccentricity error e m Cumulative deviation e of the tooth pitch of the rigid wheel and the flexible wheel p The resulting motion error has a frequency of 2ω;

[0015] Eccentricity error e as the flexible wheel rotates g The generated frequency is 2Z c / Z f The motion error of ω;

[0016] The eccentricity error e that rotates with the wave generator w The resulting motion error has a frequency of ω;

[0017] Small periodic error e caused by machining errors of rigid wheel and flexible wheel h The generated frequency is 2Z c The motion error of ω;

[0018] The total motion error of the harmonic reducer is expressed as:

[0019] (1)

[0020] In the formula, e m For fixed eccentricity error; e p e represents the cumulative pitch deviation of the rigid and flexible gears; g e is the eccentricity error caused by the rotation along with the flexible wheel; w e represents the eccentricity error caused by the wave generator rotating with it. h Small periodic errors caused by machining errors of rigid and flexible wheels; φ m The initial phase angle of the fixed eccentricity error; φ g φ is the initial phase angle of the eccentricity error due to the rotation along with the flexspline. w φ is the initial phase angle of the eccentricity error that rotates with the wave generator. p φ is the initial phase angle of the cumulative pitch deviation of the rigid and flexible gears; h Z represents the initial phase angle of the small-period error caused by the machining errors of the rigid wheel and flexible wheel. c Z represents the number of teeth on the rigid wheel. f ω is the number of teeth on the flexible gear; ω is the angular velocity of the wave generator; α n t is the pressure angle at the pitch circle of the flexible wheel; t is the time parameter.

[0021] When considering the error averaging effect of multi-tooth meshing, the modified formula for the static transmission error of the harmonic reducer is expressed as:

[0022] (2)

[0023] In the formula, To account for the static transmission error when considering the averaging effect of multi-tooth meshing error; K b z is the error influence coefficient of multi-tooth meshing transmission; t The number of teeth engaged simultaneously; The total motion error is represented by d, where d is the pitch circle diameter of the rigid wheel.

[0024] L harmonic reducer prototypes were selected for eccentricity error measurement; a coordinate measuring machine was used to measure manufacturing and assembly errors, and a gear measuring center was used to measure tooth profile errors; each error item of each prototype was measured 10 times repeatedly, and the results were verified by... Calculate and take the average; where, The average value of 10 repeated measurements for each error item of each prototype; For each error term of each prototype, perform the l-th measurement; using Calculate the mean value μ of each error index for all prototypes. j and standard deviation σ j ;

[0025] Where e i,j Let μ be the j-th error of the i-th prototype; j σ is the mean of the j-th error of all prototypes; j Let j be the standard deviation of the j-th error of all prototypes;

[0026] Based on the measured statistical parameters (μ) of each error source j , σ j 2), using the Monte Carlo method to generate 1000 independent normally distributed random samples; then these sample values ​​are substituted into formula (1) and formula (2) for batch calculation, and finally the probability distribution (μ, σ²) of the static transmission error of the whole machine is obtained through statistical analysis.

[0027] In a preferred embodiment, S2 specifically comprises:

[0028] Define the kinetic energy and potential energy of the harmonic drive system as follows:

[0029] (3)

[0030] (4)

[0031] In the formula, J in J is the moment of inertia at the input end. out θ is the moment of inertia at the output end. in θ is the input angle of the wave generator. out The output angle of the flexible gear is denoted by K; the equivalent torsional stiffness is denoted by N; and the transmission ratio is denoted by N.

[0032] The system's Lagrangian function is:

[0033] (5)

[0034] The frictional loss of the system is described using the Rayleigh dissipation function:

[0035] (6)

[0036] In the formula, B in B is the input damping coefficient; out B is the output damping coefficient; fw B is the deformation damping coefficient of the cup-shaped flexible wheel; fc is the damping coefficient of the meshing of the rigid wheel and the flexible wheel;

[0037] Therefore, the dynamic equation of the system is:

[0038] (7)

[0039] In the formula, T in Input torque;

[0040] Substituting equations (3), (4), (5), and (6) into equation (7), we obtain the second-order differential dynamic equations of the harmonic reducer as follows:

[0041] (8)

[0042] The mathematical model of the dynamic transmission error of the harmonic reducer is expressed as follows:

[0043] (9)

[0044] In a preferred embodiment, S3 specifically includes:

[0045] Constructing n-dimensional a-order interpolation points:

[0046] , (10)

[0047] Where p is the order of the Chebyshev polynomial interpolation point, p = a + 1;

[0048] Calculate the n-dimensional Chebyshev polynomial series of order a:

[0049] (11)

[0050] In the formula, Let x be a point in n-dimensional space.

[0051] The output response at the sampling point is obtained by numerically solving the mathematical model of dynamic transmission error. Construct the coefficients of the Chebyshev polynomial:

[0052] (12)

[0053] In the formula, , Let be the order of the n-dimensional Chebyshev polynomial.

[0054] Define the uncertain parameters of the harmonic reducer as follows: , corresponding to θ s , K, J in B in B fw J out B out B fc The uncertain parameter is expressed in interval form as follows: ; Indicates the lower bound of the parameter interval; Indicates the upper bound of the parameter interval; Standard interval vectors can be obtained through linear transformation. express:

[0055] (13)

[0056] The Chebyshev polynomial approximation model for dynamic transmission error is established as follows:

[0057] (14)

[0058] In the formula, h is the number of zero elements in b.

[0059] In a preferred embodiment, S4 specifically comprises:

[0060] Monte Carlo random sampling is employed: for the probabilistic parameter ξ1, the static transmission error probability distribution (μ, σ²) calculated in S1 is used. For interval type parameter ξ 2:n In the hypercube [-1,1] n-1 Internal joint sampling, And through linear transformation Mapping to actual interval ;

[0061] Construct the input matrix as

[0062] (15)

[0063] Substituting V into the Chebyshev polynomial approximation model yields the probability distribution of dynamic transmission error.

[0064] Preferably, step S5 includes the following steps:

[0065] S501. Parameter definition and initialization, constructing sub-intervals of parameter compensation quantities;

[0066] S502. Establish a Chebyshev approximation model, predict dynamic transmission errors, and calculate the fitness function;

[0067] S503. Update and optimize based on particle swarm optimization strategy;

[0068] S504, Iteration Termination and Result Output.

[0069] Preferably, step S501 includes the following steps:

[0070] Parameter definition and initialization: The dynamic uncertain parameters of the harmonic reducer are as follows Corresponding to K and J respectively in B in B fc The uncertainty interval for each parameter is set to its nominal value ξ. j,nom The range is ±20% of the center, that is, the entire parameter range is ;

[0071] Initialize the particle swarm algorithm, setting the particle swarm size to N. p The maximum number of iterations is k max Define the position vector of the i-th particle at the k-th iteration as: , where u j For parameter ξ j The normalized compensation amount (dimensionless); the target dynamic transmission error distribution is set as a truncated normal distribution, and the mean of its complete distribution is... The standard deviation is The cutoff interval is [μ target -3σ target , μ target +3σ target ], denoted as ;

[0072] Constructing subintervals of parameter compensation: in A reference point ξ is randomly selected within the area. 0,j Construct the compensation subinterval as follows And ensure that the sub-interval lies entirely within the entire interval, i.e., satisfy... and ;in For parameter ξ j The normalized compensation amount (dimensionless) at k iterations; the width of the subinterval is set to... .

[0073] Preferably, step S502 includes the following steps:

[0074] Dynamic uncertain parameters according to the entire range The remaining parameters are sampled according to their nominal values, and a Chebyshev polynomial approximation model is established according to step S3; MC sampling is applied, and the dynamic uncertainty parameters are based on the compensation sub-interval corresponding to the current particle. The remaining parameters are still sampled according to their nominal values ​​to obtain the predicted dynamic transmission error distribution, DTE. pred,k ; Calculate its mean μ pred,k and standard deviation σ pred,k ;

[0075] Define the fitness function as DTE pred,k With DTE target The weighted absolute error, i.e.

[0076] (16)

[0077] Where, ω μ Mean term The weighting coefficient, ω σ Standard deviation term The weighting coefficients.

[0078] Preferably, step S503 includes the following steps:

[0079] Update the individual's historical best position to That is, the position of particle i with the minimum fitness up to the current iteration k; update the global historical best position as follows. That is, the position with the minimum fitness among all the historical best positions of individual particles; update the velocity and position of each particle:

[0080] , (17)

[0081] in, Let be the velocity of particle i at iteration k+1; Let be the position of particle i at iteration k+1; w be the inertia weight; c1 be the learning factor for updating particle velocity using the individual's historical best position; c2 be the learning factor for updating particle velocity using the global historical best position; r1 be the random number for updating particle velocity using the individual's historical best position; and r2 be the random number for updating particle velocity using the global historical best position.

[0082] Preferably, step S504 includes the following steps:

[0083] When the number of iterations k reaches the preset maximum value k max When the time is reached, the optimization process terminates; the optimal parameter compensation amount is output as follows: Therefore, the optimal subinterval is Ultimately, the optimal predicted dynamic transmission error distribution (DTE) is obtained. pred,opt And its statistical characteristics.

[0084] Beneficial Effects: This invention incorporates static transmission error into the dynamic model of a harmonic reducer, considering the influence of harmonic reducer manufacturing and assembly on transmission error. Furthermore, it employs a probability-interval hybrid uncertainty model to mathematically describe the probability distribution characteristics of static transmission error and the interval of dynamic parameters. A Chebyshev polynomial method is used to construct an approximate model of dynamic transmission error, and Monte Carlo simulation is used for sampling calculations. A particle swarm optimization strategy is then employed to optimize the dynamic transmission error distribution characteristics. This solves the problems of difficulty in solving the dynamic transmission error of harmonic reducers under complex operating conditions and the difficulty in capturing the probability distribution characteristics. It also provides a strategy for optimizing the dynamic transmission error distribution characteristics, offering technical support for the precision design and optimization of harmonic reducers. Attached Figure Description

[0085] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort.

[0086] Figure 1 This is a schematic diagram of the method of the present invention;

[0087] Figure 2 This is a schematic diagram of the harmonic reducer structure of the present invention.

[0088] Figure 3 This is a diagram of the dynamic model of the harmonic reducer of the present invention;

[0089] Figure 4 This is a flowchart of the dynamic transmission error approximation modeling and distribution characteristic solution scheme of the present invention;

[0090] Figure 5 This is a flowchart of the dynamic transmission error distribution optimization scheme of the present invention. Detailed Implementation

[0091] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0092] In the description of this invention, it should be understood that the terms "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0093] In this invention, unless otherwise explicitly specified and limited, "above" or "below" the second feature can include direct contact between the first and second features, or contact between the first and second features through another feature between them. Furthermore, "above," "over," and "on top" of the second feature includes the first feature directly above or diagonally above the second feature, or simply indicates that the first feature is at a higher horizontal level than the second feature. "Below," "below," and "under" the second feature includes the first feature directly below or diagonally below the second feature, or simply indicates that the first feature is at a lower horizontal level than the second feature.

[0094] like Figure 1 As shown, a method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer includes the following steps:

[0095] S1. Based on the measured statistical parameters of the error sources in the processing and assembly of the harmonic reducer, a static transmission error probability model is established through Monte Carlo (MC) simulation to obtain the probability distribution of the static transmission error of the whole machine.

[0096] Specifically, S1 is:

[0097] like Figure 2 The diagram shows the structure of an existing harmonic reducer, including a wave generator (1), a rigid wheel (2), a flexible wheel (3), a small end cover (4), a large end cover (5), a deep groove ball bearing (6), a flexible bearing (7), and a cross-roller bearing (8). The manufacturing and assembly errors of these components include the cumulative error ΔF of the flexible wheel tooth pitch. p1 Flex gear tooth tangential comprehensive error Δf i1 The clearance ΔE between the flexible wheel and the mounting hole 13 Radial runout ΔE between the flexible wheel and the mounting hole 14 Cumulative error ΔF of the tooth pitch of the rigid wheel p1 Δf, the combined tangential error of the rigid gear teeth i2 Radial runout ΔE between the rigid wheel and the mounting hole 23 The clearance ΔE between the rigid wheel and the mounting hole 24 Radial runout ΔE of wave generator 31 The gap ΔE between the wave generator and the input shaft 32 Wave generator profile error ΔE33 The gap ΔE between the wave generator and the flexible bearing 34 The clearance ΔE between the flexible bearing and the flexure wheel 41 ;

[0098] The above errors are categorized into four main types based on their period:

[0099] Fixed eccentricity error e m Cumulative deviation e of the tooth pitch of the rigid wheel and the flexible wheel p The resulting motion error has a frequency of 2ω;

[0100] Eccentricity error e as the flexible wheel rotates g The generated frequency is 2Z c / Z f The motion error of ω;

[0101] The eccentricity error e that rotates with the wave generator w The resulting motion error has a frequency of ω;

[0102] Small periodic error e caused by machining errors of rigid wheel and flexible wheel h The generated frequency is 2Z c The motion error of ω;

[0103] The total motion error of the harmonic reducer is expressed as:

[0104] (1)

[0105] In the formula, e m For fixed eccentricity error; e p e represents the cumulative pitch deviation of the rigid and flexible gears; g e is the eccentricity error caused by the rotation along with the flexible wheel; w e represents the eccentricity error caused by the wave generator rotating with it. h Small periodic errors caused by machining errors of rigid and flexible wheels; φ m The initial phase angle of the fixed eccentricity error; φ g φ is the initial phase angle of the eccentricity error due to the rotation along with the flexspline. w φ is the initial phase angle of the eccentricity error that rotates with the wave generator. p φ is the initial phase angle of the cumulative pitch deviation of the rigid and flexible gears; h Z represents the initial phase angle of the small-period error caused by the machining errors of the rigid wheel and flexible wheel. c Z represents the number of teeth on the rigid wheel. f ω is the number of teeth on the flexible gear; ω is the angular velocity of the wave generator; α nThe pressure angle at the flexure pitch circle is denoted by t, which is a time parameter. When considering the error averaging effect of multi-tooth meshing, the modified static transmission error formula for the harmonic reducer is expressed as:

[0106] (2)

[0107] In the formula, To account for the static transmission error when considering the averaging effect of multi-tooth meshing error; K b z is the error influence coefficient of multi-tooth meshing transmission; t The number of teeth engaged simultaneously; The total motion error is represented by d, where d is the pitch circle diameter of the rigid wheel.

[0108] L harmonic reducer prototypes were selected for eccentricity error measurement; a coordinate measuring machine (accuracy 0.3 μm) was used to measure manufacturing and assembly errors, and a gear measuring center (accuracy 1 μm / m) was used to measure tooth profile errors; each error item of each prototype was measured 10 times repeatedly, and the results were verified by... Calculate and take the average; where, The average value of 10 repeated measurements for each error item of each prototype; For each error term of each prototype, perform the l-th measurement; using Calculate the mean value μ of each error index for all prototypes. j and standard deviation σ j ; where e i,j Let μ be the j-th error of the i-th prototype; j σ is the mean of the j-th error of all prototypes; j Let be the standard deviation of the j-th error of all prototypes. Based on the measured statistical parameters (μ) of each error source. j , σ j 2), using the Monte Carlo method to generate 1000 independent normally distributed random samples; then these sample values ​​are substituted into formula (1) and formula (2) for batch calculation, and finally the probability distribution (μ, σ²) of the static transmission error of the whole machine is obtained through statistical analysis.

[0109] S2. Based on the probability distribution of the static transmission error of the whole machine in S1, a dynamic transmission error mathematical model considering static transmission error and dynamic parameters is constructed using nonlinear differential dynamic equations.

[0110] Specifically, S2 is:

[0111] like Figure 3 The diagram shows the physical model of a harmonic reducer, which includes a wave generator 1, a rigid wheel 2, a flexible wheel 3, a flexible bearing 7, and a motor 9. According to... Figure 3The dynamic relationships between the components and the physical model of the harmonic reducer are defined as follows: the kinetic energy and potential energy of the harmonic drive system are defined as follows:

[0112] (3)

[0113] (4)

[0114] In the formula, J in J is the moment of inertia at the input end. out θ is the moment of inertia at the output end. in θ is the input angle of the wave generator. out The output angle of the flexible gear is denoted by K; the equivalent torsional stiffness is denoted by N; and the transmission ratio is denoted by N.

[0115] The system's Lagrangian function is:

[0116] (5)

[0117] The frictional loss of the system is described using the Rayleigh dissipation function:

[0118] (6)

[0119] In the formula, B in B is the input damping coefficient; out B is the output damping coefficient; fw B is the deformation damping coefficient of the cup-shaped flexible wheel; fc is the damping coefficient of the meshing of the rigid wheel and the flexible wheel;

[0120] Therefore, the dynamic equation of the system is:

[0121] (7)

[0122] In the formula, T in Input torque;

[0123] Substituting equations (3), (4), (5), and (6) into equation (7), we obtain the second-order differential dynamic equations of the harmonic reducer as follows:

[0124] (8)

[0125] The mathematical model of the dynamic transmission error of the harmonic reducer is expressed as follows:

[0126] (9)

[0127] S3. Using the Chebyshev polynomial expansion method, a high-precision surrogate model of dynamic transmission error is constructed based on the mathematical model of S2, which includes the probability distribution of static transmission error and the range of dynamic parameters.

[0128] Specifically, S3 is:

[0129] like Figure 4 As shown, construct n-dimensional a-order interpolation points:

[0130] , (10)

[0131] Where p is the order of the Chebyshev polynomial interpolation point, p = a + 1;

[0132] Calculate the n-dimensional Chebyshev polynomial series of order a:

[0133] (11)

[0134] In the formula, Let x be a point in n-dimensional space.

[0135] The output response at the sampling point is obtained by numerically solving the mathematical model of dynamic transmission error. Construct the coefficients of the Chebyshev polynomial:

[0136] (12)

[0137] In the formula, , Let be the order of the n-dimensional Chebyshev polynomial.

[0138] Define the uncertain parameters of the harmonic reducer as follows: , corresponding to θ s , K, J in B in B fw J out B out B fc The uncertain parameter is expressed in interval form as follows: ; Indicates the lower bound of the parameter interval; Indicates the upper bound of the parameter interval; Standard interval vectors can be obtained through linear transformation. express:

[0139] (13)

[0140] The Chebyshev polynomial approximation model for dynamic transmission error is established as follows:

[0141] (14)

[0142] In the formula, h is the number of zero elements in b.

[0143] S4. Using Monte Carlo sampling, the probability distribution of static transmission error and the range of dynamic parameters in the high-precision surrogate model of S3 are sampled and substituted into the model calculation in S3 to obtain the probability distribution of dynamic transmission error.

[0144] Specifically, S4 is:

[0145] Monte Carlo random sampling is employed: for the probabilistic parameter ξ1, the static transmission error probability distribution (μ, σ²) calculated in S1 is used. For interval type parameter ξ 2:n In the hypercube [-1,1] n-1 Internal joint sampling, And through linear transformation Mapping to actual interval ;

[0146] Construct the input matrix as

[0147] (15)

[0148] Substituting V into the Chebyshev polynomial approximation model yields the probability distribution of dynamic transmission error.

[0149] like Figure 5 As shown, S5 dynamically adjusts the dynamic parameter range in S3 based on the particle swarm optimization strategy, thereby finding the parameter range that makes the dynamic transmission error distribution optimal.

[0150] S501. Parameter definition and initialization, constructing sub-intervals of parameter compensation quantities;

[0151] Parameter definition and initialization: The dynamic uncertain parameters of the harmonic reducer are as follows Corresponding to K and J respectively in B in B fc The uncertainty interval for each parameter is set to its nominal value ξ. j,nom The range is ±20% of the center, that is, the entire parameter range is ;

[0152] Initialize the particle swarm algorithm, setting the particle swarm size to N. p The maximum number of iterations is k max Define the position vector of the i-th particle at the k-th iteration as: , where u j For parameter ξ j The normalized compensation amount (dimensionless); the target dynamic transmission error distribution is set as a truncated normal distribution, and the mean of its complete distribution is... The standard deviation is The cutoff interval is [μ target -3σ target , μ target +3σ target ], denoted as ;

[0153] Constructing subintervals of parameter compensation: in A reference point ξ is randomly selected within the area. 0,j Construct the compensation subinterval as follows And ensure that the sub-interval lies entirely within the entire interval, i.e., satisfy... and ;in For parameter ξ j The normalized compensation amount (dimensionless) at k iterations; the width of the subinterval is set to... .

[0154] S502. Establish a Chebyshev approximation model, predict dynamic transmission errors, and calculate the fitness function;

[0155] Dynamic uncertain parameters according to the entire range The remaining parameters are sampled according to their nominal values, and a Chebyshev polynomial approximation model is established according to step S3; MC sampling is applied, and the dynamic uncertainty parameters are based on the compensation sub-interval corresponding to the current particle. The remaining parameters are still sampled according to their nominal values ​​to obtain the predicted dynamic transmission error distribution, DTE. pred,k ; and according to Calculate its mean μ pred,k and standard deviation σ pred,k ;in The DTE during the m-th sampling in the k-th iteration. pred,k value.

[0156] Define the fitness function as DTE pred,k With DTE target The weighted absolute error, i.e.

[0157] (16)

[0158] Where, ω μ Mean term The weighting coefficient, ω σ Standard deviation term The weighting coefficients.

[0159] S503. Update and optimize based on particle swarm optimization strategy;

[0160] Update the individual's historical best position to That is, the position of particle i with the minimum fitness up to the current iteration k; update the global historical best position as follows. That is, the position with the minimum fitness among all the historical best positions of individual particles; update the velocity and position of each particle:

[0161] , (17)

[0162] in, Let be the velocity of particle i at iteration k+1; Let be the position of particle i at iteration k+1; w be the inertia weight; c1 be the learning factor for updating particle velocity using the individual's historical best position; c2 be the learning factor for updating particle velocity using the global historical best position; r1 be the random number for updating particle velocity using the individual's historical best position; and r2 be the random number for updating particle velocity using the global historical best position.

[0163] S504, Iteration Termination and Result Output;

[0164] When the number of iterations k reaches the preset maximum value k max When the time is reached, the optimization process terminates; the optimal parameter compensation amount is output as follows: Therefore, the optimal subinterval is Ultimately, the optimal predicted dynamic transmission error distribution (DTE) is obtained. pred,opt And its statistical characteristics.

[0165] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the apparatus disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple; relevant parts can be referred to the method section.

[0166] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer, characterized in that: Includes the following steps: S1. Based on the measured statistical parameters of the error sources in the processing and assembly of the harmonic reducer, a static transmission error probability model is established through Monte Carlo (MC) simulation to obtain the probability distribution of the static transmission error of the whole machine. S2. Based on the probability distribution of the static transmission error of the whole machine in S1, a dynamic transmission error mathematical model considering static transmission error and dynamic parameters is constructed using nonlinear differential dynamic equations. S3. Using the Chebyshev polynomial expansion method, a high-precision surrogate model of dynamic transmission error is constructed based on the mathematical model of S2, which includes the probability distribution of static transmission error and the range of dynamic parameters. S4. Using Monte Carlo sampling, the probability distribution of static transmission error and the range of dynamic parameters in the high-precision surrogate model of S3 are sampled and substituted into the model calculation in S3 to obtain the probability distribution of dynamic transmission error. S5. Based on the particle swarm optimization strategy, the dynamic parameter range in S3 is dynamically adjusted to find the parameter range that makes the dynamic transmission error distribution optimal.

2. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 1, characterized in that: Specifically, S1 is: The manufacturing and assembly errors of each component of the harmonic reducer include the cumulative error ΔF of the flexspline tooth pitch. p1 Flex gear tooth tangential comprehensive error Δf i1 The clearance ΔE between the flexible wheel and the mounting hole 13 Radial runout ΔE between the flexible wheel and the mounting hole 14 Cumulative error ΔF of the tooth pitch of the rigid wheel p1 Δf, the combined tangential error of the rigid gear teeth i2 Radial runout ΔE between the rigid wheel and the mounting hole 23 The clearance ΔE between the rigid wheel and the mounting hole 24 Radial runout ΔE of wave generator 31 The gap ΔE between the wave generator and the input shaft 32 Wave generator profile error ΔE 33 The gap ΔE between the wave generator and the flexible bearing 34 The clearance ΔE between the flexible bearing and the flexure wheel 41 ; The above errors are categorized into four main types based on their period: Fixed eccentricity error e m Cumulative deviation e of the tooth pitch of the rigid wheel and the flexible wheel p The resulting motion error has a frequency of 2ω; Eccentricity error e as the flexible wheel rotates g The generated frequency is 2Z c / Z f The motion error of ω; The eccentricity error e that rotates with the wave generator w The resulting motion error has a frequency of ω; Small periodic error e caused by machining errors of rigid wheel and flexible wheel h The generated frequency is 2Z c The motion error of ω; The total motion error of the harmonic reducer is expressed as: (1) In the formula, e m For fixed eccentricity error; e p e represents the cumulative pitch deviation of the rigid and flexible gears; g e is the eccentricity error caused by the rotation along with the flexible wheel; w e represents the eccentricity error caused by the wave generator rotating with it. h Small periodic errors caused by machining errors of rigid and flexible wheels; φ m The initial phase angle of the fixed eccentricity error; φ g φ is the initial phase angle of the eccentricity error due to the rotation along with the flexspline. w φ is the initial phase angle of the eccentricity error that rotates with the wave generator. p φ is the initial phase angle of the cumulative pitch deviation of the rigid and flexible gears; h Z represents the initial phase angle of the small-period error caused by the machining errors of the rigid wheel and flexible wheel. c Z represents the number of teeth on the rigid wheel. f ω is the number of teeth on the flexible gear; ω is the angular velocity of the wave generator; α n The pressure angle at the pitch circle of the flexible gear; t is the time parameter; When considering the error averaging effect of multi-tooth meshing, the modified formula for the static transmission error of the harmonic reducer is expressed as: (2) In the formula, To account for the static transmission error when considering the averaging effect of multi-tooth meshing error; K b z is the error influence coefficient of multi-tooth meshing transmission; t The number of teeth engaged simultaneously; The total motion error is represented by d, where d is the pitch circle diameter of the rigid wheel. L harmonic reducer prototypes were selected for eccentricity error measurement; a coordinate measuring machine was used to measure manufacturing and assembly errors, and a gear measuring center was used to measure tooth profile errors; each error item of each prototype was measured 10 times repeatedly, and the results were verified by... Calculate and take the average; where, The average value of 10 repeated measurements for each error item of each prototype; For each error term of each prototype, perform the l-th measurement; using Calculate the mean value μ of each error index for all prototypes. j and standard deviation σ j ; Where e i,j Let μ be the j-th error of the i-th prototype; j σ is the mean of the j-th error of all prototypes; j Let j be the standard deviation of the j-th error of all prototypes; Based on the measured statistical parameters (μ) of each error source j , σ j 2), using the Monte Carlo method to generate 1000 independent normally distributed random samples; then these sample values ​​are substituted into formula (1) and formula (2) for batch calculation, and finally the probability distribution (μ, σ²) of the static transmission error of the whole machine is obtained through statistical analysis.

3. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 1, characterized in that: Specifically, S2 is: Define the kinetic energy and potential energy of the harmonic drive system as follows: (3) (4) In the formula, J in J is the moment of inertia at the input end. out θ is the moment of inertia at the output end. in θ is the input angle of the wave generator. out The output angle of the flexible gear is denoted by K; the equivalent torsional stiffness is denoted by N; and the transmission ratio is denoted by N. The system's Lagrangian function is: (5) The frictional loss of the system is described using the Rayleigh dissipation function: (6) In the formula, B in B is the input damping coefficient; out B is the output damping coefficient; fw B is the deformation damping coefficient of the cup-shaped flexible wheel; fc is the damping coefficient of the meshing of the rigid wheel and the flexible wheel; Therefore, the dynamic equation of the system is: (7) In the formula, T in Input torque; Substituting equations (3), (4), (5), and (6) into equation (7), we obtain the second-order differential dynamic equations of the harmonic reducer as follows: (8) The mathematical model of the dynamic transmission error of the harmonic reducer is expressed as follows: (9)。 4. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 1, characterized in that: Specifically, S3 is: Constructing n-dimensional a-order interpolation points: , ）(10) Where p is the order of the Chebyshev polynomial interpolation point, p = a + 1; Calculate the n-dimensional Chebyshev polynomial series of order a: (11) In the formula, Let x be a point in n-dimensional space. ; The output response at the sampling point is obtained by numerically solving the mathematical model of dynamic transmission error. Construct the coefficients of the Chebyshev polynomial: (12) In the formula, , Let be the order of the n-dimensional Chebyshev polynomial; Define the uncertain parameters of the harmonic reducer as follows: , corresponding to θ s , K, J in B in B fw J out B out B fc The uncertain parameter is expressed in interval form as follows: ; Indicates the lower bound of the parameter interval; Indicates the upper bound of the parameter interval; Standard interval vectors can be obtained through linear transformation. express: (13) The Chebyshev polynomial approximation model for dynamic transmission error is established as follows: (14) In the formula, h is the number of zero elements in b.

5. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 1, characterized in that: Specifically, S4 is: Monte Carlo random sampling is employed: for the probabilistic parameter ξ1, the static transmission error probability distribution (μ, σ²) calculated in S1 is used. For interval type parameter ξ 2:n In the hypercube [-1,1] n-1 Internal joint sampling, And through linear transformation Mapping to actual interval ; Construct the input matrix as (15) Substituting V into the Chebyshev polynomial approximation model yields the probability distribution of dynamic transmission error.

6. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 1, characterized in that: S5 includes the following steps: S501. Parameter definition and initialization, constructing sub-intervals of parameter compensation quantities; S502. Establish a Chebyshev approximation model, predict dynamic transmission errors, and calculate the fitness function; S503. Update and optimize based on particle swarm optimization strategy; S504, Iteration Termination and Result Output.

7. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 6, characterized in that: S501 includes the following steps: Parameter definition and initialization: The dynamic uncertain parameters of the harmonic reducer are as follows Corresponding to K and J respectively in B in B fc The uncertainty interval for each parameter is set to its nominal value ξ. j,nom The range is ±20% of the center, that is, the entire parameter range is ; Initialize the particle swarm algorithm, setting the particle swarm size to N. p The maximum number of iterations is k max Define the position vector of the i-th particle at the k-th iteration as: , where u j For parameter ξ j The normalized compensation quantity is dimensionless; the target dynamic transmission error distribution is set as a truncated normal distribution, and the mean of its complete distribution is... The standard deviation is The cutoff interval is [μ target -3σ target , μ target +3σ target ], denoted as ; Constructing subintervals of parameter compensation: in A reference point ξ is randomly selected within the area. 0,j Construct the compensation subinterval as follows And ensure that the sub-interval lies entirely within the entire interval, i.e., satisfy... and ;in For parameter ξ j The normalized compensation amount at k iterations is dimensionless; the width of the subinterval is set to... .

8. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 6, characterized in that: S502 includes the following steps: Dynamic uncertain parameters according to the entire range The remaining parameters are sampled according to their nominal values, and a Chebyshev polynomial approximation model is established according to step S3; MC sampling is applied, and the dynamic uncertainty parameters are based on the compensation sub-interval corresponding to the current particle. The remaining parameters are still sampled according to their nominal values ​​to obtain the predicted dynamic transmission error distribution, DTE. pred,k ; Calculate its mean μ pred,k and standard deviation σ pred,k ; Define the fitness function as DTE pred,k With DTE target The weighted absolute error, i.e. (16) Where, ω μ Mean term The weighting coefficient, ω σ Standard deviation term The weighting coefficients.

9. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 6, characterized in that: S503 includes the following steps: Update the individual's historical best position to That is, the position of particle i with the minimum fitness up to the current iteration k; update the global historical best position as follows. That is, the position with the minimum fitness among all the historical best positions of individual particles; update the velocity and position of each particle: ， (17) in, Let be the velocity of particle i at iteration k+1; Let be the position of particle i at iteration k+1; w be the inertia weight; c1 be the learning factor for updating particle velocity using the individual's historical best position; c2 be the learning factor for updating particle velocity using the global historical best position; r1 be the random number for updating particle velocity using the individual's historical best position; and r2 be the random number for updating particle velocity using the global historical best position.

10. The method for optimizing the dynamic transmission error distribution characteristics of a harmonic reducer according to claim 6, characterized in that: S504 includes the following steps: When the number of iterations k reaches the preset maximum value k max When the time is reached, the optimization process terminates; the optimal parameter compensation amount is output as follows: Therefore, the optimal subinterval is Ultimately, the optimal predicted dynamic transmission error distribution (DTE) is obtained. pred,opt And its statistical characteristics.

Images