A method for milling pose optimization based on redundant axes milling

Abstract

The application belongs to the technical field of robot pose optimization, and particularly relates to a milling pose optimization method based on a milling redundant axis, wherein the milling spindle is vertically installed on the sixth axis of a robot, the freedom redundancy of the milling cutter pose is realized by using six degrees of freedom, and the rotation of the z-axis of the milling cutter coordinate system is described by using the redundant freedom; the normal vector of the surface point position is solved by using the point cloud model based on the PCA method, and the direction vector of the milling redundant axis is determined; then the pose optimization target of the six-degree-of-freedom joint robot under the typical five-degree-of-freedom milling operation is set, and the angle between the external load and the main stiffness direction, the stiffness performance index, the motion performance index and the motion stability are covered; finally, the angle optimization is performed in the range of -180°~180° by using the MO-MLSWOAR multi-objective optimization algorithm, and the best milling posture of the corresponding milling point position is obtained. The method can fully utilize the redundant freedom, accurately match the surface machining characteristics, realize the multi-objective collaborative optimization, and improve the milling quality and efficiency.

Description

Technical Field

[0001] This invention belongs to the field of robot pose optimization technology, and particularly relates to a milling pose optimization method based on milling redundant axes. Background Technology

[0002] In modern manufacturing, milling is an important material removal process widely used in the forming and precision machining of complex parts. With the development of industrial robot technology, six-degree-of-freedom articulated robots, due to their high flexibility and operating range, are gradually being applied to the milling field to meet the machining needs of workpieces with complex curved surfaces and irregular shapes.

[0003] However, traditional milling pose control methods still have many limitations: On the one hand, some milling equipment or robot systems use a five-degree-of-freedom configuration, which lacks redundant degrees of freedom, resulting in insufficient flexibility in adjusting the milling cutter pose. It is difficult to achieve optimal pose adaptation in complex machining scenarios, especially when dealing with complex workpieces such as free-form surfaces, which can easily lead to problems such as low machining accuracy and rapid tool wear. On the other hand, even when using a six-degree-of-freedom robot, the existing technology does not make full use of redundant degrees of freedom and fails to effectively achieve accurate description and control of the milling cutter coordinate system rotation (i.e., milling spindle deflection) through redundant degrees of freedom, making it difficult to match the dynamic requirements of different machining points.

[0004] In solving for redundant axis direction vectors, existing methods mostly rely on empirical values ​​or simplified models, lacking accurate extraction of local geometric features of the machined surface. For example, in freeform surface milling, if the normal vectors of each point on the surface cannot be accurately obtained, the redundant axis directions will not match the surface characteristics, thus affecting the effectiveness of subsequent pose optimization.

[0005] The setting of pose optimization objectives also suffers from a simplistic approach. Traditional methods often focus only on a single performance metric (such as stiffness or dexterity), neglecting the need for multi-objective collaborative optimization. In actual milling processes, the matching degree between the external loads (including milling forces and self-load) on the robot end effector and the principal stiffness directions, the volume and anisotropy of the end effector stiffness ellipsoid, the robot's dexterity, and the stability between adjacent poses all affect machining quality and equipment lifespan. Optimizing a single objective makes it difficult to consider overall performance.

[0006] In addition, existing multi-objective optimization algorithms have shortcomings in milling pose optimization: some algorithms have unreasonable leader mechanism design, which makes it difficult to balance the diversity and convergence of Pareto front solution sets; the speed update strategy lacks adaptive adjustment, which easily leads to slow convergence speed or getting trapped in local optima; at the same time, traditional algorithms have insufficient compatibility when integrating different search strategies, which affects optimization efficiency and accuracy, and makes it difficult to quickly find the best milling pose in a wide range of angles from -180° to 180°.

[0007] Therefore, in view of the shortcomings of existing technologies in robot milling pose optimization in terms of degree of freedom utilization, redundant axis solution, multi-objective cooperation and optimization algorithm performance, there is an urgent need for a milling pose optimization method that can make full use of redundant degrees of freedom, accurately extract surface features, integrate multi-objective performance and efficiently find the best, so as to improve the milling quality and processing efficiency of complex workpieces. Summary of the Invention

[0008] The purpose of this invention is to provide a milling pose optimization method based on redundant axes, in order to solve the technical problems of existing robot milling pose optimization in terms of degree of freedom utilization, redundant axis solving, multi-objective cooperation and optimization algorithm performance.

[0009] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0010] A milling pose optimization method based on redundant milling axes includes the following steps:

[0011] S1: The milling spindle is vertically mounted on the sixth axis of the robot, so that the robot uses six degrees of freedom to determine the milling pose, realizing the redundancy of the kinematics of the milling cutter pose, and using the redundant degrees of freedom to describe the rotation of the milling cutter coordinate system z-axis, that is, the deflection of the milling spindle.

[0012] S2: Solving the milling redundant axis direction vector based on PCA method;

[0013] S3: Set the pose optimization objective for a six-DOF articulated robot in a typical five-DOF milling operation;

[0014] S4: The MO-MLSWOAR multi-objective optimization algorithm is used to optimize the angle within the range of -180° to 180° to obtain the optimal milling posture for the corresponding milling point.

[0015] Preferably, in step S1, the rotation of the milling cutter coordinate system z-axis is described with redundant degrees of freedom, that is, the specific process of the milling spindle deflection is as follows:

[0016] S11: In a milling task, the pose of the milling cutter coordinate system {M} in the base coordinate system {B} is described as a homogeneous transformation matrix:

[0017]

[0018] in, Each column in the table represents the projection of the coordinate axes of the milling cutter coordinate system {M} onto the base coordinate system {B}, and the rotation axis of the milling cutter. The unit vector in the base coordinate system {B} is:

[0019]

[0020] S12: Redundancy of degrees of freedom based on the milling cutter pose, the milling cutter coordinate system {M} revolves around the milling cutter's rotation axis. For rotation, if the rotation angle is γ, according to Rodrigues' rotation formula, the rotation matrix is... Represented as:

[0021]

[0022] Where I is the identity matrix, and K is the antisymmetric matrix of the rotation axis:

[0023]

[0024] S13: New pose of the milling cutter coordinate system {M} after rotation for:

[0025]

[0026] S14: Based on the relative poses between coordinate systems and the robot's kinematic equations, establish the relationship between the robot's joint angles and the rotation angle of the milling cutter coordinate system {M}:

[0027]

[0028] S15: When the robot performs a typical milling task, in the world coordinate system {W}, the pose of the robot cutting system is described as the milling position, the milling spindle direction, and the deflection of the milling spindle:

[0029]

[0030] in, For the milling position, a m γ represents the direction of the milling spindle, and γ is the angle of rotation of the milling coordinate system around the milling spindle.

[0031] Preferably, in step S2, the solution of the milling redundant axis direction vector based on the PCA method is achieved by solving the normal vector of a point p(x, y, z) on a point cloud model obtained from a freeform surface. The specific process is as follows:

[0032] S21: Selecting a Field:

[0033] The neighborhood points p1, p2, ..., p of p are selected using Euclidean distance. n The neighboring points satisfy:

[0034]

[0035] Where r is the radius of the sphere's neighborhood;

[0036] Calculate the mean vector of the neighborhood points:

[0037]

[0038] Use the set of neighboring points as sample D:

[0039]

[0040] S22: Construct the covariance matrix:

[0041] Center sample D:

[0042]

[0043] Its covariance matrix ∑ is:

[0044]

[0045] in, It is the variance of x; Let x and y be the covariance values, and so on, to solve for the covariance matrix ∑.

[0046] S23: Solving for eigenvalues ​​and eigenvalue vectors:

[0047] The eigenvalues ​​and eigenvectors of the covariance matrix Σ satisfy:

[0048] Σv i =λ i v i ;

[0049] The eigenvalues ​​are solutions to the characteristic equation of the covariance matrix ∑.

[0050] det(∑-λ i I) = 0;

[0051] S24: For a 3×3 covariance matrix, there are three eigenvalues ​​λ1, λ2, and λ3. The magnitude of the eigenvalue represents the variance of the data in the direction of the corresponding eigenvector. For normal vector calculation, the eigenvector corresponding to the smallest eigenvalue reflects the direction of least change of the point cloud data in that local region. Therefore, the eigenvector with the smallest eigenvalue is selected as the normal vector.

[0052] λ n =min(λ1,λ2,λ3);

[0053] The eigenvector v is obtained by solving the following equation. n :

[0054] (∑-λ n I)v n =0,

[0055] Feature vector v n That is, the normal vector at point p.

[0056] Preferably, the specific process for setting the pose optimization target of the six-DOF articulated robot under a typical five-DOF milling operation in step S3 is as follows:

[0057] S31: The angle between the principal stiffness direction and the resultant force of the external load.

[0058] The external loads at the robot end effector include the milling force F. m and the original load F o This includes flanges, electric spindles, and force sensors, where the milling force F... m It is the original load F acting in the world coordinate system. o It acts on the end face, applying the external load F. ext Represented as:

[0059]

[0060] External load F ext With respect to the principal stiffness direction The included angle θ Fv1 for:

[0061]

[0062] S32: Characterizing stiffness performance indicators: end-stiffness ellipsoid volume and stiffness anisotropic elliptic area; setting the expression for end-stiffness ellipsoid volume:

[0063]

[0064] S33: Characterization of motion performance indicators: Robot dexterity:

[0065] The robot's dexterity D m Defined as condition number k m The reciprocal:

[0066]

[0067] Where, k m The condition number k(J(θ)) of the Jacobian matrix: k m =k(J(θ));

[0068] Based on the singular value analysis theory of matrices, the Jacobian matrix J(θ) of the robot in any position is analyzed by singular value decomposition:

[0069] J(θ) = U∑V;

[0070] Where, U∈R m×m , V∈R n×nU and V are both orthogonal matrices, and ∑ is:

[0071]

[0072] Among them, (σ1,σ2,…,σ m Let σ₁ ≥ σ₂ ≥ … ≥ σ₃ be the singular values ​​of J(θ). m ≥0, σ1 is the largest singular value, σ m It is the smallest singular value. The relationship between the condition number k(J) of the Jacobian matrix and the singular values ​​is:

[0073]

[0074] The condition number of the Jacobian matrix takes values ​​in the range 1 ≤ k ≤ ∞. When the condition number of the robot is 1, D m =100%, at which point the robot's dexterity is at its highest, and all singular values ​​are equal.

[0075] The motion of the robot's end effector along the x and y directions in the world coordinate system {W} was simulated using MATLAB software to obtain the dexterity during the motion.

[0076] S34: The motion capability under the current milling pose and the adjustment range between adjacent milling poses characterize the robot's motion stability.

[0077] Preferably, the specific process for characterizing the area of ​​the stiffness anisotropy ellipse in step S32 is as follows:

[0078] S321: In the world coordinate system {W}, when the robot's end effector pose is R, the three vector directions on its end effector face are:

[0079]

[0080] Let the end-effector stiffness ellipsoid coordinate system be {E}, with its origin coinciding with the end-effector coordinate system {F}. In the world coordinate system {W}, the homogeneous transformation matrix is:

[0081]

[0082] The attitude of the end-effector stiffness ellipsoid in the world coordinate system {W} is as follows:

[0083]

[0084] Under the world coordinate system {W}, the equation of the ellipsoid is:

[0085]

[0086] S322: In the end-effector stiffness ellipsoidal coordinate system {E}, the pose of the robot's end-effector coordinate system {F} is calculated as follows:

[0087]

[0088] The x and y vectors of the terminal coordinate system are expressed in the ellipsoidal coordinate system {E} as follows:

[0089]

[0090] S323: The elliptical surface cut off by the end face of the stiffness ellipsoid is obtained by calculating the intersection point of the x and y vectors of the end coordinate system with the ellipsoid.

[0091] Preferably, the specific process of obtaining the elliptical surface cut off by the end face of the stiffness ellipsoid in step S323 based on the intersection point obtained by calculating the x and y vectors of the end coordinate system with the ellipsoid is as follows:

[0092] Let the intersection points of the x and y vectors with the ellipsoid be respectively and but:

[0093]

[0094] Solve the following equations simultaneously:

[0095]

[0096]

[0097] Solving for:

[0098]

[0099] Calculate the area of ​​the stiffness anisotropic ellipse based on the compliance coefficient on the robot's end face:

[0100]

[0101] Preferably, the specific process of step S4 is as follows:

[0102] S41: Implement an improved multi-leader update mechanism: Identify the current Pareto front solution set (PF) through non-dominated sorting, and then introduce a scaling factor α to design a dynamic selection rule for the number of leaders k.

[0103] N l =min(α·card(PF),50);

[0104] S42: Setting the adaptive speed update strategy: The original speed update strategy is optimized by introducing dynamic inertia weights and adaptive acceleration constants. The optimized adaptive speed update formula is as follows:

[0105] v i (t+1)=w·v(t)+c1·r1·(p * (t)-x i (t))+c2·r2·(g * (t)-x i (t));

[0106] Among them, v i (t+1) is the velocity of the i-th particle in the (t+1)-th iteration, w is the inertia weight, controlling the influence of the particle's current velocity, c1 and c2 are acceleration constants, controlling the influence of the particle's individual cognition and social learning, respectively, r1 and r2 are random numbers in the range [0,1], used to increase the randomness of the search, p * (t) represents the historical best position of the i-th particle, g * (t) is the globally optimal position, x i (t) is the position of the i-th particle in the t-th iteration;

[0107] S43: Introduce dynamic inertia weights and adaptive speedup constants to further optimize algorithm performance;

[0108] S44: Combining the traditional PSO search strategy with the MCMLWOA position update mode, a hybrid position update mechanism is formed: the random search of MO-MLSWOAR is replaced by the position update method of the PSO algorithm.

[0109]

[0110]

[0111] S45: The optimal milling posture for the corresponding milling point is obtained based on the MO-MLSWOAR optimization algorithm.

[0112] Preferably, the specific process of introducing dynamic inertia weights and adaptive speedup constants to further optimize algorithm performance in step S43 is as follows:

[0113] The inertial weight w, adaptive acceleration constants c1 and c2 are dynamically adjusted in a linear decreasing manner to balance global and local search capabilities.

[0114]

[0115] Among them, w max and w min These represent the maximum and minimum values ​​of the inertia weight, respectively; T is the maximum number of iterations; t is the current number of iterations; and c is the maximum and minimum values ​​of the inertia weight. 1,max and c 1,min These are the maximum and minimum values ​​of the acceleration constant c1, respectively. 2，max and c2，min These are the maximum and minimum values ​​of the acceleration constant c2, respectively.

[0116] The beneficial effects of this invention include:

[0117] The milling pose optimization method based on redundant milling axes provided by this invention firstly achieves redundancy of the milling cutter pose kinematics by vertically mounting the milling spindle on the sixth axis of the robot using six degrees of freedom. The rotation of the milling cutter coordinate system z-axis (i.e., milling spindle deflection) is described by the redundant degrees of freedom, which breaks through the limitations of traditional five-degree-of-freedom milling operations. This allows the robot to adjust the milling cutter posture within a more flexible range, adapting to the needs of diverse processing scenarios such as complex curved surfaces and irregular structures.

[0118] Secondly, based on the PCA method, the normal vectors of the points on the surface are solved using a point cloud model to determine the direction vectors of the milling redundant axes. This process involves steps such as selecting neighborhood points using Euclidean distance, constructing the covariance matrix, and solving for eigenvalues ​​and eigenvectors to accurately extract the local geometric features of the surface. This ensures a high degree of matching between the direction of the redundant axes and the characteristics of the machined surface, providing reliable basic data for subsequent pose optimization.

[0119] Furthermore, the established pose optimization objectives encompass multiple key dimensions: enhancing the robot's load-bearing and adaptability by characterizing the angle between the external load and the principal stiffness direction; optimizing stiffness performance to reduce machining deformation by using the end-effector stiffness ellipsoid volume and anisotropic elliptical area; ensuring motion flexibility through robot dexterity indices; and simultaneously focusing on the stability between adjacent poses to reduce equipment wear. This multi-objective collaborative optimization effectively balances machining quality, equipment lifespan, and operational efficiency.

[0120] Finally, the MO-MLSWOAR multi-objective optimization algorithm, by improving the leader update mechanism (dynamically selecting the number of leaders), the adaptive velocity update strategy (dynamically adjusting the inertia weight and acceleration constant), and the hybrid position update mechanism (integrating different search strategies), balances the diversity and convergence of the Pareto front solution set, avoids getting trapped in local optima, and can quickly find the best milling posture in a wide range of angles from -180° to 180°, significantly improving optimization efficiency and accuracy. Attached Figure Description

[0121] Figure 1 This is a schematic diagram illustrating the calculation of the end-stiffness ellipsoid volume of the present invention.

[0122] Figure 2 This is a schematic diagram of the elliptical area of ​​the end stiffness of the present invention.

[0123] Figure 3 This is a graph showing how the dexterity of the present invention changes with movement.

[0124] Figure 4This is a thermal diagram of the planar dexterity of the present invention.

[0125] Figure 5 This is a comparison chart of the optimization algorithms of this invention.

[0126] Figure 6 This is a milling surface quality analysis diagram of the present invention.

[0127] Figure 7 This is a graph showing the change of fitness values ​​with displacement under various experimental conditions of the present invention.

[0128] Figure 8 This is a graph showing the changes in dexterity and stiffness with displacement under various experimental conditions of the present invention. Detailed Implementation

[0129] The following is in conjunction with the appendix Figures 1 to 8 The present invention will be further described in detail below:

[0130] Example 1

[0131] See appendix Figure 1 As shown, a milling pose optimization method based on redundant milling axes includes the following steps:

[0132] S1: The milling spindle is vertically mounted on the sixth axis of the robot, so that the robot uses six degrees of freedom to determine the milling pose, realizing the redundancy of the kinematics of the milling cutter pose, and using the redundant degrees of freedom to describe the rotation of the milling cutter coordinate system z-axis, that is, the deflection of the milling spindle.

[0133] S2: Solving the milling redundant axis direction vector based on PCA method;

[0134] S3: Set the pose optimization objective for a six-DOF articulated robot in a typical five-DOF milling operation;

[0135] S4: The MO-MLSWOAR multi-objective optimization algorithm is used to optimize the angle within the range of -180° to 180° to obtain the optimal milling posture for the corresponding milling point.

[0136] A typical milling task requires five degrees of freedom (DOF), with three DDFs used to locate the tool's center and two DDFs used to determine the orientation of the milling cutter axis. For a six-DOF industrial robot, five DDFs control the tool's pose. If the milling spindle coincides with the robot's sixth axis, the rotational motion of the sixth axis will not affect the milling cutter's pose, resulting in the loss of the sixth DDF's function. To solve this problem, the milling spindle can be mounted vertically on the robot's sixth axis. This mounting method allows the robot to determine the milling pose using six DDFs, thus achieving DDF redundancy in the milling cutter's kinematics. This redundant DDF can be used to describe the rotation of the milling cutter's z-axis, i.e., the deflection of the milling spindle. This invention designs multiple performance indicators and a multi-objective optimization function. By integrating motion performance indicators, stiffness performance indicators, and robot joint motion constraints, a multi-objective optimization algorithm is used to optimize the angle within the range of -180° to 180° to obtain the optimal milling posture for any milling point.

[0137] In this embodiment, the rotation of the milling cutter coordinate system z-axis in step S1 is described with redundant degrees of freedom, that is, the specific process of the milling spindle deflection is as follows:

[0138] S11: In a milling task, the pose of the milling cutter coordinate system {M} in the base coordinate system {B} is described as a homogeneous transformation matrix:

[0139]

[0140] in, Each column in the table represents the projection of the coordinate axes of the milling cutter coordinate system {M} onto the base coordinate system {B}, and the rotation axis of the milling cutter. The unit vector in the base coordinate system {B} is:

[0141]

[0142] S12: Redundancy of degrees of freedom based on the milling cutter pose, the milling cutter coordinate system {M} revolves around the milling cutter's rotation axis. For rotation, if the rotation angle is γ, according to Rodrigues' rotation formula, the rotation matrix is... Represented as:

[0143]

[0144] Where I is the identity matrix, and K is the antisymmetric matrix of the rotation axis:

[0145]

[0146] S13: New pose of the milling cutter coordinate system {M} after rotation for:

[0147]

[0148] S14: Based on the relative poses between coordinate systems and the robot's kinematic equations, establish the relationship between the robot's joint angles and the rotation angle of the milling cutter coordinate system {M}:

[0149]

[0150] S15: When the robot performs a typical milling task, in the world coordinate system {W}, the pose of the robot cutting system is described as the milling position, the milling spindle direction, and the deflection of the milling spindle:

[0151]

[0152] in, For the milling position, a m γ represents the direction of the milling spindle, and γ is the angle of rotation of the milling coordinate system around the milling spindle.

[0153] When the robot joints are not at their limits, the milling cutter coordinate system {M} can rotate around the z-axis. This characteristic allows the robot to maintain the same milling position and direction while adopting various different configurations. Since the robot's end effector stiffness is pose-dependent, changing the robot's configuration directly affects the magnitude and direction of the end effector's principal stiffness. By adjusting different milling spindle rotation angles γ, the change in the robot's end effector's principal stiffness can be actively controlled, thereby optimizing the robot's milling pose.

[0154] Example 2

[0155] Building upon Example 1, when using an end mill for cutting, the milling cutter needs to be perpendicular to the machined surface. When the cutter is perpendicular to the machined surface, all cutting edges involved in the cutting can evenly distribute the cutting load, forming a stable cutting geometry. When there is an angle between the cutter and the machined surface, it will cause excessive load on one side of the cutting edge, which will not only accelerate tool wear and shorten tool life, but also cause cutting vibration, resulting in deterioration of surface roughness and even dimensional deviations. The direction perpendicular to the machined surface is essentially the normal vector direction of the machined surface. PCA (Principal Component Analysis) is a data dimensionality reduction technique that projects data by finding the direction of maximum variance in the data. It is commonly used for data compression, feature extraction, and normal vector calculation tasks for 3D point cloud data.

[0156] In step S2, the solution for the milling redundant axis direction vector based on the PCA method is achieved by solving the normal vector of a point p(x, y, z) on a point cloud model obtained from a freeform surface. The specific process is as follows:

[0157] S21: Selecting a Field:

[0158] The neighborhood points p1, p2, ..., p of p are selected using Euclidean distance. n The neighboring points satisfy:

[0159]

[0160] Where r is the radius of the sphere's neighborhood;

[0161] Calculate the mean vector of the neighborhood points:

[0162]

[0163] Use the set of neighboring points as sample D:

[0164]

[0165] S22: Constructing the Covariance Matrix: Covariance is an important concept in statistics and multivariate analysis, reflecting the degree of variation of data points in various directions. For a dataset containing the coordinates of n three-dimensional points, the covariance matrix is ​​a 3x3 matrix.

[0166] Center sample D:

[0167]

[0168] Its covariance matrix ∑ is:

[0169]

[0170] in, It is the variance of x; Let x and y be the covariance values, and so on, to solve for the covariance matrix ∑.

[0171] S23: Solving for eigenvalues ​​and eigenvalue vectors:

[0172] The eigenvalues ​​and eigenvectors of the covariance matrix ∑ satisfy:

[0173] ∑v i =λ i v i ;

[0174] Eigenvalues ​​are solutions to the characteristic equation of the covariance matrix Σ.

[0175] det(Σ-λ i I) = 0;

[0176] S24: For a 3×3 covariance matrix, there are three eigenvalues ​​λ1, λ2, and λ3. The magnitude of the eigenvalue represents the variance of the data in the direction of the corresponding eigenvector. For normal vector calculation, the eigenvector corresponding to the smallest eigenvalue reflects the direction of least change of the point cloud data in that local region. Therefore, the eigenvector with the smallest eigenvalue is selected as the normal vector.

[0177] λ n =min(λ1, λ2, λ3);

[0178] The eigenvector v is obtained by solving the following equation. n :

[0179] (∑-λ n I)v n =0,

[0180] Feature vector v n That is, the normal vector at point p.

[0181] Example 3

[0182] Based on Example 1 or Example 2,

[0183] Because milling tasks have redundant degrees of freedom, a six-joint robot can maintain the same milling position and orientation while exhibiting multiple different configurations during milling operations. Step S1 has already summarized the expression of the robot-cutting system's milling pose as follows: Including milling positions The direction of the milling spindle a m The angle γ of the rotation of the milling coordinate system around the milling spindle. Different milling poses result in different robot configurations, which in turn causes changes in the robot's end effector stiffness.

[0184] With a fixed end effector position, the robot rotates around the z-axis of the world coordinate system {W}, resulting in dynamic changes in both the overall stiffness of the end effector and the overall stiffness of the end face. These two metrics serve as evaluation standards for end effector stiffness and are crucial for guiding milling pose planning. However, when the rotation angle θ = 180°, although the overall stiffness on the end face is relatively good, it is also at its worst. Therefore, the overall stiffness of the end effector cannot be solely judged based on the overall end effector stiffness. Regarding these two metrics, relying solely on the result of one metric to judge the overall performance of end effector stiffness has significant limitations. Therefore, when evaluating the quality of pose-dependent end effector stiffness, multiple evaluation metrics need to be considered comprehensively. This means that milling pose optimization is essentially a multi-objective optimization problem, requiring a trade-off between various stiffness and non-stiffness metrics to achieve optimal overall performance.

[0185] The specific process of setting the pose optimization target of the six-DOF articulated robot under a typical five-DOF milling operation in step S3 is as follows:

[0186] S31: The angle between the principal stiffness direction and the resultant force of the external load.

[0187] The principal stiffness direction is the direction of maximum stiffness at the robot's end effector. When the direction of the resultant external load coincides with this direction, the elastic deformation of the end effector is minimized. To ensure that the deformation generated by the robot during milling is as small as possible, the principal stiffness direction of the robot's end effector should be aligned with the direction of the resultant external load as much as possible. The external loads on the robot's end effector include the milling force F. m and the original load F o This includes flanges, electric spindles, and force sensors, where the milling force F... m It is the original load F acting in the world coordinate system. o It acts on the end face, applying the external load F. ex Represented as:

[0188]

[0189] External load F ext With respect to the principal stiffness direction The included angle θ Fv1 for:

[0190]

[0191] S32: Indicator characterizing stiffness performance:

[0192] The end effector stiffness matrix K of the robot end It is a symmetric positive definite matrix, a property that guarantees the positive definiteness of the matrix eigenvalues ​​and the orthogonality of the eigenvectors. The complete end effector stiffness matrix K of a six-DOF robot. end ∈R 6×6 The stiffness distribution of the end effector in six degrees of freedom is described:

[0193]

[0194] Among them, K trans ∈R 3×3 K is the translational stiffness matrix, describing the end effector's ability to resist deformation under external forces during pure translational motion. rot ∈R 3×3 K is the rotational stiffness matrix, describing the end effector's ability to resist torque deformation during pure rotational motion. coupling ∈R 3×3 This is the coupling stiffness matrix, describing the degree of coupling between translational and rotational stiffness. Extracting K... end The 3×3 submatrix K in the top left cornertrans :

[0195]

[0196] In this matrix, the diagonal elements represent direct stiffness, indicating the stiffness along a specific axis. The diagonal elements also represent coupled stiffness, indicating that displacement in one direction leads to a related force in another direction. For the end-effector stiffness matrix K... end Perform singular value decomposition: K end =QΛQ T ;where: Λ=diag[λ1, λ2, λ3].

[0197] Λ is a diagonal matrix containing eigenvalues ​​(λ1>λ2>λ3>0). The eigenvalues ​​λ... i This indicates the magnitude of stiffness along the principal axis. The largest eigenvalue corresponds to the direction of maximum stiffness, also known as the principal stiffness direction, while the smallest eigenvalue corresponds to the direction of minimum stiffness.

[0198] Secondly, the orthogonal matrix Q: Q = [v1, v2, v3]. Q is a three-dimensional column vector v i The orthogonal matrix formed by the column vectors v i For the eigenvalue λ i The corresponding unit eigenvector.

[0199] The end effector stiffness ellipsoid is a core tool in robotics for describing the stiffness distribution of a mechanical system within the task space. Its geometry and mathematical properties directly reflect the robot's resistance to external loads and its anisotropy. The parameters of the end effector stiffness ellipsoid are as follows:

[0200] Table 1. End Stiffness Ellipsoid Parameters

[0201]

[0202] In terms of geometric representation, the symmetry of the stiffness matrix dictates that its geometric representation must be a quadratic surface, while its positive definiteness further restricts the possible closed surfaces to ellipsoids. The linear transformation effect of the Jacobian matrix geometrically manifests as an effect on the synthetic stiffness matrix K. t Rotation and scaling operations necessarily map a sphere or ellipsoid to a new ellipsoid. Therefore, indices related to the stiffness of the robot's end effector are inherently correlated with the stiffness ellipsoid.

[0203] The end-effector stiffness ellipsoid volume is an important parameter used to measure the overall stiffness performance of a robotic arm's end effector in the task space. The end-effector stiffness ellipsoid volume and the area of ​​the stiffness anisotropy ellipse are used to define the expression for the end-effector stiffness ellipsoid volume:

[0204]

[0205] The volume of the ellipsoid is negatively correlated with the overall stiffness of the robot's end effector; a smaller ellipsoid volume indicates higher overall stiffness, while a larger ellipsoid volume indicates lower overall stiffness. If the ellipsoid is close to a sphere, the stiffness distribution is uniform. Furthermore, the shape characteristics of the ellipsoid can intuitively reflect the anisotropy of the stiffness distribution: when the ellipsoid is close to a sphere, it indicates that the stiffness distribution is uniform in all directions; if the ellipsoid shows significant elongation deformation, it indicates that there are directional differences in stiffness performance; and an ellipsoid shape that extends significantly in a specific direction indicates that the stiffness in that axis is relatively weak.

[0206] Based on the ellipsoidal volume calculation formula, the overall stiffness of the robot end effector at a constant height of z = 1000 on a plane with the same pose is calculated. Furthermore, the overall stiffness of the robot end effector is calculated when it rotates one revolution around the z-axis of the world coordinate system {w} at a position of (2000, 0, 1000). A standard external load is set as F. std =[-15, 10, -96, 0, 0, 0] T Based on the above parameters, the calculated results of the end-effector stiffness ellipsoidal volume for the two working conditions are as follows: Figure 1 As shown:

[0207] Based on the analysis of the end-stiffness ellipsoidal volumetric thermogram, from Figure 1 (a) shows that as the robot's end effector moves further away from the base, the volume of the end effector stiffness ellipsoid gradually decreases, and the robot's end effector stiffness increases. In the upper right and lower right corners of Figure (a), the heatmap color lightens, indicating that the end effector stiffness begins to decrease when the end effector extends a certain distance. Considering the robot's motion performance in this pose, the robot's dexterity gradually decreases as the end effector moves further away from the base. Therefore, when planning the working area for the cutting task, it is not simply a matter of moving closer to or further away from the base; rather, multiple factors should be considered to rationally select the working area for the cutting task.

[0208] Depend on Figure 1 (b) It is known that when the robot end effector rotates around the z-axis of the world coordinate system {w} by an angle θ = 0° (360°), the volume of the robot end effector stiffness ellipsoid is the smallest, and the robot end effector stiffness is the largest; when θ = 180°, the volume of the robot end effector stiffness ellipsoid is the largest, and the robot end effector stiffness is the smallest. During the rotations from 0 to 180° and from 180° to 360°, the volume of the robot end effector stiffness ellipsoid changes continuously and symmetrically. Obviously, when the robot end effector is in this position, selecting the milling pose of θ = 0° (360°) can obtain the best overall stiffness performance. However, in practical applications, it is necessary to further consider the stiffness change of the robot end effector stiffness ellipsoid on the surface where the end effector acts as an external load, which can be measured by the following indicators.

[0209] End stiffness ellipsoid volume V FThe order of magnitude is approximately -15. To match the consistency of the numerical magnitude between the objective function and the fitness function, and to avoid optimization bias caused by differences in magnitude, the original data is processed as follows: idx(V F ) = log V F +17.

[0210] The specific process for characterizing the area of ​​the stiffness anisotropy ellipse in step S32 is as follows:

[0211] The area of ​​the stiffness anisotropy ellipse is a geometric index used to evaluate the compliance distribution of the end effector stiffness ellipsoid within a specific plane. When an external load is applied to the robot's end effector, the end effector acts as the force-bearing surface, and the stiffness ellipsoid is truncated by this plane, forming an elliptical surface.

[0212] S321: In the world coordinate system {W}, when the robot's end effector pose is R, the three vector directions on its end effector face are:

[0213]

[0214] Let the end-effector stiffness ellipsoid coordinate system be {E}, with its origin coinciding with the end-effector coordinate system {F}. In the world coordinate system {W}, the homogeneous transformation matrix is:

[0215]

[0216] The attitude of the end-effector stiffness ellipsoid in the world coordinate system {W} is as follows:

[0217]

[0218] Under the world coordinate system {W}, the equation of the ellipsoid is:

[0219]

[0220] S322: In the end-effector stiffness ellipsoidal coordinate system {E}, the pose of the robot's end-effector coordinate system {F} is calculated as follows:

[0221]

[0222] The x and y vectors of the terminal coordinate system are expressed in the ellipsoidal coordinate system {E} as follows:

[0223]

[0224] To obtain the elliptical surface truncated by the end face of the stiffness ellipsoid, it is necessary to calculate the intersection point of the x and y vectors of the end coordinate system with the ellipsoid. The length of the line segment from the origin to the intersection point represents the compliance coefficient in that direction. and These coefficients reflect the robot end effector's ability to resist deformation in the corresponding direction.

[0225] S323: The elliptical surface cut off by the end face of the stiffness ellipsoid is obtained by calculating the intersection point of the x and y vectors of the end coordinate system with the ellipsoid.

[0226] The specific process of obtaining the elliptical surface cut off by the end face of the stiffness ellipsoid in step S323 based on the intersection point obtained by calculating the x and y vectors of the end coordinate system with the ellipsoid is as follows:

[0227] Let the intersection points of the x and y vectors with the ellipsoid be respectively and but:

[0228]

[0229] Solve the following equations simultaneously:

[0230]

[0231] Solving for:

[0232]

[0233] Calculate the area of ​​the stiffness anisotropic ellipse based on the compliance coefficient on the robot's end face:

[0234]

[0235] In the analysis of anisotropic ellipses, the size of the ellipse area directly reflects the stiffness distribution characteristics of the robot's end effector plane in different directions. A smaller ellipse area indicates better overall stiffness performance of the robot's end effector plane in that pose. Conversely, a larger ellipse area indicates worse overall stiffness performance of the robot's end effector plane in that pose. and When the values ​​are relatively close, it indicates that the robot's stiffness distribution in all directions is relatively uniform, meaning that the robot's resistance to deformation in the x and y directions tends to be consistent. The robot will not deform significantly due to forces in any particular direction. and Significant differences in the values ​​indicate that the robot's stiffness on the end effector plane exhibits obvious anisotropy, with considerable differences in its resistance to deformation in different directions. In directions with lower stiffness, the robot's end effector may be more susceptible to external forces, resulting in larger displacements or vibrations.

[0236] The area of ​​the end effector stiffness ellipse is calculated for the constant height plane at height z = 1000 and the end effector position at (2000, 0, 1000) when the end effector rotates around the z-axis of the world coordinate system {W} for one revolution.

[0237] Based on the end-stiffness elliptical area thermogram analysis, from Figure 2 (a) It is known that as the robot's end effector moves further away from the base, the area of ​​the end effector stiffness ellipse gradually decreases, while the stiffness on the end effector plane increases, exhibiting a similar trend to the volume of the stiffness ellipsoid. From Figure 2 (b) shows that the area of ​​the robot end effector stiffness ellipse is largest when θ = 90° and θ = 270°, indicating that the overall stiffness performance of the end effector surface is weakest at this time. Within the rotation angle ranges of 0–90°, 90–270°, and 270–360°, the area of ​​the end effector stiffness ellipse changes relatively smoothly without significant fluctuations. However, to further analyze the end effector surface stiffness performance within these three rotation angle ranges, it is necessary to combine the compliance coefficient value to further determine whether there is significant anisotropy in the end effector surface ellipse. Furthermore, when θ = 180°, the volume of the end effector ellipsoid is smallest, but the elliptical area on the end effector surface is largest, indicating that the size of the elliptical area on the end effector surface cannot be judged solely based on the volume of the end effector ellipsoid. Based on the results given in Figure (b), in the planning of cutting tasks, poses with rotation angles of θ = 90° and θ = 270° should be avoided as much as possible to ensure that the robot end effector has sufficient stiffness and stability during machining.

[0238] S33: Characterization of motion performance indicators: Robot dexterity:

[0239] To avoid potential singular robot configurations during task planning and to effectively evaluate a robot's multi-degree-of-freedom motion capabilities in complex tasks, the motion dexterity index D... m Manipulability Measure (D&D) is widely used to quantitatively analyze the mobility and agility of robots. Its concept is related to the singular values ​​of the Jacobian matrix. m Defined as condition number k m The reciprocal:

[0240]

[0241] Where, k m The condition number k(J(θ)) of the Jacobian matrix: k m =k(J(θ));

[0242] Based on the singular value analysis theory of matrices, the Jacobian matrix J(θ) of the robot in any position is analyzed by singular value decomposition:

[0243] J(θ) = U∑V;

[0244] Where, U∈R m×m , V∈R n×n U and V are both orthogonal matrices, and ∑ is:

[0245]

[0246] Among them, (σ1,σ2,…,σ m Let σ₁ ≥ σ₂ ≥ … ≥ σ₃ be the singular values ​​of J(θ). m ≥0, σ1 is the largest singular value, σ m It is the smallest singular value. The relationship between the condition number k(J) of the Jacobian matrix and the singular values ​​is:

[0247]

[0248] The condition number of the Jacobian matrix takes values ​​in the range 1 ≤ k ≤ ∞. When the condition number of the robot is 1, D m =100%, at which point the robot's dexterity is at its highest, and all singular values ​​are equal.

[0249] Based on the above indicators, the robot's motion performance within the workspace can be analyzed. MATLAB software is used to simulate the robot's end effector's motion along the x and y directions in the world coordinate system {W}, obtaining the dexterity during the motion process.

[0250] To align with the milling cutter coordinate system, the robot's end effector orientation was determined as follows:

[0251]

[0252] The remaining motion path parameters are shown in the table below:

[0253] Table 2 Motion Path Parameters

[0254] coordinate x-direction movement y-direction movement Starting point (x, y) (1200,0) (-1000,2000) End point (x, y) (2500,0) (1000,2000) height z 800,1000,1200,1400 800,1000,1200,1400

[0255] The changes in dexterity at different heights when moving along the x-direction and along the y-direction are shown in the figure:

[0256] Depend on Figure 3It can be seen that among the various height values ​​set, the robot's dexterity is optimal at z=800. As the height increases, the robot's dexterity gradually decreases. During movement along the x-direction, the dexterity initially increases and then decreases. This indicates that dexterity is not always directly proportional to the proximity of the robot to its base. During movement along the y-direction, the dexterity at x=0 is optimal at z=800 and z=1000. As the height increases, the dexterity at x=0 changes from a maximum to a minimum, possibly because the robot's end effector is closer to a singular shape at x=0. Heatmaps of the robot's end effector's dexterity at four planes of equal height are plotted. Each heatmap uses the x and y directions as coordinate axes, and color gradients represent the magnitude of dexterity, with warm-toned areas representing high dexterity and cool-toned areas representing low dexterity.

[0257] according to Figure 4 Analysis of the heatmap of dexterity on a plane with equal height shows that the robot's dexterity gradually decreases as the plane height increases, and the high-dexterity region on the plane gradually moves towards the left and right sides of the base. This is consistent with the conclusions of the previous analysis: when the robot's end effector is at x=0 and in the region around x=0, the robot's shape and position are closer to a singular shape and position. Therefore, when the robot is machining a workpiece with a certain height, if the surface of the part is relatively complex, it can be fixed at a position off the x-axis to ensure that the robot achieves sufficient dexterity to meet the cutting task.

[0258] S34: The motion capability under the current milling pose and the adjustment range between adjacent milling poses characterize the robot's motion stability.

[0259] In robotic milling operations, in addition to the aforementioned stiffness and kinematic performance indicators, special attention must be paid to the robot's stability during the milling process. The stability of robot operation is mainly reflected in two aspects: the motion capability in the current milling pose and the adjustment range between adjacent milling poses. First, it is essential to ensure that the robot has sufficient motion capability in the current pose. Insufficient motion capability manifests as joint velocities failing to meet the end effector's operating speed. With a constant end effector speed, insufficient motion capability will lead to joint speed exceeding limits, resulting in overload and damage to the robot's joint motors. Second, the adjustment range between adjacent milling poses must be controlled within a reasonable range. Excessive adjustment range will cause significant changes in joint velocity within a short period. Since industrial robot links typically have large mass and inertia, rapid joint movements will trigger significant mechanical vibrations, leading to instability during robot movement and ultimately damage to the tool or workpiece. Therefore, sufficient motion capability is reflected in the constraint on robot joint velocity, while the adjustment range between adjacent milling poses is reflected in the constraint on robot joint acceleration. While optimizing the above indicators, it is crucial to ensure the robot's motion stability.

[0260] If the milling speed is v, then the robot's end effector speed is:

[0261] The joint velocities can be obtained using the inverse Jacobian matrix:

[0262] The minimum singular value σ of the Jacobian matrix m It can be used as an indicator to control the upper limit of joint speed:

[0263]

[0264] In a milling task, the robot travels from point A to point B. Given a sufficiently dense network of planned points, adjacent milling points... and The distance between them is:

[0265] If the feed rate is v f Then the time t between point A and point B is:

[0266] Therefore, adjacent poses A and Pose B The average joint acceleration between them is:

[0267]

[0268] It is necessary to ensure that the average acceleration between adjacent poses is less than the robot's rated acceleration during operation.

[0269]

[0270] Based on the analysis of the above design variables, constraints, and optimization objectives, the mathematical model for pose optimization during robot milling is established as follows:

[0271]

[0272] Among them, R i This represents the robot end effector's posture when the milling spindle rotates around the redundant axis at the milling position point, and T is the homogeneous transformation matrix of the robot end effector at this time. The fitness value of the multi-objective optimization function is the angle θ between the external load and the principal stiffness. Fv1 End-stiffness ellipsoid volume index idx(V F ), terminal ellipse area index idx(S F ) and dexterity index idx (D m A linear combination of ) is used, and σ is selected based on the degree of influence of the objective on the optimization result. i The optimization objective is to minimize the fitness function fit(T). i The pose that minimizes the fitness value is the target pose T. target To ensure the robot's smooth operation during milling tasks, certain limitations are placed on joint rotational speeds and average joint accelerations. The joint rotational speeds must meet certain requirements. Joint acceleration must meet the following requirements

[0273] MO-MLSWOAR is based on MCMLWOA, modifying the original multi-leader update mechanism and incorporating the search strategy of the PSO algorithm. Furthermore, it improves the original speed update strategy of the PSO algorithm into a speed-adaptive update strategy. The modified multi-leader update mechanism and speed-adaptive update strategy are as follows:

[0274] S41: The core of the multi-leader update mechanism is a dynamic adjustment strategy for the number of leaders based on changes in the Pareto front solution set and the fitness distribution density. An improved multi-leader update mechanism is established: the current Pareto front solution set PF is identified through non-dominated sorting, and then a scaling factor α is introduced to design a dynamic selection rule for the number of leaders k.

[0275] N l =min(α·card(PF), 50);

[0276] Where card(PF) is the number of elements in the PF set. α is the leader number adjustment factor, a number ranging from (0,1), and its relationship with the fitness distribution density Cr at the current iteration number t is: Cr = 1 - α;

[0277] Fitness distribution density Cr is used to assess the degree of clustering of fitness values ​​among individuals in the t-th generation population, and its calculation formula is as follows:

[0278]

[0279] Where n is the population size, fit(X) i ) is the fitness value of the i-th individual in the population at the t-th iteration, fit max and fit min These represent the maximum and minimum fitness values ​​of the current population, respectively. If Cr is small, it indicates that the individual fitness distribution is dense, suggesting that the individual may be trapped in a local optimum. It is necessary to increase α to introduce more leaders to enhance search diversity. If Cr is large, the fitness distribution is dispersed, and the number of leaders can be appropriately reduced to decrease the computational resources consumed in calculating the leader direction and weight.

[0280] S42: Setting the Adaptive Velocity Update Strategy: In the PSO algorithm, the velocity update formula is the core mechanism for adjusting the particle's motion direction and velocity. The adaptive velocity update strategy uses the basic PSO algorithm as a reference, and performs a secondary optimization of the original velocity update strategy by introducing dynamic inertia weights and an adaptive acceleration constant. The optimized adaptive velocity update formula is as follows:

[0281] v i (t+1)=w·v(t)+c1·r1·(p * (t)-x i (t))+c2·r2·(g * (t)-x i (t));

[0282] Among them, v i (t+1) is the velocity of the i-th particle in the (t+1)-th iteration, w is the inertia weight, controlling the influence of the particle's current velocity, c1 and c2 are acceleration constants, controlling the influence of the particle's individual cognition and social learning, respectively, r1 and r2 are random numbers in the range [0,1], used to increase the randomness of the search, p * (t) represents the historical best position of the i-th particle, g * (t) is the globally optimal position, x i (t) is the position of the i-th particle in the t-th iteration;

[0283] S43: Introducing dynamic inertia weights and adaptive speedup constants to further optimize algorithm performance:

[0284] The inertial weight w, adaptive acceleration constants c1 and c2 are dynamically adjusted in a linear decreasing manner to balance global and local search capabilities.

[0285]

[0286] Among them, w max and w min These represent the maximum and minimum values ​​of the inertia weight, respectively; T is the maximum number of iterations; t is the current number of iterations; and c is the maximum and minimum values ​​of the inertia weight. 1，max and c 1,min These are the maximum and minimum values ​​of the acceleration constant c1, respectively. 2，max and c 2，min These are the maximum and minimum values ​​of the acceleration constant c2, respectively.

[0287] The MO-MLSWOAR algorithm maintains the same overall process as MCMLWOA, but differs in its position update method. MCMLWOA's position update mode follows the three basic strategies of the traditional WOA algorithm: spiral position update, encircling prey, and random search. MO-MLSWOAR combines the PSO search strategy with the MCMLWOA position update mode, forming a hybrid position update mechanism. The traditional PSO algorithm's position update method is as follows:

[0288] x i (t+1)=x i (t)+v i (t+1).

[0289] S44: Combining the traditional PSO search strategy with the MCMLWOA position update mode, a hybrid position update mechanism is formed: the random search of MO-MLSWOAR is replaced by the position update method of the PSO algorithm.

[0290]

[0291] S45: The optimal milling posture for the corresponding milling point is obtained based on the MO-MLSWOAR optimization algorithm.

[0292] To verify the performance of the MO-MLSWOAR algorithm, this invention will list the PSO algorithm and GA algorithm respectively.

[69] The algorithm, the basic WOA algorithm, and the MO-MLSWOAR algorithm have several metrics, including GD (Generational Distance), SP (Spacing), the optimal individual, and the optimal fitness. GD is used to measure the average distance between the found solution set and the true Pareto front.

[0293]

[0294] Where N is the size of the solution set, d iGD represents the minimum distance from the i-th solution to the true Pareto front, where p=2 represents the Euclidean distance. A smaller GD value indicates that the solution set is closer to the true Pareto front, reflecting the algorithm's convergence performance. On the other hand, SP is used to evaluate the uniformity of the distribution within the algorithm's solution set, and its calculation formula is:

[0295]

[0296] Where, d i The distance from the i-th solution in the solution set to its nearest neighbor solution is given by the given information. This is the average of these distances. A smaller SP value indicates a more uniform distribution of the solution set, reflecting the algorithm's performance in terms of solution set diversity. Using the above algorithm, given the design variables and multi-objective optimization conditions, the rotation angle of the redundant milling axis at the same milling position is optimized, yielding... Figure 5 Results shown:

[0297] It can be seen that the MO-MLSWOAR algorithm performs best in overall performance, with significantly better convergence and solution set uniformity than other algorithms. PSO and WOA perform similarly in convergence and uniformity but are slightly inferior to MO-MLSWOAR, while GA has the highest GD and SP values, indicating that its convergence speed is slower and the solution set distribution is not uniform enough. In terms of optimal solution quality, MO-MLSWOAR, PSO, and WOA all achieve the same optimal fitness value, but MO-MLSWOAR's optimal position is more accurate than PSO and WOA; in contrast, GA's optimal position is more deviated, and its fitness value is also slightly worse. Overall, MO-MLSWOAR demonstrates higher solution accuracy and stability in multi-objective optimization problems.

[0298] To verify the effectiveness of the multi-objective pose optimization algorithm, this invention designed and conducted multiple milling experiments. After the experiments, wire cutting was used to sample the machined surface. Subsequently, the milled surface was inspected using an Alicona InfiniteFocus G5 3D surface measuring instrument. Through the morphology reconstruction of the milled surface, the morphological features of the milling texture and parameters such as surface roughness can be accurately obtained. Table 3 shows the milling experiment conditions and parameters:

[0299] Table 3 Experimental conditions and parameters

[0300]

[0301] Based on the parameters of Experiment 3, the MO-MLSWOAR multi-objective optimization algorithm was used to optimize the milling pose of a 300mm long straight milling path with a step size of 1mm. The z-axis rotation angle and joint angle data were recorded with a node at 50mm, as shown in the table below:

[0302] Table 4. Rotation angles and joint angles (deg)

[0303] y coordinate θ J1 J2 J3 J4 J5 J6 -100 31.84543 -11.4006 -12.2068 -31.3369 -53.7798 148.13 35.78697 -150 30.34773 -12.8242 -12.3491 -31.1648 -53.7238 148.0691 35.77778 -200 28.83771 -14.2324 -12.5482 -30.9236 -53.6477 147.9846 35.76629 -250 27.32536 -15.6254 -12.8056 -30.6107 -53.5614 147.8819 35.75827 -300 25.78056 -16.9955 -13.1128 -30.236 -53.428 147.7397 35.73319 -350 24.23914 -18.3477 -13.4773 -29.7892 -53.288 147.5802 35.71402 -400 22.68954 -19.6784 -13.8952 -29.2746 -53.1266 147.3944 35.69267

[0304] See also Figure 6 The milled surface quality analysis diagram shown below:

[0305] Overall, all parameters in Experiment 1 were significantly higher than those in Experiments 2 and 3. Comparing Experiments 1 and 2, the reductions in Ra, Rq, and Rt were 76.4%, 75.4%, and 65.2%, respectively. Comparing Experiments 1 and 3, the reductions in Ra, Rq, and Rt were 78.1%, 77.2%, and 69.1%, respectively. This significant difference indicates that in milling, the cross-axis mounting method is superior to the parallel-axis mounting method. Further comparing Experiments 2 and 3, Experiment 3 was superior in several key parameters: Ra decreased by 7.1%, Rz decreased by 10.9%, and Rsm decreased by 13.7%. Except for Rp, which was slightly higher than in Experiment 2, all other parameters were lower than in Experiment 2. The results show that by optimizing the milling pose, surface roughness can be further reduced, milling line spacing can be decreased, and thus the overall machining quality can be improved.

[0306] The fitness values ​​and dexterity and stiffness performance indices under various experimental conditions were calculated using a pose optimization mathematical model, such as... Figure 7 and Figure 8 As shown:

[0307] In terms of fitness, overall, Experiment 3's fitness value is better than Experiments 1 and 2, consistently showing a decreasing trend. At the start of milling, Experiment 2's fitness value is better than Experiment 1's; however, near the end of milling, Experiment 2's fitness value is higher than Experiment 1's. Theoretically, near the end of milling, the robot's end effector stiffness and overall motion performance under Experiment 1 should be better than under Experiment 2's performance. However, based on the experimental results, Experiment 2's milled surface quality is better than Experiment 1's, mainly due to the different milling spindle mounting methods. Regarding dexterity and stiffness performance indicators, the robot's dexterity varies little across experimental conditions. In terms of stiffness performance indicators, Experiment 3 exhibits the best stiffness on the end effector surface. Since the robot's end effector load acts directly on the end effector surface, although the overall stiffness under Experiment 3 is not as good as Experiment 2, the milled surface quality under Experiment 3 is the best. In conclusion, the overall fitness analysis shows that the milling quality under the cross-axis mounting method is significantly better than that under the parallel-axis mounting method. A comparison of dexterity and stiffness performance indicators shows that the overall stiffness performance of the milled end face is crucial to the surface quality of the milled surface.

[0308] In summary, this invention first analyzes the pose optimization objective of a six-DOF articulated robot in a typical five-DOF milling operation, namely, the rotation angle around the milling redundant axis. Secondly, multiple performance indicators and a multi-objective optimization function are designed. Combining motion performance indicators, stiffness performance indicators, and robot joint motion constraints, the MO-MLSWOAR multi-objective optimization algorithm is used to optimize the angle within the range of -180° to 180°, obtaining the optimal milling pose for the milling point. The performance of MO-MLSWOAR is compared with PSO, ordinary WOA, and GA using four indicators: GD, SP, optimal fitness, and optimal individual position. Experiments demonstrate that MO-MLSWOAR exhibits higher solution accuracy and stability in multi-objective optimization problems. After pose optimization, three sets of experiments are designed to compare the milling effects of the optimized and unoptimized robots. Milling is performed using the optimized pose, and compared to the unoptimized pose, the surface roughness, average peak-to-valley distance, and waviness parameters obtained by the optimized pose milling are significantly improved. Experiments have shown that, based on the use of cross-axis mounting, this optimization algorithm has a significant effect on improving the quality of robot milling.

Claims

1. A milling pose optimization method based on redundant milling axes, characterized in that, Includes the following steps: S1: The milling spindle is vertically mounted on the robot's sixth axis, allowing the robot to determine the milling pose using six degrees of freedom. This achieves redundancy in the kinematics of the milling cutter pose, and the redundant degrees of freedom describe the milling cutter coordinate system. z The rotation of the shaft, i.e., the deflection of the milling spindle; S2: Solving the milling redundant axis direction vector based on PCA method; S3: Set the pose optimization objective for a six-DOF articulated robot in a typical five-DOF milling operation; S4: The MO-MLSWOAR multi-objective optimization algorithm is used to optimize the angle within the range of -180° to 180° to obtain the optimal milling posture for the corresponding milling point. The specific process of step S4 is as follows: S41: Implement an improved multi-leader update mechanism: Identify the current Pareto front solution set (PF) through non-dominated sorting, and then introduce a scaling factor. α Number of design leaders k Dynamic selection rules: ； S42: Setting the adaptive speed update strategy: The original speed update strategy is further optimized by introducing dynamic inertia weights and adaptive acceleration constants. The optimized adaptive speed update formula is as follows: ； in, It is the first i The particle in the first t Speed ​​in +1 iteration w It is the inertial weight, which controls the influence of the particle's current velocity. c 1 and c 2 is the acceleration constant, which controls the influence of individual particle cognition and social learning. r 1 and r 2. Random numbers in the range [0,1] are used to increase the randomness of the search. p *( t ) is the first i The historical best position of each particle It is the globally optimal position. x i ( t ) is the first i The particle in the first t Position in the next iteration; S43: Introduce dynamic inertia weights and adaptive speedup constants to further optimize algorithm performance; S44: It adopts the spiral position update, prey encirclement, and random search of the traditional WOA algorithm, and combines the traditional PSO search strategy with the MCMLWOA position update mode to form a hybrid position update mechanism: the random search of MO-MLSWOAR is replaced by the position update method of the PSO algorithm. ； S45: The optimal milling posture for the corresponding milling point is obtained based on the MO-MLSWOAR optimization algorithm.

2. The milling pose optimization method based on milling redundant axes according to claim 1, characterized in that, In step S1, the milling cutter coordinate system is described with redundant degrees of freedom. z The specific process of shaft rotation, i.e., the deflection of the milling spindle, is as follows: S11: In a milling task, the pose of the milling cutter coordinate system {M} in the base coordinate system {B} is described as a homogeneous transformation matrix: ； in, , Each column in the table represents the projection of the coordinate axes of the milling cutter coordinate system {M} onto the base coordinate system {B}, and the rotation axis of the milling cutter. The unit vector in the base coordinate system {B} is: ； S12: Redundancy of degrees of freedom based on the milling cutter pose, the milling cutter coordinate system {M} revolves around the milling cutter's rotation axis. For rotation, if the rotation angle is γ, according to Rodrigues' rotation formula, the rotation matrix is... Represented as: ； in, I It is the identity matrix. K Here, γ is the antisymmetric matrix of the rotation axis, and γ is the rotation angle: ； S13: New pose of the milling cutter coordinate system {M} after rotation for: ； S14: Based on the relative poses between coordinate systems and the robot's kinematic equations, establish the relationship between the robot's joint angles and the rotation angle of the milling cutter coordinate system {M}: ； S15: When the robot performs a typical milling task, in the world coordinate system {W}, the pose of the robot cutting system is described as the milling position, the milling spindle direction, and the deflection of the milling spindle: ； in, For milling positions, γ represents the direction of the milling spindle, and γ is the angle of rotation of the milling coordinate system around the milling spindle.

3. The milling pose optimization method based on milling redundant axes according to claim 1, characterized in that, In step S2, the milling redundant axis direction vector is solved using the PCA method by obtaining a point cloud model on a freeform surface and solving for a point on the point cloud. This is achieved using normal vectors, and the specific process is as follows: S21: Selecting a Field: Selected by Euclidean distance p neighborhood points The neighboring points satisfy: ； in, r It is the radius of the sphere's neighborhood; Calculate the mean vector of the neighborhood points: ； The set of neighboring points is used as the sample. D : ； S22: Construct the covariance matrix: Sample D Centralization: ； Its covariance matrix for: ； in, ,yes x The variance; ,yes x , y The covariance value is obtained by solving for the covariance matrix. ; S23: Solving for eigenvalues ​​and eigenvalue vectors: covariance matrix The eigenvalues ​​and eigenvectors satisfy: ； eigenvalue λ i It is the covariance matrix The solution to the characteristic equation. v i For the feature vector: ； S24: For a 3×3 covariance matrix, there are three eigenvalues ​​λ1, λ2, and λ3. The magnitude of the eigenvalue represents the variance of the data in the direction of the corresponding eigenvector. For normal vector calculation, the eigenvector corresponding to the smallest eigenvalue reflects the direction of least change of the point cloud data in that local region. Therefore, the eigenvector with the smallest eigenvalue is selected as the normal vector. ； The eigenvectors are obtained by solving the following equations. v n : ， Feature vector v n That is p The normal vector at the point.

4. The milling pose optimization method based on milling redundant axes according to claim 1, characterized in that, The specific process of setting the pose optimization target of the six-DOF articulated robot under a typical five-DOF milling operation in step S3 is as follows: S31: The angle between the principal stiffness direction and the resultant force of the external load. External loads at the robot end effector include milling forces. F m and existing load F o This includes flanges, electric spindles, and force sensors, among which the milling force... It is the existing load acting in the world coordinate system. F o It acts on the end face, transferring the external load. F ext Represented as: ； External load F ext With respect to the principal stiffness direction The included angle for: ； S32: Characterizing stiffness performance indicators: end-stiffness ellipsoid volume and stiffness anisotropic elliptic area; setting the expression for end-stiffness ellipsoid volume: ； S33: Characterization of motion performance indicators: Robot dexterity: robot dexterity D m Defined as condition number k m The reciprocal: ； in, k m It is the condition number of the Jacobian matrix. : ; Based on the singular value analysis theory of matrices, the Jacobian matrix of the robot in arbitrary positions is determined. Analysis using singular value decomposition: ； in, , , U and V They are all orthogonal matrices. for: ； in, for singular values, , It is the largest singular value. It is the smallest singular value, and the condition number of the Jacobian matrix. The relationship with singular values ​​is as follows: ； Using MATLAB software, the robot's end effector is analyzed along the world coordinate system {W}. x direction and y Simulate the directional motion process to obtain the dexterity during the motion; S34: The motion capability under the current milling pose and the adjustment range between adjacent milling poses characterize the robot's motion stability.

5. The milling pose optimization method based on milling redundant axes according to claim 4, characterized in that, The specific process for characterizing the area of ​​the stiffness anisotropy ellipse in step S32 is as follows: S321: In the world coordinate system {W}, the robot's end-effector pose is... R At that time, the three vector directions on its end face are: ； Let the end-effector stiffness ellipsoid coordinate system be {E}, with its origin coinciding with the end-effector coordinate system {F}. In the world coordinate system {W}, the homogeneous transformation matrix is: ； The attitude of the end-effector stiffness ellipsoid in the world coordinate system {W} is as follows: ； Under the world coordinate system {W}, the equation of the ellipsoid is: ； S322: In the end-effector stiffness ellipsoidal coordinate system {E}, the pose of the robot's end-effector coordinate system {F} is calculated as follows: ； End coordinate system x vectors and y The vector in the ellipsoidal coordinate system {E} is represented as: ； ； S323: Based on the calculation of the end coordinate system x vectors and y The intersection point obtained by the vector and the ellipsoid yields the elliptical surface cut off by the end face of the stiffness ellipsoid.

6. The milling pose optimization method based on milling redundant axes according to claim 5, characterized in that, In step S323, the coordinate system of the calculated end point is... x vectors and y The specific process of obtaining the elliptical surface cut off by the end face of the stiffness ellipsoid from the intersection point of the vector and the ellipsoid is as follows: make x vectors and y The intersection points obtained by the vector and the ellipsoid are respectively and ,but: ； ； Solve the following equations simultaneously: ； ； Solving for: ； ； Calculate the area of ​​the stiffness anisotropic ellipse based on the compliance coefficient on the robot's end face: 。 7. The milling pose optimization method based on milling redundant axes according to claim 1, characterized in that, The specific process of introducing dynamic inertia weights and adaptive speedup constants to further optimize algorithm performance in step S43 is as follows: Inertia weight w、 Adaptive acceleration constant c 1 and c 2. Dynamically adjust using a linear decreasing method to balance global and local search capabilities: ； ； ； in, and These are the maximum and minimum values ​​of the inertia weight, respectively. T The maximum number of iterations, t This represents the current iteration number. and These are acceleration constants. c The maximum and minimum values ​​of 1, and These are acceleration constants. c The maximum and minimum values ​​of 2.

Images