Five-axis machining cutter location linear interpolation optimization control method

Abstract

The invention discloses a five-axis machining cutter location linear interpolation optimization control method, which relates to the technical field of numerical control machining, and comprises the following steps: fitting a machining contour to generate a curve expression, determining an operation point sequence based on curvature distribution, and converting the operation point sequence into an actuator pose sequence; based on the model of a five-axis machining machine tool, a total transformation matrix from a base coordinate system to an actuator coordinate system is established through a D-H method, motion parameter changes caused by five-axis elastic deformation and inertia force are established based on the nonlinear beam theory and the Alembert principle, and a nonlinear axis body coupling motion equation is established; a correction expression is generated in combination with the Jacobian matrix, optimal unknown parameters in the correction expression are determined by introducing a quantum fluctuation simulated annealing algorithm, and a correction Jacobian matrix is generated; according to the actuator pose change of the adjacent operation points in the actuator pose sequence, the corrected Jacobian matrix is combined, the inverse kinematics algorithm is adopted to derive the motion parameter vector of the previous operation point, and higher-precision machining control is achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical control machining, and in particular to a five-axis machining tool position linear interpolation optimization control method. Background Art

[0002] With the continuous advancement of surface contouring technology, the planning algorithm of five-axis machining trajectory is becoming increasingly mature. For example, the existing invention patent with announcement number CN119179302B proposes a five-axis CNC machining trajectory interpolation method, obtains the denture contour NURBS parameter curve, derives the homogeneous coordinate expression of the N-dimensional space NURBS parameter curve and obtains the coordinates of each point on the denture contour NURBS parameter curve, uses the numerical analysis algorithm to calculate the curve length, calculates the derivative vector on the denture contour NURBS parameter, determines the curvature change of the curve at the interpolation point, locates the interpolation point coordinates according to the curve length and derivative vector, and uses the interpolation point coordinates to plan the machining path, improve machining accuracy, and smooth the machining trajectory.

[0003] After obtaining the machining trajectory, the existing technology needs to solve the Jacobian matrix based on the robot kinematic modeling and reversely deduce the motion parameters of the five-axis machining center. However, the existing technology ignores the coupling error caused by the movement of the axes of the integrated five-axis machining center, resulting in the traditional Jacobian matrix being difficult to convert the machining trajectory into accurate motion parameters of the five-axis machining center. Summary of the Invention

[0004] In view of the shortcomings of the prior art, the present invention proposes a five-axis machining tool position linear interpolation optimization control method to achieve the motion parameter determination of the five-axis machining machine tool taking coupling errors into consideration.

[0005] The technical solutions for achieving the purpose of the present invention are:

[0006] A five-axis machining tool position linear interpolation optimization control method includes the following steps:

[0007] The NURBS curve is used to fit the machining contour to generate the curve expression p(u), and the operation point sequence is adaptively generated based on the curvature distribution κ(u). Transformed into the operating point coordinate sequence in the workpiece coordinate system through recursive algorithm The actuator pose sequence is converted back into the actuator coordinate system based on the machining process M is the total number of interpolation points;

[0008] Based on the five-axis machining center model, the DH method is used to establish the total transformation matrix T from the base coordinate system to the actuator coordinate system. Based on the nonlinear beam theory and the D'Alembert principle, the changes in the motion parameters caused by the elastic deformation and inertial force of the five axes are established respectively, and the nonlinear axis-body coupling motion equation is constructed. The Jacobian matrix J is combined to generate a correction expression and the optimal unknown parameters in the correction expression are determined by introducing the simulated annealing algorithm of quantum fluctuations to generate the corrected Jacobian matrix J. * ;

[0009] According to the actuator pose sequence The actuator pose changes of adjacent operating points in the equation are combined with the modified Jacobian matrix J * The inverse kinematics algorithm is used to reversely deduce the motion parameter vector q of the earlier operating point among the adjacent operating points, and converted into the operating power vector P of the servo motor through the PID control algorithm. r and execute.

[0010] Furthermore, the operating point sequence is generated The following steps are involved:

[0011] Using N control nodes and corresponding c-order B-spline basis functions in the workpiece coordinate system, the machining contour is fitted based on the NURBS curve to generate the curve expression p(u);

[0012] Calculate the curvature distribution κ(u) corresponding to the curve expression p(u) based on the curvature formula;

[0013] The machining contour length S is calculated based on the definite integral of the square root of the dot product of the first-order curve derivative vector p′(u) in the domain [0,1] of the parameter variable u;

[0014] Divide the machining contour into L segments of equal length and determine the left segmentation point u of the lth segment on the definition domain of the parameter variable u l and right segment point u l+1 ;

[0015] Calculate the mean curvature of segment l According to the average curvature of the first segment Determine the interpolation interval Δ of segment l l u;

[0016] According to the interpolation interval Δ of the first segment l u adds interpolation points in the lth segment, and generates an operating point sequence based on the sequential statistics of the segmentation points and interpolation points in the L segment

[0017] Furthermore, the actuator pose sequence is generated based on the processing technology conversion The following steps are involved:

[0018] The coordinates of the m-th operating point in the workpiece coordinate system are p(u m ) minus the origin offset vector between the workpiece coordinate system and the actuator coordinate system to generate the actuator coordinate p of the mth operating point e (u m );

[0019] According to the parameter variable u corresponding to the mth operating point m The first-order curve derivative vector p′(u m )Calculate the tangent vector T(u) of the mth operating point m );

[0020] According to the tangent vector T(u m )’s first-order tangent derivative vector T′(u m ) Calculate the principal normal vector N(u of the mth operating point m ), the tangent vector T(u m ) and the principal normal vector N(u m ) is used as the binormal vector B(u m );

[0021] Determine the first direction of the actuator at the mth operating point based on the processing technology is the tangent vector T(u m ) or principal normal vector N(u m ) and select a reference direction is the z-axis direction of the workpiece coordinate system;

[0022] The first direction of the actuator based on the mth operating point With reference direction The cross product result is used to calculate the second direction of the actuator at the mth operating point Set the actuator of the m-th operating point in the first direction With the actuator second direction The cross product result is used as the third direction of the actuator at the mth operating point

[0023] Calculate the first direction of the actuator separately The x0 axis direction of the base coordinate system and the second direction of the actuator The y0-axis direction of the base coordinate system and the third direction of the actuator The first angle with the z0 axis of the base coordinate system The second angle and the third angle Combine and construct the actuator pose of the mth operating point Generate actuator pose sequence

[0024] Furthermore, based on the DH method, the DH parameters of the i-th axis in the five-axis machining center are established, including the torsion angle θ i , axis length b i , offset d i and the rotation angle φ i , through the homogeneous transformation matrix When describing the motion of the i-th axis, the transformation relationship of the actuator coordinate system relative to the base coordinate system is simplified according to the motion state of the i-th axis, including:

[0025] When the i-th axis is a linear axis, the angle φ i and offset d i are 0° and the translation distance of the i-th axis respectively. According to the translation direction of the i-th axis, it is the x0 axis direction of the base coordinate system or the y0 axis direction and z0 axis direction of the base coordinate system. The torsion angle θ i 0° or 90° respectively;

[0026] When the i-th axis is a rotation axis, the offset d i and the rotation angle φ i The rotation angles of the 0th and i-th axes are respectively, based on the x0 axis and z0 axis of the base coordinate system or around the y0 axis of the base coordinate system, the torsion angle θ i 0° or -90° respectively;

[0027] The homogeneous transformation matrix of each axis Multiply them in sequence to generate the total transformation matrix T.

[0028] Furthermore, in traditional robot kinematics, the actuator pose ξ is obtained by the total transformation matrix T e Regarding the functional expression of the motion parameter vector q, the actuator posture ξ e The motion parameter q of the rth pose component about the i-th axis in the motion parameter vector q i The partial derivative can be used to calculate the element in the rth row and ith column of the Jacobian matrix J, where the motion parameter q of the i-th axis is i The motion state of the i-th axis can be determined as a translation distance or a rotation angle.

[0029] Furthermore, based on the nonlinear beam theory, it is assumed that the additional force of the motion of the j-th axis on the i-th axis is F i,j And cause the structural deformation of the i-th axis, i≠j, then the deformation offset of the i-th axis caused by the structural deformation The details are as follows:

[0030]

[0031] Among them, b i is the length of the i-th axis, β i is the deformation coefficient of the i-th axis, Ei , I i and A i are the material elastic modulus, section moment of inertia and cross-sectional area of ​​the i-th axis respectively.

[0032] Furthermore, based on the D'Alembert principle, the inertial force G of the i-th axis is i Equal to the mass m of the i-th axis i Multiply by the acceleration a of the i-th axis i , the i-th axis is affected by the inertial force G i The inertial offset Δq caused i g The details are as follows:

[0033] Δq i g =γ i [exp(η i G i )-1],

[0034] Among them, γ i and η i are the coefficient of inertia and the exponential coefficient of inertia respectively, and exp() is an exponential function with a natural constant as the base.

[0035] Furthermore, the actual motion parameters of the i-th axis Equal to the motion parameter q of the i-th axis i Add the inertial offset Δq of the i-th axis i g Add the deformation offset caused by the other four axes, that is, the nonlinear axis-body coupling motion equation is The actual motion parameters of the i-th axis Replace the homogeneous transformation matrix The corresponding motion parameter q i , get the actual actuator posture ξ e* The actual function expressions of motion parameter vector q, deformation coefficient vector β, inertia coefficient vector γ and exponential inertia coefficient vector η are calculated and the partial derivatives are used to construct the modified Jacobian matrix J with unknown variables. * , that is, the corrected expression.

[0036] Furthermore, the modified Jacobian matrix J is generated * , including the following steps:

[0037] Pre-collect multiple sets of starting actuator poses ξ e,o , motion parameter vector q and the final actuator pose ξ after execution e,e , respectively define the deformation search space, inertia search space and exponential inertia search space;

[0038] Define the Hamiltonian H(s) as the product of the initial Hamiltonian H0 and the complementary annealing parameter 1-s. p The product of the annealing parameter s, where the problem Hamiltonian H p The multiple sets of predicted end actuator poses are derived and the corresponding ending actuator pose ξ e,e The mean square error of

[0039] Initialize the maximum number of iterations K, the initial temperature T0, the initial fluctuation intensity Γ(0), and the attenuation coefficient ρ. Randomly generate a set of initial deformation coefficient vectors β(0), initial inertia coefficient vectors γ(0), and initial exponential inertia coefficient vectors η(0) in the deformation search space, inertia search space, and exponential inertia search space, and splice them to generate the initial quantum state x(0).

[0040] For the k-th iteration, set the k-th step annealing parameter s(k) equal to the number of iterations k divided by the maximum number of iterations K, update the k-th step fluctuation intensity Γ(k) equal to the initial fluctuation intensity Γ(0) multiplied by the k-th step complementary annealing parameter 1-s(k), multiply the Gaussian noise by the k-th step fluctuation intensity Γ(k) and superimpose it with the k-1-th step quantum state x(k-1) to obtain the k-th step quantum state x(k) and determine the k-th step modified Jacobian matrix J * (k);

[0041] Using the modified Jacobian matrix J of step k * (k) Derive multiple sets of predicted end actuator poses Calculate the Hamiltonian H(s(k)) of the kth step and subtract the Hamiltonian H(s(k-1)) of the k-1th step to get the difference Δ of the kth step k H;

[0042] If the difference in the k-th step Δ k If H is less than the difference threshold, the modified Jacobian matrix J of the kth step is * (k) is the modified Jacobian matrix J * and stop iterating;

[0043] If the difference in the k-th step Δ k If H is greater than or equal to the quantum difference threshold and less than 0, the quantum state x(k) of the kth step is received and the iteration of the k+1th step is started;

[0044] If the difference in the k-th step Δ k H is greater than or equal to 0, according to the probability of selection in step k Randomly select and receive the quantum state x(k) of the kth step or use the quantum state x(k-1) of the kth step as the quantum state x(k) of the kth step, and start the iteration of the k+1th step, where the selection probability of the kth step is By taking the difference Δ in the kth step k The product of the opposite number of H, the initial temperature T0 and the attenuation coefficient ρ to the kth power is substituted into the exponential function with the natural constant as the base;

[0045] The iteration stops when the maximum number of iterations K is reached and the modified Jacobian matrix J of the Kth step is converted to * (K) is the modified Jacobian matrix J * .

[0046] Furthermore, the inverse kinematics algorithm includes the following steps:

[0047] The actuator pose sequence The actuator pose of the m+1th operating point in The actuator pose with respect to the mth operating point As the end actuator pose ξ e,e and the initial starting actuator pose ξ e,o (0);

[0048] In the wth iteration, the ending actuator posture ξ obtained in the w-1th iteration is e,e (w-1) is the starting actuator pose ξ for the wth iteration e,o (w) and calculate and end the actuator pose ξ e,e The difference between the two values ​​is used to obtain the actuator posture change Δξ in the wth iteration. e,o (w);

[0049] Calculate the modified Jacobian matrix J * The pseudo-inverse pinv(J * ), using the damped least squares method, combined with the modified Jacobian matrix J * The pseudo-inverse pinv(J * ) and the actuator pose change Δξ in the wth iteration e,o (w), get the single-round motion parameter vector of the m-th operating point in the w-th round iteration

[0050] The motion parameter vector q of the mth operating point determined by the w-1th round of iteration is superimposed m (w-1) and the single-round motion parameter vector of the m-th operating point in the w-th round iteration Get the motion parameter vector q of the mth operating point determined by the wth round of iteration m (w), combined with the modified Jacobian matrix J * The predicted end actuator pose of the wth iteration is derived And calculate and end the actuator pose ξ e,e The mean square error of is determined to be less than or equal to the convergence threshold;

[0051] If it is greater than the convergence threshold, the w+1th round of iteration is started. If it is less than or equal to the convergence threshold, the motion parameter vector q of the mth operating point determined by the wth round of iteration is m (w) is the motion parameter vector q of the mth operating point m Output.

[0052] Compared with the existing technology, the present invention establishes the changes in motion parameters caused by the elastic deformation and inertial force of the five axes based on nonlinear beam theory and the D'Alembert principle, and constructs a nonlinear axis-body coupling motion equation. The Jacobian matrix is ​​combined to generate a corrected expression and the optimal unknown parameters in the corrected expression are determined by introducing a simulated annealing algorithm with quantum fluctuations to generate a corrected Jacobian matrix. According to the actuator posture changes of adjacent operating points in the actuator posture sequence, combined with the corrected Jacobian matrix, the inverse kinematics algorithm is used to derive the motion parameter vector of the earlier operating point, thereby weakening the influence of the coupling error under the premise of the same machining path planning accuracy and further improving the machining accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 This is a flow chart of a five-axis machining tool position linear interpolation optimization control method;

[0054] Figure 2 Generate an actuator pose sequence flow chart for the transformation;

[0055] Figure 3 Flowchart of the simulated annealing algorithm for introducing quantum fluctuations;

[0056] Figure 4 This is the flow chart of the inverse kinematics algorithm. DETAILED DESCRIPTION

[0057] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0058] Example 1

[0059] like Figure 1 As shown, a specific embodiment of the present invention discloses a five-axis machining tool position linear interpolation optimization control method, comprising the following steps:

[0060] Extract the machining contour of the workpiece 3D model, generate the curve expression p(u) by NURBS curve fitting, calculate the curvature distribution κ(u) based on geometric analysis, and adaptively generate the operating point sequence within the definition domain of the parameter variable u. Calculate the sequence of operating points based on a recursive algorithm The corresponding operating point coordinate sequence in the workpiece coordinate system Generate actuator pose sequence in actuator coordinate system based on machining process transformation M is the total number of interpolation points, and the recursive algorithm adopts the existing De Boer algorithm;

[0061] Based on the five-axis machining center model, the DH method is used to establish the total transformation matrix T from the base coordinate system to the actuator coordinate system. Based on the nonlinear beam theory, the motion parameter changes caused by the mechanical structure deformation between the five axes are established. Based on the D'Alembert principle, the motion parameter changes caused by the inertial force of a single axis are considered. The nonlinear axis-body coupling motion equation is constructed and combined with the traditional Jacobian matrix J to generate a correction expression. The quantum fluctuation auxiliary algorithm is introduced into the simulated annealing algorithm to effectively explore and determine the optimal unknown parameters in the correction expression in the high-dimensional solution space, generating the most practical correction Jacobian matrix J. * ;

[0062] According to the actuator pose sequence The actuator pose changes of adjacent operating points in the equation are combined with the modified Jacobian matrix J * The inverse kinematics algorithm is used to reversely deduce the motion parameter vector q of the five-axis machining center at the earlier operating point among the adjacent operating points, and is converted into the operating power vector P of the servo motor based on the existing PID control algorithm. r and execute, among others, running the power vector P i r is the operating power of the servo motor driving the i-th axis.

[0063] Furthermore, the operating point sequence is generated The following steps are involved:

[0064] The NURBS curve is used to fit the machining contour and generate the curve expression p(u), as follows:

[0065]

[0066] Where N is the total number of control nodes, ω n 、 and f n,c (u) are the control weights of n control nodes, the coordinates of the control nodes in the workpiece coordinate system and the c-order B-spline basis function, where c is the order;

[0067] The curvature distribution κ(u) corresponding to the curve expression p(u) is calculated based on the curvature formula, as follows:

[0068]

[0069] Where p′(u) and p″(u) are the first-order curve derivative vector and the second-order curve derivative vector of the curve expression p(u), respectively;

[0070] The machining contour length is calculated based on the definite integral of the square root of the dot product of the first-order curve derivative vector p′(u) in the domain [0,1] of the parameter variable u. Divide the machining contour into L segments of equal length and determine the left segmentation point u of the lth segment on the definition domain of the parameter variable u l and right segment point u l+1 ;

[0071] Calculate the mean curvature of segment l According to the average curvature of the first segment Determine the interpolation interval Δ of segment l l u, specifically as follows:

[0072]

[0073] Among them, Δ min u, u0 and ε are the minimum interpolation interval, the preset interpolation interval and the minimum value to prevent the denominator from being 0, respectively;

[0074] According to the interpolation interval Δ of the first segment l u adds interpolation points in the lth segment, and generates an operating point sequence based on the sequential statistics of the segmentation points and interpolation points in the L segment

[0075] like Figure 2 As shown, further, the actuator pose sequence is generated based on the processing technology conversion The following steps are involved:

[0076] From the operating point coordinate sequence Get the coordinates of the mth operating point in the workpiece coordinate system p(u m ) and subtract the origin offset vector between the workpiece coordinate system and the actuator coordinate system to generate the actuator coordinate p of the mth operating point e (u m );

[0077] According to the parameter variable u corresponding to the mth operating point m The first-order curve derivative vector p′(u m )Calculate the tangent vector T(u) of the mth operating point m )=p′(u m ) / ||p′(u m )||;

[0078] According to the tangent vector T(u m )’s first-order tangent derivative vector T′(u m ) Calculate the principal normal vector N(u of the mth operating point m )=T′(u m ) / ||T′(u m )||, the tangent vector T(u m ) and the principal normal vector N(u m) is used as the binormal vector B(u m )=T(u m )×N(u m ), the tangent vector T(u m ), principal normal vector N(u m ) and the binormal vector B(u m ) constructs the Freynes price and constructs the local orthogonal coordinate system for the mth operating point;

[0079] Determine the first direction of the actuator at the mth operating point according to whether the machining process is tangential cutting or vertical cutting is the tangent vector T(u m ) or principal normal vector N(u m ), select the reference direction is the z-axis direction of the workpiece coordinate system;

[0080] The first direction of the actuator based on the mth operating point With reference direction The cross product result is used to calculate the second direction of the actuator at the mth operating point Set the actuator of the m-th operating point in the first direction With the actuator second direction The cross product result is used as the third direction of the actuator at the mth operating point

[0081] Calculate the first direction of the actuator separately The x0 axis direction of the base coordinate system and the second direction of the actuator The y0-axis direction of the base coordinate system and the third direction of the actuator The first angle with the z0 axis of the base coordinate system The second angle and the third angle Combine and construct the actuator pose of the mth operating point Generate actuator pose sequence

[0082] Furthermore, the five-axis machining center includes three linear axes and two rotation axes, wherein the three linear axes translate along the x0-axis, y0-axis and z0-axis directions of the base coordinate system respectively, and the two rotation axes rotate around two of the x0-axis, y0-axis and z0-axis of the base coordinate system. The DH parameters of the i-th axis are established based on the DH method, including the torsion angle θ i , axis length b i , offset d i and the rotation angle φ i , through the homogeneous transformation matrix Describes the transformation relationship between the actuator coordinate system and the base coordinate system when the i-th axis moves, the homogeneous transformation matrix The details are as follows:

[0083]

[0084] Among them, s() and c() represent the sine function and cosine function respectively. According to the motion state of the i-th axis, the homogeneous transformation matrix of the i-th axis can be further simplified. include:

[0085] When the i-th axis is a linear axis translated along the x0 axis of the base coordinate system, the torsion angle θ i and the rotation angle φ i All are 0°, offset d i is the translation distance of the i-th axis;

[0086] When the i-th axis is a linear axis translated along the y0-axis or z0-axis direction of the base coordinate system, the torsion angle θ i and the rotation angle φ i are 90° and 0° respectively, with an offset of d i is the translation distance of the i-th axis;

[0087] When the i-th axis is the rotation axis around the x0 axis or z0 axis of the base coordinate system, the torsion angle θ i and offset d i are 0° and 0 respectively, and the rotation angle φ i is the rotation angle of the i-th axis;

[0088] When the i-th axis is the rotation axis around the y0 axis of the base coordinate system, the torsion angle θ i and offset d i -90° and 0 respectively, the rotation angle φ i is the rotation angle of the i-th axis;

[0089] The homogeneous transformation matrix of each axis Multiply them in sequence to generate the total transformation matrix

[0090] Furthermore, under the premise of knowing the total transformation matrix T, the traditional robot kinematics uses the Jacobian matrix J to describe the actuator posture ξ caused by the change of the motion parameter vector q e Change the actuator posture of the m+1th operating point Subtract the actuator pose of the mth operating point And divided by the operation time to get the actuator velocity vector of the mth operation point The motion parameter vector q of the mth operating point m =[q m,1 ,q m,2 ,qm,3 ,q m,4 ,q m,5 ] divided by the operation time to obtain the motion velocity vector v of the mth operation point m , then the motion parameter vector q of the mth operating point m Equal to the Jacobian matrix J multiplied by the motion velocity vector v of the mth operating point m Specifically, the actuator pose ξ is obtained through the total transformation matrix T e Regarding the functional expression of the motion parameter vector q, the actuator posture ξ e The motion parameter q of the rth pose component about the i-th axis in the motion parameter vector q i The partial derivative can be used to calculate the element in the rth row and ith column of the Jacobian matrix J, where q m,i is the motion parameter of the i-th axis at the m-th operating point. According to the motion state of the i-th axis, it can be determined as the corresponding translation distance or rotation angle, k∈{1,2,3,4,5,6}.

[0091] Furthermore, during the operation of the five-axis machining center, the deformation of the mechanical structure is an important factor that causes nonlinear coupling between the axes. Based on the nonlinear beam theory, it is assumed that the additional force of the motion of the j-th axis on the i-th axis is F i,j And cause the structural deformation of the i-th axis, i≠j, then the deformation offset of the i-th axis caused by the structural deformation The details are as follows:

[0092]

[0093] Among them, b i is the length of the i-th axis, E i , I i and A i are the material elastic modulus, section inertia moment and cross-sectional area of ​​the i-th axis respectively, the material elastic modulus E i Reflects the ability of the i-th axis to resist elastic deformation, which is only related to the manufacturing material of the i-th axis. The moment of inertia of the section I i Reflects the ability of the i-th axis to resist bending deformation, which is only related to the cross-sectional shape and size of the i-th axis. i is the deformation coefficient of the i-th axis, the first term on the right (b i ) 3 (3E i I i ) -1 F i,j Essentially an additional force F i,j The linear change caused by the second term β on the right side i b i (E i A i )-2 (F i,j ) 2 Essentially an additional force F i,j caused by nonlinear changes.

[0094] Furthermore, based on the D'Alembert principle, the inertia force G of the i-th axis during the operation of the five-axis machining center is considered. i Equal to the mass m of the i-th axis i Multiply by the acceleration a of the i-th axis i , the i-th axis is affected by the inertial force G i The inertial offset Δq caused i g is as follows:

[0095] Δq i g =γ i [exp(η i G i )-1],

[0096] Among them, γ i and η i They are the inertia coefficient and the exponential inertia coefficient respectively. exp() is an exponential function with a natural constant as the base. It is constructed in the form of an exponential to reflect the inertia force G. i When it is small, the inertial offset Δq i g Smooth, inertial force G i When it is large, the inertial offset Δq i g A more drastic nonlinear change trend.

[0097] Furthermore, the Jacobian matrix J determined in traditional robot kinematics essentially treats the five axes of a five-axis machining center as independent entities. However, in actual machining, the additional forces between the five axes and the inertia caused by the motion of the five axes themselves will produce coupling errors. The nonlinear axis coupling motion equation is constructed. Given the motion parameter q of the i-th axis, i Under the premise of Equal to the motion parameter q of the i-th axis i Add the inertial offset Δq of the i-th axis i g Add the deformation offset caused by the other four axes, that is, the nonlinear axis-body coupling motion equation is The actual motion parameters of the i-th axis Replace the homogeneous transformation matrix The corresponding motion parameter q i , get the actual actuator posture ξ e*The actual function expressions of the motion parameter vector q, the deformation coefficient vector β = [β1, β2, β3, β4, β5], the inertia coefficient vector γ = [γ1, γ2, γ3, γ4, γ5], and the exponential inertia coefficient vector η = [η1, η2, η3, η4, η5] are derived according to the same way to construct the correction expression, which is essentially a correction Jacobian matrix J * , which reflects the coupling error caused by considering the deformation of the structure and the inertia force, and is more in line with the actual change of the actual executor pose ξ e* caused by the change of the motion parameter vector q in the actual processing process.

[0098] As shown in Figure 3 , further, the correction Jacobian matrix J * is generated, including the following steps:

[0099] A plurality of groups of initial executor poses ξ e,o , motion parameter vectors q, and final executor poses ξ e,e executed by the five-axis machining tool are collected in advance, and the deformation search space [β dw , β up ], the inertia search space [γ dw , γ up ], and the exponential inertia search space [η dw , η up ] of the deformation coefficient β, the inertia coefficient γ, and the exponential inertia coefficient η are defined, respectively, wherein β dw , γ dw , and η dw are the lower limit deformation coefficient, the lower limit inertia coefficient, and the lower limit exponential inertia coefficient, respectively, and β up , γ up , and η up are the upper limit deformation coefficient, the upper limit inertia coefficient, and the upper limit exponential inertia coefficient, respectively.

[0100] The Hamiltonian H(s) is defined as the product of the initial Hamiltonian H0 and the complementary annealing parameter 1-s plus the product of the problem Hamiltonian H p and the annealing parameter s, wherein the problem Hamiltonian H p is the mean square error of the plurality of groups of predicted final executor poses and the corresponding final executor poses ξ e,e ;

[0101] The maximum iteration step number K, the initial temperature T0, the initial fluctuation strength Γ(0), and the decay coefficient ρ are initialized, and the deformation search space [β dw , β up ], the inertia search space [γ dw , γ up ], and the exponential inertia search space [ηdw ,η up ] randomly generates a set of initial deformation coefficient vector β(0), initial inertia coefficient vector γ(0) and initial exponential inertia coefficient vector η(0), and splices them together to generate the initial quantum state x(0) = [β(0), γ(0), η(0)];

[0102] For the k-th iteration, set the k-th step annealing parameter s(k) = k / K, update the k-th step fluctuation intensity Γ(k) = Γ(0)(1-s(k)), perform quantum fluctuations to avoid falling into the optimal, and replace the Gaussian noise N(0,σ 2 ) is multiplied by the fluctuation intensity Γ(k) of the kth step and superimposed with the quantum state x(k-1) of the k-1th step, and the quantum state x(k) of the kth step is obtained. 2 ) and determine the modified Jacobian matrix J of the kth step * (k);

[0103] Using the modified Jacobian matrix J of step k * (k) Combining multiple sets of starting actuator poses ξ e,o and motion parameter vector q to derive multiple sets of predicted end actuator poses Calculate the Hamiltonian H(s(k)) of the kth step and subtract the Hamiltonian H(s(k-1)) of the k-1th step to get the difference Δ of the kth step k H;

[0104] If the difference in the k-th step Δ k If H is less than the difference threshold, the modified Jacobian matrix J of the kth step is * (k) is directly used as the modified Jacobian matrix J * and stop iterating;

[0105] If the difference in the k-th step Δ k If H is greater than or equal to the quantum difference threshold and less than 0, the quantum state x(k) of the kth step is received and the iteration of the k+1th step is started;

[0106] If the difference in the k-th step Δ k H is greater than or equal to 0, according to the probability of selection in step k Choose to receive the quantum state x(k) at step k or subtract the probability of selection at step k from 1 The quantum state x(k-1) of the k-1th step is directly used as the quantum state x(k) of the kth step, and the iteration of the k+1th step is started, where the selection probability of the kth step is By taking the difference Δ in the kth step k The product of the opposite number of H, the initial temperature T0 and the attenuation coefficient ρ to the kth power is substituted into the exponential function with the natural constant as the base;

[0107] The iteration stops when the maximum number of iterations K is reached and the modified Jacobian matrix J of the Kth step is converted to * (K) is the modified Jacobian matrix J * .

[0108] like Figure 4 As shown, further, the inverse kinematics algorithm includes the following steps:

[0109] The actuator pose sequence The actuator pose of the m+1th operating point in The actuator pose with the mth operating point As the end actuator pose ξ e,e and the initial starting actuator pose ξ e,o (0);

[0110] In the wth iteration, the ending actuator posture ξ obtained in the w-1th iteration is e,e (w-1) is the starting actuator pose ξ for the wth iteration e,o (w) and calculate and end the actuator pose ξ e,e The difference between the two values ​​is used to obtain the actuator posture change Δξ in the wth iteration. e,o (w)=Δξ e,o (0)-ξ e,o (w), since the operation time between adjacent operation points is fixed, the change in actuator posture is essentially equivalent to the change in the actuator velocity vector;

[0111] To avoid the modified Jacobian matrix J * If it is a non-square matrix or a singular situation occurs, the pseudo-inverse operation is used instead of the inverse operation to calculate the modified Jacobian matrix J * The pseudo-inverse pinv(J * )=[(J * ) T J * ] -1 (J * ) T ,in,() T and[] -1 Represent the transpose and inverse operations of the matrix respectively;

[0112] The damped least squares method is used, combined with the modified Jacobian matrix J * The pseudo-inverse pinv(J * ) and the actuator pose change Δξ in the wth iteration e,o (w), get the single-round motion parameter vector of the m-th operating point in the w-th round iteration Where ψ is the damping coefficient, which is used to avoid the divergence of singular points;

[0113] The motion parameter vector q of the m-th operating point determined in the first w-1 rounds of iterations m The single-round motion parameter vector of the m-th operating point in the w-th round of iteration is superimposed on (w-1) Get the motion parameter vector of the mth operating point determined by the wth round of iteration Combined with the modified Jacobian matrix J * and the initial starting actuator pose ξ e,o (0) Perform forward kinematic derivation to determine the predicted actuator pose at the end of the wth iteration. and the end actuator pose ξ e,e Is the mean square error of less than or equal to the convergence threshold?

[0114] If it is greater than the convergence threshold, the w+1th round of iteration is started. If it is less than or equal to the convergence threshold, the motion parameter vector q of the mth operating point determined by the wth round of iteration is m (w) is the motion parameter vector q of the mth operating point m Output.

[0115] The present invention discloses a five-axis machining tool position linear interpolation optimization control method, which adopts NURBS curve fitting to generate a curve expression for machining contour, determines an operation point sequence based on curvature distribution and converts it into an actuator posture sequence; based on the five-axis machining machine model, adopts the DH method to establish a total transformation matrix from a base coordinate system to an actuator coordinate system, establishes the motion parameter changes caused by the elastic deformation and inertial force of the five axes based on nonlinear beam theory and the D'Alembert principle, and constructs a nonlinear axis-body coupling motion equation; combines the Jacobian matrix to generate a correction expression, and determines the optimal unknown parameters in the correction expression by introducing a simulated annealing algorithm with quantum fluctuations to generate a corrected Jacobian matrix; according to the actuator posture changes of adjacent operation points in the actuator posture sequence, combined with the corrected Jacobian matrix, adopts an inverse kinematics algorithm to derive the motion parameter vector of the earlier operation point, thereby achieving higher-precision machining control.

[0116] The above description is merely a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiment. All technical solutions based on the concept of the present invention are within the scope of protection of the present invention. It should be noted that for those skilled in the art, various improvements and modifications that do not depart from the principles of the present invention should also be considered within the scope of protection of the present invention.

Claims

1. A five-axis machining tool position linear interpolation optimization control method, characterized in that: The following steps are involved: The NURBS curve is used to fit the machining contour to generate a curve expression, and the operating point sequence is adaptively generated based on the curvature distribution. The sequence is converted into the actuator pose sequence in the actuator coordinate system based on the recursive algorithm and machining process. Based on the five-axis machining center model, the DH method is used to establish the total transformation matrix from the base coordinate system to the actuator coordinate system. Based on nonlinear beam theory and the d'Alembert principle, the changes in the motion parameters caused by the elastic deformation and inertial force of the five axes are respectively established, and the nonlinear axis-body coupling motion equation is constructed. The Jacobian matrix is ​​combined to generate a correction expression. The optimal unknown parameters in the correction expression are determined by introducing a simulated annealing algorithm with quantum fluctuations to generate a corrected Jacobian matrix. According to the actuator posture changes of adjacent operating points in the actuator posture sequence, combined with the modified Jacobian matrix, the inverse kinematics algorithm is used to reversely deduce the motion parameter vector of the earlier operating point among the adjacent operating points.

2. A five-axis machining tool position linear interpolation optimization control method according to claim 1, characterized in that: Based on the nonlinear beam theory, it is assumed that the additional force acting on the ith axis due to the motion of the jth axis is F i,j And cause the structural deformation of the i-th axis, i≠j, then the deformation offset of the i-th axis caused by the structural deformation The details are as follows: Among them, b i is the length of the i-th axis, β i is the deformation coefficient of the i-th axis, E i , I i and A i are the material elastic modulus, section moment of inertia and cross-sectional area of ​​the i-th axis respectively.

3. A five-axis machining tool position linear interpolation optimization control method as claimed in claim 2, characterized in that: Based on the D'Alembert principle, the inertial force G of the i-th axis is i Equal to the mass m of the i-th axis i Multiply by the acceleration a of the i-th axis i , the i-th axis is affected by the inertial force G i The inertial offset Δq caused i g The details are as follows: Δq i g =c i [exp(η i G i )-1], Among them, γ i and η i are the coefficient of inertia and the exponential coefficient of inertia respectively, and exp() is an exponential function with a natural constant as the base.

4. A five-axis machining tool position linear interpolation optimization control method as claimed in claim 3, characterized in that: The nonlinear axis-body coupling motion equation is that the actual motion parameter of the i-th axis is equal to the motion parameter of the i-th axis plus the inertia offset of the i-th axis plus the deformation offset caused by the other four axes. The actual motion parameter of the i-th axis is replaced with the corresponding motion parameter in the homogeneous transformation matrix to obtain the actual function expression of the actual actuator posture with respect to the motion parameter vector, deformation coefficient vector, inertia coefficient vector and exponential inertia coefficient vector, and the partial derivative is calculated to construct the corrected expression.

5. The five-axis machining tool position linear interpolation optimization control method according to claim 4, characterized in that: Generating the modified Jacobian matrix involves the following steps: At the k-th iteration, the annealing parameter of the k-th step is set equal to the number of iteration steps divided by the maximum number of iteration steps, the fluctuation intensity of the k-th step is updated equal to the initial fluctuation intensity multiplied by the complementary annealing parameter of the k-th step, Gaussian noise is multiplied by the fluctuation intensity of the k-th step and superimposed with the quantum state of the k-th step to obtain the quantum state of the k-th step and determine the modified Jacobian matrix of the k-th step, where the quantum state includes a set of deformation coefficient vectors, inertia coefficient vectors, and exponential inertia coefficient vectors; Utilizing the modified Jacobian matrix of the kth step, multiple sets of predicted ending actuator poses are derived. The product of the initial Hamiltonian and the complementary annealing parameter of the kth step plus the product of the problem Hamiltonian of the kth step and the annealing parameter of the kth step is taken as the Hamiltonian of the kth step, and the Hamiltonian of the k-1th step is subtracted to obtain the quantity difference of the kth step, where the problem Hamiltonian of the kth step is the mean square error between the multiple sets of predicted ending actuator poses derived based on the modified Jacobian matrix of the kth step and the corresponding ending actuator poses. If the quantity difference of the k-th step is less than the quantity difference threshold, the modified Jacobian matrix of the k-th step is used as the modified Jacobian matrix and the iteration is stopped. In other cases, the iteration of the k+1-th step is started and further, based on whether the quantity difference of the k-th step is less than 0, a decision is made to directly receive the quantum state of the k-th step or to randomly select the quantum state x(k) of the k-th step based on the selection probability of the k-th step. The selection probability of the k-th step is obtained by substituting the product of the opposite of the quantity difference of the k-th step and the initial temperature and the attenuation coefficient to the kth power into an exponential function with a natural constant as the base; The iteration is stopped when the maximum number of iteration steps is reached and the modified Jacobian matrix of the Kth step is used as the modified Jacobian matrix.

6. A five-axis machining tool position linear interpolation optimization control method according to any one of claims 1 to 5, characterized in that: Generate an operating point sequence, including the following steps: Using N control nodes and corresponding c-order B-spline basis functions in the workpiece coordinate system, the machining contour is fitted based on the NURBS curve to generate a curve expression; The curvature distribution corresponding to the curve expression is calculated based on the curvature formula, and the machining contour length is calculated based on the definite integral of the square root of the dot product of the first-order curve derivative vector within the definition domain of the parameter variable; Divide the machining contour into L segments of equal length and determine the left segmentation point and the right segmentation point of the lth segment on the definition domain of the parameter variable, and determine the interpolation interval of the lth segment based on the average curvature of the lth segment; According to the interpolation interval of the lth segment, an interpolation point is added in the lth segment, and the segmentation points and interpolation points in the Lth segment are sequentially counted to generate an operation point sequence.

7. A five-axis machining tool position linear interpolation optimization control method as claimed in claim 6, characterized in that: The conversion into the actuator pose sequence in the actuator coordinate system based on the recursive algorithm and machining process includes the following steps: Calculate the operating point coordinate sequence corresponding to the operating point sequence in the workpiece coordinate system based on a recursive algorithm; Subtract the origin offset vector between the workpiece coordinate system and the actuator coordinate system from the operating point coordinate of the m-th operating point to generate the actuator coordinate of the m-th operating point; Calculate the tangent vector of the mth operating point based on the first-order curve derivative vector at the parameter variable corresponding to the mth operating point and further derive the principal normal vector of the mth operating point. Take the cross product of the tangent vector of the mth operating point and the principal normal vector as the binormal vector of the mth operating point. Determine the first direction of the actuator of the mth operating point as a tangent vector or a principal normal vector according to the machining process and select the reference direction as the z-axis direction of the workpiece coordinate system; Calculating the second direction of the actuator at the mth operating point based on a cross product of the first direction of the actuator at the mth operating point and the reference direction, and using the cross product of the first direction of the actuator at the mth operating point and the second direction of the actuator as the third direction of the actuator at the mth operating point; The first angle, second angle, and third angle between the actuator's first direction and the x0-axis direction of the base coordinate system, the actuator's second direction and the y0-axis direction of the base coordinate system, and the actuator's third direction and the z0-axis direction of the base coordinate system are calculated respectively. The actuator pose of the m-th operating point is constructed by combining them to generate an actuator pose sequence.

8. The five-axis machining tool position linear interpolation optimization control method according to claim 7, characterized in that: Based on the DH method, the DH parameters of the i-th axis in the five-axis machining center are established, including the torsion angle, axis length, offset and rotation angle. When the motion of the i-th axis is described by the homogeneous transformation matrix, the transformation relationship between the actuator coordinate system and the base coordinate system is simplified according to the motion state of the i-th axis, including: When the i-th axis is a linear axis, the rotation angle and offset are 0° and the translation distance of the i-th axis respectively. Depending on whether the translation direction of the i-th axis is the x0 axis direction of the base coordinate system or the y0 axis direction and the z0 axis direction of the base coordinate system, the torsion angle is 0° or 90° respectively. When the i-th axis is the rotation axis, the offset and angle are 0 and the rotation angle of the i-th axis respectively. According to the axis of rotation as the x0 axis and z0 axis of the base coordinate system or around the y0 axis of the base coordinate system, the torsion angle is 0° or -90° respectively; Multiply the homogeneous transformation matrices for each axis in sequence to generate the total transformation matrix.

9. A five-axis machining tool position linear interpolation optimization control method as claimed in claim 8, characterized in that: In traditional robot kinematics, the function expression of the actuator posture with respect to the motion parameter vector is obtained through the total transformation matrix. The element in the rth row and i-th column of the Jacobian matrix can be calculated by taking the partial derivative of the rth posture component in the actuator posture with respect to the motion parameter of the i-th axis in the motion parameter vector. Among them, the motion parameter of the i-th axis can be determined as the translation distance or rotation angle based on the motion state of the i-th axis.

10. The five-axis machining tool position linear interpolation optimization control method according to claim 9, characterized in that: The inverse kinematics algorithm consists of the following steps: The actuator pose of the m+1th operating point and the actuator pose of the mth operating point in the actuator pose sequence are respectively used as the ending actuator pose and the initial starting actuator pose; In the w-th iteration, the ending actuator pose obtained in the w-1th iteration is used as the starting actuator pose of the w-th iteration and the difference between it and the ending actuator pose is calculated to obtain the actuator pose change of the w-th iteration. Calculate the pseudo-inverse of the modified Jacobian matrix and use the damped least squares method to combine the pseudo-inverse of the modified Jacobian matrix with the actuator posture change in the w-th iteration to obtain the single-round motion parameter vector of the m-th operating point in the w-th iteration; The motion parameter vector of the mth operating point determined by the w-1th iteration is superimposed with the single-round motion parameter vector of the mth operating point in the wth iteration to obtain the motion parameter vector of the mth operating point determined by the wth iteration. The predicted ending actuator pose of the wth iteration is derived by combining the modified Jacobian matrix and the ending actuator pose ξ is calculated. e,e The mean square error of is determined to be less than or equal to the convergence threshold; If it is greater than the convergence threshold, the w+1th round of iteration is started. If it is less than or equal to the convergence threshold, the motion parameter vector of the mth operation point determined by the wth round of iteration is output as the motion parameter vector of the mth operation point.

Images