Milling pose optimization method based on milling redundant shaft

Abstract

The invention belongs to the technical field of robot pose optimization, and particularly relates to a milling pose optimization method based on a milling redundant shaft, which comprises the following steps: vertically mounting a milling main shaft on a sixth shaft of a robot, realizing degree-of-freedom redundancy of a milling cutter pose by using six degrees of freedom, and describing rotation of a z shaft of a milling cutter coordinate system by using redundant degrees of freedom; solving a normal vector of a curved surface point location through a point cloud model based on a PCA (Principal Component Analysis) method, and determining a milling redundant axis direction vector; then setting a pose optimization target of the six-degree-of-freedom joint robot under typical five-degree-of-freedom milling operation, wherein the pose optimization target covers an included angle between an external load and a main rigidity direction, a rigidity performance index, a movement performance index and movement stability; and finally, performing angle optimization in a range of-180 degrees to 180 degrees through an MO-MLSWOAR multi-objective optimization algorithm to obtain an optimal milling posture of a corresponding milling point location. According to the method, the redundant degree of freedom can be fully utilized, the curved surface machining features are accurately matched, multi-target collaborative optimization is achieved, and the milling quality and efficiency are improved.

Description

TECHNICAL FIELD

[0001] 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. BACKGROUND

[0002] In modern manufacturing, milling is an important material removal process, which is widely used in the forming and precision machining of complex parts. With the development of industrial robot technology, six-degree-of-freedom joint robots are gradually applied to the milling field due to their high flexibility and working range, in order to meet the machining needs of complex curved surfaces, special-shaped structures and other workpieces.

[0003] However, the traditional milling pose control method still has many limitations: on the one hand, some milling equipment or robot systems adopt five-degree-of-freedom configuration, which lacks redundant freedom, resulting in insufficient flexibility of milling tool pose adjustment, and it is difficult to achieve optimal pose adaptation in complex machining scenarios, especially when processing complex workpieces such as free-form surfaces, which may result in low machining precision and fast tool wear; on the other hand, even if a six-degree-of-freedom robot is used, the existing technology does not fully utilize the redundant freedom, and it is difficult to effectively achieve precise description and control of milling tool coordinate system rotation (i.e. milling spindle deflection) through redundant freedom, making it difficult to match the dynamic needs of different machining points.

[0004] In terms of solving the redundant axis direction vector, existing methods mostly rely on empirical values or simplified models, and lack accurate extraction of local geometric features of the machining surface. For example, in free-form surface milling, if the normal vector of each point on the surface cannot be accurately obtained, the redundant axis direction will not match the characteristics of the surface, which will affect the effectiveness of subsequent pose optimization.

[0005] The setting of the pose optimization target also has the problem of singleization. Traditional methods often only focus on a single performance indicator (such as stiffness or motion dexterity), while ignoring the need for multi-objective collaborative optimization. In actual milling process, the matching degree of the external load (including milling force and self-load) on the robot end to the main stiffness direction, the volume and anisotropy of the end stiffness ellipsoid, the motion dexterity of the robot and the smoothness between adjacent poses, etc. will all affect the machining quality and equipment life, and single target optimization cannot balance the overall performance.

[0006] In addition, existing multi-objective optimization algorithms have deficiencies in milling pose optimization: some algorithms have unreasonable leader mechanism design, which makes it difficult to balance the diversity and convergence of the Pareto frontier solution set; the speed update strategy lacks adaptive adjustment, which may result in slow convergence speed or falling into local optimum; at the same time, traditional algorithms have poor compatibility when integrating different search strategies, which affects the optimization efficiency and accuracy, and it is difficult to quickly find the best milling pose in a large range of angles from -180° to 180°.

[0007] Therefore, in view of the deficiencies of robot milling pose optimization in the prior art in terms of degree of freedom utilization, redundant axis solving, multi-objective coordination and optimization algorithm performance, there is an urgent need for a milling pose optimization method capable of fully utilizing redundant degrees of freedom, accurately extracting surface features, comprehensively considering multi-objective performance and efficiently searching for an optimal solution, so as to improve the milling quality and processing efficiency of complex workpieces. SUMMARY

[0008] The present application aims to provide a milling pose optimization method based on milling redundant axes to solve the technical problem of the deficiencies of robot milling pose optimization in the prior art in terms of degree of freedom utilization, redundant axis solving, multi-objective coordination and optimization algorithm performance.

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

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

[0011] S1: vertically installing a milling spindle on the sixth axis of a robot, so that the robot determines the milling pose with six degrees of freedom, achieving the kinematic redundancy of the milling tool pose and describing the rotation of the z-axis of the milling tool coordinate system, i.e. the deflection of the milling spindle, with redundant degrees of freedom;

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

[0013] S3: setting the pose optimization target of a six-degree-of-freedom joint robot under typical five-degree-of-freedom milling operations;

[0014] S4: performing angle optimization in the range of -180° to 180° through the MO-MLSWOAR multi-objective optimization algorithm to obtain the optimal milling attitude of the corresponding milling point.

[0015] Preferably, the specific process of describing the rotation of the z-axis of the milling tool coordinate system, i.e. the deflection of the milling spindle, with redundant degrees of freedom in step S1 is as follows:

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

[0017]

[0018] wherein, Each column in the matrix represents the projection of the coordinate axis of the milling tool coordinate system {M} in the base coordinate system {B}, and the rotation axis of the milling tool is The unit vector in the base coordinate system {B} is:

[0019]

[0020] S12: Based on the degree of freedom redundancy of the milling cutter pose, the milling cutter coordinate system {M} rotates around the rotation axis of the milling cutter , if the rotation angle is γ, according to the Rodrigues rotation formula, the rotation matrix is expressed as:

[0021]

[0022] Where I is the unit matrix, and K is the skew-symmetric matrix of the rotation axis:

[0023]

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

[0025]

[0026] S14: Based on the relative pose between the coordinate systems and the kinematics equation of the robot, the relationship between the joint angle of the robot and the rotation angle of the milling cutter coordinate system {M} is established:

[0027]

[0028] S15: When the robot performs a typical milling task, 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 the world coordinate system {W}:

[0029]

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

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

[0032] S21: Select the field:

[0033] Select the neighborhood points p1, p2,..., p n of p by Euclidean distance, which satisfy:

[0034]

[0035] Where r is the radius of the spherical neighborhood;

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

[0037]

[0038] Take the set of neighborhood points as sample D:

[0039]

[0040] S22: Construct the covariance matrix:

[0041] Center the sample D:

[0042]

[0043] Its covariance matrix ∑ is:

[0044]

[0045] Where, is the variance of x; is the covariance value of x, y, and so on, to solve the covariance matrix ∑;

[0046] S23: Solve the eigenvalue and eigenvector:

[0047] The eigenvalue and eigenvector of the covariance matrix Σ satisfy:

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

[0049] The eigenvalue is the solution of 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, λ3, the size of the eigenvalue represents the variance of the data in the corresponding eigenvector direction, for the normal vector calculation, the eigenvector corresponding to the smallest eigenvalue reflects the direction in which the point cloud data changes the least in the local area, therefore, select the eigenvector with the smallest eigenvalue as the normal vector:

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

[0053] Solve the eigenvector v n by the following equation:

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

[0055] The eigenvector v n is the normal vector at the point p.

[0056] Preferably, the specific process of setting the pose optimization goal of the six-degree-of-freedom joint robot under the typical five-degree-of-freedom milling operation in step S3 is as follows:

[0057] S31: Characterize the angle between the principal stiffness direction and the resultant external load:

[0058] The external load at the robot end includes the milling force F m and the original load F o , including the flange, the electric spindle, and the force sensor, wherein the milling force F m is acting in the world coordinate system, and the original load F o is acting on the end face, and the external load F ext is expressed as:

[0059]

[0060] The angle θ ext between the external load F Fv1 and the principal stiffness direction is:

[0061]

[0062] S32: Characterize the stiffness performance index: the volume of the end stiffness ellipsoid and the area of the stiffness anisotropic ellipse, and set the expression of the volume of the end stiffness ellipsoid:

[0063]

[0064] S33: Characterize the motion performance index: the robot dexterity:

[0065] The robot dexterity D m is defined as the inverse of the condition number k m :

[0066]

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

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

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

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

[0071]

[0072] where (σ1, σ2, …, σ m ) are singular values of J(θ), σ1≥σ2≥…≥σ m ≥0, σ1 is the largest singular value, and σ m is the smallest singular value, the condition number k(J) of the Jacobian matrix has the following relationship with the singular values:

[0073]

[0074] The condition number of the Jacobian matrix ranges from 1 to k, and when the condition number of the robot is 1, D m = 100%, at which time the dexterity of the robot is the highest, and each singular value is equal.

[0075] The motion process of the robot end along the x direction and the y direction in the world coordinate system {W} is simulated using MATLAB software, and the dexterity in the motion process is obtained.

[0076] S34: The motion ability of the current milling pose and the adjustment range between adjacent milling poses represent the motion stability of the robot.

[0077] Preferably, the specific process of representing the anisotropic elliptical area of stiffness in step S32 is as follows:

[0078] S321: When the robot end pose is R in the world coordinate system {W}, the three vector directions on the end face are:

[0079]

[0080] The end stiffness ellipsoid coordinate system is set as {E}, the origin of which coincides with the end coordinate system {F}, and the homogeneous transformation matrix in the world coordinate system {W} is:

[0081]

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

[0083]

[0084] The ellipsoid surface equation in the world coordinate system {W} is established as:

[0085]

[0086] S322: In the end stiffness ellipsoid coordinate system {E}, the pose of the robot end coordinate system {F} is calculated as:

[0087]

[0088] The expressions of the x vector and the y vector of the end coordinate system in the ellipsoid coordinate system {E} are:

[0089]

[0090] S323: Obtain the elliptical surface of the stiffness ellipsoid cut by the end face based on the intersection points of the x vector and the y vector of the end coordinate system and the ellipsoid.

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

[0092] Let the intersection points of the x vector and the y vector and the ellipsoid be and Then:

[0093]

[0094] Solve the following equations:

[0095]

[0096]

[0097] The solution is:

[0098]

[0099] According to the compliance coefficient on the end face of the robot, calculate the area of the stiffness anisotropic ellipse:

[0100]

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

[0102] S41: Set the improved multi-leader update mechanism: identify the current Pareto front solution set PF through non-dominated sorting, and then introduce a proportion factor a to design a dynamic selection rule for the number of leaders k:

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

[0104] S42: Set the speed adaptive update strategy: through the introduction of dynamic inertia weight and adaptive acceleration constant, the original speed update strategy is optimized twice, and the optimized speed adaptive update formula is:

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

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

[0107] S43: Introduce dynamic inertia weight and adaptive acceleration constant to further optimize the performance of the algorithm;

[0108] S44: Combine the search strategy of traditional PSO with the position update mode of MCMLWOA to form a hybrid position update mechanism: replace the random search of MO-MLSWOAR with the position update mode of PSO algorithm:

[0109]

[0110]

[0111] S45: Based on the MO-MLSWOAR optimization algorithm, the optimal milling posture of the corresponding milling point is obtained.

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

[0113] The inertia weight w, the adaptive acceleration constants c1 and c2 are dynamically adjusted by linearly decreasing to balance the global search and local search capabilities:

[0114]

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

[0116] The beneficial effects of the present application include:

[0117] The milling pose optimization method based on milling redundant axes provided by the present application firstly realizes the degree of freedom redundancy of the milling tool pose kinematics by vertically installing the milling spindle on the sixth axis of the robot, and describes the rotation of the z-axis of the milling tool coordinate system (i.e., the deflection of the milling spindle) with the redundant degrees of freedom, thereby breaking through the limitation of traditional five-degree-of-freedom milling operations, enabling the robot to adjust the milling tool pose in a more flexible range and adapt to the needs of complex curved surfaces, special-shaped structures and other diverse processing scenarios.

[0118] Secondly, the normal vector of the curved surface point is solved based on the PCA method through the point cloud model to determine the direction vector of the milling redundant axis. This process accurately extracts the local geometric features of the curved surface by steps such as filtering neighborhood points by Euclidean distance, constructing a covariance matrix, solving eigenvalues and eigenvectors, etc., to ensure that the direction of the redundant axis is highly matched with the characteristics of the processed curved surface, providing reliable basic data for subsequent pose optimization

[0119] Thirdly, the set pose optimization target covers multiple key dimensions: by representing the angle between the external load and the main stiffness direction, the load bearing and adaptability of the robot are improved; by optimizing the stiffness performance to reduce processing deformation through the volume of the end stiffness ellipsoid and the anisotropic ellipse area; by the robot dexterity index, the motion flexibility is ensured; at the same time, the smoothness between adjacent poses is concerned to reduce equipment wear and tear. Multi-objective collaborative optimization effectively balances processing quality, equipment life and operation efficiency.

[0120] Finally, the MO-MLSWOAR multi-objective optimization algorithm adopted balances the diversity and convergence of the Pareto front solution set by improving the leader update mechanism (dynamically selecting the number of leaders), the adaptive speed update strategy (dynamically adjusting the inertia weight and acceleration constant) and the hybrid position update mechanism (blending different search strategies), avoids falling into local optimum, and can quickly find the best milling pose in a large range of angles of -180° ~ 180°, significantly improving the optimization efficiency and accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0121] Figure 1 It is an end stiffness ellipsoid volume calculation schematic diagram of the present application.

[0122] Figure 2 It is an end stiffness ellipse area schematic diagram of the present application.

[0123] Figure 3 It is a dexterity change with movement diagram of the present application.

[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-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 can use six degrees of freedom to determine the milling pose. This achieves redundancy of the degrees of freedom in the kinematics of the milling cutter pose, and uses 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 pose between each coordinate system and the kinematics equation of the robot, the relationship between the joint angle of the robot and the rotation angle of the milling cutter coordinate system {M} is established:

[0149]

[0150] S15: When the robot performs a typical milling task, 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 the world coordinate system {W}:

[0151]

[0152] wherein, is the milling position, a m is 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 do not reach the limit, the milling cutter coordinate system {M} can rotate around the z-axis direction. This feature allows the robot to exhibit multiple configurations while maintaining the same milling position and direction. Since the end stiffness of the robot is pose-dependent, changing the configuration of the robot directly affects the size and direction of the end main stiffness. By adjusting the different milling spindle rotation angles γ, the change of the main stiffness of the robot end can be actively controlled, thereby optimizing the milling pose of the robot.

[0154] Embodiment 2

[0155] On the basis of embodiment 1, when using an end mill for cutting tasks, the milling cutter needs to be perpendicular to the machining surface. When the tool is perpendicular to the machining surface, all the cutting edges involved in cutting can evenly share the cutting load, forming a stable cutting geometric angle. When the tool has an inclination angle with the machining surface, the single-sided cutting edge of the tool will bear excessive load, which not only accelerates tool wear and shortens tool life, but also causes cutting vibration, resulting in surface roughness deterioration and even size out-of-tolerance. The direction perpendicular to the machining surface is essentially the direction of the normal vector of the machining surface. The PCA (Principal Component Analysis) method is a data dimensionality reduction technique that projects data by finding the direction with the maximum variance in the data. It is commonly used for data compression, feature extraction, and normal vector solving tasks for three-dimensional point cloud data.

[0156] The direction vector of the milling redundant axis in step S2 is solved based on the PCA method by obtaining a point cloud model from a free-form surface, and the normal vector of a point p(x, y, z) on the point cloud is solved. The specific process is as follows:

[0157] S21: Select the field:

[0158] Select the neighborhood points p1, p2,..., p of p by the Euclidean distance, which satisfy: n

[0159]

[0160] where r is the radius of the spherical neighborhood;

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

[0162]

[0163] Take the set of neighborhood points as the sample D:

[0164]

[0165] S22: Construct the covariance matrix: Covariance is an important concept in statistics and multivariate analysis, reflecting the degree of change of data points in each direction. For a data set containing n three-dimensional point coordinates, the covariance matrix is a 3*3 matrix.

[0166] Center the sample D:

[0167]

[0168] Its covariance matrix ∑ is:

[0169]

[0170] where, is the variance of x; is the covariance value of x, y, and so on, to solve the covariance matrix ∑;

[0171] S23: Solve the eigenvalue and eigenvector:

[0172] The eigenvalue and eigenvector of the covariance matrix ∑ satisfy:

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

[0174] The eigenvalue is the solution of the characteristic equation of the covariance matrix Σ:

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

[0176] ​S24: For a 3x3 covariance matrix, there are three eigenvalues λ1, λ2, λ3, the size of the eigenvalue represents the variance of the data in the corresponding eigenvector direction, for normal vector calculation, the eigenvector corresponding to the smallest eigenvalue reflects the direction in which the point cloud data in the local region changes the least, therefore, the eigenvector with the smallest eigenvalue is selected as the normal vector:

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

[0178] The eigenvector v is solved by the following equation n :

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

[0180] The eigenvector v is solved by the following equation n is the normal vector at p.

[0181] Example 3

[0182] On the basis of example 1 or example 2,

[0183] Due to the characteristics of redundant degrees of freedom of the milling task, the six-joint robot can present a variety of different configurations while maintaining the same milling position and direction when performing the milling task. In step S1, the expression form of the milling pose of the robot-cutting system has been summarized as including the milling position the direction a of the milling spindle m and the angle γ of the milling coordinate system rotating around the milling spindle. Due to different milling poses, the robot presents different configurations, thereby causing the change of the end stiffness of the robot.

[0184] When the end position of the robot is determined, the overall stiffness of the end and the overall stiffness of the end face of the robot rotating around the z-axis of the world coordinate system {W} both present dynamic changes. These two indicators as evaluation criteria of the end stiffness have important guiding significance for the planning of the milling pose. However, when the rotation angle θ = 180°, although the overall stiffness of the end face is good, the overall stiffness of the end is in the worst state, therefore, the size of the overall stiffness of the end face cannot be simply judged according to the overall stiffness of the end. In terms of the above two indicators, it has obvious limitations to judge the comprehensive performance of the end stiffness by only one indicator. Therefore, when evaluating the pros and cons of the end stiffness based on the pose dependence, multiple evaluation indicators need to be considered comprehensively. This means that the optimization of the milling pose is essentially a multi-objective optimization problem, which needs to balance the relationship between various stiffness and non-stiffness indicators to achieve the optimization of the 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, transferring 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] where the diagonal elements are the direct stiffness, representing the stiffness along a certain axis. The off-diagonal elements are the coupling stiffness, representing the displacement in one direction will cause a force related to it in another direction. The end stiffness matrix K end is singular value decomposed: K end = Q AQ T ; where: A = diag [λ1, λ2, λ3].

[0197] A is a diagonal matrix containing eigenvalues (λ1> λ2> λ3> 0). The eigenvalues λ i represent the magnitude of stiffness in the principal axis direction. The largest eigenvalue corresponds to the direction with the largest stiffness, also known as the principal stiffness direction, and the smallest eigenvalue corresponds to the direction with the weakest stiffness.

[0198] Secondly, the orthogonal matrix Q: Q = [v1, v2, v3]. Q is an orthogonal matrix composed of three-dimensional column vectors v i , which are the unit eigenvectors corresponding to the eigenvalues λ i . i

[0199] The end stiffness ellipsoid is a core tool in robotics to describe the stiffness distribution of a mechanical system in task space. Its geometric form and mathematical properties directly reflect the resistance of the robot to external loads and anisotropy. The parameters of the end 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 determines that its geometric representation must be a quadric surface, and the positive definiteness further limits the possible closed surface to an ellipsoid. The linear transformation of the Jacobian matrix acts geometrically on the composite stiffness matrix K t , performing rotation and scaling. This operation will inevitably map a sphere or ellipsoid to a new ellipsoid. Therefore, there is a necessary correlation between the indicators related to the end stiffness of the robot and the stiffness ellipsoid.

[0203] The end stiffness ellipsoid volume evaluation index is an important parameter for measuring the overall stiffness performance of the robot end effector in the task space. The end stiffness ellipsoid volume and the stiffness anisotropy ellipsoid area, set the end stiffness ellipsoid volume expression:

[0204]

[0205] ​The ellipsoid volume is negatively correlated with the overall stiffness of the robot end, and the smaller the ellipsoid volume, the higher the overall stiffness; otherwise, the larger the ellipsoid volume, the lower the overall stiffness. If the ellipsoid approaches a sphere, the stiffness distribution is uniform; in addition, the shape characteristics of the ellipsoid can intuitively reflect the anisotropy of the stiffness distribution: when the ellipsoid approaches a sphere, the stiffness distribution is uniform in all directions; if the ellipsoid is significantly elongated, it indicates that there is a directional difference in stiffness performance; and the ellipsoid shape significantly extended in a specific direction indicates that the axial stiffness is relatively weak.

[0206] According to the ellipsoid volume calculation formula, the overall stiffness of the robot end on the height plane z = 1000 of the robot end in the same pose is calculated. In addition, the overall stiffness of the robot end is calculated when the robot end position is (2000, 0, 1000) and the end rotates around the z-axis of the world coordinate system {w} for one revolution. The standard external load is set as F std = [-15, 10, -96, 0, 0, 0] T According to the above parameters, the end stiffness ellipsoid volume calculation results of the two working conditions are as shown in Figure 1

[0207] According to the end stiffness ellipsoid volume thermal map analysis, it can be seen from Figure 1 (a) that as the robot end gradually moves away from the base, the volume of the end stiffness ellipsoid gradually decreases, and the stiffness of the robot end increases. In the upper right and lower right corners of figure (a), the color of the thermal map becomes lighter, and the volume of the end stiffness ellipsoid becomes larger, indicating that when the end reaches a certain distance, the stiffness of the robot end begins to decrease. Considering the motion performance of the robot in this pose, the dexterity of the robot gradually decreases when the end moves away from the base, so when planning the working area for the cutting task, it is not simply close to the base or move away from the base, but should comprehensively consider multiple indicators and reasonably select the working area for the cutting task.

[0208] It can be seen from Figure 1 (b) that when the end rotates around the z-axis of the world coordinate system {w} by an angle θ = 0° (360°), the volume of the robot end stiffness ellipsoid is the smallest, and the stiffness of the robot end is the largest; when θ = 180°, the volume of the robot end stiffness ellipsoid is the largest, and the stiffness of the robot end is the smallest. During the rotation of 0-180° and 180-360°, the volume of the robot end stiffness ellipsoid changes continuously and symmetrically. Obviously, when the robot end is in this position, the best overall stiffness performance can be obtained by selecting θ = 0° (360°) milling pose. However, in practical applications, further consideration should be given to the change in stiffness of the robot end stiffness ellipsoid when the end is used as a force surface under external load. This can be measured by the following indicators.

[0209] The volume of the end stiffness ellipsoid V F ​about-15, in order to match the target function and the consistency of the fitness function in the numerical order of magnitude, to avoid the optimization deviation caused by the magnitude difference, the original data is processed as follows: idx(V F ) = log V F + 17.

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

[0211] The stiffness anisotropy ellipse area is a geometric index for evaluating the compliance distribution of the end stiffness ellipsoid in a specific plane. When an external load acts on the robot end face, the robot end face acts as a force surface, and the stiffness ellipsoid is truncated by the plane, and the cross section forms an elliptical surface.

[0212] S321: When the robot end pose is R in the world coordinate system {W}, the three vector directions on the end face are:

[0213]

[0214] The end stiffness ellipsoid coordinate system is set as {E}, the origin coincides with the end coordinate system {F}, and the homogeneous transformation matrix in the world coordinate system {W} is:

[0215]

[0216] The pose of the end stiffness ellipsoid in the world coordinate system {W} is:

[0217]

[0218] The ellipsoid surface equation in the world coordinate system {W} is established as:

[0219]

[0220] S322: In the end stiffness ellipsoid coordinate system {E}, the pose of the robot end coordinate system {F} is calculated as:

[0221]

[0222] The expressions of the x vector and the y vector of the end coordinate system in the ellipsoid coordinate system {E} are:

[0223]

[0224] If the stiffness ellipsoid is to be truncated by the end face to obtain the elliptical surface, the intersection points of the x vector and the y vector of the end coordinate system with the ellipsoid need to be calculated. The length of the line segment from the origin to the intersection point represents the compliance coefficient and These coefficients reflect the ability of the robot's end to resist deformation in the corresponding direction.

[0225] S323: Obtain the elliptical surface of the stiffness ellipsoid cut by the end face based on the intersection points of the x vector and y vector of the end coordinate system and the ellipsoid.

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

[0227] Let the intersection points of the x vector and y vector and the ellipsoid be and Then:

[0228]

[0229] Solve the following equations:

[0230]

[0231] The solution is:

[0232]

[0233] According to the compliance coefficients on the end face of the robot, calculate the stiffness anisotropic ellipse area:

[0234]

[0235] In the analysis of anisotropic ellipse, the size of the ellipse area directly reflects the stiffness distribution characteristics of the robot end plane in different directions. The smaller the ellipse area, the better the overall stiffness performance of the robot end plane in this pose. Conversely, the larger the ellipse area, the worse the overall stiffness performance of the robot end plane in this pose. If and the numerical values are close, it indicates that the stiffness distribution of the robot on the end plane is relatively uniform, indicating that the anti-deformation ability of the robot in the x and y directions tends to be consistent. The robot will not be significantly deformed due to force in a particular direction. If and the numerical values differ significantly, it indicates that the stiffness of the robot on the end plane exhibits obvious anisotropy, and the anti-deformation ability in different directions differs greatly. In the direction with lower stiffness, the robot end may be more susceptible to external force interference and may have larger displacement or vibration.

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

[0237] According to the end stiffness ellipse area thermal map analysis, it is known from Figure 2 (a) that as the robot end gradually moves away from the base, the area of the end stiffness ellipse is gradually reduced, and the stiffness on the end plane increases, which has a similar trend to the stiffness ellipsoid volume. From Figure 2 (b) it is known that when θ = 90° and θ = 270°, the robot end stiffness ellipse area is maximum, indicating that the overall stiffness performance on the end plane is weakest at this time. In the range of 0-90°, 90-270° and 270-360°, the area of the end stiffness ellipse changes relatively gently, and there is no significant fluctuation, but if you want to further analyze the stiffness performance of the end plane in these three angle intervals, you also need to combine the numerical value of the flexibility coefficient to further judge whether the end plane ellipse has significant anisotropy. In addition, when θ = 180°, the volume of the end ellipsoid is minimum, but the area of the ellipse on the end plane is maximum, which shows that the volume of the end ellipsoid cannot be simply used to judge the size of the area of the ellipse on the end plane. According to the results given in Figure (b), in the planning of cutting tasks, it is necessary to avoid selecting poses with θ = 90° and θ = 270° to ensure that the robot end has sufficient stiffness and stability during machining.

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

[0239] In order to avoid the possible robot singularity in the task planning process and effectively evaluate the multi-degree-of-freedom motion ability of the robot in complex tasks, the motion dexterity index D m (Manipulability Measure, dexterity) is widely used to quantitatively analyze the motion flexibility and mobility of the robot. Its concept is related to the singular value of the Jacobian matrix, and the robot dexterity D m is defined as the reciprocal of the condition number k m :

[0240]

[0241] where k m is the condition number k(J(θ)) of the Jacobian matrix k m = k(J(θ));

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

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

[0244] where U∈R m×m , V∈R n×n , U and V are 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] Coordinates 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 robot milling operation, in addition to the above stiffness and kinematics performance related indicators, the smoothness of the robot in the milling process needs to be paid special attention. The smoothness of the robot when running mainly reflects two aspects: the motion ability at the current milling pose and the adjustment range between adjacent milling poses. First of all, it is necessary to ensure that the robot has sufficient motion ability at the current pose, and the lack of motion ability is manifested as the joint speed cannot meet the running speed of the end. In the case of constant end speed, insufficient motion ability will lead to joint speed over-limit, and further lead to overload and damage of the robot joint motor. Secondly, the adjustment range between adjacent milling poses must be controlled within a reasonable range. Too large adjustment range will lead to a large change in joint speed in a short time, and since the connecting rod of the industrial robot usually has a large mass and inertia, sharp joint movement will cause significant mechanical vibration, which in turn leads to instability in the robot movement, and ultimately leads to damage to the tool or workpiece. Therefore, sufficient motion ability is reflected as a constraint on the joint speed of the robot, and the adjustment range between adjacent milling poses is reflected as a constraint on the joint acceleration of the robot. While optimizing the above indicators, the motion smoothness of the robot must be ensured.

[0260] If the milling speed is v, the end speed of the robot is:

[0261] The joint speed can be obtained by the inverse Jacobian matrix:

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

[0263]

[0264] In the milling task, the robot reaches point B from point A. In the case of dense enough planning points, the distance between adjacent milling points and is:

[0265] If the feed speed is v f , the time t between A and B is:

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

[0267]

[0268] The average acceleration between adjacent poses needs to be less than the rated acceleration of the robot when working:

[0269]

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

[0271]

[0272] where R i is the pose of the robot end at the milling position point when the milling spindle rotates around the redundant axis, and T is the homogeneous transformation matrix of the robot end at this time. The fitness value of the multi-objective optimization function is the linear combination of the angle θ Fv1 between the external load and the principal stiffness, the index idx(V F ) of the end stiffness ellipsoid volume, the index idx(S F ) of the end ellipsoid area, and the dexterity index idx(D m ). The value of σ i is selected according to the influence degree of the objectives on the optimization results. The optimization objective is to minimize the fitness function fit(T i ), and the pose with the minimum fitness value is the target pose T target . To ensure the smoothness of the robot in the milling task, there are certain constraints on the joint speed and average joint acceleration. The joint speed needs to satisfy The joint acceleration needs to satisfy

[0273] MO-MLSWOAR is based on MCMLWOA, modifies the original multi-leader update mechanism and integrates the search strategy of PSO algorithm. Further, the original speed update strategy of PSO algorithm is improved to adaptive speed update strategy. The modified multi-leader update mechanism and adaptive speed update strategy are as follows:

[0274] S41: The core of the multi-leader update mechanism is the dynamic adjustment strategy of the number of leaders based on the change of the Pareto front solution set and the density of the fitness distribution. Set the improved multi-leader update mechanism: identify the current Pareto front solution set PF through non-dominated sorting (Non-dominated Sorting), and then introduce a proportion factor α to design the dynamic selection rule of 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, which is a number in the range of (0, 1), and its relationship with the fitness distribution density Cr of the current iteration number t is: Cr = 1-α;

[0277] The fitness distribution density Cr is used to evaluate the degree of aggregation of the fitness values of the individuals in the tth generation, and the calculation formula is:

[0278]

[0279] Where n is the population size, fit(X i ) is the fitness value of the ith individual in the population at the tth iteration, fit max and fit min are the maximum and minimum values of the fitness of the current population, respectively. If Cr is small, it indicates that the fitness distribution of the individuals is dense, indicating that the individuals may be trapped in a local optimal solution, and more leaders need to be introduced by increasing a to enhance search diversity; if Cr is large, the fitness distribution is scattered, and the number of leaders can be appropriately reduced to reduce the computational resources consumed in calculating the direction and weight of the leader.

[0280] S42: Set the speed adaptive update strategy: In the PSO algorithm, the speed update formula is the core mechanism for adjusting the motion direction and speed of the particle. The speed adaptive update strategy takes the basic PSO algorithm as a reference, and optimizes the original speed update strategy by introducing a dynamic inertia weight and an adaptive acceleration constant, and the optimized speed adaptive update formula is:

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

[0282] Where v i (t+1) is the speed of the ith particle in the t+1th iteration, w is the inertia weight, which controls the influence of the current speed of the particle, c1 and c2 are acceleration constants, which control the influence of individual cognition and social learning of the particle, respectively, r1 and r2 are random numbers in the range of [0, 1], which are used to increase the randomness of the search, p * (t) is the historical optimal position of the ith particle, g * (t) is the global optimal position, and x i (t) is the position of the ith particle in the tth iteration.

[0283] S43: Introduce dynamic inertia weight and adaptive acceleration constant to further optimize the performance of the algorithm:

[0284] The inertia weight w, the adaptive acceleration constants c1 and c2 are dynamically adjusted by linearly decreasing to balance the global search and local search capabilities:

[0285]

[0286] where w max and w min are the maximum and minimum of the inertia weight, T is the maximum number of iterations, t is the current iteration number, c 1，max and c 1,min are the maximum and minimum of the acceleration constant c1, c 2，max and c 2，min are the maximum and minimum of the acceleration constant c2.

[0287] In the MO-MLSWOAR algorithm, the overall process is consistent with the MCMLWOA, but the position updating mode is different from it. The position updating mode of the MCMLWOA follows the three basic strategies of the traditional WOA algorithm: spiral position updating, surrounding prey and random search. The MO-MLSWOAR combines the search strategy of the PSO with the position updating mode of the MCMLWOA, forming a hybrid position updating mechanism. The position updating mode of the traditional PSO algorithm is as follows:

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

[0289] S44: The search strategy of the traditional PSO is combined with the position updating mode of the MCMLWOA to form a hybrid position updating mechanism: the random search of the MO-MLSWOAR is replaced by the position updating mode of the PSO algorithm:

[0290]

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

[0292] In order to verify the performance of the MO-MLSWOAR algorithm, the present application will list the multiple indicators of the PSO algorithm, GA

[69] algorithm, basic WOA algorithm and MO-MLSWOAR algorithm respectively, including GD (Generational Distance, intergenerational distance), SP (Spacing, spacing), optimal individual and optimal fitness. GD is used to measure the average distance between the found solution set and the true Pareto frontier:

[0293]

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

[0295]

[0296] where d i is the distance of the ith solution in the solution set to its nearest neighbor solution, and d i is the average of these distances. The smaller the SP value, the more uniform the distribution of the solution set, reflecting the performance of the algorithm in terms of solution set diversity. Using the above algorithm, the given design variables and multi-objective optimization conditions are used to optimize the redundant axis rotation angle of the same milling position, and the results are shown in FIG. 1: Figure 5

[0297] It can be seen that the MO-MLSWOAR algorithm performs best in terms of comprehensive performance, and its convergence and uniformity of solution set distribution are significantly better than those of other algorithms. PSO and WOA perform similarly in terms of convergence and uniformity, but are slightly inferior to MO-MLSWOAR, while GA has the highest GD and SP values, indicating that it has a slower convergence speed and a less uniform solution set distribution. In terms of optimal solution quality, MO-MLSWOAR, PSO and WOA all obtain the same optimal fitness value, but the optimal position of MO-MLSWOAR is more accurate than that of PSO and WOA; in contrast, the optimal position of GA deviates greatly, and the fitness value is also slightly worse. In summary, MO-MLSWOAR exhibits higher solution accuracy and stability in multi-objective optimization problems.

[0298] To verify the effectiveness of the multi-objective pose optimization algorithm, the present application designs and conducts multiple milling experiments. After the experiment is completed, the machined surface is sampled by wire cutting. Then, the Alicona InfiniteFocus G5 three-dimensional surface measuring instrument is used to detect the milling surface. Through the reconstruction of the topography of the milling surface, the topographic features and surface roughness of the milling trace can be accurately obtained. Table 3 shows the experimental conditions and parameters:

[0299] Table 3 Experimental conditions and parameters

[0300]

[0301] According to the three parameters of the experiment, the MO-MLSWOAR multi-objective optimization algorithm is used to optimize the milling pose of a straight-line milling path with a length of 300 mm, with a step size of 1 mm. With 50 mm as the node, the z-axis rotation angle and joint angle data are recorded as shown in the following table:

[0302] Table 4 Rotation angle and joint angle (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] Referring to the milling surface quality analysis diagram as shown in Figure 6

[0305] Overall, all parameters of Experiment One are significantly higher than those of Experiment Two and Experiment Three. Compared with Experiment One and Experiment Two, the reduction of Ra, Rq and Rt is 76.4%, 75.4% and 65.2% respectively. Compared with Experiment One and Experiment Three, the reduction of Ra, Rq and Rt is 78.1%, 77.2% and 69.1% respectively. This significant difference shows that the installation method of cross-axis is superior to that of parallel-axis in milling processing. Further comparison between Experiment Two and Experiment Three shows that Experiment Three is better in many key parameters, with the reduction of Ra being 7.1%, the reduction of Rz being 10.9% and the reduction of Rsm being 13.7%. Except that Rp is slightly higher than that of Experiment Two, the rest of the parameters are smaller than those of Experiment Two. The results show that the optimization of milling pose can further reduce the surface roughness and the milling trace distance, thereby comprehensively improving the processing quality.

[0306] The fitness value and the dexterity and stiffness performance indicators varying with displacement under each experimental condition are calculated by the pose optimization mathematical model, as shown in Figure 7 and Figure 8

[0307] In terms of fitness, the fitness value of Experiment Three is superior to that of Experiment One and Experiment Two, and always shows a downward trend. At the beginning of milling, the fitness value of Experiment Two is superior to that of Experiment One, but near the end of milling, the fitness value of Experiment Two is higher than that of Experiment One. In theory, the stiffness and motion performance of the robot end under the condition of Experiment One should be superior to that under the condition of Experiment Two near the end of milling. However, according to the experimental results, the milling surface quality of Experiment Two is inferior to that of Experiment One, which is mainly due to the difference in the installation method of the milling spindle. In terms of dexterity and stiffness performance indicators, the dexterity of the robot under each experimental condition does not differ much. In terms of stiffness performance indicators, the stiffness on the end face under the condition of Experiment Three is the best. Since the load of the robot end is directly applied to the end face, although the overall stiffness under the condition of Experiment Three is not as good as that under the condition of Experiment Two, the milling surface quality under the condition of Experiment Three is the best. In summary, through overall fitness analysis, it can be concluded that the milling quality under the installation method of cross-axis is significantly superior to that under the installation method of parallel-axis. Through comparison of dexterity and stiffness performance indicators, it can be seen that the overall stiffness performance of the milling end face is crucial to the milling surface quality.

[0308] ​​To sum up, the application firstly analyzes the pose optimization target of the six-degree-of-freedom joint robot under the typical five-degree-of-freedom milling operation, i.e. the rotation angle around the milling redundant axis. Secondly, multiple performance indicators and multi-objective optimization functions are designed. By comprehensively considering the motion performance indicators, stiffness performance indicators and robot joint motion limitation conditions, the angle optimization is carried out in the range of-180°-180° through the MO-MLSWOAR multi-objective optimization algorithm, and the best milling posture of the milling point is obtained. The performance of MO-MLSWOAR, PSO, ordinary WOA and GA is compared through four indicators of GD, SP, optimal fitness and optimal individual position. The experiment proves that MO-MLSWOAR shows higher solution accuracy and stability in the multi-objective optimization problem. After the pose optimization, three groups of experiments are designed to compare the milling effect of the optimized and unoptimized robots. The optimized pose is used for milling, and compared with the unoptimized pose, the surface roughness, average peak-to-valley distance and waviness of the optimized pose milling are significantly improved. The experiment proves that on the basis of the cross-axis installation, the optimization algorithm has a significant improvement effect on the robot milling quality.

Claims

1. A method of milling pose optimization based on milling redundant axes, characterized in that, Comprise the following steps: S1: the milling spindle is vertically installed on the sixth axis of the robot, so that the robot determines the pose of milling with six degrees of freedom, realizes the redundancy of the degree of freedom of the milling cutter pose kinematics, and describes the rotation of the z-axis of the milling cutter coordinate system with redundant degrees of freedom, that is, the deflection of the milling spindle; S2: solving the milling redundant axis direction vector based on the PCA method; S3: setting the pose optimization target of the six-degree-of-freedom joint robot under the typical five-degree-of-freedom milling operation; S4: through the MO-MLSWOAR multi-objective optimization algorithm, the angle is optimized in the range of-180°~180°, and the best milling posture of the corresponding milling point is obtained.

2. The method of claim 1, wherein, The specific process of describing the rotation of the milling cutter coordinate system z-axis in step S1, that is, the deflection of the milling spindle, is as follows: S11: in the milling task, the pose of the milling cutter coordinate system {M} in the base coordinate system {B} is described as a homogeneous transformation matrix: wherein, Each column in the matrix represents the projection of the coordinate axis of the milling tool coordinate system {M} in the base coordinate system {B}, the rotation axis of the milling tool The unit vector in the base coordinate system {B} is: S12: Based on the degree of freedom redundancy of the milling cutter pose, the milling cutter coordinate system {M} rotates around the rotation axis of the milling cutter Rotation, if the rotation angle is γ, according to the Rodrigues rotation formula, the rotation matrix is expressed as: Wherein, I is the unit matrix, K is the skew-symmetric matrix of the rotation axis, and γ is the rotation angle: S13: New pose of the rotary milling cutter coordinate system {M} is: S14: based on the relative pose between the coordinate systems and the kinematics equation of the robot, the relationship between the joint angle of the robot and the rotation angle of the milling cutter coordinate system {M} is established: S15: when the robot performs a typical milling task, 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 the world coordinate system {W}: wherein is the milling position, a m is the direction of the milling spindle, γ is the angle of rotation of the milling coordinate system around the milling spindle.

3. The method of claim 1, wherein, In step S2, the milling redundant axis direction vector is solved based on the PCA method. The normal vector of a point p(x, y, z) on the point cloud model obtained from a certain free-form surface is solved, and the specific process is as follows: S21: select the field: The neighborhood points p1, p2,..., p of p are selected by the Euclidean distance n , which satisfy: Wherein, r is the radius of the spherical neighborhood; Calculate the mean vector of the neighborhood points: Take the set of neighborhood points as the sample D: S22: construct the covariance matrix: Center the sample D: The covariance matrix ∑ is: wherein is the variance of x; is the covariance of x, y, and so on, to solve the covariance matrix∑; S23: solve the eigenvalue and eigenvalue vector: The eigenvalue and eigenvalue vector of the covariance matrix ∑ satisfy: ∑v i = λ i v i ; eigenvalue λ i is a solution of the eigen equation of the covariance matrix ∑, v i is an eigenvector: det(∑-λ i I) = 0; S24: for a 3×3 covariance matrix, there are three eigenvalues λ1, λ2, λ3, and the size of the eigenvalue represents the variance of the data in the corresponding eigenvector direction. For normal vector calculation, the eigenvector corresponding to the smallest eigenvalue reflects the direction in which the point cloud data changes the least in the local area, so the eigenvector with the smallest eigenvalue is selected as the normal vector: λ n = min(λ1, λ2, λ3); Solve for the eigenvector v by the following equation n : (∑ - λ n I) v n = 0, Feature vector v n That is, the normal vector at point p.

4. The method of claim 1, wherein, In step S3, the specific process of setting the pose optimization target of the six-degree-of-freedom joint robot under the typical five-degree-of-freedom milling operation is as follows: S31: represent the angle between the principal stiffness direction and the combined force of the external load: The external load at the end of the robot includes a milling force F m and an original load F o , including a flange, an electric spindle, and a force sensor, wherein the milling force F m is acted on in a world coordinate system, the original load F o is acted on in an end face, and the external load F ext is expressed as: External load F ext The angle θ of the main rigidity direction The angle θ of the main rigidity direction Fv1 is: S32: represent the stiffness performance index: the volume of the end stiffness ellipsoid and the area of the stiffness anisotropic ellipse, and set the expression of the volume of the end stiffness ellipsoid: S33: represent the motion performance index: the dexterity of the robot: The robot dexterity D m is defined as the inverse of the condition number k m : where k m is the condition number of the Jacobian matrix k(J(θ)): k m = k(J(θ)); Based on the singular value analysis theory of matrix, the Jacobian matrix J(θ) of the robot in any shape and position is analyzed through singular value decomposition: J(θ)=U∑V; where U ∈ R m×m , V ∈ R n×n , U and V are orthogonal matrices, and ∑ is: where (σ1, σ2, …, σ m ) are singular values of J(θ), σ1≥σ2≥…≥σ m ≥0, σ1 is the largest singular value, σ m is the smallest singular value, and the condition number k(J) of the Jacobian matrix has the following relationship with the singular values: Use MATLAB software to simulate the motion process of the robot end along the x direction and y direction of the world coordinate system {W}, and obtain the dexterity in the motion process; S34: The motion capability in the current milling pose and the adjustment range between adjacent milling poses represent the robot motion smoothness.

5. The method of claim 4, wherein, The specific process of representing the anisotropic stiffness ellipse area in step S32 is as follows: S321: In the world coordinate system {W}, when the robot end pose is R, the three vector directions on the end face are: Set the end stiffness ellipsoid coordinate system as {E}, the origin of which coincides with the end coordinate system {F}, and the homogeneous transformation matrix in the world coordinate system {W} is: The posture of the end stiffness ellipsoid in the world coordinate system {W} is: In the world coordinate system {W}, the ellipsoid surface equation is established as: S322: In the end stiffness ellipsoid coordinate system {E}, the pose of the robot end coordinate system {F} is calculated as: The expressions of the x vector and y vector of the end coordinate system in the ellipsoid coordinate system {E} are: S323: Based on the intersection points obtained by calculating the intersection of the x vector and y vector of the end coordinate system with the ellipsoid, the elliptical surface of the stiffness ellipsoid cut by the end face is obtained.

6. The method of claim 5, wherein, The specific process of obtaining the elliptical surface of the stiffness ellipsoid cut by the end face in step S323 is as follows: Let the intersection of the x vector and the y vector with the ellipsoid be denoted as and then: Solve the following equations: The solution is: According to the flexibility coefficient on the robot end face, the stiffness anisotropic ellipse area is calculated as:

7. The method of claim 1, wherein, The specific process of step S4 is as follows: S41: Set the improved multi-leader update mechanism: identify the current Pareto front solution set PF through non-dominated sorting, and then introduce a proportion factor α to design a dynamic selection rule for the number of leaders k: N l = min(a card(PF), 50); S42: Set the speed adaptive update strategy: introduce dynamic inertia weight and adaptive acceleration constant to optimize the original speed update strategy, and the optimized speed adaptive update formula is: v i (t + 1) = w - v(t) + c1 - r1 - (p * (t) - x i (t)) + c2 - r2 - (g * (t) - x i (t)); where v i (t + 1) is the velocity of the i-th particle in the t+1 iteration, w is the inertia weight that controls the influence of the current velocity of the particle, c1 and c2 are acceleration constants that control the influence of the individual cognition and social learning of the particle, respectively, r1 and r2 are random numbers in the range [0, 1] that are used to increase the randomness of the search, p * (t) is the historical best position of the i-th particle, g * (t) is the global best position, x i (t) is the position of the i-th particle in the t iteration. S43: Introduce dynamic inertia weight and adaptive acceleration constant to further optimize the algorithm performance; S44: Follow the traditional WOA algorithm's spiral position update, surround prey and random search, and combine the search strategy of the traditional PSO with the position update mode of MCMLWOA to form a hybrid position update mechanism: replace the random search of MO-MLSWOAR with the position update mode of PSO algorithm: S45: Based on the MO-MLSWOAR optimization algorithm, the optimal milling pose of the corresponding milling point is obtained.

8. The method of claim 7, wherein, The specific process of introducing dynamic inertia weight and adaptive acceleration constant to further optimize the algorithm performance in step S43 is as follows: The inertia weight w, adaptive acceleration constant c1 and c2 are dynamically adjusted in a linear decreasing manner to balance the global search and local search capabilities: where w max and w min are the maximum and minimum inertia weight respectively, T is the maximum number of iterations, t is the current iteration number, c 1，max and c 1,min are the maximum and minimum acceleration constant c1 respectively, c 2，max and c 2，min are the maximum and minimum acceleration constant c2 respectively.

Images