Multi-axis cooperative automatic control method and system for flange welding seam milling

Abstract

The invention provides a multi-axis cooperative automatic control method and system for flange plate welding seam milling, and relates to the technical field of automatic control, the method comprises the following steps: step 1, collecting three-dimensional contour data of a flange plate circular seam area, and detecting the pose deviation of a mobile platform in real time to obtain base error data; 2, according to the three-dimensional contour data, geometric feature points of the circular seam are extracted, two feature lines are obtained, and the included angle between the two feature lines is calculated to define a local machining area of the circular seam. According to the method, the flange plate circular seam three-dimensional contour data are collected, the pose deviation is detected, the annular closed control path is constructed, the load vibration characteristics are recognized in combination with milling force feedback and spindle motor current signals, the track is optimized, multi-axis cooperative automatic control over flange plate welding seam milling is achieved, and the machining precision, stability and efficiency are effectively improved.

Description

Technical Field

[0001] This invention relates to the field of automatic control technology, and in particular to a multi-axis collaborative automatic control method and system for milling flange weld seams. Background Technology

[0002] For milling flange circumferential seams, using industrial robots or special-purpose machines in conjunction with mobile platforms is a common solution. However, in practical applications, this method sometimes faces some technical challenges. On the one hand, due to the large size of the flange workpiece or the possibility of installation deviations, the actual position of its circumferential seam is prone to a certain degree of coordinate system deviation from the theoretical programming model. For example, in the machining of flange butt joint circumferential seams in large pressure vessel cylinders, the end face runout of the flange itself or the installation tilt of the platform may cause the entire circumferential seam trajectory to deviate to a certain extent relative to the robot base. Existing offline programming or teaching methods have limited dynamic compensation capabilities for such deviations, which may lead to a certain deviation between the milling path and the weld area, thereby affecting the machining accuracy.

[0003] On the other hand, during milling, the unevenness of the weld allowance can cause fluctuations in the milling load, which in turn leads to vibrations in the tool and robot. For example, when milling a weld with localized variations in width and height, changes in the depth of cut into the workpiece can cause fluctuations in the cutting force. Existing control strategies may not always respond smoothly enough to such load variations, which could affect the consistency of the machined surface quality. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a multi-axis collaborative automatic control method and system for milling flange weld seams, so as to avoid the problem of large-scale idle tool movement or milling beyond the weld seam area during processing, shorten the ineffective processing time, and improve the work efficiency.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: A first aspect is a multi-axis collaborative automatic control method for milling flange weld seams, the method comprising: Step 1: Collect the three-dimensional contour data of the flange circumferential seam area and detect the pose deviation of the moving platform in real time to obtain the base error data; Step 2: Based on the 3D contour data, extract the geometric feature points of the circumferential seam to obtain two feature lines, and calculate the included angle between the two feature lines to define the local processing area of ​​the circumferential seam. Step 3: Set up first and second control points inside and outside the local processing area, and construct a closed loop control path based on the time series data of the first and second control points; Step 4: Based on the geometric characteristics of the circular closed control path, generate trajectory correction parameters to correct the transformation relationship between the robot's base coordinate system and the theoretical programming coordinate system, and generate a real-time compensation signal for the base coordinate system. Step 5: Based on the real-time compensation signal of the base coordinate system, adjust the joint motion parameters of the robot to obtain the pre-compensated motion trajectory, and collect the interaction force between the milling cutter and the workpiece to obtain milling force feedback data; Step 6: Based on the milling force feedback data and the current signal of the spindle motor, identify the load changes and vibration characteristics, determine the load vibration characteristic data, adjust the robot's end-effector pose, and generate force control signals. Step 7: Correct the pre-compensated motion trajectory using force control signals to obtain the optimized motion trajectory; based on the optimized motion trajectory, control the linear motion of the mobile platform and the joint motion of the robot through a multi-axis collaborative controller to complete the milling of the flange weld.

[0006] Secondly, a multi-axis collaborative automatic control system for milling flange weld seams includes: The acquisition module is used to acquire three-dimensional contour data of the flange circumferential seam area and detect the pose deviation of the moving platform in real time to obtain the base error data. The extraction module is used to extract the geometric feature points of the circumferential seam based on the three-dimensional contour data, obtain two feature lines, and calculate the included angle between the two feature lines to define the local processing area of ​​the circumferential seam. The module is used to set first and second control points inside and outside the local processing area, and to construct a closed loop control path based on the time series data of the first and second control points. The compensation module is used to generate trajectory correction parameters based on the geometric characteristics of the circular closed control path, so as to correct the transformation relationship between the robot's base coordinate system and the theoretical programming coordinate system, and generate a real-time compensation signal for the base coordinate system. The feedback module is used to adjust the robot's joint motion parameters based on the real-time compensation signal of the base coordinate system, obtain the pre-compensated motion trajectory, and collect the force between the milling cutter and the workpiece to obtain milling force feedback data. The identification module is used to identify load changes and vibration characteristics based on milling force feedback data and spindle motor current signals, determine load vibration characteristic data, adjust the robot's end-effector pose, and generate force control signals. The optimization module is used to correct the pre-compensated motion trajectory using force control signals to obtain the optimized motion trajectory. Based on the optimized motion trajectory, the linear motion of the mobile platform and the joint motion of the robot are controlled by a multi-axis collaborative controller to complete the milling of the flange weld.

[0007] Thirdly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0008] The above-described solution of the present invention has at least the following beneficial effects: Laser scanning is used to acquire the 3D contour of the circumferential weld and detect the pose of the moving platform. Based on the geometric characteristics of the circular path, the transformation relationship between the robot's base coordinate system and the theoretical coordinate system is corrected to dynamically compensate for the base error. This avoids the shortcomings of offline programming / teaching, which cannot dynamically adjust parameters, ensuring precise alignment between the milling path and the weld and reducing the risk of trajectory misalignment. A six-dimensional force sensor is used to acquire milling force and analyze the spindle motor current signal to identify load changes and vibration characteristics. Force control signals are then generated to fine-tune the robot's end effector pose in real time, achieving a smooth response to load fluctuations, reducing the impact of vibration on the machined surface, and ensuring the smooth milling of weld seams in different sections. Consistency in quality: By extracting geometric feature points of the circumferential weld, fitting the baseline, and defining the local processing area, the system focuses only on the core section of the weld. At the same time, a closed loop path is constructed based on dual control points to avoid large-scale idle tool movement or milling beyond the weld area during processing, thus shortening ineffective processing time and improving work efficiency. The linear motion of the mobile platform and the joint motion of the robot are uniformly scheduled by a multi-axis collaborative controller to avoid the accumulation of deviations in single-axis motion. Combined with dual parameter tuning logic, the system has stronger adaptability to workpiece deviations and load changes, reducing the risk of processing interruption or rework and improving overall operational stability. Attached Figure Description

[0009] Figure 1 This is a flowchart illustrating a multi-axis collaborative automatic control method for milling flange weld seams, provided by an embodiment of the present invention.

[0010] Figure 2 This is a schematic diagram of a multi-axis collaborative automatic control system for milling flange weld seams, provided by an embodiment of the present invention. Detailed Implementation

[0011] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0012] like Figure 1 As shown in the figure, an embodiment of the present invention proposes a multi-axis collaborative automatic control method for milling flange weld seams, the method comprising the following steps: Step 1: Collect the three-dimensional contour data of the flange circumferential seam area and detect the pose deviation of the moving platform in real time to obtain the base error data; Step 2: Based on the 3D contour data, extract the geometric feature points of the circumferential seam to obtain two feature lines, and calculate the included angle between the two feature lines to define the local processing area of ​​the circumferential seam. Step 3: Set up first and second control points inside and outside the local processing area, and construct a closed loop control path based on the time series data of the first and second control points; Step 4: Based on the geometric characteristics of the circular closed control path, generate trajectory correction parameters to correct the transformation relationship between the robot's base coordinate system and the theoretical programming coordinate system, and generate a real-time compensation signal for the base coordinate system. Step 5: Based on the real-time compensation signal of the base coordinate system, adjust the joint motion parameters of the robot to obtain the pre-compensated motion trajectory, and collect the interaction force between the milling cutter and the workpiece to obtain milling force feedback data; Step 6: Based on the milling force feedback data and the current signal of the spindle motor, identify the load changes and vibration characteristics, determine the load vibration characteristic data, adjust the robot's end-effector pose, and generate force control signals. Step 7: Correct the pre-compensated motion trajectory using force control signals to obtain the optimized motion trajectory; based on the optimized motion trajectory, control the linear motion of the mobile platform and the joint motion of the robot through a multi-axis collaborative controller to complete the milling of the flange weld.

[0013] In this embodiment of the invention, laser scanning is used to acquire the three-dimensional contour of the circumferential weld and detect the pose of the moving platform. Based on the geometric characteristics of the circumferential path, the transformation relationship between the robot's base coordinate system and the theoretical coordinate system is corrected to dynamically compensate for the base error. This avoids the shortcomings of offline programming / teaching, which cannot dynamically adjust parameters, ensuring precise alignment between the milling path and the weld seam and reducing the risk of trajectory misalignment. A six-dimensional force sensor is used to acquire the milling force and analyze the spindle motor current signal to identify load changes and vibration characteristics. Force control signals are then generated to fine-tune the robot's end effector pose in real time, achieving a compliant response to load fluctuations, reducing the impact of vibration on the machined surface, and ensuring the smooth operation of different areas. The system ensures consistent milling quality of weld seams. By extracting geometric feature points of the circumferential weld, fitting baselines, and defining local processing areas, it focuses only on the core section of the weld seam. Simultaneously, it constructs a closed-loop path based on dual control points to avoid large-scale idle tool movement or milling beyond the weld seam area during processing, thus shortening ineffective processing time and improving work efficiency. Through a multi-axis collaborative controller, it uniformly schedules the linear motion of the mobile platform and the joint motion of the robot, avoiding the accumulation of deviations in single-axis motion. Combined with dual parameter tuning logic, the system has stronger adaptability to workpiece deviations and load changes, reducing the risk of processing interruptions or rework and improving overall operational stability.

[0014] In a preferred embodiment of the present invention, step 1 includes: Step 100: Acceleration and angular velocity data of the mobile platform are collected by a multi-degree-of-freedom inertial measurement unit (IMU) installed on the mobile platform. Real-time pose data of the mobile platform is obtained through attitude calculation based on the acceleration and angular velocity data. Specifically, the IMU is fixed to the base frame of the mobile platform using a custom-made aluminum alloy mounting base. The upper surface of the mounting base maintains a preset parallelism requirement with the guide reference surfaces of the X-axis and Y-axis linear guides of the mobile platform, with the parallelism error controlled within 0.02 mm / m. This ensures that the installation reference of the measurement unit is consistent with the motion reference of the mobile platform. The installation position is specifically located at the geometric center of the mobile platform in the horizontal plane, i.e., the line connecting the centers of symmetry of the two frames on both sides of the mobile platform in the X-axis direction and the line connecting the front and rear sides in the Y-axis direction. At the intersection of the lines connecting the centers of symmetry of the rear frame; in the vertical direction (Z-axis direction), the vertical distance between this installation position and the working surface of the moving platform is a preset 80mm, and the height of the mounting base is controlled within ±0.05mm to ensure the positional accuracy of the measuring unit in the Z-axis direction; the three sensitive axes of the multi-degree-of-freedom inertial measurement unit correspond one-to-one with the motion reference axes of the moving platform, that is, the X-axis sensitive axis is parallel to the central axis of the X-axis linear guide, the Y-axis sensitive axis is parallel to the central axis of the Y-axis linear guide, and the Z-axis sensitive axis is completely coincident with the central axis of the moving platform's rotation around the vertical direction (Z-axis). The coaxiality error between each sensitive axis and the corresponding motion reference axis does not exceed 0.1mm, ensuring that the measurement data can accurately reflect the motion state of the moving platform along each reference axis.

[0015] The mounting base is precisely fitted with the positioning holes of the moving platform base frame through two cylindrical locating pins with a diameter of 8 mm. The diameter tolerance of the locating pins is h6, and the tolerance of the positioning holes of the base frame is H7. The fitting clearance is controlled between 0.006 mm and 0.018 mm to ensure that the repeatability error of the installation position of the measuring unit is within 0.03 mm. At the same time, the mounting base is firmly connected to the base frame by four evenly distributed M6 hexagon socket head cap screws. The screws are tightened with a preset pre-tightening torque of 12 N·m to ensure a rigid connection between the mounting base and the base frame and avoid relative displacement during the movement of the moving platform or the vibration of milling. With a preset sampling frequency of 500 Hz to 1 kHz, the multi-degree-of-freedom inertial measurement unit synchronously collects the acceleration data (ax, ay, az) of the moving platform in the three-dimensional rectangular coordinate system and the angular velocity data (ωx, ωy, ωz) around the three-dimensional coordinate axes. For the acceleration data (ax, ay, az) collected by the multi-degree-of-freedom inertial measurement unit, where ax represents the acceleration in the X-axis direction, ay represents the acceleration in the Y-axis direction, and az represents the acceleration in the Z-axis direction, zero-drift compensation is first performed. That is, in the state where there is no external input to the measuring unit, the output values of each axis are collected as the zero-drift deviation values (the zero-drift deviation value of the X-axis is denoted as ax0, the Y-axis is denoted as ay0, and the Z-axis is denoted as az0). Subtract the zero-drift deviation value of the corresponding axis from the original acceleration data axis by axis. That is, the data after zero-drift compensation for the X-axis is , and for the Y-axis is , and for the Z-axis is , completing the zero-drift compensation. The system pre-stores a temperature-deviation calibration data table for the multi-degree-of-freedom inertial measurement unit. This table contains multiple calibration temperature points, such as -10°C, 0°C, 10°C, 20°C, 30°C, 40°C, and the temperature drift deviation values of each axis. First, obtain the current measurement ambient temperature T through the temperature sensor installed near the measuring unit. Find the two calibration temperature points adjacent to T in the temperature-deviation calibration data table, denoted as T1 (less than T) and T2 (greater than T), and the corresponding temperature drift deviation values of the X-axis are Δax1 and Δax2 respectively. Calculate the temperature drift compensation amount Δax_comp of the X-axis at the current temperature T through the linear interpolation algorithm. The calculation formula is . Subtract the temperature drift compensation amount Δax_comp from the data ax_comp after zero-drift compensation of the X-axis to obtain the corrected acceleration data of the X-axis , that is . Using the same method, calculate the temperature drift compensation amounts Δay_comp and Δaz_comp of the Y-axis and Z-axis at the current temperature respectively, and then through 、 , obtain the corrected acceleration data of the Y-axis and Z-axis .

[0016] For the angular velocity data (ωx, ωy, ωz, where ωx represents the angular velocity around the X-axis, ωy represents the angular velocity around the Y-axis, and ωz represents the angular velocity around the Z-axis) collected by the multi-degree-of-freedom inertial measurement unit, the same zero-drift compensation and temperature drift correction process as for the acceleration data is performed. First, the zero-drift deviation value of each axis is obtained (the zero-drift deviation value around the X-axis is denoted as ωx0, around the Y-axis as ωy0, and around the Z-axis as ωz0). The corresponding zero-drift deviation value is subtracted from the original angular velocity data axis by axis to obtain the zero-drift compensated angular velocity data. Then, based on the temperature-deviation calibration data table and the current ambient temperature, the temperature drift compensation for each axis is calculated using linear interpolation. The corrected angular velocity data is obtained by subtracting the corresponding temperature drift compensation amount from the zero-drift compensated angular velocity data axis by axis. (in =ωxcomplement - Δωxcomplement, =ωycomplement-Δωycomplement, =ωzcomplement-Δωzcomplement, Represents the corrected angular velocity about the X-axis. Represents the corrected angular velocity about the Y-axis. (Represents the corrected angular velocity about the Z-axis); based on corrected acceleration data. The trapezoidal integration method is used to perform two integration operations to obtain the real-time position data of the mobile platform. First, the sampling time interval Δt is determined. The value of Δt is equal to 1 divided by the sampling frequency of the multi-degree-of-freedom inertial measurement unit. For example, when the sampling frequency is 500Hz, The initial integration of the velocity components is performed with a step size of Δt, integrating the X-axis corrected acceleration data ax'. The calculation formula is as follows: Where vx(t) is the velocity component along the X-axis at the current moment, and vx(t-1) is the velocity component along the X-axis at the previous moment. This is the X-axis corrected acceleration data at the current moment. The acceleration data after X-axis correction at the previous moment; using the same formula, respectively... Integrating yields the Y-axis velocity component vy(t) and the Z-axis velocity component vz(t) (vx, vy, and vz are collectively referred to as real-time velocity components); the second integral calculates the position component with a step size of Δt, integrating the X-axis velocity component vx, with the following formula: Where x(t) is the position component in the X-axis direction at the current moment, and x(t-1) is the position component in the X-axis direction at the previous moment; using the same formula, integrate vy and vz respectively to obtain the position component y(t) in the Y-axis direction and the position component z(t) in the Z-axis direction. Integrate x(t), y(t), and z(t) to obtain the real-time position data of the mobile platform (x, y, and z represent the real-time position in the X-axis, Y-axis, and Z-axis directions, respectively).

[0017] Based on corrected angular velocity data The angular changes of the moving platform around each axis are calculated using trapezoidal integrals with a step size of Δt: angular change around the X-axis. The calculation is performed in discrete form. Where Δθx(t) is the change in angle around the X-axis at the current moment. This represents the change in angle around the X-axis at the previous moment; using the same method, calculate the change in angle around the Y-axis. and the change in angle around the Z-axis During the system initialization phase, the moving platform is calibrated using a high-precision level to obtain the initial attitude angles. This represents the attitude angle around the X-axis in the initial state. This represents the attitude angle around the Y-axis in the initial state. (This represents the initial attitude angle around the Z-axis); add the change in angle of each axis to the corresponding initial attitude angle to obtain the real-time attitude angle, i.e., the real-time attitude angle around the X-axis. Real-time attitude angle around the Y-axis Real-time attitude angle around the Z-axis ( Collectively referred to as real-time attitude angles); combining real-time position data (x, y, z) with real-time attitude angles ( The data is integrated according to a 6-DOF parameter format (including position parameters in the X, Y, and Z axes, as well as attitude angle parameters around the X, Y, and Z axes) to form the real-time pose data of the mobile platform.

[0018] Step 101: Scan the flange circumferential seam area with a laser scanner to obtain three-dimensional contour data of the seam surface. Compare the real-time pose data with the preset theoretical pose data to obtain the comparison result. Calculate the pose deviation of the moving platform based on the comparison result to generate base error data. Specifically, the laser scanner is fixed to the crossbeam at the front end of the moving platform via an L-shaped metal bracket. The bracket and the crossbeam of the moving platform are rigidly connected by bolts. The detection center of the laser scanner and the installation center of the multi-degree-of-freedom inertial measurement unit maintain a preset rigid positional relationship in three-dimensional space, i.e., a distance of 300mm in the X-axis direction, 50mm in the Y-axis direction, and 120mm in the Z-axis direction. This positional relationship is calibrated by a coordinate measuring machine. The data is pre-stored in the control system to ensure accurate correlation between their positional data. The laser scanner's scanning lens faces the flange circumferential seam area, with the lens axis maintaining a preset 45° angle with the flange end face. The scanning range covers the entire circumferential seam area and extends 50mm to the flange base material area on both sides of the circumferential seam to ensure complete acquisition of the surface contour data of the circumferential seam and surrounding base material. The circumferential seam area is 3D scanned at a preset scanning frequency of 100Hz to 200Hz to obtain 3D point cloud data containing the surface contour of the circumferential seam. The point cloud data is then denoised using a 5×5 window Gaussian filtering algorithm to filter out noise points caused by environmental interference and eliminate abnormal points that are more than 5mm away from the theoretical position of the circumferential seam, resulting in pre-processed point cloud data.

[0019] The transformation matrix between the laser scanner and the moving platform coordinate system is preset (this matrix contains 3 translation parameters). and 3 rotation parameters The preprocessed point cloud data is transformed from the laser scanner's own coordinate system to the moving platform's coordinate system to generate three-dimensional contour data of the circumferential seam surface; preset theoretical pose data of the moving platform is retrieved from the control system's parameter database, which includes theoretical position ( ) and theoretical attitude angle ( ); to obtain the real-time position from the real-time pose data ( ) and theoretical position ( Comparison along each axis, using formulas The translational deviations in the X, Y, and Z axes were calculated. ); will display the real-time attitude angle ( ) and theoretical attitude angle ( Comparison along each axis, using formulas The rotational deviations around the X, Y, and Z axes were calculated. ); the translational deviation and the rotational deviation are calculated according to ( The vector form of the data is integrated to generate basis error data.

[0020] In this embodiment, the pose of the mobile platform and the actual shape of the circumferential seam are captured through the collaborative sensing of the multi-degree-of-freedom inertial measurement unit and the laser scanner. The real-time pose data obtained after correction and integration can accurately reflect the actual motion state of the mobile platform. The three-dimensional contour data generated by the laser scan can describe the true geometric shape of the circumferential seam. By calculating the deviation between the real-time pose and the theoretical pose, the generated base error data effectively solves the problem of base coordinate system deviation caused by flange installation tilt, end face runout, etc., which helps to improve the overall accuracy of milling.

[0021] In a preferred embodiment of the present invention, step 2 includes: Step 200 involves preprocessing the 3D contour data to obtain optimized 3D contour data. Specifically, this includes: for each point in the 3D point cloud, selecting its 50 neighboring points as neighborhood points, and calculating the distance from the point to these 50 neighborhood points; calculating the average distance based on these distances, and then calculating the standard deviation; comparing the actual distance of a single point with the average distance from its neighbors; if the actual distance is greater than three times the standard deviation of the average distance, the point is identified as an outlier, such as an isolated point caused by environmental interference, and is removed from the point cloud; subsequently, a voxel mesh is performed. Grid sampling was performed, with the voxel size set to 0.5 mm × 0.5 mm × 0.5 mm. The statistically filtered point cloud space was divided into multiple cubic voxels of this size. For each voxel, the X-axis, Y-axis, and Z-axis coordinates of all points within the voxel were extracted, and the average values ​​of each axis coordinate, namely Xavg, Yavg, and Zavg, were calculated. The center point with coordinates (Xavg, Yavg, Zavg) was retained, and other redundant points within the voxel were deleted, thereby reducing the amount of point cloud data while retaining the core features of the suture contour.

[0022] The plane equation of the flange base material end face is fitted using a random sampling consensus algorithm. The plane equation is in the form of: Where A, B, and C are the three components of the plane normal vector, and D is a constant term; the specific fitting process is as follows: three non-collinear points are randomly selected from the point cloud, denoted as... , , The plane normal vector is calculated through the cross product of vectors, i.e., vectors. ,vector normal vector Then substitute point P1 into the plane equation to find the constant term. ; satisfying statistical point cloud The number of interior points (threshold distance from point to plane) is counted in millimeters. The process of randomly selecting points, calculating the plane equation, and counting the number of interior points is repeated 100 times. The plane equation with the highest number of interior points is selected as the final fitted plane equation for the flange base material end face. The normal vector direction of this final fitted plane is calculated, i.e., the direction indicated by (A, B, C). The entire point cloud data is projected onto the fitted plane along the normal vector direction, with the projection rule being the projection point of point (X, Y, Z). satisfy Parallel to (A, B, C), and satisfy This process eliminates the interference of uneven end face caused by flange installation tilt on the extraction of circumferential seam contour features, and finally obtains optimized three-dimensional contour data.

[0023] Step 201: Extract geometric feature points of the circumferential seam edge from the optimized 3D contour data. Based on the spatial coordinates of the geometric feature points, obtain the spatial position parameters of two reference feature lines by fitting using the least squares method. Specifically, this includes: calculating the gradient value of each point in the plane for the 2D point cloud data projected onto the fitting plane in step 200 (the coordinates of the projected points are marked as (X, Y), and the Z coordinate is uniformly the corresponding value on the fitting plane). The gradient value reflects the degree of coordinate change between the point and its surrounding points. The calculation method is to take the points within a 3×3 neighborhood around the point and calculate the partial derivative along the X-axis direction, i.e., the ratio of the difference in Y coordinates to the difference in X coordinates of adjacent points, and the partial derivative along the Y-axis direction, i.e., the ratio of the difference in X coordinates to the difference in Y coordinates of adjacent points. The gradient value is equal to the square root of the sum of the squares of these two partial derivatives. The preset threshold for the gradient value is set to 0.8. , gradient value The points are marked as edge candidate points, and these points are initially identified as potential points on the edge of the annular seam. Then, a region growing algorithm is used to cluster the edge candidate points. Specifically, using a candidate point as a seed point, other candidate points within a 1 mm radius are searched. If the distance between two points is less than 0.3 mm and the difference in gradient direction (the angle of the gradient vector) is less than 5°, the latter is assigned to the cluster containing the seed point. This process is repeated until all candidate points are clustered, resulting in two point sets. One set is the inner edge point set of the annular seam, containing n points, each with coordinates marked as... ,in Another type is the point set on the outer edge of the annular seam, which contains m points, each with coordinates denoted as . ,in For both types of edge point sets, the baseline feature line is fitted using the least squares method. Taking the inner edge point set of the circumferential seam as an example, let the corresponding baseline feature line be the first baseline feature line. The spatial parameters of this feature line include two parts: the fixed points through which the feature line passes, denoted as [missing information]. and the direction vector of the feature line (components denoted as) ), corresponding to the X-axis, Y-axis, and Z-axis directions respectively; construct the error function. The formula is ,in For the inner edge point set, the first The parameters corresponding to each point describe the position of that point on the feature line; the function to minimize the error is... , for parameters The partial derivatives are calculated separately, and the specific process is as follows: When finding partial derivatives, focus only on Contains item After differentiation, we get Let the partial derivative equal to 0 (at this time) (Taking the minimum value), we get the equation. ;right When finding partial derivatives, focus only on Contains item Similarly, after differentiating and setting it to 0, we obtain the equation. ;right When finding partial derivatives, focus only on Contains item Similarly, after differentiating and setting it to 0, we obtain the equation. ;right When finding partial derivatives, focus only on Contains item After differentiation, we get Setting the partial derivative to zero, we obtain the equation. ;right When finding partial derivatives, focus only on Contains item Similarly, after differentiating and setting it to 0, we obtain the equation. ;right When finding partial derivatives, focus only on Contains item Similarly, after differentiating and setting it to 0, we obtain the equation. The above six equations form a system of linear equations. Substituting the coordinates of the inner edge point set... After determining the values ​​of , the direction vector components of the first baseline feature line are obtained. and the coordinates of the points passed through That is, the spatial position parameters of the first baseline feature line; the same method is used to process the set of points on the outer edge of the circumferential seam, i.e., to construct the error function. (Form and) Consistent, only the point coordinates are replaced with the outer edge point set ( ), parameters replaced with fixed points ( ), direction vector components ( ) and parameters The formula is ,in It is the first point of the outer edge of the circumferential seam. The parameters corresponding to each point are used to describe the position of that point on the second reference feature line; by adjusting the parameters... , , , , By taking the partial derivatives and setting them to zero, a system of linear equations is constructed and solved to obtain the spatial position parameters of the second baseline characteristic line, i.e., the direction vector components. ), passing through the coordinates of the point ( ).

[0024] Step 202: Calculate the included angle between the two reference feature lines based on their spatial position parameters. Determine the processing range boundary of the circumferential seam based on the included angle value, and define the circumferential seam segment within the included angle range as the local processing area. Specifically, this includes: taking the direction vector of the first reference feature line. (The amount is) ), the direction vector of the second baseline feature line (The amount is) According to the formula Calculate the dot product of two direction vectors According to the formula respectively Calculate the magnitudes of the two direction vectors. and According to the formula Calculate the cosine of the included angle. ; for cosine value According to the formula Perform inverse cosine calculation to obtain the angle between the two reference feature lines. (Unit: degrees); then determine the local processing area based on the included angle value, that is, solve the equations of the first and second reference feature lines simultaneously. The equation of the first reference feature line is: The equation of the second baseline characteristic line is: Solving this system of equations yields the coordinates of the intersection point of the two baseline feature lines. Taking the intersection point as the vertex, extend to both sides along the circumference of the flange circumferential seam, and measure the angle between the two reference feature lines. The covered circumferential seam section is defined as the boundary of the processing range, that is, the boundary is two reference feature lines, and the length and included angle of the circumferential seam section within the boundary are... Correspondingly, the three-dimensional spatial region within this boundary is the local processing region.

[0025] This embodiment effectively filters out noise and outliers by preprocessing the three-dimensional contour data, reducing redundant data interference; it accurately captures the spatial geometry of the circumferential seam edge by fitting the baseline feature lines using the least squares method, solving the problem of feature positioning ambiguity caused by the irregularity of the weld contour; and it defines the local processing area based on the angle between the feature lines, achieving precise focusing on the core section of the circumferential seam, avoiding ineffective milling of non-weld seam areas during processing, which not only improves processing efficiency but also helps ensure the targeting and accuracy consistency of milling.

[0026] In a preferred embodiment of the present invention, step 3 includes: Step 300: Select feature points on the center line of the circumferential seam within the local processing area as the first control points, and select feature points on the reference edge of the flange base material outside the local processing area as the second control points. Specifically, this includes: based on the inner and outer edge point sets of the circumferential seam obtained in step 201, according to the positional correspondence, take each point in the inner edge point set and its corresponding point in the outer edge point set, and calculate the midpoint between the two as the point of the circumferential seam center line. The midpoint coordinates are calculated as follows: the X-axis coordinate of the midpoint is equal to the sum of the X-axis coordinates of the inner and outer edge points divided by 2; the Y-axis coordinate of the midpoint is equal to the sum of the Y-axis coordinates of the inner and outer edge points divided by 2; the Z-axis coordinate of the midpoint is equal to the sum of the Z-axis coordinates of the inner and outer edge points divided by 2. From these midpoints of the circumferential seam center lines, along the local processing area... Select points evenly along the circumference of the circumference seam, for example, one point every 30°, for a total of 12 points. These points will be used as the first control points. When selecting these points, ensure that each first control point falls within the local processing area and covers the entire circumference of the circumference seam to ensure a complete fit between the subsequent path and the circumference seam. Select control points on the flange base material reference edge outside the local processing area. The base material reference edge refers to the stable base material end face edge on the flange that is far from the weld and free from weld allowance interference. This area has a stable geometry and can serve as a reliable processing reference. Select the same number of points as the first control points evenly along the circumference of the flange, for example, 12 points as well. Each point must be at least 10 mm away from the boundary of the local processing area to avoid the influence of the weld processing process. At the same time, ensure that each point falls within a flat area of ​​the base material to ensure its positional stability. These points will be used as the second control points.

[0027] Step 301: Record the continuous position data of the first and second control points during the milling process to form a time series dataset. Based on the continuous position data in the time series dataset, generate a smooth and continuous circular closed control path. Specifically, this includes: acquiring the three-dimensional coordinates of each first and second control point in real time during the milling process at a preset sampling frequency consistent with the multi-degree-of-freedom inertial measurement unit, such as 500Hz; arranging the three-dimensional coordinates of the first and second control points at each moment into a time series dataset, which includes the X-axis, Y-axis, and Z-axis coordinates of the first control point at each moment, and the X-axis, Y-axis, and Z-axis coordinates of the second control point at that moment, ensuring that the data reflects the positional pattern of the control points over time; sorting the time series dataset according to the time axis from the first moment to the last moment to ensure the continuity of the data in the time dimension; and processing the X-axis, Y-axis, and Z-axis coordinate sequences of the first and second control points using a cubic spline interpolation algorithm, constructing an interpolation function with time as the independent variable and coordinates as the dependent variable, the function form being... ,in This refers to the position of a control point (either the first or second control point) at the current interpolation time point, calculated using a cubic spline interpolation algorithm. The coordinates on For the current time point of the interpolation calculation, This represents the start time of the current time interval. , , , These are the interpolation coefficients, which are determined as follows: for two adjacent time points, the previous time point is denoted as [the first time point]. and the current moment The time interval between the two is Using the coordinate values ​​of these two moments, i.e., the coordinate value of the previous moment, is denoted as... The current coordinate value is denoted as And the continuity of the first derivative between adjacent time intervals (the previous interval in The first derivative is denoted as The current interval is The first derivative needs to be with Equal), second derivative continuity (the previous interval is in) The second derivative is denoted as The current interval is The second derivative needs to be with (Equal) Construct a system of equations and solve for the results. The specific value, for example when =0、 Furthermore, the boundary conditions are set to natural boundaries, meaning that when the second derivative at the beginning of the first time interval is 0 and the second derivative at the end of the last time interval is 0, the coefficient values ​​can be obtained by solving the system of equations. In this way, it is ensured that the interpolated coordinate sequence is not only continuous in time, but also that its velocity (first derivative) and acceleration (second derivative) remain continuous, avoiding path fluctuations or abrupt changes and ensuring smooth robot motion. Then, combined with the circumferential geometric constraints of the flange circumferential seam, the interpolated first control point path and second control point path are fused. During fusion, the first control point path is used as the core (ensuring the path fits the circumferential seam and meets the processing position requirements), and the second control point path is used as the reference (using the base material reference to correct flange installation deviations and improve path accuracy). The final path at each moment is calculated with a weight of 70% for the first control point and 30% for the second control point. The coordinates are calculated as follows: final X-axis coordinate = 0.7 × first control point X-axis coordinate + 0.3 × second control point X-axis coordinate; final Y-axis coordinate = 0.7 × first control point Y-axis coordinate + 0.3 × second control point Y-axis coordinate; final Z-axis coordinate = 0.7 × first control point Z-axis coordinate + 0.3 × second control point Z-axis coordinate. This weighted calculation achieves complementary advantages between the two types of paths. The coordinate deviation between the first and last moments of the path is calculated by comparing the coordinate differences in the X-axis, Y-axis, and Z-axis directions. If the deviation in any direction is greater than 0.01 mm, it indicates that the path is not closed at the beginning and end, and the interpolation coefficients need to be adjusted and the interpolation and fusion calculations need to be repeated. If the deviation in all three directions is less than 0.01 mm, the path is connected at the beginning and end, and a smooth and continuous closed loop control path is finally generated.

[0028] In this embodiment, the first control point is selected from the center line of the circumferential seam, directly conforming to the core area of ​​the weld, ensuring precise correspondence between the path and the circumferential seam. The second control point is selected from the reference edge point of the base material. Relying on the stable geometry of the base material, it can offset the reference deviation caused by flange installation tilt and end face runout, providing dual positional protection for the path and reducing the impact of the base coordinate system deviation on the machining. By continuously recording the control point positions and using cubic spline interpolation to generate the path, the position, velocity, and acceleration of the path are continuous, avoiding impacts and vibrations during robot movement, reducing sudden changes in milling force caused by path fluctuations, and ensuring a smooth milling process. The circular closed control path covers the entire circumference of the circumferential seam without breaks, avoiding local missed cuts or overcuts caused by path discontinuity. At the same time, the uniform path smoothness can ensure that the milling parameters of each section of the circumferential seam, such as feed rate and depth of cut, are consistent, improving the uniformity of the machined surface quality.

[0029] In a preferred embodiment, step 4 includes: Step 400: Based on the curvature characteristics and normal vector distribution of the circular closed control path, calculate the geometric characteristic parameters of the circular closed control path. Based on these parameters, generate trajectory correction parameters. Specifically, this includes: selecting a path point every 0.5 mm along the circular closed control path to obtain a continuous sequence of path points; for any three adjacent path points in the sequence, denoted as points A, B, and C, with point B as the intermediate point, calculate vector BA (coordinates of point A minus coordinates of point B) and vector BC (coordinates of point C minus coordinates of point B); calculate the angle between the two vectors by the dot product of vectors BA and BC, and combine the lengths of vectors BA and BC with the triangle area formula (area of ​​triangle ABC = 0.5 × length of vector BA × length of vector BC × sine of the angle between the two vectors). The radius of curvature at point B is derived (radius of curvature = length of vector BA × length of vector BC × sine of the angle between the two vectors / (2 × area of ​​triangle ABC)). The curvature at point B is the reciprocal of the radius of curvature. The curvature of all path points is calculated in this way to form the path curvature distribution parameters. For each path point on the path, denoted as point P, two adjacent path points are selected (denoted as before point P and after point P). The vector before P (coordinates of point P minus coordinates of before point P) and the vector after P (coordinates of point P minus coordinates of point P) are calculated. The cross product operation is performed on the vector before P and the vector after P to obtain a vector perpendicular to the tangent direction of the path. The cross product vector is normalized (the vector length is adjusted to 1) to obtain the path unit normal vector at point P, forming the path normal vector distribution parameters.

[0030] Based on the path curvature distribution parameters, the preset curvature threshold is set to 0.05. This threshold is set based on the milling accuracy requirements of the flange circumferential seam. Path segments with curvature values ​​exceeding this threshold are prone to trajectory deviation due to motion inertia; for curvature values ​​greater than 0.05... For each path segment, the position correction is calculated as follows: Position Correction Amount = Curvature Value × Preset Compensation Coefficient. The preset compensation coefficient is set to 0.02 mm (this coefficient can effectively compensate for trajectory offset at points with large curvature). For curvature values ​​less than or equal to 0.05 mm... The path segment has a position correction value set to 0. Based on the path normal vector distribution parameters, calculate the angle between the normal vector of each path point and the reference normal vector (the preset normal vector along the flange radial direction) in the theoretical programming coordinate system. The specific process is as follows: the normal vector of the path point is a unit vector, denoted as vector N, containing three directional components: X-axis, Y-axis, and Z-axis. The reference normal vector is also a unit vector (denoted as vector B, containing three directional components: X-axis, Y-axis, and Z-axis). Calculate the dot product of the two vectors, i.e., the product of the X-axis component of vector N and the X-axis component of vector B, plus the Y-axis component of vector N. The product of the y-axis component of vector N and the y-axis component of vector B, plus the product of the z-axis component of vector N and the z-axis component of vector B; since both vectors are unit vectors (both with a magnitude of 1), the cosine of the angle between the two vectors is equal to the absolute value of this dot product; the inverse cosine operation is performed on this cosine value to obtain the angle between the two vectors, and this angle is used as the attitude correction amount to adjust the attitude of the robot end effector in the direction of the path normal vector, ensuring that the milling cutter attitude matches the path normal vector and avoiding milling direction deviation; the position correction amount and attitude correction amount of all path points are integrated to form the trajectory correction parameters.

[0031] Step 401: Fuse the trajectory correction parameters with the basis error data to correct the transformation matrix between the robot base coordinate system and the theoretical programming coordinate system, obtaining the corrected transformation matrix. Specifically, this includes: calling the initial transformation matrix (a 4×4 homogeneous transformation matrix containing rotation and translation components, where the rotation component describes the attitude relationship between the two coordinate systems, and the translation component describes the positional relationship) pre-stored in the control system. The 3×3 submatrix in the upper left corner represents the rotation component, describing the attitude correspondence between the two coordinate systems; the first three elements of the fourth column represent the translation component, describing the positional offset relationship between the two coordinate systems; the fourth row is fixed at [0, 0, 0, 1]. The standard structure of the homogeneous transformation matrix; extract the first three elements from the fourth column of the matrix, denoted as initial X translation, initial Y translation, and initial Z translation, respectively; extract the translation deviation from the basis error data generated in step 101, where ΔX is the translation deviation in the X-axis direction, ΔY is the translation deviation in the Y-axis direction, and ΔZ is the translation deviation in the Z-axis direction; extract the position correction from the trajectory correction parameters generated in step 400, where corrected X is the position correction value in the X-axis direction, corrected Y is the position correction value in the Y-axis direction, and corrected Z is the position correction value in the Z-axis direction, and calculate the corrected translation component according to the following formula: Corrected X translation = Initial X translation + ΔX + Corrected X (Initial X translation is inherent to both coordinate systems). Position offset, ΔX compensation for the pose deviation of the moving platform, and X correction compensation for the path trajectory offset are superimposed to obtain the real-time X-axis translation compensation result. Corrected Y translation = initial Y translation + ΔY + corrected Y (the calculation logic is consistent with the X-axis to obtain the real-time Y-axis translation compensation result). Corrected Z translation = initial Z translation + ΔZ + corrected Z (the calculation logic is consistent with the X-axis to obtain the real-time Z-axis translation compensation result). Based on the 3×3 rotation submatrix in the upper left corner, the specific process of calculating the initial X-axis angle, initial Y-axis angle, and initial Z-axis angle according to the XYZ Euler angle transformation rule (first rotate around the X-axis, then around the Y-axis, and finally around the Z-axis) is as follows: First, clarify the 3×3 rotation submatrix. The matrix contains 3 rows and 3 columns of elements. According to the XYZ Euler angle transformation rule, these elements are fixedly related to the three initial angles, and the angle values ​​can be derived step by step through the relationship. In the first step of calculating the initial angle around the Y-axis, according to the transformation rule, the element in the 3rd row and 1st column of the rotation submatrix and the initial angle around the Y-axis satisfy the relationship that the element in the 3rd row and 1st column is equal to the sine value of the initial angle around the Y-axis. First, take the negative value of the element in the 3rd row and 1st column to obtain the sine value of the initial angle around the Y-axis. Then, perform an arcsine operation on this sine value. Combined with the actual range of the rotation angle around the Y-axis in the flange circumferential seam processing scenario (between -90° and 90° to avoid multiple angle solutions), determine the specific value of the initial angle around the Y-axis.In the second step of calculating the initial angle around the X-axis, according to the transformation rules, the element in the 3rd row and 2nd column of the rotation submatrix satisfies the following relationship with the initial angle around the X-axis and the initial angle around the Y-axis: the element in the 3rd row and 2nd column = the sine of the initial angle around the X-axis multiplied by the cosine of the initial angle around the Y-axis. Similarly, the element in the 3rd row and 3rd column satisfies the following relationship with the two angles: the element in the 3rd row and 3rd column = the cosine of the initial angle around the X-axis multiplied by the cosine of the initial angle around the Y-axis. Since the initial angle around the Y-axis has already been calculated in the first step, its cosine value can be obtained using trigonometric functions (and the cosine value is not zero when the angle is between -90° and 90°). Dividing the element in the 3rd row and 2nd column by the element in the 3rd row and 3rd column cancels out the cosine of the initial angle around the Y-axis, yielding the tangent of the initial angle around the X-axis. The arctangent of this tangent is then calculated, and combined with the signs of the elements in the 3rd row and 2nd column and the 3rd row and 3rd column (to determine the initial angle around the X-axis). In the first quadrant, the initial X-axis angle is obtained. In the third step, when calculating the initial Z-axis angle, according to the transformation rules, the element in the first row and first column of the rotation submatrix satisfies the relationship between the initial Y-axis angle and the initial Z-axis angle: the element in the first row and first column = the cosine of the initial Y-axis angle multiplied by the cosine of the initial Z-axis angle. Similarly, the element in the second row and first column satisfies the relationship between the element in the second row and first column = the cosine of the initial Y-axis angle multiplied by the sine of the initial Z-axis angle. Using the calculated initial Y-axis angle, the cosine of the initial Y-axis angle is obtained. Dividing the element in the second row and first column by the element in the first row and first column cancels out the cosine of the initial Y-axis angle, yielding the tangent of the initial Z-axis angle. The arctangent of this tangent is then calculated, and combined with the sign of the elements in the first row and first column and the second row and first column (to determine the quadrant of the initial Z-axis angle), the specific value of the initial Z-axis angle is obtained.

[0032] Extract the rotational deviation from the basis error data generated in step 101, where ΔθX is the rotational deviation around the X-axis, ΔθY is the rotational deviation around the Y-axis, and ΔθZ is the rotational deviation around the Z-axis; extract the attitude correction from the trajectory correction parameters generated in step 400, where correction θX is the attitude correction value around the X-axis, correction θY is the attitude correction value around the Y-axis, and correction θZ is the attitude correction value around the Z-axis. Calculate the corrected rotational components using the following formulas: Corrected X-axis angle = Initial X-axis angle + ΔθX + Corrected θX (The initial X-axis angle is the inherent attitude offset between the two coordinate systems; ΔθX compensates for the attitude deviation of the moving platform; correction θX compensates for the path normal vector deviation; the three are superimposed to obtain the real-time attitude compensation result around the X-axis); Corrected Y-axis angle = Initial Y-axis angle + ΔθY + Corrected θY (Calculation logic...). Consistent with the X-axis, the real-time Y-axis attitude compensation result is obtained; the corrected Z-axis angle = initial Z-axis angle + ΔθZ + corrected θZ (the calculation logic is consistent with the X-axis, obtaining the real-time Z-axis attitude compensation result); the corrected translation components (corrected X translation, corrected Y translation, corrected Z translation) replace the first three elements of the fourth column of the initial matrix respectively; the corrected rotation components (corrected X-axis angle, corrected Y-axis angle, corrected Z-axis angle) are converted into a new 3×3 rotation submatrix according to the XYZ Euler angle transformation rule, replacing the original 3×3 rotation submatrix in the upper left corner of the initial matrix; the fourth row of the matrix remains unchanged [0, 0, 0, 1], and the final corrected transformation matrix can accurately describe the real-time correspondence between the robot base coordinate system and the theoretical programming coordinate system in the current processing state.

[0033] Step 402: Based on the corrected transformation matrix, generate a real-time compensation signal for the base coordinate system. Specifically, this includes: separating the components required for real-time compensation from the corrected transformation matrix, including real-time compensation amounts in the X-axis direction (corrected X translation - initial X translation), Y-axis direction (corrected Y translation - initial Y translation), Z-axis direction (corrected Z translation - initial Z translation), and real-time compensation angles around the X-axis (corrected X-axis angle - initial X-axis angle), Y-axis (corrected Y-axis angle - initial Y-axis angle), and Z-axis (corrected Z-axis angle - initial Z-axis angle); encoding the extracted 6 compensation amounts (3 translation compensation amounts and 3 rotation compensation angles) according to the signal protocol format of the robot controller, converting them into voltage or digital signals; adding a real-time identifier to the encoded signal, including a signal generation timestamp, to ensure that the robot controller can receive and apply the compensation amounts in chronological order, ultimately forming the real-time compensation signal for the base coordinate system.

[0034] This embodiment generates trajectory correction parameters by analyzing geometric characteristics such as path curvature and normal vectors, enabling the correction amount to adapt to the actual shape of the circumferential seam. This avoids local overcutting or undercutting caused by a fixed correction mode, ensuring the fit between the milling path and the circumferential seam contour. The trajectory correction parameters are fused with the base error data to correct the transformation matrix of the two coordinate systems, effectively offsetting the superposition effect of base error and path geometric deviation, and solving the misalignment problem between the robot's base coordinate system and the theoretically programmed coordinate system. The real-time compensation signal can dynamically adjust the robot's base coordinate system, so that the coordinate system matching relationship is updated in real time with the error changes during the machining process. This avoids the decrease in machining accuracy caused by error accumulation, improves the position and attitude control accuracy of the robot throughout the milling process, and ensures the stability of the weld milling quality.

[0035] In a preferred embodiment of the present invention, step 5 includes: Step 500: Calculate the angle compensation amount of each joint of the robot based on the real-time compensation signal of the base coordinate system; adjust the joint motion parameters of the robot based on the angle compensation amount to generate a pre-compensated motion trajectory. Specifically, this includes: acquiring the real-time compensation signal of the base coordinate system generated in step 402, which contains six components: X-axis translation compensation amount ΔX, Y-axis translation compensation amount ΔY, Z-axis translation compensation amount ΔZ, and rotation compensation angles around the X-axis Δα, Y-axis rotation compensation angle Δβ, and Z-axis rotation compensation angle Δγ. These six components are used as compensation increments in Cartesian space. Combined with the robot's link parameters, each joint corresponds to three fixed parameters: the angle α between the axes of two adjacent joints in a plane perpendicular to the link length, the vertical distance a between the two joint axes, and the distance d between the two link planes along the joint axis direction; and one variable parameter: the rotation angle θ of the joint around the axis and the joint type (rotational joint or phasing joint). The Cartesian space compensation amount is converted into the angle compensation amount of each joint through inverse kinematics. Then, adjust the motion parameters to generate a pre-compensated motion trajectory. The specific operation is as follows: First, calculate the angle compensation of each joint of the robot. Calculate the angle compensation of the rotary joints (taking a 6-axis robot as an example). To calculate the angle compensation of the rotary joints, it is necessary to first clarify the forward kinematics relationship (the forward kinematics relationship refers to the relationship of calculating the end-effector pose through joint parameters). The end-effector X-axis position = (the vertical distance a1 between the two joint axes of the first joint × the cosine of the rotation angle θ1 of the first joint around the axis) - (the distance d1 between the two link planes along the joint axis of the first joint × the sine of the rotation angle θ1 of the first joint around the axis × the sine of the angle α1 between the two adjacent joint axes of the first joint in the plane perpendicular to the link length) + (the vertical distance a2 between the two joint axes of the second joint × the cosine of the sum of the rotation angles θ1 and θ2 of the first and second joints around the axis) - ... + the fixed position C of the last link in the X-axis direction. xThe Y-axis position of the end joint = (the vertical distance a1 between the two joint axes of the first joint × the sine of the rotation angle θ1 of the first joint about the axis) + (the distance d1 between the two connecting rod planes along the joint axis of the first joint × the cosine of the rotation angle θ1 of the first joint about the axis × the sine of the angle α1 between the two adjacent joint axes of the first joint in the plane perpendicular to the length of the connecting rod) + (the vertical distance a2 between the two joint axes of the second joint × the sine of the sum of the rotation angle θ1 of the first joint about the axis and the rotation angle θ2 of the second joint about the axis) + ... + the fixed position Cᵧ of the last connecting rod in the Y-axis direction; The Z-axis position of the end joint = (the distance d1 between the two connecting rod planes along the joint axis of the first joint × the cosine of the angle α1 between the two adjacent joint axes of the first joint in the plane perpendicular to the length of the connecting rod) The formula is: (cosine of angle α1) + (distance d2 between the two connecting planes along the joint axis of the second joint × cosine of the angle α2 between the axes of the two adjacent joints of the second joint in the plane perpendicular to the length of the connecting rod) + ... + the fixed position Cz of the last connecting rod in the Z-axis direction; the rotation angle around the X-axis = (sine of the rotation angle θ3 of the third joint around the axis × cosine of the rotation angle θ4 of the fourth joint around the axis × cosine of the rotation angle θ5 of the fifth joint around the axis) + (cosine of the rotation angle θ1 of the first joint around the axis × sine of the rotation angle θ2 of the second joint around the axis × sine of the rotation angle θ6 of the sixth joint around the axis); the logic for the rotation angles around the Y-axis and Z-axis is the same as that around the X-axis, which are trigonometric function combinations of the rotation angles of each joint around the axis.

[0036] After clarifying the forward equation relationship, the angle compensation amount of each joint is inversely calculated based on the end-effector compensation increment (angle compensation amount refers to the increment Δθ of the joint's rotation angle around the axis). A system of inverse trigonometric function equations is then established. Taking the first joint as an example, the original end-effector X-axis position is the original X-axis position X0. After compensation, it needs to reach the original X-axis position X0 + X-axis translation compensation amount ΔX. Combining the relationship between the X-axis position and the rotation angle θ1 of the first joint around the axis in the forward equation, when the angle compensation amount Δθ1 of the first joint is small, the simplified result is that the X-axis translation compensation amount ΔX ≈ -1. The vertical distance between the two joint axes of the joint is a1 × the sine of the original rotation angle θ10 of the first joint × the angle compensation amount Δθ1 of the first joint + the distance d1 between the two connecting rod planes along the joint axis of the first joint × the cosine of the original rotation angle θ10 of the first joint × the sine of the angle α1 between the two adjacent joint axes of the first joint in the plane perpendicular to the length of the connecting rod × the angle compensation amount Δθ1 of the first joint; after ignoring minor terms, the X-axis translation compensation amount ΔX ≈ - the vertical distance between the two joint axes of the first joint. a1 × the sine of the original rotation angle θ10 of the first joint × the angle compensation amount Δθ1 of the first joint; Similarly, in the Y-axis direction: the Y-axis translation compensation amount ΔY ≈ the vertical distance between the two joint axes of the first joint a1 × the cosine of the original rotation angle θ10 of the first joint × the angle compensation amount Δθ1 of the first joint; The two equations constitute the system of equations for the first joint, and solving them yields the angle compensation amount Δθ1 of the first joint ≈ -X-axis translation compensation amount ΔX ÷ (the vertical distance between the two joint axes of the first joint a1 × the original rotation angle of the first joint) The sine value of θ10 (or the angle compensation amount Δθ1 of the first joint ≈ the Y-axis translation compensation amount ΔY ÷ (the vertical distance a1 between the two joint axes of the first joint × the cosine value of the original rotation angle θ10 of the first joint)). Following the same logic, establish a system of equations for the second to sixth joints, such as the angle compensation amount Δθ2 of the second joint ≈ the Z-axis translation compensation amount ΔZ ÷ (the vertical distance a2 between the two joint axes of the second joint × the cosine value of the sum of the original rotation angle θ10 of the first joint and the original rotation angle θ20 of the second joint), etc.

[0037] After establishing the system of equations, solve for the angle compensation amount, i.e., calculate the initial solution using the arcsine, arccosine, or arctangent functions. For example, the initial solution for the angle compensation amount Δθ1 of the first joint is = arcsin(-X-axis translation compensation amount ΔX ÷ (vertical distance a1 between the two joint axes of the first joint × angle compensation amount Δθ1 of the first joint)). Combine this with the joint motion range, such as the motion range of the first joint being -180° to 180° and the motion range of the second joint being -90° to 90°, to screen for valid solutions (valid solutions refer to solutions within the joint motion range). The valid solution is the angle compensation amount of the rotary joint. Substitute it into the forward equation relationship to verify and ensure that the end effector can achieve compensation. The target pose is determined by the following steps: If the robot contains translational joints, such as joints that translate along the X-axis, their joint angles are fixed (joint angles refer to the rotation angles of the translational joint, which remain constant), and the displacement is variable (displacement refers to the distance the translational joint moves along the axis). The end effector's corresponding axis position = the displacement of the translational joint + the contribution of the fixed positions of other links (the contribution of the fixed positions refers to the fixed positions of other links in the corresponding axis direction). The calculation steps are as follows: determine the Cartesian axis corresponding to the translational joint. If it only affects the X-axis, that is, the displacement of the translational joint only changes the end effector's X-axis position; obtain the transmission ratio k of the translational joint, that is, k = the number of teeth on the driven gear ÷ the number of teeth on the driving gear in the gear transmission, such as k = 1.2; the translational joint's... Displacement compensation amount = translation compensation amount of the corresponding axis × transmission ratio k. For example, when the X-axis translation compensation amount ΔX = 0.2 mm, the displacement compensation amount of the traverse joint = 0.2 mm × 1.2 = 0.24 mm, ensuring that the displacement of the traverse joint can be accurately converted into the compensation displacement of the corresponding end axis. After obtaining the angle compensation amount of each joint, the joint parameters in the robot's original motion plan are adjusted. The original target angle of the rotary joint, the original target displacement of the traverse joint, and the motion velocity of all joints are extracted from the original motion plan, i.e., the change in joint angle or displacement per unit time, and the motion acceleration a, i.e., the change in joint velocity per unit time. After adjustment, the target angle = rotary joint... The original target angle is calculated by adding the angle compensation amount of the rotary joint. For example, if the original target angle of the first joint is 30° and the angle compensation amount of the first joint is 2°, then the adjusted target angle of the first joint is 30° + 2° = 32°. The adjusted target displacement is calculated by adding the original target displacement of the locating joint to the displacement compensation amount of the locating joint. For example, if the original target displacement of the first locating joint is 50 mm and the displacement compensation amount of the first locating joint is 0.24 mm, then the adjusted target displacement of the first locating joint is 50 mm + 0.24 mm = 50.24 mm. Maintaining constant motion speed and acceleration avoids sudden parameter changes that could cause impacts or fluctuations during robot movement, ensuring smooth motion.

[0038] To verify whether the adjusted joint parameters can enable the end effector to reach the target pose, the following steps are taken: Substitute the adjusted parameters of all joints (adjusted target angle of the rotary joint, adjusted target displacement of the translating joint) into the kinematic forward kinematics relation to calculate the actual position (actual X-axis coordinate Xreal, actual Y-axis coordinate Yreal, actual Z-axis coordinate Zreal) and actual attitude (actual rotation angle α around the X-axis, actual rotation angle β around the Y-axis, actual rotation angle γ around the Z-axis) of the end effector. Extract the original trajectory pose of the end effector from the original motion plan (original X-axis position, original Y-axis position, original Z-axis position, original rotation angle around the X-axis, original rotation angle around the Y-axis, original rotation angle around the Z-axis). Calculate the theoretical position according to the following formulas: Theoretical X-axis coordinate X = Original X-axis position + X-axis translation compensation; Theoretical Y-axis coordinate Y = Original Y-axis position + Y-axis translation compensation; Theoretical Z-axis coordinate Z = Original Z-axis position + Z-axis translation compensation. Calculate the theoretical attitude according to the following formulas: Theoretical X-axis rotation angle α = Original X-axis rotation angle + X-axis rotation compensation angle; Theoretical Y-axis rotation angle β = Original Y-axis rotation angle + Y-axis rotation compensation angle; Theoretical Z-axis rotation angle γ = Original Z-axis rotation angle + Z-axis rotation compensation angle. Calculate the deviation between the actual pose and the theoretical pose, i.e., position deviation |Xactual - Xtheoretical| < 0.005 mm, |Yactual - Ytheoretical| < 0.005 mm, |Zactual - Ztheoretical| < 0.005 mm. Deviations |αactual - αtheoretical| < 0.001°, |βactual - βtheoretical| < 0.001°, |γactual - γtheoretical| < 0.001°; if not satisfied, re-examine the establishment and solution process of the inverse equation system, correct the joint angle compensation amount, and verify again until the deviation meets the requirements. Integrate the adjusted target parameters of all joints (adjusted target angle of rotary joints, adjusted target displacement of translating joints), original motion velocity, and original motion acceleration, and arrange the target parameters of each joint at each moment in chronological order (from the start time of motion to the end time of motion). For example, when the motion is 0.1 seconds, the target angle of the first joint = 30.5°, the target angle of the second joint = 15°; when the motion is 0.2 seconds, the target angle of the first joint = 31°, the target angle of the second joint = 16°, etc., to form a pre-compensated motion trajectory; this trajectory contains the motion commands of each joint of the robot at each moment.

[0039] Step 501: During the robot's movement along the pre-compensated motion trajectory, a six-dimensional force sensor is used to collect real-time three-dimensional force data between the milling cutter and the workpiece. This three-dimensional force data is then filtered to obtain milling force feedback data. Specifically, the six-dimensional force sensor is installed between the robot's end effector and the milling cutter to detect the force generated by the contact between the milling cutter and the workpiece during milling. Force data is collected at a frequency consistent with the sampling frequency of the robot control system, such as 1kHz. This data includes force components in three-dimensional space (force in the X-axis direction, force in the Y-axis direction, and force in the Z-axis direction) and three-dimensional torque components (torque around the X-axis, torque around the Y-axis, and torque around the Z-axis). Since the force signal has a more direct impact on machining quality during milling, only the three-dimensional force components are extracted to form the original three-dimensional force feedback data. Force dataset; the raw three-dimensional force data is processed using a moving average filtering algorithm. Specifically, a sliding window of 50 sampling points is set. For each three-dimensional force component (X-axis force, Y-axis force, Z-axis force), the average value of the 50 sampling points within its respective sliding window is calculated. The calculation method is: average X-axis force = sum of all X-axis force values ​​within the window / 50. The calculation methods for the average Y-axis and Z-axis forces are the same as for the average X-axis force. The calculated window average value replaces the corresponding raw force value at that moment. The window is slid sequentially over time to process the raw force values ​​of all sampling points, ultimately obtaining smoothed milling force feedback data. This data accurately reflects the stable change trend of the force during milling.

[0040] In this embodiment, by converting the base coordinate system compensation signal into joint angle compensation, the robot's motion parameters are directly corrected, enabling the pre-compensated trajectory to adapt to coordinate system deviations and path errors in real time. This avoids end effector pose shifts caused by accumulated errors and ensures that the milling cutter always moves along the target path. The collected and filtered milling force feedback data can accurately reflect the contact state between the milling cutter and the workpiece, avoiding tool wear or workpiece damage caused by sudden load changes. The pre-compensated motion trajectory reduces trajectory deviations in robot motion, and the filtered force signal eliminates noise interference. The combination of the two makes position and force control in the milling process more precise, reducing the problem of unstable machining quality caused by trajectory fluctuations or force signal jitter.

[0041] In a preferred embodiment of the present invention, step 6 includes: Step 600: Perform frequency domain analysis on the milling force feedback data to extract vibration characteristic frequencies and amplitudes; perform harmonic analysis on the current signal of the spindle motor to identify load change characteristics. Specifically, this includes: acquiring the milling force feedback data generated in step 501, including smooth force signals in the X, Y, and Z axes; determining the analysis duration to be 1 second; and extracting continuous force signal segments according to this duration; if the sampling frequency is 1kHz, 1000 time-domain sampling points will be collected within 1 second, with the force value at each point being 4... 0.2N, 4.3N, 4.1N, ..., 4.4N (a total of 1000 values, in N), ensuring the data length matches the number of points in the subsequent Discrete Fourier Transform; For the truncated single axis, such as the X-axis force signal, a Discrete Fourier Transform is performed. The Fourier coefficients of the time-domain signal are complex numbers, and their values ​​are calculated from the force values ​​at each time-domain sampling point. Taking the 250th frequency point (frequency 250Hz) as an example, the corresponding Fourier coefficients have a real part of 320 N·s and an imaginary part of -110 N·s; taking the 3rd... Taking 20 frequency points (320Hz) as an example, corresponding to a real part of 380N·s and an imaginary part of -130N·s for the Fourier coefficients, the frequency domain amplitude of each frequency point needs to be calculated as follows: First, calculate the absolute value of the Fourier coefficients, i.e., take the sum of the squares of the real and imaginary parts of the Fourier coefficients, and then take the square root of the result. Because the Discrete Fourier Transform has conjugate symmetry for the spectrum of real signals, such as milling force signals, the amplitudes of positive and negative frequencies are symmetrically distributed. In actual analysis, only the positive frequency part is considered. At this time, the amplitude of the positive frequency part is only half of the true amplitude of the signal. Therefore, the absolute value of the Fourier coefficients needs to be multiplied by 2 to restore the true amplitude corresponding to the positive frequency. Finally, divide by the number of data points involved in the transformation, 1000, which is the total number of time domain sampling points, to obtain the final frequency domain amplitude. The calculation formula is: final frequency domain amplitude = (absolute value of Fourier coefficients × 2) ÷ number of data points. Calculate the frequency domain amplitude of specific frequency points according to the above steps. For example, when calculating the 250Hz frequency point, the absolute value of the Fourier coefficients = = ≈338 N·s, multiplied by 2, we get 338 × 2 = 676 N·s, then divided by the number of data points (1000), the final frequency domain amplitude = 676 ÷ 1000 = 0.676 N; calculate the frequency domain amplitude of all 1000 frequency points (frequency range 0 to 500 Hz) on the X-axis using the above method, and establish the frequency-frequency domain amplitude correspondence; determine the maximum frequency domain amplitude in this set of data, such as 0.808 N corresponding to 320 Hz, and preset the threshold to 30% of the maximum frequency domain amplitude, i.e., 0.8 0.8 × 30% ≈ 0.242 N; Traverse all frequency points and filter out those with frequency domain amplitude greater than the threshold, such as 250 Hz, 320 Hz, etc. These frequency points are the vibration characteristic frequencies of the X-axis; The frequency domain amplitude corresponding to each vibration characteristic frequency, such as 0.676 N for 250 Hz and 0.808 N for 320 Hz, is the vibration characteristic amplitude at that frequency; Process the Y-axis and Z-axis force signals in the same way to obtain the vibration characteristic frequencies and amplitudes of the Y-axis and Z-axis respectively.

[0042] Real-time current signals from the spindle motor are acquired using a current sensor. The sampling frequency is set to match the sampling frequency (1kHz) of the milling force feedback data, and the acquisition duration is the same as the duration of the milling force frequency domain analysis (1 second), ensuring that the number of current signal data points acquired is 1000, matching the number of points for the subsequent discrete Fourier transform. A low-pass filter with a cutoff frequency of 500Hz is used to filter the acquired current signals, removing high-frequency interference components (such as high-frequency electromagnetic noise) and retaining only the fundamental and low-order harmonic components. Low-order harmonics refer to the 2nd and 3rd harmonics, as higher-order harmonics have lower sensitivity to changes in spindle motor load and cannot be effectively filtered. To effectively reflect changes in load status, these are discarded. After signal acquisition and filtering, a Discrete Fourier Transform (DFT) is performed on the processed current signal to separate the fundamental wave and harmonics. The fundamental wave frequency is equal to the rated frequency of the spindle motor, such as 50Hz. The frequencies of each harmonic are integer multiples of the fundamental wave frequency, such as 100Hz for the second harmonic and 150Hz for the third harmonic. The Fourier coefficients of the time-domain signal are as follows: the fundamental wave corresponds to a real part of 400 A·s and an imaginary part of -80 A·s; the second harmonic corresponds to a real part of 120 A·s and an imaginary part of -30 A·s; and the third harmonic corresponds to a real part of 80 A·s and an imaginary part of -20 A·s. When calculating the amplitudes of the fundamental wave and each harmonic, The same steps as those used in the frequency domain analysis of milling force are employed (since the current signal is a real signal, its spectrum also exhibits conjugate symmetry, requiring multiplication by 2 to restore the true amplitude of the positive frequency). Specifically, the steps are as follows: Calculate the absolute value of the Fourier coefficients, i.e., sum the squares of the real and imaginary parts, and then take the square root of the result; multiply the absolute value of the Fourier coefficients by 2 to restore the true amplitude corresponding to the positive frequency; divide by the number of data points (1000) to obtain the final amplitude, calculated using the formula: Final amplitude = (Absolute value of Fourier coefficients × 2) ÷ Number of data points. Calculate the specific harmonics, such as the fundamental, second, and third harmonic amplitudes, following the same steps. When calculating the amplitude proportion of each harmonic, divide the amplitude of that harmonic by... The amplitude of the fundamental wave is multiplied by 100%, i.e., the amplitude ratio of the 2nd harmonic = (2nd harmonic amplitude ÷ fundamental wave amplitude) × 100%, and the amplitude ratio of the 3rd harmonic = (3rd harmonic amplitude ÷ fundamental wave amplitude) × 100%. The amplitude ratios of the fundamental wave and the 2nd and 3rd harmonics are used as characteristic parameters of the spindle motor load change. During milling, when the milling load increases, the fundamental wave amplitude increases accordingly, and the amplitude ratios of the 2nd and 3rd harmonics also increase synchronously. When the load decreases, the fundamental wave amplitude decreases, and the amplitude ratios of the 2nd and 3rd harmonics decrease accordingly. By monitoring the changing trends of these parameters in real time, the load change characteristics of the spindle motor can be identified.

[0043] Step 601: Fuse the vibration characteristic frequency and amplitude with the load change characteristics to obtain load vibration characteristic data. Specifically, this includes: taking the historical maximum value of all vibration characteristic amplitudes for a certain axis, such as the X-axis. This value is the maximum vibration characteristic amplitude recorded during past milling processes on that axis, and is used as a standardization benchmark to ensure the comparability of vibration characteristic amplitudes under different periods and working conditions. The formula for calculating the standardized value of a vibration characteristic amplitude is: Standardized value of a vibration characteristic amplitude = (Vibration characteristic amplitude ÷ Historical maximum value) × Scaling factor 0.5; where 0.5 controls the range of the standardized value between 0 and 0.5. For example, if the vibration characteristic amplitude of the X-axis is 0.676N, its historical maximum value is 1.35N. For 2N, the standardized value = (0.676 ÷ 1.352) × 0.5 = 0.25; Take the fundamental amplitude of the spindle motor under rated load. This value is the fundamental current amplitude of the motor under rated load conditions, used as the benchmark for load characteristic standardization, eliminating the magnitude difference of the fundamental amplitude under different load levels. The formula for calculating the standardized value of the current fundamental amplitude is: Standardized value of the current fundamental amplitude = (Current fundamental amplitude ÷ Fundamental amplitude under rated load) × Scaling factor 0.5; Here, the function of 0.5 is the same as that of vibration characteristic standardization, both controlling the range of the standardized value between 0 and 0.5, matching the range of the vibration characteristic standardization value, and avoiding deviations in the fusion result due to differences in the numerical range. For example, the current fundamental amplitude... The value is 0.816A, and the fundamental amplitude under rated load is 1.632A. Therefore, the standardized value = (0.816 ÷ 1.632) × 0.5 = 0.25. Considering that vibration characteristics directly reflect the stability of the milling process (excessive vibration will lead to a decrease in machining accuracy), and load characteristics are an important cause of vibration changes (load fluctuations will directly affect the vibration state), to balance the influence of both on the milling state, the weight of the standardized vibration characteristic value is set to 0.6, and the weight of the standardized load characteristic value is set to 0.4. The standardized vibration characteristic value of each axis is weighted and summed with the corresponding standardized load characteristic value, that is, the load vibration characteristic value of a certain axis = (standardized vibration characteristic value of that axis × 0.6) + (standardized load characteristic value × 0.4). Since both have been standardized to 0 to 0.5, the range of the load vibration characteristic value after weighted summation is 0 to 0.5 (0.5×0.6+0.5×0.4=0.5). For example, if the standardized value of the X-axis vibration characteristic is 0.25 and the standardized value of the load characteristic is 0.25, then the X-axis load vibration characteristic value = (0.25×0.6) + (0.25×0.4) = 0.15 + 0.1 = 0.25. At the same time, record the vibration characteristic frequency corresponding to each axis, such as 250Hz for the X-axis and 180Hz for the Y-axis. Together with the load vibration characteristic value of each axis, they constitute the load vibration characteristic data, such as X-axis frequency 250Hz, load vibration characteristic value 0.25; Y-axis frequency 180Hz, load vibration characteristic value 0.23, etc.

[0044] Step 602: Based on the load vibration characteristic data, calculate the pose adjustment of the robot end effector using an impedance control algorithm. Based on the pose adjustment, generate a force control signal, specifically including: setting three key matrices required for impedance control (all diagonal matrices; since each axis moves independently, off-diagonal elements are 0) and the desired milling force for each axis. The stiffness matrix K represents the robot end effector's ability to resist pose changes under force; the stiffness coefficients for the X and Y axes are both 500 N / m, and the stiffness coefficient for the Z axis is 600 N / m (the Z axis is the main force direction for milling, so the stiffness coefficient is slightly higher to ensure stability). The damping matrix D represents the velocity decay characteristics of the robot end effector during movement; the damping coefficients for the X and Y axes are both 20 N·s / m, and the damping coefficient for the Z axis is 25 N·s / m (matching the stiffness coefficients to avoid overshoot). The inertia matrix M represents the inertia characteristics of the robot end effector during movement; the inertia coefficients for the X and Y axes are both 0.5 kg, and the inertia coefficient for the Z axis is 0.6kg (matching the axis load, reflecting the relationship between motion acceleration and force); the expected milling force is preset according to the milling process requirements, X-axis is 15N, Y-axis is 12N, Z-axis is 20N (ensuring that the milling force meets the machining requirements and does not exceed the tool's load limit); the core of impedance control is to establish the relationship between the pose adjustment amount (including acceleration and velocity) and the force error through equations. The equation expression is: second derivative of pose adjustment amount × inertia matrix + first derivative of pose adjustment amount × damping matrix + pose adjustment amount × stiffness matrix = expected milling force - actual milling force. The definitions and sources of each variable are as follows: the second derivative of pose adjustment amount, i.e., the acceleration of the robot's end-effector pose adjustment, reflects the rate of change of the adjustment speed; the first derivative of pose adjustment amount × damping matrix + first derivative of pose adjustment amount × stiffness matrix = expected milling force - actual milling force. The second derivative, i.e., the speed of robot end-effector pose adjustment, reflects the speed of the adjustment process; the pose adjustment amount, i.e., the displacement (translation) or rotation (attitude) that the robot end-effector needs to correct, is the final target value to be calculated; the expected milling force minus the actual milling force is defined as the force error, used to reflect the deviation between the current milling force and the target force; the actual milling force is taken from the milling force feedback data collected in step 501 (updated in real time), ensuring that the force error can dynamically reflect the current machining state; because the purpose of pose adjustment during milling is to counteract low-frequency load vibration (rather than rapid movement), the change rate of the end-effector pose adjustment amount is slow, resulting in a very small (approaching 0) value for the second derivative (acceleration) of the pose adjustment amount, and its influence on the overall equation can be ignored. The impedance equation simplifies to: Position adjustment amount × stiffness matrix = force error - first derivative of position adjustment amount × damping matrix. Since each axis moves independently, a formula for calculating the position adjustment amount can be derived for each axis (taking a specific axis as an example): Position adjustment amount of a certain axis = (force error of that axis - first derivative of the position adjustment amount of that axis × damping coefficient of that axis) ÷ stiffness coefficient of that axis. In this formula, the first derivative of the position adjustment amount of that axis, i.e., the adjustment speed of that axis, needs to be set in conjunction with the actual milling rhythm. If the speed is too fast, it is easy to induce new vibrations; if the speed is too slow, it cannot timely counteract load vibrations. It is set to 0.01 to 0.05 m / s (translation adjustment) or 0.01~0.05° / s (attitude adjustment). Taking the calculation of the X-axis translational position adjustment amount as an example, specifically... The steps are as follows: Given that the desired milling force on the X-axis is 15N, and the actual milling force on the X-axis feedback from step 501 is 13N, the X-axis force error = 15N - 13N = 2N. Considering the milling rhythm, such as the low-frequency vibration period of aluminum alloy milling (approximately 0.1 to 0.5 seconds), the first derivative of the X-axis pose adjustment (adjustment speed) is set to 0.01m / s. The X-axis pose adjustment amount = (2N - 0.01m / s × 20N・s / m) ÷ 500N / m = 0.0036m, which is 3.6mm. Following the same logic, the translation pose adjustment amounts for the Y-axis and Z-axis are calculated respectively. The desired milling force on the Y-axis is 12N, the actual milling force is 10N (feedback from step 501), and the force error is 2N; the adjustment speed is 0.0.01 m / s, damping coefficient 20 N·s / m, stiffness coefficient 500 N / m; Y-axis pose adjustment = (2 - 0.01 × 20) ÷ 500 = 0.0036 m (3.6 mm); Z-axis desired milling force 20 N, actual milling force 18 N (feedback in step 501), force error 2 N; adjustment speed 0.01 m / s, damping coefficient 25 N·s / m, stiffness coefficient 600 N / m; Z-axis pose adjustment = (2 - 0.01 × 25) ÷ 600 ≈ 0.00292 m (2.92 mm); The calculation logic for the attitude adjustment around the X-axis, Y-axis, and Z-axis is the same as that for the translation adjustment. Only the parameters of the corresponding axes need to be replaced. The stiffness coefficient, damping coefficient, and inertia coefficient of the attitude adjustment need to match the rotational motion characteristics (units are N·m / °, N·m·s / °, and kg·m, respectively). 2The desired milling force for attitude adjustment is replaced with the desired attitude angle, the actual milling force is replaced with the actual attitude angle (taken from the end effector attitude sensor feedback), and the force error is replaced with the attitude angle error. The final calculated adjustment amount is in degrees. The robot controller can only recognize signals of a specific format. The calculated translation adjustment amounts for each axis, such as 3.6mm for the X-axis, 3.6mm for the Y-axis, and 2.92mm for the Z-axis, and the attitude adjustment amounts, such as 0.2° around the X-axis, 0.15° around the Y-axis, and 0.25° around the Z-axis, need to be converted into a format supported by the controller. This is done in two specific scenarios. If the controller supports digital commands, the adjustment amount is converted into a binary digital command. For example, a 16-bit binary code is used, with a value range of 0 to 65535 corresponding to an adjustment amount of -5mm to 5mm (translation) or -0.5° to 0.5° (attitude). 3.6mm corresponds to binary code 45875 (calculated through linear mapping). If the controller supports analog voltage signals, the voltage-adjustment amount is converted according to a linear proportional relationship. 0V corresponds to 0mm (or 0°) adjustment, and 10V corresponds to 10mm (or 1°) adjustment. Therefore, 3.6mm on the X-axis corresponds to... A 3.6V analog voltage is used, with 0.2V corresponding to a 0.2° rotation around the X-axis, ensuring a one-to-one correspondence between the adjustment amount and the signal value without deviation. To avoid exceeding the robot's motion limits, which could damage the mechanism or cause abnormal accuracy, the converted signal needs to be verified. The specific standards and processing methods are as follows: the safety limit for the robot's end effector translational motion is ±5mm (determined by the mechanical structure design; exceeding this limit will cause link collisions or guide rail over-limits). If the calculated adjustment amount is 6mm (exceeding the upper limit), the signal will be truncated to the maximum allowable value corresponding to 5mm; the robot's end effector posture... The accuracy limit of the dynamic motion is ±0.5° (determined by the accuracy of the servo system; exceeding this will cause the milling surface to tilt out of tolerance). If the calculated adjustment amount is 0.6° (exceeding the upper limit), the signal is truncated to the maximum allowable value corresponding to 0.5°. After verification, the force control signal is output to the robot control system. The control system drives the motors of each joint in real time according to the signal to adjust the translational posture of the end effector, such as moving 3.6mm along the X-axis and the posture, such as rotating 0.2° around the X-axis, thereby offsetting the influence of load vibration on the milling process and ensuring the relative position stability of the tool and the workpiece.

[0045] This embodiment extracts vibration and load features from force and current signals through frequency domain analysis and harmonic analysis, avoiding the limitations of single signal analysis and providing a more comprehensive and accurate reflection of the stability and load state of the milling process. Feature fusion integrates the correlation information between vibration and load through standardization and weighting, reducing irrelevant interference and enabling the load vibration feature data to more directly reflect the core state of the milling process. By establishing the correlation between force and pose, the pose adjustment amount can be calculated in real time based on the load vibration features. The generated force control signal can dynamically offset the influence of vibration and load changes, avoiding overcutting, undercutting, or tool wear, ensuring consistent machining quality, and improving the adaptability and stability of robot motion.

[0046] In a preferred embodiment of the present invention, step 7 includes: Step 700: Based on the pose adjustment amount contained in the force control signal, the pre-compensated motion trajectory is corrected to obtain trajectory correction data. Specifically, the pre-compensated motion trajectory is an initial processing trajectory preset according to the position and shape of the flange weld. The trajectory records the pose information of the robot end effector at fixed time intervals, such as every 0.01 seconds, including the spatial coordinates at that moment, i.e., X-axis coordinates, Y-axis coordinates, and Z-axis coordinates, all in meters, and the attitude angles, i.e., the rotation angles around the X-axis, Y-axis, and Z-axis, all in degrees. For example, at a certain moment, the pre-compensated trajectory... The robot's end effector pose recorded by the trajectory is 0.5 meters on the X-axis, 0.3 meters on the Y-axis, and 0.1 meters on the Z-axis, with an attitude angle of 0 degrees around each axis. The pose adjustment amounts corresponding to the pre-compensated trajectory recording time are extracted from the force control signal. For each recording time, the translation adjustment amounts (in meters) on the X, Y, and Z axes, as well as the attitude adjustment amounts (in degrees) around each axis, are obtained. For example, the adjustment amounts corresponding to the 0.01-second time point mentioned above are 0.0036 meters for X-axis translation, 0.0025 meters for Y-axis translation, and 0.0025 meters for Z-axis translation. The translation adjustment is 0.0018 meters, the attitude adjustment around the X-axis is 0.1 degrees, the attitude adjustment around the Y-axis is 0.08 degrees, and the attitude adjustment around the Z-axis is 0.05 degrees. Following the rule that the corrected parameter = the corresponding parameter of the pre-compensation trajectory + the corresponding pose adjustment, the corrected pose parameters at each moment are calculated. The corrected X-axis coordinate = the X-axis coordinate of the pre-compensation trajectory at that moment + the X-axis translation adjustment, for example, 0.5 meters + 0.0036 meters = 0.5036 meters; the corrected Y-axis coordinate = the Y-axis coordinate of the pre-compensation trajectory at that moment + the Y-axis translation adjustment, for example, 0.0018 meters + 0.0036 meters + 0.0036 meters = 0.0036 meters. 0.3m + 0.0025m = 0.3025m; Corrected Z-axis coordinate = Z-axis coordinate of the pre-compensation trajectory at this moment + Z-axis translation adjustment, for example, 0.1m + 0.0018m = 0.1018m; Corrected X-axis attitude angle = X-axis attitude angle of the pre-compensation trajectory at this moment + X-axis attitude adjustment, for example, 0 degrees + 0.1 degrees = 0.1 degrees; Calculate the corrected Y-axis and Z-axis attitude angles in the same way, and finally arrange the corrected coordinates and attitude angles of all moments in chronological order to obtain the trajectory correction data.

[0047] Step 701: Based on the trajectory correction data, perform trajectory optimization calculations on the pre-compensated motion trajectory to generate trajectory optimization parameters; based on the trajectory optimization parameters, generate the optimized motion trajectory. Specifically, the pre-compensated motion trajectory optimization prioritizes improving trajectory smoothness as the core objective, while meeting machining accuracy requirements. The specific requirements and implementation logic are as follows: During the optimization process, the deviation between the optimized trajectory and the preset machining path of the flange weld should be controlled within the allowable range of the process. Specifically, the deviation in the translation direction (X-axis, Y-axis, Z-axis) should not exceed 0.002 meters, and the deviation in the attitude direction (around the X-axis, Y-axis, Z-axis) should not exceed 0.1 degrees. This deviation range is set to avoid... Trajectory optimization caused the milling cutter to deviate from the weld seam area; the velocity fluctuation in the translational directions (X-axis, Y-axis, Z-axis) needs to be controlled within 0.005 meters per second, and the velocity fluctuation in the attitude directions (around the X-axis, Y-axis, Z-axis) needs to be controlled within 0.02 degrees per second. The core purpose of this limitation is to avoid impact vibration caused by sudden increases or decreases in speed of the moving platform or robot joints. If the velocity fluctuation is too large, the vibration will be transmitted to the milling cutter, causing ripples on the milled surface, which will affect the machining stability and surface quality, and fail to meet the flatness requirements after weld seam milling. The least squares method is used to linearly fit the trajectory correction data, and the fitted straight line replaces the original discrete correction data. The trajectory is now smooth. The specific calculation process is carried out according to the following steps: First, select multiple consecutive time points from the trajectory correction data, such as 100 time points. Selecting consecutive time points ensures the temporal correlation of the data and avoids deviation between the fitting result and the actual motion trend due to excessively large time intervals. At the same time, extract two key data points corresponding to these 100 time points: one is the specific time value of each time point, such as 0.01 seconds for the first time point, 0.02 seconds for the second time point, and so on, up to 1 second for the 100th time point; the other is the X-axis correction coordinate corresponding to each time point, that is, the X-axis correction position calculated in step 700, such as 0.5036 meters for the first time point, 0.5036 meters for the second time point, and so on, up to 1 second for the 100th time point. The value is 0.5042 meters...; A fitting benchmark is established by calculating two basic average values. The time values ​​of 100 moments are added together to obtain the total time. The total time is then divided by the total number of moments (100) to obtain the average time. For example, if the total time of 100 moments is 50.5 seconds, then the average time = 50.5 seconds ÷ 100 = 0.505 seconds; The X-axis correction coordinates of 100 moments are added together to obtain the total X-axis correction coordinates. The total X-axis correction coordinates are then divided by the total number of moments (100) to obtain the average X-axis correction coordinates. For example, if the total X-axis correction coordinates are 50.38 meters, then the average X-axis correction coordinates = 50.38 meters ÷ 100 = 0.5038 meters; The slope is used to reflect the average velocity of the X-axis translation trajectory. The calculation is performed in the following four steps: First, for each moment, calculate the difference between the time value and the average time value, denoted as Δt, and the difference between the corrected X-axis coordinate and the average corrected X-axis coordinate, denoted as ΔX. This difference calculation eliminates the influence of absolute values ​​and focuses on the trend of data change over time. Second, for each moment, multiply the corresponding Δt and ΔX to obtain the difference product for that moment, and then add all the difference products for 100 moments to obtain the sum of the difference products. Third, for each moment, square the corresponding Δt to obtain the squared value of Δt for that moment, and then add all the squared values ​​of Δt for 100 moments to obtain the sum of the squared values ​​of Δt. Fourth, divide the sum of the difference products by the sum of the squared values ​​of Δt to obtain the slope of the fitted line. For example, if the sum of the difference products is 0.025 m / s, the sum of the squared values ​​of Δt is 0.05. Therefore, the slope = 0.025 m·s ÷ 0.05 =0.5 m / s; the intercept is used to reflect the X-axis coordinate reference of the fitted straight line at time 0. The calculation method is to use the average value of the corrected X-axis coordinates - (slope × average time value). For example, if the average value of the corrected X-axis coordinates is 0.5038 m, the slope is 0.5 m / s, and the average time value is 0.505 s, then the intercept = 0.5038 m - (0.5 m / s × 0.505 s) = 0.2513 m; the slope and intercept obtained above together constitute the optimization parameters of the X-axis translation trajectory. The slope determines the average movement speed of the X-axis, ensuring that the X-axis movement speed meets the smoothness requirements; the intercept determines the initial position reference of the X-axis trajectory, ensuring that the X-axis trajectory meets the accuracy requirements; the combination of the two can determine the linear fitting function of the X-axis translation trajectory; according to the optimization calculation logic of the X-axis translation trajectory, the optimization parameters in the following directions are calculated respectively, namely the Y-axis, Z ... The axis translation trajectory also selects continuous time points, such as 100 time points. The time values ​​and corrected coordinates of the corresponding time points are extracted, and the average time value, the average corrected coordinate value, the slope, and the intercept are calculated sequentially to obtain the optimized parameters (slope + intercept) of the Y-axis and Z-axis translation trajectories. This ensures that the Y-axis and Z-axis trajectories meet both accuracy and smoothness requirements. The optimization logic of the attitude trajectory is completely consistent with that of the translation trajectory. The time values ​​of continuous time points and the corrected attitude angles are selected, and a linear function is fitted using the least squares method to calculate the optimized parameters for each attitude direction, namely the slope and intercept. The slope reflects the average rate of change of the attitude angle, ensuring that the attitude motion meets the smoothness requirements. The intercept reflects the initial angle reference of the attitude, ensuring that the attitude trajectory meets the accuracy requirements. The optimized parameters of the X-axis, Y-axis, and Z-axis translation trajectories and the attitude trajectories around each axis are integrated to form a complete set of trajectory optimization parameters.

[0048] By setting more frequent time intervals, such as every 0.001 seconds, the continuity of the trajectory can be ensured, avoiding discontinuities in the trajectory due to excessively large time intervals, which would affect the smoothness of the motion. Then, for each frequent moment, the time value of that moment is substituted into the linear fitting function of the corresponding direction. For example, the X-axis fitting function is X-axis target coordinate = slope × time value + intercept. The target parameters in each direction at that moment are calculated, including the target coordinates of the X-axis, Y-axis, and Z-axis, as well as the target attitude angles around each axis. Finally, the target parameters of all frequent moments are arranged in chronological order to form a continuous, smooth, and optimized motion trajectory that meets the requirements of machining accuracy.

[0049] Step 702: Based on the optimized motion trajectory, generate linear motion control commands for the mobile platform and joint motion control commands for the robot through a multi-axis collaborative controller. Specifically, this includes: first, clarifying the motion responsibilities of the mobile platform and the robot to ensure they work together to complete the flange weld milling; the mobile platform is responsible for planar motion along the flange surface, i.e., adjusting its position through translation along the X and Y axes. The core purpose of this motion is to ensure the milling cutter is always aligned with the weld area, preventing deviation due to changes in weld position; the robot is responsible for adjusting the Z-axis position and orientation of the end mill, where Z-axis position adjustment controls the milling depth to prevent overcutting (milling depth exceeding process requirements) or undercutting (milling depth not meeting process requirements); orientation adjustment ensures the milling cutter axis is perpendicular to the flange surface, thus ensuring a smooth milled surface and meeting machining accuracy requirements; based on the optimized motion trajectory generated in step 701, generate linear motion control commands for the mobile platform along the X and Y axes according to the following steps, i.e., extracting the X and Y axis directions from the optimized motion trajectory. The target coordinates are the X-axis and Y-axis positions that the mobile platform needs to reach at different times (the trajectory optimization ensures continuity and smoothness). The movement speed of the mobile platform in the X-axis and Y-axis directions is calculated. The value of the movement speed is equal to the rate of change of the target coordinates with time. Since the optimized movement trajectory is smooth (optimization result of step 701), the target coordinates of a certain axis (such as the X-axis) change linearly with time. Therefore, the movement speed of that axis is the slope of the linear change function. The slope value directly reflects the average movement speed of the mobile platform in that axis direction. According to the pulse control protocol of the mobile platform drive motor, the calculated movement speed is converted into motor speed in sequence. The pulse frequency is calculated as follows: Two basic parameters are known: the distance the moving platform travels per revolution of the motor (i.e., the transmission ratio, e.g., 0.01 m / rpm), and the reduction ratio between the motor and the moving platform, e.g., 10. This means that for every 10 revolutions of the motor, the moving platform only travels the equivalent of one revolution. The motor speed is calculated as: Motor speed = Speed ​​of a certain axis of the moving platform ÷ (Distance traveled per motor revolution × Reduction ratio). For example, if the X-axis speed is 0.5 m / s, the motor travels 0.01 m per revolution, and the reduction ratio is 10, then the X-axis drive motor speed = 0.5 m / s ÷ (0.01 m / rpm × 10) = 5 rpm. Given the number of control pulses required for one revolution of the motor, i.e., the pulse equivalent parameter, such as 1000 pulses / revolution, the pulse frequency = motor speed × number of pulses required per motor revolution. This pulse frequency is directly used as the control signal parameter for driving the motor. Continuing with the previous example, if the motor requires 1000 control pulses per revolution, then the pulse frequency of the X-axis drive motor = 5 revolutions per second × 1000 pulses / revolution = 5000 Hz. Calculate the pulse frequencies of the X-axis and Y-axis drive motors of the moving platform separately using the above method, and then combine them with the motion direction signal (forward rotation corresponds to movement in the positive direction of the axis, and reverse rotation corresponds to movement in the negative direction of the axis) to form the X-axis and Y-axis linear motion control commands for the moving platform.

[0050] Taking a 6-axis robot as an example, based on the end-effector target parameters (X-axis coordinates, Y-axis coordinates, Z-axis coordinates, and attitude angles around each axis) generated in step 701, the inverse kinematics algorithm is used to calculate the target angles that each joint needs to rotate to. The specific process is as follows: Joint 1 rotates around an axis perpendicular to the flange surface to adjust the end-effector's orientation in the plane, ensuring that the end-effector is aligned with the weld direction. The input value is calculated using the arctangent function, which is the ratio of the end-effector's Y-axis coordinate to its X-axis coordinate in the optimized motion trajectory. This ratio reflects the end-effector's positional relationship in the plane. Assuming the end-effector's X-axis coordinate is 0.5 meters and its Y-axis coordinate is 0.3 meters at a certain moment in the optimized trajectory, the input value is first calculated, i.e., the ratio of the Y-axis coordinate to the X-axis coordinate = 0.3 ÷ 0. 0.5 = 0.6; then, the arctangent function is used to calculate 0.6 to obtain the target angle of joint 1. This angle ensures that after joint 1 rotates, the projection of the end effector in the plane is consistent with the weld direction; joint 3 is used to adjust the Z-axis height of the end effector to ensure that the Z-axis coordinate of the end effector meets the milling depth requirements; the input value is calculated by the inverse cosine function, which is (target Z-axis coordinate - length of link from joint 1 to joint 2 - length of link from joint 2 to joint 3) ÷ length of link from joint 3 to the end effector, where the lengths of the links from joint 1 to joint 2 (denoted as L1), from joint 2 to joint 3 (denoted as L2), and from joint 3 to the end effector (denoted as L3) are all inherent structural parameters of the robot, such as L1 = 0.2 meters, L2 = 0.3 meters, and L3 = 0.5 meters. =0.4 meters. Assuming the target Z-axis coordinate of the end point at a certain moment in the optimized trajectory is 0.6 meters, first calculate the input value: (0.6 - 0.2 - 0.3) ÷ 0.4 = 0.25; then calculate 0.25 using the inverse cosine function to obtain the target angle of joint 3. This angle ensures that after joint 3 rotates, the Z-axis coordinate of the end point reaches the milling depth requirement of 0.6 meters. The calculation logic for the target angles of joints 2, 4, 5, and 6 is similar to the above. Combining the end point attitude angle requirements (ensuring the milling cutter axis is perpendicular to the flange surface), the corresponding trigonometric functions are used for calculation. Joint 2 is used to fine-tune the pitch angle of the end point in the plane, calculated using the inverse cosine function. The input value is the ratio of the end point attitude angle around the X-axis to the link length. The end effector needs to rotate 5 degrees around the X-axis. This ratio is calculated to be 0.055. The target angle of joint 2 is obtained by calculating 0.055 using the arcsine function. Joints 4, 5, and 6 are all wrist joints used to adjust the end effector posture to ensure that the milling cutter is perpendicular to the flange surface. Joint 4 is calculated using the arctangent function. The input value is the end effector posture angle around the Y-axis. For example, if the end effector needs to rotate 5 degrees around the Y-axis, the input value is the ratio 0.087 corresponding to 5 degrees, and the target angle of joint 4 is calculated. Joints 5 and 6 are calculated using the inverse cosine function. The input value is the ratio of the end effector posture angle around the Z-axis to the length of the wrist link. For example, if the end effector needs to rotate 0 degrees around the Z-axis, the input value is 0. Finally, the target angles of joints 5 and 6 are calculated. The specific values ​​change dynamically according to the posture requirements.

[0051] Joint angular velocity reflects the speed of joint rotation, and its value is equal to the rate of change of the target angle of the corresponding joint over time. Since the optimized motion trajectory is smooth (optimization result in step 701), the target angle of the joint changes linearly with time. Therefore, the angular velocity is the slope of this linear function. Taking joint 1 as an example, if the target angle of joint 1 in the optimized trajectory is 30.96 degrees at 0.1 seconds and 32.96 degrees at 0.2 seconds, first calculate the angle change: 32.96 degrees - 30.96 degrees = 2 degrees; then calculate the time... The change in angle is 0.2 seconds - 0.1 seconds = 0.1 seconds; finally, divide the change in angle by the change in time to get the angular velocity of joint 1 = 2 degrees ÷ 0.1 seconds = 20 degrees / second. Calculate the angular velocities of the other joints in the same way. For example, joint 3 has an angle change of 0.5 degrees and a time change of 0.1 seconds, so its angular velocity is 0.5 degrees ÷ 0.1 seconds = 5 degrees / second; joint 2 has an angle change of 0.1 degrees and a time change of 0.1 seconds, so its angular velocity is 0.1 degrees ÷ 0.1 seconds = 1 degree / second, and so on. The multi-axis co-controller will process the angular velocity of each joint... Speed ​​is converted into control commands for the joint motor (taking pulse control as an example). The specific conversion process is as follows: the joint angle corresponding to one revolution of the joint motor, i.e., the reduction ratio parameter, such as 10 degrees of joint rotation for one revolution of the motor; the number of control pulses required for one revolution of the motor, i.e., the pulse equivalent, such as 1000 pulses / revolution; motor speed = joint angular velocity ÷ joint angle corresponding to one revolution of the motor. Taking joint 1 as an example, motor speed = 20 degrees / second ÷ 10 degrees / revolution = 2 revolutions / second; pulse frequency = motor speed × motor rotation rate. The number of control pulses required for one revolution, taking joint 1 as an example, is: pulse frequency = 2 revolutions / second × 1000 pulses / revolution = 2000 Hz. Calculate the pulse frequencies for joints 2 to 6 using the same method. For example, the motor speed of joint 3 is 5 degrees / second ÷ 10 degrees / revolution = 0.5 revolutions / second, and the pulse frequency is 0.5 revolutions / second × 1000 pulses / revolution = 500 Hz. Then, integrate the pulse frequencies of all joints with the rotation direction signal (forward rotation corresponds to an increase in joint angle, and reverse rotation corresponds to a decrease in joint angle) to form the robot's joint motion control command.

[0052] Step 703: Based on the linear motion control commands and joint motion control commands, control the moving platform to move along the flange surface and the robot joints to complete the milling of the flange weld. Specifically, this includes: the multi-axis collaborative controller sending linear motion control commands to the drive motors of the moving platform, such as the X-axis and Y-axis servo motors, and simultaneously sending joint motion control commands to the joint motors of the robot; using the controller's internal synchronization clock, such as synchronizing every 1 millisecond, to ensure that the moving platform and the robot start moving simultaneously, avoiding milling cutter deviation from the weld position due to start-up time difference, thus affecting machining accuracy; the moving platform using position sensors, such as grating rulers, to collect the actual X-axis and Y-axis positions in real time, and comparing the actual positions with the target positions at the corresponding moments in the optimized motion trajectory; if the deviation is less than 0.001 meters, for example, if the target X-axis position at this moment in the optimized trajectory is 0.5036 meters and the actual collected X-axis position is 0.503 meters, the deviation is 0.0006 meters, then no adjustment is needed; if the deviation reaches or exceeds 0.001 meters, the controller immediately corrects the pulses in the linear motion control commands. The frequency is adjusted to control the movement speed of the mobile platform, ensuring the actual position returns to the target position. The robot collects the actual angles of each joint in real time through joint encoders and compares the actual angles with the target joint angles at the corresponding moments in the optimized motion trajectory. If the deviation is less than 0.05 degrees, no adjustment is needed. If the deviation reaches or exceeds 0.05 degrees, the controller corrects the angular velocity in the joint motion control commands and adjusts the rotation speed of the joint motors to ensure the actual angle returns to the target angle, guaranteeing that the robot's end effector pose meets the optimized trajectory requirements. As the mobile platform moves along the flange surface according to linear motion control commands (covering the entire length of the weld seam), and the robot adjusts the Z-axis position of the end mill according to joint motion control commands (to achieve the preset milling depth, such as decreasing from 0.1 meters to 0.09 meters to achieve a milling depth of 0.01 meters) and posture, the milling cutter continuously cuts the flange weld seam. When the mobile platform completes the coverage of all weld seam areas and the robot confirms that the preset milling depth has been reached, the multi-axis collaborative controller issues a stop command, the mobile platform and the robot stop moving, and the milling of the flange weld seam is completed.

[0053] In this embodiment, the pre-compensation trajectory is corrected by adjusting the pose to offset trajectory deviations caused by load vibrations. Trajectory optimization further improves trajectory smoothness, avoiding uneven milling depths due to trajectory fluctuations and ensuring that the surface accuracy of the weld after milling meets requirements. The motion task is decomposed by a multi-axis collaborative controller, clarifying the division of labor between the mobile platform and the robot and avoiding motion interference. Synchronous start-up and real-time adjustment ensure coordinated movement between the two, eliminating the need for frequent pauses for calibration and improving overall processing efficiency. Trajectory optimization reduces speed abrupt changes, preventing vibrations in the mobile platform and robot joints due to impacts and reducing equipment wear. Real-time monitoring and deviation correction prevent the milling cutter from rigidly colliding with the flange due to trajectory deviation, protecting the tool from damage. The entire process is based on optimized trajectory control of the flange weld.

[0054] like Figure 2 As shown, embodiments of the present invention also provide a multi-axis cooperative automatic control system for milling flange weld seams, comprising: The acquisition module is used to acquire three-dimensional contour data of the flange circumferential seam area and detect the pose deviation of the moving platform in real time to obtain the base error data. The extraction module is used to extract the geometric feature points of the circumferential seam based on the three-dimensional contour data, obtain two feature lines, and calculate the included angle between the two feature lines to define the local processing area of ​​the circumferential seam. The module is used to set first and second control points inside and outside the local processing area, and to construct a closed loop control path based on the time series data of the first and second control points. The compensation module is used to generate trajectory correction parameters based on the geometric characteristics of the circular closed control path, so as to correct the transformation relationship between the robot's base coordinate system and the theoretical programming coordinate system, and generate a real-time compensation signal for the base coordinate system. The feedback module is used to adjust the robot's joint motion parameters based on the real-time compensation signal of the base coordinate system, obtain the pre-compensated motion trajectory, and collect the force between the milling cutter and the workpiece to obtain milling force feedback data. The identification module is used to identify load changes and vibration characteristics based on milling force feedback data and spindle motor current signals, determine load vibration characteristic data, adjust the robot's end position and attitude, and generate force control signals. The optimization module is used to correct the pre-compensated motion trajectory using force control signals to obtain the optimized motion trajectory. Based on the optimized motion trajectory, the linear motion of the mobile platform and the joint motion of the robot are controlled by a multi-axis collaborative controller to complete the milling of the flange weld.

[0055] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0056] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0057] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A multi-axis collaborative automatic control method for milling flange weld seams, characterized in that, The method includes: Step 1: Collect the three-dimensional contour data of the flange circumferential seam area and detect the pose deviation of the moving platform in real time to obtain the base error data; Step 2: Based on the 3D contour data, extract the geometric feature points of the circumferential seam to obtain two feature lines, and calculate the included angle between the two feature lines to define the local processing area of ​​the circumferential seam. Step 3: Set up first and second control points inside and outside the local processing area, and construct a closed loop control path based on the time series data of the first and second control points; Step 4: Based on the geometric characteristics of the circular closed control path, generate trajectory correction parameters to correct the transformation relationship between the robot's base coordinate system and the theoretical programming coordinate system, and generate a real-time compensation signal for the base coordinate system. Step 5: Based on the real-time compensation signal of the base coordinate system, adjust the joint motion parameters of the robot to obtain the pre-compensated motion trajectory, and collect the interaction force between the milling cutter and the workpiece to obtain milling force feedback data; Step 6: Based on the milling force feedback data and the current signal of the spindle motor, identify the load changes and vibration characteristics, determine the load vibration characteristic data, adjust the robot's end-effector pose, and generate force control signals. Step 7: Correct the pre-compensated motion trajectory using force control signals to obtain the optimized motion trajectory; based on the optimized motion trajectory, control the linear motion of the mobile platform and the joint motion of the robot through a multi-axis collaborative controller to complete the milling of the flange weld.

2. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 1, characterized in that, Step 1 includes: Acceleration and angular velocity data of the mobile platform are collected by a multi-degree-of-freedom inertial measurement unit installed on the mobile platform, and the real-time pose data of the mobile platform is obtained by attitude calculation based on the acceleration and angular velocity data. The flange circumferential seam area is scanned by a laser scanner to obtain three-dimensional contour data of the seam surface. The real-time pose data is compared with the preset theoretical pose data to obtain the comparison result. The pose deviation of the moving platform is calculated based on the comparison result to generate the base error data.

3. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 2, characterized in that, Step 2 includes: The 3D contour data is preprocessed to obtain optimized 3D contour data; Geometric feature points of the circumferential seam edge are extracted from the optimized 3D contour data. Based on the spatial coordinates of the geometric feature points, the spatial position parameters of the two reference feature lines are obtained by fitting using the least squares method. The included angle between the two reference feature lines is calculated based on their spatial position parameters. The processing range boundary of the circumferential seam is determined based on the included angle value, and the circumferential seam segment within the included angle range is defined as the local processing area.

4. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 3, characterized in that, Step 3 includes: Feature points on the center line of the circumferential seam within the local processing area are selected as the first control points, and feature points on the reference edge of the flange base material outside the local processing area are selected as the second control points. Record the continuous position data of the first and second control points during the milling process to form a time series dataset; generate a smooth and continuous circular closed control path based on the continuous position data in the time series dataset.

5. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 4, characterized in that, Step 4 includes: Based on the curvature characteristics and normal vector distribution of the circular closed control path, the geometric characteristic parameters of the circular closed control path are calculated, and the trajectory correction parameters are generated according to the path geometric characteristic parameters. By fusing trajectory correction parameters with basis error data, the transformation matrix between the robot's base coordinate system and the theoretical programming coordinate system is corrected, resulting in the corrected transformation matrix. Based on the corrected transformation matrix, a real-time compensation signal for the base coordinate system is generated.

6. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 5, characterized in that, Step 5 includes: Based on the real-time compensation signal of the base coordinate system, the angle compensation amount of each joint of the robot is calculated; based on the angle compensation amount, the joint motion parameters of the robot are adjusted to generate a pre-compensated motion trajectory. During the robot's movement along the pre-compensated motion trajectory, a six-dimensional force sensor collects real-time three-dimensional force data between the milling cutter and the workpiece, and filters the three-dimensional force data to obtain milling force feedback data.

7. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 6, characterized in that, Step 6 includes: Frequency domain analysis was performed on the milling force feedback data to extract the vibration characteristic frequencies and amplitudes; harmonic analysis was performed on the current signal of the spindle motor to identify load variation characteristics; By fusing the vibration characteristic frequency and amplitude with the load change characteristics, load vibration characteristic data is obtained; Based on load vibration characteristic data, the pose adjustment amount of the robot end effector is calculated through an impedance control algorithm, and a force control signal is generated according to the pose adjustment amount.

8. The multi-axis collaborative automatic control method for milling flange weld seams according to claim 7, characterized in that, Step 7 includes: Based on the pose adjustment amount contained in the force control signal, the pre-compensated motion trajectory is corrected to obtain trajectory correction data; Based on trajectory correction data, trajectory optimization calculations are performed on the pre-compensated motion trajectory to generate trajectory optimization parameters; and based on the trajectory optimization parameters, the optimized motion trajectory is generated. Based on optimized motion trajectories, linear motion control commands for the mobile platform and joint motion control commands for the robot are generated through a multi-axis collaborative controller. Based on linear motion control commands and joint motion control commands, the mobile platform is controlled to move along the flange surface and the robot joints are controlled to complete the milling of the flange weld.

9. A multi-axis collaborative automatic control system for milling flange weld seams, the system implementing the method as described in any one of claims 1 to 8, characterized in that, include: The acquisition module is used to acquire three-dimensional contour data of the flange circumferential seam area and detect the pose deviation of the moving platform in real time to obtain the base error data. The extraction module is used to extract the geometric feature points of the circumferential seam based on the three-dimensional contour data, obtain two feature lines, and calculate the included angle between the two feature lines to define the local processing area of ​​the circumferential seam. The module is used to set first and second control points inside and outside the local processing area, and to construct a closed loop control path based on the time series data of the first and second control points. The compensation module is used to generate trajectory correction parameters based on the geometric characteristics of the circular closed control path, so as to correct the transformation relationship between the robot's base coordinate system and the theoretical programming coordinate system, and generate a real-time compensation signal for the base coordinate system. The feedback module is used to adjust the robot's joint motion parameters based on the real-time compensation signal of the base coordinate system, obtain the pre-compensated motion trajectory, and collect the force between the milling cutter and the workpiece to obtain milling force feedback data. The identification module is used to identify load changes and vibration characteristics based on milling force feedback data and spindle motor current signals, determine load vibration characteristic data, adjust the robot's end-effector pose, and generate force control signals. The optimization module is used to correct the pre-compensated motion trajectory using force control signals to obtain the optimized motion trajectory. Based on the optimized motion trajectory, the linear motion of the mobile platform and the joint motion of the robot are controlled by a multi-axis collaborative controller to complete the milling of the flange weld.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 8.

Images